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