AliPhysics  ec7afe5 (ec7afe5)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskEmcalJetShapesMC.cxx
Go to the documentation of this file.
1 //
2 // Jet QG tagging analysis task.
3 //
4 // Author: D. Caffarri, L. Cunqueiro
5 
6 #include <TClonesArray.h>
7 #include <TH1F.h>
8 #include <TH2F.h>
9 #include <TH3F.h>
10 #include <THnSparse.h>
11 #include <TTree.h>
12 #include <TList.h>
13 #include <TLorentzVector.h>
14 #include <TProfile.h>
15 #include <TChain.h>
16 #include <TSystem.h>
17 #include <TFile.h>
18 #include <TKey.h>
19 #include <AliAnalysisDataSlot.h>
20 #include <AliAnalysisDataContainer.h>
21 #include "TMatrixD.h"
22 #include "TMatrixDSym.h"
23 #include "TMatrixDSymEigen.h"
24 #include "TVector3.h"
25 #include "TVector2.h"
26 #include "AliVCluster.h"
27 #include "AliVTrack.h"
28 #include "AliEmcalJet.h"
29 #include "AliRhoParameter.h"
30 #include "AliLog.h"
31 #include "AliEmcalParticle.h"
32 #include "AliMCEvent.h"
33 #include "AliGenPythiaEventHeader.h"
34 #include "AliAODMCHeader.h"
35 #include "AliMCEvent.h"
36 #include "AliAnalysisManager.h"
37 #include "AliJetContainer.h"
38 #include "AliParticleContainer.h"
39 #include "AliEmcalPythiaInfo.h"
40 #include "TRandom3.h"
41 #include "AliEmcalJetFinder.h"
42 #include "AliAODEvent.h"
44 
45 using std::cout;
46 using std::endl;
47 
49 
50 //________________________________________________________________________
53  fContainer(0),
54  fMinFractionShared(0),
55  fJetShapeType(kGenShapes),
56  fJetShapeSub(kNoSub),
57  fJetSelection(kInclusive),
58  fPtThreshold(-9999.),
59  fRMatching(0.2),
60  fJetRadius(0.4),
61  fSubjetRadius(0.2),
62  fSelectedShapes(0),
63  fminpTTrig(20.),
64  fmaxpTTrig(50.),
65  fangWindowRecoil(0.6),
66  fSemigoodCorrect(0),
67  fHolePos(0),
68  fHoleWidth(0),
69  fCentSelectOn(kTRUE),
70  fCentMin(0),
71  fCentMax(10),
72  fOneConstSelectOn(kFALSE),
73  fDerivSubtrOrder(0),
74  fPhiJetCorr6(0x0),
75  fPhiJetCorr7(0x0),
76  fEtaJetCorr6(0x0),
77  fEtaJetCorr7(0x0),
78  fPtJetCorr(0x0),
79  fPtJet(0x0),
80  fhpTjetpT(0x0),
81  fhPt(0x0),
82  fhPhi(0x0),
83  fNbOfConstvspT(0x0),
84  fTreeObservableTagging(0x0)
85 
86 {
87  for(Int_t i=0;i<33;i++){
88  fShapesVar[i]=0;}
89  SetMakeGeneralHistograms(kTRUE);
90 }
91 
92 //________________________________________________________________________
94  AliAnalysisTaskEmcalJet(name, kTRUE),
95  fContainer(0),
96  fMinFractionShared(0),
97  fJetShapeType(kGenShapes),
98  fJetShapeSub(kNoSub),
99  fJetSelection(kInclusive),
100  fPtThreshold(-9999.),
101  fRMatching(0.2),
102  fSelectedShapes(0),
103  fminpTTrig(20.),
104  fmaxpTTrig(50.),
105  fangWindowRecoil(0.6),
106  fSemigoodCorrect(0),
107  fHolePos(0),
108  fHoleWidth(0),
109  fCentSelectOn(kTRUE),
110  fCentMin(0),
111  fCentMax(10),
112  fOneConstSelectOn(kFALSE),
113  fDerivSubtrOrder(0),
114  fPhiJetCorr6(0x0),
115  fPhiJetCorr7(0x0),
116  fEtaJetCorr6(0x0),
117  fEtaJetCorr7(0x0),
118  fPtJetCorr(0x0),
119  fPtJet(0x0),
120  fhpTjetpT(0x0),
121  fhPt(0x0),
122  fhPhi(0x0),
123  fNbOfConstvspT(0x0),
124  fTreeObservableTagging(0x0)
125 
126 {
127  // Standard constructor.
128 
129 
130  for(Int_t i=0;i<33;i++){
131  fShapesVar[i]=0;}
132 
134 
135  DefineOutput(1, TList::Class());
136  DefineOutput(2, TTree::Class());
137 
138 
139 }
140 
141 //________________________________________________________________________
143 {
145  delete fTreeObservableTagging;
147  }
148 
149 }
150 
151 //________________________________________________________________________
153 {
154  // Create user output.
156 
157  Bool_t oldStatus = TH1::AddDirectoryStatus();
158  TH1::AddDirectory(oldStatus);
159 
160  //fTreeObservableTagging = new TTree("fTreeJetShape", "fTreeJetShape");
161 
162  //TH1::AddDirectory(oldStatus);
163 
164  const char* nameoutput = GetOutputSlot(2)->GetContainer()->GetName();
165  fTreeObservableTagging = new TTree(nameoutput, nameoutput);
166 
167  const Int_t nVar = 33;
168 
169  TString *fShapesVarNames = new TString [nVar];
170 
171  fShapesVarNames[0] = "partonCode";
172  fShapesVarNames[1] = "ptJet";
173  fShapesVarNames[2] = "ptDJet";
174  fShapesVarNames[3] = "mJet";
175  fShapesVarNames[4] = "nbOfConst";
176  fShapesVarNames[5] = "angularity";
177  fShapesVarNames[6] = "Nsubjet1kt";
178  fShapesVarNames[7] = "Nsubjet2kt";
179  fShapesVarNames[8] = "Nsubjet1Min";
180  fShapesVarNames[9] = "Nsubjet2Min";
181  fShapesVarNames[10] = "DeltaRkt";
182  fShapesVarNames[11] = "DeltaRMin";
183  fShapesVarNames[12] = "SDSymm";
184  fShapesVarNames[13] = "SDDeltaR";
185  fShapesVarNames[14] = "SDGroomedFrac";
186  fShapesVarNames[15] = "SDGroomedN";
187  fShapesVarNames[16] = "SDSymmkt";
188  fShapesVarNames[17] = "SDDeltaRkt";
189  fShapesVarNames[18] = "SDGroomedFrackt";
190  fShapesVarNames[19] = "SDGroomedNkt";
191  fShapesVarNames[20] = "SDSymmAkt";
192  fShapesVarNames[21] = "SDDeltaRAkt";
193  fShapesVarNames[22] = "SDGroomedFracAkt";
194  fShapesVarNames[23] = "SDGroomedNAkt";
195  fShapesVarNames[24] = "SDSymmBeta1";
196  fShapesVarNames[25] = "SDDeltaRBeta1";
197  fShapesVarNames[26] = "SDGroomedFracBeta1";
198  fShapesVarNames[27] = "SDGroomedNBeta1";
199  fShapesVarNames[28] = "SDSymmBeta2";
200  fShapesVarNames[29] = "SDDeltaRBeta2";
201  fShapesVarNames[30] = "SDGroomedFracBeta2";
202  fShapesVarNames[31] = "SDGroomedNBeta2";
203  fShapesVarNames[32] = "weightPythia";
204 
205 
206 
207  //fShapesVarNames[7] = "lesub";
208  //fShapesVarNames[8] = "CoreFraction";
209  //fShapesVarNames[9] = "Nsubjet1";
210  //fShapesVarNames[10] = "Nsubjet2";
211  //fShapesVarNames[11] = "DeltaR";
212  //fShapesVarNames[12] = "OpenAngle";
213  //fShapesVarNames[13] = "weightPythia";
214 
215  //fShapesVarNames[14] = "NT70";
216  //fShapesVarNames[15] = "nConstNT70";
217  //fShapesVarNames[16] = "NT80";
218  //fShapesVarNames[17] = "nConstNT80";
219  //fShapesVarNames[18] = "NT90";
220  //fShapesVarNames[19] = "nConstNT90";
221  //fShapesVarNames[20] = "NT95";
222  //fShapesVarNames[21] = "nConstNT95";
223 
224  //fShapesVarNames[22] = "SubjetFraction";
225 
226 
227  for(Int_t ivar=0; ivar < nVar; ivar++){
228  cout<<"looping over variables"<<endl;
229  fTreeObservableTagging->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/F", fShapesVarNames[ivar].Data()));}
230 
231 
232 
233  fPhiJetCorr6= new TH2F("fPhiJetCorr6", "fPhiJetCorr6", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
234  fOutput->Add(fPhiJetCorr6);
235  fEtaJetCorr6= new TH2F("fEtaJetCorr6", "fEtaJetCorr6", 50, -1.5, 1.5, 50, -1.5, 1.5);
236  fOutput->Add(fEtaJetCorr6);
237 
238  fPhiJetCorr7= new TH2F("fPhiJetCorr7", "fPhiJetCorr7", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
239  fOutput->Add(fPhiJetCorr7);
240  fEtaJetCorr7= new TH2F("fEtaJetCorr7", "fEtaJetCorr7", 50, -1.5, 1.5, 50, -1.5, 1.5);
241  fOutput->Add(fEtaJetCorr7);
242 
243  fPtJetCorr= new TH2F("fPtJetCorr", "fPtJetCorr", 100, 0, 200, 100, 0, 200);
244  fOutput->Add(fPtJetCorr);
245  fPtJet= new TH1F("fPtJet", "fPtJet", 100, 0, 200);
246  fOutput->Add(fPtJet);
247 
248  fhpTjetpT= new TH2F("fhpTjetpT", "fhpTjetpT", 200, 0, 200, 200, 0, 200);
249  fOutput->Add(fhpTjetpT);
250  fhPt= new TH1F("fhPt", "fhPt", 200, 0, 200);
251  fOutput->Add(fhPt);
252  fhPhi= new TH1F("fhPhi", "fhPhi", 100, -TMath::Pi(), TMath::Pi());
253  fOutput->Add(fhPhi);
254 
255  fNbOfConstvspT=new TH2F("fNbOfConstvspT", "fNbOfConstvspT", 100, 0, 100, 200, 0, 200);
256  fOutput->Add(fNbOfConstvspT);
257 
258  //fOutput->Add(fTreeObservableTagging);
259 
260 
261  PostData(1, fOutput); // Post data for ALL output slots > 0 here
262  PostData(2, fTreeObservableTagging);
263 
264  delete [] fShapesVarNames;
265 
266 }
267 
268 //________________________________________________________________________
270 {
271  // Run analysis code here, if needed. It will be executed before FillHistograms().
272 
273  return kTRUE;
274 }
275 
276 //________________________________________________________________________
278 {
279  // Fill histograms.
280  //cout<<"IntoFillHistograms"<<endl;
281  AliEmcalJet* jet1 = NULL;
282  AliJetContainer *jetCont = GetJetContainer(0);
283 
284  Float_t kWeight=1;
285  if (fCentSelectOn)
286  if ((fCent>fCentMax) || (fCent<fCentMin)) return 0;
287 
288  AliAODTrack *triggerHadron = 0x0;
289 
290  if (fJetSelection == kRecoil) {
291  //Printf("Recoil jets!!!, fminpTTrig = %f, fmaxpTTrig = %f", fminpTTrig, fmaxpTTrig);
292  Int_t triggerHadronLabel = SelectTrigger(fminpTTrig, fmaxpTTrig);
293 
294 
295  if (triggerHadronLabel==-99999) {
296  //Printf ("Trigger Hadron not found, return");
297  return 0;}
298 
299 
301  TClonesArray *trackArrayAn = partContAn->GetArray();
302  triggerHadron = static_cast<AliAODTrack*>(trackArrayAn->At(triggerHadronLabel));
303 
304  if (!triggerHadron) {
305  //Printf("No Trigger hadron with the found label!!");
306  return 0;
307  }
308 
309  if(fSemigoodCorrect){
310  Double_t disthole=RelativePhi(triggerHadron->Phi(),fHolePos);
311  if(TMath::Abs(disthole)+fHoleWidth>TMath::Pi()-fangWindowRecoil){
312  return 0;}
313  }
314 
315  fhPt->Fill(triggerHadron->Pt());
316 
317  }
318 
319  if(jetCont) {
320  jetCont->ResetCurrentID();
321  while((jet1 = jetCont->GetNextAcceptJet())) {
322  //Printf("jet1=%p", jet1);
323  if (!jet1) continue;
324  AliEmcalJet* jet2 = 0x0;
325  AliEmcalJet* jet3 = 0x0;
326  fPtJet->Fill(jet1->Pt());
327  AliEmcalJet *jetUS = NULL;
328  Int_t ifound=0, jfound=0;
329  Int_t ilab=-1, jlab=-1;
330 
332  Double_t disthole=RelativePhi(jet1->Phi(),fHolePos);
333  if(TMath::Abs(disthole)<fHoleWidth){
334  continue;
335  }
336  }
337 
338  Float_t dphiRecoil = 0.;
339  if (fJetSelection == kRecoil){
340  dphiRecoil = RelativePhi(triggerHadron->Phi(), jet1->Phi());
341  if (TMath::Abs(dphiRecoil) < (TMath::Pi() - fangWindowRecoil)) {
342  // Printf("Recoil jets back to back not found! continuing");
343  continue;
344  }
345 
346  fhpTjetpT->Fill(triggerHadron->Pt(), jet1->Pt());
347  //Printf(" ************ FILLING HISTOS****** shapeSub = %d, triggerHadron = %f, jet1 = %f", fJetShapeSub, triggerHadron->Pt(), jet1->Pt());
348  fhPhi->Fill(RelativePhi(triggerHadron->Phi(), jet1->Phi()));
349 
350  }
351 
352 
353  fShapesVar[0] = 0.;
354 
355  if (fJetShapeType == kGenShapes){
356  const AliEmcalPythiaInfo *partonsInfo = 0x0;
357  partonsInfo = GetPythiaInfo();
358  //Printf("partonsInfo=%p", partonsInfo);
359  Double_t jp1=RelativePhi(jet1->Phi(),partonsInfo->GetPartonPhi6());
360  Double_t detap1=(jet1->Eta())-(partonsInfo->GetPartonEta6());
361  kWeight=partonsInfo->GetPythiaEventWeight();
362  //Printf("kWeight=%f", kWeight);
363  fShapesVar[32] = kWeight;
364 
365  Float_t dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
366  fEtaJetCorr6->Fill(jet1->Eta(), partonsInfo->GetPartonEta6());
367  fPhiJetCorr6->Fill(jet1->Phi(), partonsInfo->GetPartonPhi6());
368  if(dRp1 < fRMatching) {
369  fShapesVar[0] = partonsInfo->GetPartonFlag6();
370  fPtJetCorr ->Fill(partonsInfo->GetPartonPt6(), jet1->Pt());
371  }
372  else {
373  jp1=RelativePhi(jet1->Phi(),partonsInfo->GetPartonPhi7());
374  detap1=(jet1->Eta())-(partonsInfo->GetPartonEta7());
375  dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
376  fEtaJetCorr7->Fill(jet1->Eta(), partonsInfo->GetPartonEta7());
377  fPhiJetCorr7->Fill(jet1->Phi(), partonsInfo->GetPartonPhi7());
378  if(dRp1 < fRMatching) {
379  fShapesVar[0] = partonsInfo->GetPartonFlag7();
380  fPtJetCorr->Fill(partonsInfo->GetPartonPt7(), jet1->Pt());
381  }
382  else fShapesVar[0]=0;
383  }
384  }
385 
386  Double_t ptSubtracted = 0;
387  ptSubtracted= jet1->Pt();
388  //Printf("ptSubtracted=%f", ptSubtracted);
389 
390 
391  if (ptSubtracted < fPtThreshold) continue;
392 
393  if ((fCentSelectOn == kFALSE) && (jet1->GetNumberOfTracks() < 2)) continue;
394 
395  //AliEmcalJetFinder *Reclusterer1; //Object containg Subjets from Subtracted Hybrid Jets
396  //Reclusterer1 = Recluster(jet1, 0, fJetRadius, fSubjetRadius, 1, 0, "SubJetFinder_1");
397 
398 
399 
400 
401 
402 
403 
404  fShapesVar[1] = ptSubtracted;
405  fShapesVar[2] = GetJetpTD(jet1,0);
406  fShapesVar[3] = GetJetMass(jet1,0);
407  fShapesVar[4] = 1.*GetJetNumberOfConstituents(jet1,0);
408  fShapesVar[5] = GetJetAngularity(jet1,0);
409  //nsub1 and nsub2 for kT
410  fShapesVar[6] =0;
411  //fjNSubJettiness(jet1,0,1,0,1,0);
412  fShapesVar[7]=0;
413  //fjNSubJettiness(jet1,0,2,0,1,0);
414  //nsub1 and nsub2 for min_axis
415  fShapesVar[8]=0;
416  //fjNSubJettiness(jet1,0,1,10,1,0);
417  fShapesVar[9]=0;
418  //fjNSubJettiness(jet1,0,2,10,1,0);
419  //deltaRkt
420  fShapesVar[10] =0;
421  //fjNSubJettiness(jet1,0,1,0,1,2);
422  //deltaRmin
423  fShapesVar[11]=0;
424  //fjNSubJettiness(jet1,0,2,10,1,2);
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  //Printf("jetfinder =%p, JetCont=%p", JetFinder, JetCont);
963  JetFinder->SetJetMaxEta(0.9-fJetRadius);
964  JetFinder->SetRadius(fJetRadius);
965  JetFinder->SetJetAlgorithm(Algorithm); //0 for anti-kt 1 for kt //this is for the JET!!!!!!!!!! Not the SubJets
966  JetFinder->SetRecombSheme(0);
967  JetFinder->SetJetMinPt(Jet->Pt());
968  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
969  //Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
970  Double_t dVtx[3]={0,0,0};
971  //Printf("JetFinder->Nsubjettiness =%f", JetFinder->Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,fSubjetRadius,Beta,Option));
972  return JetFinder->Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,fSubjetRadius,Beta,Option);
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  cout<<"HELLO"<<endl;
1088  std::vector<fastjet::PseudoJet> fInputVectors;
1089  Double_t JetInvMass=0, PseudJetInvMass=0, TrackMom = 0, TrackEnergy = 0;
1090 
1091  AliParticleContainer *fTrackCont = fJetCont->GetParticleContainer();
1092  cout<<"it is here"<<fTrackCont<<endl;
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  cout<<fJet->GetNumberOfTracks()<<" coco"<<endl;
1097  if (fTrackCont) for (Int_t i=0; i<fJet->GetNumberOfTracks(); i++) {
1098  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
1099  JetInvMass += fTrk->M();
1100  if (!fTrk) continue;
1101  fastjet::PseudoJet PseudoTracks(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
1102  TrackMom += TMath::Sqrt(TMath::Power(fTrk->Px(),2)+TMath::Power(fTrk->Py(),2)+TMath::Power(fTrk->Pz(),2));
1103  TrackEnergy += fTrk->E();
1104  PseudoTracks.set_user_index(fJet->TrackAt(i)+100);
1105  PseudJetInvMass += PseudoTracks.m();
1106  fInputVectors.push_back(PseudoTracks);
1107  EmcalJetTrackEta[i]=fTrk->Eta();
1108  EmcalJetTrackPhi[i]=fTrk->Phi();
1109  EmcalJetTrackPt[i]=fTrk->Pt();
1110  EmcalJetNTracks++;
1111  }
1112  fastjet::JetDefinition *fJetDef;
1113  fastjet::ClusterSequence *fClustSeqSA;
1114 
1115 
1116 
1117  fJetDef = new fastjet::JetDefinition(fastjet::antikt_algorithm, fJetRadius*2, static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
1118 
1119  try {
1120  fClustSeqSA = new fastjet::ClusterSequence(fInputVectors, *fJetDef);
1121  } catch (fastjet::Error) {
1122  AliError(" [w] FJ Exception caught.");
1123  //return -1;
1124  }
1125 
1126  std::vector<fastjet::PseudoJet> fOutputJets;
1127  fOutputJets.clear();
1128  fOutputJets=fClustSeqSA->inclusive_jets(0);
1129 
1130 
1131  std::vector<fastjet::PseudoJet> jet_constituents = fOutputJets[0].constituents();
1132  fastjet::contrib::SoftDrop softdrop(beta, zcut);
1133  //fastjet::contrib::SoftDrop softdrop_antikt(beta,zcut);
1134  softdrop.set_verbose_structure(kTRUE);
1135  //fastjet::JetDefinition jet_def_akt(fastjet::antikt_algorithm, 0.4);
1136  // fastjet::contrib::Recluster *antiKT_Recluster(jet_def_akt);
1137  fastjet::contrib::Recluster *recluster;
1138  if(ReclusterAlgo == 1) recluster = new fastjet::contrib::Recluster(fastjet::kt_algorithm,1,true);
1139  if(ReclusterAlgo == 2) recluster = new fastjet::contrib::Recluster(fastjet::antikt_algorithm,1,true);
1140  if(ReclusterAlgo == 0) recluster = new fastjet::contrib::Recluster(fastjet::cambridge_algorithm,1,true);
1141  softdrop.set_reclustering(true,recluster);
1142  fastjet::PseudoJet finaljet = softdrop(fOutputJets[0]);
1143  // fastjet::PseudoJet finaljet_antikt = softdrop_antikt(fOutputJets[0]);
1144  //cout<< finaljet.structure_of<fastjet::contrib::SoftDrop>().symmetry()<<endl;
1145  //cout<< finaljet_antikt.structure_of<fastjet::contrib::SoftDrop>().symmetry()<<endl;
1146 
1147 
1148  AliEmcalJet* jet = new AliEmcalJet(finaljet.perp(), finaljet.eta(), finaljet.phi(), finaljet.m());
1149  std::vector<fastjet::PseudoJet> fSDTracks=finaljet.constituents();
1150  Double_t FastjetTrackDelR,EmcalTrackDelR;
1151  for(Int_t i=0;i<fJet->GetNumberOfConstituents();i++){
1152  if(i<=finaljet.constituents().size()){
1153  FastjetTrackDelR = TMath::Sqrt(TMath::Power(fSDTracks[i].eta()-JetEta,2)+TMath::Power(fSDTracks[i].phi()-JetPhi,2));
1154  FJTrackEta[i]=fSDTracks[i].eta();
1155  FJTrackPhi[i]=fSDTracks[i].phi();
1156  FJTrackPt[i]=fSDTracks[i].perp();
1157  FJNTracks++;
1158  }
1159  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
1160  EmcalTrackDelR = TMath::Sqrt(TMath::Power(fTrk->Eta()-JetEta,2)+TMath::Power(fTrk->Phi()-JetPhi,2));
1161  }
1162  Int_t NDroppedTracks = fJet->GetNumberOfTracks()-finaljet.constituents().size();
1163  Int_t nConstituents(fClustSeqSA->constituents(finaljet).size());
1164  jet->SetNumberOfTracks(nConstituents);
1165  Double_t SymParam, Mu, DeltaR, GroomedPt;
1166  Int_t NGroomedBranches;
1167  SymParam=(finaljet.structure_of<fastjet::contrib::SoftDrop>().symmetry());
1168  Mu=(finaljet.structure_of<fastjet::contrib::SoftDrop>().mu());
1169  DeltaR=(finaljet.structure_of<fastjet::contrib::SoftDrop>().delta_R());
1170  NGroomedBranches=finaljet.structure_of<fastjet::contrib::SoftDrop>().dropped_count();
1171  GroomedPt=finaljet.perp();
1172  if(beta==0){
1173  if(ReclusterAlgo==0){
1174  fShapesVar[12]=SymParam;
1175  fShapesVar[13]=DeltaR;
1176  fShapesVar[14]=GroomedPt;
1177  fShapesVar[15]=NGroomedBranches;}
1178  if(ReclusterAlgo==1){
1179  fShapesVar[16]=SymParam;
1180  fShapesVar[17]=DeltaR;
1181  fShapesVar[18]=GroomedPt;
1182  fShapesVar[19]=NGroomedBranches;}
1183 
1184  if(ReclusterAlgo==2){
1185  fShapesVar[20]=SymParam;
1186  fShapesVar[21]=DeltaR;
1187  fShapesVar[22]=GroomedPt;
1188  fShapesVar[23]=NGroomedBranches;}}
1189  if(beta==1){
1190  fShapesVar[24]=SymParam;
1191  fShapesVar[25]=DeltaR;
1192  fShapesVar[26]=GroomedPt;
1193  fShapesVar[27]=NGroomedBranches;}
1194 
1195  if(beta==2){
1196  fShapesVar[28]=SymParam;
1197  fShapesVar[29]=DeltaR;
1198  fShapesVar[30]=GroomedPt;
1199  fShapesVar[31]=NGroomedBranches;}
1200 
1201 
1202  return;
1203 
1204 
1205 }
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
UShort_t GetNumberOfConstituents() const
Definition: AliEmcalJet.h:133
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
void SetNumberOfTracks(Int_t n)
Definition: AliEmcalJet.h:194
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:255
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