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