Bug Summary

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

Annotated Source Code

1//
2// GlueXSensitiveDetectorCCAL - class implementation
3//
4// author: richard.t.jones at uconn.edu
5// version: october 28, 2016
6
7#include "GlueXSensitiveDetectorCCAL.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 block
25int GlueXSensitiveDetectorCCAL::MAX_HITS = 100;
26
27// Geometry constants for the CCal
28//int GlueXSensitiveDetectorCCAL::CENTRAL_ROW = 8;
29//int GlueXSensitiveDetectorCCAL::CENTRAL_COLUMN = 8;
30double GlueXSensitiveDetectorCCAL::WIDTH_OF_BLOCK = 2.0*cm;
31double GlueXSensitiveDetectorCCAL::LENGTH_OF_BLOCK = 20.0*cm;
32
33// Light propagation parameters in Compton calorimeter
34double GlueXSensitiveDetectorCCAL::ATTENUATION_LENGTH = 200.*cm;
35double GlueXSensitiveDetectorCCAL::C_EFFECTIVE = 13.*cm/ns;
36
37// Minimum hit time difference for two hits on the same block
38double GlueXSensitiveDetectorCCAL::TWO_HIT_TIME_RESOL = 75.*ns;
39
40// Minimum energy deposition for a hit
41double GlueXSensitiveDetectorCCAL::THRESH_MEV = 5.;
42
43int GlueXSensitiveDetectorCCAL::instanceCount = 0;
44G4Mutex GlueXSensitiveDetectorCCAL::fMutex = G4MUTEX_INITIALIZER{ { 0, 0, 0, 0, 0, 0, 0, { 0, 0 } } };
45
46GlueXSensitiveDetectorCCAL::GlueXSensitiveDetectorCCAL(const G4String& name)
47 : G4VSensitiveDetector(name),
48 fBlocksMap(0), fPointsMap(0)
49{
50 collectionName.insert("CCALBlockHitsCollection");
51 collectionName.insert("CCALPointsCollection");
52
53 // The rest of this only needs to happen once, the first time an object
54 // of this type is instantiated for this configuration of geometry and
55 // fields. If the geometry or fields change in such a way as to modify
56 // the drift-time properties of hits in the CCAL, you must delete all old
57 // objects of this class and create new ones.
58
59 G4AutoLock barrier(&fMutex);
60 if (instanceCount++ == 0) {
61 int runno = HddmOutput::getRunNo();
62 extern jana::JApplication *japp;
63 if (japp == 0) {
64 G4cerr(*G4cerr_p) << "Error in GlueXSensitiveDetector constructor - "
65 << "jana global DApplication object not set, "
66 << "cannot continue." << G4endlstd::endl;
67 exit(-1);
68 }
69 jana::JCalibration *jcalib = japp->GetJCalibration(runno);
Value stored to 'jcalib' during its initialization is never read
70 if (japp == 0) { // dummy
71 jcalib = 0;
72 G4cout(*G4cout_p) << "CCAL: ALL parameters loaded from ccdb" << G4endlstd::endl;
73 }
74 }
75}
76
77GlueXSensitiveDetectorCCAL::GlueXSensitiveDetectorCCAL(
78 const GlueXSensitiveDetectorCCAL &src)
79 : G4VSensitiveDetector(src),
80 fBlocksMap(src.fBlocksMap), fPointsMap(src.fPointsMap)
81{
82 G4AutoLock barrier(&fMutex);
83 ++instanceCount;
84}
85
86GlueXSensitiveDetectorCCAL &GlueXSensitiveDetectorCCAL::operator=(const
87 GlueXSensitiveDetectorCCAL &src)
88{
89 G4AutoLock barrier(&fMutex);
90 *(G4VSensitiveDetector*)this = src;
91 fBlocksMap = src.fBlocksMap;
92 fPointsMap = src.fPointsMap;
93 return *this;
94}
95
96GlueXSensitiveDetectorCCAL::~GlueXSensitiveDetectorCCAL()
97{
98 G4AutoLock barrier(&fMutex);
99 --instanceCount;
100}
101
102void GlueXSensitiveDetectorCCAL::Initialize(G4HCofThisEvent* hce)
103{
104 fBlocksMap = new
105 GlueXHitsMapCCALblock(SensitiveDetectorName, collectionName[0]);
106 fPointsMap = new
107 GlueXHitsMapCCALpoint(SensitiveDetectorName, collectionName[1]);
108 G4SDManager *sdm = G4SDManager::GetSDMpointer();
109 hce->AddHitsCollection(sdm->GetCollectionID(collectionName[0]), fBlocksMap);
110 hce->AddHitsCollection(sdm->GetCollectionID(collectionName[1]), fPointsMap);
111}
112
113G4bool GlueXSensitiveDetectorCCAL::ProcessHits(G4Step* step,
114 G4TouchableHistory* ROhist)
115{
116 double dEsum = step->GetTotalEnergyDeposit();
117 if (dEsum == 0)
118 return false;
119
120 const G4ThreeVector &pin = step->GetPreStepPoint()->GetMomentum();
121 const G4ThreeVector &xin = step->GetPreStepPoint()->GetPosition();
122 const G4ThreeVector &xout = step->GetPostStepPoint()->GetPosition();
123 double Ein = step->GetPreStepPoint()->GetTotalEnergy();
124 double tin = step->GetPreStepPoint()->GetGlobalTime();
125 double tout = step->GetPostStepPoint()->GetGlobalTime();
126 G4ThreeVector x = (xin + xout) / 2;
127 double t = (tin + tout) / 2;
128
129 const G4VTouchable* touch = step->GetPreStepPoint()->GetTouchable();
130 const G4AffineTransform &local_from_global = touch->GetHistory()
131 ->GetTopTransform();
132 G4ThreeVector xlocal = local_from_global.TransformPoint(x);
133
134 // For particles that range out inside the active volume, the
135 // "out" time may sometimes be set to something enormously high.
136 // This screws up the hit. Check for this case here by looking
137 // at tout and making sure it is less than 1 second. If it's
138 // not, then just use tin for "t".
139
140 if (tout > 1.0*s)
141 t = tin;
142
143 // Post the hit to the points list in the
144 // order of appearance in the event simulation.
145
146 G4Track *track = step->GetTrack();
147 G4int trackID = track->GetTrackID();
148 GlueXUserTrackInformation *trackinfo = (GlueXUserTrackInformation*)
149 track->GetUserInformation();
150 int itrack = trackinfo->GetGlueXTrackID();
151 if (track->GetCurrentStepNumber() == 1)
152 trackinfo->SetGlueXHistory(4);
153 if (trackinfo->GetGlueXHistory() == 0 &&
154 xin.dot(pin) > 0 && Ein/MeV > THRESH_MEV)
155 {
156 G4int key = fPointsMap->entries();
157 GlueXHitCCALpoint* lastPoint = (*fPointsMap)[key - 1];
158 // Limit ccal truthPoints to one per track
159 if (lastPoint == 0 || lastPoint->track_ != trackID) {
160 int pdgtype = track->GetDynamicParticle()->GetPDGcode();
161 int g3type = GlueXPrimaryGeneratorAction::ConvertPdgToGeant3(pdgtype);
162 GlueXHitCCALpoint newPoint;
163 newPoint.ptype_G3 = g3type;
164 newPoint.track_ = trackID;
165 newPoint.trackID_ = itrack;
166 newPoint.primary_ = (track->GetParentID() == 0);
167 newPoint.t_ns = t/ns;
168 newPoint.x_cm = xin[0]/cm;
169 newPoint.y_cm = xin[1]/cm;
170 newPoint.z_cm = xin[2]/cm;
171 newPoint.px_GeV = pin[0]/GeV;
172 newPoint.py_GeV = pin[1]/GeV;
173 newPoint.pz_GeV = pin[2]/GeV;
174 newPoint.E_GeV = Ein/GeV;
175 fPointsMap->add(key, newPoint);
176 }
177 }
178
179 // Post the hit to the hits map, ordered by sector index
180
181 if (dEsum > 0) {
182 int column = GetIdent("column", touch);
183 int row = GetIdent("row", touch);
184 int key = GlueXHitCCALblock::GetKey(column, row);
185 GlueXHitCCALblock *block = (*fBlocksMap)[key];
186 if (block == 0) {
187 GlueXHitCCALblock newblock(column, row);
188 fBlocksMap->add(key, newblock);
189 block = (*fBlocksMap)[key];
190 }
191 double dist = 0.5 * LENGTH_OF_BLOCK - xlocal[2];
192 double dEcorr = dEsum * exp(-dist / ATTENUATION_LENGTH);
193 double tcorr = t + dist / C_EFFECTIVE;
194
195 // Add the hit to the hits vector, maintaining strict time ordering
196
197 int merge_hit = 0;
198 std::vector<GlueXHitCCALblock::hitinfo_t>::iterator hiter;
199 for (hiter = block->hits.begin(); hiter != block->hits.end(); ++hiter) {
200 if (fabs(hiter->t_ns*ns - tcorr) < TWO_HIT_TIME_RESOL) {
201 merge_hit = 1;
202 break;
203 }
204 else if (hiter->t_ns*ns > t) {
205 break;
206 }
207 }
208 if (merge_hit) {
209 // Use the time from the earlier hit but add the energy deposition
210 hiter->E_GeV += dEcorr/GeV;
211 if (hiter->t_ns*ns > tcorr) {
212 hiter->t_ns = tcorr/ns;
213 }
214 }
215 else {
216 // create new hit
217 hiter = block->hits.insert(hiter, GlueXHitCCALblock::hitinfo_t());
218 hiter->E_GeV = dEcorr/GeV;
219 hiter->t_ns = tcorr/ns;
220 }
221 }
222 return true;
223}
224
225void GlueXSensitiveDetectorCCAL::EndOfEvent(G4HCofThisEvent*)
226{
227 std::map<int,GlueXHitCCALblock*> *blocks = fBlocksMap->GetMap();
228 std::map<int,GlueXHitCCALpoint*> *points = fPointsMap->GetMap();
229 if (blocks->size() == 0 && points->size() == 0)
230 return;
231 std::map<int,GlueXHitCCALblock*>::iterator biter;
232 std::map<int,GlueXHitCCALpoint*>::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 CCAL: "
238 << G4endlstd::endl;
239 for (biter = blocks->begin(); biter != blocks->end(); ++biter)
240 biter->second->Print();
241
242 G4cout(*G4cout_p) << G4endlstd::endl
243 << "--------> Hits Collection: in this event there are "
244 << points->size() << " truth showers in the CCAL: "
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) << "GlueXSensitiveDetectorCCAL::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.getComptonEMcals().size() == 0)
268 hitview.addComptonEMcals();
269 hddm_s::ComptonEMcal &comptonEMcal = hitview.getComptonEMcal();
270
271 // Collect and output the ccalTruthHits
272 for (biter = blocks->begin(); biter != blocks->end(); ++biter) {
273 std::vector<GlueXHitCCALblock::hitinfo_t> &hits = biter->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 <= THRESH_MEV/1e3) {
277 hits.erase(hits.begin() + ih);
278 --ih;
279 }
280 }
281
282 if (hits.size() > 0) {
283 hddm_s::CcalBlockList block = comptonEMcal.addCcalBlocks(1);
284 block(0).setColumn(biter->second->column_);
285 block(0).setRow(biter->second->row_);
286 int hitscount = hits.size();
287 if (hitscount > MAX_HITS) {
288 hitscount = MAX_HITS;
289 G4cerr(*G4cerr_p) << "GlueXSensitiveDetectorCCAL::EndOfEvent warning: "
290 << "max hit count " << MAX_HITS << " exceeded, "
291 << hits.size() - hitscount << " hits discarded."
292 << G4endlstd::endl;
293 }
294 for (int ih=0; ih < hitscount; ++ih) {
295 hddm_s::CcalTruthHitList thit = block(0).addCcalTruthHits(1);
296 thit(0).setE(hits[ih].E_GeV);
297 thit(0).setT(hits[ih].t_ns);
298 }
299 }
300 }
301
302 // Collect and output the fcalTruthShowers
303 for (piter = points->begin(); piter != points->end(); ++piter) {
304 hddm_s::CcalTruthShowerList point = comptonEMcal.addCcalTruthShowers(1);
305 point(0).setE(piter->second->E_GeV);
306 point(0).setPrimary(piter->second->primary_);
307 point(0).setPtype(piter->second->ptype_G3);
308 point(0).setPx(piter->second->px_GeV);
309 point(0).setPy(piter->second->py_GeV);
310 point(0).setPz(piter->second->pz_GeV);
311 point(0).setX(piter->second->x_cm);
312 point(0).setY(piter->second->y_cm);
313 point(0).setZ(piter->second->z_cm);
314 point(0).setT(piter->second->t_ns);
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 GlueXSensitiveDetectorCCAL::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}