ATO-18 - a.) Removed upper limit on number of electrons. b.)Use qMax correction if...
[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 #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}]]