First commit.
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 6 Jan 2003 10:13:33 +0000 (10:13 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 6 Jan 2003 10:13:33 +0000 (10:13 +0000)
16 files changed:
EVGEN/AliCollisionGeometry.cxx [new file with mode: 0644]
EVGEN/AliCollisionGeometry.h [new file with mode: 0644]
EVGEN/AliGrayParticleModel.cxx [new file with mode: 0644]
EVGEN/AliGrayParticleModel.h [new file with mode: 0644]
EVGEN/EVGENLinkDef.h
EVGEN/libEVGEN.pkg
FASTSIM/AliFastMuonTrackingAcc.cxx [new file with mode: 0644]
FASTSIM/AliFastMuonTrackingAcc.h [new file with mode: 0644]
FASTSIM/AliFastMuonTrackingEff.cxx [new file with mode: 0644]
FASTSIM/AliFastMuonTrackingEff.h [new file with mode: 0644]
FASTSIM/AliFastMuonTrackingRes.cxx [new file with mode: 0644]
FASTSIM/AliFastMuonTrackingRes.h [new file with mode: 0644]
FASTSIM/AliFastMuonTriggerEff.cxx [new file with mode: 0644]
FASTSIM/AliFastMuonTriggerEff.h [new file with mode: 0644]
FASTSIM/AliMUONFastTracking.cxx [new file with mode: 0644]
FASTSIM/AliMUONFastTracking.h [new file with mode: 0644]

diff --git a/EVGEN/AliCollisionGeometry.cxx b/EVGEN/AliCollisionGeometry.cxx
new file mode 100644 (file)
index 0000000..f518d4d
--- /dev/null
@@ -0,0 +1,22 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/*
+$Log$
+*/
+
+
+#include "AliCollisionGeometry.h"
+ClassImp(AliCollisionGeometry)
diff --git a/EVGEN/AliCollisionGeometry.h b/EVGEN/AliCollisionGeometry.h
new file mode 100644 (file)
index 0000000..612ce78
--- /dev/null
@@ -0,0 +1,55 @@
+#ifndef ALICOLLISIONGEOMETRY_H
+#define ALICOLLISIONGEOMETRY_H
+/* Copyright(c) 198-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+#include <TClassTable.h>
+
+class AliCollisionGeometry
+{
+public:
+    AliCollisionGeometry() {;}
+    virtual ~AliCollisionGeometry(){;}
+    // Getters
+    Float_t ImpactParameter()   {return fImpactParameter;}
+    Int_t   HardScatters() {return fNHardScatters;}
+    Int_t   ProjectileParticipants()  {return fNProjectileParticipants;}
+    Int_t   TargetParticipants()      {return fNTargetParticipants;}
+    Int_t   Spectatorsn()      {return fSpecn;}
+    Int_t   Spectatorsp()      {return fSpecp;}
+    Int_t   NN()    {return fNNColl;}
+    Int_t   NNw()   {return fNNwColl;}
+    Int_t   NwN()   {return fNwNColl;}
+    Int_t   NwNw()  {return fNwNwColl;}
+    // Setters
+    void SetImpactParameter(Float_t b)     {fImpactParameter=b;}
+    void SetHardScatters(Int_t n)  {fNHardScatters=n;}
+    void SetParticipants(Int_t np, Int_t nt)
+       {fNProjectileParticipants=np, fNTargetParticipants=nt;}
+    void SetCollisions(Int_t nn, Int_t nnw, Int_t nwn, Int_t nwnw)
+       {fNNColl=nn, fNNwColl=nnw, fNwNColl=nwn,  fNwNwColl=nwnw;}
+    void SetSpectators(Int_t nspecn, Int_t nspecp)
+       {fSpecn=nspecn, fSpecp=nspecp;}
+ protected:
+    Int_t   fNHardScatters;            // Number of hard scatterings
+    Int_t   fNProjectileParticipants;  // Number of projectiles participants
+    Int_t   fNTargetParticipants;      // Number of target participants
+    Int_t   fNNColl;                   // Number of N-N collisions
+    Int_t   fNNwColl;                  // Number of N-Nwounded collisions
+    Int_t   fNwNColl;                  // Number of Nwounded-N collisons
+    Int_t   fNwNwColl;                 // Number of Nwounded-Nwounded collisions
+    Int_t   fSpecn;                    // Number of spectators neutrons
+    Int_t   fSpecp;                    // Number of spectators protons
+    Float_t fImpactParameter;          // Impact Parameter
+
+  ClassDef(AliCollisionGeometry,1)     // Collision Geometry
+};
+#endif
+
+
+
+
+
+
diff --git a/EVGEN/AliGrayParticleModel.cxx b/EVGEN/AliGrayParticleModel.cxx
new file mode 100644 (file)
index 0000000..9aef18d
--- /dev/null
@@ -0,0 +1,23 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/*
+$Log$
+*/
+
+
+#include "AliGrayParticleModel.h"
+#include "AliCollisionGeometry.h"
+ClassImp(AliGrayParticleModel)
diff --git a/EVGEN/AliGrayParticleModel.h b/EVGEN/AliGrayParticleModel.h
new file mode 100644 (file)
index 0000000..77e9b3b
--- /dev/null
@@ -0,0 +1,27 @@
+#ifndef ALIGRAYPARTICLEMODEL_H
+#define ALIGRAYPARTICLEMODEL_H
+/* Copyright(c) 198-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+#include "TObject.h"
+class AliCollisionGeometry;
+
+class AliGrayParticleModel : public TObject
+{
+public:
+    AliGrayParticleModel() {;}
+    virtual ~AliGrayParticleModel(){;}
+    virtual void GetNumberOfGrayNucleons(AliCollisionGeometry* geo, Int_t& np, Int_t& nn) {;}
+    
+ protected:
+  ClassDef(AliGrayParticleModel,1) // Gray Particle Model
+};
+#endif
+
+
+
+
+
+
index 3f9634e..07da821 100644 (file)
@@ -62,6 +62,8 @@
 #pragma link C++ class  AliStructFuncType+;
 #pragma link C++ class  AliGenGrayParticles+;
 #pragma link C++ class  AliGenEpEmv1+;
+#pragma link C++ class  AliGrayParticleModel+;
+#pragma link C++ class  AliCollisionGeometry+;
 #endif
 
 
index 33280b1..7ae6a2d 100644 (file)
@@ -21,7 +21,8 @@ SRCS          = AliGenHIJINGpara.cxx AliGenBox.cxx AliGenFixed.cxx \
                AliGenAfterBurnerFlow.cxx \
                AliPartonicEnergyLoss.cxx AliGenHerwig.cxx \
                AliStructFuncType.cxx AliGenGrayParticles.cxx \
-        AliGenGeVSimEventHeader.cxx AliGenEpEmv1.cxx
+               AliGenGeVSimEventHeader.cxx AliGenEpEmv1.cxx \
+               AliGrayParticleModel.cxx AliCollisionGeometry.cxx
 
 # Headerfiles for this particular package (Path respect to own directory)
 HDRS= $(SRCS:.cxx=.h) 
@@ -41,7 +42,7 @@ $(ROOTSYS)/include/TParticle.h
 
 DHDR:=EVGENLinkDef.h
 
-EXPORT:=AliPythia.h AliDecayer.h AliDecayerPythia.h AliGenHijingEventHeader.h
+EXPORT:=AliPythia.h AliDecayer.h AliDecayerPythia.h AliGenHijingEventHeader.h AliCollisionGeometry.h
 
 
 
diff --git a/FASTSIM/AliFastMuonTrackingAcc.cxx b/FASTSIM/AliFastMuonTrackingAcc.cxx
new file mode 100644 (file)
index 0000000..0c1d62e
--- /dev/null
@@ -0,0 +1,45 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/*
+$Log$
+*/
+
+#include "AliFastMuonTrackingAcc.h"
+#include "AliMUONFastTracking.h"
+
+ClassImp(AliFastMuonTrackingAcc)
+
+
+AliFastMuonTrackingAcc::AliFastMuonTrackingAcc() :
+    AliFastResponse("Acceptance", "Muon Tracking Acceptance")
+{
+    SetBackground();
+}
+
+void AliFastMuonTrackingAcc::Init()
+{
+    fFastTracking = AliMUONFastTracking::Instance();
+    fFastTracking->Init(fBackground);
+}
+
+
+
+Float_t AliFastMuonTrackingAcc::Evaluate(Float_t pt, Float_t theta, Float_t phi)
+{
+    Float_t p = pt / TMath::Sin(theta*TMath::Pi()/180.);
+    Float_t eff =  fFastTracking->Acceptance(p, theta, phi, Int_t(fCharge));
+    return eff;
+}
diff --git a/FASTSIM/AliFastMuonTrackingAcc.h b/FASTSIM/AliFastMuonTrackingAcc.h
new file mode 100644 (file)
index 0000000..7a14495
--- /dev/null
@@ -0,0 +1,32 @@
+#ifndef ALIFASTMUONTRACKINGACC
+#define ALIFASTMUONTRACKINGACC
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+#include "AliFastResponse.h"
+class AliMUONFastTracking;
+
+class AliFastMuonTrackingAcc :  public AliFastResponse {
+ public:
+    AliFastMuonTrackingAcc();
+    ~AliFastMuonTrackingAcc(){;}
+    void SetBackground(Float_t bg = 1.) {fBackground = bg;}
+    void SetCharge(Float_t charge = 1.) {fCharge     = charge;}
+    virtual void Init();
+    virtual Float_t Evaluate(Float_t pt, Float_t theta, Float_t phi);
+ protected:
+    Float_t              fBackground;   // Background level
+    Float_t              fCharge;       // Current charge
+    
+    AliMUONFastTracking* fFastTracking; //!Pointer to Fast Tracking Data Handler
+    ClassDef(AliFastMuonTrackingAcc,1)  // Fast MUON Tracking Acceptance
+};
+
+#endif
+
+
+
+
+
diff --git a/FASTSIM/AliFastMuonTrackingEff.cxx b/FASTSIM/AliFastMuonTrackingEff.cxx
new file mode 100644 (file)
index 0000000..55aa353
--- /dev/null
@@ -0,0 +1,45 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/*
+$Log$
+*/
+
+#include "AliFastMuonTrackingEff.h"
+#include "AliMUONFastTracking.h"
+
+ClassImp(AliFastMuonTrackingEff)
+
+
+AliFastMuonTrackingEff::AliFastMuonTrackingEff() :
+    AliFastResponse("Efficiency", "Muon Tracking Efficiency")
+{
+    SetBackground();
+}
+
+void AliFastMuonTrackingEff::Init()
+{
+    fFastTracking = AliMUONFastTracking::Instance();
+    fFastTracking->Init(fBackground);
+}
+
+
+
+Float_t AliFastMuonTrackingEff::Evaluate(Float_t pt, Float_t theta, Float_t phi)
+{
+    Float_t p = pt / TMath::Sin(theta*TMath::Pi()/180.);
+    Float_t eff =  fFastTracking->Efficiency(p, theta, phi, Int_t(fCharge));
+    return eff;
+}
diff --git a/FASTSIM/AliFastMuonTrackingEff.h b/FASTSIM/AliFastMuonTrackingEff.h
new file mode 100644 (file)
index 0000000..5facfc6
--- /dev/null
@@ -0,0 +1,32 @@
+#ifndef ALIFASTMUONTRACKINGEFF
+#define ALIFASTMUONTRACKINGEFF
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+#include "AliFastResponse.h"
+class AliMUONFastTracking;
+
+class AliFastMuonTrackingEff :  public AliFastResponse {
+ public:
+    AliFastMuonTrackingEff();
+    ~AliFastMuonTrackingEff(){;}
+    void SetBackground(Float_t bg = 1.) {fBackground = bg;}
+    void SetCharge(Float_t charge = 1.) {fCharge     = charge;}
+    virtual void Init();
+    virtual Float_t Evaluate(Float_t pt, Float_t theta, Float_t phi);
+ protected:
+    Float_t              fBackground;   // Background level
+    Float_t              fCharge;       // Current charge
+    
+    AliMUONFastTracking* fFastTracking; //!Pointer to Fast Tracking Data Handler
+    ClassDef(AliFastMuonTrackingEff,1)  // Fast MUON Tracking Efficiency
+};
+
+#endif
+
+
+
+
+
diff --git a/FASTSIM/AliFastMuonTrackingRes.cxx b/FASTSIM/AliFastMuonTrackingRes.cxx
new file mode 100644 (file)
index 0000000..e192c4e
--- /dev/null
@@ -0,0 +1,93 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/*
+$Log$
+*/
+
+#include "AliFastMuonTrackingRes.h"
+#include "AliMUONFastTracking.h"
+#include <TRandom.h>
+#include <TF1.h>
+
+ClassImp(AliFastMuonTrackingRes)
+
+
+AliFastMuonTrackingRes::AliFastMuonTrackingRes() :
+    AliFastResponse("Resolution", "Muon Tracking Resolution")
+{
+    SetBackground();
+}
+
+void AliFastMuonTrackingRes::Init()
+{
+    fFastTracking = AliMUONFastTracking::Instance();
+    fFastTracking->Init(fBackground);
+}
+
+
+
+void AliFastMuonTrackingRes::Evaluate(Float_t   p,  Float_t  theta , Float_t   phi,
+                                        Float_t& pS,  Float_t& thetaS, Float_t&  phiS)
+{
+
+    Double_t meanp    = fFastTracking->MeanP  (p, theta, phi, Int_t(fCharge));
+    Double_t sigmap   = fFastTracking->SigmaP (p, theta, phi, Int_t(fCharge));
+    Double_t sigma1p  = fFastTracking->Sigma1P(p, theta, phi, Int_t(fCharge));
+    Double_t normg2   = fFastTracking->NormG2 (p, theta, phi, Int_t(fCharge));
+    Double_t meang2   = fFastTracking->MeanG2 (p, theta, phi, Int_t(fCharge));
+    Double_t sigmag2  = fFastTracking->SigmaG2(p, theta, phi, Int_t(fCharge));
+    
+    Int_t ip,itheta,iphi;
+    TF1* fitp = fFastTracking->GetFitP();
+    
+    fFastTracking->GetIpIthetaIphi(p, theta, phi, Int_t(fCharge), ip, itheta, iphi);
+    if (sigmap == 0) {
+       printf ("WARNING!!! sigmap=0:    ");
+       printf ("ip= %d itheta = %d iphi = %d   ", ip, itheta, iphi);  
+       printf ("p= %f theta = %f phi = %f\n", p, theta, phi);  
+    }
+
+    fitp->SetParameters(meanp,sigmap,sigma1p,normg2,meang2,sigmag2);
+  
+    Double_t meantheta  = fFastTracking->MeanTheta (p, theta, phi, Int_t(fCharge));
+    Double_t sigmatheta = fFastTracking->SigmaTheta(p, theta, phi, Int_t(fCharge));
+    Double_t meanphi    = fFastTracking->MeanPhi   (p, theta, phi, Int_t(fCharge));
+    Double_t sigmaphi   = fFastTracking->SigmaPhi  (p, theta, phi, Int_t(fCharge));
+  
+    // 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 };
+
+    // Smeared momentum
+    pS   = p + fitp->GetRandom();
+    Float_t dp = pS - p; 
+
+    // Smeared phi
+    if (ip==0) sigmaphi *= pphp[0] + pphp[1] * dp + pphp[2] * dp*dp; 
+    phiS = phi + gRandom->Gaus(meanphi, sigmaphi);
+    Float_t dphi = phiS - phi; 
+    // Smeared theta
+    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; 
+    thetaS = theta + gRandom->Gaus(meantheta,sigmatheta);
+}
+
+
+
+
+
+
diff --git a/FASTSIM/AliFastMuonTrackingRes.h b/FASTSIM/AliFastMuonTrackingRes.h
new file mode 100644 (file)
index 0000000..be190c4
--- /dev/null
@@ -0,0 +1,30 @@
+#ifndef ALIFASTMUONTRACKINGRES
+#define ALIFASTMUONTRACKINGRES
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+#include "AliFastResponse.h"
+class AliMUONFastTracking;
+
+class AliFastMuonTrackingRes :  public AliFastResponse {
+ public:
+    AliFastMuonTrackingRes();
+    ~AliFastMuonTrackingRes(){;}
+    void SetBackground(Float_t bg = 1.) {fBackground = bg;}
+    void SetCharge(Float_t charge = 1.) {fCharge     = charge;}
+    virtual void Init();
+    virtual void Evaluate(Float_t   p,  Float_t  theta , Float_t   phi,
+                         Float_t& pS,  Float_t& thetaS, Float_t&  phiS);
+ protected:
+    Float_t              fBackground;   // Background level
+    Float_t              fCharge;       // Current charge
+    
+    AliMUONFastTracking* fFastTracking; //!Pointer to Fast Tracking Data Handler
+    ClassDef(AliFastMuonTrackingRes,1)  // Fast MUON Tracking 
+};
+
+#endif
+
+
diff --git a/FASTSIM/AliFastMuonTriggerEff.cxx b/FASTSIM/AliFastMuonTriggerEff.cxx
new file mode 100644 (file)
index 0000000..2a698cf
--- /dev/null
@@ -0,0 +1,370 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/*
+$Log$
+*/
+
+
+#include "AliFastMuonTriggerEff.h"
+#include <TFile.h>
+#include <TTree.h>
+#include <TROOT.h>
+#include <stdlib.h>
+
+// Debugging flag
+//#define MYTRIGDEBUG
+
+ClassImp(AliFastMuonTriggerEff);
+
+AliFastMuonTriggerEff::AliFastMuonTriggerEff():
+    AliFastResponse("Efficiency", "Muon Trigger Efficiency")
+{
+//
+// Default constructor
+//
+    SetCut(kLow);
+    fZones=0;
+}
+
+
+void AliFastMuonTriggerEff::Init()
+{
+//
+//  Initialization
+//
+    printf("Initializing %s / %s", fName.Data(), fTitle.Data());
+//
+
+    fPhiMin   = -90. ;
+    fPhiMax   =  90. ;
+    fDphi     =   9. ;
+    fThetaMin =   2. ;
+    fThetaMax =   9. ;
+    fDtheta   =   0.7;
+    fDpt      =   0.1;
+
+    InitTree();
+}
+
+void AliFastMuonTriggerEff::InitTree()
+{
+//
+//  Initialize tables from tree
+//
+    TTree          *chain;   // pointer to the analyzed TTree or TChain
+    Int_t           nSim, nZona;
+    Double_t        pmin, pmax, tmin, tmax, phimin, phimax, bkg;
+    Double_t        len50, hen50, leff, heff;
+    Double_t        vLPt[50];
+    Double_t        vHPt[50];
+    Char_t file[100]="$(ALICE_ROOT)/data/vettorpara.root";
+//
+//  Avoid memory leak in case of reinitialization
+    if(fZones!=0) {
+        printf("\nWarning: reinitialization of an object of class: AliFastMuonTriggerEff\n");
+        for (Int_t i=0; i<fZones; i++) {
+            if(fEffLow [i])delete[] fEffLow[i];
+           if(fEffHigh[i])delete[] fEffHigh[i];
+        }
+       if(fEffLow) {
+           delete[] fEffLow;
+           fEffLow=0;
+        }
+       if(fEffHigh) {
+           delete[] fEffHigh;
+           fEffHigh=0;
+        }
+    }
+    printf("AliFastMuonTriggerEff: Initialization\n");
+    TFile *f = new TFile(file);
+    if(f->IsZombie()) {
+        printf("Cannot open file: %s\n",file);
+        return;
+    }
+    f->ls();
+    
+    TTree* tree = (TTree*)gDirectory->Get("fitPar");
+
+//
+//
+    chain = tree;
+    chain->SetMakeClass(1);
+    
+    chain->SetBranchAddress("nSim",&nSim);
+    chain->SetBranchAddress("nZona",&nZona);
+    chain->SetBranchAddress("ptmin",&pmin);
+    chain->SetBranchAddress("ptmax",&pmax);
+    chain->SetBranchAddress("Thetamin",&tmin);
+    chain->SetBranchAddress("Thetamax",&tmax);
+    chain->SetBranchAddress("Phimin",&phimin);
+    chain->SetBranchAddress("Phimax",&phimax);
+    chain->SetBranchAddress("Bkg",&bkg);
+    chain->SetBranchAddress("EffLPt",vLPt);
+    chain->SetBranchAddress("EffHPt",vHPt);
+    chain->SetBranchAddress("Pt0.5Low",&len50);
+    chain->SetBranchAddress("Pt0.5High",&hen50);
+    chain->SetBranchAddress("EffLowMax",&leff);
+    chain->SetBranchAddress("EffHighMax",&heff);
+// 
+//
+    Int_t nentries = Int_t(chain->GetEntries());
+    Int_t nbytes = 0, nb = 0;
+
+// Count the number of zones of the parametrization
+    Int_t nzone0=0, nzone1=0;
+    for (Int_t jentry=0; jentry<nentries; jentry++) {
+       nb = chain->GetEntry(jentry);
+       if(bkg==0.) {
+            if(nSim==0)nzone0++;
+            if(nSim==1)nzone1++;
+       }
+    }
+    
+    printf("Trigger parametrisation for %d zones for Pt: 0. 5. GeV/c\n",nzone0);
+    printf("and %d zones extended to 10. GeV/c\n",nzone1);
+    fZones=nzone0+nzone1;
+//    printf("Ciao\n");
+    if(fZones<=0){
+        printf("Number of zones must be greater than 0\n");
+       exit(6);
+    }
+    
+    fEffLow =new Float_t*[fZones];
+    fEffHigh=new Float_t*[fZones];
+    for (Int_t i=0; i<fZones; i++) {
+        fEffLow [i]=new Float_t[50];
+       fEffHigh[i]=new Float_t[50];
+    }
+    
+//  Initialize look-up table to standard values
+    Int_t isim, itheta, iphi;
+    for (isim=0; isim<2; isim++) {
+        for (itheta=0; itheta<10; itheta++) {
+           for (iphi=0; iphi<20; iphi++) {
+               fLook[isim][itheta][iphi]=0;
+            }
+       }
+    }
+
+//  Loading Trigger efficiencies
+    Int_t myzone=0;
+#ifdef MYTRIGDEBUG
+            printf("Filling nSim nZona pmin pmax tmin tmax phimin phimax: ....\n");
+#endif
+    for (Int_t jentry=0; jentry<nentries; jentry++) {
+//     Int_t ientry = LoadTree(jentry); 
+       nb = chain->GetEntry(jentry);   
+       nbytes += nb;
+#ifdef MYTRIGDEBUG
+       printf("Getting entry %d... ",jentry);
+#endif
+// For the time being it works with background 0
+       if ((nSim == 0 || nSim == 1)&&bkg==0.) {
+#ifdef MYTRIGDEBUG
+            printf("Filling %d %d %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f: ",
+           nSim,nZona,pmin,pmax,tmin,tmax,phimin,phimax);
+#endif
+           for (Int_t k = 0; k < 50; k++) {
+               fEffLow [myzone][k] = vLPt[k];
+               fEffHigh[myzone][k] = vHPt[k];
+#ifdef MYTRIGDEBUG
+               if(k<15)printf(" %5.3f",vLPt[k]);
+#endif
+           }
+#ifdef MYTRIGDEBUG
+            printf("\n");
+#endif
+           myzone++;
+           iphi=Int_t(nZona/10)-10;
+           itheta=nZona%10;
+           if(iphi<0||iphi>19||itheta<0||itheta>9) {
+               printf("The format of data file is not consistent\nwith this version of the code\n");
+               printf("This should never happen: iphi %d, itheta: %d\n",iphi,itheta);
+               exit(7);
+           }
+           fLook[nSim][itheta][iphi]=myzone;
+       } else {
+           printf("Skipping entry with nSim=%d, bkg=%f\n",nSim,bkg);
+       }
+    }
+#ifdef MYTRIGDEBUG
+    printf("This is the content of the LUT after first step of initialization\n");
+    for(isim=0; isim<2; isim++) {
+        printf("isim=%d\n",isim);
+        for(iphi=0; iphi<20; iphi++) {
+           printf("iphi: %2d:",iphi);
+           for(itheta=0; itheta<10; itheta++) {
+               printf(" %4d",fLook[isim][itheta][iphi]);
+            }
+           printf("\n");
+       }
+    }
+#endif
+// Filling look=up table for the zones where the extended simulation does
+// not exists
+    for(iphi=0; iphi<20; iphi++) {
+       for(itheta=0; itheta<10; itheta++) {
+           if(fLook[0][itheta][iphi]==0) {
+               printf("Missing entry isim=%d itheta=%d iphi=%d\n",isim,itheta,iphi);
+               exit(8);
+           }
+           if(fLook[0][itheta][iphi]<0||fLook[0][itheta][iphi]>fZones) {
+               printf("Problem with entry isim=%d itheta=%d iphi=%d\n",isim,itheta,iphi);
+               exit(9);
+           }
+           if(fLook[1][itheta][iphi]==0) {
+                 fLook[1][itheta][iphi]=-fLook[0][itheta][iphi];
+           }
+       }
+    }
+#ifdef MYTRIGDEBUG
+    for(isim=0; isim<2; isim++) {
+        printf("isim=%d\n",isim);
+        for(iphi=0; iphi<20; iphi++) {
+           printf("iphi: %2d:",iphi);
+           for(itheta=0; itheta<10; itheta++) {
+               printf(" %4d",fLook[isim][itheta][iphi]);
+            }
+           printf("\n");
+       }
+    }
+        for(iphi=0; iphi<20; iphi++) {
+           for(itheta=0; itheta<10; itheta++) {
+    for(isim=0; isim<2; isim++) {
+               printf("%d %2d %2d %4d:",isim,iphi,itheta,fLook[isim][itheta][iphi]);
+               if(fLook[isim][itheta][iphi]>0) {
+                   myzone=fLook[isim][itheta][iphi]-1;
+                   for(Int_t ipt=0; ipt<20; ipt++) {
+                       //printf(" %5.3f",fEffLow[myzone][ipt]);
+                       printf(" %5.3f",fEffHigh[myzone][ipt]);
+                   }
+                   printf(" ...");
+                   for(Int_t ipt=40; ipt<50; ipt++) {
+                       //printf(" %5.3f",fEffLow[myzone][ipt]);
+                       printf(" %5.3f",fEffHigh[myzone][ipt]);
+                   }               
+               }
+               printf("\n");
+            }
+       }
+    }    
+#endif
+    f->Close();
+}
+
+void AliFastMuonTriggerEff::Evaluate(Float_t charge, Float_t pt,
+                Float_t theta, Float_t phi,Float_t& effLow, Float_t& effHigh)
+{
+//
+//  Trigger efficiency for pt, theta, phi (low and high cut)
+//
+    effLow=0.;
+    effHigh=0.;
+    if(fZones==0) {
+        printf("Call to uninitialized object of class: AliFastMuonTriggerEff\n");
+       return;
+    }
+    if(pt<0) {
+        printf("Warning: pt: %f < 0. GeV/c\n",pt);
+       return; 
+    }
+    Int_t iPt   = Int_t(pt/fDpt)%50;
+    Int_t iSim  = Int_t(pt/fDpt)/50;
+    Int_t iPhi  = Int_t((phi-fPhiMin)/fDphi);
+    if(phi<fPhiMin)iPhi=iPhi-1;
+    Int_t iTheta = Int_t((theta-fThetaMin)/fDtheta);
+#ifdef MYTRIGDEBUG
+    printf("iSim iPt iTheta iPhi: %d %d %d %d\n",iSim,iPt,iTheta,iPhi);
+#endif
+    iPhi=iPhi-40*(iPhi/40);
+#ifdef MYTRIGDEBUG
+    printf("1: iPhi converted to: %d for angle equivalence\n",iPhi);
+#endif
+    if(iPhi<0)iPhi=-iPhi-1;
+    if(iPhi>19)iPhi=39-iPhi;
+#ifdef MYTRIGDEBUG
+    printf("2: iPhi converted to: %d for the symmetry of the spectrometer\n",iPhi);
+#endif
+    if(charge==1.){
+    } else if(charge==-1.) {
+    iPhi=19-iPhi;
+#ifdef MYTRIGDEBUG
+    printf("3: iPhi converted to: %d for charge symmetry\n",iPhi);
+#endif
+    } else {
+        printf("Warning: not understand charge: %f\n",charge);
+        return;
+    }
+    if(iTheta<0||iTheta>9) {
+        printf("Warning: theta: %f outside acceptance\n",theta);
+        return;
+    }
+    if(iPt<0) {
+        printf("Warning: what do you mean with pt: %f <0?\n",pt);
+       return;
+    }
+    if(iSim>=fSim) {
+        iSim=fSim-1;
+       iPt=49;
+#ifdef MYTRIGDEBUG
+    printf("4: iSim iPt converted to: %d %d (last zone)\n",iSim,iPt);
+#endif
+    }
+    Int_t iLook=fLook[iSim][iTheta][iPhi];
+    if(iLook<0) {
+#ifdef MYTRIGDEBUG
+    printf("5: iLook iPt: %d %d converted to: ",iLook,iPt);
+#endif
+        iLook=-iLook-1;
+       iPt=49;
+#ifdef MYTRIGDEBUG
+    printf("%d %d from look up table contents\n",iLook,iPt);
+#endif
+    } else {
+        iLook=iLook-1;
+    }
+    effLow=fEffLow[iLook][iPt];
+    effHigh=fEffHigh[iLook][iPt];
+#ifdef MYTRIGDEBUG
+    printf("6: charge, iSim, iTheta, iPhi, iPt: %f %d %d %d %d effLow: %f, effHigh: %f\n",
+               charge,iSim,iTheta,iPhi,iPt,effLow,effHigh);
+#endif
+    
+    //fEffLow [iPhi][iTheta][iPt];
+    //fEffHigh[iPhi][iTheta][iPt];    
+}
+
+Float_t AliFastMuonTriggerEff::Evaluate(Float_t charge, Float_t pt,
+                   Float_t theta, Float_t phi)
+{
+    //
+    // Trigger efficiency for pt, theta, phi depending of fCut
+    // 
+    if(fZones==0) {
+        printf("Call to uninitialized object of class: AliFastMuonTriggerEff\n");
+       return 0.;
+    }
+    Float_t eff;
+    Float_t effLow, effHigh;
+    
+    Evaluate(charge,pt,theta,phi,effLow,effHigh);
+    if (fCut == kLow) 
+       eff  = effLow;
+    else
+       eff  = effHigh;
+
+    return eff;
+}
diff --git a/FASTSIM/AliFastMuonTriggerEff.h b/FASTSIM/AliFastMuonTriggerEff.h
new file mode 100644 (file)
index 0000000..8700c41
--- /dev/null
@@ -0,0 +1,46 @@
+#ifndef ALIFASTMUONTRIGGEREFF_H
+#define ALIFASTMUONTRIGGEREFF_H
+/*  Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+#include <AliFastResponse.h>
+enum CutTupe {kLow, kHigh};
+
+class AliFastMuonTriggerEff : public AliFastResponse {
+    
+ public:
+    AliFastMuonTriggerEff();
+    AliFastMuonTriggerEff(char* Name, char* Title) {;}    
+    virtual ~AliFastMuonTriggerEff(){;}
+    virtual void    Init();
+    virtual void    Evaluate(Float_t charge, Float_t pt, Float_t theta, Float_t phi,
+                            Float_t& effLow, Float_t& effHigh);
+    virtual Float_t Evaluate(Float_t charge, Float_t pt, Float_t theta, Float_t phi);
+    virtual void    SetCut(Int_t cut = kLow) {fCut = cut;}
+    virtual Float_t Cut() {return fCut;}
+  protected:
+    virtual void InitTree();
+  protected:
+    Int_t fLook[2][10][20];       // Look up table for bkg=0
+    Float_t fDpt;                 // Delta_pt
+    Float_t fPhiMin;              // lower limit for phi 
+    Float_t fPhiMax;              // upper limit for phi
+    Float_t fDphi;                // Delta_phi
+    Float_t fThetaMin;            // lower limit for theta
+    Float_t fThetaMax;            // upper limit for theta
+    Float_t fDtheta;              // Delta_theta
+    Int_t   fCut;                 // Cut type (low/high)
+    Int_t   fZones;               // Total number of zones
+    const static Int_t   fSim=2;  // Number of pt extentions (internal use)
+    Float_t** fEffLow;            // Table for low-pt  cut bkg=0
+    Float_t** fEffHigh;           // Table for high-pt cut bkg=0
+    
+    ClassDef(AliFastMuonTriggerEff,1)    // Fast Muon Trigger response
+};
+
+#endif 
+
+
+
diff --git a/FASTSIM/AliMUONFastTracking.cxx b/FASTSIM/AliMUONFastTracking.cxx
new file mode 100644 (file)
index 0000000..3ebd123
--- /dev/null
@@ -0,0 +1,789 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/*
+$Log$
+*/
+
+#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>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <Riostream.h>
+
+ClassImp(AliMUONFastTracking)
+
+AliMUONFastTracking* AliMUONFastTracking::fgMUONFastTracking=NULL;
+
+static Double_t FitP(Double_t *x, Double_t *par){
+    Double_t dx = x[0] - par[0];
+    Double_t dx2 = x[0] - par[4];
+    Double_t sigma = par[1] * ( 1 + par[2] * dx);
+    if (sigma == 0) {    
+
+       return 0.;
+    }
+    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);
+} 
+
+AliMUONFastTracking* AliMUONFastTracking::Instance()
+{ 
+// Set random number generator 
+    if (fgMUONFastTracking) {
+       return fgMUONFastTracking;
+    } else {
+       fgMUONFastTracking = new AliMUONFastTracking();
+       return fgMUONFastTracking;
+    }
+}
+
+AliMUONFastTracking::AliMUONFastTracking() 
+{
+//    SetBackground();
+    
+    fPrintLevel = 1;  
+    // read binning; temporarily put by hand
+    Float_t pmin = 0, pmax = 200;
+    Int_t   nbinp = 10;
+    Float_t thetamin = 2, thetamax = 9;
+    Int_t   nbintheta=10;
+    Float_t phimin = -180, phimax =180;
+    Int_t   nbinphi=10;
+    //--------------------------------------
+
+    fNbinp = nbinp;
+    fPmin  = pmin;
+    fPmax  = pmax;
+    
+    fNbintheta = nbintheta;
+    fThetamin  = thetamin;
+    fThetamax  = thetamax;
+    
+    fNbinphi = nbinphi;
+    fPhimin  = phimin;
+    fPhimax  = phimax;
+    
+    fDeltaP     = (fPmax-fPmin)/fNbinp;
+    fDeltaTheta = (fThetamax-fThetamin)/fNbintheta;
+    fDeltaPhi   = (fPhimax-fPhimin)/fNbinphi;
+}
+
+void AliMUONFastTracking::Init(Float_t bkg)
+{
+  //
+  //  Initialization
+  //
+  for (Int_t ip=0; ip< fNbinp; ip++){
+    for (Int_t itheta=0; itheta< fNbintheta; itheta++){
+      for (Int_t iphi=0; iphi< fNbinphi; iphi++){
+       fCurrentEntry[ip][itheta][iphi] = new AliMUONFastTrackingEntry;
+       for (Int_t ibkg=0; ibkg<4; ibkg++){
+         fEntry[ip][itheta][iphi][ibkg] = new AliMUONFastTrackingEntry;
+       }
+      }
+    }
+  }
+  
+  char filename [100]; 
+  sprintf (filename,"$(ALICE_ROOT)/data/MUONtrackLUT.root"); 
+  TFile *file = new TFile(filename); 
+  ReadLUT(file);
+  SetBackground(bkg);
+  UseSpline(1);
+  fFitp = new TF1("fit1",FitP,-20.,20.,6);
+  fFitp->SetNpx(200);    
+}
+
+
+void AliMUONFastTracking::ReadLUT(TFile* 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};
+  for (Int_t ibkg=0; ibkg<4; ibkg++) {
+    sprintf (tag,"BKG%g",bkg[ibkg]); 
+    file->cd(tag);
+    for (Int_t isplp = 0; isplp<kSplitP; isplp++) { 
+      for (Int_t ispltheta = 0; ispltheta<kSplitTheta; ispltheta++) { 
+       sprintf (tag2,"heff[%d][%d]",isplp,ispltheta); 
+       heff[isplp][ispltheta] = (TH3F*)gDirectory->Get(tag2);
+       sprintf (tag2,"hacc[%d][%d]",isplp,ispltheta); 
+       hacc[isplp][ispltheta] = (TH3F*)gDirectory->Get(tag2);
+      }
+    }    
+    hmeanp      = (TH3F*)gDirectory->Get("hmeanp");
+    hsigmap     = (TH3F*)gDirectory->Get("hsigmap");
+    hsigma1p    = (TH3F*)gDirectory->Get("hsigma1p");
+    hchi2p      = (TH3F*)gDirectory->Get("hchi2p");
+    hnormg2     = (TH3F*)gDirectory->Get("hnormg2");
+    hmeang2     = (TH3F*)gDirectory->Get("hmeang2");
+    hsigmag2    = (TH3F*)gDirectory->Get("hsigmag2");
+    hmeantheta  = (TH3F*)gDirectory->Get("hmeantheta");
+    hsigmatheta = (TH3F*)gDirectory->Get("hsigmatheta");
+    hchi2theta  = (TH3F*)gDirectory->Get("hchi2theta");
+    hmeanphi    = (TH3F*)gDirectory->Get("hmeanphi");
+    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++) {
+         Float_t p     = fPmin     + fDeltaP     * (ip     + 0.5);
+         Float_t theta = fThetamin + fDeltaTheta * (itheta + 0.5);
+         Float_t phi   = fPhimin   + fDeltaPhi   * (iphi   + 0.5);
+         
+         fEntry[ip][itheta][iphi][ibkg]->fP          = p;
+         fEntry[ip][itheta][iphi][ibkg]->fMeanp      = 
+           hmeanp->GetBinContent(ip+1,itheta+1,iphi+1);
+         fEntry[ip][itheta][iphi][ibkg]->fSigmap     = 
+           TMath::Abs(hsigmap->GetBinContent(ip+1,itheta+1,iphi+1));
+         fEntry[ip][itheta][iphi][ibkg]->fSigma1p    = 
+           hsigma1p->GetBinContent(ip+1,itheta+1,iphi+1);
+         fEntry[ip][itheta][iphi][ibkg]->fChi2p      = 
+           hchi2p->GetBinContent(ip+1,itheta+1,iphi+1);
+         fEntry[ip][itheta][iphi][ibkg]->fNormG2     = 
+           hnormg2->GetBinContent(ip+1,itheta+1,iphi+1);
+         fEntry[ip][itheta][iphi][ibkg]->fMeanG2     = 
+           hmeang2->GetBinContent(ip+1,itheta+1,iphi+1);
+         if (ibkg == 0)  fEntry[ip][itheta][iphi][ibkg]->fSigmaG2 = 9999;
+         else fEntry[ip][itheta][iphi][ibkg]->fSigmaG2    = 
+           hsigmag2->GetBinContent(ip+1,itheta+1,iphi+1);
+         fEntry[ip][itheta][iphi][ibkg]->fTheta      = theta;
+         fEntry[ip][itheta][iphi][ibkg]->fMeantheta  = 
+           hmeantheta->GetBinContent(ip+1,itheta+1,iphi+1);
+         fEntry[ip][itheta][iphi][ibkg]->fSigmatheta = 
+           TMath::Abs(hsigmatheta->GetBinContent(ip+1,itheta+1,iphi+1));
+         fEntry[ip][itheta][iphi][ibkg]->fChi2theta  = 
+           hchi2theta->GetBinContent(ip+1,itheta+1,iphi+1);
+         fEntry[ip][itheta][iphi][ibkg]->fPhi        = phi;
+         fEntry[ip][itheta][iphi][ibkg]->fMeanphi    = 
+           hmeanphi->GetBinContent(ip+1,itheta+1,iphi+1);
+         fEntry[ip][itheta][iphi][ibkg]->fSigmaphi   = 
+           TMath::Abs(hsigmaphi->GetBinContent(ip+1,itheta+1,iphi+1));
+         fEntry[ip][itheta][iphi][ibkg]->fChi2phi    = 
+           hchi2phi->GetBinContent(ip+1,itheta+1,iphi+1);
+         for (Int_t i=0; i<kSplitP; i++) { 
+           for (Int_t j=0; j<kSplitTheta; j++) { 
+             fEntry[ip][itheta][iphi][ibkg]->fAcc[i][j] = 
+               hacc[i][j]->GetBinContent(ip+1,itheta+1,iphi+1);
+             fEntry[ip][itheta][iphi][ibkg]->fEff[i][j] = 
+               heff[i][j]->GetBinContent(ip+1,itheta+1,iphi+1);
+           }
+         }
+       } // iphi 
+      }  // itheta
+    }   // ip
+  }    // ibkg
+
+  TGraph *graph = new TGraph(3); 
+  TF1 *f = new TF1("f","[0]+[1]*x"); 
+  
+  for (Int_t ip=0; ip< fNbinp; ip++){
+    for (Int_t itheta=0; itheta< fNbintheta; itheta++){
+      for (Int_t iphi=0; iphi< fNbinphi; iphi++){
+       graph->SetPoint(0,0.5,fEntry[ip][itheta][iphi][1]->fSigmaG2);
+       graph->SetPoint(1,1,fEntry[ip][itheta][iphi][2]->fSigmaG2);
+       graph->SetPoint(2,2,fEntry[ip][itheta][iphi][3]->fSigmaG2);
+       graph->Fit("f","q"); 
+       fEntry[ip][itheta][iphi][0]->fSigmaG2 = f->Eval(0); 
+      }
+    }
+  }
+  f->Delete();
+  graph->Delete();
+  printf ("parameters read. \n");
+}
+
+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)
+{
+    nbinp = fNbinp;
+    pmin  = fPmin;
+    pmax  = fPmax;
+    nbintheta = fNbintheta;
+    thetamin  = fThetamin;
+    thetamax  = fThetamax;
+    nbinphi = fNbinphi;
+    phimin  = fPhimin;
+    phimax  = fPhimax;
+}
+
+
+void AliMUONFastTracking::GetIpIthetaIphi(Float_t p, Float_t theta, Float_t phi, 
+                                         Int_t charge, Int_t &ip, Int_t &itheta, 
+                                         Int_t &iphi)
+{
+    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;
+    } 
+}
+
+void AliMUONFastTracking::GetSplit(Int_t ip, Int_t itheta, 
+                                  Int_t &nSplitP, Int_t &nSplitTheta) { 
+  if (ip==0) nSplitP = 5; 
+  else nSplitP = 2; 
+  if (itheta==0) nSplitTheta = 3; 
+  else nSplitTheta = 1; 
+}
+
+Float_t AliMUONFastTracking::Efficiency(Float_t p,   Float_t theta, 
+                                       Float_t phi, Int_t charge){
+  Int_t ip=0, itheta=0, iphi=0;
+  GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
+  Int_t nSplitP, nSplitTheta; 
+  GetSplit(ip,itheta,nSplitP,nSplitTheta); 
+
+  Float_t dp = p - fPmin;
+  Int_t ibinp  = Int_t(nSplitP*(dp - fDeltaP * Int_t(dp / fDeltaP))/fDeltaP); 
+  Float_t dtheta = theta - fThetamin;
+  Int_t ibintheta  = Int_t(nSplitTheta*(dtheta - fDeltaTheta * Int_t(dtheta / fDeltaTheta))/fDeltaTheta); 
+  Float_t eff = fCurrentEntry[ip][itheta][iphi]->fEff[ibinp][ibintheta];
+  return eff;   
+}
+
+Float_t AliMUONFastTracking::Acceptance(Float_t p,   Float_t theta, 
+                                       Float_t phi, Int_t charge){
+  if (theta<fThetamin || theta>fThetamax) return 0; 
+  
+  Int_t ip=0, itheta=0, iphi=0;
+  GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
+  Int_t nSplitP, nSplitTheta; 
+  GetSplit(ip,itheta,nSplitP,nSplitTheta); 
+  // central value and corrections with spline 
+  
+  Float_t dp = p - fPmin;
+  Int_t ibinp  = Int_t(nSplitP*(dp - fDeltaP * Int_t(dp / fDeltaP))/fDeltaP); 
+  Float_t dtheta = theta - fThetamin;
+  Int_t ibintheta  = Int_t(nSplitTheta*(dtheta - fDeltaTheta * Int_t(dtheta / fDeltaTheta))/fDeltaTheta); 
+  Float_t acc = fCurrentEntry[ip][itheta][iphi]->fAcc[ibinp][ibintheta];
+  return acc;   
+}
+
+Float_t AliMUONFastTracking::MeanP(Float_t p,   Float_t theta, 
+                           Float_t phi, Int_t charge)
+{
+    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)
+{
+    Int_t ip=0, itheta=0, iphi=0;
+    Int_t index;
+    GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
+    // central value and corrections with spline 
+    Float_t sigmap = fCurrentEntry[ip][itheta][iphi]->fSigmap;
+    if (!fSpline) return sigmap;
+    // corrections vs p, theta, phi 
+    index = iphi + fNbinphi * itheta;
+    Double_t xmin,ymin,xmax,ymax;
+    Float_t frac1 = fSplineSigmap[index][0]->Eval(p)/sigmap; 
+    
+    if (p>fPmax-fDeltaP/2.) {
+       Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->fSigmap;
+       Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->fSigmap;
+       Float_t s3 = fCurrentEntry[fNbinp-3][itheta][iphi]->fSigmap;
+       Float_t p1  = fDeltaP * (fNbinp - 1 + 0.5) + fPmin;
+       Float_t p2  = fDeltaP * (fNbinp - 2 + 0.5) + fPmin;
+       Float_t p3  = fDeltaP * (fNbinp - 3 + 0.5) + fPmin;
+       Float_t p12 = p1 * p1, p22 = p2 * p2, p32 = p3 * p3;
+       Float_t d = p12*p2 + p1*p32 + p22*p3 - p32*p2 - p3*p12 - p22*p1;
+       Float_t a = (s1*p2 + p1*s3 + s2*p3 - s3*p2 - p3*s1 - s2*p1) / d;
+       Float_t b = (p12*s2 + s1*p32 + p22*s3 - p32*s2 - s3*p12 - p22*s1)/d;
+       Float_t c = (p12*p2*s3 + p1*p32*s2 + p22*p3*s1 
+                    - p32*p2*s1 - p3*p12*s2 - p22*p1*s3) / d;
+       Float_t sigma = a * p * p + b * p + c;
+       frac1 = sigma/sigmap;
+    }
+    index = iphi + fNbinphi * ip;
+    fSplineEff[index][1]->GetKnot(0,xmin,ymin);
+    fSplineEff[index][1]->GetKnot(9,xmax,ymax);
+    if (theta>xmax) theta = xmax;
+    Float_t frac2 = fSplineSigmap[index][1]->Eval(theta)/sigmap; 
+    index = itheta + fNbintheta * ip;
+    fSplineEff[index][2]->GetKnot(0,xmin,ymin);
+    fSplineEff[index][2]->GetKnot(9,xmax,ymax);
+    if (phi>xmax) phi = xmax;
+    Float_t frac3 = fSplineSigmap[index][2]->Eval(phi)/sigmap;
+    Float_t sigmatot = sigmap * frac1 * frac2 * frac3;
+    if (sigmatot<0) sigmatot = sigmap;
+    return sigmatot;
+}
+
+Float_t AliMUONFastTracking::Sigma1P(Float_t p,   Float_t theta, 
+                           Float_t phi, Int_t charge)
+{
+    Int_t ip=0, itheta=0, iphi=0;
+    GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
+    if (p>fPmax) {
+       // linear extrapolation of sigmap for p out of range 
+       Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->fSigma1p;
+       Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->fSigma1p;
+       Float_t p1  = fDeltaP * (fNbinp - 1 + 0.5) + fPmin;
+       Float_t p2  = fDeltaP * (fNbinp - 2 + 0.5) + fPmin;
+       Float_t sigma = 1./(p1-p2) * ( (s1-s2)*p + (s2-s1)*p1 + s1*(p1-p2) ); 
+       return sigma;
+    }
+    else return fCurrentEntry[ip][itheta][iphi]->fSigma1p;
+}
+
+Float_t AliMUONFastTracking::NormG2(Float_t p,   Float_t theta, 
+                                   Float_t phi, Int_t charge)
+{
+    Int_t ip=0, itheta=0, iphi=0;
+    GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
+    if (p>fPmax) {
+       // linear extrapolation of sigmap for p out of range 
+       Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->fNormG2;
+       Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->fNormG2;
+       Float_t p1  = fDeltaP * (fNbinp - 1 + 0.5) + fPmin;
+       Float_t p2  = fDeltaP * (fNbinp - 2 + 0.5) + fPmin;
+       Float_t norm = 1./(p1-p2) * ( (s1-s2)*p + (s2-s1)*p1 + s1*(p1-p2) ); 
+       return norm;
+    }
+    else return fCurrentEntry[ip][itheta][iphi]->fNormG2;
+}
+
+Float_t AliMUONFastTracking::MeanG2(Float_t p,   Float_t theta, 
+                                   Float_t phi, Int_t charge)
+{
+    Int_t ip=0, itheta=0, iphi=0;
+    GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
+    if (p>fPmax) {
+       // linear extrapolation of sigmap for p out of range 
+       Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->fMeanG2;
+       Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->fMeanG2;
+       Float_t p1  = fDeltaP * (fNbinp - 1 + 0.5) + fPmin;
+       Float_t p2  = fDeltaP * (fNbinp - 2 + 0.5) + fPmin;
+       Float_t norm = 1./(p1-p2) * ( (s1-s2)*p + (s2-s1)*p1 + s1*(p1-p2) ); 
+       return norm;
+    }
+    else return fCurrentEntry[ip][itheta][iphi]->fMeanG2;
+}
+
+Float_t AliMUONFastTracking::SigmaG2(Float_t p,   Float_t theta, 
+                                    Float_t phi, Int_t charge)
+{
+    Int_t ip=0, itheta=0, iphi=0;
+    GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
+    if (p>fPmax) {
+       // linear extrapolation of sigmap for p out of range 
+       Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->fSigmaG2;
+       Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->fSigmaG2;
+       Float_t p1  = fDeltaP * (fNbinp - 1 + 0.5) + fPmin;
+       Float_t p2  = fDeltaP * (fNbinp - 2 + 0.5) + fPmin;
+       Float_t sigma = 1./(p1-p2) * ( (s1-s2)*p + (s2-s1)*p1 + s1*(p1-p2) ); 
+       return sigma;
+    }
+    else return fCurrentEntry[ip][itheta][iphi]->fSigmaG2;
+}
+
+
+Float_t AliMUONFastTracking::MeanTheta(Float_t p,   Float_t theta, 
+                                      Float_t phi, Int_t charge)
+{
+    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){
+  Int_t ip=0, itheta=0, iphi=0;
+  Int_t index;
+  GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
+  // central value and corrections with spline 
+  Float_t sigmatheta = fCurrentEntry[ip][itheta][iphi]->fSigmatheta;
+  if (!fSpline) return sigmatheta;
+  // corrections vs p, theta, phi 
+  index = iphi + fNbinphi * itheta;
+  Double_t xmin,ymin,xmax,ymax;
+  Float_t frac1 = fSplineSigmatheta[index][0]->Eval(p)/sigmatheta; 
+  if (p>fPmax-fDeltaP/2.) {
+    // linear extrapolation of sigmap for p out of range 
+    Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->fSigmatheta;
+    Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->fSigmatheta;
+    Float_t p1  = fDeltaP * (fNbinp - 1 + 0.5) + fPmin;
+    Float_t p2  = fDeltaP * (fNbinp - 2 + 0.5) + fPmin;
+    Float_t sigma = 1./(p1-p2) * ( (s1-s2)*p + (s2-s1)*p1 + s1*(p1-p2) ); 
+    frac1=sigma/sigmatheta;
+  }
+  index = iphi + fNbinphi * ip;
+  fSplineEff[index][1]->GetKnot(0,xmin,ymin);
+  fSplineEff[index][1]->GetKnot(9,xmax,ymax);
+  if (theta>xmax) theta = xmax;
+  Float_t frac2 = fSplineSigmatheta[index][1]->Eval(theta)/sigmatheta; 
+  index = itheta + fNbintheta * ip;
+  fSplineEff[index][2]->GetKnot(0,xmin,ymin);
+  fSplineEff[index][2]->GetKnot(9,xmax,ymax);
+  if (phi>xmax) phi = xmax;
+  Float_t frac3 = fSplineSigmatheta[index][2]->Eval(phi)/sigmatheta;
+  return sigmatheta * frac1 * frac2 * frac3;
+}
+
+
+Float_t AliMUONFastTracking::MeanPhi(Float_t p,   Float_t theta, 
+                           Float_t phi, Int_t charge){
+  Int_t ip=0, itheta=0, iphi=0;
+  GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
+  return fCurrentEntry[ip][itheta][iphi]->fMeanphi;
+}
+
+Float_t AliMUONFastTracking::SigmaPhi(Float_t p,   Float_t theta, 
+                           Float_t phi, Int_t charge){
+  Int_t ip=0, itheta=0, iphi=0;
+  Int_t index;
+  GetIpIthetaIphi(p,theta,phi,charge,ip,itheta,iphi);
+  // central value and corrections with spline 
+  Float_t sigmaphi = fCurrentEntry[ip][itheta][iphi]->fSigmaphi;
+  if (!fSpline) return sigmaphi;
+  // corrections vs p, theta, phi 
+  index = iphi + fNbinphi * itheta;
+  Float_t frac1 = fSplineSigmaphi[index][0]->Eval(p)/sigmaphi; 
+  Double_t xmin,ymin,xmax,ymax;
+  if (p>fPmax-fDeltaP/2.) {
+    Float_t s1 = fCurrentEntry[fNbinp-1][itheta][iphi]->fSigmaphi;
+    Float_t s2 = fCurrentEntry[fNbinp-2][itheta][iphi]->fSigmaphi;
+    Float_t p1  = fDeltaP * (fNbinp - 1 + 0.5) + fPmin;
+    Float_t p2  = fDeltaP * (fNbinp - 2 + 0.5) + fPmin;
+    Float_t sigma = 1./(p1-p2) * ( (s1-s2)*p + (s2-s1)*p1 + s1*(p1-p2) ); 
+    frac1 = sigma/sigmaphi;
+  }
+  
+  index = iphi + fNbinphi * ip;
+  fSplineEff[index][1]->GetKnot(0,xmin,ymin);
+  fSplineEff[index][1]->GetKnot(9,xmax,ymax);
+  if (theta>xmax) theta = xmax;
+  Float_t frac2 = fSplineSigmaphi[index][1]->Eval(theta)/sigmaphi; 
+  index = itheta + fNbintheta * ip;
+  fSplineEff[index][2]->GetKnot(0,xmin,ymin);
+  fSplineEff[index][2]->GetKnot(9,xmax,ymax);
+  if (phi>xmax) phi = xmax;
+  Float_t frac3 = fSplineSigmaphi[index][2]->Eval(phi)/sigmaphi;
+  return sigmaphi * frac1 * frac2 * frac3;
+}
+
+void AliMUONFastTracking::SetSpline(){
+  printf ("Setting spline functions...");
+  char splname[40];
+  Double_t x[20][3];
+  Double_t x2[50][3];
+  Int_t nbins[3] = {fNbinp, fNbintheta, fNbinphi};
+  Double_t xspl[20],yeff[50],ysigmap[20],ysigma1p[20];
+  Double_t yacc[50], ysigmatheta[20],ysigmaphi[20];
+  Double_t xsp2[50];
+  // let's calculate the x axis for p, theta, phi
+  
+  Int_t i, ispline, ivar;
+  for (i=0; i< fNbinp; i++) x[i][0] = fPmin + fDeltaP * (i + 0.5);
+  for (i=0; i< fNbintheta; i++) x[i][1] = fThetamin + fDeltaTheta * (i + 0.5);
+  for (i=0; i< fNbinphi; i++) x[i][2] = fPhimin + fDeltaPhi * (i + 0.5);
+  
+  for (i=0; i< 5 * fNbinp; i++) x2[i][0] = fPmin + fDeltaP * (i + 0.5)/5.;
+  for (i=0; i< 5 * fNbintheta; i++) x2[i][1] = fThetamin + fDeltaTheta * (i + 0.5)/5.;
+  for (i=0; i< 5 * fNbinphi; i++) x2[i][2] = fPhimin + fDeltaPhi * (i + 0.5)/5.;
+  
+  // splines in p
+  ivar = 0; 
+  for (i=0; i<nbins[ivar]; i++) xspl[i] = x[i][ivar];
+  for (i=0; i<5 * nbins[ivar]; i++) xsp2[i] = x2[i][ivar];
+  ispline=0;
+  for (Int_t itheta=0; itheta< fNbintheta; itheta++){
+    for (Int_t iphi=0; iphi< fNbinphi; iphi++){
+      for (Int_t ip=0; ip<fNbinp; ip++) {
+       //      for (Int_t i=0; i<5; i++) { 
+       //        yeff[5 * ip + i]  = fCurrentEntry[ip][itheta][iphi]->fEff[i];
+       //        yacc[5 * ip + i]  = fCurrentEntry[ip][itheta][iphi]->fAcc[i];
+       //      }
+       ysigmap[ip]     = fCurrentEntry[ip][itheta][iphi]->fSigmap;
+       ysigma1p[ip]    = fCurrentEntry[ip][itheta][iphi]->fSigma1p;
+       ysigmatheta[ip] = fCurrentEntry[ip][itheta][iphi]->fSigmatheta;
+       ysigmaphi[ip]   = fCurrentEntry[ip][itheta][iphi]->fSigmaphi;
+      }
+      if (fPrintLevel>3) cout << " creating new spline " << splname << endl;
+      sprintf (splname,"fSplineEff[%d][%d]",ispline,ivar);
+      fSplineEff[ispline][ivar] = new TSpline3(splname,xsp2,yeff,5 * nbins[ivar]);
+      sprintf (splname,"fSplineAcc[%d][%d]",ispline,ivar);
+      fSplineAcc[ispline][ivar] = new TSpline3(splname,xsp2,yacc,5 * nbins[ivar]);
+      sprintf (splname,"fSplineSigmap[%d][%d]",ispline,ivar);
+      fSplineSigmap[ispline][ivar] = new TSpline3(splname,xspl,ysigmap,nbins[ivar]);
+      sprintf (splname,"fSplineSigma1p[%d][%d]",ispline,ivar);
+      fSplineSigma1p[ispline][ivar] = new TSpline3(splname,xspl,ysigma1p,nbins[ivar]);
+      sprintf (splname,"fSplineSigmatheta[%d][%d]",ispline,ivar);
+      fSplineSigmatheta[ispline][ivar] = new TSpline3(splname,xspl,ysigmatheta,nbins[ivar]);
+      sprintf (splname,"fSplineSigmaphi[%d][%d]",ispline,ivar);
+      fSplineSigmaphi[ispline][ivar] = new TSpline3(splname,xspl,ysigmaphi,nbins[ivar]);
+      ispline++;
+    }
+  }
+
+  ivar = 1; 
+  for (i=0; i<nbins[ivar]; i++) xspl[i] = x[i][ivar];
+  ispline=0;
+  for (Int_t ip=0; ip<fNbinp; ip++) {
+    for (Int_t iphi=0; iphi< fNbinphi; iphi++){
+      for (Int_t itheta=0; itheta< fNbintheta; itheta++){
+       // for efficiency and acceptance let's take the central value
+       //      yeff[itheta]        = fCurrentEntry[ip][itheta][iphi]->fEff[2];
+       //      yacc[itheta]        = fCurrentEntry[ip][itheta][iphi]->fAcc[2];
+       ysigmap[itheta]     = fCurrentEntry[ip][itheta][iphi]->fSigmap;
+       ysigma1p[itheta]    = fCurrentEntry[ip][itheta][iphi]->fSigma1p;
+       ysigmatheta[itheta] = fCurrentEntry[ip][itheta][iphi]->fSigmatheta;
+       ysigmaphi[itheta]   = fCurrentEntry[ip][itheta][iphi]->fSigmaphi;
+      }
+      if (fPrintLevel>3) cout << " creating new spline " << splname << endl;
+      sprintf (splname,"fSplineEff[%d][%d]",ispline,ivar);
+      fSplineEff[ispline][ivar] = new TSpline3(splname,xspl,yeff, nbins[ivar]);
+      sprintf (splname,"fSplineAcc[%d][%d]",ispline,ivar);
+      fSplineAcc[ispline][ivar] = new TSpline3(splname,xspl,yacc, nbins[ivar]);
+      sprintf (splname,"fSplineSigmap[%d][%d]",ispline,ivar);
+      fSplineSigmap[ispline][ivar] = new TSpline3(splname,xspl,ysigmap,nbins[ivar]);
+      sprintf (splname,"fSplineSigma1p[%d][%d]",ispline,ivar);
+      fSplineSigma1p[ispline][ivar] = new TSpline3(splname,xspl,ysigma1p,nbins[ivar]);
+      sprintf (splname,"fSplineSigmatheta[%d][%d]",ispline,ivar);
+      fSplineSigmatheta[ispline][ivar] = new TSpline3(splname,xspl,ysigmatheta,nbins[ivar]);
+      sprintf (splname,"fSplineSigmaphi[%d][%d]",ispline,ivar);
+      fSplineSigmaphi[ispline][ivar] = new TSpline3(splname,xspl,ysigmaphi,nbins[ivar]);
+      ispline++;
+    }
+  }
+
+  ivar = 2; 
+  for (i=0; i<nbins[ivar]; i++) xspl[i] = x[i][ivar];
+  ispline=0;
+  for (Int_t ip=0; ip<fNbinp; ip++) {
+    for (Int_t itheta=0; itheta< fNbintheta; itheta++){
+      for (Int_t iphi=0; iphi< fNbinphi; iphi++){
+       // for efficiency and acceptance let's take the central value
+       //      yeff[iphi]        = fCurrentEntry[ip][itheta][iphi]->fEff[2];
+       //      yacc[iphi]        = fCurrentEntry[ip][itheta][iphi]->fAcc[2];
+       ysigmap[iphi]     = fCurrentEntry[ip][itheta][iphi]->fSigmap;
+       ysigma1p[iphi]    = fCurrentEntry[ip][itheta][iphi]->fSigma1p;
+       ysigmatheta[iphi] = fCurrentEntry[ip][itheta][iphi]->fSigmatheta;
+       ysigmaphi[iphi]   = fCurrentEntry[ip][itheta][iphi]->fSigmaphi;
+      }
+      if (fPrintLevel>3) cout << " creating new spline " << splname << endl;
+      sprintf (splname,"fSplineEff[%d][%d]",ispline,ivar);
+      fSplineEff[ispline][ivar] = new TSpline3(splname,xspl,yeff, nbins[ivar]);
+      sprintf (splname,"fSplineAcc[%d][%d]",ispline,ivar);
+      fSplineAcc[ispline][ivar] = new TSpline3(splname,xspl,yacc, nbins[ivar]);
+      sprintf (splname,"fSplineSigmap[%d][%d]",ispline,ivar);
+      fSplineSigmap[ispline][ivar] = new TSpline3(splname,xspl,ysigmap,nbins[ivar]);
+      sprintf (splname,"fSplineSigma1p[%d][%d]",ispline,ivar);
+      fSplineSigma1p[ispline][ivar] = new TSpline3(splname,xspl,ysigma1p,nbins[ivar]);
+      sprintf (splname,"fSplineSigmatheta[%d][%d]",ispline,ivar);
+      fSplineSigmatheta[ispline][ivar] = new TSpline3(splname,xspl,ysigmatheta,nbins[ivar]);
+      sprintf (splname,"fSplineSigmaphi[%d][%d]",ispline,ivar);
+      fSplineSigmaphi[ispline][ivar] = new TSpline3(splname,xspl,ysigmaphi,nbins[ivar]);
+      ispline++;
+    }
+  }
+  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
+  Int_t ibkg;
+  for (ibkg=0; ibkg<4; ibkg++) if ( bkg < BKG[ibkg]) break;
+  if (ibkg == 4) ibkg--;
+  if (ibkg == 0) ibkg++;
+  
+  Float_t x0 = BKG[ibkg-1];
+  Float_t x1 = BKG[ibkg];
+  Float_t x = (bkg - x0) / (x1 - x0); 
+  
+  Float_t y0, y1;
+
+  for (Int_t ip=0; ip< fNbinp; ip++){
+    for (Int_t itheta=0; itheta< fNbintheta; itheta++){
+      for (Int_t iphi=0; iphi< fNbinphi; iphi++){
+       fCurrentEntry[ip][itheta][iphi]->fP = fEntry[ip][itheta][iphi][ibkg]->fP;
+       fCurrentEntry[ip][itheta][iphi]->fTheta = fEntry[ip][itheta][iphi][ibkg]->fTheta;
+       fCurrentEntry[ip][itheta][iphi]->fPhi = fEntry[ip][itheta][iphi][ibkg]->fPhi;
+       fCurrentEntry[ip][itheta][iphi]->fChi2p = -1;
+       fCurrentEntry[ip][itheta][iphi]->fChi2theta = -1;
+       fCurrentEntry[ip][itheta][iphi]->fChi2phi = -1;
+
+       y0 = fEntry[ip][itheta][iphi][ibkg-1]->fMeanp;
+       y1 =   fEntry[ip][itheta][iphi][ibkg]->fMeanp;
+       fCurrentEntry[ip][itheta][iphi]      ->fMeanp = (y1 - y0) * x + y0;
+       y0 = fEntry[ip][itheta][iphi][ibkg-1]->fMeantheta;
+       y1 =   fEntry[ip][itheta][iphi][ibkg]->fMeantheta;
+       fCurrentEntry[ip][itheta][iphi]      ->fMeantheta = (y1 - y0) * x + y0;
+       y0 = fEntry[ip][itheta][iphi][ibkg-1]->fMeanphi;
+       y1 =   fEntry[ip][itheta][iphi][ibkg]->fMeanphi;
+       fCurrentEntry[ip][itheta][iphi]      ->fMeanphi = (y1 - y0) * x + y0;
+       y0 = fEntry[ip][itheta][iphi][ibkg-1]->fSigmap;
+       y1 =   fEntry[ip][itheta][iphi][ibkg]->fSigmap;
+       fCurrentEntry[ip][itheta][iphi]      ->fSigmap = (y1 - y0) * x + y0;
+       y0 = fEntry[ip][itheta][iphi][ibkg-1]->fSigmatheta;
+       y1 =   fEntry[ip][itheta][iphi][ibkg]->fSigmatheta;
+       fCurrentEntry[ip][itheta][iphi]      ->fSigmatheta = (y1 - y0) * x + y0;
+       y0 = fEntry[ip][itheta][iphi][ibkg-1]->fSigmaphi;
+       y1 =   fEntry[ip][itheta][iphi][ibkg]->fSigmaphi;
+       fCurrentEntry[ip][itheta][iphi]      ->fSigmaphi = (y1 - y0) * x + y0;
+       y0 = fEntry[ip][itheta][iphi][ibkg-1]->fSigma1p;
+       y1 =   fEntry[ip][itheta][iphi][ibkg]->fSigma1p;
+       fCurrentEntry[ip][itheta][iphi]      ->fSigma1p = (y1 - y0) * x + y0;
+       y0 = fEntry[ip][itheta][iphi][ibkg-1]->fNormG2;
+       y1 =   fEntry[ip][itheta][iphi][ibkg]->fNormG2;
+       fCurrentEntry[ip][itheta][iphi]      ->fNormG2 = (y1 - y0) * x + y0;
+       y0 = fEntry[ip][itheta][iphi][ibkg-1]->fMeanG2;
+       y1 =   fEntry[ip][itheta][iphi][ibkg]->fMeanG2;
+       fCurrentEntry[ip][itheta][iphi]      ->fMeanG2 = (y1 - y0) * x + y0;
+       
+       y0 = fEntry[ip][itheta][iphi][ibkg-1]->fSigmaG2;
+       y1 =   fEntry[ip][itheta][iphi][ibkg]->fSigmaG2;
+       fCurrentEntry[ip][itheta][iphi]      ->fSigmaG2 = (y1 - y0) * x + y0;
+       for (Int_t i=0; i<kSplitP; i++) {
+         for (Int_t j=0; j<kSplitTheta; j++) {
+           fCurrentEntry[ip][itheta][iphi]->fAcc[i][j] = fEntry[ip][itheta][iphi][ibkg]->fAcc[i][j];
+           y0 = fEntry[ip][itheta][iphi][ibkg-1]->fEff[i][j];
+           y1 =   fEntry[ip][itheta][iphi][ibkg]->fEff[i][j];
+           fCurrentEntry[ip][itheta][iphi]->fEff[i][j] = (y1 - y0) * x + y0;
+         }
+       }
+      }
+    }
+  }
+  SetSpline();
+}
+
+
+
+  // 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 
+
+
+
+
+
+
+
+
diff --git a/FASTSIM/AliMUONFastTracking.h b/FASTSIM/AliMUONFastTracking.h
new file mode 100644 (file)
index 0000000..ec56bc3
--- /dev/null
@@ -0,0 +1,84 @@
+#ifndef ALIMUONFASTTRACKING
+#define ALIMUONFASTTRACKING
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+class TF1;
+class TSpline3;
+class TFile;
+class AliMUONFastTrackingEntry;
+
+
+#include <TObject.h>
+
+class AliMUONFastTracking :  public TObject {
+ public:
+    static  AliMUONFastTracking* Instance();
+    ~AliMUONFastTracking(){;}
+    void Init(Float_t bkg);
+    void ReadLUT(TFile *file);
+    void 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);
+    void GetIpIthetaIphi(Float_t p, Float_t theta, Float_t phi, Int_t charge,
+                        Int_t &ip, Int_t &itheta, Int_t &iphi);
+    void GetSplit(Int_t ip, Int_t itheta, Int_t &nSplitP, Int_t &nSplitTheta);
+    Float_t Efficiency(Float_t p, Float_t theta, Float_t phi, Int_t charge);
+    Float_t Acceptance(Float_t p, Float_t theta, Float_t phi, Int_t charge);
+    Float_t MeanP(Float_t p, Float_t theta, Float_t phi, Int_t charge);
+    Float_t SigmaP(Float_t p, Float_t theta, Float_t phi, Int_t charge);
+    Float_t Sigma1P(Float_t p, Float_t theta, Float_t phi, Int_t charge);
+    Float_t NormG2(Float_t p, Float_t theta, Float_t phi, Int_t charge);
+    Float_t MeanG2(Float_t p, Float_t theta, Float_t phi, Int_t charge);
+    Float_t SigmaG2(Float_t p, Float_t theta, Float_t phi, Int_t charge);
+    Float_t MeanTheta(Float_t p, Float_t theta, Float_t phi, Int_t charge);
+    Float_t SigmaTheta(Float_t p, Float_t theta, Float_t phi, Int_t charge);  
+    Float_t MeanPhi(Float_t p, Float_t theta, Float_t phi, Int_t charge);
+    Float_t SigmaPhi(Float_t p, Float_t theta, Float_t phi, Int_t charge);
+
+    void SetSpline();
+    Float_t GetBackground() {return fBkg;}
+    void SetBackground(Float_t bkg);
+    void UseSpline (Int_t splineSwitch=1) {fSpline = splineSwitch;}
+    void 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);
+    TF1* GetFitP() {return fFitp;}
+ private:
+    AliMUONFastTracking();
+    AliMUONFastTracking(Float_t bkg){;}
+ protected:
+    Int_t   fNentries;
+    Int_t   fNbinp; 
+    Float_t fPmin;
+    Float_t fPmax;
+    Float_t fDeltaP;
+    Int_t   fNbintheta;
+    Float_t fThetamin;
+    Float_t fThetamax;
+    Float_t fDeltaTheta;
+    Int_t   fNbinphi;
+    Float_t fPhimin;
+    Float_t fPhimax;
+    Float_t fDeltaPhi;
+    Int_t   fPrintLevel;
+    Float_t fBkg;
+    TF1 *fFitp;                                   // func for psmear-pgen distr
+    AliMUONFastTrackingEntry *fEntry[20][20][20][4]; // array of LUT parameters
+    AliMUONFastTrackingEntry *fCurrentEntry[20][20][20]; // array of LUT parameters
+ public:
+    TSpline3 *fSplineEff[200][3];                 // spline funcs for efficiency
+    TSpline3 *fSplineAcc[200][3];                 // spline funcs for acceptance
+    TSpline3 *fSplineSigmap[200][3];              // 
+    TSpline3 *fSplineSigma1p[200][3];             //!
+    TSpline3 *fSplineSigmatheta[200][3];          //!
+    TSpline3 *fSplineSigmaphi[200][3];            //!
+ protected: 
+    Int_t fSpline;
+    static AliMUONFastTracking*    fgMUONFastTracking; //!Pointer to single instance
+    ClassDef(AliMUONFastTracking,1)                    // Fast MUON Tracking Data Handler
+};
+
+#endif
+