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