Bug Summary

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

Annotated Source Code

1//
2// GlueXSensitiveDetectorTPOL - class implementation
3//
4// author: richard.t.jones at uconn.edu
5// version: december 16, 2016
6
7#include "GlueXSensitiveDetectorTPOL.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#define sqr(x)((x)*(x)) ((x)*(x))
25
26// Cutoff on the number of allowed hits per sector
27int GlueXSensitiveDetectorTPOL::MAX_HITS = 100;
28
29// Minimum hit time difference for two hits on the same wedge
30double GlueXSensitiveDetectorTPOL::TWO_HIT_TIME_RESOL = 1000*ns;
31
32// Minimum energy deposition for a hit
33double GlueXSensitiveDetectorTPOL::THRESH_MEV = 0.050;
34
35int GlueXSensitiveDetectorTPOL::instanceCount = 0;
36G4Mutex GlueXSensitiveDetectorTPOL::fMutex = G4MUTEX_INITIALIZER{ { 0, 0, 0, 0, 0, 0, 0, { 0, 0 } } };
37
38GlueXSensitiveDetectorTPOL::GlueXSensitiveDetectorTPOL(const G4String& name)
39 : G4VSensitiveDetector(name),
40 fHitsMap(0), fPointsMap(0)
41{
42 collectionName.insert("TPOLWedgeHitsCollection");
43 collectionName.insert("TPOLPointsCollection");
44
45 // The rest of this only needs to happen once, the first time an object
46 // of this type is instantiated for this configuration of geometry and
47 // fields. If the geometry or fields change in such a way as to modify
48 // the drift-time properties of hits in the TPOL, you must delete all old
49 // objects of this class and create new ones.
50
51 G4AutoLock barrier(&fMutex);
52 if (instanceCount++ == 0) {
53 int runno = HddmOutput::getRunNo();
54 extern jana::JApplication *japp;
55 if (japp == 0) {
56 G4cerr(*G4cerr_p) << "Error in GlueXSensitiveDetector constructor - "
57 << "jana global DApplication object not set, "
58 << "cannot continue." << G4endlstd::endl;
59 exit(-1);
60 }
61 jana::JCalibration *jcalib = japp->GetJCalibration(runno);
Value stored to 'jcalib' during its initialization is never read
62 if (japp == 0) { // dummy
63 jcalib = 0;
64 G4cout(*G4cout_p) << "TPOL: ALL parameters loaded from ccdb" << G4endlstd::endl;
65 }
66 }
67}
68
69GlueXSensitiveDetectorTPOL::GlueXSensitiveDetectorTPOL(
70 const GlueXSensitiveDetectorTPOL &src)
71 : G4VSensitiveDetector(src),
72 fHitsMap(src.fHitsMap), fPointsMap(src.fPointsMap)
73{
74 G4AutoLock barrier(&fMutex);
75 ++instanceCount;
76}
77
78GlueXSensitiveDetectorTPOL &GlueXSensitiveDetectorTPOL::operator=(const
79 GlueXSensitiveDetectorTPOL &src)
80{
81 G4AutoLock barrier(&fMutex);
82 *(G4VSensitiveDetector*)this = src;
83 fHitsMap = src.fHitsMap;
84 fPointsMap = src.fPointsMap;
85 return *this;
86}
87
88GlueXSensitiveDetectorTPOL::~GlueXSensitiveDetectorTPOL()
89{
90 G4AutoLock barrier(&fMutex);
91 --instanceCount;
92}
93
94void GlueXSensitiveDetectorTPOL::Initialize(G4HCofThisEvent* hce)
95{
96 fHitsMap = new
97 GlueXHitsMapTPOLwedge(SensitiveDetectorName, collectionName[0]);
98 fPointsMap = new
99 GlueXHitsMapTPOLpoint(SensitiveDetectorName, collectionName[1]);
100 G4SDManager *sdm = G4SDManager::GetSDMpointer();
101 hce->AddHitsCollection(sdm->GetCollectionID(collectionName[0]), fHitsMap);
102 hce->AddHitsCollection(sdm->GetCollectionID(collectionName[1]), fPointsMap);
103}
104
105G4bool GlueXSensitiveDetectorTPOL::ProcessHits(G4Step* step,
106 G4TouchableHistory* ROhist)
107{
108 double dEsum = step->GetTotalEnergyDeposit();
109 if (dEsum == 0)
110 return false;
111
112 const G4ThreeVector &pin = step->GetPreStepPoint()->GetMomentum();
113 const G4ThreeVector &xin = step->GetPreStepPoint()->GetPosition();
114 const G4ThreeVector &xout = step->GetPostStepPoint()->GetPosition();
115 double Ein = step->GetPreStepPoint()->GetTotalEnergy();
116 double tin = step->GetPreStepPoint()->GetGlobalTime();
117 double tout = step->GetPostStepPoint()->GetGlobalTime();
118 G4ThreeVector x = (xin + xout) / 2;
119 G4ThreeVector dx = xout - xin;
120 double t = (tin + tout) / 2;
121
122 const G4VTouchable* touch = step->GetPreStepPoint()->GetTouchable();
123 const G4AffineTransform &local_from_global = touch->GetHistory()
124 ->GetTopTransform();
125 G4ThreeVector xlocal = local_from_global.TransformPoint(x);
126
127 // For particles that range out inside the active volume, the
128 // "out" time may sometimes be set to something enormously high.
129 // This screws up the hit. Check for this case here by looking
130 // at tout and making sure it is less than 1 second. If it's
131 // not, then just use tin for "t".
132
133 if (tout > 1.0*s)
134 t = tin;
135
136 double dr = dx.mag();
137 double dEdx = (dr > 1e-3*cm)? dEsum/dr : 0;
138
139 // Post the hit to the points list in the
140 // order of appearance in the event simulation.
141
142 int ring = 0; // GetIdent("ring", touch);
143 int sector = GetIdent("sector", touch);
144 G4Track *track = step->GetTrack();
145 G4int trackID = track->GetTrackID();
146 int pdgtype = track->GetDynamicParticle()->GetPDGcode();
147 int g3type = GlueXPrimaryGeneratorAction::ConvertPdgToGeant3(pdgtype);
148 GlueXUserTrackInformation *trackinfo = (GlueXUserTrackInformation*)
149 track->GetUserInformation();
150 int itrack = trackinfo->GetGlueXTrackID();
151 if (trackinfo->GetGlueXHistory() == 0 && itrack > 0) {
152 G4int key = fPointsMap->entries();
153 GlueXHitTPOLpoint* lastPoint = (*fPointsMap)[key - 1];
154 if (lastPoint == 0 || lastPoint->track_ != trackID ||
155 fabs(lastPoint->t_ns - t/ns) > 0.1 ||
156 (fabs(lastPoint->r_cm - x.perp()/cm) > 0.1 &&
157 fabs(lastPoint->phi_rad - x.phi()) > 0.1) )
158 {
159 GlueXHitTPOLpoint newPoint;
160 newPoint.ptype_G3 = g3type;
161 newPoint.track_ = trackID;
162 newPoint.trackID_ = itrack;
163 newPoint.primary_ = (track->GetParentID() == 0);
164 newPoint.phi_rad = x.phi();
165 newPoint.r_cm = x.perp()/cm;
166 newPoint.t_ns = t/ns;
167 newPoint.px_GeV = pin[0]/GeV;
168 newPoint.py_GeV = pin[1]/GeV;
169 newPoint.pz_GeV = pin[2]/GeV;
170 newPoint.E_GeV = Ein/GeV;
171 newPoint.dEdx_GeV_cm = dEdx/(GeV/cm);
172 fPointsMap->add(key, newPoint);
173 }
174 }
175
176 // Post the hit to the hits map, ordered by sector,ring index
177
178 if (dEsum > 0) {
179 int key = GlueXHitTPOLwedge::GetKey(sector, ring);
180 GlueXHitTPOLwedge *wedge = (*fHitsMap)[key];
181 if (wedge == 0) {
182 GlueXHitTPOLwedge newwedge(sector, ring);
183 fHitsMap->add(key, newwedge);
184 wedge = (*fHitsMap)[key];
185 }
186
187 // Add the hit to the hits vector, maintaining strict time ordering
188
189 int merge_hit = 0;
190 std::vector<GlueXHitTPOLwedge::hitinfo_t>::iterator hiter;
191 for (hiter = wedge->hits.begin(); hiter != wedge->hits.end(); ++hiter) {
192 if (fabs(hiter->t_ns*ns - t) < TWO_HIT_TIME_RESOL) {
193 merge_hit = 1;
194 break;
195 }
196 else if (hiter->t_ns*ns > t) {
197 break;
198 }
199 }
200 if (merge_hit) {
201 // Use the time from the earlier hit but add the charge
202 hiter->dE_MeV += dEsum/MeV;
203 if (hiter->t_ns*ns > t) {
204 hiter->t_ns = t/ns;
205 hiter->itrack_ = itrack;
206 hiter->ptype_G3 = g3type;
207 hiter->t0_ns = t/ns;
208 hiter->r_cm = x.perp()/cm;
209 }
210 }
211 else {
212 // create new hit
213 hiter = wedge->hits.insert(hiter, GlueXHitTPOLwedge::hitinfo_t());
214 hiter->dE_MeV = dEsum/MeV;
215 hiter->t_ns = t/ns;
216 hiter->itrack_ = itrack;
217 hiter->ptype_G3 = g3type;
218 hiter->t0_ns = t/ns;
219 hiter->r_cm = x.perp()/cm;
220 }
221 }
222 return true;
223}
224
225void GlueXSensitiveDetectorTPOL::EndOfEvent(G4HCofThisEvent*)
226{
227 std::map<int,GlueXHitTPOLwedge*> *wedges = fHitsMap->GetMap();
228 std::map<int,GlueXHitTPOLpoint*> *points = fPointsMap->GetMap();
229 if (wedges->size() == 0 && points->size() == 0)
230 return;
231 std::map<int,GlueXHitTPOLwedge*>::iterator siter;
232 std::map<int,GlueXHitTPOLpoint*>::iterator piter;
233
234 if (verboseLevel > 1) {
235 G4cout(*G4cout_p) << G4endlstd::endl
236 << "--------> Hits Collection: in this event there are "
237 << wedges->size() << " wedges with hits in the TPOL: "
238 << G4endlstd::endl;
239 for (siter = wedges->begin(); siter != wedges->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 TPOL: "
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) << "GlueXSensitiveDetectorTPOL::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.getTripletPolarimeters().size() == 0)
268 hitview.addTripletPolarimeters();
269 hddm_s::TripletPolarimeter &polarimeter = hitview.getTripletPolarimeter();
270
271 // Collect and output the tpolTruthHits
272 for (siter = wedges->begin(); siter != wedges->end(); ++siter) {
273 std::vector<GlueXHitTPOLwedge::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].dE_MeV <= THRESH_MEV) {
277 hits.erase(hits.begin() + ih);
278 --ih;
279 }
280 }
281 if (hits.size() > 0) {
282 hddm_s::TpolSectorList wedge = polarimeter.addTpolSectors(1);
283 wedge(0).setSector(siter->second->sector_);
284 wedge(0).setRing(siter->second->ring_);
285 int hitscount = hits.size();
286 if (hitscount > MAX_HITS) {
287 hitscount = MAX_HITS;
288 G4cerr(*G4cerr_p) << "GlueXSensitiveDetectorTPOL::EndOfEvent warning: "
289 << "counter max hit count " << MAX_HITS << " exceeded, "
290 << hits.size() - hitscount << " hits discarded."
291 << G4endlstd::endl;
292 }
293 for (int ih=0; ih < hitscount; ++ih) {
294 hddm_s::TpolTruthHitList thit = wedge(0).addTpolTruthHits(1);
295 thit(0).setDE(hits[ih].dE_MeV*MeV/GeV);
296 thit(0).setT(hits[ih].t_ns);
297 thit(0).setItrack(hits[ih].itrack_);
298 thit(0).setPtype(hits[ih].ptype_G3);
299 }
300 }
301 }
302
303 // Collect and output the tpolTruthPoints
304 for (piter = points->begin(); piter != points->end(); ++piter) {
305 hddm_s::TpolTruthPointList point = polarimeter.addTpolTruthPoints(1);
306 point(0).setE(piter->second->E_GeV);
307 point(0).setDEdx(piter->second->dEdx_GeV_cm);
308 point(0).setPrimary(piter->second->primary_);
309 point(0).setPtype(piter->second->ptype_G3);
310 point(0).setPx(piter->second->px_GeV);
311 point(0).setPy(piter->second->py_GeV);
312 point(0).setPz(piter->second->pz_GeV);
313 point(0).setPhi(piter->second->phi_rad);
314 point(0).setR(piter->second->r_cm);
315 point(0).setT(piter->second->t_ns);
316 point(0).setTrack(piter->second->track_);
317 hddm_s::TrackIDList tid = point(0).addTrackIDs();
318 tid(0).setItrack(piter->second->trackID_);
319 }
320}
321
322int GlueXSensitiveDetectorTPOL::GetIdent(std::string div,
323 const G4VTouchable *touch)
324{
325 const HddsG4Builder* bldr = GlueXDetectorConstruction::GetBuilder();
326 std::map<std::string, std::vector<int> >::const_iterator iter;
327 std::map<std::string, std::vector<int> > *identifiers;
328 int max_depth = touch->GetHistoryDepth();
329 for (int depth = 0; depth < max_depth; ++depth) {
330 G4VPhysicalVolume *pvol = touch->GetVolume(depth);
331 G4LogicalVolume *lvol = pvol->GetLogicalVolume();
332 int volId = fVolumeTable[lvol];
333 if (volId == 0) {
334 volId = bldr->getVolumeId(lvol);
335 fVolumeTable[lvol] = volId;
336 }
337 identifiers = &Refsys::fIdentifierTable[volId];
338 if ((iter = identifiers->find(div)) != identifiers->end()) {
339 int copyNum = touch->GetCopyNumber(depth);
340 copyNum += (dynamic_cast<G4PVPlacement*>(pvol))? -1 : 0;
341 return iter->second[copyNum];
342 }
343 }
344 return -1;
345}