]> git.uio.no Git - u/mrichter/AliRoot.git/blame - FASTSIM/AliMUONFastTracking.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Transition to NewIO
[u/mrichter/AliRoot.git] / FASTSIM / AliMUONFastTracking.cxx
CommitLineData
6255180c 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/*
17$Log$
88cb7938 18Revision 1.2.2.1 2003/04/15 15:57:31 hristov
19Merging with v3-09-08
20
21Revision 1.2 2003/01/08 10:29:33 morsch
22Path to data file changed.
23
f03af365 24Revision 1.1 2003/01/06 10:13:33 morsch
25First commit.
26
6255180c 27*/
28
29#include "AliMUONFastTracking.h"
30#include "AliMUONFastTrackingEntry.h"
31#include <TMatrixD.h>
32#include <TSpline.h>
33#include <TFile.h>
34#include <TH1.h>
35#include <TH3.h>
36#include <TF1.h>
37#include <TRandom.h>
38#include <stdlib.h>
39#include <stdio.h>
40#include <string.h>
41#include <Riostream.h>
42
43ClassImp(AliMUONFastTracking)
44
45AliMUONFastTracking* AliMUONFastTracking::fgMUONFastTracking=NULL;
46
47static Double_t FitP(Double_t *x, Double_t *par){
48 Double_t dx = x[0] - par[0];
49 Double_t dx2 = x[0] - par[4];
50 Double_t sigma = par[1] * ( 1 + par[2] * dx);
51 if (sigma == 0) {
52
53 return 0.;
54 }
55 Double_t fasymm = TMath::Exp(-0.5 * dx * dx / (sigma * sigma));
56 Double_t sigma2 = par[1] * par[5];
57 Double_t fgauss = TMath::Exp(-0.5 * dx2 * dx2 / (sigma2 * sigma2));
58 return TMath::Abs(fasymm + par[3] * fgauss);
59}
60
61AliMUONFastTracking* AliMUONFastTracking::Instance()
62{
63// Set random number generator
64 if (fgMUONFastTracking) {
65 return fgMUONFastTracking;
66 } else {
67 fgMUONFastTracking = new AliMUONFastTracking();
68 return fgMUONFastTracking;
69 }
70}
71
72AliMUONFastTracking::AliMUONFastTracking()
73{
74// SetBackground();
75
76 fPrintLevel = 1;
77 // read binning; temporarily put by hand
78 Float_t pmin = 0, pmax = 200;
79 Int_t nbinp = 10;
80 Float_t thetamin = 2, thetamax = 9;
81 Int_t nbintheta=10;
82 Float_t phimin = -180, phimax =180;
83 Int_t nbinphi=10;
84 //--------------------------------------
85
86 fNbinp = nbinp;
87 fPmin = pmin;
88 fPmax = pmax;
89
90 fNbintheta = nbintheta;
91 fThetamin = thetamin;
92 fThetamax = thetamax;
93
94 fNbinphi = nbinphi;
95 fPhimin = phimin;
96 fPhimax = phimax;
97
98 fDeltaP = (fPmax-fPmin)/fNbinp;
99 fDeltaTheta = (fThetamax-fThetamin)/fNbintheta;
100 fDeltaPhi = (fPhimax-fPhimin)/fNbinphi;
101}
102
103void AliMUONFastTracking::Init(Float_t bkg)
104{
105 //
106 // Initialization
107 //
108 for (Int_t ip=0; ip< fNbinp; ip++){
109 for (Int_t itheta=0; itheta< fNbintheta; itheta++){
110 for (Int_t iphi=0; iphi< fNbinphi; iphi++){
111 fCurrentEntry[ip][itheta][iphi] = new AliMUONFastTrackingEntry;
112 for (Int_t ibkg=0; ibkg<4; ibkg++){
113 fEntry[ip][itheta][iphi][ibkg] = new AliMUONFastTrackingEntry;
114 }
115 }
116 }
117 }
118
119 char filename [100];
f03af365 120 sprintf (filename,"$(ALICE_ROOT)/FASTSIM/data/MUONtrackLUT.root");
6255180c 121 TFile *file = new TFile(filename);
122 ReadLUT(file);
123 SetBackground(bkg);
124 UseSpline(1);
125 fFitp = new TF1("fit1",FitP,-20.,20.,6);
126 fFitp->SetNpx(200);
127}
128
129
130void AliMUONFastTracking::ReadLUT(TFile* file)
131{
132 TH3F *heff[5][3], *hacc[5][3], *hmeanp, *hsigmap, *hsigma1p, *hchi2p;
133 TH3F *hnormg2, *hmeang2, *hsigmag2, *hmeantheta, *hsigmatheta, *hchi2theta;
134 TH3F *hmeanphi, *hsigmaphi, *hchi2phi;
135 char tag[40], tag2[40];
136
137 const Float_t bkg[4] = {0, 0.5, 1, 2};
138 for (Int_t ibkg=0; ibkg<4; ibkg++) {
139 sprintf (tag,"BKG%g",bkg[ibkg]);
140 file->cd(tag);
141 for (Int_t isplp = 0; isplp<kSplitP; isplp++) {
142 for (Int_t ispltheta = 0; ispltheta<kSplitTheta; ispltheta++) {
143 sprintf (tag2,"heff[%d][%d]",isplp,ispltheta);
144 heff[isplp][ispltheta] = (TH3F*)gDirectory->Get(tag2);
145 sprintf (tag2,"hacc[%d][%d]",isplp,ispltheta);
146 hacc[isplp][ispltheta] = (TH3F*)gDirectory->Get(tag2);
147 }
148 }
149 hmeanp = (TH3F*)gDirectory->Get("hmeanp");
150 hsigmap = (TH3F*)gDirectory->Get("hsigmap");
151 hsigma1p = (TH3F*)gDirectory->Get("hsigma1p");
152 hchi2p = (TH3F*)gDirectory->Get("hchi2p");
153 hnormg2 = (TH3F*)gDirectory->Get("hnormg2");
154 hmeang2 = (TH3F*)gDirectory->Get("hmeang2");
155 hsigmag2 = (TH3F*)gDirectory->Get("hsigmag2");
156 hmeantheta = (TH3F*)gDirectory->Get("hmeantheta");
157 hsigmatheta = (TH3F*)gDirectory->Get("hsigmatheta");
158 hchi2theta = (TH3F*)gDirectory->Get("hchi2theta");
159 hmeanphi = (TH3F*)gDirectory->Get("hmeanphi");
160 hsigmaphi = (TH3F*)gDirectory->Get("hsigmaphi");
161 hchi2phi = (TH3F*)gDirectory->Get("hchi2phi");
162
163 printf ("Reading parameters from LUT file %s...\n",file->GetName());
164 for (Int_t ip=0; ip<fNbinp ;ip++) {
165 for (Int_t itheta=0; itheta<fNbintheta ;itheta++) {
166 for (Int_t iphi=0; iphi<fNbinphi ;iphi++) {
167 Float_t p = fPmin + fDeltaP * (ip + 0.5);
168 Float_t theta = fThetamin + fDeltaTheta * (itheta + 0.5);
169 Float_t phi = fPhimin + fDeltaPhi * (iphi + 0.5);
170
171 fEntry[ip][itheta][iphi][ibkg]->fP = p;
172 fEntry[ip][itheta][iphi][ibkg]->fMeanp =
173 hmeanp->GetBinContent(ip+1,itheta+1,iphi+1);
174 fEntry[ip][itheta][iphi][ibkg]->fSigmap =
175 TMath::Abs(hsigmap->GetBinContent(ip+1,itheta+1,iphi+1));
176 fEntry[ip][itheta][iphi][ibkg]->fSigma1p =
177 hsigma1p->GetBinContent(ip+1,itheta+1,iphi+1);
178 fEntry[ip][itheta][iphi][ibkg]->fChi2p =
179 hchi2p->GetBinContent(ip+1,itheta+1,iphi+1);
180 fEntry[ip][itheta][iphi][ibkg]->fNormG2 =
181 hnormg2->GetBinContent(ip+1,itheta+1,iphi+1);
182 fEntry[ip][itheta][iphi][ibkg]->fMeanG2 =
183 hmeang2->GetBinContent(ip+1,itheta+1,iphi+1);
184 if (ibkg == 0) fEntry[ip][itheta][iphi][ibkg]->fSigmaG2 = 9999;
185 else fEntry[ip][itheta][iphi][ibkg]->fSigmaG2 =
186 hsigmag2->GetBinContent(ip+1,itheta+1,iphi+1);
187 fEntry[ip][itheta][iphi][ibkg]->fTheta = theta;
188 fEntry[ip][itheta][iphi][ibkg]->fMeantheta =
189 hmeantheta->GetBinContent(ip+1,itheta+1,iphi+1);
190 fEntry[ip][itheta][iphi][ibkg]->fSigmatheta =
191 TMath::Abs(hsigmatheta->GetBinContent(ip+1,itheta+1,iphi+1));
192 fEntry[ip][itheta][iphi][ibkg]->fChi2theta =
193 hchi2theta->GetBinContent(ip+1,itheta+1,iphi+1);
194 fEntry[ip][itheta][iphi][ibkg]->fPhi = phi;
195 fEntry[ip][itheta][iphi][ibkg]->fMeanphi =
196 hmeanphi->GetBinContent(ip+1,itheta+1,iphi+1);
197 fEntry[ip][itheta][iphi][ibkg]->fSigmaphi =
198 TMath::Abs(hsigmaphi->GetBinContent(ip+1,itheta+1,iphi+1));
199 fEntry[ip][itheta][iphi][ibkg]->fChi2phi =
200 hchi2phi->GetBinContent(ip+1,itheta+1,iphi+1);
201 for (Int_t i=0; i<kSplitP; i++) {
202 for (Int_t j=0; j<kSplitTheta; j++) {
203 fEntry[ip][itheta][iphi][ibkg]->fAcc[i][j] =
204 hacc[i][j]->GetBinContent(ip+1,itheta+1,iphi+1);
205 fEntry[ip][itheta][iphi][ibkg]->fEff[i][j] =
206 heff[i][j]->GetBinContent(ip+1,itheta+1,iphi+1);
207 }
208 }
209 } // iphi
210 } // itheta
211 } // ip
212 } // ibkg
213
214 TGraph *graph = new TGraph(3);
215 TF1 *f = new TF1("f","[0]+[1]*x");
216
217 for (Int_t ip=0; ip< fNbinp; ip++){
218 for (Int_t itheta=0; itheta< fNbintheta; itheta++){
219 for (Int_t iphi=0; iphi< fNbinphi; iphi++){
220 graph->SetPoint(0,0.5,fEntry[ip][itheta][iphi][1]->fSigmaG2);
221 graph->SetPoint(1,1,fEntry[ip][itheta][iphi][2]->fSigmaG2);
222 graph->SetPoint(2,2,fEntry[ip][itheta][iphi][3]->fSigmaG2);
223 graph->Fit("f","q");
224 fEntry[ip][itheta][iphi][0]->fSigmaG2 = f->Eval(0);
225 }
226 }
227 }
228 f->Delete();
229 graph->Delete();
230 printf ("parameters read. \n");
231}
232
233void AliMUONFastTracking::GetBinning(Int_t &nbinp, Float_t &pmin, Float_t &pmax,
234 Int_t &nbintheta, Float_t &thetamin,
235 Float_t &thetamax,
236 Int_t &nbinphi, Float_t &phimin, Float_t &phimax)
237{
238 nbinp = fNbinp;
239 pmin = fPmin;
240 pmax = fPmax;
241 nbintheta = fNbintheta;
242 thetamin = fThetamin;
243 thetamax = fThetamax;
244 nbinphi = fNbinphi;
245 phimin = fPhimin;
246 phimax = fPhimax;
247}
248
249
250void AliMUONFastTracking::GetIpIthetaIphi(Float_t p, Float_t theta, Float_t phi,
251 Int_t charge, Int_t &ip, Int_t &itheta,
252 Int_t &iphi)
253{
254 if (charge < 0) phi = -phi;
255 ip = Int_t (( p - fPmin ) / fDeltaP);
256 itheta = Int_t (( theta - fThetamin ) / fDeltaTheta);
257 iphi = Int_t (( phi - fPhimin ) / fDeltaPhi);
258
259 if (ip< 0) {
260 printf ("Warning: ip= %d. Set to 0\n",ip);
261 ip = 0;
262 }
263 if (ip>= fNbinp) {
264 // printf ("Warning: ip = %d. Set to %d\n",ip,fNbinp-1);
265 ip = fNbinp-1;
266 }
267 if (itheta< 0) {
268 // printf ("Warning: itheta= %d. Set to 0\n",itheta);
269 itheta = 0;
270 }
271 if (itheta>= fNbintheta) {
272 // printf ("Warning: itheta = %d. Set to %d\n",itheta,fNbintheta-1);
273 itheta = fNbintheta-1;
274 }
275
276 if (iphi< 0) {
277 printf ("Warning: iphi= %d. Set to 0\n",iphi);
278 iphi = 0;
279 }
280 if (iphi>= fNbinphi) {
281 printf ("Warning: iphi = %d. Set to %d\n",iphi,fNbinphi-1);
282 iphi = fNbinphi-1;
283 }
284}
285
286void AliMUONFastTracking::GetSplit(Int_t ip, Int_t itheta,
287 Int_t &nSplitP, Int_t &nSplitTheta) {
288 if (ip==0) nSplitP = 5;
289 else nSplitP = 2;
290 if (itheta==0) nSplitTheta = 3;
291 else nSplitTheta = 1;
292}
293
294Float_t AliMUONFastTracking::Efficiency(Float_t p, Float_t theta,
295 Float_t phi, Int_t charge){
296 Int_t ip=0, itheta=0, iphi=0;
297 GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
298 Int_t nSplitP, nSplitTheta;
299 GetSplit(ip,itheta,nSplitP,nSplitTheta);
300
301 Float_t dp = p - fPmin;
302 Int_t ibinp = Int_t(nSplitP*(dp - fDeltaP * Int_t(dp / fDeltaP))/fDeltaP);
303 Float_t dtheta = theta - fThetamin;
304 Int_t ibintheta = Int_t(nSplitTheta*(dtheta - fDeltaTheta * Int_t(dtheta / fDeltaTheta))/fDeltaTheta);
305 Float_t eff = fCurrentEntry[ip][itheta][iphi]->fEff[ibinp][ibintheta];
306 return eff;
307}
308
309Float_t AliMUONFastTracking::Acceptance(Float_t p, Float_t theta,
310 Float_t phi, Int_t charge){
311 if (theta<fThetamin || theta>fThetamax) return 0;
312
313 Int_t ip=0, itheta=0, iphi=0;
314 GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
315 Int_t nSplitP, nSplitTheta;
316 GetSplit(ip,itheta,nSplitP,nSplitTheta);
317 // central value and corrections with spline
318
319 Float_t dp = p - fPmin;
320 Int_t ibinp = Int_t(nSplitP*(dp - fDeltaP * Int_t(dp / fDeltaP))/fDeltaP);
321 Float_t dtheta = theta - fThetamin;
322 Int_t ibintheta = Int_t(nSplitTheta*(dtheta - fDeltaTheta * Int_t(dtheta / fDeltaTheta))/fDeltaTheta);
323 Float_t acc = fCurrentEntry[ip][itheta][iphi]->fAcc[ibinp][ibintheta];
324 return acc;
325}
326
327Float_t AliMUONFastTracking::MeanP(Float_t p, Float_t theta,
328 Float_t phi, Int_t charge)
329{
330 Int_t ip=0, itheta=0, iphi=0;
331 GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
332 return fCurrentEntry[ip][itheta][iphi]->fMeanp;
333}
334
335Float_t AliMUONFastTracking::SigmaP(Float_t p, Float_t theta,
336 Float_t phi, Int_t charge)
337{
338 Int_t ip=0, itheta=0, iphi=0;
339 Int_t index;
340 GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
341 // central value and corrections with spline
342 Float_t sigmap = fCurrentEntry[ip][itheta][iphi]->fSigmap;
343 if (!fSpline) return sigmap;
344 // corrections vs p, theta, phi
345 index = iphi + fNbinphi * itheta;
346 Double_t xmin,ymin,xmax,ymax;
347 Float_t frac1 = fSplineSigmap[index][0]->Eval(p)/sigmap;
348
349 if (p>fPmax-fDeltaP/2.) {
350 Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->fSigmap;
351 Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->fSigmap;
352 Float_t s3 = fCurrentEntry[fNbinp-3][itheta][iphi]->fSigmap;
353 Float_t p1 = fDeltaP * (fNbinp - 1 + 0.5) + fPmin;
354 Float_t p2 = fDeltaP * (fNbinp - 2 + 0.5) + fPmin;
355 Float_t p3 = fDeltaP * (fNbinp - 3 + 0.5) + fPmin;
356 Float_t p12 = p1 * p1, p22 = p2 * p2, p32 = p3 * p3;
357 Float_t d = p12*p2 + p1*p32 + p22*p3 - p32*p2 - p3*p12 - p22*p1;
358 Float_t a = (s1*p2 + p1*s3 + s2*p3 - s3*p2 - p3*s1 - s2*p1) / d;
359 Float_t b = (p12*s2 + s1*p32 + p22*s3 - p32*s2 - s3*p12 - p22*s1)/d;
360 Float_t c = (p12*p2*s3 + p1*p32*s2 + p22*p3*s1
361 - p32*p2*s1 - p3*p12*s2 - p22*p1*s3) / d;
362 Float_t sigma = a * p * p + b * p + c;
363 frac1 = sigma/sigmap;
364 }
365 index = iphi + fNbinphi * ip;
366 fSplineEff[index][1]->GetKnot(0,xmin,ymin);
367 fSplineEff[index][1]->GetKnot(9,xmax,ymax);
368 if (theta>xmax) theta = xmax;
369 Float_t frac2 = fSplineSigmap[index][1]->Eval(theta)/sigmap;
370 index = itheta + fNbintheta * ip;
371 fSplineEff[index][2]->GetKnot(0,xmin,ymin);
372 fSplineEff[index][2]->GetKnot(9,xmax,ymax);
373 if (phi>xmax) phi = xmax;
374 Float_t frac3 = fSplineSigmap[index][2]->Eval(phi)/sigmap;
375 Float_t sigmatot = sigmap * frac1 * frac2 * frac3;
376 if (sigmatot<0) sigmatot = sigmap;
377 return sigmatot;
378}
379
380Float_t AliMUONFastTracking::Sigma1P(Float_t p, Float_t theta,
381 Float_t phi, Int_t charge)
382{
383 Int_t ip=0, itheta=0, iphi=0;
384 GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
385 if (p>fPmax) {
386 // linear extrapolation of sigmap for p out of range
387 Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->fSigma1p;
388 Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->fSigma1p;
389 Float_t p1 = fDeltaP * (fNbinp - 1 + 0.5) + fPmin;
390 Float_t p2 = fDeltaP * (fNbinp - 2 + 0.5) + fPmin;
391 Float_t sigma = 1./(p1-p2) * ( (s1-s2)*p + (s2-s1)*p1 + s1*(p1-p2) );
392 return sigma;
393 }
394 else return fCurrentEntry[ip][itheta][iphi]->fSigma1p;
395}
396
397Float_t AliMUONFastTracking::NormG2(Float_t p, Float_t theta,
398 Float_t phi, Int_t charge)
399{
400 Int_t ip=0, itheta=0, iphi=0;
401 GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
402 if (p>fPmax) {
403 // linear extrapolation of sigmap for p out of range
404 Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->fNormG2;
405 Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->fNormG2;
406 Float_t p1 = fDeltaP * (fNbinp - 1 + 0.5) + fPmin;
407 Float_t p2 = fDeltaP * (fNbinp - 2 + 0.5) + fPmin;
408 Float_t norm = 1./(p1-p2) * ( (s1-s2)*p + (s2-s1)*p1 + s1*(p1-p2) );
409 return norm;
410 }
411 else return fCurrentEntry[ip][itheta][iphi]->fNormG2;
412}
413
414Float_t AliMUONFastTracking::MeanG2(Float_t p, Float_t theta,
415 Float_t phi, Int_t charge)
416{
417 Int_t ip=0, itheta=0, iphi=0;
418 GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
419 if (p>fPmax) {
420 // linear extrapolation of sigmap for p out of range
421 Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->fMeanG2;
422 Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->fMeanG2;
423 Float_t p1 = fDeltaP * (fNbinp - 1 + 0.5) + fPmin;
424 Float_t p2 = fDeltaP * (fNbinp - 2 + 0.5) + fPmin;
425 Float_t norm = 1./(p1-p2) * ( (s1-s2)*p + (s2-s1)*p1 + s1*(p1-p2) );
426 return norm;
427 }
428 else return fCurrentEntry[ip][itheta][iphi]->fMeanG2;
429}
430
431Float_t AliMUONFastTracking::SigmaG2(Float_t p, Float_t theta,
432 Float_t phi, Int_t charge)
433{
434 Int_t ip=0, itheta=0, iphi=0;
435 GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
436 if (p>fPmax) {
437 // linear extrapolation of sigmap for p out of range
438 Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->fSigmaG2;
439 Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->fSigmaG2;
440 Float_t p1 = fDeltaP * (fNbinp - 1 + 0.5) + fPmin;
441 Float_t p2 = fDeltaP * (fNbinp - 2 + 0.5) + fPmin;
442 Float_t sigma = 1./(p1-p2) * ( (s1-s2)*p + (s2-s1)*p1 + s1*(p1-p2) );
443 return sigma;
444 }
445 else return fCurrentEntry[ip][itheta][iphi]->fSigmaG2;
446}
447
448
449Float_t AliMUONFastTracking::MeanTheta(Float_t p, Float_t theta,
450 Float_t phi, Int_t charge)
451{
452 Int_t ip=0, itheta=0, iphi=0;
453 GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
454 return fCurrentEntry[ip][itheta][iphi]->fMeantheta;
455}
456
457Float_t AliMUONFastTracking::SigmaTheta(Float_t p, Float_t theta,
458 Float_t phi, Int_t charge){
459 Int_t ip=0, itheta=0, iphi=0;
460 Int_t index;
461 GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
462 // central value and corrections with spline
463 Float_t sigmatheta = fCurrentEntry[ip][itheta][iphi]->fSigmatheta;
464 if (!fSpline) return sigmatheta;
465 // corrections vs p, theta, phi
466 index = iphi + fNbinphi * itheta;
467 Double_t xmin,ymin,xmax,ymax;
468 Float_t frac1 = fSplineSigmatheta[index][0]->Eval(p)/sigmatheta;
469 if (p>fPmax-fDeltaP/2.) {
470 // linear extrapolation of sigmap for p out of range
471 Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->fSigmatheta;
472 Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->fSigmatheta;
473 Float_t p1 = fDeltaP * (fNbinp - 1 + 0.5) + fPmin;
474 Float_t p2 = fDeltaP * (fNbinp - 2 + 0.5) + fPmin;
475 Float_t sigma = 1./(p1-p2) * ( (s1-s2)*p + (s2-s1)*p1 + s1*(p1-p2) );
476 frac1=sigma/sigmatheta;
477 }
478 index = iphi + fNbinphi * ip;
479 fSplineEff[index][1]->GetKnot(0,xmin,ymin);
480 fSplineEff[index][1]->GetKnot(9,xmax,ymax);
481 if (theta>xmax) theta = xmax;
482 Float_t frac2 = fSplineSigmatheta[index][1]->Eval(theta)/sigmatheta;
483 index = itheta + fNbintheta * ip;
484 fSplineEff[index][2]->GetKnot(0,xmin,ymin);
485 fSplineEff[index][2]->GetKnot(9,xmax,ymax);
486 if (phi>xmax) phi = xmax;
487 Float_t frac3 = fSplineSigmatheta[index][2]->Eval(phi)/sigmatheta;
488 return sigmatheta * frac1 * frac2 * frac3;
489}
490
491
492Float_t AliMUONFastTracking::MeanPhi(Float_t p, Float_t theta,
493 Float_t phi, Int_t charge){
494 Int_t ip=0, itheta=0, iphi=0;
495 GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
496 return fCurrentEntry[ip][itheta][iphi]->fMeanphi;
497}
498
499Float_t AliMUONFastTracking::SigmaPhi(Float_t p, Float_t theta,
500 Float_t phi, Int_t charge){
501 Int_t ip=0, itheta=0, iphi=0;
502 Int_t index;
503 GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
504 // central value and corrections with spline
505 Float_t sigmaphi = fCurrentEntry[ip][itheta][iphi]->fSigmaphi;
506 if (!fSpline) return sigmaphi;
507 // corrections vs p, theta, phi
508 index = iphi + fNbinphi * itheta;
509 Float_t frac1 = fSplineSigmaphi[index][0]->Eval(p)/sigmaphi;
510 Double_t xmin,ymin,xmax,ymax;
511 if (p>fPmax-fDeltaP/2.) {
512 Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->fSigmaphi;
513 Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->fSigmaphi;
514 Float_t p1 = fDeltaP * (fNbinp - 1 + 0.5) + fPmin;
515 Float_t p2 = fDeltaP * (fNbinp - 2 + 0.5) + fPmin;
516 Float_t sigma = 1./(p1-p2) * ( (s1-s2)*p + (s2-s1)*p1 + s1*(p1-p2) );
517 frac1 = sigma/sigmaphi;
518 }
519
520 index = iphi + fNbinphi * ip;
521 fSplineEff[index][1]->GetKnot(0,xmin,ymin);
522 fSplineEff[index][1]->GetKnot(9,xmax,ymax);
523 if (theta>xmax) theta = xmax;
524 Float_t frac2 = fSplineSigmaphi[index][1]->Eval(theta)/sigmaphi;
525 index = itheta + fNbintheta * ip;
526 fSplineEff[index][2]->GetKnot(0,xmin,ymin);
527 fSplineEff[index][2]->GetKnot(9,xmax,ymax);
528 if (phi>xmax) phi = xmax;
529 Float_t frac3 = fSplineSigmaphi[index][2]->Eval(phi)/sigmaphi;
530 return sigmaphi * frac1 * frac2 * frac3;
531}
532
533void AliMUONFastTracking::SetSpline(){
534 printf ("Setting spline functions...");
535 char splname[40];
536 Double_t x[20][3];
537 Double_t x2[50][3];
538 Int_t nbins[3] = {fNbinp, fNbintheta, fNbinphi};
539 Double_t xspl[20],yeff[50],ysigmap[20],ysigma1p[20];
540 Double_t yacc[50], ysigmatheta[20],ysigmaphi[20];
541 Double_t xsp2[50];
542 // let's calculate the x axis for p, theta, phi
543
544 Int_t i, ispline, ivar;
545 for (i=0; i< fNbinp; i++) x[i][0] = fPmin + fDeltaP * (i + 0.5);
546 for (i=0; i< fNbintheta; i++) x[i][1] = fThetamin + fDeltaTheta * (i + 0.5);
547 for (i=0; i< fNbinphi; i++) x[i][2] = fPhimin + fDeltaPhi * (i + 0.5);
548
549 for (i=0; i< 5 * fNbinp; i++) x2[i][0] = fPmin + fDeltaP * (i + 0.5)/5.;
550 for (i=0; i< 5 * fNbintheta; i++) x2[i][1] = fThetamin + fDeltaTheta * (i + 0.5)/5.;
551 for (i=0; i< 5 * fNbinphi; i++) x2[i][2] = fPhimin + fDeltaPhi * (i + 0.5)/5.;
552
553 // splines in p
554 ivar = 0;
555 for (i=0; i<nbins[ivar]; i++) xspl[i] = x[i][ivar];
556 for (i=0; i<5 * nbins[ivar]; i++) xsp2[i] = x2[i][ivar];
557 ispline=0;
558 for (Int_t itheta=0; itheta< fNbintheta; itheta++){
559 for (Int_t iphi=0; iphi< fNbinphi; iphi++){
560 for (Int_t ip=0; ip<fNbinp; ip++) {
561 // for (Int_t i=0; i<5; i++) {
562 // yeff[5 * ip + i] = fCurrentEntry[ip][itheta][iphi]->fEff[i];
563 // yacc[5 * ip + i] = fCurrentEntry[ip][itheta][iphi]->fAcc[i];
564 // }
565 ysigmap[ip] = fCurrentEntry[ip][itheta][iphi]->fSigmap;
566 ysigma1p[ip] = fCurrentEntry[ip][itheta][iphi]->fSigma1p;
567 ysigmatheta[ip] = fCurrentEntry[ip][itheta][iphi]->fSigmatheta;
568 ysigmaphi[ip] = fCurrentEntry[ip][itheta][iphi]->fSigmaphi;
569 }
570 if (fPrintLevel>3) cout << " creating new spline " << splname << endl;
571 sprintf (splname,"fSplineEff[%d][%d]",ispline,ivar);
572 fSplineEff[ispline][ivar] = new TSpline3(splname,xsp2,yeff,5 * nbins[ivar]);
573 sprintf (splname,"fSplineAcc[%d][%d]",ispline,ivar);
574 fSplineAcc[ispline][ivar] = new TSpline3(splname,xsp2,yacc,5 * nbins[ivar]);
575 sprintf (splname,"fSplineSigmap[%d][%d]",ispline,ivar);
576 fSplineSigmap[ispline][ivar] = new TSpline3(splname,xspl,ysigmap,nbins[ivar]);
577 sprintf (splname,"fSplineSigma1p[%d][%d]",ispline,ivar);
578 fSplineSigma1p[ispline][ivar] = new TSpline3(splname,xspl,ysigma1p,nbins[ivar]);
579 sprintf (splname,"fSplineSigmatheta[%d][%d]",ispline,ivar);
580 fSplineSigmatheta[ispline][ivar] = new TSpline3(splname,xspl,ysigmatheta,nbins[ivar]);
581 sprintf (splname,"fSplineSigmaphi[%d][%d]",ispline,ivar);
582 fSplineSigmaphi[ispline][ivar] = new TSpline3(splname,xspl,ysigmaphi,nbins[ivar]);
583 ispline++;
584 }
585 }
586
587 ivar = 1;
588 for (i=0; i<nbins[ivar]; i++) xspl[i] = x[i][ivar];
589 ispline=0;
590 for (Int_t ip=0; ip<fNbinp; ip++) {
591 for (Int_t iphi=0; iphi< fNbinphi; iphi++){
592 for (Int_t itheta=0; itheta< fNbintheta; itheta++){
593 // for efficiency and acceptance let's take the central value
594 // yeff[itheta] = fCurrentEntry[ip][itheta][iphi]->fEff[2];
595 // yacc[itheta] = fCurrentEntry[ip][itheta][iphi]->fAcc[2];
596 ysigmap[itheta] = fCurrentEntry[ip][itheta][iphi]->fSigmap;
597 ysigma1p[itheta] = fCurrentEntry[ip][itheta][iphi]->fSigma1p;
598 ysigmatheta[itheta] = fCurrentEntry[ip][itheta][iphi]->fSigmatheta;
599 ysigmaphi[itheta] = fCurrentEntry[ip][itheta][iphi]->fSigmaphi;
600 }
601 if (fPrintLevel>3) cout << " creating new spline " << splname << endl;
602 sprintf (splname,"fSplineEff[%d][%d]",ispline,ivar);
603 fSplineEff[ispline][ivar] = new TSpline3(splname,xspl,yeff, nbins[ivar]);
604 sprintf (splname,"fSplineAcc[%d][%d]",ispline,ivar);
605 fSplineAcc[ispline][ivar] = new TSpline3(splname,xspl,yacc, nbins[ivar]);
606 sprintf (splname,"fSplineSigmap[%d][%d]",ispline,ivar);
607 fSplineSigmap[ispline][ivar] = new TSpline3(splname,xspl,ysigmap,nbins[ivar]);
608 sprintf (splname,"fSplineSigma1p[%d][%d]",ispline,ivar);
609 fSplineSigma1p[ispline][ivar] = new TSpline3(splname,xspl,ysigma1p,nbins[ivar]);
610 sprintf (splname,"fSplineSigmatheta[%d][%d]",ispline,ivar);
611 fSplineSigmatheta[ispline][ivar] = new TSpline3(splname,xspl,ysigmatheta,nbins[ivar]);
612 sprintf (splname,"fSplineSigmaphi[%d][%d]",ispline,ivar);
613 fSplineSigmaphi[ispline][ivar] = new TSpline3(splname,xspl,ysigmaphi,nbins[ivar]);
614 ispline++;
615 }
616 }
617
618 ivar = 2;
619 for (i=0; i<nbins[ivar]; i++) xspl[i] = x[i][ivar];
620 ispline=0;
621 for (Int_t ip=0; ip<fNbinp; ip++) {
622 for (Int_t itheta=0; itheta< fNbintheta; itheta++){
623 for (Int_t iphi=0; iphi< fNbinphi; iphi++){
624 // for efficiency and acceptance let's take the central value
625 // yeff[iphi] = fCurrentEntry[ip][itheta][iphi]->fEff[2];
626 // yacc[iphi] = fCurrentEntry[ip][itheta][iphi]->fAcc[2];
627 ysigmap[iphi] = fCurrentEntry[ip][itheta][iphi]->fSigmap;
628 ysigma1p[iphi] = fCurrentEntry[ip][itheta][iphi]->fSigma1p;
629 ysigmatheta[iphi] = fCurrentEntry[ip][itheta][iphi]->fSigmatheta;
630 ysigmaphi[iphi] = fCurrentEntry[ip][itheta][iphi]->fSigmaphi;
631 }
632 if (fPrintLevel>3) cout << " creating new spline " << splname << endl;
633 sprintf (splname,"fSplineEff[%d][%d]",ispline,ivar);
634 fSplineEff[ispline][ivar] = new TSpline3(splname,xspl,yeff, nbins[ivar]);
635 sprintf (splname,"fSplineAcc[%d][%d]",ispline,ivar);
636 fSplineAcc[ispline][ivar] = new TSpline3(splname,xspl,yacc, nbins[ivar]);
637 sprintf (splname,"fSplineSigmap[%d][%d]",ispline,ivar);
638 fSplineSigmap[ispline][ivar] = new TSpline3(splname,xspl,ysigmap,nbins[ivar]);
639 sprintf (splname,"fSplineSigma1p[%d][%d]",ispline,ivar);
640 fSplineSigma1p[ispline][ivar] = new TSpline3(splname,xspl,ysigma1p,nbins[ivar]);
641 sprintf (splname,"fSplineSigmatheta[%d][%d]",ispline,ivar);
642 fSplineSigmatheta[ispline][ivar] = new TSpline3(splname,xspl,ysigmatheta,nbins[ivar]);
643 sprintf (splname,"fSplineSigmaphi[%d][%d]",ispline,ivar);
644 fSplineSigmaphi[ispline][ivar] = new TSpline3(splname,xspl,ysigmaphi,nbins[ivar]);
645 ispline++;
646 }
647 }
648 printf ("...done\n");
649}
650
651void AliMUONFastTracking::SmearMuon(Float_t pgen, Float_t thetagen, Float_t phigen,
652 Int_t charge, Float_t &psmear, Float_t &thetasmear,
653 Float_t &phismear, Float_t &eff, Float_t &acc){
654
655 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
656 // IMPORTANT NOTICE TO THE USER
657 //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
658 // THIS METHOD HAS BEEN REPLACED BY AliFastMuonTrackingEff::Evaluate()
659 // AND WILL BE DELETED SOON
660 // DO NOT USE THIS METHOD
661 //
662
663 printf ("AliMUONFastTracking::SmearMuon() THIS METHOD IS OBSOLETE ");
664 printf ("PLEASE REFER TO AliFastMuonTrackingEff::Evaluate()\n");
665 // angles are in degrees
666
667 Double_t meanp = MeanP (pgen, thetagen, phigen, charge);
668 Double_t sigmap = SigmaP (pgen, thetagen, phigen, charge);
669 Double_t sigma1p = Sigma1P(pgen, thetagen, phigen, charge);
670 Double_t normg2 = NormG2 (pgen, thetagen, phigen, charge);
671 Double_t meang2 = MeanG2 (pgen, thetagen, phigen, charge);
672 Double_t sigmag2 = SigmaG2(pgen, thetagen, phigen, charge);
673
674 printf ("fBkg = %g normg2 (100,5,0,1) = %g \n",fBkg,NormG2(100,5,0,1));
675 printf ("fBkg = %g meang2 (100,5,0,1) = %g \n",fBkg,MeanG2(100,5,0,1));
676 printf ("fBkg = %g sigmag2 (100,5,0,1) = %g \n",fBkg,SigmaG2(100,5,0,1));
677 Int_t ip,itheta,iphi;
678 GetIpIthetaIphi(pgen, thetagen, phigen, charge, ip, itheta, iphi);
679 if (sigmap == 0) {
680 if (fPrintLevel>0) {
681 printf ("WARNING!!! sigmap=0: ");
682 printf ("ip= %d itheta = %d iphi = %d ", ip, itheta, iphi);
683 printf ("p= %f theta = %f phi = %f\n", pgen, thetagen, phigen);
684 }
685 }
686
687 if (fPrintLevel>1) printf ("setting parameters: meanp = %f sigmap = %f sigma1p = %f normg2 = %f meang2 = %f sigmag2 = %f \n",meanp,sigmap,sigma1p,normg2,meang2,sigmag2);
688 fFitp->SetParameters(meanp,sigmap,sigma1p,normg2,meang2,sigmag2);
689
690 Double_t meantheta = MeanTheta (pgen, thetagen, phigen, charge);
691 Double_t sigmatheta = SigmaTheta(pgen, thetagen, phigen, charge);
692 Double_t meanphi = MeanPhi (pgen, thetagen, phigen, charge);
693 Double_t sigmaphi = SigmaPhi (pgen, thetagen, phigen, charge);
694
695 // components different from ip=0 have the RMS bigger than mean
696 Float_t ptp[3] = { 1.219576,-0.354764,-0.690117 };
697 Float_t ptph[3] = { 0.977522, 0.016269, 0.023158 };
698 Float_t pphp[3] = { 1.303256,-0.464847,-0.869322 };
699 psmear = pgen + fFitp->GetRandom();
700 Float_t dp = psmear - pgen;
701 if (ip==0) sigmaphi *= pphp[0] + pphp[1] * dp + pphp[2] * dp*dp;
702 phismear = phigen + gRandom->Gaus(meanphi, sigmaphi);
703 Float_t dphi = phismear - phigen;
704
705 if (ip==0) sigmatheta *= ptp[0] + ptp[1] * dp + ptp[2] * dp*dp;
706 if (ip==0) sigmatheta *= ptph[0] + ptph[1] * dphi + ptph[2] * dphi*dphi;
707 thetasmear = thetagen + gRandom->Gaus(meantheta,sigmatheta);
708 eff = Efficiency(pgen, thetagen, phigen, charge);
709 acc = Acceptance(pgen, thetagen, phigen, charge);
710}
711
712void AliMUONFastTracking::SetBackground(Float_t bkg){
713 // linear interpolation of the parameters in the LUT between 2 values where
714 // the background has been actually calculated
715
716 if (bkg>2) printf ("WARNING: unsafe extrapolation!\n");
717 fBkg = bkg;
718
719 Float_t BKG[4] = {0, 0.5, 1, 2}; // bkg values for which LUT is calculated
720 Int_t ibkg;
721 for (ibkg=0; ibkg<4; ibkg++) if ( bkg < BKG[ibkg]) break;
722 if (ibkg == 4) ibkg--;
723 if (ibkg == 0) ibkg++;
724
725 Float_t x0 = BKG[ibkg-1];
726 Float_t x1 = BKG[ibkg];
727 Float_t x = (bkg - x0) / (x1 - x0);
728
729 Float_t y0, y1;
730
731 for (Int_t ip=0; ip< fNbinp; ip++){
732 for (Int_t itheta=0; itheta< fNbintheta; itheta++){
733 for (Int_t iphi=0; iphi< fNbinphi; iphi++){
734 fCurrentEntry[ip][itheta][iphi]->fP = fEntry[ip][itheta][iphi][ibkg]->fP;
735 fCurrentEntry[ip][itheta][iphi]->fTheta = fEntry[ip][itheta][iphi][ibkg]->fTheta;
736 fCurrentEntry[ip][itheta][iphi]->fPhi = fEntry[ip][itheta][iphi][ibkg]->fPhi;
737 fCurrentEntry[ip][itheta][iphi]->fChi2p = -1;
738 fCurrentEntry[ip][itheta][iphi]->fChi2theta = -1;
739 fCurrentEntry[ip][itheta][iphi]->fChi2phi = -1;
740
741 y0 = fEntry[ip][itheta][iphi][ibkg-1]->fMeanp;
742 y1 = fEntry[ip][itheta][iphi][ibkg]->fMeanp;
743 fCurrentEntry[ip][itheta][iphi] ->fMeanp = (y1 - y0) * x + y0;
744 y0 = fEntry[ip][itheta][iphi][ibkg-1]->fMeantheta;
745 y1 = fEntry[ip][itheta][iphi][ibkg]->fMeantheta;
746 fCurrentEntry[ip][itheta][iphi] ->fMeantheta = (y1 - y0) * x + y0;
747 y0 = fEntry[ip][itheta][iphi][ibkg-1]->fMeanphi;
748 y1 = fEntry[ip][itheta][iphi][ibkg]->fMeanphi;
749 fCurrentEntry[ip][itheta][iphi] ->fMeanphi = (y1 - y0) * x + y0;
750 y0 = fEntry[ip][itheta][iphi][ibkg-1]->fSigmap;
751 y1 = fEntry[ip][itheta][iphi][ibkg]->fSigmap;
752 fCurrentEntry[ip][itheta][iphi] ->fSigmap = (y1 - y0) * x + y0;
753 y0 = fEntry[ip][itheta][iphi][ibkg-1]->fSigmatheta;
754 y1 = fEntry[ip][itheta][iphi][ibkg]->fSigmatheta;
755 fCurrentEntry[ip][itheta][iphi] ->fSigmatheta = (y1 - y0) * x + y0;
756 y0 = fEntry[ip][itheta][iphi][ibkg-1]->fSigmaphi;
757 y1 = fEntry[ip][itheta][iphi][ibkg]->fSigmaphi;
758 fCurrentEntry[ip][itheta][iphi] ->fSigmaphi = (y1 - y0) * x + y0;
759 y0 = fEntry[ip][itheta][iphi][ibkg-1]->fSigma1p;
760 y1 = fEntry[ip][itheta][iphi][ibkg]->fSigma1p;
761 fCurrentEntry[ip][itheta][iphi] ->fSigma1p = (y1 - y0) * x + y0;
762 y0 = fEntry[ip][itheta][iphi][ibkg-1]->fNormG2;
763 y1 = fEntry[ip][itheta][iphi][ibkg]->fNormG2;
764 fCurrentEntry[ip][itheta][iphi] ->fNormG2 = (y1 - y0) * x + y0;
765 y0 = fEntry[ip][itheta][iphi][ibkg-1]->fMeanG2;
766 y1 = fEntry[ip][itheta][iphi][ibkg]->fMeanG2;
767 fCurrentEntry[ip][itheta][iphi] ->fMeanG2 = (y1 - y0) * x + y0;
768
769 y0 = fEntry[ip][itheta][iphi][ibkg-1]->fSigmaG2;
770 y1 = fEntry[ip][itheta][iphi][ibkg]->fSigmaG2;
771 fCurrentEntry[ip][itheta][iphi] ->fSigmaG2 = (y1 - y0) * x + y0;
772 for (Int_t i=0; i<kSplitP; i++) {
773 for (Int_t j=0; j<kSplitTheta; j++) {
774 fCurrentEntry[ip][itheta][iphi]->fAcc[i][j] = fEntry[ip][itheta][iphi][ibkg]->fAcc[i][j];
775 y0 = fEntry[ip][itheta][iphi][ibkg-1]->fEff[i][j];
776 y1 = fEntry[ip][itheta][iphi][ibkg]->fEff[i][j];
777 fCurrentEntry[ip][itheta][iphi]->fEff[i][j] = (y1 - y0) * x + y0;
778 }
779 }
780 }
781 }
782 }
783 SetSpline();
784}
785
786
787
788 // to guarantee a safe extrapolation for sigmag2 to 0<bkg<0.5, let's fit
789 // with a straight line sigmag2 vs bkg for bkg=0.5, 1 and 2, and put the
790 // sigma2(BKG=0) as the extrapolation of this fit
791
792
793
794
795
796
797
798