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