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