TPC/fastSimul/AliTPCclusterFast.cxx

   1 /*
   2   gSystem->Load("libSTAT.so");
   3   .x ~/NimStyle.C
   4
   5  .L $ALICE_ROOT/TPC/fastSimul/AliTPCclusterFast.cxx+
   6  //
   7  AliTPCclusterFast::fPRF = new TF1("fprf","gausn");
   8  AliTPCclusterFast::fTRF = new TF1("ftrf","gausn");
   9  AliTPCclusterFast::fPRF->SetParameters(1,0,0.5);
  10  AliTPCclusterFast::fTRF->SetParameters(1,0,0.5);
  11  //
  12
  13  AliTPCclusterFast::Simul("aaa.root",50000);
  14  gSystem->Load("libSTAT.so");
  15
  16  TFile f("aaa.root");
  17  TTree * tree = (TTree*)f.Get("simul");
  18
  19 */
  20
  21 #include "TObject.h"
  22 #include "TF1.h"
  23 #include "TMath.h"
  24 #include "TRandom.h"
  25 #include "TVectorD.h"
  26 #include "TMatrixD.h"
  27 #include "TTreeStream.h"
  28
  29 class AliTPCclusterFast: public TObject {
  30 public:
  31   AliTPCclusterFast();
  32   virtual ~AliTPCclusterFast();
  33   void SetParam(Float_t mnprim, Float_t diff, Float_t y, Float_t z, Float_t ky, Float_t kz);
  34   void GenerElectrons();
  35   void Digitize();
  36   Double_t GetNsec();
  37   static void Simul(const char* simul, Int_t npoints);
  38 public:
  39   Float_t fMNprim;     // mean number of primary electrons
  40   //                   //electrons part input
  41   Int_t   fNprim;      // mean number of primary electrons
  42   Int_t   fNtot;       // total number of primary electrons
  43   Float_t fDiff;       // diffusion sigma
  44   Float_t fY;          // y position
  45   Float_t fZ;          // z postion
  46   Float_t fAngleY;     // y angle - tan(y)
  47   Float_t fAngleZ;     // z angle - tan z
  48   //
  49   //
  50   //                   // electron part simul
  51   TVectorD fSec;       // number of secondary electrons
  52   TVectorD fPosY;      //! position y for each electron
  53   TVectorD fPosZ;      //! position z for each electron
  54   TVectorD fGain;      //! gg for each electron
  55   //
  56   TVectorD fStatY;     // stat Y
  57   TVectorD fStatZ;     // stat Y
  58   //
  59   // digitization part
  60   //
  61   TMatrixD fDigits;    // response matrix
  62   static TF1* fPRF;    // Pad response
  63   static TF1* fTRF;    // Time response function
  64   ClassDef(AliTPCclusterFast,1)  // container for
  65 };
  66 ClassImp(AliTPCclusterFast)
  67
  68 TF1 *AliTPCclusterFast::fPRF=0;
  69 TF1 *AliTPCclusterFast::fTRF=0;
  70
  71
  72 AliTPCclusterFast::AliTPCclusterFast(){
  73   //
  74   //
  75   fDigits.ResizeTo(5,7);
  76 }
  77
  78 AliTPCclusterFast::~AliTPCclusterFast(){
  79 }
  80
  81
  82 void AliTPCclusterFast::SetParam(Float_t mnprim, Float_t diff, Float_t y, Float_t z, Float_t ky, Float_t kz){
  83   //
  84   //
  85   fMNprim = mnprim; fDiff = diff;
  86   fY=y; fZ=z;
  87   fAngleY=ky; fAngleZ=kz;
  88 }
  89 Double_t AliTPCclusterFast::GetNsec(){
  90   //
  91   // Generate number of secondary electrons
  92   // copy of procedure implemented in geant
  93   //
  94   const Double_t FPOT=20.77E-9, EEND=10E-6, EEXPO=2.2, EEND1=1E-6;
  95   const Double_t XEXPO=-EEXPO+1, YEXPO=1/XEXPO;
  96   const Double_t W=20.77E-9;
  97   Float_t RAN = gRandom->Rndm();
  98   return TMath::Nint(TMath::Power((TMath::Power(FPOT,XEXPO)*(1-RAN)+TMath::Power(EEND,XEXPO)*RAN),YEXPO)/W);
  99 }
 100
 101 void AliTPCclusterFast::GenerElectrons(){
 102   //
 103   //
 104   //
 105   //
 106   const Int_t knMax=1000;
 107   if (fPosY.GetNrows()<knMax){
 108     fPosY.ResizeTo(knMax);
 109     fPosZ.ResizeTo(knMax);
 110     fGain.ResizeTo(knMax);
 111     fSec.ResizeTo(knMax);
 112     fStatY.ResizeTo(3);
 113     fStatZ.ResizeTo(3);
 114   }
 115   fNprim = gRandom->Poisson(fMNprim);  //number of primary electrons
 116   fNtot=0;
 117   //
 118   Double_t sumQ=0;
 119   Double_t sumYQ=0;
 120   Double_t sumZQ=0;
 121   Double_t sumY2Q=0;
 122   Double_t sumZ2Q=0;
 123   for (Int_t i=0;i<knMax;i++){
 124     fSec[i]=0;
 125   }
 126   for (Int_t iprim=0; iprim<fNprim;iprim++){
 127     Float_t dN   =  GetNsec();
 128     fSec[iprim]=dN;
 129     Double_t yc = fY+(gRandom->Rndm()-0.5)*fAngleY;
 130     Double_t zc = fZ+(gRandom->Rndm()-0.5)*fAngleZ;
 131     for (Int_t isec=0;isec<=dN;isec++){
 132       //
 133       //
 134       Double_t y = gRandom->Gaus(0,fDiff)+yc;
 135       Double_t z = gRandom->Gaus(0,fDiff)+zc;
 136       Double_t gg = -TMath::Log(gRandom->Rndm());
 137       fPosY[fNtot]=y;
 138       fPosZ[fNtot]=z;
 139       fGain[fNtot]=gg;
 140       fNtot++;
 141       sumQ+=gg;
 142       sumYQ+=gg*y;
 143       sumY2Q+=gg*y*y;
 144       sumZQ+=gg*z;
 145       sumZ2Q+=gg*z*z;
 146       if (fNtot>=knMax) break;
 147     }
 148     if (fNtot>=knMax) break;
 149   }
 150   if (sumQ>0){
 151     fStatY[0]=sumQ;
 152     fStatY[1]=sumYQ/sumQ;
 153     fStatY[2]=sumY2Q/sumQ-fStatY[1]*fStatY[1];
 154     fStatZ[0]=sumQ;
 155     fStatZ[1]=sumZQ/sumQ;
 156     fStatZ[2]=sumZ2Q/sumQ-fStatZ[1]*fStatZ[1];
 157   }
 158 }
 159
 160 void AliTPCclusterFast::Digitize(){
 161   //
 162   //
 163   //
 164   // 1. Clear digits
 165   for (Int_t i=0; i<5;i++)
 166     for (Int_t j=0; j<7;j++){
 167       fDigits(i,j)=0;
 168     }
 169   //
 170   // Fill digits
 171   for (Int_t iel = 0; iel<fNtot; iel++){
 172     for (Int_t di=-2; di<=2;di++)
 173       for (Int_t dj=-2; dj<=2;dj++){
 174         Float_t fac = fPRF->Eval(di-fPosY[iel])*fTRF->Eval(dj-fPosZ[iel]);
 175         fac*=fGain[iel];
 176         fDigits(2+di,3+dj)+=fac;
 177       }
 178   }
 179 }
 180
 181
 182
 183 void AliTPCclusterFast::Simul(const char* fname, Int_t npoints){
 184   AliTPCclusterFast fast;
 185   TTreeSRedirector cstream(fname);
 186   for (Int_t icl=0; icl<npoints; icl++){
 187     Float_t nprim=(10+40*gRandom->Rndm());
 188     Float_t diff =0.01 +0.3*gRandom->Rndm();
 189     Float_t posY = gRandom->Rndm()-0.5;
 190     Float_t posZ = gRandom->Rndm()-0.5;
 191     Float_t ky   = 1.*(gRandom->Rndm()-0.5);
 192     Float_t kz   = 1.*(gRandom->Rndm()-0.5);
 193     fast.SetParam(nprim,diff,posY,posZ,ky,kz);
 194     fast.GenerElectrons();
 195     fast.Digitize();
 196     cstream<<"simul"<<
 197       "s.="<<&fast<<
 198       "\n";
 199   }
 200 }
 201
 202
 203 /*
 204   TH2F *hisL = new TH2F("hisL","hisL",10,10,50,100,0,10)
 205
 206 */