AliPhysics  41af4b0 (41af4b0)
AliAnalysisTaskJetCoreEmcal.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2016, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 #include <TClonesArray.h>
17 #include <TList.h>
18 #include "TChain.h"
19 #include "TTree.h"
20 #include "TMath.h"
21 #include "TH1F.h"
22 #include "TH2F.h"
23 #include "TH3F.h"
24 #include "THnSparse.h"
25 #include "TCanvas.h"
26 #include "TRandom3.h"
27 #include "AliLog.h"
28 
29 #include <AliAnalysisManager.h>
30 #include <AliVEventHandler.h>
31 #include <AliVEvent.h>
32 #include <AliVCluster.h>
33 #include <AliVParticle.h>
34 #include <AliLog.h>
35 #include <AliAnalysisDataSlot.h>
36 #include <AliAnalysisDataContainer.h>
37 
38 #include "AliTLorentzVector.h"
39 #include "AliEmcalJet.h"
40 #include "AliRhoParameter.h"
41 #include "AliJetContainer.h"
42 #include "AliParticleContainer.h"
43 #include "AliClusterContainer.h"
44 #include "AliEmcalPythiaInfo.h"
45 
47 
51 
57  fHistManager(),
58  fJetShapeType(AliAnalysisTaskJetCoreEmcal::kData),
59  fCentMin(0.),
60  fCentMax(100.),
61  fTTLowRef(8.),
62  fTTUpRef(9.),
63  fTTLowSig(20.),
64  fTTUpSig(50.),
65  fNRPBins(50),
66  fFrac(0.8),
67  fJetEtaMin(-.5),
68  fJetEtaMax(.5),
69  fJetHadronDeltaPhi(0.6),
70  fMinFractionSharedPt(0.5),
71  fMinEmbJetPt(15.),
72  fJetContName(""),
73  fJetContTrueName(""),
74  fJetContPartName(""),
75  fFillTrackHistograms(kTRUE),
76  fFillJetHistograms(kTRUE),
77  fFillRecoilTHnSparse(kTRUE),
78  fFillInclusiveTree(kFALSE),
79  fFillRecoilTree(kFALSE),
80  fPtHardBin(0.),
81  fRejectionFactorInclusiveJets(1),
82  fRandom(0),
83  fHistEvtSelection(0x0),
84  fHJetSpec(0x0),
85  fh1TrigRef(0x0),
86  fh1TrigSig(0x0),
87  fh2Ntriggers(0x0),
88  fh2RPJetsC10(0x0),
89  fh2RPJetsC20(0x0),
90  fh2RPTC10(0x0),
91  fh2RPTC20(0x0),
92  fhDphiPtSig(0x0),
93  fhDphiPtRef(0x0),
94  fhPtDetPart(0x0),
95  fhPtHybrDet(0x0),
96  fhPtHybrPart(0x0),
97  fhPtHybrPartCor(0x0),
98  fhPhiHybrPartCor(0x0),
99  fhPtDet(0x0),
100  fhPtDetMatchedToPart(0x0),
101  fhResidual(0x0),
102  fhPtResidual(0x0),
103  fhPhiResidual(0x0),
104  fhPhiPhiResidual(0x0),
105  fhPtDetPartRecoil(0x0),
106  fhPtHybrDetRecoil(0x0),
107  fhPtHybrPartRecoil(0x0),
108  fhPtHybrPartCorRecoil(0x0),
109  fhPtDetRecoil(0x0),
110  fhPtDetMatchedToPartRecoil(0x0),
111  fhResidualRecoil(0x0),
112  fhPtResidualRecoil(0x0),
113  fhDphiResidualRecoil(0x0),
114  fhDphiphiResidualRecoil(0x0),
115  fhTTPtDetMatchedToPart(0x0),
116  fhTTPhiDetMatchedToPart(0x0),
117  fhDPhiHybrPartCorRecoil(0x0),
118  fTreeEmbInclusive(0x0),
119  fTreeEmbRecoil(0x0)
120 {
122 
123  DefineOutput(1, TList::Class());
124  //if(fJetShapeType == AliAnalysisTaskJetCoreEmcal::kDetEmbPart && fFillInclusiveTree)
125  DefineOutput(2, TTree::Class());
126 // if(fJetShapeType == AliAnalysisTaskJetCoreEmcal::kDetEmbPart && fFillRecoilTree)
127  DefineOutput(3, TTree::Class());
128 }
129 
136  AliAnalysisTaskEmcalJet(name, kTRUE),
137  fHistManager(name),
139  fCentMin(0.),
140  fCentMax(100.),
141  fTTLowRef(8),
142  fTTUpRef(9.),
143  fTTLowSig(20.),
144  fTTUpSig(50.),
145  fNRPBins(50),
146  fFrac(0.8),
147  fJetEtaMin(-.5),
148  fJetEtaMax(.5),
149  fJetHadronDeltaPhi(0.6),
151  fMinEmbJetPt(15.),
152  fJetContName(""),
153  fJetContTrueName(""),
154  fJetContPartName(""),
155  fFillTrackHistograms(kTRUE),
156  fFillJetHistograms(kTRUE),
157  fFillRecoilTHnSparse(kTRUE),
158  fFillInclusiveTree(kFALSE),
159  fFillRecoilTree(kFALSE),
160  fPtHardBin(0.),
162  fRandom(0),
163  fHistEvtSelection(0x0),
164  fHJetSpec(0x0),
165  fh1TrigRef(0x0),
166  fh1TrigSig(0x0),
167  fh2Ntriggers(0x0),
168  fh2RPJetsC10(0x0),
169  fh2RPJetsC20(0x0),
170  fh2RPTC10(0x0),
171  fh2RPTC20(0x0),
172  fhDphiPtSig(0x0),
173  fhDphiPtRef(0x0),
174  fhPtDetPart(0x0),
175  fhPtHybrDet(0x0),
176  fhPtHybrPart(0x0),
177  fhPtHybrPartCor(0x0),
178  fhPhiHybrPartCor(0x0),
179  fhPtDet(0x0),
181  fhResidual(0x0),
182  fhPtResidual(0x0),
183  fhPhiResidual(0x0),
184  fhPhiPhiResidual(0x0),
185  fhPtDetPartRecoil(0x0),
186  fhPtHybrDetRecoil(0x0),
187  fhPtHybrPartRecoil(0x0),
189  fhPtDetRecoil(0x0),
191  fhResidualRecoil(0x0),
192  fhPtResidualRecoil(0x0),
198  fTreeEmbInclusive(0x0),
199  fTreeEmbRecoil(0x0)
200 {
202 
203  DefineOutput(1, TList::Class());
204  DefineOutput(2, TTree::Class());
205  DefineOutput(3, TTree::Class());
206 }
207 
212 {
213 }
214 
220 {
222 
225  //AllocateClusterHistograms();
226  //AllocateCellHistograms();
228 
229  TIter next(fHistManager.GetListOfHistograms());
230  TObject* obj = 0;
231  while ((obj = next())) {
232  fOutput->Add(obj);
233  }
234 
235  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
236  if(fJetShapeType == AliAnalysisTaskJetCoreEmcal::kDetEmbPart && fFillInclusiveTree) PostData(2, fTreeEmbInclusive); // Post data for ALL output slots > 0 here.
237  if(fJetShapeType == AliAnalysisTaskJetCoreEmcal::kDetEmbPart && fFillRecoilTree) PostData(3, fTreeEmbRecoil); // Post data for ALL output slots > 0 here.
238 }
239 
240 /*
241  * This function allocates the histograms for basic EMCal cluster QA.
242  * A set of histograms (energy, eta, phi, number of cluster) is allocated
243  * per each cluster container and per each centrality bin.
244  */
246 {
247  TString histname;
248  TString histtitle;
249  TString groupname;
250  AliClusterContainer* clusCont = 0;
251  TIter next(&fClusterCollArray);
252  while ((clusCont = static_cast<AliClusterContainer*>(next()))) {
253  groupname = clusCont->GetName();
254  // Protect against creating the histograms twice
255  if (fHistManager.FindObject(groupname)) {
256  AliWarning(TString::Format("%s: Found groupname %s in hist manager. The cluster containers will be filled into the same histograms.", GetName(), groupname.Data()));
257  continue;
258  }
259  fHistManager.CreateHistoGroup(groupname);
260  for (Int_t cent = 0; cent < fNcentBins; cent++) {
261  histname = TString::Format("%s/histClusterEnergy_%d", groupname.Data(), cent);
262  histtitle = TString::Format("%s;#it{E}_{cluster} (GeV);counts", histname.Data());
263  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
264 
265  histname = TString::Format("%s/histClusterEnergyExotic_%d", groupname.Data(), cent);
266  histtitle = TString::Format("%s;#it{E}_{cluster}^{exotic} (GeV);counts", histname.Data());
267  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
268 
269  histname = TString::Format("%s/histClusterNonLinCorrEnergy_%d", groupname.Data(), cent);
270  histtitle = TString::Format("%s;#it{E}_{cluster}^{non-lin.corr.} (GeV);counts", histname.Data());
271  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
272 
273  histname = TString::Format("%s/histClusterHadCorrEnergy_%d", groupname.Data(), cent);
274  histtitle = TString::Format("%s;#it{E}_{cluster}^{had.corr.} (GeV);counts", histname.Data());
275  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
276 
277  histname = TString::Format("%s/histClusterPhi_%d", groupname.Data(), cent);
278  histtitle = TString::Format("%s;#it{#phi}_{custer};counts", histname.Data());
279  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, 0, TMath::TwoPi());
280 
281  histname = TString::Format("%s/histClusterEta_%d", groupname.Data(), cent);
282  histtitle = TString::Format("%s;#it{#eta}_{custer};counts", histname.Data());
283  fHistManager.CreateTH1(histname, histtitle, fNbins / 6, -1, 1);
284 
285  histname = TString::Format("%s/histNClusters_%d", groupname.Data(), cent);
286  histtitle = TString::Format("%s;number of clusters;events", histname.Data());
287  if (fForceBeamType != kpp) {
288  fHistManager.CreateTH1(histname, histtitle, 500, 0, 3000);
289  }
290  else {
291  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
292  }
293  }
294  }
295 
296  histname = "fHistSumNClusters";
297  histtitle = TString::Format("%s;Sum of n clusters;events", histname.Data());
298  if (fForceBeamType != kpp) {
299  fHistManager.CreateTH1(histname, histtitle, 500, 0, 3000);
300  }
301  else {
302  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
303  }
304 }
305 
306 /*
307  * This function allocates the histograms for basic EMCal QA.
308  * One 2D histogram with the cell energy spectra and the number of cells
309  * per event is allocated per each centrality bin.
310  */
312 {
313  TString histname;
314  TString histtitle;
315  TString groupname(fCaloCellsName);
316 
317  fHistManager.CreateHistoGroup(groupname);
318  for (Int_t cent = 0; cent < fNcentBins; cent++) {
319  histname = TString::Format("%s/histCellEnergy_%d", groupname.Data(), cent);
320  histtitle = TString::Format("%s;#it{E}_{cell} (GeV);counts", histname.Data());
321  fHistManager.CreateTH1(histname, histtitle, 300, 0, 150);
322 
323  histname = TString::Format("%s/histNCells_%d", groupname.Data(), cent);
324  histtitle = TString::Format("%s;number of cells;events", histname.Data());
325  if (fForceBeamType != kpp) {
326  fHistManager.CreateTH1(histname, histtitle, 500, 0, 6000);
327  }
328  else {
329  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
330  }
331  }
332 }
333 
334 /*
335  * This function allocates the histograms for basic tracking QA.
336  * A set of histograms (pT, eta, phi, difference between kinematic properties
337  * at the vertex and at the EMCal surface, number of tracks) is allocated
338  * per each particle container and per each centrality bin.
339  */
341 {
342  TString histname;
343  TString histtitle;
344  TString groupname;
345  AliParticleContainer* partCont = 0;
346  TIter next(&fParticleCollArray);
347  while ((partCont = static_cast<AliParticleContainer*>(next()))) {
348  groupname = partCont->GetName();
349 
350  // Protect against creating the histograms twice
351  if (fHistManager.FindObject(groupname)) {
352  AliWarning(TString::Format("%s: Found groupname %s in hist manager. The track containers will be filled into the same histograms.", GetName(), groupname.Data()));
353  continue;
354  }
355  fHistManager.CreateHistoGroup(groupname);
356  for (Int_t cent = 0; cent < fNcentBins; cent++) {
357  histname = TString::Format("%s/histTrackPt_%d", groupname.Data(), cent);
358  histtitle = TString::Format("%s;#it{p}_{T,track} (GeV/#it{c});counts", histname.Data());
359  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
360 
361  histname = TString::Format("%s/histTrackPhi_%d", groupname.Data(), cent);
362  histtitle = TString::Format("%s;#it{#phi}_{track};counts", histname.Data());
363  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, 0, TMath::TwoPi());
364 
365  histname = TString::Format("%s/histTrackEta_%d", groupname.Data(), cent);
366  histtitle = TString::Format("%s;#it{#eta}_{track};counts", histname.Data());
367  fHistManager.CreateTH1(histname, histtitle, fNbins / 6, -1, 1);
368 
369  if (TClass(partCont->GetClassName()).InheritsFrom("AliVTrack")) {
370  histname = TString::Format("%s/fHistDeltaEtaPt_%d", groupname.Data(), cent);
371  histtitle = TString::Format("%s;#it{p}_{T,track}^{vertex} (GeV/#it{c});#it{#eta}_{track}^{vertex} - #it{#eta}_{track}^{EMCal};counts", histname.Data());
372  fHistManager.CreateTH2(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt, 50, -0.5, 0.5);
373 
374  histname = TString::Format("%s/fHistDeltaPhiPt_%d", groupname.Data(), cent);
375  histtitle = TString::Format("%s;#it{p}_{T,track}^{vertex} (GeV/#it{c});#it{#phi}_{track}^{vertex} - #it{#phi}_{track}^{EMCal};counts", histname.Data());
376  fHistManager.CreateTH2(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt, 200, -2, 2);
377 
378  histname = TString::Format("%s/fHistDeltaPtvsPt_%d", groupname.Data(), cent);
379  histtitle = TString::Format("%s;#it{p}_{T,track}^{vertex} (GeV/#it{c});#it{p}_{T,track}^{vertex} - #it{p}_{T,track}^{EMCal} (GeV/#it{c});counts", histname.Data());
380  fHistManager.CreateTH2(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt, fNbins / 2, -fMaxBinPt/2, fMaxBinPt/2);
381 
382  histname = TString::Format("%s/fHistEoverPvsP_%d", groupname.Data(), cent);
383  histtitle = TString::Format("%s;#it{P}_{track} (GeV/#it{c});#it{E}_{cluster} / #it{P}_{track} #it{c};counts", histname.Data());
384  fHistManager.CreateTH2(histname, histtitle, fNbins / 2, fMinBinPt, fMaxBinPt, fNbins / 2, 0, 4);
385  }
386 
387  histname = TString::Format("%s/histNTracks_%d", groupname.Data(), cent);
388  histtitle = TString::Format("%s;number of tracks;events", histname.Data());
389  if (fForceBeamType != kpp) {
390  fHistManager.CreateTH1(histname, histtitle, 500, 0, 5000);
391  }
392  else {
393  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
394  }
395  }
396  }
397 
398  histname = "fHistSumNTracks";
399  histtitle = TString::Format("%s;Sum of n tracks;events", histname.Data());
400  if (fForceBeamType != kpp) {
401  fHistManager.CreateTH1(histname, histtitle, 500, 0, 5000);
402  }
403  else {
404  fHistManager.CreateTH1(histname, histtitle, 200, 0, 200);
405  }
406 }
407 
408 /*
409  * This function allocates the histograms for basic jet QA.
410  * A set of histograms (pT, eta, phi, area, number of jets, corrected pT) is allocated
411  * per each jet container and per each centrality bin.
412  */
414 {
415  TString histname;
416  TString histtitle;
417  TString groupname;
418  AliJetContainer* jetCont = 0;
419  TIter next(&fJetCollArray);
420  while ((jetCont = static_cast<AliJetContainer*>(next()))) {
421  groupname = jetCont->GetName();
422  // Protect against creating the histograms twice
423  if (fHistManager.FindObject(groupname)) {
424  AliWarning(TString::Format("%s: Found groupname %s in hist manager. The jet containers will be filled into the same histograms.", GetName(), groupname.Data()));
425  continue;
426  }
427  fHistManager.CreateHistoGroup(groupname);
428  for (Int_t cent = 0; cent < fNcentBins; cent++) {
429  histname = TString::Format("%s/histJetPt_%d", groupname.Data(), cent);
430  histtitle = TString::Format("%s;#it{p}_{T,jet} (GeV/#it{c});counts", histname.Data());
431  fHistManager.CreateTH1(histname, histtitle, fNbins, fMinBinPt, fMaxBinPt);
432 
433  histname = TString::Format("%s/histJetArea_%d", groupname.Data(), cent);
434  histtitle = TString::Format("%s;#it{A}_{jet};counts", histname.Data());
435  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, 0, 3);
436 
437  histname = TString::Format("%s/histJetPhi_%d", groupname.Data(), cent);
438  histtitle = TString::Format("%s;#it{#phi}_{jet};counts", histname.Data());
439  fHistManager.CreateTH1(histname, histtitle, fNbins / 2, 0, TMath::TwoPi());
440 
441  histname = TString::Format("%s/histJetEta_%d", groupname.Data(), cent);
442  histtitle = TString::Format("%s;#it{#eta}_{jet};counts", histname.Data());
443  fHistManager.CreateTH1(histname, histtitle, fNbins / 6, -1, 1);
444 
445  histname = TString::Format("%s/histNJets_%d", groupname.Data(), cent);
446  histtitle = TString::Format("%s;number of jets;events", histname.Data());
447  if (fForceBeamType != kpp) {
448  fHistManager.CreateTH1(histname, histtitle, 500, 0, 500);
449  }
450  else {
451  fHistManager.CreateTH1(histname, histtitle, 100, 0, 100);
452  }
453 
454  if (!jetCont->GetRhoName().IsNull()) {
455  histname = TString::Format("%s/histJetCorrPt_%d", groupname.Data(), cent);
456  histtitle = TString::Format("%s;#it{p}_{T,jet}^{corr} (GeV/#it{c});counts", histname.Data());
457  fHistManager.CreateTH1(histname, histtitle, fNbins, -fMaxBinPt / 2, fMaxBinPt / 2);
458  }
459  }
460  }
461 }
462 
464 {
465 
466  Bool_t oldStatus = TH1::AddDirectoryStatus();
467  TH1::AddDirectory(kFALSE);
468 
469  // set seed
470  fRandom = new TRandom3(0);
471 
472  fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
473  fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
474  fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
475  fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
476  fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
477  fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
478  fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
479 
480  fh1TrigRef=new TH1D("Trig Ref","",10,0.,10);
481  fh1TrigSig=new TH1D("Trig Sig","",10,0.,10);
482  fh2Ntriggers=new TH2F("# of triggers","",100,0.,100.,50,0.,50.);
483  fh2RPJetsC10=new TH2F("RPJetC10","",35,0.,3.5,100,0.,100.);
484  fh2RPJetsC20=new TH2F("RPJetC20","",35,0.,3.5,100,0.,100.);
485  fh2RPTC10=new TH2F("RPTriggerC10","",35,0.,3.5,50,0.,50.);
486  fh2RPTC20=new TH2F("RPTriggerC20","",35,0.,3.5,50,0.,50.);
487 
489 
490  fOutput->Add(fh1TrigRef);
491  fOutput->Add(fh1TrigSig);
492  fOutput->Add(fh2Ntriggers);
493  fOutput->Add(fh2RPJetsC10);
494  fOutput->Add(fh2RPJetsC20);
495  fOutput->Add(fh2RPTC10);
496  fOutput->Add(fh2RPTC20);
497 
499  const Int_t dimSpec = 6;
500  const Int_t nBinsSpec[dimSpec] = {100,100, 280, 50,200, fNRPBins};
501  const Double_t lowBinSpec[dimSpec] = {0,0,-80, 0,-0.5*TMath::Pi(), 0};
502  const Double_t hiBinSpec[dimSpec] = {100,1, 200, 50,1.5*TMath::Pi(), static_cast<Double_t>(fNRPBins)};
503  fHJetSpec = new THnSparseF("fHJetSpec","Recoil jet spectrum",dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
504  }
505 
506  // comment out since I want finer binning in jet area, to make it easier
507  // to change selection on jet area (Leticia used 0.8*R^2*Pi whereas 0.6 is used
508  // for inclusive jets)
509  //change binning in jet area
510 // Double_t *xPt6 = new Double_t[7];
511 // xPt6[0] = 0.;
512 // xPt6[1]=0.07;
513 // xPt6[2]=0.2;
514 // xPt6[3]=0.4;
515 // xPt6[4]=0.6;
516 // xPt6[5]=0.8;
517 // xPt6[6]=1;
518 // fHJetSpec->SetBinEdges(1,xPt6);
519 // delete [] xPt6;
520 
521  fOutput->Add(fHJetSpec);
522 
523  // azimuthal correlation
524 
525  fhDphiPtSig = new TH2F("hDphiPtS","recoil #Delta #phi vs jet pT signal",100,-2,5,250,-50,200);
526  fhDphiPtSig->GetXaxis()->SetTitle("#Delta #phi");
527  fhDphiPtSig->GetYaxis()->SetTitle("p^{reco,ch}_{T,jet} (GeV/c)");
528  fhDphiPtRef = new TH2F("hDphiPtR","recoil #Delta #phi vs jet pT reference",100,-2,5,250,-50,200);
529  fhDphiPtRef->GetXaxis()->SetTitle("#Delta #phi");
530  fhDphiPtRef->GetYaxis()->SetTitle("p^{reco,ch}_{T,jet} (GeV/c)");
531 
532  fOutput->Add(fhDphiPtRef);
533  fOutput->Add(fhDphiPtSig);
534 
535 
536  fhPtHybrDet= new TH2F("hPtHybrDet","pT response Pb-Pb+PYTHIA vs PYTHIA",200,0,200,200,0,200);
537  fhPtHybrDet->GetXaxis()->SetTitle("p^{Pb-Pb+PYTHIA,ch}_{T} (GeV/c)");
538  fhPtHybrDet->GetYaxis()->SetTitle("p^{reco,PYTHIA,ch}_{T} (GeV/c)");
539  fhPtDetPart= new TH2F("hPtDetPart","pT response PYTHIA vs part",200,0,200,200,0,200);
540  fhPtDetPart->GetXaxis()->SetTitle("p^{reco,PYTHIA,ch}_{T} (GeV/c)");
541  fhPtDetPart->GetYaxis()->SetTitle("p^{true}_{T} (GeV/c)");
542  fhPtHybrPart = new TH2F("hPtHybrPart",Form("pT response Pb-Pb+PYTHIA vs part, min shared pT > %.0f",fMinFractionSharedPt*100),200,0,200,200,0,200);
543  fhPtHybrPart->GetXaxis()->SetTitle("p^{Pb-Pb+PYTHIA,ch}_{T} (GeV/c)");
544  fhPtHybrPart->GetYaxis()->SetTitle("p^{true}_{T} (GeV/c)");
545  fhPtHybrPartCor = new TH2F("hPtHybrPartCor",Form("pT response Pb-Pb+PYTHIA corrected vs part, min shared pT > %.0f",fMinFractionSharedPt*100),200,0,200,200,0,200);
546  fhPtHybrPartCor->GetXaxis()->SetTitle("p^{Pb-Pb+PYTHIA,ch}_{T} (GeV/c)");
547  fhPtHybrPartCor->GetYaxis()->SetTitle("p^{true}_{T} (GeV/c)");
548 
549  fhPhiHybrPartCor = new TH2F("hPhiHybrPartCor",Form("phi response Pb-Pb+PYTHIA corrected vs part, min shared pT > %.0f",fMinFractionSharedPt*100),200,-0.5*TMath::Pi(),1.5*TMath::Pi(),200,-0.5*TMath::Pi(),1.5*TMath::Pi());
550  fhPhiHybrPartCor->GetXaxis()->SetTitle("#phi^{Pb-Pb+PYTHIA,ch}");
551  fhPhiHybrPartCor->GetYaxis()->SetTitle("#phi^{true}");
552  fhResidual = new TH1F("hResidual","residual",50,-1,1);
553  fhResidual->GetXaxis()->SetTitle("p^{reco}_{T} - p^{true}_{T} / p^{true}_{T}");
554  fhPtResidual= new TH2F("hPtResidual","pT vs residual",200,0,200,50,-1,1);
555  fhPtResidual->GetXaxis()->SetTitle("p^{true}_{T} (GeV/c)");
556  fhPtResidual->GetYaxis()->SetTitle("p^{reco}_{T} - p^{true}_{T} / p^{true}_{T}");
557  fhPhiResidual = new TH1F("hPhiResidual","residual phi",50,-1,1);
558  fhPhiResidual->GetXaxis()->SetTitle("#phi^{reco} - #phi^{true} / #phi^{true}");
559  fhPhiPhiResidual= new TH2F("hPhiPhiResidual","pT vs residual",600,-3,3,200,-0.5*TMath::Pi(),1.5*TMath::Pi());
560  fhPhiPhiResidual->GetXaxis()->SetTitle("#phi^{true}");
561  fhPhiPhiResidual->GetYaxis()->SetTitle("#phi^{reco} - #phi^{true} / #phi^{true}");
562 
563  fhPtDet= new TH1F("hPtDet","pT detector level",200,0,200);
564  fhPtDet->GetXaxis()->SetTitle("p^{reco,PYTHIA,ch}_{T} (GeV/c)");
565  fhPtDetMatchedToPart = new TH1F("hPtDetMatchedToPart","pT detector level matched to particle level jet",200,0,200);
566  fhPtDetMatchedToPart->GetXaxis()->SetTitle("p^{reco,PYTHIA,ch}_{T} (GeV/c)");
567 
568  fhPtHybrDetRecoil= new TH2F("hPtHybrDetRecoil","pT response Pb-Pb+PYTHIA vs PYTHIA",200,0,200,200,0,200);
569  fhPtHybrDetRecoil->GetXaxis()->SetTitle("p^{Pb-Pb+PYTHIA,ch}_{T} (GeV/c)");
570  fhPtHybrDetRecoil->GetYaxis()->SetTitle("p^{reco,PYTHIA,ch}_{T} (GeV/c)");
571  fhPtDetPartRecoil= new TH2F("hPtDetPartRecoil","pT response PYTHIA vs part",200,0,200,200,0,200);
572  fhPtDetPartRecoil->GetXaxis()->SetTitle("p^{reco,PYTHIA,ch}_{T} (GeV/c)");
573  fhPtDetPartRecoil->GetYaxis()->SetTitle("p^{true}_{T} (GeV/c)");
574  fhPtHybrPartRecoil = new TH2F("hPtHybrPartRecoil",Form("pT response Pb-Pb+PYTHIA vs part, min shared pT > %.0f",fMinFractionSharedPt*100),200,0,200,200,0,200);
575  fhPtHybrPartRecoil->GetXaxis()->SetTitle("p^{Pb-Pb+PYTHIA,ch}_{T} (GeV/c)");
576  fhPtHybrPartRecoil->GetYaxis()->SetTitle("p^{true}_{T} (GeV/c)");
577  fhPtHybrPartCorRecoil = new TH2F("hPtHybrPartCorRecoil",Form("pT response Pb-Pb+PYTHIA corrected vs part, min shared pT > %.0f",fMinFractionSharedPt*100),200,0,200,200,0,200);
578  fhPtHybrPartCorRecoil->GetXaxis()->SetTitle("p^{Pb-Pb+PYTHIA,ch}_{T} (GeV/c)");
579  fhPtHybrPartCorRecoil->GetYaxis()->SetTitle("p^{true}_{T} (GeV/c)");
580 
581  fhDPhiHybrPartCorRecoil = new TH2F("hDPhiHybrPartCorRecoil",Form("#Delta#phi response Pb-Pb+PYTHIA corrected vs part, min shared pT > %.0f",fMinFractionSharedPt*100),200,-0.5*TMath::Pi(),1.5*TMath::Pi(),200,-0.5*TMath::Pi(),1.5*TMath::Pi());
582  fhDPhiHybrPartCorRecoil->GetXaxis()->SetTitle("#Delta#phi^{Pb-Pb+PYTHIA,ch}");
583  fhDPhiHybrPartCorRecoil->GetYaxis()->SetTitle("#Delta#phi^{true}");
584  fhResidualRecoil = new TH1F("hResidualRecoil","residual",50,-1,1);
585  fhResidualRecoil->GetXaxis()->SetTitle("p^{reco}_{T} - p^{true}_{T} / p^{true}_{T}");
586  fhPtResidualRecoil= new TH2F("hPtResidualRecoil","pT vs residual",200,0,200,50,-1,1);
587  fhPtResidualRecoil->GetXaxis()->SetTitle("p^{true}_{T} (GeV/c)");
588  fhPtResidualRecoil->GetYaxis()->SetTitle("p^{reco}_{T} - p^{true}_{T} / p^{true}_{T}");
589 
590  fhDphiResidualRecoil = new TH1F("hDphiResidualRecoil","residual",50,-1,1);
591  fhDphiResidualRecoil->GetXaxis()->SetTitle("#Delta#phi^{reco} - #Delta#phi^{true} / #Delta#phi^{true}");
592  fhDphiphiResidualRecoil= new TH2F("hDphiphiResidualRecoil","#Delta#phi vs residual in phi",400,-3,3,200,-0.5*TMath::Pi(),1.5*TMath::Pi());
593  fhDphiphiResidualRecoil->GetXaxis()->SetTitle("#Delta#phi^{true}");
594  fhDphiphiResidualRecoil->GetYaxis()->SetTitle("#Delta#phi^{reco} - #Delta#phi^{true} / #Delta#phi^{true}");
595  fhPtDetRecoil= new TH1F("hPtDetRecoil","pT detector level",200,0,200);
596  fhPtDetRecoil->GetXaxis()->SetTitle("p^{reco,PYTHIA,ch}_{T} (GeV/c)");
597  fhPtDetMatchedToPartRecoil = new TH1F("hPtDetMatchedToPartRecoil","pT detector level matched to particle level jet",200,0,200);
598  fhPtDetMatchedToPartRecoil->GetXaxis()->SetTitle("p^{reco,PYTHIA,ch}_{T} (GeV/c)");
599 
600  fhTTPtDetMatchedToPart = new TH2F("hTTPtDetMatchedToPart","trigger track pT response reco vs partice",140,0,70,140,0,70);
601  fhTTPtDetMatchedToPart->GetXaxis()->SetTitle("p^{TT,reco}_{T} (GeV/c)");
602  fhTTPtDetMatchedToPart->GetYaxis()->SetTitle("p^{TT,part}_{T} (GeV/c)");
603  fhTTPhiDetMatchedToPart = new TH2F("hTTPhiDetMatchedToPart","trigger track #phi response reco vs partice",200,-0.5*TMath::Pi(),1.5*TMath::Pi(),200,-0.5*TMath::Pi(),1.5*TMath::Pi());
604  fhTTPhiDetMatchedToPart->GetXaxis()->SetTitle("p^{TT,reco}_{T} (GeV/c)");
605  fhTTPhiDetMatchedToPart->GetYaxis()->SetTitle("p^{TT,part}_{T} (GeV/c)");
606 
607  fOutput->Add(fhPtHybrDet);
608  fOutput->Add(fhPtDetPart);
609  fOutput->Add(fhPtHybrPart);
610  fOutput->Add(fhPtHybrPartCor);
612  fOutput->Add(fhPtDet);
614  fOutput->Add(fhResidual);
615  fOutput->Add(fhPtResidual);
616  fOutput->Add(fhPhiResidual);
618 
624  fOutput->Add(fhPtDetRecoil);
630 
633 
634  TString varNamesInclusive[8]={"centrality","ptRawRec","areaRec","ptCorrRec","phiRec","ptPart","phiPart","binPtHard"};
635  TString varNamesRecoil[11]={"centrality","ptTT","ptRawRec","areaRec","ptCorrRec","phiRec","DPhiRec","ptPart","phiPart","DPhiPart","binPtHard"};
637  const char* nameEmbInclusive = GetOutputSlot(2)->GetContainer()->GetName();
638  fTreeEmbInclusive = new TTree(nameEmbInclusive, nameEmbInclusive);
639  for(Int_t ivar=0; ivar < 8; ivar++){
640  fTreeEmbInclusive->Branch(varNamesInclusive[ivar].Data(), &fTreeVarsInclusive[ivar], Form("%s/F", varNamesInclusive[ivar].Data()));
641  }
642  }
643 
645  const char* nameEmbRecoil= GetOutputSlot(3)->GetContainer()->GetName();
646  fTreeEmbRecoil = new TTree(nameEmbRecoil, nameEmbRecoil);
647  for(Int_t ivar=0; ivar < 11; ivar++){
648  fTreeEmbRecoil->Branch(varNamesRecoil[ivar].Data(), &fTreeVarsRecoil[ivar], Form("%s/F", varNamesRecoil[ivar].Data()));
649  }
650  }
651 
652  // =========== Switch on Sumw2 for all histos ===========
653  for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
654  TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
655  if (h1){
656  h1->Sumw2();
657  continue;
658  }
659  THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
660  if (hn){
661  hn->Sumw2();
662  }
663  }
664 
665  // add QA plots from fEventCuts
666  fEventCuts.AddQAplotsToList(fOutput);
667 
668  TH1::AddDirectory(oldStatus);
669 
670  PostData(1, fOutput);
673 }
674 
682 {
683 
684  fHistEvtSelection->Fill(1); // number of events before event selection
685  AliVEvent *ev = InputEvent();
686  if (!fEventCuts.AcceptEvent(ev)) {
687  fHistEvtSelection->Fill(2);
688  return kTRUE;
689  }
690 
693  //DoClusterLoop();
694  //DoCellLoop();
695  DoJetCoreLoop();
697 
698  return kTRUE;
699 }
700 
702 {
703  // Do jet core analysis and fill histograms.
704 
706  AliJetContainer *jetContPart = 0x0;
707  AliJetContainer *jetContTrue = 0x0;
709  jetContTrue = GetJetContainer(fJetContTrueName);
710  jetContPart = GetJetContainer(fJetContPartName);
711  }
712 
713  if(!jetCont ||
715  AliError(Form("jet container not found - check name %s",fJetContName.Data()));
716  TIter next(&fJetCollArray);
717  while ((jetCont = static_cast<AliJetContainer*>(next())))
718  AliError(Form("%s",jetCont->GetName()));
719  AliFatal("Exit...");
720  return;
721  }
722 
723  // centrality selection
724  if(fDebug) Printf("centrality: %f\n", fCent);
726  ((fCent>fCentMax) || (fCent<fCentMin))) {
727  fHistEvtSelection->Fill(4);
728  return;
729  }
730  fHistEvtSelection->Fill(0);
731 
732  // Background
733  Double_t rho = 0;
734  if (jetCont->GetRhoParameter()) rho = jetCont->GetRhoVal();
735  if(fDebug) Printf("rho = %f",rho);
736 
737  // get MC particle container in case running embedding, to match
738  // reconstructed and MC-level trigger tracks
739  AliParticleContainer *partCont = 0x0;
741 // AliParticleContainer *partCont = 0x0;
742 // TIter next(&fParticleCollArray);
743 // while ((partCont = static_cast<AliParticleContainer*>(next()))) {
744 // TString groupname = partCont->GetName();
745 // Printf("particle name = %s",groupname.Data());
746 // }
747 
748  // Choose trigger track
749  Int_t nT=0;
750  TList ParticleList;
751  Double_t minT=0;
752  Double_t maxT=0;
753  Int_t number=0;
754  Double_t dice=fRandom->Uniform(0,1);
755  Bool_t isSignal = kFALSE;
756  if(dice>fFrac){
757  minT=fTTLowRef;
758  maxT=fTTUpRef;
759  }
760  if(dice<=fFrac){
761  isSignal = kTRUE;
762  minT=fTTLowSig;
763  maxT=fTTUpSig;
764  }
765  nT=SelectTrigger(&ParticleList,minT,maxT,number);
766  if(fDebug) Printf("%s class ---> n triggers between %f and %f = %i, index of trigger chosen = %i",dice>fFrac?"ref.":"sig.",minT,maxT,number,nT);
767  if(nT<0) return;
768 
769  if(dice>fFrac) fh1TrigRef->Fill(number);
770  if(dice<=fFrac)fh1TrigSig->Fill(number);
771 
772 
773  // particle loop -
774  for(Int_t tt=0;tt<ParticleList.GetEntries();tt++){
775  // histogram
776  //if(fHardest==0||fHardest==1){if(tt!=nT) continue;}
777  if(tt!=nT) continue;
778  AliVParticle *partback = (AliVParticle*)ParticleList.At(tt);
779  if(!partback) continue;
780  if(fDebug) Printf("trigger particle pt = %f \teta = %f \t phi = %f",partback->Pt(),partback->Eta(),partback->Phi());
781  // if(partback->Pt()<8) continue;
782 
783  Int_t injet4=0;
784  Int_t injet=0;
785 
786  fh2Ntriggers->Fill(fCent,partback->Pt());
787  Double_t phiBinT = RelativePhi(partback->Phi(),fEPV0);
788  if(fCent<20.) fh2RPTC20->Fill(TMath::Abs(phiBinT),partback->Pt());
789  if(fCent<10.) fh2RPTC10->Fill(TMath::Abs(phiBinT),partback->Pt());
790 
791  Double_t etabig=0;
792  Double_t ptbig=0;
793  Double_t areabig=0;
794  Double_t phibig=0.;
795  // Double_t areasmall=0;
796 
797  TString histname;
798  TString groupname;
799  groupname = jetCont->GetName();
800  UInt_t count = 0;
801  for(auto jetbig : jetCont->accepted()) {
802  if (!jetbig) continue;
803  count++;
804  ptbig = jetbig->Pt();
805  etabig = jetbig->Eta();
806  phibig = jetbig->Phi();
807  if(ptbig==0) continue;
808  Double_t phiBin = RelativePhi(phibig,fEPV0); //relative phi between jet and ev. plane
809  areabig = jetbig->Area();
810  Double_t ptcorr=ptbig-rho*areabig;
811  //JJJ - perhaps should change eta selection if implemented in jet container
812  if((etabig<fJetEtaMin)||(etabig>fJetEtaMax)) continue;
813  if(areabig>=0.07) injet=injet+1;
814  if(areabig>=0.4) injet4=injet4+1;
815  Double_t dphi=RelativePhi(partback->Phi(),phibig);
816  if(fDebug) Printf("jet properties...\n\teta = %f \t phi = %f \t pt = %f \t relativephi = %f\t area = %f\t rho = %f",etabig,phibig,ptbig,dphi,areabig,rho);
817 
818  // do azimuthal correlation analysis
819  // dPhi between -0.5 < dPhi < 1.5
820  Double_t dPhiShift=phibig-partback->Phi();
821  if(dPhiShift>2*TMath::Pi()) dPhiShift -= 2*TMath::Pi();
822  if(dPhiShift<-2*TMath::Pi()) dPhiShift += 2*TMath::Pi();
823  if(dPhiShift<-0.5*TMath::Pi()) dPhiShift += 2*TMath::Pi();
824  if(dPhiShift>1.5*TMath::Pi()) dPhiShift -= 2*TMath::Pi();
825  if(isSignal) fhDphiPtSig->Fill(dPhiShift,ptcorr);
826  else fhDphiPtRef->Fill(dPhiShift,ptcorr);
827 
828  // selection on relative phi
829  if(fJetHadronDeltaPhi>0. &&
830  TMath::Abs(dphi)<TMath::Pi()-fJetHadronDeltaPhi) continue;
831 
832  if(fCent<10.) fh2RPJetsC10->Fill(TMath::Abs(phiBin), ptcorr);
833  if(fCent<20.) fh2RPJetsC20->Fill(TMath::Abs(phiBin), ptcorr);
834 
836  Float_t phitt=partback->Phi();
837  if(phitt<0)phitt+=TMath::Pi()*2.;
838  Int_t phiBintt = GetPhiBin(phitt-fEPV0);
839 
840  Double_t fillspec[] = {fCent,areabig,ptcorr,partback->Pt(),dPhiShift, static_cast<Double_t>(phiBintt)};
841  fHJetSpec->Fill(fillspec);
842  }
843 
845  //
846  // embedding for recoil jets
847  // get MC info
848  //
849  if(ptcorr<fMinEmbJetPt) continue;
850  Double_t ptTTMC = 0;
851  Double_t phiTTMC = 0;
852  Int_t TTmatched = 0;
853  for(auto partMC : partCont->accepted()) {
854  Int_t labtr = partback->GetLabel();
855  Int_t labpa = partMC->GetLabel();
856  if(labtr==labpa) {
857  ptTTMC = partMC->Pt();
858  phiTTMC = partMC->Phi();
859  TTmatched++;
860  break;
861  }
862  }
863  if(TTmatched!=1) continue;
864  Double_t ptTTreco = partback->Pt();
865  Double_t phiTTreco = partback->Phi();
866  if(fDebug) Printf("found corresponding truth-level particle, pt reco = %f pt MC = %f, phi reco = %f phi MC = %f",ptTTreco,ptTTMC,phiTTreco,phiTTMC);
867  fhTTPtDetMatchedToPart->Fill(ptTTreco,ptTTMC);
868  fhTTPhiDetMatchedToPart->Fill(phiTTreco,phiTTMC);
869 
870  auto jet2 = jetbig->ClosestJet();
871  if(!jet2) {
872  //Printf("jet 2 cant be found");
873  continue;}
874  Double_t ptJet2 = jet2->Pt();
875  fhPtDetRecoil->Fill(ptJet2);
876  auto jet3 = jet2->ClosestJet();
877  if(!jet3) {
878  //Printf("jet3 can't be found");
879  continue;
880  }
881  fhPtDetMatchedToPartRecoil->Fill(ptJet2);
882  Double_t ptJet3 = jet3->Pt();
883  Double_t phiJet3 = jet3->Phi();
884 
885 
886  Double_t dPhiPart=phiJet3-phiTTMC;
887  if(dPhiPart>2*TMath::Pi()) dPhiPart -= 2*TMath::Pi();
888  if(dPhiPart<-2*TMath::Pi()) dPhiPart += 2*TMath::Pi();
889  if(dPhiPart<-0.5*TMath::Pi()) dPhiPart += 2*TMath::Pi();
890  if(dPhiPart>1.5*TMath::Pi()) dPhiPart -= 2*TMath::Pi();
891 
892  if(fDebug) Printf("--- recoil - jet pt = jet hybrid pt = %f\t jet matched det pt = %f\t jet matched particle level pt = %f\t\n\tjet reco phi = %f\t jet particle phi = %f",ptbig,ptJet2,ptJet3,phibig,phiJet3);
893 
894  fhPtDetPartRecoil->Fill(ptJet2,ptJet3);
895  Double_t fraction = jetCont->GetFractionSharedPt(jetbig);
896  if(fraction < fMinFractionSharedPt) continue;
897 
898  Double_t residual = (ptcorr - ptJet3) / ptJet3;
899  Double_t residualDphi = (dPhiShift - dPhiPart) / dPhiPart;
900 
901  fhPtHybrDetRecoil->Fill(ptbig,ptJet2);
902  fhPtHybrPartRecoil->Fill(ptbig,ptJet3);
903  fhPtHybrPartCorRecoil->Fill(ptcorr,ptJet3);
904  fhDPhiHybrPartCorRecoil->Fill(dPhiShift,dPhiPart);
905  fhResidualRecoil->Fill(residual);
906  fhPtResidualRecoil->Fill(ptJet3,residual);
907  fhDphiResidualRecoil->Fill(residualDphi);
908  fhDphiphiResidualRecoil->Fill(dPhiPart,residualDphi);
909 
910  if(fFillRecoilTree) {
911  fTreeVarsRecoil[0] = fCent;
912  fTreeVarsRecoil[1] = partback->Pt();
913  fTreeVarsRecoil[2] = ptbig;
914  fTreeVarsRecoil[3] = areabig;
915  fTreeVarsRecoil[4] = ptcorr;
916  fTreeVarsRecoil[5] = phibig;
917  fTreeVarsRecoil[6] = dPhiShift;
918  fTreeVarsRecoil[7] = ptJet3;
919  fTreeVarsRecoil[8] = phiJet3;
920  fTreeVarsRecoil[9] = dPhiPart;
922  fTreeEmbRecoil->Fill();
923  }
924  }
925  }
926  }
927 }
928 
930 
934  if(fDebug) Printf("particle container 0 entries = %i \t1 entries = %i\t 2 entries = %i",partCont0->GetNParticles(),partCont1->GetNParticles(),partCont2->GetNParticles());
938 
939  if(!jetCont || !jetContPart || !jetContTrue)
940  {
941  AliError(Form("jet container not found - check name %s",fJetContName.Data()));
942  TIter next(&fJetCollArray);
943  while ((jetCont = static_cast<AliJetContainer*>(next())))
944  AliError(Form("%s",jetCont->GetName()));
945  AliFatal("Exit...");
946  return;
947  }
948 
949  if(fDebug) {
950  Printf("n particle jets = %i",jetContPart->GetNJets());
951  Printf("n reco jets = %i",jetCont->GetNJets());
952  Printf("n PYTHIA jets = %i",jetContTrue->GetNJets());
953  }
954 
955  // Background
956  Double_t rho = 0;
957  if (jetCont->GetRhoParameter()) rho = jetCont->GetRhoVal();
958 
959  // PYTHIA event weight
960  // note - not used
961  // AliGenPythiaEventHeader *pyHeader = 0x0; //!<! Pythia header of the current external event
962  // AliAODEvent *ev = dynamic_cast<AliAODEvent*>(InputEvent());
963  // AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(ev->FindListObject(AliAODMCHeader::StdBranchName()));
964  // if (aodMCH) {
965  // for (UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) {
966  // pyHeader= dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
967  // if (pyHeader) break;
968  // }
969  // }
970  //
971  // Double_t pythiaCrossSection = 0;
972  // Double_t pythiaTrials = 0;
973  // Double_t pythiaWeight = 0;
974  // if (pyHeader)
975  // {
976  // if(fDebug) Printf("have pythia header - get weight");
977  // pythiaCrossSection = pyHeader->GetXsection();
978  // pythiaTrials = pyHeader->Trials();
979  // //fPythiaPtHard = fPythiaHeader->GetPtHard();
980  // pythiaWeight = pythiaCrossSection / pythiaTrials;
981  // }
982  // if(fDebug) Printf("pythia weight = %f",pythiaWeight);
983 
984  for(auto jet1 : jetCont->accepted()) { // loop over hybrid jets
985 
986  Double_t ptJet1 = jet1->Pt();
987  Double_t phiJet1 = jet1->Phi();
988  Double_t area = jet1->Area();
989  Double_t ptCorr = ptJet1-rho*area;
990  if(ptCorr<fMinEmbJetPt) continue;
991  auto jet2 = jet1->ClosestJet();
992  if(!jet2) {
993  //Printf("jet 2 cant be found");
994  continue;}
995  Double_t ptJet2 = jet2->Pt();
996  fhPtDet->Fill(ptJet2);
997  auto jet3 = jet2->ClosestJet();
998  if(!jet3) {
999  //Printf("jet3 can't be found");
1000  continue;
1001  }
1002  fhPtDetMatchedToPart->Fill(ptJet2);
1003  Double_t ptJet3 = jet3->Pt();
1004  Double_t phiJet3 = jet3->Phi();
1005 
1006  if(phiJet3>2*TMath::Pi()) phiJet3 -= 2*TMath::Pi();
1007  if(phiJet3<-2*TMath::Pi()) phiJet3 += 2*TMath::Pi();
1008  if(phiJet3<-0.5*TMath::Pi()) phiJet3 += 2*TMath::Pi();
1009  if(phiJet3>1.5*TMath::Pi()) phiJet3 -= 2*TMath::Pi();
1010 
1011  if(phiJet1>2*TMath::Pi()) phiJet1 -= 2*TMath::Pi();
1012  if(phiJet1<-2*TMath::Pi()) phiJet1 += 2*TMath::Pi();
1013  if(phiJet1<-0.5*TMath::Pi()) phiJet1 += 2*TMath::Pi();
1014  if(phiJet1>1.5*TMath::Pi()) phiJet1 -= 2*TMath::Pi();
1015 
1016  if(fDebug) Printf("--- jet pt = jet hybrid pt = %f\t jet matched det pt = %f\t jet matched particle level pt = %f\t",ptJet1,ptJet2,ptJet3);
1017 
1018  fhPtDetPart->Fill(ptJet2,ptJet3);
1019  Double_t fraction = jetCont->GetFractionSharedPt(jet1);
1020  if(fraction < fMinFractionSharedPt) continue;
1021 
1022 
1023  Double_t residual = (ptCorr - ptJet3) / ptJet3;
1024  Double_t residualPhi = (phiJet1 - phiJet3) / phiJet3;
1025 
1026  fhPtHybrDet->Fill(ptJet1,ptJet2);
1027  fhPtHybrPart->Fill(ptJet1,ptJet3);
1028  fhPtHybrPartCor->Fill(ptCorr,ptJet3);
1029  fhPhiHybrPartCor->Fill(phiJet1,phiJet3);
1030 
1031  fhResidual->Fill(residual);
1032  fhPtResidual->Fill(ptJet3,residual);
1033  fhPhiResidual->Fill(residualPhi);
1034  fhPhiPhiResidual->Fill(phiJet3,residualPhi);
1035 
1038  fTreeVarsInclusive[1] = ptJet1;
1039  fTreeVarsInclusive[2] = area;
1040  fTreeVarsInclusive[3] = ptCorr;
1041  fTreeVarsInclusive[4] = phiJet1;
1042  fTreeVarsInclusive[5] = ptJet3;
1043  fTreeVarsInclusive[6] = phiJet3;
1045  fTreeEmbInclusive->Fill();
1046  }
1047  }
1048 // for(auto jettrue : jetContPart->accepted()) {
1049 // // jettrue
1050 // Double_t ptTrue = jettrue->Pt();
1051 // Double_t phiTrue = jettrue->Phi();
1052 // Double_t etaTrue= jettrue->Eta();
1053 // auto jetmatched = jettrue->ClosestJet();
1054 // if(!jetmatched) continue;
1055 // Double_t ptMatched = jetmatched->Pt();
1056 // Double_t phiMatched = jetmatched->Phi();
1057 // Double_t etaMatched = jetmatched->Eta();
1058 // Double_t residual = (ptMatched - ptTrue) / ptTrue;
1059 // fhPtHybrTrue->Fill(ptTrue,ptMatched);
1060 // fhResidual->Fill(residual);
1061 // fhPtResidual->Fill(ptTrue,residual);
1062 // }
1063 }
1064 
1065 
1066 
1072 {
1073  TString histname;
1074  TString groupname;
1075  AliJetContainer* jetCont = 0;
1076  TIter next(&fJetCollArray);
1077  while ((jetCont = static_cast<AliJetContainer*>(next()))) {
1078  groupname = jetCont->GetName();
1079  UInt_t count = 0;
1080  for(auto jet : jetCont->accepted()) {
1081  if (!jet) continue;
1082  count++;
1083 
1084  histname = TString::Format("%s/histJetPt_%d", groupname.Data(), fCentBin);
1085  fHistManager.FillTH1(histname, jet->Pt());
1086 
1087  histname = TString::Format("%s/histJetArea_%d", groupname.Data(), fCentBin);
1088  fHistManager.FillTH1(histname, jet->Area());
1089 
1090  histname = TString::Format("%s/histJetPhi_%d", groupname.Data(), fCentBin);
1091  fHistManager.FillTH1(histname, jet->Phi());
1092 
1093  histname = TString::Format("%s/histJetEta_%d", groupname.Data(), fCentBin);
1094  fHistManager.FillTH1(histname, jet->Eta());
1095 
1096  if (jetCont->GetRhoParameter()) {
1097  histname = TString::Format("%s/histJetCorrPt_%d", groupname.Data(), fCentBin);
1098  fHistManager.FillTH1(histname, jet->Pt() - jetCont->GetRhoVal() * jet->Area());
1099  }
1100  }
1101  histname = TString::Format("%s/histNJets_%d", groupname.Data(), fCentBin);
1102  fHistManager.FillTH1(histname, count);
1103  }
1104 }
1105 
1111 {
1113 
1114  TString histname;
1115  TString groupname;
1116  UInt_t sumAcceptedTracks = 0;
1117  AliParticleContainer* partCont = 0;
1118  TIter next(&fParticleCollArray);
1119  while ((partCont = static_cast<AliParticleContainer*>(next()))) {
1120  groupname = partCont->GetName();
1121  UInt_t count = 0;
1122  for(auto part : partCont->accepted()) {
1123  if (!part) continue;
1124  count++;
1125 
1126  histname = TString::Format("%s/histTrackPt_%d", groupname.Data(), fCentBin);
1127  fHistManager.FillTH1(histname, part->Pt());
1128 
1129  histname = TString::Format("%s/histTrackPhi_%d", groupname.Data(), fCentBin);
1130  fHistManager.FillTH1(histname, part->Phi());
1131 
1132  histname = TString::Format("%s/histTrackEta_%d", groupname.Data(), fCentBin);
1133  fHistManager.FillTH1(histname, part->Eta());
1134 
1135  if (partCont->GetLoadedClass()->InheritsFrom("AliVTrack")) {
1136  const AliVTrack* track = static_cast<const AliVTrack*>(part);
1137 
1138  histname = TString::Format("%s/fHistDeltaEtaPt_%d", groupname.Data(), fCentBin);
1139  fHistManager.FillTH1(histname, track->Pt(), track->Eta() - track->GetTrackEtaOnEMCal());
1140 
1141  histname = TString::Format("%s/fHistDeltaPhiPt_%d", groupname.Data(), fCentBin);
1142  fHistManager.FillTH1(histname, track->Pt(), track->Phi() - track->GetTrackPhiOnEMCal());
1143 
1144  histname = TString::Format("%s/fHistDeltaPtvsPt_%d", groupname.Data(), fCentBin);
1145  fHistManager.FillTH1(histname, track->Pt(), track->Pt() - track->GetTrackPtOnEMCal());
1146 
1147  if (clusCont) {
1148  Int_t iCluster = track->GetEMCALcluster();
1149  if (iCluster >= 0) {
1150  AliVCluster* cluster = clusCont->GetAcceptCluster(iCluster);
1151  if (cluster) {
1152  histname = TString::Format("%s/fHistEoverPvsP_%d", groupname.Data(), fCentBin);
1153  fHistManager.FillTH2(histname, track->P(), cluster->GetNonLinCorrEnergy() / track->P());
1154  }
1155  }
1156  }
1157  }
1158  }
1159  sumAcceptedTracks += count;
1160 
1161  histname = TString::Format("%s/histNTracks_%d", groupname.Data(), fCentBin);
1162  fHistManager.FillTH1(histname, count);
1163  }
1164 
1165  histname = "fHistSumNTracks";
1166  fHistManager.FillTH1(histname, sumAcceptedTracks);
1167 }
1168 
1174 {
1175  TString histname;
1176  TString groupname;
1177  UInt_t sumAcceptedClusters = 0;
1178  AliClusterContainer* clusCont = 0;
1179  TIter next(&fClusterCollArray);
1180  while ((clusCont = static_cast<AliClusterContainer*>(next()))) {
1181  groupname = clusCont->GetName();
1182 
1183  for(auto cluster : clusCont->all()) {
1184  if (!cluster) continue;
1185 
1186  if (cluster->GetIsExotic()) {
1187  histname = TString::Format("%s/histClusterEnergyExotic_%d", groupname.Data(), fCentBin);
1188  fHistManager.FillTH1(histname, cluster->E());
1189  }
1190  }
1191 
1192  UInt_t count = 0;
1193  for(auto cluster : clusCont->accepted()) {
1194  if (!cluster) continue;
1195  count++;
1196 
1197  AliTLorentzVector nPart;
1198  cluster->GetMomentum(nPart, fVertex);
1199 
1200  histname = TString::Format("%s/histClusterEnergy_%d", groupname.Data(), fCentBin);
1201  fHistManager.FillTH1(histname, cluster->E());
1202 
1203  histname = TString::Format("%s/histClusterNonLinCorrEnergy_%d", groupname.Data(), fCentBin);
1204  fHistManager.FillTH1(histname, cluster->GetNonLinCorrEnergy());
1205 
1206  histname = TString::Format("%s/histClusterHadCorrEnergy_%d", groupname.Data(), fCentBin);
1207  fHistManager.FillTH1(histname, cluster->GetHadCorrEnergy());
1208 
1209  histname = TString::Format("%s/histClusterPhi_%d", groupname.Data(), fCentBin);
1210  fHistManager.FillTH1(histname, nPart.Phi_0_2pi());
1211 
1212  histname = TString::Format("%s/histClusterEta_%d", groupname.Data(), fCentBin);
1213  fHistManager.FillTH1(histname, nPart.Eta());
1214  }
1215  sumAcceptedClusters += count;
1216 
1217  histname = TString::Format("%s/histNClusters_%d", groupname.Data(), fCentBin);
1218  fHistManager.FillTH1(histname, count);
1219  }
1220 
1221  histname = "fHistSumNClusters";
1222  fHistManager.FillTH1(histname, sumAcceptedClusters);
1223 }
1224 
1230 {
1231  if (!fCaloCells) return;
1232 
1233  TString histname;
1234 
1235  const Short_t ncells = fCaloCells->GetNumberOfCells();
1236 
1237  histname = TString::Format("%s/histNCells_%d", fCaloCellsName.Data(), fCentBin);
1238  fHistManager.FillTH1(histname, ncells);
1239 
1240  histname = TString::Format("%s/histCellEnergy_%d", fCaloCellsName.Data(), fCentBin);
1241  for (Short_t pos = 0; pos < ncells; pos++) {
1242  Double_t amp = fCaloCells->GetAmplitude(pos);
1243 
1244  fHistManager.FillTH1(histname, amp);
1245  }
1246 }
1247 
1253 {
1255 }
1256 
1265 {
1266  return kTRUE;
1267 }
1268 
1273 {
1274 }
1275 
1277 
1278  Int_t index=-1;
1279  Int_t triggers[100];
1280 
1281  for(Int_t cr=0;cr<100;cr++) triggers[cr]=-1;
1282 
1283  Int_t im=0;
1284 
1285  TString groupname = "";
1286  AliParticleContainer* partCont = 0x0;
1288  else partCont = GetParticleContainer(0);
1289  groupname = partCont->GetName();
1290  UInt_t iCount = 0;
1291  for(auto part : partCont->accepted()) {
1292  if (!part) continue;
1293  list->Add(part);
1294  iCount++;
1295  if(part->Pt()>=minT && part->Pt()<maxT){
1296  triggers[im]=iCount-1;
1297  im=im+1;
1298  }
1299  }
1300  number=im;
1301  Int_t rd=0;
1302  if(im>0) rd=fRandom->Integer(im);
1303  index=triggers[rd];
1304  return index;
1305 }
1306 
1308 
1309  // get relative DeltaPhi from -pi < DeltaPhi < pi
1310 
1311  if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1312  else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1313  if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1314  else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1315  double dphi = mphi-vphi;
1316  if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1317  else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1318  return dphi;//dphi in [-Pi, Pi]
1319 }
1320 
1322 {
1323  Int_t phibin=-1;
1324  if(!(TMath::Abs(phi)<=2*TMath::Pi())) return -1;
1325  Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
1326  phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
1327  //if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
1328  return phibin;
1329 }
Float_t fJetHadronDeltaPhi
max angle from pi (set <0 for no selection)
THashList * CreateHistoGroup(const char *groupname)
Create a new group of histograms within a parent group.
TObjArray fClusterCollArray
cluster collection array
Float_t fTTLowRef
minimum reference trigger track pt
Double_t GetRhoVal() const
const TString & GetRhoName() const
double Double_t
Definition: External.C:58
Float_t fFrac
fraction of events that are used to fill signal recoil jet population
Bool_t fFillRecoilTHnSparse
switch to fill recoil THnSparse for main analysis
Definition: External.C:236
Int_t SelectTrigger(TList *list, Double_t minT, Double_t maxT, Int_t &number)
AliJetContainer * GetJetContainer(Int_t i=0) const
Float_t fTTUpRef
maximum reference trigger track pt
Int_t GetNParticles() const
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
Fill a 2D histogram within the container.
Declaration of class AliTLorentzVector.
Bool_t fFillInclusiveTree
switch to fill embedding tree with inclusive jet info
Double_t fMinBinPt
min pt in histograms
Double_t fEPV0
!event plane V0
Int_t fCentBin
!event centrality bin
Float_t fMinEmbJetPt
min corrected jet pt to use in embedding
TString fJetContName
Base level jet container name.
Float_t fMinFractionSharedPt
min fraction of pt between hybrid / detector jets
Container for particles within the EMCAL framework.
TString fJetContTrueName
True pp (detector) level jet container name.
TString kData
Declare data MC or deltaAOD.
TObjArray fParticleCollArray
particle/track collection array
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
const AliClusterIterableContainer all() const
TH2 * CreateTH2(const char *name, const char *title, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax, Option_t *opt="")
Create a new TH2 within the container.
TObject * FindObject(const char *name) const
Find an object inside the container.
int Int_t
Definition: External.C:63
Float_t fTTUpSig
maximum signal trigger track pt
Definition: External.C:204
unsigned int UInt_t
Definition: External.C:33
THashList * GetListOfHistograms() const
Get the list of histograms.
Definition: THistManager.h:671
float Float_t
Definition: External.C:68
Double_t Phi_0_2pi() const
Float_t fTTLowSig
minimum signal trigger track pt
Int_t GetNJets() const
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
Create a new TH1 within the container.
BeamType fForceBeamType
forced beam type
Int_t fNcentBins
how many centrality bins
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Get cluster container attached to this task.
Definition: External.C:212
AliVCluster * GetAcceptCluster(Int_t i) const
TString fJetContPartName
Particle(MC) level jet container name.
const AliClusterIterableContainer accepted() const
Double_t fCent
!event centrality
Float_t fCentMax
maximum centrality
TString fCaloCellsName
name of calo cell collection
Bool_t fFillJetHistograms
switch to fill jet histograms
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
Fill a 1D histogram within the container.
TObjArray fJetCollArray
jet collection array
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
Double_t RelativePhi(Double_t mphi, Double_t vphi)
short Short_t
Definition: External.C:23
AliVCaloCells * fCaloCells
!cells
AliRhoParameter * GetRhoParameter()
AliEmcalList * fOutput
!output list
Double_t fMaxBinPt
max pt in histograms
Int_t fPtHardBin
pt hard bin if running embedding
THistManager fHistManager
Histogram manager.
Double_t fVertex[3]
!event vertex
void SetMakeGeneralHistograms(Bool_t g)
Base task in the EMCAL jet framework.
Double_t GetFractionSharedPt(const AliEmcalJet *jet, AliParticleContainer *cont2=0x0) const
Declaration of class AliEmcalPythiaInfo.
const AliParticleIterableContainer accepted() const
const char Option_t
Definition: External.C:48
void UserCreateOutputObjects()
Main initialization function on the worker.
const AliJetIterableContainer accepted() const
bool Bool_t
Definition: External.C:53
Bool_t fFillRecoilTree
switch to fill embedding tree with recoil jet info
Container structure for EMCAL clusters.
Float_t fCentMin
minimum centrality
Container for jet within the EMCAL jet framework.
Definition: External.C:196
Bool_t fFillTrackHistograms
switch to fill track histograms
Int_t fNbins
no. of pt bins
Int_t fRejectionFactorInclusiveJets
factor to reject inclusive jets, to reduce size of ttree