Bug Summary

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