HarpoSimReadout.cxx

Go to the documentation of this file.

#include "HarpoSimReadout.h"
#include "HarpoSimElectronics.h"
#include "HarpoConfig.h"
#include "HarpoDebug.h"

#include "TMath.h"
#include "TRandom3.h"

ClassImp(HarpoSimReadout);



HarpoSimReadout::HarpoSimReadout()
{

  fStrips = (HarpoSimStrips*)NULL;//new HarpoSimStrips();

  //  fGain = 500;
  fGainMM = 20;
  fGainGEM1 = 5;
  fGainGEM2 = 5;
  fRatioYX = 1;
  fSigmaPrf = 0.35;

  fStripWidth = 1.;
  fPolya = 1.;

  fDigi = (HarpoSimElectronics*)NULL;



  Double_t ratioYX;
  if ( ! gHConfig->Lookup("Sim.RatioYX",ratioYX) )
    Info("Constructor","Use default GainMM %g",fRatioYX);
  else
    fRatioYX = ratioYX;

  Double_t gainMM;
  if ( ! gHConfig->Lookup("Sim.GainMM",gainMM) )
    Info("Constructor","Use default GainMM %g",fGainMM);
  else
    fGainMM = gainMM;

  Double_t gainGEM1;
  if ( ! gHConfig->Lookup("Sim.GainGEM1",gainGEM1) )
    Info("Constructor","Use default GainGEM1 %g",fGainGEM1);
  else
    fGainGEM1 = gainGEM1;

  Double_t gainGEM2;
  if ( ! gHConfig->Lookup("Sim.GainGEM2",gainGEM2) )
    Info("Constructor","Use default GainGEM2 %g",fGainGEM2);
  else
    fGainGEM2 = gainGEM2;

  Double_t polya;
  if ( ! gHConfig->Lookup("Sim.Polya",polya) )
    Info("Constructor","Use default Polya %g",fPolya);
  else
    fPolya = polya;

  Double_t gain;
  if ( ! gHConfig->Lookup("Sim.Gain",gain) )
    Info("Constructor","Use default Gain %g/%g/%g",fGainMM,fGainGEM1,fGainGEM2);
  else{
    fGainMM = gain;
    fGainGEM1 = 1;
    fGainGEM2 = 1;
  }
 
  Double_t sigmaPrf;
  if ( ! gHConfig->Lookup("Sim.SigmaPrf",sigmaPrf) )
    Info("Constructor","Use default SigmaPrf %g",fSigmaPrf);
  else
    fSigmaPrf = sigmaPrf;
 
  Double_t stripWidth;
  if ( ! gHConfig->Lookup("Sim.StripWidth",stripWidth) )
    Info("Constructor","Use default StripWidth %g",fStripWidth);
  else
    fStripWidth = stripWidth;
      

  Info("HarpoSimReadout","SigmaPrf: %g",fSigmaPrf);


}


HarpoSimReadout::~HarpoSimReadout()
{

  fStrips->Delete();

}

Int_t HarpoSimReadout::Reset()
{

  //  fStrips->Reset();
  // fStrips = new HarpoSimStrips();


  return 1;
}

HarpoSimStrips* HarpoSimReadout::ProcessElectron(Double_t x, Double_t y, Double_t t, hSimReadoutRegion region)
{


  // Double_t gain = 0;
  // if(region == hSimReadoutRegionDrift) gain = fGainGEM2*fGainGEM1*fGainMM;
  // if(region == hSimReadoutRegionTransfer2) gain = fGainGEM1*fGainMM;
  // if(region == hSimReadoutRegionTransfer1) gain = fGainMM;

  Double_t charge = 1;
  if(region == hSimReadoutRegionDrift){
    charge *= fGainGEM2;
    charge *= fGainGEM1;
    charge *= fGainMM;
  }
  if(region == hSimReadoutRegionTransfer2){
    charge *= fGainGEM1;
    charge *= fGainMM;
  }
  if(region == hSimReadoutRegionTransfer1){
    charge *= fGainMM;
  }

  // charge *= gRandom->Exp(1);
  // charge *= gRandom->Exp(1)*gRandom->Exp(1)gRandom->Exp(1);
  Double_t rndm = gRandom->Rndm();
  charge *= 4*(1.-fPolya)*rndm*rndm*rndm + gRandom->Exp(fPolya);


  if(charge<=0) return 0;
  //  Double_t charge = gain*gRandom->Exp(1);

  Double_t sigma = fSigmaPrf;

  //HarpoSimStrips* strips = new HarpoSimStrips();

  if(gHarpoDebug>4) Info("ProcessElectron","AddCharge(%f %f %f %f %f)",charge,sigma,x,y,t);


  if(isBlocked(x,y)) return 0;

  //fStrips->
  AddCharge(charge,sigma,x,y,t);

  return fStrips;
}




Int_t HarpoSimReadout::AddCharge(Double_t charge, Double_t sigma, Double_t x, Double_t y, Double_t t)
{

  //Info("AddCharge","%d %g %g %g %g %g",fNhitsX,charge,sigma,x,y,t);
  
  if(charge<0.1) return 0;

  // Find Strip number (X and Y coordinate)
  Int_t iX = Int_t((x + 150.)/fStripWidth);
  Int_t iY = Int_t((y + 150.)/fStripWidth);
  
  if(iX<0 || iX >= NCHAN) return 0;
  if(iY<0 || iY >= NCHAN) return 0;

  Double_t chargeSum = 0;

  Int_t nbins = Int_t(5*sigma/fStripWidth);
  for(Int_t i = iX - nbins; i<iX+nbins+1; i++){ // One layer (fStripWidth mm) in the X direction
    Double_t Xi1 = fStripWidth*i-150.;
    Double_t Xi2 = Xi1 + fStripWidth*0.379;
    Double_t Xi3 = Xi2 + fStripWidth*0.1;
    Double_t Xi4 = Xi3 + fStripWidth*0.421;
    Double_t chXstrip = 0, chXpads = 0;
    if(sigma>0){
      // Charge fraction chX on the strip number i (between x = Xi1 and x = Xi2)
      chXstrip = (TMath::Abs(TMath::Erf((x-Xi1)/sigma) - TMath::Erf((x-Xi2)/sigma)))/2;
      // Charge fraction chY on the neighbour pads  (between x = Xi3 and x = Xi4)
      chXpads = (TMath::Abs(TMath::Erf((x-Xi3)/sigma) - TMath::Erf((x-Xi4)/sigma)))/2;
    }else{ // punctual charge deposition (sigma = 0)
      if(x>Xi3 && x<Xi4) chXstrip = 1; 
      if(x>Xi1 && x<Xi2) chXpads = 1; 
    }
    chXstrip *= charge;
    chXpads *= charge;
    if(chXstrip>0.1){ // Add the charge to the X strip i
      if(fStrips != NULL) fStrips->AddHit(1,i,t,chXstrip);
      fDigi->AddCharge(XPLANE,i,t,chXstrip);
      chargeSum += chXstrip;
    }
    if(chXpads<0.1) continue;
    for(Int_t j = iY - nbins; j<iY+nbins+1; j++){ // One pad on neighbour row to the strip i
      Double_t Yj1 = fStripWidth*j-150.;
      Double_t Yj2 = Yj1 + fStripWidth*0.9;
      // charge fraction chY on the pad  (between x = Xi3 and x = Xi4 and y = Yj1 and y = Yj2)
      Double_t chY = 0;
      if(sigma>0) 
        chY = (TMath::Abs(TMath::Erf((y-Yj1)/sigma) - TMath::Erf((y-Yj2)/sigma)))/2;
      else{ // punctual charge deposition (sigma = 0)
        if(y>Yj1 && y<Yj2) chY = 1;
      }
      chY *= chXpads;
      chY *= fRatioYX;
      if(chY>0.1){ // Add charge to Y strip j
        if(fStrips != NULL) fStrips->AddHit(0,i,t,chY);
        fDigi->AddCharge(YPLANE,j,t,chY);
        chargeSum += chY;
      }
    }
      
  }
  
  if(gHarpoDebug>3) Info("AddCharge","### In: %.3g, OUT: %.3g   =>   %.3g",charge, chargeSum, chargeSum/charge);

  return 1;

}


Bool_t HarpoSimReadout::isBlocked(Double_t x, Double_t y)
{

  if(TMath::Abs(fmod(TMath::Abs(x),33)-16.5)<0.5) return kTRUE; // GEM sector

  if(TMath::Abs(fmod(TMath::Abs(x),100)-50)<1) return kTRUE; // separator
  if(TMath::Abs(fmod(TMath::Abs(y),100)-50)<1) return kTRUE; // separator

  return kFALSE;
}