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