1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | #include "GlueXSensitiveDetectorCDC.hh" |
8 | #include "GlueXDetectorConstruction.hh" |
9 | #include "GlueXPrimaryGeneratorAction.hh" |
10 | #include "GlueXUserEventInformation.hh" |
11 | #include "GlueXUserTrackInformation.hh" |
12 | #include "GlueXUserOptions.hh" |
13 | #include "HddmOutput.hh" |
14 | |
15 | #include <CLHEP/Random/RandPoisson.h> |
16 | #include <Randomize.hh> |
17 | |
18 | #include "G4VPhysicalVolume.hh" |
19 | #include "G4PVPlacement.hh" |
20 | #include "G4EventManager.hh" |
21 | #include "G4HCofThisEvent.hh" |
22 | #include "G4Step.hh" |
23 | #include "G4SDManager.hh" |
24 | #include "G4ios.hh" |
25 | |
26 | #include <JANA/JApplication.h> |
27 | |
28 | #include <stdlib.h> |
29 | |
30 | const double fC = 1e-15 * coulomb; |
31 | const double GlueXSensitiveDetectorCDC::ELECTRON_CHARGE = 1.6022e-4*fC; |
32 | |
33 | |
34 | double GlueXSensitiveDetectorCDC::DRIFT_SPEED = 0.0055*cm/ns; |
35 | |
36 | |
37 | double GlueXSensitiveDetectorCDC::TWO_HIT_TIME_RESOL = 25*ns; |
38 | |
39 | |
40 | int GlueXSensitiveDetectorCDC::MAX_HITS = 1000; |
41 | |
42 | |
43 | double GlueXSensitiveDetectorCDC::THRESH_KEV = 1.; |
44 | double GlueXSensitiveDetectorCDC::THRESH_MV = 1.; |
45 | |
46 | |
47 | double GlueXSensitiveDetectorCDC::STRAW_RADIUS = 0.776*cm; |
48 | |
49 | |
50 | double GlueXSensitiveDetectorCDC::CDC_TIME_WINDOW = 1000*ns; |
51 | |
52 | |
53 | double GlueXSensitiveDetectorCDC::GAS_GAIN = 1e5; |
54 | |
55 | |
56 | int GlueXSensitiveDetectorCDC::N_SECOND_PER_PRIMARY = 1.94; |
57 | |
58 | |
59 | double GlueXSensitiveDetectorCDC::W_EFF_PER_ION = 29.5*eV; |
60 | |
61 | |
62 | |
63 | |
64 | |
65 | |
66 | int GlueXSensitiveDetectorCDC::instanceCount = 0; |
67 | int GlueXSensitiveDetectorCDC::fDrift_clusters = 0; |
68 | double GlueXSensitiveDetectorCDC::fDrift_time[CDC_DRIFT_TABLE_LEN]; |
69 | double GlueXSensitiveDetectorCDC::fDrift_distance[CDC_DRIFT_TABLE_LEN]; |
70 | double GlueXSensitiveDetectorCDC::fBscale_par1; |
71 | double GlueXSensitiveDetectorCDC::fBscale_par2; |
72 | |
73 | G4Mutex GlueXSensitiveDetectorCDC::fMutex = G4MUTEX_INITIALIZER{ { 0, 0, 0, 0, 0, 0, 0, { 0, 0 } } }; |
74 | |
75 | GlueXSensitiveDetectorCDC::GlueXSensitiveDetectorCDC(const G4String& name) |
76 | : G4VSensitiveDetector(name), |
77 | fStrawsMap(0), fPointsMap(0) |
78 | { |
79 | collectionName.insert("CDCStrawHitsCollection"); |
80 | collectionName.insert("CDCPointsCollection"); |
81 | |
82 | |
83 | |
84 | |
85 | |
86 | |
87 | |
88 | G4AutoLock barrier(&fMutex); |
89 | if (instanceCount++ == 0) { |
90 | int runno = HddmOutput::getRunNo(); |
91 | extern jana::JApplication *japp; |
92 | if (japp == 0) { |
93 | G4cerr(*G4cerr_p) << "Error in GlueXSensitiveDetector constructor - " |
94 | << "jana global DApplication object not set, " |
95 | << "cannot continue." << G4endlstd::endl; |
96 | exit(-1); |
97 | } |
98 | jana::JCalibration *jcalib = japp->GetJCalibration(runno); |
99 | std::map<string, float> cdc_parms; |
100 | jcalib->Get("CDC/cdc_parms", cdc_parms); |
101 | DRIFT_SPEED = cdc_parms.at("CDC_DRIFT_SPEED")*cm/ns; |
102 | TWO_HIT_TIME_RESOL = cdc_parms.at("CDC_TWO_HIT_RESOL")*ns; |
103 | MAX_HITS = cdc_parms.at("CDC_MAX_HITS"); |
104 | THRESH_KEV = cdc_parms.at("CDC_THRESH_KEV"); |
105 | |
106 | |
107 | |
108 | |
109 | G4ThreeVector B = GlueXDetectorConstruction::GetInstance() |
110 | ->GetMagneticField(G4ThreeVector(0, 0, 65), tesla); |
111 | if (B.mag() > 1e-3) { |
112 | int nvalues = CDC_DRIFT_TABLE_LEN; |
113 | std::vector< std::map<std::string, float> > values; |
114 | jcalib->Get("CDC/cdc_drift_table", values); |
115 | for (int k=0; k < nvalues; ++k) { |
116 | fDrift_distance[k] = 0.01 * k; |
117 | fDrift_time[k] = values[k]["t"] * 1000; |
118 | } |
119 | |
120 | |
121 | |
122 | |
123 | std::vector< std::vector<float> > cdc_drift_parms; |
124 | jcalib->Get("/CDC/drift_parameters", cdc_drift_parms); |
125 | fBscale_par1 = cdc_drift_parms.at(0).at(9); |
126 | fBscale_par2 = cdc_drift_parms.at(0).at(10); |
127 | } |
128 | else { |
129 | int nvalues = CDC_DRIFT_TABLE_LEN; |
130 | std::vector< std::map<std::string, float> > values; |
131 | jcalib->Get("CDC/cdc_drift_table::NoBField", values); |
132 | for (int k=0; k < nvalues; ++k) { |
133 | fDrift_distance[k] = 0.01 * k; |
134 | fDrift_time[k] = values[k]["t"] * 1000; |
135 | } |
136 | fBscale_par1 = 1.; |
137 | fBscale_par2 = 0; |
138 | } |
139 | G4cout(*G4cout_p) << "CDC: ALL parameters loaded from ccdb" << G4endlstd::endl; |
140 | |
141 | |
142 | |
143 | GlueXUserOptions *opts = GlueXUserOptions::GetInstance(); |
144 | if (opts) { |
145 | std::map<int, int> driftclusters_opts; |
146 | if (opts->Find("driftclusters", driftclusters_opts)) |
147 | fDrift_clusters = driftclusters_opts[1]; |
148 | } |
149 | } |
150 | } |
151 | |
152 | GlueXSensitiveDetectorCDC::GlueXSensitiveDetectorCDC( |
153 | const GlueXSensitiveDetectorCDC &src) |
154 | : G4VSensitiveDetector(src), |
155 | fStrawsMap(src.fStrawsMap), fPointsMap(src.fPointsMap) |
156 | { |
157 | G4AutoLock barrier(&fMutex); |
158 | ++instanceCount; |
159 | } |
160 | |
161 | GlueXSensitiveDetectorCDC &GlueXSensitiveDetectorCDC::operator=(const |
162 | GlueXSensitiveDetectorCDC &src) |
163 | { |
164 | G4AutoLock barrier(&fMutex); |
165 | *(G4VSensitiveDetector*)this = src; |
166 | fStrawsMap = src.fStrawsMap; |
167 | fPointsMap = src.fPointsMap; |
168 | return *this; |
169 | } |
170 | |
171 | GlueXSensitiveDetectorCDC::~GlueXSensitiveDetectorCDC() |
172 | { |
173 | G4AutoLock barrier(&fMutex); |
174 | --instanceCount; |
175 | } |
176 | |
177 | void GlueXSensitiveDetectorCDC::Initialize(G4HCofThisEvent* hce) |
178 | { |
179 | fStrawsMap = new |
180 | GlueXHitsMapCDCstraw(SensitiveDetectorName, collectionName[0]); |
181 | fPointsMap = new |
182 | GlueXHitsMapCDCpoint(SensitiveDetectorName, collectionName[1]); |
183 | G4SDManager *sdm = G4SDManager::GetSDMpointer(); |
184 | hce->AddHitsCollection(sdm->GetCollectionID(collectionName[0]), fStrawsMap); |
185 | hce->AddHitsCollection(sdm->GetCollectionID(collectionName[1]), fPointsMap); |
186 | } |
187 | |
188 | G4bool GlueXSensitiveDetectorCDC::ProcessHits(G4Step* step, |
189 | G4TouchableHistory* ROhist) |
190 | { |
191 | double dEsum = step->GetTotalEnergyDeposit(); |
192 | if (dEsum == 0) |
193 | return false; |
194 | |
195 | const G4ThreeVector &pin = step->GetPreStepPoint()->GetMomentum(); |
196 | const G4ThreeVector &xin = step->GetPreStepPoint()->GetPosition(); |
197 | const G4ThreeVector &xout = step->GetPostStepPoint()->GetPosition(); |
198 | double tin = step->GetPreStepPoint()->GetGlobalTime(); |
199 | double tout = step->GetPostStepPoint()->GetGlobalTime(); |
200 | G4ThreeVector x = (xin + xout) / 2; |
201 | G4ThreeVector dx = xout - xin; |
202 | double t = (tin + tout) / 2; |
203 | |
204 | const G4VTouchable* touch = step->GetPreStepPoint()->GetTouchable(); |
205 | const G4AffineTransform &local_from_global = touch->GetHistory() |
206 | ->GetTopTransform(); |
207 | G4ThreeVector xinlocal = local_from_global.TransformPoint(xin); |
208 | G4ThreeVector xoutlocal = local_from_global.TransformPoint(xout); |
209 | |
210 | |
211 | |
212 | |
213 | |
214 | |
215 | |
216 | if (tout > 1.0*s) |
217 | t = tin; |
218 | |
219 | double drin = xinlocal.perp(); |
220 | double drout = xoutlocal.perp(); |
221 | |
222 | |
223 | G4ThreeVector xlocal = (xinlocal + xoutlocal) / 2; |
224 | G4ThreeVector dxlocal = xoutlocal - xinlocal; |
225 | double alpha = -(xinlocal[0] * dxlocal[0] + |
226 | xinlocal[1] * dxlocal[1]) / (dxlocal.perp2() + 1e-30); |
227 | G4ThreeVector xpoca = xinlocal + alpha * dxlocal; |
228 | double alpha01 = (alpha < 0)? 0 : (alpha > 1)? 1 : alpha; |
229 | G4ThreeVector x2local = xinlocal + alpha01 * dxlocal; |
230 | |
231 | |
232 | if (fabs(x2local[2]) >= 75.45*cm) { |
233 | int sign = (xoutlocal[2] > 0)? 1 : -1; |
234 | int ring = GetIdent("ring", touch); |
235 | if (ring <= 4 || (ring >= 13 && ring <= 16) || ring >= 25) { |
236 | alpha = (sign * 75.45*cm - xinlocal[2]) / (dxlocal[2] + 1e-30); |
237 | xpoca = xinlocal + alpha * dxlocal; |
238 | alpha01 = (alpha < 0)? 0 : (alpha > 1)? 1 : alpha; |
239 | x2local = xinlocal + alpha01 * dxlocal; |
240 | } |
241 | else if (fabs(x2local[2]) >= 75.575*cm) { |
242 | alpha = (sign * 75.575*cm - xinlocal[2]) / (dxlocal[2] + 1e-30); |
243 | xpoca = xinlocal + alpha * dxlocal; |
244 | alpha01 = (alpha < 0)? 0 : (alpha > 1)? 1 : alpha; |
245 | x2local = xinlocal + alpha01 * dxlocal; |
246 | } |
247 | } |
248 | |
249 | |
250 | |
251 | |
252 | |
253 | |
254 | |
255 | |
256 | |
257 | |
258 | |
259 | if (drin < 0.0050*cm) { |
260 | return false; |
261 | } |
262 | |
263 | if ((drin > (STRAW_RADIUS - 0.0200*cm) && drout < 0.0050*cm) || |
264 | (drin < 0.274*cm && drin > 0.234*cm && drout < 0.0050*cm)) |
265 | { |
266 | |
267 | |
268 | |
269 | |
270 | |
271 | |
272 | x = xout; |
273 | t = tout; |
274 | xlocal = xoutlocal; |
275 | |
276 | |
277 | |
278 | |
279 | |
280 | dx *= 2; |
281 | dEsum *= 2; |
282 | } |
283 | |
284 | double dradius = xpoca.perp(); |
285 | double dr = dx.mag(); |
286 | double dEdx = (dr > 1e-3*cm)? dEsum/dr : 0; |
287 | |
288 | |
289 | |
290 | |
291 | int ring = GetIdent("ring", touch); |
292 | int sector = GetIdent("sector", touch); |
293 | G4Track *track = step->GetTrack(); |
294 | G4int trackID = track->GetTrackID(); |
295 | GlueXUserTrackInformation *trackinfo = (GlueXUserTrackInformation*) |
296 | track->GetUserInformation(); |
297 | G4int itrack = trackinfo->GetGlueXTrackID(); |
298 | int pdgtype = track->GetDynamicParticle()->GetPDGcode(); |
299 | int g3type = GlueXPrimaryGeneratorAction::ConvertPdgToGeant3(pdgtype); |
300 | if (trackinfo->GetGlueXHistory() == 0 && itrack > 0) { |
301 | G4int key = fPointsMap->entries(); |
302 | GlueXHitCDCpoint* lastPoint = (*fPointsMap)[key - 1]; |
303 | |
304 | int ring = GetIdent("ring", touch); |
305 | if (lastPoint == 0 || lastPoint->track_ != trackID || |
306 | lastPoint->ring_ != ring) |
307 | { |
308 | GlueXHitCDCpoint newPoint; |
309 | newPoint.ptype_G3 = g3type; |
310 | newPoint.track_ = trackID; |
311 | newPoint.trackID_ = itrack; |
312 | newPoint.primary_ = (track->GetParentID() == 0); |
313 | newPoint.t_ns = t/ns; |
314 | newPoint.z_cm = x[2]/cm; |
315 | newPoint.r_cm = x.perp()/cm; |
316 | newPoint.phi_rad = x.phi(); |
317 | newPoint.dradius_cm = dradius/cm; |
318 | newPoint.px_GeV = pin[0]/GeV; |
319 | newPoint.py_GeV = pin[1]/GeV; |
320 | newPoint.pz_GeV = pin[2]/GeV; |
321 | newPoint.dEdx_GeV_cm = dEdx/(GeV/cm); |
322 | newPoint.sector_ = sector; |
323 | newPoint.ring_ = ring; |
324 | fPointsMap->add(key, newPoint); |
325 | } |
326 | } |
327 | |
328 | |
329 | |
330 | if (dEsum > 0) { |
331 | int key = GlueXHitCDCstraw::GetKey(ring, sector); |
332 | GlueXHitCDCstraw *straw = (*fStrawsMap)[key]; |
333 | if (straw == 0) { |
334 | GlueXHitCDCstraw newstraw(ring, sector); |
335 | fStrawsMap->add(key, newstraw); |
336 | straw = (*fStrawsMap)[key]; |
337 | } |
338 | |
339 | |
340 | |
341 | |
342 | int merge_hit = 0; |
343 | hit_vector_t::iterator hiter; |
344 | for (hiter = straw->hits.begin(); hiter != straw->hits.end(); ++hiter) { |
345 | if (itrack == hiter->itrack_ && fabs(tin - hiter->t1_ns*ns) < 0.01) { |
346 | merge_hit = 1; |
347 | break; |
348 | } |
349 | else if (hiter->t0_ns*ns > tin) { |
350 | break; |
351 | } |
352 | } |
353 | if (merge_hit) { |
354 | double d_cm = x2local.perp()/cm; |
355 | if (d_cm < hiter->d_cm) |
356 | hiter->d_cm = d_cm; |
357 | hiter->q_fC += dEsum; |
358 | hiter->t1_ns = tout/ns; |
359 | hiter->x1_g = xout; |
360 | } |
361 | else if ((int)straw->hits.size() < MAX_HITS) { |
362 | hiter = straw->hits.insert(hiter, GlueXHitCDCstraw::hitinfo_t()); |
363 | hiter->track_ = trackID; |
364 | hiter->q_fC = dEsum; |
365 | hiter->d_cm = x2local.perp()/cm; |
366 | hiter->itrack_ = itrack; |
367 | hiter->ptype_G3 = g3type; |
368 | hiter->t0_ns = tin/ns; |
369 | hiter->t1_ns = tout/ns; |
370 | hiter->x0_g = xin; |
371 | hiter->x1_g = xout; |
372 | } |
373 | else { |
374 | G4cerr(*G4cerr_p) << "GlueXSensitiveDetectorCDC::ProcessHits error: " |
375 | << "max hit count " << MAX_HITS << " exceeded, truncating!" |
376 | << G4endlstd::endl; |
377 | } |
378 | } |
379 | return true; |
380 | } |
381 | |
382 | void GlueXSensitiveDetectorCDC::EndOfEvent(G4HCofThisEvent*) |
383 | { |
384 | std::map<int, GlueXHitCDCstraw*> *straws = fStrawsMap->GetMap(); |
385 | std::map<int, GlueXHitCDCpoint*> *points = fPointsMap->GetMap(); |
386 | if (straws->size() == 0 && points->size() == 0) |
387 | return; |
388 | std::map<int, GlueXHitCDCstraw*>::iterator siter; |
389 | std::map<int, GlueXHitCDCpoint*>::iterator piter; |
390 | |
391 | if (verboseLevel > 1) { |
392 | G4cout(*G4cout_p) << G4endlstd::endl |
393 | << "--------> Hits Collection: in this event there are " |
394 | << straws->size() << " straws with hits in the CDC: " |
395 | << G4endlstd::endl; |
396 | for (siter = straws->begin(); siter != straws->end(); ++siter) |
397 | siter->second->Print(); |
398 | |
399 | G4cout(*G4cout_p) << G4endlstd::endl |
400 | << "--------> Hits Collection: in this event there are " |
401 | << points->size() << " truth points in the CDC: " |
402 | << G4endlstd::endl; |
403 | for (piter = points->begin(); piter != points->end(); ++piter) |
404 | piter->second->Print(); |
405 | } |
406 | |
407 | |
408 | |
409 | G4EventManager* mgr = G4EventManager::GetEventManager(); |
410 | G4VUserEventInformation* info = mgr->GetUserInformation(); |
411 | hddm_s::HDDM *record = ((GlueXUserEventInformation*)info)->getOutputRecord(); |
412 | if (record == 0) { |
413 | G4cerr(*G4cerr_p) << "GlueXSensitiveDetectorCDC::EndOfEvent error - " |
414 | << "hits seen but no output hddm record to save them into, " |
415 | << "cannot continue!" << G4endlstd::endl; |
416 | exit(1); |
417 | } |
418 | |
419 | if (record->getPhysicsEvents().size() == 0) |
420 | record->addPhysicsEvents(); |
421 | if (record->getHitViews().size() == 0) |
422 | record->getPhysicsEvent().addHitViews(); |
423 | hddm_s::HitView &hitview = record->getPhysicsEvent().getHitView(); |
424 | if (hitview.getCentralDCs().size() == 0) |
425 | hitview.addCentralDCs(); |
426 | hddm_s::CentralDC ¢ralDC = hitview.getCentralDC(); |
427 | |
428 | |
429 | |
430 | for (siter = straws->begin(); siter != straws->end(); ++siter) { |
431 | |
432 | |
433 | |
434 | |
435 | hit_vector_t &splits = siter->second->hits; |
436 | hit_vector_t hits; |
437 | while (splits.size() > 0) { |
438 | for (unsigned int ih=1; ih < splits.size(); ++ih) { |
439 | if (fabs(splits[ih].t0_ns - splits[0].t1_ns) < 0.5*ns) { |
440 | splits[0].t1_ns = splits[ih].t1_ns; |
441 | splits[0].q_fC += splits[ih].q_fC; |
442 | if (splits[ih].d_cm < splits[0].d_cm) |
443 | splits[0].d_cm = splits[ih].d_cm; |
444 | splits[0].x1_g = splits[ih].x1_g; |
445 | splits.erase(splits.begin() + ih); |
446 | --ih; |
447 | } |
448 | } |
449 | if (splits[0].q_fC/keV > THRESH_KEV) { |
450 | |
451 | |
452 | |
453 | |
454 | |
455 | |
456 | double t = (splits[0].t0_ns + splits[0].t1_ns)*ns / 2; |
457 | double n_p_mean = (splits[0].q_fC / W_EFF_PER_ION) / |
458 | (1 + N_SECOND_PER_PRIMARY); |
459 | |
460 | int n_p = CLHEP::RandPoisson::shoot(n_p_mean); |
461 | if (fDrift_clusters == 0) { |
462 | add_cluster(hits, splits[0], n_p, t, splits[0].x0_g); |
463 | } |
464 | else { |
465 | |
466 | |
467 | |
468 | G4ThreeVector dx = splits[0].x1_g - splits[0].x0_g; |
469 | for (int n=0; n < n_p; n++) { |
470 | double u = G4RandFlatG4MTRandFlat::shoot(); |
471 | G4ThreeVector x = splits[0].x0_g + u * dx; |
472 | add_cluster(hits, splits[0], 1, t, x); |
473 | } |
474 | } |
475 | } |
476 | splits.erase(splits.begin()); |
477 | } |
478 | |
479 | |
480 | |
481 | |
482 | if (fDrift_clusters) { |
483 | |
484 | int num_samples = int(CDC_TIME_WINDOW/ns); |
485 | double *samples = new double[num_samples]; |
486 | for (int i=0; i < num_samples; i++) { |
487 | samples[i] = cdc_wire_signal_mV(double(i), hits); |
488 | } |
489 | |
490 | |
491 | double dradius_cm = hits[0].d_cm; |
492 | int ptype_G3 = hits[0].ptype_G3; |
493 | int track = hits[0].track_; |
494 | int itrack = hits[0].itrack_; |
495 | |
496 | hits.clear(); |
497 | double q_mV_ns = 0.; |
498 | int over_threshold = 0; |
499 | for (int i=0; i < num_samples; ++i) { |
500 | if (samples[i] > THRESH_MV) { |
501 | if (!over_threshold) { |
502 | hits.push_back(GlueXHitCDCstraw::hitinfo_t()); |
503 | hits.back().track_ = track; |
504 | hits.back().t_ns = double(i); |
505 | hits.back().d_cm = dradius_cm; |
506 | hits.back().itrack_ = itrack; |
507 | hits.back().ptype_G3 = ptype_G3; |
508 | over_threshold = 1; |
509 | } |
510 | q_mV_ns += samples[i]; |
511 | } |
512 | else if (over_threshold) { |
513 | hits.back().q_fC = q_mV_ns; |
514 | over_threshold = 0; |
515 | q_mV_ns = 0; |
516 | } |
517 | } |
518 | delete [] samples; |
519 | } |
520 | else { |
521 | |
522 | |
523 | |
524 | hit_vector_t::iterator hiter; |
525 | for (hiter = hits.begin(); hiter != hits.end(); ++hiter) { |
526 | int merge_hits = 0; |
527 | hit_vector_t::iterator iter; |
528 | for (iter = splits.begin(); iter != splits.end(); ++iter) { |
529 | if (fabs(hiter->t_ns - iter->t_ns) < TWO_HIT_TIME_RESOL) { |
530 | merge_hits = 1; |
531 | break; |
532 | } |
533 | else if (hiter->t_ns < iter->t_ns) { |
534 | break; |
535 | } |
536 | } |
537 | if (merge_hits) { |
538 | if (hiter->t_ns < iter->t_ns) { |
539 | double q_fC = iter->q_fC; |
540 | *iter = *hiter; |
541 | iter->q_fC += q_fC; |
542 | } |
543 | else { |
544 | iter->q_fC += hiter->q_fC; |
545 | } |
546 | } |
547 | else { |
548 | splits.insert(iter, *hiter); |
549 | } |
550 | } |
551 | } |
552 | |
553 | if (hits.size() > 0) { |
554 | hddm_s::CdcStrawList straw = centralDC.addCdcStraws(1); |
555 | straw(0).setRing(siter->second->ring_); |
556 | straw(0).setStraw(siter->second->sector_); |
557 | for (int ih=0; ih < (int)hits.size(); ++ih) { |
558 | hddm_s::CdcStrawTruthHitList thit = straw(0).addCdcStrawTruthHits(1); |
559 | thit(0).setQ(hits[ih].q_fC); |
560 | thit(0).setT(hits[ih].t_ns); |
561 | thit(0).setD(hits[ih].d_cm); |
562 | thit(0).setItrack(hits[ih].itrack_); |
563 | thit(0).setPtype(hits[ih].ptype_G3); |
564 | } |
565 | } |
566 | } |
567 | |
568 | |
569 | for (piter = points->begin(); piter != points->end(); ++piter) { |
570 | hddm_s::CdcTruthPointList point = centralDC.addCdcTruthPoints(1); |
571 | point(0).setDEdx(piter->second->dEdx_GeV_cm); |
572 | point(0).setDradius(piter->second->dradius_cm); |
573 | point(0).setPrimary(piter->second->primary_); |
574 | point(0).setPtype(piter->second->ptype_G3); |
575 | point(0).setPx(piter->second->px_GeV); |
576 | point(0).setPy(piter->second->py_GeV); |
577 | point(0).setPz(piter->second->pz_GeV); |
578 | point(0).setR(piter->second->r_cm); |
579 | point(0).setPhi(piter->second->phi_rad); |
580 | point(0).setZ(piter->second->z_cm); |
581 | point(0).setT(piter->second->t_ns); |
582 | point(0).setTrack(piter->second->track_); |
583 | hddm_s::TrackIDList tid = point(0).addTrackIDs(); |
584 | tid(0).setItrack(piter->second->trackID_); |
585 | } |
586 | } |
587 | |
588 | double GlueXSensitiveDetectorCDC::asic_response(double t_ns) |
589 | { |
590 | |
591 | |
592 | double par[11] = {-0.01986, 0.01802, -0.001097, 10.3, 11.72, |
593 | -0.03701, 35.84, 15.93, 0.006141, 80.95, 24.77}; |
594 | if (t_ns < par[3]) |
595 | return par[0] * t_ns + par[1] * t_ns*t_ns + par[2] * t_ns*t_ns*t_ns; |
596 | else |
597 | return (par[0] * par[3] + |
598 | par[1] * par[3]*par[3] + |
599 | par[2] * par[3]*par[3]*par[3]) * |
600 | exp(-pow((t_ns - par[3])/par[4], 2)) + |
601 | par[5] * exp(-pow((t_ns - par[6])/par[7], 2)) + |
602 | par[8] * exp(-pow((t_ns - par[9])/par[10], 2)); |
603 | } |
604 | |
605 | double GlueXSensitiveDetectorCDC::cdc_wire_signal_mV(double t_ns, |
606 | hit_vector_t &hits) |
607 | { |
608 | |
609 | |
610 | double asic_gain = 0.5; |
611 | double signal_mV = 0; |
612 | hit_vector_t::iterator hiter; |
613 | for (hiter = hits.begin(); hiter != hits.end(); ++hiter) { |
614 | if (t_ns > hiter->t_ns) { |
615 | double my_time_ns = t_ns - hiter->t_ns; |
616 | signal_mV += asic_gain * hiter->q_fC * asic_response(my_time_ns); |
617 | } |
618 | } |
619 | return signal_mV; |
620 | } |
621 | |
622 | void GlueXSensitiveDetectorCDC::add_cluster(hit_vector_t &hits, |
623 | GlueXHitCDCstraw::hitinfo_t &hit, |
624 | int n_p, |
625 | double t, |
626 | G4ThreeVector &x) |
627 | { |
628 | |
629 | double q_fC = 0; |
630 | |
631 | |
632 | double dradius_cm = hit.d_cm; |
633 | double d2 = dradius_cm * dradius_cm; |
634 | double d3 = dradius_cm * d2; |
| Value stored to 'd3' during its initialization is never read |
635 | |
636 | |
637 | |
638 | |
639 | G4ThreeVector B = GlueXDetectorConstruction::GetInstance() |
640 | ->GetMagneticField(x, tesla); |
641 | double BmagT = B.mag(); |
642 | |
643 | |
644 | |
645 | double my_t_ns; |
646 | double my_t_err; |
647 | int i = (int)(dradius_cm * 100); |
648 | if (i >= CDC_DRIFT_TABLE_LEN - 3) { |
649 | |
650 | my_t_ns = fDrift_time[CDC_DRIFT_TABLE_LEN - 3] + |
651 | (dradius_cm - fDrift_distance[CDC_DRIFT_TABLE_LEN - 3]) * |
652 | (fDrift_time[CDC_DRIFT_TABLE_LEN - 1] - |
653 | fDrift_time[CDC_DRIFT_TABLE_LEN - 3]) / 0.02; |
654 | } |
655 | else { |
656 | int index = (i < 1)? 0 : i - 1; |
657 | |
658 | |
659 | |
660 | polint(&fDrift_distance[index], |
661 | &fDrift_time[index], 4, dradius_cm, &my_t_ns, &my_t_err); |
662 | } |
663 | double tdrift_ns = my_t_ns / (fBscale_par1 + fBscale_par2 * BmagT); |
664 | |
665 | |
666 | |
667 | double v_max = 0.08; |
668 | double tmin_ns = dradius_cm / v_max; |
669 | if (tdrift_ns < tmin_ns) { |
670 | tdrift_ns = tmin_ns; |
671 | } |
672 | double total_time = t + tdrift_ns*ns; |
673 | |
674 | |
675 | if (total_time > CDC_TIME_WINDOW) { |
676 | hit.q_fC = 0; |
677 | hit.t_ns = 0; |
678 | return; |
679 | } |
680 | |
681 | if (fDrift_clusters == 0) { |
682 | |
683 | |
684 | |
685 | |
686 | double n_s_mean = n_p * N_SECOND_PER_PRIMARY; |
687 | int n_s = CLHEP::RandPoisson::shoot(n_s_mean); |
688 | int n_t = n_s + n_p; |
689 | q_fC = n_t * GAS_GAIN * ELECTRON_CHARGE/fC; |
690 | } |
691 | else { |
692 | |
693 | |
694 | |
695 | |
696 | int n_s = CLHEP::RandPoisson::shoot(N_SECOND_PER_PRIMARY); |
697 | |
698 | q_fC = GAS_GAIN * ELECTRON_CHARGE/fC * (1 + n_s); |
699 | } |
700 | hit.q_fC = q_fC; |
701 | hit.t_ns = total_time/ns; |
702 | hits.push_back(hit); |
703 | } |
704 | |
705 | void GlueXSensitiveDetectorCDC::polint(double *xa, double *ya, int n, |
706 | double x, double *y, double *dy) |
707 | { |
708 | |
709 | |
710 | |
711 | |
712 | |
713 | |
714 | double *c = NULL__null; |
715 | double *d = NULL__null; |
716 | double den; |
717 | double dif; |
718 | double dift; |
719 | double ho; |
720 | double hp; |
721 | double w; |
722 | |
723 | int i; |
724 | int m; |
725 | int ns; |
726 | |
727 | if ((c = (double*)malloc(n*sizeof(double))) == NULL__null || |
728 | (d = (double*)malloc(n*sizeof(double))) == NULL__null ) |
729 | { |
730 | fprintf(stderrstderr, "polint error: allocating workspace\n" ); |
731 | fprintf(stderrstderr, "polint error: setting y = 0 and dy = 1e9\n" ); |
732 | *y = 0.0; |
733 | *dy = 1.e9; |
734 | if (c != NULL__null) |
735 | free(c); |
736 | if (d != NULL__null ) |
737 | free(d); |
738 | return; |
739 | } |
740 | |
741 | ns = 0; |
742 | dif = fabs(x-xa[0]); |
743 | for (i = 0; i < n; ++i) { |
744 | dift = fabs(x-xa[i]); |
745 | if (dift < dif) { |
746 | ns = i; |
747 | dif = dift; |
748 | } |
749 | c[i] = ya[i]; |
750 | d[i] = ya[i]; |
751 | } |
752 | *y = ya[ns]; |
753 | ns = ns-1; |
754 | for (m = 0; m < n-1; ++m) { |
755 | for (i = 0; i < n-m-1; ++i) { |
756 | ho = xa[i]-x; |
757 | hp = xa[i+m+1]-x; |
758 | w = c[i+1]-d[i]; |
759 | den = ho-hp; |
760 | if (den == 0) { |
761 | fprintf( stderrstderr, "polint error: den = 0\n" ); |
762 | fprintf( stderrstderr, "polint error: setting y = 0 and dy = 1e9\n" ); |
763 | *y = 0.0; |
764 | *dy = 1.e9; |
765 | if (c != NULL__null) |
766 | free(c); |
767 | if (d != NULL__null) |
768 | free(d); |
769 | return; |
770 | } |
771 | den = w/den; |
772 | d[i] = hp*den; |
773 | c[i] = ho*den; |
774 | } |
775 | if (2*(ns+1) < n-m-1) { |
776 | *dy = c[ns+1]; |
777 | } |
778 | else { |
779 | *dy = d[ns]; |
780 | ns = ns-1; |
781 | } |
782 | *y = (*y)+(*dy); |
783 | } |
784 | |
785 | if (c != NULL__null) |
786 | free(c); |
787 | if (d != NULL__null) |
788 | free(d); |
789 | } |
790 | |
791 | int GlueXSensitiveDetectorCDC::GetIdent(std::string div, |
792 | const G4VTouchable *touch) |
793 | { |
794 | const HddsG4Builder* bldr = GlueXDetectorConstruction::GetBuilder(); |
795 | std::map<std::string, std::vector<int> >::const_iterator iter; |
796 | std::map<std::string, std::vector<int> > *identifiers; |
797 | int max_depth = touch->GetHistoryDepth(); |
798 | for (int depth = 0; depth < max_depth; ++depth) { |
799 | G4VPhysicalVolume *pvol = touch->GetVolume(depth); |
800 | G4LogicalVolume *lvol = pvol->GetLogicalVolume(); |
801 | int volId = fVolumeTable[lvol]; |
802 | if (volId == 0) { |
803 | volId = bldr->getVolumeId(lvol); |
804 | fVolumeTable[lvol] = volId; |
805 | } |
806 | identifiers = &Refsys::fIdentifierTable[volId]; |
807 | if ((iter = identifiers->find(div)) != identifiers->end()) { |
808 | int copyNum = touch->GetCopyNumber(depth); |
809 | copyNum += (dynamic_cast<G4PVPlacement*>(pvol))? -1 : 0; |
810 | return iter->second[copyNum]; |
811 | } |
812 | } |
813 | return -1; |
814 | } |