Corrections implemented into AliAnalysisHadEtReconstructed, simulated hadronic and...
[u/mrichter/AliRoot.git] / PWG4 / totEt / AliAnalysisHadEt.cxx
CommitLineData
cf6522d1 1//Created by Christine Nattrass, Rebecca Scott, Irakli Martashvili
641e1e0c 2//University of Tennessee at Knoxville
cf6522d1 3//
4// This class is designed for the analysis of the hadronic component of
5// transverse energy. It is used by AliAnalysisTaskHadEt.
6// This gets information about the hadronic component of the transverse energy
7// from tracks reconstructed in an event
8// it has daughters, AliAnalysisHadEtMonteCarlo and
9// AliAnalysisHadEtReconstructed which loop over either Monte Carlo data or
10// real data to get Et
3ce6b879 11
641e1e0c 12#include "AliAnalysisHadEt.h"
13#include "TMath.h"
14#include "TList.h"
15#include "TH1F.h"
16#include "TH2F.h"
17#include <iostream>
18#include "AliAnalysisEtCuts.h"
020abff0 19#include "AliMCEvent.h"
641e1e0c 20#include "AliVEvent.h"
020abff0 21#include "AliStack.h"
22#include "AliESDtrackCuts.h"
23#include "TDatabasePDG.h"
24#include "TParticle.h"
25#include "Rtypes.h"
641e1e0c 26
27using namespace std;
28
29ClassImp(AliAnalysisHadEt);
30
31
020abff0 32Int_t AliAnalysisHadEt::fgnumOfEtaBins = 46;
33Float_t AliAnalysisHadEt::fgEtaAxis[47]={-0.78, -0.74, -0.7, -0.66, -0.62, -0.58, -0.54, -0.5, -0.46, -0.42, -0.38, -0.34, -0.3, -0.26, -0.22, -0.18, -0.14, -0.12, -0.1, -0.08, -0.06, -0.04, -0.02, -0.0, 0.02, 0.04, 0.06, 0.08, 0.1, 0.12, 0.14, 0.18, 0.22, 0.26, 0.3, 0.34, 0.38, 0.42, 0.46, 0.5, 0.54, 0.58, 0.62, 0.66, 0.7, 0.74, 0.78};
34Int_t AliAnalysisHadEt::fgNumOfPtBins = 111;
35Float_t AliAnalysisHadEt::fgPtAxis[117]=
641e1e0c 36 {0.0,0.01,0.02,0.03,0.04, 0.05, 0.06,0.07,0.08,0.09, 0.10,0.11, .12,0.13, .14,0.15, .16,0.17, .18,0.19,
37 0.2, .22, .24, .26, .28, 0.30, 0.32, .34, .36, .38, 0.40, .42, .44, .46, .48,
38 0.5, .52, .54, .56, .58, 0.60, 0.62, .64, .66, .68, 0.70, .72, .74, .76, .78,
39 .80, .82, .84, .86, .88, 0.90, 0.92, .94, .96, .98, 1.00,1.05, 1.1,1.15, 1.2,
40 1.25, 1.3,1.35,1.40,1.45, 1.50, 1.55, 1.6,1.65, 1.7, 1.75, 1.8,1.85, 1.9,1.95,
41 2.0, 2.2, 2.4, 2.6, 2.8, 3.00, 3.20, 3.4, 3.6, 3.8, 4.00, 4.2, 4.4, 4.6, 4.8,
42 5.0, 5.5, 6.0, 6.5, 7.0, 7.50, 8.00, 8.5, 9.0, 9.5, 10.0,12.0,14.0,16.0,18.0,
43 20.0,25.0,30.0,35.0,40.0, 45.0, 50.0};
44
d6214a64 45Float_t AliAnalysisHadEt::fgPtTPCCutOff = 0.15;
46Float_t AliAnalysisHadEt::fgPtITSCutOff = 0.10;
47
641e1e0c 48AliAnalysisHadEt::AliAnalysisHadEt() :
49 fHistogramNameSuffix("")
4998becf 50 ,fCuts(0)
641e1e0c 51 ,fPdgDB(0)
020abff0 52 ,fPiPlusCode(0)
53 ,fPiMinusCode(0)
54 ,fKPlusCode(0)
55 ,fKMinusCode(0)
56 ,fProtonCode(0)
57 ,fAntiProtonCode(0)
58 ,fLambdaCode(0)
59 ,fAntiLambdaCode(0)
60 ,fK0SCode(0)
61 ,fOmegaCode(0)
62 ,fAntiOmegaCode(0)
63 ,fXi0Code(0)
64 ,fAntiXi0Code(0)
65 ,fXiCode(0)
66 ,fAntiXiCode(0)
67 ,fSigmaCode(0)
68 ,fAntiSigmaCode(0)
69 ,fK0LCode(0)
70 ,fNeutronCode(0)
71 ,fAntiNeutronCode(0)
72 ,fEPlusCode(0)
73 ,fEMinusCode(0)
f427cbed 74 ,fGammaCode(0)
75 ,fPi0Code(0)
76 ,fEtaCode(0)
b64de20c 77 ,fOmega0Code(0)
020abff0 78 ,fPionMass(0)
3ce6b879 79 ,fKaonMass(0)
80 ,fProtonMass(0)
81 ,fElectronMass(0)
641e1e0c 82 ,fSumEt(0)
83 ,fSumEtAcc(0)
84 ,fTotEt(0)
85 ,fTotEtAcc(0)
86 ,fTotNeutralEt(0)
87 ,fTotNeutralEtAcc(0)
88 ,fTotChargedEt(0)
89 ,fTotChargedEtAcc(0)
90 ,fMultiplicity(0)
91 ,fChargedMultiplicity(0)
92 ,fNeutralMultiplicity(0)
cf6522d1 93 ,fEsdtrackCutsITSTPC(0)
94 ,fEsdtrackCutsTPC(0)
95 ,fEsdtrackCutsITS(0)
020abff0 96 ,fhistoList(0)
97{//default constructor
641e1e0c 98
99}
100
101AliAnalysisHadEt::~AliAnalysisHadEt()
020abff0 102{//destructor
641e1e0c 103
104}
105
106Int_t AliAnalysisHadEt::AnalyseEvent(AliVEvent *event)
020abff0 107{ //this line is basically here to eliminate a compiler warning that event is not used. Making it a virtual function did not work with the plugin.
0065531a 108 cout << "This event has " << event->GetNumberOfTracks() << " tracks" << endl;
109 ResetEventValues();
641e1e0c 110 return 0;
111}
112
113void AliAnalysisHadEt::FillOutputList()
020abff0 114{//fill the output histogram list with histograms in all AliAnalysisHadEt's
641e1e0c 115}
116
117void AliAnalysisHadEt::Init()
0065531a 118{// clear variables, set up cuts and PDG info
119 ResetEventValues();
3ce6b879 120
641e1e0c 121}
122
123void AliAnalysisHadEt::CreateHistograms()
020abff0 124{//creates histograms included in all AliAnalysisHadEt's
641e1e0c 125}
126
127void AliAnalysisHadEt::FillHistograms()
020abff0 128{//Fills histograms filled for all AliAnalysisHadEt's
641e1e0c 129}
130
131void AliAnalysisHadEt::ResetEventValues()
020abff0 132{//Resets event values of et to zero
0065531a 133 fTotEt = 0;
134 fTotEtAcc = 0;
135 fTotNeutralEt = 0;
136 fTotNeutralEtAcc = 0;
137 fTotChargedEt = 0;
138 fTotChargedEtAcc = 0;
139 fMultiplicity = 0;
140 fChargedMultiplicity = 0;
141 fNeutralMultiplicity = 0;
142
143 if (!fCuts || !fPdgDB || fPiPlusCode==0) { // some Init's needed
144 cout << __FILE__ << ":" << __LINE__ << " : Init " << endl;
145 if (!fCuts) {
146 cout << " setting up Cuts " << endl;
147 fCuts = new AliAnalysisEtCuts();
148 }
149 if(!fPdgDB) {
150 cout << " setting up PdgDB " << endl;
151 fPdgDB = new TDatabasePDG();
152 }
153
154 if (fPiPlusCode==0) {
155 SetParticleCodes();
156 }
157 }
641e1e0c 158}
4998becf 159
160void AliAnalysisHadEt::SetParticleCodes()
161{ //the codes are defined in $ROOTSYS/etc/pdg_table.txt
f427cbed 162 fPionMass = fPdgDB->GetParticle("pi+")->Mass();
3ce6b879 163 fKaonMass = fPdgDB->GetParticle("K+")->Mass();
164 fProtonMass = fPdgDB->GetParticle("proton")->Mass();
165 fElectronMass = fPdgDB->GetParticle("e+")->Mass();
f427cbed 166 fPiPlusCode = fPdgDB->GetParticle("pi+")->PdgCode();
167 fPiMinusCode = fPdgDB->GetParticle("pi-")->PdgCode();
168 fKPlusCode = fPdgDB->GetParticle("K+")->PdgCode();
169 fKMinusCode = fPdgDB->GetParticle("K-")->PdgCode();
170 fProtonCode = fPdgDB->GetParticle("proton")->PdgCode();
171 fAntiProtonCode = fPdgDB->GetParticle("antiproton")->PdgCode();
172 fLambdaCode = fPdgDB->GetParticle("Lambda0")->PdgCode();
173 fAntiLambdaCode = fPdgDB->GetParticle("Lambda0_bar")->PdgCode();
174 fK0SCode = fPdgDB->GetParticle("K_S0")->PdgCode();
175 fOmegaCode = fPdgDB->GetParticle("Omega-")->PdgCode();
176 fAntiOmegaCode = fPdgDB->GetParticle("Omega+")->PdgCode();
177 fXi0Code = fPdgDB->GetParticle("Xi0")->PdgCode();
178 fAntiXi0Code = fPdgDB->GetParticle("Xi0_bar")->PdgCode();
179 fXiCode = fPdgDB->GetParticle("Xi-")->PdgCode();
180 fAntiXiCode = fPdgDB->GetParticle("Xi-_bar")->PdgCode();
181 fSigmaCode = fPdgDB->GetParticle("Sigma-")->PdgCode();
182 fAntiSigmaCode = fPdgDB->GetParticle("Sigma+")->PdgCode();
183 fK0LCode = fPdgDB->GetParticle("K_L0")->PdgCode();
184 fNeutronCode = fPdgDB->GetParticle("neutron")->PdgCode();
185 fAntiNeutronCode = fPdgDB->GetParticle("antineutron")->PdgCode();
186 fEPlusCode = fPdgDB->GetParticle("e+")->PdgCode();
187 fEMinusCode = fPdgDB->GetParticle("e-")->PdgCode();
188 fGammaCode = fPdgDB->GetParticle("gamma")->PdgCode();
189 fPi0Code = fPdgDB->GetParticle("pi0")->PdgCode();
190 fEtaCode = fPdgDB->GetParticle("eta")->PdgCode();
b64de20c 191 fOmega0Code = fPdgDB->GetParticle("omega")->PdgCode();
f427cbed 192 cout << "Resetting Codes: Pion " << fPiPlusCode
193 << "," << fPiMinusCode
194 << " Kaon " << fKPlusCode
195 << "," << fKMinusCode
196 << " Gamma "<<fGammaCode
197 << " Pi0 "<<fPi0Code
198 << " Eta "<<fEtaCode
b64de20c 199 << " omega "<<fOmega0Code
f427cbed 200 << endl;
4998becf 201}
202
641e1e0c 203void AliAnalysisHadEt::CreateEtaPtHisto2D(TString name, TString title)
020abff0 204{ //creates a 2-d histogram in eta and phi and adds it to the list of histograms to be saved
641e1e0c 205 TString *histoname = new TString();
206 TString *histotitle = new TString();
207
208 histoname->Append(name);
209 histotitle->Append(title);
020abff0 210 //TH2F *h1 = new TH2F("h1", "Histogram with Gaussian random distribution", fgNumOfPtBins, ptBinsArray, fgnumOfEtaBins, etaBinsArray);
641e1e0c 211
020abff0 212 TH2F *histo = new TH2F(histoname->Data(),histotitle->Data(),fgNumOfPtBins, fgPtAxis, fgnumOfEtaBins, fgEtaAxis);
641e1e0c 213 histo->SetYTitle("#eta");
214 histo->SetXTitle("p_{T}");
215 histo->SetZTitle("E_{T}");
216 histo->Sumw2();
020abff0 217 fhistoList->Add(histo);
641e1e0c 218 delete histoname;
219 delete histotitle;
220
221}
222
223void AliAnalysisHadEt::CreateHisto1D(TString name, TString title, TString xtitle, TString ytitle,Int_t xbins, Float_t xlow,Float_t xhigh)
020abff0 224{ //creates a 1d histogram of the given dimensions and adds it to the list of histograms to be saved
641e1e0c 225 TString *histoname = new TString();
226 TString *histotitle = new TString();
227
228 //cout<<"creating "<<name<<endl;
229
230 histoname->Append(name);
231 histotitle->Append(title);
232 // printf("%s \n ",histoname->Data());
233 TH1F *histo = new TH1F(histoname->Data(),histotitle->Data(),xbins,xlow,xhigh);
234 histo->SetYTitle(ytitle);
235 histo->SetXTitle(xtitle);
236 histo->Sumw2();
020abff0 237 fhistoList->Add(histo);
641e1e0c 238 delete histoname;
239 delete histotitle;
240
241}
242void AliAnalysisHadEt::CreateIntHisto1D(TString name, TString title, TString xtitle, TString ytitle,Int_t xbins, Int_t xlow,Int_t xhigh)
020abff0 243{ //creates a 1d integer histogram and adds it to the list of histograms to be saved
641e1e0c 244 TString *histoname = new TString();
245 TString *histotitle = new TString();
246
247 //cout<<"creating "<<name<<endl;
248
249 histoname->Append(name);
250 histotitle->Append(title);
251 // printf("%s \n ",histoname->Data());
252 TH1I *histo = new TH1I(histoname->Data(),histotitle->Data(),xbins,xlow,xhigh);
253 histo->SetYTitle(ytitle);
254 histo->SetXTitle(xtitle);
255 histo->Sumw2();
020abff0 256 fhistoList->Add(histo);
641e1e0c 257 delete histoname;
258 delete histotitle;
259
260}
261void AliAnalysisHadEt::CreateHisto2D(TString name, TString title, TString xtitle, TString ytitle,Int_t xbins, Float_t xlow,Float_t xhigh,Int_t ybins,Float_t ylow,Float_t yhigh)
020abff0 262{ //creates a 2d histogram and adds it to the list of histograms to be saved
641e1e0c 263 TString *histoname = new TString();
264 TString *histotitle = new TString();
265
266 //cout<<"creating "<<name<<endl;
267
268 histoname->Append(name);
269 histotitle->Append(title);
270 // printf("%s \n ",histoname->Data());
271 TH2F *histo = new TH2F(histoname->Data(),histotitle->Data(),xbins,xlow,xhigh,ybins,ylow,yhigh);
272 histo->SetYTitle(ytitle);
273 histo->SetXTitle(xtitle);
274 histo->Sumw2();
020abff0 275 fhistoList->Add(histo);
641e1e0c 276 delete histoname;
277 delete histotitle;
278
279}
280void AliAnalysisHadEt::CreateIntHisto2D(TString name, TString title, TString xtitle, TString ytitle,Int_t xbins, Int_t xlow,Int_t xhigh,Int_t ybins,Int_t ylow,Int_t yhigh)
020abff0 281{ //creates a 2-d integer histogram and adds it to the list of histograms to be saved
641e1e0c 282 TString *histoname = new TString();
283 TString *histotitle = new TString();
284
285 //cout<<"creating "<<name<<endl;
286
287 histoname->Append(name);
288 histotitle->Append(title);
289 // printf("%s \n ",histoname->Data());
290 TH2I *histo = new TH2I(histoname->Data(),histotitle->Data(),xbins,xlow,xhigh,ybins,ylow,yhigh);
291 histo->SetYTitle(ytitle);
292 histo->SetXTitle(xtitle);
293 histo->Sumw2();
020abff0 294 fhistoList->Add(histo);
641e1e0c 295 delete histoname;
296 delete histotitle;
297
298}
299
300void AliAnalysisHadEt::CreateEtaHisto1D(TString name, TString title)
020abff0 301{ //creates 1d histogram in eta and adds it to the list of histograms to be saved
641e1e0c 302 TString *histoname = new TString();
303 TString *histotitle = new TString();
304
305
306 histoname->Append(name);
307 histotitle->Append(title);
020abff0 308 TH1F *histo = new TH1F(histoname->Data(),histotitle->Data(),fgnumOfEtaBins, fgEtaAxis);
641e1e0c 309 histo->SetYTitle("E_{T}");
310 histo->SetXTitle("#eta");
311 histo->Sumw2();
020abff0 312 fhistoList->Add(histo);
641e1e0c 313 delete histoname;
314 delete histotitle;
315
316}
317void AliAnalysisHadEt::FillHisto1D(TString histname, Float_t x, Float_t weight)
020abff0 318{//fills a 1d histogram with the name histoname with the value x and the weight "weight"
641e1e0c 319 TH1F *histo;
320 TString *name = new TString();
321
322 name->Append(histname);
020abff0 323 histo = (TH1F *)fhistoList->FindObject(name->Data());
641e1e0c 324 if(histo){
325 histo->Fill((Double_t)x, weight);
326 }
327 else{cerr<<"CorrelationMaker::FillHisto1D: no histogram "<< name->Data()<<endl;}
328 delete name;
329}
330void AliAnalysisHadEt::FillHisto2D(TString histname, Float_t x, Float_t y, Float_t weight)
020abff0 331{//fills a 2d histogram with the name histoname with the value x and the weight "weight"
641e1e0c 332 TH2F *histo;
333 TString *name = new TString();
334
335 name->Append(histname);
020abff0 336 histo = (TH2F *)fhistoList->FindObject(name->Data());
641e1e0c 337 if(histo){
338 histo->Fill((Double_t)x,(Double_t)y, weight);
339 }
340 else{cerr<<"CorrelationMaker::FillHisto2D: no histogram "<< name->Data()<<endl;}
341 delete name;
342}
343
344
020abff0 345Float_t AliAnalysisHadEt::Et(TParticle *part, float mass){//function to calculate et in the same way as it would be calculated in a calorimeter
346 if(mass+1000<0.01){//if no mass given return default. The default argument is -1000
641e1e0c 347 if(TMath::Abs(part->GetPDG(0)->PdgCode())==2212 || TMath::Abs(part->GetPDG(0)->PdgCode())==2112){
348 if(part->GetPDG(0)->PdgCode()==-2212 || part->GetPDG(0)->PdgCode()==-2112){//antiproton or antineutron
349 //for antinucleons we specifically want to return the kinetic energy plus twice the rest mass
350 return (part->Energy()+part->GetMass())*TMath::Sin(part->Theta());
351 }
352 if(part->GetPDG(0)->PdgCode()==2212 || part->GetPDG(0)->PdgCode()==2112){//antiproton or antineutron
353 //for nucleons we specifically want to return the kinetic energy only
354 return (part->Energy()-part->GetMass())*TMath::Sin(part->Theta());
355 }
356 }
357 else{//otherwise go to the default
358 return part->Energy()*TMath::Sin(part->Theta());
359 }
360 }
361 else{//otherwise use the mass that was given
362 return (TMath::Sqrt(TMath::Power(part->P(),2.0)+TMath::Power(mass,2.0)))*TMath::Sin(part->Theta());
363 }
364 return 0.0;
365}
3ce6b879 366Float_t AliAnalysisHadEt::Et(Float_t p, Float_t theta, Int_t pid, Short_t charge){//function to calculate et in the same way as it would be calculated in a calorimeter
367 if(pid==fPiPlusCode || pid==fPiMinusCode){//Nothing special for pions
368 return TMath::Sqrt(p*p + fPionMass*fPionMass) * TMath::Sin(theta);
369 }
370 if(pid==fKPlusCode || pid==fKMinusCode){//Nothing special for kaons
371 return TMath::Sqrt(p*p + fKaonMass*fKaonMass) * TMath::Sin(theta);
372 }
373 if(pid==fEPlusCode || pid==fEMinusCode){//Nothing special for electrons
374 return TMath::Sqrt(p*p + fElectronMass*fElectronMass) * TMath::Sin(theta);
375 }
376 if(pid==fProtonCode || pid==fAntiProtonCode){//But for protons we must be careful...
377 if(charge<0.0){//antiprotns: kinetic energy plus twice the rest mass
378 return (TMath::Sqrt(p*p + fProtonMass*fProtonMass) + fProtonMass) * TMath::Sin(theta);
379 }
380 if(charge>0.0){//antiprotns: kinetic energy only
381 return (TMath::Sqrt(p*p + fProtonMass*fProtonMass) - fProtonMass) * TMath::Sin(theta);
382 }
383 }
384 cerr<<"Uh-oh! Et not set properly!"<<endl;
385 return 0.0;
386}