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