6 gSystem->Load("libSTAT.so");
9 .L $ALICE_ROOT/TPC/fastSimul/AliTPCclusterFast.cxx+
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);
17 AliTPCtrackFast::Simul("trackerSimul.root",100);
18 // AliTPCclusterFast::Simul("cluterSimul.root",20000);
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.
44 #include "TClonesArray.h"
45 #include "TTreeStream.h"
48 class AliTPCclusterFast: public TObject {
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);
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);
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);
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
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
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
89 TVectorD fStatY; //!stat Y
90 TVectorD fStatZ; //!stat Y
94 TMatrixD fDigits; // response matrix
95 static TF1* fPRF; // Pad response
96 static TF1* fTRF; // Time response function
97 ClassDef(AliTPCclusterFast,1) // container for
101 class AliTPCtrackFast: public TObject {
104 void Add(AliTPCtrackFast &track2);
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);
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);
115 Double_t CookdEdx(Int_t npoints, Double_t *amp, Double_t f0,Float_t f1, Int_t dEdxMode);
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
125 Bool_t fInit; // initialization flag
128 ClassDef(AliTPCtrackFast,2) // container for
133 ClassImp(AliTPCclusterFast)
134 ClassImp(AliTPCtrackFast)
140 TF1 *AliTPCclusterFast::fPRF=0;
141 TF1 *AliTPCclusterFast::fTRF=0;
144 AliTPCtrackFast::AliTPCtrackFast():
158 void AliTPCtrackFast::Add(AliTPCtrackFast &track2){
159 if (!track2.fInit) return;
165 void AliTPCtrackFast::MakeTrack(){
169 if (!fCl) fCl = new TClonesArray("AliTPCclusterFast",160);
171 // 0.) Init data structure
173 for (Int_t i=0;i<fN;i++){
174 AliTPCclusterFast * cluster = (AliTPCclusterFast*) fCl->UncheckedAt(i);
175 if (!cluster) cluster = new ((*fCl)[i]) AliTPCclusterFast;
179 // 1.) Create hits - with crosstalk diffusion
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;
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);
193 cluster->GenerElectrons(cluster, clusterm, clusterp);
196 // 2.) make digitization
198 for (Int_t i=0;i<fN;i++){
199 AliTPCclusterFast * cluster = (AliTPCclusterFast*) fCl->UncheckedAt(i);
205 Double_t AliTPCtrackFast::CookdEdxNtot(Double_t f0,Float_t f1){
207 // Double_t CookdEdxNtot(Double_t f0,Float_t f1); // dEdx_{hit} reconstructed meen number of electrons
210 for (Int_t i=0;i<fN;i++){
211 AliTPCclusterFast * cluster = ( AliTPCclusterFast *)((*fCl)[i]);
212 amp[i]=cluster->fNtot;
214 return CookdEdx(fN,amp,f0,f1,0);
217 Double_t AliTPCtrackFast::CookdEdxQtot(Double_t f0,Float_t f1){
219 // dEdx_{Q} reconstructed mean number of electronsxGain
222 for (Int_t i=0;i<fN;i++){
223 AliTPCclusterFast * cluster = ( AliTPCclusterFast *)((*fCl)[i]);
224 amp[i]=cluster->fQtot;
226 return CookdEdx(fN,amp,f0,f1,0);
230 Double_t AliTPCtrackFast::CookdEdxNtotThr(Double_t f0,Float_t f1, Double_t thr, Int_t dEdxMode){
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
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;
247 if (minAbove<0) minAbove=clQ;
248 if (minAbove>clQ) minAbove=clQ;
251 if (dEdxMode==-1) return Double_t(nBellow)/Double_t(fN);
253 for (Int_t i=0;i<fN;i++){
254 AliTPCclusterFast * cluster = ( AliTPCclusterFast *)((*fCl)[i]);
255 Double_t clQ= cluster->fNtot;
257 if (dEdxMode==0) amp[i]=clQ; // dEdxMode0 - not threshold - keep default
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
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
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
271 return CookdEdx(fN,amp,f0,f1, dEdxMode);
276 Double_t AliTPCtrackFast::CookdEdxQtotThr(Double_t f0,Float_t f1, Double_t thr, Int_t dEdxMode){
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
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;
296 if (minAbove<0) minAbove=clQ;
297 if (minAbove>clQ) minAbove=clQ;
300 if (dEdxMode==-1) return Double_t(nBellow)/Double_t(fN);
302 for (Int_t i=0;i<fN;i++){
303 AliTPCclusterFast * cluster = ( AliTPCclusterFast *)((*fCl)[i]);
304 Double_t clQ= cluster->fQtot;
306 if (dEdxMode==0) amp[i]=clQ; // dEdxMode0 - not threshold - keep default
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
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
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
320 return CookdEdx(fN,amp,f0,f1, dEdxMode);
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){
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}
335 for (Int_t i=0;i<fN;i++){
336 AliTPCclusterFast * cluster = ( AliTPCclusterFast *)((*fCl)[i]);
338 if (dEdxMode==0) camp = cluster->GetQtot(gain,0,noise);
340 camp = cluster->GetQtot(gain,thr,noise);
342 if (doCorr) corr = cluster->GetQtotCorr(0.5,0.5,gain,thr);
346 if (minAmp <0) minAmp=camp;
347 if (minAmp >camp) minAmp=camp;
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);
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){
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
365 for (Int_t i=0;i<fN;i++){
366 AliTPCclusterFast * cluster = ( AliTPCclusterFast *)((*fCl)[i]);
368 if (dEdxMode==0) camp = cluster->GetQmax(gain,0,noise);
370 camp = cluster->GetQmax(gain,thr,noise);
372 if (doCorr) corr = cluster->GetQmaxCorr(0.5,0.5);
376 if (minAmp <0) minAmp=camp;
377 if (minAmp >camp) minAmp=camp;
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);
386 Double_t AliTPCtrackFast::CookdEdx(Int_t npoints, Double_t *amp,Double_t f0,Float_t f1, Int_t dEdxMode){
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)
401 // 0. sorted the array of amplitudes
404 TMath::Sort(npoints,amp,index,kFALSE);
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;
411 for (Int_t i=0;i<npoints;i++) if (amp[index[i]]>0) above++;
413 for (Int_t i=0;i<npoints;i++){
415 if (dEdxMode==1 && amp[index[i]]==0) {
418 if (accepted<npoints*f0) continue;
419 if (accepted>npoints*f1) continue;
421 sum1+= amp[index[i]];
422 sum2+= amp[index[i]];
425 if (dEdxMode==-1) return 1-Double_t(above)/Double_t(npoints);
426 if (sum0<=0) return 0;
430 void AliTPCtrackFast::Simul(const char* fname, Int_t ntracks, Double_t diffFactor){
434 AliTPCtrackFast fast;
435 TTreeSRedirector *pcstream = new TTreeSRedirector(fname,"recreate");
436 for (Int_t itr=0; itr<ntracks; itr++){
438 fast.fMNprim=(10.+100*gRandom->Rndm());
439 if (gRandom->Rndm()>0.5) fast.fMNprim=1./(0.00001+gRandom->Rndm()*0.1);
441 fast.fDiff =0.01 +0.35*gRandom->Rndm();
442 // fast.fDiffLong = fast.fDiff*0.6/1.;
443 fast.fDiffLong = fast.fDiff*diffFactor/1.;
445 fast.fAngleY = 4.0*(gRandom->Rndm()-0.5);
446 fast.fAngleZ = 4.0*(gRandom->Rndm()-0.5);
449 if (itr%100==0) printf("%d\n",itr);
450 (*pcstream)<<"simulTrack"<<
460 AliTPCclusterFast::AliTPCclusterFast(){
463 fDigits.ResizeTo(5,7);
466 void AliTPCclusterFast::Init(){
468 // reset all counters
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
477 fPosY.ResizeTo(knMax);
478 fPosZ.ResizeTo(knMax);
479 fGain.ResizeTo(knMax);
480 fSec.ResizeTo(knMax);
483 for (Int_t i=0; i<knMax; i++){
493 AliTPCclusterFast::~AliTPCclusterFast(){
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){
500 fMNprim = mnprim; fDiff = diff; fDiffLong=diffL;
502 fAngleY=ky; fAngleZ=kz;
504 Double_t AliTPCclusterFast::GetNsec(){
506 // Generate number of secondary electrons
507 // copy of procedure implemented in geant
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);
519 void AliTPCclusterFast::GenerElectrons(AliTPCclusterFast *cl0, AliTPCclusterFast *clm, AliTPCclusterFast *clp){
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
534 // for (Int_t i=0;i<knMax;i++){
537 for (Int_t iprim=0; iprim<cl0->fNprim;iprim++){
538 Float_t dN = cl0->GetNsec();
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);
544 for (Int_t isec=0;isec<=dN;isec++){
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;
551 AliTPCclusterFast *cl=cl0;
552 if (r<-0.5 &&cl) cl=clm;
553 if (r>0.5 &&cl) cl=clp;
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;
564 // cl->sumY2Q+=gg*y*y;
566 // cl->sumZ2Q+=gg*z*z;
567 if (cl->fNtot>=knMax) continue;
569 if (cl0->fNtot>=knMax) break;
574 // fStatY[1]=sumYQ/sumQ;
575 // fStatY[2]=sumY2Q/sumQ-fStatY[1]*fStatY[1];
577 // fStatZ[1]=sumZQ/sumQ;
578 // fStatZ[2]=sumZ2Q/sumQ-fStatZ[1]*fStatZ[1];
582 void AliTPCclusterFast::Digitize(){
587 for (Int_t i=0; i<5;i++)
588 for (Int_t j=0; j<7;j++){
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]);
598 fDigits(2+di,3+dj)+=fac;
608 // void AliTPCclusterFast::Simul(const char* fname, Int_t npoints){
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;
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();
625 // if (icl%10000==0) printf("%d\n",icl);
626 // cstream<<"simul"<<
633 Double_t AliTPCclusterFast::GetQtot(Float_t gain, Float_t thr, Float_t noise, Bool_t brounding, Bool_t baddPedestal){
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;
650 Double_t AliTPCclusterFast::GetQmax(Float_t gain, Float_t thr, Float_t noise, Bool_t brounding, Bool_t baddPedestal){
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;
669 Double_t AliTPCclusterFast::GetQmaxCorr(Float_t rmsy0, Float_t rmsz0){
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
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);
684 Double_t AliTPCclusterFast::GetQtotCorr(Float_t rmsy0, Float_t rmsz0, Float_t gain, Float_t thr){
686 // Calculates the fraction of the charge over threshol to total charge
687 // The response function
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
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;
705 if (sumAll>0&&sumThr>0) corr=(sumThr)/sumAll;
706 if (sumThr==0) corr=GetQmaxCorr(0.5,0.5);
708 if (corr>0) qnorm=qtot/corr;
718 Double_t AliTPCclusterFast::GaussConvolution(Double_t x0, Double_t x1, Double_t k0, Double_t k1, Double_t s0, Double_t s1){
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}]]
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)
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());
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));
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);
748 Double_t AliTPCclusterFast::GaussExpConvolution(Double_t x0, Double_t s0,Double_t t1){
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;
765 Double_t AliTPCclusterFast::Gamma4(Double_t x, Double_t p0, Double_t p1){
767 // Gamma 4 Time response function of ALTRO
770 Double_t g1 = TMath::Exp(-4.*x/p1);
771 Double_t g2 = TMath::Power(x/p1,4);
777 Double_t AliTPCclusterFast::GaussGamma4(Double_t x, Double_t s0, Double_t p1){
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)
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);
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}]]
796 // No analytical solution
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}]]