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, Double_t diff);
 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, Double_t diffFactor){
 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     fast.fDiffLong =  fast.fDiff*diffFactor/1.;
 431     //
 432     fast.fAngleY   = 4.0*(gRandom->Rndm()-0.5);
 433     fast.fAngleZ   = 4.0*(gRandom->Rndm()-0.5);
 434     fast.fN  = 160;
 435     fast.MakeTrack();
 436     if (itr%100==0) printf("%d\n",itr);
 437     (*pcstream)<<"simulTrack"<<
 438       "tr.="<<&fast<<
 439       "\n";
 440   }
 441   fast.Write("track");
 442   delete pcstream;
 443 }
 444
 445
 446
 447 AliTPCclusterFast::AliTPCclusterFast(){
 448   //
 449   //
 450   fDigits.ResizeTo(5,7);
 451 }
 452
 453 void AliTPCclusterFast::Init(){
 454   //
 455   // reset all counters
 456   //
 457   const Int_t knMax=10000;
 458   fMNprim=0;     // mean number of primary electrons
 459   //                   //electrons part input
 460   fNprim=0;      // mean number of primary electrons
 461   fNtot=0;       // total number of  electrons
 462   fQtot=0;       // total charge - Gas gain flucuation taken into account
 463   //
 464   fPosY.ResizeTo(knMax);
 465   fPosZ.ResizeTo(knMax);
 466   fGain.ResizeTo(knMax);
 467   fSec.ResizeTo(knMax);
 468   fStatY.ResizeTo(3);
 469   fStatZ.ResizeTo(3);
 470   for (Int_t i=0; i<knMax; i++){
 471     fPosY[i]=0;
 472     fPosZ[i]=0;
 473     fGain[i]=0;
 474   }
 475 }
 476
 477
 478
 479 AliTPCclusterFast::~AliTPCclusterFast(){
 480 }
 481
 482
 483 void AliTPCclusterFast::SetParam(Float_t mnprim, Float_t diff,  Float_t diffL,Float_t y, Float_t z, Float_t ky, Float_t kz){
 484   //
 485   //
 486   fMNprim = mnprim; fDiff = diff; fDiffLong=diffL;
 487   fY=y; fZ=z;
 488   fAngleY=ky; fAngleZ=kz;
 489 }
 490 Double_t AliTPCclusterFast::GetNsec(){
 491   //
 492   // Generate number of secondary electrons
 493   // copy of procedure implemented in geant
 494   //
 495   const Double_t FPOT=20.77E-9, EEND=10E-6, EEXPO=2.2, EEND1=1E-6;
 496   const Double_t XEXPO=-EEXPO+1, YEXPO=1/XEXPO;
 497   const Double_t W=20.77E-9;
 498   Float_t RAN = gRandom->Rndm();
 499   //Double_t edep = TMath::Power((TMath::Power(FPOT,XEXPO)*(1-RAN)+TMath::Power(EEND,XEXPO)*RAN),YEXPO);
 500   //edep = TMath::Min(edep, EEND);
 501   //return TMath::Nint(edep/W);
 502   return TMath::Nint(TMath::Power((TMath::Power(FPOT,XEXPO)*(1-RAN)+TMath::Power(EEND,XEXPO)*RAN),YEXPO)/W);
 503 }
 504
 505 void AliTPCclusterFast::GenerElectrons(AliTPCclusterFast *cl0, AliTPCclusterFast *clm, AliTPCclusterFast *clp){
 506   //
 507   //
 508   //
 509   //
 510   const Int_t knMax=1000;
 511   cl0->fNprim = gRandom->Poisson(cl0->fMNprim);  //number of primary electrons
 512   cl0->fNtot=0; //total number of electrons
 513   cl0->fQtot=0; //total number of electrons after gain multiplification
 514   //
 515   Double_t sumQ=0;
 516   Double_t sumYQ=0;
 517   Double_t sumZQ=0;
 518   Double_t sumY2Q=0;
 519   Double_t sumZ2Q=0;
 520   for (Int_t i=0;i<knMax;i++){
 521     cl0->fSec[i]=0;
 522   }
 523   for (Int_t iprim=0; iprim<cl0->fNprim;iprim++){
 524     Float_t dN   =  cl0->GetNsec();
 525     cl0->fSec[iprim]=dN;
 526     Double_t yc = cl0->fY+(gRandom->Rndm()-0.5)*cl0->fAngleY;
 527     Double_t zc = cl0->fZ+(gRandom->Rndm()-0.5)*cl0->fAngleZ;
 528     Double_t rc = (gRandom->Rndm()-0.5);
 529
 530     for (Int_t isec=0;isec<=dN;isec++){
 531       //
 532       //
 533       Double_t y = gRandom->Gaus(0,cl0->fDiff)+yc;
 534       Double_t z = gRandom->Gaus(0,cl0->fDiff)+zc;
 535       Double_t r = gRandom->Gaus(0,cl0->fDiffLong)+rc;
 536       // choose pad row
 537       AliTPCclusterFast *cl=cl0;
 538       if (r<-0.5 &&cl) cl=clm;
 539       if (r>0.5 &&cl)  cl=clp;
 540       //
 541       Double_t gg = -TMath::Log(gRandom->Rndm());
 542       cl->fPosY[cl->fNtot]=y;
 543       cl->fPosZ[cl->fNtot]=z;
 544       cl->fGain[cl->fNtot]=gg;
 545       cl->fQtot+=gg;
 546       cl->fNtot++;
 547       //
 548   //     cl->sumQ+=gg;
 549 //       cl->sumYQ+=gg*y;
 550 //       cl->sumY2Q+=gg*y*y;
 551 //       cl->sumZQ+=gg*z;
 552 //       cl->sumZ2Q+=gg*z*z;
 553       if (cl->fNtot>=knMax) continue;
 554     }
 555     if (cl0->fNtot>=knMax) break;
 556   }
 557
 558  //  if (sumQ>0){
 559 //     fStatY[0]=sumQ;
 560 //     fStatY[1]=sumYQ/sumQ;
 561 //     fStatY[2]=sumY2Q/sumQ-fStatY[1]*fStatY[1];
 562 //     fStatZ[0]=sumQ;
 563 //     fStatZ[1]=sumZQ/sumQ;
 564 //     fStatZ[2]=sumZ2Q/sumQ-fStatZ[1]*fStatZ[1];
 565 //   }
 566 }
 567
 568 void AliTPCclusterFast::Digitize(){
 569   //
 570   //
 571   //
 572   // 1. Clear digits
 573   for (Int_t i=0; i<5;i++)
 574     for (Int_t j=0; j<7;j++){
 575       fDigits(i,j)=0;
 576     }
 577   //
 578   // Fill digits
 579   for (Int_t iel = 0; iel<fNtot; iel++){
 580     for (Int_t di=-2; di<=2;di++)
 581       for (Int_t dj=-3; dj<=3;dj++){
 582         Float_t fac = fPRF->Eval(di-fPosY[iel])*fTRF->Eval(dj-fPosZ[iel]);
 583         fac*=fGain[iel];
 584         fDigits(2+di,3+dj)+=fac;
 585       }
 586   }
 587   //
 588   //
 589   //
 590 }
 591
 592
 593
 594 // void AliTPCclusterFast::Simul(const char* fname, Int_t npoints){
 595 //   //
 596 //   // Calc rms
 597 //   //
 598 //   AliTPCclusterFast fast;
 599 //   TTreeSRedirector cstream(fname);
 600 //   for (Int_t icl=0; icl<npoints; icl++){
 601 //     Float_t nprim=(10+20*gRandom->Rndm());
 602 //     Float_t diff =0.01 +0.35*gRandom->Rndm();
 603 //     Float_t posY = gRandom->Rndm()-0.5;
 604 //     Float_t posZ = gRandom->Rndm()-0.5;
 605 //     //
 606 //     Float_t ky   = 4.0*(gRandom->Rndm()-0.5);
 607 //     Float_t kz   = 4.0*(gRandom->Rndm()-0.5);
 608 //     fast.SetParam(nprim,diff,posY,posZ,ky,kz);
 609 //     fast.GenerElectrons();
 610 //     fast.Digitize();
 611 //     if (icl%10000==0) printf("%d\n",icl);
 612 //     cstream<<"simul"<<
 613 //       "s.="<<&fast<<
 614 //       "\n";
 615 //   }
 616 // }
 617
 618
 619 Double_t AliTPCclusterFast::GetQtot(Float_t gain, Float_t thr, Float_t noise, Bool_t brounding, Bool_t baddPedestal){
 620   //
 621   //
 622   //
 623   Float_t sum =0;
 624   for (Int_t ip=0;ip<5;ip++){
 625     Float_t pedestal=gRandom->Rndm()-0.5; //pedestal offset different for each pad
 626     for (Int_t it=0;it<7;it++){
 627       Float_t amp = gain*fDigits(ip,it)+gRandom->Gaus()*noise;
 628       if (baddPedestal) amp+=pedestal;
 629       if (brounding) amp=TMath::Nint(amp);
 630       if (amp>thr) sum+=amp;
 631     }
 632   }
 633   return sum;
 634 }
 635
 636 Double_t AliTPCclusterFast::GetQmax(Float_t gain, Float_t thr, Float_t noise, Bool_t brounding, Bool_t baddPedestal){
 637   //
 638   //
 639   //
 640   Float_t max =0;
 641   for (Int_t ip=0;ip<5;ip++){
 642     Float_t pedestal=gRandom->Rndm()-0.5; //pedestal offset different for each pad
 643     for (Int_t it=0;it<7;it++){
 644       Float_t amp = gain*fDigits(ip,it)+gRandom->Gaus()*noise;
 645       if (baddPedestal) amp+=pedestal;
 646       if (brounding) amp=TMath::Nint(amp);
 647       if (amp>max && amp>thr) max=amp;
 648     }
 649   }
 650   return max;
 651 }
 652
 653
 654
 655 Double_t  AliTPCclusterFast::GetQmaxCorr(Float_t rmsy0, Float_t rmsz0){
 656   //
 657   // Gaus distribution convolueted with rectangular
 658   // Gaus width sy and sz is determined by RF width and diffusion
 659   // Integral of Q is equal 1
 660   // Q max is calculated at position fY,fX
 661   //
 662   //
 663   //
 664   Double_t sy = TMath::Sqrt(rmsy0*rmsy0+fDiff*fDiff);
 665   Double_t sz = TMath::Sqrt(rmsz0*rmsz0+fDiff*fDiff);
 666   return GaussConvolution(fY,fZ, fAngleY,fAngleZ,sy,sz);
 667 }
 668
 669
 670 Double_t  AliTPCclusterFast::GetQtotCorr(Float_t rmsy0, Float_t rmsz0, Float_t gain, Float_t thr){
 671   //
 672   //  Calculates the fraction of the charge over threshol to total charge
 673   //  The response function
 674   //
 675   Double_t sy = TMath::Sqrt(rmsy0*rmsy0+fDiff*fDiff);
 676   Double_t sz = TMath::Sqrt(rmsz0*rmsz0+fDiff*fDiff);
 677   Double_t sumAll=0,sumThr=0;
 678   Double_t qtot = GetQtot(gain,thr,0); // sum of signal over threshold
 679   Double_t qmax = GetQmax(gain,thr,0); // qmax
 680   //
 681   Double_t corr =1;
 682   Double_t qnorm=qtot;
 683   for (Int_t iter=0;iter<2;iter++){
 684     for (Int_t iy=-2;iy<=2;iy++)
 685       for (Int_t iz=-2;iz<=2;iz++){
 686         Double_t val = GaussConvolution(fY-iy,fZ-iz, fAngleY,fAngleZ,sy,sz);
 687         Double_t qlocal =TMath::Nint(qnorm*val);
 688         if (qlocal>thr) sumThr+=qlocal;
 689         sumAll+=qlocal;
 690       }
 691     if (sumAll>0&&sumThr>0) corr=(sumThr)/sumAll;
 692     if (sumThr==0) corr=GetQmaxCorr(0.5,0.5);
 693     //corr = sumThr;
 694     if (corr>0) qnorm=qtot/corr;
 695
 696   }
 697   return corr;
 698 }
 699
 700
 701
 702
 703
 704 Double_t  AliTPCclusterFast::GaussConvolution(Double_t x0, Double_t x1, Double_t k0, Double_t k1, Double_t s0, Double_t s1){
 705   //
 706   // 2 D gaus convoluted with angular effect
 707   // See in mathematica:
 708   //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}]]
 709   //
 710   //TF1 f1("f1","AliTPCclusterFast::GaussConvolution(x,0,1,0,0.1,0.1)",-2,2)
 711   //TF2 f2("f2","AliTPCclusterFast::GaussConvolution(x,y,1,1,0.1,0.1)",-2,2,-2,2)
 712   //
 713   const Float_t kEpsilon = 0.0001;
 714   if ((TMath::Abs(k0)+TMath::Abs(k1))<kEpsilon*(s0+s1)){
 715     // small angular effect
 716     Double_t val = (TMath::Gaus(x0,0,s0)*TMath::Gaus(x1,0,s1))/(s0*s1*2.*TMath::Pi());
 717     return val;
 718   }
 719
 720   Double_t sigma2 = k1*k1*s0*s0+k0*k0*s1*s1;
 721   Double_t exp0 = TMath::Exp(-(k1*x0-k0*x1)*(k1*x0-k0*x1)/(2*sigma2));
 722   //
 723   Double_t sigmaErf =  2*s0*s1*TMath::Sqrt(2*sigma2);
 724   Double_t erf0 = TMath::Erf( (k0*s1*s1*(k0-2*x0)+k1*s0*s0*(k1-2*x1))/sigmaErf);
 725   Double_t erf1 = TMath::Erf( (k0*s1*s1*(k0+2*x0)+k1*s0*s0*(k1+2*x1))/sigmaErf);
 726   Double_t norm = 1./TMath::Sqrt(sigma2);
 727   norm/=2.*TMath::Sqrt(2.*TMath::Pi());
 728   Double_t val  = norm*exp0*(erf0+erf1);
 729   return val;
 730
 731 }
 732
 733
 734 Double_t  AliTPCclusterFast::GaussExpConvolution(Double_t x0, Double_t s0,Double_t t1){
 735  //
 736   // 2 D gaus convoluted with exponential
 737   // Integral nomalized to 1
 738   // See in mathematica:
 739   //Simplify[Integrate[Exp[-(x0-x1)*(x0-x1)/(2*s0*s0)]*Exp[-x1*t1],{x1,0,Infinity}]]
 740   // TF1 fgexp("fgexp","AliTPCclusterFast::GaussExpConvolution(x,0.5,1)",-2,2)
 741   Double_t exp1 = (s0*s0*t1-2*x0)*t1/2.;
 742   exp1 = TMath::Exp(exp1);
 743   Double_t erf = 1+TMath::Erf((-s0*s0*t1+x0)/(s0*TMath::Sqrt(2.)));
 744   Double_t val = exp1*erf;
 745   val *=t1/(2.);
 746   return val;
 747
 748 }
 749
 750
 751 Double_t  AliTPCclusterFast::Gamma4(Double_t x, Double_t p0, Double_t p1){
 752   //
 753   // Gamma 4 Time response function of ALTRO
 754   //
 755   if (x<0) return 0;
 756   Double_t g1 = TMath::Exp(-4.*x/p1);
 757   Double_t g2 = TMath::Power(x/p1,4);
 758   return p0*g1*g2;
 759 }
 760
 761
 762
 763 Double_t  AliTPCclusterFast::GaussGamma4(Double_t x, Double_t s0, Double_t p1){
 764   //
 765   // Gamma 4 Time response function of ALTRO convoluted with Gauss
 766   // Simplify[Integrate[Exp[-(x0-x1)*(x0-x1)/(2*s0*s0)]*Exp[-4*x1/p1]*(x/p1)^4/s0,{x1,0,Infinity}]]
 767   //TF1 fgg4("fgg4","AliTPCclusterFast::GaussGamma4(x,0.5,0.5)",-2,2)
 768
 769   Double_t exp1 = (8*s0*s0-4.*p1*x)/(p1*p1);
 770   exp1 = TMath::Exp(exp1);
 771   Double_t erf1 = 1+TMath::Erf((-4*s0/p1+x/s0)/TMath::Sqrt(2));
 772   //  Double_t xp14 = TMath::Power(TMath::Abs((x/p1)),4);
 773   return exp1*erf1;
 774
 775
 776 }
 777
 778 // Analytical sollution only in 1D - too long expression
 779 // Simplify[Integrate[Exp[-(x0-(x1-k*x2))*(x0-(x1-k*x2))/(2*s0*s0)]*Exp[-(x1*t1-k*x2)],{x2,-1,1}]]
 780 //
 781 //
 782 // No analytical solution
 783 //
 784 //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}]]