]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/fastSimul/AliTPCclusterFast.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fixes for #86059: Install data when ALICE_ROOT does not point to source (Christian)
[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"
46a0eaea 24#include "THnSparse.h"
13f45dd8 25#include "TClonesArray.h"
0716cb63 26#include "TTreeStream.h"
27
28class AliTPCclusterFast: public TObject {
29public:
30 AliTPCclusterFast();
31 virtual ~AliTPCclusterFast();
32 void SetParam(Float_t mnprim, Float_t diff, Float_t y, Float_t z, Float_t ky, Float_t kz);
33 void GenerElectrons();
37045aeb 34 void Digitize();
13f45dd8 35 Double_t GetQtot(Float_t gain,Float_t thr, Float_t noise, Bool_t rounding=kTRUE, Bool_t addPedestal=kTRUE);
36 Double_t GetQmax(Float_t gain,Float_t thr, Float_t noise, Bool_t rounding=kTRUE, Bool_t addPedestal=kTRUE);
37 Double_t GetQmaxCorr(Float_t rmsy0, Float_t rmsz0);
38 Double_t GetQtotCorr(Float_t rmsy0, Float_t rmsz0, Float_t gain, Float_t thr);
39
37045aeb 40 Double_t GetNsec();
0716cb63 41 static void Simul(const char* simul, Int_t npoints);
13f45dd8 42 static Double_t GaussConvolution(Double_t x0, Double_t x1, Double_t k0, Double_t k1, Double_t s0, Double_t s1);
43 static Double_t GaussExpConvolution(Double_t x0, Double_t s0,Double_t t1);
44 static Double_t GaussGamma4(Double_t x, Double_t s0, Double_t p1);
45 static Double_t Gamma4(Double_t x, Double_t p0, Double_t p1);
0716cb63 46public:
47 Float_t fMNprim; // mean number of primary electrons
48 // //electrons part input
37045aeb 49 Int_t fNprim; // mean number of primary electrons
13f45dd8 50 Int_t fNtot; // total number of electrons
51 Float_t fQtot; // total charge - Gas gain flucuation taken into account
52 //
0716cb63 53 Float_t fDiff; // diffusion sigma
54 Float_t fY; // y position
55 Float_t fZ; // z postion
56 Float_t fAngleY; // y angle - tan(y)
57 Float_t fAngleZ; // z angle - tan z
37045aeb 58 //
59 //
0716cb63 60 // // electron part simul
13f45dd8 61 TVectorD fSec; //! number of secondary electrons
0716cb63 62 TVectorD fPosY; //! position y for each electron
63 TVectorD fPosZ; //! position z for each electron
64 TVectorD fGain; //! gg for each electron
37045aeb 65 //
13f45dd8 66 TVectorD fStatY; //!stat Y
67 TVectorD fStatZ; //!stat Y
0716cb63 68 //
37045aeb 69 // digitization part
70 //
71 TMatrixD fDigits; // response matrix
72 static TF1* fPRF; // Pad response
73 static TF1* fTRF; // Time response function
0716cb63 74 ClassDef(AliTPCclusterFast,1) // container for
75};
13f45dd8 76
77
78class AliTPCtrackFast: public TObject {
79public:
80 AliTPCtrackFast();
46a0eaea 81 void Add(AliTPCtrackFast &track2);
13f45dd8 82 void MakeTrack();
83 void UpdatedEdxHisto();
84 void MakeHisto();
85 static void Simul(const char* simul, Int_t ntracks);
86 Double_t CookdEdxNtot(Double_t f0,Float_t f1);
87 Double_t CookdEdxQtot(Double_t f0,Float_t f1);
46a0eaea 88 //
89 Double_t CookdEdxDtot(Double_t f0,Float_t f1, Float_t gain,Float_t thr, Float_t noise, Bool_t corr = kTRUE);
90 Double_t CookdEdxDmax(Double_t f0,Float_t f1,Float_t gain,Float_t thr, Float_t noise, Bool_t corr=kTRUE);
91 //
92 Double_t CookdEdx(Int_t npoints, Double_t *amp, Double_t f0,Float_t f1);
13f45dd8 93 //
94 Float_t fMNprim; // mean number of primary electrons
95 Float_t fAngleY; // y angle - tan(y)
96 Float_t fAngleZ; // z angle - tan z
97 Float_t fDiff; // diffusion
98 Int_t fN; // number of clusters
99 TClonesArray *fCl; // array of clusters
100 //
101 Bool_t fInit; // initialization flag
46a0eaea 102 THnSparse *fHistoNtot; // histograms of trunc mean Ntot
103 THnSparse *fHistoQtot; // histograms of trunc mean Qtot
104 THnSparse *fHistoQNtot; // histograms of trunc mean Qtot/Ntot
13f45dd8 105 //
46a0eaea 106 THnSparse *fHistoDtot; // histograms of trunc mean digit tot
107 THnSparse *fHistoDmax; // histograms of trunc mean digit max
108 THnSparse *fHistoDtotRaw; // histograms of trunc mean digit tot
109 THnSparse *fHistoDmaxRaw; // histograms of trunc mean digit max
13f45dd8 110
46a0eaea 111 //
112 //
13f45dd8 113 ClassDef(AliTPCtrackFast,2) // container for
114};
115
116
117
0716cb63 118ClassImp(AliTPCclusterFast)
13f45dd8 119ClassImp(AliTPCtrackFast)
120
121
122
123
0716cb63 124
37045aeb 125TF1 *AliTPCclusterFast::fPRF=0;
126TF1 *AliTPCclusterFast::fTRF=0;
0716cb63 127
128
13f45dd8 129AliTPCtrackFast::AliTPCtrackFast():
130 TObject(),
131 fMNprim(0),
132 fAngleY(0),
133 fAngleZ(0),
134 fN(0),
135 fCl(0),
46a0eaea 136 fInit(kFALSE),
137 fHistoNtot(0),
138 fHistoQtot(0),
139 fHistoQNtot(0),
140 fHistoDtot(0),
141 fHistoDmax(0),
142 fHistoDtotRaw(0),
143 fHistoDmaxRaw(0)
13f45dd8 144{
145 //
146 //
147 //
46a0eaea 148}
149
150void AliTPCtrackFast::Add(AliTPCtrackFast &track2){
151 if (!track2.fInit) return;
152
153 fHistoNtot->Add(track2.fHistoNtot); // histograms of trunc mean Ntot
154 fHistoQtot->Add(track2.fHistoQtot); // histograms of trunc mean Qtot
155 fHistoQNtot->Add(track2.fHistoQNtot); // histograms of trunc mean Qtot/Ntot
156 //
157 fHistoDtot->Add(track2.fHistoDtot); // histograms of trunc mean digit tot
158 fHistoDmax->Add(track2.fHistoDmax); // histograms of trunc mean digit max
159 fHistoDtotRaw->Add(track2.fHistoDtotRaw); // histograms of trunc mean digit tot
160 fHistoDmaxRaw->Add(track2.fHistoDmaxRaw); // histograms of trunc mean digit max
13f45dd8 161}
162
163void AliTPCtrackFast::MakeHisto(){
164 //
165 // make default histo
166 //
46a0eaea 167 // dEdx histogram THnSparse
168 // 0 - value
169 // 1 - fMNprim - number of generated primaries
170 // 2 - fNpoints
171 // 3 - fFraction
172 // 4 - fDiff
173 // 5 - fAngleY
174 // 6 - fAngleZ
175
176 Double_t xmin[7], xmax[7];
177 Int_t nbins[7];
13f45dd8 178 if (fInit) return;
46a0eaea 179 //
180 nbins[1] = 10; xmin[1]=10; xmax[1]=30; // fMNprim
181 nbins[2] = 8; xmin[2]=80; xmax[2]=160; // fNPoints
182 nbins[3] = 6; xmin[3]=0.45; xmax[3]=1.05; // trunc mean fraction
183
184 nbins[4] = 5; xmin[4]=0.0; xmax[4]=0.4; // fDiff
185 nbins[5] = 10; xmin[5]=0; xmax[5]=2; // fAngleY
186 nbins[6] = 10; xmin[6]=0; xmax[6]=2; // fAngleZ
187 //
188 nbins[0] =100; xmin[0]=2; xmax[0]=8;
189 fHistoNtot = new THnSparseF("dNdxall/dNdxprim","dNdxall/dNdxprim", 4, nbins, xmin,xmax);
190 nbins[0] =100; xmin[0]=2; xmax[0]=8;
191 fHistoQtot = new THnSparseF("dQdx/dNdxprim","dQdxall/dNdxprim", 4, nbins, xmin,xmax);
192 nbins[0] =100; xmin[0]=0.5; xmax[0]=1.5;
193 fHistoQNtot = new THnSparseF("dQdx/dNdxprim","dQdxprim/dNdxprim", 4, nbins, xmin,xmax);
194 //
195 nbins[0] =100; xmin[0]=0.05; xmax[0]=8;
196 fHistoDtot = new THnSparseF("dQtotdx/dNdxprim","dQtotdx/dNdx", 7, nbins, xmin,xmax);
197 fHistoDmax = new THnSparseF("dQmaxdx/dNdxprim","dQmaxdx/dNdx", 7, nbins, xmin,xmax);
198 fHistoDtotRaw = new THnSparseF("raw dQtotdx/dNdxprim","raw dQtotdx/dNdx", 7, nbins, xmin,xmax);
199 fHistoDmaxRaw = new THnSparseF("raw dQmaxdx/dNdxprim","raw dQmaxdx/dNdx", 7, nbins, xmin,xmax);
13f45dd8 200 fInit=kTRUE;
201}
202
203void AliTPCtrackFast::UpdatedEdxHisto(){
204 //
46a0eaea 205 //fill default histo
13f45dd8 206 //
207 if (!fInit) MakeHisto();
46a0eaea 208 Double_t x[7];
209 x[1] = fMNprim;
210 x[2] = fN;
211 //
212 x[4] = fDiff;
213 x[5] = TMath::Abs(fAngleY);
214 x[6] = TMath::Abs(fAngleZ);
215
216 for (Int_t i=0;i<7;i++){
217 Float_t frac = 0.5+Float_t(i)*0.1;
218 x[3] = frac;
219 Double_t cNtot = CookdEdxNtot(0.01,frac);
220 Double_t cQtot = CookdEdxQtot(0.01,frac);
221 // MC -using hits
222 x[0] = cNtot/fMNprim;
223 fHistoNtot->Fill(x);
224 x[0] = cQtot/fMNprim;
225 fHistoQtot->Fill(x);
226 x[0] = cQtot/cNtot;
227 fHistoQNtot->Fill(x);
228 // MC - using digits
229 Double_t dQtot = CookdEdxDtot(0.01,frac,1,2.5,1,kTRUE);
230 Double_t dQmax = CookdEdxDmax(0.01,frac,1,2.5,1,kTRUE);
231 Double_t dQrawtot = CookdEdxDtot(0.01,frac,1,2.5,1,kFALSE);
232 Double_t dQrawmax = CookdEdxDmax(0.01,frac,1,2.5,1,kFALSE);
233 x[0] = dQtot/fMNprim;
234 fHistoDtot->Fill(x);
235 x[0] = dQmax/fMNprim;
236 fHistoDmax->Fill(x);
237 x[0] = dQrawtot/fMNprim;
238 fHistoDtotRaw->Fill(x);
239 x[0] = dQrawmax/fMNprim;
240 fHistoDmaxRaw->Fill(x);
13f45dd8 241 }
242}
243
244void AliTPCtrackFast::MakeTrack(){
245 //
246 //
247 //
46a0eaea 248 if (!fCl) fCl = new TClonesArray("AliTPCclusterFast",160);
13f45dd8 249 for (Int_t i=0;i<fN;i++){
46a0eaea 250 Double_t tY = i*fAngleY;
251 Double_t tZ = i*fAngleZ;
252 AliTPCclusterFast * cluster = (AliTPCclusterFast*) fCl->UncheckedAt(i);
253 if (!cluster) cluster = new ((*fCl)[i]) AliTPCclusterFast;
254 //
255 Double_t posY = tY-TMath::Nint(tY);
256 Double_t posZ = tZ-TMath::Nint(tZ);
257 cluster->SetParam(fMNprim,fDiff,posY,posZ,fAngleY,fAngleZ);
258 //
13f45dd8 259 cluster->GenerElectrons();
260 cluster->Digitize();
261 }
262 UpdatedEdxHisto();
263}
264
265Double_t AliTPCtrackFast::CookdEdxNtot(Double_t f0,Float_t f1){
266 //
267 Double_t amp[160];
13f45dd8 268 for (Int_t i=0;i<fN;i++){
269 AliTPCclusterFast * cluster = ( AliTPCclusterFast *)((*fCl)[i]);
270 amp[i]=cluster->fNtot;
271 }
46a0eaea 272 return CookdEdx(fN,amp,f0,f1);
13f45dd8 273}
274
275Double_t AliTPCtrackFast::CookdEdxQtot(Double_t f0,Float_t f1){
276 //
277 Double_t amp[160];
13f45dd8 278 for (Int_t i=0;i<fN;i++){
279 AliTPCclusterFast * cluster = ( AliTPCclusterFast *)((*fCl)[i]);
280 amp[i]=cluster->fQtot;
281 }
46a0eaea 282 return CookdEdx(fN,amp,f0,f1);
13f45dd8 283}
284
46a0eaea 285Double_t AliTPCtrackFast::CookdEdxDtot(Double_t f0,Float_t f1, Float_t gain,Float_t thr, Float_t noise, Bool_t doCorr){
286 //
287 //
288 //
289 Double_t amp[160];
290 Int_t over=0;
291 for (Int_t i=0;i<fN;i++){
292 AliTPCclusterFast * cluster = ( AliTPCclusterFast *)((*fCl)[i]);
293 Float_t camp = cluster->GetQtot(gain,thr,noise);
294 if (camp==0) continue;
295 Float_t corr = 1;
296 if (doCorr) corr = cluster->GetQtotCorr(0.5,0.5,gain,thr);
297 amp[over]=camp/corr;
298 over++;
299 }
300 return CookdEdx(over,amp,f0,f1);
301
302}
303
304Double_t AliTPCtrackFast::CookdEdxDmax(Double_t f0,Float_t f1, Float_t gain,Float_t thr, Float_t noise, Bool_t doCorr){
305 //
306 //
307 //
308 Double_t amp[160];
309 Int_t over=0;
310 for (Int_t i=0;i<fN;i++){
311 AliTPCclusterFast * cluster = ( AliTPCclusterFast *)((*fCl)[i]);
312 Float_t camp = cluster->GetQmax(gain,thr,noise);
313 if (camp==0) continue;
314 Float_t corr = 1;
315 if (doCorr) corr = cluster->GetQmaxCorr(0.5,0.5);
316 amp[over]=camp/corr;
317 over++;
318 }
319 return CookdEdx(over,amp,f0,f1);
320
321}
322
323
324Double_t AliTPCtrackFast::CookdEdx(Int_t npoints, Double_t *amp,Double_t f0,Float_t f1){
13f45dd8 325 //
326 //
327 //
328 Int_t index[160];
46a0eaea 329 TMath::Sort(npoints,amp,index,kFALSE);
13f45dd8 330 Float_t sum0=0, sum1=0,sum2=0;
46a0eaea 331 for (Int_t i=0;i<npoints;i++){
332 if (i<npoints*f0) continue;
333 if (i>npoints*f1) continue;
13f45dd8 334 sum0++;
335 sum1+= amp[index[i]];
336 sum2+= amp[index[i]];
337 }
338 if (sum0<=0) return 0;
339 return sum1/sum0;
340}
341
342void AliTPCtrackFast::Simul(const char* fname, Int_t ntracks){
343 //
344 //
345 //
346 AliTPCtrackFast fast;
347 TTreeSRedirector cstream(fname);
348 for (Int_t itr=0; itr<ntracks; itr++){
349 //
c9e8d4ec 350 fast.fMNprim=(5+50*gRandom->Rndm());
13f45dd8 351 fast.fDiff =0.01 +0.35*gRandom->Rndm();
352 //
353 fast.fAngleY = 4.0*(gRandom->Rndm()-0.5);
354 fast.fAngleZ = 4.0*(gRandom->Rndm()-0.5);
46a0eaea 355 fast.fN = TMath::Nint(80.+gRandom->Rndm()*80.);
13f45dd8 356 fast.MakeTrack();
357 if (itr%100==0) printf("%d\n",itr);
358 cstream<<"simulTrack"<<
359 "tr.="<<&fast<<
360 "\n";
361 }
362 fast.Write("track");
363}
364
365
366
0716cb63 367AliTPCclusterFast::AliTPCclusterFast(){
37045aeb 368 //
369 //
370 fDigits.ResizeTo(5,7);
0716cb63 371}
37045aeb 372
0716cb63 373AliTPCclusterFast::~AliTPCclusterFast(){
374}
375
376
377void AliTPCclusterFast::SetParam(Float_t mnprim, Float_t diff, Float_t y, Float_t z, Float_t ky, Float_t kz){
378 //
379 //
380 fMNprim = mnprim; fDiff = diff;
381 fY=y; fZ=z;
382 fAngleY=ky; fAngleZ=kz;
383}
37045aeb 384Double_t AliTPCclusterFast::GetNsec(){
0716cb63 385 //
37045aeb 386 // Generate number of secondary electrons
387 // copy of procedure implemented in geant
0716cb63 388 //
0716cb63 389 const Double_t FPOT=20.77E-9, EEND=10E-6, EEXPO=2.2, EEND1=1E-6;
390 const Double_t XEXPO=-EEXPO+1, YEXPO=1/XEXPO;
391 const Double_t W=20.77E-9;
37045aeb 392 Float_t RAN = gRandom->Rndm();
46a0eaea 393 //Double_t edep = TMath::Power((TMath::Power(FPOT,XEXPO)*(1-RAN)+TMath::Power(EEND,XEXPO)*RAN),YEXPO);
394 //edep = TMath::Min(edep, EEND);
395 //return TMath::Nint(edep/W);
37045aeb 396 return TMath::Nint(TMath::Power((TMath::Power(FPOT,XEXPO)*(1-RAN)+TMath::Power(EEND,XEXPO)*RAN),YEXPO)/W);
397}
398
399void AliTPCclusterFast::GenerElectrons(){
0716cb63 400 //
37045aeb 401 //
402 //
403 //
404 const Int_t knMax=1000;
0716cb63 405 if (fPosY.GetNrows()<knMax){
406 fPosY.ResizeTo(knMax);
407 fPosZ.ResizeTo(knMax);
408 fGain.ResizeTo(knMax);
37045aeb 409 fSec.ResizeTo(knMax);
0716cb63 410 fStatY.ResizeTo(3);
411 fStatZ.ResizeTo(3);
412 }
37045aeb 413 fNprim = gRandom->Poisson(fMNprim); //number of primary electrons
13f45dd8 414 fNtot=0; //total number of electrons
415 fQtot=0; //total number of electrons after gain multiplification
0716cb63 416 //
417 Double_t sumQ=0;
418 Double_t sumYQ=0;
419 Double_t sumZQ=0;
420 Double_t sumY2Q=0;
421 Double_t sumZ2Q=0;
37045aeb 422 for (Int_t i=0;i<knMax;i++){
423 fSec[i]=0;
424 }
0716cb63 425 for (Int_t iprim=0; iprim<fNprim;iprim++){
37045aeb 426 Float_t dN = GetNsec();
427 fSec[iprim]=dN;
0716cb63 428 Double_t yc = fY+(gRandom->Rndm()-0.5)*fAngleY;
429 Double_t zc = fZ+(gRandom->Rndm()-0.5)*fAngleZ;
37045aeb 430 for (Int_t isec=0;isec<=dN;isec++){
0716cb63 431 //
432 //
433 Double_t y = gRandom->Gaus(0,fDiff)+yc;
434 Double_t z = gRandom->Gaus(0,fDiff)+zc;
435 Double_t gg = -TMath::Log(gRandom->Rndm());
436 fPosY[fNtot]=y;
437 fPosZ[fNtot]=z;
438 fGain[fNtot]=gg;
13f45dd8 439 fQtot+=gg;
0716cb63 440 fNtot++;
441 sumQ+=gg;
442 sumYQ+=gg*y;
443 sumY2Q+=gg*y*y;
444 sumZQ+=gg*z;
445 sumZ2Q+=gg*z*z;
446 if (fNtot>=knMax) break;
447 }
448 if (fNtot>=knMax) break;
449 }
450 if (sumQ>0){
451 fStatY[0]=sumQ;
452 fStatY[1]=sumYQ/sumQ;
453 fStatY[2]=sumY2Q/sumQ-fStatY[1]*fStatY[1];
454 fStatZ[0]=sumQ;
455 fStatZ[1]=sumZQ/sumQ;
456 fStatZ[2]=sumZ2Q/sumQ-fStatZ[1]*fStatZ[1];
457 }
458}
459
37045aeb 460void AliTPCclusterFast::Digitize(){
461 //
462 //
463 //
464 // 1. Clear digits
465 for (Int_t i=0; i<5;i++)
466 for (Int_t j=0; j<7;j++){
467 fDigits(i,j)=0;
468 }
469 //
470 // Fill digits
471 for (Int_t iel = 0; iel<fNtot; iel++){
472 for (Int_t di=-2; di<=2;di++)
13f45dd8 473 for (Int_t dj=-3; dj<=3;dj++){
37045aeb 474 Float_t fac = fPRF->Eval(di-fPosY[iel])*fTRF->Eval(dj-fPosZ[iel]);
475 fac*=fGain[iel];
476 fDigits(2+di,3+dj)+=fac;
477 }
478 }
13f45dd8 479
37045aeb 480}
481
482
0716cb63 483
484void AliTPCclusterFast::Simul(const char* fname, Int_t npoints){
13f45dd8 485 //
486 // Calc rms
487 //
0716cb63 488 AliTPCclusterFast fast;
489 TTreeSRedirector cstream(fname);
490 for (Int_t icl=0; icl<npoints; icl++){
13f45dd8 491 Float_t nprim=(10+20*gRandom->Rndm());
a314e83c 492 Float_t diff =0.01 +0.35*gRandom->Rndm();
0716cb63 493 Float_t posY = gRandom->Rndm()-0.5;
494 Float_t posZ = gRandom->Rndm()-0.5;
13f45dd8 495 //
496 Float_t ky = 4.0*(gRandom->Rndm()-0.5);
497 Float_t kz = 4.0*(gRandom->Rndm()-0.5);
0716cb63 498 fast.SetParam(nprim,diff,posY,posZ,ky,kz);
499 fast.GenerElectrons();
37045aeb 500 fast.Digitize();
13f45dd8 501 if (icl%10000==0) printf("%d\n",icl);
0716cb63 502 cstream<<"simul"<<
503 "s.="<<&fast<<
504 "\n";
505 }
506}
507
508
13f45dd8 509Double_t AliTPCclusterFast::GetQtot(Float_t gain, Float_t thr, Float_t noise, Bool_t brounding, Bool_t baddPedestal){
a314e83c 510 //
511 //
512 //
513 Float_t sum =0;
13f45dd8 514 for (Int_t ip=0;ip<5;ip++){
515 Float_t pedestal=gRandom->Rndm()-0.5; //pedestal offset different for each pad
516 for (Int_t it=0;it<7;it++){
517 Float_t amp = gain*fDigits(ip,it)+gRandom->Gaus()*noise;
518 if (baddPedestal) amp+=pedestal;
519 if (brounding) amp=TMath::Nint(amp);
a314e83c 520 if (amp>thr) sum+=amp;
521 }
13f45dd8 522 }
a314e83c 523 return sum;
524}
525
13f45dd8 526Double_t AliTPCclusterFast::GetQmax(Float_t gain, Float_t thr, Float_t noise, Bool_t brounding, Bool_t baddPedestal){
a314e83c 527 //
528 //
529 //
530 Float_t max =0;
13f45dd8 531 for (Int_t ip=0;ip<5;ip++){
532 Float_t pedestal=gRandom->Rndm()-0.5; //pedestal offset different for each pad
533 for (Int_t it=0;it<7;it++){
534 Float_t amp = gain*fDigits(ip,it)+gRandom->Gaus()*noise;
535 if (baddPedestal) amp+=pedestal;
536 if (brounding) amp=TMath::Nint(amp);
537 if (amp>max && amp>thr) max=amp;
a314e83c 538 }
13f45dd8 539 }
a314e83c 540 return max;
541}
542
543
a314e83c 544
13f45dd8 545Double_t AliTPCclusterFast::GetQmaxCorr(Float_t rmsy0, Float_t rmsz0){
546 //
547 // Gaus distribution convolueted with rectangular
548 // Gaus width sy and sz is determined by RF width and diffusion
549 // Integral of Q is equal 1
550 // Q max is calculated at position fY,fX
551 //
552 //
553 //
554 Double_t sy = TMath::Sqrt(rmsy0*rmsy0+fDiff*fDiff);
555 Double_t sz = TMath::Sqrt(rmsz0*rmsz0+fDiff*fDiff);
556 return GaussConvolution(fY,fZ, fAngleY,fAngleZ,sy,sz);
557}
558
559
560Double_t AliTPCclusterFast::GetQtotCorr(Float_t rmsy0, Float_t rmsz0, Float_t gain, Float_t thr){
561 //
562 // Calculates the fraction of the charge over threshol to total charge
563 // The response function
564 //
565 Double_t sy = TMath::Sqrt(rmsy0*rmsy0+fDiff*fDiff);
566 Double_t sz = TMath::Sqrt(rmsz0*rmsz0+fDiff*fDiff);
567 Double_t sumAll=0,sumThr=0;
568 Double_t qtot = GetQtot(gain,thr,0); // sum of signal over threshold
569 //
570 Double_t corr =1;
571 Double_t qnorm=qtot;
572 for (Int_t iter=0;iter<2;iter++){
573 for (Int_t iy=-2;iy<=2;iy++)
574 for (Int_t iz=-2;iz<=2;iz++){
575 Double_t val = GaussConvolution(fY-iy,fZ-iz, fAngleY,fAngleZ,sy,sz);
576 Double_t qlocal =TMath::Nint(qnorm*val);
577 if (qlocal>thr) sumThr+=qlocal;
578 sumAll+=qlocal;
579 }
580 if (sumAll>0&&sumThr>0) corr=(sumThr)/sumAll;
581 //corr = sumThr;
582 if (corr>0) qnorm=qtot/corr;
583
584 }
585 return corr;
586}
a314e83c 587
588
0716cb63 589
13f45dd8 590
591
592Double_t AliTPCclusterFast::GaussConvolution(Double_t x0, Double_t x1, Double_t k0, Double_t k1, Double_t s0, Double_t s1){
593 //
594 // 2 D gaus convoluted with angular effect
595 // See in mathematica:
596 //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}]]
597 //
598 //TF1 f1("f1","AliTPCclusterFast::GaussConvolution(x,0,1,0,0.1,0.1)",-2,2)
599 //TF2 f2("f2","AliTPCclusterFast::GaussConvolution(x,y,1,1,0.1,0.1)",-2,2,-2,2)
600 //
601 const Float_t kEpsilon = 0.0001;
602 if ((TMath::Abs(k0)+TMath::Abs(k1))<kEpsilon*(s0+s1)){
603 // small angular effect
604 Double_t val = (TMath::Gaus(x0,0,s0)*TMath::Gaus(x1,0,s1))/(s0*s1*2.*TMath::Pi());
605 return val;
606 }
607
608 Double_t sigma2 = k1*k1*s0*s0+k0*k0*s1*s1;
609 Double_t exp0 = TMath::Exp(-(k1*x0-k0*x1)*(k1*x0-k0*x1)/(2*sigma2));
610 //
611 Double_t sigmaErf = 2*s0*s1*TMath::Sqrt(2*sigma2);
612 Double_t erf0 = TMath::Erf( (k0*s1*s1*(k0-2*x0)+k1*s0*s0*(k1-2*x1))/sigmaErf);
613 Double_t erf1 = TMath::Erf( (k0*s1*s1*(k0+2*x0)+k1*s0*s0*(k1+2*x1))/sigmaErf);
614 Double_t norm = 1./TMath::Sqrt(sigma2);
615 norm/=2.*TMath::Sqrt(2.*TMath::Pi());
616 Double_t val = norm*exp0*(erf0+erf1);
617 return val;
618
619}
620
621
622Double_t AliTPCclusterFast::GaussExpConvolution(Double_t x0, Double_t s0,Double_t t1){
623 //
624 // 2 D gaus convoluted with exponential
625 // Integral nomalized to 1
626 // See in mathematica:
627 //Simplify[Integrate[Exp[-(x0-x1)*(x0-x1)/(2*s0*s0)]*Exp[-x1*t1],{x1,0,Infinity}]]
628 // TF1 fgexp("fgexp","AliTPCclusterFast::GaussExpConvolution(x,0.5,1)",-2,2)
629 Double_t exp1 = (s0*s0*t1-2*x0)*t1/2.;
630 exp1 = TMath::Exp(exp1);
631 Double_t erf = 1+TMath::Erf((-s0*s0*t1+x0)/(s0*TMath::Sqrt(2.)));
632 Double_t val = exp1*erf;
633 val *=t1/(2.);
634 return val;
635
636}
637
638
639Double_t AliTPCclusterFast::Gamma4(Double_t x, Double_t p0, Double_t p1){
640 //
641 // Gamma 4 Time response function of ALTRO
642 //
643 if (x<0) return 0;
644 Double_t g1 = TMath::Exp(-4.*x/p1);
645 Double_t g2 = TMath::Power(x/p1,4);
646 return p0*g1*g2;
647}
648
649
650
651Double_t AliTPCclusterFast::GaussGamma4(Double_t x, Double_t s0, Double_t p1){
652 //
653 // Gamma 4 Time response function of ALTRO convoluted with Gauss
654 // Simplify[Integrate[Exp[-(x0-x1)*(x0-x1)/(2*s0*s0)]*Exp[-4*x1/p1]*(x/p1)^4/s0,{x1,0,Infinity}]]
655 //TF1 fgg4("fgg4","AliTPCclusterFast::GaussGamma4(x,0.5,0.5)",-2,2)
656
657 Double_t exp1 = (8*s0*s0-4.*p1*x)/(p1*p1);
658 exp1 = TMath::Exp(exp1);
659 Double_t erf1 = 1+TMath::Erf((-4*s0/p1+x/s0)/TMath::Sqrt(2));
46a0eaea 660 // Double_t xp14 = TMath::Power(TMath::Abs((x/p1)),4);
13f45dd8 661 return exp1*erf1;
662
663
664}
665
666// Analytical sollution only in 1D - too long expression
667// Simplify[Integrate[Exp[-(x0-(x1-k*x2))*(x0-(x1-k*x2))/(2*s0*s0)]*Exp[-(x1*t1-k*x2)],{x2,-1,1}]]
668//
669//
670// No analytical solution
671//
672//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}]]