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