a55a6748118828b058208a790d42280cbbd7da9a
[u/mrichter/AliRoot.git] / 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}]]