#include <BesTofDigitizerEcV2.hh>
Inheritance diagram for BesTofDigitizerEcV2:
Public Member Functions | |
void | AccuSignal (G4double, G4int) |
void | AccuSignal (G4double, G4int) |
BesTofDigitizerEcV2 () | |
BesTofDigitizerEcV2 () | |
virtual void | Digitize (ScintSingle *, BesTofDigitsCollection *) |
virtual void | Digitize (ScintSingle *, BesTofDigitsCollection *) |
void | DirectPh (G4int, G4ThreeVector, G4double &, G4int &) |
void | DirectPh (G4int, G4ThreeVector, G4double &, G4int &) |
void | Initialize () |
void | Initialize () |
void | ReadData () |
void | ReadData () |
G4double | Scintillation (G4int) |
G4double | Scintillation (G4int) |
void | TofPmtAccum (BesTofHit *) |
void | TofPmtAccum (BesTofHit *) |
void | TofPmtInit () |
void | TofPmtInit () |
void | TofPmtRspns (G4int) |
void | TofPmtRspns (G4int) |
G4double | TransitTime () |
G4double | TransitTime () |
~BesTofDigitizerEcV2 () | |
~BesTofDigitizerEcV2 () | |
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 | |
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_profBinNEcV2][2] |
G4double | m_peCorFacEc |
G4double | m_phNConstEc |
G4double | m_PMTgainEc |
G4double | m_preGainEc |
G4double | m_QEEc |
G4double | m_rAngleEc |
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 |
|
00027 { 00028 ReadData(); 00029 m_timeBinSize=0.005; 00030 00031 //retrieve G4Svc 00032 ISvcLocator* svcLocator = Gaudi::svcLocator(); 00033 IG4Svc* tmpSvc; 00034 StatusCode sc = svcLocator->service("G4Svc", tmpSvc); 00035 m_G4Svc = dynamic_cast<G4Svc *>(tmpSvc); 00036 00037 }
|
|
00069 { 00070 ; 00071 }
|
|
|
|
|
|
|
|
00301 { 00302 G4int ihst; 00303 ihst=G4int(endTime/m_timeBinSize); 00304 if (ihst>0 &&ihst<m_profBinNEcV2) 00305 { 00306 m_nPhot[ihst][forb]=m_nPhot[ihst][forb]+1; 00307 m_totalPhot[forb]=m_totalPhot[forb]+1; 00308 } 00309 }
|
|
Reimplemented from BesTofDigitizerV. |
|
Reimplemented from BesTofDigitizerV. 00074 { 00075 m_beamTime = m_G4Svc->GetBeamTime() * ns; 00076 m_besTofDigitsCollection = DC; 00077 00078 G4DigiManager* digiManager = G4DigiManager::GetDMpointer(); 00079 00080 G4int THCID = digiManager->GetHitsCollectionID("BesTofHitsCollection"); 00081 m_THC = (BesTofHitsCollection*) (digiManager->GetHitsCollection(THCID)); 00082 00083 if (m_THC) 00084 { 00085 //for each digi, compute TDC and ADC 00086 G4int partId, scinNb, nHits; 00087 G4double edep; 00088 BesTofHit* hit; 00089 partId=scint->GetPartId(); 00090 scinNb=scint->GetScinNb(); 00091 edep = scint->GetEdep(); 00092 nHits=scint->GetHitIndexes()->size(); 00093 00094 TofPmtInit(); 00095 00096 if (edep>0.01) 00097 { 00098 for (G4int j=0;j<nHits;j++) 00099 { 00100 hit= (*m_THC)[( *(scint->GetHitIndexes()) )[j]]; 00101 TofPmtAccum(hit); 00102 } 00103 00104 //get final tdc and adc 00105 TofPmtRspns(partId); 00106 00107 G4double temp0 = m_ADC[0]+m_TDC[0]; 00108 G4double temp1 = m_ADC[1]+m_TDC[1]; 00109 if ( (partId==1&&temp0>0&&temp1>0) || ((partId!=1)&&temp0>0) ) 00110 { 00111 if (m_ADC[0]>1000) m_ADC[0]=1000; 00112 if (m_ADC[1]>1000) m_ADC[1]=1000; 00113 BesTofDigi* digi = new BesTofDigi; 00114 digi->SetTrackIndex(m_trackIndex); 00115 digi->SetPartId(partId); 00116 digi->SetScinNb(scinNb); 00117 digi->SetForwADC( m_ADC[0]) ; 00118 digi->SetBackADC( m_ADC[1]) ; 00119 if (m_TDC[0]!=-999) 00120 m_TDC[0] = m_TDC[0]+m_beamTime; 00121 if (m_TDC[1]!=-999) 00122 m_TDC[1] = m_TDC[1]+m_beamTime; 00123 digi->SetForwTDC( m_TDC[0]) ; 00124 digi->SetBackTDC( m_TDC[1]) ; 00125 m_besTofDigitsCollection->insert(digi); 00126 } 00127 } 00128 } 00129 }
|
|
|
|
00223 { 00224 //generation photon have random direction 00225 //optical parameters of scintillation simulation 00226 //G4double cos_spanEc = 1-1/m_refIndex; 00227 G4double cos_spanEc = 1; 00228 G4double ran; 00229 //dcos=cos_span*(ran*2.0-1.0); 00230 ran=G4UniformRand(); 00231 //assuming uniform phi distribution, simulate cos distr only 00232 G4double dcosEc = cos_spanEc*(ran*2.0-1.0); 00233 G4double r2=sqrt(emtPos.x()*emtPos.x()+emtPos.y()*emtPos.y()); 00234 if (dcosEc > 0.0 && dcosEc != 1) 00235 { 00236 forb=0; //forward 00237 pathL = (r2-m_ecR1)/(1.0-dcosEc); 00238 } 00239 else if (dcosEc < 0 && dcosEc != -1) 00240 { 00241 forb=0; 00242 pathL=(2*838-r2-m_ecR1)/(1.0+dcosEc); 00243 } 00244 }
|
|
|
|
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 }
|
|
|
|
00040 { 00041 BesTofGeoParameter* tofPara = BesTofGeoParameter::GetInstance(); 00042 00043 m_ecR1 = tofPara->GetEcR1(); 00044 m_tau1Ec = tofPara->GetTau1Ec(); 00045 m_tau2Ec = tofPara->GetTau2Ec(); 00046 m_tau3Ec = tofPara->GetTau3Ec(); 00047 m_tauRatioEc = tofPara->GetTauRatioEc(); 00048 m_refIndexEc = tofPara->GetRefIndexEc(); 00049 m_phNConstEc = tofPara->GetPhNConstEc(); 00050 m_Cpe2pmtEc = tofPara->GetCpe2pmtEc(); 00051 m_rAngleEc = tofPara->GetRAngleEc(); 00052 m_QEEc = tofPara->GetQEEc(); 00053 m_CEEc = tofPara->GetCEEc(); 00054 m_peCorFacEc = tofPara->GetPeCorFacEc(); 00055 m_attenEc = tofPara->GetAttenEc(); 00056 00057 m_ttsMeanEc = tofPara->GetTTSmeanEc(); 00058 m_ttsSigmaEc = tofPara->GetTTSsigmaEc(); 00059 m_PMTgainEc = tofPara->GetPMTgainEc(); 00060 m_CeEc = tofPara->GetCeEc(); 00061 m_riseTimeEc = tofPara->GetRiseTimeEc(); 00062 m_LLthreshEc = tofPara->GetLLthreshEc(); 00063 m_HLthreshEc = tofPara->GetHLthreshEc(); 00064 m_preGainEc = tofPara->GetPreGainEc(); 00065 m_noiseSigmaEc = tofPara->GetNoiseSigmaEc(); 00066 }
|
|
|
|
00247 { 00248 G4double tmp_tauRatio,tmp_tau1,tmp_tau2,tmp_tau3; 00249 tmp_tauRatio = m_tauRatioEc; 00250 tmp_tau1 = m_tau1Ec; 00251 tmp_tau2 = m_tau2Ec; 00252 tmp_tau3 = m_tau3Ec; 00253 00254 G4double UniformR = tmp_tauRatio/(1+tmp_tauRatio); 00255 G4double EmissionTime; 00256 if (G4UniformRand()>UniformR) { 00257 while (1) { 00258 EmissionTime = -tmp_tau2*log( G4UniformRand() ); 00259 if (G4UniformRand()-exp(EmissionTime/tmp_tau2-EmissionTime/tmp_tau1)>1.E-8) 00260 break; 00261 } 00262 } 00263 else EmissionTime = -tmp_tau3*log( G4UniformRand() ); 00264 return EmissionTime; 00265 }
|
|
|
|
00148 { 00149 G4double cvelScint=c_light/m_refIndexEc; 00150 //Get information of this step 00151 G4ThreeVector pos=hit->GetPos(); 00152 G4int trackIndex = hit->GetTrackIndex(); 00153 G4int partId =hit->GetPartId(); 00154 G4double edep=hit->GetEdep(); 00155 G4double stepL=hit->GetStepL(); 00156 G4double deltaT=hit->GetDeltaT(); 00157 G4double timeFlight=hit->GetTime()-m_beamTime; 00158 if (timeFlight < m_globalTime) 00159 { 00160 m_globalTime = timeFlight; 00161 m_trackIndex = trackIndex; 00162 } 00163 00164 G4ThreeVector pDirection=hit->GetPDirection(); 00165 G4double nx=pDirection.x(); 00166 G4double ny=pDirection.y(); 00167 G4double nz=pDirection.z(); 00168 00169 //Cpe2pmt=cathode area/scint area 00170 //peCorFac=correction factor for eff. of available Npe 00171 00172 //number of photons generated in this step 00173 G4int NphStep; 00174 NphStep=G4int(edep*m_phNConstEc*m_Cpe2pmtEc*0.66*m_QEEc*m_CEEc*m_peCorFacEc); 00175 00176 G4double ddS, ddT; 00177 if (NphStep>0) 00178 { 00179 for (G4int i=0;i<NphStep;i++) 00180 { 00181 //uniform scintilation in each step 00182 ddS=stepL*G4UniformRand(); 00183 ddT=deltaT*G4UniformRand(); 00184 G4ThreeVector emtPos; 00185 emtPos.setX(pos.x() + nx*ddS); 00186 emtPos.setY(pos.y() + ny*ddS); 00187 emtPos.setZ(pos.z() + nz*ddS); 00188 00189 //check scinillation light whether to hit the pmt or not 00190 //forb=0/1 for forward(z>0, east) and backward(z<0, west) 00191 G4double pathL=0; 00192 G4int forb; 00193 DirectPh(partId, emtPos, pathL, forb); 00194 00195 //check if photon can reach PMT or not, after attenuation 00196 G4double ran=G4UniformRand(); 00197 if (pathL>0 && exp(-pathL/m_attenEc) > ran) 00198 { 00199 //propagation time in scintillator 00200 G4double scinSwim=pathL/cvelScint; 00201 //scintillation timing 00202 //G4double scinTime=GenPhoton(partId); 00203 G4double scinTime=Scintillation(partId); 00204 00205 //PMT transit time 00206 G4double transitTime=TransitTime(); 00207 //sum up all time components 00208 G4double endTime= timeFlight + ddT + scinSwim + scinTime + transitTime; 00209 00210 //store timing into binned buffer 00211 AccuSignal(endTime, forb); 00212 00213 //update 1st and last timings here 00214 if (m_t1st[forb]>endTime) m_t1st[forb] = endTime; 00215 if (m_tLast[forb]<endTime) m_tLast[forb]= endTime; 00216 } 00217 } 00218 } 00219 }
|
|
|
|
00132 { 00133 m_trackIndex=-999; 00134 m_globalTime = 9999; 00135 m_t1st[0]=100; 00136 m_t1st[1]=100; 00137 m_tLast[0]=0.; 00138 m_tLast[1]=0; 00139 m_totalPhot[0]=0; 00140 m_totalPhot[1]=0; 00141 for (G4int i=0;i<2;i++) 00142 for (G4int j=0;j<m_profBinNEcV2;j++) 00143 m_nPhot[j][i]=0; 00144 00145 }
|
|
|
|
00312 { 00313 //to generate PMT signal shape for single photoelectron. 00314 //only one time for a job. 00315 static G4double snpe[m_snpeBinNEcV2]; 00316 static G4int istore_snpe=-1; 00317 00318 //Model: f(t)=Gain*mv_1pe* t**2 * exp-(t/tau)**2/normal-const 00319 //normalization const =sqrt(pi)*tau*tau*tau/4.0 00320 G4double tau = m_riseTimeEc; 00321 G4double norma_const=sqrt(M_PI)*tau*tau*tau/4.0; //in unit of ns**3 00322 G4double echarge=1.6e-7; //in unit of PC 00323 00324 //time profile of pmt signals for Back and Forward PMT. 00325 G4double profPmt[m_profBinNEcV2][2]; 00326 00327 G4double t; 00328 G4int n1, n2, ii; 00329 G4int phtn; 00330 00331 if (istore_snpe<0) 00332 { 00333 istore_snpe = 1; 00334 for (G4int i=0;i<m_snpeBinNEcV2;i++) 00335 { 00336 t=(i+1)*m_timeBinSize; 00337 snpe[i]=m_PMTgainEc*m_CeEc*t*t*exp(- (t/tau) * (t/tau) )/norma_const; 00338 } 00339 } 00340 //for barrel and endcap tof: fb=2 or 1 00341 G4int fb=1; 00342 00343 for (G4int j=0;j<2;j++) 00344 { 00345 m_TDC[j]=-999.0; 00346 m_ADC[j]=-999.0; 00347 } 00348 for (G4int j=0; j<fb; j++) 00349 { 00350 if (m_totalPhot[j] > 0) 00351 { 00352 n1=G4int(m_t1st[j]/m_timeBinSize); 00353 n2=G4int(m_tLast[j]/m_timeBinSize); 00354 00355 for (G4int i=0;i<m_profBinNEcV2;i++) 00356 profPmt[i][j]=0.0; 00357 00358 //generate PMT pulse 00359 n2 = n2<m_profBinNEcV2 ? n2:m_profBinNEcV2; 00360 for (G4int i=n1;i<n2;i++) 00361 { 00362 phtn=m_nPhot[i][j]; 00363 if (phtn>0) 00364 { 00365 for (G4int ihst=0; ihst<m_snpeBinNEcV2; ihst++) 00366 { 00367 ii=i+ihst; 00368 if (ii<m_profBinNEcV2) 00369 profPmt[ii][j] += phtn*snpe[ihst]; 00370 } 00371 } 00372 } 00373 00374 //add preamplifier and noise 00375 for (int i=0;i<m_profBinNEcV2;i++) 00376 { 00377 if (profPmt[i][j]>0) 00378 profPmt[i][j] = m_preGainEc*profPmt[i][j]+G4RandGauss::shoot(0,m_noiseSigmaEc); 00379 } 00380 00381 //get pulse height 00382 G4double max=0; 00383 for (int i=0;i<m_profBinNEcV2;i++) 00384 { 00385 if (profPmt[i][j]>max) 00386 max=profPmt[i][j]; 00387 } 00388 00389 G4double tmp_HLthresh, tmp_LLthresh; 00390 if (partId==1) { 00391 tmp_HLthresh = m_HLthreshEc; 00392 tmp_LLthresh = m_LLthreshEc; 00393 } 00394 else { 00395 tmp_HLthresh = m_HLthreshEc; 00396 tmp_LLthresh = m_LLthreshEc; 00397 } 00398 00399 //get final tdc and adc 00400 if (max>=tmp_HLthresh) 00401 { 00402 for (int i=0;i<m_profBinNEcV2;i++) 00403 { 00404 if ( profPmt[i][j] >= tmp_LLthresh ) 00405 { 00406 m_TDC[j] = i*m_timeBinSize; 00407 m_ADC[j] = m_preGainEc*m_totalPhot[j]*echarge*m_PMTgainEc; 00408 break; 00409 } 00410 } 00411 } 00412 } 00413 } 00414 }
|
|
|
|
00295 { 00296 //get time transit spread 00297 return G4RandGauss::shoot(m_ttsMeanEc,m_ttsSigmaEc); 00298 }
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Reimplemented from BesTofDigitizerV. |
|
Reimplemented from BesTofDigitizerV. |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|