Bug Summary

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

Annotated Source Code

1//
2// GlueXSensitiveDetectorFDC - class implementation
3//
4// author: richard.t.jones at uconn.edu
5// version: october 11, 2016
6
7#include "GlueXSensitiveDetectorFDC.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#include <math.h>
30
31const double fC = 1e-15 * coulomb;
32const double pC = 1e-12 * coulomb;
33const double GlueXSensitiveDetectorFDC::ELECTRON_CHARGE = 1.6022e-4*fC;
34
35// Drift speed 2.2cm/us is appropriate for a 90/10 Argon/Methane mixture
36double GlueXSensitiveDetectorFDC::DRIFT_SPEED = 0.0055*cm/ns;
37
38// Minimum hit time difference for two hits on the same wire
39double GlueXSensitiveDetectorFDC::TWO_HIT_TIME_RESOL = 25*ns;
40
41// Cutoff on the number of allowed hits per anode wire or cathode strip
42int GlueXSensitiveDetectorFDC::MAX_HITS = 1000;
43
44// Minimum energy deposition for a wire hit (keV, pC)
45double GlueXSensitiveDetectorFDC::THRESH_KEV = 1.;
46double GlueXSensitiveDetectorFDC::THRESH_ANODE = 1.*pC;
47double GlueXSensitiveDetectorFDC::THRESH_STRIPS = 5.*pC;
48
49// Straw hits are accepted from t=0 up to this maximum time
50double GlueXSensitiveDetectorFDC::FDC_TIME_WINDOW = 1000*ns;
51
52// Parameters for setting signal pulse height
53double GlueXSensitiveDetectorFDC::GAS_GAIN = 8e4;
54
55// Average number of secondary ion pairs for 40/60 Ar/CO2 mixture
56double GlueXSensitiveDetectorFDC::N_SECOND_PER_PRIMARY = 1.89;
57
58// Average energy needed to produce an ion pair for 40/60 mixture
59double GlueXSensitiveDetectorFDC::W_EFF_PER_ION = 30.2*eV;
60
61// Geometry parameters in the FDC chambers
62int GlueXSensitiveDetectorFDC::WIRES_PER_PLANE = 96;
63int GlueXSensitiveDetectorFDC::STRIPS_PER_PLANE = 192;
64double GlueXSensitiveDetectorFDC::ACTIVE_AREA_OUTER_RADIUS = 48.5*cm;
65double GlueXSensitiveDetectorFDC::ANODE_CATHODE_SPACING = 0.5*cm;
66double GlueXSensitiveDetectorFDC::WIRE_SPACING = 1.0*cm;
67double GlueXSensitiveDetectorFDC::STRIP_SPACING = 0.5*cm;
68double GlueXSensitiveDetectorFDC::U_OF_WIRE_ONE = -47.5*cm; //-(WIRES_PER_PLANE-1)*WIRE_SPACING/2
69double GlueXSensitiveDetectorFDC::U_OF_STRIP_ONE = -47.75*cm; //-(STRIPS_PER_PLANE-1)*STRIP_SPACING/2
70double GlueXSensitiveDetectorFDC::CATHODE_ROT_ANGLE = 1.309; // radians (75 degrees)
71double GlueXSensitiveDetectorFDC::STRIP_GAP = 0.1*cm;
72double GlueXSensitiveDetectorFDC::STRIP_NODES = 3;
73
74// Parameters for calculating the drift time-distance relation
75double GlueXSensitiveDetectorFDC::LORENTZ_NR_PAR1;
76double GlueXSensitiveDetectorFDC::LORENTZ_NR_PAR2;
77double GlueXSensitiveDetectorFDC::LORENTZ_NZ_PAR1;
78double GlueXSensitiveDetectorFDC::LORENTZ_NZ_PAR2;
79double GlueXSensitiveDetectorFDC::DRIFT_RES_PARMS[3];
80double GlueXSensitiveDetectorFDC::DRIFT_FUNC_PARMS[6];
81double GlueXSensitiveDetectorFDC::DRIFT_BSCALE_PAR1;
82double GlueXSensitiveDetectorFDC::DRIFT_BSCALE_PAR2;
83
84// Parameters for estimating magnetic field drift effects
85double GlueXSensitiveDetectorFDC::DIFFUSION_COEFF = 1.1e-6*cm*cm/s; // cm^2/s --> 200 microns at 1 cm
86double GlueXSensitiveDetectorFDC::K2 = 1.15;
87
88// Radius of deadened region around the beam
89double GlueXSensitiveDetectorFDC::wire_dead_zone_radius[4] =
90 {3.0*cm, 3.0*cm, 3.9*cm, 3.9*cm};
91double GlueXSensitiveDetectorFDC::strip_dead_zone_radius[4] =
92 {1.3*cm, 1.3*cm, 1.3*cm, 1.3*cm};
93
94// Drift time - distance lookup table
95int GlueXSensitiveDetectorFDC::drift_table_len;
96double *GlueXSensitiveDetectorFDC::drift_table_t_ns;
97double *GlueXSensitiveDetectorFDC::drift_table_d_cm;
98
99int GlueXSensitiveDetectorFDC::instanceCount = 0;
100G4Mutex GlueXSensitiveDetectorFDC::fMutex = G4MUTEX_INITIALIZER{ { 0, 0, 0, 0, 0, 0, 0, { 0, 0 } } };
101int GlueXSensitiveDetectorFDC::fDrift_clusters = 0;
102
103GlueXSensitiveDetectorFDC::GlueXSensitiveDetectorFDC(const G4String& name)
104 : G4VSensitiveDetector(name),
105 fWiresMap(0), fCathodesMap(0), fPointsMap(0)
106{
107 collectionName.insert("FDCWireHitsCollection");
108 collectionName.insert("FDCCathodeHitsCollection");
109 collectionName.insert("FDCPointsCollection");
110
111 // The rest of this only needs to happen once, the first time an object
112 // of this type is instantiated for this configuration of geometry and
113 // fields. If the geometry or fields change in such a way as to modify
114 // the drift-time properties of hits in the FDC, you must delete all old
115 // objects of this class and create new ones.
116
117 G4AutoLock barrier(&fMutex);
118 if (instanceCount++ == 0) {
119 int runno = HddmOutput::getRunNo();
120 extern jana::JApplication *japp;
121 if (japp == 0) {
122 G4cerr(*G4cerr_p) << "Error in GlueXSensitiveDetector constructor - "
123 << "jana global DApplication object not set, "
124 << "cannot continue." << G4endlstd::endl;
125 exit(-1);
126 }
127 jana::JCalibration *jcalib = japp->GetJCalibration(runno);
128 std::map<string, float> fdc_parms;
129 jcalib->Get("FDC/fdc_parms", fdc_parms);
130 DRIFT_SPEED = fdc_parms.at("FDC_DRIFT_SPEED")*cm/ns;
131 ACTIVE_AREA_OUTER_RADIUS = fdc_parms.at("FDC_ACTIVE_AREA_OUTER_RADIUS")*cm;
132 ANODE_CATHODE_SPACING = fdc_parms.at("FDC_ANODE_CATHODE_SPACING")*cm;
133 TWO_HIT_TIME_RESOL = fdc_parms.at("FDC_TWO_HIT_RESOL")*ns;
134 WIRES_PER_PLANE = fdc_parms.at("FDC_WIRES_PER_PLANE");
135 WIRE_SPACING = fdc_parms.at("FDC_WIRE_SPACING")*cm;
136 STRIPS_PER_PLANE = fdc_parms.at("FDC_STRIPS_PER_PLANE");
137 STRIP_SPACING = fdc_parms.at("FDC_STRIP_SPACING")*cm;
138 STRIP_GAP = fdc_parms.at("FDC_STRIP_GAP")*cm;
139 MAX_HITS = fdc_parms.at("FDC_MAX_HITS");
140 K2 = fdc_parms.at("FDC_K2");
141 STRIP_NODES = fdc_parms.at("FDC_STRIP_NODES");
142 THRESH_KEV = fdc_parms.at("FDC_THRESH_KEV");
143 THRESH_STRIPS = fdc_parms.at("FDC_THRESH_STRIPS");
144 DIFFUSION_COEFF = fdc_parms.at("FDC_DIFFUSION_COEFF")*cm*cm/s;
145 U_OF_WIRE_ONE = -(WIRES_PER_PLANE -1) * WIRE_SPACING / 2;
146 U_OF_STRIP_ONE = -(STRIPS_PER_PLANE -1) * STRIP_SPACING / 2;
147
148 // Parameters for correcting for deflection due to Lorentz force
149 std::map<string, double> lorentz_parms;
150 jcalib->Get("FDC/lorentz_deflection_parms", lorentz_parms);
151 LORENTZ_NR_PAR1 = lorentz_parms["nr_par1"];
152 LORENTZ_NR_PAR2 = lorentz_parms["nr_par2"];
153 LORENTZ_NZ_PAR1 = lorentz_parms["nz_par1"];
154 LORENTZ_NZ_PAR2 = lorentz_parms["nz_par2"];
155
156 // Parameters for accounting for variation in drift distance from FDC
157 std::map<string, double> drift_res_parms;
158 jcalib->Get("FDC/drift_resolution_parms", drift_res_parms);
159 DRIFT_RES_PARMS[0] = drift_res_parms["p0"];
160 DRIFT_RES_PARMS[1] = drift_res_parms["p1"];
161 DRIFT_RES_PARMS[2] = drift_res_parms["p2"];
162
163 // Time-to-distance function parameters for FDC
164 std::map<string, double> drift_func_parms;
165 jcalib->Get("FDC/drift_function_parms", drift_func_parms);
166 DRIFT_FUNC_PARMS[0] = drift_func_parms["p0"];
167 DRIFT_FUNC_PARMS[1] = drift_func_parms["p1"];
168 DRIFT_FUNC_PARMS[2] = drift_func_parms["p2"];
169 DRIFT_FUNC_PARMS[3] = drift_func_parms["p3"];
170 DRIFT_FUNC_PARMS[4] = 1000.;
171 DRIFT_FUNC_PARMS[5] = 0.;
172 std::map<string, double> drift_func_ext;
173 if (jcalib->Get("FDC/drift_function_ext", drift_func_ext) == false) {
174 DRIFT_FUNC_PARMS[4] = drift_func_ext["p4"];
175 DRIFT_FUNC_PARMS[5] = drift_func_ext["p5"];
176 }
177
178 // Factors for taking care of B-dependence of drift time for FDC
179 std::map<string, double> fdc_drift_parms;
180 jcalib->Get("FDC/fdc_drift_parms", fdc_drift_parms);
181 DRIFT_BSCALE_PAR1 = fdc_drift_parms["bscale_par1"];
182 DRIFT_BSCALE_PAR2 = fdc_drift_parms["bscale_par2"];
183
184 // Build a lookup table of drift time->distance for the FDC,
185 // used in the code to build an efficient reverse-map function.
186 drift_table_len = 1000;
187 drift_table_t_ns = new double[drift_table_len];
188 drift_table_d_cm = new double[drift_table_len];
189 double thigh = DRIFT_FUNC_PARMS[4];
190 double tstep = 0.5; //ns
191 for (int j=0; j < drift_table_len; j++) {
192 double t = j * tstep;
193 if (t < thigh) {
194 double t2 = t*t;
195 drift_table_t_ns[j] = t;
196 drift_table_d_cm[j] = DRIFT_FUNC_PARMS[0] * sqrt(t) +
197 DRIFT_FUNC_PARMS[1] * t +
198 DRIFT_FUNC_PARMS[2] * t2 +
199 DRIFT_FUNC_PARMS[3] * t*t2;
200 }
201 else {
202 double thigh2 = thigh * thigh;
203 drift_table_t_ns[j] = t;
204 drift_table_d_cm[j] = DRIFT_FUNC_PARMS[0] * sqrt(thigh) +
205 DRIFT_FUNC_PARMS[1] * thigh +
206 DRIFT_FUNC_PARMS[2] * thigh2 +
207 DRIFT_FUNC_PARMS[3] * thigh2 * thigh +
208 DRIFT_FUNC_PARMS[5] * (t - thigh);
209 }
210 }
211
212 G4cout(*G4cout_p) << "FDC: ALL parameters loaded from ccdb" << G4endlstd::endl;
213
214 // Check for "driftclusters" option in control.in
215
216 GlueXUserOptions *opts = GlueXUserOptions::GetInstance();
217 if (opts) {
218 std::map<int, int> driftclusters_opts;
219 if (opts->Find("driftclusters", driftclusters_opts))
220 fDrift_clusters = driftclusters_opts[1];
221 }
222 }
223}
224
225GlueXSensitiveDetectorFDC::GlueXSensitiveDetectorFDC(
226 const GlueXSensitiveDetectorFDC &src)
227 : G4VSensitiveDetector(src),
228 fWiresMap(src.fWiresMap),
229 fCathodesMap(src.fCathodesMap),
230 fPointsMap(src.fPointsMap)
231{
232 G4AutoLock barrier(&fMutex);
233 ++instanceCount;
234}
235
236GlueXSensitiveDetectorFDC &GlueXSensitiveDetectorFDC::operator=(const
237 GlueXSensitiveDetectorFDC &src)
238{
239 G4AutoLock barrier(&fMutex);
240 *(G4VSensitiveDetector*)this = src;
241 fWiresMap = src.fWiresMap;
242 fCathodesMap = src.fCathodesMap;
243 fPointsMap = src.fPointsMap;
244 return *this;
245}
246
247GlueXSensitiveDetectorFDC::~GlueXSensitiveDetectorFDC()
248{
249 G4AutoLock barrier(&fMutex);
250 --instanceCount;
251}
252
253void GlueXSensitiveDetectorFDC::Initialize(G4HCofThisEvent* hce)
254{
255 fWiresMap = new
256 GlueXHitsMapFDCwire(SensitiveDetectorName, collectionName[0]);
257 fCathodesMap = new
258 GlueXHitsMapFDCcathode(SensitiveDetectorName, collectionName[1]);
259 fPointsMap = new
260 GlueXHitsMapFDCpoint(SensitiveDetectorName, collectionName[2]);
261 G4SDManager *sdm = G4SDManager::GetSDMpointer();
262 hce->AddHitsCollection(sdm->GetCollectionID(collectionName[0]), fWiresMap);
263 hce->AddHitsCollection(sdm->GetCollectionID(collectionName[1]), fCathodesMap);
264 hce->AddHitsCollection(sdm->GetCollectionID(collectionName[2]), fPointsMap);
265}
266
267G4bool GlueXSensitiveDetectorFDC::ProcessHits(G4Step* step,
268 G4TouchableHistory* ROhist)
269{
270 double dEsum = step->GetTotalEnergyDeposit();
271 if (dEsum == 0)
272 return false;
273
274 G4ThreeVector pin = step->GetPreStepPoint()->GetMomentum();
275 G4ThreeVector xin = step->GetPreStepPoint()->GetPosition();
276 G4ThreeVector xout = step->GetPostStepPoint()->GetPosition();
277 double Ein = step->GetPreStepPoint()->GetTotalEnergy();
278 double tin = step->GetPreStepPoint()->GetGlobalTime();
279 double tout = step->GetPostStepPoint()->GetGlobalTime();
280
281 G4Track *track = step->GetTrack();
282 int trackID = track->GetTrackID();
283 const G4VTouchable* touch = step->GetPreStepPoint()->GetTouchable();
284 int package = GetIdent("package", touch);
285 int layer = GetIdent("layer", touch);
286 if (layer == 0) {
287 fprintf(stderrstderr, "hitFDC error: FDC layer number evaluates to zero! "
288 "THIS SHOULD NEVER HAPPEN! drop this particle.\n");
289 return false;
290 }
291 else if (package == 0) {
292 fprintf(stderrstderr, "hitFDC error: FDC package number evaluates to zero! "
293 "THIS SHOULD NEVER HAPPEN! drop this particle.\n");
294 return false;
295 }
296 // Normally numeric identifiers start at 1, eg. layer, package, module
297 // but if it is an index counting from zero, add the "No" suffix.
298 int packNo = package - 1;
299 int module = 2 * packNo + ((layer - 1) / 3) + 1;
300 int chamber = (module * 10) + ((layer - 1) % 3) + 1;
301
302#ifdef MERGE_STEPS_BEFORE_HITS_GENERATION
303
304 // This section forces hdgeant4 to conform to the behavior of hdgeant
305 // in merging together all of the steps taken by a single track as it
306 // moves through a single fdc chamber (plane) before generating hits.
307 // For this purpose, it borrows existing class GlueXHitFDCpoint to
308 // save the information from previous tracking steps by this track
309 // inside the current fdc active volume, and merging the saved info
310 // into the current step when the track either stops inside the fdc
311 // active volume or when it exits through an exterior boundary.
312 //
313 // NOTE
314 // This algorithm has a blind spot for tracks that stop inside the
315 // wire layer because the track disappears without ever crossing
316 // through an exterior boundary. This condition produces a warning
317 // message if verboseLevel > 0, but it is an inherent inefficiency
318 // in this model, and should not be considered to be a bug.
319
320 G4int rkey = 1 << 24; // reserved key for this thread
321 std::map<int,GlueXHitFDCpoint*>::iterator saved_segment =
322 fPointsMap->GetMap()->find(rkey);
323 G4String invol = step->GetPostStepPoint()->GetTouchable()
324 ->GetVolume()->GetName();
325 if (track->GetTrackStatus() == fAlive &&
326 (invol.index("FDA") == 0 || invol.index("FDX") == 0))
327 {
328 if (saved_segment == fPointsMap->GetMap()->end()) {
329 GlueXHitFDCpoint segment(chamber);
330 segment.track_ = trackID;
331 segment.x_cm = xin[0]; // this is borrowed storage for
332 segment.y_cm = xin[1]; // local track step information,
333 segment.z_cm = xin[2]; // don't worry about the units!
334 segment.t_ns = tin;
335 segment.px_GeV = pin[0];
336 segment.py_GeV = pin[1];
337 segment.pz_GeV = pin[2];
338 segment.E_GeV = Ein;
339 segment.dEdx_GeV_cm = dEsum;
340 fPointsMap->add(rkey, segment);
341 return true;
342 }
343 else if (saved_segment->second->track_ == trackID) {
344 saved_segment->second->dEdx_GeV_cm += dEsum;
345 return true;
346 }
347 else {
348 if (verboseLevel > 0)
349 fprintf(stderrstderr, "hitFDC warning: FDC saved track segment "
350 "was lost, drop this step.\n");
351 fPointsMap->GetMap()->erase(saved_segment);
352 return ProcessHits(step, ROhist);
353 }
354 }
355 else if (saved_segment != fPointsMap->GetMap()->end()) {
356 if (saved_segment->second->track_ == trackID) {
357 xin.setX(saved_segment->second->x_cm);
358 xin.setY(saved_segment->second->y_cm);
359 xin.setZ(saved_segment->second->z_cm);
360 tin = saved_segment->second->t_ns;
361 pin.setX(saved_segment->second->px_GeV);
362 pin.setY(saved_segment->second->py_GeV);
363 pin.setZ(saved_segment->second->pz_GeV);
364 Ein = saved_segment->second->E_GeV;
365 dEsum += saved_segment->second->dEdx_GeV_cm;
366 fPointsMap->GetMap()->erase(saved_segment);
367 }
368 else {
369 if (verboseLevel > 0)
370 fprintf(stderrstderr, "hitFDC warning: FDC saved track segment "
371 "was orphaned, drop this step.\n");
372 fPointsMap->GetMap()->erase(saved_segment);
373 return ProcessHits(step, ROhist);
374 }
375 }
376
377#endif
378
379 G4ThreeVector x = (xin + xout) / 2;
380 G4ThreeVector dx = xout - xin;
381 double t = (tin + tout) / 2;
382 double dr = dx.mag();
383 double dEdx = (dr > 1e-3*cm)? dEsum/dr : 0;
384
385 const G4AffineTransform &local_from_global = touch->GetHistory()
386 ->GetTopTransform();
387 G4ThreeVector xinlocal = local_from_global.TransformPoint(xin);
388 G4ThreeVector xoutlocal = local_from_global.TransformPoint(xout);
389
390 // For particles that range out inside the active volume, the
391 // "out" time may sometimes be set to something enormously high.
392 // This screws up the hit. Check for this case here by looking
393 // at tout and making sure it is less than 1 second. If it's
394 // not, then just use tin for "t".
395
396 if (tout > 1.0*s)
397 t = tin;
398
399 double alpha = atan2(xoutlocal[0] - xinlocal[0], xoutlocal[2] - xinlocal[2]);
400 double sinalpha = sin(alpha);
401 double cosalpha = cos(alpha);
402 G4ThreeVector xlocal = (xinlocal + xoutlocal) / 2;
403
404 // Make a fuzzy boundary around the forward dead region
405 // by killing any track segment whose midpoint is within the boundary
406
407 if (xlocal.perp() < wire_dead_zone_radius[packNo])
408 return false;
409
410 int wire = ceil((xlocal[0] - U_OF_WIRE_ONE) / WIRE_SPACING + 0.5);
411 double xwire = U_OF_WIRE_ONE + (wire - 1) * WIRE_SPACING;
412 double uwire = xinlocal[2];
413 double vwire = xinlocal[0] - xwire;
414 double dradius = fabs(vwire * cosalpha - uwire * sinalpha);
415 int pdgtype = track->GetDynamicParticle()->GetPDGcode();
416 int g3type = GlueXPrimaryGeneratorAction::ConvertPdgToGeant3(pdgtype);
417
418#if VERBOSE_PRINT_POINTS
419 std::cout << "fdc in u,v=" << uwire << "," << vwire
420 << " out u,v=" << xoutlocal[2] << "," << xoutlocal[0] - xwire
421 << " dradius=" << dradius
422 << " in=" << step->GetPreStepPoint()->GetTouchable()->GetVolume()->GetName()
423 << " out=" << step->GetPostStepPoint()->GetTouchable()->GetVolume()->GetName()
424 << std::endl;
425#endif
426
427 // Post the hit to the points list in the
428 // order of appearance in the event simulation.
429
430 GlueXUserTrackInformation *trackinfo = (GlueXUserTrackInformation*)
431 track->GetUserInformation();
432 int itrack = trackinfo->GetGlueXTrackID();
433 if (trackinfo->GetGlueXHistory() == 0 && itrack > 0) {
434 G4int key = (chamber << 20) + fPointsMap->entries();
435 GlueXHitFDCpoint* lastPoint = (*fPointsMap)[key - 1];
436 // Limit fdc truthPoints to one per chamber
437 if (lastPoint == 0 || lastPoint->track_ != trackID ||
438 lastPoint->chamber_ != chamber)
439 {
440 GlueXHitFDCpoint newPoint(chamber);
441 newPoint.primary_ = (track->GetParentID() == 0);
442 newPoint.track_ = trackID;
443 newPoint.x_cm = xout[0]/cm;
444 newPoint.y_cm = xout[1]/cm;
445 newPoint.z_cm = xout[2]/cm;
446 newPoint.t_ns = tout/ns;
447 newPoint.px_GeV = pin[0]/GeV;
448 newPoint.py_GeV = pin[1]/GeV;
449 newPoint.pz_GeV = pin[2]/GeV;
450 newPoint.E_GeV = Ein/GeV;
451 newPoint.dradius_cm = dradius/cm;
452 newPoint.dEdx_GeV_cm = dEdx/(GeV/cm);
453 newPoint.ptype_G3 = g3type;
454 newPoint.trackID_ = itrack;
455 fPointsMap->add(key, newPoint);
456 }
457 }
458
459 // Post the hit to the hits tree, ordered by wire, strip number
460
461 if (dEsum > 0) {
462 double u0 = xinlocal[0];
463 double u1 = xoutlocal[0];
464 int wire1 = ceil((u0 - U_OF_WIRE_ONE) / WIRE_SPACING + 0.5);
465 int wire2 = ceil((u1 - U_OF_WIRE_ONE) / WIRE_SPACING + 0.5);
466
467 // Check that wire numbers are not out of range,
468 // making sure at least one wire number is valid
469 if (wire1 > WIRES_PER_PLANE && wire2 > WIRES_PER_PLANE)
470 return false;
471 else if (wire1 < 1 && wire2 < 1)
472 return false;
473 wire1 = (wire1 > WIRES_PER_PLANE)? WIRES_PER_PLANE :
474 (wire1 < 1)? 1 : wire1;
475 wire2 = (wire2 > WIRES_PER_PLANE)? WIRES_PER_PLANE :
476 (wire2 < 1)? 1 : wire2;
477 int dwire = (wire1 < wire2)? 1 : -1;
478
479 // deal with the case of tracks crossing two cells
480 for (int wire = wire1; wire != wire2 + dwire; wire += dwire) {
481 double xwire = U_OF_WIRE_ONE + (wire - 1) * WIRE_SPACING;
482 G4ThreeVector x0;
483 G4ThreeVector x1;
484 double dE;
485 if (wire1 == wire2) {
486 dE = dEsum;
487 x0 = xinlocal;
488 x1 = xoutlocal;
489 }
490 else {
491 x0[0] = xwire - 0.5 * dwire * WIRE_SPACING;
492 x0[1] = xinlocal[1] + (x0[0] - xinlocal[0] + 1e-20) *
493 (xoutlocal[1] - xinlocal[1]) /
494 (xoutlocal[0] - xinlocal[0] + 1e-20);
495 x0[2] = xinlocal[2] + (x0[0] - xinlocal[0] + 1e-20) *
496 (xoutlocal[2] - xinlocal[2]) /
497 (xoutlocal[0] - xinlocal[0] + 1e-20);
498 if (fabs(x0[2] - xoutlocal[2]) > fabs(xinlocal[2] - xoutlocal[2]))
499 x0 = xinlocal;
500
501 x1[0] = xwire + 0.5 * dwire * WIRE_SPACING;
502 x1[1] = xinlocal[1] + (x1[0] - xinlocal[0] + 1e-20) *
503 (xoutlocal[1] - xinlocal[1]) /
504 (xoutlocal[0] - xinlocal[0] + 1e-20);
505 x1[2] = xinlocal[2] + (x1[0] - xinlocal[0] + 1e-20) *
506 (xoutlocal[2] - xinlocal[2]) /
507 (xoutlocal[0] - xinlocal[0] + 1e-20);
508 if (fabs(x1[2] - xinlocal[2]) > fabs(xoutlocal[2] - xinlocal[2]))
509 x1 = xoutlocal;
510
511 dE = dEsum * (x1[2] - x0[2]) /
512 (xoutlocal[2] - xinlocal[2] + 1e-20);
513 }
514
515 int key = GlueXHitFDCwire::GetKey(chamber, wire);
516 GlueXHitFDCwire *anode = (*fWiresMap)[key];
517 if (anode == 0) {
518 GlueXHitFDCwire newanode(chamber, wire);
519 fWiresMap->add(key, newanode);
520 anode = (*fWiresMap)[key];
521 }
522
523 // Add the hit to the hits vector, maintaining track time ordering,
524 // re-ordering according to hit times will take place at end of event.
525
526 int merge_hits = 0;
527 std::vector<GlueXHitFDCwire::hitinfo_t>::iterator hiter;
528 for (hiter = anode->hits.begin(); hiter != anode->hits.end(); ++hiter) {
529 if (itrack == hiter->itrack_ && fabs(tin - hiter->t1_ns*ns) < 0.01) {
530 merge_hits = 1;
531 break;
532 }
533 else if (hiter->t0_ns*ns > tin) {
534 break;
535 }
536 }
537 if (merge_hits) {
538 hiter->dE_keV += dE/keV;
539 hiter->t1_ns = tout/ns;
540 hiter->x1_g = xout;
541 hiter->x1_l = x1;
542 }
543 else {
544 // create new hit
545 hiter = anode->hits.insert(hiter, GlueXHitFDCwire::hitinfo_t());
546 hiter->dE_keV = dE/keV;
547 hiter->itrack_ = itrack;
548 hiter->ptype_G3 = g3type;
549 hiter->t0_ns = tin/ns;
550 hiter->t1_ns = tout/ns;
551 hiter->x0_g = xin;
552 hiter->x1_g = xout;
553 hiter->x0_l = x0;
554 hiter->x1_l = x1;
555 }
556 }
557 }
558 return true;
559}
560
561void GlueXSensitiveDetectorFDC::EndOfEvent(G4HCofThisEvent*)
562{
563 std::map<int,GlueXHitFDCwire*> *wires = fWiresMap->GetMap();
564 std::map<int,GlueXHitFDCcathode*> *strips = fCathodesMap->GetMap();
565 std::map<int,GlueXHitFDCpoint*> *points = fPointsMap->GetMap();
566 if (wires->size() == 0 && strips->size() == 0 && points->size() == 0)
567 return;
568 std::map<int,GlueXHitFDCwire*>::iterator witer;
569 std::map<int,GlueXHitFDCcathode*>::iterator siter;
570 std::map<int,GlueXHitFDCpoint*>::iterator piter;
571
572 if (verboseLevel > 1) {
573 G4cout(*G4cout_p) << G4endlstd::endl
574 << "--------> Hits Collection: in this event there are "
575 << wires->size() << " anode wires with hits in the FDC: "
576 << G4endlstd::endl;
577 for (witer = wires->begin(); witer != wires->end(); ++witer)
578 witer->second->Print();
579
580 G4cout(*G4cout_p) << G4endlstd::endl
581 << "--------> Hits Collection: in this event there are "
582 << strips->size() << " cathode strips with hits in the FDC: "
583 << G4endlstd::endl;
584 for (siter = strips->begin(); siter != strips->end(); ++siter)
585 siter->second->Print();
586
587 G4cout(*G4cout_p) << G4endlstd::endl
588 << "--------> Hits Collection: in this event there are "
589 << points->size() << " truth points in the FDC: "
590 << G4endlstd::endl;
591 for (piter = points->begin(); piter != points->end(); ++piter)
592 piter->second->Print();
593 }
594
595 // pack hits into ouptut hddm record
596
597 G4EventManager* mgr = G4EventManager::GetEventManager();
598 G4VUserEventInformation* info = mgr->GetUserInformation();
599 hddm_s::HDDM *record = ((GlueXUserEventInformation*)info)->getOutputRecord();
600 if (record == 0) {
601 G4cerr(*G4cerr_p) << "GlueXSensitiveDetectorFDC::EndOfEvent error - "
602 << "hits seen but no output hddm record to save them into, "
603 << "cannot continue!" << G4endlstd::endl;
604 exit(1);
605 }
606
607 if (record->getPhysicsEvents().size() == 0)
608 record->addPhysicsEvents();
609 if (record->getHitViews().size() == 0)
610 record->getPhysicsEvent().addHitViews();
611 hddm_s::HitView &hitview = record->getPhysicsEvent().getHitView();
612 if (hitview.getForwardDCs().size() == 0)
613 hitview.addForwardDCs();
614 hddm_s::ForwardDC &forwardDC = hitview.getForwardDC();
615
616 // Collect and output the wireTruthHits
617
618 for (witer = wires->begin(); witer != wires->end(); ++witer) {
619
620 // Merge multiple segments from a single track into one, and
621 // apply the drift time algorithm to get a single hit time for each.
622
623 int chamber = witer->second->chamber_;
624 int wire = witer->second->wire_;
625 int module = chamber / 10;
626 int packNo = (module - 1) / 2;
627 int layer = chamber % 10;
628 int glayerNo = 3*(layer - 1) + module-1;
629 int global_wire_number = 96 * glayerNo + wire - 1;
630 std::vector<GlueXHitFDCwire::hitinfo_t> &splits = witer->second->hits;
631 std::vector<GlueXHitFDCwire::hitinfo_t> hits;
632 while (splits.size() > 0) {
633 for (unsigned int ih=1; ih < splits.size(); ++ih) {
634 if (fabs(splits[ih].t0_ns - splits[0].t1_ns) < 0.5*ns) {
635 splits[0].dE_keV += splits[ih].dE_keV;
636 splits[0].t1_ns = splits[ih].t1_ns;
637 splits[0].x1_g = splits[ih].x1_g;
638 splits[0].x1_l = splits[ih].x1_l;
639 splits.erase(splits.begin() + ih);
640 --ih;
641 }
642 }
643
644 // Simulate the number of primary ion pairs.
645 // The total number of ion pairs depends on the energy deposition
646 // and the effective average energy to produce a pair, w_eff.
647 // On average for each primary ion pair produced there are n_s_per_p
648 // secondary ion pairs produced.
649
650 double xwire = U_OF_WIRE_ONE + (wire - 1) * WIRE_SPACING;
651 double dE = splits[0].dE_keV*keV;
652 if (dE > THRESH_KEV*keV) {
653 // Average number of primary ion pairs
654 double n_p_mean = dE / W_EFF_PER_ION / (1 + N_SECOND_PER_PRIMARY);
655 // number of primary ion pairs
656 int n_p = CLHEP::RandPoisson::shoot(n_p_mean);
657 G4ThreeVector &x0 = splits[0].x0_l;
658 G4ThreeVector &x1 = splits[0].x1_l;
659
660 if (fDrift_clusters == 0) {
661 G4ThreeVector dx(x1 - x0);
662 double alpha = -((x0[0] - xwire) * dx[0] + x0[2] * dx[2]) /
663 (dx[0] * dx[0] + dx[2] * dx[2] + 1e-99);
664 alpha = (alpha < 0)? 0 : (alpha > 1)? 1 : alpha;
665 double t = (splits[0].t0_ns + splits[0].t1_ns)*ns / 2;
666 G4ThreeVector x((splits[0].x0_g + splits[0].x1_g) / 2);
667 G4ThreeVector xlocal(x0 + alpha * dx);
668 double tdrift;
669 int wire_fired = add_anode_hit(hits, splits[0], layer,
670 xwire, x, xlocal, dE,
671 t, tdrift);
672 if (wire_fired) {
673 add_cathode_hit(splits[0], packNo, xwire, xlocal[1],
674 tdrift, n_p, chamber, module, layer,
675 global_wire_number);
676 }
677 }
678 else {
679 G4ThreeVector xlocal;
680 G4ThreeVector x((splits[0].x0_g + splits[0].x1_g) / 2);
681 double t = (splits[0].t0_ns + splits[0].t1_ns)*ns / 2;
682 // Loop over the number of primary ion pairs
683 for (int n=0; n < n_p; n++) {
684 // Generate a cluster at a random position
685 // along the path within the cell
686 double u = G4UniformRand()G4MTHepRandom::getTheEngine()->flat();
687 xlocal = x0 + u * (x1 - x0);
688 double tdrift;
689 int wire_fired = add_anode_hit(hits, splits[0], layer,
690 xwire, x, xlocal, dE,
691 t, tdrift);
692 if (wire_fired) {
693 add_cathode_hit(splits[0], packNo, xwire, xlocal[1],
694 tdrift, n_p, chamber, module, layer,
695 global_wire_number);
696 }
697 }
698 }
699 }
700 splits.erase(splits.begin());
701 }
702
703 if (fDrift_clusters) {
704 // store waveform data in sampled sequence with 1 ns bins
705 int num_samples = (int)FDC_TIME_WINDOW;
706 double *samples = new double[num_samples];
707 for (int i=0; i < num_samples; i++) {
708 samples[i] = fdc_wire_signal_mV(double(i), hits);
709 }
710
711 // take the earliest hit to identify the track parameters
712 double dradius_cm = hits[0].d_cm;
713 int ptype_G3 = hits[0].ptype_G3;
714 int itrack = hits[0].itrack_;
715 double t0_ns = hits[0].t0_ns;
716 double t1_ns = hits[0].t1_ns;
717 G4ThreeVector x0_g = hits[0].x0_g;
718 G4ThreeVector x0_l = hits[0].x0_l;
719 G4ThreeVector x1_g = hits[0].x1_g;
720 G4ThreeVector x1_l = hits[0].x1_l;
721
722 double dE_keV = 0;
723 int over_threshold = 0;
724 for (int i=0; i < num_samples; i++) {
725 if (samples[i] > THRESH_ANODE) {
726 if (!over_threshold) {
727 splits.push_back(GlueXHitFDCwire::hitinfo_t());
728 splits.back().d_cm = dradius_cm;
729 splits.back().ptype_G3 = ptype_G3;
730 splits.back().itrack_ = itrack;
731 splits.back().t_ns = double(i);
732 splits.back().t0_ns = t0_ns;
733 splits.back().t1_ns = t1_ns;
734 splits.back().x0_g = x0_g;
735 splits.back().x0_l = x0_l;
736 splits.back().x1_g = x1_g;
737 splits.back().x1_l = x1_l;
738
739 // Do an interpolation to find the time
740 // at which the threshold was crossed.
741 double t_array[4];
742 double t_ns, t_err;
743 for (int k=0; k < 4; k++)
744 t_array[k] = i - 1 + k;
745 polint(&samples[i-1], t_array, 4,
746 THRESH_ANODE, &t_ns, &t_err);
747 splits.back().t_ns = t_ns;
748 over_threshold = 1;
749 }
750 dE_keV += samples[i];
751 }
752 else if (over_threshold) {
753 splits.back().dE_keV = dE_keV;
754 over_threshold = 0;
755 dE_keV = 0;
756 }
757 }
758 delete [] samples;
759 }
760 else {
761
762 // Build new reduced hit list ordered by hit time
763
764 std::vector<GlueXHitFDCwire::hitinfo_t>::iterator hiter;
765 for (hiter = hits.begin(); hiter != hits.end(); ++hiter) {
766 int merge_hits = 0;
767 std::vector<GlueXHitFDCwire::hitinfo_t>::iterator siter;
768 for (siter = splits.begin(); siter != splits.end(); ++siter) {
769 if (fabs(hiter->t_ns - siter->t_ns) < TWO_HIT_TIME_RESOL) {
770 merge_hits = 1;
771 break;
772 }
773 else if (hiter->t_ns < siter->t_ns) {
774 break;
775 }
776 }
777 if (merge_hits) {
778 if (hiter->t_ns < siter->t_ns) {
779 double dE_keV = siter->dE_keV;
780 *siter = *hiter;
781 siter->dE_keV += dE_keV;
782 }
783 else {
784 siter->dE_keV += hiter->dE_keV;
785 }
786 }
787 else {
788 splits.insert(siter, *hiter);
789 }
790 }
791 }
792
793 if (splits.size() > 0) {
794 hddm_s::FdcChamberList chambers = forwardDC.getFdcChambers();
795 hddm_s::FdcChamberList::iterator citer;
796 for (citer = chambers.begin(); citer != chambers.end(); ++citer) {
797 if (citer->getModule() == module && citer->getLayer() == layer)
798 break;
799 }
800 if (citer == chambers.end()) {
801 chambers = forwardDC.addFdcChambers(1);
802 chambers(0).setModule(module);
803 chambers(0).setLayer(layer);
804 citer = chambers.begin();
805 }
806 hddm_s::FdcAnodeWireList anodes = citer->getFdcAnodeWires();
807 hddm_s::FdcAnodeWireList::iterator aiter;
808 for (aiter = anodes.begin(); aiter != anodes.end(); ++aiter) {
809 if (aiter->getWire() == wire)
810 break;
811 }
812 if (aiter == anodes.end()) {
813 anodes = citer->addFdcAnodeWires(1);
814 anodes(0).setWire(wire);
815 aiter = anodes.begin();
816 }
817 int hitscount = splits.size();
818 if (hitscount + aiter->getFdcAnodeTruthHits().size() > MAX_HITS) {
819 hitscount = MAX_HITS - aiter->getFdcAnodeTruthHits().size();
820 G4cerr(*G4cerr_p) << "GlueXSensitiveDetectorFDC::EndOfEvent warning: "
821 << "wire hit count exceeds max hit count " << MAX_HITS
822 << ", " << splits.size() - hitscount << " hits discarded."
823 << G4endlstd::endl;
824 }
825 for (int ih=0; ih < hitscount; ++ih) {
826 hddm_s::FdcAnodeTruthHitList thit = aiter->addFdcAnodeTruthHits(1);
827 thit(0).setDE(splits[ih].dE_keV * 1e-6);
828 thit(0).setT(splits[ih].t_ns);
829 thit(0).setD(splits[ih].d_cm);
830 thit(0).setItrack(splits[ih].itrack_);
831 thit(0).setPtype(splits[ih].ptype_G3);
832 thit(0).setT_unsmeared(splits[ih].t_unsmeared_ns);
833 }
834 }
835 }
836
837 // Collect and output the cathodeTruthHits
838
839 for (siter = strips->begin(); siter != strips->end(); ++siter) {
840 int chamber = siter->second->chamber_;
841 int planeNo = siter->second->plane_;
842 int stripNo = siter->second->strip_;
843 int module = chamber / 10;
844 int layer = chamber % 10;
845 std::vector<GlueXHitFDCcathode::hitinfo_t> &hits = siter->second->hits;
846 std::vector<GlueXHitFDCcathode::hitinfo_t>::iterator hiter;
847 if (fDrift_clusters) {
848 // store waveform data in sampled sequence with 1 ns bins
849 int num_samples = (int)FDC_TIME_WINDOW;
850 double *samples = new double[num_samples];
851 for (int i=0; i < num_samples; i++) {
852 samples[i] = fdc_cathode_signal_mV(double(i), hits);
853 }
854
855 // take the earliest hit to identify the track parameters
856 int ptype_G3 = hits[0].ptype_G3;
857 int itrack = hits[0].itrack_;
858 double u_cm = hits[0].u_cm;
859 double v_cm = hits[0].v_cm;
860
861 hits.clear();
862 int istart = 0;
863 int threshold_toggle = 0;
864 int FADC_BIN_SIZE = 1;
865 for (int i=0; i < num_samples; i += FADC_BIN_SIZE) {
866 if (samples[i] > THRESH_STRIPS) {
867 if (threshold_toggle == 0) {
868 hits.push_back(GlueXHitFDCcathode::hitinfo_t());
869 hits.back().t_ns = i;
870 hits.back().ptype_G3 = ptype_G3;
871 hits.back().itrack_ = itrack;
872 hits.back().v_cm = v_cm;
873 hits.back().u_cm = u_cm;
874 istart = (i > 0)? i-1 : 0;
875 threshold_toggle=1;
876 }
877 }
878 else if (threshold_toggle) {
879 // Find the first peak
880 for (int j = istart + 1; j < i - 1; j++) {
881 if (samples[j] > samples[j-1] && samples[j] > samples[j+1]) {
882 hits.back().q_fC = samples[j];
883 break;
884 }
885 }
886 threshold_toggle=0;
887 }
888 }
889 int i = num_samples - 1;
890 if (samples[i] > THRESH_STRIPS && threshold_toggle) {
891 for (int j = istart + 1; j < i - 1; j++) {
892 if (samples[j] > samples[j-1] && samples[j] > samples[j+1]) {
893 hits.back().q_fC = samples[j];
894 break;
895 }
896 }
897 }
898 delete [] samples;
899 }
900 else {
901 double t_ns = -1e9;
902 for (hiter = hits.begin(); hiter != hits.end(); ++hiter) {
903 if (hiter == hits.begin())
904 continue;
905 // combine separate hits that are too close together in time
906 if (fabs(hiter->t_ns - t_ns) < TWO_HIT_TIME_RESOL ||
907 hiter->q_fC == 0)
908 {
909 // Use the time from the earlier hit but add the charge
910 (hiter - 1)->q_fC += hiter->q_fC;
911 hits.erase(hiter);
912 hiter = hits.begin();
913 }
914 else {
915 t_ns = hiter->t_ns;
916 }
917 }
918 }
919
920 if (hits.size() > 0) {
921 hddm_s::FdcChamberList chambers = forwardDC.getFdcChambers();
922 hddm_s::FdcChamberList::iterator citer;
923 for (citer = chambers.begin(); citer != chambers.end(); ++citer) {
924 if (citer->getModule() == module && citer->getLayer() == layer)
925 break;
926 }
927 if (citer == chambers.end()) {
928 chambers = forwardDC.addFdcChambers(1);
929 chambers(0).setModule(module);
930 chambers(0).setLayer(layer);
931 citer = chambers.begin();
932 }
933 hddm_s::FdcCathodeStripList cathodes = citer->getFdcCathodeStrips();
934 hddm_s::FdcCathodeStripList::iterator kiter;
935 for (kiter = cathodes.begin(); kiter != cathodes.end(); ++kiter) {
936 if (kiter->getPlane() == planeNo && kiter->getStrip() == stripNo)
937 break;
938 }
939 if (kiter == cathodes.end()) {
940 cathodes = citer->addFdcCathodeStrips(1);
941 cathodes(0).setPlane(planeNo);
942 cathodes(0).setStrip(stripNo);
943 kiter = cathodes.begin();
944 }
945 int hitscount = hits.size();
946 if (hitscount + kiter->getFdcCathodeTruthHits().size() > MAX_HITS) {
947 hitscount = MAX_HITS - kiter->getFdcCathodeTruthHits().size();
948 G4cerr(*G4cerr_p) << "GlueXSensitiveDetectorFDC::EndOfEvent warning: "
949 << "cathode hit count exceeds max hit count " << MAX_HITS
950 << ", " << hits.size() - hitscount << " hits discarded."
951 << G4endlstd::endl;
952 }
953 for (int ih=0; ih < hitscount; ++ih) {
954 hddm_s::FdcCathodeTruthHitList thit = kiter->addFdcCathodeTruthHits(1);
955 thit(0).setQ(hits[ih].q_fC);
956 thit(0).setT(hits[ih].t_ns);
957 thit(0).setItrack(hits[ih].itrack_);
958 thit(0).setPtype(hits[ih].ptype_G3);
959 }
960 }
961 }
962
963 // Collect and output the fdcTruthPoints
964
965 int last_chamber = -1;
966 hddm_s::FdcTruthPoint *last_point = 0;
967 for (piter = points->begin(); piter != points->end(); ++piter) {
968 if (piter->second->chamber_ == last_chamber &&
969 piter->second->track_ == last_point->getTrack() &&
970 fabs(piter->second->t_ns - last_point->getT()) < 0.1)
971 {
972 if (piter->second->dradius_cm < last_point->getDradius())
973 last_point->setDradius(piter->second->dradius_cm);
974 else
975 piter->second->dradius_cm = last_point->getDradius();
976 if (piter->second->t_ns > last_point->getT())
977 continue;
978 }
979 int module = piter->second->chamber_ / 10;
980 int layer = piter->second->chamber_ % 10;
981 hddm_s::FdcChamberList chambers = forwardDC.getFdcChambers();
982 hddm_s::FdcChamberList::iterator citer;
983 for (citer = chambers.begin(); citer != chambers.end(); ++citer) {
984 if (citer->getModule() == module && citer->getLayer() == layer)
985 break;
986 }
987 if (citer == chambers.end()) {
988 chambers = forwardDC.addFdcChambers(1);
989 chambers(0).setModule(module);
990 chambers(0).setLayer(layer);
991 citer = chambers.begin();
992 }
993 hddm_s::FdcTruthPointList point = citer->addFdcTruthPoints(1);
994 point(0).setE(piter->second->E_GeV);
995 point(0).setDEdx(piter->second->dEdx_GeV_cm);
996 point(0).setDradius(piter->second->dradius_cm);
997 point(0).setPrimary(piter->second->primary_);
998 point(0).setPtype(piter->second->ptype_G3);
999 point(0).setPx(piter->second->px_GeV);
1000 point(0).setPy(piter->second->py_GeV);
1001 point(0).setPz(piter->second->pz_GeV);
1002 point(0).setT(piter->second->t_ns);
1003 point(0).setX(piter->second->x_cm);
1004 point(0).setY(piter->second->y_cm);
1005 point(0).setZ(piter->second->z_cm);
1006 point(0).setTrack(piter->second->track_);
1007 hddm_s::TrackIDList tid = point(0).addTrackIDs();
1008 tid(0).setItrack(piter->second->trackID_);
1009 last_chamber = piter->second->chamber_;
1010 last_point = &point(0);
1011 }
1012}
1013
1014double GlueXSensitiveDetectorFDC::asic_response(double t_ns)
1015{
1016 // Simulation of the ASIC response to a pulse due to a cluster
1017
1018 double par[11] = {-0.01986, 0.01802, -0.001097, 10.3, 11.72,
1019 -0.03701, 35.84, 15.93, 0.006141, 80.95, 24.77};
1020 if (t_ns < par[3])
1021 return par[0] * t_ns + par[1] * t_ns*t_ns + par[2] * t_ns*t_ns*t_ns;
1022 else
1023 return (par[0] * par[3] +
1024 par[1] * par[3]*par[3] +
1025 par[2] * par[3]*par[3]*par[3]) *
1026 exp(-pow((t_ns - par[3])/par[4], 2)) +
1027 par[5] * exp(-pow((t_ns - par[6])/par[7], 2)) +
1028 par[8] * exp(-pow((t_ns - par[9])/par[10], 2));
1029}
1030
1031double GlueXSensitiveDetectorFDC::fdc_wire_signal_mV(
1032 double t_ns, std::vector<GlueXHitFDCwire::hitinfo_t> &hits)
1033{
1034 // Simulation of signal on a wire
1035
1036 double asic_gain = 0.76; // mV/fC
1037 double signal_mV = 0;
1038 std::vector<GlueXHitFDCwire::hitinfo_t>::iterator hiter;
1039 for (hiter = hits.begin(); hiter != hits.end(); ++hiter) {
1040 if (t_ns > hiter->t_ns) {
1041 double my_time_ns = t_ns - hiter->t_ns;
1042 signal_mV += asic_gain * hiter->dE_keV * asic_response(my_time_ns);
1043 }
1044 }
1045 return signal_mV;
1046}
1047
1048double GlueXSensitiveDetectorFDC::fdc_cathode_signal_mV(
1049 double t_ns, std::vector<GlueXHitFDCcathode::hitinfo_t> &hits)
1050{
1051 // Simulation of signal on a cathode strip (ASIC output)
1052
1053 double asic_gain = 2.3; // mv/fC
1054 double signal_mV = 0;
1055 std::vector<GlueXHitFDCcathode::hitinfo_t>::iterator hiter;
1056 for (hiter =hits.begin(); hiter != hits.end(); ++hiter) {
1057 if (t_ns > hiter->t_ns) {
1058 double my_time_ns = t_ns - hiter->t_ns;
1059 signal_mV += asic_gain * hiter->q_fC * asic_response(my_time_ns);
1060 }
1061 }
1062 return signal_mV;
1063}
1064
1065void GlueXSensitiveDetectorFDC::add_cathode_hit(
1066 GlueXHitFDCwire::hitinfo_t &wirehit,
1067 int packageNo,
1068 double xwire,
1069 double yavalanche,
1070 double tdrift,
1071 int n_p,
1072 int chamber,
1073 int module,
1074 int layer,
1075 int global_wire_number)
1076{
1077 double q_anode;
1078 if (!fDrift_clusters) {
1079 // Total number of ion pairs. On average for each primary ion
1080 // pair produced there are n_s secondary ion pairs produced. The
1081 // probability distribution is a compound poisson distribution
1082 // that requires generating two Poisson variables.
1083 double n_s_mean = n_p * N_SECOND_PER_PRIMARY;
1084 int n_s = CLHEP::RandPoisson::shoot(n_s_mean);
1085 q_anode = (n_s + n_p) * GAS_GAIN * ELECTRON_CHARGE;
1086 }
1087 else {
1088 // Distribute the number of secondary ionizations for this primary
1089 // ionization according to a Poisson distribution with mean n_s_over_p.
1090 // For simplicity we assume these secondary electrons and the primary
1091 // electron stay together as a cluster.
1092 int n_s = CLHEP::RandPoisson::shoot(N_SECOND_PER_PRIMARY);
1093 q_anode = (1 + n_s) * GAS_GAIN * ELECTRON_CHARGE;
1094 }
1095
1096 // Mock-up of cathode strip charge distribution
1097 for (int plane=1; plane < 4; plane += 2) {
1098 double theta = (plane == 1)? M_PI3.14159265358979323846-CATHODE_ROT_ANGLE : CATHODE_ROT_ANGLE;
1099 double cathode_u = -xwire * cos(theta) - yavalanche * sin(theta);
1100 int strip1 = ceil((cathode_u - U_OF_STRIP_ONE) / STRIP_SPACING + 0.5);
1101 double cathode_u1 = (strip1 - 1) * STRIP_SPACING + U_OF_STRIP_ONE;
1102 double delta_u = cathode_u - cathode_u1;
1103 for (int node = -STRIP_NODES; node <= STRIP_NODES; node++) {
1104 // Induce charge on the strips according to the Mathieson
1105 // function tuned to results from FDC prototype
1106 double lambda1 = ((node - 0.5) * STRIP_SPACING +
1107 STRIP_GAP / 2. - delta_u) / ANODE_CATHODE_SPACING;
1108 double lambda2 = ((node + 0.5) * STRIP_SPACING -
1109 STRIP_GAP / 2. - delta_u) / ANODE_CATHODE_SPACING;
1110 double factor = 0.25 * M_PI3.14159265358979323846 * K2;
1111 double q = 0.25 * q_anode * (tanh(factor * lambda2) -
1112 tanh(factor * lambda1));
1113 int strip = strip1 + node;
1114 // Throw away hits on strips falling within a certain dead-zone radius
1115 double strip_outer_u = cathode_u1;
1116 strip_outer_u += node * (STRIP_SPACING + STRIP_GAP / 2.);
1117 double cathode_v = -xwire * sin(theta) + yavalanche * cos(theta);
1118 double check_radius = sqrt(strip_outer_u * strip_outer_u +
1119 cathode_v * cathode_v);
1120 if (strip > 0 && strip <= STRIPS_PER_PLANE &&
1121 check_radius > strip_dead_zone_radius[packageNo])
1122 {
1123 int key = GlueXHitFDCcathode::GetKey(chamber, plane, strip);
1124 GlueXHitFDCcathode *cathode = (*fCathodesMap)[key];
1125 if (cathode == 0) {
1126 GlueXHitFDCcathode newcathode(chamber, plane, strip);
1127 fCathodesMap->add(key, newcathode);
1128 cathode = (*fCathodesMap)[key];
1129 }
1130 std::vector<GlueXHitFDCcathode::hitinfo_t>::iterator hiter;
1131 for (hiter = cathode->hits.begin();
1132 hiter != cathode->hits.end(); ++hiter)
1133 {
1134 if (hiter->t_ns*ns > tdrift)
1135 break;
1136 }
1137 // create new hit
1138 hiter = cathode->hits.insert(hiter, GlueXHitFDCcathode::hitinfo_t());
1139 hiter->t_ns = tdrift/ns;
1140 hiter->q_fC = q/fC;
1141 hiter->itrack_ = wirehit.itrack_;
1142 hiter->ptype_G3 = wirehit.ptype_G3;
1143 hiter->u_cm = xwire/cm;
1144 hiter->v_cm = yavalanche/cm;
1145 }
1146 } // loop over cathode strips
1147 } // loop over cathode planes
1148}
1149
1150int GlueXSensitiveDetectorFDC::add_anode_hit(
1151 std::vector<GlueXHitFDCwire::hitinfo_t> &hits,
1152 GlueXHitFDCwire::hitinfo_t &wirehit,
1153 int layer,
1154 double xwire,
1155 G4ThreeVector &xglobal,
1156 G4ThreeVector &xlocal,
1157 double dE,
1158 double t,
1159 double &tdrift)
1160{
1161 // Get the magnetic field at this cluster position
1162 G4ThreeVector B = GlueXDetectorConstruction::GetInstance()
1163 ->GetMagneticField(xglobal, tesla);
1164 double BrhoT = B.perp();
1165
1166 // Find the angle between the wire direction and the direction of the
1167 // magnetic field in the x-y plane
1168 double wire_theta = 1.0472 * ((layer % 3) - 1);
1169 double wire_dir[2] = {sin(wire_theta), cos(wire_theta)};
1170 double phi = 0;
1171 if (BrhoT > 0)
1172 phi = acos((B[0] * wire_dir[0] + B[1] * wire_dir[1]) / BrhoT);
1173
1174 // useful combinations of dx and dz
1175 double dx = xlocal[0] - xwire;
1176 double dx2 = dx * dx;
1177 double dx4 = dx2 * dx2;
1178 double dz = xlocal[2];
1179 double dz2 = dz * dz;
1180 double dz4 = dz2 * dz2;
1181
1182 // Next compute the avalanche position along wire.
1183 // Correct avalanche position with deflection along wire
1184 // due to the Lorentz force.
1185 double cm2 = cm * cm;
1186 double cm4 = cm2 * cm2;
1187 xlocal[1] += (LORENTZ_NR_PAR1 * B[2] * (1 + LORENTZ_NR_PAR2 * BrhoT)) * dx +
1188 (LORENTZ_NZ_PAR1 + LORENTZ_NZ_PAR2 * B[2]) * BrhoT * cos(phi) * xlocal[2] +
1189 (-0.000176 * dx * dx2 / (dz2 + 0.001*cm2));
1190 // Add transverse diffusion
1191 xlocal[1] += G4RandGaussG4MTRandGaussQ::shoot() *
1192 (0.01*cm * pow((dx2 + dz2)/cm2, 0.125) + 0.0061*cm * dx2/cm2);
1193
1194 // Do not use this cluster if the Lorentz force would deflect
1195 // the electrons outside the active region of the detector
1196 if (sqrt(xlocal[1] * xlocal[1] + xwire * xwire) > ACTIVE_AREA_OUTER_RADIUS)
1197 return 0;
1198
1199 // Model the drift time and longitudinal diffusion as a function of
1200 // position of the cluster within the cell
1201
1202#if OLD_FDC_DRIFT_TIME_MODEL
1203 double tdrift_unsmeared = 1086.0*ns * (1 + 0.039 * B.mag()) * dx2/cm2 +
1204 1068.0*ns * dz2/cm2 +
1205 (-2.675*ns / (dz2/cm2 + 0.001) +
1206 2.4e4*ns * dz2/cm2) * dx4/cm4;
1207#else
1208 double dradius = sqrt(dx2 + dz2);
1209 int index = locate(drift_table_d_cm, drift_table_len, dradius/cm);
1210 index = (index < drift_table_len - 3)? index : drift_table_len - 3;
1211 double *dd = &drift_table_d_cm[index];
1212 double tt = 0.5; //ns
1213 double dd10 = dd[1] - dd[0];
1214 double dd20 = dd[2] - dd[0];
1215 double dd21 = dd[2] - dd[1];
1216 double qa = tt*index;
1217 double qb = (dd20/dd10 - 2*dd10/dd20) * tt/dd21;
1218 double qc = (2/dd20 - 1/dd10) * tt/dd21;
1219 double d0 = dradius/cm - dd[0];
1220 double tdrift_unsmeared = qa + qb*d0 + qc*d0*d0;
1221#endif
1222
1223 // Apply small B-field dependence on the drift time
1224 tdrift_unsmeared *= 1. + DRIFT_BSCALE_PAR1 + DRIFT_BSCALE_PAR2*B[2]*B[2];
1225
1226 // Minimum drift time for docas near wire (very crude approximation)
1227 double v_max = 0.08*cm/ns; // guess for now based on Garfield, near wire
1228 double tmin = dradius / v_max;
1229
1230 // longitidinal diffusion, derived from Garfield calculations
1231 double dt = (G4RandGaussG4MTRandGaussQ::shoot() - 0.5) *
1232 (39.44*ns * dx4/cm4 / (0.5 - dz2/cm2) +
1233 56.00*ns * dz4/cm4 / (0.5 - dx2/cm2) +
1234 0.01566*ns * dx4/cm4 / (dz4/cm4 + 0.002) / (0.251 - dx2/cm2));
1235
1236 double tdrift_smeared = tdrift_unsmeared + dt;
1237 if (tdrift_smeared < tmin) {
1238 tdrift_smeared = tmin;
1239 }
1240
1241 // Avalanche time
1242 tdrift = t + tdrift_smeared;
1243
1244 // Skip cluster if the time would go beyond readout window
1245 if (tdrift > FDC_TIME_WINDOW)
1246 return 0;
1247
1248 // Record the anode hit
1249 std::vector<GlueXHitFDCwire::hitinfo_t>::iterator hiter;
1250 for (hiter = hits.begin(); hiter != hits.end(); ++hiter) {
1251 if (hiter->t_ns*ns > tdrift) {
1252 break;
1253 }
1254 }
1255 hiter = hits.insert(hiter, wirehit);
1256 hiter->dE_keV = dE/keV;
1257 hiter->t_ns = tdrift/ns;
1258 hiter->t_unsmeared_ns = tdrift_unsmeared/ns;
1259 hiter->d_cm = dradius/cm;
1260 return 1;
1261}
1262
1263void GlueXSensitiveDetectorFDC::polint(double *xa, double *ya, int n,
1264 double x, double *y, double *dy)
1265{
1266 // Slightly modified versions of the "polint" routine from
1267 // Press, William H., Brian P. Flannery, Saul A. Teukolsky and
1268 // William T. Vetterling, 1986, "Numerical Recipes: The Art of
1269 // Scientific Computing" (Fortran), Cambrigde University Press,
1270 // pp. 80-82.
1271
1272 double *c = NULL__null;
1273 double *d = NULL__null;
1274 double den;
1275 double dif;
1276 double dift;
1277 double ho;
1278 double hp;
1279 double w;
1280
1281 int i;
1282 int m;
1283 int ns;
1284
1285 if ((c = (double*)malloc(n*sizeof(double))) == NULL__null ||
1
Taking false branch
→
1286 (d = (double*)malloc(n*sizeof(double))) == NULL__null )
1287 {
1288 fprintf(stderrstderr, "polint error: allocating workspace\n" );
1289 fprintf(stderrstderr, "polint error: setting y = 0 and dy = 1e9\n" );
1290 *y = 0.0;
1291 *dy = 1.e9;
1292 if (c != NULL__null)
1293 free(c);
1294 if (d != NULL__null )
1295 free(d);
1296 return;
1297 }
1298
1299 ns = 0;
1300 dif = fabs(x-xa[0]);
1301 for (i = 0; i < n; ++i) {
2
←
Assuming 'i' is >= 'n'
→
3
←
Loop condition is false. Execution continues on line 1310
→
1302 dift = fabs(x-xa[i]);
1303 if (dift < dif) {
1304 ns = i;
1305 dif = dift;
1306 }
1307 c[i] = ya[i];
1308 d[i] = ya[i];
1309 }
1310 *y = ya[ns];
1311 ns = ns-1;
1312 for (m = 0; m < n-1; ++m) {
4
←
Loop condition is true. Entering loop body
→
1313 for (i = 0; i < n-m-1; ++i) {
5
←
Loop condition is true. Entering loop body
→
1314 ho = xa[i]-x;
1315 hp = xa[i+m+1]-x;
1316 w = c[i+1]-d[i];
6
←
The left operand of '-' is a garbage value
1317 den = ho-hp;
1318 if (den == 0) {
1319 fprintf( stderrstderr, "polint error: den = 0\n" );
1320 fprintf( stderrstderr, "polint error: setting y = 0 and dy = 1e9\n" );
1321 *y = 0.0;
1322 *dy = 1.e9;
1323 if (c != NULL__null)
1324 free(c);
1325 if (d != NULL__null)
1326 free(d);
1327 return;
1328 }
1329 den = w/den;
1330 d[i] = hp*den;
1331 c[i] = ho*den;
1332 }
1333 if (2*(ns+1) < n-m-1) {
1334 *dy = c[ns+1];
1335 }
1336 else {
1337 *dy = d[ns];
1338 ns = ns-1;
1339 }
1340 *y = (*y)+(*dy);
1341 }
1342
1343 if (c != NULL__null)
1344 free(c);
1345 if (d != NULL__null)
1346 free(d);
1347}
1348
1349int GlueXSensitiveDetectorFDC::locate(double *xx, int n, double x)
1350{
1351 // Locate a position in array xx of value x
1352
1353 int j;
1354 int ju;
1355 int jm;
1356 int jl;
1357 int ascnd;
1358
1359 jl = -1;
1360 ju = n;
1361 ascnd = (xx[n-1] >= xx[0]);
1362 while (ju - jl > 1) {
1363 jm = (ju + jl) >> 1;
1364 if ((x >= xx[jm]) == ascnd)
1365 jl = jm;
1366 else
1367 ju = jm;
1368 }
1369 if (x == xx[0])
1370 j = 0;
1371 else if (x == xx[n-1])
1372 j = n-2;
1373 else
1374 j = jl;
1375 return j;
1376}
1377
1378int GlueXSensitiveDetectorFDC::GetIdent(std::string div,
1379 const G4VTouchable *touch)
1380{
1381 const HddsG4Builder* bldr = GlueXDetectorConstruction::GetBuilder();
1382 std::map<std::string, std::vector<int> >::const_iterator iter;
1383 std::map<std::string, std::vector<int> > *identifiers;
1384 int max_depth = touch->GetHistoryDepth();
1385 for (int depth = 0; depth < max_depth; ++depth) {
1386 G4VPhysicalVolume *pvol = touch->GetVolume(depth);
1387 G4LogicalVolume *lvol = pvol->GetLogicalVolume();
1388 int volId = fVolumeTable[lvol];
1389 if (volId == 0) {
1390 volId = bldr->getVolumeId(lvol);
1391 fVolumeTable[lvol] = volId;
1392 }
1393 identifiers = &Refsys::fIdentifierTable[volId];
1394 if ((iter = identifiers->find(div)) != identifiers->end()) {
1395 int copyNum = touch->GetCopyNumber(depth);
1396 copyNum += (dynamic_cast<G4PVPlacement*>(pvol))? -1 : 0;
1397 return iter->second[copyNum];
1398 }
1399 }
1400 return -1;
1401}