FASTSIM/AliMUONFastTracking.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /*
  17 $Log$
  18 Revision 1.2  2003/01/08 10:29:33  morsch
  19 Path to data file changed.
  20
  21 Revision 1.1  2003/01/06 10:13:33  morsch
  22 First commit.
  23
  24 */
  25
  26 #include "AliMUONFastTracking.h"
  27 #include "AliMUONFastTrackingEntry.h"
  28 #include <TMatrixD.h>
  29 #include <TSpline.h>
  30 #include <TFile.h>
  31 #include <TH1.h>
  32 #include <TH3.h>
  33 #include <TF1.h>
  34 #include <TRandom.h>
  35 #include <stdlib.h>
  36 #include <stdio.h>
  37 #include <string.h>
  38 #include <Riostream.h>
  39
  40 ClassImp(AliMUONFastTracking)
  41
  42 AliMUONFastTracking* AliMUONFastTracking::fgMUONFastTracking=NULL;
  43
  44 static Double_t FitP(Double_t *x, Double_t *par){
  45     Double_t dx = x[0] - par[0];
  46     Double_t dx2 = x[0] - par[4];
  47     Double_t sigma = par[1] * ( 1 + par[2] * dx);
  48     if (sigma == 0) {
  49
  50         return 0.;
  51     }
  52     Double_t fasymm = TMath::Exp(-0.5 * dx * dx / (sigma * sigma));
  53     Double_t sigma2 = par[1] * par[5];
  54     Double_t fgauss = TMath::Exp(-0.5 * dx2 * dx2 / (sigma2 * sigma2));
  55     Double_t value = fasymm + par[3] * fgauss;
  56     return value;
  57 }
  58
  59 AliMUONFastTracking* AliMUONFastTracking::Instance()
  60 {
  61 // Set random number generator
  62     if (fgMUONFastTracking) {
  63         return fgMUONFastTracking;
  64     } else {
  65         fgMUONFastTracking = new AliMUONFastTracking();
  66         return fgMUONFastTracking;
  67     }
  68 }
  69
  70 AliMUONFastTracking::AliMUONFastTracking()
  71 {
  72     fClusterFinder = kOld;
  73     fPrintLevel = 1;
  74     // read binning; temporarily put by hand
  75     Float_t pmin = 0, pmax = 200;
  76     Int_t   nbinp = 10;
  77     Float_t thetamin = 2, thetamax = 9;
  78     Int_t   nbintheta=10;
  79     Float_t phimin = -180, phimax =180;
  80     Int_t   nbinphi=10;
  81     //--------------------------------------
  82
  83     fNbinp = nbinp;
  84     fPmin  = pmin;
  85     fPmax  = pmax;
  86
  87     fNbintheta = nbintheta;
  88     fThetamin  = thetamin;
  89     fThetamax  = thetamax;
  90
  91     fNbinphi = nbinphi;
  92     fPhimin  = phimin;
  93     fPhimax  = phimax;
  94
  95     fDeltaP     = (fPmax-fPmin)/fNbinp;
  96     fDeltaTheta = (fThetamax-fThetamin)/fNbintheta;
  97     fDeltaPhi   = (fPhimax-fPhimin)/fNbinphi;
  98 }
  99
 100 void AliMUONFastTracking::Init(Float_t bkg)
 101 {
 102   //
 103   //  Initialization
 104   //
 105   for (Int_t ip=0; ip< fNbinp; ip++){
 106     for (Int_t itheta=0; itheta< fNbintheta; itheta++){
 107       for (Int_t iphi=0; iphi< fNbinphi; iphi++){
 108         fCurrentEntry[ip][itheta][iphi] = new AliMUONFastTrackingEntry;
 109         for (Int_t ibkg=0; ibkg<4; ibkg++){
 110           fEntry[ip][itheta][iphi][ibkg] = new AliMUONFastTrackingEntry;
 111         }
 112       }
 113     }
 114   }
 115
 116   char filename [100];
 117   if (fClusterFinder==kOld) sprintf (filename,"$(ALICE_ROOT)/FASTSIM/data/MUONtrackLUT.root");
 118   else sprintf (filename,"$(ALICE_ROOT)/FASTSIM/data/MUONtrackLUT.root");
 119
 120   TFile *file = new TFile(filename);
 121   ReadLUT(file);
 122   SetBackground(bkg);
 123   UseSpline(0);
 124 }
 125
 126
 127 void AliMUONFastTracking::ReadLUT(TFile* file)
 128 {
 129   TH3F *heff[5][3], *hacc[5][3], *hmeanp, *hsigmap, *hsigma1p, *hchi2p;
 130   TH3F *hnormg2, *hmeang2, *hsigmag2, *hmeantheta, *hsigmatheta, *hchi2theta;
 131   TH3F *hmeanphi, *hsigmaphi, *hchi2phi;
 132   char tag[40], tag2[40];
 133
 134   printf ("Reading parameters from LUT file %s...\n",file->GetName());
 135
 136   const Float_t bkg[4] = {0, 0.5, 1, 2};
 137   for (Int_t ibkg=0; ibkg<4; ibkg++) {
 138     sprintf (tag,"BKG%g",bkg[ibkg]);
 139     file->cd(tag);
 140     for (Int_t isplp = 0; isplp<kSplitP; isplp++) {
 141       for (Int_t ispltheta = 0; ispltheta<kSplitTheta; ispltheta++) {
 142         sprintf (tag2,"heff[%d][%d]",isplp,ispltheta);
 143         heff[isplp][ispltheta] = (TH3F*)gDirectory->Get(tag2);
 144         sprintf (tag2,"hacc[%d][%d]",isplp,ispltheta);
 145         hacc[isplp][ispltheta] = (TH3F*)gDirectory->Get(tag2);
 146       }
 147     }
 148     hmeanp      = (TH3F*)gDirectory->Get("hmeanp");
 149     hsigmap     = (TH3F*)gDirectory->Get("hsigmap");
 150     hsigma1p    = (TH3F*)gDirectory->Get("hsigma1p");
 151     hchi2p      = (TH3F*)gDirectory->Get("hchi2p");
 152     hnormg2     = (TH3F*)gDirectory->Get("hnormg2");
 153     hmeang2     = (TH3F*)gDirectory->Get("hmeang2");
 154     hsigmag2    = (TH3F*)gDirectory->Get("hsigmag2");
 155     hmeantheta  = (TH3F*)gDirectory->Get("hmeantheta");
 156     hsigmatheta = (TH3F*)gDirectory->Get("hsigmatheta");
 157     hchi2theta  = (TH3F*)gDirectory->Get("hchi2theta");
 158     hmeanphi    = (TH3F*)gDirectory->Get("hmeanphi");
 159     hsigmaphi   = (TH3F*)gDirectory->Get("hsigmaphi");
 160     hchi2phi    = (TH3F*)gDirectory->Get("hchi2phi");
 161
 162     for (Int_t ip=0; ip<fNbinp ;ip++) {
 163       for (Int_t itheta=0; itheta<fNbintheta ;itheta++) {
 164         for (Int_t iphi=0; iphi<fNbinphi ;iphi++) {
 165           Float_t p     = fPmin     + fDeltaP     * (ip     + 0.5);
 166           Float_t theta = fThetamin + fDeltaTheta * (itheta + 0.5);
 167           Float_t phi   = fPhimin   + fDeltaPhi   * (iphi   + 0.5);
 168
 169           fEntry[ip][itheta][iphi][ibkg]->fP          = p;
 170           fEntry[ip][itheta][iphi][ibkg]->fMeanp      =
 171             hmeanp->GetBinContent(ip+1,itheta+1,iphi+1);
 172           fEntry[ip][itheta][iphi][ibkg]->fSigmap     =
 173             TMath::Abs(hsigmap->GetBinContent(ip+1,itheta+1,iphi+1));
 174           fEntry[ip][itheta][iphi][ibkg]->fSigma1p    =
 175             hsigma1p->GetBinContent(ip+1,itheta+1,iphi+1);
 176           fEntry[ip][itheta][iphi][ibkg]->fChi2p      =
 177             hchi2p->GetBinContent(ip+1,itheta+1,iphi+1);
 178           fEntry[ip][itheta][iphi][ibkg]->fNormG2     =
 179             hnormg2->GetBinContent(ip+1,itheta+1,iphi+1);
 180           fEntry[ip][itheta][iphi][ibkg]->fMeanG2     =
 181             hmeang2->GetBinContent(ip+1,itheta+1,iphi+1);
 182           if (ibkg == 0)  fEntry[ip][itheta][iphi][ibkg]->fSigmaG2 = 9999;
 183           else fEntry[ip][itheta][iphi][ibkg]->fSigmaG2    =
 184             hsigmag2->GetBinContent(ip+1,itheta+1,iphi+1);
 185           fEntry[ip][itheta][iphi][ibkg]->fTheta      = theta;
 186           fEntry[ip][itheta][iphi][ibkg]->fMeantheta  =
 187             hmeantheta->GetBinContent(ip+1,itheta+1,iphi+1);
 188           fEntry[ip][itheta][iphi][ibkg]->fSigmatheta =
 189             TMath::Abs(hsigmatheta->GetBinContent(ip+1,itheta+1,iphi+1));
 190           fEntry[ip][itheta][iphi][ibkg]->fChi2theta  =
 191             hchi2theta->GetBinContent(ip+1,itheta+1,iphi+1);
 192           fEntry[ip][itheta][iphi][ibkg]->fPhi        = phi;
 193           fEntry[ip][itheta][iphi][ibkg]->fMeanphi    =
 194             hmeanphi->GetBinContent(ip+1,itheta+1,iphi+1);
 195           fEntry[ip][itheta][iphi][ibkg]->fSigmaphi   =
 196             TMath::Abs(hsigmaphi->GetBinContent(ip+1,itheta+1,iphi+1));
 197           fEntry[ip][itheta][iphi][ibkg]->fChi2phi    =
 198             hchi2phi->GetBinContent(ip+1,itheta+1,iphi+1);
 199           for (Int_t i=0; i<kSplitP; i++) {
 200             for (Int_t j=0; j<kSplitTheta; j++) {
 201               fEntry[ip][itheta][iphi][ibkg]->fAcc[i][j] =
 202                 hacc[i][j]->GetBinContent(ip+1,itheta+1,iphi+1);
 203               fEntry[ip][itheta][iphi][ibkg]->fEff[i][j] =
 204                 heff[i][j]->GetBinContent(ip+1,itheta+1,iphi+1);
 205             }
 206           }
 207         } // iphi
 208       }  // itheta
 209     }   // ip
 210   }    // ibkg
 211
 212   TGraph *graph = new TGraph(3);
 213   TF1 *f = new TF1("f","[0]+[1]*x");
 214
 215   for (Int_t ip=0; ip< fNbinp; ip++){
 216     for (Int_t itheta=0; itheta< fNbintheta; itheta++){
 217       for (Int_t iphi=0; iphi< fNbinphi; iphi++){
 218         graph->SetPoint(0,0.5,fEntry[ip][itheta][iphi][1]->fSigmaG2);
 219         graph->SetPoint(1,1,fEntry[ip][itheta][iphi][2]->fSigmaG2);
 220         graph->SetPoint(2,2,fEntry[ip][itheta][iphi][3]->fSigmaG2);
 221         graph->Fit("f","q");
 222         fEntry[ip][itheta][iphi][0]->fSigmaG2 = f->Eval(0);
 223       }
 224     }
 225   }
 226   f->Delete();
 227   graph->Delete();
 228   printf ("parameters read. \n");
 229 }
 230
 231 void AliMUONFastTracking::GetBinning(Int_t &nbinp, Float_t &pmin, Float_t &pmax,
 232                                      Int_t &nbintheta, Float_t &thetamin,
 233                                      Float_t &thetamax,
 234                                      Int_t &nbinphi, Float_t &phimin, Float_t &phimax)
 235 {
 236     nbinp = fNbinp;
 237     pmin  = fPmin;
 238     pmax  = fPmax;
 239     nbintheta = fNbintheta;
 240     thetamin  = fThetamin;
 241     thetamax  = fThetamax;
 242     nbinphi = fNbinphi;
 243     phimin  = fPhimin;
 244     phimax  = fPhimax;
 245 }
 246
 247
 248 void AliMUONFastTracking::GetIpIthetaIphi(Float_t p, Float_t theta, Float_t phi,
 249                                           Int_t charge, Int_t &ip, Int_t &itheta,
 250                                           Int_t &iphi)
 251 {
 252     if (charge < 0) phi = -phi;
 253     ip           = Int_t (( p - fPmin ) / fDeltaP);
 254     itheta       = Int_t (( theta - fThetamin ) / fDeltaTheta);
 255     iphi         = Int_t (( phi - fPhimin ) / fDeltaPhi);
 256
 257     if (ip< 0)  ip = 0;
 258     if (ip>= fNbinp) ip = fNbinp-1;
 259     if (itheta< 0) itheta = 0;
 260     if (itheta>= fNbintheta) itheta = fNbintheta-1;
 261
 262     if (iphi< 0) iphi = 0;
 263     if (iphi>= fNbinphi) iphi = fNbinphi-1;
 264 }
 265
 266 void AliMUONFastTracking::GetSplit(Int_t ip, Int_t itheta,
 267                                    Int_t &nSplitP, Int_t &nSplitTheta) {
 268   if (ip==0) nSplitP = 5;
 269   else nSplitP = 2;
 270   if (itheta==0) nSplitTheta = 3;
 271   else nSplitTheta = 1;
 272 }
 273
 274 Float_t AliMUONFastTracking::Efficiency(Float_t p,   Float_t theta,
 275                                         Float_t phi, Int_t charge){
 276   Int_t ip=0, itheta=0, iphi=0;
 277   GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
 278   Int_t nSplitP, nSplitTheta;
 279   GetSplit(ip,itheta,nSplitP,nSplitTheta);
 280
 281   Float_t dp = p - fPmin;
 282   Int_t ibinp  = Int_t(nSplitP*(dp - fDeltaP * Int_t(dp / fDeltaP))/fDeltaP);
 283   Float_t dtheta = theta - fThetamin;
 284   Int_t ibintheta  = Int_t(nSplitTheta*(dtheta - fDeltaTheta * Int_t(dtheta / fDeltaTheta))/fDeltaTheta);
 285   Float_t eff = fCurrentEntry[ip][itheta][iphi]->fEff[ibinp][ibintheta];
 286   return eff;
 287 }
 288
 289 Float_t AliMUONFastTracking::Acceptance(Float_t p,   Float_t theta,
 290                                         Float_t phi, Int_t charge){
 291   if (theta<fThetamin || theta>fThetamax) return 0;
 292
 293   Int_t ip=0, itheta=0, iphi=0;
 294   GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
 295   Int_t nSplitP, nSplitTheta;
 296   GetSplit(ip,itheta,nSplitP,nSplitTheta);
 297   // central value and corrections with spline
 298
 299   Float_t dp = p - fPmin;
 300   Int_t ibinp  = Int_t(nSplitP*(dp - fDeltaP * Int_t(dp / fDeltaP))/fDeltaP);
 301   Float_t dtheta = theta - fThetamin;
 302   Int_t ibintheta  = Int_t(nSplitTheta*(dtheta - fDeltaTheta * Int_t(dtheta / fDeltaTheta))/fDeltaTheta);
 303   Float_t acc = fCurrentEntry[ip][itheta][iphi]->fAcc[ibinp][ibintheta];
 304   return acc;
 305 }
 306
 307 Float_t AliMUONFastTracking::MeanP(Float_t p,   Float_t theta,
 308                             Float_t phi, Int_t charge)
 309 {
 310     Int_t ip=0, itheta=0, iphi=0;
 311     GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
 312     return fCurrentEntry[ip][itheta][iphi]->fMeanp;
 313 }
 314
 315 Float_t AliMUONFastTracking::SigmaP(Float_t p,   Float_t theta,
 316                                     Float_t phi, Int_t charge)
 317 {
 318     Int_t ip=0, itheta=0, iphi=0;
 319     Int_t index;
 320     GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
 321     // central value and corrections with spline
 322     Float_t sigmap = fCurrentEntry[ip][itheta][iphi]->fSigmap;
 323     if (!fSpline) return sigmap;
 324     // corrections vs p, theta, phi
 325     index = iphi + fNbinphi * itheta;
 326     Double_t xmin,ymin,xmax,ymax;
 327     Float_t frac1 = fSplineSigmap[index][0]->Eval(p)/sigmap;
 328
 329     if (p>fPmax-fDeltaP/2.) {
 330         Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->fSigmap;
 331         Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->fSigmap;
 332         Float_t s3 = fCurrentEntry[fNbinp-3][itheta][iphi]->fSigmap;
 333         Float_t p1  = fDeltaP * (fNbinp - 1 + 0.5) + fPmin;
 334         Float_t p2  = fDeltaP * (fNbinp - 2 + 0.5) + fPmin;
 335         Float_t p3  = fDeltaP * (fNbinp - 3 + 0.5) + fPmin;
 336         Float_t p12 = p1 * p1, p22 = p2 * p2, p32 = p3 * p3;
 337         Float_t d = p12*p2 + p1*p32 + p22*p3 - p32*p2 - p3*p12 - p22*p1;
 338         Float_t a = (s1*p2 + p1*s3 + s2*p3 - s3*p2 - p3*s1 - s2*p1) / d;
 339         Float_t b = (p12*s2 + s1*p32 + p22*s3 - p32*s2 - s3*p12 - p22*s1)/d;
 340         Float_t c = (p12*p2*s3 + p1*p32*s2 + p22*p3*s1
 341                      - p32*p2*s1 - p3*p12*s2 - p22*p1*s3) / d;
 342         Float_t sigma = a * p * p + b * p + c;
 343         frac1 = sigma/sigmap;
 344     }
 345     index = iphi + fNbinphi * ip;
 346     fSplineEff[index][1]->GetKnot(0,xmin,ymin);
 347     fSplineEff[index][1]->GetKnot(9,xmax,ymax);
 348     if (theta>xmax) theta = xmax;
 349     Float_t frac2 = fSplineSigmap[index][1]->Eval(theta)/sigmap;
 350     index = itheta + fNbintheta * ip;
 351     fSplineEff[index][2]->GetKnot(0,xmin,ymin);
 352     fSplineEff[index][2]->GetKnot(9,xmax,ymax);
 353     if (phi>xmax) phi = xmax;
 354     Float_t frac3 = fSplineSigmap[index][2]->Eval(phi)/sigmap;
 355     Float_t sigmatot = sigmap * frac1 * frac2 * frac3;
 356     if (sigmatot<0) sigmatot = sigmap;
 357     return sigmatot;
 358 }
 359
 360 Float_t AliMUONFastTracking::Sigma1P(Float_t p,   Float_t theta,
 361                             Float_t phi, Int_t charge)
 362 {
 363     Int_t ip=0, itheta=0, iphi=0;
 364     GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
 365     if (p>fPmax) {
 366         // linear extrapolation of sigmap for p out of range
 367         Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->fSigma1p;
 368         Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->fSigma1p;
 369         Float_t p1  = fDeltaP * (fNbinp - 1 + 0.5) + fPmin;
 370         Float_t p2  = fDeltaP * (fNbinp - 2 + 0.5) + fPmin;
 371         Float_t sigma = 1./(p1-p2) * ( (s1-s2)*p + (s2-s1)*p1 + s1*(p1-p2) );
 372         return sigma;
 373     }
 374     else return fCurrentEntry[ip][itheta][iphi]->fSigma1p;
 375 }
 376
 377 Float_t AliMUONFastTracking::NormG2(Float_t p,   Float_t theta,
 378                                     Float_t phi, Int_t charge)
 379 {
 380     Int_t ip=0, itheta=0, iphi=0;
 381     GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
 382     if (p>fPmax) {
 383         // linear extrapolation of sigmap for p out of range
 384         Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->fNormG2;
 385         Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->fNormG2;
 386         Float_t p1  = fDeltaP * (fNbinp - 1 + 0.5) + fPmin;
 387         Float_t p2  = fDeltaP * (fNbinp - 2 + 0.5) + fPmin;
 388         Float_t norm = 1./(p1-p2) * ( (s1-s2)*p + (s2-s1)*p1 + s1*(p1-p2) );
 389         return norm;
 390     }
 391     else return fCurrentEntry[ip][itheta][iphi]->fNormG2;
 392 }
 393
 394 Float_t AliMUONFastTracking::MeanG2(Float_t p,   Float_t theta,
 395                                     Float_t phi, Int_t charge)
 396 {
 397     Int_t ip=0, itheta=0, iphi=0;
 398     GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
 399     if (p>fPmax) {
 400         // linear extrapolation of sigmap for p out of range
 401         Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->fMeanG2;
 402         Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->fMeanG2;
 403         Float_t p1  = fDeltaP * (fNbinp - 1 + 0.5) + fPmin;
 404         Float_t p2  = fDeltaP * (fNbinp - 2 + 0.5) + fPmin;
 405         Float_t norm = 1./(p1-p2) * ( (s1-s2)*p + (s2-s1)*p1 + s1*(p1-p2) );
 406         return norm;
 407     }
 408     else return fCurrentEntry[ip][itheta][iphi]->fMeanG2;
 409 }
 410
 411 Float_t AliMUONFastTracking::SigmaG2(Float_t p,   Float_t theta,
 412                                      Float_t phi, Int_t charge)
 413 {
 414     Int_t ip=0, itheta=0, iphi=0;
 415     GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
 416     if (p>fPmax) {
 417         // linear extrapolation of sigmap for p out of range
 418         Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->fSigmaG2;
 419         Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->fSigmaG2;
 420         Float_t p1  = fDeltaP * (fNbinp - 1 + 0.5) + fPmin;
 421         Float_t p2  = fDeltaP * (fNbinp - 2 + 0.5) + fPmin;
 422         Float_t sigma = 1./(p1-p2) * ( (s1-s2)*p + (s2-s1)*p1 + s1*(p1-p2) );
 423         return sigma;
 424     }
 425     else return fCurrentEntry[ip][itheta][iphi]->fSigmaG2;
 426 }
 427
 428
 429 Float_t AliMUONFastTracking::MeanTheta(Float_t p,   Float_t theta,
 430                                        Float_t phi, Int_t charge)
 431 {
 432     Int_t ip=0, itheta=0, iphi=0;
 433     GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
 434     return fCurrentEntry[ip][itheta][iphi]->fMeantheta;
 435 }
 436
 437 Float_t AliMUONFastTracking::SigmaTheta(Float_t p,   Float_t theta,
 438                             Float_t phi, Int_t charge){
 439   Int_t ip=0, itheta=0, iphi=0;
 440   Int_t index;
 441   GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
 442   // central value and corrections with spline
 443   Float_t sigmatheta = fCurrentEntry[ip][itheta][iphi]->fSigmatheta;
 444   if (!fSpline) return sigmatheta;
 445   // corrections vs p, theta, phi
 446   index = iphi + fNbinphi * itheta;
 447   Double_t xmin,ymin,xmax,ymax;
 448   Float_t frac1 = fSplineSigmatheta[index][0]->Eval(p)/sigmatheta;
 449   if (p>fPmax-fDeltaP/2.) {
 450     // linear extrapolation of sigmap for p out of range
 451     Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->fSigmatheta;
 452     Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->fSigmatheta;
 453     Float_t p1  = fDeltaP * (fNbinp - 1 + 0.5) + fPmin;
 454     Float_t p2  = fDeltaP * (fNbinp - 2 + 0.5) + fPmin;
 455     Float_t sigma = 1./(p1-p2) * ( (s1-s2)*p + (s2-s1)*p1 + s1*(p1-p2) );
 456     frac1=sigma/sigmatheta;
 457   }
 458   index = iphi + fNbinphi * ip;
 459   fSplineEff[index][1]->GetKnot(0,xmin,ymin);
 460   fSplineEff[index][1]->GetKnot(9,xmax,ymax);
 461   if (theta>xmax) theta = xmax;
 462   Float_t frac2 = fSplineSigmatheta[index][1]->Eval(theta)/sigmatheta;
 463   index = itheta + fNbintheta * ip;
 464   fSplineEff[index][2]->GetKnot(0,xmin,ymin);
 465   fSplineEff[index][2]->GetKnot(9,xmax,ymax);
 466   if (phi>xmax) phi = xmax;
 467   Float_t frac3 = fSplineSigmatheta[index][2]->Eval(phi)/sigmatheta;
 468   return sigmatheta * frac1 * frac2 * frac3;
 469 }
 470
 471
 472 Float_t AliMUONFastTracking::MeanPhi(Float_t p,   Float_t theta,
 473                             Float_t phi, Int_t charge){
 474   Int_t ip=0, itheta=0, iphi=0;
 475   GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
 476   return fCurrentEntry[ip][itheta][iphi]->fMeanphi;
 477 }
 478
 479 Float_t AliMUONFastTracking::SigmaPhi(Float_t p,   Float_t theta,
 480                             Float_t phi, Int_t charge){
 481   Int_t ip=0, itheta=0, iphi=0;
 482   Int_t index;
 483   GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
 484   // central value and corrections with spline
 485   Float_t sigmaphi = fCurrentEntry[ip][itheta][iphi]->fSigmaphi;
 486   if (!fSpline) return sigmaphi;
 487   // corrections vs p, theta, phi
 488   index = iphi + fNbinphi * itheta;
 489   Float_t frac1 = fSplineSigmaphi[index][0]->Eval(p)/sigmaphi;
 490   Double_t xmin,ymin,xmax,ymax;
 491   if (p>fPmax-fDeltaP/2.) {
 492     Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->fSigmaphi;
 493     Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->fSigmaphi;
 494     Float_t p1  = fDeltaP * (fNbinp - 1 + 0.5) + fPmin;
 495     Float_t p2  = fDeltaP * (fNbinp - 2 + 0.5) + fPmin;
 496     Float_t sigma = 1./(p1-p2) * ( (s1-s2)*p + (s2-s1)*p1 + s1*(p1-p2) );
 497     frac1 = sigma/sigmaphi;
 498   }
 499
 500   index = iphi + fNbinphi * ip;
 501   fSplineEff[index][1]->GetKnot(0,xmin,ymin);
 502   fSplineEff[index][1]->GetKnot(9,xmax,ymax);
 503   if (theta>xmax) theta = xmax;
 504   Float_t frac2 = fSplineSigmaphi[index][1]->Eval(theta)/sigmaphi;
 505   index = itheta + fNbintheta * ip;
 506   fSplineEff[index][2]->GetKnot(0,xmin,ymin);
 507   fSplineEff[index][2]->GetKnot(9,xmax,ymax);
 508   if (phi>xmax) phi = xmax;
 509   Float_t frac3 = fSplineSigmaphi[index][2]->Eval(phi)/sigmaphi;
 510   return sigmaphi * frac1 * frac2 * frac3;
 511 }
 512
 513 void AliMUONFastTracking::SetSpline(){
 514   printf ("Setting spline functions...");
 515   char splname[40];
 516   Double_t x[20][3];
 517   Double_t x2[50][3];
 518   Int_t nbins[3] = {fNbinp, fNbintheta, fNbinphi};
 519   Double_t xspl[20],yeff[50],ysigmap[20],ysigma1p[20];
 520   Double_t yacc[50], ysigmatheta[20],ysigmaphi[20];
 521   Double_t xsp2[50];
 522   // let's calculate the x axis for p, theta, phi
 523
 524   Int_t i, ispline, ivar;
 525   for (i=0; i< fNbinp; i++) x[i][0] = fPmin + fDeltaP * (i + 0.5);
 526   for (i=0; i< fNbintheta; i++) x[i][1] = fThetamin + fDeltaTheta * (i + 0.5);
 527   for (i=0; i< fNbinphi; i++) x[i][2] = fPhimin + fDeltaPhi * (i + 0.5);
 528
 529   for (i=0; i< 5 * fNbinp; i++) x2[i][0] = fPmin + fDeltaP * (i + 0.5)/5.;
 530   for (i=0; i< 5 * fNbintheta; i++) x2[i][1] = fThetamin + fDeltaTheta * (i + 0.5)/5.;
 531   for (i=0; i< 5 * fNbinphi; i++) x2[i][2] = fPhimin + fDeltaPhi * (i + 0.5)/5.;
 532
 533   // splines in p
 534   ivar = 0;
 535   for (i=0; i<nbins[ivar]; i++) xspl[i] = x[i][ivar];
 536   for (i=0; i<5 * nbins[ivar]; i++) xsp2[i] = x2[i][ivar];
 537   ispline=0;
 538   for (Int_t itheta=0; itheta< fNbintheta; itheta++){
 539     for (Int_t iphi=0; iphi< fNbinphi; iphi++){
 540       for (Int_t ip=0; ip<fNbinp; ip++) {
 541         //      for (Int_t i=0; i<5; i++) {
 542         //        yeff[5 * ip + i]  = fCurrentEntry[ip][itheta][iphi]->fEff[i];
 543         //        yacc[5 * ip + i]  = fCurrentEntry[ip][itheta][iphi]->fAcc[i];
 544         //      }
 545         ysigmap[ip]     = fCurrentEntry[ip][itheta][iphi]->fSigmap;
 546         ysigma1p[ip]    = fCurrentEntry[ip][itheta][iphi]->fSigma1p;
 547         ysigmatheta[ip] = fCurrentEntry[ip][itheta][iphi]->fSigmatheta;
 548         ysigmaphi[ip]   = fCurrentEntry[ip][itheta][iphi]->fSigmaphi;
 549       }
 550       if (fPrintLevel>3) cout << " creating new spline " << splname << endl;
 551       sprintf (splname,"fSplineEff[%d][%d]",ispline,ivar);
 552       fSplineEff[ispline][ivar] = new TSpline3(splname,xsp2,yeff,5 * nbins[ivar]);
 553       sprintf (splname,"fSplineAcc[%d][%d]",ispline,ivar);
 554       fSplineAcc[ispline][ivar] = new TSpline3(splname,xsp2,yacc,5 * nbins[ivar]);
 555       sprintf (splname,"fSplineSigmap[%d][%d]",ispline,ivar);
 556       fSplineSigmap[ispline][ivar] = new TSpline3(splname,xspl,ysigmap,nbins[ivar]);
 557       sprintf (splname,"fSplineSigma1p[%d][%d]",ispline,ivar);
 558       fSplineSigma1p[ispline][ivar] = new TSpline3(splname,xspl,ysigma1p,nbins[ivar]);
 559       sprintf (splname,"fSplineSigmatheta[%d][%d]",ispline,ivar);
 560       fSplineSigmatheta[ispline][ivar] = new TSpline3(splname,xspl,ysigmatheta,nbins[ivar]);
 561       sprintf (splname,"fSplineSigmaphi[%d][%d]",ispline,ivar);
 562       fSplineSigmaphi[ispline][ivar] = new TSpline3(splname,xspl,ysigmaphi,nbins[ivar]);
 563       ispline++;
 564     }
 565   }
 566
 567   ivar = 1;
 568   for (i=0; i<nbins[ivar]; i++) xspl[i] = x[i][ivar];
 569   ispline=0;
 570   for (Int_t ip=0; ip<fNbinp; ip++) {
 571     for (Int_t iphi=0; iphi< fNbinphi; iphi++){
 572       for (Int_t itheta=0; itheta< fNbintheta; itheta++){
 573         // for efficiency and acceptance let's take the central value
 574         //      yeff[itheta]        = fCurrentEntry[ip][itheta][iphi]->fEff[2];
 575         //      yacc[itheta]        = fCurrentEntry[ip][itheta][iphi]->fAcc[2];
 576         ysigmap[itheta]     = fCurrentEntry[ip][itheta][iphi]->fSigmap;
 577         ysigma1p[itheta]    = fCurrentEntry[ip][itheta][iphi]->fSigma1p;
 578         ysigmatheta[itheta] = fCurrentEntry[ip][itheta][iphi]->fSigmatheta;
 579         ysigmaphi[itheta]   = fCurrentEntry[ip][itheta][iphi]->fSigmaphi;
 580       }
 581       if (fPrintLevel>3) cout << " creating new spline " << splname << endl;
 582       sprintf (splname,"fSplineEff[%d][%d]",ispline,ivar);
 583       fSplineEff[ispline][ivar] = new TSpline3(splname,xspl,yeff, nbins[ivar]);
 584       sprintf (splname,"fSplineAcc[%d][%d]",ispline,ivar);
 585       fSplineAcc[ispline][ivar] = new TSpline3(splname,xspl,yacc, nbins[ivar]);
 586       sprintf (splname,"fSplineSigmap[%d][%d]",ispline,ivar);
 587       fSplineSigmap[ispline][ivar] = new TSpline3(splname,xspl,ysigmap,nbins[ivar]);
 588       sprintf (splname,"fSplineSigma1p[%d][%d]",ispline,ivar);
 589       fSplineSigma1p[ispline][ivar] = new TSpline3(splname,xspl,ysigma1p,nbins[ivar]);
 590       sprintf (splname,"fSplineSigmatheta[%d][%d]",ispline,ivar);
 591       fSplineSigmatheta[ispline][ivar] = new TSpline3(splname,xspl,ysigmatheta,nbins[ivar]);
 592       sprintf (splname,"fSplineSigmaphi[%d][%d]",ispline,ivar);
 593       fSplineSigmaphi[ispline][ivar] = new TSpline3(splname,xspl,ysigmaphi,nbins[ivar]);
 594       ispline++;
 595     }
 596   }
 597
 598   ivar = 2;
 599   for (i=0; i<nbins[ivar]; i++) xspl[i] = x[i][ivar];
 600   ispline=0;
 601   for (Int_t ip=0; ip<fNbinp; ip++) {
 602     for (Int_t itheta=0; itheta< fNbintheta; itheta++){
 603       for (Int_t iphi=0; iphi< fNbinphi; iphi++){
 604         // for efficiency and acceptance let's take the central value
 605         //      yeff[iphi]        = fCurrentEntry[ip][itheta][iphi]->fEff[2];
 606         //      yacc[iphi]        = fCurrentEntry[ip][itheta][iphi]->fAcc[2];
 607         ysigmap[iphi]     = fCurrentEntry[ip][itheta][iphi]->fSigmap;
 608         ysigma1p[iphi]    = fCurrentEntry[ip][itheta][iphi]->fSigma1p;
 609         ysigmatheta[iphi] = fCurrentEntry[ip][itheta][iphi]->fSigmatheta;
 610         ysigmaphi[iphi]   = fCurrentEntry[ip][itheta][iphi]->fSigmaphi;
 611       }
 612       if (fPrintLevel>3) cout << " creating new spline " << splname << endl;
 613       sprintf (splname,"fSplineEff[%d][%d]",ispline,ivar);
 614       fSplineEff[ispline][ivar] = new TSpline3(splname,xspl,yeff, nbins[ivar]);
 615       sprintf (splname,"fSplineAcc[%d][%d]",ispline,ivar);
 616       fSplineAcc[ispline][ivar] = new TSpline3(splname,xspl,yacc, nbins[ivar]);
 617       sprintf (splname,"fSplineSigmap[%d][%d]",ispline,ivar);
 618       fSplineSigmap[ispline][ivar] = new TSpline3(splname,xspl,ysigmap,nbins[ivar]);
 619       sprintf (splname,"fSplineSigma1p[%d][%d]",ispline,ivar);
 620       fSplineSigma1p[ispline][ivar] = new TSpline3(splname,xspl,ysigma1p,nbins[ivar]);
 621       sprintf (splname,"fSplineSigmatheta[%d][%d]",ispline,ivar);
 622       fSplineSigmatheta[ispline][ivar] = new TSpline3(splname,xspl,ysigmatheta,nbins[ivar]);
 623       sprintf (splname,"fSplineSigmaphi[%d][%d]",ispline,ivar);
 624       fSplineSigmaphi[ispline][ivar] = new TSpline3(splname,xspl,ysigmaphi,nbins[ivar]);
 625       ispline++;
 626     }
 627   }
 628   printf ("...done\n");
 629 }
 630
 631 void AliMUONFastTracking::SetBackground(Float_t bkg){
 632   // linear interpolation of the parameters in the LUT between 2 values where
 633   // the background has been actually calculated
 634
 635   if (bkg>2) printf ("WARNING: unsafe extrapolation!\n");
 636   fBkg = bkg;
 637
 638   Float_t BKG[4] = {0, 0.5, 1, 2}; // bkg values for which LUT is calculated
 639   Int_t ibkg;
 640   for (ibkg=0; ibkg<4; ibkg++) if ( bkg < BKG[ibkg]) break;
 641   if (ibkg == 4) ibkg--;
 642   if (ibkg == 0) ibkg++;
 643
 644   Float_t x0 = BKG[ibkg-1];
 645   Float_t x1 = BKG[ibkg];
 646   Float_t x = (bkg - x0) / (x1 - x0);
 647
 648   Float_t y0, y1;
 649
 650   for (Int_t ip=0; ip< fNbinp; ip++){
 651     for (Int_t itheta=0; itheta< fNbintheta; itheta++){
 652       for (Int_t iphi=0; iphi< fNbinphi; iphi++){
 653         fCurrentEntry[ip][itheta][iphi]->fP = fEntry[ip][itheta][iphi][ibkg]->fP;
 654         fCurrentEntry[ip][itheta][iphi]->fTheta = fEntry[ip][itheta][iphi][ibkg]->fTheta;
 655         fCurrentEntry[ip][itheta][iphi]->fPhi = fEntry[ip][itheta][iphi][ibkg]->fPhi;
 656         fCurrentEntry[ip][itheta][iphi]->fChi2p = -1;
 657         fCurrentEntry[ip][itheta][iphi]->fChi2theta = -1;
 658         fCurrentEntry[ip][itheta][iphi]->fChi2phi = -1;
 659
 660         y0 = fEntry[ip][itheta][iphi][ibkg-1]->fMeanp;
 661         y1 =   fEntry[ip][itheta][iphi][ibkg]->fMeanp;
 662         fCurrentEntry[ip][itheta][iphi]      ->fMeanp = (y1 - y0) * x + y0;
 663         y0 = fEntry[ip][itheta][iphi][ibkg-1]->fMeantheta;
 664         y1 =   fEntry[ip][itheta][iphi][ibkg]->fMeantheta;
 665         fCurrentEntry[ip][itheta][iphi]      ->fMeantheta = (y1 - y0) * x + y0;
 666         y0 = fEntry[ip][itheta][iphi][ibkg-1]->fMeanphi;
 667         y1 =   fEntry[ip][itheta][iphi][ibkg]->fMeanphi;
 668         fCurrentEntry[ip][itheta][iphi]      ->fMeanphi = (y1 - y0) * x + y0;
 669         y0 = fEntry[ip][itheta][iphi][ibkg-1]->fSigmap;
 670         y1 =   fEntry[ip][itheta][iphi][ibkg]->fSigmap;
 671         fCurrentEntry[ip][itheta][iphi]      ->fSigmap = (y1 - y0) * x + y0;
 672         y0 = fEntry[ip][itheta][iphi][ibkg-1]->fSigmatheta;
 673         y1 =   fEntry[ip][itheta][iphi][ibkg]->fSigmatheta;
 674         fCurrentEntry[ip][itheta][iphi]      ->fSigmatheta = (y1 - y0) * x + y0;
 675         y0 = fEntry[ip][itheta][iphi][ibkg-1]->fSigmaphi;
 676         y1 =   fEntry[ip][itheta][iphi][ibkg]->fSigmaphi;
 677         fCurrentEntry[ip][itheta][iphi]      ->fSigmaphi = (y1 - y0) * x + y0;
 678         y0 = fEntry[ip][itheta][iphi][ibkg-1]->fSigma1p;
 679         y1 =   fEntry[ip][itheta][iphi][ibkg]->fSigma1p;
 680         fCurrentEntry[ip][itheta][iphi]      ->fSigma1p = (y1 - y0) * x + y0;
 681         y0 = fEntry[ip][itheta][iphi][ibkg-1]->fNormG2;
 682         y1 =   fEntry[ip][itheta][iphi][ibkg]->fNormG2;
 683         fCurrentEntry[ip][itheta][iphi]      ->fNormG2 = (y1 - y0) * x + y0;
 684         y0 = fEntry[ip][itheta][iphi][ibkg-1]->fMeanG2;
 685         y1 =   fEntry[ip][itheta][iphi][ibkg]->fMeanG2;
 686         fCurrentEntry[ip][itheta][iphi]      ->fMeanG2 = (y1 - y0) * x + y0;
 687
 688         y0 = fEntry[ip][itheta][iphi][ibkg-1]->fSigmaG2;
 689         y1 =   fEntry[ip][itheta][iphi][ibkg]->fSigmaG2;
 690         fCurrentEntry[ip][itheta][iphi]      ->fSigmaG2 = (y1 - y0) * x + y0;
 691         for (Int_t i=0; i<kSplitP; i++) {
 692           for (Int_t j=0; j<kSplitTheta; j++) {
 693             fCurrentEntry[ip][itheta][iphi]->fAcc[i][j] = fEntry[ip][itheta][iphi][ibkg]->fAcc[i][j];
 694             y0 = fEntry[ip][itheta][iphi][ibkg-1]->fEff[i][j];
 695             y1 =   fEntry[ip][itheta][iphi][ibkg]->fEff[i][j];
 696             fCurrentEntry[ip][itheta][iphi]->fEff[i][j] = (y1 - y0) * x + y0;
 697           }
 698         }
 699       }
 700     }
 701   }
 702   SetSpline();
 703 }
 704
 705 TF1* AliMUONFastTracking::GetFitP(Int_t ip,Int_t itheta,Int_t iphi) {
 706   if (!fFitp[ip][itheta][iphi]) {
 707     fFitp[ip][itheta][iphi] = new TF1("fit1",FitP,-20.,20.,6);
 708     fFitp[ip][itheta][iphi]->SetNpx(500);
 709     fFitp[ip][itheta][iphi]->SetParameters(0,0,0,0,0,0);
 710   }
 711   return fFitp[ip][itheta][iphi];
 712 }
 713
 714   // to guarantee a safe extrapolation for sigmag2 to 0<bkg<0.5, let's fit
 715   // with a straight line sigmag2 vs bkg for bkg=0.5, 1 and 2, and put the
 716   // sigma2(BKG=0) as the extrapolation of this fit
 717
 718
 719
 720
 721
 722
 723
 724