moving AliAnalysisTaskDiHadron (Jason) to PWGCF
[u/mrichter/AliRoot.git] / PWGDQ / base / AliQuarkoniaEfficiency.cxx
1
2 /* $Id$ */
3
4 //===================================================================
5 //  Class AliQuarkoniaEfficiency                               
6 //
7 //  This class will provide the quarkonia reconstruction efficiency 
8 //  in ALICE without acceptance consideration.
9 //
10 //
11 //  Reconstruction efficiency has been evaluated by means of a flat
12 //  y and pt distribution of quarkonia in -4 < y < -2.5, 
13 //  0 < pt < 20 GeV/c. Weights have been used to evaluate the
14 //  reconstruction efficiency in different parameterizations.
15 //
16 //  Beware that efficiency histos are preliminary.
17 //  Just Jpsi, Dimuon, UnlikePair trigger
18 //  efficiencies should be considered.
19 //
20 //
21 //  Example:
22 //   Double_t eff,error;
23 //   AliQuarkoniaEfficiency * JPsiEff = new AliQuarkoniaEfficiency();
24 //   JPsiEff->Init();
25 //   JPsiEff->GetEfficiency(Rapidity,Pt,eff,error);
26 //   printf(" eff = %2.2e  error %2.2e \n",eff,error);
27 //                                                          
28 //                                                              
29 //  Subatech 2006
30 //===================================================================
31
32 // Root 
33 #include "TFile.h"
34 #include "TH2.h"
35
36 // AliRoot includes
37 #include "AliQuarkoniaEfficiency.h"
38 #include "AliLog.h"
39
40
41 //_______________________________________________________________________
42 AliQuarkoniaEfficiency::AliQuarkoniaEfficiency(Int_t quarkoniaResonance, Int_t decayChannel,
43                                                Int_t simParameterization):
44   fEfficiencyFileName("QuarkoniaEfficiency.root"),
45   fQuarkoniaResonance(quarkoniaResonance),     
46   fDecayChannel(decayChannel),           
47   fParameterization(simParameterization),  
48   fTriggerType(kPairUnlikeApt),
49   fTrigger(kFALSE),
50   fEfficiency(0x0)
51 {
52   // Constructor
53 }
54
55 //_______________________________________________________________________
56 AliQuarkoniaEfficiency::~AliQuarkoniaEfficiency()
57 {
58   // Destructor
59   delete fEfficiency;
60 }
61
62 //_______________________________________________________________________
63 void AliQuarkoniaEfficiency::Init()
64 {
65   // Initialize method
66   switch (fQuarkoniaResonance)  {  
67   case kJpsi:
68     SetTitle("Jpsi");
69     break;
70   case kPsiP:
71     SetTitle("PsiP");
72     break;
73   case kUpsilon:
74     SetTitle("Upsilon");
75     break;
76   case kUpsilonP:
77     SetTitle("UpsilonP");
78     break;
79   case kUpsilonPP:
80     SetTitle("UpsilonPP");
81     break;
82   case kOmega:
83     SetTitle("Omega");
84     break;
85   case kPhi:
86     SetTitle("Phi");
87     break;
88   }
89
90   switch (fDecayChannel) {
91   case kDimuon:
92     SetName("Dimuon");
93     break;
94   case kDielectron:
95     SetName("Dielectron");
96     break;
97   }
98
99   const char *param=0;
100   switch (fParameterization){
101   case kFlat:
102     param = "Flat";
103     break;
104   case kCDFscaled:
105     param = "CDFscaled";
106     break;
107   case kCDFscaledPP:
108     param = "CDFscaledPP";
109     break;
110   }
111   
112   const char *trig=0;
113   switch (fTriggerType){
114   case kSinglePlusLpt:
115     trig = "SinglePlusLpt";
116     break;
117   case kSinglePlusHpt: 
118     trig = "SinglePlusHpt";
119     break;
120   case kSinglePlusApt:
121     trig = "SinglePlusApt";
122     break;    
123   case kSingleMinusLpt:
124     trig = "SingleMinusLpt";
125     break;    
126   case kSingleMinusHpt:
127     trig = "SingleMinusHpt";
128     break;    
129   case kSingleMinusApt:
130     trig = "SingleMinusApt";
131     break;    
132   case kSingleUndefLpt:
133     trig = "SingleUndefLpt";
134     break;    
135   case kSingleUndefHpt:
136     trig = "SingleUndefHpt";
137     break;    
138   case kSingleUndefApt:
139     trig = "SingleUndefApt";
140     break;    
141   case kPairUnlikeLpt:
142     trig = "PairUnlikeLpt";
143     break;    
144   case kPairUnlikeHpt:
145     trig = "PairUnlikeHpt";
146     break;    
147   case kPairUnlikeApt:
148     trig = "PairUnlikeApt";
149     break;    
150   case kPairLikeLpt:
151     trig = "PairLikeLpt";
152     break;    
153   case kPairLikeHpt:
154     trig = "PairLikeHpt";
155     break;    
156   case kPairLikeApt:
157     trig = "PairLikeApt";
158     break;    
159   }
160
161
162   if(!fEfficiency) delete fEfficiency; 
163
164   TFile efficiencyFile(fEfficiencyFileName);
165   if ( efficiencyFile.IsOpen() ) {
166
167     char quarkoniaDir[15];
168     snprintf(quarkoniaDir,15,"%s",GetTitle());
169     if (! efficiencyFile.cd(quarkoniaDir) ){
170       AliError(Form("Directory %s not found in file %s \n Efficiency data for quarkonia %s and channel %s not found ",
171                     quarkoniaDir,fEfficiencyFileName.Data(),GetTitle(),GetName() ));
172       return;
173     }
174     
175     char histosDir[30];
176     snprintf(histosDir,30,"%s/%s_%s_%s",quarkoniaDir,GetTitle(),GetName(),param);
177     if(! efficiencyFile.cd(histosDir) ){
178       AliError(Form("Subdirectory %s/%s not found in file %s \n Efficiency data for quarkonia %s and channel %s not found ",
179                     quarkoniaDir,histosDir,fEfficiencyFileName.Data(),GetTitle(),GetName() ));
180       return;
181     }
182
183     char histoname[50];
184     if(fTrigger) snprintf(histoname,50,"h%sEfficiencyPtRap_%s",GetTitle(),trig);
185     else snprintf(histoname,50,"h%sEfficiencyPtRap",GetTitle());
186     char histonameposition[99];
187     snprintf(histonameposition,99,"%s/%s",histosDir,histoname);
188     fEfficiency = (TH2F*)efficiencyFile.Get(histonameposition);
189
190     if ( !fEfficiency ) {
191       AliError(Form("Histo %s not found in file %s \n Efficiency data for quarkonia %s and channel %s not found",
192                     histoname, fEfficiencyFileName.Data(), GetTitle(), GetName() ));
193     }
194     else {
195       fEfficiency->SetDirectory(0);
196     }
197     efficiencyFile.Close();
198
199   }
200   else {
201     AliError(Form("File %s not found",fEfficiencyFileName.Data()));
202   }
203
204 }
205
206 //_______________________________________________________________________
207 TH2F*  AliQuarkoniaEfficiency::GetEfficiencyHisto() const
208 {
209   // Returns the efficiency histogram
210   if (fEfficiency) return fEfficiency;
211   else {
212     AliError(Form("Efficiency data for quarkonia %s and channel %s not found",GetTitle(),GetName()));
213     return 0x0;
214   }
215 }
216
217 //_______________________________________________________________________
218 void  AliQuarkoniaEfficiency::GetEfficiency(Float_t rap, Float_t pT, Double_t &eff, Double_t &error)
219 {
220   // Evaluates the efficiency for a given (y,Pt) of the quarkonia
221   Int_t binx=0;
222   Int_t biny=0;
223   
224   if (!fEfficiency) {
225     AliError(Form("Efficiency data for quarkonia %s and channel %s not found",GetTitle(),GetName()));
226   }
227   else {
228     if ( rap < (fEfficiency->GetXaxis())->GetXmin()  ||  
229          rap > (fEfficiency->GetXaxis())->GetXmax()  ||
230          pT  < (fEfficiency->GetYaxis())->GetXmin()  ||  
231          pT  > (fEfficiency->GetYaxis())->GetXmax()   ) {
232       AliInfo("Values out of range");
233       eff   = 0.;
234       error = 0.;
235     }
236     else  { 
237       binx  = fEfficiency->GetXaxis()->FindBin(rap);  
238       biny  = fEfficiency->GetYaxis()->FindBin(pT);
239       eff   = fEfficiency->GetBinContent(binx,biny);
240       error = fEfficiency->GetBinError(binx,biny);
241     }
242   } 
243 }