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