PndMdtParamDigi.h

Go to the documentation of this file.

//==================================================================
//description
//parameterized digitizaiton model from garfield simulation
//author, Jifeng Hu, hu@to.infn.it, Torino University
//

#ifndef PNDMDTPARAMDIGI_H
#define PNDMDTPARAMDIGI_H 1

#include "TVector3.h"
#include "TGraphErrors.h"
#include <vector>
#include "TF1.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TRandom.h"
#include <math.h>
#include <map>

//value and its error
typedef std::pair<Double_t, Double_t> ValueErrorType;

namespace CLHEP{
class RandGeneral;
}

class PndMdtParamDigi : public TNamed
{

  public:

    PndMdtParamDigi();

    ~PndMdtParamDigi();
    void SetVerbose(Int_t v) { fVerbose = v; }
    //call Init before invoking the other interfaces.
    Bool_t Init() ;
    //signal on Anode wire
    //==========================================================================
    // Int_t particleType: electron, muon, pion, kaon+-, proton, indiced as 0, 1, 2, 3, 4
    // TVector3 iniP     : initial 3-momentum of this charged particle
    // TVector3 iniPos   : initial position when entering into a tube
    // TVector3 fnlPos   : final position when exiting a tube
    //==========================================================================
    //by default, strip lenght 100 cm
    //this inteface should be invoked first when parameters change
    PndMdtParamDigi& SetParams(Int_t ptlType, TVector3 iniP, TVector3 iniPos, TVector3 finalPos, Double_t stripLen=100.);

    void UseNoise(Bool_t swith)             { fUseNoise = swith; }
    void UseDetailedSim(Bool_t swith=kTRUE) { fDetailedSim = swith; }
    void UsePlot(Bool_t swith=kTRUE)        { fUsePlot = swith; }
    void UseGaussianAmp(Bool_t swith)       { fGaussianAmp = swith; }
    void SetOptimization(Int_t val)         { fNumofTruncation = val; }
    void SetNoiseWidth(Double_t anode, Double_t strip) { fNoiseSigmaAnode = anode; fNoiseSigmaStrip = strip; }
    //===============================================================
    //readout interfaces
    //only fired information is provided 
    //for wire, <-1, fired time>, for strips, <strip number from 0, fired time>
    std::vector<std::pair<Int_t, Double_t> > GetFiredInfo();
    //both anode and strip signals
    void Compute(Bool_t useConvolution=kTRUE) ;
    //access data of induced current
    const std::vector<Double_t>& GetWireSignal() const { return fSignalDataAnode; }
    const std::map<Int_t, std::vector<Double_t> > & GetStripSignals() const { return fSignalDataStripM; }
    //===============================================================

    //plot a charged track and its produced signal 
    void Draw();

  private:
    //anode signal only
    void GetSignal(Bool_t useConvolution=kTRUE);
    //get threshold-crossing time and amplitude
    // if tube is fired, return kTRUE;
    // in case of kTRUE, time and amplitude will be set.
    Bool_t Digitize(Double_t& time, Double_t& amp);//anode
    Bool_t Digitize(Int_t stripNo, Double_t& time, Double_t& amp);//strip
    //
    void GetRawSignalbySimAvalanche(Double_t fNoiseLevel=1.);
    void GetRawSignalbyWeightingAvalanche(Double_t fNoiseLevel = 1.);
    void AddNoise(Double_t fNoiseLevel = 1., Int_t isAnode=1);
    void ApplyTransferFunction(Double_t* fSignalData, Int_t nSize=fSamplingSize);

    //by default, sampling rate is set to 100 MHz, 200 sampling points
    static const Int_t fSamplingSize = 20;
    Int_t fSamplingRate;
    Int_t fSamplingInterval;
    std::vector<Double_t> fSignalDataAnode;
    std::map<Int_t, std::vector<Double_t> > fSignalDataStripM;//strip no, and its corresponding signal
    TF1* fTrF;//transfer function
    Bool_t fUseNoise;
    Bool_t fGaussianAmp;//amplication model
    Double_t fNoiseSigmaAnode ;//0.01 mu Ampere
    Double_t fNoiseSigmaStrip ;//0.01 mu Ampere
    Int_t fNumofTruncation;


  private:
    //drift time of primary ionized electrons to wire surface
    inline Double_t GetElectronDriftTime(TVector2 iniPos)
    {
      const Double_t  p0 =  5.35513e+03;
      const Double_t  p1 =  3.92749e+02;
      const Double_t  p2 = -5.41433e+03;
      const Double_t  p3 = -3.31810e+03;
      Double_t x = iniPos.Mod() - cWIRERADIUS;
      Double_t x1 = log(1.+x);
      return (x*(p0 + p1*pow(x, 0.25)) + x1*(p2 + p3*x1));
    }
    //drift time of secondary ions from wire surface to a given point
    inline Double_t GetShiftTime(TVector2 fWireSurfacePos, TVector2 fProductionPos)
    {
      const Double_t p0 = 1.81940e-8;//unit: micro sencond
      const Double_t p1 = 1.80248;
      const Double_t p2 = 3.68652e+2;
      const Double_t p3 = 1.67353e+3;
      //const Double_t mean = 0; //[R.K. 01/2017] unused variable?
      //const Double_t sigma = 2.73477e-3;//micro second //[R.K. 01/2017] unused variable?
      Double_t dr = (fProductionPos-fWireSurfacePos).Mod();
      return (((p3*dr+p2)*dr+p1)*dr+p0 /*+ gRandom->Gaus(mean, sigma)*/)*1.e3;//nano sencod
    }

  private:
    Int_t    fParticleType;//electron 0, muon 1, pion 2, 
    TVector3 fParticleMomentum;// momentum of the incident charged particle
    TVector3 fInitPostion;
    TVector3 fExitPosition;
    Bool_t   fParamsChanged;
    Double_t fStripLength;
    Double_t fRestMass[5];


    void GetMPVofPrimaryIonization(Int_t particleType, const TVector3& momentum, ValueErrorType& val) const;
    Int_t GetAmplicationFactor(Int_t particleType, Double_t momentum) const ;
    //sampling position of secondary ion/electron for a given center 
    Bool_t fDetailedSim;
    void SamplingPosition(TVector2 fDirection, TVector2& fIonProductionPos);
    //mean free path function as velocity
    TGraphErrors* gFreePath[5];//

    inline Double_t GetMeanFreePath(Int_t ptlType, Double_t mom) const 
    {
      Double_t mMass = fRestMass[ptlType];
      //Double_t mEnergy = sqrt(mom*mom + mMass*mMass);
      Double_t mBeta = mom/mMass;
      Double_t mPath = gFreePath[ptlType]->Eval(mBeta);
      return gRandom->Exp(mPath);
    }

    //avalanche size function
    Bool_t fUseAvaHist;
    TH1F* hAvaSize;
    //TH1F* hAnodeI;
    //TH1F* hCathodeIx;
    //TH1F* hCathodeIy;
    //
    static const Int_t NRAD= 100;//radial bins, not changeable
    static const Int_t NPHI= 100;//phi bins, not changeable
    //geometry constant
    Double_t cWIRERADIUS ;//unit: cm
    Double_t cCELLSIZE   ;//unit: cm
    Double_t cSTRIPWIDTH ;//unit: cm
    Double_t cMAXRADIUS  ;//not changeable
    Double_t cMAPFACTOR  ;//not changeable
    Double_t cUNITLOGDR ;//= TMath::Log(1. + cMAPFACTOR*(cMAXRADIUS - cWIRERADIUS))/NRAD;
    Double_t cUNITPHI ;//= TMath::Pi()*2./NPHI;
    TH1F* hAnodeI1d[NRAD];
    TH1F* hCathodeI2d[NRAD][NPHI];
    Int_t fNr[NRAD];
    Int_t fNrNphi[NRAD*NPHI];

    //cluster initialization
    struct ClusInfo 
    {
      explicit ClusInfo(TVector3 _v=TVector3()): fPosition(_v) {}
      TVector3 fPosition;
      Int_t    fNumofPrimaryIonization;
      Int_t    fNumofIonsofThisCluster;
      Double_t fStartTime;
    };


    //histograms 
    Bool_t fUsePlot;
    TH1F* hIndex;
    TH2F* h2dIndex;
    TH2F* hTrack;
    TH1F* hIonNum;
    TH2F* hAvaPos;
    TH1F* hAvaNum;
    TH1F* hAva1D;
    TH2F* hAva2D;
    //
    TH1F* hExp1;
    TH1F* hExp2;
    TH1F* hGen;
    TH1F* hLan;
    Int_t fVerbose;
    //
    public:
    struct AvaBinType{
      AvaBinType(Int_t _i=0, Double_t _v=0.) : Index(_i), Probabilty(_v){}
      Int_t Index;
      Double_t Probabilty;
    };//index, density
    private:
    std::vector<AvaBinType> fProbFunc1D;
    std::vector<AvaBinType> fProbFunc2D;

    //
    Double_t SPEEDOFLIGHT ;
    CLHEP::RandGeneral* fRandAva;

    ClassDef(PndMdtParamDigi,1);
};

#endif