BesTofDigitizerEcV3 Class Reference

#include <BesTofDigitizerEcV3.hh>

Inheritance diagram for BesTofDigitizerEcV3:


Public Member Functions
	BesTofDigitizerEcV3 ()
	~BesTofDigitizerEcV3 ()
virtual void	Digitize (ScintSingle , BesTofDigitsCollection )
void	ReadData ()
void	ReadEffTree ()
void	TofPmtInit ()
void	TofPmtAccum (BesTofHit *)
void	DirectPh (G4int, G4int, G4int, G4double &)
G4double	Scintillation (G4int)
G4double	TransitTime ()
void	AccuSignal (G4double, G4int)
void	TofPmtRspns (G4int, G4int)
G4double	BirksLaw (BesTofHit *hit)
void	Initialize ()
Protected Attributes
BesTofDigitsCollection *	m_besTofDigitsCollection
BesTofHitsCollection *	m_THC
ITofCaliSvc *	m_tofCaliSvc
ITofSimSvc *	m_tofSimSvc
ITofQElecSvc *	m_tofQElecSvc
G4double	m_ADC [2]
G4double	m_TDC [2]
G4int	m_trackIndex
G4double	m_globalTime
Static Protected Attributes
static bool	m_booked = false
static NTuple::Tuple *	m_tupleTof1 = 0
static NTuple::Item< double >	m_partId
static NTuple::Item< double >	m_scinNb
static NTuple::Item< double >	m_edep
static NTuple::Item< double >	m_nHits
static NTuple::Item< double >	m_time1st0
static NTuple::Item< double >	m_time1st1
static NTuple::Item< double >	m_timelast0
static NTuple::Item< double >	m_timelast1
static NTuple::Item< double >	m_totalPhot0
static NTuple::Item< double >	m_totalPhot1
static NTuple::Item< double >	m_NphAllSteps
static NTuple::Item< double >	m_max0
static NTuple::Item< double >	m_max1
static NTuple::Item< double >	m_tdc0
static NTuple::Item< double >	m_adc0
static NTuple::Item< double >	m_tdc1
static NTuple::Item< double >	m_adc1
static NTuple::Tuple *	m_tupleTof2 = 0
static NTuple::Item< double >	m_eTotal
static NTuple::Item< double >	m_nDigi
static NTuple::Item< double >	m_partIdMPV
static NTuple::Item< double >	m_scinNbMPV
static NTuple::Item< double >	m_edepMPV
static NTuple::Item< double >	m_nDigiOut
static NTuple::Tuple *	m_tupleTof3 = 0
static NTuple::Item< double >	m_forb
static NTuple::Item< double >	m_timeFlight
static NTuple::Item< double >	m_ddT
static NTuple::Item< double >	m_scinSwim
static NTuple::Item< double >	m_scinTime
static NTuple::Item< double >	m_transitTime
static NTuple::Item< double >	m_endTime
static NTuple::Item< double >	m_edepHit
Private Attributes
G4double	m_ecR1
G4double	m_tau1Ec
G4double	m_tau2Ec
G4double	m_tau3Ec
G4double	m_tauRatioEc
G4double	m_refIndexEc
G4double	m_phNConstEc
G4double	m_Cpe2pmtEc
G4double	m_rAngleEc
G4double	m_QEEc
G4double	m_CEEc
G4double	m_peCorFacEc
G4double	m_attenEc
G4double	m_ttsMeanEc
G4double	m_ttsSigmaEc
G4double	m_PMTgainEc
G4double	m_CeEc
G4double	m_riseTimeEc
G4double	m_LLthreshEc
G4double	m_HLthreshEc
G4double	m_preGainEc
G4double	m_noiseSigmaEc
G4double	m_t1st [2]
G4double	m_tLast [2]
G4double	m_timeBinSize
G4int	m_totalPhot [2]
G4int	m_nPhot [m_profBinNEcV3][2]
G4Svc *	m_G4Svc
G4double	m_beamTime
RealizationSvc *	m_RealizationSvc
float	propTime [50][10][10][num1]
float	prob [50][10][10][num1]
float	eff [50][10][10]

Detailed Description

Definition at line 40 of file BesTofDigitizerEcV3.hh.

Constructor & Destructor Documentation

BesTofDigitizerEcV3::BesTofDigitizerEcV3 ( )

Definition at line 31 of file BesTofDigitizerEcV3.cc.

References eff, genRecEmupikp::i, ganga-rec::j, m_G4Svc, m_RealizationSvc, m_timeBinSize, num1, prob, propTime, ReadData(), and ReadEffTree().

00032 {
00033   ReadData(); //Get some basic data (not Hit data, but information about smearing from electronics,... )
00034   m_timeBinSize = 0.005;
00035 
00036 
00037 
00038 
00039   // Simulation/G4Svc/G4Svc-00-01-47/src/G4Svc.cpp
00040   //Get the basic parameters from the event: Starttime, Beamposition,...
00041 
00042   //retrieve G4Svc   
00043   ISvcLocator* svcLocator = Gaudi::svcLocator();
00044   IG4Svc* tmpSvc;
00045   StatusCode sc = svcLocator->service("G4Svc", tmpSvc);
00046   if(!sc.isSuccess())
00047   {
00048     std::cout << " Could not initialize Realization Service in BesTofDigitizerBrV2" << std::endl;
00049   }
00050   else
00051   {
00052     m_G4Svc = dynamic_cast<G4Svc *>(tmpSvc);
00053   }
00054 
00055   //retrieve RealizationSvc
00056   IRealizationSvc *tmpReal;
00057   StatusCode scReal = svcLocator->service("RealizationSvc",tmpReal);
00058   if (!scReal.isSuccess())
00059   {
00060     std::cout << " Could not initialize Realization Service in BesTofDigitizerEcV3" << std::endl;
00061   } 
00062   else 
00063   {
00064     m_RealizationSvc = dynamic_cast<RealizationSvc*>(tmpReal);
00065   }
00066 
00067 
00068 
00069   for(int i=0;i<50;i++)
00070   {
00071     for(int j=0;j<10;j++)
00072     {
00073       for(int k=0;k<10;k++)
00074       {
00075         for(int m=0;m<num1;m++) // num1 bin number in histogram is 400 at the moment
00076         {  
00077           //G4cout << "time:" << propTime[i][j][k][m] << "; prob:" << prob[i][j][k][m] << "; eff:" << eff[i][j][k] << G4endl;
00078           propTime[i][j][k][m] = 0;
00079           prob[i][j][k][m] = 0;
00080           eff[i][j][k] = 0;
00081         }
00082       }
00083     }
00084   }
00085 
00086   ReadEffTree();
00087   G4cout << "ETofSim: Reading nTuples of is completed." << G4endl;
00088 }

BesTofDigitizerEcV3::~BesTofDigitizerEcV3 ( )

Definition at line 177 of file BesTofDigitizerEcV3.cc.

00178 {;}

Member Function Documentation

void BesTofDigitizerEcV3::AccuSignal	(	G4double	,
		G4int
	)

Definition at line 556 of file BesTofDigitizerEcV3.cc.

References m_nPhot, m_profBinNEcV3, m_timeBinSize, and m_totalPhot.

Referenced by TofPmtAccum().

00557 {
00558   G4int ihst;
00559   ihst=G4int(endTime/m_timeBinSize);
00560   if (ihst>0 &&ihst<m_profBinNEcV3)
00561   {
00562     m_nPhot[ihst][forb]=m_nPhot[ihst][forb]+1;
00563     m_totalPhot[forb]=m_totalPhot[forb]+1;
00564   }
00565 }

G4double BesTofDigitizerEcV3::BirksLaw ( BesTofHit * hit )

Definition at line 476 of file BesTofDigitizerEcV3.cc.

References BesTofHit::GetCharge(), BesTofHit::GetEdep(), and BesTofHit::GetStepL().

Referenced by TofPmtAccum().

00477 {
00478   const G4double kappa = 0.015*cm/MeV;
00479   const G4String brMaterial = "BC404";
00480   G4double dE = hit->GetEdep();
00481   //G4cout << "The edep is "<< dE << G4endl;
00482   G4double dX = hit->GetStepL();
00483   //G4Material* materiral = hit->GetMaterial();
00484   G4double charge = hit->GetCharge();
00485   G4double cor_dE = dE;
00486   //if((materiral->GetName()==brMaterial) && charge!=0.&& dX!=0.)
00487   if(charge!=0.&& dX!=0.)
00488   {
00489     cor_dE = dE/(1+kappa*dE/dX);
00490     //if(dE>20)
00491     //{
00492     //  G4cout << "\n dE > 20. Details are below:" << G4endl;
00493     //  G4cout << "dE/dx:" << dE/dX << G4endl;
00494     //  G4cout << "dE:" << dE << "; dX:" << dX << G4endl;
00495     //  G4cout << "It is BC404. cor_dE is " << cor_dE << G4endl;
00496     //  G4double ratio = cor_dE/dE;
00497     //  G4cout << "The ratio cor_dE/edep is "<< ratio << G4endl;
00498     //}
00499     //G4cout << "It is BC408. cor_dE is " << cor_dE << G4endl;
00500     //G4double ratio = cor_dE/dE;
00501     //G4cout << "The ratio cor_dE/edep is "<< ratio << G4endl;
00502   }
00503   return cor_dE;
00504 
00505 }

void BesTofDigitizerEcV3::Digitize	(	ScintSingle *	,
		BesTofDigitsCollection *
	)			`[virtual]`

Reimplemented from BesTofDigitizerV.

Definition at line 180 of file BesTofDigitizerEcV3.cc.

Referenced by BesTofDigitizer::Digitize().

00181 {
00182   m_beamTime = m_G4Svc->GetBeamTime() * ns;
00183   m_besTofDigitsCollection = DC;
00184 
00185   G4DigiManager* digiManager = G4DigiManager::GetDMpointer();
00186 
00187   G4int THCID = digiManager->GetHitsCollectionID("BesTofHitsCollection");
00188   m_THC = (BesTofHitsCollection*) (digiManager->GetHitsCollection(THCID));
00189 
00190   if (m_G4Svc->TofRootFlag())
00191   {
00192     m_eTotal = 0;
00193     m_nDigi = 0;
00194     m_partIdMPV = -9;
00195     m_scinNbMPV = -9;
00196     m_edepMPV = 0;
00197     m_nDigiOut = 0;
00198   }
00199 
00200   if (m_THC)
00201   {
00202     //for each digi, compute TDC and ADC
00203     G4int partId, scinNb, nHits;
00204     G4double edep;
00205     BesTofHit* hit;
00206     partId=scint->GetPartId();
00207     scinNb=scint->GetScinNb();
00208     edep = scint->GetEdep();
00209     nHits=scint->GetHitIndexes()->size();
00210 
00211     TofPmtInit();
00212 
00213     //fill tof Ntuple
00214     if (m_G4Svc->TofRootFlag())
00215     {
00216       if (edep>m_edepMPV)
00217       {
00218         m_partIdMPV = partId;
00219         m_scinNbMPV = scinNb;
00220         m_edepMPV = edep;
00221       }
00222       m_eTotal += edep;
00223       m_nDigi ++;
00224 
00225       m_partId = partId;
00226       m_scinNb = scinNb;
00227       m_edep = edep;
00228       m_nHits = nHits;
00229     }
00230 
00231     if (edep>0.01)
00232     {
00233       for (G4int j=0;j<nHits;j++)
00234       {
00235         hit= (*m_THC)[( *(scint->GetHitIndexes()) )[j]];
00236         TofPmtAccum(hit);
00237       }
00238 
00239       if (m_G4Svc->TofRootFlag())
00240       {
00241         m_time1st0=m_t1st[0];
00242         m_time1st1=m_t1st[1];
00243         m_timelast0=m_tLast[0];
00244         m_timelast1=m_tLast[1];
00245         m_totalPhot0=m_totalPhot[0];
00246         m_totalPhot1=m_totalPhot[1];
00247       }
00248 
00249       //get final tdc and adc
00250       TofPmtRspns(partId,scinNb);
00251 
00252       G4double temp0 = m_ADC[0]+m_TDC[0];
00253       G4double temp1 = m_ADC[1]+m_TDC[1];
00254       //const double MAX_ADC = 8191*0.3; // channel set up to 8192 will lead to overflow.
00255       if ( (partId!=1) && temp0>0. )
00256       {
00257         BesTofDigi* digi = new BesTofDigi;
00258         digi->SetTrackIndex(m_trackIndex);
00259         digi->SetPartId(partId);
00260         digi->SetScinNb(scinNb);
00261         digi->SetForwADC( m_ADC[0]) ;
00262         digi->SetBackADC( m_ADC[1]) ;
00263         if (m_TDC[0]>0.)
00264           m_TDC[0] = m_TDC[0]+m_beamTime;
00265         digi->SetForwTDC( m_TDC[0]) ;
00266         digi->SetBackTDC( m_TDC[1]) ;
00267   
00268         //G4cout << "BesTofDigitizerEcV3:     ForwTDC   |   BackTDC  "  << m_TDC[0] << "  |  " << m_TDC[1] << G4endl;
00269         //G4cout<<"endcap\nadc0:"<<m_ADC[0]<<"; adc1:"<<m_ADC[1]<<"; tdc0:"<<m_TDC[0]<<"; tdc1:"<<m_TDC[1]<<G4endl;
00270         m_besTofDigitsCollection->insert(digi);
00271         if (m_G4Svc->TofRootFlag() )
00272           m_nDigiOut++;
00273       }
00274       if (m_G4Svc->TofRootFlag()  )
00275         m_tupleTof1->write();
00276       //cout << "m_tupleTof1->write()" << endl;
00277     }
00278 
00279   }
00280   if (m_G4Svc->TofRootFlag())
00281     m_tupleTof2->write();
00282   //cout << "m_tupleTof2->write()" << endl;
00283 
00284 }

void BesTofDigitizerEcV3::DirectPh	(	G4int	,
		G4int	,
		G4int	,
		G4double &
	)

Definition at line 507 of file BesTofDigitizerEcV3.cc.

References key, prob, and propTime.

Referenced by TofPmtAccum().

00508 {
00509   G4double ran = G4UniformRand();
00510   G4double p = 0; 
00511   G4int nth = 1;
00512   G4int key = 0;
00513   t = 0;
00514   while(1)
00515   {
00516     if(p>ran||nth==400)
00517     {
00518       key = nth;
00519       //G4cout << "Value found!" << G4endl;
00520       break;
00521     }
00522     p = p + prob[rBin][phiBin][zBin][nth];
00523     nth++;
00524   }
00525   t = propTime[rBin][phiBin][zBin][key-1];
00526 }

void BesTofDigitizerV::Initialize ( ) [inherited]

Definition at line 189 of file BesTofDigitizerV.cc.

References genRecEmupikp::i, BesTofDigitizerV::m_ADC, BesTofDigitizerV::m_globalTime, BesTofDigitizerV::m_TDC, and BesTofDigitizerV::m_trackIndex.

Referenced by BesTofDigitizerEcV1::TofPmtInit(), BesTofDigitizerBrV2::TofPmtInit(), and BesTofDigitizerBrV1::TofPmtInit().

00190 {
00191     for (G4int i=0;i<2;i++)
00192     {
00193         m_ADC[i] = -999;
00194         m_TDC[i] = -999;
00195     }
00196     m_trackIndex = -999;
00197     m_globalTime = 9999;
00198 }

void BesTofDigitizerEcV3::ReadData ( )

Definition at line 90 of file BesTofDigitizerEcV3.cc.

Referenced by BesTofDigitizerEcV3().

00091 {
00092   BesTofGeoParameter* tofPara = BesTofGeoParameter::GetInstance();
00093 
00094   m_ecR1         = tofPara->GetEcR1();
00095   m_tau1Ec       = tofPara->GetTau1Ec();
00096   m_tau2Ec       = tofPara->GetTau2Ec();
00097   m_tau3Ec       = tofPara->GetTau3Ec();
00098   m_tauRatioEc   = tofPara->GetTauRatioEc();
00099   m_refIndexEc   = tofPara->GetRefIndexEc();
00100   m_phNConstEc   = tofPara->GetPhNConstEc();
00101   m_Cpe2pmtEc    = tofPara->GetCpe2pmtEc();
00102   m_rAngleEc     = tofPara->GetRAngleEc();
00103   m_QEEc         = tofPara->GetQEEc();
00104   m_CEEc         = tofPara->GetCEEc();
00105   m_peCorFacEc   = tofPara->GetPeCorFacEc();
00106   m_attenEc      = tofPara->GetAttenEc();
00107 
00108   m_ttsMeanEc    = tofPara->GetTTSmeanEc();
00109   m_ttsSigmaEc   = tofPara->GetTTSsigmaEc();
00110   m_PMTgainEc    = tofPara->GetPMTgainEc();
00111   m_CeEc         = tofPara->GetCeEc();
00112   m_riseTimeEc   = tofPara->GetRiseTimeEc();
00113   m_LLthreshEc   = tofPara->GetLLthreshEc();
00114   m_HLthreshEc   = tofPara->GetHLthreshEc();
00115   m_preGainEc    = tofPara->GetPreGainEc();
00116   m_noiseSigmaEc = tofPara->GetNoiseSigmaEc();
00117 
00118   //G4cout << "m_LLthreshEc:" << m_LLthreshEc << "; m_HLthreshEc:" << m_HLthreshEc << G4endl;
00119 
00120 }

void BesTofDigitizerEcV3::ReadEffTree ( )

Definition at line 123 of file BesTofDigitizerEcV3.cc.

References eff, RealizationSvc::getRunId(), genRecEmupikp::i, ganga-rec::j, m_RealizationSvc, nPhi, nR, nZ, prob, propTime, t(), and x.

Referenced by BesTofDigitizerEcV3().

00124 {
00125   int rBin,phiBin,zBin;
00126   const int nR = 43;
00127   const int nPhi = 6;
00128   const int nZ = 6;
00129   float efficiency0,x[400],y[400]; 
00130 
00131 
00132   G4String dataPath = getenv("TOFSIMROOT");
00133   if(!dataPath)
00134   {
00135     G4Exception("Boss environment is not set!");
00136   }
00137 
00138   char treePath[200];
00139   G4int runId = m_RealizationSvc->getRunId();
00140   if(runId>=-80000 && runId<=-9484)
00141   {
00142     // After TOF HV adjustment, endcap attenL was set to 5000mm.
00143     sprintf(treePath,"%s/dat/effTree_1600mm.root",dataPath.c_str());
00144   }
00145   else
00146   {
00147     //Before TOF HV adjustment, endcap attenL was set to 1600mm.
00148     sprintf(treePath,"%s/dat/effTree_1600mm.root",dataPath.c_str());
00149   }
00150 
00151   TFile *f = new TFile(treePath, "read");
00152   TTree *t = (TTree*)f->Get("effTree");
00153 
00154   t->SetBranchAddress("rBin", &rBin);
00155   t->SetBranchAddress("phiBin", &phiBin);
00156   t->SetBranchAddress("zBin", &zBin);
00157   t->SetBranchAddress("efficiency0", &efficiency0);
00158   t->SetBranchAddress("x", x);
00159   t->SetBranchAddress("y", y);
00160 
00161   int r,phi,z;
00162   for (Int_t i = 0; i < nR*nPhi*nZ; i++){
00163     t->GetEntry(i);
00164     r = rBin;
00165     phi = phiBin;
00166     z = zBin;
00167     eff[r][phi][z] = efficiency0;
00168     for (Int_t j = 0; j < 400; j++){
00169       propTime[r][phi][z][j] = x[j];
00170       prob[r][phi][z][j] = y[j];
00171       //cout << "\n" << propTime[r][phi][z][j] << "   " << prob[r][phi][z][j] << endl;
00172     }
00173   }
00174 
00175 }

G4double BesTofDigitizerEcV3::Scintillation ( G4int )

Definition at line 528 of file BesTofDigitizerEcV3.cc.

References exp(), m_tau1Ec, m_tau2Ec, m_tau3Ec, and m_tauRatioEc.

Referenced by TofPmtAccum().

00529 {
00530   G4double tmp_tauRatio,tmp_tau1,tmp_tau2,tmp_tau3;
00531   tmp_tauRatio = m_tauRatioEc;
00532   tmp_tau1 = m_tau1Ec;
00533   tmp_tau2 = m_tau2Ec;
00534   tmp_tau3 = m_tau3Ec;
00535 
00536   G4double UniformR = tmp_tauRatio/(1+tmp_tauRatio);
00537   G4double EmissionTime;
00538   if (G4UniformRand()>UniformR) {
00539     while (1) {
00540       EmissionTime = -tmp_tau2*log( G4UniformRand() );
00541       if (G4UniformRand()-exp(EmissionTime/tmp_tau2-EmissionTime/tmp_tau1)>1.E-8)
00542         break;
00543     }
00544   }
00545   else EmissionTime = -tmp_tau3*log( G4UniformRand() );
00546   return EmissionTime;
00547 }

void BesTofDigitizerEcV3::TofPmtAccum ( BesTofHit * )

Definition at line 328 of file BesTofDigitizerEcV3.cc.

Referenced by Digitize().

00329 {
00330   G4double cvelScint = c_light/m_refIndexEc;
00331   //Get information of this step
00332   G4ThreeVector pos = hit->GetPos();
00333   G4int trackIndex = hit->GetTrackIndex();
00334   G4int partId = hit->GetPartId();
00335   G4double edep = hit->GetEdep();
00336   G4double stepL = hit->GetStepL();
00337   //G4String particleName = hit->GetPName();
00338   G4double deltaT=hit->GetDeltaT();
00339   G4double timeFlight=hit->GetTime()-m_beamTime;
00340   if (timeFlight < m_globalTime)
00341   {
00342     m_globalTime = timeFlight;
00343     m_trackIndex = trackIndex;
00344   }
00345 
00346   G4ThreeVector pDirection=hit->GetPDirection();
00347   G4double nx=pDirection.x();
00348   G4double ny=pDirection.y();
00349   G4double nz=pDirection.z();
00350 
00351 
00352   //number of photons generated in this step
00353   G4int NphStep;
00354   G4double nMean, nPhoton;
00355   nMean = m_phNConstEc*BirksLaw(hit);
00356 
00357   if(nMean>10)
00358   {
00359     G4double resolutionScale=1.;
00360     G4double sigma=resolutionScale*sqrt(nMean);
00361     nPhoton=G4int(G4RandGauss::shoot(nMean,sigma)+0.5);
00362   } 
00363   else
00364     nPhoton=G4int(G4Poisson(nMean));
00365 
00366 
00367   NphStep=G4int(nPhoton*0.66*m_QEEc*m_CEEc);
00368   //G4cout << "\nradius:" << radius << "; phi:" << phi << "; z:" << z << G4endl;
00369   //G4cout << "\nrBin:" << rBin << ";phiBin:" << phiBin << ";zBin" << zBin << G4endl;
00370 
00371   if (m_G4Svc->TofRootFlag())
00372     m_NphAllSteps += G4int(nPhoton*0.66*m_QEEc*m_CEEc);
00373 
00374   if (NphStep>0)
00375   {
00376     for (G4int i=0;i<NphStep;i++)
00377     {
00378       //uniform scintilation in each step
00379       G4double ddS, ddT;
00380       ddS=stepL*G4UniformRand();
00381       ddT=deltaT*G4UniformRand();
00382       G4ThreeVector emtPos;
00383       emtPos.setX(pos.x() + nx*ddS);
00384       emtPos.setY(pos.y() + ny*ddS);
00385       emtPos.setZ(pos.z() + nz*ddS);
00386 
00387       //retrieve the histogram info
00388       G4double radius = sqrt(emtPos.x()*emtPos.x()+emtPos.y()*emtPos.y())-m_ecR1;
00389       const G4double pie = 2.*asin(1.);
00390       G4double phi;
00391       if(emtPos.x()>0 && emtPos.y()>0) 
00392         phi = atan(emtPos.y()/emtPos.x());
00393       else if(emtPos.x()==0 && emtPos.y()>0)
00394         phi = pie/2.;
00395       else if(emtPos.x()<0)
00396         phi = atan(emtPos.y()/emtPos.x())+pie;
00397       else if(emtPos.x()==0 && emtPos.y()<0)
00398         phi = 1.5*pie;
00399       else if(emtPos.x()>0 && emtPos.y()<0)
00400         phi = 2.*pie+atan(emtPos.y()/emtPos.x());
00401       phi = phi*180./pie; // in degrees
00402       G4double z = fabs(emtPos.z());
00403       // Warning: Should obtain absolute value of z to determine zBinNum
00404 
00405       G4int rBin = G4int(radius/10.); // radius bin no.
00406       G4double resPhi = phi-(G4int(phi/7.5)*7.5); // residual of phi in period of 7.5
00407       G4int phiBin = G4int(resPhi/1.25);
00408       G4int zBin = G4int((z-1332.)/8.);
00409 
00410       //check scinillation light whether to hit the pmt or not
00411       //forb=0/1 for forward(z>0, east) and backward(z<0, west)
00412       G4int forb = 0;
00413       G4double transpTime = 0;
00414       G4double pathL = 0;
00415       G4double efficiency1;
00416       G4double efficiency2;
00417       efficiency1 = G4RandGauss::shoot(0,0.004);
00418       if(rBin>=0&&rBin<=nR && phiBin>=0&& phiBin<=nPhi && zBin>=0&&zBin<=nZ)
00419         efficiency1 += eff[rBin][phiBin][zBin];
00420       else
00421         efficiency1 = 0;
00422       //G4cout << "FATAL: The collection efficiency does NOT exist!" << G4endl;
00423       if(m_attenEc==0)
00424       {
00425         G4cout <<" ERROR: Attenuation Length is null!" << G4endl;
00426         break;
00427       }
00428       //efficiency2 = pow(efficiency1,(1600/m_attenEc));
00429       if(G4UniformRand() <= efficiency1)
00430       {
00431         DirectPh(rBin, phiBin, zBin, transpTime);
00432         //cout << "transpTime:" << transpTime << endl;
00433       }
00434 
00435       //check if photon can reach PMT or not, after attenuation
00436       //G4double ran = G4UniformRand();
00437       //pathL = transpTime*cvelScint;
00438       //if (pathL>0 && exp(-pathL/m_attenEc) > ran)  // Note: Do NOT double count attuation!
00439       if(transpTime>0)
00440       {
00441         //propagation time in scintillator
00442         G4double scinSwim = transpTime;
00443         //scintillation timing
00444         G4double scinTime = Scintillation(partId);
00445 
00446         //PMT transit time
00447         G4double transitTime = TransitTime();
00448         //sum up all time components
00449         G4double endTime = timeFlight + ddT + scinSwim + scinTime + transitTime;
00450 
00451         if (m_G4Svc->TofRootFlag())
00452         {
00453           //m_forb = forb;
00454           m_timeFlight = timeFlight+ddT;
00455           m_ddT = ddT;
00456           m_scinSwim = scinSwim;
00457           m_scinTime = scinTime;
00458           m_transitTime = transitTime;
00459           m_endTime = endTime;
00460           m_tupleTof3->write();
00461         }
00462 
00463         //store timing into binned buffer
00464         AccuSignal(endTime, forb);
00465 
00466         //update 1st and last timings here
00467         if (m_t1st[forb]>endTime)   m_t1st[forb] = endTime;
00468         if (m_tLast[forb]<endTime)  m_tLast[forb]= endTime;
00469         //if(m_tLast[0]>100)
00470         //std::cout<<"endTime: "<<endTime<<std::endl;
00471       }
00472     }
00473   }
00474 }

void BesTofDigitizerEcV3::TofPmtInit ( )

Definition at line 286 of file BesTofDigitizerEcV3.cc.

Referenced by Digitize().

00287 {
00288   m_ADC[0] = -999.;
00289   m_ADC[1] = -999.;
00290   m_TDC[0] = -999.;
00291   m_TDC[1] = -999.;
00292   m_trackIndex = -999;
00293   m_globalTime = 9999;
00294 
00295   m_t1st[0]=100;
00296   m_t1st[1]=100;
00297   m_tLast[0]=0.;
00298   m_tLast[1]=0;
00299   m_totalPhot[0]=0;
00300   m_totalPhot[1]=0;
00301   for (G4int i=0;i<2;i++)
00302     for (G4int j=0;j<m_profBinNEcV3;j++)
00303       m_nPhot[j][i]=0;
00304 
00305   if (m_G4Svc->TofRootFlag())
00306   {
00307     m_partId = -9;
00308     m_scinNb = -9;
00309     m_edep = 0;
00310     m_nHits = 0;
00311     m_time1st0 = 100;
00312     m_time1st1 = 100;
00313     m_timelast0 = 0;
00314     m_timelast1 = 0;
00315     m_totalPhot0 = 0;
00316     m_totalPhot1 = 0;
00317     m_NphAllSteps = 0;
00318     m_max0 = 0;
00319     m_max1 = 0;
00320     m_tdc0 = -999;
00321     m_adc0 = -999;
00322     m_tdc1 = -999;
00323     m_adc1 = -999;
00324   }
00325 
00326 }

void BesTofDigitizerEcV3::TofPmtRspns	(	G4int	,
		G4int
	)

Definition at line 567 of file BesTofDigitizerEcV3.cc.

Referenced by Digitize().

00568 {
00569   //to generate PMT signal shape for single photoelectron.
00570   //only one time for a job.
00571   static G4double snpe[m_snpeBinNEcV3];
00572   static G4int istore_snpe=-1;
00573 
00574   //Model: f(t)=Gain*mv_1pe* t**2 * exp-(t/tau)**2/normal-const
00575   //normalization const =sqrt(pi)*tau*tau*tau/4.0
00576   G4double tau = m_riseTimeEc;
00577   G4double norma_const=sqrt(M_PI)*tau*tau*tau/4.0;    //in unit of ns**3
00578   G4double echarge=1.6e-7;   //in unit of PC
00579 
00580   //time profile of pmt signals for Back and Forward PMT.
00581   G4double profPmt[m_profBinNEcV3][2];
00582 
00583   G4double t;
00584   G4int n1, n2, ii;
00585   G4int phtn;
00586 
00587   if (istore_snpe<0)
00588   {
00589     istore_snpe = 1;
00590     for (G4int i=0;i<m_snpeBinNEcV3;i++)
00591     {
00592       t=(i+1)*m_timeBinSize;
00593       snpe[i]=m_CeEc*t*t*exp(- (t/tau) * (t/tau) )/norma_const;
00594     }
00595   }
00596   //for barrel and endcap tof: fb=2 or 1
00597   G4int fb=1;   // for endcap part
00598 
00599   G4double tmpADC[2] = {0,0};
00600 
00601   for (G4int j=0; j<fb; j++)
00602   {
00603     if (m_totalPhot[j] > 0)
00604     {
00605       n1=G4int(m_t1st[j]/m_timeBinSize);
00606       n2=G4int(m_tLast[j]/m_timeBinSize);
00607 
00608       for (G4int i=0;i<m_profBinNEcV3;i++)
00609         profPmt[i][j]=0.0;
00610 
00611       //generate PMT pulse
00612       n2 = n2<m_profBinNEcV3 ? n2:m_profBinNEcV3;
00613       for (G4int i=n1;i<n2;i++)
00614       {
00615         phtn=m_nPhot[i][j];
00616         if (phtn>0)
00617         {
00618           G4double Npoisson;
00619           while(1) {
00620             Npoisson=G4Poisson(10.0);
00621             if(Npoisson>0) break;
00622           }
00623           G4double tmpPMTgain;
00624           while(1) {
00625             m_PMTgainEc = m_tofSimSvc->EndPMTGain(); 
00626             tmpPMTgain=G4RandGauss::shoot(m_PMTgainEc,m_PMTgainEc/sqrt(Npoisson));
00627             //tmpPMTgain = m_PMTgainEc;
00628             if(tmpPMTgain>0) break;
00629           }
00630           tmpADC[j]+=phtn*tmpPMTgain;
00631 
00632           for (G4int ihst=0; ihst<m_snpeBinNEcV3; ihst++)
00633           {
00634             ii=i+ihst;
00635             if (ii<m_profBinNEcV3)
00636               profPmt[ii][j] += tmpPMTgain*phtn*snpe[ihst];
00637             else
00638               break;
00639           }
00640         }
00641       }
00642 
00643       //add preamplifier and noise
00644       for (int i=0;i<m_profBinNEcV3;i++)
00645       {
00646         if (profPmt[i][j]>0)
00647           profPmt[i][j] = m_preGainEc*profPmt[i][j]+G4RandGauss::shoot(0,m_noiseSigmaEc);
00648       }
00649 
00650       //get pulse height
00651       G4double max=0;
00652       for (int i=n1;i<m_profBinNEcV3;i++)
00653       {
00654         if (profPmt[i][j]>max)
00655           max=profPmt[i][j];
00656       }
00657       if (m_G4Svc->TofRootFlag())
00658       {
00659         if (j==0)  m_max0=max;
00660         else      m_max1=max;
00661       }
00662 
00663 
00664       G4double tmp_HLthresh, tmp_LLthresh, ratio;
00665 
00666 //      G4int runId = m_RealizationSvc->getRunId();
00667 //      if(runId>=-80000 && runId<=-9484)
00668 //      {
00669 //        // After TOF HV adjustment, PMT Gain of barrel TOF in MC is set to 5.E5
00670 //        // High/Low Threshold for barrel: 500/125 mV
00671 //        tmp_HLthresh = m_HLthreshEc;
00672 //        tmp_LLthresh = m_LLthreshEc;
00673 //      }
00674 //      else
00675 //      {
00676 //        // Before TOF HV adjustment, PMT Gain of barrel TOF in MC is set to 2.5E5
00677 //        // High/Low Threshold for barrel: 600/150 mV
00678 //        tmp_HLthresh = 800;
00679 //        tmp_LLthresh = 200;
00680 //      }
00681 
00682       G4double adcFactor = 3.35;
00683 //      double ratio;
00684 //      if (runId>=-80000 && runId<=-9484) 
00685 //      {
00686 //        ratio = 1.22*1354.4/1282.1;
00687 //      }
00688 //      else 
00689 //      {
00690 //        ratio = 1.77*2028.8/1931.4;
00691 //      }
00692  
00693       tmp_HLthresh = m_tofSimSvc->EndHighThres();
00694       tmp_LLthresh = m_tofSimSvc->EndLowThres();
00695       ratio = m_tofSimSvc->EndConstant();
00696      
00697       //get final tdc and adc
00698       if (max>=tmp_HLthresh)
00699       {
00700         for (int i=0;i<m_profBinNEcV3;i++)
00701         {
00702           if ( profPmt[i][j] >= tmp_LLthresh )
00703           {
00704             m_TDC[j] = i*m_timeBinSize + G4RandGauss::shoot(0,0.025); // Adding Electronical Uncertainty of 25ps
00705             G4double NoiseSigma;
00706             G4int EndNoiseSwitch = int(m_tofSimSvc->EndNoiseSwitch());
00707             //std::cout << " EndNoiseSwitch = " << EndNoiseSwitch << std::endl;
00708             switch (EndNoiseSwitch) {
00709               case 0:
00710                 NoiseSigma = 0.;                
00711                 break;
00712               case 1:
00713                 if (partId==0) { NoiseSigma = m_tofSimSvc->EndNoiseSmear(scinNb); }
00714                 if (partId==2) { NoiseSigma = 0.; }
00715                 break;
00716               case 2: 
00717                 if (partId==0) { NoiseSigma = m_tofSimSvc->EndNoiseSmear(scinNb); }
00718                 if (partId==2) { NoiseSigma = m_tofSimSvc->EndNoiseSmear(scinNb+48.); }
00719                 break;      
00720             }
00721             //std::cout << "ID : " << scinNb + 48.*partId/2. << " Sigma : " << NoiseSigma << std::endl;
00722             m_TDC[j] = m_TDC[j] + G4RandGauss::shoot(0,NoiseSigma);   // Adding Endcap Noise Smear
00723             
00724             //use the saturation curve
00725 
00726             if( m_G4Svc->TofSaturationFlag())
00727             {
00728               double x = tmpADC[j]*m_preGainEc*echarge*ratio;
00729               int id;         
00730               if (partId==0) { id = scinNb;}
00731               if (partId==2) { id = scinNb+48;}
00732 
00733               m_ADC[j] = m_tofQElecSvc->EQChannel(id,x);
00734             }
00735             else
00736               m_ADC[j] = tmpADC[j]*m_preGainEc*echarge*adcFactor;
00737 
00738             if (m_G4Svc->TofRootFlag())
00739             {
00740               if (j==0) {
00741                 m_tdc0 = m_TDC[0];
00742                 m_adc0 = m_ADC[0];
00743               }
00744               else     {
00745                 m_tdc1 = m_TDC[1];
00746                 m_adc1 = m_ADC[1];
00747               }
00748             }
00749             break;
00750           }
00751         }
00752       }
00753     }
00754   }
00755 }

G4double BesTofDigitizerEcV3::TransitTime ( )

Definition at line 549 of file BesTofDigitizerEcV3.cc.

References m_ttsMeanEc, and m_ttsSigmaEc.

Referenced by TofPmtAccum().

00550 {
00551   //get time transit spread
00552   //G4cout << "TTS mean:" << m_ttsMeanEc << "; " << m_ttsSigmaEc << G4endl;
00553   return G4RandGauss::shoot(m_ttsMeanEc,m_ttsSigmaEc);
00554 }