B2 analysis code
[u/mrichter/AliRoot.git] / 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 {
39 //
40 // constructor
41 //
42 }
43
44 AliLnEfficiency::~AliLnEfficiency()
45 {
46 //
47 // destructor
48 //
49 }
50
51 Int_t AliLnEfficiency::Exec()
52 {
53 //
54 // extract the efficiency for the given simulation
55 //
56         using namespace std;
57         
58         TFile* fsimu = new TFile(fSimuFilename.Data(),"read");
59         if (fsimu->IsZombie()) exit(1);
60         
61         TFile* foutput = new TFile(fOutputFilename.Data(),"recreate");
62         if(fOutputTag != "")
63         {
64                 foutput->mkdir(fOutputTag.Data());
65                 foutput->cd(fOutputTag.Data());
66         }
67         
68         TH1D* hGenPhSpPt = (TH1D*)FindObj(fsimu, fParticle + "_Gen_PhS_Prim_Pt");
69         TH1D* hGenTrigPt = (TH1D*)FindObj(fsimu, fParticle + "_Gen_Trig_Prim_Pt");
70         TH1D* hGenVtxPt  = (TH1D*)FindObj(fsimu, fParticle + "_Gen_Vtx_Prim_Pt");
71         TH1D* hGenAccPt  = (TH1D*)FindObj(fsimu, fParticle + "_Gen_Acc_Prim_Pt");
72         TH1D* hPrimRecPt = (TH1D*)FindObj(fsimu, fParticle + "_Sim_Prim_Pt");
73         
74         TH1D* hEffTrigPt   = Divide(hGenTrigPt, hGenPhSpPt, fParticle + "_Eff_Trig_Pt");
75         TH1D* hEffVtxPt    = Divide(hGenVtxPt, hGenTrigPt, fParticle + "_Eff_Vtx_Pt");
76         TH1D* hEffAccPt    = Divide(hGenAccPt, hGenVtxPt, fParticle + "_Eff_Acc_Pt");
77         TH1D* hEffTrkPt    = Divide(hPrimRecPt, hGenAccPt, fParticle + "_Eff_Trk_Pt");
78         TH1D* hEffAccTrkPt = Divide(hPrimRecPt, hGenVtxPt, fParticle + "_Eff_AccTrk_Pt");
79         
80         if(fG3Fluka)
81         {
82                 TF1* fncG3Fluka = 0;
83                 if(fParticle == "Proton") fncG3Fluka = GetProtonG3FlukaCor("Proton_G3FLUKA_Corr_Pt");
84                 if(fParticle == "AntiProton") fncG3Fluka = GetAntiProtonG3FlukaCor("AntiProton_G3FLUKA_Corr_Pt");
85                 
86                 if(fncG3Fluka != 0)
87                 {
88                         hEffTrkPt->Divide(fncG3Fluka);
89                         hEffAccTrkPt->Divide(fncG3Fluka);
90                 }
91                 
92                 delete fncG3Fluka;
93         }
94         
95         hEffAccTrkPt->SetTitle(fParticle.Data());
96         hEffTrkPt->SetTitle(fParticle.Data());
97         hEffAccPt->SetTitle(fParticle.Data());
98         hEffVtxPt->SetTitle(fParticle.Data());
99         hEffTrigPt->SetTitle(fParticle.Data());
100         
101         hEffAccTrkPt->SetYTitle("Acceptance #times Track Efficiency");
102         hEffAccPt->SetYTitle("Acceptance Efficiency");
103         hEffTrkPt->SetYTitle("Track Efficiency");
104         hEffVtxPt->SetYTitle("Vertex Efficiency");
105         hEffTrigPt->SetYTitle("Trigger Efficiency");
106         
107         hEffAccTrkPt->Write();
108         hEffAccPt->Write();
109         hEffTrkPt->Write();
110         hEffVtxPt->Write();
111         hEffTrigPt->Write();
112         
113         delete hEffAccTrkPt;
114         delete hEffAccPt;
115         delete hEffTrkPt;
116         delete hEffVtxPt;
117         delete hEffTrigPt;
118         
119         delete foutput;
120         delete fsimu;
121         
122         return 0;
123 }
124
125 TF1* AliLnEfficiency::GetProtonG3FlukaCor(const TString& name) const
126 {
127 //
128 // GEANT3/FLUKA correction for protons (Baryon Twiki)
129 // (TPC tracks)
130 //
131         TF1* fnc = new TF1(name.Data(), "1-[0]*exp([1]*x)+[2]", 0.001, 10);
132         fnc->SetParameters(4113, -32.9, -0.009383);
133         Double_t err[] = {1640.4, 1.5, 0.000917};
134         fnc->SetParErrors(err);
135         
136         return fnc;
137 }
138
139 TF1* AliLnEfficiency::GetAntiProtonG3FlukaCor(const TString& name) const
140 {
141 //
142 // GEANT3/FLUKA correction for antiprotons (Baryon Twiki)
143 // (TPC tracks)
144 //
145         TF1* fnc = new TF1(name.Data(), "1-[0]*exp([1]*x)+[2]+[3]*log(x)/pow(x,0.2)", 0.001, 10);
146         fnc->SetParameters(177.4, -22.03, -0.06538, 0.04431 );
147         Double_t err[] = {73., 1.5, 0.00221, 0.00555};
148         fnc->SetParErrors(err);
149         
150         return fnc;
151 }
152