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