- adding histograms for energy deposited by identified charged particles
[u/mrichter/AliRoot.git] / PWG4 / totEt / AliAnalysisEtReconstructed.cxx
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 //_________________________________________________________________________
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"
18 #include "TDatabasePDG.h"
19 #include "TList.h"
20 #include <iostream>
21 #include "TH2F.h"
22
23 AliAnalysisEtReconstructed::AliAnalysisEtReconstructed() :
24         AliAnalysisEt()
25         ,fNTpcClustersCut(0)
26         ,fNItsClustersCut(0)
27         ,fTrackDistanceCut(0)
28         ,fPidCut(0)
29         ,fClusterType(0)
30         ,fHistChargedPionEnergyDeposit(0)
31         ,fHistProtonEnergyDeposit(0)
32         ,fHistAntiProtonEnergyDeposit(0)
33         ,fHistChargedKaonEnergyDeposit(0)
34         ,fHistMuonEnergyDeposit(0)
35 {
36
37 }
38
39 AliAnalysisEtReconstructed::~AliAnalysisEtReconstructed() 
40 {
41 }
42
43 Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
44 { // analyse ESD event
45     ResetEventValues();
46     AliESDEvent *event = dynamic_cast<AliESDEvent*>(ev);
47
48     Double_t protonMass = fPdgDB->GetParticle("proton")->Mass();
49
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();
67         Int_t maxpid = -1;
68         Double_t maxpidweight = 0;
69             
70         if (pidWeights)
71         {
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             {
82               //by definition of ET
83                 massPart = -protonMass*track->Charge();
84             }
85
86         }
87
88         Double_t et = track->E() * TMath::Sin(track->Theta()) + massPart;
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
90
91         if (TMath::Abs(track->Eta()) < fCuts->GetCommonEtaCut() && CheckGoodVertex(track) && nItsClusters > fCuts->GetReconstructedNItsClustersCut() && nTPCClusters > fCuts->GetReconstructedNTpcClustersCut() )
92         {
93             fTotChargedEt +=  et;
94             fChargedMultiplicity++;
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             }
114
115             if (TMath::Abs(track->Eta()) < fEtaCutAcc && track->Phi() < fPhiCutAccMax && track->Phi() > fPhiCutAccMin)
116             {
117                 fTotChargedEtAcc += track->E()*TMath::Sin(track->Theta()) + massPart;
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            
138             }
139         }
140
141         if (TrackHitsCalorimeter(track, event->GetMagneticField()))
142         {
143           Double_t phi = track->Phi();
144           Double_t pt = track->Pt();
145           // printf("Rec track hit: iTrack %03d phi %4.3f pt %4.3f\n", iTrack, phi, pt); // tmp/debug printout
146           if (track->Charge() > 0) fHistPhivsPtPos->Fill(phi, pt);
147           else fHistPhivsPtNeg->Fill(phi, pt);
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
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         
165         if (cluster->E() < fClusterEnergyCut) continue;
166         Float_t pos[3];
167         TVector3 cp(pos);
168         cluster->GetPosition(pos);
169
170         fHistTMDeltaR->Fill(cluster->GetEmcCpvDistance());
171         if (cluster->GetEmcCpvDistance() < fTrackDistanceCut)
172         {
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
220             continue;
221         }
222
223         if (cluster->E() >  fSingleCellEnergyCut && cluster->GetNCells() == fCuts->GetCommonSingleCell()) continue;
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);
234         fNeutralMultiplicity++;
235
236         fMultiplicity++;
237     }
238
239     fTotNeutralEtAcc = fTotNeutralEt;
240     fTotEt = fTotChargedEt + fTotNeutralEt;
241     fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc;
242
243     // Fill the histograms...
244     FillHistograms();
245
246     return 0;
247 }
248
249 bool AliAnalysisEtReconstructed::CheckGoodVertex(AliVParticle* track)
250 { // check vertex
251
252     Float_t bxy = 999.;
253     Float_t bz = 999.;
254     dynamic_cast<AliESDtrack*>(track)->GetImpactParametersTPC(bxy,bz);
255
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()); 
261
262     return status;
263 }
264
265 void AliAnalysisEtReconstructed::Init()
266 { // Init
267     AliAnalysisEt::Init();
268     fNItsClustersCut = fCuts->GetReconstructedNItsClustersCut();
269     fNTpcClustersCut = fCuts->GetReconstructedNTpcClustersCut();
270     fPidCut = fCuts->GetCommonPidCut();
271 }
272
273 bool AliAnalysisEtReconstructed::TrackHitsCalorimeter(AliVParticle* track, Double_t magField)
274 { // propagate track to detector radius
275
276    AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
277    // Printf("Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
278
279     Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);
280
281     // if (prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
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
288 void 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
299 void 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 }