AliPhysics  8bb951a (8bb951a)
 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 "AliEmcalPythiaInfo.h"
39 #include "TRandom3.h"
40 #include "AliPicoTrack.h"
41 #include "AliEmcalJetFinder.h"
42 #include "AliAODEvent.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
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 
69 
70 //________________________________________________________________________
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  fhJetPt(0x0),
97  fhJetPt_1(0x0),
98  fhJetPt_2(0x0),
99  fhJetPhi(0x0),
100  fhJetPhi_1(0x0),
101  fhJetPhi_2(0x0),
102  fhJetEta(0x0),
103  fhJetEta_1(0x0),
104  fhJetEta_2(0x0),
105  fhJetMass(0x0),
106  fhJetMass_1(0x0),
107  fhJetMass_2(0x0),
108  fhJetRadius(0x0),
109  fhJetRadius_1(0x0),
110  fhJetRadius_2(0x0),
111  fhJetAngularity(0x0),
112  fhJetAngularity_1(0x0),
113  fhJetAngularity_2(0x0),
114  fhJetAngularityJetPt(0x0),
115  fhJetAngularityJetPt_1(0x0),
116  fhJetAngularityJetPt_2(0x0),
117  fhJetPTD(0x0),
118  fhJetPTD_1(0x0),
119  fhJetPTD_2(0x0),
120  fhJetPTDJetPt(0x0),
121  fhJetPTDJetPt_1(0x0),
122  fhJetPTDJetPt_2(0x0),
123  fhJetCounter(0x0),
124  fhJetCounter_1(0x0),
125  fhJetCounter_2(0x0),
126  fhNumberOfJetTracks(0x0),
127  fhNumberOfJetTracks_1(0x0),
128  fhNumberOfJetTracks_2(0x0),
129  fhSubJetPt(0x0),
130  fhSubJetPt_1(0x0),
131  fhSubJetPt_2(0x0),
132  fhSubJetMass(0x0),
133  fhSubJetMass_1(0x0),
134  fhSubJetMass_2(0x0),
135  fhSubJetRadius(0x0),
136  fhSubJetRadius_1(0x0),
137  fhSubJetRadius_2(0x0),
138  fhSubJetCounter(0x0),
139  fhSubJetCounter_1(0x0),
140  fhSubJetCounter_2(0x0),
141  fhNumberOfSubJetTracks(0x0),
142  fhNumberOfSubJetTracks_1(0x0),
143  fhNumberOfSubJetTracks_2(0x0),
144  fhSubJetPtFrac(0x0),
145  fhSubJetPtFrac_1(0x0),
146  fhSubJetPtFrac_2(0x0),
147  fhSubJetPtFrac2(0x0),
148  fhSubJetPtFrac2_1(0x0),
149  fhSubJetPtFrac2_2(0x0),
150  fhSubJetPtLoss(0x0),
151  fhSubJetPtLoss_1(0x0),
152  fhSubJetPtLoss_2(0x0),
153  fhSubJetPtLoss2(0x0),
154  fhSubJetPtLoss2_1(0x0),
155  fhSubJetPtLoss2_2(0x0),
156  fhSubJetEnergyFrac(0x0),
157  fhSubJetEnergyFrac_1(0x0),
158  fhSubJetEnergyFrac_2(0x0),
159  fhSubJetEnergyFrac2(0x0),
160  fhSubJetEnergyFrac2_1(0x0),
161  fhSubJetEnergyFrac2_2(0x0),
162  fhSubJetEnergyLoss(0x0),
163  fhSubJetEnergyLoss_1(0x0),
164  fhSubJetEnergyLoss_2(0x0),
165  fhSubJetEnergyLoss2(0x0),
166  fhSubJetEnergyLoss2_1(0x0),
167  fhSubJetEnergyLoss2_2(0x0),
168  fh2PtRatio(0x0),
169  fhSubJetiness1(0x0),
170  fhSubJetiness1_1(0x0),
171  fhSubJetiness1_2(0x0),
172  fhSubJetiness1JetPt(0x0),
173  fhSubJetiness1JetPt_1(0x0),
174  fhSubJetiness1JetPt_2(0x0),
175  fhSubJetiness2(0x0),
176  fhSubJetiness2_1(0x0),
177  fhSubJetiness2_2(0x0),
178  fhSubJetiness2JetPt(0x0),
179  fhSubJetiness2JetPt_1(0x0),
180  fhSubJetiness2JetPt_2(0x0),
181  fh2to1SubJetinessRatio(0x0),
182  fh2to1SubJetinessRatio_1(0x0),
183  fh2to1SubJetinessRatio_2(0x0),
184  fh2to1SubJetinessRatioJetPt(0x0),
185  fh2to1SubJetinessRatioJetPt_1(0x0),
186  fh2to1SubJetinessRatioJetPt_2(0x0),
187  fhEventCounter(0x0),
188  fhEventCounter_1(0x0),
189  fhEventCounter_2(0x0),
190  fhJetPtJetEta(0x0),
191  fhJetPtJetEta_1(0x0),
192  fhJetPtJetEta_2(0x0),
193  fhSubJetiness2Distance(0x0),
194  fhSubJetiness2Distance_1(0x0),
195  fhSubJetiness2Distance_2(0x0),
196  fTreeResponseMatrixAxis(0)
197 
198 {
199  SetMakeGeneralHistograms(kTRUE);
200 }
201 
202 //________________________________________________________________________
204  AliAnalysisTaskEmcalJet(name, kTRUE),
205  fContainer(0),
206  fMinFractionShared(0),
207  fJetShapeType(kData),
208  fJetShapeSub(kNoSub),
209  fJetSelection(kInclusive),
210  fShapesVar(0),
211  fPtThreshold(-9999.),
212  fRMatching(0.2),
213  fminpTTrig(20.),
214  fmaxpTTrig(50.),
215  fangWindowRecoil(0.6),
216  fSemigoodCorrect(0),
217  fHolePos(0),
218  fHoleWidth(0),
219  fCentSelectOn(kTRUE),
220  fCentMin(0),
221  fCentMax(10),
222  fSubJetAlgorithm(0),
223  fSubJetRadius(0.1),
224  fSubJetMinPt(1),
225  fJetRadius(0.4),
226  fRMatched(0.2),
227  fSharedFractionPtMin(0.5),
228  fhJetPt(0x0),
229  fhJetPt_1(0x0),
230  fhJetPt_2(0x0),
231  fhJetPhi(0x0),
232  fhJetPhi_1(0x0),
233  fhJetPhi_2(0x0),
234  fhJetEta(0x0),
235  fhJetEta_1(0x0),
236  fhJetEta_2(0x0),
237  fhJetMass(0x0),
238  fhJetMass_1(0x0),
239  fhJetMass_2(0x0),
240  fhJetRadius(0x0),
241  fhJetRadius_1(0x0),
242  fhJetRadius_2(0x0),
243  fhJetAngularity(0x0),
244  fhJetAngularity_1(0x0),
245  fhJetAngularity_2(0x0),
246  fhJetAngularityJetPt(0x0),
247  fhJetAngularityJetPt_1(0x0),
248  fhJetAngularityJetPt_2(0x0),
249  fhJetPTD(0x0),
250  fhJetPTD_1(0x0),
251  fhJetPTD_2(0x0),
252  fhJetPTDJetPt(0x0),
253  fhJetPTDJetPt_1(0x0),
254  fhJetPTDJetPt_2(0x0),
255  fhJetCounter(0x0),
256  fhJetCounter_1(0x0),
257  fhJetCounter_2(0x0),
258  fhNumberOfJetTracks(0x0),
259  fhNumberOfJetTracks_1(0x0),
260  fhNumberOfJetTracks_2(0x0),
261  fhSubJetPt(0x0),
262  fhSubJetPt_1(0x0),
263  fhSubJetPt_2(0x0),
264  fhSubJetMass(0x0),
265  fhSubJetMass_1(0x0),
266  fhSubJetMass_2(0x0),
267  fhSubJetRadius(0x0),
268  fhSubJetRadius_1(0x0),
269  fhSubJetRadius_2(0x0),
270  fhSubJetCounter(0x0),
271  fhSubJetCounter_1(0x0),
272  fhSubJetCounter_2(0x0),
273  fhNumberOfSubJetTracks(0x0),
274  fhNumberOfSubJetTracks_1(0x0),
275  fhNumberOfSubJetTracks_2(0x0),
276  fhSubJetPtFrac(0x0),
277  fhSubJetPtFrac_1(0x0),
278  fhSubJetPtFrac_2(0x0),
279  fhSubJetPtFrac2(0x0),
280  fhSubJetPtFrac2_1(0x0),
281  fhSubJetPtFrac2_2(0x0),
282  fhSubJetPtLoss(0x0),
283  fhSubJetPtLoss_1(0x0),
284  fhSubJetPtLoss_2(0x0),
285  fhSubJetPtLoss2(0x0),
286  fhSubJetPtLoss2_1(0x0),
287  fhSubJetPtLoss2_2(0x0),
288  fhSubJetEnergyFrac(0x0),
289  fhSubJetEnergyFrac_1(0x0),
290  fhSubJetEnergyFrac_2(0x0),
291  fhSubJetEnergyFrac2(0x0),
292  fhSubJetEnergyFrac2_1(0x0),
293  fhSubJetEnergyFrac2_2(0x0),
294  fhSubJetEnergyLoss(0x0),
295  fhSubJetEnergyLoss_1(0x0),
296  fhSubJetEnergyLoss_2(0x0),
297  fhSubJetEnergyLoss2(0x0),
298  fhSubJetEnergyLoss2_1(0x0),
299  fhSubJetEnergyLoss2_2(0x0),
300  fh2PtRatio(0x0),
301  fhSubJetiness1(0x0),
302  fhSubJetiness1_1(0x0),
303  fhSubJetiness1_2(0x0),
304  fhSubJetiness1JetPt(0x0),
305  fhSubJetiness1JetPt_1(0x0),
306  fhSubJetiness1JetPt_2(0x0),
307  fhSubJetiness2(0x0),
308  fhSubJetiness2_1(0x0),
309  fhSubJetiness2_2(0x0),
310  fhSubJetiness2JetPt(0x0),
311  fhSubJetiness2JetPt_1(0x0),
312  fhSubJetiness2JetPt_2(0x0),
313  fh2to1SubJetinessRatio(0x0),
314  fh2to1SubJetinessRatio_1(0x0),
315  fh2to1SubJetinessRatio_2(0x0),
316  fh2to1SubJetinessRatioJetPt(0x0),
317  fh2to1SubJetinessRatioJetPt_1(0x0),
318  fh2to1SubJetinessRatioJetPt_2(0x0),
319  fhEventCounter(0x0),
320  fhEventCounter_1(0x0),
321  fhEventCounter_2(0x0),
322  fhJetPtJetEta(0x0),
323  fhJetPtJetEta_1(0x0),
324  fhJetPtJetEta_2(0x0),
325  fhSubJetiness2Distance(0x0),
326  fhSubJetiness2Distance_1(0x0),
327  fhSubJetiness2Distance_2(0x0),
328  fTreeResponseMatrixAxis(0)
329 
330 {
331  // Standard constructor.
332 
334 
335  DefineOutput(1, TTree::Class());
336 
337 }
338 
339 //________________________________________________________________________
341 {
342  // Destructor.
343 }
344 
345 //________________________________________________________________________
347 {
348  // Create user output.
349 
351 
352  Bool_t oldStatus = TH1::AddDirectoryStatus();
353  TH1::AddDirectory(kFALSE);
354 
355  //create a tree used for the MC data and making a 4D response matrix
356  const Int_t nVar = 14;
357  fShapesVar = new Double_t [nVar]; //shapes used for tagging
358  fTreeResponseMatrixAxis = new TTree("fTreeJetShape", "fTreeJetShape");
359  TString *fShapesVarNames = new TString [nVar];
360 
361  fShapesVarNames[0] = "Pt_Data_or_Detector_Level_or_Embeded";
362  fShapesVarNames[1] = "Pt_Particle_Level";
363  fShapesVarNames[2] = "Frac_Data_or_Detector_Level_or_Embeded";
364  fShapesVarNames[3] = "Frac_Particle_Level";
365  fShapesVarNames[4] = "Frac2_Data_or_Detector_Level_or_Embeded";
366  fShapesVarNames[5] = "Frac2_Particle_Level";
367  fShapesVarNames[6] = "1SubJettiness_Data_or_Detector_Level_or_Embeded";
368  fShapesVarNames[7] = "1SubJettiness_Particle_Level";
369  fShapesVarNames[8] = "2SubJettiness_Data_or_Detector_Level_or_Embeded";
370  fShapesVarNames[9] = "2SubJettiness_Particle_Level";
371  fShapesVarNames[10] = "Distance_Between_Two_Hardest_Subjets_Axis_Data_or_Detector_Level_or_Embeded";
372  fShapesVarNames[11] = "Distance_Between_Two_Hardest_Subjets_Axis_Particle_Level";
373  fShapesVarNames[12]= "Parton_Flavour_Data_or_Detector_Level_or_Embeded";
374  fShapesVarNames[13]= "Parton_Flavour_Particle_Level";
375  for(Int_t ivar=0; ivar < nVar; ivar++){
376  cout<<"looping over variables"<<endl;
377  fTreeResponseMatrixAxis->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/D", fShapesVarNames[ivar].Data()));
378  }
379 
380 
381 
383 
384  // fhJetPt= new TH1F("fhJetPt", "Jet Pt", (XBinsJetPtSize)-1, XBinsJetPt);
385  fhJetPt= new TH1F("fhJetPt", "Jet Pt",1500,-0.5,149.5 );
386  fOutput->Add(fhJetPt);
387  //fhJetPhi= new TH1F("fhJetPhi", "Jet Phi", Phi_Bins, Phi_Low, Phi_Up);
388  fhJetPhi= new TH1F("fhJetPhi", "Jet Phi",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
389  fOutput->Add(fhJetPhi);
390  fhJetEta= new TH1F("fhJetEta", "Jet Eta", Eta_Bins, Eta_Low, Eta_Up);
391  fOutput->Add(fhJetEta);
392  fhJetMass= new TH1F("fhJetMass", "Jet Mass", 4000,-0.5, 39.5);
393  fOutput->Add(fhJetMass);
394  fhJetRadius= new TH1F("fhJetRadius", "Jet Radius", 100, -0.05,0.995);
395  fOutput->Add(fhJetRadius);
396  fhJetAngularity= new TH1F("fhJetAngularity", "Jet Angularity", 400, -2,2);
397  fOutput->Add(fhJetAngularity);
398  fhJetAngularityJetPt= new TH2F("fhJetAngularityJetPt", "Jet Angularity vs Jet Pt", 1500, -0.5, 149.50, 400, -2,2);
400  fhJetPTD= new TH1F("fhJetPTD", "Jet PTD", 400, -2,2);
401  fOutput->Add(fhJetPTD);
402  fhJetPTDJetPt= new TH2F("fhJetPTDJetPt", "Jet PTD vs Jet Pt", 1500, -0.5, 149.50, 400, -2,2);
403  fOutput->Add(fhJetPTDJetPt);
404  fhNumberOfJetTracks= new TH1F("fhNumberOfJetTracks", "Number of Tracks within a Jet", 30, -0.5,29.5);
406  fhSubJetRadius= new TH1F("fhSubJetRadius", "SubJet Radius", 100, -0.05,0.995);
407  fOutput->Add(fhSubJetRadius);
408  fhSubJetPt= new TH1F("fhSubJetPt", "SubJet Pt", 1500, -0.5,149.5);
409  fOutput->Add(fhSubJetPt);
410  fhSubJetMass= new TH1F("fhSubJetMass", "Sub Jet Mass", 4000,-0.5, 39.5);
411  fOutput->Add(fhSubJetMass);
412  fhNumberOfSubJetTracks= new TH1F("fhNumberOfSubJetTracks", "Number of Tracks within a Sub Jet", 30, -0.5,29.5);
414  fhSubJetPtFrac= new TH1F("fhSubJetPtFrac", "Pt Fraction of Highest Pt Subjet compared to original Jet",101, -0.05,1.05);
415  fOutput->Add(fhSubJetPtFrac);
416  fhSubJetPtFrac2= new TH1F("fhSubJetPtFrac2", "Pt Fraction of Two Highest Pt Subjets compared to original Jet",101, -0.05,1.05);
417  fOutput->Add(fhSubJetPtFrac2);
418  fhSubJetPtLoss= new TH1F("fhSubJetPtLoss", "Pt Difference of Highest Pt Subjet compared to original Jet",101, -0.05,1.05);
419  fOutput->Add(fhSubJetPtLoss);
420  fhSubJetPtLoss2= new TH1F("fhSubJetPtLoss2", "Pt Difference of Two Highest Pt Subjets compared to original Jet",101, -0.05,1.05);
421  fOutput->Add(fhSubJetPtLoss2);
422  fhSubJetEnergyFrac= new TH1F("fhSubJetEnergyFrac", "Energy Fraction of Most Energetic Subjet compared to original Jet",101, -0.05,1.05);
424  fhSubJetEnergyFrac2= new TH1F("fhSubJetEnergyFrac2", "Energy Fraction of Two Most Energetic Subjets compared to original Jet",101, -0.05,1.05);
426  fhSubJetEnergyLoss= new TH1F("fhSubJetEnergyLoss", "Energy Difference of Most Energetic Subjet compared to original Jet",101, -0.05,1.05);
428  fhSubJetEnergyLoss2= new TH1F("fhSubJetEnergyLoss2", "Pt Difference of Two Most Energetic Subjets compared to original Jet",101, -0.05,1.05);
430  fhSubJetiness1= new TH1F("fhSubjetiness1", "Tau 1 value",101, -0.05,1.05);
431  fOutput->Add(fhSubJetiness1);
432  fhSubJetiness1JetPt= new TH2F("fhSubjetiness1JetPt", "Tau 1 value vs Jet Pt",1500, -0.5, 149.5, 101, -0.05,1.05);
434  fhSubJetiness2= new TH1F("fhSubjetiness2", "Tau 2 value",101, -0.05,1.05);
435  fOutput->Add(fhSubJetiness2);
436  fhSubJetiness2JetPt= new TH2F("fhSubjetiness2JetPt", "Tau 2 value vs Jet Pt",1500, -0.5, 149.5, 101, -0.05,1.05);
438  fhJetCounter= new TH1F("fhJetCounter", "Jet Counter", 100, -0.5, 99.5);
439  fOutput->Add(fhJetCounter);
440  fhSubJetCounter = new TH1F("fhSubJetCounter", "SubJet Counter",50, -0.5,49.5);
441  fOutput->Add(fhSubJetCounter);
442  fh2to1SubJetinessRatio= new TH1F("fh2to1SubJetinessRatio", "Ratio of #tau 1 to #tau 2",200, -0.5,1.5);
444  fh2to1SubJetinessRatioJetPt= new TH2F("fh2to1SubJetinessRatioJetPt", "Ratio of #tau 1 to #tau 2 vs Jet Pt", 1500, -0.5, 149.5, 200, -0.5,1.5);
446  fhJetPtJetEta = new TH2F("fhJetPtJetEta", "Jet Pt vs Jet Eta", 1500,-0.5,150,Eta_Bins, Eta_Low, Eta_Up);
447  fOutput->Add(fhJetPtJetEta);
448  fhSubJetiness2Distance = new TH2F("fhSubJetiness2Distance", "#tau 2 as a function of distance between subject axis",100,0,10,101,-0.05,1.05);
450  fhEventCounter= new TH1F("fhEventCounter", "Event Counter", 15,0.5,15.5);
451  fOutput->Add(fhEventCounter);
452  }
454  fhJetPt_1= new TH1F("fhJetPt_1", "Jet Pt Detector Level",1500,-0.5,149.5 );
455  fOutput->Add(fhJetPt_1);
456  fhJetPt_2= new TH1F("fhJetPt_2", "Jet Pt Particle Level",1500,-0.5,149.5 );
457  fOutput->Add(fhJetPt_2);
458  fhJetPhi_1= new TH1F("fhJetPhi_1", "Jet Phi Detector Level",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
459  fOutput->Add(fhJetPhi_1);
460  fhJetPhi_2= new TH1F("fhJetPhi_2", "Jet Phi Particle Level",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
461  fOutput->Add(fhJetPhi_2);
462  fhJetEta_1= new TH1F("fhJetEta_1", "Jet Eta Detector Level", Eta_Bins, Eta_Low, Eta_Up);
463  fOutput->Add(fhJetEta_1);
464  fhJetEta_2= new TH1F("fhJetEta_2", "Jet Eta Particle Level", Eta_Bins, Eta_Low, Eta_Up);
465  fOutput->Add(fhJetEta_2);
466  fhJetMass_1= new TH1F("fhJetMass_1", "Jet Mass Detector Level", 4000,-0.5, 39.5);
467  fOutput->Add(fhJetMass_1);
468  fhJetMass_2= new TH1F("fhJetMass_2", "Jet Mass Particle Level", 4000,-0.5, 39.5);
469  fOutput->Add(fhJetMass_2);
470  fhJetRadius_1= new TH1F("fhJetRadius_1", "Jet Radius Detector Level", 100, -0.05,0.995);
471  fOutput->Add(fhJetRadius_1);
472  fhJetRadius_2= new TH1F("fhJetRadius_2", "Jet Radius Particle Level", 100, -0.05,0.995);
473  fOutput->Add(fhJetRadius_2);
474  fhJetAngularity_1= new TH1F("fhJetAngularity_1", "Jet Angularity Detector Level", 400, -2,2);
476  fhJetAngularity_2= new TH1F("fhJetAngularity_2", "Jet Angularity Particle Level", 400, -2,2);
478  fhJetAngularityJetPt_1= new TH2F("fhJetAngularityJetPt_1", "Jet Angularity vs Jet Pt Detector Level", 1500, -0.5, 149.50, 400, -2,2);
480  fhJetAngularityJetPt_2= new TH2F("fhJetAngularityJetPt_2", "Jet Angularity vs Jet Pt Particle Level", 1500, -0.5, 149.50, 400, -2,2);
482  fhJetPTD_1= new TH1F("fhJetPTD_1", "Jet PTD Detector Level", 400, -2,2);
483  fOutput->Add(fhJetPTD_1);
484  fhJetPTD_2= new TH1F("fhJetPTD_2", "Jet PTD Particle Level", 400, -2,2);
485  fOutput->Add(fhJetPTD_2);
486  fhJetPTDJetPt_1= new TH2F("fhJetPTDJetPt_1", "Jet PTD vs Jet Pt Detector Level", 1500, -0.5, 149.50, 400, -2,2);
487  fOutput->Add(fhJetPTDJetPt_1);
488  fhJetPTDJetPt_2= new TH2F("fhJetPTDJetPt_2", "Jet PTD vs Jet Pt Particle Level", 1500, -0.5, 149.50, 400, -2,2);
489  fOutput->Add(fhJetPTDJetPt_2);
490  fhNumberOfJetTracks_1= new TH1F("fhNumberOfJetTracks_1", "Number of Tracks within a Jet Detector Level", 30, -0.5,29.5);
492  fhNumberOfJetTracks_2= new TH1F("fhNumberOfJetTracks_2", "Number of Tracks within a Jet Particle Level", 30, -0.5,29.5);
494  fhSubJetRadius_1= new TH1F("fhSubJetRadius_1", "SubJet Radius Detector Level", 100, -0.05,0.995);
496  fhSubJetRadius_2= new TH1F("fhSubJetRadius_2", "SubJet Radius Particle Level", 100, -0.05,0.995);
498  fhSubJetPt_1= new TH1F("fhSubJetPt_1", "SubJet Pt Detector Level", 1500, -0.5,149.5);
499  fOutput->Add(fhSubJetPt_1);
500  fhSubJetPt_2= new TH1F("fhSubJetPt_2", "SubJet Pt Particle Level", 1500, -0.5,149.5);
501  fOutput->Add(fhSubJetPt_2);
502  fhSubJetMass_1= new TH1F("fhSubJetMass_1", "Sub Jet Mass Detector Level", 4000,-0.5, 39.5);
503  fOutput->Add(fhSubJetMass_1);
504  fhSubJetMass_2= new TH1F("fhSubJetMass_2", "Sub Jet Mass Particle Level", 4000,-0.5, 39.5);
505  fOutput->Add(fhSubJetMass_2);
506  fhNumberOfSubJetTracks_1= new TH1F("fhNumberOfSubJetTracks_1", "Number of Tracks within a Sub Jet Detector Level", 30, -0.5,29.5);
508  fhNumberOfSubJetTracks_2= new TH1F("fhNumberOfSubJetTracks_2", "Number of Tracks within a Sub Jet Particle Level", 30, -0.5,29.5);
510  fhSubJetPtFrac_1= new TH1F("fhSubJetPtFrac_1", "Pt Fraction of Highest Pt Subjet compared to original Jet Detector Level",101, -0.05,1.05);
512  fhSubJetPtFrac_2= new TH1F("fhSubJetPtFrac_2", "Pt Fraction of Highest Pt Subjet compared to original Jet Particle Level",101, -0.05,1.05);
514  fhSubJetPtFrac2_1= new TH1F("fhSubJetPtFrac2_1", "Pt Fraction of Two Highest Pt Subjets compared to original Jet Detector Level",101, -0.05,1.05);
516  fhSubJetPtFrac2_2= new TH1F("fhSubJetPtFrac2_2", "Pt Fraction of Two Highest Pt Subjets compared to original Jet Particle Level",101, -0.05,1.05);
518  fhSubJetPtLoss_1= new TH1F("fhSubJetPtLoss_1", "Pt Difference of Highest Pt Subjet compared to original Jet Detector Level",101, -0.05,1.05);
520  fhSubJetPtLoss_2= new TH1F("fhSubJetPtLoss_2", "Pt Difference of Highest Pt Subjet compared to original Jet Particle Level",101, -0.05,1.05);
522  fhSubJetPtLoss2_1= new TH1F("fhSubJetPtLoss2_1", "Pt Difference of Two Highest Pt Subjets compared to original Jet Detector Level",101, -0.05,1.05);
524  fhSubJetPtLoss2_2= new TH1F("fhSubJetPtLoss2_2", "Pt Difference of Two Highest Pt Subjets compared to original Jet Particle Level",101, -0.05,1.05);
526  fhSubJetEnergyFrac_1= new TH1F("fhSubJetEnergyFrac_1", "Energy Fraction of Most Energetic Subjet compared to original Jet Detector Level",101, -0.05,1.05);
528  fhSubJetEnergyFrac_2= new TH1F("fhSubJetEnergyFrac_2", "Energy Fraction of Most Energetic Subjet compared to original Jet Particle Level",101, -0.05,1.05);
530  fhSubJetEnergyFrac2_1= new TH1F("fhSubJetEnergyFrac2_1", "Energy Fraction of Two Most Energetic Subjets compared to original Jet Detector Level",101, -0.05,1.05);
532  fhSubJetEnergyFrac2_2= new TH1F("fhSubJetEnergyFrac2_2", "Energy Fraction of Two Most Energetic Subjets compared to original Jet Particle Level",101, -0.05,1.05);
534  fhSubJetEnergyLoss_1= new TH1F("fhSubJetEnergyLoss_1", "Energy Difference of Most Energetic Subjet compared to original Jet Detector Level",101, -0.05,1.05);
536  fhSubJetEnergyLoss_2= new TH1F("fhSubJetEnergyLoss_2", "Energy Difference of Most Energetic Subjet compared to original Jet Particle Level",101, -0.05,1.05);
538  fhSubJetEnergyLoss2_1= new TH1F("fhSubJetEnergyLoss2_1", "Pt Difference of Two Most Energetic Subjets compared to original Jet Detector Level",101, -0.05,1.05);
540  fhSubJetEnergyLoss2_2= new TH1F("fhSubJetEnergyLoss2_2", "Pt Difference of Two Most Energetic Subjets compared to original Jet Particle Level",101, -0.05,1.05);
542  fhSubJetiness1_1= new TH1F("fhSubjetiness1_1", "#Tau 1 value Detector Level",101, -0.05,1.05);
544  fhSubJetiness1_2= new TH1F("fhSubjetiness1_2", "#Tau 1 value Particle Level",101, -0.05,1.05);
546  fhSubJetiness1JetPt_1= new TH2F("fhSubjetiness1JetPt_1", "#Tau 1 value vs Jet Pt Detector Level",1500, -0.5, 149.5, 101, -0.05,1.05);
548  fhSubJetiness1JetPt_2= new TH2F("fhSubjetiness1JetPt_2", "#Tau 1 value vs Jet Pt Particle Level",1500, -0.5, 149.5, 101, -0.05,1.05);
550  fhSubJetiness2_1= new TH1F("fhSubjetiness2_1", "#Tau 2 value Detector Level",101, -0.05,1.05);
552  fhSubJetiness2_2= new TH1F("fhSubjetiness2_2", "#Tau 2 value Particle Level",101, -0.05,1.05);
554  fhSubJetiness2JetPt_1= new TH2F("fhSubjetiness2JetPt_1", "#Tau 2 value vs Jet Pt Detector Level",1500, -0.5, 149.5, 101, -0.05,1.05);
556  fhSubJetiness2JetPt_2= new TH2F("fhSubjetiness2JetPt_2", "#Tau 2 value vs Jet Pt Particle Level",1500, -0.5, 149.5, 101, -0.05,1.05);
558  fhJetCounter_1= new TH1F("fhJetCounter_1", "Jet Counter Detector Level", 100, -0.5, 99.5);
559  fOutput->Add(fhJetCounter_1);
560  fhJetCounter_2= new TH1F("fhJetCounter_2", "Jet Counter Particle Level", 100, -0.5, 99.5);
561  fOutput->Add(fhJetCounter_2);
562  fhSubJetCounter_1 = new TH1F("fhSubJetCounter_1", "SubJet Counter Detector Level",50, -0.5,49.5);
564  fhSubJetCounter_2 = new TH1F("fhSubJetCounter_2", "SubJet Counter Particle Level",50, -0.5,49.5);
566  fh2to1SubJetinessRatio_1= new TH1F("fh2to1SubJetinessRatio_1", "Ratio of #tau 1 to #tau 2 Detector Level",200, -0.5,1.5);
568  fh2to1SubJetinessRatio_2= new TH1F("fh2to1SubJetinessRatio_2", "Ratio of #tau 1 to #tau 2 Particle Level",200, -0.5,1.5);
570  fh2to1SubJetinessRatioJetPt_1= new TH2F("fh2to1SubJetinessRatioJetPt_1", "Ratio of #tau 1 to #tau 2 vs Jet Pt Detector Level", 1500, -0.5, 149.5, 200, -0.5,1.5);
572  fh2to1SubJetinessRatioJetPt_2= new TH2F("fh2to1SubJetinessRatioJetPt_2", "Ratio of #tau 1 to #tau 2 vs Jet Pt Particle Level", 1500, -0.5, 149.5, 200, -0.5,1.5);
574  fhJetPtJetEta_1 = new TH2F("fhJetPtJetEta_1", "Jet Pt vs Jet Eta Detector Level", 1500,-0.5,150,Eta_Bins, Eta_Low, Eta_Up);
575  fOutput->Add(fhJetPtJetEta_1);
576  fhJetPtJetEta_2 = new TH2F("fhJetPtJetEta_2", "Jet Pt vs Jet Eta Particle Level", 1500,-0.5,150,Eta_Bins, Eta_Low, Eta_Up);
577  fOutput->Add(fhJetPtJetEta_2);
578  fhSubJetiness2Distance_1 = new TH2F("fhSubJetiness2Distance_1", "#tau 2 as a function of distance between subject axis Detector Level",100,0,10,101,-0.05,1.05);
580  fhSubJetiness2Distance_2 = new TH2F("fhSubJetiness2Distance_2", "#tau 2 as a function of distance between subject axis Particle Level",100,0,10,101,-0.05,1.05);
582  fh2PtRatio= new TH2F("fhPtRatio", "MC pt for detector level vs particle level jets",1500,-0.5,149.5,1500,-0.5,149.5);
583  fOutput->Add(fh2PtRatio);
584  fhEventCounter_1= new TH1F("fhEventCounter_1", "Event Counter Detector Level", 15,0.5,15.5);
586  fhEventCounter_2= new TH1F("fhEventCounter_2", "Event Counter Particle Level", 15,0.5,15.5);
588  }
590  fhEventCounter= new TH1F("fhEventCounter", "Event Counter", 15,0.5,15.5);
591  fOutput->Add(fhEventCounter);
592  }
594  TH1::AddDirectory(oldStatus);
595  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
596 }
597 
598 //________________________________________________________________________
600 {
601  // Run analysis code here, if needed. It will be executed before FillHistograms().
602 
603 
604  return kTRUE;
605 }
606 
607 //________________________________________________________________________
609 {
610 
611  //fhEventCounter Key:
612  // 1: Number of events with a Jet Container
613  // 2: Number of Jets without a Jet Container
614  // 3:
615  // 4: Number of Jets found in all events
616  // 5: Number of Jets that were reclustered in all events
617  // 6: Number of SubJets found in all events
618  // 7: Number of events were primary vertext was found
619  // 8: Number of Jets with more than one SubJet
620  // 9:Number of Jets with more than two SubJets
621  // 12:Number of SubJetinessEvents in kData
622 
623 
624  // const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
625  // Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
626  // if(vert) fhEventCounter->Fill(7);
627  if (fCentSelectOn){
628  if ((fCent>fCentMax) || (fCent<fCentMin)) return 0;
629  }
630 
632  AliEmcalJet *Jet1 = NULL; //Subtracted hybrid Jet
633  AliEmcalJet *Jet2 = NULL; //Unsubtracted Hybrid Jet
634  AliEmcalJet *Jet3 = NULL; //Detector Level Pythia Jet
635  AliEmcalJet *Jet4 = NULL; //Particle Level Pyhtia Jet
636  AliJetContainer *JetCont1= GetJetContainer(0); //Jet Container for Subtracted Hybrid Jets
637  AliJetContainer *JetCont2= GetJetContainer(1); //Jet Container for Unsubtracted Hybrid Jets
638  AliJetContainer *JetCont3= GetJetContainer(2); //Jet Container for Detector Level Pyhtia Jets
639  AliJetContainer *JetCont4= GetJetContainer(3); //Jet Container for Particle Level Pythia Jets
640  AliEmcalJetFinder *Reclusterer1; //Object containg Subjets from Subtracted Hybrid Jets
641  AliEmcalJetFinder *Reclusterer2; //Object containg Subjets from Unsubtracted Hybrid Jets
642  AliEmcalJetFinder *Reclusterer3; //Object containg Subjets from Detector Level Pythia Jets
643  AliEmcalJetFinder *Reclusterer4; //Object containg Subjets from Particle Level Pythia Jets
644  Double_t Two_Hardest_SubJet_Distance1=0;
645  Double_t Two_Hardest_SubJet_Distance4=0;
646  Bool_t JetsMatched=kFALSE;
647  Int_t PartonFlag1=0;
648  Int_t PartonFlag4=0;
649  Bool_t EventCounter=kFALSE;
650  Int_t JetNumber=-1;
651  fhEventCounter->Fill(1);
652  if(JetCont1) {
653  fhEventCounter->Fill(2);
654  JetCont1->ResetCurrentID();
655  while((Jet1=JetCont1->GetNextAcceptJet())) {
656  if( (!Jet1) || (Jet1->Pt()<fPtThreshold)) {
657  continue;
658  }
659  else {
660  if (!EventCounter){
661  fhEventCounter->Fill(3);
662  EventCounter=kTRUE;
663  }
664  Reclusterer1 = Recluster(Jet1, 0, fSubJetRadius, fSubJetMinPt, fSubJetAlgorithm, "SubJetFinder_1");
665  if (Reclusterer1->GetNumberOfJets()>=2){
666  Two_Hardest_SubJet_Distance1=TMath::Sqrt(TMath::Power(Reclusterer1->GetJet(SubJetOrdering(Jet1, Reclusterer1, 1, 0, kTRUE))->Phi()-Reclusterer1->GetJet(SubJetOrdering(Jet1,Reclusterer1, 2, 0, kTRUE))->Phi(),2)+TMath::Power(Reclusterer1->GetJet(SubJetOrdering(Jet1, Reclusterer1, 1, 0, kTRUE))->Eta()-Reclusterer1->GetJet(SubJetOrdering(Jet1, Reclusterer1, 2, 0, kTRUE))->Eta(),2));
667  if (Two_Hardest_SubJet_Distance1>(2*fJetRadius)){
668  Two_Hardest_SubJet_Distance1=TMath::Sqrt(TMath::Power((2*TMath::Pi())-TMath::Abs(Reclusterer1->GetJet(SubJetOrdering(Jet1, Reclusterer1, 1, 0, kTRUE))->Phi()-Reclusterer1->GetJet(SubJetOrdering(Jet1,Reclusterer1, 2, 0, kTRUE))->Phi()),2)+TMath::Power(Reclusterer1->GetJet(SubJetOrdering(Jet1, Reclusterer1, 1, 0, kTRUE))->Eta()-Reclusterer1->GetJet(SubJetOrdering(Jet1, Reclusterer1, 2,0, kTRUE))->Eta(),2));
669  }
670  }
671  else Two_Hardest_SubJet_Distance1=-2;
672  if(fJetShapeSub==kConstSub){
673  JetNumber=-1;
674  for(Int_t i = 0; i<JetCont2->GetNJets(); i++) {
675  Jet2 = JetCont2->GetJet(i);
676  if(Jet2->GetLabel()==Jet1->GetLabel()) {
677  JetNumber=i;
678  }
679  }
680  if(JetNumber==-1) continue;
681  Jet2=JetCont2->GetJet(JetNumber);
682  if (JetCont2->AliJetContainer::GetFractionSharedPt(Jet2)<fSharedFractionPtMin) continue;
683  Jet3=Jet2->ClosestJet();
684  }
685  if(!(fJetShapeSub==kConstSub)){
686  if (!(JetCont1->AliJetContainer::GetFractionSharedPt(Jet1)<fSharedFractionPtMin)) continue;
687  Jet3 = Jet1->ClosestJet();
688  }
689  if (!Jet3) continue;
690  Jet4=Jet3->ClosestJet();
691  if(!Jet4) continue;
692  JetsMatched=kTRUE;
693  Reclusterer4 = Recluster(Jet4, 3, fSubJetRadius, 0, fSubJetAlgorithm, "SubJetFinder_4");
694  if (Reclusterer4->GetNumberOfJets()>=2){
695  Two_Hardest_SubJet_Distance4=TMath::Sqrt(TMath::Power(Reclusterer4->GetJet(SubJetOrdering(Jet4, Reclusterer4, 1, 0, kTRUE))->Phi()-Reclusterer4->GetJet(SubJetOrdering(Jet4,Reclusterer4, 2, 0, kTRUE))->Phi(),2)+TMath::Power(Reclusterer4->GetJet(SubJetOrdering(Jet4, Reclusterer4, 1, 0, kTRUE))->Eta()-Reclusterer4->GetJet(SubJetOrdering(Jet4, Reclusterer4, 2, 0, kTRUE))->Eta(),2));
696  if (Two_Hardest_SubJet_Distance4>(2*fJetRadius)){
697  Two_Hardest_SubJet_Distance4=TMath::Sqrt(TMath::Power((2*TMath::Pi())-TMath::Abs(Reclusterer4->GetJet(SubJetOrdering(Jet4, Reclusterer4, 1, 0, kTRUE))->Phi()-Reclusterer4->GetJet(SubJetOrdering(Jet4,Reclusterer4, 2, 0, kTRUE))->Phi()),2)+TMath::Power(Reclusterer4->GetJet(SubJetOrdering(Jet4, Reclusterer4, 1, 0, kTRUE))->Eta()-Reclusterer4->GetJet(SubJetOrdering(Jet4, Reclusterer4, 2,0, kTRUE))->Eta(),2));
698  }
699  }
700  else Two_Hardest_SubJet_Distance4=-2;
702  fShapesVar[0]=Jet1->Pt();
703  fShapesVar[2]=SubJetFraction(Jet1, Reclusterer1, 1, 0, kTRUE, kFALSE);
704  fShapesVar[4]=SubJetFraction(Jet1, Reclusterer1, 2, 0, kTRUE, kFALSE);
705  fShapesVar[6]=NSubJettiness(Jet1, 0, fJetRadius, Reclusterer1, 1, 0, 1);
706  fShapesVar[8]=NSubJettiness(Jet1, 0, fJetRadius, Reclusterer1, 2, 0, 1);
707  fShapesVar[10]=Two_Hardest_SubJet_Distance1;
708  fShapesVar[12]=PartonFlag1;
709  if (JetsMatched){ //even needed? Not now but might be if you want to fill trees when jets aren't matched too
710  fShapesVar[1]=Jet4->Pt();
711  fShapesVar[3]=SubJetFraction(Jet4, Reclusterer4, 1, 0, kTRUE, kFALSE);
712  fShapesVar[5]=SubJetFraction(Jet4, Reclusterer4, 2, 0, kTRUE, kFALSE);
713  fShapesVar[7]=NSubJettiness(Jet4, 3, fJetRadius, Reclusterer4, 1, 0, 1);
714  fShapesVar[9]=NSubJettiness(Jet4, 3, fJetRadius, Reclusterer4, 2, 0, 1);
715  fShapesVar[11]=Two_Hardest_SubJet_Distance4;
716  fShapesVar[13]=PartonFlag4;
717  }
718  else{
719  fShapesVar[1]=-990;
720  fShapesVar[3]=-2;
721  fShapesVar[5]=-2;
722  fShapesVar[7]=-2;
723  fShapesVar[9]=-2;
724  fShapesVar[11]=-2;
725  fShapesVar[13]=-2;
726  }
727  fTreeResponseMatrixAxis->Fill();
728  JetsMatched=kFALSE;
729  }
730  }
731  }
732  }
733 
734 
735 
736 
737 
738 
739 
740 
741 
742 
743 
744 
745 
746 
747 
748 
749 
750 
751 
752 
753 
754 
755 
756 
757 
758 
759 
760 
761 
762 
763 
764 
765 
767  AliEmcalJet *Jet1 = NULL; //Detector Level Jet
768  AliEmcalJet *Jet2 = NULL; //Particle Level Jet
769  AliJetContainer *JetCont1= GetJetContainer(0); //Jet Container for Detector Level Pythia
770  AliJetContainer *JetCont2= GetJetContainer(1); //Jet Container for Particle Level Pythia
771  //const AliEmcalPythiaInfo *PartonInfo1 = 0x0;
772  //const AliEmcalPythiaInfo *PartonInfo2 = 0x0;
773  AliEmcalJetFinder *Reclusterer1; //Object containg Subjets Detector Level
774  AliEmcalJetFinder *Reclusterer2; //Object containg Subjets Particle Level
775  Int_t JetCounter1=0; //Counts number of jets in event
776  Int_t JetCounter2=0; //Counts number of jets in event
777  Double_t JetPhi1=0;
778  Double_t JetPhi2=0;
779  Double_t Two_Hardest_SubJet_Distance1=0;
780  Double_t Two_Hardest_SubJet_Distance2=0;
781  Bool_t JetsMatched=kFALSE;
782  Double_t SubJettiness1_1 =0; //1SubJettiness Detector Level
783  Double_t SubJettiness2_1 =0; //2SubJettiness Detector Level
784  Double_t SubJettiness1_2 =0; //1SubJettiness Particle Level
785  Double_t SubJettiness2_2 =0; //2SubJettiness Particle Level
786  Double_t Pythia_Event_Weight=1;
787  Double_t PartonFlag1=0;
788  Double_t PartonFlag2=0;
789  Bool_t EventCounter=kFALSE;
790  fhEventCounter_1->Fill(1);
791  if(JetCont1) {
792  fhEventCounter_1->Fill(2); //Number of events with a jet container
793  JetCont1->ResetCurrentID();
794  while((Jet1=JetCont1->GetNextAcceptJet())) {
795  if( (!Jet1) || ((Jet1->Pt())<fPtThreshold)) {
796  // fhEventCounter_1->Fill(3); //events where the jet had a problem
797  continue;
798  }
799  else {
800  if (!EventCounter){
801  fhEventCounter_1->Fill(3);
802  EventCounter=kTRUE;
803  }
804  if((Jet1->GetNumberOfTracks())==0){
805  fhEventCounter_1->Fill(10); //zero track jets
806  }
807  if((Jet1->GetNumberOfTracks())==1){
808  fhEventCounter_1->Fill(11); //one track jets
809  }
810  fhEventCounter_1->Fill(4); //Number of Jets found in all events
811  JetCounter1++;
812  fhJetPt_1->Fill(Jet1->Pt());
813  JetPhi1=Jet1->Phi();
814  if(JetPhi1 < -1*TMath::Pi()) JetPhi1 += (2*TMath::Pi());
815  else if (JetPhi1 > TMath::Pi()) JetPhi1 -= (2*TMath::Pi());
816  fhJetPhi_1->Fill(JetPhi1);
817  fhJetPtJetEta_1->Fill(Jet1->Pt(),Jet1->Eta());
818  fhJetEta_1->Fill(Jet1->Eta());
819  fhJetMass_1->Fill(Jet1->M());
820  fhJetRadius_1->Fill(TMath::Sqrt((Jet1->Area()/TMath::Pi()))); //Radius of Jets per event
822  fhJetAngularity_1->Fill(Angularity(Jet1,0)); //Angularity Jet Shape
823  fhJetAngularityJetPt_1->Fill(Jet1->Pt(),Angularity(Jet1,0));
824  fhJetPTD_1->Fill(PTD(Jet1,0)); //PTD Jet Shape
825  fhJetPTDJetPt_1->Fill(Jet1->Pt(),PTD(Jet1,0));
826  Reclusterer1 = Recluster(Jet1, 0, fSubJetRadius, fSubJetMinPt, fSubJetAlgorithm, "SubJetFinder_1");
827  fhSubJetCounter_1->Fill(Reclusterer1->GetNumberOfJets());
828  for (Int_t i= 0; i<Reclusterer1->GetNumberOfJets(); i++){
829  fhEventCounter_1->Fill(6); //Number of overall subjets in all events
830  fhSubJetPt_1->Fill(Reclusterer1->GetJet(i)->Pt());
831  fhSubJetMass_1->Fill(Reclusterer1->GetJet(i)->M());
832  fhNumberOfSubJetTracks_1->Fill(Reclusterer1->GetJet(i)->GetNumberOfTracks());
833  fhSubJetRadius_1->Fill(TMath::Sqrt((Reclusterer1->GetJet(i)->Area()/TMath::Pi()))); //Radius of SubJets per event
834  }
835  fhSubJetPtFrac_1->Fill(SubJetFraction(Jet1, Reclusterer1, 1, 0, kTRUE, kFALSE));
836  fhSubJetPtLoss_1->Fill(SubJetFraction(Jet1, Reclusterer1, 1, 0, kTRUE, kTRUE));
837  fhSubJetEnergyFrac_1->Fill(SubJetFraction(Jet1, Reclusterer1, 1, 1, kTRUE, kFALSE));
838  fhSubJetEnergyLoss_1->Fill(SubJetFraction(Jet1, Reclusterer1, 1, 1, kTRUE, kTRUE));
839  fhSubJetPtFrac2_1->Fill(SubJetFraction(Jet1, Reclusterer1, 2, 0, kTRUE, kFALSE));
840  fhSubJetPtLoss2_1->Fill(SubJetFraction(Jet1, Reclusterer1, 2, 0, kTRUE, kTRUE));
841  fhSubJetEnergyFrac2_1->Fill(SubJetFraction(Jet1, Reclusterer1, 2, 1, kTRUE, kFALSE));
842  fhSubJetEnergyLoss2_1->Fill(SubJetFraction(Jet1, Reclusterer1, 2, 1, kTRUE, kTRUE));
843  SubJettiness1_1=NSubJettiness(Jet1, 0, fJetRadius, Reclusterer1, 1, 0, 1);
844  SubJettiness2_1=NSubJettiness(Jet1, 0, fJetRadius, Reclusterer1, 2, 0, 1);
845  fhSubJetiness1_1->Fill(SubJettiness1_1);
846  fhSubJetiness1JetPt_1->Fill(Jet1->Pt(),SubJettiness1_1);
847  fhSubJetiness2_1->Fill(SubJettiness2_1);
848  fhSubJetiness2JetPt_1->Fill(Jet1->Pt(),SubJettiness2_1);
849  if(SubJettiness1_1!=-2 && SubJettiness2_1!=-2){ //have to be careful on ratios
850  fh2to1SubJetinessRatio_1->Fill(SubJettiness2_1/SubJettiness1_1);
851  fh2to1SubJetinessRatioJetPt_1->Fill(Jet1->Pt(),SubJettiness2_1/SubJettiness1_1);
852  }
853  if (Reclusterer1->GetNumberOfJets()>=2){
854  Two_Hardest_SubJet_Distance1=TMath::Sqrt(TMath::Power(Reclusterer1->GetJet(SubJetOrdering(Jet1, Reclusterer1, 1, 0, kTRUE))->Phi()-Reclusterer1->GetJet(SubJetOrdering(Jet1,Reclusterer1, 2, 0, kTRUE))->Phi(),2)+TMath::Power(Reclusterer1->GetJet(SubJetOrdering(Jet1, Reclusterer1, 1, 0, kTRUE))->Eta()-Reclusterer1->GetJet(SubJetOrdering(Jet1, Reclusterer1,2, 0, kTRUE))->Eta(),2));
855  if (Two_Hardest_SubJet_Distance1>(2*fJetRadius)){//this can happen when the two subjet axis fall either side of the Phi boundary (i.e one close to zero the other close to 2*Pi IS THIS CORRECT?????????????????????????????????????????????????????????
856  Two_Hardest_SubJet_Distance1=TMath::Sqrt(TMath::Power((2*TMath::Pi())-TMath::Abs(Reclusterer1->GetJet(SubJetOrdering(Jet1, Reclusterer1, 1, 0, kTRUE))->Phi()-Reclusterer1->GetJet(SubJetOrdering(Jet1,Reclusterer1, 2, 0, kTRUE))->Phi()),2)+TMath::Power(Reclusterer1->GetJet(SubJetOrdering(Jet1, Reclusterer1, 1, 0, kTRUE))->Eta()-Reclusterer1->GetJet(SubJetOrdering(Jet1, Reclusterer1, 2,0, kTRUE))->Eta(),2));
857  }
858  }
859  else Two_Hardest_SubJet_Distance1=-2;
860  /*
861  PartonInfo1 = (AliPythiaInfo*) JetCont1->GetPythiaInfo();
862  if(!PartonInfo1) {
863  return 0;
864  }
865  Pythia_Event_Weight=PartonInfo1->GetPythiaEventWeight();
866  if(TMath::Sqrt(TMath::Power(RelativePhi(Jet1->Phi(),PartonInfo1->GetPartonPhi6()),2)+TMath::Power(Jet1->Eta()-PartonInfo1->GetPartonEta6(),2))<fRMatching){
867  PartonFlag1=PartonInfo1->GetPartonFlag6(); //6 and 7 refer to the back to back scattered partons. The flag returns a number corresponding to the parton flavour
868  }
869  else if(TMath::Sqrt(TMath::Power(RelativePhi(Jet1->Phi(),PartonInfo1->GetPartonPhi7()),2)+TMath::Power(Jet1->Eta()-PartonInfo1->GetPartonEta7(),2))<fRMatching){
870  PartonFlag1=PartonInfo1->GetPartonFlag7();
871  }
872  else PartonFlag1=0;
873  */
874  if((Jet2 = Jet1->ClosestJet())){
875  JetsMatched=kTRUE;
876  if((Jet2->GetNumberOfTracks())==0){
877  fhEventCounter_2->Fill(10); //zero track jets
878  }
879  if((Jet2->GetNumberOfTracks())==1){
880  fhEventCounter_2->Fill(11); //one track jets
881  }
882  fhEventCounter_2->Fill(4); //Number of Jets found in all events
883  JetCounter2++;
884  fhJetPt_2->Fill(Jet2->Pt());
885  JetPhi2=Jet2->Phi();
886  if(JetPhi2 < -1*TMath::Pi()) JetPhi2 += (2*TMath::Pi());
887  else if (JetPhi2 > TMath::Pi()) JetPhi2 -= (2*TMath::Pi());
888  fhJetPhi_2->Fill(JetPhi2);
889  fhJetPtJetEta_2->Fill(Jet2->Pt(),Jet2->Eta());
890  fhJetEta_2->Fill(Jet2->Eta());
891  fhJetMass_2->Fill(Jet2->M());
892  fhJetRadius_2->Fill(TMath::Sqrt((Jet2->Area()/TMath::Pi()))); //Radius of Jets per event
894  fhJetAngularity_2->Fill(Angularity(Jet2,1)); //Angularity Jet Shape
895  fhJetAngularityJetPt_2->Fill(Jet2->Pt(),Angularity(Jet2,1));
896  fhJetPTD_2->Fill(PTD(Jet2,1)); //PTD Jet Shape
897  fhJetPTDJetPt_2->Fill(Jet2->Pt(),PTD(Jet2,1));
898  Reclusterer2 = Recluster(Jet2, 1, fSubJetRadius, 0, fSubJetAlgorithm, "SubJetFinder_2");
899  fhSubJetCounter_2->Fill(Reclusterer2->GetNumberOfJets());
900  for (Int_t i= 0; i<Reclusterer2->GetNumberOfJets(); i++){
901  fhEventCounter_2->Fill(6); //Number of overall subjets in all events
902  fhSubJetPt_2->Fill(Reclusterer2->GetJet(i)->Pt());
903  fhSubJetMass_2->Fill(Reclusterer2->GetJet(i)->M());
904  fhNumberOfSubJetTracks_2->Fill(Reclusterer2->GetJet(i)->GetNumberOfTracks());
905  fhSubJetRadius_2->Fill(TMath::Sqrt((Reclusterer2->GetJet(i)->Area()/TMath::Pi()))); //Radius of SubJets per event
906  }
907  fhSubJetPtFrac_2->Fill(SubJetFraction(Jet2, Reclusterer2, 1, 0, kTRUE, kFALSE));
908  fhSubJetPtLoss_2->Fill(SubJetFraction(Jet2, Reclusterer2, 1, 0, kTRUE, kTRUE));
909  fhSubJetEnergyFrac_2->Fill(SubJetFraction(Jet2, Reclusterer2, 1, 1, kTRUE, kFALSE));
910  fhSubJetEnergyLoss_2->Fill(SubJetFraction(Jet2, Reclusterer2, 1, 1, kTRUE, kTRUE));
911  fhSubJetPtFrac2_2->Fill(SubJetFraction(Jet2, Reclusterer2, 2, 0, kTRUE, kFALSE));
912  fhSubJetPtLoss2_2->Fill(SubJetFraction(Jet2, Reclusterer2, 2, 0, kTRUE, kTRUE));
913  fhSubJetEnergyFrac2_2->Fill(SubJetFraction(Jet2, Reclusterer2, 2, 1, kTRUE, kFALSE));
914  fhSubJetEnergyLoss2_2->Fill(SubJetFraction(Jet2, Reclusterer2, 2, 1, kTRUE, kTRUE));
915  SubJettiness1_2=NSubJettiness(Jet2, 1, fJetRadius, Reclusterer2, 1, 0, 1);
916  SubJettiness2_2=NSubJettiness(Jet2, 1, fJetRadius, Reclusterer2, 2, 0, 1);
917  fhSubJetiness1_2->Fill(SubJettiness1_2);
918  fhSubJetiness1JetPt_2->Fill(Jet2->Pt(),SubJettiness1_2);
919  fhSubJetiness2_2->Fill(SubJettiness2_2);
920  fhSubJetiness2JetPt_2->Fill(Jet2->Pt(),SubJettiness2_2);
921  if(SubJettiness1_2!=-2 && SubJettiness2_2!=-2){ //have to be careful on ratios
922  fh2to1SubJetinessRatio_2->Fill(SubJettiness2_2/SubJettiness1_2);
923  fh2to1SubJetinessRatioJetPt_2->Fill(Jet2->Pt(),SubJettiness2_2/SubJettiness1_2);
924  }
925  if (Reclusterer2->GetNumberOfJets()>=2){
926  Two_Hardest_SubJet_Distance2=TMath::Sqrt(TMath::Power(Reclusterer2->GetJet(SubJetOrdering(Jet2, Reclusterer2, 1, 0, kTRUE))->Phi()-Reclusterer2->GetJet(SubJetOrdering(Jet2, Reclusterer2, 2, 0, kTRUE))->Phi(),2)+TMath::Power(Reclusterer2->GetJet(SubJetOrdering(Jet2, Reclusterer2, 1, 0, kTRUE))->Eta()-Reclusterer2->GetJet(SubJetOrdering(Jet2, Reclusterer2,2, 0, kTRUE))->Eta(),2));
927  if (Two_Hardest_SubJet_Distance2>(2*fJetRadius)){//this can happen when the two subjet axis fall either side of the Phi boundary (i.e one close to zero the other closeto 2*Pi IS THIS CORRECT?????????????????????????????????????????????????????????
928  Two_Hardest_SubJet_Distance2=TMath::Sqrt(TMath::Power((2*TMath::Pi())-TMath::Abs(Reclusterer2->GetJet(SubJetOrdering(Jet2, Reclusterer2, 1, 0, kTRUE))->Phi()-Reclusterer2->GetJet(SubJetOrdering(Jet2,Reclusterer2, 2, 0, kTRUE))->Phi()),2)+TMath::Power(Reclusterer2->GetJet(SubJetOrdering(Jet2, Reclusterer2, 1, 0, kTRUE))->Eta()-Reclusterer2->GetJet(SubJetOrdering(Jet2, Reclusterer2, 2,0, kTRUE))->Eta(),2));
929  }
930  }
931  else Two_Hardest_SubJet_Distance2=-2;
932 /*
933  PartonInfo2 = (AliPythiaInfo*) JetCont2->GetPythiaInfo();
934  if(!PartonInfo2) return 0;
935  if(TMath::Sqrt(TMath::Power(RelativePhi(Jet2->Phi(),PartonInfo2->GetPartonPhi6()),2)+TMath::Power(Jet2->Eta()-PartonInfo2->GetPartonEta6(),2))<fRMatching){
936  PartonFlag2=PartonInfo2->GetPartonFlag6();
937  }
938  else if(TMath::Sqrt(TMath::Power(RelativePhi(Jet2->Phi(),PartonInfo2->GetPartonPhi7()),2)+TMath::Power(Jet2->Eta()-PartonInfo2->GetPartonEta7(),2))<fRMatching){
939  PartonFlag2=PartonInfo2->GetPartonFlag7();
940  }
941  else PartonFlag2=0;
942 */
943  fh2PtRatio->Fill(Jet1->Pt(),Jet2->Pt(),Pythia_Event_Weight);
944  }
945  else {
946  fhEventCounter_2->Fill(1);
947  //continue;
948  }
950 
951  fhSubJetiness2Distance_1->Fill(Two_Hardest_SubJet_Distance1,SubJettiness2_1);
952  fShapesVar[0]=Jet1->Pt();
953  fShapesVar[2]=SubJetFraction(Jet1, Reclusterer1, 1, 0, kTRUE, kFALSE);
954  fShapesVar[4]=SubJetFraction(Jet1, Reclusterer1, 2, 0, kTRUE, kFALSE);
955  fShapesVar[6]=SubJettiness1_1;
956  fShapesVar[8]=SubJettiness2_1;
957  fShapesVar[10]=Two_Hardest_SubJet_Distance1;
958  fShapesVar[12]=PartonFlag1;
959  if (JetsMatched){
960  fhSubJetiness2Distance_2->Fill(Two_Hardest_SubJet_Distance2,SubJettiness2_2);
961  fShapesVar[1]=Jet2->Pt();
962  fShapesVar[3]=SubJetFraction(Jet2, Reclusterer2, 1, 0, kTRUE, kFALSE);
963  fShapesVar[5]=SubJetFraction(Jet2, Reclusterer2, 2, 0, kTRUE, kFALSE);
964  fShapesVar[7]=SubJettiness1_2;
965  fShapesVar[9]=SubJettiness2_2;
966  fShapesVar[11]=Two_Hardest_SubJet_Distance2;
967  fShapesVar[13]=PartonFlag2;
968  }
969  else{
970  fShapesVar[1]=-990;
971  fShapesVar[3]=-2;
972  fShapesVar[5]=-2;
973  fShapesVar[7]=-2;
974  fShapesVar[9]=-2;
975  fShapesVar[11]=-2;
976  fShapesVar[13]=-2;
977  }
978  fTreeResponseMatrixAxis->Fill();
979  JetsMatched=kFALSE;
980  }
981  }
982  fhJetCounter_1->Fill(JetCounter1); //Number of Jets in Each Event Particle Level
983  fhJetCounter_2->Fill(JetCounter2); //Number of Jets in Each Event Detector Level
984  }
985  //else {fhEventCounter_1->Fill(3);} //Events with no jet container
986  }
987 
989 
990  AliEmcalJet *Jet1 = NULL; //Original Jet in the event
991  AliJetContainer *JetCont= GetJetContainer(0); //Jet Container for event
992  Int_t JetCounter=0; //Counts number of jets in event
993  Double_t JetPhi=0;
994  Double_t Two_Hardest_SubJet_Distance=0;
995  Double_t SubJettiness1 =0;
996  Double_t SubJettiness2 =0;
997  Bool_t EventCounter=kFALSE;
998  Double_t JetPt_ForThreshold=0;
999  fhEventCounter->Fill(1);
1000  if(JetCont) {
1001  fhEventCounter->Fill(2); //Number of events with a jet container
1002  JetCont->ResetCurrentID();
1003  while((Jet1=JetCont->GetNextAcceptJet())) {
1004  if(!Jet1) continue;
1005  if(fJetShapeSub==kNoSub) JetPt_ForThreshold = Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
1006  else JetPt_ForThreshold = Jet1->Pt();
1007  if(JetPt_ForThreshold<fPtThreshold) {
1008  //fhEventCounter->Fill(3); //events where the jet had a problem
1009  continue;
1010  }
1011  else {
1012  if (!EventCounter){
1013  fhEventCounter->Fill(3);
1014  EventCounter=kTRUE;
1015  }
1016  if((Jet1->GetNumberOfTracks())==0){
1017  fhEventCounter->Fill(10); //zero track jets
1018  }
1019  if((Jet1->GetNumberOfTracks())==1){
1020  fhEventCounter->Fill(11); //one track jets
1021  }
1022  fhEventCounter->Fill(4); //Number of Jets found in all events
1023  JetCounter++;
1024  fhJetPt->Fill(Jet1->Pt());
1025  JetPhi=Jet1->Phi();
1026  if(JetPhi < -1*TMath::Pi()) JetPhi += (2*TMath::Pi());
1027  else if (JetPhi > TMath::Pi()) JetPhi -= (2*TMath::Pi());
1028  fhJetPhi->Fill(JetPhi);
1029  fhJetPtJetEta->Fill(Jet1->Pt(),Jet1->Eta());
1030  fhJetEta->Fill(Jet1->Eta());
1031  fhJetMass->Fill(Jet1->M());
1032  fhJetRadius->Fill(TMath::Sqrt((Jet1->Area()/TMath::Pi()))); //Radius of Jets per event
1033  fhNumberOfJetTracks->Fill(Jet1->GetNumberOfTracks());
1034  fhJetAngularity->Fill(Angularity(Jet1,0)); //Angularity Jet Shape
1035  fhJetAngularityJetPt->Fill(Jet1->Pt(),Angularity(Jet1,0));
1036  fhJetPTD->Fill(PTD(Jet1,0)); //PTD Jet Shape
1037  fhJetPTDJetPt->Fill(Jet1->Pt(),PTD(Jet1,0));
1038  AliEmcalJetFinder *Reclusterer = Recluster(Jet1, 0, fSubJetRadius, fSubJetMinPt, fSubJetAlgorithm, "SubJetFinder");
1039  fhSubJetCounter->Fill(Reclusterer->GetNumberOfJets());
1040  for (Int_t i= 0; i<Reclusterer->GetNumberOfJets(); i++){
1041  fhEventCounter->Fill(6); //Number of overall subjets in all events
1042  fhSubJetPt->Fill(Reclusterer->GetJet(i)->Pt());
1043  fhSubJetMass->Fill(Reclusterer->GetJet(i)->M());
1044  fhNumberOfSubJetTracks->Fill(Reclusterer->GetJet(i)->GetNumberOfTracks());
1045  fhSubJetRadius->Fill(TMath::Sqrt((Reclusterer->GetJet(i)->Area()/TMath::Pi()))); //Radius of SubJets per event
1046  }
1047  fhSubJetPtFrac->Fill(SubJetFraction(Jet1, Reclusterer, 1, 0, kTRUE, kFALSE));
1048  fhSubJetPtLoss->Fill(SubJetFraction(Jet1, Reclusterer, 1, 0, kTRUE, kTRUE));
1049  fhSubJetEnergyFrac->Fill(SubJetFraction(Jet1, Reclusterer, 1, 1, kTRUE, kFALSE));
1050  fhSubJetEnergyLoss->Fill(SubJetFraction(Jet1, Reclusterer, 1, 1, kTRUE, kTRUE));
1051  fhSubJetPtFrac2->Fill(SubJetFraction(Jet1, Reclusterer, 2, 0, kTRUE, kFALSE));
1052  fhSubJetPtLoss2->Fill(SubJetFraction(Jet1, Reclusterer, 2, 0, kTRUE, kTRUE));
1053  fhSubJetEnergyFrac2->Fill(SubJetFraction(Jet1, Reclusterer, 2, 1, kTRUE, kFALSE));
1054  fhSubJetEnergyLoss2->Fill(SubJetFraction(Jet1, Reclusterer, 2, 1, kTRUE, kTRUE));
1055  SubJettiness1=NSubJettiness(Jet1, 0, fJetRadius, Reclusterer, 1, 0, 1);
1056  SubJettiness2=NSubJettiness(Jet1, 0, fJetRadius, Reclusterer, 2, 0, 1);
1057  fhSubJetiness1->Fill(SubJettiness1);
1058  fhSubJetiness1JetPt->Fill(Jet1->Pt(),SubJettiness1);
1059  fhSubJetiness2->Fill(SubJettiness2);
1060  fhSubJetiness2JetPt->Fill(Jet1->Pt(),SubJettiness2);
1061  if(SubJettiness1!=-2 && SubJettiness2!=-2){ //have to be careful on ratios
1062  fh2to1SubJetinessRatio->Fill(SubJettiness2/SubJettiness1);
1063  fh2to1SubJetinessRatioJetPt->Fill(Jet1->Pt(),SubJettiness2/SubJettiness1);
1064  }
1065  if (Reclusterer->GetNumberOfJets()>=2){
1066  Two_Hardest_SubJet_Distance=TMath::Sqrt(TMath::Power(Reclusterer->GetJet(SubJetOrdering(Jet1, Reclusterer, 1, 0, kTRUE))->Phi()-Reclusterer->GetJet(SubJetOrdering(Jet1, Reclusterer, 2, 0, kTRUE))->Phi(),2)+TMath::Power(Reclusterer->GetJet(SubJetOrdering(Jet1, Reclusterer, 1, 0, kTRUE))->Eta()-Reclusterer->GetJet(SubJetOrdering(Jet1, Reclusterer, 2, 0, kTRUE))->Eta(),2));
1067  if (Two_Hardest_SubJet_Distance>(2*fJetRadius)){//this can happen when the two subjet axis fall either side of the Phi boundary (i.e one close to zero the other close to 2*Pi IS THIS CORRECT?????????????????????????????????????????????????????????
1068  Two_Hardest_SubJet_Distance=TMath::Sqrt(TMath::Power((2*TMath::Pi())-TMath::Abs(Reclusterer->GetJet(SubJetOrdering(Jet1, Reclusterer, 1, 0, kTRUE))->Phi()-Reclusterer->GetJet(SubJetOrdering(Jet1,Reclusterer, 2, 0, kTRUE))->Phi()),2)+TMath::Power(Reclusterer->GetJet(SubJetOrdering(Jet1, Reclusterer, 1, 0, kTRUE))->Eta()-Reclusterer->GetJet(SubJetOrdering(Jet1, Reclusterer, 2,0, kTRUE))->Eta(),2));
1069  }
1070  }
1071  else Two_Hardest_SubJet_Distance=-2;
1072  fhSubJetiness2Distance->Fill(Two_Hardest_SubJet_Distance,SubJettiness2);
1073  if(fJetShapeSub==kNoSub) fShapesVar[0]= Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
1074  else fShapesVar[0]=Jet1->Pt();
1075  fShapesVar[1]=0;
1076  fShapesVar[2]=SubJetFraction(Jet1, Reclusterer, 1, 0, kTRUE, kFALSE);
1077  fShapesVar[3]=0;
1078  fhEventCounter->Fill(12);
1079  fShapesVar[4]=SubJetFraction(Jet1, Reclusterer, 2, 0, kTRUE, kFALSE);
1080  fShapesVar[5]=0;
1081  fShapesVar[6]=SubJettiness1;
1082  fShapesVar[7]=0;
1083  fShapesVar[8]=SubJettiness2;
1084  fShapesVar[9]=0;
1085  fShapesVar[10]=Two_Hardest_SubJet_Distance;
1086  fShapesVar[11]=0;
1087  fShapesVar[12]=0;
1088  fShapesVar[13]=0;
1089  fTreeResponseMatrixAxis->Fill();
1090  }
1091  }
1092  fhJetCounter->Fill(JetCounter); //Number of Jets in Each Event
1093  }
1094  //else {fhEventCounter->Fill(2);} //Events with no jet container
1095  }
1096  // if (JetCounter>40) cout << "Too Big!" << JetCounter <<endl;
1097 
1098  return kTRUE;
1099 }
1100 
1101 //________________________________________________________________________
1102 Double_t AliAnalysisTaskSubJetFraction::RelativePhi(Double_t Phi1, Double_t Phi2){
1103 
1104  if(Phi1 < -1*TMath::Pi()) Phi1 += (2*TMath::Pi()); // Turns the range of 0to2Pi into -PitoPi ???????????
1105  else if (Phi1 > TMath::Pi()) Phi1 -= (2*TMath::Pi());
1106  if(Phi2 < -1*TMath::Pi()) Phi2 += (2*TMath::Pi());
1107  else if (Phi2 > TMath::Pi()) Phi2 -= (2*TMath::Pi());
1108  Double_t DeltaPhi=Phi2-Phi1;
1109  if(DeltaPhi < -1*TMath::Pi()) DeltaPhi += (2*TMath::Pi());
1110  else if (DeltaPhi > TMath::Pi()) DeltaPhi -= (2*TMath::Pi());
1111  return DeltaPhi;
1112 }
1113 
1114 
1115 
1116 //--------------------------------------------------------------------------
1118 
1119  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1120  Double_t Angularity_Numerator=0; //Reset these values
1121  Double_t Angularity_Denominator=0;
1122  AliVParticle *Particle=0x0;
1123  Double_t DeltaPhi=0;
1124  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks (particles in the jet
1125  Particle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
1126  if(!Particle) continue;
1127  DeltaPhi=RelativePhi(Jet->Phi(),Particle->Phi());
1128  Angularity_Numerator=Angularity_Numerator+(Particle->Pt()*TMath::Sqrt(((Particle->Eta()-Jet->Eta())*(Particle->Eta()-Jet->Eta()))+(DeltaPhi*DeltaPhi)));
1129  Angularity_Denominator= Angularity_Denominator+Particle->Pt();
1130  }
1131  if(Angularity_Denominator!=0) return Angularity_Numerator/Angularity_Denominator;
1132  else return -1;
1133 }
1134 
1135 
1136 
1137 //--------------------------------------------------------------------------
1138 Double_t AliAnalysisTaskSubJetFraction::PTD(AliEmcalJet *Jet, Int_t JetContNb){
1139 
1140  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1141  Double_t PTD_Numerator=0; //Reset these values
1142  Double_t PTD_Denominator=0;
1143  AliVParticle *Particle=0x0;
1144  Double_t DeltaPhi=0;
1145  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks (particles in the jet
1146  Particle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
1147  if(!Particle) continue;
1148  PTD_Numerator=PTD_Numerator+(Particle->Pt()*Particle->Pt());
1149  PTD_Denominator=PTD_Denominator+Particle->Pt();
1150  }
1151  if(PTD_Denominator!=0) return TMath::Sqrt(PTD_Numerator)/PTD_Denominator;
1152  else return -1;
1153 }
1154 
1155 
1156 //----------------------------------------------------------------------
1158 AliEmcalJetFinder *AliAnalysisTaskSubJetFraction::Recluster(AliEmcalJet *Jet, Int_t JetContNb, Double_t SubJetRadius, Double_t SubJetMinPt, Int_t Algorithm, const char* Name){
1159 
1160  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1161  AliEmcalJetFinder *Reclusterer = new AliEmcalJetFinder(Name); //JetFinder Object for reclustered jets
1162  Reclusterer->SetRadius(SubJetRadius);
1163  Reclusterer->SetJetMinPt(SubJetMinPt);
1164  Reclusterer->SetJetAlgorithm(Algorithm); //0 for anti-kt 1 for kt
1165  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
1166  Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
1167  if(Reclusterer->AliEmcalJetFinder::Filter(Jet, JetCont, dVtx)){;} //reclustering jet1 using the jetfinderobject Reclusterer
1168  return Reclusterer;
1169 }
1170 
1171 
1172 
1173 
1174 
1175 
1176 
1177 
1178 
1179 //----------------------------------------------------------------------
1180 Double_t AliAnalysisTaskSubJetFraction::SubJetOrdering(AliEmcalJet *Jet, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t Type, Bool_t Index){
1181  AliEmcalJet *SubJet=NULL;
1182  Double_t SortingVariable;
1183  Int_t ArraySize =N+1;
1184  TArrayD *JetSorter = new TArrayD(ArraySize);
1185  TArrayD *JetIndexSorter = new TArrayD(ArraySize);
1186  for (Int_t i=0; i<ArraySize; i++){
1187  JetSorter->SetAt(0,i);
1188  }
1189  for (Int_t i=0; i<ArraySize; i++){
1190  JetIndexSorter->SetAt(0,i);
1191  }
1192  if(Reclusterer->GetNumberOfJets()<N) return -999;
1193  for (Int_t i=0; i<Reclusterer->GetNumberOfJets(); i++){
1194  SubJet=Reclusterer->GetJet(i);
1195  if (Type==0) SortingVariable=SubJet->Pt();
1196  else if (Type==1) SortingVariable=SubJet->E();
1197  else if (Type==2) SortingVariable=SubJet->M();
1198  for (Int_t j=0; j<N; j++){
1199  if (SortingVariable>JetSorter->GetAt(j)){
1200  for (Int_t k=N-1; k>=j; k--){
1201  JetSorter->SetAt(JetSorter->GetAt(k),k+1);
1202  JetIndexSorter->SetAt(JetIndexSorter->GetAt(k),k+1);
1203  }
1204  JetSorter->SetAt(SortingVariable,j);
1205  JetIndexSorter->SetAt(i,j);
1206  break;
1207  }
1208  }
1209  }
1210  if (!Index) return JetSorter->GetAt(N-1);
1211  else return JetIndexSorter->GetAt(N-1);
1212 }
1213 
1214 
1215 
1216 //returns -1 if the Nth hardest jet is requested where N>number of available jets
1217 //type: 0=Pt 1=E 2=M
1218 //Index TRUE=returns index FALSE=returns value of quantatiy in question
1219 
1220 
1221 
1222 
1223 
1224 
1225 
1226 //----------------------------------------------------------------------
1227 Double_t AliAnalysisTaskSubJetFraction::SubJetFraction(AliEmcalJet *Jet, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t Type, Bool_t Add, Bool_t Loss){
1228  AliEmcalJet *SubJet=NULL;
1229  Double_t Observable=0;
1230  Double_t Fraction_Numerator=0;
1231  Bool_t Error=kFALSE;
1232  if (!Jet->GetNumberOfTracks()) return -2;
1233  if (Add){
1234  for (Int_t i=1; i<=N; i++){
1235  Observable=SubJetOrdering(Jet,Reclusterer,i,Type,kFALSE);
1236  if(Observable==-999){
1237  Error = kTRUE;
1238  return -2;
1239  }
1240  Fraction_Numerator=Fraction_Numerator+Observable;
1241  }
1242  }
1243  else {
1244  Fraction_Numerator=SubJetOrdering(Jet,Reclusterer,N,Type,kFALSE);
1245  if (Fraction_Numerator==-999) return -2;
1246  }
1247  if (Type==0){
1248  if(Loss) return (Jet->Pt()-Fraction_Numerator)/Jet->Pt();
1249  else return Fraction_Numerator/Jet->Pt();
1250  }
1251  else if (Type==1){
1252  if(Loss) return (Jet->E()-Fraction_Numerator)/Jet->E();
1253  else return Fraction_Numerator/Jet->E();
1254  }
1255  else { //change to else if if you want to add more later
1256  if(Loss) return (Jet->M()-Fraction_Numerator)/Jet->M();
1257  else return Fraction_Numerator/Jet->M();
1258  }
1259 }
1260 //N number of hardest subjets involved
1261 //Type 0=Pt 1=E 2=M
1262 // Add TRUE: Add 1 to Nth hardest subjets together/total jet False: Nth hardest Jet/toal jet
1263 //Loss TRUE: Jet-Subjet(s)/Jet FALSE: Subjet(s)/Jet
1264 
1265 
1266 //----------------------------------------------------------------------
1267 Double_t AliAnalysisTaskSubJetFraction::NSubJettiness(AliEmcalJet *Jet, Int_t JetContNb, Double_t JetRadius, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t A, Int_t B){
1268  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1269  AliEmcalJet *SubJet=NULL;
1270  Double_t DeltaR1=0;
1271  Double_t DeltaR2=0;
1272  AliVParticle *JetParticle=0x0;
1273  Double_t SubJetiness_Numerator = 0;
1274  Double_t SubJetiness_Denominator = 0;
1275  Double_t Index=-2;
1276  Bool_t Error=kFALSE;
1277  if (!Jet->GetNumberOfTracks()) return -2;
1278  if (Reclusterer->GetNumberOfJets() < N) return -2;
1279  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks (particles in the jet
1280  JetParticle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
1281  for (Int_t j=1; j<=N; j++){
1282  Index=SubJetOrdering(Jet,Reclusterer,j,0,kTRUE);
1283  if(Index==-999){
1284  Error = kTRUE;
1285  i=Jet->GetNumberOfTracks();
1286  break;
1287  }
1288  if(j==1){
1289  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);
1290  }
1291  else{
1292  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);
1293  if (DeltaR2<DeltaR1) DeltaR1=DeltaR2;
1294  }
1295  }
1296  SubJetiness_Numerator=SubJetiness_Numerator+(JetParticle->Pt()*DeltaR1);
1297  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));
1298  else SubJetiness_Denominator=SubJetiness_Denominator+(TMath::Power((Reclusterer->GetJet(SubJetOrdering(Jet,Reclusterer,N,0,kTRUE))->Pt()),A)*JetParticle->Pt()*TMath::Power(JetRadius,B));
1299  }
1300  if (SubJetiness_Denominator!=0 && !Error) return SubJetiness_Numerator/SubJetiness_Denominator;
1301  else return -2;
1302 }
1303 
1304 
1305 //________________________________________________________________________
1307  //
1308  // retrieve event objects
1309  //
1311  return kFALSE;
1312 
1313  return kTRUE;
1314 }
1315 
1316 //_______________________________________________________________________
1318 {
1319  // Called once at the end of the analysis.
1321 
1322  /*
1323  for (int i=1; i<=(XBinsJetPtSize)-1; i++){ //Rescales the fhJetPt Histograms according to the with of the variable bins and number of events
1324  fhJetPt->SetBinContent(i,((1.0/(XBinsJetPt[i]-XBinsJetPt[i-1]))*((fhJetPt->GetBinContent(i))*(1.0/(fhEventCounter->GetBinContent(1))))));
1325  fhSubJetPt->SetBinContent(i,((1.0/(XBinsJetPt[i]-XBinsJetPt[i-1]))*((fhSubJetPt->GetBinContent(i))*(1.0/(fhEventCounter->GetBinContent(5))))));
1326 
1327  //fhJetPt->SetBinContent(i,((1.0/(fhPt->GetBinWidth(i)))*((fhJetPt->GetBinContent(i))*(fhEventCounter->GetBinContent(1)))));
1328  // fhJetPt->SetBinContent(i,((1.0/((fhPt->GetLowEdge(i+1))-(fhJetPt->GetLowEdge(i))))*((fhPt->GetBinContent(i))*(fhEventCounter->GetBinContent(1)))));
1329  }
1330  for (int i=1; i<=(XBinsJetMassSize)-1; i++){ //Rescales the fhPt Histograms according to the with of the variable bins and number of events
1331  fhJetMass->SetBinContent(i,((1.0/(XBinsJetMass[i]-XBinsJetMass[i-1]))*((fhJetMass->GetBinContent(i))*(1.0/(fhEventCounter->GetBinContent(1))))));
1332  }
1333 
1334  fhJetPhi->Scale(Phi_Bins/((Phi_Up-Phi_Low)*((fhEventCounter->GetBinContent(1)))));
1335  fhJetEta->Scale(Eta_Bins/((Eta_Up-Eta_Low)*((fhEventCounter->GetBinContent(1)))));
1336  fhJetRadius->Scale(100/(fhEventCounter->GetBinContent(4))); //should this and JetAngularity be divided by Bin 1 or 4????
1337  fhJetAngularity->Scale(100/(fhEventCounter->GetBinContent(4)));
1338  fhJetPTD->Scale(100/(fhEventCounter->GetBinContent(4)));
1339  fhNumberOfJetTracks->Scale(1.0/(fhEventCounter->GetBinContent(4)));
1340  fhSubJetPtFrac->Scale(100/(fhEventCounter->GetBinContent(8)));
1341  fhSubJetPtFrac2->Scale(100/(fhEventCounter->GetBinContent(9)));
1342  fhSubJetEnergyFrac->Scale(100/(fhEventCounter->GetBinContent(8)));
1343  fhSubJetEnergyFrac2->Scale(100/(fhEventCounter->GetBinContent(9)));
1344  fhJetCounter->Scale(1.0/(fhEventCounter->GetBinContent(1))); //is the first bin the correct one to look at?
1345 
1346 
1347  fhSubJetCounter->Scale(1.0/(fhEventCounter->GetBinContent(5)));
1348  fhSubJetiness1->Scale(100/(fhEventCounter->GetBinContent(12)));
1349  fhSubJetiness2->Scale(100/(fhEventCounter->GetBinContent(12)));
1350  fh2to1SubJetinessRatio->Scale(100/(fhEventCounter->GetBinContent(12)));
1351 
1352  */
1353  /*
1354 
1355  fhJetPt->Scale(1.0/(fhEventCounter->GetBinContent(1)));
1356  fhSubJetPt->Scale(1.0/(fhEventCounter->GetBinContent(5)));
1357  fhJetMass->Scale(1.0/(fhEventCounter->GetBinContent(1)));
1358  fhJetPhi->Scale(1.0/(fhEventCounter->GetBinContent(1)));
1359  fhJetEta->Scale(1.0/(fhEventCounter->GetBinContent(1)));
1360  fhJetRadius->Scale(1.0/(fhEventCounter->GetBinContent(4))); //should this and JetAngularity be divided by Bin 1 or 4????
1361  fhJetAngularity->Scale(1.0/(fhEventCounter->GetBinContent(4)));
1362  fhJetPTD->Scale(1.0/(fhEventCounter->GetBinContent(4)));
1363  fhNumberOfJetTracks->Scale(1.0/(fhEventCounter->GetBinContent(4)));
1364  fhSubJetPtFrac->Scale(1.0/(fhEventCounter->GetBinContent(8)));
1365  fhSubJetPtFrac2->Scale(1.0/(fhEventCounter->GetBinContent(9)));
1366  fhSubJetEnergyFrac->Scale(1.0/(fhEventCounter->GetBinContent(8)));
1367  fhSubJetEnergyFrac2->Scale(1.0/(fhEventCounter->GetBinContent(9)));
1368  fhJetCounter->Scale(1.0/(fhEventCounter->GetBinContent(1))); //is the first bin the correct one to look at?
1369  fhSubJetCounter->Scale(1.0/(fhEventCounter->GetBinContent(5)));
1370  fhSubJetiness1->Scale(1.0/(fhEventCounter->GetBinContent(12)));
1371  fhSubJetiness2->Scale(1.0/(fhEventCounter->GetBinContent(12)));
1372  fh2to1SubJetinessRatio->Scale(1.0/(fhEventCounter->GetBinContent(12)));
1373  */
1374  }
1375 
1376 }
void SetRadius(Double_t val)
Double_t Area() const
Definition: AliEmcalJet.h:69
Double_t XBinsPtSize
void SetJetMinPt(Double_t val)
static Double_t DeltaPhi(Double_t phia, Double_t phib, Double_t rMin=-TMath::Pi()/2, Double_t rMax=3 *TMath::Pi()/2)
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:158
AliJetContainer * GetJetContainer(Int_t i=0) const
Double_t Eta() const
Definition: AliEmcalJet.h:60
Double_t Phi() const
Definition: AliEmcalJet.h:55
Double_t PTD(AliEmcalJet *Jet, Int_t JetContNb)
Int_t GetLabel() const
Definition: AliEmcalJet.h:63
Double_t NSubJettiness(AliEmcalJet *Jet, Int_t JetContNb, Double_t JetRadius, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t A, Int_t B)
Double_t RelativePhi(Double_t Phi1, Double_t Phi2)
Double_t E() const
Definition: AliEmcalJet.h:58
TList * fOutput
!output list
Double_t XBinsJetMass[28]
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:83
TString kData
Declare data MC or deltaAOD.
AliParticleContainer * GetParticleContainer() const
void SetJetAlgorithm(Int_t val)
ClassImp(AliAnalysisTaskSubJetFraction) AliAnalysisTaskSubJetFraction
Int_t GetNJets() const
Double_t fCent
!event centrality
TClonesArray * GetArray() const
Double_t SubJetFraction(AliEmcalJet *Jet, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t Type, Bool_t Add, Bool_t Loss)
AliEmcalJet * GetNextAcceptJet()
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
AliEmcalJetFinder * Recluster(AliEmcalJet *Jet, Int_t JetContNb, Double_t SubJetRadius, Double_t SubJetMinPt, Int_t Algorithm, const char *Name)
AliEmcalJet * GetJet(Int_t index)
Double_t Pt() const
Definition: AliEmcalJet.h:47
Double_t GetRhoVal(Int_t i=0) const
Double_t XBinsPt[66]
Short_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:107
void SetMakeGeneralHistograms(Bool_t g)
Declaration of class AliEmcalPythiaInfo.
Double_t Angularity(AliEmcalJet *Jet, Int_t JetContNb)
Double_t SubJetOrdering(AliEmcalJet *Jet, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t Type, Bool_t Index)
Double_t XBinsJetPt[35]
Double_t M() const
Definition: AliEmcalJet.h:59
void ResetCurrentID(Int_t i=-1)
AliEmcalJet * GetJet(Int_t i) const