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