TPC/fastSimul/AliTPCclusterFast.cxx

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