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