AliPhysics  648edd6 (648edd6)
 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] =fjNSubJettiness(jet1,0,1,10,1,0);
419  fShapesVar[9] =fjNSubJettiness(jet1,0,2,10,1,0);
420  //nsub1 and nsub2 for akt
421  fShapesVar[10] = 0;
422  //fjNSubJettiness(jet1,0,1,10,1,0);
423  fShapesVar[11] =0;
424  //fjNSubJettiness(jet1,0,2,10,1,0);
425 
426  //SoftDropParameters for different reclustering strategies and beta values
427  SoftDrop(jet1,jetCont,0.1,0,0);
428  SoftDrop(jet1,jetCont,0.1,0,1);
429  SoftDrop(jet1,jetCont,0.1,0,2);
430  SoftDrop(jet1,jetCont,0.1,1,0);
431  SoftDrop(jet1,jetCont,0.1,2,0);
432 
433 
434  // Float_t nTFractions[8]={0.,0.,0.,0.,0.,0.,0.,0.};
435  //NTValues(jet1, 0, nTFractions);
436  //shape 13 is pythia weight!
437  //for (Int_t ishape=14; ishape<22; ishape++) fShapesVar[ishape] = nTFractions[ishape-14];
438 
439  //fShapesVar[22]= GetSubjetFraction(jet1,0,fJetRadius,Reclusterer1);
440 
441  fTreeObservableTagging->Fill();
442 
443  }
444 
445  }
446 
447  return kTRUE;
448 }
449 
450 //________________________________________________________________________
452  //calc subtracted jet mass
453  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
455  else return jet->GetShapeProperties()->GetSecondOrderSubtracted();
456  else
457  return jet->M();
458 }
459 
460 //________________________________________________________________________
462 
463  AliJetContainer *jetCont = GetJetContainer(jetContNb);
464  if (!jet->GetNumberOfTracks())
465  return 0;
466  Double_t den=0.;
467  Double_t num = 0.;
468  AliVParticle *vp1 = 0x0;
469  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
470  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
471 
472  if (!vp1){
473  Printf("AliVParticle associated to constituent not found");
474  continue;
475  }
476 
477  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
478  Double_t dr2 = (vp1->Eta()-jet->Eta())*(vp1->Eta()-jet->Eta()) + dphi*dphi;
479  Double_t dr = TMath::Sqrt(dr2);
480  num=num+vp1->Pt()*dr;
481  den=den+vp1->Pt();
482  }
483  return num/den;
484 }
485 
486 //________________________________________________________________________
488 
489  if((fJetShapeSub==kDerivSub) && (jetContNb==0))
492  else
493  return Angularity(jet, jetContNb);
494 
495 }
496 
497 
498 //________________________________________________________________________
500 
501  AliJetContainer *jetCont = GetJetContainer(jetContNb);
502  if (!jet->GetNumberOfTracks())
503  return 0;
504  Double_t den=0.;
505  Double_t num = 0.;
506  AliVParticle *vp1 = 0x0;
507  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
508  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
509 
510  if (!vp1){
511  Printf("AliVParticle associated to constituent not found");
512  continue;
513  }
514 
515  num=num+vp1->Pt()*vp1->Pt();
516  den=den+vp1->Pt();
517  }
518  return TMath::Sqrt(num)/den;
519 }
520 
521 //________________________________________________________________________
523  //calc subtracted jet mass
524  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
526  else return jet->GetShapeProperties()->GetSecondOrderSubtractedpTD();
527  else
528  return PTD(jet, jetContNb);
529 
530 }
531 
532 //_____________________________________________________________________________
534 
535  AliJetContainer *jetCont = GetJetContainer(jetContNb);
536  if (!jet->GetNumberOfTracks())
537  return 0;
538  Double_t mxx = 0.;
539  Double_t myy = 0.;
540  Double_t mxy = 0.;
541  int nc = 0;
542  Double_t sump2 = 0.;
543  Double_t pxjet=jet->Px();
544  Double_t pyjet=jet->Py();
545  Double_t pzjet=jet->Pz();
546 
547 
548  //2 general normalized vectors perpendicular to the jet
549  TVector3 ppJ1(pxjet, pyjet, pzjet);
550  TVector3 ppJ3(- pxjet* pzjet, - pyjet * pzjet, pxjet * pxjet + pyjet * pyjet);
551  ppJ3.SetMag(1.);
552  TVector3 ppJ2(-pyjet, pxjet, 0);
553  ppJ2.SetMag(1.);
554  AliVParticle *vp1 = 0x0;
555  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
556  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
557 
558  if (!vp1){
559  Printf("AliVParticle associated to constituent not found");
560  continue;
561  }
562 
563  TVector3 pp(vp1->Px(), vp1->Py(), vp1->Pz());
564 
565  //local frame
566  TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
567  TVector3 pPerp = pp - pLong;
568  //projection onto the two perpendicular vectors defined above
569 
570  Float_t ppjX = pPerp.Dot(ppJ2);
571  Float_t ppjY = pPerp.Dot(ppJ3);
572  Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
573  if(ppjT<=0) return 0;
574 
575  mxx += (ppjX * ppjX / ppjT);
576  myy += (ppjY * ppjY / ppjT);
577  mxy += (ppjX * ppjY / ppjT);
578  nc++;
579  sump2 += ppjT;}
580 
581  if(nc<2) return 0;
582  if(sump2==0) return 0;
583  // Sphericity Matrix
584  Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
585  TMatrixDSym m0(2,ele);
586 
587  // Find eigenvectors
588  TMatrixDSymEigen m(m0);
589  TVectorD eval(2);
590  TMatrixD evecm = m.GetEigenVectors();
591  eval = m.GetEigenValues();
592  // Largest eigenvector
593  int jev = 0;
594  // cout<<eval[0]<<" "<<eval[1]<<endl;
595  if (eval[0] < eval[1]) jev = 1;
596  TVectorD evec0(2);
597  // Principle axis
598  evec0 = TMatrixDColumn(evecm, jev);
599  Double_t compx=evec0[0];
600  Double_t compy=evec0[1];
601  TVector2 evec(compx, compy);
602  Double_t circ=0;
603  if(jev==1) circ=2*eval[0];
604  if(jev==0) circ=2*eval[1];
605 
606  return circ;
607 
608 
609 
610 }
611 
612 
613 
614 
615 //________________________________________________________________________
617  //calc subtracted jet mass
618 
619  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
622  else
623  return Circularity(jet, jetContNb);
624 
625 }
626 
627 //________________________________________________________________________
629 
630  AliJetContainer *jetCont = GetJetContainer(jetContNb);
631  if (!jet->GetNumberOfTracks())
632  return 0;
633  Double_t den=0.;
634  Double_t num = 0.;
635  AliVParticle *vp1 = 0x0;
636  AliVParticle *vp2 = 0x0;
637  std::vector<int> ordindex;
638  ordindex=jet->GetPtSortedTrackConstituentIndexes(jetCont->GetParticleContainer()->GetArray());
639  //Printf("Nbof const = %d", jet->GetNumberOfTracks());
640  //Printf("ordindex[0] = %d, ordindex[1] = %d", ordindex[0], ordindex[1]);
641 
642  if(ordindex.size()<2) return -1;
643 
644  vp1 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[0], jetCont->GetParticleContainer()->GetArray()));
645  if (!vp1){
646  Printf("AliVParticle associated to Leading constituent not found");
647  return -1;
648  }
649 
650  vp2 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[1], jetCont->GetParticleContainer()->GetArray()));
651  if (!vp2){
652  Printf("AliVParticle associated to Subleading constituent not found");
653  return -1;
654  }
655 
656 
657  num=vp1->Pt();
658  den=vp2->Pt();
659  //Printf("vp1->Pt() =%f, vp2->Pt() =%f", vp1->Pt(), vp2->Pt());
660 
661 return num-den;
662 }
663 
664 //________________________________________________________________________
666  //calc subtracted jet mass
667 
668  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
671  else
672  return LeSub(jet, jetContNb);
673 
674 }
675 
676 //________________________________________________________________________
678  //calc subtracted jet mass
679 
680  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
683  else
684  return jet->GetNumberOfTracks();
685 
686  }
687 
688 
689 //________________________________________________________________________
690 AliEmcalJetFinder *AliAnalysisTaskEmcalJetShapesMC::Recluster(AliEmcalJet *Jet, Int_t JetContNb, Double_t JetRadius, Double_t SubJetRadius, Double_t SubJetMinPt, Int_t Algorithm, const char* Name){
691 
692  AliJetContainer *JetCont = GetJetContainer(JetContNb);
693  AliEmcalJetFinder *Reclusterer = new AliEmcalJetFinder(Name); //JetFinder Object for reclustered jets
694  Reclusterer->SetRadius(SubJetRadius);
695  Reclusterer->SetJetMinPt(SubJetMinPt);
696  Reclusterer->SetJetAlgorithm(Algorithm); //0 for anti-kt 1 for kt
697  Reclusterer->SetJetMaxEta(0.9-JetRadius);
698  Reclusterer->SetRecombSheme(0);
699  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
700 
701  //Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
702  Double_t dVtx[3]={0.,0.,0.};
703  if(Reclusterer->AliEmcalJetFinder::Filter(Jet, JetCont, dVtx)){;} //reclustering jet1 using the jetfinderobject Reclusterer
704  return Reclusterer;
705 }
706 
707 
708 
709 
710 //----------------------------------------------------------------------
712  AliJetContainer *JetCont = GetJetContainer(JetContNb);
713  AliEmcalJet *SubJet=NULL;
714  Double_t DeltaR1=0;
715  Double_t DeltaR2=0;
716  AliVParticle *JetParticle=0x0;
717  Double_t SubJetiness_Numerator = 0;
718  Double_t SubJetiness_Denominator = 0;
719  Double_t Index=-2;
720  if (Reclusterer->GetNumberOfJets() < 1) return -2;
721  Index=SubJetOrdering(Jet,Reclusterer,1,0,kTRUE);
722  if(Index==-999) return -2;
723  SubJetiness_Numerator=(Reclusterer->GetJet(Index)->Pt());
724  SubJetiness_Denominator=Jet->Pt();
725  return SubJetiness_Numerator/SubJetiness_Denominator;
726 
727 
728 }
729 //__________________________________________________________________________________
731 
732  AliJetContainer *jetCont = GetJetContainer(jetContNb);
733  if (!jet->GetNumberOfTracks())
734  return 0;
735  Double_t den=0.;
736  Double_t num = 0.;
737  AliVParticle *vp1 = 0x0;
738  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
739  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
740 
741  if (!vp1){
742  Printf("AliVParticle associated to constituent not found");
743  continue;
744  }
745 
746  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
747  Double_t dr2 = (vp1->Eta()-jet->Eta())*(vp1->Eta()-jet->Eta()) + dphi*dphi;
748  Double_t dr = TMath::Sqrt(dr2);
749  if(dr<=fSubjetRadius) num=num+vp1->Pt();
750 
751  }
752  return num/jet->Pt();
753 }
754 
755 
756 
757 
758 //________________________________________________________________________
760  //calc subtracted jet mass
761 
762  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
765  else
766  return CoreFrac(jet, jetContNb);
767 
768 }
769 
770 
771 
772 
773 //----------------------------------------------------------------------
775  AliEmcalJet *SubJet=NULL;
776  Double_t SortingVariable;
777  Int_t ArraySize =N+1;
778  TArrayD *JetSorter = new TArrayD(ArraySize);
779  TArrayD *JetIndexSorter = new TArrayD(ArraySize);
780  for (Int_t i=0; i<ArraySize; i++){
781  JetSorter->SetAt(0,i);
782  }
783  for (Int_t i=0; i<ArraySize; i++){
784  JetIndexSorter->SetAt(0,i);
785  }
786  if(Reclusterer->GetNumberOfJets()<N) return -999;
787  for (Int_t i=0; i<Reclusterer->GetNumberOfJets(); i++){
788  SubJet=Reclusterer->GetJet(i);
789  if (Type==0) SortingVariable=SubJet->Pt();
790  else if (Type==1) SortingVariable=SubJet->E();
791  else if (Type==2) SortingVariable=SubJet->M();
792  for (Int_t j=0; j<N; j++){
793  if (SortingVariable>JetSorter->GetAt(j)){
794  for (Int_t k=N-1; k>=j; k--){
795  JetSorter->SetAt(JetSorter->GetAt(k),k+1);
796  JetIndexSorter->SetAt(JetIndexSorter->GetAt(k),k+1);
797  }
798  JetSorter->SetAt(SortingVariable,j);
799  JetIndexSorter->SetAt(i,j);
800  break;
801  }
802  }
803  }
804  if (!Index) return JetSorter->GetAt(N-1);
805  else return JetIndexSorter->GetAt(N-1);
806 }
807 
808 
809 
810 //returns -1 if the Nth hardest jet is requested where N>number of available jets
811 //type: 0=Pt 1=E 2=M
812 //Index TRUE=returns index FALSE=returns value of quantatiy in question
813 
814 
815 
816 
817 
818 //______________________________________________________________________________
820 
821  AliJetContainer *jetCont = GetJetContainer(jetContNb);
822  if (!jet->GetNumberOfTracks())
823  return 0;
824  Double_t mxx = 0.;
825  Double_t myy = 0.;
826  Double_t mxy = 0.;
827  int nc = 0;
828  Double_t sump2 = 0.;
829 
830  AliVParticle *vp1 = 0x0;
831  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
832  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
833 
834  if (!vp1){
835  Printf("AliVParticle associated to constituent not found");
836  continue;
837  }
838 
839  Double_t ppt=vp1->Pt();
840  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
841 
842  Double_t deta = vp1->Eta()-jet->Eta();
843  mxx += ppt*ppt*deta*deta;
844  myy += ppt*ppt*dphi*dphi;
845  mxy -= ppt*ppt*deta*TMath::Abs(dphi);
846  nc++;
847  sump2 += ppt*ppt;
848 
849  }
850  if(nc<2) return 0;
851  if(sump2==0) return 0;
852  // Sphericity Matrix
853  Double_t ele[4] = {mxx , mxy , mxy , myy };
854  TMatrixDSym m0(2,ele);
855 
856  // Find eigenvectors
857  TMatrixDSymEigen m(m0);
858  TVectorD eval(2);
859  TMatrixD evecm = m.GetEigenVectors();
860  eval = m.GetEigenValues();
861  // Largest eigenvector
862  int jev = 0;
863  // cout<<eval[0]<<" "<<eval[1]<<endl;
864  if (eval[0] < eval[1]) jev = 1;
865  TVectorD evec0(2);
866  // Principle axis
867  evec0 = TMatrixDColumn(evecm, jev);
868  Double_t compx=evec0[0];
869  Double_t compy=evec0[1];
870  TVector2 evec(compx, compy);
871  Double_t sig=0;
872  if(jev==1) sig=TMath::Sqrt(TMath::Abs(eval[0])/sump2);
873  if(jev==0) sig=TMath::Sqrt(TMath::Abs(eval[1])/sump2);
874 
875  return sig;
876 
877 }
878 
879 //________________________________________________________________________
881  //calc subtracted jet mass
882 
883  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
886  else
887  return Sigma2(jet, jetContNb);
888 
889 }
890 
891 
892 //_________________________________________________________________________
894 
895  AliJetContainer *jetCont = GetJetContainer(jetContNb);
896  if (!jet->GetNumberOfTracks())
897  return;
898 
899  Double_t ptJet = jet->Pt();
900 
901  AliVParticle *vp1 = 0x0;
902  std::vector<int> ordindex;
903  ordindex=jet->GetPtSortedTrackConstituentIndexes(jetCont->GetParticleContainer()->GetArray());
904  //Printf("Nbof const = %d", jet->GetNumberOfTracks());
905  //Printf("ordindex[0] = %d, ordindex[1] = %d", ordindex[0], ordindex[1]);
906  //if(ordindex.size()<2) return -1;
907 
908  for(Int_t iconst =0; iconst<jet->GetNumberOfTracks(); iconst++){
909 
910  vp1 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[iconst], jetCont->GetParticleContainer()->GetArray()));
911  if (!vp1){
912  Printf("AliVParticle associated to Leading constituent not found");
913  return;
914  }
915 
916  if (nTFractions[0] <= 0.7*ptJet){
917  nTFractions[0] += vp1->Pt();
918  nTFractions[1] +=1;
919  }
920 
921  if (nTFractions[2] <= 0.8*ptJet){
922  nTFractions[2] += vp1->Pt();
923  nTFractions[3] +=1;
924  }
925 
926  if (nTFractions[4] <= 0.9*ptJet){
927  nTFractions[4] += vp1->Pt();
928  nTFractions[5] +=1;
929  }
930 
931  if (nTFractions[6] <= 0.95*ptJet){
932  nTFractions[6] += vp1->Pt();
933  nTFractions[7] +=1;
934  }
935  }
936 }
937 //_________________________________________________________________________________________________
939 
940  //WARNING!!! Only works for parent jets that are clustered with Anti-Kt! To change go to AliEmcalJetFinder.cxx and look at the Nsubjettiness() function
941 
942  //Algorithm==0 -> kt_axes;
943  // Algorithm==1 -> ca_axes;
944  //Algorithm==2 -> antikt_0p2_axes;
945  //Algorithm==3 -> wta_kt_axes;
946  //Algorithm==4 -> wta_ca_axes;
947  //Algorithm==5 -> onepass_kt_axes;
948  //Algorithm==6 -> onepass_ca_axes;
949  //Algorithm==7 -> onepass_antikt_0p2_axes;
950  //Algorithm==8 -> onepass_wta_kt_axes;
951  //Algorithm==9 -> onepass_wta_ca_axes;
952  //Algorithm==10 -> min_axes;
953 
954 
955  //Option==0 returns Nsubjettiness Value
956  //Option==1 && N==2 returns opening angle between two subjet axes(Delta R?)
957  //Option==2 && N==2 returns Delta R
958 
959  if (Jet->GetNumberOfTracks()>=N){
960  AliJetContainer *JetCont = GetJetContainer(JetContNb);
961  AliEmcalJetFinder *JetFinder=new AliEmcalJetFinder("Nsubjettiness");
962  JetFinder->SetJetMaxEta(0.9-fJetRadius);
963  JetFinder->SetRadius(fJetRadius);
964  JetFinder->SetJetAlgorithm(0); //0 for anti-kt 1 for kt //this is for the JET!!!!!!!!!! Not the SubJets
965  JetFinder->SetRecombSheme(0);
966  JetFinder->SetJetMinPt(Jet->Pt());
967  //Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
968  Double_t dVtx[3]={1,1,1};
969  //Printf("JetFinder->Nsubjettiness =%f", JetFinder->Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,fSubjetRadius,Beta,Option));
970  return JetFinder->Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,0.2,Beta,Option,0,0,0);
971 
972  }
973  else return -2;
974 }
975 
976 
977 
978 //________________________________________________________________________
980 
982  TClonesArray *tracksArray = partCont->GetArray();
983 
984  if(!partCont || !tracksArray) return -99999;
985  AliAODTrack *track = 0x0;
986  AliEmcalParticle *emcPart = 0x0;
987 
988 
989  TList *trackList = new TList();
990  Int_t triggers[100];
991  for (Int_t iTrigger=0; iTrigger<100; iTrigger++) triggers[iTrigger] = 0;
992  Int_t iTT = 0;
993 
994  for(Int_t iTrack=0; iTrack <= tracksArray->GetEntriesFast(); iTrack++){
995 
996 
997  if (fJetShapeSub == kConstSub){
998  emcPart = static_cast<AliEmcalParticle*>(tracksArray->At(iTrack));
999  if (!emcPart) continue;
1000  if(TMath::Abs(emcPart->Eta())>0.9) continue;
1001  if (emcPart->Pt()<0.15) continue;
1002 
1003  if ((emcPart->Pt() >= minpT) && (emcPart->Pt()< maxpT)) {
1004  trackList->Add(emcPart);
1005  triggers[iTT] = iTrack;
1006  iTT++;
1007  }
1008  }
1009  else{
1010  track = static_cast<AliAODTrack*>(tracksArray->At(iTrack));
1011  if (!track) continue;
1012  if(TMath::Abs(track->Eta())>0.9) continue;
1013  if (track->Pt()<0.15) continue;
1014  if (!(track->TestFilterBit(768))) continue;
1015 
1016  if ((track->Pt() >= minpT) && (track->Pt()< maxpT)) {
1017  trackList->Add(track);
1018  triggers[iTT] = iTrack;
1019  iTT++;
1020 
1021  }
1022  }
1023  }
1024 
1025  if (iTT == 0) return -99999;
1026  Int_t nbRn = 0, index = 0 ;
1027  TRandom3* random = new TRandom3(0);
1028  nbRn = random->Integer(iTT);
1029 
1030  index = triggers[nbRn];
1031  //Printf("iTT Total= %d, nbRn = %d, Index = %d",iTT, nbRn, index );
1032  return index;
1033 
1034 }
1035 
1036 //__________________________________________________________________________________
1038 
1039  if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1040  else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1041  if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1042  else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1043  double dphi = mphi-vphi;
1044  if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1045  else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1046  return dphi;//dphi in [-Pi, Pi]
1047 }
1048 
1049 
1050 //________________________________________________________________________
1052  //
1053  // retrieve event objects
1054  //
1056  return kFALSE;
1057 
1058  return kTRUE;
1059 }
1060 
1061 //_______________________________________________________________________
1063 {
1064  // Called once at the end of the analysis.
1065 
1066  AliInfo("Terminate");
1067  AliAnalysisTaskSE::Terminate();
1068 
1069  fOutput = dynamic_cast<AliEmcalList*> (GetOutputData(1));
1070  if (!fOutput) {
1071  AliError("fOutput not available");
1072  return;
1073  }
1074 
1075  fTreeObservableTagging = dynamic_cast<TTree*>(GetOutputData(2));
1076  if (!fTreeObservableTagging){
1077  Printf("ERROR: fTreeObservableTagging not available");
1078  return;
1079  }
1080 
1081 }
1082 
1083 //_________________________________________________________________________
1084 void AliAnalysisTaskEmcalJetShapesMC::SoftDrop(AliEmcalJet *fJet,AliJetContainer *fJetCont, double zcut, double beta, int ReclusterAlgo){
1085 
1086  std::vector<fastjet::PseudoJet> fInputVectors;
1087  Double_t JetInvMass=0, PseudJetInvMass=0, TrackMom = 0, TrackEnergy = 0;
1088 
1089  AliParticleContainer *fTrackCont = fJetCont->GetParticleContainer();
1090 
1091  Double_t JetEta=fJet->Eta(),JetPhi=fJet->Phi();
1092  Double_t FJTrackEta[9999],FJTrackPhi[9999],FJTrackPt[9999],EmcalJetTrackEta[9999],EmcalJetTrackPhi[9999],EmcalJetTrackPt[9999];
1093  UShort_t FJNTracks=0,EmcalJetNTracks=0;
1094 
1095  if (fTrackCont) for (Int_t i=0; i<fJet->GetNumberOfTracks(); i++) {
1096  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
1097  if (!fTrk) continue;
1098  JetInvMass += fTrk->M();
1099 
1100  fastjet::PseudoJet PseudoTracks(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
1101  TrackMom += TMath::Sqrt(TMath::Power(fTrk->Px(),2)+TMath::Power(fTrk->Py(),2)+TMath::Power(fTrk->Pz(),2));
1102  TrackEnergy += fTrk->E();
1103  PseudoTracks.set_user_index(fJet->TrackAt(i)+100);
1104  PseudJetInvMass += PseudoTracks.m();
1105  fInputVectors.push_back(PseudoTracks);
1106  EmcalJetTrackEta[i]=fTrk->Eta();
1107  EmcalJetTrackPhi[i]=fTrk->Phi();
1108  EmcalJetTrackPt[i]=fTrk->Pt();
1109  EmcalJetNTracks++;
1110  }
1111  fastjet::JetDefinition *fJetDef;
1112  fastjet::ClusterSequence *fClustSeqSA;
1113 
1114 
1115 
1116  fJetDef = new fastjet::JetDefinition(fastjet::antikt_algorithm, fJetRadius*2, static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
1117 
1118  try {
1119  fClustSeqSA = new fastjet::ClusterSequence(fInputVectors, *fJetDef);
1120  } catch (fastjet::Error) {
1121  AliError(" [w] FJ Exception caught.");
1122  //return -1;
1123  }
1124 
1125  std::vector<fastjet::PseudoJet> fOutputJets;
1126  fOutputJets.clear();
1127  fOutputJets=fClustSeqSA->inclusive_jets(0);
1128 
1129 
1130  std::vector<fastjet::PseudoJet> jet_constituents = fOutputJets[0].constituents();
1131  fastjet::contrib::SoftDrop softdrop(beta, zcut);
1132  //fastjet::contrib::SoftDrop softdrop_antikt(beta,zcut);
1133  softdrop.set_verbose_structure(kTRUE);
1134  //fastjet::JetDefinition jet_def_akt(fastjet::antikt_algorithm, 0.4);
1135  // fastjet::contrib::Recluster *antiKT_Recluster(jet_def_akt);
1136  fastjet::contrib::Recluster *recluster;
1137  if(ReclusterAlgo == 1) recluster = new fastjet::contrib::Recluster(fastjet::kt_algorithm,1,true);
1138  if(ReclusterAlgo == 2) recluster = new fastjet::contrib::Recluster(fastjet::antikt_algorithm,1,true);
1139  if(ReclusterAlgo == 0) recluster = new fastjet::contrib::Recluster(fastjet::cambridge_algorithm,1,true);
1140  softdrop.set_reclustering(true,recluster);
1141  fastjet::PseudoJet finaljet = softdrop(fOutputJets[0]);
1142  // fastjet::PseudoJet finaljet_antikt = softdrop_antikt(fOutputJets[0]);
1143  //cout<< finaljet.structure_of<fastjet::contrib::SoftDrop>().symmetry()<<endl;
1144  //cout<< finaljet_antikt.structure_of<fastjet::contrib::SoftDrop>().symmetry()<<endl;
1145 
1146 
1147  //AliEmcalJet* jet = new AliEmcalJet(finaljet.perp(), finaljet.eta(), finaljet.phi(), finaljet.m());
1148  //std::vector<fastjet::PseudoJet> fSDTracks=finaljet.constituents();
1149  //Double_t FastjetTrackDelR,EmcalTrackDelR;
1150  //for(Int_t i=0;i<fJet->GetNumberOfConstituents();i++){
1151  // if(i<=finaljet.constituents().size()){
1152  // FastjetTrackDelR = TMath::Sqrt(TMath::Power(fSDTracks[i].eta()-JetEta,2)+TMath::Power(fSDTracks[i].phi()-JetPhi,2));
1153  // FJTrackEta[i]=fSDTracks[i].eta();
1154  // FJTrackPhi[i]=fSDTracks[i].phi();
1155  // FJTrackPt[i]=fSDTracks[i].perp();
1156  // FJNTracks++;
1157  // }
1158  // AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
1159  // EmcalTrackDelR = TMath::Sqrt(TMath::Power(fTrk->Eta()-JetEta,2)+TMath::Power(fTrk->Phi()-JetPhi,2));
1160  //}
1161  Int_t NDroppedTracks = fJet->GetNumberOfTracks()-finaljet.constituents().size();
1162  //Int_t nConstituents(fClustSeqSA->constituents(finaljet).size());
1163  //jet->SetNumberOfTracks(nConstituents);
1164  Double_t SymParam, Mu, DeltaR, GroomedPt;
1165  Int_t NGroomedBranches;
1166  SymParam=(finaljet.structure_of<fastjet::contrib::SoftDrop>().symmetry());
1167  Mu=(finaljet.structure_of<fastjet::contrib::SoftDrop>().mu());
1168  DeltaR=(finaljet.structure_of<fastjet::contrib::SoftDrop>().delta_R());
1169  NGroomedBranches=finaljet.structure_of<fastjet::contrib::SoftDrop>().dropped_count();
1170  GroomedPt=finaljet.perp();
1171  if(beta==0){
1172  if(ReclusterAlgo==0){
1173  fShapesVar[12]=SymParam;
1174  fShapesVar[13]=DeltaR;
1175  fShapesVar[14]=GroomedPt;
1176  fShapesVar[15]=NGroomedBranches;
1177  fShapesVar[16]=Mu;}
1178  if(ReclusterAlgo==1){
1179  fShapesVar[17]=SymParam;
1180  fShapesVar[18]=DeltaR;
1181  fShapesVar[19]=GroomedPt;
1182  fShapesVar[20]=NGroomedBranches;
1183  fShapesVar[21]=Mu; }
1184 
1185  if(ReclusterAlgo==2){
1186  fShapesVar[22]=SymParam;
1187  fShapesVar[23]=DeltaR;
1188  fShapesVar[24]=GroomedPt;
1189  fShapesVar[25]=NGroomedBranches;
1190  fShapesVar[26]=Mu;
1191  }}
1192  if(beta==1){
1193  fShapesVar[27]=SymParam;
1194  fShapesVar[28]=DeltaR;
1195  fShapesVar[29]=GroomedPt;
1196  fShapesVar[30]=NGroomedBranches;
1197  fShapesVar[31]=Mu;
1198  }
1199 
1200  if(beta==2){
1201  fShapesVar[32]=SymParam;
1202  fShapesVar[33]=DeltaR;
1203  fShapesVar[34]=GroomedPt;
1204  fShapesVar[35]=NGroomedBranches;
1205  fShapesVar[36]=Mu; }
1206 
1207 
1208  return;
1209 
1210 
1211 }
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