00001
00011
00012 #include "BesTofDigitizerEcV2.hh"
00013 #include "BesTofHit.hh"
00014 #include "G4DigiManager.hh"
00015 #include "G4RunManager.hh"
00016 #include "Randomize.hh"
00017 #include "G4Poisson.hh"
00018 #include "BesTofDigi.hh"
00019 #include "BesTofGeoParameter.hh"
00020 #include "GaudiKernel/ISvcLocator.h"
00021 #include "GaudiKernel/Bootstrap.h"
00022 #include "GaudiKernel/IDataProviderSvc.h"
00023 #include "G4Svc/IG4Svc.h"
00024 #include "G4Svc/G4Svc.h"
00025
00026 BesTofDigitizerEcV2::BesTofDigitizerEcV2()
00027 {
00028 ReadData();
00029 m_timeBinSize=0.005;
00030
00031
00032 ISvcLocator* svcLocator = Gaudi::svcLocator();
00033 IG4Svc* tmpSvc;
00034 StatusCode sc = svcLocator->service("G4Svc", tmpSvc);
00035 m_G4Svc = dynamic_cast<G4Svc *>(tmpSvc);
00036
00037 }
00038
00039 void BesTofDigitizerEcV2::ReadData()
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 }
00067
00068 BesTofDigitizerEcV2::~BesTofDigitizerEcV2()
00069 {
00070 ;
00071 }
00072
00073 void BesTofDigitizerEcV2::Digitize(ScintSingle* scint, BesTofDigitsCollection* DC)
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
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
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 }
00130
00131 void BesTofDigitizerEcV2::TofPmtInit()
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 }
00146
00147 void BesTofDigitizerEcV2::TofPmtAccum(BesTofHit* hit)
00148 {
00149 G4double cvelScint=c_light/m_refIndexEc;
00150
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
00170
00171
00172
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
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
00190
00191 G4double pathL=0;
00192 G4int forb;
00193 DirectPh(partId, emtPos, pathL, forb);
00194
00195
00196 G4double ran=G4UniformRand();
00197 if (pathL>0 && exp(-pathL/m_attenEc) > ran)
00198 {
00199
00200 G4double scinSwim=pathL/cvelScint;
00201
00202
00203 G4double scinTime=Scintillation(partId);
00204
00205
00206 G4double transitTime=TransitTime();
00207
00208 G4double endTime= timeFlight + ddT + scinSwim + scinTime + transitTime;
00209
00210
00211 AccuSignal(endTime, forb);
00212
00213
00214 if (m_t1st[forb]>endTime) m_t1st[forb] = endTime;
00215 if (m_tLast[forb]<endTime) m_tLast[forb]= endTime;
00216 }
00217 }
00218 }
00219 }
00220
00221
00222 void BesTofDigitizerEcV2::DirectPh(G4int partId, G4ThreeVector emtPos, G4double& pathL, G4int& forb)
00223 {
00224
00225
00226
00227 G4double cos_spanEc = 1;
00228 G4double ran;
00229
00230 ran=G4UniformRand();
00231
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;
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 }
00245
00246 G4double BesTofDigitizerEcV2::Scintillation(G4int partId)
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 }
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294 G4double BesTofDigitizerEcV2::TransitTime()
00295 {
00296
00297 return G4RandGauss::shoot(m_ttsMeanEc,m_ttsSigmaEc);
00298 }
00299
00300 void BesTofDigitizerEcV2::AccuSignal(G4double endTime, G4int forb)
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 }
00310
00311 void BesTofDigitizerEcV2::TofPmtRspns(G4int partId)
00312 {
00313
00314
00315 static G4double snpe[m_snpeBinNEcV2];
00316 static G4int istore_snpe=-1;
00317
00318
00319
00320 G4double tau = m_riseTimeEc;
00321 G4double norma_const=sqrt(M_PI)*tau*tau*tau/4.0;
00322 G4double echarge=1.6e-7;
00323
00324
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
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
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
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
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
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 }
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425