TPC/fastSimul/AliTPCclusterFast.cxx

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