BesMcTruthWriter.cc

Go to the documentation of this file.

//---------------------------------------------------------------------------//
////      BOOST --- BESIII Object_Oriented Simulation Tool                     //
////---------------------------------------------------------------------------//
////Description:
////Author:   Dengzy
////Created:  Mar, 2004
////Modified:
////Comment:
//
#include "G4DigiManager.hh"
#include "BesMdcHit.hh"
#include "BesCgemHit.hh"
#include "BesTofHit.hh"
#include "BesEmcHit.hh"
#include "BesMucHit.hh"
#include "BesTruthTrack.hh"
#include "BesTruthVertex.hh"
#include "BesSensitiveManager.hh"
//#include "G4HCofThisEvent.hh"
//#include "G4SDManager.hh"
//#include "G4PrimaryVertex.hh"
//#include "G4PrimaryParticle.hh"
 
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/RegistryEntry.h"
#include "GaudiKernel/MsgStream.h"
#include "CLHEP/Vector/LorentzVector.h"
#include "CLHEP/Geometry/Point3D.h"
 
#include "McTruth/DecayMode.h"
#include "McTruth/MdcMcHit.h"
#include "McTruth/CgemMcHit.h"
#include "McTruth/TofMcHit.h"
#include "McTruth/EmcMcHit.h"
#include "McTruth/MucMcHit.h"
 
#include "Identifier/Identifier.h"
#include "Identifier/MdcID.h"
#include "Identifier/CgemID.h"
#include "Identifier/TofID.h"
#include "Identifier/EmcID.h"
#include "Identifier/MucID.h"
 
 
#include "TofSim/BesTofDigitizerEcV4.hh"
 
#include "McTruth/McEvent.h"
#include "EventModel/EventModel.h"
 
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/IJobOptionsSvc.h"
#include "BesMcTruthWriter.hh"
 
#include "TArrayI.h"
 
BesMcTruthWriter::BesMcTruthWriter()
{
  m_DigiMan = G4DigiManager::GetDMpointer(); 
  //mdcGeoPointer=BesMdcGeoParameter::GetGeo(); 
  
  // set m_cgem_misAligned
  m_cgem_misAligned=0;
  ISvcLocator* svcLocator = Gaudi::svcLocator();
  static IJobOptionsSvc* jobSvc = 0;
  StatusCode sc = svcLocator->service("JobOptionsSvc", jobSvc);
  if ( sc.isFailure() ) {
      std::cout << "BesMcTruthWriter: Can't get the JobOptionsSvc " << std::endl;
      //return sc;
  }
  else 
  {
      const vector<const Property*>* properties = jobSvc->getProperties("BesSim");
      if ( properties != NULL ) {
          for ( unsigned int i = 0; i < properties->size(); ++i ) {
              if ( properties->at(i)->name() == "CgemMisAligned" ) {
                  string strRnd = properties->at(i)->toString();
                  sscanf(strRnd.c_str(), "%i", &m_cgem_misAligned);
                  cout<<"BesMcTruthWriter: cgem_misAligned is set to "<<m_cgem_misAligned<<" by jobOptions"<<endl;
                  break;
              }
          }
      }
  }
 
 
}
 
BesMcTruthWriter::~BesMcTruthWriter()
{
}
 
void BesMcTruthWriter::SaveMcTruth()
{
  //interface to event data service
  ISvcLocator* svcLocator = Gaudi::svcLocator();
  StatusCode sc=svcLocator->service("EventDataSvc", m_evtSvc);
  if (sc.isFailure())
    G4cout<<"Could not accesss EventDataSvc!"<<G4endl;
  
  DataObject *aMcEvent;
  m_evtSvc->findObject("/Event/MC",aMcEvent);
  if(aMcEvent==NULL) {       
    McEvent* aMcEvent = new McEvent;
    sc = m_evtSvc->registerObject("/Event/MC",aMcEvent);
    if(sc!=StatusCode::SUCCESS) 
      G4cout<< "Could not register McEvent" <<G4endl;
  }
 
  SaveMcParticle();
  SaveDecayMode();
  SaveMdcTruth();
  SaveCgemTruth();
  SaveTofTruth();
  SaveEmcTruth();
  SaveMucTruth();
}
 
void BesMcTruthWriter::SaveMcParticle()
{
  BesSensitiveManager* sensitiveManager = BesSensitiveManager::GetSensitiveManager(); 
  vector<BesTruthTrack*>* trackList =  sensitiveManager->GetTrackList();
  G4int nTrack = trackList->size();
  BesTruthTrack* track;
 
  vector<BesTruthVertex*>* vertexList =  sensitiveManager->GetVertexList();
  G4int nVertex = vertexList->size();
  BesTruthVertex* vertex;
  
  //arrange TruthTrack in trackList in order of trackIndex
  for(int i=0;i<nTrack-1;i++)
  for(int j=i+1;j<nTrack;j++)
    if((*trackList)[i]->GetIndex()>(*trackList)[j]->GetIndex())
    {
      track=(*trackList)[i];
      (*trackList)[i]=(*trackList)[j];
      (*trackList)[j]=track;
    }
  
  Event::McParticleCol *particleCol = new Event::McParticleCol;
  
  //loop over tracks
  for(int i=0;i<nTrack;i++)
  {
    track = (*trackList) [i];
    
    // find out the start point
    bool isPrimary = false;
    bool startPointFound = false;
    BesTruthVertex* startPoint;
    
    for (int i=0;i<nVertex;i++) 
    {
      vertex = (*vertexList) [i];
      if ( track->GetVertex()->GetIndex() == vertex->GetIndex() ) 
      {
        if (!vertex->GetParentTrack()) isPrimary = true;
        startPointFound = true;
        startPoint = vertex;
        break;
      }
    }
 
    if (!startPointFound)  
      std::cout << "error in finding the start point of a track" <<std::endl;
 
 
    bool endPointFound = false;
    // find out the end point
    for (int i=0;i<nVertex;i++)
    {
      vertex = (*vertexList) [i];
      if( track->GetTerminalVertex() )
      {
        if (track->GetTerminalVertex()->GetIndex() == vertex->GetIndex() ) 
        {
            //create the mc particle
          Event::McParticle* mcParticle = new Event::McParticle;
        
            // initialization 
          HepLorentzVector initialMomentum(track->GetP4().x()/1000., track->GetP4().y()/1000., track->GetP4().z()/1000., track->GetP4().t()/1000.);
                 
          HepLorentzVector  initialPosition(startPoint->GetPosition().x()/10., startPoint->GetPosition().y()/10., startPoint->GetPosition().z()/10.,startPoint->GetTime());
            
            mcParticle->initialize(track->GetPDGCode(), 0, initialMomentum, initialPosition); 
          
          // Set track index
          mcParticle->setTrackIndex( track->GetIndex() );
          
          // Adding status flag  
            if (isPrimary ) mcParticle->addStatusFlag(Event::McParticle::PRIMARY);
 
          if (track->GetDaughterIndexes().size()==0)
            mcParticle->addStatusFlag(Event::McParticle::LEAF);
        
          //std::cout<<"pdg:"<<track->GetPDGCode()<<" "
          //         <<"source:"<<track->GetSource()<<" "
          //         <<std::endl;
          if (track->GetSource()=="FromGenerator")
            mcParticle->addStatusFlag(Event::McParticle::DECAYED);
          else if(track->GetSource()=="FromG4")
            mcParticle->addStatusFlag(Event::McParticle::DECAYFLT); 
          else
            mcParticle->addStatusFlag(Event::McParticle::ERROR);
 
          //std::cout<<"status:"<<mcParticle->statusFlags()<<" "
          //         <<"FGN:"<<mcParticle->decayFromGenerator()<<" "
          //         <<"FG4:"<<mcParticle->decayInFlight()<<std::endl;
            // Adding vertex index 
            mcParticle->setVertexIndex0(startPoint->GetIndex());
          mcParticle->setVertexIndex1(vertex->GetIndex());
          
          // Set the final position 
          HepLorentzVector  finalPosition(vertex->GetPosition().x()/10., vertex->GetPosition().y()/10., vertex->GetPosition().z()/10., vertex->GetTime());
          mcParticle->finalize(finalPosition);
 
          // push back 
            particleCol->push_back(mcParticle);
 
          endPointFound = true; 
        }
      }
    }
   
    if (!endPointFound)
    {
        //std::cout << "warning in finding the end point of a track" <<std::endl;
        if (track->GetDaughterIndexes().size()==0) 
      {
          // create the mc particle
          Event::McParticle* mcParticle = new Event::McParticle;
          // initialization 
          HepLorentzVector initialMomentum( track->GetP4().x()/1000., track->GetP4().y()/1000., track->GetP4().z()/1000., track->GetP4().t()/1000.);
          HepLorentzVector  initialPosition(startPoint->GetPosition().x()/10., startPoint->GetPosition().y()/10., startPoint->GetPosition().z()/10., startPoint->GetTime());
 
          mcParticle->initialize(track->GetPDGCode(), 0, initialMomentum, initialPosition); 
 
        // Set track index
          mcParticle->setTrackIndex( track->GetIndex() );
 
        // Adding status flag  
          mcParticle->addStatusFlag(Event::McParticle::LEAF);
          if (isPrimary ) mcParticle->addStatusFlag(Event::McParticle::PRIMARY);
 
          // Adding vertex index 
          mcParticle->setVertexIndex0( startPoint->GetIndex() );
          mcParticle->setVertexIndex1( -99 );
 
        //std::cout<<"pdg:"<<track->GetPDGCode()<<" "
        //           <<"source:"<<track->GetSource()<<" "
        //           <<std::endl;
        if (track->GetSource()=="FromGenerator")
          mcParticle->addStatusFlag(Event::McParticle::DECAYED);
        else if(track->GetSource()=="FromG4")
          mcParticle->addStatusFlag(Event::McParticle::DECAYFLT); 
        else
          mcParticle->addStatusFlag(Event::McParticle::ERROR);
 
        //std::cout<<"status:"<<mcParticle->statusFlags()<<" "
        //         <<"FGN:"<<mcParticle->decayFromGenerator()<<" "
        //         <<"FG4:"<<mcParticle->decayInFlight()<<std::endl;
 
          // push back 
          particleCol->push_back(mcParticle);
          continue;
      }
    }
    }
 
  // Get primary McParticles 
  SmartRefVector<Event::McParticle> primaryParticleCol;
  Event::McParticleCol::iterator iter_particle;
  for ( iter_particle = particleCol->begin(); 
        iter_particle != particleCol->end(); iter_particle++) {
 
    if ((*iter_particle)->primaryParticle()) {
       Event::McParticle* mcPart =  (Event::McParticle*)(*iter_particle); 
       primaryParticleCol.push_back(mcPart);
     }
  }
 
  if (primaryParticleCol.empty())
    std::cout << "no primary particle found!" << std::endl;
    
  // Add mother particle recursively 
  SmartRefVector<Event::McParticle>::iterator iter_primary; 
  for ( iter_primary = primaryParticleCol.begin(); 
        iter_primary != primaryParticleCol.end(); iter_primary++) {
    if ( !((*iter_primary)->leafParticle()) ) 
      AddMother((*iter_primary), particleCol);
  }
 
  //register McParticle collection to TDS
  StatusCode sc = m_evtSvc->registerObject("/Event/MC/McParticleCol", particleCol);
  if(sc!=StatusCode::SUCCESS) 
    G4cout<< "Could not register McParticle collection" <<G4endl;
  
  //retrive McParticle from TDS
  /*SmartDataPtr<Event::McParticleCol> mcParticleCol(m_evtSvc,"/Event/MC/McParticleCol");
  if(!mcParticleCol)
    G4cout<<"Could not retrieve McParticelCol"<<G4endl;  
  Event::McParticleCol::iterator iter_mc = mcParticleCol->begin();
  for (;iter_mc != mcParticleCol->end(); iter_mc++)
  {
    G4cout<<(*iter_mc)->getTrackIndex();
    G4cout<<" "<<(*iter_mc)->particleProperty();
    G4cout<<" "<<(*iter_mc)->getVertexIndex0();
    G4cout<<" "<<(*iter_mc)->getVertexIndex1();
    G4cout<<" "<<(*iter_mc)->initialFourMomentum().x();
    G4cout<<" "<<(*iter_mc)->initialFourMomentum().y();
    G4cout<<" "<<(*iter_mc)->initialFourMomentum().z();
    G4cout<<" "<<(*iter_mc)->initialFourMomentum().t();
    G4cout<<" "<<(*iter_mc)->initialPosition().x();
    G4cout<<" "<<(*iter_mc)->initialPosition().y();
    G4cout<<" "<<(*iter_mc)->initialPosition().z();
    G4cout<<G4endl;
  }
  G4cout<<"end of retrieve McParticle from TDS"<<G4endl;
  */
 
}
 
void BesMcTruthWriter::AddMother(Event::McParticle* currentParticle, Event::McParticleCol *particleCol)
{
  if (currentParticle->leafParticle()) {
    return;
  }
  Event::McParticleCol::iterator iter;
  bool found = false;
  for ( iter = particleCol->begin(); iter != particleCol->end(); iter++) {
    if (currentParticle->vertexIndex1() == (*iter)->vertexIndex0()) {
      found = true;
      (*iter)->setMother(currentParticle);
      currentParticle->addDaughter(*iter);
      AddMother((*iter), particleCol);
    }
  }
  if (!found) std::cout << "AddMother: inconsistency was found!" << std::endl;
 
}
 
void BesMcTruthWriter::SaveDecayMode()
{
  SmartDataPtr<DecayMode> decayMode(m_evtSvc,"/Event/MC/DecayMode");
  if(!decayMode)
  {
    //G4cout<<"Could not retrieve DecayMode"<<G4endl;
    DecayMode*  aDecayMode = new DecayMode;
    //register Decay Mode to TDS
    int dm[10]={0,0,0,0,0,0,0,0,0,0};
    aDecayMode->putData(dm,10);
    StatusCode scDM = m_evtSvc->registerObject("/Event/MC/DecayMode", aDecayMode);
    if(scDM!=StatusCode::SUCCESS)
      G4cout<< "Could not register DecayMode" <<G4endl;
  }
}
  
void BesMcTruthWriter::SaveMdcTruth()
{
    //Mdc McHit collection defined in BOSS
    Event::MdcMcHitCol*  aMdcMcHitCol = new Event::MdcMcHitCol;
 
    G4int HCID = -1;
    //HCID = m_DigiMan->GetHitsCollectionID("BesMdcTruthCollection");
    HCID = m_DigiMan->GetHitsCollectionID("BesMdcHitsCollection");
    //G4cout<< "SaveMdcTruth(): BesMdcHitsCollection ID="<<HCID <<G4endl;
    if(HCID>=0)
    {
        BesMdcHitsCollection* HC = 0;
        HC = (BesMdcHitsCollection*) (m_DigiMan->GetHitsCollection(HCID));
        G4int n_hit = HC->entries();
        //G4cout<< "SaveMdcTruth(): BesMdcHitsCollection n_hit="<<n_hit <<G4endl;
        if(n_hit>0)
        {
            //arrange hits in hits collection in order of trackIndex
            BesMdcHit* hit;
            /*
            vector<BesMdcHit*>* vecHC = HC->GetVector();
            for(int i=0;i<n_hit-1;i++)
                for(int j=i+1;j<n_hit;j++)
                    if((*vecHC)[i]->GetTrackID()>(*vecHC)[j]->GetTrackID())
                    {
                        hit = (*vecHC)[i];
                        (*vecHC)[i] = (*vecHC)[j];
                        (*vecHC)[j] = hit;
                    }
            */
 
            //push back MdcMcHits to MdcMcHitCol in TDS
            for(G4int i=0;i<n_hit;i++)
            {
                hit = (*HC)[i];
                const Identifier ident =  MdcID::wire_id ( hit->GetLayerNo(), hit->GetCellNo() );
                Event::MdcMcHit* mdcMcHit = new Event::MdcMcHit(ident, hit->GetTrackID(), 
                        hit->GetPos().x(),  hit->GetPos().y(), hit->GetPos().z(), 
                        hit->GetDriftD(), hit->GetEdep(),hit->GetPosFlag() );
                G4ThreeVector p3=hit->GetMomentum();
                mdcMcHit->setMomentumX(p3.x());
                mdcMcHit->setMomentumY(p3.y());
                mdcMcHit->setMomentumZ(p3.z());
                mdcMcHit->setIsSecondary(hit->GetIsSecondary());
                mdcMcHit->setCurrentTrackPID(hit->GetPDGCode());
                mdcMcHit->setFlightLength(hit->GetFlightLength());
                mdcMcHit->setCreatorProcess(hit->GetCreatorProcess());
                mdcMcHit->setDigiIdx(hit->GetDigiId());
                //mdcMcHit->();
                aMdcMcHitCol->push_back(mdcMcHit);
            }      
 
        }
    }
 
    //register MDC McTruth collection to TDS
    StatusCode scMdc = m_evtSvc->registerObject("/Event/MC/MdcMcHitCol", aMdcMcHitCol);
    if(scMdc!=StatusCode::SUCCESS) 
        G4cout<< "Could not register MDC McTruth collection" <<G4endl;
 
    //retrieve MDC McTruth from TDS
    /*SmartDataPtr<Event::MdcMcHitCol> aMcHitCol(m_evtSvc,"/Event/MC/MdcMcHitCol");
      if(!aMcHitCol)
      G4cout<<"Could not retrieve MDC McTruth collection"<<G4endl;
 
      Event::MdcMcHitCol::iterator iMcHitCol;
      for(iMcHitCol=aMcHitCol->begin(); iMcHitCol!=aMcHitCol->end(); iMcHitCol++)
      {
      const Identifier ident = (*iMcHitCol)->identify();
      G4cout<<(*iMcHitCol)->getTrackIndex();
      G4cout<<" "<<MdcID::layer(ident);
      G4cout<<" "<<MdcID::wire(ident);
      G4cout<<" "<<(*iMcHitCol)->getDepositEnergy();
      G4cout<<" "<<(*iMcHitCol)->getDriftDistance();
      G4cout<<" "<<(*iMcHitCol)->getPositionX();
      G4cout<<" "<<(*iMcHitCol)->getPositionY();
      G4cout<<" "<<(*iMcHitCol)->getPositionZ();
      G4cout<<G4endl;
      }
      G4cout<<"end of retrieve MDC McTruth collection"<<G4endl;
      */
}
 
void BesMcTruthWriter::SaveCgemTruth()
{
    //G4cout << "|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||" << G4endl;
    //G4cout << "INFO : BesSim::BesMcTruthWriter::SaveCgemTruth(), Begin! " << G4endl;
    /* Cgem McHit collection defined in BOSS */
    Event::CgemMcHitCol*  aCgemMcHitCol = new Event::CgemMcHitCol;
 
    G4int HCID = -1;
    G4int BESID = -1;
    HCID  = m_DigiMan->GetHitsCollectionID("BesCgemTruthCollection");
    BESID = m_DigiMan->GetHitsCollectionID("BesCgemHitsCollection");
    if (HCID > 0)
    {
        BesCgemHitsCollection *HC = 0;
        HC = (BesCgemHitsCollection*)(m_DigiMan->GetHitsCollection(HCID));
        BesCgemHitsCollection *BES = 0;
        /*Add BEShit for Identifier by sunxh*/
        BES = (BesCgemHitsCollection*)(m_DigiMan->GetHitsCollection(BESID));
        G4int n_hit = HC->entries();
        G4int n_BEShit = BES->entries();
        if (n_hit > 0)
        {
            /* arrange hits in hits collection in order of trackIndex */
            BesCgemHit* hit;
            BesCgemHit* BEShit;
            vector<BesCgemHit*>* vecHC = HC->GetVector();
            for(int i=0; i<n_hit-1; i++)
                for(int j=i+1; j<n_hit; j++)
                    if((*vecHC)[i]->GetTrackID()>(*vecHC)[j]->GetTrackID())
                    {
                        hit = (*vecHC)[i];
                        (*vecHC)[i] = (*vecHC)[j];
                        (*vecHC)[j] = hit;
                    }
 
            ///* push back CgemMcHits to CgemMcHitCol in TDS */
            //for(G4int i=0;i<n_hit;i++)
            //{
            //hit = (*HC)[i];
            //const Identifier ident =  CgemID::chip_id ( hit->GetLayerID(), hit->GetSheetID(),
            //    hit->GetChipID(), hit->GetRow(),
            //    hit->GetCol() );
 
            //ident.show();
            //G4cout<<" "<<ident.get_valueMSB()<<" "<<ident.get_valueLSB()<<" "<<CgemID::layer(ident)<<" "<<CgemID::ladder(ident)
            //<<"   "<<CgemID::chip(ident) <<"  "<<CgemID::row(ident)<<" "<<CgemID::column(ident)<<G4endl;
 
            //Event::CgemMcHit* cgemMcHit = new Event::CgemMcHit(ident, hit->GetTrackID(), 
            //    hit->GetPos().x(),  hit->GetPos().y(), hit->GetPos().z(), hit->GetEdep() );
            //aCgemMcHitCol->push_back(cgemMcHit);
            /*for(G4int j=0;j<n_BEShit;j++)
              {
              BEShit = (*BES)[j];
 
              G4cout << "McWriter track_ID1 = "<< BEShit->GetTrackID() << G4endl;
              G4cout << "McWriter HitID1 = "<< BEShit->GetHitID() << G4endl;
              TArrayI BES_Ident = BEShit->GetIdentifier();
              for(int jj = 0; jj < BES_Ident.GetSize(); jj++){
              G4cout << "McWriter ident1["<<jj<<"] = "<< BES_Ident.GetAt(jj) << G4endl;
              }
              }*/
            /* push back CgemMcHits to CgemMcHitCol in TDS */
            for(G4int i=0;i<n_hit;i++)
            {
                hit = (*HC)[i];
 
                //       cout<<"hitTrack_ID = "<<hit->GetTrackID()<<"\thitPDG = "<<hit->GetPDGCode()<<"\thitParent ID = "<<hit->GetParentID()<<endl;
 
                /* Since no ID in Hit, create McHit without Identifier object */
                Event::CgemMcHit *cgemMcHit = new Event::CgemMcHit(
                        hit->GetTrackID(), hit->GetLayerID() , hit->GetPDGCode(),
                        hit->GetParentID(), hit->GetTotalEnergyDeposit()      ,
                        hit->GetPositionOfPrePoint().x()     , hit->GetPositionOfPrePoint().y()  ,
                        hit->GetPositionOfPrePoint().z()     , hit->GetPositionOfPostPoint().x() ,
                        hit->GetPositionOfPostPoint().y()    , hit->GetPositionOfPostPoint().z() ,
                        hit->GetMomentumOfPrePoint().x()     , hit->GetMomentumOfPrePoint().y()  ,
                        hit->GetMomentumOfPrePoint().z()     , hit->GetMomentumOfPostPoint().x() ,
                        hit->GetMomentumOfPostPoint().y()    , hit->GetMomentumOfPostPoint().z() );
                /*
                   hit->GetPrePositionInCu().x()        , hit->GetPrePositionInCu().y()     ,
                   hit->GetPrePositionInCu().z()        , hit->GetPostPositionInCu().x()    ,
                   hit->GetPostPositionInCu().y()       , hit->GetPostPositionInCu().z()    ,
                   hit->GetMomentumOfCuPre().x()        , hit->GetMomentumOfCuPre().y()     ,
                   hit->GetMomentumOfCuPre().z()        , hit->GetMomentumOfCuPost().x()    ,
                   hit->GetMomentumOfCuPost().y()       , hit->GetMomentumOfCuPost().z()    ); 
                   */
 
                if(m_cgem_misAligned==1) 
                {
                    cgemMcHit->SetPositionXOfPrePoint(hit->GetPositionOfPrePointAlign().x());
                    cgemMcHit->SetPositionYOfPrePoint(hit->GetPositionOfPrePointAlign().y());
                    cgemMcHit->SetPositionZOfPrePoint(hit->GetPositionOfPrePointAlign().z());
                    cgemMcHit->SetPositionXOfPostPoint(hit->GetPositionOfPostPointAlign().x());
                    cgemMcHit->SetPositionYOfPostPoint(hit->GetPositionOfPostPointAlign().y());
                    cgemMcHit->SetPositionZOfPostPoint(hit->GetPositionOfPostPointAlign().z());
                }
                cgemMcHit->SetIsSecondary(hit->GetIsSecondary());
 
                //   G4cout << "McWriter track_ID = "<< hit->GetTrackID() << G4endl;
                int tmp[2000];
                int size(0);
                for(G4int j=0;j<n_BEShit;j++)
                {
                    BEShit = (*BES)[j];
                    //      cout<<"BEShitTrack_ID = "<<BEShit->GetTrackID()<<"\tBEShitPDG = "<<BEShit->GetPDGCode()<<"\tBEShitParent ID = "<<BEShit->GetParentID()<<endl;
                    if(BEShit->GetTrackID() != hit->GetTrackID())continue;
                    TArrayI hit_ID = hit->GetIdentifier();
                    for(int ii = 0; ii < hit_ID.GetSize(); ii++){
                        if(hit_ID.GetAt(ii) != BEShit->GetHitID())continue;
                        //              G4cout << "McWriter EneDp = "<< BEShit->GetTotalEnergyDeposit() << G4endl;
 
                        TArrayI BES_Ident = BEShit->GetIdentifier();
                        for(int jj = 0; jj < BES_Ident.GetSize(); jj++){
                            tmp[size] = BES_Ident.GetAt(jj);
                            size++;
                        }
                    } 
                }
                cgemMcHit->AddIdentifier(tmp,size);
                cgemMcHit->SetCreatorProcess(hit->GetCreatorProcess());
 
                aCgemMcHitCol->push_back(cgemMcHit);
 
                /*End Add*/
 
            } /* End of 'for(G4int i=0;i<n_hit;i++)' */     
        } /* End of 'if (n_hit > 0)' */
        } /* End of 'if (HCID > 0)' */
        //register CGEM McTruth collection to TDS
        StatusCode scCgem = m_evtSvc->registerObject("/Event/MC/CgemMcHitCol", aCgemMcHitCol);
        if(scCgem!=StatusCode::SUCCESS) 
        {
            G4cout << "ERROR : BesSim::BesMcTruthWriter::SaveCgemTruth(), Could not register CGEM McTruth collection!" << G4endl;
        }
 
        /*
        //retrieve CGEM McTruth from TDS
        SmartDataPtr<Event::CgemMcHitCol> aMcHitCol(m_evtSvc,"/Event/MC/CgemMcHitCol");
        if(!aMcHitCol)
        G4cout<<"Could not retrieve CGEM McTruth collection"<<G4endl;
 
        Event::CgemMcHitCol::iterator iMcHitCol;
        for(iMcHitCol=aMcHitCol->begin(); iMcHitCol!=aMcHitCol->end(); iMcHitCol++)
        {
        cout<<"Track_ID = "<<(*iMcHitCol)->GetTrackID()<<"\tPDG = "<<(*iMcHitCol)->GetPDGCode()<<"\tParent ID = "<<(*iMcHitCol)->GetParentID()<<endl;
        //       const Identifier ident = (*iMcHitCol)->identify();
        //       G4cout<<(*iMcHitCol)->getTrackIndex();
        //    G4cout<<" "<<CgemID::layer(ident);
        //    G4cout<<" "<<CgemID::wire(ident);
        //    G4cout<<" "<<(*iMcHitCol)->getDepositEnergy();
        //    G4cout<<" "<<(*iMcHitCol)->getDriftDistance();
        //    G4cout<<" "<<(*iMcHitCol)->getPositionX();
        //    G4cout<<" "<<(*iMcHitCol)->getPositionY();
        //    G4cout<<" "<<(*iMcHitCol)->getPositionZ();
        //    G4cout<<G4endl;
        }
        G4cout<<"end of retrieve CGEM McTruth collection"<<G4endl;
        */
}
 
void BesMcTruthWriter::SaveTofTruth()
{ 
  //Tof McHit collection defined in BOSS
  Event::TofMcHitCol*  aTofMcHitCol = new Event::TofMcHitCol;
 
  G4int HCID = -1;
  HCID = m_DigiMan->GetHitsCollectionID("BesTofHitsList");
  if(HCID>0)
  {
    BesTofHitsCollection* HC = 0;
    HC = (BesTofHitsCollection*) (m_DigiMan->GetHitsCollection(HCID));
    G4int n_hit = HC->entries();
    if(n_hit>0)
    {
      //arrange hits in hits collection in order of trackIndex
      BesTofHit* hit;
      vector<BesTofHit*>* vecHC = HC->GetVector();
      for(int i=0;i<n_hit-1;i++)
       for(int j=i+1;j<n_hit;j++)
         if((*vecHC)[i]->GetTrackIndex()>(*vecHC)[j]->GetTrackIndex())
         {
           hit = (*vecHC)[i];
           (*vecHC)[i] = (*vecHC)[j];
           (*vecHC)[j] = hit;
         }
      
      //push back TofMcHit collection to TDS
      for(G4int i=0;i<n_hit;i++)
      {
        hit = (*HC)[i];
        unsigned int scinNum, barrel_ec, layer = 0;
        scinNum = hit->GetScinNb();
        barrel_ec = hit->GetPartId();    
 
    //std::cout << "BesMcTruthWriter        PartID | scinNum  "  <<  barrel_ec  << "  |  "   << scinNum  << std::endl;
    //std::cout << "BesMcTruthWriter        is_barrel(barrel_ec)    "  <<  TofID::is_barrel(barrel_ec)  << std::endl;
 
    
        if (TofID::is_barrel(barrel_ec) && scinNum > TofID::getPHI_BARREL_MAX()) {
          layer = 1;
          scinNum = scinNum -  TofID::getPHI_BARREL_MAX() - 1;
        }
 
 
 
    Identifier ident;
    
    if(barrel_ec==0 || barrel_ec==1 || barrel_ec==2 )   //Old-Tof 
      {
        ident =  TofID::cell_id ( barrel_ec, layer, scinNum, 0);
      }
    else  //MRPC Part
      {
        int help = BesTofDigitizerEcV4::Calculate_Readoutstrip_number_continuum(hit->GetPos().x()*mm,hit->GetPos().y()*mm, barrel_ec, scinNum);
        scinNum = BesTofDigitizerEcV4::Produce_unique_identifier(scinNum, help);
        //std::cout << "BesMcTruthWriter        Readoutstrip "  <<  help  << std::endl;
        //std::cout << "BesMcTruthWriter        unique Identifier    "  <<  scinNum  << std::endl;
        ident =  TofID::cell_id_mrpc_mc ( barrel_ec, scinNum);
      }
 
   
    
        Event::TofMcHit* tofMcHit = new Event::TofMcHit( ident, hit->GetTrackIndex(),
            hit->GetPos().x(), hit->GetPos().y(), hit->GetPos().z(), hit->GetMomentum().x(),
            hit->GetMomentum().y(), hit->GetMomentum().z(), hit->GetTrackL(), hit->GetTime() );
 
        aTofMcHitCol->push_back(tofMcHit);
      }
    }
  }
 
  //register TOF McTruth collection to TDS
  StatusCode scTof = m_evtSvc->registerObject("/Event/MC/TofMcHitCol", aTofMcHitCol);
  if(scTof!=StatusCode::SUCCESS) 
    G4cout<< "Could not register TOF McTruth collection" <<G4endl;
 
  //retrieve TOF McTruth from TDS
  /*SmartDataPtr<Event::TofMcHitCol> aMcHitCol(m_evtSvc,"/Event/MC/TofMcHitCol");
  if(!aMcHitCol)
    G4cout<<"Could not retrieve TOF McTruth collection"<<G4endl;
 
  Event::TofMcHitCol::iterator iMcHitCol;
  for(iMcHitCol=aMcHitCol->begin(); iMcHitCol!=aMcHitCol->end(); iMcHitCol++)
  {
    const Identifier ident = (*iMcHitCol)->identify();
    G4cout<<(*iMcHitCol)->getTrackIndex();
    G4cout<<" "<<TofID::barrel_ec(ident);;
    G4cout<<" "<<TofID::layer(ident);
    G4cout<<" "<<TofID::phi_module(ident);
    G4cout<<" "<<(*iMcHitCol)->getPositionX();
    G4cout<<" "<<(*iMcHitCol)->getPositionY();
    G4cout<<" "<<(*iMcHitCol)->getPositionZ();
    G4cout<<" "<<(*iMcHitCol)->getPx();
    G4cout<<" "<<(*iMcHitCol)->getPy();
    G4cout<<" "<<(*iMcHitCol)->getPz();
    G4cout<<" "<<(*iMcHitCol)->getTrackLength();
    G4cout<<" "<<(*iMcHitCol)->getFlightTime();
    G4cout<<G4endl;
  }
  G4cout<<"end of retrieve TOF McTruth"<<G4endl;
  */
}  
 
void BesMcTruthWriter::SaveEmcTruth()
{
  //Emc McHit collection defined in BOSS
  Event::EmcMcHitCol*  aEmcMcHitCol = new Event::EmcMcHitCol;
 
  G4int fullMc = 1;
  if(fullMc==1) {     //full Mc Truth
    G4int HCID = -1;
    HCID = m_DigiMan->GetHitsCollectionID("BesEmcTruthHitsList");
    if(HCID>0)
    {
      BesEmcTruthHitsCollection* HC = 0;
      HC = (BesEmcTruthHitsCollection*) (m_DigiMan->GetHitsCollection(HCID));
      G4int n_hit = HC->entries();
      if(n_hit>0)
      {
        //arrange hits in hits collection in order of trackIndex
        BesEmcTruthHit* hit;
        vector<BesEmcTruthHit*>* vecHC = HC->GetVector();
        for(int i=0;i<n_hit-1;i++) {
          for(int j=i+1;j<n_hit;j++) {
            if((*vecHC)[i]->GetTrackIndex()>(*vecHC)[j]->GetTrackIndex()) {
              hit = (*vecHC)[i];
              (*vecHC)[i] = (*vecHC)[j];
              (*vecHC)[j] = hit;
            }
          }
        }
 
        for(G4int i=0;i<n_hit;i++)
        {
          BesEmcTruthHit* hit;
          hit = (*HC)[i];
 
          std::map<Identifier,double> hitMap;
          std::map<Identifier,double>::const_iterator iHitMap;
          hitMap.clear();
          for(iHitMap=hit->Begin();iHitMap!=hit->End();iHitMap++) {
            hitMap[iHitMap->first]=iHitMap->second;
          }
          
          Event::EmcMcHit* emcMcHit = new Event::EmcMcHit(hit->GetIdentify(), hit->GetTrackIndex(),
              hit->GetPosition().x(), hit->GetPosition().y(), hit->GetPosition().z(),
              hit->GetMomentum().x(), hit->GetMomentum().y(), hit->GetMomentum().z(),
              hit->GetEDep() );
 
          emcMcHit->setHitEmc(hit->GetHitEmc());
          emcMcHit->setPDGCode(hit->GetPDGCode());
          emcMcHit->setPDGCharge(hit->GetPDGCharge());
          emcMcHit->setTime(hit->GetTime());
          emcMcHit->setHitMap(hitMap);
 
          aEmcMcHitCol->push_back(emcMcHit);
        }
      }
    }
  } else {    //simple Mc Truth
    G4int HCID = -1;
    HCID = m_DigiMan->GetHitsCollectionID("BesEmcHitsList");
    if(HCID>0)
    {
      BesEmcHitsCollection* HC = 0;
      HC = (BesEmcHitsCollection*) (m_DigiMan->GetHitsCollection(HCID));
      G4int n_hit = HC->entries();
      if(n_hit>0)
      {
        //arrange hits in hits collection in order of trackIndex
        BesEmcHit* hit;
        vector<BesEmcHit*>* vecHC = HC->GetVector();
        for(int i=0;i<n_hit-1;i++)
          for(int j=i+1;j<n_hit;j++)
            if((*vecHC)[i]->GetTrackIndex()>(*vecHC)[j]->GetTrackIndex())
            {
              hit = (*vecHC)[i];
              (*vecHC)[i] = (*vecHC)[j];
              (*vecHC)[j] = hit;
            }
            
        //Emc McHit collection defined in BOSS
        Event::EmcMcHitCol*  aEmcMcHitCol = new Event::EmcMcHitCol;
        
        for(G4int i=0;i<n_hit;i++)
        {
          hit = (*HC)[i];
          Identifier ident =  EmcID::crystal_id ( hit->GetPartId(), hit->GetNumThetaCrystal(), hit->GetNumPhiCrystal() );
          
          std::map<Identifier,double> hitMap;
          Event::EmcMcHit* emcMcHit = new Event::EmcMcHit(ident, hit->GetTrackIndex(),
              hit->GetPosCrystal().x(), hit->GetPosCrystal().y(), hit->GetPosCrystal().z(),
              hit->GetMomentum().x(), hit->GetMomentum().y(), hit->GetMomentum().z(),
              hit->GetEdepCrystal() );
          
          aEmcMcHitCol->push_back(emcMcHit);
        }
      }
    }           
  }
 
  //register EMC McTruth collection to TDS
  StatusCode scEmc = m_evtSvc->registerObject("/Event/MC/EmcMcHitCol", aEmcMcHitCol);
  if(scEmc!=StatusCode::SUCCESS)
    G4cout<< "Could not register EMC McTruth collection" <<G4endl;
 
  //retrieve EMC McTruth from TDS
  /*SmartDataPtr<Event::EmcMcHitCol> aMcHitCol(m_evtSvc,"/Event/MC/EmcMcHitCol");
  if(!aMcHitCol)
    G4cout<<"Could not retrieve EMC McTruth collection"<<G4endl;
 
  Event::EmcMcHitCol::iterator iMcHitCol;
  for(iMcHitCol=aMcHitCol->begin(); iMcHitCol!=aMcHitCol->end(); iMcHitCol++)
  {
    const Identifier ident = (*iMcHitCol)->identify();
    G4cout<<(*iMcHitCol)->getTrackIndex();
    G4cout<<" "<<EmcID::barrel_ec(ident);
    G4cout<<" "<<EmcID::theta_module(ident);
    G4cout<<" "<<EmcID::phi_module(ident);
    G4cout<<" "<<(*iMcHitCol)->getPositionX();
    G4cout<<" "<<(*iMcHitCol)->getPositionY();
    G4cout<<" "<<(*iMcHitCol)->getPositionZ();
    G4cout<<" "<<(*iMcHitCol)->getPx();
    G4cout<<" "<<(*iMcHitCol)->getPy();
    G4cout<<" "<<(*iMcHitCol)->getPz();
    G4cout<<" "<<(*iMcHitCol)->getDepositEnergy();
    G4cout<<G4endl;
  }
  G4cout<<"end of retrieve EMC McTruth"<<G4endl;
  */
}
 
void BesMcTruthWriter::SaveMucTruth()
{
  //Muc McHit collection defined in BOSS
  Event::MucMcHitCol*  aMucMcHitCol = new Event::MucMcHitCol;
 
  G4int HCID = -1;
  HCID = m_DigiMan->GetHitsCollectionID("BesMucHitsList");
  if(HCID>0)
  {
    BesMucHitsCollection* HC = 0;
    HC = (BesMucHitsCollection*) (m_DigiMan->GetHitsCollection(HCID));
    G4int n_hit = HC->entries();
    if(n_hit>0)
    {
      //arrange hits in hits collection in order of trackIndex
      BesMucHit* hit;
      vector<BesMucHit*>* vecHC = HC->GetVector();
      for(int i=0;i<n_hit-1;i++)
        for(int j=i+1;j<n_hit;j++)
          if((*vecHC)[i]->GetTrackIndex()>(*vecHC)[j]->GetTrackIndex())
          {
            hit = (*vecHC)[i];
            (*vecHC)[i] = (*vecHC)[j];
            (*vecHC)[j] = hit;
          }
      
      for(G4int i=0;i<n_hit;i++)
      {
        hit = (*HC)[i];
        Identifier ident =  MucID::channel_id ( hit->GetPart(), hit->GetSeg(), hit->GetGap(),
          hit->GetStrip() );
        Event::MucMcHit* mucMcHit = new Event::MucMcHit(ident, hit->GetTrackIndex(), 
            hit->GetPos().x(), hit->GetPos().y(), hit->GetPos().z(), hit->GetMomentum().x(),
            hit->GetMomentum().y(), hit->GetMomentum().z() );
        aMucMcHitCol->push_back(mucMcHit);
      }
      
    }
  }
 
  //register MUC McTruth collection to TDS
  StatusCode scMuc = m_evtSvc->registerObject("/Event/MC/MucMcHitCol", aMucMcHitCol);
  if(scMuc!=StatusCode::SUCCESS) 
    G4cout<< "Could not register MUC McTruth collection" <<G4endl;
 
  //retrieve MUC McTruth from TDS
  /*SmartDataPtr<Event::MucMcHitCol> aMcHitCol(m_evtSvc,"/Event/MC/MucMcHitCol");
  if(!aMcHitCol)
    G4cout<<"Could not retrieve MUC McTruth collection"<<G4endl;
 
  Event::MucMcHitCol::iterator iMcHitCol;
  for(iMcHitCol=aMcHitCol->begin(); iMcHitCol!=aMcHitCol->end(); iMcHitCol++)
  {
    const Identifier ident = (*iMcHitCol)->identify();
    G4cout<<(*iMcHitCol)->getTrackIndex();
    G4cout<<" "<<MucID::part(ident);
    G4cout<<" "<<MucID::seg(ident);
    G4cout<<" "<<MucID::gap(ident);
    G4cout<<" "<<MucID::strip(ident);
    G4cout<<" "<<(*iMcHitCol)->getPositionX();
    G4cout<<" "<<(*iMcHitCol)->getPositionY();
    G4cout<<" "<<(*iMcHitCol)->getPositionZ();
    G4cout<<" "<<(*iMcHitCol)->getPx();
    G4cout<<" "<<(*iMcHitCol)->getPy();
    G4cout<<" "<<(*iMcHitCol)->getPz();
    G4cout<<G4endl;
  }
  G4cout<<"end of retrieve MUC McTruth"<<G4endl;
  */
}