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