TestHit.cxx

Go to the documentation of this file.

// **************************************************************************/
// authors: L. Lavezzi (univ. of Torino & INFN, Italy)
//
 
//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 "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/TestHit.h"
#include "TMath.h"
//using namespace std;
 
TestHit::TestHit(const std::string& name, ISvcLocator* pSvcLocator):
    Algorithm(name,pSvcLocator){   
 
  declareProperty("LUTfile", lutfile = "/bes3fs/cgemCosmic/data/CGEM_cosmic_look_up_table_from_17_to_17.root");
  declareProperty("TimeMin", timemin = -8900); // noise or signal
  declareProperty("TimeMax", timemax = -8500); // 
 
}
 
TestHit::~TestHit(){
 
  // delete tree;
  // delete output;
  // delete lutreader;
}
 
StatusCode TestHit::initialize(){
  
  MsgStream log(msgSvc(), name());
  log << MSG::INFO << "TestHit initialize()" << endreq;
 
  // CgemGeomSvc 
  StatusCode sc;       
 
  sc = service("CgemGeomSvc", m_geoSvc); 
  if(sc != StatusCode::SUCCESS)  { 
    log << MSG::ERROR << "can not use CgemGeomSvc" << endreq; 
    return StatusCode::FAILURE;
  }  
 
  event = 0;
  
  output = new TFile("hit.root", "RECREATE");
  tree = new TTree("tree","hit info tree");
  tree->Branch("event", &event, "event/I");
  tree->Branch("nhit", &nhit, "nhit/I");
  tree->Branch("nhit_L1_S1_x", &nhit_L1_S1_x,"nhit_L1_S1_x/I");
  tree->Branch("nhit_L2_S1_x", &nhit_L2_S1_x,"nhit_L2_S1_x/I");
  tree->Branch("nhit_L2_S2_x", &nhit_L2_S2_x,"nhit_L2_S2_x/I");
  tree->Branch("nhit_L1_S1_v", &nhit_L1_S1_v,"nhit_L1_S1_v/I");
  tree->Branch("nhit_L2_S1_v", &nhit_L2_S1_v,"nhit_L2_S1_v/I");
  tree->Branch("nhit_L2_S2_v", &nhit_L2_S2_v,"nhit_L2_S2_v/I");
  //
  tree->Branch("ntwin_L1_S1_x", &ntwin_L1_S1_x,"ntwin_L1_S1_x/I");
  tree->Branch("ntwin_L2_S1_x", &ntwin_L2_S1_x,"ntwin_L2_S1_x/I");
  tree->Branch("ntwin_L2_S2_x", &ntwin_L2_S2_x,"ntwin_L2_S2_x/I");
  tree->Branch("ntwin_L1_S1_v", &ntwin_L1_S1_v,"ntwin_L1_S1_v/I");
  tree->Branch("ntwin_L2_S1_v", &ntwin_L2_S1_v,"ntwin_L2_S1_v/I");
  tree->Branch("ntwin_L2_S2_v", &ntwin_L2_S2_v,"ntwin_L2_S2_v/I");
  // from geometry
  tree->Branch("hit_strip",   &hit_strip,  "hit_strip[nhit]/I");
  tree->Branch("hit_view",  &hit_view,  "hit_view[nhit]/I");
  tree->Branch("hit_layer", &hit_layer, "hit_layer[nhit]/I");
  tree->Branch("hit_sheet", &hit_sheet, "hit_sheet[nhit]/I");
  tree->Branch("hit_length", &hit_length, "hit_length[nhit]/D");
  // from LUT
  tree->Branch("hit_channel", &hit_channel,"hit_channel[nhit]/I");
  tree->Branch("hit_roc",     &hit_roc,    "hit_roc[nhit]/I");
  tree->Branch("hit_feb",     &hit_feb,    "hit_feb[nhit]/I");
  tree->Branch("hit_tiger",   &hit_tiger,  "hit_tiger[nhit]/I");
  tree->Branch("hit_chip",    &hit_chip,   "hit_chip[nhit]/I");
  // from signal
  tree->Branch("hit_t",     &hit_t,     "hit_t[nhit]/D");
  tree->Branch("hit_q",     &hit_q,     "hit_q[nhit]/D");
  tree->Branch("hit_saturated",  &hit_saturated,  "hit_saturated[nhit]/I");
  tree->Branch("hit_quality", &hit_quality, "hit_quality[nhit]/I");
 
  // LUT
  lutreader = new CgemLUTReader(lutfile);
  lutreader->ReadLUT();
 
  return StatusCode::SUCCESS;
}
 
void TestHit::reset() {
  
  event = 0;
  nhit = 0;
  nhit_L1_S1_x= 0 ;
  nhit_L2_S1_x= 0 ;
  nhit_L2_S2_x= 0 ;
  nhit_L1_S1_v= 0 ;
  nhit_L2_S1_v= 0 ;
  nhit_L2_S2_v= 0 ;
 
  ntwin_L1_S1_x= 0 ;
  ntwin_L2_S1_x= 0 ;
  ntwin_L2_S2_x= 0 ;
  ntwin_L1_S1_v= 0 ;
  ntwin_L2_S1_v= 0 ;
  ntwin_L2_S2_v= 0 ;
 
  for(int ihit = 0; ihit < MAXNOFHIT; ihit++) {
    hit_strip[ihit] = -1;
    hit_view[ihit]  = -1;
    hit_layer[ihit] = -1;
    hit_sheet[ihit] = -1;
    hit_length[ihit] = 0.;
 
    hit_channel[ihit] = -1;
    hit_roc[ihit] = -1;
    hit_feb[ihit] = -1;
    hit_tiger[ihit] = -1;
    hit_chip[ihit] = -1;
    hit_quality[ihit] = -1;
 
    hit_t[ihit] = 0 ;
    hit_q[ihit] = 0 ;
    hit_saturated[ihit] = 0;
  }
  
}
 
StatusCode TestHit::execute(){
  
  MsgStream log(msgSvc(), name());
    //    cout << "->TestHit::execute" << endl;
 
    //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;
 
    // -------------------
                 
 
    // CgemGeoReadoutPlane* anode = m_SvcCgem->getReadoutPlane(0, 0);
    // anode_radius_L1_x = anode->getRX();
    // anode_radius_L1_v = anode->getRV();
    // anode_mid_gap_L1 = anode->getMidRAtGap();
    // anode = m_SvcCgem->getReadoutPlane(1, 0);
    // anode_radius_L2_x = anode->getRX();
    // anode_radius_L2_v = anode->getRV();
    // anode_mid_gap_L2 = anode->getMidRAtGap();
 
    // -------------------
   
 
    reset();
    // -------------
 
    // retrieve CGEM digits from TDS
    SmartDataPtr<CgemDigiCol> aDigiCol(m_evtSvc,"/Event/Digi/CgemDigiCol");
    if(!aDigiCol)
      cout<<"Could not retrieve CGEM digi collection"<<endl;
    
    nhit = aDigiCol->size();
    CgemDigiCol::iterator iDigiCol;
    
    //    event = m_evtSvc->GetEventID(); // CHECK
 
    // loop on hits
    int ihit=0;
    for(iDigiCol=aDigiCol->begin(); iDigiCol!=aDigiCol->end(); iDigiCol++)
      {
    const Identifier ident = (*iDigiCol)->identify();
 
    int layer = CgemID::layer(ident);
    int sheet = CgemID::sheet(ident);
    int strip = CgemID::strip(ident);
    int view = 1;
    bool is_xstrip = CgemID::is_xstrip(ident);
    if(is_xstrip == true) view = 0;
    double charge = (*iDigiCol)->getCharge_fc();
    double time = (*iDigiCol)->getTime_ns();
    
    hit_strip[ihit] = strip;
    hit_layer[ihit] = layer;
    hit_sheet[ihit] = sheet;
    hit_view[ihit] = view;
 
    anode = m_geoSvc->getReadoutPlane(layer, sheet);
    if(view == 1)  hit_length[ihit] = anode->getVStripLength(strip);
    else hit_length[ihit] = anode->getLength();
 
    hit_channel[ihit] = lutreader->GetChannel(layer, sheet, view, strip);
    hit_roc[ihit]     = lutreader->GetROC(layer, sheet, view, strip);
    hit_feb[ihit]     = lutreader->GetFEB(layer, sheet, view, strip);
    hit_tiger[ihit]   = lutreader->GetTIGER(layer, sheet, view, strip);
    hit_chip[ihit]    = lutreader->GetChip(layer, sheet, view, strip);
    hit_quality[ihit] = lutreader->GetQuality(layer, sheet, view, strip);
    
    //  std::cout << ihit << " strip " << strip << " ch " << hit_channel[ihit] << " roc " << hit_roc[ihit] << " feb " << hit_feb[ihit] << " tg " << hit_tiger[ihit] << std::endl;
 
    hit_t[ihit] = time;
    hit_q[ihit] = charge;
    hit_saturated[ihit] = 0; // CHECK
 
    if(is_xstrip == true) {
      if(layer==0 && sheet==0)  nhit_L1_S1_x++;
      else if(layer==1 && sheet==0)  nhit_L2_S1_x++;
          else if(layer==1 && sheet==1)  nhit_L2_S2_x++;
    }
    else {
      if(layer==0 && sheet==0)  nhit_L1_S1_v++;
      else if(layer==1 && sheet==0)  nhit_L2_S1_v++;
      else if(layer==1 && sheet==1)  nhit_L2_S2_v++;
    }
 
    // --- time window
    if(time > timemin && time < timemax) {
      if(is_xstrip == true) {
        if(layer==0 && sheet==0)  ntwin_L1_S1_x++;
        else if(layer==1 && sheet==0)  ntwin_L2_S1_x++;
        else if(layer==1 && sheet==1)  ntwin_L2_S2_x++;
      }
      else {
        if(layer==0 && sheet==0)  ntwin_L1_S1_v++;
        else if(layer==1 && sheet==0)  ntwin_L2_S1_v++;
        else if(layer==1 && sheet==1)  ntwin_L2_S2_v++;
      }
    }
    // -----------
    ihit++;
      }
 
    tree->Fill();
    event++;
   
    return StatusCode::SUCCESS;
}
 
StatusCode TestHit::finalize(){
  
    MsgStream log(msgSvc(),name());
    log << MSG::INFO << "TestHit finalize()" << endreq;
    output->Write();
    output->Close();
  
    return StatusCode::SUCCESS;
}