BesMdcSD Class Reference

#include <BesMdcSD.hh>

Inheritance diagram for BesMdcSD:

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)
	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)
	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)
Public Member Functions inherited from BesSensitiveDetector
	BesSensitiveDetector (const G4String name)
virtual	~BesSensitiveDetector ()
virtual void	BeginOfTrack (const G4Track *)
virtual void	EndOfTrack (const G4Track *)
	BesSensitiveDetector (const G4String name)
virtual	~BesSensitiveDetector ()
virtual void	BeginOfTrack (const G4Track *)
virtual void	EndOfTrack (const G4Track *)
	BesSensitiveDetector (const G4String name)
virtual	~BesSensitiveDetector ()
virtual void	BeginOfTrack (const G4Track *)
virtual void	EndOfTrack (const G4Track *)

Additional Inherited Members
Protected Member Functions inherited from BesSensitiveDetector
void	GetCurrentTrackIndex (G4int &trackIndex, G4int &g4TrackId) const
void	GetCurrentTrackIndex (G4int &trackIndex, G4int &g4TrackId) const
void	GetCurrentTrackIndex (G4int &trackIndex, G4int &g4TrackId) const

Detailed Description

Definition at line 37 of file InstallArea/x86_64-el9-gcc13-dbg/include/MdcSim/BesMdcSD.hh.

Constructor & Destructor Documentation

BesMdcSD() [1/3]

BesMdcSD::BesMdcSD

(

G4String

name

)

Definition at line 25 of file BesMdcSD.cc.

                                  : BesSensitiveDetector( name ) {
  collectionName.insert( "BesMdcHitsCollection" );
  collectionName.insert( "BesMdcTruthCollection" );
 
  mdcGeoPointer = BesMdcGeoParameter::GetGeo();
  mdcCalPointer = new BesMdcCalTransfer;
 
  StatusCode sc = Gaudi::svcLocator()->service( "MdcGeomSvc", mdcGeomSvc );
  if ( !sc.isSuccess() ) std::cout << "BesMdcSD::Could not open MdcGeomSvc" << std::endl;
 
  sc = Gaudi::svcLocator()->service( "G4Svc", m_G4Svc );
  if ( !sc.isSuccess() ) G4cout << " MdcSD::Error,could not open G4Svc" << G4endl;
 
  if ( m_G4Svc->GetMdcDedxFlag() == 1 )
  {
    G4cout << " MdcSD: Use sampled dedx instead of Geant4 value" << G4endl;
    dedxFuncInti();
  }
 
  // get DedxSimSvc
  sc = Gaudi::svcLocator()->service( "DedxSimSvc", m_pDedxSimSvc, true );
  if ( !sc.isSuccess() ) { std::cout << "Could not get DedxSimSvc" << std::endl; }
 
  // get CalibCurSvc
  IDedxCurSvc* tmp_dedxCurSvc;
  sc = Gaudi::svcLocator()->service( "DedxCurSvc", m_pDedxCurSvc, true );
  if ( !sc.isSuccess() ) { std::cout << "Could not get DedxCurSvc" << std::endl; }
 
  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() [1/3]

BesMdcSD::~BesMdcSD

(

)

Definition at line 67 of file BesMdcSD.cc.

{}

BesMdcSD() [2/3]

BesMdcSD::BesMdcSD

(

G4String

)

~BesMdcSD() [2/3]

BesMdcSD::~BesMdcSD

(

)

BesMdcSD() [3/3]

BesMdcSD::BesMdcSD

(

G4String

)

~BesMdcSD() [3/3]

BesMdcSD::~BesMdcSD

(

)

Member Function Documentation

BeginOfTruthEvent() [1/3]

void BesMdcSD::BeginOfTruthEvent ( const G4Event * evt )

virtual

Reimplemented from BesSensitiveDetector.

Definition at line 87 of file BesMdcSD.cc.

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

BeginOfTruthEvent() [2/3]

void BesMdcSD::BeginOfTruthEvent ( const G4Event * )

virtual

Reimplemented from BesSensitiveDetector.

BeginOfTruthEvent() [3/3]

void BesMdcSD::BeginOfTruthEvent ( const G4Event * )

virtual

Reimplemented from BesSensitiveDetector.

dedxFuncInti() [1/3]

void BesMdcSD::dedxFuncInti

(

void

)

Definition at line 405 of file BesMdcSD.cc.

                                  {
  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" );
}

Referenced by BesMdcSD().

dedxFuncInti() [2/3]

void BesMdcSD::dedxFuncInti

(

void

)

dedxFuncInti() [3/3]

void BesMdcSD::dedxFuncInti

(

void

)

Distance() [1/3]

G4double BesMdcSD::Distance	(	G4int	layerId,
		G4int	cellId,
		G4ThreeVector	pointIn,
		G4ThreeVector	pointOut,
		G4ThreeVector &	hitPosition,
		G4double &	transferT )

Definition at line 336 of file BesMdcSD.cc.

                                                   {
  // 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;
}

Referenced by ProcessHits().

Distance() [2/3]

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

Distance() [3/3]

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

EndOfEvent() [1/3]

void BesMdcSD::EndOfEvent

(

G4HCofThisEvent *

)

Definition at line 323 of file BesMdcSD.cc.

                                            {
  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();
    */
  }
}

EndOfEvent() [2/3]

void BesMdcSD::EndOfEvent

(

G4HCofThisEvent *

)

EndOfEvent() [3/3]

void BesMdcSD::EndOfEvent

(

G4HCofThisEvent *

)

EndOfTruthEvent() [1/3]

void BesMdcSD::EndOfTruthEvent ( const G4Event * evt )

virtual

Reimplemented from BesSensitiveDetector.

Definition at line 92 of file BesMdcSD.cc.

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

EndOfTruthEvent() [2/3]

void BesMdcSD::EndOfTruthEvent ( const G4Event * )

virtual

Reimplemented from BesSensitiveDetector.

EndOfTruthEvent() [3/3]

void BesMdcSD::EndOfTruthEvent ( const G4Event * )

virtual

Reimplemented from BesSensitiveDetector.

Initialize() [1/3]

void BesMdcSD::Initialize

(

G4HCofThisEvent *

HCE

)

Definition at line 69 of file BesMdcSD.cc.

                                                {
  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;
    }
  }
}

Initialize() [2/3]

void BesMdcSD::Initialize

(

G4HCofThisEvent *

)

Initialize() [3/3]

void BesMdcSD::Initialize

(

G4HCofThisEvent *

)

ProcessHits() [1/3]

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

Definition at line 100 of file BesMdcSD.cc.

                                                                 {
  G4Track* gTrack = aStep->GetTrack();
 
  G4double globalT =
      gTrack->GetGlobalTime(); // Time since the event in which the track belongs is created
  if ( std::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, 0, 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;
}

ProcessHits() [2/3]

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

ProcessHits() [3/3]

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

---

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

• InstallArea/x86_64-el9-gcc13-dbg/include/MdcSim/BesMdcSD.hh
• InstallArea/x86_64-el9-gcc13-opt/include/MdcSim/BesMdcSD.hh
• Simulation/BOOST/MdcSim/include/MdcSim/BesMdcSD.hh
• BesMdcSD.cc