Adding phenix points to p/pi plot and implementing a method in AliBWTools to divide...
[u/mrichter/AliRoot.git] / PWG4 / totEt / AliAnalysisHadEt.cxx
CommitLineData
641e1e0c 1//Create by Christine Nattrass, Rebecca Scott, Irakli Martashvili
2//University of Tennessee at Knoxville
020abff0 3//This class is designed for the analysis of the hadronic component of transverse energy. It is used by AliAnalysisTaskHadEt.
4//This gets information about the hadronic component of the transverse energy from tracks reconstructed in an event
5//it has daughters, AliAnalysisHadEtMonteCarlo and AliAnalysisHadEtReconstructed, which loop over either Monte Carlo data or real data to get Et
641e1e0c 6#include "AliAnalysisHadEt.h"
7#include "TMath.h"
8#include "TList.h"
9#include "TH1F.h"
10#include "TH2F.h"
11#include <iostream>
12#include "AliAnalysisEtCuts.h"
020abff0 13#include "AliMCEvent.h"
641e1e0c 14#include "AliVEvent.h"
020abff0 15#include "AliStack.h"
16#include "AliESDtrackCuts.h"
17#include "TDatabasePDG.h"
18#include "TParticle.h"
19#include "Rtypes.h"
641e1e0c 20
21using namespace std;
22
23ClassImp(AliAnalysisHadEt);
24
25
020abff0 26Int_t AliAnalysisHadEt::fgnumOfEtaBins = 46;
27Float_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};
28Int_t AliAnalysisHadEt::fgNumOfPtBins = 111;
29Float_t AliAnalysisHadEt::fgPtAxis[117]=
641e1e0c 30 {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,
31 0.2, .22, .24, .26, .28, 0.30, 0.32, .34, .36, .38, 0.40, .42, .44, .46, .48,
32 0.5, .52, .54, .56, .58, 0.60, 0.62, .64, .66, .68, 0.70, .72, .74, .76, .78,
33 .80, .82, .84, .86, .88, 0.90, 0.92, .94, .96, .98, 1.00,1.05, 1.1,1.15, 1.2,
34 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,
35 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,
36 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,
37 20.0,25.0,30.0,35.0,40.0, 45.0, 50.0};
38
39AliAnalysisHadEt::AliAnalysisHadEt() :
40 fHistogramNameSuffix("")
41 ,fPdgDB(0)
020abff0 42 ,fPiPlusCode(0)
43 ,fPiMinusCode(0)
44 ,fKPlusCode(0)
45 ,fKMinusCode(0)
46 ,fProtonCode(0)
47 ,fAntiProtonCode(0)
48 ,fLambdaCode(0)
49 ,fAntiLambdaCode(0)
50 ,fK0SCode(0)
51 ,fOmegaCode(0)
52 ,fAntiOmegaCode(0)
53 ,fXi0Code(0)
54 ,fAntiXi0Code(0)
55 ,fXiCode(0)
56 ,fAntiXiCode(0)
57 ,fSigmaCode(0)
58 ,fAntiSigmaCode(0)
59 ,fK0LCode(0)
60 ,fNeutronCode(0)
61 ,fAntiNeutronCode(0)
62 ,fEPlusCode(0)
63 ,fEMinusCode(0)
64 ,fPionMass(0)
641e1e0c 65 ,fSumEt(0)
66 ,fSumEtAcc(0)
67 ,fTotEt(0)
68 ,fTotEtAcc(0)
69 ,fTotNeutralEt(0)
70 ,fTotNeutralEtAcc(0)
71 ,fTotChargedEt(0)
72 ,fTotChargedEtAcc(0)
73 ,fMultiplicity(0)
74 ,fChargedMultiplicity(0)
75 ,fNeutralMultiplicity(0)
76 ,fEtaCut(EtCommonCuts::kEtaCut)
77 ,fEtaCutAcc(0)
78 //,fPhiCutAccMin(0)
79 //,fPhiCutAccMax(360.)
80 ,fVertexXCut(0)
81 ,fVertexYCut(0)
82 ,fVertexZCut(0)
83 ,fIPxyCut(0)
84 ,fIPzCut(0)
85 //,fSingleCellEnergyCut(0)
86 //,fClusterEnergyCut(EtCommonCuts::kClusterEnergyCut)
87 //,fTrackPtCut(EtCommonCuts::kTrackPtCut)
020abff0 88 ,ffesdtrackCutsITSTPC(0)
89 ,fesdtrackCutsTPC(0)
90 ,fesdtrackCutsITS(0)
91 ,fhistoList(0)
92{//default constructor
641e1e0c 93
94}
95
96AliAnalysisHadEt::~AliAnalysisHadEt()
020abff0 97{//destructor
641e1e0c 98
99}
100
101Int_t AliAnalysisHadEt::AnalyseEvent(AliVEvent *event)
020abff0 102{ //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 103 cout<<"This event has "<<event->GetNumberOfTracks()<<" tracks"<<endl;
020abff0 104 ResetEventValues();
105 AliMCEvent *mcEvent = dynamic_cast<AliMCEvent*>(event);
106
107 // Let's play with the stack!
108 AliStack *stack = mcEvent->Stack();
109
110 Int_t nPrim = stack->GetNtrack();
111
112 for (Int_t iPart = 0; iPart < nPrim; iPart++)
113 {
114
115 TParticle *part = stack->Particle(iPart);
116
117 if (!part)
118 {
119 Printf("ERROR: Could not get particle %d", iPart);
120 continue;
121 }
122 }
641e1e0c 123 return 0;
124}
125
126void AliAnalysisHadEt::FillOutputList()
020abff0 127{//fill the output histogram list with histograms in all AliAnalysisHadEt's
641e1e0c 128}
129
130void AliAnalysisHadEt::Init()
020abff0 131{//Initiate member vaiables to reasonable values
641e1e0c 132 if(!fPdgDB) fPdgDB = new TDatabasePDG();
133 //the codes are defined in $ROOTSYS/etc/pdg_table.txt
020abff0 134 fPionMass = fPdgDB->GetParticle("pi+")->Mass();
135 fPiPlusCode = fPdgDB->GetParticle("pi+")->PdgCode();
136 fPiMinusCode = fPdgDB->GetParticle("pi-")->PdgCode();
137 fKPlusCode = fPdgDB->GetParticle("K+")->PdgCode();
138 fKMinusCode = fPdgDB->GetParticle("K-")->PdgCode();
139 fProtonCode = fPdgDB->GetParticle("proton")->PdgCode();
140 fAntiProtonCode = fPdgDB->GetParticle("antiproton")->PdgCode();
141 fLambdaCode = fPdgDB->GetParticle("Lambda0")->PdgCode();
142 fAntiLambdaCode = fPdgDB->GetParticle("Lambda0_bar")->PdgCode();
143 fK0SCode = fPdgDB->GetParticle("K_S0")->PdgCode();
144 fOmegaCode = fPdgDB->GetParticle("Omega-")->PdgCode();
145 fAntiOmegaCode = fPdgDB->GetParticle("Omega+")->PdgCode();
146 fXi0Code = fPdgDB->GetParticle("Xi0")->PdgCode();
147 fAntiXi0Code = fPdgDB->GetParticle("Xi0_bar")->PdgCode();
148 fXiCode = fPdgDB->GetParticle("Xi-")->PdgCode();
149 fAntiXiCode = fPdgDB->GetParticle("Xi-_bar")->PdgCode();
150 fSigmaCode = fPdgDB->GetParticle("Sigma-")->PdgCode();
151 fAntiSigmaCode = fPdgDB->GetParticle("Sigma+")->PdgCode();
152 fK0LCode = fPdgDB->GetParticle("K_L0")->PdgCode();
153 fNeutronCode = fPdgDB->GetParticle("neutron")->PdgCode();
154 fAntiNeutronCode = fPdgDB->GetParticle("antineutron")->PdgCode();
155 fEPlusCode = fPdgDB->GetParticle("e+")->PdgCode();
156 fEMinusCode = fPdgDB->GetParticle("e-")->PdgCode();
641e1e0c 157}
158
159void AliAnalysisHadEt::CreateHistograms()
020abff0 160{//creates histograms included in all AliAnalysisHadEt's
641e1e0c 161}
162
163void AliAnalysisHadEt::FillHistograms()
020abff0 164{//Fills histograms filled for all AliAnalysisHadEt's
641e1e0c 165}
166
167void AliAnalysisHadEt::ResetEventValues()
020abff0 168{//Resets event values of et to zero
641e1e0c 169 fTotEt = 0;
170 fTotEtAcc = 0;
171 fTotNeutralEt = 0;
172 fTotNeutralEtAcc = 0;
173 fTotChargedEt = 0;
174 fTotChargedEtAcc = 0;
175 fMultiplicity = 0;
176 fChargedMultiplicity = 0;
177 fNeutralMultiplicity = 0;
178}
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}