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