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