CgemDigitizerSvc.cxx

Go to the documentation of this file.

#include "CgemDigitizerSvc/CgemDigitizerSvc.h"
#include "GaudiKernel/Kernel.h"
#include "GaudiKernel/IInterface.h"
#include "GaudiKernel/StatusCode.h"
#include "GaudiKernel/SvcFactory.h"
#include "GaudiKernel/IJobOptionsSvc.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/PropertyMgr.h"
#include "GaudiKernel/IIncidentSvc.h"
#include "GaudiKernel/Incident.h"
#include "GaudiKernel/IIncidentListener.h"
#include "GaudiKernel/INTupleSvc.h"
#include "GaudiKernel/NTuple.h"
 
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/Bootstrap.h"
 
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/DataSvc.h"
 
#include "CalibData/CalibModel.h"
// #include "CalibData/Cgem/CgemCalibData.h"
// #include "CalibData/Cgem/CgemDataConst.h"
#include "EventModel/EventHeader.h"
#include "GaudiKernel/SmartDataPtr.h"
 
#include <iomanip>
#include <iostream>
#include <fstream>
#include <cmath>
 
#include "TObjArray.h"
 
using namespace std;
 
CgemDigitizerSvc::CgemDigitizerSvc( const string& name, ISvcLocator* svcloc) :Service (name, svcloc){
     // declare properties
     declareProperty("IonizationModel", m_ionModel);
     declareProperty("DriftAvalancheModel", m_driftAvaModel);
     declareProperty("InductionModel", m_inductionModel);
     declareProperty("GarfieldDebugging", m_garfDebugging);
     declareProperty("InductionDebugging", m_debugInduction);
     declareProperty("MagConfig", m_magConfig);
     //cout << "CgemDigitizerSvc::m_ionModel = " << m_ionModel << endl;
}
 
CgemDigitizerSvc::~CgemDigitizerSvc(){
}
 
StatusCode CgemDigitizerSvc::queryInterface(const InterfaceID& riid, void** ppvInterface){
     if( IID_ICgemDigitizerSvc.versionMatch(riid) ){
      *ppvInterface = static_cast<ICgemDigitizerSvc*> (this);
     } else{
      return Service::queryInterface(riid, ppvInterface);
     }
     return StatusCode::SUCCESS;
}
 
StatusCode CgemDigitizerSvc::initialize(){
     MsgStream log(messageService(), name());
     log << MSG::INFO << "CgemDigitizerSvc::initialize()" << endreq;
 
     StatusCode sc = Service::initialize();
     if( sc.isFailure() ) return sc;
 
     static IJobOptionsSvc* jobSvc = 0;
     if ( jobSvc == 0 ) {
      sc = service("JobOptionsSvc", jobSvc);
      if ( sc.isFailure() ) {
           std::cout << "Can't get the JobOptionsSvc @ DistBoss::GetPropertyValue()" << std::endl;
           return sc;
      }
     }
 
     // const vector<const Property*>* properties = jobSvc->getProperties(client);
     const vector<const Property*>* properties = jobSvc->getProperties("BesRndmGenSvc");
     // cout << "client of jobSvc: " << jobSvc->getClients() << endl;
     if ( properties == NULL ) {
      std::cout << "In CgemDigitizerSvc::initialize(), can't get client: " << std::endl;
      return StatusCode::FAILURE;
     }
 
     unsigned int randSeed;
     for ( unsigned int i = 0; i < properties->size(); ++i ) {
      if ( properties->at(i)->name() == "RndmSeed" ) {
           cout << "name of property: " << properties->at(i)->name() << endl;
           string strRnd = properties->at(i)->toString();
           sscanf(strRnd.c_str(), "%u", &randSeed);
           cout << "random seed from jobOption: " << randSeed << endl;
           break;
      }
     }
 
     sc = service("CgemGeomSvc", m_geomSvc); // initial pointer of CgemGeomSvc
     if(sc != StatusCode::SUCCESS)
     {
          log << MSG::ERROR << "can not use CgemGeomSvc" << endreq;
          return StatusCode::FAILURE;
     }
 
     gRandom->SetSeed(randSeed); // set random seed
 
     cout << "ionModel = " << m_ionModel << endl;
     if(1 == m_ionModel){
          m_pIon = new IonizationGar();
     } else if(2 == m_ionModel){
      m_pIon = new IonizationGTS();
     }
     // m_pIon = new IonizationGar();
     m_pIon->init(randSeed, m_geomSvc, m_magConfig);
     m_pIon->setDebugging(m_garfDebugging);
 
     if(1 == m_driftAvaModel){
          m_pDriftAndAva = new SamplingGar();
     } else if(2 == m_driftAvaModel){
       m_pDriftAndAva = new SamplingGTS();
     }
     m_pDriftAndAva->init(m_geomSvc, m_magConfig);
 
     if(1 == m_inductionModel){
          m_pInduction = new InductionGar();
     } else if(2 == m_inductionModel){
          m_pInduction = new InductionGTS();
     }
     m_pInduction->init(m_geomSvc, m_magConfig);
     m_pInduction->setDebugOutput(m_debugInduction);
 
     // ntuple
     INTupleSvc* ntupleSvc;
     Gaudi::svcLocator() -> service("NTupleSvc", ntupleSvc);
     NTuplePtr nt(ntupleSvc, "CgemDigitizerSvcNT/digi");
     // NTuplePtr nt(ntupleSvc(), "CgemDigitizerSvcNT/digi");
     if(nt) m_tuple = nt;
     else{
          m_tuple = ntupleSvc->book("CgemDigitizerSvcNT/digi", CLID_ColumnWiseTuple, "noise");
          if(m_tuple){
               m_tuple->addItem("nIonE",          m_ntNIonE);
           m_tuple->addItem("nMultiE",        m_ntNMultiE);
           m_tuple->addItem("nXstrips",       m_ntNxstrips,0,100);
           m_tuple->addItem("nVstrips",       m_ntNvstrips,0,100);
               m_tuple->addIndexedItem("XstripQ",   m_ntNxstrips,m_ntxstripQ);
               m_tuple->addIndexedItem("XstripT",   m_ntNxstrips,m_ntxstripT);
           m_tuple->addIndexedItem("VstripQ",   m_ntNvstrips,m_ntvstripQ);
               m_tuple->addIndexedItem("VstripT",   m_ntNvstrips,m_ntvstripT);
               m_tuple->addIndexedItem("XstripID",  m_ntNxstrips,m_XstripID);
               m_tuple->addIndexedItem("VstripID",  m_ntNvstrips,m_VstripID);
      }
     }
 
     return StatusCode::SUCCESS;
}
 
StatusCode CgemDigitizerSvc::finalize(){
     MsgStream log(messageService(), name());
     log << MSG::INFO << "CgemDigitizerSvc::finalize()" << endreq;
     delete m_pIon;
     delete m_pDriftAndAva;
     delete m_pInduction;
 
     return StatusCode::SUCCESS;
}
 
StatusCode CgemDigitizerSvc::setTrack(int layer, int particle, int charge, double p, double trkPosIn[], double trkPosOut[]){
 
  double time_spent = 0.;
  clock_t begin = clock();
 
     clear();
     m_layer = layer;
 
     m_pIon->setTrack(particle, charge, p, trkPosIn, trkPosOut);
 
     vector<double> ionEx;
     vector<double> ionEy;
     vector<double> ionEz;
     vector<double> ionEt;
     int nIonE = m_pIon->getNumberIonE();
     for(int i=0; i<nIonE; i++){
          ionEx.push_back(m_pIon->getEx(i));
          ionEy.push_back(m_pIon->getEy(i));
          ionEz.push_back(m_pIon->getEz(i));
          ionEt.push_back(m_pIon->getEt(i));
     }
     clock_t end = clock();
     time_spent += (double)(end - begin) / CLOCKS_PER_SEC;
      // cout << "CgemDigitizerSvc::ionization " << time_spent << " seconds" << endl;
     // ---------------------------
     time_spent = 0.;
     begin = clock();
 
     m_pDriftAndAva->setIonElectrons(layer, nIonE, ionEx, ionEy, ionEz, ionEt);
 
     m_nMultiE = m_pDriftAndAva->getNelectrons();
     for(int i=0; i<m_nMultiE; i++){
      double x = m_pDriftAndAva->getX(i);
      double y = m_pDriftAndAva->getY(i);
      double z = m_pDriftAndAva->getZ(i);
      double t = m_pDriftAndAva->getT(i);
          m_multiEx.push_back(x);
          m_multiEy.push_back(y);
          m_multiEz.push_back(z);
          m_multiEt.push_back(t);
     }
     end = clock();
     time_spent += (double)(end - begin) / CLOCKS_PER_SEC;
     // cout << "CgemDigitizerSvc::diffusion " << time_spent << " seconds" << endl;
     // --------------------
     time_spent = 0.;
     begin = clock();
 
     // Induction process
     // positionToStrips();
     // cout << "nXstrips = " << m_nXstrips << " ,   m_nVstrips = " << m_nVstrips << endl;
     m_pInduction->setMultiElectrons(layer, m_nMultiE, m_multiEx, m_multiEy, m_multiEz, m_multiEt);
     m_nXstrips = m_pInduction->getNXstrips();
     for(int i=0; i<m_nXstrips; i++){
      m_xstripSheet.push_back(m_pInduction->getXstripSheet(i));
      m_xstripID.push_back(m_pInduction->getXstripID(i));
      m_xstripQ.push_back(m_pInduction->getXstripQ(i));
      m_xstripT.push_back(m_pInduction->getXstripT(i));
     }
     m_nVstrips = m_pInduction->getNVstrips();
     for(int i=0; i<m_nVstrips; i++){
      m_vstripSheet.push_back(m_pInduction->getVstripSheet(i));
      m_vstripID.push_back(m_pInduction->getVstripID(i));
      m_vstripQ.push_back(m_pInduction->getVstripQ(i));
      m_vstripT.push_back(m_pInduction->getVstripT(i));
     }
     end = clock();
     time_spent += (double)(end - begin) / CLOCKS_PER_SEC;
     // cout << "CgemDigitizerSvc::induction " << time_spent << " seconds" << endl;
     // ---------------
     time_spent = 0.;
     begin = clock();
 
     //m_ntNIonE = nIonE;
     //m_ntNMultiE = m_nMultiE;
     //m_ntNxstrips = m_nXstrips;
     //m_ntNvstrips = m_nVstrips;
     //m_tuple->write();
 
     end = clock();
     time_spent += (double)(end - begin) / CLOCKS_PER_SEC;
     // cout << "CgemDigitizerSvc::write " << time_spent << " seconds" << endl;
 
 
     return StatusCode::SUCCESS;
}
 
void CgemDigitizerSvc::positionToStrips(){
     m_cgemLayer = m_geomSvc->getCgemLayer(m_layer);
     m_nSheet    = m_cgemLayer->getNumberOfSheet();
     int i_sheet;
 
     char hname[200];
     TObjArray* hlist = new TObjArray(0);
     // cout << "nMultiE :" << m_nMultiE << endl;
     for(int i=0; i<m_nMultiE; i++){
     // ofstream fout("elecPos.txt");
     // for(int i=0; i<1000; i++){
      // cout << "nMultiE :" << i << endl;
      double phi = atan2(m_multiEy[i], m_multiEx[i]);
      double z   = m_multiEz[i];
      G4ThreeVector pos(m_multiEx[i], m_multiEy[i], z);
      // fout << setw(5) << i << setw(20) << m_multiEx[i] << setw(20) << m_multiEy[i] << setw(20) << m_multiEz[i] << endl;
      for(int j=0; j<m_nSheet; j++)
      {
           m_cgemReadoutPlane = m_geomSvc->getReadoutPlane(m_layer, j);
           if(m_cgemReadoutPlane->OnThePlane(phi,z))
           {
            int iView = 1;
            int ID    = -1;
            double distV;
            int vID = m_cgemReadoutPlane->getClosestVStripID(pos, distV);
            ID = vID;
            if(vID>=0) { 
             if(distV>0) {
                  double distX; 
                  int xID = m_cgemReadoutPlane->getClosestXStripID(phi, distX);
                  if(xID>=0&&distX<distV) {
                   ID = xID;
                   iView = 0;
                  }
             }
             map<int, int>::iterator it = m_mapQ[j][iView].find(ID);
             if(it==m_mapQ[j][iView].end()) 
             {
                  m_mapQ[j][iView][ID]=1;
                  m_mapT[j][iView][ID]=m_multiEt[i];
             }
             else {
                  int Q = (it->second) +1;
                  m_mapQ[j][iView][ID]=Q;
                  if(m_multiEt[i] < m_mapT[j][iView][ID]) m_mapT[j][iView][ID] = m_multiEt[i];
             }
             // T info
             sprintf(hname, "time%d_%d_%d", j, iView, ID);
             TH1F* hist = (TH1F*)hlist->FindObject(hname);
             if(NULL == hist){
                  hist = new TH1F(hname, "", 1000, 0, 1000);
                  hlist->Add(hist);
             }
             hist->Fill(m_multiEt[i]);
            }
            break;
           }
      }
     }
     m_nXstrips = m_mapQ[0][0].size() + m_mapQ[1][0].size();
     m_nVstrips = m_mapQ[0][1].size() + m_mapQ[1][1].size();
     // fout.close();
 
     int nhist = hlist->GetEntries();
     // cout << setw(10) << m_nXstrips+m_nVstrips << setw(10) << nhist
     //       << "nXstrips = " << m_nXstrips << " ,   m_nVstrips = " << m_nVstrips << endl;
     for(int i=0; i<nhist; i++){
      TH1F* hist;
      TH1F* hlast;
      if(0==i){
           hist = (TH1F*)hlist->First();
      } else{
           hist = (TH1F*)hlist->After(hlast);
      }
      hlast = hist;
      const char* histname = hist->GetName();
      int iSheet, iView, ID;
      sscanf(histname, "%*4s%1d%*1s%1d%*1s%d", &iSheet, &iView, &ID);
 
      if(0==iView){
           m_xstripSheet.push_back(iSheet);
           m_xstripID.push_back(ID);
           m_xstripQ.push_back(m_mapQ[iSheet][iView][ID]);
           m_xstripT.push_back(m_mapT[iSheet][iView][ID]);
      } else{
           m_vstripSheet.push_back(iSheet);
           m_vstripID.push_back(ID);
           m_vstripQ.push_back(m_mapQ[iSheet][iView][ID]);
           m_vstripT.push_back(m_mapT[iSheet][iView][ID]);
      }
     }
 
     for(int i=0; i<nhist; i++){
      TH1F* hist = (TH1F*)hlist->Last();
      hlist->Remove(hist);
      delete hist;
     }
     hlist->Clear();
     delete hlist;
}
 
void CgemDigitizerSvc::clear(){
     m_nXstrips = 0;
     m_nVstrips = 0;
     m_xstripSheet.clear();
     m_xstripID.clear();
     m_vstripSheet.clear();
     m_vstripID.clear();
     m_xstripQ.clear();
     m_vstripQ.clear();
     m_xstripT.clear();
     m_vstripT.clear();
 
     m_nMultiE = 0;
     m_multiEx.clear();
     m_multiEy.clear();
     m_multiEz.clear();
     m_multiEt.clear();
 
     for(int i=0; i<2; i++){
      for(int k=0; k<2; k++){
           m_mapQ[i][k].clear();
           m_mapT[i][k].clear();
      }
     }
}