ReadCosmicRayData.cxx

Go to the documentation of this file.

/**************************************************************************
 * CgemBOSS (BESIII Offline Software System)                              *
 *                                                                        *
 * Author: The BESII Collaboration                                        *
 * Contributors: Aiqiang Guo, Liangliang Wang, Lia Lavezzi                *
 * 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 = "input_file.txt");
  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");
  declareProperty("runNo",                 m_runNo = 1);
 
  //  mapper = new StripMapper(CosmicRayDataSetID);
}
 
ReadCosmicRayData::~ReadCosmicRayData(){
}
 
StatusCode ReadCosmicRayData::initialize(){
  MsgStream log(msgSvc(), name());
  log << MSG::INFO  << "ReadCosmicRayData initialize()" << endreq;
  
  //  bool ismap = mapper->FillMap();
 
  TString TDir_file(Dir_file);
  if(ReadDigi)
  {
      TString TTreeDigi(TreeDigi);
      Tdigi = new TChain(TTreeDigi); 
      TFileCollection* fc = new TFileCollection("mylist", "mylist",TDir_file);
      Tdigi->AddFileInfoList((TCollection*)fc->GetList());
 
      // 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_sheet", m_sheet);     // sheet
      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();
      cout<<"total Entry is "<<No_Entries_D<<endl;
      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;
    }
 
    //cout<<"layer "<<m_LayerID[ihit]<<", StripType "<<m_StripType[ihit]<<", original strip id: "<< m_strip[ihit]<<endl;
    // sheet ID/strip ID
    if(m_LayerID[ihit] == 0) {               // LAYER 1
        m_SheetID[ihit] = 0; // LAYER 1
        m_StripID[ihit] = m_strip[ihit]-1;// strip ID starts from 0
        if(m_StripType[ihit] == 0 && (m_StripID[ihit]<0||m_StripID[ihit]>=CgemID::getXSTRIP_MAX(m_LayerID[ihit]))) cout<<"ReadCosmicRayData::ConvertHitToDigi: X-strip ID out of range on layer1, sheet "<<m_SheetID[ihit]<<" : "<<m_StripID[ihit]<<endl;
        if(m_StripType[ihit] == 1 && (m_StripID[ihit]<0||m_StripID[ihit]>=CgemID::getVSTRIP_MAX(m_LayerID[ihit]))) cout<<"ReadCosmicRayData::ConvertHitToDigi: V-strip ID out of range on layer1, sheet "<<m_SheetID[ihit]<<" : "<<m_StripID[ihit]<<endl;
    }
    if(m_LayerID[ihit] == 1) {               // LAYER 2
        m_SheetID[ihit] = 0; 
        if(m_StripType[ihit] == 0) // X strip
        {
 
            // sheet down --> sheet 1
            if(m_strip[ihit] > CgemID::getXSTRIP_MAX(m_LayerID[ihit])) 
            {
                m_SheetID[ihit] = 0;
                m_StripID[ihit] = 2*CgemID::getXSTRIP_MAX(m_LayerID[ihit]) - m_strip[ihit];// strip ID starts from 0
            }
            else {
                // sheet up --> sheet 2
                m_StripID[ihit] = CgemID::getXSTRIP_MAX(m_LayerID[ihit]) - m_strip[ihit];
                m_SheetID[ihit] = 1;
            }
            if(m_StripID[ihit]<0||m_StripID[ihit]>=CgemID::getXSTRIP_MAX(m_LayerID[ihit])) cout<<"ReadCosmicRayData::ConvertHitToDigi: X-strip ID out of range on layer2, sheet "<<m_SheetID[ihit]<<" : "<<m_StripID[ihit]<<endl;
        }
        else if(m_StripType[ihit]  == 1)// V strip with correction for the wrong mapping at strip 538 and 1617 in integration
        {
            if(m_strip[ihit] == 539 || m_strip[ihit] == 1617) {
                cout<<"ReadCosmicRayData::ConvertHitToDigi: unexpated firing strip"<<m_strip[ihit]<<endl;
                return false;
            }
            else if(m_strip[ihit] < 539) {
                // sheet up --> sheet 2 
                m_StripID[ihit] = CgemID::getVSTRIP_MAX(m_LayerID[ihit]) - m_strip[ihit];
                m_SheetID[ihit] = 1;
            }
            else if(m_strip[ihit] < CgemID::getVSTRIP_MAX(m_LayerID[ihit]) + 2) {
                // sheet up --> sheet 2
                m_StripID[ihit] = CgemID::getVSTRIP_MAX(m_LayerID[ihit]) - m_strip[ihit] + 1;
                m_SheetID[ihit] = 1;
            }
            else if(m_strip[ihit] < 1617) {
                // sheet down --> sheet 1
                m_StripID[ihit] = 2*CgemID::getVSTRIP_MAX(m_LayerID[ihit]) - m_strip[ihit] + 1;
                m_SheetID[ihit] = 0;
            }
            else {
                // sheet down --> sheet 1
                m_StripID[ihit] = 2*CgemID::getVSTRIP_MAX(m_LayerID[ihit]) - m_strip[ihit] + 2;
                m_SheetID[ihit] = 0;
            }
            if(m_StripID[ihit]<0||m_StripID[ihit]>=CgemID::getVSTRIP_MAX(m_LayerID[ihit])) cout<<"ReadCosmicRayData::ConvertHitToDigi: V-strip ID out of range on layer2, sheet "<<m_SheetID[ihit]<<" : "<<m_StripID[ihit]<<endl;
        }
    }
    return true;
  }
  else if(CosmicRayDataSetID == "CR202001") {
    //    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;
    }
    
    if(m_strip[ihit] < 0) return false;
 
    // sheet ID/strip ID
    m_SheetID[ihit] = 0;
 
    if(m_StripType[ihit] == 0) // X strip
      {
    if(m_LayerID[ihit] == 0) {               // LAYER 1  
      m_StripID[ihit] = CgemID::getXSTRIP_MAX(m_LayerID[ihit]) - m_strip[ihit];
      m_SheetID[ihit] = 0;
    }
    else if(m_LayerID[ihit] == 1) {               // LAYER 2  
      // sheet down --> sheet 1
      if(m_strip[ihit] > CgemID::getXSTRIP_MAX(m_LayerID[ihit])) 
        {
          m_SheetID[ihit] = 0;
          m_StripID[ihit] = 2*CgemID::getXSTRIP_MAX(m_LayerID[ihit]) - m_strip[ihit];// strip ID starts from 0
        }
      else {
        // sheet up --> sheet 2
        m_StripID[ihit] = CgemID::getXSTRIP_MAX(m_LayerID[ihit]) - m_strip[ihit];
        m_SheetID[ihit] = 1;
      }
    }
    if(m_StripID[ihit]<0||m_StripID[ihit]>=CgemID::getXSTRIP_MAX(m_LayerID[ihit])) cout<<"ReadCosmicRayData::ConvertHitToDigi: X-strip ID out of range on layer2, sheet "<<m_SheetID[ihit]<<" : "<<m_StripID[ihit]<<endl;
      }
    else if(m_StripType[ihit]  == 1)// V strip with correction for the wrong mapping at strip 538 and 1617 in integration
      {
    if(m_LayerID[ihit] == 0) {               // LAYER 1
          m_StripID[ihit] = CgemID::getVSTRIP_MAX(m_LayerID[ihit]) - m_strip[ihit];
          m_SheetID[ihit] = 0;
        }
        else if(m_LayerID[ihit] == 1) {               // LAYER 2       
      // sheet down --> sheet 1
      if(m_strip[ihit] > CgemID::getVSTRIP_MAX(m_LayerID[ihit]))
        {
          m_SheetID[ihit] = 0;
          m_StripID[ihit] = 2*CgemID::getVSTRIP_MAX(m_LayerID[ihit]) - m_strip[ihit];// strip ID starts from 0
        }
      else {
        // sheet up --> sheet 2
        m_StripID[ihit] = CgemID::getVSTRIP_MAX(m_LayerID[ihit]) - m_strip[ihit];
        m_SheetID[ihit] = 1;
      }
    }
    if(m_StripID[ihit]<0||m_StripID[ihit]>=CgemID::getVSTRIP_MAX(m_LayerID[ihit])) cout<<"ReadCosmicRayData::ConvertHitToDigi: V-strip ID out of range on layer2, sheet "<<m_SheetID[ihit]<<" : "<<m_StripID[ihit]<<endl;
      }
    
    // 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;
  }
  else if(CosmicRayDataSetID == "NEW202001") {
 
    charge_channel = 1;
    time_channel = 1;
 
    // from plane to LayerID
    if(m_plane[ihit] != 0 && m_plane[ihit]  != 1) {
      cout<<"ReadCosmicRayData::ConvertHitToDigi: unknown plane "<<m_plane[ihit]<<endl;
      return false;
    }
    else m_LayerID[ihit] = m_plane[ihit] ;
    
    // 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;
    }
    
    if(m_strip[ihit] < 0) return false;
 
    // sheet
    m_SheetID[ihit] = m_sheet[ihit] ;
 
    // strip
    m_StripID[ihit] = m_strip[ihit];
 
    // some checks
    if(m_StripType[ihit] == 0 && m_StripID[ihit] > CgemID::getXSTRIP_MAX(m_LayerID[ihit])) {
    cout<<"ReadCosmicRayData::ConvertHitToDigi: X-strip ID out of range on layer " << m_LayerID[ihit] << ", sheet " << m_SheetID[ihit] << ": " << m_StripID[ihit] << endl;
    }
    else if(m_StripType[ihit] == 1 && m_StripID[ihit] > CgemID::getVSTRIP_MAX(m_LayerID[ihit])) {
      cout<<"ReadCosmicRayData::ConvertHitToDigi: V-strip ID out of range on layer " << m_LayerID[ihit] << ", sheet " << m_SheetID[ihit] << ": " << m_StripID[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;
 
    // set /Event/EventHeader
    SmartDataPtr<Event::EventHeader> eventHeader(m_evtSvc,"/Event/EventHeader");
    if (!eventHeader) {
        Event::EventHeader *eventHeader = new Event::EventHeader;
        StatusCode sc = m_evtSvc->registerObject("/Event/EventHeader",eventHeader);    
    }
    eventHeader->setEventNumber(Ind_Entry_D); 
    eventHeader->setRunNumber( m_runNo);
 
 
    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;
 
    // DEBUG
    /**
     const int nlayer = 3;                 // CHECK hardcoded 
     const int nsheet[nlayer] = {1, 2, 2}; // CHECK hardcoded
     const int nview = 2;                  // CHECK hardcoded 
     int nstrip[nlayer][nview] = {{856, 1173}, {630, 1077}, {832, 1395}}; // CHECK hardcoded 
     
     int layer=1;
     int sheet=0;
     int type=0;
     for(int strip=0; strip<nstrip[layer][type]; strip++) {
     cout << strip << " " 
     << mapper->GetGEMROC(strip, type, layer, sheet) << " "  
     << mapper->GetFEB(strip, type, layer, sheet) << " "
     << mapper->GetTIGER(strip, type, layer, sheet)
     << endl;
     }
  **/
  // --------
 
  return StatusCode::SUCCESS;
}