MrpcDBSRec.cxx

Go to the documentation of this file.

//  Package: MrpcDBSRec
//  BESIII Tof Reconstruction Algorithm. 
//  Created by Sun Shengsen (EPC IHEP)
//
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/AlgFactory.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/SmartDataLocator.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/PropertyMgr.h"
#include "EventModel/EventHeader.h"
#include "EventModel/EventModel.h"
#include "ReconEvent/ReconEvent.h"
#include "McTruth/McParticle.h"
#include "McTruth/TofMcHit.h"
#include "McTruth/McParticle.h"
#include "EvTimeEvent/RecEsTime.h"
#include "ExtEvent/RecExtTrack.h"
#include "DstEvent/TofHitStatus.h"
#include "TofRecEvent/RecTofTrack.h"
 
//#include "TofRecEvent/RecBTofCalHit.h"
//#include "TofRecEvent/RecETofCalHit.h"
 
#include "TofGeomSvc/ITofGeomSvc.h"
#include "TofCaliSvc/ITofCaliSvc.h"
#include "RawDataProviderSvc/IRawDataProviderSvc.h"
#include "RawDataProviderSvc/TofData.h"
#include "MrpcRecDBS/MrpcDBSRec.h"
#include "MrpcRecDBS/MrpcDBSTrack.h"
#include "MrpcRecDBS/MrpcDBSCount.h"
#include "MrpcRecDBS/MrpcDBSRecTDS.h"
#include <iostream>
#include "ReconEvent/ReconEvent.h"
#include "McTruth/McParticle.h"
#include "McTruth/TofMcHit.h"
#include "McTruth/McParticle.h"
#include "MdcRecEvent/RecMdcTrack.h"
#include "EmcRecEventModel/RecEmcShower.h"
 
#include "TofSim/BesTofDigitizerEcV4_dbs.hh"
#include "Identifier/TofID.h"
 
 
using namespace std;
using namespace Event;
 
/////////////////////////////////////////////////////////////////////////////
 
ITofGeomSvc* tofGeomSvc;
ITofCaliSvc* tofCaliSvc;
// ITofQCorrSvc* tofQCorrSvc;
// ITofQElecSvc* tofQElecSvc;
IRawDataProviderSvc* tofDigiSvc;
 
MrpcDBSRec::MrpcDBSRec(const std::string& name, ISvcLocator* pSvcLocator) :
  Algorithm(name, pSvcLocator)
{   
  declareProperty( "AcceleratorStatus", m_acceleratorStatus );
  declareProperty( "MagneticField",     m_magneticField     );
  declareProperty( "FirstIteration",    m_firstIteration    );
  declareProperty( "PrintOutInfo",      m_printOutInfo      );
  declareProperty( "neighborhood",      m_neighborhood=4    );//For the time being: 5 is the standardvalue -> It yield the highest reconstruction efficiency! 
  declareProperty( "store_trans_time",  m_store_trans_time=0);//Usually the root file is not produced!
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
 
StatusCode MrpcDBSRec::initialize(){
 
 
 
 
 
  MsgStream log(msgSvc(), name());
  log << MSG::INFO << "MrpcRecDBS in initialize()" << endreq;   
 
  //Get TOF Geomertry Service
  StatusCode scg = service("TofGeomSvc", tofGeomSvc);
  if (scg==  StatusCode::SUCCESS) {
    log << MSG::INFO << "MrpcRecDBS Get Geometry Service Sucessfully!" << endreq;
  }else {
    log << MSG::ERROR << "MrpcRecDBS Get Geometry Service Failed !" << endreq;
    return StatusCode::FAILURE;
  }
 
  //Get TOF Calibtration Service
  StatusCode scc = service("TofCaliSvc", tofCaliSvc);
  if (scc ==  StatusCode::SUCCESS) { 
    log << MSG::INFO << "MrpcRecDBS Get Calibration Service Sucessfully!" << endreq;
  } else {
    log << MSG::ERROR << "MrpcRecDBS Get Calibration Service Failed !" << endreq;
    return StatusCode::FAILURE;
  }
 
  //Get TOF Raw Data Provider Service
  StatusCode scd = service("RawDataProviderSvc", tofDigiSvc);
  if (scd ==  StatusCode::SUCCESS) { 
    log << MSG::INFO << "MrpcRecDBS Get Tof Raw Data Service Sucessfully!" << endreq;
  } else {
    log << MSG::ERROR << "MrpcRecDBS Get Tof Raw Data Service Failed !" << endreq;
    return StatusCode::FAILURE;
  }
 
  if( m_printOutInfo ) { m_printOut = new MrpcDBSCount; }
 
 
 
 
 
  //Store data of MRPC RECONSTRUCTION                                                                                                                                               
  if(m_store_trans_time==1)
    {
      FileName = "output_mrpc_reconstructed.root";
      m_TreeFile=new TFile(TString(FileName),"RECREATE","TreeFile");
      m_tree_mrpc = new TTree("tree_mrpc_rec", "");   //Create Tree
      m_tree_mrpc->Branch("m_partId",&m_partID,"m_partId/D");
      m_tree_mrpc->Branch("m_stripidentifier",&m_stripidentifier,"m_stripidentifier/D");
      m_tree_mrpc->Branch("m_time_1",&m_time_1,"m_time_1/D");
      m_tree_mrpc->Branch("m_time_2",&m_time_2,"m_time_2/D");
      m_tree_mrpc->Branch("m_trans_time",&m_trans_time,"m_trans_time/D");
      m_tree_mrpc->Branch("m_r_ext",&m_r_ext,"m_r_ext/D");
      m_tree_mrpc->Branch("m_phi_ext",&m_phi_ext,"m_phi_ext/D");
      m_tree_mrpc->Branch("m_mrpc_dbs",&m_mrpc_dbs,"m_mrpc_dbs/D");
      m_tree_mrpc->Branch("m_mrpc_extrap",&m_mrpc_extrap,"m_extrap/D");
      m_tree_mrpc->Branch("m_mrpc_aver",&m_mrpc_aver,"m_mrpc_aver/D");
      m_tree_mrpc->Branch("m_transition_time_simulation",&m_transition_time_simulation,"m_transition_time_simulation/D");
 
    }
 
 
 
 
 
 
  return StatusCode::SUCCESS;
}
 
StatusCode MrpcDBSRec::beginRun(){
  MsgStream log(msgSvc(), name());
  log << MSG::INFO <<"MrpcRecDBS begin_run!"<< endreq;
  return StatusCode::SUCCESS;
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
 
StatusCode MrpcDBSRec::execute() {
 
 
 
  MsgStream log(msgSvc(), name());
  log << MSG::INFO << "MrpcRecDBS in execute()!" << endreq;
 
  SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
  if( !eventHeader ) {
    log << MSG::FATAL << "MrpcRecDBS could not find Event Header!" << endreq;
    return StatusCode::FAILURE;
  }
  int run = eventHeader->runNumber();
  int event = eventHeader->eventNumber();
  if( ( event % 1000 == 0 ) && m_printOutInfo ) {
    std::cout << "run:" << run << "  event: " << event << std::endl;
  }
  log << MSG::INFO << "run= " << run << "  event= " << event << endreq;
 
  MrpcDBSRecTDS tds;
  tds.RegisterNullRecTofTrackCol();
 
  // Magnetic Field: Colliding data and Cosmic Ray
  if( m_acceleratorStatus == "Colliding" ) {
    // Retrieve Event Start Time
    SmartDataPtr<RecEsTimeCol> estimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
    if( !estimeCol || ( estimeCol->size() == 0 ) ) { 
      log << MSG::WARNING << "MrpcRecDBS Could not find RecEsTimeCol!   Run = " << run << " Event = " << event << endreq;
      return StatusCode::SUCCESS;
    } 
    RecEsTimeCol::iterator iter_ESTime=estimeCol->begin();
    double t0  = (*iter_ESTime)->getTest();
    int t0Stat = (*iter_ESTime)->getStat();
    log << MSG::INFO << "t0= " << t0 << "  t0Stat= " << t0Stat << endreq;
 
 
 
 
    // Kalman Filter Track
    SmartDataPtr<RecMdcKalTrackCol> mdcKalTrackCol(eventSvc(),"/Event/Recon/RecMdcKalTrackCol");
    if( !mdcKalTrackCol ) {
      log << MSG::WARNING << "No MdcKalTrackCol in TDS!   Run = " << run << "  Event = " << event << endreq;
      return StatusCode::SUCCESS;
    }
 
 
 
    // Tof Get Extrapolated Track
    SmartDataPtr<RecExtTrackCol> extTrackCol(eventSvc(),"/Event/Recon/RecExtTrackCol");
    if( !extTrackCol ) {
      log << MSG::WARNING << "No ExtTrackCol in TDS!   Run = " << run << "  Event = " << event << endreq;
      return StatusCode::SUCCESS;
    }
    else 
      {
    if( m_printOutInfo ) { m_printOut->setExtTrackNum( extTrackCol->size() ); }
      }
 
    
 
 
    TofTrackVec* tofTrackVec = new TofTrackVec;
    RecExtTrackCol::iterator iter_track = extTrackCol->begin();
    for( ; iter_track < extTrackCol->end(); iter_track++ ) {
      RecMdcKalTrackCol::iterator iter_kal = mdcKalTrackCol->begin();
      for( ; iter_kal != mdcKalTrackCol->end(); iter_kal++ ) {
    if( (*iter_kal)->trackId() == (*iter_track)->trackId() ) break;
      }
      int kal[5] = {-1};
      if( iter_kal != mdcKalTrackCol->end() ) {
    for( unsigned int i=0; i<5; i++ ) {
      kal[i] = (*iter_kal)->getStat( 1, i );
    }
      }
      MrpcDBSTrack* tof = new MrpcDBSTrack;
      tof->setExtTrack( (*iter_track), kal, t0, t0Stat );
 
      if( tofTrackVec->size()>0 ) {
    std::vector<MrpcDBSTrack*>::iterator iterExt = tofTrackVec->begin();
    for( ; iterExt < tofTrackVec->end(); iterExt++ ) {
      if( (*iterExt)->isNoHit() ) continue;
      tof->getMultiHit( *iterExt );
    }
      }
 
      tofTrackVec->push_back( tof );
    }
 
 
    if( m_printOutInfo ) 
      {
      m_printOut->setTrack1Col( tofTrackVec );
      }
 
    // Retrieve Tof Digi Data
    TofDataMap tofDataMap = tofDigiSvc->tofDataMapTof( t0 );
    if( tofDataMap.empty() ) {
      log << MSG::WARNING << "No Tof Data Map in RawDataProviderSvc!   Run=" << run << " Event=" << event << endreq;
    }
    
 
    std::vector<MrpcDBSTrack*>::iterator iter = tofTrackVec->begin();
    for( ; iter < tofTrackVec->end(); iter++ ) {
      if( (*iter)->isNoHit() ) continue;
      (*iter)->setTofData( tofDataMap,m_neighborhood );
      if( m_printOutInfo ) { m_printOut->setTrack2( (*iter) ); }
      if( (*iter)->isNoHit() ) continue;
      (*iter)->match(tofTrackVec );
      if( m_printOutInfo ) { m_printOut->setTrack3( (*iter) ); }
    }
 
    iter = tofTrackVec->begin();
 
    for( ; iter < tofTrackVec->end(); iter++ ) 
      {
    (*iter)->setCalibration();
    if( m_printOutInfo ) { m_printOut->setTrack4( (*iter) ); }
 
 
 
 
    //Store MRPC Rootfile! //Variables are set in the Calibration Process!                                                                                                        
      if(m_store_trans_time==1)
      {
 
 
        //How was the transmission time calculated
        m_mrpc_dbs=(*iter)->mrpc_dbs();
        m_mrpc_extrap=(*iter)->mrpc_extra();
        m_mrpc_aver=(*iter)->mrpc_aver();
 
        m_partID=(*iter)->mrpc_partid();
        m_stripidentifier=(*iter)->mrpc_stripidentifier();
        m_time_1=(*iter)->mrpc_time_1();
        m_time_2=(*iter)->mrpc_time_2();
        m_trans_time=(*iter)->mrpc_trans_time();
 
        m_r_ext=(*iter)->mrpc_r_ext();
        m_phi_ext=(*iter)->mrpc_phi_ext();
 
 
        //Retrieve the simulated transitiontime                                                                                                                                   
        SmartDataPtr<Event::TofMcHitCol> mcParticleCol(eventSvc(), "/Event/MC/TofMcHitCol");
        double mc_truth_x,mc_truth_y;
        double mc_truth_trans_time;
        bool write_root_data=false;
        bool help=false;
        if(!mcParticleCol) help==true ;
 
 
 
 
        if (!mcParticleCol)
          {
        std::cout << "Could not retrieve McParticelCol" << std::endl;
          }
        else
          {
        
        Event::TofMcHitCol::iterator iter_mc = mcParticleCol->begin();
                for (; iter_mc != mcParticleCol->end(); iter_mc++)
                  {
 
 
 
                    mc_truth_x=(*iter_mc)->getPositionX();
                    mc_truth_y=(*iter_mc)->getPositionY();
                    const Identifier ident = (*iter_mc)->identify();
                    int mc_truth_partid=TofID::barrel_ec(ident);
                    int mc_truth_tof_id=TofID::phi_module(ident);
                    int mc_truth_tofmod=BesTofDigitizerEcV4_dbs::Get_module_mrpc_from_unique_identifier(mc_truth_tof_id);
                    int mc_truth_tofstr=BesTofDigitizerEcV4_dbs::Get_stripnumber_from_unique_identifier(mc_truth_tof_id);
                        
            mc_truth_trans_time = BesTofDigitizerEcV4_dbs::GetTransitionTime_extrap_track(mc_truth_x,mc_truth_y,mc_truth_partid,mc_truth_tofmod,mc_truth_tofstr);
              
 
 
 
 
              if(m_partID==mc_truth_partid && mc_truth_tofmod==((int)m_stripidentifier/25) && ((int)m_stripidentifier%25)==(mc_truth_tofstr))
                      {
            write_root_data=true;
                        break;
                      }
 
                  }
 
          }
        if(mc_truth_trans_time==0) m_transition_time_simulation=-100;
        else m_transition_time_simulation=mc_truth_trans_time;
 
        if(write_root_data) m_tree_mrpc->Fill();
      }
 
      
 
    }//close for loop
 
    tds.RegisterTDS(tofTrackVec);
 
    clearTofTrackVec( tofTrackVec );
 
// Check RecTofTrackCol Registered
    SmartDataPtr<RecTofTrackCol> tofTrackCol(eventSvc(),"/Event/Recon/RecTofTrackCol");  
    if (!tofTrackCol) { 
      log << MSG::FATAL << "MrpcRecDBS could not find RecTofTrackCol!" << endreq;
      return StatusCode::FAILURE;
    }
 
 
  }
  else {
    log << MSG::FATAL << "In MrpcRecDBS: AcceleratorStatus is NOT correct! m_acceleratorStatus = " << m_acceleratorStatus << endreq;
    return StatusCode::FAILURE;
  }
 
  return StatusCode::SUCCESS;
 
}
 
//   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
 
StatusCode MrpcDBSRec::finalize() {
  if( m_printOutInfo ) {
 
    m_printOut->final();
    delete m_printOut;
  }
 
  //Store data for MRPC                                                                                                                                                             
  if(m_store_trans_time==1) 
   {
      m_TreeFile->Write();
      m_TreeFile->Close();
    }
 
 
 
  return StatusCode::SUCCESS;
}
 
 
//----------------------------------------------------------------------------
 
void MrpcDBSRec::clearTofTrackVec( std::vector<MrpcDBSTrack*>*& tofTrackVec) {
  if( tofTrackVec ) {
    std::vector<MrpcDBSTrack*>::iterator iter = tofTrackVec->begin();
    for( ; iter < tofTrackVec->end(); iter++ ) {
      delete (*iter);
    }
    tofTrackVec->clear();
    delete tofTrackVec;
  }
  return;
}