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