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",-5,5);
   8  AliTPCclusterFast::fTRF = new TF1("ftrf","gausn",-5,5);
   9  AliTPCclusterFast::fPRF->SetParameters(1,0,0.5);
  10  AliTPCclusterFast::fTRF->SetParameters(1,0,0.5);
  11
  12  //
  13  AliTPCtrackFast::Simul("trackerSimul.root",100);
  14 // AliTPCclusterFast::Simul("cluterSimul.root",20000);
  15 */
  16
  17 #include "TObject.h"
  18 #include "TF1.h"
  19 #include "TMath.h"
  20 #include "TRandom.h"
  21 #include "TVectorD.h"
  22 #include "TMatrixD.h"
  23 #include "TH1.h"
  24 #include "THnSparse.h"
  25 #include "TClonesArray.h"
  26 #include "TTreeStream.h"
  27
  28 class AliTPCclusterFast: public TObject {
  29 public:
  30   AliTPCclusterFast();
  31   virtual ~AliTPCclusterFast();
  32   void SetParam(Float_t mnprim, Float_t diff, Float_t y, Float_t z, Float_t ky, Float_t kz);
  33   void GenerElectrons();
  34   void Digitize();
  35   Double_t GetQtot(Float_t gain,Float_t thr, Float_t noise, Bool_t rounding=kTRUE, Bool_t addPedestal=kTRUE);
  36   Double_t GetQmax(Float_t gain,Float_t thr, Float_t noise, Bool_t rounding=kTRUE, Bool_t addPedestal=kTRUE);
  37   Double_t GetQmaxCorr(Float_t rmsy0, Float_t rmsz0);
  38   Double_t GetQtotCorr(Float_t rmsy0, Float_t rmsz0, Float_t gain, Float_t thr);
  39
  40   Double_t GetNsec();
  41   static void Simul(const char* simul, Int_t npoints);
  42   static Double_t GaussConvolution(Double_t x0, Double_t x1, Double_t k0, Double_t k1, Double_t s0, Double_t s1);
  43   static Double_t GaussExpConvolution(Double_t x0, Double_t s0,Double_t t1);
  44   static Double_t GaussGamma4(Double_t x, Double_t s0, Double_t p1);
  45   static Double_t Gamma4(Double_t x, Double_t p0, Double_t p1);
  46 public:
  47   Float_t fMNprim;     // mean number of primary electrons
  48   //                   //electrons part input
  49   Int_t   fNprim;      // mean number of primary electrons
  50   Int_t   fNtot;       // total number of  electrons
  51   Float_t fQtot;       // total charge - Gas gain flucuation taken into account
  52   //
  53   Float_t fDiff;       // diffusion sigma
  54   Float_t fY;          // y position
  55   Float_t fZ;          // z postion
  56   Float_t fAngleY;     // y angle - tan(y)
  57   Float_t fAngleZ;     // z angle - tan z
  58   //
  59   //
  60   //                   // electron part simul
  61   TVectorD fSec;       //! number of secondary electrons
  62   TVectorD fPosY;      //! position y for each electron
  63   TVectorD fPosZ;      //! position z for each electron
  64   TVectorD fGain;      //! gg for each electron
  65   //
  66   TVectorD fStatY;     //!stat Y
  67   TVectorD fStatZ;     //!stat Y
  68   //
  69   // digitization part
  70   //
  71   TMatrixD fDigits;    // response matrix
  72   static TF1* fPRF;    // Pad response
  73   static TF1* fTRF;    // Time response function
  74   ClassDef(AliTPCclusterFast,1)  // container for
  75 };
  76
  77
  78 class AliTPCtrackFast: public TObject {
  79 public:
  80   AliTPCtrackFast();
  81   void Add(AliTPCtrackFast &track2);
  82   void MakeTrack();
  83   void UpdatedEdxHisto();
  84   void MakeHisto();
  85   static void Simul(const char* simul, Int_t ntracks);
  86   Double_t  CookdEdxNtot(Double_t f0,Float_t f1);
  87   Double_t  CookdEdxQtot(Double_t f0,Float_t f1);
  88   //
  89   Double_t  CookdEdxDtot(Double_t f0,Float_t f1, Float_t gain,Float_t thr, Float_t noise, Bool_t corr = kTRUE);
  90   Double_t  CookdEdxDmax(Double_t f0,Float_t f1,Float_t gain,Float_t thr, Float_t noise, Bool_t corr=kTRUE);
  91   //
  92   Double_t  CookdEdx(Int_t npoints, Double_t *amp, Double_t f0,Float_t f1);
  93   //
  94   Float_t fMNprim;     // mean number of primary electrons
  95   Float_t fAngleY;     // y angle - tan(y)
  96   Float_t fAngleZ;     // z angle - tan z
  97   Float_t fDiff;       // diffusion
  98   Int_t   fN;          // number of clusters
  99   TClonesArray *fCl;   // array of clusters
 100   //
 101   Bool_t   fInit;      // initialization flag
 102   THnSparse    *fHistoNtot;    // histograms of trunc mean Ntot
 103   THnSparse    *fHistoQtot;    // histograms of trunc mean Qtot
 104   THnSparse    *fHistoQNtot;   // histograms of trunc mean Qtot/Ntot
 105   //
 106   THnSparse    *fHistoDtot;    // histograms of trunc mean digit tot
 107   THnSparse    *fHistoDmax;    // histograms of trunc mean digit max
 108   THnSparse    *fHistoDtotRaw;    // histograms of trunc mean digit tot
 109   THnSparse    *fHistoDmaxRaw;    // histograms of trunc mean digit max
 110
 111   //
 112   //
 113   ClassDef(AliTPCtrackFast,2)  // container for
 114 };
 115
 116
 117
 118 ClassImp(AliTPCclusterFast)
 119 ClassImp(AliTPCtrackFast)
 120
 121
 122
 123
 124
 125 TF1 *AliTPCclusterFast::fPRF=0;
 126 TF1 *AliTPCclusterFast::fTRF=0;
 127
 128
 129 AliTPCtrackFast::AliTPCtrackFast():
 130   TObject(),
 131   fMNprim(0),
 132   fAngleY(0),
 133   fAngleZ(0),
 134   fN(0),
 135   fCl(0),
 136   fInit(kFALSE),
 137   fHistoNtot(0),
 138   fHistoQtot(0),
 139   fHistoQNtot(0),
 140   fHistoDtot(0),
 141   fHistoDmax(0),
 142   fHistoDtotRaw(0),
 143   fHistoDmaxRaw(0)
 144 {
 145   //
 146   //
 147   //
 148 }
 149
 150 void AliTPCtrackFast::Add(AliTPCtrackFast &track2){
 151   if (!track2.fInit) return;
 152
 153   fHistoNtot->Add(track2.fHistoNtot);    // histograms of trunc mean Ntot
 154   fHistoQtot->Add(track2.fHistoQtot);    // histograms of trunc mean Qtot
 155   fHistoQNtot->Add(track2.fHistoQNtot);   // histograms of trunc mean Qtot/Ntot
 156   //
 157   fHistoDtot->Add(track2.fHistoDtot);    // histograms of trunc mean digit tot
 158   fHistoDmax->Add(track2.fHistoDmax);    // histograms of trunc mean digit max
 159   fHistoDtotRaw->Add(track2.fHistoDtotRaw);    // histograms of trunc mean digit tot
 160   fHistoDmaxRaw->Add(track2.fHistoDmaxRaw);    // histograms of trunc mean digit max
 161 }
 162
 163 void AliTPCtrackFast::MakeHisto(){
 164   //
 165   // make default histo
 166   //
 167   // dEdx histogram THnSparse
 168   // 0 - value
 169   // 1 - fMNprim - number of generated primaries
 170   // 2 - fNpoints
 171   // 3 - fFraction
 172   // 4 - fDiff
 173   // 5 - fAngleY
 174   // 6 - fAngleZ
 175
 176   Double_t xmin[7],  xmax[7];
 177   Int_t    nbins[7];
 178   if (fInit) return;
 179   //
 180   nbins[1] = 10; xmin[1]=10;  xmax[1]=30;    // fMNprim
 181   nbins[2] = 8;  xmin[2]=80;  xmax[2]=160;   // fNPoints
 182   nbins[3] = 6;  xmin[3]=0.45; xmax[3]=1.05;     // trunc mean fraction
 183
 184   nbins[4] = 5;  xmin[4]=0.0; xmax[4]=0.4;   // fDiff
 185   nbins[5] = 10; xmin[5]=0;   xmax[5]=2;     // fAngleY
 186   nbins[6] = 10; xmin[6]=0;   xmax[6]=2;     // fAngleZ
 187   //
 188   nbins[0] =100; xmin[0]=2; xmax[0]=8;
 189   fHistoNtot = new THnSparseF("dNdxall/dNdxprim","dNdxall/dNdxprim", 4, nbins, xmin,xmax);
 190   nbins[0] =100; xmin[0]=2; xmax[0]=8;
 191   fHistoQtot = new THnSparseF("dQdx/dNdxprim","dQdxall/dNdxprim", 4, nbins, xmin,xmax);
 192   nbins[0] =100; xmin[0]=0.5; xmax[0]=1.5;
 193   fHistoQNtot = new THnSparseF("dQdx/dNdxprim","dQdxprim/dNdxprim", 4, nbins, xmin,xmax);
 194   //
 195   nbins[0] =100; xmin[0]=0.05; xmax[0]=8;
 196   fHistoDtot = new THnSparseF("dQtotdx/dNdxprim","dQtotdx/dNdx", 7, nbins, xmin,xmax);
 197   fHistoDmax = new THnSparseF("dQmaxdx/dNdxprim","dQmaxdx/dNdx", 7, nbins, xmin,xmax);
 198   fHistoDtotRaw = new THnSparseF("raw dQtotdx/dNdxprim","raw dQtotdx/dNdx", 7, nbins, xmin,xmax);
 199   fHistoDmaxRaw = new THnSparseF("raw dQmaxdx/dNdxprim","raw dQmaxdx/dNdx", 7, nbins, xmin,xmax);
 200   fInit=kTRUE;
 201 }
 202
 203 void  AliTPCtrackFast::UpdatedEdxHisto(){
 204   //
 205   //fill default histo
 206   //
 207   if (!fInit) MakeHisto();
 208   Double_t x[7];
 209   x[1] = fMNprim;
 210   x[2] = fN;
 211   //
 212   x[4] = fDiff;
 213   x[5] = TMath::Abs(fAngleY);
 214   x[6] = TMath::Abs(fAngleZ);
 215
 216   for (Int_t i=0;i<7;i++){
 217     Float_t frac = 0.5+Float_t(i)*0.1;
 218     x[3] = frac;
 219     Double_t cNtot = CookdEdxNtot(0.01,frac);
 220     Double_t cQtot = CookdEdxQtot(0.01,frac);
 221     // MC -using hits
 222     x[0] = cNtot/fMNprim;
 223     fHistoNtot->Fill(x);
 224     x[0] = cQtot/fMNprim;
 225     fHistoQtot->Fill(x);
 226     x[0] = cQtot/cNtot;
 227     fHistoQNtot->Fill(x);
 228     // MC - using digits
 229     Double_t dQtot = CookdEdxDtot(0.01,frac,1,2.5,1,kTRUE);
 230     Double_t dQmax = CookdEdxDmax(0.01,frac,1,2.5,1,kTRUE);
 231     Double_t dQrawtot = CookdEdxDtot(0.01,frac,1,2.5,1,kFALSE);
 232     Double_t dQrawmax = CookdEdxDmax(0.01,frac,1,2.5,1,kFALSE);
 233     x[0] = dQtot/fMNprim;
 234     fHistoDtot->Fill(x);
 235     x[0] = dQmax/fMNprim;
 236     fHistoDmax->Fill(x);
 237     x[0] = dQrawtot/fMNprim;
 238     fHistoDtotRaw->Fill(x);
 239     x[0] = dQrawmax/fMNprim;
 240     fHistoDmaxRaw->Fill(x);
 241   }
 242 }
 243
 244 void AliTPCtrackFast::MakeTrack(){
 245   //
 246   //
 247   //
 248   if (!fCl) fCl = new TClonesArray("AliTPCclusterFast",160);
 249   for (Int_t i=0;i<fN;i++){
 250     Double_t tY = i*fAngleY;
 251     Double_t tZ = i*fAngleZ;
 252     AliTPCclusterFast * cluster = (AliTPCclusterFast*) fCl->UncheckedAt(i);
 253     if (!cluster) cluster =   new ((*fCl)[i]) AliTPCclusterFast;
 254     //
 255     Double_t posY = tY-TMath::Nint(tY);
 256     Double_t posZ = tZ-TMath::Nint(tZ);
 257     cluster->SetParam(fMNprim,fDiff,posY,posZ,fAngleY,fAngleZ);
 258     //
 259     cluster->GenerElectrons();
 260     cluster->Digitize();
 261   }
 262   UpdatedEdxHisto();
 263 }
 264
 265 Double_t  AliTPCtrackFast::CookdEdxNtot(Double_t f0,Float_t f1){
 266   //
 267   Double_t amp[160];
 268   for (Int_t i=0;i<fN;i++){
 269     AliTPCclusterFast * cluster = ( AliTPCclusterFast *)((*fCl)[i]);
 270     amp[i]=cluster->fNtot;
 271   }
 272   return CookdEdx(fN,amp,f0,f1);
 273 }
 274
 275 Double_t  AliTPCtrackFast::CookdEdxQtot(Double_t f0,Float_t f1){
 276   //
 277   Double_t amp[160];
 278   for (Int_t i=0;i<fN;i++){
 279     AliTPCclusterFast * cluster = ( AliTPCclusterFast *)((*fCl)[i]);
 280     amp[i]=cluster->fQtot;
 281   }
 282   return CookdEdx(fN,amp,f0,f1);
 283 }
 284
 285 Double_t   AliTPCtrackFast::CookdEdxDtot(Double_t f0,Float_t f1, Float_t gain,Float_t thr, Float_t noise, Bool_t doCorr){
 286   //
 287   //
 288   //
 289   Double_t amp[160];
 290   Int_t over=0;
 291   for (Int_t i=0;i<fN;i++){
 292     AliTPCclusterFast * cluster = ( AliTPCclusterFast *)((*fCl)[i]);
 293     Float_t camp = cluster->GetQtot(gain,thr,noise);
 294     if (camp==0) continue;
 295     Float_t corr =  1;
 296     if (doCorr) corr = cluster->GetQtotCorr(0.5,0.5,gain,thr);
 297     amp[over]=camp/corr;
 298     over++;
 299   }
 300   return CookdEdx(over,amp,f0,f1);
 301
 302 }
 303
 304 Double_t   AliTPCtrackFast::CookdEdxDmax(Double_t f0,Float_t f1, Float_t gain,Float_t thr, Float_t noise, Bool_t doCorr){
 305   //
 306   //
 307   //
 308   Double_t amp[160];
 309   Int_t over=0;
 310   for (Int_t i=0;i<fN;i++){
 311     AliTPCclusterFast * cluster = ( AliTPCclusterFast *)((*fCl)[i]);
 312     Float_t camp = cluster->GetQmax(gain,thr,noise);
 313     if (camp==0) continue;
 314     Float_t corr =  1;
 315     if (doCorr) corr = cluster->GetQmaxCorr(0.5,0.5);
 316     amp[over]=camp/corr;
 317     over++;
 318   }
 319   return CookdEdx(over,amp,f0,f1);
 320
 321 }
 322
 323
 324 Double_t  AliTPCtrackFast::CookdEdx(Int_t npoints, Double_t *amp,Double_t f0,Float_t f1){
 325   //
 326   //
 327   //
 328   Int_t index[160];
 329   TMath::Sort(npoints,amp,index,kFALSE);
 330   Float_t sum0=0, sum1=0,sum2=0;
 331   for (Int_t i=0;i<npoints;i++){
 332     if (i<npoints*f0) continue;
 333     if (i>npoints*f1) continue;
 334     sum0++;
 335     sum1+= amp[index[i]];
 336     sum2+= amp[index[i]];
 337   }
 338   if (sum0<=0) return 0;
 339   return sum1/sum0;
 340 }
 341
 342 void AliTPCtrackFast::Simul(const char* fname, Int_t ntracks){
 343   //
 344   //
 345   //
 346   AliTPCtrackFast fast;
 347   TTreeSRedirector cstream(fname);
 348   for (Int_t itr=0; itr<ntracks; itr++){
 349     //
 350     fast.fMNprim=(5+50*gRandom->Rndm());
 351     fast.fDiff =0.01 +0.35*gRandom->Rndm();
 352     //
 353     fast.fAngleY   = 4.0*(gRandom->Rndm()-0.5);
 354     fast.fAngleZ   = 4.0*(gRandom->Rndm()-0.5);
 355     fast.fN  = TMath::Nint(80.+gRandom->Rndm()*80.);
 356     fast.MakeTrack();
 357     if (itr%100==0) printf("%d\n",itr);
 358     cstream<<"simulTrack"<<
 359       "tr.="<<&fast<<
 360       "\n";
 361   }
 362   fast.Write("track");
 363 }
 364
 365
 366
 367 AliTPCclusterFast::AliTPCclusterFast(){
 368   //
 369   //
 370   fDigits.ResizeTo(5,7);
 371 }
 372
 373 AliTPCclusterFast::~AliTPCclusterFast(){
 374 }
 375
 376
 377 void AliTPCclusterFast::SetParam(Float_t mnprim, Float_t diff, Float_t y, Float_t z, Float_t ky, Float_t kz){
 378   //
 379   //
 380   fMNprim = mnprim; fDiff = diff;
 381   fY=y; fZ=z;
 382   fAngleY=ky; fAngleZ=kz;
 383 }
 384 Double_t AliTPCclusterFast::GetNsec(){
 385   //
 386   // Generate number of secondary electrons
 387   // copy of procedure implemented in geant
 388   //
 389   const Double_t FPOT=20.77E-9, EEND=10E-6, EEXPO=2.2, EEND1=1E-6;
 390   const Double_t XEXPO=-EEXPO+1, YEXPO=1/XEXPO;
 391   const Double_t W=20.77E-9;
 392   Float_t RAN = gRandom->Rndm();
 393   //Double_t edep = TMath::Power((TMath::Power(FPOT,XEXPO)*(1-RAN)+TMath::Power(EEND,XEXPO)*RAN),YEXPO);
 394   //edep = TMath::Min(edep, EEND);
 395   //return TMath::Nint(edep/W);
 396   return TMath::Nint(TMath::Power((TMath::Power(FPOT,XEXPO)*(1-RAN)+TMath::Power(EEND,XEXPO)*RAN),YEXPO)/W);
 397 }
 398
 399 void AliTPCclusterFast::GenerElectrons(){
 400   //
 401   //
 402   //
 403   //
 404   const Int_t knMax=1000;
 405   if (fPosY.GetNrows()<knMax){
 406     fPosY.ResizeTo(knMax);
 407     fPosZ.ResizeTo(knMax);
 408     fGain.ResizeTo(knMax);
 409     fSec.ResizeTo(knMax);
 410     fStatY.ResizeTo(3);
 411     fStatZ.ResizeTo(3);
 412   }
 413   fNprim = gRandom->Poisson(fMNprim);  //number of primary electrons
 414   fNtot=0; //total number of electrons
 415   fQtot=0; //total number of electrons after gain multiplification
 416   //
 417   Double_t sumQ=0;
 418   Double_t sumYQ=0;
 419   Double_t sumZQ=0;
 420   Double_t sumY2Q=0;
 421   Double_t sumZ2Q=0;
 422   for (Int_t i=0;i<knMax;i++){
 423     fSec[i]=0;
 424   }
 425   for (Int_t iprim=0; iprim<fNprim;iprim++){
 426     Float_t dN   =  GetNsec();
 427     fSec[iprim]=dN;
 428     Double_t yc = fY+(gRandom->Rndm()-0.5)*fAngleY;
 429     Double_t zc = fZ+(gRandom->Rndm()-0.5)*fAngleZ;
 430     for (Int_t isec=0;isec<=dN;isec++){
 431       //
 432       //
 433       Double_t y = gRandom->Gaus(0,fDiff)+yc;
 434       Double_t z = gRandom->Gaus(0,fDiff)+zc;
 435       Double_t gg = -TMath::Log(gRandom->Rndm());
 436       fPosY[fNtot]=y;
 437       fPosZ[fNtot]=z;
 438       fGain[fNtot]=gg;
 439       fQtot+=gg;
 440       fNtot++;
 441       sumQ+=gg;
 442       sumYQ+=gg*y;
 443       sumY2Q+=gg*y*y;
 444       sumZQ+=gg*z;
 445       sumZ2Q+=gg*z*z;
 446       if (fNtot>=knMax) break;
 447     }
 448     if (fNtot>=knMax) break;
 449   }
 450   if (sumQ>0){
 451     fStatY[0]=sumQ;
 452     fStatY[1]=sumYQ/sumQ;
 453     fStatY[2]=sumY2Q/sumQ-fStatY[1]*fStatY[1];
 454     fStatZ[0]=sumQ;
 455     fStatZ[1]=sumZQ/sumQ;
 456     fStatZ[2]=sumZ2Q/sumQ-fStatZ[1]*fStatZ[1];
 457   }
 458 }
 459
 460 void AliTPCclusterFast::Digitize(){
 461   //
 462   //
 463   //
 464   // 1. Clear digits
 465   for (Int_t i=0; i<5;i++)
 466     for (Int_t j=0; j<7;j++){
 467       fDigits(i,j)=0;
 468     }
 469   //
 470   // Fill digits
 471   for (Int_t iel = 0; iel<fNtot; iel++){
 472     for (Int_t di=-2; di<=2;di++)
 473       for (Int_t dj=-3; dj<=3;dj++){
 474         Float_t fac = fPRF->Eval(di-fPosY[iel])*fTRF->Eval(dj-fPosZ[iel]);
 475         fac*=fGain[iel];
 476         fDigits(2+di,3+dj)+=fac;
 477       }
 478   }
 479
 480 }
 481
 482
 483
 484 void AliTPCclusterFast::Simul(const char* fname, Int_t npoints){
 485   //
 486   // Calc rms
 487   //
 488   AliTPCclusterFast fast;
 489   TTreeSRedirector cstream(fname);
 490   for (Int_t icl=0; icl<npoints; icl++){
 491     Float_t nprim=(10+20*gRandom->Rndm());
 492     Float_t diff =0.01 +0.35*gRandom->Rndm();
 493     Float_t posY = gRandom->Rndm()-0.5;
 494     Float_t posZ = gRandom->Rndm()-0.5;
 495     //
 496     Float_t ky   = 4.0*(gRandom->Rndm()-0.5);
 497     Float_t kz   = 4.0*(gRandom->Rndm()-0.5);
 498     fast.SetParam(nprim,diff,posY,posZ,ky,kz);
 499     fast.GenerElectrons();
 500     fast.Digitize();
 501     if (icl%10000==0) printf("%d\n",icl);
 502     cstream<<"simul"<<
 503       "s.="<<&fast<<
 504       "\n";
 505   }
 506 }
 507
 508
 509 Double_t AliTPCclusterFast::GetQtot(Float_t gain, Float_t thr, Float_t noise, Bool_t brounding, Bool_t baddPedestal){
 510   //
 511   //
 512   //
 513   Float_t sum =0;
 514   for (Int_t ip=0;ip<5;ip++){
 515     Float_t pedestal=gRandom->Rndm()-0.5; //pedestal offset different for each pad
 516     for (Int_t it=0;it<7;it++){
 517       Float_t amp = gain*fDigits(ip,it)+gRandom->Gaus()*noise;
 518       if (baddPedestal) amp+=pedestal;
 519       if (brounding) amp=TMath::Nint(amp);
 520       if (amp>thr) sum+=amp;
 521     }
 522   }
 523   return sum;
 524 }
 525
 526 Double_t AliTPCclusterFast::GetQmax(Float_t gain, Float_t thr, Float_t noise, Bool_t brounding, Bool_t baddPedestal){
 527   //
 528   //
 529   //
 530   Float_t max =0;
 531   for (Int_t ip=0;ip<5;ip++){
 532     Float_t pedestal=gRandom->Rndm()-0.5; //pedestal offset different for each pad
 533     for (Int_t it=0;it<7;it++){
 534       Float_t amp = gain*fDigits(ip,it)+gRandom->Gaus()*noise;
 535       if (baddPedestal) amp+=pedestal;
 536       if (brounding) amp=TMath::Nint(amp);
 537       if (amp>max && amp>thr) max=amp;
 538     }
 539   }
 540   return max;
 541 }
 542
 543
 544
 545 Double_t  AliTPCclusterFast::GetQmaxCorr(Float_t rmsy0, Float_t rmsz0){
 546   //
 547   // Gaus distribution convolueted with rectangular
 548   // Gaus width sy and sz is determined by RF width and diffusion
 549   // Integral of Q is equal 1
 550   // Q max is calculated at position fY,fX
 551   //
 552   //
 553   //
 554   Double_t sy = TMath::Sqrt(rmsy0*rmsy0+fDiff*fDiff);
 555   Double_t sz = TMath::Sqrt(rmsz0*rmsz0+fDiff*fDiff);
 556   return GaussConvolution(fY,fZ, fAngleY,fAngleZ,sy,sz);
 557 }
 558
 559
 560 Double_t  AliTPCclusterFast::GetQtotCorr(Float_t rmsy0, Float_t rmsz0, Float_t gain, Float_t thr){
 561   //
 562   //  Calculates the fraction of the charge over threshol to total charge
 563   //  The response function
 564   //
 565   Double_t sy = TMath::Sqrt(rmsy0*rmsy0+fDiff*fDiff);
 566   Double_t sz = TMath::Sqrt(rmsz0*rmsz0+fDiff*fDiff);
 567   Double_t sumAll=0,sumThr=0;
 568   Double_t qtot = GetQtot(gain,thr,0); // sum of signal over threshold
 569   //
 570   Double_t corr =1;
 571   Double_t qnorm=qtot;
 572   for (Int_t iter=0;iter<2;iter++){
 573     for (Int_t iy=-2;iy<=2;iy++)
 574       for (Int_t iz=-2;iz<=2;iz++){
 575         Double_t val = GaussConvolution(fY-iy,fZ-iz, fAngleY,fAngleZ,sy,sz);
 576         Double_t qlocal =TMath::Nint(qnorm*val);
 577         if (qlocal>thr) sumThr+=qlocal;
 578         sumAll+=qlocal;
 579       }
 580     if (sumAll>0&&sumThr>0) corr=(sumThr)/sumAll;
 581     //corr = sumThr;
 582     if (corr>0) qnorm=qtot/corr;
 583
 584   }
 585   return corr;
 586 }
 587
 588
 589
 590
 591
 592 Double_t  AliTPCclusterFast::GaussConvolution(Double_t x0, Double_t x1, Double_t k0, Double_t k1, Double_t s0, Double_t s1){
 593   //
 594   // 2 D gaus convoluted with angular effect
 595   // See in mathematica:
 596   //Simplify[Integrate[Exp[-(x0-k0*xd)*(x0-k0*xd)/(2*s0*s0)-(x1-k1*xd)*(x1-k1*xd)/(2*s1*s1)]/(s0*s1),{xd,-1/2,1/2}]]
 597   //
 598   //TF1 f1("f1","AliTPCclusterFast::GaussConvolution(x,0,1,0,0.1,0.1)",-2,2)
 599   //TF2 f2("f2","AliTPCclusterFast::GaussConvolution(x,y,1,1,0.1,0.1)",-2,2,-2,2)
 600   //
 601   const Float_t kEpsilon = 0.0001;
 602   if ((TMath::Abs(k0)+TMath::Abs(k1))<kEpsilon*(s0+s1)){
 603     // small angular effect
 604     Double_t val = (TMath::Gaus(x0,0,s0)*TMath::Gaus(x1,0,s1))/(s0*s1*2.*TMath::Pi());
 605     return val;
 606   }
 607
 608   Double_t sigma2 = k1*k1*s0*s0+k0*k0*s1*s1;
 609   Double_t exp0 = TMath::Exp(-(k1*x0-k0*x1)*(k1*x0-k0*x1)/(2*sigma2));
 610   //
 611   Double_t sigmaErf =  2*s0*s1*TMath::Sqrt(2*sigma2);
 612   Double_t erf0 = TMath::Erf( (k0*s1*s1*(k0-2*x0)+k1*s0*s0*(k1-2*x1))/sigmaErf);
 613   Double_t erf1 = TMath::Erf( (k0*s1*s1*(k0+2*x0)+k1*s0*s0*(k1+2*x1))/sigmaErf);
 614   Double_t norm = 1./TMath::Sqrt(sigma2);
 615   norm/=2.*TMath::Sqrt(2.*TMath::Pi());
 616   Double_t val  = norm*exp0*(erf0+erf1);
 617   return val;
 618
 619 }
 620
 621
 622 Double_t  AliTPCclusterFast::GaussExpConvolution(Double_t x0, Double_t s0,Double_t t1){
 623  //
 624   // 2 D gaus convoluted with exponential
 625   // Integral nomalized to 1
 626   // See in mathematica:
 627   //Simplify[Integrate[Exp[-(x0-x1)*(x0-x1)/(2*s0*s0)]*Exp[-x1*t1],{x1,0,Infinity}]]
 628   // TF1 fgexp("fgexp","AliTPCclusterFast::GaussExpConvolution(x,0.5,1)",-2,2)
 629   Double_t exp1 = (s0*s0*t1-2*x0)*t1/2.;
 630   exp1 = TMath::Exp(exp1);
 631   Double_t erf = 1+TMath::Erf((-s0*s0*t1+x0)/(s0*TMath::Sqrt(2.)));
 632   Double_t val = exp1*erf;
 633   val *=t1/(2.);
 634   return val;
 635
 636 }
 637
 638
 639 Double_t  AliTPCclusterFast::Gamma4(Double_t x, Double_t p0, Double_t p1){
 640   //
 641   // Gamma 4 Time response function of ALTRO
 642   //
 643   if (x<0) return 0;
 644   Double_t g1 = TMath::Exp(-4.*x/p1);
 645   Double_t g2 = TMath::Power(x/p1,4);
 646   return p0*g1*g2;
 647 }
 648
 649
 650
 651 Double_t  AliTPCclusterFast::GaussGamma4(Double_t x, Double_t s0, Double_t p1){
 652   //
 653   // Gamma 4 Time response function of ALTRO convoluted with Gauss
 654   // Simplify[Integrate[Exp[-(x0-x1)*(x0-x1)/(2*s0*s0)]*Exp[-4*x1/p1]*(x/p1)^4/s0,{x1,0,Infinity}]]
 655   //TF1 fgg4("fgg4","AliTPCclusterFast::GaussGamma4(x,0.5,0.5)",-2,2)
 656
 657   Double_t exp1 = (8*s0*s0-4.*p1*x)/(p1*p1);
 658   exp1 = TMath::Exp(exp1);
 659   Double_t erf1 = 1+TMath::Erf((-4*s0/p1+x/s0)/TMath::Sqrt(2));
 660   //  Double_t xp14 = TMath::Power(TMath::Abs((x/p1)),4);
 661   return exp1*erf1;
 662
 663
 664 }
 665
 666 // Analytical sollution only in 1D - too long expression
 667 // Simplify[Integrate[Exp[-(x0-(x1-k*x2))*(x0-(x1-k*x2))/(2*s0*s0)]*Exp[-(x1*t1-k*x2)],{x2,-1,1}]]
 668 //
 669 //
 670 // No analytical solution
 671 //
 672 //Simplify[Integrate[Exp[-(x0-k0*xd)*(x0-k0*xd)/(2*s0*s0)-(x1-xt-k1*xd)*(x1-xt-k1*xd)/(2*s1*s1)]*Exp[-kt*xt]/(s0*s1),{xd,-1/2,1/2},{xt,0,Infinity}]]