- adding histograms for energy deposited by identified charged particles
[u/mrichter/AliRoot.git] / PWG4 / totEt / AliAnalysisEtReconstructed.cxx
CommitLineData
cf6522d1 1//_________________________________________________________________________
2// Utility Class for transverse energy studies
3// Base class for ESD analysis
4// - reconstruction output
5// implementation file
6//
7//*-- Authors: Oystein Djuvsland (Bergen), David Silvermyr (ORNL)
8//_________________________________________________________________________
2fbf38ac 9
10#include "AliAnalysisEtReconstructed.h"
11#include "AliAnalysisEtCuts.h"
12#include "AliESDtrack.h"
13#include "AliESDCaloCluster.h"
14#include "TVector3.h"
15#include "AliVEvent.h"
16#include "AliESDEvent.h"
17#include "AliVParticle.h"
87efb15c 18#include "TDatabasePDG.h"
19#include "TList.h"
2fbf38ac 20#include <iostream>
21#include "TH2F.h"
22
23AliAnalysisEtReconstructed::AliAnalysisEtReconstructed() :
24 AliAnalysisEt()
87efb15c 25 ,fNTpcClustersCut(0)
26 ,fNItsClustersCut(0)
2fbf38ac 27 ,fTrackDistanceCut(0)
87efb15c 28 ,fPidCut(0)
2fbf38ac 29 ,fClusterType(0)
87efb15c 30 ,fHistChargedPionEnergyDeposit(0)
31 ,fHistProtonEnergyDeposit(0)
32 ,fHistAntiProtonEnergyDeposit(0)
33 ,fHistChargedKaonEnergyDeposit(0)
34 ,fHistMuonEnergyDeposit(0)
2fbf38ac 35{
36
37}
38
cf6522d1 39AliAnalysisEtReconstructed::~AliAnalysisEtReconstructed()
2fbf38ac 40{
cf6522d1 41}
42
43Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
44{ // analyse ESD event
2fbf38ac 45 ResetEventValues();
46 AliESDEvent *event = dynamic_cast<AliESDEvent*>(ev);
47
87efb15c 48 Double_t protonMass = fPdgDB->GetParticle("proton")->Mass();
49
2fbf38ac 50 for (Int_t iTrack = 0; iTrack < event->GetNumberOfTracks(); iTrack++)
51 {
52 AliVParticle *track = event->GetTrack(iTrack);
53 if (!track)
54 {
55 Printf("ERROR: Could not get track %d", iTrack);
56 continue;
57 }
58
59 fMultiplicity++;
60
61 Int_t nItsClusters = dynamic_cast<AliESDtrack*>(track)->GetNcls(0);
62 Int_t nTPCClusters = dynamic_cast<AliESDtrack*>(track)->GetNcls(1);
63
64 Float_t massPart = 0;
65
66 const Double_t *pidWeights = track->PID();
87efb15c 67 Int_t maxpid = -1;
68 Double_t maxpidweight = 0;
69
2fbf38ac 70 if (pidWeights)
71 {
2fbf38ac 72 for (Int_t p =0; p < AliPID::kSPECIES; p++)
73 {
74 if (pidWeights[p] > maxpidweight)
75 {
76 maxpidweight = pidWeights[p];
77 maxpid = p;
78 }
79 }
80 if (maxpid == AliPID::kProton)
81 {
87efb15c 82 //by definition of ET
83 massPart = -protonMass*track->Charge();
2fbf38ac 84 }
85
86 }
87
88 Double_t et = track->E() * TMath::Sin(track->Theta()) + massPart;
cf6522d1 89 // printf("Rec track: iTrack %03d eta %4.3f phi %4.3f nITSCl %d nTPCCl %d\n", iTrack, track->Eta(), track->Phi(), nItsClusters, nTPCClusters); // tmp/debug printout
2fbf38ac 90
4998becf 91 if (TMath::Abs(track->Eta()) < fCuts->GetCommonEtaCut() && CheckGoodVertex(track) && nItsClusters > fCuts->GetReconstructedNItsClustersCut() && nTPCClusters > fCuts->GetReconstructedNTpcClustersCut() )
2fbf38ac 92 {
93 fTotChargedEt += et;
94 fChargedMultiplicity++;
87efb15c 95 if (maxpid != -1)
96 {
97 if (maxpid == AliPID::kPion)
98 {
99 fProtonEt += et;
100 }
101 if (maxpid == AliPID::kKaon)
102 {
103 fChargedKaonEt += et;
104 }
105 if (maxpid == AliPID::kMuon)
106 {
107 fMuonEt += et;
108 }
109 if (maxpid == AliPID::kElectron)
110 {
111 fElectronEt += et;
112 }
113 }
2fbf38ac 114
115 if (TMath::Abs(track->Eta()) < fEtaCutAcc && track->Phi() < fPhiCutAccMax && track->Phi() > fPhiCutAccMin)
116 {
117 fTotChargedEtAcc += track->E()*TMath::Sin(track->Theta()) + massPart;
87efb15c 118 if (maxpid != -1)
119 {
120 if (maxpid == AliPID::kPion)
121 {
122 fProtonEtAcc += et;
123 }
124 if (maxpid == AliPID::kKaon)
125 {
126 fChargedKaonEtAcc += et;
127 }
128 if (maxpid == AliPID::kMuon)
129 {
130 fMuonEtAcc += et;
131 }
132 if (maxpid == AliPID::kElectron)
133 {
134 fElectronEtAcc += et;
135 }
136 }
137
2fbf38ac 138 }
139 }
140
2fbf38ac 141 if (TrackHitsCalorimeter(track, event->GetMagneticField()))
142 {
13b0d3c1 143 Double_t phi = track->Phi();
144 Double_t pt = track->Pt();
cf6522d1 145 // printf("Rec track hit: iTrack %03d phi %4.3f pt %4.3f\n", iTrack, phi, pt); // tmp/debug printout
13b0d3c1 146 if (track->Charge() > 0) fHistPhivsPtPos->Fill(phi, pt);
147 else fHistPhivsPtNeg->Fill(phi, pt);
2fbf38ac 148 }
149 }
150
151 for (Int_t iCluster = 0; iCluster < event->GetNumberOfCaloClusters(); iCluster++)
152 {
153 AliESDCaloCluster* cluster = event->GetCaloCluster(iCluster);
154 if (!cluster)
155 {
156 Printf("ERROR: Could not get cluster %d", iCluster);
157 continue;
158 }
159
87efb15c 160 if (cluster->GetType() != fClusterType) continue;
161 if(cluster->GetTracksMatched() > 0)
162 //printf("Rec Cluster: iCluster %03d E %4.3f type %d NCells %d, nmatched: %d, distance to closest: %f\n", iCluster, cluster->E(), (int)(cluster->GetType()), cluster->GetNCells(), cluster->GetNTracksMatched(), cluster->GetEmcCpvDistance()); // tmp/debug printout
163
164
2fbf38ac 165 if (cluster->E() < fClusterEnergyCut) continue;
166 Float_t pos[3];
167 TVector3 cp(pos);
168 cluster->GetPosition(pos);
87efb15c 169
170 fHistTMDeltaR->Fill(cluster->GetEmcCpvDistance());
2fbf38ac 171 if (cluster->GetEmcCpvDistance() < fTrackDistanceCut)
172 {
87efb15c 173 if (cluster->GetNTracksMatched() == 1)
174 {
175 AliVTrack *track = event->GetTrack(cluster->GetTrackMatchedIndex());
176 const Double_t *pidWeights = track->PID();
177
178 Double_t maxpidweight = 0;
179 Int_t maxpid = 0;
180
181 if (pidWeights)
182 {
183 for (Int_t p =0; p < AliPID::kSPECIES; p++)
184 {
185 if (pidWeights[p] > maxpidweight)
186 {
187 maxpidweight = pidWeights[p];
188 maxpid = p;
189 }
190 }
191 if(maxpidweight > fPidCut)
192 {
193 if(maxpid == AliPID::kProton)
194 {
195 if(track->Charge() == 1)
196 {
197 fHistProtonEnergyDeposit->Fill(cluster->E(), track->E());
198 }
199 else if(track->Charge() == -1)
200 {
201 fHistAntiProtonEnergyDeposit->Fill(cluster->E(), track->E());
202 }
203 }
204 else if(maxpid == AliPID::kPion)
205 {
206 fHistChargedPionEnergyDeposit->Fill(cluster->E(), track->E());
207 }
208 else if(maxpid == AliPID::kKaon)
209 {
210 fHistChargedKaonEnergyDeposit->Fill(cluster->E(), track->E());
211 }
212 else if(maxpid == AliPID::kMuon)
213 {
214 fHistMuonEnergyDeposit->Fill(cluster->E(), track->E());
215 }
216 }
217 }
218 }
219
2fbf38ac 220 continue;
2fbf38ac 221 }
222
4998becf 223 if (cluster->E() > fSingleCellEnergyCut && cluster->GetNCells() == fCuts->GetCommonSingleCell()) continue;
2fbf38ac 224
225 cluster->GetPosition(pos);
226
227 // TODO: replace with TVector3, too lazy now...
228
229 float dist = TMath::Sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
230
231 float theta = TMath::ATan(pos[2]/dist)+TMath::Pi()/2;
232 // float eta = TMath::Log(TMath::Abs( TMath::Tan( 0.5 * theta ) ) );
233 fTotNeutralEt += cluster->E() * TMath::Sin(theta);
13b0d3c1 234 fNeutralMultiplicity++;
2fbf38ac 235
236 fMultiplicity++;
2fbf38ac 237 }
238
239 fTotNeutralEtAcc = fTotNeutralEt;
240 fTotEt = fTotChargedEt + fTotNeutralEt;
241 fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;
242
2fbf38ac 243 // Fill the histograms...
244 FillHistograms();
245
246 return 0;
247}
248
249bool AliAnalysisEtReconstructed::CheckGoodVertex(AliVParticle* track)
cf6522d1 250{ // check vertex
2fbf38ac 251
252 Float_t bxy = 999.;
253 Float_t bz = 999.;
254 dynamic_cast<AliESDtrack*>(track)->GetImpactParametersTPC(bxy,bz);
255
4998becf 256 bool status = (TMath::Abs(track->Xv()) < fCuts->GetReconstructedVertexXCut()) &&
257 (TMath::Abs(track->Yv()) < fCuts->GetReconstructedVertexYCut()) &&
258 (TMath::Abs(track->Zv()) < fCuts->GetReconstructedVertexZCut()) &&
259 (TMath::Abs(bxy) < fCuts->GetReconstructedIPxyCut()) &&
260 (TMath::Abs(bz) < fCuts->GetReconstructedIPzCut());
2fbf38ac 261
4998becf 262 return status;
2fbf38ac 263}
264
265void AliAnalysisEtReconstructed::Init()
cf6522d1 266{ // Init
2fbf38ac 267 AliAnalysisEt::Init();
87efb15c 268 fNItsClustersCut = fCuts->GetReconstructedNItsClustersCut();
269 fNTpcClustersCut = fCuts->GetReconstructedNTpcClustersCut();
270 fPidCut = fCuts->GetCommonPidCut();
2fbf38ac 271}
272
273bool AliAnalysisEtReconstructed::TrackHitsCalorimeter(AliVParticle* track, Double_t magField)
cf6522d1 274{ // propagate track to detector radius
2fbf38ac 275
276 AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
cf6522d1 277 // Printf("Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
2fbf38ac 278
279 Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);
280
cf6522d1 281 // if (prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
2fbf38ac 282 return prop &&
283 TMath::Abs(esdTrack->Eta()) < fEtaCutAcc &&
284 esdTrack->Phi() > fPhiCutAccMin*TMath::Pi()/180. &&
285 esdTrack->Phi() < fPhiCutAccMax*TMath::Pi()/180.;
286}
287
87efb15c 288void AliAnalysisEtReconstructed::FillOutputList(TList* list)
289{
290 AliAnalysisEt::FillOutputList(list);
291
292 list->Add(fHistChargedPionEnergyDeposit);
293 list->Add(fHistProtonEnergyDeposit);
294 list->Add(fHistAntiProtonEnergyDeposit);
295 list->Add(fHistChargedKaonEnergyDeposit);
296 list->Add(fHistMuonEnergyDeposit);
297}
298
299void AliAnalysisEtReconstructed::CreateHistograms()
300{
301 AliAnalysisEt::CreateHistograms();
302
303 TString histname;
304 histname = "fHistChargedPionEnergyDeposit" + fHistogramNameSuffix;
305 fHistChargedPionEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #pi^{+/-}", 1000, 0, 10, 1000, 0, 10);
306 fHistChargedPionEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
307 fHistChargedPionEnergyDeposit->SetYTitle("Energy of track");
308
309 histname = "fHistProtonEnergyDeposit" + fHistogramNameSuffix;
310 fHistProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by protons", 1000, 0, 10, 1000, 0, 10);
311 fHistProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
312 fHistProtonEnergyDeposit->SetYTitle("Energy of track");
313
314 histname = "fHistAntiProtonEnergyDeposit" + fHistogramNameSuffix;
315 fHistAntiProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by anti-protons", 1000, 0, 10, 1000, 0, 10);
316 fHistAntiProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
317 fHistAntiProtonEnergyDeposit->SetYTitle("Energy of track");
318
319 histname = "fHistChargedKaonEnergyDeposit" + fHistogramNameSuffix;
320 fHistChargedKaonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by K^{+/-}", 1000, 0, 10, 1000, 0, 10);
321 fHistChargedKaonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
322 fHistChargedKaonEnergyDeposit->SetYTitle("Energy of track");
323
324 histname = "fHistMuonEnergyDeposit" + fHistogramNameSuffix;
325 fHistMuonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #mu^{+/-}", 1000, 0, 10, 1000, 0, 10);
326 fHistMuonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
327 fHistMuonEnergyDeposit->SetYTitle("Energy of track");
328
329
330}