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"
#include "TSystem.h"
 
using namespace std;
 
CgemDigitizerSvc::CgemDigitizerSvc( const string& name, ISvcLocator* svcloc) :Service (name, svcloc){
     // declare properties
     declareProperty("IonizationModel", m_ionModel=1);
     declareProperty("DriftAvalancheModel", m_driftAvaModel=1); // 3 mean SamplingGar2, 
     declareProperty("InductionModel", m_inductionModel=1);     // 3 mean InductionGar2,
                                                                // this 2 doesnot extend from base normally and have to be used in conjunction.
     declareProperty("GarfieldDebugging", m_garfDebugging=false);
     declareProperty("SamplingDebugging", m_samplingDebugging=false);
     declareProperty("InductionDebugging", m_debugInduction=false);
     declareProperty("MagConfig", m_magConfig=1);
     declareProperty("LUTFilePath",m_LUTFilePath = "/bes3fs/cgemCosmic/data/CGEM_cosmic_look_up_table_17.root");
     declareProperty("QBranchVsampleDelay",sample_delay = 162.5);
     declareProperty("StoreUnderThrethold",m_storeFlag = false);
     declareProperty("Saturation",m_saturation = true);
     declareProperty("SaveNt",m_saveNt= false);
     declareProperty("SamplingGarGainMultiplier",m_GainMultiplier);
     declareProperty("SamplingGarTransMultiplier",m_TransMultiplier= 1.0);
     declareProperty("SamplingGarDiffuMultiplier",m_DiffuMultiplier= 1.0);
     //cout << "CgemDigitizerSvc::m_ionModel = " << m_ionModel << endl;
     declareProperty("Ngaps_microSector",m_Ngaps_microSector= 40);
     declareProperty("gapShift_microSector",m_gapShift_microSector);
     declareProperty("gap_microSector", m_gap_microSector= 1.1);// in mm
     declareProperty("microSector_width", m_microSector_width);// in rad
     declareProperty("QinGausSigma", m_QinGausSigma);
     declareProperty("ScaleSignalX", m_ScaleSignalX=1.0);
}
 
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->setDebugging(m_garfDebugging);
    m_pIon->init(randSeed, m_geomSvc, m_magConfig);
 
    if ((3==m_driftAvaModel && 3!=m_inductionModel )||(3!=m_driftAvaModel && 3==m_inductionModel )) {
        cout<<"CgemDigitizerSvc::initialize() Error: DriftAvalancheModel type3 and InductionModel type3 must be used together!"<<endl;
        return StatusCode::FAILURE;
    }
 
    if(1 == m_driftAvaModel){
        m_pDriftAndAva = new SamplingGar();
        SamplingGar *pDriftAndAvaGar=(SamplingGar *) m_pDriftAndAva;
        pDriftAndAvaGar->setGainMultiplier(m_GainMultiplier);
        pDriftAndAvaGar->setTransMultiplier(m_TransMultiplier);
        pDriftAndAvaGar->setDiffuMultiplier(m_DiffuMultiplier);
    } else if(2 == m_driftAvaModel){
        m_pDriftAndAva = new SamplingGTS();
    } else if(3 == m_driftAvaModel){
        m_pDriftAndAva = new SamplingGar2();
        SamplingGar2* pDriftAndAva=(SamplingGar2*)m_pDriftAndAva;
        pDriftAndAva->setMultiElectronMapAddr(&m_hitmap);
        
        // migrate of samGar1;
        pDriftAndAva->setGainMultiplier(m_GainMultiplier);
        pDriftAndAva->setTransMultiplier(m_TransMultiplier);
        pDriftAndAva->setDiffuMultiplier(m_DiffuMultiplier);
        // migrate of indGar1
        pDriftAndAva->setNgapsSect(m_Ngaps_microSector);
        pDriftAndAva->setGapSizeSect(m_gap_microSector);
        pDriftAndAva->setGapShiftSect(m_gapShift_microSector);
        pDriftAndAva->setMicroSectorWidthRad(m_microSector_width);
        
    }
    m_pDriftAndAva->init(m_geomSvc, m_magConfig);
    m_pDriftAndAva->setDebugging(m_samplingDebugging);
 
    if(1 == m_inductionModel){
        InductionGar* indGar = new InductionGar();
        indGar->setNgapsSect(m_Ngaps_microSector);
        indGar->setGapSizeSect(m_gap_microSector);
        indGar->setGapShiftSect(m_gapShift_microSector);
        indGar->setMicroSectorWidthRad(m_microSector_width);
        indGar->setQinGausSigma(m_QinGausSigma);
        indGar->setScaleSignalX(m_ScaleSignalX);
        m_pInduction = indGar;
        //m_pInduction = new InductionGar();
        m_pInduction->setVsampleDelay(sample_delay);
        m_pInduction->setLUTFilePath(m_LUTFilePath);
        m_pInduction->setStoreFlag(m_storeFlag);
        m_pInduction->setSaturation(m_saturation);
    } else if(2 == m_inductionModel){
        m_pInduction = new InductionGTS();
        //m_pInduction = new InductionGar();
    }else if (3 == m_inductionModel){
        m_pInduction = new InductionGar2();
        InductionGar2* indGar=(InductionGar2*)m_pInduction;
        indGar->setQinGausSigma(m_QinGausSigma);
        indGar->setScaleSignalX(m_ScaleSignalX);
 
        m_pInduction->setVsampleDelay(sample_delay);
        m_pInduction->setLUTFilePath(m_LUTFilePath);
        m_pInduction->setStoreFlag(m_storeFlag);
        m_pInduction->setSaturation(m_saturation);
    }
    m_pInduction->init(m_geomSvc, m_magConfig);
    m_pInduction->setDebugOutput(m_debugInduction);
 
    // --- check gaps
    cout<<"CgemDigitizerSvc:"<<endl;
    cout<<"gapShift_microSector="<<m_gapShift_microSector[0]
        <<", "<<m_gapShift_microSector[1]
        <<", "<<m_gapShift_microSector[2]
        <<", "<<m_gapShift_microSector[3]
        <<", "<<m_gapShift_microSector[4]
        <<", "<<m_gapShift_microSector[5]
        <<endl;
    cout<<"microSector_width="<<m_microSector_width[0]
        <<", "<<m_microSector_width[1]
        <<", "<<m_microSector_width[2]
        <<", "<<m_microSector_width[3]
        <<", "<<m_microSector_width[4]
        <<", "<<m_microSector_width[5]
        <<endl;
    cout<<"QinGausSigma="<<m_QinGausSigma[0]
        <<", "<<m_QinGausSigma[1]
        <<endl;
    cout<<"Ngaps_microSector="<<m_Ngaps_microSector<<endl;
    cout<<"gap_microSector="<<m_gap_microSector<<endl;
    
    // --- read saturations
    TFile *LUTFile = TFile::Open(m_LUTFilePath.c_str(),"read");
    TTree *tree = (TTree*)LUTFile->Get("tree");
    Float_t QDC_saturation;
    Int_t layer,sheet,strip_v,strip_x;
    tree->SetBranchAddress("strip_x_boss", &strip_x);
    tree->SetBranchAddress("strip_v_boss", &strip_v);
    tree->SetBranchAddress("layer", &layer);
    tree->SetBranchAddress("sheet", &sheet);
    tree->SetBranchAddress("calib_QDC_saturation", &QDC_saturation);
    for(int i=0;i<tree->GetEntries();i++) {
        tree->GetEntry(i);
        if(strip_x!=-1) {
            Qsaturation[layer][sheet][0][strip_x] = QDC_saturation;
        }
        if(strip_v!=-1) {
            Qsaturation[layer][sheet][1][strip_v] = QDC_saturation;
        }
    }
    for(int i=0;i<832;i++) {
        Qsaturation[2][0][0][i] = 45.;
        Qsaturation[2][1][0][i] = 45.;
    }
    for(int i=0;i<1395;i++) {
        Qsaturation[2][0][1][i] = 45.;
        Qsaturation[2][1][1][i] = 45.;
    }
    LUTFile->Close();
 
    // ntuple
    if(m_saveNt) {
        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[]){
    vector<int> Vparticle             (1,particle);
    vector<int> Vcharge               (1,charge);
    vector<double> Vp                 (1,p);
 
    // awkward way of saying 
    //vector<vector<double> >VtrkPosIn={{trkPosIn[0], trkPosIn[1], trkPosIn[2]}};
    vector<double> posin;
    posin.push_back(trkPosIn[0]);
    posin.push_back(trkPosIn[1]);
    posin.push_back(trkPosIn[2]);
    vector<vector<double> > VtrkPosIn;
    VtrkPosIn.push_back(posin);
 
    vector<double> posout;
    posout.push_back(trkPosOut[0]);
    posout.push_back(trkPosOut[1]);
    posout.push_back(trkPosOut[2]);
    vector<vector<double> > VtrkPosOut;
    VtrkPosOut.push_back(posout);//= {{trkPosOut[0], trkPosOut[1], trkPosOut[2]}};
 
    StatusCode sc=setTrack(layer,Vparticle, Vcharge, Vp, VtrkPosIn, VtrkPosOut);
 
    return sc;
 
    //*******below are original contents of this function.**********
 
    //  //ProcInfo_t info;
    //  //gSystem->GetProcInfo(&info);
    //  //double currentMemory = info.fMemResident/1024;// MB
    //  //cout<<"CgemDigitizerSvc::setTrack() memory 1 "<<currentMemory<<" MB"<<endl;
    //
    //  double time_spent = 0.;
    //  clock_t begin = clock();
    //  //cout<<"setTrack"<<endl;
    //  clear();
    //  m_layer = layer;
    //
    //  //cout<<"before m_pIon->setTrack()"<<endl;
    //  //cout<<"set particle "<<particle<<", p="<<p<<endl;
    //  m_pIon->setTrack(particle, charge, p, trkPosIn, trkPosOut);
    //  //cout<<"after m_pIon->setTrack()"<<endl;
    //  //gSystem->GetProcInfo(&info); currentMemory = info.fMemResident/1024; cout<<"CgemDigitizerSvc::setTrack() memory 2 "<<currentMemory<<" MB"<<endl;
    //
    //  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));
    //  }
    //  //cout<<"after pushing back ionE "<<endl;
    //  if(m_saveNt) m_ntNIonE=nIonE;
    //  clock_t end = clock();
    //  //cout<<" CLOCKS_PER_SEC "<<endl;
    //  //cout<<CLOCKS_PER_SEC<<endl;
    //  //cout<<" calculate time_spent"<<endl;
    //  time_spent += (double)(end - begin) / CLOCKS_PER_SEC;
    //  //cout << "CgemDigitizerSvc::ionization " << time_spent << " seconds" << endl;
    //  // ---------------------------
    //  time_spent = 0.;
    //  begin = clock();
    //
    //  //gSystem->GetProcInfo(&info); currentMemory = info.fMemResident/1024; cout<<"CgemDigitizerSvc::setTrack() memory 3 "<<currentMemory<<" MB"<<endl;
    //  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_pmultiEx).push_back(x);
    //  //  (*m_pmultiEy).push_back(y);
    //  //  (*m_pmultiEz).push_back(z);
    //  //  (*m_pmultiEt).push_back(t);
    //  //}
    //  m_pmultiEx=& (m_pDriftAndAva->getXContainer());
    //  m_pmultiEy=& (m_pDriftAndAva->getYContainer());
    //  m_pmultiEz=& (m_pDriftAndAva->getZContainer());
    //  m_pmultiEt=& (m_pDriftAndAva->getTContainer());
    //  //std::cout<<m_nMultiE<<"\t";
    //  if(m_saveNt) m_ntNMultiE=m_nMultiE;
    //  end = clock();
    //  time_spent += (double)(end - begin) / CLOCKS_PER_SEC;
    //  //cout << "CgemDigitizerSvc::diffusion " << time_spent << " seconds" << endl;
    //  // --------------------
    //  time_spent = 0.;
    //  begin = clock();
    //
    //  //gSystem->GetProcInfo(&info); currentMemory = info.fMemResident/1024; cout<<"CgemDigitizerSvc::setTrack() memory 4 "<<currentMemory<<" MB"<<endl;
    //  // Induction process
    //  // positionToStrips();
    //  //cout<<" before m_pInduction->setMultiElectrons "<<endl;
    //  m_pInduction->setMultiElectrons(layer, m_nMultiE, *m_pmultiEx, *m_pmultiEy, *m_pmultiEz, *m_pmultiEt);
    //  //cout<<" after  m_pInduction->setMultiElectrons "<<endl;
    //  m_nXstrips = m_pInduction->getNXstrips();
    //  //gSystem->GetProcInfo(&info); currentMemory = info.fMemResident/1024; cout<<"CgemDigitizerSvc::setTrack() memory 5 "<<currentMemory<<" MB"<<endl;
    //  //m_ntNxstrips = m_nXstrips;
    //  //m_ntNvstrips = m_nVstrips;
    //  //  cout<<"read data"<<endl;
    //  for(int i=0; i<m_nXstrips; i++){
    //      //m_ntNxstrips=i;
    //      if(m_saveNt && i<100) {
    //          m_ntxstripQ[i]=m_pInduction->getXstripQ(i);
    //          m_ntxstripT[i]=m_pInduction->getXstripT(i);
    //          m_XstripID[i]=m_pInduction->getXstripID(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_xstripQ2.push_back(m_pInduction->getXstripQ_Branch(i));
    //      m_xstripT.push_back(m_pInduction->getXstripT(i));
    //      m_xfirstT.push_back(m_pInduction->getXfirstT(i));
    //
    //  }
    //
    //  m_nVstrips = m_pInduction->getNVstrips();
    //  for(int i=0; i<m_nVstrips; i++){
    //      //m_ntNvstrips=i;
    //      if(m_saveNt && i<100) {
    //          m_ntvstripQ[i]=m_pInduction->getVstripQ(i);
    //          m_ntvstripT[i]=m_pInduction->getVstripT(i);
    //          m_VstripID[i]=m_pInduction->getVstripID(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_vstripQ2.push_back(m_pInduction->getVstripQ_Branch(i));
    //      m_vstripT.push_back(m_pInduction->getVstripT(i));
    //      m_vfirstT.push_back(m_pInduction->getVfirstT(i));
    //  }
    //  end = clock();
    //  time_spent += (double)(end - begin) / CLOCKS_PER_SEC;
    //  // cout << "CgemDigitizerSvc::induction " << time_spent << " seconds" << endl;
    //  // cout << "nXstrips = " << m_nXstrips << " ,   m_nVstrips = " << m_nVstrips << endl;
    //  // ---------------
    //  time_spent = 0.;
    //  begin = clock();
    //
    //  if(m_saveNt) {
    //      //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;
    //
    //  //std::cout<<"wow test2!\t";
    //  //checkMemorySize();
    //  //gSystem->GetProcInfo(&info);
    //  //double currentMemory2 = info.fMemResident/1024;// MB
    //  //cout<<"CgemDigitizerSvc::setTrack() memory from "<<currentMemory<<" to "<<currentMemory2<<" MB"<<endl;
    //
    //  return StatusCode::SUCCESS;
}
 
StatusCode CgemDigitizerSvc::setTrack(int layer, vector<int>particle, vector<int> charge, vector<double> p, vector<vector<double> > trkPosIn, vector<vector<double> > trkPosOut){
 
    //ProcInfo_t info;
    //gSystem->GetProcInfo(&info);
    //double currentMemory = info.fMemResident/1024;// MB
    //cout<<"CgemDigitizerSvc::setTrack() memory 1 "<<currentMemory<<" MB"<<endl;
    
    double time_spent = 0.;
    clock_t begin = clock();
    int nTotIonE = 0;
    //cout<<"setTrack"<<endl;
    clear();
    m_layer = layer;
    vector<double> ionEx;
    vector<double> ionEy;
    vector<double> ionEz;
    vector<double> ionEt;
    for(unsigned TrackSeg=0;TrackSeg<particle.size();TrackSeg++){
 
        double SegTrkPosIn[3];
        double SegTrkPosOut[3];
        for(int i=0;i<3;i++){
            SegTrkPosIn[i]=trkPosIn[TrackSeg][i];
            SegTrkPosOut[i]=trkPosOut[TrackSeg][i];
        }
        //cout<<"trrpos    "<<SegTrkPosIn[0]<<"    "<<SegTrkPosIn[1]<<"    "<<SegTrkPosIn[2]<<"    "<<SegTrkPosOut[0]<<"    "<<SegTrkPosOut[1]<<"    "<<SegTrkPosOut[2]<<endl;
        m_pIon->setTrack(particle[TrackSeg], charge[TrackSeg], p[TrackSeg], SegTrkPosIn, SegTrkPosOut);
 
        int nIonE = m_pIon->getNumberIonE();
        nTotIonE=nIonE+nTotIonE;
        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));
        }
    }
    //cout<<"after pushing back ionE "<<endl;  
    if(m_saveNt) m_ntNIonE=nTotIonE;
    clock_t end = clock();
    time_spent += (double)(end - begin) / CLOCKS_PER_SEC;
     //cout << "CgemDigitizerSvc::ionization " << time_spent << " seconds" << endl;
    // ---------------------------
    time_spent = 0.;
    
 
    if (m_driftAvaModel !=3){
        begin = clock();
        m_pDriftAndAva->setIonElectrons(layer, nTotIonE, 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_pmultiEx).push_back(x);
        //  (*m_pmultiEy).push_back(y);
        //  (*m_pmultiEz).push_back(z);
        //  (*m_pmultiEt).push_back(t);
        //}
        m_pmultiEx=& (m_pDriftAndAva->getXContainer());
        m_pmultiEy=& (m_pDriftAndAva->getYContainer());
        m_pmultiEz=& (m_pDriftAndAva->getZContainer());
        m_pmultiEt=& (m_pDriftAndAva->getTContainer());
 
        if(m_saveNt) {
            m_ntNMultiE=m_nMultiE;
        }
        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_pmultiEx, *m_pmultiEy, *m_pmultiEz, *m_pmultiEt);
 
    }
    else{ //m_driftAvaModel ==3
        SamplingGar2 *pSampling2=(SamplingGar2 *) m_pDriftAndAva;
        InductionGar2 *pInduction2=(InductionGar2 *) m_pInduction;
 
        pSampling2->setIonElectrons(layer, nTotIonE, ionEx, ionEy, ionEz, ionEt);
        //yields m_hitmap
        // as test;
        end = clock();
        time_spent += (double)(end - begin) / CLOCKS_PER_SEC;
 
    //  cout << "CgemDigitizerSvc::diffusion " << time_spent << " seconds" << endl;
        // --------------------
        m_nMultiE = m_pDriftAndAva->getNelectrons();
 
        time_spent = 0.;
        begin = clock();
        if(m_saveNt) {
            m_ntNMultiE=m_nMultiE;
        }
 
        InductionGar2::PVectorXYZT vMultiElectrons;
        vMultiElectrons.x=& (m_pDriftAndAva->getXContainer());
        vMultiElectrons.y=& (m_pDriftAndAva->getYContainer());
        vMultiElectrons.z=& (m_pDriftAndAva->getZContainer());
        vMultiElectrons.t=& (m_pDriftAndAva->getTContainer());
        InductionGar2::PVectorXYZT *pMultiElectronsInfo=&vMultiElectrons;
        if (!m_samplingDebugging){
            pMultiElectronsInfo=0;
        }
        pInduction2->setMultiElectrons(layer, m_hitmap,pMultiElectronsInfo);
 
    }
    m_nXstrips = m_pInduction->getNXstrips();
 
    for(int i=0; i<m_nXstrips; i++){
        // m_ntNxstrips=i;
        if(m_saveNt && i<100) {
            m_ntxstripQ[i]=m_pInduction->getXstripQ(i);
            m_ntxstripT[i]=m_pInduction->getXstripT(i);
            m_XstripID[i]=m_pInduction->getXstripID(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_xstripQ2.push_back(m_pInduction->getXstripQ_Branch(i));
        m_xstripT.push_back(m_pInduction->getXstripT(i));
        m_xfirstT.push_back(m_pInduction->getXfirstT(i));
 
 
    }
 
    m_nVstrips = m_pInduction->getNVstrips();
    for(int i=0; i<m_nVstrips; i++){
        // m_ntNvstrips=i;
        if(m_saveNt && i<100) {
            m_ntvstripQ[i]=m_pInduction->getVstripQ(i);
            m_ntvstripT[i]=m_pInduction->getVstripT(i);
            m_VstripID[i]=m_pInduction->getVstripID(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_vstripQ2.push_back(m_pInduction->getVstripQ_Branch(i));
        m_vstripT.push_back(m_pInduction->getVstripT(i));
        m_vfirstT.push_back(m_pInduction->getVfirstT(i));
    }
    end = clock();
 
 
    //cout<<"read finish"<<endl;
    
    time_spent += (double)(end - begin) / CLOCKS_PER_SEC;
    //cout << "CgemDigitizerSvc::induction " << time_spent << " seconds" << endl;
    //cout<< "CgemDigitizerSvc multiplied electrons " << m_nMultiE <<  endl;
    // ---------------
    time_spent = 0.;
    begin = clock();
 
    if(m_saveNt) {
        //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;
    //checkMemorySize();
 
    return StatusCode::SUCCESS;
}
// ZZH i think  this function should have been depreciated?
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_pmultiEy)[i], (*m_pmultiEx)[i]);
        double z   = (*m_pmultiEz)[i];
        G4ThreeVector pos((*m_pmultiEx)[i], (*m_pmultiEy)[i], z);
        // fout << setw(5) << i << setw(20) << (*m_pmultiEx)[i] << setw(20) << (*m_pmultiEy)[i] << setw(20) << (*m_pmultiEz)[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_pmultiEt)[i];
                    }
                    else {
                        int Q = (it->second) +1;
                        m_mapQ[j][iView][ID]=Q;
                        if((*m_pmultiEt)[i] < m_mapT[j][iView][ID]) m_mapT[j][iView][ID] = (*m_pmultiEt)[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_pmultiEt)[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_xstripQ2.clear();
    m_vstripQ2.clear();
    m_xfirstT.clear();
    m_vfirstT.clear();
    m_nMultiE = 0;
    // now that we dont own this. they are cleared by Sampling Obj.
    //(*m_pmultiEx).clear(); 
    //(*m_pmultiEy).clear(); 
    //(*m_pmultiEz).clear(); 
    //(*m_pmultiEt).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();
        }
    }
}
 
 
void CgemDigitizerSvc::checkMemorySize()
{
    ProcInfo_t info;
    gSystem->GetProcInfo(&info);
    double currentMemory = info.fMemResident/1024;// MB
    cout<<"CgemDigitizerSvc: current memory size "<<currentMemory<<" MB"<<endl;
}