AliPhysics  f1cc956 (f1cc956)
 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  AliParticleContainer *PartCont =NULL;
713  if (fJetShapeSub==kConstSub) PartCont = GetParticleContainer(1);
714  else PartCont = GetParticleContainer(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  AliParticleContainer *PartCont = NULL;
1274  if (fJetShapeSub==kConstSub) PartCont = GetParticleContainer(1);
1275  else PartCont = GetParticleContainer(0);
1276  TClonesArray *TracksArray = PartCont->GetArray();
1277 
1278  if(!PartCont || !TracksArray) return -99999;
1279  AliAODTrack *Track = 0x0;
1280  AliEmcalParticle *EmcalParticle = 0x0;
1281  Int_t Trigger_Index[100];
1282  for (Int_t i=0; i<100; i++) Trigger_Index[i] = 0;
1283  Int_t Trigger_Counter = 0;
1284  for(Int_t i=0; i <= TracksArray->GetEntriesFast(); i++){
1285  if (fJetShapeSub == kConstSub){
1286  EmcalParticle = static_cast<AliEmcalParticle*>(TracksArray->At(i));
1287  if (!EmcalParticle) continue;
1288  if(TMath::Abs(EmcalParticle->Eta())>0.9) continue;
1289  if (EmcalParticle->Pt()<0.15) continue;
1290  if ((EmcalParticle->Pt() >= PtMin) && (EmcalParticle->Pt()< PtMax)) {
1291  Trigger_Index[Trigger_Counter] = i;
1292  Trigger_Counter++;
1293  }
1294  }
1295  else{
1296  Track = static_cast<AliAODTrack*>(TracksArray->At(i));
1297  if (!Track) continue;
1298  if(TMath::Abs(Track->Eta())>0.9) continue;
1299  if (Track->Pt()<0.15) continue;
1300  if ((Track->Pt() >= PtMin) && (Track->Pt()< PtMax)) {
1301  Trigger_Index[Trigger_Counter] = i;
1302  Trigger_Counter++;
1303  }
1304  }
1305  }
1306  if (Trigger_Counter == 0) return -99999;
1307  Int_t RandomNumber = 0, Index = 0 ;
1308  TRandom3* Random = new TRandom3(0);
1309  RandomNumber = Random->Integer(Trigger_Counter);
1310  Index = Trigger_Index[RandomNumber];
1311  return Index;
1312 }
1313 
1314 
1315 //--------------------------------------------------------------------------
1317 
1318  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1319  Double_t Angularity_Numerator=0; //Reset these values
1320  Double_t Angularity_Denominator=0;
1321  AliVParticle *Particle=0x0;
1322  Double_t DeltaPhi=0;
1323 
1324 
1325  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks (particles in the jet
1326  Particle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
1327  if(!Particle) continue;
1328  DeltaPhi=RelativePhi(Jet->Phi(),Particle->Phi());
1329  Angularity_Numerator=Angularity_Numerator+(Particle->Pt()*TMath::Sqrt(((Particle->Eta()-Jet->Eta())*(Particle->Eta()-Jet->Eta()))+(DeltaPhi*DeltaPhi)));
1330  Angularity_Denominator= Angularity_Denominator+Particle->Pt();
1331  }
1332  if(Angularity_Denominator!=0) return Angularity_Numerator/Angularity_Denominator;
1333  else return -1;
1334 
1335 }
1336 
1337 
1338 
1339 //--------------------------------------------------------------------------
1340 Double_t AliAnalysisTaskSubJetFraction::PTD(AliEmcalJet *Jet, Int_t JetContNb){
1341 
1342  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1343  Double_t PTD_Numerator=0; //Reset these values
1344  Double_t PTD_Denominator=0;
1345  AliVParticle *Particle=0x0;
1346  Double_t DeltaPhi=0;
1347  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks (particles in the jet
1348  Particle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
1349  if(!Particle) continue;
1350  PTD_Numerator=PTD_Numerator+(Particle->Pt()*Particle->Pt());
1351  PTD_Denominator=PTD_Denominator+Particle->Pt();
1352  }
1353  if(PTD_Denominator!=0) return TMath::Sqrt(PTD_Numerator)/PTD_Denominator;
1354  else return -1;
1355 
1356 }
1357 
1358 
1359 //----------------------------------------------------------------------
1361 AliEmcalJetFinder *AliAnalysisTaskSubJetFraction::Recluster(AliEmcalJet *Jet, Int_t JetContNb, Double_t SubJetRadius, Double_t SubJetMinPt, Int_t Algorithm, const char* Name){
1362 
1363  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1364  AliEmcalJetFinder *Reclusterer = new AliEmcalJetFinder(Name); //JetFinder Object for reclustered jets
1365  Reclusterer->SetRadius(SubJetRadius);
1366  Reclusterer->SetJetMinPt(SubJetMinPt);
1367  Reclusterer->SetJetAlgorithm(Algorithm); //0 for anti-kt 1 for kt
1368  Reclusterer->SetJetMaxEta(0.9);
1369  Reclusterer->SetRecombSheme(0);
1371  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
1372  Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
1373  if(Reclusterer->AliEmcalJetFinder::Filter(Jet, JetCont, dVtx)){;} //reclustering jet1 using the jetfinderobject Reclusterer
1374  }
1375  else{
1376  Double_t dVtx[3]={0,0,0};
1377  if(Reclusterer->AliEmcalJetFinder::Filter(Jet, JetCont, dVtx)){;} //reclustering jet1 using the jetfinderobject Reclusterer
1378  }
1379  return Reclusterer;
1380 
1381 }
1382 
1383 
1384 
1385 
1386 
1387 
1388 
1389 
1390 
1391 //----------------------------------------------------------------------
1392 Double_t AliAnalysisTaskSubJetFraction::SubJetOrdering(AliEmcalJet *Jet, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t Type, Bool_t Index){
1393  AliEmcalJet *SubJet=NULL;
1394  Double_t SortingVariable;
1395  Int_t ArraySize =N+1;
1396  TArrayD *JetSorter = new TArrayD(ArraySize);
1397  TArrayD *JetIndexSorter = new TArrayD(ArraySize);
1398  for (Int_t i=0; i<ArraySize; i++){
1399  JetSorter->SetAt(0,i);
1400  }
1401  for (Int_t i=0; i<ArraySize; i++){
1402  JetIndexSorter->SetAt(0,i);
1403  }
1404  if(Reclusterer->GetNumberOfJets()<N) return -999;
1405  for (Int_t i=0; i<Reclusterer->GetNumberOfJets(); i++){
1406  SubJet=Reclusterer->GetJet(i);
1407  if (Type==0) SortingVariable=SubJet->Pt();
1408  else if (Type==1) SortingVariable=SubJet->E();
1409  else if (Type==2) SortingVariable=SubJet->M();
1410  for (Int_t j=0; j<N; j++){
1411  if (SortingVariable>JetSorter->GetAt(j)){
1412  for (Int_t k=N-1; k>=j; k--){
1413  JetSorter->SetAt(JetSorter->GetAt(k),k+1);
1414  JetIndexSorter->SetAt(JetIndexSorter->GetAt(k),k+1);
1415  }
1416  JetSorter->SetAt(SortingVariable,j);
1417  JetIndexSorter->SetAt(i,j);
1418  break;
1419  }
1420  }
1421  }
1422  if (!Index) return JetSorter->GetAt(N-1);
1423  else return JetIndexSorter->GetAt(N-1);
1424 }
1425 
1426 
1427 
1428 //returns -1 if the Nth hardest jet is requested where N>number of available jets
1429 //type: 0=Pt 1=E 2=M
1430 //Index TRUE=returns index FALSE=returns value of quantatiy in question
1431 
1432 
1433 
1434 
1435 
1436 
1437 
1438 //----------------------------------------------------------------------
1439 Double_t AliAnalysisTaskSubJetFraction::SubJetFraction(AliEmcalJet *Jet, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t Type, Bool_t Add, Bool_t Loss){
1440  AliEmcalJet *SubJet=NULL;
1441  Double_t Observable=0;
1442  Double_t Fraction_Numerator=0;
1443  Bool_t Error=kFALSE;
1444  if (!Jet->GetNumberOfTracks()) return -2;
1445  if (Add){
1446  for (Int_t i=1; i<=N; i++){
1447  Observable=SubJetOrdering(Jet,Reclusterer,i,Type,kFALSE);
1448  if(Observable==-999){
1449  Error = kTRUE;
1450  return -2;
1451  }
1452  Fraction_Numerator=Fraction_Numerator+Observable;
1453  }
1454  }
1455  else {
1456  Fraction_Numerator=SubJetOrdering(Jet,Reclusterer,N,Type,kFALSE);
1457  if (Fraction_Numerator==-999) return -2;
1458  }
1459  if (Type==0){
1460  if(Loss) return (Jet->Pt()-Fraction_Numerator)/Jet->Pt();
1461  else return Fraction_Numerator/Jet->Pt();
1462  }
1463  else if (Type==1){
1464  if(Loss) return (Jet->E()-Fraction_Numerator)/Jet->E();
1465  else return Fraction_Numerator/Jet->E();
1466  }
1467  else { //change to else if if you want to add more later
1468  if(Loss) return (Jet->M()-Fraction_Numerator)/Jet->M();
1469  else return Fraction_Numerator/Jet->M();
1470  }
1471 }
1472 //N number of hardest subjets involved
1473 //Type 0=Pt 1=E 2=M
1474 // Add TRUE: Add 1 to Nth hardest subjets together/total jet False: Nth hardest Jet/toal jet
1475 //Loss TRUE: Jet-Subjet(s)/Jet FALSE: Subjet(s)/Jet
1476 
1477 
1478 //----------------------------------------------------------------------
1479 Double_t AliAnalysisTaskSubJetFraction::NSubJettiness(AliEmcalJet *Jet, Int_t JetContNb, Double_t JetRadius, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t A, Int_t B){
1480  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1481  AliEmcalJet *SubJet=NULL;
1482  Double_t DeltaR1=0;
1483  Double_t DeltaR2=0;
1484  AliVParticle *JetParticle=0x0;
1485  Double_t SubJetiness_Numerator = 0;
1486  Double_t SubJetiness_Denominator = 0;
1487  Double_t Index=-2;
1488  Bool_t Error=kFALSE;
1489  // JetRadius=TMath::Sqrt((Jet->Area()/TMath::Pi())); //comment out later
1490  if (!Jet->GetNumberOfTracks()) return -2;
1491  if (Reclusterer->GetNumberOfJets() < N) return -2;
1492  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks (particles in the jet
1493  JetParticle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
1494  for (Int_t j=1; j<=N; j++){
1495  Index=SubJetOrdering(Jet,Reclusterer,j,0,kTRUE);
1496  if(Index==-999){
1497  Error = kTRUE;
1498  i=Jet->GetNumberOfTracks();
1499  break;
1500  }
1501  if(j==1){
1502  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);
1503  }
1504  else{
1505  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);
1506  if (DeltaR2<DeltaR1) DeltaR1=DeltaR2;
1507  }
1508  }
1509  SubJetiness_Numerator=SubJetiness_Numerator+(JetParticle->Pt()*DeltaR1);
1510  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));
1511  else SubJetiness_Denominator=SubJetiness_Denominator+(TMath::Power((Reclusterer->GetJet(SubJetOrdering(Jet,Reclusterer,N,0,kTRUE))->Pt()),A)*JetParticle->Pt()*TMath::Power(JetRadius,B));
1512  }
1513  if (SubJetiness_Denominator!=0 && !Error){
1514  // if (SubJetiness_Numerator/SubJetiness_Denominator>1.0) cout << JetRadius<<endl;
1515  return SubJetiness_Numerator/SubJetiness_Denominator; }
1516  else return -2;
1517 }
1518 
1519 
1520 //______________________________________________________________________________________
1521 Double_t AliAnalysisTaskSubJetFraction::fjNSubJettiness(AliEmcalJet *Jet, Int_t JetContNb,Int_t N, Int_t Algorithm, Double_t Beta, Int_t Option){
1522 
1523  //WARNING!!! Only works for parent jets that are clustered with Anti-Kt! To change go to AliEmcalJetFinder.cxx and look at the Nsubjettiness() function
1524 
1525  //Algorithm==0 -> kt_axes;
1526  // Algorithm==1 -> ca_axes;
1527  //Algorithm==2 -> antikt_0p2_axes;
1528  //Algorithm==3 -> wta_kt_axes;
1529  //Algorithm==4 -> wta_ca_axes;
1530  //Algorithm==5 -> onepass_kt_axes;
1531  //Algorithm==6 -> onepass_ca_axes;
1532  //Algorithm==7 -> onepass_antikt_0p2_axes;
1533  //Algorithm==8 -> onepass_wta_kt_axes;
1534  //Algorithm==9 -> onepass_wta_ca_axes;
1535  //Algorithm==10 -> min_axes;
1536 
1537 
1538  //Option==0 returns Nsubjettiness Value
1539  //Option==1 && N==2 returns opening angle between two subjet axes(Delta R?)
1540  //Option==2 && N==2 returns Delta R
1541  if (Jet->GetNumberOfTracks()>=N){
1542  if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==1) && (Algorithm==0) && (Beta==1.0) && (Option==0)){
1545  }
1546  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==0) && (Beta==1.0) && (Option==0)){
1549  }
1550  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==3) && (Algorithm==0) && (Beta==1.0) && (Option==0)){
1553  }
1554  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==0) && (Beta==1.0) && (Option==1)){
1557  }
1558  else{
1559  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1560  AliEmcalJetFinder *JetFinder=new AliEmcalJetFinder("Nsubjettiness");
1561  JetFinder->SetJetMaxEta(0.9-fJetRadius);
1562  JetFinder->SetRadius(fJetRadius);
1563  JetFinder->SetJetAlgorithm(0); //0 for anti-kt 1 for kt //this is for the JET!!!!!!!!!! Not the SubJets
1564  JetFinder->SetRecombSheme(0);
1565  JetFinder->SetJetMinPt(Jet->Pt());
1567  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
1568  Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
1569  return JetFinder->Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,fSubJetRadius,Beta,Option);
1570  }
1571  else{
1572  Double_t dVtx[3]={0,0,0};
1573  return JetFinder->Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,fSubJetRadius,Beta,Option);
1574  }
1575  }
1576  }
1577  else return -2;
1578 }
1579 
1580 
1581 //________________________________________________________________________
1583  //
1584  // retrieve event objects
1585  //
1587  return kFALSE;
1588 
1589  return kTRUE;
1590 }
1591 
1592 
1593 //_______________________________________________________________________
1595 {
1596  // Called once at the end of the analysis.
1598 
1599  /*
1600  for (int i=1; i<=(XBinsJetPtSize)-1; i++){ //Rescales the fhJetPt Histograms according to the with of the variable bins and number of events
1601  fhJetPt->SetBinContent(i,((1.0/(XBinsJetPt[i]-XBinsJetPt[i-1]))*((fhJetPt->GetBinContent(i))*(1.0/(fhEventCounter->GetBinContent(1))))));
1602  fhSubJetPt->SetBinContent(i,((1.0/(XBinsJetPt[i]-XBinsJetPt[i-1]))*((fhSubJetPt->GetBinContent(i))*(1.0/(fhEventCounter->GetBinContent(5))))));
1603 
1604  //fhJetPt->SetBinContent(i,((1.0/(fhPt->GetBinWidth(i)))*((fhJetPt->GetBinContent(i))*(fhEventCounter->GetBinContent(1)))));
1605  // fhJetPt->SetBinContent(i,((1.0/((fhPt->GetLowEdge(i+1))-(fhJetPt->GetLowEdge(i))))*((fhPt->GetBinContent(i))*(fhEventCounter->GetBinContent(1)))));
1606  }
1607  for (int i=1; i<=(XBinsJetMassSize)-1; i++){ //Rescales the fhPt Histograms according to the with of the variable bins and number of events
1608  fhJetMass->SetBinContent(i,((1.0/(XBinsJetMass[i]-XBinsJetMass[i-1]))*((fhJetMass->GetBinContent(i))*(1.0/(fhEventCounter->GetBinContent(1))))));
1609  }
1610 
1611  fhJetPhi->Scale(Phi_Bins/((Phi_Up-Phi_Low)*((fhEventCounter->GetBinContent(1)))));
1612  fhJetEta->Scale(Eta_Bins/((Eta_Up-Eta_Low)*((fhEventCounter->GetBinContent(1)))));
1613  fhJetRadius->Scale(100/(fhEventCounter->GetBinContent(4))); //should this and JetAngularity be divided by Bin 1 or 4????
1614  fhNumberOfJetTracks->Scale(1.0/(fhEventCounter->GetBinContent(4)));
1615 
1616  fhJetCounter->Scale(1.0/(fhEventCounter->GetBinContent(1))); //is the first bin the correct one to look at?
1617 
1618 
1619  fhSubJetCounter->Scale(1.0/(fhEventCounter->GetBinContent(5)));
1620 
1621  */
1622  /*
1623 
1624  fhJetPt->Scale(1.0/(fhEventCounter->GetBinContent(1)));
1625  fhSubJetPt->Scale(1.0/(fhEventCounter->GetBinContent(5)));
1626  fhJetMass->Scale(1.0/(fhEventCounter->GetBinContent(1)));
1627  fhJetPhi->Scale(1.0/(fhEventCounter->GetBinContent(1)));
1628  fhJetEta->Scale(1.0/(fhEventCounter->GetBinContent(1)));
1629  fhJetRadius->Scale(1.0/(fhEventCounter->GetBinContent(4))); //should this and JetAngularity be divided by Bin 1 or 4????
1630  fhNumberOfJetTracks->Scale(1.0/(fhEventCounter->GetBinContent(4)));
1631  fhJetCounter->Scale(1.0/(fhEventCounter->GetBinContent(1))); //is the first bin the correct one to look at?
1632  fhSubJetCounter->Scale(1.0/(fhEventCounter->GetBinContent(5)));
1633  */
1634  }
1635 
1636 }
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
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
Container for particles within the EMCAL framework.
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(Int_t i=0) const
AliParticleContainer * GetParticleContainer() const
void SetRecombSheme(Int_t val)
void SetJetAlgorithm(Int_t val)
Double_t Eta() const
Double_t Pt() const
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
void SetMakeGeneralHistograms(Bool_t g)
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