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