AliPhysics  45476af (45476af)
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  fhJetPhi= new TH1F("fhJetPhi", "Jet Phi",780 , -7, 7);
512  fOutput->Add(fhJetPhi);
513  fhJetEta= new TH1F("fhJetEta", "Jet Eta", Eta_Bins, Eta_Low, Eta_Up);
514  fOutput->Add(fhJetEta);
515  fhJetMass= new TH1F("fhJetMass", "Jet Mass", 4000,-0.5, 39.5);
516  fOutput->Add(fhJetMass);
517  fhJetRadius= new TH1F("fhJetRadius", "Jet Radius", 100, -0.05,0.995);
518  fOutput->Add(fhJetRadius);
519  fhNumberOfJetTracks= new TH1F("fhNumberOfJetTracks", "Number of Tracks within a Jet", 300, -0.5,299.5);
521  fhSubJetRadius= new TH1F("fhSubJetRadius", "SubJet Radius", 100, -0.05,0.995);
522  fOutput->Add(fhSubJetRadius);
523  fhSubJetPt= new TH1F("fhSubJetPt", "SubJet Pt", 1500, -0.5,149.5);
524  fOutput->Add(fhSubJetPt);
525  fhSubJetMass= new TH1F("fhSubJetMass", "Sub Jet Mass", 4000,-0.5, 39.5);
526  fOutput->Add(fhSubJetMass);
527  fhNumberOfSubJetTracks= new TH1F("fhNumberOfSubJetTracks", "Number of Tracks within a Sub Jet", 300, -0.5,299.5);
529  fhJetCounter= new TH1F("fhJetCounter", "Jet Counter", 150, -0.5, 149.5);
530  fOutput->Add(fhJetCounter);
531  fhSubJetCounter = new TH1F("fhSubJetCounter", "SubJet Counter",50, -0.5,49.5);
532  fOutput->Add(fhSubJetCounter);
533  fhEventCounter= new TH1F("fhEventCounter", "Event Counter", 15,0.5,15.5);
534  fOutput->Add(fhEventCounter);
535  }
537  fhSubJettiness1_FJ_KT= new TH1D("fhSubJettiness1_FJ_KT","fhSubJettiness1_FJ_KT",400,-2,2);
539  fhSubJettiness1_FJ_MIN= new TH1D("fhSubJettiness1_FJ_MIN","fhSubJettiness1_FJ_MIN",400,-2,2);
541  fhSubJettiness2_FJ_KT= new TH1D("fhSubJettiness2_FJ_KT","fhSubJettiness2_FJ_KT",400,-2,2);
543  fhSubJettiness2_FJ_MIN= new TH1D("fhSubJettiness2_FJ_MIN","fhSubJettiness2_FJ_MIN",400,-2,2);
545  fhSubJettiness1CheckRatio_FJ_AKT = new TH2D("fhSubJettiness1CheckRatio_FJ_AKT","fhSubJettiness1CheckRatio_FJ_AKT",400,-2,2,300,-1,2);
547  fhSubJettiness1CheckRatio_FJ_KT= new TH2D("fhSubJettiness1CheckRatio_FJ_KT","fhSubJettiness1CheckRatio_FJ_KT",400,-2,2,300,-1,2);
549  fhSubJettiness1CheckRatio_FJ_CA= new TH2D("fhSubJettiness1CheckRatio_FJ_CA","fhSubJettiness1CheckRatio_FJ_CA",400,-2,2,300,-1,2);
551  fhSubJettiness1CheckRatio_FJ_WTA_KT= new TH2D("fhSubJettiness1CheckRatio_FJ_WTA_KT","fhSubJettiness1CheckRatio_FJ_WTA_KT",400,-2,2,300,-1,2);
553  fhSubJettiness1CheckRatio_FJ_WTA_CA= new TH2D("fhSubJettiness1CheckRatio_FJ_WTA_CA","fhSubJettiness1CheckRatio_FJ_WTA_CA",400,-2,2,300,-1,2);
555  fhSubJettiness1CheckRatio_FJ_OP_AKT= new TH2D("fhSubJettiness1CheckRatio_FJ_OP_AKT","fhSubJettiness1CheckRatio_FJ_OP_AKT",400,-2,2,300,-1,2);
557  fhSubJettiness1CheckRatio_FJ_OP_KT= new TH2D("fhSubJettiness1CheckRatio_FJ_OP_KT","fhSubJettiness1CheckRatio_FJ_OP_KT",400,-2,2,300,-1,2);
559  fhSubJettiness1CheckRatio_FJ_OP_CA= new TH2D("fhSubJettiness1CheckRatio_FJ_OP_CA","fhSubJettiness1CheckRatio_FJ_OP_CA",400,-2,2,300,-1,2);
561  fhSubJettiness1CheckRatio_FJ_OP_WTA_KT= new TH2D("fhSubJettiness1CheckRatio_FJ_OP_WTA_KT","fhSubJettiness1CheckRatio_FJ_OP_WTA_KT",400,-2,2,300,-1,2);
563  fhSubJettiness1CheckRatio_FJ_OP_WTA_CA= new TH2D("fhSubJettiness1CheckRatio_FJ_OP_WTA_CA","fhSubJettiness1CheckRatio_FJ_OP_WTA_CA",400,-2,2,300,-1,2);
565  fhSubJettiness1CheckRatio_FJ_MIN= new TH2D("fhSubJettiness1CheckRatio_FJ_MIN","fhSubJettiness1CheckRatio_FJ_MIN",400,-2,2,300,-1,2);
567  fhSubJettiness2CheckRatio_FJ_AKT = new TH2D("fhSubJettiness2CheckRatio_FJ_AKT","fhSubJettiness2CheckRatio_FJ_AKT",400,-2,2,300,-1,2);
569  fhSubJettiness2CheckRatio_FJ_KT= new TH2D("fhSubJettiness2CheckRatio_FJ_KT","fhSubJettiness2CheckRatio_FJ_KT",400,-2,2,300,-1,2);
571  fhSubJettiness2CheckRatio_FJ_CA= new TH2D("fhSubJettiness2CheckRatio_FJ_CA","fhSubJettiness2CheckRatio_FJ_CA",400,-2,2,300,-1,2);
573  fhSubJettiness2CheckRatio_FJ_WTA_KT= new TH2D("fhSubJettiness2CheckRatio_FJ_WTA_KT","fhSubJettiness2CheckRatio_FJ_WTA_KT",400,-2,2,300,-1,2);
575  fhSubJettiness2CheckRatio_FJ_WTA_CA= new TH2D("fhSubJettiness2CheckRatio_FJ_WTA_CA","fhSubJettiness2CheckRatio_FJ_WTA_CA",400,-2,2,300,-1,2);
577  fhSubJettiness2CheckRatio_FJ_OP_AKT= new TH2D("fhSubJettiness2CheckRatio_FJ_OP_AKT","fhSubJettiness2CheckRatio_FJ_OP_AKT",400,-2,2,300,-1,2);
579  fhSubJettiness2CheckRatio_FJ_OP_KT= new TH2D("fhSubJettiness2CheckRatio_FJ_OP_KT","fhSubJettiness2CheckRatio_FJ_OP_KT",400,-2,2,300,-1,2);
581  fhSubJettiness2CheckRatio_FJ_OP_CA= new TH2D("fhSubJettiness2CheckRatio_FJ_OP_CA","fhSubJettiness2CheckRatio_FJ_OP_CA",400,-2,2,300,-1,2);
583  fhSubJettiness2CheckRatio_FJ_OP_WTA_KT= new TH2D("fhSubJettiness2CheckRatio_FJ_OP_WTA_KT","fhSubJettiness2CheckRatio_FJ_OP_WTA_KT",400,-2,2,300,-1,2);
585  fhSubJettiness2CheckRatio_FJ_OP_WTA_CA= new TH2D("fhSubJettiness2CheckRatio_FJ_OP_WTA_CA","fhSubJettiness2CheckRatio_FJ_OP_WTA_CA",400,-2,2,300,-1,2);
587  fhSubJettiness2CheckRatio_FJ_MIN= new TH2D("fhSubJettiness2CheckRatio_FJ_MIN","fhSubJettiness2CheckRatio_FJ_MIN",400,-2,2,300,-1,2);
589  fhSubJettiness2to1CheckRatio_FJ_AKT = new TH2D("fhSubJettiness2to1CheckRatio_FJ_AKT","fhSubJettiness2to1CheckRatio_FJ_AKT",400,-2,2,300,-1,2);
591  fhSubJettiness2to1CheckRatio_FJ_KT= new TH2D("fhSubJettiness2to1CheckRatio_FJ_KT","fhSubJettiness2to1CheckRatio_FJ_KT",400,-2,2,300,-1,2);
593  fhSubJettiness2to1CheckRatio_FJ_CA= new TH2D("fhSubJettiness2to1CheckRatio_FJ_CA","fhSubJettiness2to1CheckRatio_FJ_CA",400,-2,2,300,-1,2);
595  fhSubJettiness2to1CheckRatio_FJ_WTA_KT= new TH2D("fhSubJettiness2to1CheckRatio_FJ_WTA_KT","fhSubJettiness2to1CheckRatio_FJ_WTA_KT",400,-2,2,300,-1,2);
597  fhSubJettiness2to1CheckRatio_FJ_WTA_CA= new TH2D("fhSubJettiness2to1CheckRatio_FJ_WTA_CA","fhSubJettiness2to1CheckRatio_FJ_WTA_CA",400,-2,2,300,-1,2);
599  fhSubJettiness2to1CheckRatio_FJ_OP_AKT= new TH2D("fhSubJettiness2to1CheckRatio_FJ_OP_AKT","fhSubJettiness2to1CheckRatio_FJ_OP_AKT",400,-2,2,300,-1,2);
601  fhSubJettiness2to1CheckRatio_FJ_OP_KT= new TH2D("fhSubJettiness2to1CheckRatio_FJ_OP_KT","fhSubJettiness2to1CheckRatio_FJ_OP_KT",400,-2,2,300,-1,2);
603  fhSubJettiness2to1CheckRatio_FJ_OP_CA= new TH2D("fhSubJettiness2to1CheckRatio_FJ_OP_CA","fhSubJettiness2to1CheckRatio_FJ_OP_CA",400,-2,2,300,-1,2);
605  fhSubJettiness2to1CheckRatio_FJ_OP_WTA_KT= new TH2D("fhSubJettiness2to1CheckRatio_FJ_OP_WTA_KT","fhSubJettiness2to1CheckRatio_FJ_OP_WTA_KT",400,-2,2,300,-1,2);
607  fhSubJettiness2to1CheckRatio_FJ_OP_WTA_CA= new TH2D("fhSubJettiness2to1CheckRatio_FJ_OP_WTA_CA","fhSubJettiness2to1CheckRatio_FJ_OP_WTA_CA",400,-2,2,300,-1,2);
609  fhSubJettiness2to1CheckRatio_FJ_MIN= new TH2D("fhSubJettiness2to1CheckRatio_FJ_MIN","fhSubJettiness2to1CheckRatio_FJ_MIN",400,-2,2,300,-1,2);
611  fhSubJettiness1_FJ_AKT= new TH1D("fhSubJettiness1_FJ_AKT","fhSubJettiness1_FJ_AKT",400,-2,2);
613  fhSubJettiness1_FJ_CA= new TH1D("fhSubJettiness1_FJ_CA","fhSubJettiness1_FJ_CA",400,-2,2);
615  fhSubJettiness1_FJ_WTA_KT= new TH1D("fhSubJettiness1_FJ_WTA_KT","fhSubJettiness1_FJ_WTA_KT",400,-2,2);
617  fhSubJettiness1_FJ_WTA_CA= new TH1D("fhSubJettiness1_FJ_WTA_CA","fhSubJettiness1_FJ_WTA_CA",400,-2,2);
619  fhSubJettiness1_FJ_OP_AKT= new TH1D("fhSubJettiness1_FJ_OP_AKT","fhSubJettiness1_FJ_OP_AKT",400,-2,2);
621  fhSubJettiness1_FJ_OP_KT= new TH1D("fhSubJettiness1_FJ_OP_KT","fhSubJettiness1_FJ_OP_KT",400,-2,2);
623  fhSubJettiness1_FJ_OP_CA= new TH1D("fhSubJettiness1_FJ_OP_CA","fhSubJettiness1_FJ_OP_CA",400,-2,2);
625  fhSubJettiness1_FJ_OP_WTA_KT= new TH1D("fhSubJettiness1_FJ_OP_WTA_KT","fhSubJettiness1_FJ_OP_WTA_KT",400,-2,2);
627  fhSubJettiness1_FJ_OP_WTA_CA= new TH1D("fhSubJettiness1_FJ_OP_WTA_CA","fhSubJettiness1_FJ_OP_WTA_CA",400,-2,2);
629  fhSubJettiness2_FJ_AKT= new TH1D("fhSubJettiness2_FJ_AKT","fhSubJettiness2_FJ_AKT",400,-2,2);
631  fhSubJettiness2_FJ_CA= new TH1D("fhSubJettiness2_FJ_CA","fhSubJettiness2_FJ_CA",400,-2,2);
633  fhSubJettiness2_FJ_WTA_KT= new TH1D("fhSubJettiness2_FJ_WTA_KT","fhSubJettiness2_FJ_WTA_KT",400,-2,2);
635  fhSubJettiness2_FJ_WTA_CA= new TH1D("fhSubJettiness2_FJ_WTA_CA","fhSubJettiness2_FJ_WTA_CA",400,-2,2);
637  fhSubJettiness2_FJ_OP_AKT= new TH1D("fhSubJettiness2_FJ_OP_AKT","fhSubJettiness2_FJ_OP_AKT",400,-2,2);
639  fhSubJettiness2_FJ_OP_KT= new TH1D("fhSubJettiness2_FJ_OP_KT","fhSubJettiness2_FJ_OP_KT",400,-2,2);
641  fhSubJettiness2_FJ_OP_CA= new TH1D("fhSubJettiness2_FJ_OP_CA","fhSubJettiness2_FJ_OP_CA",400,-2,2);
643  fhSubJettiness2_FJ_OP_WTA_KT= new TH1D("fhSubJettiness2_FJ_OP_WTA_KT","fhSubJettiness2_FJ_OP_WTA_KT",400,-2,2);
645  fhSubJettiness2_FJ_OP_WTA_CA= new TH1D("fhSubJettiness2_FJ_OP_WTA_CA","fhSubJettiness2_FJ_OP_WTA_CA",400,-2,2);
647  fhSubJettiness2to1_FJ_KT= new TH1D("fhSubJettiness2to1_FJ_KT","fhSubJettiness2to1_FJ_KT",400,-2,2);
649  fhSubJettiness2to1_FJ_AKT= new TH1D("fhSubJettiness2to1_FJ_AKT","fhSubJettiness2to1_FJ_AKT",400,-2,2);
651  fhSubJettiness2to1_FJ_CA= new TH1D("fhSubJettiness2to1_FJ_CA","fhSubJettiness2to1_FJ_CA",400,-2,2);
653  fhSubJettiness2to1_FJ_WTA_KT= new TH1D("fhSubJettiness2to1_FJ_WTA_KT","fhSubJettiness2to1_FJ_WTA_KT",400,-2,2);
655  fhSubJettiness2to1_FJ_WTA_CA= new TH1D("fhSubJettiness2to1_FJ_WTA_CA","fhSubJettiness2to1_FJ_WTA_CA",400,-2,2);
657  fhSubJettiness2to1_FJ_OP_AKT= new TH1D("fhSubJettiness2to1_FJ_OP_AKT","fhSubJettiness2to1_FJ_OP_AKT",400,-2,2);
659  fhSubJettiness2to1_FJ_OP_KT= new TH1D("fhSubJettiness2to1_FJ_OP_KT","fhSubJettiness2to1_FJ_OP_KT",400,-2,2);
661  fhSubJettiness2to1_FJ_OP_CA= new TH1D("fhSubJettiness2to1_FJ_OP_CA","fhSubJettiness2to1_FJ_OP_CA",400,-2,2);
663  fhSubJettiness2to1_FJ_OP_WTA_KT= new TH1D("fhSubJettiness2to1_FJ_OP_WTA_KT","fhSubJettiness2to1_FJ_OP_WTA_KT",400,-2,2);
665  fhSubJettiness2to1_FJ_OP_WTA_CA= new TH1D("fhSubJettiness2to1_FJ_OP_WTA_CA","fhSubJettiness2to1_FJ_OP_WTA_CA",400,-2,2);
667  fhSubJettiness2to1_FJ_MIN= new TH1D("fhSubJettiness2to1_FJ_MIN","fhSubJettiness2to1_FJ_MIN",400,-2,2);
669  }
671  fhJetPt_1= new TH1F("fhJetPt_1", "Jet Pt Detector Level",1500,-0.5,149.5 );
672  fOutput->Add(fhJetPt_1);
673  fhJetPt_2= new TH1F("fhJetPt_2", "Jet Pt Particle Level",1500,-0.5,149.5 );
674  fOutput->Add(fhJetPt_2);
675  fhJetPhi_1= new TH1F("fhJetPhi_1", "Jet Phi Detector Level",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
676  fOutput->Add(fhJetPhi_1);
677  fhJetPhi_2= new TH1F("fhJetPhi_2", "Jet Phi Particle Level",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
678  fOutput->Add(fhJetPhi_2);
679  fhJetEta_1= new TH1F("fhJetEta_1", "Jet Eta Detector Level", Eta_Bins, Eta_Low, Eta_Up);
680  fOutput->Add(fhJetEta_1);
681  fhJetEta_2= new TH1F("fhJetEta_2", "Jet Eta Particle Level", Eta_Bins, Eta_Low, Eta_Up);
682  fOutput->Add(fhJetEta_2);
683  fhJetMass_1= new TH1F("fhJetMass_1", "Jet Mass Detector Level", 4000,-0.5, 39.5);
684  fOutput->Add(fhJetMass_1);
685  fhJetMass_2= new TH1F("fhJetMass_2", "Jet Mass Particle Level", 4000,-0.5, 39.5);
686  fOutput->Add(fhJetMass_2);
687  fhJetRadius_1= new TH1F("fhJetRadius_1", "Jet Radius Detector Level", 100, -0.05,0.995);
688  fOutput->Add(fhJetRadius_1);
689  fhJetRadius_2= new TH1F("fhJetRadius_2", "Jet Radius Particle Level", 100, -0.05,0.995);
690  fOutput->Add(fhJetRadius_2);
691  fhNumberOfJetTracks_1= new TH1F("fhNumberOfJetTracks_1", "Number of Tracks within a Jet Detector Level", 300, -0.5,299.5);
693  fhNumberOfJetTracks_2= new TH1F("fhNumberOfJetTracks_2", "Number of Tracks within a Jet Particle Level", 300, -0.5,299.5);
695  fhSubJetRadius_1= new TH1F("fhSubJetRadius_1", "SubJet Radius Detector Level", 100, -0.05,0.995);
697  fhSubJetRadius_2= new TH1F("fhSubJetRadius_2", "SubJet Radius Particle Level", 100, -0.05,0.995);
699  fhSubJetPt_1= new TH1F("fhSubJetPt_1", "SubJet Pt Detector Level", 1500, -0.5,149.5);
700  fOutput->Add(fhSubJetPt_1);
701  fhSubJetPt_2= new TH1F("fhSubJetPt_2", "SubJet Pt Particle Level", 1500, -0.5,149.5);
702  fOutput->Add(fhSubJetPt_2);
703  fhSubJetMass_1= new TH1F("fhSubJetMass_1", "Sub Jet Mass Detector Level", 4000,-0.5, 39.5);
704  fOutput->Add(fhSubJetMass_1);
705  fhSubJetMass_2= new TH1F("fhSubJetMass_2", "Sub Jet Mass Particle Level", 4000,-0.5, 39.5);
706  fOutput->Add(fhSubJetMass_2);
707  fhNumberOfSubJetTracks_1= new TH1F("fhNumberOfSubJetTracks_1", "Number of Tracks within a Sub Jet Detector Level", 300, -0.5,299.5);
709  fhNumberOfSubJetTracks_2= new TH1F("fhNumberOfSubJetTracks_2", "Number of Tracks within a Sub Jet Particle Level", 300, -0.5,299.5);
711  fhJetCounter_1= new TH1F("fhJetCounter_1", "Jet Counter Detector Level", 150, -0.5, 149.5);
712  fOutput->Add(fhJetCounter_1);
713  fhJetCounter_2= new TH1F("fhJetCounter_2", "Jet Counter Particle Level", 150, -0.5, 149.5);
714  fOutput->Add(fhJetCounter_2);
715  fhSubJetCounter_1 = new TH1F("fhSubJetCounter_1", "SubJet Counter Detector Level",50, -0.5,49.5);
717  fhSubJetCounter_2 = new TH1F("fhSubJetCounter_2", "SubJet Counter Particle Level",50, -0.5,49.5);
719  fh2PtRatio= new TH2F("fhPtRatio", "MC pt for detector level vs particle level jets",1500,-0.5,149.5,1500,-0.5,149.5);
720  fOutput->Add(fh2PtRatio);
721  fhEventCounter_1= new TH1F("fhEventCounter_1", "Event Counter Detector Level", 15,0.5,15.5);
723  fhEventCounter_2= new TH1F("fhEventCounter_2", "Event Counter Particle Level", 15,0.5,15.5);
725  }
727  fhEventCounter= new TH1F("fhEventCounter", "Event Counter", 15,0.5,15.5);
728  fOutput->Add(fhEventCounter);
729  }
730 
731  PostData(1,fOutput);
732  PostData(2,fTreeResponseMatrixAxis);
733  // delete [] fShapesVarNames;
734 }
735 
736 //________________________________________________________________________
738 {
739  // Run analysis code here, if needed. It will be executed before FillHistograms().
740 
741 
742  return kTRUE;
743 }
744 
745 //________________________________________________________________________
747 {
748 
749  //fhEventCounter Key:
750  // 1: Number of events with a Jet Container
751  // 2: Number of Jets without a Jet Container
752  // 3:
753  // 4: Number of Jets found in all events
754  // 5: Number of Jets that were reclustered in all events
755  // 6: Number of SubJets found in all events
756  // 7: Number of events were primary vertext was found
757  // 8: Number of Jets with more than one SubJet
758  // 9:Number of Jets with more than two SubJets
759  // 12:Number of SubJetinessEvents in kData
760 
761 
762  // const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
763  // Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
764  // if(vert) fhEventCounter->Fill(7);
765 
766  //cout << ((AliVAODHeader*)(InputEvent()->GetHeader()))->GetEventplane()<< " "<<fEPV0<<endl;
767  // cout << InputEvent()->GetEventplane()->GetEventplane("Q")<< " "<<fEPV0<<endl;
768  if (fCentSelectOn){
769  if ((fCent>fCentMax) || (fCent<fCentMin)) return 0;
770  }
771 
772  AliAODTrack *TriggerHadron = 0x0;
773  if (fJetSelection == kRecoil) {
774  //you have to set a flag and the limits of the pT interval for your trigger
776  if (TriggerHadronLabel==-99999) return 0; //Trigger Hadron Not Found
777  AliTrackContainer *PartCont =NULL;
778  AliParticleContainer *PartContMC=NULL;
779  if (fJetShapeSub==kConstSub){
781  else PartCont = GetTrackContainer(1);
782  }
783  else{
785  else PartCont = GetTrackContainer(0);
786  }
787  TClonesArray *TrackArray = NULL;
788  TClonesArray *TrackArrayMC = NULL;
789  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) TrackArrayMC = PartContMC->GetArray();
790  else TrackArray = PartCont->GetArray();
791  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) TriggerHadron = static_cast<AliAODTrack*>(TrackArrayMC->At(TriggerHadronLabel));
792  else TriggerHadron = static_cast<AliAODTrack*>(TrackArray->At(TriggerHadronLabel));
793  if (!TriggerHadron) return 0;//No trigger hadron with label found
794  if(fSemigoodCorrect){
795  Double_t HoleDistance=RelativePhi(TriggerHadron->Phi(),fHolePos);
796  if(TMath::Abs(HoleDistance)+fHoleWidth+fJetRadius>TMath::Pi()-fRecoilAngularWindow) return 0;
797  }
798  fhPtTriggerHadron->Fill(TriggerHadron->Pt()); //Needed for per trigger Normalisation
799  if (fJetShapeType != AliAnalysisTaskSubJetFraction::kGenOnTheFly) fhPhiTriggerHadronEventPlane->Fill(TMath::Abs(RelativePhiEventPlane(fEPV0,TriggerHadron->Phi()))); //fEPV0 is the event plane from AliAnalysisTaskEmcal
800  if (fJetShapeType != AliAnalysisTaskSubJetFraction::kGenOnTheFly) fhPhiTriggerHadronEventPlaneTPC->Fill(TMath::Abs(RelativePhiEventPlane(((AliVAODHeader*)(InputEvent()->GetHeader()))->GetEventplane(),TriggerHadron->Phi()))); //TPC event plane
801  }
802 
803 
806  AliEmcalJet *Jet1 = NULL; //Subtracted hybrid Jet
807  AliEmcalJet *Jet2 = NULL; //Unsubtracted Hybrid Jet
808  AliEmcalJet *Jet3 = NULL; //Detector Level Pythia Jet
809  AliEmcalJet *Jet4 = NULL; //Particle Level Pyhtia Jet
810  AliJetContainer *JetCont1= GetJetContainer(0); //Jet Container for Subtracted Hybrid Jets
811  AliJetContainer *JetCont2= GetJetContainer(1); //Jet Container for Unsubtracted Hybrid Jets
812  AliJetContainer *JetCont3= GetJetContainer(2); //Jet Container for Detector Level Pyhtia Jets
813  AliJetContainer *JetCont4= GetJetContainer(3); //Jet Container for Particle Level Pythia Jets
814  AliEmcalJetFinder *Reclusterer1; //Object containg Subjets from Subtracted Hybrid Jets
815  AliEmcalJetFinder *Reclusterer4; //Object containg Subjets from Particle Level Pythia Jets
816  Bool_t JetsMatched=kFALSE;
817  Bool_t EventCounter=kFALSE;
818  Int_t JetNumber=-1;
819  Double_t JetPtThreshold=-2;
820  fhEventCounter->Fill(1);
821  if(JetCont1) {
822  fhEventCounter->Fill(2);
823  JetCont1->ResetCurrentID();
824  while((Jet1=JetCont1->GetNextAcceptJet())) {
825  if (fJetShapeSub==kConstSub) JetPtThreshold=Jet1->Pt();
826  else JetPtThreshold=Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
827  if ( (!Jet1) || (JetPtThreshold<fPtThreshold)) continue;
828  else {
829  /* if(fSemigoodCorrect){
830  Double_t HoleDistance=RelativePhi(Jet1->Phi(),fHolePos);
831  if(TMath::Abs(HoleDistance)<fHoleWidth) continue;
832  }*/
833  Float_t RecoilDeltaPhi = 0.;
834  if (fJetSelection == kRecoil){
835  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), Jet1->Phi());
836  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
837  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), Jet1->Pt());
838  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), Jet1->Phi()));
839  }
840  if (!EventCounter){
841  fhEventCounter->Fill(3);
842  EventCounter=kTRUE;
843  }
844  if(fJetShapeSub==kConstSub){
845  JetNumber=-1;
846  for(Int_t i = 0; i<JetCont2->GetNJets(); i++) {
847  Jet2 = JetCont2->GetJet(i);
848  if(Jet2->GetLabel()==Jet1->GetLabel()) {
849  JetNumber=i;
850  }
851  }
852  if(JetNumber==-1) continue;
853  Jet2=JetCont2->GetJet(JetNumber);
854  if (JetCont2->AliJetContainer::GetFractionSharedPt(Jet2)<fSharedFractionPtMin) continue;
855  Jet3=Jet2->ClosestJet();
856  }
857  if(!(fJetShapeSub==kConstSub)){
858  if (JetCont1->AliJetContainer::GetFractionSharedPt(Jet1)<fSharedFractionPtMin) continue;
859  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
860  }
861  if (!Jet3) continue;
862  Jet4=Jet3->ClosestJet();
863  if(!Jet4) continue;
864  JetsMatched=kTRUE;
866  if (fJetShapeSub==kConstSub) fShapesVar[0]=Jet1->Pt();
867  else fShapesVar[0]=Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
868  fShapesVar[2]=FjNSubJettiness(Jet1,0,1,0,1,0);
869  fShapesVar[4]=FjNSubJettiness(Jet1,0,2,0,1,0);
870  fShapesVar[6]=FjNSubJettiness(Jet1,0,2,0,1,8);
871  fShapesVar[8]=FjNSubJettiness(Jet1,0,2,0,1,1);
872  // fShapesVar[10]=Jet1->GetNumberOfTracks();
873  fShapesVar[10]=FjNSubJettiness(Jet1,0,2,0,1,9);
874  fShapesVar[12]=FjNSubJettiness(Jet1,0,2,0,1,3,fBeta_SD,fZCut);
875  fShapesVar[14]=FjNSubJettiness(Jet1,0,2,0,1,5,fBeta_SD,fZCut);
876  fShapesVar[16]=Jet1->GetLeadingTrack(JetCont1->GetParticleContainer()->GetArray())->Pt();
878  //fShapesVar[20]=RelativePhiEventPlane(((AliVAODHeader*)(InputEvent()->GetHeader()))->GetEventplane(),Jet1->Phi());
879  fShapesVar[20]=FjNSubJettiness(Jet1,0,2,0,1,6,fBeta_SD,fZCut);
880  if (fFullTree){
881  fShapesVar[22]=FjNSubJettiness(Jet1,0,2,0,1,2);
882  Reclusterer1 = Recluster(Jet1, 0, fSubJetRadius, fSubJetMinPt, fSubJetAlgorithm, "SubJetFinder_1");
883  fShapesVar[24]=SubJetFraction(Jet1, Reclusterer1, 1, 0, kTRUE, kFALSE);
884  fShapesVar[26]=SubJetFraction(Jet1, Reclusterer1, 2, 0, kTRUE, kFALSE);
885  }
886  if (JetsMatched){ //even needed? Not now but might be if you want to fill trees when jets aren't matched too
887  fShapesVar[1]=Jet4->Pt();
888  fShapesVar[3]=FjNSubJettiness(Jet4,3,1,0,1,0);
889  fShapesVar[5]=FjNSubJettiness(Jet4,3,2,0,1,0);
890  fShapesVar[7]=FjNSubJettiness(Jet4,3,2,0,1,8);
891  fShapesVar[9]=FjNSubJettiness(Jet4,3,2,0,1,1);
892  //fShapesVar[11]=Jet4->GetNumberOfTracks();
893  fShapesVar[11]=FjNSubJettiness(Jet4,3,2,0,1,9);
894  fShapesVar[13]=FjNSubJettiness(Jet4,3,2,0,1,3,fBeta_SD,fZCut);
895  fShapesVar[15]=FjNSubJettiness(Jet4,3,2,0,1,5,fBeta_SD,fZCut);
896  fShapesVar[17]=Jet4->GetLeadingTrack(JetCont4->GetParticleContainer()->GetArray())->Pt();
897  fShapesVar[19]=RelativePhiEventPlane(fEPV0,Jet4->Phi());
898  //fShapesVar[21]=RelativePhiEventPlane(((AliVAODHeader*)(InputEvent()->GetHeader()))->GetEventplane(),Jet4->Phi());
899  fShapesVar[21]=FjNSubJettiness(Jet4,3,2,0,1,6,fBeta_SD,fZCut);
900  if (fFullTree){
901  fShapesVar[23]=FjNSubJettiness(Jet4,3,2,0,1,2);
902  Reclusterer4=Recluster(Jet4, 3, fSubJetRadius, 0, fSubJetAlgorithm, "SubJetFinder_4");
903  fShapesVar[25]=SubJetFraction(Jet4, Reclusterer4, 1, 0, kTRUE, kFALSE);
904  fShapesVar[27]=SubJetFraction(Jet4, Reclusterer4, 2, 0, kTRUE, kFALSE);
905  }
906  }
907  else{
908  fShapesVar[1]=-2;
909  fShapesVar[3]=-2;
910  fShapesVar[5]=-2;
911  fShapesVar[7]=-2;
912  fShapesVar[9]=-2;
913  fShapesVar[11]=-2;
914  fShapesVar[13]=-2;
915  fShapesVar[15]=-2;
916  fShapesVar[17]=-2;
917  fShapesVar[19]=-2;
918  fShapesVar[21]=-2;
919  if (fFullTree){
920  fShapesVar[23]=-2;
921  fShapesVar[25]=-2;
922  fShapesVar[27]=-2;
923  }
924  }
925  fTreeResponseMatrixAxis->Fill();
926  JetsMatched=kFALSE;
927  }
928  }
929  }
930  }
931 
934  AliEmcalJet *Jet1 = NULL; //Detector Level Jet
935  AliEmcalJet *Jet2 = NULL; //Particle Level Jet
936  AliJetContainer *JetCont1= GetJetContainer(0); //Jet Container for Detector Level Pythia
937  AliJetContainer *JetCont2= GetJetContainer(1); //Jet Container for Particle Level Pythia
938  AliEmcalJetFinder *Reclusterer1; //Object containg Subjets Detector Level
939  AliEmcalJetFinder *Reclusterer2; //Object containg Subjets Particle Level
940  Int_t JetCounter1=0; //Counts number of jets in event
941  Int_t JetCounter2=0; //Counts number of jets in event
942  Double_t JetPhi1=0;
943  Double_t JetPhi2=0;
944  Bool_t JetsMatched=kFALSE;
945  Double_t Pythia_Event_Weight=1;
946  Bool_t EventCounter=kFALSE;
947  fhEventCounter_1->Fill(1);
948  if(JetCont1) {
949  fhEventCounter_1->Fill(2); //Number of events with a jet container
950  JetCont1->ResetCurrentID();
951  while((Jet1=JetCont1->GetNextAcceptJet())) {
952  if( (!Jet1) || ((Jet1->Pt())<fPtThreshold)) {
953  // fhEventCounter_1->Fill(3); //events where the jet had a problem
954  continue;
955  }
956  else {
957  /* if(fSemigoodCorrect){
958  Double_t HoleDistance=RelativePhi(Jet1->Phi(),fHolePos);
959  if(TMath::Abs(HoleDistance)<fHoleWidth) continue;
960  }*/
961  Float_t RecoilDeltaPhi = 0.;
962  if (fJetSelection == kRecoil){
963  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), Jet1->Phi());
964  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
965  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), Jet1->Pt());
966  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), Jet1->Phi()));
967  }
968  if (!EventCounter){
969  fhEventCounter_1->Fill(3);
970  EventCounter=kTRUE;
971  }
972  if((Jet1->GetNumberOfTracks())==0){
973  fhEventCounter_1->Fill(10); //zero track jets
974  }
975  if((Jet1->GetNumberOfTracks())==1){
976  fhEventCounter_1->Fill(11); //one track jets
977  }
978  fhEventCounter_1->Fill(4); //Number of Jets found in all events
979  JetCounter1++;
980  fhJetPt_1->Fill(Jet1->Pt());
981  JetPhi1=Jet1->Phi();
982  if(JetPhi1 < -1*TMath::Pi()) JetPhi1 += (2*TMath::Pi());
983  else if (JetPhi1 > TMath::Pi()) JetPhi1 -= (2*TMath::Pi());
984  fhJetPhi_1->Fill(JetPhi1);
985  fhJetEta_1->Fill(Jet1->Eta());
986  fhJetMass_1->Fill(Jet1->M());
987  fhJetRadius_1->Fill(TMath::Sqrt((Jet1->Area()/TMath::Pi()))); //Radius of Jets per event
989  if((Jet2 = Jet1->ClosestJet())){
990  JetsMatched=kTRUE;
991  if((Jet2->GetNumberOfTracks())==0){
992  fhEventCounter_2->Fill(10); //zero track jets
993  }
994  if((Jet2->GetNumberOfTracks())==1){
995  fhEventCounter_2->Fill(11); //one track jets
996  }
997  fhEventCounter_2->Fill(4); //Number of Jets found in all events
998  JetCounter2++;
999  fhJetPt_2->Fill(Jet2->Pt());
1000  JetPhi2=Jet2->Phi();
1001  if(JetPhi2 < -1*TMath::Pi()) JetPhi2 += (2*TMath::Pi());
1002  else if (JetPhi2 > TMath::Pi()) JetPhi2 -= (2*TMath::Pi());
1003  fhJetPhi_2->Fill(JetPhi2);
1004  fhJetEta_2->Fill(Jet2->Eta());
1005  fhJetMass_2->Fill(Jet2->M());
1006  fhJetRadius_2->Fill(TMath::Sqrt((Jet2->Area()/TMath::Pi()))); //Radius of Jets per event
1008  fh2PtRatio->Fill(Jet1->Pt(),Jet2->Pt(),Pythia_Event_Weight);
1009  }
1010  else {
1011  fhEventCounter_2->Fill(1);
1012  //continue;
1013  }
1015 
1016 
1017  fShapesVar[0]=Jet1->Pt();
1018  fShapesVar[2]=FjNSubJettiness(Jet1,0,1,0,1,0);
1019  fShapesVar[4]=FjNSubJettiness(Jet1,0,2,0,1,0);
1020  fShapesVar[6]=FjNSubJettiness(Jet1,0,2,0,1,8);
1021  fShapesVar[8]=FjNSubJettiness(Jet1,0,2,0,1,1);
1022  // fShapesVar[10]=Jet1->GetNumberOfTracks();
1023  fShapesVar[10]=FjNSubJettiness(Jet1,0,2,0,1,9);
1024  fShapesVar[12]=FjNSubJettiness(Jet1,0,2,0,1,3,fBeta_SD,fZCut);
1025  fShapesVar[14]=FjNSubJettiness(Jet1,0,2,0,1,5,fBeta_SD,fZCut);
1026  fShapesVar[16]=Jet1->GetLeadingTrack(JetCont1->GetParticleContainer()->GetArray())->Pt();
1027  fShapesVar[18]=-2; //event plane calculation only needed for PbPb recoils
1028  fShapesVar[20]=FjNSubJettiness(Jet1,0,2,0,1,6,fBeta_SD,fZCut);
1029  Reclusterer1 = Recluster(Jet1, 0, fSubJetRadius, fSubJetMinPt, fSubJetAlgorithm, "SubJetFinder_1");
1030  if (fFullTree){
1031  fShapesVar[22]=FjNSubJettiness(Jet1,0,2,0,1,2);
1032  fShapesVar[24]=SubJetFraction(Jet1, Reclusterer1, 1, 0, kTRUE, kFALSE);
1033  fShapesVar[26]=SubJetFraction(Jet1, Reclusterer1, 2, 0, kTRUE, kFALSE);
1034  }
1035  if (JetsMatched){ //even needed? Not now but might be if you want to fill trees when jets aren't matched too
1036  fShapesVar[1]=Jet2->Pt();
1037  fShapesVar[3]=FjNSubJettiness(Jet2,1,1,0,1,0);
1038  fShapesVar[5]=FjNSubJettiness(Jet2,1,2,0,1,0);
1039  fShapesVar[7]=FjNSubJettiness(Jet2,1,2,0,1,8);
1040  fShapesVar[9]=FjNSubJettiness(Jet2,1,2,0,1,1);
1041  // fShapesVar[11]=Jet2->GetNumberOfTracks();
1042  fShapesVar[11]=FjNSubJettiness(Jet2,1,2,0,1,9);
1043  fShapesVar[13]=FjNSubJettiness(Jet2,1,2,0,1,3,fBeta_SD,fZCut);
1044  fShapesVar[15]=FjNSubJettiness(Jet2,1,2,0,1,5,fBeta_SD,fZCut);
1045  fShapesVar[17]=Jet2->GetLeadingTrack(JetCont2->GetParticleContainer()->GetArray())->Pt();
1046  fShapesVar[19]=-2;
1047  fShapesVar[21]=FjNSubJettiness(Jet2,1,2,0,1,6,fBeta_SD,fZCut);
1048  Reclusterer2 = Recluster(Jet2, 1, fSubJetRadius, 0, fSubJetAlgorithm, "SubJetFinder_2");
1049  if (fFullTree){
1050  fShapesVar[23]=FjNSubJettiness(Jet2,1,2,0,1,2);
1051  fShapesVar[25]=SubJetFraction(Jet2, Reclusterer2, 1, 0, kTRUE, kFALSE);
1052  fShapesVar[27]=SubJetFraction(Jet2, Reclusterer2, 2, 0, kTRUE, kFALSE);
1053  }
1054  }
1055  else{
1056  fShapesVar[1]=-2;
1057  fShapesVar[3]=-2;
1058  fShapesVar[5]=-2;
1059  fShapesVar[7]=-2;
1060  fShapesVar[9]=-2;
1061  fShapesVar[11]=-2;
1062  fShapesVar[13]=-2;
1063  fShapesVar[15]=-2;
1064  fShapesVar[17]=-2;
1065  fShapesVar[19]=-2;
1066  fShapesVar[21]=-2;
1067  if (fFullTree){
1068  fShapesVar[23]=-2;
1069  fShapesVar[25]=-2;
1070  fShapesVar[27]=-2;
1071  }
1072  }
1073  fTreeResponseMatrixAxis->Fill();
1074 
1075  fhSubJetCounter_1->Fill(Reclusterer1->GetNumberOfJets());
1076  for (Int_t i= 0; i<Reclusterer1->GetNumberOfJets(); i++){
1077  fhEventCounter_1->Fill(6); //Number of overall subjets in all events
1078  fhSubJetPt_1->Fill(Reclusterer1->GetJet(i)->Pt());
1079  fhSubJetMass_1->Fill(Reclusterer1->GetJet(i)->M());
1080  fhNumberOfSubJetTracks_1->Fill(Reclusterer1->GetJet(i)->GetNumberOfTracks());
1081  fhSubJetRadius_1->Fill(TMath::Sqrt((Reclusterer1->GetJet(i)->Area()/TMath::Pi()))); //Radius of SubJets per event
1082  }
1083  if(JetsMatched){
1084  fhSubJetCounter_2->Fill(Reclusterer2->GetNumberOfJets());
1085  for (Int_t i= 0; i<Reclusterer2->GetNumberOfJets(); i++){
1086  fhEventCounter_2->Fill(6); //Number of overall subjets in all events
1087  fhSubJetPt_2->Fill(Reclusterer2->GetJet(i)->Pt());
1088  fhSubJetMass_2->Fill(Reclusterer2->GetJet(i)->M());
1089  fhNumberOfSubJetTracks_2->Fill(Reclusterer2->GetJet(i)->GetNumberOfTracks());
1090  fhSubJetRadius_2->Fill(TMath::Sqrt((Reclusterer2->GetJet(i)->Area()/TMath::Pi()))); //Radius of SubJets per event
1091  }
1092  }
1093  JetsMatched=kFALSE;
1094  }
1095  }
1096  fhJetCounter_1->Fill(JetCounter1); //Number of Jets in Each Event Particle Level
1097  fhJetCounter_2->Fill(JetCounter2); //Number of Jets in Each Event Detector Level
1098  }
1099  //else {fhEventCounter_1->Fill(3);} //Events with no jet container
1100  }
1101 
1102 
1103 
1105  AliEmcalJet *Jet1 = NULL; //Original Jet in the event
1106  AliJetContainer *JetCont= GetJetContainer(0); //Jet Container for event
1107  Int_t JetCounter=0; //Counts number of jets in event
1108  Double_t JetPhi=0;
1109  Bool_t EventCounter=kFALSE;
1110  Double_t JetPt_ForThreshold=0;
1111  fhEventCounter->Fill(1);
1112  if(JetCont) {
1113  fhEventCounter->Fill(2); //Number of events with a jet container
1114  JetCont->ResetCurrentID();
1115  while((Jet1=JetCont->GetNextAcceptJet())) {
1116  if(!Jet1) continue;
1117  if (fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) JetPt_ForThreshold = Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
1118  else JetPt_ForThreshold = Jet1->Pt();
1119  if(JetPt_ForThreshold<fPtThreshold) {
1120  //fhEventCounter->Fill(3); //events where the jet had a problem
1121  continue;
1122  }
1123  else {
1124  /* if(fSemigoodCorrect){
1125  Double_t HoleDistance=RelativePhi(Jet1->Phi(),fHolePos);
1126  if(TMath::Abs(HoleDistance)<fHoleWidth) continue;
1127  }*/
1128  Float_t RecoilDeltaPhi = 0.;
1129  if (fJetSelection == kRecoil){
1130  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), Jet1->Phi());
1131  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
1132  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), Jet1->Pt());
1133  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), Jet1->Phi()));
1134  }
1135  if (!EventCounter){
1136  fhEventCounter->Fill(3);
1137  EventCounter=kTRUE;
1138  }
1139  if((Jet1->GetNumberOfTracks())==0){
1140  fhEventCounter->Fill(10); //zero track jets
1141  }
1142  if((Jet1->GetNumberOfTracks())==1){
1143  fhEventCounter->Fill(11); //one track jets
1144  }
1145  fhEventCounter->Fill(4); //Number of Jets found in all events
1146  JetCounter++;
1147  fhJetPt->Fill(Jet1->Pt());
1148  JetPhi=Jet1->Phi();
1149  //if(JetPhi < -1*TMath::Pi()) JetPhi += (2*TMath::Pi());
1150  //else if (JetPhi > TMath::Pi()) JetPhi -= (2*TMath::Pi());
1151  fhJetPhi->Fill(JetPhi);
1152  fhJetEta->Fill(Jet1->Eta());
1153  fhJetMass->Fill(Jet1->M());
1154  fhJetRadius->Fill(TMath::Sqrt((Jet1->Area()/TMath::Pi()))); //Radius of Jets per event
1155  fhNumberOfJetTracks->Fill(Jet1->GetNumberOfTracks());
1156  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) fShapesVar[0]= Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
1157  else fShapesVar[0]=Jet1->Pt();
1158  fShapesVar[2]=FjNSubJettiness(Jet1,0,1,0,1,0);
1159  fShapesVar[4]=FjNSubJettiness(Jet1,0,2,0,1,0);
1160  fShapesVar[6]=FjNSubJettiness(Jet1,0,2,0,1,8);
1161  fShapesVar[8]=FjNSubJettiness(Jet1,0,2,0,1,1);
1162  //fShapesVar[10]=Jet1->GetNumberOfTracks();
1163  fShapesVar[10]=FjNSubJettiness(Jet1,0,2,0,1,9);
1164  fShapesVar[12]=FjNSubJettiness(Jet1,0,2,0,1,3,fBeta_SD,fZCut);
1165  fShapesVar[14]=FjNSubJettiness(Jet1,0,2,0,1,5,fBeta_SD,fZCut);
1166  fShapesVar[16]=Jet1->GetLeadingTrack(JetCont->GetParticleContainer()->GetArray())->Pt();
1167  fShapesVar[18]=-2; //event plane calculation not needed for data
1168  fShapesVar[20]=FjNSubJettiness(Jet1,0,2,0,1,6,fBeta_SD,fZCut);
1169  AliEmcalJetFinder *Reclusterer1 = Recluster(Jet1, 0, fSubJetRadius, fSubJetMinPt, fSubJetAlgorithm, "SubJetFinder");
1170  if (fFullTree){
1171  fShapesVar[22]=FjNSubJettiness(Jet1,0,2,0,1,2);
1172  fShapesVar[24]=SubJetFraction(Jet1, Reclusterer1, 1, 0, kTRUE, kFALSE);
1173  fShapesVar[26]=SubJetFraction(Jet1, Reclusterer1, 2, 0, kTRUE, kFALSE);
1174  }
1175  fShapesVar[1]=-2;
1176  fShapesVar[3]=-2;
1177  fShapesVar[5]=-2;
1178  fShapesVar[7]=-2;
1179  fShapesVar[9]=-2;
1180  fShapesVar[11]=-2;
1181  fShapesVar[13]=-2;
1182  fShapesVar[15]=-2;
1183  fShapesVar[17]=-2;
1184  fShapesVar[19]=-2;
1185  fShapesVar[21]=-2;
1186  if (fFullTree){
1187  fShapesVar[23]=-2;
1188  fShapesVar[25]=-2;
1189  fShapesVar[27]=-2;
1190  }
1191  fTreeResponseMatrixAxis->Fill();
1192  fhSubJetCounter->Fill(Reclusterer1->GetNumberOfJets());
1193  for (Int_t i= 0; i<Reclusterer1->GetNumberOfJets(); i++){
1194  fhEventCounter->Fill(6); //Number of overall subjets in all events
1195  fhSubJetPt->Fill(Reclusterer1->GetJet(i)->Pt());
1196  fhSubJetMass->Fill(Reclusterer1->GetJet(i)->M());
1197  fhNumberOfSubJetTracks->Fill(Reclusterer1->GetJet(i)->GetNumberOfTracks());
1198  fhSubJetRadius->Fill(TMath::Sqrt((Reclusterer1->GetJet(i)->Area()/TMath::Pi()))); //Radius of SubJets per event
1199  }
1200  }
1201  }
1202  fhJetCounter->Fill(JetCounter); //Number of Jets in Each Event
1203  }
1204  //else {fhEventCounter->Fill(2);} //Events with no jet container
1205  }
1206 
1207 
1209 
1210  AliEmcalJet *Jet1 = NULL; //Original Jet in the event
1211  AliJetContainer *JetCont= GetJetContainer(0); //Jet Container for event
1212  Int_t JetCounter=0; //Counts number of jets in event
1213  Double_t JetPhi=0;
1214  Bool_t EventCounter=kFALSE;
1215  Double_t JetPt_ForThreshold=0;
1216  fhEventCounter->Fill(1);
1217  if(JetCont) {
1218  fhEventCounter->Fill(2); //Number of events with a jet container
1219  JetCont->ResetCurrentID();
1220  while((Jet1=JetCont->GetNextAcceptJet())) {
1221  // if (((Jet1=JetCont->GetLeadingJet("rho")) && (fPtThreshold<=Jet1->Pt()-(GetRhoVal(0)*Jet1->Area())))){
1222  if(!Jet1) continue;
1223  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) JetPt_ForThreshold = Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
1224  else JetPt_ForThreshold = Jet1->Pt();
1225  if(JetPt_ForThreshold<fPtThreshold) {
1226  //fhEventCounter->Fill(3); //events where the jet had a problem
1227  continue;
1228  }
1229  else {
1230  /* if(fSemigoodCorrect){
1231  Double_t HoleDistance=RelativePhi(Jet1->Phi(),fHolePos);
1232  if(TMath::Abs(HoleDistance)<fHoleWidth) continue;
1233  }*/
1234  Float_t RecoilDeltaPhi = 0.;
1235  if (fJetSelection == kRecoil){
1236  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), Jet1->Phi());
1237  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
1238  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), Jet1->Pt());
1239  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), Jet1->Phi()));
1240  }
1241  if (!EventCounter){
1242  fhEventCounter->Fill(3);
1243  EventCounter=kTRUE;
1244  }
1245  if((Jet1->GetNumberOfTracks())==0){
1246  fhEventCounter->Fill(10); //zero track jets
1247  }
1248  if((Jet1->GetNumberOfTracks())==1){
1249  fhEventCounter->Fill(11); //one track jets
1250  }
1251  fhEventCounter->Fill(4); //Number of Jets found in all events
1252  JetCounter++;
1253  fhJetPt->Fill(Jet1->Pt());
1254  JetPhi=Jet1->Phi();
1255  if(JetPhi < -1*TMath::Pi()) JetPhi += (2*TMath::Pi());
1256  else if (JetPhi > TMath::Pi()) JetPhi -= (2*TMath::Pi());
1257  fhJetPhi->Fill(JetPhi);
1258  fhJetEta->Fill(Jet1->Eta());
1259  fhJetMass->Fill(Jet1->M());
1260  fhJetRadius->Fill(TMath::Sqrt((Jet1->Area()/TMath::Pi()))); //Radius of Jets per event
1261  fhNumberOfJetTracks->Fill(Jet1->GetNumberOfTracks());
1262  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> Randomised_Jet=ModifyJet(Jet1,0,"Randomise");
1263  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_02_15=ModifyJet(Jet1,0,"AddExtraProng_02_15");
1264  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_02_30=ModifyJet(Jet1,0,"AddExtraProng_02_30");
1265  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_02_45=ModifyJet(Jet1,0,"AddExtraProng_02_45");
1266  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_02_60=ModifyJet(Jet1,0,"AddExtraProng_02_60");
1267  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_01_30=ModifyJet(Jet1,0,"AddExtraProng_01_30");
1268  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_01_45=ModifyJet(Jet1,0,"AddExtraProng_01_45");
1269  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_03_30=ModifyJet(Jet1,0,"AddExtraProng_03_30");
1270  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_03_45=ModifyJet(Jet1,0,"AddExtraProng_03_45");
1271  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) fShapesVar[0]= Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
1272  else fShapesVar[0]=Jet1->Pt();
1273  fShapesVar[2]=FjNSubJettiness(Jet1,0,1,0,1,0);
1274  fShapesVar[4]=FjNSubJettiness(Jet1,0,2,0,1,0);
1275  fShapesVar[6]=FjNSubJettinessFastJet(ExtraProng_Jet_02_45,0,1,0,1,0);
1276  fShapesVar[8]=FjNSubJettiness(Jet1,0,2,0,1,1);
1277  // fShapesVar[10]=Jet1->GetNumberOfTracks();
1278  fShapesVar[10]=FjNSubJettinessFastJet(ExtraProng_Jet_01_45,0,1,0,1,0);
1279  fShapesVar[12]=FjNSubJettinessFastJet(ExtraProng_Jet_03_45,0,1,0,1,0);
1280  fShapesVar[14]=FjNSubJettinessFastJet(ExtraProng_Jet_02_15,0,2,0,1,0);
1281  fShapesVar[16]=FjNSubJettinessFastJet(ExtraProng_Jet_02_45,0,2,0,1,0);
1282  fShapesVar[18]=FjNSubJettinessFastJet(ExtraProng_Jet_01_30,0,2,0,1,0);
1283  fShapesVar[20]=FjNSubJettinessFastJet(ExtraProng_Jet_03_30,0,2,0,1,0);
1284  AliEmcalJetFinder *Reclusterer1 = Recluster(Jet1, 0, fSubJetRadius, fSubJetMinPt, fSubJetAlgorithm, "SubJetFinder");
1285  if (fFullTree){
1286  fShapesVar[22]=FjNSubJettiness(Jet1,0,2,0,1,2);
1287  fShapesVar[24]=SubJetFraction(Jet1, Reclusterer1, 1, 0, kTRUE, kFALSE);
1288  fShapesVar[26]=SubJetFraction(Jet1, Reclusterer1, 2, 0, kTRUE, kFALSE);
1289  }
1290  fShapesVar[1]=FjNSubJettinessFastJet(Randomised_Jet,0,1,0,1,0);
1291  // fShapesVar[3]=FjNSubJettiness(std::unique_ptr<AliEmcalJet>(ModifyJet(Jet1,0,"Randomise")).get(),0,1,0,1,0);
1292  //fShapesVar[5]=FjNSubJettiness(std::unique_ptr<AliEmcalJet>(ModifyJet(Jet1,0,"AddExtraProng_02_30")).get(),0,1,0,1,0);
1293  fShapesVar[3]=FjNSubJettinessFastJet(ExtraProng_Jet_02_15,0,1,0,1,0);
1294  fShapesVar[5]=FjNSubJettinessFastJet(ExtraProng_Jet_02_30,0,1,0,1,0);
1295  fShapesVar[7]=FjNSubJettinessFastJet(ExtraProng_Jet_02_60,0,1,0,1,0);
1296  fShapesVar[9]=FjNSubJettinessFastJet(ExtraProng_Jet_01_30,0,1,0,1,0);
1297  fShapesVar[11]=FjNSubJettinessFastJet(ExtraProng_Jet_03_30,0,1,0,1,0);
1298  fShapesVar[13]=FjNSubJettinessFastJet(Randomised_Jet,0,2,0,1,0);
1299  fShapesVar[15]=FjNSubJettinessFastJet(ExtraProng_Jet_02_30,0,2,0,1,0);
1300  fShapesVar[17]=FjNSubJettinessFastJet(ExtraProng_Jet_02_60,0,2,0,1,0);
1301  fShapesVar[19]=FjNSubJettinessFastJet(ExtraProng_Jet_01_45,0,2,0,1,0);
1302  fShapesVar[21]=FjNSubJettinessFastJet(ExtraProng_Jet_03_45,0,2,0,1,0);
1303  if (fFullTree){
1304  fShapesVar[23]=-2;
1305  fShapesVar[25]=-2;
1306  fShapesVar[27]=-2;
1307  }
1308  fTreeResponseMatrixAxis->Fill();
1309 
1310  fhSubJetCounter->Fill(Reclusterer1->GetNumberOfJets());
1311  for (Int_t i= 0; i<Reclusterer1->GetNumberOfJets(); i++){
1312  fhEventCounter->Fill(6); //Number of overall subjets in all events
1313  fhSubJetPt->Fill(Reclusterer1->GetJet(i)->Pt());
1314  fhSubJetMass->Fill(Reclusterer1->GetJet(i)->M());
1315  fhNumberOfSubJetTracks->Fill(Reclusterer1->GetJet(i)->GetNumberOfTracks());
1316  fhSubJetRadius->Fill(TMath::Sqrt((Reclusterer1->GetJet(i)->Area()/TMath::Pi()))); //Radius of SubJets per event
1317  }
1318  delete Randomised_Jet.second;
1319  delete ExtraProng_Jet_02_15.second;
1320  delete ExtraProng_Jet_02_30.second;
1321  delete ExtraProng_Jet_02_45.second;
1322  delete ExtraProng_Jet_02_60.second;
1323  delete ExtraProng_Jet_01_30.second;
1324  delete ExtraProng_Jet_01_45.second;
1325  delete ExtraProng_Jet_03_30.second;
1326  delete ExtraProng_Jet_03_45.second;
1327  }
1328  }
1329  fhJetCounter->Fill(JetCounter); //Number of Jets in Each Event
1330  }
1331  //else {fhEventCounter->Fill(2);} //Events with no jet container
1332  }
1333  return kTRUE;
1334 }
1335 //________________________________________________________________________
1337 
1338  if(Phi < -1*TMath::Pi()) Phi += (2*TMath::Pi());
1339  else if (Phi > TMath::Pi()) Phi -= (2*TMath::Pi());
1340  Double_t DeltaPhi=Phi-EventPlane;
1341  if(DeltaPhi < -1*TMath::Pi()) DeltaPhi += (2*TMath::Pi());
1342  else if (DeltaPhi > TMath::Pi()) DeltaPhi -= (2*TMath::Pi());
1343  return DeltaPhi;
1344 }
1345 //________________________________________________________________________
1347 
1348  if(Phi1 < -1*TMath::Pi()) Phi1 += (2*TMath::Pi()); // Turns the range of 0to2Pi into -PitoPi ???????????
1349  else if (Phi1 > TMath::Pi()) Phi1 -= (2*TMath::Pi());
1350  if(Phi2 < -1*TMath::Pi()) Phi2 += (2*TMath::Pi());
1351  else if (Phi2 > TMath::Pi()) Phi2 -= (2*TMath::Pi());
1352  Double_t DeltaPhi=Phi2-Phi1;
1353  if(DeltaPhi < -1*TMath::Pi()) DeltaPhi += (2*TMath::Pi());
1354  else if (DeltaPhi > TMath::Pi()) DeltaPhi -= (2*TMath::Pi());
1355  return DeltaPhi;
1356 }
1357 
1358 
1359 //--------------------------------------------------------------------------
1361 
1362  AliTrackContainer *PartCont = NULL;
1363  AliParticleContainer *PartContMC = NULL;
1364 
1365 
1366  if (fJetShapeSub==kConstSub){
1368  else PartCont = GetTrackContainer(1);
1369  }
1370  else{
1372  else PartCont = GetTrackContainer(0);
1373  }
1374 
1375  TClonesArray *TracksArray = NULL;
1376  TClonesArray *TracksArrayMC = NULL;
1377  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) TracksArrayMC = PartContMC->GetArray();
1378  else TracksArray = PartCont->GetArray();
1379 
1381  if(!PartContMC || !TracksArrayMC) return -99999;
1382  }
1383  else {
1384  if(!PartCont || !TracksArray) return -99999;
1385  }
1386 
1387  AliAODTrack *Track = 0x0;
1388  Int_t Trigger_Index[100];
1389  for (Int_t i=0; i<100; i++) Trigger_Index[i] = 0;
1390  Int_t Trigger_Counter = 0;
1391  Int_t NTracks=0;
1392  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) NTracks = TracksArrayMC->GetEntriesFast();
1393  else NTracks = TracksArray->GetEntriesFast();
1394  for(Int_t i=0; i < NTracks; i++){
1396  if((Track = static_cast<AliAODTrack*>(PartContMC->GetAcceptParticle(i)))){
1397  if (!Track) continue;
1398  if(TMath::Abs(Track->Eta())>0.9) continue;
1399  if (Track->Pt()<0.15) continue;
1400  if ((Track->Pt() >= PtMin) && (Track->Pt()< PtMax)) {
1401  Trigger_Index[Trigger_Counter] = i;
1402  Trigger_Counter++;
1403  }
1404  }
1405  }
1406  else{
1407  if((Track = static_cast<AliAODTrack*>(PartCont->GetAcceptTrack(i)))){
1408  if (!Track) continue;
1409  if(TMath::Abs(Track->Eta())>0.9) continue;
1410  if (Track->Pt()<0.15) continue;
1411  if ((Track->Pt() >= PtMin) && (Track->Pt()< PtMax)) {
1412  Trigger_Index[Trigger_Counter] = i;
1413  Trigger_Counter++;
1414  }
1415  }
1416  }
1417  }
1418  if (Trigger_Counter == 0) return -99999;
1419  Int_t RandomNumber = 0, Index = 0 ;
1420  TRandom3* Random = new TRandom3(0);
1421  RandomNumber = Random->Integer(Trigger_Counter);
1422  Index = Trigger_Index[RandomNumber];
1423  return Index;
1424 }
1425 
1426 
1427 //--------------------------------------------------------------------------
1429 
1430  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1431  Double_t Angularity_Numerator=0; //Reset these values
1432  Double_t Angularity_Denominator=0;
1433  AliVParticle *Particle=0x0;
1434  Double_t DeltaPhi=0;
1435 
1436 
1437  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks (particles in the jet
1438  Particle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
1439  if(!Particle) continue;
1440  DeltaPhi=RelativePhi(Jet->Phi(),Particle->Phi());
1441  Angularity_Numerator=Angularity_Numerator+(Particle->Pt()*TMath::Sqrt(((Particle->Eta()-Jet->Eta())*(Particle->Eta()-Jet->Eta()))+(DeltaPhi*DeltaPhi)));
1442  Angularity_Denominator= Angularity_Denominator+Particle->Pt();
1443  }
1444  if(Angularity_Denominator!=0) return Angularity_Numerator/Angularity_Denominator;
1445  else return -1;
1446 
1447 }
1448 
1449 
1450 
1451 //--------------------------------------------------------------------------
1453 
1454  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1455  Double_t PTD_Numerator=0; //Reset these values
1456  Double_t PTD_Denominator=0;
1457  AliVParticle *Particle=0x0;
1458  Double_t DeltaPhi=0;
1459  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks (particles in the jet
1460  Particle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
1461  if(!Particle) continue;
1462  PTD_Numerator=PTD_Numerator+(Particle->Pt()*Particle->Pt());
1463  PTD_Denominator=PTD_Denominator+Particle->Pt();
1464  }
1465  if(PTD_Denominator!=0) return TMath::Sqrt(PTD_Numerator)/PTD_Denominator;
1466  else return -1;
1467 
1468 }
1469 
1470 
1471 //----------------------------------------------------------------------
1473 AliEmcalJetFinder *AliAnalysisTaskSubJetFraction::Recluster(AliEmcalJet *Jet, Int_t JetContNb, Double_t SubJetRadius, Double_t SubJetMinPt, Int_t Algorithm, const char* Name){
1474 
1475  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1476  AliEmcalJetFinder *Reclusterer = new AliEmcalJetFinder(Name); //JetFinder Object for reclustered jets
1477  Reclusterer->SetRadius(SubJetRadius);
1478  Reclusterer->SetJetMinPt(SubJetMinPt);
1479  Reclusterer->SetJetAlgorithm(Algorithm); //0 for anti-kt 1 for kt
1480  Reclusterer->SetJetMaxEta(0.9);
1481  Reclusterer->SetRecombSheme(0);
1483  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
1484  Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
1485  if(Reclusterer->AliEmcalJetFinder::Filter(Jet, JetCont, dVtx)){;} //reclustering jet1 using the jetfinderobject Reclusterer
1486  }
1487  else{
1488  Double_t dVtx[3]={1,1,1};
1489  if(Reclusterer->AliEmcalJetFinder::Filter(Jet, JetCont, dVtx)){;} //reclustering jet1 using the jetfinderobject Reclusterer
1490  }
1491  return Reclusterer;
1492 
1493 }
1494 
1495 
1496 
1497 //----------------------------------------------------------------------
1499  AliEmcalJet *SubJet=NULL;
1500  Double_t SortingVariable;
1501  Int_t ArraySize =N+1;
1502  TArrayD JetSorter(ArraySize);
1503  TArrayD JetIndexSorter(ArraySize);
1504  for (Int_t i=0; i<ArraySize; i++){
1505  JetSorter[i]=0;
1506  }
1507  for (Int_t i=0; i<ArraySize; i++){
1508  JetIndexSorter[i]=0;
1509  }
1510  if(Reclusterer->GetNumberOfJets()<N) return -999;
1511  for (Int_t i=0; i<Reclusterer->GetNumberOfJets(); i++){
1512  SubJet=Reclusterer->GetJet(i);
1513  if (Type==0) SortingVariable=SubJet->Pt();
1514  else if (Type==1) SortingVariable=SubJet->E();
1515  else if (Type==2) SortingVariable=SubJet->M();
1516  for (Int_t j=0; j<N; j++){
1517  if (SortingVariable>JetSorter[j]){
1518  for (Int_t k=N-1; k>=j; k--){
1519  JetSorter[k+1]=JetSorter[k];
1520  JetIndexSorter[k+1]=JetIndexSorter[k];
1521  }
1522  JetSorter[j]=SortingVariable;
1523  JetIndexSorter[j]=i;
1524  break;
1525  }
1526  }
1527  }
1528  if (!Index) return JetSorter[N-1];
1529  else return JetIndexSorter[N-1];
1530 }
1531 
1532 
1533 
1534 //returns -1 if the Nth hardest jet is requested where N>number of available jets
1535 //type: 0=Pt 1=E 2=M
1536 //Index TRUE=returns index FALSE=returns value of quantatiy in question
1537 
1538 
1539 
1540 
1541 
1542 
1543 
1544 //----------------------------------------------------------------------
1546  AliEmcalJet *SubJet=NULL;
1547  Double_t Observable=0;
1548  Double_t Fraction_Numerator=0;
1549  Bool_t Error=kFALSE;
1550  if (!Jet->GetNumberOfTracks()) return -2;
1551  if (Add){
1552  for (Int_t i=1; i<=N; i++){
1553  Observable=SubJetOrdering(Jet,Reclusterer,i,Type,kFALSE);
1554  if(Observable==-999){
1555  Error = kTRUE;
1556  return -2;
1557  }
1558  Fraction_Numerator=Fraction_Numerator+Observable;
1559  }
1560  }
1561  else {
1562  Fraction_Numerator=SubJetOrdering(Jet,Reclusterer,N,Type,kFALSE);
1563  if (Fraction_Numerator==-999) return -2;
1564  }
1565  if (Type==0){
1566  if(Loss) return (Jet->Pt()-Fraction_Numerator)/Jet->Pt();
1567  else return Fraction_Numerator/Jet->Pt();
1568  }
1569  else if (Type==1){
1570  if(Loss) return (Jet->E()-Fraction_Numerator)/Jet->E();
1571  else return Fraction_Numerator/Jet->E();
1572  }
1573  else { //change to else if if you want to add more later
1574  if(Loss) return (Jet->M()-Fraction_Numerator)/Jet->M();
1575  else return Fraction_Numerator/Jet->M();
1576  }
1577 }
1578 //N number of hardest subjets involved
1579 //Type 0=Pt 1=E 2=M
1580 // Add TRUE: Add 1 to Nth hardest subjets together/total jet False: Nth hardest Jet/toal jet
1581 //Loss TRUE: Jet-Subjet(s)/Jet FALSE: Subjet(s)/Jet
1582 
1583 
1584 //----------------------------------------------------------------------
1586  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1587  AliEmcalJet *SubJet=NULL;
1588  Double_t DeltaR1=0;
1589  Double_t DeltaR2=0;
1590  AliVParticle *JetParticle=0x0;
1591  Double_t SubJetiness_Numerator = 0;
1592  Double_t SubJetiness_Denominator = 0;
1593  Double_t Index=-2;
1594  Bool_t Error=kFALSE;
1595  // JetRadius=TMath::Sqrt((Jet->Area()/TMath::Pi())); //comment out later
1596  if (!Jet->GetNumberOfTracks()) return -2;
1597  if (Reclusterer->GetNumberOfJets() < N) return -2;
1598  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks (particles in the jet
1599  JetParticle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
1600  for (Int_t j=1; j<=N; j++){
1601  Index=SubJetOrdering(Jet,Reclusterer,j,0,kTRUE);
1602  if(Index==-999){
1603  Error = kTRUE;
1604  i=Jet->GetNumberOfTracks();
1605  break;
1606  }
1607  if(j==1){
1608  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);
1609  }
1610  else{
1611  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);
1612  if (DeltaR2<DeltaR1) DeltaR1=DeltaR2;
1613  }
1614  }
1615  SubJetiness_Numerator=SubJetiness_Numerator+(JetParticle->Pt()*DeltaR1);
1616  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));
1617  else SubJetiness_Denominator=SubJetiness_Denominator+(TMath::Power((Reclusterer->GetJet(SubJetOrdering(Jet,Reclusterer,N,0,kTRUE))->Pt()),A)*JetParticle->Pt()*TMath::Power(JetRadius,B));
1618  }
1619  if (SubJetiness_Denominator!=0 && !Error){
1620  // if (SubJetiness_Numerator/SubJetiness_Denominator>1.0) cout << JetRadius<<endl;
1621  return SubJetiness_Numerator/SubJetiness_Denominator; }
1622  else return -2;
1623 }
1624 
1625 
1626 //______________________________________________________________________________________
1628 
1629  //WARNING!!! Only works for parent jets that are clustered with Anti-Kt! To change go to AliEmcalJetFinder.cxx and look at the Nsubjettiness() function
1630 
1631  //Algorithm==0 -> kt_axes;
1632  // Algorithm==1 -> ca_axes;
1633  //Algorithm==2 -> antikt_0p2_axes;
1634  //Algorithm==3 -> wta_kt_axes;
1635  //Algorithm==4 -> wta_ca_axes;
1636  //Algorithm==5 -> onepass_kt_axes;
1637  //Algorithm==6 -> onepass_ca_axes;
1638  //Algorithm==7 -> onepass_antikt_0p2_axes;
1639  //Algorithm==8 -> onepass_wta_kt_axes;
1640  //Algorithm==9 -> onepass_wta_ca_axes;
1641  //Algorithm==10 -> min_axes;
1642 
1643 
1644  //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.
1645 
1646 
1647  //Option==0 returns Nsubjettiness Value
1648  //Option==1 && N==2 returns opening angle between two subjet axes(Delta R?)
1649  //Option==2 && N==2 returns Delta R
1650  //Option==3 returns first splitting distance for soft dropped jet (CA)
1651  //Option==4 returns Symmetry measure (Zg) for soft dropped jet (CA)
1652  //Option==5 returns Pt of Subjet1
1653  //Option==6 returns Pt of Subjet2
1654  //Options==7 trutns deltaR of subjets...Is this different to before??
1655  //Option==8 Subjet1 Leading track Pt
1656  //Option==9 Subjet1 Leading track Pt
1657 
1658 
1659 
1660  if (Jet->GetNumberOfTracks()>=N){
1661  if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==1) && (Algorithm==0) && (Beta==1.0) && (Option==0)){
1664  }
1665  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==0) && (Beta==1.0) && (Option==0)){
1668  }
1669  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==3) && (Algorithm==0) && (Beta==1.0) && (Option==0)){
1672  }
1673  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==0) && (Beta==1.0) && (Option==1)){
1676  }
1677  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==1) && (Algorithm==1) && (Beta==1.0) && (Option==0)){
1680  }
1681  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==1) && (Beta==1.0) && (Option==0)){
1684  }
1685  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==1) && (Beta==1.0) && (Option==1)){
1688  }
1689  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==1) && (Algorithm==2) && (Beta==1.0) && (Option==0)){
1692  }
1693  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==2) && (Beta==1.0) && (Option==0)){
1696  }
1697  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==2) && (Beta==1.0) && (Option==1)){
1700  }
1701  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==1) && (Algorithm==6) && (Beta==1.0) && (Option==0)){
1704  }
1705  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==6) && (Beta==1.0) && (Option==0)){
1708  }
1709  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==6) && (Beta==1.0) && (Option==1)){
1712  }
1713  else{
1714  Algorithm=fReclusteringAlgorithm;
1715  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1716  AliEmcalJetFinder JetFinder("Nsubjettiness");
1717  JetFinder.SetJetMaxEta(0.9-fJetRadius);
1718  JetFinder.SetRadius(fJetRadius);
1719  JetFinder.SetJetAlgorithm(0); //0 for anti-kt 1 for kt //this is for the JET!!!!!!!!!! Not the SubJets
1720  JetFinder.SetRecombSheme(0);
1721  JetFinder.SetJetMinPt(Jet->Pt());
1723  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
1724  Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
1725  return JetFinder.Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,fSubJetRadius,Beta,Option,0,Beta_SD,ZCut,fSoftDropOn);
1726  }
1727  else{
1728  Double_t dVtx[3]={1,1,1};
1729  if (!fNsubMeasure){
1730  return JetFinder.Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,fSubJetRadius,Beta,Option,0,Beta_SD,ZCut,fSoftDropOn);
1731  }
1732  else{
1733  if (Option==3) return JetFinder.Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,fSubJetRadius,Beta,Option,0,Beta_SD,ZCut,fSoftDropOn);
1734  else return JetFinder.Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,fSubJetRadius,Beta,Option,1);
1735  }
1736  }
1737  }
1738  }
1739  else return -2;
1740 }
1741 
1742 
1743 //________________________________________________________________________
1745  //
1746  // retrieve event objects
1747  //
1749  return kFALSE;
1750 
1751  return kTRUE;
1752 }
1753 
1754 
1755 //_______________________________________________________________________
1757 {
1758  // Called once at the end of the analysis.
1760 
1761  /*
1762  for (int i=1; i<=(XBinsJetPtSize)-1; i++){ //Rescales the fhJetPt Histograms according to the with of the variable bins and number of events
1763  fhJetPt->SetBinContent(i,((1.0/(XBinsJetPt[i]-XBinsJetPt[i-1]))*((fhJetPt->GetBinContent(i))*(1.0/(fhEventCounter->GetBinContent(1))))));
1764  fhSubJetPt->SetBinContent(i,((1.0/(XBinsJetPt[i]-XBinsJetPt[i-1]))*((fhSubJetPt->GetBinContent(i))*(1.0/(fhEventCounter->GetBinContent(5))))));
1765 
1766  //fhJetPt->SetBinContent(i,((1.0/(fhPt->GetBinWidth(i)))*((fhJetPt->GetBinContent(i))*(fhEventCounter->GetBinContent(1)))));
1767  // fhJetPt->SetBinContent(i,((1.0/((fhPt->GetLowEdge(i+1))-(fhJetPt->GetLowEdge(i))))*((fhPt->GetBinContent(i))*(fhEventCounter->GetBinContent(1)))));
1768  }
1769  for (int i=1; i<=(XBinsJetMassSize)-1; i++){ //Rescales the fhPt Histograms according to the with of the variable bins and number of events
1770  fhJetMass->SetBinContent(i,((1.0/(XBinsJetMass[i]-XBinsJetMass[i-1]))*((fhJetMass->GetBinContent(i))*(1.0/(fhEventCounter->GetBinContent(1))))));
1771  }
1772 
1773  fhJetPhi->Scale(Phi_Bins/((Phi_Up-Phi_Low)*((fhEventCounter->GetBinContent(1)))));
1774  fhJetEta->Scale(Eta_Bins/((Eta_Up-Eta_Low)*((fhEventCounter->GetBinContent(1)))));
1775  fhJetRadius->Scale(100/(fhEventCounter->GetBinContent(4))); //should this and JetAngularity be divided by Bin 1 or 4????
1776  fhNumberOfJetTracks->Scale(1.0/(fhEventCounter->GetBinContent(4)));
1777 
1778  fhJetCounter->Scale(1.0/(fhEventCounter->GetBinContent(1))); //is the first bin the correct one to look at?
1779 
1780 
1781  fhSubJetCounter->Scale(1.0/(fhEventCounter->GetBinContent(5)));
1782 
1783  */
1784  /*
1785 
1786  fhJetPt->Scale(1.0/(fhEventCounter->GetBinContent(1)));
1787  fhSubJetPt->Scale(1.0/(fhEventCounter->GetBinContent(5)));
1788  fhJetMass->Scale(1.0/(fhEventCounter->GetBinContent(1)));
1789  fhJetPhi->Scale(1.0/(fhEventCounter->GetBinContent(1)));
1790  fhJetEta->Scale(1.0/(fhEventCounter->GetBinContent(1)));
1791  fhJetRadius->Scale(1.0/(fhEventCounter->GetBinContent(4))); //should this and JetAngularity be divided by Bin 1 or 4????
1792  fhNumberOfJetTracks->Scale(1.0/(fhEventCounter->GetBinContent(4)));
1793  fhJetCounter->Scale(1.0/(fhEventCounter->GetBinContent(1))); //is the first bin the correct one to look at?
1794  fhSubJetCounter->Scale(1.0/(fhEventCounter->GetBinContent(5)));
1795  */
1796  }
1797 
1798 }
1799 
1801 std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> AliAnalysisTaskSubJetFraction::ModifyJet(AliEmcalJet* Jet, Int_t JetContNb, TString Modification){
1802  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> Jet_ClusterSequence;
1803  std::vector<fastjet::PseudoJet> fInputVectors;
1804  AliJetContainer *fJetCont = GetJetContainer(JetContNb);
1805  AliParticleContainer *fTrackCont = fJetCont->GetParticleContainer();
1806  std::vector<fastjet::PseudoJet> Modified_Jet;
1807  fastjet::ClusterSequence *fClustSeqSA;
1808  Int_t NJetTracksTest = Jet->GetNumberOfTracks();
1809  if (fTrackCont) for (Int_t i=0; i<Jet->GetNumberOfTracks(); i++) {
1810  AliVParticle *fTrk = Jet->TrackAt(i, fTrackCont->GetArray());
1811  if (!fTrk) continue;
1812  fastjet::PseudoJet PseudoTracks(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
1813  PseudoTracks.set_user_index(Jet->TrackAt(i)+100); //100 is very arbitary....why is it here anyway?
1814  fInputVectors.push_back(PseudoTracks);
1815  }
1816  if (Modification=="Randomise") fInputVectors=RandomiseTracks(Jet,fInputVectors);
1817  if (Modification=="AddExtraProng_02_15") fInputVectors=AddExtraProng(fInputVectors,0.2,0.15);
1818  if (Modification=="AddExtraProng_02_30") fInputVectors=AddExtraProng(fInputVectors,0.2,0.30);
1819  if (Modification=="AddExtraProng_02_45") fInputVectors=AddExtraProng(fInputVectors,0.2,0.45);
1820  if (Modification=="AddExtraProng_02_60") fInputVectors=AddExtraProng(fInputVectors,0.2,0.60);
1821  if (Modification=="AddExtraProng_01_30") fInputVectors=AddExtraProng(fInputVectors,0.1,0.30);
1822  if (Modification=="AddExtraProng_01_45") fInputVectors=AddExtraProng(fInputVectors,0.1,0.45);
1823  if (Modification=="AddExtraProng_03_30") fInputVectors=AddExtraProng(fInputVectors,0.3,0.30);
1824  if (Modification=="AddExtraProng_03_45") fInputVectors=AddExtraProng(fInputVectors,0.3,0.45);
1825  fJetRadius=fJetRadius*2.0;
1826  try {
1827  fastjet::JetDefinition fJetDef(fastjet::antikt_algorithm, fJetRadius, static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
1828  fClustSeqSA=new fastjet::ClusterSequence(fInputVectors, fJetDef);
1829  Modified_Jet=fClustSeqSA->inclusive_jets(0);
1830  Jet_ClusterSequence.second=fClustSeqSA;
1831  } catch (fastjet::Error) {
1832  AliError(" [w] FJ Exception caught.");
1833  //return -1;
1834  }
1835  fJetRadius=fJetRadius/2.0;
1836  Jet_ClusterSequence.first=Modified_Jet[0];
1837  return Jet_ClusterSequence;
1838 }
1839 
1840 
1841 std::vector<fastjet::PseudoJet> AliAnalysisTaskSubJetFraction::RandomiseTracks(AliEmcalJet *Jet,std::vector<fastjet::PseudoJet> fInputVectors){
1842  TRandom3 fRandom1;
1843  fRandom1.SetSeed(0);
1844  Double_t Random_Phi=0.0;
1845  Double_t Random_Phi_Change=0.0;
1846  Double_t Random_Eta=0.0;
1847  Double_t Random_Eta_Change=0.0;
1848  Double_t Track_Pt=0.0;
1849  std::vector<fastjet::PseudoJet> Random_Track_Vector;
1850  Random_Track_Vector.clear();
1851  for (Int_t i=0; i< fInputVectors.size(); i++){
1852  Random_Phi_Change=fRandom1.Uniform(-1*fJetRadius,fJetRadius);
1853  Random_Phi=Jet->Phi()+Random_Phi_Change;
1854  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)));
1855  Random_Eta=Jet->Eta()+Random_Eta_Change;
1856  Track_Pt=fInputVectors[i].perp();
1857  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))));
1858  Random_Track.set_user_index(i);
1859  Random_Track_Vector.push_back(Random_Track);
1860  }
1861  fInputVectors.clear();
1862  return Random_Track_Vector;
1863 }
1864 
1865 std::vector<fastjet::PseudoJet> AliAnalysisTaskSubJetFraction::AddExtraProng(std::vector<fastjet::PseudoJet> fInputVectors, Double_t Distance, Double_t PtFrac){
1866  TRandom3 fRandom1;
1867  fRandom1.SetSeed(0);
1868  std::vector<fastjet::PseudoJet> Jet;
1869  Jet.clear();
1870  try {
1871  fastjet::JetDefinition fJetDef(fastjet::antikt_algorithm, fJetRadius*2, static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
1872  fastjet::ClusterSequence fClustSeqSA_Prong(fInputVectors, fJetDef);
1873  Jet= fClustSeqSA_Prong.inclusive_jets(0);
1874  } catch (fastjet::Error) {
1875  AliError(" [w] FJ Exception caught.");
1876  }
1877  Double_t Extra_Track_Phi_Change=fRandom1.Uniform(-1*Distance,Distance);
1878  Double_t Extra_Track_Phi=Jet[0].phi()+Extra_Track_Phi_Change;
1879  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)));
1880  Double_t Extra_Track_Eta=Jet[0].pseudorapidity()+Extra_Track_Eta_Change;
1881  Double_t Extra_Track_Pt=Jet[0].perp()*PtFrac;
1882  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))));
1883  ExtraTrack.set_user_index(150);
1884  fInputVectors.push_back(ExtraTrack);
1885  return fInputVectors;
1886 }
1888 
1889 //______________________________________________________________________________________
1890 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){
1891 
1892  //Only Works for kGenOnTheFly
1893 
1894  //WARNING!!! Only works for parent jets that are clustered with Anti-Kt! To change go to AliEmcalJetFinder.cxx and look at the Nsubjettiness() function
1895 
1896  //Algorithm==0 -> kt_axes;
1897  // Algorithm==1 -> ca_axes;
1898  //Algorithm==2 -> antikt_0p2_axes;
1899  //Algorithm==3 -> wta_kt_axes;
1900  //Algorithm==4 -> wta_ca_axes;
1901  //Algorithm==5 -> onepass_kt_axes;
1902  //Algorithm==6 -> onepass_ca_axes;
1903  //Algorithm==7 -> onepass_antikt_0p2_axes;
1904  //Algorithm==8 -> onepass_wta_kt_axes;
1905  //Algorithm==9 -> onepass_wta_ca_axes;
1906  //Algorithm==10 -> min_axes;
1907 
1908 
1909  //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.
1910 
1911 
1912  //Option==0 returns Nsubjettiness Value
1913  //Option==1 && N==2 returns opening angle between two subjet axes(Delta R?)
1914  //Option==2 && N==2 returns Delta R
1915  //Option==3 returns first splitting distance for soft dropped jet (CA)
1916  //Option==4 returns Symmetry measure (Zg) for soft dropped jet (CA)
1917  //Option==5 returns Pt of Subjet1
1918  //Option==6 returns Pt of Subjet2
1919  //Options==7 trutns deltaR of subjets...Is this different to before??
1920  //Option==8 Subjet1 Leading track Pt
1921  //Option==9 Subjet1 Leading track Pt
1922 
1923  fastjet::PseudoJet Jet=Jet_ClusterSequence.first;
1924  AliFJWrapper FJWrapper("FJWrapper", "FJWrapper");
1925  FJWrapper.SetR(fJetRadius);
1926  FJWrapper.SetMinJetPt(Jet.perp());
1927  FJWrapper.SetAlgorithm(fastjet::antikt_algorithm); //this is for the jet clustering not the subjet reclustering.
1928  FJWrapper.SetRecombScheme(static_cast<fastjet::RecombinationScheme>(0));
1929  if (Jet.constituents().size()>=N){
1930  Algorithm=fReclusteringAlgorithm;
1931  Double_t dVtx[3]={1,1,1};
1932  return FJWrapper.NSubjettinessDerivativeSub(N,Algorithm,fSubJetRadius,Beta,fJetRadius,Jet,Option,0,Beta_SD,ZCut,fSoftDropOn); //change this
1933  }
1934  else return -2;
1935 }
1936 
void SetRadius(Double_t val)
Bool_t fCentSelectOn
Random number generator.
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 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 GetFirstOrderSubtractedOpeningAngle_onepassca() const
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)
Double_t GetSecondOrderSubtracted1subjettiness_onepassca() const
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 GetFirstOrderSubtracted2subjettiness_onepassca() 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 GetSecondOrderSubtractedOpeningAngle_onepassca() 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 GetSecondOrderSubtracted2subjettiness_onepassca() const
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)
Double_t GetFirstOrderSubtracted1subjettiness_onepassca() const
const Int_t nVar
AliTrackContainer * GetTrackContainer(Int_t i=0) const
void SetMakeGeneralHistograms(Bool_t g)
virtual AliVTrack * GetAcceptTrack(Int_t i=-1) const
Base task in the EMCAL jet framework.
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 GetSecondOrderSubtractedOpeningAngle_ca() const
Double_t GetSecondOrderSubtractedOpeningAngle_kt() const
AliEmcalJet * GetJet(Int_t i) const