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