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