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