BesRawDataWriter.cc

Go to the documentation of this file.

//---------------------------------------------------------------------------//
////      BOOST --- BESIII Object_Oriented Simulation Tool                     //
////---------------------------------------------------------------------------//
////Description:
////Author:   Dengzy
////Created:  Mar, 2004
////Modified:
////Comment:
//
#include "BesMdcDigi.hh"
#include "BesCgemDigi.hh"
#include "BesTofDigi.hh"
#include "BesEmcDigi.hh"
#include "BesMucDigi.hh"
#include "BesRawDataWriter.hh"
#include "G4DigiManager.hh"
#include "G4Svc/IG4Svc.h"
#include "G4Svc/G4Svc.h"
 
#include "Randomize.hh"
 
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/RegistryEntry.h"
#include "GaudiKernel/MsgStream.h"
 
#include "MdcRawEvent/MdcDigi.h"
#include "CgemRawEvent/CgemDigi.h"
#include "TofRawEvent/TofDigi.h"
#include "EmcRawEvent/EmcDigi.h"
#include "MucRawEvent/MucDigi.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 "RawEvent/RawDataUtil.h"
#include "RawEvent/DigiEvent.h"
#include "GaudiKernel/SmartDataPtr.h"
 
 
BesRawDataWriter::BesRawDataWriter()
{
  m_DigiMan = G4DigiManager::GetDMpointer(); 
  //yuany 2005-2-28
  mdcGeoPointer=BesMdcGeoParameter::GetGeo();
 
  ICgemGeomSvc* ISvc;
  StatusCode sc=Gaudi::svcLocator()->service("CgemGeomSvc", ISvc);
  if (!sc.isSuccess())
    G4cout <<" BesRawDataWriter::Error,could not open CgemGeomSvc"<<G4endl;
  cgemGeoPointer=dynamic_cast<CgemGeomSvc *>(ISvc);
 
  IG4Svc* tmpSvc;
  sc=Gaudi::svcLocator()->service("G4Svc", tmpSvc);
  if (!sc.isSuccess())
    G4cout <<" BesRawDataWriter::Error,could not open G4Svc"<<G4endl;
  m_G4Svc=dynamic_cast<G4Svc *>(tmpSvc); 
 
  IRealizationSvc *tmpReal;
  StatusCode scReal = Gaudi::svcLocator()->service("RealizationSvc",tmpReal);
  if (!scReal.isSuccess())
  {
    std::cout << " Could not initialize Realization Service in BesRawDataWriter" << std::endl;
  } 
  else 
  {
    m_RealizationSvc=dynamic_cast<RealizationSvc*>(tmpReal);
  }
 
}
 
BesRawDataWriter::~BesRawDataWriter()
{
}
 
void BesRawDataWriter::SaveDigits()
{
  //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;
 
  DigiEvent* aDigiEvent = new DigiEvent;
  sc = m_evtSvc->registerObject("/Event/Digi",aDigiEvent);
  if(sc!=StatusCode::SUCCESS) {
    G4cout<< "Could not register DigiEvent" <<G4endl;
  }
 
  SaveMdcDigits();
  SaveCgemDigits();
  SaveTofDigits();
  SaveEmcDigits();
  SaveMucDigits();
}
 
void BesRawDataWriter::SaveMdcDigits()
{
  //mdc digis collection defined in BOSS
  MdcDigiCol* aMdcDigiCol = new MdcDigiCol;
 
  G4int mdcDigiCollID = -1;
  mdcDigiCollID = m_DigiMan->GetDigiCollectionID("BesMdcDigisCollection");
  if(mdcDigiCollID>=0)
  {  
    BesMdcDigisCollection* mdcDC = (BesMdcDigisCollection*)m_DigiMan->GetDigiCollection(mdcDigiCollID);
    G4int nDigi = mdcDC->entries();
    if(nDigi>0)
    {
      //push back mdc digits to MdcDigiCol in BOSS
      BesMdcDigi* mdcDigi;
      for(int i=0;i<nDigi;i++)
      {
        mdcDigi=(*mdcDC)[i];
        unsigned int charge;
        if(m_G4Svc->GetMdcDedxFlag()==1){
          charge = int(mdcDigi->GetEdep());
        }else{
          charge = RawDataUtil::MdcChargeChannel(mdcDigi->GetEdep());
        }
        unsigned int time = RawDataUtil::MdcTimeChannel(mdcDigi->GetDriftT());
        const Identifier ident =  MdcID::wire_id ( mdcDigi->GetLayerNo(), mdcDigi->GetCellNo() );
 
        MdcDigi* aMdcDigi = new MdcDigi(ident, time, charge);  
        aMdcDigi->setTrackIndex(mdcDigi->GetTrackID());
        aMdcDigiCol->push_back(aMdcDigi);      
      }
    }
  }
 
  //register MDC digits collection to TDS
  StatusCode scMdc = m_evtSvc->registerObject("/Event/Digi/MdcDigiCol", aMdcDigiCol);
  if(scMdc!=StatusCode::SUCCESS) 
    G4cout<< "Could not register MDC digi collection" <<G4endl;
 
  //retrieve MDC digits from TDS
  /*SmartDataPtr<MdcDigiCol> aDigiCol(m_evtSvc,"/Event/Digi/MdcDigiCol");
    if(!aDigiCol)
    G4cout<<"Could not retrieve MDC digi collection"<<G4endl;
 
    MdcDigiCol::iterator iDigiCol;
    for(iDigiCol=aDigiCol->begin(); iDigiCol!=aDigiCol->end(); iDigiCol++)
    {
    const Identifier ident = (*iDigiCol)->identify();
    G4cout<<"layer: "<<MdcID::layer(ident);
    G4cout<<"  cell: "<<MdcID::wire(ident);
    G4cout<<"  charge: "<<(*iDigiCol)->getChargeChannel();
    G4cout<<"  time: "<<(*iDigiCol)->getTimeChannel()<<G4endl;
    }
    G4cout<<"end of retrieve MDC digi collection"<<G4endl;
   */
}
 
void BesRawDataWriter::SaveCgemDigits()
{
     bool printFlag=false;
     if(printFlag) {
      G4cout << "|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||" << G4endl;
      G4cout << "INFO : BesSim::BesRawDataWriter::SaveCgemDigits(), Begin!" << G4endl;
     }
 
     /* cgem digis collection defined in BOSS */
     CgemDigiCol* aCgemDigiCol = new CgemDigiCol;
 
     G4int cgemDigiCollID = -1;
     cgemDigiCollID = m_DigiMan->GetDigiCollectionID("BesCgemDigisCollection");
     if (cgemDigiCollID >= 0)
     {  
      BesCgemDigisCollection *cgemDC=(BesCgemDigisCollection*)m_DigiMan->GetDigiCollection(cgemDigiCollID);
      G4int nDigi = cgemDC->entries();
      if (nDigi > 0)
      {
           /* push back cgem digits to CgemDigiCol in BOSS */
           BesCgemDigi* cgemDigi;
           for(int i=0;i<nDigi;i++)
           {
            cgemDigi=(*cgemDC)[i];
            unsigned int charge;
 
            double stripEff = m_G4Svc->GetCgemStripEff();
            G4double frandom = G4UniformRand();
            // cout << setw(10) << cgemDigi->GetLayerID() << setw(10) << cgemDigi->GetStripID();
            // cout << "stripEff = " << stripEff << "  ,  frandom = " << frandom << endl;
            if(frandom > stripEff) continue;
 
            if(m_G4Svc->GetCgemDedxFlag()==1)
            {
             //charge = int(cgemDigi->GetEnergyDeposit());
             charge = int(cgemDigi->GetEnergyDeposit()*1e10); /* Change MeV to eV */
             //cout<<"use dE/dx flag!"<<endl;
            }
            else
            {
             charge = RawDataUtil::CgemChargeChannel(cgemDigi->GetEnergyDeposit());
            }
 
            unsigned int time = RawDataUtil::CgemTimeChannel(cgemDigi->GetGlobalTime());
            const Identifier ident =  CgemID::strip_id(cgemDigi->GetLayerID   (),
                                   cgemDigi->GetSheetID   (),
                                   cgemDigi->GetStripType (),
                                   cgemDigi->GetStripID   () );
            double Q_in_fC=cgemDigi->GetEnergyDeposit();
            double T_in_ns=cgemDigi->GetGlobalTime();
            if(printFlag) 
             G4cout  << "digiID="     << ident.get_value()
                 << "   layerID=" << CgemID::layer(ident)
                 << "   sheetID=" << CgemID::sheet(ident)
                 << "   isXstrip=" << CgemID::is_xstrip(ident)
                 << "   stripID=" << CgemID::strip(ident)
                 << "   time="    << time                 
                 << "   charge="  << charge                    //<< G4endl;
                 << "   time_ns="    << T_in_ns                 
                 << "   charge_fc="  << Q_in_fC                   << G4endl;
 
            CgemDigi* aCgemDigi = new CgemDigi(ident, time, charge);  
            aCgemDigi->setTrackIndex(cgemDigi->GetTrackID());
            
            aCgemDigi->setCharge_fc(Q_in_fC);
            aCgemDigi->setTime_ns(T_in_ns);
 
            aCgemDigiCol->push_back(aCgemDigi);      
 
           } /* End of 'for(int i=0;i<nDigi;i++)' */
      } /* End of 'if (nDigi > 0)' */
     } /* End of 'if (cgemDigiCollID >= 0)' */
 
     /* register CGEM digits collection to TDS */
     StatusCode scCgem = m_evtSvc->registerObject("/Event/Digi/CgemDigiCol", aCgemDigiCol);
     if(scCgem!=StatusCode::SUCCESS) 
     {
      G4cout << "ERROR : BesSim::BesRawDataWriter::SaveCgemDigits(), Could not register CGEM digi collection! " << G4endl;
     }
 
 
     //std::cout<<__FILE__<<" "<<__LINE__<<std::endl;
     //CgemDigiCol::iterator iter = aCgemDigiCol->begin();
     //for ( ; iter != aCgemDigiCol->end(); iter++) 
     //{
     //  std::cout << " layer "  << CgemID::layer((*iter)->identify())
     //            << " sheet "  << CgemID::sheet((*iter)->identify())
     //            << " time "   << (*iter)->getTimeChannel()
     //            << " charge " << (*iter)->getChargeChannel()        << std::endl;
     //}
 
     //retrieve CGEM digits from TDS
     /*SmartDataPtr<CgemDigiCol> aDigiCol(m_evtSvc,"/Event/Digi/CgemDigiCol");
       if(!aDigiCol)
       G4cout<<"Could not retrieve CGEM digi collection"<<G4endl;
 
       CgemDigiCol::iterator iDigiCol;
       for(iDigiCol=aDigiCol->begin(); iDigiCol!=aDigiCol->end(); iDigiCol++)
       {
       const Identifier ident = (*iDigiCol)->identify();
       G4cout<<"layer: "<<CgemID::layer(ident);
       G4cout<<"  cell: "<<CgemID::wire(ident);
       G4cout<<"  charge: "<<(*iDigiCol)->getChargeChannel();
       G4cout<<"  time: "<<(*iDigiCol)->getTimeChannel()<<G4endl;
       }
       G4cout<<"end of retrieve CGEM digi collection"<<G4endl;
     */
}
 
void BesRawDataWriter::SaveTofDigits()
{  
  //Tof digits collection defined in BOSS
  TofDigiCol* aTofDigiCol = new TofDigiCol;
 
  G4int tofDigiCollID = -1;
  tofDigiCollID = m_DigiMan->GetDigiCollectionID("BesTofDigitsCollection");
 
  if(tofDigiCollID>=0)
  {
    BesTofDigitsCollection* tofDC = (BesTofDigitsCollection*)m_DigiMan->GetDigiCollection(tofDigiCollID);
    G4int nDigi = tofDC->entries();
    if(nDigi>0)
    {
      //arrange digis in digitsCollection in order of trackIndex
 
      //std::cout << "BesRawDataWriter      Arange Digis     " << std::endl;
      BesTofDigi* digi;
      vector<BesTofDigi*>* vecDC = tofDC->GetVector();
      for(int i=0;i<nDigi-1;i++)
        for(int j=i+1;j<nDigi;j++)
          if((*vecDC)[i]->GetTrackIndex()>(*vecDC)[j]->GetTrackIndex())
          {
            digi = (*vecDC)[i];
            (*vecDC)[i] = (*vecDC)[j];
            (*vecDC)[j] = digi;
          }
      //std::cout << "BesRawDataWriter      Digis are arranged     " << std::endl;
 
 
      //push back tof digits to TofDigiCol in BOSS
      for(int i=0;i<nDigi;i++)
      {
    //std::cout << "BesRawDataWriter      Try:   digi = (*tofDC)[i]     " << std::endl;
        digi = (*tofDC)[i];
    //std::cout << "BesRawDataWriter      OK:   digi = (*tofDC)[i]     " << std::endl;
        G4int scinNum = digi->GetScinNb();
        G4int barrel_ec = digi->GetPartId();
 
    
    //std::cout << "BesRawDataWriter      scinNum:     "<< scinNum << std::endl;
    //std::cout << "BesRawDataWriter      barrel_ec:     "<< barrel_ec << std::endl;
    //std::cout << "BesRawDataWriter       is_barrel(barrel_ec):     "<<TofID::is_barrel(barrel_ec)<<std::endl;
 
 
    //#Matthias. scinNum:     154
    //#Matthias. barrel_ec:     1
    //#Matthias is_barrel(barrel_ec):     1
 
 
    unsigned int layer = 0;
    Identifier ident;
    unsigned int time =0;
 
 
 
 
    if(barrel_ec == 0 || barrel_ec == 1 || barrel_ec == 2)  //These numbers stand for the old TOF
    {
      //std::cout << "BesRawDataWriter:    You are using the old TOF as TOF Detector. PartID = " << barrel_ec << std::endl;
      
 
      if (TofID::is_barrel(barrel_ec) && scinNum > TofID::getPHI_BARREL_MAX()) {
        layer = 1;
        scinNum = scinNum -  TofID::getPHI_BARREL_MAX() - 1;
      }
 
      // Store tofDigi for end==0
      ident =  TofID::cell_id ( barrel_ec, layer, scinNum, 0); 
          
      if(digi->GetForwTDC() > 0 )
          time = RawDataUtil::TofTimeChannel(digi->GetForwTDC());
      else
          time = 0;
      
 
        unsigned int charge;
        
        if(digi->GetForwADC() > 0 )
        {
          // ADC linear conversion
          charge = digi->GetForwADC(); // convert ADC(pC) of bTof into Q channel in Data
               // Treat ADC overflow
          if(charge>=8192)
          {
            // The overlow flag is the 20th bit of the ADC_channel.
            charge = (charge | 0x080000);// if ADC_channel>=8192, set overflow flag = 1.  
        //  10000000000000000000  //The 20th bit is the overflow flag
        //         1111111111111  // =8191   -> 13 bits are required to store the information
        //   1111110000000000000  //The time information is stored in bit 14-19, however the time is stored in an additonal indentifier, too.
        //std::cout << "BesRawDataWriter        Charge overflow ForwardADC" << std::endl;
          }
          charge =  ( (time & 0x07e000) | charge) ; // Store tclock into ADC_channel as qclock
        }
        else
          charge =0;
 
    
        TofDigi* tofDigi = new TofDigi(ident, time, charge);
 
        if( charge & 0x80000 ) // Treat overflow flag(20th bit)
    {tofDigi->setOverflow(0x2);}
        
    tofDigi->setTrackIndex(digi->GetTrackIndex());
        aTofDigiCol->push_back(tofDigi);
 
    
    if(barrel_ec == 0 || barrel_ec == 1 || barrel_ec == 2)  //These numbers stand for the old TOF    
      {
 
    ident =  TofID::cell_id ( barrel_ec, layer, scinNum, 1); 
        if(digi->GetBackTDC()>0)
          time = RawDataUtil::TofTimeChannel(digi->GetBackTDC());
        else
          time = 0;
    /*
    G4cout << G4endl << "#Matthias charge:  " << charge << G4endl;
    G4cout <<  "#Matthias time  " << time << G4endl;
    G4cout <<  "#Matthias time /fac  " << time*0.00239 << G4endl;
    G4cout <<  "#Matthias (time & 0x07e000) :  " << (time & 0x07e000)  << G4endl;
    G4cout <<  "#Matthias ( (time & 0x07e000) | charge)  " << ( (time & 0x07e000) | charge)  << G4endl;
    if( charge & 0x80000 ) G4cout <<  "#Matthias: overflow" << G4endl; 
    if( ( (time & 0x07e000) | charge) != charge )  G4cout <<  "#Matthias: ungleich" << G4endl;
    */
        if(digi->GetBackADC()>0 )
        {
 
          // ADC linear conversion
          charge = digi->GetBackADC();
      
          // Treat ADC overflow
          if(charge>=8192)
          {
            // The overlow flag is the 20th bit of the ADC_channel.
            charge = (charge | 0x080000);// if ADC_channel>=8192, set overflow flag = 1.
          }
          charge =  ( (time & 0x07e000) | charge) ; // Store tclock into ADC_channel as qclock
        }
        else
          charge = 0;
 
 
        
        tofDigi = new TofDigi(ident, time, charge);
 
        if( charge & 0x80000 ) // Treat overflow flag(20th bit)
          tofDigi->setOverflow(0x2);
        tofDigi->setTrackIndex(digi->GetTrackIndex());
    
        aTofDigiCol->push_back(tofDigi);
 
 
      }//close if(barrel_ec == 0 || ...)
 
  }//close if(barrel_ec == 0 || barrel_ec == 1 || barrel_ec == 2) 
 
 
 
  ////////////////
  /// NewMRPC ////
  ////////////////
 
  else if(barrel_ec == 3 || barrel_ec == 4 || barrel_ec == 5 || barrel_ec == 6) //MRPC-Detector
    {
      //std::cout << "BesRawDataWriter:    You are using the new MRPC as TOF Detector. PartID | scinNum   " << barrel_ec << "   |   " << scinNum << std::endl;
      ident =  TofID::cell_id_mrpc ( barrel_ec, scinNum);
  
 
      //static double TofTime(unsigned int timeChannel) { return ((double)timeChannel)*TOF_TIME_FACTOR; }
      //static unsigned int TofTimeChannel (double time) { return (unsigned int)(time/TOF_TIME_FACTOR); } //Konvertiere TofTime in Zeitchannel: factor = 0.00293
 
          if(digi->GetForwTDC() > 0 )
        time = RawDataUtil::TofTimeChannel(digi->GetForwTDC());
          else
        time = 0;
 
 
 
      unsigned int time2;
 
      if(digi->GetBackTDC() > 0 ) time2 = RawDataUtil::TofTimeChannel(digi->GetBackTDC());
      else time2 = 0;
 
 
      
      /*
      std::cout << "Time1  ns  " << digi->GetForwTDC() << std::endl;
      std::cout << "Time2  ns  " << digi->GetBackTDC() << std::endl;
      std::cout << "Time1  channel  " << time << std::endl;
      std::cout << "Time2  channel  " << time2 << std::endl;
      std::cout << "charge from time2    " << exp(RawDataUtil::TofTime(time2-time)-10.)/1000. << " pC"<< std::endl;
      */
 
 
 
    //std::cout << "BesRawDataWriter       : ident     "<<ident<<std::endl;
    //std::cout << "BesRawDataWriter       : partId    |time     "<<  barrel_ec <<"     "<<time<<std::endl;
 
 
 
 
        /*  //Here I did store the charge directly, this is not done anymore as the real raw data will be to times. The difference of those time decodes the input charge!
        unsigned int charge;
 
        if(digi->GetForwADC() > 0 )
        {
          // ADC linear conversion
          charge = digi->GetForwADC(); // convert ADC(pC) of bTof into Q channel in Data
      //std::cout << "BesRawDataWriter      :partID    |  charge   "<<  barrel_ec<<"     "<<charge<<std::endl;
 
          // Treat ADC overflow
          if(charge>=8192)
          {
            // The overlow flag is the 20th bit of the ADC_channel.
            charge = (charge | 0x080000);// if ADC_channel>=8192, set overflow flag = 1.  
        //  10000000000000000000  //The 20th bit is the overflow flag
        //         1111111111111  // =8191   -> 13 bits are required to store the information
        //   1111110000000000000  //The time information is stored in bit 14-19, however the time is stored in an additonal indentifier, too.
        //std::cout << "BesRawDataWriter        Charge overflow ForwardADC" << std::endl;
          }
          charge =  ( (time & 0x07e000) | charge) ; // Store tclock into ADC_channel as qclock
        }
        else
          charge =0;
 
 
    //std::cout << "BesRawDataWriter      : Ende - ident    "<<ident<<std::endl;
    
    //  std::cout << "BesRawDataWriter      : real - time     "<<digi->GetForwTDC()<<std::endl;
    //std::cout << "BesRawDataWriter      :  charge mit Zeit   " <<charge<<std::endl;
    //std::cout << "BesRawDataWriter      :  charge rekons     " <<(charge & 0x1fff)<<std::endl;
    //std::cout << "BesRawDataWriter      : TDC - time     "<<RawDataUtil::TofTimeChannel(digi->GetForwTDC())<<std::endl;
    //std::cout << "BesRawDataWriter      : ADC - charge   "<<digi->GetForwADC()<<std::endl;
    */
 
 
        TofDigi* tofDigi = new TofDigi(ident, time, time2);
    //std::cout << "BesRawDataWriter      : TofDigi is stored   "<<std::endl;
            
    tofDigi->setTrackIndex(digi->GetTrackIndex());
    //std::cout << "BesRawDataWriter      : TofDigi: SetTrackIndex   "<<std::endl;
        aTofDigiCol->push_back(tofDigi);
    //std::cout << "BesRawDataWriter      : TofDigi: Pushback   "<<std::endl;
 
 
 
    }
 
 
 
 
 
 
 
    
 
 
 
      }//close for
    } //close if(nDigi>0)
  } //close if(tofDigiCollID>=0)
 
  //register TOF digits collection to TDS
  StatusCode scTof = m_evtSvc->registerObject("/Event/Digi/TofDigiCol", aTofDigiCol);
  if(scTof!=StatusCode::SUCCESS) 
    G4cout<< "Could not register TOF digi collection" <<G4endl;
 
  //std::cout <<  "BesRawDataWriter      : TOF digi colleceton is registered!  "<<std::endl; 
 
 
  //retrieve TOF digits from TDS
  /*SmartDataPtr<TofDigiCol> aDigiCol(m_evtSvc,"/Event/Digi/TofDigiCol");
    if(!aDigiCol)
    G4cout<<"Could not retrieve TOF digi collection"<<G4endl;
 
    TofDigiCol::iterator iDigiCol;
    for(iDigiCol=aDigiCol->begin(); iDigiCol!=aDigiCol->end(); iDigiCol++)
    {
    const Identifier ident = (*iDigiCol)->identify();
    G4cout<<"partId: "<<TofID::barrel_ec(ident);
    G4cout<<"  layer: "<<TofID::layer(ident);
    G4cout<<"  scinNb: "<<TofID::phi_module(ident);
  // G4cout<<"  end: "<<TofID::end(ident);
  G4cout<<"  charge: "<<(*iDigiCol)->getChargeChannel();
  G4cout<<"  time: "<<(*iDigiCol)->getTimeChannel()<<G4endl;
  }
  G4cout<<"end of retrieve TOF digits"<<G4endl;
   */
}  
 
void BesRawDataWriter::SaveEmcDigits()
{
  //Emc digits collection defined in BOSS
  EmcDigiCol* aEmcDigiCol = new EmcDigiCol;
 
  G4int emcDigiCollID = -1;
  emcDigiCollID = m_DigiMan->GetDigiCollectionID("BesEmcDigitsCollection");
  if(emcDigiCollID>=0)
  {
    BesEmcDigitsCollection* emcDC = (BesEmcDigitsCollection*)m_DigiMan->GetDigiCollection(emcDigiCollID);
    G4int nDigi = emcDC->entries();
    if(nDigi>0)
    {
      //arrange digis in digitsCollection in order of trackIndex
      BesEmcDigi* digi;
      vector<BesEmcDigi*>* vecDC = emcDC->GetVector();
      for(int i=0;i<nDigi-1;i++)
        for(int j=i+1;j<nDigi;j++)
          if((*vecDC)[i]->GetTrackIndex()>(*vecDC)[j]->GetTrackIndex())
          {
            digi = (*vecDC)[i];
            (*vecDC)[i] = (*vecDC)[j];
            (*vecDC)[j] = digi;
          }
      //push back EmcDigi to EmcDigiCol in BOSS
      for(int i=0;i<nDigi;i++)
      {
        digi = (*emcDC)[i];
 
        int charge = RawDataUtil::EmcChargeChannel(digi->GetEnergy());
        int measure = RawDataUtil::EmcChargeMeasure(digi->GetEnergy());
        int time = RawDataUtil::EmcTimeChannel(digi->GetTime());
        Identifier ident =  EmcID::crystal_id (digi->GetPartId(), digi->GetThetaNb(),  digi->GetPhiNb());
        EmcDigi* emcDigi = new EmcDigi(ident, time, charge);
        emcDigi->setTrackIndex(digi->GetTrackIndex());
        emcDigi->setMeasure(measure);
        aEmcDigiCol->push_back(emcDigi);
      }
    }
  }
 
  //register EMC digits collection to TDS
  StatusCode scEmc = m_evtSvc->registerObject("/Event/Digi/EmcDigiCol", aEmcDigiCol);
  if(scEmc!=StatusCode::SUCCESS)
    G4cout<< "Could not register EMC digi collection" <<G4endl;
 
  //retrieve EMC digits from TDS
  /*SmartDataPtr<EmcDigiCol> aDigiCol(m_evtSvc,"/Event/Digi/EmcDigiCol");
    if(!aDigiCol)
    G4cout<<"Could not retrieve EMC digi collection"<<G4endl;
 
    EmcDigiCol::iterator iDigiCol;
    for(iDigiCol=aDigiCol->begin(); iDigiCol!=aDigiCol->end(); iDigiCol++)
    {
    const Identifier ident = (*iDigiCol)->identify();
    G4cout<<"barrel_ec: "<<EmcID::barrel_ec(ident);
    G4cout<<"  theta: "<<EmcID::theta_module(ident);
    G4cout<<"  phi: "<<EmcID::phi_module(ident);
    G4cout<<"  charge: "<<(*iDigiCol)->getChargeChannel();
    G4cout<<"  time: "<<(*iDigiCol)->getTimeChannel()<<G4endl;
    }
    G4cout<<"end of retrieve EMC digits"<<G4endl;
   */
}
 
void BesRawDataWriter::SaveMucDigits()
{
  //Muc digits collection defined in BOSS
  MucDigiCol* aMucDigiCol = new MucDigiCol;
 
  G4int mucDigiCollID =-1;
  mucDigiCollID = m_DigiMan->GetDigiCollectionID("BesMucDigisCollection");
  if(mucDigiCollID>=0)
  {
    BesMucDigisCollection* mucDC = (BesMucDigisCollection*)m_DigiMan->GetDigiCollection(mucDigiCollID);
    G4int nDigi = mucDC->entries();
 
    if(nDigi > 0) {
      BesMucDigi* digi;
      for(int i = 0; i < nDigi; i++) {
        digi = (*mucDC)[i];
        Identifier ident =  MucID::channel_id( digi->GetPartId(), digi->GetSegId(),digi->GetGapId(),digi->GetStripId() );
        MucDigi* mucDigi = new MucDigi(ident);
        mucDigi->setTrackIndex(digi->GetTrackIndex());
        aMucDigiCol->push_back(mucDigi);
      }
    }
  }
 
  //register MUC digits collection to TDS
  StatusCode scMuc = m_evtSvc->registerObject("/Event/Digi/MucDigiCol", aMucDigiCol);
  if(scMuc!=StatusCode::SUCCESS)
    G4cout<< "Could not register MUC digi collection" <<G4endl;
 
  //retrieve MUC digits from TDS
  /*SmartDataPtr<MucDigiCol> aDigiCol(m_evtSvc,"/Event/Digi/MucDigiCol");
    if(!aDigiCol)
    G4cout<<"Could not retrieve MUC digi collection"<<G4endl;
 
    MucDigiCol::iterator iDigiCol;
    for(iDigiCol=aDigiCol->begin(); iDigiCol!=aDigiCol->end(); iDigiCol++)
    {
    const Identifier ident = (*iDigiCol)->identify();
    G4cout<<"Part: "<<MucID::part(ident);
    G4cout<<"  Seg: "<<MucID::seg(ident);
    G4cout<<"  Gap: "<<MucID::gap(ident);
    G4cout<<"  Strip: "<<MucID::strip(ident)<<G4endl;
    }
    G4cout<<"end of retrieve MUC digits"<<G4endl;
   */
}