Some problems with too big fluctuations corrected. (A. de Falco)
[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.2  2003/01/08 10:29:33  morsch
19 Path to data file changed.
20
21 Revision 1.1  2003/01/06 10:13:33  morsch
22 First commit.
23
24 */
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
40 ClassImp(AliMUONFastTracking)
41
42 AliMUONFastTracking* AliMUONFastTracking::fgMUONFastTracking=NULL;
43
44 static 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));
55     Double_t value = fasymm + par[3] * fgauss; 
56     return value;
57
58
59 AliMUONFastTracking* 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
70 AliMUONFastTracking::AliMUONFastTracking() 
71 {
72     fClusterFinder = kOld; 
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
100 void 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]; 
117   if (fClusterFinder==kOld) sprintf (filename,"$(ALICE_ROOT)/FASTSIM/data/MUONtrackLUT.root"); 
118   else sprintf (filename,"$(ALICE_ROOT)/FASTSIM/data/MUONtrackLUT.root"); 
119
120   TFile *file = new TFile(filename); 
121   ReadLUT(file);
122   SetBackground(bkg);
123   UseSpline(0);
124 }
125
126
127 void 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   
134   printf ("Reading parameters from LUT file %s...\n",file->GetName());
135
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     
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
231 void 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
248 void 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     
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;
261     
262     if (iphi< 0) iphi = 0;
263     if (iphi>= fNbinphi) iphi = fNbinphi-1;
264 }
265
266 void 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
274 Float_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
289 Float_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
307 Float_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
315 Float_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
360 Float_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
377 Float_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
394 Float_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
411 Float_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
429 Float_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
437 Float_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
472 Float_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
479 Float_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
513 void 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   
631 void 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
705 TF1* 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 }
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