kBoth + R_photon > 35 cm no weights
[u/mrichter/AliRoot.git] / TPC / fastSimul / AliTPCclusterFast.cxx
1 /*
2   gSystem->Load("libSTAT.so");
3   .x ~/NimStyle.C
4  
5  .L $ALICE_ROOT/TPC/fastSimul/AliTPCclusterFast.cxx+
6  //
7  AliTPCclusterFast::fPRF = new TF1("fprf","gausn",-5,5);
8  AliTPCclusterFast::fTRF = new TF1("ftrf","gausn",-5,5);
9  AliTPCclusterFast::fPRF->SetParameters(1,0,0.5);
10  AliTPCclusterFast::fTRF->SetParameters(1,0,0.5);
11
12  //
13  AliTPCtrackFast::Simul("trackerSimul.root",100); 
14 // AliTPCclusterFast::Simul("cluterSimul.root",20000); 
15 */
16
17 #include "TObject.h"
18 #include "TF1.h"
19 #include "TMath.h"
20 #include "TRandom.h"
21 #include "TVectorD.h"
22 #include "TMatrixD.h"
23 #include "TH1.h"
24 #include "TClonesArray.h"
25 #include "TTreeStream.h"
26
27 class AliTPCclusterFast: public TObject {
28 public:
29   AliTPCclusterFast();
30   virtual ~AliTPCclusterFast();
31   void SetParam(Float_t mnprim, Float_t diff, Float_t y, Float_t z, Float_t ky, Float_t kz);
32   void GenerElectrons();
33   void Digitize();
34   Double_t GetQtot(Float_t gain,Float_t thr, Float_t noise, Bool_t rounding=kTRUE, Bool_t addPedestal=kTRUE);
35   Double_t GetQmax(Float_t gain,Float_t thr, Float_t noise, Bool_t rounding=kTRUE, Bool_t addPedestal=kTRUE);
36   Double_t GetQmaxCorr(Float_t rmsy0, Float_t rmsz0);
37   Double_t GetQtotCorr(Float_t rmsy0, Float_t rmsz0, Float_t gain, Float_t thr);
38   
39   Double_t GetNsec();
40   static void Simul(const char* simul, Int_t npoints);
41   static Double_t GaussConvolution(Double_t x0, Double_t x1, Double_t k0, Double_t k1, Double_t s0, Double_t s1);
42   static Double_t GaussExpConvolution(Double_t x0, Double_t s0,Double_t t1);
43   static Double_t GaussGamma4(Double_t x, Double_t s0, Double_t p1);
44   static Double_t Gamma4(Double_t x, Double_t p0, Double_t p1);
45 public:
46   Float_t fMNprim;     // mean number of primary electrons
47   //                   //electrons part input
48   Int_t   fNprim;      // mean number of primary electrons
49   Int_t   fNtot;       // total number of  electrons
50   Float_t fQtot;       // total charge - Gas gain flucuation taken into account
51   //
52   Float_t fDiff;       // diffusion sigma
53   Float_t fY;          // y position 
54   Float_t fZ;          // z postion 
55   Float_t fAngleY;     // y angle - tan(y)
56   Float_t fAngleZ;     // z angle - tan z
57   //
58   //
59   //                   // electron part simul
60   TVectorD fSec;       //! number of secondary electrons
61   TVectorD fPosY;      //! position y for each electron
62   TVectorD fPosZ;      //! position z for each electron
63   TVectorD fGain;      //! gg for each electron
64   //
65   TVectorD fStatY;     //!stat Y  
66   TVectorD fStatZ;     //!stat Y
67   //
68   // digitization part
69   //
70   TMatrixD fDigits;    // response matrix
71   static TF1* fPRF;    // Pad response
72   static TF1* fTRF;    // Time response function 
73   ClassDef(AliTPCclusterFast,1)  // container for
74 };
75
76
77 class AliTPCtrackFast: public TObject {
78 public:
79   AliTPCtrackFast();
80   void Add(AliTPCtrackFast &track2);
81   void MakeTrack();
82   static void Simul(const char* simul, Int_t ntracks);
83   Double_t  CookdEdxNtot(Double_t f0,Float_t f1);
84   Double_t  CookdEdxQtot(Double_t f0,Float_t f1);
85   //
86   Double_t  CookdEdxDtot(Double_t f0,Float_t f1, Float_t gain,Float_t thr, Float_t noise, Bool_t corr = kTRUE);
87   Double_t  CookdEdxDmax(Double_t f0,Float_t f1,Float_t gain,Float_t thr, Float_t noise, Bool_t corr=kTRUE);
88   //
89   Double_t  CookdEdx(Int_t npoints, Double_t *amp, Double_t f0,Float_t f1);
90   //
91   Float_t fMNprim;     // mean number of primary electrons
92   Float_t fAngleY;     // y angle - tan(y)
93   Float_t fAngleZ;     // z angle - tan z
94   Float_t fDiff;       // diffusion
95   Int_t   fN;          // number of clusters
96   TClonesArray *fCl;   // array of clusters  
97   //
98   Bool_t   fInit;      // initialization flag
99   //
100   //
101   ClassDef(AliTPCtrackFast,2)  // container for
102 };
103
104
105
106 ClassImp(AliTPCclusterFast)
107 ClassImp(AliTPCtrackFast)
108
109
110
111
112
113 TF1 *AliTPCclusterFast::fPRF=0;
114 TF1 *AliTPCclusterFast::fTRF=0;
115
116
117 AliTPCtrackFast::AliTPCtrackFast():
118   TObject(),
119   fMNprim(0),
120   fAngleY(0),
121   fAngleZ(0),
122   fN(0),
123   fCl(0),
124   fInit(kFALSE)
125 {
126   //
127   //
128   //
129 }
130
131 void AliTPCtrackFast::Add(AliTPCtrackFast &track2){
132   if (!track2.fInit) return;
133   
134 }
135
136
137
138 void AliTPCtrackFast::MakeTrack(){
139   //
140   //
141   //
142   if (!fCl) fCl = new TClonesArray("AliTPCclusterFast",160);
143   for (Int_t i=0;i<fN;i++){
144     Double_t tY = i*fAngleY;
145     Double_t tZ = i*fAngleZ;
146     AliTPCclusterFast * cluster = (AliTPCclusterFast*) fCl->UncheckedAt(i);
147     if (!cluster) cluster =   new ((*fCl)[i]) AliTPCclusterFast;
148     //
149     Double_t posY = tY-TMath::Nint(tY);
150     Double_t posZ = tZ-TMath::Nint(tZ);
151     cluster->SetParam(fMNprim,fDiff,posY,posZ,fAngleY,fAngleZ); 
152     //
153     cluster->GenerElectrons();
154     cluster->Digitize();
155   }
156 }
157
158 Double_t  AliTPCtrackFast::CookdEdxNtot(Double_t f0,Float_t f1){
159   //
160   Double_t amp[160];
161   for (Int_t i=0;i<fN;i++){ 
162     AliTPCclusterFast * cluster = ( AliTPCclusterFast *)((*fCl)[i]);
163     amp[i]=cluster->fNtot;
164   }
165   return CookdEdx(fN,amp,f0,f1);
166 }
167
168 Double_t  AliTPCtrackFast::CookdEdxQtot(Double_t f0,Float_t f1){
169   //
170   Double_t amp[160];
171   for (Int_t i=0;i<fN;i++){ 
172     AliTPCclusterFast * cluster = ( AliTPCclusterFast *)((*fCl)[i]);
173     amp[i]=cluster->fQtot;
174   }
175   return CookdEdx(fN,amp,f0,f1);
176 }
177
178 Double_t   AliTPCtrackFast::CookdEdxDtot(Double_t f0,Float_t f1, Float_t gain,Float_t thr, Float_t noise, Bool_t doCorr){
179   //
180   //
181   //
182   Double_t amp[160];
183   Int_t over=0;
184   for (Int_t i=0;i<fN;i++){ 
185     AliTPCclusterFast * cluster = ( AliTPCclusterFast *)((*fCl)[i]);
186     Float_t camp = cluster->GetQtot(gain,thr,noise);
187     if (camp==0) continue;
188     Float_t corr =  1;
189     if (doCorr) corr = cluster->GetQtotCorr(0.5,0.5,gain,thr);
190     amp[over]=camp/corr;
191     over++;
192   }
193   return CookdEdx(over,amp,f0,f1);
194
195 }
196
197 Double_t   AliTPCtrackFast::CookdEdxDmax(Double_t f0,Float_t f1, Float_t gain,Float_t thr, Float_t noise, Bool_t doCorr){
198   //
199   //
200   //
201   Double_t amp[160];
202   Int_t over=0;
203   for (Int_t i=0;i<fN;i++){ 
204     AliTPCclusterFast * cluster = ( AliTPCclusterFast *)((*fCl)[i]);
205     Float_t camp = cluster->GetQmax(gain,thr,noise);
206     if (camp==0) continue;    
207     Float_t corr =  1;
208     if (doCorr) corr = cluster->GetQmaxCorr(0.5,0.5);
209     amp[over]=camp/corr;
210     over++;
211   }
212   return CookdEdx(over,amp,f0,f1);
213
214 }
215
216
217 Double_t  AliTPCtrackFast::CookdEdx(Int_t npoints, Double_t *amp,Double_t f0,Float_t f1){
218   //
219   //
220   //
221   Int_t index[160];
222   TMath::Sort(npoints,amp,index,kFALSE);
223   Float_t sum0=0, sum1=0,sum2=0;
224   for (Int_t i=0;i<npoints;i++){
225     if (i<npoints*f0) continue;
226     if (i>npoints*f1) continue;
227     sum0++;
228     sum1+= amp[index[i]];
229     sum2+= amp[index[i]];
230   }
231   if (sum0<=0) return 0;
232   return sum1/sum0;
233 }
234
235 void AliTPCtrackFast::Simul(const char* fname, Int_t ntracks){
236   //
237   // 
238   //
239   AliTPCtrackFast fast;
240   TTreeSRedirector cstream(fname,"recreate");
241   for (Int_t itr=0; itr<ntracks; itr++){
242     //
243     fast.fMNprim=(5+50*gRandom->Rndm());
244     fast.fDiff =0.01 +0.35*gRandom->Rndm();
245     //
246     fast.fAngleY   = 4.0*(gRandom->Rndm()-0.5);
247     fast.fAngleZ   = 4.0*(gRandom->Rndm()-0.5);
248     fast.fN  = TMath::Nint(80.+gRandom->Rndm()*80.);
249     fast.MakeTrack();
250     if (itr%100==0) printf("%d\n",itr);
251     cstream<<"simulTrack"<<
252       "tr.="<<&fast<<
253       "\n";
254   }
255   fast.Write("track");
256 }
257
258
259
260 AliTPCclusterFast::AliTPCclusterFast(){
261   //
262   //
263   fDigits.ResizeTo(5,7);
264 }
265
266 AliTPCclusterFast::~AliTPCclusterFast(){
267 }
268
269
270 void AliTPCclusterFast::SetParam(Float_t mnprim, Float_t diff, Float_t y, Float_t z, Float_t ky, Float_t kz){
271   //
272   //
273   fMNprim = mnprim; fDiff = diff;
274   fY=y; fZ=z; 
275   fAngleY=ky; fAngleZ=kz;
276 }
277 Double_t AliTPCclusterFast::GetNsec(){
278   //
279   // Generate number of secondary electrons
280   // copy of procedure implemented in geant
281   //
282   const Double_t FPOT=20.77E-9, EEND=10E-6, EEXPO=2.2, EEND1=1E-6;
283   const Double_t XEXPO=-EEXPO+1, YEXPO=1/XEXPO;
284   const Double_t W=20.77E-9;
285   Float_t RAN = gRandom->Rndm();
286   //Double_t edep = TMath::Power((TMath::Power(FPOT,XEXPO)*(1-RAN)+TMath::Power(EEND,XEXPO)*RAN),YEXPO);
287   //edep = TMath::Min(edep, EEND);
288   //return TMath::Nint(edep/W);
289   return TMath::Nint(TMath::Power((TMath::Power(FPOT,XEXPO)*(1-RAN)+TMath::Power(EEND,XEXPO)*RAN),YEXPO)/W);
290 }
291
292 void AliTPCclusterFast::GenerElectrons(){
293   //
294   //
295   //
296   //
297   const Int_t knMax=1000;
298   if (fPosY.GetNrows()<knMax){
299     fPosY.ResizeTo(knMax);
300     fPosZ.ResizeTo(knMax);
301     fGain.ResizeTo(knMax);
302     fSec.ResizeTo(knMax);
303     fStatY.ResizeTo(3);
304     fStatZ.ResizeTo(3);
305   }
306   fNprim = gRandom->Poisson(fMNprim);  //number of primary electrons
307   fNtot=0; //total number of electrons
308   fQtot=0; //total number of electrons after gain multiplification
309   //
310   Double_t sumQ=0;
311   Double_t sumYQ=0;
312   Double_t sumZQ=0;
313   Double_t sumY2Q=0;
314   Double_t sumZ2Q=0;
315   for (Int_t i=0;i<knMax;i++){ 
316     fSec[i]=0;
317   }
318   for (Int_t iprim=0; iprim<fNprim;iprim++){
319     Float_t dN   =  GetNsec();
320     fSec[iprim]=dN;
321     Double_t yc = fY+(gRandom->Rndm()-0.5)*fAngleY;
322     Double_t zc = fZ+(gRandom->Rndm()-0.5)*fAngleZ;
323     for (Int_t isec=0;isec<=dN;isec++){
324       //
325       //
326       Double_t y = gRandom->Gaus(0,fDiff)+yc;
327       Double_t z = gRandom->Gaus(0,fDiff)+zc;
328       Double_t gg = -TMath::Log(gRandom->Rndm());
329       fPosY[fNtot]=y;
330       fPosZ[fNtot]=z;
331       fGain[fNtot]=gg;
332       fQtot+=gg;
333       fNtot++;
334       sumQ+=gg;
335       sumYQ+=gg*y;
336       sumY2Q+=gg*y*y;
337       sumZQ+=gg*z;
338       sumZ2Q+=gg*z*z;
339       if (fNtot>=knMax) break;
340     }
341     if (fNtot>=knMax) break;
342   }
343   if (sumQ>0){
344     fStatY[0]=sumQ;
345     fStatY[1]=sumYQ/sumQ;
346     fStatY[2]=sumY2Q/sumQ-fStatY[1]*fStatY[1];
347     fStatZ[0]=sumQ;
348     fStatZ[1]=sumZQ/sumQ;
349     fStatZ[2]=sumZ2Q/sumQ-fStatZ[1]*fStatZ[1];
350   }
351 }
352
353 void AliTPCclusterFast::Digitize(){
354   //
355   //
356   //
357   // 1. Clear digits
358   for (Int_t i=0; i<5;i++)
359     for (Int_t j=0; j<7;j++){
360       fDigits(i,j)=0;
361     }
362   //
363   // Fill digits
364   for (Int_t iel = 0; iel<fNtot; iel++){
365     for (Int_t di=-2; di<=2;di++)
366       for (Int_t dj=-3; dj<=3;dj++){
367         Float_t fac = fPRF->Eval(di-fPosY[iel])*fTRF->Eval(dj-fPosZ[iel]);
368         fac*=fGain[iel];
369         fDigits(2+di,3+dj)+=fac;
370       }
371   }
372   
373 }
374
375
376
377 void AliTPCclusterFast::Simul(const char* fname, Int_t npoints){
378   //
379   // Calc rms
380   //
381   AliTPCclusterFast fast;
382   TTreeSRedirector cstream(fname);
383   for (Int_t icl=0; icl<npoints; icl++){
384     Float_t nprim=(10+20*gRandom->Rndm());
385     Float_t diff =0.01 +0.35*gRandom->Rndm();
386     Float_t posY = gRandom->Rndm()-0.5;
387     Float_t posZ = gRandom->Rndm()-0.5;
388     //
389     Float_t ky   = 4.0*(gRandom->Rndm()-0.5);
390     Float_t kz   = 4.0*(gRandom->Rndm()-0.5);
391     fast.SetParam(nprim,diff,posY,posZ,ky,kz);
392     fast.GenerElectrons();
393     fast.Digitize();
394     if (icl%10000==0) printf("%d\n",icl);
395     cstream<<"simul"<<
396       "s.="<<&fast<<
397       "\n";
398   }
399 }
400
401
402 Double_t AliTPCclusterFast::GetQtot(Float_t gain, Float_t thr, Float_t noise, Bool_t brounding, Bool_t baddPedestal){
403   //
404   //
405   //
406   Float_t sum =0;
407   for (Int_t ip=0;ip<5;ip++){
408     Float_t pedestal=gRandom->Rndm()-0.5; //pedestal offset different for each pad
409     for (Int_t it=0;it<7;it++){
410       Float_t amp = gain*fDigits(ip,it)+gRandom->Gaus()*noise;
411       if (baddPedestal) amp+=pedestal;
412       if (brounding) amp=TMath::Nint(amp);
413       if (amp>thr) sum+=amp;
414     } 
415   }
416   return sum;
417 }
418
419 Double_t AliTPCclusterFast::GetQmax(Float_t gain, Float_t thr, Float_t noise, Bool_t brounding, Bool_t baddPedestal){
420   //
421   //
422   //
423   Float_t max =0;
424   for (Int_t ip=0;ip<5;ip++){
425     Float_t pedestal=gRandom->Rndm()-0.5; //pedestal offset different for each pad
426     for (Int_t it=0;it<7;it++){
427       Float_t amp = gain*fDigits(ip,it)+gRandom->Gaus()*noise;
428       if (baddPedestal) amp+=pedestal;
429       if (brounding) amp=TMath::Nint(amp);
430       if (amp>max && amp>thr) max=amp;
431     } 
432   }
433   return max;
434 }
435
436
437
438 Double_t  AliTPCclusterFast::GetQmaxCorr(Float_t rmsy0, Float_t rmsz0){
439   //
440   // Gaus distribution convolueted with rectangular
441   // Gaus width sy and sz is determined by RF width and diffusion 
442   // Integral of Q is equal 1
443   // Q max is calculated at position fY,fX 
444   //          
445   //  
446   //  
447   Double_t sy = TMath::Sqrt(rmsy0*rmsy0+fDiff*fDiff);
448   Double_t sz = TMath::Sqrt(rmsz0*rmsz0+fDiff*fDiff); 
449   return GaussConvolution(fY,fZ, fAngleY,fAngleZ,sy,sz);
450 }
451
452
453 Double_t  AliTPCclusterFast::GetQtotCorr(Float_t rmsy0, Float_t rmsz0, Float_t gain, Float_t thr){
454   //
455   //  Calculates the fraction of the charge over threshol to total charge
456   //  The response function
457   //
458   Double_t sy = TMath::Sqrt(rmsy0*rmsy0+fDiff*fDiff);
459   Double_t sz = TMath::Sqrt(rmsz0*rmsz0+fDiff*fDiff);
460   Double_t sumAll=0,sumThr=0;
461   Double_t qtot = GetQtot(gain,thr,0); // sum of signal over threshold
462   //
463   Double_t corr =1;
464   Double_t qnorm=qtot;
465   for (Int_t iter=0;iter<2;iter++){
466     for (Int_t iy=-2;iy<=2;iy++)
467       for (Int_t iz=-2;iz<=2;iz++){      
468         Double_t val = GaussConvolution(fY-iy,fZ-iz, fAngleY,fAngleZ,sy,sz);
469         Double_t qlocal =TMath::Nint(qnorm*val);
470         if (qlocal>thr) sumThr+=qlocal;
471         sumAll+=qlocal;
472       }
473     if (sumAll>0&&sumThr>0) corr=(sumThr)/sumAll;
474     //corr = sumThr;
475     if (corr>0) qnorm=qtot/corr;
476     
477   }
478   return corr;
479 }
480
481
482
483
484
485 Double_t  AliTPCclusterFast::GaussConvolution(Double_t x0, Double_t x1, Double_t k0, Double_t k1, Double_t s0, Double_t s1){
486   //
487   // 2 D gaus convoluted with angular effect
488   // See in mathematica: 
489   //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}]]
490   // 
491   //TF1 f1("f1","AliTPCclusterFast::GaussConvolution(x,0,1,0,0.1,0.1)",-2,2)
492   //TF2 f2("f2","AliTPCclusterFast::GaussConvolution(x,y,1,1,0.1,0.1)",-2,2,-2,2)
493   //
494   const Float_t kEpsilon = 0.0001;
495   if ((TMath::Abs(k0)+TMath::Abs(k1))<kEpsilon*(s0+s1)){
496     // small angular effect
497     Double_t val = (TMath::Gaus(x0,0,s0)*TMath::Gaus(x1,0,s1))/(s0*s1*2.*TMath::Pi());
498     return val;
499   }
500   
501   Double_t sigma2 = k1*k1*s0*s0+k0*k0*s1*s1;
502   Double_t exp0 = TMath::Exp(-(k1*x0-k0*x1)*(k1*x0-k0*x1)/(2*sigma2));
503   //
504   Double_t sigmaErf =  2*s0*s1*TMath::Sqrt(2*sigma2);                        
505   Double_t erf0 = TMath::Erf( (k0*s1*s1*(k0-2*x0)+k1*s0*s0*(k1-2*x1))/sigmaErf);
506   Double_t erf1 = TMath::Erf( (k0*s1*s1*(k0+2*x0)+k1*s0*s0*(k1+2*x1))/sigmaErf);
507   Double_t norm = 1./TMath::Sqrt(sigma2);
508   norm/=2.*TMath::Sqrt(2.*TMath::Pi());
509   Double_t val  = norm*exp0*(erf0+erf1);
510   return val;
511
512 }
513
514
515 Double_t  AliTPCclusterFast::GaussExpConvolution(Double_t x0, Double_t s0,Double_t t1){
516  //
517   // 2 D gaus convoluted with exponential
518   // Integral nomalized to 1
519   // See in mathematica: 
520   //Simplify[Integrate[Exp[-(x0-x1)*(x0-x1)/(2*s0*s0)]*Exp[-x1*t1],{x1,0,Infinity}]]
521   // TF1 fgexp("fgexp","AliTPCclusterFast::GaussExpConvolution(x,0.5,1)",-2,2)
522   Double_t exp1 = (s0*s0*t1-2*x0)*t1/2.;
523   exp1 = TMath::Exp(exp1);
524   Double_t erf = 1+TMath::Erf((-s0*s0*t1+x0)/(s0*TMath::Sqrt(2.)));
525   Double_t val = exp1*erf;
526   val *=t1/(2.);
527   return val;
528
529 }
530
531
532 Double_t  AliTPCclusterFast::Gamma4(Double_t x, Double_t p0, Double_t p1){
533   //
534   // Gamma 4 Time response function of ALTRO
535   //
536   if (x<0) return 0;
537   Double_t g1 = TMath::Exp(-4.*x/p1);
538   Double_t g2 = TMath::Power(x/p1,4);
539   return p0*g1*g2;
540 }
541  
542
543
544 Double_t  AliTPCclusterFast::GaussGamma4(Double_t x, Double_t s0, Double_t p1){
545   //
546   // Gamma 4 Time response function of ALTRO convoluted with Gauss
547   // Simplify[Integrate[Exp[-(x0-x1)*(x0-x1)/(2*s0*s0)]*Exp[-4*x1/p1]*(x/p1)^4/s0,{x1,0,Infinity}]]
548   //TF1 fgg4("fgg4","AliTPCclusterFast::GaussGamma4(x,0.5,0.5)",-2,2)
549
550   Double_t exp1 = (8*s0*s0-4.*p1*x)/(p1*p1);
551   exp1 = TMath::Exp(exp1);
552   Double_t erf1 = 1+TMath::Erf((-4*s0/p1+x/s0)/TMath::Sqrt(2));
553   //  Double_t xp14 = TMath::Power(TMath::Abs((x/p1)),4);
554   return exp1*erf1;
555
556  
557 }
558  
559 // Analytical sollution only in 1D - too long expression
560 // Simplify[Integrate[Exp[-(x0-(x1-k*x2))*(x0-(x1-k*x2))/(2*s0*s0)]*Exp[-(x1*t1-k*x2)],{x2,-1,1}]] 
561 //
562 //
563 // No analytical solution
564 // 
565 //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}]]