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