BesMdcSD Class Reference

#include <BesMdcSD.hh>

Inheritance diagram for BesMdcSD:

List of all members.

Public Member Functions
	BesMdcSD (G4String)
	~BesMdcSD ()
void	Initialize (G4HCofThisEvent *)
G4bool	ProcessHits (G4Step *, G4TouchableHistory *)
void	EndOfEvent (G4HCofThisEvent *)
void	BeginOfTruthEvent (const G4Event *)
void	EndOfTruthEvent (const G4Event *)
G4double	Distance (G4int, G4int, G4ThreeVector, G4ThreeVector, G4ThreeVector &, G4double &)
void	dedxFuncInti (void)
virtual void	BeginOfTrack (const G4Track *)
virtual void	EndOfTrack (const G4Track *)
Protected Member Functions
void	GetCurrentTrackIndex (G4int &trackIndex, G4int &g4TrackId) const
Private Member Functions
G4int	GetBetagammaIndex (G4double bg)
G4int	GetAngleIndex (G4double)
G4int	GetChargeIndex (G4int)
G4double	GetValDedxCurve (G4double bg, G4double charge)
G4double	dedxSample (G4double betagamma, G4double length, G4double theta)
Private Attributes
G4int	hitPointer [43][288]
G4int	truthPointer [43][288]
BesMdcHitsCollection *	hitsCollection
BesMdcHitsCollection *	truthCollection
BesMdcGeoParameter *	mdcGeoPointer
BesMdcCalTransfer *	mdcCalPointer
MdcGeomSvc *	mdcGeomSvc
G4Svc *	m_G4Svc
TF1 *	dEdE_mylanfunc
CalibDataSvc *	m_calibDataSvc
	dedx sim ---------------------------
IDedxCurSvc *	m_pDedxCurSvc
TH1F *	m_dedx_hists
G4int	m_version
G4int	m_numDedxHists
G4int	m_numBg
G4int	m_numTheta
std::vector< G4double >	m_bgRange
NTuple::Tuple *	m_tupleMdc
NTuple::Item< double >	m_betaGamma
NTuple::Item< double >	m_fitval
NTuple::Item< double >	m_random
NTuple::Item< double >	m_dedx
NTuple::Item< double >	m_de
NTuple::Item< double >	m_charge
NTuple::Item< double >	m_costheta

---

Detailed Description

Definition at line 35 of file BesMdcSD.hh.

---

Constructor & Destructor Documentation

BesMdcSD::BesMdcSD

(

G4String

)

Definition at line 29 of file BesMdcSD.cc.

References dedxFuncInti(), BesMdcGeoParameter::GetGeo(), G4Svc::GetMdcDedxFlag(), G4Svc::GetTupleMdc(), genRecEmupikp::i, m_betaGamma, m_bgRange, m_calibDataSvc, m_charge, m_costheta, m_de, m_dedx, m_dedx_hists, m_fitval, m_G4Svc, m_numBg, m_numDedxHists, m_numTheta, m_pDedxCurSvc, m_tupleMdc, m_version, mdcCalPointer, mdcGeomSvc, mdcGeoPointer, G4Svc::MdcRootFlag(), and deljobs::string.

:BesSensitiveDetector(name)
{
    collectionName.insert("BesMdcHitsCollection");
    collectionName.insert("BesMdcTruthCollection");
  
    mdcGeoPointer=BesMdcGeoParameter::GetGeo();
    mdcCalPointer=new BesMdcCalTransfer;
  
    IMdcGeomSvc* ISvc;
    StatusCode sc=Gaudi::svcLocator()->service("MdcGeomSvc", ISvc);
    if (!sc.isSuccess()) 
      std::cout<<"BesMdcSD::Could not open MdcGeomSvc"<<std::endl;
    mdcGeomSvc=dynamic_cast<MdcGeomSvc *>(ISvc);
  
    IG4Svc* tmpSvc;
    sc=Gaudi::svcLocator()->service("G4Svc", tmpSvc);
    if (!sc.isSuccess())
      G4cout <<" MdcSD::Error,could not open G4Svc"<<G4endl;
    m_G4Svc=dynamic_cast<G4Svc *>(tmpSvc); 
  
    if(m_G4Svc->GetMdcDedxFlag()==1){
      G4cout <<" MdcSD: Use sampled dedx instead of Geant4 value"<<G4endl;
      dedxFuncInti();
    }

    
    //get DedxSimData
    std::string dedxTDSPath = "/Calib/DedxSim";
    IDataProviderSvc* pCalibDataSvc;
    sc = Gaudi::svcLocator()->service("CalibDataSvc", pCalibDataSvc, true);
    if (!sc.isSuccess())
        std::cout << "BesMdcSD::Could not open CalibDataSvc" << std::endl;
    m_calibDataSvc = dynamic_cast<CalibDataSvc*>(pCalibDataSvc);
    if (!sc.isSuccess())
    {
        std::cout << "Could not get CalibDataSvc"
                  << std::endl;
    }
    SmartDataPtr<CalibData::DedxSimData> pDedxSimData(m_calibDataSvc, dedxTDSPath);
    m_version = pDedxSimData->getVersion();
    m_numDedxHists = pDedxSimData->gethistNo();
    m_numBg = pDedxSimData->getRangeNo();
    if (m_version == 0) m_numTheta = 10;
    else m_numTheta = pDedxSimData->getThetaNo();
    m_dedx_hists = new TH1F[m_numDedxHists];
    for (G4int i = 0; i < m_numBg; i++)
    {
        m_bgRange.push_back(pDedxSimData->getRange(i));
    }
    for (G4int i = 0; i < m_numDedxHists; i++)
    {
        m_dedx_hists[i] = pDedxSimData->getHist(i);
    }

    //get CalibCurSvc
    IDedxCurSvc* tmp_dedxCurSvc;
    sc = Gaudi::svcLocator()->service("DedxCurSvc", tmp_dedxCurSvc, true);
    if (!sc.isSuccess())
    {
        std::cout << "Could not get DedxCurSvc"
                  << std::endl;
    }
    m_pDedxCurSvc = dynamic_cast<DedxCurSvc*>(tmp_dedxCurSvc);

    if(m_G4Svc->MdcRootFlag())
    {
      m_tupleMdc = m_G4Svc->GetTupleMdc();
      sc = m_tupleMdc->addItem("betaGamma",m_betaGamma);
      sc = m_tupleMdc->addItem("fitval",m_fitval);
      sc = m_tupleMdc->addItem("dedx",m_dedx);
      sc = m_tupleMdc->addItem("de",m_de);
      //sc = m_tupleMdc->addItem("length",m_length);
      //sc = m_tupleMdc->addItem("random",m_random);
      sc = m_tupleMdc->addItem("charge", m_charge);
      sc = m_tupleMdc->addItem("costheta", m_costheta);
    }
}

BesMdcSD::~BesMdcSD

(

)

Definition at line 109 of file BesMdcSD.cc.

References m_dedx_hists.

                   {
  delete []m_dedx_hists;
}

---

Member Function Documentation

virtual void BesSensitiveDetector::BeginOfTrack

(

const G4Track *

)

[inline, virtual, inherited]

Definition at line 72 of file BesSensitiveDetector.hh.

{;}

void BesMdcSD::BeginOfTruthEvent

(

const G4Event *

)

[virtual]

Reimplemented from BesSensitiveDetector.

Definition at line 131 of file BesMdcSD.cc.

References truthCollection.

{  
  truthCollection = new BesMdcHitsCollection
                          (SensitiveDetectorName,collectionName[1]); 
  //  G4cout<<" begin event "<<evt->GetEventID()<<G4endl;
}

void BesMdcSD::dedxFuncInti

(

void

)

Definition at line 436 of file BesMdcSD.cc.

References dEdE_mylanfunc.

Referenced by BesMdcSD().

{
    dEdE_mylanfunc = new TF1("dEdE_mylanfunc",
                             "[3]*TMath::Exp([2]*((x[0]-[0])/[1]+TMath::Exp(-1*((x[0]-[0])/[1]))))",
                             0,
                             7500);
    //dEdE_mylanfunc->SetParameters(2009.35,559.776,-1.0932,6327.38);
    dEdE_mylanfunc->SetParNames("MPV","Sigma","constant1","constant2");
}

G4double BesMdcSD::dedxSample	(	G4double	betagamma,
		G4double	length,
		G4double	theta
	)			[private]

Definition at line 446 of file BesMdcSD.cc.

References cos(), GetAngleIndex(), GetBetagammaIndex(), GetChargeIndex(), GetValDedxCurve(), m_betaGamma, m_bgRange, m_charge, m_costheta, m_de, m_dedx, m_dedx_hists, m_fitval, m_G4Svc, m_numBg, m_numTheta, m_tupleMdc, G4Svc::MdcRootFlag(), and x.

Referenced by ProcessHits().

{
    G4double x = betagamma;
    G4double fitval = GetValDedxCurve(x, charge);
    if(fitval <= 0)return 0;
   
    G4double random1, random2, dedx1, dedx2, de;
    G4double standard1, standard2, beta_temp1, beta_temp2;
    G4double dedx = -1;

    G4int range_idx, bg_idx, angle_idx, charge_idx, hist_idx;
    range_idx = GetBetagammaIndex(betagamma);
    angle_idx = GetAngleIndex(theta);
    charge_idx = GetChargeIndex(charge);
    
    if (range_idx == -1)
    {
        while (dedx <= 0)
        {
            bg_idx = 0;
            hist_idx = bg_idx * 2 * m_numTheta + angle_idx * 2 + charge_idx;
            random1 = m_dedx_hists[hist_idx].GetRandom();
            beta_temp1 = (m_bgRange[0] + m_bgRange[1]) / 2;
            standard1 = GetValDedxCurve(beta_temp1, charge); 
            dedx = random1 + fitval - standard1;
        }
    }
    else if (range_idx == m_numBg - 1)
    {
        while (dedx <= 0)
        {
            bg_idx = (G4int)(range_idx / 2);
            hist_idx = bg_idx * 2 * m_numTheta + angle_idx * 2 + charge_idx;
            random1 = m_dedx_hists[hist_idx].GetRandom();
            beta_temp1 = (m_bgRange[m_numBg - 2] +
                          m_bgRange[m_numBg - 1]) / 2;
            standard1 = GetValDedxCurve(beta_temp1, charge);
            dedx = random1 + fitval - standard1;
        }     
    }
    else
    {
        //case 1: Given betagamma fall in one histograph range
        if (range_idx % 2 == 0)
        {
            while (dedx <= 0)
            {
                bg_idx = (G4int)(range_idx / 2);
                hist_idx = bg_idx * 2 * m_numTheta + angle_idx * 2 + charge_idx;
                random1 = m_dedx_hists[hist_idx].GetRandom();
                beta_temp1 = (m_bgRange[range_idx] + 
                              m_bgRange[range_idx + 1]) / 2;
                standard1 = GetValDedxCurve(beta_temp1, charge); 
                dedx1 = random1 + fitval - standard1;
                dedx = dedx1;
            }
        }
        //case 2: Given betagamma fall in interval between
        //        two histographs.
        else
        {
            while (dedx <= 0)
            {
                //standard1
                beta_temp1 = (m_bgRange[range_idx - 1] + 
                              m_bgRange[range_idx]) / 2;
                standard1 = GetValDedxCurve(beta_temp1, charge);

                //stardard2
                beta_temp2 = (m_bgRange[range_idx + 1] + 
                              m_bgRange[range_idx + 2]) / 2;
                standard2 = GetValDedxCurve(beta_temp2, charge);
               
                //random1 
                bg_idx = (G4int)(range_idx / 2);
                hist_idx = bg_idx * 2 * m_numTheta + angle_idx * 2 + charge_idx;
                random1 = m_dedx_hists[hist_idx].GetRandom();

                //random2
                bg_idx++;
                hist_idx = bg_idx * 2 * m_numTheta + angle_idx * 2 + charge_idx;
                random2 = m_dedx_hists[hist_idx].GetRandom();

                //combine dedx1 and dedx2
                dedx1 = random1 + fitval - standard1;
                dedx2 = random2 + fitval - standard2;
                dedx = (dedx2 * (x - m_bgRange[range_idx]) + 
                        dedx1 * (m_bgRange[range_idx + 1] - x)) / 
                       (m_bgRange[range_idx + 1] - m_bgRange[range_idx]);
            }
        }
    }
   
   
    de= dedx;// * y/10./1.5;// stepLength unit is mm in Geant4, cm in Rec.& Cal. software 
                         // dedx counts *1.5 in dedx Cal.
   
    if(m_G4Svc->MdcRootFlag())
    {
       m_betaGamma= x;
       m_fitval= fitval;
       //m_trackId = trackId; 
       //m_layer = layerId; 
       //m_wire  = cellId; 
       //m_random= random;
       m_dedx= dedx;
       m_de= de;
       //m_length=y;
       m_charge = charge;
       m_costheta = cos(theta);
       m_tupleMdc->write(); 
    }
    return de;
}

G4double BesMdcSD::Distance	(	G4int	,
		G4int	,
		G4ThreeVector	,
		G4ThreeVector	,
		G4ThreeVector &	,
		G4double &
	)

Definition at line 375 of file BesMdcSD.cc.

References MdcGeoWire::Backward(), MdcGeoWire::Forward(), MdcGeomSvc::Layer(), MdcGeoLayer::Length(), mdcGeomSvc, and MdcGeomSvc::Wire().

Referenced by ProcessHits().

{
  //For two lines r=r1+t1.v1 & r=r2+t2.v2 
  //the closest approach is d=|(r2-r1).(v1 x v2)|/|v1 x v2|
  //the point where closest approach are
  //t1=(v1 x v2).[(r2-r1) x v2]/[(v1 x v2).(v1 x v2)]
  //t2=(v1 x v2).[(r2-r1) x v1]/[(v1 x v2).(v1 x v2)]
  //if v1 x v2=0 means two lines are parallel
  //d=|(r2-r1) x v1|/|v1| 

  G4double t1,distance,dInOut,dHitIn,dHitOut;
  //Get wirepoint @ endplate
  G4ThreeVector east=mdcGeomSvc->Wire(layerId,cellId)->Backward();
  G4ThreeVector west=mdcGeomSvc->Wire(layerId,cellId)->Forward();
  G4ThreeVector wireLine=east-west;
  G4ThreeVector hitLine=pointOut-pointIn;

  G4ThreeVector hitXwire=hitLine.cross(wireLine);
  G4ThreeVector wire2hit=east-pointOut;
  //Hitposition is the position on hit line where closest approach 
  //of two lines, but it may out the area from pointIn to pointOut
  if(hitXwire.mag()==0){
    distance=wireLine.cross(wire2hit).mag()/wireLine.mag();
    hitPosition=pointIn;
  }else{
    t1=hitXwire.dot(wire2hit.cross(wireLine))/hitXwire.mag2();
    hitPosition=pointOut+t1*hitLine;
    
    dInOut=(pointOut-pointIn).mag();
    dHitIn=(hitPosition-pointIn).mag();
    dHitOut=(hitPosition-pointOut).mag();
    if(dHitIn<=dInOut && dHitOut<=dInOut){ //Between point in & out
      distance=fabs(wire2hit.dot(hitXwire)/hitXwire.mag());
    }else if(dHitOut>dHitIn){ // out pointIn
      distance=wireLine.cross(pointIn-east).mag()/wireLine.mag();
      hitPosition=pointIn;
    }else{ // out pointOut
      distance=wireLine.cross(pointOut-east).mag()/wireLine.mag(); 
      hitPosition=pointOut;
    }
  }

  //Calculate signal transferT on wire
  G4double halfLayerLength=mdcGeomSvc->Layer(layerId)->Length()/2.;
  G4double halfWireLength=wireLine.mag()/2.;
  G4double transferZ=0;
  if(layerId%2==0){
    transferZ=halfLayerLength+hitPosition.z(); //West readout
  }else{
    transferZ=halfLayerLength-hitPosition.z(); //East readout
  }
  if(layerId<8){
    transferT=transferZ*halfWireLength/halfLayerLength/220;
  }else{
    transferT=transferZ*halfWireLength/halfLayerLength/240;
  }

  return distance;

}

void BesMdcSD::EndOfEvent

(

G4HCofThisEvent *

)

Definition at line 362 of file BesMdcSD.cc.

References hitsCollection.

{
  if (verboseLevel>0) { 
  hitsCollection->PrintAllHits();
  /*
     G4int NbHits = hitsCollection->entries();
     G4cout << "\n-------->Hits Collection: in this event they are " << NbHits 
            << " hits in the MDC chambers: " << G4endl;
     for (G4int i=0;i<NbHits;i++) (*hitsCollection)[i]->Print();
  */
  } 
}

virtual void BesSensitiveDetector::EndOfTrack

(

const G4Track *

)

[inline, virtual, inherited]

Definition at line 74 of file BesSensitiveDetector.hh.

{;}

void BesMdcSD::EndOfTruthEvent

(

const G4Event *

)

[virtual]

Reimplemented from BesSensitiveDetector.

Definition at line 138 of file BesMdcSD.cc.

References truthCollection.

{  static G4int HLID=-1;
  if(HLID<0)
  {
    HLID = G4SDManager::GetSDMpointer()->
                GetCollectionID(collectionName[1]); 
  }
  G4HCofThisEvent* HCE = evt->GetHCofThisEvent();
  HCE->AddHitsCollection(HLID,truthCollection);
}

G4int BesMdcSD::GetAngleIndex

(

G4double

)

[private]

Definition at line 628 of file BesMdcSD.cc.

References cos(), m_numTheta, and m_version.

Referenced by dedxSample().

{   
    if (m_version == 0) {
        if (fabs(cos(theta)) >= 0.93) return 9;
        return (G4int)(fabs(cos(theta)) * 10 / 0.93);
    }
    else {
        G4double cos_min = -0.93;
        G4double cos_max = 0.93;
        G4int nbin = m_numTheta;
        G4double cos_step = (cos_max - cos_min) / nbin;

        if (cos(theta) < cos_min) return 0;
        else if (cos_min < cos(theta) && cos(theta) < cos_max) {
            return (G4int)((cos(theta) - cos_min) / cos_step);
        }
        else return nbin - 1;
    }
}

G4int BesMdcSD::GetBetagammaIndex

(

G4double

bg

)

[private]

Definition at line 607 of file BesMdcSD.cc.

References genRecEmupikp::i, m_bgRange, and m_numBg.

Referenced by dedxSample().

{
    if (bg < m_bgRange[0]) return -1;
    for (G4int i = 0; i < m_numBg - 1; i++)
    {
        if (bg > m_bgRange[i] && bg < m_bgRange[i + 1])
        {
            return i;
        }
    }
    if (bg > m_bgRange[m_numBg - 1])
        return (m_numBg - 1);
    return -1;
}

G4int BesMdcSD::GetChargeIndex

(

G4int

)

[private]

Definition at line 653 of file BesMdcSD.cc.

Referenced by dedxSample().

{
    if (charge > 0)  return 0;
    if (charge == 0) return -1;  // warning: -1 is illegal, for charged tracks are expected.
    if (charge < 0)  return 1;
    return -1;
}

void BesSensitiveDetector::GetCurrentTrackIndex	(	G4int &	trackIndex,
		G4int &	g4TrackId
	)			const [protected, inherited]

Definition at line 30 of file BesSensitiveDetector.cc.

References BesSensitiveManager::GetCurrentTrackIndex(), BesSensitiveManager::GetSensitiveManager(), BesSensitiveManager::GetTrackList(), genRecEmupikp::i, and delete_small_size::size.

Referenced by BesTofSD::ProcessHits(), BesMucSD::ProcessHits(), ProcessHits(), and BesEmcSD::ProcessHits().

{
  BesSensitiveManager* manager = BesSensitiveManager::GetSensitiveManager();
  if(manager)
  {
    trackIndex = manager->GetCurrentTrackIndex();
    std::vector<BesTruthTrack*>* trackList = manager->GetTrackList();
    if(trackList)
    {
      G4int size = trackList->size();
      if(size>0)
      {
        for(G4int i=0;i<size;i++)
        {
          if( (*trackList)[i]->GetIndex() == trackIndex )
          {
            g4TrackId = (*trackList)[i]->GetG4TrackId();
            break;
          }
        }
      } 
    }
  }
}

G4double BesMdcSD::GetValDedxCurve	(	G4double	bg,
		G4double	charge
	)			[private]

Definition at line 566 of file BesMdcSD.cc.

References EvtCyclic3::A, EvtCyclic3::B, EvtCyclic3::C, exp(), IDedxCurSvc::getCurve(), IDedxCurSvc::getCurveSize(), genRecEmupikp::i, m_pDedxCurSvc, and delete_small_size::size.

Referenced by dedxSample().

{
    G4int dedxflag = -1;
    G4int size = -1;
    G4double A = 0.;
    G4double B = 0.;
    G4double C = 0.;
    std::vector<G4double> par;
    G4double val;
    G4int index = -1;

    par.clear();
    dedxflag = m_pDedxCurSvc->getCurve(0);
    size = m_pDedxCurSvc->getCurveSize();
    for (G4int i = 1; i < size; i++) {
        par.push_back(m_pDedxCurSvc->getCurve(i));
    }

    if (x < 4.5)
        A = 1;
    else if(x < 10)
        B = 1;
    else
        C = 1;
    
    G4double partA = par[0] * pow(sqrt(x * x + 1), par[2]) / pow(x, par[2]) *
                  (par[1] - par[5] * log(pow(1 / x, par[3]))) -
                   par[4] + exp(par[6] + par[7] * x);
    G4double partB = par[8] * pow(x, 3) + par[9] * pow(x, 2) + par[10] * x + par[11];
    G4double partC = - par[12] * log(par[15] + pow(1 / x, par[13])) + par[14];
    
    val = 550 * (A * partA + B * partB + C * partC);
    return val;
    
}

void BesMdcSD::Initialize

(

G4HCofThisEvent *

)

Definition at line 113 of file BesMdcSD.cc.

References hitPointer, hitsCollection, genRecEmupikp::i, ganga-rec::j, and truthPointer.

{
  hitsCollection = new BesMdcHitsCollection
                          (SensitiveDetectorName,collectionName[0]); 
  static G4int HCID = -1;
  if(HCID<0)
  { HCID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]); }
  HCE->AddHitsCollection( HCID, hitsCollection ); 
  G4int i,j;
  for(i=0; i<43;i++){
    for(j=0;j<288;j++){
      hitPointer[i][j]=-1;
      truthPointer[i][j]=-1;
    }
  }
}

G4bool BesMdcSD::ProcessHits	(	G4Step *	,
		G4TouchableHistory *
	)

Definition at line 149 of file BesMdcSD.cc.

References BesMdcCalTransfer::D2T(), dedxSample(), Distance(), BesSensitiveDetector::GetCurrentTrackIndex(), ReadBoostRoot::GetMdc(), G4Svc::GetMdcDedxFlag(), hitPointer, hitsCollection, MdcGeomSvc::Layer(), MdcGeoLayer::Length(), m_G4Svc, mdcCalPointer, mdcGeomSvc, mdcGeoPointer, BesMdcWire::Phi(), pi, BesMdcHit::SetCellNo(), BesMdcHit::SetDriftD(), BesMdcHit::SetDriftT(), BesMdcHit::SetEdep(), BesMdcHit::SetEnterAngle(), BesMdcHit::SetGlobalT(), BesMdcCalTransfer::SetHitPointer(), BesMdcHit::SetLayerNo(), BesMdcHit::SetPos(), BesMdcHit::SetPosFlag(), BesMdcHit::SetTheta(), BesMdcHit::SetTrackID(), BesMdcGeoParameter::SignalWire(), truthCollection, and truthPointer.

{
  G4Track* gTrack = aStep->GetTrack() ;

  G4double globalT=gTrack->GetGlobalTime();//Time since the event in which the track belongs is created
  if(isnan(globalT)){
    G4cout<<"MdcSD:error, globalT is nan "<<G4endl;
    return false;
  }
  if(globalT > 2000)return false; //MDC T window is 2 microsecond 

  //skip neutral tracks
  G4int charge = gTrack->GetDefinition()->GetPDGCharge();
  if (charge == 0) return false;

  //skip no energy deposit tracks
  G4double stepLength=aStep->GetStepLength();
  if(stepLength==0){
    //    G4cout<<"step length equal 0!!"<<G4endl;
    return false;
  }

  G4double edep = aStep->GetTotalEnergyDeposit() / stepLength;
  if(edep==0.) return false;

  // get position of the track at the beginning and at the end of step
  G4StepPoint* prePoint  = aStep->GetPreStepPoint() ;
  G4StepPoint* postPoint = aStep->GetPostStepPoint() ;

  //get position coordinate
  G4ThreeVector pointIn  = prePoint->GetPosition();
  G4ThreeVector pointOut = postPoint->GetPosition();

  // get physical volumes
  const G4VTouchable *touchable = prePoint->GetTouchable();
  G4VPhysicalVolume *volume = touchable->GetVolume(0);

  G4double driftD = 0;
  G4double driftT = 0;
  G4double edepTemp = 0;
  G4double lengthTemp = 0;
  G4int cellId=0;
  G4int layerId = touchable->GetVolume(1)->GetCopyNo();
  if(volume->IsReplicated()){
    cellId  = touchable->GetReplicaNumber();
  }else{
    cellId=touchable->GetVolume(0)->GetCopyNo();
  }
  if(layerId==18&&(cellId==27||cellId==42))return false; // no sense wire 

  if(ReadBoostRoot::GetMdc() == 2) { //gdml
    //layerId  0-35 -> CopyNo 0-35 in gdml 
    //layerId 36-42 -> CopyNo (36,37),(38,39),...(48,49)    
    if(layerId >= 36) layerId = 36 + (layerId-36)/2;
  }

  G4double halfLayerLength=mdcGeomSvc->Layer(layerId)->Length()/2.;
  if(((fabs(pointIn.z())-halfLayerLength)>-7.)
     &&((fabs(pointOut.z())-halfLayerLength)>-7.))return false;//Out sensitive area

  G4int trackID= gTrack->GetTrackID(); //G4 track ID of current track.
  G4int truthID, g4TrackID;
  GetCurrentTrackIndex(truthID, g4TrackID); //ID of current primary track.
  
  G4double theta=gTrack->GetMomentumDirection().theta();

  G4ThreeVector hitPosition=0;
  G4double transferT=0;
  driftD  = Distance(layerId,cellId,pointIn,pointOut,hitPosition,transferT);

  G4double posPhi, wirePhi;
  posPhi=hitPosition.phi();//from -pi to pi
  if(posPhi<0)posPhi += 2*pi;
  wirePhi=mdcGeoPointer->SignalWire(layerId,cellId).Phi(hitPosition.z());//from 0 to 2pi

  G4int posFlag;
  if(posPhi<=wirePhi){
    posFlag    = 0;
  }else{
    posFlag    = 1;
  }
  // if x axis is between pos and wire, phi will has a jump of one of them.
  if(posPhi < 1. && wirePhi > 5.)posFlag = 1;
  if(posPhi > 5. && wirePhi < 1.)posFlag = 0;

  G4ThreeVector hitLine=pointOut-pointIn;
  G4double enterAngle=hitLine.phi()-wirePhi;
  while(enterAngle<-pi/2.)enterAngle+=pi;
  while(enterAngle>pi/2.)enterAngle-=pi;

  if(m_G4Svc->GetMdcDedxFlag()==1){
     G4double betaGamma=aStep->GetPreStepPoint()->GetBeta() * aStep->GetPreStepPoint()->GetGamma();
     if(betaGamma<0.01)return false;//too low momentum
     //if (betaGamma < 10.0) betaGamma = 10.0;

     G4double eCount=dedxSample(betaGamma, charge, theta);
     edep=eCount;
  }

  BesMdcHit* newHit = new BesMdcHit();
  newHit->SetTrackID(truthID);
  //newHit->SetTrkID(trackID); 
  newHit->SetLayerNo(layerId);
  newHit->SetCellNo(cellId);
  newHit->SetEdep(edep);
  newHit->SetPos(hitPosition);
  newHit->SetDriftD(driftD);
  newHit->SetTheta(theta);
  newHit->SetPosFlag(posFlag);
  newHit->SetEnterAngle(enterAngle);

  //Transfer hit pointer to BesMdcCalTransfer
  mdcCalPointer->SetHitPointer(newHit);

  driftT=mdcCalPointer->D2T(driftD);
  globalT+=transferT;
  driftT+=globalT;
  
  newHit->SetDriftT (driftT);
  newHit->SetGlobalT(globalT);

  if (hitPointer[layerId][cellId] == -1) {
      hitsCollection->insert(newHit);
      G4int NbHits = hitsCollection->entries();
      hitPointer[layerId][cellId]=NbHits-1;
  } 
  else
  {
      G4int pointer=hitPointer[layerId][cellId]; 
      if (g4TrackID == trackID) {
          G4double preDriftT=(*hitsCollection)[pointer]->GetDriftT();
      }
      
      G4double preDriftT = (*hitsCollection)[pointer]->GetDriftT();
      if (driftT < preDriftT) {   
          (*hitsCollection)[pointer]->SetTrackID(truthID);
          //(*hitsCollection)[pointer]->SetTrkID(trackID);
          (*hitsCollection)[pointer]->SetDriftD(driftD);
          (*hitsCollection)[pointer]->SetDriftT(driftT);
          (*hitsCollection)[pointer]->SetPos(hitPosition);
          (*hitsCollection)[pointer]->SetGlobalT(globalT);
          (*hitsCollection)[pointer]->SetTheta(theta);
          (*hitsCollection)[pointer]->SetPosFlag(posFlag);
          (*hitsCollection)[pointer]->SetEnterAngle(enterAngle);
      }    
   
      delete newHit;
  }

  //for mc truth 
  if(truthCollection){
    if(g4TrackID==trackID){ //This track is the primary track & will appear in MC truth
      G4int pointer=truthPointer[layerId][cellId];  
      if(pointer==-1){
        BesMdcHit* truthHit = new BesMdcHit();
        truthHit->SetTrackID  (truthID);
        truthHit->SetLayerNo(layerId);
        truthHit->SetCellNo(cellId);
        truthHit->SetEdep     (edep);
        truthHit->SetPos      (hitPosition);
        truthHit->SetDriftD (driftD);
        truthHit->SetDriftT (driftT);
        truthHit->SetGlobalT(globalT);
        truthHit->SetTheta(theta);
        truthHit->SetPosFlag(posFlag);
        truthHit->SetEnterAngle(enterAngle);

        truthCollection->insert(truthHit);
        G4int NbHits = truthCollection->entries();
        truthPointer[layerId][cellId]=NbHits-1;
      } 
      else {
        if(truthID == (*truthCollection)[pointer]->GetTrackID()){
            G4double preDriftT=(*truthCollection)[pointer]->GetDriftT();
            if(driftT<preDriftT){
                (*truthCollection)[pointer]->SetDriftD(driftD);
                (*truthCollection)[pointer]->SetDriftT(driftT);
                (*truthCollection)[pointer]->SetPos(hitPosition);
                (*truthCollection)[pointer]->SetGlobalT(globalT);
                (*truthCollection)[pointer]->SetTheta(theta);
                (*truthCollection)[pointer]->SetPosFlag(posFlag);
                (*truthCollection)[pointer]->SetEnterAngle(enterAngle);
                }
            edepTemp=(*truthCollection)[pointer]->GetEdep();
            (*truthCollection)[pointer]->SetEdep(edepTemp+edep);
            } else {
                BesMdcHit* truthHit = new BesMdcHit();
                truthHit->SetTrackID  (truthID);
                truthHit->SetLayerNo(layerId);
                truthHit->SetCellNo(cellId);
                truthHit->SetEdep(edep);
                truthHit->SetPos(hitPosition);
                truthHit->SetDriftD (driftD);
                truthHit->SetDriftT (driftT);
                truthHit->SetGlobalT(globalT);
                truthHit->SetTheta(theta);
                truthHit->SetPosFlag(posFlag);
                truthHit->SetEnterAngle(enterAngle);

                truthCollection->insert(truthHit);
                G4int NbHits = truthCollection->entries();
                truthPointer[layerId][cellId]=NbHits-1;
            }
      }
    }
  }

  //newHit->Print();
//  newHit->Draw();

  return true;
}

---

Member Data Documentation

TF1* BesMdcSD::dEdE_mylanfunc [private]

Definition at line 60 of file BesMdcSD.hh.

Referenced by dedxFuncInti().

G4int BesMdcSD::hitPointer[43][288] [private]

Definition at line 53 of file BesMdcSD.hh.

Referenced by Initialize(), and ProcessHits().

BesMdcHitsCollection* BesMdcSD::hitsCollection [private]

Definition at line 54 of file BesMdcSD.hh.

Referenced by EndOfEvent(), Initialize(), and ProcessHits().

NTuple::Item<double> BesMdcSD::m_betaGamma [private]

Definition at line 79 of file BesMdcSD.hh.

Referenced by BesMdcSD(), and dedxSample().

std::vector<G4double> BesMdcSD::m_bgRange [private]

Definition at line 70 of file BesMdcSD.hh.

Referenced by BesMdcSD(), dedxSample(), and GetBetagammaIndex().

CalibDataSvc* BesMdcSD::m_calibDataSvc [private]

dedx sim ---------------------------

Definition at line 63 of file BesMdcSD.hh.

Referenced by BesMdcSD().

NTuple::Item<double> BesMdcSD::m_charge [private]

Definition at line 85 of file BesMdcSD.hh.

Referenced by BesMdcSD(), and dedxSample().

NTuple::Item<double> BesMdcSD::m_costheta [private]

Definition at line 86 of file BesMdcSD.hh.

Referenced by BesMdcSD(), and dedxSample().

NTuple::Item<double> BesMdcSD::m_de [private]

Definition at line 83 of file BesMdcSD.hh.

Referenced by BesMdcSD(), and dedxSample().

NTuple::Item<double> BesMdcSD::m_dedx [private]

Definition at line 82 of file BesMdcSD.hh.

Referenced by BesMdcSD(), and dedxSample().

TH1F* BesMdcSD::m_dedx_hists [private]

Definition at line 65 of file BesMdcSD.hh.

Referenced by BesMdcSD(), dedxSample(), and ~BesMdcSD().

NTuple::Item<double> BesMdcSD::m_fitval [private]

Definition at line 80 of file BesMdcSD.hh.

Referenced by BesMdcSD(), and dedxSample().

G4Svc* BesMdcSD::m_G4Svc [private]

Definition at line 59 of file BesMdcSD.hh.

Referenced by BesMdcSD(), dedxSample(), and ProcessHits().

G4int BesMdcSD::m_numBg [private]

Definition at line 68 of file BesMdcSD.hh.

Referenced by BesMdcSD(), dedxSample(), and GetBetagammaIndex().

G4int BesMdcSD::m_numDedxHists [private]

Definition at line 67 of file BesMdcSD.hh.

Referenced by BesMdcSD().

G4int BesMdcSD::m_numTheta [private]

Definition at line 69 of file BesMdcSD.hh.

Referenced by BesMdcSD(), dedxSample(), and GetAngleIndex().

IDedxCurSvc* BesMdcSD::m_pDedxCurSvc [private]

Definition at line 64 of file BesMdcSD.hh.

Referenced by BesMdcSD(), and GetValDedxCurve().

NTuple::Item<double> BesMdcSD::m_random [private]

Definition at line 81 of file BesMdcSD.hh.

NTuple::Tuple* BesMdcSD::m_tupleMdc [private]

Definition at line 78 of file BesMdcSD.hh.

Referenced by BesMdcSD(), and dedxSample().

G4int BesMdcSD::m_version [private]

Definition at line 66 of file BesMdcSD.hh.

Referenced by BesMdcSD(), and GetAngleIndex().

BesMdcCalTransfer* BesMdcSD::mdcCalPointer [private]

Definition at line 57 of file BesMdcSD.hh.

Referenced by BesMdcSD(), and ProcessHits().

MdcGeomSvc* BesMdcSD::mdcGeomSvc [private]

Definition at line 58 of file BesMdcSD.hh.

Referenced by BesMdcSD(), Distance(), and ProcessHits().

BesMdcGeoParameter* BesMdcSD::mdcGeoPointer [private]

Definition at line 56 of file BesMdcSD.hh.

Referenced by BesMdcSD(), and ProcessHits().

BesMdcHitsCollection* BesMdcSD::truthCollection [private]

Definition at line 55 of file BesMdcSD.hh.

Referenced by BeginOfTruthEvent(), EndOfTruthEvent(), and ProcessHits().

G4int BesMdcSD::truthPointer[43][288] [private]

Definition at line 53 of file BesMdcSD.hh.

Referenced by Initialize(), and ProcessHits().

---