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 number of allowed hits per chamber
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 {
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 }
274 return true;
275}
276
277void GlueXSensitiveDetectorFMWPC::EndOfEvent(G4HCofThisEvent*)
278{
279 std::map<int,GlueXHitFMWPCwire*> *wires = fWireHitsMap->GetMap();
280 std::map<int,GlueXHitFMWPCpoint*> *points = fPointsMap->GetMap();
281 if (wires->size() == 0 && points->size() == 0)
282 return;
283 std::map<int,GlueXHitFMWPCwire*>::iterator siter;
284 std::map<int,GlueXHitFMWPCpoint*>::iterator piter;
285
286 if (verboseLevel > 1) {
287 G4cout(*G4cout_p) << G4endlstd::endl
288 << "--------> Hits Collection: in this event there are "
289 << wires->size() << " wires with hits in the FMWPC: "
290 << G4endlstd::endl;
291 for (siter = wires->begin(); siter != wires->end(); ++siter)
292 siter->second->Print();
293
294 G4cout(*G4cout_p) << G4endlstd::endl
295 << "--------> Hits Collection: in this event there are "
296 << points->size() << " truth points in the FWMPC: "
297 << G4endlstd::endl;
298 for (piter = points->begin(); piter != points->end(); ++piter)
299 piter->second->Print();
300 }
301
302 // pack hits into ouptut hddm record
303
304 G4EventManager* mgr = G4EventManager::GetEventManager();
305 G4VUserEventInformation* info = mgr->GetUserInformation();
306 hddm_s::HDDM *record = ((GlueXUserEventInformation*)info)->getOutputRecord();
307 if (record == 0) {
308 G4cerr(*G4cerr_p) << "GlueXSensitiveDetectorFMWPC::EndOfEvent error - "
309 << "hits seen but no output hddm record to save them into, "
310 << "cannot continue!" << G4endlstd::endl;
311 exit(1);
312 }
313
314 if (record->getPhysicsEvents().size() == 0)
315 record->addPhysicsEvents();
316 if (record->getHitViews().size() == 0)
317 record->getPhysicsEvent().addHitViews();
318 hddm_s::HitView &hitview = record->getPhysicsEvent().getHitView();
319 if (hitview.getForwardMWPCs().size() == 0)
320 hitview.addForwardMWPCs();
321 hddm_s::ForwardMWPC &fmwpc = hitview.getForwardMWPC();
322
323 // Collect and output the fmwpcTruthHits
324 for (siter = wires->begin(); siter != wires->end(); ++siter) {
325 std::vector<GlueXHitFMWPCwire::hitinfo_t> &hits = siter->second->hits;
326 // apply a pulse height threshold cut
327 for (unsigned int ih=0; ih < hits.size(); ++ih) {
328 if (hits[ih].dE_keV <= THRESH_KEV) {
329 hits.erase(hits.begin() + ih);
330 --ih;
331 }
332 }
333 if (hits.size() > 0) {
334 hddm_s::FmwpcChamberList wire = fmwpc.addFmwpcChambers(1);
335 wire(0).setLayer(siter->second->layer_);
336 wire(0).setWire(siter->second->wire_);
337 int hitscount = hits.size();
338 if (hitscount > MAX_HITS) {
339 hitscount = MAX_HITS;
340 G4cerr(*G4cerr_p) << "GlueXSensitiveDetectorFMWPC::EndOfEvent warning: "
341 << "max hit count " << MAX_HITS << "exceeded, "
342 << hits.size() - hitscount << " hits discarded."
343 << G4endlstd::endl;
344 }
345 for (int ih=0; ih < (int)hits.size(); ++ih) {
346 double dE=hits[ih].dE_keV*keV;
347 double n_p_mean=dE/W_EFF_PER_ION/(1.+N_SECOND_PER_PRIMARY);
348 int n_p=CLHEP::RandPoisson::shoot(n_p_mean);
349 double n_s_mean=n_p*N_SECOND_PER_PRIMARY;
350 int n_s=CLHEP::RandPoisson::shoot(n_s_mean);
351 int n_t=n_p+n_s;
352 const double pC=1e-12*coulomb;
353 double q_pC=n_t*GAS_GAIN*ELECTRON_CHARGE/pC;
354
355 hddm_s::FmwpcTruthHitList thit = wire(0).addFmwpcTruthHits(1);
356 thit(0).setDE(dE);
357 thit(0).setT(hits[ih].t_ns);
358 thit(0).setDx(0.); // not used any more (SJT 2/28/22)
359 hddm_s::FmwpcTruthHitQList charges=thit(0).addFmwpcTruthHitQs(1);
360 charges(0).setD(hits[ih].d_cm);
361 charges(0).setQ(q_pC);
362 }
363 }
364 }
365
366 // Collect and output the fmwpcTruthPoints
367 for (piter = points->begin(); piter != points->end(); ++piter) {
368 hddm_s::FmwpcTruthPointList point = fmwpc.addFmwpcTruthPoints(1);
369 point(0).setPrimary(piter->second->primary_);
370 point(0).setPtype(piter->second->ptype_G3);
371 point(0).setPx(piter->second->px_GeV);
372 point(0).setPy(piter->second->py_GeV);
373 point(0).setPz(piter->second->pz_GeV);
374 point(0).setE(piter->second->E_GeV);
375 point(0).setX(piter->second->x_cm);
376 point(0).setY(piter->second->y_cm);
377 point(0).setZ(piter->second->z_cm);
378 point(0).setT(piter->second->t_ns);
379 point(0).setTrack(piter->second->track_);
380 hddm_s::TrackIDList tid = point(0).addTrackIDs();
381 tid(0).setItrack(piter->second->trackID_);
382 }
383}
384
385int GlueXSensitiveDetectorFMWPC::GetIdent(std::string div,
386 const G4VTouchable *touch)
387{
388 const HddsG4Builder* bldr = GlueXDetectorConstruction::GetBuilder();
389 std::map<std::string, std::vector<int> >::const_iterator iter;
390 std::map<std::string, std::vector<int> > *identifiers;
391 int max_depth = touch->GetHistoryDepth();
392 for (int depth = 0; depth < max_depth; ++depth) {
393 G4VPhysicalVolume *pvol = touch->GetVolume(depth);
394 G4LogicalVolume *lvol = pvol->GetLogicalVolume();
395 int volId = fVolumeTable[lvol];
396 if (volId == 0) {
397 volId = bldr->getVolumeId(lvol);
398 fVolumeTable[lvol] = volId;
399 }
400 identifiers = &Refsys::fIdentifierTable[volId];
401 if ((iter = identifiers->find(div)) != identifiers->end()) {
402 int copyNum = touch->GetCopyNumber(depth);
403 copyNum += (dynamic_cast<G4PVPlacement*>(pvol))? -1 : 0;
404 return iter->second[copyNum];
405 }
406 }
407 return -1;
408}