#include <BesTofDigitizerEcV3.hh>
Inheritance diagram for BesTofDigitizerEcV3:
Public Member Functions | |
| void | AccuSignal (G4double, G4int) |
| void | AccuSignal (G4double, G4int) |
| BesTofDigitizerEcV3 () | |
| BesTofDigitizerEcV3 () | |
| G4double | BirksLaw (BesTofHit *hit) |
| G4double | BirksLaw (BesTofHit *hit) |
| virtual void | Digitize (ScintSingle *, BesTofDigitsCollection *) |
| virtual void | Digitize (ScintSingle *, BesTofDigitsCollection *) |
| void | DirectPh (G4int, G4int, G4int, G4double &) |
| void | DirectPh (G4int, G4int, G4int, G4double &) |
| void | Initialize () |
| void | Initialize () |
| void | ReadData () |
| void | ReadData () |
| void | ReadEffTree () |
| void | ReadEffTree () |
| G4double | Scintillation (G4int) |
| G4double | Scintillation (G4int) |
| void | TofPmtAccum (BesTofHit *) |
| void | TofPmtAccum (BesTofHit *) |
| void | TofPmtInit () |
| void | TofPmtInit () |
| void | TofPmtRspns (G4int, G4int) |
| void | TofPmtRspns (G4int, G4int) |
| G4double | TransitTime () |
| G4double | TransitTime () |
| ~BesTofDigitizerEcV3 () | |
| ~BesTofDigitizerEcV3 () | |
Protected Attributes | |
| G4double | m_ADC [2] |
| BesTofDigitsCollection * | m_besTofDigitsCollection |
| BesTofDigitsCollection * | m_besTofDigitsCollection |
| G4double | m_globalTime |
| G4double | m_TDC [2] |
| BesTofHitsCollection * | m_THC |
| BesTofHitsCollection * | m_THC |
| ITofCaliSvc * | m_tofCaliSvc |
| ITofCaliSvc * | m_tofCaliSvc |
| ITofQElecSvc * | m_tofQElecSvc |
| ITofQElecSvc * | m_tofQElecSvc |
| ITofSimSvc * | m_tofSimSvc |
| ITofSimSvc * | m_tofSimSvc |
| G4int | m_trackIndex |
Static Protected Attributes | |
| NTuple::Item< double > | m_adc0 |
| NTuple::Item< double > | m_adc0 |
| NTuple::Item< double > | m_adc1 |
| NTuple::Item< double > | m_adc1 |
| bool | m_booked = false |
| NTuple::Item< double > | m_ddT |
| NTuple::Item< double > | m_ddT |
| NTuple::Item< double > | m_edep |
| NTuple::Item< double > | m_edep |
| NTuple::Item< double > | m_edepHit |
| NTuple::Item< double > | m_edepHit |
| NTuple::Item< double > | m_edepMPV |
| NTuple::Item< double > | m_edepMPV |
| NTuple::Item< double > | m_endTime |
| NTuple::Item< double > | m_endTime |
| NTuple::Item< double > | m_eTotal |
| NTuple::Item< double > | m_eTotal |
| NTuple::Item< double > | m_forb |
| NTuple::Item< double > | m_forb |
| NTuple::Item< double > | m_max0 |
| NTuple::Item< double > | m_max0 |
| NTuple::Item< double > | m_max1 |
| NTuple::Item< double > | m_max1 |
| NTuple::Item< double > | m_nDigi |
| NTuple::Item< double > | m_nDigi |
| NTuple::Item< double > | m_nDigiOut |
| NTuple::Item< double > | m_nDigiOut |
| NTuple::Item< double > | m_nHits |
| NTuple::Item< double > | m_nHits |
| NTuple::Item< double > | m_NphAllSteps |
| NTuple::Item< double > | m_NphAllSteps |
| NTuple::Item< double > | m_partId |
| NTuple::Item< double > | m_partId |
| NTuple::Item< double > | m_partIdMPV |
| NTuple::Item< double > | m_partIdMPV |
| NTuple::Item< double > | m_scinNb |
| NTuple::Item< double > | m_scinNb |
| NTuple::Item< double > | m_scinNbMPV |
| NTuple::Item< double > | m_scinNbMPV |
| NTuple::Item< double > | m_scinSwim |
| NTuple::Item< double > | m_scinSwim |
| NTuple::Item< double > | m_scinTime |
| NTuple::Item< double > | m_scinTime |
| NTuple::Item< double > | m_tdc0 |
| NTuple::Item< double > | m_tdc0 |
| NTuple::Item< double > | m_tdc1 |
| NTuple::Item< double > | m_tdc1 |
| NTuple::Item< double > | m_time1st0 |
| NTuple::Item< double > | m_time1st0 |
| NTuple::Item< double > | m_time1st1 |
| NTuple::Item< double > | m_time1st1 |
| NTuple::Item< double > | m_timeFlight |
| NTuple::Item< double > | m_timeFlight |
| NTuple::Item< double > | m_timelast0 |
| NTuple::Item< double > | m_timelast0 |
| NTuple::Item< double > | m_timelast1 |
| NTuple::Item< double > | m_timelast1 |
| NTuple::Item< double > | m_totalPhot0 |
| NTuple::Item< double > | m_totalPhot0 |
| NTuple::Item< double > | m_totalPhot1 |
| NTuple::Item< double > | m_totalPhot1 |
| NTuple::Item< double > | m_transitTime |
| NTuple::Item< double > | m_transitTime |
| NTuple::Tuple * | m_tupleTof1 |
| NTuple::Tuple * | m_tupleTof1 = 0 |
| NTuple::Tuple * | m_tupleTof2 |
| NTuple::Tuple * | m_tupleTof2 = 0 |
| NTuple::Tuple * | m_tupleTof3 |
| NTuple::Tuple * | m_tupleTof3 = 0 |
Private Attributes | |
| float | eff [50][10][10] |
| G4double | m_attenEc |
| G4double | m_beamTime |
| G4double | m_CeEc |
| G4double | m_CEEc |
| G4double | m_Cpe2pmtEc |
| G4double | m_ecR1 |
| G4Svc * | m_G4Svc |
| G4Svc * | m_G4Svc |
| G4double | m_HLthreshEc |
| G4double | m_LLthreshEc |
| G4double | m_noiseSigmaEc |
| G4int | m_nPhot [m_profBinNEcV3][2] |
| G4double | m_peCorFacEc |
| G4double | m_phNConstEc |
| G4double | m_PMTgainEc |
| G4double | m_preGainEc |
| G4double | m_QEEc |
| G4double | m_rAngleEc |
| RealizationSvc * | m_RealizationSvc |
| RealizationSvc * | m_RealizationSvc |
| G4double | m_refIndexEc |
| G4double | m_riseTimeEc |
| G4double | m_t1st [2] |
| G4double | m_tau1Ec |
| G4double | m_tau2Ec |
| G4double | m_tau3Ec |
| G4double | m_tauRatioEc |
| G4double | m_timeBinSize |
| G4double | m_tLast [2] |
| G4int | m_totalPhot [2] |
| G4double | m_ttsMeanEc |
| G4double | m_ttsSigmaEc |
| float | prob [50][10][10][num1] |
| float | propTime [50][10][10][num1] |
|
|
00032 {
00033 ReadData();
00034 m_timeBinSize = 0.005;
00035
00036
00037
00038 //retrieve G4Svc
00039 ISvcLocator* svcLocator = Gaudi::svcLocator();
00040 IG4Svc* tmpSvc;
00041 StatusCode sc = svcLocator->service("G4Svc", tmpSvc);
00042 if(!sc.isSuccess())
00043 {
00044 std::cout << " Could not initialize Realization Service in BesTofDigitizerBrV2" << std::endl;
00045 }
00046 else
00047 {
00048 m_G4Svc = dynamic_cast<G4Svc *>(tmpSvc);
00049 }
00050
00051 //retrieve RealizationSvc
00052 IRealizationSvc *tmpReal;
00053 StatusCode scReal = svcLocator->service("RealizationSvc",tmpReal);
00054 if (!scReal.isSuccess())
00055 {
00056 std::cout << " Could not initialize Realization Service in BesTofDigitizerEcV3" << std::endl;
00057 }
00058 else
00059 {
00060 m_RealizationSvc = dynamic_cast<RealizationSvc*>(tmpReal);
00061 }
00062
00063
00064
00065 for(int i=0;i<50;i++)
00066 {
00067 for(int j=0;j<10;j++)
00068 {
00069 for(int k=0;k<10;k++)
00070 {
00071 for(int m=0;m<num1;m++)
00072 {
00073 //G4cout << "time:" << propTime[i][j][k][m] << "; prob:" << prob[i][j][k][m] << "; eff:" << eff[i][j][k] << G4endl;
00074 propTime[i][j][k][m] = 0;
00075 prob[i][j][k][m] = 0;
00076 eff[i][j][k] = 0;
00077 }
00078 }
00079 }
00080 }
00081
00082 ReadEffTree();
00083 G4cout << "ETofSim: Reading nTuples of is completed." << G4endl;
00084 }
|
|
|
00173 {;}
|
|
|
|
|
|
|
|
||||||||||||
|
|
|
||||||||||||
|
00550 {
00551 G4int ihst;
00552 ihst=G4int(endTime/m_timeBinSize);
00553 if (ihst>0 &&ihst<m_profBinNEcV3)
00554 {
00555 m_nPhot[ihst][forb]=m_nPhot[ihst][forb]+1;
00556 m_totalPhot[forb]=m_totalPhot[forb]+1;
00557 }
00558 }
|
|
|
|
|
|
00470 {
00471 const G4double kappa = 0.015*cm/MeV;
00472 const G4String brMaterial = "BC404";
00473 G4double dE = hit->GetEdep();
00474 //G4cout << "The edep is "<< dE << G4endl;
00475 G4double dX = hit->GetStepL();
00476 //G4Material* materiral = hit->GetMaterial();
00477 G4double charge = hit->GetCharge();
00478 G4double cor_dE = dE;
00479 //if((materiral->GetName()==brMaterial) && charge!=0.&& dX!=0.)
00480 if(charge!=0.&& dX!=0.)
00481 {
00482 cor_dE = dE/(1+kappa*dE/dX);
00483 //if(dE>20)
00484 //{
00485 // G4cout << "\n dE > 20. Details are below:" << G4endl;
00486 // G4cout << "dE/dx:" << dE/dX << G4endl;
00487 // G4cout << "dE:" << dE << "; dX:" << dX << G4endl;
00488 // G4cout << "It is BC404. cor_dE is " << cor_dE << G4endl;
00489 // G4double ratio = cor_dE/dE;
00490 // G4cout << "The ratio cor_dE/edep is "<< ratio << G4endl;
00491 //}
00492 //G4cout << "It is BC408. cor_dE is " << cor_dE << G4endl;
00493 //G4double ratio = cor_dE/dE;
00494 //G4cout << "The ratio cor_dE/edep is "<< ratio << G4endl;
00495 }
00496 return cor_dE;
00497
00498 }
|
|
||||||||||||
|
Reimplemented from BesTofDigitizerV. |
|
||||||||||||
|
Reimplemented from BesTofDigitizerV. 00176 {
00177 m_beamTime = m_G4Svc->GetBeamTime() * ns;
00178 m_besTofDigitsCollection = DC;
00179
00180 G4DigiManager* digiManager = G4DigiManager::GetDMpointer();
00181
00182 G4int THCID = digiManager->GetHitsCollectionID("BesTofHitsCollection");
00183 m_THC = (BesTofHitsCollection*) (digiManager->GetHitsCollection(THCID));
00184
00185 if (m_G4Svc->TofRootFlag())
00186 {
00187 m_eTotal = 0;
00188 m_nDigi = 0;
00189 m_partIdMPV = -9;
00190 m_scinNbMPV = -9;
00191 m_edepMPV = 0;
00192 m_nDigiOut = 0;
00193 }
00194
00195 if (m_THC)
00196 {
00197 //for each digi, compute TDC and ADC
00198 G4int partId, scinNb, nHits;
00199 G4double edep;
00200 BesTofHit* hit;
00201 partId=scint->GetPartId();
00202 scinNb=scint->GetScinNb();
00203 edep = scint->GetEdep();
00204 nHits=scint->GetHitIndexes()->size();
00205
00206 TofPmtInit();
00207
00208 //fill tof Ntuple
00209 if (m_G4Svc->TofRootFlag())
00210 {
00211 if (edep>m_edepMPV)
00212 {
00213 m_partIdMPV = partId;
00214 m_scinNbMPV = scinNb;
00215 m_edepMPV = edep;
00216 }
00217 m_eTotal += edep;
00218 m_nDigi ++;
00219
00220 m_partId = partId;
00221 m_scinNb = scinNb;
00222 m_edep = edep;
00223 m_nHits = nHits;
00224 }
00225
00226 if (edep>0.01)
00227 {
00228 for (G4int j=0;j<nHits;j++)
00229 {
00230 hit= (*m_THC)[( *(scint->GetHitIndexes()) )[j]];
00231 TofPmtAccum(hit);
00232 }
00233
00234 if (m_G4Svc->TofRootFlag())
00235 {
00236 m_time1st0=m_t1st[0];
00237 m_time1st1=m_t1st[1];
00238 m_timelast0=m_tLast[0];
00239 m_timelast1=m_tLast[1];
00240 m_totalPhot0=m_totalPhot[0];
00241 m_totalPhot1=m_totalPhot[1];
00242 }
00243
00244 //get final tdc and adc
00245 TofPmtRspns(partId,scinNb);
00246
00247 G4double temp0 = m_ADC[0]+m_TDC[0];
00248 G4double temp1 = m_ADC[1]+m_TDC[1];
00249 //const double MAX_ADC = 8191*0.3; // channel set up to 8192 will lead to overflow.
00250 if ( (partId!=1) && temp0>0. )
00251 {
00252 BesTofDigi* digi = new BesTofDigi;
00253 digi->SetTrackIndex(m_trackIndex);
00254 digi->SetPartId(partId);
00255 digi->SetScinNb(scinNb);
00256 digi->SetForwADC( m_ADC[0]) ;
00257 digi->SetBackADC( m_ADC[1]) ;
00258 if (m_TDC[0]>0.)
00259 m_TDC[0] = m_TDC[0]+m_beamTime;
00260 digi->SetForwTDC( m_TDC[0]) ;
00261 digi->SetBackTDC( m_TDC[1]) ;
00262 //G4cout<<"endcap\nadc0:"<<m_ADC[0]<<"; adc1:"<<m_ADC[1]<<"; tdc0:"<<m_TDC[0]<<"; tdc1:"<<m_TDC[1]<<G4endl;
00263 m_besTofDigitsCollection->insert(digi);
00264 if (m_G4Svc->TofRootFlag() )
00265 m_nDigiOut++;
00266 }
00267 if (m_G4Svc->TofRootFlag() )
00268 m_tupleTof1->write();
00269 //cout << "m_tupleTof1->write()" << endl;
00270 }
00271
00272 }
00273 if (m_G4Svc->TofRootFlag())
00274 m_tupleTof2->write();
00275 //cout << "m_tupleTof2->write()" << endl;
00276
00277 }
|
|
||||||||||||||||||||
|
|
|
||||||||||||||||||||
|
00501 {
00502 G4double ran = G4UniformRand();
00503 G4double p = 0;
00504 G4int nth = 1;
00505 G4int key = 0;
00506 t = 0;
00507 while(1)
00508 {
00509 if(p>ran||nth==400)
00510 {
00511 key = nth;
00512 //G4cout << "Value found!" << G4endl;
00513 break;
00514 }
00515 p = p + prob[rBin][phiBin][zBin][nth];
00516 nth++;
00517 }
00518 t = propTime[rBin][phiBin][zBin][key-1];
00519 }
|
|
|
|
|
|
00149 {
00150 for (G4int i=0;i<2;i++)
00151 {
00152 m_ADC[i] = -999;
00153 m_TDC[i] = -999;
00154 }
00155 m_trackIndex = -999;
00156 m_globalTime = 9999;
00157 }
|
|
|
|
|
|
00087 {
00088 BesTofGeoParameter* tofPara = BesTofGeoParameter::GetInstance();
00089
00090 m_ecR1 = tofPara->GetEcR1();
00091 m_tau1Ec = tofPara->GetTau1Ec();
00092 m_tau2Ec = tofPara->GetTau2Ec();
00093 m_tau3Ec = tofPara->GetTau3Ec();
00094 m_tauRatioEc = tofPara->GetTauRatioEc();
00095 m_refIndexEc = tofPara->GetRefIndexEc();
00096 m_phNConstEc = tofPara->GetPhNConstEc();
00097 m_Cpe2pmtEc = tofPara->GetCpe2pmtEc();
00098 m_rAngleEc = tofPara->GetRAngleEc();
00099 m_QEEc = tofPara->GetQEEc();
00100 m_CEEc = tofPara->GetCEEc();
00101 m_peCorFacEc = tofPara->GetPeCorFacEc();
00102 m_attenEc = tofPara->GetAttenEc();
00103
00104 m_ttsMeanEc = tofPara->GetTTSmeanEc();
00105 m_ttsSigmaEc = tofPara->GetTTSsigmaEc();
00106 m_PMTgainEc = tofPara->GetPMTgainEc();
00107 m_CeEc = tofPara->GetCeEc();
00108 m_riseTimeEc = tofPara->GetRiseTimeEc();
00109 m_LLthreshEc = tofPara->GetLLthreshEc();
00110 m_HLthreshEc = tofPara->GetHLthreshEc();
00111 m_preGainEc = tofPara->GetPreGainEc();
00112 m_noiseSigmaEc = tofPara->GetNoiseSigmaEc();
00113
00114 //G4cout << "m_LLthreshEc:" << m_LLthreshEc << "; m_HLthreshEc:" << m_HLthreshEc << G4endl;
00115
00116 }
|
|
|
|
|
|
00119 {
00120 int rBin,phiBin,zBin;
00121 const int nR = 43;
00122 const int nPhi = 6;
00123 const int nZ = 6;
00124 float efficiency0,x[400],y[400];
00125
00126
00127 G4String dataPath = getenv("TOFSIMROOT");
00128 if(!dataPath)
00129 {
00130 G4Exception("Boss environment is not set!");
00131 }
00132
00133 char treePath[200];
00134 G4int runId = m_RealizationSvc->getRunId();
00135 if(runId>=-80000 && runId<=-9484)
00136 {
00137 // After TOF HV adjustment, endcap attenL was set to 5000mm.
00138 sprintf(treePath,"%s/dat/effTree_1600mm.root",dataPath.c_str());
00139 }
00140 else
00141 {
00142 //Before TOF HV adjustment, endcap attenL was set to 1600mm.
00143 sprintf(treePath,"%s/dat/effTree_1600mm.root",dataPath.c_str());
00144 }
00145
00146 TFile *f = new TFile(treePath, "read");
00147 TTree *t = (TTree*)f->Get("effTree");
00148
00149 t->SetBranchAddress("rBin", &rBin);
00150 t->SetBranchAddress("phiBin", &phiBin);
00151 t->SetBranchAddress("zBin", &zBin);
00152 t->SetBranchAddress("efficiency0", &efficiency0);
00153 t->SetBranchAddress("x", x);
00154 t->SetBranchAddress("y", y);
00155
00156 int r,phi,z;
00157 for (Int_t i = 0; i < nR*nPhi*nZ; i++){
00158 t->GetEntry(i);
00159 r = rBin;
00160 phi = phiBin;
00161 z = zBin;
00162 eff[r][phi][z] = efficiency0;
00163 for (Int_t j = 0; j < 400; j++){
00164 propTime[r][phi][z][j] = x[j];
00165 prob[r][phi][z][j] = y[j];
00166 //cout << "\n" << propTime[r][phi][z][j] << " " << prob[r][phi][z][j] << endl;
00167 }
00168 }
00169
00170 }
|
|
|
|
|
|
00522 {
00523 G4double tmp_tauRatio,tmp_tau1,tmp_tau2,tmp_tau3;
00524 tmp_tauRatio = m_tauRatioEc;
00525 tmp_tau1 = m_tau1Ec;
00526 tmp_tau2 = m_tau2Ec;
00527 tmp_tau3 = m_tau3Ec;
00528
00529 G4double UniformR = tmp_tauRatio/(1+tmp_tauRatio);
00530 G4double EmissionTime;
00531 if (G4UniformRand()>UniformR) {
00532 while (1) {
00533 EmissionTime = -tmp_tau2*log( G4UniformRand() );
00534 if (G4UniformRand()-exp(EmissionTime/tmp_tau2-EmissionTime/tmp_tau1)>1.E-8)
00535 break;
00536 }
00537 }
00538 else EmissionTime = -tmp_tau3*log( G4UniformRand() );
00539 return EmissionTime;
00540 }
|
|
|
|
|
|
00322 {
00323 G4double cvelScint = c_light/m_refIndexEc;
00324 //Get information of this step
00325 G4ThreeVector pos = hit->GetPos();
00326 G4int trackIndex = hit->GetTrackIndex();
00327 G4int partId = hit->GetPartId();
00328 G4double edep = hit->GetEdep();
00329 G4double stepL = hit->GetStepL();
00330 //G4String particleName = hit->GetPName();
00331 G4double deltaT=hit->GetDeltaT();
00332 G4double timeFlight=hit->GetTime()-m_beamTime;
00333 if (timeFlight < m_globalTime)
00334 {
00335 m_globalTime = timeFlight;
00336 m_trackIndex = trackIndex;
00337 }
00338
00339 G4ThreeVector pDirection=hit->GetPDirection();
00340 G4double nx=pDirection.x();
00341 G4double ny=pDirection.y();
00342 G4double nz=pDirection.z();
00343
00344
00345 //number of photons generated in this step
00346 G4int NphStep;
00347 G4double nMean, nPhoton;
00348 nMean = m_phNConstEc*BirksLaw(hit);
00349
00350 if(nMean>10)
00351 {
00352 G4double resolutionScale=1.;
00353 G4double sigma=resolutionScale*sqrt(nMean);
00354 nPhoton=G4int(G4RandGauss::shoot(nMean,sigma)+0.5);
00355 }
00356 else
00357 nPhoton=G4int(G4Poisson(nMean));
00358
00359
00360 NphStep=G4int(nPhoton*0.66*m_QEEc*m_CEEc);
00361 //G4cout << "\nradius:" << radius << "; phi:" << phi << "; z:" << z << G4endl;
00362 //G4cout << "\nrBin:" << rBin << ";phiBin:" << phiBin << ";zBin" << zBin << G4endl;
00363
00364 if (m_G4Svc->TofRootFlag())
00365 m_NphAllSteps += G4int(nPhoton*0.66*m_QEEc*m_CEEc);
00366
00367 if (NphStep>0)
00368 {
00369 for (G4int i=0;i<NphStep;i++)
00370 {
00371 //uniform scintilation in each step
00372 G4double ddS, ddT;
00373 ddS=stepL*G4UniformRand();
00374 ddT=deltaT*G4UniformRand();
00375 G4ThreeVector emtPos;
00376 emtPos.setX(pos.x() + nx*ddS);
00377 emtPos.setY(pos.y() + ny*ddS);
00378 emtPos.setZ(pos.z() + nz*ddS);
00379
00380 //retrieve the histogram info
00381 G4double radius = sqrt(emtPos.x()*emtPos.x()+emtPos.y()*emtPos.y())-m_ecR1;
00382 const G4double pie = 2.*asin(1.);
00383 G4double phi;
00384 if(emtPos.x()>0 && emtPos.y()>0)
00385 phi = atan(emtPos.y()/emtPos.x());
00386 else if(emtPos.x()==0 && emtPos.y()>0)
00387 phi = pie/2.;
00388 else if(emtPos.x()<0)
00389 phi = atan(emtPos.y()/emtPos.x())+pie;
00390 else if(emtPos.x()==0 && emtPos.y()<0)
00391 phi = 1.5*pie;
00392 else if(emtPos.x()>0 && emtPos.y()<0)
00393 phi = 2.*pie+atan(emtPos.y()/emtPos.x());
00394 phi = phi*180./pie; // in degrees
00395 G4double z = fabs(emtPos.z());
00396 // Warning: Should obtain absolute value of z to determine zBinNum
00397
00398 G4int rBin = G4int(radius/10.); // radius bin no.
00399 G4double resPhi = phi-(G4int(phi/7.5)*7.5); // residual of phi in period of 7.5
00400 G4int phiBin = G4int(resPhi/1.25);
00401 G4int zBin = G4int((z-1332.)/8.);
00402
00403 //check scinillation light whether to hit the pmt or not
00404 //forb=0/1 for forward(z>0, east) and backward(z<0, west)
00405 G4int forb = 0;
00406 G4double transpTime = 0;
00407 G4double pathL = 0;
00408 G4double efficiency1;
00409 G4double efficiency2;
00410 efficiency1 = G4RandGauss::shoot(0,0.004);
00411 if(rBin>=0&&rBin<=nR && phiBin>=0&& phiBin<=nPhi && zBin>=0&&zBin<=nZ)
00412 efficiency1 += eff[rBin][phiBin][zBin];
00413 else
00414 efficiency1 = 0;
00415 //G4cout << "FATAL: The collection efficiency does NOT exist!" << G4endl;
00416 if(m_attenEc==0)
00417 {
00418 G4cout <<" ERROR: Attenuation Length is null!" << G4endl;
00419 break;
00420 }
00421 //efficiency2 = pow(efficiency1,(1600/m_attenEc));
00422 if(G4UniformRand() <= efficiency1)
00423 {
00424 DirectPh(rBin, phiBin, zBin, transpTime);
00425 //cout << "transpTime:" << transpTime << endl;
00426 }
00427
00428 //check if photon can reach PMT or not, after attenuation
00429 //G4double ran = G4UniformRand();
00430 //pathL = transpTime*cvelScint;
00431 //if (pathL>0 && exp(-pathL/m_attenEc) > ran) // Note: Do NOT double count attuation!
00432 if(transpTime>0)
00433 {
00434 //propagation time in scintillator
00435 G4double scinSwim = transpTime;
00436 //scintillation timing
00437 G4double scinTime = Scintillation(partId);
00438
00439 //PMT transit time
00440 G4double transitTime = TransitTime();
00441 //sum up all time components
00442 G4double endTime = timeFlight + ddT + scinSwim + scinTime + transitTime;
00443
00444 if (m_G4Svc->TofRootFlag())
00445 {
00446 //m_forb = forb;
00447 m_timeFlight = timeFlight+ddT;
00448 m_ddT = ddT;
00449 m_scinSwim = scinSwim;
00450 m_scinTime = scinTime;
00451 m_transitTime = transitTime;
00452 m_endTime = endTime;
00453 m_tupleTof3->write();
00454 }
00455
00456 //store timing into binned buffer
00457 AccuSignal(endTime, forb);
00458
00459 //update 1st and last timings here
00460 if (m_t1st[forb]>endTime) m_t1st[forb] = endTime;
00461 if (m_tLast[forb]<endTime) m_tLast[forb]= endTime;
00462 //if(m_tLast[0]>100)
00463 //std::cout<<"endTime: "<<endTime<<std::endl;
00464 }
00465 }
00466 }
00467 }
|
|
|
|
|
|
00280 {
00281 m_ADC[0] = -999.;
00282 m_ADC[1] = -999.;
00283 m_TDC[0] = -999.;
00284 m_TDC[1] = -999.;
00285 m_trackIndex = -999;
00286 m_globalTime = 9999;
00287
00288 m_t1st[0]=100;
00289 m_t1st[1]=100;
00290 m_tLast[0]=0.;
00291 m_tLast[1]=0;
00292 m_totalPhot[0]=0;
00293 m_totalPhot[1]=0;
00294 for (G4int i=0;i<2;i++)
00295 for (G4int j=0;j<m_profBinNEcV3;j++)
00296 m_nPhot[j][i]=0;
00297
00298 if (m_G4Svc->TofRootFlag())
00299 {
00300 m_partId = -9;
00301 m_scinNb = -9;
00302 m_edep = 0;
00303 m_nHits = 0;
00304 m_time1st0 = 100;
00305 m_time1st1 = 100;
00306 m_timelast0 = 0;
00307 m_timelast1 = 0;
00308 m_totalPhot0 = 0;
00309 m_totalPhot1 = 0;
00310 m_NphAllSteps = 0;
00311 m_max0 = 0;
00312 m_max1 = 0;
00313 m_tdc0 = -999;
00314 m_adc0 = -999;
00315 m_tdc1 = -999;
00316 m_adc1 = -999;
00317 }
00318
00319 }
|
|
||||||||||||
|
|
|
||||||||||||
|
00561 {
00562 //to generate PMT signal shape for single photoelectron.
00563 //only one time for a job.
00564 static G4double snpe[m_snpeBinNEcV3];
00565 static G4int istore_snpe=-1;
00566
00567 //Model: f(t)=Gain*mv_1pe* t**2 * exp-(t/tau)**2/normal-const
00568 //normalization const =sqrt(pi)*tau*tau*tau/4.0
00569 G4double tau = m_riseTimeEc;
00570 G4double norma_const=sqrt(M_PI)*tau*tau*tau/4.0; //in unit of ns**3
00571 G4double echarge=1.6e-7; //in unit of PC
00572
00573 //time profile of pmt signals for Back and Forward PMT.
00574 G4double profPmt[m_profBinNEcV3][2];
00575
00576 G4double t;
00577 G4int n1, n2, ii;
00578 G4int phtn;
00579
00580 if (istore_snpe<0)
00581 {
00582 istore_snpe = 1;
00583 for (G4int i=0;i<m_snpeBinNEcV3;i++)
00584 {
00585 t=(i+1)*m_timeBinSize;
00586 snpe[i]=m_CeEc*t*t*exp(- (t/tau) * (t/tau) )/norma_const;
00587 }
00588 }
00589 //for barrel and endcap tof: fb=2 or 1
00590 G4int fb=1; // for endcap part
00591
00592 G4double tmpADC[2] = {0,0};
00593
00594 for (G4int j=0; j<fb; j++)
00595 {
00596 if (m_totalPhot[j] > 0)
00597 {
00598 n1=G4int(m_t1st[j]/m_timeBinSize);
00599 n2=G4int(m_tLast[j]/m_timeBinSize);
00600
00601 for (G4int i=0;i<m_profBinNEcV3;i++)
00602 profPmt[i][j]=0.0;
00603
00604 //generate PMT pulse
00605 n2 = n2<m_profBinNEcV3 ? n2:m_profBinNEcV3;
00606 for (G4int i=n1;i<n2;i++)
00607 {
00608 phtn=m_nPhot[i][j];
00609 if (phtn>0)
00610 {
00611 G4double Npoisson;
00612 while(1) {
00613 Npoisson=G4Poisson(10.0);
00614 if(Npoisson>0) break;
00615 }
00616 G4double tmpPMTgain;
00617 while(1) {
00618 m_PMTgainEc = m_tofSimSvc->EndPMTGain();
00619 tmpPMTgain=G4RandGauss::shoot(m_PMTgainEc,m_PMTgainEc/sqrt(Npoisson));
00620 //tmpPMTgain = m_PMTgainEc;
00621 if(tmpPMTgain>0) break;
00622 }
00623 tmpADC[j]+=phtn*tmpPMTgain;
00624
00625 for (G4int ihst=0; ihst<m_snpeBinNEcV3; ihst++)
00626 {
00627 ii=i+ihst;
00628 if (ii<m_profBinNEcV3)
00629 profPmt[ii][j] += tmpPMTgain*phtn*snpe[ihst];
00630 else
00631 break;
00632 }
00633 }
00634 }
00635
00636 //add preamplifier and noise
00637 for (int i=0;i<m_profBinNEcV3;i++)
00638 {
00639 if (profPmt[i][j]>0)
00640 profPmt[i][j] = m_preGainEc*profPmt[i][j]+G4RandGauss::shoot(0,m_noiseSigmaEc);
00641 }
00642
00643 //get pulse height
00644 G4double max=0;
00645 for (int i=n1;i<m_profBinNEcV3;i++)
00646 {
00647 if (profPmt[i][j]>max)
00648 max=profPmt[i][j];
00649 }
00650 if (m_G4Svc->TofRootFlag())
00651 {
00652 if (j==0) m_max0=max;
00653 else m_max1=max;
00654 }
00655
00656
00657 G4double tmp_HLthresh, tmp_LLthresh, ratio;
00658
00659 // G4int runId = m_RealizationSvc->getRunId();
00660 // if(runId>=-80000 && runId<=-9484)
00661 // {
00662 // // After TOF HV adjustment, PMT Gain of barrel TOF in MC is set to 5.E5
00663 // // High/Low Threshold for barrel: 500/125 mV
00664 // tmp_HLthresh = m_HLthreshEc;
00665 // tmp_LLthresh = m_LLthreshEc;
00666 // }
00667 // else
00668 // {
00669 // // Before TOF HV adjustment, PMT Gain of barrel TOF in MC is set to 2.5E5
00670 // // High/Low Threshold for barrel: 600/150 mV
00671 // tmp_HLthresh = 800;
00672 // tmp_LLthresh = 200;
00673 // }
00674
00675 G4double adcFactor = 3.35;
00676 // double ratio;
00677 // if (runId>=-80000 && runId<=-9484)
00678 // {
00679 // ratio = 1.22*1354.4/1282.1;
00680 // }
00681 // else
00682 // {
00683 // ratio = 1.77*2028.8/1931.4;
00684 // }
00685
00686 tmp_HLthresh = m_tofSimSvc->EndHighThres();
00687 tmp_LLthresh = m_tofSimSvc->EndLowThres();
00688 ratio = m_tofSimSvc->EndConstant();
00689
00690 //get final tdc and adc
00691 if (max>=tmp_HLthresh)
00692 {
00693 for (int i=0;i<m_profBinNEcV3;i++)
00694 {
00695 if ( profPmt[i][j] >= tmp_LLthresh )
00696 {
00697 m_TDC[j] = i*m_timeBinSize + G4RandGauss::shoot(0,0.025); // Adding Electronical Uncertainty of 25ps
00698 G4double NoiseSigma;
00699 G4int EndNoiseSwitch = int(m_tofSimSvc->EndNoiseSwitch());
00700 //std::cout << " EndNoiseSwitch = " << EndNoiseSwitch << std::endl;
00701 switch (EndNoiseSwitch) {
00702 case 0:
00703 NoiseSigma = 0.;
00704 break;
00705 case 1:
00706 if (partId==0) { NoiseSigma = m_tofSimSvc->EndNoiseSmear(scinNb); }
00707 if (partId==2) { NoiseSigma = 0.; }
00708 break;
00709 case 2:
00710 if (partId==0) { NoiseSigma = m_tofSimSvc->EndNoiseSmear(scinNb); }
00711 if (partId==2) { NoiseSigma = m_tofSimSvc->EndNoiseSmear(scinNb+48.); }
00712 break;
00713 }
00714 //std::cout << "ID : " << scinNb + 48.*partId/2. << " Sigma : " << NoiseSigma << std::endl;
00715 m_TDC[j] = m_TDC[j] + G4RandGauss::shoot(0,NoiseSigma); // Adding Endcap Noise Smear
00716
00717 //use the saturation curve
00718
00719 if( m_G4Svc->TofSaturationFlag())
00720 {
00721 double x = tmpADC[j]*m_preGainEc*echarge*ratio;
00722 int id;
00723 if (partId==0) { id = scinNb;}
00724 if (partId==2) { id = scinNb+48;}
00725
00726 m_ADC[j] = m_tofQElecSvc->EQChannel(id,x);
00727 }
00728 else
00729 m_ADC[j] = tmpADC[j]*m_preGainEc*echarge*adcFactor;
00730
00731 if (m_G4Svc->TofRootFlag())
00732 {
00733 if (j==0) {
00734 m_tdc0 = m_TDC[0];
00735 m_adc0 = m_ADC[0];
00736 }
00737 else {
00738 m_tdc1 = m_TDC[1];
00739 m_adc1 = m_ADC[1];
00740 }
00741 }
00742 break;
00743 }
00744 }
00745 }
00746 }
00747 }
00748 }
|
|
|
|
|
|
00543 {
00544 //get time transit spread
00545 //G4cout << "TTS mean:" << m_ttsMeanEc << "; " << m_ttsSigmaEc << G4endl;
00546 return G4RandGauss::shoot(m_ttsMeanEc,m_ttsSigmaEc);
00547 }
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Reimplemented from BesTofDigitizerV. |
|
|
Reimplemented from BesTofDigitizerV. |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
1.3.9.1