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