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