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