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 total number of allowed hits
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 if ((int)wedge->hits.size() < MAX_HITS) {
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 else {
222 G4cerr(*G4cerr_p) << "GlueXSensitiveDetectorTPOL::ProcessHits error: "
223 << "max hit count " << MAX_HITS << " exceeded, truncating!"
224 << G4endlstd::endl;
225 }
226 }
227 return true;
228}
229
230void GlueXSensitiveDetectorTPOL::EndOfEvent(G4HCofThisEvent*)
231{
232 std::map<int,GlueXHitTPOLwedge*> *wedges = fHitsMap->GetMap();
233 std::map<int,GlueXHitTPOLpoint*> *points = fPointsMap->GetMap();
234 if (wedges->size() == 0 && points->size() == 0)
235 return;
236 std::map<int,GlueXHitTPOLwedge*>::iterator siter;
237 std::map<int,GlueXHitTPOLpoint*>::iterator piter;
238
239 if (verboseLevel > 1) {
240 G4cout(*G4cout_p) << G4endlstd::endl
241 << "--------> Hits Collection: in this event there are "
242 << wedges->size() << " wedges with hits in the TPOL: "
243 << G4endlstd::endl;
244 for (siter = wedges->begin(); siter != wedges->end(); ++siter)
245 siter->second->Print();
246
247 G4cout(*G4cout_p) << G4endlstd::endl
248 << "--------> Hits Collection: in this event there are "
249 << points->size() << " truth points in the TPOL: "
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) << "GlueXSensitiveDetectorTPOL::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.getTripletPolarimeters().size() == 0)
273 hitview.addTripletPolarimeters();
274 hddm_s::TripletPolarimeter &polarimeter = hitview.getTripletPolarimeter();
275
276 // Collect and output the tpolTruthHits
277 for (siter = wedges->begin(); siter != wedges->end(); ++siter) {
278 std::vector<GlueXHitTPOLwedge::hitinfo_t> &hits = siter->second->hits;
279 // apply a pulse height threshold cut
280 for (unsigned int ih=0; ih < hits.size(); ++ih) {
281 if (hits[ih].dE_MeV <= THRESH_MEV) {
282 hits.erase(hits.begin() + ih);
283 --ih;
284 }
285 }
286 if (hits.size() > 0) {
287 hddm_s::TpolSectorList wedge = polarimeter.addTpolSectors(1);
288 wedge(0).setSector(siter->second->sector_);
289 wedge(0).setRing(siter->second->ring_);
290 for (int ih=0; ih < (int)hits.size(); ++ih) {
291 hddm_s::TpolTruthHitList thit = wedge(0).addTpolTruthHits(1);
292 thit(0).setDE(hits[ih].dE_MeV*MeV/GeV);
293 thit(0).setT(hits[ih].t_ns);
294 thit(0).setItrack(hits[ih].itrack_);
295 thit(0).setPtype(hits[ih].ptype_G3);
296 }
297 }
298 }
299
300 // Collect and output the tpolTruthPoints
301 for (piter = points->begin(); piter != points->end(); ++piter) {
302 hddm_s::TpolTruthPointList point = polarimeter.addTpolTruthPoints(1);
303 point(0).setE(piter->second->E_GeV);
304 point(0).setDEdx(piter->second->dEdx_GeV_cm);
305 point(0).setPrimary(piter->second->primary_);
306 point(0).setPtype(piter->second->ptype_G3);
307 point(0).setPx(piter->second->px_GeV);
308 point(0).setPy(piter->second->py_GeV);
309 point(0).setPz(piter->second->pz_GeV);
310 point(0).setPhi(piter->second->phi_rad);
311 point(0).setR(piter->second->r_cm);
312 point(0).setT(piter->second->t_ns);
313 point(0).setTrack(piter->second->track_);
314 hddm_s::TrackIDList tid = point(0).addTrackIDs();
315 tid(0).setItrack(piter->second->trackID_);
316 }
317}
318
319int GlueXSensitiveDetectorTPOL::GetIdent(std::string div,
320 const G4VTouchable *touch)
321{
322 const HddsG4Builder* bldr = GlueXDetectorConstruction::GetBuilder();
323 std::map<std::string, std::vector<int> >::const_iterator iter;
324 std::map<std::string, std::vector<int> > *identifiers;
325 int max_depth = touch->GetHistoryDepth();
326 for (int depth = 0; depth < max_depth; ++depth) {
327 G4VPhysicalVolume *pvol = touch->GetVolume(depth);
328 G4LogicalVolume *lvol = pvol->GetLogicalVolume();
329 int volId = fVolumeTable[lvol];
330 if (volId == 0) {
331 volId = bldr->getVolumeId(lvol);
332 fVolumeTable[lvol] = volId;
333 }
334 identifiers = &Refsys::fIdentifierTable[volId];
335 if ((iter = identifiers->find(div)) != identifiers->end()) {
336 int copyNum = touch->GetCopyNumber(depth);
337 copyNum += (dynamic_cast<G4PVPlacement*>(pvol))? -1 : 0;
338 return iter->second[copyNum];
339 }
340 }
341 return -1;
342}