AliPhysics  d497afb (d497afb)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskSubJetFraction.cxx
Go to the documentation of this file.
1 //
2 // Subjet fraction analysis task.
3 //
4 //Nima Zardoshti
5 
6 #include <TClonesArray.h>
7 #include <TH1F.h>
8 #include <TH2F.h>
9 #include <TH3F.h>
10 #include <TCanvas.h>
11 #include <THnSparse.h>
12 #include <TTree.h>
13 #include <TList.h>
14 #include <TLorentzVector.h>
15 #include <TProfile.h>
16 #include <TChain.h>
17 #include <TSystem.h>
18 #include <TFile.h>
19 #include <TKey.h>
20 #include <AliAnalysisDataSlot.h>
21 #include <AliAnalysisDataContainer.h>
22 #include "TMatrixD.h"
23 #include "TMatrixDSym.h"
24 #include "TMatrixDSymEigen.h"
25 #include "TVector3.h"
26 #include "TVector2.h"
27 #include "AliVCluster.h"
28 #include "AliVTrack.h"
29 #include "AliEmcalJet.h"
30 #include "AliRhoParameter.h"
31 #include "AliLog.h"
32 #include "AliEmcalParticle.h"
33 #include "AliMCEvent.h"
34 #include "AliGenPythiaEventHeader.h"
35 #include "AliAODMCHeader.h"
36 #include "AliMCEvent.h"
37 #include "AliAnalysisManager.h"
38 #include "AliJetContainer.h"
39 #include "AliParticleContainer.h"
40 //#include "AliPythiaInfo.h"
41 #include "TRandom3.h"
42 #include "AliPicoTrack.h"
43 #include "AliEmcalJetFinder.h"
44 #include "AliAODEvent.h"
46 
47 #include "FJ_includes.h"
48 
49 //Globals
50 
51 
52 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
54 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
56 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
61 Double_t Phi_Up=2*(TMath::Pi());
62 Double_t Phi_Low=(-1*(TMath::Pi()));
64 
65 
66 
67 
68 
69 using std::cout;
70 using std::endl;
71 
73 
74 //________________________________________________________________________
77  fContainer(0),
78  fMinFractionShared(0),
79  fJetShapeType(kData),
80  fJetShapeSub(kNoSub),
81  fJetSelection(kInclusive),
82  fPtThreshold(-9999.),
83  fRMatching(0.2),
84  fPtMinTriggerHadron(20.),
85  fPtMaxTriggerHadron(50.),
86  fRecoilAngularWindow(0.6),
87  fSemigoodCorrect(0),
88  fHolePos(0),
89  fHoleWidth(0),
90  fRandom(0x0),
91  fCentSelectOn(kTRUE),
92  fCentMin(0),
93  fCentMax(10),
94  fSubJetAlgorithm(0),
95  fSubJetRadius(0.1),
96  fSubJetMinPt(1),
97  fJetRadius(0.4),
98  fRMatched(0.2),
99  fSharedFractionPtMin(0.5),
100  fDerivSubtrOrder(0),
101  fFullTree(kFALSE),
102  fBeta_SD(0),
103  fZCut(0.1),
104  fReclusteringAlgorithm(0),
105  fSoftDropOn(0),
106  fNsubMeasure(kFALSE),
107  fhPtTriggerHadron(0x0),
108  fhJetPt(0x0),
109  fhJetPt_1(0x0),
110  fhJetPt_2(0x0),
111  fhJetPhi(0x0),
112  fhJetPhi_1(0x0),
113  fhJetPhi_2(0x0),
114  fhJetEta(0x0),
115  fhJetEta_1(0x0),
116  fhJetEta_2(0x0),
117  fhJetMass(0x0),
118  fhJetMass_1(0x0),
119  fhJetMass_2(0x0),
120  fhJetRadius(0x0),
121  fhJetRadius_1(0x0),
122  fhJetRadius_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  fh2PtTriggerHadronJet(0x0),
145  fhPhiTriggerHadronJet(0x0),
146  fhPhiTriggerHadronEventPlane(0x0),
147  fhPhiTriggerHadronEventPlaneTPC(0x0),
148  fh2PtRatio(0x0),
149  fhEventCounter(0x0),
150  fhEventCounter_1(0x0),
151  fhEventCounter_2(0x0),
152  fhSubJettiness1CheckRatio_FJ_AKT(0x0),
153  fhSubJettiness1CheckRatio_FJ_KT(0x0),
154  fhSubJettiness1CheckRatio_FJ_CA(0x0),
155  fhSubJettiness1CheckRatio_FJ_WTA_KT(0x0),
156  fhSubJettiness1CheckRatio_FJ_WTA_CA(0x0),
157  fhSubJettiness1CheckRatio_FJ_OP_AKT(0x0),
158  fhSubJettiness1CheckRatio_FJ_OP_KT(0x0),
159  fhSubJettiness1CheckRatio_FJ_OP_CA(0x0),
160  fhSubJettiness1CheckRatio_FJ_OP_WTA_KT(0x0),
161  fhSubJettiness1CheckRatio_FJ_OP_WTA_CA(0x0),
162  fhSubJettiness1CheckRatio_FJ_MIN(0x0),
163  fhSubJettiness2CheckRatio_FJ_AKT(0x0),
164  fhSubJettiness2CheckRatio_FJ_KT(0x0),
165  fhSubJettiness2CheckRatio_FJ_CA(0x0),
166  fhSubJettiness2CheckRatio_FJ_WTA_KT(0x0),
167  fhSubJettiness2CheckRatio_FJ_WTA_CA(0x0),
168  fhSubJettiness2CheckRatio_FJ_OP_AKT(0x0),
169  fhSubJettiness2CheckRatio_FJ_OP_KT(0x0),
170  fhSubJettiness2CheckRatio_FJ_OP_CA(0x0),
171  fhSubJettiness2CheckRatio_FJ_OP_WTA_KT(0x0),
172  fhSubJettiness2CheckRatio_FJ_OP_WTA_CA(0x0),
173  fhSubJettiness2CheckRatio_FJ_MIN(0x0),
174  fhSubJettiness2to1CheckRatio_FJ_AKT(0x0),
175  fhSubJettiness2to1CheckRatio_FJ_KT(0x0),
176  fhSubJettiness2to1CheckRatio_FJ_CA(0x0),
177  fhSubJettiness2to1CheckRatio_FJ_WTA_KT(0x0),
178  fhSubJettiness2to1CheckRatio_FJ_WTA_CA(0x0),
179  fhSubJettiness2to1CheckRatio_FJ_OP_AKT(0x0),
180  fhSubJettiness2to1CheckRatio_FJ_OP_KT(0x0),
181  fhSubJettiness2to1CheckRatio_FJ_OP_CA(0x0),
182  fhSubJettiness2to1CheckRatio_FJ_OP_WTA_KT(0x0),
183  fhSubJettiness2to1CheckRatio_FJ_OP_WTA_CA(0x0),
184  fhSubJettiness2to1CheckRatio_FJ_MIN(0x0),
185  fhSubJettiness1_FJ_AKT(0x0),
186  fhSubJettiness1_FJ_KT(0x0),
187  fhSubJettiness1_FJ_CA(0x0),
188  fhSubJettiness1_FJ_WTA_KT(0x0),
189  fhSubJettiness1_FJ_WTA_CA(0x0),
190  fhSubJettiness1_FJ_OP_AKT(0x0),
191  fhSubJettiness1_FJ_OP_KT(0x0),
192  fhSubJettiness1_FJ_OP_CA(0x0),
193  fhSubJettiness1_FJ_OP_WTA_KT(0x0),
194  fhSubJettiness1_FJ_OP_WTA_CA(0x0),
195  fhSubJettiness1_FJ_MIN(0x0),
196  fhSubJettiness2_FJ_AKT(0x0),
197  fhSubJettiness2_FJ_KT(0x0),
198  fhSubJettiness2_FJ_CA(0x0),
199  fhSubJettiness2_FJ_WTA_KT(0x0),
200  fhSubJettiness2_FJ_WTA_CA(0x0),
201  fhSubJettiness2_FJ_OP_AKT(0x0),
202  fhSubJettiness2_FJ_OP_KT(0x0),
203  fhSubJettiness2_FJ_OP_CA(0x0),
204  fhSubJettiness2_FJ_OP_WTA_KT(0x0),
205  fhSubJettiness2_FJ_OP_WTA_CA(0x0),
206  fhSubJettiness2_FJ_MIN(0x0),
207  fhSubJettiness2to1_FJ_AKT(0x0),
208  fhSubJettiness2to1_FJ_KT(0x0),
209  fhSubJettiness2to1_FJ_CA(0x0),
210  fhSubJettiness2to1_FJ_WTA_KT(0x0),
211  fhSubJettiness2to1_FJ_WTA_CA(0x0),
212  fhSubJettiness2to1_FJ_OP_AKT(0x0),
213  fhSubJettiness2to1_FJ_OP_KT(0x0),
214  fhSubJettiness2to1_FJ_OP_CA(0x0),
215  fhSubJettiness2to1_FJ_OP_WTA_KT(0x0),
216  fhSubJettiness2to1_FJ_OP_WTA_CA(0x0),
217  fhSubJettiness2to1_FJ_MIN(0x0),
218  fTreeResponseMatrixAxis(0)
219 
220 {
221  for(Int_t i=0;i<nVar;i++){
222  fShapesVar[i]=0;
223  }
224  SetMakeGeneralHistograms(kTRUE);
225  DefineOutput(1, TList::Class());
226  DefineOutput(2, TTree::Class());
227 }
228 
229 //________________________________________________________________________
231  AliAnalysisTaskEmcalJet(name, kTRUE),
232  fContainer(0),
233  fMinFractionShared(0),
234  fJetShapeType(kData),
235  fJetShapeSub(kNoSub),
236  fJetSelection(kInclusive),
237  fPtThreshold(-9999.),
238  fRMatching(0.2),
239  fPtMinTriggerHadron(20.),
240  fPtMaxTriggerHadron(50.),
241  fRecoilAngularWindow(0.6),
242  fSemigoodCorrect(0),
243  fHolePos(0),
244  fHoleWidth(0),
245  fRandom(0x0),
246  fCentSelectOn(kTRUE),
247  fCentMin(0),
248  fCentMax(10),
249  fSubJetAlgorithm(0),
250  fSubJetRadius(0.1),
251  fSubJetMinPt(1),
252  fJetRadius(0.4),
253  fRMatched(0.2),
254  fSharedFractionPtMin(0.5),
255  fDerivSubtrOrder(0),
256  fFullTree(kFALSE),
257  fBeta_SD(0),
258  fZCut(0.1),
259  fReclusteringAlgorithm(0),
260  fSoftDropOn(0),
261  fNsubMeasure(kFALSE),
262  fhPtTriggerHadron(0x0),
263  fhJetPt(0x0),
264  fhJetPt_1(0x0),
265  fhJetPt_2(0x0),
266  fhJetPhi(0x0),
267  fhJetPhi_1(0x0),
268  fhJetPhi_2(0x0),
269  fhJetEta(0x0),
270  fhJetEta_1(0x0),
271  fhJetEta_2(0x0),
272  fhJetMass(0x0),
273  fhJetMass_1(0x0),
274  fhJetMass_2(0x0),
275  fhJetRadius(0x0),
276  fhJetRadius_1(0x0),
277  fhJetRadius_2(0x0),
278  fhJetCounter(0x0),
279  fhJetCounter_1(0x0),
280  fhJetCounter_2(0x0),
281  fhNumberOfJetTracks(0x0),
282  fhNumberOfJetTracks_1(0x0),
283  fhNumberOfJetTracks_2(0x0),
284  fhSubJetPt(0x0),
285  fhSubJetPt_1(0x0),
286  fhSubJetPt_2(0x0),
287  fhSubJetMass(0x0),
288  fhSubJetMass_1(0x0),
289  fhSubJetMass_2(0x0),
290  fhSubJetRadius(0x0),
291  fhSubJetRadius_1(0x0),
292  fhSubJetRadius_2(0x0),
293  fhSubJetCounter(0x0),
294  fhSubJetCounter_1(0x0),
295  fhSubJetCounter_2(0x0),
296  fhNumberOfSubJetTracks(0x0),
297  fhNumberOfSubJetTracks_1(0x0),
298  fhNumberOfSubJetTracks_2(0x0),
299  fh2PtTriggerHadronJet(0x0),
300  fhPhiTriggerHadronJet(0x0),
301  fhPhiTriggerHadronEventPlane(0x0),
302  fhPhiTriggerHadronEventPlaneTPC(0x0),
303  fh2PtRatio(0x0),
304  fhEventCounter(0x0),
305  fhEventCounter_1(0x0),
306  fhEventCounter_2(0x0),
307  fhSubJettiness1CheckRatio_FJ_AKT(0x0),
308  fhSubJettiness1CheckRatio_FJ_KT(0x0),
309  fhSubJettiness1CheckRatio_FJ_CA(0x0),
310  fhSubJettiness1CheckRatio_FJ_WTA_KT(0x0),
311  fhSubJettiness1CheckRatio_FJ_WTA_CA(0x0),
312  fhSubJettiness1CheckRatio_FJ_OP_AKT(0x0),
313  fhSubJettiness1CheckRatio_FJ_OP_KT(0x0),
314  fhSubJettiness1CheckRatio_FJ_OP_CA(0x0),
315  fhSubJettiness1CheckRatio_FJ_OP_WTA_KT(0x0),
316  fhSubJettiness1CheckRatio_FJ_OP_WTA_CA(0x0),
317  fhSubJettiness1CheckRatio_FJ_MIN(0x0),
318  fhSubJettiness2CheckRatio_FJ_AKT(0x0),
319  fhSubJettiness2CheckRatio_FJ_KT(0x0),
320  fhSubJettiness2CheckRatio_FJ_CA(0x0),
321  fhSubJettiness2CheckRatio_FJ_WTA_KT(0x0),
322  fhSubJettiness2CheckRatio_FJ_WTA_CA(0x0),
323  fhSubJettiness2CheckRatio_FJ_OP_AKT(0x0),
324  fhSubJettiness2CheckRatio_FJ_OP_KT(0x0),
325  fhSubJettiness2CheckRatio_FJ_OP_CA(0x0),
326  fhSubJettiness2CheckRatio_FJ_OP_WTA_KT(0x0),
327  fhSubJettiness2CheckRatio_FJ_OP_WTA_CA(0x0),
328  fhSubJettiness2CheckRatio_FJ_MIN(0x0),
329  fhSubJettiness2to1CheckRatio_FJ_AKT(0x0),
330  fhSubJettiness2to1CheckRatio_FJ_KT(0x0),
331  fhSubJettiness2to1CheckRatio_FJ_CA(0x0),
332  fhSubJettiness2to1CheckRatio_FJ_WTA_KT(0x0),
333  fhSubJettiness2to1CheckRatio_FJ_WTA_CA(0x0),
334  fhSubJettiness2to1CheckRatio_FJ_OP_AKT(0x0),
335  fhSubJettiness2to1CheckRatio_FJ_OP_KT(0x0),
336  fhSubJettiness2to1CheckRatio_FJ_OP_CA(0x0),
337  fhSubJettiness2to1CheckRatio_FJ_OP_WTA_KT(0x0),
338  fhSubJettiness2to1CheckRatio_FJ_OP_WTA_CA(0x0),
339  fhSubJettiness2to1CheckRatio_FJ_MIN(0x0),
340  fhSubJettiness1_FJ_AKT(0x0),
341  fhSubJettiness1_FJ_KT(0x0),
342  fhSubJettiness1_FJ_CA(0x0),
343  fhSubJettiness1_FJ_WTA_KT(0x0),
344  fhSubJettiness1_FJ_WTA_CA(0x0),
345  fhSubJettiness1_FJ_OP_AKT(0x0),
346  fhSubJettiness1_FJ_OP_KT(0x0),
347  fhSubJettiness1_FJ_OP_CA(0x0),
348  fhSubJettiness1_FJ_OP_WTA_KT(0x0),
349  fhSubJettiness1_FJ_OP_WTA_CA(0x0),
350  fhSubJettiness1_FJ_MIN(0x0),
351  fhSubJettiness2_FJ_AKT(0x0),
352  fhSubJettiness2_FJ_KT(0x0),
353  fhSubJettiness2_FJ_CA(0x0),
354  fhSubJettiness2_FJ_WTA_KT(0x0),
355  fhSubJettiness2_FJ_WTA_CA(0x0),
356  fhSubJettiness2_FJ_OP_AKT(0x0),
357  fhSubJettiness2_FJ_OP_KT(0x0),
358  fhSubJettiness2_FJ_OP_CA(0x0),
359  fhSubJettiness2_FJ_OP_WTA_KT(0x0),
360  fhSubJettiness2_FJ_OP_WTA_CA(0x0),
361  fhSubJettiness2_FJ_MIN(0x0),
362  fhSubJettiness2to1_FJ_AKT(0x0),
363  fhSubJettiness2to1_FJ_KT(0x0),
364  fhSubJettiness2to1_FJ_CA(0x0),
365  fhSubJettiness2to1_FJ_WTA_KT(0x0),
366  fhSubJettiness2to1_FJ_WTA_CA(0x0),
367  fhSubJettiness2to1_FJ_OP_AKT(0x0),
368  fhSubJettiness2to1_FJ_OP_KT(0x0),
369  fhSubJettiness2to1_FJ_OP_CA(0x0),
370  fhSubJettiness2to1_FJ_OP_WTA_KT(0x0),
371  fhSubJettiness2to1_FJ_OP_WTA_CA(0x0),
372  fhSubJettiness2to1_FJ_MIN(0x0),
373  fTreeResponseMatrixAxis(0)
374 
375 {
376  // Standard constructor.
377  for(Int_t i=0;i<nVar;i++){
378  fShapesVar[i]=0;
379  }
381  DefineOutput(1, TList::Class());
382  DefineOutput(2, TTree::Class());
383 }
384 
385 //________________________________________________________________________
387 {
388  // Destructor.
389 }
390 
391 //________________________________________________________________________
393 {
394  // Create user output.
395 
397 
398  Bool_t oldStatus = TH1::AddDirectoryStatus();
399  TH1::AddDirectory(kFALSE);
400  TH1::AddDirectory(oldStatus);
401  //create a tree used for the MC data and making a 4D response matrix
402  const char* nameoutput = GetOutputSlot(2)->GetContainer()->GetName();
403  fTreeResponseMatrixAxis = new TTree(nameoutput, nameoutput);
404  if (fFullTree){
405  TString *fShapesVarNames = new TString [nVar];
406 
407  fShapesVarNames[0] = "Pt";
408  fShapesVarNames[1] = "Pt_Truth";
409  fShapesVarNames[2] = "Tau1";
410  fShapesVarNames[3] = "Tau1_Truth";
411  fShapesVarNames[4] = "Tau2";
412  fShapesVarNames[5] = "Tau2_Truth";
413  fShapesVarNames[6] = "Tau3";
414  fShapesVarNames[7] = "Tau3_Truth";
415  fShapesVarNames[8] = "OpeningAngle";
416  fShapesVarNames[9] = "OpeningAngle_Truth";
417  fShapesVarNames[10] = "JetMultiplicity";
418  fShapesVarNames[11] = "JetMultiplicity_Truth";
419  fShapesVarNames[12] = "OpeningAngleSD";
420  fShapesVarNames[13] = "OpeningAngleSD_Truth";
421  fShapesVarNames[14] = "Zg";
422  fShapesVarNames[15] = "Zg_Truth";
423  fShapesVarNames[16] = "LeadingTrackPt";
424  fShapesVarNames[17] = "LeadingTrackPt_Truth";
425  fShapesVarNames[18] = "EventPlaneTriggerHadron";
426  fShapesVarNames[19] = "EventPlaneTriggerHadron_Truth";
427  fShapesVarNames[20] = "EventPlaneTPCTriggerHadron";
428  fShapesVarNames[21] = "EventPlaneTPCTriggerHadron_Truth";
429  fShapesVarNames[22] = "DeltaR";
430  fShapesVarNames[23] = "DeltaR_Truth";
431  fShapesVarNames[24] = "Frac1";
432  fShapesVarNames[25] = "Frac1_Truth";
433  fShapesVarNames[26] = "Frac2";
434  fShapesVarNames[27] = "Frac2_Truth";
435  for(Int_t ivar=0; ivar < nVar; ivar++){
436  cout<<"looping over variables"<<endl;
437  fTreeResponseMatrixAxis->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/D", fShapesVarNames[ivar].Data()));
438  }
439  }
440 
441  if (!fFullTree){
442  const Int_t nVarMin = 22;
443  TString *fShapesVarNames = new TString [nVarMin];
444 
445  fShapesVarNames[0] = "Pt";
446  fShapesVarNames[1] = "Pt_Truth";
447  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[1] = "Tau1_Randomisation";
448  fShapesVarNames[2] = "Tau1";
449  fShapesVarNames[3] = "Tau1_Truth";
450  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[3] = "Tau1_ExtraProng_02_15";
451  fShapesVarNames[4] = "Tau2";
452  fShapesVarNames[5] = "Tau2_Truth";
453  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[5] = "Tau1_ExtraProng_02_30";
454  fShapesVarNames[6] = "SubJet1LeadingTrackPt";
455  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[6] = "Tau1_ExtraProng_02_45";
456  fShapesVarNames[7] = "SubJet1LeadingTrackPt_Truth";
457  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[7] = "Tau1_ExtraProng_02_60";
458  fShapesVarNames[8] = "OpeningAngle";
459  fShapesVarNames[9] = "OpeningAngle_Truth";
460  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[9] = "Tau1_ExtraProng_01_30";
461  fShapesVarNames[10] = "SubJet2LeadingTrackPt";
462  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[10] = "Tau1_ExtraProng_01_45";
463  fShapesVarNames[11] = "SubJet2LeadingTrackPt_Truth";
464  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[11] = "Tau1_ExtraProng_03_30";
465  fShapesVarNames[12] = "OpeningAngleSD";
466  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[12] = "Tau1_ExtraProng_03_45";
467  fShapesVarNames[13] = "OpeningAngleSD_Truth";
468  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[13] = "Tau2_Randomisation";
469  fShapesVarNames[14] = "SubJet1Pt";
470  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[14] = "Tau2_ExtraProng_02_15";
471  fShapesVarNames[15] = "SubJet1Pt_Truth";
472  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[15] = "Tau2_ExtraProng_02_30";
473  fShapesVarNames[16] = "LeadingTrackPt";
474  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[16] = "Tau2_ExtraProng_02_45";
475  fShapesVarNames[17] = "LeadingTrackPt_Truth";
476  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[17] = "Tau2_ExtraProng_02_60";
477  fShapesVarNames[18] = "EventPlaneTriggerHadron";
478  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[18] = "Tau2_ExtraProng_01_30";
479  fShapesVarNames[19] = "EventPlaneTriggerHadron_Truth";
480  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[19] = "Tau2_ExtraProng_01_45";
481  fShapesVarNames[20] = "SubJet2Pt";
482  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[20] = "Tau2_ExtraProng_03_30";
483  fShapesVarNames[21] = "SubJet2Pt_Truth";
484  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[21] = "Tau2_ExtraProng_03_45";
485 
486  for(Int_t ivar=0; ivar < nVarMin; ivar++){
487  cout<<"looping over variables"<<endl;
488  fTreeResponseMatrixAxis->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/D", fShapesVarNames[ivar].Data()));
489  }
490  }
491 
493  fhPtTriggerHadron= new TH1F("fhPtTriggerHadron", "fhPtTriggerHadron",1500,-0.5,149.5);
495  fh2PtTriggerHadronJet= new TH2F("fh2PtTriggerHadronJet", "fh2PtTriggerHadronJet",1500,-0.5,149.5,1500,-0.5,149.5);
497  fhPhiTriggerHadronJet= new TH1F("fhPhiTriggerHadronJet", "fhPhiTriggerHadronJet",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
499  fhPhiTriggerHadronEventPlane= new TH1F("fhPhiTriggerHadronEventPlane", "fhPhiTriggerHadronEventPlane",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
501  fhPhiTriggerHadronEventPlaneTPC= new TH1F("fhPhiTriggerHadronEventPlaneTPC", "fhPhiTriggerHadronEventPlaneTPC",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
503  }
505 
506  // fhJetPt= new TH1F("fhJetPt", "Jet Pt", (XBinsJetPtSize)-1, XBinsJetPt);
507  fhJetPt= new TH1F("fhJetPt", "Jet Pt",1500,-0.5,149.5 );
508  fOutput->Add(fhJetPt);
509  //fhJetPhi= new TH1F("fhJetPhi", "Jet Phi", Phi_Bins, Phi_Low, Phi_Up);
510  fhJetPhi= new TH1F("fhJetPhi", "Jet Phi",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
511  fOutput->Add(fhJetPhi);
512  fhJetEta= new TH1F("fhJetEta", "Jet Eta", Eta_Bins, Eta_Low, Eta_Up);
513  fOutput->Add(fhJetEta);
514  fhJetMass= new TH1F("fhJetMass", "Jet Mass", 4000,-0.5, 39.5);
515  fOutput->Add(fhJetMass);
516  fhJetRadius= new TH1F("fhJetRadius", "Jet Radius", 100, -0.05,0.995);
517  fOutput->Add(fhJetRadius);
518  fhNumberOfJetTracks= new TH1F("fhNumberOfJetTracks", "Number of Tracks within a Jet", 300, -0.5,299.5);
520  fhSubJetRadius= new TH1F("fhSubJetRadius", "SubJet Radius", 100, -0.05,0.995);
521  fOutput->Add(fhSubJetRadius);
522  fhSubJetPt= new TH1F("fhSubJetPt", "SubJet Pt", 1500, -0.5,149.5);
523  fOutput->Add(fhSubJetPt);
524  fhSubJetMass= new TH1F("fhSubJetMass", "Sub Jet Mass", 4000,-0.5, 39.5);
525  fOutput->Add(fhSubJetMass);
526  fhNumberOfSubJetTracks= new TH1F("fhNumberOfSubJetTracks", "Number of Tracks within a Sub Jet", 300, -0.5,299.5);
528  fhJetCounter= new TH1F("fhJetCounter", "Jet Counter", 150, -0.5, 149.5);
529  fOutput->Add(fhJetCounter);
530  fhSubJetCounter = new TH1F("fhSubJetCounter", "SubJet Counter",50, -0.5,49.5);
531  fOutput->Add(fhSubJetCounter);
532  fhEventCounter= new TH1F("fhEventCounter", "Event Counter", 15,0.5,15.5);
533  fOutput->Add(fhEventCounter);
534  }
536  fhSubJettiness1_FJ_KT= new TH1D("fhSubJettiness1_FJ_KT","fhSubJettiness1_FJ_KT",400,-2,2);
538  fhSubJettiness1_FJ_MIN= new TH1D("fhSubJettiness1_FJ_MIN","fhSubJettiness1_FJ_MIN",400,-2,2);
540  fhSubJettiness2_FJ_KT= new TH1D("fhSubJettiness2_FJ_KT","fhSubJettiness2_FJ_KT",400,-2,2);
542  fhSubJettiness2_FJ_MIN= new TH1D("fhSubJettiness2_FJ_MIN","fhSubJettiness2_FJ_MIN",400,-2,2);
544  fhSubJettiness1CheckRatio_FJ_AKT = new TH2D("fhSubJettiness1CheckRatio_FJ_AKT","fhSubJettiness1CheckRatio_FJ_AKT",400,-2,2,300,-1,2);
546  fhSubJettiness1CheckRatio_FJ_KT= new TH2D("fhSubJettiness1CheckRatio_FJ_KT","fhSubJettiness1CheckRatio_FJ_KT",400,-2,2,300,-1,2);
548  fhSubJettiness1CheckRatio_FJ_CA= new TH2D("fhSubJettiness1CheckRatio_FJ_CA","fhSubJettiness1CheckRatio_FJ_CA",400,-2,2,300,-1,2);
550  fhSubJettiness1CheckRatio_FJ_WTA_KT= new TH2D("fhSubJettiness1CheckRatio_FJ_WTA_KT","fhSubJettiness1CheckRatio_FJ_WTA_KT",400,-2,2,300,-1,2);
552  fhSubJettiness1CheckRatio_FJ_WTA_CA= new TH2D("fhSubJettiness1CheckRatio_FJ_WTA_CA","fhSubJettiness1CheckRatio_FJ_WTA_CA",400,-2,2,300,-1,2);
554  fhSubJettiness1CheckRatio_FJ_OP_AKT= new TH2D("fhSubJettiness1CheckRatio_FJ_OP_AKT","fhSubJettiness1CheckRatio_FJ_OP_AKT",400,-2,2,300,-1,2);
556  fhSubJettiness1CheckRatio_FJ_OP_KT= new TH2D("fhSubJettiness1CheckRatio_FJ_OP_KT","fhSubJettiness1CheckRatio_FJ_OP_KT",400,-2,2,300,-1,2);
558  fhSubJettiness1CheckRatio_FJ_OP_CA= new TH2D("fhSubJettiness1CheckRatio_FJ_OP_CA","fhSubJettiness1CheckRatio_FJ_OP_CA",400,-2,2,300,-1,2);
560  fhSubJettiness1CheckRatio_FJ_OP_WTA_KT= new TH2D("fhSubJettiness1CheckRatio_FJ_OP_WTA_KT","fhSubJettiness1CheckRatio_FJ_OP_WTA_KT",400,-2,2,300,-1,2);
562  fhSubJettiness1CheckRatio_FJ_OP_WTA_CA= new TH2D("fhSubJettiness1CheckRatio_FJ_OP_WTA_CA","fhSubJettiness1CheckRatio_FJ_OP_WTA_CA",400,-2,2,300,-1,2);
564  fhSubJettiness1CheckRatio_FJ_MIN= new TH2D("fhSubJettiness1CheckRatio_FJ_MIN","fhSubJettiness1CheckRatio_FJ_MIN",400,-2,2,300,-1,2);
566  fhSubJettiness2CheckRatio_FJ_AKT = new TH2D("fhSubJettiness2CheckRatio_FJ_AKT","fhSubJettiness2CheckRatio_FJ_AKT",400,-2,2,300,-1,2);
568  fhSubJettiness2CheckRatio_FJ_KT= new TH2D("fhSubJettiness2CheckRatio_FJ_KT","fhSubJettiness2CheckRatio_FJ_KT",400,-2,2,300,-1,2);
570  fhSubJettiness2CheckRatio_FJ_CA= new TH2D("fhSubJettiness2CheckRatio_FJ_CA","fhSubJettiness2CheckRatio_FJ_CA",400,-2,2,300,-1,2);
572  fhSubJettiness2CheckRatio_FJ_WTA_KT= new TH2D("fhSubJettiness2CheckRatio_FJ_WTA_KT","fhSubJettiness2CheckRatio_FJ_WTA_KT",400,-2,2,300,-1,2);
574  fhSubJettiness2CheckRatio_FJ_WTA_CA= new TH2D("fhSubJettiness2CheckRatio_FJ_WTA_CA","fhSubJettiness2CheckRatio_FJ_WTA_CA",400,-2,2,300,-1,2);
576  fhSubJettiness2CheckRatio_FJ_OP_AKT= new TH2D("fhSubJettiness2CheckRatio_FJ_OP_AKT","fhSubJettiness2CheckRatio_FJ_OP_AKT",400,-2,2,300,-1,2);
578  fhSubJettiness2CheckRatio_FJ_OP_KT= new TH2D("fhSubJettiness2CheckRatio_FJ_OP_KT","fhSubJettiness2CheckRatio_FJ_OP_KT",400,-2,2,300,-1,2);
580  fhSubJettiness2CheckRatio_FJ_OP_CA= new TH2D("fhSubJettiness2CheckRatio_FJ_OP_CA","fhSubJettiness2CheckRatio_FJ_OP_CA",400,-2,2,300,-1,2);
582  fhSubJettiness2CheckRatio_FJ_OP_WTA_KT= new TH2D("fhSubJettiness2CheckRatio_FJ_OP_WTA_KT","fhSubJettiness2CheckRatio_FJ_OP_WTA_KT",400,-2,2,300,-1,2);
584  fhSubJettiness2CheckRatio_FJ_OP_WTA_CA= new TH2D("fhSubJettiness2CheckRatio_FJ_OP_WTA_CA","fhSubJettiness2CheckRatio_FJ_OP_WTA_CA",400,-2,2,300,-1,2);
586  fhSubJettiness2CheckRatio_FJ_MIN= new TH2D("fhSubJettiness2CheckRatio_FJ_MIN","fhSubJettiness2CheckRatio_FJ_MIN",400,-2,2,300,-1,2);
588  fhSubJettiness2to1CheckRatio_FJ_AKT = new TH2D("fhSubJettiness2to1CheckRatio_FJ_AKT","fhSubJettiness2to1CheckRatio_FJ_AKT",400,-2,2,300,-1,2);
590  fhSubJettiness2to1CheckRatio_FJ_KT= new TH2D("fhSubJettiness2to1CheckRatio_FJ_KT","fhSubJettiness2to1CheckRatio_FJ_KT",400,-2,2,300,-1,2);
592  fhSubJettiness2to1CheckRatio_FJ_CA= new TH2D("fhSubJettiness2to1CheckRatio_FJ_CA","fhSubJettiness2to1CheckRatio_FJ_CA",400,-2,2,300,-1,2);
594  fhSubJettiness2to1CheckRatio_FJ_WTA_KT= new TH2D("fhSubJettiness2to1CheckRatio_FJ_WTA_KT","fhSubJettiness2to1CheckRatio_FJ_WTA_KT",400,-2,2,300,-1,2);
596  fhSubJettiness2to1CheckRatio_FJ_WTA_CA= new TH2D("fhSubJettiness2to1CheckRatio_FJ_WTA_CA","fhSubJettiness2to1CheckRatio_FJ_WTA_CA",400,-2,2,300,-1,2);
598  fhSubJettiness2to1CheckRatio_FJ_OP_AKT= new TH2D("fhSubJettiness2to1CheckRatio_FJ_OP_AKT","fhSubJettiness2to1CheckRatio_FJ_OP_AKT",400,-2,2,300,-1,2);
600  fhSubJettiness2to1CheckRatio_FJ_OP_KT= new TH2D("fhSubJettiness2to1CheckRatio_FJ_OP_KT","fhSubJettiness2to1CheckRatio_FJ_OP_KT",400,-2,2,300,-1,2);
602  fhSubJettiness2to1CheckRatio_FJ_OP_CA= new TH2D("fhSubJettiness2to1CheckRatio_FJ_OP_CA","fhSubJettiness2to1CheckRatio_FJ_OP_CA",400,-2,2,300,-1,2);
604  fhSubJettiness2to1CheckRatio_FJ_OP_WTA_KT= new TH2D("fhSubJettiness2to1CheckRatio_FJ_OP_WTA_KT","fhSubJettiness2to1CheckRatio_FJ_OP_WTA_KT",400,-2,2,300,-1,2);
606  fhSubJettiness2to1CheckRatio_FJ_OP_WTA_CA= new TH2D("fhSubJettiness2to1CheckRatio_FJ_OP_WTA_CA","fhSubJettiness2to1CheckRatio_FJ_OP_WTA_CA",400,-2,2,300,-1,2);
608  fhSubJettiness2to1CheckRatio_FJ_MIN= new TH2D("fhSubJettiness2to1CheckRatio_FJ_MIN","fhSubJettiness2to1CheckRatio_FJ_MIN",400,-2,2,300,-1,2);
610  fhSubJettiness1_FJ_AKT= new TH1D("fhSubJettiness1_FJ_AKT","fhSubJettiness1_FJ_AKT",400,-2,2);
612  fhSubJettiness1_FJ_CA= new TH1D("fhSubJettiness1_FJ_CA","fhSubJettiness1_FJ_CA",400,-2,2);
614  fhSubJettiness1_FJ_WTA_KT= new TH1D("fhSubJettiness1_FJ_WTA_KT","fhSubJettiness1_FJ_WTA_KT",400,-2,2);
616  fhSubJettiness1_FJ_WTA_CA= new TH1D("fhSubJettiness1_FJ_WTA_CA","fhSubJettiness1_FJ_WTA_CA",400,-2,2);
618  fhSubJettiness1_FJ_OP_AKT= new TH1D("fhSubJettiness1_FJ_OP_AKT","fhSubJettiness1_FJ_OP_AKT",400,-2,2);
620  fhSubJettiness1_FJ_OP_KT= new TH1D("fhSubJettiness1_FJ_OP_KT","fhSubJettiness1_FJ_OP_KT",400,-2,2);
622  fhSubJettiness1_FJ_OP_CA= new TH1D("fhSubJettiness1_FJ_OP_CA","fhSubJettiness1_FJ_OP_CA",400,-2,2);
624  fhSubJettiness1_FJ_OP_WTA_KT= new TH1D("fhSubJettiness1_FJ_OP_WTA_KT","fhSubJettiness1_FJ_OP_WTA_KT",400,-2,2);
626  fhSubJettiness1_FJ_OP_WTA_CA= new TH1D("fhSubJettiness1_FJ_OP_WTA_CA","fhSubJettiness1_FJ_OP_WTA_CA",400,-2,2);
628  fhSubJettiness2_FJ_AKT= new TH1D("fhSubJettiness2_FJ_AKT","fhSubJettiness2_FJ_AKT",400,-2,2);
630  fhSubJettiness2_FJ_CA= new TH1D("fhSubJettiness2_FJ_CA","fhSubJettiness2_FJ_CA",400,-2,2);
632  fhSubJettiness2_FJ_WTA_KT= new TH1D("fhSubJettiness2_FJ_WTA_KT","fhSubJettiness2_FJ_WTA_KT",400,-2,2);
634  fhSubJettiness2_FJ_WTA_CA= new TH1D("fhSubJettiness2_FJ_WTA_CA","fhSubJettiness2_FJ_WTA_CA",400,-2,2);
636  fhSubJettiness2_FJ_OP_AKT= new TH1D("fhSubJettiness2_FJ_OP_AKT","fhSubJettiness2_FJ_OP_AKT",400,-2,2);
638  fhSubJettiness2_FJ_OP_KT= new TH1D("fhSubJettiness2_FJ_OP_KT","fhSubJettiness2_FJ_OP_KT",400,-2,2);
640  fhSubJettiness2_FJ_OP_CA= new TH1D("fhSubJettiness2_FJ_OP_CA","fhSubJettiness2_FJ_OP_CA",400,-2,2);
642  fhSubJettiness2_FJ_OP_WTA_KT= new TH1D("fhSubJettiness2_FJ_OP_WTA_KT","fhSubJettiness2_FJ_OP_WTA_KT",400,-2,2);
644  fhSubJettiness2_FJ_OP_WTA_CA= new TH1D("fhSubJettiness2_FJ_OP_WTA_CA","fhSubJettiness2_FJ_OP_WTA_CA",400,-2,2);
646  fhSubJettiness2to1_FJ_KT= new TH1D("fhSubJettiness2to1_FJ_KT","fhSubJettiness2to1_FJ_KT",400,-2,2);
648  fhSubJettiness2to1_FJ_AKT= new TH1D("fhSubJettiness2to1_FJ_AKT","fhSubJettiness2to1_FJ_AKT",400,-2,2);
650  fhSubJettiness2to1_FJ_CA= new TH1D("fhSubJettiness2to1_FJ_CA","fhSubJettiness2to1_FJ_CA",400,-2,2);
652  fhSubJettiness2to1_FJ_WTA_KT= new TH1D("fhSubJettiness2to1_FJ_WTA_KT","fhSubJettiness2to1_FJ_WTA_KT",400,-2,2);
654  fhSubJettiness2to1_FJ_WTA_CA= new TH1D("fhSubJettiness2to1_FJ_WTA_CA","fhSubJettiness2to1_FJ_WTA_CA",400,-2,2);
656  fhSubJettiness2to1_FJ_OP_AKT= new TH1D("fhSubJettiness2to1_FJ_OP_AKT","fhSubJettiness2to1_FJ_OP_AKT",400,-2,2);
658  fhSubJettiness2to1_FJ_OP_KT= new TH1D("fhSubJettiness2to1_FJ_OP_KT","fhSubJettiness2to1_FJ_OP_KT",400,-2,2);
660  fhSubJettiness2to1_FJ_OP_CA= new TH1D("fhSubJettiness2to1_FJ_OP_CA","fhSubJettiness2to1_FJ_OP_CA",400,-2,2);
662  fhSubJettiness2to1_FJ_OP_WTA_KT= new TH1D("fhSubJettiness2to1_FJ_OP_WTA_KT","fhSubJettiness2to1_FJ_OP_WTA_KT",400,-2,2);
664  fhSubJettiness2to1_FJ_OP_WTA_CA= new TH1D("fhSubJettiness2to1_FJ_OP_WTA_CA","fhSubJettiness2to1_FJ_OP_WTA_CA",400,-2,2);
666  fhSubJettiness2to1_FJ_MIN= new TH1D("fhSubJettiness2to1_FJ_MIN","fhSubJettiness2to1_FJ_MIN",400,-2,2);
668  }
670  fhJetPt_1= new TH1F("fhJetPt_1", "Jet Pt Detector Level",1500,-0.5,149.5 );
671  fOutput->Add(fhJetPt_1);
672  fhJetPt_2= new TH1F("fhJetPt_2", "Jet Pt Particle Level",1500,-0.5,149.5 );
673  fOutput->Add(fhJetPt_2);
674  fhJetPhi_1= new TH1F("fhJetPhi_1", "Jet Phi Detector Level",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
675  fOutput->Add(fhJetPhi_1);
676  fhJetPhi_2= new TH1F("fhJetPhi_2", "Jet Phi Particle Level",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
677  fOutput->Add(fhJetPhi_2);
678  fhJetEta_1= new TH1F("fhJetEta_1", "Jet Eta Detector Level", Eta_Bins, Eta_Low, Eta_Up);
679  fOutput->Add(fhJetEta_1);
680  fhJetEta_2= new TH1F("fhJetEta_2", "Jet Eta Particle Level", Eta_Bins, Eta_Low, Eta_Up);
681  fOutput->Add(fhJetEta_2);
682  fhJetMass_1= new TH1F("fhJetMass_1", "Jet Mass Detector Level", 4000,-0.5, 39.5);
683  fOutput->Add(fhJetMass_1);
684  fhJetMass_2= new TH1F("fhJetMass_2", "Jet Mass Particle Level", 4000,-0.5, 39.5);
685  fOutput->Add(fhJetMass_2);
686  fhJetRadius_1= new TH1F("fhJetRadius_1", "Jet Radius Detector Level", 100, -0.05,0.995);
687  fOutput->Add(fhJetRadius_1);
688  fhJetRadius_2= new TH1F("fhJetRadius_2", "Jet Radius Particle Level", 100, -0.05,0.995);
689  fOutput->Add(fhJetRadius_2);
690  fhNumberOfJetTracks_1= new TH1F("fhNumberOfJetTracks_1", "Number of Tracks within a Jet Detector Level", 300, -0.5,299.5);
692  fhNumberOfJetTracks_2= new TH1F("fhNumberOfJetTracks_2", "Number of Tracks within a Jet Particle Level", 300, -0.5,299.5);
694  fhSubJetRadius_1= new TH1F("fhSubJetRadius_1", "SubJet Radius Detector Level", 100, -0.05,0.995);
696  fhSubJetRadius_2= new TH1F("fhSubJetRadius_2", "SubJet Radius Particle Level", 100, -0.05,0.995);
698  fhSubJetPt_1= new TH1F("fhSubJetPt_1", "SubJet Pt Detector Level", 1500, -0.5,149.5);
699  fOutput->Add(fhSubJetPt_1);
700  fhSubJetPt_2= new TH1F("fhSubJetPt_2", "SubJet Pt Particle Level", 1500, -0.5,149.5);
701  fOutput->Add(fhSubJetPt_2);
702  fhSubJetMass_1= new TH1F("fhSubJetMass_1", "Sub Jet Mass Detector Level", 4000,-0.5, 39.5);
703  fOutput->Add(fhSubJetMass_1);
704  fhSubJetMass_2= new TH1F("fhSubJetMass_2", "Sub Jet Mass Particle Level", 4000,-0.5, 39.5);
705  fOutput->Add(fhSubJetMass_2);
706  fhNumberOfSubJetTracks_1= new TH1F("fhNumberOfSubJetTracks_1", "Number of Tracks within a Sub Jet Detector Level", 300, -0.5,299.5);
708  fhNumberOfSubJetTracks_2= new TH1F("fhNumberOfSubJetTracks_2", "Number of Tracks within a Sub Jet Particle Level", 300, -0.5,299.5);
710  fhJetCounter_1= new TH1F("fhJetCounter_1", "Jet Counter Detector Level", 150, -0.5, 149.5);
711  fOutput->Add(fhJetCounter_1);
712  fhJetCounter_2= new TH1F("fhJetCounter_2", "Jet Counter Particle Level", 150, -0.5, 149.5);
713  fOutput->Add(fhJetCounter_2);
714  fhSubJetCounter_1 = new TH1F("fhSubJetCounter_1", "SubJet Counter Detector Level",50, -0.5,49.5);
716  fhSubJetCounter_2 = new TH1F("fhSubJetCounter_2", "SubJet Counter Particle Level",50, -0.5,49.5);
718  fh2PtRatio= new TH2F("fhPtRatio", "MC pt for detector level vs particle level jets",1500,-0.5,149.5,1500,-0.5,149.5);
719  fOutput->Add(fh2PtRatio);
720  fhEventCounter_1= new TH1F("fhEventCounter_1", "Event Counter Detector Level", 15,0.5,15.5);
722  fhEventCounter_2= new TH1F("fhEventCounter_2", "Event Counter Particle Level", 15,0.5,15.5);
724  }
726  fhEventCounter= new TH1F("fhEventCounter", "Event Counter", 15,0.5,15.5);
727  fOutput->Add(fhEventCounter);
728  }
729 
730  PostData(1,fOutput);
731  PostData(2,fTreeResponseMatrixAxis);
732  // delete [] fShapesVarNames;
733 }
734 
735 //________________________________________________________________________
737 {
738  // Run analysis code here, if needed. It will be executed before FillHistograms().
739 
740 
741  return kTRUE;
742 }
743 
744 //________________________________________________________________________
746 {
747 
748  //fhEventCounter Key:
749  // 1: Number of events with a Jet Container
750  // 2: Number of Jets without a Jet Container
751  // 3:
752  // 4: Number of Jets found in all events
753  // 5: Number of Jets that were reclustered in all events
754  // 6: Number of SubJets found in all events
755  // 7: Number of events were primary vertext was found
756  // 8: Number of Jets with more than one SubJet
757  // 9:Number of Jets with more than two SubJets
758  // 12:Number of SubJetinessEvents in kData
759 
760 
761  // const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
762  // Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
763  // if(vert) fhEventCounter->Fill(7);
764 
765  //cout << ((AliVAODHeader*)(InputEvent()->GetHeader()))->GetEventplane()<< " "<<fEPV0<<endl;
766  // cout << InputEvent()->GetEventplane()->GetEventplane("Q")<< " "<<fEPV0<<endl;
767  if (fCentSelectOn){
768  if ((fCent>fCentMax) || (fCent<fCentMin)) return 0;
769  }
770 
771  AliAODTrack *TriggerHadron = 0x0;
772  if (fJetSelection == kRecoil) {
773  //you have to set a flag and the limits of the pT interval for your trigger
775  if (TriggerHadronLabel==-99999) return 0; //Trigger Hadron Not Found
776  AliTrackContainer *PartCont =NULL;
777  AliParticleContainer *PartContMC=NULL;
778  if (fJetShapeSub==kConstSub){
780  else PartCont = GetTrackContainer(1);
781  }
782  else{
784  else PartCont = GetTrackContainer(0);
785  }
786  TClonesArray *TrackArray = NULL;
787  TClonesArray *TrackArrayMC = NULL;
788  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) TrackArrayMC = PartContMC->GetArray();
789  else TrackArray = PartCont->GetArray();
790  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) TriggerHadron = static_cast<AliAODTrack*>(TrackArrayMC->At(TriggerHadronLabel));
791  else TriggerHadron = static_cast<AliAODTrack*>(TrackArray->At(TriggerHadronLabel));
792  if (!TriggerHadron) return 0;//No trigger hadron with label found
793  if(fSemigoodCorrect){
794  Double_t HoleDistance=RelativePhi(TriggerHadron->Phi(),fHolePos);
795  if(TMath::Abs(HoleDistance)+fHoleWidth+fJetRadius>TMath::Pi()-fRecoilAngularWindow) return 0;
796  }
797  fhPtTriggerHadron->Fill(TriggerHadron->Pt()); //Needed for per trigger Normalisation
798  if (fJetShapeType != AliAnalysisTaskSubJetFraction::kGenOnTheFly) fhPhiTriggerHadronEventPlane->Fill(TMath::Abs(RelativePhiEventPlane(fEPV0,TriggerHadron->Phi()))); //fEPV0 is the event plane from AliAnalysisTaskEmcal
799  if (fJetShapeType != AliAnalysisTaskSubJetFraction::kGenOnTheFly) fhPhiTriggerHadronEventPlaneTPC->Fill(TMath::Abs(RelativePhiEventPlane(((AliVAODHeader*)(InputEvent()->GetHeader()))->GetEventplane(),TriggerHadron->Phi()))); //TPC event plane
800  }
801 
802 
805  AliEmcalJet *Jet1 = NULL; //Subtracted hybrid Jet
806  AliEmcalJet *Jet2 = NULL; //Unsubtracted Hybrid Jet
807  AliEmcalJet *Jet3 = NULL; //Detector Level Pythia Jet
808  AliEmcalJet *Jet4 = NULL; //Particle Level Pyhtia Jet
809  AliJetContainer *JetCont1= GetJetContainer(0); //Jet Container for Subtracted Hybrid Jets
810  AliJetContainer *JetCont2= GetJetContainer(1); //Jet Container for Unsubtracted Hybrid Jets
811  AliJetContainer *JetCont3= GetJetContainer(2); //Jet Container for Detector Level Pyhtia Jets
812  AliJetContainer *JetCont4= GetJetContainer(3); //Jet Container for Particle Level Pythia Jets
813  AliEmcalJetFinder *Reclusterer1; //Object containg Subjets from Subtracted Hybrid Jets
814  AliEmcalJetFinder *Reclusterer4; //Object containg Subjets from Particle Level Pythia Jets
815  Bool_t JetsMatched=kFALSE;
816  Bool_t EventCounter=kFALSE;
817  Int_t JetNumber=-1;
818  Double_t JetPtThreshold=-2;
819  fhEventCounter->Fill(1);
820  if(JetCont1) {
821  fhEventCounter->Fill(2);
822  JetCont1->ResetCurrentID();
823  while((Jet1=JetCont1->GetNextAcceptJet())) {
824  if (fJetShapeSub==kConstSub) JetPtThreshold=Jet1->Pt();
825  else JetPtThreshold=Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
826  if ( (!Jet1) || (JetPtThreshold<fPtThreshold)) continue;
827  else {
828  /* if(fSemigoodCorrect){
829  Double_t HoleDistance=RelativePhi(Jet1->Phi(),fHolePos);
830  if(TMath::Abs(HoleDistance)<fHoleWidth) continue;
831  }*/
832  Float_t RecoilDeltaPhi = 0.;
833  if (fJetSelection == kRecoil){
834  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), Jet1->Phi());
835  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
836  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), Jet1->Pt());
837  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), Jet1->Phi()));
838  }
839  if (!EventCounter){
840  fhEventCounter->Fill(3);
841  EventCounter=kTRUE;
842  }
843  if(fJetShapeSub==kConstSub){
844  JetNumber=-1;
845  for(Int_t i = 0; i<JetCont2->GetNJets(); i++) {
846  Jet2 = JetCont2->GetJet(i);
847  if(Jet2->GetLabel()==Jet1->GetLabel()) {
848  JetNumber=i;
849  }
850  }
851  if(JetNumber==-1) continue;
852  Jet2=JetCont2->GetJet(JetNumber);
853  if (JetCont2->AliJetContainer::GetFractionSharedPt(Jet2)<fSharedFractionPtMin) continue;
854  Jet3=Jet2->ClosestJet();
855  }
856  if(!(fJetShapeSub==kConstSub)){
857  if (JetCont1->AliJetContainer::GetFractionSharedPt(Jet1)<fSharedFractionPtMin) continue;
858  Jet3 = Jet1->ClosestJet(); //Note for NoSub and Deriv Sub cases you must fill both the Unsubtracted and Subtracted Hybrid jet containers with the same jet branch
859  }
860  if (!Jet3) continue;
861  Jet4=Jet3->ClosestJet();
862  if(!Jet4) continue;
863  JetsMatched=kTRUE;
865  if (fJetShapeSub==kConstSub) fShapesVar[0]=Jet1->Pt();
866  else fShapesVar[0]=Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
867  fShapesVar[2]=FjNSubJettiness(Jet1,0,1,0,1,0);
868  fShapesVar[4]=FjNSubJettiness(Jet1,0,2,0,1,0);
869  fShapesVar[6]=FjNSubJettiness(Jet1,0,2,0,1,8);
870  fShapesVar[8]=FjNSubJettiness(Jet1,0,2,0,1,1);
871  // fShapesVar[10]=Jet1->GetNumberOfTracks();
872  fShapesVar[10]=FjNSubJettiness(Jet1,0,2,0,1,9);
873  fShapesVar[12]=FjNSubJettiness(Jet1,0,2,0,1,3,fBeta_SD,fZCut);
874  fShapesVar[14]=FjNSubJettiness(Jet1,0,2,0,1,5,fBeta_SD,fZCut);
875  fShapesVar[16]=Jet1->GetLeadingTrack(JetCont1->GetParticleContainer()->GetArray())->Pt();
877  //fShapesVar[20]=RelativePhiEventPlane(((AliVAODHeader*)(InputEvent()->GetHeader()))->GetEventplane(),Jet1->Phi());
878  fShapesVar[20]=FjNSubJettiness(Jet1,0,2,0,1,6,fBeta_SD,fZCut);
879  if (fFullTree){
880  fShapesVar[22]=FjNSubJettiness(Jet1,0,2,0,1,2);
881  Reclusterer1 = Recluster(Jet1, 0, fSubJetRadius, fSubJetMinPt, fSubJetAlgorithm, "SubJetFinder_1");
882  fShapesVar[24]=SubJetFraction(Jet1, Reclusterer1, 1, 0, kTRUE, kFALSE);
883  fShapesVar[26]=SubJetFraction(Jet1, Reclusterer1, 2, 0, kTRUE, kFALSE);
884  }
885  if (JetsMatched){ //even needed? Not now but might be if you want to fill trees when jets aren't matched too
886  fShapesVar[1]=Jet4->Pt();
887  fShapesVar[3]=FjNSubJettiness(Jet4,3,1,0,1,0);
888  fShapesVar[5]=FjNSubJettiness(Jet4,3,2,0,1,0);
889  fShapesVar[7]=FjNSubJettiness(Jet4,3,2,0,1,8);
890  fShapesVar[9]=FjNSubJettiness(Jet4,3,2,0,1,1);
891  //fShapesVar[11]=Jet4->GetNumberOfTracks();
892  fShapesVar[11]=FjNSubJettiness(Jet4,3,2,0,1,9);
893  fShapesVar[13]=FjNSubJettiness(Jet4,3,2,0,1,3,fBeta_SD,fZCut);
894  fShapesVar[15]=FjNSubJettiness(Jet4,3,2,0,1,5,fBeta_SD,fZCut);
895  fShapesVar[17]=Jet4->GetLeadingTrack(JetCont4->GetParticleContainer()->GetArray())->Pt();
896  fShapesVar[19]=RelativePhiEventPlane(fEPV0,Jet4->Phi());
897  //fShapesVar[21]=RelativePhiEventPlane(((AliVAODHeader*)(InputEvent()->GetHeader()))->GetEventplane(),Jet4->Phi());
898  fShapesVar[21]=FjNSubJettiness(Jet4,3,2,0,1,6,fBeta_SD,fZCut);
899  if (fFullTree){
900  fShapesVar[23]=FjNSubJettiness(Jet4,3,2,0,1,2);
901  Reclusterer4=Recluster(Jet4, 3, fSubJetRadius, 0, fSubJetAlgorithm, "SubJetFinder_4");
902  fShapesVar[25]=SubJetFraction(Jet4, Reclusterer4, 1, 0, kTRUE, kFALSE);
903  fShapesVar[27]=SubJetFraction(Jet4, Reclusterer4, 2, 0, kTRUE, kFALSE);
904  }
905  }
906  else{
907  fShapesVar[1]=-2;
908  fShapesVar[3]=-2;
909  fShapesVar[5]=-2;
910  fShapesVar[7]=-2;
911  fShapesVar[9]=-2;
912  fShapesVar[11]=-2;
913  fShapesVar[13]=-2;
914  fShapesVar[15]=-2;
915  fShapesVar[17]=-2;
916  fShapesVar[19]=-2;
917  fShapesVar[21]=-2;
918  if (fFullTree){
919  fShapesVar[23]=-2;
920  fShapesVar[25]=-2;
921  fShapesVar[27]=-2;
922  }
923  }
924  fTreeResponseMatrixAxis->Fill();
925  JetsMatched=kFALSE;
926  }
927  }
928  }
929  }
930 
933  AliEmcalJet *Jet1 = NULL; //Detector Level Jet
934  AliEmcalJet *Jet2 = NULL; //Particle Level Jet
935  AliJetContainer *JetCont1= GetJetContainer(0); //Jet Container for Detector Level Pythia
936  AliJetContainer *JetCont2= GetJetContainer(1); //Jet Container for Particle Level Pythia
937  AliEmcalJetFinder *Reclusterer1; //Object containg Subjets Detector Level
938  AliEmcalJetFinder *Reclusterer2; //Object containg Subjets Particle Level
939  Int_t JetCounter1=0; //Counts number of jets in event
940  Int_t JetCounter2=0; //Counts number of jets in event
941  Double_t JetPhi1=0;
942  Double_t JetPhi2=0;
943  Bool_t JetsMatched=kFALSE;
944  Double_t Pythia_Event_Weight=1;
945  Bool_t EventCounter=kFALSE;
946  fhEventCounter_1->Fill(1);
947  if(JetCont1) {
948  fhEventCounter_1->Fill(2); //Number of events with a jet container
949  JetCont1->ResetCurrentID();
950  while((Jet1=JetCont1->GetNextAcceptJet())) {
951  if( (!Jet1) || ((Jet1->Pt())<fPtThreshold)) {
952  // fhEventCounter_1->Fill(3); //events where the jet had a problem
953  continue;
954  }
955  else {
956  /* if(fSemigoodCorrect){
957  Double_t HoleDistance=RelativePhi(Jet1->Phi(),fHolePos);
958  if(TMath::Abs(HoleDistance)<fHoleWidth) continue;
959  }*/
960  Float_t RecoilDeltaPhi = 0.;
961  if (fJetSelection == kRecoil){
962  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), Jet1->Phi());
963  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
964  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), Jet1->Pt());
965  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), Jet1->Phi()));
966  }
967  if (!EventCounter){
968  fhEventCounter_1->Fill(3);
969  EventCounter=kTRUE;
970  }
971  if((Jet1->GetNumberOfTracks())==0){
972  fhEventCounter_1->Fill(10); //zero track jets
973  }
974  if((Jet1->GetNumberOfTracks())==1){
975  fhEventCounter_1->Fill(11); //one track jets
976  }
977  fhEventCounter_1->Fill(4); //Number of Jets found in all events
978  JetCounter1++;
979  fhJetPt_1->Fill(Jet1->Pt());
980  JetPhi1=Jet1->Phi();
981  if(JetPhi1 < -1*TMath::Pi()) JetPhi1 += (2*TMath::Pi());
982  else if (JetPhi1 > TMath::Pi()) JetPhi1 -= (2*TMath::Pi());
983  fhJetPhi_1->Fill(JetPhi1);
984  fhJetEta_1->Fill(Jet1->Eta());
985  fhJetMass_1->Fill(Jet1->M());
986  fhJetRadius_1->Fill(TMath::Sqrt((Jet1->Area()/TMath::Pi()))); //Radius of Jets per event
988  if((Jet2 = Jet1->ClosestJet())){
989  JetsMatched=kTRUE;
990  if((Jet2->GetNumberOfTracks())==0){
991  fhEventCounter_2->Fill(10); //zero track jets
992  }
993  if((Jet2->GetNumberOfTracks())==1){
994  fhEventCounter_2->Fill(11); //one track jets
995  }
996  fhEventCounter_2->Fill(4); //Number of Jets found in all events
997  JetCounter2++;
998  fhJetPt_2->Fill(Jet2->Pt());
999  JetPhi2=Jet2->Phi();
1000  if(JetPhi2 < -1*TMath::Pi()) JetPhi2 += (2*TMath::Pi());
1001  else if (JetPhi2 > TMath::Pi()) JetPhi2 -= (2*TMath::Pi());
1002  fhJetPhi_2->Fill(JetPhi2);
1003  fhJetEta_2->Fill(Jet2->Eta());
1004  fhJetMass_2->Fill(Jet2->M());
1005  fhJetRadius_2->Fill(TMath::Sqrt((Jet2->Area()/TMath::Pi()))); //Radius of Jets per event
1007  fh2PtRatio->Fill(Jet1->Pt(),Jet2->Pt(),Pythia_Event_Weight);
1008  }
1009  else {
1010  fhEventCounter_2->Fill(1);
1011  //continue;
1012  }
1014 
1015 
1016  fShapesVar[0]=Jet1->Pt();
1017  fShapesVar[2]=FjNSubJettiness(Jet1,0,1,0,1,0);
1018  fShapesVar[4]=FjNSubJettiness(Jet1,0,2,0,1,0);
1019  fShapesVar[6]=FjNSubJettiness(Jet1,0,2,0,1,8);
1020  fShapesVar[8]=FjNSubJettiness(Jet1,0,2,0,1,1);
1021  // fShapesVar[10]=Jet1->GetNumberOfTracks();
1022  fShapesVar[10]=FjNSubJettiness(Jet1,0,2,0,1,9);
1023  fShapesVar[12]=FjNSubJettiness(Jet1,0,2,0,1,3,fBeta_SD,fZCut);
1024  fShapesVar[14]=FjNSubJettiness(Jet1,0,2,0,1,5,fBeta_SD,fZCut);
1025  fShapesVar[16]=Jet1->GetLeadingTrack(JetCont1->GetParticleContainer()->GetArray())->Pt();
1026  fShapesVar[18]=-2; //event plane calculation only needed for PbPb recoils
1027  fShapesVar[20]=FjNSubJettiness(Jet1,0,2,0,1,6,fBeta_SD,fZCut);
1028  Reclusterer1 = Recluster(Jet1, 0, fSubJetRadius, fSubJetMinPt, fSubJetAlgorithm, "SubJetFinder_1");
1029  if (fFullTree){
1030  fShapesVar[22]=FjNSubJettiness(Jet1,0,2,0,1,2);
1031  fShapesVar[24]=SubJetFraction(Jet1, Reclusterer1, 1, 0, kTRUE, kFALSE);
1032  fShapesVar[26]=SubJetFraction(Jet1, Reclusterer1, 2, 0, kTRUE, kFALSE);
1033  }
1034  if (JetsMatched){ //even needed? Not now but might be if you want to fill trees when jets aren't matched too
1035  fShapesVar[1]=Jet2->Pt();
1036  fShapesVar[3]=FjNSubJettiness(Jet2,1,1,0,1,0);
1037  fShapesVar[5]=FjNSubJettiness(Jet2,1,2,0,1,0);
1038  fShapesVar[7]=FjNSubJettiness(Jet2,1,2,0,1,8);
1039  fShapesVar[9]=FjNSubJettiness(Jet2,1,2,0,1,1);
1040  // fShapesVar[11]=Jet2->GetNumberOfTracks();
1041  fShapesVar[11]=FjNSubJettiness(Jet2,1,2,0,1,9);
1042  fShapesVar[13]=FjNSubJettiness(Jet2,1,2,0,1,3,fBeta_SD,fZCut);
1043  fShapesVar[15]=FjNSubJettiness(Jet2,1,2,0,1,5,fBeta_SD,fZCut);
1044  fShapesVar[17]=Jet2->GetLeadingTrack(JetCont2->GetParticleContainer()->GetArray())->Pt();
1045  fShapesVar[19]=-2;
1046  fShapesVar[21]=FjNSubJettiness(Jet2,1,2,0,1,6,fBeta_SD,fZCut);
1047  Reclusterer2 = Recluster(Jet2, 1, fSubJetRadius, 0, fSubJetAlgorithm, "SubJetFinder_2");
1048  if (fFullTree){
1049  fShapesVar[23]=FjNSubJettiness(Jet2,1,2,0,1,2);
1050  fShapesVar[25]=SubJetFraction(Jet2, Reclusterer2, 1, 0, kTRUE, kFALSE);
1051  fShapesVar[27]=SubJetFraction(Jet2, Reclusterer2, 2, 0, kTRUE, kFALSE);
1052  }
1053  }
1054  else{
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  fShapesVar[13]=-2;
1062  fShapesVar[15]=-2;
1063  fShapesVar[17]=-2;
1064  fShapesVar[19]=-2;
1065  fShapesVar[21]=-2;
1066  if (fFullTree){
1067  fShapesVar[23]=-2;
1068  fShapesVar[25]=-2;
1069  fShapesVar[27]=-2;
1070  }
1071  }
1072  fTreeResponseMatrixAxis->Fill();
1073 
1074  fhSubJetCounter_1->Fill(Reclusterer1->GetNumberOfJets());
1075  for (Int_t i= 0; i<Reclusterer1->GetNumberOfJets(); i++){
1076  fhEventCounter_1->Fill(6); //Number of overall subjets in all events
1077  fhSubJetPt_1->Fill(Reclusterer1->GetJet(i)->Pt());
1078  fhSubJetMass_1->Fill(Reclusterer1->GetJet(i)->M());
1079  fhNumberOfSubJetTracks_1->Fill(Reclusterer1->GetJet(i)->GetNumberOfTracks());
1080  fhSubJetRadius_1->Fill(TMath::Sqrt((Reclusterer1->GetJet(i)->Area()/TMath::Pi()))); //Radius of SubJets per event
1081  }
1082  if(JetsMatched){
1083  fhSubJetCounter_2->Fill(Reclusterer2->GetNumberOfJets());
1084  for (Int_t i= 0; i<Reclusterer2->GetNumberOfJets(); i++){
1085  fhEventCounter_2->Fill(6); //Number of overall subjets in all events
1086  fhSubJetPt_2->Fill(Reclusterer2->GetJet(i)->Pt());
1087  fhSubJetMass_2->Fill(Reclusterer2->GetJet(i)->M());
1088  fhNumberOfSubJetTracks_2->Fill(Reclusterer2->GetJet(i)->GetNumberOfTracks());
1089  fhSubJetRadius_2->Fill(TMath::Sqrt((Reclusterer2->GetJet(i)->Area()/TMath::Pi()))); //Radius of SubJets per event
1090  }
1091  }
1092  JetsMatched=kFALSE;
1093  }
1094  }
1095  fhJetCounter_1->Fill(JetCounter1); //Number of Jets in Each Event Particle Level
1096  fhJetCounter_2->Fill(JetCounter2); //Number of Jets in Each Event Detector Level
1097  }
1098  //else {fhEventCounter_1->Fill(3);} //Events with no jet container
1099  }
1100 
1101 
1102 
1104  AliEmcalJet *Jet1 = NULL; //Original Jet in the event
1105  AliJetContainer *JetCont= GetJetContainer(0); //Jet Container for event
1106  Int_t JetCounter=0; //Counts number of jets in event
1107  Double_t JetPhi=0;
1108  Bool_t EventCounter=kFALSE;
1109  Double_t JetPt_ForThreshold=0;
1110  fhEventCounter->Fill(1);
1111  if(JetCont) {
1112  fhEventCounter->Fill(2); //Number of events with a jet container
1113  JetCont->ResetCurrentID();
1114  while((Jet1=JetCont->GetNextAcceptJet())) {
1115  if(!Jet1) continue;
1116  if (fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) JetPt_ForThreshold = Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
1117  else JetPt_ForThreshold = Jet1->Pt();
1118  if(JetPt_ForThreshold<fPtThreshold) {
1119  //fhEventCounter->Fill(3); //events where the jet had a problem
1120  continue;
1121  }
1122  else {
1123  /* if(fSemigoodCorrect){
1124  Double_t HoleDistance=RelativePhi(Jet1->Phi(),fHolePos);
1125  if(TMath::Abs(HoleDistance)<fHoleWidth) continue;
1126  }*/
1127  Float_t RecoilDeltaPhi = 0.;
1128  if (fJetSelection == kRecoil){
1129  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), Jet1->Phi());
1130  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
1131  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), Jet1->Pt());
1132  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), Jet1->Phi()));
1133  }
1134  if (!EventCounter){
1135  fhEventCounter->Fill(3);
1136  EventCounter=kTRUE;
1137  }
1138  if((Jet1->GetNumberOfTracks())==0){
1139  fhEventCounter->Fill(10); //zero track jets
1140  }
1141  if((Jet1->GetNumberOfTracks())==1){
1142  fhEventCounter->Fill(11); //one track jets
1143  }
1144  fhEventCounter->Fill(4); //Number of Jets found in all events
1145  JetCounter++;
1146  fhJetPt->Fill(Jet1->Pt());
1147  JetPhi=Jet1->Phi();
1148  if(JetPhi < -1*TMath::Pi()) JetPhi += (2*TMath::Pi());
1149  else if (JetPhi > TMath::Pi()) JetPhi -= (2*TMath::Pi());
1150  fhJetPhi->Fill(JetPhi);
1151  fhJetEta->Fill(Jet1->Eta());
1152  fhJetMass->Fill(Jet1->M());
1153  fhJetRadius->Fill(TMath::Sqrt((Jet1->Area()/TMath::Pi()))); //Radius of Jets per event
1154  fhNumberOfJetTracks->Fill(Jet1->GetNumberOfTracks());
1155  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) fShapesVar[0]= Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
1156  else fShapesVar[0]=Jet1->Pt();
1157  fShapesVar[2]=FjNSubJettiness(Jet1,0,1,0,1,0);
1158  fShapesVar[4]=FjNSubJettiness(Jet1,0,2,0,1,0);
1159  fShapesVar[6]=FjNSubJettiness(Jet1,0,2,0,1,8);
1160  fShapesVar[8]=FjNSubJettiness(Jet1,0,2,0,1,1);
1161  //fShapesVar[10]=Jet1->GetNumberOfTracks();
1162  fShapesVar[10]=FjNSubJettiness(Jet1,0,2,0,1,9);
1163  fShapesVar[12]=FjNSubJettiness(Jet1,0,2,0,1,3,fBeta_SD,fZCut);
1164  fShapesVar[14]=FjNSubJettiness(Jet1,0,2,0,1,5,fBeta_SD,fZCut);
1165  fShapesVar[16]=Jet1->GetLeadingTrack(JetCont->GetParticleContainer()->GetArray())->Pt();
1166  fShapesVar[18]=-2; //event plane calculation not needed for data
1167  fShapesVar[20]=FjNSubJettiness(Jet1,0,2,0,1,6,fBeta_SD,fZCut);
1168  AliEmcalJetFinder *Reclusterer1 = Recluster(Jet1, 0, fSubJetRadius, fSubJetMinPt, fSubJetAlgorithm, "SubJetFinder");
1169  if (fFullTree){
1170  fShapesVar[22]=FjNSubJettiness(Jet1,0,2,0,1,2);
1171  fShapesVar[24]=SubJetFraction(Jet1, Reclusterer1, 1, 0, kTRUE, kFALSE);
1172  fShapesVar[26]=SubJetFraction(Jet1, Reclusterer1, 2, 0, kTRUE, kFALSE);
1173  }
1174  fShapesVar[1]=-2;
1175  fShapesVar[3]=-2;
1176  fShapesVar[5]=-2;
1177  fShapesVar[7]=-2;
1178  fShapesVar[9]=-2;
1179  fShapesVar[11]=-2;
1180  fShapesVar[13]=-2;
1181  fShapesVar[15]=-2;
1182  fShapesVar[17]=-2;
1183  fShapesVar[19]=-2;
1184  fShapesVar[21]=-2;
1185  if (fFullTree){
1186  fShapesVar[23]=-2;
1187  fShapesVar[25]=-2;
1188  fShapesVar[27]=-2;
1189  }
1190  fTreeResponseMatrixAxis->Fill();
1191  fhSubJetCounter->Fill(Reclusterer1->GetNumberOfJets());
1192  for (Int_t i= 0; i<Reclusterer1->GetNumberOfJets(); i++){
1193  fhEventCounter->Fill(6); //Number of overall subjets in all events
1194  fhSubJetPt->Fill(Reclusterer1->GetJet(i)->Pt());
1195  fhSubJetMass->Fill(Reclusterer1->GetJet(i)->M());
1196  fhNumberOfSubJetTracks->Fill(Reclusterer1->GetJet(i)->GetNumberOfTracks());
1197  fhSubJetRadius->Fill(TMath::Sqrt((Reclusterer1->GetJet(i)->Area()/TMath::Pi()))); //Radius of SubJets per event
1198  }
1199  }
1200  }
1201  fhJetCounter->Fill(JetCounter); //Number of Jets in Each Event
1202  }
1203  //else {fhEventCounter->Fill(2);} //Events with no jet container
1204  }
1205 
1206 
1208 
1209  AliEmcalJet *Jet1 = NULL; //Original Jet in the event
1210  AliJetContainer *JetCont= GetJetContainer(0); //Jet Container for event
1211  Int_t JetCounter=0; //Counts number of jets in event
1212  Double_t JetPhi=0;
1213  Bool_t EventCounter=kFALSE;
1214  Double_t JetPt_ForThreshold=0;
1215  fhEventCounter->Fill(1);
1216  if(JetCont) {
1217  fhEventCounter->Fill(2); //Number of events with a jet container
1218  JetCont->ResetCurrentID();
1219  while((Jet1=JetCont->GetNextAcceptJet())) {
1220  // if (((Jet1=JetCont->GetLeadingJet("rho")) && (fPtThreshold<=Jet1->Pt()-(GetRhoVal(0)*Jet1->Area())))){
1221  if(!Jet1) continue;
1222  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) JetPt_ForThreshold = Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
1223  else JetPt_ForThreshold = Jet1->Pt();
1224  if(JetPt_ForThreshold<fPtThreshold) {
1225  //fhEventCounter->Fill(3); //events where the jet had a problem
1226  continue;
1227  }
1228  else {
1229  /* if(fSemigoodCorrect){
1230  Double_t HoleDistance=RelativePhi(Jet1->Phi(),fHolePos);
1231  if(TMath::Abs(HoleDistance)<fHoleWidth) continue;
1232  }*/
1233  Float_t RecoilDeltaPhi = 0.;
1234  if (fJetSelection == kRecoil){
1235  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), Jet1->Phi());
1236  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
1237  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), Jet1->Pt());
1238  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), Jet1->Phi()));
1239  }
1240  if (!EventCounter){
1241  fhEventCounter->Fill(3);
1242  EventCounter=kTRUE;
1243  }
1244  if((Jet1->GetNumberOfTracks())==0){
1245  fhEventCounter->Fill(10); //zero track jets
1246  }
1247  if((Jet1->GetNumberOfTracks())==1){
1248  fhEventCounter->Fill(11); //one track jets
1249  }
1250  fhEventCounter->Fill(4); //Number of Jets found in all events
1251  JetCounter++;
1252  fhJetPt->Fill(Jet1->Pt());
1253  JetPhi=Jet1->Phi();
1254  if(JetPhi < -1*TMath::Pi()) JetPhi += (2*TMath::Pi());
1255  else if (JetPhi > TMath::Pi()) JetPhi -= (2*TMath::Pi());
1256  fhJetPhi->Fill(JetPhi);
1257  fhJetEta->Fill(Jet1->Eta());
1258  fhJetMass->Fill(Jet1->M());
1259  fhJetRadius->Fill(TMath::Sqrt((Jet1->Area()/TMath::Pi()))); //Radius of Jets per event
1260  fhNumberOfJetTracks->Fill(Jet1->GetNumberOfTracks());
1261  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> Randomised_Jet=ModifyJet(Jet1,0,"Randomise");
1262  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_02_15=ModifyJet(Jet1,0,"AddExtraProng_02_15");
1263  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_02_30=ModifyJet(Jet1,0,"AddExtraProng_02_30");
1264  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_02_45=ModifyJet(Jet1,0,"AddExtraProng_02_45");
1265  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_02_60=ModifyJet(Jet1,0,"AddExtraProng_02_60");
1266  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_01_30=ModifyJet(Jet1,0,"AddExtraProng_01_30");
1267  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_01_45=ModifyJet(Jet1,0,"AddExtraProng_01_45");
1268  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_03_30=ModifyJet(Jet1,0,"AddExtraProng_03_30");
1269  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_03_45=ModifyJet(Jet1,0,"AddExtraProng_03_45");
1270  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) fShapesVar[0]= Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
1271  else fShapesVar[0]=Jet1->Pt();
1272  fShapesVar[2]=FjNSubJettiness(Jet1,0,1,0,1,0);
1273  fShapesVar[4]=FjNSubJettiness(Jet1,0,2,0,1,0);
1274  fShapesVar[6]=FjNSubJettinessFastJet(ExtraProng_Jet_02_45,0,1,0,1,0);
1275  fShapesVar[8]=FjNSubJettiness(Jet1,0,2,0,1,1);
1276  // fShapesVar[10]=Jet1->GetNumberOfTracks();
1277  fShapesVar[10]=FjNSubJettinessFastJet(ExtraProng_Jet_01_45,0,1,0,1,0);
1278  fShapesVar[12]=FjNSubJettinessFastJet(ExtraProng_Jet_03_45,0,1,0,1,0);
1279  fShapesVar[14]=FjNSubJettinessFastJet(ExtraProng_Jet_02_15,0,2,0,1,0);
1280  fShapesVar[16]=FjNSubJettinessFastJet(ExtraProng_Jet_02_45,0,2,0,1,0);
1281  fShapesVar[18]=FjNSubJettinessFastJet(ExtraProng_Jet_01_30,0,2,0,1,0);
1282  fShapesVar[20]=FjNSubJettinessFastJet(ExtraProng_Jet_03_30,0,2,0,1,0);
1283  AliEmcalJetFinder *Reclusterer1 = Recluster(Jet1, 0, fSubJetRadius, fSubJetMinPt, fSubJetAlgorithm, "SubJetFinder");
1284  if (fFullTree){
1285  fShapesVar[22]=FjNSubJettiness(Jet1,0,2,0,1,2);
1286  fShapesVar[24]=SubJetFraction(Jet1, Reclusterer1, 1, 0, kTRUE, kFALSE);
1287  fShapesVar[26]=SubJetFraction(Jet1, Reclusterer1, 2, 0, kTRUE, kFALSE);
1288  }
1289  fShapesVar[1]=FjNSubJettinessFastJet(Randomised_Jet,0,1,0,1,0);
1290  // fShapesVar[3]=FjNSubJettiness(std::unique_ptr<AliEmcalJet>(ModifyJet(Jet1,0,"Randomise")).get(),0,1,0,1,0);
1291  //fShapesVar[5]=FjNSubJettiness(std::unique_ptr<AliEmcalJet>(ModifyJet(Jet1,0,"AddExtraProng_02_30")).get(),0,1,0,1,0);
1292  fShapesVar[3]=FjNSubJettinessFastJet(ExtraProng_Jet_02_15,0,1,0,1,0);
1293  fShapesVar[5]=FjNSubJettinessFastJet(ExtraProng_Jet_02_30,0,1,0,1,0);
1294  fShapesVar[7]=FjNSubJettinessFastJet(ExtraProng_Jet_02_60,0,1,0,1,0);
1295  fShapesVar[9]=FjNSubJettinessFastJet(ExtraProng_Jet_01_30,0,1,0,1,0);
1296  fShapesVar[11]=FjNSubJettinessFastJet(ExtraProng_Jet_03_30,0,1,0,1,0);
1297  fShapesVar[13]=FjNSubJettinessFastJet(Randomised_Jet,0,2,0,1,0);
1298  fShapesVar[15]=FjNSubJettinessFastJet(ExtraProng_Jet_02_30,0,2,0,1,0);
1299  fShapesVar[17]=FjNSubJettinessFastJet(ExtraProng_Jet_02_60,0,2,0,1,0);
1300  fShapesVar[19]=FjNSubJettinessFastJet(ExtraProng_Jet_01_45,0,2,0,1,0);
1301  fShapesVar[21]=FjNSubJettinessFastJet(ExtraProng_Jet_03_45,0,2,0,1,0);
1302  if (fFullTree){
1303  fShapesVar[23]=-2;
1304  fShapesVar[25]=-2;
1305  fShapesVar[27]=-2;
1306  }
1307  fTreeResponseMatrixAxis->Fill();
1308 
1309  fhSubJetCounter->Fill(Reclusterer1->GetNumberOfJets());
1310  for (Int_t i= 0; i<Reclusterer1->GetNumberOfJets(); i++){
1311  fhEventCounter->Fill(6); //Number of overall subjets in all events
1312  fhSubJetPt->Fill(Reclusterer1->GetJet(i)->Pt());
1313  fhSubJetMass->Fill(Reclusterer1->GetJet(i)->M());
1314  fhNumberOfSubJetTracks->Fill(Reclusterer1->GetJet(i)->GetNumberOfTracks());
1315  fhSubJetRadius->Fill(TMath::Sqrt((Reclusterer1->GetJet(i)->Area()/TMath::Pi()))); //Radius of SubJets per event
1316  }
1317  delete Randomised_Jet.second;
1318  delete ExtraProng_Jet_02_15.second;
1319  delete ExtraProng_Jet_02_30.second;
1320  delete ExtraProng_Jet_02_45.second;
1321  delete ExtraProng_Jet_02_60.second;
1322  delete ExtraProng_Jet_01_30.second;
1323  delete ExtraProng_Jet_01_45.second;
1324  delete ExtraProng_Jet_03_30.second;
1325  delete ExtraProng_Jet_03_45.second;
1326  }
1327  }
1328  fhJetCounter->Fill(JetCounter); //Number of Jets in Each Event
1329  }
1330  //else {fhEventCounter->Fill(2);} //Events with no jet container
1331  }
1332  return kTRUE;
1333 }
1334 //________________________________________________________________________
1336 
1337  if(Phi < -1*TMath::Pi()) Phi += (2*TMath::Pi());
1338  else if (Phi > TMath::Pi()) Phi -= (2*TMath::Pi());
1339  Double_t DeltaPhi=Phi-EventPlane;
1340  if(DeltaPhi < -1*TMath::Pi()) DeltaPhi += (2*TMath::Pi());
1341  else if (DeltaPhi > TMath::Pi()) DeltaPhi -= (2*TMath::Pi());
1342  return DeltaPhi;
1343 }
1344 //________________________________________________________________________
1346 
1347  if(Phi1 < -1*TMath::Pi()) Phi1 += (2*TMath::Pi()); // Turns the range of 0to2Pi into -PitoPi ???????????
1348  else if (Phi1 > TMath::Pi()) Phi1 -= (2*TMath::Pi());
1349  if(Phi2 < -1*TMath::Pi()) Phi2 += (2*TMath::Pi());
1350  else if (Phi2 > TMath::Pi()) Phi2 -= (2*TMath::Pi());
1351  Double_t DeltaPhi=Phi2-Phi1;
1352  if(DeltaPhi < -1*TMath::Pi()) DeltaPhi += (2*TMath::Pi());
1353  else if (DeltaPhi > TMath::Pi()) DeltaPhi -= (2*TMath::Pi());
1354  return DeltaPhi;
1355 }
1356 
1357 
1358 //--------------------------------------------------------------------------
1360 
1361  AliTrackContainer *PartCont = NULL;
1362  AliParticleContainer *PartContMC = NULL;
1363 
1364 
1365  if (fJetShapeSub==kConstSub){
1367  else PartCont = GetTrackContainer(1);
1368  }
1369  else{
1371  else PartCont = GetTrackContainer(0);
1372  }
1373 
1374  TClonesArray *TracksArray = NULL;
1375  TClonesArray *TracksArrayMC = NULL;
1376  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) TracksArrayMC = PartContMC->GetArray();
1377  else TracksArray = PartCont->GetArray();
1378 
1380  if(!PartContMC || !TracksArrayMC) return -99999;
1381  }
1382  else {
1383  if(!PartCont || !TracksArray) return -99999;
1384  }
1385 
1386  AliAODTrack *Track = 0x0;
1387  Int_t Trigger_Index[100];
1388  for (Int_t i=0; i<100; i++) Trigger_Index[i] = 0;
1389  Int_t Trigger_Counter = 0;
1390  Int_t NTracks=0;
1391  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) NTracks = TracksArrayMC->GetEntriesFast();
1392  else NTracks = TracksArray->GetEntriesFast();
1393  for(Int_t i=0; i < NTracks; i++){
1395  if((Track = static_cast<AliAODTrack*>(PartContMC->GetAcceptParticle(i)))){
1396  if (!Track) continue;
1397  if(TMath::Abs(Track->Eta())>0.9) continue;
1398  if (Track->Pt()<0.15) continue;
1399  if ((Track->Pt() >= PtMin) && (Track->Pt()< PtMax)) {
1400  Trigger_Index[Trigger_Counter] = i;
1401  Trigger_Counter++;
1402  }
1403  }
1404  }
1405  else{
1406  if((Track = static_cast<AliAODTrack*>(PartCont->GetAcceptTrack(i)))){
1407  if (!Track) continue;
1408  if(TMath::Abs(Track->Eta())>0.9) continue;
1409  if (Track->Pt()<0.15) continue;
1410  if ((Track->Pt() >= PtMin) && (Track->Pt()< PtMax)) {
1411  Trigger_Index[Trigger_Counter] = i;
1412  Trigger_Counter++;
1413  }
1414  }
1415  }
1416  }
1417  if (Trigger_Counter == 0) return -99999;
1418  Int_t RandomNumber = 0, Index = 0 ;
1419  TRandom3* Random = new TRandom3(0);
1420  RandomNumber = Random->Integer(Trigger_Counter);
1421  Index = Trigger_Index[RandomNumber];
1422  return Index;
1423 }
1424 
1425 
1426 //--------------------------------------------------------------------------
1428 
1429  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1430  Double_t Angularity_Numerator=0; //Reset these values
1431  Double_t Angularity_Denominator=0;
1432  AliVParticle *Particle=0x0;
1433  Double_t DeltaPhi=0;
1434 
1435 
1436  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks (particles in the jet
1437  Particle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
1438  if(!Particle) continue;
1439  DeltaPhi=RelativePhi(Jet->Phi(),Particle->Phi());
1440  Angularity_Numerator=Angularity_Numerator+(Particle->Pt()*TMath::Sqrt(((Particle->Eta()-Jet->Eta())*(Particle->Eta()-Jet->Eta()))+(DeltaPhi*DeltaPhi)));
1441  Angularity_Denominator= Angularity_Denominator+Particle->Pt();
1442  }
1443  if(Angularity_Denominator!=0) return Angularity_Numerator/Angularity_Denominator;
1444  else return -1;
1445 
1446 }
1447 
1448 
1449 
1450 //--------------------------------------------------------------------------
1452 
1453  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1454  Double_t PTD_Numerator=0; //Reset these values
1455  Double_t PTD_Denominator=0;
1456  AliVParticle *Particle=0x0;
1457  Double_t DeltaPhi=0;
1458  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks (particles in the jet
1459  Particle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
1460  if(!Particle) continue;
1461  PTD_Numerator=PTD_Numerator+(Particle->Pt()*Particle->Pt());
1462  PTD_Denominator=PTD_Denominator+Particle->Pt();
1463  }
1464  if(PTD_Denominator!=0) return TMath::Sqrt(PTD_Numerator)/PTD_Denominator;
1465  else return -1;
1466 
1467 }
1468 
1469 
1470 //----------------------------------------------------------------------
1472 AliEmcalJetFinder *AliAnalysisTaskSubJetFraction::Recluster(AliEmcalJet *Jet, Int_t JetContNb, Double_t SubJetRadius, Double_t SubJetMinPt, Int_t Algorithm, const char* Name){
1473 
1474  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1475  AliEmcalJetFinder *Reclusterer = new AliEmcalJetFinder(Name); //JetFinder Object for reclustered jets
1476  Reclusterer->SetRadius(SubJetRadius);
1477  Reclusterer->SetJetMinPt(SubJetMinPt);
1478  Reclusterer->SetJetAlgorithm(Algorithm); //0 for anti-kt 1 for kt
1479  Reclusterer->SetJetMaxEta(0.9);
1480  Reclusterer->SetRecombSheme(0);
1482  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
1483  Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
1484  if(Reclusterer->AliEmcalJetFinder::Filter(Jet, JetCont, dVtx)){;} //reclustering jet1 using the jetfinderobject Reclusterer
1485  }
1486  else{
1487  Double_t dVtx[3]={1,1,1};
1488  if(Reclusterer->AliEmcalJetFinder::Filter(Jet, JetCont, dVtx)){;} //reclustering jet1 using the jetfinderobject Reclusterer
1489  }
1490  return Reclusterer;
1491 
1492 }
1493 
1494 
1495 
1496 //----------------------------------------------------------------------
1498  AliEmcalJet *SubJet=NULL;
1499  Double_t SortingVariable;
1500  Int_t ArraySize =N+1;
1501  TArrayD JetSorter(ArraySize);
1502  TArrayD JetIndexSorter(ArraySize);
1503  for (Int_t i=0; i<ArraySize; i++){
1504  JetSorter[i]=0;
1505  }
1506  for (Int_t i=0; i<ArraySize; i++){
1507  JetIndexSorter[i]=0;
1508  }
1509  if(Reclusterer->GetNumberOfJets()<N) return -999;
1510  for (Int_t i=0; i<Reclusterer->GetNumberOfJets(); i++){
1511  SubJet=Reclusterer->GetJet(i);
1512  if (Type==0) SortingVariable=SubJet->Pt();
1513  else if (Type==1) SortingVariable=SubJet->E();
1514  else if (Type==2) SortingVariable=SubJet->M();
1515  for (Int_t j=0; j<N; j++){
1516  if (SortingVariable>JetSorter[j]){
1517  for (Int_t k=N-1; k>=j; k--){
1518  JetSorter[k+1]=JetSorter[k];
1519  JetIndexSorter[k+1]=JetIndexSorter[k];
1520  }
1521  JetSorter[j]=SortingVariable;
1522  JetIndexSorter[j]=i;
1523  break;
1524  }
1525  }
1526  }
1527  if (!Index) return JetSorter[N-1];
1528  else return JetIndexSorter[N-1];
1529 }
1530 
1531 
1532 
1533 //returns -1 if the Nth hardest jet is requested where N>number of available jets
1534 //type: 0=Pt 1=E 2=M
1535 //Index TRUE=returns index FALSE=returns value of quantatiy in question
1536 
1537 
1538 
1539 
1540 
1541 
1542 
1543 //----------------------------------------------------------------------
1545  AliEmcalJet *SubJet=NULL;
1546  Double_t Observable=0;
1547  Double_t Fraction_Numerator=0;
1548  Bool_t Error=kFALSE;
1549  if (!Jet->GetNumberOfTracks()) return -2;
1550  if (Add){
1551  for (Int_t i=1; i<=N; i++){
1552  Observable=SubJetOrdering(Jet,Reclusterer,i,Type,kFALSE);
1553  if(Observable==-999){
1554  Error = kTRUE;
1555  return -2;
1556  }
1557  Fraction_Numerator=Fraction_Numerator+Observable;
1558  }
1559  }
1560  else {
1561  Fraction_Numerator=SubJetOrdering(Jet,Reclusterer,N,Type,kFALSE);
1562  if (Fraction_Numerator==-999) return -2;
1563  }
1564  if (Type==0){
1565  if(Loss) return (Jet->Pt()-Fraction_Numerator)/Jet->Pt();
1566  else return Fraction_Numerator/Jet->Pt();
1567  }
1568  else if (Type==1){
1569  if(Loss) return (Jet->E()-Fraction_Numerator)/Jet->E();
1570  else return Fraction_Numerator/Jet->E();
1571  }
1572  else { //change to else if if you want to add more later
1573  if(Loss) return (Jet->M()-Fraction_Numerator)/Jet->M();
1574  else return Fraction_Numerator/Jet->M();
1575  }
1576 }
1577 //N number of hardest subjets involved
1578 //Type 0=Pt 1=E 2=M
1579 // Add TRUE: Add 1 to Nth hardest subjets together/total jet False: Nth hardest Jet/toal jet
1580 //Loss TRUE: Jet-Subjet(s)/Jet FALSE: Subjet(s)/Jet
1581 
1582 
1583 //----------------------------------------------------------------------
1585  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1586  AliEmcalJet *SubJet=NULL;
1587  Double_t DeltaR1=0;
1588  Double_t DeltaR2=0;
1589  AliVParticle *JetParticle=0x0;
1590  Double_t SubJetiness_Numerator = 0;
1591  Double_t SubJetiness_Denominator = 0;
1592  Double_t Index=-2;
1593  Bool_t Error=kFALSE;
1594  // JetRadius=TMath::Sqrt((Jet->Area()/TMath::Pi())); //comment out later
1595  if (!Jet->GetNumberOfTracks()) return -2;
1596  if (Reclusterer->GetNumberOfJets() < N) return -2;
1597  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks (particles in the jet
1598  JetParticle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
1599  for (Int_t j=1; j<=N; j++){
1600  Index=SubJetOrdering(Jet,Reclusterer,j,0,kTRUE);
1601  if(Index==-999){
1602  Error = kTRUE;
1603  i=Jet->GetNumberOfTracks();
1604  break;
1605  }
1606  if(j==1){
1607  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);
1608  }
1609  else{
1610  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);
1611  if (DeltaR2<DeltaR1) DeltaR1=DeltaR2;
1612  }
1613  }
1614  SubJetiness_Numerator=SubJetiness_Numerator+(JetParticle->Pt()*DeltaR1);
1615  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));
1616  else SubJetiness_Denominator=SubJetiness_Denominator+(TMath::Power((Reclusterer->GetJet(SubJetOrdering(Jet,Reclusterer,N,0,kTRUE))->Pt()),A)*JetParticle->Pt()*TMath::Power(JetRadius,B));
1617  }
1618  if (SubJetiness_Denominator!=0 && !Error){
1619  // if (SubJetiness_Numerator/SubJetiness_Denominator>1.0) cout << JetRadius<<endl;
1620  return SubJetiness_Numerator/SubJetiness_Denominator; }
1621  else return -2;
1622 }
1623 
1624 
1625 //______________________________________________________________________________________
1627 
1628  //WARNING!!! Only works for parent jets that are clustered with Anti-Kt! To change go to AliEmcalJetFinder.cxx and look at the Nsubjettiness() function
1629 
1630  //Algorithm==0 -> kt_axes;
1631  // Algorithm==1 -> ca_axes;
1632  //Algorithm==2 -> antikt_0p2_axes;
1633  //Algorithm==3 -> wta_kt_axes;
1634  //Algorithm==4 -> wta_ca_axes;
1635  //Algorithm==5 -> onepass_kt_axes;
1636  //Algorithm==6 -> onepass_ca_axes;
1637  //Algorithm==7 -> onepass_antikt_0p2_axes;
1638  //Algorithm==8 -> onepass_wta_kt_axes;
1639  //Algorithm==9 -> onepass_wta_ca_axes;
1640  //Algorithm==10 -> min_axes;
1641 
1642 
1643  //fSoftDropOn==1 All the below are done on the jet after it has undergone softdrop with CA. Option 3 and 4 are done whether this is turned on or not.
1644 
1645 
1646  //Option==0 returns Nsubjettiness Value
1647  //Option==1 && N==2 returns opening angle between two subjet axes(Delta R?)
1648  //Option==2 && N==2 returns Delta R
1649  //Option==3 returns first splitting distance for soft dropped jet (CA)
1650  //Option==4 returns Symmetry measure (Zg) for soft dropped jet (CA)
1651  //Option==5 returns Pt of Subjet1
1652  //Option==6 returns Pt of Subjet2
1653  //Options==7 trutns deltaR of subjets...Is this different to before??
1654  //Option==8 Subjet1 Leading track Pt
1655  //Option==9 Subjet1 Leading track Pt
1656 
1657 
1658 
1659  if (Jet->GetNumberOfTracks()>=N){
1660  if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==1) && (Algorithm==0) && (Beta==1.0) && (Option==0)){
1663  }
1664  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==0) && (Beta==1.0) && (Option==0)){
1667  }
1668  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==3) && (Algorithm==0) && (Beta==1.0) && (Option==0)){
1671  }
1672  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==0) && (Beta==1.0) && (Option==1)){
1675  }
1676  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==1) && (Algorithm==1) && (Beta==1.0) && (Option==0)){
1679  }
1680  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==1) && (Beta==1.0) && (Option==0)){
1683  }
1684  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==1) && (Beta==1.0) && (Option==1)){
1687  }
1688  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==1) && (Algorithm==2) && (Beta==1.0) && (Option==0)){
1691  }
1692  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==2) && (Beta==1.0) && (Option==0)){
1695  }
1696  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==2) && (Beta==1.0) && (Option==1)){
1699  }
1700  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==1) && (Algorithm==1) && (Beta==1.0) && (Option==0) && (Beta_SD==0.0) && (ZCut==0.1) && (fSoftDropOn==1)){
1703  }
1704  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==1) && (Beta==1.0) && (Option==0) && (Beta_SD==0.0) && (ZCut==0.1) && (fSoftDropOn==1)){
1707  }
1708  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==1) && (Beta==1.0) && (Option==1) && (Beta_SD==0.0) && (ZCut==0.1) && (fSoftDropOn==1)){
1711  }
1712  else{
1713  Algorithm=fReclusteringAlgorithm;
1714  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1715  AliEmcalJetFinder JetFinder("Nsubjettiness");
1716  JetFinder.SetJetMaxEta(0.9-fJetRadius);
1717  JetFinder.SetRadius(fJetRadius);
1718  JetFinder.SetJetAlgorithm(0); //0 for anti-kt 1 for kt //this is for the JET!!!!!!!!!! Not the SubJets
1719  JetFinder.SetRecombSheme(0);
1720  JetFinder.SetJetMinPt(Jet->Pt());
1722  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
1723  Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
1724  return JetFinder.Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,fSubJetRadius,Beta,Option,0,Beta_SD,ZCut,fSoftDropOn);
1725  }
1726  else{
1727  Double_t dVtx[3]={1,1,1};
1728  if (!fNsubMeasure){
1729  return JetFinder.Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,fSubJetRadius,Beta,Option,0,Beta_SD,ZCut,fSoftDropOn);
1730  }
1731  else{
1732  if (Option==3) return JetFinder.Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,fSubJetRadius,Beta,Option,0,Beta_SD,ZCut,fSoftDropOn);
1733  else return JetFinder.Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,fSubJetRadius,Beta,Option,1);
1734  }
1735  }
1736  }
1737  }
1738  else return -2;
1739 }
1740 
1741 
1742 //________________________________________________________________________
1744  //
1745  // retrieve event objects
1746  //
1748  return kFALSE;
1749 
1750  return kTRUE;
1751 }
1752 
1753 
1754 //_______________________________________________________________________
1756 {
1757  // Called once at the end of the analysis.
1759 
1760  /*
1761  for (int i=1; i<=(XBinsJetPtSize)-1; i++){ //Rescales the fhJetPt Histograms according to the with of the variable bins and number of events
1762  fhJetPt->SetBinContent(i,((1.0/(XBinsJetPt[i]-XBinsJetPt[i-1]))*((fhJetPt->GetBinContent(i))*(1.0/(fhEventCounter->GetBinContent(1))))));
1763  fhSubJetPt->SetBinContent(i,((1.0/(XBinsJetPt[i]-XBinsJetPt[i-1]))*((fhSubJetPt->GetBinContent(i))*(1.0/(fhEventCounter->GetBinContent(5))))));
1764 
1765  //fhJetPt->SetBinContent(i,((1.0/(fhPt->GetBinWidth(i)))*((fhJetPt->GetBinContent(i))*(fhEventCounter->GetBinContent(1)))));
1766  // fhJetPt->SetBinContent(i,((1.0/((fhPt->GetLowEdge(i+1))-(fhJetPt->GetLowEdge(i))))*((fhPt->GetBinContent(i))*(fhEventCounter->GetBinContent(1)))));
1767  }
1768  for (int i=1; i<=(XBinsJetMassSize)-1; i++){ //Rescales the fhPt Histograms according to the with of the variable bins and number of events
1769  fhJetMass->SetBinContent(i,((1.0/(XBinsJetMass[i]-XBinsJetMass[i-1]))*((fhJetMass->GetBinContent(i))*(1.0/(fhEventCounter->GetBinContent(1))))));
1770  }
1771 
1772  fhJetPhi->Scale(Phi_Bins/((Phi_Up-Phi_Low)*((fhEventCounter->GetBinContent(1)))));
1773  fhJetEta->Scale(Eta_Bins/((Eta_Up-Eta_Low)*((fhEventCounter->GetBinContent(1)))));
1774  fhJetRadius->Scale(100/(fhEventCounter->GetBinContent(4))); //should this and JetAngularity be divided by Bin 1 or 4????
1775  fhNumberOfJetTracks->Scale(1.0/(fhEventCounter->GetBinContent(4)));
1776 
1777  fhJetCounter->Scale(1.0/(fhEventCounter->GetBinContent(1))); //is the first bin the correct one to look at?
1778 
1779 
1780  fhSubJetCounter->Scale(1.0/(fhEventCounter->GetBinContent(5)));
1781 
1782  */
1783  /*
1784 
1785  fhJetPt->Scale(1.0/(fhEventCounter->GetBinContent(1)));
1786  fhSubJetPt->Scale(1.0/(fhEventCounter->GetBinContent(5)));
1787  fhJetMass->Scale(1.0/(fhEventCounter->GetBinContent(1)));
1788  fhJetPhi->Scale(1.0/(fhEventCounter->GetBinContent(1)));
1789  fhJetEta->Scale(1.0/(fhEventCounter->GetBinContent(1)));
1790  fhJetRadius->Scale(1.0/(fhEventCounter->GetBinContent(4))); //should this and JetAngularity be divided by Bin 1 or 4????
1791  fhNumberOfJetTracks->Scale(1.0/(fhEventCounter->GetBinContent(4)));
1792  fhJetCounter->Scale(1.0/(fhEventCounter->GetBinContent(1))); //is the first bin the correct one to look at?
1793  fhSubJetCounter->Scale(1.0/(fhEventCounter->GetBinContent(5)));
1794  */
1795  }
1796 
1797 }
1798 
1800 std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> AliAnalysisTaskSubJetFraction::ModifyJet(AliEmcalJet* Jet, Int_t JetContNb, TString Modification){
1801  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> Jet_ClusterSequence;
1802  std::vector<fastjet::PseudoJet> fInputVectors;
1803  AliJetContainer *fJetCont = GetJetContainer(JetContNb);
1804  AliParticleContainer *fTrackCont = fJetCont->GetParticleContainer();
1805  std::vector<fastjet::PseudoJet> Modified_Jet;
1806  fastjet::ClusterSequence *fClustSeqSA;
1807  Int_t NJetTracksTest = Jet->GetNumberOfTracks();
1808  if (fTrackCont) for (Int_t i=0; i<Jet->GetNumberOfTracks(); i++) {
1809  AliVParticle *fTrk = Jet->TrackAt(i, fTrackCont->GetArray());
1810  if (!fTrk) continue;
1811  fastjet::PseudoJet PseudoTracks(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
1812  PseudoTracks.set_user_index(Jet->TrackAt(i)+100); //100 is very arbitary....why is it here anyway?
1813  fInputVectors.push_back(PseudoTracks);
1814  }
1815  if (Modification=="Randomise") fInputVectors=RandomiseTracks(Jet,fInputVectors);
1816  if (Modification=="AddExtraProng_02_15") fInputVectors=AddExtraProng(fInputVectors,0.2,0.15);
1817  if (Modification=="AddExtraProng_02_30") fInputVectors=AddExtraProng(fInputVectors,0.2,0.30);
1818  if (Modification=="AddExtraProng_02_45") fInputVectors=AddExtraProng(fInputVectors,0.2,0.45);
1819  if (Modification=="AddExtraProng_02_60") fInputVectors=AddExtraProng(fInputVectors,0.2,0.60);
1820  if (Modification=="AddExtraProng_01_30") fInputVectors=AddExtraProng(fInputVectors,0.1,0.30);
1821  if (Modification=="AddExtraProng_01_45") fInputVectors=AddExtraProng(fInputVectors,0.1,0.45);
1822  if (Modification=="AddExtraProng_03_30") fInputVectors=AddExtraProng(fInputVectors,0.3,0.30);
1823  if (Modification=="AddExtraProng_03_45") fInputVectors=AddExtraProng(fInputVectors,0.3,0.45);
1824  fJetRadius=fJetRadius*2.0;
1825  try {
1826  fastjet::JetDefinition fJetDef(fastjet::antikt_algorithm, fJetRadius, static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
1827  fClustSeqSA=new fastjet::ClusterSequence(fInputVectors, fJetDef);
1828  Modified_Jet=fClustSeqSA->inclusive_jets(0);
1829  Jet_ClusterSequence.second=fClustSeqSA;
1830  } catch (fastjet::Error) {
1831  AliError(" [w] FJ Exception caught.");
1832  //return -1;
1833  }
1834  fJetRadius=fJetRadius/2.0;
1835  Jet_ClusterSequence.first=Modified_Jet[0];
1836  return Jet_ClusterSequence;
1837 }
1838 
1839 
1840 std::vector<fastjet::PseudoJet> AliAnalysisTaskSubJetFraction::RandomiseTracks(AliEmcalJet *Jet,std::vector<fastjet::PseudoJet> fInputVectors){
1841  TRandom3 fRandom1;
1842  fRandom1.SetSeed(0);
1843  Double_t Random_Phi=0.0;
1844  Double_t Random_Phi_Change=0.0;
1845  Double_t Random_Eta=0.0;
1846  Double_t Random_Eta_Change=0.0;
1847  Double_t Track_Pt=0.0;
1848  std::vector<fastjet::PseudoJet> Random_Track_Vector;
1849  Random_Track_Vector.clear();
1850  for (Int_t i=0; i< fInputVectors.size(); i++){
1851  Random_Phi_Change=fRandom1.Uniform(-1*fJetRadius,fJetRadius);
1852  Random_Phi=Jet->Phi()+Random_Phi_Change;
1853  Random_Eta_Change=fRandom1.Uniform(-1*(TMath::Sqrt((fJetRadius*fJetRadius)-(Random_Phi_Change*Random_Phi_Change))),TMath::Sqrt((fJetRadius*fJetRadius)-(Random_Phi_Change*Random_Phi_Change)));
1854  Random_Eta=Jet->Eta()+Random_Eta_Change;
1855  Track_Pt=fInputVectors[i].perp();
1856  fastjet::PseudoJet Random_Track(Track_Pt*TMath::Cos(Random_Phi),Track_Pt*TMath::Sin(Random_Phi),Track_Pt*TMath::SinH(Random_Eta),TMath::Sqrt((Track_Pt*Track_Pt)+(Track_Pt*TMath::SinH(Random_Eta)*Track_Pt*TMath::SinH(Random_Eta))));
1857  Random_Track.set_user_index(i);
1858  Random_Track_Vector.push_back(Random_Track);
1859  }
1860  fInputVectors.clear();
1861  return Random_Track_Vector;
1862 }
1863 
1864 std::vector<fastjet::PseudoJet> AliAnalysisTaskSubJetFraction::AddExtraProng(std::vector<fastjet::PseudoJet> fInputVectors, Double_t Distance, Double_t PtFrac){
1865  TRandom3 fRandom1;
1866  fRandom1.SetSeed(0);
1867  std::vector<fastjet::PseudoJet> Jet;
1868  Jet.clear();
1869  try {
1870  fastjet::JetDefinition fJetDef(fastjet::antikt_algorithm, fJetRadius*2, static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
1871  fastjet::ClusterSequence fClustSeqSA_Prong(fInputVectors, fJetDef);
1872  Jet= fClustSeqSA_Prong.inclusive_jets(0);
1873  } catch (fastjet::Error) {
1874  AliError(" [w] FJ Exception caught.");
1875  }
1876  Double_t Extra_Track_Phi_Change=fRandom1.Uniform(-1*Distance,Distance);
1877  Double_t Extra_Track_Phi=Jet[0].phi()+Extra_Track_Phi_Change;
1878  Double_t Extra_Track_Eta_Change=fRandom1.Uniform(-1*(TMath::Sqrt((Distance*Distance)-(Extra_Track_Phi_Change*Extra_Track_Phi_Change))),TMath::Sqrt((Distance*Distance)-(Extra_Track_Phi_Change*Extra_Track_Phi_Change)));
1879  Double_t Extra_Track_Eta=Jet[0].pseudorapidity()+Extra_Track_Eta_Change;
1880  Double_t Extra_Track_Pt=Jet[0].perp()*PtFrac;
1881  fastjet::PseudoJet ExtraTrack(Extra_Track_Pt*TMath::Cos(Extra_Track_Phi),Extra_Track_Pt*TMath::Sin(Extra_Track_Phi),Extra_Track_Pt*TMath::SinH(Extra_Track_Eta),TMath::Sqrt((Extra_Track_Pt*Extra_Track_Pt)+(Extra_Track_Pt*TMath::SinH(Extra_Track_Eta)*Extra_Track_Pt*TMath::SinH(Extra_Track_Eta))));
1882  ExtraTrack.set_user_index(150);
1883  fInputVectors.push_back(ExtraTrack);
1884  return fInputVectors;
1885 }
1887 
1888 //______________________________________________________________________________________
1889 Double_t AliAnalysisTaskSubJetFraction::FjNSubJettinessFastJet(std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> Jet_ClusterSequence, Int_t JetContNb,Int_t N, Int_t Algorithm, Double_t Beta, Int_t Option, Double_t Beta_SD, Double_t ZCut){
1890 
1891  //Only Works for kGenOnTheFly
1892 
1893  //WARNING!!! Only works for parent jets that are clustered with Anti-Kt! To change go to AliEmcalJetFinder.cxx and look at the Nsubjettiness() function
1894 
1895  //Algorithm==0 -> kt_axes;
1896  // Algorithm==1 -> ca_axes;
1897  //Algorithm==2 -> antikt_0p2_axes;
1898  //Algorithm==3 -> wta_kt_axes;
1899  //Algorithm==4 -> wta_ca_axes;
1900  //Algorithm==5 -> onepass_kt_axes;
1901  //Algorithm==6 -> onepass_ca_axes;
1902  //Algorithm==7 -> onepass_antikt_0p2_axes;
1903  //Algorithm==8 -> onepass_wta_kt_axes;
1904  //Algorithm==9 -> onepass_wta_ca_axes;
1905  //Algorithm==10 -> min_axes;
1906 
1907 
1908  //fSoftDropOn==1 All the below are done on the jet after it has undergone softdrop with CA. Option 3 and 4 are done whether this is turned on or not.
1909 
1910 
1911  //Option==0 returns Nsubjettiness Value
1912  //Option==1 && N==2 returns opening angle between two subjet axes(Delta R?)
1913  //Option==2 && N==2 returns Delta R
1914  //Option==3 returns first splitting distance for soft dropped jet (CA)
1915  //Option==4 returns Symmetry measure (Zg) for soft dropped jet (CA)
1916  //Option==5 returns Pt of Subjet1
1917  //Option==6 returns Pt of Subjet2
1918  //Options==7 trutns deltaR of subjets...Is this different to before??
1919  //Option==8 Subjet1 Leading track Pt
1920  //Option==9 Subjet1 Leading track Pt
1921 
1922  fastjet::PseudoJet Jet=Jet_ClusterSequence.first;
1923  AliFJWrapper FJWrapper("FJWrapper", "FJWrapper");
1924  FJWrapper.SetR(fJetRadius);
1925  FJWrapper.SetMinJetPt(Jet.perp());
1926  FJWrapper.SetAlgorithm(fastjet::antikt_algorithm); //this is for the jet clustering not the subjet reclustering.
1927  FJWrapper.SetRecombScheme(static_cast<fastjet::RecombinationScheme>(0));
1928  if (Jet.constituents().size()>=N){
1929  Algorithm=fReclusteringAlgorithm;
1930  Double_t dVtx[3]={1,1,1};
1931  return FJWrapper.NSubjettinessDerivativeSub(N,Algorithm,fSubJetRadius,Beta,fJetRadius,Jet,Option,0,Beta_SD,ZCut,fSoftDropOn); //change this
1932  }
1933  else return -2;
1934 }
1935 
void SetRadius(Double_t val)
Bool_t fCentSelectOn
Random number generator.
Double_t GetSecondOrderSubtractedOpeningAngle_casd() const
Double_t Area() const
Definition: AliEmcalJet.h:130
double Double_t
Definition: External.C:58
Double_t XBinsPtSize
void SetJetMinPt(Double_t val)
Double_t GetSecondOrderSubtracted1subjettiness_akt02() const
static Double_t DeltaPhi(Double_t phia, Double_t phib, Double_t rMin=-TMath::Pi()/2, Double_t rMax=3 *TMath::Pi()/2)
Definition: External.C:236
AliVParticle * GetLeadingTrack(TClonesArray *tracks=0) const
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:327
Bool_t Run()
Run function. This is the core function of the analysis and contains the user code. Therefore users have to implement this function.
AliJetContainer * GetJetContainer(Int_t i=0) const
Double_t GetFirstOrderSubtracted2subjettiness_casd() const
Double_t Eta() const
Definition: AliEmcalJet.h:121
Double_t Phi() const
Definition: AliEmcalJet.h:117
Container with name, TClonesArray and cuts for particles.
std::pair< fastjet::PseudoJet, fastjet::ClusterSequence * > ModifyJet(AliEmcalJet *Jet, Int_t JetContNb, TString Modification)
Double_t PTD(AliEmcalJet *Jet, Int_t JetContNb)
Double_t GetSecondOrderSubtracted2subjettiness_akt02() const
Bool_t RetrieveEventObjects()
Retrieve common objects from event.
Double_t fEPV0
!event plane V0
Int_t GetLabel() const
Definition: AliEmcalJet.h:124
std::vector< fastjet::PseudoJet > AddExtraProng(std::vector< fastjet::PseudoJet > fInputVectors, Double_t Distance, Double_t PtFrac)
Double_t GetFirstOrderSubtracted2subjettiness_kt() const
void SetMinJetPt(Double_t MinPt)
Definition: AliFJWrapper.h:142
Double_t NSubJettiness(AliEmcalJet *Jet, Int_t JetContNb, Double_t JetRadius, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t A, Int_t B)
Double_t GetFirstOrderSubtractedOpeningAngle_kt() const
void SetRecombScheme(const fastjet::RecombinationScheme &scheme)
Definition: AliFJWrapper.h:122
Double_t GetSecondOrderSubtractedOpeningAngle_akt02() const
Double_t RelativePhi(Double_t Phi1, Double_t Phi2)
Double_t E() const
Definition: AliEmcalJet.h:119
Bool_t FillHistograms()
Function filling histograms.
Double_t GetFirstOrderSubtracted2subjettiness_ca() const
Container for particles within the EMCAL framework.
Double_t XBinsJetMass[28]
Double_t GetFirstOrderSubtractedOpeningAngle_ca() const
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:160
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:139
void SetR(Double_t r)
Definition: AliFJWrapper.h:127
Double_t GetFirstOrderSubtracted3subjettiness_kt() const
TString kData
Declare data MC or deltaAOD.
Double_t GetSecondOrderSubtracted2subjettiness_kt() const
Double_t Nsubjettiness(AliEmcalJet *pJet, AliJetContainer *pContJets, Double_t dVtx[3], Int_t N, Int_t Algorithm, Double_t Radius, Double_t Beta, Int_t Option=0, Int_t Measure=0, Double_t Beta_SD=0, Double_t ZCut=0.1, Int_t SoftDropOn=0)
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
AliParticleContainer * GetParticleContainer() const
void SetRecombSheme(Int_t val)
void SetJetAlgorithm(Int_t val)
int Int_t
Definition: External.C:63
Double_t GetSecondOrderSubtracted1subjettiness_ca() const
float Float_t
Definition: External.C:68
Double_t GetFirstOrderSubtracted2subjettiness_akt02() const
Int_t GetNJets() const
void SetAlgorithm(const fastjet::JetAlgorithm &algor)
Definition: AliFJWrapper.h:121
Int_t SelectTriggerHadron(Float_t PtMin, Float_t PtMax)
Definition: External.C:228
Definition: External.C:212
virtual AliVParticle * GetAcceptParticle(Int_t i=-1) const
Definition: Option.C:68
Double_t GetSecondOrderSubtracted2subjettiness_ca() const
Double_t fCent
!event centrality
Double_t GetSecondOrderSubtracted2subjettiness_casd() const
Double_t SubJetFraction(AliEmcalJet *Jet, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t Type, Bool_t Add, Bool_t Loss)
void SetJetMaxEta(Double_t val)
Double_t GetSecondOrderSubtracted1subjettiness_casd() const
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 FjNSubJettinessFastJet(std::pair< fastjet::PseudoJet, fastjet::ClusterSequence * > Jet_ClusterSequence, Int_t JetContNb, Int_t N, Int_t Algorithm, Double_t Beta, Int_t Option=0, Double_t Beta_SD=0, Double_t ZCut=0.1)
Double_t Pt() const
Definition: AliEmcalJet.h:109
Double_t GetRhoVal(Int_t i=0) const
Double_t GetFirstOrderSubtractedOpeningAngle_akt02() const
Double_t GetFirstOrderSubtracted1subjettiness_akt02() const
Double_t FjNSubJettiness(AliEmcalJet *Jet, Int_t JetContNb, Int_t N, Int_t Algorithm, Double_t Beta, Int_t Option=0, Double_t Beta_SD=0, Double_t ZCut=0.1)
Double_t XBinsPt[66]
AliEmcalList * fOutput
!output list
std::vector< fastjet::PseudoJet > RandomiseTracks(AliEmcalJet *Jet, std::vector< fastjet::PseudoJet > fInputVectors)
const Int_t nVar
Double_t GetFirstOrderSubtracted1subjettiness_casd() const
AliTrackContainer * GetTrackContainer(Int_t i=0) const
void SetMakeGeneralHistograms(Bool_t g)
virtual AliVTrack * GetAcceptTrack(Int_t i=-1) const
Base task in the EMCAL jet framework.
Double_t GetFirstOrderSubtracted1subjettiness_ca() const
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
Double_t GetFirstOrderSubtracted1subjettiness_kt() const
Double_t Angularity(AliEmcalJet *Jet, Int_t JetContNb)
const char Option_t
Definition: External.C:48
Double32_t NSubjettinessDerivativeSub(Int_t N, Int_t Algorithm, Double_t Radius, Double_t Beta, Double_t JetR, fastjet::PseudoJet jet, Int_t Option=0, Int_t Measure=0, Double_t Beta_SD=0.0, Double_t ZCut=0.1, Int_t SoftDropOn=0)
Double_t RelativePhiEventPlane(Double_t EventPlane, Double_t Phi)
void UserCreateOutputObjects()
Main initialization function on the worker.
bool Bool_t
Definition: External.C:53
Double_t SubJetOrdering(AliEmcalJet *Jet, AliEmcalJetFinder *Reclusterer, Int_t N, Int_t Type, Bool_t Index)
AliEmcalJetShapeProperties * GetShapeProperties() const
Definition: AliEmcalJet.h:361
Double_t GetSecondOrderSubtracted3subjettiness_kt() const
Double_t XBinsJetPt[35]
Double_t M() const
Definition: AliEmcalJet.h:120
Double_t GetSecondOrderSubtracted1subjettiness_kt() const
Container for jet within the EMCAL jet framework.
Double_t GetFirstOrderSubtractedOpeningAngle_casd() const
Double_t GetSecondOrderSubtractedOpeningAngle_ca() const
Double_t GetSecondOrderSubtractedOpeningAngle_kt() const
AliEmcalJet * GetJet(Int_t i) const