1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | #include "GlueXSensitiveDetectorUPV.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 GlueXSensitiveDetectorUPV::MAX_HITS = 100; |
26 | |
27 | |
28 | double GlueXSensitiveDetectorUPV::ATTENUATION_LENGTH = 150.*cm; |
29 | double GlueXSensitiveDetectorUPV::C_EFFECTIVE = 19*cm/ns; |
30 | |
31 | |
32 | double GlueXSensitiveDetectorUPV::TWO_HIT_TIME_RESOL = 50*ns; |
33 | |
34 | |
35 | double GlueXSensitiveDetectorUPV::THRESH_MEV = 5.; |
36 | |
37 | int GlueXSensitiveDetectorUPV::instanceCount = 0; |
38 | G4Mutex GlueXSensitiveDetectorUPV::fMutex = G4MUTEX_INITIALIZER{ { 0, 0, 0, 0, 0, 0, 0, { 0, 0 } } }; |
39 | |
40 | GlueXSensitiveDetectorUPV::GlueXSensitiveDetectorUPV(const G4String& name) |
41 | : G4VSensitiveDetector(name), |
42 | fBarHitsMap(0), fPointsMap(0) |
43 | { |
44 | collectionName.insert("UPVBarHitsCollection"); |
45 | collectionName.insert("UPVPointsCollection"); |
46 | |
47 | |
48 | |
49 | |
50 | |
51 | |
52 | |
53 | G4AutoLock barrier(&fMutex); |
54 | if (instanceCount++ == 0) { |
55 | int runno = HddmOutput::getRunNo(); |
56 | extern jana::JApplication *japp; |
57 | if (japp == 0) { |
58 | G4cerr(*G4cerr_p) << "Error in GlueXSensitiveDetector constructor - " |
59 | << "jana global DApplication object not set, " |
60 | << "cannot continue." << G4endlstd::endl; |
61 | exit(-1); |
62 | } |
63 | jana::JCalibration *jcalib = japp->GetJCalibration(runno); |
| Value stored to 'jcalib' during its initialization is never read |
64 | if (japp == 0) { |
65 | jcalib = 0; |
66 | G4cout(*G4cout_p) << "UPV: ALL parameters loaded from ccdb" << G4endlstd::endl; |
67 | } |
68 | } |
69 | } |
70 | |
71 | GlueXSensitiveDetectorUPV::GlueXSensitiveDetectorUPV( |
72 | const GlueXSensitiveDetectorUPV &src) |
73 | : G4VSensitiveDetector(src), |
74 | fBarHitsMap(src.fBarHitsMap), fPointsMap(src.fPointsMap) |
75 | { |
76 | G4AutoLock barrier(&fMutex); |
77 | ++instanceCount; |
78 | } |
79 | |
80 | GlueXSensitiveDetectorUPV &GlueXSensitiveDetectorUPV::operator=(const |
81 | GlueXSensitiveDetectorUPV &src) |
82 | { |
83 | G4AutoLock barrier(&fMutex); |
84 | *(G4VSensitiveDetector*)this = src; |
85 | fBarHitsMap = src.fBarHitsMap; |
86 | fPointsMap = src.fPointsMap; |
87 | return *this; |
88 | } |
89 | |
90 | GlueXSensitiveDetectorUPV::~GlueXSensitiveDetectorUPV() |
91 | { |
92 | G4AutoLock barrier(&fMutex); |
93 | --instanceCount; |
94 | } |
95 | |
96 | void GlueXSensitiveDetectorUPV::Initialize(G4HCofThisEvent* hce) |
97 | { |
98 | fBarHitsMap = new |
99 | GlueXHitsMapUPVbar(SensitiveDetectorName, collectionName[0]); |
100 | fPointsMap = new |
101 | GlueXHitsMapUPVpoint(SensitiveDetectorName, collectionName[1]); |
102 | G4SDManager *sdm = G4SDManager::GetSDMpointer(); |
103 | hce->AddHitsCollection(sdm->GetCollectionID(collectionName[0]), fBarHitsMap); |
104 | hce->AddHitsCollection(sdm->GetCollectionID(collectionName[1]), fPointsMap); |
105 | } |
106 | |
107 | G4bool GlueXSensitiveDetectorUPV::ProcessHits(G4Step* step, |
108 | G4TouchableHistory* ROhist) |
109 | { |
110 | double dEsum = step->GetTotalEnergyDeposit(); |
111 | if (dEsum == 0) |
112 | return false; |
113 | |
114 | const G4ThreeVector &pin = step->GetPreStepPoint()->GetMomentum(); |
115 | const G4ThreeVector &xin = step->GetPreStepPoint()->GetPosition(); |
116 | const G4ThreeVector &xout = step->GetPostStepPoint()->GetPosition(); |
117 | double Ein = step->GetPreStepPoint()->GetTotalEnergy(); |
118 | double tin = step->GetPreStepPoint()->GetGlobalTime(); |
119 | double tout = step->GetPostStepPoint()->GetGlobalTime(); |
120 | G4ThreeVector x = (xin + xout) / 2; |
121 | G4ThreeVector dx = xout - xin; |
122 | double t = (tin + tout) / 2; |
123 | |
124 | const G4VTouchable* touch = step->GetPreStepPoint()->GetTouchable(); |
125 | const G4AffineTransform &local_from_global = touch->GetHistory() |
126 | ->GetTopTransform(); |
127 | G4ThreeVector xlocal = local_from_global.TransformPoint(x); |
128 | |
129 | |
130 | |
131 | |
132 | |
133 | |
134 | |
135 | if (tout > 1.0*s) |
136 | t = tin; |
137 | |
138 | |
139 | |
140 | |
141 | int layer = GetIdent("layer", touch); |
142 | int row = GetIdent("row", touch); |
143 | G4Track *track = step->GetTrack(); |
144 | G4int trackID = track->GetTrackID(); |
145 | GlueXUserTrackInformation *trackinfo = (GlueXUserTrackInformation*) |
146 | track->GetUserInformation(); |
147 | if (trackinfo->GetGlueXHistory() == 0 && Ein/MeV > THRESH_MEV) { |
148 | G4int key = fPointsMap->entries(); |
149 | GlueXHitUPVpoint* lastPoint = (*fPointsMap)[key - 1]; |
150 | if (lastPoint == 0 || lastPoint->track_ != trackID || |
151 | fabs(lastPoint->t_ns - t/ns) > 0.1 || |
152 | fabs(lastPoint->x_cm - x[0]/cm) > 2. || |
153 | fabs(lastPoint->y_cm - x[1]/cm) > 2. || |
154 | fabs(lastPoint->z_cm - x[2]/cm) > 2.) |
155 | { |
156 | int pdgtype = track->GetDynamicParticle()->GetPDGcode(); |
157 | int g3type = GlueXPrimaryGeneratorAction::ConvertPdgToGeant3(pdgtype); |
158 | GlueXHitUPVpoint newPoint; |
159 | newPoint.ptype_G3 = g3type; |
160 | newPoint.track_ = trackID; |
161 | newPoint.trackID_ = trackinfo->GetGlueXTrackID(); |
162 | newPoint.primary_ = (track->GetParentID() == 0); |
163 | newPoint.t_ns = t/ns; |
164 | newPoint.x_cm = x[0]/cm; |
165 | newPoint.y_cm = x[1]/cm; |
166 | newPoint.z_cm = x[2]/cm; |
167 | newPoint.px_GeV = pin[0]/GeV; |
168 | newPoint.py_GeV = pin[1]/GeV; |
169 | newPoint.pz_GeV = pin[2]/GeV; |
170 | newPoint.E_GeV = Ein/GeV; |
171 | fPointsMap->add(key, newPoint); |
172 | } |
173 | } |
174 | |
175 | |
176 | |
177 | if (dEsum > 0) { |
178 | int key = GlueXHitUPVbar::GetKey(layer, row); |
179 | GlueXHitUPVbar *counter = (*fBarHitsMap)[key]; |
180 | if (counter == 0) { |
181 | GlueXHitUPVbar newcounter(layer, row); |
182 | fBarHitsMap->add(key, newcounter); |
183 | counter = (*fBarHitsMap)[key]; |
184 | } |
185 | |
186 | double dxleft = xlocal[0]; |
187 | double dxright = -xlocal[0]; |
188 | |
189 | |
190 | |
191 | double tleft = t + dxleft / C_EFFECTIVE; |
192 | double tright = t + dxright / C_EFFECTIVE; |
193 | |
194 | |
195 | double dEleft = dEsum * exp(-dxleft / ATTENUATION_LENGTH); |
196 | double dEright = dEsum * exp(-dxright / ATTENUATION_LENGTH); |
197 | |
198 | |
199 | |
200 | if (dEleft > 0) { |
201 | |
202 | int merge_hit = 0; |
203 | std::vector<GlueXHitUPVbar::hitinfo_t>::iterator hiter; |
204 | for (hiter = counter->hits.begin(); hiter != counter->hits.end(); ++hiter) { |
205 | if (hiter->end_ == 0) { |
206 | if (fabs(hiter->t_ns*ns - tleft) < TWO_HIT_TIME_RESOL) { |
207 | merge_hit = 1; |
208 | break; |
209 | } |
210 | else if (hiter->t_ns*ns > tleft) { |
211 | break; |
212 | } |
213 | } |
214 | } |
215 | if (merge_hit) { |
216 | |
217 | hiter->E_GeV += dEleft/GeV; |
218 | if (hiter->t_ns*ns > tleft) { |
219 | hiter->t_ns = tleft/ns; |
220 | } |
221 | } |
222 | else if ((int)counter->hits.size() < MAX_HITS) { |
223 | |
224 | hiter = counter->hits.insert(hiter, GlueXHitUPVbar::hitinfo_t()); |
225 | hiter->end_ = 0; |
226 | hiter->E_GeV = dEleft/GeV; |
227 | hiter->t_ns = tleft/ns; |
228 | } |
229 | else { |
230 | G4cerr(*G4cerr_p) << "GlueXSensitiveDetectorUPV::ProcessHits error: " |
231 | << "max hit count " << MAX_HITS << " exceeded, truncating!" |
232 | << G4endlstd::endl; |
233 | } |
234 | } |
235 | |
236 | if (dEright/MeV > 0) { |
237 | |
238 | int merge_hit = 0; |
239 | std::vector<GlueXHitUPVbar::hitinfo_t>::iterator hiter; |
240 | for (hiter = counter->hits.begin(); hiter != counter->hits.end(); ++hiter) { |
241 | if (hiter->end_ == 1) { |
242 | if (fabs(hiter->t_ns*ns - tright) < TWO_HIT_TIME_RESOL) { |
243 | merge_hit = 1; |
244 | break; |
245 | } |
246 | else if (hiter->t_ns*ns > tright) { |
247 | break; |
248 | } |
249 | } |
250 | } |
251 | if (merge_hit) { |
252 | |
253 | hiter->E_GeV += dEright/GeV; |
254 | if (hiter->t_ns*ns > tright) { |
255 | hiter->t_ns = tright/ns; |
256 | } |
257 | } |
258 | else if ((int)counter->hits.size() < MAX_HITS) { |
259 | |
260 | hiter = counter->hits.insert(hiter, GlueXHitUPVbar::hitinfo_t()); |
261 | hiter->end_ = 1; |
262 | hiter->E_GeV = dEright/GeV; |
263 | hiter->t_ns = tright/ns; |
264 | } |
265 | else { |
266 | G4cerr(*G4cerr_p) << "GlueXSensitiveDetectorUPV::ProcessHits error: " |
267 | << "max hit count " << MAX_HITS << " exceeded, truncating!" |
268 | << G4endlstd::endl; |
269 | } |
270 | } |
271 | } |
272 | return true; |
273 | } |
274 | |
275 | void GlueXSensitiveDetectorUPV::EndOfEvent(G4HCofThisEvent*) |
276 | { |
277 | std::map<int,GlueXHitUPVbar*> *bars = fBarHitsMap->GetMap(); |
278 | std::map<int,GlueXHitUPVpoint*> *points = fPointsMap->GetMap(); |
279 | if (bars->size() == 0 && points->size() == 0) |
280 | return; |
281 | std::map<int,GlueXHitUPVbar*>::iterator siter; |
282 | std::map<int,GlueXHitUPVpoint*>::iterator piter; |
283 | |
284 | if (verboseLevel > 1) { |
285 | G4cout(*G4cout_p) << G4endlstd::endl |
286 | << "--------> Hits Collection: in this event there are " |
287 | << bars->size() << " bars with hits in the UPV: " |
288 | << G4endlstd::endl; |
289 | for (siter = bars->begin(); siter != bars->end(); ++siter) |
290 | siter->second->Print(); |
291 | |
292 | G4cout(*G4cout_p) << G4endlstd::endl |
293 | << "--------> Hits Collection: in this event there are " |
294 | << points->size() << " truth points in the UPV: " |
295 | << G4endlstd::endl; |
296 | for (piter = points->begin(); piter != points->end(); ++piter) |
297 | piter->second->Print(); |
298 | } |
299 | |
300 | |
301 | |
302 | G4EventManager* mgr = G4EventManager::GetEventManager(); |
303 | G4VUserEventInformation* info = mgr->GetUserInformation(); |
304 | hddm_s::HDDM *record = ((GlueXUserEventInformation*)info)->getOutputRecord(); |
305 | if (record == 0) { |
306 | G4cerr(*G4cerr_p) << "GlueXSensitiveDetectorUPV::EndOfEvent error - " |
307 | << "hits seen but no output hddm record to save them into, " |
308 | << "cannot continue!" << G4endlstd::endl; |
309 | exit(1); |
310 | } |
311 | |
312 | if (record->getPhysicsEvents().size() == 0) |
313 | record->addPhysicsEvents(); |
314 | if (record->getHitViews().size() == 0) |
315 | record->getPhysicsEvent().addHitViews(); |
316 | hddm_s::HitView &hitview = record->getPhysicsEvent().getHitView(); |
317 | if (hitview.getUpstreamEMvetos().size() == 0) |
318 | hitview.addUpstreamEMvetos(); |
319 | hddm_s::UpstreamEMveto &upv = hitview.getUpstreamEMveto(); |
320 | |
321 | |
322 | for (siter = bars->begin(); siter != bars->end(); ++siter) { |
323 | std::vector<GlueXHitUPVbar::hitinfo_t> &hits = siter->second->hits; |
324 | |
325 | for (unsigned int ih=0; ih < hits.size(); ++ih) { |
326 | if (hits[ih].E_GeV*1e3 <= THRESH_MEV) { |
327 | hits.erase(hits.begin() + ih); |
328 | --ih; |
329 | } |
330 | } |
331 | if (hits.size() > 0) { |
332 | hddm_s::UpvPaddleList counter = upv.addUpvPaddles(1); |
333 | counter(0).setLayer(siter->second->layer_); |
334 | counter(0).setRow(siter->second->row_); |
335 | |
336 | for (int ih=0; ih < (int)hits.size(); ++ih) { |
337 | if (hits[ih].end_ == 0) { |
338 | hddm_s::UpvTruthHitList thit = counter(0).addUpvTruthHits(1); |
339 | thit(0).setEnd(hits[ih].end_); |
340 | thit(0).setE(hits[ih].E_GeV); |
341 | thit(0).setT(hits[ih].t_ns); |
342 | } |
343 | } |
344 | |
345 | for (int ih=0; ih < (int)hits.size(); ++ih) { |
346 | if (hits[ih].end_ == 1) { |
347 | hddm_s::UpvTruthHitList thit = counter(0).addUpvTruthHits(1); |
348 | thit(0).setEnd(hits[ih].end_); |
349 | thit(0).setE(hits[ih].E_GeV); |
350 | thit(0).setT(hits[ih].t_ns); |
351 | } |
352 | } |
353 | } |
354 | } |
355 | |
356 | |
357 | for (piter = points->begin(); piter != points->end(); ++piter) { |
358 | hddm_s::UpvTruthShowerList point = upv.addUpvTruthShowers(1); |
359 | point(0).setPrimary(piter->second->primary_); |
360 | point(0).setPtype(piter->second->ptype_G3); |
361 | point(0).setPx(piter->second->px_GeV); |
362 | point(0).setPy(piter->second->py_GeV); |
363 | point(0).setPz(piter->second->pz_GeV); |
364 | point(0).setE(piter->second->E_GeV); |
365 | point(0).setX(piter->second->x_cm); |
366 | point(0).setY(piter->second->y_cm); |
367 | point(0).setZ(piter->second->z_cm); |
368 | point(0).setT(piter->second->t_ns); |
369 | point(0).setTrack(piter->second->track_); |
370 | hddm_s::TrackIDList tid = point(0).addTrackIDs(); |
371 | tid(0).setItrack(piter->second->trackID_); |
372 | } |
373 | } |
374 | |
375 | int GlueXSensitiveDetectorUPV::GetIdent(std::string div, |
376 | const G4VTouchable *touch) |
377 | { |
378 | const HddsG4Builder* bldr = GlueXDetectorConstruction::GetBuilder(); |
379 | std::map<std::string, std::vector<int> >::const_iterator iter; |
380 | std::map<std::string, std::vector<int> > *identifiers; |
381 | int max_depth = touch->GetHistoryDepth(); |
382 | for (int depth = 0; depth < max_depth; ++depth) { |
383 | G4VPhysicalVolume *pvol = touch->GetVolume(depth); |
384 | G4LogicalVolume *lvol = pvol->GetLogicalVolume(); |
385 | int volId = fVolumeTable[lvol]; |
386 | if (volId == 0) { |
387 | volId = bldr->getVolumeId(lvol); |
388 | fVolumeTable[lvol] = volId; |
389 | } |
390 | identifiers = &Refsys::fIdentifierTable[volId]; |
391 | if ((iter = identifiers->find(div)) != identifiers->end()) { |
392 | int copyNum = touch->GetCopyNumber(depth); |
393 | copyNum += (dynamic_cast<G4PVPlacement*>(pvol))? -1 : 0; |
394 | return iter->second[copyNum]; |
395 | } |
396 | } |
397 | return -1; |
398 | } |