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