d7060435904b7f3f02b7076985705a05c699bdcd
[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         ,fEtaCutAcc(0)
64         ,fPhiCutAccMin(0)
65         ,fPhiCutAccMax(0)
66         ,fDetectorRadius(0)
67         ,fClusterEnergyCut(0) 
68         ,fSingleCellEnergyCut(0)
69         ,fHistEt(0)
70         ,fHistChargedEt(0)
71         ,fHistNeutralEt(0)
72         ,fHistEtAcc(0)
73         ,fHistChargedEtAcc(0)
74         ,fHistNeutralEtAcc(0)
75         ,fHistMult(0)
76         ,fHistChargedMult(0)
77         ,fHistNeutralMult(0)
78         ,fHistPhivsPtPos(0)
79         ,fHistPhivsPtNeg(0)
80         ,fHistBaryonEt(0)
81         ,fHistAntiBaryonEt(0)
82         ,fHistMesonEt(0)
83         ,fHistBaryonEtAcc(0)
84         ,fHistAntiBaryonEtAcc(0)
85         ,fHistMesonEtAcc(0)
86         ,fHistTMDeltaR(0)
87 {
88 }
89
90 AliAnalysisEt::~AliAnalysisEt()
91 {
92
93 }
94
95 void AliAnalysisEt::FillOutputList(TList *list)
96 { // histograms to be added to output
97     list->Add(fHistEt);
98     list->Add(fHistChargedEt);
99     list->Add(fHistNeutralEt);
100
101     list->Add(fHistEtAcc);
102     list->Add(fHistChargedEtAcc);
103     list->Add(fHistNeutralEtAcc);
104
105     list->Add(fHistMult);
106     list->Add(fHistChargedMult);
107     list->Add(fHistNeutralMult);
108
109     list->Add(fHistPhivsPtPos);
110     list->Add(fHistPhivsPtNeg);
111
112     list->Add(fHistBaryonEt);
113     list->Add(fHistAntiBaryonEt);
114     list->Add(fHistMesonEt);
115
116     list->Add(fHistBaryonEtAcc);
117     list->Add(fHistAntiBaryonEtAcc);
118     list->Add(fHistMesonEtAcc);
119
120     list->Add(fHistTMDeltaR);
121 }
122
123 void AliAnalysisEt::Init()
124 {// clear variables, set up cuts and PDG info
125   ResetEventValues();
126 }
127
128 void AliAnalysisEt::CreateHistograms()
129 { // create histograms..
130   // histogram binning for E_T, p_T and Multiplicity: defaults for p+p
131   Int_t nbinsEt = 1000;
132   Double_t minEt = 0.0001;
133   Double_t maxEt = 100;
134   Int_t nbinsPt = 200;
135   Double_t minPt = 0;
136   Double_t maxPt = 20;
137   Int_t nbinsMult = 200;
138   Double_t minMult = -0.5; // offset -0.5 to have integer bins centered around 0
139   Double_t maxMult = nbinsMult + minMult; // 1 bin per integer value
140
141     TString histname = "fHistEt" + fHistogramNameSuffix;
142     fHistEt = new TH1F(histname.Data(), "Total E_{T} Distribution", nbinsEt, minEt, maxEt);
143     fHistEt->GetXaxis()->SetTitle("E_{T} (GeV/c^{2})");
144     fHistEt->GetYaxis()->SetTitle("dN/dE_{T} (c^{2}/GeV)");
145
146     histname = "fHistChargedEt" + fHistogramNameSuffix;
147     fHistChargedEt = new TH1F(histname.Data(), "Total Charged E_{T} Distribution", nbinsEt, minEt, maxEt);
148     fHistChargedEt->GetXaxis()->SetTitle("E_{T} (GeV/c^{2})");
149     fHistChargedEt->GetYaxis()->SetTitle("dN/dE_{T} (c^{2}/GeV)");
150
151     histname = "fHistNeutralEt" + fHistogramNameSuffix;
152     fHistNeutralEt = new TH1F(histname.Data(), "Total Neutral E_{T} Distribution", nbinsEt, minEt, maxEt);
153     fHistNeutralEt->GetXaxis()->SetTitle("E_{T} (GeV/c^{2})");
154     fHistNeutralEt->GetYaxis()->SetTitle("dN/dE_{T} (c^{2}/GeV)");
155
156     histname = "fHistEtAcc" + fHistogramNameSuffix;
157     fHistEtAcc = new TH1F(histname.Data(), "Total E_{T} Distribution in Acceptance", nbinsEt, minEt, maxEt);
158     fHistEtAcc->GetXaxis()->SetTitle("E_{T} (GeV/c^{2})");
159     fHistEtAcc->GetYaxis()->SetTitle("dN/dE_{T} (c^{2}/GeV)");
160
161     histname = "fHistChargedEtAcc" + fHistogramNameSuffix;
162     fHistChargedEtAcc = new TH1F(histname.Data(), "Total Charged E_{T} Distribution in Acceptance", nbinsEt, minEt, maxEt);
163     fHistChargedEtAcc->GetXaxis()->SetTitle("E_{T} (GeV/c^{2})");
164     fHistChargedEtAcc->GetYaxis()->SetTitle("dN/dE_{T} (c^{2}/GeV)");
165
166     histname = "fHistNeutralEtAcc" + fHistogramNameSuffix;
167     fHistNeutralEtAcc = new TH1F(histname.Data(), "Total Neutral E_{T} Distribution in Acceptance", nbinsEt, minEt, maxEt);
168     fHistNeutralEtAcc->GetXaxis()->SetTitle("E_{T} (GeV/c^{2})");
169     fHistNeutralEtAcc->GetYaxis()->SetTitle("dN/dE_{T} (c^{2}/GeV)");
170     std::cout << histname << std::endl;
171     histname = "fHistMult" + fHistogramNameSuffix;
172     fHistMult = new TH1F(histname.Data(), "Total Multiplicity", nbinsMult, minMult, maxMult);
173     fHistMult->GetXaxis()->SetTitle("N");
174     fHistMult->GetYaxis()->SetTitle("Multiplicity");
175
176     histname = "fHistChargedMult" + fHistogramNameSuffix;
177     fHistChargedMult = new TH1F(histname.Data(), "Charged Multiplicity", nbinsMult, minMult, maxMult);
178     fHistChargedMult->GetXaxis()->SetTitle("N");
179     fHistChargedMult->GetYaxis()->SetTitle("Multiplicity");
180
181     histname = "fHistNeutralMult" + fHistogramNameSuffix;
182     fHistNeutralMult = new TH1F(histname.Data(), "Neutral Multiplicity", nbinsMult, minMult, maxMult);
183     fHistNeutralMult->GetXaxis()->SetTitle("N");
184     fHistNeutralMult->GetYaxis()->SetTitle("Multiplicity");
185
186     histname = "fHistPhivsPtPos" + fHistogramNameSuffix;
187     fHistPhivsPtPos = new TH2F(histname.Data(), "Phi vs pT of positively charged tracks hitting the calorimeter",       200, 0, 2*TMath::Pi(), nbinsPt, minPt, maxPt);
188
189     histname = "fHistPhivsPtNeg" + fHistogramNameSuffix;
190     fHistPhivsPtNeg = new TH2F(histname.Data(), "Phi vs pT of negatively charged tracks hitting the calorimeter",       200, 0, 2*TMath::Pi(), nbinsPt, minPt, maxPt);
191
192     histname = "fHistBaryonEt" + fHistogramNameSuffix;
193     fHistBaryonEt = new TH1F(histname.Data(), "E_{T} for baryons",  nbinsEt, minEt, maxEt);
194
195     histname = "fHistAntiBaryonEt" + fHistogramNameSuffix;
196     fHistAntiBaryonEt = new TH1F(histname.Data(), "E_{T} for anti baryons",  nbinsEt, minEt, maxEt);
197
198     histname = "fHistMesonEt" + fHistogramNameSuffix;
199     fHistMesonEt = new TH1F(histname.Data(), "E_{T} for mesons",  nbinsEt, minEt, maxEt);
200
201     histname = "fHistBaryonEtAcc" + fHistogramNameSuffix;
202     fHistBaryonEtAcc = new TH1F(histname.Data(), "E_{T} for baryons in calorimeter acceptance",  nbinsEt, minEt, maxEt);
203
204     histname = "fHistAntiBaryonEtAcc" + fHistogramNameSuffix;
205     fHistAntiBaryonEtAcc = new TH1F(histname.Data(), "E_{T} for anti baryons in calorimeter acceptance",  nbinsEt, minEt, maxEt);
206
207     histname = "fHistMesonEtAcc" + fHistogramNameSuffix;
208     fHistMesonEtAcc = new TH1F(histname.Data(), "E_{T} for mesons in calorimeter acceptance",  nbinsEt, minEt, maxEt);
209
210     //
211     histname = "fHistTMDeltaR" + fHistogramNameSuffix;
212     fHistTMDeltaR = new TH1F(histname.Data(), "#Delta R for calorimeter clusters", 200, 0, 50);
213
214 }
215
216 void AliAnalysisEt::FillHistograms()
217 { // fill histograms..
218     fHistEt->Fill(fTotEt);
219     fHistChargedEt->Fill(fTotChargedEt);
220     fHistNeutralEt->Fill(fTotNeutralEt);
221
222     fHistEtAcc->Fill(fTotEtAcc);
223     fHistChargedEtAcc->Fill(fTotChargedEtAcc);
224     fHistNeutralEtAcc->Fill(fTotNeutralEtAcc);
225
226     fHistMult->Fill(fMultiplicity);
227     fHistChargedMult->Fill(fChargedMultiplicity);
228     fHistNeutralMult->Fill(fNeutralMultiplicity);
229
230     /* // DS commented out non-fills to prevent compilation warnings
231     fHistPhivsPtPos;
232     fHistPhivsPtNeg;
233
234     fHistBaryonEt;
235     fHistAntiBaryonEt;
236     fHistMesonEt;
237
238     fHistBaryonEtAcc;
239     fHistAntiBaryonEtAcc;
240     fHistMesonEtAcc;
241
242     fHistTMDeltaR;
243     */
244 }
245
246 Int_t AliAnalysisEt::AnalyseEvent(AliVEvent *event)
247 { //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.
248   cout << "This event has " << event->GetNumberOfTracks() << " tracks" << endl;
249   ResetEventValues();
250   return 0;
251 }
252
253 void AliAnalysisEt::ResetEventValues()
254 { // clear
255   fTotEt = 0;
256   fTotEtAcc = 0;
257   fTotNeutralEt = 0;
258   fTotNeutralEtAcc = 0;
259   fTotChargedEt  = 0;
260   fTotChargedEtAcc = 0;
261   fMultiplicity = 0;
262   fChargedMultiplicity = 0;
263   fNeutralMultiplicity = 0;
264   
265   if (!fCuts || !fPdgDB || fPiPlusCode==0) { // some Init's needed
266     cout << __FILE__ << ":" << __LINE__ << " : Init " << endl;
267     if (!fCuts) {
268       cout << " setting up Cuts " << endl;
269       fCuts = new AliAnalysisEtCuts();
270     }
271     if(!fPdgDB) {
272       cout << " setting up PdgDB " << endl;
273       fPdgDB = new TDatabasePDG();
274     }
275     
276     if (fPiPlusCode==0) {
277       SetParticleCodes();
278     }
279   }
280   return;
281 }
282
283
284 void AliAnalysisEt::SetParticleCodes()
285 { // set PDG info    
286   fPionMass = fPdgDB->GetParticle("pi+")->Mass();
287   fPiPlusCode = fPdgDB->GetParticle("pi+")->PdgCode();
288   fPiMinusCode = fPdgDB->GetParticle("pi-")->PdgCode();
289   fKPlusCode = fPdgDB->GetParticle("K+")->PdgCode();
290   fKMinusCode = fPdgDB->GetParticle("K-")->PdgCode();
291   fProtonCode = fPdgDB->GetParticle("proton")->PdgCode();
292   fAntiProtonCode = fPdgDB->GetParticle("antiproton")->PdgCode();
293   fLambdaCode = fPdgDB->GetParticle("Lambda0")->PdgCode();
294   fAntiLambdaCode = fPdgDB->GetParticle("Lambda0_bar")->PdgCode();
295   fK0SCode = fPdgDB->GetParticle("K_S0")->PdgCode();
296   fOmegaCode = fPdgDB->GetParticle("Omega-")->PdgCode();
297   fAntiOmegaCode = fPdgDB->GetParticle("Omega+")->PdgCode();
298   fXi0Code = fPdgDB->GetParticle("Xi0")->PdgCode();
299   fAntiXi0Code = fPdgDB->GetParticle("Xi0_bar")->PdgCode();
300   fXiCode = fPdgDB->GetParticle("Xi-")->PdgCode();
301   fAntiXiCode = fPdgDB->GetParticle("Xi-_bar")->PdgCode();
302   fSigmaCode = fPdgDB->GetParticle("Sigma-")->PdgCode();
303   fAntiSigmaCode = fPdgDB->GetParticle("Sigma+")->PdgCode();
304   fK0LCode = fPdgDB->GetParticle("K_L0")->PdgCode();
305   fNeutronCode = fPdgDB->GetParticle("neutron")->PdgCode();
306   fAntiNeutronCode = fPdgDB->GetParticle("antineutron")->PdgCode();
307   fEPlusCode = fPdgDB->GetParticle("e+")->PdgCode();
308   fEMinusCode = fPdgDB->GetParticle("e-")->PdgCode();
309   cout << "Resetting Codes: Pion " << fPiPlusCode
310        << "," << fPiMinusCode 
311        << " Kaon " << fKPlusCode 
312        << "," << fKMinusCode << endl;
313 }
314