00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include "GaudiKernel/MsgStream.h"
00012 #include "GaudiKernel/Bootstrap.h"
00013 #include "GaudiKernel/PropertyMgr.h"
00014 #include "GaudiKernel/IJobOptionsSvc.h"
00015 #include "GaudiKernel/ISvcLocator.h"
00016 #include "GaudiKernel/IDataProviderSvc.h"
00017
00018 #include "BesTofDigitizerEcV4.hh"
00019 #include "BesTofDigi.hh"
00020 #include "BesTofHit.hh"
00021 #include "G4DigiManager.hh"
00022 #include "Randomize.hh"
00023 #include "TMath.h"
00024 #include <math.h>
00025 #include "TNtuple.h"
00026 #include "TFile.h"
00027 #include "TH1F.h"
00028 #include "TROOT.h"
00029 #include "TSpline.h"
00030
00031
00032 BesTofDigitizerEcV4::BesTofDigitizerEcV4()
00033 {
00034 PropertyMgr m_propMgr;
00035 m_propMgr.declareProperty("FileName", m_fileName = "mrpc.root");
00036 m_propMgr.declareProperty("RootFlag", m_rootFlag = false);
00037 m_propMgr.declareProperty("E", m_V = 7000);
00038 m_propMgr.declareProperty("Threshold", m_threshold = 5.5e+08);
00039
00040 m_propMgr.declareProperty("nstep", m_nstep = -1);
00041 m_propMgr.declareProperty("E_weight", m_E_weight = -1);
00042 m_propMgr.declareProperty("saturationFlag", m_saturationFlag = true);
00043 m_propMgr.declareProperty("tdcRes_const", tdcRes_const = -1);
00044 m_propMgr.declareProperty("adcRes_const", adcRes_const = -1);
00045 m_propMgr.declareProperty("calTdcRes_charge_flag", m_calTdcRes_charge_flag=0);
00046 m_propMgr.declareProperty("charge2Time_flag", m_charge2Time_flag=0);
00047 m_propMgr.declareProperty("calAdcRes_charge_flag", m_calAdcRes_charge_flag=0);
00048 cout<<"Property:"<<endl
00049 <<" FileName= "<<m_fileName
00050 <<" E= "<<m_V
00051 <<" Threshold= "<<m_threshold
00052 <<" nstep= "<<m_nstep
00053 <<" E_weight= "<<m_E_weight
00054 <<" saturationFlag= "<<m_saturationFlag
00055 <<" tdcRes_const= "<<tdcRes_const
00056 <<" adcRes_const= "<<adcRes_const
00057 <<" calTdcRes_charge_flag= "<<m_calTdcRes_charge_flag
00058 <<" charge2Time_flag= "<<m_charge2Time_flag
00059 <<" calAdcRes_charge_flag= "<<m_calAdcRes_charge_flag
00060 <<endl;
00061
00062
00063 IJobOptionsSvc* jobSvc;
00064 Gaudi::svcLocator()->service("JobOptionsSvc", jobSvc);
00065 jobSvc->setMyProperties("BesTofDigitizerEcV4", &m_propMgr);
00066
00067 initial();
00068
00069 if(m_rootFlag)
00070 {
00071 m_file = new TFile(m_fileName.c_str(), "RECREATE");
00072 m_tree = new TTree("mrpc", "mrpc");
00073
00074 m_tree->Branch("event", &m_event, "event/D");
00075 m_tree->Branch("partId", &m_partId, "partId/D");
00076 m_tree->Branch("module", &m_module, "module/D");
00077 m_tree->Branch("time_leading_sphi", &m_time_leading_sphi, "time_leading_sphi/D");
00078 m_tree->Branch("time_leading_xphi", &m_time_leading_xphi, "time_leading_xphi/D");
00079 m_tree->Branch("time_trailing_sphi", &m_time_trailing_sphi, "time_trailing_sphi/D");
00080 m_tree->Branch("time_trailing_xphi", &m_time_trailing_xphi, "time_trailing_xphi/D");
00081 m_tree->Branch("tdcRes", &m_tdcRes, "tdcRes/D");
00082 m_tree->Branch("tdcRes_charge", &m_tdcRes_charge, "tdcRes_charge/D");
00083 m_tree->Branch("adc", &m_adc, "adc/D");
00084 m_tree->Branch("adcRes", &m_adcRes, "adcRes/D");
00085 m_tree->Branch("adcRes_charge", &m_adcRes_charge, "adcRes_charge/D");
00086 m_tree->Branch("strip", &m_strip, "strip/D");
00087 m_tree->Branch("trkIndex", &m_trkIndex, "trkIndex/D");
00088 m_tree->Branch("tStart", &m_tStart, "tStart/D");
00089 m_tree->Branch("tPropagate_sphi", &m_tPropagate_sphi, "tPropagate_sphi/D");
00090 m_tree->Branch("tPropagate_xphi", &m_tPropagate_xphi, "tPropagate_xphi/D");
00091 m_tree->Branch("tThreshold", &m_tThreshold, "tThreshold/D");
00092 m_tree->Branch("charge", &m_charge, "charge/D");
00093 m_tree->Branch("nhit", &m_nhit, "nhit/I");
00094 m_tree->Branch("ions_hit", m_ions_hit, "ions_hit[nhit]/D");
00095 m_tree->Branch("trkIndex_hit", m_trkIndex_hit, "trkIndex_hit[nhit]/D");
00096 m_tree->Branch("pdgCode_hit", m_pdgCode_hit, "pdgCode_hit[nhit]/D");
00097 m_tree->Branch("gap_hit", m_gap_hit, "gap_hit[nhit]/D");
00098 m_tree->Branch("underStrip_hit", m_underStrip_hit, "underStrip_hit[nhit]/D");
00099 m_tree->Branch("locx_hit", m_locx_hit, "locx_hit[nhit]/D");
00100 m_tree->Branch("locy_hit", m_locy_hit, "locy_hit[nhit]/D");
00101 m_tree->Branch("locz_hit", m_locz_hit, "locz_hit[nhit]/D");
00102 m_tree->Branch("x_hit", m_x_hit, "x_hit[nhit]/D");
00103 m_tree->Branch("y_hit", m_y_hit, "y_hit[nhit]/D");
00104 m_tree->Branch("z_hit", m_z_hit, "z_hit[nhit]/D");
00105 m_tree->Branch("px_hit", m_px_hit, "px_hit[nhit]/D");
00106 m_tree->Branch("py_hit", m_py_hit, "py_hit[nhit]/D");
00107 m_tree->Branch("pz_hit", m_pz_hit, "pz_hit[nhit]/D");
00108 }
00109 }
00110
00111
00112 BesTofDigitizerEcV4::~BesTofDigitizerEcV4()
00113 {
00114 if(m_rootFlag)
00115 {
00116 m_file->Write();
00117 m_file->Close();
00118 }
00119 }
00120
00121 void BesTofDigitizerEcV4::initial()
00122 {
00123 m_param = Param();
00124 m_param.setPar(m_nstep, m_E_weight);
00125 m_event=-1;
00126 partId = -999;
00127 module = -999;
00128 tdc_sphi = -999;
00129 tdc_xphi = -999;
00130 if(tdcRes_const<0) tdcRes_const = 38;
00131 tdcRes = -999;
00132
00133 adc = -999;
00134 if(adcRes_const<0) adcRes_const = 27;
00135 adcRes = -999;
00136
00137 time_leading_sphi = -999;
00138 time_leading_xphi = -999;
00139 time_trailing_sphi = -999;
00140 time_trailing_xphi = -999;
00141 }
00142
00143 void BesTofDigitizerEcV4::Digitize( ScintSingle* single_module, BesTofDigitsCollection* DC )
00144 {
00145 m_besTofDigitsCollection = DC;
00146 G4DigiManager* digiManager = G4DigiManager::GetDMpointer();
00147 G4int THCID = digiManager->GetHitsCollectionID("BesTofHitsCollection");
00148 m_THC = (BesTofHitsCollection*) (digiManager->GetHitsCollection(THCID));
00149 if( !m_THC ) return;
00150
00151 partId = single_module->GetPartId();
00152 module = single_module->GetModule_mrpc();
00153
00154
00155
00156 int nstrip = m_param.nstrip;
00157 StripStruct stripStruct[12];
00158 for(int i=0; i<nstrip; i++)
00159 {
00160 stripStruct[i].m_param = m_param;
00161 stripStruct[i].setPar(m_V, m_threshold, m_saturationFlag);
00162 }
00163
00164 BesTofHit* hit;
00165 for(unsigned int i=0; i<single_module->GetHitIndexes_mrpc()->size(); i++)
00166 {
00167 hit = (*m_THC)[ (*(single_module->GetHitIndexes_mrpc()))[i] ];
00168 m_event = hit->GetEvent();
00169
00170 HitStruct hitStruct;
00171 hitStruct.m_param = m_param;
00172 hitStruct.trkIndex = hit->GetG4Index();
00173 hitStruct.pdgCode = hit->GetPDGcode();
00174 hitStruct.ions = hit->GetIons();
00175 hitStruct.strip = calStrip(hit->GetLocPos().z()/mm);
00176 hitStruct.underStrip = underStrip(hit->GetLocPos().x()/mm, hit->GetLocPos().z()/mm);
00177 hitStruct.gap = hit->GetGapNb();
00178 hitStruct.glbTime = hit->GetTime()/ns;
00179 hitStruct.locx = hit->GetLocPos().x()/mm;
00180 hitStruct.locy = hit->GetLocPos().y()/mm;
00181 hitStruct.locz = hit->GetLocPos().z()/mm;
00182 hitStruct.x = hit->GetPos().x()/mm;
00183 hitStruct.y = hit->GetPos().y()/mm;
00184 hitStruct.z = hit->GetPos().z()/mm;
00185 hitStruct.px = hit->GetMomentum().x()/(GeV/(3e+08*m/s));
00186 hitStruct.py = hit->GetMomentum().y()/(GeV/(3e+08*m/s));
00187 hitStruct.pz = hit->GetMomentum().z()/(GeV/(3e+08*m/s));
00188
00189
00190 if(hitStruct.ions>0 && hitStruct.glbTime>0)
00191 {
00192 stripStruct[hitStruct.strip].hitStructCol.push_back(hitStruct);
00193 }
00194 }
00195
00196
00197 int hitSize=0;
00198
00199 for(int i=0; i<nstrip; i++)
00200 {
00201 if(stripStruct[i].hitStructCol.size()==0) continue;
00202 stripStruct[i].strip = i;
00203 stripStruct[i].calFirstHit();
00204 stripStruct[i].avalanche();
00205
00206
00207 if(stripStruct[i].tThreshold<=0) continue;
00208
00209 tdc_sphi = stripStruct[i].tStart + stripStruct[i].tThreshold + stripStruct[i].tPropagate_sphi;
00210 tdc_xphi = stripStruct[i].tStart + stripStruct[i].tThreshold + stripStruct[i].tPropagate_xphi;
00211
00212 double tdcRes_charge;
00213 if(m_calTdcRes_charge_flag==0)
00214 {
00215 tdcRes_charge = calTdcRes_charge(stripStruct[i].charge*1000);
00216 }
00217 else if(m_calTdcRes_charge_flag==1)
00218 {
00219 tdcRes_charge = calTdcRes_charge1(stripStruct[i].charge*1000);
00220 }
00221 else if(m_calTdcRes_charge_flag==2)
00222 {
00223 tdcRes_charge = 0;
00224 }
00225
00226 tdcRes = sqrt(tdcRes_charge*tdcRes_charge+tdcRes_const*tdcRes_const)/1000;
00227
00228 tdc_sphi = G4RandGauss::shoot(tdc_sphi, tdcRes);
00229 tdc_xphi = G4RandGauss::shoot(tdc_xphi, tdcRes);
00230
00231 if(m_charge2Time_flag==0)
00232 {
00233 adc = charge2Time(stripStruct[i].charge*1000);
00234 }
00235 else if(m_charge2Time_flag==1)
00236 {
00237 adc = charge2Time1(stripStruct[i].charge*1000);
00238 }
00239
00240 double adcRes_charge;
00241 if(m_calAdcRes_charge_flag==0)
00242 {
00243 adcRes_charge = calAdcRes_charge(stripStruct[i].charge*1000);
00244 }
00245 else if(m_calAdcRes_charge_flag==1)
00246 {
00247 adcRes_charge = calAdcRes_charge1(stripStruct[i].charge*1000);
00248 }
00249 else if(m_calAdcRes_charge_flag==2)
00250 {
00251 adcRes_charge = 0;
00252 }
00253
00254 adcRes = sqrt(adcRes_charge*adcRes_charge+adcRes_const*adcRes_const)/1000;
00255 adc = G4RandGauss::shoot(adc, adcRes);
00256 if(adc<0) adc=0;
00257
00258 time_leading_sphi = tdc_sphi;
00259 time_leading_xphi = tdc_xphi;
00260 time_trailing_sphi = tdc_sphi+adc;
00261 time_trailing_xphi = tdc_xphi+adc;
00262
00263
00264
00265 BesTofDigi *digi = new BesTofDigi;
00266 digi->SetTrackIndex(stripStruct[i].trkIndex);
00267 digi->SetPartId(partId);
00268 digi->SetModule(module);
00269 digi->SetStrip(stripStruct[i].strip);
00270 int mo = (partId-3)*36+module;
00271 int st = stripStruct[i].strip;
00272 if(m_param.deadChannel[mo][st]==0 || m_param.deadChannel[mo][st]==2)
00273 {
00274
00275 digi->SetForwT1(-999);
00276 digi->SetForwT2(-999);
00277 }
00278 else
00279 {
00280 digi->SetForwT1(time_leading_sphi);
00281 digi->SetForwT2(time_trailing_sphi);
00282 }
00283 if(m_param.deadChannel[mo][st]==1 || m_param.deadChannel[mo][st]==2)
00284 {
00285 digi->SetBackT1(-999);
00286 digi->SetBackT2(-999);
00287 }
00288 else
00289 {
00290 digi->SetBackT1(time_leading_xphi);
00291 digi->SetBackT2(time_trailing_xphi);
00292 }
00293 m_besTofDigitsCollection->insert(digi);
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306 if(m_rootFlag)
00307 {
00308 m_partId = partId;
00309 m_module = module;
00310 m_time_leading_sphi = time_leading_sphi;
00311 m_time_leading_xphi = time_leading_xphi;
00312 m_time_trailing_sphi = time_trailing_sphi;
00313 m_time_trailing_xphi = time_trailing_xphi;
00314 m_tdcRes = tdcRes;
00315 m_tdcRes_charge = tdcRes_charge;
00316 m_adc = adc;
00317 m_adcRes = adcRes;
00318 m_adcRes_charge = adcRes_charge;
00319
00320 m_strip = stripStruct[i].strip;
00321 m_trkIndex = stripStruct[i].trkIndex;
00322 m_tStart = stripStruct[i].tStart;
00323 m_tPropagate_sphi = stripStruct[i].tPropagate_sphi;
00324 m_tPropagate_xphi = stripStruct[i].tPropagate_xphi;
00325 m_tThreshold = stripStruct[i].tThreshold;
00326 m_charge = stripStruct[i].charge;
00327
00328 m_nhit = stripStruct[i].hitStructCol.size();
00329
00330 for(int j=0; j<m_nhit; j++)
00331 {
00332 m_ions_hit[j] = stripStruct[i].hitStructCol[j].ions;
00333 m_trkIndex_hit[j] = stripStruct[i].hitStructCol[j].trkIndex;
00334 m_pdgCode_hit[j] = stripStruct[i].hitStructCol[j].pdgCode;
00335 m_gap_hit[j] = stripStruct[i].hitStructCol[j].gap;
00336 m_underStrip_hit[j] = stripStruct[i].hitStructCol[j].underStrip;
00337 m_locx_hit[j] = stripStruct[i].hitStructCol[j].locx;
00338 m_locy_hit[j] = stripStruct[i].hitStructCol[j].locy;
00339 m_locz_hit[j] = stripStruct[i].hitStructCol[j].locz;
00340 m_x_hit[j] = stripStruct[i].hitStructCol[j].x;
00341 m_y_hit[j] = stripStruct[i].hitStructCol[j].y;
00342 m_z_hit[j] = stripStruct[i].hitStructCol[j].z;
00343 m_px_hit[j] = stripStruct[i].hitStructCol[j].px;
00344 m_py_hit[j] = stripStruct[i].hitStructCol[j].py;
00345 m_pz_hit[j] = stripStruct[i].hitStructCol[j].pz;
00346 }
00347 m_tree->Fill();
00348 }
00349 }
00350 }
00351
00352 int BesTofDigitizerEcV4::calStrip(double locZ)
00353 {
00354 int strip=-1;
00355 double stripWidth = m_param.strip_z+m_param.strip_gap;
00356 int nstrip = m_param.nstrip;
00357
00358 double length = locZ+stripWidth*nstrip/2-0.5;
00359 if(length<=0)
00360 {
00361 strip=0;
00362 }
00363 else if(length<stripWidth*nstrip)
00364 {
00365 for(int i=0; i<nstrip; i++)
00366 {
00367 if(length>i*stripWidth && length<(i+1)*stripWidth)
00368 {
00369 strip = i;
00370 break;
00371 }
00372 }
00373 }
00374 else
00375 {
00376 strip=nstrip-1;
00377 }
00378 if(strip<0) strip=0;
00379 if(strip>nstrip-1) strip=nstrip-1;
00380
00381 return strip;
00382 }
00383
00384 bool BesTofDigitizerEcV4::underStrip(double locX, double locZ)
00385 {
00386 bool flag = 0;
00387 int strip=-1;
00388 double stripWidth = m_param.strip_z+m_param.strip_gap;
00389 int nstrip = m_param.nstrip;
00390
00391 double length = locZ+stripWidth*nstrip/2-0.5;
00392 if(length<stripWidth*nstrip)
00393 {
00394 for(int i=0; i<nstrip; i++)
00395 {
00396 if(length>i*stripWidth && length<(i+1)*stripWidth)
00397 {
00398 strip = i;
00399 if(length>i*stripWidth+m_param.strip_gap/2 && length<(i+1)*stripWidth-m_param.strip_gap/2
00400 && locX>-m_param.strip_x[strip]/2 && locX<m_param.strip_x[strip]/2) flag=1;
00401 break;
00402 }
00403 }
00404 }
00405
00406 return flag;
00407 }
00408
00409
00410 void BesTofDigitizerEcV4::StripStruct::calFirstHit()
00411 {
00412 for(unsigned int i=0; i<hitStructCol.size(); i++)
00413 {
00414 if(hitStructCol[i].glbTime<tStart)
00415 {
00416 tStart = hitStructCol[i].glbTime;
00417 trkIndex = hitStructCol[i].trkIndex;
00418 hitStructCol[i].calTPropagate();
00419 tPropagate_sphi = hitStructCol[i].tPropagate_sphi;
00420 tPropagate_xphi = hitStructCol[i].tPropagate_xphi;
00421 }
00422 }
00423 }
00424
00425 void BesTofDigitizerEcV4::HitStruct::calTPropagate()
00426 {
00427 if(strip<0 || strip>m_param.nstrip-1)
00428 {
00429 cout<<"!! BesTofDigitizerEcV4::HitStruct::calTPropagate Wrong Strip !!!"<<endl;
00430 return;
00431 }
00432
00433
00434 double length_sphi = m_param.strip_x[strip]/2-locx;
00435 tPropagate_sphi = abs(length_sphi)/v_propagate;
00436
00437 double length_xphi = m_param.strip_x[strip]/2+locx;
00438 tPropagate_xphi = abs(length_xphi)/v_propagate;
00439 }
00440
00441 double BesTofDigitizerEcV4::HitStruct::calAvaLength()
00442 {
00443
00444
00445
00446 double length=0;
00447 if(gap>=0 && gap<m_param.ngap/2) length = m_param.gapWidth/2+locy;
00448 else if(gap<m_param.ngap) length = m_param.gapWidth/2-locy;
00449 else
00450 {
00451 cout<<"BesTofDigitizerEcV4::StripStruct::calAvaLength Wrong gap calculation !!!"<<endl;
00452 return -999.0;
00453 }
00454
00455 return length;
00456 }
00457
00458 void BesTofDigitizerEcV4::StripStruct::avalanche()
00459 {
00460
00461 for(unsigned int i=0; i<hitStructCol.size(); i++)
00462 {
00463 hitStructCol[i].ava_pos.clear();
00464 hitStructCol[i].ava_num.clear();
00465
00466 hitStructCol[i].ava_pos[0] = hitStructCol[i].calAvaLength();
00467 hitStructCol[i].ava_num[0] = hitStructCol[i].ions;
00468
00469 for(int j=1; j<m_param.nstep; j++)
00470 {
00471 hitStructCol[i].ava_pos[j] = hitStructCol[i].ava_pos[j-1] + m_param.stepWidth;
00472 if(saturationFlag==true && hitStructCol[i].ava_num[j-1]>1.5e+07)
00473 {
00474 hitStructCol[i].ava_num[j] = hitStructCol[i].ava_num[j-1];
00475 }
00476 else
00477 {
00478 hitStructCol[i].ava_num[j] = calNextN(hitStructCol[i].ava_num[j-1]);
00479 }
00480 if(hitStructCol[i].ava_pos[j]>m_param.gapWidth) break;
00481 }
00482 }
00483
00484
00485 bool over_threshold = false;
00486 long int sum = 0;
00487 for(int i=0; i<m_param.nstep; i++)
00488 {
00489 for(unsigned int j=0; j<hitStructCol.size(); j++)
00490 {
00491 if(i<hitStructCol[j].ava_pos.size() && hitStructCol[j].ava_pos[i]<m_param.gapWidth)
00492 {
00493 sum += hitStructCol[j].ava_num[i];
00494 }
00495 }
00496
00497
00498 if(over_threshold==false)
00499 {
00500 if(sum>threshold)
00501 {
00502 over_threshold = true;
00503 tThreshold = (m_param.gapWidth)/(m_param.nstep)/v_drift*(i+1);
00504 }
00505 }
00506 }
00507
00508 charge = sum*(m_param.E_weight)*(v_drift)*(m_param.eCharge)*(m_param.gapWidth)/(m_param.nstep)/v_drift;
00509 }
00510
00511
00512 long int BesTofDigitizerEcV4::StripStruct::calNextN(int num)
00513 {
00514 if(num<150)
00515 {
00516 long int nextN = 0;
00517 double rdm;
00518 for(int i=0; i<num; i++)
00519 {
00520 rdm = G4UniformRand();
00521 nextN += multiply(rdm);
00522 }
00523 return nextN;
00524 }
00525 else
00526 {
00527 double nbar = exp((alpha-eta)*m_param.stepWidth);
00528 double sigma = calSigma();
00529 double mean = num*nbar;
00530 double resolution = G4RandGauss::shoot(0,(sqrt(num)*sigma));
00531 long int nextN = mean+resolution;
00532 return nextN;
00533 }
00534 }
00535
00536 long int BesTofDigitizerEcV4::StripStruct::multiply(double rdm)
00537 {
00538 double nbar = exp((alpha-eta)*m_param.stepWidth);
00539 double k = eta/alpha;
00540 double rdm_border = k*(nbar-1)/(nbar-k);
00541 if(rdm<rdm_border)
00542 {
00543 return 0;
00544 }
00545 else
00546 {
00547 long int number = 1.+1./log((nbar-1.)/(nbar-k))*log((nbar-k)*(rdm-1)/(k-1)/nbar);
00548 return number;
00549 }
00550 }
00551
00552 double BesTofDigitizerEcV4::StripStruct::calSigma()
00553 {
00554 double nbar = exp((alpha-eta)*m_param.stepWidth);
00555 double k = eta/alpha;
00556 double sigma = sqrt((1+k)/(1-k)*nbar*(nbar-1));
00557 return sigma;
00558 }
00559
00560 double BesTofDigitizerEcV4::calTdcRes_charge(double charge_fC)
00561 {
00562 double time=0;
00563 if( charge_fC<250)
00564 {
00565 time = 100.764*exp(-charge_fC*0.0413966+0.377154)+ 13.814;
00566 }
00567 else
00568 {
00569 time = 12.8562+0.000507299*charge_fC;
00570 }
00571 if(time<0) time=0;
00572 return time;
00573 }
00574
00575 double BesTofDigitizerEcV4::charge2Time(double charge_fC)
00576 {
00577 double time=0;
00578 time=-120.808/log(charge_fC*30.1306)+26.6024;
00579 if(time<0) time=0;
00580 return time;
00581 }
00582
00583 double BesTofDigitizerEcV4::calAdcRes_charge(double charge_fC)
00584 {
00585 double time=0;
00586 if(charge_fC<250)
00587 {
00588 time = 72.6005*exp(-charge_fC*0.0302626 + 1.49059) + 40.8627;
00589 }
00590 else
00591 {
00592 time = 32.6233+0.00404149*charge_fC;
00593 }
00594 if(time<0) time=0;
00595 return time;
00596 }
00597
00598 double BesTofDigitizerEcV4::calTdcRes_charge1(double charge_fC)
00599 {
00600 double result =0;
00601 if( charge_fC > 50.)
00602 {
00603 result = 67.6737*exp(-charge_fC/50.9995-0.27755)+9.06223;
00604 }
00605 else
00606 {
00607 result = 176.322-2.98345*charge_fC;
00608 }
00609 if(result<0) result=0;
00610 return result;
00611 }
00612
00613 double BesTofDigitizerEcV4::charge2Time1(double charge_fC)
00614 {
00615 double time=0;
00616 time=-4.50565/log(charge_fC*0.0812208)+16.6653;
00617 if(time<0) time=0;
00618 return time;
00619 }
00620
00621 double BesTofDigitizerEcV4::calAdcRes_charge1(double charge_fC)
00622 {
00623 double time = 64.3326*exp(-charge_fC/25.4638 + 0.944184)+19.4846;
00624 if(time<0) time=0;
00625 return time;
00626 }
00627
00628
00629 BesTofDigitizerEcV4::HitStruct::HitStruct()
00630 {
00631 initial();
00632 }
00633
00634
00635 void BesTofDigitizerEcV4::HitStruct::initial()
00636 {
00637 trkIndex = -999.0;
00638 pdgCode = -999.0;
00639 ions = -999.0;
00640 strip = -999.0;
00641 gap = -999.0;
00642 glbTime = -999.0;
00643 locx = -999.0;
00644 locy = -999.0;
00645 locz = -999.0;
00646 x = -999.0;
00647 y = -999.0;
00648 z = -999.0;
00649 px = -999.0;
00650 py = -999.0;
00651 pz = -999.0;
00652 v_propagate = 0.5*0.299792458e+3;
00653 tPropagate_sphi = -999.0;
00654 tPropagate_xphi = -999.0;
00655 }
00656
00657 BesTofDigitizerEcV4::StripStruct::StripStruct()
00658 {
00659 initial();
00660 }
00661
00662 void BesTofDigitizerEcV4::StripStruct::initial()
00663 {
00664
00665 strip = -999.0;
00666 trkIndex = -999.0;
00667 tStart = 99999.0;
00668 tPropagate_sphi = -999.0;
00669 tPropagate_xphi = -999.0;
00670 tThreshold = -999.0;
00671 charge = -999.0;
00672
00673
00674 E = 106;
00675 alpha = 144800./1000;
00676 eta = 5013./1000;
00677 threshold = 1.5e+08;
00678 v_drift = 210.9e-3;
00679
00680 hitStructCol.clear();
00681 }
00682
00683 void BesTofDigitizerEcV4::StripStruct::setPar(double E_V, double threshold_n, bool saturationFlag_n)
00684 {
00685 threshold = threshold_n;
00686 E = E_V/1000*2/6/(m_param.gapWidth/10);
00687 alpha = getAlpha(E);
00688 eta = getEta(E);
00689 v_drift = getV(E);
00690
00691 saturationFlag = saturationFlag_n;
00692 }
00693
00694 void BesTofDigitizerEcV4::Param::setPar(int nstep_n, double E_weight_n)
00695 {
00696 if(nstep_n>0) nstep = nstep_n;
00697 if(E_weight_n>0) E_weight = E_weight_n;
00698 }
00699
00700 BesTofDigitizerEcV4::Param::Param()
00701 {
00702
00703 tofPara = BesTofGeoParameter::GetInstance();
00704 nstrip = 12;
00705 nmodule = 72;
00706 std::stringstream ss;
00707 for(int i=0; i<nstrip; i++)
00708 {
00709 ss.str("");
00710 ss<<"strip_x["<<i<<"]";
00711 strip_x[i] = tofPara->Get(ss.str().c_str());
00712 }
00713 strip_z = tofPara->Get("strip_z");
00714 strip_gap = tofPara->Get("strip_gap");
00715
00716 ngap = 12;
00717 gapWidth = 0.22;
00718 nstep = 200;
00719 stepWidth = gapWidth/nstep;
00720 E_weight = 1./(6.*0.22+(5.*0.4+2*0.55)/3.7);
00721 eCharge = 1.60217733e-7;
00722 tofPara->Get_deadChannel(deadChannel);
00723
00724
00725 }
00726
00727 double BesTofDigitizerEcV4::StripStruct::getAlpha(double E)
00728 {
00729
00730
00731 double e[22] =
00732 {
00733 65,
00734 70 ,
00735 75 ,
00736 80 ,
00737 85 ,
00738 90 ,
00739 95 ,
00740 100 ,
00741 102 ,
00742 104 ,
00743 106 ,
00744 108 ,
00745 110 ,
00746 112 ,
00747 114 ,
00748 116 ,
00749 118 ,
00750 120 ,
00751 125 ,
00752 130 ,
00753 135 ,
00754 140
00755 };
00756
00757
00758 double alpha[22]=
00759 {
00760 383.5/10 ,
00761 471 /10 ,
00762 564.5/10 ,
00763 663.6/10 ,
00764 777.1/10 ,
00765 877 /10 ,
00766 990.8/10 ,
00767 1106 /10 ,
00768 1154 /10 ,
00769 1199 /10 ,
00770 1253 /10 ,
00771 1296 /10 ,
00772 1344 /10 ,
00773 1396 /10 ,
00774 1448 /10 ,
00775 1502 /10 ,
00776 1545 /10 ,
00777 1597 /10 ,
00778 1726 /10 ,
00779 1858 /10 ,
00780 1992 /10 ,
00781 2124 /10 ,
00782 };
00783
00784 TSpline3* sp_alpha = new TSpline3("sp_alpha", e, alpha, 22);
00785 double alphaVal = sp_alpha->Eval(E);
00786 return alphaVal;
00787 }
00788
00789 double BesTofDigitizerEcV4::StripStruct::getEta(double E)
00790 {
00791
00792
00793 double e[22] =
00794 {
00795 65,
00796 70 ,
00797 75 ,
00798 80 ,
00799 85 ,
00800 90 ,
00801 95 ,
00802 100 ,
00803 102 ,
00804 104 ,
00805 106 ,
00806 108 ,
00807 110 ,
00808 112 ,
00809 114 ,
00810 116 ,
00811 118 ,
00812 120 ,
00813 125 ,
00814 130 ,
00815 135 ,
00816 140
00817 };
00818
00819
00820 double eta[22]=
00821 {
00822 132.6/10 ,
00823 117.2/10 ,
00824 102.6/10 ,
00825 88.26/10 ,
00826 79.81/10 ,
00827 74.0 /10 ,
00828 66.7 /10 ,
00829 62.7 /10 ,
00830 61.4 /10 ,
00831 57.4 /10 ,
00832 55.45/10 ,
00833 54.35/10 ,
00834 52.48/10 ,
00835 51.3 /10 ,
00836 50.1 /10 ,
00837 48.3 /10 ,
00838 48.28/10 ,
00839 46.00/10 ,
00840 44.08/10 ,
00841 41.67/10 ,
00842 39.97/10 ,
00843 38.04/10
00844 };
00845
00846 TSpline3* sp_eta = new TSpline3("sp_eta", e, eta, 22);
00847 double etaVal = sp_eta->Eval(E);
00848 return etaVal;
00849 }
00850
00851 double BesTofDigitizerEcV4::StripStruct::getV(double E)
00852 {
00853
00854
00855 double e[22] =
00856 {
00857 65,
00858 70 ,
00859 75 ,
00860 80 ,
00861 85 ,
00862 90 ,
00863 95 ,
00864 100 ,
00865 102 ,
00866 104 ,
00867 106 ,
00868 108 ,
00869 110 ,
00870 112 ,
00871 114 ,
00872 116 ,
00873 118 ,
00874 120 ,
00875 125 ,
00876 130 ,
00877 135 ,
00878 140
00879 };
00880
00881
00882 double v[22]=
00883 {
00884 130.2/1000 ,
00885 138.5/1000 ,
00886 146.7/1000 ,
00887 155.0/1000 ,
00888 163.3/1000 ,
00889 171.4/1000 ,
00890 179.7/1000 ,
00891 187.7/1000 ,
00892 191.2/1000 ,
00893 194.5/1000 ,
00894 197.9/1000 ,
00895 201.2/1000 ,
00896 204.5/1000 ,
00897 207.6/1000 ,
00898 210.9/1000 ,
00899 214.4/1000 ,
00900 217.5/1000 ,
00901 220.9/1000 ,
00902 228.8/1000 ,
00903 237.0/1000 ,
00904 244.7/1000 ,
00905 252.9/1000
00906 };
00907
00908 TSpline3* sp_v = new TSpline3("sp_v", e, v, 22);
00909 double vVal = sp_v->Eval(E);
00910 return vVal;
00911 }
00912
00913 void BesTofDigitizerEcV4::Param::print()
00914 {
00915 cout<<"Fixed parameters: "<<endl;
00916 for(int i=0; i<nstrip; i++)
00917 {
00918 cout<<" strip_x["<<i<<"]= "<<strip_x[i];
00919 }
00920 for(int i=0; i<nmodule; i++)
00921 {
00922 for(int j=0; j<nstrip; j++)
00923 {
00924 cout<<" deadChannel["<<i<<"]["<<j<<"]= "<<deadChannel[i][j];
00925 }
00926 }
00927
00928 cout<<" strip_z= "<<strip_z
00929 <<" strip_gap= "<<strip_gap
00930 <<" ngap= "<<ngap
00931 <<" gapWidth= "<<gapWidth
00932 <<" nstep= "<<nstep
00933 <<" stepWidth= "<<stepWidth
00934 <<" E_weight= "<<E_weight
00935 <<" eCharge= "<<eCharge
00936 <<endl;
00937 }
00938
00939 void BesTofDigitizerEcV4::HitStruct::print()
00940 {
00941 cout<<"Hit information: "<<endl;
00942 cout<<" trkIndex= "<<trkIndex
00943 <<" pdgCode= "<<pdgCode
00944 <<" ions= "<<pdgCode
00945 <<" strip= "<<strip
00946 <<" gap= "<<gap
00947 <<" glbTime= "<<glbTime
00948 <<" locx= "<<locx
00949 <<" locy= "<<locy
00950 <<" locz= "<<locz
00951 <<" x= "<<x
00952 <<" y= "<<y
00953 <<" z= "<<z
00954 <<" px= "<<px
00955 <<" py= "<<py
00956 <<" pz= "<<pz
00957 <<" v_propagate= "<<v_propagate
00958 <<" tPropagate_sphi= "<<tPropagate_sphi
00959 <<" tPropagate_xphi= "<<tPropagate_xphi
00960 <<endl;
00961 }
00962
00963 void BesTofDigitizerEcV4::StripStruct::print()
00964 {
00965 cout<<"Strip information: "<<endl;
00966 cout<<" strip= "<<strip
00967 <<" trkIndex= "<<trkIndex
00968 <<" tStart= "<<tStart
00969 <<" tPropagate_sphi= "<<tPropagate_sphi
00970 <<" tPropagate_xphi= "<<tPropagate_xphi
00971 <<" tThreshold "<<tThreshold
00972 <<" charge= "<<charge
00973 <<" E= "<<E
00974 <<" alpha= "<<alpha
00975 <<" eta= "<<eta
00976 <<" threshold= "<<threshold
00977 <<" v_drift= "<<v_drift
00978 <<endl;
00979 }