setter to assume pion mass for clusters
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliAnalysisTaskEmcalJetSample.cxx
1 // $Id$
2 //
3 // Jet sample analysis task.
4 //
5 // Author: S.Aiola
6
7 #include <TClonesArray.h>
8 #include <TH1F.h>
9 #include <TH2F.h>
10 #include <TH3F.h>
11 #include <TList.h>
12 #include <TLorentzVector.h>
13
14 #include "AliVCluster.h"
15 #include "AliAODCaloCluster.h"
16 #include "AliESDCaloCluster.h"
17 #include "AliVTrack.h"
18 #include "AliEmcalJet.h"
19 #include "AliRhoParameter.h"
20 #include "AliLog.h"
21 #include "AliJetContainer.h"
22 #include "AliParticleContainer.h"
23 #include "AliClusterContainer.h"
24 #include "AliPicoTrack.h"
25
26 #include "AliAnalysisTaskEmcalJetSample.h"
27
28 ClassImp(AliAnalysisTaskEmcalJetSample)
29
30 //________________________________________________________________________
31 AliAnalysisTaskEmcalJetSample::AliAnalysisTaskEmcalJetSample() : 
32   AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetSample", kTRUE),
33   fHistTracksPt(0),
34   fHistClustersPt(0),
35   fHistLeadingJetPt(0),
36   fHistJetsPhiEta(0),
37   fHistJetsPtArea(0),
38   fHistJetsPtLeadHad(0),
39   fHistJetsCorrPtArea(0),
40   fHistPtDEtaDPhiTrackClus(0),
41   fHistPtDEtaDPhiClusTrack(0),
42   fJetsCont(0),
43   fTracksCont(0),
44   fCaloClustersCont(0)
45
46 {
47   // Default constructor.
48
49   fHistTracksPt       = new TH1*[fNcentBins];
50   fHistClustersPt     = new TH1*[fNcentBins];
51   fHistLeadingJetPt   = new TH1*[fNcentBins];
52   fHistJetsPhiEta     = new TH2*[fNcentBins];
53   fHistJetsPtArea     = new TH2*[fNcentBins];
54   fHistJetsPtLeadHad  = new TH2*[fNcentBins];
55   fHistJetsCorrPtArea = new TH2*[fNcentBins];
56
57   for (Int_t i = 0; i < fNcentBins; i++) {
58     fHistTracksPt[i] = 0;
59     fHistClustersPt[i] = 0;
60     fHistLeadingJetPt[i] = 0;
61     fHistJetsPhiEta[i] = 0;
62     fHistJetsPtArea[i] = 0;
63     fHistJetsPtLeadHad[i] = 0;
64     fHistJetsCorrPtArea[i] = 0;
65   }
66
67   SetMakeGeneralHistograms(kTRUE);
68 }
69
70 //________________________________________________________________________
71 AliAnalysisTaskEmcalJetSample::AliAnalysisTaskEmcalJetSample(const char *name) : 
72   AliAnalysisTaskEmcalJet(name, kTRUE),
73   fHistTracksPt(0),
74   fHistClustersPt(0),
75   fHistLeadingJetPt(0),
76   fHistJetsPhiEta(0),
77   fHistJetsPtArea(0),
78   fHistJetsPtLeadHad(0),
79   fHistJetsCorrPtArea(0),
80   fHistPtDEtaDPhiTrackClus(0),
81   fHistPtDEtaDPhiClusTrack(0),
82   fJetsCont(0),
83   fTracksCont(0),
84   fCaloClustersCont(0)
85 {
86   // Standard constructor.
87
88   fHistTracksPt       = new TH1*[fNcentBins];
89   fHistClustersPt     = new TH1*[fNcentBins];
90   fHistLeadingJetPt   = new TH1*[fNcentBins];
91   fHistJetsPhiEta     = new TH2*[fNcentBins];
92   fHistJetsPtArea     = new TH2*[fNcentBins];
93   fHistJetsPtLeadHad  = new TH2*[fNcentBins];
94   fHistJetsCorrPtArea = new TH2*[fNcentBins];
95
96   for (Int_t i = 0; i < fNcentBins; i++) {
97     fHistTracksPt[i] = 0;
98     fHistClustersPt[i] = 0;
99     fHistLeadingJetPt[i] = 0;
100     fHistJetsPhiEta[i] = 0;
101     fHistJetsPtArea[i] = 0;
102     fHistJetsPtLeadHad[i] = 0;
103     fHistJetsCorrPtArea[i] = 0;
104   }
105
106   SetMakeGeneralHistograms(kTRUE);
107 }
108
109 //________________________________________________________________________
110 AliAnalysisTaskEmcalJetSample::~AliAnalysisTaskEmcalJetSample()
111 {
112   // Destructor.
113 }
114
115 //________________________________________________________________________
116 void AliAnalysisTaskEmcalJetSample::UserCreateOutputObjects()
117 {
118   // Create user output.
119
120   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
121
122   fJetsCont           = GetJetContainer(0);
123   if(fJetsCont) { //get particles and clusters connected to jets
124     fTracksCont       = fJetsCont->GetParticleContainer();
125     fCaloClustersCont = fJetsCont->GetClusterContainer();
126   } else {        //no jets, just analysis tracks and clusters
127     fTracksCont       = GetParticleContainer(0);
128     fCaloClustersCont = GetClusterContainer(0);
129   }
130   fTracksCont->SetClassName("AliVTrack");
131   fCaloClustersCont->SetClassName("AliVCluster");
132
133   TString histname;
134
135   for (Int_t i = 0; i < fNcentBins; i++) {
136     if (fParticleCollArray.GetEntriesFast()>0) {
137       histname = "fHistTracksPt_";
138       histname += i;
139       fHistTracksPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2);
140       fHistTracksPt[i]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
141       fHistTracksPt[i]->GetYaxis()->SetTitle("counts");
142       fOutput->Add(fHistTracksPt[i]);
143     }
144
145     if (fClusterCollArray.GetEntriesFast()>0) {
146       histname = "fHistClustersPt_";
147       histname += i;
148       fHistClustersPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2);
149       fHistClustersPt[i]->GetXaxis()->SetTitle("p_{T,clus} (GeV/c)");
150       fHistClustersPt[i]->GetYaxis()->SetTitle("counts");
151       fOutput->Add(fHistClustersPt[i]);
152     }
153
154     if (fJetCollArray.GetEntriesFast()>0) {
155       histname = "fHistLeadingJetPt_";
156       histname += i;
157       fHistLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
158       fHistLeadingJetPt[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
159       fHistLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
160       fOutput->Add(fHistLeadingJetPt[i]);
161       
162       histname = "fHistJetsPhiEta_";
163       histname += i;
164       fHistJetsPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 50, -1, 1, 101, 0, TMath::Pi()*2 + TMath::Pi()/200);
165       fHistJetsPhiEta[i]->GetXaxis()->SetTitle("#eta");
166       fHistJetsPhiEta[i]->GetYaxis()->SetTitle("#phi");
167       fOutput->Add(fHistJetsPhiEta[i]);
168       
169       histname = "fHistJetsPtArea_";
170       histname += i;
171       fHistJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, 30, 0, 3);
172       fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
173       fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
174       fOutput->Add(fHistJetsPtArea[i]);
175
176       histname = "fHistJetsPtLeadHad_";
177       histname += i;
178       fHistJetsPtLeadHad[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
179       fHistJetsPtLeadHad[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
180       fHistJetsPtLeadHad[i]->GetYaxis()->SetTitle("p_{T,lead} (GeV/c)");
181       fHistJetsPtLeadHad[i]->GetZaxis()->SetTitle("counts");
182       fOutput->Add(fHistJetsPtLeadHad[i]);
183     
184       if (!(GetJetContainer()->GetRhoName().IsNull())) {
185         histname = "fHistJetsCorrPtArea_";
186         histname += i;
187         fHistJetsCorrPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins*2, -fMaxBinPt, fMaxBinPt, 30, 0, 3);
188         fHistJetsCorrPtArea[i]->GetXaxis()->SetTitle("p_{T}^{corr} [GeV/c]");
189         fHistJetsCorrPtArea[i]->GetYaxis()->SetTitle("area");
190         fOutput->Add(fHistJetsCorrPtArea[i]);
191       }
192     }
193   }
194
195   histname = "fHistPtDEtaDPhiTrackClus";
196   fHistPtDEtaDPhiTrackClus = new TH3F(histname.Data(),Form("%s;#it{p}_{T}^{track};#Delta#eta;#Delta#varphi",histname.Data()),100,0.,100.,100,-0.1,0.1,100,-0.1,0.1);
197   fOutput->Add(fHistPtDEtaDPhiTrackClus);
198
199   histname = "fHistPtDEtaDPhiClusTrack";
200   fHistPtDEtaDPhiClusTrack = new TH3F(histname.Data(),Form("%s;#it{p}_{T}^{clus};#Delta#eta;#Delta#varphi",histname.Data()),100,0.,100.,100,-0.1,0.1,100,-0.1,0.1);
201   fOutput->Add(fHistPtDEtaDPhiClusTrack);
202
203   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
204 }
205
206 //________________________________________________________________________
207 Bool_t AliAnalysisTaskEmcalJetSample::FillHistograms()
208 {
209   // Fill histograms.
210
211   if (fTracksCont) {
212     AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle(0)); 
213     while(track) {
214       fHistTracksPt[fCentBin]->Fill(track->Pt()); 
215       track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
216     }
217   }
218   
219   if (fCaloClustersCont) {
220     AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster(0); 
221     while(cluster) {
222       TLorentzVector nPart;
223       cluster->GetMomentum(nPart, fVertex);
224       fHistClustersPt[fCentBin]->Fill(nPart.Pt());
225
226       cluster = fCaloClustersCont->GetNextAcceptCluster();
227     }
228   }
229
230   if (fJetsCont) {
231     AliEmcalJet *jet = fJetsCont->GetNextAcceptJet(0); 
232     while(jet) {
233
234       fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());
235       fHistJetsPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi());
236
237       Float_t ptLeading = fJetsCont->GetLeadingHadronPt(jet);
238       fHistJetsPtLeadHad[fCentBin]->Fill(jet->Pt(), ptLeading);
239
240       if (fHistJetsCorrPtArea[fCentBin]) {
241         Float_t corrPt = jet->Pt() - fJetsCont->GetRhoVal() * jet->Area();
242         fHistJetsCorrPtArea[fCentBin]->Fill(corrPt, jet->Area());
243       }
244       jet = fJetsCont->GetNextAcceptJet(); 
245     }
246     
247     jet = fJetsCont->GetLeadingJet();
248     if(jet) fHistLeadingJetPt[fCentBin]->Fill(jet->Pt());
249   }
250
251   CheckClusTrackMatching();
252
253   return kTRUE;
254 }
255
256 //________________________________________________________________________
257 void AliAnalysisTaskEmcalJetSample::CheckClusTrackMatching()
258 {
259   
260   if(!fTracksCont || !fCaloClustersCont)
261     return;
262
263   Double_t deta = 999;
264   Double_t dphi = 999;
265
266   //Get closest cluster to track
267   AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle(0)); 
268   while(track) {
269     //Get matched cluster
270     Int_t emc1 = track->GetEMCALcluster();
271     if(fCaloClustersCont && emc1>=0) {
272       AliVCluster *clusMatch = fCaloClustersCont->GetCluster(emc1);
273       if(clusMatch) {
274         AliPicoTrack::GetEtaPhiDiff(track, clusMatch, dphi, deta);
275         fHistPtDEtaDPhiTrackClus->Fill(track->Pt(),deta,dphi);
276       }
277     }
278     track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
279   }
280   
281   //Get closest track to cluster
282   AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster(0); 
283   while(cluster) {
284     TLorentzVector nPart;
285     cluster->GetMomentum(nPart, fVertex);
286     fHistClustersPt[fCentBin]->Fill(nPart.Pt());
287     
288     //Get matched track
289     AliVTrack *mt = NULL;      
290     AliAODCaloCluster *acl = dynamic_cast<AliAODCaloCluster*>(cluster);
291     if(acl) {
292       if(acl->GetNTracksMatched()>1)
293         mt = static_cast<AliVTrack*>(acl->GetTrackMatched(0));
294     }
295     else {
296       AliESDCaloCluster *ecl = dynamic_cast<AliESDCaloCluster*>(cluster);
297       Int_t im = ecl->GetTrackMatchedIndex();
298       if(fTracksCont && im>=0) {
299         mt = static_cast<AliVTrack*>(fTracksCont->GetParticle(im));
300       }
301     }
302     if(mt) {
303       AliPicoTrack::GetEtaPhiDiff(mt, cluster, dphi, deta);
304       fHistPtDEtaDPhiClusTrack->Fill(nPart.Pt(),deta,dphi);
305       
306       /* //debugging
307          if(mt->IsEMCAL()) {
308          Int_t emc1 = mt->GetEMCALcluster();
309          Printf("current id: %d  emc1: %d",fCaloClustersCont->GetCurrentID(),emc1);
310          AliVCluster *clm = fCaloClustersCont->GetCluster(emc1);
311          AliPicoTrack::GetEtaPhiDiff(mt, clm, dphi, deta);
312          Printf("deta: %f dphi: %f",deta,dphi);
313          }
314       */
315     }
316     cluster = fCaloClustersCont->GetNextAcceptCluster();
317   }
318 }
319
320 //________________________________________________________________________
321 void AliAnalysisTaskEmcalJetSample::ExecOnce() {
322
323   AliAnalysisTaskEmcalJet::ExecOnce();
324
325   if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
326   if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
327   if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
328
329 }
330
331 //________________________________________________________________________
332 Bool_t AliAnalysisTaskEmcalJetSample::Run()
333 {
334   // Run analysis code here, if needed. It will be executed before FillHistograms().
335
336   return kTRUE;  // If return kFALSE FillHistogram() will NOT be executed.
337 }
338
339 //________________________________________________________________________
340 void AliAnalysisTaskEmcalJetSample::Terminate(Option_t *) 
341 {
342   // Called once at the end of the analysis.
343 }