Bug Summary

File:scratch/gluex/scan-build-work/hdgeant4/src/GlueXSensitiveDetectorFDC.cc
Location:line 397, column 7
Description:Value stored to 't' is never read

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 total number of allowed anode and cathode hits
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
56int 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;
Value stored to 't' is never read
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 if ((int)anode->hits.size() < MAX_HITS) {
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 else {
557 G4cerr(*G4cerr_p) << "GlueXSensitiveDetectorFDC::ProcessHits error: "
558 << "max hit count " << MAX_HITS << " exceeded, truncating!"
559 << G4endlstd::endl;
560 }
561 }
562 }
563 return true;
564}
565
566void GlueXSensitiveDetectorFDC::EndOfEvent(G4HCofThisEvent*)
567{
568 std::map<int,GlueXHitFDCwire*> *wires = fWiresMap->GetMap();
569 std::map<int,GlueXHitFDCcathode*> *strips = fCathodesMap->GetMap();
570 std::map<int,GlueXHitFDCpoint*> *points = fPointsMap->GetMap();
571 if (wires->size() == 0 && strips->size() == 0 && points->size() == 0)
572 return;
573 std::map<int,GlueXHitFDCwire*>::iterator witer;
574 std::map<int,GlueXHitFDCcathode*>::iterator siter;
575 std::map<int,GlueXHitFDCpoint*>::iterator piter;
576
577 if (verboseLevel > 1) {
578 G4cout(*G4cout_p) << G4endlstd::endl
579 << "--------> Hits Collection: in this event there are "
580 << wires->size() << " anode wires with hits in the FDC: "
581 << G4endlstd::endl;
582 for (witer = wires->begin(); witer != wires->end(); ++witer)
583 witer->second->Print();
584
585 G4cout(*G4cout_p) << G4endlstd::endl
586 << "--------> Hits Collection: in this event there are "
587 << strips->size() << " cathode strips with hits in the FDC: "
588 << G4endlstd::endl;
589 for (siter = strips->begin(); siter != strips->end(); ++siter)
590 siter->second->Print();
591
592 G4cout(*G4cout_p) << G4endlstd::endl
593 << "--------> Hits Collection: in this event there are "
594 << points->size() << " truth points in the FDC: "
595 << G4endlstd::endl;
596 for (piter = points->begin(); piter != points->end(); ++piter)
597 piter->second->Print();
598 }
599
600 // pack hits into ouptut hddm record
601
602 G4EventManager* mgr = G4EventManager::GetEventManager();
603 G4VUserEventInformation* info = mgr->GetUserInformation();
604 hddm_s::HDDM *record = ((GlueXUserEventInformation*)info)->getOutputRecord();
605 if (record == 0) {
606 G4cerr(*G4cerr_p) << "GlueXSensitiveDetectorFDC::EndOfEvent error - "
607 << "hits seen but no output hddm record to save them into, "
608 << "cannot continue!" << G4endlstd::endl;
609 exit(1);
610 }
611
612 if (record->getPhysicsEvents().size() == 0)
613 record->addPhysicsEvents();
614 if (record->getHitViews().size() == 0)
615 record->getPhysicsEvent().addHitViews();
616 hddm_s::HitView &hitview = record->getPhysicsEvent().getHitView();
617 if (hitview.getForwardDCs().size() == 0)
618 hitview.addForwardDCs();
619 hddm_s::ForwardDC &forwardDC = hitview.getForwardDC();
620
621 // Collect and output the wireTruthHits
622
623 for (witer = wires->begin(); witer != wires->end(); ++witer) {
624
625 // Merge multiple segments from a single track into one, and
626 // apply the drift time algorithm to get a single hit time for each.
627
628 int chamber = witer->second->chamber_;
629 int wire = witer->second->wire_;
630 int module = chamber / 10;
631 int packNo = (module - 1) / 2;
632 int layer = chamber % 10;
633 int glayerNo = 3*(layer - 1) + module-1;
634 int global_wire_number = 96 * glayerNo + wire - 1;
635 std::vector<GlueXHitFDCwire::hitinfo_t> &splits = witer->second->hits;
636 std::vector<GlueXHitFDCwire::hitinfo_t> hits;
637 while (splits.size() > 0) {
638 for (unsigned int ih=1; ih < splits.size(); ++ih) {
639 if (fabs(splits[ih].t0_ns - splits[0].t1_ns) < 0.5*ns) {
640 splits[0].dE_keV += splits[ih].dE_keV;
641 splits[0].t1_ns = splits[ih].t1_ns;
642 splits[0].x1_g = splits[ih].x1_g;
643 splits[0].x1_l = splits[ih].x1_l;
644 splits.erase(splits.begin() + ih);
645 --ih;
646 }
647 }
648
649 // Simulate the number of primary ion pairs.
650 // The total number of ion pairs depends on the energy deposition
651 // and the effective average energy to produce a pair, w_eff.
652 // On average for each primary ion pair produced there are n_s_per_p
653 // secondary ion pairs produced.
654
655 double xwire = U_OF_WIRE_ONE + (wire - 1) * WIRE_SPACING;
656 double dE = splits[0].dE_keV*keV;
657 if (dE > THRESH_KEV*keV) {
658 // Average number of primary ion pairs
659 double n_p_mean = dE / W_EFF_PER_ION / (1 + N_SECOND_PER_PRIMARY);
660 // number of primary ion pairs
661 int n_p = CLHEP::RandPoisson::shoot(n_p_mean);
662 G4ThreeVector &x0 = splits[0].x0_l;
663 G4ThreeVector &x1 = splits[0].x1_l;
664
665 if (fDrift_clusters == 0) {
666 G4ThreeVector dx(x1 - x0);
667 double alpha = -((x0[0] - xwire) * dx[0] + x0[2] * dx[2]) /
668 (dx[0] * dx[0] + dx[2] * dx[2] + 1e-99);
669 alpha = (alpha < 0)? 0 : (alpha > 1)? 1 : alpha;
670 double t = (splits[0].t0_ns + splits[0].t1_ns)*ns / 2;
671 G4ThreeVector x((splits[0].x0_g + splits[0].x1_g) / 2);
672 G4ThreeVector xlocal(x0 + alpha * dx);
673 double tdrift;
674 int wire_fired = add_anode_hit(hits, splits[0], layer,
675 xwire, x, xlocal, dE,
676 t, tdrift);
677 if (wire_fired) {
678 add_cathode_hit(splits[0], packNo, xwire, xlocal[1],
679 tdrift, n_p, chamber, module, layer,
680 global_wire_number);
681 }
682 }
683 else {
684 G4ThreeVector xlocal;
685 G4ThreeVector x((splits[0].x0_g + splits[0].x1_g) / 2);
686 double t = (splits[0].t0_ns + splits[0].t1_ns)*ns / 2;
687 // Loop over the number of primary ion pairs
688 for (int n=0; n < n_p; n++) {
689 // Generate a cluster at a random position
690 // along the path within the cell
691 double u = G4UniformRand()G4MTHepRandom::getTheEngine()->flat();
692 xlocal = x0 + u * (x1 - x0);
693 double tdrift;
694 int wire_fired = add_anode_hit(hits, splits[0], layer,
695 xwire, x, xlocal, dE,
696 t, tdrift);
697 if (wire_fired) {
698 add_cathode_hit(splits[0], packNo, xwire, xlocal[1],
699 tdrift, n_p, chamber, module, layer,
700 global_wire_number);
701 }
702 }
703 }
704 }
705 splits.erase(splits.begin());
706 }
707
708 if (fDrift_clusters) {
709 // store waveform data in sampled sequence with 1 ns bins
710 int num_samples = (int)FDC_TIME_WINDOW;
711 double *samples = new double[num_samples];
712 for (int i=0; i < num_samples; i++) {
713 samples[i] = fdc_wire_signal_mV(double(i), hits);
714 }
715
716 // take the earliest hit to identify the track parameters
717 double dradius_cm = hits[0].d_cm;
718 int ptype_G3 = hits[0].ptype_G3;
719 int itrack = hits[0].itrack_;
720 double t0_ns = hits[0].t0_ns;
721 double t1_ns = hits[0].t1_ns;
722 G4ThreeVector x0_g = hits[0].x0_g;
723 G4ThreeVector x0_l = hits[0].x0_l;
724 G4ThreeVector x1_g = hits[0].x1_g;
725 G4ThreeVector x1_l = hits[0].x1_l;
726
727 double dE_keV = 0;
728 int over_threshold = 0;
729 for (int i=0; i < num_samples; i++) {
730 if (samples[i] > THRESH_ANODE) {
731 if (!over_threshold) {
732 splits.push_back(GlueXHitFDCwire::hitinfo_t());
733 splits.back().d_cm = dradius_cm;
734 splits.back().ptype_G3 = ptype_G3;
735 splits.back().itrack_ = itrack;
736 splits.back().t_ns = double(i);
737 splits.back().t0_ns = t0_ns;
738 splits.back().t1_ns = t1_ns;
739 splits.back().x0_g = x0_g;
740 splits.back().x0_l = x0_l;
741 splits.back().x1_g = x1_g;
742 splits.back().x1_l = x1_l;
743
744 // Do an interpolation to find the time
745 // at which the threshold was crossed.
746 double t_array[4];
747 double t_ns, t_err;
748 for (int k=0; k < 4; k++)
749 t_array[k] = i - 1 + k;
750 polint(&samples[i-1], t_array, 4,
751 THRESH_ANODE, &t_ns, &t_err);
752 splits.back().t_ns = t_ns;
753 over_threshold = 1;
754 }
755 dE_keV += samples[i];
756 }
757 else if (over_threshold) {
758 splits.back().dE_keV = dE_keV;
759 over_threshold = 0;
760 dE_keV = 0;
761 }
762 }
763 delete [] samples;
764 }
765 else {
766
767 // Build new reduced hit list ordered by hit time
768
769 std::vector<GlueXHitFDCwire::hitinfo_t>::iterator hiter;
770 for (hiter = hits.begin(); hiter != hits.end(); ++hiter) {
771 int merge_hits = 0;
772 std::vector<GlueXHitFDCwire::hitinfo_t>::iterator siter;
773 for (siter = splits.begin(); siter != splits.end(); ++siter) {
774 if (fabs(hiter->t_ns - siter->t_ns) < TWO_HIT_TIME_RESOL) {
775 merge_hits = 1;
776 break;
777 }
778 else if (hiter->t_ns < siter->t_ns) {
779 break;
780 }
781 }
782 if (merge_hits) {
783 if (hiter->t_ns < siter->t_ns) {
784 double dE_keV = siter->dE_keV;
785 *siter = *hiter;
786 siter->dE_keV += dE_keV;
787 }
788 else {
789 siter->dE_keV += hiter->dE_keV;
790 }
791 }
792 else {
793 splits.insert(siter, *hiter);
794 }
795 }
796 }
797
798 if (splits.size() > 0) {
799 hddm_s::FdcChamberList chambers = forwardDC.getFdcChambers();
800 hddm_s::FdcChamberList::iterator citer;
801 for (citer = chambers.begin(); citer != chambers.end(); ++citer) {
802 if (citer->getModule() == module && citer->getLayer() == layer)
803 break;
804 }
805 if (citer == chambers.end()) {
806 chambers = forwardDC.addFdcChambers(1);
807 chambers(0).setModule(module);
808 chambers(0).setLayer(layer);
809 citer = chambers.begin();
810 }
811 hddm_s::FdcAnodeWireList anodes = citer->getFdcAnodeWires();
812 hddm_s::FdcAnodeWireList::iterator aiter;
813 for (aiter = anodes.begin(); aiter != anodes.end(); ++aiter) {
814 if (aiter->getWire() == wire)
815 break;
816 }
817 if (aiter == anodes.end()) {
818 anodes = citer->addFdcAnodeWires(1);
819 anodes(0).setWire(wire);
820 aiter = anodes.begin();
821 }
822 for (int ih=0; ih < (int)splits.size(); ++ih) {
823 hddm_s::FdcAnodeTruthHitList thit = aiter->addFdcAnodeTruthHits(1);
824 thit(0).setDE(splits[ih].dE_keV * 1e-6);
825 thit(0).setT(splits[ih].t_ns);
826 thit(0).setD(splits[ih].d_cm);
827 thit(0).setItrack(splits[ih].itrack_);
828 thit(0).setPtype(splits[ih].ptype_G3);
829 thit(0).setT_unsmeared(splits[ih].t_unsmeared_ns);
830 }
831 }
832 }
833
834 // Collect and output the cathodeTruthHits
835
836 for (siter = strips->begin(); siter != strips->end(); ++siter) {
837 int chamber = siter->second->chamber_;
838 int planeNo = siter->second->plane_;
839 int stripNo = siter->second->strip_;
840 int module = chamber / 10;
841 int layer = chamber % 10;
842 std::vector<GlueXHitFDCcathode::hitinfo_t> &hits = siter->second->hits;
843 std::vector<GlueXHitFDCcathode::hitinfo_t>::iterator hiter;
844 if (fDrift_clusters) {
845 // store waveform data in sampled sequence with 1 ns bins
846 int num_samples = (int)FDC_TIME_WINDOW;
847 double *samples = new double[num_samples];
848 for (int i=0; i < num_samples; i++) {
849 samples[i] = fdc_cathode_signal_mV(double(i), hits);
850 }
851
852 // take the earliest hit to identify the track parameters
853 int ptype_G3 = hits[0].ptype_G3;
854 int itrack = hits[0].itrack_;
855 double u_cm = hits[0].u_cm;
856 double v_cm = hits[0].v_cm;
857
858 hits.clear();
859 int istart = 0;
860 int threshold_toggle = 0;
861 int FADC_BIN_SIZE = 1;
862 for (int i=0; i < num_samples; i += FADC_BIN_SIZE) {
863 if (samples[i] > THRESH_STRIPS) {
864 if (threshold_toggle == 0) {
865 hits.push_back(GlueXHitFDCcathode::hitinfo_t());
866 hits.back().t_ns = i;
867 hits.back().ptype_G3 = ptype_G3;
868 hits.back().itrack_ = itrack;
869 hits.back().v_cm = v_cm;
870 hits.back().u_cm = u_cm;
871 istart = (i > 0)? i-1 : 0;
872 threshold_toggle=1;
873 }
874 }
875 else if (threshold_toggle) {
876 // Find the first peak
877 for (int j = istart + 1; j < i - 1; j++) {
878 if (samples[j] > samples[j-1] && samples[j] > samples[j+1]) {
879 hits.back().q_fC = samples[j];
880 break;
881 }
882 }
883 threshold_toggle=0;
884 }
885 }
886 int i = num_samples - 1;
887 if (samples[i] > THRESH_STRIPS && threshold_toggle) {
888 for (int j = istart + 1; j < i - 1; j++) {
889 if (samples[j] > samples[j-1] && samples[j] > samples[j+1]) {
890 hits.back().q_fC = samples[j];
891 break;
892 }
893 }
894 }
895 delete [] samples;
896 }
897 else {
898 double t_ns = -1e9;
899 for (hiter = hits.begin(); hiter != hits.end(); ++hiter) {
900 if (hiter == hits.begin())
901 continue;
902 // combine separate hits that are too close together in time
903 if (fabs(hiter->t_ns - t_ns) < TWO_HIT_TIME_RESOL ||
904 hiter->q_fC == 0)
905 {
906 // Use the time from the earlier hit but add the charge
907 (hiter - 1)->q_fC += hiter->q_fC;
908 hits.erase(hiter);
909 hiter = hits.begin();
910 }
911 else {
912 t_ns = hiter->t_ns;
913 }
914 }
915 }
916
917 if (hits.size() > 0) {
918 hddm_s::FdcChamberList chambers = forwardDC.getFdcChambers();
919 hddm_s::FdcChamberList::iterator citer;
920 for (citer = chambers.begin(); citer != chambers.end(); ++citer) {
921 if (citer->getModule() == module && citer->getLayer() == layer)
922 break;
923 }
924 if (citer == chambers.end()) {
925 chambers = forwardDC.addFdcChambers(1);
926 chambers(0).setModule(module);
927 chambers(0).setLayer(layer);
928 citer = chambers.begin();
929 }
930 hddm_s::FdcCathodeStripList cathodes = citer->getFdcCathodeStrips();
931 hddm_s::FdcCathodeStripList::iterator kiter;
932 for (kiter = cathodes.begin(); kiter != cathodes.end(); ++kiter) {
933 if (kiter->getPlane() == planeNo && kiter->getStrip() == stripNo)
934 break;
935 }
936 if (kiter == cathodes.end()) {
937 cathodes = citer->addFdcCathodeStrips(1);
938 cathodes(0).setPlane(planeNo);
939 cathodes(0).setStrip(stripNo);
940 kiter = cathodes.begin();
941 }
942 for (int ih=0; ih < (int)hits.size(); ++ih) {
943 hddm_s::FdcCathodeTruthHitList thit = kiter->addFdcCathodeTruthHits(1);
944 thit(0).setQ(hits[ih].q_fC);
945 thit(0).setT(hits[ih].t_ns);
946 thit(0).setItrack(hits[ih].itrack_);
947 thit(0).setPtype(hits[ih].ptype_G3);
948 }
949 }
950 }
951
952 // Collect and output the fdcTruthPoints
953
954 int last_chamber = -1;
955 hddm_s::FdcTruthPoint *last_point = 0;
956 for (piter = points->begin(); piter != points->end(); ++piter) {
957 if (piter->second->chamber_ == last_chamber &&
958 piter->second->track_ == last_point->getTrack() &&
959 fabs(piter->second->t_ns - last_point->getT()) < 0.1)
960 {
961 if (piter->second->dradius_cm < last_point->getDradius())
962 last_point->setDradius(piter->second->dradius_cm);
963 else
964 piter->second->dradius_cm = last_point->getDradius();
965 if (piter->second->t_ns > last_point->getT())
966 continue;
967 }
968 int module = piter->second->chamber_ / 10;
969 int layer = piter->second->chamber_ % 10;
970 hddm_s::FdcChamberList chambers = forwardDC.getFdcChambers();
971 hddm_s::FdcChamberList::iterator citer;
972 for (citer = chambers.begin(); citer != chambers.end(); ++citer) {
973 if (citer->getModule() == module && citer->getLayer() == layer)
974 break;
975 }
976 if (citer == chambers.end()) {
977 chambers = forwardDC.addFdcChambers(1);
978 chambers(0).setModule(module);
979 chambers(0).setLayer(layer);
980 citer = chambers.begin();
981 }
982 hddm_s::FdcTruthPointList point = citer->addFdcTruthPoints(1);
983 point(0).setE(piter->second->E_GeV);
984 point(0).setDEdx(piter->second->dEdx_GeV_cm);
985 point(0).setDradius(piter->second->dradius_cm);
986 point(0).setPrimary(piter->second->primary_);
987 point(0).setPtype(piter->second->ptype_G3);
988 point(0).setPx(piter->second->px_GeV);
989 point(0).setPy(piter->second->py_GeV);
990 point(0).setPz(piter->second->pz_GeV);
991 point(0).setT(piter->second->t_ns);
992 point(0).setX(piter->second->x_cm);
993 point(0).setY(piter->second->y_cm);
994 point(0).setZ(piter->second->z_cm);
995 point(0).setTrack(piter->second->track_);
996 hddm_s::TrackIDList tid = point(0).addTrackIDs();
997 tid(0).setItrack(piter->second->trackID_);
998 last_chamber = piter->second->chamber_;
999 last_point = &point(0);
1000 }
1001}
1002
1003double GlueXSensitiveDetectorFDC::asic_response(double t_ns)
1004{
1005 // Simulation of the ASIC response to a pulse due to a cluster
1006
1007 double par[11] = {-0.01986, 0.01802, -0.001097, 10.3, 11.72,
1008 -0.03701, 35.84, 15.93, 0.006141, 80.95, 24.77};
1009 if (t_ns < par[3])
1010 return par[0] * t_ns + par[1] * t_ns*t_ns + par[2] * t_ns*t_ns*t_ns;
1011 else
1012 return (par[0] * par[3] +
1013 par[1] * par[3]*par[3] +
1014 par[2] * par[3]*par[3]*par[3]) *
1015 exp(-pow((t_ns - par[3])/par[4], 2)) +
1016 par[5] * exp(-pow((t_ns - par[6])/par[7], 2)) +
1017 par[8] * exp(-pow((t_ns - par[9])/par[10], 2));
1018}
1019
1020double GlueXSensitiveDetectorFDC::fdc_wire_signal_mV(
1021 double t_ns, std::vector<GlueXHitFDCwire::hitinfo_t> &hits)
1022{
1023 // Simulation of signal on a wire
1024
1025 double asic_gain = 0.76; // mV/fC
1026 double signal_mV = 0;
1027 std::vector<GlueXHitFDCwire::hitinfo_t>::iterator hiter;
1028 for (hiter = hits.begin(); hiter != hits.end(); ++hiter) {
1029 if (t_ns > hiter->t_ns) {
1030 double my_time_ns = t_ns - hiter->t_ns;
1031 signal_mV += asic_gain * hiter->dE_keV * asic_response(my_time_ns);
1032 }
1033 }
1034 return signal_mV;
1035}
1036
1037double GlueXSensitiveDetectorFDC::fdc_cathode_signal_mV(
1038 double t_ns, std::vector<GlueXHitFDCcathode::hitinfo_t> &hits)
1039{
1040 // Simulation of signal on a cathode strip (ASIC output)
1041
1042 double asic_gain = 2.3; // mv/fC
1043 double signal_mV = 0;
1044 std::vector<GlueXHitFDCcathode::hitinfo_t>::iterator hiter;
1045 for (hiter =hits.begin(); hiter != hits.end(); ++hiter) {
1046 if (t_ns > hiter->t_ns) {
1047 double my_time_ns = t_ns - hiter->t_ns;
1048 signal_mV += asic_gain * hiter->q_fC * asic_response(my_time_ns);
1049 }
1050 }
1051 return signal_mV;
1052}
1053
1054void GlueXSensitiveDetectorFDC::add_cathode_hit(
1055 GlueXHitFDCwire::hitinfo_t &wirehit,
1056 int packageNo,
1057 double xwire,
1058 double yavalanche,
1059 double tdrift,
1060 int n_p,
1061 int chamber,
1062 int module,
1063 int layer,
1064 int global_wire_number)
1065{
1066 double q_anode;
1067 if (!fDrift_clusters) {
1068 // Total number of ion pairs. On average for each primary ion
1069 // pair produced there are n_s secondary ion pairs produced. The
1070 // probability distribution is a compound poisson distribution
1071 // that requires generating two Poisson variables.
1072 double n_s_mean = n_p * N_SECOND_PER_PRIMARY;
1073 int n_s = CLHEP::RandPoisson::shoot(n_s_mean);
1074 q_anode = (n_s + n_p) * GAS_GAIN * ELECTRON_CHARGE;
1075 }
1076 else {
1077 // Distribute the number of secondary ionizations for this primary
1078 // ionization according to a Poisson distribution with mean n_s_over_p.
1079 // For simplicity we assume these secondary electrons and the primary
1080 // electron stay together as a cluster.
1081 int n_s = CLHEP::RandPoisson::shoot(N_SECOND_PER_PRIMARY);
1082 q_anode = (1 + n_s) * GAS_GAIN * ELECTRON_CHARGE;
1083 }
1084
1085 // Mock-up of cathode strip charge distribution
1086 for (int plane=1; plane < 4; plane += 2) {
1087 double theta = (plane == 1)? M_PI3.14159265358979323846-CATHODE_ROT_ANGLE : CATHODE_ROT_ANGLE;
1088 double cathode_u = -xwire * cos(theta) - yavalanche * sin(theta);
1089 int strip1 = ceil((cathode_u - U_OF_STRIP_ONE) / STRIP_SPACING + 0.5);
1090 double cathode_u1 = (strip1 - 1) * STRIP_SPACING + U_OF_STRIP_ONE;
1091 double delta_u = cathode_u - cathode_u1;
1092 for (int node = -STRIP_NODES; node <= STRIP_NODES; node++) {
1093 // Induce charge on the strips according to the Mathieson
1094 // function tuned to results from FDC prototype
1095 double lambda1 = ((node - 0.5) * STRIP_SPACING +
1096 STRIP_GAP / 2. - delta_u) / ANODE_CATHODE_SPACING;
1097 double lambda2 = ((node + 0.5) * STRIP_SPACING -
1098 STRIP_GAP / 2. - delta_u) / ANODE_CATHODE_SPACING;
1099 double factor = 0.25 * M_PI3.14159265358979323846 * K2;
1100 double q = 0.25 * q_anode * (tanh(factor * lambda2) -
1101 tanh(factor * lambda1));
1102 int strip = strip1 + node;
1103 // Throw away hits on strips falling within a certain dead-zone radius
1104 double strip_outer_u = cathode_u1;
1105 strip_outer_u += node * (STRIP_SPACING + STRIP_GAP / 2.);
1106 double cathode_v = -xwire * sin(theta) + yavalanche * cos(theta);
1107 double check_radius = sqrt(strip_outer_u * strip_outer_u +
1108 cathode_v * cathode_v);
1109 if (strip > 0 && strip <= STRIPS_PER_PLANE &&
1110 check_radius > strip_dead_zone_radius[packageNo])
1111 {
1112 int key = GlueXHitFDCcathode::GetKey(chamber, plane, strip);
1113 GlueXHitFDCcathode *cathode = (*fCathodesMap)[key];
1114 if (cathode == 0) {
1115 GlueXHitFDCcathode newcathode(chamber, plane, strip);
1116 fCathodesMap->add(key, newcathode);
1117 cathode = (*fCathodesMap)[key];
1118 }
1119 std::vector<GlueXHitFDCcathode::hitinfo_t>::iterator hiter;
1120 for (hiter = cathode->hits.begin();
1121 hiter != cathode->hits.end(); ++hiter)
1122 {
1123 if (hiter->t_ns*ns > tdrift)
1124 break;
1125 }
1126 // create new hit
1127 hiter = cathode->hits.insert(hiter, GlueXHitFDCcathode::hitinfo_t());
1128 hiter->t_ns = tdrift/ns;
1129 hiter->q_fC = q/fC;
1130 hiter->itrack_ = wirehit.itrack_;
1131 hiter->ptype_G3 = wirehit.ptype_G3;
1132 hiter->u_cm = xwire/cm;
1133 hiter->v_cm = yavalanche/cm;
1134 }
1135 } // loop over cathode strips
1136 } // loop over cathode planes
1137}
1138
1139int GlueXSensitiveDetectorFDC::add_anode_hit(
1140 std::vector<GlueXHitFDCwire::hitinfo_t> &hits,
1141 GlueXHitFDCwire::hitinfo_t &wirehit,
1142 int layer,
1143 double xwire,
1144 G4ThreeVector &xglobal,
1145 G4ThreeVector &xlocal,
1146 double dE,
1147 double t,
1148 double &tdrift)
1149{
1150 // Get the magnetic field at this cluster position
1151 G4ThreeVector B = GlueXDetectorConstruction::GetInstance()
1152 ->GetMagneticField(xglobal, tesla);
1153 double BrhoT = B.perp();
1154
1155 // Find the angle between the wire direction and the direction of the
1156 // magnetic field in the x-y plane
1157 double wire_theta = 1.0472 * ((layer % 3) - 1);
1158 double wire_dir[2] = {sin(wire_theta), cos(wire_theta)};
1159 double phi = 0;
1160 if (BrhoT > 0)
1161 phi = acos((B[0] * wire_dir[0] + B[1] * wire_dir[1]) / BrhoT);
1162
1163 // useful combinations of dx and dz
1164 double dx = xlocal[0] - xwire;
1165 double dx2 = dx * dx;
1166 double dx4 = dx2 * dx2;
1167 double dz = xlocal[2];
1168 double dz2 = dz * dz;
1169 double dz4 = dz2 * dz2;
1170
1171 // Next compute the avalanche position along wire.
1172 // Correct avalanche position with deflection along wire
1173 // due to the Lorentz force.
1174 double cm2 = cm * cm;
1175 double cm4 = cm2 * cm2;
1176 xlocal[1] += (LORENTZ_NR_PAR1 * B[2] * (1 + LORENTZ_NR_PAR2 * BrhoT)) * dx +
1177 (LORENTZ_NZ_PAR1 + LORENTZ_NZ_PAR2 * B[2]) * BrhoT * cos(phi) * xlocal[2] +
1178 (-0.000176 * dx * dx2 / (dz2 + 0.001*cm2));
1179 // Add transverse diffusion
1180 xlocal[1] += G4RandGaussG4MTRandGaussQ::shoot() *
1181 (0.01*cm * pow((dx2 + dz2)/cm2, 0.125) + 0.0061*cm * dx2/cm2);
1182
1183 // Do not use this cluster if the Lorentz force would deflect
1184 // the electrons outside the active region of the detector
1185 if (sqrt(xlocal[1] * xlocal[1] + xwire * xwire) > ACTIVE_AREA_OUTER_RADIUS)
1186 return 0;
1187
1188 // Model the drift time and longitudinal diffusion as a function of
1189 // position of the cluster within the cell
1190
1191#if OLD_FDC_DRIFT_TIME_MODEL
1192 double tdrift_unsmeared = 1086.0*ns * (1 + 0.039 * B.mag()) * dx2/cm2 +
1193 1068.0*ns * dz2/cm2 +
1194 (-2.675*ns / (dz2/cm2 + 0.001) +
1195 2.4e4*ns * dz2/cm2) * dx4/cm4;
1196#else
1197 double dradius = sqrt(dx2 + dz2);
1198 int index = locate(drift_table_d_cm, drift_table_len, dradius/cm);
1199 index = (index < drift_table_len - 3)? index : drift_table_len - 3;
1200 double *dd = &drift_table_d_cm[index];
1201 double tt = 0.5; //ns
1202 double dd10 = dd[1] - dd[0];
1203 double dd20 = dd[2] - dd[0];
1204 double dd21 = dd[2] - dd[1];
1205 double qa = tt*index;
1206 double qb = (dd20/dd10 - 2*dd10/dd20) * tt/dd21;
1207 double qc = (2/dd20 - 1/dd10) * tt/dd21;
1208 double d0 = dradius/cm - dd[0];
1209 double tdrift_unsmeared = qa + qb*d0 + qc*d0*d0;
1210#endif
1211
1212 // Apply small B-field dependence on the drift time
1213 tdrift_unsmeared *= 1. + DRIFT_BSCALE_PAR1 + DRIFT_BSCALE_PAR2*B[2]*B[2];
1214
1215 // Minimum drift time for docas near wire (very crude approximation)
1216 double v_max = 0.08*cm/ns; // guess for now based on Garfield, near wire
1217 double tmin = dradius / v_max;
1218
1219 // longitidinal diffusion, derived from Garfield calculations
1220 double dt = (G4RandGaussG4MTRandGaussQ::shoot() - 0.5) *
1221 (39.44*ns * dx4/cm4 / (0.5 - dz2/cm2) +
1222 56.00*ns * dz4/cm4 / (0.5 - dx2/cm2) +
1223 0.01566*ns * dx4/cm4 / (dz4/cm4 + 0.002) / (0.251 - dx2/cm2));
1224
1225 double tdrift_smeared = tdrift_unsmeared + dt;
1226 if (tdrift_smeared < tmin) {
1227 tdrift_smeared = tmin;
1228 }
1229
1230 // Avalanche time
1231 tdrift = t + tdrift_smeared;
1232
1233 // Skip cluster if the time would go beyond readout window
1234 if (tdrift > FDC_TIME_WINDOW)
1235 return 0;
1236
1237 // Record the anode hit
1238 std::vector<GlueXHitFDCwire::hitinfo_t>::iterator hiter;
1239 for (hiter = hits.begin(); hiter != hits.end(); ++hiter) {
1240 if (hiter->t_ns*ns > tdrift) {
1241 break;
1242 }
1243 }
1244 hiter = hits.insert(hiter, wirehit);
1245 hiter->dE_keV = dE/keV;
1246 hiter->t_ns = tdrift/ns;
1247 hiter->t_unsmeared_ns = tdrift_unsmeared/ns;
1248 hiter->d_cm = dradius/cm;
1249 return 1;
1250}
1251
1252void GlueXSensitiveDetectorFDC::polint(double *xa, double *ya, int n,
1253 double x, double *y, double *dy)
1254{
1255 // Slightly modified versions of the "polint" routine from
1256 // Press, William H., Brian P. Flannery, Saul A. Teukolsky and
1257 // William T. Vetterling, 1986, "Numerical Recipes: The Art of
1258 // Scientific Computing" (Fortran), Cambrigde University Press,
1259 // pp. 80-82.
1260
1261 double *c = NULL__null;
1262 double *d = NULL__null;
1263 double den;
1264 double dif;
1265 double dift;
1266 double ho;
1267 double hp;
1268 double w;
1269
1270 int i;
1271 int m;
1272 int ns;
1273
1274 if ((c = (double*)malloc(n*sizeof(double))) == NULL__null ||
1275 (d = (double*)malloc(n*sizeof(double))) == NULL__null )
1276 {
1277 fprintf(stderrstderr, "polint error: allocating workspace\n" );
1278 fprintf(stderrstderr, "polint error: setting y = 0 and dy = 1e9\n" );
1279 *y = 0.0;
1280 *dy = 1.e9;
1281 if (c != NULL__null)
1282 free(c);
1283 if (d != NULL__null )
1284 free(d);
1285 return;
1286 }
1287
1288 ns = 0;
1289 dif = fabs(x-xa[0]);
1290 for (i = 0; i < n; ++i) {
1291 dift = fabs(x-xa[i]);
1292 if (dift < dif) {
1293 ns = i;
1294 dif = dift;
1295 }
1296 c[i] = ya[i];
1297 d[i] = ya[i];
1298 }
1299 *y = ya[ns];
1300 ns = ns-1;
1301 for (m = 0; m < n-1; ++m) {
1302 for (i = 0; i < n-m-1; ++i) {
1303 ho = xa[i]-x;
1304 hp = xa[i+m+1]-x;
1305 w = c[i+1]-d[i];
1306 den = ho-hp;
1307 if (den == 0) {
1308 fprintf( stderrstderr, "polint error: den = 0\n" );
1309 fprintf( stderrstderr, "polint error: setting y = 0 and dy = 1e9\n" );
1310 *y = 0.0;
1311 *dy = 1.e9;
1312 if (c != NULL__null)
1313 free(c);
1314 if (d != NULL__null)
1315 free(d);
1316 return;
1317 }
1318 den = w/den;
1319 d[i] = hp*den;
1320 c[i] = ho*den;
1321 }
1322 if (2*(ns+1) < n-m-1) {
1323 *dy = c[ns+1];
1324 }
1325 else {
1326 *dy = d[ns];
1327 ns = ns-1;
1328 }
1329 *y = (*y)+(*dy);
1330 }
1331
1332 if (c != NULL__null)
1333 free(c);
1334 if (d != NULL__null)
1335 free(d);
1336}
1337
1338int GlueXSensitiveDetectorFDC::locate(double *xx, int n, double x)
1339{
1340 // Locate a position in array xx of value x
1341
1342 int j;
1343 int ju;
1344 int jm;
1345 int jl;
1346 int ascnd;
1347
1348 jl = -1;
1349 ju = n;
1350 ascnd = (xx[n-1] >= xx[0]);
1351 while (ju - jl > 1) {
1352 jm = (ju + jl) >> 1;
1353 if ((x >= xx[jm]) == ascnd)
1354 jl = jm;
1355 else
1356 ju = jm;
1357 }
1358 if (x == xx[0])
1359 j = 0;
1360 else if (x == xx[n-1])
1361 j = n-2;
1362 else
1363 j = jl;
1364 return j;
1365}
1366
1367int GlueXSensitiveDetectorFDC::GetIdent(std::string div,
1368 const G4VTouchable *touch)
1369{
1370 const HddsG4Builder* bldr = GlueXDetectorConstruction::GetBuilder();
1371 std::map<std::string, std::vector<int> >::const_iterator iter;
1372 std::map<std::string, std::vector<int> > *identifiers;
1373 int max_depth = touch->GetHistoryDepth();
1374 for (int depth = 0; depth < max_depth; ++depth) {
1375 G4VPhysicalVolume *pvol = touch->GetVolume(depth);
1376 G4LogicalVolume *lvol = pvol->GetLogicalVolume();
1377 int volId = fVolumeTable[lvol];
1378 if (volId == 0) {
1379 volId = bldr->getVolumeId(lvol);
1380 fVolumeTable[lvol] = volId;
1381 }
1382 identifiers = &Refsys::fIdentifierTable[volId];
1383 if ((iter = identifiers->find(div)) != identifiers->end()) {
1384 int copyNum = touch->GetCopyNumber(depth);
1385 copyNum += (dynamic_cast<G4PVPlacement*>(pvol))? -1 : 0;
1386 return iter->second[copyNum];
1387 }
1388 }
1389 return -1;
1390}