AliPhysics  7140ed4 (7140ed4)
 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 "TMatrixD.h"
21 #include "TMatrixDSym.h"
22 #include "TMatrixDSymEigen.h"
23 #include "TVector3.h"
24 #include "TVector2.h"
25 #include "AliVCluster.h"
26 #include "AliVTrack.h"
27 #include "AliEmcalJet.h"
28 #include "AliRhoParameter.h"
29 #include "AliLog.h"
30 #include "AliEmcalParticle.h"
31 #include "AliMCEvent.h"
32 #include "AliGenPythiaEventHeader.h"
33 #include "AliAODMCHeader.h"
34 #include "AliMCEvent.h"
35 #include "AliAnalysisManager.h"
36 #include "AliJetContainer.h"
37 #include "AliParticleContainer.h"
38 //#include "AliPythiaInfo.h"
39 #include "TRandom3.h"
40 #include "AliPicoTrack.h"
41 #include "AliEmcalJetFinder.h"
42 #include "AliAODEvent.h"
43 #include "AliAnalysisTaskSubJetFraction.h"
44 
45 //Globals
46 
47 
48 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
49 Double_t XBinsPtSize=66;
50 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
51 Int_t XBinsJetPtSize=35;
52 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
53 Int_t XBinsJetMassSize=28;
54 Double_t Eta_Up=1.00;
55 Double_t Eta_Low=-1.00;
56 Int_t Eta_Bins=100;
57 Double_t Phi_Up=2*(TMath::Pi());
58 Double_t Phi_Low=(-1*(TMath::Pi()));
59 Int_t Phi_Bins=100;
60 
61 
62 
63 
64 
65 using std::cout;
66 using std::endl;
67 
68 ClassImp(AliAnalysisTaskSubJetFraction)
69 
70 //________________________________________________________________________
71 AliAnalysisTaskSubJetFraction::AliAnalysisTaskSubJetFraction() :
72  AliAnalysisTaskEmcalJet("AliAnalysisTaskSubJetFraction", kTRUE),
73  fContainer(0),
74  fMinFractionShared(0),
75  fJetShapeType(kData),
76  fJetShapeSub(kNoSub),
77  fJetSelection(kInclusive),
78  fShapesVar(0),
79  fPtThreshold(-9999.),
80  fRMatching(0.2),
81  fminpTTrig(20.),
82  fmaxpTTrig(50.),
83  fangWindowRecoil(0.6),
84  fSemigoodCorrect(0),
85  fHolePos(0),
86  fHoleWidth(0),
87  fCentSelectOn(kTRUE),
88  fCentMin(0),
89  fCentMax(10),
90  fSubJetAlgorithm(0),
91  fSubJetRadius(0.1),
92  fSubJetMinPt(1),
93  fJetRadius(0.4),
94  fRMatched(0.2),
95  fSharedFractionPtMin(0.5),
96  fDerivSubtrOrder(0),
97  fFullTree(kFALSE),
98  fhJetPt(0x0),
99  fhJetPt_1(0x0),
100  fhJetPt_2(0x0),
101  fhJetPhi(0x0),
102  fhJetPhi_1(0x0),
103  fhJetPhi_2(0x0),
104  fhJetEta(0x0),
105  fhJetEta_1(0x0),
106  fhJetEta_2(0x0),
107  fhJetMass(0x0),
108  fhJetMass_1(0x0),
109  fhJetMass_2(0x0),
110  fhJetRadius(0x0),
111  fhJetRadius_1(0x0),
112  fhJetRadius_2(0x0),
113  fhJetAngularity(0x0),
114  fhJetAngularityJetPt(0x0),
115  fhJetPTD(0x0),
116  fhJetPTDJetPt(0x0),
117  fhJetCounter(0x0),
118  fhJetCounter_1(0x0),
119  fhJetCounter_2(0x0),
120  fhNumberOfJetTracks(0x0),
121  fhNumberOfJetTracks_1(0x0),
122  fhNumberOfJetTracks_2(0x0),
123  fhSubJetPt(0x0),
124  fhSubJetPt_1(0x0),
125  fhSubJetPt_2(0x0),
126  fhSubJetMass(0x0),
127  fhSubJetMass_1(0x0),
128  fhSubJetMass_2(0x0),
129  fhSubJetRadius(0x0),
130  fhSubJetRadius_1(0x0),
131  fhSubJetRadius_2(0x0),
132  fhSubJetCounter(0x0),
133  fhSubJetCounter_1(0x0),
134  fhSubJetCounter_2(0x0),
135  fhNumberOfSubJetTracks(0x0),
136  fhNumberOfSubJetTracks_1(0x0),
137  fhNumberOfSubJetTracks_2(0x0),
138  fhSubJetPtFrac(0x0),
139  fhSubJetPtFrac2(0x0),
140  fhSubJetPtLoss(0x0),
141  fhSubJetPtLoss2(0x0),
142  fhSubJetEnergyFrac(0x0),
143  fhSubJetEnergyFrac2(0x0),
144  fhSubJetEnergyLoss(0x0),
145  fhSubJetEnergyLoss2(0x0),
146  fh2PtRatio(0x0),
147  fhSubJetiness1(0x0),
148  fhSubJetiness1JetPt(0x0),
149  fhSubJetiness2(0x0),
150  fhSubJetiness2JetPt(0x0),
151  fh2to1SubJetinessRatio(0x0),
152  fh2to1SubJetinessRatioJetPt(0x0),
153  fhEventCounter(0x0),
154  fhEventCounter_1(0x0),
155  fhEventCounter_2(0x0),
156  fhJetPtJetEta(0x0),
157  fhSubJetiness2Distance(0x0),
158  fh2JetTracksEtaPhiPt(0x0),
159  fhSubJettiness1CheckRatio_FJ_AKT(0x0),
160  fhSubJettiness1CheckRatio_FJ_KT(0x0),
161  fhSubJettiness1CheckRatio_FJ_CA(0x0),
162  fhSubJettiness1CheckRatio_FJ_WTA_KT(0x0),
163  fhSubJettiness1CheckRatio_FJ_WTA_CA(0x0),
164  fhSubJettiness1CheckRatio_FJ_OP_AKT(0x0),
165  fhSubJettiness1CheckRatio_FJ_OP_KT(0x0),
166  fhSubJettiness1CheckRatio_FJ_OP_CA(0x0),
167  fhSubJettiness1CheckRatio_FJ_OP_WTA_KT(0x0),
168  fhSubJettiness1CheckRatio_FJ_OP_WTA_CA(0x0),
169  fhSubJettiness1CheckRatio_FJ_MIN(0x0),
170  fhSubJettiness2CheckRatio_FJ_AKT(0x0),
171  fhSubJettiness2CheckRatio_FJ_KT(0x0),
172  fhSubJettiness2CheckRatio_FJ_CA(0x0),
173  fhSubJettiness2CheckRatio_FJ_WTA_KT(0x0),
174  fhSubJettiness2CheckRatio_FJ_WTA_CA(0x0),
175  fhSubJettiness2CheckRatio_FJ_OP_AKT(0x0),
176  fhSubJettiness2CheckRatio_FJ_OP_KT(0x0),
177  fhSubJettiness2CheckRatio_FJ_OP_CA(0x0),
178  fhSubJettiness2CheckRatio_FJ_OP_WTA_KT(0x0),
179  fhSubJettiness2CheckRatio_FJ_OP_WTA_CA(0x0),
180  fhSubJettiness2CheckRatio_FJ_MIN(0x0),
181  fhSubJettiness1_FJ_AKT(0x0),
182  fhSubJettiness1_FJ_KT(0x0),
183  fhSubJettiness1_FJ_CA(0x0),
184  fhSubJettiness1_FJ_WTA_KT(0x0),
185  fhSubJettiness1_FJ_WTA_CA(0x0),
186  fhSubJettiness1_FJ_OP_AKT(0x0),
187  fhSubJettiness1_FJ_OP_KT(0x0),
188  fhSubJettiness1_FJ_OP_CA(0x0),
189  fhSubJettiness1_FJ_OP_WTA_KT(0x0),
190  fhSubJettiness1_FJ_OP_WTA_CA(0x0),
191  fhSubJettiness1_FJ_MIN(0x0),
192  fhSubJettiness2_FJ_AKT(0x0),
193  fhSubJettiness2_FJ_KT(0x0),
194  fhSubJettiness2_FJ_CA(0x0),
195  fhSubJettiness2_FJ_WTA_KT(0x0),
196  fhSubJettiness2_FJ_WTA_CA(0x0),
197  fhSubJettiness2_FJ_OP_AKT(0x0),
198  fhSubJettiness2_FJ_OP_KT(0x0),
199  fhSubJettiness2_FJ_OP_CA(0x0),
200  fhSubJettiness2_FJ_OP_WTA_KT(0x0),
201  fhSubJettiness2_FJ_OP_WTA_CA(0x0),
202  fhSubJettiness2_FJ_MIN(0x0),
203  fTreeResponseMatrixAxis(0)
204 
205 {
206  SetMakeGeneralHistograms(kTRUE);
207 }
208 
209 //________________________________________________________________________
210 AliAnalysisTaskSubJetFraction::AliAnalysisTaskSubJetFraction(const char *name) :
211  AliAnalysisTaskEmcalJet(name, kTRUE),
212  fContainer(0),
213  fMinFractionShared(0),
214  fJetShapeType(kData),
215  fJetShapeSub(kNoSub),
216  fJetSelection(kInclusive),
217  fShapesVar(0),
218  fPtThreshold(-9999.),
219  fRMatching(0.2),
220  fminpTTrig(20.),
221  fmaxpTTrig(50.),
222  fangWindowRecoil(0.6),
223  fSemigoodCorrect(0),
224  fHolePos(0),
225  fHoleWidth(0),
226  fCentSelectOn(kTRUE),
227  fCentMin(0),
228  fCentMax(10),
229  fSubJetAlgorithm(0),
230  fSubJetRadius(0.1),
231  fSubJetMinPt(1),
232  fJetRadius(0.4),
233  fRMatched(0.2),
234  fSharedFractionPtMin(0.5),
235  fDerivSubtrOrder(0),
236  fFullTree(kFALSE),
237  fhJetPt(0x0),
238  fhJetPt_1(0x0),
239  fhJetPt_2(0x0),
240  fhJetPhi(0x0),
241  fhJetPhi_1(0x0),
242  fhJetPhi_2(0x0),
243  fhJetEta(0x0),
244  fhJetEta_1(0x0),
245  fhJetEta_2(0x0),
246  fhJetMass(0x0),
247  fhJetMass_1(0x0),
248  fhJetMass_2(0x0),
249  fhJetRadius(0x0),
250  fhJetRadius_1(0x0),
251  fhJetRadius_2(0x0),
252  fhJetAngularity(0x0),
253  fhJetAngularityJetPt(0x0),
254  fhJetPTD(0x0),
255  fhJetPTDJetPt(0x0),
256  fhJetCounter(0x0),
257  fhJetCounter_1(0x0),
258  fhJetCounter_2(0x0),
259  fhNumberOfJetTracks(0x0),
260  fhNumberOfJetTracks_1(0x0),
261  fhNumberOfJetTracks_2(0x0),
262  fhSubJetPt(0x0),
263  fhSubJetPt_1(0x0),
264  fhSubJetPt_2(0x0),
265  fhSubJetMass(0x0),
266  fhSubJetMass_1(0x0),
267  fhSubJetMass_2(0x0),
268  fhSubJetRadius(0x0),
269  fhSubJetRadius_1(0x0),
270  fhSubJetRadius_2(0x0),
271  fhSubJetCounter(0x0),
272  fhSubJetCounter_1(0x0),
273  fhSubJetCounter_2(0x0),
274  fhNumberOfSubJetTracks(0x0),
275  fhNumberOfSubJetTracks_1(0x0),
276  fhNumberOfSubJetTracks_2(0x0),
277  fhSubJetPtFrac(0x0),
278  fhSubJetPtFrac2(0x0),
279  fhSubJetPtLoss(0x0),
280  fhSubJetPtLoss2(0x0),
281  fhSubJetEnergyFrac(0x0),
282  fhSubJetEnergyFrac2(0x0),
283  fhSubJetEnergyLoss(0x0),
284  fhSubJetEnergyLoss2(0x0),
285  fh2PtRatio(0x0),
286  fhSubJetiness1(0x0),
287  fhSubJetiness1JetPt(0x0),
288  fhSubJetiness2(0x0),
289  fhSubJetiness2JetPt(0x0),
290  fh2to1SubJetinessRatio(0x0),
291  fh2to1SubJetinessRatioJetPt(0x0),
292  fhEventCounter(0x0),
293  fhEventCounter_1(0x0),
294  fhEventCounter_2(0x0),
295  fhJetPtJetEta(0x0),
296  fhSubJetiness2Distance(0x0),
297  fh2JetTracksEtaPhiPt(0x0),
298  fhSubJettiness1CheckRatio_FJ_AKT(0x0),
299  fhSubJettiness1CheckRatio_FJ_KT(0x0),
300  fhSubJettiness1CheckRatio_FJ_CA(0x0),
301  fhSubJettiness1CheckRatio_FJ_WTA_KT(0x0),
302  fhSubJettiness1CheckRatio_FJ_WTA_CA(0x0),
303  fhSubJettiness1CheckRatio_FJ_OP_AKT(0x0),
304  fhSubJettiness1CheckRatio_FJ_OP_KT(0x0),
305  fhSubJettiness1CheckRatio_FJ_OP_CA(0x0),
306  fhSubJettiness1CheckRatio_FJ_OP_WTA_KT(0x0),
307  fhSubJettiness1CheckRatio_FJ_OP_WTA_CA(0x0),
308  fhSubJettiness1CheckRatio_FJ_MIN(0x0),
309  fhSubJettiness2CheckRatio_FJ_AKT(0x0),
310  fhSubJettiness2CheckRatio_FJ_KT(0x0),
311  fhSubJettiness2CheckRatio_FJ_CA(0x0),
312  fhSubJettiness2CheckRatio_FJ_WTA_KT(0x0),
313  fhSubJettiness2CheckRatio_FJ_WTA_CA(0x0),
314  fhSubJettiness2CheckRatio_FJ_OP_AKT(0x0),
315  fhSubJettiness2CheckRatio_FJ_OP_KT(0x0),
316  fhSubJettiness2CheckRatio_FJ_OP_CA(0x0),
317  fhSubJettiness2CheckRatio_FJ_OP_WTA_KT(0x0),
318  fhSubJettiness2CheckRatio_FJ_OP_WTA_CA(0x0),
319  fhSubJettiness2CheckRatio_FJ_MIN(0x0),
320  fhSubJettiness1_FJ_AKT(0x0),
321  fhSubJettiness1_FJ_KT(0x0),
322  fhSubJettiness1_FJ_CA(0x0),
323  fhSubJettiness1_FJ_WTA_KT(0x0),
324  fhSubJettiness1_FJ_WTA_CA(0x0),
325  fhSubJettiness1_FJ_OP_AKT(0x0),
326  fhSubJettiness1_FJ_OP_KT(0x0),
327  fhSubJettiness1_FJ_OP_CA(0x0),
328  fhSubJettiness1_FJ_OP_WTA_KT(0x0),
329  fhSubJettiness1_FJ_OP_WTA_CA(0x0),
330  fhSubJettiness1_FJ_MIN(0x0),
331  fhSubJettiness2_FJ_AKT(0x0),
332  fhSubJettiness2_FJ_KT(0x0),
333  fhSubJettiness2_FJ_CA(0x0),
334  fhSubJettiness2_FJ_WTA_KT(0x0),
335  fhSubJettiness2_FJ_WTA_CA(0x0),
336  fhSubJettiness2_FJ_OP_AKT(0x0),
337  fhSubJettiness2_FJ_OP_KT(0x0),
338  fhSubJettiness2_FJ_OP_CA(0x0),
339  fhSubJettiness2_FJ_OP_WTA_KT(0x0),
340  fhSubJettiness2_FJ_OP_WTA_CA(0x0),
341  fhSubJettiness2_FJ_MIN(0x0),
342  fTreeResponseMatrixAxis(0)
343 
344 {
345  // Standard constructor.
346 
347  SetMakeGeneralHistograms(kTRUE);
348 
349  DefineOutput(1, TTree::Class());
350 
351 }
352 
353 //________________________________________________________________________
354 AliAnalysisTaskSubJetFraction::~AliAnalysisTaskSubJetFraction()
355 {
356  // Destructor.
357 }
358 
359 //________________________________________________________________________
360  void AliAnalysisTaskSubJetFraction::UserCreateOutputObjects()
361 {
362  // Create user output.
363 
364  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
365 
366  Bool_t oldStatus = TH1::AddDirectoryStatus();
367  TH1::AddDirectory(kFALSE);
368 
369  //create a tree used for the MC data and making a 4D response matrix
370  fTreeResponseMatrixAxis = new TTree("fTreeJetShape", "fTreeJetShape");
371  if (fFullTree){
372  const Int_t nVar = 18;
373  fShapesVar = new Double_t [nVar]; //shapes used for tagging
374  TString *fShapesVarNames = new TString [nVar];
375 
376  fShapesVarNames[0] = "Pt";
377  fShapesVarNames[1] = "Pt_Truth";
378  fShapesVarNames[2] = "Tau1";
379  fShapesVarNames[3] = "Tau1_Truth";
380  fShapesVarNames[4] = "Tau2";
381  fShapesVarNames[5] = "Tau2_Truth";
382  fShapesVarNames[6] = "Tau3";
383  fShapesVarNames[7] = "Tau3_Truth";
384  fShapesVarNames[8] = "OpeningAngle";
385  fShapesVarNames[9] = "OpeningAngle_Truth";
386  fShapesVarNames[10] = "JetMultiplicity";
387  fShapesVarNames[11] = "JetMultiplicity_Truth";
388  fShapesVarNames[12] = "DeltaR";
389  fShapesVarNames[13] = "DeltaR_Truth";
390  fShapesVarNames[14] = "Frac1";
391  fShapesVarNames[15] = "Frac1_Truth";
392  fShapesVarNames[16] = "Frac2";
393  fShapesVarNames[17] = "Frac2_Truth";
394  for(Int_t ivar=0; ivar < nVar; ivar++){
395  cout<<"looping over variables"<<endl;
396  fTreeResponseMatrixAxis->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/D", fShapesVarNames[ivar].Data()));
397  }
398  }
399 
400  if (!fFullTree){
401  const Int_t nVar = 12;
402  fShapesVar = new Double_t [nVar]; //shapes used for tagging
403  TString *fShapesVarNames = new TString [nVar];
404 
405  fShapesVarNames[0] = "Pt";
406  fShapesVarNames[1] = "Pt_Truth";
407  fShapesVarNames[2] = "Tau1";
408  fShapesVarNames[3] = "Tau1_Truth";
409  fShapesVarNames[4] = "Tau2";
410  fShapesVarNames[5] = "Tau2_Truth";
411  fShapesVarNames[6] = "Tau3";
412  fShapesVarNames[7] = "Tau3_Truth";
413  fShapesVarNames[8] = "OpeningAngle";
414  fShapesVarNames[9] = "OpeningAngle_Truth";
415  fShapesVarNames[10] = "JetMultiplicity";
416  fShapesVarNames[11] = "JetMultiplicity_Truth";
417  for(Int_t ivar=0; ivar < nVar; ivar++){
418  cout<<"looping over variables"<<endl;
419  fTreeResponseMatrixAxis->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/D", fShapesVarNames[ivar].Data()));
420  }
421  }
422 
423 
424  if (fJetShapeType==AliAnalysisTaskSubJetFraction::kData || fJetShapeType==AliAnalysisTaskSubJetFraction::kSim || fJetShapeType==AliAnalysisTaskSubJetFraction::kGenOnTheFly){
425 
426  // fhJetPt= new TH1F("fhJetPt", "Jet Pt", (XBinsJetPtSize)-1, XBinsJetPt);
427  fhJetPt= new TH1F("fhJetPt", "Jet Pt",1500,-0.5,149.5 );
428  fOutput->Add(fhJetPt);
429  //fhJetPhi= new TH1F("fhJetPhi", "Jet Phi", Phi_Bins, Phi_Low, Phi_Up);
430  fhJetPhi= new TH1F("fhJetPhi", "Jet Phi",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
431  fOutput->Add(fhJetPhi);
432  fhJetEta= new TH1F("fhJetEta", "Jet Eta", Eta_Bins, Eta_Low, Eta_Up);
433  fOutput->Add(fhJetEta);
434  fhJetMass= new TH1F("fhJetMass", "Jet Mass", 4000,-0.5, 39.5);
435  fOutput->Add(fhJetMass);
436  fhJetRadius= new TH1F("fhJetRadius", "Jet Radius", 100, -0.05,0.995);
437  fOutput->Add(fhJetRadius);
438  fhNumberOfJetTracks= new TH1F("fhNumberOfJetTracks", "Number of Tracks within a Jet", 300, -0.5,299.5);
439  fOutput->Add(fhNumberOfJetTracks);
440  fhSubJetRadius= new TH1F("fhSubJetRadius", "SubJet Radius", 100, -0.05,0.995);
441  fOutput->Add(fhSubJetRadius);
442  fhSubJetPt= new TH1F("fhSubJetPt", "SubJet Pt", 1500, -0.5,149.5);
443  fOutput->Add(fhSubJetPt);
444  fhSubJetMass= new TH1F("fhSubJetMass", "Sub Jet Mass", 4000,-0.5, 39.5);
445  fOutput->Add(fhSubJetMass);
446  fhNumberOfSubJetTracks= new TH1F("fhNumberOfSubJetTracks", "Number of Tracks within a Sub Jet", 300, -0.5,299.5);
447  fOutput->Add(fhNumberOfSubJetTracks);
448  fhJetCounter= new TH1F("fhJetCounter", "Jet Counter", 150, -0.5, 149.5);
449  fOutput->Add(fhJetCounter);
450  fhSubJetCounter = new TH1F("fhSubJetCounter", "SubJet Counter",50, -0.5,49.5);
451  fOutput->Add(fhSubJetCounter);
452  fhEventCounter= new TH1F("fhEventCounter", "Event Counter", 15,0.5,15.5);
453  fOutput->Add(fhEventCounter);
454  }
455  if(fJetShapeType==AliAnalysisTaskSubJetFraction::kSim || fJetShapeType==AliAnalysisTaskSubJetFraction::kGenOnTheFly){
456  fhJetAngularity= new TH1F("fhJetAngularity", "Jet Angularity", 400, -2,2);
457  fOutput->Add(fhJetAngularity);
458  fhJetAngularityJetPt= new TH2F("fhJetAngularityJetPt", "Jet Angularity vs Jet Pt", 1500, -0.5, 149.50, 400, -2,2);
459  fOutput->Add(fhJetAngularityJetPt);
460  fhJetPTD= new TH1F("fhJetPTD", "Jet PTD", 400, -2,2);
461  fOutput->Add(fhJetPTD);
462  fhJetPTDJetPt= new TH2F("fhJetPTDJetPt", "Jet PTD vs Jet Pt", 1500, -0.5, 149.50, 400, -2,2);
463  fOutput->Add(fhJetPTDJetPt);
464  fhSubJetPtFrac= new TH1F("fhSubJetPtFrac", "Pt Fraction of Highest Pt Subjet compared to original Jet",101, -0.05,1.05);
465  fOutput->Add(fhSubJetPtFrac);
466  fhSubJetPtFrac2= new TH1F("fhSubJetPtFrac2", "Pt Fraction of Two Highest Pt Subjets compared to original Jet",101, -0.05,1.05);
467  fOutput->Add(fhSubJetPtFrac2);
468  fhSubJetPtLoss= new TH1F("fhSubJetPtLoss", "Pt Difference of Highest Pt Subjet compared to original Jet",101, -0.05,1.05);
469  fOutput->Add(fhSubJetPtLoss);
470  fhSubJetPtLoss2= new TH1F("fhSubJetPtLoss2", "Pt Difference of Two Highest Pt Subjets compared to original Jet",101, -0.05,1.05);
471  fOutput->Add(fhSubJetPtLoss2);
472  fhSubJetEnergyFrac= new TH1F("fhSubJetEnergyFrac", "Energy Fraction of Most Energetic Subjet compared to original Jet",101, -0.05,1.05);
473  fOutput->Add(fhSubJetEnergyFrac);
474  fhSubJetEnergyFrac2= new TH1F("fhSubJetEnergyFrac2", "Energy Fraction of Two Most Energetic Subjets compared to original Jet",101, -0.05,1.05);
475  fOutput->Add(fhSubJetEnergyFrac2);
476  fhSubJetEnergyLoss= new TH1F("fhSubJetEnergyLoss", "Energy Difference of Most Energetic Subjet compared to original Jet",101, -0.05,1.05);
477  fOutput->Add(fhSubJetEnergyLoss);
478  fhSubJetEnergyLoss2= new TH1F("fhSubJetEnergyLoss2", "Pt Difference of Two Most Energetic Subjets compared to original Jet",101, -0.05,1.05);
479  fOutput->Add(fhSubJetEnergyLoss2);
480  fhSubJetiness1= new TH1F("fhSubjetiness1", "Tau 1 value",101, -0.05,1.05);
481  fOutput->Add(fhSubJetiness1);
482  fhSubJetiness1JetPt= new TH2F("fhSubjetiness1JetPt", "Tau 1 value vs Jet Pt",1500, -0.5, 149.5, 101, -0.05,1.05);
483  fOutput->Add(fhSubJetiness1JetPt);
484  fhSubJetiness2= new TH1F("fhSubjetiness2", "Tau 2 value",101, -0.05,1.05);
485  fOutput->Add(fhSubJetiness2);
486  fhSubJetiness2JetPt= new TH2F("fhSubjetiness2JetPt", "Tau 2 value vs Jet Pt",1500, -0.5, 149.5, 101, -0.05,1.05);
487  fOutput->Add(fhSubJetiness2JetPt);
488  fh2to1SubJetinessRatio= new TH1F("fh2to1SubJetinessRatio", "Ratio of #tau 1 to #tau 2",200, -0.5,1.5);
489  fOutput->Add(fh2to1SubJetinessRatio);
490  fh2to1SubJetinessRatioJetPt= new TH2F("fh2to1SubJetinessRatioJetPt", "Ratio of #tau 1 to #tau 2 vs Jet Pt", 1500, -0.5, 149.5, 200, -0.5,1.5);
491  fOutput->Add(fh2to1SubJetinessRatioJetPt);
492  fhJetPtJetEta = new TH2F("fhJetPtJetEta", "Jet Pt vs Jet Eta", 1500,-0.5,150,Eta_Bins, Eta_Low, Eta_Up);
493  fOutput->Add(fhJetPtJetEta);
494  fhSubJetiness2Distance = new TH2F("fhSubJetiness2Distance", "#tau 2 as a function of distance between subject axis",100,0,10,101,-0.05,1.05);
495  fOutput->Add(fhSubJetiness2Distance);
496  fhSubJettiness1_FJ_KT= new TH1D("fhSubJettiness1_FJ_KT","fhSubJettiness1_FJ_KT",400,-2,2);
497  fOutput->Add(fhSubJettiness1_FJ_KT);
498  fhSubJettiness1_FJ_MIN= new TH1D("fhSubJettiness1_FJ_MIN","fhSubJettiness1_FJ_MIN",400,-2,2);
499  fOutput->Add(fhSubJettiness1_FJ_MIN);
500  fhSubJettiness2_FJ_KT= new TH1D("fhSubJettiness2_FJ_KT","fhSubJettiness2_FJ_KT",400,-2,2);
501  fOutput->Add(fhSubJettiness2_FJ_KT);
502  fhSubJettiness2_FJ_MIN= new TH1D("fhSubJettiness2_FJ_MIN","fhSubJettiness2_FJ_MIN",400,-2,2);
503  fOutput->Add(fhSubJettiness2_FJ_MIN);
504  fh2JetTracksEtaPhiPt = new TH2D("fh3JetTracksEtaPhiPt","fh3JetTracksEtaPgiPt",Eta_Bins, Eta_Low, Eta_Up,360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
505  fOutput->Add(fh2JetTracksEtaPhiPt);
506  fhSubJettiness1CheckRatio_FJ_AKT = new TH2D("fhSubJettiness1CheckRatio_FJ_AKT","fhSubJettiness1CheckRatio_FJ_AKT",400,-2,2,300,-1,2);
507  fOutput->Add(fhSubJettiness1CheckRatio_FJ_AKT);
508  fhSubJettiness1CheckRatio_FJ_KT= new TH2D("fhSubJettiness1CheckRatio_FJ_KT","fhSubJettiness1CheckRatio_FJ_KT",400,-2,2,300,-1,2);
509  fOutput->Add(fhSubJettiness1CheckRatio_FJ_KT);
510  fhSubJettiness1CheckRatio_FJ_CA= new TH2D("fhSubJettiness1CheckRatio_FJ_CA","fhSubJettiness1CheckRatio_FJ_CA",400,-2,2,300,-1,2);
511  fOutput->Add(fhSubJettiness1CheckRatio_FJ_CA);
512  fhSubJettiness1CheckRatio_FJ_WTA_KT= new TH2D("fhSubJettiness1CheckRatio_FJ_WTA_KT","fhSubJettiness1CheckRatio_FJ_WTA_KT",400,-2,2,300,-1,2);
513  fOutput->Add(fhSubJettiness1CheckRatio_FJ_WTA_KT);
514  fhSubJettiness1CheckRatio_FJ_WTA_CA= new TH2D("fhSubJettiness1CheckRatio_FJ_WTA_CA","fhSubJettiness1CheckRatio_FJ_WTA_CA",400,-2,2,300,-1,2);
515  fOutput->Add(fhSubJettiness1CheckRatio_FJ_WTA_CA);
516  fhSubJettiness1CheckRatio_FJ_OP_AKT= new TH2D("fhSubJettiness1CheckRatio_FJ_OP_AKT","fhSubJettiness1CheckRatio_FJ_OP_AKT",400,-2,2,300,-1,2);
517  fOutput->Add(fhSubJettiness1CheckRatio_FJ_OP_AKT);
518  fhSubJettiness1CheckRatio_FJ_OP_KT= new TH2D("fhSubJettiness1CheckRatio_FJ_OP_KT","fhSubJettiness1CheckRatio_FJ_OP_KT",400,-2,2,300,-1,2);
519  fOutput->Add(fhSubJettiness1CheckRatio_FJ_OP_KT);
520  fhSubJettiness1CheckRatio_FJ_OP_CA= new TH2D("fhSubJettiness1CheckRatio_FJ_OP_CA","fhSubJettiness1CheckRatio_FJ_OP_CA",400,-2,2,300,-1,2);
521  fOutput->Add(fhSubJettiness1CheckRatio_FJ_OP_CA);
522  fhSubJettiness1CheckRatio_FJ_OP_WTA_KT= new TH2D("fhSubJettiness1CheckRatio_FJ_OP_WTA_KT","fhSubJettiness1CheckRatio_FJ_OP_WTA_KT",400,-2,2,300,-1,2);
523  fOutput->Add(fhSubJettiness1CheckRatio_FJ_OP_WTA_KT);
524  fhSubJettiness1CheckRatio_FJ_OP_WTA_CA= new TH2D("fhSubJettiness1CheckRatio_FJ_OP_WTA_CA","fhSubJettiness1CheckRatio_FJ_OP_WTA_CA",400,-2,2,300,-1,2);
525  fOutput->Add(fhSubJettiness1CheckRatio_FJ_OP_WTA_CA);
526  fhSubJettiness1CheckRatio_FJ_MIN= new TH2D("fhSubJettiness1CheckRatio_FJ_MIN","fhSubJettiness1CheckRatio_FJ_MIN",400,-2,2,300,-1,2);
527  fOutput->Add(fhSubJettiness1CheckRatio_FJ_MIN);
528  fhSubJettiness2CheckRatio_FJ_AKT = new TH2D("fhSubJettiness2CheckRatio_FJ_AKT","fhSubJettiness2CheckRatio_FJ_AKT",400,-2,2,300,-1,2);
529  fOutput->Add(fhSubJettiness2CheckRatio_FJ_AKT);
530  fhSubJettiness2CheckRatio_FJ_KT= new TH2D("fhSubJettiness2CheckRatio_FJ_KT","fhSubJettiness2CheckRatio_FJ_KT",400,-2,2,300,-1,2);
531  fOutput->Add(fhSubJettiness2CheckRatio_FJ_KT);
532  fhSubJettiness2CheckRatio_FJ_CA= new TH2D("fhSubJettiness2CheckRatio_FJ_CA","fhSubJettiness2CheckRatio_FJ_CA",400,-2,2,300,-1,2);
533  fOutput->Add(fhSubJettiness2CheckRatio_FJ_CA);
534  fhSubJettiness2CheckRatio_FJ_WTA_KT= new TH2D("fhSubJettiness2CheckRatio_FJ_WTA_KT","fhSubJettiness2CheckRatio_FJ_WTA_KT",400,-2,2,300,-1,2);
535  fOutput->Add(fhSubJettiness2CheckRatio_FJ_WTA_KT);
536  fhSubJettiness2CheckRatio_FJ_WTA_CA= new TH2D("fhSubJettiness2CheckRatio_FJ_WTA_CA","fhSubJettiness2CheckRatio_FJ_WTA_CA",400,-2,2,300,-1,2);
537  fOutput->Add(fhSubJettiness2CheckRatio_FJ_WTA_CA);
538  fhSubJettiness2CheckRatio_FJ_OP_AKT= new TH2D("fhSubJettiness2CheckRatio_FJ_OP_AKT","fhSubJettiness2CheckRatio_FJ_OP_AKT",400,-2,2,300,-1,2);
539  fOutput->Add(fhSubJettiness2CheckRatio_FJ_OP_AKT);
540  fhSubJettiness2CheckRatio_FJ_OP_KT= new TH2D("fhSubJettiness2CheckRatio_FJ_OP_KT","fhSubJettiness2CheckRatio_FJ_OP_KT",400,-2,2,300,-1,2);
541  fOutput->Add(fhSubJettiness2CheckRatio_FJ_OP_KT);
542  fhSubJettiness2CheckRatio_FJ_OP_CA= new TH2D("fhSubJettiness2CheckRatio_FJ_OP_CA","fhSubJettiness2CheckRatio_FJ_OP_CA",400,-2,2,300,-1,2);
543  fOutput->Add(fhSubJettiness2CheckRatio_FJ_OP_CA);
544  fhSubJettiness2CheckRatio_FJ_OP_WTA_KT= new TH2D("fhSubJettiness2CheckRatio_FJ_OP_WTA_KT","fhSubJettiness2CheckRatio_FJ_OP_WTA_KT",400,-2,2,300,-1,2);
545  fOutput->Add(fhSubJettiness2CheckRatio_FJ_OP_WTA_KT);
546  fhSubJettiness2CheckRatio_FJ_OP_WTA_CA= new TH2D("fhSubJettiness2CheckRatio_FJ_OP_WTA_CA","fhSubJettiness2CheckRatio_FJ_OP_WTA_CA",400,-2,2,300,-1,2);
547  fOutput->Add(fhSubJettiness2CheckRatio_FJ_OP_WTA_CA);
548  fhSubJettiness2CheckRatio_FJ_MIN= new TH2D("fhSubJettiness2CheckRatio_FJ_MIN","fhSubJettiness2CheckRatio_FJ_MIN",400,-2,2,300,-1,2);
549  fOutput->Add(fhSubJettiness2CheckRatio_FJ_MIN);
550  fhSubJettiness1_FJ_AKT= new TH1D("fhSubJettiness1_FJ_AKT","fhSubJettiness1_FJ_AKT",400,-2,2);
551  fOutput->Add(fhSubJettiness1_FJ_AKT);
552  fhSubJettiness1_FJ_CA= new TH1D("fhSubJettiness1_FJ_CA","fhSubJettiness1_FJ_CA",400,-2,2);
553  fOutput->Add(fhSubJettiness1_FJ_CA);
554  fhSubJettiness1_FJ_WTA_KT= new TH1D("fhSubJettiness1_FJ_WTA_KT","fhSubJettiness1_FJ_WTA_KT",400,-2,2);
555  fOutput->Add(fhSubJettiness1_FJ_WTA_KT);
556  fhSubJettiness1_FJ_WTA_CA= new TH1D("fhSubJettiness1_FJ_WTA_CA","fhSubJettiness1_FJ_WTA_CA",400,-2,2);
557  fOutput->Add(fhSubJettiness1_FJ_WTA_CA);
558  fhSubJettiness1_FJ_OP_AKT= new TH1D("fhSubJettiness1_FJ_OP_AKT","fhSubJettiness1_FJ_OP_AKT",400,-2,2);
559  fOutput->Add(fhSubJettiness1_FJ_OP_AKT);
560  fhSubJettiness1_FJ_OP_KT= new TH1D("fhSubJettiness1_FJ_OP_KT","fhSubJettiness1_FJ_OP_KT",400,-2,2);
561  fOutput->Add(fhSubJettiness1_FJ_OP_KT);
562  fhSubJettiness1_FJ_OP_CA= new TH1D("fhSubJettiness1_FJ_OP_CA","fhSubJettiness1_FJ_OP_CA",400,-2,2);
563  fOutput->Add(fhSubJettiness1_FJ_OP_CA);
564  fhSubJettiness1_FJ_OP_WTA_KT= new TH1D("fhSubJettiness1_FJ_OP_WTA_KT","fhSubJettiness1_FJ_OP_WTA_KT",400,-2,2);
565  fOutput->Add(fhSubJettiness1_FJ_OP_WTA_KT);
566  fhSubJettiness1_FJ_OP_WTA_CA= new TH1D("fhSubJettiness1_FJ_OP_WTA_CA","fhSubJettiness1_FJ_OP_WTA_CA",400,-2,2);
567  fOutput->Add(fhSubJettiness1_FJ_OP_WTA_CA);
568  fhSubJettiness2_FJ_AKT= new TH1D("fhSubJettiness2_FJ_AKT","fhSubJettiness2_FJ_AKT",400,-2,2);
569  fOutput->Add(fhSubJettiness2_FJ_AKT);
570  fhSubJettiness2_FJ_CA= new TH1D("fhSubJettiness2_FJ_CA","fhSubJettiness2_FJ_CA",400,-2,2);
571  fOutput->Add(fhSubJettiness2_FJ_CA);
572  fhSubJettiness2_FJ_WTA_KT= new TH1D("fhSubJettiness2_FJ_WTA_KT","fhSubJettiness2_FJ_WTA_KT",400,-2,2);
573  fOutput->Add(fhSubJettiness2_FJ_WTA_KT);
574  fhSubJettiness2_FJ_WTA_CA= new TH1D("fhSubJettiness2_FJ_WTA_CA","fhSubJettiness2_FJ_WTA_CA",400,-2,2);
575  fOutput->Add(fhSubJettiness2_FJ_WTA_CA);
576  fhSubJettiness2_FJ_OP_AKT= new TH1D("fhSubJettiness2_FJ_OP_AKT","fhSubJettiness2_FJ_OP_AKT",400,-2,2);
577  fOutput->Add(fhSubJettiness2_FJ_OP_AKT);
578  fhSubJettiness2_FJ_OP_KT= new TH1D("fhSubJettiness2_FJ_OP_KT","fhSubJettiness2_FJ_OP_KT",400,-2,2);
579  fOutput->Add(fhSubJettiness2_FJ_OP_KT);
580  fhSubJettiness2_FJ_OP_CA= new TH1D("fhSubJettiness2_FJ_OP_CA","fhSubJettiness2_FJ_OP_CA",400,-2,2);
581  fOutput->Add(fhSubJettiness2_FJ_OP_CA);
582  fhSubJettiness2_FJ_OP_WTA_KT= new TH1D("fhSubJettiness2_FJ_OP_WTA_KT","fhSubJettiness2_FJ_OP_WTA_KT",400,-2,2);
583  fOutput->Add(fhSubJettiness2_FJ_OP_WTA_KT);
584  fhSubJettiness2_FJ_OP_WTA_CA= new TH1D("fhSubJettiness2_FJ_OP_WTA_CA","fhSubJettiness2_FJ_OP_WTA_CA",400,-2,2);
585  fOutput->Add(fhSubJettiness2_FJ_OP_WTA_CA);
586  }
587  if(fJetShapeType==AliAnalysisTaskSubJetFraction::kTrueDet){
588  fhJetPt_1= new TH1F("fhJetPt_1", "Jet Pt Detector Level",1500,-0.5,149.5 );
589  fOutput->Add(fhJetPt_1);
590  fhJetPt_2= new TH1F("fhJetPt_2", "Jet Pt Particle Level",1500,-0.5,149.5 );
591  fOutput->Add(fhJetPt_2);
592  fhJetPhi_1= new TH1F("fhJetPhi_1", "Jet Phi Detector Level",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
593  fOutput->Add(fhJetPhi_1);
594  fhJetPhi_2= new TH1F("fhJetPhi_2", "Jet Phi Particle Level",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
595  fOutput->Add(fhJetPhi_2);
596  fhJetEta_1= new TH1F("fhJetEta_1", "Jet Eta Detector Level", Eta_Bins, Eta_Low, Eta_Up);
597  fOutput->Add(fhJetEta_1);
598  fhJetEta_2= new TH1F("fhJetEta_2", "Jet Eta Particle Level", Eta_Bins, Eta_Low, Eta_Up);
599  fOutput->Add(fhJetEta_2);
600  fhJetMass_1= new TH1F("fhJetMass_1", "Jet Mass Detector Level", 4000,-0.5, 39.5);
601  fOutput->Add(fhJetMass_1);
602  fhJetMass_2= new TH1F("fhJetMass_2", "Jet Mass Particle Level", 4000,-0.5, 39.5);
603  fOutput->Add(fhJetMass_2);
604  fhJetRadius_1= new TH1F("fhJetRadius_1", "Jet Radius Detector Level", 100, -0.05,0.995);
605  fOutput->Add(fhJetRadius_1);
606  fhJetRadius_2= new TH1F("fhJetRadius_2", "Jet Radius Particle Level", 100, -0.05,0.995);
607  fOutput->Add(fhJetRadius_2);
608  fhNumberOfJetTracks_1= new TH1F("fhNumberOfJetTracks_1", "Number of Tracks within a Jet Detector Level", 300, -0.5,299.5);
609  fOutput->Add(fhNumberOfJetTracks_1);
610  fhNumberOfJetTracks_2= new TH1F("fhNumberOfJetTracks_2", "Number of Tracks within a Jet Particle Level", 300, -0.5,299.5);
611  fOutput->Add(fhNumberOfJetTracks_2);
612  fhSubJetRadius_1= new TH1F("fhSubJetRadius_1", "SubJet Radius Detector Level", 100, -0.05,0.995);
613  fOutput->Add(fhSubJetRadius_1);
614  fhSubJetRadius_2= new TH1F("fhSubJetRadius_2", "SubJet Radius Particle Level", 100, -0.05,0.995);
615  fOutput->Add(fhSubJetRadius_2);
616  fhSubJetPt_1= new TH1F("fhSubJetPt_1", "SubJet Pt Detector Level", 1500, -0.5,149.5);
617  fOutput->Add(fhSubJetPt_1);
618  fhSubJetPt_2= new TH1F("fhSubJetPt_2", "SubJet Pt Particle Level", 1500, -0.5,149.5);
619  fOutput->Add(fhSubJetPt_2);
620  fhSubJetMass_1= new TH1F("fhSubJetMass_1", "Sub Jet Mass Detector Level", 4000,-0.5, 39.5);
621  fOutput->Add(fhSubJetMass_1);
622  fhSubJetMass_2= new TH1F("fhSubJetMass_2", "Sub Jet Mass Particle Level", 4000,-0.5, 39.5);
623  fOutput->Add(fhSubJetMass_2);
624  fhNumberOfSubJetTracks_1= new TH1F("fhNumberOfSubJetTracks_1", "Number of Tracks within a Sub Jet Detector Level", 300, -0.5,299.5);
625  fOutput->Add(fhNumberOfSubJetTracks_1);
626  fhNumberOfSubJetTracks_2= new TH1F("fhNumberOfSubJetTracks_2", "Number of Tracks within a Sub Jet Particle Level", 300, -0.5,299.5);
627  fOutput->Add(fhNumberOfSubJetTracks_2);
628  fhJetCounter_1= new TH1F("fhJetCounter_1", "Jet Counter Detector Level", 150, -0.5, 149.5);
629  fOutput->Add(fhJetCounter_1);
630  fhJetCounter_2= new TH1F("fhJetCounter_2", "Jet Counter Particle Level", 150, -0.5, 149.5);
631  fOutput->Add(fhJetCounter_2);
632  fhSubJetCounter_1 = new TH1F("fhSubJetCounter_1", "SubJet Counter Detector Level",50, -0.5,49.5);
633  fOutput->Add(fhSubJetCounter_1);
634  fhSubJetCounter_2 = new TH1F("fhSubJetCounter_2", "SubJet Counter Particle Level",50, -0.5,49.5);
635  fOutput->Add(fhSubJetCounter_2);
636  fh2PtRatio= new TH2F("fhPtRatio", "MC pt for detector level vs particle level jets",1500,-0.5,149.5,1500,-0.5,149.5);
637  fOutput->Add(fh2PtRatio);
638  fhEventCounter_1= new TH1F("fhEventCounter_1", "Event Counter Detector Level", 15,0.5,15.5);
639  fOutput->Add(fhEventCounter_1);
640  fhEventCounter_2= new TH1F("fhEventCounter_2", "Event Counter Particle Level", 15,0.5,15.5);
641  fOutput->Add(fhEventCounter_2);
642  }
643  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kDetEmbPart){
644  fhEventCounter= new TH1F("fhEventCounter", "Event Counter", 15,0.5,15.5);
645  fOutput->Add(fhEventCounter);
646  }
647  fOutput->Add(fTreeResponseMatrixAxis);
648  TH1::AddDirectory(oldStatus);
649  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
650 }
651 
652 //________________________________________________________________________
653 Bool_t AliAnalysisTaskSubJetFraction::Run()
654 {
655  // Run analysis code here, if needed. It will be executed before FillHistograms().
656 
657 
658  return kTRUE;
659 }
660 
661 //________________________________________________________________________
662 Bool_t AliAnalysisTaskSubJetFraction::FillHistograms()
663 {
664 
665  //fhEventCounter Key:
666  // 1: Number of events with a Jet Container
667  // 2: Number of Jets without a Jet Container
668  // 3:
669  // 4: Number of Jets found in all events
670  // 5: Number of Jets that were reclustered in all events
671  // 6: Number of SubJets found in all events
672  // 7: Number of events were primary vertext was found
673  // 8: Number of Jets with more than one SubJet
674  // 9:Number of Jets with more than two SubJets
675  // 12:Number of SubJetinessEvents in kData
676 
677 
678  // const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
679  // Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
680  // if(vert) fhEventCounter->Fill(7);
681  if (fCentSelectOn){
682  if ((fCent>fCentMax) || (fCent<fCentMin)) return 0;
683  }
684 
686  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kDetEmbPart){
687  AliEmcalJet *Jet1 = NULL; //Subtracted hybrid Jet
688  AliEmcalJet *Jet2 = NULL; //Unsubtracted Hybrid Jet
689  AliEmcalJet *Jet3 = NULL; //Detector Level Pythia Jet
690  AliEmcalJet *Jet4 = NULL; //Particle Level Pyhtia Jet
691  AliJetContainer *JetCont1= GetJetContainer(0); //Jet Container for Subtracted Hybrid Jets
692  AliJetContainer *JetCont2= GetJetContainer(1); //Jet Container for Unsubtracted Hybrid Jets
693  AliJetContainer *JetCont3= GetJetContainer(2); //Jet Container for Detector Level Pyhtia Jets
694  AliJetContainer *JetCont4= GetJetContainer(3); //Jet Container for Particle Level Pythia Jets
695  AliEmcalJetFinder *Reclusterer1; //Object containg Subjets from Subtracted Hybrid Jets
696  AliEmcalJetFinder *Reclusterer4; //Object containg Subjets from Particle Level Pythia Jets
697  Bool_t JetsMatched=kFALSE;
698  Bool_t EventCounter=kFALSE;
699  Int_t JetNumber=-1;
700  Double_t JetPtThreshold=-2;
701  fhEventCounter->Fill(1);
702  if(JetCont1) {
703  fhEventCounter->Fill(2);
704  JetCont1->ResetCurrentID();
705  while((Jet1=JetCont1->GetNextAcceptJet())) {
706  if (fJetShapeSub==kConstSub) JetPtThreshold=Jet1->Pt();
707  else JetPtThreshold=Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
708  if ( (!Jet1) || (JetPtThreshold<fPtThreshold)) continue;
709  else {
710  if(fSemigoodCorrect){
711  Double_t disthole=RelativePhi(Jet1->Phi(),fHolePos);
712  if(TMath::Abs(disthole)<fHoleWidth){
713  continue;}
714  }
715  if (!EventCounter){
716  fhEventCounter->Fill(3);
717  EventCounter=kTRUE;
718  }
719  if(fJetShapeSub==kConstSub){
720  JetNumber=-1;
721  for(Int_t i = 0; i<JetCont2->GetNJets(); i++) {
722  Jet2 = JetCont2->GetJet(i);
723  if(Jet2->GetLabel()==Jet1->GetLabel()) {
724  JetNumber=i;
725  }
726  }
727  if(JetNumber==-1) continue;
728  Jet2=JetCont2->GetJet(JetNumber);
729  if (JetCont2->AliJetContainer::GetFractionSharedPt(Jet2)<fSharedFractionPtMin) continue;
730  Jet3=Jet2->ClosestJet();
731  }
732  if(!(fJetShapeSub==kConstSub)){
733  if (!(JetCont1->AliJetContainer::GetFractionSharedPt(Jet1)<fSharedFractionPtMin)) continue;
734  Jet3 = Jet1->ClosestJet();
735  }
736  if (!Jet3) continue;
737  Jet4=Jet3->ClosestJet();
738  if(!Jet4) continue;
739  JetsMatched=kTRUE;
741  if (fJetShapeSub==kConstSub) fShapesVar[0]=Jet1->Pt();
742  else fShapesVar[0]=Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
743  fShapesVar[2]=fjNSubJettiness(Jet1,0,1,0,1,0);
744  fShapesVar[4]=fjNSubJettiness(Jet1,0,2,0,1,0);
745  fShapesVar[6]=fjNSubJettiness(Jet1,0,3,0,1,0);
746  fShapesVar[8]=fjNSubJettiness(Jet1,0,2,0,1,1);
747  fShapesVar[10]=Jet1->GetNumberOfTracks();
748  if (fFullTree){
749  fShapesVar[12]=fjNSubJettiness(Jet1,0,2,0,1,2);
750  Reclusterer1 = Recluster(Jet1, 0, fSubJetRadius, fSubJetMinPt, fSubJetAlgorithm, "SubJetFinder_1");
751  fShapesVar[14]=SubJetFraction(Jet1, Reclusterer1, 1, 0, kTRUE, kFALSE);
752  fShapesVar[16]=SubJetFraction(Jet1, Reclusterer1, 2, 0, kTRUE, kFALSE);
753  }
754  if (JetsMatched){ //even needed? Not now but might be if you want to fill trees when jets aren't matched too
755  fShapesVar[1]=Jet4->Pt();
756  fShapesVar[3]=fjNSubJettiness(Jet4,3,1,0,1,0);
757  fShapesVar[5]=fjNSubJettiness(Jet4,3,2,0,1,0);
758  fShapesVar[7]=fjNSubJettiness(Jet4,3,3,0,1,0);
759  fShapesVar[9]=fjNSubJettiness(Jet4,3,2,0,1,1);
760  fShapesVar[11]=Jet4->GetNumberOfTracks();
761  if (fFullTree){
762  fShapesVar[13]=fjNSubJettiness(Jet4,3,2,0,1,2);
763  Reclusterer4=Recluster(Jet4, 3, fSubJetRadius, 0, fSubJetAlgorithm, "SubJetFinder_4");
764  fShapesVar[15]=SubJetFraction(Jet4, Reclusterer4, 1, 0, kTRUE, kFALSE);
765  fShapesVar[17]=SubJetFraction(Jet4, Reclusterer4, 2, 0, kTRUE, kFALSE);
766  }
767  }
768  else{
769  fShapesVar[1]=-2;
770  fShapesVar[3]=-2;
771  fShapesVar[5]=-2;
772  fShapesVar[7]=-2;
773  fShapesVar[9]=-2;
774  fShapesVar[11]=-2;
775  if (fFullTree){
776  fShapesVar[13]=-2;
777  fShapesVar[15]=-2;
778  fShapesVar[17]=-2;
779  }
780  }
781  fTreeResponseMatrixAxis->Fill();
782  JetsMatched=kFALSE;
783  }
784  }
785  }
786  }
787 
789  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kTrueDet){
790  AliEmcalJet *Jet1 = NULL; //Detector Level Jet
791  AliEmcalJet *Jet2 = NULL; //Particle Level Jet
792  AliJetContainer *JetCont1= GetJetContainer(0); //Jet Container for Detector Level Pythia
793  AliJetContainer *JetCont2= GetJetContainer(1); //Jet Container for Particle Level Pythia
794  AliEmcalJetFinder *Reclusterer1; //Object containg Subjets Detector Level
795  AliEmcalJetFinder *Reclusterer2; //Object containg Subjets Particle Level
796  Int_t JetCounter1=0; //Counts number of jets in event
797  Int_t JetCounter2=0; //Counts number of jets in event
798  Double_t JetPhi1=0;
799  Double_t JetPhi2=0;
800  Bool_t JetsMatched=kFALSE;
801  Double_t Pythia_Event_Weight=1;
802  Bool_t EventCounter=kFALSE;
803  fhEventCounter_1->Fill(1);
804  if(JetCont1) {
805  fhEventCounter_1->Fill(2); //Number of events with a jet container
806  JetCont1->ResetCurrentID();
807  while((Jet1=JetCont1->GetNextAcceptJet())) {
808  if( (!Jet1) || ((Jet1->Pt())<fPtThreshold)) {
809  // fhEventCounter_1->Fill(3); //events where the jet had a problem
810  continue;
811  }
812  else {
813  if(fSemigoodCorrect){
814  Double_t disthole=RelativePhi(Jet1->Phi(),fHolePos);
815  if(TMath::Abs(disthole)<fHoleWidth){
816  continue;}
817  }
818  if (!EventCounter){
819  fhEventCounter_1->Fill(3);
820  EventCounter=kTRUE;
821  }
822  if((Jet1->GetNumberOfTracks())==0){
823  fhEventCounter_1->Fill(10); //zero track jets
824  }
825  if((Jet1->GetNumberOfTracks())==1){
826  fhEventCounter_1->Fill(11); //one track jets
827  }
828  fhEventCounter_1->Fill(4); //Number of Jets found in all events
829  JetCounter1++;
830  fhJetPt_1->Fill(Jet1->Pt());
831  JetPhi1=Jet1->Phi();
832  if(JetPhi1 < -1*TMath::Pi()) JetPhi1 += (2*TMath::Pi());
833  else if (JetPhi1 > TMath::Pi()) JetPhi1 -= (2*TMath::Pi());
834  fhJetPhi_1->Fill(JetPhi1);
835  fhJetEta_1->Fill(Jet1->Eta());
836  fhJetMass_1->Fill(Jet1->M());
837  fhJetRadius_1->Fill(TMath::Sqrt((Jet1->Area()/TMath::Pi()))); //Radius of Jets per event
838  fhNumberOfJetTracks_1->Fill(Jet1->GetNumberOfTracks());
839  if((Jet2 = Jet1->ClosestJet())){
840  JetsMatched=kTRUE;
841  if((Jet2->GetNumberOfTracks())==0){
842  fhEventCounter_2->Fill(10); //zero track jets
843  }
844  if((Jet2->GetNumberOfTracks())==1){
845  fhEventCounter_2->Fill(11); //one track jets
846  }
847  fhEventCounter_2->Fill(4); //Number of Jets found in all events
848  JetCounter2++;
849  fhJetPt_2->Fill(Jet2->Pt());
850  JetPhi2=Jet2->Phi();
851  if(JetPhi2 < -1*TMath::Pi()) JetPhi2 += (2*TMath::Pi());
852  else if (JetPhi2 > TMath::Pi()) JetPhi2 -= (2*TMath::Pi());
853  fhJetPhi_2->Fill(JetPhi2);
854  fhJetEta_2->Fill(Jet2->Eta());
855  fhJetMass_2->Fill(Jet2->M());
856  fhJetRadius_2->Fill(TMath::Sqrt((Jet2->Area()/TMath::Pi()))); //Radius of Jets per event
857  fhNumberOfJetTracks_2->Fill(Jet2->GetNumberOfTracks());
858  fh2PtRatio->Fill(Jet1->Pt(),Jet2->Pt(),Pythia_Event_Weight);
859  }
860  else {
861  fhEventCounter_2->Fill(1);
862  //continue;
863  }
865 
866 
867  fShapesVar[0]=Jet1->Pt();
868  fShapesVar[2]=fjNSubJettiness(Jet1,0,1,0,1,0);
869  fShapesVar[4]=fjNSubJettiness(Jet1,0,2,0,1,0);
870  fShapesVar[6]=fjNSubJettiness(Jet1,0,3,0,1,0);
871  fShapesVar[8]=fjNSubJettiness(Jet1,0,2,0,1,1);
872  fShapesVar[10]=Jet1->GetNumberOfTracks();
873  Reclusterer1 = Recluster(Jet1, 0, fSubJetRadius, fSubJetMinPt, fSubJetAlgorithm, "SubJetFinder_1");
874  if (fFullTree){
875  fShapesVar[12]=fjNSubJettiness(Jet1,0,2,0,1,2);
876  fShapesVar[14]=SubJetFraction(Jet1, Reclusterer1, 1, 0, kTRUE, kFALSE);
877  fShapesVar[16]=SubJetFraction(Jet1, Reclusterer1, 2, 0, kTRUE, kFALSE);
878  }
879  if (JetsMatched){ //even needed? Not now but might be if you want to fill trees when jets aren't matched too
880  fShapesVar[1]=Jet2->Pt();
881  fShapesVar[3]=fjNSubJettiness(Jet2,1,1,0,1,0);
882  fShapesVar[5]=fjNSubJettiness(Jet2,1,2,0,1,0);
883  fShapesVar[7]=fjNSubJettiness(Jet2,1,3,0,1,0);
884  fShapesVar[9]=fjNSubJettiness(Jet2,1,2,0,1,1);
885  fShapesVar[11]=Jet2->GetNumberOfTracks();
886  Reclusterer2 = Recluster(Jet2, 1, fSubJetRadius, 0, fSubJetAlgorithm, "SubJetFinder_2");
887  if (fFullTree){
888  fShapesVar[13]=fjNSubJettiness(Jet2,1,2,0,1,2);
889  fShapesVar[15]=SubJetFraction(Jet2, Reclusterer2, 1, 0, kTRUE, kFALSE);
890  fShapesVar[17]=SubJetFraction(Jet2, Reclusterer2, 2, 0, kTRUE, kFALSE);
891  }
892  }
893  else{
894  fShapesVar[1]=-2;
895  fShapesVar[3]=-2;
896  fShapesVar[5]=-2;
897  fShapesVar[7]=-2;
898  fShapesVar[9]=-2;
899  fShapesVar[11]=-2;
900  if (fFullTree){
901  fShapesVar[13]=-2;
902  fShapesVar[15]=-2;
903  fShapesVar[17]=-2;
904  }
905  }
906  fTreeResponseMatrixAxis->Fill();
907 
908  fhSubJetCounter_1->Fill(Reclusterer1->GetNumberOfJets());
909  for (Int_t i= 0; i<Reclusterer1->GetNumberOfJets(); i++){
910  fhEventCounter_1->Fill(6); //Number of overall subjets in all events
911  fhSubJetPt_1->Fill(Reclusterer1->GetJet(i)->Pt());
912  fhSubJetMass_1->Fill(Reclusterer1->GetJet(i)->M());
913  fhNumberOfSubJetTracks_1->Fill(Reclusterer1->GetJet(i)->GetNumberOfTracks());
914  fhSubJetRadius_1->Fill(TMath::Sqrt((Reclusterer1->GetJet(i)->Area()/TMath::Pi()))); //Radius of SubJets per event
915  }
916  if(JetsMatched){
917  fhSubJetCounter_2->Fill(Reclusterer2->GetNumberOfJets());
918  for (Int_t i= 0; i<Reclusterer2->GetNumberOfJets(); i++){
919  fhEventCounter_2->Fill(6); //Number of overall subjets in all events
920  fhSubJetPt_2->Fill(Reclusterer2->GetJet(i)->Pt());
921  fhSubJetMass_2->Fill(Reclusterer2->GetJet(i)->M());
922  fhNumberOfSubJetTracks_2->Fill(Reclusterer2->GetJet(i)->GetNumberOfTracks());
923  fhSubJetRadius_2->Fill(TMath::Sqrt((Reclusterer2->GetJet(i)->Area()/TMath::Pi()))); //Radius of SubJets per event
924  }
925  }
926  JetsMatched=kFALSE;
927  }
928  }
929  fhJetCounter_1->Fill(JetCounter1); //Number of Jets in Each Event Particle Level
930  fhJetCounter_2->Fill(JetCounter2); //Number of Jets in Each Event Detector Level
931  }
932  //else {fhEventCounter_1->Fill(3);} //Events with no jet container
933  }
934 
935 
936 
937  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kData){
938  AliEmcalJet *Jet1 = NULL; //Original Jet in the event
939  AliJetContainer *JetCont= GetJetContainer(0); //Jet Container for event
940  Int_t JetCounter=0; //Counts number of jets in event
941  Double_t JetPhi=0;
942  Bool_t EventCounter=kFALSE;
943  Double_t JetPt_ForThreshold=0;
944  fhEventCounter->Fill(1);
945  if(JetCont) {
946  fhEventCounter->Fill(2); //Number of events with a jet container
947  JetCont->ResetCurrentID();
948  while((Jet1=JetCont->GetNextAcceptJet())) {
949  if(!Jet1) continue;
950  if (fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) JetPt_ForThreshold = Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
951  else JetPt_ForThreshold = Jet1->Pt();
952  if(JetPt_ForThreshold<fPtThreshold) {
953  //fhEventCounter->Fill(3); //events where the jet had a problem
954  continue;
955  }
956  else {
957  if(fSemigoodCorrect){
958  Double_t disthole=RelativePhi(Jet1->Phi(),fHolePos);
959  if(TMath::Abs(disthole)<fHoleWidth){
960  continue;}
961  }
962  if (!EventCounter){
963  fhEventCounter->Fill(3);
964  EventCounter=kTRUE;
965  }
966  if((Jet1->GetNumberOfTracks())==0){
967  fhEventCounter->Fill(10); //zero track jets
968  }
969  if((Jet1->GetNumberOfTracks())==1){
970  fhEventCounter->Fill(11); //one track jets
971  }
972  fhEventCounter->Fill(4); //Number of Jets found in all events
973  JetCounter++;
974  fhJetPt->Fill(Jet1->Pt());
975  JetPhi=Jet1->Phi();
976  if(JetPhi < -1*TMath::Pi()) JetPhi += (2*TMath::Pi());
977  else if (JetPhi > TMath::Pi()) JetPhi -= (2*TMath::Pi());
978  fhJetPhi->Fill(JetPhi);
979  fhJetEta->Fill(Jet1->Eta());
980  fhJetMass->Fill(Jet1->M());
981  fhJetRadius->Fill(TMath::Sqrt((Jet1->Area()/TMath::Pi()))); //Radius of Jets per event
982  fhNumberOfJetTracks->Fill(Jet1->GetNumberOfTracks());
983  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) fShapesVar[0]= Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
984  else fShapesVar[0]=Jet1->Pt();
985  fShapesVar[2]=fjNSubJettiness(Jet1,0,1,0,1,0);
986  fShapesVar[4]=fjNSubJettiness(Jet1,0,2,0,1,0);
987  fShapesVar[6]=fjNSubJettiness(Jet1,0,3,0,1,0);
988  fShapesVar[8]=fjNSubJettiness(Jet1,0,2,0,1,1);
989  fShapesVar[10]=Jet1->GetNumberOfTracks();
990  AliEmcalJetFinder *Reclusterer1 = Recluster(Jet1, 0, fSubJetRadius, fSubJetMinPt, fSubJetAlgorithm, "SubJetFinder");
991  if (fFullTree){
992  fShapesVar[12]=fjNSubJettiness(Jet1,0,2,0,1,2);
993  fShapesVar[14]=SubJetFraction(Jet1, Reclusterer1, 1, 0, kTRUE, kFALSE);
994  fShapesVar[16]=SubJetFraction(Jet1, Reclusterer1, 2, 0, kTRUE, kFALSE);
995  }
996  fShapesVar[1]=-2;
997  fShapesVar[3]=-2;
998  fShapesVar[5]=-2;
999  fShapesVar[7]=-2;
1000  fShapesVar[9]=-2;
1001  fShapesVar[11]=-2;
1002  if (fFullTree){
1003  fShapesVar[13]=-2;
1004  fShapesVar[15]=-2;
1005  fShapesVar[17]=-2;
1006  }
1007  fTreeResponseMatrixAxis->Fill();
1008  fhSubJetCounter->Fill(Reclusterer1->GetNumberOfJets());
1009  for (Int_t i= 0; i<Reclusterer1->GetNumberOfJets(); i++){
1010  fhEventCounter->Fill(6); //Number of overall subjets in all events
1011  fhSubJetPt->Fill(Reclusterer1->GetJet(i)->Pt());
1012  fhSubJetMass->Fill(Reclusterer1->GetJet(i)->M());
1013  fhNumberOfSubJetTracks->Fill(Reclusterer1->GetJet(i)->GetNumberOfTracks());
1014  fhSubJetRadius->Fill(TMath::Sqrt((Reclusterer1->GetJet(i)->Area()/TMath::Pi()))); //Radius of SubJets per event
1015  }
1016  }
1017  }
1018  fhJetCounter->Fill(JetCounter); //Number of Jets in Each Event
1019  }
1020  //else {fhEventCounter->Fill(2);} //Events with no jet container
1021  }
1022 
1023 
1024  if (fJetShapeType==AliAnalysisTaskSubJetFraction::kSim || fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly){
1025 
1026  AliEmcalJet *Jet1 = NULL; //Original Jet in the event
1027  AliJetContainer *JetCont= GetJetContainer(0); //Jet Container for event
1028  Int_t JetCounter=0; //Counts number of jets in event
1029  Double_t JetPhi=0;
1030  Bool_t EventCounter=kFALSE;
1031  Double_t JetPt_ForThreshold=0;
1032  fhEventCounter->Fill(1);
1033  if(JetCont) {
1034  fhEventCounter->Fill(2); //Number of events with a jet container
1035  JetCont->ResetCurrentID();
1036  while((Jet1=JetCont->GetNextAcceptJet())) {
1037  // if (((Jet1=JetCont->GetLeadingJet("rho")) && (fPtThreshold<=Jet1->Pt()-(GetRhoVal(0)*Jet1->Area())))){
1038  if(!Jet1) continue;
1039  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) JetPt_ForThreshold = Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
1040  else JetPt_ForThreshold = Jet1->Pt();
1041  if(JetPt_ForThreshold<fPtThreshold) {
1042  //fhEventCounter->Fill(3); //events where the jet had a problem
1043  continue;
1044  }
1045  else {
1046  if(fSemigoodCorrect){
1047  Double_t disthole=RelativePhi(Jet1->Phi(),fHolePos);
1048  if(TMath::Abs(disthole)<fHoleWidth){
1049  continue;}
1050  }
1051  if (!EventCounter){
1052  fhEventCounter->Fill(3);
1053  EventCounter=kTRUE;
1054  }
1055  if((Jet1->GetNumberOfTracks())==0){
1056  fhEventCounter->Fill(10); //zero track jets
1057  }
1058  if((Jet1->GetNumberOfTracks())==1){
1059  fhEventCounter->Fill(11); //one track jets
1060  }
1061  fhEventCounter->Fill(4); //Number of Jets found in all events
1062  JetCounter++;
1063  fhJetPt->Fill(Jet1->Pt());
1064  JetPhi=Jet1->Phi();
1065  if(JetPhi < -1*TMath::Pi()) JetPhi += (2*TMath::Pi());
1066  else if (JetPhi > TMath::Pi()) JetPhi -= (2*TMath::Pi());
1067  fhJetPhi->Fill(JetPhi);
1068  fhJetEta->Fill(Jet1->Eta());
1069  fhJetMass->Fill(Jet1->M());
1070  fhJetRadius->Fill(TMath::Sqrt((Jet1->Area()/TMath::Pi()))); //Radius of Jets per event
1071  fhNumberOfJetTracks->Fill(Jet1->GetNumberOfTracks());
1072 
1073  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
1074  Double_t FJ_Beta=1.0;
1075  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));
1076  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));
1077  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));
1078  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));
1079  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));
1080  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));
1081  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));
1082  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));
1083  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));
1084  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));
1085  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));
1086  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));
1087  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));
1088  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));
1089  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));
1090  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));
1091  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));
1092  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));
1093  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));
1094  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));
1095  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));
1096  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));
1097 
1098  fhSubJettiness1_FJ_KT->Fill(fjNSubJettiness(Jet1,0,1,0,FJ_Beta,0));
1099  fhSubJettiness1_FJ_CA->Fill(fjNSubJettiness(Jet1,0,1,1,FJ_Beta,0));
1100  fhSubJettiness1_FJ_AKT->Fill(fjNSubJettiness(Jet1,0,1,2,FJ_Beta,0));
1101  fhSubJettiness1_FJ_WTA_KT->Fill(fjNSubJettiness(Jet1,0,1,3,FJ_Beta,0));
1102  fhSubJettiness1_FJ_WTA_CA->Fill(fjNSubJettiness(Jet1,0,1,4,FJ_Beta,0));
1103  fhSubJettiness1_FJ_OP_KT->Fill(fjNSubJettiness(Jet1,0,1,5,FJ_Beta,0));
1104  fhSubJettiness1_FJ_OP_CA->Fill(fjNSubJettiness(Jet1,0,1,6,FJ_Beta,0));
1105  fhSubJettiness1_FJ_OP_AKT->Fill(fjNSubJettiness(Jet1,0,1,7,FJ_Beta,0));
1106  fhSubJettiness1_FJ_OP_WTA_KT->Fill(fjNSubJettiness(Jet1,0,1,8,FJ_Beta,0));
1107  fhSubJettiness1_FJ_OP_WTA_CA->Fill(fjNSubJettiness(Jet1,0,1,9,FJ_Beta,0));
1108  fhSubJettiness1_FJ_MIN->Fill(fjNSubJettiness(Jet1,0,1,10,FJ_Beta,0));
1109  fhSubJettiness2_FJ_KT->Fill(fjNSubJettiness(Jet1,0,2,0,FJ_Beta,0));
1110  fhSubJettiness2_FJ_CA->Fill(fjNSubJettiness(Jet1,0,2,1,FJ_Beta,0));
1111  fhSubJettiness2_FJ_AKT->Fill(fjNSubJettiness(Jet1,0,2,2,FJ_Beta,0)); //because in this case we aren't garantueed to get 2 subjets
1112  fhSubJettiness2_FJ_WTA_KT->Fill(fjNSubJettiness(Jet1,0,2,3,FJ_Beta,0));
1113  fhSubJettiness2_FJ_WTA_CA->Fill(fjNSubJettiness(Jet1,0,2,4,FJ_Beta,0));
1114  fhSubJettiness2_FJ_OP_KT->Fill(fjNSubJettiness(Jet1,0,2,5,FJ_Beta,0));
1115  fhSubJettiness2_FJ_OP_CA->Fill(fjNSubJettiness(Jet1,0,2,6,FJ_Beta,0));
1116  fhSubJettiness2_FJ_OP_AKT->Fill(fjNSubJettiness(Jet1,0,2,7,FJ_Beta,0));
1117  fhSubJettiness2_FJ_OP_WTA_KT->Fill(fjNSubJettiness(Jet1,0,2,8,FJ_Beta,0));
1118  fhSubJettiness2_FJ_OP_WTA_CA->Fill(fjNSubJettiness(Jet1,0,2,9,FJ_Beta,0));
1119  fhSubJettiness2_FJ_MIN->Fill(fjNSubJettiness(Jet1,0,2,10,FJ_Beta,0));
1120 
1121  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) fShapesVar[0]= Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
1122  else fShapesVar[0]=Jet1->Pt();
1123  fShapesVar[2]=fjNSubJettiness(Jet1,0,1,0,1,0);
1124  fShapesVar[4]=fjNSubJettiness(Jet1,0,2,0,1,0);
1125  fShapesVar[6]=fjNSubJettiness(Jet1,0,3,0,1,0);
1126  fShapesVar[8]=fjNSubJettiness(Jet1,0,2,0,1,1);
1127  fShapesVar[10]=Jet1->GetNumberOfTracks();
1128  AliEmcalJetFinder *Reclusterer1 = Recluster(Jet1, 0, fSubJetRadius, fSubJetMinPt, fSubJetAlgorithm, "SubJetFinder");
1129  if (fFullTree){
1130  fShapesVar[12]=fjNSubJettiness(Jet1,0,2,0,1,2);
1131  fShapesVar[14]=SubJetFraction(Jet1, Reclusterer1, 1, 0, kTRUE, kFALSE);
1132  fShapesVar[16]=SubJetFraction(Jet1, Reclusterer1, 2, 0, kTRUE, kFALSE);
1133  }
1134  fShapesVar[1]=-2;
1135  fShapesVar[3]=-2;
1136  fShapesVar[5]=-2;
1137  fShapesVar[7]=-2;
1138  fShapesVar[9]=-2;
1139  fShapesVar[11]=-2;
1140  if (fFullTree){
1141  fShapesVar[13]=-2;
1142  fShapesVar[15]=-2;
1143  fShapesVar[17]=-2;
1144  }
1145  fTreeResponseMatrixAxis->Fill();
1146 
1147  fhSubJetCounter->Fill(Reclusterer1->GetNumberOfJets());
1148  for (Int_t i= 0; i<Reclusterer1->GetNumberOfJets(); i++){
1149  fhEventCounter->Fill(6); //Number of overall subjets in all events
1150  fhSubJetPt->Fill(Reclusterer1->GetJet(i)->Pt());
1151  fhSubJetMass->Fill(Reclusterer1->GetJet(i)->M());
1152  fhNumberOfSubJetTracks->Fill(Reclusterer1->GetJet(i)->GetNumberOfTracks());
1153  fhSubJetRadius->Fill(TMath::Sqrt((Reclusterer1->GetJet(i)->Area()/TMath::Pi()))); //Radius of SubJets per event
1154  }
1155  }
1156  }
1157  fhJetCounter->Fill(JetCounter); //Number of Jets in Each Event
1158  }
1159  //else {fhEventCounter->Fill(2);} //Events with no jet container
1160  }
1161  return kTRUE;
1162 }
1163 
1164 //________________________________________________________________________
1165 Double_t AliAnalysisTaskSubJetFraction::RelativePhi(Double_t Phi1, Double_t Phi2){
1166 
1167  if(Phi1 < -1*TMath::Pi()) Phi1 += (2*TMath::Pi()); // Turns the range of 0to2Pi into -PitoPi ???????????
1168  else if (Phi1 > TMath::Pi()) Phi1 -= (2*TMath::Pi());
1169  if(Phi2 < -1*TMath::Pi()) Phi2 += (2*TMath::Pi());
1170  else if (Phi2 > TMath::Pi()) Phi2 -= (2*TMath::Pi());
1171  Double_t DeltaPhi=Phi2-Phi1;
1172  if(DeltaPhi < -1*TMath::Pi()) DeltaPhi += (2*TMath::Pi());
1173  else if (DeltaPhi > TMath::Pi()) DeltaPhi -= (2*TMath::Pi());
1174  return DeltaPhi;
1175 }
1176 
1177 
1178 
1179 //--------------------------------------------------------------------------
1180 Double_t AliAnalysisTaskSubJetFraction::Angularity(AliEmcalJet *Jet, Int_t JetContNb){
1181 
1182  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1183  Double_t Angularity_Numerator=0; //Reset these values
1184  Double_t Angularity_Denominator=0;
1185  AliVParticle *Particle=0x0;
1186  Double_t DeltaPhi=0;
1187 
1188 
1189  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks (particles in the jet
1190  Particle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
1191  if(!Particle) continue;
1192  DeltaPhi=RelativePhi(Jet->Phi(),Particle->Phi());
1193  Angularity_Numerator=Angularity_Numerator+(Particle->Pt()*TMath::Sqrt(((Particle->Eta()-Jet->Eta())*(Particle->Eta()-Jet->Eta()))+(DeltaPhi*DeltaPhi)));
1194  Angularity_Denominator= Angularity_Denominator+Particle->Pt();
1195  }
1196  if(Angularity_Denominator!=0) return Angularity_Numerator/Angularity_Denominator;
1197  else return -1;
1198 
1199 }
1200 
1201 
1202 
1203 //--------------------------------------------------------------------------
1204 Double_t AliAnalysisTaskSubJetFraction::PTD(AliEmcalJet *Jet, Int_t JetContNb){
1205 
1206  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1207  Double_t PTD_Numerator=0; //Reset these values
1208  Double_t PTD_Denominator=0;
1209  AliVParticle *Particle=0x0;
1210  Double_t DeltaPhi=0;
1211  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks (particles in the jet
1212  Particle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
1213  if(!Particle) continue;
1214  PTD_Numerator=PTD_Numerator+(Particle->Pt()*Particle->Pt());
1215  PTD_Denominator=PTD_Denominator+Particle->Pt();
1216  }
1217  if(PTD_Denominator!=0) return TMath::Sqrt(PTD_Numerator)/PTD_Denominator;
1218  else return -1;
1219 
1220 }
1221 
1222 
1223 //----------------------------------------------------------------------
1225 AliEmcalJetFinder *AliAnalysisTaskSubJetFraction::Recluster(AliEmcalJet *Jet, Int_t JetContNb, Double_t SubJetRadius, Double_t SubJetMinPt, Int_t Algorithm, const char* Name){
1226 
1227  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1228  AliEmcalJetFinder *Reclusterer = new AliEmcalJetFinder(Name); //JetFinder Object for reclustered jets
1229  Reclusterer->SetRadius(SubJetRadius);
1230  Reclusterer->SetJetMinPt(SubJetMinPt);
1231  Reclusterer->SetJetAlgorithm(Algorithm); //0 for anti-kt 1 for kt
1232  Reclusterer->SetJetMaxEta(0.9);
1233  Reclusterer->SetRecombSheme(0);
1234  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
1235  Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
1236  if(Reclusterer->AliEmcalJetFinder::Filter(Jet, JetCont, dVtx)){;} //reclustering jet1 using the jetfinderobject Reclusterer
1237  return Reclusterer;
1238 }
1239 
1240 
1241 
1242 
1243 
1244 
1245 
1246 
1247 
1248 //----------------------------------------------------------------------
1249 Double_t AliAnalysisTaskSubJetFraction::SubJetOrdering(AliEmcalJet *Jet, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t Type, Bool_t Index){
1250  AliEmcalJet *SubJet=NULL;
1251  Double_t SortingVariable;
1252  Int_t ArraySize =N+1;
1253  TArrayD *JetSorter = new TArrayD(ArraySize);
1254  TArrayD *JetIndexSorter = new TArrayD(ArraySize);
1255  for (Int_t i=0; i<ArraySize; i++){
1256  JetSorter->SetAt(0,i);
1257  }
1258  for (Int_t i=0; i<ArraySize; i++){
1259  JetIndexSorter->SetAt(0,i);
1260  }
1261  if(Reclusterer->GetNumberOfJets()<N) return -999;
1262  for (Int_t i=0; i<Reclusterer->GetNumberOfJets(); i++){
1263  SubJet=Reclusterer->GetJet(i);
1264  if (Type==0) SortingVariable=SubJet->Pt();
1265  else if (Type==1) SortingVariable=SubJet->E();
1266  else if (Type==2) SortingVariable=SubJet->M();
1267  for (Int_t j=0; j<N; j++){
1268  if (SortingVariable>JetSorter->GetAt(j)){
1269  for (Int_t k=N-1; k>=j; k--){
1270  JetSorter->SetAt(JetSorter->GetAt(k),k+1);
1271  JetIndexSorter->SetAt(JetIndexSorter->GetAt(k),k+1);
1272  }
1273  JetSorter->SetAt(SortingVariable,j);
1274  JetIndexSorter->SetAt(i,j);
1275  break;
1276  }
1277  }
1278  }
1279  if (!Index) return JetSorter->GetAt(N-1);
1280  else return JetIndexSorter->GetAt(N-1);
1281 }
1282 
1283 
1284 
1285 //returns -1 if the Nth hardest jet is requested where N>number of available jets
1286 //type: 0=Pt 1=E 2=M
1287 //Index TRUE=returns index FALSE=returns value of quantatiy in question
1288 
1289 
1290 
1291 
1292 
1293 
1294 
1295 //----------------------------------------------------------------------
1296 Double_t AliAnalysisTaskSubJetFraction::SubJetFraction(AliEmcalJet *Jet, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t Type, Bool_t Add, Bool_t Loss){
1297  AliEmcalJet *SubJet=NULL;
1298  Double_t Observable=0;
1299  Double_t Fraction_Numerator=0;
1300  Bool_t Error=kFALSE;
1301  if (!Jet->GetNumberOfTracks()) return -2;
1302  if (Add){
1303  for (Int_t i=1; i<=N; i++){
1304  Observable=SubJetOrdering(Jet,Reclusterer,i,Type,kFALSE);
1305  if(Observable==-999){
1306  Error = kTRUE;
1307  return -2;
1308  }
1309  Fraction_Numerator=Fraction_Numerator+Observable;
1310  }
1311  }
1312  else {
1313  Fraction_Numerator=SubJetOrdering(Jet,Reclusterer,N,Type,kFALSE);
1314  if (Fraction_Numerator==-999) return -2;
1315  }
1316  if (Type==0){
1317  if(Loss) return (Jet->Pt()-Fraction_Numerator)/Jet->Pt();
1318  else return Fraction_Numerator/Jet->Pt();
1319  }
1320  else if (Type==1){
1321  if(Loss) return (Jet->E()-Fraction_Numerator)/Jet->E();
1322  else return Fraction_Numerator/Jet->E();
1323  }
1324  else { //change to else if if you want to add more later
1325  if(Loss) return (Jet->M()-Fraction_Numerator)/Jet->M();
1326  else return Fraction_Numerator/Jet->M();
1327  }
1328 }
1329 //N number of hardest subjets involved
1330 //Type 0=Pt 1=E 2=M
1331 // Add TRUE: Add 1 to Nth hardest subjets together/total jet False: Nth hardest Jet/toal jet
1332 //Loss TRUE: Jet-Subjet(s)/Jet FALSE: Subjet(s)/Jet
1333 
1334 
1335 //----------------------------------------------------------------------
1336 Double_t AliAnalysisTaskSubJetFraction::NSubJettiness(AliEmcalJet *Jet, Int_t JetContNb, Double_t JetRadius, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t A, Int_t B){
1337  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1338  AliEmcalJet *SubJet=NULL;
1339  Double_t DeltaR1=0;
1340  Double_t DeltaR2=0;
1341  AliVParticle *JetParticle=0x0;
1342  Double_t SubJetiness_Numerator = 0;
1343  Double_t SubJetiness_Denominator = 0;
1344  Double_t Index=-2;
1345  Bool_t Error=kFALSE;
1346  // JetRadius=TMath::Sqrt((Jet->Area()/TMath::Pi())); //comment out later
1347  if (!Jet->GetNumberOfTracks()) return -2;
1348  if (Reclusterer->GetNumberOfJets() < N) return -2;
1349  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks (particles in the jet
1350  JetParticle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
1351  for (Int_t j=1; j<=N; j++){
1352  Index=SubJetOrdering(Jet,Reclusterer,j,0,kTRUE);
1353  if(Index==-999){
1354  Error = kTRUE;
1355  i=Jet->GetNumberOfTracks();
1356  break;
1357  }
1358  if(j==1){
1359  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);
1360  }
1361  else{
1362  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);
1363  if (DeltaR2<DeltaR1) DeltaR1=DeltaR2;
1364  }
1365  }
1366  SubJetiness_Numerator=SubJetiness_Numerator+(JetParticle->Pt()*DeltaR1);
1367  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));
1368  else SubJetiness_Denominator=SubJetiness_Denominator+(TMath::Power((Reclusterer->GetJet(SubJetOrdering(Jet,Reclusterer,N,0,kTRUE))->Pt()),A)*JetParticle->Pt()*TMath::Power(JetRadius,B));
1369  }
1370  if (SubJetiness_Denominator!=0 && !Error){
1371  // if (SubJetiness_Numerator/SubJetiness_Denominator>1.0) cout << JetRadius<<endl;
1372  return SubJetiness_Numerator/SubJetiness_Denominator; }
1373  else return -2;
1374 }
1375 
1376 
1377 //______________________________________________________________________________________
1378 Double_t AliAnalysisTaskSubJetFraction::fjNSubJettiness(AliEmcalJet *Jet, Int_t JetContNb,Int_t N, Int_t Algorithm, Double_t Beta, Int_t Option){
1379 
1380  //WARNING!!! Only works for parent jets that are clustered with Anti-Kt! To change go to AliEmcalJetFinder.cxx and look at the Nsubjettiness() function
1381 
1382  //Algorithm==0 -> kt_axes;
1383  // Algorithm==1 -> ca_axes;
1384  //Algorithm==2 -> antikt_0p2_axes;
1385  //Algorithm==3 -> wta_kt_axes;
1386  //Algorithm==4 -> wta_ca_axes;
1387  //Algorithm==5 -> onepass_kt_axes;
1388  //Algorithm==6 -> onepass_ca_axes;
1389  //Algorithm==7 -> onepass_antikt_0p2_axes;
1390  //Algorithm==8 -> onepass_wta_kt_axes;
1391  //Algorithm==9 -> onepass_wta_ca_axes;
1392  //Algorithm==10 -> min_axes;
1393 
1394 
1395  //Option==0 returns Nsubjettiness Value
1396  //Option==1 && N==2 returns opening angle between two subjet axes(Delta R?)
1397  //Option==2 && N==2 returns Delta R
1398  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1399  AliEmcalJetFinder *JetFinder=new AliEmcalJetFinder("Nsubjettiness");
1400  if (Jet->GetNumberOfTracks()>=N){
1401  if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==1) && (Algorithm==0) && (Beta==1.0) && (Option==0)){
1402  if (fDerivSubtrOrder == kFirstOrder) return Jet->GetFirstOrderSubtracted1subjettiness_kt();
1403  else return Jet->GetSecondOrderSubtracted1subjettiness_kt();
1404  }
1405  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==0) && (Beta==1.0) && (Option==0)){
1406  if (fDerivSubtrOrder == kFirstOrder) return Jet->GetFirstOrderSubtracted2subjettiness_kt();
1407  else return Jet->GetSecondOrderSubtracted2subjettiness_kt();
1408  }
1409  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==3) && (Algorithm==0) && (Beta==1.0) && (Option==0)){
1410  if (fDerivSubtrOrder == kFirstOrder) return Jet->GetFirstOrderSubtracted3subjettiness_kt();
1411  else return Jet->GetSecondOrderSubtracted3subjettiness_kt();
1412  }
1413  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==0) && (Beta==1.0) && (Option==1)){
1414  if (fDerivSubtrOrder == kFirstOrder) return Jet->GetFirstOrderSubtractedOpeningAngle_kt();
1415  else return Jet->GetSecondOrderSubtractedOpeningAngle_kt();
1416  }
1417  else{
1418  JetFinder->SetJetMaxEta(0.9-fJetRadius);
1419  JetFinder->SetRadius(fJetRadius);
1420  JetFinder->SetJetAlgorithm(0); //0 for anti-kt 1 for kt //this is for the JET!!!!!!!!!! Not the SubJets
1421  JetFinder->SetRecombSheme(0);
1422  JetFinder->SetJetMinPt(Jet->Pt());
1423  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
1424  Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
1425  return JetFinder->Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,fSubJetRadius,Beta,Option);
1426  }
1427  }
1428  else return -2;
1429 }
1430 
1431 
1432 //________________________________________________________________________
1433 Bool_t AliAnalysisTaskSubJetFraction::RetrieveEventObjects() {
1434  //
1435  // retrieve event objects
1436  //
1437  if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
1438  return kFALSE;
1439 
1440  return kTRUE;
1441 }
1442 
1443 
1444 //_______________________________________________________________________
1445 void AliAnalysisTaskSubJetFraction::Terminate(Option_t *)
1446 {
1447  // Called once at the end of the analysis.
1448  if (fJetShapeType==AliAnalysisTaskSubJetFraction::kData){
1449 
1450  /*
1451  for (int i=1; i<=(XBinsJetPtSize)-1; i++){ //Rescales the fhJetPt Histograms according to the with of the variable bins and number of events
1452  fhJetPt->SetBinContent(i,((1.0/(XBinsJetPt[i]-XBinsJetPt[i-1]))*((fhJetPt->GetBinContent(i))*(1.0/(fhEventCounter->GetBinContent(1))))));
1453  fhSubJetPt->SetBinContent(i,((1.0/(XBinsJetPt[i]-XBinsJetPt[i-1]))*((fhSubJetPt->GetBinContent(i))*(1.0/(fhEventCounter->GetBinContent(5))))));
1454 
1455  //fhJetPt->SetBinContent(i,((1.0/(fhPt->GetBinWidth(i)))*((fhJetPt->GetBinContent(i))*(fhEventCounter->GetBinContent(1)))));
1456  // fhJetPt->SetBinContent(i,((1.0/((fhPt->GetLowEdge(i+1))-(fhJetPt->GetLowEdge(i))))*((fhPt->GetBinContent(i))*(fhEventCounter->GetBinContent(1)))));
1457  }
1458  for (int i=1; i<=(XBinsJetMassSize)-1; i++){ //Rescales the fhPt Histograms according to the with of the variable bins and number of events
1459  fhJetMass->SetBinContent(i,((1.0/(XBinsJetMass[i]-XBinsJetMass[i-1]))*((fhJetMass->GetBinContent(i))*(1.0/(fhEventCounter->GetBinContent(1))))));
1460  }
1461 
1462  fhJetPhi->Scale(Phi_Bins/((Phi_Up-Phi_Low)*((fhEventCounter->GetBinContent(1)))));
1463  fhJetEta->Scale(Eta_Bins/((Eta_Up-Eta_Low)*((fhEventCounter->GetBinContent(1)))));
1464  fhJetRadius->Scale(100/(fhEventCounter->GetBinContent(4))); //should this and JetAngularity be divided by Bin 1 or 4????
1465  fhJetAngularity->Scale(100/(fhEventCounter->GetBinContent(4)));
1466  fhJetPTD->Scale(100/(fhEventCounter->GetBinContent(4)));
1467  fhNumberOfJetTracks->Scale(1.0/(fhEventCounter->GetBinContent(4)));
1468  fhSubJetPtFrac->Scale(100/(fhEventCounter->GetBinContent(8)));
1469  fhSubJetPtFrac2->Scale(100/(fhEventCounter->GetBinContent(9)));
1470  fhSubJetEnergyFrac->Scale(100/(fhEventCounter->GetBinContent(8)));
1471  fhSubJetEnergyFrac2->Scale(100/(fhEventCounter->GetBinContent(9)));
1472  fhJetCounter->Scale(1.0/(fhEventCounter->GetBinContent(1))); //is the first bin the correct one to look at?
1473 
1474 
1475  fhSubJetCounter->Scale(1.0/(fhEventCounter->GetBinContent(5)));
1476  fhSubJetiness1->Scale(100/(fhEventCounter->GetBinContent(12)));
1477  fhSubJetiness2->Scale(100/(fhEventCounter->GetBinContent(12)));
1478  fh2to1SubJetinessRatio->Scale(100/(fhEventCounter->GetBinContent(12)));
1479 
1480  */
1481  /*
1482 
1483  fhJetPt->Scale(1.0/(fhEventCounter->GetBinContent(1)));
1484  fhSubJetPt->Scale(1.0/(fhEventCounter->GetBinContent(5)));
1485  fhJetMass->Scale(1.0/(fhEventCounter->GetBinContent(1)));
1486  fhJetPhi->Scale(1.0/(fhEventCounter->GetBinContent(1)));
1487  fhJetEta->Scale(1.0/(fhEventCounter->GetBinContent(1)));
1488  fhJetRadius->Scale(1.0/(fhEventCounter->GetBinContent(4))); //should this and JetAngularity be divided by Bin 1 or 4????
1489  fhJetAngularity->Scale(1.0/(fhEventCounter->GetBinContent(4)));
1490  fhJetPTD->Scale(1.0/(fhEventCounter->GetBinContent(4)));
1491  fhNumberOfJetTracks->Scale(1.0/(fhEventCounter->GetBinContent(4)));
1492  fhSubJetPtFrac->Scale(1.0/(fhEventCounter->GetBinContent(8)));
1493  fhSubJetPtFrac2->Scale(1.0/(fhEventCounter->GetBinContent(9)));
1494  fhSubJetEnergyFrac->Scale(1.0/(fhEventCounter->GetBinContent(8)));
1495  fhSubJetEnergyFrac2->Scale(1.0/(fhEventCounter->GetBinContent(9)));
1496  fhJetCounter->Scale(1.0/(fhEventCounter->GetBinContent(1))); //is the first bin the correct one to look at?
1497  fhSubJetCounter->Scale(1.0/(fhEventCounter->GetBinContent(5)));
1498  fhSubJetiness1->Scale(1.0/(fhEventCounter->GetBinContent(12)));
1499  fhSubJetiness2->Scale(1.0/(fhEventCounter->GetBinContent(12)));
1500  fh2to1SubJetinessRatio->Scale(1.0/(fhEventCounter->GetBinContent(12)));
1501  */
1502  }
1503 
1504 }
AliEmcalJetFinder::SetRadius
void SetRadius(Double_t val)
Definition: AliEmcalJetFinder.h:40
AliEmcalJet::Area
Double_t Area() const
Definition: AliEmcalJet.h:69
AliAnalysisTaskSubJetFraction::fjNSubJettiness
Double_t fjNSubJettiness(AliEmcalJet *Jet, Int_t JetContNb, Int_t N, Int_t Algorithm, Double_t Beta, Int_t Option=0)
Definition: AliAnalysisTaskSubJetFraction.cxx:1378
AliPicoTrack.h
XBinsPtSize
Double_t XBinsPtSize
Definition: AliAnalysisTaskSubJetFraction.cxx:49
AliEmcalJetFinder::SetJetMinPt
void SetJetMinPt(Double_t val)
Definition: AliEmcalJetFinder.h:45
AliEmcalJetFinder.h
AliAnalysisTaskEmcal::DeltaPhi
static Double_t DeltaPhi(Double_t phia, Double_t phib, Double_t rMin=-TMath::Pi()/2, Double_t rMax=3 *TMath::Pi()/2)
Definition: AliAnalysisTaskEmcal.h:336
AliEmcalJet::ClosestJet
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:158
AliAnalysisTaskEmcalJet::GetJetContainer
AliJetContainer * GetJetContainer(Int_t i=0) const
Definition: AliAnalysisTaskEmcalJet.cxx:315
AliEmcalJet::Eta
Double_t Eta() const
Definition: AliEmcalJet.h:60
AliEmcalJet::Phi
Double_t Phi() const
Definition: AliEmcalJet.h:55
AliAnalysisTaskSubJetFraction::PTD
Double_t PTD(AliEmcalJet *Jet, Int_t JetContNb)
Definition: AliAnalysisTaskSubJetFraction.cxx:1204
AliEmcalJet::GetLabel
Int_t GetLabel() const
Definition: AliEmcalJet.h:63
AliEmcalParticle.h
AliEmcalJet::GetSecondOrderSubtracted3subjettiness_kt
Double_t GetSecondOrderSubtracted3subjettiness_kt() const
Definition: AliEmcalJet.h:286
AliAnalysisTaskSubJetFraction::NSubJettiness
Double_t NSubJettiness(AliEmcalJet *Jet, Int_t JetContNb, Double_t JetRadius, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t A, Int_t B)
Definition: AliAnalysisTaskSubJetFraction.cxx:1336
AliEmcalJet::GetFirstOrderSubtractedOpeningAngle_kt
Double_t GetFirstOrderSubtractedOpeningAngle_kt() const
Definition: AliEmcalJet.h:295
AliAnalysisTaskSubJetFraction::RelativePhi
Double_t RelativePhi(Double_t Phi1, Double_t Phi2)
Definition: AliAnalysisTaskSubJetFraction.cxx:1165
AliEmcalJet::E
Double_t E() const
Definition: AliEmcalJet.h:58
AliAnalysisTaskEmcal::fOutput
TList * fOutput
!output list
Definition: AliAnalysisTaskEmcal.h:304
Eta_Bins
Int_t Eta_Bins
Definition: AliAnalysisTaskSubJetFraction.cxx:56
XBinsJetMass
Double_t XBinsJetMass[28]
Definition: AliAnalysisTaskSubJetFraction.cxx:52
AliEmcalJet::GetFirstOrderSubtracted2subjettiness_kt
Double_t GetFirstOrderSubtracted2subjettiness_kt() const
Definition: AliEmcalJet.h:275
AliEmcalJet::GetNumberOfTracks
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:83
kData
TString kData
Declare data MC or deltaAOD.
Definition: AddTaskCaloTrackCorr.C:21
AliJetContainer::GetParticleContainer
AliParticleContainer * GetParticleContainer() const
Definition: AliJetContainer.h:156
AliEmcalJetFinder::SetRecombSheme
void SetRecombSheme(Int_t val)
Definition: AliEmcalJetFinder.h:42
Eta_Up
Double_t Eta_Up
Definition: AliAnalysisTaskSubJetFraction.cxx:54
AliEmcalJetFinder::SetJetAlgorithm
void SetJetAlgorithm(Int_t val)
Definition: AliEmcalJetFinder.h:41
ClassImp
ClassImp(AliAnalysisTaskSubJetFraction) AliAnalysisTaskSubJetFraction
Definition: AliAnalysisTaskSubJetFraction.cxx:68
Eta_Low
Double_t Eta_Low
Definition: AliAnalysisTaskSubJetFraction.cxx:55
AliEmcalJetFinder::Nsubjettiness
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)
Definition: AliEmcalJetFinder.cxx:286
XBinsJetPtSize
Int_t XBinsJetPtSize
Definition: AliAnalysisTaskSubJetFraction.cxx:51
AliEmcalJet::GetFirstOrderSubtracted3subjettiness_kt
Double_t GetFirstOrderSubtracted3subjettiness_kt() const
Definition: AliEmcalJet.h:285
AliParticleContainer.h
AliJetContainer::GetNJets
Int_t GetNJets() const
Definition: AliJetContainer.h:129
Phi_Up
Double_t Phi_Up
Definition: AliAnalysisTaskSubJetFraction.cxx:57
AliJetContainer.h
AliAnalysisTaskEmcal::fCent
Double_t fCent
!event centrality
Definition: AliAnalysisTaskEmcal.h:286
AliEmcalJet::GetSecondOrderSubtracted2subjettiness_kt
Double_t GetSecondOrderSubtracted2subjettiness_kt() const
Definition: AliEmcalJet.h:276
AliEmcalContainer::GetArray
TClonesArray * GetArray() const
Definition: AliEmcalContainer.h:119
XBinsJetMassSize
Int_t XBinsJetMassSize
Definition: AliAnalysisTaskSubJetFraction.cxx:53
AliAnalysisTaskSubJetFraction::SubJetFraction
Double_t SubJetFraction(AliEmcalJet *Jet, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t Type, Bool_t Add, Bool_t Loss)
Definition: AliAnalysisTaskSubJetFraction.cxx:1296
AliRhoParameter.h
AliEmcalJetFinder::SetJetMaxEta
void SetJetMaxEta(Double_t val)
Definition: AliEmcalJetFinder.h:44
Phi_Low
Double_t Phi_Low
Definition: AliAnalysisTaskSubJetFraction.cxx:58
AliJetContainer::GetNextAcceptJet
AliEmcalJet * GetNextAcceptJet()
Definition: AliJetContainer.cxx:307
Data
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)
Definition: charmCutsOptimization.C:371
AliAnalysisTaskSubJetFraction::Recluster
AliEmcalJetFinder * Recluster(AliEmcalJet *Jet, Int_t JetContNb, Double_t SubJetRadius, Double_t SubJetMinPt, Int_t Algorithm, const char *Name)
Definition: AliAnalysisTaskSubJetFraction.cxx:1225
AliEmcalJet::GetFirstOrderSubtracted1subjettiness_kt
Double_t GetFirstOrderSubtracted1subjettiness_kt() const
Definition: AliEmcalJet.h:265
AliEmcalJetFinder::GetJet
AliEmcalJet * GetJet(Int_t index)
Definition: AliEmcalJetFinder.h:37
AliAnalysisTaskSubJetFraction.h
AliEmcalJet::Pt
Double_t Pt() const
Definition: AliEmcalJet.h:47
AliAnalysisTaskEmcalJet::GetRhoVal
Double_t GetRhoVal(Int_t i=0) const
Definition: AliAnalysisTaskEmcalJet.cxx:598
XBinsPt
Double_t XBinsPt[66]
Definition: AliAnalysisTaskSubJetFraction.cxx:48
AliEmcalJet::TrackAt
Short_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:107
AliAnalysisTaskEmcal::SetMakeGeneralHistograms
void SetMakeGeneralHistograms(Bool_t g)
Definition: AliAnalysisTaskEmcal.h:149
AliEmcalJet::GetSecondOrderSubtracted1subjettiness_kt
Double_t GetSecondOrderSubtracted1subjettiness_kt() const
Definition: AliEmcalJet.h:266
AliAnalysisTaskSubJetFraction::Angularity
Double_t Angularity(AliEmcalJet *Jet, Int_t JetContNb)
Definition: AliAnalysisTaskSubJetFraction.cxx:1180
AliAnalysisTaskSubJetFraction::SubJetOrdering
Double_t SubJetOrdering(AliEmcalJet *Jet, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t Type, Bool_t Index)
Definition: AliAnalysisTaskSubJetFraction.cxx:1249
AliEmcalJet::GetSecondOrderSubtractedOpeningAngle_kt
Double_t GetSecondOrderSubtractedOpeningAngle_kt() const
Definition: AliEmcalJet.h:296
XBinsJetPt
Double_t XBinsJetPt[35]
Definition: AliAnalysisTaskSubJetFraction.cxx:50
Phi_Bins
Int_t Phi_Bins
Definition: AliAnalysisTaskSubJetFraction.cxx:59
AliEmcalJet::M
Double_t M() const
Definition: AliEmcalJet.h:59
AliEmcalContainer::ResetCurrentID
void ResetCurrentID(Int_t i=-1)
Definition: AliEmcalContainer.h:141
AliEmcalJet.h
AliJetContainer::GetJet
AliEmcalJet * GetJet(Int_t i) const
Definition: AliJetContainer.cxx:266