AliPhysics  b752f14 (b752f14)
 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 
95  fHistNTracks = new TH1*[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  fhCorrPtRg = new TH2F("fhCorrPtRg", "R_{g}; p_{T}^{corr} [GeV/c]; R_{g}", 16, 0, 160, 40, 0., 0.5);
232  fOutput->Add(fhCorrPtRg);
233 
234  fhCorrPtPtfrac = new TH2F("fhCorrPtPtfrac", "#deltap_{T}; p_{T}^{corr} [GeV/c]; #deltap_{T}", 16, 0, 160, 80, 0., 1.0);
235  fOutput->Add(fhCorrPtPtfrac);
236 
237  fhCorrPtDropCount = new TH2F("fhCorrPtDropCount", "fhCorrPtDropCount; p_{T}^{corr} [GeV/c]; Counts", 16, 0, 160, 50, 0., 50);
239 
240  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
241 }
242 
243 //________________________________________________________________________
245 {
246  // Fill histograms.
247 
248  if (fTracksCont) {
250  fTracksCont->ResetCurrentID();
251  AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
252  while(track) {
253  fHistTracksPt[fCentBin]->Fill(track->Pt());
254  track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
255 
256  }
257  }
258 
259  if (fCaloClustersCont) {
260  AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster();
261  while(cluster) {
262  TLorentzVector nPart;
263  cluster->GetMomentum(nPart, fVertex);
264  fHistClustersPt[fCentBin]->Fill(nPart.Pt());
265  Double_t dx = cluster->GetTrackDx();
266  Double_t dz = cluster->GetTrackDz();
267  fHistClustDx->Fill(dx);
268  fHistClustDz->Fill(dz);
270  }
271  }
272 
273  if (fJetsCont) {
274  Int_t count = 0;
275  for (auto jet : fJetsCont->accepted() ) {
276  count++;
277  fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());
278  fHistJetsPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi());
279 
280  Float_t ptLeading = fJetsCont->GetLeadingHadronPt(jet);
281  fHistJetsPtLeadHad[fCentBin]->Fill(jet->Pt(), ptLeading);
282 
284  Float_t corrPt = jet->Pt() - fJetsCont->GetRhoVal() * jet->Area();
285  fHistJetsCorrPtArea[fCentBin]->Fill(corrPt, jet->Area());
286  }
287 
288  Double_t jetpt_ungrmd = jet->Pt() / ( jet->GetShapeProperties()->GetSoftDropPtfrac() );
289 
290  fhZg->Fill(jet->GetShapeProperties()->GetSoftDropZg());
291  fhCorrPtZg->Fill(jetpt_ungrmd - fJetsCont->GetRhoVal() * jet->Area(), jet->GetShapeProperties()->GetSoftDropZg() );
292  fhCorrPtRg->Fill(jetpt_ungrmd - fJetsCont->GetRhoVal() * jet->Area(), jet->GetShapeProperties()->GetSoftDropdR() );
293  fhCorrPtPtfrac->Fill(jetpt_ungrmd - fJetsCont->GetRhoVal() * jet->Area(), jet->GetShapeProperties()->GetSoftDropPtfrac() );
294  fhCorrPtDropCount->Fill(jetpt_ungrmd - fJetsCont->GetRhoVal() * jet->Area(), jet->GetShapeProperties()->GetSoftDropDropCount() );
295  }
296  fNAccJets->Fill(count);
297  auto jet = (AliEmcalJet*)fJetsCont->GetLeadingJet();
298  if(jet) fHistLeadingJetPt[fCentBin]->Fill(jet->Pt());
299  }
300 
302 
303  return kTRUE;
304 }
305 
306 //________________________________________________________________________
308 {
309 
311  return;
312 
313  Double_t deta = 999;
314  Double_t dphi = 999;
315 
316  //Get closest cluster to track
317  AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
318  while(track) {
319  //Get matched cluster
320  Int_t emc1 = track->GetEMCALcluster();
321  if(fCaloClustersCont && emc1>=0) {
322  AliVCluster *clusMatch = fCaloClustersCont->GetCluster(emc1);
323  if(clusMatch) {
324  AliPicoTrack::GetEtaPhiDiff(track, clusMatch, dphi, deta);
325  fHistPtDEtaDPhiTrackClus->Fill(track->Pt(),deta,dphi);
326  }
327  }
328  track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
329  }
330 
331  //Get closest track to cluster
332  AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster();
333  while(cluster) {
334  TLorentzVector nPart;
335  cluster->GetMomentum(nPart, fVertex);
336  fHistClustersPt[fCentBin]->Fill(nPart.Pt());
337 
338  //Get matched track
339  AliVTrack *mt = NULL;
340  AliAODCaloCluster *acl = dynamic_cast<AliAODCaloCluster*>(cluster);
341  if(acl) {
342  if(acl->GetNTracksMatched()>1)
343  mt = static_cast<AliVTrack*>(acl->GetTrackMatched(0));
344  }
345  else {
346  AliESDCaloCluster *ecl = dynamic_cast<AliESDCaloCluster*>(cluster);
347  Int_t im = ecl->GetTrackMatchedIndex();
348  if(fTracksCont && im>=0) {
349  mt = static_cast<AliVTrack*>(fTracksCont->GetParticle(im));
350  }
351  }
352  if(mt) {
353  AliPicoTrack::GetEtaPhiDiff(mt, cluster, dphi, deta);
354  fHistPtDEtaDPhiClusTrack->Fill(nPart.Pt(),deta,dphi);
355 
356  /* //debugging
357  if(mt->IsEMCAL()) {
358  Int_t emc1 = mt->GetEMCALcluster();
359  Printf("current id: %d emc1: %d",fCaloClustersCont->GetCurrentID(),emc1);
360  AliVCluster *clm = fCaloClustersCont->GetCluster(emc1);
361  AliPicoTrack::GetEtaPhiDiff(mt, clm, dphi, deta);
362  Printf("deta: %f dphi: %f",deta,dphi);
363  }
364  */
365  }
367  }
368 }
369 
370 //________________________________________________________________________
372 
374 
375  if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
376  if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
377  if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
378 
379 }
380 
381 //________________________________________________________________________
383 {
384  // Run analysis code here, if needed. It will be executed before FillHistograms().
385 
386  return kTRUE; // If return kFALSE FillHistogram() will NOT be executed.
387 }
388 
389 //________________________________________________________________________
391 {
392  // Called once at the end of the analysis.
393 }
TObjArray fClusterCollArray
cluster collection array
AliClusterContainer * fCaloClustersCont
Tracks.
virtual AliVParticle * GetNextAcceptParticle()
Double_t GetRhoVal() const
double Double_t
Definition: External.C:58
Definition: External.C:260
Definition: External.C:236
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
Get particle container attached to this task.
AliParticleContainer * GetParticleContainer() const
TH2 ** fHistJetsPtArea
Phi-Eta distribution of jets.
AliEmcalJet * GetLeadingJet(const char *opt="")
TH2 * fhCorrPtRg
! distribution of groomed jet radius, jet pt-diff
int Int_t
Definition: External.C:63
Double_t GetLeadingHadronPt(const AliEmcalJet *jet) const
float Float_t
Definition: External.C:68
void Terminate(Option_t *option)
ClassImp(AliAnalysisTaskSoftDrop) AliAnalysisTaskSoftDrop
virtual AliVParticle * GetParticle(Int_t i=-1) const
TH2 * fhCorrPtDropCount
! distribution of dropped branches number, jet pt-diff
Int_t fNcentBins
how many centrality bins
TH2 * fhCorrPtZg
! distribution of zg, jet pt-diff
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Get cluster container attached to this task.
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
Bool_t FillHistograms()
Function filling histograms.
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.
AliEmcalList * fOutput
!output list
Double_t fMaxBinPt
max pt in histograms
TH1 * fHistClustDx
number of tracks per event
Definition: External.C:220
TH2 ** fHistJetsPhiEta
Leading jet pt spectrum.
Double_t fVertex[3]
!event vertex
TH2 * fhCorrPtPtfrac
! distribution of ptfrac, jet pt-diff
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
const char Option_t
Definition: External.C:48
TH1 ** fHistClustersPt
Track pt spectrum.
void UserCreateOutputObjects()
Main initialization function on the worker.
static void GetEtaPhiDiff(const AliVTrack *t, const AliVCluster *v, Double_t &phidiff, Double_t &etadiff)
const AliJetIterableContainer accepted() const
bool Bool_t
Definition: External.C:53
TH3 * fHistPtDEtaDPhiTrackClus
Jet pt - bkg vs. area.
Int_t GetNAcceptedParticles() const
void ExecOnce()
Perform steps needed to initialize the analysis.
AliVCluster * GetNextAcceptCluster()
TH2 ** fHistJetsPtLeadHad
Jet pt vs. area.
Definition: External.C:196
Int_t fNbins
no. of pt bins