BOSS_7.0.2: /home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Simulation/BOOST/MucSim/MucSim-00-01-03/src/BesMucSD.cc Source File

00001 //---------------------------------------------------------------------------//
00002 //      BOOST --- BESIII Object_Oriented Simulation Tool                     //
00003 //---------------------------------------------------------------------------//
00004 //Description:
00005 //Author:  Youzy      Peking University      mail: youzy@hep.pku.cn
00006 //Created: Nov, 2003
00007 //Modified:
00008 //Comment:
00009 //  2006.10.17 update to new gdml, structure changed
00010 //  new gdml structure: gap<---aluminumbox<---gaschamber  
00011 //  but still use the former simulation structure  
00012 //---------------------------------------------------------------------------//
00013 
00014 //
00015 // $Id: BesMucSD.cc,v 1.13 2009/08/25 13:33:54 xieyg Exp $
00016 // GEANT4 tag $Name: MucSim-00-01-03 $
00017 //
00018 
00019 #include "GaudiKernel/Bootstrap.h"
00020 #include "GaudiKernel/ISvcLocator.h"
00021 #include "G4Svc/G4Svc.h"
00022 
00023 #include "BesMucSD.hh"
00024 #include "BesMucDigit.hh"
00025 #include "G4HCofThisEvent.hh"
00026 #include "G4Step.hh"
00027 #include "G4ThreeVector.hh"
00028 #include "G4SDManager.hh"
00029 #include "G4EventManager.hh"
00030 #include "G4ios.hh"
00031 #include "G4UnitsTable.hh"
00032 #include "BesMucEfficiency.hh"
00033 #include "BesMucNoise.hh"
00034 #include "ReadBoostRoot.hh"
00035 #include "Randomize.hh"
00036 #include "strstream"
00037 #include "G4Box.hh"
00038 #include "stdlib.h"
00039 
00040   BesMucSD::BesMucSD(G4String name,BesMucConstruction* Det)
00041 :BesSensitiveDetector(name),detector(Det)
00042 {
00043   collectionName.insert("BesMucHitsCollection");
00044   collectionName.insert("BesMucHitsList");
00045   HitID = new G4int[500];
00046 
00047   //liangyt 2006.3.1 initialize muc-efficiency;
00048   G4String GeometryPath = ReadBoostRoot::GetBoostRoot(); 
00049   G4String GeometryPath2 = GeometryPath;
00050   if(!GeometryPath){ 
00051     G4Exception("BOOST environment not set!");
00052   }
00053   GeometryPath += "/dat/muc-effi.dat";
00054   GeometryPath2 += "/dat/muc-noise.dat";
00055 
00056   m_effi=BesMucEfficiency::Instance();
00057   //  m_effi->Initialize(GeometryPath);  
00058   
00059   //retrieve G4Svc
00060   ISvcLocator* svcLocator = Gaudi::svcLocator();
00061   IG4Svc* iG4Svc;
00062   StatusCode sc=svcLocator->service("G4Svc", iG4Svc);
00063   m_G4Svc = dynamic_cast<G4Svc *>(iG4Svc);
00064   m_noiseMode = m_G4Svc->MucNoiseMode();  
00065   G4cout << "MucNoiseMode:\t"<<m_noiseMode<<G4endl;
00066   
00067   if( m_noiseMode != 0 ) 
00068   {
00069     m_noise=BesMucNoise::Instance();
00070     G4LogicalVolume* logicalMuc = detector->GetPhysicalMuc()->GetLogicalVolume();
00071     m_noise->Initialize(GeometryPath2,logicalMuc);
00072     //m_noise->Initialize(GeometryPath2,logicalMuc,"ROOT");
00073   }
00074 
00075   m_PreviousPrimaryTrackG4Id = 0;
00076   m_CurEvent = 0;
00077   m_TrackCon = 0;
00078 
00079 }
00080 
00081 BesMucSD::~BesMucSD(){ 
00082   delete [] HitID;
00083 }
00084 
00085 void BesMucSD::Initialize(G4HCofThisEvent*)
00086 {
00087   MucHitCollection = new BesMucHitsCollection(SensitiveDetectorName, collectionName[0]); 
00088   // for (G4int j=0;j<detector->GetBesMucNbOfTraps();j++) {HitID[j] = -1;};  
00089   TotEneDeposit = 0;
00090   // G4cout<<"----------in SD::Init()---  "<<detector->GetPhysicalMuc()->GetLogicalVolume()->GetName()<<G4endl;  
00091 }
00092 
00093 void BesMucSD::BeginOfTruthEvent(const G4Event* )
00094 {
00095   MucHitList = new BesMucHitsCollection(SensitiveDetectorName, collectionName[1]);
00096   m_trackIndex = -99;
00097   m_trackIndexes.clear();
00098   m_prePart  = -99;
00099   m_preSeg   = -99;
00100   m_preGap   = -99;
00101   m_preStrip = -99;
00102   //G4cout<<"---in BesMucSD::BeginOfTruthEvent()---"<<MucHitCollection->entries()<<" "<<MucHitList->entries()<<G4endl;
00103   
00104   //add oise
00105   if(m_noiseMode != 0) m_noise->AddNoise(m_noiseMode, MucHitCollection, MucHitList);  //commen out 2006.10.21
00106 }
00107 
00108 void BesMucSD::EndOfTruthEvent(const G4Event* evt)
00109 {
00110   static G4int HLID=-1;
00111   if(HLID<0)
00112     HLID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[1]);
00113   G4HCofThisEvent* HCE = evt->GetHCofThisEvent();
00114   HCE->AddHitsCollection(HLID,MucHitList);
00115 }
00116 
00117 G4bool BesMucSD::ProcessHits(G4Step* aStep,G4TouchableHistory*)
00118 {
00119   G4Track *curTrack = aStep->GetTrack();
00120 
00121   m_CurEvent = G4EventManager::GetEventManager()->GetConstCurrentEvent();
00122   //  if(m_CurEvent->GetEventID()<19565)return true;  // ---del
00123   if (m_CurEvent) {
00124     m_TrackCon = m_CurEvent->GetTrajectoryContainer();
00125 
00126     //  if (!m_TrackCon) G4cout << "m_TrackCon does not exist " << G4endl;
00127     //  else {
00128     //  G4cout << "m_TrackCon size " << m_TrackCon->size() << G4endl;
00129     //  for (int i = 0; i < (int)m_TrackCon->size(); i++) {
00130     //  G4cout << "track " << i << " g4TrackId " << (*m_TrackCon)[i]->GetTrackID() << " parent ID " << (*m_TrackCon)[i]->GetParentID() << G4endl;
00131     //  }
00132     //  }
00133   }
00134 
00135 
00136   if (curTrack->GetDefinition()->GetPDGCharge() == 0.) return false;
00137 
00138   G4double edep = aStep->GetTotalEnergyDeposit();
00139   //if(edep == 0.) return false;
00140 
00141   //TotEneDeposit+=edep;
00142   G4int trackIndex = -99, g4TrackId = -99;
00143   GetCurrentTrackIndex(trackIndex, g4TrackId);
00144 
00145   G4TouchableHistory* theTouchable
00146     = (G4TouchableHistory*)(aStep->GetPreStepPoint()->GetTouchable());
00147 
00148   BesMucHit *newHit = new BesMucHit();
00149 
00150   G4int trackID = curTrack->GetTrackID(); //G4 track ID of current track.
00151   G4int parentID = curTrack->GetParentID(); //G4 track ID of parent track.
00152   newHit->SetTrackID(trackID);
00153   newHit->SetTrackIndex(trackIndex); // MC truth track index
00154 
00155   G4int pdg = curTrack->GetDefinition()->GetPDGEncoding();
00156   newHit->SetPDGCode(pdg);
00157 
00158   newHit->SetEdep(edep);
00159 
00160   G4ThreeVector pos = 0.5*( aStep->GetPostStepPoint()->GetPosition()
00161       + aStep->GetPreStepPoint()->GetPosition() );
00162   newHit->SetPos(pos);
00163 
00164   G4ThreeVector posInGas = theTouchable->GetHistory()->GetTopTransform().TransformPoint(pos);
00165   G4int       stackDepth = theTouchable->GetHistory()->GetDepth();
00166   G4ThreeVector posInBox = theTouchable->GetHistory()->GetTransform(stackDepth-1).TransformPoint(pos);
00167   newHit->SetPosLocal(posInBox);
00168 
00169   G4ThreeVector posInGap = theTouchable->GetHistory()->GetTransform(stackDepth-2).TransformPoint(pos);
00170   // cout<<"  pos  "<<pos.x()<<" "<<pos.y()<<"  "<<pos.z()<<"depth= "<<stackDepth<<endl
00171   //       <<" posInBox(3)"<<posInBox.x()<<"  "<<posInBox.y()<<" "<<posInBox.z()<<endl
00172   //       <<" posInGap(2)"<<posInGap.x()<<"  "<<posInGap.y()<<" "<<posInGap.z()<<endl;
00173 
00174   G4double energy = aStep->GetPreStepPoint()->GetKineticEnergy();
00175   newHit->SetEnergy(energy);
00176 
00177   G4ThreeVector dir = aStep->GetPreStepPoint()->GetMomentumDirection();
00178   newHit->SetDir(dir);
00179 
00180   G4ThreeVector momentum = aStep->GetPreStepPoint()->GetMomentum();
00181   newHit->SetMomentum(momentum);
00182 
00183   G4double GlobalTime = aStep->GetPostStepPoint()->GetGlobalTime();
00184   newHit->SetTime(GlobalTime);
00185 
00186   G4VPhysicalVolume* vl = theTouchable->GetVolume(0);
00187   newHit->SetVolume(vl);
00188 
00189   //G4cout<<"in SD "<<newHit->GetPart()<<" "<<newHit->GetSeg()<<" "<<newHit->GetGap()<<" "<<newHit->GetGasChamber()<<" "<<newHit->GetPanel()<<endl;
00190 
00191   newHit->Draw();
00192   //newHit->Print();
00193   //MucHitCollection->insert(newHit);
00194   m_PreviousPrimaryTrackG4Id = g4TrackId;
00195 
00196   //-----------add by liangyt for efficiency
00197 
00198   //for MC Truth
00199   if (MucHitList) {
00200 
00201     // cout <<"g4TrackId " << g4TrackId << " trackID " << trackID << " parentID " << parentID << "trackIndex " << trackIndex << endl; 
00202     // g4TrackId:  G4 track ID of previous primary track in trackList,
00203     // trackID:    G4 track ID of current track.
00204     // parentID:   G4 track ID of parent track of current track.
00205     // trackIndex: track ID in MC truth
00206 
00207     // no longer necessary to find the track's ancestor primary track,
00208     // it is kept in trackIndex untill next primary track generated. So IsChildof is uncessary
00209     //if ( g4TrackId == trackID ) { //|| IsChildOf(curTrack, g4TrackId) ) {     // This track is the primary track & will appear in MC truth
00210     G4int newTrackFlag = 0;
00211     newHit->SetTrackIndex(trackIndex);
00212     if(m_trackIndex != trackIndex) {
00213       m_trackIndex = trackIndex;
00214       G4int size = m_trackIndexes.size();
00215       newTrackFlag = 1;
00216       if (size > 0) {
00217         for(G4int i=0;i<size;i++)
00218           if(m_trackIndexes[i] == trackIndex ) {
00219             newTrackFlag = 0; 
00220             break;
00221           }
00222       }
00223     }
00224 
00225     if (newTrackFlag) {
00226       m_trackIndexes.push_back(trackIndex);
00227       m_prePart  = -99;
00228       m_preSeg   = -99;
00229       m_preGap   = -99;
00230       m_preStrip = -99;
00231     }
00232     BesMucHit* truHit = new BesMucHit();
00233     *truHit = *newHit;
00234     if (g4TrackId != trackID) {
00235       //        trackIndex += 1000; // a sencondary track
00236       trackIndex += 0; //2006.12.22 do not indicate secondary track now
00237       truHit->SetTrackIndex(trackIndex);
00238     }
00239     //cout << "trackIndex " << trackIndex << endl;
00240 
00241     BesMucDigit aDigit;
00242     aDigit.SetHit(truHit);
00243     G4int curPart, curSeg, curGap, curStrip;
00244     curPart  = aDigit.GetPart();
00245     curSeg   = aDigit.GetSeg();
00246     curGap   = aDigit.GetGap();
00247     curStrip = aDigit.GetNearestStripNo();
00248 
00249     m_effi->CheckCalibSvc();
00250     m_effi->SetHit(truHit);
00251     G4double need_eff = m_effi->GetEfficiency();
00252 
00253     //G4cout<<"in SD effi= "<<need_eff<<endl;
00254     //need_eff = 1.0;  //2006.12.28
00255     if (curPart == m_prePart && curSeg == m_preSeg &&
00256         curGap  == m_preGap && curStrip == m_preStrip) {
00257       //cout<<MucHitList->entries()<<" "<<MucHitCollection->entries()<<" "<<need_eff<<"---if----curPart "<<curPart<<"curSeg "<<curSeg<<"curGap "<<curGap<<"curStrip "<<curStrip<<" momentum "<<momentum.x()<<" "<<momentum.y()<<" "<<momentum.z()<<" "<<endl;
00258       delete truHit;delete newHit;
00259     }
00260     else {
00261       //cout<<MucHitList->entries()<<"----else--Part "<<curPart<<" Seg "<<curSeg<<" Gap "<<curGap<<" Strip "<<curStrip<<endl;
00262       truHit->SetPart(curPart);
00263       truHit->SetSeg(curSeg);
00264       truHit->SetGap(curGap);
00265       truHit->SetStrip(curStrip);
00266 
00267       // if a truHit with the same id(part, seg, gap, strip) and trackIndex(%1000) exist,
00268       // they belong to the same primary track,(maybe one is primary, the other is secondary), dont add.
00269       bool truHitExist = false;
00270       G4int n_hit = MucHitList->entries();
00271       for(G4int iTru=0;iTru<n_hit;iTru++) {
00272         BesMucHit* aTruHit = (*MucHitList)[iTru];
00273         if ( aTruHit->GetTrackIndex()%1000 == truHit->GetTrackIndex()%1000 &&
00274             aTruHit->GetPart()  == truHit->GetPart() &&
00275             aTruHit->GetSeg()   == truHit->GetSeg()  &&
00276             aTruHit->GetGap()   == truHit->GetGap()  &&
00277             aTruHit->GetStrip() == truHit->GetStrip() ) 
00278         {
00279           truHitExist = true;
00280           break;
00281         }
00282       }
00283       G4float random=G4UniformRand();   //*** use other random
00284       // G4float random=(rand()%100000)/100000.0;
00285       // G4cout<<"---in SD---"<<random<<endl;
00286       if (random<=need_eff){ MucHitCollection->insert(newHit);}
00287       else delete newHit;
00288       if (!truHitExist&&random<=need_eff)
00289       { MucHitList->insert(truHit);}
00290       else delete truHit;
00291       m_prePart  = curPart;
00292       m_preSeg   = curSeg;
00293       m_preGap   = curGap;
00294       m_preStrip = curStrip;
00295     }
00296     //}
00297   }
00298 
00299   return true;
00300   }
00301 
00302   bool BesMucSD::IsChildOf(G4Track *curTrack, G4int primaryG4TrackID) 
00303   {
00304     G4cout << "IsChildof " << "curTrackID " <<  curTrack->GetTrackID() << G4endl;
00305 
00306     G4VTrajectory* aTraj = GetTrajFromID(curTrack->GetParentID());
00307     G4cout << "IsChildof " << "parentTrackID " <<  aTraj->GetTrackID() << G4endl;
00308 
00309     while ( aTraj->GetTrackID() != 0 && aTraj->GetTrackID() != primaryG4TrackID ) {
00310       //G4cout << "loop: parentID " << aTraj->GetParentID() << G4endl;
00311       aTraj = GetTrajFromID(aTraj->GetParentID());
00312     }
00313 
00314     if (aTraj->GetTrackID() == primaryG4TrackID) return true;
00315     else return false;
00316   }
00317 
00318   G4VTrajectory* BesMucSD::GetTrajFromID(G4int id) 
00319   {
00320     // TrajContainer does not contain current track
00321     G4cout << "begin of GetTrajFromID, id : " << id << G4endl;
00322 
00323     if (!m_TrackCon) G4cout << "Trajectory not saved? Set /tracking/storeTrajectory 1" << G4endl;
00324 
00325     if (m_TrackCon->size() == 0) {
00326       G4cout << "BesMucSD::GetTrajFromID, TrackCon size 0" << G4endl;
00327       return 0;
00328     }
00329 
00330     int k = 0;
00331     while( k < (int)m_TrackCon->size() && (*m_TrackCon)[k]->GetTrackID() != id ) {
00332       G4cout << "GetTrajFromID " << k << " : ID " << (*m_TrackCon)[k]->GetTrackID() << G4endl;
00333       k++;
00334       if (!(*m_TrackCon)[k]) {
00335         G4cout << "G4Track ID " << (*m_TrackCon)[k]->GetTrackID() << " doesnt exist in TrajContainer of this event! " << G4endl; 
00336         return 0;
00337       }
00338     }
00339 
00340     if ( k == (int)m_TrackCon->size() ) {
00341       G4cout << "BesMucSD::GetTrajFromID, track with ID " << id << " not found" << G4endl;
00342       return 0;
00343     }
00344 
00345     return (*m_TrackCon)[k];
00346   }
00347 
00348   void BesMucSD::EndOfEvent(G4HCofThisEvent* HCE)
00349   {
00350     static G4int HCID=-1;
00351     if(HCID<0)
00352     {HCID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]); }
00353     HCE->AddHitsCollection(HCID, MucHitCollection);
00354   }
00355