PWGLF/SPECTRA/Nuclei/B2/AliLnEfficiency.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 // efficiency correction
  17 // author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
  18
  19 #include <Riostream.h>
  20 #include <TFile.h>
  21 #include <TH1D.h>
  22 #include <TH2D.h>
  23 #include <TString.h>
  24 #include <TF1.h>
  25
  26 #include "AliLnEfficiency.h"
  27 #include "B2.h"
  28
  29 ClassImp(AliLnEfficiency)
  30
  31 AliLnEfficiency::AliLnEfficiency(const TString& particle, const TString& simuFilename, const TString& outputFilename, const TString& otag)
  32 : TObject()
  33 , fParticle(particle)
  34 , fSimuFilename(simuFilename)
  35 , fOutputFilename(outputFilename)
  36 , fOutputTag(otag)
  37 , fG3Fluka(0)
  38 , fAddFakeTracks(1)
  39 {
  40 //
  41 // constructor
  42 //
  43         TH1::SetDefaultSumw2();
  44 }
  45
  46 AliLnEfficiency::~AliLnEfficiency()
  47 {
  48 //
  49 // destructor
  50 //
  51 }
  52
  53 Int_t AliLnEfficiency::Exec()
  54 {
  55 //
  56 // extract the efficiency for the given simulation
  57 //
  58         using namespace std;
  59
  60         TFile* fsimu = new TFile(fSimuFilename.Data(),"read");
  61         if (fsimu->IsZombie()) exit(1);
  62
  63         TFile* foutput = new TFile(fOutputFilename.Data(),"recreate");
  64         if(fOutputTag != "")
  65         {
  66                 foutput->mkdir(fOutputTag.Data());
  67                 foutput->cd(fOutputTag.Data());
  68         }
  69
  70         TH1D* hGenPhSpPt = FindObj<TH1D>(fsimu, fParticle + "_Gen_PhS_Prim_Pt");
  71         TH1D* hGenTrigPt = FindObj<TH1D>(fsimu, fParticle + "_Gen_Trig_Prim_Pt");
  72         TH1D* hGenVtxPt  = FindObj<TH1D>(fsimu, fParticle + "_Gen_Vtx_Prim_Pt");
  73         TH1D* hGenAccPt  = FindObj<TH1D>(fsimu, fParticle + "_Gen_Acc_Prim_Pt");
  74         TH1D* hPrimRecPt = FindObj<TH1D>(fsimu, fParticle + "_Sim_Prim_Pt");
  75         TH1D* hFakePrimRecPt = FindObj<TH1D>(fsimu, fParticle + "_Sim_Fake_Prim_Pt");
  76
  77         if(fAddFakeTracks) hPrimRecPt->Add(hFakePrimRecPt);
  78
  79         TH1D* hEffTrigPt   = Divide(hGenTrigPt, hGenPhSpPt, fParticle + "_Eff_Trig_Pt");
  80         TH1D* hEffVtxPt    = Divide(hGenVtxPt, hGenTrigPt, fParticle + "_Eff_Vtx_Pt");
  81         TH1D* hEffAccPt    = Divide(hGenAccPt, hGenVtxPt, fParticle + "_Eff_Acc_Pt");
  82         TH1D* hEffTrkPt    = Divide(hPrimRecPt, hGenAccPt, fParticle + "_Eff_Trk_Pt");
  83         TH1D* hEffAccTrkPt = Divide(hPrimRecPt, hGenVtxPt, fParticle + "_Eff_AccTrk_Pt");
  84
  85         if(fG3Fluka)
  86         {
  87                 TF1* fncG3Fluka = 0;
  88                 if(fParticle == "Proton") fncG3Fluka = GetProtonG3FlukaCor("Proton_G3FLUKA_Corr_Pt");
  89                 if(fParticle == "AntiProton") fncG3Fluka = GetAntiProtonG3FlukaCor("AntiProton_G3FLUKA_Corr_Pt");
  90
  91                 if(fncG3Fluka != 0)
  92                 {
  93                         hEffTrkPt->Divide(fncG3Fluka);
  94                         hEffAccTrkPt->Divide(fncG3Fluka);
  95                 }
  96
  97                 delete fncG3Fluka;
  98         }
  99
 100         if(!fParticle.Contains("Anti"))
 101         {
 102                 TH1D* hStats = FindObj<TH1D>(fsimu, fParticle + "_Stats");
 103                 hStats->Write();
 104         }
 105
 106         hEffAccTrkPt->SetTitle(fParticle.Data());
 107         hEffTrkPt->SetTitle(fParticle.Data());
 108         hEffAccPt->SetTitle(fParticle.Data());
 109         hEffVtxPt->SetTitle(fParticle.Data());
 110         hEffTrigPt->SetTitle(fParticle.Data());
 111
 112         hEffAccTrkPt->SetYTitle("Acceptance #times Track Efficiency");
 113         hEffAccPt->SetYTitle("Acceptance Efficiency");
 114         hEffTrkPt->SetYTitle("Track Efficiency");
 115         hEffVtxPt->SetYTitle("Vertex Efficiency");
 116         hEffTrigPt->SetYTitle("Trigger Efficiency");
 117
 118         hEffAccTrkPt->Write();
 119         hEffAccPt->Write();
 120         hEffTrkPt->Write();
 121         hEffVtxPt->Write();
 122         hEffTrigPt->Write();
 123
 124         delete hEffAccTrkPt;
 125         delete hEffAccPt;
 126         delete hEffTrkPt;
 127         delete hEffVtxPt;
 128         delete hEffTrigPt;
 129
 130         delete foutput;
 131         delete fsimu;
 132
 133         return 0;
 134 }
 135
 136 TF1* AliLnEfficiency::GetProtonG3FlukaCor(const TString& name) const
 137 {
 138 //
 139 // GEANT3/FLUKA correction for protons (Baryon Twiki)
 140 // (TPC tracks)
 141 //
 142         TF1* fnc = new TF1(name.Data(), "1-[0]*exp([1]*x)+[2]", 0.001, 10);
 143         fnc->SetParameters(4113, -32.9, -0.009383);
 144         Double_t err[] = {1640.4, 1.5, 0.000917};
 145         fnc->SetParErrors(err);
 146
 147         return fnc;
 148 }
 149
 150 TF1* AliLnEfficiency::GetAntiProtonG3FlukaCor(const TString& name) const
 151 {
 152 //
 153 // GEANT3/FLUKA correction for antiprotons (Baryon Twiki)
 154 // (TPC tracks)
 155 //
 156         TF1* fnc = new TF1(name.Data(), "1-[0]*exp([1]*x)+[2]+[3]*log(x)/pow(x,0.2)", 0.001, 10);
 157         fnc->SetParameters(177.4, -22.03, -0.06538, 0.04431 );
 158         Double_t err[] = {73., 1.5, 0.00221, 0.00555};
 159         fnc->SetParErrors(err);
 160
 161         return fnc;
 162 }
 163