SamplingGTS.cxx

Go to the documentation of this file.

// implementation:
//          R. Farinelli (University of Ferrara & INFN of Ferrara) 
//          L. Lavezzi   (IHEP & INFN of Torino)   
#include "CgemDigitizerSvc/SamplingGTS.h"
 
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/DataSvc.h"
 
#include "G4DigiManager.hh"
#include "Randomize.hh"
#include "G4ios.hh"
 
#include "CLHEP/Units/PhysicalConstants.h"
 
#include <iomanip>
#include <iostream>
#include <fstream>
#include <sstream>
#include <cmath>
#include <time.h>
 
#include "TFile.h"
#include "TRandom.h"
#include "TMath.h"
#include "TVector3.h"
 
 
using namespace std;
 
SamplingGTS::SamplingGTS(){
}
 
SamplingGTS::~SamplingGTS(){
  if(m_testing_sam == true) {
    output->Write();
    output->Close();
  }
}
 
void SamplingGTS::init(ICgemGeomSvc* geomSvc, double magConfig){
  // std::cout << "SamplingGTS::init" << std::endl;
  double time_spent = 0.;
  clock_t begin = clock();
  m_geomSvc = geomSvc;
  m_magConfig = magConfig;
 
     // SETTINGS -------------------------------------------------
     string filePath = getenv("CGEMDIGITIZERSVCROOT");
     string fileName;
 
     fileName = filePath + "/dat/par_settings_GTS.txt";
     ifstream fin(fileName.c_str(), ios::in);
 
     string strline;
     string strpar;
     if( ! fin.is_open() ){
       std::cout << "SamplingGTS::init ERROR: can not open file " << filePath << endl;
     }
     else  std::cout << "SamplingGTS::init: open file " << fileName << endl;
 
 
     while(fin >> strpar){
 
       if(strpar[0] == '#'){
     getline(fin, strline);
       }
       else if(strpar == "high_voltage") {
     fin >> m_hv;
       }
       else if(strpar == "magnetic_field") {
     fin >> m_field;
       }
       else if(strpar == "tuning_factor_gain") {
     fin >> m_tuning_factor_gain;
       }
       else if(strpar == "tuning_factor_diff_perp") {
     fin >> m_tuning_factor_diff_perp;
       }
       else if(strpar == "tuning_factor_diff_paral") {
     fin >> m_tuning_factor_diff_paral;
       }
     }
 
     cout << "SamplingGTS::init: drift and avalanche parameters" << endl;
     cout << "high voltage " << m_hv << "V" << endl;
     cout << "field        " << m_field <<endl;
     cout << "tuning factors:" << endl;
     cout << "1) for gain " << m_tuning_factor_gain << endl;
     cout << "2) for diffusion " << m_tuning_factor_diff_perp << " (for orthogonal) " << m_tuning_factor_diff_paral << " (for parallel)" << endl;
 
 
     // GAIN
     bool gempar = readGemParameters();
     if(!gempar) return; // CHECK ??  
   
     // DIFFUSION & MAG FIELD
     bool gaspar = readGasParameters();
     if(!gaspar) return; // CHECK ??  
 
     m_nMulElec = 0;     
 
     if(m_testing_sam == true) {
       output = new TFile("gem.root", "RECREATE");
       tree = new TNtuple("tree", "tree", "dphi:darc:dx:dy:dz:z_local");
     }
 
 
     clock_t end = clock();
     time_spent += (double)(end - begin) / CLOCKS_PER_SEC;
     cout << "Sampling::init " << time_spent << " seconds" << endl;
 
 
}
 
void SamplingGTS::setIonElectrons(int layer, int nElectrons, std::vector<double> x, std::vector<double> y, std::vector<double> z, std::vector<double> t){
  cout << "SamplingGTS::setIonElectrons " << nElectrons << endl;
  clear();
 
  double time_spent = 0.;
  clock_t begin = clock();
 
  //  cout << "drift setIonElectrons " << nElectrons << endl;
  // cout << endl;
  // loop on ions from ionization
  for(int iele = 0; iele < nElectrons; iele++){
    //  cout << "iele " << iele << endl;
    // does the electron survive to GEM1?
    if(!is_survived()) continue;
 
    // in MARS: 
    // 1) the z_local of the GEM is the radial direction
    // --> the x_local shift and smearing are on the arc
    // 2) the y_local is the direction parallel to the mag field
    // --> the y_local shift and smearing are on z_global
 
    // generate effective gain from POLYA
    int tot_gain = compute_gain();
 
    //  double time_spent = 0.;
    // clock_t begin = clock();
  
    // cout << "tot_gain " << tot_gain << endl;
    for(int jele = 0; jele < tot_gain; jele++) {
      double xf, yf, zf, tf;
 
      compute_diffusion_on_GEM3(x[iele], y[iele], z[iele], t[iele], xf, yf, zf, tf);
      // cout << "jele " << jele << endl;
 
      // set output
      m_eX.push_back(xf);
      m_eY.push_back(yf);
      m_eZ.push_back(zf);
      m_eT.push_back(tf);
      m_nMulElec++;
      /**
     cout << iele << " " << jele << " diffused x " << x[iele] << " to " << xf << endl;
     cout << iele << " " << jele << " diffused y " << y[iele] << " to  " << yf << endl;
     cout << iele << " " << jele << " diffused z " << z[iele] << " to " << zf << endl;
     cout << iele << " " << jele << " diffused t " << t[iele] << " to " << tf << endl;
      **/
 
    }
  }
  clock_t end = clock();
  time_spent += (double)(end - begin) / CLOCKS_PER_SEC;
  cout << "Sampling::diffusion " << m_nMulElec << " time " << time_spent << " seconds" << endl;
    
 
  
}
 
void SamplingGTS::clear(){
  m_nMulElec = 0;
  m_eX.clear();
  m_eY.clear();
  m_eZ.clear();
  m_eT.clear();
}
 
// ------------------------------------------------------- set parameters
bool SamplingGTS::readGasParameters() {
 
  string filePath = getenv("CGEMDIGITIZERSVCROOT");
  // diffusion orthogonal to mag field
  string fileName;
  if(m_field) fileName = filePath + "/dat/par_ArIso_1T_GTS.txt";
  else        fileName = filePath + "/dat/par_ArIso_0T_GTS.txt";
  // diffusion parallel to mag field
  string fileName2 = filePath + "/dat/par_ArIso_0T_GTS.txt";
 
  diffusion = new DiffusionGTS(m_geomSvc);
  diffusion->readGasPerpParameters(fileName);
  diffusion->readGasParalParameters(fileName2);
 
  return true;
 
}
 
bool SamplingGTS::readGemParameters() {
 
  // TRANSPARENCY -------------------------------------------------
  string filePath = getenv("CGEMDIGITIZERSVCROOT");
  string fileName;
 
  fileName = filePath + "/dat/par_ArIso_GEM_GTS.txt";
  ifstream fin(fileName.c_str(), ios::in);
 
  string strline;
  string strpar;
  if( ! fin.is_open() ){
    std::cout << "SamplingGTS::readGemParameters ERROR: can not open file " << filePath << endl;
    return false;
  }
  else  std::cout << "SamplingGTS::readGemParameters: open file " << fileName << endl;
 
  bool is_tra = false;
  while(fin >> strpar){
 
    if(strpar[0] == '#'){
      getline(fin, strline);
    } 
    else if( atoi(strpar.c_str()) == m_hv) {
      fin >> m_eff_col[0] >> m_eff_col[1] >> m_eff_col[2] >> m_eff_ext[0] >> m_eff_ext[1] >> m_eff_ext[2];
      is_tra = true;   
      break;
    }
  }
 
  cout << "SamplingGTS::readGemParameters: parameters of GEM @ HV = " << m_hv << " V" << endl;
  cout << " collection eff GEM1: " << m_eff_col[0] 
       << " collection eff GEM2: " << m_eff_col[1] 
       << " collection eff GEM3: " << m_eff_col[2] << endl;
  cout << " extraction eff GEM1: " << m_eff_ext[0] 
       << " extraction eff GEM2: " << m_eff_ext[1] 
       << " extraction eff GEM3: " << m_eff_ext[2] << endl;
  cout << " is transparent " << is_tra << endl;
 
  if(!is_tra) return false;
 
  // GAIN ------------------------------
  fileName = filePath + "/dat/CHECK_ArIso_EffGain_10M_GTS.root"; // to be substituted (this obtained by gain_c1.C)
  TFile *file_gain_eff = new TFile(fileName.c_str(), "READ");
  TString histoname = "h_gain_eff_0_"; histoname += m_hv;
  h_gain_eff = (TH1F*) file_gain_eff->Get(histoname);
  return true;
}
 
// ----------------------------------------------------------- transparency
bool SamplingGTS::is_survived() {
  return (gRandom->Uniform() < m_eff_col[0]);  // CHECK 14) TRandom, as in ionization
}
 
double SamplingGTS::compute_gain() {
  return m_tuning_factor_gain * h_gain_eff->GetRandom();
}
 
// ----------------------------------------------------------- drifting
void SamplingGTS::compute_diffusion_on_GEM3(double xi, double yi, double zi, double ti, double &xf, double &yf, double &zf, double &tf) {
 
 
  double R_catout[3] = {m_geomSvc->getCgemLayer(0)->getInnerROfCgemLayer() + m_geomSvc->getThicknessOfCathode(0),
                        m_geomSvc->getCgemLayer(1)->getInnerROfCgemLayer() + m_geomSvc->getThicknessOfCathode(1), 
                        m_geomSvc->getCgemLayer(2)->getInnerROfCgemLayer() + m_geomSvc->getThicknessOfCathode(2)}; 
 
  double thick[3] = {m_geomSvc->getThicknessOfGapD(0), 
                     m_geomSvc->getThicknessOfGapD(1), 
                     m_geomSvc->getThicknessOfGapD(2)}; 
 
  double R = TMath::Sqrt(xi*xi + yi*yi);
  // cout << "xi " << xi << " yi " << yi << " R " << R << endl;
  int ilayer = 0;
  if(R < R_catout[1])      ilayer = 0;
  else if(R < R_catout[2]) ilayer = 1; 
  else                     ilayer = 2; 
  
  // cout << "xi " << xi << " yi " << yi << " R " << R << " ilayer " << ilayer <<  endl; 
 
  double drift_thick = thick[ilayer];
  double cat_rad_out = R_catout[ilayer];
  double gem3_rad_in = m_geomSvc->getCgemLayer(0)->getCgemFoil(2)->getInnerROfCgemFoil();
 
  // cout << "drift_thick " << drift_thick << endl;
  // cout << "cat_rad_out " << cat_rad_out << endl;
  // cout << "gem3_rad_in " << gem3_rad_in << endl;
 
 
  // in MARS:
  // 1) the z_local of the GEM is the radial direction
  // --> the x_local shift and smearing are on the arc
  // 2) the y_local is the direction parallel to the mag field
  // --> the y_local shift and smearing are on z_global
 
  // 1)
  double phi = TMath::ATan2(yi, xi);
  // cout << "phi " << phi << "[rad] -> [deg]" << phi*TMath::RadToDeg() << endl;
 
  double z_local = R - cat_rad_out;
  if(z_local < 0 || z_local > drift_thick) cout << "ERRRRRRORE " << z_local << endl;
  // cout<< "z_local " << z_local << endl;
 
  double shift_x_drift  = diffusion->shift_x_drift(z_local);
  double shift_x_transf = diffusion->shift_x_transf();
  double sigma_x_drift  = diffusion->sigma_x_drift(z_local);
  double sigma_x_transf = diffusion->sigma_x_transf();
 
  double shift_x = shift_x_drift + 2 * shift_x_transf;
  double sigma_x = m_tuning_factor_diff_perp * TMath::Sqrt(pow(sigma_x_drift, 2) + 2 * pow(sigma_x_transf, 2));
  // cout << "m_tuning_factor_diff_perp " << m_tuning_factor_diff_perp << endl;
  // cout << diffusion->sigma_x_drift(z_local) << " " << diffusion->sigma_x_transf() << endl;
  // cout << "shift_x " << shift_x << " sigma_x " << sigma_x << endl;
 
  double shift_phi = shift_x/gem3_rad_in;
  double sigma_phi = sigma_x/gem3_rad_in;
  // cout << "shift_phi " << shift_phi << " sigma_phi " << sigma_phi << endl;                                                                   
  double phif = phi + gRandom->Gaus(shift_phi, sigma_phi);
  // cout << "phif " << phif << endl;
  xf = gem3_rad_in * TMath::Cos(phif);
  yf = gem3_rad_in * TMath::Sin(phif);
  // cout << "xf " << xf << " yf " << yf << endl;
 
  // 2)
  double shift_y = diffusion->shift_y_drift(z_local) + 2 * diffusion->shift_y_transf();
  double sigma_y = m_tuning_factor_diff_paral * TMath::Sqrt(pow(diffusion->sigma_y_drift(z_local), 2) + 2 * pow(diffusion->sigma_y_transf(), 2));
  zf = zi + gRandom->Gaus(shift_y, sigma_y);
 
  // cout<< "shift_y " << shift_y << "sigma_y " << sigma_y << endl;
  // cout << "zf " << zf << endl;
 
  // time
  double shift_t = diffusion->shift_t_drift(z_local) + 2 * diffusion->shift_t_transf();
  double sigma_t = TMath::Sqrt(pow(diffusion->sigma_t_drift(z_local), 2) + 2 * pow(diffusion->sigma_t_transf(), 2));
  tf = ti + gRandom->Gaus(shift_t, sigma_t);
 
  if(m_testing_sam==true) {
    
    //      cout << "z_local " << z_local << " phif " << phif << " xf " << xf << " yf " << yf << " zf " << zf << endl; 
    float dphi = phif-phi;
    float darc = 1000 * gem3_rad_in * dphi;
    float dx = 1000*(xf-xi);
    float dy = 1000*(yf-yi);
    float dz = 1000*(zf-zi);
 
    tree->Fill(dphi, darc, dx, dy, dz, z_local);
 
    // cout << phif-phi << " " << 100*(xf-xi) << " " << 100*(yf-yi) << " " << 100*(zf-zi) << endl;
    // cout << "delta x " << xf << " " << xi << endl;
    // cout << "delta y " << yf << " " << yi << endl;
    // cout << "delta z " << zf << " " << zi << endl;
  }
 
 
 
 
 
 
}