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