BesTofSD.cc

Go to the documentation of this file.

//---------------------------------------------------------------------------//
//      BOOST --- BESIII Object_Oriented Simulation Tool                     //
//---------------------------------------------------------------------------//
//Description:
//Author: Dengzy
//Created: Mar, 2004
//Modified:
//Comment:
//---------------------------------------------------------------------------//
//$Id: BesTofSD.cc
 
#include "BesTofSD.hh"
#include "ReadBoostRoot.hh"
#include "G4HCofThisEvent.hh"
#include "G4Step.hh"
#include "G4ThreeVector.hh"
#include "G4SDManager.hh"
#include "G4ios.hh"
#include "G4UnitsTable.hh"
#include "G4Event.hh"
 
 
 
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/MsgStream.h"
 
 
#include "G4LossTableManager.hh"
#include "G4ElectronIonPair.hh"
 
 
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 
 
BesTofSD::BesTofSD(G4String name)
        :BesSensitiveDetector(name),m_besTofList(0)
{
    collectionName.insert("BesTofHitsCollection");
    collectionName.insert("BesTofHitsList");
    elIonPair = new G4ElectronIonPair();
 
}
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
BesTofSD::~BesTofSD()
{ }
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
void BesTofSD::Initialize(G4HCofThisEvent* )
{
    m_besTofCollection = new BesTofHitsCollection
    (SensitiveDetectorName,collectionName[0]);
 
    //elIonPair = G4LossTableManager::Instance()->ElectronIonPair();
}
 
//for MC Truth
void BesTofSD::BeginOfTruthEvent(const G4Event* )
{
    m_besTofList = new BesTofHitsCollection(SensitiveDetectorName, collectionName[1]);
    m_trackIndexes.clear();
    m_trackIndex = -99;
}
 
void BesTofSD::EndOfTruthEvent(const G4Event* evt)
{
    static G4int HLID=-1;
    if (HLID<0)
    {
        HLID = G4SDManager::GetSDMpointer()->
               GetCollectionID(collectionName[1]);
    }
 
    G4HCofThisEvent* HCE = evt->GetHCofThisEvent();
    HCE->AddHitsCollection(HLID,m_besTofList);
 
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
G4bool BesTofSD::ProcessHits(G4Step* aStep,G4TouchableHistory*) 
{
   
   G4double chg=aStep->GetTrack()->GetDefinition()->GetPDGCharge();
    G4double edep = aStep->GetTotalEnergyDeposit();
    G4double stepL=aStep->GetStepLength();
    G4double deltaT=aStep->GetDeltaTime();
    G4StepPoint* preStep = aStep->GetPreStepPoint();
    G4ThreeVector pDirection=preStep->GetMomentumDirection();
    G4String particleName = aStep->GetTrack()->GetDefinition()->GetParticleName();
    G4Material* scinMaterial = aStep->GetTrack()->GetMaterial();
    G4double charge = aStep->GetTrack()->GetDefinition()->GetPDGCharge();
    G4int pdgcode = aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
 
    if (chg==0&&edep==0&&stepL==0)
        return false;
    BesTofHit *newHit = new BesTofHit();
 
    
    G4int trackId = aStep->GetTrack()->GetTrackID();
 
    newHit->SetTrackIndex(-99);
    newHit->SetG4Index(aStep->GetTrack()->GetTrackID());
    newHit->SetEdep(edep);
    newHit->SetStepL(stepL);
    //newHit->SetPName(particleName);
    newHit->SetTrackL(aStep->GetTrack()->GetTrackLength());
    G4ThreeVector pos=preStep->GetPosition();
    newHit->SetPos(pos);
    G4double globalTime=preStep->GetGlobalTime();
    newHit->SetTime(globalTime);
    newHit->SetDeltaT(deltaT);
    newHit->SetPDirection(pDirection);
    newHit->SetMomentum(preStep->GetMomentum());
    //newHit->SetMaterial(scinMaterial);
    newHit->SetCharge(charge);
    newHit->SetPDGcode(pdgcode);
    //barrel segment number: 0-175
    //endcap west:           0-47      z<0
    //endcap east:           0-47      z>0
 
    G4TouchableHistory* theTouchable
      = (G4TouchableHistory*)(preStep->GetTouchable()); 
 
    //partId: barrel(1), east endcap(0), west endcap(2);
    //partId:mrpc: east_layer_2  (3), east_layer_1  (4), west_layer_1  (5), west_layer_2  (6), 
    //scinNb:  barrel(0-175),east endcap(0-47), west endcap(0-47)
    G4String name = theTouchable->GetVolume(0)->GetName();
 
    G4int partId, scinNb, number,help_mrpc,module_mrpc ;
 
   number = theTouchable->GetReplicaNumber(2);
   help_mrpc = theTouchable->GetReplicaNumber(3);
 
 
   //MRPC ENDCAPS construct by geant4 code
     if(name=="physical_sensitive_detector_west_1")
       {
     partId =5;
     scinNb = help_mrpc;
      }
     else if(name=="physical_sensitive_detector_west_2")
       {
     partId =6;
     scinNb = help_mrpc;
       }
     else if(name=="physical_sensitive_detector_east_1")
       {
     partId =4;
     scinNb = help_mrpc;
       }
     else if(name=="physical_sensitive_detector_east_2")
       {
     partId =3;
     scinNb = help_mrpc;
       }
     //MRPC ENDCAPS construct by GDML
     else if(name=="logical_sensitive_detector_west_1")
       {
         partId =5;
         scinNb = (help_mrpc)*(-0.5)+18.5;
       }
     else if(name=="logical_sensitive_detector_west_2")
       {
         partId =6;
         scinNb = (help_mrpc)*(-0.5)+18;
       }
     else if(name=="logical_sensitive_detector_east_1")
       {
         partId =4;
         scinNb = (help_mrpc)*(0.5)+0.5;
       }
 
     else if(name=="logical_sensitive_detector_east_2")
       {
         partId =3;
         scinNb =( help_mrpc)*(0.5)+1;
       }
     //old TOF
     else if (name=="physicalScinBr1")
       {
     partId = 1;
     scinNb = number;
       }
     else if (name=="physicalScinBr2")
       {
     partId = 1;
     scinNb = number+88;
       }
     else if (name=="physicalScinEcWest")
       {
     partId = 2;
     scinNb = number;
       }
     else if (name=="physicalScinEcEast")
       {
     partId = 0;
     scinNb = number;
       }
 
     //if construct by gdml files
     else if (name=="logicalScinBr1" ||name=="logicalScinBr2")
       {
     partId = 1;
     scinNb = (527-number)/3;
       }
     else if (name=="logicalScinEcEast")
       { //copy number: east:1-95(scinNb:47-0), west:1-95(scinNb:47-0)
     //east: partId=0  west: partId=2
     partId = 0;
     scinNb = (95-number)/2;
       }
     else if (name=="logicalScinEcWest")
       {
     partId  = 2;
     scinNb = (95-number)/2;
       }
     else
       return false;
 
 
 
 
 
 
 
 
    newHit->SetPartId(partId);
    newHit->SetScinNb(scinNb);
    newHit->SetModule_mrpc(scinNb);
    //newHit->SetIons(elIonPair->SampleNumberOfIonsAlongStep(aStep));
    newHit->SetIons(SampleNumberOfIonsAlongStep(aStep,elIonPair));
    
 
    
 
  //newHit->Draw();
    G4int trackIndex, g4TrackId;
    GetCurrentTrackIndex(trackIndex, g4TrackId);
    newHit->SetTrackIndex(trackIndex);
    //newHit->Print();
    if (edep>0)
      {
        m_besTofCollection->insert(newHit);
    
      }
 
    //for mc truth
    if (m_besTofList)
    {
        G4int trackIndex, g4TrackId;
        GetCurrentTrackIndex(trackIndex, g4TrackId);
        newHit->SetTrackIndex(trackIndex);
        if (m_trackIndex != trackIndex)
        {
            m_trackIndex = trackIndex;
            //G4int size = m_trackIndexes.size();
            //G4int flag=1;
            //if(size>0)
            //{
            //  for(G4int i=0;i<size;i++)
            //    if(m_trackIndexes[i] == trackIndex )
            //    {flag =0; break;}
            //}
            //if(flag && aStep->GetTrack()->GetTrackID()==g4TrackId )
            G4int flag=1;
            G4int pdg = abs(aStep->GetTrack()->GetDefinition()->GetPDGEncoding());
            if (pdg==12 || pdg==14 || pdg==16)
                flag=0;
            if (flag && aStep->GetTrack()->GetTrackID()==g4TrackId)
            {
                m_trackIndexes.push_back(trackIndex);
                BesTofHit* truHit = new BesTofHit();
                *truHit = *newHit;
                m_besTofList->insert(truHit);
            }
        }
    }
    if (edep<=0)  delete newHit;
 
 
 
return true;
 
}//Close Process Hits
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
void BesTofSD::EndOfEvent(G4HCofThisEvent* HCE)
{
    static G4int HCID=-1;
    if (HCID<0)
    {
        HCID = G4SDManager::GetSDMpointer()->
               GetCollectionID(collectionName[0]);
    }
    HCE->AddHitsCollection(HCID,m_besTofCollection);
 
}
 
 
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 
 
G4int BesTofSD::SampleNumberOfIonsAlongStep(const G4Step* step, G4ElectronIonPair* eipair)
{
 
  G4double FanoFactor = 0.2;  //Is this value appropiate?
  G4double meanion = eipair->MeanNumberOfIonsAlongStep(step);
  G4double sig = std::sqrt(FanoFactor*meanion);  //Be carefull : In originally Geant4 code this calculatin is wrong! 
  G4int nion = G4int(G4RandGauss::shoot(meanion,sig) + 0.5);
  return nion;
 
}