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