Bug Summary

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

Annotated Source Code

1//
2// GlueXSensitiveDetectorPSC - class implementation
3//
4// author: richard.t.jones at uconn.edu
5// version: october 26, 2016
6
7#include "GlueXSensitiveDetectorPSC.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 GlueXSensitiveDetectorPSC::MAX_HITS = 100;
26
27// Minimum hit time difference for two hits on the same paddle
28double GlueXSensitiveDetectorPSC::TWO_HIT_TIME_RESOL = 25*ns;
29
30// Minimum energy deposition for a hit
31double GlueXSensitiveDetectorPSC::THRESH_MEV = 0.010;
32
33// Geometric parameters
34int GlueXSensitiveDetectorPSC::NUM_MODULES_PER_ARM = 8;
35
36int GlueXSensitiveDetectorPSC::instanceCount = 0;
37G4Mutex GlueXSensitiveDetectorPSC::fMutex = G4MUTEX_INITIALIZER{ { 0, 0, 0, 0, 0, 0, 0, { 0, 0 } } };
38
39GlueXSensitiveDetectorPSC::GlueXSensitiveDetectorPSC(const G4String& name)
40 : G4VSensitiveDetector(name),
41 fCounterHitsMap(0), fPointsMap(0)
42{
43 collectionName.insert("PSCPaddleHitsCollection");
44 collectionName.insert("PSCPointsCollection");
45
46 // The rest of this only needs to happen once, the first time an object
47 // of this type is instantiated for this configuration of geometry and
48 // fields. If the geometry or fields change in such a way as to modify
49 // the drift-time properties of hits in the PSC, you must delete all old
50 // objects of this class and create new ones.
51
52 G4AutoLock barrier(&fMutex);
53 if (instanceCount++ == 0) {
54 int runno = HddmOutput::getRunNo();
55 extern jana::JApplication *japp;
56 if (japp == 0) {
57 G4cerr(*G4cerr_p) << "Error in GlueXSensitiveDetector constructor - "
58 << "jana global DApplication object not set, "
59 << "cannot continue." << G4endlstd::endl;
60 exit(-1);
61 }
62 jana::JCalibration *jcalib = japp->GetJCalibration(runno);
Value stored to 'jcalib' during its initialization is never read
63 if (japp == 0) { // dummy
64 jcalib = 0;
65 G4cout(*G4cout_p) << "PSC: ALL parameters loaded from ccdb" << G4endlstd::endl;
66 }
67 }
68}
69
70GlueXSensitiveDetectorPSC::GlueXSensitiveDetectorPSC(
71 const GlueXSensitiveDetectorPSC &src)
72 : G4VSensitiveDetector(src),
73 fCounterHitsMap(src.fCounterHitsMap), fPointsMap(src.fPointsMap)
74{
75 G4AutoLock barrier(&fMutex);
76 ++instanceCount;
77}
78
79GlueXSensitiveDetectorPSC &GlueXSensitiveDetectorPSC::operator=(const
80 GlueXSensitiveDetectorPSC &src)
81{
82 G4AutoLock barrier(&fMutex);
83 *(G4VSensitiveDetector*)this = src;
84 fCounterHitsMap = src.fCounterHitsMap;
85 fPointsMap = src.fPointsMap;
86 return *this;
87}
88
89GlueXSensitiveDetectorPSC::~GlueXSensitiveDetectorPSC()
90{
91 G4AutoLock barrier(&fMutex);
92 --instanceCount;
93}
94
95void GlueXSensitiveDetectorPSC::Initialize(G4HCofThisEvent* hce)
96{
97 fCounterHitsMap = new
98 GlueXHitsMapPSCpaddle(SensitiveDetectorName, collectionName[0]);
99 fPointsMap = new
100 GlueXHitsMapPSCpoint(SensitiveDetectorName, collectionName[1]);
101 G4SDManager *sdm = G4SDManager::GetSDMpointer();
102 hce->AddHitsCollection(sdm->GetCollectionID(collectionName[0]), fCounterHitsMap);
103 hce->AddHitsCollection(sdm->GetCollectionID(collectionName[1]), fPointsMap);
104}
105
106G4bool GlueXSensitiveDetectorPSC::ProcessHits(G4Step* step,
107 G4TouchableHistory* ROhist)
108{
109 double dEsum = step->GetTotalEnergyDeposit();
110 if (dEsum == 0)
111 return false;
112
113 const G4ThreeVector &pin = step->GetPreStepPoint()->GetMomentum();
114 const G4ThreeVector &xin = step->GetPreStepPoint()->GetPosition();
115 const G4ThreeVector &xout = step->GetPostStepPoint()->GetPosition();
116 double Ein = step->GetPreStepPoint()->GetTotalEnergy();
117 double tin = step->GetPreStepPoint()->GetGlobalTime();
118 double tout = step->GetPostStepPoint()->GetGlobalTime();
119 G4ThreeVector x = (xin + xout) / 2;
120 G4ThreeVector dx = xout - xin;
121 double t = (tin + tout) / 2;
122
123 const G4VTouchable* touch = step->GetPreStepPoint()->GetTouchable();
124 const G4AffineTransform &local_from_global = touch->GetHistory()
125 ->GetTopTransform();
126 G4ThreeVector xlocal = local_from_global.TransformPoint(x);
127
128 // For particles that range out inside the active volume, the
129 // "out" time may sometimes be set to something enormously high.
130 // This screws up the hit. Check for this case here by looking
131 // at tout and making sure it is less than 1 second. If it's
132 // not, then just use tin for "t".
133
134 if (tout > 1.0*s)
135 t = tin;
136
137 double dr = dx.mag();
138 double dEdx = (dr > 1e-3*cm)? dEsum/dr : 0;
139
140 // Post the hit to the points list in the
141 // order of appearance in the event simulation.
142
143 int module = GetIdent("module", touch);
144 int arm = (module - 1) / NUM_MODULES_PER_ARM;
145 module = (module - 1) % NUM_MODULES_PER_ARM + 1;
146 G4Track *track = step->GetTrack();
147 G4int trackID = track->GetTrackID();
148 int pdgtype = track->GetDynamicParticle()->GetPDGcode();
149 int g3type = GlueXPrimaryGeneratorAction::ConvertPdgToGeant3(pdgtype);
150 GlueXUserTrackInformation *trackinfo = (GlueXUserTrackInformation*)
151 track->GetUserInformation();
152 int itrack = trackinfo->GetGlueXTrackID();
153 if (trackinfo->GetGlueXHistory() == 0 && itrack > 0) {
154 G4int key = fPointsMap->entries();
155 GlueXHitPSCpoint* lastPoint = (*fPointsMap)[key - 1];
156 if (lastPoint == 0 || lastPoint->track_ != trackID ||
157 fabs(lastPoint->t_ns - t/ns) > 0.2 ||
158 fabs(lastPoint->x_cm - x[0]/cm) > 5.0 ||
159 fabs(lastPoint->y_cm - x[1]/cm) > 5.0 ||
160 fabs(lastPoint->z_cm - x[2]/cm) > 5.0)
161 {
162 GlueXHitPSCpoint newPoint;
163 newPoint.arm_ = arm;
164 newPoint.module_ = module;
165 newPoint.ptype_G3 = g3type;
166 newPoint.track_ = trackID;
167 newPoint.trackID_ = itrack;
168 newPoint.primary_ = (track->GetParentID() == 0);
169 newPoint.t_ns = t/ns;
170 newPoint.x_cm = x[0]/cm;
171 newPoint.y_cm = x[1]/cm;
172 newPoint.z_cm = x[2]/cm;
173 newPoint.px_GeV = pin[0]/GeV;
174 newPoint.py_GeV = pin[1]/GeV;
175 newPoint.pz_GeV = pin[2]/GeV;
176 newPoint.E_GeV = Ein/GeV;
177 newPoint.dEdx_GeV_cm = dEdx/(GeV/cm);
178 fPointsMap->add(key, newPoint);
179 }
180 }
181
182 // Post the hit to the hits map, ordered by arm,module index
183
184 if (dEsum > 0) {
185 int key = GlueXHitPSCpaddle::GetKey(arm, module);
186 GlueXHitPSCpaddle *paddle = (*fCounterHitsMap)[key];
187 if (paddle == 0) {
188 GlueXHitPSCpaddle newpaddle(arm, module);
189 fCounterHitsMap->add(key, newpaddle);
190 paddle = (*fCounterHitsMap)[key];
191 }
192
193 // Add the hit to the hits vector, maintaining strict time ordering
194
195 int merge_hit = 0;
196 std::vector<GlueXHitPSCpaddle::hitinfo_t>::iterator hiter;
197 for (hiter = paddle->hits.begin(); hiter != paddle->hits.end(); ++hiter) {
198 if (fabs(hiter->t_ns*ns - t) < TWO_HIT_TIME_RESOL) {
199 merge_hit = 1;
200 break;
201 }
202 else if (hiter->t_ns*ns > t) {
203 break;
204 }
205 }
206 if (merge_hit) {
207 // Add the charge, do energy-weighted time averaging
208 hiter->t_ns = (hiter->t_ns * hiter->dE_GeV + t/ns * dEsum/GeV) /
209 (hiter->dE_GeV + dEsum/GeV);
210 hiter->dE_GeV += dEsum/GeV;
211 }
212 else if ((int)paddle->hits.size() < MAX_HITS) {
213 // create new hit
214 hiter = paddle->hits.insert(hiter, GlueXHitPSCpaddle::hitinfo_t());
215 hiter->dE_GeV = dEsum/GeV;
216 hiter->t_ns = t/ns;
217 hiter->itrack_ = itrack;
218 hiter->ptype_G3 = g3type;
219 }
220 else {
221 G4cerr(*G4cerr_p) << "GlueXSensitiveDetectorPSC::ProcessHits error: "
222 << "max hit count " << MAX_HITS << " exceeded, truncating!"
223 << G4endlstd::endl;
224 }
225 }
226 return true;
227}
228
229void GlueXSensitiveDetectorPSC::EndOfEvent(G4HCofThisEvent*)
230{
231 std::map<int,GlueXHitPSCpaddle*> *paddles = fCounterHitsMap->GetMap();
232 std::map<int,GlueXHitPSCpoint*> *points = fPointsMap->GetMap();
233 if (paddles->size() == 0 && points->size() == 0)
234 return;
235 std::map<int,GlueXHitPSCpaddle*>::iterator siter;
236 std::map<int,GlueXHitPSCpoint*>::iterator piter;
237
238 if (verboseLevel > 1) {
239 G4cout(*G4cout_p) << G4endlstd::endl
240 << "--------> Hits Collection: in this event there are "
241 << paddles->size() << " paddles with hits in the PSC: "
242 << G4endlstd::endl;
243 for (siter = paddles->begin(); siter != paddles->end(); ++siter)
244 siter->second->Print();
245
246 G4cout(*G4cout_p) << G4endlstd::endl
247 << "--------> Hits Collection: in this event there are "
248 << points->size() << " truth points in the PSC: "
249 << G4endlstd::endl;
250 for (piter = points->begin(); piter != points->end(); ++piter)
251 piter->second->Print();
252 }
253
254 // pack hits into ouptut hddm record
255
256 G4EventManager* mgr = G4EventManager::GetEventManager();
257 G4VUserEventInformation* info = mgr->GetUserInformation();
258 hddm_s::HDDM *record = ((GlueXUserEventInformation*)info)->getOutputRecord();
259 if (record == 0) {
260 G4cerr(*G4cerr_p) << "GlueXSensitiveDetectorPSC::EndOfEvent error - "
261 << "hits seen but no output hddm record to save them into, "
262 << "cannot continue!" << G4endlstd::endl;
263 exit(1);
264 }
265
266 if (record->getPhysicsEvents().size() == 0)
267 record->addPhysicsEvents();
268 if (record->getHitViews().size() == 0)
269 record->getPhysicsEvent().addHitViews();
270 hddm_s::HitView &hitview = record->getPhysicsEvent().getHitView();
271 if (hitview.getPairSpectrometerCoarses().size() == 0)
272 hitview.addPairSpectrometerCoarses();
273 hddm_s::PairSpectrometerCoarse &psc = hitview.getPairSpectrometerCoarse();
274
275 // Collect and output the PscTruthHits
276 for (siter = paddles->begin(); siter != paddles->end(); ++siter) {
277 std::vector<GlueXHitPSCpaddle::hitinfo_t> &hits = siter->second->hits;
278 // apply a pulse height threshold cut
279 for (unsigned int ih=0; ih < hits.size(); ++ih) {
280 if (hits[ih].dE_GeV*1e3 <= THRESH_MEV) {
281 hits.erase(hits.begin() + ih);
282 --ih;
283 }
284 }
285 if (hits.size() > 0) {
286 hddm_s::PscPaddleList paddle = psc.addPscPaddles(1);
287 paddle(0).setArm(siter->second->arm_);
288 paddle(0).setModule(siter->second->module_);
289 for (int ih=0; ih < (int)hits.size(); ++ih) {
290 hddm_s::PscTruthHitList thit = paddle(0).addPscTruthHits(1);
291 thit(0).setDE(hits[ih].dE_GeV);
292 thit(0).setT(hits[ih].t_ns);
293 thit(0).setItrack(hits[ih].itrack_);
294 thit(0).setPtype(hits[ih].ptype_G3);
295 }
296 }
297 }
298
299 // Collect and output the pscTruthPoints
300 for (piter = points->begin(); piter != points->end(); ++piter) {
301 hddm_s::PscTruthPointList point = psc.addPscTruthPoints(1);
302 point(0).setE(piter->second->E_GeV);
303 point(0).setDEdx(piter->second->dEdx_GeV_cm);
304 point(0).setPrimary(piter->second->primary_);
305 point(0).setPtype(piter->second->ptype_G3);
306 point(0).setPx(piter->second->px_GeV);
307 point(0).setPy(piter->second->py_GeV);
308 point(0).setPz(piter->second->pz_GeV);
309 point(0).setX(piter->second->x_cm);
310 point(0).setY(piter->second->y_cm);
311 point(0).setZ(piter->second->z_cm);
312 point(0).setT(piter->second->t_ns);
313 point(0).setArm(piter->second->arm_);
314 point(0).setModule(piter->second->module_);
315 point(0).setTrack(piter->second->track_);
316 hddm_s::TrackIDList tid = point(0).addTrackIDs();
317 tid(0).setItrack(piter->second->trackID_);
318 }
319}
320
321int GlueXSensitiveDetectorPSC::GetIdent(std::string div,
322 const G4VTouchable *touch)
323{
324 const HddsG4Builder* bldr = GlueXDetectorConstruction::GetBuilder();
325 std::map<std::string, std::vector<int> >::const_iterator iter;
326 std::map<std::string, std::vector<int> > *identifiers;
327 int max_depth = touch->GetHistoryDepth();
328 for (int depth = 0; depth < max_depth; ++depth) {
329 G4VPhysicalVolume *pvol = touch->GetVolume(depth);
330 G4LogicalVolume *lvol = pvol->GetLogicalVolume();
331 int volId = fVolumeTable[lvol];
332 if (volId == 0) {
333 volId = bldr->getVolumeId(lvol);
334 fVolumeTable[lvol] = volId;
335 }
336 identifiers = &Refsys::fIdentifierTable[volId];
337 if ((iter = identifiers->find(div)) != identifiers->end()) {
338 int copyNum = touch->GetCopyNumber(depth);
339 copyNum += (dynamic_cast<G4PVPlacement*>(pvol))? -1 : 0;
340 return iter->second[copyNum];
341 }
342 }
343 return -1;
344}