AliPhysics  master (3d17d9d)
AliAnalysisTaskHardestBranch.cxx
Go to the documentation of this file.
1 //
2 // Checking the hardest branch
3 // Leticia Cunqueiro
4 //
5 //
6 #include "AliAODMCHeader.h"
7 #include "AliAnalysisManager.h"
8 #include "AliEmcalJet.h"
9 #include "AliEmcalParticle.h"
10 #include "AliEmcalPythiaInfo.h"
11 #include "AliGenPythiaEventHeader.h"
12 #include "AliJetContainer.h"
13 #include "AliLog.h"
14 #include "AliMCEvent.h"
15 #include "AliParticleContainer.h"
16 #include "AliRhoParameter.h"
17 #include "AliVCluster.h"
18 #include "AliVTrack.h"
19 #include "TMatrixD.h"
20 #include "TMatrixDSym.h"
21 #include "TMatrixDSymEigen.h"
22 #include "TRandom3.h"
23 #include "TVector2.h"
24 #include "TVector3.h"
25 #include <AliAnalysisDataContainer.h>
26 #include <AliAnalysisDataSlot.h>
27 #include <vector>
28 #include <TChain.h>
29 #include <TClonesArray.h>
30 #include <TFile.h>
31 #include <TH1F.h>
32 #include <TH2F.h>
33 #include <TH3F.h>
34 #include <THnSparse.h>
35 #include <TKey.h>
36 #include <TList.h>
37 #include <TLorentzVector.h>
38 #include <TProfile.h>
39 #include <TSystem.h>
40 #include <TTree.h>
41 #include <iostream> // std::cout
42 #include <algorithm>
43 
44 #include "AliAODEvent.h"
46 
47 using std::cout;
48 using std::endl;
49 
51 
52  //________________________________________________________________________
55  fContainer(0), fMinFractionShared(0), fJetShapeType(kData),
56  fJetShapeSub(kNoSub), fJetSelection(kInclusive), fPtThreshold(-9999.),
57  fRMatching(0.2), fCentSelectOn(kTRUE), fCentMin(0), fCentMax(10),
58  fOneConstSelectOn(kFALSE), fTrackCheckPlots(kFALSE),
59  fDoFillMCLund(kFALSE), fCheckResolution(kFALSE), fSubjetCutoff(0.1),
60  fMinPtConst(1), fHardCutoff(0), fDoTwoTrack(kFALSE),
61  fDoAreaIterative(kTRUE), fPowerAlgo(1), fPhiCutValue(0.02),
62  fEtaCutValue(0.02), fMagFieldPolarity(1), fDerivSubtrOrder(0),
63  fPtJet(0x0),fTreeSubstructure(0)
64 
65 {
66  for (Int_t i = 0; i < 10; i++) {
67  fShapesVar[i] = 0;
68  }
69  SetMakeGeneralHistograms(kTRUE);
70  DefineOutput(1, TList::Class());
71  DefineOutput(2, TTree::Class());
72 }
73 
74 //________________________________________________________________________
76  const char *name)
77  : AliAnalysisTaskEmcalJet(name, kTRUE), fContainer(0),
78  fMinFractionShared(0), fJetShapeType(kData), fJetShapeSub(kNoSub),
79  fJetSelection(kInclusive), fPtThreshold(-9999.), fRMatching(0.2),
80  fCentSelectOn(kTRUE), fCentMin(0), fCentMax(10),
81  fOneConstSelectOn(kFALSE), fTrackCheckPlots(kFALSE),
82  fDoFillMCLund(kFALSE), fCheckResolution(kFALSE), fSubjetCutoff(0.1),
83  fMinPtConst(1), fHardCutoff(0), fDoTwoTrack(kFALSE),
84  fDoAreaIterative(kTRUE), fPowerAlgo(1), fPhiCutValue(0.02),
85  fEtaCutValue(0.02), fMagFieldPolarity(1), fDerivSubtrOrder(0),
86  fPtJet(0x0),fTreeSubstructure(0)
87 
88 {
89  // Standard constructor.
90  for (Int_t i = 0; i < 10; i++) {
91  fShapesVar[i] = 0;
92  }
94 
95  DefineOutput(1, TList::Class());
96  DefineOutput(2, TTree::Class());
97 }
98 
99 //________________________________________________________________________
101  // Destructor.
102 }
103 
104 //________________________________________________________________________
106  // Create user output.
107 
109 
110  Bool_t oldStatus = TH1::AddDirectoryStatus();
111  TH1::AddDirectory(kFALSE);
112 
113  fPtJet = new TH1F("fPtJet", "fPtJet", 100, 0, 200);
114  fOutput->Add(fPtJet);
115 
116 
117 
118  // =========== Switch on Sumw2 for all histos ===========
119  for (Int_t i = 0; i < fOutput->GetEntries(); ++i) {
120  TH1 *h1 = dynamic_cast<TH1 *>(fOutput->At(i));
121  if (h1) {
122  h1->Sumw2();
123  continue;
124  }
125  THnSparse *hn = dynamic_cast<THnSparse *>(fOutput->At(i));
126  if (hn)
127  hn->Sumw2();
128  }
129 
130  TH1::AddDirectory(oldStatus);
131  const Int_t nVar = 13;
132  const char *nameoutput = GetOutputSlot(2)->GetContainer()->GetName();
133  fTreeSubstructure = new TTree(nameoutput, nameoutput);
134  TString *fShapesVarNames = new TString[nVar];
135 
136  fShapesVarNames[0] = "ptJet";
137  fShapesVarNames[1] = "ktg";
138  fShapesVarNames[2] = "tfg";
139  fShapesVarNames[3] = "zg";
140  fShapesVarNames[4] = "rg";
141  fShapesVarNames[5] = "ng";
142  fShapesVarNames[6] = "ptJetMatch";
143  fShapesVarNames[7] = "ktgMatch";
144  fShapesVarNames[8] = "tfgMatch";
145  fShapesVarNames[9] = "zgMatch";
146  fShapesVarNames[10] = "rgMatch";
147  fShapesVarNames[11] = "ngMatch";
148  fShapesVarNames[12] = "LeadingTrackPtMatch";
149 
150 
151  for (Int_t ivar = 0; ivar < nVar; ivar++) {
152  cout << "looping over variables" << endl;
153  fTreeSubstructure->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar],
154  Form("%s/F", fShapesVarNames[ivar].Data()));
155  }
156 
157  PostData(1, fOutput);
158  PostData(2, fTreeSubstructure);
159 
160  delete[] fShapesVarNames;
161 }
162 
163 //________________________________________________________________________
165  // Run analysis code here, if needed. It will be executed before
166  // FillHistograms().
167 
168  return kTRUE;
169 }
170 
171 //________________________________________________________________________
173 
174  AliEmcalJet *jet1 = NULL;
175  AliJetContainer *jetCont = GetJetContainer(0);
176  // container zero is always the base containe: the data container, the
177  // embedded subtracted in the case of embedding or the detector level in case
178  // of pythia
179 
180  if (fCentSelectOn)
181  if ((fCent > fCentMax) || (fCent < fCentMin))
182  return 0;
183 
184  Float_t rhoVal = 0, rhoMassVal = 0.;
185  if (jetCont) {
186 
187  jetCont->ResetCurrentID();
188  if ((fJetShapeSub == kConstSub) || (fJetShapeSub == kDerivSub)) {
189  // rho
190  AliRhoParameter *rhoParam = dynamic_cast<AliRhoParameter *>(
191  InputEvent()->FindListObject("RhoSparseR020"));
192  if (!rhoParam) {
193  Printf("%s: Could not retrieve rho %s (some histograms will be filled "
194  "with zero)!",
195  GetName(), jetCont->GetRhoName().Data());
196  } else
197  rhoVal = rhoParam->GetVal();
198  // rhom
199  AliRhoParameter *rhomParam = dynamic_cast<AliRhoParameter *>(
200  InputEvent()->FindListObject("RhoMassSparseR020"));
201  if (!rhomParam) {
202  Printf("%s: Could not retrieve rho_m %s (some histograms will be "
203  "filled with zero)!",
204  GetName(), jetCont->GetRhoMassName().Data());
205  } else
206  rhoMassVal = rhomParam->GetVal();
207  }
208 
209  while ((jet1 = jetCont->GetNextAcceptJet())) {
210  if (!jet1)
211  continue;
212  AliEmcalJet *jet2 = 0x0;
213  AliEmcalJet *jet3 = 0x0;
214  fPtJet->Fill(jet1->Pt());
215  AliEmcalJet *jetUS = NULL;
216  Int_t ifound = 0, jfound = 0;
217  Int_t ilab = -1, jlab = -1;
218 
219  // The embedding mode
220  // the matching is done between unsubtracted embedded jets and detector
221  // level jets unsubtracted and subtracted jets share the label. Once we
222  // identify the corresponding unsubtracted jet, jetUS, then we fetch jet2,
223  // which is the matched detector level jet In the case we are not
224  // considering constituent subtraction, then the detector-level matched jet
225  // is the one that was directly matched to the base jet1. Then, the
226  // particle-level jet jet3 is obtained as the matched one to jet2 In short,
227  // there are 2 consecutive matchinges, between particle-level (jet3) and
228  // detector-level (jet2) pythia jets and between jet2 and the embedding
229  // unsubtracted jet. Note that the matching obtained via ClosestJet is
230  // purely geometrical. So below we require a fraction of the probe momentum
231  // to be reconstructed in the embedded jet.
233 
234  AliJetContainer *jetContUS = GetJetContainer(2);
235 
236  if (fJetShapeSub == kConstSub) {
237 
238  for (Int_t i = 0; i < jetContUS->GetNJets(); i++) {
239  jetUS = jetContUS->GetJet(i);
240  if (jetUS->GetLabel() == jet1->GetLabel()) {
241  ifound++;
242  if (ifound == 1)
243  ilab = i;
244  }
245  }
246  if (ilab == -1)
247  continue;
248  jetUS = jetContUS->GetJet(ilab);
249  jet2 = jetUS->ClosestJet();
250  }
251 
252  if (!(fJetShapeSub == kConstSub))
253  jet2 = jet1->ClosestJet();
254  if (!jet2) {
255  Printf("jet2 does not exist, returning");
256  continue;
257  }
258 
259  // AliJetContainer *jetContPart=GetJetContainer(3);
260  jet3 = jet2->ClosestJet();
261 
262  if (!jet3) {
263  Printf("jet3 does not exist, returning");
264  continue;
265  }
266  cout << "jet 3 exists" << jet3->Pt() << endl;
267 
268  Double_t fraction = 0;
269  if (!(fJetShapeSub == kConstSub))
270  fraction = jetCont->GetFractionSharedPt(jet1);
271  if (fJetShapeSub == kConstSub)
272  fraction = jetContUS->GetFractionSharedPt(jetUS);
273 
274  if (fraction < fMinFractionShared)
275  continue;
276  }
277 
278  // this is the mode to run over pythia to produce a det-part response
279  // here we have also added the constituent-subtraction case, but we don't
280  // use it normally in pp the matching is purely geometrical
281  if (fJetShapeType == kPythiaDef) {
282 
283  AliJetContainer *jetContTrue = GetJetContainer(1);
284  AliJetContainer *jetContUS = GetJetContainer(2);
285  AliJetContainer *jetContPart = GetJetContainer(3);
286 
287  if (fJetShapeSub == kConstSub) {
288 
289  for (Int_t i = 0; i < jetContUS->GetNJets(); i++) {
290  jetUS = jetContUS->GetJet(i);
291  if (jetUS->GetLabel() == jet1->GetLabel()) {
292  ifound++;
293  if (ifound == 1)
294  ilab = i;
295  }
296  }
297  if (ilab == -1)
298  continue;
299  jetUS = jetContUS->GetJet(ilab);
300  jet2 = jetUS->ClosestJet();
301 
302  if (!jet2) {
303  Printf("jet2 does not exist, returning");
304  continue;
305  }
306 
307  for (Int_t j = 0; j < jetContPart->GetNJets(); j++) {
308 
309  jet3 = jetContPart->GetJet(j);
310  if (!jet3)
311  continue;
312  if (jet3->GetLabel() == jet2->GetLabel()) {
313  jfound++;
314  if (jfound == 1)
315  jlab = j;
316  }
317  }
318  if (jlab == -1)
319  continue;
320  jet3 = jetContPart->GetJet(jlab);
321  if (!jet3) {
322  Printf("jet3 does not exist, returning");
323  continue;
324  }
325  }
326  if (!(fJetShapeSub == kConstSub))
327  jet3 = jet1->ClosestJet();
328  if (!jet3) {
329  Printf("jet3 does not exist, returning");
330  continue;
331  }
332 
333  if (fCheckResolution)
334  CheckSubjetResolution(jet1, jetCont, jet3, jetContTrue);
335  }
336 
337  Double_t ptSubtracted = 0;
338  if (fJetShapeSub == kConstSub)
339  ptSubtracted = jet1->Pt();
340 
341  else if (fJetShapeSub == kDerivSub) {
342  ptSubtracted = jet1->Pt() - GetRhoVal(0) * jet1->Area();
343  }
344 
345  else if (fJetShapeSub == kNoSub) {
347  ptSubtracted = jet1->Pt() - GetRhoVal(0) * jet1->Area();
348  else if ((fJetShapeType == kPythiaDef) || (fJetShapeType == kMCTrue) ||
350  ptSubtracted = jet1->Pt();
351  }
352 
353  if (ptSubtracted < fPtThreshold)
354  continue;
355 
356  if ((fCentSelectOn == kFALSE) && (jet1->GetNumberOfTracks() <= 1))
357  continue;
358 
359  fShapesVar[0] = ptSubtracted;
360  fShapesVar[10] = jet1->MaxTrackPt();
361  IterativeParents(jet1, jetCont);
362 
363  Float_t ptMatch = 0.;
364  Float_t leadTrackMatch = 0.;
365  Double_t ktgMatch = 0;
366  Double_t tfgMatch =0;
367  Double_t ngMatch = 0;
368  Double_t zgMatch = 0;
369  Double_t rgMatch = 0;
370  Float_t ptDet = 0.;
371 
372 
373  Double_t aver1 = 0;
374  Double_t aver2 = 0;
375  Double_t aver3 = 0;
376  Double_t aver4 = 0;
377  Double_t aver5 = 0;
378  Int_t kMatched = 0;
379  if (fJetShapeType == kPythiaDef) {
380  kMatched = 1;
381  if (fJetShapeSub == kConstSub)
382  kMatched = 3;
383 
384  ptMatch = jet3->Pt();
385  leadTrackMatch = jet3->MaxTrackPt();
386  IterativeParentsMCAverage(jet3, kMatched, aver1, aver2, aver3, aver4, aver5);
387  ktgMatch = aver1;
388  tfgMatch = aver2;
389  zgMatch = aver3;
390  rgMatch = aver4;
391  ngMatch = aver5;
392  }
393 
395  if (fJetShapeSub == kConstSub)
396  kMatched = 3;
397  if (fJetShapeSub == kDerivSub)
398  kMatched = 2;
399  ptMatch = jet3->Pt();
400  leadTrackMatch = jet3->MaxTrackPt();
401  IterativeParentsMCAverage(jet3, kMatched, aver1, aver2, aver3, aver4, aver5);
402  ktgMatch = aver1;
403  tfgMatch = aver2;
404  zgMatch = aver3;
405  rgMatch = aver4;
406  ngMatch = aver5;
407  }
408 
409  if (fJetShapeType == kMCTrue || fJetShapeType == kData ||
411 
412  ptMatch = 0.;
413  leadTrackMatch = 0.;
414  ktgMatch = 0.;
415  tfgMatch = 0.;
416  zgMatch = 0;
417  rgMatch = 0;
418  ngMatch = 0;
419  }
420 
421  fShapesVar[6] = ptMatch;
422  fShapesVar[7] = ktgMatch;
423  fShapesVar[8] = tfgMatch;
424  fShapesVar[9] = zgMatch;
425  fShapesVar[10] = rgMatch;
426  fShapesVar[11] = ngMatch;
427  fShapesVar[12] = leadTrackMatch;
428 
429 
430  fTreeSubstructure->Fill();
431  }
432  }
433 
434  return kTRUE;
435 }
436 
437 //________________________________________________________________________
439  Int_t jetContNb = 0) {
440  // calc subtracted jet mass
441  if ((fJetShapeSub == kDerivSub) && (jetContNb == 0))
442  if (fDerivSubtrOrder == 1)
444  else
446  else
447  return jet->M();
448 }
449 
450 //________________________________________________________________________
452  Int_t jetContNb = 0) {
453 
454  AliJetContainer *jetCont = GetJetContainer(jetContNb);
455  if (!jet->GetNumberOfTracks())
456  return 0;
457  Double_t den = 0.;
458  Double_t num = 0.;
459  AliVParticle *vp1 = 0x0;
460  for (UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
461  vp1 = static_cast<AliVParticle *>(
462  jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
463 
464  if (!vp1) {
465  Printf("AliVParticle associated to constituent not found");
466  continue;
467  }
468 
469  Double_t dphi = RelativePhi(vp1->Phi(), jet->Phi());
470  Double_t dr2 =
471  (vp1->Eta() - jet->Eta()) * (vp1->Eta() - jet->Eta()) + dphi * dphi;
472  Double_t dr = TMath::Sqrt(dr2);
473  num = num + vp1->Pt() * dr;
474  den = den + vp1->Pt();
475  }
476  return num / den;
477 }
478 
479 //________________________________________________________________________
480 Float_t
482  Int_t jetContNb = 0) {
483 
484  if ((fJetShapeSub == kDerivSub) && (jetContNb == 0))
485  if (fDerivSubtrOrder == 1)
487  else
489  else
490  return Angularity(jet, jetContNb);
491 }
492 
493 //__________________________________________________________________________________
495  Double_t vphi) {
496 
497  if (vphi < -1 * TMath::Pi())
498  vphi += (2 * TMath::Pi());
499  else if (vphi > TMath::Pi())
500  vphi -= (2 * TMath::Pi());
501  if (mphi < -1 * TMath::Pi())
502  mphi += (2 * TMath::Pi());
503  else if (mphi > TMath::Pi())
504  mphi -= (2 * TMath::Pi());
505  double dphi = mphi - vphi;
506  if (dphi < -1 * TMath::Pi())
507  dphi += (2 * TMath::Pi());
508  else if (dphi > TMath::Pi())
509  dphi -= (2 * TMath::Pi());
510  return dphi; // dphi in [-Pi, Pi]
511 }
512 
513 
514 //_________________________________________________________________________
516 
517  std::vector<fastjet::PseudoJet> fInputVectors;
518  fInputVectors.clear();
519  fastjet::PseudoJet PseudoTracks;
520 
521  AliParticleContainer *fTrackCont = fJetCont->GetParticleContainer();
522 
523  if (fTrackCont)
524  for (Int_t i = 0; i < fJet->GetNumberOfTracks(); i++) {
525  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
526  if (!fTrk)
527  continue;
528  if (fDoTwoTrack == kTRUE && CheckClosePartner(i, fJet, fTrk, fTrackCont))
529  continue;
530  PseudoTracks.reset(fTrk->Px(), fTrk->Py(), fTrk->Pz(), fTrk->E());
531  PseudoTracks.set_user_index(fJet->TrackAt(i) + 100);
532  fInputVectors.push_back(PseudoTracks);
533  }
534  fastjet::JetAlgorithm jetalgo(fastjet::genkt_algorithm);
535  fastjet::GhostedAreaSpec ghost_spec(1, 1, 0.05);
536 
537  fastjet::JetDefinition fJetDef(jetalgo, 1., fPowerAlgo,
538  static_cast<fastjet::RecombinationScheme>(0),
539  fastjet::BestFJ30);
540  fastjet::AreaDefinition fAreaDef(fastjet::passive_area, ghost_spec);
541  try {
542  fastjet::ClusterSequenceArea fClustSeqSA(fInputVectors, fJetDef, fAreaDef);
543  std::vector<fastjet::PseudoJet> fOutputJets;
544  fOutputJets.clear();
545  fOutputJets = fClustSeqSA.inclusive_jets(0);
546 
547  fastjet::PseudoJet jj;
548  fastjet::PseudoJet j1;
549  fastjet::PseudoJet j2;
550  jj = fOutputJets[0];
551 
552  double nall = 0;
553  double nsd = 0;
554  std::vector<Double_t> zvec;
555  std::vector<Double_t> ktvec;
556  std::vector<Double_t> thetavec;
557  std::vector<Double_t> tformvec;
558  std::vector<Double_t> nvec;
559  int indmax;
560  while (jj.has_parents(j1, j2)) {
561  nall = nall + 1;
562 
563  if (j1.perp() < j2.perp())
564  swap(j1, j2);
565 
566  double delta_R = j1.delta_R(j2);
567  double xkt = j2.perp() * sin(delta_R);
568 
569  double z = j2.perp() / (j2.perp() + j1.perp());
570  double form = 2 * 0.197 * j2.e() / ((1.-z)*xkt * xkt);
571 
572  if(z>0.1) nsd=nsd+1;
573  nvec.push_back(nsd);
574  zvec.push_back(z);
575  ktvec.push_back(xkt);
576  thetavec.push_back(delta_R);
577  tformvec.push_back(form);
578  jj = j1;
579  }
580  if(nall>0){
581 
582  auto result = std::max_element(ktvec.begin(), ktvec.end());
583  indmax=std::distance(ktvec.begin(), result);
584 
585 
586 
587  fShapesVar[1] = ktvec[indmax];
588  fShapesVar[2] = tformvec[indmax];
589  fShapesVar[3] = zvec[indmax];
590  fShapesVar[4] = thetavec[indmax];
591  fShapesVar[5] = nvec[indmax];}
592 
593  if(nall==0){
594 
595  fShapesVar[1] = -1;
596  fShapesVar[2] = -1;;
597  fShapesVar[3] = -1;;
598  fShapesVar[4] = -1;;
599  fShapesVar[5] = -1;}
600 
601  } catch (fastjet::Error) {
602  AliError(" [w] FJ Exception caught.");
603  // return -1;
604  }
605 
606  return;
607 
608 }
609 //_________________________________________________________________________
610 void AliAnalysisTaskHardestBranch::IterativeParentsMCAverage(AliEmcalJet *fJet, Int_t km, Double_t &average1, Double_t &average2, Double_t &average3, Double_t &average4, Double_t &average5) {
611  AliJetContainer *jetCont = GetJetContainer(km);
612  std::vector<fastjet::PseudoJet> fInputVectors;
613  fInputVectors.clear();
614  fastjet::PseudoJet PseudoTracks;
615 
616  AliParticleContainer *fTrackCont = jetCont->GetParticleContainer();
617 
618  if (fTrackCont)
619  for (Int_t i = 0; i < fJet->GetNumberOfTracks(); i++) {
620  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
621  if (!fTrk)
622  continue;
623 
624  PseudoTracks.reset(fTrk->Px(), fTrk->Py(), fTrk->Pz(), fTrk->E());
625  PseudoTracks.set_user_index(fJet->TrackAt(i) + 100);
626  fInputVectors.push_back(PseudoTracks);
627  }
628  fastjet::JetAlgorithm jetalgo(fastjet::cambridge_algorithm);
629 
630  fastjet::JetDefinition fJetDef(jetalgo, 1.,
631  static_cast<fastjet::RecombinationScheme>(0),
632  fastjet::BestFJ30);
633 
634  try {
635  fastjet::ClusterSequence fClustSeqSA(fInputVectors, fJetDef);
636  std::vector<fastjet::PseudoJet> fOutputJets;
637  fOutputJets.clear();
638  fOutputJets = fClustSeqSA.inclusive_jets(0);
639 
640  fastjet::PseudoJet jj;
641  fastjet::PseudoJet j1;
642  fastjet::PseudoJet j2;
643  jj = fOutputJets[0];
644  int flagSubjet = 0;
645  double nall = 0;
646  double nsd = 0;
647 
648  std::vector<Double_t> zvec;
649  std::vector<Double_t> ktvec;
650  std::vector<Double_t> thetavec;
651  std::vector<Double_t> tformvec;
652  std::vector<Double_t> nvec;
653  int indmax;
654 
655 
656  while (jj.has_parents(j1, j2)) {
657  nall = nall + 1;
658 
659  if (j1.perp() < j2.perp())
660  swap(j1, j2);
661  double delta_R = j1.delta_R(j2);
662  double xkt = j2.perp() * sin(delta_R);
663 
664  double z = j2.perp() / (j2.perp() + j1.perp());
665  double form = 2 * 0.197 * j2.e() / ((1.-z)*xkt * xkt);
666  double rad = j2.e();
667 
668  if(z>0.1) nsd=nsd+1;
669 
670  zvec.push_back(z);
671  ktvec.push_back(xkt);
672  thetavec.push_back(delta_R);
673  tformvec.push_back(form);
674  nvec.push_back(nsd);
675  jj = j1;
676  }
677 
678  if(nall>0){
679  auto result = std::max_element(ktvec.begin(), ktvec.end());
680  indmax=std::distance(ktvec.begin(), result);
681 
682 
683  average1 = ktvec[indmax];
684  average2 = tformvec[indmax];
685  average3 = zvec[indmax];
686  average4 = thetavec[indmax];
687  average5 = nvec[indmax];}
688 
689  if(nall==0){
690  average1 = -1;
691  average2 = -1;
692  average3 = -1;
693  average4 = -1;
694  average5 = -1;}
695 
696 
697 
698  } catch (fastjet::Error) {
699  AliError(" [w] FJ Exception caught.");
700  // return -1;
701  }
702 
703  return;
704 }
705 
706 
707 //________________________________________________________________________
709  //
710  // retrieve event objects
711  //
713  return kFALSE;
714 
715  return kTRUE;
716 }
717 
718 //_______________________________________________________________________
720  // Called once at the end of the analysis.
721 
722  // fTreeObservableTagging = dynamic_cast<TTree*>(GetOutputData(1));
723  // if (!fTreeObservableTagging){
724  // Printf("ERROR: fTreeObservableTagging not available");
725  // return;
726  // }
727 }
Double_t GetFirstOrderSubtractedAngularity() const
Double_t Area() const
Definition: AliEmcalJet.h:130
const TString & GetRhoName() const
double Double_t
Definition: External.C:58
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:341
AliJetContainer * GetJetContainer(Int_t i=0) const
void IterativeParentsMCAverage(AliEmcalJet *fJet, Int_t km, Double_t &aver1, Double_t &aver2, Double_t &aver3, Double_t &aver4, Double_t &aver5)
Double_t Eta() const
Definition: AliEmcalJet.h:121
Bool_t FillHistograms()
Function filling histograms.
Double_t Phi() const
Definition: AliEmcalJet.h:117
void CheckSubjetResolution(AliEmcalJet *fJet, AliJetContainer *fJetCont, AliEmcalJet *fJetM, AliJetContainer *fJetContM)
Int_t GetLabel() const
Definition: AliEmcalJet.h:124
Float_t Angularity(AliEmcalJet *jet, Int_t jetContNb)
Container for particles within the EMCAL framework.
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:160
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:139
TString kData
Declare data MC or deltaAOD.
Bool_t RetrieveEventObjects()
Retrieve common objects from event.
const TString & GetRhoMassName() const
AliParticleContainer * GetParticleContainer() const
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
Int_t GetNJets() const
Double_t GetSecondOrderSubtractedAngularity() const
Double_t RelativePhi(Double_t mphi, Double_t vphi)
Double_t MaxTrackPt() const
Definition: AliEmcalJet.h:155
Double_t fCent
!event centrality
TClonesArray * GetArray() const
Float_t GetJetMass(AliEmcalJet *jet, Int_t jetContNb)
AliEmcalJet * GetNextAcceptJet()
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)
void IterativeParents(AliEmcalJet *fJet, AliJetContainer *fJetCont)
Double_t Pt() const
Definition: AliEmcalJet.h:109
Double_t GetRhoVal(Int_t i=0) const
Bool_t CheckClosePartner(Int_t index, AliEmcalJet *fJet, AliVParticle *fTrack, AliParticleContainer *fTrackCont)
AliEmcalList * fOutput
!output list
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.
void SetMakeGeneralHistograms(Bool_t g)
Enable general histograms.
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
Double_t GetFractionSharedPt(const AliEmcalJet *jet, AliParticleContainer *cont2=0x0) const
Declaration of class AliEmcalPythiaInfo.
const char Option_t
Definition: External.C:48
void UserCreateOutputObjects()
Main initialization function on the worker.
Float_t GetJetAngularity(AliEmcalJet *jet, Int_t jetContNb)
bool Bool_t
Definition: External.C:53
AliEmcalJetShapeProperties * GetShapeProperties() const
Definition: AliEmcalJet.h:375
void swap(PWGJE::EMCALJetTasks::AliAnalysisTaskEmcalJetHPerformance &first, PWGJE::EMCALJetTasks::AliAnalysisTaskEmcalJetHPerformance &second)
Double_t M() const
Definition: AliEmcalJet.h:120
void ResetCurrentID(Int_t i=-1)
Reset the iterator to a given index.
Container for jet within the EMCAL jet framework.
Definition: External.C:196
AliEmcalJet * GetJet(Int_t i) const