kBoth + R_photon > 35 cm no weights
[u/mrichter/AliRoot.git] / TPC / fastSimul / AliTPCclusterFast.cxx
CommitLineData
0716cb63 1/*
2 gSystem->Load("libSTAT.so");
37045aeb 3 .x ~/NimStyle.C
0716cb63 4
5 .L $ALICE_ROOT/TPC/fastSimul/AliTPCclusterFast.cxx+
37045aeb 6 //
13f45dd8 7 AliTPCclusterFast::fPRF = new TF1("fprf","gausn",-5,5);
8 AliTPCclusterFast::fTRF = new TF1("ftrf","gausn",-5,5);
37045aeb 9 AliTPCclusterFast::fPRF->SetParameters(1,0,0.5);
10 AliTPCclusterFast::fTRF->SetParameters(1,0,0.5);
0716cb63 11
13f45dd8 12 //
46a0eaea 13 AliTPCtrackFast::Simul("trackerSimul.root",100);
14// AliTPCclusterFast::Simul("cluterSimul.root",20000);
0716cb63 15*/
16
17#include "TObject.h"
37045aeb 18#include "TF1.h"
0716cb63 19#include "TMath.h"
20#include "TRandom.h"
21#include "TVectorD.h"
37045aeb 22#include "TMatrixD.h"
13f45dd8 23#include "TH1.h"
13f45dd8 24#include "TClonesArray.h"
0716cb63 25#include "TTreeStream.h"
26
27class AliTPCclusterFast: public TObject {
28public:
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();
37045aeb 33 void Digitize();
13f45dd8 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
37045aeb 39 Double_t GetNsec();
0716cb63 40 static void Simul(const char* simul, Int_t npoints);
13f45dd8 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);
0716cb63 45public:
46 Float_t fMNprim; // mean number of primary electrons
47 // //electrons part input
37045aeb 48 Int_t fNprim; // mean number of primary electrons
13f45dd8 49 Int_t fNtot; // total number of electrons
50 Float_t fQtot; // total charge - Gas gain flucuation taken into account
51 //
0716cb63 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
37045aeb 57 //
58 //
0716cb63 59 // // electron part simul
13f45dd8 60 TVectorD fSec; //! number of secondary electrons
0716cb63 61 TVectorD fPosY; //! position y for each electron
62 TVectorD fPosZ; //! position z for each electron
63 TVectorD fGain; //! gg for each electron
37045aeb 64 //
13f45dd8 65 TVectorD fStatY; //!stat Y
66 TVectorD fStatZ; //!stat Y
0716cb63 67 //
37045aeb 68 // digitization part
69 //
70 TMatrixD fDigits; // response matrix
71 static TF1* fPRF; // Pad response
72 static TF1* fTRF; // Time response function
0716cb63 73 ClassDef(AliTPCclusterFast,1) // container for
74};
13f45dd8 75
76
77class AliTPCtrackFast: public TObject {
78public:
79 AliTPCtrackFast();
46a0eaea 80 void Add(AliTPCtrackFast &track2);
13f45dd8 81 void MakeTrack();
13f45dd8 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);
46a0eaea 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);
13f45dd8 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
46a0eaea 99 //
100 //
13f45dd8 101 ClassDef(AliTPCtrackFast,2) // container for
102};
103
104
105
0716cb63 106ClassImp(AliTPCclusterFast)
13f45dd8 107ClassImp(AliTPCtrackFast)
108
109
110
111
0716cb63 112
37045aeb 113TF1 *AliTPCclusterFast::fPRF=0;
114TF1 *AliTPCclusterFast::fTRF=0;
0716cb63 115
116
13f45dd8 117AliTPCtrackFast::AliTPCtrackFast():
118 TObject(),
119 fMNprim(0),
120 fAngleY(0),
121 fAngleZ(0),
122 fN(0),
123 fCl(0),
44a0fef0 124 fInit(kFALSE)
13f45dd8 125{
126 //
127 //
128 //
46a0eaea 129}
130
131void AliTPCtrackFast::Add(AliTPCtrackFast &track2){
132 if (!track2.fInit) return;
133
13f45dd8 134}
135
13f45dd8 136
13f45dd8 137
138void AliTPCtrackFast::MakeTrack(){
139 //
140 //
141 //
46a0eaea 142 if (!fCl) fCl = new TClonesArray("AliTPCclusterFast",160);
13f45dd8 143 for (Int_t i=0;i<fN;i++){
46a0eaea 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 //
13f45dd8 153 cluster->GenerElectrons();
154 cluster->Digitize();
155 }
13f45dd8 156}
157
158Double_t AliTPCtrackFast::CookdEdxNtot(Double_t f0,Float_t f1){
159 //
160 Double_t amp[160];
13f45dd8 161 for (Int_t i=0;i<fN;i++){
162 AliTPCclusterFast * cluster = ( AliTPCclusterFast *)((*fCl)[i]);
163 amp[i]=cluster->fNtot;
164 }
46a0eaea 165 return CookdEdx(fN,amp,f0,f1);
13f45dd8 166}
167
168Double_t AliTPCtrackFast::CookdEdxQtot(Double_t f0,Float_t f1){
169 //
170 Double_t amp[160];
13f45dd8 171 for (Int_t i=0;i<fN;i++){
172 AliTPCclusterFast * cluster = ( AliTPCclusterFast *)((*fCl)[i]);
173 amp[i]=cluster->fQtot;
174 }
46a0eaea 175 return CookdEdx(fN,amp,f0,f1);
13f45dd8 176}
177
46a0eaea 178Double_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
197Double_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
217Double_t AliTPCtrackFast::CookdEdx(Int_t npoints, Double_t *amp,Double_t f0,Float_t f1){
13f45dd8 218 //
219 //
220 //
221 Int_t index[160];
46a0eaea 222 TMath::Sort(npoints,amp,index,kFALSE);
13f45dd8 223 Float_t sum0=0, sum1=0,sum2=0;
46a0eaea 224 for (Int_t i=0;i<npoints;i++){
225 if (i<npoints*f0) continue;
226 if (i>npoints*f1) continue;
13f45dd8 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
235void AliTPCtrackFast::Simul(const char* fname, Int_t ntracks){
236 //
237 //
238 //
239 AliTPCtrackFast fast;
44a0fef0 240 TTreeSRedirector cstream(fname,"recreate");
13f45dd8 241 for (Int_t itr=0; itr<ntracks; itr++){
242 //
c9e8d4ec 243 fast.fMNprim=(5+50*gRandom->Rndm());
13f45dd8 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);
46a0eaea 248 fast.fN = TMath::Nint(80.+gRandom->Rndm()*80.);
13f45dd8 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
0716cb63 260AliTPCclusterFast::AliTPCclusterFast(){
37045aeb 261 //
262 //
263 fDigits.ResizeTo(5,7);
0716cb63 264}
37045aeb 265
0716cb63 266AliTPCclusterFast::~AliTPCclusterFast(){
267}
268
269
270void 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}
37045aeb 277Double_t AliTPCclusterFast::GetNsec(){
0716cb63 278 //
37045aeb 279 // Generate number of secondary electrons
280 // copy of procedure implemented in geant
0716cb63 281 //
0716cb63 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;
37045aeb 285 Float_t RAN = gRandom->Rndm();
46a0eaea 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);
37045aeb 289 return TMath::Nint(TMath::Power((TMath::Power(FPOT,XEXPO)*(1-RAN)+TMath::Power(EEND,XEXPO)*RAN),YEXPO)/W);
290}
291
292void AliTPCclusterFast::GenerElectrons(){
0716cb63 293 //
37045aeb 294 //
295 //
296 //
297 const Int_t knMax=1000;
0716cb63 298 if (fPosY.GetNrows()<knMax){
299 fPosY.ResizeTo(knMax);
300 fPosZ.ResizeTo(knMax);
301 fGain.ResizeTo(knMax);
37045aeb 302 fSec.ResizeTo(knMax);
0716cb63 303 fStatY.ResizeTo(3);
304 fStatZ.ResizeTo(3);
305 }
37045aeb 306 fNprim = gRandom->Poisson(fMNprim); //number of primary electrons
13f45dd8 307 fNtot=0; //total number of electrons
308 fQtot=0; //total number of electrons after gain multiplification
0716cb63 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;
37045aeb 315 for (Int_t i=0;i<knMax;i++){
316 fSec[i]=0;
317 }
0716cb63 318 for (Int_t iprim=0; iprim<fNprim;iprim++){
37045aeb 319 Float_t dN = GetNsec();
320 fSec[iprim]=dN;
0716cb63 321 Double_t yc = fY+(gRandom->Rndm()-0.5)*fAngleY;
322 Double_t zc = fZ+(gRandom->Rndm()-0.5)*fAngleZ;
37045aeb 323 for (Int_t isec=0;isec<=dN;isec++){
0716cb63 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;
13f45dd8 332 fQtot+=gg;
0716cb63 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
37045aeb 353void 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++)
13f45dd8 366 for (Int_t dj=-3; dj<=3;dj++){
37045aeb 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 }
13f45dd8 372
37045aeb 373}
374
375
0716cb63 376
377void AliTPCclusterFast::Simul(const char* fname, Int_t npoints){
13f45dd8 378 //
379 // Calc rms
380 //
0716cb63 381 AliTPCclusterFast fast;
382 TTreeSRedirector cstream(fname);
383 for (Int_t icl=0; icl<npoints; icl++){
13f45dd8 384 Float_t nprim=(10+20*gRandom->Rndm());
a314e83c 385 Float_t diff =0.01 +0.35*gRandom->Rndm();
0716cb63 386 Float_t posY = gRandom->Rndm()-0.5;
387 Float_t posZ = gRandom->Rndm()-0.5;
13f45dd8 388 //
389 Float_t ky = 4.0*(gRandom->Rndm()-0.5);
390 Float_t kz = 4.0*(gRandom->Rndm()-0.5);
0716cb63 391 fast.SetParam(nprim,diff,posY,posZ,ky,kz);
392 fast.GenerElectrons();
37045aeb 393 fast.Digitize();
13f45dd8 394 if (icl%10000==0) printf("%d\n",icl);
0716cb63 395 cstream<<"simul"<<
396 "s.="<<&fast<<
397 "\n";
398 }
399}
400
401
13f45dd8 402Double_t AliTPCclusterFast::GetQtot(Float_t gain, Float_t thr, Float_t noise, Bool_t brounding, Bool_t baddPedestal){
a314e83c 403 //
404 //
405 //
406 Float_t sum =0;
13f45dd8 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);
a314e83c 413 if (amp>thr) sum+=amp;
414 }
13f45dd8 415 }
a314e83c 416 return sum;
417}
418
13f45dd8 419Double_t AliTPCclusterFast::GetQmax(Float_t gain, Float_t thr, Float_t noise, Bool_t brounding, Bool_t baddPedestal){
a314e83c 420 //
421 //
422 //
423 Float_t max =0;
13f45dd8 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;
a314e83c 431 }
13f45dd8 432 }
a314e83c 433 return max;
434}
435
436
a314e83c 437
13f45dd8 438Double_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
453Double_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}
a314e83c 480
481
0716cb63 482
13f45dd8 483
484
485Double_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
515Double_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
532Double_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
544Double_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));
46a0eaea 553 // Double_t xp14 = TMath::Power(TMath::Abs((x/p1)),4);
13f45dd8 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}]]