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