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