- adding histograms for energy deposited by identified charged particles
[u/mrichter/AliRoot.git] / PWG4 / totEt / AliAnalysisEt.cxx
1 //_________________________________________________________________________
2 //  Utility Class for transverse energy studies
3 //  Base class for ESD & MC analysis
4 //  - reconstruction and MonteCarlo output
5 // implementation file
6 //
7 //*-- Authors: Oystein Djuvsland (Bergen), David Silvermyr (ORNL)
8 //_________________________________________________________________________
9
10 #include "AliAnalysisEt.h"
11 #include "TMath.h"
12 #include "TList.h"
13 #include "TH1F.h"
14 #include "TH2F.h"
15 #include <iostream>
16 #include "AliAnalysisEtCuts.h"
17 #include "AliVEvent.h"
18 #include "TDatabasePDG.h"
19 #include "Rtypes.h"
20
21 using namespace std;
22 ClassImp(AliAnalysisEt);
23
24
25 AliAnalysisEt::AliAnalysisEt() :
26         fHistogramNameSuffix("")
27         ,fCuts(0)
28         ,fPdgDB(0)
29         ,fPiPlusCode(0)
30         ,fPiMinusCode(0)
31         ,fKPlusCode(0)
32         ,fKMinusCode(0)
33         ,fProtonCode(0)
34         ,fAntiProtonCode(0)
35         ,fLambdaCode(0)
36         ,fAntiLambdaCode(0)
37         ,fK0SCode(0)
38         ,fOmegaCode(0)
39         ,fAntiOmegaCode(0)
40         ,fXi0Code(0)
41         ,fAntiXi0Code(0)
42         ,fXiCode(0)
43         ,fAntiXiCode(0)
44         ,fSigmaCode(0)
45         ,fAntiSigmaCode(0)
46         ,fK0LCode(0)
47         ,fNeutronCode(0)
48         ,fAntiNeutronCode(0)
49         ,fEPlusCode(0)
50         ,fEMinusCode(0)
51         ,fPionMass(0)
52         ,fSumEt(0)
53         ,fSumEtAcc(0)
54         ,fTotEt(0)
55         ,fTotEtAcc(0)
56         ,fTotNeutralEt(0)
57         ,fTotNeutralEtAcc(0)
58         ,fTotChargedEt(0)
59         ,fTotChargedEtAcc(0)
60         ,fMultiplicity(0)
61         ,fChargedMultiplicity(0)
62         ,fNeutralMultiplicity(0)
63         ,fBaryonEt(0)
64         ,fAntiBaryonEt(0)
65         ,fMesonEt(0)
66         ,fBaryonEtAcc(0)
67         ,fAntiBaryonEtAcc(0)
68         ,fMesonEtAcc(0)
69         ,fProtonEt(0)
70         ,fChargedKaonEt(0)
71         ,fMuonEt(0)
72         ,fElectronEt(0)
73         ,fProtonEtAcc(0)
74         ,fChargedKaonEtAcc(0)
75         ,fMuonEtAcc(0)
76         ,fElectronEtAcc(0)
77         ,fEtaCut(0)
78         ,fEtaCutAcc(0)
79         ,fPhiCutAccMin(0)
80         ,fPhiCutAccMax(0)
81         ,fDetectorRadius(0)
82         ,fClusterEnergyCut(0) 
83         ,fSingleCellEnergyCut(0)
84         ,fHistEt(0)
85         ,fHistChargedEt(0)
86         ,fHistNeutralEt(0)
87         ,fHistEtAcc(0)
88         ,fHistChargedEtAcc(0)
89         ,fHistNeutralEtAcc(0)
90         ,fHistMult(0)
91         ,fHistChargedMult(0)
92         ,fHistNeutralMult(0)
93         ,fHistPhivsPtPos(0)
94         ,fHistPhivsPtNeg(0)
95         ,fHistBaryonEt(0)
96         ,fHistAntiBaryonEt(0)
97         ,fHistMesonEt(0)
98         ,fHistBaryonEtAcc(0)
99         ,fHistAntiBaryonEtAcc(0)
100         ,fHistMesonEtAcc(0)
101         ,fHistProtonEt(0)
102         ,fHistChargedKaonEt(0)
103         ,fHistMuonEt(0)
104         ,fHistElectronEt(0)
105         ,fHistProtonEtAcc(0)
106         ,fHistChargedKaonEtAcc(0)
107         ,fHistMuonEtAcc(0)
108         ,fHistElectronEtAcc(0)
109         ,fHistEtRecvsEtMC(0)
110         ,fHistTMDeltaR(0)
111 {
112 }
113
114 AliAnalysisEt::~AliAnalysisEt()
115 {
116
117 }
118
119 void AliAnalysisEt::FillOutputList(TList *list)
120 { // histograms to be added to output
121     list->Add(fHistEt);
122     list->Add(fHistChargedEt);
123     list->Add(fHistNeutralEt);
124
125     list->Add(fHistEtAcc);
126     list->Add(fHistChargedEtAcc);
127     list->Add(fHistNeutralEtAcc);
128
129     list->Add(fHistMult);
130     list->Add(fHistChargedMult);
131     list->Add(fHistNeutralMult);
132
133     list->Add(fHistPhivsPtPos);
134     list->Add(fHistPhivsPtNeg);
135
136     list->Add(fHistBaryonEt);
137     list->Add(fHistAntiBaryonEt);
138     list->Add(fHistMesonEt);
139
140     list->Add(fHistBaryonEtAcc);
141     list->Add(fHistAntiBaryonEtAcc);
142     list->Add(fHistMesonEtAcc);
143
144     list->Add(fHistProtonEtAcc);
145     list->Add(fHistChargedKaonEtAcc);
146     list->Add(fHistMuonEtAcc);
147     list->Add(fHistElectronEtAcc);
148
149     list->Add(fHistEtRecvsEtMC);
150
151     list->Add(fHistTMDeltaR);
152 }
153
154 void AliAnalysisEt::Init()
155 {// clear variables, set up cuts and PDG info
156   ResetEventValues();
157 }
158
159 void AliAnalysisEt::CreateHistograms()
160 { // create histograms..
161   // histogram binning for E_T, p_T and Multiplicity: defaults for p+p
162   Int_t nbinsEt = 1000;
163   Double_t minEt = 0.0001;
164   Double_t maxEt = 100;
165   Int_t nbinsPt = 200;
166   Double_t minPt = 0;
167   Double_t maxPt = 20;
168   Int_t nbinsMult = 200;
169   Double_t minMult = -0.5; // offset -0.5 to have integer bins centered around 0
170   Double_t maxMult = nbinsMult + minMult; // 1 bin per integer value
171
172     TString histname = "fHistEt" + fHistogramNameSuffix;
173     fHistEt = new TH1F(histname.Data(), "Total E_{T} Distribution", nbinsEt, minEt, maxEt);
174     fHistEt->GetXaxis()->SetTitle("E_{T} (GeV/c^{2})");
175     fHistEt->GetYaxis()->SetTitle("dN/dE_{T} (c^{2}/GeV)");
176
177     histname = "fHistChargedEt" + fHistogramNameSuffix;
178     fHistChargedEt = new TH1F(histname.Data(), "Total Charged E_{T} Distribution", nbinsEt, minEt, maxEt);
179     fHistChargedEt->GetXaxis()->SetTitle("E_{T} (GeV/c^{2})");
180     fHistChargedEt->GetYaxis()->SetTitle("dN/dE_{T} (c^{2}/GeV)");
181
182     histname = "fHistNeutralEt" + fHistogramNameSuffix;
183     fHistNeutralEt = new TH1F(histname.Data(), "Total Neutral E_{T} Distribution", nbinsEt, minEt, maxEt);
184     fHistNeutralEt->GetXaxis()->SetTitle("E_{T} (GeV/c^{2})");
185     fHistNeutralEt->GetYaxis()->SetTitle("dN/dE_{T} (c^{2}/GeV)");
186
187     histname = "fHistEtAcc" + fHistogramNameSuffix;
188     fHistEtAcc = new TH1F(histname.Data(), "Total E_{T} Distribution in Acceptance", nbinsEt, minEt, maxEt);
189     fHistEtAcc->GetXaxis()->SetTitle("E_{T} (GeV/c^{2})");
190     fHistEtAcc->GetYaxis()->SetTitle("dN/dE_{T} (c^{2}/GeV)");
191
192     histname = "fHistChargedEtAcc" + fHistogramNameSuffix;
193     fHistChargedEtAcc = new TH1F(histname.Data(), "Total Charged E_{T} Distribution in Acceptance", nbinsEt, minEt, maxEt);
194     fHistChargedEtAcc->GetXaxis()->SetTitle("E_{T} (GeV/c^{2})");
195     fHistChargedEtAcc->GetYaxis()->SetTitle("dN/dE_{T} (c^{2}/GeV)");
196
197     histname = "fHistNeutralEtAcc" + fHistogramNameSuffix;
198     fHistNeutralEtAcc = new TH1F(histname.Data(), "Total Neutral E_{T} Distribution in Acceptance", nbinsEt, minEt, maxEt);
199     fHistNeutralEtAcc->GetXaxis()->SetTitle("E_{T} (GeV/c^{2})");
200     fHistNeutralEtAcc->GetYaxis()->SetTitle("dN/dE_{T} (c^{2}/GeV)");
201     std::cout << histname << std::endl;
202     histname = "fHistMult" + fHistogramNameSuffix;
203     fHistMult = new TH1F(histname.Data(), "Total Multiplicity", nbinsMult, minMult, maxMult);
204     fHistMult->GetXaxis()->SetTitle("N");
205     fHistMult->GetYaxis()->SetTitle("Multiplicity");
206
207     histname = "fHistChargedMult" + fHistogramNameSuffix;
208     fHistChargedMult = new TH1F(histname.Data(), "Charged Multiplicity", nbinsMult, minMult, maxMult);
209     fHistChargedMult->GetXaxis()->SetTitle("N");
210     fHistChargedMult->GetYaxis()->SetTitle("Multiplicity");
211
212     histname = "fHistNeutralMult" + fHistogramNameSuffix;
213     fHistNeutralMult = new TH1F(histname.Data(), "Neutral Multiplicity", nbinsMult, minMult, maxMult);
214     fHistNeutralMult->GetXaxis()->SetTitle("N");
215     fHistNeutralMult->GetYaxis()->SetTitle("Multiplicity");
216
217     histname = "fHistPhivsPtPos" + fHistogramNameSuffix;
218     fHistPhivsPtPos = new TH2F(histname.Data(), "Phi vs pT of positively charged tracks hitting the calorimeter",       200, 0, 2*TMath::Pi(), nbinsPt, minPt, maxPt);
219
220     histname = "fHistPhivsPtNeg" + fHistogramNameSuffix;
221     fHistPhivsPtNeg = new TH2F(histname.Data(), "Phi vs pT of negatively charged tracks hitting the calorimeter",       200, 0, 2*TMath::Pi(), nbinsPt, minPt, maxPt);
222
223     histname = "fHistBaryonEt" + fHistogramNameSuffix;
224     fHistBaryonEt = new TH1F(histname.Data(), "E_{T} for baryons",  nbinsEt, minEt, maxEt);
225
226     histname = "fHistAntiBaryonEt" + fHistogramNameSuffix;
227     fHistAntiBaryonEt = new TH1F(histname.Data(), "E_{T} for anti baryons",  nbinsEt, minEt, maxEt);
228
229     histname = "fHistMesonEt" + fHistogramNameSuffix;
230     fHistMesonEt = new TH1F(histname.Data(), "E_{T} for mesons",  nbinsEt, minEt, maxEt);
231
232     histname = "fHistBaryonEtAcc" + fHistogramNameSuffix;
233     fHistBaryonEtAcc = new TH1F(histname.Data(), "E_{T} for baryons in calorimeter acceptance",  nbinsEt, minEt, maxEt);
234
235     histname = "fHistAntiBaryonEtAcc" + fHistogramNameSuffix;
236     fHistAntiBaryonEtAcc = new TH1F(histname.Data(), "E_{T} for anti baryons in calorimeter acceptance",  nbinsEt, minEt, maxEt);
237
238     histname = "fHistMesonEtAcc" + fHistogramNameSuffix;
239     fHistMesonEtAcc = new TH1F(histname.Data(), "E_{T} for mesons in calorimeter acceptance",  nbinsEt, minEt, maxEt);
240
241     histname = "fHistProtonEt" + fHistogramNameSuffix;
242     fHistProtonEt = new TH1F(histname.Data(), "E_{T} for (anti-)protons", nbinsEt, minEt, maxEt);
243
244     histname = "fHistKaonEt" + fHistogramNameSuffix;
245     fHistChargedKaonEt = new TH1F(histname.Data(), "E_{T} for charged kaons", nbinsEt, minEt, maxEt);
246
247     histname = "fHistMuonEt" + fHistogramNameSuffix;
248     fHistMuonEt = new TH1F(histname.Data(), "E_{T} for muons", nbinsEt, minEt, maxEt);
249
250     histname = "fHistElectronEt" + fHistogramNameSuffix;
251     fHistElectronEt = new TH1F(histname.Data(), "E_{T} for electrons/positrons", nbinsEt, minEt, maxEt);
252
253     histname = "fHistProtonEtAcc" + fHistogramNameSuffix;
254     fHistProtonEtAcc = new TH1F(histname.Data(), "E_{T} for (anti-)protons in calorimeter acceptance", nbinsEt, minEt, maxEt);
255
256     histname = "fHistKaonEtAcc" + fHistogramNameSuffix;
257     fHistChargedKaonEtAcc = new TH1F(histname.Data(), "E_{T} for charged kaons in calorimeter acceptance", nbinsEt, minEt, maxEt);
258
259     histname = "fHistMuonEtAcc" + fHistogramNameSuffix;
260     fHistMuonEtAcc = new TH1F(histname.Data(), "E_{T} for muons in calorimeter acceptance", nbinsEt, minEt, maxEt);
261
262     histname = "fHistElectronEtAcc" + fHistogramNameSuffix;
263     fHistElectronEtAcc = new TH1F(histname.Data(), "E_{T} for electrons/positrons in calorimeter acceptance", nbinsEt, minEt, maxEt);
264
265     histname = "fHistEtRecvsEtMC" + fHistogramNameSuffix;
266     fHistEtRecvsEtMC = new TH2F(histname.Data(), "Reconstructed E_{t} vs MC E_{t}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
267
268     //
269     histname = "fHistTMDeltaR" + fHistogramNameSuffix;
270     fHistTMDeltaR = new TH1F(histname.Data(), "#Delta R for calorimeter clusters", 200, 0, 50);
271
272 }
273
274 void AliAnalysisEt::FillHistograms()
275 { // fill histograms..
276     fHistEt->Fill(fTotEt);
277     fHistChargedEt->Fill(fTotChargedEt);
278     fHistNeutralEt->Fill(fTotNeutralEt);
279
280     fHistEtAcc->Fill(fTotEtAcc);
281     fHistChargedEtAcc->Fill(fTotChargedEtAcc);
282     fHistNeutralEtAcc->Fill(fTotNeutralEtAcc);
283
284     fHistMult->Fill(fMultiplicity);
285     fHistChargedMult->Fill(fChargedMultiplicity);
286     fHistNeutralMult->Fill(fNeutralMultiplicity);
287
288     /* // DS commented out non-fills to prevent compilation warnings
289     fHistPhivsPtPos;
290     fHistPhivsPtNeg;
291
292     fHistBaryonEt;
293     fHistAntiBaryonEt;
294     fHistMesonEt;
295
296     fHistBaryonEtAcc;
297     fHistAntiBaryonEtAcc;
298     fHistMesonEtAcc;
299
300     fHistTMDeltaR;
301     */
302     fHistProtonEt->Fill(fProtonEt);
303     fHistChargedKaonEt->Fill(fChargedKaonEt);
304     fHistMuonEt->Fill(fMuonEt);
305     fHistElectronEt->Fill(fElectronEt);
306     
307     fHistProtonEtAcc->Fill(fProtonEtAcc);
308     fHistChargedKaonEtAcc->Fill(fChargedKaonEtAcc);
309     fHistMuonEtAcc->Fill(fMuonEtAcc);
310     fHistElectronEtAcc->Fill(fElectronEtAcc);
311 }
312
313 Int_t AliAnalysisEt::AnalyseEvent(AliVEvent *event)
314 { //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.
315   cout << "This event has " << event->GetNumberOfTracks() << " tracks" << endl;
316   ResetEventValues();
317   return 0;
318 }
319
320 void AliAnalysisEt::ResetEventValues()
321 { // clear
322   fTotEt = 0;
323   fTotEtAcc = 0;
324   fTotNeutralEt = 0;
325   fTotNeutralEtAcc = 0;
326   fTotChargedEt  = 0;
327   fTotChargedEtAcc = 0;
328   fMultiplicity = 0;
329   fChargedMultiplicity = 0;
330   fNeutralMultiplicity = 0;
331   fBaryonEt = 0;
332   fAntiBaryonEt = 0;
333   fMesonEt = 0;
334   fBaryonEtAcc = 0;
335   fAntiBaryonEtAcc = 0;
336   fMesonEtAcc = 0;
337   fProtonEt = 0;
338   fChargedKaonEt = 0;
339   fMuonEt = 0;
340   fElectronEt = 0;
341   fProtonEtAcc = 0;
342   fChargedKaonEtAcc = 0;
343   fMuonEtAcc = 0;
344     
345   if (!fCuts || !fPdgDB || fPiPlusCode==0) { // some Init's needed
346     cout << __FILE__ << ":" << __LINE__ << " : Init " << endl;
347     if (!fCuts) {
348       cout << " setting up Cuts " << endl;
349       fCuts = new AliAnalysisEtCuts();
350     }
351     if(!fPdgDB) {
352       cout << " setting up PdgDB " << endl;
353       fPdgDB = new TDatabasePDG();
354     }
355     
356     if (fPiPlusCode==0) {
357       SetParticleCodes();
358     }
359   }
360   return;
361 }
362
363
364 void AliAnalysisEt::SetParticleCodes()
365 { // set PDG info    
366   fPionMass = fPdgDB->GetParticle("pi+")->Mass();
367   fPiPlusCode = fPdgDB->GetParticle("pi+")->PdgCode();
368   fPiMinusCode = fPdgDB->GetParticle("pi-")->PdgCode();
369   fKPlusCode = fPdgDB->GetParticle("K+")->PdgCode();
370   fKMinusCode = fPdgDB->GetParticle("K-")->PdgCode();
371   fProtonCode = fPdgDB->GetParticle("proton")->PdgCode();
372   fAntiProtonCode = fPdgDB->GetParticle("antiproton")->PdgCode();
373   fLambdaCode = fPdgDB->GetParticle("Lambda0")->PdgCode();
374   fAntiLambdaCode = fPdgDB->GetParticle("Lambda0_bar")->PdgCode();
375   fK0SCode = fPdgDB->GetParticle("K_S0")->PdgCode();
376   fOmegaCode = fPdgDB->GetParticle("Omega-")->PdgCode();
377   fAntiOmegaCode = fPdgDB->GetParticle("Omega+")->PdgCode();
378   fXi0Code = fPdgDB->GetParticle("Xi0")->PdgCode();
379   fAntiXi0Code = fPdgDB->GetParticle("Xi0_bar")->PdgCode();
380   fXiCode = fPdgDB->GetParticle("Xi-")->PdgCode();
381   fAntiXiCode = fPdgDB->GetParticle("Xi-_bar")->PdgCode();
382   fSigmaCode = fPdgDB->GetParticle("Sigma-")->PdgCode();
383   fAntiSigmaCode = fPdgDB->GetParticle("Sigma+")->PdgCode();
384   fK0LCode = fPdgDB->GetParticle("K_L0")->PdgCode();
385   fNeutronCode = fPdgDB->GetParticle("neutron")->PdgCode();
386   fAntiNeutronCode = fPdgDB->GetParticle("antineutron")->PdgCode();
387   fEPlusCode = fPdgDB->GetParticle("e+")->PdgCode();
388   fEMinusCode = fPdgDB->GetParticle("e-")->PdgCode();
389   cout << "Resetting Codes: Pion " << fPiPlusCode
390        << "," << fPiMinusCode 
391        << " Kaon " << fKPlusCode 
392        << "," << fKMinusCode << endl;
393 }
394