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