git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
ATO-58 digitization and hit creation in anohter loop (to take into account properly...
[u/mrichter/AliRoot.git] / TPC / fastSimul / AliTPCclusterFast.cxx
CommitLineData
fa9e5d19 1
0716cb63 2/*
99957137 3
4
5
0716cb63 6 gSystem->Load("libSTAT.so");
37045aeb 7 .x ~/NimStyle.C
0716cb63 8
9 .L $ALICE_ROOT/TPC/fastSimul/AliTPCclusterFast.cxx+
37045aeb 10 //
13f45dd8 11 AliTPCclusterFast::fPRF = new TF1("fprf","gausn",-5,5);
12 AliTPCclusterFast::fTRF = new TF1("ftrf","gausn",-5,5);
37045aeb 13 AliTPCclusterFast::fPRF->SetParameters(1,0,0.5);
14 AliTPCclusterFast::fTRF->SetParameters(1,0,0.5);
0716cb63 15
13f45dd8 16 //
46a0eaea 17 AliTPCtrackFast::Simul("trackerSimul.root",100);
18// AliTPCclusterFast::Simul("cluterSimul.root",20000);
0716cb63 19*/
20
994b9770 21
22
23
24/*
25 Modifications to add:
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.
31 hwMode
32
33 */
34
35
36
0716cb63 37#include "TObject.h"
37045aeb 38#include "TF1.h"
0716cb63 39#include "TMath.h"
40#include "TRandom.h"
41#include "TVectorD.h"
37045aeb 42#include "TMatrixD.h"
13f45dd8 43#include "TH1.h"
13f45dd8 44#include "TClonesArray.h"
0716cb63 45#include "TTreeStream.h"
e05accd3 46#include "TGrid.h"
0716cb63 47
48class AliTPCclusterFast: public TObject {
49public:
50 AliTPCclusterFast();
f0266161 51 void Init();
52
0716cb63 53 virtual ~AliTPCclusterFast();
f0266161 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);
37045aeb 56 void Digitize();
13f45dd8 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);
61
37045aeb 62 Double_t GetNsec();
f0266161 63 //static void Simul(const char* simul, Int_t npoints);
13f45dd8 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);
0716cb63 68public:
69 Float_t fMNprim; // mean number of primary electrons
70 // //electrons part input
37045aeb 71 Int_t fNprim; // mean number of primary electrons
13f45dd8 72 Int_t fNtot; // total number of electrons
73 Float_t fQtot; // total charge - Gas gain flucuation taken into account
74 //
0716cb63 75 Float_t fDiff; // diffusion sigma
f0266161 76 Float_t fDiffLong; // diffusion sigma longitudinal direction
0716cb63 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
37045aeb 81 //
82 //
0716cb63 83 // // electron part simul
13f45dd8 84 TVectorD fSec; //! number of secondary electrons
0716cb63 85 TVectorD fPosY; //! position y for each electron
86 TVectorD fPosZ; //! position z for each electron
87 TVectorD fGain; //! gg for each electron
37045aeb 88 //
13f45dd8 89 TVectorD fStatY; //!stat Y
90 TVectorD fStatZ; //!stat Y
0716cb63 91 //
37045aeb 92 // digitization part
93 //
94 TMatrixD fDigits; // response matrix
95 static TF1* fPRF; // Pad response
96 static TF1* fTRF; // Time response function
0716cb63 97 ClassDef(AliTPCclusterFast,1) // container for
98};
13f45dd8 99
100
101class AliTPCtrackFast: public TObject {
102public:
103 AliTPCtrackFast();
46a0eaea 104 void Add(AliTPCtrackFast &track2);
13f45dd8 105 void MakeTrack();
52560b2e 106 static void Simul(const char* simul, Int_t ntracks, Double_t diff);
13f45dd8 107 Double_t CookdEdxNtot(Double_t f0,Float_t f1);
108 Double_t CookdEdxQtot(Double_t f0,Float_t f1);
994b9770 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);
46a0eaea 111 //
994b9770 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);
46a0eaea 114 //
994b9770 115 Double_t CookdEdx(Int_t npoints, Double_t *amp, Double_t f0,Float_t f1, Int_t dEdxMode);
13f45dd8 116 //
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
f0266161 121 Float_t fDiffLong; // diffusion sigma longitudinal direction
13f45dd8 122 Int_t fN; // number of clusters
123 TClonesArray *fCl; // array of clusters
124 //
125 Bool_t fInit; // initialization flag
46a0eaea 126 //
127 //
13f45dd8 128 ClassDef(AliTPCtrackFast,2) // container for
129};
130
131
132
0716cb63 133ClassImp(AliTPCclusterFast)
13f45dd8 134ClassImp(AliTPCtrackFast)
135
136
137
138
0716cb63 139
37045aeb 140TF1 *AliTPCclusterFast::fPRF=0;
141TF1 *AliTPCclusterFast::fTRF=0;
0716cb63 142
143
13f45dd8 144AliTPCtrackFast::AliTPCtrackFast():
145 TObject(),
146 fMNprim(0),
147 fAngleY(0),
148 fAngleZ(0),
149 fN(0),
150 fCl(0),
44a0fef0 151 fInit(kFALSE)
13f45dd8 152{
153 //
154 //
155 //
46a0eaea 156}
157
158void AliTPCtrackFast::Add(AliTPCtrackFast &track2){
159 if (!track2.fInit) return;
160
13f45dd8 161}
162
13f45dd8 163
13f45dd8 164
165void AliTPCtrackFast::MakeTrack(){
166 //
167 //
168 //
46a0eaea 169 if (!fCl) fCl = new TClonesArray("AliTPCclusterFast",160);
e05accd3 170 //
171 // 0.) Init data structure
172 //
f0266161 173 for (Int_t i=0;i<fN;i++){
174 AliTPCclusterFast * cluster = (AliTPCclusterFast*) fCl->UncheckedAt(i);
175 if (!cluster) cluster = new ((*fCl)[i]) AliTPCclusterFast;
176 cluster->Init();
177 }
e05accd3 178 //
179 // 1.) Create hits - with crosstalk diffusion
180 //
13f45dd8 181 for (Int_t i=0;i<fN;i++){
46a0eaea 182 Double_t tY = i*fAngleY;
183 Double_t tZ = i*fAngleZ;
184 AliTPCclusterFast * cluster = (AliTPCclusterFast*) fCl->UncheckedAt(i);
f0266161 185 AliTPCclusterFast * clusterm = (AliTPCclusterFast*) fCl->UncheckedAt(TMath::Max(i-1,0));
186 AliTPCclusterFast * clusterp = (AliTPCclusterFast*) fCl->UncheckedAt(TMath::Min(i+1,159));
46a0eaea 187 if (!cluster) cluster = new ((*fCl)[i]) AliTPCclusterFast;
188 //
189 Double_t posY = tY-TMath::Nint(tY);
190 Double_t posZ = tZ-TMath::Nint(tZ);
f0266161 191 cluster->SetParam(fMNprim,fDiff, fDiffLong, posY,posZ,fAngleY,fAngleZ);
46a0eaea 192 //
f0266161 193 cluster->GenerElectrons(cluster, clusterm, clusterp);
e05accd3 194 }
195 //
196 // 2.) make digitization
197 //
198 for (Int_t i=0;i<fN;i++){
199 AliTPCclusterFast * cluster = (AliTPCclusterFast*) fCl->UncheckedAt(i);
13f45dd8 200 cluster->Digitize();
201 }
e05accd3 202
13f45dd8 203}
204
205Double_t AliTPCtrackFast::CookdEdxNtot(Double_t f0,Float_t f1){
99957137 206 //
207 // Double_t CookdEdxNtot(Double_t f0,Float_t f1); // dEdx_{hit} reconstructed meen number of electrons
13f45dd8 208 //
209 Double_t amp[160];
13f45dd8 210 for (Int_t i=0;i<fN;i++){
211 AliTPCclusterFast * cluster = ( AliTPCclusterFast *)((*fCl)[i]);
212 amp[i]=cluster->fNtot;
213 }
99957137 214 return CookdEdx(fN,amp,f0,f1,0);
13f45dd8 215}
216
217Double_t AliTPCtrackFast::CookdEdxQtot(Double_t f0,Float_t f1){
99957137 218 //
219 // dEdx_{Q} reconstructed mean number of electronsxGain
13f45dd8 220 //
221 Double_t amp[160];
13f45dd8 222 for (Int_t i=0;i<fN;i++){
223 AliTPCclusterFast * cluster = ( AliTPCclusterFast *)((*fCl)[i]);
224 amp[i]=cluster->fQtot;
225 }
99957137 226 return CookdEdx(fN,amp,f0,f1,0);
227}
228
229
994b9770 230Double_t AliTPCtrackFast::CookdEdxNtotThr(Double_t f0,Float_t f1, Double_t thr, Int_t dEdxMode){
99957137 231 //
232 // dEdx_{hit} reconstructed mean number of electrons
233 // thr = threshold in terms of the number of electrons
994b9770 234 // dEdxMode = algorithm to deal with trhesold values replacing
99957137 235 //
236 Double_t amp[160];
237 Int_t nBellow=0;
238 //
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;
243 if (clQ<thr) {
244 nBellow++;
245 continue;
246 }
247 if (minAbove<0) minAbove=clQ;
248 if (minAbove>clQ) minAbove=clQ;
249 }
250 //
994b9770 251 if (dEdxMode==-1) return Double_t(nBellow)/Double_t(fN);
99957137 252
253 for (Int_t i=0;i<fN;i++){
254 AliTPCclusterFast * cluster = ( AliTPCclusterFast *)((*fCl)[i]);
255 Double_t clQ= cluster->fNtot;
256 //
994b9770 257 if (dEdxMode==0) amp[i]=clQ; // dEdxMode0 - not threshold - keep default
99957137 258 //
259 //
994b9770 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
99957137 262 //
263 //
994b9770 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
99957137 266 //
267 //
994b9770 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
99957137 270 }
994b9770 271 return CookdEdx(fN,amp,f0,f1, dEdxMode);
13f45dd8 272}
273
99957137 274
275
994b9770 276Double_t AliTPCtrackFast::CookdEdxQtotThr(Double_t f0,Float_t f1, Double_t thr, Int_t dEdxMode){
fa9e5d19 277 //
278 //
99957137 279 // dEdx_{Q} reconstructed mean number of electrons xgain
280 // thr = threshold in terms of the number of electrons
994b9770 281 // dEdxMode = algorithm to deal with trhesold values replacing
99957137 282 //
283
fa9e5d19 284 //
285 Double_t amp[160];
fa9e5d19 286 Int_t nBellow=0;
287 //
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;
292 if (clQ<thr) {
293 nBellow++;
294 continue;
295 }
296 if (minAbove<0) minAbove=clQ;
297 if (minAbove>clQ) minAbove=clQ;
298 }
99957137 299 //
994b9770 300 if (dEdxMode==-1) return Double_t(nBellow)/Double_t(fN);
fa9e5d19 301
302 for (Int_t i=0;i<fN;i++){
303 AliTPCclusterFast * cluster = ( AliTPCclusterFast *)((*fCl)[i]);
304 Double_t clQ= cluster->fQtot;
99957137 305 //
994b9770 306 if (dEdxMode==0) amp[i]=clQ; // dEdxMode0 - not threshold - keep default
99957137 307 //
308 //
994b9770 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
99957137 311 //
312 //
994b9770 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
99957137 315 //
316 //
994b9770 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
fa9e5d19 319 }
994b9770 320 return CookdEdx(fN,amp,f0,f1, dEdxMode);
fa9e5d19 321}
322
323
324
325
326
994b9770 327Double_t AliTPCtrackFast::CookdEdxDtot(Double_t f0,Float_t f1, Float_t gain,Float_t thr, Float_t noise, Bool_t doCorr, Int_t dEdxMode){
46a0eaea 328 //
99957137 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}
46a0eaea 331 //
332 Double_t amp[160];
99957137 333 Double_t minAmp=-1;
334 //
46a0eaea 335 for (Int_t i=0;i<fN;i++){
336 AliTPCclusterFast * cluster = ( AliTPCclusterFast *)((*fCl)[i]);
99957137 337 Float_t camp = 0;
994b9770 338 if (dEdxMode==0) camp = cluster->GetQtot(gain,0,noise);
99957137 339 else
340 camp = cluster->GetQtot(gain,thr,noise);
46a0eaea 341 Float_t corr = 1;
342 if (doCorr) corr = cluster->GetQtotCorr(0.5,0.5,gain,thr);
99957137 343 camp/=corr;
344 amp[i]=camp;
345 if (camp>0){
346 if (minAmp <0) minAmp=camp;
347 if (minAmp >camp) minAmp=camp;
348 }
46a0eaea 349 }
994b9770 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);
46a0eaea 353}
354
99957137 355
356
994b9770 357Double_t AliTPCtrackFast::CookdEdxDmax(Double_t f0,Float_t f1, Float_t gain,Float_t thr, Float_t noise, Bool_t doCorr, Int_t dEdxMode){
46a0eaea 358 //
99957137 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
46a0eaea 361 //
362 Double_t amp[160];
99957137 363 Double_t minAmp=-1;
364 //
46a0eaea 365 for (Int_t i=0;i<fN;i++){
366 AliTPCclusterFast * cluster = ( AliTPCclusterFast *)((*fCl)[i]);
99957137 367 Float_t camp = 0;
994b9770 368 if (dEdxMode==0) camp = cluster->GetQmax(gain,0,noise);
99957137 369 else
370 camp = cluster->GetQmax(gain,thr,noise);
46a0eaea 371 Float_t corr = 1;
372 if (doCorr) corr = cluster->GetQmaxCorr(0.5,0.5);
99957137 373 camp/=corr;
374 amp[i]=camp;
375 if (camp>0){
376 if (minAmp <0) minAmp=camp;
377 if (minAmp >camp) minAmp=camp;
378 }
46a0eaea 379 }
994b9770 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);
46a0eaea 383}
384
385
994b9770 386Double_t AliTPCtrackFast::CookdEdx(Int_t npoints, Double_t *amp,Double_t f0,Float_t f1, Int_t dEdxMode){
99957137 387 //
388 // Calculate truncated mean
389 // npoints - number of points in array
390 // amp - array with points
391 // f0-f1 - truncation range
994b9770 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)
13f45dd8 398 //
99957137 399
13f45dd8 400 //
99957137 401 // 0. sorted the array of amplitudes
13f45dd8 402 //
403 Int_t index[160];
46a0eaea 404 TMath::Sort(npoints,amp,index,kFALSE);
99957137 405 //
406 // 1.) Calculate truncated mean from the selected range of the array (ranking statistic )
994b9770 407 // dependening on the dEdxMode 0 amplitude can be skipped
13f45dd8 408 Float_t sum0=0, sum1=0,sum2=0;
99957137 409 Int_t accepted=0;
c40cdd78 410 Int_t above=0;
411 for (Int_t i=0;i<npoints;i++) if (amp[index[i]]>0) above++;
99957137 412
46a0eaea 413 for (Int_t i=0;i<npoints;i++){
99957137 414 //
994b9770 415 if (dEdxMode==1 && amp[index[i]]==0) {
99957137 416 continue;
417 }
418 if (accepted<npoints*f0) continue;
419 if (accepted>npoints*f1) continue;
13f45dd8 420 sum0++;
421 sum1+= amp[index[i]];
422 sum2+= amp[index[i]];
99957137 423 accepted++;
13f45dd8 424 }
994b9770 425 if (dEdxMode==-1) return 1-Double_t(above)/Double_t(npoints);
13f45dd8 426 if (sum0<=0) return 0;
427 return sum1/sum0;
428}
429
52560b2e 430void AliTPCtrackFast::Simul(const char* fname, Int_t ntracks, Double_t diffFactor){
13f45dd8 431 //
432 //
433 //
434 AliTPCtrackFast fast;
994b9770 435 TTreeSRedirector *pcstream = new TTreeSRedirector(fname,"recreate");
13f45dd8 436 for (Int_t itr=0; itr<ntracks; itr++){
437 //
f0266161 438 fast.fMNprim=(10.+100*gRandom->Rndm());
439 if (gRandom->Rndm()>0.5) fast.fMNprim=1./(0.00001+gRandom->Rndm()*0.1);
c40cdd78 440
13f45dd8 441 fast.fDiff =0.01 +0.35*gRandom->Rndm();
52560b2e 442 // fast.fDiffLong = fast.fDiff*0.6/1.;
443 fast.fDiffLong = fast.fDiff*diffFactor/1.;
13f45dd8 444 //
445 fast.fAngleY = 4.0*(gRandom->Rndm()-0.5);
446 fast.fAngleZ = 4.0*(gRandom->Rndm()-0.5);
c40cdd78 447 fast.fN = 160;
13f45dd8 448 fast.MakeTrack();
449 if (itr%100==0) printf("%d\n",itr);
994b9770 450 (*pcstream)<<"simulTrack"<<
13f45dd8 451 "tr.="<<&fast<<
452 "\n";
453 }
454 fast.Write("track");
994b9770 455 delete pcstream;
13f45dd8 456}
457
458
459
0716cb63 460AliTPCclusterFast::AliTPCclusterFast(){
37045aeb 461 //
462 //
463 fDigits.ResizeTo(5,7);
0716cb63 464}
37045aeb 465
f0266161 466void AliTPCclusterFast::Init(){
467 //
468 // reset all counters
469 //
5fc70f07 470 const Int_t knMax=10000;
f0266161 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
476 //
477 fPosY.ResizeTo(knMax);
478 fPosZ.ResizeTo(knMax);
479 fGain.ResizeTo(knMax);
480 fSec.ResizeTo(knMax);
481 fStatY.ResizeTo(3);
482 fStatZ.ResizeTo(3);
483 for (Int_t i=0; i<knMax; i++){
484 fPosY[i]=0;
485 fPosZ[i]=0;
486 fGain[i]=0;
e05accd3 487 fSec[i]=0;
f0266161 488 }
489}
490
491
492
0716cb63 493AliTPCclusterFast::~AliTPCclusterFast(){
494}
495
496
f0266161 497void AliTPCclusterFast::SetParam(Float_t mnprim, Float_t diff, Float_t diffL,Float_t y, Float_t z, Float_t ky, Float_t kz){
0716cb63 498 //
499 //
f0266161 500 fMNprim = mnprim; fDiff = diff; fDiffLong=diffL;
0716cb63 501 fY=y; fZ=z;
502 fAngleY=ky; fAngleZ=kz;
503}
37045aeb 504Double_t AliTPCclusterFast::GetNsec(){
0716cb63 505 //
37045aeb 506 // Generate number of secondary electrons
507 // copy of procedure implemented in geant
0716cb63 508 //
e05accd3 509 const Double_t FPOT=20.77E-9, EEND=10E-6, EEXPO=2.2; // EEND1=1E-6;
0716cb63 510 const Double_t XEXPO=-EEXPO+1, YEXPO=1/XEXPO;
511 const Double_t W=20.77E-9;
37045aeb 512 Float_t RAN = gRandom->Rndm();
46a0eaea 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);
37045aeb 516 return TMath::Nint(TMath::Power((TMath::Power(FPOT,XEXPO)*(1-RAN)+TMath::Power(EEND,XEXPO)*RAN),YEXPO)/W);
517}
518
f0266161 519void AliTPCclusterFast::GenerElectrons(AliTPCclusterFast *cl0, AliTPCclusterFast *clm, AliTPCclusterFast *clp){
0716cb63 520 //
37045aeb 521 //
522 //
523 //
524 const Int_t knMax=1000;
f0266161 525 cl0->fNprim = gRandom->Poisson(cl0->fMNprim); //number of primary electrons
e05accd3 526 // cl0->fNtot=0; //total number of electrons
527 // cl0->fQtot=0; //total number of electrons after gain multiplification
0716cb63 528 //
529 Double_t sumQ=0;
530 Double_t sumYQ=0;
531 Double_t sumZQ=0;
532 Double_t sumY2Q=0;
533 Double_t sumZ2Q=0;
e05accd3 534 // for (Int_t i=0;i<knMax;i++){
535 // cl0->fSec[i]=0;
536 //}
f0266161 537 for (Int_t iprim=0; iprim<cl0->fNprim;iprim++){
538 Float_t dN = cl0->GetNsec();
539 cl0->fSec[iprim]=dN;
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);
543
37045aeb 544 for (Int_t isec=0;isec<=dN;isec++){
0716cb63 545 //
546 //
f0266161 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;
550 // choose pad row
551 AliTPCclusterFast *cl=cl0;
552 if (r<-0.5 &&cl) cl=clm;
553 if (r>0.5 &&cl) cl=clp;
554 //
0716cb63 555 Double_t gg = -TMath::Log(gRandom->Rndm());
f0266161 556 cl->fPosY[cl->fNtot]=y;
557 cl->fPosZ[cl->fNtot]=z;
558 cl->fGain[cl->fNtot]=gg;
559 cl->fQtot+=gg;
560 cl->fNtot++;
561 //
e05accd3 562 // cl->sumQ+=gg;
563 // cl->sumYQ+=gg*y;
564 // cl->sumY2Q+=gg*y*y;
565 // cl->sumZQ+=gg*z;
566 // cl->sumZ2Q+=gg*z*z;
f0266161 567 if (cl->fNtot>=knMax) continue;
0716cb63 568 }
f0266161 569 if (cl0->fNtot>=knMax) break;
0716cb63 570 }
f0266161 571
572 // if (sumQ>0){
573// fStatY[0]=sumQ;
574// fStatY[1]=sumYQ/sumQ;
575// fStatY[2]=sumY2Q/sumQ-fStatY[1]*fStatY[1];
576// fStatZ[0]=sumQ;
577// fStatZ[1]=sumZQ/sumQ;
578// fStatZ[2]=sumZ2Q/sumQ-fStatZ[1]*fStatZ[1];
579// }
0716cb63 580}
581
37045aeb 582void AliTPCclusterFast::Digitize(){
583 //
584 //
585 //
586 // 1. Clear digits
587 for (Int_t i=0; i<5;i++)
588 for (Int_t j=0; j<7;j++){
589 fDigits(i,j)=0;
590 }
591 //
592 // Fill digits
593 for (Int_t iel = 0; iel<fNtot; iel++){
594 for (Int_t di=-2; di<=2;di++)
13f45dd8 595 for (Int_t dj=-3; dj<=3;dj++){
37045aeb 596 Float_t fac = fPRF->Eval(di-fPosY[iel])*fTRF->Eval(dj-fPosZ[iel]);
597 fac*=fGain[iel];
598 fDigits(2+di,3+dj)+=fac;
599 }
600 }
fa9e5d19 601 //
602 //
603 //
37045aeb 604}
605
606
0716cb63 607
f0266161 608// void AliTPCclusterFast::Simul(const char* fname, Int_t npoints){
609// //
610// // Calc rms
611// //
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;
619// //
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();
624// fast.Digitize();
625// if (icl%10000==0) printf("%d\n",icl);
626// cstream<<"simul"<<
627// "s.="<<&fast<<
628// "\n";
629// }
630// }
0716cb63 631
632
13f45dd8 633Double_t AliTPCclusterFast::GetQtot(Float_t gain, Float_t thr, Float_t noise, Bool_t brounding, Bool_t baddPedestal){
a314e83c 634 //
635 //
636 //
637 Float_t sum =0;
13f45dd8 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);
a314e83c 644 if (amp>thr) sum+=amp;
645 }
13f45dd8 646 }
a314e83c 647 return sum;
648}
649
13f45dd8 650Double_t AliTPCclusterFast::GetQmax(Float_t gain, Float_t thr, Float_t noise, Bool_t brounding, Bool_t baddPedestal){
a314e83c 651 //
652 //
653 //
654 Float_t max =0;
13f45dd8 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;
a314e83c 662 }
13f45dd8 663 }
a314e83c 664 return max;
665}
666
667
a314e83c 668
13f45dd8 669Double_t AliTPCclusterFast::GetQmaxCorr(Float_t rmsy0, Float_t rmsz0){
670 //
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
675 //
676 //
677 //
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);
681}
682
683
684Double_t AliTPCclusterFast::GetQtotCorr(Float_t rmsy0, Float_t rmsz0, Float_t gain, Float_t thr){
685 //
686 // Calculates the fraction of the charge over threshol to total charge
687 // The response function
688 //
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
5fc70f07 693 Double_t qmax = GetQmax(gain,thr,0); // qmax
13f45dd8 694 //
695 Double_t corr =1;
696 Double_t qnorm=qtot;
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;
703 sumAll+=qlocal;
704 }
705 if (sumAll>0&&sumThr>0) corr=(sumThr)/sumAll;
5fc70f07 706 if (sumThr==0) corr=GetQmaxCorr(0.5,0.5);
13f45dd8 707 //corr = sumThr;
708 if (corr>0) qnorm=qtot/corr;
709
710 }
711 return corr;
712}
a314e83c 713
714
0716cb63 715
13f45dd8 716
717
718Double_t AliTPCclusterFast::GaussConvolution(Double_t x0, Double_t x1, Double_t k0, Double_t k1, Double_t s0, Double_t s1){
719 //
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}]]
723 //
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)
726 //
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());
731 return val;
732 }
733
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));
736 //
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);
743 return val;
744
745}
746
747
748Double_t AliTPCclusterFast::GaussExpConvolution(Double_t x0, Double_t s0,Double_t t1){
749 //
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;
759 val *=t1/(2.);
760 return val;
761
762}
763
764
765Double_t AliTPCclusterFast::Gamma4(Double_t x, Double_t p0, Double_t p1){
766 //
767 // Gamma 4 Time response function of ALTRO
768 //
769 if (x<0) return 0;
770 Double_t g1 = TMath::Exp(-4.*x/p1);
771 Double_t g2 = TMath::Power(x/p1,4);
772 return p0*g1*g2;
773}
774
775
776
777Double_t AliTPCclusterFast::GaussGamma4(Double_t x, Double_t s0, Double_t p1){
778 //
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)
782
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));
46a0eaea 786 // Double_t xp14 = TMath::Power(TMath::Abs((x/p1)),4);
13f45dd8 787 return exp1*erf1;
788
789
790}
791
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}]]
794//
795//
796// No analytical solution
797//
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}]]