bak_ReadCosmicRayData-00-00-06/src/ReadCosmicRayData.cxx

Go to the documentation of this file.

/**************************************************************************
 * CgemBOSS (BESIII Offline Software System)                              *
 *                                                                        *
 * Author: The BESII Collaboration                                        *
 * Contributors: Aiqiang Guo, Liangliang Wang                             *
 * Date: Dec 10 2018                                                      *
 *                                                                        *
 **************************************************************************/
/**************************************************************************
 * Note for 01
 * Change the reading method for vectorized data
 * By: Aiqiang Guo
 * Date: Dec 17 2018
 **************************************************************************/
/**************************************************************************
 * Note for 02
 * Add CC and mTPC for cluster object
 * By: Aiqiang Guo
 * Date: Dec 27 2018
 **************************************************************************/
 
// include system lib
#include <iostream>
#include <iomanip>
#include <string>
#include <cmath>
 
// Include files
#include "GaudiKernel/AlgFactory.h"
#include "GaudiKernel/DataObject.h"
#include "GaudiKernel/IEventProcessor.h"
 
#include "GaudiKernel/Incident.h"
#include "GaudiKernel/IIncidentSvc.h"
#include "GaudiKernel/Memory.h"
 
#include <csignal>
 
// for save digit & cluster
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/RegistryEntry.h"
#include "GaudiKernel/MsgStream.h"
 
#include "CgemRawEvent/CgemDigi.h"
#include "CgemRecEvent/RecCgemCluster.h"
 
#include "Identifier/CgemID.h"
 
#include "RawEvent/RawDataUtil.h"
#include "RawEvent/DigiEvent.h"
#include "ReconEvent/ReconEvent.h"
#include "EventModel/EventHeader.h"
#include "GaudiKernel/SmartDataPtr.h"
 
#include "ReadCosmicRayData/ReadCosmicRayData.h"
 
//using namespace std;
ReadCosmicRayData::ReadCosmicRayData(const std::string& name, ISvcLocator* pSvcLocator):
  Algorithm(name,pSvcLocator){
 
  declareProperty("Dir_file",              Dir_file = "Cosmic_data_01.root");
  declareProperty("TreeDigi",              TreeDigi = "t1");
  declareProperty("TreeCluster",           TreeCluster = "t1");
  declareProperty("ReadDigi",              ReadDigi = true);
  declareProperty("ReadCluster",           ReadCluster = true);
  declareProperty("DigiSheetID",           DigiSheetID = 0);
  declareProperty("Cut_on_tpc",            Cut_on_tpc = false);
  declareProperty("ClusterSheetID",        ClusterSheetID = 0);
  declareProperty("ClusterRecZ",           ClusterRecZ = 0);
  declareProperty("R_Cluster",             R_Cluster = 1.0);
  declareProperty("Shift_DigitLayerID",    Shift_DigitLayerID = 0);
  declareProperty("Shift_DigitSheetID",    Shift_DigitSheetID = 0);
  declareProperty("Shift_DigitXStripID",   Shift_DigitXStripID = 0);
  declareProperty("Shift_DigitVStripID",   Shift_DigitVStripID = 0);
  declareProperty("Shift_ClusterLayerID",  Shift_ClusterLayerID = 0);
  declareProperty("Shift_ClusterSheetID",  Shift_ClusterSheetID = 0);
  declareProperty("Shift_RecPhi",          Shift_RecPhi = 0);
  declareProperty("Shift_RecV",            Shift_RecV = 0);
  declareProperty("Shift_RecZ",            Shift_RecZ = 0);
  declareProperty("CosmicRayDataSetID",    CosmicRayDataSetID = "CR201909");
 
}
 
ReadCosmicRayData::~ReadCosmicRayData(){
}
 
StatusCode ReadCosmicRayData::initialize(){
  MsgStream log(msgSvc(), name());
  log << MSG::INFO  << "ReadCosmicRayData initialize()" << endreq;
  
  TString TDir_file(Dir_file);
  f  = new TFile(TDir_file);
  if(ReadDigi)
    {
      TString TTreeDigi(TreeDigi);
      Tdigi  = (TTree*)f->Get(TTreeDigi);
      
      // Get Cgem digi tree 
      Tdigi->SetBranchAddress("Event", &m_Event_D);        // event ID
      Tdigi->SetBranchAddress("nGemHit", &m_nGemHit);      // nof GEM hits
      // Tdigi->SetBranchAddress("GemHit_nHit", &m_nGemHit); // nof GEM hits it is the same thing as before 
      
      // information on the IDs
      Tdigi->SetBranchAddress("GemHit_channel", m_channel); // channel no. [0, 63]
      Tdigi->SetBranchAddress("GemHit_ROC", m_ROC);         // ROC no.
      Tdigi->SetBranchAddress("GemHit_chip", m_chip);       // chip no.
      Tdigi->SetBranchAddress("GemHit_FEB", m_FEB);         // FEB no.
      Tdigi->SetBranchAddress("GemHit_plane", m_plane);     // plane
      Tdigi->SetBranchAddress("GemHit_view", m_view);       // view (axial or stereo strips)
      Tdigi->SetBranchAddress("GemHit_strip", m_strip);     // strip no.
 
      // physical information
      Tdigi->SetBranchAddress("GemHit_saturated", m_saturated); // is the ASIC channel saturated
      Tdigi->SetBranchAddress("GemHit_q", m_charge);            // charge (fC)
      Tdigi->SetBranchAddress("GemHit_time", m_time);           // time (ns) 
      //        Tdigi->SetBranchAddress("GemHit_is_tpc", m_GemHit_is_tpc);
 
      No_Entries_D = Tdigi->GetEntries();
      Ind_Entry_D = 0;
    }
 
  /**
     if(ReadCluster)
     {
     //Get Cgem cluster tree 
     TString TTreeCluster(TreeCluster);
     Tcluster  = (TTree*)f->Get(TTreeCluster);
 
     Tcluster->SetBranchAddress("Event", &m_Event_C);
     Tcluster->SetBranchAddress("GemCluster1d_nCluster", &m_nGemCluster);
     Tcluster->SetBranchAddress("GemCluster1d_nHit", m_ClusternHit);
     Tcluster->SetBranchAddress("GemCluster1d_HitIndex", m_ClusterHitIndex);
     Tcluster->SetBranchAddress("GemCluster1d_plane", m_ClusterLayerID);
     Tcluster->SetBranchAddress("GemCluster1d_view", m_Flag);
     Tcluster->SetBranchAddress("GemCluster1d_q", m_EnergyDeposit);
     Tcluster->SetBranchAddress("GemCluster1d_x", m_Cluster_x);
     Tcluster->SetBranchAddress("GemCluster1d_z", m_Cluster_z);
     Tcluster->SetBranchAddress("GemCluster1d_x_cc", m_Cluster_x_cc);
     Tcluster->SetBranchAddress("GemCluster1d_x_tpc", m_Cluster_x_tpc);
     Tcluster->SetBranchAddress("GemCluster1d_z_cc", m_Cluster_z_cc);
     Tcluster->SetBranchAddress("GemCluster1d_z_tpc", m_Cluster_z_tpc);
     //Tcluster->SetBranchAddress("ClusterSheetID", m_ClusterSheetID);
     //Tcluster->SetBranchAddress("ClusterFlagB", m_ClusterFlagB);
     //Tcluster->SetBranchAddress("ClusterFlagE", m_ClusterFlagE);
     //Tcluster->SetBranchAddress("RecV", m_RecV);
     //Tcluster->SetBranchAddress("RecZ", m_RecZ);
 
     No_Entries_C = Tcluster->GetEntries();
     Ind_Entry_C = 0;
     }
  **/
 
  return StatusCode::SUCCESS;
}
 
 
int ReadCosmicRayData::TranslateDigitLayerID(int Input_LayerID)
{
  int ShiftValue = Shift_DigitLayerID;
  return Input_LayerID+ShiftValue;
}
 
int ReadCosmicRayData::TranslateDigitSheetID(int Input_SheetID)
{
  int ShiftValue = Shift_DigitSheetID;
  return Input_SheetID+ShiftValue;
}
 
int ReadCosmicRayData::TranslateDigitXStripID(int Input_StripID)
{
  int ShiftValue = Shift_DigitXStripID;
  return Input_StripID+ShiftValue;
}
 
int ReadCosmicRayData::TranslateDigitVStripID(int Input_StripID)
{
  int ShiftValue = Shift_DigitVStripID;
  return Input_StripID+ShiftValue;
}
 
int ReadCosmicRayData::TranslateDigitStripID(int Input_StripID, int StripType)
{
  int Output_StripID = -1;
  if(StripType==0) Output_StripID =  TranslateDigitXStripID(Input_StripID);
  if(StripType==1) Output_StripID =  TranslateDigitVStripID(Input_StripID);
  return Output_StripID;
}
 
int ReadCosmicRayData::TranslateDigitStripType(int Input_StripType)
{
  return Input_StripType;
}
 
int ReadCosmicRayData::TranslateClusterLayerID(int Input_LayerID)
{
  int ShiftValue = Shift_ClusterLayerID;
  return Input_LayerID+ShiftValue;
}
 
int ReadCosmicRayData::TranslateClusterSheetID(int Input_SheetID)
{
  int ShiftValue = Shift_ClusterSheetID;
  return Input_SheetID+ShiftValue;
}
 
int ReadCosmicRayData::TranslateClusterFlag(int Input_Flag)
{
  return Input_Flag-2;
}
 
double ReadCosmicRayData::TranslateRecPhi(double Input_RecPhi)
{
  double ShiftValue = Shift_RecPhi;
  return Input_RecPhi+ShiftValue;
}
 
double ReadCosmicRayData::TranslateRecV(double Input_RecV)
{
  double ShiftValue = Shift_RecV;
  return Input_RecV+ShiftValue;
}
 
double ReadCosmicRayData::TranslateRecZ(double Input_RecZ)
{
  double ShiftValue = Shift_RecZ;
  return Input_RecZ+ShiftValue;
}
 
void ReadCosmicRayData::ReadCgemDigits()
{
  //cout<<"Start to read the CgemDigits"<<endl;
  Tdigi->GetEntry(Ind_Entry_D);
  Ind_Entry_D++;
}
 
void ReadCosmicRayData::ReadCgemClusters()
{
  //cout<<"Start to read the CgemClusters"<<endl;
  Tcluster->GetEntry(Ind_Entry_C);
  Ind_Entry_C++;
}
 
bool ReadCosmicRayData::ConvertHitToDigi(int ihit, unsigned int &charge_channel, unsigned int &time_channel)
{
 
  if(CosmicRayDataSetID == "CR201909") {
    //    cout << "ReadCosmicRayData::ConvertGRAALToCgemBoss(), converting set " << CosmicRayDataSetID << endl;
 
    //    if(!m_GemHit_is_tpc[i]&&Cut_on_tpc) return false; // ??
    charge_channel = 1;
    time_channel = 1;
 
    // from plane to LayerID
    if(m_plane[ihit] == 0) m_LayerID[ihit] = 0; // LAYER1
    else if(m_plane[ihit] == 1) m_LayerID[ihit] = 1; // LAYER2
    else if(m_plane[ihit] == 2) m_LayerID[ihit] = 1; // LAYER2
    else {
        cout<<"ReadCosmicRayData::ConvertHitToDigi: unknown plane "<<m_plane[ihit]<<endl;
        return false;
    }
                                                                     
    // from view to StripType
    if(m_view[ihit] == 2) m_StripType[ihit] = 0;       // x strip (phi)
    else if(m_view[ihit] == 3) m_StripType[ihit] = 1;  // v strip (stereo)
    else {
        cout<<"ReadCosmicRayData::ConvertHitToDigi: unknown view "<<m_view[ihit]<<endl;
        return false;
    }
 
    // sheet ID/strip ID
    m_SheetID[ihit] = 0; // LAYER 1, or sheet 1 of layer 2
    m_StripID[ihit] = m_strip[ihit]-1;// strip ID starts from 0
    if(m_LayerID[ihit] == 1) {               // LAYER 2
      if(m_StripType[ihit] == 0) // X strip
      {
          if(m_strip[ihit] > CgemID::getXSTRIP_MAX(m_LayerID[ihit])) // sheet 2
          {
              m_SheetID[ihit] = 1;
              m_StripID[ihit] = m_strip[ihit] - 1 - CgemID::getXSTRIP_MAX(m_LayerID[ihit]);// strip ID starts from 0
          }
      }
      else if(m_StripType[ihit]  == 1)// V strip
      {
          if( m_strip[ihit] >= 539 && m_strip[ihit]-1<=CgemID::getVSTRIP_MAX(m_LayerID[ihit]) )
          {
              m_StripID[ihit] = m_strip[ihit] - 2; // correction for the wrong mapping at strip 539 in integration
          }
          else if( m_strip[ihit]-1 > CgemID::getVSTRIP_MAX(m_LayerID[ihit]) ) // sheet 2
          {
              m_SheetID[ihit] = 1;
              m_StripID[ihit] = m_strip[ihit] - 2 - CgemID::getVSTRIP_MAX(m_LayerID[ihit]);
              if(m_strip[ihit] >= 1617) m_StripID[ihit] = m_strip[ihit] - 3 - CgemID::getVSTRIP_MAX(m_LayerID[ihit]);// correction for the wrong mapping at strip 1617 in integration
          }
      }
    }
 
    // cout << "GRAAL " << m_strip[ihit] << " view " << m_view[ihit] << " plane " << m_plane[ihit] << endl;
    // cout << "CGEMB " << m_StripID[ihit] << " view " << m_StripType[ihit] << " plane " << m_LayerID[ihit] << " sheet " << m_SheetID[ihit] << endl;
 
 
    return true;
  }
 
  cout << "ERROR : ReadCosmicRayData::ConvertGRAALToCgemBoss(), the data set " << CosmicRayDataSetID << " is unknown! " << endl;
  return false;
 
}
 
void ReadCosmicRayData::SaveCgemDigits()
{
 
  bool printFlag=false;
  //cout<<"Start to save CgemDigi"<<endl;
  // cgem digis collection defined in BOSS 
  CgemDigiCol* aCgemDigiCol = new CgemDigiCol;
    
  //cout<<"nDigi : "<<m_nGemHit<<endl;
  if (m_nGemHit > 0)
    {
      // push back cgem digits to CgemDigiCol in BOSS 
      for(int i=0;i<m_nGemHit;i++)
    {
      unsigned int charge_channel;
      unsigned int time_channel;
      bool is_converted = ConvertHitToDigi(i, charge_channel, time_channel);
      if(!is_converted) {
          cout<<"ReadCosmicRayData::SaveCgemDigits failed to convert hit "<<i<<" to digi in event "<<m_Event_D<<endl;
          continue;
      }
      const Identifier ident = CgemID::strip_id(
                            TranslateDigitLayerID(m_LayerID[i]),
                            TranslateDigitSheetID(m_SheetID[i]),
                            TranslateDigitStripType(m_StripType[i]),
                            TranslateDigitStripID(m_StripID[i],TranslateDigitStripType(m_StripType[i])));
 
      /**
         cout << "graal strip " << m_strip[i] 
         << " graal view " << m_view[i]
         << " stripid " << m_StripID[i] 
         << " striptype " << m_StripType[i] 
         << " strip " << TranslateDigitStripID(m_StripID[i],TranslateDigitStripType(m_StripType[i]))
         << " identifier " << CgemID::getIntID(TranslateDigitLayerID(m_LayerID[i]),
         TranslateDigitSheetID(m_SheetID[i]),
         TranslateDigitStripType(m_StripType[i]),
         TranslateDigitStripID(m_StripID[i],TranslateDigitStripType(m_StripType[i]))) << endl;;
      **/
      
      int layerid = CgemID::layer(ident);
      int stripid = CgemID::strip(ident);
          int sheetid = CgemID::sheet(ident);
      bool is_xstrip = CgemID::is_xstrip(ident);
 
      if(printFlag) 
        {
          cout  
        << "   digiID="  << ident.get_value()
        << "   layerID=" << CgemID::layer(ident)
        << "   sheetID=" << CgemID::sheet(ident)
        << "   isXstrip="<< CgemID::is_xstrip(ident)
        << "   stripID=" << CgemID::strip(ident)
        << "   time channel="    << time_channel                 
        << "   charge channel="  << charge_channel           << endl;
        }
      
      CgemDigi* aCgemDigi = new CgemDigi(ident, time_channel, charge_channel);  
      aCgemDigi->setTime_ns(m_time[i]);
      aCgemDigi->setCharge_fc(m_charge[i]);
      aCgemDigiCol->push_back(aCgemDigi);      
        
    } // End of 'for(int i=0;i<nDigi;i++)' 
    } // End of 'if (nDigi > 0)' 
 
  // register CGEM digits collection to TDS 
  StatusCode scCgem = m_evtSvc->registerObject("/Event/Digi/CgemDigiCol", aCgemDigiCol);
  if(scCgem!=StatusCode::SUCCESS) 
    {
      cout << "ERROR : ReadCosmicRayData::SaveCgemDigits(), Could not register CGEM digi collection! " << endl;
    }
 
  //retrieve CGEM digits from TDS
  /**
     SmartDataPtr<CgemDigiCol> aDigiCol(m_evtSvc,"/Event/Digi/CgemDigiCol");
     if(!aDigiCol)
     cout<<"Could not retrieve CGEM digi collection"<<endl;
       
     CgemDigiCol::iterator iDigiCol;
     for(iDigiCol=aDigiCol->begin(); iDigiCol!=aDigiCol->end(); iDigiCol++)
     {
     const Identifier ident = (*iDigiCol)->identify();
     cout<<" layer: "<<CgemID::layer(ident)
     <<" sheet: "<<CgemID::sheet(ident) 
     <<" strip: "<<CgemID::strip(ident) 
     <<" charge: "<<(*iDigiCol)->getCharge_fc()
     <<" time: "<<(*iDigiCol)->getTime_ns()<<endl;
     }
     cout<<"end of retrieve CGEM digi collection"<<endl;
  **/
}
 
 
 
void ReadCosmicRayData::SaveCgemClusters()
{
 
  bool printFlag=false;
  //cout<<"Start to save CgemCluster"<<endl;
 
  // cgem digis collection defined in BOSS 
  RecCgemClusterCol* aRecCgemClusterCol = new RecCgemClusterCol;
    
  int nCluster = m_nGemCluster; // the number of Clusters in one event
  if (nCluster > 0)
    {
      // push back cgem clusters to RecCgemClusterCol in BOSS 
      for(int i=0;i<nCluster;i++)
    {
      RecCgemCluster* aRecCgemCluster = new RecCgemCluster();  
 
      aRecCgemCluster->setclusterid(i);
      aRecCgemCluster->setlayerid(TranslateClusterLayerID(m_ClusterLayerID[i]));
      aRecCgemCluster->setsheetid(ClusterSheetID);
      aRecCgemCluster->setflag(TranslateClusterFlag(m_Flag[i]));
      aRecCgemCluster->setenergydeposit(m_EnergyDeposit[i]);
      if(TranslateClusterFlag(m_Flag[i])==0) 
        {
          aRecCgemCluster->setrecphi(m_Cluster_x[i]);
          aRecCgemCluster->setrecphi_CC(m_Cluster_x_cc[i]);
          aRecCgemCluster->setrecphi_mTPC(m_Cluster_x_tpc[i]);
            }
      if(TranslateClusterFlag(m_Flag[i])==1) 
        {
          aRecCgemCluster->setrecv(m_Cluster_x[i]);
          aRecCgemCluster->setrecv_CC(m_Cluster_x_cc[i]);
          aRecCgemCluster->setrecv_mTPC(m_Cluster_x_tpc[i]);
            }
      aRecCgemCluster->setRecZ(ClusterRecZ); //currently, no Z information is available
      aRecCgemCluster->setRecZ_CC(ClusterRecZ);
      aRecCgemCluster->setRecZ_mTPC(ClusterRecZ);
      //aRecCgemCluster->setRecZ(m_Cluster_z[i]);
      //aRecCgemCluster->setRecZ_CC(m_Cluster_z_cc[i]);
      //aRecCgemCluster->setRecZ_mTPC(m_Cluster_z_tpc[i]);
      aRecCgemCluster->setclusterflag(m_ClusterHitIndex[i][0],m_ClusterHitIndex[i][m_ClusternHit[i]-1]);
 
      if(printFlag)
        {
          cout  
        << "   clusterID="            << aRecCgemCluster->getclusterid()
        << "   clusterlayerID="       << aRecCgemCluster->getlayerid()
        << "   clustersheetID="       << aRecCgemCluster->getsheetid()
        << "   flag="                 << aRecCgemCluster->getflag()
        << "   energydeposit="        << aRecCgemCluster->getenergydeposit()
        << "   recphi="               << aRecCgemCluster->getrecphi()
        << "   recv="                 << aRecCgemCluster->getrecv()
        << "   recZ="                 << aRecCgemCluster->getRecZ()
        << "   clusterflagb="         << aRecCgemCluster->getclusterflagb()
        << "   clusterflage="         << aRecCgemCluster->getclusterflage()<< endl;
        }
 
      aRecCgemClusterCol->push_back(aRecCgemCluster);      
 
    }  
    } 
 
 
  // register CGEM clusters collection to TDS 
  StatusCode scCgem = m_evtSvc->registerObject("/Event/Recon/RecCgemClusterCol", aRecCgemClusterCol); //Where should I put the vector?
  if(scCgem!=StatusCode::SUCCESS) 
    {
      cout << "ERROR : ReadCosmicRayData:::SaveCgemClusters(), Could not register CGEM cluster collection! " << endl;
    }
 
}
 
 
StatusCode ReadCosmicRayData::execute(){
  
  MsgStream log(msgSvc(), name());
  if(ReadDigi&&!ReadCluster) log << MSG::INFO << "ReadCosmicRayData execute(): "<<Ind_Entry_D+1<<"/"<<No_Entries_D<<" events are finished !" << endreq;
  if(!ReadDigi&&ReadCluster) log << MSG::INFO << "ReadCosmicRayData execute(): "<<Ind_Entry_C+1<<"/"<<No_Entries_C<<" events are finished !" << endreq;
  if(ReadDigi&&ReadCluster) log << MSG::INFO << "ReadCosmicRayData execute(): "<<Ind_Entry_C+1<<"/"<<No_Entries_C<<" events are finished !" << endreq;
 
  //interface to event data service
  ISvcLocator* svcLocator = Gaudi::svcLocator();
  StatusCode sc=svcLocator->service("EventDataSvc", m_evtSvc);
  if (sc.isFailure())
    cout<<"Could not accesss EventDataSvc!"<<endl;
 
  if(ReadDigi)
    {
      DigiEvent* aDigiEvent = new DigiEvent;
      sc = m_evtSvc->registerObject("/Event/Digi",aDigiEvent);
      if(sc!=StatusCode::SUCCESS) {
    cout<< "Could not register DigiEvent" <<endl;
      }
 
      ReadCgemDigits();
      SaveCgemDigits();
      //cout<<"Ind_Entry_D "<<Ind_Entry_D<<", Max_Ind_Entry_D "<<No_Entries_D<<endl;
 
        if(Ind_Entry_D==No_Entries_D)
        {
                log << MSG::INFO << "scheduling a event processing stop...." << endreq;
                SmartIF<IEventProcessor> ep(serviceLocator());
                if (ep) ep->stopRun();
        }
 
    }
  if(ReadCluster)
    {
        ReconEvent* aReconEvent = new ReconEvent;
        sc = m_evtSvc->registerObject("/Event/Recon",aReconEvent);
        if(sc!=StatusCode::SUCCESS) {
            cout<< "Could not register ReconEvent" <<endl;
        }
 
        ReadCgemClusters();
        SaveCgemClusters();
        //cout<<"Ind_Entry_C "<<Ind_Entry_C<<", Max_Ind_Entry_C "<<No_Entries_C<<endl;
        if(Ind_Entry_C==No_Entries_C)
        {
            log << MSG::INFO << "scheduling a event processing stop...." << endreq;
            SmartIF<IEventProcessor> ep(serviceLocator());
            if (ep) ep->stopRun();
        }
    }
  return StatusCode::SUCCESS;
}
 
StatusCode ReadCosmicRayData::finalize(){
  MsgStream log(msgSvc(),name());
  log << MSG::INFO << "ReadCosmicRayData finalize()" << endreq;
 
  return StatusCode::SUCCESS;
}