BesTofDigitizerBrV2 Class Reference

#include <BesTofDigitizerBrV2.hh>

Inheritance diagram for BesTofDigitizerBrV2:

Public Member Functions
	BesTofDigitizerBrV2 ()
	~BesTofDigitizerBrV2 ()
virtual void	Digitize (ScintSingle *, BesTofDigitsCollection *)
void	ReadData ()
void	TofPmtInit ()
void	TofPmtAccum (BesTofHit *, G4int)
void	DirectPh (G4int, G4ThreeVector, G4double &, G4int &)
G4double	Scintillation (G4int)
G4double	TransitTime ()
void	AccuSignal (G4double, G4int)
void	TofPmtRspns (G4int)
G4double	Reflectivity (G4double n1, G4double n2, G4double theta)
G4double	Reflectivity (G4double n1, G4double n2, G4double n3, G4double theta)
G4double	BirksLaw (BesTofHit *hit)
Public Member Functions inherited from BesTofDigitizerV
	BesTofDigitizerV ()
	~BesTofDigitizerV ()
void	Initialize ()
virtual void	Digitize (ScintSingle *, BesTofDigitsCollection *)

Additional Inherited Members
Protected Attributes inherited from BesTofDigitizerV
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 inherited from BesTofDigitizerV
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

Detailed Description

Definition at line 34 of file BesTofDigitizerBrV2.hh.

Constructor & Destructor Documentation

BesTofDigitizerBrV2()

BesTofDigitizerBrV2::BesTofDigitizerBrV2

(

)

Definition at line 27 of file BesTofDigitizerBrV2.cc.

{
  ReadData();
  m_timeBinSize=0.005;
 
  //retrieve G4Svc
  ISvcLocator* svcLocator = Gaudi::svcLocator();
  IG4Svc* tmpSvc;
  StatusCode sc = svcLocator->service("G4Svc", tmpSvc);
  if(!sc.isSuccess())
  {
    std::cout << " Could not initialize Realization Service in BesTofDigitizerBrV2" << std::endl;
  }
  else
  {
    m_G4Svc = dynamic_cast<G4Svc *>(tmpSvc);
  }
 
  //retrieve RealizationSvc
  IRealizationSvc *tmpReal;
  StatusCode scReal = svcLocator->service("RealizationSvc",tmpReal);
  if (!scReal.isSuccess())
  {
    std::cout << " Could not initialize Realization Service in BesTofDigitizerBrV2" << std::endl;
  } 
  else 
  {
    m_RealizationSvc = dynamic_cast<RealizationSvc*>(tmpReal);
  }
 
 
}

~BesTofDigitizerBrV2()

BesTofDigitizerBrV2::~BesTofDigitizerBrV2

(

)

Definition at line 88 of file BesTofDigitizerBrV2.cc.

{
  ;
}

Member Function Documentation

AccuSignal()

void BesTofDigitizerBrV2::AccuSignal	(	G4double	endTime,
		G4int	forb
	)

Definition at line 746 of file BesTofDigitizerBrV2.cc.

{
  G4int ihst;
  ihst=G4int(endTime/m_timeBinSize);
  if (ihst>0 &&ihst<m_profBinN)
  {
    m_nPhot[ihst][forb]=m_nPhot[ihst][forb]+1;
    m_totalPhot[forb]=m_totalPhot[forb]+1;
  }
}

Referenced by TofPmtAccum().

BirksLaw()

G4double BesTofDigitizerBrV2::BirksLaw

(

BesTofHit *

hit

)

Definition at line 379 of file BesTofDigitizerBrV2.cc.

{
  const G4double kappa = 0.015*cm/MeV;
  const G4String brMaterial = "BC408";
  G4double dE = hit->GetEdep();
  //G4cout << "The edep is "<< dE << G4endl;
  G4double dX = hit->GetStepL();
  //G4Material* materiral = hit->GetMaterial();
  G4double charge = hit->GetCharge();
  G4double cor_dE = dE;
  //if((materiral->GetName()==brMaterial) && charge!=0.&& dX!=0.)
  if(charge!=0.&& dX!=0.)
  {
    cor_dE = dE/(1+kappa*dE/dX);
    //if(dE>20)
    //{
    //  G4cout << "\n dE > 20. Details are below:" << G4endl;
    //  G4cout << "dE/dx:" << dE/dX << G4endl;
    //  G4cout << "dE:" << dE << "; dX:" << dX << G4endl;
    //  G4cout << "It is BC408. cor_dE is " << cor_dE << G4endl;
    //  G4double ratio = cor_dE/dE;
    //  G4cout << "The ratio cor_dE/edep is "<< ratio << G4endl;
    //}
    //G4cout << "It is BC408. cor_dE is " << cor_dE << G4endl;
    //G4double ratio = cor_dE/dE;
    //G4cout << "The ratio cor_dE/edep is "<< ratio << G4endl;
  }
  return cor_dE;
 
}

Referenced by TofPmtAccum().

Digitize()

void BesTofDigitizerBrV2::Digitize ( ScintSingle *  scint, BesTofDigitsCollection *  DC  )

virtual

Reimplemented from BesTofDigitizerV.

Definition at line 93 of file BesTofDigitizerBrV2.cc.

{
  m_beamTime = m_G4Svc->GetBeamTime() * ns;
  m_besTofDigitsCollection = DC;
 
  G4DigiManager* digiManager = G4DigiManager::GetDMpointer();
  G4int THCID = digiManager->GetHitsCollectionID("BesTofHitsCollection");
  m_THC = (BesTofHitsCollection*) (digiManager->GetHitsCollection(THCID));
 
  if (m_G4Svc->TofRootFlag())
  {
    m_eTotal = 0;
    m_nDigi = 0;
    m_partIdMPV = -9;
    m_scinNbMPV = -9;
    m_edepMPV = 0;
    m_nDigiOut = 0;
  }
 
  if (m_THC)
  {
    //for each digi, compute TDC and ADC
    G4int partId, scinNb, nHits;
    G4double edep;
    BesTofHit* hit;
    partId=scint->GetPartId();
    scinNb=scint->GetScinNb();
    edep = scint->GetEdep();
    nHits=scint->GetHitIndexes()->size();
 
 
    //std::cout << "BesTofDigitizerBrV2      Partid   scinNb   " << partId  << "   "  << scinNb << std::endl; 
    //cout << "*** scinNb:" << scinNb << " *** " << m_tofCaliSvc->BAtten(scinNb) << "***" << endl;
    //cout << "*****scinNb:"<< scinNb << " ***** A2:" << m_tofCaliSvc->BGainBackward(scinNb) 
    //  << " ***** A1:" << m_tofCaliSvc->BGainForward(scinNb) << endl;
 
    TofPmtInit();
 
    //fill tof Ntuple
    if (m_G4Svc->TofRootFlag())
    {
      if (edep>m_edepMPV)
      {
        m_partIdMPV = partId;
        m_scinNbMPV = scinNb;
        m_edepMPV = edep;
      }
      m_eTotal += edep;
      m_nDigi ++;
 
      m_partId = partId;
      m_scinNb = scinNb;
      m_edep = edep;
      m_nHits = nHits;
    }
 
    if (edep>0.01)
    {
      for (G4int j=0;j<nHits;j++)
      {
        hit= (*m_THC)[( *(scint->GetHitIndexes()) )[j]];
        TofPmtAccum(hit, scinNb);
      }
      if (m_G4Svc->TofRootFlag())
      {
        m_time1st0=m_t1st[0];
        m_time1st1=m_t1st[1];
        m_timelast0=m_tLast[0];
        m_timelast1=m_tLast[1];
        m_totalPhot0=m_totalPhot[0];
        m_totalPhot1=m_totalPhot[1];
      }
      //get final tdc and adc
      TofPmtRspns(scinNb);
      //G4cout<<"pre-cut " << partId << "\nadc0:"<<m_ADC[0]<<"; adc1:"
      //  <<m_ADC[1]<<"; tdc0:"<<m_TDC[0]<<"; tdc1:"
      //  <<m_TDC[1]<<G4endl;
 
      G4double temp0 = m_ADC[0]+m_TDC[0];
      G4double temp1 = m_ADC[1]+m_TDC[1];
      //cout << "partid: " <<  partId << " temp0: " <<  temp0 << " temp1: " <<  temp1 << endl;
      //if ( partId==1 && m_ADC[0]>255 && m_ADC[1]>255 && m_TDC[0]>0. && m_TDC[1]>0.)
      if ( partId==1 && (temp0 > -1998. || temp1 > -1998.))
      {
        //const double MAX_ADC = 8191 * 0.3;  // channel set up to 8192 will lead to overflow.
        BesTofDigi* digi = new BesTofDigi;
        digi->SetTrackIndex(m_trackIndex);
        digi->SetPartId(partId);
        digi->SetScinNb(scinNb);
        digi->SetForwADC( m_ADC[0]) ;
        digi->SetBackADC( m_ADC[1]) ;
        if (m_TDC[0]>0)
          m_TDC[0] = m_TDC[0]+m_beamTime;
        if (m_TDC[1]>0)
          m_TDC[1] = m_TDC[1]+m_beamTime;
        digi->SetForwTDC( m_TDC[0]) ;
        digi->SetBackTDC( m_TDC[1]) ;
        //G4cout<<"+++++++++++++++++++++++++++++barrel\nadc0:"<<m_ADC[0]<<"; adc1:"
        //  <<m_ADC[1]<<"; tdc0:"<<m_TDC[0]<<"; tdc1:"
        //  <<m_TDC[1]<<G4endl;
 
        m_besTofDigitsCollection->insert(digi);
 
        if (m_G4Svc->TofRootFlag())
          m_nDigiOut++;
 
      }
      if (m_G4Svc->TofRootFlag())
        m_tupleTof1->write();
 
    }
 
    if (m_G4Svc->TofRootFlag())
      m_tupleTof2->write();
  }
}

Referenced by BesTofDigitizer::Digitize().

DirectPh()

void BesTofDigitizerBrV2::DirectPh	(	G4int	partId,
		G4ThreeVector	emtPos,
		G4double &	pathL,
		G4int &	forb
	)

Definition at line 457 of file BesTofDigitizerBrV2.cc.

{
  //std::cout << "BesTofDigitizerBrV2::DirectPh()" << std::endl;
  const static G4double silicon_oil_index = 1.43; 
  const static G4double glass_index = 1.532;
  //generation photon have random direction
  //optical parameters of scintillation simulation
  G4double cos_span=1-cos( asin(silicon_oil_index/m_refIndex) );
  //G4double cos_spanEc = 1;
  G4double dcos, ran;
  ran=G4UniformRand();
  //assuming uniform phi distribution, simulate cos distr only
  dcos=cos_span*(ran*2.0-1.0);
  //G4double dcosEc = cos_spanEc*(ran*2.0-1.0);
  G4double r2=sqrt(emtPos.x()*emtPos.x()+emtPos.y()*emtPos.y());
  G4int nRef=0;
  G4double costheta,sintheta;
  G4double theta,thetaR; // thetaR is scin to air full ref angle. about 39.265 degree.
  costheta=dcos>0?(1-dcos):(1+dcos);
  theta=acos(costheta);
  sintheta=sin(theta);
  thetaR=asin(1/m_refIndex);
  G4double R1;
  R1=Reflectivity(m_refIndex,1.0,theta);
  G4double ratio1Mean=(CLHEP::pi)*25.5*25.5/(57.1+60.7)/25.0; //0.693657
  G4double ratio2Mean=(19.5/25.5)*(19.5/25.5); //0.584775
  G4double ratio1=G4RandGauss::shoot(ratio1Mean,0.015);
  G4double ratio2=G4RandGauss::shoot(ratio2Mean,0.015);
  G4double ratio3Mean=0.945*ratio1Mean*ratio2Mean;
  G4double ratio3=G4RandGauss::shoot(ratio3Mean,0.015);
 
  G4double R2=1-Reflectivity(m_refIndex,silicon_oil_index, glass_index, theta);
  if (dcos > 0 && dcos != 1)
  {
    if (theta < thetaR) // theta < 39.265 degree
    {
      if (G4UniformRand()<ratio1) //coup1
      {
        if (G4UniformRand()<R2)
        {
          if (G4UniformRand()<ratio2) //PMT 39mm
          {
            forb=0;
            pathL=(m_scinLength/2-emtPos.z())/costheta;
          }
        }
        else
        {
          if (G4UniformRand()<ratio3)
          {
            forb=1;
            pathL=(1.5*m_scinLength-emtPos.z())/costheta;
          }
        }
      }
      else //Air
      {
        if (G4UniformRand()<R1)
        {
          G4double tempran=G4UniformRand();
          if (tempran<ratio3)
          {
            forb=1;
            pathL=(1.5*m_scinLength-emtPos.z())/costheta;
          }
          else if (tempran>ratio1&&G4UniformRand()<R1&&G4UniformRand()<ratio3)
          {
            forb=0;
            pathL=(2.5*m_scinLength-emtPos.z())/costheta;
          }
        }
      }
    }
    else  // 39.265 <= theta < 64.832
    {
      if (G4UniformRand()<ratio1) //coup1
      {
        if (G4UniformRand()<R2)
        {
          if (G4UniformRand()<ratio2) //PMT 39mm
          {
            forb=0;
            pathL=(m_scinLength/2-emtPos.z())/costheta;
          }
        }
        else
        {
          if (G4UniformRand()<ratio3)
          {
            forb=1;
            pathL=(1.5*m_scinLength-emtPos.z())/costheta;
          }
        }
      }
      else //Air
      {
        G4double tempran=G4UniformRand();
        if (tempran<ratio3)
        {
          forb=1;
          pathL=(1.5*m_scinLength-emtPos.z())/costheta;
        }
        else if (tempran>ratio1 && G4UniformRand()<ratio3)
        {
          forb=0;
          pathL=(2.5*m_scinLength-emtPos.z())/costheta;
        }
      }
    }
  }
  else if (dcos < 0 && dcos != -1)
  {
    if (theta < thetaR) // theta < 39.265 degree
    {
      if (G4UniformRand()<ratio1) //coup1
      {
        if (G4UniformRand()<R2)
        {
          if (G4UniformRand()<ratio2) //PMT 39mm
          {
            forb=1;
            pathL=(m_scinLength/2+emtPos.z())/costheta;
          }
        }
        else
        {
          if (G4UniformRand()<ratio3)
          {
            forb=0;
            pathL=(1.5*m_scinLength+emtPos.z())/costheta;
          }
        }
      }
      else //Air
      {
        if (G4UniformRand()<R1)
        {
          G4double tempran=G4UniformRand();
          if (tempran<ratio3)
          {
            forb=0;
            pathL=(1.5*m_scinLength+emtPos.z())/costheta;
          }
          else if (tempran>ratio1&&G4UniformRand()<R1&&G4UniformRand()<ratio3)
          {
            forb=1;
            pathL=(2.5*m_scinLength+emtPos.z())/costheta;
          }
        }
      }
    }
    else  // 39.265 <= theta < 64.832
    {
      if (G4UniformRand()<ratio1) //coup1
      {
        if (G4UniformRand()<R2)
        {
          if (G4UniformRand()<ratio2) //PMT 39mm
          {
            forb=1;
            pathL=(m_scinLength/2+emtPos.z())/costheta;
          }
        }
        else
        {
          if (G4UniformRand()<ratio3)
          {
            forb=0;
            pathL=(1.5*m_scinLength+emtPos.z())/costheta;
          }
        }
      }
      else //Air
      {
        G4double tempran=G4UniformRand();
        if (tempran<ratio3)
        {
          forb=0;
          pathL=(1.5*m_scinLength+emtPos.z())/costheta;
        }
        else if (tempran>ratio1 && G4UniformRand()<ratio3)
        {
          forb=1;
          pathL=(2.5*m_scinLength+emtPos.z())/costheta;
        }
      }
    }
  }
 
  G4double convFactor=180./3.1415926;
  if (theta>asin(1)-thetaR)
  {
    G4double alpha = G4UniformRand()*asin(1.);
    G4int nRef1 = pathL*sintheta*cos(alpha)/50.0+0.5;
    G4int nRef2 = pathL*sintheta*sin(alpha)/58.9+0.5;
    G4double beta1=acos(sintheta*cos(alpha));
    G4double beta2=acos(sintheta*sin(alpha));
    beta2 -= 3.75/convFactor;
    G4double R21,R22;
    R21=Reflectivity(m_refIndex,1.0,beta1);
    R22=Reflectivity(m_refIndex,1.0,beta2);
    for (G4int i=0;i<nRef1;i++)
    {
      if (G4UniformRand()<(1-R21) && G4UniformRand()<0.15)
        pathL=-9;
    }
    for (G4int i=0;i<nRef2;i++)
    {
      if (G4UniformRand()<(1-R22) && G4UniformRand()<0.15)
        pathL=-9;
    }
  }
}

Referenced by TofPmtAccum().

ReadData()

void BesTofDigitizerBrV2::ReadData

(

)

Definition at line 60 of file BesTofDigitizerBrV2.cc.

{
  BesTofGeoParameter* tofPara = BesTofGeoParameter::GetInstance();
 
  m_scinLength = tofPara->GetBr1L();
  m_tau1       = tofPara->GetTau1();
  m_tau2       = tofPara->GetTau2();
  m_tau3       = tofPara->GetTau3();
  m_tauRatio   = tofPara->GetTauRatio();
  m_refIndex   = tofPara->GetRefIndex();
  m_phNConst   = tofPara->GetPhNConst();
  m_Cpe2pmt    = tofPara->GetCpe2pmt();
  m_rAngle     = tofPara->GetRAngle();
  m_QE         = tofPara->GetQE();
  m_CE         = tofPara->GetCE();
  m_peCorFac   = tofPara->GetPeCorFac();
 
  m_ttsMean    = tofPara->GetTTSmean();
  m_ttsSigma   = tofPara->GetTTSsigma();
  m_Ce         = tofPara->GetCe();
  m_LLthresh   = tofPara->GetLLthresh();
  m_HLthresh   = tofPara->GetHLthresh();
  m_preGain    = tofPara->GetPreGain();
  m_noiseSigma = tofPara->GetNoiseSigma();
 
 
}

Referenced by BesTofDigitizerBrV2().

Reflectivity() [1/2]

G4double BesTofDigitizerBrV2::Reflectivity	(	G4double	n1,
		G4double	n2,
		G4double	n3,
		G4double	theta
	)

Definition at line 410 of file BesTofDigitizerBrV2.cc.

{
  G4double I1,I2,I3,rp1,rs1,rp2,rs2,Rp1,Rs1,Rp2,Rs2,Rp,Rs;
  G4double R=1.0;
  //n1=m_refIndex;
  //n2=1.0;
  I1=theta;
  if(I1<asin(n2/n1))
  {
    I2=asin(sin(I1)*(n1/n2));
    rp1=(n1/cos(I1)-n2/cos(I2))/(n1/cos(I1)+n2/cos(I2));
    rs1=(n1*cos(I1)-n2*cos(I2))/(n1*cos(I1)+n2*cos(I2));
    Rp1=rp1*rp1;
    Rs1=rs1*rs1;
 
    I3=asin(sin(I1)*(n1/n3));
    rp2=(n2/cos(I2)-n3/cos(I3))/(n2/cos(I2)+n3/cos(I3));
    rs2=(n2*cos(I2)-n3*cos(I3))/(n2*cos(I2)+n3*cos(I3));
    Rp2=rp2*rp2;
    Rs2=rs2*rs2;
    Rp=(Rp1+Rp2-2*Rp1*Rp2)/(1-Rp1*Rp2);
    Rs=(Rs1+Rs2-2*Rs1*Rs2)/(1-Rs1*Rs2);
    R=(Rp+Rs)/2.;
  }
  return R;
}

Reflectivity() [2/2]

G4double BesTofDigitizerBrV2::Reflectivity	(	G4double	n1,
		G4double	n2,
		G4double	theta
	)

Definition at line 438 of file BesTofDigitizerBrV2.cc.

{
  G4double I1,I2,rp,rs,Rp,Rs;
  G4double R=1.0;
  //n1=m_refIndex;
  //n2=1.0;
  I1=theta;
  if (I1<asin(n2/n1))
  {
    I2=asin(sin(I1)*(n1/n2));
    rp=(n1/cos(I1)-n2/cos(I2))/(n1/cos(I1)+n2/cos(I2));
    rs=(n1*cos(I1)-n2*cos(I2))/(n1*cos(I1)+n2*cos(I2));
    Rp=rp*rp;
    Rs=rs*rs;
    R=(Rp+Rs)/2.;
  }
  return R;
}

Referenced by DirectPh().

Scintillation()

G4double BesTofDigitizerBrV2::Scintillation

(

G4int

partId

)

Definition at line 693 of file BesTofDigitizerBrV2.cc.

{
  G4double tmp_tauRatio,tmp_tau1,tmp_tau2,tmp_tau3;
  tmp_tauRatio = m_tauRatio;
  tmp_tau1 = m_tau1;
  tmp_tau2 = m_tau2;
  tmp_tau3 = m_tau3;
  G4double UniformR = tmp_tauRatio/(1+tmp_tauRatio);
  G4double EmissionTime;
  if (G4UniformRand()>UniformR) {
    while (1) {
      EmissionTime = -tmp_tau2*log( G4UniformRand() );
      if (G4UniformRand()-exp(EmissionTime/tmp_tau2-EmissionTime/tmp_tau1)>1.E-8)
        break;
    }
  }
  else EmissionTime = -tmp_tau3*log( G4UniformRand() );
  return EmissionTime;
}

Referenced by TofPmtAccum().

TofPmtAccum()

void BesTofDigitizerBrV2::TofPmtAccum	(	BesTofHit *	hit,
		G4int	scinNb
	)

Definition at line 246 of file BesTofDigitizerBrV2.cc.

{
  BesTofGeoParameter* tofPara = BesTofGeoParameter::GetInstance();
  //std::cout << "BesTofDigitizerBrV2::TofPmtAccum()" << std::endl;
  G4double cvelScint=c_light/m_refIndex/1.07;
  //Get information of this step
  G4ThreeVector pos=hit->GetPos();
  G4int trackIndex = hit->GetTrackIndex();
  G4int partId =hit->GetPartId();
  G4double edep=hit->GetEdep();
  G4double stepL=hit->GetStepL();
  G4double deltaT=hit->GetDeltaT();
  G4double timeFlight=hit->GetTime()-m_beamTime;
  //std::cout << "timeFlight: " << timeFlight << std::endl;
  if (timeFlight < m_globalTime)
  {
    m_globalTime = timeFlight;
    m_trackIndex = trackIndex;
  }
 
  G4ThreeVector pDirection=hit->GetPDirection();
  G4double nx=pDirection.x();
  G4double ny=pDirection.y();
  G4double nz=pDirection.z();
 
  //phNConst=(Nph/MeV)*(QE)*(CE)*(1-cos(crit));
  //        =10000 * 0.2 * 0.6 * (1-cos(39.265))=270.931
  //asin(1/1.58) = 39.265248
  //only the light has theta=0---39.265 can go out to PMT, the probability is computed with solid angle
  //solid angle = phi*(1-cos(theta)), phi=2pi
 
  //Cpe2pmt=cathode area/scint area
  //peCorFac=correction factor for eff. of available Npe
 
  //G4double thetaProb = 1-sqrt(m_refIndex*m_refIndex-1)/m_refIndex;
  G4double thetaProb=1-cos( asin(1.43/m_refIndex));
  //G4double thetaProbEc = 1-1/m_refIndex;
 
  //number of photons generated in this step 
  G4double nMean, nPhoton;
  //std::cout << "0 BirksLaw(): " << std::endl;
  nMean = m_phNConst*BirksLaw(hit);
  //std::cout << "1 BirksLaw(): " << std::endl;
 
  if(nMean>10)
  {
    G4double resolutionScale=1.;
    G4double sigma=resolutionScale*sqrt(nMean);
    nPhoton=G4int(G4RandGauss::shoot(nMean,sigma)+0.5);
  }
  else
    nPhoton=G4int(G4Poisson(nMean));
  //G4cout<<"nPhoton:"<<nPhoton<<G4endl;
 
  G4int NphStep;
  if (partId==1)
    NphStep=G4int(nPhoton*thetaProb*m_rAngle*m_QE*m_CE*m_peCorFac);
  //else
  //NphStep=G4int(edep*m_phNConstEc*m_Cpe2pmtEc*0.66*m_QEEc*m_CE*m_peCorFac);
  //introduce poission distribution of Npe
  //G4double Navr = NphStep;
  //G4double NpePoiss = G4double(G4Poisson(Navr));
  //G4double adcW_corr =5.0;
  //NphStep=G4int( (NpePoiss - Navr)*adcW_corr + Navr );
 
  if (m_G4Svc->TofRootFlag())
    m_NphAllSteps +=  NphStep;
  //std::cout << "m_G4Svc->TofRootFlag(): " << m_G4Svc->TofRootFlag() << std::endl;
 
  G4double ddS, ddT;
  if (NphStep>0)
  {
    for (G4int i=0;i<NphStep;i++)
    {
      //uniform scintilation in each step
      ddS=stepL*G4UniformRand();
      ddT=deltaT*G4UniformRand();
      G4ThreeVector emtPos;
      emtPos.setX(pos.x() + nx*ddS);
      emtPos.setY(pos.y() + ny*ddS);
      emtPos.setZ(pos.z() + nz*ddS);
 
      //check scinillation light whether to hit the pmt or not
      //forb=0/1 for forward(z>0, east) and backward(z<0, west)
      G4double pathL=0;
      G4int forb;
      DirectPh(partId, emtPos, pathL, forb);
 
 
      //check if photon can reach PMT or not, after attenuation
 
      G4double ran = G4UniformRand();
      //G4double attenL = tofPara->GetAtten(scinNb);
      //G4double attenL = 10.*(m_tofCaliSvc->BAtten(scinNb))/0.75; // cm into mm
      G4double attenL = m_tofSimSvc->BarAttenLength(scinNb); 
      attenL = 10.*attenL/0.75; // cm into mm
 
      if (pathL>0 && exp(-pathL/attenL) > ran)
      {
        //propagation time in scintillator
        G4double scinSwim=pathL/cvelScint;
        //scintillation timing
        //G4double scinTime=GenPhoton(partId);
        G4double scinTime=Scintillation(partId);
 
        //PMT transit time
        G4double transitTime=TransitTime();
        //sum up all time components
        G4double endTime= timeFlight + ddT + scinSwim + scinTime + transitTime;
        //std::cout << "endtime: " << endTime << std::endl;
 
        if (m_G4Svc->TofRootFlag())
        {
          //m_forb = forb;
          //m_timeFlight = timeFlight;
          //m_ddT = ddT;
          m_scinSwim = scinSwim;
          //m_scinTime = scinTime;
          //m_transitTime = transitTime;
          m_endTime = endTime;
          m_tupleTof3->write();
        }
        //store timing into binned buffer
        AccuSignal(endTime, forb);
 
        //update 1st and last timings here
        if (m_t1st[forb]>endTime)   m_t1st[forb] = endTime;
        if (m_tLast[forb]<endTime)  m_tLast[forb]= endTime;
      }
    }
  }
}

Referenced by Digitize().

TofPmtInit()

void BesTofDigitizerBrV2::TofPmtInit

(

)

Definition at line 210 of file BesTofDigitizerBrV2.cc.

{
  Initialize();
 
  m_t1st[0]=100;
  m_t1st[1]=100;
  m_tLast[0]=0.;
  m_tLast[1]=0;
  m_totalPhot[0]=0;
  m_totalPhot[1]=0;
  for (G4int i=0;i<2;i++)
    for (G4int j=0;j<m_profBinN;j++)
      m_nPhot[j][i]=0;
 
  if (m_G4Svc->TofRootFlag())
  {
    m_partId = -9;
    m_scinNb = -9;
    m_edep = 0;
    m_nHits = 0;
    m_time1st0 = 100;
    m_time1st1 = 100;
    m_timelast0 = 0;
    m_timelast1 = 0;
    m_totalPhot0 = 0;
    m_totalPhot1 = 0;
    m_NphAllSteps = 0;
    m_max0 = 0;
    m_max1 = 0;
    m_tdc0 = -999;
    m_adc0 = -999;
    m_tdc1 = -999;
    m_adc1 = -999;
  }
}

Referenced by Digitize().

TofPmtRspns()

void BesTofDigitizerBrV2::TofPmtRspns

(

G4int

scinNb

)

Definition at line 757 of file BesTofDigitizerBrV2.cc.

{
  //std::cout << "BesTofDigitizerBrV2::TofPmtRspns()" << std::endl;
  BesTofGeoParameter* tofPara = BesTofGeoParameter::GetInstance();
  //to generate PMT signal shape for single photoelectron.
  //only one time for a job.
  G4double snpe[m_snpeBinN][2];
 
  //Model: f(t)=Gain*mv_1pe* t**2 * exp-(t/tau)**2/normal-const
  //normalization const =sqrt(pi)*tau*tau*tau/4.0
  //G4double tau = m_riseTime;
  G4double norma_const;
  G4double echarge=1.6e-7;   //in unit of PC
 
  //time profile of pmt signals for Back and Forward PMT.
  G4double profPmt[m_profBinN][2];
 
  G4double t;
  G4double tau;
  G4int n1, n2, ii;
  G4int phtn;
  // forb = 0, east
  for (G4int i=0;i<m_snpeBinN;i++)
  {
    tau = tofPara->GetBrERiseTime(scinNb);
    norma_const = sqrt(CLHEP::pi)*tau*tau*tau/4.0;    //in unit of ns**3
    t = (i+1)*m_timeBinSize;
    snpe[i][0] = m_Ce*t*t*exp(- (t/tau) * (t/tau) )/norma_const;// Pulse of one single photoelectron
  }
  // forb = 1, west
  for (G4int i=0;i<m_snpeBinN;i++)
  {
    tau = tofPara->GetBrWRiseTime(scinNb);
    norma_const = sqrt(CLHEP::pi)*tau*tau*tau/4.0;    //in unit of ns**3
    t = (i+1)*m_timeBinSize;
    snpe[i][1] = m_Ce*t*t*exp(- (t/tau) * (t/tau) )/norma_const;
  }
  //for barrel and endcap tof: fb=2 or 1
  G4int fb=2;
  G4double Npoisson;
  G4double pmtGain0,pmtGain,relativeGain,smearPMTgain;
  G4double smearADC[2] = {0.,0.};
  G4int runId = m_RealizationSvc->getRunId();
  pmtGain0 = m_tofSimSvc->BarPMTGain();
 
//  if(runId>=-80000 && runId<=-9484)
//  {
//    // After TOF HV adjustment, PMT Gain of barrel TOF in MC is set to 5.E5
//    // High/Low Threshold for barrel: 500/125 mV
//    pmtGain0 = 6.E5;
//  }
//  else
//  {
//    // Before TOF HV adjustment, PMT Gain of barrel TOF in MC is set to 2.5E5
//    // High/Low Threshold for barrel: 600/150 mV
//    pmtGain0 = 5.E5;
//  }
 
  G4double timeSmear = G4RandGauss::shoot(0,0.020);
  if (runId>=-22913 && runId<=-20448) {//for 2011-psipp(20448-22913), smear barrel TOF resolution to ~78ps
     timeSmear = G4RandGauss::shoot(0,0.040);
  }
 
  for (G4int j=0; j<fb; j++)
  {
    if (m_totalPhot[j] > 0)
    {
      n1=G4int(m_t1st[j]/m_timeBinSize);
      n2=G4int(m_tLast[j]/m_timeBinSize);
      //std::cout << "n1: " << n1 << std::endl;
 
      for (G4int i=0;i<m_profBinN;i++)
        profPmt[i][j]=0.0;
 
      //generate PMT pulse
      n2 = n2<m_profBinN ? n2:m_profBinN;
      for (G4int i=n1;i<n2;i++)
      {
        phtn=m_nPhot[i][j];
        if (phtn>0)
        {
          while(1) {
            Npoisson = G4Poisson(10.0);
            if(Npoisson>0.) break;
          }
          while(1) {
            //pmtGain = j ? tofPara->GetBrWPMTgain(scinNb) : tofPara->GetBrEPMTgain(scinNb);
            //relativeGain = j ? m_tofCaliSvc->BGainBackward(scinNb) : m_tofCaliSvc->BGainForward(scinNb);
            relativeGain = j ? m_tofSimSvc->BarGain2(scinNb) : m_tofSimSvc->BarGain1(scinNb);
        pmtGain = pmtGain0*relativeGain;
            smearPMTgain = G4RandGauss::shoot(pmtGain,pmtGain/sqrt(Npoisson));
            //smearPMTgain = pmtGain;
            if(smearPMTgain>0) break;
          }
          smearADC[j] += phtn*smearPMTgain;
 
          for (G4int ihst=0; ihst<m_snpeBinN; ihst++)
          {
            ii=i+ihst;
            if (ii<m_profBinN)
              profPmt[ii][j] += smearPMTgain*phtn*snpe[ihst][j];
            else
              break;
          }
        }
      }
 
      //add preamplifier and noise
      for (int i=0;i<m_profBinN;i++)
      {
        if (profPmt[i][j]>0)
          profPmt[i][j] = m_preGain*profPmt[i][j]+G4RandGauss::shoot(0,m_noiseSigma);
      }
 
      //get pulse height
      G4double max=0;
      for (int i=n1;i<m_profBinN;i++)
      {
        if (profPmt[i][j]>max)
          max=profPmt[i][j];
      }
      if (m_G4Svc->TofRootFlag())
      {
        if (j==0)  m_max0=max;
        else      m_max1=max;
      }
 
      G4double tmp_HLthresh, tmp_LLthresh, adcFactor;
      G4double ratio;
 
      if(runId>=-80000 && runId<=-9484)
      {
        // After TOF HV adjustment, PMT Gain of barrel TOF in MC is set to 5.E5
        // High/Low Threshold for barrel: 500/125 mV
        tmp_HLthresh = m_HLthresh;
        tmp_LLthresh = m_LLthresh;
        adcFactor = 5.89;
      }
      else
      {
        // Before TOF HV adjustment, PMT Gain of barrel TOF in MC is set to 2.5E5
        // High/Low Threshold for barrel: 600/150 mV
        tmp_HLthresh = 600;
        tmp_LLthresh = 150;
        adcFactor = 4.8;
      }
//      if(runId>=-80000 && runId<=-9484)
//      {
//        ratio=2.16*1923.8/2197.8;
//      }
//      else
//      {
//        ratio = 2.11*2437.0/3102.9;
//      }
      tmp_HLthresh = m_tofSimSvc->BarHighThres();
      tmp_LLthresh = m_tofSimSvc->BarLowThres();
      ratio = m_tofSimSvc->BarConstant();
 
      //get final tdc and adc
      if (max>=tmp_HLthresh)
      {
        for (int i=0;i<m_profBinN;i++)
        {
          if ( profPmt[i][j] >= tmp_LLthresh )
          {
            m_TDC[j] = i*m_timeBinSize + timeSmear +G4RandGauss::shoot(0,0.025); // Adding Electronical Uncertainty of 25ps
            
            if( m_G4Svc->TofSaturationFlag())
            {
              //get ADC[j] using tofQElecSvc
              double x = m_preGain*smearADC[j]*echarge*ratio;
              if (j==0) 
              { 
                m_ADC[j] = m_tofQElecSvc->BQChannel1(scinNb,x);
              }
              else 
              {
                m_ADC[j] = m_tofQElecSvc->BQChannel2(scinNb,x);
              }
            }
            else
              m_ADC[j] = m_preGain*smearADC[j]*echarge*adcFactor;
 
 
            if (m_G4Svc->TofRootFlag())
            {
              if (j==0) {
                m_tdc0 = m_TDC[0];
                m_adc0 = m_ADC[0];
              }
              else     {
                m_tdc1 = m_TDC[1];
                m_adc1 = m_ADC[1];
              }
            }
            break;
          }
        }
      }
    }
  }
}

Referenced by Digitize().

TransitTime()

G4double BesTofDigitizerBrV2::TransitTime

(

)

Definition at line 740 of file BesTofDigitizerBrV2.cc.

{
  //get time transit spread
  return G4RandGauss::shoot(m_ttsMean,m_ttsSigma);
}

Referenced by TofPmtAccum().

---

The documentation for this class was generated from the following files:

• BesTofDigitizerBrV2.hh
• BesTofDigitizerBrV2.cc