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 total number of allowed hits
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 if ((int)block->hits.size() < MAX_HITS) {
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 else {
217 G4cerr(*G4cerr_p) << "GlueXSensitiveDetectorGCAL::ProcessHits error: "
218 << "max hit count " << MAX_HITS << " exceeded, truncating!"
219 << G4endlstd::endl;
220 }
221 }
222 return true;
223}
224
225void GlueXSensitiveDetectorGCAL::EndOfEvent(G4HCofThisEvent*)
226{
227 std::map<int,GlueXHitGCALblock*> *blocks = fBlockHitsMap->GetMap();
228 std::map<int,GlueXHitGCALpoint*> *points = fPointsMap->GetMap();
229 if (blocks->size() == 0 && points->size() == 0)
230 return;
231 std::map<int,GlueXHitGCALblock*>::iterator siter;
232 std::map<int,GlueXHitGCALpoint*>::iterator piter;
233
234 if (verboseLevel > 1) {
235 G4cout(*G4cout_p) << G4endlstd::endl
236 << "--------> Hits Collection: in this event there are "
237 << blocks->size() << " blocks with hits in the GCAL: "
238 << G4endlstd::endl;
239 for (siter = blocks->begin(); siter != blocks->end(); ++siter)
240 siter->second->Print();
241
242 G4cout(*G4cout_p) << G4endlstd::endl
243 << "--------> Hits Collection: in this event there are "
244 << points->size() << " truth points in the GCAL: "
245 << G4endlstd::endl;
246 for (piter = points->begin(); piter != points->end(); ++piter)
247 piter->second->Print();
248 }
249
250 // pack hits into ouptut hddm record
251
252 G4EventManager* mgr = G4EventManager::GetEventManager();
253 G4VUserEventInformation* info = mgr->GetUserInformation();
254 hddm_s::HDDM *record = ((GlueXUserEventInformation*)info)->getOutputRecord();
255 if (record == 0) {
256 G4cerr(*G4cerr_p) << "GlueXSensitiveDetectorGCAL::EndOfEvent error - "
257 << "hits seen but no output hddm record to save them into, "
258 << "cannot continue!" << G4endlstd::endl;
259 exit(1);
260 }
261
262 if (record->getPhysicsEvents().size() == 0)
263 record->addPhysicsEvents();
264 if (record->getHitViews().size() == 0)
265 record->getPhysicsEvent().addHitViews();
266 hddm_s::HitView &hitview = record->getPhysicsEvent().getHitView();
267 if (hitview.getGapEMcals().size() == 0)
268 hitview.addGapEMcals();
269 hddm_s::GapEMcal &gcal = hitview.getGapEMcal();
270
271 // Collect and output the gcalTruthHits
272 for (siter = blocks->begin(); siter != blocks->end(); ++siter) {
273 std::vector<GlueXHitGCALblock::hitinfo_t> &hits = siter->second->hits;
274 // apply a pulse height threshold cut
275 for (unsigned int ih=0; ih < hits.size(); ++ih) {
276 if (hits[ih].E_GeV*1e3 <= THRESH_MEV) {
277 hits.erase(hits.begin() + ih);
278 --ih;
279 }
280 }
281 if (hits.size() > 0) {
282 hddm_s::GcalCellList block = gcal.addGcalCells(1);
283 block(0).setModule(siter->second->module_);
284 for (int ih=0; ih < (int)hits.size(); ++ih) {
285 hddm_s::GcalTruthHitList thit = block(0).addGcalTruthHits(1);
286 thit(0).setE(hits[ih].E_GeV);
287 thit(0).setT(hits[ih].t_ns);
288 }
289 }
290 }
291
292 // Collect and output the gcalTruthShowers
293 for (piter = points->begin(); piter != points->end(); ++piter) {
294 hddm_s::GcalTruthShowerList point = gcal.addGcalTruthShowers(1);
295 point(0).setE(piter->second->E_GeV);
296 point(0).setPrimary(piter->second->primary_);
297 point(0).setPtype(piter->second->ptype_G3);
298 point(0).setPx(piter->second->px_GeV);
299 point(0).setPy(piter->second->py_GeV);
300 point(0).setPz(piter->second->pz_GeV);
301 point(0).setR(piter->second->r_cm);
302 point(0).setPhi(piter->second->phi_rad);
303 point(0).setZ(piter->second->z_cm);
304 point(0).setT(piter->second->t_ns);
305 point(0).setTrack(piter->second->track_);
306 hddm_s::TrackIDList tid = point(0).addTrackIDs();
307 tid(0).setItrack(piter->second->trackID_);
308 }
309}
310
311int GlueXSensitiveDetectorGCAL::GetIdent(std::string div,
312 const G4VTouchable *touch)
313{
314 const HddsG4Builder* bldr = GlueXDetectorConstruction::GetBuilder();
315 std::map<std::string, std::vector<int> >::const_iterator iter;
316 std::map<std::string, std::vector<int> > *identifiers;
317 int max_depth = touch->GetHistoryDepth();
318 for (int depth = 0; depth < max_depth; ++depth) {
319 G4VPhysicalVolume *pvol = touch->GetVolume(depth);
320 G4LogicalVolume *lvol = pvol->GetLogicalVolume();
321 int volId = fVolumeTable[lvol];
322 if (volId == 0) {
323 volId = bldr->getVolumeId(lvol);
324 fVolumeTable[lvol] = volId;
325 }
326 identifiers = &Refsys::fIdentifierTable[volId];
327 if ((iter = identifiers->find(div)) != identifiers->end()) {
328 int copyNum = touch->GetCopyNumber(depth);
329 copyNum += (dynamic_cast<G4PVPlacement*>(pvol))? -1 : 0;
330 return iter->second[copyNum];
331 }
332 }
333 return -1;
334}