]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - FASTSIM/AliMUONFastTracking.cxx
Some small corrections to avoid infinite loops at high momenta.
[u/mrichter/AliRoot.git] / FASTSIM / AliMUONFastTracking.cxx
index ce9f021ac45ef294151079a2027c448c71b39ac1..3f979a64c44e7dd2f84dcce22aa3380969358fc1 100644 (file)
 
 /*
 $Log$
+Revision 1.9  2004/02/03 16:45:20  morsch
+Unique name for functions (TF1).
+
+Revision 1.8  2004/01/27 18:02:38  hristov
+Removing some warning (Sun)
+
+Revision 1.7  2003/11/13 14:21:57  morsch
+Coding Rule violation corrections.
+
+Revision 1.6  2003/08/12 15:16:25  morsch
+Saver initialisation of fFitp  array. (Lenaic COUEDEL)
+
+Revision 1.5  2003/08/05 16:14:20  morsch
+Some problems with too big fluctuations corrected. (A. de Falco)
+
+Revision 1.2  2003/01/08 10:29:33  morsch
+Path to data file changed.
+
 Revision 1.1  2003/01/06 10:13:33  morsch
 First commit.
 
 */
 
+//-------------------------------------------------------------------------
+//        Class AliMUONFastTracking 
+//
+//  Manager for the fast simulation of tracking in the muon spectrometer
+//  This class reads the lookup tables containing the parameterization 
+//  of the deltap, deltatheta, deltaphi for different background levels
+//  and provides the related smeared parameters. 
+//  Used by AliFastMuonTrackingEff, AliFastMuonTrackingAcc, 
+//  AliFastMuonTrackingRes. 
+//-------------------------------------------------------------------------
+
 #include "AliMUONFastTracking.h"
 #include "AliMUONFastTrackingEntry.h"
-#include <TMatrixD.h>
 #include <TSpline.h>
 #include <TFile.h>
-#include <TH1.h>
 #include <TH3.h>
 #include <TF1.h>
 #include <TRandom.h>
@@ -36,6 +63,7 @@ First commit.
 
 ClassImp(AliMUONFastTracking)
 
+
 AliMUONFastTracking* AliMUONFastTracking::fgMUONFastTracking=NULL;
 
 static Double_t FitP(Double_t *x, Double_t *par){
@@ -49,9 +77,17 @@ static Double_t FitP(Double_t *x, Double_t *par){
     Double_t fasymm = TMath::Exp(-0.5 * dx * dx / (sigma * sigma));
     Double_t sigma2 = par[1] * par[5];
     Double_t fgauss = TMath::Exp(-0.5 * dx2 * dx2 / (sigma2 * sigma2));
-    return TMath::Abs(fasymm + par[3] * fgauss);
+    Double_t value = fasymm + par[3] * fgauss; 
+    return TMath::Abs(value);
 } 
 
+AliMUONFastTracking::AliMUONFastTracking(const AliMUONFastTracking & ft):TObject()
+{
+// Copy constructor
+    ft.Copy(*this);
+}
+
+
 AliMUONFastTracking* AliMUONFastTracking::Instance()
 { 
 // Set random number generator 
@@ -65,8 +101,18 @@ AliMUONFastTracking* AliMUONFastTracking::Instance()
 
 AliMUONFastTracking::AliMUONFastTracking() 
 {
-//    SetBackground();
-    
+//
+// constructor
+//
+  for (Int_t i = 0; i<20;i++) {
+    for (Int_t j = 0; j<20; j++) {
+      for (Int_t k = 0; k<20; k++) {
+       fFitp[i][j][k] = 0x0;
+      }
+    }
+  }
+
+    fClusterFinder = kOld; 
     fPrintLevel = 1;  
     // read binning; temporarily put by hand
     Float_t pmin = 0, pmax = 200;
@@ -111,26 +157,31 @@ void AliMUONFastTracking::Init(Float_t bkg)
   }
   
   char filename [100]; 
-  sprintf (filename,"$(ALICE_ROOT)/FASTSIM/data/MUONtrackLUT.root"); 
+  if (fClusterFinder==kOld) sprintf (filename,"$(ALICE_ROOT)/FASTSIM/data/MUONtrackLUT.root"); 
+  else sprintf (filename,"$(ALICE_ROOT)/FASTSIM/data/MUONtrackLUT-AZ.root"); 
+
   TFile *file = new TFile(filename); 
   ReadLUT(file);
   SetBackground(bkg);
-  UseSpline(1);
-  fFitp = new TF1("fit1",FitP,-20.,20.,6);
-  fFitp->SetNpx(200);    
+  UseSpline(0);
 }
 
 
 void AliMUONFastTracking::ReadLUT(TFile* file)
 {
+  //
+  // read the lookup tables from file
+  //
   TH3F *heff[5][3], *hacc[5][3], *hmeanp, *hsigmap, *hsigma1p, *hchi2p;
   TH3F *hnormg2, *hmeang2, *hsigmag2, *hmeantheta, *hsigmatheta, *hchi2theta;
   TH3F *hmeanphi, *hsigmaphi, *hchi2phi;
   char tag[40], tag2[40]; 
   
-  const Float_t bkg[4] = {0, 0.5, 1, 2};
+  printf ("Reading parameters from LUT file %s...\n",file->GetName());
+
+  const Float_t kBkg[4] = {0, 0.5, 1, 2};
   for (Int_t ibkg=0; ibkg<4; ibkg++) {
-    sprintf (tag,"BKG%g",bkg[ibkg]); 
+    sprintf (tag,"BKG%g",kBkg[ibkg]); 
     file->cd(tag);
     for (Int_t isplp = 0; isplp<kSplitP; isplp++) { 
       for (Int_t ispltheta = 0; ispltheta<kSplitTheta; ispltheta++) { 
@@ -154,7 +205,6 @@ void AliMUONFastTracking::ReadLUT(TFile* file)
     hsigmaphi   = (TH3F*)gDirectory->Get("hsigmaphi");
     hchi2phi    = (TH3F*)gDirectory->Get("hchi2phi");
     
-    printf ("Reading parameters from LUT file %s...\n",file->GetName());
     for (Int_t ip=0; ip<fNbinp ;ip++) {
       for (Int_t itheta=0; itheta<fNbintheta ;itheta++) {
        for (Int_t iphi=0; iphi<fNbinphi ;iphi++) {
@@ -227,8 +277,11 @@ void AliMUONFastTracking::ReadLUT(TFile* file)
 void AliMUONFastTracking::GetBinning(Int_t &nbinp, Float_t &pmin, Float_t &pmax,
                                     Int_t &nbintheta, Float_t &thetamin, 
                                     Float_t &thetamax,
-                                    Int_t &nbinphi, Float_t &phimin, Float_t &phimax)
+                                    Int_t &nbinphi, Float_t &phimin, Float_t &phimax) const 
 {
+  //
+  // gets the binning for the discrete parametrizations in the lookup table
+  //
     nbinp = fNbinp;
     pmin  = fPmin;
     pmax  = fPmax;
@@ -243,42 +296,33 @@ void AliMUONFastTracking::GetBinning(Int_t &nbinp, Float_t &pmin, Float_t &pmax,
 
 void AliMUONFastTracking::GetIpIthetaIphi(Float_t p, Float_t theta, Float_t phi, 
                                          Int_t charge, Int_t &ip, Int_t &itheta, 
-                                         Int_t &iphi)
+                                         Int_t &iphi) const
 {
+  //
+  // gets the id of the cells in the LUT for a given (p,theta,phi, charge)
+  //
     if (charge < 0) phi = -phi;
     ip           = Int_t (( p - fPmin ) / fDeltaP);
     itheta       = Int_t (( theta - fThetamin ) / fDeltaTheta);
     iphi         = Int_t (( phi - fPhimin ) / fDeltaPhi);
     
-    if (ip< 0) {
-       printf ("Warning: ip= %d. Set to 0\n",ip);
-       ip = 0;
-    } 
-    if (ip>= fNbinp) {
-      //       printf ("Warning: ip = %d. Set to %d\n",ip,fNbinp-1);
-       ip = fNbinp-1;
-    } 
-    if (itheta< 0) {
-      //       printf ("Warning: itheta= %d. Set to 0\n",itheta);
-       itheta = 0;
-    } 
-    if (itheta>= fNbintheta) {
-      //       printf ("Warning: itheta = %d. Set to %d\n",itheta,fNbintheta-1);
-       itheta = fNbintheta-1;
-    } 
     
-    if (iphi< 0) {
-       printf ("Warning: iphi= %d. Set to 0\n",iphi);
-       iphi = 0;
-    } 
-    if (iphi>= fNbinphi) {
-       printf ("Warning: iphi = %d. Set to %d\n",iphi,fNbinphi-1);
-       iphi = fNbinphi-1;
-    } 
+    if (ip< 0)         ip = 0;
+    if (ip>= fNbinp) ip = fNbinp-1;
+    if (itheta< 0) itheta = 0;
+    if (itheta>= fNbintheta) itheta = fNbintheta-1;
+    
+    if (iphi< 0) iphi = 0;
+    if (iphi>= fNbinphi) iphi = fNbinphi-1;
 }
 
 void AliMUONFastTracking::GetSplit(Int_t ip, Int_t itheta, 
-                                  Int_t &nSplitP, Int_t &nSplitTheta) { 
+                                  Int_t &nSplitP, Int_t &nSplitTheta) const 
+{ 
+  //
+  // the first cell is splitted in more bins for theta and momentum
+  // parameterizations. Get the number of divisions for the splitted bins
+  //
   if (ip==0) nSplitP = 5; 
   else nSplitP = 2; 
   if (itheta==0) nSplitTheta = 3; 
@@ -287,6 +331,9 @@ void AliMUONFastTracking::GetSplit(Int_t ip, Int_t itheta,
 
 Float_t AliMUONFastTracking::Efficiency(Float_t p,   Float_t theta, 
                                        Float_t phi, Int_t charge){
+  //
+  // gets the tracking efficiency
+  //
   Int_t ip=0, itheta=0, iphi=0;
   GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
   Int_t nSplitP, nSplitTheta; 
@@ -302,6 +349,9 @@ Float_t AliMUONFastTracking::Efficiency(Float_t p,   Float_t theta,
 
 Float_t AliMUONFastTracking::Acceptance(Float_t p,   Float_t theta, 
                                        Float_t phi, Int_t charge){
+  //
+  // gets the geometrical acceptance
+  //
   if (theta<fThetamin || theta>fThetamax) return 0; 
   
   Int_t ip=0, itheta=0, iphi=0;
@@ -319,16 +369,22 @@ Float_t AliMUONFastTracking::Acceptance(Float_t p,   Float_t theta,
 }
 
 Float_t AliMUONFastTracking::MeanP(Float_t p,   Float_t theta, 
-                           Float_t phi, Int_t charge)
+                           Float_t phi, Int_t charge) const
 {
+  //
+  // gets the mean value of the prec-pgen distribution
+  //
     Int_t ip=0, itheta=0, iphi=0;
     GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
     return fCurrentEntry[ip][itheta][iphi]->fMeanp;
 }
 
 Float_t AliMUONFastTracking::SigmaP(Float_t p,   Float_t theta, 
-                                   Float_t phi, Int_t charge)
+                                   Float_t phi, Int_t charge) const
 {
+  //
+  // gets the width of the prec-pgen distribution
+  //
     Int_t ip=0, itheta=0, iphi=0;
     Int_t index;
     GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
@@ -372,8 +428,11 @@ Float_t AliMUONFastTracking::SigmaP(Float_t p,   Float_t theta,
 }
 
 Float_t AliMUONFastTracking::Sigma1P(Float_t p,   Float_t theta, 
-                           Float_t phi, Int_t charge)
+                           Float_t phi, Int_t charge) const
 {
+  //
+  // gets the width correction of the prec-pgen distribution (see FitP)
+  //
     Int_t ip=0, itheta=0, iphi=0;
     GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
     if (p>fPmax) {
@@ -389,8 +448,12 @@ Float_t AliMUONFastTracking::Sigma1P(Float_t p,   Float_t theta,
 }
 
 Float_t AliMUONFastTracking::NormG2(Float_t p,   Float_t theta, 
-                                   Float_t phi, Int_t charge)
+                                   Float_t phi, Int_t charge) const
 {
+  //
+  // gets the relative normalization of the background
+  // (gaussian) component in the prec-pgen distribution
+  //
     Int_t ip=0, itheta=0, iphi=0;
     GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
     if (p>fPmax) {
@@ -406,8 +469,12 @@ Float_t AliMUONFastTracking::NormG2(Float_t p,   Float_t theta,
 }
 
 Float_t AliMUONFastTracking::MeanG2(Float_t p,   Float_t theta, 
-                                   Float_t phi, Int_t charge)
+                                   Float_t phi, Int_t charge) const
 {
+  //
+  // gets the mean value of the background
+  // (gaussian) component in the prec-pgen distribution
+  //
     Int_t ip=0, itheta=0, iphi=0;
     GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
     if (p>fPmax) {
@@ -423,8 +490,12 @@ Float_t AliMUONFastTracking::MeanG2(Float_t p,   Float_t theta,
 }
 
 Float_t AliMUONFastTracking::SigmaG2(Float_t p,   Float_t theta, 
-                                    Float_t phi, Int_t charge)
+                                    Float_t phi, Int_t charge) const
 {
+  //
+  // gets the width of the background
+  // (gaussian) component in the prec-pgen distribution
+  //
     Int_t ip=0, itheta=0, iphi=0;
     GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
     if (p>fPmax) {
@@ -441,15 +512,22 @@ Float_t AliMUONFastTracking::SigmaG2(Float_t p,   Float_t theta,
 
 
 Float_t AliMUONFastTracking::MeanTheta(Float_t p,   Float_t theta, 
-                                      Float_t phi, Int_t charge)
+                                      Float_t phi, Int_t charge) const
 {
+  //
+  // gets the mean value of the thetarec-thetagen distribution
+  //
     Int_t ip=0, itheta=0, iphi=0;
     GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
     return fCurrentEntry[ip][itheta][iphi]->fMeantheta;
 }
 
-Float_t AliMUONFastTracking::SigmaTheta(Float_t p,   Float_t theta, 
-                           Float_t phi, Int_t charge){
+Float_t AliMUONFastTracking::SigmaTheta(Float_t p,   Float_t theta,  
+                           Float_t phi, Int_t charge) const
+{
+  //
+  // gets the width of the thetarec-thetagen distribution
+  //
   Int_t ip=0, itheta=0, iphi=0;
   Int_t index;
   GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
@@ -484,7 +562,11 @@ Float_t AliMUONFastTracking::SigmaTheta(Float_t p,   Float_t theta,
 
 
 Float_t AliMUONFastTracking::MeanPhi(Float_t p,   Float_t theta, 
-                           Float_t phi, Int_t charge){
+                           Float_t phi, Int_t charge) const
+{
+  //
+  // gets the mean value of the phirec-phigen distribution
+  //
   Int_t ip=0, itheta=0, iphi=0;
   GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
   return fCurrentEntry[ip][itheta][iphi]->fMeanphi;
@@ -492,6 +574,9 @@ Float_t AliMUONFastTracking::MeanPhi(Float_t p,   Float_t theta,
 
 Float_t AliMUONFastTracking::SigmaPhi(Float_t p,   Float_t theta, 
                            Float_t phi, Int_t charge){
+  //
+  // gets the width of the phirec-phigen distribution
+  //
   Int_t ip=0, itheta=0, iphi=0;
   Int_t index;
   GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
@@ -525,6 +610,10 @@ Float_t AliMUONFastTracking::SigmaPhi(Float_t p,   Float_t theta,
 }
 
 void AliMUONFastTracking::SetSpline(){
+  //
+  // sets the spline functions for a smooth behaviour of the parameters
+  // when going from one cell to another
+  //
   printf ("Setting spline functions...");
   char splname[40];
   Double_t x[20][3];
@@ -642,82 +731,22 @@ void AliMUONFastTracking::SetSpline(){
   printf ("...done\n");
 }
   
-void AliMUONFastTracking::SmearMuon(Float_t pgen, Float_t thetagen, Float_t phigen, 
-                           Int_t charge, Float_t &psmear, Float_t &thetasmear,
-                           Float_t &phismear, Float_t &eff, Float_t &acc){
-
-  // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-  // IMPORTANT NOTICE TO THE USER
-  //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-  // THIS METHOD HAS BEEN REPLACED BY AliFastMuonTrackingEff::Evaluate()
-  // AND WILL BE DELETED SOON
-  // DO NOT USE THIS METHOD
-  // 
-
-  printf ("AliMUONFastTracking::SmearMuon()   THIS METHOD IS OBSOLETE ");
-  printf ("PLEASE REFER TO AliFastMuonTrackingEff::Evaluate()\n");
-  // angles are in degrees   
-  
-  Double_t meanp    = MeanP  (pgen, thetagen, phigen, charge);
-  Double_t sigmap   = SigmaP (pgen, thetagen, phigen, charge);
-  Double_t sigma1p  = Sigma1P(pgen, thetagen, phigen, charge);
-  Double_t normg2   = NormG2 (pgen, thetagen, phigen, charge);
-  Double_t meang2   = MeanG2 (pgen, thetagen, phigen, charge);
-  Double_t sigmag2  = SigmaG2(pgen, thetagen, phigen, charge);
-
-  printf ("fBkg = %g    normg2 (100,5,0,1) = %g \n",fBkg,NormG2(100,5,0,1));
-  printf ("fBkg = %g    meang2 (100,5,0,1) = %g \n",fBkg,MeanG2(100,5,0,1));
-  printf ("fBkg = %g    sigmag2 (100,5,0,1) = %g \n",fBkg,SigmaG2(100,5,0,1));
-  Int_t ip,itheta,iphi;
-  GetIpIthetaIphi(pgen, thetagen, phigen, charge, ip, itheta, iphi);
-  if (sigmap == 0) {
-    if (fPrintLevel>0) {
-      printf ("WARNING!!! sigmap=0:    ");
-      printf ("ip= %d itheta = %d iphi = %d   ", ip, itheta, iphi);  
-      printf ("p= %f theta = %f phi = %f\n", pgen, thetagen, phigen);  
-    }
-  }
-  
-  if (fPrintLevel>1) printf ("setting parameters: meanp = %f sigmap = %f sigma1p = %f normg2 = %f meang2 = %f sigmag2 = %f \n",meanp,sigmap,sigma1p,normg2,meang2,sigmag2);
-  fFitp->SetParameters(meanp,sigmap,sigma1p,normg2,meang2,sigmag2);
-  
-  Double_t meantheta  = MeanTheta (pgen, thetagen, phigen, charge);
-  Double_t sigmatheta = SigmaTheta(pgen, thetagen, phigen, charge);
-  Double_t meanphi    = MeanPhi   (pgen, thetagen, phigen, charge);
-  Double_t sigmaphi   = SigmaPhi  (pgen, thetagen, phigen, charge);
-  
-  // components different from ip=0 have the RMS bigger than mean
-  Float_t ptp[3]  =  { 1.219576,-0.354764,-0.690117 };
-  Float_t ptph[3] =  { 0.977522, 0.016269, 0.023158 }; 
-  Float_t pphp[3] =  { 1.303256,-0.464847,-0.869322 };
-  psmear   = pgen + fFitp->GetRandom();
-  Float_t dp = psmear - pgen; 
-  if (ip==0) sigmaphi *= pphp[0] + pphp[1] * dp + pphp[2] * dp*dp; 
-  phismear = phigen + gRandom->Gaus(meanphi, sigmaphi);
-  Float_t dphi = phismear - phigen; 
-
-  if (ip==0) sigmatheta *= ptp[0] + ptp[1] * dp + ptp[2] * dp*dp; 
-  if (ip==0) sigmatheta *= ptph[0] + ptph[1] * dphi + ptph[2] * dphi*dphi; 
-  thetasmear = thetagen + gRandom->Gaus(meantheta,sigmatheta);
-  eff      = Efficiency(pgen, thetagen, phigen, charge);
-  acc      = Acceptance(pgen, thetagen, phigen, charge);
-}
-
 void AliMUONFastTracking::SetBackground(Float_t bkg){
+  //
   // linear interpolation of the parameters in the LUT between 2 values where
   // the background has been actually calculated
-
+  //
   if (bkg>2) printf ("WARNING: unsafe extrapolation!\n");
   fBkg = bkg;
 
-  Float_t BKG[4] = {0, 0.5, 1, 2}; // bkg values for which LUT is calculated
+  Float_t bkgLevel[4] = {0, 0.5, 1, 2}; // bkg values for which LUT is calculated
   Int_t ibkg;
-  for (ibkg=0; ibkg<4; ibkg++) if ( bkg < BKG[ibkg]) break;
+  for (ibkg=0; ibkg<4; ibkg++) if ( bkg < bkgLevel[ibkg]) break;
   if (ibkg == 4) ibkg--;
   if (ibkg == 0) ibkg++;
   
-  Float_t x0 = BKG[ibkg-1];
-  Float_t x1 = BKG[ibkg];
+  Float_t x0 = bkgLevel[ibkg-1];
+  Float_t x1 = bkgLevel[ibkg];
   Float_t x = (bkg - x0) / (x1 - x0); 
   
   Float_t y0, y1;
@@ -777,11 +806,32 @@ void AliMUONFastTracking::SetBackground(Float_t bkg){
   SetSpline();
 }
 
+TF1* AliMUONFastTracking::GetFitP(Int_t ip,Int_t itheta,Int_t iphi) { 
+  // gets the correct prec-pgen distribution for a given LUT cell 
+  if (!fFitp[ip][itheta][iphi]) {
+    char name[256];
+    sprintf(name, "fit_%d_%d_%d", ip, itheta, iphi);
+    fFitp[ip][itheta][iphi] = new TF1(name ,FitP,-20.,20.,6);
+    fFitp[ip][itheta][iphi]->SetNpx(500);    
+    fFitp[ip][itheta][iphi]->SetParameters(0.,0.,0.,0.,0.,0.);    
+  }
+  return fFitp[ip][itheta][iphi]; 
+}
 
+AliMUONFastTracking& AliMUONFastTracking::operator=(const  AliMUONFastTracking& rhs)
+{
+// Assignment operator
+    rhs.Copy(*this);
+    return *this;
+}
 
-  // to guarantee a safe extrapolation for sigmag2 to 0<bkg<0.5, let's fit 
-  // with a straight line sigmag2 vs bkg for bkg=0.5, 1 and 2, and put the 
-  // sigma2(BKG=0) as the extrapolation of this fit 
+void AliMUONFastTracking::Copy(TObject&) const
+{
+    //
+    // Copy 
+    //
+    Fatal("Copy","Not implemented!\n");
+}