IonizationGTS.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/IonizationGTS.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 <iomanip>
#include <iostream>
#include <fstream>
#include <cmath>
 
#include "TRandom.h"
#include "TMath.h"
 
// using namespace std;
 
// mean number of clusters in 5 mm
#define N_PRI_AVE   18.0016 // CHECK 2) deve essere letto da un file di settaggi
 
// n_electron per cluster ArIso CHECK 2)
const double int_prob_electron_per_cluster[100]={0.902864,0.920309,0.93045,0.936283,0.940138,0.943808,0.948277,0.954485,0.96098,0.966479,
                         0.970936,0.974523,0.977442,0.979843,0.981959,0.983657,0.985152,0.986421,0.987404,0.988233,
                         0.990002,0.990783,0.991482,0.992107,0.992666,0.993168,0.99362,0.994028,0.994397,0.994731,
                         0.995036,0.995313,0.995567,0.9958,0.996015,0.996212,0.996395,0.996564,0.996721,0.996867,
                         0.997003,0.99713,0.997248,0.99736,0.997464,0.997562,0.997655,0.997742,0.997825,0.997903,
                         0.997976,0.998046,0.998113,0.998176,0.998236,0.998293,0.998348,0.9984,0.998449,0.998497,
                         0.998542,0.998586,0.998628,0.998668,0.998706,0.998743,0.998779,0.998813,0.998846,0.998877,
                         0.998908,0.998938,0.998966,0.998993,0.99902,0.999046,0.999071,0.999095,0.999118,0.99914,
                         0.999162,0.999183,0.999204,0.999224,0.999243,0.999262,0.99928,0.999298,0.999315,0.999332,
                         0.999348,0.999364,0.99938,0.999395,0.99941,0.999424,0.999438,0.999452,0.999465,0.999478};
 
 
 
IonizationGTS::IonizationGTS(){
}
 
IonizationGTS::~IonizationGTS(){
  if(m_testing == true) {
    output->Write();
    output->Close();
  }
}
 
void IonizationGTS::init(unsigned int random, ICgemGeomSvc* geomSvc, double magConfig){
  double time_spent = 0.;
  clock_t begin = clock();
  
  //  std::cout << "IonizationGTS::init" << std::endl;
  m_random = random;
  gRandom->SetSeed(m_random); // CHECK 6) TRandom 
  m_geomSvc = geomSvc;
  m_magConfig = magConfig;
  
  cout << "IonizationGTS::init: m_testing " << m_testing << endl;
  
  
  if(m_testing==true) {
    output = new TFile("ionization.root", "RECREATE");
    h_distance_cluster = new TH1F("h_distance_cluster", "distance between two consecutive clusters", 100, 0, 5);
    h_nof_cluster = new TH1F("h_nof_cluster", "number of clusters", 100, 0, 100);
  }
 
  clock_t end = clock();
  time_spent += (double)(end - begin) / CLOCKS_PER_SEC;
 }
 
// CHECK 7) not using particle now
// CHECK 8) not using momentum now
void IonizationGTS::setTrack(int particle, int charge, double p, double trkPosIn[], double trkPosOut[]) {
   
  double time_spent = 0.;
  clock_t begin = clock();
 
  clear();
  // cout << "IonizationGTS setTrack" << endl;
  
  // CHECK 1) 3D
  m_track_length_limit = TMath::Sqrt((trkPosOut[0] - trkPosIn[0])*(trkPosOut[0] - trkPosIn[0]) + (trkPosOut[1] - trkPosIn[1])*(trkPosOut[1] - trkPosIn[1]) + (trkPosOut[2] - trkPosIn[2])*(trkPosOut[2] - trkPosIn[2]));
  m_track_length = 0;
 
  m_n_ion_mm = N_PRI_AVE/m_track_length_limit; // CHECK 2) N_PRI_AVE
 
  /**
  cout << "m_n_ion_mm " << m_n_ion_mm << " N_PRI_AVE " << N_PRI_AVE << " m_track_length_limit " << m_track_length_limit << endl;
  cout << "xIN  " << trkPosIn[0]  << " yIN  " << trkPosIn[1]  << " zIN  " << trkPosIn[2] << endl;    
  cout << "xOUT " << trkPosOut[0] << " yOUT " << trkPosOut[1] << " zOUT " << trkPosOut[2] << endl;    
  **/
 
  // loop over ionization clusters
  int counter = 0;
 
  double last_x;
  double last_y;
  double last_z;
  double last_t;
 
  while(1) {
    
    // primary ionization
    bool is_inside = generate_primary_ele();
    if(!is_inside) break;
    
    // secondary ionization
    int nof_sec_ele = generate_secondary_ele();
    m_nIonE += nof_sec_ele;
    
    // compute position
    double x = 0;
    double y = 0; 
    double z = 0;
    double t = 0; // <-track_time;  // CHECK 3) add the track time 
    compute_pos(trkPosIn, trkPosOut, x, y, z);
    //    cout << "ionization electron pos " << x << " " << y << " " << z << endl;    
 
    // set cluster pos
    m_cx.push_back(x); // CHECK 5) cluster info
    m_cy.push_back(y); // CHECK 5) cluster info
    m_cz.push_back(z); // CHECK 5) cluster info
    m_ct.push_back(t); // CHECK 5) cluster info
 
    //    cout << "x " << x << " y " << y << " z " << z << " t " << t << endl;
    
    if(m_testing==true && counter > 0) {
      double distance = TMath::Sqrt((last_x - x)*(last_x - x) + (last_y - y)*(last_y - y) + (last_z - z)*(last_z - z));
      h_distance_cluster->Fill(distance);
      last_x = x;
      last_y = y;
      last_z = z;
      last_t = t;
    }
 
    // set ele pos
    for(int iele = 0; iele < nof_sec_ele + 1; iele++)
      {
    m_ex.push_back(x);
    m_ey.push_back(y);
    m_ez.push_back(z);
    m_et.push_back(t);
    
      }
    counter++;
    //    cout << "ion position " << x << " " << y << " " << z << endl;
 
  }
  
  if(m_nIonE != m_ex.size()) std::cout << "IonizationGTS::setTrack ERRRRRRRRRRRRRROR" << std::endl; // CHECK delete this
 
  if(m_testing==true) h_nof_cluster->Fill(counter-1);
 
  clock_t end = clock();
  time_spent += (double)(end - begin) / CLOCKS_PER_SEC;
  cout << "Sampling::ionization " << m_nIonE << " time " << time_spent << " seconds" << endl;
 
 
 
 
}
 
// computation of the ionization clusters:
// the number of clusters is Poissonian -> the inder-distance is neg-exponential
bool IonizationGTS::generate_primary_ele() {
  
  double u = gRandom->Uniform(0, 1); // CHECK 6) TRandom
  // step along the trajectory, that for now is a straight line [x0, z0] -> [x1, z1]
  double dl_extracted = -(1./m_n_ion_mm) * TMath::Log(1 - u);                                                                                
  m_track_length += dl_extracted;
 
  //  cout << "m_track_length " <<   m_track_length << " dl_extracted " << dl_extracted << endl;
 
 
  if(m_track_length > m_track_length_limit) return false;
  
  m_nIonC++;
  m_nIonE++;
  return true;
}
 
// computation of the number of secondary elecrons per primary ionization:
// they are all generated in the same position of the prim cluster
int IonizationGTS::generate_secondary_ele() {
  int n_elec=0;
  double prob = gRandom->Uniform(0, 1);  // CHECK 6) TRandom
 
  for(int iele = 0; iele < 100; iele++) { // CHECK 4) cut at 100
    n_elec = iele + 1;
    if(prob < int_prob_electron_per_cluster[iele]) break;
  }
 
  if(n_elec>100) { 
    std::cout << "IonizationGTS::generate_secondary_ele: increase the size above 100"<<endl; // CHECK 4)
    n_elec=100;
  }
 
  return n_elec;
}
 
void IonizationGTS::compute_pos(double posin[], double posout[], double &x, double &y, double &z) {
  
  double t = (m_track_length_limit - m_track_length)/m_track_length_limit;
  x = posin[0] * t + posout[0] * (1 - t);
  y = posin[1] * t + posout[1] * (1 - t); 
  z = posin[2] * t + posout[2] * (1 - t);
 
}
 
/**
// from global to local frame
// xl = coordinate along the cylinder circunference (parallel to the electrodes) = (yg/xg) * sqrt(xg**2 + yg**2)
// yl = coordinate along the cylinder axis (coincident with zg)
// zl = coordinate along the cylinder radius (orthogonal to the electrodes - [inner_radius_layer + cathode_thick]) = sqrt(xg**2 + yg**2) - [inner_radius_layer + cathode_thick])
void IonizationGTS::port_track_to_loc(double xg_in, double yg_in, double zg_in, double xg_out, double yg_out, double zg_out, 
                      double &xl_in, double &yl_in, double &zl_in, double &xl_out, double &yl_out, double &zl_out) 
{
 
  xl_in = 0;
  yl_in = 0;
  zl_in = 0;
 
  // ---------------------------------------------------
  double R_catout[3] = {m_geomSvc->getCgemLayer(0)->getInnerROfCgemLayer() + m_geomSvc->getThicknessOfCathode(), 
            m_geomSvc->getCgemLayer(1)->getInnerROfCgemLayer() + m_geomSvc->getThicknessOfCathode(),
            m_geomSvc->getCgemLayer(2)->getInnerROfCgemLayer() + m_geomSvc->getThicknessOfCathode()};
  
  double thick[3] = {m_geomSvc->getCgemLayer(0)->getThicknessOfGapD(), 
             m_geomSvc->getCgemLayer(1)->getThicknessOfGapD(),
             m_geomSvc->getCgemLayer(2)->getThicknessOfGapD()};
 
  TVector3 vg_in_xy(xg_in, yg_in, 0);
  TVector3 vg_out_xy(xg_out, yg_out, 0);
  // radial distance
  double radg_in = vg_in_xy.Perp();
  if(radg_in < R_catout[1])      zl_out = thick[0];
  else if(radg_in < R_catout[2]) zl_out = thick[1];
  else                           zl_out = thick[2];
 
  // direction xy
  TVector3 d_xy = vg_out_xy - vg_in_xy;
  // direction vin, with magnitute=thick
  TVector3 r_xy = vg_in_xy; r_xy.SetMag(zl_out);
  double alpha = TMath::ACos(d_xy.Dot(r_xy)/(d_xy.Mag()*r_xy.Mag()));
  double sign = d_xy.Cross(r_xy).Z()/fabs(d_xy.Cross(r_xy).Z());
  if(fabs(yg_out) < 1e-7) xl_out = 0;
  else xl_out = sign * zl_out * TMath::Tan(alpha);
 
  // much easier procedure for yl_out                                                                        
  double delta_z = zg_out - zg_in;
  yl_out = delta_z;
 
}
**/
 
void IonizationGTS::clear(){
  m_track_length_limit = 0;
  m_track_length = 0;
 // cluster
  m_nIonC = 0;
  m_cx.clear();
  m_cy.clear();
  m_cz.clear();
  m_ct.clear();
  // electron
  m_nIonE = 0;
  m_ex.clear();
  m_ey.clear();
  m_ez.clear();
  m_et.clear();
}