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