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 (trackinfo->GetGlueXHistory() == 0 && itrack > 0 &&
149 xin.dot(pin) > 0 && Ein/MeV > THRESH_MEV)
150 {
151 int pdgtype = track->GetDynamicParticle()->GetPDGcode();
152 int g3type = GlueXPrimaryGeneratorAction::ConvertPdgToGeant3(pdgtype);
153 GlueXHitGCALpoint newPoint;
154 newPoint.ptype_G3 = g3type;
155 newPoint.track_ = trackID;
156 newPoint.trackID_ = itrack;
157 newPoint.primary_ = (track->GetParentID() == 0);
158 newPoint.t_ns = t/ns;
159 newPoint.z_cm = x[2]/cm;
160 newPoint.r_cm = x.perp()/cm;
161 newPoint.phi_rad = x.phi();
162 newPoint.px_GeV = pin[0]/GeV;
163 newPoint.py_GeV = pin[1]/GeV;
164 newPoint.pz_GeV = pin[2]/GeV;
165 newPoint.E_GeV = Ein/GeV;
166 G4int key = fPointsMap->entries();
167 fPointsMap->add(key, newPoint);
168 trackinfo->SetGlueXHistory(3);
169 }
170
171 // Post the hit to the hits map, ordered by module index
172
173 if (dEsum > 0) {
174 int key = GlueXHitGCALblock::GetKey(module);
175 GlueXHitGCALblock *block = (*fBlockHitsMap)[key];
176 if (block == 0) {
177 GlueXHitGCALblock newblock(module);
178 fBlockHitsMap->add(key, newblock);
179 block = (*fBlockHitsMap)[key];
180 }
181 double dist = 0.5 * LENGTH_OF_BLOCK - xlocal[2];
182 double dEcorr = dEsum * exp(-dist / ATTENUATION_LENGTH);
183 double tcorr = t + dist / C_EFFECTIVE;
184
185 // Add the hit to the hits vector, maintaining strict time ordering
186
187 int merge_hit = 0;
188 std::vector<GlueXHitGCALblock::hitinfo_t>::iterator hiter;
189 for (hiter = block->hits.begin(); hiter != block->hits.end(); ++hiter) {
190 if (fabs(hiter->t_ns*ns - tcorr) < TWO_HIT_TIME_RESOL) {
191 merge_hit = 1;
192 break;
193 }
194 else if (hiter->t_ns*ns > tcorr) {
195 break;
196 }
197 }
198 if (merge_hit) {
199 // Use the time from the earlier hit but add the charge
200 hiter->E_GeV += dEcorr/GeV;
201 if (hiter->t_ns*ns > tcorr) {
202 hiter->t_ns = tcorr/ns;
203 }
204 }
205 else if ((int)block->hits.size() < MAX_HITS) {
206 // create new hit
207 hiter = block->hits.insert(hiter, GlueXHitGCALblock::hitinfo_t());
208 hiter->E_GeV = dEcorr/GeV;
209 hiter->t_ns = tcorr/ns;
210 }
211 else {
212 G4cerr(*G4cerr_p) << "GlueXSensitiveDetectorGCAL::ProcessHits error: "
213 << "max hit count " << MAX_HITS << " exceeded, truncating!"
214 << G4endlstd::endl;
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 for (int ih=0; ih < (int)hits.size(); ++ih) {
280 hddm_s::GcalTruthHitList thit = block(0).addGcalTruthHits(1);
281 thit(0).setE(hits[ih].E_GeV);
282 thit(0).setT(hits[ih].t_ns);
283 }
284 }
285 }
286
287 // Collect and output the gcalTruthShowers
288 for (piter = points->begin(); piter != points->end(); ++piter) {
289 hddm_s::GcalTruthShowerList point = gcal.addGcalTruthShowers(1);
290 point(0).setE(piter->second->E_GeV);
291 point(0).setPrimary(piter->second->primary_);
292 point(0).setPtype(piter->second->ptype_G3);
293 point(0).setPx(piter->second->px_GeV);
294 point(0).setPy(piter->second->py_GeV);
295 point(0).setPz(piter->second->pz_GeV);
296 point(0).setR(piter->second->r_cm);
297 point(0).setPhi(piter->second->phi_rad);
298 point(0).setZ(piter->second->z_cm);
299 point(0).setT(piter->second->t_ns);
300 point(0).setTrack(piter->second->track_);
301 hddm_s::TrackIDList tid = point(0).addTrackIDs();
302 tid(0).setItrack(piter->second->trackID_);
303 }
304}
305
306int GlueXSensitiveDetectorGCAL::GetIdent(std::string div,
307 const G4VTouchable *touch)
308{
309 const HddsG4Builder* bldr = GlueXDetectorConstruction::GetBuilder();
310 std::map<std::string, std::vector<int> >::const_iterator iter;
311 std::map<std::string, std::vector<int> > *identifiers;
312 int max_depth = touch->GetHistoryDepth();
313 for (int depth = 0; depth < max_depth; ++depth) {
314 G4VPhysicalVolume *pvol = touch->GetVolume(depth);
315 G4LogicalVolume *lvol = pvol->GetLogicalVolume();
316 int volId = fVolumeTable[lvol];
317 if (volId == 0) {
318 volId = bldr->getVolumeId(lvol);
319 fVolumeTable[lvol] = volId;
320 }
321 identifiers = &Refsys::fIdentifierTable[volId];
322 if ((iter = identifiers->find(div)) != identifiers->end()) {
323 int copyNum = touch->GetCopyNumber(depth);
324 copyNum += (dynamic_cast<G4PVPlacement*>(pvol))? -1 : 0;
325 return iter->second[copyNum];
326 }
327 }
328 return -1;
329}