Bug Summary

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

Annotated Source Code

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