InductionGTS.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/InductionGTS.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 "TMath.h"
#include "TRandom.h"
#include "TCanvas.h"
 
#include <iomanip>
#include <iostream>
#include <fstream>
#include <cmath>
 
using namespace std;
 
InductionGTS::InductionGTS(){
}
 
InductionGTS::~InductionGTS()
{
 
  if(m_testing_ind == true) {
        output->Write();
    output->Close();
  }
 
  double tminbin = hcollected_charge->FindBin(apv_tstart);
  //  for(int it = tminbin; it < hnbin; it++) delete f[it];
  delete hintegratore;
  delete hcollected_charge;
  delete hcharge;
  delete f_FD;
 
 
}
 
void InductionGTS::init(ICgemGeomSvc* geomSvc, double magConfig){
  m_geomSvc = geomSvc;
  m_magConfig = magConfig;
  evt = 0;
  // 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 << "InductionGTS::init ERROR: can not open file " << filePath << endl;
  }
  else  std::cout << "InductionGTS::init: open file " << fileName << endl;
 
 
  while(fin >> strpar){
 
    if(strpar[0] == '#'){
      getline(fin, strline);
    }
    else if(strpar == "magnetic_field") {
      fin >> m_field;
    }
    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 << "InductionGTS::init: drift and avalanche parameters" << endl;
  cout << "field        " << m_field <<endl;
  cout << "tuning factors for diffusion " << m_tuning_factor_diff_perp << " (for orthogonal) " << m_tuning_factor_diff_paral << " (for parallel)" << endl;
 
  // diffusion orthogonal to mag field
  if(m_field) fileName = filePath + "/dat/par_ArIso_1T_GTS.txt";
  else        fileName = filePath + "/dat/par_ArIso_0T_GTS.txt";
 
 
  if(m_testing_ind == true) {
    output = new TFile("induction.root", "RECREATE");
    tree = new TNtuple("tree", "tree", "dphi:darc");
    tree_strip = new TNtuple("tree_strip", "tree for strip", "evt:id:view:qind:q:t");
  }
 
  // diffusion parallel to mag field
  string fileName2 = filePath + "/dat/par_ArIso_0T_GTS.txt";
  hcollected_charge = new TH1F("hcollected_charge", "hcollected_charge", hnbin, hxmin, hxmax);
  hintegratore = new TH1F("hintegratore", "", hnbin, hxmin, hxmax);
  hcharge = new TH1F("hcharge", "ASIC measured charge", 28, hxmin, hxmax);
  double tminbin = hcollected_charge->FindBin(apv_tstart);
  double dt = hcollected_charge->GetXaxis()->GetBinWidth(1);
 
  for(int it = tminbin; it < hnbin; it++) {
    double t = apv_tstart + it * dt;
    TString name = "f"; name += it;
    f[it] = new TF1(name, "[0] * ((x - [1]) / [2]) * TMath::Exp(-(x - [1]) / [2])", t, hxmax);
  }
 
  f_FD = new TF1("f_FD", "[0] + [1]/(1+TMath::Exp(-(x - [2])/[3]))", hxmin, hxmax);
 
 
 
  diffusion = new DiffusionGTS(m_geomSvc);
  diffusion->readGasPerpParameters(fileName);
  diffusion->readGasParalParameters(fileName2);
}
 
void InductionGTS::clear(){
 
  // temp                                                                                                                                   
  stripid_x.clear();                                                                                                             
  sheetid_x.clear();                                                                                                             
  charge_x.clear();                                                                                                             
  time_x.clear();                                                                                                               
  stripid_v.clear();                                                                                                             
  sheetid_v.clear();                                                                                                             
  charge_v.clear();                                                                                                             
  time_v.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();
}
 
 
void InductionGTS::setMultiElectrons(int layer, int nElectrons, std::vector<double> x, std::vector<double> y, std::vector<double> z, std::vector<double> t){
  double time_spent = 0.;
  clock_t begin = clock();
 
  clear();
  
  std::vector < double > stripid_x;
  std::vector < double > ind_charge_x;
  std::vector < double > ind_time_x;
 
  std::vector < double > stripid_v;
  std::vector < double > ind_charge_v;
  std::vector < double > ind_time_v;
 
  //  cout << "NELETTRONI " << nElectrons << endl;
  // ---------------------------------------------
  // loop on ions in avalanche
  for(int iele = 0; iele < nElectrons; iele++){
    
    double x_start = x.at(iele);
    double y_start = y.at(iele);
    double z_start = z.at(iele);
    double t_start = t.at(iele); // CHECK time here already contains the trigger time or not?? must be track_t0 + t_drift + t_diffusion
 
    // diffusion is local
    double x_end, y_end, z_end, t_end;
    bool isanode = drive_to_anode(layer, x_start, y_start, z_start, t_start, x_end, y_end, z_end, t_end);
 
    int sheet = 0; // CHECK THIS!!
    double phi = TMath::ATan2(y_start, x_start);
    if(layer > 0 && phi > TMath::Pi()) sheet = 1;
 
    CgemGeoReadoutPlane* anode = m_geomSvc->getReadoutPlane(layer, sheet);
    int X_ID;
    int V_ID;
    G4ThreeVector pos(x_end, y_end, z_end);
 
    anode->getStripID(pos, X_ID, V_ID); //  get the closest XID/VID
 
    // cout << endl;
    // cout << "start "<< x_start << " " << y_start << " " << z_start << " " << t_start << endl;
    // cout << "end "<< x_end << " " << y_end << " " << z_end << " " << t_end << endl;
    // cout << "X_ID " << X_ID << " V_ID " << V_ID << endl;
 
    // cout << "->radius " << TMath::Sqrt(x_end*x_end + y_end*y_end) 
    //   << " arco " << TMath::ATan2(y_end,x_end) * TMath::Sqrt(x_end*x_end + y_end*y_end)
    //   << " X_ID " << X_ID << endl;
 
 
    // charge sharing // CHECK TO BE CHANGED
    double percX = 0.5;
    double percV = 0.5;
    double qele = 1.6e-4;
 
    // save all the vectors in parallel
    stripid_x.push_back(X_ID);
    ind_charge_x.push_back(percX * qele);
    ind_time_x.push_back(t_end);
    sheetid_x.push_back(sheet);
    stripid_v.push_back(V_ID);
    ind_charge_v.push_back(percV * qele);
    ind_time_v.push_back(t_end);
    sheetid_v.push_back(sheet);
 
  }
 
  // final output
  
  // count and fill X stripIDs -----------------------------------------------------------
  for (std::vector< double >::iterator it = stripid_x.begin() ; it != stripid_x.end(); it++) {
    double stripid = *it;
    int ipos = it - stripid_x.begin();
    double sheetid = sheetid_x.at(ipos);
 
    std::vector< double >::iterator it2 = std::find(m_XstripID.begin(), m_XstripID.end(), stripid);
    std::vector< double >::iterator it3 = std::find(m_XstripSheet.begin(), m_XstripSheet.end(), sheetid);
 
    if(it2 == m_XstripID.end() || it3 == m_XstripSheet.end()) { // if the strip is not present, then add stripid and sheetid
      m_XstripID.push_back(stripid);
      m_XstripSheet.push_back(sheetid);
    } 
    else {                        // if the strip is present, then check if it belongs to the same sheetid
      double sheetid3 = *it3;
      if(sheetid != sheetid3) { // if the sheetid is not present, then add stripid and sheetid 
     m_XstripID.push_back(stripid);
     m_XstripSheet.push_back(sheetid);
      }
    }
  }
  m_NXstrips = m_XstripID.size();
  // cout << "nstrip X " << m_NXstrips << endl;
  
  // measure charge, after electronics, now APV-25 ---------------------------------------- 
  for(int istrip = 0; istrip < m_NXstrips; istrip++) {
    double xID = m_XstripID.at(istrip);  
    int view = 0;
    double strip_charge;
    double strip_time;
    double strip_dtime; // CHECK this!  
    bool isASIC = useAPV(xID, view, stripid_x, ind_charge_x, ind_time_x, strip_charge, strip_time, strip_dtime);
    m_XstripQ.push_back(strip_charge); 
    m_XstripT.push_back(strip_time);
  }
 
 
  // count and fill V stripIDs -----------------------------------------------------------
  for (std::vector< double >::iterator it = stripid_v.begin() ; it != stripid_v.end(); it++) {
    double stripid = *it;
    int ipos = it - stripid_v.begin();
    double sheetid = sheetid_v.at(ipos);
 
    std::vector< double >::iterator it2 = std::find(m_VstripID.begin(), m_VstripID.end(), stripid);
    std::vector< double >::iterator it3 = std::find(m_VstripSheet.begin(), m_VstripSheet.end(), sheetid);
 
    if(it2 == m_VstripID.end() || it3 == m_VstripSheet.end()) { // if the strip is not present, then add stripid and sheetid
      m_VstripID.push_back(stripid);
      m_VstripSheet.push_back(sheetid);
    } 
    else {                        // if the strip is present, then check if it belongs to the same sheetid
      double sheetid3 = *it3;
      if(sheetid != sheetid3) { // if the sheetid is not present, then add stripid and sheetid 
     m_VstripID.push_back(stripid);
     m_VstripSheet.push_back(sheetid);
      }
    }
  }
  m_NVstrips = m_VstripID.size();
  //  cout << "nstripv " << m_NVstrips << endl;
 
 
  // count and fill V stripIDs -----------------------------------------------------------
  for(int istrip = 0; istrip < m_NVstrips; istrip++) {
    double vID = m_VstripID.at(istrip);
    int view = 1;  
    //  cout << "V stripid " << vID << " sheetid " << m_VstripSheet.at(istrip) << endl;
    double strip_charge;
    double strip_time;
    double strip_dtime; // CHECK this!
    bool isASIC = useAPV(vID, view, stripid_v, ind_charge_v, ind_time_v, strip_charge, strip_time, strip_dtime);
    m_VstripQ.push_back(strip_charge); 
    m_VstripT.push_back(strip_time);
  }
 
  clock_t end = clock();
  time_spent += (double)(end - begin) / CLOCKS_PER_SEC;
  cout << "InductionGTS::induzione " << time_spent << " seconds" << endl;
 
  evt++;
 
}
 
 
 
bool InductionGTS::useAPV(int stripid, int view, std::vector< double > stripidvec,  std::vector< double > indchargevec,  std::vector< double > indtimevec, double &charge, double &time, double &dtime)
{
  
  // cout << "--> x stripid " << stripid << " sheetid " << m_XstripSheet.at(istrip) << endl;
   
    //    cout << "APV for strip " << stripid << " no. " << istrip << "/" << m_NXstrips << endl;
    // ----------- APV25 FAST SIMULATION ------------------
  
    double tau_APV = 50; // ns
    int ntime = 27; // time bin
    double tstep = 25; // ns
 
    hcollected_charge->Reset();
    TString name;
    if(m_testing_ind) {
      name = "hcollected_charge";
      name += stripid;
      hcollected_charge->SetName(name);
    }
 
   std::vector< double >::iterator iter = stripidvec.begin();
    while(1) {
      if(iter != stripidvec.end()+1) iter = std::find(iter, stripidvec.end(), stripid);
      else break;
      int ipos = iter - stripidvec.begin();
      std::vector< double >::iterator iter2 = indtimevec.begin() + ipos;
      std::vector< double >::iterator iter3 = indchargevec.begin() + ipos;
      hcollected_charge->Fill(*iter2, *iter3);
      //      cout << ipos << " time " << *iter2 << " charge " << *iter3 << endl;
      iter++;
    }
    double ind_charge = hcollected_charge->Integral();
    
    double dt = hcollected_charge->GetXaxis()->GetBinWidth(1);
    // integrator
    double tminbin = hcollected_charge->FindBin(apv_tstart);
    for(int it = tminbin; it < hnbin; it++) {
      double t = apv_tstart + it * dt;
      double timebin = hcollected_charge->FindBin(t);
      double integral = hcollected_charge->GetBinContent(timebin);
      
      name = "f"; name += it; name += "_strip"; name += stripid;
      f[it]->SetName(name);
      //      f[it] = new TF1(name, "[0] * ((x - [1]) / [2]) * TMath::Exp(-(x - [1]) / [2])", t, xmax);
      f[it]->SetParameter(0, TMath::Exp(1.)*integral);
      f[it]->SetParameter(1, t);
      f[it]->SetParameter(2, tau_APV);
    }
 
    hintegratore->Reset();
    if(m_testing_ind) {
      name = "hintegratore";
      name += stripid;
      hintegratore->SetName(name);
    }
 
    for(int it = tminbin; it < hnbin; it++) {
      double t = apv_tstart + it * dt;
      double qmeas = 0;
      for(int jt = tminbin; jt < hnbin; jt++) {
    double t2 = apv_tstart + jt * dt;
    if(t2 < t)  qmeas += f[jt]->Eval(t);
      }
      hintegratore->SetBinContent(it, qmeas);
    }
    
    hcharge->Reset();
    if(m_testing_ind) {
      name = "hcharge";
      name += stripid;
      hcharge->SetName(name);
    }
 
    double jitter = 0.; // CHECK APV
    for(int iapv = 0; iapv < ntime; iapv++) {
      double t = jitter + apv_tstart + iapv * tstep;
      double timebin = hintegratore->FindBin(t);
      double qmeas = hintegratore->GetBinContent(timebin);
      hcharge->AddBinContent(iapv+1, qmeas);
    }
 
    //    output->Write();
    // hcollected_charge->Write();
    // hintegratore->Write();
    // hcharge->Write();
    // -----------------
 
    // ------------FERMI DIRAC-------
    // CHARGE
    charge = hcharge->GetMaximum();
    //  cout << "charge " << charge << endl;  
    // TIME from FERMI DIRAC %%%%%%%%%%%%%%%%%%%%%%%%%%
    // start parameters                                
      double maxbin    = hcharge->GetMaximumBin();
      double QMaxHisto = hcharge->GetMaximum();
 
      //find the lower edge of the rising @ 10% of the maximum 
      double minbin = 0;
      for(int ibin = maxbin; ibin > 1; ibin--){
        if(hcharge->GetBinContent(ibin) < 0.1 * QMaxHisto) {
          minbin = ibin;
          break;
        }
      }
      double startvalues[5] = {0};
      double fitlimup[5]    = {0};
      double fitlimlow[5]   = {0};
 
      // average of the three bins before minbin  
      double meanfirstbin = (hcharge->GetBinContent(minbin-1) + hcharge->GetBinContent(minbin-2) + hcharge->GetBinContent(minbin-3))/3.;
      double sigma_MFB = TMath::Sqrt( ( pow(hcharge->GetBinContent(minbin-1) - meanfirstbin, 2) + pow(hcharge->GetBinContent(minbin-2) - meanfirstbin,2) + pow(hcharge->GetBinContent(minbin-3)- meanfirstbin,2)) / 3);
      if(sigma_MFB < TMath::Abs(meanfirstbin)) sigma_MFB = TMath::Abs(meanfirstbin);
 
      startvalues[0] = meanfirstbin;
      startvalues[1] = QMaxHisto;
      startvalues[2] = 0.5* (minbin + maxbin) * tstep;
      startvalues[3] = 0.5*tstep;
 
      // put limits
      fitlimlow[0] = startvalues[0] - 2 * sigma_MFB;
      fitlimlow[1] = startvalues[1] - 0.3 * startvalues[1];
      fitlimlow[2] = minbin * tstep;
      fitlimlow[3] = 0.1*tstep;
 
      fitlimup[0] = startvalues[0] + 2 * sigma_MFB;
      fitlimup[1] = startvalues[1] + 0.3 * startvalues[1];
      fitlimup[2] = maxbin * tstep;
      fitlimup[3] = 0.9*tstep;
 
      // ..................
      // the real Fermi Dirac
      //
      // FD(t) = K ----------------------------- + B 
      //            1 + exp( -(t - tFD)/sigFD ) 
      // [0] = B
      // [1] = K   
      // [2] = tFD
      // [3] = sigFD
      double minFD = minbin - 4;
      double maxFD = maxbin + 0;
      // cout<<"FERMI DIRAC FIT: \n"; 
      
      f_FD->SetRange(minFD*tstep, maxFD*tstep);
      name = "f_FD";
      name += stripid; 
      f_FD->SetName(name);
      
      f_FD->SetParameters(startvalues[0], startvalues[1], startvalues[2], startvalues[3]);
      f_FD->SetParLimits(0, fitlimlow[0], fitlimup[0]);
      f_FD->SetParLimits(1, fitlimlow[1], fitlimup[1]);
      f_FD->SetParLimits(2, fitlimlow[2], fitlimup[2]);
      f_FD->SetParLimits(3, fitlimlow[3], fitlimup[3]);
  
      hcharge->Fit(f_FD,"WQRB0");  // quiet/range/params
      time = 0;
      time = gRandom->Gaus(f_FD->GetParameter(2), 8); //CHECK ??
      // time = f_FD->GetParameter(2); // CHECK??
      dtime  = 0.;
      dtime = f_FD->GetParError(2);
      
      if(m_testing_ind) tree_strip->Fill(evt, stripid, view, ind_charge, charge, time);
    
 
    //-------------
}
 
bool InductionGTS::drive_to_anode(int layer, double xi, double yi, double zi, double ti, double &xf, double &yf, double &zf, double &tf) {
  
  double R = TMath::Sqrt(xi*xi + yi*yi);
  double phi = TMath::ATan2(yi, xi);
 
  double shift_x_induct = diffusion->shift_x_induct();
  double sigma_x_induct = diffusion->sigma_x_induct();
 
  double gem3_rad_out = m_geomSvc->getCgemLayer(0)->getCgemFoil(2)->getOuterROfCgemFoil();
  double shift_phi = shift_x_induct/gem3_rad_out;
  double sigma_phi = sigma_x_induct/gem3_rad_out;
  double phif = phi + gRandom->Gaus(shift_phi, sigma_phi);
 
  double anode_rad_in = m_geomSvc->getCgemLayer(layer)->getInnerROfAnode();
  xf = anode_rad_in * TMath::Cos(phif);
  yf = anode_rad_in * TMath::Sin(phif);
 
  // y
  double shift_y_induct = diffusion->shift_y_induct();
  double sigma_y_induct = diffusion->sigma_y_induct();
  zf = zi + gRandom->Gaus(shift_y_induct, sigma_y_induct);
 
  // cout << "xi " << xi << " yi " << yi << " zi " << zi << endl;
  // cout << "R " << R << " phi " << phi << " phif " << phif << endl;
  // cout << "xf " << xf << " yf " << yf << " zf " << zf << endl;
 
  // t                      
  double shift_t_induct = diffusion->shift_t_induct();
  double sigma_t_induct = diffusion->sigma_t_induct();
  tf = ti + gRandom->Gaus(shift_t_induct, sigma_t_induct);
 
  double dphi = phif - phi;
  double darc = anode_rad_in * dphi;
 
 
  if(m_testing_ind) tree->Fill(dphi, darc);
 
  return true;
}
 
 
/* output info of fired strips */
int    InductionGTS::getNXstrips()         const { return m_NXstrips; }
int    InductionGTS::getNVstrips()         const { return m_NVstrips; }
int    InductionGTS::getXstripSheet(int n) const { return m_XstripSheet.at(n); }
int    InductionGTS::getXstripID(int n)    const { return m_XstripID.at(n); }
int    InductionGTS::getVstripSheet(int n) const { return m_VstripSheet.at(n); }
int    InductionGTS::getVstripID(int n)    const { return m_VstripID.at(n); }
double InductionGTS::getXstripQ(int n)     const { return m_XstripQ.at(n); }
double InductionGTS::getVstripQ(int n)     const { return m_VstripQ.at(n); }
double InductionGTS::getXstripT(int n)     const { return m_XstripT.at(n); }
double InductionGTS::getVstripT(int n)     const { return m_VstripT.at(n); }
 
//  LocalWords:  tstep