]>
Commit | Line | Data |
---|---|---|
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" |
75bbe354 | 46 | #include "TGrid.h" |
0716cb63 | 47 | |
48 | class AliTPCclusterFast: public TObject { | |
49 | public: | |
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 | 68 | public: |
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 | ||
101 | class AliTPCtrackFast: public TObject { | |
102 | public: | |
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 | 133 | ClassImp(AliTPCclusterFast) |
13f45dd8 | 134 | ClassImp(AliTPCtrackFast) |
135 | ||
136 | ||
137 | ||
138 | ||
0716cb63 | 139 | |
37045aeb | 140 | TF1 *AliTPCclusterFast::fPRF=0; |
141 | TF1 *AliTPCclusterFast::fTRF=0; | |
0716cb63 | 142 | |
143 | ||
13f45dd8 | 144 | AliTPCtrackFast::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 | ||
158 | void AliTPCtrackFast::Add(AliTPCtrackFast &track2){ | |
159 | if (!track2.fInit) return; | |
160 | ||
13f45dd8 | 161 | } |
162 | ||
13f45dd8 | 163 | |
13f45dd8 | 164 | |
165 | void AliTPCtrackFast::MakeTrack(){ | |
166 | // | |
167 | // | |
168 | // | |
46a0eaea | 169 | if (!fCl) fCl = new TClonesArray("AliTPCclusterFast",160); |
75bbe354 | 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 | } | |
75bbe354 | 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); |
75bbe354 | 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 | } | |
75bbe354 | 202 | |
13f45dd8 | 203 | } |
204 | ||
205 | Double_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 | ||
217 | Double_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 | 230 | Double_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 | 276 | Double_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 | 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){ |
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 | 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){ |
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 | 386 | Double_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 | 430 | void 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 | 460 | AliTPCclusterFast::AliTPCclusterFast(){ |
37045aeb | 461 | // |
462 | // | |
463 | fDigits.ResizeTo(5,7); | |
0716cb63 | 464 | } |
37045aeb | 465 | |
f0266161 | 466 | void AliTPCclusterFast::Init(){ |
467 | // | |
468 | // reset all counters | |
469 | // | |
ceb30540 | 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; | |
75bbe354 | 487 | fSec[i]=0; |
f0266161 | 488 | } |
489 | } | |
490 | ||
491 | ||
492 | ||
0716cb63 | 493 | AliTPCclusterFast::~AliTPCclusterFast(){ |
494 | } | |
495 | ||
496 | ||
f0266161 | 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){ |
0716cb63 | 498 | // |
499 | // | |
f0266161 | 500 | fMNprim = mnprim; fDiff = diff; fDiffLong=diffL; |
0716cb63 | 501 | fY=y; fZ=z; |
502 | fAngleY=ky; fAngleZ=kz; | |
503 | } | |
37045aeb | 504 | Double_t AliTPCclusterFast::GetNsec(){ |
0716cb63 | 505 | // |
37045aeb | 506 | // Generate number of secondary electrons |
507 | // copy of procedure implemented in geant | |
0716cb63 | 508 | // |
75bbe354 | 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 | 519 | void 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 |
75bbe354 | 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; | |
75bbe354 | 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 | // | |
75bbe354 | 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 | 582 | void 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 | 633 | Double_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 | 650 | Double_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 | 669 | Double_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 | ||
684 | Double_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 | |
ceb30540 | 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; | |
ceb30540 | 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 | ||
718 | Double_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 | ||
748 | Double_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 | ||
765 | Double_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 | ||
777 | Double_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}]] |