AliPhysics  608b256 (608b256)
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 "TF1.h"
42 #include "AliEmcalJetFinder.h"
43 #include "AliAODEvent.h"
45 
46 #include <fastjet/config.h>
47 #if FASJET_VERSION_NUMBER >= 30302
48 #include <fastjet/tools/Recluster.hh>
49 #else
50 #include <fastjet/contrib/Recluster.hh>
51 #endif
52 
53 using std::cout;
54 using std::endl;
55 
57 
58 //________________________________________________________________________
61  fContainer(0),
62  fMinFractionShared(0),
63  fJetShapeType(kGenShapes),
64  fJetShapeSub(kNoSub),
65  fJetSelection(kInclusive),
66  fPtThreshold(-9999.),
67  fRMatching(0.2),
68  fJetRadius(0.4),
69  fSubjetRadius(0.2),
70  fSelectedShapes(0),
71  fSwitchKtNSub(0),
72  fSwitchMinNSub(0),
73  fSwitchAktNSub(0),
74  fSwitchSDKtNSub(0),
75  fSwitchSDMinNSub(0),
76  fAdditionalTracks(0),
77  fHardCutoff(0),
78  fOptionalPartonInfo(0),
79  fminpTTrig(20.),
80  fmaxpTTrig(50.),
81  fangWindowRecoil(0.6),
82  fSemigoodCorrect(0),
83  fHolePos(0),
84  fHoleWidth(0),
85  fRandom(0),
86  fqhat(1),
87  fxlength(2),
88  fCentSelectOn(kTRUE),
89  fCentMin(0),
90  fCentMax(10),
91  fOneConstSelectOn(kFALSE),
92  fDerivSubtrOrder(0),
93  fPhiJetCorr6(0x0),
94  fPhiJetCorr7(0x0),
95  fEtaJetCorr6(0x0),
96  fEtaJetCorr7(0x0),
97  fPtJetCorr(0x0),
98  fPtJet(0x0),
99  fhpTjetpT(0x0),
100  fhPt(0x0),
101  fhPhi(0x0),
102  fHLundIterative(0x0),
103  fHLundIterative_ktaxis(0x0),
104  fHLundIterativeInject(0x0),
105  fNbOfConstvspT(0x0),
106  fTreeObservableTagging(0x0),
107  fTf1SoftOmega(0x0),
108  fTf1SoftKt(0x0),
109  fTf1Omega(0x0),
110  fTf1Kt(0x0),
111  fScaleELoss(kFALSE),
112  xfraction(1),
113  fAddMedScat(kFALSE),
114  fAddMedScatPtFrac(1),
115  fAddMedScatN(100)
116 
117 
118 {
119  for(Int_t i=0;i<13;i++){
120  fShapesVar[i]=0;}
121  SetMakeGeneralHistograms(kTRUE);
122  DefineOutput(1, TList::Class());
123  DefineOutput(2, TTree::Class());
124 }
125 
126 //________________________________________________________________________
128  AliAnalysisTaskEmcalJet(name, kTRUE),
129  fContainer(0),
130  fMinFractionShared(0),
131  fJetShapeType(kGenShapes),
132  fJetShapeSub(kNoSub),
133  fJetSelection(kInclusive),
134  fPtThreshold(-9999.),
135  fRMatching(0.2),
136  fSelectedShapes(0),
137  fSwitchKtNSub(0),
138  fSwitchMinNSub(0),
139  fSwitchAktNSub(0),
140  fSwitchSDKtNSub(0),
141  fSwitchSDMinNSub(0),
142  fAdditionalTracks(0),
143  fHardCutoff(0),
144  fOptionalPartonInfo(0),
145  fminpTTrig(20.),
146  fmaxpTTrig(50.),
147  fangWindowRecoil(0.6),
148  fSemigoodCorrect(0),
149  fHolePos(0),
150  fHoleWidth(0),
151  fRandom(0),
152  fqhat(1),
153  fxlength(2),
154  fCentSelectOn(kTRUE),
155  fCentMin(0),
156  fCentMax(10),
157  fOneConstSelectOn(kFALSE),
158  fDerivSubtrOrder(0),
159  fPhiJetCorr6(0x0),
160  fPhiJetCorr7(0x0),
161  fEtaJetCorr6(0x0),
162  fEtaJetCorr7(0x0),
163  fPtJetCorr(0x0),
164  fPtJet(0x0),
165  fhpTjetpT(0x0),
166  fhPt(0x0),
167  fhPhi(0x0),
168  fHLundIterative(0x0),
169  fHLundIterative_ktaxis(0x0),
170  fHLundIterativeInject(0x0),
171  fNbOfConstvspT(0x0),
172  fTreeObservableTagging(0x0),
173  fTf1SoftOmega(0x0),
174  fTf1SoftKt(0x0),
175  fTf1Omega(0x0),
176  fTf1Kt(0x0),
177  fScaleELoss(kFALSE),
178  xfraction(1),
179  fAddMedScat(kFALSE),
180  fAddMedScatPtFrac(1),
181  fAddMedScatN(100)
182 {
183  // Standard constructor.
184 
185 
186  for(Int_t i=0;i<13;i++){
187  fShapesVar[i]=0;}
188 
190 
191  DefineOutput(1, TList::Class());
192  DefineOutput(2, TTree::Class());
193 
194 
195 }
196 
197 //________________________________________________________________________
199 {
201  delete fTreeObservableTagging;
203  }
204 
205  if(fRandom) delete fRandom;
206  if(fTf1SoftOmega) delete fTf1SoftOmega;
207  if(fTf1SoftKt) delete fTf1SoftKt;
208  if(fTf1Omega) delete fTf1Omega;
209  if(fTf1Kt) delete fTf1Kt;
210 
211 }
212 
213 //________________________________________________________________________
215 {
216  // Create user output.
218 
219  Bool_t oldStatus = TH1::AddDirectoryStatus();
220  TH1::AddDirectory(oldStatus);
221 
222  //fTreeObservableTagging = new TTree("fTreeJetShape", "fTreeJetShape");
223 
224  //TH1::AddDirectory(oldStatus);
225 
226  const char* nameoutput = GetOutputSlot(2)->GetContainer()->GetName();
227  fTreeObservableTagging = new TTree(nameoutput, nameoutput);
228 
229  const Int_t nVar = 13;
230 
231  TString *fShapesVarNames = new TString [nVar];
232 
233  fShapesVarNames[0] = "partonCode";
234  fShapesVarNames[1] = "ptJet";
235  fShapesVarNames[2] = "ptDJet";
236  fShapesVarNames[3] = "mJet";
237  fShapesVarNames[4] = "nbOfConst";
238  fShapesVarNames[5] = "angularity";
239  fShapesVarNames[6] = "nitersd";
240  fShapesVarNames[7] = "niterall";
241  fShapesVarNames[8] = "weightPythia";
242  //fShapesVarNames[6] = "Nsubjet1kt";
243  // fShapesVarNames[7] = "Nsubjet2kt";
244  // fShapesVarNames[8] = "Nsubjet1Min";
245  // fShapesVarNames[9] = "Nsubjet2Min";
246  // fShapesVarNames[10] = "DeltaRkt";
247  // fShapesVarNames[11] = "DeltaRMin";
248  fShapesVarNames[9] = "SDSymm";
249  fShapesVarNames[10] = "scaledptJet";
250  fShapesVarNames[11] = "zg";
251  fShapesVarNames[12] = "rg";
252  // fShapesVarNames[15] = "SDGroomedN";
253  // fShapesVarNames[16] = "SDMass";
254  // fShapesVarNames[17] = "SDSymmkt";
255  // fShapesVarNames[18] = "SDDeltaRkt";
256  // fShapesVarNames[19] = "SDGroomedFrackt";
257  // fShapesVarNames[20] = "SDGroomedNkt";
258  // fShapesVarNames[21] = "SDMasskt";
259  // fShapesVarNames[22] = "SDSymmAkt";
260  // fShapesVarNames[23] = "SDDeltaRAkt";
261  // fShapesVarNames[24] = "SDGroomedFracAkt";
262  // fShapesVarNames[25] = "SDGroomedNAkt";
263  // fShapesVarNames[26] = "SDMassAkt";
264  // fShapesVarNames[27] = "SDSymmktForm";
265  // fShapesVarNames[28] = "SDDeltaRktForm";
266  // fShapesVarNames[29] = "SDGroomedFracktForm";
267  // fShapesVarNames[30] = "SDGroomedNktForm";
268  // fShapesVarNames[31] = "SDMassktForm";
269  // fShapesVarNames[32] = "SDSymmDemo";
270  // fShapesVarNames[33] = "SDDeltaRDemo";
271  // fShapesVarNames[34] = "SDGroomedFracDemo";
272  // fShapesVarNames[35] = "SDGroomedNDemo";
273  // fShapesVarNames[36] = "SDMassDemo";
274  // fShapesVarNames[42] = "SDSymmForm";
275  // fShapesVarNames[43] = "SDDeltaRForm";
276  // fShapesVarNames[44] = "SDGroomedFracForm";
277  // fShapesVarNames[45] = "SDGroomedNForm";
278  // fShapesVarNames[46] = "SDMassForm";
279  // fShapesVarNames[47] = "weightPythia";
280  // fShapesVarNames[42] = "SDSymmNoCut";
281  // fShapesVarNames[43] = "SDDeltaRNoCut";
282  // fShapesVarNames[44] = "SDGroomedFracNoCut";
283  // fShapesVarNames[45] = "SDGroomedNNoCut";
284  // fShapesVarNames[46] = "SDMassNoCut";
285 
286 
287  //fShapesVarNames[7] = "lesub";
288  //fShapesVarNames[8] = "CoreFraction";
289  //fShapesVarNames[9] = "Nsubjet1";
290  //fShapesVarNames[10] = "Nsubjet2";
291  //fShapesVarNames[11] = "DeltaR";
292  //fShapesVarNames[12] = "OpenAngle";
293  //fShapesVarNames[13] = "weightPythia";
294 
295  //fShapesVarNames[14] = "NT70";
296  //fShapesVarNames[15] = "nConstNT70";
297  //fShapesVarNames[16] = "NT80";
298  //fShapesVarNames[17] = "nConstNT80";
299  //fShapesVarNames[18] = "NT90";
300  //fShapesVarNames[19] = "nConstNT90";
301  //fShapesVarNames[20] = "NT95";
302  //fShapesVarNames[21] = "nConstNT95";
303 
304  //fShapesVarNames[22] = "SubjetFraction";
305 
306 
307  for(Int_t ivar=0; ivar < nVar; ivar++){
308  cout<<"looping over variables"<<endl;
309  fTreeObservableTagging->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/F", fShapesVarNames[ivar].Data()));}
310 
311 
312 
313  fPhiJetCorr6= new TH2F("fPhiJetCorr6", "fPhiJetCorr6", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
314  fOutput->Add(fPhiJetCorr6);
315  fEtaJetCorr6= new TH2F("fEtaJetCorr6", "fEtaJetCorr6", 50, -1.5, 1.5, 50, -1.5, 1.5);
316  fOutput->Add(fEtaJetCorr6);
317 
318  fPhiJetCorr7= new TH2F("fPhiJetCorr7", "fPhiJetCorr7", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
319  fOutput->Add(fPhiJetCorr7);
320  fEtaJetCorr7= new TH2F("fEtaJetCorr7", "fEtaJetCorr7", 50, -1.5, 1.5, 50, -1.5, 1.5);
321  fOutput->Add(fEtaJetCorr7);
322 
323  fPtJetCorr= new TH2F("fPtJetCorr", "fPtJetCorr", 100, 0, 200, 100, 0, 200);
324  fOutput->Add(fPtJetCorr);
325  fPtJet= new TH1F("fPtJet", "fPtJet", 100, 0, 200);
326  fOutput->Add(fPtJet);
327 
328  fhpTjetpT= new TH2F("fhpTjetpT", "fhpTjetpT", 200, 0, 200, 200, 0, 200);
329  fOutput->Add(fhpTjetpT);
330  fhPt= new TH1F("fhPt", "fhPt", 200, 0, 200);
331  fOutput->Add(fhPt);
332  fhPhi= new TH1F("fhPhi", "fhPhi", 100, -TMath::Pi(), TMath::Pi());
333  fOutput->Add(fhPhi);
334 
335 
336 
337  //log(1/theta),log(z*theta),jetpT,algo//
338  const Int_t dimSpec = 6;
339  const Int_t nBinsSpec[6] = {50,100,60,3,22,10};
340  const Double_t lowBinSpec[6] = {0.9,-10,0.0,0.0,0.0,0.0};
341  const Double_t hiBinSpec[6] = {5.0,2.0,600.0,3.0,22.0,10.0};
342  fHLundIterative = new THnSparseF("fHLundIterative",
343  "LundIterativePlot [log(1/theta),log(z*theta),pTjet,algo,partonFlavor,depth]",
344  dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
345  fOutput->Add(fHLundIterative);
346  const Int_t nBinsSpec_ktaxis[6] = {50,100,60,3,22,10};
347  const Double_t lowBinSpec_ktaxis[6] = {0.9,-4.0,0.0,0.0,0.0,0.0};
348  const Double_t hiBinSpec_ktaxis[6] = {5.0,8.0,600.0,3.0,22.0,10.0};
349  fHLundIterative_ktaxis = new THnSparseF("fHLundIterative_ktaxis",
350  "LundIterativePlot [log(1/theta),log(z*theta),pTjet,algo,partonFlavor,depth]",
351  dimSpec,nBinsSpec_ktaxis,lowBinSpec_ktaxis,hiBinSpec_ktaxis);
353 
354  if(fAdditionalTracks>0){
355  //log(1/theta),log(z*theta),jetpT of added tracks
356  const Int_t dimSpecb = 4;
357  const Int_t nBinsSpecb[4] = {50,50,20,20};
358  const Double_t lowBinSpecb[4] = {0.0,-10, 0,0};
359  const Double_t hiBinSpecb[4] = {5.0, 0,200,200};
360  fHLundIterativeInject = new THnSparseF("fHLundIterativeInject",
361  "LundIterativePlotInject [log(1/theta),log(z*theta),pTjet,algo]",
362  dimSpecb,nBinsSpecb,lowBinSpecb,hiBinSpecb);
364 
365 
366  fNbOfConstvspT=new TH2F("fNbOfConstvspT", "fNbOfConstvspT", 100, 0, 100, 200, 0, 200);
367  fOutput->Add(fNbOfConstvspT);
368 
369  //fOutput->Add(fTreeObservableTagging);
370 
371  fRandom = new TRandom3(0);
372  PostData(1, fOutput); // Post data for ALL output slots > 0 here
373  PostData(2, fTreeObservableTagging);
374 
375  delete [] fShapesVarNames;
376 
377  }
378 
379 //________________________________________________________________________
381 {
382  // Run analysis code here, if needed. It will be executed before FillHistograms().
383  // if (gRandom) delete gRandom;
384  // gRandom = new TRandom3(0);
385  return kTRUE;
386 }
387 
388 //________________________________________________________________________
390 {
391  // Fill histograms.
392  //cout<<"IntoFillHistograms"<<endl;
393  AliEmcalJet* jet1 = NULL;
394  AliJetContainer *jetCont = GetJetContainer(0);
395 
396  Float_t kWeight=1;
397  if (fCentSelectOn)
398  if ((fCent>fCentMax) || (fCent<fCentMin)) return 0;
399 
400  AliAODTrack *triggerHadron = 0x0;
401 
402  if (fJetSelection == kRecoil) {
403  //Printf("Recoil jets!!!, fminpTTrig = %f, fmaxpTTrig = %f", fminpTTrig, fmaxpTTrig);
404  Int_t triggerHadronLabel = SelectTrigger(fminpTTrig, fmaxpTTrig);
405 
406 
407  if (triggerHadronLabel==-99999) {
408  //Printf ("Trigger Hadron not found, return");
409  return 0;}
410 
411 
413  TClonesArray *trackArrayAn = partContAn->GetArray();
414  triggerHadron = static_cast<AliAODTrack*>(trackArrayAn->At(triggerHadronLabel));
415 
416  if (!triggerHadron) {
417  //Printf("No Trigger hadron with the found label!!");
418  return 0;
419  }
420 
421  if(fSemigoodCorrect){
422  Double_t disthole=RelativePhi(triggerHadron->Phi(),fHolePos);
423  if(TMath::Abs(disthole)+fHoleWidth>TMath::Pi()-fangWindowRecoil){
424  return 0;}
425  }
426 
427  fhPt->Fill(triggerHadron->Pt());
428 
429  }
430 
431  if(jetCont) {
432  jetCont->ResetCurrentID();
433  while((jet1 = jetCont->GetNextAcceptJet())) {
434  //Printf("jet1=%p", jet1);
435  if (!jet1) continue;
436 
437  fPtJet->Fill(jet1->Pt());
438 
439 
440 
441 
443  Double_t disthole=RelativePhi(jet1->Phi(),fHolePos);
444  if(TMath::Abs(disthole)<fHoleWidth){
445  continue;
446  }
447  }
448 
449  Float_t dphiRecoil = 0.;
450  if (fJetSelection == kRecoil){
451  dphiRecoil = RelativePhi(triggerHadron->Phi(), jet1->Phi());
452  if (TMath::Abs(dphiRecoil) < (TMath::Pi() - fangWindowRecoil)) {
453  // Printf("Recoil jets back to back not found! continuing");
454  continue;
455  }
456 
457  fhpTjetpT->Fill(triggerHadron->Pt(), jet1->Pt());
458  //Printf(" ************ FILLING HISTOS****** shapeSub = %d, triggerHadron = %f, jet1 = %f", fJetShapeSub, triggerHadron->Pt(), jet1->Pt());
459  fhPhi->Fill(RelativePhi(triggerHadron->Phi(), jet1->Phi()));
460 
461  }
462 
463 
464  fShapesVar[0] = 0.;
465 
466  if (fOptionalPartonInfo==1){
467  const AliEmcalPythiaInfo *partonsInfo = 0x0;
468  partonsInfo = GetPythiaInfo();
469  //Printf("partonsInfo=%p", partonsInfo);
470  Double_t jp1=RelativePhi(jet1->Phi(),partonsInfo->GetPartonPhi6());
471  Double_t detap1=(jet1->Eta())-(partonsInfo->GetPartonEta6());
472  kWeight=partonsInfo->GetPythiaEventWeight();
473  //Printf("kWeight=%f", kWeight);
474  fShapesVar[8] = kWeight;
475 
476  Float_t dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
477  fEtaJetCorr6->Fill(jet1->Eta(), partonsInfo->GetPartonEta6());
478  fPhiJetCorr6->Fill(jet1->Phi(), partonsInfo->GetPartonPhi6());
479  if(dRp1 < fRMatching) {
480  fShapesVar[0] = partonsInfo->GetPartonFlag6();
481  fPtJetCorr ->Fill(partonsInfo->GetPartonPt6(), jet1->Pt());
482  }
483  else {
484  jp1=RelativePhi(jet1->Phi(),partonsInfo->GetPartonPhi7());
485  detap1=(jet1->Eta())-(partonsInfo->GetPartonEta7());
486  dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
487  fEtaJetCorr7->Fill(jet1->Eta(), partonsInfo->GetPartonEta7());
488  fPhiJetCorr7->Fill(jet1->Phi(), partonsInfo->GetPartonPhi7());
489  if(dRp1 < fRMatching) {
490  fShapesVar[0] = partonsInfo->GetPartonFlag7();
491  fPtJetCorr->Fill(partonsInfo->GetPartonPt7(), jet1->Pt());
492  }
493  else fShapesVar[0]=0;
494  }
495  }
496 
497  Double_t ptSubtracted = 0;
498  ptSubtracted= jet1->Pt();
499  //Printf("ptSubtracted=%f", ptSubtracted);
500 
501 
502  if (ptSubtracted < fPtThreshold) continue;
503 
504  if ((fCentSelectOn == kFALSE) && (jet1->GetNumberOfTracks() < 2)) continue;
505 
506  //AliEmcalJetFinder *Reclusterer1; //Object containg Subjets from Subtracted Hybrid Jets
507  //Reclusterer1 = Recluster(jet1, 0, fJetRadius, fSubjetRadius, 1, 0, "SubJetFinder_1");
508 
509 
510 
511 
512 
513 
514 
515  fShapesVar[1] = ptSubtracted;
516  fShapesVar[2] = GetJetpTD(jet1,0);
517  fShapesVar[3] = GetJetMass(jet1,0);
518  fShapesVar[4] = 1.*GetJetNumberOfConstituents(jet1,0);
519  fShapesVar[5] = GetJetAngularity(jet1,0);
520  SoftDrop(jet1,jetCont,0.1,0,0);
521 
522  RecursiveParents(jet1,jetCont,0,fShapesVar[0]);
523  RecursiveParents(jet1,jetCont,1,fShapesVar[0]);
524  RecursiveParents(jet1,jetCont,2,fShapesVar[0]);
525 
526 
527  fTreeObservableTagging->Fill();
528 
529  }
530 
531  }
532 
533  return kTRUE;
534 }
535 
536 //________________________________________________________________________
538  //calc subtracted jet mass
539  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
541  else return jet->GetShapeProperties()->GetSecondOrderSubtracted();
542  else
543  return jet->M();
544 }
545 
546 //________________________________________________________________________
548 
549  AliJetContainer *jetCont = GetJetContainer(jetContNb);
550  if (!jet->GetNumberOfTracks())
551  return 0;
552  Double_t den=0.;
553  Double_t num = 0.;
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  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
564  Double_t dr2 = (vp1->Eta()-jet->Eta())*(vp1->Eta()-jet->Eta()) + dphi*dphi;
565  Double_t dr = TMath::Sqrt(dr2);
566  num=num+vp1->Pt()*dr;
567  den=den+vp1->Pt();
568  }
569  return num/den;
570 }
571 
572 //________________________________________________________________________
574 
575  if((fJetShapeSub==kDerivSub) && (jetContNb==0))
578  else
579  return Angularity(jet, jetContNb);
580 
581 }
582 
583 
584 //________________________________________________________________________
586 
587  AliJetContainer *jetCont = GetJetContainer(jetContNb);
588  if (!jet->GetNumberOfTracks())
589  return 0;
590  Double_t den=0.;
591  Double_t num = 0.;
592  AliVParticle *vp1 = 0x0;
593  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
594  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
595 
596  if (!vp1){
597  Printf("AliVParticle associated to constituent not found");
598  continue;
599  }
600 
601  num=num+vp1->Pt()*vp1->Pt();
602  den=den+vp1->Pt();
603  }
604  return TMath::Sqrt(num)/den;
605 }
606 
607 //________________________________________________________________________
609  //calc subtracted jet mass
610  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
612  else return jet->GetShapeProperties()->GetSecondOrderSubtractedpTD();
613  else
614  return PTD(jet, jetContNb);
615 
616 }
617 
618 //_____________________________________________________________________________
620 
621  AliJetContainer *jetCont = GetJetContainer(jetContNb);
622  if (!jet->GetNumberOfTracks())
623  return 0;
624  Double_t mxx = 0.;
625  Double_t myy = 0.;
626  Double_t mxy = 0.;
627  int nc = 0;
628  Double_t sump2 = 0.;
629  Double_t pxjet=jet->Px();
630  Double_t pyjet=jet->Py();
631  Double_t pzjet=jet->Pz();
632 
633 
634  //2 general normalized vectors perpendicular to the jet
635  TVector3 ppJ1(pxjet, pyjet, pzjet);
636  TVector3 ppJ3(- pxjet* pzjet, - pyjet * pzjet, pxjet * pxjet + pyjet * pyjet);
637  ppJ3.SetMag(1.);
638  TVector3 ppJ2(-pyjet, pxjet, 0);
639  ppJ2.SetMag(1.);
640  AliVParticle *vp1 = 0x0;
641  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
642  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
643 
644  if (!vp1){
645  Printf("AliVParticle associated to constituent not found");
646  continue;
647  }
648 
649  TVector3 pp(vp1->Px(), vp1->Py(), vp1->Pz());
650 
651  //local frame
652  TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
653  TVector3 pPerp = pp - pLong;
654  //projection onto the two perpendicular vectors defined above
655 
656  Float_t ppjX = pPerp.Dot(ppJ2);
657  Float_t ppjY = pPerp.Dot(ppJ3);
658  Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
659  if(ppjT<=0) return 0;
660 
661  mxx += (ppjX * ppjX / ppjT);
662  myy += (ppjY * ppjY / ppjT);
663  mxy += (ppjX * ppjY / ppjT);
664  nc++;
665  sump2 += ppjT;}
666 
667  if(nc<2) return 0;
668  if(sump2==0) return 0;
669  // Sphericity Matrix
670  Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
671  TMatrixDSym m0(2,ele);
672 
673  // Find eigenvectors
674  TMatrixDSymEigen m(m0);
675  TVectorD eval(2);
676  TMatrixD evecm = m.GetEigenVectors();
677  eval = m.GetEigenValues();
678  // Largest eigenvector
679  int jev = 0;
680  // cout<<eval[0]<<" "<<eval[1]<<endl;
681  if (eval[0] < eval[1]) jev = 1;
682  TVectorD evec0(2);
683  // Principle axis
684  evec0 = TMatrixDColumn(evecm, jev);
685  Double_t compx=evec0[0];
686  Double_t compy=evec0[1];
687  TVector2 evec(compx, compy);
688  Double_t circ=0;
689  if(jev==1) circ=2*eval[0];
690  if(jev==0) circ=2*eval[1];
691 
692  return circ;
693 
694 
695 
696 }
697 
698 
699 
700 
701 //________________________________________________________________________
703  //calc subtracted jet mass
704 
705  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
708  else
709  return Circularity(jet, jetContNb);
710 
711 }
712 
713 //________________________________________________________________________
715 
716  AliJetContainer *jetCont = GetJetContainer(jetContNb);
717  if (!jet->GetNumberOfTracks())
718  return 0;
719  Double_t den=0.;
720  Double_t num = 0.;
721  AliVParticle *vp1 = 0x0;
722  AliVParticle *vp2 = 0x0;
723  std::vector<int> ordindex;
725  //Printf("Nbof const = %d", jet->GetNumberOfTracks());
726  //Printf("ordindex[0] = %d, ordindex[1] = %d", ordindex[0], ordindex[1]);
727 
728  if(ordindex.size()<2) return -1;
729 
730  vp1 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[0], jetCont->GetParticleContainer()->GetArray()));
731  if (!vp1){
732  Printf("AliVParticle associated to Leading constituent not found");
733  return -1;
734  }
735 
736  vp2 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[1], jetCont->GetParticleContainer()->GetArray()));
737  if (!vp2){
738  Printf("AliVParticle associated to Subleading constituent not found");
739  return -1;
740  }
741 
742 
743  num=vp1->Pt();
744  den=vp2->Pt();
745  //Printf("vp1->Pt() =%f, vp2->Pt() =%f", vp1->Pt(), vp2->Pt());
746 
747 return num-den;
748 }
749 
750 //________________________________________________________________________
752  //calc subtracted jet mass
753 
754  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
757  else
758  return LeSub(jet, jetContNb);
759 
760 }
761 
762 //________________________________________________________________________
764  //calc subtracted jet mass
765 
766  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
769  else
770  return jet->GetNumberOfTracks();
771 
772  }
773 
774 
775 //________________________________________________________________________
776 AliEmcalJetFinder *AliAnalysisTaskEmcalJetShapesMC::Recluster(AliEmcalJet *Jet, Int_t JetContNb, Double_t JetRadius, Double_t SubJetRadius, Double_t SubJetMinPt, Int_t Algorithm, const char* Name){
777 
778  AliJetContainer *JetCont = GetJetContainer(JetContNb);
779  AliEmcalJetFinder *Reclusterer = new AliEmcalJetFinder(Name); //JetFinder Object for reclustered jets
780  Reclusterer->SetRadius(SubJetRadius);
781  Reclusterer->SetJetMinPt(SubJetMinPt);
782  Reclusterer->SetJetAlgorithm(Algorithm); //0 for anti-kt 1 for kt
783  Reclusterer->SetJetMaxEta(0.9-JetRadius);
784  Reclusterer->SetRecombSheme(0);
785 
786 
787  //Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
788  Double_t dVtx[3]={0.,0.,0.};
789  if(Reclusterer->AliEmcalJetFinder::Filter(Jet, JetCont, dVtx)){;} //reclustering jet1 using the jetfinderobject Reclusterer
790  return Reclusterer;
791 }
792 
793 
794 
795 
796 //----------------------------------------------------------------------
798  AliJetContainer *JetCont = GetJetContainer(JetContNb);
799 
800 
801 
802  Double_t SubJetiness_Numerator = 0;
803  Double_t SubJetiness_Denominator = 0;
804  Double_t Index=-2;
805  if (Reclusterer->GetNumberOfJets() < 1) return -2;
806  Index=SubJetOrdering(Jet,Reclusterer,1,0,kTRUE);
807  if(Index==-999) return -2;
808  SubJetiness_Numerator=(Reclusterer->GetJet(Index)->Pt());
809  SubJetiness_Denominator=Jet->Pt();
810  return SubJetiness_Numerator/SubJetiness_Denominator;
811 
812 
813 }
814 //__________________________________________________________________________________
816 
817  AliJetContainer *jetCont = GetJetContainer(jetContNb);
818  if (!jet->GetNumberOfTracks())
819  return 0;
820  Double_t den=0.;
821  Double_t num = 0.;
822  AliVParticle *vp1 = 0x0;
823  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
824  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
825 
826  if (!vp1){
827  Printf("AliVParticle associated to constituent not found");
828  continue;
829  }
830 
831  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
832  Double_t dr2 = (vp1->Eta()-jet->Eta())*(vp1->Eta()-jet->Eta()) + dphi*dphi;
833  Double_t dr = TMath::Sqrt(dr2);
834  if(dr<=fSubjetRadius) num=num+vp1->Pt();
835 
836  }
837  return num/jet->Pt();
838 }
839 
840 
841 
842 
843 //________________________________________________________________________
845  //calc subtracted jet mass
846 
847  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
850  else
851  return CoreFrac(jet, jetContNb);
852 
853 }
854 
855 
856 
857 
858 //----------------------------------------------------------------------
860  AliEmcalJet *SubJet=NULL;
861  Double_t SortingVariable;
862  Int_t ArraySize =N+1;
863  TArrayD *JetSorter = new TArrayD(ArraySize);
864  TArrayD *JetIndexSorter = new TArrayD(ArraySize);
865  for (Int_t i=0; i<ArraySize; i++){
866  JetSorter->SetAt(0,i);
867  }
868  for (Int_t i=0; i<ArraySize; i++){
869  JetIndexSorter->SetAt(0,i);
870  }
871  if(Reclusterer->GetNumberOfJets()<N) return -999;
872  for (Int_t i=0; i<Reclusterer->GetNumberOfJets(); i++){
873  SubJet=Reclusterer->GetJet(i);
874  if (Type==0) SortingVariable=SubJet->Pt();
875  else if (Type==1) SortingVariable=SubJet->E();
876  else if (Type==2) SortingVariable=SubJet->M();
877  for (Int_t j=0; j<N; j++){
878  if (SortingVariable>JetSorter->GetAt(j)){
879  for (Int_t k=N-1; k>=j; k--){
880  JetSorter->SetAt(JetSorter->GetAt(k),k+1);
881  JetIndexSorter->SetAt(JetIndexSorter->GetAt(k),k+1);
882  }
883  JetSorter->SetAt(SortingVariable,j);
884  JetIndexSorter->SetAt(i,j);
885  break;
886  }
887  }
888  }
889  if (!Index) return JetSorter->GetAt(N-1);
890  else return JetIndexSorter->GetAt(N-1);
891 }
892 
893 
894 
895 //returns -1 if the Nth hardest jet is requested where N>number of available jets
896 //type: 0=Pt 1=E 2=M
897 //Index TRUE=returns index FALSE=returns value of quantatiy in question
898 
899 
900 
901 
902 
903 //______________________________________________________________________________
905 
906  AliJetContainer *jetCont = GetJetContainer(jetContNb);
907  if (!jet->GetNumberOfTracks())
908  return 0;
909  Double_t mxx = 0.;
910  Double_t myy = 0.;
911  Double_t mxy = 0.;
912  int nc = 0;
913  Double_t sump2 = 0.;
914 
915  AliVParticle *vp1 = 0x0;
916  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
917  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
918 
919  if (!vp1){
920  Printf("AliVParticle associated to constituent not found");
921  continue;
922  }
923 
924  Double_t ppt=vp1->Pt();
925  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
926 
927  Double_t deta = vp1->Eta()-jet->Eta();
928  mxx += ppt*ppt*deta*deta;
929  myy += ppt*ppt*dphi*dphi;
930  mxy -= ppt*ppt*deta*TMath::Abs(dphi);
931  nc++;
932  sump2 += ppt*ppt;
933 
934  }
935  if(nc<2) return 0;
936  if(sump2==0) return 0;
937  // Sphericity Matrix
938  Double_t ele[4] = {mxx , mxy , mxy , myy };
939  TMatrixDSym m0(2,ele);
940 
941  // Find eigenvectors
942  TMatrixDSymEigen m(m0);
943  TVectorD eval(2);
944  TMatrixD evecm = m.GetEigenVectors();
945  eval = m.GetEigenValues();
946  // Largest eigenvector
947  int jev = 0;
948  // cout<<eval[0]<<" "<<eval[1]<<endl;
949  if (eval[0] < eval[1]) jev = 1;
950  TVectorD evec0(2);
951  // Principle axis
952  evec0 = TMatrixDColumn(evecm, jev);
953  Double_t compx=evec0[0];
954  Double_t compy=evec0[1];
955  TVector2 evec(compx, compy);
956  Double_t sig=0;
957  if(jev==1) sig=TMath::Sqrt(TMath::Abs(eval[0])/sump2);
958  if(jev==0) sig=TMath::Sqrt(TMath::Abs(eval[1])/sump2);
959 
960  return sig;
961 
962 }
963 
964 //________________________________________________________________________
966  //calc subtracted jet mass
967 
968  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
971  else
972  return Sigma2(jet, jetContNb);
973 
974 }
975 
976 
977 //_________________________________________________________________________
979 
980  AliJetContainer *jetCont = GetJetContainer(jetContNb);
981  if (!jet->GetNumberOfTracks())
982  return;
983 
984  Double_t ptJet = jet->Pt();
985 
986  AliVParticle *vp1 = 0x0;
987  std::vector<int> ordindex;
989  //Printf("Nbof const = %d", jet->GetNumberOfTracks());
990  //Printf("ordindex[0] = %d, ordindex[1] = %d", ordindex[0], ordindex[1]);
991  //if(ordindex.size()<2) return -1;
992 
993  for(Int_t iconst =0; iconst<jet->GetNumberOfTracks(); iconst++){
994 
995  vp1 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[iconst], jetCont->GetParticleContainer()->GetArray()));
996  if (!vp1){
997  Printf("AliVParticle associated to Leading constituent not found");
998  return;
999  }
1000 
1001  if (nTFractions[0] <= 0.7*ptJet){
1002  nTFractions[0] += vp1->Pt();
1003  nTFractions[1] +=1;
1004  }
1005 
1006  if (nTFractions[2] <= 0.8*ptJet){
1007  nTFractions[2] += vp1->Pt();
1008  nTFractions[3] +=1;
1009  }
1010 
1011  if (nTFractions[4] <= 0.9*ptJet){
1012  nTFractions[4] += vp1->Pt();
1013  nTFractions[5] +=1;
1014  }
1015 
1016  if (nTFractions[6] <= 0.95*ptJet){
1017  nTFractions[6] += vp1->Pt();
1018  nTFractions[7] +=1;
1019  }
1020  }
1021 }
1022 //_________________________________________________________________________________________________
1024 
1025  //WARNING!!! Only works for parent jets that are clustered with Anti-Kt! To change go to AliEmcalJetFinder.cxx and look at the Nsubjettiness() function
1026 
1027  //Algorithm==0 -> kt_axes;
1028  // Algorithm==1 -> ca_axes;
1029  //Algorithm==2 -> antikt_0p2_axes;
1030  //Algorithm==3 -> wta_kt_axes;
1031  //Algorithm==4 -> wta_ca_axes;
1032  //Algorithm==5 -> onepass_kt_axes;
1033  //Algorithm==6 -> onepass_ca_axes;
1034  //Algorithm==7 -> onepass_antikt_0p2_axes;
1035  //Algorithm==8 -> onepass_wta_kt_axes;
1036  //Algorithm==9 -> onepass_wta_ca_axes;
1037  //Algorithm==10 -> min_axes;
1038 
1039 
1040  //Option==0 returns Nsubjettiness Value
1041  //Option==1 && N==2 returns opening angle between two subjet axes(Delta R?)
1042  //Option==2 && N==2 returns Delta R
1043 
1044  if (Jet->GetNumberOfTracks()>=N){
1045  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1046  AliEmcalJetFinder *JetFinder=new AliEmcalJetFinder("Nsubjettiness");
1047  JetFinder->SetJetMaxEta(0.9-fJetRadius);
1048  JetFinder->SetRadius(fJetRadius);
1049  JetFinder->SetJetAlgorithm(0); //0 for anti-kt 1 for kt //this is for the JET!!!!!!!!!! Not the SubJets
1050  JetFinder->SetRecombSheme(0);
1051  JetFinder->SetJetMinPt(Jet->Pt());
1052  //Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
1053  Double_t dVtx[3]={1,1,1};
1054  //Printf("JetFinder->Nsubjettiness =%f", JetFinder->Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,fSubjetRadius,Beta,Option));
1055  return JetFinder->Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,0.2,Beta,Option,0,Beta_SD,ZCut,SoftDropOn);
1056 
1057  }
1058  else return -2;
1059 }
1060 
1061 
1062 
1063 //________________________________________________________________________
1065 
1067  TClonesArray *tracksArray = partCont->GetArray();
1068 
1069  if(!partCont || !tracksArray) return -99999;
1070  AliAODTrack *track = 0x0;
1071  AliEmcalParticle *emcPart = 0x0;
1072 
1073 
1074  TList *trackList = new TList();
1075  Int_t triggers[100];
1076  for (Int_t iTrigger=0; iTrigger<100; iTrigger++) triggers[iTrigger] = 0;
1077  Int_t iTT = 0;
1078 
1079  for(Int_t iTrack=0; iTrack <= tracksArray->GetEntriesFast(); iTrack++){
1080 
1081 
1082  if (fJetShapeSub == kConstSub){
1083  emcPart = static_cast<AliEmcalParticle*>(tracksArray->At(iTrack));
1084  if (!emcPart) continue;
1085  if(TMath::Abs(emcPart->Eta())>0.9) continue;
1086  if (emcPart->Pt()<0.15) continue;
1087 
1088  if ((emcPart->Pt() >= minpT) && (emcPart->Pt()< maxpT)) {
1089  trackList->Add(emcPart);
1090  triggers[iTT] = iTrack;
1091  iTT++;
1092  }
1093  }
1094  else{
1095  track = static_cast<AliAODTrack*>(tracksArray->At(iTrack));
1096  if (!track) continue;
1097  if(TMath::Abs(track->Eta())>0.9) continue;
1098  if (track->Pt()<0.15) continue;
1099  if (!(track->TestFilterBit(768))) continue;
1100 
1101  if ((track->Pt() >= minpT) && (track->Pt()< maxpT)) {
1102  trackList->Add(track);
1103  triggers[iTT] = iTrack;
1104  iTT++;
1105 
1106  }
1107  }
1108  }
1109 
1110  if (iTT == 0) return -99999;
1111  Int_t nbRn = 0, index = 0 ;
1112  TRandom3* random = new TRandom3(0);
1113  nbRn = random->Integer(iTT);
1114 
1115  index = triggers[nbRn];
1116  //Printf("iTT Total= %d, nbRn = %d, Index = %d",iTT, nbRn, index );
1117  return index;
1118 
1119 }
1120 
1121 //__________________________________________________________________________________
1123 
1124  if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1125  else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1126  if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1127  else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1128  double dphi = mphi-vphi;
1129  if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1130  else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1131  return dphi;//dphi in [-Pi, Pi]
1132 }
1133 
1134 
1135 //________________________________________________________________________
1137  //
1138  // retrieve event objects
1139  //
1141  return kFALSE;
1142 
1143  return kTRUE;
1144 }
1145 
1146 //_______________________________________________________________________
1148 {
1149  // Called once at the end of the analysis.
1150 
1151  AliInfo("Terminate");
1152  AliAnalysisTaskSE::Terminate();
1153 
1154  fOutput = dynamic_cast<AliEmcalList*> (GetOutputData(1));
1155  if (!fOutput) {
1156  AliError("fOutput not available");
1157  return;
1158  }
1159 
1160  fTreeObservableTagging = dynamic_cast<TTree*>(GetOutputData(2));
1161  if (!fTreeObservableTagging){
1162  Printf("ERROR: fTreeObservableTagging not available");
1163  return;
1164  }
1165 
1166 }
1167 
1168 //_________________________________________________________________________
1170 
1171  std::vector<fastjet::PseudoJet> fInputVectors;
1172  fInputVectors.clear();
1173  Double_t JetInvMass=0, PseudJetInvMass=0, TrackMom = 0, TrackEnergy = 0;
1174  fastjet::PseudoJet PseudoTracks;
1175  fastjet::PseudoJet MyJet;
1176  fastjet::PseudoJet PseudoTracksCMS;
1177  AliParticleContainer *fTrackCont = fJetCont->GetParticleContainer();
1178  //cout<<"CALL TO SOFTDROP"<<endl;
1179 
1180  Double_t zeta=0;
1181  Double_t angle=0;
1182  if (fTrackCont) for (Int_t i=0; i<fJet->GetNumberOfTracks(); i++) {
1183  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
1184  if (!fTrk) continue;
1185  JetInvMass += fTrk->M();
1186 
1187  PseudoTracks.reset(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
1188  TrackMom += TMath::Sqrt(TMath::Power(fTrk->Px(),2)+TMath::Power(fTrk->Py(),2)+TMath::Power(fTrk->Pz(),2));
1189  TrackEnergy += fTrk->E();
1190  PseudoTracks.set_user_index(fJet->TrackAt(i)+100);
1191  PseudJetInvMass += PseudoTracks.m();
1192  fInputVectors.push_back(PseudoTracks);
1193 
1194  }
1195 
1196  //fRandom->SetSeed(0);
1197  //here add N tracks with random phi and eta and theta according to bdmps distrib.
1198 
1199  MyJet.reset(fJet->Px(),fJet->Py(),fJet->Pz(),fJet->E());
1200  Double_t omegac=0.5*fqhat*fxlength*fxlength/0.2;
1201  Double_t thetac=TMath::Sqrt(12*0.2/(fqhat*TMath::Power(fxlength,3)));
1202  Double_t xQs=TMath::Sqrt(fqhat*fxlength);
1203  //cout<<"medium parameters "<<omegac<<" "<<thetac<<" "<<xQs<<endl;
1204 
1205  for(Int_t i=0;i<fAdditionalTracks;i++){
1206 
1207  Double_t ppx,ppy,ppz,kTscale,lim2o,lim1o;
1208  Double_t lim2=xQs;
1209  Double_t lim1=10000;
1210 
1211  //generation of kT according to 1/kT^4, with minimum QS=2 GeV and maximum ~sqrt(ptjet*T)
1212  fTf1Kt= new TF1("fTf1Kt","1/(x*x*x*x)",lim2,lim1);
1213  kTscale=fTf1Kt->GetRandom();
1214  //generation within the jet cone
1215 
1216  //generation of w according to 1/w, with minimum wc
1217  //omega needs to be larger than kT so to have well defined angles
1218  lim2o=kTscale;
1219  lim1o=kTscale/TMath::Sin(0.1);
1220  fTf1Omega= new TF1("fTf1Omega","1/x",lim2o,lim1o);
1221  Double_t omega=fTf1Omega->GetRandom();
1222 
1223  Double_t sinpptheta=kTscale/omega;
1224  Double_t pptheta=TMath::ASin(sinpptheta);
1225  //cout<<"angle_omega_kt"<<pptheta<<" "<<omega<<" "<<kTscale<<endl;
1226  if(pptheta>fJetRadius) continue;
1227 
1228  PseudoTracksCMS.reset(kTscale/TMath::Sqrt(2),kTscale/TMath::Sqrt(2),omega*TMath::Cos(pptheta),omega);
1229  //boost the particle in the rest frame of the jet to the lab frame
1230  fastjet::PseudoJet PseudoTracksLab=PseudoTracksCMS.boost(MyJet);
1231  PseudoTracksLab.set_user_index(i+fJet->GetNumberOfTracks()+100);
1232  fInputVectors.push_back(PseudoTracksLab);
1233  //in the frame of the jet
1234  zeta=omega/fJet->E();
1235  angle=pptheta;
1236 
1237 
1238  }
1239 
1240 
1241 
1242 
1243 
1244 
1245 
1246 
1247 
1248  fastjet::JetDefinition fJetDef(fastjet::antikt_algorithm, fJetRadius*2, static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
1249  try {
1250  fastjet::ClusterSequence fClustSeqSA(fInputVectors, fJetDef);
1251  std::vector<fastjet::PseudoJet> fOutputJets;
1252  fOutputJets.clear();
1253  fOutputJets=fClustSeqSA.inclusive_jets(0);
1254 
1255  //cout<<fOutputJets[0].perp()<<" "<<fJet->Pt()<<endl;
1256 
1257  fastjet::contrib::SoftDrop softdrop(beta, zcut);
1258 
1259  softdrop.set_verbose_structure(kTRUE);
1260  fastjet::JetAlgorithm jetalgo(fastjet::cambridge_algorithm);
1261  if(ReclusterAlgo==2) jetalgo=fastjet::antikt_algorithm;
1262  if(ReclusterAlgo==1) jetalgo=fastjet::kt_algorithm;
1263  if(ReclusterAlgo==0) jetalgo=fastjet::cambridge_algorithm;
1264 
1265 #if FASTJET_VERSION_NUMBER >= 30302
1266  fastjet::Recluster recluster(jetalgo,1,fastjet::Recluster::keep_only_hardest);
1267 #else
1268  fastjet::contrib::Recluster recluster(jetalgo,1,true);
1269 #endif
1270  softdrop.set_reclustering(true, &recluster);
1271  fastjet::PseudoJet finaljet = softdrop(fOutputJets[0]);
1272  Int_t NDroppedTracks = fJet->GetNumberOfTracks()-finaljet.constituents().size();
1273 
1274  Double_t SymParam, Mu, DeltaR, GroomedPt,GroomedMass;
1275  Int_t NGroomedBranches;
1276  SymParam=(finaljet.structure_of<fastjet::contrib::SoftDrop>().symmetry());
1277  Mu=(finaljet.structure_of<fastjet::contrib::SoftDrop>().mu());
1278  DeltaR=(finaljet.structure_of<fastjet::contrib::SoftDrop>().delta_R());
1279  NGroomedBranches=finaljet.structure_of<fastjet::contrib::SoftDrop>().dropped_count();
1280  GroomedPt=finaljet.perp();
1281  GroomedMass=finaljet.m();
1282 
1283 
1284  fShapesVar[11]=SymParam;
1285  fShapesVar[12]=DeltaR;
1286  // fShapesVar[14]=zeta;
1287  // fShapesVar[15]=angle;
1288  // fShapesVar[16]=GroomedMass;}
1289  // if(ReclusterAlgo==1){
1290  // fShapesVar[17]=SymParam;
1291  // fShapesVar[18]=DeltaR;
1292  // fShapesVar[19]=zeta;
1293  // fShapesVar[20]=angle;
1294  // fShapesVar[21]=GroomedMass; }
1295 
1296  // if(ReclusterAlgo==2){
1297  // fShapesVar[22]=SymParam;
1298  // fShapesVar[23]=DeltaR;
1299  // fShapesVar[24]=zeta;
1300  // fShapesVar[25]=angle;
1301  // fShapesVar[26]=GroomedMass;
1302  // }}
1303  // if(beta==1){
1304  // fShapesVar[27]=SymParam;
1305  // fShapesVar[28]=DeltaR;
1306  // fShapesVar[29]=zeta;
1307  // fShapesVar[30]=angle;
1308  // fShapesVar[31]=GroomedMass;
1309  // }
1310  // //this one kills soft and large angle radiation
1311  // if((beta==1.5) && (zcut==0.5)){
1312  // fShapesVar[32]=SymParam;
1313  // fShapesVar[33]=DeltaR;
1314  // fShapesVar[34]=zeta;
1315  // fShapesVar[35]=angle;
1316  // fShapesVar[36]=GroomedMass; }
1317  // //this option favour democratic branches at large kt
1318  // if((beta==-1) && (zcut==0.005)){
1319  // fShapesVar[37]=SymParam;
1320  // fShapesVar[38]=DeltaR;
1321  // fShapesVar[39]=zeta;
1322  // fShapesVar[40]=angle;
1323  // fShapesVar[41]=GroomedMass; }
1324 
1325  // if((beta==-2) && (zcut==0.005)){
1326  // fShapesVar[42]=SymParam;
1327  // fShapesVar[43]=DeltaR;
1328  // fShapesVar[44]=zeta;
1329  // fShapesVar[45]=angle;
1330  // fShapesVar[46]=GroomedMass; }
1331 
1332 
1333 
1334 
1335 
1336 
1337 
1338 
1339 
1340  } catch (fastjet::Error) {
1341  AliError(" [w] FJ Exception caught.");
1342  //return -1;
1343  }
1344  if(fTf1Kt){ delete fTf1Kt;}
1345  if(fTf1Omega){ delete fTf1Omega;}
1346 }
1347 
1348 
1349 //_________________________________________________________________________
1351 
1352  std::vector<fastjet::PseudoJet> fInputVectors;
1353  fInputVectors.clear();
1354  fastjet::PseudoJet PseudoTracks;
1355  fastjet::PseudoJet PseudoTracksLab;
1356  double xflagalgo=0;
1357  double lnpt_relinject=0;
1358  double yinject=0;
1359  int xflagAdded=0;
1360  double zinject,angleinject,pptheta,sinpptheta,omega,omega2,angle2;
1361  AliParticleContainer *fTrackCont = fJetCont->GetParticleContainer();
1362  Float_t pTscale=0., phiscale=0., thetascale=0., pXscale=0., pYscale=0., pZscale=0., pscale=0.;
1363 
1364  if (fTrackCont) for (Int_t i=0; i<fJet->GetNumberOfTracks(); i++) {
1365  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
1366  if (!fTrk) continue;
1367  if (fScaleELoss){
1368  pTscale = xfraction*sqrt(pow(fTrk->Px(),2)+pow(fTrk->Py(),2));
1369  phiscale = fTrk->Phi();
1370  thetascale = 2.*TMath::ATan(TMath::Exp(-1.*(fTrk->Eta())));
1371  pXscale = pTscale * TMath::Cos(phiscale);
1372  pYscale = pTscale * TMath::Sin(phiscale);
1373  pZscale = pTscale/TMath::Tan(thetascale);
1374  pscale = TMath::Sqrt(pTscale*pTscale+pZscale*pZscale);
1375  PseudoTracks.reset(pXscale, pYscale, pZscale, pscale);
1376  }
1377  else PseudoTracks.reset(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
1378  PseudoTracks.set_user_index(fJet->TrackAt(i)+100);
1379  fInputVectors.push_back(PseudoTracks);
1380 
1381  }
1382  /*if(fAddMedScat){
1383  for(int i = 0; i < fAddMedScatN; i++){
1384  TRandom3 rand1(0),rand2(0),rand3(0); //set range +- jet R
1385  Double_t randN1 = 0.4*0.4*rand1.Rndm();
1386  Double_t randN2 = 2*TMath::Pi()*rand2.Rndm();
1387  Double_t phi_rand = (fJet->Phi())+TMath::Sqrt(randN1)*TMath::Sin(randN2);
1388  Double_t eta_rand = (fJet->Eta())+TMath::Sqrt(randN1)*TMath::Cos(randN2);
1389  Double_t fAddMedScatPt = (fAddMedScatPtFrac*fJet->Pt())/fAddMedScatN;
1390  PseudoTracks.reset(fAddMedScatPt*TMath::Cos(phi_rand),fAddMedScatPt*TMath::Sin(phi_rand),fAddMedScatPt/TMath::Tan(eta_rand),fAddMedScatPt);
1391  PseudoTracks.set_user_index(i+fJet->GetNumberOfTracks()+100);
1392  fInputVectors.push_back(PseudoTracks);
1393  }
1394  }*/
1395  if(fAddMedScat){
1396  for(int i = 0; i < fAddMedScatN; i++){
1397  Double_t ppx,ppy,ppz,SoftkTscale,lim2o,lim1o;
1398  Double_t lim1=0.1;
1399  Double_t lim2=0.5;
1400  fTf1SoftKt= new TF1("fTf1SoftKt","1/(x)",lim1,lim2);
1401  SoftkTscale=fTf1SoftKt->GetRandom();
1402 
1403  lim2o=SoftkTscale;
1404  lim1o=SoftkTscale/TMath::Sin(0.1);
1405  fTf1SoftOmega= new TF1("fTf1SoftOmega","1/x",lim2o,lim1o);
1406  omega=fTf1SoftOmega->GetRandom();
1407  sinpptheta=SoftkTscale/omega;
1408  pptheta=TMath::ASin(sinpptheta);
1409  if(pptheta>fJetRadius) continue;
1410 
1411  TLorentzVector pTrackCMS(SoftkTscale/TMath::Sqrt(2),SoftkTscale/TMath::Sqrt(2),omega*TMath::Cos(pptheta),omega);
1412  TVector3 MyJet(fJet->Px(),fJet->Py(),fJet->Pz());
1413  TVector3 direction = MyJet.Unit();
1414  //rotate the track to the jet frame
1415  pTrackCMS.RotateUz(direction);
1416 
1417  //add the rotated track to the jet
1418  PseudoTracksLab.reset(pTrackCMS.Px(),pTrackCMS.Py(),pTrackCMS.Pz(),pTrackCMS.E());
1419 
1420  PseudoTracksLab.set_user_index(i+fJet->GetNumberOfTracks()+100);
1421 
1422  omega2=PseudoTracksLab.perp();
1423  angle2=pTrackCMS.Angle(MyJet);
1424 
1425  fInputVectors.push_back(PseudoTracksLab);
1426  }
1427  }
1428 
1429 
1430  //add tracks to the jet prior to the reclusterer in case of iterative mapping of splittings
1431 
1432  Double_t omegac=0.5*fqhat*fxlength*fxlength/0.2;
1433  Double_t thetac=TMath::Sqrt(12*0.2/(fqhat*TMath::Power(fxlength,3)));
1434  Double_t xQs=TMath::Sqrt(fqhat*fxlength);
1435 
1436 
1437  for(Int_t i=0;i<fAdditionalTracks;i++){
1438 
1439  Double_t ppx,ppy,ppz,kTscale,lim2o,lim1o;
1440  Double_t lim2=xQs;
1441  Double_t lim1=10000;
1442 
1443  //generation of kT according to 1/kT^4, with minimum QS=2 GeV and maximum ~sqrt(ptjet*T)
1444  fTf1Kt= new TF1("fTf1Kt","1/(x*x*x*x)",lim2,lim1);
1445  kTscale=fTf1Kt->GetRandom();
1446  //generation within the jet cone
1447 
1448  //generation of w according to 1/w, with minimum wc
1449  //omega needs to be larger than kT so to have well defined angles
1450  lim2o=kTscale;
1451  lim1o=kTscale/TMath::Sin(0.1);
1452  fTf1Omega= new TF1("fTf1Omega","1/x",lim2o,lim1o);
1453  omega=fTf1Omega->GetRandom();
1454  sinpptheta=kTscale/omega;
1455  pptheta=TMath::ASin(sinpptheta);
1456  if(pptheta>fJetRadius) continue;
1457 
1458  //Lorentz vector in the frame where the jet moves along z axis
1459  TLorentzVector pTrackCMS(kTscale/TMath::Sqrt(2),kTscale/TMath::Sqrt(2),omega*TMath::Cos(pptheta),omega);
1460  TVector3 MyJet(fJet->Px(),fJet->Py(),fJet->Pz());
1461  TVector3 direction = MyJet.Unit();
1462  //rotate the track to the jet frame
1463  pTrackCMS.RotateUz(direction);
1464 
1465  //add the rotated track to the jet
1466  PseudoTracksLab.reset(pTrackCMS.Px(),pTrackCMS.Py(),pTrackCMS.Pz(),pTrackCMS.E());
1467 
1468  PseudoTracksLab.set_user_index(i+fJet->GetNumberOfTracks()+100);
1469 
1470  omega2=PseudoTracksLab.perp();
1471  angle2=pTrackCMS.Angle(MyJet);
1472 
1473  fInputVectors.push_back(PseudoTracksLab);
1474  xflagAdded=1;
1475  }
1476 
1477 
1478 
1479 
1480  fastjet::JetAlgorithm jetalgo(fastjet::antikt_algorithm);
1481  if(ReclusterAlgo==0){ xflagalgo=0.5;
1482  jetalgo=fastjet::kt_algorithm ;}
1483 
1484  if(ReclusterAlgo==1){ xflagalgo=1.5;
1485  jetalgo=fastjet::cambridge_algorithm;}
1486  if(ReclusterAlgo==2){ xflagalgo=2.5;
1487  jetalgo=fastjet::antikt_algorithm;}
1488 
1489  fastjet::JetDefinition fJetDef(jetalgo, 1., static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
1490 
1491  try {
1492  fastjet::ClusterSequence fClustSeqSA(fInputVectors, fJetDef);
1493  std::vector<fastjet::PseudoJet> fOutputJets;
1494  fOutputJets.clear();
1495  fOutputJets=fClustSeqSA.inclusive_jets(0);
1496 
1497  fastjet::PseudoJet jj;
1498  fastjet::PseudoJet j1;
1499  fastjet::PseudoJet j2;
1500  jj=fOutputJets[0];
1501  double ndepth=0;
1502  double nsd=0;
1503  double nall=0;
1504  fShapesVar[10]=fOutputJets[0].perp();
1505  while(jj.has_parents(j1,j2)){
1506 
1507  if(j1.perp() < j2.perp()) swap(j1,j2);
1508  double delta_R=j1.delta_R(j2);
1509  double z=j2.perp()/(j1.perp()+j2.perp());
1510  double y =log(1.0/delta_R);
1511  double lnpt_rel=log(z*delta_R);
1512  double lnkt=log(jj.perp()*z*delta_R);
1513 
1514  nall=nall+1;
1515  if(z>fHardCutoff){
1516  ndepth=ndepth+1;
1517  nsd=nsd+1;
1518  if(nsd==1 && ReclusterAlgo==1) fShapesVar[9]=z;
1519  Double_t LundEntries[6] = {y,lnpt_rel,fOutputJets[0].perp(),xflagalgo,partonFlavor,ndepth};
1520  Double_t LundEntries_kt[6] = {y,lnkt,fOutputJets[0].perp(),xflagalgo,partonFlavor,ndepth};
1521  fHLundIterative->Fill(LundEntries);
1522  fHLundIterative_ktaxis->Fill(LundEntries_kt);
1523  }
1524  jj=j1;}
1525 
1526  if(ReclusterAlgo==1){
1527  fShapesVar[6]=nsd;
1528  fShapesVar[7]=nall;}
1529 
1530  if(fAdditionalTracks>0 && xflagAdded>0){
1531  zinject=omega2/fOutputJets[0].perp();
1532  angleinject=angle2;
1533 
1534  yinject =log(1.0/angleinject);
1535  lnpt_relinject=log(zinject*angleinject);
1536  Double_t LundEntriesInject[4] = {yinject,lnpt_relinject,fOutputJets[0].perp(),fJet->Pt()};
1537  fHLundIterativeInject->Fill(LundEntriesInject);}
1538 
1539 
1540  } catch (fastjet::Error) {
1541  AliError(" [w] FJ Exception caught.");
1542  //return -1;
1543  }
1544 
1545 
1546  if(fTf1SoftKt){ delete fTf1SoftKt;}
1547  if(fTf1SoftOmega){ delete fTf1SoftOmega;}
1548  if(fTf1Kt){ delete fTf1Kt;}
1549  if(fTf1Omega){ delete fTf1Omega;}
1550 
1551 
1552  return;
1553 
1554 
1555 }
1556 
1557 
1559  const Double_t jetradius,
1560  const Double_t subjetradius,
1561  const char *ntracksPartLevel,
1562  const char *type,
1563  const char *CentEst,
1564  Int_t pSel,
1565  TString trigClass,
1566  TString kEmcalTriggers,
1567  TString tag,
1568  const char *rhoName,
1572  Float_t minpTHTrigger, Float_t maxpTHTrigger,
1574 
1575 
1576 
1577  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1578  if (!mgr)
1579  {
1580  ::Error("AliAnalysisTaskEmcalJetShapesMC","No analysis manager found.");
1581  return 0;
1582  }
1583  Bool_t ismc=kFALSE;
1584  ismc = (mgr->GetMCtruthEventHandler())?kTRUE:kFALSE;
1585 
1586  // Check the analysis type using the event handlers connected to the analysis manager.
1587  //==============================================================================
1588  if (!mgr->GetInputEventHandler())
1589  {
1590  ::Error("AddTaskJetShapesMC", "This task requires an input event handler");
1591  return NULL;
1592  }
1593 
1594  TString wagonName1 = Form("JetShapesMC_%s_Histos%s%s",njetsBase,trigClass.Data(),tag.Data());
1595  TString wagonName2 = Form("JetShapesMC_%s_Tree%s%s",njetsBase,trigClass.Data(),tag.Data());
1596 
1597  //Configure jet tagger task
1599 
1600  //task->SetNCentBins(4);
1601  task->SetJetShapeType(jetShapeType);
1602  task->SetJetShapeSub(jetShapeSub);
1603  task->SetJetSelection(jetSelection);
1604  task->SetDerivativeSubtractionOrder(derivSubtrOrder);
1605  task->SetJetRadius(jetradius);
1606  task->SetSubjetRadius(subjetradius);
1607 
1608  if (jetSelection == AliAnalysisTaskEmcalJetShapesMC::kRecoil) task->SetPtTriggerSelections(minpTHTrigger, maxpTHTrigger);
1609 
1610  TString thename(njetsBase);
1611 
1612  AliParticleContainer *trackContPartLevel = task->AddMCParticleContainer(ntracksPartLevel);
1613 
1614  AliJetContainer *jetContBase=0x0;
1615  TString strType(type);
1616 
1617  jetContBase = task->AddJetContainer(njetsBase,strType,jetradius);
1618  if(jetContBase) {
1619  jetContBase->SetRhoName(rhoName);
1620  jetContBase->ConnectParticleContainer(trackContPartLevel);
1621  //jetContBase->ConnectClusterContainer(clusterCont);
1622  jetContBase->SetPercAreaCut(0.6);
1623  }
1624 
1625  task->SetCaloTriggerPatchInfoName(kEmcalTriggers.Data());
1626  task->SetCentralityEstimator(CentEst);
1627  task->SelectCollisionCandidates(pSel);
1628  task->SetUseAliAnaUtils(kFALSE);
1629 
1630  mgr->AddTask(task);
1631 
1632  //Connnect input
1633  mgr->ConnectInput (task, 0, mgr->GetCommonInputContainer() );
1634 
1635  //Connect output
1636  TString contName1(wagonName1);
1637  TString contName2(wagonName2);
1638 
1639  if (jetShapeType == AliAnalysisTaskEmcalJetShapesMC::kGenShapes) {
1640  contName1 += "_GenShapes";
1641  contName2 += "_GenShapes";
1642  }
1643 
1644  switch (jetShapeSub) {
1645 
1647  contName1 += "_NoSub";
1648  contName2 += "_NoSub";
1649  break;
1650 
1652  contName1 += "_ConstSub";
1653  contName2 += "_ConstSub";
1654  break;
1655 
1657  contName1 += "_DerivSub";
1658  contName2 += "_DerivSub";
1659  break;
1660 
1661 
1662  }
1663 
1664  switch (jetSelection) {
1666  contName1 += "_Incl";
1667  contName2 += "_Incl";
1668  break;
1669 
1670 
1672  TString recoilTriggerString = Form("_Recoil_%.0f_%0.f", minpTHTrigger, maxpTHTrigger);
1673  contName1 += recoilTriggerString;
1674  contName2 += recoilTriggerString;
1675 
1676  break;
1677 
1678  }
1679 
1680 
1681  TString outputfile = Form("%s",AliAnalysisManager::GetCommonFileName());
1682 
1683 
1684  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contName1.Data(),TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
1685  AliAnalysisDataContainer *coutput2 = mgr->CreateContainer(contName2.Data(),TTree::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
1686 
1687  mgr->ConnectOutput(task,1,coutput1);
1688  mgr->ConnectOutput(task,2,coutput2);
1689 
1690  return task;
1691 
1692 }
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:121
Double_t FjNSubJettiness(AliEmcalJet *Jet, Int_t JetContNb, Int_t N, Int_t Algorithm, Double_t Beta, Int_t Option, Double_t Beta_SD=0.0, Double_t ZCut=0.1, Int_t SoftDropOn=0)
Double_t Py() const
Definition: AliEmcalJet.h:107
Double_t Phi() const
Definition: AliEmcalJet.h:117
Float_t GetPythiaEventWeight() const
void RecursiveParents(AliEmcalJet *fJet, AliJetContainer *fJetCont, Int_t ReclusterAlgo, Float_t PartonFlavor)
Float_t GetPartonPhi7() const
void SetPercAreaCut(Float_t p)
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:119
TF1 * fTf1SoftOmega
to generate kT according to BDMS tail
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:160
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:139
Double_t Px() const
Definition: AliEmcalJet.h:106
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.
Int_t GetPartonFlag6() const
void SetRhoName(const char *n)
AliParticleContainer * GetParticleContainer() const
Double_t GetFirstOrderSubtractedConstituent() const
void SetRecombSheme(Int_t val)
static AliAnalysisTaskEmcalJetShapesMC * AddTaskJetShapesMC(const char *njetsBase, const Double_t jetradius, const Double_t subjetradius, const char *ntracksPartLevel, const char *type, const char *CentEst, Int_t pSel, TString trigClass="", TString kEmcalTriggers="", TString tag="", const char *rhoName="", AliAnalysisTaskEmcalJetShapesMC::JetShapeType jetShapeType=AliAnalysisTaskEmcalJetShapesMC::kGenShapes, AliAnalysisTaskEmcalJetShapesMC::JetShapeSub jetShapeSub=AliAnalysisTaskEmcalJetShapesMC::kNoSub, AliAnalysisTaskEmcalJetShapesMC::JetSelectionType jetSelection=AliAnalysisTaskEmcalJetShapesMC::kInclusive, Float_t minpTHTrigger=0., Float_t maxpTHTrigger=0., AliAnalysisTaskEmcalJetShapesMC::DerivSubtrOrder derivSubtrOrder=AliAnalysisTaskEmcalJetShapesMC::kSecondOrder)
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
Float_t fqhat
Random number generator.
Double_t DeltaR(const AliVParticle *part1, const AliVParticle *part2)
Task to store and correlate the MC shapes.
TF1 * fTf1SoftKt
to generate omega for soft background
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)
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
TClonesArray * GetArray() const
Float_t PTD(AliEmcalJet *jet, Int_t jetContNb=0)
void SetJetMaxEta(Double_t val)
AliEmcalJet * GetNextAcceptJet()
TF1 * fTf1Kt
to generate omega according to BDMPS tail
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)
void ConnectParticleContainer(AliParticleContainer *c)
Double_t Pt() const
Definition: AliEmcalJet.h:109
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)
Double_t GetSecondOrderSubtractedCircularity() const
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)
Enable general histograms.
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
const AliEmcalPythiaInfo * GetPythiaInfo() const
Double_t Pz() const
Definition: AliEmcalJet.h:108
Float_t Circularity(AliEmcalJet *jet, Int_t jetContNb=0)
void SoftDrop(AliEmcalJet *fJet, AliJetContainer *fJetCont, Double_t zcut, Double_t beta, Int_t ReclusterAlgo)
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:375
Float_t GetPartonPt7() const
Float_t LeSub(AliEmcalJet *jet, Int_t jetContNb=0)
void swap(PWGJE::EMCALJetTasks::AliAnalysisTaskEmcalJetHPerformance &first, PWGJE::EMCALJetTasks::AliAnalysisTaskEmcalJetHPerformance &second)
Double_t M() const
Definition: AliEmcalJet.h:120
void ResetCurrentID(Int_t i=-1)
Reset the iterator to a given index.
Container for jet within the EMCAL jet framework.
std::vector< int > GetPtSortedTrackConstituentIndexes(TClonesArray *tracks) const
Float_t GetPartonPt6() const