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