]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/totEt/AliAnalysisEtReconstructed.cxx
Important bug fixes
[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"
ba136eb4 13#include "AliEMCALTrack.h"
2fbf38ac 14#include "AliESDCaloCluster.h"
15#include "TVector3.h"
ba136eb4 16#include "TGeoGlobalMagField.h"
17#include "AliMagF.h"
2fbf38ac 18#include "AliVEvent.h"
19#include "AliESDEvent.h"
20#include "AliVParticle.h"
87efb15c 21#include "TDatabasePDG.h"
22#include "TList.h"
2fbf38ac 23#include <iostream>
24#include "TH2F.h"
25
26AliAnalysisEtReconstructed::AliAnalysisEtReconstructed() :
27 AliAnalysisEt()
87efb15c 28 ,fNTpcClustersCut(0)
29 ,fNItsClustersCut(0)
2fbf38ac 30 ,fTrackDistanceCut(0)
87efb15c 31 ,fPidCut(0)
2fbf38ac 32 ,fClusterType(0)
87efb15c 33 ,fHistChargedPionEnergyDeposit(0)
34 ,fHistProtonEnergyDeposit(0)
35 ,fHistAntiProtonEnergyDeposit(0)
36 ,fHistChargedKaonEnergyDeposit(0)
37 ,fHistMuonEnergyDeposit(0)
2fbf38ac 38{
39
40}
41
cf6522d1 42AliAnalysisEtReconstructed::~AliAnalysisEtReconstructed()
2fbf38ac 43{
cf6522d1 44}
45
46Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
47{ // analyse ESD event
2fbf38ac 48 ResetEventValues();
49 AliESDEvent *event = dynamic_cast<AliESDEvent*>(ev);
50
87efb15c 51 Double_t protonMass = fPdgDB->GetParticle("proton")->Mass();
52
2fbf38ac 53 for (Int_t iTrack = 0; iTrack < event->GetNumberOfTracks(); iTrack++)
54 {
55 AliVParticle *track = event->GetTrack(iTrack);
56 if (!track)
57 {
58 Printf("ERROR: Could not get track %d", iTrack);
59 continue;
60 }
61
62 fMultiplicity++;
63
64 Int_t nItsClusters = dynamic_cast<AliESDtrack*>(track)->GetNcls(0);
65 Int_t nTPCClusters = dynamic_cast<AliESDtrack*>(track)->GetNcls(1);
66
67 Float_t massPart = 0;
68
69 const Double_t *pidWeights = track->PID();
87efb15c 70 Int_t maxpid = -1;
71 Double_t maxpidweight = 0;
72
2fbf38ac 73 if (pidWeights)
74 {
2fbf38ac 75 for (Int_t p =0; p < AliPID::kSPECIES; p++)
76 {
77 if (pidWeights[p] > maxpidweight)
78 {
79 maxpidweight = pidWeights[p];
80 maxpid = p;
81 }
82 }
83 if (maxpid == AliPID::kProton)
84 {
87efb15c 85 //by definition of ET
86 massPart = -protonMass*track->Charge();
2fbf38ac 87 }
88
89 }
90
91 Double_t et = track->E() * TMath::Sin(track->Theta()) + massPart;
cf6522d1 92 // 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 93
4998becf 94 if (TMath::Abs(track->Eta()) < fCuts->GetCommonEtaCut() && CheckGoodVertex(track) && nItsClusters > fCuts->GetReconstructedNItsClustersCut() && nTPCClusters > fCuts->GetReconstructedNTpcClustersCut() )
2fbf38ac 95 {
96 fTotChargedEt += et;
97 fChargedMultiplicity++;
87efb15c 98 if (maxpid != -1)
99 {
100 if (maxpid == AliPID::kPion)
101 {
102 fProtonEt += et;
103 }
104 if (maxpid == AliPID::kKaon)
105 {
106 fChargedKaonEt += et;
107 }
108 if (maxpid == AliPID::kMuon)
109 {
110 fMuonEt += et;
111 }
112 if (maxpid == AliPID::kElectron)
113 {
114 fElectronEt += et;
115 }
116 }
2fbf38ac 117
118 if (TMath::Abs(track->Eta()) < fEtaCutAcc && track->Phi() < fPhiCutAccMax && track->Phi() > fPhiCutAccMin)
119 {
120 fTotChargedEtAcc += track->E()*TMath::Sin(track->Theta()) + massPart;
87efb15c 121 if (maxpid != -1)
122 {
123 if (maxpid == AliPID::kPion)
124 {
125 fProtonEtAcc += et;
126 }
127 if (maxpid == AliPID::kKaon)
128 {
129 fChargedKaonEtAcc += et;
130 }
131 if (maxpid == AliPID::kMuon)
132 {
133 fMuonEtAcc += et;
134 }
135 if (maxpid == AliPID::kElectron)
136 {
137 fElectronEtAcc += et;
138 }
139 }
140
2fbf38ac 141 }
142 }
143
2fbf38ac 144 if (TrackHitsCalorimeter(track, event->GetMagneticField()))
145 {
13b0d3c1 146 Double_t phi = track->Phi();
147 Double_t pt = track->Pt();
cf6522d1 148 // printf("Rec track hit: iTrack %03d phi %4.3f pt %4.3f\n", iTrack, phi, pt); // tmp/debug printout
13b0d3c1 149 if (track->Charge() > 0) fHistPhivsPtPos->Fill(phi, pt);
150 else fHistPhivsPtNeg->Fill(phi, pt);
2fbf38ac 151 }
152 }
153
154 for (Int_t iCluster = 0; iCluster < event->GetNumberOfCaloClusters(); iCluster++)
155 {
156 AliESDCaloCluster* cluster = event->GetCaloCluster(iCluster);
157 if (!cluster)
158 {
159 Printf("ERROR: Could not get cluster %d", iCluster);
160 continue;
161 }
162
87efb15c 163 if (cluster->GetType() != fClusterType) continue;
ba136eb4 164 //if(cluster->GetTracksMatched() > 0)
87efb15c 165 //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
166
167
2fbf38ac 168 if (cluster->E() < fClusterEnergyCut) continue;
169 Float_t pos[3];
170 TVector3 cp(pos);
171 cluster->GetPosition(pos);
87efb15c 172
ba136eb4 173 Double_t distance = cluster->GetEmcCpvDistance();
174 Int_t trackMatchedIndex = cluster->GetTrackMatchedIndex();
175 if ( cluster->IsEMCAL() ) {
176 distance = CalcTrackClusterDistance(pos, &trackMatchedIndex, event);
177 }
178
179 fHistTMDeltaR->Fill(distance);
180 if (distance < fTrackDistanceCut)
2fbf38ac 181 {
ba136eb4 182 if (cluster->GetNTracksMatched() == 1 && trackMatchedIndex>0)
87efb15c 183 {
ba136eb4 184 AliVTrack *track = event->GetTrack(trackMatchedIndex);
87efb15c 185 const Double_t *pidWeights = track->PID();
186
187 Double_t maxpidweight = 0;
188 Int_t maxpid = 0;
189
190 if (pidWeights)
191 {
192 for (Int_t p =0; p < AliPID::kSPECIES; p++)
193 {
194 if (pidWeights[p] > maxpidweight)
195 {
196 maxpidweight = pidWeights[p];
197 maxpid = p;
198 }
199 }
200 if(maxpidweight > fPidCut)
201 {
202 if(maxpid == AliPID::kProton)
203 {
204 if(track->Charge() == 1)
205 {
206 fHistProtonEnergyDeposit->Fill(cluster->E(), track->E());
207 }
208 else if(track->Charge() == -1)
209 {
210 fHistAntiProtonEnergyDeposit->Fill(cluster->E(), track->E());
211 }
212 }
213 else if(maxpid == AliPID::kPion)
214 {
215 fHistChargedPionEnergyDeposit->Fill(cluster->E(), track->E());
216 }
217 else if(maxpid == AliPID::kKaon)
218 {
219 fHistChargedKaonEnergyDeposit->Fill(cluster->E(), track->E());
220 }
221 else if(maxpid == AliPID::kMuon)
222 {
223 fHistMuonEnergyDeposit->Fill(cluster->E(), track->E());
224 }
225 }
226 }
227 }
228
2fbf38ac 229 continue;
ba136eb4 230 } // distance
2fbf38ac 231
4998becf 232 if (cluster->E() > fSingleCellEnergyCut && cluster->GetNCells() == fCuts->GetCommonSingleCell()) continue;
2fbf38ac 233
234 cluster->GetPosition(pos);
235
236 // TODO: replace with TVector3, too lazy now...
237
238 float dist = TMath::Sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
239
240 float theta = TMath::ATan(pos[2]/dist)+TMath::Pi()/2;
241 // float eta = TMath::Log(TMath::Abs( TMath::Tan( 0.5 * theta ) ) );
242 fTotNeutralEt += cluster->E() * TMath::Sin(theta);
13b0d3c1 243 fNeutralMultiplicity++;
2fbf38ac 244
245 fMultiplicity++;
2fbf38ac 246 }
247
248 fTotNeutralEtAcc = fTotNeutralEt;
249 fTotEt = fTotChargedEt + fTotNeutralEt;
250 fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;
251
2fbf38ac 252 // Fill the histograms...
253 FillHistograms();
254
255 return 0;
256}
257
258bool AliAnalysisEtReconstructed::CheckGoodVertex(AliVParticle* track)
cf6522d1 259{ // check vertex
2fbf38ac 260
261 Float_t bxy = 999.;
262 Float_t bz = 999.;
263 dynamic_cast<AliESDtrack*>(track)->GetImpactParametersTPC(bxy,bz);
264
4998becf 265 bool status = (TMath::Abs(track->Xv()) < fCuts->GetReconstructedVertexXCut()) &&
266 (TMath::Abs(track->Yv()) < fCuts->GetReconstructedVertexYCut()) &&
267 (TMath::Abs(track->Zv()) < fCuts->GetReconstructedVertexZCut()) &&
268 (TMath::Abs(bxy) < fCuts->GetReconstructedIPxyCut()) &&
269 (TMath::Abs(bz) < fCuts->GetReconstructedIPzCut());
2fbf38ac 270
4998becf 271 return status;
2fbf38ac 272}
273
274void AliAnalysisEtReconstructed::Init()
cf6522d1 275{ // Init
2fbf38ac 276 AliAnalysisEt::Init();
87efb15c 277 fNItsClustersCut = fCuts->GetReconstructedNItsClustersCut();
278 fNTpcClustersCut = fCuts->GetReconstructedNTpcClustersCut();
83d0f02c 279 fPidCut = fCuts->GetReconstructedPidCut();
ba136eb4 280 TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
281
2fbf38ac 282}
283
284bool AliAnalysisEtReconstructed::TrackHitsCalorimeter(AliVParticle* track, Double_t magField)
cf6522d1 285{ // propagate track to detector radius
2fbf38ac 286
287 AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
cf6522d1 288 // Printf("Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
2fbf38ac 289
290 Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);
291
cf6522d1 292 // if (prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
2fbf38ac 293 return prop &&
294 TMath::Abs(esdTrack->Eta()) < fEtaCutAcc &&
295 esdTrack->Phi() > fPhiCutAccMin*TMath::Pi()/180. &&
296 esdTrack->Phi() < fPhiCutAccMax*TMath::Pi()/180.;
297}
298
87efb15c 299void AliAnalysisEtReconstructed::FillOutputList(TList* list)
300{
301 AliAnalysisEt::FillOutputList(list);
302
303 list->Add(fHistChargedPionEnergyDeposit);
304 list->Add(fHistProtonEnergyDeposit);
305 list->Add(fHistAntiProtonEnergyDeposit);
306 list->Add(fHistChargedKaonEnergyDeposit);
307 list->Add(fHistMuonEnergyDeposit);
308}
309
310void AliAnalysisEtReconstructed::CreateHistograms()
311{
312 AliAnalysisEt::CreateHistograms();
313
314 TString histname;
315 histname = "fHistChargedPionEnergyDeposit" + fHistogramNameSuffix;
316 fHistChargedPionEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #pi^{+/-}", 1000, 0, 10, 1000, 0, 10);
317 fHistChargedPionEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
318 fHistChargedPionEnergyDeposit->SetYTitle("Energy of track");
319
320 histname = "fHistProtonEnergyDeposit" + fHistogramNameSuffix;
321 fHistProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by protons", 1000, 0, 10, 1000, 0, 10);
322 fHistProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
323 fHistProtonEnergyDeposit->SetYTitle("Energy of track");
324
325 histname = "fHistAntiProtonEnergyDeposit" + fHistogramNameSuffix;
326 fHistAntiProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by anti-protons", 1000, 0, 10, 1000, 0, 10);
327 fHistAntiProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
328 fHistAntiProtonEnergyDeposit->SetYTitle("Energy of track");
329
330 histname = "fHistChargedKaonEnergyDeposit" + fHistogramNameSuffix;
331 fHistChargedKaonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by K^{+/-}", 1000, 0, 10, 1000, 0, 10);
332 fHistChargedKaonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
333 fHistChargedKaonEnergyDeposit->SetYTitle("Energy of track");
334
335 histname = "fHistMuonEnergyDeposit" + fHistogramNameSuffix;
336 fHistMuonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #mu^{+/-}", 1000, 0, 10, 1000, 0, 10);
337 fHistMuonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
338 fHistMuonEnergyDeposit->SetYTitle("Energy of track");
339
340
341}
ba136eb4 342
343Double_t
344AliAnalysisEtReconstructed::CalcTrackClusterDistance(Float_t clsPos[3],
345 Int_t *trkMatchId,
346 AliESDEvent *event)
347{ // calculate distance between cluster and closest track
348
349 Double_t trkPos[3] = {0,0,0};
350
351 Int_t bestTrkMatchId = -1;
352 Double_t distance = 9999; // init to a big number
353
354 Double_t dist = 0;
355 Double_t distX = 0, distY = 0, distZ = 0;
356
357 for (Int_t iTrack = 0; iTrack < event->GetNumberOfTracks(); iTrack++) {
358 AliESDtrack *track = event->GetTrack(iTrack);
359 if (!track) {
360 Printf("ERROR: Could not get track %d", iTrack);
361 continue;
362 }
363
364 // check for approx. eta and phi range before we propagate..
365 // TBD
366
367 AliEMCALTrack *emctrack = new AliEMCALTrack(*track);
368 if (!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) ) {
369 continue;
370 }
371 emctrack->GetXYZ(trkPos);
372 if (emctrack) delete emctrack;
373
374 distX = clsPos[0]-trkPos[0];
375 distY = clsPos[1]-trkPos[1];
376 distZ = clsPos[2]-trkPos[2];
377 dist = TMath::Sqrt(distX*distX + distY*distY + distZ*distZ);
378
379 if (dist < distance) {
380 distance = dist;
381 bestTrkMatchId = iTrack;
382 }
383 } // iTrack
384
385 // printf("CalcTrackClusterDistance: bestTrkMatch %d origTrkMatch %d distance %f\n", bestTrkMatchId, *trkMatchId, distance);
386 *trkMatchId = bestTrkMatchId;
387 return distance;
388}