AliPhysics  1168478 (1168478)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 "AliEmcalJetFinder.h"
42 #include "AliAODEvent.h"
44 
45 using std::cout;
46 using std::endl;
47 
49 
50 //________________________________________________________________________
53  fContainer(0),
54  fMinFractionShared(0),
55  fJetShapeType(kGenShapes),
56  fJetShapeSub(kNoSub),
57  fJetSelection(kInclusive),
58  fPtThreshold(-9999.),
59  fRMatching(0.2),
60  fJetRadius(0.4),
61  fSubjetRadius(0.2),
62  fSelectedShapes(0),
63  fminpTTrig(20.),
64  fmaxpTTrig(50.),
65  fangWindowRecoil(0.6),
66  fSemigoodCorrect(0),
67  fHolePos(0),
68  fHoleWidth(0),
69  fCentSelectOn(kTRUE),
70  fCentMin(0),
71  fCentMax(10),
72  fOneConstSelectOn(kFALSE),
73  fDerivSubtrOrder(0),
74  fPhiJetCorr6(0x0),
75  fPhiJetCorr7(0x0),
76  fEtaJetCorr6(0x0),
77  fEtaJetCorr7(0x0),
78  fPtJetCorr(0x0),
79  fPtJet(0x0),
80  fhpTjetpT(0x0),
81  fhPt(0x0),
82  fhPhi(0x0),
83  fNbOfConstvspT(0x0),
84  fTreeObservableTagging(0x0)
85 
86 {
87  for(Int_t i=0;i<33;i++){
88  fShapesVar[i]=0;}
89  SetMakeGeneralHistograms(kTRUE);
90 }
91 
92 //________________________________________________________________________
94  AliAnalysisTaskEmcalJet(name, kTRUE),
95  fContainer(0),
96  fMinFractionShared(0),
97  fJetShapeType(kGenShapes),
98  fJetShapeSub(kNoSub),
99  fJetSelection(kInclusive),
100  fPtThreshold(-9999.),
101  fRMatching(0.2),
102  fSelectedShapes(0),
103  fminpTTrig(20.),
104  fmaxpTTrig(50.),
105  fangWindowRecoil(0.6),
106  fSemigoodCorrect(0),
107  fHolePos(0),
108  fHoleWidth(0),
109  fCentSelectOn(kTRUE),
110  fCentMin(0),
111  fCentMax(10),
112  fOneConstSelectOn(kFALSE),
113  fDerivSubtrOrder(0),
114  fPhiJetCorr6(0x0),
115  fPhiJetCorr7(0x0),
116  fEtaJetCorr6(0x0),
117  fEtaJetCorr7(0x0),
118  fPtJetCorr(0x0),
119  fPtJet(0x0),
120  fhpTjetpT(0x0),
121  fhPt(0x0),
122  fhPhi(0x0),
123  fNbOfConstvspT(0x0),
124  fTreeObservableTagging(0x0)
125 
126 {
127  // Standard constructor.
128 
129 
130  for(Int_t i=0;i<33;i++){
131  fShapesVar[i]=0;}
132 
134 
135  DefineOutput(1, TList::Class());
136  DefineOutput(2, TTree::Class());
137 
138 
139 }
140 
141 //________________________________________________________________________
143 {
145  delete fTreeObservableTagging;
147  }
148 
149 }
150 
151 //________________________________________________________________________
153 {
154  // Create user output.
156 
157  Bool_t oldStatus = TH1::AddDirectoryStatus();
158  TH1::AddDirectory(oldStatus);
159 
160  //fTreeObservableTagging = new TTree("fTreeJetShape", "fTreeJetShape");
161 
162  //TH1::AddDirectory(oldStatus);
163 
164  const char* nameoutput = GetOutputSlot(2)->GetContainer()->GetName();
165  fTreeObservableTagging = new TTree(nameoutput, nameoutput);
166 
167  const Int_t nVar = 33;
168 
169  TString *fShapesVarNames = new TString [nVar];
170 
171  fShapesVarNames[0] = "partonCode";
172  fShapesVarNames[1] = "ptJet";
173  fShapesVarNames[2] = "ptDJet";
174  fShapesVarNames[3] = "mJet";
175  fShapesVarNames[4] = "nbOfConst";
176  fShapesVarNames[5] = "angularity";
177  fShapesVarNames[6] = "Nsubjet1kt";
178  fShapesVarNames[7] = "Nsubjet2kt";
179  fShapesVarNames[8] = "Nsubjet1Min";
180  fShapesVarNames[9] = "Nsubjet2Min";
181  fShapesVarNames[10] = "DeltaRkt";
182  fShapesVarNames[11] = "DeltaRMin";
183  fShapesVarNames[12] = "SDSymm";
184  fShapesVarNames[13] = "SDDeltaR";
185  fShapesVarNames[14] = "SDGroomedFrac";
186  fShapesVarNames[15] = "SDGroomedN";
187  fShapesVarNames[16] = "SDSymmkt";
188  fShapesVarNames[17] = "SDDeltaRkt";
189  fShapesVarNames[18] = "SDGroomedFrackt";
190  fShapesVarNames[19] = "SDGroomedNkt";
191  fShapesVarNames[20] = "SDSymmAkt";
192  fShapesVarNames[21] = "SDDeltaRAkt";
193  fShapesVarNames[22] = "SDGroomedFracAkt";
194  fShapesVarNames[23] = "SDGroomedNAkt";
195  fShapesVarNames[24] = "SDSymmBeta1";
196  fShapesVarNames[25] = "SDDeltaRBeta1";
197  fShapesVarNames[26] = "SDGroomedFracBeta1";
198  fShapesVarNames[27] = "SDGroomedNBeta1";
199  fShapesVarNames[28] = "SDSymmBeta2";
200  fShapesVarNames[29] = "SDDeltaRBeta2";
201  fShapesVarNames[30] = "SDGroomedFracBeta2";
202  fShapesVarNames[31] = "SDGroomedNBeta2";
203  fShapesVarNames[32] = "weightPythia";
204 
205 
206 
207  //fShapesVarNames[7] = "lesub";
208  //fShapesVarNames[8] = "CoreFraction";
209  //fShapesVarNames[9] = "Nsubjet1";
210  //fShapesVarNames[10] = "Nsubjet2";
211  //fShapesVarNames[11] = "DeltaR";
212  //fShapesVarNames[12] = "OpenAngle";
213  //fShapesVarNames[13] = "weightPythia";
214 
215  //fShapesVarNames[14] = "NT70";
216  //fShapesVarNames[15] = "nConstNT70";
217  //fShapesVarNames[16] = "NT80";
218  //fShapesVarNames[17] = "nConstNT80";
219  //fShapesVarNames[18] = "NT90";
220  //fShapesVarNames[19] = "nConstNT90";
221  //fShapesVarNames[20] = "NT95";
222  //fShapesVarNames[21] = "nConstNT95";
223 
224  //fShapesVarNames[22] = "SubjetFraction";
225 
226 
227  for(Int_t ivar=0; ivar < nVar; ivar++){
228  cout<<"looping over variables"<<endl;
229  fTreeObservableTagging->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/F", fShapesVarNames[ivar].Data()));}
230 
231 
232 
233  fPhiJetCorr6= new TH2F("fPhiJetCorr6", "fPhiJetCorr6", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
234  fOutput->Add(fPhiJetCorr6);
235  fEtaJetCorr6= new TH2F("fEtaJetCorr6", "fEtaJetCorr6", 50, -1.5, 1.5, 50, -1.5, 1.5);
236  fOutput->Add(fEtaJetCorr6);
237 
238  fPhiJetCorr7= new TH2F("fPhiJetCorr7", "fPhiJetCorr7", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
239  fOutput->Add(fPhiJetCorr7);
240  fEtaJetCorr7= new TH2F("fEtaJetCorr7", "fEtaJetCorr7", 50, -1.5, 1.5, 50, -1.5, 1.5);
241  fOutput->Add(fEtaJetCorr7);
242 
243  fPtJetCorr= new TH2F("fPtJetCorr", "fPtJetCorr", 100, 0, 200, 100, 0, 200);
244  fOutput->Add(fPtJetCorr);
245  fPtJet= new TH1F("fPtJet", "fPtJet", 100, 0, 200);
246  fOutput->Add(fPtJet);
247 
248  fhpTjetpT= new TH2F("fhpTjetpT", "fhpTjetpT", 200, 0, 200, 200, 0, 200);
249  fOutput->Add(fhpTjetpT);
250  fhPt= new TH1F("fhPt", "fhPt", 200, 0, 200);
251  fOutput->Add(fhPt);
252  fhPhi= new TH1F("fhPhi", "fhPhi", 100, -TMath::Pi(), TMath::Pi());
253  fOutput->Add(fhPhi);
254 
255  fNbOfConstvspT=new TH2F("fNbOfConstvspT", "fNbOfConstvspT", 100, 0, 100, 200, 0, 200);
256  fOutput->Add(fNbOfConstvspT);
257 
258  //fOutput->Add(fTreeObservableTagging);
259 
260 
261  PostData(1, fOutput); // Post data for ALL output slots > 0 here
262  PostData(2, fTreeObservableTagging);
263 
264  delete [] fShapesVarNames;
265 
266 }
267 
268 //________________________________________________________________________
270 {
271  // Run analysis code here, if needed. It will be executed before FillHistograms().
272 
273  return kTRUE;
274 }
275 
276 //________________________________________________________________________
278 {
279  // Fill histograms.
280  //cout<<"IntoFillHistograms"<<endl;
281  AliEmcalJet* jet1 = NULL;
282  AliJetContainer *jetCont = GetJetContainer(0);
283 
284  Float_t kWeight=1;
285  if (fCentSelectOn)
286  if ((fCent>fCentMax) || (fCent<fCentMin)) return 0;
287 
288  AliAODTrack *triggerHadron = 0x0;
289 
290  if (fJetSelection == kRecoil) {
291  //Printf("Recoil jets!!!, fminpTTrig = %f, fmaxpTTrig = %f", fminpTTrig, fmaxpTTrig);
292  Int_t triggerHadronLabel = SelectTrigger(fminpTTrig, fmaxpTTrig);
293 
294 
295  if (triggerHadronLabel==-99999) {
296  //Printf ("Trigger Hadron not found, return");
297  return 0;}
298 
299 
301  TClonesArray *trackArrayAn = partContAn->GetArray();
302  triggerHadron = static_cast<AliAODTrack*>(trackArrayAn->At(triggerHadronLabel));
303 
304  if (!triggerHadron) {
305  //Printf("No Trigger hadron with the found label!!");
306  return 0;
307  }
308 
309  if(fSemigoodCorrect){
310  Double_t disthole=RelativePhi(triggerHadron->Phi(),fHolePos);
311  if(TMath::Abs(disthole)+fHoleWidth>TMath::Pi()-fangWindowRecoil){
312  return 0;}
313  }
314 
315  fhPt->Fill(triggerHadron->Pt());
316 
317  }
318 
319  if(jetCont) {
320  jetCont->ResetCurrentID();
321  while((jet1 = jetCont->GetNextAcceptJet())) {
322  //Printf("jet1=%p", jet1);
323  if (!jet1) continue;
324  AliEmcalJet* jet2 = 0x0;
325  AliEmcalJet* jet3 = 0x0;
326  fPtJet->Fill(jet1->Pt());
327  AliEmcalJet *jetUS = NULL;
328  Int_t ifound=0, jfound=0;
329  Int_t ilab=-1, jlab=-1;
330 
332  Double_t disthole=RelativePhi(jet1->Phi(),fHolePos);
333  if(TMath::Abs(disthole)<fHoleWidth){
334  continue;
335  }
336  }
337 
338  Float_t dphiRecoil = 0.;
339  if (fJetSelection == kRecoil){
340  dphiRecoil = RelativePhi(triggerHadron->Phi(), jet1->Phi());
341  if (TMath::Abs(dphiRecoil) < (TMath::Pi() - fangWindowRecoil)) {
342  // Printf("Recoil jets back to back not found! continuing");
343  continue;
344  }
345 
346  fhpTjetpT->Fill(triggerHadron->Pt(), jet1->Pt());
347  //Printf(" ************ FILLING HISTOS****** shapeSub = %d, triggerHadron = %f, jet1 = %f", fJetShapeSub, triggerHadron->Pt(), jet1->Pt());
348  fhPhi->Fill(RelativePhi(triggerHadron->Phi(), jet1->Phi()));
349 
350  }
351 
352 
353  fShapesVar[0] = 0.;
354 
355  if (fJetShapeType == kGenShapes){
356  const AliEmcalPythiaInfo *partonsInfo = 0x0;
357  partonsInfo = GetPythiaInfo();
358  //Printf("partonsInfo=%p", partonsInfo);
359  Double_t jp1=RelativePhi(jet1->Phi(),partonsInfo->GetPartonPhi6());
360  Double_t detap1=(jet1->Eta())-(partonsInfo->GetPartonEta6());
361  kWeight=partonsInfo->GetPythiaEventWeight();
362  //Printf("kWeight=%f", kWeight);
363  fShapesVar[32] = kWeight;
364 
365  Float_t dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
366  fEtaJetCorr6->Fill(jet1->Eta(), partonsInfo->GetPartonEta6());
367  fPhiJetCorr6->Fill(jet1->Phi(), partonsInfo->GetPartonPhi6());
368  if(dRp1 < fRMatching) {
369  fShapesVar[0] = partonsInfo->GetPartonFlag6();
370  fPtJetCorr ->Fill(partonsInfo->GetPartonPt6(), jet1->Pt());
371  }
372  else {
373  jp1=RelativePhi(jet1->Phi(),partonsInfo->GetPartonPhi7());
374  detap1=(jet1->Eta())-(partonsInfo->GetPartonEta7());
375  dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
376  fEtaJetCorr7->Fill(jet1->Eta(), partonsInfo->GetPartonEta7());
377  fPhiJetCorr7->Fill(jet1->Phi(), partonsInfo->GetPartonPhi7());
378  if(dRp1 < fRMatching) {
379  fShapesVar[0] = partonsInfo->GetPartonFlag7();
380  fPtJetCorr->Fill(partonsInfo->GetPartonPt7(), jet1->Pt());
381  }
382  else fShapesVar[0]=0;
383  }
384  }
385 
386  Double_t ptSubtracted = 0;
387  ptSubtracted= jet1->Pt();
388  //Printf("ptSubtracted=%f", ptSubtracted);
389 
390 
391  if (ptSubtracted < fPtThreshold) continue;
392 
393  if ((fCentSelectOn == kFALSE) && (jet1->GetNumberOfTracks() < 2)) continue;
394 
395  //AliEmcalJetFinder *Reclusterer1; //Object containg Subjets from Subtracted Hybrid Jets
396  //Reclusterer1 = Recluster(jet1, 0, fJetRadius, fSubjetRadius, 1, 0, "SubJetFinder_1");
397 
398 
399 
400 
401 
402 
403 
404  fShapesVar[1] = ptSubtracted;
405  fShapesVar[2] = GetJetpTD(jet1,0);
406  fShapesVar[3] = GetJetMass(jet1,0);
407  fShapesVar[4] = 1.*GetJetNumberOfConstituents(jet1,0);
408  fShapesVar[5] = GetJetAngularity(jet1,0);
409  //nsub1 and nsub2 for kT
410  fShapesVar[6] = fjNSubJettiness(jet1,0,1,0,1,0);
411  fShapesVar[7] = fjNSubJettiness(jet1,0,2,0,1,0);
412  //nsub1 and nsub2 for min_axis
413  fShapesVar[8] =0;
414  // fjNSubJettiness(jet1,0,1,10,1,0);
415  fShapesVar[9] =0;
416  // fjNSubJettiness(jet1,0,2,10,1,0);
417  //deltaRkt
418  fShapesVar[10] = 0;
419  //fjNSubJettiness(jet1,0,2,0,1,2);
420  //deltaRmin
421  fShapesVar[11] = 0;
422  //fjNSubJettiness(jet1,0,2,10,1,2);
423 
424  //SoftDropParameters for different reclustering strategies and beta values
425  SoftDrop(jet1,jetCont,0.1,0,0);
426  SoftDrop(jet1,jetCont,0.1,0,1);
427  SoftDrop(jet1,jetCont,0.1,0,2);
428  SoftDrop(jet1,jetCont,0.1,1,0);
429  SoftDrop(jet1,jetCont,0.1,2,0);
430 
431 
432  // Float_t nTFractions[8]={0.,0.,0.,0.,0.,0.,0.,0.};
433  //NTValues(jet1, 0, nTFractions);
434  //shape 13 is pythia weight!
435  //for (Int_t ishape=14; ishape<22; ishape++) fShapesVar[ishape] = nTFractions[ishape-14];
436 
437  //fShapesVar[22]= GetSubjetFraction(jet1,0,fJetRadius,Reclusterer1);
438 
439  fTreeObservableTagging->Fill();
440 
441  }
442 
443  }
444 
445  return kTRUE;
446 }
447 
448 //________________________________________________________________________
450  //calc subtracted jet mass
451  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
453  else return jet->GetShapeProperties()->GetSecondOrderSubtracted();
454  else
455  return jet->M();
456 }
457 
458 //________________________________________________________________________
460 
461  AliJetContainer *jetCont = GetJetContainer(jetContNb);
462  if (!jet->GetNumberOfTracks())
463  return 0;
464  Double_t den=0.;
465  Double_t num = 0.;
466  AliVParticle *vp1 = 0x0;
467  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
468  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
469 
470  if (!vp1){
471  Printf("AliVParticle associated to constituent not found");
472  continue;
473  }
474 
475  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
476  Double_t dr2 = (vp1->Eta()-jet->Eta())*(vp1->Eta()-jet->Eta()) + dphi*dphi;
477  Double_t dr = TMath::Sqrt(dr2);
478  num=num+vp1->Pt()*dr;
479  den=den+vp1->Pt();
480  }
481  return num/den;
482 }
483 
484 //________________________________________________________________________
486 
487  if((fJetShapeSub==kDerivSub) && (jetContNb==0))
490  else
491  return Angularity(jet, jetContNb);
492 
493 }
494 
495 
496 //________________________________________________________________________
498 
499  AliJetContainer *jetCont = GetJetContainer(jetContNb);
500  if (!jet->GetNumberOfTracks())
501  return 0;
502  Double_t den=0.;
503  Double_t num = 0.;
504  AliVParticle *vp1 = 0x0;
505  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
506  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
507 
508  if (!vp1){
509  Printf("AliVParticle associated to constituent not found");
510  continue;
511  }
512 
513  num=num+vp1->Pt()*vp1->Pt();
514  den=den+vp1->Pt();
515  }
516  return TMath::Sqrt(num)/den;
517 }
518 
519 //________________________________________________________________________
521  //calc subtracted jet mass
522  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
524  else return jet->GetShapeProperties()->GetSecondOrderSubtractedpTD();
525  else
526  return PTD(jet, jetContNb);
527 
528 }
529 
530 //_____________________________________________________________________________
532 
533  AliJetContainer *jetCont = GetJetContainer(jetContNb);
534  if (!jet->GetNumberOfTracks())
535  return 0;
536  Double_t mxx = 0.;
537  Double_t myy = 0.;
538  Double_t mxy = 0.;
539  int nc = 0;
540  Double_t sump2 = 0.;
541  Double_t pxjet=jet->Px();
542  Double_t pyjet=jet->Py();
543  Double_t pzjet=jet->Pz();
544 
545 
546  //2 general normalized vectors perpendicular to the jet
547  TVector3 ppJ1(pxjet, pyjet, pzjet);
548  TVector3 ppJ3(- pxjet* pzjet, - pyjet * pzjet, pxjet * pxjet + pyjet * pyjet);
549  ppJ3.SetMag(1.);
550  TVector3 ppJ2(-pyjet, pxjet, 0);
551  ppJ2.SetMag(1.);
552  AliVParticle *vp1 = 0x0;
553  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
554  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
555 
556  if (!vp1){
557  Printf("AliVParticle associated to constituent not found");
558  continue;
559  }
560 
561  TVector3 pp(vp1->Px(), vp1->Py(), vp1->Pz());
562 
563  //local frame
564  TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
565  TVector3 pPerp = pp - pLong;
566  //projection onto the two perpendicular vectors defined above
567 
568  Float_t ppjX = pPerp.Dot(ppJ2);
569  Float_t ppjY = pPerp.Dot(ppJ3);
570  Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
571  if(ppjT<=0) return 0;
572 
573  mxx += (ppjX * ppjX / ppjT);
574  myy += (ppjY * ppjY / ppjT);
575  mxy += (ppjX * ppjY / ppjT);
576  nc++;
577  sump2 += ppjT;}
578 
579  if(nc<2) return 0;
580  if(sump2==0) return 0;
581  // Sphericity Matrix
582  Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
583  TMatrixDSym m0(2,ele);
584 
585  // Find eigenvectors
586  TMatrixDSymEigen m(m0);
587  TVectorD eval(2);
588  TMatrixD evecm = m.GetEigenVectors();
589  eval = m.GetEigenValues();
590  // Largest eigenvector
591  int jev = 0;
592  // cout<<eval[0]<<" "<<eval[1]<<endl;
593  if (eval[0] < eval[1]) jev = 1;
594  TVectorD evec0(2);
595  // Principle axis
596  evec0 = TMatrixDColumn(evecm, jev);
597  Double_t compx=evec0[0];
598  Double_t compy=evec0[1];
599  TVector2 evec(compx, compy);
600  Double_t circ=0;
601  if(jev==1) circ=2*eval[0];
602  if(jev==0) circ=2*eval[1];
603 
604  return circ;
605 
606 
607 
608 }
609 
610 
611 
612 
613 //________________________________________________________________________
615  //calc subtracted jet mass
616 
617  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
620  else
621  return Circularity(jet, jetContNb);
622 
623 }
624 
625 //________________________________________________________________________
627 
628  AliJetContainer *jetCont = GetJetContainer(jetContNb);
629  if (!jet->GetNumberOfTracks())
630  return 0;
631  Double_t den=0.;
632  Double_t num = 0.;
633  AliVParticle *vp1 = 0x0;
634  AliVParticle *vp2 = 0x0;
635  std::vector<int> ordindex;
636  ordindex=jet->GetPtSortedTrackConstituentIndexes(jetCont->GetParticleContainer()->GetArray());
637  //Printf("Nbof const = %d", jet->GetNumberOfTracks());
638  //Printf("ordindex[0] = %d, ordindex[1] = %d", ordindex[0], ordindex[1]);
639 
640  if(ordindex.size()<2) return -1;
641 
642  vp1 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[0], jetCont->GetParticleContainer()->GetArray()));
643  if (!vp1){
644  Printf("AliVParticle associated to Leading constituent not found");
645  return -1;
646  }
647 
648  vp2 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[1], jetCont->GetParticleContainer()->GetArray()));
649  if (!vp2){
650  Printf("AliVParticle associated to Subleading constituent not found");
651  return -1;
652  }
653 
654 
655  num=vp1->Pt();
656  den=vp2->Pt();
657  //Printf("vp1->Pt() =%f, vp2->Pt() =%f", vp1->Pt(), vp2->Pt());
658 
659 return num-den;
660 }
661 
662 //________________________________________________________________________
664  //calc subtracted jet mass
665 
666  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
669  else
670  return LeSub(jet, jetContNb);
671 
672 }
673 
674 //________________________________________________________________________
676  //calc subtracted jet mass
677 
678  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
681  else
682  return jet->GetNumberOfTracks();
683 
684  }
685 
686 
687 //________________________________________________________________________
688 AliEmcalJetFinder *AliAnalysisTaskEmcalJetShapesMC::Recluster(AliEmcalJet *Jet, Int_t JetContNb, Double_t JetRadius, Double_t SubJetRadius, Double_t SubJetMinPt, Int_t Algorithm, const char* Name){
689 
690  AliJetContainer *JetCont = GetJetContainer(JetContNb);
691  AliEmcalJetFinder *Reclusterer = new AliEmcalJetFinder(Name); //JetFinder Object for reclustered jets
692  Reclusterer->SetRadius(SubJetRadius);
693  Reclusterer->SetJetMinPt(SubJetMinPt);
694  Reclusterer->SetJetAlgorithm(Algorithm); //0 for anti-kt 1 for kt
695  Reclusterer->SetJetMaxEta(0.9-JetRadius);
696  Reclusterer->SetRecombSheme(0);
697  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
698 
699  //Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
700  Double_t dVtx[3]={0.,0.,0.};
701  if(Reclusterer->AliEmcalJetFinder::Filter(Jet, JetCont, dVtx)){;} //reclustering jet1 using the jetfinderobject Reclusterer
702  return Reclusterer;
703 }
704 
705 
706 
707 
708 //----------------------------------------------------------------------
710  AliJetContainer *JetCont = GetJetContainer(JetContNb);
711  AliEmcalJet *SubJet=NULL;
712  Double_t DeltaR1=0;
713  Double_t DeltaR2=0;
714  AliVParticle *JetParticle=0x0;
715  Double_t SubJetiness_Numerator = 0;
716  Double_t SubJetiness_Denominator = 0;
717  Double_t Index=-2;
718  if (Reclusterer->GetNumberOfJets() < 1) return -2;
719  Index=SubJetOrdering(Jet,Reclusterer,1,0,kTRUE);
720  if(Index==-999) return -2;
721  SubJetiness_Numerator=(Reclusterer->GetJet(Index)->Pt());
722  SubJetiness_Denominator=Jet->Pt();
723  return SubJetiness_Numerator/SubJetiness_Denominator;
724 
725 
726 }
727 //__________________________________________________________________________________
729 
730  AliJetContainer *jetCont = GetJetContainer(jetContNb);
731  if (!jet->GetNumberOfTracks())
732  return 0;
733  Double_t den=0.;
734  Double_t num = 0.;
735  AliVParticle *vp1 = 0x0;
736  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
737  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
738 
739  if (!vp1){
740  Printf("AliVParticle associated to constituent not found");
741  continue;
742  }
743 
744  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
745  Double_t dr2 = (vp1->Eta()-jet->Eta())*(vp1->Eta()-jet->Eta()) + dphi*dphi;
746  Double_t dr = TMath::Sqrt(dr2);
747  if(dr<=fSubjetRadius) num=num+vp1->Pt();
748 
749  }
750  return num/jet->Pt();
751 }
752 
753 
754 
755 
756 //________________________________________________________________________
758  //calc subtracted jet mass
759 
760  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
763  else
764  return CoreFrac(jet, jetContNb);
765 
766 }
767 
768 
769 
770 
771 //----------------------------------------------------------------------
773  AliEmcalJet *SubJet=NULL;
774  Double_t SortingVariable;
775  Int_t ArraySize =N+1;
776  TArrayD *JetSorter = new TArrayD(ArraySize);
777  TArrayD *JetIndexSorter = new TArrayD(ArraySize);
778  for (Int_t i=0; i<ArraySize; i++){
779  JetSorter->SetAt(0,i);
780  }
781  for (Int_t i=0; i<ArraySize; i++){
782  JetIndexSorter->SetAt(0,i);
783  }
784  if(Reclusterer->GetNumberOfJets()<N) return -999;
785  for (Int_t i=0; i<Reclusterer->GetNumberOfJets(); i++){
786  SubJet=Reclusterer->GetJet(i);
787  if (Type==0) SortingVariable=SubJet->Pt();
788  else if (Type==1) SortingVariable=SubJet->E();
789  else if (Type==2) SortingVariable=SubJet->M();
790  for (Int_t j=0; j<N; j++){
791  if (SortingVariable>JetSorter->GetAt(j)){
792  for (Int_t k=N-1; k>=j; k--){
793  JetSorter->SetAt(JetSorter->GetAt(k),k+1);
794  JetIndexSorter->SetAt(JetIndexSorter->GetAt(k),k+1);
795  }
796  JetSorter->SetAt(SortingVariable,j);
797  JetIndexSorter->SetAt(i,j);
798  break;
799  }
800  }
801  }
802  if (!Index) return JetSorter->GetAt(N-1);
803  else return JetIndexSorter->GetAt(N-1);
804 }
805 
806 
807 
808 //returns -1 if the Nth hardest jet is requested where N>number of available jets
809 //type: 0=Pt 1=E 2=M
810 //Index TRUE=returns index FALSE=returns value of quantatiy in question
811 
812 
813 
814 
815 
816 //______________________________________________________________________________
818 
819  AliJetContainer *jetCont = GetJetContainer(jetContNb);
820  if (!jet->GetNumberOfTracks())
821  return 0;
822  Double_t mxx = 0.;
823  Double_t myy = 0.;
824  Double_t mxy = 0.;
825  int nc = 0;
826  Double_t sump2 = 0.;
827 
828  AliVParticle *vp1 = 0x0;
829  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
830  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
831 
832  if (!vp1){
833  Printf("AliVParticle associated to constituent not found");
834  continue;
835  }
836 
837  Double_t ppt=vp1->Pt();
838  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
839 
840  Double_t deta = vp1->Eta()-jet->Eta();
841  mxx += ppt*ppt*deta*deta;
842  myy += ppt*ppt*dphi*dphi;
843  mxy -= ppt*ppt*deta*TMath::Abs(dphi);
844  nc++;
845  sump2 += ppt*ppt;
846 
847  }
848  if(nc<2) return 0;
849  if(sump2==0) return 0;
850  // Sphericity Matrix
851  Double_t ele[4] = {mxx , mxy , mxy , myy };
852  TMatrixDSym m0(2,ele);
853 
854  // Find eigenvectors
855  TMatrixDSymEigen m(m0);
856  TVectorD eval(2);
857  TMatrixD evecm = m.GetEigenVectors();
858  eval = m.GetEigenValues();
859  // Largest eigenvector
860  int jev = 0;
861  // cout<<eval[0]<<" "<<eval[1]<<endl;
862  if (eval[0] < eval[1]) jev = 1;
863  TVectorD evec0(2);
864  // Principle axis
865  evec0 = TMatrixDColumn(evecm, jev);
866  Double_t compx=evec0[0];
867  Double_t compy=evec0[1];
868  TVector2 evec(compx, compy);
869  Double_t sig=0;
870  if(jev==1) sig=TMath::Sqrt(TMath::Abs(eval[0])/sump2);
871  if(jev==0) sig=TMath::Sqrt(TMath::Abs(eval[1])/sump2);
872 
873  return sig;
874 
875 }
876 
877 //________________________________________________________________________
879  //calc subtracted jet mass
880 
881  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
884  else
885  return Sigma2(jet, jetContNb);
886 
887 }
888 
889 
890 //_________________________________________________________________________
892 
893  AliJetContainer *jetCont = GetJetContainer(jetContNb);
894  if (!jet->GetNumberOfTracks())
895  return;
896 
897  Double_t ptJet = jet->Pt();
898 
899  AliVParticle *vp1 = 0x0;
900  std::vector<int> ordindex;
901  ordindex=jet->GetPtSortedTrackConstituentIndexes(jetCont->GetParticleContainer()->GetArray());
902  //Printf("Nbof const = %d", jet->GetNumberOfTracks());
903  //Printf("ordindex[0] = %d, ordindex[1] = %d", ordindex[0], ordindex[1]);
904  //if(ordindex.size()<2) return -1;
905 
906  for(Int_t iconst =0; iconst<jet->GetNumberOfTracks(); iconst++){
907 
908  vp1 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[iconst], jetCont->GetParticleContainer()->GetArray()));
909  if (!vp1){
910  Printf("AliVParticle associated to Leading constituent not found");
911  return;
912  }
913 
914  if (nTFractions[0] <= 0.7*ptJet){
915  nTFractions[0] += vp1->Pt();
916  nTFractions[1] +=1;
917  }
918 
919  if (nTFractions[2] <= 0.8*ptJet){
920  nTFractions[2] += vp1->Pt();
921  nTFractions[3] +=1;
922  }
923 
924  if (nTFractions[4] <= 0.9*ptJet){
925  nTFractions[4] += vp1->Pt();
926  nTFractions[5] +=1;
927  }
928 
929  if (nTFractions[6] <= 0.95*ptJet){
930  nTFractions[6] += vp1->Pt();
931  nTFractions[7] +=1;
932  }
933  }
934 }
935 //_________________________________________________________________________________________________
937 
938  //WARNING!!! Only works for parent jets that are clustered with Anti-Kt! To change go to AliEmcalJetFinder.cxx and look at the Nsubjettiness() function
939 
940  //Algorithm==0 -> kt_axes;
941  // Algorithm==1 -> ca_axes;
942  //Algorithm==2 -> antikt_0p2_axes;
943  //Algorithm==3 -> wta_kt_axes;
944  //Algorithm==4 -> wta_ca_axes;
945  //Algorithm==5 -> onepass_kt_axes;
946  //Algorithm==6 -> onepass_ca_axes;
947  //Algorithm==7 -> onepass_antikt_0p2_axes;
948  //Algorithm==8 -> onepass_wta_kt_axes;
949  //Algorithm==9 -> onepass_wta_ca_axes;
950  //Algorithm==10 -> min_axes;
951 
952 
953  //Option==0 returns Nsubjettiness Value
954  //Option==1 && N==2 returns opening angle between two subjet axes(Delta R?)
955  //Option==2 && N==2 returns Delta R
956 
957  if (Jet->GetNumberOfTracks()>=N){
958  AliJetContainer *JetCont = GetJetContainer(JetContNb);
959  AliEmcalJetFinder *JetFinder=new AliEmcalJetFinder("Nsubjettiness");
960  JetFinder->SetJetMaxEta(0.9-fJetRadius);
961  JetFinder->SetRadius(fJetRadius);
962  JetFinder->SetJetAlgorithm(0); //0 for anti-kt 1 for kt //this is for the JET!!!!!!!!!! Not the SubJets
963  JetFinder->SetRecombSheme(0);
964  JetFinder->SetJetMinPt(Jet->Pt());
965  //Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
966  Double_t dVtx[3]={1,1,1};
967  //Printf("JetFinder->Nsubjettiness =%f", JetFinder->Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,fSubjetRadius,Beta,Option));
968  return JetFinder->Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,0.2,Beta,Option,0,0,0);
969 
970  }
971  else return -2;
972 }
973 
974 
975 
976 //________________________________________________________________________
978 
980  TClonesArray *tracksArray = partCont->GetArray();
981 
982  if(!partCont || !tracksArray) return -99999;
983  AliAODTrack *track = 0x0;
984  AliEmcalParticle *emcPart = 0x0;
985 
986 
987  TList *trackList = new TList();
988  Int_t triggers[100];
989  for (Int_t iTrigger=0; iTrigger<100; iTrigger++) triggers[iTrigger] = 0;
990  Int_t iTT = 0;
991 
992  for(Int_t iTrack=0; iTrack <= tracksArray->GetEntriesFast(); iTrack++){
993 
994 
995  if (fJetShapeSub == kConstSub){
996  emcPart = static_cast<AliEmcalParticle*>(tracksArray->At(iTrack));
997  if (!emcPart) continue;
998  if(TMath::Abs(emcPart->Eta())>0.9) continue;
999  if (emcPart->Pt()<0.15) continue;
1000 
1001  if ((emcPart->Pt() >= minpT) && (emcPart->Pt()< maxpT)) {
1002  trackList->Add(emcPart);
1003  triggers[iTT] = iTrack;
1004  iTT++;
1005  }
1006  }
1007  else{
1008  track = static_cast<AliAODTrack*>(tracksArray->At(iTrack));
1009  if (!track) continue;
1010  if(TMath::Abs(track->Eta())>0.9) continue;
1011  if (track->Pt()<0.15) continue;
1012  if (!(track->TestFilterBit(768))) continue;
1013 
1014  if ((track->Pt() >= minpT) && (track->Pt()< maxpT)) {
1015  trackList->Add(track);
1016  triggers[iTT] = iTrack;
1017  iTT++;
1018 
1019  }
1020  }
1021  }
1022 
1023  if (iTT == 0) return -99999;
1024  Int_t nbRn = 0, index = 0 ;
1025  TRandom3* random = new TRandom3(0);
1026  nbRn = random->Integer(iTT);
1027 
1028  index = triggers[nbRn];
1029  //Printf("iTT Total= %d, nbRn = %d, Index = %d",iTT, nbRn, index );
1030  return index;
1031 
1032 }
1033 
1034 //__________________________________________________________________________________
1036 
1037  if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1038  else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1039  if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1040  else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1041  double dphi = mphi-vphi;
1042  if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1043  else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1044  return dphi;//dphi in [-Pi, Pi]
1045 }
1046 
1047 
1048 //________________________________________________________________________
1050  //
1051  // retrieve event objects
1052  //
1054  return kFALSE;
1055 
1056  return kTRUE;
1057 }
1058 
1059 //_______________________________________________________________________
1061 {
1062  // Called once at the end of the analysis.
1063 
1064  AliInfo("Terminate");
1065  AliAnalysisTaskSE::Terminate();
1066 
1067  fOutput = dynamic_cast<AliEmcalList*> (GetOutputData(1));
1068  if (!fOutput) {
1069  AliError("fOutput not available");
1070  return;
1071  }
1072 
1073  fTreeObservableTagging = dynamic_cast<TTree*>(GetOutputData(2));
1074  if (!fTreeObservableTagging){
1075  Printf("ERROR: fTreeObservableTagging not available");
1076  return;
1077  }
1078 
1079 }
1080 
1081 //_________________________________________________________________________
1082 void AliAnalysisTaskEmcalJetShapesMC::SoftDrop(AliEmcalJet *fJet,AliJetContainer *fJetCont, double zcut, double beta, int ReclusterAlgo){
1083 
1084  std::vector<fastjet::PseudoJet> fInputVectors;
1085  Double_t JetInvMass=0, PseudJetInvMass=0, TrackMom = 0, TrackEnergy = 0;
1086 
1087  AliParticleContainer *fTrackCont = fJetCont->GetParticleContainer();
1088 
1089  Double_t JetEta=fJet->Eta(),JetPhi=fJet->Phi();
1090  Double_t FJTrackEta[9999],FJTrackPhi[9999],FJTrackPt[9999],EmcalJetTrackEta[9999],EmcalJetTrackPhi[9999],EmcalJetTrackPt[9999];
1091  UShort_t FJNTracks=0,EmcalJetNTracks=0;
1092 
1093  if (fTrackCont) for (Int_t i=0; i<fJet->GetNumberOfTracks(); i++) {
1094  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
1095  if (!fTrk) continue;
1096  JetInvMass += fTrk->M();
1097 
1098  fastjet::PseudoJet PseudoTracks(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
1099  TrackMom += TMath::Sqrt(TMath::Power(fTrk->Px(),2)+TMath::Power(fTrk->Py(),2)+TMath::Power(fTrk->Pz(),2));
1100  TrackEnergy += fTrk->E();
1101  PseudoTracks.set_user_index(fJet->TrackAt(i)+100);
1102  PseudJetInvMass += PseudoTracks.m();
1103  fInputVectors.push_back(PseudoTracks);
1104  EmcalJetTrackEta[i]=fTrk->Eta();
1105  EmcalJetTrackPhi[i]=fTrk->Phi();
1106  EmcalJetTrackPt[i]=fTrk->Pt();
1107  EmcalJetNTracks++;
1108  }
1109  fastjet::JetDefinition *fJetDef;
1110  fastjet::ClusterSequence *fClustSeqSA;
1111 
1112 
1113 
1114  fJetDef = new fastjet::JetDefinition(fastjet::antikt_algorithm, fJetRadius*2, static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
1115 
1116  try {
1117  fClustSeqSA = new fastjet::ClusterSequence(fInputVectors, *fJetDef);
1118  } catch (fastjet::Error) {
1119  AliError(" [w] FJ Exception caught.");
1120  //return -1;
1121  }
1122 
1123  std::vector<fastjet::PseudoJet> fOutputJets;
1124  fOutputJets.clear();
1125  fOutputJets=fClustSeqSA->inclusive_jets(0);
1126 
1127 
1128  std::vector<fastjet::PseudoJet> jet_constituents = fOutputJets[0].constituents();
1129  fastjet::contrib::SoftDrop softdrop(beta, zcut);
1130  //fastjet::contrib::SoftDrop softdrop_antikt(beta,zcut);
1131  softdrop.set_verbose_structure(kTRUE);
1132  //fastjet::JetDefinition jet_def_akt(fastjet::antikt_algorithm, 0.4);
1133  // fastjet::contrib::Recluster *antiKT_Recluster(jet_def_akt);
1134  fastjet::contrib::Recluster *recluster;
1135  if(ReclusterAlgo == 1) recluster = new fastjet::contrib::Recluster(fastjet::kt_algorithm,1,true);
1136  if(ReclusterAlgo == 2) recluster = new fastjet::contrib::Recluster(fastjet::antikt_algorithm,1,true);
1137  if(ReclusterAlgo == 0) recluster = new fastjet::contrib::Recluster(fastjet::cambridge_algorithm,1,true);
1138  softdrop.set_reclustering(true,recluster);
1139  fastjet::PseudoJet finaljet = softdrop(fOutputJets[0]);
1140  // fastjet::PseudoJet finaljet_antikt = softdrop_antikt(fOutputJets[0]);
1141  //cout<< finaljet.structure_of<fastjet::contrib::SoftDrop>().symmetry()<<endl;
1142  //cout<< finaljet_antikt.structure_of<fastjet::contrib::SoftDrop>().symmetry()<<endl;
1143 
1144 
1145  //AliEmcalJet* jet = new AliEmcalJet(finaljet.perp(), finaljet.eta(), finaljet.phi(), finaljet.m());
1146  //std::vector<fastjet::PseudoJet> fSDTracks=finaljet.constituents();
1147  //Double_t FastjetTrackDelR,EmcalTrackDelR;
1148  //for(Int_t i=0;i<fJet->GetNumberOfConstituents();i++){
1149  // if(i<=finaljet.constituents().size()){
1150  // FastjetTrackDelR = TMath::Sqrt(TMath::Power(fSDTracks[i].eta()-JetEta,2)+TMath::Power(fSDTracks[i].phi()-JetPhi,2));
1151  // FJTrackEta[i]=fSDTracks[i].eta();
1152  // FJTrackPhi[i]=fSDTracks[i].phi();
1153  // FJTrackPt[i]=fSDTracks[i].perp();
1154  // FJNTracks++;
1155  // }
1156  // AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
1157  // EmcalTrackDelR = TMath::Sqrt(TMath::Power(fTrk->Eta()-JetEta,2)+TMath::Power(fTrk->Phi()-JetPhi,2));
1158  //}
1159  Int_t NDroppedTracks = fJet->GetNumberOfTracks()-finaljet.constituents().size();
1160  //Int_t nConstituents(fClustSeqSA->constituents(finaljet).size());
1161  //jet->SetNumberOfTracks(nConstituents);
1162  Double_t SymParam, Mu, DeltaR, GroomedPt;
1163  Int_t NGroomedBranches;
1164  SymParam=(finaljet.structure_of<fastjet::contrib::SoftDrop>().symmetry());
1165  Mu=(finaljet.structure_of<fastjet::contrib::SoftDrop>().mu());
1166  DeltaR=(finaljet.structure_of<fastjet::contrib::SoftDrop>().delta_R());
1167  NGroomedBranches=finaljet.structure_of<fastjet::contrib::SoftDrop>().dropped_count();
1168  GroomedPt=finaljet.perp();
1169  if(beta==0){
1170  if(ReclusterAlgo==0){
1171  fShapesVar[12]=SymParam;
1172  fShapesVar[13]=DeltaR;
1173  fShapesVar[14]=GroomedPt;
1174  fShapesVar[15]=NGroomedBranches;}
1175  if(ReclusterAlgo==1){
1176  fShapesVar[16]=SymParam;
1177  fShapesVar[17]=DeltaR;
1178  fShapesVar[18]=GroomedPt;
1179  fShapesVar[19]=NGroomedBranches;}
1180 
1181  if(ReclusterAlgo==2){
1182  fShapesVar[20]=SymParam;
1183  fShapesVar[21]=DeltaR;
1184  fShapesVar[22]=GroomedPt;
1185  fShapesVar[23]=NGroomedBranches;}}
1186  if(beta==1){
1187  fShapesVar[24]=SymParam;
1188  fShapesVar[25]=DeltaR;
1189  fShapesVar[26]=GroomedPt;
1190  fShapesVar[27]=NGroomedBranches;}
1191 
1192  if(beta==2){
1193  fShapesVar[28]=SymParam;
1194  fShapesVar[29]=DeltaR;
1195  fShapesVar[30]=GroomedPt;
1196  fShapesVar[31]=NGroomedBranches;}
1197 
1198 
1199  return;
1200 
1201 
1202 }
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:114
ClassImp(AliAnalysisTaskEmcalJetShapesMC) AliAnalysisTaskEmcalJetShapesMC
Double_t Py() const
Definition: AliEmcalJet.h:100
Double_t Phi() const
Definition: AliEmcalJet.h:110
Float_t GetPythiaEventWeight() const
Float_t GetPartonPhi7() const
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)
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:112
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:153
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:132
Double_t Px() const
Definition: AliEmcalJet.h:99
Float_t GetJetCoreFrac(AliEmcalJet *jet, Int_t jetContNb=0)
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
Double_t fjNSubJettiness(AliEmcalJet *Jet, Int_t JetContNb, Int_t N, Int_t Algorithm, Double_t Beta, Int_t Option)
Int_t GetPartonFlag6() const
AliParticleContainer * GetParticleContainer() const
Double_t GetFirstOrderSubtractedConstituent() const
void SetRecombSheme(Int_t val)
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
Task to store and correlate the MC shapes.
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)
void SoftDrop(AliEmcalJet *fJet, AliJetContainer *fJetCont, double zcut, double beta, Int_t ReclusterAlgo)
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()
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)
Double_t Pt() const
Definition: AliEmcalJet.h:102
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:44
const AliEmcalPythiaInfo * GetPythiaInfo() const
unsigned short UShort_t
Definition: External.C:28
Double_t Pz() const
Definition: AliEmcalJet.h:101
Float_t Circularity(AliEmcalJet *jet, Int_t jetContNb=0)
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:256
Float_t GetPartonPt7() const
Float_t LeSub(AliEmcalJet *jet, Int_t jetContNb=0)
Double_t M() const
Definition: AliEmcalJet.h:113
Container for jet within the EMCAL jet framework.
std::vector< int > GetPtSortedTrackConstituentIndexes(TClonesArray *tracks) const
Float_t GetPartonPt6() const