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