Bug fix in PID step (ChiaraZ)
[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
641e1e0c 11#include "AliAnalysisHadEt.h"
12#include "TMath.h"
13#include "TList.h"
14#include "TH1F.h"
15#include "TH2F.h"
16#include <iostream>
17#include "AliAnalysisEtCuts.h"
020abff0 18#include "AliMCEvent.h"
641e1e0c 19#include "AliVEvent.h"
020abff0 20#include "AliStack.h"
21#include "AliESDtrackCuts.h"
22#include "TDatabasePDG.h"
23#include "TParticle.h"
24#include "Rtypes.h"
641e1e0c 25
26using namespace std;
27
28ClassImp(AliAnalysisHadEt);
29
30
020abff0 31Int_t AliAnalysisHadEt::fgnumOfEtaBins = 46;
32Float_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};
33Int_t AliAnalysisHadEt::fgNumOfPtBins = 111;
34Float_t AliAnalysisHadEt::fgPtAxis[117]=
641e1e0c 35 {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,
36 0.2, .22, .24, .26, .28, 0.30, 0.32, .34, .36, .38, 0.40, .42, .44, .46, .48,
37 0.5, .52, .54, .56, .58, 0.60, 0.62, .64, .66, .68, 0.70, .72, .74, .76, .78,
38 .80, .82, .84, .86, .88, 0.90, 0.92, .94, .96, .98, 1.00,1.05, 1.1,1.15, 1.2,
39 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,
40 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,
41 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,
42 20.0,25.0,30.0,35.0,40.0, 45.0, 50.0};
43
44AliAnalysisHadEt::AliAnalysisHadEt() :
45 fHistogramNameSuffix("")
46 ,fPdgDB(0)
020abff0 47 ,fPiPlusCode(0)
48 ,fPiMinusCode(0)
49 ,fKPlusCode(0)
50 ,fKMinusCode(0)
51 ,fProtonCode(0)
52 ,fAntiProtonCode(0)
53 ,fLambdaCode(0)
54 ,fAntiLambdaCode(0)
55 ,fK0SCode(0)
56 ,fOmegaCode(0)
57 ,fAntiOmegaCode(0)
58 ,fXi0Code(0)
59 ,fAntiXi0Code(0)
60 ,fXiCode(0)
61 ,fAntiXiCode(0)
62 ,fSigmaCode(0)
63 ,fAntiSigmaCode(0)
64 ,fK0LCode(0)
65 ,fNeutronCode(0)
66 ,fAntiNeutronCode(0)
67 ,fEPlusCode(0)
68 ,fEMinusCode(0)
69 ,fPionMass(0)
641e1e0c 70 ,fSumEt(0)
71 ,fSumEtAcc(0)
72 ,fTotEt(0)
73 ,fTotEtAcc(0)
74 ,fTotNeutralEt(0)
75 ,fTotNeutralEtAcc(0)
76 ,fTotChargedEt(0)
77 ,fTotChargedEtAcc(0)
78 ,fMultiplicity(0)
79 ,fChargedMultiplicity(0)
80 ,fNeutralMultiplicity(0)
81 ,fEtaCut(EtCommonCuts::kEtaCut)
82 ,fEtaCutAcc(0)
83 //,fPhiCutAccMin(0)
84 //,fPhiCutAccMax(360.)
85 ,fVertexXCut(0)
86 ,fVertexYCut(0)
87 ,fVertexZCut(0)
88 ,fIPxyCut(0)
89 ,fIPzCut(0)
90 //,fSingleCellEnergyCut(0)
91 //,fClusterEnergyCut(EtCommonCuts::kClusterEnergyCut)
92 //,fTrackPtCut(EtCommonCuts::kTrackPtCut)
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.
641e1e0c 108 cout<<"This event has "<<event->GetNumberOfTracks()<<" tracks"<<endl;
020abff0 109 ResetEventValues();
110 AliMCEvent *mcEvent = dynamic_cast<AliMCEvent*>(event);
111
112 // Let's play with the stack!
113 AliStack *stack = mcEvent->Stack();
114
115 Int_t nPrim = stack->GetNtrack();
116
117 for (Int_t iPart = 0; iPart < nPrim; iPart++)
118 {
119
120 TParticle *part = stack->Particle(iPart);
121
122 if (!part)
123 {
124 Printf("ERROR: Could not get particle %d", iPart);
125 continue;
126 }
127 }
641e1e0c 128 return 0;
129}
130
131void AliAnalysisHadEt::FillOutputList()
020abff0 132{//fill the output histogram list with histograms in all AliAnalysisHadEt's
641e1e0c 133}
134
135void AliAnalysisHadEt::Init()
020abff0 136{//Initiate member vaiables to reasonable values
cf6522d1 137 fVertexXCut = EtReconstructedCuts::kVertexXCut;
138 fVertexYCut = EtReconstructedCuts::kVertexYCut;
139 fVertexZCut = EtReconstructedCuts::kVertexZCut;
140 fIPxyCut = EtReconstructedCuts::kIPxyCut;
141 fIPzCut = EtReconstructedCuts::kIPzCut;
142
641e1e0c 143 if(!fPdgDB) fPdgDB = new TDatabasePDG();
144 //the codes are defined in $ROOTSYS/etc/pdg_table.txt
020abff0 145 fPionMass = fPdgDB->GetParticle("pi+")->Mass();
146 fPiPlusCode = fPdgDB->GetParticle("pi+")->PdgCode();
147 fPiMinusCode = fPdgDB->GetParticle("pi-")->PdgCode();
148 fKPlusCode = fPdgDB->GetParticle("K+")->PdgCode();
149 fKMinusCode = fPdgDB->GetParticle("K-")->PdgCode();
150 fProtonCode = fPdgDB->GetParticle("proton")->PdgCode();
151 fAntiProtonCode = fPdgDB->GetParticle("antiproton")->PdgCode();
152 fLambdaCode = fPdgDB->GetParticle("Lambda0")->PdgCode();
153 fAntiLambdaCode = fPdgDB->GetParticle("Lambda0_bar")->PdgCode();
154 fK0SCode = fPdgDB->GetParticle("K_S0")->PdgCode();
155 fOmegaCode = fPdgDB->GetParticle("Omega-")->PdgCode();
156 fAntiOmegaCode = fPdgDB->GetParticle("Omega+")->PdgCode();
157 fXi0Code = fPdgDB->GetParticle("Xi0")->PdgCode();
158 fAntiXi0Code = fPdgDB->GetParticle("Xi0_bar")->PdgCode();
159 fXiCode = fPdgDB->GetParticle("Xi-")->PdgCode();
160 fAntiXiCode = fPdgDB->GetParticle("Xi-_bar")->PdgCode();
161 fSigmaCode = fPdgDB->GetParticle("Sigma-")->PdgCode();
162 fAntiSigmaCode = fPdgDB->GetParticle("Sigma+")->PdgCode();
163 fK0LCode = fPdgDB->GetParticle("K_L0")->PdgCode();
164 fNeutronCode = fPdgDB->GetParticle("neutron")->PdgCode();
165 fAntiNeutronCode = fPdgDB->GetParticle("antineutron")->PdgCode();
166 fEPlusCode = fPdgDB->GetParticle("e+")->PdgCode();
167 fEMinusCode = fPdgDB->GetParticle("e-")->PdgCode();
641e1e0c 168}
169
170void AliAnalysisHadEt::CreateHistograms()
020abff0 171{//creates histograms included in all AliAnalysisHadEt's
641e1e0c 172}
173
174void AliAnalysisHadEt::FillHistograms()
020abff0 175{//Fills histograms filled for all AliAnalysisHadEt's
641e1e0c 176}
177
178void AliAnalysisHadEt::ResetEventValues()
020abff0 179{//Resets event values of et to zero
641e1e0c 180 fTotEt = 0;
181 fTotEtAcc = 0;
182 fTotNeutralEt = 0;
183 fTotNeutralEtAcc = 0;
184 fTotChargedEt = 0;
185 fTotChargedEtAcc = 0;
186 fMultiplicity = 0;
187 fChargedMultiplicity = 0;
188 fNeutralMultiplicity = 0;
189}
190void AliAnalysisHadEt::CreateEtaPtHisto2D(TString name, TString title)
020abff0 191{ //creates a 2-d histogram in eta and phi and adds it to the list of histograms to be saved
641e1e0c 192 TString *histoname = new TString();
193 TString *histotitle = new TString();
194
195 histoname->Append(name);
196 histotitle->Append(title);
020abff0 197 //TH2F *h1 = new TH2F("h1", "Histogram with Gaussian random distribution", fgNumOfPtBins, ptBinsArray, fgnumOfEtaBins, etaBinsArray);
641e1e0c 198
020abff0 199 TH2F *histo = new TH2F(histoname->Data(),histotitle->Data(),fgNumOfPtBins, fgPtAxis, fgnumOfEtaBins, fgEtaAxis);
641e1e0c 200 histo->SetYTitle("#eta");
201 histo->SetXTitle("p_{T}");
202 histo->SetZTitle("E_{T}");
203 histo->Sumw2();
020abff0 204 fhistoList->Add(histo);
641e1e0c 205 delete histoname;
206 delete histotitle;
207
208}
209
210void AliAnalysisHadEt::CreateHisto1D(TString name, TString title, TString xtitle, TString ytitle,Int_t xbins, Float_t xlow,Float_t xhigh)
020abff0 211{ //creates a 1d histogram of the given dimensions and adds it to the list of histograms to be saved
641e1e0c 212 TString *histoname = new TString();
213 TString *histotitle = new TString();
214
215 //cout<<"creating "<<name<<endl;
216
217 histoname->Append(name);
218 histotitle->Append(title);
219 // printf("%s \n ",histoname->Data());
220 TH1F *histo = new TH1F(histoname->Data(),histotitle->Data(),xbins,xlow,xhigh);
221 histo->SetYTitle(ytitle);
222 histo->SetXTitle(xtitle);
223 histo->Sumw2();
020abff0 224 fhistoList->Add(histo);
641e1e0c 225 delete histoname;
226 delete histotitle;
227
228}
229void AliAnalysisHadEt::CreateIntHisto1D(TString name, TString title, TString xtitle, TString ytitle,Int_t xbins, Int_t xlow,Int_t xhigh)
020abff0 230{ //creates a 1d integer histogram and adds it to the list of histograms to be saved
641e1e0c 231 TString *histoname = new TString();
232 TString *histotitle = new TString();
233
234 //cout<<"creating "<<name<<endl;
235
236 histoname->Append(name);
237 histotitle->Append(title);
238 // printf("%s \n ",histoname->Data());
239 TH1I *histo = new TH1I(histoname->Data(),histotitle->Data(),xbins,xlow,xhigh);
240 histo->SetYTitle(ytitle);
241 histo->SetXTitle(xtitle);
242 histo->Sumw2();
020abff0 243 fhistoList->Add(histo);
641e1e0c 244 delete histoname;
245 delete histotitle;
246
247}
248void 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 249{ //creates a 2d histogram and adds it to the list of histograms to be saved
641e1e0c 250 TString *histoname = new TString();
251 TString *histotitle = new TString();
252
253 //cout<<"creating "<<name<<endl;
254
255 histoname->Append(name);
256 histotitle->Append(title);
257 // printf("%s \n ",histoname->Data());
258 TH2F *histo = new TH2F(histoname->Data(),histotitle->Data(),xbins,xlow,xhigh,ybins,ylow,yhigh);
259 histo->SetYTitle(ytitle);
260 histo->SetXTitle(xtitle);
261 histo->Sumw2();
020abff0 262 fhistoList->Add(histo);
641e1e0c 263 delete histoname;
264 delete histotitle;
265
266}
267void 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 268{ //creates a 2-d integer histogram and adds it to the list of histograms to be saved
641e1e0c 269 TString *histoname = new TString();
270 TString *histotitle = new TString();
271
272 //cout<<"creating "<<name<<endl;
273
274 histoname->Append(name);
275 histotitle->Append(title);
276 // printf("%s \n ",histoname->Data());
277 TH2I *histo = new TH2I(histoname->Data(),histotitle->Data(),xbins,xlow,xhigh,ybins,ylow,yhigh);
278 histo->SetYTitle(ytitle);
279 histo->SetXTitle(xtitle);
280 histo->Sumw2();
020abff0 281 fhistoList->Add(histo);
641e1e0c 282 delete histoname;
283 delete histotitle;
284
285}
286
287void AliAnalysisHadEt::CreateEtaHisto1D(TString name, TString title)
020abff0 288{ //creates 1d histogram in eta and adds it to the list of histograms to be saved
641e1e0c 289 TString *histoname = new TString();
290 TString *histotitle = new TString();
291
292
293 histoname->Append(name);
294 histotitle->Append(title);
020abff0 295 TH1F *histo = new TH1F(histoname->Data(),histotitle->Data(),fgnumOfEtaBins, fgEtaAxis);
641e1e0c 296 histo->SetYTitle("E_{T}");
297 histo->SetXTitle("#eta");
298 histo->Sumw2();
020abff0 299 fhistoList->Add(histo);
641e1e0c 300 delete histoname;
301 delete histotitle;
302
303}
304void AliAnalysisHadEt::FillHisto1D(TString histname, Float_t x, Float_t weight)
020abff0 305{//fills a 1d histogram with the name histoname with the value x and the weight "weight"
641e1e0c 306 TH1F *histo;
307 TString *name = new TString();
308
309 name->Append(histname);
020abff0 310 histo = (TH1F *)fhistoList->FindObject(name->Data());
641e1e0c 311 if(histo){
312 histo->Fill((Double_t)x, weight);
313 }
314 else{cerr<<"CorrelationMaker::FillHisto1D: no histogram "<< name->Data()<<endl;}
315 delete name;
316}
317void AliAnalysisHadEt::FillHisto2D(TString histname, Float_t x, Float_t y, Float_t weight)
020abff0 318{//fills a 2d histogram with the name histoname with the value x and the weight "weight"
641e1e0c 319 TH2F *histo;
320 TString *name = new TString();
321
322 name->Append(histname);
020abff0 323 histo = (TH2F *)fhistoList->FindObject(name->Data());
641e1e0c 324 if(histo){
325 histo->Fill((Double_t)x,(Double_t)y, weight);
326 }
327 else{cerr<<"CorrelationMaker::FillHisto2D: no histogram "<< name->Data()<<endl;}
328 delete name;
329}
330
331
020abff0 332Float_t AliAnalysisHadEt::Et(TParticle *part, float mass){//function to calculate et in the same way as it would be calculated in a calorimeter
333 if(mass+1000<0.01){//if no mass given return default. The default argument is -1000
641e1e0c 334 if(TMath::Abs(part->GetPDG(0)->PdgCode())==2212 || TMath::Abs(part->GetPDG(0)->PdgCode())==2112){
335 if(part->GetPDG(0)->PdgCode()==-2212 || part->GetPDG(0)->PdgCode()==-2112){//antiproton or antineutron
336 //for antinucleons we specifically want to return the kinetic energy plus twice the rest mass
337 return (part->Energy()+part->GetMass())*TMath::Sin(part->Theta());
338 }
339 if(part->GetPDG(0)->PdgCode()==2212 || part->GetPDG(0)->PdgCode()==2112){//antiproton or antineutron
340 //for nucleons we specifically want to return the kinetic energy only
341 return (part->Energy()-part->GetMass())*TMath::Sin(part->Theta());
342 }
343 }
344 else{//otherwise go to the default
345 return part->Energy()*TMath::Sin(part->Theta());
346 }
347 }
348 else{//otherwise use the mass that was given
349 return (TMath::Sqrt(TMath::Power(part->P(),2.0)+TMath::Power(mass,2.0)))*TMath::Sin(part->Theta());
350 }
351 return 0.0;
352}