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