BOSS_7.0.2: /home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Simulation/BOOST/TofSim/TofSim-00-02-33/src/BesTofSD.cc Source File

00001 //---------------------------------------------------------------------------//
00002 //      BOOST --- BESIII Object_Oriented Simulation Tool                     //
00003 //---------------------------------------------------------------------------//
00004 //Description:
00005 //Author: Dengzy
00006 //Created: Mar, 2004
00007 //Modified:
00008 //Comment:
00009 //---------------------------------------------------------------------------//
00010 //$Id: BesTofSD.cc
00011 
00012 #include "BesTofSD.hh"
00013 #include "ReadBoostRoot.hh"
00014 #include "G4HCofThisEvent.hh"
00015 #include "G4Step.hh"
00016 #include "G4ThreeVector.hh"
00017 #include "G4SDManager.hh"
00018 #include "G4ios.hh"
00019 #include "G4UnitsTable.hh"
00020 #include "G4Event.hh"
00021 #include "GaudiKernel/ISvcLocator.h"
00022 #include "GaudiKernel/Bootstrap.h"
00023 #include "GaudiKernel/MsgStream.h"
00024 #include "G4LossTableManager.hh"
00025 #include "G4ElectronIonPair.hh"
00026 
00027 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00028 
00029 BesTofSD::BesTofSD( G4String name ) : BesSensitiveDetector(name), m_besTofList(0) {
00030     collectionName.insert("BesTofHitsCollection");
00031     collectionName.insert("BesTofHitsList");
00032     elIonPair = new G4ElectronIonPair();
00033 }
00034 
00035 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00036 
00037 BesTofSD::~BesTofSD() {
00038 }
00039 
00040 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00041 
00042 void BesTofSD::Initialize( G4HCofThisEvent* ) {
00043     m_besTofCollection = new BesTofHitsCollection( SensitiveDetectorName, collectionName[0] );
00044     //elIonPair = G4LossTableManager::Instance()->ElectronIonPair();
00045 }
00046 
00047 //for MC Truth
00048 void BesTofSD::BeginOfTruthEvent( const G4Event* evt) {
00049     m_event = evt->GetEventID();
00050     m_besTofList = new BesTofHitsCollection( SensitiveDetectorName, collectionName[1] );
00051     m_trackIndexes.clear();
00052     m_trackIndex = -99;
00053 }
00054 
00055 void BesTofSD::EndOfTruthEvent( const G4Event* evt ) {
00056     static G4int HLID=-1;
00057     if( HLID<0 ) {
00058         HLID = G4SDManager::GetSDMpointer()->GetCollectionID( collectionName[1] );
00059     }
00060     G4HCofThisEvent* HCE = evt->GetHCofThisEvent();
00061     HCE->AddHitsCollection( HLID, m_besTofList );
00062 }
00063 
00064 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00065 
00066 G4bool BesTofSD::ProcessHits( G4Step* aStep, G4TouchableHistory* ) {
00067 
00068     G4double chg=aStep->GetTrack()->GetDefinition()->GetPDGCharge();
00069     G4double edep = aStep->GetTotalEnergyDeposit();
00070     G4double stepL=aStep->GetStepLength();
00071     G4double deltaT=aStep->GetDeltaTime();
00072     G4StepPoint* preStep = aStep->GetPreStepPoint();
00073     G4ThreeVector pDirection=preStep->GetMomentumDirection();
00074     G4String particleName = aStep->GetTrack()->GetDefinition()->GetParticleName();
00075     G4Material* scinMaterial = aStep->GetTrack()->GetMaterial();
00076     G4double charge = aStep->GetTrack()->GetDefinition()->GetPDGCharge();
00077     G4int pdgcode = aStep->GetTrack()->GetDefinition()->GetPDGEncoding();
00078 
00079     if( chg==0 && edep==0 && stepL==0 ) { return false; }
00080 
00081     BesTofHit *newHit = new BesTofHit();
00082     G4int trackId = aStep->GetTrack()->GetTrackID();
00083 
00084     newHit->SetEvent(m_event);
00085     newHit->SetTrackIndex(-99);
00086     newHit->SetG4Index(aStep->GetTrack()->GetTrackID());
00087     newHit->SetEdep(edep);
00088     newHit->SetStepL(stepL);
00089     //newHit->SetPName(particleName);
00090     newHit->SetTrackL(aStep->GetTrack()->GetTrackLength());
00091     G4ThreeVector pos=preStep->GetPosition();
00092     newHit->SetPos(pos);
00093     G4double globalTime=preStep->GetGlobalTime();
00094     newHit->SetTime(globalTime);
00095     newHit->SetDeltaT(deltaT);
00096     newHit->SetPDirection(pDirection);
00097     newHit->SetMomentum(preStep->GetMomentum());
00098     //newHit->SetMaterial(scinMaterial);
00099     newHit->SetCharge(charge);
00100     newHit->SetPDGcode(pdgcode);
00101 
00102     //Get local position in sensitive detector, add by anff
00103     G4ThreeVector locPos(0,0,0);
00104     G4TouchableHistory* theTouchable = (G4TouchableHistory*)( preStep->GetTouchable() );
00105     //partId: barrel(1), east endcap(0), west endcap(2);
00106     //scinNb: barrel(0-175),east endcap(0-47), west endcap(0-47)
00107     //partId: mrpc: east(3), west(4), 
00108 
00109     G4String name;
00110     //for(G4int i=0; i<theTouchable->GetHistoryDepth(); i++)
00111     //{
00112     //    name = theTouchable->GetVolume(i)->GetName();    
00113     //    G4cout<<"********* name of the "<<i<<" volume:  "<<name<<" ***************"<<G4endl;
00114     //}
00115     name = theTouchable->GetVolume(0)->GetName();
00116 
00117     G4int partId=-1, scinNb=-1, gapNb=-1, number=-1;
00118     G4int strip = -1;
00119     gapNb = theTouchable->GetReplicaNumber(0);
00120     number = theTouchable->GetReplicaNumber(2); //This is valid only for scintillator
00121 
00122     //MRPC ENDCAPS construct by geant4 code
00123     if(name.contains("physical_gasLayer"))
00124     {
00125         locPos = theTouchable->GetHistory()->GetTopTransform().TransformPoint(pos);
00126         number = theTouchable->GetReplicaNumber(3);
00127         scinNb = number;
00128 
00129         G4String name1 = theTouchable->GetVolume(4)->GetName();
00130         if(name1 == "physicalEcTofEast") partId=3;
00131         else if(name1 == "physicalEcTofWest") partId=4;
00132         //G4cout<<"scinNb= "<<scinNb<<"  gapNb= "<<gapNb<<"  partId="<<partId<<G4endl;
00133     }
00134     //MRPC ENDCAPS construct by GDML
00135     else if(name=="logical_gasLayer")
00136     {
00137         locPos = theTouchable->GetHistory()->GetTopTransform().TransformPoint(pos);
00138         number = theTouchable->GetReplicaNumber(3);
00139         scinNb = 35-number;
00140 
00141         G4String name1 = theTouchable->GetVolume(4)->GetName();
00142         if(name1 == "logicalEcTofEast") partId=3;
00143         else if(name1 == "logicalEcTofWest") partId=4;
00144         //G4cout<<"scinNb= "<<scinNb<<"  gapNb= "<<gapNb<<"  partId="<<partId<<G4endl;
00145     }
00146 
00147 
00148     // Scintillator
00149     else if( name=="physicalScinBr1" ) {
00150         partId = 1;
00151         scinNb = number;
00152         //G4cout<<"scinNb= "<<scinNb<<"  partId="<<partId<<G4endl;
00153     }
00154     else if( name=="physicalScinBr2" ) {
00155         partId = 1;
00156         scinNb = number+88;
00157         //G4cout<<"scinNb= "<<scinNb<<"  partId="<<partId<<G4endl;
00158     }
00159     else if( name=="physicalScinEcWest" ) {
00160         partId = 2;
00161         scinNb = number;
00162         //G4cout<<"scinNb= "<<scinNb<<"  partId="<<partId<<G4endl;
00163     }
00164     else if( name=="physicalScinEcEast" ) {
00165         partId = 0;
00166         scinNb = number;
00167         //G4cout<<"scinNb= "<<scinNb<<"  partId="<<partId<<G4endl;
00168     }
00169 
00170 
00171     //if construct by gdml files
00172     else if( name=="logicalScinBr1" || name=="logicalScinBr2" ) {
00173         partId = 1;
00174         scinNb = (527-number)/3;
00175     }
00176     //copy number: east:1-95(scinNb:47-0), west:1-95(scinNb:47-0)
00177     else if( name=="logicalScinEcEast" ) {
00178         partId = 0;
00179         scinNb = (95-number)/2;
00180     }
00181     else if( name=="logicalScinEcWest" ) {
00182         partId  = 2;
00183         scinNb = (95-number)/2;
00184     }
00185     else { return false; }
00186 
00187 
00188     if(name.contains("physical_gasLayer") || name.contains("logical_gasLayer"))
00189     {
00190         G4double zz = locPos.z()-0.5*mm+(24+3)*mm*6; //0.5mm is the offset between strip and gasLayer
00191         if(zz<=0)
00192         {
00193             strip=0;
00194         }
00195         else if(zz>0 && zz<12*27*mm)
00196         {
00197             for(G4int i=0; i<12; i++)
00198             {
00199                 if(zz>i*27*mm && zz<=(i+1)*27*mm)
00200                 {
00201                     strip = i;
00202                     //G4cout<<"Local z is "<<locPos.z()<<"  strip number is: "<<strip<<" !!!"<<G4endl;        
00203                     break;
00204                 }
00205             }
00206         }
00207         else
00208         {
00209             strip=11;
00210         }
00211         if(strip<0) strip=0;
00212         if(strip>11) strip=11;
00213     }
00214 
00215 
00216     newHit->SetPartId(partId);
00217     newHit->SetScinNb(scinNb);
00218     newHit->SetGapNb(gapNb);
00219     newHit->SetLocPos(locPos);
00220     newHit->SetModule(scinNb);
00221     newHit->SetStrip(strip);
00222     //newHit->SetIons(elIonPair->SampleNumberOfIonsAlongStep(aStep));
00223     newHit->SetIons(SampleNumberOfIonsAlongStep(aStep,elIonPair));
00224     //newHit->Draw();
00225 
00226     G4int trackIndex, g4TrackId;
00227     GetCurrentTrackIndex( trackIndex, g4TrackId );
00228     newHit->SetTrackIndex( trackIndex );
00229     //newHit->Print();
00230     if( edep>0 ) {
00231         m_besTofCollection->insert( newHit );
00232     }
00233 
00234     //for mc truth
00235     if( m_besTofList ) {
00236         G4int trackIndex, g4TrackId;
00237         GetCurrentTrackIndex(trackIndex, g4TrackId);
00238         newHit->SetTrackIndex(trackIndex);
00239         if( m_trackIndex != trackIndex ) {
00240             m_trackIndex = trackIndex;
00241             //G4int size = m_trackIndexes.size();
00242             //G4int flag=1;
00243             //if(size>0)
00244             //{
00245             //  for(G4int i=0;i<size;i++)
00246             //    if(m_trackIndexes[i] == trackIndex )
00247             //    {flag =0; break;}
00248             //}
00249             //if(flag && aStep->GetTrack()->GetTrackID()==g4TrackId )
00250             G4int flag = 1;
00251             G4int pdg = abs(aStep->GetTrack()->GetDefinition()->GetPDGEncoding());
00252             if( pdg==12 || pdg==14 || pdg==16 ) { flag=0; }
00253             if( flag && aStep->GetTrack()->GetTrackID()==g4TrackId ) {
00254                 m_trackIndexes.push_back(trackIndex);
00255                 BesTofHit* truHit = new BesTofHit();
00256                 *truHit = *newHit;
00257                 m_besTofList->insert(truHit);
00258             }
00259         }
00260     }
00261     if( edep<=0 ) { delete newHit; }
00262 
00263     return true;
00264 
00265 }//Close Process Hits
00266 
00267 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00268 
00269 void BesTofSD::EndOfEvent( G4HCofThisEvent* HCE ) {
00270     static G4int HCID=-1;
00271     if( HCID<0 ) {
00272         HCID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);
00273     }
00274     HCE->AddHitsCollection(HCID,m_besTofCollection);
00275 }
00276 
00277 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00278 
00279 G4int BesTofSD::SampleNumberOfIonsAlongStep( const G4Step* step, G4ElectronIonPair* eipair ) {
00280     G4double FanoFactor = 0.2;  //Is this value appropiate?
00281     G4double meanion = eipair->MeanNumberOfIonsAlongStep(step);
00282     G4double sig = std::sqrt(FanoFactor*meanion);  //Be carefull : In originally Geant4 code this calculatin is wrong! 
00283     G4int nion = G4int(G4RandGauss::shoot(meanion,sig) + 0.5);
00284     return nion;
00285 }
00286 
00287