1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | #include "GlueXSensitiveDetectorCCAL.hh" |
8 | #include "GlueXDetectorConstruction.hh" |
9 | #include "GlueXPrimaryGeneratorAction.hh" |
10 | #include "GlueXUserEventInformation.hh" |
11 | #include "GlueXUserTrackInformation.hh" |
12 | #include "HddmOutput.hh" |
13 | |
14 | #include "G4VPhysicalVolume.hh" |
15 | #include "G4PVPlacement.hh" |
16 | #include "G4EventManager.hh" |
17 | #include "G4HCofThisEvent.hh" |
18 | #include "G4Step.hh" |
19 | #include "G4SDManager.hh" |
20 | #include "G4ios.hh" |
21 | |
22 | #include <JANA/JApplication.h> |
23 | |
24 | |
25 | int GlueXSensitiveDetectorCCAL::MAX_HITS = 100; |
26 | |
27 | |
28 | |
29 | |
30 | double GlueXSensitiveDetectorCCAL::WIDTH_OF_BLOCK = 2.0*cm; |
31 | double GlueXSensitiveDetectorCCAL::LENGTH_OF_BLOCK = 20.0*cm; |
32 | |
33 | |
34 | double GlueXSensitiveDetectorCCAL::ATTENUATION_LENGTH = 200.*cm; |
35 | double GlueXSensitiveDetectorCCAL::C_EFFECTIVE = 13.*cm/ns; |
36 | |
37 | |
38 | double GlueXSensitiveDetectorCCAL::TWO_HIT_TIME_RESOL = 75.*ns; |
39 | |
40 | |
41 | double GlueXSensitiveDetectorCCAL::THRESH_MEV = 5.; |
42 | |
43 | int GlueXSensitiveDetectorCCAL::instanceCount = 0; |
44 | G4Mutex GlueXSensitiveDetectorCCAL::fMutex = G4MUTEX_INITIALIZER{ { 0, 0, 0, 0, 0, 0, 0, { 0, 0 } } }; |
45 | |
46 | GlueXSensitiveDetectorCCAL::GlueXSensitiveDetectorCCAL(const G4String& name) |
47 | : G4VSensitiveDetector(name), |
48 | fBlocksMap(0), fPointsMap(0) |
49 | { |
50 | collectionName.insert("CCALBlockHitsCollection"); |
51 | collectionName.insert("CCALPointsCollection"); |
52 | |
53 | |
54 | |
55 | |
56 | |
57 | |
58 | |
59 | G4AutoLock barrier(&fMutex); |
60 | if (instanceCount++ == 0) { |
61 | int runno = HddmOutput::getRunNo(); |
62 | extern jana::JApplication *japp; |
63 | if (japp == 0) { |
64 | G4cerr(*G4cerr_p) << "Error in GlueXSensitiveDetector constructor - " |
65 | << "jana global DApplication object not set, " |
66 | << "cannot continue." << G4endlstd::endl; |
67 | exit(-1); |
68 | } |
69 | jana::JCalibration *jcalib = japp->GetJCalibration(runno); |
| Value stored to 'jcalib' during its initialization is never read |
70 | if (japp == 0) { |
71 | jcalib = 0; |
72 | G4cout(*G4cout_p) << "CCAL: ALL parameters loaded from ccdb" << G4endlstd::endl; |
73 | } |
74 | } |
75 | } |
76 | |
77 | GlueXSensitiveDetectorCCAL::GlueXSensitiveDetectorCCAL( |
78 | const GlueXSensitiveDetectorCCAL &src) |
79 | : G4VSensitiveDetector(src), |
80 | fBlocksMap(src.fBlocksMap), fPointsMap(src.fPointsMap) |
81 | { |
82 | G4AutoLock barrier(&fMutex); |
83 | ++instanceCount; |
84 | } |
85 | |
86 | GlueXSensitiveDetectorCCAL &GlueXSensitiveDetectorCCAL::operator=(const |
87 | GlueXSensitiveDetectorCCAL &src) |
88 | { |
89 | G4AutoLock barrier(&fMutex); |
90 | *(G4VSensitiveDetector*)this = src; |
91 | fBlocksMap = src.fBlocksMap; |
92 | fPointsMap = src.fPointsMap; |
93 | return *this; |
94 | } |
95 | |
96 | GlueXSensitiveDetectorCCAL::~GlueXSensitiveDetectorCCAL() |
97 | { |
98 | G4AutoLock barrier(&fMutex); |
99 | --instanceCount; |
100 | } |
101 | |
102 | void GlueXSensitiveDetectorCCAL::Initialize(G4HCofThisEvent* hce) |
103 | { |
104 | fBlocksMap = new |
105 | GlueXHitsMapCCALblock(SensitiveDetectorName, collectionName[0]); |
106 | fPointsMap = new |
107 | GlueXHitsMapCCALpoint(SensitiveDetectorName, collectionName[1]); |
108 | G4SDManager *sdm = G4SDManager::GetSDMpointer(); |
109 | hce->AddHitsCollection(sdm->GetCollectionID(collectionName[0]), fBlocksMap); |
110 | hce->AddHitsCollection(sdm->GetCollectionID(collectionName[1]), fPointsMap); |
111 | } |
112 | |
113 | G4bool GlueXSensitiveDetectorCCAL::ProcessHits(G4Step* step, |
114 | G4TouchableHistory* ROhist) |
115 | { |
116 | double dEsum = step->GetTotalEnergyDeposit(); |
117 | if (dEsum == 0) |
118 | return false; |
119 | |
120 | const G4ThreeVector &pin = step->GetPreStepPoint()->GetMomentum(); |
121 | const G4ThreeVector &xin = step->GetPreStepPoint()->GetPosition(); |
122 | const G4ThreeVector &xout = step->GetPostStepPoint()->GetPosition(); |
123 | double Ein = step->GetPreStepPoint()->GetTotalEnergy(); |
124 | double tin = step->GetPreStepPoint()->GetGlobalTime(); |
125 | double tout = step->GetPostStepPoint()->GetGlobalTime(); |
126 | G4ThreeVector x = (xin + xout) / 2; |
127 | double t = (tin + tout) / 2; |
128 | |
129 | const G4VTouchable* touch = step->GetPreStepPoint()->GetTouchable(); |
130 | const G4AffineTransform &local_from_global = touch->GetHistory() |
131 | ->GetTopTransform(); |
132 | G4ThreeVector xlocal = local_from_global.TransformPoint(x); |
133 | |
134 | |
135 | |
136 | |
137 | |
138 | |
139 | |
140 | if (tout > 1.0*s) |
141 | t = tin; |
142 | |
143 | |
144 | |
145 | |
146 | G4Track *track = step->GetTrack(); |
147 | G4int trackID = track->GetTrackID(); |
148 | GlueXUserTrackInformation *trackinfo = (GlueXUserTrackInformation*) |
149 | track->GetUserInformation(); |
150 | int itrack = trackinfo->GetGlueXTrackID(); |
151 | if (track->GetCurrentStepNumber() == 1) |
152 | trackinfo->SetGlueXHistory(4); |
153 | if (trackinfo->GetGlueXHistory() == 0 && |
154 | xin.dot(pin) > 0 && Ein/MeV > THRESH_MEV) |
155 | { |
156 | G4int key = fPointsMap->entries(); |
157 | GlueXHitCCALpoint* lastPoint = (*fPointsMap)[key - 1]; |
158 | |
159 | if (lastPoint == 0 || lastPoint->track_ != trackID) { |
160 | int pdgtype = track->GetDynamicParticle()->GetPDGcode(); |
161 | int g3type = GlueXPrimaryGeneratorAction::ConvertPdgToGeant3(pdgtype); |
162 | GlueXHitCCALpoint newPoint; |
163 | newPoint.ptype_G3 = g3type; |
164 | newPoint.track_ = trackID; |
165 | newPoint.trackID_ = itrack; |
166 | newPoint.primary_ = (track->GetParentID() == 0); |
167 | newPoint.t_ns = t/ns; |
168 | newPoint.x_cm = xin[0]/cm; |
169 | newPoint.y_cm = xin[1]/cm; |
170 | newPoint.z_cm = xin[2]/cm; |
171 | newPoint.px_GeV = pin[0]/GeV; |
172 | newPoint.py_GeV = pin[1]/GeV; |
173 | newPoint.pz_GeV = pin[2]/GeV; |
174 | newPoint.E_GeV = Ein/GeV; |
175 | fPointsMap->add(key, newPoint); |
176 | } |
177 | } |
178 | |
179 | |
180 | |
181 | if (dEsum > 0) { |
182 | int column = GetIdent("column", touch); |
183 | int row = GetIdent("row", touch); |
184 | int key = GlueXHitCCALblock::GetKey(column, row); |
185 | GlueXHitCCALblock *block = (*fBlocksMap)[key]; |
186 | if (block == 0) { |
187 | GlueXHitCCALblock newblock(column, row); |
188 | fBlocksMap->add(key, newblock); |
189 | block = (*fBlocksMap)[key]; |
190 | } |
191 | double dist = 0.5 * LENGTH_OF_BLOCK - xlocal[2]; |
192 | double dEcorr = dEsum * exp(-dist / ATTENUATION_LENGTH); |
193 | double tcorr = t + dist / C_EFFECTIVE; |
194 | |
195 | |
196 | |
197 | int merge_hit = 0; |
198 | std::vector<GlueXHitCCALblock::hitinfo_t>::iterator hiter; |
199 | for (hiter = block->hits.begin(); hiter != block->hits.end(); ++hiter) { |
200 | if (fabs(hiter->t_ns*ns - tcorr) < TWO_HIT_TIME_RESOL) { |
201 | merge_hit = 1; |
202 | break; |
203 | } |
204 | else if (hiter->t_ns*ns > t) { |
205 | break; |
206 | } |
207 | } |
208 | if (merge_hit) { |
209 | |
210 | hiter->E_GeV += dEcorr/GeV; |
211 | if (hiter->t_ns*ns > tcorr) { |
212 | hiter->t_ns = tcorr/ns; |
213 | } |
214 | } |
215 | else if ((int)block->hits.size() < MAX_HITS) { |
216 | |
217 | hiter = block->hits.insert(hiter, GlueXHitCCALblock::hitinfo_t()); |
218 | hiter->E_GeV = dEcorr/GeV; |
219 | hiter->t_ns = tcorr/ns; |
220 | } |
221 | else { |
222 | G4cerr(*G4cerr_p) << "GlueXSensitiveDetectorCCAL::ProcessHits error: " |
223 | << "max hit count " << MAX_HITS << " exceeded, truncating!" |
224 | << G4endlstd::endl; |
225 | } |
226 | } |
227 | return true; |
228 | } |
229 | |
230 | void GlueXSensitiveDetectorCCAL::EndOfEvent(G4HCofThisEvent*) |
231 | { |
232 | std::map<int,GlueXHitCCALblock*> *blocks = fBlocksMap->GetMap(); |
233 | std::map<int,GlueXHitCCALpoint*> *points = fPointsMap->GetMap(); |
234 | if (blocks->size() == 0 && points->size() == 0) |
235 | return; |
236 | std::map<int,GlueXHitCCALblock*>::iterator biter; |
237 | std::map<int,GlueXHitCCALpoint*>::iterator piter; |
238 | |
239 | if (verboseLevel > 1) { |
240 | G4cout(*G4cout_p) << G4endlstd::endl |
241 | << "--------> Hits Collection: in this event there are " |
242 | << blocks->size() << " blocks with hits in the CCAL: " |
243 | << G4endlstd::endl; |
244 | for (biter = blocks->begin(); biter != blocks->end(); ++biter) |
245 | biter->second->Print(); |
246 | |
247 | G4cout(*G4cout_p) << G4endlstd::endl |
248 | << "--------> Hits Collection: in this event there are " |
249 | << points->size() << " truth showers in the CCAL: " |
250 | << G4endlstd::endl; |
251 | for (piter = points->begin(); piter != points->end(); ++piter) |
252 | piter->second->Print(); |
253 | } |
254 | |
255 | |
256 | |
257 | G4EventManager* mgr = G4EventManager::GetEventManager(); |
258 | G4VUserEventInformation* info = mgr->GetUserInformation(); |
259 | hddm_s::HDDM *record = ((GlueXUserEventInformation*)info)->getOutputRecord(); |
260 | if (record == 0) { |
261 | G4cerr(*G4cerr_p) << "GlueXSensitiveDetectorCCAL::EndOfEvent error - " |
262 | << "hits seen but no output hddm record to save them into, " |
263 | << "cannot continue!" << G4endlstd::endl; |
264 | exit(1); |
265 | } |
266 | |
267 | if (record->getPhysicsEvents().size() == 0) |
268 | record->addPhysicsEvents(); |
269 | if (record->getHitViews().size() == 0) |
270 | record->getPhysicsEvent().addHitViews(); |
271 | hddm_s::HitView &hitview = record->getPhysicsEvent().getHitView(); |
272 | if (hitview.getComptonEMcals().size() == 0) |
273 | hitview.addComptonEMcals(); |
274 | hddm_s::ComptonEMcal &comptonEMcal = hitview.getComptonEMcal(); |
275 | |
276 | |
277 | for (biter = blocks->begin(); biter != blocks->end(); ++biter) { |
278 | std::vector<GlueXHitCCALblock::hitinfo_t> &hits = biter->second->hits; |
279 | |
280 | for (unsigned int ih=0; ih < hits.size(); ++ih) { |
281 | if (hits[ih].E_GeV <= THRESH_MEV/1e3) { |
282 | hits.erase(hits.begin() + ih); |
283 | --ih; |
284 | } |
285 | } |
286 | |
287 | if (hits.size() > 0) { |
288 | hddm_s::CcalBlockList block = comptonEMcal.addCcalBlocks(1); |
289 | block(0).setColumn(biter->second->column_); |
290 | block(0).setRow(biter->second->row_); |
291 | for (int ih=0; ih < (int)hits.size(); ++ih) { |
292 | hddm_s::CcalTruthHitList thit = block(0).addCcalTruthHits(1); |
293 | thit(0).setE(hits[ih].E_GeV); |
294 | thit(0).setT(hits[ih].t_ns); |
295 | } |
296 | } |
297 | } |
298 | |
299 | |
300 | for (piter = points->begin(); piter != points->end(); ++piter) { |
301 | hddm_s::CcalTruthShowerList point = comptonEMcal.addCcalTruthShowers(1); |
302 | point(0).setE(piter->second->E_GeV); |
303 | point(0).setPrimary(piter->second->primary_); |
304 | point(0).setPtype(piter->second->ptype_G3); |
305 | point(0).setPx(piter->second->px_GeV); |
306 | point(0).setPy(piter->second->py_GeV); |
307 | point(0).setPz(piter->second->pz_GeV); |
308 | point(0).setX(piter->second->x_cm); |
309 | point(0).setY(piter->second->y_cm); |
310 | point(0).setZ(piter->second->z_cm); |
311 | point(0).setT(piter->second->t_ns); |
312 | point(0).setTrack(piter->second->track_); |
313 | hddm_s::TrackIDList tid = point(0).addTrackIDs(); |
314 | tid(0).setItrack(piter->second->trackID_); |
315 | } |
316 | } |
317 | |
318 | int GlueXSensitiveDetectorCCAL::GetIdent(std::string div, |
319 | const G4VTouchable *touch) |
320 | { |
321 | const HddsG4Builder* bldr = GlueXDetectorConstruction::GetBuilder(); |
322 | std::map<std::string, std::vector<int> >::const_iterator iter; |
323 | std::map<std::string, std::vector<int> > *identifiers; |
324 | int max_depth = touch->GetHistoryDepth(); |
325 | for (int depth = 0; depth < max_depth; ++depth) { |
326 | G4VPhysicalVolume *pvol = touch->GetVolume(depth); |
327 | G4LogicalVolume *lvol = pvol->GetLogicalVolume(); |
328 | int volId = fVolumeTable[lvol]; |
329 | if (volId == 0) { |
330 | volId = bldr->getVolumeId(lvol); |
331 | fVolumeTable[lvol] = volId; |
332 | } |
333 | identifiers = &Refsys::fIdentifierTable[volId]; |
334 | if ((iter = identifiers->find(div)) != identifiers->end()) { |
335 | int copyNum = touch->GetCopyNumber(depth); |
336 | copyNum += (dynamic_cast<G4PVPlacement*>(pvol))? -1 : 0; |
337 | return iter->second[copyNum]; |
338 | } |
339 | } |
340 | return -1; |
341 | } |