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