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