SamplingGar.cxx

Go to the documentation of this file.

#include "CgemDigitizerSvc/SamplingGar.h"
 
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/DataSvc.h"
 
#include "CLHEP/Units/PhysicalConstants.h"
 
#include "G4DigiManager.hh"
#include "Randomize.hh"
#include "G4ios.hh"
 
#include "CLHEP/Units/PhysicalConstants.h"
 
#include <iomanip>
#include <iostream>
#include <fstream>
#include <cmath>
 
#include "TRandom.h"
#include "TMath.h"
 
using namespace std;
 
SamplingGar::SamplingGar(){
    m_GainMultiplier[0]=1.60;
    m_GainMultiplier[1]=1.60;
    m_GainMultiplier[2]=1.60;
}
 
SamplingGar::~SamplingGar(){
}
 
void SamplingGar::init(ICgemGeomSvc* geomSvc, double magConfig){
     m_geomSvc = geomSvc;
     m_magConfig = magConfig;
 
     readSmearPar();
 
     // initialize polya functions
     char funname[10];
     for(int l=0; l<3;l++){
         for(int i=0; i<3; i++){
             sprintf(funname, "funPolya%d%d", l, i);
             m_funPolya[l][i] = new TF1(funname, "[0]*pow(1+[2],1+[2])*pow(x/[1],[2])*exp(-(1+[2])*x/[1])/TMath::Gamma(1+[2])", 1, 500);
             m_funPolya[l][i]->SetParameter(0, m_gain_gem[i][0]);
             m_funPolya[l][i]->SetParameter(1, m_gain_gem[i][1]*m_GainMultiplier[l]);
             m_funPolya[l][i]->SetParameter(2, m_gain_gem[i][2]);
         }
     }
}
 
void SamplingGar::setIonElectrons(int layer, int nElectrons, std::vector<double> x, std::vector<double> y, std::vector<double> z, std::vector<double> t){
     clear();
     m_nIonE = nElectrons;
     // cout << "nIonE = " << m_nIonE << endl;
     for(int i=0; i<nElectrons; i++){
      double r = sqrt(x[i]*x[i] + y[i]*y[i]);
      double phi;
      if(y[i] >= 0) phi = acos(x[i] / r);
      else phi = CLHEP::twopi - acos(x[i] / r);
      m_vecIonR.push_back(r);
      m_vecIonPhi.push_back(phi);
      m_vecIonZ.push_back(z[i]);
      m_vecIonT.push_back(t[i]);
 
      // cout << "check geom layer" << setw(15) << layer << setw(15) << m_geomSvc->getCgemLayer(layer)->getInnerROfGapD()
      //      << setw(15) << m_geomSvc->getCgemLayer(layer)->getOuterROfGapD()
      //      << setw(15) << m_geomSvc->getCgemLayer(layer)->getCgemFoil(0)->getInnerROfCgemFoil() << endl;
 
      double radii_gem1 = m_geomSvc->getCgemLayer(layer)->getCgemFoil(0)->getInnerROfCgemFoil();
      double driftD = radii_gem1 - r;
      // cout << "driftD = " << setw(15) << driftD << setw(15) << radii_gem1 << setw(15) << r << endl;
      samplingDrift(layer,driftD, i);
     }
     // cout << "nMultiGem1 = " << m_nEgem1 << endl;
     for(int i=0; i<m_nEgem1; i++) samplingTransfer(layer,1, i);
     // cout << "nMultiGem2 = " << m_nEgem2 << endl;
     for(int i=0; i<m_nEgem2; i++) samplingTransfer(layer,2, i);
     // cout << "nMultiGem3 = " << m_nEgem3 << endl;
     calMultiE(layer);
}
 
void SamplingGar::clear(){
     m_nIonE = 0;
     m_vecIonR.clear();
     m_vecIonPhi.clear();
     m_vecIonZ.clear();
     m_vecIonT.clear();
 
     m_nEgem1 = 0;
     m_idIon.clear();
     m_vecGem1dX.clear();
     m_vecGem1dZ.clear();
     m_vecGem1dT.clear();
 
     m_nEgem2 = 0;
     m_idGem1.clear();
     m_vecGem2dX.clear();
     m_vecGem2dZ.clear();
     m_vecGem2dT.clear();
 
     m_nEgem3 = 0;
     m_idGem2.clear();
     m_vecGem3dX.clear();
     m_vecGem3dZ.clear();
     m_vecGem3dT.clear();
 
     m_nMulElec = 0;
     m_eX.clear();
     m_eY.clear();
     m_eZ.clear();
     m_eT.clear();
}
 
bool SamplingGar::readSmearPar(){
     string filePath = getenv("CGEMDIGITIZERSVCROOT");
     string fileName;
     if(m_magConfig<0.05) fileName = filePath + "/dat/par_SamplingGar_0T.txt";
     else fileName = filePath + "/dat/par_SamplingGar_1T.txt";
     ifstream fin(fileName.c_str(), ios::in);
     cout << "fileName: " << fileName << endl;
 
     string strline;
     string strpar;
     if( ! fin.is_open() ){
      cout << "ERROR: can not open file " << fileName << endl;
      return false;
     }
     while(fin >> strpar){
      if('#' == strpar[0]){
           getline(fin, strline);
      } else if("Mean_X_function" == strpar){
           fin >> m_meanXpar[0] >> m_meanXpar[1];
      } else if("Sigma_X_function" == strpar){
           fin >> m_sigmaXpar[0] >> m_sigmaXpar[1] >> m_sigmaXpar[2] >> m_sigmaXpar[3];
      } else if("Mean_Z_function" == strpar){
           fin >> m_meanZpar;
      } else if("Sigma_Z_function" == strpar){
           fin >> m_sigmaZpar[0] >> m_sigmaZpar[1] >> m_sigmaZpar[2] >> m_sigmaZpar[3];
      } else if("Mean_T_function" == strpar){
           fin >> m_meanTpar[0] >> m_meanTpar[1];
      } else if("Sigma_T_function" == strpar){
           fin >> m_sigmaTpar[0] >> m_sigmaTpar[1] >> m_sigmaTpar[2] >> m_sigmaTpar[3];
      } else if("Transparency_Gem1" == strpar){
           fin >> m_transparency_gem1;
           m_transparency_gem1=1-(1-m_transparency_gem1)*m_TransMultiplier;
      } else if("Gain_Gem1" == strpar){
           fin >> m_gain_gem[0][0] >> m_gain_gem[0][1] >> m_gain_gem[0][2];
           //m_gain_gem[0][1]=m_gain_gem[0][1]*m_GainMultiplier;
 
      } else if("MeanX_Transf1" == strpar){
           fin >> m_meanX_transf[0];
      } else if("SigmaX_Transf1" == strpar){
           fin >> m_sigmaX_transf[0];
      } else if("MeanZ_Transf1" == strpar){
           fin >> m_meanZ_transf[0];
      } else if("SigmaZ_Transf1" == strpar){
           fin >> m_sigmaZ_transf[0];
      } else if("MeanT_Transf1" == strpar){
           fin >> m_meanT_transf[0];
      } else if("SigmaT_Transf1" == strpar){
           fin >> m_sigmaT_transf[0];
      } else if("Transparency_Gem2" == strpar){
           fin >> m_transparency_gem2;
           m_transparency_gem2=1-(1-m_transparency_gem2)*m_TransMultiplier;
      } else if("Gain_Gem2" == strpar){
           fin >> m_gain_gem[1][0] >> m_gain_gem[1][1] >> m_gain_gem[1][2];
           //m_gain_gem[1][1]=m_gain_gem[1][1]*m_GainMultiplier;
 
      } else if("MeanX_Transf2" == strpar){
           fin >> m_meanX_transf[1];
      } else if("SigmaX_Transf2" == strpar){
           fin >> m_sigmaX_transf[1];
      } else if("MeanZ_Transf2" == strpar){
           fin >> m_meanZ_transf[1];
      } else if("SigmaZ_Transf2" == strpar){
           fin >> m_sigmaZ_transf[1];
      } else if("MeanT_Transf2" == strpar){
           fin >> m_meanT_transf[1];
      } else if("SigmaT_Transf2" == strpar){
           fin >> m_sigmaT_transf[1];
      } else if("Transparency_Gem3" == strpar){
           fin >> m_transparency_gem3;
           m_transparency_gem3=1-(1-m_transparency_gem3)*m_TransMultiplier;
      } else if("Gain_Gem3" == strpar){
           fin >> m_gain_gem[2][0] >> m_gain_gem[2][1] >> m_gain_gem[2][2];
           //m_gain_gem[2][1]=m_gain_gem[2][1]*m_GainMultiplier;
 
      } else if("MeanX_Induc" == strpar){
           fin >> m_meanX_induc;
      } else if("SigmaX_Induc" == strpar){
           fin >> m_sigmaX_induc;
      } else if("MeanZ_Induc" == strpar){
           fin >> m_meanZ_induc;
      } else if("SigmaZ_Induc" == strpar){
           fin >> m_sigmaZ_induc;
      } else if("MeanT_Induc" == strpar){
           fin >> m_meanT_induc;
      } else if("SigmaT_Induc" == strpar){
           fin >> m_sigmaT_induc;
      }
     }
     fin.close();
     return true;
}
 
bool SamplingGar::samplingDrift(int layer, double driftD, int idIon){
     G4double frandom = G4UniformRand();
     // cout << "frandom = " << frandom << ",  trans = " << m_transparency_gem1 << endl;
     if(frandom > m_transparency_gem1) return false;
     int n = samplingGain(layer,0);
     double meanPhi = m_meanXpar[0] + m_meanXpar[1]*driftD;
     double sigmaPhi = m_sigmaXpar[0] + m_sigmaXpar[1]*driftD + m_sigmaXpar[2]*driftD*driftD + m_sigmaXpar[3]*driftD*driftD*driftD;
     double meanZ = 0.0;
     double sigmaZ = m_sigmaZpar[0] + m_sigmaZpar[1]*driftD + m_sigmaZpar[2]*driftD*driftD + m_sigmaZpar[3]*driftD*driftD*driftD;
     double meanT = m_meanTpar[0] + m_meanTpar[1]*driftD;
     double sigmaT = m_sigmaTpar[0] + m_sigmaTpar[1]*driftD + m_sigmaTpar[2]*driftD*driftD + m_sigmaTpar[3]*driftD*driftD*driftD;
     // cout << "debug in samplingDrift: " << setw(15) << driftD << setw(15) << meanPhi << setw(15) << sigmaPhi
     //       << setw(15) << sigmaZ << setw(15) << meanT << setw(15) << sigmaT << endl;
     for(int i=0; i<n; i++){
          double dphi = G4RandGauss::shoot(meanPhi, sigmaPhi*m_DiffuMultiplier);
      double dz = G4RandGauss::shoot(meanZ, sigmaZ*m_DiffuMultiplier);
      double dt = G4RandGauss::shoot(meanT, sigmaT);
 
      m_idIon.push_back(idIon);
      m_vecGem1dX.push_back(-1.0*dphi);
      m_vecGem1dZ.push_back(dz);
      m_vecGem1dT.push_back(dt);
     }
     m_nEgem1 += n;
     return true;
}
 
bool SamplingGar::samplingTransfer(int layer, int region, int idLastStep){
     double transp;
     if(1==region) transp = m_transparency_gem2;
     else transp = m_transparency_gem3;
     G4double frandom = G4UniformRand();
     if(frandom > transp) return false;
     int n = samplingGain(layer,region);
 
     for(int i=0; i<n; i++){
          double dphi = G4RandGauss::shoot(m_meanX_transf[region-1], m_sigmaX_transf[region-1]*m_DiffuMultiplier);
      double dz = G4RandGauss::shoot(m_meanZ_transf[region-1], m_sigmaZ_transf[region-1]*m_DiffuMultiplier);
      double dt = G4RandGauss::shoot(m_meanT_transf[region-1], m_sigmaT_transf[region-1]);
 
      if(1==region){
           m_idGem1.push_back(idLastStep);
           m_vecGem2dX.push_back(-1.0*dphi);
           m_vecGem2dZ.push_back(dz);
           m_vecGem2dT.push_back(dt);
      } else{
           m_idGem2.push_back(idLastStep);
           m_vecGem3dX.push_back(-1.0*dphi);
           m_vecGem3dZ.push_back(dz);
           m_vecGem3dT.push_back(dt);
      }
     }
     if(1==region) m_nEgem2 += n;
     else m_nEgem3 += n;
     return true;
}
 
int SamplingGar::samplingGain(int layer ,int region){
     double xRand = m_funPolya[layer][region]->GetRandom(1, 500);
     int gain = (int)xRand;
     return gain;
}
 
bool SamplingGar::calMultiE(int layer){
     for(int i=0; i<m_nEgem3; i++){
      //double dphiInduc = -1.0 * G4RandGauss::shoot(m_meanX_induc, m_sigmaX_induc);
      //double dzInduc = G4RandGauss::shoot(m_meanZ_induc, m_sigmaZ_induc);
      //double dtInduc = G4RandGauss::shoot(m_meanT_induc, m_sigmaT_induc);
 
      double dphiGem3 = m_vecGem3dX[i];
      double dzGem3 = m_vecGem3dZ[i];
      double dtGem3 = m_vecGem3dT[i];
 
      int idGem2 = m_idGem2[i];
      double dphiGem2 = m_vecGem2dX[idGem2];
      double dzGem2 = m_vecGem2dZ[idGem2];
      double dtGem2 = m_vecGem2dT[idGem2];
 
      int idGem1 = m_idGem1[idGem2];
      double dphiGem1 = m_vecGem1dX[idGem1];
      double dzGem1 = m_vecGem1dZ[idGem1];
      double dtGem1 = m_vecGem1dT[idGem1];
 
      int idIon = m_idIon[idGem1];
      double rini = m_vecIonR[idIon];
      double phiini = m_vecIonPhi[idIon];
      double zini = m_vecIonZ[idIon];
      double tini = m_vecIonT[idIon];
 
      double rNew = m_geomSvc->getCgemLayer(layer)->getInnerROfGapI();
      double dphiTot = dphiGem1 + dphiGem2 + dphiGem3;
      double xNew = rNew * phiini + dphiTot;
      double phiOut;    // phi angle
      if(xNew < 0){
           phiOut = xNew/rNew + CLHEP::twopi;
      } else if(xNew > (rNew * CLHEP::twopi)){
           phiOut = xNew/rNew - CLHEP::twopi;
      } else{
           phiOut = xNew / rNew;
      }
 
      m_eX.push_back(rNew * cos(phiOut));
      m_eY.push_back(rNew * sin(phiOut));
      m_eZ.push_back(zini + dzGem1 + dzGem2 + dzGem3);
      m_eT.push_back(tini + dtGem1 + dtGem2 + dtGem3);
      m_nMulElec++;
 
      // double rNew = m_geomSvc->getCgemLayer(layer)->getInnerROfAnode();
      // double dphiTot = dphiGem1 + dphiGem2 + dphiGem3 + dphiInduc;
      // double xNew = rNew * phiini + dphiTot;
      // double phiOut; // phi angle
      // if(xNew < 0){
      //      phiOut = xNew/rNew + CLHEP::twopi;
      // } else if(xNew > (rNew * CLHEP::twopi)){
      //      phiOut = xNew/rNew - CLHEP::twopi;
      // } else{
      //      phiOut = xNew / rNew;
      // }
 
      // m_eX.push_back(rNew * cos(phiOut));
      // m_eY.push_back(rNew * sin(phiOut));
      // m_eZ.push_back(zini + dzGem1 + dzGem2 + dzGem3 + dzInduc);
      // m_eT.push_back(tini + dtGem1 + dtGem2 + dtGem3 + dtInduc);
      // m_nMulElec++;
     }
     return true;
}