Saver initialisation of fFitp array. (Lenaic COUEDEL)
[u/mrichter/AliRoot.git] / FASTSIM / AliMUONFastTracking.cxx
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$
18 Revision 1.5  2003/08/05 16:14:20  morsch
19 Some problems with too big fluctuations corrected. (A. de Falco)
20
21 Revision 1.2  2003/01/08 10:29:33  morsch
22 Path to data file changed.
23
24 Revision 1.1  2003/01/06 10:13:33  morsch
25 First 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
43 ClassImp(AliMUONFastTracking)
44
45 AliMUONFastTracking* AliMUONFastTracking::fgMUONFastTracking=NULL;
46
47 static 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
62 AliMUONFastTracking* 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
73 AliMUONFastTracking::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
111 void 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
138 void 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
242 void 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
259 void 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
278 void 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
286 Float_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
301 Float_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
319 Float_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
327 Float_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
372 Float_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
389 Float_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
406 Float_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
423 Float_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
441 Float_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
449 Float_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
484 Float_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
491 Float_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
525 void 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   
643 void 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
717 TF1* 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