Bug Summary

File:scratch/gluex/scan-build-work/hdgeant4/src/GlueXSensitiveDetectorCDC.cc
Location:line 756, column 20
Description:The left operand of '-' is a garbage value

Annotated Source Code

1//
2// GlueXSensitiveDetectorCDC - class implementation
3//
4// author: richard.t.jones at uconn.edu
5// version: may 28, 2015
6
7#include "GlueXSensitiveDetectorCDC.hh"
8#include "GlueXDetectorConstruction.hh"
9#include "GlueXPrimaryGeneratorAction.hh"
10#include "GlueXUserEventInformation.hh"
11#include "GlueXUserTrackInformation.hh"
12#include "GlueXUserOptions.hh"
13#include "HddmOutput.hh"
14
15#include <CLHEP/Random/RandPoisson.h>
16#include <Randomize.hh>
17
18#include "G4VPhysicalVolume.hh"
19#include "G4PVPlacement.hh"
20#include "G4EventManager.hh"
21#include "G4HCofThisEvent.hh"
22#include "G4Step.hh"
23#include "G4SDManager.hh"
24#include "G4ios.hh"
25
26#include <JANA/JApplication.h>
27
28#include <stdlib.h>
29
30const double fC = 1e-15 * coulomb;
31const double GlueXSensitiveDetectorCDC::ELECTRON_CHARGE = 1.6022e-4*fC;
32
33// Drift speed 2.2cm/us is appropriate for a 90/10 Argon/Methane mixture
34double GlueXSensitiveDetectorCDC::DRIFT_SPEED = 0.0055*cm/ns;
35
36// Minimum hit time difference for two hits on the same wire
37double GlueXSensitiveDetectorCDC::TWO_HIT_TIME_RESOL = 25*ns;
38
39// Cutoff on the total number of allowed straw hits
40int GlueXSensitiveDetectorCDC::MAX_HITS = 1000;
41
42// Minimum energy deposition for a straw hit (keV and mV)
43double GlueXSensitiveDetectorCDC::THRESH_KEV = 1.;
44double GlueXSensitiveDetectorCDC::THRESH_MV = 1.;
45
46// Drift distance can be somewhat larger than this in a field
47double GlueXSensitiveDetectorCDC::STRAW_RADIUS = 0.776*cm;
48
49// Straw hits are accepted from t=0 up to this maximum time
50double GlueXSensitiveDetectorCDC::CDC_TIME_WINDOW = 1000*ns;
51
52// Parameters for setting signal pulse height
53double GlueXSensitiveDetectorCDC::GAS_GAIN = 1e5;
54
55// Average number of secondary ion pairs for 50/50 Ar/CO2 mixture
56double GlueXSensitiveDetectorCDC::N_SECOND_PER_PRIMARY = 1.94;
57
58// Average energy needed to produce an ion pair for 50/50 mixture
59double GlueXSensitiveDetectorCDC::W_EFF_PER_ION = 29.5*eV;
60
61// These shared constants/tables initialized once from ccdb at run-time.
62// It is assumed that all threads share common values of all of these
63// lookup tables, and as long as all threads are working on events
64// belonging to the same run, this should be true.
65
66int GlueXSensitiveDetectorCDC::instanceCount = 0;
67int GlueXSensitiveDetectorCDC::fDrift_clusters = 0;
68double GlueXSensitiveDetectorCDC::fDrift_time[CDC_DRIFT_TABLE_LEN];
69double GlueXSensitiveDetectorCDC::fDrift_distance[CDC_DRIFT_TABLE_LEN];
70double GlueXSensitiveDetectorCDC::fBscale_par1;
71double GlueXSensitiveDetectorCDC::fBscale_par2;
72
73G4Mutex GlueXSensitiveDetectorCDC::fMutex = G4MUTEX_INITIALIZER{ { 0, 0, 0, 0, 0, 0, 0, { 0, 0 } } };
74
75GlueXSensitiveDetectorCDC::GlueXSensitiveDetectorCDC(const G4String& name)
76 : G4VSensitiveDetector(name),
77 fStrawsMap(0), fPointsMap(0)
78{
79 collectionName.insert("CDCStrawHitsCollection");
80 collectionName.insert("CDCPointsCollection");
81
82 // The rest of this only needs to happen once, the first time an object
83 // of this type is instantiated for this configuration of geometry and
84 // fields. If the geometry or fields change in such a way as to modify
85 // the drift-time properties of hits in the CDC, you must delete all old
86 // objects of this class and create new ones.
87
88 G4AutoLock barrier(&fMutex);
89 if (instanceCount++ == 0) {
90 int runno = HddmOutput::getRunNo();
91 extern jana::JApplication *japp;
92 if (japp == 0) {
93 G4cerr(*G4cerr_p) << "Error in GlueXSensitiveDetector constructor - "
94 << "jana global DApplication object not set, "
95 << "cannot continue." << G4endlstd::endl;
96 exit(-1);
97 }
98 jana::JCalibration *jcalib = japp->GetJCalibration(runno);
99 std::map<string, float> cdc_parms;
100 jcalib->Get("CDC/cdc_parms", cdc_parms);
101 DRIFT_SPEED = cdc_parms.at("CDC_DRIFT_SPEED")*cm/ns;
102 TWO_HIT_TIME_RESOL = cdc_parms.at("CDC_TWO_HIT_RESOL")*ns;
103 MAX_HITS = cdc_parms.at("CDC_MAX_HITS");
104 THRESH_KEV = cdc_parms.at("CDC_THRESH_KEV");
105
106 // Decide which drift time tables to load based on
107 // whether or not a magnetic field is present.
108
109 G4ThreeVector B = GlueXDetectorConstruction::GetInstance()
110 ->GetMagneticField(G4ThreeVector(0, 0, 65), tesla);
111 if (B.mag() > 1e-3) {
112 int nvalues = CDC_DRIFT_TABLE_LEN;
113 std::vector< std::map<std::string, float> > values;
114 jcalib->Get("CDC/cdc_drift_table", values);
115 for (int k=0; k < nvalues; ++k) {
116 fDrift_distance[k] = 0.01 * k; // 100 micron increments;
117 fDrift_time[k] = values[k]["t"] * 1000; // from us to ns
118 }
119
120 // The /CDC/drift_parameters consists of 2 rows and 10 columns
121 // The B-field scale parameters B1 and B2 are in columns 9 and 10
122 // They should be the same for both rows, so we get them from the first
123 std::vector< std::vector<float> > cdc_drift_parms;
124 jcalib->Get("/CDC/drift_parameters", cdc_drift_parms);
125 fBscale_par1 = cdc_drift_parms.at(0).at(9);
126 fBscale_par2 = cdc_drift_parms.at(0).at(10);
127 }
128 else {
129 int nvalues = CDC_DRIFT_TABLE_LEN;
130 std::vector< std::map<std::string, float> > values;
131 jcalib->Get("CDC/cdc_drift_table::NoBField", values);
132 for (int k=0; k < nvalues; ++k) {
133 fDrift_distance[k] = 0.01 * k; // 100 micron increments;
134 fDrift_time[k] = values[k]["t"] * 1000; // from us to ns
135 }
136 fBscale_par1 = 1.;
137 fBscale_par2 = 0;
138 }
139 G4cout(*G4cout_p) << "CDC: ALL parameters loaded from ccdb" << G4endlstd::endl;
140
141 // Check for "driftclusters" option in control.in
142
143 GlueXUserOptions *opts = GlueXUserOptions::GetInstance();
144 if (opts) {
145 std::map<int, int> driftclusters_opts;
146 if (opts->Find("driftclusters", driftclusters_opts))
147 fDrift_clusters = driftclusters_opts[1];
148 }
149 }
150}
151
152GlueXSensitiveDetectorCDC::GlueXSensitiveDetectorCDC(
153 const GlueXSensitiveDetectorCDC &src)
154 : G4VSensitiveDetector(src),
155 fStrawsMap(src.fStrawsMap), fPointsMap(src.fPointsMap)
156{
157 G4AutoLock barrier(&fMutex);
158 ++instanceCount;
159}
160
161GlueXSensitiveDetectorCDC &GlueXSensitiveDetectorCDC::operator=(const
162 GlueXSensitiveDetectorCDC &src)
163{
164 G4AutoLock barrier(&fMutex);
165 *(G4VSensitiveDetector*)this = src;
166 fStrawsMap = src.fStrawsMap;
167 fPointsMap = src.fPointsMap;
168 return *this;
169}
170
171GlueXSensitiveDetectorCDC::~GlueXSensitiveDetectorCDC()
172{
173 G4AutoLock barrier(&fMutex);
174 --instanceCount;
175}
176
177void GlueXSensitiveDetectorCDC::Initialize(G4HCofThisEvent* hce)
178{
179 fStrawsMap = new
180 GlueXHitsMapCDCstraw(SensitiveDetectorName, collectionName[0]);
181 fPointsMap = new
182 GlueXHitsMapCDCpoint(SensitiveDetectorName, collectionName[1]);
183 G4SDManager *sdm = G4SDManager::GetSDMpointer();
184 hce->AddHitsCollection(sdm->GetCollectionID(collectionName[0]), fStrawsMap);
185 hce->AddHitsCollection(sdm->GetCollectionID(collectionName[1]), fPointsMap);
186}
187
188G4bool GlueXSensitiveDetectorCDC::ProcessHits(G4Step* step,
189 G4TouchableHistory* ROhist)
190{
191 double dEsum = step->GetTotalEnergyDeposit();
192 if (dEsum == 0)
193 return false;
194
195 const G4ThreeVector &pin = step->GetPreStepPoint()->GetMomentum();
196 const G4ThreeVector &xin = step->GetPreStepPoint()->GetPosition();
197 const G4ThreeVector &xout = step->GetPostStepPoint()->GetPosition();
198 double tin = step->GetPreStepPoint()->GetGlobalTime();
199 double tout = step->GetPostStepPoint()->GetGlobalTime();
200 G4ThreeVector x = (xin + xout) / 2;
201 G4ThreeVector dx = xout - xin;
202 double t = (tin + tout) / 2;
203
204 const G4VTouchable* touch = step->GetPreStepPoint()->GetTouchable();
205 const G4AffineTransform &local_from_global = touch->GetHistory()
206 ->GetTopTransform();
207 G4ThreeVector xinlocal = local_from_global.TransformPoint(xin);
208 G4ThreeVector xoutlocal = local_from_global.TransformPoint(xout);
209
210 // For particles that range out inside the active volume, the
211 // "out" time may sometimes be set to something enormously high.
212 // This screws up the hit. Check for this case here by looking
213 // at tout and making sure it is less than 1 second. If it's
214 // not, then just use tin for "t".
215
216 if (tout > 1.0*s)
217 t = tin;
218
219 double drin = xinlocal.perp();
220 double drout = xoutlocal.perp();
221
222 // Find x2local = the point of closest approach to the z axis
223 G4ThreeVector xlocal = (xinlocal + xoutlocal) / 2;
224 G4ThreeVector dxlocal = xoutlocal - xinlocal;
225 double alpha = -(xinlocal[0] * dxlocal[0] +
226 xinlocal[1] * dxlocal[1]) / (dxlocal.perp2() + 1e-30);
227 G4ThreeVector xpoca = xinlocal + alpha * dxlocal;
228 double alpha01 = (alpha < 0)? 0 : (alpha > 1)? 1 : alpha;
229 G4ThreeVector x2local = xinlocal + alpha01 * dxlocal;
230
231 // Deal with tracks exiting the ends of the straws
232 if (fabs(x2local[2]) >= 75.45*cm) {
233 int sign = (xoutlocal[2] > 0)? 1 : -1;
234 int ring = GetIdent("ring", touch);
235 if (ring <= 4 || (ring >= 13 && ring <= 16) || ring >= 25) {
236 alpha = (sign * 75.45*cm - xinlocal[2]) / (dxlocal[2] + 1e-30);
237 xpoca = xinlocal + alpha * dxlocal;
238 alpha01 = (alpha < 0)? 0 : (alpha > 1)? 1 : alpha;
239 x2local = xinlocal + alpha01 * dxlocal;
240 }
241 else if (fabs(x2local[2]) >= 75.575*cm) {
242 alpha = (sign * 75.575*cm - xinlocal[2]) / (dxlocal[2] + 1e-30);
243 xpoca = xinlocal + alpha * dxlocal;
244 alpha01 = (alpha < 0)? 0 : (alpha > 1)? 1 : alpha;
245 x2local = xinlocal + alpha01 * dxlocal;
246 }
247 }
248
249 // Handle the case when the particle actually passes through
250 // the wire volume itself. For these cases, we should set the
251 // location of the hit to be the point on the wire itself. To
252 // determine if this is what is happening, we check drout to
253 // see if it is very close to the wire and drin to see if it is
254 // close to the tube. For the other case, when drin is close to
255 // the wire, we assume it is because it is emerging from the
256 // wire volume and automatically ignore those hits by returning
257 // immediately.
258
259 if (drin < 0.0050*cm) {
260 return false; /* entering straw within 50 microns of wire. ignore */
261 }
262
263 if ((drin > (STRAW_RADIUS - 0.0200*cm) && drout < 0.0050*cm) ||
264 (drin < 0.274*cm && drin > 0.234*cm && drout < 0.0050*cm))
265 {
266 // Either we entered within 200 microns of the straw tube and left
267 // within 50 microns of the wire or we entered the stub region near the
268 // donuts at either end of the straw (the inner radius of the feedthrough
269 // region is 0.254 cm) and passed near the wire. Assume the track passed
270 // through the wire volume.
271
272 x = xout;
273 t = tout;
274 xlocal = xoutlocal;
275
276 // For dx, we will just assume it is twice the distance from
277 // the straw to wire. Approximate the energy loss in the straw to
278 // be twice the energy loss in the first half of the straw.
279
280 dx *= 2;
281 dEsum *= 2;
282 }
283
284 double dradius = xpoca.perp();
285 double dr = dx.mag();
286 double dEdx = (dr > 1e-3*cm)? dEsum/dr : 0;
287
288 // Post the hit to the points list in the
289 // order of appearance in the event simulation.
290
291 int ring = GetIdent("ring", touch);
292 int sector = GetIdent("sector", touch);
293 G4Track *track = step->GetTrack();
294 G4int trackID = track->GetTrackID();
295 GlueXUserTrackInformation *trackinfo = (GlueXUserTrackInformation*)
296 track->GetUserInformation();
297 G4int itrack = trackinfo->GetGlueXTrackID();
298 int pdgtype = track->GetDynamicParticle()->GetPDGcode();
299 int g3type = GlueXPrimaryGeneratorAction::ConvertPdgToGeant3(pdgtype);
300 if (trackinfo->GetGlueXHistory() == 0 && itrack > 0) {
301 G4int key = fPointsMap->entries();
302 GlueXHitCDCpoint* lastPoint = (*fPointsMap)[key - 1];
303 // No more than one truth point per ring in the cdc
304 int ring = GetIdent("ring", touch);
305 if (lastPoint == 0 || lastPoint->track_ != trackID ||
306 lastPoint->ring_ != ring)
307 {
308 GlueXHitCDCpoint newPoint;
309 newPoint.ptype_G3 = g3type;
310 newPoint.track_ = trackID;
311 newPoint.trackID_ = itrack;
312 newPoint.primary_ = (track->GetParentID() == 0);
313 newPoint.t_ns = t/ns;
314 newPoint.z_cm = x[2]/cm;
315 newPoint.r_cm = x.perp()/cm;
316 newPoint.phi_rad = x.phi();
317 newPoint.dradius_cm = dradius/cm;
318 newPoint.px_GeV = pin[0]/GeV;
319 newPoint.py_GeV = pin[1]/GeV;
320 newPoint.pz_GeV = pin[2]/GeV;
321 newPoint.dEdx_GeV_cm = dEdx/(GeV/cm);
322 newPoint.sector_ = sector;
323 newPoint.ring_ = ring;
324 fPointsMap->add(key, newPoint);
325 }
326 }
327
328 // Post the hit to the straw hits map, ordered by straw index
329
330 if (dEsum > 0) {
331 int key = GlueXHitCDCstraw::GetKey(ring, sector);
332 GlueXHitCDCstraw *straw = (*fStrawsMap)[key];
333 if (straw == 0) {
334 GlueXHitCDCstraw newstraw(ring, sector);
335 fStrawsMap->add(key, newstraw);
336 straw = (*fStrawsMap)[key];
337 }
338
339 // Add the hit to the hits vector, maintaining track time ordering,
340 // re-ordering according to hit times will take place and end of event.
341
342 int merge_hit = 0;
343 hit_vector_t::iterator hiter;
344 for (hiter = straw->hits.begin(); hiter != straw->hits.end(); ++hiter) {
345 if (itrack == hiter->itrack_ && fabs(tin - hiter->t1_ns*ns) < 0.01) {
346 merge_hit = 1;
347 break;
348 }
349 else if (hiter->t0_ns*ns > tin) {
350 break;
351 }
352 }
353 if (merge_hit) {
354 double d_cm = x2local.perp()/cm;
355 if (d_cm < hiter->d_cm)
356 hiter->d_cm = d_cm;
357 hiter->q_fC += dEsum;
358 hiter->t1_ns = tout/ns;
359 hiter->x1_g = xout;
360 }
361 else if ((int)straw->hits.size() < MAX_HITS) { // create new hit
362 hiter = straw->hits.insert(hiter, GlueXHitCDCstraw::hitinfo_t());
363 hiter->track_ = trackID;
364 hiter->q_fC = dEsum;
365 hiter->d_cm = x2local.perp()/cm;
366 hiter->itrack_ = itrack;
367 hiter->ptype_G3 = g3type;
368 hiter->t0_ns = tin/ns;
369 hiter->t1_ns = tout/ns;
370 hiter->x0_g = xin;
371 hiter->x1_g = xout;
372 }
373 else {
374 G4cerr(*G4cerr_p) << "GlueXSensitiveDetectorCDC::ProcessHits error: "
375 << "max hit count " << MAX_HITS << " exceeded, truncating!"
376 << G4endlstd::endl;
377 }
378 }
379 return true;
380}
381
382void GlueXSensitiveDetectorCDC::EndOfEvent(G4HCofThisEvent*)
383{
384 std::map<int, GlueXHitCDCstraw*> *straws = fStrawsMap->GetMap();
385 std::map<int, GlueXHitCDCpoint*> *points = fPointsMap->GetMap();
386 if (straws->size() == 0 && points->size() == 0)
387 return;
388 std::map<int, GlueXHitCDCstraw*>::iterator siter;
389 std::map<int, GlueXHitCDCpoint*>::iterator piter;
390
391 if (verboseLevel > 1) {
392 G4cout(*G4cout_p) << G4endlstd::endl
393 << "--------> Hits Collection: in this event there are "
394 << straws->size() << " straws with hits in the CDC: "
395 << G4endlstd::endl;
396 for (siter = straws->begin(); siter != straws->end(); ++siter)
397 siter->second->Print();
398
399 G4cout(*G4cout_p) << G4endlstd::endl
400 << "--------> Hits Collection: in this event there are "
401 << points->size() << " truth points in the CDC: "
402 << G4endlstd::endl;
403 for (piter = points->begin(); piter != points->end(); ++piter)
404 piter->second->Print();
405 }
406
407 // pack hits into ouptut hddm record
408
409 G4EventManager* mgr = G4EventManager::GetEventManager();
410 G4VUserEventInformation* info = mgr->GetUserInformation();
411 hddm_s::HDDM *record = ((GlueXUserEventInformation*)info)->getOutputRecord();
412 if (record == 0) {
413 G4cerr(*G4cerr_p) << "GlueXSensitiveDetectorCDC::EndOfEvent error - "
414 << "hits seen but no output hddm record to save them into, "
415 << "cannot continue!" << G4endlstd::endl;
416 exit(1);
417 }
418
419 if (record->getPhysicsEvents().size() == 0)
420 record->addPhysicsEvents();
421 if (record->getHitViews().size() == 0)
422 record->getPhysicsEvent().addHitViews();
423 hddm_s::HitView &hitview = record->getPhysicsEvent().getHitView();
424 if (hitview.getCentralDCs().size() == 0)
425 hitview.addCentralDCs();
426 hddm_s::CentralDC &centralDC = hitview.getCentralDC();
427
428 // Collect and output the cdcTruthHits
429
430 for (siter = straws->begin(); siter != straws->end(); ++siter) {
431
432 // Merge multiple segments from a single track into one, and
433 // apply the drift time algorithm to get a single hit time for each.
434
435 hit_vector_t &splits = siter->second->hits;
436 hit_vector_t hits;
437 while (splits.size() > 0) {
438 for (unsigned int ih=1; ih < splits.size(); ++ih) {
439 if (fabs(splits[ih].t0_ns - splits[0].t1_ns) < 0.5*ns) {
440 splits[0].t1_ns = splits[ih].t1_ns;
441 splits[0].q_fC += splits[ih].q_fC;
442 if (splits[ih].d_cm < splits[0].d_cm)
443 splits[0].d_cm = splits[ih].d_cm;
444 splits[0].x1_g = splits[ih].x1_g;
445 splits.erase(splits.begin() + ih);
446 --ih;
447 }
448 }
449 if (splits[0].q_fC/keV > THRESH_KEV) {
450
451 // Simulate number of primary ion pairs.
452 // The total number of ion pairs depends on the energy deposition
453 // and the effective average energy to produce a pair.
454
455 // Average number of primary ion pairs
456 double t = (splits[0].t0_ns + splits[0].t1_ns)*ns / 2;
457 double n_p_mean = (splits[0].q_fC / W_EFF_PER_ION) /
458 (1 + N_SECOND_PER_PRIMARY);
459 // Number of generated primary ion pairs
460 int n_p = CLHEP::RandPoisson::shoot(n_p_mean);
461 if (fDrift_clusters == 0) {
462 add_cluster(hits, splits[0], n_p, t, splits[0].x0_g);
463 }
464 else {
465 // Loop over the number of primary ion pairs,
466 // generating a cluster at a random position
467 // along the track within the straw.
468 G4ThreeVector dx = splits[0].x1_g - splits[0].x0_g;
469 for (int n=0; n < n_p; n++) {
470 double u = G4RandFlatG4MTRandFlat::shoot();
471 G4ThreeVector x = splits[0].x0_g + u * dx;
472 add_cluster(hits, splits[0], 1, t, x);
473 }
474 }
475 }
476 splits.erase(splits.begin());
477 }
478
479 // If doing driftclusters generate a sampled waveform and
480 // analyze it to rebuild the hits list from scratch.
481
482 if (fDrift_clusters) {
483 // store waveform data in sampled sequence with 1 ns bins
484 int num_samples = int(CDC_TIME_WINDOW/ns);
485 double *samples = new double[num_samples];
486 for (int i=0; i < num_samples; i++) {
487 samples[i] = cdc_wire_signal_mV(double(i), hits);
488 }
489
490 // take the earliest hit to identify the track parameters
491 double dradius_cm = hits[0].d_cm;
492 int ptype_G3 = hits[0].ptype_G3;
493 int track = hits[0].track_;
494 int itrack = hits[0].itrack_;
495
496 hits.clear();
497 double q_mV_ns = 0.;
498 int over_threshold = 0;
499 for (int i=0; i < num_samples; ++i) {
500 if (samples[i] > THRESH_MV) {
501 if (!over_threshold) {
502 hits.push_back(GlueXHitCDCstraw::hitinfo_t());
503 hits.back().track_ = track;
504 hits.back().t_ns = double(i);
505 hits.back().d_cm = dradius_cm;
506 hits.back().itrack_ = itrack;
507 hits.back().ptype_G3 = ptype_G3;
508 over_threshold = 1;
509 }
510 q_mV_ns += samples[i];
511 }
512 else if (over_threshold) {
513 hits.back().q_fC = q_mV_ns; // warning -- faking the units
514 over_threshold = 0;
515 q_mV_ns = 0;
516 }
517 }
518 delete [] samples;
519 }
520 else {
521
522 // New reduced hit list is ordered by hit time
523
524 hit_vector_t::iterator hiter;
525 for (hiter = hits.begin(); hiter != hits.end(); ++hiter) {
526 int merge_hits = 0;
527 hit_vector_t::iterator iter;
528 for (iter = splits.begin(); iter != splits.end(); ++iter) {
529 if (fabs(hiter->t_ns - iter->t_ns) < TWO_HIT_TIME_RESOL) {
530 merge_hits = 1;
531 break;
532 }
533 else if (hiter->t_ns < iter->t_ns) {
534 break;
535 }
536 }
537 if (merge_hits) {
538 if (hiter->t_ns < iter->t_ns) {
539 double q_fC = iter->q_fC;
540 *iter = *hiter;
541 iter->q_fC += q_fC;
542 }
543 else {
544 iter->q_fC += hiter->q_fC;
545 }
546 }
547 else {
548 splits.insert(iter, *hiter);
549 }
550 }
551 }
552
553 if (hits.size() > 0) {
554 hddm_s::CdcStrawList straw = centralDC.addCdcStraws(1);
555 straw(0).setRing(siter->second->ring_);
556 straw(0).setStraw(siter->second->sector_);
557 for (int ih=0; ih < (int)hits.size(); ++ih) {
558 hddm_s::CdcStrawTruthHitList thit = straw(0).addCdcStrawTruthHits(1);
559 thit(0).setQ(hits[ih].q_fC);
560 thit(0).setT(hits[ih].t_ns);
561 thit(0).setD(hits[ih].d_cm);
562 thit(0).setItrack(hits[ih].itrack_);
563 thit(0).setPtype(hits[ih].ptype_G3);
564 }
565 }
566 }
567
568 // Collect and output the cdcTruthPoints
569 for (piter = points->begin(); piter != points->end(); ++piter) {
570 hddm_s::CdcTruthPointList point = centralDC.addCdcTruthPoints(1);
571 point(0).setDEdx(piter->second->dEdx_GeV_cm);
572 point(0).setDradius(piter->second->dradius_cm);
573 point(0).setPrimary(piter->second->primary_);
574 point(0).setPtype(piter->second->ptype_G3);
575 point(0).setPx(piter->second->px_GeV);
576 point(0).setPy(piter->second->py_GeV);
577 point(0).setPz(piter->second->pz_GeV);
578 point(0).setR(piter->second->r_cm);
579 point(0).setPhi(piter->second->phi_rad);
580 point(0).setZ(piter->second->z_cm);
581 point(0).setT(piter->second->t_ns);
582 point(0).setTrack(piter->second->track_);
583 hddm_s::TrackIDList tid = point(0).addTrackIDs();
584 tid(0).setItrack(piter->second->trackID_);
585 }
586}
587
588double GlueXSensitiveDetectorCDC::asic_response(double t_ns)
589{
590 // Simulation of the ASIC response to a pulse due to a cluster
591
592 double par[11] = {-0.01986, 0.01802, -0.001097, 10.3, 11.72,
593 -0.03701, 35.84, 15.93, 0.006141, 80.95, 24.77};
594 if (t_ns < par[3])
595 return par[0] * t_ns + par[1] * t_ns*t_ns + par[2] * t_ns*t_ns*t_ns;
596 else
597 return (par[0] * par[3] +
598 par[1] * par[3]*par[3] +
599 par[2] * par[3]*par[3]*par[3]) *
600 exp(-pow((t_ns - par[3])/par[4], 2)) +
601 par[5] * exp(-pow((t_ns - par[6])/par[7], 2)) +
602 par[8] * exp(-pow((t_ns - par[9])/par[10], 2));
603}
604
605double GlueXSensitiveDetectorCDC::cdc_wire_signal_mV(double t_ns,
606 hit_vector_t &hits)
607{
608 // Simulation of signal on a wire
609
610 double asic_gain = 0.5; // mV/fC
611 double signal_mV = 0;
612 hit_vector_t::iterator hiter;
613 for (hiter = hits.begin(); hiter != hits.end(); ++hiter) {
614 if (t_ns > hiter->t_ns) {
615 double my_time_ns = t_ns - hiter->t_ns;
616 signal_mV += asic_gain * hiter->q_fC * asic_response(my_time_ns);
617 }
618 }
619 return signal_mV;
620}
621
622void GlueXSensitiveDetectorCDC::add_cluster(hit_vector_t &hits,
623 GlueXHitCDCstraw::hitinfo_t &hit,
624 int n_p,
625 double t,
626 G4ThreeVector &x)
627{
628 // measured charge
629 double q_fC = 0;
630
631 // drift radius
632 double dradius_cm = hit.d_cm;
633
634 // Find the drift time for this cluster. Drift time depends on B:
635 // (dependence derived from Garfield calculations)
636
637 G4ThreeVector B = GlueXDetectorConstruction::GetInstance()
638 ->GetMagneticField(x, tesla);
639 double BmagT = B.mag();
640
641 // Check for closeness to boundaries of the drift table
642
643 double my_t_ns;
644 double my_t_err;
645 int i = (int)(dradius_cm * 100);
646 if (i >= CDC_DRIFT_TABLE_LEN - 3) {
647 // Do a crude linear extrapolation
648 my_t_ns = fDrift_time[CDC_DRIFT_TABLE_LEN - 3] +
649 (dradius_cm - fDrift_distance[CDC_DRIFT_TABLE_LEN - 3]) *
650 (fDrift_time[CDC_DRIFT_TABLE_LEN - 1] -
651 fDrift_time[CDC_DRIFT_TABLE_LEN - 3]) / 0.02;
652 }
653 else {
654 int index = (i < 1)? 0 : i - 1;
655
656 // Interpolate over the drift table to find
657 // an approximation for the drift time
658 polint(&fDrift_distance[index],
659 &fDrift_time[index], 4, dradius_cm, &my_t_ns, &my_t_err);
660 }
661 double tdrift_ns = my_t_ns / (fBscale_par1 + fBscale_par2 * BmagT);
662
663 // Prevent unphysical times (drift electrons arriving
664 // at wire before particle passes the doca to the wire)
665 double v_max = 0.08; // guess for now based on Garfield, near wire
666 double tmin_ns = dradius_cm / v_max;
667 if (tdrift_ns < tmin_ns) {
668 tdrift_ns = tmin_ns;
669 }
670 double total_time = t + tdrift_ns*ns;
671
672 // Skip cluster if the time would go beyond readout window
673 if (total_time > CDC_TIME_WINDOW) {
674 hit.q_fC = 0;
675 hit.t_ns = 0;
676 return;
677 }
678
679 if (fDrift_clusters == 0) {
680 // Total number of ion pairs. On average for each primary ion
681 // pair produced there are n_s secondary ion pairs produced. The
682 // probability distribution is a compound poisson distribution
683 // that requires generating two Poisson variables.
684 double n_s_mean = n_p * N_SECOND_PER_PRIMARY;
685 int n_s = CLHEP::RandPoisson::shoot(n_s_mean);
686 int n_t = n_s + n_p;
687 q_fC = n_t * GAS_GAIN * ELECTRON_CHARGE/fC;
688 }
689 else {
690 // Distribute the number of secondary ionizations for this primary
691 // ionization according to a Poisson distribution with mean n_s_over_p.
692 // For simplicity we assume these secondary electrons and the primary
693 // electron stay together as a cluster.
694 int n_s = CLHEP::RandPoisson::shoot(N_SECOND_PER_PRIMARY);
695 // Energy deposition, equivalent to anode charge, in units of fC
696 q_fC = GAS_GAIN * ELECTRON_CHARGE/fC * (1 + n_s);
697 }
698 hit.q_fC = q_fC;
699 hit.t_ns = total_time/ns;
700 hits.push_back(hit);
701}
702
703void GlueXSensitiveDetectorCDC::polint(double *xa, double *ya, int n,
704 double x, double *y, double *dy)
705{
706 // Slightly modified versions of the "polint" routine from
707 // Press, William H., Brian P. Flannery, Saul A. Teukolsky and
708 // William T. Vetterling, 1986, "Numerical Recipes: The Art of
709 // Scientific Computing" (Fortran), Cambrigde University Press,
710 // pp. 80-82.
711
712 double *c = NULL__null;
713 double *d = NULL__null;
714 double den;
715 double dif;
716 double dift;
717 double ho;
718 double hp;
719 double w;
720
721 int i;
722 int m;
723 int ns;
724
725 if ((c = (double*)malloc(n*sizeof(double))) == NULL__null ||
1
Taking false branch
→
726 (d = (double*)malloc(n*sizeof(double))) == NULL__null )
727 {
728 fprintf(stderrstderr, "polint error: allocating workspace\n" );
729 fprintf(stderrstderr, "polint error: setting y = 0 and dy = 1e9\n" );
730 *y = 0.0;
731 *dy = 1.e9;
732 if (c != NULL__null)
733 free(c);
734 if (d != NULL__null )
735 free(d);
736 return;
737 }
738
739 ns = 0;
740 dif = fabs(x-xa[0]);
741 for (i = 0; i < n; ++i) {
2
←
Assuming 'i' is >= 'n'
→
3
←
Loop condition is false. Execution continues on line 750
→
742 dift = fabs(x-xa[i]);
743 if (dift < dif) {
744 ns = i;
745 dif = dift;
746 }
747 c[i] = ya[i];
748 d[i] = ya[i];
749 }
750 *y = ya[ns];
751 ns = ns-1;
752 for (m = 0; m < n-1; ++m) {
4
←
Loop condition is true. Entering loop body
→
753 for (i = 0; i < n-m-1; ++i) {
5
←
Loop condition is true. Entering loop body
→
754 ho = xa[i]-x;
755 hp = xa[i+m+1]-x;
756 w = c[i+1]-d[i];
6
←
The left operand of '-' is a garbage value
757 den = ho-hp;
758 if (den == 0) {
759 fprintf( stderrstderr, "polint error: den = 0\n" );
760 fprintf( stderrstderr, "polint error: setting y = 0 and dy = 1e9\n" );
761 *y = 0.0;
762 *dy = 1.e9;
763 if (c != NULL__null)
764 free(c);
765 if (d != NULL__null)
766 free(d);
767 return;
768 }
769 den = w/den;
770 d[i] = hp*den;
771 c[i] = ho*den;
772 }
773 if (2*(ns+1) < n-m-1) {
774 *dy = c[ns+1];
775 }
776 else {
777 *dy = d[ns];
778 ns = ns-1;
779 }
780 *y = (*y)+(*dy);
781 }
782
783 if (c != NULL__null)
784 free(c);
785 if (d != NULL__null)
786 free(d);
787}
788
789int GlueXSensitiveDetectorCDC::GetIdent(std::string div,
790 const G4VTouchable *touch)
791{
792 const HddsG4Builder* bldr = GlueXDetectorConstruction::GetBuilder();
793 std::map<std::string, std::vector<int> >::const_iterator iter;
794 std::map<std::string, std::vector<int> > *identifiers;
795 int max_depth = touch->GetHistoryDepth();
796 for (int depth = 0; depth < max_depth; ++depth) {
797 G4VPhysicalVolume *pvol = touch->GetVolume(depth);
798 G4LogicalVolume *lvol = pvol->GetLogicalVolume();
799 int volId = fVolumeTable[lvol];
800 if (volId == 0) {
801 volId = bldr->getVolumeId(lvol);
802 fVolumeTable[lvol] = volId;
803 }
804 identifiers = &Refsys::fIdentifierTable[volId];
805 if ((iter = identifiers->find(div)) != identifiers->end()) {
806 int copyNum = touch->GetCopyNumber(depth);
807 copyNum += (dynamic_cast<G4PVPlacement*>(pvol))? -1 : 0;
808 return iter->second[copyNum];
809 }
810 }
811 return -1;
812}