AliPhysics  a8afd6c (a8afd6c)
AliAnalysisTaskSubJetFraction.cxx
Go to the documentation of this file.
1 //
2 // Subjet fraction analysis task.
3 //
4 //Nima Zardoshti
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"
45 #include "AliAnalysisTaskSubJetFraction.h"
46 #include "AliEmcalPythiaInfo.h"
47 
48 #include "FJ_includes.h"
49 
50 //Globals
51 
52 
53 Double_t XBinsPt[66]={0.15,0.20,0.25,0.30,0.35,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.9,0.95,1.00,1.10,1.20,1.30,1.40,1.50,1.60,1.70,1.80,1.90,2.00,2.20,2.40,2.60,2.80,3.00,3.20,3.40,3.60,3.80,4.00,4.50,5.00,5.50,6.00,6.50,7.00,8.00,9.00,10.00,11.00,12.00,13.00,14.00,15.00,16.00,18.00,20.00,22.00,24.00,26.00,28.00,30.00,32.00,34.00,36.00,40.00,45.00,50.00}; //Lower edges of new bins for rebbing purposes for track pt
54 Double_t XBinsPtSize=66;
55 Double_t XBinsJetPt[35]={0,0.50,1.00,2.00,3.00,4.00,5.00,6.00,7.00,8.00,9.00,10.00,12.00,14.00,16.00,18.00,20.00,25.00,30.00,35.00,40.00,45.00,50.00,60.00,70.00,80.00,90.00,100.00,120.00,140.00,160.00,180.00,200.00,250.00,300.00}; //for jet pt
56 Int_t XBinsJetPtSize=35;
57 Double_t XBinsJetMass[28]={0,0.50,1.00,2.00,3.00,4.00,5.00,6.00,7.00,8.00,9.00,10.00,12.00,14.00,16.00,18.00,20.00,25.00,30.00,35.00,40.00,45.00,50.00,60.00,70.00,80.00,90.00,100.00}; //for jet mass
58 Int_t XBinsJetMassSize=28;
59 Double_t Eta_Up=1.00;
60 Double_t Eta_Low=-1.00;
61 Int_t Eta_Bins=100;
62 Double_t Phi_Up=2*(TMath::Pi());
63 Double_t Phi_Low=(-1*(TMath::Pi()));
64 Int_t Phi_Bins=100;
65 
66 
67 
68 
69 
70 using std::cout;
71 using std::endl;
72 
73 ClassImp(AliAnalysisTaskSubJetFraction)
74 
75 //________________________________________________________________________
76 AliAnalysisTaskSubJetFraction::AliAnalysisTaskSubJetFraction() :
77  AliAnalysisTaskEmcalJet("AliAnalysisTaskSubJetFraction", kTRUE),
78  fContainer(0),
79  fMinFractionShared(0),
80  fJetShapeType(kData),
81  fJetShapeSub(kNoSub),
82  fJetSelection(kInclusive),
83  fPtThreshold(-9999.),
84  fRMatching(0.2),
85  fPtMinTriggerHadron(20.),
86  fPtMaxTriggerHadron(50.),
87  fRecoilAngularWindow(0.6),
88  fSemigoodCorrect(0),
89  fHolePos(0),
90  fHoleWidth(0),
91  fRandom(0x0),
92  fCentSelectOn(kTRUE),
93  fCentMin(0),
94  fCentMax(10),
95  fSubJetAlgorithm(0),
96  fSubJetRadius(0.1),
97  fSubJetMinPt(1),
98  fJetRadius(0.4),
99  fRMatched(0.2),
100  fSharedFractionPtMin(0.5),
101  fDerivSubtrOrder(0),
102  fFullTree(kFALSE),
103  fBeta_SD(0),
104  fZCut(0.1),
105  fReclusteringAlgorithm(0),
106  fSoftDropOn(0),
107  fMLOn(0),
108  fNsubMeasure(kFALSE),
109  fRandomisationEqualPt(kFALSE),
110  fhPtTriggerHadron(0x0),
111  fhJetPt(0x0),
112  fhJetPt_1(0x0),
113  fhJetPt_2(0x0),
114  fhJetPhi(0x0),
115  fhJetPhi_1(0x0),
116  fhJetPhi_2(0x0),
117  fhJetEta(0x0),
118  fhJetEta_1(0x0),
119  fhJetEta_2(0x0),
120  fhJetMass(0x0),
121  fhJetMass_1(0x0),
122  fhJetMass_2(0x0),
123  fhJetRadius(0x0),
124  fhJetRadius_1(0x0),
125  fhJetRadius_2(0x0),
126  fhJetCounter(0x0),
127  fhJetCounter_1(0x0),
128  fhJetCounter_2(0x0),
129  fhNumberOfJetTracks(0x0),
130  fhNumberOfJetTracks_1(0x0),
131  fhNumberOfJetTracks_2(0x0),
132  fhSubJetPt(0x0),
133  fhSubJetPt_1(0x0),
134  fhSubJetPt_2(0x0),
135  fhSubJetMass(0x0),
136  fhSubJetMass_1(0x0),
137  fhSubJetMass_2(0x0),
138  fhSubJetRadius(0x0),
139  fhSubJetRadius_1(0x0),
140  fhSubJetRadius_2(0x0),
141  fhSubJetCounter(0x0),
142  fhSubJetCounter_1(0x0),
143  fhSubJetCounter_2(0x0),
144  fhNumberOfSubJetTracks(0x0),
145  fhNumberOfSubJetTracks_1(0x0),
146  fhNumberOfSubJetTracks_2(0x0),
147  fh2PtTriggerHadronJet(0x0),
148  fhPhiTriggerHadronJet(0x0),
149  fhPhiTriggerHadronEventPlane(0x0),
150  fhPhiTriggerHadronEventPlaneTPC(0x0),
151  fhTrackPhi(0x0),
152  fhTrackPhi_Cut(0x0),
153  fh2PtRatio(0x0),
154  fhEventCounter(0x0),
155  fhEventCounter_1(0x0),
156  fhEventCounter_2(0x0),
157  fhSubJettiness1CheckRatio_FJ_AKT(0x0),
158  fhSubJettiness1CheckRatio_FJ_KT(0x0),
159  fhSubJettiness1CheckRatio_FJ_CA(0x0),
160  fhSubJettiness1CheckRatio_FJ_WTA_KT(0x0),
161  fhSubJettiness1CheckRatio_FJ_WTA_CA(0x0),
162  fhSubJettiness1CheckRatio_FJ_OP_AKT(0x0),
163  fhSubJettiness1CheckRatio_FJ_OP_KT(0x0),
164  fhSubJettiness1CheckRatio_FJ_OP_CA(0x0),
165  fhSubJettiness1CheckRatio_FJ_OP_WTA_KT(0x0),
166  fhSubJettiness1CheckRatio_FJ_OP_WTA_CA(0x0),
167  fhSubJettiness1CheckRatio_FJ_MIN(0x0),
168  fhSubJettiness2CheckRatio_FJ_AKT(0x0),
169  fhSubJettiness2CheckRatio_FJ_KT(0x0),
170  fhSubJettiness2CheckRatio_FJ_CA(0x0),
171  fhSubJettiness2CheckRatio_FJ_WTA_KT(0x0),
172  fhSubJettiness2CheckRatio_FJ_WTA_CA(0x0),
173  fhSubJettiness2CheckRatio_FJ_OP_AKT(0x0),
174  fhSubJettiness2CheckRatio_FJ_OP_KT(0x0),
175  fhSubJettiness2CheckRatio_FJ_OP_CA(0x0),
176  fhSubJettiness2CheckRatio_FJ_OP_WTA_KT(0x0),
177  fhSubJettiness2CheckRatio_FJ_OP_WTA_CA(0x0),
178  fhSubJettiness2CheckRatio_FJ_MIN(0x0),
179  fhSubJettiness2to1CheckRatio_FJ_AKT(0x0),
180  fhSubJettiness2to1CheckRatio_FJ_KT(0x0),
181  fhSubJettiness2to1CheckRatio_FJ_CA(0x0),
182  fhSubJettiness2to1CheckRatio_FJ_WTA_KT(0x0),
183  fhSubJettiness2to1CheckRatio_FJ_WTA_CA(0x0),
184  fhSubJettiness2to1CheckRatio_FJ_OP_AKT(0x0),
185  fhSubJettiness2to1CheckRatio_FJ_OP_KT(0x0),
186  fhSubJettiness2to1CheckRatio_FJ_OP_CA(0x0),
187  fhSubJettiness2to1CheckRatio_FJ_OP_WTA_KT(0x0),
188  fhSubJettiness2to1CheckRatio_FJ_OP_WTA_CA(0x0),
189  fhSubJettiness2to1CheckRatio_FJ_MIN(0x0),
190  fhSubJettiness1_FJ_AKT(0x0),
191  fhSubJettiness1_FJ_KT(0x0),
192  fhSubJettiness1_FJ_CA(0x0),
193  fhSubJettiness1_FJ_WTA_KT(0x0),
194  fhSubJettiness1_FJ_WTA_CA(0x0),
195  fhSubJettiness1_FJ_OP_AKT(0x0),
196  fhSubJettiness1_FJ_OP_KT(0x0),
197  fhSubJettiness1_FJ_OP_CA(0x0),
198  fhSubJettiness1_FJ_OP_WTA_KT(0x0),
199  fhSubJettiness1_FJ_OP_WTA_CA(0x0),
200  fhSubJettiness1_FJ_MIN(0x0),
201  fhSubJettiness2_FJ_AKT(0x0),
202  fhSubJettiness2_FJ_KT(0x0),
203  fhSubJettiness2_FJ_CA(0x0),
204  fhSubJettiness2_FJ_WTA_KT(0x0),
205  fhSubJettiness2_FJ_WTA_CA(0x0),
206  fhSubJettiness2_FJ_OP_AKT(0x0),
207  fhSubJettiness2_FJ_OP_KT(0x0),
208  fhSubJettiness2_FJ_OP_CA(0x0),
209  fhSubJettiness2_FJ_OP_WTA_KT(0x0),
210  fhSubJettiness2_FJ_OP_WTA_CA(0x0),
211  fhSubJettiness2_FJ_MIN(0x0),
212  fhSubJettiness2to1_FJ_AKT(0x0),
213  fhSubJettiness2to1_FJ_KT(0x0),
214  fhSubJettiness2to1_FJ_CA(0x0),
215  fhSubJettiness2to1_FJ_WTA_KT(0x0),
216  fhSubJettiness2to1_FJ_WTA_CA(0x0),
217  fhSubJettiness2to1_FJ_OP_AKT(0x0),
218  fhSubJettiness2to1_FJ_OP_KT(0x0),
219  fhSubJettiness2to1_FJ_OP_CA(0x0),
220  fhSubJettiness2to1_FJ_OP_WTA_KT(0x0),
221  fhSubJettiness2to1_FJ_OP_WTA_CA(0x0),
222  fhSubJettiness2to1_FJ_MIN(0x0),
223  fShapesVar_Tracks_Rec(0),
224  fShapesVar_Tracks_Truth(0),
225  fTreeResponseMatrixAxis(0),
226  fTreeTracks(0)
227 
228 {
229  for(Int_t i=0;i<nVar;i++){
230  fShapesVar[i]=0;
231  }
232  SetMakeGeneralHistograms(kTRUE);
233  DefineOutput(1, TList::Class());
234  DefineOutput(2, TTree::Class());
235  DefineOutput(3, TTree::Class());
236 }
237 
238 //________________________________________________________________________
239 AliAnalysisTaskSubJetFraction::AliAnalysisTaskSubJetFraction(const char *name) :
240  AliAnalysisTaskEmcalJet(name, kTRUE),
241  fContainer(0),
242  fMinFractionShared(0),
243  fJetShapeType(kData),
244  fJetShapeSub(kNoSub),
245  fJetSelection(kInclusive),
246  fPtThreshold(-9999.),
247  fRMatching(0.2),
248  fPtMinTriggerHadron(20.),
249  fPtMaxTriggerHadron(50.),
250  fRecoilAngularWindow(0.6),
251  fSemigoodCorrect(0),
252  fHolePos(0),
253  fHoleWidth(0),
254  fRandom(0x0),
255  fCentSelectOn(kTRUE),
256  fCentMin(0),
257  fCentMax(10),
258  fSubJetAlgorithm(0),
259  fSubJetRadius(0.1),
260  fSubJetMinPt(1),
261  fJetRadius(0.4),
262  fRMatched(0.2),
263  fSharedFractionPtMin(0.5),
264  fDerivSubtrOrder(0),
265  fFullTree(kFALSE),
266  fBeta_SD(0),
267  fZCut(0.1),
268  fReclusteringAlgorithm(0),
269  fSoftDropOn(0),
270  fMLOn(0),
271  fNsubMeasure(kFALSE),
272  fRandomisationEqualPt(kFALSE),
273  fhPtTriggerHadron(0x0),
274  fhJetPt(0x0),
275  fhJetPt_1(0x0),
276  fhJetPt_2(0x0),
277  fhJetPhi(0x0),
278  fhJetPhi_1(0x0),
279  fhJetPhi_2(0x0),
280  fhJetEta(0x0),
281  fhJetEta_1(0x0),
282  fhJetEta_2(0x0),
283  fhJetMass(0x0),
284  fhJetMass_1(0x0),
285  fhJetMass_2(0x0),
286  fhJetRadius(0x0),
287  fhJetRadius_1(0x0),
288  fhJetRadius_2(0x0),
289  fhJetCounter(0x0),
290  fhJetCounter_1(0x0),
291  fhJetCounter_2(0x0),
292  fhNumberOfJetTracks(0x0),
293  fhNumberOfJetTracks_1(0x0),
294  fhNumberOfJetTracks_2(0x0),
295  fhSubJetPt(0x0),
296  fhSubJetPt_1(0x0),
297  fhSubJetPt_2(0x0),
298  fhSubJetMass(0x0),
299  fhSubJetMass_1(0x0),
300  fhSubJetMass_2(0x0),
301  fhSubJetRadius(0x0),
302  fhSubJetRadius_1(0x0),
303  fhSubJetRadius_2(0x0),
304  fhSubJetCounter(0x0),
305  fhSubJetCounter_1(0x0),
306  fhSubJetCounter_2(0x0),
307  fhNumberOfSubJetTracks(0x0),
308  fhNumberOfSubJetTracks_1(0x0),
309  fhNumberOfSubJetTracks_2(0x0),
310  fh2PtTriggerHadronJet(0x0),
311  fhPhiTriggerHadronJet(0x0),
312  fhPhiTriggerHadronEventPlane(0x0),
313  fhPhiTriggerHadronEventPlaneTPC(0x0),
314  fhTrackPhi(0x0),
315  fhTrackPhi_Cut(0x0),
316  fh2PtRatio(0x0),
317  fhEventCounter(0x0),
318  fhEventCounter_1(0x0),
319  fhEventCounter_2(0x0),
320  fhSubJettiness1CheckRatio_FJ_AKT(0x0),
321  fhSubJettiness1CheckRatio_FJ_KT(0x0),
322  fhSubJettiness1CheckRatio_FJ_CA(0x0),
323  fhSubJettiness1CheckRatio_FJ_WTA_KT(0x0),
324  fhSubJettiness1CheckRatio_FJ_WTA_CA(0x0),
325  fhSubJettiness1CheckRatio_FJ_OP_AKT(0x0),
326  fhSubJettiness1CheckRatio_FJ_OP_KT(0x0),
327  fhSubJettiness1CheckRatio_FJ_OP_CA(0x0),
328  fhSubJettiness1CheckRatio_FJ_OP_WTA_KT(0x0),
329  fhSubJettiness1CheckRatio_FJ_OP_WTA_CA(0x0),
330  fhSubJettiness1CheckRatio_FJ_MIN(0x0),
331  fhSubJettiness2CheckRatio_FJ_AKT(0x0),
332  fhSubJettiness2CheckRatio_FJ_KT(0x0),
333  fhSubJettiness2CheckRatio_FJ_CA(0x0),
334  fhSubJettiness2CheckRatio_FJ_WTA_KT(0x0),
335  fhSubJettiness2CheckRatio_FJ_WTA_CA(0x0),
336  fhSubJettiness2CheckRatio_FJ_OP_AKT(0x0),
337  fhSubJettiness2CheckRatio_FJ_OP_KT(0x0),
338  fhSubJettiness2CheckRatio_FJ_OP_CA(0x0),
339  fhSubJettiness2CheckRatio_FJ_OP_WTA_KT(0x0),
340  fhSubJettiness2CheckRatio_FJ_OP_WTA_CA(0x0),
341  fhSubJettiness2CheckRatio_FJ_MIN(0x0),
342  fhSubJettiness2to1CheckRatio_FJ_AKT(0x0),
343  fhSubJettiness2to1CheckRatio_FJ_KT(0x0),
344  fhSubJettiness2to1CheckRatio_FJ_CA(0x0),
345  fhSubJettiness2to1CheckRatio_FJ_WTA_KT(0x0),
346  fhSubJettiness2to1CheckRatio_FJ_WTA_CA(0x0),
347  fhSubJettiness2to1CheckRatio_FJ_OP_AKT(0x0),
348  fhSubJettiness2to1CheckRatio_FJ_OP_KT(0x0),
349  fhSubJettiness2to1CheckRatio_FJ_OP_CA(0x0),
350  fhSubJettiness2to1CheckRatio_FJ_OP_WTA_KT(0x0),
351  fhSubJettiness2to1CheckRatio_FJ_OP_WTA_CA(0x0),
352  fhSubJettiness2to1CheckRatio_FJ_MIN(0x0),
353  fhSubJettiness1_FJ_AKT(0x0),
354  fhSubJettiness1_FJ_KT(0x0),
355  fhSubJettiness1_FJ_CA(0x0),
356  fhSubJettiness1_FJ_WTA_KT(0x0),
357  fhSubJettiness1_FJ_WTA_CA(0x0),
358  fhSubJettiness1_FJ_OP_AKT(0x0),
359  fhSubJettiness1_FJ_OP_KT(0x0),
360  fhSubJettiness1_FJ_OP_CA(0x0),
361  fhSubJettiness1_FJ_OP_WTA_KT(0x0),
362  fhSubJettiness1_FJ_OP_WTA_CA(0x0),
363  fhSubJettiness1_FJ_MIN(0x0),
364  fhSubJettiness2_FJ_AKT(0x0),
365  fhSubJettiness2_FJ_KT(0x0),
366  fhSubJettiness2_FJ_CA(0x0),
367  fhSubJettiness2_FJ_WTA_KT(0x0),
368  fhSubJettiness2_FJ_WTA_CA(0x0),
369  fhSubJettiness2_FJ_OP_AKT(0x0),
370  fhSubJettiness2_FJ_OP_KT(0x0),
371  fhSubJettiness2_FJ_OP_CA(0x0),
372  fhSubJettiness2_FJ_OP_WTA_KT(0x0),
373  fhSubJettiness2_FJ_OP_WTA_CA(0x0),
374  fhSubJettiness2_FJ_MIN(0x0),
375  fhSubJettiness2to1_FJ_AKT(0x0),
376  fhSubJettiness2to1_FJ_KT(0x0),
377  fhSubJettiness2to1_FJ_CA(0x0),
378  fhSubJettiness2to1_FJ_WTA_KT(0x0),
379  fhSubJettiness2to1_FJ_WTA_CA(0x0),
380  fhSubJettiness2to1_FJ_OP_AKT(0x0),
381  fhSubJettiness2to1_FJ_OP_KT(0x0),
382  fhSubJettiness2to1_FJ_OP_CA(0x0),
383  fhSubJettiness2to1_FJ_OP_WTA_KT(0x0),
384  fhSubJettiness2to1_FJ_OP_WTA_CA(0x0),
385  fhSubJettiness2to1_FJ_MIN(0x0),
386  fShapesVar_Tracks_Rec(0),
387  fShapesVar_Tracks_Truth(0),
388  fTreeResponseMatrixAxis(0),
389  fTreeTracks(0)
390 {
391  // Standard constructor.
392  for(Int_t i=0;i<nVar;i++){
393  fShapesVar[i]=0;
394  }
395  SetMakeGeneralHistograms(kTRUE);
396  DefineOutput(1, TList::Class());
397  DefineOutput(2, TTree::Class());
398  DefineOutput(3, TTree::Class());
399 }
400 
401 //________________________________________________________________________
402 AliAnalysisTaskSubJetFraction::~AliAnalysisTaskSubJetFraction()
403 {
404  // Destructor.
405 }
406 
407 //________________________________________________________________________
408  void AliAnalysisTaskSubJetFraction::UserCreateOutputObjects()
409 {
410  // Create user output.
411 
412  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
413 
414  Bool_t oldStatus = TH1::AddDirectoryStatus();
415  TH1::AddDirectory(kFALSE);
416  TH1::AddDirectory(oldStatus);
417  //create a tree used for the MC data and making a 4D response matrix
418 
419 
420 
421 
422  const char* nameoutput_1 = GetOutputSlot(3)->GetContainer()->GetName();
423  fTreeTracks = new TTree(nameoutput_1, nameoutput_1);
424  TString *fShapesVarNames_Tracks=new TString[2];
425  fShapesVarNames_Tracks[0] = "Jet_Constituents_Det";
426  fShapesVarNames_Tracks[1] = "Jet_Constituents_Truth";
427  fTreeTracks->Branch(fShapesVarNames_Tracks[0].Data(), &fShapesVar_Tracks_Rec, 0,1);
428  fTreeTracks->Branch(fShapesVarNames_Tracks[1].Data(), &fShapesVar_Tracks_Truth, 0,1);
429  //fTreeTracks->Branch(fShapesVarNames_Tracks[0].Data(), &fShapesVar_Tracks[0], "dfdf");
430  //fTreeTracks->Branch(fShapesVarNames_Tracks[1].Data(), &fShapesVar_Tracks[1], "fdss");
431 
432 
433 
434  const char* nameoutput = GetOutputSlot(2)->GetContainer()->GetName();
435  fTreeResponseMatrixAxis = new TTree(nameoutput, nameoutput);
436  if (fFullTree){
437  TString *fShapesVarNames = new TString [nVar];
438  if (fMLOn==0){
439  fShapesVarNames[0] = "Pt";
440  fShapesVarNames[1] = "Pt_Truth";
441  fShapesVarNames[2] = "Tau1";
442  fShapesVarNames[3] = "Tau1_Truth";
443  fShapesVarNames[4] = "Tau2";
444  fShapesVarNames[5] = "Tau2_Truth";
445  fShapesVarNames[6] = "Tau3";
446  fShapesVarNames[7] = "Tau3_Truth";
447  fShapesVarNames[8] = "OpeningAngle";
448  fShapesVarNames[9] = "OpeningAngle_Truth";
449  fShapesVarNames[10] = "JetMultiplicity";
450  fShapesVarNames[11] = "JetMultiplicity_Truth";
451  fShapesVarNames[12] = "OpeningAngleSD";
452  fShapesVarNames[13] = "OpeningAngleSD_Truth";
453  fShapesVarNames[14] = "Zg";
454  fShapesVarNames[15] = "Zg_Truth";
455  fShapesVarNames[16] = "LeadingTrackPt";
456  fShapesVarNames[17] = "LeadingTrackPt_Truth";
457  fShapesVarNames[18] = "EventPlaneTriggerHadron";
458  fShapesVarNames[19] = "EventPlaneTriggerHadron_Truth";
459  fShapesVarNames[20] = "EventPlaneTPCTriggerHadron";
460  fShapesVarNames[21] = "EventPlaneTPCTriggerHadron_Truth";
461  fShapesVarNames[22] = "DeltaR";
462  fShapesVarNames[23] = "DeltaR_Truth";
463  fShapesVarNames[24] = "Frac1";
464  fShapesVarNames[25] = "Frac1_Truth";
465  fShapesVarNames[26] = "Frac2";
466  fShapesVarNames[27] = "Frac2_Truth";
467  }
468  if (fMLOn==1){
469  fShapesVarNames[0] = "Pt_Rec";
470  fShapesVarNames[1] = "Pt_Truth";
471  fShapesVarNames[2] = "Eta_Rec";
472  fShapesVarNames[3] = "Eta_Truth";
473  fShapesVarNames[4] = "Phi_Rec";
474  fShapesVarNames[5] = "Phi_Truth";
475  fShapesVarNames[6] = "Mass_Rec";
476  fShapesVarNames[7] = "Mass_Truth";
477  fShapesVarNames[8] = "JetMultiplicity_Rec";
478  fShapesVarNames[9] = "JetMultiplicity_Truth";
479  fShapesVarNames[10] = "Parton_1_Flag";
480  fShapesVarNames[11] = "Blank_2";
481  fShapesVarNames[12] = "Parton_1_Eta";
482  fShapesVarNames[13] = "Blank_4";
483  fShapesVarNames[14] = "Parton_1_Phi";
484  fShapesVarNames[15] = "Blank_6";
485  fShapesVarNames[16] = "Parton_2_Flag";
486  fShapesVarNames[17] = "Blank_8";
487  fShapesVarNames[18] = "Parton_2_Eta";
488  fShapesVarNames[19] = "Blank_10";
489  fShapesVarNames[20] = "Parton_2_Phi";
490  fShapesVarNames[21] = "Blank_12";
491  fShapesVarNames[22] = "Blank_13";
492  fShapesVarNames[23] = "Blank_14";
493  fShapesVarNames[24] = "Blank_15";
494  fShapesVarNames[25] = "Blank_16";
495  fShapesVarNames[26] = "Blank_17";
496  fShapesVarNames[27] = "Blank_18";
497  }
498  for(Int_t ivar=0; ivar < nVar; ivar++){
499  cout<<"looping over variables"<<endl;
500  fTreeResponseMatrixAxis->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/D", fShapesVarNames[ivar].Data()));
501  }
502  }
503 
504  if (!fFullTree){
505  const Int_t nVarMin = 22;
506  TString *fShapesVarNames = new TString [nVarMin];
507  if (fMLOn==0){
508  fShapesVarNames[0] = "Pt";
509  fShapesVarNames[1] = "Pt_Truth";
510  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[1] = "Tau1_Randomisation";
511  fShapesVarNames[2] = "Tau1";
512  fShapesVarNames[3] = "Tau1_Truth";
513  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[3] = "Tau1_ExtraProng_01_10";
514  fShapesVarNames[4] = "Tau2";
515  fShapesVarNames[5] = "Tau2_Truth";
516  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[5] = "Tau1_ExtraProng_01_15";
517  fShapesVarNames[6] = "SubJet1LeadingTrackPt";
518  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[6] = "Tau1_ExtraProng_01_20";
519  fShapesVarNames[7] = "SubJet1LeadingTrackPt_Truth";
520  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[7] = "Tau1_ExtraProng_01_25";
521  fShapesVarNames[8] = "OpeningAngle";
522  fShapesVarNames[9] = "OpeningAngle_Truth";
523  //if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[9] = "Tau1_ExtraProng_01_30";
524  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[9] = "Tau1_kTTracks_1_2_1";
525  fShapesVarNames[10] = "SubJet2LeadingTrackPt";
526  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[10] = "Tau1_ExtraProng_03_10";
527  fShapesVarNames[11] = "SubJet2LeadingTrackPt_Truth";
528  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[11] = "Tau1_ExtraProng_03_15";
529  fShapesVarNames[12] = "OpeningAngleSD";
530  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[12] = "Tau1_ExtraProng_03_20";
531  fShapesVarNames[13] = "OpeningAngleSD_Truth";
532  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[13] = "Tau2_Randomisation";
533  fShapesVarNames[14] = "SubJet1Pt";
534  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[14] = "Tau2_ExtraProng_01_10";
535  fShapesVarNames[15] = "SubJet1Pt_Truth";
536  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[15] = "Tau2_ExtraProng_01_15";
537  fShapesVarNames[16] = "LeadingTrackPt";
538  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[16] = "Tau2_ExtraProng_01_20";
539  fShapesVarNames[17] = "LeadingTrackPt_Truth";
540  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[17] = "Tau2_ExtraProng_01_25";
541  fShapesVarNames[18] = "EventPlaneTriggerHadron";
542  // if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[18] = "Tau2_ExtraProng_01_30";
543  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[18] = "Tau2_kTTracks_1_2_1";
544  fShapesVarNames[19] = "EventPlaneTriggerHadron_Truth";
545  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[19] = "Tau2_ExtraProng_03_10";
546  fShapesVarNames[20] = "SubJet2Pt";
547  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[20] = "Tau2_ExtraProng_03_15";
548  fShapesVarNames[21] = "SubJet2Pt_Truth";
549  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[21] = "Tau2_ExtraProng_03_20";
550  }
551  if (fMLOn==1){
552  fShapesVarNames[0] = "Pt_Rec";
553  fShapesVarNames[1] = "Pt_Truth";
554  fShapesVarNames[2] = "Eta_Rec";
555  fShapesVarNames[3] = "Eta_Truth";
556  fShapesVarNames[4] = "Phi_Rec";
557  fShapesVarNames[5] = "Phi_Truth";
558  fShapesVarNames[6] = "Mass_Rec";
559  fShapesVarNames[7] = "Mass_Truth";
560  fShapesVarNames[8] = "JetMultiplicity_Rec";
561  fShapesVarNames[9] = "JetMultiplicity_Truth";
562  fShapesVarNames[10] = "Blank_1";
563  fShapesVarNames[11] = "Blank_2";
564  fShapesVarNames[12] = "Blank_3";
565  fShapesVarNames[13] = "Blank_4";
566  fShapesVarNames[14] = "Blank_5";
567  fShapesVarNames[15] = "Blank_6";
568  fShapesVarNames[16] = "Blank_7";
569  fShapesVarNames[17] = "Blank_8";
570  fShapesVarNames[18] = "Blank_9";
571  fShapesVarNames[19] = "Blank_10";
572  fShapesVarNames[20] = "Blank_11";
573  fShapesVarNames[21] = "Blank_12";
574  }
575 
576  for(Int_t ivar=0; ivar < nVarMin; ivar++){
577  cout<<"looping over variables"<<endl;
578  fTreeResponseMatrixAxis->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/D", fShapesVarNames[ivar].Data()));
579  }
580  }
581 
582  if (fJetSelection == AliAnalysisTaskSubJetFraction::kRecoil){
583  fhPtTriggerHadron= new TH1F("fhPtTriggerHadron", "fhPtTriggerHadron",1500,-0.5,149.5);
584  fOutput->Add(fhPtTriggerHadron);
585  fh2PtTriggerHadronJet= new TH2F("fh2PtTriggerHadronJet", "fh2PtTriggerHadronJet",1500,-0.5,149.5,1500,-0.5,149.5);
586  fOutput->Add(fh2PtTriggerHadronJet);
587  fhPhiTriggerHadronJet= new TH1F("fhPhiTriggerHadronJet", "fhPhiTriggerHadronJet",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
588  fOutput->Add(fhPhiTriggerHadronJet);
589  fhPhiTriggerHadronEventPlane= new TH1F("fhPhiTriggerHadronEventPlane", "fhPhiTriggerHadronEventPlane",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
590  fOutput->Add(fhPhiTriggerHadronEventPlane);
591  fhPhiTriggerHadronEventPlaneTPC= new TH1F("fhPhiTriggerHadronEventPlaneTPC", "fhPhiTriggerHadronEventPlaneTPC",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
592  fOutput->Add(fhPhiTriggerHadronEventPlaneTPC);
593  }
594  if (fJetShapeType==AliAnalysisTaskSubJetFraction::kData || fJetShapeType==AliAnalysisTaskSubJetFraction::kSim || fJetShapeType==AliAnalysisTaskSubJetFraction::kGenOnTheFly){
595 
596  // fhJetPt= new TH1F("fhJetPt", "Jet Pt", (XBinsJetPtSize)-1, XBinsJetPt);
597  fhJetPt= new TH1F("fhJetPt", "Jet Pt",1500,-0.5,149.5 );
598  fOutput->Add(fhJetPt);
599  //fhJetPhi= new TH1F("fhJetPhi", "Jet Phi", Phi_Bins, Phi_Low, Phi_Up);
600  // fhJetPhi= new TH1F("fhJetPhi", "Jet Phi",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
601  fhJetPhi= new TH1F("fhJetPhi", "Jet Phi",780 , -7, 7);
602  fOutput->Add(fhJetPhi);
603  fhJetEta= new TH1F("fhJetEta", "Jet Eta", Eta_Bins, Eta_Low, Eta_Up);
604  fOutput->Add(fhJetEta);
605  fhJetMass= new TH1F("fhJetMass", "Jet Mass", 4000,-0.5, 39.5);
606  fOutput->Add(fhJetMass);
607  fhJetRadius= new TH1F("fhJetRadius", "Jet Radius", 100, -0.05,0.995);
608  fOutput->Add(fhJetRadius);
609  fhNumberOfJetTracks= new TH1F("fhNumberOfJetTracks", "Number of Tracks within a Jet", 300, -0.5,299.5);
610  fOutput->Add(fhNumberOfJetTracks);
611  fhSubJetRadius= new TH1F("fhSubJetRadius", "SubJet Radius", 100, -0.05,0.995);
612  fOutput->Add(fhSubJetRadius);
613  fhSubJetPt= new TH1F("fhSubJetPt", "SubJet Pt", 1500, -0.5,149.5);
614  fOutput->Add(fhSubJetPt);
615  fhSubJetMass= new TH1F("fhSubJetMass", "Sub Jet Mass", 4000,-0.5, 39.5);
616  fOutput->Add(fhSubJetMass);
617  fhNumberOfSubJetTracks= new TH1F("fhNumberOfSubJetTracks", "Number of Tracks within a Sub Jet", 300, -0.5,299.5);
618  fOutput->Add(fhNumberOfSubJetTracks);
619  fhJetCounter= new TH1F("fhJetCounter", "Jet Counter", 150, -0.5, 149.5);
620  fOutput->Add(fhJetCounter);
621  fhSubJetCounter = new TH1F("fhSubJetCounter", "SubJet Counter",50, -0.5,49.5);
622  fOutput->Add(fhSubJetCounter);
623  fhEventCounter= new TH1F("fhEventCounter", "Event Counter", 15,0.5,15.5);
624  fOutput->Add(fhEventCounter);
625  fhTrackPhi= new TH1F("fhTrackPhi", "fhTrackPhi",780 , -7, 7);
626  fOutput->Add(fhTrackPhi);
627  fhTrackPhi_Cut= new TH1F("fhTrackPhi_Cut", "fhTrackPhi_Cut",780 , -7, 7);
628  fOutput->Add(fhTrackPhi_Cut);
629  }
630  if(fJetShapeType==AliAnalysisTaskSubJetFraction::kSim || fJetShapeType==AliAnalysisTaskSubJetFraction::kGenOnTheFly){
631  fhSubJettiness1_FJ_KT= new TH1D("fhSubJettiness1_FJ_KT","fhSubJettiness1_FJ_KT",400,-2,2);
632  fOutput->Add(fhSubJettiness1_FJ_KT);
633  fhSubJettiness1_FJ_MIN= new TH1D("fhSubJettiness1_FJ_MIN","fhSubJettiness1_FJ_MIN",400,-2,2);
634  fOutput->Add(fhSubJettiness1_FJ_MIN);
635  fhSubJettiness2_FJ_KT= new TH1D("fhSubJettiness2_FJ_KT","fhSubJettiness2_FJ_KT",400,-2,2);
636  fOutput->Add(fhSubJettiness2_FJ_KT);
637  fhSubJettiness2_FJ_MIN= new TH1D("fhSubJettiness2_FJ_MIN","fhSubJettiness2_FJ_MIN",400,-2,2);
638  fOutput->Add(fhSubJettiness2_FJ_MIN);
639  fhSubJettiness1CheckRatio_FJ_AKT = new TH2D("fhSubJettiness1CheckRatio_FJ_AKT","fhSubJettiness1CheckRatio_FJ_AKT",400,-2,2,300,-1,2);
640  fOutput->Add(fhSubJettiness1CheckRatio_FJ_AKT);
641  fhSubJettiness1CheckRatio_FJ_KT= new TH2D("fhSubJettiness1CheckRatio_FJ_KT","fhSubJettiness1CheckRatio_FJ_KT",400,-2,2,300,-1,2);
642  fOutput->Add(fhSubJettiness1CheckRatio_FJ_KT);
643  fhSubJettiness1CheckRatio_FJ_CA= new TH2D("fhSubJettiness1CheckRatio_FJ_CA","fhSubJettiness1CheckRatio_FJ_CA",400,-2,2,300,-1,2);
644  fOutput->Add(fhSubJettiness1CheckRatio_FJ_CA);
645  fhSubJettiness1CheckRatio_FJ_WTA_KT= new TH2D("fhSubJettiness1CheckRatio_FJ_WTA_KT","fhSubJettiness1CheckRatio_FJ_WTA_KT",400,-2,2,300,-1,2);
646  fOutput->Add(fhSubJettiness1CheckRatio_FJ_WTA_KT);
647  fhSubJettiness1CheckRatio_FJ_WTA_CA= new TH2D("fhSubJettiness1CheckRatio_FJ_WTA_CA","fhSubJettiness1CheckRatio_FJ_WTA_CA",400,-2,2,300,-1,2);
648  fOutput->Add(fhSubJettiness1CheckRatio_FJ_WTA_CA);
649  fhSubJettiness1CheckRatio_FJ_OP_AKT= new TH2D("fhSubJettiness1CheckRatio_FJ_OP_AKT","fhSubJettiness1CheckRatio_FJ_OP_AKT",400,-2,2,300,-1,2);
650  fOutput->Add(fhSubJettiness1CheckRatio_FJ_OP_AKT);
651  fhSubJettiness1CheckRatio_FJ_OP_KT= new TH2D("fhSubJettiness1CheckRatio_FJ_OP_KT","fhSubJettiness1CheckRatio_FJ_OP_KT",400,-2,2,300,-1,2);
652  fOutput->Add(fhSubJettiness1CheckRatio_FJ_OP_KT);
653  fhSubJettiness1CheckRatio_FJ_OP_CA= new TH2D("fhSubJettiness1CheckRatio_FJ_OP_CA","fhSubJettiness1CheckRatio_FJ_OP_CA",400,-2,2,300,-1,2);
654  fOutput->Add(fhSubJettiness1CheckRatio_FJ_OP_CA);
655  fhSubJettiness1CheckRatio_FJ_OP_WTA_KT= new TH2D("fhSubJettiness1CheckRatio_FJ_OP_WTA_KT","fhSubJettiness1CheckRatio_FJ_OP_WTA_KT",400,-2,2,300,-1,2);
656  fOutput->Add(fhSubJettiness1CheckRatio_FJ_OP_WTA_KT);
657  fhSubJettiness1CheckRatio_FJ_OP_WTA_CA= new TH2D("fhSubJettiness1CheckRatio_FJ_OP_WTA_CA","fhSubJettiness1CheckRatio_FJ_OP_WTA_CA",400,-2,2,300,-1,2);
658  fOutput->Add(fhSubJettiness1CheckRatio_FJ_OP_WTA_CA);
659  fhSubJettiness1CheckRatio_FJ_MIN= new TH2D("fhSubJettiness1CheckRatio_FJ_MIN","fhSubJettiness1CheckRatio_FJ_MIN",400,-2,2,300,-1,2);
660  fOutput->Add(fhSubJettiness1CheckRatio_FJ_MIN);
661  fhSubJettiness2CheckRatio_FJ_AKT = new TH2D("fhSubJettiness2CheckRatio_FJ_AKT","fhSubJettiness2CheckRatio_FJ_AKT",400,-2,2,300,-1,2);
662  fOutput->Add(fhSubJettiness2CheckRatio_FJ_AKT);
663  fhSubJettiness2CheckRatio_FJ_KT= new TH2D("fhSubJettiness2CheckRatio_FJ_KT","fhSubJettiness2CheckRatio_FJ_KT",400,-2,2,300,-1,2);
664  fOutput->Add(fhSubJettiness2CheckRatio_FJ_KT);
665  fhSubJettiness2CheckRatio_FJ_CA= new TH2D("fhSubJettiness2CheckRatio_FJ_CA","fhSubJettiness2CheckRatio_FJ_CA",400,-2,2,300,-1,2);
666  fOutput->Add(fhSubJettiness2CheckRatio_FJ_CA);
667  fhSubJettiness2CheckRatio_FJ_WTA_KT= new TH2D("fhSubJettiness2CheckRatio_FJ_WTA_KT","fhSubJettiness2CheckRatio_FJ_WTA_KT",400,-2,2,300,-1,2);
668  fOutput->Add(fhSubJettiness2CheckRatio_FJ_WTA_KT);
669  fhSubJettiness2CheckRatio_FJ_WTA_CA= new TH2D("fhSubJettiness2CheckRatio_FJ_WTA_CA","fhSubJettiness2CheckRatio_FJ_WTA_CA",400,-2,2,300,-1,2);
670  fOutput->Add(fhSubJettiness2CheckRatio_FJ_WTA_CA);
671  fhSubJettiness2CheckRatio_FJ_OP_AKT= new TH2D("fhSubJettiness2CheckRatio_FJ_OP_AKT","fhSubJettiness2CheckRatio_FJ_OP_AKT",400,-2,2,300,-1,2);
672  fOutput->Add(fhSubJettiness2CheckRatio_FJ_OP_AKT);
673  fhSubJettiness2CheckRatio_FJ_OP_KT= new TH2D("fhSubJettiness2CheckRatio_FJ_OP_KT","fhSubJettiness2CheckRatio_FJ_OP_KT",400,-2,2,300,-1,2);
674  fOutput->Add(fhSubJettiness2CheckRatio_FJ_OP_KT);
675  fhSubJettiness2CheckRatio_FJ_OP_CA= new TH2D("fhSubJettiness2CheckRatio_FJ_OP_CA","fhSubJettiness2CheckRatio_FJ_OP_CA",400,-2,2,300,-1,2);
676  fOutput->Add(fhSubJettiness2CheckRatio_FJ_OP_CA);
677  fhSubJettiness2CheckRatio_FJ_OP_WTA_KT= new TH2D("fhSubJettiness2CheckRatio_FJ_OP_WTA_KT","fhSubJettiness2CheckRatio_FJ_OP_WTA_KT",400,-2,2,300,-1,2);
678  fOutput->Add(fhSubJettiness2CheckRatio_FJ_OP_WTA_KT);
679  fhSubJettiness2CheckRatio_FJ_OP_WTA_CA= new TH2D("fhSubJettiness2CheckRatio_FJ_OP_WTA_CA","fhSubJettiness2CheckRatio_FJ_OP_WTA_CA",400,-2,2,300,-1,2);
680  fOutput->Add(fhSubJettiness2CheckRatio_FJ_OP_WTA_CA);
681  fhSubJettiness2CheckRatio_FJ_MIN= new TH2D("fhSubJettiness2CheckRatio_FJ_MIN","fhSubJettiness2CheckRatio_FJ_MIN",400,-2,2,300,-1,2);
682  fOutput->Add(fhSubJettiness2CheckRatio_FJ_MIN);
683  fhSubJettiness2to1CheckRatio_FJ_AKT = new TH2D("fhSubJettiness2to1CheckRatio_FJ_AKT","fhSubJettiness2to1CheckRatio_FJ_AKT",400,-2,2,300,-1,2);
684  fOutput->Add(fhSubJettiness2to1CheckRatio_FJ_AKT);
685  fhSubJettiness2to1CheckRatio_FJ_KT= new TH2D("fhSubJettiness2to1CheckRatio_FJ_KT","fhSubJettiness2to1CheckRatio_FJ_KT",400,-2,2,300,-1,2);
686  fOutput->Add(fhSubJettiness2to1CheckRatio_FJ_KT);
687  fhSubJettiness2to1CheckRatio_FJ_CA= new TH2D("fhSubJettiness2to1CheckRatio_FJ_CA","fhSubJettiness2to1CheckRatio_FJ_CA",400,-2,2,300,-1,2);
688  fOutput->Add(fhSubJettiness2to1CheckRatio_FJ_CA);
689  fhSubJettiness2to1CheckRatio_FJ_WTA_KT= new TH2D("fhSubJettiness2to1CheckRatio_FJ_WTA_KT","fhSubJettiness2to1CheckRatio_FJ_WTA_KT",400,-2,2,300,-1,2);
690  fOutput->Add(fhSubJettiness2to1CheckRatio_FJ_WTA_KT);
691  fhSubJettiness2to1CheckRatio_FJ_WTA_CA= new TH2D("fhSubJettiness2to1CheckRatio_FJ_WTA_CA","fhSubJettiness2to1CheckRatio_FJ_WTA_CA",400,-2,2,300,-1,2);
692  fOutput->Add(fhSubJettiness2to1CheckRatio_FJ_WTA_CA);
693  fhSubJettiness2to1CheckRatio_FJ_OP_AKT= new TH2D("fhSubJettiness2to1CheckRatio_FJ_OP_AKT","fhSubJettiness2to1CheckRatio_FJ_OP_AKT",400,-2,2,300,-1,2);
694  fOutput->Add(fhSubJettiness2to1CheckRatio_FJ_OP_AKT);
695  fhSubJettiness2to1CheckRatio_FJ_OP_KT= new TH2D("fhSubJettiness2to1CheckRatio_FJ_OP_KT","fhSubJettiness2to1CheckRatio_FJ_OP_KT",400,-2,2,300,-1,2);
696  fOutput->Add(fhSubJettiness2to1CheckRatio_FJ_OP_KT);
697  fhSubJettiness2to1CheckRatio_FJ_OP_CA= new TH2D("fhSubJettiness2to1CheckRatio_FJ_OP_CA","fhSubJettiness2to1CheckRatio_FJ_OP_CA",400,-2,2,300,-1,2);
698  fOutput->Add(fhSubJettiness2to1CheckRatio_FJ_OP_CA);
699  fhSubJettiness2to1CheckRatio_FJ_OP_WTA_KT= new TH2D("fhSubJettiness2to1CheckRatio_FJ_OP_WTA_KT","fhSubJettiness2to1CheckRatio_FJ_OP_WTA_KT",400,-2,2,300,-1,2);
700  fOutput->Add(fhSubJettiness2to1CheckRatio_FJ_OP_WTA_KT);
701  fhSubJettiness2to1CheckRatio_FJ_OP_WTA_CA= new TH2D("fhSubJettiness2to1CheckRatio_FJ_OP_WTA_CA","fhSubJettiness2to1CheckRatio_FJ_OP_WTA_CA",400,-2,2,300,-1,2);
702  fOutput->Add(fhSubJettiness2to1CheckRatio_FJ_OP_WTA_CA);
703  fhSubJettiness2to1CheckRatio_FJ_MIN= new TH2D("fhSubJettiness2to1CheckRatio_FJ_MIN","fhSubJettiness2to1CheckRatio_FJ_MIN",400,-2,2,300,-1,2);
704  fOutput->Add(fhSubJettiness2to1CheckRatio_FJ_MIN);
705  fhSubJettiness1_FJ_AKT= new TH1D("fhSubJettiness1_FJ_AKT","fhSubJettiness1_FJ_AKT",400,-2,2);
706  fOutput->Add(fhSubJettiness1_FJ_AKT);
707  fhSubJettiness1_FJ_CA= new TH1D("fhSubJettiness1_FJ_CA","fhSubJettiness1_FJ_CA",400,-2,2);
708  fOutput->Add(fhSubJettiness1_FJ_CA);
709  fhSubJettiness1_FJ_WTA_KT= new TH1D("fhSubJettiness1_FJ_WTA_KT","fhSubJettiness1_FJ_WTA_KT",400,-2,2);
710  fOutput->Add(fhSubJettiness1_FJ_WTA_KT);
711  fhSubJettiness1_FJ_WTA_CA= new TH1D("fhSubJettiness1_FJ_WTA_CA","fhSubJettiness1_FJ_WTA_CA",400,-2,2);
712  fOutput->Add(fhSubJettiness1_FJ_WTA_CA);
713  fhSubJettiness1_FJ_OP_AKT= new TH1D("fhSubJettiness1_FJ_OP_AKT","fhSubJettiness1_FJ_OP_AKT",400,-2,2);
714  fOutput->Add(fhSubJettiness1_FJ_OP_AKT);
715  fhSubJettiness1_FJ_OP_KT= new TH1D("fhSubJettiness1_FJ_OP_KT","fhSubJettiness1_FJ_OP_KT",400,-2,2);
716  fOutput->Add(fhSubJettiness1_FJ_OP_KT);
717  fhSubJettiness1_FJ_OP_CA= new TH1D("fhSubJettiness1_FJ_OP_CA","fhSubJettiness1_FJ_OP_CA",400,-2,2);
718  fOutput->Add(fhSubJettiness1_FJ_OP_CA);
719  fhSubJettiness1_FJ_OP_WTA_KT= new TH1D("fhSubJettiness1_FJ_OP_WTA_KT","fhSubJettiness1_FJ_OP_WTA_KT",400,-2,2);
720  fOutput->Add(fhSubJettiness1_FJ_OP_WTA_KT);
721  fhSubJettiness1_FJ_OP_WTA_CA= new TH1D("fhSubJettiness1_FJ_OP_WTA_CA","fhSubJettiness1_FJ_OP_WTA_CA",400,-2,2);
722  fOutput->Add(fhSubJettiness1_FJ_OP_WTA_CA);
723  fhSubJettiness2_FJ_AKT= new TH1D("fhSubJettiness2_FJ_AKT","fhSubJettiness2_FJ_AKT",400,-2,2);
724  fOutput->Add(fhSubJettiness2_FJ_AKT);
725  fhSubJettiness2_FJ_CA= new TH1D("fhSubJettiness2_FJ_CA","fhSubJettiness2_FJ_CA",400,-2,2);
726  fOutput->Add(fhSubJettiness2_FJ_CA);
727  fhSubJettiness2_FJ_WTA_KT= new TH1D("fhSubJettiness2_FJ_WTA_KT","fhSubJettiness2_FJ_WTA_KT",400,-2,2);
728  fOutput->Add(fhSubJettiness2_FJ_WTA_KT);
729  fhSubJettiness2_FJ_WTA_CA= new TH1D("fhSubJettiness2_FJ_WTA_CA","fhSubJettiness2_FJ_WTA_CA",400,-2,2);
730  fOutput->Add(fhSubJettiness2_FJ_WTA_CA);
731  fhSubJettiness2_FJ_OP_AKT= new TH1D("fhSubJettiness2_FJ_OP_AKT","fhSubJettiness2_FJ_OP_AKT",400,-2,2);
732  fOutput->Add(fhSubJettiness2_FJ_OP_AKT);
733  fhSubJettiness2_FJ_OP_KT= new TH1D("fhSubJettiness2_FJ_OP_KT","fhSubJettiness2_FJ_OP_KT",400,-2,2);
734  fOutput->Add(fhSubJettiness2_FJ_OP_KT);
735  fhSubJettiness2_FJ_OP_CA= new TH1D("fhSubJettiness2_FJ_OP_CA","fhSubJettiness2_FJ_OP_CA",400,-2,2);
736  fOutput->Add(fhSubJettiness2_FJ_OP_CA);
737  fhSubJettiness2_FJ_OP_WTA_KT= new TH1D("fhSubJettiness2_FJ_OP_WTA_KT","fhSubJettiness2_FJ_OP_WTA_KT",400,-2,2);
738  fOutput->Add(fhSubJettiness2_FJ_OP_WTA_KT);
739  fhSubJettiness2_FJ_OP_WTA_CA= new TH1D("fhSubJettiness2_FJ_OP_WTA_CA","fhSubJettiness2_FJ_OP_WTA_CA",400,-2,2);
740  fOutput->Add(fhSubJettiness2_FJ_OP_WTA_CA);
741  fhSubJettiness2to1_FJ_KT= new TH1D("fhSubJettiness2to1_FJ_KT","fhSubJettiness2to1_FJ_KT",400,-2,2);
742  fOutput->Add(fhSubJettiness2to1_FJ_KT);
743  fhSubJettiness2to1_FJ_AKT= new TH1D("fhSubJettiness2to1_FJ_AKT","fhSubJettiness2to1_FJ_AKT",400,-2,2);
744  fOutput->Add(fhSubJettiness2to1_FJ_AKT);
745  fhSubJettiness2to1_FJ_CA= new TH1D("fhSubJettiness2to1_FJ_CA","fhSubJettiness2to1_FJ_CA",400,-2,2);
746  fOutput->Add(fhSubJettiness2to1_FJ_CA);
747  fhSubJettiness2to1_FJ_WTA_KT= new TH1D("fhSubJettiness2to1_FJ_WTA_KT","fhSubJettiness2to1_FJ_WTA_KT",400,-2,2);
748  fOutput->Add(fhSubJettiness2to1_FJ_WTA_KT);
749  fhSubJettiness2to1_FJ_WTA_CA= new TH1D("fhSubJettiness2to1_FJ_WTA_CA","fhSubJettiness2to1_FJ_WTA_CA",400,-2,2);
750  fOutput->Add(fhSubJettiness2to1_FJ_WTA_CA);
751  fhSubJettiness2to1_FJ_OP_AKT= new TH1D("fhSubJettiness2to1_FJ_OP_AKT","fhSubJettiness2to1_FJ_OP_AKT",400,-2,2);
752  fOutput->Add(fhSubJettiness2to1_FJ_OP_AKT);
753  fhSubJettiness2to1_FJ_OP_KT= new TH1D("fhSubJettiness2to1_FJ_OP_KT","fhSubJettiness2to1_FJ_OP_KT",400,-2,2);
754  fOutput->Add(fhSubJettiness2to1_FJ_OP_KT);
755  fhSubJettiness2to1_FJ_OP_CA= new TH1D("fhSubJettiness2to1_FJ_OP_CA","fhSubJettiness2to1_FJ_OP_CA",400,-2,2);
756  fOutput->Add(fhSubJettiness2to1_FJ_OP_CA);
757  fhSubJettiness2to1_FJ_OP_WTA_KT= new TH1D("fhSubJettiness2to1_FJ_OP_WTA_KT","fhSubJettiness2to1_FJ_OP_WTA_KT",400,-2,2);
758  fOutput->Add(fhSubJettiness2to1_FJ_OP_WTA_KT);
759  fhSubJettiness2to1_FJ_OP_WTA_CA= new TH1D("fhSubJettiness2to1_FJ_OP_WTA_CA","fhSubJettiness2to1_FJ_OP_WTA_CA",400,-2,2);
760  fOutput->Add(fhSubJettiness2to1_FJ_OP_WTA_CA);
761  fhSubJettiness2to1_FJ_MIN= new TH1D("fhSubJettiness2to1_FJ_MIN","fhSubJettiness2to1_FJ_MIN",400,-2,2);
762  fOutput->Add(fhSubJettiness2to1_FJ_MIN);
763  }
764  if(fJetShapeType==AliAnalysisTaskSubJetFraction::kTrueDet){
765  fhJetPt_1= new TH1F("fhJetPt_1", "Jet Pt Detector Level",1500,-0.5,149.5 );
766  fOutput->Add(fhJetPt_1);
767  fhJetPt_2= new TH1F("fhJetPt_2", "Jet Pt Particle Level",1500,-0.5,149.5 );
768  fOutput->Add(fhJetPt_2);
769  fhJetPhi_1= new TH1F("fhJetPhi_1", "Jet Phi Detector Level",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
770  fOutput->Add(fhJetPhi_1);
771  fhJetPhi_2= new TH1F("fhJetPhi_2", "Jet Phi Particle Level",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
772  fOutput->Add(fhJetPhi_2);
773  fhJetEta_1= new TH1F("fhJetEta_1", "Jet Eta Detector Level", Eta_Bins, Eta_Low, Eta_Up);
774  fOutput->Add(fhJetEta_1);
775  fhJetEta_2= new TH1F("fhJetEta_2", "Jet Eta Particle Level", Eta_Bins, Eta_Low, Eta_Up);
776  fOutput->Add(fhJetEta_2);
777  fhJetMass_1= new TH1F("fhJetMass_1", "Jet Mass Detector Level", 4000,-0.5, 39.5);
778  fOutput->Add(fhJetMass_1);
779  fhJetMass_2= new TH1F("fhJetMass_2", "Jet Mass Particle Level", 4000,-0.5, 39.5);
780  fOutput->Add(fhJetMass_2);
781  fhJetRadius_1= new TH1F("fhJetRadius_1", "Jet Radius Detector Level", 100, -0.05,0.995);
782  fOutput->Add(fhJetRadius_1);
783  fhJetRadius_2= new TH1F("fhJetRadius_2", "Jet Radius Particle Level", 100, -0.05,0.995);
784  fOutput->Add(fhJetRadius_2);
785  fhNumberOfJetTracks_1= new TH1F("fhNumberOfJetTracks_1", "Number of Tracks within a Jet Detector Level", 300, -0.5,299.5);
786  fOutput->Add(fhNumberOfJetTracks_1);
787  fhNumberOfJetTracks_2= new TH1F("fhNumberOfJetTracks_2", "Number of Tracks within a Jet Particle Level", 300, -0.5,299.5);
788  fOutput->Add(fhNumberOfJetTracks_2);
789  fhSubJetRadius_1= new TH1F("fhSubJetRadius_1", "SubJet Radius Detector Level", 100, -0.05,0.995);
790  fOutput->Add(fhSubJetRadius_1);
791  fhSubJetRadius_2= new TH1F("fhSubJetRadius_2", "SubJet Radius Particle Level", 100, -0.05,0.995);
792  fOutput->Add(fhSubJetRadius_2);
793  fhSubJetPt_1= new TH1F("fhSubJetPt_1", "SubJet Pt Detector Level", 1500, -0.5,149.5);
794  fOutput->Add(fhSubJetPt_1);
795  fhSubJetPt_2= new TH1F("fhSubJetPt_2", "SubJet Pt Particle Level", 1500, -0.5,149.5);
796  fOutput->Add(fhSubJetPt_2);
797  fhSubJetMass_1= new TH1F("fhSubJetMass_1", "Sub Jet Mass Detector Level", 4000,-0.5, 39.5);
798  fOutput->Add(fhSubJetMass_1);
799  fhSubJetMass_2= new TH1F("fhSubJetMass_2", "Sub Jet Mass Particle Level", 4000,-0.5, 39.5);
800  fOutput->Add(fhSubJetMass_2);
801  fhNumberOfSubJetTracks_1= new TH1F("fhNumberOfSubJetTracks_1", "Number of Tracks within a Sub Jet Detector Level", 300, -0.5,299.5);
802  fOutput->Add(fhNumberOfSubJetTracks_1);
803  fhNumberOfSubJetTracks_2= new TH1F("fhNumberOfSubJetTracks_2", "Number of Tracks within a Sub Jet Particle Level", 300, -0.5,299.5);
804  fOutput->Add(fhNumberOfSubJetTracks_2);
805  fhJetCounter_1= new TH1F("fhJetCounter_1", "Jet Counter Detector Level", 150, -0.5, 149.5);
806  fOutput->Add(fhJetCounter_1);
807  fhJetCounter_2= new TH1F("fhJetCounter_2", "Jet Counter Particle Level", 150, -0.5, 149.5);
808  fOutput->Add(fhJetCounter_2);
809  fhSubJetCounter_1 = new TH1F("fhSubJetCounter_1", "SubJet Counter Detector Level",50, -0.5,49.5);
810  fOutput->Add(fhSubJetCounter_1);
811  fhSubJetCounter_2 = new TH1F("fhSubJetCounter_2", "SubJet Counter Particle Level",50, -0.5,49.5);
812  fOutput->Add(fhSubJetCounter_2);
813  fh2PtRatio= new TH2F("fhPtRatio", "MC pt for detector level vs particle level jets",1500,-0.5,149.5,1500,-0.5,149.5);
814  fOutput->Add(fh2PtRatio);
815  fhEventCounter_1= new TH1F("fhEventCounter_1", "Event Counter Detector Level", 15,0.5,15.5);
816  fOutput->Add(fhEventCounter_1);
817  fhEventCounter_2= new TH1F("fhEventCounter_2", "Event Counter Particle Level", 15,0.5,15.5);
818  fOutput->Add(fhEventCounter_2);
819  fhTrackPhi= new TH1F("fhTrackPhi", "fhTrackPhi",780 , -7, 7);
820  fOutput->Add(fhTrackPhi);
821  fhTrackPhi_Cut= new TH1F("fhTrackPhi_Cut", "fhTrackPhi_Cut",780 , -7, 7);
822  fOutput->Add(fhTrackPhi_Cut);
823  }
824  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kDetEmbPart){
825  fhEventCounter= new TH1F("fhEventCounter", "Event Counter", 15,0.5,15.5);
826  fOutput->Add(fhEventCounter);
827  }
828 
829  PostData(1,fOutput);
830  PostData(2,fTreeResponseMatrixAxis);
831  PostData(3,fTreeTracks);
832  // delete [] fShapesVarNames;
833 }
834 
835 //________________________________________________________________________
836 Bool_t AliAnalysisTaskSubJetFraction::Run()
837 {
838  // Run analysis code here, if needed. It will be executed before FillHistograms().
839 
840 
841  return kTRUE;
842 }
843 
844 //________________________________________________________________________
845 Bool_t AliAnalysisTaskSubJetFraction::FillHistograms()
846 {
847 
848  //fhEventCounter Key:
849  // 1: Number of events with a Jet Container
850  // 2: Number of Jets without a Jet Container
851  // 3:
852  // 4: Number of Jets found in all events
853  // 5: Number of Jets that were reclustered in all events
854  // 6: Number of SubJets found in all events
855  // 7: Number of events were primary vertext was found
856  // 8: Number of Jets with more than one SubJet
857  // 9:Number of Jets with more than two SubJets
858  // 12:Number of SubJetinessEvents in kData
859 
860 
861  // const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
862  // Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
863  // if(vert) fhEventCounter->Fill(7);
864 
865  //cout << ((AliVAODHeader*)(InputEvent()->GetHeader()))->GetEventplane()<< " "<<fEPV0<<endl;
866  // cout << InputEvent()->GetEventplane()->GetEventplane("Q")<< " "<<fEPV0<<endl;
867  if (fCentSelectOn){
868  if ((fCent>fCentMax) || (fCent<fCentMin)) return 0;
869  }
870 
872  AliTrackContainer *PartCont_Particles = NULL;
873  AliParticleContainer *PartContMC_Particles = NULL;
874 
875  if (fJetShapeSub==kConstSub){
876  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) PartContMC_Particles = GetParticleContainer(1);
877  else PartCont_Particles = GetTrackContainer(1);
878  }
879  else{
880  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) PartContMC_Particles = GetParticleContainer(0);
881  else PartCont_Particles = GetTrackContainer(0);
882  }
883 
884  TClonesArray *TracksArray_Particles = NULL;
885  TClonesArray *TracksArrayMC_Particles = NULL;
886  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly && PartContMC_Particles) TracksArrayMC_Particles = PartContMC_Particles->GetArray();
887  else if(PartCont_Particles) TracksArray_Particles = PartCont_Particles->GetArray();
888 
889  AliAODTrack *Track_Particles = 0x0;
890  Int_t NTracks_Particles=0;
891  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly && TracksArrayMC_Particles) NTracks_Particles = TracksArrayMC_Particles->GetEntriesFast();
892  else if(TracksArray_Particles) NTracks_Particles = TracksArray_Particles->GetEntriesFast();
893 
894  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly){
895  if (PartContMC_Particles && TracksArrayMC_Particles){
896  for(Int_t i=0; i < NTracks_Particles; i++){
897  if((Track_Particles = static_cast<AliAODTrack*>(PartContMC_Particles->GetAcceptParticle(i)))){
898  if (!Track_Particles) continue;
899  if(TMath::Abs(Track_Particles->Eta())>0.9) continue;
900  if (Track_Particles->Pt()<0.15) continue;
901  fhTrackPhi->Fill(Track_Particles->Phi());
902  if (Track_Particles->Pt()>=4.0) fhTrackPhi_Cut->Fill(Track_Particles->Phi());
903  }
904  }
905  }
906  }
907  else{
908  if (PartCont_Particles && TracksArray_Particles){
909  for(Int_t i=0; i < NTracks_Particles; i++){
910  if((Track_Particles = static_cast<AliAODTrack*>(PartCont_Particles->GetAcceptTrack(i)))){
911  if (!Track_Particles) continue;
912  if(TMath::Abs(Track_Particles->Eta())>0.9) continue;
913  if (Track_Particles->Pt()<0.15) continue;
914  fhTrackPhi->Fill(Track_Particles->Phi());
915  if (Track_Particles->Pt()>=4.0) fhTrackPhi_Cut->Fill(Track_Particles->Phi());
916  }
917  }
918  }
919  }
921 
922  AliAODTrack *TriggerHadron = 0x0;
923  if (fJetSelection == kRecoil) {
924  //you have to set a flag and the limits of the pT interval for your trigger
925  Int_t TriggerHadronLabel = SelectTriggerHadron(fPtMinTriggerHadron, fPtMaxTriggerHadron);
926  if (TriggerHadronLabel==-99999) return 0; //Trigger Hadron Not Found
927  AliTrackContainer *PartCont =NULL;
928  AliParticleContainer *PartContMC=NULL;
929  if (fJetShapeSub==kConstSub){
930  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) PartContMC = GetParticleContainer(1);
931  else PartCont = GetTrackContainer(1);
932  }
933  else{
934  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) PartContMC = GetParticleContainer(0);
935  else PartCont = GetTrackContainer(0);
936  }
937  TClonesArray *TrackArray = NULL;
938  TClonesArray *TrackArrayMC = NULL;
939  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) TrackArrayMC = PartContMC->GetArray();
940  else TrackArray = PartCont->GetArray();
941  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) TriggerHadron = static_cast<AliAODTrack*>(TrackArrayMC->At(TriggerHadronLabel));
942  else TriggerHadron = static_cast<AliAODTrack*>(TrackArray->At(TriggerHadronLabel));
943  if (!TriggerHadron) return 0;//No trigger hadron with label found
944  if(fSemigoodCorrect){
945  Double_t HoleDistance=RelativePhi(TriggerHadron->Phi(),fHolePos);
946  if(TMath::Abs(HoleDistance)+fHoleWidth+fJetRadius>TMath::Pi()-fRecoilAngularWindow) return 0;
947  }
948  fhPtTriggerHadron->Fill(TriggerHadron->Pt()); //Needed for per trigger Normalisation
949  if (fJetShapeType != AliAnalysisTaskSubJetFraction::kGenOnTheFly) fhPhiTriggerHadronEventPlane->Fill(TMath::Abs(RelativePhiEventPlane(fEPV0,TriggerHadron->Phi()))); //fEPV0 is the event plane from AliAnalysisTaskEmcal
950  if (fJetShapeType != AliAnalysisTaskSubJetFraction::kGenOnTheFly) fhPhiTriggerHadronEventPlaneTPC->Fill(TMath::Abs(RelativePhiEventPlane(((AliVAODHeader*)(InputEvent()->GetHeader()))->GetEventplane(),TriggerHadron->Phi()))); //TPC event plane
951  }
952 
953 
955  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kDetEmbPart){
956  AliEmcalJet *Jet1 = NULL; //Subtracted hybrid Jet
957  AliEmcalJet *Jet2 = NULL; //Unsubtracted Hybrid Jet
958  AliEmcalJet *Jet3 = NULL; //Detector Level Pythia Jet
959  AliEmcalJet *Jet4 = NULL; //Particle Level Pyhtia Jet
960  AliJetContainer *JetCont1= GetJetContainer(0); //Jet Container for Subtracted Hybrid Jets
961  AliJetContainer *JetCont2= GetJetContainer(1); //Jet Container for Unsubtracted Hybrid Jets
962  AliJetContainer *JetCont3= GetJetContainer(2); //Jet Container for Detector Level Pyhtia Jets
963  AliJetContainer *JetCont4= GetJetContainer(3); //Jet Container for Particle Level Pythia Jets
964  AliEmcalJetFinder *Reclusterer1; //Object containg Subjets from Subtracted Hybrid Jets
965  AliEmcalJetFinder *Reclusterer4; //Object containg Subjets from Particle Level Pythia Jets
966  Bool_t JetsMatched=kFALSE;
967  Bool_t EventCounter=kFALSE;
968  Int_t JetNumber=-1;
969  Double_t JetPtThreshold=-2;
970  fhEventCounter->Fill(1);
971  if(JetCont1) {
972  fhEventCounter->Fill(2);
973  JetCont1->ResetCurrentID();
974  while((Jet1=JetCont1->GetNextAcceptJet())) {
975  if (fJetShapeSub==kConstSub) JetPtThreshold=Jet1->Pt();
976  else JetPtThreshold=Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
977  if ( (!Jet1) || (JetPtThreshold<fPtThreshold)) continue;
978  else {
979  /* if(fSemigoodCorrect){
980  Double_t HoleDistance=RelativePhi(Jet1->Phi(),fHolePos);
981  if(TMath::Abs(HoleDistance)<fHoleWidth) continue;
982  }*/
983  Float_t RecoilDeltaPhi = 0.;
984  if (fJetSelection == kRecoil){
985  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), Jet1->Phi());
986  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
987  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), Jet1->Pt());
988  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), Jet1->Phi()));
989  }
990  if (!EventCounter){
991  fhEventCounter->Fill(3);
992  EventCounter=kTRUE;
993  }
994  if(fJetShapeSub==kConstSub){
995  JetNumber=-1;
996  for(Int_t i = 0; i<JetCont2->GetNJets(); i++) {
997  Jet2 = JetCont2->GetJet(i);
998  if(Jet2->GetLabel()==Jet1->GetLabel()) {
999  JetNumber=i;
1000  }
1001  }
1002  if(JetNumber==-1) continue;
1003  Jet2=JetCont2->GetJet(JetNumber);
1004  if (JetCont2->AliJetContainer::GetFractionSharedPt(Jet2)<fSharedFractionPtMin) continue;
1005  Jet3=Jet2->ClosestJet();
1006  }
1007  if(!(fJetShapeSub==kConstSub)){
1008  if (JetCont1->AliJetContainer::GetFractionSharedPt(Jet1)<fSharedFractionPtMin) continue;
1009  Jet3 = Jet1->ClosestJet(); //Note for NoSub and Deriv Sub cases you must fill both the Unsubtracted and Subtracted Hybrid jet containers with the same jet branch
1010  }
1011  if (!Jet3) continue;
1012  Jet4=Jet3->ClosestJet();
1013  if(!Jet4) continue;
1014  JetsMatched=kTRUE;
1015 
1016  std::vector <Double_t> Jet_Constituents_Emb;
1017  std::vector <Double_t> Jet_Constituents_Truth;
1018  AliVParticle *JetConstituent_Emb=NULL;
1019  AliVParticle *JetConstituent_Truth=NULL;
1020  if (fMLOn==1){
1021  for (Int_t Tracks_Emb=0; Tracks_Emb<Jet1->GetNumberOfTracks(); Tracks_Emb++){
1022  JetConstituent_Emb = static_cast<AliVParticle*>(Jet1->TrackAt(Tracks_Emb, JetCont1->GetParticleContainer()->GetArray()));
1023  if(!JetConstituent_Emb) continue;
1024  Jet_Constituents_Emb.push_back(JetConstituent_Emb->Px());
1025  Jet_Constituents_Emb.push_back(JetConstituent_Emb->Py());
1026  Jet_Constituents_Emb.push_back(JetConstituent_Emb->Pz());
1027  Jet_Constituents_Emb.push_back(JetConstituent_Emb->E());
1028  }
1029  for (Int_t Tracks_Truth=0; Tracks_Truth<Jet4->GetNumberOfTracks(); Tracks_Truth++){
1030  JetConstituent_Truth = static_cast<AliVParticle*>(Jet4->TrackAt(Tracks_Truth, JetCont4->GetParticleContainer()->GetArray()));
1031  if(!JetConstituent_Truth) continue;
1032  Jet_Constituents_Truth.push_back(JetConstituent_Truth->Px());
1033  Jet_Constituents_Truth.push_back(JetConstituent_Truth->Py());
1034  Jet_Constituents_Truth.push_back(JetConstituent_Truth->Pz());
1035  Jet_Constituents_Truth.push_back(JetConstituent_Truth->E());
1036  }
1037  }
1038 
1039  if (fMLOn==1) fShapesVar_Tracks_Rec.push_back(Jet_Constituents_Emb);
1040  if (fMLOn==1)fShapesVar_Tracks_Truth.push_back(Jet_Constituents_Truth);
1041 
1043  if (fJetShapeSub==kConstSub) fShapesVar[0]=Jet1->Pt();
1044  else fShapesVar[0]=Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
1045  if (fMLOn==0) fShapesVar[2]=FjNSubJettiness(Jet1,0,1,0,1,0);
1046  else fShapesVar[2]=Jet1->Eta();
1047  if (fMLOn==0) fShapesVar[4]=FjNSubJettiness(Jet1,0,2,0,1,0);
1048  else fShapesVar[4]=Jet1->Phi();
1049  if (fMLOn==0) fShapesVar[6]=FjNSubJettiness(Jet1,0,2,0,1,8);
1050  else fShapesVar[6]=Jet1->M();
1051  if (fMLOn==0) fShapesVar[8]=FjNSubJettiness(Jet1,0,2,0,1,1);
1052  else fShapesVar[8]=Jet1->GetNumberOfTracks();
1053  if (fMLOn==0) fShapesVar[10]=FjNSubJettiness(Jet1,0,2,0,1,9);
1054  else fShapesVar[10]=0.0;
1055  if (fMLOn==0) fShapesVar[12]=FjNSubJettiness(Jet1,0,2,0,1,3,fBeta_SD,fZCut);
1056  else fShapesVar[12]=0.0;
1057  if (fMLOn==0) fShapesVar[14]=FjNSubJettiness(Jet1,0,2,0,1,5,fBeta_SD,fZCut);
1058  else fShapesVar[14]=0.0;
1059  if (fMLOn==0) fShapesVar[16]=Jet1->GetLeadingTrack(JetCont1->GetParticleContainer()->GetArray())->Pt();
1060  else fShapesVar[16]=0.0;
1061  if (fMLOn==0) fShapesVar[18]=RelativePhiEventPlane(fEPV0,Jet1->Phi());
1062  else fShapesVar[18]=0.0;
1063  //fShapesVar[20]=RelativePhiEventPlane(((AliVAODHeader*)(InputEvent()->GetHeader()))->GetEventplane(),Jet1->Phi());
1064  if (fMLOn==0) fShapesVar[20]=FjNSubJettiness(Jet1,0,2,0,1,6,fBeta_SD,fZCut);
1065  else fShapesVar[20]=0.0;
1066  if (fFullTree){
1067  if (fMLOn==0) fShapesVar[22]=FjNSubJettiness(Jet1,0,2,0,1,2);
1068  else fShapesVar[22]=0.0;
1069  Reclusterer1 = Recluster(Jet1, 0, fSubJetRadius, fSubJetMinPt, fSubJetAlgorithm, "SubJetFinder_1");
1070  if (fMLOn==0) fShapesVar[24]=SubJetFraction(Jet1, Reclusterer1, 1, 0, kTRUE, kFALSE);
1071  else fShapesVar[24]=0.0;
1072  if (fMLOn==0) fShapesVar[26]=SubJetFraction(Jet1, Reclusterer1, 2, 0, kTRUE, kFALSE);
1073  else fShapesVar[26]=0.0;
1074  }
1075  if (JetsMatched){ //even needed? Not now but might be if you want to fill trees when jets aren't matched too
1076  fShapesVar[1]=Jet4->Pt();
1077  if (fMLOn==0) fShapesVar[3]=FjNSubJettiness(Jet4,3,1,0,1,0);
1078  else fShapesVar[3]=Jet4->Eta();
1079  if (fMLOn==0) fShapesVar[5]=FjNSubJettiness(Jet4,3,2,0,1,0);
1080  else fShapesVar[5]=Jet4->Phi();
1081  if (fMLOn==0) fShapesVar[7]=FjNSubJettiness(Jet4,3,2,0,1,8);
1082  else fShapesVar[7]=Jet4->M();
1083  if (fMLOn==0) fShapesVar[9]=FjNSubJettiness(Jet4,3,2,0,1,1);
1084  else fShapesVar[9]=Jet4->GetNumberOfTracks();
1085  if (fMLOn==0) fShapesVar[11]=FjNSubJettiness(Jet4,3,2,0,1,9);
1086  else fShapesVar[11]=0.0;
1087  if (fMLOn==0) fShapesVar[13]=FjNSubJettiness(Jet4,3,2,0,1,3,fBeta_SD,fZCut);
1088  else fShapesVar[13]=0.0;
1089  if (fMLOn==0) fShapesVar[15]=FjNSubJettiness(Jet4,3,2,0,1,5,fBeta_SD,fZCut);
1090  else fShapesVar[15]=0.0;
1091  if (fMLOn==0) fShapesVar[17]=Jet4->GetLeadingTrack(JetCont4->GetParticleContainer()->GetArray())->Pt();
1092  else fShapesVar[17]=0.0;
1093  if (fMLOn==0) fShapesVar[19]=RelativePhiEventPlane(fEPV0,Jet4->Phi());
1094  else fShapesVar[19]=0.0;
1095  //fShapesVar[21]=RelativePhiEventPlane(((AliVAODHeader*)(InputEvent()->GetHeader()))->GetEventplane(),Jet4->Phi());
1096  if (fMLOn==0) fShapesVar[21]=FjNSubJettiness(Jet4,3,2,0,1,6,fBeta_SD,fZCut);
1097  else fShapesVar[21]=0.0;
1098  if (fFullTree){
1099  if (fMLOn==0) fShapesVar[23]=FjNSubJettiness(Jet4,3,2,0,1,2);
1100  else fShapesVar[23]=0.0;
1101  Reclusterer4=Recluster(Jet4, 3, fSubJetRadius, 0, fSubJetAlgorithm, "SubJetFinder_4");
1102  if (fMLOn==0) fShapesVar[25]=SubJetFraction(Jet4, Reclusterer4, 1, 0, kTRUE, kFALSE);
1103  else fShapesVar[25]=0.0;
1104  if (fMLOn==0) fShapesVar[27]=SubJetFraction(Jet4, Reclusterer4, 2, 0, kTRUE, kFALSE);
1105  else fShapesVar[27]=0.0;
1106  }
1107  }
1108  else{
1109  fShapesVar[1]=-2;
1110  fShapesVar[3]=-2;
1111  fShapesVar[5]=-2;
1112  fShapesVar[7]=-2;
1113  fShapesVar[9]=-2;
1114  fShapesVar[11]=-2;
1115  fShapesVar[13]=-2;
1116  fShapesVar[15]=-2;
1117  fShapesVar[17]=-2;
1118  fShapesVar[19]=-2;
1119  fShapesVar[21]=-2;
1120  if (fFullTree){
1121  fShapesVar[23]=-2;
1122  fShapesVar[25]=-2;
1123  fShapesVar[27]=-2;
1124  }
1125  }
1126  fTreeResponseMatrixAxis->Fill();
1127  if (fMLOn==1) fTreeTracks->Fill();
1128  JetsMatched=kFALSE;
1129  }
1130  }
1131  }
1132  }
1133 
1135  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kTrueDet){
1136  AliEmcalJet *Jet1 = NULL; //Detector Level Jet
1137  AliEmcalJet *Jet2 = NULL; //Particle Level Jet
1138  AliJetContainer *JetCont1= GetJetContainer(0); //Jet Container for Detector Level Pythia
1139  AliJetContainer *JetCont2= GetJetContainer(1); //Jet Container for Particle Level Pythia
1140  AliEmcalJetFinder *Reclusterer1; //Object containg Subjets Detector Level
1141  AliEmcalJetFinder *Reclusterer2; //Object containg Subjets Particle Level
1142  Int_t JetCounter1=0; //Counts number of jets in event
1143  Int_t JetCounter2=0; //Counts number of jets in event
1144  Double_t JetPhi1=0;
1145  Double_t JetPhi2=0;
1146  Bool_t JetsMatched=kFALSE;
1147  Double_t Pythia_Event_Weight=1;
1148  Bool_t EventCounter=kFALSE;
1149  fhEventCounter_1->Fill(1);
1150  if(JetCont1) {
1151  fhEventCounter_1->Fill(2); //Number of events with a jet container
1152  JetCont1->ResetCurrentID();
1153  while((Jet1=JetCont1->GetNextAcceptJet())) {
1154  std::vector <Double_t> Jet_Constituents_Det;
1155  std::vector <Double_t> Jet_Constituents_Truth;
1156  if( (!Jet1) || ((Jet1->Pt())<fPtThreshold)) {
1157  // fhEventCounter_1->Fill(3); //events where the jet had a problem
1158  continue;
1159  }
1160  else {
1161  /* if(fSemigoodCorrect){
1162  Double_t HoleDistance=RelativePhi(Jet1->Phi(),fHolePos);
1163  if(TMath::Abs(HoleDistance)<fHoleWidth) continue;
1164  }*/
1165  Float_t RecoilDeltaPhi = 0.;
1166  if (fJetSelection == kRecoil){
1167  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), Jet1->Phi());
1168  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
1169  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), Jet1->Pt());
1170  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), Jet1->Phi()));
1171  }
1172  if (!EventCounter){
1173  fhEventCounter_1->Fill(3);
1174  EventCounter=kTRUE;
1175  }
1176  if((Jet1->GetNumberOfTracks())==0){
1177  fhEventCounter_1->Fill(10); //zero track jets
1178  }
1179  if((Jet1->GetNumberOfTracks())==1){
1180  fhEventCounter_1->Fill(11); //one track jets
1181  }
1182  fhEventCounter_1->Fill(4); //Number of Jets found in all events
1183  JetCounter1++;
1184  fhJetPt_1->Fill(Jet1->Pt());
1185  JetPhi1=Jet1->Phi();
1186  if(JetPhi1 < -1*TMath::Pi()) JetPhi1 += (2*TMath::Pi());
1187  else if (JetPhi1 > TMath::Pi()) JetPhi1 -= (2*TMath::Pi());
1188  fhJetPhi_1->Fill(JetPhi1);
1189  fhJetEta_1->Fill(Jet1->Eta());
1190  fhJetMass_1->Fill(Jet1->M());
1191  fhJetRadius_1->Fill(TMath::Sqrt((Jet1->Area()/TMath::Pi()))); //Radius of Jets per event
1192  fhNumberOfJetTracks_1->Fill(Jet1->GetNumberOfTracks());
1193  if((Jet2 = Jet1->ClosestJet())){
1194  JetsMatched=kTRUE;
1195  if((Jet2->GetNumberOfTracks())==0){
1196  fhEventCounter_2->Fill(10); //zero track jets
1197  }
1198  if((Jet2->GetNumberOfTracks())==1){
1199  fhEventCounter_2->Fill(11); //one track jets
1200  }
1201  fhEventCounter_2->Fill(4); //Number of Jets found in all events
1202  JetCounter2++;
1203  fhJetPt_2->Fill(Jet2->Pt());
1204  JetPhi2=Jet2->Phi();
1205  if(JetPhi2 < -1*TMath::Pi()) JetPhi2 += (2*TMath::Pi());
1206  else if (JetPhi2 > TMath::Pi()) JetPhi2 -= (2*TMath::Pi());
1207  fhJetPhi_2->Fill(JetPhi2);
1208  fhJetEta_2->Fill(Jet2->Eta());
1209  fhJetMass_2->Fill(Jet2->M());
1210  fhJetRadius_2->Fill(TMath::Sqrt((Jet2->Area()/TMath::Pi()))); //Radius of Jets per event
1211  fhNumberOfJetTracks_2->Fill(Jet2->GetNumberOfTracks());
1212  fh2PtRatio->Fill(Jet1->Pt(),Jet2->Pt(),Pythia_Event_Weight);
1213 
1214 
1215  AliVParticle *JetConstituent_Det=NULL;
1216  AliVParticle *JetConstituent_Truth=NULL;
1217  if (fMLOn==1){
1218  for (Int_t Tracks_Det=0; Tracks_Det<Jet1->GetNumberOfTracks(); Tracks_Det++){
1219  JetConstituent_Det = static_cast<AliVParticle*>(Jet1->TrackAt(Tracks_Det, JetCont1->GetParticleContainer()->GetArray()));
1220  if(!JetConstituent_Det) continue;
1221  Jet_Constituents_Det.push_back(JetConstituent_Det->Px());
1222  Jet_Constituents_Det.push_back(JetConstituent_Det->Py());
1223  Jet_Constituents_Det.push_back(JetConstituent_Det->Pz());
1224  Jet_Constituents_Det.push_back(JetConstituent_Det->E());
1225  }
1226  for (Int_t Tracks_Truth=0; Tracks_Truth<Jet2->GetNumberOfTracks(); Tracks_Truth++){
1227  JetConstituent_Truth = static_cast<AliVParticle*>(Jet2->TrackAt(Tracks_Truth, JetCont2->GetParticleContainer()->GetArray()));
1228  if(!JetConstituent_Truth) continue;
1229  Jet_Constituents_Truth.push_back(JetConstituent_Truth->Px());
1230  Jet_Constituents_Truth.push_back(JetConstituent_Truth->Py());
1231  Jet_Constituents_Truth.push_back(JetConstituent_Truth->Pz());
1232  Jet_Constituents_Truth.push_back(JetConstituent_Truth->E());
1233  }
1234  }
1235 
1236  }
1237  else {
1238  fhEventCounter_2->Fill(1);
1239  //continue;
1240  }
1242 
1243  if (fMLOn==1) fShapesVar_Tracks_Rec.push_back(Jet_Constituents_Det);
1244  if (fMLOn==1)fShapesVar_Tracks_Truth.push_back(Jet_Constituents_Truth);
1245 
1246  fShapesVar[0]=Jet1->Pt();
1247  if (fMLOn==0) fShapesVar[2]=FjNSubJettiness(Jet1,0,1,0,1,0);
1248  else fShapesVar[2]=Jet1->Eta();
1249  if (fMLOn==0) fShapesVar[4]=FjNSubJettiness(Jet1,0,2,0,1,0);
1250  else fShapesVar[4]=Jet1->Phi();
1251  if (fMLOn==0) fShapesVar[6]=FjNSubJettiness(Jet1,0,2,0,1,8);
1252  else fShapesVar[6]=Jet1->M();
1253  if (fMLOn==0) fShapesVar[8]=FjNSubJettiness(Jet1,0,2,0,1,1);
1254  else fShapesVar[8]=Jet1->GetNumberOfTracks();
1255  if (fMLOn==0) fShapesVar[10]=FjNSubJettiness(Jet1,0,2,0,1,9);
1256  else fShapesVar[10]=0.0;
1257  if (fMLOn==0) fShapesVar[12]=FjNSubJettiness(Jet1,0,2,0,1,3,fBeta_SD,fZCut);
1258  else fShapesVar[12]=0.0;
1259  if (fMLOn==0) fShapesVar[14]=FjNSubJettiness(Jet1,0,2,0,1,5,fBeta_SD,fZCut);
1260  else fShapesVar[14]=0.0;
1261  if (fMLOn==0) fShapesVar[16]=Jet1->GetLeadingTrack(JetCont1->GetParticleContainer()->GetArray())->Pt();
1262  else fShapesVar[16]=0.0;
1263  if (fMLOn==0) fShapesVar[18]=-2; //event plane calculation only needed for PbPb recoils
1264  else fShapesVar[18]=0.0;
1265  if (fMLOn==0) fShapesVar[20]=FjNSubJettiness(Jet1,0,2,0,1,6,fBeta_SD,fZCut);
1266  else fShapesVar[20]=0.0;
1267  Reclusterer1 = Recluster(Jet1, 0, fSubJetRadius, fSubJetMinPt, fSubJetAlgorithm, "SubJetFinder_1");
1268  if (fFullTree){
1269  if (fMLOn==0) fShapesVar[22]=FjNSubJettiness(Jet1,0,2,0,1,2);
1270  else fShapesVar[22]=0.0;
1271  if (fMLOn==0) fShapesVar[24]=SubJetFraction(Jet1, Reclusterer1, 1, 0, kTRUE, kFALSE);
1272  else fShapesVar[24]=0.0;
1273  if (fMLOn==0) fShapesVar[26]=SubJetFraction(Jet1, Reclusterer1, 2, 0, kTRUE, kFALSE);
1274  else fShapesVar[26]=0.0;
1275  }
1276  if (JetsMatched){ //even needed? Not now but might be if you want to fill trees when jets aren't matched too
1277  fShapesVar[1]=Jet2->Pt();
1278  if (fMLOn==0) fShapesVar[3]=FjNSubJettiness(Jet2,1,1,0,1,0);
1279  else fShapesVar[3]=Jet2->Eta();
1280  if (fMLOn==0) fShapesVar[5]=FjNSubJettiness(Jet2,1,2,0,1,0);
1281  else fShapesVar[5]=Jet2->Phi();
1282  if (fMLOn==0) fShapesVar[7]=FjNSubJettiness(Jet2,1,2,0,1,8);
1283  else fShapesVar[7]=Jet2->M();
1284  if (fMLOn==0) fShapesVar[9]=FjNSubJettiness(Jet2,1,2,0,1,1);
1285  else fShapesVar[9]=Jet2->GetNumberOfTracks();
1286  if (fMLOn==0) fShapesVar[11]=FjNSubJettiness(Jet2,1,2,0,1,9);
1287  else fShapesVar[11]=0.0;
1288  if (fMLOn==0) fShapesVar[13]=FjNSubJettiness(Jet2,1,2,0,1,3,fBeta_SD,fZCut);
1289  else fShapesVar[13]=0.0;
1290  if (fMLOn==0) fShapesVar[15]=FjNSubJettiness(Jet2,1,2,0,1,5,fBeta_SD,fZCut);
1291  else fShapesVar[15]=0.0;
1292  if (fMLOn==0) fShapesVar[17]=Jet2->GetLeadingTrack(JetCont2->GetParticleContainer()->GetArray())->Pt();
1293  else fShapesVar[17]=0.0;
1294  if (fMLOn==0) fShapesVar[19]=-2;
1295  else fShapesVar[19]=0.0;
1296  if (fMLOn==0) fShapesVar[21]=FjNSubJettiness(Jet2,1,2,0,1,6,fBeta_SD,fZCut);
1297  else fShapesVar[21]=0.0;
1298  Reclusterer2 = Recluster(Jet2, 1, fSubJetRadius, 0, fSubJetAlgorithm, "SubJetFinder_2");
1299  if (fFullTree){
1300  if (fMLOn==0) fShapesVar[23]=FjNSubJettiness(Jet2,1,2,0,1,2);
1301  else fShapesVar[23]=0.0;
1302  if (fMLOn==0) fShapesVar[25]=SubJetFraction(Jet2, Reclusterer2, 1, 0, kTRUE, kFALSE);
1303  else fShapesVar[25]=0.0;
1304  if (fMLOn==0) fShapesVar[27]=SubJetFraction(Jet2, Reclusterer2, 2, 0, kTRUE, kFALSE);
1305  else fShapesVar[27]=0.0;
1306  }
1307  }
1308  else{
1309  fShapesVar[1]=-2;
1310  fShapesVar[3]=-2;
1311  fShapesVar[5]=-2;
1312  fShapesVar[7]=-2;
1313  fShapesVar[9]=-2;
1314  fShapesVar[11]=-2;
1315  fShapesVar[13]=-2;
1316  fShapesVar[15]=-2;
1317  fShapesVar[17]=-2;
1318  fShapesVar[19]=-2;
1319  fShapesVar[21]=-2;
1320  if (fFullTree){
1321  fShapesVar[23]=-2;
1322  fShapesVar[25]=-2;
1323  fShapesVar[27]=-2;
1324  }
1325  }
1326  fTreeResponseMatrixAxis->Fill();
1327  if (fMLOn==1)fTreeTracks->Fill();
1328 
1329  fhSubJetCounter_1->Fill(Reclusterer1->GetNumberOfJets());
1330  for (Int_t i= 0; i<Reclusterer1->GetNumberOfJets(); i++){
1331  fhEventCounter_1->Fill(6); //Number of overall subjets in all events
1332  fhSubJetPt_1->Fill(Reclusterer1->GetJet(i)->Pt());
1333  fhSubJetMass_1->Fill(Reclusterer1->GetJet(i)->M());
1334  fhNumberOfSubJetTracks_1->Fill(Reclusterer1->GetJet(i)->GetNumberOfTracks());
1335  fhSubJetRadius_1->Fill(TMath::Sqrt((Reclusterer1->GetJet(i)->Area()/TMath::Pi()))); //Radius of SubJets per event
1336  }
1337  if(JetsMatched){
1338  fhSubJetCounter_2->Fill(Reclusterer2->GetNumberOfJets());
1339  for (Int_t i= 0; i<Reclusterer2->GetNumberOfJets(); i++){
1340  fhEventCounter_2->Fill(6); //Number of overall subjets in all events
1341  fhSubJetPt_2->Fill(Reclusterer2->GetJet(i)->Pt());
1342  fhSubJetMass_2->Fill(Reclusterer2->GetJet(i)->M());
1343  fhNumberOfSubJetTracks_2->Fill(Reclusterer2->GetJet(i)->GetNumberOfTracks());
1344  fhSubJetRadius_2->Fill(TMath::Sqrt((Reclusterer2->GetJet(i)->Area()/TMath::Pi()))); //Radius of SubJets per event
1345  }
1346  }
1347  JetsMatched=kFALSE;
1348  }
1349  }
1350  fhJetCounter_1->Fill(JetCounter1); //Number of Jets in Each Event Particle Level
1351  fhJetCounter_2->Fill(JetCounter2); //Number of Jets in Each Event Detector Level
1352  }
1353  //else {fhEventCounter_1->Fill(3);} //Events with no jet container
1354  }
1355 
1356 
1357 
1358  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kData){
1359  AliEmcalJet *Jet1 = NULL; //Original Jet in the event
1360  AliJetContainer *JetCont= GetJetContainer(0); //Jet Container for event
1361  Int_t JetCounter=0; //Counts number of jets in event
1362  Double_t JetPhi=0;
1363  Bool_t EventCounter=kFALSE;
1364  Double_t JetPt_ForThreshold=0;
1365  fhEventCounter->Fill(1);
1366  if(JetCont) {
1367  fhEventCounter->Fill(2); //Number of events with a jet container
1368  JetCont->ResetCurrentID();
1369  while((Jet1=JetCont->GetNextAcceptJet())) {
1370  if(!Jet1) continue;
1371  if (fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) JetPt_ForThreshold = Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
1372  else JetPt_ForThreshold = Jet1->Pt();
1373  if(JetPt_ForThreshold<fPtThreshold) {
1374  //fhEventCounter->Fill(3); //events where the jet had a problem
1375  continue;
1376  }
1377  else {
1378  /* if(fSemigoodCorrect){
1379  Double_t HoleDistance=RelativePhi(Jet1->Phi(),fHolePos);
1380  if(TMath::Abs(HoleDistance)<fHoleWidth) continue;
1381  }*/
1382  Float_t RecoilDeltaPhi = 0.;
1383  if (fJetSelection == kRecoil){
1384  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), Jet1->Phi());
1385  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
1386  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), Jet1->Pt());
1387  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), Jet1->Phi()));
1388  }
1389  if (!EventCounter){
1390  fhEventCounter->Fill(3);
1391  EventCounter=kTRUE;
1392  }
1393  if((Jet1->GetNumberOfTracks())==0){
1394  fhEventCounter->Fill(10); //zero track jets
1395  }
1396  if((Jet1->GetNumberOfTracks())==1){
1397  fhEventCounter->Fill(11); //one track jets
1398  }
1399  fhEventCounter->Fill(4); //Number of Jets found in all events
1400  JetCounter++;
1401  fhJetPt->Fill(Jet1->Pt());
1402  JetPhi=Jet1->Phi();
1403  //if(JetPhi < -1*TMath::Pi()) JetPhi += (2*TMath::Pi());
1404  //else if (JetPhi > TMath::Pi()) JetPhi -= (2*TMath::Pi());
1405  fhJetPhi->Fill(JetPhi);
1406  fhJetEta->Fill(Jet1->Eta());
1407  fhJetMass->Fill(Jet1->M());
1408  fhJetRadius->Fill(TMath::Sqrt((Jet1->Area()/TMath::Pi()))); //Radius of Jets per event
1409  fhNumberOfJetTracks->Fill(Jet1->GetNumberOfTracks());
1410  std::vector <Double_t> Jet_Constituents_Data;
1411  AliVParticle *JetConstituent_Data=NULL;
1412  AliVParticle *JetConstituent_Truth=NULL;
1413  if (fMLOn==1){
1414  for (Int_t Tracks_Data=0; Tracks_Data<Jet1->GetNumberOfTracks(); Tracks_Data++){
1415  JetConstituent_Data = static_cast<AliVParticle*>(Jet1->TrackAt(Tracks_Data, JetCont->GetParticleContainer()->GetArray()));
1416  if(!JetConstituent_Data) continue;
1417  Jet_Constituents_Data.push_back(JetConstituent_Data->Px());
1418  Jet_Constituents_Data.push_back(JetConstituent_Data->Py());
1419  Jet_Constituents_Data.push_back(JetConstituent_Data->Pz());
1420  Jet_Constituents_Data.push_back(JetConstituent_Data->E());
1421  }
1422  fShapesVar_Tracks_Rec.push_back(Jet_Constituents_Data);
1423  fShapesVar_Tracks_Truth.push_back(Jet_Constituents_Data);
1424  }
1425  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) fShapesVar[0]= Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
1426  else fShapesVar[0]=Jet1->Pt();
1427  if (fMLOn==0) fShapesVar[2]=FjNSubJettiness(Jet1,0,1,0,1,0);
1428  else fShapesVar[2]=Jet1->Eta();
1429  if (fMLOn==0) fShapesVar[4]=FjNSubJettiness(Jet1,0,2,0,1,0);
1430  else fShapesVar[4]=Jet1->Phi();
1431  if (fMLOn==0) fShapesVar[6]=FjNSubJettiness(Jet1,0,2,0,1,8);
1432  else fShapesVar[6]=Jet1->M();
1433  if (fMLOn==0) fShapesVar[8]=FjNSubJettiness(Jet1,0,2,0,1,1);
1434  else fShapesVar[8]=Jet1->GetNumberOfTracks();
1435  if (fMLOn==0) fShapesVar[10]=FjNSubJettiness(Jet1,0,2,0,1,9);
1436  else fShapesVar[10]=0.0;
1437  if (fMLOn==0) fShapesVar[12]=FjNSubJettiness(Jet1,0,2,0,1,3,fBeta_SD,fZCut);
1438  else fShapesVar[12]=0.0;
1439  if (fMLOn==0) fShapesVar[14]=FjNSubJettiness(Jet1,0,2,0,1,5,fBeta_SD,fZCut);
1440  else fShapesVar[14]=0.0;
1441  if (fMLOn==0) fShapesVar[16]=Jet1->GetLeadingTrack(JetCont->GetParticleContainer()->GetArray())->Pt();
1442  else fShapesVar[16]=0.0;
1443  if (fMLOn==0) fShapesVar[18]=-2; //event plane calculation not needed for data
1444  else fShapesVar[18]=0.0;
1445  if (fMLOn==0) fShapesVar[20]=FjNSubJettiness(Jet1,0,2,0,1,6,fBeta_SD,fZCut);
1446  else fShapesVar[20]=0.0;
1447  AliEmcalJetFinder *Reclusterer1 = Recluster(Jet1, 0, fSubJetRadius, fSubJetMinPt, fSubJetAlgorithm, "SubJetFinder");
1448  if (fFullTree){
1449  if (fMLOn==0) fShapesVar[22]=FjNSubJettiness(Jet1,0,2,0,1,2);
1450  else fShapesVar[22]=0.0;
1451  if (fMLOn==0) fShapesVar[24]=SubJetFraction(Jet1, Reclusterer1, 1, 0, kTRUE, kFALSE);
1452  else fShapesVar[24]=0.0;
1453  if (fMLOn==0) fShapesVar[26]=SubJetFraction(Jet1, Reclusterer1, 2, 0, kTRUE, kFALSE);
1454  else fShapesVar[26]=0.0;
1455  }
1456  fShapesVar[1]=-2;
1457  fShapesVar[3]=-2;
1458  fShapesVar[5]=-2;
1459  fShapesVar[7]=-2;
1460  fShapesVar[9]=-2;
1461  fShapesVar[11]=-2;
1462  fShapesVar[13]=-2;
1463  fShapesVar[15]=-2;
1464  fShapesVar[17]=-2;
1465  fShapesVar[19]=-2;
1466  fShapesVar[21]=-2;
1467  if (fFullTree){
1468  fShapesVar[23]=-2;
1469  fShapesVar[25]=-2;
1470  fShapesVar[27]=-2;
1471  }
1472  fTreeResponseMatrixAxis->Fill();
1473  if (fMLOn==1) fTreeTracks->Fill();
1474  fhSubJetCounter->Fill(Reclusterer1->GetNumberOfJets());
1475  for (Int_t i= 0; i<Reclusterer1->GetNumberOfJets(); i++){
1476  fhEventCounter->Fill(6); //Number of overall subjets in all events
1477  fhSubJetPt->Fill(Reclusterer1->GetJet(i)->Pt());
1478  fhSubJetMass->Fill(Reclusterer1->GetJet(i)->M());
1479  fhNumberOfSubJetTracks->Fill(Reclusterer1->GetJet(i)->GetNumberOfTracks());
1480  fhSubJetRadius->Fill(TMath::Sqrt((Reclusterer1->GetJet(i)->Area()/TMath::Pi()))); //Radius of SubJets per event
1481  }
1482  }
1483  }
1484  fhJetCounter->Fill(JetCounter); //Number of Jets in Each Event
1485  }
1486  //else {fhEventCounter->Fill(2);} //Events with no jet container
1487  }
1488 
1489 
1490  if (fJetShapeType==AliAnalysisTaskSubJetFraction::kSim || fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly){
1491 
1492  AliEmcalJet *Jet1 = NULL; //Original Jet in the event
1493  AliJetContainer *JetCont= GetJetContainer(0); //Jet Container for event
1494  Int_t JetCounter=0; //Counts number of jets in event
1495  Double_t JetPhi=0;
1496  Bool_t EventCounter=kFALSE;
1497  Double_t JetPt_ForThreshold=0;
1498  const AliEmcalPythiaInfo *Parton_Info = 0x0;
1499  Parton_Info=GetPythiaInfo();
1500  fhEventCounter->Fill(1);
1501  if(JetCont) {
1502  fhEventCounter->Fill(2); //Number of events with a jet container
1503  JetCont->ResetCurrentID();
1504  while((Jet1=JetCont->GetNextAcceptJet())) {
1505  // if (((Jet1=JetCont->GetLeadingJet("rho")) && (fPtThreshold<=Jet1->Pt()-(GetRhoVal(0)*Jet1->Area())))){
1506  if(!Jet1) continue;
1507  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) JetPt_ForThreshold = Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
1508  else JetPt_ForThreshold = Jet1->Pt();
1509  if(JetPt_ForThreshold<fPtThreshold) {
1510  //fhEventCounter->Fill(3); //events where the jet had a problem
1511  continue;
1512  }
1513  else {
1514  /* if(fSemigoodCorrect){
1515  Double_t HoleDistance=RelativePhi(Jet1->Phi(),fHolePos);
1516  if(TMath::Abs(HoleDistance)<fHoleWidth) continue;
1517  }*/
1518  Float_t RecoilDeltaPhi = 0.;
1519  if (fJetSelection == kRecoil){
1520  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), Jet1->Phi());
1521  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
1522  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), Jet1->Pt());
1523  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), Jet1->Phi()));
1524  }
1525  if (!EventCounter){
1526  fhEventCounter->Fill(3);
1527  EventCounter=kTRUE;
1528  }
1529  if((Jet1->GetNumberOfTracks())==0){
1530  fhEventCounter->Fill(10); //zero track jets
1531  }
1532  if((Jet1->GetNumberOfTracks())==1){
1533  fhEventCounter->Fill(11); //one track jets
1534  }
1535  fhEventCounter->Fill(4); //Number of Jets found in all events
1536  JetCounter++;
1537  fhJetPt->Fill(Jet1->Pt());
1538  JetPhi=Jet1->Phi();
1539  if(JetPhi < -1*TMath::Pi()) JetPhi += (2*TMath::Pi());
1540  else if (JetPhi > TMath::Pi()) JetPhi -= (2*TMath::Pi());
1541  fhJetPhi->Fill(JetPhi);
1542  fhJetEta->Fill(Jet1->Eta());
1543  fhJetMass->Fill(Jet1->M());
1544  fhJetRadius->Fill(TMath::Sqrt((Jet1->Area()/TMath::Pi()))); //Radius of Jets per event
1545  fhNumberOfJetTracks->Fill(Jet1->GetNumberOfTracks());
1546  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> Randomised_Jet=ModifyJet(Jet1,0,"Randomise");
1547  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_02_15=ModifyJet(Jet1,0,"AddExtraProng_02_15");
1548  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_02_30=ModifyJet(Jet1,0,"AddExtraProng_02_30");
1549  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_02_45=ModifyJet(Jet1,0,"AddExtraProng_02_45");
1550  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_02_60=ModifyJet(Jet1,0,"AddExtraProng_02_60");
1551  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_01_10=ModifyJet(Jet1,0,"AddExtraProng_01_10");
1552  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_01_15=ModifyJet(Jet1,0,"AddExtraProng_01_15");
1553  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_01_20=ModifyJet(Jet1,0,"AddExtraProng_01_20");
1554  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_01_25=ModifyJet(Jet1,0,"AddExtraProng_01_25");
1555  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_01_30=ModifyJet(Jet1,0,"AddExtraProng_01_30");
1556  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_01_45=ModifyJet(Jet1,0,"AddExtraProng_01_45");
1557  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_03_10=ModifyJet(Jet1,0,"AddExtraProng_03_10");
1558  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_03_15=ModifyJet(Jet1,0,"AddExtraProng_03_15");
1559  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_03_20=ModifyJet(Jet1,0,"AddExtraProng_03_20");
1560  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_03_30=ModifyJet(Jet1,0,"AddExtraProng_03_30");
1561  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_03_45=ModifyJet(Jet1,0,"AddExtraProng_03_45");
1562 
1563 
1564  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> kTTrack_1_2_1=ModifyJet(Jet1,0,"AddkTTrack_1_2_1");
1565 
1566  std::vector <Double_t> Jet_Constituents_Data;
1567  AliVParticle *JetConstituent_Data=NULL;
1568  AliVParticle *JetConstituent_Truth=NULL;
1569  if (fMLOn==1){
1570  for (Int_t Tracks_Data=0; Tracks_Data<Jet1->GetNumberOfTracks(); Tracks_Data++){
1571  JetConstituent_Data = static_cast<AliVParticle*>(Jet1->TrackAt(Tracks_Data, JetCont->GetParticleContainer()->GetArray()));
1572  if(!JetConstituent_Data) continue;
1573  Jet_Constituents_Data.push_back(JetConstituent_Data->Px());
1574  Jet_Constituents_Data.push_back(JetConstituent_Data->Py());
1575  Jet_Constituents_Data.push_back(JetConstituent_Data->Pz());
1576  Jet_Constituents_Data.push_back(JetConstituent_Data->E());
1577  }
1578  fShapesVar_Tracks_Rec.push_back(Jet_Constituents_Data);
1579  fShapesVar_Tracks_Truth.push_back(Jet_Constituents_Data);
1580  }
1581 
1582 
1583  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) fShapesVar[0]= Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
1584  else fShapesVar[0]=Jet1->Pt();
1585  if (fMLOn==0) fShapesVar[2]=FjNSubJettiness(Jet1,0,1,0,1,0);
1586  else fShapesVar[2]=Jet1->Eta();
1587  if (fMLOn==0) fShapesVar[4]=FjNSubJettiness(Jet1,0,2,0,1,0);
1588  else fShapesVar[4]=Jet1->Phi();
1589  if (fMLOn==0) fShapesVar[6]=FjNSubJettinessFastJet(ExtraProng_Jet_01_20,0,1,0,1,0);
1590  else fShapesVar[6]=Jet1->M();
1591  if (fMLOn==0) fShapesVar[8]=FjNSubJettiness(Jet1,0,2,0,1,1);
1592  else fShapesVar[8]=Jet1->GetNumberOfTracks();
1593  if (fMLOn==0) fShapesVar[10]=FjNSubJettinessFastJet(ExtraProng_Jet_03_10,0,1,0,1,0);
1594  else fShapesVar[10]=Parton_Info->GetPartonFlag6();
1595  if (fMLOn==0) fShapesVar[12]=FjNSubJettinessFastJet(ExtraProng_Jet_03_20,0,1,0,1,0);
1596  else fShapesVar[12]=Parton_Info->GetPartonEta6();
1597  if (fMLOn==0) fShapesVar[14]=FjNSubJettinessFastJet(ExtraProng_Jet_01_10,0,2,0,1,0);
1598  else fShapesVar[14]=Parton_Info->GetPartonPhi6();
1599  if (fMLOn==0) fShapesVar[16]=FjNSubJettinessFastJet(ExtraProng_Jet_01_20,0,2,0,1,0);
1600  else fShapesVar[16]=Parton_Info->GetPartonFlag7();
1601  // fShapesVar[18]=FjNSubJettinessFastJet(ExtraProng_Jet_01_30,0,2,0,1,0);
1602  if (fMLOn==0) fShapesVar[18]=FjNSubJettinessFastJet(kTTrack_1_2_1,0,2,0,1,0);
1603  else fShapesVar[18]=Parton_Info->GetPartonEta7();
1604  if (fMLOn==0) fShapesVar[20]=FjNSubJettinessFastJet(ExtraProng_Jet_03_15,0,2,0,1,0);
1605  else fShapesVar[20]=Parton_Info->GetPartonPhi7();
1606  AliEmcalJetFinder *Reclusterer1 = Recluster(Jet1, 0, fSubJetRadius, fSubJetMinPt, fSubJetAlgorithm, "SubJetFinder");
1607  if (fFullTree){
1608  if (fMLOn==0) fShapesVar[22]=FjNSubJettiness(Jet1,0,2,0,1,2);
1609  else fShapesVar[22]=0.0;
1610  if (fMLOn==0) fShapesVar[24]=SubJetFraction(Jet1, Reclusterer1, 1, 0, kTRUE, kFALSE);
1611  else fShapesVar[24]=0.0;
1612  if (fMLOn==0) fShapesVar[26]=SubJetFraction(Jet1, Reclusterer1, 2, 0, kTRUE, kFALSE);
1613  else fShapesVar[26]=0.0;
1614  }
1615  if (fMLOn==0) fShapesVar[1]=FjNSubJettinessFastJet(Randomised_Jet,0,1,0,1,0);
1616  else fShapesVar[1]=0.0;
1617  // fShapesVar[3]=FjNSubJettiness(std::unique_ptr<AliEmcalJet>(ModifyJet(Jet1,0,"Randomise")).get(),0,1,0,1,0);
1618  //fShapesVar[5]=FjNSubJettiness(std::unique_ptr<AliEmcalJet>(ModifyJet(Jet1,0,"AddExtraProng_02_30")).get(),0,1,0,1,0);
1619  if (fMLOn==0) fShapesVar[3]=FjNSubJettinessFastJet(ExtraProng_Jet_01_10,0,1,0,1,0);
1620  else fShapesVar[3]=0.0;
1621  if (fMLOn==0) fShapesVar[5]=FjNSubJettinessFastJet(ExtraProng_Jet_01_15,0,1,0,1,0);
1622  else fShapesVar[5]=0.0;
1623  if (fMLOn==0) fShapesVar[7]=FjNSubJettinessFastJet(ExtraProng_Jet_01_25,0,1,0,1,0);
1624  else fShapesVar[7]=0.0;
1625  //fShapesVar[9]=FjNSubJettinessFastJet(ExtraProng_Jet_01_30,0,1,0,1,0);
1626  if (fMLOn==0) fShapesVar[9]=FjNSubJettinessFastJet(kTTrack_1_2_1,0,1,0,1,0);
1627  else fShapesVar[9]=0.0;
1628  if (fMLOn==0) fShapesVar[11]=FjNSubJettinessFastJet(ExtraProng_Jet_03_15,0,1,0,1,0);
1629  else fShapesVar[11]=0.0;
1630  if (fMLOn==0) fShapesVar[13]=FjNSubJettinessFastJet(Randomised_Jet,0,2,0,1,0);
1631  else fShapesVar[13]=0.0;
1632  if (fMLOn==0) fShapesVar[15]=FjNSubJettinessFastJet(ExtraProng_Jet_01_15,0,2,0,1,0);
1633  else fShapesVar[15]=0.0;
1634  if (fMLOn==0) fShapesVar[17]=FjNSubJettinessFastJet(ExtraProng_Jet_01_25,0,2,0,1,0);
1635  else fShapesVar[17]=0.0;
1636  if (fMLOn==0) fShapesVar[19]=FjNSubJettinessFastJet(ExtraProng_Jet_03_10,0,2,0,1,0);
1637  else fShapesVar[19]=0.0;
1638  if (fMLOn==0) fShapesVar[21]=FjNSubJettinessFastJet(ExtraProng_Jet_03_20,0,2,0,1,0);
1639  else fShapesVar[21]=0.0;
1640  if (fFullTree){
1641  if (fMLOn==0) fShapesVar[23]=-2;
1642  else fShapesVar[23]=0.0;
1643  if (fMLOn==0) fShapesVar[25]=-2;
1644  else fShapesVar[25]=0.0;
1645  if (fMLOn==0) fShapesVar[27]=-2;
1646  else fShapesVar[27]=0.0;
1647  }
1648  fTreeResponseMatrixAxis->Fill();
1649  if (fMLOn==0) fTreeTracks->Fill();
1650 
1651  fhSubJetCounter->Fill(Reclusterer1->GetNumberOfJets());
1652  for (Int_t i= 0; i<Reclusterer1->GetNumberOfJets(); i++){
1653  fhEventCounter->Fill(6); //Number of overall subjets in all events
1654  fhSubJetPt->Fill(Reclusterer1->GetJet(i)->Pt());
1655  fhSubJetMass->Fill(Reclusterer1->GetJet(i)->M());
1656  fhNumberOfSubJetTracks->Fill(Reclusterer1->GetJet(i)->GetNumberOfTracks());
1657  fhSubJetRadius->Fill(TMath::Sqrt((Reclusterer1->GetJet(i)->Area()/TMath::Pi()))); //Radius of SubJets per event
1658  }
1659  delete Randomised_Jet.second;
1660  delete ExtraProng_Jet_02_15.second;
1661  delete ExtraProng_Jet_02_30.second;
1662  delete ExtraProng_Jet_02_45.second;
1663  delete ExtraProng_Jet_02_60.second;
1664  delete ExtraProng_Jet_01_30.second;
1665  delete ExtraProng_Jet_01_45.second;
1666  delete ExtraProng_Jet_03_30.second;
1667  delete ExtraProng_Jet_03_45.second;
1668  }
1669  }
1670  fhJetCounter->Fill(JetCounter); //Number of Jets in Each Event
1671  }
1672  //else {fhEventCounter->Fill(2);} //Events with no jet container
1673  }
1674  return kTRUE;
1675 }
1676 //________________________________________________________________________
1677 Double_t AliAnalysisTaskSubJetFraction::RelativePhiEventPlane(Double_t EventPlane, Double_t Phi){
1678 
1679  if(Phi < -1*TMath::Pi()) Phi += (2*TMath::Pi());
1680  else if (Phi > TMath::Pi()) Phi -= (2*TMath::Pi());
1681  Double_t DeltaPhi=Phi-EventPlane;
1682  if(DeltaPhi < -1*TMath::Pi()) DeltaPhi += (2*TMath::Pi());
1683  else if (DeltaPhi > TMath::Pi()) DeltaPhi -= (2*TMath::Pi());
1684  return DeltaPhi;
1685 }
1686 //________________________________________________________________________
1687 Double_t AliAnalysisTaskSubJetFraction::RelativePhi(Double_t Phi1, Double_t Phi2){
1688 
1689  if(Phi1 < -1*TMath::Pi()) Phi1 += (2*TMath::Pi()); // Turns the range of 0to2Pi into -PitoPi ???????????
1690  else if (Phi1 > TMath::Pi()) Phi1 -= (2*TMath::Pi());
1691  if(Phi2 < -1*TMath::Pi()) Phi2 += (2*TMath::Pi());
1692  else if (Phi2 > TMath::Pi()) Phi2 -= (2*TMath::Pi());
1693  Double_t DeltaPhi=Phi2-Phi1;
1694  if(DeltaPhi < -1*TMath::Pi()) DeltaPhi += (2*TMath::Pi());
1695  else if (DeltaPhi > TMath::Pi()) DeltaPhi -= (2*TMath::Pi());
1696  return DeltaPhi;
1697 }
1698 
1699 
1700 //--------------------------------------------------------------------------
1701 Int_t AliAnalysisTaskSubJetFraction::SelectTriggerHadron(Float_t PtMin, Float_t PtMax){
1702 
1703  AliTrackContainer *PartCont = NULL;
1704  AliParticleContainer *PartContMC = NULL;
1705 
1706 
1707  if (fJetShapeSub==kConstSub){
1708  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) PartContMC = GetParticleContainer(1);
1709  else PartCont = GetTrackContainer(1);
1710  }
1711  else{
1712  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) PartContMC = GetParticleContainer(0);
1713  else PartCont = GetTrackContainer(0);
1714  }
1715 
1716  TClonesArray *TracksArray = NULL;
1717  TClonesArray *TracksArrayMC = NULL;
1718  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) TracksArrayMC = PartContMC->GetArray();
1719  else TracksArray = PartCont->GetArray();
1720 
1721  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly){
1722  if(!PartContMC || !TracksArrayMC) return -99999;
1723  }
1724  else {
1725  if(!PartCont || !TracksArray) return -99999;
1726  }
1727 
1728  AliAODTrack *Track = 0x0;
1729  Int_t Trigger_Index[100];
1730  for (Int_t i=0; i<100; i++) Trigger_Index[i] = 0;
1731  Int_t Trigger_Counter = 0;
1732  Int_t NTracks=0;
1733  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) NTracks = TracksArrayMC->GetEntriesFast();
1734  else NTracks = TracksArray->GetEntriesFast();
1735  for(Int_t i=0; i < NTracks; i++){
1736  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly){
1737  if((Track = static_cast<AliAODTrack*>(PartContMC->GetAcceptParticle(i)))){
1738  if (!Track) continue;
1739  if(TMath::Abs(Track->Eta())>0.9) continue;
1740  if (Track->Pt()<0.15) continue;
1741  if ((Track->Pt() >= PtMin) && (Track->Pt()< PtMax)) {
1742  Trigger_Index[Trigger_Counter] = i;
1743  Trigger_Counter++;
1744  }
1745  }
1746  }
1747  else{
1748  if((Track = static_cast<AliAODTrack*>(PartCont->GetAcceptTrack(i)))){
1749  if (!Track) continue;
1750  if(TMath::Abs(Track->Eta())>0.9) continue;
1751  if (Track->Pt()<0.15) continue;
1752  if ((Track->Pt() >= PtMin) && (Track->Pt()< PtMax)) {
1753  Trigger_Index[Trigger_Counter] = i;
1754  Trigger_Counter++;
1755  }
1756  }
1757  }
1758  }
1759  if (Trigger_Counter == 0) return -99999;
1760  Int_t RandomNumber = 0, Index = 0 ;
1761  TRandom3* Random = new TRandom3(0);
1762  RandomNumber = Random->Integer(Trigger_Counter);
1763  Index = Trigger_Index[RandomNumber];
1764  return Index;
1765 }
1766 
1767 
1768 //--------------------------------------------------------------------------
1769 Double_t AliAnalysisTaskSubJetFraction::Angularity(AliEmcalJet *Jet, Int_t JetContNb){
1770 
1771  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1772  Double_t Angularity_Numerator=0; //Reset these values
1773  Double_t Angularity_Denominator=0;
1774  AliVParticle *Particle=0x0;
1775  Double_t DeltaPhi=0;
1776 
1777 
1778  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks (particles in the jet
1779  Particle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
1780  if(!Particle) continue;
1781  DeltaPhi=RelativePhi(Jet->Phi(),Particle->Phi());
1782  Angularity_Numerator=Angularity_Numerator+(Particle->Pt()*TMath::Sqrt(((Particle->Eta()-Jet->Eta())*(Particle->Eta()-Jet->Eta()))+(DeltaPhi*DeltaPhi)));
1783  Angularity_Denominator= Angularity_Denominator+Particle->Pt();
1784  }
1785  if(Angularity_Denominator!=0) return Angularity_Numerator/Angularity_Denominator;
1786  else return -1;
1787 
1788 }
1789 
1790 
1791 
1792 //--------------------------------------------------------------------------
1793 Double_t AliAnalysisTaskSubJetFraction::PTD(AliEmcalJet *Jet, Int_t JetContNb){
1794 
1795  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1796  Double_t PTD_Numerator=0; //Reset these values
1797  Double_t PTD_Denominator=0;
1798  AliVParticle *Particle=0x0;
1799  Double_t DeltaPhi=0;
1800  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks (particles in the jet
1801  Particle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
1802  if(!Particle) continue;
1803  PTD_Numerator=PTD_Numerator+(Particle->Pt()*Particle->Pt());
1804  PTD_Denominator=PTD_Denominator+Particle->Pt();
1805  }
1806  if(PTD_Denominator!=0) return TMath::Sqrt(PTD_Numerator)/PTD_Denominator;
1807  else return -1;
1808 
1809 }
1810 
1811 
1812 //----------------------------------------------------------------------
1814 AliEmcalJetFinder *AliAnalysisTaskSubJetFraction::Recluster(AliEmcalJet *Jet, Int_t JetContNb, Double_t SubJetRadius, Double_t SubJetMinPt, Int_t Algorithm, const char* Name){
1815 
1816  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1817  AliEmcalJetFinder *Reclusterer = new AliEmcalJetFinder(Name); //JetFinder Object for reclustered jets
1818  Reclusterer->SetRadius(SubJetRadius);
1819  Reclusterer->SetJetMinPt(SubJetMinPt);
1820  Reclusterer->SetJetAlgorithm(Algorithm); //0 for anti-kt 1 for kt
1821  Reclusterer->SetJetMaxEta(0.9);
1822  Reclusterer->SetRecombSheme(0);
1823  if(fJetShapeType != AliAnalysisTaskSubJetFraction::kGenOnTheFly){
1824  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
1825  Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
1826  if(Reclusterer->AliEmcalJetFinder::Filter(Jet, JetCont, dVtx)){;} //reclustering jet1 using the jetfinderobject Reclusterer
1827  }
1828  else{
1829  Double_t dVtx[3]={1,1,1};
1830  if(Reclusterer->AliEmcalJetFinder::Filter(Jet, JetCont, dVtx)){;} //reclustering jet1 using the jetfinderobject Reclusterer
1831  }
1832  return Reclusterer;
1833 
1834 }
1835 
1836 
1837 
1838 //----------------------------------------------------------------------
1839 Double_t AliAnalysisTaskSubJetFraction::SubJetOrdering(AliEmcalJet *Jet, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t Type, Bool_t Index){
1840  AliEmcalJet *SubJet=NULL;
1841  Double_t SortingVariable;
1842  Int_t ArraySize =N+1;
1843  TArrayD JetSorter(ArraySize);
1844  TArrayD JetIndexSorter(ArraySize);
1845  for (Int_t i=0; i<ArraySize; i++){
1846  JetSorter[i]=0;
1847  }
1848  for (Int_t i=0; i<ArraySize; i++){
1849  JetIndexSorter[i]=0;
1850  }
1851  if(Reclusterer->GetNumberOfJets()<N) return -999;
1852  for (Int_t i=0; i<Reclusterer->GetNumberOfJets(); i++){
1853  SubJet=Reclusterer->GetJet(i);
1854  if (Type==0) SortingVariable=SubJet->Pt();
1855  else if (Type==1) SortingVariable=SubJet->E();
1856  else if (Type==2) SortingVariable=SubJet->M();
1857  for (Int_t j=0; j<N; j++){
1858  if (SortingVariable>JetSorter[j]){
1859  for (Int_t k=N-1; k>=j; k--){
1860  JetSorter[k+1]=JetSorter[k];
1861  JetIndexSorter[k+1]=JetIndexSorter[k];
1862  }
1863  JetSorter[j]=SortingVariable;
1864  JetIndexSorter[j]=i;
1865  break;
1866  }
1867  }
1868  }
1869  if (!Index) return JetSorter[N-1];
1870  else return JetIndexSorter[N-1];
1871 }
1872 
1873 
1874 
1875 //returns -1 if the Nth hardest jet is requested where N>number of available jets
1876 //type: 0=Pt 1=E 2=M
1877 //Index TRUE=returns index FALSE=returns value of quantatiy in question
1878 
1879 
1880 
1881 
1882 
1883 
1884 
1885 //----------------------------------------------------------------------
1886 Double_t AliAnalysisTaskSubJetFraction::SubJetFraction(AliEmcalJet *Jet, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t Type, Bool_t Add, Bool_t Loss){
1887  AliEmcalJet *SubJet=NULL;
1888  Double_t Observable=0;
1889  Double_t Fraction_Numerator=0;
1890  Bool_t Error=kFALSE;
1891  if (!Jet->GetNumberOfTracks()) return -2;
1892  if (Add){
1893  for (Int_t i=1; i<=N; i++){
1894  Observable=SubJetOrdering(Jet,Reclusterer,i,Type,kFALSE);
1895  if(Observable==-999){
1896  Error = kTRUE;
1897  return -2;
1898  }
1899  Fraction_Numerator=Fraction_Numerator+Observable;
1900  }
1901  }
1902  else {
1903  Fraction_Numerator=SubJetOrdering(Jet,Reclusterer,N,Type,kFALSE);
1904  if (Fraction_Numerator==-999) return -2;
1905  }
1906  if (Type==0){
1907  if(Loss) return (Jet->Pt()-Fraction_Numerator)/Jet->Pt();
1908  else return Fraction_Numerator/Jet->Pt();
1909  }
1910  else if (Type==1){
1911  if(Loss) return (Jet->E()-Fraction_Numerator)/Jet->E();
1912  else return Fraction_Numerator/Jet->E();
1913  }
1914  else { //change to else if if you want to add more later
1915  if(Loss) return (Jet->M()-Fraction_Numerator)/Jet->M();
1916  else return Fraction_Numerator/Jet->M();
1917  }
1918 }
1919 //N number of hardest subjets involved
1920 //Type 0=Pt 1=E 2=M
1921 // Add TRUE: Add 1 to Nth hardest subjets together/total jet False: Nth hardest Jet/toal jet
1922 //Loss TRUE: Jet-Subjet(s)/Jet FALSE: Subjet(s)/Jet
1923 
1924 
1925 //----------------------------------------------------------------------
1926 Double_t AliAnalysisTaskSubJetFraction::NSubJettiness(AliEmcalJet *Jet, Int_t JetContNb, Double_t JetRadius, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t A, Int_t B){
1927  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1928  AliEmcalJet *SubJet=NULL;
1929  Double_t DeltaR1=0;
1930  Double_t DeltaR2=0;
1931  AliVParticle *JetParticle=0x0;
1932  Double_t SubJetiness_Numerator = 0;
1933  Double_t SubJetiness_Denominator = 0;
1934  Double_t Index=-2;
1935  Bool_t Error=kFALSE;
1936  // JetRadius=TMath::Sqrt((Jet->Area()/TMath::Pi())); //comment out later
1937  if (!Jet->GetNumberOfTracks()) return -2;
1938  if (Reclusterer->GetNumberOfJets() < N) return -2;
1939  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks (particles in the jet
1940  JetParticle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
1941  for (Int_t j=1; j<=N; j++){
1942  Index=SubJetOrdering(Jet,Reclusterer,j,0,kTRUE);
1943  if(Index==-999){
1944  Error = kTRUE;
1945  i=Jet->GetNumberOfTracks();
1946  break;
1947  }
1948  if(j==1){
1949  DeltaR1=TMath::Power((Reclusterer->GetJet(Index)->Pt()),A)*TMath::Power((TMath::Sqrt((((JetParticle->Eta())-(Reclusterer->GetJet(Index)->Eta()))*((JetParticle->Eta())- (Reclusterer->GetJet(Index)->Eta())))+((RelativePhi((Reclusterer->GetJet(Index)->Phi()),JetParticle->Phi()))*(RelativePhi((Reclusterer->GetJet(Index)->Phi()),JetParticle->Phi()))))),B);
1950  }
1951  else{
1952  DeltaR2=TMath::Power((Reclusterer->GetJet(Index)->Pt()),A)*TMath::Power((TMath::Sqrt((((JetParticle->Eta())-(Reclusterer->GetJet(Index)->Eta()))*((JetParticle->Eta())- (Reclusterer->GetJet(Index)->Eta())))+((RelativePhi((Reclusterer->GetJet(Index)->Phi()),JetParticle->Phi()))*(RelativePhi((Reclusterer->GetJet(Index)->Phi()),JetParticle->Phi()))))),B);
1953  if (DeltaR2<DeltaR1) DeltaR1=DeltaR2;
1954  }
1955  }
1956  SubJetiness_Numerator=SubJetiness_Numerator+(JetParticle->Pt()*DeltaR1);
1957  if (A>=0) SubJetiness_Denominator=SubJetiness_Denominator+(TMath::Power((Reclusterer->GetJet(SubJetOrdering(Jet,Reclusterer,1,0,kTRUE))->Pt()),A)*JetParticle->Pt()*TMath::Power(JetRadius,B));
1958  else SubJetiness_Denominator=SubJetiness_Denominator+(TMath::Power((Reclusterer->GetJet(SubJetOrdering(Jet,Reclusterer,N,0,kTRUE))->Pt()),A)*JetParticle->Pt()*TMath::Power(JetRadius,B));
1959  }
1960  if (SubJetiness_Denominator!=0 && !Error){
1961  // if (SubJetiness_Numerator/SubJetiness_Denominator>1.0) cout << JetRadius<<endl;
1962  return SubJetiness_Numerator/SubJetiness_Denominator; }
1963  else return -2;
1964 }
1965 
1966 
1967 //______________________________________________________________________________________
1968 Double_t AliAnalysisTaskSubJetFraction::FjNSubJettiness(AliEmcalJet *Jet, Int_t JetContNb,Int_t N, Int_t Algorithm, Double_t Beta, Int_t Option, Double_t Beta_SD, Double_t ZCut){
1969 
1970  //WARNING!!! Only works for parent jets that are clustered with Anti-Kt! To change go to AliEmcalJetFinder.cxx and look at the Nsubjettiness() function
1971 
1972  //Algorithm==0 -> kt_axes;
1973  // Algorithm==1 -> ca_axes;
1974  //Algorithm==2 -> antikt_0p2_axes;
1975  //Algorithm==3 -> wta_kt_axes;
1976  //Algorithm==4 -> wta_ca_axes;
1977  //Algorithm==5 -> onepass_kt_axes;
1978  //Algorithm==6 -> onepass_ca_axes;
1979  //Algorithm==7 -> onepass_antikt_0p2_axes;
1980  //Algorithm==8 -> onepass_wta_kt_axes;
1981  //Algorithm==9 -> onepass_wta_ca_axes;
1982  //Algorithm==10 -> min_axes;
1983 
1984 
1985  //fSoftDropOn==1 All the below are done on the jet after it has undergone softdrop with CA. Option 3 and 4 are done whether this is turned on or not.
1986 
1987 
1988  //Option==0 returns Nsubjettiness Value
1989  //Option==1 && N==2 returns opening angle between two subjet axes(Delta R?)
1990  //Option==2 && N==2 returns Delta R
1991  //Option==3 returns first splitting distance for soft dropped jet (CA)
1992  //Option==4 returns Symmetry measure (Zg) for soft dropped jet (CA)
1993  //Option==5 returns Pt of Subjet1
1994  //Option==6 returns Pt of Subjet2
1995  //Options==7 trutns deltaR of subjets...Is this different to before??
1996  //Option==8 Subjet1 Leading track Pt
1997  //Option==9 Subjet1 Leading track Pt
1998 
1999 
2000  Algorithm=fReclusteringAlgorithm;
2001  if (Jet->GetNumberOfTracks()>=N){
2002  if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==1) && (Algorithm==0) && (Beta==1.0) && (Option==0)){
2003  if (fDerivSubtrOrder == kFirstOrder) return Jet->GetShapeProperties()->GetFirstOrderSubtracted1subjettiness_kt();
2004  else return Jet->GetShapeProperties()->GetSecondOrderSubtracted1subjettiness_kt();
2005  }
2006  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==0) && (Beta==1.0) && (Option==0)){
2007  if (fDerivSubtrOrder == kFirstOrder) return Jet->GetShapeProperties()->GetFirstOrderSubtracted2subjettiness_kt();
2008  else return Jet->GetShapeProperties()->GetSecondOrderSubtracted2subjettiness_kt();
2009  }
2010  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==3) && (Algorithm==0) && (Beta==1.0) && (Option==0)){
2011  if (fDerivSubtrOrder == kFirstOrder) return Jet->GetShapeProperties()->GetFirstOrderSubtracted3subjettiness_kt();
2012  else return Jet->GetShapeProperties()->GetSecondOrderSubtracted3subjettiness_kt();
2013  }
2014  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==0) && (Beta==1.0) && (Option==1)){
2015  if (fDerivSubtrOrder == kFirstOrder) return Jet->GetShapeProperties()->GetFirstOrderSubtractedOpeningAngle_kt();
2016  else return Jet->GetShapeProperties()->GetSecondOrderSubtractedOpeningAngle_kt();
2017  }
2018  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==1) && (Algorithm==1) && (Beta==1.0) && (Option==0)){
2019  if (fDerivSubtrOrder == kFirstOrder) return Jet->GetShapeProperties()->GetFirstOrderSubtracted1subjettiness_ca();
2020  else return Jet->GetShapeProperties()->GetSecondOrderSubtracted1subjettiness_ca();
2021  }
2022  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==1) && (Beta==1.0) && (Option==0)){
2023  if (fDerivSubtrOrder == kFirstOrder) return Jet->GetShapeProperties()->GetFirstOrderSubtracted2subjettiness_ca();
2024  else return Jet->GetShapeProperties()->GetSecondOrderSubtracted2subjettiness_ca();
2025  }
2026  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==1) && (Beta==1.0) && (Option==1)){
2027  if (fDerivSubtrOrder == kFirstOrder) return Jet->GetShapeProperties()->GetFirstOrderSubtractedOpeningAngle_ca();
2028  else return Jet->GetShapeProperties()->GetSecondOrderSubtractedOpeningAngle_ca();
2029  }
2030  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==1) && (Algorithm==2) && (Beta==1.0) && (Option==0)){
2031  if (fDerivSubtrOrder == kFirstOrder) return Jet->GetShapeProperties()->GetFirstOrderSubtracted1subjettiness_akt02();
2032  else return Jet->GetShapeProperties()->GetSecondOrderSubtracted1subjettiness_akt02();
2033  }
2034  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==2) && (Beta==1.0) && (Option==0)){
2035  if (fDerivSubtrOrder == kFirstOrder) return Jet->GetShapeProperties()->GetFirstOrderSubtracted2subjettiness_akt02();
2036  else return Jet->GetShapeProperties()->GetSecondOrderSubtracted2subjettiness_akt02();
2037  }
2038  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==2) && (Beta==1.0) && (Option==1)){
2039  if (fDerivSubtrOrder == kFirstOrder) return Jet->GetShapeProperties()->GetFirstOrderSubtractedOpeningAngle_akt02();
2040  else return Jet->GetShapeProperties()->GetSecondOrderSubtractedOpeningAngle_akt02();
2041  }
2042  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==1) && (Algorithm==6) && (Beta==1.0) && (Option==0)){
2043  if (fDerivSubtrOrder == kFirstOrder) return Jet->GetShapeProperties()->GetFirstOrderSubtracted1subjettiness_onepassca();
2044  else return Jet->GetShapeProperties()->GetSecondOrderSubtracted1subjettiness_onepassca();
2045  }
2046  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==6) && (Beta==1.0) && (Option==0)){
2047  if (fDerivSubtrOrder == kFirstOrder) return Jet->GetShapeProperties()->GetFirstOrderSubtracted2subjettiness_onepassca();
2048  else return Jet->GetShapeProperties()->GetSecondOrderSubtracted2subjettiness_onepassca();
2049  }
2050  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==6) && (Beta==1.0) && (Option==1)){
2051  if (fDerivSubtrOrder == kFirstOrder) return Jet->GetShapeProperties()->GetFirstOrderSubtractedOpeningAngle_onepassca();
2052  else return Jet->GetShapeProperties()->GetSecondOrderSubtractedOpeningAngle_onepassca();
2053  }
2054  else{
2055  //Algorithm=fReclusteringAlgorithm;
2056  AliJetContainer *JetCont = GetJetContainer(JetContNb);
2057  AliEmcalJetFinder JetFinder("Nsubjettiness");
2058  JetFinder.SetJetMaxEta(0.9-fJetRadius);
2059  JetFinder.SetRadius(fJetRadius);
2060  JetFinder.SetJetAlgorithm(0); //0 for anti-kt 1 for kt //this is for the JET!!!!!!!!!! Not the SubJets
2061  JetFinder.SetRecombSheme(0);
2062  JetFinder.SetJetMinPt(Jet->Pt());
2063  if(fJetShapeType != AliAnalysisTaskSubJetFraction::kGenOnTheFly){
2064  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
2065  Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
2066  return JetFinder.Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,fSubJetRadius,Beta,Option,0,Beta_SD,ZCut,fSoftDropOn);
2067  }
2068  else{
2069  Double_t dVtx[3]={1,1,1};
2070  if (!fNsubMeasure){
2071  return JetFinder.Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,fSubJetRadius,Beta,Option,0,Beta_SD,ZCut,fSoftDropOn);
2072  }
2073  else{
2074  if (Option==3) return JetFinder.Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,fSubJetRadius,Beta,Option,0,Beta_SD,ZCut,fSoftDropOn);
2075  else return JetFinder.Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,fSubJetRadius,Beta,Option,1);
2076  }
2077  }
2078  }
2079  }
2080  else return -2;
2081 }
2082 
2083 
2084 //________________________________________________________________________
2085 Bool_t AliAnalysisTaskSubJetFraction::RetrieveEventObjects() {
2086  //
2087  // retrieve event objects
2088  //
2089  if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
2090  return kFALSE;
2091 
2092  return kTRUE;
2093 }
2094 
2095 
2096 //_______________________________________________________________________
2097 void AliAnalysisTaskSubJetFraction::Terminate(Option_t *)
2098 {
2099  // Called once at the end of the analysis.
2100  if (fJetShapeType==AliAnalysisTaskSubJetFraction::kData){
2101 
2102  /*
2103  for (int i=1; i<=(XBinsJetPtSize)-1; i++){ //Rescales the fhJetPt Histograms according to the with of the variable bins and number of events
2104  fhJetPt->SetBinContent(i,((1.0/(XBinsJetPt[i]-XBinsJetPt[i-1]))*((fhJetPt->GetBinContent(i))*(1.0/(fhEventCounter->GetBinContent(1))))));
2105  fhSubJetPt->SetBinContent(i,((1.0/(XBinsJetPt[i]-XBinsJetPt[i-1]))*((fhSubJetPt->GetBinContent(i))*(1.0/(fhEventCounter->GetBinContent(5))))));
2106 
2107  //fhJetPt->SetBinContent(i,((1.0/(fhPt->GetBinWidth(i)))*((fhJetPt->GetBinContent(i))*(fhEventCounter->GetBinContent(1)))));
2108  // fhJetPt->SetBinContent(i,((1.0/((fhPt->GetLowEdge(i+1))-(fhJetPt->GetLowEdge(i))))*((fhPt->GetBinContent(i))*(fhEventCounter->GetBinContent(1)))));
2109  }
2110  for (int i=1; i<=(XBinsJetMassSize)-1; i++){ //Rescales the fhPt Histograms according to the with of the variable bins and number of events
2111  fhJetMass->SetBinContent(i,((1.0/(XBinsJetMass[i]-XBinsJetMass[i-1]))*((fhJetMass->GetBinContent(i))*(1.0/(fhEventCounter->GetBinContent(1))))));
2112  }
2113 
2114  fhJetPhi->Scale(Phi_Bins/((Phi_Up-Phi_Low)*((fhEventCounter->GetBinContent(1)))));
2115  fhJetEta->Scale(Eta_Bins/((Eta_Up-Eta_Low)*((fhEventCounter->GetBinContent(1)))));
2116  fhJetRadius->Scale(100/(fhEventCounter->GetBinContent(4))); //should this and JetAngularity be divided by Bin 1 or 4????
2117  fhNumberOfJetTracks->Scale(1.0/(fhEventCounter->GetBinContent(4)));
2118  fhJetCounter->Scale(1.0/(fhEventCounter->GetBinContent(1))); //is the first bin the correct one to look at?
2119 
2120  fhSubJetCounter->Scale(1.0/(fhEventCounter->GetBinContent(5)));
2121  */
2122  /*
2123  fhJetPt->Scale(1.0/(fhEventCounter->GetBinContent(1)));
2124  fhSubJetPt->Scale(1.0/(fhEventCounter->GetBinContent(5)));
2125  fhJetMass->Scale(1.0/(fhEventCounter->GetBinContent(1)));
2126  fhJetPhi->Scale(1.0/(fhEventCounter->GetBinContent(1)));
2127  fhJetEta->Scale(1.0/(fhEventCounter->GetBinContent(1)));
2128  fhJetRadius->Scale(1.0/(fhEventCounter->GetBinContent(4))); //should this and JetAngularity be divided by Bin 1 or 4????
2129  fhNumberOfJetTracks->Scale(1.0/(fhEventCounter->GetBinContent(4)));
2130  fhJetCounter->Scale(1.0/(fhEventCounter->GetBinContent(1))); //is the first bin the correct one to look at?
2131  fhSubJetCounter->Scale(1.0/(fhEventCounter->GetBinContent(5)));
2132  */
2133  }
2134 
2135 }
2136 
2138 std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> AliAnalysisTaskSubJetFraction::ModifyJet(AliEmcalJet* Jet, Int_t JetContNb, TString Modification){
2139  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> Jet_ClusterSequence;
2140  std::vector<fastjet::PseudoJet> fInputVectors;
2141  AliJetContainer *fJetCont = GetJetContainer(JetContNb);
2142  AliParticleContainer *fTrackCont = fJetCont->GetParticleContainer();
2143  std::vector<fastjet::PseudoJet> Modified_Jet;
2144  fastjet::ClusterSequence *fClustSeqSA;
2145  Int_t NJetTracksTest = Jet->GetNumberOfTracks();
2146  if (fTrackCont) for (Int_t i=0; i<Jet->GetNumberOfTracks(); i++) {
2147  AliVParticle *fTrk = Jet->TrackAt(i, fTrackCont->GetArray());
2148  if (!fTrk) continue;
2149  fastjet::PseudoJet PseudoTracks(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
2150  PseudoTracks.set_user_index(Jet->TrackAt(i)+100); //100 is very arbitary....why is it here anyway?
2151  fInputVectors.push_back(PseudoTracks);
2152  }
2153  if (Modification=="Randomise") fInputVectors=RandomiseTracks(Jet,fInputVectors);
2154  if (Modification=="AddExtraProng_02_15") fInputVectors=AddExtraProng(fInputVectors,0.2,0.15);
2155  if (Modification=="AddExtraProng_02_30") fInputVectors=AddExtraProng(fInputVectors,0.2,0.30);
2156  if (Modification=="AddExtraProng_02_45") fInputVectors=AddExtraProng(fInputVectors,0.2,0.45);
2157  if (Modification=="AddExtraProng_02_60") fInputVectors=AddExtraProng(fInputVectors,0.2,0.60);
2158  if (Modification=="AddExtraProng_01_10") fInputVectors=AddExtraProng(fInputVectors,0.1,0.10);
2159  if (Modification=="AddExtraProng_01_15") fInputVectors=AddExtraProng(fInputVectors,0.1,0.15);
2160  if (Modification=="AddExtraProng_01_20") fInputVectors=AddExtraProng(fInputVectors,0.1,0.20);
2161  if (Modification=="AddExtraProng_01_25") fInputVectors=AddExtraProng(fInputVectors,0.1,0.25);
2162  if (Modification=="AddExtraProng_01_30") fInputVectors=AddExtraProng(fInputVectors,0.1,0.30);
2163  if (Modification=="AddExtraProng_01_45") fInputVectors=AddExtraProng(fInputVectors,0.1,0.45);
2164  if (Modification=="AddExtraProng_03_10") fInputVectors=AddExtraProng(fInputVectors,0.3,0.10);
2165  if (Modification=="AddExtraProng_03_15") fInputVectors=AddExtraProng(fInputVectors,0.3,0.15);
2166  if (Modification=="AddExtraProng_03_20") fInputVectors=AddExtraProng(fInputVectors,0.3,0.20);
2167  if (Modification=="AddExtraProng_03_30") fInputVectors=AddExtraProng(fInputVectors,0.3,0.30);
2168  if (Modification=="AddExtraProng_03_45") fInputVectors=AddExtraProng(fInputVectors,0.3,0.45);
2169  if (Modification=="AddkTTrack_1_2_1") fInputVectors=AddkTTracks(Jet,fInputVectors,1.0,2.0,1);
2170  fJetRadius=fJetRadius*2.0;
2171  try {
2172  fastjet::JetDefinition fJetDef(fastjet::antikt_algorithm, fJetRadius, static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
2173  fClustSeqSA=new fastjet::ClusterSequence(fInputVectors, fJetDef);
2174  Modified_Jet=fClustSeqSA->inclusive_jets(0);
2175  Jet_ClusterSequence.second=fClustSeqSA;
2176  } catch (fastjet::Error) {
2177  AliError(" [w] FJ Exception caught.");
2178  //return -1;
2179  }
2180  fJetRadius=fJetRadius/2.0;
2181  Jet_ClusterSequence.first=Modified_Jet[0];
2182  return Jet_ClusterSequence;
2183 }
2184 
2185 
2186 std::vector<fastjet::PseudoJet> AliAnalysisTaskSubJetFraction::RandomiseTracks(AliEmcalJet *Jet,std::vector<fastjet::PseudoJet> fInputVectors){
2187  TRandom3 fRandom1;
2188  fRandom1.SetSeed(0);
2189  Double_t Random_Phi=0.0;
2190  // Double_t Random_Phi_Change=0.0;
2191  Double_t Random_Eta=0.0;
2192  // Double_t Random_Eta_Change=0.0;
2193  Double_t Track_Pt=0.0;
2194  std::vector<fastjet::PseudoJet> Random_Track_Vector;
2195  Random_Track_Vector.clear();
2196  for (Int_t i=0; i< fInputVectors.size(); i++){
2197  do{
2198  Random_Phi=fRandom1.Uniform(Jet->Phi()-fJetRadius,Jet->Phi()+fJetRadius);
2199  Random_Eta=fRandom1.Uniform(Jet->Eta()-fJetRadius,Jet->Eta()+fJetRadius);
2200  }while(TMath::Sqrt(TMath::Power(RelativePhi(Jet->Phi(),Random_Phi),2)+TMath::Power(Jet->Eta()-Random_Eta,2))>fJetRadius);
2201  /* Random_Phi_Change=fRandom1.Uniform(-1*fJetRadius,fJetRadius);
2202  Random_Phi=Jet->Phi()+Random_Phi_Change;
2203  Random_Eta_Change=fRandom1.Uniform(-1*(TMath::Sqrt((fJetRadius*fJetRadius)-(Random_Phi_Change*Random_Phi_Change))),TMath::Sqrt((fJetRadius*fJetRadius)-(Random_Phi_Change*Random_Phi_Change)));
2204  Random_Eta=Jet->Eta()+Random_Eta_Change;
2205  */
2206  if (fRandomisationEqualPt) Track_Pt=Jet->Pt()/fInputVectors.size();
2207  else Track_Pt=fInputVectors[i].perp();
2208  fastjet::PseudoJet Random_Track(Track_Pt*TMath::Cos(Random_Phi),Track_Pt*TMath::Sin(Random_Phi),Track_Pt*TMath::SinH(Random_Eta),TMath::Sqrt((Track_Pt*Track_Pt)+(Track_Pt*TMath::SinH(Random_Eta)*Track_Pt*TMath::SinH(Random_Eta))));
2209  Random_Track.set_user_index(i);
2210  Random_Track_Vector.push_back(Random_Track);
2211  }
2212  fInputVectors.clear();
2213  return Random_Track_Vector;
2214 }
2215 
2216 std::vector<fastjet::PseudoJet> AliAnalysisTaskSubJetFraction::AddExtraProng(std::vector<fastjet::PseudoJet> fInputVectors, Double_t Distance, Double_t PtFrac){
2217  TRandom3 fRandom1;
2218  fRandom1.SetSeed(0);
2219  std::vector<fastjet::PseudoJet> Jet;
2220  Jet.clear();
2221  try {
2222  fastjet::JetDefinition fJetDef(fastjet::antikt_algorithm, fJetRadius*2, static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
2223  fastjet::ClusterSequence fClustSeqSA_Prong(fInputVectors, fJetDef);
2224  Jet= fClustSeqSA_Prong.inclusive_jets(0);
2225  } catch (fastjet::Error) {
2226  AliError(" [w] FJ Exception caught.");
2227  }
2228  Double_t Extra_Track_Phi_Change=fRandom1.Uniform(-1*Distance,Distance);
2229  Double_t Extra_Track_Phi=Jet[0].phi()+Extra_Track_Phi_Change;
2230  Double_t Extra_Track_Eta_Change=fRandom1.Uniform(-1*(TMath::Sqrt((Distance*Distance)-(Extra_Track_Phi_Change*Extra_Track_Phi_Change))),TMath::Sqrt((Distance*Distance)-(Extra_Track_Phi_Change*Extra_Track_Phi_Change)));
2231  Double_t Extra_Track_Eta=Jet[0].pseudorapidity()+Extra_Track_Eta_Change;
2232  Double_t Extra_Track_Pt=Jet[0].perp()*PtFrac;
2233  fastjet::PseudoJet ExtraTrack(Extra_Track_Pt*TMath::Cos(Extra_Track_Phi),Extra_Track_Pt*TMath::Sin(Extra_Track_Phi),Extra_Track_Pt*TMath::SinH(Extra_Track_Eta),TMath::Sqrt((Extra_Track_Pt*Extra_Track_Pt)+(Extra_Track_Pt*TMath::SinH(Extra_Track_Eta)*Extra_Track_Pt*TMath::SinH(Extra_Track_Eta))));
2234  ExtraTrack.set_user_index(150);
2235  fInputVectors.push_back(ExtraTrack);
2236  return fInputVectors;
2237 }
2238 
2239 std::vector<fastjet::PseudoJet> AliAnalysisTaskSubJetFraction::AddkTTracks(AliEmcalJet *Jet,std::vector<fastjet::PseudoJet> fInputVectors, Double_t QHat,Double_t XLength, Int_t NAdditionalTracks)
2240 {
2241 //here add N tracks with random phi and eta and theta according to bdmps distrib.
2242  fastjet::PseudoJet Lab_Jet;
2243  Lab_Jet.reset(Jet->Px(),Jet->Py(),Jet->Pz(),Jet->E());
2244  // Double_t Omega_C=0.5*QHat*XLength*XLength/0.2;
2245  // Double_t Theta_C=TMath::Sqrt(12*0.2/(QHat*TMath::Power(XLength,3)));
2246  Double_t xQs=TMath::Sqrt(QHat*XLength); //QHat=1,XLength=2
2247  //cout<<"medium parameters "<<Omega_C<<" "<<Theta_C<<" "<<xQs<<endl;
2248 
2249  for(Int_t i=0;i<NAdditionalTracks;i++){
2250 
2251  Double_t kT_Scale,Limit_Min,Limit_Max;
2252 
2253  //generation of kT according to 1/kT^4, with minimum QS^2=2 GeV and maximum ~sqrt(ptjet*T)
2254  TF1 *kT_Func= new TF1("kT_Func","1/(x*x*x*x)",xQs*xQs,10000);
2255  kT_Scale=kT_Func->GetRandom();
2256  //generation within the jet cone
2257 
2258  //generation of w according to 1/w, with minimum wc
2259  //omega needs to be larger than kT so to have well defined angles
2260  Limit_Min=kT_Scale;
2261  Limit_Max=kT_Scale/TMath::Sin(0.01);
2262  TF1 *Omega_Func= new TF1("Omega_Func","1/x",Limit_Min,Limit_Max);
2263  Double_t Omega=Omega_Func->GetRandom();
2264 
2265  Double_t Theta=TMath::ASin(kT_Scale/Omega);
2266  //cout<<"angle_omega_kt"<<Theta<<" "<<Omega<<" "<<kT_Scale<<endl;
2267  if(Theta>fJetRadius) continue;
2268 
2269  fastjet::PseudoJet ExtraTrack(kT_Scale/TMath::Sqrt(2),kT_Scale/TMath::Sqrt(2),Omega*TMath::Cos(Theta),Omega);
2270  //boost the particle in the rest frame of the jet to the lab frame
2271  fastjet::PseudoJet ExtraTrack_LabFrame=ExtraTrack.boost(Lab_Jet);
2272  ExtraTrack_LabFrame.set_user_index(i+Jet->GetNumberOfTracks()+100);
2273  fInputVectors.push_back(ExtraTrack_LabFrame);
2274  delete kT_Func;
2275  delete Omega_Func;
2276  }
2277  return fInputVectors;
2278 
2279 }
2280 
2281 
2283 
2284 //______________________________________________________________________________________
2285 Double_t AliAnalysisTaskSubJetFraction::FjNSubJettinessFastJet(std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> Jet_ClusterSequence, Int_t JetContNb,Int_t N, Int_t Algorithm, Double_t Beta, Int_t Option, Double_t Beta_SD, Double_t ZCut){
2286 
2287  //Only Works for kGenOnTheFly
2288 
2289  //WARNING!!! Only works for parent jets that are clustered with Anti-Kt! To change go to AliEmcalJetFinder.cxx and look at the Nsubjettiness() function
2290 
2291  //Algorithm==0 -> kt_axes;
2292  // Algorithm==1 -> ca_axes;
2293  //Algorithm==2 -> antikt_0p2_axes;
2294  //Algorithm==3 -> wta_kt_axes;
2295  //Algorithm==4 -> wta_ca_axes;
2296  //Algorithm==5 -> onepass_kt_axes;
2297  //Algorithm==6 -> onepass_ca_axes;
2298  //Algorithm==7 -> onepass_antikt_0p2_axes;
2299  //Algorithm==8 -> onepass_wta_kt_axes;
2300  //Algorithm==9 -> onepass_wta_ca_axes;
2301  //Algorithm==10 -> min_axes;
2302 
2303 
2304  //fSoftDropOn==1 All the below are done on the jet after it has undergone softdrop with CA. Option 3 and 4 are done whether this is turned on or not.
2305 
2306 
2307  //Option==0 returns Nsubjettiness Value
2308  //Option==1 && N==2 returns opening angle between two subjet axes(Delta R?)
2309  //Option==2 && N==2 returns Delta R
2310  //Option==3 returns first splitting distance for soft dropped jet (CA)
2311  //Option==4 returns Symmetry measure (Zg) for soft dropped jet (CA)
2312  //Option==5 returns Pt of Subjet1
2313  //Option==6 returns Pt of Subjet2
2314  //Options==7 trutns deltaR of subjets...Is this different to before??
2315  //Option==8 Subjet1 Leading track Pt
2316  //Option==9 Subjet1 Leading track Pt
2317 
2318  fastjet::PseudoJet Jet=Jet_ClusterSequence.first;
2319  AliFJWrapper FJWrapper("FJWrapper", "FJWrapper");
2320  FJWrapper.SetR(fJetRadius);
2321  FJWrapper.SetMinJetPt(Jet.perp());
2322  FJWrapper.SetAlgorithm(fastjet::antikt_algorithm); //this is for the jet clustering not the subjet reclustering.
2323  FJWrapper.SetRecombScheme(static_cast<fastjet::RecombinationScheme>(0));
2324  if (Jet.constituents().size()>=N){
2325  Algorithm=fReclusteringAlgorithm;
2326  Double_t dVtx[3]={1,1,1};
2327  return FJWrapper.NSubjettinessDerivativeSub(N,Algorithm,fSubJetRadius,Beta,fJetRadius,Jet,Option,0,Beta_SD,ZCut,fSoftDropOn); //change this
2328  }
2329  else return -2;
2330 }
AliEmcalJetFinder::SetRadius
void SetRadius(Double_t val)
Definition: AliEmcalJetFinder.h:40
AliAnalysisTaskSubJetFraction::fCentSelectOn
Bool_t fCentSelectOn
Random number generator.
Definition: AliAnalysisTaskSubJetFraction.h:166
AliEmcalJet::Area
Double_t Area() const
Definition: AliEmcalJet.h:130
AliPicoTrack.h
Double_t
double Double_t
Definition: External.C:58
XBinsPtSize
Double_t XBinsPtSize
Definition: AliAnalysisTaskSubJetFraction.cxx:54
AliEmcalJetFinder::SetJetMinPt
void SetJetMinPt(Double_t val)
Definition: AliEmcalJetFinder.h:45
AliEmcalJetFinder.h
AliEmcalJetShapeProperties::GetSecondOrderSubtracted1subjettiness_akt02
Double_t GetSecondOrderSubtracted1subjettiness_akt02() const
Definition: AliEmcalJetShapeProperties.h:194
AliAnalysisTaskEmcal::DeltaPhi
static Double_t DeltaPhi(Double_t phia, Double_t phib, Double_t rMin=-TMath::Pi()/2, Double_t rMax=3 *TMath::Pi()/2)
Definition: AliAnalysisTaskEmcal.h:960
TH2F
Definition: External.C:236
AliEmcalJet::GetLeadingTrack
AliVParticle * GetLeadingTrack(TClonesArray *tracks=0) const
Definition: AliEmcalJet.cxx:487
AliEmcalPythiaInfo::GetPartonEta6
Float_t GetPartonEta6() const
Definition: AliEmcalPythiaInfo.h:45
AliEmcalJet::ClosestJet
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:327
AliAnalysisTaskSubJetFraction::Run
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.
Definition: AliAnalysisTaskSubJetFraction.cxx:836
AliAnalysisTaskEmcalJet::GetJetContainer
AliJetContainer * GetJetContainer(Int_t i=0) const
Definition: AliAnalysisTaskEmcalJet.cxx:377
AliEmcalPythiaInfo::GetPartonEta7
Float_t GetPartonEta7() const
Definition: AliEmcalPythiaInfo.h:50
AliAnalysisTaskSubJetFraction::fShapesVar_Tracks_Truth
std::vector< std::vector< Double_t > > fShapesVar_Tracks_Truth
Definition: AliAnalysisTaskSubJetFraction.h:155
AliEmcalJet::Eta
Double_t Eta() const
Definition: AliEmcalJet.h:121
AliEmcalJet::Py
Double_t Py() const
Definition: AliEmcalJet.h:107
AliEmcalJet::Phi
Double_t Phi() const
Definition: AliEmcalJet.h:117
AliTrackContainer
Container with name, TClonesArray and cuts for particles.
Definition: AliTrackContainer.h:56
AliAnalysisTaskSubJetFraction::ModifyJet
std::pair< fastjet::PseudoJet, fastjet::ClusterSequence * > ModifyJet(AliEmcalJet *Jet, Int_t JetContNb, TString Modification)
Definition: AliAnalysisTaskSubJetFraction.cxx:2138
AliEmcalJetShapeProperties::GetFirstOrderSubtractedOpeningAngle_onepassca
Double_t GetFirstOrderSubtractedOpeningAngle_onepassca() const
Definition: AliEmcalJetShapeProperties.h:243
AliAnalysisTaskSubJetFraction::PTD
Double_t PTD(AliEmcalJet *Jet, Int_t JetContNb)
Definition: AliAnalysisTaskSubJetFraction.cxx:1793
AliEmcalJetShapeProperties::GetSecondOrderSubtracted2subjettiness_akt02
Double_t GetSecondOrderSubtracted2subjettiness_akt02() const
Definition: AliEmcalJetShapeProperties.h:204
AliAnalysisTaskSubJetFraction::RetrieveEventObjects
Bool_t RetrieveEventObjects()
Retrieve common objects from event.
Definition: AliAnalysisTaskSubJetFraction.cxx:2085
AliAnalysisTaskEmcal::fEPV0
Double_t fEPV0
!event plane V0
Definition: AliAnalysisTaskEmcal.h:903
AliEmcalJet::GetLabel
Int_t GetLabel() const
Definition: AliEmcalJet.h:124
AliAnalysisTaskSubJetFraction::AddExtraProng
std::vector< fastjet::PseudoJet > AddExtraProng(std::vector< fastjet::PseudoJet > fInputVectors, Double_t Distance, Double_t PtFrac)
Definition: AliAnalysisTaskSubJetFraction.cxx:2216
AliEmcalJetShapeProperties::GetFirstOrderSubtracted2subjettiness_kt
Double_t GetFirstOrderSubtracted2subjettiness_kt() const
Definition: AliEmcalJetShapeProperties.h:133
AliEmcalParticle.h
AliEmcalPythiaInfo::GetPartonPhi7
Float_t GetPartonPhi7() const
Definition: AliEmcalPythiaInfo.h:51
AliFJWrapper::SetMinJetPt
void SetMinJetPt(Double_t MinPt)
Definition: AliFJWrapper.h:142
AliAnalysisTaskSubJetFraction::NSubJettiness
Double_t NSubJettiness(AliEmcalJet *Jet, Int_t JetContNb, Double_t JetRadius, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t A, Int_t B)
Definition: AliAnalysisTaskSubJetFraction.cxx:1926
AliEmcalJetShapeProperties::GetFirstOrderSubtractedOpeningAngle_kt
Double_t GetFirstOrderSubtractedOpeningAngle_kt() const
Definition: AliEmcalJetShapeProperties.h:153
AliFJWrapper::SetRecombScheme
void SetRecombScheme(const fastjet::RecombinationScheme &scheme)
Definition: AliFJWrapper.h:122
AliEmcalJetShapeProperties::GetSecondOrderSubtractedOpeningAngle_akt02
Double_t GetSecondOrderSubtractedOpeningAngle_akt02() const
Definition: AliEmcalJetShapeProperties.h:214
AliAnalysisTaskSubJetFraction::RelativePhi
Double_t RelativePhi(Double_t Phi1, Double_t Phi2)
Definition: AliAnalysisTaskSubJetFraction.cxx:1687
AliEmcalJet::E
Double_t E() const
Definition: AliEmcalJet.h:119
AliAnalysisTaskSubJetFraction::FillHistograms
Bool_t FillHistograms()
Function filling histograms.
Definition: AliAnalysisTaskSubJetFraction.cxx:845
Eta_Bins
Int_t Eta_Bins
Definition: AliAnalysisTaskSubJetFraction.cxx:61
AliEmcalJetShapeProperties::GetFirstOrderSubtracted2subjettiness_ca
Double_t GetFirstOrderSubtracted2subjettiness_ca() const
Definition: AliEmcalJetShapeProperties.h:173
AliEmcalPythiaInfo::GetPartonFlag7
Int_t GetPartonFlag7() const
Definition: AliEmcalPythiaInfo.h:48
AliParticleContainer
Container for particles within the EMCAL framework.
Definition: AliParticleContainer.h:31
XBinsJetMass
Double_t XBinsJetMass[28]
Definition: AliAnalysisTaskSubJetFraction.cxx:57
AliEmcalPythiaInfo::GetPartonPhi6
Float_t GetPartonPhi6() const
Definition: AliEmcalPythiaInfo.h:46
AliEmcalJetShapeProperties::GetFirstOrderSubtractedOpeningAngle_ca
Double_t GetFirstOrderSubtractedOpeningAngle_ca() const
Definition: AliEmcalJetShapeProperties.h:183
AliEmcalJet::TrackAt
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:160
AliEmcalJet::GetNumberOfTracks
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:139
FJ_includes.h
AliFJWrapper::SetR
void SetR(Double_t r)
Definition: AliFJWrapper.h:127
AliEmcalJetShapeProperties::GetFirstOrderSubtracted3subjettiness_kt
Double_t GetFirstOrderSubtracted3subjettiness_kt() const
Definition: AliEmcalJetShapeProperties.h:143
kData
TString kData
Declare data MC or deltaAOD.
Definition: AddTaskCaloTrackCorr.C:63
AliEmcalJet::Px
Double_t Px() const
Definition: AliEmcalJet.h:106
AliEmcalJetShapeProperties::GetSecondOrderSubtracted2subjettiness_kt
Double_t GetSecondOrderSubtracted2subjettiness_kt() const
Definition: AliEmcalJetShapeProperties.h:134
AliEmcalJetFinder::Nsubjettiness
Double_t Nsubjettiness(AliEmcalJet *pJet, AliJetContainer *pContJets, Double_t dVtx[3], Int_t N, Int_t Algorithm, Double_t Radius, Double_t Beta, Int_t Option=0, Int_t Measure=0, Double_t Beta_SD=0, Double_t ZCut=0.1, Int_t SoftDropOn=0)
Definition: AliEmcalJetFinder.cxx:286
AliAnalysisTaskEmcal::GetParticleContainer
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
Definition: AliAnalysisTaskEmcal.cxx:1639
AliEmcalPythiaInfo::GetPartonFlag6
Int_t GetPartonFlag6() const
Definition: AliEmcalPythiaInfo.h:43
AliJetContainer::GetParticleContainer
AliParticleContainer * GetParticleContainer() const
Definition: AliJetContainer.h:186
AliEmcalJetFinder::SetRecombSheme
void SetRecombSheme(Int_t val)
Definition: AliEmcalJetFinder.h:42
AliEmcalJetShapeProperties::GetSecondOrderSubtracted1subjettiness_onepassca
Double_t GetSecondOrderSubtracted1subjettiness_onepassca() const
Definition: AliEmcalJetShapeProperties.h:224
Eta_Up
Double_t Eta_Up
Definition: AliAnalysisTaskSubJetFraction.cxx:59
AliEmcalJetFinder::SetJetAlgorithm
void SetJetAlgorithm(Int_t val)
Definition: AliEmcalJetFinder.h:41
Int_t
int Int_t
Definition: External.C:63
AliAnalysisTaskSubJetFraction::fShapesVar_Tracks_Rec
std::vector< std::vector< Double_t > > fShapesVar_Tracks_Rec
Definition: AliAnalysisTaskSubJetFraction.h:154
AliEmcalJetShapeProperties::GetSecondOrderSubtracted1subjettiness_ca
Double_t GetSecondOrderSubtracted1subjettiness_ca() const
Definition: AliEmcalJetShapeProperties.h:164
Eta_Low
Double_t Eta_Low
Definition: AliAnalysisTaskSubJetFraction.cxx:60
Float_t
float Float_t
Definition: External.C:68
XBinsJetPtSize
Int_t XBinsJetPtSize
Definition: AliAnalysisTaskSubJetFraction.cxx:56
AliEmcalJetShapeProperties::GetFirstOrderSubtracted2subjettiness_akt02
Double_t GetFirstOrderSubtracted2subjettiness_akt02() const
Definition: AliEmcalJetShapeProperties.h:203
AliParticleContainer.h
AliJetContainer::GetNJets
Int_t GetNJets() const
Definition: AliJetContainer.h:155
AliFJWrapper::SetAlgorithm
void SetAlgorithm(const fastjet::JetAlgorithm &algor)
Definition: AliFJWrapper.h:121
AliAnalysisTaskSubJetFraction::SelectTriggerHadron
Int_t SelectTriggerHadron(Float_t PtMin, Float_t PtMax)
Definition: AliAnalysisTaskSubJetFraction.cxx:1701
Phi_Up
Double_t Phi_Up
Definition: AliAnalysisTaskSubJetFraction.cxx:62
TH2D
Definition: External.C:228
TH1D
Definition: External.C:212
AliJetContainer.h
AliParticleContainer::GetAcceptParticle
virtual AliVParticle * GetAcceptParticle(Int_t i=-1) const
Definition: AliParticleContainer.cxx:113
Option
Definition: Option.C:68
AliEmcalJetShapeProperties::GetSecondOrderSubtracted2subjettiness_ca
Double_t GetSecondOrderSubtracted2subjettiness_ca() const
Definition: AliEmcalJetShapeProperties.h:174
AliAnalysisTaskEmcal::fCent
Double_t fCent
!event centrality
Definition: AliAnalysisTaskEmcal.h:901
AliEmcalJetShapeProperties::GetFirstOrderSubtracted2subjettiness_onepassca
Double_t GetFirstOrderSubtracted2subjettiness_onepassca() const
Definition: AliEmcalJetShapeProperties.h:233
XBinsJetMassSize
Int_t XBinsJetMassSize
Definition: AliAnalysisTaskSubJetFraction.cxx:58
AliAnalysisTaskSubJetFraction::SubJetFraction
Double_t SubJetFraction(AliEmcalJet *Jet, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t Type, Bool_t Add, Bool_t Loss)
Definition: AliAnalysisTaskSubJetFraction.cxx:1886
AliRhoParameter.h
AliEmcalJetFinder::SetJetMaxEta
void SetJetMaxEta(Double_t val)
Definition: AliEmcalJetFinder.h:44
Phi_Low
Double_t Phi_Low
Definition: AliAnalysisTaskSubJetFraction.cxx:63
AliEmcalJetShapeProperties::GetSecondOrderSubtractedOpeningAngle_onepassca
Double_t GetSecondOrderSubtractedOpeningAngle_onepassca() const
Definition: AliEmcalJetShapeProperties.h:244
AliJetContainer::GetNextAcceptJet
AliEmcalJet * GetNextAcceptJet()
Definition: AliJetContainer.cxx:308
Data
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)
Definition: charmCutsOptimization.C:371
AliAnalysisTaskSubJetFraction::Recluster
AliEmcalJetFinder * Recluster(AliEmcalJet *Jet, Int_t JetContNb, Double_t SubJetRadius, Double_t SubJetMinPt, Int_t Algorithm, const char *Name)
Definition: AliAnalysisTaskSubJetFraction.cxx:1814
AliEmcalJetFinder::GetJet
AliEmcalJet * GetJet(Int_t index)
Definition: AliEmcalJetFinder.h:37
AliAnalysisTaskSubJetFraction.h
AliAnalysisTaskSubJetFraction::FjNSubJettinessFastJet
Double_t FjNSubJettinessFastJet(std::pair< fastjet::PseudoJet, fastjet::ClusterSequence * > Jet_ClusterSequence, Int_t JetContNb, Int_t N, Int_t Algorithm, Double_t Beta, Int_t Option=0, Double_t Beta_SD=0, Double_t ZCut=0.1)
Definition: AliAnalysisTaskSubJetFraction.cxx:2285
AliEmcalJet::Pt
Double_t Pt() const
Definition: AliEmcalJet.h:109
AliEmcalJetShapeProperties::GetSecondOrderSubtracted2subjettiness_onepassca
Double_t GetSecondOrderSubtracted2subjettiness_onepassca() const
Definition: AliEmcalJetShapeProperties.h:234
AliAnalysisTaskEmcalJet::GetRhoVal
Double_t GetRhoVal(Int_t i=0) const
Definition: AliAnalysisTaskEmcalJet.cxx:626
Track
void Track()
Definition: DrawProductionComparison.C:1494
AliEmcalJetShapeProperties::GetFirstOrderSubtractedOpeningAngle_akt02
Double_t GetFirstOrderSubtractedOpeningAngle_akt02() const
Definition: AliEmcalJetShapeProperties.h:213
AliEmcalJetShapeProperties::GetFirstOrderSubtracted1subjettiness_akt02
Double_t GetFirstOrderSubtracted1subjettiness_akt02() const
Definition: AliEmcalJetShapeProperties.h:193
AliAnalysisTaskSubJetFraction::FjNSubJettiness
Double_t FjNSubJettiness(AliEmcalJet *Jet, Int_t JetContNb, Int_t N, Int_t Algorithm, Double_t Beta, Int_t Option=0, Double_t Beta_SD=0, Double_t ZCut=0.1)
Definition: AliAnalysisTaskSubJetFraction.cxx:1968
XBinsPt
Double_t XBinsPt[66]
Definition: AliAnalysisTaskSubJetFraction.cxx:53
AliAnalysisTaskEmcal::fOutput
AliEmcalList * fOutput
!output list
Definition: AliAnalysisTaskEmcal.h:924
AliAnalysisTaskSubJetFraction::RandomiseTracks
std::vector< fastjet::PseudoJet > RandomiseTracks(AliEmcalJet *Jet, std::vector< fastjet::PseudoJet > fInputVectors)
Definition: AliAnalysisTaskSubJetFraction.cxx:2186
AliEmcalJetShapeProperties::GetFirstOrderSubtracted1subjettiness_onepassca
Double_t GetFirstOrderSubtracted1subjettiness_onepassca() const
Definition: AliEmcalJetShapeProperties.h:223
nVar
const Int_t nVar
Definition: AliAnalysisTaskSubJetFraction.h:22
AliAnalysisTaskEmcal::GetTrackContainer
AliTrackContainer * GetTrackContainer(Int_t i=0) const
Definition: AliAnalysisTaskEmcal.h:223
AliEmcalPythiaInfo
Store some informaion about a Pythia eventThis class is used to store some information about a Pythia...
Definition: AliEmcalPythiaInfo.h:31
AliAnalysisTaskEmcal::SetMakeGeneralHistograms
void SetMakeGeneralHistograms(Bool_t g)
Definition: AliAnalysisTaskEmcal.h:257
AliTrackContainer::GetAcceptTrack
virtual AliVTrack * GetAcceptTrack(Int_t i=-1) const
Definition: AliTrackContainer.cxx:269
AliAnalysisTaskEmcalJet
Base task in the EMCAL jet framework.
Definition: AliAnalysisTaskEmcalJet.h:30
AliEmcalJetShapeProperties::GetFirstOrderSubtracted1subjettiness_ca
Double_t GetFirstOrderSubtracted1subjettiness_ca() const
Definition: AliEmcalJetShapeProperties.h:163
AliEmcalJet
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
AliAnalysisTaskEmcal::GetPythiaInfo
const AliEmcalPythiaInfo * GetPythiaInfo() const
Definition: AliAnalysisTaskEmcal.h:334
AliEmcalJet::Pz
Double_t Pz() const
Definition: AliEmcalJet.h:108
AliEmcalJetShapeProperties::GetFirstOrderSubtracted1subjettiness_kt
Double_t GetFirstOrderSubtracted1subjettiness_kt() const
Definition: AliEmcalJetShapeProperties.h:123
AliEmcalPythiaInfo.h
Declaration of class AliEmcalPythiaInfo.
AliAnalysisTaskSubJetFraction::Angularity
Double_t Angularity(AliEmcalJet *Jet, Int_t JetContNb)
Definition: AliAnalysisTaskSubJetFraction.cxx:1769
Option_t
const char Option_t
Definition: External.C:48
AliFJWrapper::NSubjettinessDerivativeSub
Double32_t NSubjettinessDerivativeSub(Int_t N, Int_t Algorithm, Double_t Radius, Double_t Beta, Double_t JetR, fastjet::PseudoJet jet, Int_t Option=0, Int_t Measure=0, Double_t Beta_SD=0.0, Double_t ZCut=0.1, Int_t SoftDropOn=0)
AliAnalysisTaskSubJetFraction::RelativePhiEventPlane
Double_t RelativePhiEventPlane(Double_t EventPlane, Double_t Phi)
Definition: AliAnalysisTaskSubJetFraction.cxx:1677
AliAnalysisTaskEmcal::UserCreateOutputObjects
void UserCreateOutputObjects()
Main initialization function on the worker.
Definition: AliAnalysisTaskEmcal.cxx:342
AliAnalysisTaskSubJetFraction::AddkTTracks
std::vector< fastjet::PseudoJet > AddkTTracks(AliEmcalJet *Jet, std::vector< fastjet::PseudoJet > fInputVectors, Double_t QHat, Double_t Xlength, Int_t NAdditionalTracks)
Definition: AliAnalysisTaskSubJetFraction.cxx:2239
Bool_t
bool Bool_t
Definition: External.C:53
AliAnalysisTaskSubJetFraction::SubJetOrdering
Double_t SubJetOrdering(AliEmcalJet *Jet, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t Type, Bool_t Index)
Definition: AliAnalysisTaskSubJetFraction.cxx:1839
AliEmcalJet::GetShapeProperties
AliEmcalJetShapeProperties * GetShapeProperties() const
Definition: AliEmcalJet.h:361
AliEmcalJetShapeProperties::GetSecondOrderSubtracted3subjettiness_kt
Double_t GetSecondOrderSubtracted3subjettiness_kt() const
Definition: AliEmcalJetShapeProperties.h:144
XBinsJetPt
Double_t XBinsJetPt[35]
Definition: AliAnalysisTaskSubJetFraction.cxx:55
Phi_Bins
Int_t Phi_Bins
Definition: AliAnalysisTaskSubJetFraction.cxx:64
AliEmcalJet::M
Double_t M() const
Definition: AliEmcalJet.h:120
AliEmcalJetShapeProperties::GetSecondOrderSubtracted1subjettiness_kt
Double_t GetSecondOrderSubtracted1subjettiness_kt() const
Definition: AliEmcalJetShapeProperties.h:124
AliEmcalJet.h
AliJetContainer
Container for jet within the EMCAL jet framework.
Definition: AliJetContainer.h:37
AliEmcalJetShapeProperties::GetSecondOrderSubtractedOpeningAngle_ca
Double_t GetSecondOrderSubtractedOpeningAngle_ca() const
Definition: AliEmcalJetShapeProperties.h:184
AliEmcalJetShapeProperties::GetSecondOrderSubtractedOpeningAngle_kt
Double_t GetSecondOrderSubtractedOpeningAngle_kt() const
Definition: AliEmcalJetShapeProperties.h:154
AliJetContainer::GetJet
AliEmcalJet * GetJet(Int_t i) const
Definition: AliJetContainer.cxx:279