BOSS 7.0.1
BESIII Offline Software System
Loading...
Searching...
No Matches
BesTofDigitizerEcV4.cc
Go to the documentation of this file.
1//---------------------------------------------------------------------------//
2// BOOST --- BESIII Object_Oriented Simulation Tool //
3//---------------------------------------------------------------------------//
4//Description: This is the Digitizer for the MRPC with doubled sided readout
5//Author: An Fenfen
6//Created: Nov, 2015
7
8//---------------------------------------------------------------------------//
9//$Id: BesTofDigitizerEcV4.cc
10
11#include "GaudiKernel/MsgStream.h"
12#include "GaudiKernel/Bootstrap.h"
13#include "GaudiKernel/PropertyMgr.h"
14#include "GaudiKernel/IJobOptionsSvc.h"
15#include "GaudiKernel/ISvcLocator.h"
16#include "GaudiKernel/IDataProviderSvc.h"
17
18#include "BesTofDigitizerEcV4.hh"
19#include "BesTofDigi.hh"
20#include "BesTofHit.hh"
21#include "G4DigiManager.hh"
22#include "Randomize.hh"
23#include "TMath.h"
24#include <math.h>
25#include "TNtuple.h"
26#include "TFile.h"
27#include "TH1F.h"
28#include "TROOT.h"
29#include "TSpline.h"
30
31
32BesTofDigitizerEcV4::BesTofDigitizerEcV4()
33{
34 PropertyMgr m_propMgr;
35 m_propMgr.declareProperty("FileName", m_fileName = "mrpc.root");
36 m_propMgr.declareProperty("RootFlag", m_rootFlag = false);
37 m_propMgr.declareProperty("E", m_V = 7000);
38 m_propMgr.declareProperty("Threshold", m_threshold = 6e+08);
39
40 m_propMgr.declareProperty("nstep", m_nstep = -1);
41 m_propMgr.declareProperty("E_weight", m_E_weight = -1);
42 m_propMgr.declareProperty("saturationFlag", m_saturationFlag = true);
43 m_propMgr.declareProperty("tdcRes_const", tdcRes_const = -1);
44 m_propMgr.declareProperty("adcRes_const", adcRes_const = -1);
45 m_propMgr.declareProperty("calTdcRes_charge_flag", m_calTdcRes_charge_flag=0);
46 m_propMgr.declareProperty("charge2Time_flag", m_charge2Time_flag=0);
47 m_propMgr.declareProperty("calAdcRes_charge_flag", m_calAdcRes_charge_flag=0);
48
49
50 IJobOptionsSvc* jobSvc;
51 Gaudi::svcLocator()->service("JobOptionsSvc", jobSvc);
52 jobSvc->setMyProperties("BesTofDigitizerEcV4", &m_propMgr);
53
54 initial();
55
56 if(m_rootFlag)
57 {
58 m_file = new TFile(m_fileName.c_str(), "RECREATE");
59 m_tree = new TTree("mrpc", "mrpc");
60
61 m_tree->Branch("partId", &m_partId, "partId/D");
62 m_tree->Branch("module", &m_module, "module/D");
63 m_tree->Branch("time_leading_sphi", &m_time_leading_sphi, "time_leading_sphi/D");
64 m_tree->Branch("time_leading_xphi", &m_time_leading_xphi, "time_leading_xphi/D");
65 m_tree->Branch("time_trailing_sphi", &m_time_trailing_sphi, "time_trailing_sphi/D");
66 m_tree->Branch("time_trailing_xphi", &m_time_trailing_xphi, "time_trailing_xphi/D");
67 m_tree->Branch("tdcRes", &m_tdcRes, "tdcRes/D");
68 m_tree->Branch("tdcRes_charge", &m_tdcRes_charge, "tdcRes_charge/D");
69 m_tree->Branch("adc", &m_adc, "adc/D");
70 m_tree->Branch("adcRes", &m_adcRes, "adcRes/D");
71 m_tree->Branch("adcRes_charge", &m_adcRes_charge, "adcRes_charge/D");
72 m_tree->Branch("strip", &m_strip, "strip/D");
73 m_tree->Branch("trkIndex", &m_trkIndex, "trkIndex/D");
74 m_tree->Branch("tStart", &m_tStart, "tStart/D");
75 m_tree->Branch("tPropagate_sphi", &m_tPropagate_sphi, "tPropagate_sphi/D");
76 m_tree->Branch("tPropagate_xphi", &m_tPropagate_xphi, "tPropagate_xphi/D");
77 m_tree->Branch("tThreshold", &m_tThreshold, "tThreshold/D");
78 m_tree->Branch("charge", &m_charge, "charge/D");
79 m_tree->Branch("nhit", &m_nhit, "nhit/I");
80 m_tree->Branch("ions_hit", m_ions_hit, "ions_hit[nhit]/D");
81 m_tree->Branch("trkIndex_hit", m_trkIndex_hit, "trkIndex_hit[nhit]/D");
82 m_tree->Branch("pdgCode_hit", m_pdgCode_hit, "pdgCode_hit[nhit]/D");
83 m_tree->Branch("gap_hit", m_gap_hit, "gap_hit[nhit]/D");
84 m_tree->Branch("underStrip_hit", m_underStrip_hit, "underStrip_hit[nhit]/D");
85 m_tree->Branch("locx_hit", m_locx_hit, "locx_hit[nhit]/D");
86 m_tree->Branch("locy_hit", m_locy_hit, "locy_hit[nhit]/D");
87 m_tree->Branch("locz_hit", m_locz_hit, "locz_hit[nhit]/D");
88 m_tree->Branch("x_hit", m_x_hit, "x_hit[nhit]/D");
89 m_tree->Branch("y_hit", m_y_hit, "y_hit[nhit]/D");
90 m_tree->Branch("z_hit", m_z_hit, "z_hit[nhit]/D");
91 m_tree->Branch("px_hit", m_px_hit, "px_hit[nhit]/D");
92 m_tree->Branch("py_hit", m_py_hit, "py_hit[nhit]/D");
93 m_tree->Branch("pz_hit", m_pz_hit, "pz_hit[nhit]/D");
94 }
95}
96
97
98BesTofDigitizerEcV4::~BesTofDigitizerEcV4()
99{
100 if(m_rootFlag)
101 {
102 m_file->Write();
103 m_file->Close();
104 }
105}
106
107void BesTofDigitizerEcV4::initial()
108{
109 m_param = Param();
110 m_param.setPar(m_nstep, m_E_weight);
111 partId = -999;
112 module = -999;
113 tdc_sphi = -999;
114 tdc_xphi = -999;
115 if(tdcRes_const<0) tdcRes_const = sqrt(27*27.+10.*10+20.*20)+18.; //ps, 27:TDC; 10:gapNb; 20:cables..
116 tdcRes = -999;
117
118 adc = -999;
119 if(adcRes_const<0) adcRes_const = 27; //ps TDC
120 adcRes = -999;
121
122 time_leading_sphi = -999;
123 time_leading_xphi = -999;
124 time_trailing_sphi = -999;
125 time_trailing_xphi = -999;
126}
127
128void BesTofDigitizerEcV4::Digitize( ScintSingle* single_module, BesTofDigitsCollection* DC )
129{
130 m_besTofDigitsCollection = DC;
131 G4DigiManager* digiManager = G4DigiManager::GetDMpointer();
132 G4int THCID = digiManager->GetHitsCollectionID("BesTofHitsCollection");
133 m_THC = (BesTofHitsCollection*) (digiManager->GetHitsCollection(THCID));
134 if( !m_THC ) return;
135
136 partId = single_module->GetPartId();
137 module = single_module->GetModule_mrpc();
138 //cout<<"partId= "<<partId<<" module= "<<module<<endl;
139
140 //Process the hits
141 int nstrip = m_param.nstrip;
142 StripStruct stripStruct[12];
143 for(int i=0; i<nstrip; i++)
144 {
145 stripStruct[i].m_param = m_param;
146 stripStruct[i].setPar(m_V, m_threshold, m_saturationFlag);
147 }
148
149 BesTofHit* hit;
150 for(unsigned int i=0; i<single_module->GetHitIndexes_mrpc()->size(); i++)
151 {
152 hit = (*m_THC)[ (*(single_module->GetHitIndexes_mrpc()))[i] ];
153
154 HitStruct hitStruct;
155 hitStruct.m_param = m_param;
156 hitStruct.trkIndex = hit->GetG4Index();
157 hitStruct.pdgCode = hit->GetPDGcode();
158 hitStruct.ions = hit->GetIons();
159 hitStruct.strip = calStrip(hit->GetLocPos().z()/mm);
160 hitStruct.underStrip = underStrip(hit->GetLocPos().x()/mm, hit->GetLocPos().z()/mm);
161 hitStruct.gap = hit->GetGapNb();
162 hitStruct.glbTime = hit->GetTime()/ns;
163 hitStruct.locx = hit->GetLocPos().x()/mm;
164 hitStruct.locy = hit->GetLocPos().y()/mm;
165 hitStruct.locz = hit->GetLocPos().z()/mm;
166 hitStruct.x = hit->GetPos().x()/mm;
167 hitStruct.y = hit->GetPos().y()/mm;
168 hitStruct.z = hit->GetPos().z()/mm;
169 hitStruct.px = hit->GetMomentum().x()/(GeV/(3e+08*m/s));
170 hitStruct.py = hit->GetMomentum().y()/(GeV/(3e+08*m/s));
171 hitStruct.pz = hit->GetMomentum().z()/(GeV/(3e+08*m/s));
172 //hitStruct.print();
173
174 if(hitStruct.ions>0 && hitStruct.glbTime>0)
175 {
176 stripStruct[hitStruct.strip].hitStructCol.push_back(hitStruct);
177 }
178 }
179
180 for(int i=0; i<nstrip; i++)
181 {
182 if(stripStruct[i].hitStructCol.size()==0) continue;
183 stripStruct[i].strip = i;
184 stripStruct[i].calFirstHit();
185 stripStruct[i].avalanche();
186 //stripStruct[i].print();
187
188 if(stripStruct[i].tThreshold<=0) continue;
189
190 tdc_sphi = stripStruct[i].tStart + stripStruct[i].tThreshold + stripStruct[i].tPropagate_sphi;
191 tdc_xphi = stripStruct[i].tStart + stripStruct[i].tThreshold + stripStruct[i].tPropagate_xphi;
192
193 double tdcRes_charge;
194 if(m_calTdcRes_charge_flag==0)
195 {
196 tdcRes_charge = calTdcRes_charge(stripStruct[i].charge*1000); //ps, charge in fC
197 }
198 else if(m_calTdcRes_charge_flag==1)
199 {
200 tdcRes_charge = calTdcRes_charge1(stripStruct[i].charge*1000);
201 }
202 else if(m_calTdcRes_charge_flag==2)
203 {
204 tdcRes_charge = 0;
205 }
206
207 tdcRes = sqrt(tdcRes_charge*tdcRes_charge+tdcRes_const*tdcRes_const)/1000; //ns
208
209 tdc_sphi = G4RandGauss::shoot(tdc_sphi, tdcRes);
210 tdc_xphi = G4RandGauss::shoot(tdc_xphi, tdcRes);
211
212 if(m_charge2Time_flag==0)
213 {
214 adc = charge2Time(stripStruct[i].charge*1000); //ns, charge in fC
215 }
216 else if(m_charge2Time_flag==1)
217 {
218 adc = charge2Time1(stripStruct[i].charge*1000);
219 }
220
221 double adcRes_charge;
222 if(m_calAdcRes_charge_flag==0)
223 {
224 adcRes_charge = calAdcRes_charge(stripStruct[i].charge*1000); //ps, charge in fC
225 }
226 else if(m_calAdcRes_charge_flag==1)
227 {
228 adcRes_charge = calAdcRes_charge1(stripStruct[i].charge*1000);
229 }
230 else if(m_calAdcRes_charge_flag==2)
231 {
232 adcRes_charge = 0;
233 }
234
235 adcRes = sqrt(adcRes_charge*adcRes_charge+adcRes_const*adcRes_const)/1000;
236 adc = G4RandGauss::shoot(adc, adcRes);
237
238 time_leading_sphi = tdc_sphi;
239 time_leading_xphi = tdc_xphi;
240 time_trailing_sphi = tdc_sphi+adc;
241 time_trailing_xphi = tdc_xphi+adc;
242
243
244 //save digi information
245 BesTofDigi *digi = new BesTofDigi;
246 digi->SetTrackIndex(stripStruct[i].trkIndex);
247 digi->SetPartId(partId);
248 digi->SetModule(module);
249 digi->SetStrip(stripStruct[i].strip);
250 digi->SetForwT1(time_leading_sphi);
251 digi->SetForwT2(time_trailing_sphi);
252 digi->SetBackT1(time_leading_xphi);
253 digi->SetBackT2(time_trailing_xphi);
254 m_besTofDigitsCollection->insert(digi);
255 //cout<<"Print digi info: "
256 // <<" partId= "<<partId
257 // <<" module= "<<module
258 // <<" strip= "<<stripStruct[i].strip
259 // <<" time_leading_sphi= "<<time_leading_sphi
260 // <<" time_leading_xphi= "<<time_leading_xphi
261 // <<" time_trailing_sphi= "<<time_trailing_sphi
262 // <<" time_trailing_xphi= "<<time_trailing_xphi
263 // <<endl;
264
265
266 //save digi information
267 if(m_rootFlag)
268 {
269 m_partId = partId;
270 m_module = module;
271 m_time_leading_sphi = time_leading_sphi;
272 m_time_leading_xphi = time_leading_xphi;
273 m_time_trailing_sphi = time_trailing_sphi;
274 m_time_trailing_xphi = time_trailing_xphi;
275 m_tdcRes = tdcRes;
276 m_tdcRes_charge = tdcRes_charge;
277 m_adc = adc;
278 m_adcRes = adcRes;
279 m_adcRes_charge = adcRes_charge;
280
281 m_strip = stripStruct[i].strip;
282 m_trkIndex = stripStruct[i].trkIndex;
283 m_tStart = stripStruct[i].tStart;
284 m_tPropagate_sphi = stripStruct[i].tPropagate_sphi;
285 m_tPropagate_xphi = stripStruct[i].tPropagate_xphi;
286 m_tThreshold = stripStruct[i].tThreshold;
287 m_charge = stripStruct[i].charge;
288
289 m_nhit = stripStruct[i].hitStructCol.size();
290 //cout<<"m_nhit= "<<m_nhit<<endl;
291 for(int j=0; j<m_nhit; j++)
292 {
293 m_ions_hit[j] = stripStruct[i].hitStructCol[j].ions;
294 m_trkIndex_hit[j] = stripStruct[i].hitStructCol[j].trkIndex;
295 m_pdgCode_hit[j] = stripStruct[i].hitStructCol[j].pdgCode;
296 m_gap_hit[j] = stripStruct[i].hitStructCol[j].gap;
297 m_underStrip_hit[j] = stripStruct[i].hitStructCol[j].underStrip;
298 m_locx_hit[j] = stripStruct[i].hitStructCol[j].locx;
299 m_locy_hit[j] = stripStruct[i].hitStructCol[j].locy;
300 m_locz_hit[j] = stripStruct[i].hitStructCol[j].locz;
301 m_x_hit[j] = stripStruct[i].hitStructCol[j].x;
302 m_y_hit[j] = stripStruct[i].hitStructCol[j].y;
303 m_z_hit[j] = stripStruct[i].hitStructCol[j].z;
304 m_px_hit[j] = stripStruct[i].hitStructCol[j].px;
305 m_py_hit[j] = stripStruct[i].hitStructCol[j].py;
306 m_pz_hit[j] = stripStruct[i].hitStructCol[j].pz;
307 }
308 m_tree->Fill();
309 }
310 }
311}
312
313int BesTofDigitizerEcV4::calStrip(double locZ)
314{
315 int strip=-1;
316 double stripWidth = m_param.strip_z+m_param.strip_gap; //Strip spread: (24+3)mm
317 int nstrip = m_param.nstrip;
318 //the offset between strip coordinate and gas SD: 0.5
319 double length = locZ+stripWidth*nstrip/2-0.5;
320 if(length<=0)
321 {
322 strip=0;
323 }
324 else if(length<stripWidth*nstrip)
325 {
326 for(int i=0; i<nstrip; i++)
327 {
328 if(length>i*stripWidth && length<(i+1)*stripWidth)
329 {
330 strip = i;
331 break;
332 }
333 }
334 }
335 else
336 {
337 strip=nstrip-1;
338 }
339 if(strip<0) strip=0;
340 if(strip>nstrip-1) strip=nstrip-1;
341
342 return strip;
343}
344
345bool BesTofDigitizerEcV4::underStrip(double locX, double locZ)
346{
347 bool flag = 0;
348 int strip=-1;
349 double stripWidth = m_param.strip_z+m_param.strip_gap; //Strip spread: (24+3)mm
350 int nstrip = m_param.nstrip;
351 //the offset between strip coordinate and gas SD: 0.5
352 double length = locZ+stripWidth*nstrip/2-0.5;
353 if(length<stripWidth*nstrip)
354 {
355 for(int i=0; i<nstrip; i++)
356 {
357 if(length>i*stripWidth && length<(i+1)*stripWidth)
358 {
359 strip = i;
360 if(length>i*stripWidth+m_param.strip_gap/2 && length<(i+1)*stripWidth-m_param.strip_gap/2
361 && locX>-m_param.strip_x[strip]/2 && locX<m_param.strip_x[strip]/2) flag=1;
362 break;
363 }
364 }
365 }
366
367 return flag;
368}
369
370
371void BesTofDigitizerEcV4::StripStruct::calFirstHit()
372{
373 for(unsigned int i=0; i<hitStructCol.size(); i++)
374 {
375 if(hitStructCol[i].glbTime<tStart)
376 {
377 tStart = hitStructCol[i].glbTime;
378 trkIndex = hitStructCol[i].trkIndex;
379 hitStructCol[i].calTPropagate();
380 tPropagate_sphi = hitStructCol[i].tPropagate_sphi;
381 tPropagate_xphi = hitStructCol[i].tPropagate_xphi;
382 }
383 }
384}
385
386void BesTofDigitizerEcV4::HitStruct::calTPropagate()
387{
388 if(strip<0 || strip>m_param.nstrip-1)
389 {
390 cout<<"!! BesTofDigitizerEcV4::HitStruct::calTPropagate Wrong Strip !!!"<<endl;
391 return;
392 }
393
394 //It can be minus, consistent with calibration
395 double length_sphi = m_param.strip_x[strip]/2-locx; //mm
396 tPropagate_sphi = abs(length_sphi)/v_propagate;
397
398 double length_xphi = m_param.strip_x[strip]/2+locx; //mm
399 tPropagate_xphi = abs(length_xphi)/v_propagate;
400}
401
402double BesTofDigitizerEcV4::HitStruct::calAvaLength()
403{
404 //This calculation depends on the arangements of the gasLayer order and the turnover of gasContainer.
405 //all modules have the same local y trends: y larger, 11->0
406 //In units of mm
407 double length=0;
408 if(gap>=0 && gap<m_param.ngap/2) length = m_param.gapWidth/2+locy;
409 else if(gap<m_param.ngap) length = m_param.gapWidth/2-locy;
410 else
411 {
412 cout<<"BesTofDigitizerEcV4::StripStruct::calAvaLength Wrong gap calculation !!!"<<endl;
413 return -999.0;
414 }
415
416 return length;
417}
418
419void BesTofDigitizerEcV4::StripStruct::avalanche()
420{
421 //process each hit
422 for(unsigned int i=0; i<hitStructCol.size(); i++)
423 {
424 hitStructCol[i].ava_pos.clear();
425 hitStructCol[i].ava_num.clear();
426
427 hitStructCol[i].ava_pos[0] = hitStructCol[i].calAvaLength();
428 hitStructCol[i].ava_num[0] = hitStructCol[i].ions;
429 //cout<<"i= "<<i<<" gap= "<<hitStructCol[i].gap<<" initial pos= "<<hitStructCol[i].ava_pos[0]<<endl;
430 for(int j=1; j<m_param.nstep; j++)
431 {
432 hitStructCol[i].ava_pos[j] = hitStructCol[i].ava_pos[j-1] + m_param.stepWidth;
433 if(saturationFlag==true && hitStructCol[i].ava_num[j-1]>1.5e+07) //saturation e+07~e+08, ~2pC, Reather limit
434 {
435 hitStructCol[i].ava_num[j] = hitStructCol[i].ava_num[j-1];
436 }
437 else
438 {
439 hitStructCol[i].ava_num[j] = calNextN(hitStructCol[i].ava_num[j-1]);
440 }
441 if(hitStructCol[i].ava_pos[j]>m_param.gapWidth) break;
442 }
443 }
444
445 //decide threshold and charge
446 bool over_threshold = false;
447 long int sum = 0;
448 for(int i=0; i<m_param.nstep; i++)
449 {
450 for(unsigned int j=0; j<hitStructCol.size(); j++)
451 {
452 if(i<hitStructCol[j].ava_pos.size() && hitStructCol[j].ava_pos[i]<m_param.gapWidth)
453 {
454 sum += hitStructCol[j].ava_num[i];
455 }
456 }
457 //cout<<"sum= "<<sum<<" avaSize= "<<hitStructCol.size()<<endl;
458
459 if(over_threshold==false)
460 {
461 if(sum>threshold)
462 {
463 over_threshold = true;
464 tThreshold = (m_param.gapWidth)/(m_param.nstep)/v_drift*(i+1);
465 }
466 }
467 }
468
469 charge = sum*(m_param.E_weight)*(v_drift)*(m_param.eCharge)*(m_param.gapWidth)/(m_param.nstep)/v_drift;
470}
471
472
473long int BesTofDigitizerEcV4::StripStruct::calNextN(int num)
474{
475 if(num<150)
476 {
477 long int nextN = 0;
478 double rdm;
479 for(int i=0; i<num; i++)
480 {
481 rdm = G4UniformRand();
482 nextN += multiply(rdm);
483 }
484 return nextN;
485 }
486 else
487 {
488 double nbar = exp((alpha-eta)*m_param.stepWidth);
489 double sigma = calSigma();
490 double mean = num*nbar;
491 double resolution = G4RandGauss::shoot(0,(sqrt(num)*sigma));
492 long int nextN = mean+resolution;
493 return nextN;
494 }
495}
496
497long int BesTofDigitizerEcV4::StripStruct::multiply(double rdm)
498{
499 double nbar = exp((alpha-eta)*m_param.stepWidth);
500 double k = eta/alpha;
501 double rdm_border = k*(nbar-1)/(nbar-k);
502 if(rdm<rdm_border)
503 {
504 return 0;
505 }
506 else
507 {
508 long int number = 1.+1./log((nbar-1.)/(nbar-k))*log((nbar-k)*(rdm-1)/(k-1)/nbar);
509 return number;
510 }
511}
512
513double BesTofDigitizerEcV4::StripStruct::calSigma()
514{
515 double nbar = exp((alpha-eta)*m_param.stepWidth);
516 double k = eta/alpha;
517 double sigma = sqrt((1+k)/(1-k)*nbar*(nbar-1));
518 return sigma;
519}
520
521double BesTofDigitizerEcV4::calTdcRes_charge(double charge_fC)
522{
523 double time =0;
524 if( charge_fC>=13. && charge_fC<250)
525 {
526 time = 100.764*exp(-charge_fC*0.0413966+0.377154)+ 13.814;
527 }
528 else if(charge_fC>=250)
529 {
530 time = 12.8562+0.000507299*charge_fC;
531 }
532 return time; //ps
533}
534
535double BesTofDigitizerEcV4::charge2Time(double charge_fC)
536{
537 double time=0;
538 if(charge_fC>=13) time=-120.808/log(charge_fC*30.1306)+26.6024; //ns
539 return time;
540}
541
542double BesTofDigitizerEcV4::calAdcRes_charge(double charge_fC)
543{
544 double time = 0;
545 if( charge_fC>=13. && charge_fC<250)
546 {
547 time = 72.6005*exp(-charge_fC*0.0302626 + 1.49059) + 40.8627;
548 }
549 else if(charge_fC>=250)
550 {
551 time = 32.6233+0.00404149*charge_fC;
552 }
553 return time; //ps
554}
555
556double BesTofDigitizerEcV4::calTdcRes_charge1(double charge_fC)
557{
558 double result =0;
559 if( charge_fC > 50.)
560 {
561 result = 67.6737*exp(-charge_fC/50.9995-0.27755)+9.06223;
562 }
563 else
564 {
565 result = 176.322-2.98345*charge_fC;
566 }
567 return result; //ps
568}
569
570double BesTofDigitizerEcV4::charge2Time1(double charge_fC)
571{
572 double time=0;
573 if(charge_fC>=13) time=-4.50565/log(charge_fC*0.0812208)+16.6653; //ns
574 return time;
575}
576
577double BesTofDigitizerEcV4::calAdcRes_charge1(double charge_fC)
578{
579 double time = 64.3326*exp(-charge_fC/25.4638 + 0.944184)+19.4846;
580 return time; //ps
581}
582
583
584BesTofDigitizerEcV4::HitStruct::HitStruct()
585{
586 initial();
587}
588
589//in units of mm or ns
590void BesTofDigitizerEcV4::HitStruct::initial()
591{
592 trkIndex = -999.0;
593 pdgCode = -999.0;
594 ions = -999.0;
595 strip = -999.0;
596 gap = -999.0;
597 glbTime = -999.0;
598 locx = -999.0;
599 locy = -999.0;
600 locz = -999.0;
601 x = -999.0;
602 y = -999.0;
603 z = -999.0;
604 px = -999.0;
605 py = -999.0;
606 pz = -999.0;
607 v_propagate = 0.5*0.299792458e+3; //mm/ns
608 tPropagate_sphi = -999.0;
609 tPropagate_xphi = -999.0;
610}
611
612BesTofDigitizerEcV4::StripStruct::StripStruct()
613{
614 initial();
615}
616
617void BesTofDigitizerEcV4::StripStruct::initial()
618{
619 //properties to get
620 strip = -999.0;
621 trkIndex = -999.0;
622 tStart = 99999.0;
623 tPropagate_sphi = -999.0;
624 tPropagate_xphi = -999.0;
625 tThreshold = -999.0;
626 charge = -999.0;
627
628 //parameters to tune
629 E = 106;
630 alpha = 144800./1000; //-999.0; /mm^-1
631 eta = 5013./1000; //-999.0; /mm^-1
632 threshold = 1.5e+08; //Correspond to induced charge of 15 fC
633 v_drift = 210.9e-3; // mm/ns
634
635 hitStructCol.clear();
636}
637
638void BesTofDigitizerEcV4::StripStruct::setPar(double E_V, double threshold_n, bool saturationFlag_n)
639{
640 threshold = threshold_n;
641 E = E_V/1000*2/6/(m_param.gapWidth/10); //kV/cm
642 alpha = getAlpha(E); //mm^-1
643 eta = getEta(E); //mm^-1
644 v_drift = getV(E); // mm/ns
645
646 saturationFlag = saturationFlag_n;
647}
648
649void BesTofDigitizerEcV4::Param::setPar(int nstep_n, double E_weight_n)
650{
651 if(nstep_n>0) nstep = nstep_n;
652 if(E_weight_n>0) E_weight = E_weight_n;
653}
654
655BesTofDigitizerEcV4::Param::Param()
656{
657 //parameters fixed
658 tofPara = BesTofGeoParameter::GetInstance();
659 nstrip = 12;
660 std::stringstream ss;
661 for(int i=0; i<nstrip; i++)
662 {
663 ss.str("");
664 ss<<"strip_x["<<i<<"]";
665 strip_x[i] = tofPara->Get(ss.str().c_str()); //mm
666 }
667 strip_z = tofPara->Get("strip_z");
668 strip_gap = tofPara->Get("strip_gap");
669
670 ngap = 12;
671 gapWidth = 0.22; //mm
672 nstep = 200;
673 stepWidth = gapWidth/nstep;
674 E_weight = 1./(6.*0.22+(5.*0.4+2*0.55)/3.7); //V/mm
675 eCharge = 1.60217733e-7; //pC
676
677 //print();
678}
679
680double BesTofDigitizerEcV4::StripStruct::getAlpha(double E)
681{
682 //electric field: kV/cm; alpha: mm-1
683 //kV/cm
684 double e[22] =
685 {
686 65,
687 70 ,
688 75 ,
689 80 ,
690 85 ,
691 90 ,
692 95 ,
693 100 ,
694 102 ,
695 104 ,
696 106 ,
697 108 ,
698 110 ,
699 112 ,
700 114 ,
701 116 ,
702 118 ,
703 120 ,
704 125 ,
705 130 ,
706 135 ,
707 140
708 };
709
710 //mm-1
711 double alpha[22]=
712 {
713 383.5/10 ,
714 471 /10 ,
715 564.5/10 ,
716 663.6/10 ,
717 777.1/10 ,
718 877 /10 ,
719 990.8/10 ,
720 1106 /10 ,
721 1154 /10 ,
722 1199 /10 ,
723 1253 /10 ,
724 1296 /10 ,
725 1344 /10 ,
726 1396 /10 ,
727 1448 /10 ,
728 1502 /10 ,
729 1545 /10 ,
730 1597 /10 ,
731 1726 /10 ,
732 1858 /10 ,
733 1992 /10 ,
734 2124 /10 ,
735 };
736
737 TSpline3* sp_alpha = new TSpline3("sp_alpha", e, alpha, 22);
738 double alphaVal = sp_alpha->Eval(E);
739 return alphaVal;
740}
741
742double BesTofDigitizerEcV4::StripStruct::getEta(double E)
743{
744 //electric field: kV/cm; eta: mm-1
745 //kV/cm
746 double e[22] =
747 {
748 65,
749 70 ,
750 75 ,
751 80 ,
752 85 ,
753 90 ,
754 95 ,
755 100 ,
756 102 ,
757 104 ,
758 106 ,
759 108 ,
760 110 ,
761 112 ,
762 114 ,
763 116 ,
764 118 ,
765 120 ,
766 125 ,
767 130 ,
768 135 ,
769 140
770 };
771
772 //mm-1
773 double eta[22]=
774 {
775 132.6/10 ,
776 117.2/10 ,
777 102.6/10 ,
778 88.26/10 ,
779 79.81/10 ,
780 74.0 /10 ,
781 66.7 /10 ,
782 62.7 /10 ,
783 61.4 /10 ,
784 57.4 /10 ,
785 55.45/10 ,
786 54.35/10 ,
787 52.48/10 ,
788 51.3 /10 ,
789 50.1 /10 ,
790 48.3 /10 ,
791 48.28/10 ,
792 46.00/10 ,
793 44.08/10 ,
794 41.67/10 ,
795 39.97/10 ,
796 38.04/10
797 };
798
799 TSpline3* sp_eta = new TSpline3("sp_eta", e, eta, 22);
800 double etaVal = sp_eta->Eval(E);
801 return etaVal;
802}
803
804double BesTofDigitizerEcV4::StripStruct::getV(double E)
805{
806 //electric field: kV/cm; velocity: mm/ns
807 //kV/cm
808 double e[22] =
809 {
810 65,
811 70 ,
812 75 ,
813 80 ,
814 85 ,
815 90 ,
816 95 ,
817 100 ,
818 102 ,
819 104 ,
820 106 ,
821 108 ,
822 110 ,
823 112 ,
824 114 ,
825 116 ,
826 118 ,
827 120 ,
828 125 ,
829 130 ,
830 135 ,
831 140
832 };
833
834 //mm/ns
835 double v[22]=
836 {
837 130.2/1000 ,
838 138.5/1000 ,
839 146.7/1000 ,
840 155.0/1000 ,
841 163.3/1000 ,
842 171.4/1000 ,
843 179.7/1000 ,
844 187.7/1000 ,
845 191.2/1000 ,
846 194.5/1000 ,
847 197.9/1000 ,
848 201.2/1000 ,
849 204.5/1000 ,
850 207.6/1000 ,
851 210.9/1000 ,
852 214.4/1000 ,
853 217.5/1000 ,
854 220.9/1000 ,
855 228.8/1000 ,
856 237.0/1000 ,
857 244.7/1000 ,
858 252.9/1000
859 };
860
861 TSpline3* sp_v = new TSpline3("sp_v", e, v, 22);
862 double vVal = sp_v->Eval(E);
863 return vVal;
864}
865
866void BesTofDigitizerEcV4::Param::print()
867{
868 cout<<"Fixed parameters: "<<endl;
869 for(int i=0; i<nstrip; i++)
870 {
871 cout<<" strip_x["<<i<<"]= "<<strip_x[i];
872 }
873 cout<<" strip_z= "<<strip_z
874 <<" strip_gap= "<<strip_gap
875 <<" ngap= "<<ngap
876 <<" gapWidth= "<<gapWidth
877 <<" nstep= "<<nstep
878 <<" stepWidth= "<<stepWidth
879 <<" E_weight= "<<E_weight
880 <<" eCharge= "<<eCharge
881 <<endl;
882}
883
884void BesTofDigitizerEcV4::HitStruct::print()
885{
886 cout<<"Hit information: "<<endl;
887 cout<<" trkIndex= "<<trkIndex
888 <<" pdgCode= "<<pdgCode
889 <<" ions= "<<pdgCode
890 <<" strip= "<<strip
891 <<" gap= "<<gap
892 <<" glbTime= "<<glbTime
893 <<" locx= "<<locx
894 <<" locy= "<<locy
895 <<" locz= "<<locz
896 <<" x= "<<x
897 <<" y= "<<y
898 <<" z= "<<z
899 <<" px= "<<px
900 <<" py= "<<py
901 <<" pz= "<<pz
902 <<" v_propagate= "<<v_propagate
903 <<" tPropagate_sphi= "<<tPropagate_sphi
904 <<" tPropagate_xphi= "<<tPropagate_xphi
905 <<endl;
906}
907
908void BesTofDigitizerEcV4::StripStruct::print()
909{
910 cout<<"Strip information: "<<endl;
911 cout<<" strip= "<<strip
912 <<" trkIndex= "<<trkIndex
913 <<" tStart= "<<tStart
914 <<" tPropagate_sphi= "<<tPropagate_sphi
915 <<" tPropagate_xphi= "<<tPropagate_xphi
916 <<" tThreshold "<<tThreshold
917 <<" charge= "<<charge
918 <<" E= "<<E
919 <<" alpha= "<<alpha
920 <<" eta= "<<eta
921 <<" threshold= "<<threshold
922 <<" v_drift= "<<v_drift
923 <<endl;
924}
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
time
Double_t time
Definition: DataBase/tau_mode.c:7
exp
EvtComplex exp(const EvtComplex &c)
Definition: EvtComplex.hh:252
alpha
const double alpha
Definition: FastVertexFit.cxx:4
s
XmlRpcServer s
Definition: HelloServer.cpp:11
BesTofDigitsCollection
G4TDigiCollection< BesTofDigi > BesTofDigitsCollection
Definition: InstallArea/include/TofSim/TofSim/BesTofDigi.hh:79
BesTofHitsCollection
G4THitsCollection< BesTofHit > BesTofHitsCollection
Definition: InstallArea/include/TofSim/TofSim/BesTofHit.hh:113
v
**********Class see also m_nmax DOUBLE PRECISION m_amel DOUBLE PRECISION m_x2 DOUBLE PRECISION m_alfinv DOUBLE PRECISION m_Xenph INTEGER m_KeyWtm INTEGER m_idyfs DOUBLE PRECISION m_zini DOUBLE PRECISION m_q2 DOUBLE PRECISION m_Wt_KF DOUBLE PRECISION m_WtCut INTEGER m_KFfin *COMMON c_KarLud $ !Input CMS energy[GeV] $ !CMS energy after beam spread beam strahlung[GeV] $ !Beam energy spread[GeV] $ !z boost due to beam spread $ !electron beam mass *ff pair spectrum $ !minimum v
Definition: KarLud.h:35
BesTofDigitizerEcV4::calTdcRes_charge1
double calTdcRes_charge1(double charge_fC)
Definition: BesTofDigitizerEcV4.cc:556
BesTofDigitizerEcV4::calAdcRes_charge
double calAdcRes_charge(double charge_fC)
Definition: BesTofDigitizerEcV4.cc:542
BesTofDigitizerEcV4::calAdcRes_charge1
double calAdcRes_charge1(double charge_fC)
Definition: BesTofDigitizerEcV4.cc:577
BesTofDigitizerEcV4::underStrip
bool underStrip(double locX, double locZ)
Definition: BesTofDigitizerEcV4.cc:345
BesTofDigitizerEcV4::calStrip
int calStrip(double locZ)
Definition: BesTofDigitizerEcV4.cc:313
BesTofDigitizerEcV4::Digitize
virtual void Digitize(ScintSingle *, BesTofDigitsCollection *)
Definition: BesTofDigitizerEcV4.cc:128
BesTofDigitizerEcV4::charge2Time1
double charge2Time1(double charge_fC)
Definition: BesTofDigitizerEcV4.cc:570
BesTofDigitizerEcV4::charge2Time
double charge2Time(double charge_fC)
Definition: BesTofDigitizerEcV4.cc:535
BesTofDigitizerEcV4::calTdcRes_charge
double calTdcRes_charge(double charge_fC)
Definition: BesTofDigitizerEcV4.cc:521
BesTofGeoParameter::GetInstance
static BesTofGeoParameter * GetInstance()
Definition: BesTofGeoParameter.cc:21
num
int num[96]
Definition: ranlxd.c:373
BesTofDigitizerEcV4::Param::setPar
void setPar(int nstep, double E_weight)
Definition: BesTofDigitizerEcV4.cc:649
BesTofDigitizerEcV4::StripStruct::setPar
void setPar(double V, double threshold, bool saturationFlag=true)
Definition: BesTofDigitizerEcV4.cc:638
ns
#define ns(x)
Definition: xmltok.c:1504