Bug Summary

File:scratch/gluex/scan-build-work/hdgeant4/src/GlueXSensitiveDetectorUPV.cc
Location:line 63, column 27
Description:Value stored to 'jcalib' during its initialization is never read

Annotated Source Code

1//
2// GlueXSensitiveDetectorUPV - class implementation
3//
4// author: richard.t.jones at uconn.edu
5// version: october 26, 2016
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// Cutoff on the total number of allowed hits
25int GlueXSensitiveDetectorUPV::MAX_HITS = 100;
26
27// Light propagation parameters in tof bars
28double GlueXSensitiveDetectorUPV::ATTENUATION_LENGTH = 150.*cm;
29double GlueXSensitiveDetectorUPV::C_EFFECTIVE = 19*cm/ns;
30
31// Minimum hit time difference for two hits on the same bar
32double GlueXSensitiveDetectorUPV::TWO_HIT_TIME_RESOL = 50*ns;
33
34// Minimum energy deposition for a hit
35double GlueXSensitiveDetectorUPV::THRESH_MEV = 5.;
36
37int GlueXSensitiveDetectorUPV::instanceCount = 0;
38G4Mutex GlueXSensitiveDetectorUPV::fMutex = G4MUTEX_INITIALIZER{ { 0, 0, 0, 0, 0, 0, 0, { 0, 0 } } };
39
40GlueXSensitiveDetectorUPV::GlueXSensitiveDetectorUPV(const G4String& name)
41 : G4VSensitiveDetector(name),
42 fBarHitsMap(0), fPointsMap(0)
43{
44 collectionName.insert("UPVBarHitsCollection");
45 collectionName.insert("UPVPointsCollection");
46
47 // The rest of this only needs to happen once, the first time an object
48 // of this type is instantiated for this configuration of geometry and
49 // fields. If the geometry or fields change in such a way as to modify
50 // the drift-time properties of hits in the UPV, you must delete all old
51 // objects of this class and create new ones.
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) { // dummy
65 jcalib = 0;
66 G4cout(*G4cout_p) << "UPV: ALL parameters loaded from ccdb" << G4endlstd::endl;
67 }
68 }
69}
70
71GlueXSensitiveDetectorUPV::GlueXSensitiveDetectorUPV(
72 const GlueXSensitiveDetectorUPV &src)
73 : G4VSensitiveDetector(src),
74 fBarHitsMap(src.fBarHitsMap), fPointsMap(src.fPointsMap)
75{
76 G4AutoLock barrier(&fMutex);
77 ++instanceCount;
78}
79
80GlueXSensitiveDetectorUPV &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
90GlueXSensitiveDetectorUPV::~GlueXSensitiveDetectorUPV()
91{
92 G4AutoLock barrier(&fMutex);
93 --instanceCount;
94}
95
96void 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
107G4bool 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 // For particles that range out inside the active volume, the
130 // "out" time may sometimes be set to something enormously high.
131 // This screws up the hit. Check for this case here by looking
132 // at tout and making sure it is less than 1 second. If it's
133 // not, then just use tin for "t".
134
135 if (tout > 1.0*s)
136 t = tin;
137
138 // Post the hit to the points list in the
139 // order of appearance in the event simulation.
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 // Post the hit to the hits map, ordered by layer,row,end index
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 // Calculate time at the PMT "normalized" to the center,
190 // so a hit in the center will have time "t" at both PMTs.
191 double tleft = t + dxleft / C_EFFECTIVE;
192 double tright = t + dxright / C_EFFECTIVE;
193
194 // calculate energy seen by PM for this track step using attenuation factor
195 double dEleft = dEsum * exp(-dxleft / ATTENUATION_LENGTH);
196 double dEright = dEsum * exp(-dxright / ATTENUATION_LENGTH);
197
198 // Add the hit to the hits vector, maintaining strict time ordering
199
200 if (dEleft > 0) {
201 // add the hit on end=0 (north/left end of the bar)
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 // Use the time from the earlier hit but add the charge
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 // create new hit
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 // add the hit on end=1 (south/bottom end of the bar)
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 // Use the time from the earlier hit but add the charge
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 // create new hit
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
275void 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 // pack hits into ouptut hddm record
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 // Collect and output the tofTruthHits
322 for (siter = bars->begin(); siter != bars->end(); ++siter) {
323 std::vector<GlueXHitUPVbar::hitinfo_t> &hits = siter->second->hits;
324 // apply a pulse height threshold cut
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 // first the end=0 hits
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 // followed by the end=1 hits
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 // Collect and output the barTruthShowers
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
375int 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}