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