Merge remote-tracking branch 'remotes/origin/master' into TPCdev
[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 #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}]]