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