Bug Summary

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

Annotated Source Code

1//
2// GlueXSensitiveDetectorFMWPC - class implementation
3//
4// author: richard.t.jones at uconn.edu
5// version: november 29, 2016
6
7#include "GlueXSensitiveDetectorFMWPC.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 GlueXSensitiveDetectorFMWPC::MAX_HITS = 100;
26
27// Minimum hit time difference for two hits on the same wire
28double GlueXSensitiveDetectorFMWPC::TWO_HIT_TIME_RESOL = 400*ns;
29
30// Minimum photoelectron count for a hit
31double GlueXSensitiveDetectorFMWPC::THRESH_KEV = 0.;
32
33// Coordinate of wire 0, transverse to wire direction
34double GlueXSensitiveDetectorFMWPC::WIRE_OFFSET = -73*cm;
35
36// Minimum photoelectron count for a hit
37double GlueXSensitiveDetectorFMWPC::WIRE_PITCH = 1.0*cm;
38
39int GlueXSensitiveDetectorFMWPC::instanceCount = 0;
40G4Mutex GlueXSensitiveDetectorFMWPC::fMutex = G4MUTEX_INITIALIZER{ { 0, 0, 0, 0, 0, 0, 0, { 0, 0 } } };
41
42GlueXSensitiveDetectorFMWPC::GlueXSensitiveDetectorFMWPC(const G4String& name)
43 : G4VSensitiveDetector(name),
44 fWireHitsMap(0), fPointsMap(0)
45{
46 collectionName.insert("FMWPCWireHitsCollection");
47 collectionName.insert("FMWPCPointsCollection");
48
49 // The rest of this only needs to happen once, the first time an object
50 // of this type is instantiated for this configuration of geometry and
51 // fields. If the geometry or fields change in such a way as to modify
52 // the drift-time properties of hits in the FMWPC, you must delete all old
53 // objects of this class and create new ones.
54
55 G4AutoLock wirerier(&fMutex);
56 if (instanceCount++ == 0) {
57 int runno = HddmOutput::getRunNo();
58 extern jana::JApplication *japp;
59 if (japp == 0) {
60 G4cerr(*G4cerr_p) << "Error in GlueXSensitiveDetector constructor - "
61 << "jana global DApplication object not set, "
62 << "cannot continue." << G4endlstd::endl;
63 exit(-1);
64 }
65 jana::JCalibration *jcalib = japp->GetJCalibration(runno);
Value stored to 'jcalib' during its initialization is never read
66 if (japp == 0) { // dummy
67 jcalib = 0;
68 G4cout(*G4cout_p) << "FMWPC: ALL parameters loaded from ccdb" << G4endlstd::endl;
69 }
70 }
71}
72
73GlueXSensitiveDetectorFMWPC::GlueXSensitiveDetectorFMWPC(
74 const GlueXSensitiveDetectorFMWPC &src)
75 : G4VSensitiveDetector(src),
76 fWireHitsMap(src.fWireHitsMap), fPointsMap(src.fPointsMap)
77{
78 G4AutoLock wirerier(&fMutex);
79 ++instanceCount;
80}
81
82GlueXSensitiveDetectorFMWPC &GlueXSensitiveDetectorFMWPC::operator=(const
83 GlueXSensitiveDetectorFMWPC &src)
84{
85 G4AutoLock wirerier(&fMutex);
86 *(G4VSensitiveDetector*)this = src;
87 fWireHitsMap = src.fWireHitsMap;
88 fPointsMap = src.fPointsMap;
89 return *this;
90}
91
92GlueXSensitiveDetectorFMWPC::~GlueXSensitiveDetectorFMWPC()
93{
94 G4AutoLock wirerier(&fMutex);
95 --instanceCount;
96}
97
98void GlueXSensitiveDetectorFMWPC::Initialize(G4HCofThisEvent* hce)
99{
100 fWireHitsMap = new
101 GlueXHitsMapFMWPCwire(SensitiveDetectorName, collectionName[0]);
102 fPointsMap = new
103 GlueXHitsMapFMWPCpoint(SensitiveDetectorName, collectionName[1]);
104 G4SDManager *sdm = G4SDManager::GetSDMpointer();
105 hce->AddHitsCollection(sdm->GetCollectionID(collectionName[0]), fWireHitsMap);
106 hce->AddHitsCollection(sdm->GetCollectionID(collectionName[1]), fPointsMap);
107}
108
109G4bool GlueXSensitiveDetectorFMWPC::ProcessHits(G4Step* step,
110 G4TouchableHistory* ROhist)
111{
112 double dEsum = step->GetTotalEnergyDeposit();
113 if (dEsum == 0)
114 return false;
115
116 const G4ThreeVector &pin = step->GetPreStepPoint()->GetMomentum();
117 const G4ThreeVector &xin = step->GetPreStepPoint()->GetPosition();
118 const G4ThreeVector &xout = step->GetPostStepPoint()->GetPosition();
119 double Ein = step->GetPreStepPoint()->GetTotalEnergy();
120 double tin = step->GetPreStepPoint()->GetGlobalTime();
121 double tout = step->GetPostStepPoint()->GetGlobalTime();
122 G4ThreeVector x = (xin + xout) / 2;
123 G4ThreeVector dx = xout - xin;
124 double t = (tin + tout) / 2;
125
126 const G4VTouchable* touch = step->GetPreStepPoint()->GetTouchable();
127 const G4AffineTransform &local_from_global = touch->GetHistory()
128 ->GetTopTransform();
129 G4ThreeVector xlocal = local_from_global.TransformPoint(x);
130
131 // For particles that range out inside the active volume, the
132 // "out" time may sometimes be set to something enormously high.
133 // This screws up the hit. Check for this case here by looking
134 // at tout and making sure it is less than 1 second. If it's
135 // not, then just use tin for "t".
136
137 if (tout > 1.0*s)
138 t = tin;
139
140 // Post the hit to the points list in the
141 // order of appearance in the event simulation.
142
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 && xin.dot(pin) > 0) {
149 G4int key = fPointsMap->entries();
150 GlueXHitFMWPCpoint* lastPoint = (*fPointsMap)[key - 1];
151 if (lastPoint == 0 || lastPoint->track_ != trackID ||
152 fabs(lastPoint->t_ns - t/ns) > 0.1 ||
153 fabs(lastPoint->x_cm - x[0]/cm) > 2. ||
154 fabs(lastPoint->y_cm - x[1]/cm) > 2. ||
155 fabs(lastPoint->z_cm - x[2]/cm) > 2.)
156 {
157 int pdgtype = track->GetDynamicParticle()->GetPDGcode();
158 int g3type = GlueXPrimaryGeneratorAction::ConvertPdgToGeant3(pdgtype);
159 GlueXHitFMWPCpoint newPoint;
160 newPoint.ptype_G3 = g3type;
161 newPoint.track_ = trackID;
162 newPoint.trackID_ = itrack;
163 newPoint.primary_ = (track->GetParentID() == 0);
164 newPoint.t_ns = t/ns;
165 newPoint.x_cm = x[0]/cm;
166 newPoint.y_cm = x[1]/cm;
167 newPoint.z_cm = x[2]/cm;
168 newPoint.px_GeV = pin[0]/GeV;
169 newPoint.py_GeV = pin[1]/GeV;
170 newPoint.pz_GeV = pin[2]/GeV;
171 newPoint.E_GeV = Ein/GeV;
172 fPointsMap->add(key, newPoint);
173 }
174 }
175
176 // Post the hit to the hits map, ordered by plane,wire,end index
177
178 if (dEsum > 0) {
179
180 // HDDS geometry does not include wires nor plane rotations. Assume 1 cm
181 // wire spacing with wires at 0.5cm on either side of the beamline at
182 // x,y = 0,0. Also assume odd number planes have wires in the vertical
183 // direction and even numbered planes have wires in the horizontal direction.
184 // Vertical wires start with wire 1 at x=-71.5 and wire 144 at x=+71.5
185 // (the gas volume ends at x=+/-72.0). Horizontal wires start with wire 1
186 // at y=-71.5 (i.e. closest to the ground) and wire 144 at y=+71.5
187 // (i.e. closest to the sky).
188
189 int layer = GetIdent("layer", touch);
190 int wire = 0;
191 if (layer % 2 != 0) {
192 // Vertical wires
193 wire = floor(x[0] - WIRE_OFFSET)/WIRE_PITCH;
194 }
195 else {
196 // Horizontal wires
197 wire = floor(x[1] - WIRE_OFFSET)/WIRE_PITCH;
198 }
199 if (layer == 6) {
200 // Vertical wires
201 wire = floor(x[0] - WIRE_OFFSET)/WIRE_PITCH;
202 }
203
204 //cout<<"MWPC: layer/wire = "<<layer<<" / "<<wire<<endl;
205
206 if (wire < 1 || wire > 144)
207 return false;
208
209 int key = GlueXHitFMWPCwire::GetKey(layer, wire);
210 GlueXHitFMWPCwire *counter = (*fWireHitsMap)[key];
211 if (counter == 0) {
212 GlueXHitFMWPCwire newwire(layer, wire);
213 fWireHitsMap->add(key, newwire);
214 counter = (*fWireHitsMap)[key];
215 }
216
217 // Add the hit to the hits vector, maintaining strict time ordering
218
219 int merge_hit = 0;
220 std::vector<GlueXHitFMWPCwire::hitinfo_t>::iterator hiter;
221 for (hiter = counter->hits.begin(); hiter != counter->hits.end(); ++hiter) {
222 if (fabs(hiter->t_ns*ns - t) < TWO_HIT_TIME_RESOL) {
223 merge_hit = 1;
224 break;
225 }
226 else if (hiter->t_ns*ns > t) {
227 break;
228 }
229 }
230 if (merge_hit) {
231 // Use the time from the earlier hit but add the charge
232 hiter->dE_keV += dEsum/keV;
233 hiter->dx_cm += dx.mag()/cm;
234 if (hiter->t_ns*ns > t) {
235 hiter->t_ns = t/ns;
236 }
237 }
238 else if ((int)counter->hits.size() < MAX_HITS) {
239 // create new hit
240 hiter = counter->hits.insert(hiter, GlueXHitFMWPCwire::hitinfo_t());
241 hiter->dE_keV = dEsum/keV;
242 hiter->dx_cm = dx.mag()/cm;
243 hiter->t_ns = t/ns;
244 }
245 else {
246 G4cerr(*G4cerr_p) << "GlueXSensitiveDetectorFMWPC::ProcessHits error: "
247 << "max hit count " << MAX_HITS
248 << " exceeded, truncating!"
249 << G4endlstd::endl;
250 }
251 }
252 return true;
253}
254
255void GlueXSensitiveDetectorFMWPC::EndOfEvent(G4HCofThisEvent*)
256{
257 std::map<int,GlueXHitFMWPCwire*> *wires = fWireHitsMap->GetMap();
258 std::map<int,GlueXHitFMWPCpoint*> *points = fPointsMap->GetMap();
259 if (wires->size() == 0 && points->size() == 0)
260 return;
261 std::map<int,GlueXHitFMWPCwire*>::iterator siter;
262 std::map<int,GlueXHitFMWPCpoint*>::iterator piter;
263
264 if (verboseLevel > 1) {
265 G4cout(*G4cout_p) << G4endlstd::endl
266 << "--------> Hits Collection: in this event there are "
267 << wires->size() << " wires with hits in the FMWPC: "
268 << G4endlstd::endl;
269 for (siter = wires->begin(); siter != wires->end(); ++siter)
270 siter->second->Print();
271
272 G4cout(*G4cout_p) << G4endlstd::endl
273 << "--------> Hits Collection: in this event there are "
274 << points->size() << " truth points in the FWMPC: "
275 << G4endlstd::endl;
276 for (piter = points->begin(); piter != points->end(); ++piter)
277 piter->second->Print();
278 }
279
280 // pack hits into ouptut hddm record
281
282 G4EventManager* mgr = G4EventManager::GetEventManager();
283 G4VUserEventInformation* info = mgr->GetUserInformation();
284 hddm_s::HDDM *record = ((GlueXUserEventInformation*)info)->getOutputRecord();
285 if (record == 0) {
286 G4cerr(*G4cerr_p) << "GlueXSensitiveDetectorFMWPC::EndOfEvent error - "
287 << "hits seen but no output hddm record to save them into, "
288 << "cannot continue!" << G4endlstd::endl;
289 exit(1);
290 }
291
292 if (record->getPhysicsEvents().size() == 0)
293 record->addPhysicsEvents();
294 if (record->getHitViews().size() == 0)
295 record->getPhysicsEvent().addHitViews();
296 hddm_s::HitView &hitview = record->getPhysicsEvent().getHitView();
297 if (hitview.getForwardMWPCs().size() == 0)
298 hitview.addForwardMWPCs();
299 hddm_s::ForwardMWPC &fmwpc = hitview.getForwardMWPC();
300
301 // Collect and output the fmwpcTruthHits
302 for (siter = wires->begin(); siter != wires->end(); ++siter) {
303 std::vector<GlueXHitFMWPCwire::hitinfo_t> &hits = siter->second->hits;
304 // apply a pulse height threshold cut
305 for (unsigned int ih=0; ih < hits.size(); ++ih) {
306 if (hits[ih].dE_keV <= THRESH_KEV) {
307 hits.erase(hits.begin() + ih);
308 --ih;
309 }
310 }
311 if (hits.size() > 0) {
312 hddm_s::FmwpcChamberList wire = fmwpc.addFmwpcChambers(1);
313 wire(0).setLayer(siter->second->layer_);
314 wire(0).setWire(siter->second->wire_);
315 for (int ih=0; ih < (int)hits.size(); ++ih) {
316 hddm_s::FmwpcTruthHitList thit = wire(0).addFmwpcTruthHits(1);
317 thit(0).setDE(hits[ih].dE_keV);
318 thit(0).setDx(hits[ih].dx_cm);
319 thit(0).setT(hits[ih].t_ns);
320 }
321 }
322 }
323
324 // Collect and output the fmwpcTruthPoints
325 for (piter = points->begin(); piter != points->end(); ++piter) {
326 hddm_s::FmwpcTruthPointList point = fmwpc.addFmwpcTruthPoints(1);
327 point(0).setPrimary(piter->second->primary_);
328 point(0).setPtype(piter->second->ptype_G3);
329 point(0).setPx(piter->second->px_GeV);
330 point(0).setPy(piter->second->py_GeV);
331 point(0).setPz(piter->second->pz_GeV);
332 point(0).setE(piter->second->E_GeV);
333 point(0).setX(piter->second->x_cm);
334 point(0).setY(piter->second->y_cm);
335 point(0).setZ(piter->second->z_cm);
336 point(0).setT(piter->second->t_ns);
337 point(0).setTrack(piter->second->track_);
338 hddm_s::TrackIDList tid = point(0).addTrackIDs();
339 tid(0).setItrack(piter->second->trackID_);
340 }
341}
342
343int GlueXSensitiveDetectorFMWPC::GetIdent(std::string div,
344 const G4VTouchable *touch)
345{
346 const HddsG4Builder* bldr = GlueXDetectorConstruction::GetBuilder();
347 std::map<std::string, std::vector<int> >::const_iterator iter;
348 std::map<std::string, std::vector<int> > *identifiers;
349 int max_depth = touch->GetHistoryDepth();
350 for (int depth = 0; depth < max_depth; ++depth) {
351 G4VPhysicalVolume *pvol = touch->GetVolume(depth);
352 G4LogicalVolume *lvol = pvol->GetLogicalVolume();
353 int volId = fVolumeTable[lvol];
354 if (volId == 0) {
355 volId = bldr->getVolumeId(lvol);
356 fVolumeTable[lvol] = volId;
357 }
358 identifiers = &Refsys::fIdentifierTable[volId];
359 if ((iter = identifiers->find(div)) != identifiers->end()) {
360 int copyNum = touch->GetCopyNumber(depth);
361 copyNum += (dynamic_cast<G4PVPlacement*>(pvol))? -1 : 0;
362 return iter->second[copyNum];
363 }
364 }
365 return -1;
366}