From 7aa9b805fdd7e2adf82e353f909329da58f5c12a Mon Sep 17 00:00:00 2001 From: Alexander Somov Date: Tue, 23 Apr 2024 18:50:04 -0400 Subject: [PATCH 1/2] Initial ECAL libs --- src/GlueXDetectorConstruction.cc | 9 + src/GlueXHitECALblock.cc | 95 +++++++++ src/GlueXHitECALblock.hh | 66 ++++++ src/GlueXHitECALpoint.cc | 86 ++++++++ src/GlueXHitECALpoint.hh | 64 ++++++ src/GlueXSensitiveDetectorECAL.cc | 344 ++++++++++++++++++++++++++++++ src/GlueXSensitiveDetectorECAL.hh | 54 +++++ 7 files changed, 718 insertions(+) create mode 100644 src/GlueXHitECALblock.cc create mode 100644 src/GlueXHitECALblock.hh create mode 100644 src/GlueXHitECALpoint.cc create mode 100644 src/GlueXHitECALpoint.hh create mode 100644 src/GlueXSensitiveDetectorECAL.cc create mode 100644 src/GlueXSensitiveDetectorECAL.hh diff --git a/src/GlueXDetectorConstruction.cc b/src/GlueXDetectorConstruction.cc index a2848e8..2aa2c5d 100644 --- a/src/GlueXDetectorConstruction.cc +++ b/src/GlueXDetectorConstruction.cc @@ -17,6 +17,7 @@ #include "GlueXSensitiveDetectorFCALinsert.hh" #include "GlueXSensitiveDetectorGCAL.hh" #include "GlueXSensitiveDetectorCCAL.hh" +#include "GlueXSensitiveDetectorECAL.hh" #include "GlueXSensitiveDetectorFTOF.hh" #include "GlueXSensitiveDetectorDIRC.hh" #include "GlueXSensitiveDetectorCERE.hh" @@ -280,6 +281,7 @@ void GlueXDetectorConstruction::ConstructSDandField() GlueXSensitiveDetectorFCALinsert* fcalInsertHandler = 0; GlueXSensitiveDetectorGCAL* gcalHandler = 0; GlueXSensitiveDetectorCCAL* ccalHandler = 0; + GlueXSensitiveDetectorECAL* ecalHandler = 0; GlueXSensitiveDetectorFTOF* ftofHandler = 0; GlueXSensitiveDetectorDIRC* dircHandler = 0; GlueXSensitiveDetectorCERE* cereHandler = 0; @@ -377,6 +379,13 @@ void GlueXDetectorConstruction::ConstructSDandField() } iter->second->SetSensitiveDetector(ccalHandler); } + else if (volname == "XTBL") { + if (ecalHandler == 0) { + ecalHandler = new GlueXSensitiveDetectorECAL("ecal"); + SDman->AddNewDetector(ecalHandler); + } + iter->second->SetSensitiveDetector(ecalHandler); + } else if (volname == "FTOC" || volname == "FTOX" || volname == "FTOH" || volname == "FTOL") { diff --git a/src/GlueXHitECALblock.cc b/src/GlueXHitECALblock.cc new file mode 100644 index 0000000..de31053 --- /dev/null +++ b/src/GlueXHitECALblock.cc @@ -0,0 +1,95 @@ +// +// GlueXHitECALblock - class implementation +// + +#include "GlueXHitECALblock.hh" + +G4ThreadLocal G4Allocator* GlueXHitECALblockAllocator = 0; + +GlueXHitECALblock::GlueXHitECALblock(G4int column, G4int row) + : G4VHit(), + column_(column), + row_(row) +{} + +GlueXHitECALblock::GlueXHitECALblock(const GlueXHitECALblock &src) +{ + column_ = src.column_; + row_ = src.row_; + hits = src.hits; +} + +int GlueXHitECALblock::operator==(const GlueXHitECALblock &right) const +{ + if (column_ != right.column_ || row_ != right.row_) { + return 0; + } + else if (hits.size() != right.hits.size()) { + return 0; + } + + for (int ih=0; ih < (int)hits.size(); ++ih) { + if (hits[ih].E_GeV != right.hits[ih].E_GeV || + hits[ih].t_ns != right.hits[ih].t_ns || + hits[ih].dE_lightguide_GeV != right.hits[ih].dE_lightguide_GeV || + hits[ih].t_lightguide_ns != right.hits[ih].t_lightguide_ns) + { + return 0; + } + } + return 1; +} + +GlueXHitECALblock &GlueXHitECALblock::operator+=(const GlueXHitECALblock &right) +{ + if (column_ != right.column_ || row_ != right.row_) { + G4cerr << "Error in GlueXHitECALblock::operator+=() - " + << "illegal attempt to merge hits from two different blocks!" + << G4endl; + return *this; + } + std::vector::iterator hiter = hits.begin(); + std::vector::const_iterator hitsrc; + for (hitsrc = right.hits.begin(); hitsrc != right.hits.end(); ++hitsrc) { + for (; hiter != hits.end(); ++hiter) { + if (hiter->t_ns > hitsrc->t_ns) + break; + } + hiter = hits.insert(hiter, *hitsrc); + } + return *this; +} + +void GlueXHitECALblock::Draw() const +{ + // not yet implemented +} + +void GlueXHitECALblock::Print() const +{ + G4cout << "GlueXHitECALblock: " + << " column = " << column_ + << ", row = " << row_ << G4endl; + std::vector::const_iterator hiter; + for (hiter = hits.begin(); hiter != hits.end(); ++hiter) { + G4cout << " E = " << hiter->E_GeV << " GeV" << G4endl + << " t = " << hiter->t_ns << " ns" << G4endl + << " E(lightguide) = " + << hiter->dE_lightguide_GeV << " GeV" << G4endl + << " t(lightguide) = " + << hiter->t_lightguide_ns << " ns" << G4endl + << G4endl; + } +} + +void printallhits(GlueXHitsMapECALblock *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/GlueXHitECALblock.hh b/src/GlueXHitECALblock.hh new file mode 100644 index 0000000..b5ffa8f --- /dev/null +++ b/src/GlueXHitECALblock.hh @@ -0,0 +1,66 @@ +// +// GlueXHitECALblock - class header +// +// 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 GlueXHitECALblock_h +#define GlueXHitECALblock_h 1 + +#include "G4VHit.hh" +#include "G4THitsMap.hh" +#include "G4Allocator.hh" + +class GlueXHitECALblock : public G4VHit +{ + public: + GlueXHitECALblock() {} + GlueXHitECALblock(G4int column, G4int row); + GlueXHitECALblock(const GlueXHitECALblock &src); + int operator==(const GlueXHitECALblock &right) const; + GlueXHitECALblock &operator+=(const GlueXHitECALblock &right); + + void *operator new(size_t); + void operator delete(void *aHit); + + void Draw() const; + void Print() const; + + // no reason to hide hit data + + G4int column_; // ECAL block column, from 1 increasing x + G4int row_; // ECAL block row, from 1 increasing y + + struct hitinfo_t { + G4double E_GeV; // energy deposition (GeV) + G4double t_ns; // pulse leading-edge time (ns) + G4double dE_lightguide_GeV; // light guide energy deposition (GeV) + G4double t_lightguide_ns; // light guide pulse leading-edge time (ns) + }; + std::vector hits; + + G4int GetKey() const { return GetKey(column_, row_); } + static G4int GetKey(G4int column, G4int row) { + return ((row + 1) << 16) + column + 1; + } +}; + +typedef G4THitsMap GlueXHitsMapECALblock; + +extern G4ThreadLocal G4Allocator* GlueXHitECALblockAllocator; + +inline void* GlueXHitECALblock::operator new(size_t) +{ + if (!GlueXHitECALblockAllocator) + GlueXHitECALblockAllocator = new G4Allocator; + return (void *) GlueXHitECALblockAllocator->MallocSingle(); +} + +inline void GlueXHitECALblock::operator delete(void *aHit) +{ + GlueXHitECALblockAllocator->FreeSingle((GlueXHitECALblock*) aHit); +} + +#endif diff --git a/src/GlueXHitECALpoint.cc b/src/GlueXHitECALpoint.cc new file mode 100644 index 0000000..954dd6f --- /dev/null +++ b/src/GlueXHitECALpoint.cc @@ -0,0 +1,86 @@ +// +// GlueXHitECALpoint - class implementation +// + +#include "GlueXHitECALpoint.hh" + +G4ThreadLocal G4Allocator* GlueXHitECALpointAllocator = 0; + +GlueXHitECALpoint::GlueXHitECALpoint(const GlueXHitECALpoint &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 GlueXHitECALpoint::operator==(const GlueXHitECALpoint &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; +} + +GlueXHitECALpoint &GlueXHitECALpoint::operator+=(const GlueXHitECALpoint &right) +{ + G4cerr << "Error in GlueXHitECALpoint::operator+= - " + << "illegal attempt to merge two TruthShower objects in the ecal!" + << G4endl; + return *this; +} + +void GlueXHitECALpoint::Draw() const +{ + // not yet implemented +} + +void GlueXHitECALpoint::Print() const +{ + G4cout << "GlueXHitECALpoint:" << 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(GlueXHitsMapECALpoint *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/GlueXHitECALpoint.hh b/src/GlueXHitECALpoint.hh new file mode 100644 index 0000000..24e2d07 --- /dev/null +++ b/src/GlueXHitECALpoint.hh @@ -0,0 +1,64 @@ +// +// GlueXHitECALpoint - class header +// +// 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 GlueXHitECALpoint_h +#define GlueXHitECALpoint_h 1 + +#include "G4VHit.hh" +#include "G4THitsMap.hh" +#include "G4Allocator.hh" + +class GlueXHitECALpoint : public G4VHit +{ + public: + GlueXHitECALpoint() {} + GlueXHitECALpoint(const GlueXHitECALpoint &src); + int operator==(const GlueXHitECALpoint &right) const; + GlueXHitECALpoint &operator+=(const GlueXHitECALpoint &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 shower particle at this point + G4bool primary_; // true if track belongs to from a primary particle + G4int ptype_G3; // G3 type of particle making this shower + G4double px_GeV; // momentum (GeV/c) of shower particle at point, x component + G4double py_GeV; // momentum (GeV/c) of shower particle at point, y component + G4double pz_GeV; // momentum (GeV/c) of shower particle at point, z component + G4double x_cm; // global x coordinate of shower particle at point (cm) + G4double y_cm; // global y coordinate of shower particle at point (cm) + G4double z_cm; // global z coordinate of shower particle at point (cm) + G4double t_ns; // time of shower particle crossing at point (ns) + G4int track_; // Geant4 track ID of particle making this shower + G4int trackID_; // GlueX-assigned track ID of particle making this shower + + G4int GetKey() const { return (track_ << 20) + int(t_ns * 100); } +}; + +typedef G4THitsMap GlueXHitsMapECALpoint; + +extern G4ThreadLocal G4Allocator* GlueXHitECALpointAllocator; + +inline void* GlueXHitECALpoint::operator new(size_t) +{ + if (!GlueXHitECALpointAllocator) + GlueXHitECALpointAllocator = new G4Allocator; + return (void *) GlueXHitECALpointAllocator->MallocSingle(); +} + +inline void GlueXHitECALpoint::operator delete(void *aHit) +{ + GlueXHitECALpointAllocator->FreeSingle((GlueXHitECALpoint*) aHit); +} + +#endif diff --git a/src/GlueXSensitiveDetectorECAL.cc b/src/GlueXSensitiveDetectorECAL.cc new file mode 100644 index 0000000..2cba17f --- /dev/null +++ b/src/GlueXSensitiveDetectorECAL.cc @@ -0,0 +1,344 @@ +// +// GlueXSensitiveDetectorECAL - class implementation +// + +#include "GlueXSensitiveDetectorECAL.hh" +#include "GlueXDetectorConstruction.hh" +#include "GlueXPrimaryGeneratorAction.hh" +#include "GlueXUserEventInformation.hh" +#include "GlueXUserTrackInformation.hh" +#include "HddmOutput.hh" + +#include "G4VPhysicalVolume.hh" +#include "G4PVPlacement.hh" +#include "G4EventManager.hh" +#include "G4HCofThisEvent.hh" +#include "G4Step.hh" +#include "G4SDManager.hh" +#include "G4ios.hh" + +#include + +// Cutoff on the number of allowed hits per block +int GlueXSensitiveDetectorECAL::MAX_HITS = 100; + +// Geometry constants for the ECAL +double GlueXSensitiveDetectorECAL::WIDTH_OF_BLOCK = 2.09*cm; +double GlueXSensitiveDetectorECAL::LENGTH_OF_BLOCK = 20.0*cm; +int GlueXSensitiveDetectorECAL::CENTRAL_COLUMN = 24; +int GlueXSensitiveDetectorECAL::CENTRAL_ROW = 24; + +// Light propagation parameters in forward calorimeter +double GlueXSensitiveDetectorECAL::ATTENUATION_LENGTH = 100.*cm; +double GlueXSensitiveDetectorECAL::C_EFFECTIVE = 15.*cm/ns; + +// Minimum hit time difference for two hits on the same block +double GlueXSensitiveDetectorECAL::TWO_HIT_TIME_RESOL = 75.*ns; + +// Minimum energy deposition for a hit +double GlueXSensitiveDetectorECAL::THRESH_MEV = 5.; + +int GlueXSensitiveDetectorECAL::instanceCount = 0; +G4Mutex GlueXSensitiveDetectorECAL::fMutex = G4MUTEX_INITIALIZER; + +GlueXSensitiveDetectorECAL::GlueXSensitiveDetectorECAL(const G4String& name) + : G4VSensitiveDetector(name), + fBlocksMap(0), fPointsMap(0) +{ + collectionName.insert("ECALBlockHitsCollection"); + collectionName.insert("ECALPointsCollection"); + + // The rest of this only needs to happen once, the first time an object + // of this type is instantiated for this configuration of geometry and + // fields. If the geometry or fields change in such a way as to modify + // the drift-time properties of hits in the ECAL, you must delete all old + // objects of this class and create new ones. + + G4AutoLock barrier(&fMutex); + if (instanceCount++ == 0) { + int runno = HddmOutput::getRunNo(); + extern jana::JApplication *japp; + if (japp == 0) { + G4cerr << "Error in GlueXSensitiveDetector constructor - " + << "jana global DApplication object not set, " + << "cannot continue." << G4endl; + exit(-1); + } + jana::JCalibration *jcalib = japp->GetJCalibration(runno); + std::map fcal_parms; + jcalib->Get("FCAL/fcal_parms", fcal_parms); + ATTENUATION_LENGTH = fcal_parms.at("FCAL_ATTEN_LENGTH")*cm; + C_EFFECTIVE = fcal_parms.at("FCAL_C_EFFECTIVE")*cm/ns; + TWO_HIT_TIME_RESOL = fcal_parms.at("FCAL_TWO_HIT_RESOL")*ns; + MAX_HITS = fcal_parms.at("FCAL_MAX_HITS"); + THRESH_MEV = fcal_parms.at("FCAL_THRESH_MEV"); + + G4cout << "ECAL: ALL parameters loaded from ccdb" << G4endl; + } +} + +GlueXSensitiveDetectorECAL::GlueXSensitiveDetectorECAL( + const GlueXSensitiveDetectorECAL &src) + : G4VSensitiveDetector(src), + fBlocksMap(src.fBlocksMap), fPointsMap(src.fPointsMap) +{ + G4AutoLock barrier(&fMutex); + ++instanceCount; +} + +GlueXSensitiveDetectorECAL &GlueXSensitiveDetectorECAL::operator=(const + GlueXSensitiveDetectorECAL &src) +{ + G4AutoLock barrier(&fMutex); + *(G4VSensitiveDetector*)this = src; + fBlocksMap = src.fBlocksMap; + fPointsMap = src.fPointsMap; + return *this; +} + +GlueXSensitiveDetectorECAL::~GlueXSensitiveDetectorECAL() +{ + G4AutoLock barrier(&fMutex); + --instanceCount; +} + +void GlueXSensitiveDetectorECAL::Initialize(G4HCofThisEvent* hce) +{ + fBlocksMap = new + GlueXHitsMapECALblock(SensitiveDetectorName, collectionName[0]); + fPointsMap = new + GlueXHitsMapECALpoint(SensitiveDetectorName, collectionName[1]); + G4SDManager *sdm = G4SDManager::GetSDMpointer(); + hce->AddHitsCollection(sdm->GetCollectionID(collectionName[0]), fBlocksMap); + hce->AddHitsCollection(sdm->GetCollectionID(collectionName[1]), fPointsMap); +} + +G4bool GlueXSensitiveDetectorECAL::ProcessHits(G4Step* step, + G4TouchableHistory* ROhist) +{ + double dEsum = step->GetTotalEnergyDeposit(); + const G4ThreeVector &pin = step->GetPreStepPoint()->GetMomentum(); + const G4ThreeVector &xin = step->GetPreStepPoint()->GetPosition(); + const G4ThreeVector &xout = step->GetPostStepPoint()->GetPosition(); + double Ein = step->GetPreStepPoint()->GetTotalEnergy(); + double tin = step->GetPreStepPoint()->GetGlobalTime(); + double tout = step->GetPostStepPoint()->GetGlobalTime(); + G4ThreeVector x = (xin + xout) / 2; + double t = (tin + tout) / 2; + + const G4VTouchable* touch = step->GetPreStepPoint()->GetTouchable(); + const G4AffineTransform &local_from_global = touch->GetHistory() + ->GetTopTransform(); + G4ThreeVector xlocal = local_from_global.TransformPoint(x); + + // For particles that range out inside the active volume, the + // "out" time may sometimes be set to something enormously high. + // This screws up the hit. Check for this case here by looking + // at tout and making sure it is less than 1 second. If it's + // not, then just use tin for "t". + + if (tout > 1.0*s) + t = tin; + + // Post the hit to the points list in the + // order of appearance in the event simulation. + + G4Track *track = step->GetTrack(); + G4int trackID = track->GetTrackID(); + GlueXUserTrackInformation *trackinfo = (GlueXUserTrackInformation*) + track->GetUserInformation(); + int itrack = trackinfo->GetGlueXTrackID(); + if (trackinfo->GetGlueXHistory() == 0 && + xin.dot(pin) > 0 && Ein/MeV > THRESH_MEV) + { + int pdgtype = track->GetDynamicParticle()->GetPDGcode(); + int g3type = GlueXPrimaryGeneratorAction::ConvertPdgToGeant3(pdgtype); + GlueXHitECALpoint newPoint; + newPoint.ptype_G3 = g3type; + newPoint.track_ = trackID; + newPoint.trackID_ = itrack; + newPoint.primary_ = (track->GetParentID() == 0); + newPoint.t_ns = t/ns; + newPoint.x_cm = xin[0]/cm; + newPoint.y_cm = xin[1]/cm; + newPoint.z_cm = xin[2]/cm; + newPoint.px_GeV = pin[0]/GeV; + newPoint.py_GeV = pin[1]/GeV; + newPoint.pz_GeV = pin[2]/GeV; + newPoint.E_GeV = Ein/GeV; + G4int key = fPointsMap->entries(); + fPointsMap->add(key, newPoint); + trackinfo->SetGlueXHistory(2); + } + + // Post the hit to the hits map, ordered by sector index + + if (dEsum > 0) { + int column = GetIdent("column", touch); + int row = GetIdent("row", touch); + int key = GlueXHitECALblock::GetKey(column, row); + GlueXHitECALblock *block = (*fBlocksMap)[key]; + if (block == 0) { + GlueXHitECALblock newblock(column, row); + fBlocksMap->add(key, newblock); + block = (*fBlocksMap)[key]; + } + + // Handle hits in the lead tungtate crystal + + if (touch->GetVolume()->GetName() == "XTBL") { + double dist = 0.5 * LENGTH_OF_BLOCK - xlocal[2]; + double dEcorr = dEsum * exp(-dist / ATTENUATION_LENGTH); + double tcorr = t + dist / C_EFFECTIVE; + + // Add the hit to the hits vector, maintaining strict time ordering + + int merge_hit = 0; + std::vector::iterator hiter; + for (hiter = block->hits.begin(); hiter != block->hits.end(); ++hiter) { + if (fabs(hiter->t_ns*ns - tcorr) < TWO_HIT_TIME_RESOL) { + merge_hit = 1; + break; + } + else if (hiter->t_ns*ns > tcorr) { + break; + } + } + if (merge_hit) { + // Merge the time with the existing hit, add the energy deposition + hiter->t_ns = hiter->t_ns * hiter->E_GeV + dEcorr/GeV * tcorr/ns; + hiter->E_GeV += dEcorr/GeV; + hiter->t_ns /= hiter->E_GeV; + } + else { + // create new hit + hiter = block->hits.insert(hiter, GlueXHitECALblock::hitinfo_t()); + hiter->E_GeV = dEcorr/GeV; + hiter->t_ns = tcorr/ns; + hiter->dE_lightguide_GeV = 0; + hiter->t_lightguide_ns = 0; + } + } + } + return true; +} + +void GlueXSensitiveDetectorECAL::EndOfEvent(G4HCofThisEvent*) +{ + std::map *blocks = fBlocksMap->GetMap(); + std::map *points = fPointsMap->GetMap(); + if (blocks->size() == 0 && points->size() == 0) + return; + std::map::iterator biter; + std::map::iterator piter; + + if (verboseLevel > 1) { + G4cout << G4endl + << "--------> Hits Collection: in this event there are " + << blocks->size() << " blocks with hits in the ECAL: " + << G4endl; + for (biter = blocks->begin(); biter != blocks->end(); ++biter) + biter->second->Print(); + + G4cout << G4endl + << "--------> Hits Collection: in this event there are " + << points->size() << " truth showers in the ECAL: " + << G4endl; + for (piter = points->begin(); piter != points->end(); ++piter) + piter->second->Print(); + } + + // pack hits into ouptut hddm record + + G4EventManager* mgr = G4EventManager::GetEventManager(); + G4VUserEventInformation* info = mgr->GetUserInformation(); + hddm_s::HDDM *record = ((GlueXUserEventInformation*)info)->getOutputRecord(); + if (record == 0) { + G4cerr << "GlueXSensitiveDetectorECAL::EndOfEvent error - " + << "hits seen but no output hddm record to save them into, " + << "cannot continue!" << G4endl; + exit(1); + } + + if (record->getPhysicsEvents().size() == 0) + record->addPhysicsEvents(); + if (record->getHitViews().size() == 0) + record->getPhysicsEvent().addHitViews(); + hddm_s::HitView &hitview = record->getPhysicsEvent().getHitView(); + if (hitview.getCrystalEcals().size() == 0) + hitview.addCrystalEcals(); + + hddm_s::CrystalEcal &CrystalEcal = hitview.getCrystalEcal(); + + // Collect and output the ecalTruthHits + for (biter = blocks->begin(); biter != blocks->end(); ++biter) { + std::vector &hits = biter->second->hits; + // apply a pulse height threshold cut + for (unsigned int ih=0; ih < hits.size(); ++ih) { + if (hits[ih].E_GeV <= THRESH_MEV/1e3) { + hits.erase(hits.begin() + ih); + --ih; + } + } + + hddm_s::EcalBlockList block = CrystalEcal.addEcalBlocks(1); + block(0).setColumn(biter->second->column_); + block(0).setRow(biter->second->row_); + int hitscount = hits.size(); + if (hitscount > MAX_HITS) { + hitscount = MAX_HITS; + G4cerr << "GlueXSensitiveDetectorECAL::EndOfEvent warning:" + << "max hit count " << MAX_HITS << " exceeded, " + << hits.size() - hitscount << " hits discarded." + << G4endl; + } + for (int ih=0; ih < hitscount; ++ih) { + hddm_s::EcalTruthHitList thit = block(0).addEcalTruthHits(1); + thit(0).setE(hits[ih].E_GeV); + thit(0).setT(hits[ih].t_ns); + } + } + + // Collect and output the ecalTruthShowers + for (piter = points->begin(); piter != points->end(); ++piter) { + hddm_s::EcalTruthShowerList point = CrystalEcal.addEcalTruthShowers(1); + point(0).setE(piter->second->E_GeV); + point(0).setPrimary(piter->second->primary_); + point(0).setPtype(piter->second->ptype_G3); + point(0).setPx(piter->second->px_GeV); + point(0).setPy(piter->second->py_GeV); + point(0).setPz(piter->second->pz_GeV); + point(0).setX(piter->second->x_cm); + point(0).setY(piter->second->y_cm); + point(0).setZ(piter->second->z_cm); + point(0).setT(piter->second->t_ns); + point(0).setTrack(piter->second->track_); + hddm_s::TrackIDList tid = point(0).addTrackIDs(); + tid(0).setItrack(piter->second->trackID_); + } +} + +int GlueXSensitiveDetectorECAL::GetIdent(std::string div, + const G4VTouchable *touch) +{ + const HddsG4Builder* bldr = GlueXDetectorConstruction::GetBuilder(); + std::map >::const_iterator iter; + std::map > *identifiers; + int max_depth = touch->GetHistoryDepth(); + for (int depth = 0; depth < max_depth; ++depth) { + G4VPhysicalVolume *pvol = touch->GetVolume(depth); + G4LogicalVolume *lvol = pvol->GetLogicalVolume(); + int volId = fVolumeTable[lvol]; + if (volId == 0) { + volId = bldr->getVolumeId(lvol); + fVolumeTable[lvol] = volId; + } + identifiers = &Refsys::fIdentifierTable[volId]; + if ((iter = identifiers->find(div)) != identifiers->end()) { + int copyNum = touch->GetCopyNumber(depth); + copyNum += (dynamic_cast(pvol))? -1 : 0; + return iter->second[copyNum]; + } + } + return -1; +} diff --git a/src/GlueXSensitiveDetectorECAL.hh b/src/GlueXSensitiveDetectorECAL.hh new file mode 100644 index 0000000..4cda004 --- /dev/null +++ b/src/GlueXSensitiveDetectorECAL.hh @@ -0,0 +1,54 @@ +// +// GlueXSensitiveDetectorECAL - class header +// +// 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. + +#ifndef GlueXSensitiveDetectorECAL_h +#define GlueXSensitiveDetectorECAL_h 1 + +#include "G4VSensitiveDetector.hh" +#include "G4AutoLock.hh" + +#include "GlueXHitECALblock.hh" +#include "GlueXHitECALpoint.hh" + +class G4Step; +class G4HCofThisEvent; + +class GlueXSensitiveDetectorECAL : public G4VSensitiveDetector +{ + public: + GlueXSensitiveDetectorECAL(const G4String& name); + GlueXSensitiveDetectorECAL(const GlueXSensitiveDetectorECAL &right); + GlueXSensitiveDetectorECAL &operator=(const GlueXSensitiveDetectorECAL &right); + virtual ~GlueXSensitiveDetectorECAL(); + + virtual void Initialize(G4HCofThisEvent* hitCollection); + virtual G4bool ProcessHits(G4Step* step, G4TouchableHistory* ROhist); + virtual void EndOfEvent(G4HCofThisEvent* hitCollection); + + int GetIdent(std::string div, const G4VTouchable *touch); + + private: + GlueXHitsMapECALblock* fBlocksMap; + GlueXHitsMapECALpoint* fPointsMap; + + std::map fVolumeTable; + + static int MAX_HITS; + static int CENTRAL_COLUMN; + static int CENTRAL_ROW; + static double WIDTH_OF_BLOCK; + static double LENGTH_OF_BLOCK; + static double ATTENUATION_LENGTH; + static double C_EFFECTIVE; + static double TWO_HIT_TIME_RESOL; + static double THRESH_MEV; + + static int instanceCount; + static G4Mutex fMutex; +}; + +#endif From 2a18fc240b0b4ed96d034188daf0af287ad6b114 Mon Sep 17 00:00:00 2001 From: Alexander Somov Date: Tue, 23 Apr 2024 18:50:43 -0400 Subject: [PATCH 2/2] Initial ECAL libs --- GNUmakefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GNUmakefile b/GNUmakefile index 1a1e302..481929f 100644 --- a/GNUmakefile +++ b/GNUmakefile @@ -87,7 +87,7 @@ else endif DANALIBS = -L$(HALLD_RECON_HOME)/$(BMS_OSNAME)/lib -lHDGEOMETRY -lDANA \ - -lANALYSIS -lBCAL -lCCAL -lCDC -lCERE -lTRD -lDIRC -lFCAL \ + -lANALYSIS -lBCAL -lCCAL -lECAL -lCDC -lCERE -lTRD -lDIRC -lFCAL \ -lFDC -lFMWPC -lHDDM -lPAIR_SPECTROMETER -lPID -lRF \ -lSTART_COUNTER -lTAGGER -lTOF -lTPOL -lTRACKING \ -lTRIGGER -lDAQ -lTTAB -lEVENTSTORE -lKINFITTER -lTAC \