AliPhysics  695988a (695988a)
 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
713  TClonesArray *TrackArray = PartCont->GetArray();
714  TriggerHadron = static_cast<AliAODTrack*>(TrackArray->At(TriggerHadronLabel));
715  if (!TriggerHadron) return 0;//No trigger hadron with label found
716  if(fSemigoodCorrect){
717  Double_t HoleDistance=RelativePhi(TriggerHadron->Phi(),fHolePos);
718  if(TMath::Abs(HoleDistance)+fHoleWidth>TMath::Pi()-fRecoilAngularWindow){
719  return 0;}
720  }
721  fhPtTriggerHadron->Fill(TriggerHadron->Pt()); //Needed for per trigger Normalisation
722  }
723 
724 
725 
728  AliEmcalJet *Jet1 = NULL; //Subtracted hybrid Jet
729  AliEmcalJet *Jet2 = NULL; //Unsubtracted Hybrid Jet
730  AliEmcalJet *Jet3 = NULL; //Detector Level Pythia Jet
731  AliEmcalJet *Jet4 = NULL; //Particle Level Pyhtia Jet
732  AliJetContainer *JetCont1= GetJetContainer(0); //Jet Container for Subtracted Hybrid Jets
733  AliJetContainer *JetCont2= GetJetContainer(1); //Jet Container for Unsubtracted Hybrid Jets
734  AliJetContainer *JetCont3= GetJetContainer(2); //Jet Container for Detector Level Pyhtia Jets
735  AliJetContainer *JetCont4= GetJetContainer(3); //Jet Container for Particle Level Pythia Jets
736  AliEmcalJetFinder *Reclusterer1; //Object containg Subjets from Subtracted Hybrid Jets
737  AliEmcalJetFinder *Reclusterer4; //Object containg Subjets from Particle Level Pythia Jets
738  Bool_t JetsMatched=kFALSE;
739  Bool_t EventCounter=kFALSE;
740  Int_t JetNumber=-1;
741  Double_t JetPtThreshold=-2;
742  fhEventCounter->Fill(1);
743  if(JetCont1) {
744  fhEventCounter->Fill(2);
745  JetCont1->ResetCurrentID();
746  while((Jet1=JetCont1->GetNextAcceptJet())) {
747  if (fJetShapeSub==kConstSub) JetPtThreshold=Jet1->Pt();
748  else JetPtThreshold=Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
749  if ( (!Jet1) || (JetPtThreshold<fPtThreshold)) continue;
750  else {
751  if(fSemigoodCorrect){
752  Double_t HoleDistance=RelativePhi(Jet1->Phi(),fHolePos);
753  if(TMath::Abs(HoleDistance)<fHoleWidth) continue;
754  }
755  Float_t RecoilDeltaPhi = 0.;
756  if (fJetSelection == kRecoil){
757  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), Jet1->Phi());
758  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
759  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), Jet1->Pt());
760  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), Jet1->Phi()));
761  }
762  if (!EventCounter){
763  fhEventCounter->Fill(3);
764  EventCounter=kTRUE;
765  }
766  if(fJetShapeSub==kConstSub){
767  JetNumber=-1;
768  for(Int_t i = 0; i<JetCont2->GetNJets(); i++) {
769  Jet2 = JetCont2->GetJet(i);
770  if(Jet2->GetLabel()==Jet1->GetLabel()) {
771  JetNumber=i;
772  }
773  }
774  if(JetNumber==-1) continue;
775  Jet2=JetCont2->GetJet(JetNumber);
776  if (JetCont2->AliJetContainer::GetFractionSharedPt(Jet2)<fSharedFractionPtMin) continue;
777  Jet3=Jet2->ClosestJet();
778  }
779  if(!(fJetShapeSub==kConstSub)){
780  if (!(JetCont1->AliJetContainer::GetFractionSharedPt(Jet1)<fSharedFractionPtMin)) continue;
781  Jet3 = Jet1->ClosestJet();
782  }
783  if (!Jet3) continue;
784  Jet4=Jet3->ClosestJet();
785  if(!Jet4) continue;
786  JetsMatched=kTRUE;
788  if (fJetShapeSub==kConstSub) fShapesVar[0]=Jet1->Pt();
789  else fShapesVar[0]=Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
790  fShapesVar[2]=fjNSubJettiness(Jet1,0,1,0,1,0);
791  fShapesVar[4]=fjNSubJettiness(Jet1,0,2,0,1,0);
792  fShapesVar[6]=fjNSubJettiness(Jet1,0,3,0,1,0);
793  fShapesVar[8]=fjNSubJettiness(Jet1,0,2,0,1,1);
794  fShapesVar[10]=Jet1->GetNumberOfTracks();
795  if (fFullTree){
796  fShapesVar[12]=fjNSubJettiness(Jet1,0,2,0,1,2);
797  Reclusterer1 = Recluster(Jet1, 0, fSubJetRadius, fSubJetMinPt, fSubJetAlgorithm, "SubJetFinder_1");
798  fShapesVar[14]=SubJetFraction(Jet1, Reclusterer1, 1, 0, kTRUE, kFALSE);
799  fShapesVar[16]=SubJetFraction(Jet1, Reclusterer1, 2, 0, kTRUE, kFALSE);
800  }
801  if (JetsMatched){ //even needed? Not now but might be if you want to fill trees when jets aren't matched too
802  fShapesVar[1]=Jet4->Pt();
803  fShapesVar[3]=fjNSubJettiness(Jet4,3,1,0,1,0);
804  fShapesVar[5]=fjNSubJettiness(Jet4,3,2,0,1,0);
805  fShapesVar[7]=fjNSubJettiness(Jet4,3,3,0,1,0);
806  fShapesVar[9]=fjNSubJettiness(Jet4,3,2,0,1,1);
807  fShapesVar[11]=Jet4->GetNumberOfTracks();
808  if (fFullTree){
809  fShapesVar[13]=fjNSubJettiness(Jet4,3,2,0,1,2);
810  Reclusterer4=Recluster(Jet4, 3, fSubJetRadius, 0, fSubJetAlgorithm, "SubJetFinder_4");
811  fShapesVar[15]=SubJetFraction(Jet4, Reclusterer4, 1, 0, kTRUE, kFALSE);
812  fShapesVar[17]=SubJetFraction(Jet4, Reclusterer4, 2, 0, kTRUE, kFALSE);
813  }
814  }
815  else{
816  fShapesVar[1]=-2;
817  fShapesVar[3]=-2;
818  fShapesVar[5]=-2;
819  fShapesVar[7]=-2;
820  fShapesVar[9]=-2;
821  fShapesVar[11]=-2;
822  if (fFullTree){
823  fShapesVar[13]=-2;
824  fShapesVar[15]=-2;
825  fShapesVar[17]=-2;
826  }
827  }
828  fTreeResponseMatrixAxis->Fill();
829  JetsMatched=kFALSE;
830  }
831  }
832  }
833  }
834 
837  AliEmcalJet *Jet1 = NULL; //Detector Level Jet
838  AliEmcalJet *Jet2 = NULL; //Particle Level Jet
839  AliJetContainer *JetCont1= GetJetContainer(0); //Jet Container for Detector Level Pythia
840  AliJetContainer *JetCont2= GetJetContainer(1); //Jet Container for Particle Level Pythia
841  AliEmcalJetFinder *Reclusterer1; //Object containg Subjets Detector Level
842  AliEmcalJetFinder *Reclusterer2; //Object containg Subjets Particle Level
843  Int_t JetCounter1=0; //Counts number of jets in event
844  Int_t JetCounter2=0; //Counts number of jets in event
845  Double_t JetPhi1=0;
846  Double_t JetPhi2=0;
847  Bool_t JetsMatched=kFALSE;
848  Double_t Pythia_Event_Weight=1;
849  Bool_t EventCounter=kFALSE;
850  fhEventCounter_1->Fill(1);
851  if(JetCont1) {
852  fhEventCounter_1->Fill(2); //Number of events with a jet container
853  JetCont1->ResetCurrentID();
854  while((Jet1=JetCont1->GetNextAcceptJet())) {
855  if( (!Jet1) || ((Jet1->Pt())<fPtThreshold)) {
856  // fhEventCounter_1->Fill(3); //events where the jet had a problem
857  continue;
858  }
859  else {
860  if(fSemigoodCorrect){
861  Double_t HoleDistance=RelativePhi(Jet1->Phi(),fHolePos);
862  if(TMath::Abs(HoleDistance)<fHoleWidth) continue;
863  }
864  Float_t RecoilDeltaPhi = 0.;
865  if (fJetSelection == kRecoil){
866  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), Jet1->Phi());
867  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
868  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), Jet1->Pt());
869  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), Jet1->Phi()));
870  }
871  if (!EventCounter){
872  fhEventCounter_1->Fill(3);
873  EventCounter=kTRUE;
874  }
875  if((Jet1->GetNumberOfTracks())==0){
876  fhEventCounter_1->Fill(10); //zero track jets
877  }
878  if((Jet1->GetNumberOfTracks())==1){
879  fhEventCounter_1->Fill(11); //one track jets
880  }
881  fhEventCounter_1->Fill(4); //Number of Jets found in all events
882  JetCounter1++;
883  fhJetPt_1->Fill(Jet1->Pt());
884  JetPhi1=Jet1->Phi();
885  if(JetPhi1 < -1*TMath::Pi()) JetPhi1 += (2*TMath::Pi());
886  else if (JetPhi1 > TMath::Pi()) JetPhi1 -= (2*TMath::Pi());
887  fhJetPhi_1->Fill(JetPhi1);
888  fhJetEta_1->Fill(Jet1->Eta());
889  fhJetMass_1->Fill(Jet1->M());
890  fhJetRadius_1->Fill(TMath::Sqrt((Jet1->Area()/TMath::Pi()))); //Radius of Jets per event
892  if((Jet2 = Jet1->ClosestJet())){
893  JetsMatched=kTRUE;
894  if((Jet2->GetNumberOfTracks())==0){
895  fhEventCounter_2->Fill(10); //zero track jets
896  }
897  if((Jet2->GetNumberOfTracks())==1){
898  fhEventCounter_2->Fill(11); //one track jets
899  }
900  fhEventCounter_2->Fill(4); //Number of Jets found in all events
901  JetCounter2++;
902  fhJetPt_2->Fill(Jet2->Pt());
903  JetPhi2=Jet2->Phi();
904  if(JetPhi2 < -1*TMath::Pi()) JetPhi2 += (2*TMath::Pi());
905  else if (JetPhi2 > TMath::Pi()) JetPhi2 -= (2*TMath::Pi());
906  fhJetPhi_2->Fill(JetPhi2);
907  fhJetEta_2->Fill(Jet2->Eta());
908  fhJetMass_2->Fill(Jet2->M());
909  fhJetRadius_2->Fill(TMath::Sqrt((Jet2->Area()/TMath::Pi()))); //Radius of Jets per event
911  fh2PtRatio->Fill(Jet1->Pt(),Jet2->Pt(),Pythia_Event_Weight);
912  }
913  else {
914  fhEventCounter_2->Fill(1);
915  //continue;
916  }
918 
919 
920  fShapesVar[0]=Jet1->Pt();
921  fShapesVar[2]=fjNSubJettiness(Jet1,0,1,0,1,0);
922  fShapesVar[4]=fjNSubJettiness(Jet1,0,2,0,1,0);
923  fShapesVar[6]=fjNSubJettiness(Jet1,0,3,0,1,0);
924  fShapesVar[8]=fjNSubJettiness(Jet1,0,2,0,1,1);
925  fShapesVar[10]=Jet1->GetNumberOfTracks();
926  Reclusterer1 = Recluster(Jet1, 0, fSubJetRadius, fSubJetMinPt, fSubJetAlgorithm, "SubJetFinder_1");
927  if (fFullTree){
928  fShapesVar[12]=fjNSubJettiness(Jet1,0,2,0,1,2);
929  fShapesVar[14]=SubJetFraction(Jet1, Reclusterer1, 1, 0, kTRUE, kFALSE);
930  fShapesVar[16]=SubJetFraction(Jet1, Reclusterer1, 2, 0, kTRUE, kFALSE);
931  }
932  if (JetsMatched){ //even needed? Not now but might be if you want to fill trees when jets aren't matched too
933  fShapesVar[1]=Jet2->Pt();
934  fShapesVar[3]=fjNSubJettiness(Jet2,1,1,0,1,0);
935  fShapesVar[5]=fjNSubJettiness(Jet2,1,2,0,1,0);
936  fShapesVar[7]=fjNSubJettiness(Jet2,1,3,0,1,0);
937  fShapesVar[9]=fjNSubJettiness(Jet2,1,2,0,1,1);
938  fShapesVar[11]=Jet2->GetNumberOfTracks();
939  Reclusterer2 = Recluster(Jet2, 1, fSubJetRadius, 0, fSubJetAlgorithm, "SubJetFinder_2");
940  if (fFullTree){
941  fShapesVar[13]=fjNSubJettiness(Jet2,1,2,0,1,2);
942  fShapesVar[15]=SubJetFraction(Jet2, Reclusterer2, 1, 0, kTRUE, kFALSE);
943  fShapesVar[17]=SubJetFraction(Jet2, Reclusterer2, 2, 0, kTRUE, kFALSE);
944  }
945  }
946  else{
947  fShapesVar[1]=-2;
948  fShapesVar[3]=-2;
949  fShapesVar[5]=-2;
950  fShapesVar[7]=-2;
951  fShapesVar[9]=-2;
952  fShapesVar[11]=-2;
953  if (fFullTree){
954  fShapesVar[13]=-2;
955  fShapesVar[15]=-2;
956  fShapesVar[17]=-2;
957  }
958  }
959  fTreeResponseMatrixAxis->Fill();
960 
961  fhSubJetCounter_1->Fill(Reclusterer1->GetNumberOfJets());
962  for (Int_t i= 0; i<Reclusterer1->GetNumberOfJets(); i++){
963  fhEventCounter_1->Fill(6); //Number of overall subjets in all events
964  fhSubJetPt_1->Fill(Reclusterer1->GetJet(i)->Pt());
965  fhSubJetMass_1->Fill(Reclusterer1->GetJet(i)->M());
966  fhNumberOfSubJetTracks_1->Fill(Reclusterer1->GetJet(i)->GetNumberOfTracks());
967  fhSubJetRadius_1->Fill(TMath::Sqrt((Reclusterer1->GetJet(i)->Area()/TMath::Pi()))); //Radius of SubJets per event
968  }
969  if(JetsMatched){
970  fhSubJetCounter_2->Fill(Reclusterer2->GetNumberOfJets());
971  for (Int_t i= 0; i<Reclusterer2->GetNumberOfJets(); i++){
972  fhEventCounter_2->Fill(6); //Number of overall subjets in all events
973  fhSubJetPt_2->Fill(Reclusterer2->GetJet(i)->Pt());
974  fhSubJetMass_2->Fill(Reclusterer2->GetJet(i)->M());
975  fhNumberOfSubJetTracks_2->Fill(Reclusterer2->GetJet(i)->GetNumberOfTracks());
976  fhSubJetRadius_2->Fill(TMath::Sqrt((Reclusterer2->GetJet(i)->Area()/TMath::Pi()))); //Radius of SubJets per event
977  }
978  }
979  JetsMatched=kFALSE;
980  }
981  }
982  fhJetCounter_1->Fill(JetCounter1); //Number of Jets in Each Event Particle Level
983  fhJetCounter_2->Fill(JetCounter2); //Number of Jets in Each Event Detector Level
984  }
985  //else {fhEventCounter_1->Fill(3);} //Events with no jet container
986  }
987 
988 
989 
991  AliEmcalJet *Jet1 = NULL; //Original Jet in the event
992  AliJetContainer *JetCont= GetJetContainer(0); //Jet Container for event
993  Int_t JetCounter=0; //Counts number of jets in event
994  Double_t JetPhi=0;
995  Bool_t EventCounter=kFALSE;
996  Double_t JetPt_ForThreshold=0;
997  fhEventCounter->Fill(1);
998  if(JetCont) {
999  fhEventCounter->Fill(2); //Number of events with a jet container
1000  JetCont->ResetCurrentID();
1001  while((Jet1=JetCont->GetNextAcceptJet())) {
1002  if(!Jet1) continue;
1003  if (fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) JetPt_ForThreshold = Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
1004  else JetPt_ForThreshold = Jet1->Pt();
1005  if(JetPt_ForThreshold<fPtThreshold) {
1006  //fhEventCounter->Fill(3); //events where the jet had a problem
1007  continue;
1008  }
1009  else {
1010  if(fSemigoodCorrect){
1011  Double_t HoleDistance=RelativePhi(Jet1->Phi(),fHolePos);
1012  if(TMath::Abs(HoleDistance)<fHoleWidth) continue;
1013  }
1014  Float_t RecoilDeltaPhi = 0.;
1015  if (fJetSelection == kRecoil){
1016  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), Jet1->Phi());
1017  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
1018  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), Jet1->Pt());
1019  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), Jet1->Phi()));
1020  }
1021  if (!EventCounter){
1022  fhEventCounter->Fill(3);
1023  EventCounter=kTRUE;
1024  }
1025  if((Jet1->GetNumberOfTracks())==0){
1026  fhEventCounter->Fill(10); //zero track jets
1027  }
1028  if((Jet1->GetNumberOfTracks())==1){
1029  fhEventCounter->Fill(11); //one track jets
1030  }
1031  fhEventCounter->Fill(4); //Number of Jets found in all events
1032  JetCounter++;
1033  fhJetPt->Fill(Jet1->Pt());
1034  JetPhi=Jet1->Phi();
1035  if(JetPhi < -1*TMath::Pi()) JetPhi += (2*TMath::Pi());
1036  else if (JetPhi > TMath::Pi()) JetPhi -= (2*TMath::Pi());
1037  fhJetPhi->Fill(JetPhi);
1038  fhJetEta->Fill(Jet1->Eta());
1039  fhJetMass->Fill(Jet1->M());
1040  fhJetRadius->Fill(TMath::Sqrt((Jet1->Area()/TMath::Pi()))); //Radius of Jets per event
1041  fhNumberOfJetTracks->Fill(Jet1->GetNumberOfTracks());
1042  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) fShapesVar[0]= Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
1043  else fShapesVar[0]=Jet1->Pt();
1044  fShapesVar[2]=fjNSubJettiness(Jet1,0,1,0,1,0);
1045  fShapesVar[4]=fjNSubJettiness(Jet1,0,2,0,1,0);
1046  fShapesVar[6]=fjNSubJettiness(Jet1,0,3,0,1,0);
1047  fShapesVar[8]=fjNSubJettiness(Jet1,0,2,0,1,1);
1048  fShapesVar[10]=Jet1->GetNumberOfTracks();
1049  AliEmcalJetFinder *Reclusterer1 = Recluster(Jet1, 0, fSubJetRadius, fSubJetMinPt, fSubJetAlgorithm, "SubJetFinder");
1050  if (fFullTree){
1051  fShapesVar[12]=fjNSubJettiness(Jet1,0,2,0,1,2);
1052  fShapesVar[14]=SubJetFraction(Jet1, Reclusterer1, 1, 0, kTRUE, kFALSE);
1053  fShapesVar[16]=SubJetFraction(Jet1, Reclusterer1, 2, 0, kTRUE, kFALSE);
1054  }
1055  fShapesVar[1]=-2;
1056  fShapesVar[3]=-2;
1057  fShapesVar[5]=-2;
1058  fShapesVar[7]=-2;
1059  fShapesVar[9]=-2;
1060  fShapesVar[11]=-2;
1061  if (fFullTree){
1062  fShapesVar[13]=-2;
1063  fShapesVar[15]=-2;
1064  fShapesVar[17]=-2;
1065  }
1066  fTreeResponseMatrixAxis->Fill();
1067  fhSubJetCounter->Fill(Reclusterer1->GetNumberOfJets());
1068  for (Int_t i= 0; i<Reclusterer1->GetNumberOfJets(); i++){
1069  fhEventCounter->Fill(6); //Number of overall subjets in all events
1070  fhSubJetPt->Fill(Reclusterer1->GetJet(i)->Pt());
1071  fhSubJetMass->Fill(Reclusterer1->GetJet(i)->M());
1072  fhNumberOfSubJetTracks->Fill(Reclusterer1->GetJet(i)->GetNumberOfTracks());
1073  fhSubJetRadius->Fill(TMath::Sqrt((Reclusterer1->GetJet(i)->Area()/TMath::Pi()))); //Radius of SubJets per event
1074  }
1075  }
1076  }
1077  fhJetCounter->Fill(JetCounter); //Number of Jets in Each Event
1078  }
1079  //else {fhEventCounter->Fill(2);} //Events with no jet container
1080  }
1081 
1082 
1084 
1085  AliEmcalJet *Jet1 = NULL; //Original Jet in the event
1086  AliJetContainer *JetCont= GetJetContainer(0); //Jet Container for event
1087  Int_t JetCounter=0; //Counts number of jets in event
1088  Double_t JetPhi=0;
1089  Bool_t EventCounter=kFALSE;
1090  Double_t JetPt_ForThreshold=0;
1091  fhEventCounter->Fill(1);
1092  if(JetCont) {
1093  fhEventCounter->Fill(2); //Number of events with a jet container
1094  JetCont->ResetCurrentID();
1095  while((Jet1=JetCont->GetNextAcceptJet())) {
1096  // if (((Jet1=JetCont->GetLeadingJet("rho")) && (fPtThreshold<=Jet1->Pt()-(GetRhoVal(0)*Jet1->Area())))){
1097  if(!Jet1) continue;
1098  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) JetPt_ForThreshold = Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
1099  else JetPt_ForThreshold = Jet1->Pt();
1100  if(JetPt_ForThreshold<fPtThreshold) {
1101  //fhEventCounter->Fill(3); //events where the jet had a problem
1102  continue;
1103  }
1104  else {
1105  if(fSemigoodCorrect){
1106  Double_t HoleDistance=RelativePhi(Jet1->Phi(),fHolePos);
1107  if(TMath::Abs(HoleDistance)<fHoleWidth) continue;
1108  }
1109  Float_t RecoilDeltaPhi = 0.;
1110  if (fJetSelection == kRecoil){
1111  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), Jet1->Phi());
1112  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
1113  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), Jet1->Pt());
1114  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), Jet1->Phi()));
1115  }
1116  if (!EventCounter){
1117  fhEventCounter->Fill(3);
1118  EventCounter=kTRUE;
1119  }
1120  if((Jet1->GetNumberOfTracks())==0){
1121  fhEventCounter->Fill(10); //zero track jets
1122  }
1123  if((Jet1->GetNumberOfTracks())==1){
1124  fhEventCounter->Fill(11); //one track jets
1125  }
1126  fhEventCounter->Fill(4); //Number of Jets found in all events
1127  JetCounter++;
1128  fhJetPt->Fill(Jet1->Pt());
1129  JetPhi=Jet1->Phi();
1130  if(JetPhi < -1*TMath::Pi()) JetPhi += (2*TMath::Pi());
1131  else if (JetPhi > TMath::Pi()) JetPhi -= (2*TMath::Pi());
1132  fhJetPhi->Fill(JetPhi);
1133  fhJetEta->Fill(Jet1->Eta());
1134  fhJetMass->Fill(Jet1->M());
1135  fhJetRadius->Fill(TMath::Sqrt((Jet1->Area()/TMath::Pi()))); //Radius of Jets per event
1136  fhNumberOfJetTracks->Fill(Jet1->GetNumberOfTracks());
1137 
1138  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
1139  Double_t FJ_Beta=1.0;
1140  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));
1141  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));
1142  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));
1143  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));
1144  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));
1145  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));
1146  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));
1147  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));
1148  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));
1149  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));
1150  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));
1151  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));
1152  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));
1153  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));
1154  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));
1155  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));
1156  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));
1157  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));
1158  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));
1159  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));
1160  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));
1161  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));
1162 
1163  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));
1164  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));
1165  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));
1166  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));
1167  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));
1168  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));
1169  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));
1170  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));
1171  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));
1172  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));
1173  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));
1174 
1175 
1176  fhSubJettiness1_FJ_KT->Fill(fjNSubJettiness(Jet1,0,1,0,FJ_Beta,0));
1177  fhSubJettiness1_FJ_CA->Fill(fjNSubJettiness(Jet1,0,1,1,FJ_Beta,0));
1178  fhSubJettiness1_FJ_AKT->Fill(fjNSubJettiness(Jet1,0,1,2,FJ_Beta,0));
1179  fhSubJettiness1_FJ_WTA_KT->Fill(fjNSubJettiness(Jet1,0,1,3,FJ_Beta,0));
1180  fhSubJettiness1_FJ_WTA_CA->Fill(fjNSubJettiness(Jet1,0,1,4,FJ_Beta,0));
1181  fhSubJettiness1_FJ_OP_KT->Fill(fjNSubJettiness(Jet1,0,1,5,FJ_Beta,0));
1182  fhSubJettiness1_FJ_OP_CA->Fill(fjNSubJettiness(Jet1,0,1,6,FJ_Beta,0));
1183  fhSubJettiness1_FJ_OP_AKT->Fill(fjNSubJettiness(Jet1,0,1,7,FJ_Beta,0));
1184  fhSubJettiness1_FJ_OP_WTA_KT->Fill(fjNSubJettiness(Jet1,0,1,8,FJ_Beta,0));
1185  fhSubJettiness1_FJ_OP_WTA_CA->Fill(fjNSubJettiness(Jet1,0,1,9,FJ_Beta,0));
1186  fhSubJettiness1_FJ_MIN->Fill(fjNSubJettiness(Jet1,0,1,10,FJ_Beta,0));
1187  fhSubJettiness2_FJ_KT->Fill(fjNSubJettiness(Jet1,0,2,0,FJ_Beta,0));
1188  fhSubJettiness2_FJ_CA->Fill(fjNSubJettiness(Jet1,0,2,1,FJ_Beta,0));
1189  fhSubJettiness2_FJ_AKT->Fill(fjNSubJettiness(Jet1,0,2,2,FJ_Beta,0)); //because in this case we aren't garantueed to get 2 subjets
1190  fhSubJettiness2_FJ_WTA_KT->Fill(fjNSubJettiness(Jet1,0,2,3,FJ_Beta,0));
1191  fhSubJettiness2_FJ_WTA_CA->Fill(fjNSubJettiness(Jet1,0,2,4,FJ_Beta,0));
1192  fhSubJettiness2_FJ_OP_KT->Fill(fjNSubJettiness(Jet1,0,2,5,FJ_Beta,0));
1193  fhSubJettiness2_FJ_OP_CA->Fill(fjNSubJettiness(Jet1,0,2,6,FJ_Beta,0));
1194  fhSubJettiness2_FJ_OP_AKT->Fill(fjNSubJettiness(Jet1,0,2,7,FJ_Beta,0));
1195  fhSubJettiness2_FJ_OP_WTA_KT->Fill(fjNSubJettiness(Jet1,0,2,8,FJ_Beta,0));
1196  fhSubJettiness2_FJ_OP_WTA_CA->Fill(fjNSubJettiness(Jet1,0,2,9,FJ_Beta,0));
1197  fhSubJettiness2_FJ_MIN->Fill(fjNSubJettiness(Jet1,0,2,10,FJ_Beta,0));
1198 
1199  fhSubJettiness2to1_FJ_KT->Fill(fjNSubJettiness(Jet1,0,2,0,FJ_Beta,0)/fjNSubJettiness(Jet1,0,1,0,FJ_Beta,0));
1200  fhSubJettiness2to1_FJ_CA->Fill(fjNSubJettiness(Jet1,0,2,1,FJ_Beta,0)/fjNSubJettiness(Jet1,0,2,1,FJ_Beta,0));
1201  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
1202  fhSubJettiness2to1_FJ_WTA_KT->Fill(fjNSubJettiness(Jet1,0,2,3,FJ_Beta,0)/fjNSubJettiness(Jet1,0,2,3,FJ_Beta,0));
1203  fhSubJettiness2to1_FJ_WTA_CA->Fill(fjNSubJettiness(Jet1,0,2,4,FJ_Beta,0)/fjNSubJettiness(Jet1,0,2,4,FJ_Beta,0));
1204  fhSubJettiness2to1_FJ_OP_KT->Fill(fjNSubJettiness(Jet1,0,2,5,FJ_Beta,0)/fjNSubJettiness(Jet1,0,2,5,FJ_Beta,0));
1205  fhSubJettiness2to1_FJ_OP_CA->Fill(fjNSubJettiness(Jet1,0,2,6,FJ_Beta,0)/fjNSubJettiness(Jet1,0,2,6,FJ_Beta,0));
1206  fhSubJettiness2to1_FJ_OP_AKT->Fill(fjNSubJettiness(Jet1,0,2,7,FJ_Beta,0)/fjNSubJettiness(Jet1,0,2,7,FJ_Beta,0));
1207  fhSubJettiness2to1_FJ_OP_WTA_KT->Fill(fjNSubJettiness(Jet1,0,2,8,FJ_Beta,0)/fjNSubJettiness(Jet1,0,2,8,FJ_Beta,0));
1208  fhSubJettiness2to1_FJ_OP_WTA_CA->Fill(fjNSubJettiness(Jet1,0,2,9,FJ_Beta,0)/fjNSubJettiness(Jet1,0,2,9,FJ_Beta,0));
1209  fhSubJettiness2to1_FJ_MIN->Fill(fjNSubJettiness(Jet1,0,2,10,FJ_Beta,0)/fjNSubJettiness(Jet1,0,2,10,FJ_Beta,0));
1210 
1211  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) fShapesVar[0]= Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
1212  else fShapesVar[0]=Jet1->Pt();
1213  fShapesVar[2]=fjNSubJettiness(Jet1,0,1,0,1,0);
1214  fShapesVar[4]=fjNSubJettiness(Jet1,0,2,0,1,0);
1215  fShapesVar[6]=fjNSubJettiness(Jet1,0,3,0,1,0);
1216  fShapesVar[8]=fjNSubJettiness(Jet1,0,2,0,1,1);
1217  fShapesVar[10]=Jet1->GetNumberOfTracks();
1218  AliEmcalJetFinder *Reclusterer1 = Recluster(Jet1, 0, fSubJetRadius, fSubJetMinPt, fSubJetAlgorithm, "SubJetFinder");
1219  if (fFullTree){
1220  fShapesVar[12]=fjNSubJettiness(Jet1,0,2,0,1,2);
1221  fShapesVar[14]=SubJetFraction(Jet1, Reclusterer1, 1, 0, kTRUE, kFALSE);
1222  fShapesVar[16]=SubJetFraction(Jet1, Reclusterer1, 2, 0, kTRUE, kFALSE);
1223  }
1224  fShapesVar[1]=-2;
1225  fShapesVar[3]=-2;
1226  fShapesVar[5]=-2;
1227  fShapesVar[7]=-2;
1228  fShapesVar[9]=-2;
1229  fShapesVar[11]=-2;
1230  if (fFullTree){
1231  fShapesVar[13]=-2;
1232  fShapesVar[15]=-2;
1233  fShapesVar[17]=-2;
1234  }
1235  fTreeResponseMatrixAxis->Fill();
1236 
1237  fhSubJetCounter->Fill(Reclusterer1->GetNumberOfJets());
1238  for (Int_t i= 0; i<Reclusterer1->GetNumberOfJets(); i++){
1239  fhEventCounter->Fill(6); //Number of overall subjets in all events
1240  fhSubJetPt->Fill(Reclusterer1->GetJet(i)->Pt());
1241  fhSubJetMass->Fill(Reclusterer1->GetJet(i)->M());
1242  fhNumberOfSubJetTracks->Fill(Reclusterer1->GetJet(i)->GetNumberOfTracks());
1243  fhSubJetRadius->Fill(TMath::Sqrt((Reclusterer1->GetJet(i)->Area()/TMath::Pi()))); //Radius of SubJets per event
1244  }
1245  }
1246  }
1247  fhJetCounter->Fill(JetCounter); //Number of Jets in Each Event
1248  }
1249  //else {fhEventCounter->Fill(2);} //Events with no jet container
1250  }
1251  return kTRUE;
1252 }
1253 
1254 //________________________________________________________________________
1255 Double_t AliAnalysisTaskSubJetFraction::RelativePhi(Double_t Phi1, Double_t Phi2){
1256 
1257  if(Phi1 < -1*TMath::Pi()) Phi1 += (2*TMath::Pi()); // Turns the range of 0to2Pi into -PitoPi ???????????
1258  else if (Phi1 > TMath::Pi()) Phi1 -= (2*TMath::Pi());
1259  if(Phi2 < -1*TMath::Pi()) Phi2 += (2*TMath::Pi());
1260  else if (Phi2 > TMath::Pi()) Phi2 -= (2*TMath::Pi());
1261  Double_t DeltaPhi=Phi2-Phi1;
1262  if(DeltaPhi < -1*TMath::Pi()) DeltaPhi += (2*TMath::Pi());
1263  else if (DeltaPhi > TMath::Pi()) DeltaPhi -= (2*TMath::Pi());
1264  return DeltaPhi;
1265 }
1266 
1267 
1268 //--------------------------------------------------------------------------
1269 Int_t AliAnalysisTaskSubJetFraction::SelectTriggerHadron(Float_t PtMin, Float_t PtMax){
1270 
1272  TClonesArray *TracksArray = PartCont->GetArray();
1273 
1274  if(!PartCont || !TracksArray) return -99999;
1275  AliAODTrack *Track = 0x0;
1276  AliEmcalParticle *EmcalParticle = 0x0;
1277  TList *TrackList = new TList();
1278  Int_t Trigger_Index[100];
1279  for (Int_t i=0; i<100; i++) Trigger_Index[i] = 0;
1280  Int_t Trigger_Counter = 0;
1281  for(Int_t i=0; i <= TracksArray->GetEntriesFast(); i++){
1282  if (fJetShapeSub == kConstSub){
1283  EmcalParticle = static_cast<AliEmcalParticle*>(TracksArray->At(i));
1284  if (!EmcalParticle) continue;
1285  if(TMath::Abs(EmcalParticle->Eta())>0.9) continue;
1286  if (EmcalParticle->Pt()<0.15) continue;
1287  if ((EmcalParticle->Pt() >= PtMin) && (EmcalParticle->Pt()< PtMax)) {
1288  TrackList->Add(EmcalParticle);
1289  Trigger_Index[Trigger_Counter] = i;
1290  Trigger_Counter++;
1291  }
1292  }
1293  else{
1294  Track = static_cast<AliAODTrack*>(TracksArray->At(i));
1295  if (!Track) continue;
1296  if(TMath::Abs(Track->Eta())>0.9) continue;
1297  if (Track->Pt()<0.15) continue;
1298  if ((Track->Pt() >= PtMin) && (Track->Pt()< PtMax)) {
1299  TrackList->Add(Track);
1300  Trigger_Index[Trigger_Counter] = i;
1301  Trigger_Counter++;
1302  }
1303  }
1304  }
1305  if (Trigger_Counter == 0) return -99999;
1306  Int_t RandomNumber = 0, Index = 0 ;
1307  TRandom3* Random = new TRandom3(0);
1308  RandomNumber = Random->Integer(Trigger_Counter);
1309  Index = Trigger_Index[RandomNumber];
1310  return Index;
1311 }
1312 
1313 
1314 //--------------------------------------------------------------------------
1316 
1317  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1318  Double_t Angularity_Numerator=0; //Reset these values
1319  Double_t Angularity_Denominator=0;
1320  AliVParticle *Particle=0x0;
1321  Double_t DeltaPhi=0;
1322 
1323 
1324  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks (particles in the jet
1325  Particle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
1326  if(!Particle) continue;
1327  DeltaPhi=RelativePhi(Jet->Phi(),Particle->Phi());
1328  Angularity_Numerator=Angularity_Numerator+(Particle->Pt()*TMath::Sqrt(((Particle->Eta()-Jet->Eta())*(Particle->Eta()-Jet->Eta()))+(DeltaPhi*DeltaPhi)));
1329  Angularity_Denominator= Angularity_Denominator+Particle->Pt();
1330  }
1331  if(Angularity_Denominator!=0) return Angularity_Numerator/Angularity_Denominator;
1332  else return -1;
1333 
1334 }
1335 
1336 
1337 
1338 //--------------------------------------------------------------------------
1339 Double_t AliAnalysisTaskSubJetFraction::PTD(AliEmcalJet *Jet, Int_t JetContNb){
1340 
1341  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1342  Double_t PTD_Numerator=0; //Reset these values
1343  Double_t PTD_Denominator=0;
1344  AliVParticle *Particle=0x0;
1345  Double_t DeltaPhi=0;
1346  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks (particles in the jet
1347  Particle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
1348  if(!Particle) continue;
1349  PTD_Numerator=PTD_Numerator+(Particle->Pt()*Particle->Pt());
1350  PTD_Denominator=PTD_Denominator+Particle->Pt();
1351  }
1352  if(PTD_Denominator!=0) return TMath::Sqrt(PTD_Numerator)/PTD_Denominator;
1353  else return -1;
1354 
1355 }
1356 
1357 
1358 //----------------------------------------------------------------------
1360 AliEmcalJetFinder *AliAnalysisTaskSubJetFraction::Recluster(AliEmcalJet *Jet, Int_t JetContNb, Double_t SubJetRadius, Double_t SubJetMinPt, Int_t Algorithm, const char* Name){
1361 
1362  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1363  AliEmcalJetFinder *Reclusterer = new AliEmcalJetFinder(Name); //JetFinder Object for reclustered jets
1364  Reclusterer->SetRadius(SubJetRadius);
1365  Reclusterer->SetJetMinPt(SubJetMinPt);
1366  Reclusterer->SetJetAlgorithm(Algorithm); //0 for anti-kt 1 for kt
1367  Reclusterer->SetJetMaxEta(0.9);
1368  Reclusterer->SetRecombSheme(0);
1370  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
1371  Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
1372  if(Reclusterer->AliEmcalJetFinder::Filter(Jet, JetCont, dVtx)){;} //reclustering jet1 using the jetfinderobject Reclusterer
1373  }
1374  else{
1375  Double_t dVtx[3]={0,0,0};
1376  if(Reclusterer->AliEmcalJetFinder::Filter(Jet, JetCont, dVtx)){;} //reclustering jet1 using the jetfinderobject Reclusterer
1377  }
1378  return Reclusterer;
1379 
1380 }
1381 
1382 
1383 
1384 
1385 
1386 
1387 
1388 
1389 
1390 //----------------------------------------------------------------------
1391 Double_t AliAnalysisTaskSubJetFraction::SubJetOrdering(AliEmcalJet *Jet, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t Type, Bool_t Index){
1392  AliEmcalJet *SubJet=NULL;
1393  Double_t SortingVariable;
1394  Int_t ArraySize =N+1;
1395  TArrayD *JetSorter = new TArrayD(ArraySize);
1396  TArrayD *JetIndexSorter = new TArrayD(ArraySize);
1397  for (Int_t i=0; i<ArraySize; i++){
1398  JetSorter->SetAt(0,i);
1399  }
1400  for (Int_t i=0; i<ArraySize; i++){
1401  JetIndexSorter->SetAt(0,i);
1402  }
1403  if(Reclusterer->GetNumberOfJets()<N) return -999;
1404  for (Int_t i=0; i<Reclusterer->GetNumberOfJets(); i++){
1405  SubJet=Reclusterer->GetJet(i);
1406  if (Type==0) SortingVariable=SubJet->Pt();
1407  else if (Type==1) SortingVariable=SubJet->E();
1408  else if (Type==2) SortingVariable=SubJet->M();
1409  for (Int_t j=0; j<N; j++){
1410  if (SortingVariable>JetSorter->GetAt(j)){
1411  for (Int_t k=N-1; k>=j; k--){
1412  JetSorter->SetAt(JetSorter->GetAt(k),k+1);
1413  JetIndexSorter->SetAt(JetIndexSorter->GetAt(k),k+1);
1414  }
1415  JetSorter->SetAt(SortingVariable,j);
1416  JetIndexSorter->SetAt(i,j);
1417  break;
1418  }
1419  }
1420  }
1421  if (!Index) return JetSorter->GetAt(N-1);
1422  else return JetIndexSorter->GetAt(N-1);
1423 }
1424 
1425 
1426 
1427 //returns -1 if the Nth hardest jet is requested where N>number of available jets
1428 //type: 0=Pt 1=E 2=M
1429 //Index TRUE=returns index FALSE=returns value of quantatiy in question
1430 
1431 
1432 
1433 
1434 
1435 
1436 
1437 //----------------------------------------------------------------------
1438 Double_t AliAnalysisTaskSubJetFraction::SubJetFraction(AliEmcalJet *Jet, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t Type, Bool_t Add, Bool_t Loss){
1439  AliEmcalJet *SubJet=NULL;
1440  Double_t Observable=0;
1441  Double_t Fraction_Numerator=0;
1442  Bool_t Error=kFALSE;
1443  if (!Jet->GetNumberOfTracks()) return -2;
1444  if (Add){
1445  for (Int_t i=1; i<=N; i++){
1446  Observable=SubJetOrdering(Jet,Reclusterer,i,Type,kFALSE);
1447  if(Observable==-999){
1448  Error = kTRUE;
1449  return -2;
1450  }
1451  Fraction_Numerator=Fraction_Numerator+Observable;
1452  }
1453  }
1454  else {
1455  Fraction_Numerator=SubJetOrdering(Jet,Reclusterer,N,Type,kFALSE);
1456  if (Fraction_Numerator==-999) return -2;
1457  }
1458  if (Type==0){
1459  if(Loss) return (Jet->Pt()-Fraction_Numerator)/Jet->Pt();
1460  else return Fraction_Numerator/Jet->Pt();
1461  }
1462  else if (Type==1){
1463  if(Loss) return (Jet->E()-Fraction_Numerator)/Jet->E();
1464  else return Fraction_Numerator/Jet->E();
1465  }
1466  else { //change to else if if you want to add more later
1467  if(Loss) return (Jet->M()-Fraction_Numerator)/Jet->M();
1468  else return Fraction_Numerator/Jet->M();
1469  }
1470 }
1471 //N number of hardest subjets involved
1472 //Type 0=Pt 1=E 2=M
1473 // Add TRUE: Add 1 to Nth hardest subjets together/total jet False: Nth hardest Jet/toal jet
1474 //Loss TRUE: Jet-Subjet(s)/Jet FALSE: Subjet(s)/Jet
1475 
1476 
1477 //----------------------------------------------------------------------
1478 Double_t AliAnalysisTaskSubJetFraction::NSubJettiness(AliEmcalJet *Jet, Int_t JetContNb, Double_t JetRadius, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t A, Int_t B){
1479  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1480  AliEmcalJet *SubJet=NULL;
1481  Double_t DeltaR1=0;
1482  Double_t DeltaR2=0;
1483  AliVParticle *JetParticle=0x0;
1484  Double_t SubJetiness_Numerator = 0;
1485  Double_t SubJetiness_Denominator = 0;
1486  Double_t Index=-2;
1487  Bool_t Error=kFALSE;
1488  // JetRadius=TMath::Sqrt((Jet->Area()/TMath::Pi())); //comment out later
1489  if (!Jet->GetNumberOfTracks()) return -2;
1490  if (Reclusterer->GetNumberOfJets() < N) return -2;
1491  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks (particles in the jet
1492  JetParticle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
1493  for (Int_t j=1; j<=N; j++){
1494  Index=SubJetOrdering(Jet,Reclusterer,j,0,kTRUE);
1495  if(Index==-999){
1496  Error = kTRUE;
1497  i=Jet->GetNumberOfTracks();
1498  break;
1499  }
1500  if(j==1){
1501  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);
1502  }
1503  else{
1504  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);
1505  if (DeltaR2<DeltaR1) DeltaR1=DeltaR2;
1506  }
1507  }
1508  SubJetiness_Numerator=SubJetiness_Numerator+(JetParticle->Pt()*DeltaR1);
1509  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));
1510  else SubJetiness_Denominator=SubJetiness_Denominator+(TMath::Power((Reclusterer->GetJet(SubJetOrdering(Jet,Reclusterer,N,0,kTRUE))->Pt()),A)*JetParticle->Pt()*TMath::Power(JetRadius,B));
1511  }
1512  if (SubJetiness_Denominator!=0 && !Error){
1513  // if (SubJetiness_Numerator/SubJetiness_Denominator>1.0) cout << JetRadius<<endl;
1514  return SubJetiness_Numerator/SubJetiness_Denominator; }
1515  else return -2;
1516 }
1517 
1518 
1519 //______________________________________________________________________________________
1520 Double_t AliAnalysisTaskSubJetFraction::fjNSubJettiness(AliEmcalJet *Jet, Int_t JetContNb,Int_t N, Int_t Algorithm, Double_t Beta, Int_t Option){
1521 
1522  //WARNING!!! Only works for parent jets that are clustered with Anti-Kt! To change go to AliEmcalJetFinder.cxx and look at the Nsubjettiness() function
1523 
1524  //Algorithm==0 -> kt_axes;
1525  // Algorithm==1 -> ca_axes;
1526  //Algorithm==2 -> antikt_0p2_axes;
1527  //Algorithm==3 -> wta_kt_axes;
1528  //Algorithm==4 -> wta_ca_axes;
1529  //Algorithm==5 -> onepass_kt_axes;
1530  //Algorithm==6 -> onepass_ca_axes;
1531  //Algorithm==7 -> onepass_antikt_0p2_axes;
1532  //Algorithm==8 -> onepass_wta_kt_axes;
1533  //Algorithm==9 -> onepass_wta_ca_axes;
1534  //Algorithm==10 -> min_axes;
1535 
1536 
1537  //Option==0 returns Nsubjettiness Value
1538  //Option==1 && N==2 returns opening angle between two subjet axes(Delta R?)
1539  //Option==2 && N==2 returns Delta R
1540  if (Jet->GetNumberOfTracks()>=N){
1541  if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==1) && (Algorithm==0) && (Beta==1.0) && (Option==0)){
1544  }
1545  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==0) && (Beta==1.0) && (Option==0)){
1548  }
1549  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==3) && (Algorithm==0) && (Beta==1.0) && (Option==0)){
1552  }
1553  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==0) && (Beta==1.0) && (Option==1)){
1556  }
1557  else{
1558  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1559  AliEmcalJetFinder *JetFinder=new AliEmcalJetFinder("Nsubjettiness");
1560  JetFinder->SetJetMaxEta(0.9-fJetRadius);
1561  JetFinder->SetRadius(fJetRadius);
1562  JetFinder->SetJetAlgorithm(0); //0 for anti-kt 1 for kt //this is for the JET!!!!!!!!!! Not the SubJets
1563  JetFinder->SetRecombSheme(0);
1564  JetFinder->SetJetMinPt(Jet->Pt());
1566  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
1567  Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
1568  return JetFinder->Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,fSubJetRadius,Beta,Option);
1569  }
1570  else{
1571  Double_t dVtx[3]={0,0,0};
1572  return JetFinder->Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,fSubJetRadius,Beta,Option);
1573  }
1574  }
1575  }
1576  else return -2;
1577 }
1578 
1579 
1580 //________________________________________________________________________
1582  //
1583  // retrieve event objects
1584  //
1586  return kFALSE;
1587 
1588  return kTRUE;
1589 }
1590 
1591 
1592 //_______________________________________________________________________
1594 {
1595  // Called once at the end of the analysis.
1597 
1598  /*
1599  for (int i=1; i<=(XBinsJetPtSize)-1; i++){ //Rescales the fhJetPt Histograms according to the with of the variable bins and number of events
1600  fhJetPt->SetBinContent(i,((1.0/(XBinsJetPt[i]-XBinsJetPt[i-1]))*((fhJetPt->GetBinContent(i))*(1.0/(fhEventCounter->GetBinContent(1))))));
1601  fhSubJetPt->SetBinContent(i,((1.0/(XBinsJetPt[i]-XBinsJetPt[i-1]))*((fhSubJetPt->GetBinContent(i))*(1.0/(fhEventCounter->GetBinContent(5))))));
1602 
1603  //fhJetPt->SetBinContent(i,((1.0/(fhPt->GetBinWidth(i)))*((fhJetPt->GetBinContent(i))*(fhEventCounter->GetBinContent(1)))));
1604  // fhJetPt->SetBinContent(i,((1.0/((fhPt->GetLowEdge(i+1))-(fhJetPt->GetLowEdge(i))))*((fhPt->GetBinContent(i))*(fhEventCounter->GetBinContent(1)))));
1605  }
1606  for (int i=1; i<=(XBinsJetMassSize)-1; i++){ //Rescales the fhPt Histograms according to the with of the variable bins and number of events
1607  fhJetMass->SetBinContent(i,((1.0/(XBinsJetMass[i]-XBinsJetMass[i-1]))*((fhJetMass->GetBinContent(i))*(1.0/(fhEventCounter->GetBinContent(1))))));
1608  }
1609 
1610  fhJetPhi->Scale(Phi_Bins/((Phi_Up-Phi_Low)*((fhEventCounter->GetBinContent(1)))));
1611  fhJetEta->Scale(Eta_Bins/((Eta_Up-Eta_Low)*((fhEventCounter->GetBinContent(1)))));
1612  fhJetRadius->Scale(100/(fhEventCounter->GetBinContent(4))); //should this and JetAngularity be divided by Bin 1 or 4????
1613  fhNumberOfJetTracks->Scale(1.0/(fhEventCounter->GetBinContent(4)));
1614 
1615  fhJetCounter->Scale(1.0/(fhEventCounter->GetBinContent(1))); //is the first bin the correct one to look at?
1616 
1617 
1618  fhSubJetCounter->Scale(1.0/(fhEventCounter->GetBinContent(5)));
1619 
1620  */
1621  /*
1622 
1623  fhJetPt->Scale(1.0/(fhEventCounter->GetBinContent(1)));
1624  fhSubJetPt->Scale(1.0/(fhEventCounter->GetBinContent(5)));
1625  fhJetMass->Scale(1.0/(fhEventCounter->GetBinContent(1)));
1626  fhJetPhi->Scale(1.0/(fhEventCounter->GetBinContent(1)));
1627  fhJetEta->Scale(1.0/(fhEventCounter->GetBinContent(1)));
1628  fhJetRadius->Scale(1.0/(fhEventCounter->GetBinContent(4))); //should this and JetAngularity be divided by Bin 1 or 4????
1629  fhNumberOfJetTracks->Scale(1.0/(fhEventCounter->GetBinContent(4)));
1630  fhJetCounter->Scale(1.0/(fhEventCounter->GetBinContent(1))); //is the first bin the correct one to look at?
1631  fhSubJetCounter->Scale(1.0/(fhEventCounter->GetBinContent(5)));
1632  */
1633  }
1634 
1635 }
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