TPC/fastSimul/AliTPCclusterFast.cxx

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