AliPhysics  9c66e61 (9c66e61)
AliAnalysisTaskEmcalJetShapesMC.cxx
Go to the documentation of this file.
1 //
2 // Jet QG tagging analysis task.
3 //
4 // Author: D. Caffarri, L. Cunqueiro
5 
6 #include <TClonesArray.h>
7 #include <TH1F.h>
8 #include <TH2F.h>
9 #include <TH3F.h>
10 #include <THnSparse.h>
11 #include <TTree.h>
12 #include <TList.h>
13 #include <TLorentzVector.h>
14 #include <TProfile.h>
15 #include <TChain.h>
16 #include <TSystem.h>
17 #include <TFile.h>
18 #include <TKey.h>
19 #include <AliAnalysisDataSlot.h>
20 #include <AliAnalysisDataContainer.h>
21 #include "TMatrixD.h"
22 #include "TMatrixDSym.h"
23 #include "TMatrixDSymEigen.h"
24 #include "TVector3.h"
25 #include "TVector2.h"
26 #include "AliVCluster.h"
27 #include "AliVTrack.h"
28 #include "AliEmcalJet.h"
29 #include "AliRhoParameter.h"
30 #include "AliLog.h"
31 #include "AliEmcalParticle.h"
32 #include "AliMCEvent.h"
33 #include "AliGenPythiaEventHeader.h"
34 #include "AliAODMCHeader.h"
35 #include "AliMCEvent.h"
36 #include "AliAnalysisManager.h"
37 #include "AliJetContainer.h"
38 #include "AliParticleContainer.h"
39 #include "AliEmcalPythiaInfo.h"
40 #include "TRandom3.h"
41 #include "TF1.h"
42 #include "AliEmcalJetFinder.h"
43 #include "AliAODEvent.h"
45 
46 using std::cout;
47 using std::endl;
48 
50 
51 //________________________________________________________________________
54  fContainer(0),
55  fMinFractionShared(0),
56  fJetShapeType(kGenShapes),
57  fJetShapeSub(kNoSub),
58  fJetSelection(kInclusive),
59  fPtThreshold(-9999.),
60  fRMatching(0.2),
61  fJetRadius(0.4),
62  fSubjetRadius(0.2),
63  fSelectedShapes(0),
64  fSwitchKtNSub(0),
65  fSwitchMinNSub(0),
66  fSwitchAktNSub(0),
67  fSwitchSDKtNSub(0),
68  fSwitchSDMinNSub(0),
69  fAdditionalTracks(0),
70  fHardCutoff(0),
71  fminpTTrig(20.),
72  fmaxpTTrig(50.),
73  fangWindowRecoil(0.6),
74  fSemigoodCorrect(0),
75  fHolePos(0),
76  fHoleWidth(0),
77  fRandom(0),
78  fqhat(1),
79  fxlength(2),
80  fCentSelectOn(kTRUE),
81  fCentMin(0),
82  fCentMax(10),
83  fOneConstSelectOn(kFALSE),
84  fDerivSubtrOrder(0),
85  fPhiJetCorr6(0x0),
86  fPhiJetCorr7(0x0),
87  fEtaJetCorr6(0x0),
88  fEtaJetCorr7(0x0),
89  fPtJetCorr(0x0),
90  fPtJet(0x0),
91  fhpTjetpT(0x0),
92  fhPt(0x0),
93  fhPhi(0x0),
94  fHLundIterative(0x0),
95  fHLundIterativeInject(0x0),
96  fNbOfConstvspT(0x0),
97  fTreeObservableTagging(0x0),
98  fTf1SoftOmega(0x0),
99  fTf1SoftKt(0x0),
100  fTf1Omega(0x0),
101  fTf1Kt(0x0),
102  fScaleELoss(kFALSE),
103  xfraction(1),
104  fAddMedScat(kFALSE),
105  fAddMedScatPtFrac(1),
106  fAddMedScatN(100)
107 
108 
109 {
110  for(Int_t i=0;i<13;i++){
111  fShapesVar[i]=0;}
112  SetMakeGeneralHistograms(kTRUE);
113  DefineOutput(1, TList::Class());
114  DefineOutput(2, TTree::Class());
115 }
116 
117 //________________________________________________________________________
119  AliAnalysisTaskEmcalJet(name, kTRUE),
120  fContainer(0),
121  fMinFractionShared(0),
122  fJetShapeType(kGenShapes),
123  fJetShapeSub(kNoSub),
124  fJetSelection(kInclusive),
125  fPtThreshold(-9999.),
126  fRMatching(0.2),
127  fSelectedShapes(0),
128  fSwitchKtNSub(0),
129  fSwitchMinNSub(0),
130  fSwitchAktNSub(0),
131  fSwitchSDKtNSub(0),
132  fSwitchSDMinNSub(0),
133  fAdditionalTracks(0),
134  fHardCutoff(0),
135  fminpTTrig(20.),
136  fmaxpTTrig(50.),
137  fangWindowRecoil(0.6),
138  fSemigoodCorrect(0),
139  fHolePos(0),
140  fHoleWidth(0),
141  fRandom(0),
142  fqhat(1),
143  fxlength(2),
144  fCentSelectOn(kTRUE),
145  fCentMin(0),
146  fCentMax(10),
147  fOneConstSelectOn(kFALSE),
148  fDerivSubtrOrder(0),
149  fPhiJetCorr6(0x0),
150  fPhiJetCorr7(0x0),
151  fEtaJetCorr6(0x0),
152  fEtaJetCorr7(0x0),
153  fPtJetCorr(0x0),
154  fPtJet(0x0),
155  fhpTjetpT(0x0),
156  fhPt(0x0),
157  fhPhi(0x0),
158  fHLundIterative(0x0),
159  fHLundIterativeInject(0x0),
160  fNbOfConstvspT(0x0),
161  fTreeObservableTagging(0x0),
162  fTf1SoftOmega(0x0),
163  fTf1SoftKt(0x0),
164  fTf1Omega(0x0),
165  fTf1Kt(0x0),
166  fScaleELoss(kFALSE),
167  xfraction(1),
168  fAddMedScat(kFALSE),
169  fAddMedScatPtFrac(1),
170  fAddMedScatN(100)
171 {
172  // Standard constructor.
173 
174 
175  for(Int_t i=0;i<13;i++){
176  fShapesVar[i]=0;}
177 
179 
180  DefineOutput(1, TList::Class());
181  DefineOutput(2, TTree::Class());
182 
183 
184 }
185 
186 //________________________________________________________________________
188 {
190  delete fTreeObservableTagging;
192  }
193 
194  if(fRandom) delete fRandom;
195  if(fTf1SoftOmega) delete fTf1SoftOmega;
196  if(fTf1SoftKt) delete fTf1SoftKt;
197  if(fTf1Omega) delete fTf1Omega;
198  if(fTf1Kt) delete fTf1Kt;
199 
200 }
201 
202 //________________________________________________________________________
204 {
205  // Create user output.
207 
208  Bool_t oldStatus = TH1::AddDirectoryStatus();
209  TH1::AddDirectory(oldStatus);
210 
211  //fTreeObservableTagging = new TTree("fTreeJetShape", "fTreeJetShape");
212 
213  //TH1::AddDirectory(oldStatus);
214 
215  const char* nameoutput = GetOutputSlot(2)->GetContainer()->GetName();
216  fTreeObservableTagging = new TTree(nameoutput, nameoutput);
217 
218  const Int_t nVar = 13;
219 
220  TString *fShapesVarNames = new TString [nVar];
221 
222  fShapesVarNames[0] = "partonCode";
223  fShapesVarNames[1] = "ptJet";
224  fShapesVarNames[2] = "ptDJet";
225  fShapesVarNames[3] = "mJet";
226  fShapesVarNames[4] = "nbOfConst";
227  fShapesVarNames[5] = "angularity";
228  fShapesVarNames[6] = "nitersd";
229  fShapesVarNames[7] = "niterall";
230  fShapesVarNames[8] = "weightPythia";
231  //fShapesVarNames[6] = "Nsubjet1kt";
232  // fShapesVarNames[7] = "Nsubjet2kt";
233  // fShapesVarNames[8] = "Nsubjet1Min";
234  // fShapesVarNames[9] = "Nsubjet2Min";
235  // fShapesVarNames[10] = "DeltaRkt";
236  // fShapesVarNames[11] = "DeltaRMin";
237  fShapesVarNames[9] = "SDSymm";
238  fShapesVarNames[10] = "scaledptJet";
239  fShapesVarNames[11] = "zg";
240  fShapesVarNames[12] = "rg";
241  // fShapesVarNames[15] = "SDGroomedN";
242  // fShapesVarNames[16] = "SDMass";
243  // fShapesVarNames[17] = "SDSymmkt";
244  // fShapesVarNames[18] = "SDDeltaRkt";
245  // fShapesVarNames[19] = "SDGroomedFrackt";
246  // fShapesVarNames[20] = "SDGroomedNkt";
247  // fShapesVarNames[21] = "SDMasskt";
248  // fShapesVarNames[22] = "SDSymmAkt";
249  // fShapesVarNames[23] = "SDDeltaRAkt";
250  // fShapesVarNames[24] = "SDGroomedFracAkt";
251  // fShapesVarNames[25] = "SDGroomedNAkt";
252  // fShapesVarNames[26] = "SDMassAkt";
253  // fShapesVarNames[27] = "SDSymmktForm";
254  // fShapesVarNames[28] = "SDDeltaRktForm";
255  // fShapesVarNames[29] = "SDGroomedFracktForm";
256  // fShapesVarNames[30] = "SDGroomedNktForm";
257  // fShapesVarNames[31] = "SDMassktForm";
258  // fShapesVarNames[32] = "SDSymmDemo";
259  // fShapesVarNames[33] = "SDDeltaRDemo";
260  // fShapesVarNames[34] = "SDGroomedFracDemo";
261  // fShapesVarNames[35] = "SDGroomedNDemo";
262  // fShapesVarNames[36] = "SDMassDemo";
263  // fShapesVarNames[42] = "SDSymmForm";
264  // fShapesVarNames[43] = "SDDeltaRForm";
265  // fShapesVarNames[44] = "SDGroomedFracForm";
266  // fShapesVarNames[45] = "SDGroomedNForm";
267  // fShapesVarNames[46] = "SDMassForm";
268  // fShapesVarNames[47] = "weightPythia";
269  // fShapesVarNames[42] = "SDSymmNoCut";
270  // fShapesVarNames[43] = "SDDeltaRNoCut";
271  // fShapesVarNames[44] = "SDGroomedFracNoCut";
272  // fShapesVarNames[45] = "SDGroomedNNoCut";
273  // fShapesVarNames[46] = "SDMassNoCut";
274 
275 
276  //fShapesVarNames[7] = "lesub";
277  //fShapesVarNames[8] = "CoreFraction";
278  //fShapesVarNames[9] = "Nsubjet1";
279  //fShapesVarNames[10] = "Nsubjet2";
280  //fShapesVarNames[11] = "DeltaR";
281  //fShapesVarNames[12] = "OpenAngle";
282  //fShapesVarNames[13] = "weightPythia";
283 
284  //fShapesVarNames[14] = "NT70";
285  //fShapesVarNames[15] = "nConstNT70";
286  //fShapesVarNames[16] = "NT80";
287  //fShapesVarNames[17] = "nConstNT80";
288  //fShapesVarNames[18] = "NT90";
289  //fShapesVarNames[19] = "nConstNT90";
290  //fShapesVarNames[20] = "NT95";
291  //fShapesVarNames[21] = "nConstNT95";
292 
293  //fShapesVarNames[22] = "SubjetFraction";
294 
295 
296  for(Int_t ivar=0; ivar < nVar; ivar++){
297  cout<<"looping over variables"<<endl;
298  fTreeObservableTagging->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/F", fShapesVarNames[ivar].Data()));}
299 
300 
301 
302  fPhiJetCorr6= new TH2F("fPhiJetCorr6", "fPhiJetCorr6", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
303  fOutput->Add(fPhiJetCorr6);
304  fEtaJetCorr6= new TH2F("fEtaJetCorr6", "fEtaJetCorr6", 50, -1.5, 1.5, 50, -1.5, 1.5);
305  fOutput->Add(fEtaJetCorr6);
306 
307  fPhiJetCorr7= new TH2F("fPhiJetCorr7", "fPhiJetCorr7", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
308  fOutput->Add(fPhiJetCorr7);
309  fEtaJetCorr7= new TH2F("fEtaJetCorr7", "fEtaJetCorr7", 50, -1.5, 1.5, 50, -1.5, 1.5);
310  fOutput->Add(fEtaJetCorr7);
311 
312  fPtJetCorr= new TH2F("fPtJetCorr", "fPtJetCorr", 100, 0, 200, 100, 0, 200);
313  fOutput->Add(fPtJetCorr);
314  fPtJet= new TH1F("fPtJet", "fPtJet", 100, 0, 200);
315  fOutput->Add(fPtJet);
316 
317  fhpTjetpT= new TH2F("fhpTjetpT", "fhpTjetpT", 200, 0, 200, 200, 0, 200);
318  fOutput->Add(fhpTjetpT);
319  fhPt= new TH1F("fhPt", "fhPt", 200, 0, 200);
320  fOutput->Add(fhPt);
321  fhPhi= new TH1F("fhPhi", "fhPhi", 100, -TMath::Pi(), TMath::Pi());
322  fOutput->Add(fhPhi);
323 
324 
325 
326  //log(1/theta),log(z*theta),jetpT,algo//
327  const Int_t dimSpec = 6;
328  const Int_t nBinsSpec[6] = {50,50,10,3,22,10};
329  const Double_t lowBinSpec[6] = {0.0,-10, 0,0,0,0};
330  const Double_t hiBinSpec[6] = {5.0, 0,200,3,22,10};
331  fHLundIterative = new THnSparseF("fHLundIterative",
332  "LundIterativePlot [log(1/theta),log(z*theta),pTjet,algo,partonFlavor,depth]",
333  dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
334  fOutput->Add(fHLundIterative);
335 
336  if(fAdditionalTracks>0){
337  //log(1/theta),log(z*theta),jetpT of added tracks
338  const Int_t dimSpecb = 4;
339  const Int_t nBinsSpecb[4] = {50,50,20,20};
340  const Double_t lowBinSpecb[4] = {0.0,-10, 0,0};
341  const Double_t hiBinSpecb[4] = {5.0, 0,200,200};
342  fHLundIterativeInject = new THnSparseF("fHLundIterativeInject",
343  "LundIterativePlotInject [log(1/theta),log(z*theta),pTjet,algo]",
344  dimSpecb,nBinsSpecb,lowBinSpecb,hiBinSpecb);
346 
347 
348  fNbOfConstvspT=new TH2F("fNbOfConstvspT", "fNbOfConstvspT", 100, 0, 100, 200, 0, 200);
349  fOutput->Add(fNbOfConstvspT);
350 
351  //fOutput->Add(fTreeObservableTagging);
352 
353  fRandom = new TRandom3(0);
354  PostData(1, fOutput); // Post data for ALL output slots > 0 here
355  PostData(2, fTreeObservableTagging);
356 
357  delete [] fShapesVarNames;
358 
359  }
360 
361 //________________________________________________________________________
363 {
364  // Run analysis code here, if needed. It will be executed before FillHistograms().
365  // if (gRandom) delete gRandom;
366  // gRandom = new TRandom3(0);
367  return kTRUE;
368 }
369 
370 //________________________________________________________________________
372 {
373  // Fill histograms.
374  //cout<<"IntoFillHistograms"<<endl;
375  AliEmcalJet* jet1 = NULL;
376  AliJetContainer *jetCont = GetJetContainer(0);
377 
378  Float_t kWeight=1;
379  if (fCentSelectOn)
380  if ((fCent>fCentMax) || (fCent<fCentMin)) return 0;
381 
382  AliAODTrack *triggerHadron = 0x0;
383 
384  if (fJetSelection == kRecoil) {
385  //Printf("Recoil jets!!!, fminpTTrig = %f, fmaxpTTrig = %f", fminpTTrig, fmaxpTTrig);
386  Int_t triggerHadronLabel = SelectTrigger(fminpTTrig, fmaxpTTrig);
387 
388 
389  if (triggerHadronLabel==-99999) {
390  //Printf ("Trigger Hadron not found, return");
391  return 0;}
392 
393 
395  TClonesArray *trackArrayAn = partContAn->GetArray();
396  triggerHadron = static_cast<AliAODTrack*>(trackArrayAn->At(triggerHadronLabel));
397 
398  if (!triggerHadron) {
399  //Printf("No Trigger hadron with the found label!!");
400  return 0;
401  }
402 
403  if(fSemigoodCorrect){
404  Double_t disthole=RelativePhi(triggerHadron->Phi(),fHolePos);
405  if(TMath::Abs(disthole)+fHoleWidth>TMath::Pi()-fangWindowRecoil){
406  return 0;}
407  }
408 
409  fhPt->Fill(triggerHadron->Pt());
410 
411  }
412 
413  if(jetCont) {
414  jetCont->ResetCurrentID();
415  while((jet1 = jetCont->GetNextAcceptJet())) {
416  //Printf("jet1=%p", jet1);
417  if (!jet1) continue;
418 
419  fPtJet->Fill(jet1->Pt());
420 
421 
422 
423 
425  Double_t disthole=RelativePhi(jet1->Phi(),fHolePos);
426  if(TMath::Abs(disthole)<fHoleWidth){
427  continue;
428  }
429  }
430 
431  Float_t dphiRecoil = 0.;
432  if (fJetSelection == kRecoil){
433  dphiRecoil = RelativePhi(triggerHadron->Phi(), jet1->Phi());
434  if (TMath::Abs(dphiRecoil) < (TMath::Pi() - fangWindowRecoil)) {
435  // Printf("Recoil jets back to back not found! continuing");
436  continue;
437  }
438 
439  fhpTjetpT->Fill(triggerHadron->Pt(), jet1->Pt());
440  //Printf(" ************ FILLING HISTOS****** shapeSub = %d, triggerHadron = %f, jet1 = %f", fJetShapeSub, triggerHadron->Pt(), jet1->Pt());
441  fhPhi->Fill(RelativePhi(triggerHadron->Phi(), jet1->Phi()));
442 
443  }
444 
445 
446  fShapesVar[0] = 0.;
447 
448  if (fJetShapeType == kGenShapes){
449  const AliEmcalPythiaInfo *partonsInfo = 0x0;
450  partonsInfo = GetPythiaInfo();
451  //Printf("partonsInfo=%p", partonsInfo);
452  Double_t jp1=RelativePhi(jet1->Phi(),partonsInfo->GetPartonPhi6());
453  Double_t detap1=(jet1->Eta())-(partonsInfo->GetPartonEta6());
454  kWeight=partonsInfo->GetPythiaEventWeight();
455  //Printf("kWeight=%f", kWeight);
456  fShapesVar[8] = kWeight;
457 
458  Float_t dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
459  fEtaJetCorr6->Fill(jet1->Eta(), partonsInfo->GetPartonEta6());
460  fPhiJetCorr6->Fill(jet1->Phi(), partonsInfo->GetPartonPhi6());
461  if(dRp1 < fRMatching) {
462  fShapesVar[0] = partonsInfo->GetPartonFlag6();
463  fPtJetCorr ->Fill(partonsInfo->GetPartonPt6(), jet1->Pt());
464  }
465  else {
466  jp1=RelativePhi(jet1->Phi(),partonsInfo->GetPartonPhi7());
467  detap1=(jet1->Eta())-(partonsInfo->GetPartonEta7());
468  dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
469  fEtaJetCorr7->Fill(jet1->Eta(), partonsInfo->GetPartonEta7());
470  fPhiJetCorr7->Fill(jet1->Phi(), partonsInfo->GetPartonPhi7());
471  if(dRp1 < fRMatching) {
472  fShapesVar[0] = partonsInfo->GetPartonFlag7();
473  fPtJetCorr->Fill(partonsInfo->GetPartonPt7(), jet1->Pt());
474  }
475  else fShapesVar[0]=0;
476  }
477  }
478 
479  Double_t ptSubtracted = 0;
480  ptSubtracted= jet1->Pt();
481  //Printf("ptSubtracted=%f", ptSubtracted);
482 
483 
484  if (ptSubtracted < fPtThreshold) continue;
485 
486  if ((fCentSelectOn == kFALSE) && (jet1->GetNumberOfTracks() < 2)) continue;
487 
488  //AliEmcalJetFinder *Reclusterer1; //Object containg Subjets from Subtracted Hybrid Jets
489  //Reclusterer1 = Recluster(jet1, 0, fJetRadius, fSubjetRadius, 1, 0, "SubJetFinder_1");
490 
491 
492 
493 
494 
495 
496 
497  fShapesVar[1] = ptSubtracted;
498  fShapesVar[2] = GetJetpTD(jet1,0);
499  fShapesVar[3] = GetJetMass(jet1,0);
500  fShapesVar[4] = 1.*GetJetNumberOfConstituents(jet1,0);
501  fShapesVar[5] = GetJetAngularity(jet1,0);
502  SoftDrop(jet1,jetCont,0.1,0,0);
503 
504  RecursiveParents(jet1,jetCont,0,fShapesVar[0]);
505  RecursiveParents(jet1,jetCont,1,fShapesVar[0]);
506  RecursiveParents(jet1,jetCont,2,fShapesVar[0]);
507 
508 
509  fTreeObservableTagging->Fill();
510 
511  }
512 
513  }
514 
515  return kTRUE;
516 }
517 
518 //________________________________________________________________________
520  //calc subtracted jet mass
521  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
523  else return jet->GetShapeProperties()->GetSecondOrderSubtracted();
524  else
525  return jet->M();
526 }
527 
528 //________________________________________________________________________
530 
531  AliJetContainer *jetCont = GetJetContainer(jetContNb);
532  if (!jet->GetNumberOfTracks())
533  return 0;
534  Double_t den=0.;
535  Double_t num = 0.;
536  AliVParticle *vp1 = 0x0;
537  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
538  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
539 
540  if (!vp1){
541  Printf("AliVParticle associated to constituent not found");
542  continue;
543  }
544 
545  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
546  Double_t dr2 = (vp1->Eta()-jet->Eta())*(vp1->Eta()-jet->Eta()) + dphi*dphi;
547  Double_t dr = TMath::Sqrt(dr2);
548  num=num+vp1->Pt()*dr;
549  den=den+vp1->Pt();
550  }
551  return num/den;
552 }
553 
554 //________________________________________________________________________
556 
557  if((fJetShapeSub==kDerivSub) && (jetContNb==0))
560  else
561  return Angularity(jet, jetContNb);
562 
563 }
564 
565 
566 //________________________________________________________________________
568 
569  AliJetContainer *jetCont = GetJetContainer(jetContNb);
570  if (!jet->GetNumberOfTracks())
571  return 0;
572  Double_t den=0.;
573  Double_t num = 0.;
574  AliVParticle *vp1 = 0x0;
575  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
576  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
577 
578  if (!vp1){
579  Printf("AliVParticle associated to constituent not found");
580  continue;
581  }
582 
583  num=num+vp1->Pt()*vp1->Pt();
584  den=den+vp1->Pt();
585  }
586  return TMath::Sqrt(num)/den;
587 }
588 
589 //________________________________________________________________________
591  //calc subtracted jet mass
592  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
594  else return jet->GetShapeProperties()->GetSecondOrderSubtractedpTD();
595  else
596  return PTD(jet, jetContNb);
597 
598 }
599 
600 //_____________________________________________________________________________
602 
603  AliJetContainer *jetCont = GetJetContainer(jetContNb);
604  if (!jet->GetNumberOfTracks())
605  return 0;
606  Double_t mxx = 0.;
607  Double_t myy = 0.;
608  Double_t mxy = 0.;
609  int nc = 0;
610  Double_t sump2 = 0.;
611  Double_t pxjet=jet->Px();
612  Double_t pyjet=jet->Py();
613  Double_t pzjet=jet->Pz();
614 
615 
616  //2 general normalized vectors perpendicular to the jet
617  TVector3 ppJ1(pxjet, pyjet, pzjet);
618  TVector3 ppJ3(- pxjet* pzjet, - pyjet * pzjet, pxjet * pxjet + pyjet * pyjet);
619  ppJ3.SetMag(1.);
620  TVector3 ppJ2(-pyjet, pxjet, 0);
621  ppJ2.SetMag(1.);
622  AliVParticle *vp1 = 0x0;
623  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
624  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
625 
626  if (!vp1){
627  Printf("AliVParticle associated to constituent not found");
628  continue;
629  }
630 
631  TVector3 pp(vp1->Px(), vp1->Py(), vp1->Pz());
632 
633  //local frame
634  TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
635  TVector3 pPerp = pp - pLong;
636  //projection onto the two perpendicular vectors defined above
637 
638  Float_t ppjX = pPerp.Dot(ppJ2);
639  Float_t ppjY = pPerp.Dot(ppJ3);
640  Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
641  if(ppjT<=0) return 0;
642 
643  mxx += (ppjX * ppjX / ppjT);
644  myy += (ppjY * ppjY / ppjT);
645  mxy += (ppjX * ppjY / ppjT);
646  nc++;
647  sump2 += ppjT;}
648 
649  if(nc<2) return 0;
650  if(sump2==0) return 0;
651  // Sphericity Matrix
652  Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
653  TMatrixDSym m0(2,ele);
654 
655  // Find eigenvectors
656  TMatrixDSymEigen m(m0);
657  TVectorD eval(2);
658  TMatrixD evecm = m.GetEigenVectors();
659  eval = m.GetEigenValues();
660  // Largest eigenvector
661  int jev = 0;
662  // cout<<eval[0]<<" "<<eval[1]<<endl;
663  if (eval[0] < eval[1]) jev = 1;
664  TVectorD evec0(2);
665  // Principle axis
666  evec0 = TMatrixDColumn(evecm, jev);
667  Double_t compx=evec0[0];
668  Double_t compy=evec0[1];
669  TVector2 evec(compx, compy);
670  Double_t circ=0;
671  if(jev==1) circ=2*eval[0];
672  if(jev==0) circ=2*eval[1];
673 
674  return circ;
675 
676 
677 
678 }
679 
680 
681 
682 
683 //________________________________________________________________________
685  //calc subtracted jet mass
686 
687  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
690  else
691  return Circularity(jet, jetContNb);
692 
693 }
694 
695 //________________________________________________________________________
697 
698  AliJetContainer *jetCont = GetJetContainer(jetContNb);
699  if (!jet->GetNumberOfTracks())
700  return 0;
701  Double_t den=0.;
702  Double_t num = 0.;
703  AliVParticle *vp1 = 0x0;
704  AliVParticle *vp2 = 0x0;
705  std::vector<int> ordindex;
706  ordindex=jet->GetPtSortedTrackConstituentIndexes(jetCont->GetParticleContainer()->GetArray());
707  //Printf("Nbof const = %d", jet->GetNumberOfTracks());
708  //Printf("ordindex[0] = %d, ordindex[1] = %d", ordindex[0], ordindex[1]);
709 
710  if(ordindex.size()<2) return -1;
711 
712  vp1 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[0], jetCont->GetParticleContainer()->GetArray()));
713  if (!vp1){
714  Printf("AliVParticle associated to Leading constituent not found");
715  return -1;
716  }
717 
718  vp2 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[1], jetCont->GetParticleContainer()->GetArray()));
719  if (!vp2){
720  Printf("AliVParticle associated to Subleading constituent not found");
721  return -1;
722  }
723 
724 
725  num=vp1->Pt();
726  den=vp2->Pt();
727  //Printf("vp1->Pt() =%f, vp2->Pt() =%f", vp1->Pt(), vp2->Pt());
728 
729 return num-den;
730 }
731 
732 //________________________________________________________________________
734  //calc subtracted jet mass
735 
736  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
739  else
740  return LeSub(jet, jetContNb);
741 
742 }
743 
744 //________________________________________________________________________
746  //calc subtracted jet mass
747 
748  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
751  else
752  return jet->GetNumberOfTracks();
753 
754  }
755 
756 
757 //________________________________________________________________________
758 AliEmcalJetFinder *AliAnalysisTaskEmcalJetShapesMC::Recluster(AliEmcalJet *Jet, Int_t JetContNb, Double_t JetRadius, Double_t SubJetRadius, Double_t SubJetMinPt, Int_t Algorithm, const char* Name){
759 
760  AliJetContainer *JetCont = GetJetContainer(JetContNb);
761  AliEmcalJetFinder *Reclusterer = new AliEmcalJetFinder(Name); //JetFinder Object for reclustered jets
762  Reclusterer->SetRadius(SubJetRadius);
763  Reclusterer->SetJetMinPt(SubJetMinPt);
764  Reclusterer->SetJetAlgorithm(Algorithm); //0 for anti-kt 1 for kt
765  Reclusterer->SetJetMaxEta(0.9-JetRadius);
766  Reclusterer->SetRecombSheme(0);
767 
768 
769  //Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
770  Double_t dVtx[3]={0.,0.,0.};
771  if(Reclusterer->AliEmcalJetFinder::Filter(Jet, JetCont, dVtx)){;} //reclustering jet1 using the jetfinderobject Reclusterer
772  return Reclusterer;
773 }
774 
775 
776 
777 
778 //----------------------------------------------------------------------
780  AliJetContainer *JetCont = GetJetContainer(JetContNb);
781 
782 
783 
784  Double_t SubJetiness_Numerator = 0;
785  Double_t SubJetiness_Denominator = 0;
786  Double_t Index=-2;
787  if (Reclusterer->GetNumberOfJets() < 1) return -2;
788  Index=SubJetOrdering(Jet,Reclusterer,1,0,kTRUE);
789  if(Index==-999) return -2;
790  SubJetiness_Numerator=(Reclusterer->GetJet(Index)->Pt());
791  SubJetiness_Denominator=Jet->Pt();
792  return SubJetiness_Numerator/SubJetiness_Denominator;
793 
794 
795 }
796 //__________________________________________________________________________________
798 
799  AliJetContainer *jetCont = GetJetContainer(jetContNb);
800  if (!jet->GetNumberOfTracks())
801  return 0;
802  Double_t den=0.;
803  Double_t num = 0.;
804  AliVParticle *vp1 = 0x0;
805  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
806  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
807 
808  if (!vp1){
809  Printf("AliVParticle associated to constituent not found");
810  continue;
811  }
812 
813  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
814  Double_t dr2 = (vp1->Eta()-jet->Eta())*(vp1->Eta()-jet->Eta()) + dphi*dphi;
815  Double_t dr = TMath::Sqrt(dr2);
816  if(dr<=fSubjetRadius) num=num+vp1->Pt();
817 
818  }
819  return num/jet->Pt();
820 }
821 
822 
823 
824 
825 //________________________________________________________________________
827  //calc subtracted jet mass
828 
829  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
832  else
833  return CoreFrac(jet, jetContNb);
834 
835 }
836 
837 
838 
839 
840 //----------------------------------------------------------------------
842  AliEmcalJet *SubJet=NULL;
843  Double_t SortingVariable;
844  Int_t ArraySize =N+1;
845  TArrayD *JetSorter = new TArrayD(ArraySize);
846  TArrayD *JetIndexSorter = new TArrayD(ArraySize);
847  for (Int_t i=0; i<ArraySize; i++){
848  JetSorter->SetAt(0,i);
849  }
850  for (Int_t i=0; i<ArraySize; i++){
851  JetIndexSorter->SetAt(0,i);
852  }
853  if(Reclusterer->GetNumberOfJets()<N) return -999;
854  for (Int_t i=0; i<Reclusterer->GetNumberOfJets(); i++){
855  SubJet=Reclusterer->GetJet(i);
856  if (Type==0) SortingVariable=SubJet->Pt();
857  else if (Type==1) SortingVariable=SubJet->E();
858  else if (Type==2) SortingVariable=SubJet->M();
859  for (Int_t j=0; j<N; j++){
860  if (SortingVariable>JetSorter->GetAt(j)){
861  for (Int_t k=N-1; k>=j; k--){
862  JetSorter->SetAt(JetSorter->GetAt(k),k+1);
863  JetIndexSorter->SetAt(JetIndexSorter->GetAt(k),k+1);
864  }
865  JetSorter->SetAt(SortingVariable,j);
866  JetIndexSorter->SetAt(i,j);
867  break;
868  }
869  }
870  }
871  if (!Index) return JetSorter->GetAt(N-1);
872  else return JetIndexSorter->GetAt(N-1);
873 }
874 
875 
876 
877 //returns -1 if the Nth hardest jet is requested where N>number of available jets
878 //type: 0=Pt 1=E 2=M
879 //Index TRUE=returns index FALSE=returns value of quantatiy in question
880 
881 
882 
883 
884 
885 //______________________________________________________________________________
887 
888  AliJetContainer *jetCont = GetJetContainer(jetContNb);
889  if (!jet->GetNumberOfTracks())
890  return 0;
891  Double_t mxx = 0.;
892  Double_t myy = 0.;
893  Double_t mxy = 0.;
894  int nc = 0;
895  Double_t sump2 = 0.;
896 
897  AliVParticle *vp1 = 0x0;
898  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
899  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
900 
901  if (!vp1){
902  Printf("AliVParticle associated to constituent not found");
903  continue;
904  }
905 
906  Double_t ppt=vp1->Pt();
907  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
908 
909  Double_t deta = vp1->Eta()-jet->Eta();
910  mxx += ppt*ppt*deta*deta;
911  myy += ppt*ppt*dphi*dphi;
912  mxy -= ppt*ppt*deta*TMath::Abs(dphi);
913  nc++;
914  sump2 += ppt*ppt;
915 
916  }
917  if(nc<2) return 0;
918  if(sump2==0) return 0;
919  // Sphericity Matrix
920  Double_t ele[4] = {mxx , mxy , mxy , myy };
921  TMatrixDSym m0(2,ele);
922 
923  // Find eigenvectors
924  TMatrixDSymEigen m(m0);
925  TVectorD eval(2);
926  TMatrixD evecm = m.GetEigenVectors();
927  eval = m.GetEigenValues();
928  // Largest eigenvector
929  int jev = 0;
930  // cout<<eval[0]<<" "<<eval[1]<<endl;
931  if (eval[0] < eval[1]) jev = 1;
932  TVectorD evec0(2);
933  // Principle axis
934  evec0 = TMatrixDColumn(evecm, jev);
935  Double_t compx=evec0[0];
936  Double_t compy=evec0[1];
937  TVector2 evec(compx, compy);
938  Double_t sig=0;
939  if(jev==1) sig=TMath::Sqrt(TMath::Abs(eval[0])/sump2);
940  if(jev==0) sig=TMath::Sqrt(TMath::Abs(eval[1])/sump2);
941 
942  return sig;
943 
944 }
945 
946 //________________________________________________________________________
948  //calc subtracted jet mass
949 
950  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
953  else
954  return Sigma2(jet, jetContNb);
955 
956 }
957 
958 
959 //_________________________________________________________________________
961 
962  AliJetContainer *jetCont = GetJetContainer(jetContNb);
963  if (!jet->GetNumberOfTracks())
964  return;
965 
966  Double_t ptJet = jet->Pt();
967 
968  AliVParticle *vp1 = 0x0;
969  std::vector<int> ordindex;
970  ordindex=jet->GetPtSortedTrackConstituentIndexes(jetCont->GetParticleContainer()->GetArray());
971  //Printf("Nbof const = %d", jet->GetNumberOfTracks());
972  //Printf("ordindex[0] = %d, ordindex[1] = %d", ordindex[0], ordindex[1]);
973  //if(ordindex.size()<2) return -1;
974 
975  for(Int_t iconst =0; iconst<jet->GetNumberOfTracks(); iconst++){
976 
977  vp1 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[iconst], jetCont->GetParticleContainer()->GetArray()));
978  if (!vp1){
979  Printf("AliVParticle associated to Leading constituent not found");
980  return;
981  }
982 
983  if (nTFractions[0] <= 0.7*ptJet){
984  nTFractions[0] += vp1->Pt();
985  nTFractions[1] +=1;
986  }
987 
988  if (nTFractions[2] <= 0.8*ptJet){
989  nTFractions[2] += vp1->Pt();
990  nTFractions[3] +=1;
991  }
992 
993  if (nTFractions[4] <= 0.9*ptJet){
994  nTFractions[4] += vp1->Pt();
995  nTFractions[5] +=1;
996  }
997 
998  if (nTFractions[6] <= 0.95*ptJet){
999  nTFractions[6] += vp1->Pt();
1000  nTFractions[7] +=1;
1001  }
1002  }
1003 }
1004 //_________________________________________________________________________________________________
1006 
1007  //WARNING!!! Only works for parent jets that are clustered with Anti-Kt! To change go to AliEmcalJetFinder.cxx and look at the Nsubjettiness() function
1008 
1009  //Algorithm==0 -> kt_axes;
1010  // Algorithm==1 -> ca_axes;
1011  //Algorithm==2 -> antikt_0p2_axes;
1012  //Algorithm==3 -> wta_kt_axes;
1013  //Algorithm==4 -> wta_ca_axes;
1014  //Algorithm==5 -> onepass_kt_axes;
1015  //Algorithm==6 -> onepass_ca_axes;
1016  //Algorithm==7 -> onepass_antikt_0p2_axes;
1017  //Algorithm==8 -> onepass_wta_kt_axes;
1018  //Algorithm==9 -> onepass_wta_ca_axes;
1019  //Algorithm==10 -> min_axes;
1020 
1021 
1022  //Option==0 returns Nsubjettiness Value
1023  //Option==1 && N==2 returns opening angle between two subjet axes(Delta R?)
1024  //Option==2 && N==2 returns Delta R
1025 
1026  if (Jet->GetNumberOfTracks()>=N){
1027  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1028  AliEmcalJetFinder *JetFinder=new AliEmcalJetFinder("Nsubjettiness");
1029  JetFinder->SetJetMaxEta(0.9-fJetRadius);
1030  JetFinder->SetRadius(fJetRadius);
1031  JetFinder->SetJetAlgorithm(0); //0 for anti-kt 1 for kt //this is for the JET!!!!!!!!!! Not the SubJets
1032  JetFinder->SetRecombSheme(0);
1033  JetFinder->SetJetMinPt(Jet->Pt());
1034  //Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
1035  Double_t dVtx[3]={1,1,1};
1036  //Printf("JetFinder->Nsubjettiness =%f", JetFinder->Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,fSubjetRadius,Beta,Option));
1037  return JetFinder->Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,0.2,Beta,Option,0,Beta_SD,ZCut,SoftDropOn);
1038 
1039  }
1040  else return -2;
1041 }
1042 
1043 
1044 
1045 //________________________________________________________________________
1047 
1049  TClonesArray *tracksArray = partCont->GetArray();
1050 
1051  if(!partCont || !tracksArray) return -99999;
1052  AliAODTrack *track = 0x0;
1053  AliEmcalParticle *emcPart = 0x0;
1054 
1055 
1056  TList *trackList = new TList();
1057  Int_t triggers[100];
1058  for (Int_t iTrigger=0; iTrigger<100; iTrigger++) triggers[iTrigger] = 0;
1059  Int_t iTT = 0;
1060 
1061  for(Int_t iTrack=0; iTrack <= tracksArray->GetEntriesFast(); iTrack++){
1062 
1063 
1064  if (fJetShapeSub == kConstSub){
1065  emcPart = static_cast<AliEmcalParticle*>(tracksArray->At(iTrack));
1066  if (!emcPart) continue;
1067  if(TMath::Abs(emcPart->Eta())>0.9) continue;
1068  if (emcPart->Pt()<0.15) continue;
1069 
1070  if ((emcPart->Pt() >= minpT) && (emcPart->Pt()< maxpT)) {
1071  trackList->Add(emcPart);
1072  triggers[iTT] = iTrack;
1073  iTT++;
1074  }
1075  }
1076  else{
1077  track = static_cast<AliAODTrack*>(tracksArray->At(iTrack));
1078  if (!track) continue;
1079  if(TMath::Abs(track->Eta())>0.9) continue;
1080  if (track->Pt()<0.15) continue;
1081  if (!(track->TestFilterBit(768))) continue;
1082 
1083  if ((track->Pt() >= minpT) && (track->Pt()< maxpT)) {
1084  trackList->Add(track);
1085  triggers[iTT] = iTrack;
1086  iTT++;
1087 
1088  }
1089  }
1090  }
1091 
1092  if (iTT == 0) return -99999;
1093  Int_t nbRn = 0, index = 0 ;
1094  TRandom3* random = new TRandom3(0);
1095  nbRn = random->Integer(iTT);
1096 
1097  index = triggers[nbRn];
1098  //Printf("iTT Total= %d, nbRn = %d, Index = %d",iTT, nbRn, index );
1099  return index;
1100 
1101 }
1102 
1103 //__________________________________________________________________________________
1105 
1106  if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1107  else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1108  if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1109  else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1110  double dphi = mphi-vphi;
1111  if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1112  else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1113  return dphi;//dphi in [-Pi, Pi]
1114 }
1115 
1116 
1117 //________________________________________________________________________
1119  //
1120  // retrieve event objects
1121  //
1123  return kFALSE;
1124 
1125  return kTRUE;
1126 }
1127 
1128 //_______________________________________________________________________
1130 {
1131  // Called once at the end of the analysis.
1132 
1133  AliInfo("Terminate");
1134  AliAnalysisTaskSE::Terminate();
1135 
1136  fOutput = dynamic_cast<AliEmcalList*> (GetOutputData(1));
1137  if (!fOutput) {
1138  AliError("fOutput not available");
1139  return;
1140  }
1141 
1142  fTreeObservableTagging = dynamic_cast<TTree*>(GetOutputData(2));
1143  if (!fTreeObservableTagging){
1144  Printf("ERROR: fTreeObservableTagging not available");
1145  return;
1146  }
1147 
1148 }
1149 
1150 //_________________________________________________________________________
1152 
1153  std::vector<fastjet::PseudoJet> fInputVectors;
1154  fInputVectors.clear();
1155  Double_t JetInvMass=0, PseudJetInvMass=0, TrackMom = 0, TrackEnergy = 0;
1156  fastjet::PseudoJet PseudoTracks;
1157  fastjet::PseudoJet MyJet;
1158  fastjet::PseudoJet PseudoTracksCMS;
1159  AliParticleContainer *fTrackCont = fJetCont->GetParticleContainer();
1160  //cout<<"CALL TO SOFTDROP"<<endl;
1161 
1162  Double_t zeta=0;
1163  Double_t angle=0;
1164  if (fTrackCont) for (Int_t i=0; i<fJet->GetNumberOfTracks(); i++) {
1165  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
1166  if (!fTrk) continue;
1167  JetInvMass += fTrk->M();
1168 
1169  PseudoTracks.reset(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
1170  TrackMom += TMath::Sqrt(TMath::Power(fTrk->Px(),2)+TMath::Power(fTrk->Py(),2)+TMath::Power(fTrk->Pz(),2));
1171  TrackEnergy += fTrk->E();
1172  PseudoTracks.set_user_index(fJet->TrackAt(i)+100);
1173  PseudJetInvMass += PseudoTracks.m();
1174  fInputVectors.push_back(PseudoTracks);
1175 
1176  }
1177 
1178  //fRandom->SetSeed(0);
1179  //here add N tracks with random phi and eta and theta according to bdmps distrib.
1180 
1181  MyJet.reset(fJet->Px(),fJet->Py(),fJet->Pz(),fJet->E());
1182  Double_t omegac=0.5*fqhat*fxlength*fxlength/0.2;
1183  Double_t thetac=TMath::Sqrt(12*0.2/(fqhat*TMath::Power(fxlength,3)));
1184  Double_t xQs=TMath::Sqrt(fqhat*fxlength);
1185  //cout<<"medium parameters "<<omegac<<" "<<thetac<<" "<<xQs<<endl;
1186 
1187  for(Int_t i=0;i<fAdditionalTracks;i++){
1188 
1189  Double_t ppx,ppy,ppz,kTscale,lim2o,lim1o;
1190  Double_t lim2=xQs;
1191  Double_t lim1=10000;
1192 
1193  //generation of kT according to 1/kT^4, with minimum QS=2 GeV and maximum ~sqrt(ptjet*T)
1194  fTf1Kt= new TF1("fTf1Kt","1/(x*x*x*x)",lim2,lim1);
1195  kTscale=fTf1Kt->GetRandom();
1196  //generation within the jet cone
1197 
1198  //generation of w according to 1/w, with minimum wc
1199  //omega needs to be larger than kT so to have well defined angles
1200  lim2o=kTscale;
1201  lim1o=kTscale/TMath::Sin(0.1);
1202  fTf1Omega= new TF1("fTf1Omega","1/x",lim2o,lim1o);
1203  Double_t omega=fTf1Omega->GetRandom();
1204 
1205  Double_t sinpptheta=kTscale/omega;
1206  Double_t pptheta=TMath::ASin(sinpptheta);
1207  //cout<<"angle_omega_kt"<<pptheta<<" "<<omega<<" "<<kTscale<<endl;
1208  if(pptheta>fJetRadius) continue;
1209 
1210  PseudoTracksCMS.reset(kTscale/TMath::Sqrt(2),kTscale/TMath::Sqrt(2),omega*TMath::Cos(pptheta),omega);
1211  //boost the particle in the rest frame of the jet to the lab frame
1212  fastjet::PseudoJet PseudoTracksLab=PseudoTracksCMS.boost(MyJet);
1213  PseudoTracksLab.set_user_index(i+fJet->GetNumberOfTracks()+100);
1214  fInputVectors.push_back(PseudoTracksLab);
1215  //in the frame of the jet
1216  zeta=omega/fJet->E();
1217  angle=pptheta;
1218 
1219 
1220  }
1221 
1222 
1223 
1224 
1225 
1226 
1227 
1228 
1229 
1230  fastjet::JetDefinition fJetDef(fastjet::antikt_algorithm, fJetRadius*2, static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
1231  fastjet::contrib::Recluster *recluster;
1232  try {
1233  fastjet::ClusterSequence fClustSeqSA(fInputVectors, fJetDef);
1234  std::vector<fastjet::PseudoJet> fOutputJets;
1235  fOutputJets.clear();
1236  fOutputJets=fClustSeqSA.inclusive_jets(0);
1237 
1238  //cout<<fOutputJets[0].perp()<<" "<<fJet->Pt()<<endl;
1239 
1240  fastjet::contrib::SoftDrop softdrop(beta, zcut);
1241 
1242  softdrop.set_verbose_structure(kTRUE);
1243  fastjet::JetAlgorithm jetalgo(fastjet::cambridge_algorithm);
1244  if(ReclusterAlgo==2) jetalgo=fastjet::antikt_algorithm;
1245  if(ReclusterAlgo==1) jetalgo=fastjet::kt_algorithm;
1246  if(ReclusterAlgo==0) jetalgo=fastjet::cambridge_algorithm;
1247 
1248  recluster = new fastjet::contrib::Recluster(jetalgo,1,true);
1249  softdrop.set_reclustering(true,recluster);
1250  fastjet::PseudoJet finaljet = softdrop(fOutputJets[0]);
1251  Int_t NDroppedTracks = fJet->GetNumberOfTracks()-finaljet.constituents().size();
1252 
1253  Double_t SymParam, Mu, DeltaR, GroomedPt,GroomedMass;
1254  Int_t NGroomedBranches;
1255  SymParam=(finaljet.structure_of<fastjet::contrib::SoftDrop>().symmetry());
1256  Mu=(finaljet.structure_of<fastjet::contrib::SoftDrop>().mu());
1257  DeltaR=(finaljet.structure_of<fastjet::contrib::SoftDrop>().delta_R());
1258  NGroomedBranches=finaljet.structure_of<fastjet::contrib::SoftDrop>().dropped_count();
1259  GroomedPt=finaljet.perp();
1260  GroomedMass=finaljet.m();
1261 
1262 
1263  fShapesVar[11]=SymParam;
1264  fShapesVar[12]=DeltaR;
1265  // fShapesVar[14]=zeta;
1266  // fShapesVar[15]=angle;
1267  // fShapesVar[16]=GroomedMass;}
1268  // if(ReclusterAlgo==1){
1269  // fShapesVar[17]=SymParam;
1270  // fShapesVar[18]=DeltaR;
1271  // fShapesVar[19]=zeta;
1272  // fShapesVar[20]=angle;
1273  // fShapesVar[21]=GroomedMass; }
1274 
1275  // if(ReclusterAlgo==2){
1276  // fShapesVar[22]=SymParam;
1277  // fShapesVar[23]=DeltaR;
1278  // fShapesVar[24]=zeta;
1279  // fShapesVar[25]=angle;
1280  // fShapesVar[26]=GroomedMass;
1281  // }}
1282  // if(beta==1){
1283  // fShapesVar[27]=SymParam;
1284  // fShapesVar[28]=DeltaR;
1285  // fShapesVar[29]=zeta;
1286  // fShapesVar[30]=angle;
1287  // fShapesVar[31]=GroomedMass;
1288  // }
1289  // //this one kills soft and large angle radiation
1290  // if((beta==1.5) && (zcut==0.5)){
1291  // fShapesVar[32]=SymParam;
1292  // fShapesVar[33]=DeltaR;
1293  // fShapesVar[34]=zeta;
1294  // fShapesVar[35]=angle;
1295  // fShapesVar[36]=GroomedMass; }
1296  // //this option favour democratic branches at large kt
1297  // if((beta==-1) && (zcut==0.005)){
1298  // fShapesVar[37]=SymParam;
1299  // fShapesVar[38]=DeltaR;
1300  // fShapesVar[39]=zeta;
1301  // fShapesVar[40]=angle;
1302  // fShapesVar[41]=GroomedMass; }
1303 
1304  // if((beta==-2) && (zcut==0.005)){
1305  // fShapesVar[42]=SymParam;
1306  // fShapesVar[43]=DeltaR;
1307  // fShapesVar[44]=zeta;
1308  // fShapesVar[45]=angle;
1309  // fShapesVar[46]=GroomedMass; }
1310 
1311 
1312 
1313 
1314 
1315 
1316 
1317 
1318 
1319  } catch (fastjet::Error) {
1320  AliError(" [w] FJ Exception caught.");
1321  //return -1;
1322  }
1323 
1324 
1325 
1326 
1327  if(recluster) { delete recluster; recluster = NULL; }
1328 
1329  if(fTf1Kt){ delete fTf1Kt;}
1330  if(fTf1Omega){ delete fTf1Omega;}
1331 
1332 
1333 
1334  return;
1335 
1336 
1337 }
1338 
1339 
1340 //_________________________________________________________________________
1342 
1343  std::vector<fastjet::PseudoJet> fInputVectors;
1344  fInputVectors.clear();
1345  fastjet::PseudoJet PseudoTracks;
1346  fastjet::PseudoJet PseudoTracksLab;
1347  double xflagalgo=0;
1348  double lnpt_relinject=0;
1349  double yinject=0;
1350  int xflagAdded=0;
1351  double zinject,angleinject,pptheta,sinpptheta,omega,omega2,angle2;
1352  AliParticleContainer *fTrackCont = fJetCont->GetParticleContainer();
1353  Float_t pTscale=0., phiscale=0., thetascale=0., pXscale=0., pYscale=0., pZscale=0., pscale=0.;
1354 
1355  if (fTrackCont) for (Int_t i=0; i<fJet->GetNumberOfTracks(); i++) {
1356  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
1357  if (!fTrk) continue;
1358  if (fScaleELoss){
1359  pTscale = xfraction*sqrt(pow(fTrk->Px(),2)+pow(fTrk->Py(),2));
1360  phiscale = fTrk->Phi();
1361  thetascale = 2.*TMath::ATan(TMath::Exp(-1.*(fTrk->Eta())));
1362  pXscale = pTscale * TMath::Cos(phiscale);
1363  pYscale = pTscale * TMath::Sin(phiscale);
1364  pZscale = pTscale/TMath::Tan(thetascale);
1365  pscale = TMath::Sqrt(pTscale*pTscale+pZscale*pZscale);
1366  PseudoTracks.reset(pXscale, pYscale, pZscale, pscale);
1367  }
1368  else PseudoTracks.reset(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
1369  PseudoTracks.set_user_index(fJet->TrackAt(i)+100);
1370  fInputVectors.push_back(PseudoTracks);
1371 
1372  }
1373  /*if(fAddMedScat){
1374  for(int i = 0; i < fAddMedScatN; i++){
1375  TRandom3 rand1(0),rand2(0),rand3(0); //set range +- jet R
1376  Double_t randN1 = 0.4*0.4*rand1.Rndm();
1377  Double_t randN2 = 2*TMath::Pi()*rand2.Rndm();
1378  Double_t phi_rand = (fJet->Phi())+TMath::Sqrt(randN1)*TMath::Sin(randN2);
1379  Double_t eta_rand = (fJet->Eta())+TMath::Sqrt(randN1)*TMath::Cos(randN2);
1380  Double_t fAddMedScatPt = (fAddMedScatPtFrac*fJet->Pt())/fAddMedScatN;
1381  PseudoTracks.reset(fAddMedScatPt*TMath::Cos(phi_rand),fAddMedScatPt*TMath::Sin(phi_rand),fAddMedScatPt/TMath::Tan(eta_rand),fAddMedScatPt);
1382  PseudoTracks.set_user_index(i+fJet->GetNumberOfTracks()+100);
1383  fInputVectors.push_back(PseudoTracks);
1384  }
1385  }*/
1386  if(fAddMedScat){
1387  for(int i = 0; i < fAddMedScatN; i++){
1388  Double_t ppx,ppy,ppz,SoftkTscale,lim2o,lim1o;
1389  Double_t lim1=0.1;
1390  Double_t lim2=0.5;
1391  fTf1SoftKt= new TF1("fTf1SoftKt","1/(x)",lim1,lim2);
1392  SoftkTscale=fTf1SoftKt->GetRandom();
1393 
1394  lim2o=SoftkTscale;
1395  lim1o=SoftkTscale/TMath::Sin(0.1);
1396  fTf1SoftOmega= new TF1("fTf1SoftOmega","1/x",lim2o,lim1o);
1397  omega=fTf1SoftOmega->GetRandom();
1398  sinpptheta=SoftkTscale/omega;
1399  pptheta=TMath::ASin(sinpptheta);
1400  if(pptheta>fJetRadius) continue;
1401 
1402  TLorentzVector pTrackCMS(SoftkTscale/TMath::Sqrt(2),SoftkTscale/TMath::Sqrt(2),omega*TMath::Cos(pptheta),omega);
1403  TVector3 MyJet(fJet->Px(),fJet->Py(),fJet->Pz());
1404  TVector3 direction = MyJet.Unit();
1405  //rotate the track to the jet frame
1406  pTrackCMS.RotateUz(direction);
1407 
1408  //add the rotated track to the jet
1409  PseudoTracksLab.reset(pTrackCMS.Px(),pTrackCMS.Py(),pTrackCMS.Pz(),pTrackCMS.E());
1410 
1411  PseudoTracksLab.set_user_index(i+fJet->GetNumberOfTracks()+100);
1412 
1413  omega2=PseudoTracksLab.perp();
1414  angle2=pTrackCMS.Angle(MyJet);
1415 
1416  fInputVectors.push_back(PseudoTracksLab);
1417  }
1418  }
1419 
1420 
1421  //add tracks to the jet prior to the reclusterer in case of iterative mapping of splittings
1422 
1423  Double_t omegac=0.5*fqhat*fxlength*fxlength/0.2;
1424  Double_t thetac=TMath::Sqrt(12*0.2/(fqhat*TMath::Power(fxlength,3)));
1425  Double_t xQs=TMath::Sqrt(fqhat*fxlength);
1426 
1427 
1428  for(Int_t i=0;i<fAdditionalTracks;i++){
1429 
1430  Double_t ppx,ppy,ppz,kTscale,lim2o,lim1o;
1431  Double_t lim2=xQs;
1432  Double_t lim1=10000;
1433 
1434  //generation of kT according to 1/kT^4, with minimum QS=2 GeV and maximum ~sqrt(ptjet*T)
1435  fTf1Kt= new TF1("fTf1Kt","1/(x*x*x*x)",lim2,lim1);
1436  kTscale=fTf1Kt->GetRandom();
1437  //generation within the jet cone
1438 
1439  //generation of w according to 1/w, with minimum wc
1440  //omega needs to be larger than kT so to have well defined angles
1441  lim2o=kTscale;
1442  lim1o=kTscale/TMath::Sin(0.1);
1443  fTf1Omega= new TF1("fTf1Omega","1/x",lim2o,lim1o);
1444  omega=fTf1Omega->GetRandom();
1445  sinpptheta=kTscale/omega;
1446  pptheta=TMath::ASin(sinpptheta);
1447  if(pptheta>fJetRadius) continue;
1448 
1449  //Lorentz vector in the frame where the jet moves along z axis
1450  TLorentzVector pTrackCMS(kTscale/TMath::Sqrt(2),kTscale/TMath::Sqrt(2),omega*TMath::Cos(pptheta),omega);
1451  TVector3 MyJet(fJet->Px(),fJet->Py(),fJet->Pz());
1452  TVector3 direction = MyJet.Unit();
1453  //rotate the track to the jet frame
1454  pTrackCMS.RotateUz(direction);
1455 
1456  //add the rotated track to the jet
1457  PseudoTracksLab.reset(pTrackCMS.Px(),pTrackCMS.Py(),pTrackCMS.Pz(),pTrackCMS.E());
1458 
1459  PseudoTracksLab.set_user_index(i+fJet->GetNumberOfTracks()+100);
1460 
1461  omega2=PseudoTracksLab.perp();
1462  angle2=pTrackCMS.Angle(MyJet);
1463 
1464  fInputVectors.push_back(PseudoTracksLab);
1465  xflagAdded=1;
1466  }
1467 
1468 
1469 
1470 
1471  fastjet::JetAlgorithm jetalgo(fastjet::antikt_algorithm);
1472  if(ReclusterAlgo==0){ xflagalgo=0.5;
1473  jetalgo=fastjet::kt_algorithm ;}
1474 
1475  if(ReclusterAlgo==1){ xflagalgo=1.5;
1476  jetalgo=fastjet::cambridge_algorithm;}
1477  if(ReclusterAlgo==2){ xflagalgo=2.5;
1478  jetalgo=fastjet::antikt_algorithm;}
1479 
1480  fastjet::JetDefinition fJetDef(jetalgo, 1., static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
1481 
1482  try {
1483  fastjet::ClusterSequence fClustSeqSA(fInputVectors, fJetDef);
1484  std::vector<fastjet::PseudoJet> fOutputJets;
1485  fOutputJets.clear();
1486  fOutputJets=fClustSeqSA.inclusive_jets(0);
1487 
1488  fastjet::PseudoJet jj;
1489  fastjet::PseudoJet j1;
1490  fastjet::PseudoJet j2;
1491  jj=fOutputJets[0];
1492  double ndepth=0;
1493  double nsd=0;
1494  double nall=0;
1495  fShapesVar[10]=fOutputJets[0].perp();
1496  while(jj.has_parents(j1,j2)){
1497 
1498  if(j1.perp() < j2.perp()) swap(j1,j2);
1499  double delta_R=j1.delta_R(j2);
1500  double z=j2.perp()/(j1.perp()+j2.perp());
1501  double y =log(1.0/delta_R);
1502  double lnpt_rel=log(z*delta_R);
1503  nall=nall+1;
1504  if(z>fHardCutoff){
1505  ndepth=ndepth+1;
1506  nsd=nsd+1;
1507  if(nsd==1 && ReclusterAlgo==1) fShapesVar[9]=z;
1508  Double_t LundEntries[6] = {y,lnpt_rel,fOutputJets[0].perp(),xflagalgo,partonFlavor,ndepth};
1509  fHLundIterative->Fill(LundEntries);}
1510  jj=j1;}
1511 
1512  if(ReclusterAlgo==1){
1513  fShapesVar[6]=nsd;
1514  fShapesVar[7]=nall;}
1515 
1516  if(fAdditionalTracks>0 && xflagAdded>0){
1517  zinject=omega2/fOutputJets[0].perp();
1518  angleinject=angle2;
1519 
1520  yinject =log(1.0/angleinject);
1521  lnpt_relinject=log(zinject*angleinject);
1522  Double_t LundEntriesInject[4] = {yinject,lnpt_relinject,fOutputJets[0].perp(),fJet->Pt()};
1523  fHLundIterativeInject->Fill(LundEntriesInject);}
1524 
1525 
1526  } catch (fastjet::Error) {
1527  AliError(" [w] FJ Exception caught.");
1528  //return -1;
1529  }
1530 
1531 
1532  if(fTf1SoftKt){ delete fTf1SoftKt;}
1533  if(fTf1SoftOmega){ delete fTf1SoftOmega;}
1534  if(fTf1Kt){ delete fTf1Kt;}
1535  if(fTf1Omega){ delete fTf1Omega;}
1536 
1537 
1538  return;
1539 
1540 
1541 }
1542 
1543 
1545  const Double_t jetradius,
1546  const Double_t subjetradius,
1547  const char *ntracksPartLevel,
1548  const char *type,
1549  const char *CentEst,
1550  Int_t pSel,
1551  TString trigClass,
1552  TString kEmcalTriggers,
1553  TString tag,
1554  const char *rhoName,
1558  Float_t minpTHTrigger, Float_t maxpTHTrigger,
1560 
1561 
1562 
1563  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1564  if (!mgr)
1565  {
1566  ::Error("AliAnalysisTaskEmcalJetShapesMC","No analysis manager found.");
1567  return 0;
1568  }
1569  Bool_t ismc=kFALSE;
1570  ismc = (mgr->GetMCtruthEventHandler())?kTRUE:kFALSE;
1571 
1572  // Check the analysis type using the event handlers connected to the analysis manager.
1573  //==============================================================================
1574  if (!mgr->GetInputEventHandler())
1575  {
1576  ::Error("AddTaskJetShapesMC", "This task requires an input event handler");
1577  return NULL;
1578  }
1579 
1580  TString wagonName1 = Form("JetShapesMC_%s_Histos%s%s",njetsBase,trigClass.Data(),tag.Data());
1581  TString wagonName2 = Form("JetShapesMC_%s_Tree%s%s",njetsBase,trigClass.Data(),tag.Data());
1582 
1583  //Configure jet tagger task
1585 
1586  //task->SetNCentBins(4);
1587  task->SetJetShapeType(jetShapeType);
1588  task->SetJetShapeSub(jetShapeSub);
1589  task->SetJetSelection(jetSelection);
1590  task->SetDerivativeSubtractionOrder(derivSubtrOrder);
1591  task->SetJetRadius(jetradius);
1592  task->SetSubjetRadius(subjetradius);
1593 
1594  if (jetSelection == AliAnalysisTaskEmcalJetShapesMC::kRecoil) task->SetPtTriggerSelections(minpTHTrigger, maxpTHTrigger);
1595 
1596  TString thename(njetsBase);
1597 
1598  AliParticleContainer *trackContPartLevel = task->AddMCParticleContainer(ntracksPartLevel);
1599 
1600  AliJetContainer *jetContBase=0x0;
1601  TString strType(type);
1602 
1603  jetContBase = task->AddJetContainer(njetsBase,strType,jetradius);
1604  if(jetContBase) {
1605  jetContBase->SetRhoName(rhoName);
1606  jetContBase->ConnectParticleContainer(trackContPartLevel);
1607  //jetContBase->ConnectClusterContainer(clusterCont);
1608  jetContBase->SetPercAreaCut(0.6);
1609  }
1610 
1611  task->SetCaloTriggerPatchInfoName(kEmcalTriggers.Data());
1612  task->SetCentralityEstimator(CentEst);
1613  task->SelectCollisionCandidates(pSel);
1614  task->SetUseAliAnaUtils(kFALSE);
1615 
1616  mgr->AddTask(task);
1617 
1618  //Connnect input
1619  mgr->ConnectInput (task, 0, mgr->GetCommonInputContainer() );
1620 
1621  //Connect output
1622  TString contName1(wagonName1);
1623  TString contName2(wagonName2);
1624 
1625  if (jetShapeType == AliAnalysisTaskEmcalJetShapesMC::kGenShapes) {
1626  contName1 += "_GenShapes";
1627  contName2 += "_GenShapes";
1628  }
1629 
1630  switch (jetShapeSub) {
1631 
1633  contName1 += "_NoSub";
1634  contName2 += "_NoSub";
1635  break;
1636 
1638  contName1 += "_ConstSub";
1639  contName2 += "_ConstSub";
1640  break;
1641 
1643  contName1 += "_DerivSub";
1644  contName2 += "_DerivSub";
1645  break;
1646 
1647 
1648  }
1649 
1650  switch (jetSelection) {
1652  contName1 += "_Incl";
1653  contName2 += "_Incl";
1654  break;
1655 
1656 
1658  TString recoilTriggerString = Form("_Recoil_%.0f_%0.f", minpTHTrigger, maxpTHTrigger);
1659  contName1 += recoilTriggerString;
1660  contName2 += recoilTriggerString;
1661 
1662  break;
1663 
1664  }
1665 
1666 
1667  TString outputfile = Form("%s",AliAnalysisManager::GetCommonFileName());
1668 
1669 
1670  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contName1.Data(),TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
1671  AliAnalysisDataContainer *coutput2 = mgr->CreateContainer(contName2.Data(),TTree::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
1672 
1673  mgr->ConnectOutput(task,1,coutput1);
1674  mgr->ConnectOutput(task,2,coutput2);
1675 
1676  return task;
1677 
1678 }
1679 
1680 
1681 
void SetRadius(Double_t val)
Float_t GetSigma2(AliEmcalJet *jet, Int_t jetContNb=0)
Double_t GetFirstOrderSubtractedAngularity() const
Float_t GetJetMass(AliEmcalJet *jet, Int_t jetContNb=0)
double Double_t
Definition: External.C:58
void SetJetMinPt(Double_t val)
Double_t GetSecondOrderSubtractedSigma2() const
Definition: External.C:236
Float_t GetPartonEta6() const
AliJetContainer * GetJetContainer(Int_t i=0) const
Float_t GetPartonEta7() const
Double_t GetSecondOrderSubtractedConstituent() const
Double_t Eta() const
Definition: AliEmcalJet.h:121
Double_t FjNSubJettiness(AliEmcalJet *Jet, Int_t JetContNb, Int_t N, Int_t Algorithm, Double_t Beta, Int_t Option, Double_t Beta_SD=0.0, Double_t ZCut=0.1, Int_t SoftDropOn=0)
Double_t Py() const
Definition: AliEmcalJet.h:107
Double_t Phi() const
Definition: AliEmcalJet.h:117
Float_t GetPythiaEventWeight() const
void RecursiveParents(AliEmcalJet *fJet, AliJetContainer *fJetCont, Int_t ReclusterAlgo, Float_t PartonFlavor)
Float_t GetPartonPhi7() const
void SetPercAreaCut(Float_t p)
Float_t CoreFrac(AliEmcalJet *jet, Int_t jetContNb=0)
Float_t GetJetCircularity(AliEmcalJet *jet, Int_t jetContNb=0)
Double_t E() const
Definition: AliEmcalJet.h:119
TF1 * fTf1SoftOmega
to generate kT according to BDMS tail
Float_t GetJetpTD(AliEmcalJet *jet, Int_t jetContNb=0)
Int_t GetPartonFlag7() const
Container for particles within the EMCAL framework.
Float_t GetPartonPhi6() const
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:160
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:139
Double_t Px() const
Definition: AliEmcalJet.h:106
Float_t GetJetCoreFrac(AliEmcalJet *jet, Int_t jetContNb=0)
Double_t Nsubjettiness(AliEmcalJet *pJet, AliJetContainer *pContJets, Double_t dVtx[3], Int_t N, Int_t Algorithm, Double_t Radius, Double_t Beta, Int_t Option=0, Int_t Measure=0, Double_t Beta_SD=0, Double_t ZCut=0.1, Int_t SoftDropOn=0)
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
Int_t GetPartonFlag6() const
void SetRhoName(const char *n)
AliParticleContainer * GetParticleContainer() const
Double_t GetFirstOrderSubtractedConstituent() const
void SetRecombSheme(Int_t val)
static AliAnalysisTaskEmcalJetShapesMC * AddTaskJetShapesMC(const char *njetsBase, const Double_t jetradius, const Double_t subjetradius, const char *ntracksPartLevel, const char *type, const char *CentEst, Int_t pSel, TString trigClass="", TString kEmcalTriggers="", TString tag="", const char *rhoName="", AliAnalysisTaskEmcalJetShapesMC::JetShapeType jetShapeType=AliAnalysisTaskEmcalJetShapesMC::kGenShapes, AliAnalysisTaskEmcalJetShapesMC::JetShapeSub jetShapeSub=AliAnalysisTaskEmcalJetShapesMC::kNoSub, AliAnalysisTaskEmcalJetShapesMC::JetSelectionType jetSelection=AliAnalysisTaskEmcalJetShapesMC::kInclusive, Float_t minpTHTrigger=0., Float_t maxpTHTrigger=0., AliAnalysisTaskEmcalJetShapesMC::DerivSubtrOrder derivSubtrOrder=AliAnalysisTaskEmcalJetShapesMC::kSecondOrder)
void SetJetAlgorithm(Int_t val)
int Int_t
Definition: External.C:63
Double_t Eta() const
Double_t Pt() const
unsigned int UInt_t
Definition: External.C:33
Float_t GetJetNumberOfConstituents(AliEmcalJet *jet, Int_t jetContNb=0)
float Float_t
Definition: External.C:68
Float_t fqhat
Random number generator.
Task to store and correlate the MC shapes.
TF1 * fTf1SoftKt
to generate omega for soft background
TTree * fTreeObservableTagging
! Tree with tagging variables subtracted MC or true MC or raw
Double_t RelativePhi(Double_t mphi, Double_t vphi)
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.
Double_t GetSecondOrderSubtractedAngularity() const
Double_t GetSubjetFraction(AliEmcalJet *Jet, Int_t JetContNb, Double_t JetRadius, AliEmcalJetFinder *Reclusterer)
Definition: Option.C:68
AliEmcalJetFinder * Recluster(AliEmcalJet *Jet, Int_t JetContNb, Double_t JetRadius, Double_t SubJetRadius, Double_t SubJetMinPt, Int_t Algorithm, const char *Name)
Double_t fCent
!event centrality
Float_t PTD(AliEmcalJet *jet, Int_t jetContNb=0)
void SetJetMaxEta(Double_t val)
AliEmcalJet * GetNextAcceptJet()
TF1 * fTf1Kt
to generate omega according to BDMPS tail
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)
Float_t GetJetAngularity(AliEmcalJet *jet, Int_t jetContNb=0)
Int_t SelectTrigger(Float_t minpT, Float_t maxpT)
AliEmcalJet * GetJet(Int_t index)
void ConnectParticleContainer(AliParticleContainer *c)
Double_t Pt() const
Definition: AliEmcalJet.h:109
Enhanced TList-derived class that implements correct merging for pt_hard binned production.
Definition: AliEmcalList.h:25
Float_t Angularity(AliEmcalJet *jet, Int_t jetContNb=0)
Bool_t FillHistograms()
Function filling histograms.
Double_t GetFirstOrderSubtractedCircularity() const
Bool_t RetrieveEventObjects()
Retrieve common objects from event.
AliEmcalList * fOutput
!output list
Double_t SubJetOrdering(AliEmcalJet *Jet, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t Type, Bool_t Index)
void NTValues(AliEmcalJet *jet, Int_t jetContNb, Float_t *nTFractions)
const Int_t nVar
Double_t GetSecondOrderSubtractedCircularity() const
Double_t DeltaR(const AliVParticle *part1, const AliVParticle *part2)
Float_t GetJetLeSub(AliEmcalJet *jet, Int_t jetContNb=0)
Store some informaion about a Pythia eventThis class is used to store some information about a Pythia...
void SetMakeGeneralHistograms(Bool_t g)
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
const AliEmcalPythiaInfo * GetPythiaInfo() const
Double_t Pz() const
Definition: AliEmcalJet.h:108
Float_t Circularity(AliEmcalJet *jet, Int_t jetContNb=0)
void SoftDrop(AliEmcalJet *fJet, AliJetContainer *fJetCont, Double_t zcut, Double_t beta, Int_t ReclusterAlgo)
Declaration of class AliEmcalPythiaInfo.
const char Option_t
Definition: External.C:48
void UserCreateOutputObjects()
Main initialization function on the worker.
Double_t GetFirstOrderSubtractedSigma2() const
Float_t Sigma2(AliEmcalJet *jet, Int_t jetContNb=0)
bool Bool_t
Definition: External.C:53
AliEmcalJetShapeProperties * GetShapeProperties() const
Definition: AliEmcalJet.h:361
Float_t GetPartonPt7() const
Float_t LeSub(AliEmcalJet *jet, Int_t jetContNb=0)
void swap(PWGJE::EMCALJetTasks::AliAnalysisTaskEmcalJetHPerformance &first, PWGJE::EMCALJetTasks::AliAnalysisTaskEmcalJetHPerformance &second)
Double_t M() const
Definition: AliEmcalJet.h:120
Container for jet within the EMCAL jet framework.
std::vector< int > GetPtSortedTrackConstituentIndexes(TClonesArray *tracks) const
Float_t GetPartonPt6() const