AliPhysics  fde8a9f (fde8a9f)
 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  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
238 }
239 
240 //________________________________________________________________________
242 {
243  // Fill histograms.
244 
245  if (fTracksCont) {
247  fTracksCont->ResetCurrentID();
248  AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
249  while(track) {
250  fHistTracksPt[fCentBin]->Fill(track->Pt());
251  track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
252 
253  }
254  }
255 
256  if (fCaloClustersCont) {
257  AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster();
258  while(cluster) {
259  TLorentzVector nPart;
260  cluster->GetMomentum(nPart, fVertex);
261  fHistClustersPt[fCentBin]->Fill(nPart.Pt());
262  Double_t dx = cluster->GetTrackDx();
263  Double_t dz = cluster->GetTrackDz();
264  fHistClustDx->Fill(dx);
265  fHistClustDz->Fill(dz);
267  }
268  }
269 
270  if (fJetsCont) {
271  Int_t count = 0;
272  for (auto jet : fJetsCont->accepted() ) {
273  count++;
274  fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());
275  fHistJetsPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi());
276 
277  Float_t ptLeading = fJetsCont->GetLeadingHadronPt(jet);
278  fHistJetsPtLeadHad[fCentBin]->Fill(jet->Pt(), ptLeading);
279 
281  Float_t corrPt = jet->Pt() - fJetsCont->GetRhoVal() * jet->Area();
282  fHistJetsCorrPtArea[fCentBin]->Fill(corrPt, jet->Area());
283  }
284 
285  Double_t jetpt_ungrmd = jet->Pt() / ( jet->GetShapeProperties()->GetSoftDropPtfrac() );
286 
287  fhZg->Fill(jet->GetShapeProperties()->GetSoftDropZg());
288  fhCorrPtZg->Fill(jetpt_ungrmd - fJetsCont->GetRhoVal() * jet->Area(), jet->GetShapeProperties()->GetSoftDropZg() );
289  fhCorrPtRg->Fill(jetpt_ungrmd - fJetsCont->GetRhoVal() * jet->Area(), jet->GetShapeProperties()->GetSoftDropdR() );
290  fhCorrPtPtfrac->Fill(jetpt_ungrmd - fJetsCont->GetRhoVal() * jet->Area(), jet->GetShapeProperties()->GetSoftDropPtfrac() );
291  }
292  fNAccJets->Fill(count);
293  auto jet = (AliEmcalJet*)fJetsCont->GetLeadingJet();
294  if(jet) fHistLeadingJetPt[fCentBin]->Fill(jet->Pt());
295  }
296 
298 
299  return kTRUE;
300 }
301 
302 //________________________________________________________________________
304 {
305 
307  return;
308 
309  Double_t deta = 999;
310  Double_t dphi = 999;
311 
312  //Get closest cluster to track
313  AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
314  while(track) {
315  //Get matched cluster
316  Int_t emc1 = track->GetEMCALcluster();
317  if(fCaloClustersCont && emc1>=0) {
318  AliVCluster *clusMatch = fCaloClustersCont->GetCluster(emc1);
319  if(clusMatch) {
320  AliPicoTrack::GetEtaPhiDiff(track, clusMatch, dphi, deta);
321  fHistPtDEtaDPhiTrackClus->Fill(track->Pt(),deta,dphi);
322  }
323  }
324  track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
325  }
326 
327  //Get closest track to cluster
328  AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster();
329  while(cluster) {
330  TLorentzVector nPart;
331  cluster->GetMomentum(nPart, fVertex);
332  fHistClustersPt[fCentBin]->Fill(nPart.Pt());
333 
334  //Get matched track
335  AliVTrack *mt = NULL;
336  AliAODCaloCluster *acl = dynamic_cast<AliAODCaloCluster*>(cluster);
337  if(acl) {
338  if(acl->GetNTracksMatched()>1)
339  mt = static_cast<AliVTrack*>(acl->GetTrackMatched(0));
340  }
341  else {
342  AliESDCaloCluster *ecl = dynamic_cast<AliESDCaloCluster*>(cluster);
343  Int_t im = ecl->GetTrackMatchedIndex();
344  if(fTracksCont && im>=0) {
345  mt = static_cast<AliVTrack*>(fTracksCont->GetParticle(im));
346  }
347  }
348  if(mt) {
349  AliPicoTrack::GetEtaPhiDiff(mt, cluster, dphi, deta);
350  fHistPtDEtaDPhiClusTrack->Fill(nPart.Pt(),deta,dphi);
351 
352  /* //debugging
353  if(mt->IsEMCAL()) {
354  Int_t emc1 = mt->GetEMCALcluster();
355  Printf("current id: %d emc1: %d",fCaloClustersCont->GetCurrentID(),emc1);
356  AliVCluster *clm = fCaloClustersCont->GetCluster(emc1);
357  AliPicoTrack::GetEtaPhiDiff(mt, clm, dphi, deta);
358  Printf("deta: %f dphi: %f",deta,dphi);
359  }
360  */
361  }
363  }
364 }
365 
366 //________________________________________________________________________
368 
370 
371  if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
372  if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
373  if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
374 
375 }
376 
377 //________________________________________________________________________
379 {
380  // Run analysis code here, if needed. It will be executed before FillHistograms().
381 
382  return kTRUE; // If return kFALSE FillHistogram() will NOT be executed.
383 }
384 
385 //________________________________________________________________________
387 {
388  // Called once at the end of the analysis.
389 }
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
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
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
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.
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
AliVCluster * GetNextAcceptCluster()
TH2 ** fHistJetsPtLeadHad
Jet pt vs. area.
Definition: External.C:196
Int_t fNbins
no. of pt bins