AliPhysics  e6c8d43 (e6c8d43)
AliAnalysisTaskDcalDijetPerf.cxx
Go to the documentation of this file.
1 // $Id$
2 //
3 // Dcal dijet performance task
4 //
5 // Author: R. Reed
6 
7 #include <TClonesArray.h>
8 #include <TH1F.h>
9 #include <TH2F.h>
10 #include <TH3F.h>
11 #include <THnSparse.h>
12 #include <TList.h>
13 #include <TLorentzVector.h>
14 
15 #include "AliVCluster.h"
16 #include "AliAODCaloCluster.h"
17 #include "AliESDCaloCluster.h"
18 #include "AliVTrack.h"
19 #include "AliEmcalJet.h"
20 #include "AliRhoParameter.h"
21 #include "AliLog.h"
22 #include "AliJetContainer.h"
23 #include "AliParticleContainer.h"
24 #include "AliClusterContainer.h"
25 #include "AliPicoTrack.h"
26 
28 
30 
31 //________________________________________________________________________
34  fHistTracksPt(0),
35  fHistTracksEtaPhi(0),
36  fHistClustersPt(0),
37  fHistLeadingJetPt(0),
38  fHistJetsPhiEta(0),
39  fHistJetsPtArea(0),
40  fHistJetsPtLeadHad(0),
41  fHistJetsCorrPtArea(0),
42  fHistJet1(0),
43  fHistJet1m(0),
44  fHistJet1nm(0),
45  fHistJet2(0),
46  fHistJet1to2(0),
47  fHistDiJet1(0),
48  fHistDiJet1m(0),
49  fJetsCont(0),
50  fJetsCont2(0),
51  fJetsCont3(0),
52  fTracksCont(0),
53  fCaloClustersCont(0)
54 
55 {
56  // Default constructor.
57 
58  fHistTracksPt = new TH1*[fNcentBins];
59  fHistTracksEtaPhi = new TH2*[fNcentBins];
60  fHistClustersPt = new TH1*[fNcentBins];
61  fHistLeadingJetPt = new TH1*[fNcentBins];
62  fHistJetsPhiEta = new TH2*[fNcentBins];
63  fHistJetsPtArea = new TH2*[fNcentBins];
64  fHistJetsPtLeadHad = new TH2*[fNcentBins];
65  fHistJetsCorrPtArea = new TH2*[fNcentBins];
66 
67  for (Int_t i = 0; i < fNcentBins; i++) {
68  fHistTracksPt[i] = 0;
69  fHistTracksEtaPhi[i] = 0;
70  fHistClustersPt[i] = 0;
71  fHistLeadingJetPt[i] = 0;
72  fHistJetsPhiEta[i] = 0;
73  fHistJetsPtArea[i] = 0;
74  fHistJetsPtLeadHad[i] = 0;
75  fHistJetsCorrPtArea[i] = 0;
76  }
77  fHistJet1 = 0;
78  fHistJet1m = 0;
79  fHistJet1nm = 0;
80  fHistJet2 = 0;
81  fHistJet1to2 = 0;
82  fHistDiJet1 = 0;
83  fHistDiJet1m = 0;
84  SetMakeGeneralHistograms(kTRUE);
85 }
86 
87 //________________________________________________________________________
89  AliAnalysisTaskEmcalJet(name, kTRUE),
90  fHistTracksPt(0),
91  fHistTracksEtaPhi(0),
92  fHistClustersPt(0),
93  fHistLeadingJetPt(0),
94  fHistJetsPhiEta(0),
95  fHistJetsPtArea(0),
96  fHistJetsPtLeadHad(0),
97  fHistJetsCorrPtArea(0),
98  fHistJet1(0),
99  fHistJet1m(0),
100  fHistJet1nm(0),
101  fHistJet2(0),
102  fHistJet1to2(0),
103  fHistDiJet1(0),
104  fHistDiJet1m(0),
105  fJetsCont(0),
106  fJetsCont2(0),
107  fJetsCont3(0),
108  fTracksCont(0),
109  fCaloClustersCont(0)
110 {
111  // Standard constructor.
112 
113  fHistTracksPt = new TH1*[fNcentBins];
121 
122  for (Int_t i = 0; i < fNcentBins; i++) {
123  fHistTracksPt[i] = 0;
124  fHistTracksEtaPhi[i] = 0;
125  fHistClustersPt[i] = 0;
126  fHistLeadingJetPt[i] = 0;
127  fHistJetsPhiEta[i] = 0;
128  fHistJetsPtArea[i] = 0;
129  fHistJetsPtLeadHad[i] = 0;
130  fHistJetsCorrPtArea[i] = 0;
131  }
132  fHistJet1 = 0;
133  fHistJet1m = 0;
134  fHistJet1nm = 0;
135  fHistJet2 = 0;
136  fHistJet1to2 = 0;
137  fHistDiJet1 = 0;
138  fHistDiJet1m = 0;
139 
141 }
142 
143 //________________________________________________________________________
145 {
146  // Destructor.
147 }
148 
149 //________________________________________________________________________
151 {
152  // Create user output.
153 
155 
159 
160  if(fJetsCont) { //get particles and clusters connected to jets
163  } else { //no jets, just analysis tracks and clusters
166  }
167  fTracksCont->SetClassName("AliVTrack");
168  //fCaloClustersCont->SetClassName("AliVCluster");
169  TString histname;
170 
171  for (Int_t i = 0; i < fNcentBins; i++) {
172  if (fParticleCollArray.GetEntriesFast()>0) {
173  histname = "fHistTracksPt_";
174  histname += i;
175  fHistTracksPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2);
176  fHistTracksPt[i]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
177  fHistTracksPt[i]->GetYaxis()->SetTitle("counts");
178  fOutput->Add(fHistTracksPt[i]);
179  histname = "fHistTracksEtaPhi_";
180  histname += i;
181  fHistTracksEtaPhi[i] = new TH2F(histname.Data(), histname.Data(), fNbins, -0.7, 0.7, fNbins, 0, TMath::Pi()*2);
182  fHistTracksEtaPhi[i]->GetXaxis()->SetTitle("eta");
183  fHistTracksEtaPhi[i]->GetYaxis()->SetTitle("phi");
184  fOutput->Add(fHistTracksEtaPhi[i]);
185  }
186 
187  if (fClusterCollArray.GetEntriesFast()>0) {
188  histname = "fHistClustersPt_";
189  histname += i;
190  fHistClustersPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2);
191  fHistClustersPt[i]->GetXaxis()->SetTitle("p_{T,clus} (GeV/c)");
192  fHistClustersPt[i]->GetYaxis()->SetTitle("counts");
193  fOutput->Add(fHistClustersPt[i]);
194  }
195 
196  if (fJetCollArray.GetEntriesFast()>0) {
197  histname = "fHistLeadingJetPt_";
198  histname += i;
199  fHistLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
200  fHistLeadingJetPt[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
201  fHistLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
202  fOutput->Add(fHistLeadingJetPt[i]);
203 
204  histname = "fHistJetsPhiEta_";
205  histname += i;
206  fHistJetsPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 50, -1, 1, 101, 0, TMath::Pi()*2 + TMath::Pi()/200);
207  fHistJetsPhiEta[i]->GetXaxis()->SetTitle("#eta");
208  fHistJetsPhiEta[i]->GetYaxis()->SetTitle("#phi");
209  fOutput->Add(fHistJetsPhiEta[i]);
210 
211  histname = "fHistJetsPtArea_";
212  histname += i;
213  fHistJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, 30, 0, 3);
214  fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
215  fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
216  fOutput->Add(fHistJetsPtArea[i]);
217 
218  histname = "fHistJetsPtLeadHad_";
219  histname += i;
220  fHistJetsPtLeadHad[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
221  fHistJetsPtLeadHad[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
222  fHistJetsPtLeadHad[i]->GetYaxis()->SetTitle("p_{T,lead} (GeV/c)");
223  fHistJetsPtLeadHad[i]->GetZaxis()->SetTitle("counts");
224  fOutput->Add(fHistJetsPtLeadHad[i]);
225 
226  if (!(GetJetContainer()->GetRhoName().IsNull())) {
227  histname = "fHistJetsCorrPtArea_";
228  histname += i;
229  fHistJetsCorrPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins*2, -fMaxBinPt, fMaxBinPt, 30, 0, 3);
230  fHistJetsCorrPtArea[i]->GetXaxis()->SetTitle("p_{T}^{corr} [GeV/c]");
231  fHistJetsCorrPtArea[i]->GetYaxis()->SetTitle("area");
232  fOutput->Add(fHistJetsCorrPtArea[i]);
233  }
234  }
235  }
236 
237  Int_t nbins[] = {150,100,100,100,150};
238  Double_t xmin[] = { 0,-0.7, 0, 0, 0};
239  Double_t xmax[] = {150, 0.7,TMath::TwoPi(), 1, 150};
240  fHistJet1 = new THnSparseF("Jets1Collection","Jets1Collection",5,nbins,xmin,xmax);
241  fOutput->Add(fHistJet1);
242  fHistJet1m = new THnSparseF("Jets1CollectionMatched","Jets1Collection",5,nbins,xmin,xmax);
243  fOutput->Add(fHistJet1m);
244  fHistJet1nm = new THnSparseF("Jets1CollectionNotMatched","Jets1Collection",5,nbins,xmin,xmax);
245  fOutput->Add(fHistJet1nm);
246  fHistJet2 = new THnSparseF("Jets2Collection","Jets2Collection",5,nbins,xmin,xmax);
247  fOutput->Add(fHistJet2);
248 
249  Int_t nbins2[] = {150,100,100,100,150,100,100,100,100};
250  Double_t xmin2[] = {0,-0.7,0,0,0,-0.7,0,0,0};
251  Double_t xmax2[] = {150,0.7,6.28,1,150,0.7,6.28,1,0.2};
252  fHistJet1to2 = new THnSparseF("Jets1to2Collection","Jets1to2Collection",9,nbins2,xmin2,xmax2);
253  fOutput->Add(fHistJet1to2);
254 
255  Int_t nbins3[] = {150,100,100,100,150,100,100,100,100};
256  Double_t xmin3[] = {0,-0.7,0,0,0,-0.7,0,0,0};
257  Double_t xmax3[] = {150,0.7,6.28,1,150,0.7,6.28,1,1};
258  fHistDiJet1 = new THnSparseF("fHistDiJet1","fHistDiJet1",9,nbins3,xmin3,xmax3);
259  fOutput->Add(fHistDiJet1);
260 
261  Int_t nbins4[] = {150,100,100,100,150,100,100,100,100,150,100,100,100};
262  Double_t xmin4[] = {0,-0.7,0,0,0,-0.7,0,0,0,0,-0.7,0,0}; //pt1 eta1 phi1 NEF1 pt2 eta2 phi2 NEF2 AJ pt3 eta3 phi3 R
263  Double_t xmax4[] = {150,0.7,6.28,1,150,0.7,6.28,1,1,150,0.7,6.28,0.2};
264  fHistDiJet1m = new THnSparseF("fHistDiJet1m","fHistDiJet1m",13,nbins4,xmin4,xmax4);
265  fOutput->Add(fHistDiJet1m);
266 
267 
268 
269  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
270 }
271 
272 //________________________________________________________________________
274 {
275  // Fill histograms.
276 
277  if (fTracksCont) {
278  fTracksCont->ResetCurrentID();
279  AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
280  while(track) {
281  fHistTracksPt[fCentBin]->Fill(track->Pt());
282  fHistTracksEtaPhi[fCentBin]->Fill(track->Eta(),track->Phi());
283  track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
284  }
285  }
286 
287  if (fCaloClustersCont) {
288  fCaloClustersCont->ResetCurrentID();
289  AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster();
290  while(cluster) {
291  TLorentzVector nPart;
292  cluster->GetMomentum(nPart, fVertex);
293  fHistClustersPt[fCentBin]->Fill(nPart.Pt());
294 
296  }
297  }
298 
299  int N1 = 0;
300  if (fJetsCont) {
301  fJetsCont->ResetCurrentID();
303  while(jet) {
304  Float_t NEFpT = jet->Pt()*jet->NEF();
305  N1++;
306  fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());
307  fHistJetsPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi());
308  Float_t ptLeading = fJetsCont->GetLeadingHadronPt(jet);
309  fHistJetsPtLeadHad[fCentBin]->Fill(jet->Pt(), ptLeading);
310  Double_t jetarray[] = {jet->Pt(),jet->Eta(),jet->Phi(),jet->NEF(),NEFpT};
311  fHistJet1->Fill(jetarray);
313  Float_t corrPt = jet->Pt() - fJetsCont->GetRhoVal() * jet->Area();
314  fHistJetsCorrPtArea[fCentBin]->Fill(corrPt, jet->Area());
315  }
316  jet = fJetsCont->GetNextAcceptJet();
317  }
318 
319  jet = fJetsCont->GetLeadingJet();
320  if(jet) fHistLeadingJetPt[fCentBin]->Fill(jet->Pt());
321  }//Loop over the first collection of jets
322 
323  int N2 = 0;
324  if(fJetsCont2){
325  fJetsCont2->ResetCurrentID();
327  while(jet){
328  Float_t NEFpT = jet->Pt()*jet->NEF();
329  N2++;
330  Double_t jetarray[] = {jet->Pt(),jet->Eta(),jet->Phi(),jet->NEF(),NEFpT};
331  fHistJet2->Fill(jetarray);
332  jet = fJetsCont2->GetNextAcceptJet();
333  }
334  } // loop over the trigger jerts.
335  int N1N2 = 0;
336  int N1N2m = 0;
337  if (fJetsCont&&fJetsCont2) {
338  fJetsCont->ResetCurrentID();
339  fJetsCont2->ResetCurrentID();
342  while(jet1){
343  bool ismatched = false;
344  Float_t NEFpT1 = jet1->Pt()*jet1->NEF();
345  Double_t jetarray1[] = {jet1->Pt(),jet1->Eta(),jet1->Phi(),jet1->NEF(),NEFpT1};
346  while(jet2){
347  N1N2++;
348  Double_t deta = jet1->Eta()-jet2->Eta();
349  Double_t dphi = RelativePhi(jet1->Phi(),jet2->Phi());
350  Double_t deta2 = deta*deta;
351  Double_t dphi2 = dphi*dphi;
352  Double_t dR = pow(deta2+dphi2,0.5);
353  Double_t jetarray[] = {jet1->Pt(),jet1->Eta(),jet1->Phi(),jet1->NEF(),jet2->Pt(),jet2->Eta(),jet2->Phi(),jet2->NEF(),dR};
354  if (dR<0.2){
355  N1N2m++;
356  fHistJet1to2->Fill(jetarray);
357  ismatched = true;
358  }
359  jet2 = fJetsCont2->GetNextAcceptJet();
360  }
361  if (ismatched)
362  fHistJet1m->Fill(jetarray1);
363  else
364  fHistJet1nm->Fill(jetarray1);
365  jet1 = fJetsCont->GetNextAcceptJet();
366  fJetsCont2->ResetCurrentID();
367  jet2 = fJetsCont2->GetNextAcceptJet();
368  }
369  }
370 
371 
372  if (fJetsCont&&fJetsCont3) {
373  fJetsCont->ResetCurrentID();
375  fJetsCont3->ResetCurrentID();
377  while(jet1){
378  while(jet3){
379  Double_t deta = jet1->Eta()-jet3->Eta();
380  Double_t dphi = RelativePhi(jet1->Phi(),jet3->Phi());
381  Double_t deta2 = deta*deta;
382  Double_t dphi2 = dphi*dphi;
383  Double_t Aj = (jet1->Pt()-jet3->Pt())/(jet1->Pt()+jet3->Pt());
384  Double_t jetarray[] = {jet1->Pt(),jet1->Eta(),jet1->Phi(),jet1->NEF(),jet3->Pt(),jet3->Eta(),jet3->Phi(),jet3->NEF(),Aj};
385  //Marta used |dphi - pi|<pi/3
386  if (fabs(fabs(dphi)-TMath::Pi())< TMath::Pi()/3.0){//dijet
387  fHistDiJet1->Fill(jetarray);
388  //we have a dijet, lets see if there is also a matched trigger
389  if (fJetsCont2) {
390  fJetsCont2->ResetCurrentID();
392  while(jet2){
393  Double_t tdeta = jet1->Eta()-jet2->Eta();
394  Double_t tdphi = RelativePhi(jet1->Phi(),jet2->Phi());
395  Double_t tdeta2 = tdeta*tdeta;
396  Double_t tdphi2 = tdphi*tdphi;
397  Double_t dR = pow(tdeta2+tdphi2,0.5);
398 
399  if (dR<0.2){
400  Double_t jetarray3[] = {jet1->Pt(),jet1->Eta(),jet1->Phi(),jet1->NEF(),jet3->Pt(),jet3->Eta(),jet3->Phi(),jet3->NEF(),Aj,jet2->Pt(),jet2->Eta(),jet2->Phi(),jet2->NEF(),dR};
401  //this dijet is triggered on!
402  fHistDiJet1m->Fill(jetarray3);
403  }
404  jet2 = fJetsCont2->GetNextAcceptJet();
405  } //while jet2
406  } // if jetscont2
407  }// if dijet
408  jet3 = fJetsCont3->GetNextAcceptJet();
409  }//while jet 3
410  jet1 = fJetsCont->GetNextAcceptJet();
411  fJetsCont3->ResetCurrentID();
412  jet3 = fJetsCont3->GetNextAcceptJet();
413  }//while jet 1
414  } //if jet cont 1 and 3
415 
416  return kTRUE;
417 }
418 
419 //________________________________________________________________________
421 {
422  if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
423  else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
424  if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
425  else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
426  double dphi = mphi-vphi;
427  if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
428  else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
429 
430  return dphi;//dphi in [-Pi, Pi]
431 }
432 
433 
434 //________________________________________________________________________
436 
438 
439  if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
440  if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
441  if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
442 
443 }
444 
445 //________________________________________________________________________
447 {
448  // Run analysis code here, if needed. It will be executed before FillHistograms().
449 
450  return kTRUE; // If return kFALSE FillHistogram() will NOT be executed.
451 }
452 
453 //________________________________________________________________________
455 {
456  // Called once at the end of the analysis.
457 }
Double_t Area() const
Definition: AliEmcalJet.h:130
TObjArray fClusterCollArray
cluster collection array
virtual AliVParticle * GetNextAcceptParticle()
Double_t GetRhoVal() const
double Double_t
Definition: External.C:58
TH2 ** fHistTracksEtaPhi
Track pt spectrum.
Bool_t FillHistograms()
Function filling histograms.
THnSparse * fHistJet1to2
jet collection 2
Definition: External.C:236
Bool_t Run()
Run function. This is the core function of the analysis and contains the user code. Therefore users have to implement this function.
AliJetContainer * GetJetContainer(Int_t i=0) const
Double_t Eta() const
Definition: AliEmcalJet.h:121
const TString & GetRhoName(Int_t c=0) const
Double_t Phi() const
Definition: AliEmcalJet.h:117
Double_t fMinBinPt
min pt in histograms
void ExecOnce()
Perform steps needed to initialize the analysis.
AliClusterContainer * GetClusterContainer() const
Int_t fCentBin
!event centrality bin
THnSparse * fHistJet2
jet collection 1 unmatched
AliJetContainer * fJetsCont2
Jets Jet 1.
TObjArray fParticleCollArray
particle/track collection array
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
AliParticleContainer * GetParticleContainer() const
AliEmcalJet * GetLeadingJet(const char *opt="")
AliJetContainer * fJetsCont
Dijet collection 1 and 3 matched.
AliClusterContainer * fCaloClustersCont
Tracks.
int Int_t
Definition: External.C:63
Double_t GetLeadingHadronPt(const AliEmcalJet *jet) const
float Float_t
Definition: External.C:68
Int_t fNcentBins
how many centrality bins
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Get cluster container attached to this task.
THnSparse * fHistJet1m
jet collection 1
TH2 ** fHistJetsPtLeadHad
Jet pt vs. area.
TH2 ** fHistJetsPhiEta
Leading jet pt spectrum.
TH2 ** fHistJetsPtArea
Phi-Eta distribution of jets.
AliEmcalJet * GetNextAcceptJet()
THnSparse * fHistDiJet1
jet collection 1 and 2
TObjArray fJetCollArray
jet collection array
TH2 ** fHistJetsCorrPtArea
Jet pt vs. leading hadron.
Double_t Pt() const
Definition: AliEmcalJet.h:109
AliParticleContainer * fTracksCont
Jets DiJet.
AliEmcalList * fOutput
!output list
Double_t fMaxBinPt
max pt in histograms
Float_t RelativePhi(Double_t mphi, Double_t vphi) const
Definition: External.C:220
Double_t fVertex[3]
!event vertex
void SetMakeGeneralHistograms(Bool_t g)
THnSparse * fHistJet1nm
jet collection 1 matched
Base task in the EMCAL jet framework.
THnSparse * fHistDiJet1m
Dijet collection 1 and 3.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
THnSparse * fHistJet1
Jet pt - bkg vs. area.
const char Option_t
Definition: External.C:48
void UserCreateOutputObjects()
Main initialization function on the worker.
const Int_t nbins
bool Bool_t
Definition: External.C:53
Double_t NEF() const
Definition: AliEmcalJet.h:148
TH1 ** fHistLeadingJetPt
Cluster pt spectrum.
AliJetContainer * fJetsCont3
Jets Trigger Jer.
AliVCluster * GetNextAcceptCluster()
Definition: External.C:196
Int_t fNbins
no. of pt bins