AliPhysics  e59a9ba (e59a9ba)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskSoftDrop.cxx
Go to the documentation of this file.
1 #include <TClonesArray.h>
2 #include <TH1F.h>
3 #include <TH2F.h>
4 #include <TH3F.h>
5 #include <TList.h>
6 #include <TLorentzVector.h>
7 
8 #include "AliVCluster.h"
9 #include "AliAODCaloCluster.h"
10 #include "AliESDCaloCluster.h"
11 #include "AliVTrack.h"
12 #include "AliEmcalJet.h"
13 #include "AliRhoParameter.h"
14 #include "AliLog.h"
15 #include "AliJetContainer.h"
16 #include "AliParticleContainer.h"
17 #include "AliClusterContainer.h"
18 #include "AliPicoTrack.h"
19 
21 
23 
24 //________________________________________________________________________
27  fHistTracksPt(0),
28  fHistNTracks(0),
29  fHistClustersPt(0),
30  fHistLeadingJetPt(0),
31  fHistJetsPhiEta(0),
32  fHistJetsPtArea(0),
33  fHistJetsPtLeadHad(0),
34  fHistJetsCorrPtArea(0),
35  fHistPtDEtaDPhiTrackClus(0),
36  fHistPtDEtaDPhiClusTrack(0),
37  fHistClustDx(0),
38  fHistClustDz(0),
39  fNAccJets(0),
40  fhZg(0),
41  fJetsCont(0),
42  fTracksCont(0),
43  fCaloClustersCont(0)
44 
45 {
46  // Default constructor.
47 
48  fHistTracksPt = new TH1*[fNcentBins];
49  fHistNTracks = 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  fHistNTracks[i] = 0;
60  fHistClustersPt[i] = 0;
61  fHistLeadingJetPt[i] = 0;
62  fHistJetsPhiEta[i] = 0;
63  fHistJetsPtArea[i] = 0;
64  fHistJetsPtLeadHad[i] = 0;
65  fHistJetsCorrPtArea[i] = 0;
66  }
67 
68  SetMakeGeneralHistograms(kTRUE);
69 }
70 
71 //________________________________________________________________________
73  AliAnalysisTaskEmcalJet(name, kTRUE),
74  fHistTracksPt(0),
75  fHistNTracks(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  fNAccJets(0),
87  fhZg(0),
88  fJetsCont(0),
89  fTracksCont(0),
90  fCaloClustersCont(0)
91 {
92  // Standard constructor.
93 
94  fHistTracksPt = new TH1*[fNcentBins];
95  fHistNTracks = new TH1*[fNcentBins];
96  fHistClustersPt = new TH1*[fNcentBins];
97  fHistLeadingJetPt = new TH1*[fNcentBins];
98  fHistJetsPhiEta = new TH2*[fNcentBins];
99  fHistJetsPtArea = new TH2*[fNcentBins];
100  fHistJetsPtLeadHad = new TH2*[fNcentBins];
101  fHistJetsCorrPtArea = new TH2*[fNcentBins];
102 
103  for (Int_t i = 0; i < fNcentBins; i++) {
104  fHistTracksPt[i] = 0;
105  fHistNTracks[i] = 0;
106  fHistClustersPt[i] = 0;
107  fHistLeadingJetPt[i] = 0;
108  fHistJetsPhiEta[i] = 0;
109  fHistJetsPtArea[i] = 0;
110  fHistJetsPtLeadHad[i] = 0;
111  fHistJetsCorrPtArea[i] = 0;
112  }
113 
115 }
116 
117 //________________________________________________________________________
119 {
120  // Destructor.
121 }
122 
123 //________________________________________________________________________
125 {
126  // Create user output.
127 
129 
131  if(fJetsCont) { //get particles and clusters connected to jets
134  } else { //no jets, just analysis tracks and clusters
137  }
138  if(fTracksCont) fTracksCont->SetClassName("AliVTrack");
139  if(fCaloClustersCont) fCaloClustersCont->SetClassName("AliVCluster");
140 
141  TString histname;
142 
143  for (Int_t i = 0; i < fNcentBins; i++) {
144  if (fParticleCollArray.GetEntriesFast()>0) {
145  histname = "fHistTracksPt_";
146  histname += i;
147  fHistTracksPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2);
148  fHistTracksPt[i]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
149  fHistTracksPt[i]->GetYaxis()->SetTitle("counts");
150  fOutput->Add(fHistTracksPt[i]);
151 
152  histname = "fHistNTracks_";
153  histname += i;
154  fHistNTracks[i] = new TH1F(histname.Data(), histname.Data(), 200, 0., 199.);
155  fOutput->Add(fHistNTracks[i]);
156  }
157 
158  if (fClusterCollArray.GetEntriesFast()>0) {
159  histname = "fHistClustersPt_";
160  histname += i;
161  fHistClustersPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2);
162  fHistClustersPt[i]->GetXaxis()->SetTitle("p_{T,clus} (GeV/c)");
163  fHistClustersPt[i]->GetYaxis()->SetTitle("counts");
164  fOutput->Add(fHistClustersPt[i]);
165  }
166 
167  if (fJetCollArray.GetEntriesFast()>0) {
168  histname = "fHistLeadingJetPt_";
169  histname += i;
170  fHistLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
171  fHistLeadingJetPt[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
172  fHistLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
173  fOutput->Add(fHistLeadingJetPt[i]);
174 
175  histname = "fHistJetsPhiEta_";
176  histname += i;
177  fHistJetsPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 50, -1, 1, 101, 0, TMath::Pi()*2 + TMath::Pi()/200);
178  fHistJetsPhiEta[i]->GetXaxis()->SetTitle("#eta");
179  fHistJetsPhiEta[i]->GetYaxis()->SetTitle("#phi");
180  fOutput->Add(fHistJetsPhiEta[i]);
181 
182  histname = "fHistJetsPtArea_";
183  histname += i;
184  fHistJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, 50, 0, 1);
185  fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
186  fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
187  fOutput->Add(fHistJetsPtArea[i]);
188 
189  histname = "fHistJetsPtLeadHad_";
190  histname += i;
191  fHistJetsPtLeadHad[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
192  fHistJetsPtLeadHad[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
193  fHistJetsPtLeadHad[i]->GetYaxis()->SetTitle("p_{T,lead} (GeV/c)");
194  fHistJetsPtLeadHad[i]->GetZaxis()->SetTitle("counts");
195  fOutput->Add(fHistJetsPtLeadHad[i]);
196 
197  if (!(GetJetContainer()->GetRhoName().IsNull())) {
198  histname = "fHistJetsCorrPtArea_";
199  histname += i;
200  fHistJetsCorrPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins*2, -fMaxBinPt, fMaxBinPt, 50, 0, 1);
201  fHistJetsCorrPtArea[i]->GetXaxis()->SetTitle("p_{T}^{corr} [GeV/c]");
202  fHistJetsCorrPtArea[i]->GetYaxis()->SetTitle("area");
203  fOutput->Add(fHistJetsCorrPtArea[i]);
204  }
205  }
206  }
207 
208  histname = "fHistPtDEtaDPhiTrackClus";
209  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);
211 
212  histname = "fHistPtDEtaDPhiClusTrack";
213  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);
215 
216  fHistClustDx = new TH1F("fHistClustDx","fHistClustDx;Dx",1000,0.,1.);
217  fOutput->Add(fHistClustDx);
218 
219  fHistClustDz = new TH1F("fHistClustDz","fHistClustDz;Dz",1000,0.,1.);
220  fOutput->Add(fHistClustDz);
221 
222  fNAccJets = new TH1F("fNAccJets","fNAccJets;N/ev",11,-0.5, 9.5);
223  fOutput->Add(fNAccJets);
224 
225  fhZg = new TH1F("fhZg", "#it{Z}_{g}; #it{Z}_{g}; Entries", 200, 0., 0.5);
226  fOutput->Add(fhZg);
227 
228  fhCorrPtZg = new TH2F("fhCorrPtZg", "#it{Z}_{g}; p_{T}^{corr} [GeV/c]; #it{Z}_{g}", 16, 0, 160, 20, 0., 0.5);
229  fOutput->Add(fhCorrPtZg);
230 
231  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
232 }
233 
234 //________________________________________________________________________
236 {
237  // Fill histograms.
238 
239  if (fTracksCont) {
241  fTracksCont->ResetCurrentID();
242  AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
243  while(track) {
244  fHistTracksPt[fCentBin]->Fill(track->Pt());
245  track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
246 
247  }
248  }
249 
250  if (fCaloClustersCont) {
251  AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster();
252  while(cluster) {
253  TLorentzVector nPart;
254  cluster->GetMomentum(nPart, fVertex);
255  fHistClustersPt[fCentBin]->Fill(nPart.Pt());
256  Double_t dx = cluster->GetTrackDx();
257  Double_t dz = cluster->GetTrackDz();
258  fHistClustDx->Fill(dx);
259  fHistClustDz->Fill(dz);
261  }
262  }
263 
264  if (fJetsCont) {
265  Int_t count = 0;
266  for (auto jet : fJetsCont->accepted() ) {
267  count++;
268  fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());
269  fHistJetsPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi());
270 
271  Float_t ptLeading = fJetsCont->GetLeadingHadronPt(jet);
272  fHistJetsPtLeadHad[fCentBin]->Fill(jet->Pt(), ptLeading);
273 
275  Float_t corrPt = jet->Pt() - fJetsCont->GetRhoVal() * jet->Area();
276  fHistJetsCorrPtArea[fCentBin]->Fill(corrPt, jet->Area());
277  }
278  fhZg->Fill(jet->GetShapeProperties()->GetSoftDropZg());
279  fhCorrPtZg->Fill(jet->Pt() - fJetsCont->GetRhoVal() * jet->Area(), jet->GetShapeProperties()->GetSoftDropZg());
280  }
281  fNAccJets->Fill(count);
282  auto jet = (AliEmcalJet*)fJetsCont->GetLeadingJet();
283  if(jet) fHistLeadingJetPt[fCentBin]->Fill(jet->Pt());
284  }
285 
287 
288  return kTRUE;
289 }
290 
291 //________________________________________________________________________
293 {
294 
296  return;
297 
298  Double_t deta = 999;
299  Double_t dphi = 999;
300 
301  //Get closest cluster to track
302  AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
303  while(track) {
304  //Get matched cluster
305  Int_t emc1 = track->GetEMCALcluster();
306  if(fCaloClustersCont && emc1>=0) {
307  AliVCluster *clusMatch = fCaloClustersCont->GetCluster(emc1);
308  if(clusMatch) {
309  AliPicoTrack::GetEtaPhiDiff(track, clusMatch, dphi, deta);
310  fHistPtDEtaDPhiTrackClus->Fill(track->Pt(),deta,dphi);
311  }
312  }
313  track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
314  }
315 
316  //Get closest track to cluster
317  AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster();
318  while(cluster) {
319  TLorentzVector nPart;
320  cluster->GetMomentum(nPart, fVertex);
321  fHistClustersPt[fCentBin]->Fill(nPart.Pt());
322 
323  //Get matched track
324  AliVTrack *mt = NULL;
325  AliAODCaloCluster *acl = dynamic_cast<AliAODCaloCluster*>(cluster);
326  if(acl) {
327  if(acl->GetNTracksMatched()>1)
328  mt = static_cast<AliVTrack*>(acl->GetTrackMatched(0));
329  }
330  else {
331  AliESDCaloCluster *ecl = dynamic_cast<AliESDCaloCluster*>(cluster);
332  Int_t im = ecl->GetTrackMatchedIndex();
333  if(fTracksCont && im>=0) {
334  mt = static_cast<AliVTrack*>(fTracksCont->GetParticle(im));
335  }
336  }
337  if(mt) {
338  AliPicoTrack::GetEtaPhiDiff(mt, cluster, dphi, deta);
339  fHistPtDEtaDPhiClusTrack->Fill(nPart.Pt(),deta,dphi);
340 
341  /* //debugging
342  if(mt->IsEMCAL()) {
343  Int_t emc1 = mt->GetEMCALcluster();
344  Printf("current id: %d emc1: %d",fCaloClustersCont->GetCurrentID(),emc1);
345  AliVCluster *clm = fCaloClustersCont->GetCluster(emc1);
346  AliPicoTrack::GetEtaPhiDiff(mt, clm, dphi, deta);
347  Printf("deta: %f dphi: %f",deta,dphi);
348  }
349  */
350  }
352  }
353 }
354 
355 //________________________________________________________________________
357 
359 
360  if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
361  if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
362  if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
363 
364 }
365 
366 //________________________________________________________________________
368 {
369  // Run analysis code here, if needed. It will be executed before FillHistograms().
370 
371  return kTRUE; // If return kFALSE FillHistogram() will NOT be executed.
372 }
373 
374 //________________________________________________________________________
376 {
377  // Called once at the end of the analysis.
378 }
TObjArray fClusterCollArray
cluster collection array
AliClusterContainer * fCaloClustersCont
Tracks.
virtual AliVParticle * GetNextAcceptParticle()
Double_t GetRhoVal() const
AliJetContainer * GetJetContainer(Int_t i=0) const
const TString & GetRhoName(Int_t c=0) const
TH3 * fHistPtDEtaDPhiClusTrack
track pt, delta eta, delta phi to matched cluster
Double_t fMinBinPt
min pt in histograms
AliClusterContainer * GetClusterContainer() const
Int_t fCentBin
!event centrality bin
TH1 ** fHistLeadingJetPt
Cluster pt spectrum.
TH1 ** fHistNTracks
cluster pt, delta eta, delta phi to matched track
TObjArray fParticleCollArray
particle/track collection array
AliParticleContainer * GetParticleContainer(Int_t i=0) const
AliParticleContainer * GetParticleContainer() const
TH2 ** fHistJetsPtArea
Phi-Eta distribution of jets.
AliEmcalJet * GetLeadingJet(const char *opt="")
Double_t GetLeadingHadronPt(const AliEmcalJet *jet) const
void Terminate(Option_t *option)
ClassImp(AliAnalysisTaskSoftDrop) AliAnalysisTaskSoftDrop
virtual AliVParticle * GetParticle(Int_t i=-1) const
Int_t fNcentBins
how many centrality bins
TH2 * fhCorrPtZg
! distribution of zg, jet pt-diff
AliClusterContainer * GetClusterContainer(Int_t i=0) const
TH2 ** fHistJetsCorrPtArea
Jet pt vs. leading hadron.
TH1 * fhZg
number of jets per event
AliParticleContainer * fTracksCont
Jets.
AliVCluster * GetCluster(Int_t i) const
TObjArray fJetCollArray
jet collection array
AliEmcalList * fOutput
!output list
Double_t fMaxBinPt
max pt in histograms
TH1 * fHistClustDx
number of tracks per event
TH2 ** fHistJetsPhiEta
Leading jet pt spectrum.
Double_t fVertex[3]
!event vertex
void SetMakeGeneralHistograms(Bool_t g)
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
TH1 ** fHistClustersPt
Track pt spectrum.
static void GetEtaPhiDiff(const AliVTrack *t, const AliVCluster *v, Double_t &phidiff, Double_t &etadiff)
const AliJetIterableContainer accepted() const
TH3 * fHistPtDEtaDPhiTrackClus
Jet pt - bkg vs. area.
Int_t GetNAcceptedParticles() const
AliVCluster * GetNextAcceptCluster()
TH2 ** fHistJetsPtLeadHad
Jet pt vs. area.
Int_t fNbins
no. of pt bins