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