1 | |
2 | |
3 | |
4 | |
5 | |
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 | |
31 | const double fC = 1e-15 * coulomb; |
32 | const double pC = 1e-12 * coulomb; |
33 | const double GlueXSensitiveDetectorFDC::ELECTRON_CHARGE = 1.6022e-4*fC; |
34 | |
35 | |
36 | double GlueXSensitiveDetectorFDC::DRIFT_SPEED = 0.0055*cm/ns; |
37 | |
38 | |
39 | double GlueXSensitiveDetectorFDC::TWO_HIT_TIME_RESOL = 25*ns; |
40 | |
41 | |
42 | int GlueXSensitiveDetectorFDC::MAX_HITS = 1000; |
43 | |
44 | |
45 | double GlueXSensitiveDetectorFDC::THRESH_KEV = 1.; |
46 | double GlueXSensitiveDetectorFDC::THRESH_ANODE = 1.*pC; |
47 | double GlueXSensitiveDetectorFDC::THRESH_STRIPS = 5.*pC; |
48 | |
49 | |
50 | double GlueXSensitiveDetectorFDC::FDC_TIME_WINDOW = 1000*ns; |
51 | |
52 | |
53 | double GlueXSensitiveDetectorFDC::GAS_GAIN = 8e4; |
54 | |
55 | |
56 | double GlueXSensitiveDetectorFDC::N_SECOND_PER_PRIMARY = 1.89; |
57 | |
58 | |
59 | double GlueXSensitiveDetectorFDC::W_EFF_PER_ION = 30.2*eV; |
60 | |
61 | |
62 | int GlueXSensitiveDetectorFDC::WIRES_PER_PLANE = 96; |
63 | int GlueXSensitiveDetectorFDC::STRIPS_PER_PLANE = 192; |
64 | double GlueXSensitiveDetectorFDC::ACTIVE_AREA_OUTER_RADIUS = 48.5*cm; |
65 | double GlueXSensitiveDetectorFDC::ANODE_CATHODE_SPACING = 0.5*cm; |
66 | double GlueXSensitiveDetectorFDC::WIRE_SPACING = 1.0*cm; |
67 | double GlueXSensitiveDetectorFDC::STRIP_SPACING = 0.5*cm; |
68 | double GlueXSensitiveDetectorFDC::U_OF_WIRE_ONE = -47.5*cm; |
69 | double GlueXSensitiveDetectorFDC::U_OF_STRIP_ONE = -47.75*cm; |
70 | double GlueXSensitiveDetectorFDC::CATHODE_ROT_ANGLE = 1.309; |
71 | double GlueXSensitiveDetectorFDC::STRIP_GAP = 0.1*cm; |
72 | double GlueXSensitiveDetectorFDC::STRIP_NODES = 3; |
73 | |
74 | |
75 | double GlueXSensitiveDetectorFDC::LORENTZ_NR_PAR1; |
76 | double GlueXSensitiveDetectorFDC::LORENTZ_NR_PAR2; |
77 | double GlueXSensitiveDetectorFDC::LORENTZ_NZ_PAR1; |
78 | double GlueXSensitiveDetectorFDC::LORENTZ_NZ_PAR2; |
79 | double GlueXSensitiveDetectorFDC::DRIFT_RES_PARMS[3]; |
80 | double GlueXSensitiveDetectorFDC::DRIFT_FUNC_PARMS[6]; |
81 | double GlueXSensitiveDetectorFDC::DRIFT_BSCALE_PAR1; |
82 | double GlueXSensitiveDetectorFDC::DRIFT_BSCALE_PAR2; |
83 | |
84 | |
85 | double GlueXSensitiveDetectorFDC::DIFFUSION_COEFF = 1.1e-6*cm*cm/s; |
86 | double GlueXSensitiveDetectorFDC::K2 = 1.15; |
87 | |
88 | |
89 | double GlueXSensitiveDetectorFDC::wire_dead_zone_radius[4] = |
90 | {3.0*cm, 3.0*cm, 3.9*cm, 3.9*cm}; |
91 | double GlueXSensitiveDetectorFDC::strip_dead_zone_radius[4] = |
92 | {1.3*cm, 1.3*cm, 1.3*cm, 1.3*cm}; |
93 | |
94 | |
95 | int GlueXSensitiveDetectorFDC::drift_table_len; |
96 | double *GlueXSensitiveDetectorFDC::drift_table_t_ns; |
97 | double *GlueXSensitiveDetectorFDC::drift_table_d_cm; |
98 | |
99 | int GlueXSensitiveDetectorFDC::instanceCount = 0; |
100 | G4Mutex GlueXSensitiveDetectorFDC::fMutex = G4MUTEX_INITIALIZER{ { 0, 0, 0, 0, 0, 0, 0, { 0, 0 } } }; |
101 | int GlueXSensitiveDetectorFDC::fDrift_clusters = 0; |
102 | |
103 | GlueXSensitiveDetectorFDC::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 | |
112 | |
113 | |
114 | |
115 | |
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 | |
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 | |
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 | |
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 | |
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 | |
185 | |
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; |
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 | |
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 | |
225 | GlueXSensitiveDetectorFDC::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 | |
236 | GlueXSensitiveDetectorFDC &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 | |
247 | GlueXSensitiveDetectorFDC::~GlueXSensitiveDetectorFDC() |
248 | { |
249 | G4AutoLock barrier(&fMutex); |
250 | --instanceCount; |
251 | } |
252 | |
253 | void 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 | |
267 | G4bool 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 | |
297 | |
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 | |
305 | |
306 | |
307 | |
308 | |
309 | |
310 | |
311 | |
312 | |
313 | |
314 | |
315 | |
316 | |
317 | |
318 | |
319 | |
320 | G4int rkey = 1 << 24; |
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]; |
332 | segment.y_cm = xin[1]; |
333 | segment.z_cm = xin[2]; |
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; |
| Value stored to 't' during its initialization is never read |
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 | |
391 | |
392 | |
393 | |
394 | |
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 | |
405 | |
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 | |
428 | |
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 | |
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 | |
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 | |
468 | |
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 | |
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 | |
524 | |
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 | |
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 | |
566 | void 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 | |
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 | |
622 | |
623 | for (witer = wires->begin(); witer != wires->end(); ++witer) { |
624 | |
625 | |
626 | |
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 | |
650 | |
651 | |
652 | |
653 | |
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 | |
659 | double n_p_mean = dE / W_EFF_PER_ION / (1 + N_SECOND_PER_PRIMARY); |
660 | |
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 | |
688 | for (int n=0; n < n_p; n++) { |
689 | |
690 | |
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 | |
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 | |
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 | |
745 | |
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 | |
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 | |
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 | |
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 | |
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 | |
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 | |
903 | if (fabs(hiter->t_ns - t_ns) < TWO_HIT_TIME_RESOL || |
904 | hiter->q_fC == 0) |
905 | { |
906 | |
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 | |
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 | |
1003 | double GlueXSensitiveDetectorFDC::asic_response(double t_ns) |
1004 | { |
1005 | |
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 | |
1020 | double GlueXSensitiveDetectorFDC::fdc_wire_signal_mV( |
1021 | double t_ns, std::vector<GlueXHitFDCwire::hitinfo_t> &hits) |
1022 | { |
1023 | |
1024 | |
1025 | double asic_gain = 0.76; |
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 | |
1037 | double GlueXSensitiveDetectorFDC::fdc_cathode_signal_mV( |
1038 | double t_ns, std::vector<GlueXHitFDCcathode::hitinfo_t> &hits) |
1039 | { |
1040 | |
1041 | |
1042 | double asic_gain = 2.3; |
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 | |
1054 | void 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 | |
1069 | |
1070 | |
1071 | |
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 | |
1078 | |
1079 | |
1080 | |
1081 | int n_s = CLHEP::RandPoisson::shoot(N_SECOND_PER_PRIMARY); |
1082 | q_anode = (1 + n_s) * GAS_GAIN * ELECTRON_CHARGE; |
1083 | } |
1084 | |
1085 | |
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 | |
1094 | |
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 | |
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 | |
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 | } |
1136 | } |
1137 | } |
1138 | |
1139 | int 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 | |
1151 | G4ThreeVector B = GlueXDetectorConstruction::GetInstance() |
1152 | ->GetMagneticField(xglobal, tesla); |
1153 | double BrhoT = B.perp(); |
1154 | |
1155 | |
1156 | |
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 | |
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 | |
1172 | |
1173 | |
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 | |
1180 | xlocal[1] += G4RandGaussG4MTRandGaussQ::shoot() * |
1181 | (0.01*cm * pow((dx2 + dz2)/cm2, 0.125) + 0.0061*cm * dx2/cm2); |
1182 | |
1183 | |
1184 | |
1185 | if (sqrt(xlocal[1] * xlocal[1] + xwire * xwire) > ACTIVE_AREA_OUTER_RADIUS) |
1186 | return 0; |
1187 | |
1188 | |
1189 | |
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; |
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 | |
1213 | tdrift_unsmeared *= 1. + DRIFT_BSCALE_PAR1 + DRIFT_BSCALE_PAR2*B[2]*B[2]; |
1214 | |
1215 | |
1216 | double v_max = 0.08*cm/ns; |
1217 | double tmin = dradius / v_max; |
1218 | |
1219 | |
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 | |
1231 | tdrift = t + tdrift_smeared; |
1232 | |
1233 | |
1234 | if (tdrift > FDC_TIME_WINDOW) |
1235 | return 0; |
1236 | |
1237 | |
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 | |
1252 | void GlueXSensitiveDetectorFDC::polint(double *xa, double *ya, int n, |
1253 | double x, double *y, double *dy) |
1254 | { |
1255 | |
1256 | |
1257 | |
1258 | |
1259 | |
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 | |
1338 | int GlueXSensitiveDetectorFDC::locate(double *xx, int n, double x) |
1339 | { |
1340 | |
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 | |
1367 | int 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 | } |