]>
Commit | Line | Data |
---|---|---|
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 | ||
28 | class AliTPCclusterFast: public TObject { | |
29 | public: | |
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 | 46 | public: |
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 | ||
78 | class AliTPCtrackFast: public TObject { | |
79 | public: | |
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 | 118 | ClassImp(AliTPCclusterFast) |
13f45dd8 | 119 | ClassImp(AliTPCtrackFast) |
120 | ||
121 | ||
122 | ||
123 | ||
0716cb63 | 124 | |
37045aeb | 125 | TF1 *AliTPCclusterFast::fPRF=0; |
126 | TF1 *AliTPCclusterFast::fTRF=0; | |
0716cb63 | 127 | |
128 | ||
13f45dd8 | 129 | AliTPCtrackFast::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 | ||
150 | void 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 | ||
163 | void 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 | ||
203 | void 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 | ||
244 | void 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 | ||
265 | Double_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 | ||
275 | Double_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 | 285 | Double_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 | ||
304 | Double_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 | ||
324 | Double_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 | ||
342 | void 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 | 367 | AliTPCclusterFast::AliTPCclusterFast(){ |
37045aeb | 368 | // |
369 | // | |
370 | fDigits.ResizeTo(5,7); | |
0716cb63 | 371 | } |
37045aeb | 372 | |
0716cb63 | 373 | AliTPCclusterFast::~AliTPCclusterFast(){ |
374 | } | |
375 | ||
376 | ||
377 | void 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 | 384 | Double_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 | ||
399 | void 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 | 460 | void 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 | |
484 | void 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 | 509 | Double_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 | 526 | Double_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 | 545 | Double_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 | ||
560 | Double_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 | ||
592 | Double_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 | ||
622 | Double_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 | ||
639 | Double_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 | ||
651 | Double_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}]] |