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