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