AliPhysics  b0c77bb (b0c77bb)
AliAnalysisTaskRecoilJetYield.cxx
Go to the documentation of this file.
1 //
2 // My analysis task
3 //
4 // Author: Harry Andrews
5 
6 #include <TClonesArray.h>
7 #include <TH1F.h>
8 #include <TH2F.h>
9 #include <TH3F.h>
10 #include <TCanvas.h>
11 #include <THnSparse.h>
12 #include <TTree.h>
13 #include <TList.h>
14 #include <TLorentzVector.h>
15 #include <TProfile.h>
16 #include <TChain.h>
17 #include <TSystem.h>
18 #include <TFile.h>
19 #include <TKey.h>
20 #include <AliAnalysisDataSlot.h>
21 #include <AliAnalysisDataContainer.h>
22 #include "TMatrixD.h"
23 #include "TMatrixDSym.h"
24 #include "TMatrixDSymEigen.h"
25 #include "TVector3.h"
26 #include "TVector2.h"
27 #include "AliVCluster.h"
28 #include "AliVTrack.h"
29 #include "AliEmcalJet.h"
30 #include "AliRhoParameter.h"
31 #include "AliLog.h"
32 #include "AliEmcalParticle.h"
33 #include "AliMCEvent.h"
34 #include "AliGenPythiaEventHeader.h"
35 #include "AliAODMCHeader.h"
36 #include "AliMCEvent.h"
37 #include "AliAnalysisManager.h"
38 #include "AliJetContainer.h"
39 #include "AliParticleContainer.h"
40 //#include "AliPythiaInfo.h"
41 #include "TRandom3.h"
42 #include "AliPicoTrack.h"
43 #include "AliEmcalJetFinder.h"
44 #include "AliAODEvent.h"
46 #include "AliMultSelection.h"
47 
48 //Globals
49 
50 #include <fastjet/config.h>
51 #if FASJET_VERSION_NUMBER >= 30302
52 #include <fastjet/tools/Recluster.hh>
53 #else
54 #include <fastjet/contrib/Recluster.hh>
55 #endif
56 
57 
61 Double_t Phi_Upper=2*(TMath::Pi());
62 Double_t Phi_Lower=(-1*(TMath::Pi()));
65 
66 
67 
68 
69 using std::cout;
70 using std::endl;
71 
73 
74 //________________________________________________________________________
77  fContainer(0),
78  fMinFractionShared(0),
79  fJetShapeType(kData),
80  fJetShapeSub(kNoSub),
81  fJetSelection(kInclusive),
82  fPtThreshold(-9999.),
83  fRMatching(0.2),
84  fPtMinTriggerHadron(20.),
85  fPtMaxTriggerHadron(50.),
86  fRecoilAngularWindow(0.6),
87  fSemigoodCorrect(0),
88  fHolePos(0),
89  fHoleWidth(0),
90  fCentSelectOn(kTRUE),
91  fCentMin(0),
92  fCentMax(10),
93  fSubJetAlgorithm(0),
94  fSubJetRadius(0.1),
95  fSubJetMinPt(1),
96  fJetRadius(0.4),
97  fRMatched(0.2),
98  fSharedFractionPtMin(0.5),
99  fFullTree(kFALSE),
100  fBeta_SD(0),
101  fZCut(0.1),
102  fRho(1e-6),
103  fRhoParam(0),
104  fNsubMeasure(kFALSE),
105  fDoSoftDrop(kFALSE),
106  fhJetPt(0x0),
107  fhJetPhi(0x0),
108  fhJetEta(0x0),
109  fhJetMass(0x0),
110  fhJetRadius(0x0),
111  fhJetCounter(0x0),
112  fhNumberOfJetTracks(0x0),
113  fTreeJetInfo(0),
114  fhJetArea(0x0),
115  fhTrackPt(0x0),
116  fhTrackEta(0x0),
117  fhTrackPhi(0x0),
118  fhGroomedPtvJetPt(0x0),
119  fhDroppedBranches(0x0),
120  fhPtTriggerHadron(0x0),
121  fh2PtTriggerHadronJet(0x0),
122  fhPhiTriggerHadronJet(0x0),
123  fhPhiTriggerHadronEventPlane(0x0),
124  fhPhiTriggerHadronEventPlaneTPC(0x0),
125  fhDetJetPt_Incl(0x0),
126  fhDetJetPt_Matched(0x0),
127  fhJetPt_Det(0x0),
128  fhJetPt_True(0x0),
129  fhJetPhi_Det(0x0),
130  fhJetPhi_True(0x0),
131  fhJetEta_Det(0x0),
132  fhJetEta_True(0x0),
133  fhJetMass_Det(0x0),
134  fhJetMass_True(0x0),
135  fhJetRadius_Det(0x0),
136  fhJetRadius_True(0x0),
137  fhJetCounter_Det(0x0),
138  fhJetCounter_True(0x0),
139  fhNumberOfJetTracks_Det(0x0),
140  fhNumberOfJetTracks_True(0x0),
141  fh2PtRatio(0x0),
142  fhLundIterative(0x0),
143  fhLundIterativeTrue(0x0),
144  fReclusterAlgo(0),
145  fSubMatching(kFALSE),
146  bMinSubjetPt(kFALSE),
147  fMinSubjetPt(0)
148 
149 
150 
151 {
152  for(Int_t i=0;i<nBranch;i++){
153  fJetInfoVar[i]=0;
154  }
155  SetMakeGeneralHistograms(kTRUE);
156  DefineOutput(1, TList::Class());
157  DefineOutput(2, TTree::Class());
158 }
159 
160 //________________________________________________________________________
162  AliAnalysisTaskEmcalJet(name, kTRUE),
163  fContainer(0),
164  fMinFractionShared(0),
165  fJetShapeType(kData),
166  fJetShapeSub(kNoSub),
167  fJetSelection(kInclusive),
168  fPtThreshold(-9999.),
169  fRMatching(0.2),
170  fPtMinTriggerHadron(20.),
171  fPtMaxTriggerHadron(50.),
172  fRecoilAngularWindow(0.6),
173  fSemigoodCorrect(0),
174  fHolePos(0),
175  fHoleWidth(0),
176  fCentSelectOn(kTRUE),
177  fCentMin(0),
178  fCentMax(10),
179  fSubJetAlgorithm(0),
180  fSubJetRadius(0.1),
181  fSubJetMinPt(1),
182  fJetRadius(0.4),
183  fRMatched(0.2),
184  fSharedFractionPtMin(0.5),
185  fFullTree(kFALSE),
186  fBeta_SD(0),
187  fZCut(0.1),
188  fRho(1e-6),
189  fRhoParam(0),
190  fNsubMeasure(kFALSE),
191  fDoSoftDrop(kFALSE),
192  fhJetPt(0x0),
193  fhJetPhi(0x0),
194  fhJetEta(0x0),
195  fhJetMass(0x0),
196  fhJetRadius(0x0),
197  fhJetCounter(0x0),
198  fhNumberOfJetTracks(0x0),
199  fTreeJetInfo(0),
200  fhJetArea(0x0),
201  fhTrackPt(0x0),
202  fhTrackEta(0x0),
203  fhTrackPhi(0x0),
204  fhGroomedPtvJetPt(0x0),
205  fhDroppedBranches(0x0),
206  fhPtTriggerHadron(0x0),
207  fh2PtTriggerHadronJet(0x0),
208  fhPhiTriggerHadronJet(0x0),
209  fhPhiTriggerHadronEventPlane(0x0),
210  fhPhiTriggerHadronEventPlaneTPC(0x0),
211  fhDetJetPt_Incl(0x0),
212  fhDetJetPt_Matched(0x0),
213  fhJetPt_Det(0x0),
214  fhJetPt_True(0x0),
215  fhJetPhi_Det(0x0),
216  fhJetPhi_True(0x0),
217  fhJetEta_Det(0x0),
218  fhJetEta_True(0x0),
219  fhJetMass_Det(0x0),
220  fhJetMass_True(0x0),
221  fhJetRadius_Det(0x0),
222  fhJetRadius_True(0x0),
223  fhJetCounter_Det(0x0),
224  fhJetCounter_True(0x0),
225  fhNumberOfJetTracks_Det(0x0),
226  fhNumberOfJetTracks_True(0x0),
227  fh2PtRatio(0x0),
228  fhLundIterative(0x0),
229  fhLundIterativeTrue(0x0),
230  fReclusterAlgo(0),
231  fSubMatching(kFALSE),
232  bMinSubjetPt(kFALSE),
233  fMinSubjetPt(0)
234 {
235  // Standard constructor.
236  for(Int_t i=0;i<nBranch;i++){
237  fJetInfoVar[i]=0;
238  }
240  DefineOutput(1, TList::Class());
241  DefineOutput(2, TTree::Class());
242 }
243 
244 //________________________________________________________________________
246 {
247  // Destructor.
248 }
249 
250 //________________________________________________________________________
252 {
253  // Create user output.
254 
256 
257  Bool_t oldStatus = TH1::AddDirectoryStatus();
258  TH1::AddDirectory(kFALSE);
259  TH1::AddDirectory(oldStatus);
260  const char* nameoutput = GetOutputSlot(2)->GetContainer()->GetName();
261  fTreeJetInfo = new TTree(nameoutput, nameoutput);
262 
263  if (!fFullTree){
264  const Int_t nVarMin = 22;
265  TString *fJetInfoVarNames = new TString [nVarMin];
266 
267  fJetInfoVarNames[0] = "Pt";
268  fJetInfoVarNames[1] = "Pt_Truth";
269  fJetInfoVarNames[2] = "SymParam";
270  fJetInfoVarNames[3] = "SymParam_Truth";
271  fJetInfoVarNames[4] = "Mass";
272  fJetInfoVarNames[5] = "Mass_Truth";
273  fJetInfoVarNames[6] = "SLSubJetMass";
274  fJetInfoVarNames[7] = "SLSubJetMass_Truth";
275  fJetInfoVarNames[8] = "LeadingParton";
276  fJetInfoVarNames[9] = "LeadingParton_Truth";
277  fJetInfoVarNames[10] = "SLSubJetPt";
278  fJetInfoVarNames[11] = "SLSubJetPt_Truth";
279  fJetInfoVarNames[12] = "DelR";
280  fJetInfoVarNames[13] = "DelR_Truth";
281  fJetInfoVarNames[14] = "N_SplittingsHard";
282  fJetInfoVarNames[15] = "N_SplittingsHard_Truth";
283  fJetInfoVarNames[16] = "Groomed_Jet_Pt";
284  fJetInfoVarNames[17] = "Groomed_Jet_Pt_Truth";
285  fJetInfoVarNames[18] = "Groomed_Mass";
286  fJetInfoVarNames[19] = "Groomed_Mass_Truth";
287  fJetInfoVarNames[20] = "N_SplittingsAll";
288  fJetInfoVarNames[21] = "N_SplittingsAll_Truth";
289 
290 
291  for(Int_t ivar=0; ivar < nVarMin; ivar++){
292  cout<<"looping over variables"<<endl;
293  fTreeJetInfo->Branch(fJetInfoVarNames[ivar].Data(), &fJetInfoVar[ivar], Form("%s/D", fJetInfoVarNames[ivar].Data()));
294  }
295  }
296  if(fJetSelection == kRecoil){
297  fhPtTriggerHadron= new TH1F("fhPtTriggerHadron", "fhPtTriggerHadron",1500,-0.5,149.5);
299  fh2PtTriggerHadronJet= new TH2F("fh2PtTriggerHadronJet", "fh2PtTriggerHadronJet",1500,-0.5,149.5,1500,-0.5,149.5);
301  fhPhiTriggerHadronJet= new TH1F("fhPhiTriggerHadronJet", "fhPhiTriggerHadronJet",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
303  fhPhiTriggerHadronEventPlane= new TH1F("fhPhiTriggerHadronEventPlane", "fhPhiTriggerHadronEventPlane",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
305  fhPhiTriggerHadronEventPlaneTPC= new TH1F("fhPhiTriggerHadronEventPlaneTPC", "fhPhiTriggerHadronEventPlaneTPC",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
307 
308 
309  }
310 
311 
312 
314 
315  fhJetPt= new TH1F("fhJetPt", "Jet Pt",150,-0.5,149.5 );
316  fOutput->Add(fhJetPt);
317  fhJetPhi= new TH1F("fhJetPhi", "Jet Phi",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
318  fOutput->Add(fhJetPhi);
319  fhJetEta= new TH1F("fhJetEta", "Jet Eta", Eta_Bins_1, Eta_Lower, Eta_Upper);
320  fOutput->Add(fhJetEta);
321  fhJetMass= new TH1F("fhJetMass", "Jet Mass", 4000,-0.5, 39.5);
322  fOutput->Add(fhJetMass);
323  fhJetRadius= new TH1F("fhJetRadius", "Jet Radius", 100, -0.05,0.995);
324  fOutput->Add(fhJetRadius);
325  fhNumberOfJetTracks= new TH1F("fhNumberOfJetTracks", "Number of Tracks within a Jet", 300, -0.5,299.5);
327  fhJetCounter= new TH1F("fhJetCounter", "Jet Counter", 150, -0.5, 149.5);
328  fOutput->Add(fhJetCounter);
329  fhJetArea= new TH1F("fhJetArea", "Jet Area", 400,-0.5, 1.95);
330  fOutput->Add(fhJetArea);
331  fhTrackPt= new TH1F("fhTrackPt", "Track Pt",600,-0.5,59.5);
332  fOutput->Add(fhTrackPt);
333  fhTrackEta = new TH1F("fhTrackEta","Track Eta ",600,-1.5,1.5);
334  fOutput->Add(fhTrackEta);
335  fhTrackPhi = new TH1F("fhTrackPhi","Track Phi ",600,-1.5,1.5);
336  fOutput->Add(fhTrackPhi);
337  fhGroomedPtvJetPt= new TH2F("fhGroomedPtvJetPt","Groomed Jet p_{T} v Original Jet p_{T}",150,0,150,150,0,150);
339  fhDroppedBranches= new TH1F("fhDroppedBranches","Number of Softdropped branches",50,0,50);
341  fhDetJetPt_Incl= new TH1F("fhDetJetPt_Incl", "Jet Pt",200,-0.5,199.5 );
342  fOutput->Add(fhDetJetPt_Incl);
343  fhDetJetPt_Matched= new TH1F("fhDetJetPt_Matched", "Jet Pt",200,-0.5,199.5 );
345 
346  const Int_t dimSpec = 5;
347  const Int_t nBinsSpec[5] = {50,50,10,3,10};
348  const Double_t lowBinSpec[5] = {0.0,-10, 0,0,0};
349  const Double_t hiBinSpec[5] = {5.0, 0,200,3,10};
350  fhLundIterative = new THnSparseF("fHLundIterative",
351  "LundIterativePlot [log(1/theta),log(z*theta),pTjet,algo,ndepth]",
352  dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
353  fOutput->Add(fhLundIterative);
355 
356  fhLundIterativeTrue = new THnSparseF("fHLundIterativeTrue",
357  "LundIterativePlot [log(1/theta),log(z*theta),pTjet,algo,ndepth]",
358  dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
360  }
361 
362  }
364  fhJetPt_Det= new TH1F("fhJetPt_Det", "Jet Pt Detector Level",1500,-0.5,149.5 );
365  fOutput->Add(fhJetPt_Det);
366  fhJetPt_True= new TH1F("fhJetPt_True", "Jet Pt Particle Level",1500,-0.5,149.5 );
367  fOutput->Add(fhJetPt_True);
368  fhJetPhi_Det= new TH1F("fhJetPhi_Det", "Jet Phi Detector Level",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
369  fOutput->Add(fhJetPhi_Det);
370  fhJetPhi_True= new TH1F("fhJetPhi_True", "Jet Phi Particle Level",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
371  fOutput->Add(fhJetPhi_True);
372  fhJetEta_Det= new TH1F("fhJetEta_Det", "Jet Eta Detector Level", Eta_Bins_1, Eta_Lower, Eta_Upper);
373  fOutput->Add(fhJetEta_Det);
374  fhJetEta_True= new TH1F("fhJetEta_True", "Jet Eta Particle Level", Eta_Bins_1, Eta_Lower, Eta_Upper);
375  fOutput->Add(fhJetEta_True);
376  fhJetMass_Det= new TH1F("fhJetMass_Det", "Jet Mass Detector Level", 4000,-0.5, 39.5);
377  fOutput->Add(fhJetMass_Det);
378  fhJetMass_True= new TH1F("fhJetMass_True", "Jet Mass Particle Level", 4000,-0.5, 39.5);
379  fOutput->Add(fhJetMass_True);
380  fhJetRadius_Det= new TH1F("fhJetRadius_Det", "Jet Radius Detector Level", 100, -0.05,0.995);
381  fOutput->Add(fhJetRadius_Det);
382  fhJetRadius_True= new TH1F("fhJetRadius_True", "Jet Radius Particle Level", 100, -0.05,0.995);
384  fhNumberOfJetTracks_Det= new TH1F("fhNumberOfJetTracks_Det", "Number of Tracks within a Jet Detector Level", 300, -0.5,299.5);
386  fhNumberOfJetTracks_True= new TH1F("fhNumberOfJetTracks_True", "Number of Tracks within a Jet Particle Level", 300, -0.5,299.5);
388 
389  fhJetCounter_Det= new TH1F("fhJetCounter_Det", "Jet Counter Detector Level", 150, -0.5, 149.5);
391  fhJetCounter_True= new TH1F("fhJetCounter_True", "Jet Counter Particle Level", 150, -0.5, 149.5);
393  fh2PtRatio= new TH2F("fhPtRatio", "MC pt for detector level vs particle level jets",1500,-0.5,149.5,1500,-0.5,149.5);
394  fOutput->Add(fh2PtRatio);
395 
396  const Int_t dimSpec = 5;
397  const Int_t nBinsSpec[5] = {50,50,10,3,10};
398  const Double_t lowBinSpec[5] = {0.0,-10, 0,0,0};
399  const Double_t hiBinSpec[5] = {5.0, 0,200,3,10};
400  fhLundIterative = new THnSparseF("fHLundIterative",
401  "LundIterativePlot [log(1/theta),log(z*theta),pTjet,algo,ndepth]",
402  dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
403  fOutput->Add(fhLundIterative);
404 
405  fhLundIterativeTrue = new THnSparseF("fHLundIterativeTrue",
406  "LundIterativePlot [log(1/theta),log(z*theta),pTjet,algo,ndepth]",
407  dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
409 
410  }
411  PostData(1,fOutput);
412  PostData(2,fTreeJetInfo);
413  //cout<<"End of UserCreateOutputObjects"<<endl;
414  // delete [] fShapesVarNames;
415 }
416 
417 //________________________________________________________________________
419 {
420  // Run analysis code here, if needed. It will be executed before FillHistograms().
421  return kTRUE;
422 }
423 
424 //________________________________________________________________________
426 {
427  //AliMultSelection *MultSelection = 0x0;
428  AliAODEvent *fEvent = dynamic_cast<AliAODEvent*>(InputEvent());
429  if (!fEvent) {
430  Error("UserExec","AOD not available");
431  return 0;
432  }
433  /*MultSelection = (AliMultSelection * ) fEvent->FindListObject("MultSelection");
434  if(!MultSelection) {
435  AliWarning("AliMultSelection object not found!");
436  }
437  else{
438  Percentile_1 = MultSelection->GetMultiplicityPercentile("V0M");
439  }
440  */
441  if (fCentSelectOn){
442  if ((fCent>fCentMax) || (fCent<fCentMin)) return 0;
443  }
444 
446  AliAODTrack *TriggerHadron = 0x0;
447  if (fJetSelection == kRecoil) {
448  //you have to set a flag and the limits of the pT interval for your trigger
450  if (TriggerHadronLabel==-99999) return 0; //Trigger Hadron Not Found
451  AliTrackContainer *PartCont =NULL;
452  AliParticleContainer *PartContMC = NULL;
453  if (fJetShapeSub==kConstSub){
455  else PartCont = GetTrackContainer(1);
456  }
457  else {
459  else PartCont = GetTrackContainer(0);
460  }
461  TClonesArray *TrackArray = NULL;
462  TClonesArray *TrackArrayMC = NULL;
463  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kGenOnTheFly) TrackArrayMC = PartContMC->GetArray();
464  else TrackArray = PartCont->GetArray();
465  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kGenOnTheFly) TriggerHadron = static_cast<AliAODTrack*>(TrackArrayMC->At(TriggerHadronLabel));
466  else TriggerHadron = static_cast<AliAODTrack*>(TrackArray->At(TriggerHadronLabel));
467  if (!TriggerHadron) return 0;//No trigger hadron with label found
468  if(fSemigoodCorrect){
469  Double_t HoleDistance=RelativePhi(TriggerHadron->Phi(),fHolePos);
470  if(TMath::Abs(HoleDistance)+fHoleWidth+fJetRadius>TMath::Pi()-fRecoilAngularWindow) return 0;
471  }
472  fhPtTriggerHadron->Fill(TriggerHadron->Pt()); //Needed for per trigger Normalisation
473  if (fJetShapeType != AliAnalysisTaskRecoilJetYield::kGenOnTheFly) fhPhiTriggerHadronEventPlane->Fill(TMath::Abs(RelativePhiEventPlane(fEPV0,TriggerHadron->Phi()))); //fEPV0 is the event plane from AliAnalysisTaskEmcal
474  if (fJetShapeType != AliAnalysisTaskRecoilJetYield::kGenOnTheFly) fhPhiTriggerHadronEventPlaneTPC->Fill(TMath::Abs(RelativePhiEventPlane(((AliVAODHeader*)(InputEvent()->GetHeader()))->GetEventplane(),TriggerHadron->Phi()))); //TPC event plane
475  }
476 
477 
478 
479 
480 
483  // cout<<"entering kData"<<endl;
484 
485  AliEmcalJet *Jet1 = NULL; //Original Jet in the event
486  AliJetContainer *JetCont= GetJetContainer(0); //Jet Container for event
487  Int_t JetCounter=0; //Counts number of jets in event
488  Double_t JetPhi=0,TrackPhi=0;
489  Bool_t EventCounter=kFALSE;
490  Double_t JetPt_ForThreshold=0;
491  AliEmcalJet *Jet2= NULL;
492  if(JetCont) {
493  JetCont->ResetCurrentID();
494  while((Jet1=JetCont->GetNextAcceptJet())) {
495  if(!Jet1) continue;
496  if (fJetShapeSub==kNoSub) JetPt_ForThreshold = Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
497  else JetPt_ForThreshold = Jet1->Pt();
498  if(JetPt_ForThreshold<fPtThreshold) {
499  continue;
500  }
501  else {
502  Float_t RecoilDeltaPhi = 0.;
503  if (fJetSelection == kRecoil){
504  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), Jet1->Phi());
505  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
506  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), Jet1->Pt());
507  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), Jet1->Phi()));
508  }
509  if (!EventCounter){
510  EventCounter=kTRUE;
511  }
512  JetCounter++;
513  if(fJetShapeSub==kNoSub) fhJetPt->Fill(Jet1->Pt()-(GetRhoVal(0)*Jet1->Area()));
514  else fhJetPt->Fill(Jet1->Pt());
515  AliTrackContainer *TrackCont = NULL;
516  if(fJetShapeSub==kConstSub) TrackCont = GetTrackContainer(1);
517  else TrackCont = GetTrackContainer(0);
518  TClonesArray *TracksArray = NULL;
519  TracksArray = TrackCont->GetArray();
520  Double_t NTracks = TracksArray->GetEntriesFast();
521  AliAODTrack *Track = 0x0;
522  for(Int_t i=0; i < NTracks; i++){
523  Track = static_cast<AliAODTrack*>(TrackCont->GetAcceptTrack(i));
524  if (!Track) continue;
525  TrackPhi=Track->Phi();
526  if(TrackPhi < -1*TMath::Pi()) TrackPhi += (2*TMath::Pi());
527  else if (TrackPhi > TMath::Pi()) TrackPhi -= (2*TMath::Pi());
528  fhTrackPt->Fill(Track->Pt());
529  fhTrackEta->Fill(Track->Eta());
530  fhTrackPhi->Fill(TrackPhi);
531  }
532  fhJetArea->Fill(Jet1->Area());
533  JetPhi=Jet1->Phi();
534  if(JetPhi < -1*TMath::Pi()) JetPhi += (2*TMath::Pi());
535  else if (JetPhi > TMath::Pi()) JetPhi -= (2*TMath::Pi());
536  fhJetPhi->Fill(JetPhi);
537  fhJetEta->Fill(Jet1->Eta());
538  fhJetMass->Fill(Jet1->M());
539  fJetInfoVar[4]=Jet1->M();
540  fhJetRadius->Fill(TMath::Sqrt((Jet1->Area()/TMath::Pi()))); //Radius of Jets per event
541  fhNumberOfJetTracks->Fill(Jet1->GetNumberOfTracks());
542  if(fJetShapeSub==kNoSub) fJetInfoVar[0]= Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
543  else fJetInfoVar[0]=Jet1->Pt();
544 
545  fJetInfoVar[1]=0;
546  if(fDoSoftDrop) {
547  SoftDrop(Jet1,JetCont,fZCut,fBeta_SD,kFALSE);
548  RecursiveParents(Jet1,JetCont,0,kFALSE);
549  RecursiveParents(Jet1,JetCont,1,kFALSE);
550  RecursiveParents(Jet1,JetCont,2,kFALSE);
551  //SoftDrop(Jet1,JetCont,fZCut,fBeta_SD,kTRUE);
552 
553  fJetInfoVar[3]=0;
554  fJetInfoVar[5]=0;
555  fJetInfoVar[7]=0;
556  fJetInfoVar[9]=0;
557  fJetInfoVar[11]=0;
558  fJetInfoVar[13]=0;
559  fJetInfoVar[15]=0;
560  fJetInfoVar[17]=0;
561  fJetInfoVar[19]=0;
562  fJetInfoVar[21]=0;
563  }
564  else{
565  fJetInfoVar[2]=0;
566  fJetInfoVar[3]=0;
567  fJetInfoVar[4]=0;
568  fJetInfoVar[5]=0;
569  fJetInfoVar[6]=0;
570  fJetInfoVar[7]=0;
571  fJetInfoVar[8]=0;
572  fJetInfoVar[9]=0;
573  fJetInfoVar[10]=0;
574  fJetInfoVar[11]=0;
575  fJetInfoVar[12]=0;
576  fJetInfoVar[13]=0;
577  fJetInfoVar[14]=0;
578  fJetInfoVar[15]=0;
579  fJetInfoVar[16]=0;
580  fJetInfoVar[17]=0;
581  fJetInfoVar[18]=0;
582  fJetInfoVar[19]=0;
583  fJetInfoVar[20]=0;
584  fJetInfoVar[21]=0;
585 
586  }
587  fTreeJetInfo->Fill();
588  }
589  }
590  fhJetCounter->Fill(JetCounter); //Number of Jets in Each Event
591  }
592  }
593 
595  if(fJetShapeType == kDetEmbPart){
596  AliEmcalJet *JetHybridS = NULL; //Subtracted hybrid Jet
597  AliEmcalJet *JetHybridUS = NULL; //Unsubtracted Hybrid Jet //For matching SubtractedHybrid->DetPythia this jet container is also Subtracted Hybrid
598  AliEmcalJet *JetPythDet = NULL; //Detector Level Pythia Jet
599  AliEmcalJet *JetPythTrue = NULL; //Particle Level Pyhtia Jet
600  AliJetContainer *JetContHybridS= GetJetContainer(0); //Jet Container for Subtracted Hybrid Jets
601  AliJetContainer *JetContHybridUS= GetJetContainer(1); //Jet Container for Unsubtracted Hybrid Jets
602  AliJetContainer *JetContPythDet= GetJetContainer(2); //Jet Container for Detector Level Pyhtia Jets
603  AliJetContainer *JetContPythTrue= GetJetContainer(3); //Jet Container for Particle Level Pythia Jets
604 
605 
606 
607  Bool_t JetsMatched = kFALSE;
608  Double_t JetPtThreshold;
609  JetContHybridS->ResetCurrentID();
610  JetContHybridUS->ResetCurrentID();
611  JetContPythDet->ResetCurrentID();
612  JetContPythTrue->ResetCurrentID();
613 
614  while((JetPythDet = JetContPythDet->GetNextAcceptJet())){
615  fhDetJetPt_Incl->Fill(JetPythDet->Pt()); //Fill histogram with all Detector level jets
616  }
617 
618 
619  while((JetHybridS = JetContHybridS->GetNextAcceptJet())){
620  if (fJetShapeSub==kConstSub) JetPtThreshold=JetHybridS->Pt();
621  else JetPtThreshold=JetHybridS->Pt()-(GetRhoVal(0)*JetHybridS->Area());
622  if ( (!JetHybridS) || (JetPtThreshold<fPtThreshold)) continue;
623  Int_t JetNumber=-1;
624  for(Int_t i = 0; i<JetContHybridUS->GetNJets(); i++) {
625  JetHybridUS = JetContHybridUS->GetJet(i);
626  if (!JetHybridUS) continue;
627 
628  if(JetHybridUS->GetLabel()==JetHybridS->GetLabel()) {
629  JetNumber=i;
630  }
631  }
632  if(JetNumber==-1) continue;
633  JetHybridUS=JetContHybridUS->GetJet(JetNumber);
634  if (JetContHybridUS->AliJetContainer::GetFractionSharedPt(JetHybridUS)<fSharedFractionPtMin && !fSubMatching) {
635  continue;
636  }
638  continue;
639  }
640  JetPythDet=JetHybridUS->ClosestJet();
641  if (!JetHybridUS) {
642  Printf("Unsubtracted embedded jet does not exist, returning");
643  continue;
644  }
645  if (!JetPythDet) continue;
646  UInt_t rejectionReason = 0;
647  if (!(JetContPythDet->AcceptJet(JetPythDet,rejectionReason))) continue;
648  fhDetJetPt_Matched->Fill(JetPythDet->Pt()); //Fill only matched detector level jets for tagging efficiency comparison
649  JetPythTrue=JetPythDet->ClosestJet();
650  if(!JetPythTrue) continue;
651  JetsMatched=kTRUE;
652 
653  if(fJetShapeSub==kConstSub) fJetInfoVar[0]=JetHybridS->Pt();
654  else fJetInfoVar[0]=JetHybridS->Pt()-(GetRhoVal(0)*JetHybridS->Area());
655  fJetInfoVar[1]=JetPythTrue->Pt();
656  if(fDoSoftDrop) {
657  SoftDrop(JetHybridS,JetContHybridS,fZCut,fBeta_SD,kFALSE);
658  RecursiveParents(JetHybridS,JetContHybridS,0,kFALSE);
659  RecursiveParents(JetHybridS,JetContHybridS,1,kFALSE);
660  RecursiveParents(JetHybridS,JetContHybridS,2,kFALSE);
661  SoftDrop(JetPythTrue,JetContPythTrue,fZCut,fBeta_SD,kTRUE);
662  RecursiveParents(JetPythTrue,JetContPythTrue,0,kTRUE);
663  RecursiveParents(JetPythTrue,JetContPythTrue,1,kTRUE);
664  RecursiveParents(JetPythTrue,JetContPythTrue,2,kTRUE);
665  }
666  else{
667  fJetInfoVar[2]=0;
668  fJetInfoVar[3]=0;
669  fJetInfoVar[4]=0;
670  fJetInfoVar[5]=0;
671  fJetInfoVar[6]=0;
672  fJetInfoVar[7]=0;
673  fJetInfoVar[8]=0;
674  fJetInfoVar[9]=0;
675  fJetInfoVar[10]=0;
676  fJetInfoVar[11]=0;
677  fJetInfoVar[12]=0;
678  fJetInfoVar[13]=0;
679  fJetInfoVar[14]=0;
680  fJetInfoVar[15]=0;
681  fJetInfoVar[16]=0;
682  fJetInfoVar[17]=0;
683  fJetInfoVar[18]=0;
684  fJetInfoVar[19]=0;
685  fJetInfoVar[20]=0;
686  fJetInfoVar[21]=0;
687 
688  }
689  fJetInfoVar[4]=JetHybridS->M();
690  fJetInfoVar[5]=JetPythTrue->M();
691  fTreeJetInfo->Fill();
692  }
693 
694  }
695 
697  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kTrueDet){ //Truth->Detector response
698  AliEmcalJet *JetDet = NULL; //Detector Level Jet
699  AliEmcalJet *JetTrue = NULL; //Particle Level Jet
700  AliJetContainer *JetContDet= GetJetContainer(0); //Jet Container for Detector Level Pythia
701  AliJetContainer *JetContTrue= GetJetContainer(1); //Jet Container for Particle Level Pythia
702  Double_t JetPhiDet=0;
703  Double_t JetPhiTrue=0;
704  Bool_t JetsMatched=kFALSE;
705  Double_t Pythia_Event_Weight=1;
706  Bool_t EventCounter=kFALSE;
707  Int_t JetCounter_Det=0,JetCounter_True=0;
708  JetContDet->ResetCurrentID();
709  JetContTrue->ResetCurrentID();
710  if(JetContDet) {
711  while((JetDet=JetContDet->GetNextAcceptJet())) {
712  if( (!JetDet) || ((JetDet->Pt())<fPtThreshold)) {
713  // fhEventCounter_1->Fill(3); //events where the jet had a problem
714  continue;
715  }
716  else {
717  /* if(fSemigoodCorrect){
718  Double_t HoleDistance=RelativePhi(Jet1->Phi(),fHolePos);
719  if(TMath::Abs(HoleDistance)<fHoleWidth) continue;
720  }*/
721  Float_t RecoilDeltaPhi = 0.;
722  if (fJetSelection == kRecoil){
723  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), JetDet->Phi());
724  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
725  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), JetDet->Pt());
726  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), JetDet->Phi()));
727  }
728  JetCounter_Det++;
729  fhJetPt_Det->Fill(JetDet->Pt());
730  JetPhiDet=JetDet->Phi();
731  if(JetPhiDet < -1*TMath::Pi()) JetPhiDet += (2*TMath::Pi());
732  else if (JetPhiDet > TMath::Pi()) JetPhiDet -= (2*TMath::Pi());
733  fhJetPhi_Det->Fill(JetPhiDet);
734  fhJetEta_Det->Fill(JetDet->Eta());
735  fhJetMass_Det->Fill(JetDet->M());
736  fhJetRadius_Det->Fill(TMath::Sqrt((JetDet->Area()/TMath::Pi())));
738  if((JetTrue = JetDet->ClosestJet())){
739  JetsMatched=kTRUE;
740  JetCounter_True++;
741  fhJetPt_True->Fill(JetTrue->Pt());
742  JetPhiTrue=JetTrue->Phi();
743  if(JetPhiTrue < -1*TMath::Pi()) JetPhiTrue += (2*TMath::Pi());
744  else if (JetPhiTrue > TMath::Pi()) JetPhiTrue -= (2*TMath::Pi());
745  fhJetPhi_True->Fill(JetPhiTrue);
746  fhJetEta_True->Fill(JetTrue->Eta());
747  fhJetMass_True->Fill(JetTrue->M());
748  fhJetRadius_True->Fill(TMath::Sqrt((JetTrue->Area()/TMath::Pi())));
750  fh2PtRatio->Fill(JetDet->Pt(),JetTrue->Pt(),Pythia_Event_Weight);
751  }
752  else continue;
753  if(fJetShapeSub==kNoSub) fJetInfoVar[0]= JetDet->Pt()-(GetRhoVal(0)*JetDet->Area());
754  else fJetInfoVar[0]=JetDet->Pt();
755  if(fJetShapeSub==kNoSub) fJetInfoVar[1]= JetTrue->Pt()-(GetRhoVal(0)*JetTrue->Area());
756  else fJetInfoVar[1]=JetTrue->Pt();
757  if(fDoSoftDrop) {
758  SoftDrop(JetDet,JetContDet,fZCut,fBeta_SD,kFALSE);
759  RecursiveParents(JetDet,JetContDet,0,kFALSE);
760  RecursiveParents(JetDet,JetContDet,1,kFALSE);
761  RecursiveParents(JetDet,JetContDet,2,kFALSE);
762  SoftDrop(JetTrue,JetContTrue,fZCut,fBeta_SD,kTRUE);
763  RecursiveParents(JetTrue,JetContTrue,0,kTRUE);
764  RecursiveParents(JetTrue,JetContTrue,1,kTRUE);
765  RecursiveParents(JetTrue,JetContTrue,2,kTRUE);
766  }
767  else{
768  fJetInfoVar[2]=0;
769  fJetInfoVar[3]=0;
770  fJetInfoVar[4]=0;
771  fJetInfoVar[5]=0;
772  fJetInfoVar[6]=0;
773  fJetInfoVar[7]=0;
774  fJetInfoVar[8]=0;
775  fJetInfoVar[9]=0;
776  fJetInfoVar[10]=0;
777  fJetInfoVar[11]=0;
778  fJetInfoVar[12]=0;
779  fJetInfoVar[13]=0;
780  fJetInfoVar[14]=0;
781  fJetInfoVar[15]=0;
782  fJetInfoVar[16]=0;
783  fJetInfoVar[17]=0;
784  fJetInfoVar[18]=0;
785  fJetInfoVar[19]=0;
786  fJetInfoVar[20]=0;
787  fJetInfoVar[21]=0;
788 
789  }
790  fJetInfoVar[4]=JetDet->M();
791  fJetInfoVar[5]=JetTrue->M();
792  fTreeJetInfo->Fill();
793 
794  JetsMatched=kFALSE;
795  }
796  }
797  fhJetCounter_Det->Fill(JetCounter_Det); //Number of Jets in Each Event Particle Level
798  fhJetCounter_True->Fill(JetCounter_True); //Number of Jets in Each Event Detector Level
799  }
800  }
801 
802 
805  AliEmcalJet *Jet1 = NULL; //Original Jet in the event
806  AliJetContainer *JetCont= GetJetContainer(0); //Jet Container for event
807  Int_t JetCounter=0; //Counts number of jets in event
808  Double_t JetPhi=0;
809  Bool_t EventCounter=kFALSE;
810  Double_t JetPt_ForThreshold=0;
811  AliEmcalJet *Jet2= NULL;
812  if(JetCont) {
813  JetCont->ResetCurrentID();
814  while((Jet1=JetCont->GetNextAcceptJet())) {
815  if(!Jet1) continue;
816  if (fJetShapeSub==kNoSub) JetPt_ForThreshold = Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
817  else JetPt_ForThreshold = Jet1->Pt();
818  if(JetPt_ForThreshold<fPtThreshold) {
819  continue;
820  }
821  else {
822  Float_t RecoilDeltaPhi = 0.;
823  if (fJetSelection == kRecoil){
824  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), Jet1->Phi());
825  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
826  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), Jet1->Pt());
827  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), Jet1->Phi()));
828  }
829  if (!EventCounter){
830  EventCounter=kTRUE;
831  }
832  JetCounter++;
833  fhJetPt->Fill(Jet1->Pt());
834  fhJetArea->Fill(Jet1->Area());
835  JetPhi=Jet1->Phi();
836  if(JetPhi < -1*TMath::Pi()) JetPhi += (2*TMath::Pi());
837  else if (JetPhi > TMath::Pi()) JetPhi -= (2*TMath::Pi());
838  fhJetPhi->Fill(JetPhi);
839  fhJetEta->Fill(Jet1->Eta());
840  fhJetMass->Fill(Jet1->M());
841  fhJetRadius->Fill(TMath::Sqrt((Jet1->Area()/TMath::Pi()))); //Radius of Jets per event
842  fhNumberOfJetTracks->Fill(Jet1->GetNumberOfTracks());
843  if(fJetShapeSub==kNoSub) fJetInfoVar[0]= Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
844  else fJetInfoVar[0]=Jet1->Pt();
845  fJetInfoVar[1]=0;
846  if(fDoSoftDrop) {
847  SoftDrop(Jet1,JetCont,fZCut,fBeta_SD,kFALSE);
848  SoftDrop(Jet1,JetCont,fZCut,fBeta_SD,kTRUE);
849  }
850  else{
851  fJetInfoVar[2]=0;
852  fJetInfoVar[3]=0;
853  fJetInfoVar[4]=0;
854  fJetInfoVar[5]=0;
855  fJetInfoVar[6]=0;
856  fJetInfoVar[7]=0;
857  fJetInfoVar[8]=0;
858  fJetInfoVar[9]=0;
859  fJetInfoVar[10]=0;
860  fJetInfoVar[11]=0;
861  fJetInfoVar[12]=0;
862  fJetInfoVar[13]=0;
863  fJetInfoVar[14]=0;
864  fJetInfoVar[15]=0;
865  fJetInfoVar[16]=0;
866  fJetInfoVar[17]=0;
867  fJetInfoVar[18]=0;
868  fJetInfoVar[19]=0;
869  fJetInfoVar[20]=0;
870  fJetInfoVar[21]=0;
871  }
872  fJetInfoVar[4]=Jet1->M();
873  fJetInfoVar[5]=0;
874  fTreeJetInfo->Fill();
875  }
876  }
877  fhJetCounter->Fill(JetCounter); //Number of Jets in Each Event
878  }
879  }
880 
881 
882 
883  return kTRUE;
884 }
885 
886 //________________________________________________________________________
888 
889  if(Phi < -1*TMath::Pi()) Phi += (2*TMath::Pi());
890  else if (Phi > TMath::Pi()) Phi -= (2*TMath::Pi());
891  Double_t DeltaPhi=Phi-EventPlane;
892  if(DeltaPhi < -1*TMath::Pi()) DeltaPhi += (2*TMath::Pi());
893  else if (DeltaPhi > TMath::Pi()) DeltaPhi -= (2*TMath::Pi());
894  return DeltaPhi;
895 }
896 //________________________________________________________________________
898 
899  if(Phi1 < -1*TMath::Pi()) Phi1 += (2*TMath::Pi());
900  else if (Phi1 > TMath::Pi()) Phi1 -= (2*TMath::Pi());
901  if(Phi2 < -1*TMath::Pi()) Phi2 += (2*TMath::Pi());
902  else if (Phi2 > TMath::Pi()) Phi2 -= (2*TMath::Pi());
903  Double_t DeltaPhi=Phi2-Phi1;
904  if(DeltaPhi < -1*TMath::Pi()) DeltaPhi += (2*TMath::Pi());
905  else if (DeltaPhi > TMath::Pi()) DeltaPhi -= (2*TMath::Pi());
906  return DeltaPhi;
907 }
908 
909 
910 
911 
912 //--------------------------------------------------------------------------
914 
915  AliJetContainer *JetCont = GetJetContainer(JetContNb);
916  Double_t Angularity_Numerator=0;
917  Double_t Angularity_Denominator=0;
918  AliVParticle *Particle=0x0;
919  Double_t DeltaPhi=0;
920 
921 
922  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks
923  Particle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
924  if(!Particle) continue;
925  DeltaPhi=RelativePhi(Jet->Phi(),Particle->Phi());
926  Angularity_Numerator=Angularity_Numerator+(Particle->Pt()*TMath::Sqrt(((Particle->Eta()-Jet->Eta())*(Particle->Eta()-Jet->Eta()))+(DeltaPhi*DeltaPhi)));
927  Angularity_Denominator= Angularity_Denominator+Particle->Pt();
928  }
929  if(Angularity_Denominator!=0) return Angularity_Numerator/Angularity_Denominator;
930  else return -1;
931 
932 }
933 
934 
935 
936 //--------------------------------------------------------------------------
938 
939  AliJetContainer *JetCont = GetJetContainer(JetContNb);
940  Double_t PTD_Numerator=0; //Reset these values
941  Double_t PTD_Denominator=0;
942  AliVParticle *Particle=0x0;
943  Double_t DeltaPhi=0;
944  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks
945  Particle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
946  if(!Particle) continue;
947  PTD_Numerator=PTD_Numerator+(Particle->Pt()*Particle->Pt());
948  PTD_Denominator=PTD_Denominator+Particle->Pt();
949  }
950  if(PTD_Denominator!=0) return TMath::Sqrt(PTD_Numerator)/PTD_Denominator;
951  else return -1;
952 
953 }
954 
955 //_________________________________________________________________________
956 void AliAnalysisTaskRecoilJetYield::SoftDrop(AliEmcalJet *fJet,AliJetContainer *fJetCont, double zcut, double beta, Bool_t fTruthJet){
957  try {
958  std::vector<fastjet::PseudoJet> fInputVectors;
959  Double_t JetInvMass=0, PseudJetInvMass=0, TrackMom = 0, TrackEnergy = 0;
960  AliParticleContainer *fTrackCont = fJetCont->GetParticleContainer();
961  Double_t JetEta=fJet->Eta(),JetPhi=fJet->Phi();
962  Double_t FJTrackEta[9999],FJTrackPhi[9999],FJTrackPt[9999],EmcalJetTrackEta[9999],EmcalJetTrackPhi[9999],EmcalJetTrackPt[9999];
963  UShort_t FJNTracks=0,EmcalJetNTracks=0;
964  if (fTrackCont) for (Int_t i=0; i<fJet->GetNumberOfTracks(); i++) {
965  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
966  JetInvMass += fTrk->M();
967  if (!fTrk) continue;
968  fastjet::PseudoJet PseudoTracks(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
969  TrackMom += TMath::Sqrt(TMath::Power(fTrk->Px(),2)+TMath::Power(fTrk->Py(),2)+TMath::Power(fTrk->Pz(),2));
970  TrackEnergy += fTrk->E();
971  PseudoTracks.set_user_index(fJet->TrackAt(i)+100);
972  PseudJetInvMass += PseudoTracks.m();
973  fInputVectors.push_back(PseudoTracks);
974  EmcalJetTrackEta[i]=fTrk->Eta();
975  EmcalJetTrackPhi[i]=fTrk->Phi();
976  EmcalJetTrackPt[i]=fTrk->Pt();
977  EmcalJetNTracks++;
978  }
979 
980 
981  fastjet::JetDefinition fJetDef = fastjet::JetDefinition(fastjet::antikt_algorithm, fJetRadius*2, static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
982  fastjet::ClusterSequence fClustSeqSA = fastjet::ClusterSequence(fInputVectors,fJetDef);
983 
984 
985  std::vector<fastjet::PseudoJet> fOutputJets;
986  fOutputJets.clear();
987  fOutputJets=fClustSeqSA.inclusive_jets(0);
988 
989  std::vector<fastjet::PseudoJet> jet_constituents = fOutputJets[0].constituents();
990  Float_t NSubjettinessResult[3], NSubBeta = 1, R0=0.4;
991  std::vector<fastjet::PseudoJet> Subjet_Axes;
992  fastjet::PseudoJet SubJet1_Axis,SubJet2_Axis;
993  Double_t DelR=-5;
994 
995 
996  for(Int_t j=1; j<3; j++){
997  if(jet_constituents.size() < j){
998  if(!fTruthJet){
999  //if(j==1) fJetInfoVar[4]=-5;
1000  // if(j==2) fJetInfoVar[6]=-5;
1001  }
1002  else {
1003  //if(j==1) fJetInfoVar[5]=-5;
1004  // if(j==2) fJetInfoVar[7]=-5;
1005  }
1006  continue;
1007  }
1008  else{
1009  fastjet::contrib::Nsubjettiness nSub(j,fastjet::contrib::KT_Axes(),fastjet::contrib::NormalizedMeasure(NSubBeta,R0));
1010  NSubjettinessResult[j] = nSub.result(fOutputJets[0]);
1011  if(j==2){ //find subjet axis to calculate delR
1012  Subjet_Axes = nSub.currentAxes();
1013  SubJet1_Axis = Subjet_Axes[0];
1014  SubJet2_Axis = Subjet_Axes[1];
1015 
1016  Double_t SubJet1_Eta=SubJet1_Axis.pseudorapidity();
1017  Double_t SubJet2_Eta=SubJet2_Axis.pseudorapidity();
1018  Double_t SubJet1_Phi=SubJet1_Axis.phi();
1019  if(SubJet1_Phi < -1*TMath::Pi()) SubJet1_Phi += (2*TMath::Pi());
1020  else if (SubJet1_Phi > TMath::Pi()) SubJet1_Phi -= (2*TMath::Pi());
1021  Double_t SubJet2_Phi=SubJet2_Axis.phi();
1022  if(SubJet1_Phi < -1*TMath::Pi()) SubJet1_Phi += (2*TMath::Pi());
1023  else if (SubJet1_Phi > TMath::Pi()) SubJet1_Phi -= (2*TMath::Pi());
1024 
1025  Double_t DeltaPhiSubJets,DeltaEtaSubJets;
1026  DeltaPhiSubJets=SubJet1_Phi-SubJet2_Phi;
1027  DeltaEtaSubJets=SubJet1_Eta-SubJet2_Eta;
1028  if(DeltaPhiSubJets < -1*TMath::Pi()) DeltaPhiSubJets += (2*TMath::Pi());
1029  else if (DeltaPhiSubJets > TMath::Pi()) DeltaPhiSubJets -= (2*TMath::Pi());
1030 
1031  DelR = TMath::Power(TMath::Power(DeltaPhiSubJets,2)+TMath::Power(DeltaEtaSubJets,2),0.5);
1032  }
1033  }
1034  }
1035  if(!fTruthJet){
1036  //fJetInfoVar[4]=NSubjettinessResult[1];
1037  //fJetInfoVar[6]=NSubjettinessResult[2];
1038  // fJetInfoVar[12]=DelR;
1039  }
1040  else {
1041  //fJetInfoVar[5]=NSubjettinessResult[1];
1042  //fJetInfoVar[7]=NSubjettinessResult[2];
1043  // fJetInfoVar[13]=DelR;
1044  }
1045 
1046 
1047 
1048  fastjet::contrib::SoftDrop softdrop(beta, zcut);
1049  //fastjet::contrib::SoftDrop softdrop_antikt(beta,zcut);
1050  softdrop.set_verbose_structure(kTRUE);
1051  //fastjet::JetDefinition jet_def_akt(fastjet::antikt_algorithm, 0.4);
1052  // fastjet::contrib::Recluster *antiKT_Recluster(jet_def_akt);
1053  fastjet::JetAlgorithm algo;
1054  switch(fReclusterAlgo){
1055  case 0: algo = fastjet::cambridge_aachen_algorithm; break;
1056  case 1: algo = fastjet::antikt_algorithm; break;
1057  case 2: algo = fastjet::kt_algorithm; break;
1058  };
1059 #if FASTJET_VERSION_NUMBER >= 30302
1060  fastjet::Recluster recluster(algo, 1, fastjet::Recluster::keep_only_hardest);
1061 #else
1062  fastjet::contrib::Recluster recluster(algo, 1, true);
1063 #endif
1064  softdrop.set_reclustering(true,&recluster);
1065  fastjet::PseudoJet finaljet = softdrop(fOutputJets[0]);
1066  // fastjet::PseudoJet finaljet_antikt = softdrop_antikt(fOutputJets[0]);
1067  //cout<< finaljet.structure_of<fastjet::contrib::SoftDrop>().symmetry()<<endl;
1068  //cout<< finaljet_antikt.structure_of<fastjet::contrib::SoftDrop>().symmetry()<<endl;
1069 
1070 
1071  AliEmcalJet* jet = new AliEmcalJet(finaljet.perp(), finaljet.eta(), finaljet.phi(), finaljet.m());
1072  std::vector<fastjet::PseudoJet> fSDTracks=finaljet.constituents();
1073  Double_t FastjetTrackDelR,EmcalTrackDelR;
1074  for(Int_t i=0;i<fJet->GetNumberOfConstituents();i++){
1075  if(i<=finaljet.constituents().size()){
1076  FastjetTrackDelR = TMath::Sqrt(TMath::Power(fSDTracks[i].eta()-JetEta,2)+TMath::Power(fSDTracks[i].phi()-JetPhi,2));
1077  FJTrackEta[i]=fSDTracks[i].eta();
1078  FJTrackPhi[i]=fSDTracks[i].phi();
1079  FJTrackPt[i]=fSDTracks[i].perp();
1080  FJNTracks++;
1081  }
1082  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
1083  EmcalTrackDelR = TMath::Sqrt(TMath::Power(fTrk->Eta()-JetEta,2)+TMath::Power(fTrk->Phi()-JetPhi,2));
1084  }
1085  Int_t NDroppedTracks = fJet->GetNumberOfTracks()-finaljet.constituents().size();
1086  Int_t nConstituents(fClustSeqSA.constituents(finaljet).size());
1087  jet->SetNumberOfTracks(nConstituents);
1088  Double_t SymParam, Mu, DeltaR;
1089  SymParam=(finaljet.structure_of<fastjet::contrib::SoftDrop>().symmetry());
1090  Mu=(finaljet.structure_of<fastjet::contrib::SoftDrop>().mu());
1091  DeltaR=(finaljet.structure_of<fastjet::contrib::SoftDrop>().delta_R());
1092  //fhGroomedPtvJetPt->Fill(finaljet.perp(),fJet->Pt());
1093  //fhDroppedBranches->Fill(finaljet.structure_of<fastjet::contrib::SoftDrop>().dropped_count());
1094  std::vector<fastjet::PseudoJet> subjets;
1095  if ( finaljet.has_pieces() ) {
1096  subjets = finaljet.pieces();
1097  fastjet::PseudoJet subjet1 = subjets[0];
1098  fastjet::PseudoJet subjet2 = subjets[1];
1099  if(!fTruthJet){
1100  if(subjets[0].perp() > subjets[1].perp()){
1101  fJetInfoVar[6]=subjets[1].m();
1102  fJetInfoVar[8]=LeadingTrackPt(subjets[1]);
1103  fJetInfoVar[10]=subjets[1].perp();
1104  }
1105  else {
1106  fJetInfoVar[6]= subjets[0].m();
1107  fJetInfoVar[8]=LeadingTrackPt(subjets[0]);
1108  fJetInfoVar[10]= subjets[0].perp();
1109  }
1110  }
1111  else {
1112  if(subjets[0].perp() > subjets[1].perp()){
1113  fJetInfoVar[7]=subjets[1].m();
1114  fJetInfoVar[9]=LeadingTrackPt(subjets[1]);
1115  fJetInfoVar[11]=subjets[1].perp();
1116  }
1117  else{
1118  fJetInfoVar[7]= subjets[0].m();
1119  fJetInfoVar[9]=LeadingTrackPt(subjets[0]);
1120  fJetInfoVar[11]=subjets[0].perp();
1121  }
1122  }
1123 
1124  }
1125  else {
1126  if(!fTruthJet){
1127  fJetInfoVar[6]=0;
1128  fJetInfoVar[10]=-1;
1129  }
1130  else {
1131  fJetInfoVar[7]=0;
1132  fJetInfoVar[11]=-1;
1133  }
1134  }
1135 
1136  if(!fTruthJet) fJetInfoVar[2]=SymParam;
1137  else fJetInfoVar[3]=SymParam;
1138  if(!fTruthJet) fJetInfoVar[12] = DeltaR;
1139  else fJetInfoVar[13] = DeltaR;
1140  //if(!fTruthJet) fJetInfoVar[14]=finaljet.structure_of<fastjet::contrib::SoftDrop>().dropped_count();
1141  //else fJetInfoVar[15]=finaljet.structure_of<fastjet::contrib::SoftDrop>().dropped_count();
1142  if(!(fJetShapeSub==kConstSub)){
1143  if(!fTruthJet) fJetInfoVar[16]=(finaljet.perp()-(GetRhoVal(0)*fJet->Area()));
1144  else fJetInfoVar[17]=finaljet.perp();
1145  }
1146  else{
1147  if(!fTruthJet) fJetInfoVar[16]=(finaljet.perp());
1148  else fJetInfoVar[17]=(finaljet.perp());
1149  }
1150  if(!fTruthJet) fJetInfoVar[18]=(finaljet.m());
1151  else fJetInfoVar[19]=(finaljet.m());
1152 
1153 } catch (fastjet::Error) {
1154  AliError(" [w] FJ Exception caught.");
1155  //return -1;
1156 }
1157  return;
1158 
1159 
1160 }
1161 
1162 //_________________________________________________________________________
1164  try {
1165  std::vector<fastjet::PseudoJet> fInputVectors;
1166  fInputVectors.clear();
1167  fastjet::PseudoJet PseudoTracks;
1168  double xflagalgo=0;
1169  AliParticleContainer *fTrackCont = fJetCont->GetParticleContainer();
1170 
1171  if (fTrackCont) for (Int_t i=0; i<fJet->GetNumberOfTracks(); i++) {
1172  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
1173  if (!fTrk) continue;
1174  PseudoTracks.reset(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
1175  PseudoTracks.set_user_index(fJet->TrackAt(i)+100);
1176  fInputVectors.push_back(PseudoTracks);
1177 
1178  }
1179  fastjet::JetAlgorithm jetalgo(fastjet::antikt_algorithm);
1180 
1181  if(ReclusterAlgo==0){ xflagalgo=0.5;
1182  jetalgo=fastjet::kt_algorithm ;}
1183  if(ReclusterAlgo==1){ xflagalgo=1.5;
1184  jetalgo=fastjet::cambridge_algorithm;}
1185  if(ReclusterAlgo==2){ xflagalgo=2.5;
1186  jetalgo=fastjet::antikt_algorithm;}
1187 
1188  double ndepth=0;
1189  double nhard=0;
1190  double nall=0;
1191 
1192  fastjet::JetDefinition fJetDef(jetalgo, 1., static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
1193  fastjet::ClusterSequence fClustSeqSA(fInputVectors, fJetDef);
1194  std::vector<fastjet::PseudoJet> fOutputJets;
1195  fOutputJets.clear();
1196  fOutputJets=fClustSeqSA.inclusive_jets(0);
1197 
1198  fastjet::PseudoJet jj;
1199  fastjet::PseudoJet j1;
1200  fastjet::PseudoJet j2;
1201  jj=fOutputJets[0];
1202 
1203  while(jj.has_parents(j1,j2)){
1204  ndepth=ndepth+1;
1205  nall=nall+1;
1206  if(j1.perp() < j2.perp()) swap(j1,j2);
1207  double delta_R=j1.delta_R(j2);
1208  double z=j2.perp()/(j1.perp()+j2.perp());
1209  if(z > 0.1) nhard=nhard+1;
1210  double y =log(1.0/delta_R);
1211  double lnpt_rel=log(z*delta_R);
1212  Double_t LundEntries[5] = {y,lnpt_rel,fOutputJets[0].perp(),xflagalgo,ndepth};
1213  if(!bTruth) fhLundIterative->Fill(LundEntries);
1214  else if(bTruth) fhLundIterativeTrue->Fill(LundEntries);
1215  jj=j1;
1216  if(bMinSubjetPt && jj.perp() < fMinSubjetPt) break;
1217  }
1218 
1219 
1220 
1221 
1222 
1223  if(ReclusterAlgo==1){
1224  if(!bTruth) fJetInfoVar[14]=nhard;
1225  else if(bTruth) fJetInfoVar[15]=nhard;
1226  if(!bTruth) fJetInfoVar[20] = nall;
1227  else if(bTruth) fJetInfoVar[21] = nall;
1228  }
1229 
1230 
1231 }
1232 catch (fastjet::Error) {
1233  AliError(" [w] FJ Exception caught.");
1234  //return -1;
1235 }
1236 
1237  return;
1238 
1239 
1240 }
1241 
1242 
1243 //--------------------------------------------------------------------------
1245  std::vector< fastjet::PseudoJet > constituents = jet.constituents();
1246  fastjet::PseudoJet leadingtrack = constituents[0];
1247  Double_t track_pt;
1248  for(size_t i=0; i<constituents.size(); i++){
1249  track_pt=constituents[i].perp();
1250  if (track_pt > leadingtrack.perp()){
1251  leadingtrack = constituents[i];
1252  }
1253  }
1254  return leadingtrack.perp();
1255 }
1256 
1257 
1258 //--------------------------------------------------------------------------
1260 
1261  AliTrackContainer *PartCont = NULL;
1262  AliParticleContainer *PartContMC = NULL;
1263 
1264 
1265  if (fJetShapeSub==kConstSub){
1267  else PartCont = GetTrackContainer(1);
1268  }
1269  else{
1271  else PartCont = GetTrackContainer(0);
1272  }
1273 
1274  TClonesArray *TracksArray = NULL;
1275  TClonesArray *TracksArrayMC = NULL;
1276  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kGenOnTheFly) TracksArrayMC = PartContMC->GetArray();
1277  else TracksArray = PartCont->GetArray();
1278 
1280  if(!PartContMC || !TracksArrayMC) return -99999;
1281  }
1282  else {
1283  if(!PartCont || !TracksArray) return -99999;
1284  }
1285 
1286  AliAODTrack *Track = 0x0;
1287  Int_t Trigger_Index[100];
1288  for (Int_t i=0; i<100; i++) Trigger_Index[i] = 0;
1289  Int_t Trigger_Counter = 0;
1290  Int_t NTracks=0;
1291  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kGenOnTheFly) NTracks = TracksArrayMC->GetEntriesFast();
1292  else NTracks = TracksArray->GetEntriesFast();
1293  for(Int_t i=0; i < NTracks; i++){
1295  if((Track = static_cast<AliAODTrack*>(PartContMC->GetAcceptParticle(i)))){
1296  if (!Track) continue;
1297  if(TMath::Abs(Track->Eta())>0.9) continue;
1298  if (Track->Pt()<0.15) continue;
1299  if ((Track->Pt() >= PtMin) && (Track->Pt()< PtMax)) {
1300  Trigger_Index[Trigger_Counter] = i;
1301  Trigger_Counter++;
1302  }
1303  }
1304  }
1305  else{
1306  if((Track = static_cast<AliAODTrack*>(PartCont->GetAcceptTrack(i)))){
1307  if (!Track) continue;
1308  if(TMath::Abs(Track->Eta())>0.9) continue;
1309  if (Track->Pt()<0.15) continue;
1310  if ((Track->Pt() >= PtMin) && (Track->Pt()< PtMax)) {
1311  Trigger_Index[Trigger_Counter] = i;
1312  Trigger_Counter++;
1313  }
1314  }
1315  }
1316  }
1317  if (Trigger_Counter == 0) return -99999;
1318  Int_t RandomNumber = 0, Index = 0 ;
1319  TRandom3 Random (0);
1320  RandomNumber = Random.Integer(Trigger_Counter);
1321  Index = Trigger_Index[RandomNumber];
1322  return Index;
1323 }
1324 
1325 
1326 //________________________________________________________________________
1328 {
1329  AliEmcalJet *jet2 = jet1->ClosestJet();
1330  if (!jet2) return -1;
1331 
1332  Double_t jetPt2 = jet2->Pt();
1333  if (jetPt2 <= 0) return -1;
1334 
1335  Int_t bgeom = kTRUE;
1336  if (!cont2) bgeom = kFALSE;
1337  Double_t sumPt = 0.;
1338  AliVParticle *vpf = 0x0;
1339  Int_t iFound = 0;
1340  for (Int_t icc = 0; icc < jet2->GetNumberOfTracks(); icc++) {
1341  Int_t idx = (Int_t)jet2->TrackAt(icc);
1342  //get particle
1343  AliVParticle *p2 = 0x0;
1344  if (bgeom) p2 = static_cast<AliVParticle*>(jet2->TrackAt(icc, cont2->GetArray()));
1345  iFound = 0;
1346  for (Int_t icf = 0; icf < jet1->GetNumberOfTracks(); icf++) {
1347  if (!bgeom && idx == jet1->TrackAt(icf) && iFound == 0 ) {
1348  iFound = 1;
1349  vpf = jet1->Track(icf);
1350  if (vpf) sumPt += vpf->Pt();
1351  continue;
1352  }
1353  if (bgeom){
1354  vpf = jet1->Track(icf);
1355  if (!vpf) continue;
1356  if (!SameParticle(vpf, p2, 1.e-4)) continue; //not the same particle
1357  sumPt += vpf->Pt();
1358  }
1359  }
1360  }
1361 
1362  Double_t fraction = sumPt / jetPt2;
1363 
1364  return fraction;
1365 }
1366 
1367 //________________________________________________________________________
1368 Bool_t AliAnalysisTaskRecoilJetYield::SameParticle(const AliVParticle* part1, const AliVParticle* part2, Double_t dist)
1369 {
1370  if(!part1) return kFALSE;
1371  if(!part2) return kFALSE;
1372  Double_t dPhi = TMath::Abs(part1->Phi() - part2->Phi());
1373  dPhi = TVector2::Phi_mpi_pi(dPhi);
1374  if (dPhi > dist) return kFALSE;
1375  return kTRUE;
1376 }
1377 
1378 //________________________________________________________________________
1380  //
1381  // retrieve event objects
1382  //
1384  return kFALSE;
1385 
1386  return kTRUE;
1387 }
1388 
1389 
1390 //_______________________________________________________________________
1392 {
1393  // Called once at the end of the analysis.
1394 
1395 
1396 }
Bool_t RetrieveEventObjects()
Retrieve common objects from event.
Double_t Area() const
Definition: AliEmcalJet.h:130
double Double_t
Definition: External.C:58
Double_t PTD(AliEmcalJet *Jet, Int_t JetContNb)
static Double_t DeltaPhi(Double_t phia, Double_t phib, Double_t rMin=-TMath::Pi()/2, Double_t rMax=3 *TMath::Pi()/2)
Definition: External.C:236
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:341
AliJetContainer * GetJetContainer(Int_t i=0) const
Bool_t FillHistograms()
Function filling histograms.
Double_t Eta() const
Definition: AliEmcalJet.h:121
Double_t Phi_Upper
Double_t Phi() const
Definition: AliEmcalJet.h:117
Container with name, TClonesArray and cuts for particles.
Double_t fEPV0
!event plane V0
Int_t GetLabel() const
Definition: AliEmcalJet.h:124
AliVParticle * Track(Int_t idx) const
virtual Bool_t AcceptJet(Int_t i, UInt_t &rejectionReason) const
Double_t Eta_Lower
Float_t Percentile_1
Double_t Eta_Upper
void SoftDrop(AliEmcalJet *fJet, AliJetContainer *fJetCont, double zcut, double beta, Bool_t fTruthJet)
Double_t Phi_Lower
UShort_t GetNumberOfConstituents() const
Definition: AliEmcalJet.h:140
Container for particles within the EMCAL framework.
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:160
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:139
TString kData
Declare data MC or deltaAOD.
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
AliParticleContainer * GetParticleContainer() const
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
Int_t SelectTriggerHadron(Float_t PtMin, Float_t PtMax)
Double_t GetFractionSharedPt_SubMatching(const AliEmcalJet *jet, AliParticleContainer *cont2=0x0) const
Int_t GetNJets() const
virtual AliVParticle * GetAcceptParticle(Int_t i=-1) const
Double_t fCent
!event centrality
const Int_t nBranch
Double_t Angularity(AliEmcalJet *Jet, Int_t JetContNb)
AliEmcalJet * GetNextAcceptJet()
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)
Double_t RelativePhi(Double_t Phi1, Double_t Phi2)
static Bool_t SameParticle(const AliVParticle *part1, const AliVParticle *part2, Double_t dist=1.e-4)
Double_t Pt() const
Definition: AliEmcalJet.h:109
Double_t GetRhoVal(Int_t i=0) const
void Track()
AliEmcalList * fOutput
!output list
AliTrackContainer * GetTrackContainer(Int_t i=0) const
Double_t DeltaR(const AliVParticle *part1, const AliVParticle *part2)
void SetMakeGeneralHistograms(Bool_t g)
void RecursiveParents(AliEmcalJet *fJet, AliJetContainer *fJetCont, Int_t ReclusterAlgo, Bool_t bTrue)
virtual AliVTrack * GetAcceptTrack(Int_t i=-1) const
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
void SetNumberOfTracks(Int_t n)
Definition: AliEmcalJet.h:266
unsigned short UShort_t
Definition: External.C:28
const char Option_t
Definition: External.C:48
void UserCreateOutputObjects()
Main initialization function on the worker.
bool Bool_t
Definition: External.C:53
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.
void swap(PWGJE::EMCALJetTasks::AliAnalysisTaskEmcalJetHPerformance &first, PWGJE::EMCALJetTasks::AliAnalysisTaskEmcalJetHPerformance &second)
Double_t RelativePhiEventPlane(Double_t EventPlane, Double_t Phi)
Double_t LeadingTrackPt(fastjet::PseudoJet jet)
Double_t M() const
Definition: AliEmcalJet.h:120
Container for jet within the EMCAL jet framework.
AliEmcalJet * GetJet(Int_t i) const