diff --git a/GNUmakefile b/GNUmakefile index fe93eb5..ff0cdac 100644 --- a/GNUmakefile +++ b/GNUmakefile @@ -46,6 +46,9 @@ CPPFLAGS += -DG4MULTITHREADED #CPPFLAGS += -DDEBUG_SECTIONPLANE #CPPFLAGS += -DDEBUG_SECTIONPLANE_ZAVE #CPPFLAGS += -DG4VERSION_10_04_OR_LATER=1 +ifneq (, $(wildcard $(HALLD_RECON_HOME)/src/libraries/DIRC/DDIRCPmtHit.*)) + CPPFLAGS += -DDIRCTRUTHEXTRA +endif # If you want to build against Geant4.10.03 or greater, you will need this line uncommented #CPPFLAGS += -DG4VUSERPHYSICSLIST_HAS_GETPARTICLEITERATOR diff --git a/src/GlueXDetectorConstruction.cc b/src/GlueXDetectorConstruction.cc index 084d077..8176da1 100644 --- a/src/GlueXDetectorConstruction.cc +++ b/src/GlueXDetectorConstruction.cc @@ -368,13 +368,29 @@ void GlueXDetectorConstruction::ConstructSDandField() } iter->second->SetSensitiveDetector(ftofHandler); } - else if (volname == "QZBL" || volname == "PIXV") { + // radiator volume: BNNM (NN = bar number 0-47 and M is sub-bar character A-D) + else if (volname == "PIXV" || + (volname(0,1)(0) == 'B' && + 10*((int)volname(1,1)(0)-48)+(int)volname(2,1)(0)-48 >= 0 && + 10*((int)volname(1,1)(0)-48)+(int)volname(2,1)(0)-48 < 48)) + { // this is nasty, but it works if (dircHandler == 0) { dircHandler = new GlueXSensitiveDetectorDIRC("dirc"); SDman->AddNewDetector(dircHandler); } iter->second->SetSensitiveDetector(dircHandler); } + else if (volname == "FWM1" || volname == "FWM2" || volname == "FTMR" || + volname == "TSM1" || volname == "TSM2" || volname == "TSM3" || + volname == "FSM1" || volname == "FSM2" || volname == "OWDG" || + (volname(0,1)(0) == 'A' && volname(0,1)(1) == 'G') ) + { + if (dircHandler == 0) { + dircHandler = new GlueXSensitiveDetectorDIRC("dirc"); + SDman->AddNewDetector(dircHandler); + } + iter->second->SetSensitiveDetector(dircHandler); + } else if (volname == "CERW" || volname == "CPPC") { if (cereHandler == 0) { cereHandler = new GlueXSensitiveDetectorCERE("cere"); diff --git a/src/GlueXHitDIRCPmt.cc b/src/GlueXHitDIRCPmt.cc index 3f41164..7ac4878 100644 --- a/src/GlueXHitDIRCPmt.cc +++ b/src/GlueXHitDIRCPmt.cc @@ -13,11 +13,15 @@ GlueXHitDIRCPmt::GlueXHitDIRCPmt(const GlueXHitDIRCPmt &src) { E_GeV = src.E_GeV; t_ns = src.t_ns; + t_fixed_ns = src.t_fixed_ns; x_cm = src.x_cm; y_cm = src.y_cm; z_cm = src.z_cm; ch = src.ch; key_bar = src.key_bar; + path = src.path; + refl = src.refl; + bbrefl = src.bbrefl; } void GlueXHitDIRCPmt::Draw() const @@ -29,9 +33,13 @@ void GlueXHitDIRCPmt::Print() const { G4cout << " E = " << E_GeV << " GeV" << G4endl << " t = " << t_ns << " ns" << G4endl + << " t_fixed = " << t_fixed_ns << " ns" << G4endl << " x = " << x_cm << " cm" << G4endl << " y = " << y_cm << " cm" << G4endl << " z = " << z_cm << " cm" << G4endl + << " path = " << path << " " << G4endl + << " refl = " << refl << " " << G4endl + << " bbrefl = " << bbrefl << " " << G4endl << " ch = " << ch << " " << G4endl << " key_bar = " << key_bar << " " << G4endl << G4endl; diff --git a/src/GlueXHitDIRCPmt.hh b/src/GlueXHitDIRCPmt.hh index 8422460..afd9c47 100644 --- a/src/GlueXHitDIRCPmt.hh +++ b/src/GlueXHitDIRCPmt.hh @@ -25,11 +25,16 @@ public: G4double E_GeV; // photon energy [GeV] G4double t_ns; // detection time [ns] + G4double t_fixed_ns; // fixed pathlength time [ns] G4double x_cm; // x coordinate where hit was created G4double y_cm; // y coordinate where hit was created G4double z_cm; // z coordinate where hit was created + int64_t path; // photon's path id in the optical box + G4int refl; // number of reflections in the oprical box + G4bool bbrefl; // reflected off far end mirror of bar box G4int ch; // PMT channel of the hit G4int key_bar; // key of the corresponding bar hit + G4int track; // index of the MC track }; typedef G4THitsMap GlueXHitsMapDIRCPmt; diff --git a/src/GlueXHitDIRCWob.cc b/src/GlueXHitDIRCWob.cc new file mode 100644 index 0000000..2269660 --- /dev/null +++ b/src/GlueXHitDIRCWob.cc @@ -0,0 +1,26 @@ +// +// GlueXHitDIRCWob - class implementation +// +// created on: 05.04.2017 +// original author: r.dzhygadlo at gsi.de + +#include "GlueXHitDIRCWob.hh" + +G4ThreadLocal G4Allocator* GlueXHitDIRCWobAllocator = 0; + +GlueXHitDIRCWob::GlueXHitDIRCWob():G4VHit() +{ +} + +void GlueXHitDIRCWob::Draw() const +{ + // not yet implemented +} + +void GlueXHitDIRCWob::Print() const +{ + G4cout << "GlueXHitDIRCWob:" << G4endl + << " track = " << track << G4endl + << " normalId = " << normalId << " " << G4endl + << G4endl; +} diff --git a/src/GlueXHitDIRCWob.hh b/src/GlueXHitDIRCWob.hh new file mode 100644 index 0000000..3e71bea --- /dev/null +++ b/src/GlueXHitDIRCWob.hh @@ -0,0 +1,45 @@ +// +// GlueXHitDIRCWob - class implementation +// +// created on: 05.04.2017 +// original author: r.dzhygadlo at gsi.de + +#ifndef GlueXHitDIRCWob_h +#define GlueXHitDIRCWob_h 1 + +#include "G4VHit.hh" +#include "G4THitsMap.hh" +#include "G4Allocator.hh" + +class GlueXHitDIRCWob : public G4VHit +{ +public: + GlueXHitDIRCWob(); + + void *operator new(size_t); + void operator delete(void *aHit); + + void Draw() const; + void Print() const; + + G4double normalId; // uniq identifier of the normal + G4int track; // index of the MC track +}; + +typedef G4THitsMap GlueXHitsMapDIRCWob; + +extern G4ThreadLocal G4Allocator* GlueXHitDIRCWobAllocator; + +inline void* GlueXHitDIRCWob::operator new(size_t) +{ + if (!GlueXHitDIRCWobAllocator) + GlueXHitDIRCWobAllocator = new G4Allocator; + return (void *) GlueXHitDIRCWobAllocator->MallocSingle(); +} + +inline void GlueXHitDIRCWob::operator delete(void *aHit) +{ + GlueXHitDIRCWobAllocator->FreeSingle((GlueXHitDIRCWob*) aHit); +} + +#endif diff --git a/src/GlueXHitDIRCpoint.cc b/src/GlueXHitDIRCpoint.cc deleted file mode 100644 index c72f682..0000000 --- a/src/GlueXHitDIRCpoint.cc +++ /dev/null @@ -1,88 +0,0 @@ -// -// GlueXHitDIRCpoint - class implementation -// -// author: richard.t.jones at uconn.edu -// version: november 26, 2016 - -#include "GlueXHitDIRCpoint.hh" - -G4ThreadLocal G4Allocator* GlueXHitDIRCpointAllocator = 0; - -GlueXHitDIRCpoint::GlueXHitDIRCpoint(const GlueXHitDIRCpoint &src) -{ - E_GeV = src.E_GeV; - primary_ = src.primary_; - ptype_G3 = src.ptype_G3; - px_GeV = src.px_GeV; - py_GeV = src.py_GeV; - pz_GeV = src.pz_GeV; - x_cm = src.x_cm; - y_cm = src.y_cm; - z_cm = src.z_cm; - t_ns = src.t_ns; - track_ = src.track_; - trackID_ = src.trackID_; -} - -int GlueXHitDIRCpoint::operator==(const GlueXHitDIRCpoint &right) const -{ - if (E_GeV != right.E_GeV || - primary_ != right.primary_ || - ptype_G3 != right.ptype_G3 || - px_GeV != right.px_GeV || - py_GeV != right.py_GeV || - pz_GeV != right.pz_GeV || - x_cm != right.x_cm || - y_cm != right.y_cm || - z_cm != right.z_cm || - t_ns != right.t_ns || - track_ != right.track_ || - trackID_ != right.trackID_ ) - { - return 0; - } - return 1; -} - -GlueXHitDIRCpoint &GlueXHitDIRCpoint::operator+=(const GlueXHitDIRCpoint &right) -{ - G4cerr << "Error in GlueXHitDIRCpoint::operator+= - " - << "illegal attempt to merge two TruthPoint objects in the DIRC!" - << G4endl; - return *this; -} - -void GlueXHitDIRCpoint::Draw() const -{ - // not yet implemented -} - -void GlueXHitDIRCpoint::Print() const -{ - G4cout << "GlueXHitDIRCpoint:" << G4endl - << " track = " << track_ << G4endl - << " trackID = " << trackID_ << G4endl - << " E = " << E_GeV << " GeV" << G4endl - << " primary = " << primary_ << G4endl - << " ptype = " << ptype_G3 << G4endl - << " px = " << px_GeV << " GeV/c" << G4endl - << " py = " << py_GeV << " GeV/c" << G4endl - << " pz = " << pz_GeV << " GeV/c" << G4endl - << " x = " << x_cm << " cm" << G4endl - << " y = " << y_cm << " cm" << G4endl - << " z = " << z_cm << " cm" << G4endl - << " t = " << t_ns << " ns" << G4endl - << G4endl; -} - -void printallhits(GlueXHitsMapDIRCpoint *hitsmap) -{ - std::map *map = hitsmap->GetMap(); - std::map::const_iterator iter; - G4cout << "G4THitsMap " << hitsmap->GetName() << " with " << hitsmap->entries() - << " entries:" << G4endl; - for (iter = map->begin(); iter != map->end(); ++iter) { - G4cout << " key=" << iter->first << " "; - iter->second->Print(); - } -} diff --git a/src/GlueXHitDIRCpoint.hh b/src/GlueXHitDIRCpoint.hh deleted file mode 100644 index 52ad32c..0000000 --- a/src/GlueXHitDIRCpoint.hh +++ /dev/null @@ -1,67 +0,0 @@ -// -// GlueXHitDIRCpoint - class header -// -// author: richard.t.jones at uconn.edu -// version: november 26, 2016 -// -// In the context of the Geant4 event-level multithreading model, -// this class is "thread-local", ie. has thread-local state. Its -// allocator is designed to run within a worker thread context. -// This class is final, do NOT try to derive another class from it. - -#ifndef GlueXHitDIRCpoint_h -#define GlueXHitDIRCpoint_h 1 - -#include "G4VHit.hh" -#include "G4THitsMap.hh" -#include "G4Allocator.hh" - -class GlueXHitDIRCpoint : public G4VHit -{ - public: - GlueXHitDIRCpoint() {} - GlueXHitDIRCpoint(const GlueXHitDIRCpoint &src); - int operator==(const GlueXHitDIRCpoint &right) const; - GlueXHitDIRCpoint &operator+=(const GlueXHitDIRCpoint &right); - - void *operator new(size_t); - void operator delete(void *aHit); - - void Draw() const; - void Print() const; - - // no reason to hide hit data - - G4double E_GeV; // total energy (GeV) of this track at this point - G4bool primary_; // true if track belongs to from a primary particle - G4int ptype_G3; // G3 type of particle making this track - G4double px_GeV; // momentum (GeV/c) of track at point, x component - G4double py_GeV; // momentum (GeV/c) of track at point, y component - G4double pz_GeV; // momentum (GeV/c) of track at point, z component - G4double x_cm; // global x coordinate of track at point (cm) - G4double y_cm; // global y coordinate of track at point (cm) - G4double z_cm; // global z coordinate of track at point (cm) - G4double t_ns; // time of track crossing at point (ns) - G4int track_; // Geant4 track ID of particle making this track - G4int trackID_; // GlueX-assigned track ID of particle making this track - - G4int GetKey() const { return (track_ << 20) + int(t_ns * 100); } -}; - -typedef G4THitsMap GlueXHitsMapDIRCpoint; - -extern G4ThreadLocal G4Allocator* GlueXHitDIRCpointAllocator; - -inline void* GlueXHitDIRCpoint::operator new(size_t) -{ - if (!GlueXHitDIRCpointAllocator) - GlueXHitDIRCpointAllocator = new G4Allocator; - return (void *) GlueXHitDIRCpointAllocator->MallocSingle(); -} - -inline void GlueXHitDIRCpoint::operator delete(void *aHit) -{ - GlueXHitDIRCpointAllocator->FreeSingle((GlueXHitDIRCpoint*) aHit); -} - -#endif diff --git a/src/GlueXPrimaryGeneratorAction.cc b/src/GlueXPrimaryGeneratorAction.cc index 8d981c4..6427f59 100644 --- a/src/GlueXPrimaryGeneratorAction.cc +++ b/src/GlueXPrimaryGeneratorAction.cc @@ -8,6 +8,7 @@ #include "GlueXPrimaryGenerator.hh" #include "GlueXUserEventInformation.hh" #include "GlueXUserOptions.hh" +#include "G4OpticalPhoton.hh" #include "G4Event.hh" #include "G4ParticleTable.hh" @@ -72,6 +73,93 @@ GlueXPrimaryGeneratorAction::GlueXPrimaryGeneratorAction() exit(-1); } + // get positions for LUT from XML geometry + std::map dirclutpars; + if (instanceCount == 1) { + + if (user_opts->Find("DIRCLUT", dirclutpars)) { + + extern int run_number; + extern jana::JApplication *japp; + if (japp == 0) { + G4cerr << "Error in GlueXPrimaryGeneratorAction constructor - " + << "jana global DApplication object not set, " + << "cannot continue." << G4endl; + exit(-1); + } + jana::JGeometry *jgeom = japp->GetJGeometry(run_number); + if (japp == 0) { // dummy + jgeom = 0; + G4cout << "DIRC: ALL parameters loaded from ccdb" << G4endl; + } + + vectorDIRC; + vectorDRCC; + vectorDCML00_XYZ; + vectorDCML01_XYZ; + vectorDCML10_XYZ; + vectorDCML11_XYZ; + vectorWNGL00_XYZ; + vectorWNGL01_XYZ; + vectorWNGL10_XYZ; + vectorWNGL11_XYZ; + vectorOWDG_XYZ; + jgeom->Get("//section/composition/posXYZ[@volume='DIRC']/@X_Y_Z", DIRC); + jgeom->Get("//composition[@name='DIRC']/posXYZ[@volume='DRCC']/@X_Y_Z", DRCC); + jgeom->Get("//composition[@name='DRCC']/posXYZ[@volume='DCML00']/@X_Y_Z", DCML00_XYZ); + jgeom->Get("//composition[@name='DRCC']/posXYZ[@volume='DCML01']/@X_Y_Z", DCML01_XYZ); + jgeom->Get("//composition[@name='DRCC']/posXYZ[@volume='DCML10']/@X_Y_Z", DCML10_XYZ); + jgeom->Get("//composition[@name='DRCC']/posXYZ[@volume='DCML11']/@X_Y_Z", DCML11_XYZ); + jgeom->Get("//composition[@name='DCML00']/posXYZ[@volume='WNGL']/@X_Y_Z", WNGL00_XYZ); + jgeom->Get("//composition[@name='DCML01']/posXYZ[@volume='WNGL']/@X_Y_Z", WNGL01_XYZ); + jgeom->Get("//composition[@name='DCML10']/posXYZ[@volume='WNGL']/@X_Y_Z", WNGL10_XYZ); + jgeom->Get("//composition[@name='DCML11']/posXYZ[@volume='WNGL']/@X_Y_Z", WNGL11_XYZ); + jgeom->Get("//composition[@name='DCML11']/posXYZ[@volume='WNGL']/@X_Y_Z", WNGL11_XYZ); + jgeom->Get("//trd[@name='OWDG']/@Xmp_Ymp_Z", OWDG_XYZ); + DIRC_LUT_Z = (DIRC[2] + DRCC[2] + DCML01_XYZ[2] + 0.8625) * cm; + DIRC_QZBL_DY = 3.5 * cm; // nominal width to generate LUT + DIRC_QZBL_DZ = 1.725 * cm; // nominal thickness to generate LUT + DIRC_OWDG_DZ = OWDG_XYZ[4]; + + // set array of bar positions + for (int i=0; i<48; i++) { + vectorDCBR_XYZ; + if (i<12) { + std::stringstream geomDCML10; + geomDCML10 << "//composition[@name='DCML10']/posXYZ[@volume='DCBR" + << std::setfill('0') << std::setw(2) << i << "']/@X_Y_Z"; + jgeom->Get(geomDCML10.str(), DCBR_XYZ); + DIRC_BAR_Y[i] = (DCML10_XYZ[1] - DCBR_XYZ[1]) * cm; + DIRC_LUT_X[i] = (DIRC[0] + DRCC[0] + DCML10_XYZ[0] - WNGL10_XYZ[0] + DIRC_OWDG_DZ) * cm; + } + else if (i<24) { + std::stringstream geomDCML11; + geomDCML11 << "//composition[@name='DCML11']/posXYZ[@volume='DCBR" + << std::setfill('0') << std::setw(2) << i << "']/@X_Y_Z"; + jgeom->Get(geomDCML11.str(), DCBR_XYZ); + DIRC_BAR_Y[i] = (DCML11_XYZ[1] - DCBR_XYZ[1]) * cm; + DIRC_LUT_X[i] = (DIRC[0] + DRCC[0] + DCML11_XYZ[0] - WNGL11_XYZ[0] + DIRC_OWDG_DZ) * cm; + } + else if (i<36) { + std::stringstream geomDCML01; + geomDCML01 << "//composition[@name='DCML01']/posXYZ[@volume='DCBR" + << std::setfill('0') << std::setw(2) << i << "']/@X_Y_Z"; + jgeom->Get(geomDCML01.str(), DCBR_XYZ); + DIRC_BAR_Y[i] = (DCML01_XYZ[1] + DCBR_XYZ[1]) * cm; + DIRC_LUT_X[i] = (DIRC[0] + DRCC[0] + DCML01_XYZ[0] + WNGL01_XYZ[0] - DIRC_OWDG_DZ) * cm; + } + else if (i<48) { + std::stringstream geomDCML00; + geomDCML00 << "//composition[@name='DCML00']/posXYZ[@volume='DCBR" + << std::setfill('0') << std::setw(2) << i << "']/@X_Y_Z"; + jgeom->Get(geomDCML00.str(), DCBR_XYZ); + DIRC_BAR_Y[i] = (DCML00_XYZ[1] + DCBR_XYZ[1]) * cm; + DIRC_LUT_X[i] = (DIRC[0] + DRCC[0] + DCML00_XYZ[0] + WNGL00_XYZ[0] - DIRC_OWDG_DZ) * cm; + } + } + } + } + std::map infile; std::map beampars; std::map kinepars; @@ -102,6 +190,23 @@ GlueXPrimaryGeneratorAction::GlueXPrimaryGeneratorAction() fSourceType = SOURCE_TYPE_COBREMS_GEN; } + else if (user_opts->Find("DIRCLUT", dirclutpars)) + { + fGunParticle.geantType = 0; + fGunParticle.pdgType = 999999; + fGunParticle.partDef = fParticleTable->FindParticle("opticalphoton"); + fGunParticle.deltaR = 0; + fGunParticle.deltaZ = 0; + fGunParticle.mom = 3.18 * eV; + + fGunParticle.deltaMom = 0; + fGunParticle.deltaTheta = 0; + fGunParticle.deltaPhi = 0; + fParticleGun->SetParticleDefinition(fGunParticle.partDef); + + fSourceType = SOURCE_TYPE_PARTICLE_GUN; + } + else if (user_opts->Find("KINE", kinepars)) { if (kinepars[1] == 1000) { @@ -329,6 +434,11 @@ void GlueXPrimaryGeneratorAction::GeneratePrimariesParticleGun(G4Event* anEvent) // our own derived class. (Sheesh!!) fParticleGun->Reset(); + // std::cout<<"GlueXPrimaryGeneratorAction:: GENERATE PRIMARIES PARTICLE GUN"< dirclutpars; + // place and smear the particle gun origin G4ThreeVector pos(fGunParticle.pos); if (fGunParticle.deltaR > 0) { @@ -378,6 +488,35 @@ void GlueXPrimaryGeneratorAction::GeneratePrimariesParticleGun(G4Event* anEvent) } if (fGunParticle.deltaPhi > 0) phip += (G4UniformRand() - 0.5) * fGunParticle.deltaPhi; + + // Special case of Cherenkov photon gun for DIRC Look Up Tables (LUT) + if (user_opts->Find("DIRCLUT", dirclutpars)) { + + // array of bar y-positions for LUT from JGeometry + double y = 0.; // no shift + double x = DIRC_LUT_X[dirclutpars[1]]; + double z = DIRC_LUT_Z; + + G4ThreeVector vec(0,0,1); + double rand1 = G4UniformRand(); + double rand2 = G4UniformRand(); + vec.setTheta(acos(rand1)); + vec.setPhi(2*M_PI*rand2); + vec.rotateY(M_PI/2.); + y = DIRC_BAR_Y[dirclutpars[1]]; + if (dirclutpars[1] < 24) { + vec.rotateY(M_PI); + } + + // spread over end of bar in y and z + y += DIRC_QZBL_DY/2.0 - DIRC_QZBL_DY*G4UniformRand(); + z += DIRC_QZBL_DZ/2.0 - DIRC_QZBL_DZ*G4UniformRand(); + + thetap = vec.theta(); + phip = vec.phi(); + fParticleGun->SetParticlePosition(G4ThreeVector(x,y,z)); + } + G4ThreeVector mom(p * sin(thetap) * cos(phip), p * sin(thetap) * sin(phip), p * cos(thetap)); diff --git a/src/GlueXPrimaryGeneratorAction.hh b/src/GlueXPrimaryGeneratorAction.hh index f45e67b..f0de3da 100644 --- a/src/GlueXPrimaryGeneratorAction.hh +++ b/src/GlueXPrimaryGeneratorAction.hh @@ -26,6 +26,7 @@ #include "G4Event.hh" #include +#include #include @@ -142,6 +143,9 @@ class GlueXPrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction private: static G4Mutex fMutex; static std::list fInstance; + double DIRC_LUT_X[48], DIRC_BAR_Y[48]; + double DIRC_LUT_Z; + double DIRC_QZBL_DY, DIRC_QZBL_DZ, DIRC_OWDG_DZ; }; inline G4ParticleDefinition *GlueXPrimaryGeneratorAction::GetParticle(int PDGtype) diff --git a/src/GlueXSensitiveDetectorDIRC.cc b/src/GlueXSensitiveDetectorDIRC.cc index 89d5460..f4f0d24 100644 --- a/src/GlueXSensitiveDetectorDIRC.cc +++ b/src/GlueXSensitiveDetectorDIRC.cc @@ -9,6 +9,7 @@ #include "GlueXPrimaryGeneratorAction.hh" #include "GlueXUserEventInformation.hh" #include "GlueXUserTrackInformation.hh" +#include "GlueXUserOptions.hh" #include "G4VPhysicalVolume.hh" #include "G4PVPlacement.hh" @@ -17,6 +18,8 @@ #include "G4Step.hh" #include "G4SDManager.hh" #include "G4ios.hh" +#include "G4TransportationManager.hh" +#include "G4ParallelWorldProcess.hh" #include @@ -24,7 +27,8 @@ #define OPTICAL_PHOTON 50 // Cutoff on the total number of allowed hits -int GlueXSensitiveDetectorDIRC::MAX_HITS = 500; +int GlueXSensitiveDetectorDIRC::MAX_HITS = 1000; +int GlueXSensitiveDetectorDIRC::MAX_PIXELS = 6912; // Minimum hit time difference for two hits on the same tube double GlueXSensitiveDetectorDIRC::TWO_HIT_TIME_RESOL = 50*ns; @@ -49,8 +53,8 @@ GlueXSensitiveDetectorDIRC::GlueXSensitiveDetectorDIRC(const G4String& name) extern jana::JApplication *japp; if (japp == 0) { G4cerr << "Error in GlueXSensitiveDetector constructor - " - << "jana global DApplication object not set, " - << "cannot continue." << G4endl; + << "jana global DApplication object not set, " + << "cannot continue." << G4endl; exit(-1); } jana::JCalibration *jcalib = japp->GetJCalibration(run_number); @@ -59,6 +63,15 @@ GlueXSensitiveDetectorDIRC::GlueXSensitiveDetectorDIRC(const G4String& name) G4cout << "DIRC: ALL parameters loaded from ccdb" << G4endl; } } + + GlueXUserOptions *user_opts = GlueXUserOptions::GetInstance(); + std::map dirclutpars; + if (user_opts->Find("DIRCLUT", dirclutpars)) { + fLutId = dirclutpars[1]; + } + else { + fLutId = 100; + } } GlueXSensitiveDetectorDIRC::GlueXSensitiveDetectorDIRC(const GlueXSensitiveDetectorDIRC &src) @@ -69,7 +82,7 @@ GlueXSensitiveDetectorDIRC::GlueXSensitiveDetectorDIRC(const G4String& name) } GlueXSensitiveDetectorDIRC &GlueXSensitiveDetectorDIRC::operator=(const - GlueXSensitiveDetectorDIRC &src) + GlueXSensitiveDetectorDIRC &src) { G4AutoLock barrier(&fMutex); *(G4VSensitiveDetector*)this = src; @@ -117,13 +130,15 @@ G4bool GlueXSensitiveDetectorDIRC::ProcessHits(G4Step* step, // order of appearance in the event simulation. G4Track *track = step->GetTrack(); - // G4int trackID = track->GetTrackID(); GlueXUserTrackInformation *trackinfo = (GlueXUserTrackInformation*) track->GetUserInformation(); int itrack = trackinfo->GetGlueXTrackID(); + G4String volname = touch->GetVolume()->GetName(); + + // radiator volume: BNNM (NN = bar number 0-47 and M is sub-bar character A-D) + int ibar = 10*((int)volname(1,1)(0)-48)+(int)volname(2,1)(0)-48; // this is nasty, but it works + if (volname(0,1)(0) == 'B' && ibar >= 0 && ibar < 48) { - // radiator volume - if (touch->GetVolume()->GetName() == "QZBL") { if (trackinfo->GetGlueXHistory() == 0 && itrack > 0 && xin.dot(pin) > 0) { int pdgtype = track->GetDynamicParticle()->GetPDGcode(); int g3type = GlueXPrimaryGeneratorAction::ConvertPdgToGeant3(pdgtype); @@ -138,33 +153,156 @@ G4bool GlueXSensitiveDetectorDIRC::ProcessHits(G4Step* step, barhit.py_GeV = pin[1]/GeV; barhit.pz_GeV = pin[2]/GeV; barhit.pdg = g3type; - barhit.bar = touch_hist->GetReplicaNumber(0)/4; // each bar is glued from 4 pieces + barhit.bar = ibar; // from HDDS geometry barhit.track = itrack; // track id of the charged particle fHitsBar.push_back(barhit); - } return true; } + // wedge and mirrors volumes + if (volname == "FWM1" || volname == "FWM2" || volname == "FTMR" || + volname == "TSM1" || volname == "TSM2" || volname == "TSM3" || + volname == "FSM1" || volname == "FSM2" || volname == "OWDG" || + (volname(0,1)(0) == 'A' && volname(0,1)(1) == 'G') ) + { + + GlueXHitDIRCWob wobhit; + wobhit.track = track->GetTrackID(); + + // store normal to the closest boundary + G4int hNavId = G4ParallelWorldProcess::GetHypNavigatorID(); + std::vector::iterator iNav = + G4TransportationManager::GetTransportationManager()->GetActiveNavigatorsIterator(); + + G4bool valid; + G4ThreeVector localNormal = (iNav[hNavId])->GetLocalExitNormal(&valid); + if (valid){ + int mid=-1; + if (volname == "OWDG") { + if (localNormal.y()<-0.999) + mid=1; + else if (localNormal.y()>0.999) + mid=2; + else if(localNormal.z()>0.999) + mid=3; + else if(fabs(localNormal.z()+0.86)<0.01) + mid=4; + } + if (volname == "FSM1") + mid = 5; + if (volname == "FSM2") + mid = 6; + if (volname == "FWM1") + mid = 7; + if (volname == "FWM2") + mid = 8; + if (volname == "FTMR") + mid = 0; + if (volname == "TSM1") + mid=91; + if (volname == "TSM2") + mid=92; + if (volname == "TSM3") + mid=93; + if ((volname(0,1)(0) == 'A' && volname(0,1)(1) == 'G')) + mid=100; + + if (mid!=-1) { + G4double normalId = mid;// localNormal.x() + 10*localNormal.y() + 100*localNormal.z(); + wobhit.normalId = normalId; + fHitsWob.push_back(wobhit); + } + } + return true; + } + // PMT's pixel volume - if (touch->GetVolume()->GetName() == "PIXV"){ - if ((int)fHitsPmt.size() < MAX_HITS){ + if (volname == "PIXV") { + if ((int)fHitsPmt.size() < MAX_HITS) { + + // fix propagation speed for op + double tracklen=track->GetTrackLength()/cm; + double en=Ein/GeV; + double refindex= 1.43603+0.0132404*en-0.00225287*en*en+0.000500109*en*en*en; + double time_fixed=tracklen/(29.9792458/refindex); + +#ifdef DIRC_CHECK_PROPAGATION_TIME + double l_QZBL = 9.1 + 0.96 + 122.5; // 2 * (4*122.5) // 2 * bar length + wedge + window + double l_EPOTEK = 0.005 + 8 * 0.005; // window+wedge glue + 6 * bar joing glue + double l_AIR = 2 * 0.01; // air gap to mirror + double l_H2O = tracklen - l_QZBL - l_EPOTEK - l_AIR; + + // hard coded propagation time for 3.5 eV OpticalPhoton + double angle = 45/180*3.14159; + double time_propagated = l_QZBL/(29.9792458/1.476)/cos(angle); + time_propagated += l_EPOTEK/(29.9792458/1.616)/cos(angle); + time_propagated += l_AIR/(29.9792458)/cos(angle); + time_propagated += l_H2O/(29.9792458/1.343); + G4cout<<"Propagated time = "< 1) { G4cout << G4endl - << "--------> Hits Collection: in this event there are " - << fHitsBar.size() << " bar hits:" - << G4endl; + << "--------> Hits Collection: in this event there are " + << fHitsBar.size() << " bar hits:" + << G4endl; for(unsigned int h=0; h Hits Collection: in this event there are " - << fHitsPmt.size() << " PMT hits: " - << G4endl; + << "--------> Hits Collection: in this event there are " + << fHitsPmt.size() << " PMT hits: " + << G4endl; for(unsigned int h=0; hgetOutputRecord(); if (record == 0) { G4cerr << "GlueXSensitiveDetectorDIRC::EndOfEvent error - " - << "hits seen but no output hddm record to save them into, " - << "cannot continue!" << G4endl; + << "hits seen but no output hddm record to save them into, " + << "cannot continue!" << G4endl; exit(1); } @@ -228,7 +370,7 @@ void GlueXSensitiveDetectorDIRC::EndOfEvent(G4HCofThisEvent*) bhit(0).setTrack(fHitsBar[h].track); } - // Collect and output the DircTruthPoints + // Collect and output the DircTruthPmtHit for(unsigned int h=0; h @@ -44,11 +44,14 @@ class GlueXSensitiveDetectorDIRC : public G4VSensitiveDetector private: std::vector fHitsBar; + std::vector fHitsWob; std::vector fHitsPmt; + int fLutId; std::map fVolumeTable; static int MAX_HITS; + static int MAX_PIXELS; // put all other detector response parameters here static double TWO_HIT_TIME_RESOL; diff --git a/src/GlueXStackingAction.cc b/src/GlueXStackingAction.cc index 990678b..b348cf7 100644 --- a/src/GlueXStackingAction.cc +++ b/src/GlueXStackingAction.cc @@ -13,6 +13,15 @@ #include "G4ios.hh" #include "G4ParticleTable.hh" #include "TMath.h" +#include "G4TransportationManager.hh" +#include "G4PhysicalVolumeStore.hh" + +#ifdef DIRC_MONITORING_HISTOS +#include "TH1.h" +#include "TCanvas.h" +TH1F *hbouncez = new TH1F("hbouncez",";bounces along z [#];entries [#]",1000,0,2000); +TH1F *hbouncey = new TH1F("hbouncey",";bounces along y [#];entries [#]",1000,0,2000); +#endif GlueXStackingAction::GlueXStackingAction() { @@ -43,10 +52,28 @@ GlueXStackingAction::GlueXStackingAction() if (user_opts->Find("NOSECONDARIES", opt)) { nosecondaries = (opt[1] != 0); } + + fBarEnd = 0; + G4PhysicalVolumeStore* pvStore = G4PhysicalVolumeStore::GetInstance(); + for (size_t i=0; isize(); i++) { + if ((*pvStore)[i]->GetName()=="QZWN") { + fBarEnd = fabs((*pvStore)[i]->GetTranslation().x()); + } + } } GlueXStackingAction::~GlueXStackingAction() -{} +{ +#ifdef DIRC_MONITORING_HISTOS + TCanvas *c = new TCanvas("c","c",800,400); + hbouncez->Draw(); + c->Print("cbounces_z.png"); + c->Print("cbounces_z.C"); + hbouncey->Draw(); + c->Print("cbounces_y.png"); + c->Print("cbounces_y.C"); +#endif +} G4ClassificationOfNewTrack GlueXStackingAction::ClassifyNewTrack( const G4Track *aTrack) @@ -75,15 +102,60 @@ G4ClassificationOfNewTrack GlueXStackingAction::ClassifyNewTrack( // apply detection efficiency for the DIRC at production stage: G4String ParticleName = aTrack->GetDynamicParticle()->GetParticleDefinition()->GetParticleName(); - if (aTrack->GetParentID() > 0) { // particle is secondary + if (aTrack->GetParentID() != 0) { // particle is secondary if (ParticleName == "opticalphoton") { Double_t Ephoton = aTrack->GetMomentum().mag(); Double_t ra = G4UniformRand(); if (ra > GlueXSensitiveDetectorDIRC::GetDetectionEfficiency(Ephoton)) - return fKill; - } - } - + return fKill; + + G4ThreeVector v = aTrack->GetPosition(); + G4ThreeVector n = aTrack->GetMomentumDirection().unit(); + + Double_t bary = 35; // bar width + Double_t barz = 17.25; // bar height + Double_t barx = 4*1225; // bar length + + Double_t lenx; + if (v.y()<0) { + lenx = fabs(fBarEnd+v.x()); + if (n.x()>0) + lenx = 2*barx - lenx; + } + else { + lenx = fabs(v.x()-fBarEnd); + if (n.x()<0) + lenx = 2*barx - lenx; + } + + Double_t lenz = lenx*n.z()/fabs(n.x()); + Double_t leny = lenx*n.y()/fabs(n.x()); + + int bouncesz = fabs(lenz/barz); + int bouncesy = fabs(leny/bary); + +#ifdef DIRC_MONITORING_HISTOS + hbouncez->Fill(bouncesz); + hbouncey->Fill(bouncesy); +#endif + + Double_t anglez = fabs(n.getTheta()-CLHEP::pi/2.); + Double_t angley = fabs(n.angle(G4ThreeVector(0,1,0))-CLHEP::pi/2.); + + Double_t lambda = 197.0*2.0*CLHEP::pi/(aTrack->GetMomentum().mag()*1.0E6); + Double_t lambda2 = lambda*lambda; + + // calculate bounce probability + Double_t n_quartz = sqrt(1 + (0.696*lambda2/(lambda2-pow(0.068,2))) + (0.407*lambda2/(lambda2-pow(0.116,2))) + 0.897*lambda2/(lambda2-pow(9.896,2))); + Double_t bounce_probz = 1 - pow(4*CLHEP::pi*cos(anglez)*0.5*n_quartz/lambda,2);// 0.5 [nm] - roughness + Double_t bounce_proby = 1 - pow(4*CLHEP::pi*cos(angley)*0.5*n_quartz/lambda,2); + Double_t prob = pow(bounce_probz,bouncesz) * pow(bounce_proby,bouncesy); + + // transport efficiency + if (G4UniformRand() > prob) + return fKill; + } //else {return fKill;} // remove this condition!!! + } // if particle is secondary return fUrgent; } diff --git a/src/GlueXStackingAction.hh b/src/GlueXStackingAction.hh index 9dd21fb..5e3da41 100644 --- a/src/GlueXStackingAction.hh +++ b/src/GlueXStackingAction.hh @@ -33,8 +33,7 @@ class GlueXStackingAction : public G4UserStackingAction int nosecondaries; private: - TRandom *fRand; - TGraph* fDetEff; + double fBarEnd; }; #endif