AliPhysics  d565ceb (d565ceb)
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  fhTrackPhi(0x0),
149  fhTrackPhi_Cut(0x0),
150  fh2PtRatio(0x0),
151  fhEventCounter(0x0),
152  fhEventCounter_1(0x0),
153  fhEventCounter_2(0x0),
154  fhSubJettiness1CheckRatio_FJ_AKT(0x0),
155  fhSubJettiness1CheckRatio_FJ_KT(0x0),
156  fhSubJettiness1CheckRatio_FJ_CA(0x0),
157  fhSubJettiness1CheckRatio_FJ_WTA_KT(0x0),
158  fhSubJettiness1CheckRatio_FJ_WTA_CA(0x0),
159  fhSubJettiness1CheckRatio_FJ_OP_AKT(0x0),
160  fhSubJettiness1CheckRatio_FJ_OP_KT(0x0),
161  fhSubJettiness1CheckRatio_FJ_OP_CA(0x0),
162  fhSubJettiness1CheckRatio_FJ_OP_WTA_KT(0x0),
163  fhSubJettiness1CheckRatio_FJ_OP_WTA_CA(0x0),
164  fhSubJettiness1CheckRatio_FJ_MIN(0x0),
165  fhSubJettiness2CheckRatio_FJ_AKT(0x0),
166  fhSubJettiness2CheckRatio_FJ_KT(0x0),
167  fhSubJettiness2CheckRatio_FJ_CA(0x0),
168  fhSubJettiness2CheckRatio_FJ_WTA_KT(0x0),
169  fhSubJettiness2CheckRatio_FJ_WTA_CA(0x0),
170  fhSubJettiness2CheckRatio_FJ_OP_AKT(0x0),
171  fhSubJettiness2CheckRatio_FJ_OP_KT(0x0),
172  fhSubJettiness2CheckRatio_FJ_OP_CA(0x0),
173  fhSubJettiness2CheckRatio_FJ_OP_WTA_KT(0x0),
174  fhSubJettiness2CheckRatio_FJ_OP_WTA_CA(0x0),
175  fhSubJettiness2CheckRatio_FJ_MIN(0x0),
176  fhSubJettiness2to1CheckRatio_FJ_AKT(0x0),
177  fhSubJettiness2to1CheckRatio_FJ_KT(0x0),
178  fhSubJettiness2to1CheckRatio_FJ_CA(0x0),
179  fhSubJettiness2to1CheckRatio_FJ_WTA_KT(0x0),
180  fhSubJettiness2to1CheckRatio_FJ_WTA_CA(0x0),
181  fhSubJettiness2to1CheckRatio_FJ_OP_AKT(0x0),
182  fhSubJettiness2to1CheckRatio_FJ_OP_KT(0x0),
183  fhSubJettiness2to1CheckRatio_FJ_OP_CA(0x0),
184  fhSubJettiness2to1CheckRatio_FJ_OP_WTA_KT(0x0),
185  fhSubJettiness2to1CheckRatio_FJ_OP_WTA_CA(0x0),
186  fhSubJettiness2to1CheckRatio_FJ_MIN(0x0),
187  fhSubJettiness1_FJ_AKT(0x0),
188  fhSubJettiness1_FJ_KT(0x0),
189  fhSubJettiness1_FJ_CA(0x0),
190  fhSubJettiness1_FJ_WTA_KT(0x0),
191  fhSubJettiness1_FJ_WTA_CA(0x0),
192  fhSubJettiness1_FJ_OP_AKT(0x0),
193  fhSubJettiness1_FJ_OP_KT(0x0),
194  fhSubJettiness1_FJ_OP_CA(0x0),
195  fhSubJettiness1_FJ_OP_WTA_KT(0x0),
196  fhSubJettiness1_FJ_OP_WTA_CA(0x0),
197  fhSubJettiness1_FJ_MIN(0x0),
198  fhSubJettiness2_FJ_AKT(0x0),
199  fhSubJettiness2_FJ_KT(0x0),
200  fhSubJettiness2_FJ_CA(0x0),
201  fhSubJettiness2_FJ_WTA_KT(0x0),
202  fhSubJettiness2_FJ_WTA_CA(0x0),
203  fhSubJettiness2_FJ_OP_AKT(0x0),
204  fhSubJettiness2_FJ_OP_KT(0x0),
205  fhSubJettiness2_FJ_OP_CA(0x0),
206  fhSubJettiness2_FJ_OP_WTA_KT(0x0),
207  fhSubJettiness2_FJ_OP_WTA_CA(0x0),
208  fhSubJettiness2_FJ_MIN(0x0),
209  fhSubJettiness2to1_FJ_AKT(0x0),
210  fhSubJettiness2to1_FJ_KT(0x0),
211  fhSubJettiness2to1_FJ_CA(0x0),
212  fhSubJettiness2to1_FJ_WTA_KT(0x0),
213  fhSubJettiness2to1_FJ_WTA_CA(0x0),
214  fhSubJettiness2to1_FJ_OP_AKT(0x0),
215  fhSubJettiness2to1_FJ_OP_KT(0x0),
216  fhSubJettiness2to1_FJ_OP_CA(0x0),
217  fhSubJettiness2to1_FJ_OP_WTA_KT(0x0),
218  fhSubJettiness2to1_FJ_OP_WTA_CA(0x0),
219  fhSubJettiness2to1_FJ_MIN(0x0),
220  fTreeResponseMatrixAxis(0)
221 
222 {
223  for(Int_t i=0;i<nVar;i++){
224  fShapesVar[i]=0;
225  }
226  SetMakeGeneralHistograms(kTRUE);
227  DefineOutput(1, TList::Class());
228  DefineOutput(2, TTree::Class());
229 }
230 
231 //________________________________________________________________________
233  AliAnalysisTaskEmcalJet(name, kTRUE),
234  fContainer(0),
235  fMinFractionShared(0),
236  fJetShapeType(kData),
237  fJetShapeSub(kNoSub),
238  fJetSelection(kInclusive),
239  fPtThreshold(-9999.),
240  fRMatching(0.2),
241  fPtMinTriggerHadron(20.),
242  fPtMaxTriggerHadron(50.),
243  fRecoilAngularWindow(0.6),
244  fSemigoodCorrect(0),
245  fHolePos(0),
246  fHoleWidth(0),
247  fRandom(0x0),
248  fCentSelectOn(kTRUE),
249  fCentMin(0),
250  fCentMax(10),
251  fSubJetAlgorithm(0),
252  fSubJetRadius(0.1),
253  fSubJetMinPt(1),
254  fJetRadius(0.4),
255  fRMatched(0.2),
256  fSharedFractionPtMin(0.5),
257  fDerivSubtrOrder(0),
258  fFullTree(kFALSE),
259  fBeta_SD(0),
260  fZCut(0.1),
261  fReclusteringAlgorithm(0),
262  fSoftDropOn(0),
263  fNsubMeasure(kFALSE),
264  fhPtTriggerHadron(0x0),
265  fhJetPt(0x0),
266  fhJetPt_1(0x0),
267  fhJetPt_2(0x0),
268  fhJetPhi(0x0),
269  fhJetPhi_1(0x0),
270  fhJetPhi_2(0x0),
271  fhJetEta(0x0),
272  fhJetEta_1(0x0),
273  fhJetEta_2(0x0),
274  fhJetMass(0x0),
275  fhJetMass_1(0x0),
276  fhJetMass_2(0x0),
277  fhJetRadius(0x0),
278  fhJetRadius_1(0x0),
279  fhJetRadius_2(0x0),
280  fhJetCounter(0x0),
281  fhJetCounter_1(0x0),
282  fhJetCounter_2(0x0),
283  fhNumberOfJetTracks(0x0),
284  fhNumberOfJetTracks_1(0x0),
285  fhNumberOfJetTracks_2(0x0),
286  fhSubJetPt(0x0),
287  fhSubJetPt_1(0x0),
288  fhSubJetPt_2(0x0),
289  fhSubJetMass(0x0),
290  fhSubJetMass_1(0x0),
291  fhSubJetMass_2(0x0),
292  fhSubJetRadius(0x0),
293  fhSubJetRadius_1(0x0),
294  fhSubJetRadius_2(0x0),
295  fhSubJetCounter(0x0),
296  fhSubJetCounter_1(0x0),
297  fhSubJetCounter_2(0x0),
298  fhNumberOfSubJetTracks(0x0),
299  fhNumberOfSubJetTracks_1(0x0),
300  fhNumberOfSubJetTracks_2(0x0),
301  fh2PtTriggerHadronJet(0x0),
302  fhPhiTriggerHadronJet(0x0),
303  fhPhiTriggerHadronEventPlane(0x0),
304  fhPhiTriggerHadronEventPlaneTPC(0x0),
305  fhTrackPhi(0x0),
306  fhTrackPhi_Cut(0x0),
307  fh2PtRatio(0x0),
308  fhEventCounter(0x0),
309  fhEventCounter_1(0x0),
310  fhEventCounter_2(0x0),
311  fhSubJettiness1CheckRatio_FJ_AKT(0x0),
312  fhSubJettiness1CheckRatio_FJ_KT(0x0),
313  fhSubJettiness1CheckRatio_FJ_CA(0x0),
314  fhSubJettiness1CheckRatio_FJ_WTA_KT(0x0),
315  fhSubJettiness1CheckRatio_FJ_WTA_CA(0x0),
316  fhSubJettiness1CheckRatio_FJ_OP_AKT(0x0),
317  fhSubJettiness1CheckRatio_FJ_OP_KT(0x0),
318  fhSubJettiness1CheckRatio_FJ_OP_CA(0x0),
319  fhSubJettiness1CheckRatio_FJ_OP_WTA_KT(0x0),
320  fhSubJettiness1CheckRatio_FJ_OP_WTA_CA(0x0),
321  fhSubJettiness1CheckRatio_FJ_MIN(0x0),
322  fhSubJettiness2CheckRatio_FJ_AKT(0x0),
323  fhSubJettiness2CheckRatio_FJ_KT(0x0),
324  fhSubJettiness2CheckRatio_FJ_CA(0x0),
325  fhSubJettiness2CheckRatio_FJ_WTA_KT(0x0),
326  fhSubJettiness2CheckRatio_FJ_WTA_CA(0x0),
327  fhSubJettiness2CheckRatio_FJ_OP_AKT(0x0),
328  fhSubJettiness2CheckRatio_FJ_OP_KT(0x0),
329  fhSubJettiness2CheckRatio_FJ_OP_CA(0x0),
330  fhSubJettiness2CheckRatio_FJ_OP_WTA_KT(0x0),
331  fhSubJettiness2CheckRatio_FJ_OP_WTA_CA(0x0),
332  fhSubJettiness2CheckRatio_FJ_MIN(0x0),
333  fhSubJettiness2to1CheckRatio_FJ_AKT(0x0),
334  fhSubJettiness2to1CheckRatio_FJ_KT(0x0),
335  fhSubJettiness2to1CheckRatio_FJ_CA(0x0),
336  fhSubJettiness2to1CheckRatio_FJ_WTA_KT(0x0),
337  fhSubJettiness2to1CheckRatio_FJ_WTA_CA(0x0),
338  fhSubJettiness2to1CheckRatio_FJ_OP_AKT(0x0),
339  fhSubJettiness2to1CheckRatio_FJ_OP_KT(0x0),
340  fhSubJettiness2to1CheckRatio_FJ_OP_CA(0x0),
341  fhSubJettiness2to1CheckRatio_FJ_OP_WTA_KT(0x0),
342  fhSubJettiness2to1CheckRatio_FJ_OP_WTA_CA(0x0),
343  fhSubJettiness2to1CheckRatio_FJ_MIN(0x0),
344  fhSubJettiness1_FJ_AKT(0x0),
345  fhSubJettiness1_FJ_KT(0x0),
346  fhSubJettiness1_FJ_CA(0x0),
347  fhSubJettiness1_FJ_WTA_KT(0x0),
348  fhSubJettiness1_FJ_WTA_CA(0x0),
349  fhSubJettiness1_FJ_OP_AKT(0x0),
350  fhSubJettiness1_FJ_OP_KT(0x0),
351  fhSubJettiness1_FJ_OP_CA(0x0),
352  fhSubJettiness1_FJ_OP_WTA_KT(0x0),
353  fhSubJettiness1_FJ_OP_WTA_CA(0x0),
354  fhSubJettiness1_FJ_MIN(0x0),
355  fhSubJettiness2_FJ_AKT(0x0),
356  fhSubJettiness2_FJ_KT(0x0),
357  fhSubJettiness2_FJ_CA(0x0),
358  fhSubJettiness2_FJ_WTA_KT(0x0),
359  fhSubJettiness2_FJ_WTA_CA(0x0),
360  fhSubJettiness2_FJ_OP_AKT(0x0),
361  fhSubJettiness2_FJ_OP_KT(0x0),
362  fhSubJettiness2_FJ_OP_CA(0x0),
363  fhSubJettiness2_FJ_OP_WTA_KT(0x0),
364  fhSubJettiness2_FJ_OP_WTA_CA(0x0),
365  fhSubJettiness2_FJ_MIN(0x0),
366  fhSubJettiness2to1_FJ_AKT(0x0),
367  fhSubJettiness2to1_FJ_KT(0x0),
368  fhSubJettiness2to1_FJ_CA(0x0),
369  fhSubJettiness2to1_FJ_WTA_KT(0x0),
370  fhSubJettiness2to1_FJ_WTA_CA(0x0),
371  fhSubJettiness2to1_FJ_OP_AKT(0x0),
372  fhSubJettiness2to1_FJ_OP_KT(0x0),
373  fhSubJettiness2to1_FJ_OP_CA(0x0),
374  fhSubJettiness2to1_FJ_OP_WTA_KT(0x0),
375  fhSubJettiness2to1_FJ_OP_WTA_CA(0x0),
376  fhSubJettiness2to1_FJ_MIN(0x0),
377  fTreeResponseMatrixAxis(0)
378 
379 {
380  // Standard constructor.
381  for(Int_t i=0;i<nVar;i++){
382  fShapesVar[i]=0;
383  }
385  DefineOutput(1, TList::Class());
386  DefineOutput(2, TTree::Class());
387 }
388 
389 //________________________________________________________________________
391 {
392  // Destructor.
393 }
394 
395 //________________________________________________________________________
397 {
398  // Create user output.
399 
401 
402  Bool_t oldStatus = TH1::AddDirectoryStatus();
403  TH1::AddDirectory(kFALSE);
404  TH1::AddDirectory(oldStatus);
405  //create a tree used for the MC data and making a 4D response matrix
406  const char* nameoutput = GetOutputSlot(2)->GetContainer()->GetName();
407  fTreeResponseMatrixAxis = new TTree(nameoutput, nameoutput);
408  if (fFullTree){
409  TString *fShapesVarNames = new TString [nVar];
410 
411  fShapesVarNames[0] = "Pt";
412  fShapesVarNames[1] = "Pt_Truth";
413  fShapesVarNames[2] = "Tau1";
414  fShapesVarNames[3] = "Tau1_Truth";
415  fShapesVarNames[4] = "Tau2";
416  fShapesVarNames[5] = "Tau2_Truth";
417  fShapesVarNames[6] = "Tau3";
418  fShapesVarNames[7] = "Tau3_Truth";
419  fShapesVarNames[8] = "OpeningAngle";
420  fShapesVarNames[9] = "OpeningAngle_Truth";
421  fShapesVarNames[10] = "JetMultiplicity";
422  fShapesVarNames[11] = "JetMultiplicity_Truth";
423  fShapesVarNames[12] = "OpeningAngleSD";
424  fShapesVarNames[13] = "OpeningAngleSD_Truth";
425  fShapesVarNames[14] = "Zg";
426  fShapesVarNames[15] = "Zg_Truth";
427  fShapesVarNames[16] = "LeadingTrackPt";
428  fShapesVarNames[17] = "LeadingTrackPt_Truth";
429  fShapesVarNames[18] = "EventPlaneTriggerHadron";
430  fShapesVarNames[19] = "EventPlaneTriggerHadron_Truth";
431  fShapesVarNames[20] = "EventPlaneTPCTriggerHadron";
432  fShapesVarNames[21] = "EventPlaneTPCTriggerHadron_Truth";
433  fShapesVarNames[22] = "DeltaR";
434  fShapesVarNames[23] = "DeltaR_Truth";
435  fShapesVarNames[24] = "Frac1";
436  fShapesVarNames[25] = "Frac1_Truth";
437  fShapesVarNames[26] = "Frac2";
438  fShapesVarNames[27] = "Frac2_Truth";
439  for(Int_t ivar=0; ivar < nVar; ivar++){
440  cout<<"looping over variables"<<endl;
441  fTreeResponseMatrixAxis->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/D", fShapesVarNames[ivar].Data()));
442  }
443  }
444 
445  if (!fFullTree){
446  const Int_t nVarMin = 22;
447  TString *fShapesVarNames = new TString [nVarMin];
448 
449  fShapesVarNames[0] = "Pt";
450  fShapesVarNames[1] = "Pt_Truth";
451  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[1] = "Tau1_Randomisation";
452  fShapesVarNames[2] = "Tau1";
453  fShapesVarNames[3] = "Tau1_Truth";
454  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[3] = "Tau1_ExtraProng_02_15";
455  fShapesVarNames[4] = "Tau2";
456  fShapesVarNames[5] = "Tau2_Truth";
457  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[5] = "Tau1_ExtraProng_02_30";
458  fShapesVarNames[6] = "SubJet1LeadingTrackPt";
459  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[6] = "Tau1_ExtraProng_02_45";
460  fShapesVarNames[7] = "SubJet1LeadingTrackPt_Truth";
461  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[7] = "Tau1_ExtraProng_02_60";
462  fShapesVarNames[8] = "OpeningAngle";
463  fShapesVarNames[9] = "OpeningAngle_Truth";
464  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[9] = "Tau1_ExtraProng_01_30";
465  fShapesVarNames[10] = "SubJet2LeadingTrackPt";
466  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[10] = "Tau1_ExtraProng_01_45";
467  fShapesVarNames[11] = "SubJet2LeadingTrackPt_Truth";
468  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[11] = "Tau1_ExtraProng_03_30";
469  fShapesVarNames[12] = "OpeningAngleSD";
470  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[12] = "Tau1_ExtraProng_03_45";
471  fShapesVarNames[13] = "OpeningAngleSD_Truth";
472  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[13] = "Tau2_Randomisation";
473  fShapesVarNames[14] = "SubJet1Pt";
474  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[14] = "Tau2_ExtraProng_02_15";
475  fShapesVarNames[15] = "SubJet1Pt_Truth";
476  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[15] = "Tau2_ExtraProng_02_30";
477  fShapesVarNames[16] = "LeadingTrackPt";
478  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[16] = "Tau2_ExtraProng_02_45";
479  fShapesVarNames[17] = "LeadingTrackPt_Truth";
480  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[17] = "Tau2_ExtraProng_02_60";
481  fShapesVarNames[18] = "EventPlaneTriggerHadron";
482  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[18] = "Tau2_ExtraProng_01_30";
483  fShapesVarNames[19] = "EventPlaneTriggerHadron_Truth";
484  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[19] = "Tau2_ExtraProng_01_45";
485  fShapesVarNames[20] = "SubJet2Pt";
486  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[20] = "Tau2_ExtraProng_03_30";
487  fShapesVarNames[21] = "SubJet2Pt_Truth";
488  if(fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) fShapesVarNames[21] = "Tau2_ExtraProng_03_45";
489 
490  for(Int_t ivar=0; ivar < nVarMin; ivar++){
491  cout<<"looping over variables"<<endl;
492  fTreeResponseMatrixAxis->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/D", fShapesVarNames[ivar].Data()));
493  }
494  }
495 
497  fhPtTriggerHadron= new TH1F("fhPtTriggerHadron", "fhPtTriggerHadron",1500,-0.5,149.5);
499  fh2PtTriggerHadronJet= new TH2F("fh2PtTriggerHadronJet", "fh2PtTriggerHadronJet",1500,-0.5,149.5,1500,-0.5,149.5);
501  fhPhiTriggerHadronJet= new TH1F("fhPhiTriggerHadronJet", "fhPhiTriggerHadronJet",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
503  fhPhiTriggerHadronEventPlane= new TH1F("fhPhiTriggerHadronEventPlane", "fhPhiTriggerHadronEventPlane",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
505  fhPhiTriggerHadronEventPlaneTPC= new TH1F("fhPhiTriggerHadronEventPlaneTPC", "fhPhiTriggerHadronEventPlaneTPC",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
507  }
509 
510  // fhJetPt= new TH1F("fhJetPt", "Jet Pt", (XBinsJetPtSize)-1, XBinsJetPt);
511  fhJetPt= new TH1F("fhJetPt", "Jet Pt",1500,-0.5,149.5 );
512  fOutput->Add(fhJetPt);
513  //fhJetPhi= new TH1F("fhJetPhi", "Jet Phi", Phi_Bins, Phi_Low, Phi_Up);
514  // fhJetPhi= new TH1F("fhJetPhi", "Jet Phi",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
515  fhJetPhi= new TH1F("fhJetPhi", "Jet Phi",780 , -7, 7);
516  fOutput->Add(fhJetPhi);
517  fhJetEta= new TH1F("fhJetEta", "Jet Eta", Eta_Bins, Eta_Low, Eta_Up);
518  fOutput->Add(fhJetEta);
519  fhJetMass= new TH1F("fhJetMass", "Jet Mass", 4000,-0.5, 39.5);
520  fOutput->Add(fhJetMass);
521  fhJetRadius= new TH1F("fhJetRadius", "Jet Radius", 100, -0.05,0.995);
522  fOutput->Add(fhJetRadius);
523  fhNumberOfJetTracks= new TH1F("fhNumberOfJetTracks", "Number of Tracks within a Jet", 300, -0.5,299.5);
525  fhSubJetRadius= new TH1F("fhSubJetRadius", "SubJet Radius", 100, -0.05,0.995);
526  fOutput->Add(fhSubJetRadius);
527  fhSubJetPt= new TH1F("fhSubJetPt", "SubJet Pt", 1500, -0.5,149.5);
528  fOutput->Add(fhSubJetPt);
529  fhSubJetMass= new TH1F("fhSubJetMass", "Sub Jet Mass", 4000,-0.5, 39.5);
530  fOutput->Add(fhSubJetMass);
531  fhNumberOfSubJetTracks= new TH1F("fhNumberOfSubJetTracks", "Number of Tracks within a Sub Jet", 300, -0.5,299.5);
533  fhJetCounter= new TH1F("fhJetCounter", "Jet Counter", 150, -0.5, 149.5);
534  fOutput->Add(fhJetCounter);
535  fhSubJetCounter = new TH1F("fhSubJetCounter", "SubJet Counter",50, -0.5,49.5);
536  fOutput->Add(fhSubJetCounter);
537  fhEventCounter= new TH1F("fhEventCounter", "Event Counter", 15,0.5,15.5);
538  fOutput->Add(fhEventCounter);
539  fhTrackPhi= new TH1F("fhTrackPhi", "fhTrackPhi",780 , -7, 7);
540  fOutput->Add(fhTrackPhi);
541  fhTrackPhi_Cut= new TH1F("fhTrackPhi_Cut", "fhTrackPhi_Cut",780 , -7, 7);
542  fOutput->Add(fhTrackPhi_Cut);
543  }
545  fhSubJettiness1_FJ_KT= new TH1D("fhSubJettiness1_FJ_KT","fhSubJettiness1_FJ_KT",400,-2,2);
547  fhSubJettiness1_FJ_MIN= new TH1D("fhSubJettiness1_FJ_MIN","fhSubJettiness1_FJ_MIN",400,-2,2);
549  fhSubJettiness2_FJ_KT= new TH1D("fhSubJettiness2_FJ_KT","fhSubJettiness2_FJ_KT",400,-2,2);
551  fhSubJettiness2_FJ_MIN= new TH1D("fhSubJettiness2_FJ_MIN","fhSubJettiness2_FJ_MIN",400,-2,2);
553  fhSubJettiness1CheckRatio_FJ_AKT = new TH2D("fhSubJettiness1CheckRatio_FJ_AKT","fhSubJettiness1CheckRatio_FJ_AKT",400,-2,2,300,-1,2);
555  fhSubJettiness1CheckRatio_FJ_KT= new TH2D("fhSubJettiness1CheckRatio_FJ_KT","fhSubJettiness1CheckRatio_FJ_KT",400,-2,2,300,-1,2);
557  fhSubJettiness1CheckRatio_FJ_CA= new TH2D("fhSubJettiness1CheckRatio_FJ_CA","fhSubJettiness1CheckRatio_FJ_CA",400,-2,2,300,-1,2);
559  fhSubJettiness1CheckRatio_FJ_WTA_KT= new TH2D("fhSubJettiness1CheckRatio_FJ_WTA_KT","fhSubJettiness1CheckRatio_FJ_WTA_KT",400,-2,2,300,-1,2);
561  fhSubJettiness1CheckRatio_FJ_WTA_CA= new TH2D("fhSubJettiness1CheckRatio_FJ_WTA_CA","fhSubJettiness1CheckRatio_FJ_WTA_CA",400,-2,2,300,-1,2);
563  fhSubJettiness1CheckRatio_FJ_OP_AKT= new TH2D("fhSubJettiness1CheckRatio_FJ_OP_AKT","fhSubJettiness1CheckRatio_FJ_OP_AKT",400,-2,2,300,-1,2);
565  fhSubJettiness1CheckRatio_FJ_OP_KT= new TH2D("fhSubJettiness1CheckRatio_FJ_OP_KT","fhSubJettiness1CheckRatio_FJ_OP_KT",400,-2,2,300,-1,2);
567  fhSubJettiness1CheckRatio_FJ_OP_CA= new TH2D("fhSubJettiness1CheckRatio_FJ_OP_CA","fhSubJettiness1CheckRatio_FJ_OP_CA",400,-2,2,300,-1,2);
569  fhSubJettiness1CheckRatio_FJ_OP_WTA_KT= new TH2D("fhSubJettiness1CheckRatio_FJ_OP_WTA_KT","fhSubJettiness1CheckRatio_FJ_OP_WTA_KT",400,-2,2,300,-1,2);
571  fhSubJettiness1CheckRatio_FJ_OP_WTA_CA= new TH2D("fhSubJettiness1CheckRatio_FJ_OP_WTA_CA","fhSubJettiness1CheckRatio_FJ_OP_WTA_CA",400,-2,2,300,-1,2);
573  fhSubJettiness1CheckRatio_FJ_MIN= new TH2D("fhSubJettiness1CheckRatio_FJ_MIN","fhSubJettiness1CheckRatio_FJ_MIN",400,-2,2,300,-1,2);
575  fhSubJettiness2CheckRatio_FJ_AKT = new TH2D("fhSubJettiness2CheckRatio_FJ_AKT","fhSubJettiness2CheckRatio_FJ_AKT",400,-2,2,300,-1,2);
577  fhSubJettiness2CheckRatio_FJ_KT= new TH2D("fhSubJettiness2CheckRatio_FJ_KT","fhSubJettiness2CheckRatio_FJ_KT",400,-2,2,300,-1,2);
579  fhSubJettiness2CheckRatio_FJ_CA= new TH2D("fhSubJettiness2CheckRatio_FJ_CA","fhSubJettiness2CheckRatio_FJ_CA",400,-2,2,300,-1,2);
581  fhSubJettiness2CheckRatio_FJ_WTA_KT= new TH2D("fhSubJettiness2CheckRatio_FJ_WTA_KT","fhSubJettiness2CheckRatio_FJ_WTA_KT",400,-2,2,300,-1,2);
583  fhSubJettiness2CheckRatio_FJ_WTA_CA= new TH2D("fhSubJettiness2CheckRatio_FJ_WTA_CA","fhSubJettiness2CheckRatio_FJ_WTA_CA",400,-2,2,300,-1,2);
585  fhSubJettiness2CheckRatio_FJ_OP_AKT= new TH2D("fhSubJettiness2CheckRatio_FJ_OP_AKT","fhSubJettiness2CheckRatio_FJ_OP_AKT",400,-2,2,300,-1,2);
587  fhSubJettiness2CheckRatio_FJ_OP_KT= new TH2D("fhSubJettiness2CheckRatio_FJ_OP_KT","fhSubJettiness2CheckRatio_FJ_OP_KT",400,-2,2,300,-1,2);
589  fhSubJettiness2CheckRatio_FJ_OP_CA= new TH2D("fhSubJettiness2CheckRatio_FJ_OP_CA","fhSubJettiness2CheckRatio_FJ_OP_CA",400,-2,2,300,-1,2);
591  fhSubJettiness2CheckRatio_FJ_OP_WTA_KT= new TH2D("fhSubJettiness2CheckRatio_FJ_OP_WTA_KT","fhSubJettiness2CheckRatio_FJ_OP_WTA_KT",400,-2,2,300,-1,2);
593  fhSubJettiness2CheckRatio_FJ_OP_WTA_CA= new TH2D("fhSubJettiness2CheckRatio_FJ_OP_WTA_CA","fhSubJettiness2CheckRatio_FJ_OP_WTA_CA",400,-2,2,300,-1,2);
595  fhSubJettiness2CheckRatio_FJ_MIN= new TH2D("fhSubJettiness2CheckRatio_FJ_MIN","fhSubJettiness2CheckRatio_FJ_MIN",400,-2,2,300,-1,2);
597  fhSubJettiness2to1CheckRatio_FJ_AKT = new TH2D("fhSubJettiness2to1CheckRatio_FJ_AKT","fhSubJettiness2to1CheckRatio_FJ_AKT",400,-2,2,300,-1,2);
599  fhSubJettiness2to1CheckRatio_FJ_KT= new TH2D("fhSubJettiness2to1CheckRatio_FJ_KT","fhSubJettiness2to1CheckRatio_FJ_KT",400,-2,2,300,-1,2);
601  fhSubJettiness2to1CheckRatio_FJ_CA= new TH2D("fhSubJettiness2to1CheckRatio_FJ_CA","fhSubJettiness2to1CheckRatio_FJ_CA",400,-2,2,300,-1,2);
603  fhSubJettiness2to1CheckRatio_FJ_WTA_KT= new TH2D("fhSubJettiness2to1CheckRatio_FJ_WTA_KT","fhSubJettiness2to1CheckRatio_FJ_WTA_KT",400,-2,2,300,-1,2);
605  fhSubJettiness2to1CheckRatio_FJ_WTA_CA= new TH2D("fhSubJettiness2to1CheckRatio_FJ_WTA_CA","fhSubJettiness2to1CheckRatio_FJ_WTA_CA",400,-2,2,300,-1,2);
607  fhSubJettiness2to1CheckRatio_FJ_OP_AKT= new TH2D("fhSubJettiness2to1CheckRatio_FJ_OP_AKT","fhSubJettiness2to1CheckRatio_FJ_OP_AKT",400,-2,2,300,-1,2);
609  fhSubJettiness2to1CheckRatio_FJ_OP_KT= new TH2D("fhSubJettiness2to1CheckRatio_FJ_OP_KT","fhSubJettiness2to1CheckRatio_FJ_OP_KT",400,-2,2,300,-1,2);
611  fhSubJettiness2to1CheckRatio_FJ_OP_CA= new TH2D("fhSubJettiness2to1CheckRatio_FJ_OP_CA","fhSubJettiness2to1CheckRatio_FJ_OP_CA",400,-2,2,300,-1,2);
613  fhSubJettiness2to1CheckRatio_FJ_OP_WTA_KT= new TH2D("fhSubJettiness2to1CheckRatio_FJ_OP_WTA_KT","fhSubJettiness2to1CheckRatio_FJ_OP_WTA_KT",400,-2,2,300,-1,2);
615  fhSubJettiness2to1CheckRatio_FJ_OP_WTA_CA= new TH2D("fhSubJettiness2to1CheckRatio_FJ_OP_WTA_CA","fhSubJettiness2to1CheckRatio_FJ_OP_WTA_CA",400,-2,2,300,-1,2);
617  fhSubJettiness2to1CheckRatio_FJ_MIN= new TH2D("fhSubJettiness2to1CheckRatio_FJ_MIN","fhSubJettiness2to1CheckRatio_FJ_MIN",400,-2,2,300,-1,2);
619  fhSubJettiness1_FJ_AKT= new TH1D("fhSubJettiness1_FJ_AKT","fhSubJettiness1_FJ_AKT",400,-2,2);
621  fhSubJettiness1_FJ_CA= new TH1D("fhSubJettiness1_FJ_CA","fhSubJettiness1_FJ_CA",400,-2,2);
623  fhSubJettiness1_FJ_WTA_KT= new TH1D("fhSubJettiness1_FJ_WTA_KT","fhSubJettiness1_FJ_WTA_KT",400,-2,2);
625  fhSubJettiness1_FJ_WTA_CA= new TH1D("fhSubJettiness1_FJ_WTA_CA","fhSubJettiness1_FJ_WTA_CA",400,-2,2);
627  fhSubJettiness1_FJ_OP_AKT= new TH1D("fhSubJettiness1_FJ_OP_AKT","fhSubJettiness1_FJ_OP_AKT",400,-2,2);
629  fhSubJettiness1_FJ_OP_KT= new TH1D("fhSubJettiness1_FJ_OP_KT","fhSubJettiness1_FJ_OP_KT",400,-2,2);
631  fhSubJettiness1_FJ_OP_CA= new TH1D("fhSubJettiness1_FJ_OP_CA","fhSubJettiness1_FJ_OP_CA",400,-2,2);
633  fhSubJettiness1_FJ_OP_WTA_KT= new TH1D("fhSubJettiness1_FJ_OP_WTA_KT","fhSubJettiness1_FJ_OP_WTA_KT",400,-2,2);
635  fhSubJettiness1_FJ_OP_WTA_CA= new TH1D("fhSubJettiness1_FJ_OP_WTA_CA","fhSubJettiness1_FJ_OP_WTA_CA",400,-2,2);
637  fhSubJettiness2_FJ_AKT= new TH1D("fhSubJettiness2_FJ_AKT","fhSubJettiness2_FJ_AKT",400,-2,2);
639  fhSubJettiness2_FJ_CA= new TH1D("fhSubJettiness2_FJ_CA","fhSubJettiness2_FJ_CA",400,-2,2);
641  fhSubJettiness2_FJ_WTA_KT= new TH1D("fhSubJettiness2_FJ_WTA_KT","fhSubJettiness2_FJ_WTA_KT",400,-2,2);
643  fhSubJettiness2_FJ_WTA_CA= new TH1D("fhSubJettiness2_FJ_WTA_CA","fhSubJettiness2_FJ_WTA_CA",400,-2,2);
645  fhSubJettiness2_FJ_OP_AKT= new TH1D("fhSubJettiness2_FJ_OP_AKT","fhSubJettiness2_FJ_OP_AKT",400,-2,2);
647  fhSubJettiness2_FJ_OP_KT= new TH1D("fhSubJettiness2_FJ_OP_KT","fhSubJettiness2_FJ_OP_KT",400,-2,2);
649  fhSubJettiness2_FJ_OP_CA= new TH1D("fhSubJettiness2_FJ_OP_CA","fhSubJettiness2_FJ_OP_CA",400,-2,2);
651  fhSubJettiness2_FJ_OP_WTA_KT= new TH1D("fhSubJettiness2_FJ_OP_WTA_KT","fhSubJettiness2_FJ_OP_WTA_KT",400,-2,2);
653  fhSubJettiness2_FJ_OP_WTA_CA= new TH1D("fhSubJettiness2_FJ_OP_WTA_CA","fhSubJettiness2_FJ_OP_WTA_CA",400,-2,2);
655  fhSubJettiness2to1_FJ_KT= new TH1D("fhSubJettiness2to1_FJ_KT","fhSubJettiness2to1_FJ_KT",400,-2,2);
657  fhSubJettiness2to1_FJ_AKT= new TH1D("fhSubJettiness2to1_FJ_AKT","fhSubJettiness2to1_FJ_AKT",400,-2,2);
659  fhSubJettiness2to1_FJ_CA= new TH1D("fhSubJettiness2to1_FJ_CA","fhSubJettiness2to1_FJ_CA",400,-2,2);
661  fhSubJettiness2to1_FJ_WTA_KT= new TH1D("fhSubJettiness2to1_FJ_WTA_KT","fhSubJettiness2to1_FJ_WTA_KT",400,-2,2);
663  fhSubJettiness2to1_FJ_WTA_CA= new TH1D("fhSubJettiness2to1_FJ_WTA_CA","fhSubJettiness2to1_FJ_WTA_CA",400,-2,2);
665  fhSubJettiness2to1_FJ_OP_AKT= new TH1D("fhSubJettiness2to1_FJ_OP_AKT","fhSubJettiness2to1_FJ_OP_AKT",400,-2,2);
667  fhSubJettiness2to1_FJ_OP_KT= new TH1D("fhSubJettiness2to1_FJ_OP_KT","fhSubJettiness2to1_FJ_OP_KT",400,-2,2);
669  fhSubJettiness2to1_FJ_OP_CA= new TH1D("fhSubJettiness2to1_FJ_OP_CA","fhSubJettiness2to1_FJ_OP_CA",400,-2,2);
671  fhSubJettiness2to1_FJ_OP_WTA_KT= new TH1D("fhSubJettiness2to1_FJ_OP_WTA_KT","fhSubJettiness2to1_FJ_OP_WTA_KT",400,-2,2);
673  fhSubJettiness2to1_FJ_OP_WTA_CA= new TH1D("fhSubJettiness2to1_FJ_OP_WTA_CA","fhSubJettiness2to1_FJ_OP_WTA_CA",400,-2,2);
675  fhSubJettiness2to1_FJ_MIN= new TH1D("fhSubJettiness2to1_FJ_MIN","fhSubJettiness2to1_FJ_MIN",400,-2,2);
677  }
679  fhJetPt_1= new TH1F("fhJetPt_1", "Jet Pt Detector Level",1500,-0.5,149.5 );
680  fOutput->Add(fhJetPt_1);
681  fhJetPt_2= new TH1F("fhJetPt_2", "Jet Pt Particle Level",1500,-0.5,149.5 );
682  fOutput->Add(fhJetPt_2);
683  fhJetPhi_1= new TH1F("fhJetPhi_1", "Jet Phi Detector Level",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
684  fOutput->Add(fhJetPhi_1);
685  fhJetPhi_2= new TH1F("fhJetPhi_2", "Jet Phi Particle Level",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
686  fOutput->Add(fhJetPhi_2);
687  fhJetEta_1= new TH1F("fhJetEta_1", "Jet Eta Detector Level", Eta_Bins, Eta_Low, Eta_Up);
688  fOutput->Add(fhJetEta_1);
689  fhJetEta_2= new TH1F("fhJetEta_2", "Jet Eta Particle Level", Eta_Bins, Eta_Low, Eta_Up);
690  fOutput->Add(fhJetEta_2);
691  fhJetMass_1= new TH1F("fhJetMass_1", "Jet Mass Detector Level", 4000,-0.5, 39.5);
692  fOutput->Add(fhJetMass_1);
693  fhJetMass_2= new TH1F("fhJetMass_2", "Jet Mass Particle Level", 4000,-0.5, 39.5);
694  fOutput->Add(fhJetMass_2);
695  fhJetRadius_1= new TH1F("fhJetRadius_1", "Jet Radius Detector Level", 100, -0.05,0.995);
696  fOutput->Add(fhJetRadius_1);
697  fhJetRadius_2= new TH1F("fhJetRadius_2", "Jet Radius Particle Level", 100, -0.05,0.995);
698  fOutput->Add(fhJetRadius_2);
699  fhNumberOfJetTracks_1= new TH1F("fhNumberOfJetTracks_1", "Number of Tracks within a Jet Detector Level", 300, -0.5,299.5);
701  fhNumberOfJetTracks_2= new TH1F("fhNumberOfJetTracks_2", "Number of Tracks within a Jet Particle Level", 300, -0.5,299.5);
703  fhSubJetRadius_1= new TH1F("fhSubJetRadius_1", "SubJet Radius Detector Level", 100, -0.05,0.995);
705  fhSubJetRadius_2= new TH1F("fhSubJetRadius_2", "SubJet Radius Particle Level", 100, -0.05,0.995);
707  fhSubJetPt_1= new TH1F("fhSubJetPt_1", "SubJet Pt Detector Level", 1500, -0.5,149.5);
708  fOutput->Add(fhSubJetPt_1);
709  fhSubJetPt_2= new TH1F("fhSubJetPt_2", "SubJet Pt Particle Level", 1500, -0.5,149.5);
710  fOutput->Add(fhSubJetPt_2);
711  fhSubJetMass_1= new TH1F("fhSubJetMass_1", "Sub Jet Mass Detector Level", 4000,-0.5, 39.5);
712  fOutput->Add(fhSubJetMass_1);
713  fhSubJetMass_2= new TH1F("fhSubJetMass_2", "Sub Jet Mass Particle Level", 4000,-0.5, 39.5);
714  fOutput->Add(fhSubJetMass_2);
715  fhNumberOfSubJetTracks_1= new TH1F("fhNumberOfSubJetTracks_1", "Number of Tracks within a Sub Jet Detector Level", 300, -0.5,299.5);
717  fhNumberOfSubJetTracks_2= new TH1F("fhNumberOfSubJetTracks_2", "Number of Tracks within a Sub Jet Particle Level", 300, -0.5,299.5);
719  fhJetCounter_1= new TH1F("fhJetCounter_1", "Jet Counter Detector Level", 150, -0.5, 149.5);
720  fOutput->Add(fhJetCounter_1);
721  fhJetCounter_2= new TH1F("fhJetCounter_2", "Jet Counter Particle Level", 150, -0.5, 149.5);
722  fOutput->Add(fhJetCounter_2);
723  fhSubJetCounter_1 = new TH1F("fhSubJetCounter_1", "SubJet Counter Detector Level",50, -0.5,49.5);
725  fhSubJetCounter_2 = new TH1F("fhSubJetCounter_2", "SubJet Counter Particle Level",50, -0.5,49.5);
727  fh2PtRatio= new TH2F("fhPtRatio", "MC pt for detector level vs particle level jets",1500,-0.5,149.5,1500,-0.5,149.5);
728  fOutput->Add(fh2PtRatio);
729  fhEventCounter_1= new TH1F("fhEventCounter_1", "Event Counter Detector Level", 15,0.5,15.5);
731  fhEventCounter_2= new TH1F("fhEventCounter_2", "Event Counter Particle Level", 15,0.5,15.5);
733  fhTrackPhi= new TH1F("fhTrackPhi", "fhTrackPhi",780 , -7, 7);
734  fOutput->Add(fhTrackPhi);
735  fhTrackPhi_Cut= new TH1F("fhTrackPhi_Cut", "fhTrackPhi_Cut",780 , -7, 7);
736  fOutput->Add(fhTrackPhi_Cut);
737  }
739  fhEventCounter= new TH1F("fhEventCounter", "Event Counter", 15,0.5,15.5);
740  fOutput->Add(fhEventCounter);
741  }
742 
743  PostData(1,fOutput);
744  PostData(2,fTreeResponseMatrixAxis);
745  // delete [] fShapesVarNames;
746 }
747 
748 //________________________________________________________________________
750 {
751  // Run analysis code here, if needed. It will be executed before FillHistograms().
752 
753 
754  return kTRUE;
755 }
756 
757 //________________________________________________________________________
759 {
760 
761  //fhEventCounter Key:
762  // 1: Number of events with a Jet Container
763  // 2: Number of Jets without a Jet Container
764  // 3:
765  // 4: Number of Jets found in all events
766  // 5: Number of Jets that were reclustered in all events
767  // 6: Number of SubJets found in all events
768  // 7: Number of events were primary vertext was found
769  // 8: Number of Jets with more than one SubJet
770  // 9:Number of Jets with more than two SubJets
771  // 12:Number of SubJetinessEvents in kData
772 
773 
774  // const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
775  // Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
776  // if(vert) fhEventCounter->Fill(7);
777 
778  //cout << ((AliVAODHeader*)(InputEvent()->GetHeader()))->GetEventplane()<< " "<<fEPV0<<endl;
779  // cout << InputEvent()->GetEventplane()->GetEventplane("Q")<< " "<<fEPV0<<endl;
780  if (fCentSelectOn){
781  if ((fCent>fCentMax) || (fCent<fCentMin)) return 0;
782  }
783 
785  AliTrackContainer *PartCont_Particles = NULL;
786  AliParticleContainer *PartContMC_Particles = NULL;
787 
788  if (fJetShapeSub==kConstSub){
790  else PartCont_Particles = GetTrackContainer(1);
791  }
792  else{
794  else PartCont_Particles = GetTrackContainer(0);
795  }
796 
797  TClonesArray *TracksArray_Particles = NULL;
798  TClonesArray *TracksArrayMC_Particles = NULL;
799  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly && PartContMC_Particles) TracksArrayMC_Particles = PartContMC_Particles->GetArray();
800  else if(PartCont_Particles) TracksArray_Particles = PartCont_Particles->GetArray();
801 
802  AliAODTrack *Track_Particles = 0x0;
803  Int_t NTracks_Particles=0;
804  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly && TracksArrayMC_Particles) NTracks_Particles = TracksArrayMC_Particles->GetEntriesFast();
805  else if(TracksArray_Particles) NTracks_Particles = TracksArray_Particles->GetEntriesFast();
806 
808  if (PartContMC_Particles && TracksArrayMC_Particles){
809  for(Int_t i=0; i < NTracks_Particles; i++){
810  if((Track_Particles = static_cast<AliAODTrack*>(PartContMC_Particles->GetAcceptParticle(i)))){
811  if (!Track_Particles) continue;
812  if(TMath::Abs(Track_Particles->Eta())>0.9) continue;
813  if (Track_Particles->Pt()<0.15) continue;
814  fhTrackPhi->Fill(Track_Particles->Phi());
815  if (Track_Particles->Pt()>=4.0) fhTrackPhi_Cut->Fill(Track_Particles->Phi());
816  }
817  }
818  }
819  }
820  else{
821  if (PartCont_Particles && TracksArray_Particles){
822  for(Int_t i=0; i < NTracks_Particles; i++){
823  if((Track_Particles = static_cast<AliAODTrack*>(PartCont_Particles->GetAcceptTrack(i)))){
824  if (!Track_Particles) continue;
825  if(TMath::Abs(Track_Particles->Eta())>0.9) continue;
826  if (Track_Particles->Pt()<0.15) continue;
827  fhTrackPhi->Fill(Track_Particles->Phi());
828  if (Track_Particles->Pt()>=4.0) fhTrackPhi_Cut->Fill(Track_Particles->Phi());
829  }
830  }
831  }
832  }
834 
835  AliAODTrack *TriggerHadron = 0x0;
836  if (fJetSelection == kRecoil) {
837  //you have to set a flag and the limits of the pT interval for your trigger
839  if (TriggerHadronLabel==-99999) return 0; //Trigger Hadron Not Found
840  AliTrackContainer *PartCont =NULL;
841  AliParticleContainer *PartContMC=NULL;
842  if (fJetShapeSub==kConstSub){
844  else PartCont = GetTrackContainer(1);
845  }
846  else{
848  else PartCont = GetTrackContainer(0);
849  }
850  TClonesArray *TrackArray = NULL;
851  TClonesArray *TrackArrayMC = NULL;
852  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) TrackArrayMC = PartContMC->GetArray();
853  else TrackArray = PartCont->GetArray();
854  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) TriggerHadron = static_cast<AliAODTrack*>(TrackArrayMC->At(TriggerHadronLabel));
855  else TriggerHadron = static_cast<AliAODTrack*>(TrackArray->At(TriggerHadronLabel));
856  if (!TriggerHadron) return 0;//No trigger hadron with label found
857  if(fSemigoodCorrect){
858  Double_t HoleDistance=RelativePhi(TriggerHadron->Phi(),fHolePos);
859  if(TMath::Abs(HoleDistance)+fHoleWidth+fJetRadius>TMath::Pi()-fRecoilAngularWindow) return 0;
860  }
861  fhPtTriggerHadron->Fill(TriggerHadron->Pt()); //Needed for per trigger Normalisation
862  if (fJetShapeType != AliAnalysisTaskSubJetFraction::kGenOnTheFly) fhPhiTriggerHadronEventPlane->Fill(TMath::Abs(RelativePhiEventPlane(fEPV0,TriggerHadron->Phi()))); //fEPV0 is the event plane from AliAnalysisTaskEmcal
863  if (fJetShapeType != AliAnalysisTaskSubJetFraction::kGenOnTheFly) fhPhiTriggerHadronEventPlaneTPC->Fill(TMath::Abs(RelativePhiEventPlane(((AliVAODHeader*)(InputEvent()->GetHeader()))->GetEventplane(),TriggerHadron->Phi()))); //TPC event plane
864  }
865 
866 
869  AliEmcalJet *Jet1 = NULL; //Subtracted hybrid Jet
870  AliEmcalJet *Jet2 = NULL; //Unsubtracted Hybrid Jet
871  AliEmcalJet *Jet3 = NULL; //Detector Level Pythia Jet
872  AliEmcalJet *Jet4 = NULL; //Particle Level Pyhtia Jet
873  AliJetContainer *JetCont1= GetJetContainer(0); //Jet Container for Subtracted Hybrid Jets
874  AliJetContainer *JetCont2= GetJetContainer(1); //Jet Container for Unsubtracted Hybrid Jets
875  AliJetContainer *JetCont3= GetJetContainer(2); //Jet Container for Detector Level Pyhtia Jets
876  AliJetContainer *JetCont4= GetJetContainer(3); //Jet Container for Particle Level Pythia Jets
877  AliEmcalJetFinder *Reclusterer1; //Object containg Subjets from Subtracted Hybrid Jets
878  AliEmcalJetFinder *Reclusterer4; //Object containg Subjets from Particle Level Pythia Jets
879  Bool_t JetsMatched=kFALSE;
880  Bool_t EventCounter=kFALSE;
881  Int_t JetNumber=-1;
882  Double_t JetPtThreshold=-2;
883  fhEventCounter->Fill(1);
884  if(JetCont1) {
885  fhEventCounter->Fill(2);
886  JetCont1->ResetCurrentID();
887  while((Jet1=JetCont1->GetNextAcceptJet())) {
888  if (fJetShapeSub==kConstSub) JetPtThreshold=Jet1->Pt();
889  else JetPtThreshold=Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
890  if ( (!Jet1) || (JetPtThreshold<fPtThreshold)) continue;
891  else {
892  /* if(fSemigoodCorrect){
893  Double_t HoleDistance=RelativePhi(Jet1->Phi(),fHolePos);
894  if(TMath::Abs(HoleDistance)<fHoleWidth) continue;
895  }*/
896  Float_t RecoilDeltaPhi = 0.;
897  if (fJetSelection == kRecoil){
898  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), Jet1->Phi());
899  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
900  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), Jet1->Pt());
901  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), Jet1->Phi()));
902  }
903  if (!EventCounter){
904  fhEventCounter->Fill(3);
905  EventCounter=kTRUE;
906  }
907  if(fJetShapeSub==kConstSub){
908  JetNumber=-1;
909  for(Int_t i = 0; i<JetCont2->GetNJets(); i++) {
910  Jet2 = JetCont2->GetJet(i);
911  if(Jet2->GetLabel()==Jet1->GetLabel()) {
912  JetNumber=i;
913  }
914  }
915  if(JetNumber==-1) continue;
916  Jet2=JetCont2->GetJet(JetNumber);
917  if (JetCont2->AliJetContainer::GetFractionSharedPt(Jet2)<fSharedFractionPtMin) continue;
918  Jet3=Jet2->ClosestJet();
919  }
920  if(!(fJetShapeSub==kConstSub)){
921  if (JetCont1->AliJetContainer::GetFractionSharedPt(Jet1)<fSharedFractionPtMin) continue;
922  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
923  }
924  if (!Jet3) continue;
925  Jet4=Jet3->ClosestJet();
926  if(!Jet4) continue;
927  JetsMatched=kTRUE;
929  if (fJetShapeSub==kConstSub) fShapesVar[0]=Jet1->Pt();
930  else fShapesVar[0]=Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
931  fShapesVar[2]=FjNSubJettiness(Jet1,0,1,0,1,0);
932  fShapesVar[4]=FjNSubJettiness(Jet1,0,2,0,1,0);
933  fShapesVar[6]=FjNSubJettiness(Jet1,0,2,0,1,8);
934  fShapesVar[8]=FjNSubJettiness(Jet1,0,2,0,1,1);
935  // fShapesVar[10]=Jet1->GetNumberOfTracks();
936  fShapesVar[10]=FjNSubJettiness(Jet1,0,2,0,1,9);
937  fShapesVar[12]=FjNSubJettiness(Jet1,0,2,0,1,3,fBeta_SD,fZCut);
938  fShapesVar[14]=FjNSubJettiness(Jet1,0,2,0,1,5,fBeta_SD,fZCut);
939  fShapesVar[16]=Jet1->GetLeadingTrack(JetCont1->GetParticleContainer()->GetArray())->Pt();
941  //fShapesVar[20]=RelativePhiEventPlane(((AliVAODHeader*)(InputEvent()->GetHeader()))->GetEventplane(),Jet1->Phi());
942  fShapesVar[20]=FjNSubJettiness(Jet1,0,2,0,1,6,fBeta_SD,fZCut);
943  if (fFullTree){
944  fShapesVar[22]=FjNSubJettiness(Jet1,0,2,0,1,2);
945  Reclusterer1 = Recluster(Jet1, 0, fSubJetRadius, fSubJetMinPt, fSubJetAlgorithm, "SubJetFinder_1");
946  fShapesVar[24]=SubJetFraction(Jet1, Reclusterer1, 1, 0, kTRUE, kFALSE);
947  fShapesVar[26]=SubJetFraction(Jet1, Reclusterer1, 2, 0, kTRUE, kFALSE);
948  }
949  if (JetsMatched){ //even needed? Not now but might be if you want to fill trees when jets aren't matched too
950  fShapesVar[1]=Jet4->Pt();
951  fShapesVar[3]=FjNSubJettiness(Jet4,3,1,0,1,0);
952  fShapesVar[5]=FjNSubJettiness(Jet4,3,2,0,1,0);
953  fShapesVar[7]=FjNSubJettiness(Jet4,3,2,0,1,8);
954  fShapesVar[9]=FjNSubJettiness(Jet4,3,2,0,1,1);
955  //fShapesVar[11]=Jet4->GetNumberOfTracks();
956  fShapesVar[11]=FjNSubJettiness(Jet4,3,2,0,1,9);
957  fShapesVar[13]=FjNSubJettiness(Jet4,3,2,0,1,3,fBeta_SD,fZCut);
958  fShapesVar[15]=FjNSubJettiness(Jet4,3,2,0,1,5,fBeta_SD,fZCut);
959  fShapesVar[17]=Jet4->GetLeadingTrack(JetCont4->GetParticleContainer()->GetArray())->Pt();
960  fShapesVar[19]=RelativePhiEventPlane(fEPV0,Jet4->Phi());
961  //fShapesVar[21]=RelativePhiEventPlane(((AliVAODHeader*)(InputEvent()->GetHeader()))->GetEventplane(),Jet4->Phi());
962  fShapesVar[21]=FjNSubJettiness(Jet4,3,2,0,1,6,fBeta_SD,fZCut);
963  if (fFullTree){
964  fShapesVar[23]=FjNSubJettiness(Jet4,3,2,0,1,2);
965  Reclusterer4=Recluster(Jet4, 3, fSubJetRadius, 0, fSubJetAlgorithm, "SubJetFinder_4");
966  fShapesVar[25]=SubJetFraction(Jet4, Reclusterer4, 1, 0, kTRUE, kFALSE);
967  fShapesVar[27]=SubJetFraction(Jet4, Reclusterer4, 2, 0, kTRUE, kFALSE);
968  }
969  }
970  else{
971  fShapesVar[1]=-2;
972  fShapesVar[3]=-2;
973  fShapesVar[5]=-2;
974  fShapesVar[7]=-2;
975  fShapesVar[9]=-2;
976  fShapesVar[11]=-2;
977  fShapesVar[13]=-2;
978  fShapesVar[15]=-2;
979  fShapesVar[17]=-2;
980  fShapesVar[19]=-2;
981  fShapesVar[21]=-2;
982  if (fFullTree){
983  fShapesVar[23]=-2;
984  fShapesVar[25]=-2;
985  fShapesVar[27]=-2;
986  }
987  }
988  fTreeResponseMatrixAxis->Fill();
989  JetsMatched=kFALSE;
990  }
991  }
992  }
993  }
994 
997  AliEmcalJet *Jet1 = NULL; //Detector Level Jet
998  AliEmcalJet *Jet2 = NULL; //Particle Level Jet
999  AliJetContainer *JetCont1= GetJetContainer(0); //Jet Container for Detector Level Pythia
1000  AliJetContainer *JetCont2= GetJetContainer(1); //Jet Container for Particle Level Pythia
1001  AliEmcalJetFinder *Reclusterer1; //Object containg Subjets Detector Level
1002  AliEmcalJetFinder *Reclusterer2; //Object containg Subjets Particle Level
1003  Int_t JetCounter1=0; //Counts number of jets in event
1004  Int_t JetCounter2=0; //Counts number of jets in event
1005  Double_t JetPhi1=0;
1006  Double_t JetPhi2=0;
1007  Bool_t JetsMatched=kFALSE;
1008  Double_t Pythia_Event_Weight=1;
1009  Bool_t EventCounter=kFALSE;
1010  fhEventCounter_1->Fill(1);
1011  if(JetCont1) {
1012  fhEventCounter_1->Fill(2); //Number of events with a jet container
1013  JetCont1->ResetCurrentID();
1014  while((Jet1=JetCont1->GetNextAcceptJet())) {
1015  if( (!Jet1) || ((Jet1->Pt())<fPtThreshold)) {
1016  // fhEventCounter_1->Fill(3); //events where the jet had a problem
1017  continue;
1018  }
1019  else {
1020  /* if(fSemigoodCorrect){
1021  Double_t HoleDistance=RelativePhi(Jet1->Phi(),fHolePos);
1022  if(TMath::Abs(HoleDistance)<fHoleWidth) continue;
1023  }*/
1024  Float_t RecoilDeltaPhi = 0.;
1025  if (fJetSelection == kRecoil){
1026  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), Jet1->Phi());
1027  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
1028  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), Jet1->Pt());
1029  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), Jet1->Phi()));
1030  }
1031  if (!EventCounter){
1032  fhEventCounter_1->Fill(3);
1033  EventCounter=kTRUE;
1034  }
1035  if((Jet1->GetNumberOfTracks())==0){
1036  fhEventCounter_1->Fill(10); //zero track jets
1037  }
1038  if((Jet1->GetNumberOfTracks())==1){
1039  fhEventCounter_1->Fill(11); //one track jets
1040  }
1041  fhEventCounter_1->Fill(4); //Number of Jets found in all events
1042  JetCounter1++;
1043  fhJetPt_1->Fill(Jet1->Pt());
1044  JetPhi1=Jet1->Phi();
1045  if(JetPhi1 < -1*TMath::Pi()) JetPhi1 += (2*TMath::Pi());
1046  else if (JetPhi1 > TMath::Pi()) JetPhi1 -= (2*TMath::Pi());
1047  fhJetPhi_1->Fill(JetPhi1);
1048  fhJetEta_1->Fill(Jet1->Eta());
1049  fhJetMass_1->Fill(Jet1->M());
1050  fhJetRadius_1->Fill(TMath::Sqrt((Jet1->Area()/TMath::Pi()))); //Radius of Jets per event
1052  if((Jet2 = Jet1->ClosestJet())){
1053  JetsMatched=kTRUE;
1054  if((Jet2->GetNumberOfTracks())==0){
1055  fhEventCounter_2->Fill(10); //zero track jets
1056  }
1057  if((Jet2->GetNumberOfTracks())==1){
1058  fhEventCounter_2->Fill(11); //one track jets
1059  }
1060  fhEventCounter_2->Fill(4); //Number of Jets found in all events
1061  JetCounter2++;
1062  fhJetPt_2->Fill(Jet2->Pt());
1063  JetPhi2=Jet2->Phi();
1064  if(JetPhi2 < -1*TMath::Pi()) JetPhi2 += (2*TMath::Pi());
1065  else if (JetPhi2 > TMath::Pi()) JetPhi2 -= (2*TMath::Pi());
1066  fhJetPhi_2->Fill(JetPhi2);
1067  fhJetEta_2->Fill(Jet2->Eta());
1068  fhJetMass_2->Fill(Jet2->M());
1069  fhJetRadius_2->Fill(TMath::Sqrt((Jet2->Area()/TMath::Pi()))); //Radius of Jets per event
1071  fh2PtRatio->Fill(Jet1->Pt(),Jet2->Pt(),Pythia_Event_Weight);
1072  }
1073  else {
1074  fhEventCounter_2->Fill(1);
1075  //continue;
1076  }
1078 
1079 
1080  fShapesVar[0]=Jet1->Pt();
1081  fShapesVar[2]=FjNSubJettiness(Jet1,0,1,0,1,0);
1082  fShapesVar[4]=FjNSubJettiness(Jet1,0,2,0,1,0);
1083  fShapesVar[6]=FjNSubJettiness(Jet1,0,2,0,1,8);
1084  fShapesVar[8]=FjNSubJettiness(Jet1,0,2,0,1,1);
1085  // fShapesVar[10]=Jet1->GetNumberOfTracks();
1086  fShapesVar[10]=FjNSubJettiness(Jet1,0,2,0,1,9);
1087  fShapesVar[12]=FjNSubJettiness(Jet1,0,2,0,1,3,fBeta_SD,fZCut);
1088  fShapesVar[14]=FjNSubJettiness(Jet1,0,2,0,1,5,fBeta_SD,fZCut);
1089  fShapesVar[16]=Jet1->GetLeadingTrack(JetCont1->GetParticleContainer()->GetArray())->Pt();
1090  fShapesVar[18]=-2; //event plane calculation only needed for PbPb recoils
1091  fShapesVar[20]=FjNSubJettiness(Jet1,0,2,0,1,6,fBeta_SD,fZCut);
1092  Reclusterer1 = Recluster(Jet1, 0, fSubJetRadius, fSubJetMinPt, fSubJetAlgorithm, "SubJetFinder_1");
1093  if (fFullTree){
1094  fShapesVar[22]=FjNSubJettiness(Jet1,0,2,0,1,2);
1095  fShapesVar[24]=SubJetFraction(Jet1, Reclusterer1, 1, 0, kTRUE, kFALSE);
1096  fShapesVar[26]=SubJetFraction(Jet1, Reclusterer1, 2, 0, kTRUE, kFALSE);
1097  }
1098  if (JetsMatched){ //even needed? Not now but might be if you want to fill trees when jets aren't matched too
1099  fShapesVar[1]=Jet2->Pt();
1100  fShapesVar[3]=FjNSubJettiness(Jet2,1,1,0,1,0);
1101  fShapesVar[5]=FjNSubJettiness(Jet2,1,2,0,1,0);
1102  fShapesVar[7]=FjNSubJettiness(Jet2,1,2,0,1,8);
1103  fShapesVar[9]=FjNSubJettiness(Jet2,1,2,0,1,1);
1104  // fShapesVar[11]=Jet2->GetNumberOfTracks();
1105  fShapesVar[11]=FjNSubJettiness(Jet2,1,2,0,1,9);
1106  fShapesVar[13]=FjNSubJettiness(Jet2,1,2,0,1,3,fBeta_SD,fZCut);
1107  fShapesVar[15]=FjNSubJettiness(Jet2,1,2,0,1,5,fBeta_SD,fZCut);
1108  fShapesVar[17]=Jet2->GetLeadingTrack(JetCont2->GetParticleContainer()->GetArray())->Pt();
1109  fShapesVar[19]=-2;
1110  fShapesVar[21]=FjNSubJettiness(Jet2,1,2,0,1,6,fBeta_SD,fZCut);
1111  Reclusterer2 = Recluster(Jet2, 1, fSubJetRadius, 0, fSubJetAlgorithm, "SubJetFinder_2");
1112  if (fFullTree){
1113  fShapesVar[23]=FjNSubJettiness(Jet2,1,2,0,1,2);
1114  fShapesVar[25]=SubJetFraction(Jet2, Reclusterer2, 1, 0, kTRUE, kFALSE);
1115  fShapesVar[27]=SubJetFraction(Jet2, Reclusterer2, 2, 0, kTRUE, kFALSE);
1116  }
1117  }
1118  else{
1119  fShapesVar[1]=-2;
1120  fShapesVar[3]=-2;
1121  fShapesVar[5]=-2;
1122  fShapesVar[7]=-2;
1123  fShapesVar[9]=-2;
1124  fShapesVar[11]=-2;
1125  fShapesVar[13]=-2;
1126  fShapesVar[15]=-2;
1127  fShapesVar[17]=-2;
1128  fShapesVar[19]=-2;
1129  fShapesVar[21]=-2;
1130  if (fFullTree){
1131  fShapesVar[23]=-2;
1132  fShapesVar[25]=-2;
1133  fShapesVar[27]=-2;
1134  }
1135  }
1136  fTreeResponseMatrixAxis->Fill();
1137 
1138  fhSubJetCounter_1->Fill(Reclusterer1->GetNumberOfJets());
1139  for (Int_t i= 0; i<Reclusterer1->GetNumberOfJets(); i++){
1140  fhEventCounter_1->Fill(6); //Number of overall subjets in all events
1141  fhSubJetPt_1->Fill(Reclusterer1->GetJet(i)->Pt());
1142  fhSubJetMass_1->Fill(Reclusterer1->GetJet(i)->M());
1143  fhNumberOfSubJetTracks_1->Fill(Reclusterer1->GetJet(i)->GetNumberOfTracks());
1144  fhSubJetRadius_1->Fill(TMath::Sqrt((Reclusterer1->GetJet(i)->Area()/TMath::Pi()))); //Radius of SubJets per event
1145  }
1146  if(JetsMatched){
1147  fhSubJetCounter_2->Fill(Reclusterer2->GetNumberOfJets());
1148  for (Int_t i= 0; i<Reclusterer2->GetNumberOfJets(); i++){
1149  fhEventCounter_2->Fill(6); //Number of overall subjets in all events
1150  fhSubJetPt_2->Fill(Reclusterer2->GetJet(i)->Pt());
1151  fhSubJetMass_2->Fill(Reclusterer2->GetJet(i)->M());
1152  fhNumberOfSubJetTracks_2->Fill(Reclusterer2->GetJet(i)->GetNumberOfTracks());
1153  fhSubJetRadius_2->Fill(TMath::Sqrt((Reclusterer2->GetJet(i)->Area()/TMath::Pi()))); //Radius of SubJets per event
1154  }
1155  }
1156  JetsMatched=kFALSE;
1157  }
1158  }
1159  fhJetCounter_1->Fill(JetCounter1); //Number of Jets in Each Event Particle Level
1160  fhJetCounter_2->Fill(JetCounter2); //Number of Jets in Each Event Detector Level
1161  }
1162  //else {fhEventCounter_1->Fill(3);} //Events with no jet container
1163  }
1164 
1165 
1166 
1168  AliEmcalJet *Jet1 = NULL; //Original Jet in the event
1169  AliJetContainer *JetCont= GetJetContainer(0); //Jet Container for event
1170  Int_t JetCounter=0; //Counts number of jets in event
1171  Double_t JetPhi=0;
1172  Bool_t EventCounter=kFALSE;
1173  Double_t JetPt_ForThreshold=0;
1174  fhEventCounter->Fill(1);
1175  if(JetCont) {
1176  fhEventCounter->Fill(2); //Number of events with a jet container
1177  JetCont->ResetCurrentID();
1178  while((Jet1=JetCont->GetNextAcceptJet())) {
1179  if(!Jet1) continue;
1180  if (fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) JetPt_ForThreshold = Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
1181  else JetPt_ForThreshold = Jet1->Pt();
1182  if(JetPt_ForThreshold<fPtThreshold) {
1183  //fhEventCounter->Fill(3); //events where the jet had a problem
1184  continue;
1185  }
1186  else {
1187  /* if(fSemigoodCorrect){
1188  Double_t HoleDistance=RelativePhi(Jet1->Phi(),fHolePos);
1189  if(TMath::Abs(HoleDistance)<fHoleWidth) continue;
1190  }*/
1191  Float_t RecoilDeltaPhi = 0.;
1192  if (fJetSelection == kRecoil){
1193  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), Jet1->Phi());
1194  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
1195  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), Jet1->Pt());
1196  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), Jet1->Phi()));
1197  }
1198  if (!EventCounter){
1199  fhEventCounter->Fill(3);
1200  EventCounter=kTRUE;
1201  }
1202  if((Jet1->GetNumberOfTracks())==0){
1203  fhEventCounter->Fill(10); //zero track jets
1204  }
1205  if((Jet1->GetNumberOfTracks())==1){
1206  fhEventCounter->Fill(11); //one track jets
1207  }
1208  fhEventCounter->Fill(4); //Number of Jets found in all events
1209  JetCounter++;
1210  fhJetPt->Fill(Jet1->Pt());
1211  JetPhi=Jet1->Phi();
1212  //if(JetPhi < -1*TMath::Pi()) JetPhi += (2*TMath::Pi());
1213  //else if (JetPhi > TMath::Pi()) JetPhi -= (2*TMath::Pi());
1214  fhJetPhi->Fill(JetPhi);
1215  fhJetEta->Fill(Jet1->Eta());
1216  fhJetMass->Fill(Jet1->M());
1217  fhJetRadius->Fill(TMath::Sqrt((Jet1->Area()/TMath::Pi()))); //Radius of Jets per event
1218  fhNumberOfJetTracks->Fill(Jet1->GetNumberOfTracks());
1219  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) fShapesVar[0]= Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
1220  else fShapesVar[0]=Jet1->Pt();
1221  fShapesVar[2]=FjNSubJettiness(Jet1,0,1,0,1,0);
1222  fShapesVar[4]=FjNSubJettiness(Jet1,0,2,0,1,0);
1223  fShapesVar[6]=FjNSubJettiness(Jet1,0,2,0,1,8);
1224  fShapesVar[8]=FjNSubJettiness(Jet1,0,2,0,1,1);
1225  //fShapesVar[10]=Jet1->GetNumberOfTracks();
1226  fShapesVar[10]=FjNSubJettiness(Jet1,0,2,0,1,9);
1227  fShapesVar[12]=FjNSubJettiness(Jet1,0,2,0,1,3,fBeta_SD,fZCut);
1228  fShapesVar[14]=FjNSubJettiness(Jet1,0,2,0,1,5,fBeta_SD,fZCut);
1229  fShapesVar[16]=Jet1->GetLeadingTrack(JetCont->GetParticleContainer()->GetArray())->Pt();
1230  fShapesVar[18]=-2; //event plane calculation not needed for data
1231  fShapesVar[20]=FjNSubJettiness(Jet1,0,2,0,1,6,fBeta_SD,fZCut);
1232  AliEmcalJetFinder *Reclusterer1 = Recluster(Jet1, 0, fSubJetRadius, fSubJetMinPt, fSubJetAlgorithm, "SubJetFinder");
1233  if (fFullTree){
1234  fShapesVar[22]=FjNSubJettiness(Jet1,0,2,0,1,2);
1235  fShapesVar[24]=SubJetFraction(Jet1, Reclusterer1, 1, 0, kTRUE, kFALSE);
1236  fShapesVar[26]=SubJetFraction(Jet1, Reclusterer1, 2, 0, kTRUE, kFALSE);
1237  }
1238  fShapesVar[1]=-2;
1239  fShapesVar[3]=-2;
1240  fShapesVar[5]=-2;
1241  fShapesVar[7]=-2;
1242  fShapesVar[9]=-2;
1243  fShapesVar[11]=-2;
1244  fShapesVar[13]=-2;
1245  fShapesVar[15]=-2;
1246  fShapesVar[17]=-2;
1247  fShapesVar[19]=-2;
1248  fShapesVar[21]=-2;
1249  if (fFullTree){
1250  fShapesVar[23]=-2;
1251  fShapesVar[25]=-2;
1252  fShapesVar[27]=-2;
1253  }
1254  fTreeResponseMatrixAxis->Fill();
1255  fhSubJetCounter->Fill(Reclusterer1->GetNumberOfJets());
1256  for (Int_t i= 0; i<Reclusterer1->GetNumberOfJets(); i++){
1257  fhEventCounter->Fill(6); //Number of overall subjets in all events
1258  fhSubJetPt->Fill(Reclusterer1->GetJet(i)->Pt());
1259  fhSubJetMass->Fill(Reclusterer1->GetJet(i)->M());
1260  fhNumberOfSubJetTracks->Fill(Reclusterer1->GetJet(i)->GetNumberOfTracks());
1261  fhSubJetRadius->Fill(TMath::Sqrt((Reclusterer1->GetJet(i)->Area()/TMath::Pi()))); //Radius of SubJets per event
1262  }
1263  }
1264  }
1265  fhJetCounter->Fill(JetCounter); //Number of Jets in Each Event
1266  }
1267  //else {fhEventCounter->Fill(2);} //Events with no jet container
1268  }
1269 
1270 
1272 
1273  AliEmcalJet *Jet1 = NULL; //Original Jet in the event
1274  AliJetContainer *JetCont= GetJetContainer(0); //Jet Container for event
1275  Int_t JetCounter=0; //Counts number of jets in event
1276  Double_t JetPhi=0;
1277  Bool_t EventCounter=kFALSE;
1278  Double_t JetPt_ForThreshold=0;
1279  fhEventCounter->Fill(1);
1280  if(JetCont) {
1281  fhEventCounter->Fill(2); //Number of events with a jet container
1282  JetCont->ResetCurrentID();
1283  while((Jet1=JetCont->GetNextAcceptJet())) {
1284  // if (((Jet1=JetCont->GetLeadingJet("rho")) && (fPtThreshold<=Jet1->Pt()-(GetRhoVal(0)*Jet1->Area())))){
1285  if(!Jet1) continue;
1286  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) JetPt_ForThreshold = Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
1287  else JetPt_ForThreshold = Jet1->Pt();
1288  if(JetPt_ForThreshold<fPtThreshold) {
1289  //fhEventCounter->Fill(3); //events where the jet had a problem
1290  continue;
1291  }
1292  else {
1293  /* if(fSemigoodCorrect){
1294  Double_t HoleDistance=RelativePhi(Jet1->Phi(),fHolePos);
1295  if(TMath::Abs(HoleDistance)<fHoleWidth) continue;
1296  }*/
1297  Float_t RecoilDeltaPhi = 0.;
1298  if (fJetSelection == kRecoil){
1299  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), Jet1->Phi());
1300  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
1301  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), Jet1->Pt());
1302  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), Jet1->Phi()));
1303  }
1304  if (!EventCounter){
1305  fhEventCounter->Fill(3);
1306  EventCounter=kTRUE;
1307  }
1308  if((Jet1->GetNumberOfTracks())==0){
1309  fhEventCounter->Fill(10); //zero track jets
1310  }
1311  if((Jet1->GetNumberOfTracks())==1){
1312  fhEventCounter->Fill(11); //one track jets
1313  }
1314  fhEventCounter->Fill(4); //Number of Jets found in all events
1315  JetCounter++;
1316  fhJetPt->Fill(Jet1->Pt());
1317  JetPhi=Jet1->Phi();
1318  if(JetPhi < -1*TMath::Pi()) JetPhi += (2*TMath::Pi());
1319  else if (JetPhi > TMath::Pi()) JetPhi -= (2*TMath::Pi());
1320  fhJetPhi->Fill(JetPhi);
1321  fhJetEta->Fill(Jet1->Eta());
1322  fhJetMass->Fill(Jet1->M());
1323  fhJetRadius->Fill(TMath::Sqrt((Jet1->Area()/TMath::Pi()))); //Radius of Jets per event
1324  fhNumberOfJetTracks->Fill(Jet1->GetNumberOfTracks());
1325  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> Randomised_Jet=ModifyJet(Jet1,0,"Randomise");
1326  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_02_15=ModifyJet(Jet1,0,"AddExtraProng_02_15");
1327  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_02_30=ModifyJet(Jet1,0,"AddExtraProng_02_30");
1328  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_02_45=ModifyJet(Jet1,0,"AddExtraProng_02_45");
1329  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_02_60=ModifyJet(Jet1,0,"AddExtraProng_02_60");
1330  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_01_30=ModifyJet(Jet1,0,"AddExtraProng_01_30");
1331  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_01_45=ModifyJet(Jet1,0,"AddExtraProng_01_45");
1332  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_03_30=ModifyJet(Jet1,0,"AddExtraProng_03_30");
1333  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> ExtraProng_Jet_03_45=ModifyJet(Jet1,0,"AddExtraProng_03_45");
1334  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) fShapesVar[0]= Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
1335  else fShapesVar[0]=Jet1->Pt();
1336  fShapesVar[2]=FjNSubJettiness(Jet1,0,1,0,1,0);
1337  fShapesVar[4]=FjNSubJettiness(Jet1,0,2,0,1,0);
1338  fShapesVar[6]=FjNSubJettinessFastJet(ExtraProng_Jet_02_45,0,1,0,1,0);
1339  fShapesVar[8]=FjNSubJettiness(Jet1,0,2,0,1,1);
1340  // fShapesVar[10]=Jet1->GetNumberOfTracks();
1341  fShapesVar[10]=FjNSubJettinessFastJet(ExtraProng_Jet_01_45,0,1,0,1,0);
1342  fShapesVar[12]=FjNSubJettinessFastJet(ExtraProng_Jet_03_45,0,1,0,1,0);
1343  fShapesVar[14]=FjNSubJettinessFastJet(ExtraProng_Jet_02_15,0,2,0,1,0);
1344  fShapesVar[16]=FjNSubJettinessFastJet(ExtraProng_Jet_02_45,0,2,0,1,0);
1345  fShapesVar[18]=FjNSubJettinessFastJet(ExtraProng_Jet_01_30,0,2,0,1,0);
1346  fShapesVar[20]=FjNSubJettinessFastJet(ExtraProng_Jet_03_30,0,2,0,1,0);
1347  AliEmcalJetFinder *Reclusterer1 = Recluster(Jet1, 0, fSubJetRadius, fSubJetMinPt, fSubJetAlgorithm, "SubJetFinder");
1348  if (fFullTree){
1349  fShapesVar[22]=FjNSubJettiness(Jet1,0,2,0,1,2);
1350  fShapesVar[24]=SubJetFraction(Jet1, Reclusterer1, 1, 0, kTRUE, kFALSE);
1351  fShapesVar[26]=SubJetFraction(Jet1, Reclusterer1, 2, 0, kTRUE, kFALSE);
1352  }
1353  fShapesVar[1]=FjNSubJettinessFastJet(Randomised_Jet,0,1,0,1,0);
1354  // fShapesVar[3]=FjNSubJettiness(std::unique_ptr<AliEmcalJet>(ModifyJet(Jet1,0,"Randomise")).get(),0,1,0,1,0);
1355  //fShapesVar[5]=FjNSubJettiness(std::unique_ptr<AliEmcalJet>(ModifyJet(Jet1,0,"AddExtraProng_02_30")).get(),0,1,0,1,0);
1356  fShapesVar[3]=FjNSubJettinessFastJet(ExtraProng_Jet_02_15,0,1,0,1,0);
1357  fShapesVar[5]=FjNSubJettinessFastJet(ExtraProng_Jet_02_30,0,1,0,1,0);
1358  fShapesVar[7]=FjNSubJettinessFastJet(ExtraProng_Jet_02_60,0,1,0,1,0);
1359  fShapesVar[9]=FjNSubJettinessFastJet(ExtraProng_Jet_01_30,0,1,0,1,0);
1360  fShapesVar[11]=FjNSubJettinessFastJet(ExtraProng_Jet_03_30,0,1,0,1,0);
1361  fShapesVar[13]=FjNSubJettinessFastJet(Randomised_Jet,0,2,0,1,0);
1362  fShapesVar[15]=FjNSubJettinessFastJet(ExtraProng_Jet_02_30,0,2,0,1,0);
1363  fShapesVar[17]=FjNSubJettinessFastJet(ExtraProng_Jet_02_60,0,2,0,1,0);
1364  fShapesVar[19]=FjNSubJettinessFastJet(ExtraProng_Jet_01_45,0,2,0,1,0);
1365  fShapesVar[21]=FjNSubJettinessFastJet(ExtraProng_Jet_03_45,0,2,0,1,0);
1366  if (fFullTree){
1367  fShapesVar[23]=-2;
1368  fShapesVar[25]=-2;
1369  fShapesVar[27]=-2;
1370  }
1371  fTreeResponseMatrixAxis->Fill();
1372 
1373  fhSubJetCounter->Fill(Reclusterer1->GetNumberOfJets());
1374  for (Int_t i= 0; i<Reclusterer1->GetNumberOfJets(); i++){
1375  fhEventCounter->Fill(6); //Number of overall subjets in all events
1376  fhSubJetPt->Fill(Reclusterer1->GetJet(i)->Pt());
1377  fhSubJetMass->Fill(Reclusterer1->GetJet(i)->M());
1378  fhNumberOfSubJetTracks->Fill(Reclusterer1->GetJet(i)->GetNumberOfTracks());
1379  fhSubJetRadius->Fill(TMath::Sqrt((Reclusterer1->GetJet(i)->Area()/TMath::Pi()))); //Radius of SubJets per event
1380  }
1381  delete Randomised_Jet.second;
1382  delete ExtraProng_Jet_02_15.second;
1383  delete ExtraProng_Jet_02_30.second;
1384  delete ExtraProng_Jet_02_45.second;
1385  delete ExtraProng_Jet_02_60.second;
1386  delete ExtraProng_Jet_01_30.second;
1387  delete ExtraProng_Jet_01_45.second;
1388  delete ExtraProng_Jet_03_30.second;
1389  delete ExtraProng_Jet_03_45.second;
1390  }
1391  }
1392  fhJetCounter->Fill(JetCounter); //Number of Jets in Each Event
1393  }
1394  //else {fhEventCounter->Fill(2);} //Events with no jet container
1395  }
1396  return kTRUE;
1397 }
1398 //________________________________________________________________________
1400 
1401  if(Phi < -1*TMath::Pi()) Phi += (2*TMath::Pi());
1402  else if (Phi > TMath::Pi()) Phi -= (2*TMath::Pi());
1403  Double_t DeltaPhi=Phi-EventPlane;
1404  if(DeltaPhi < -1*TMath::Pi()) DeltaPhi += (2*TMath::Pi());
1405  else if (DeltaPhi > TMath::Pi()) DeltaPhi -= (2*TMath::Pi());
1406  return DeltaPhi;
1407 }
1408 //________________________________________________________________________
1410 
1411  if(Phi1 < -1*TMath::Pi()) Phi1 += (2*TMath::Pi()); // Turns the range of 0to2Pi into -PitoPi ???????????
1412  else if (Phi1 > TMath::Pi()) Phi1 -= (2*TMath::Pi());
1413  if(Phi2 < -1*TMath::Pi()) Phi2 += (2*TMath::Pi());
1414  else if (Phi2 > TMath::Pi()) Phi2 -= (2*TMath::Pi());
1415  Double_t DeltaPhi=Phi2-Phi1;
1416  if(DeltaPhi < -1*TMath::Pi()) DeltaPhi += (2*TMath::Pi());
1417  else if (DeltaPhi > TMath::Pi()) DeltaPhi -= (2*TMath::Pi());
1418  return DeltaPhi;
1419 }
1420 
1421 
1422 //--------------------------------------------------------------------------
1424 
1425  AliTrackContainer *PartCont = NULL;
1426  AliParticleContainer *PartContMC = NULL;
1427 
1428 
1429  if (fJetShapeSub==kConstSub){
1431  else PartCont = GetTrackContainer(1);
1432  }
1433  else{
1435  else PartCont = GetTrackContainer(0);
1436  }
1437 
1438  TClonesArray *TracksArray = NULL;
1439  TClonesArray *TracksArrayMC = NULL;
1440  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) TracksArrayMC = PartContMC->GetArray();
1441  else TracksArray = PartCont->GetArray();
1442 
1444  if(!PartContMC || !TracksArrayMC) return -99999;
1445  }
1446  else {
1447  if(!PartCont || !TracksArray) return -99999;
1448  }
1449 
1450  AliAODTrack *Track = 0x0;
1451  Int_t Trigger_Index[100];
1452  for (Int_t i=0; i<100; i++) Trigger_Index[i] = 0;
1453  Int_t Trigger_Counter = 0;
1454  Int_t NTracks=0;
1455  if (fJetShapeType == AliAnalysisTaskSubJetFraction::kGenOnTheFly) NTracks = TracksArrayMC->GetEntriesFast();
1456  else NTracks = TracksArray->GetEntriesFast();
1457  for(Int_t i=0; i < NTracks; i++){
1459  if((Track = static_cast<AliAODTrack*>(PartContMC->GetAcceptParticle(i)))){
1460  if (!Track) continue;
1461  if(TMath::Abs(Track->Eta())>0.9) continue;
1462  if (Track->Pt()<0.15) continue;
1463  if ((Track->Pt() >= PtMin) && (Track->Pt()< PtMax)) {
1464  Trigger_Index[Trigger_Counter] = i;
1465  Trigger_Counter++;
1466  }
1467  }
1468  }
1469  else{
1470  if((Track = static_cast<AliAODTrack*>(PartCont->GetAcceptTrack(i)))){
1471  if (!Track) continue;
1472  if(TMath::Abs(Track->Eta())>0.9) continue;
1473  if (Track->Pt()<0.15) continue;
1474  if ((Track->Pt() >= PtMin) && (Track->Pt()< PtMax)) {
1475  Trigger_Index[Trigger_Counter] = i;
1476  Trigger_Counter++;
1477  }
1478  }
1479  }
1480  }
1481  if (Trigger_Counter == 0) return -99999;
1482  Int_t RandomNumber = 0, Index = 0 ;
1483  TRandom3* Random = new TRandom3(0);
1484  RandomNumber = Random->Integer(Trigger_Counter);
1485  Index = Trigger_Index[RandomNumber];
1486  return Index;
1487 }
1488 
1489 
1490 //--------------------------------------------------------------------------
1492 
1493  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1494  Double_t Angularity_Numerator=0; //Reset these values
1495  Double_t Angularity_Denominator=0;
1496  AliVParticle *Particle=0x0;
1497  Double_t DeltaPhi=0;
1498 
1499 
1500  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks (particles in the jet
1501  Particle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
1502  if(!Particle) continue;
1503  DeltaPhi=RelativePhi(Jet->Phi(),Particle->Phi());
1504  Angularity_Numerator=Angularity_Numerator+(Particle->Pt()*TMath::Sqrt(((Particle->Eta()-Jet->Eta())*(Particle->Eta()-Jet->Eta()))+(DeltaPhi*DeltaPhi)));
1505  Angularity_Denominator= Angularity_Denominator+Particle->Pt();
1506  }
1507  if(Angularity_Denominator!=0) return Angularity_Numerator/Angularity_Denominator;
1508  else return -1;
1509 
1510 }
1511 
1512 
1513 
1514 //--------------------------------------------------------------------------
1516 
1517  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1518  Double_t PTD_Numerator=0; //Reset these values
1519  Double_t PTD_Denominator=0;
1520  AliVParticle *Particle=0x0;
1521  Double_t DeltaPhi=0;
1522  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks (particles in the jet
1523  Particle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
1524  if(!Particle) continue;
1525  PTD_Numerator=PTD_Numerator+(Particle->Pt()*Particle->Pt());
1526  PTD_Denominator=PTD_Denominator+Particle->Pt();
1527  }
1528  if(PTD_Denominator!=0) return TMath::Sqrt(PTD_Numerator)/PTD_Denominator;
1529  else return -1;
1530 
1531 }
1532 
1533 
1534 //----------------------------------------------------------------------
1536 AliEmcalJetFinder *AliAnalysisTaskSubJetFraction::Recluster(AliEmcalJet *Jet, Int_t JetContNb, Double_t SubJetRadius, Double_t SubJetMinPt, Int_t Algorithm, const char* Name){
1537 
1538  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1539  AliEmcalJetFinder *Reclusterer = new AliEmcalJetFinder(Name); //JetFinder Object for reclustered jets
1540  Reclusterer->SetRadius(SubJetRadius);
1541  Reclusterer->SetJetMinPt(SubJetMinPt);
1542  Reclusterer->SetJetAlgorithm(Algorithm); //0 for anti-kt 1 for kt
1543  Reclusterer->SetJetMaxEta(0.9);
1544  Reclusterer->SetRecombSheme(0);
1546  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
1547  Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
1548  if(Reclusterer->AliEmcalJetFinder::Filter(Jet, JetCont, dVtx)){;} //reclustering jet1 using the jetfinderobject Reclusterer
1549  }
1550  else{
1551  Double_t dVtx[3]={1,1,1};
1552  if(Reclusterer->AliEmcalJetFinder::Filter(Jet, JetCont, dVtx)){;} //reclustering jet1 using the jetfinderobject Reclusterer
1553  }
1554  return Reclusterer;
1555 
1556 }
1557 
1558 
1559 
1560 //----------------------------------------------------------------------
1562  AliEmcalJet *SubJet=NULL;
1563  Double_t SortingVariable;
1564  Int_t ArraySize =N+1;
1565  TArrayD JetSorter(ArraySize);
1566  TArrayD JetIndexSorter(ArraySize);
1567  for (Int_t i=0; i<ArraySize; i++){
1568  JetSorter[i]=0;
1569  }
1570  for (Int_t i=0; i<ArraySize; i++){
1571  JetIndexSorter[i]=0;
1572  }
1573  if(Reclusterer->GetNumberOfJets()<N) return -999;
1574  for (Int_t i=0; i<Reclusterer->GetNumberOfJets(); i++){
1575  SubJet=Reclusterer->GetJet(i);
1576  if (Type==0) SortingVariable=SubJet->Pt();
1577  else if (Type==1) SortingVariable=SubJet->E();
1578  else if (Type==2) SortingVariable=SubJet->M();
1579  for (Int_t j=0; j<N; j++){
1580  if (SortingVariable>JetSorter[j]){
1581  for (Int_t k=N-1; k>=j; k--){
1582  JetSorter[k+1]=JetSorter[k];
1583  JetIndexSorter[k+1]=JetIndexSorter[k];
1584  }
1585  JetSorter[j]=SortingVariable;
1586  JetIndexSorter[j]=i;
1587  break;
1588  }
1589  }
1590  }
1591  if (!Index) return JetSorter[N-1];
1592  else return JetIndexSorter[N-1];
1593 }
1594 
1595 
1596 
1597 //returns -1 if the Nth hardest jet is requested where N>number of available jets
1598 //type: 0=Pt 1=E 2=M
1599 //Index TRUE=returns index FALSE=returns value of quantatiy in question
1600 
1601 
1602 
1603 
1604 
1605 
1606 
1607 //----------------------------------------------------------------------
1609  AliEmcalJet *SubJet=NULL;
1610  Double_t Observable=0;
1611  Double_t Fraction_Numerator=0;
1612  Bool_t Error=kFALSE;
1613  if (!Jet->GetNumberOfTracks()) return -2;
1614  if (Add){
1615  for (Int_t i=1; i<=N; i++){
1616  Observable=SubJetOrdering(Jet,Reclusterer,i,Type,kFALSE);
1617  if(Observable==-999){
1618  Error = kTRUE;
1619  return -2;
1620  }
1621  Fraction_Numerator=Fraction_Numerator+Observable;
1622  }
1623  }
1624  else {
1625  Fraction_Numerator=SubJetOrdering(Jet,Reclusterer,N,Type,kFALSE);
1626  if (Fraction_Numerator==-999) return -2;
1627  }
1628  if (Type==0){
1629  if(Loss) return (Jet->Pt()-Fraction_Numerator)/Jet->Pt();
1630  else return Fraction_Numerator/Jet->Pt();
1631  }
1632  else if (Type==1){
1633  if(Loss) return (Jet->E()-Fraction_Numerator)/Jet->E();
1634  else return Fraction_Numerator/Jet->E();
1635  }
1636  else { //change to else if if you want to add more later
1637  if(Loss) return (Jet->M()-Fraction_Numerator)/Jet->M();
1638  else return Fraction_Numerator/Jet->M();
1639  }
1640 }
1641 //N number of hardest subjets involved
1642 //Type 0=Pt 1=E 2=M
1643 // Add TRUE: Add 1 to Nth hardest subjets together/total jet False: Nth hardest Jet/toal jet
1644 //Loss TRUE: Jet-Subjet(s)/Jet FALSE: Subjet(s)/Jet
1645 
1646 
1647 //----------------------------------------------------------------------
1649  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1650  AliEmcalJet *SubJet=NULL;
1651  Double_t DeltaR1=0;
1652  Double_t DeltaR2=0;
1653  AliVParticle *JetParticle=0x0;
1654  Double_t SubJetiness_Numerator = 0;
1655  Double_t SubJetiness_Denominator = 0;
1656  Double_t Index=-2;
1657  Bool_t Error=kFALSE;
1658  // JetRadius=TMath::Sqrt((Jet->Area()/TMath::Pi())); //comment out later
1659  if (!Jet->GetNumberOfTracks()) return -2;
1660  if (Reclusterer->GetNumberOfJets() < N) return -2;
1661  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks (particles in the jet
1662  JetParticle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
1663  for (Int_t j=1; j<=N; j++){
1664  Index=SubJetOrdering(Jet,Reclusterer,j,0,kTRUE);
1665  if(Index==-999){
1666  Error = kTRUE;
1667  i=Jet->GetNumberOfTracks();
1668  break;
1669  }
1670  if(j==1){
1671  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);
1672  }
1673  else{
1674  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);
1675  if (DeltaR2<DeltaR1) DeltaR1=DeltaR2;
1676  }
1677  }
1678  SubJetiness_Numerator=SubJetiness_Numerator+(JetParticle->Pt()*DeltaR1);
1679  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));
1680  else SubJetiness_Denominator=SubJetiness_Denominator+(TMath::Power((Reclusterer->GetJet(SubJetOrdering(Jet,Reclusterer,N,0,kTRUE))->Pt()),A)*JetParticle->Pt()*TMath::Power(JetRadius,B));
1681  }
1682  if (SubJetiness_Denominator!=0 && !Error){
1683  // if (SubJetiness_Numerator/SubJetiness_Denominator>1.0) cout << JetRadius<<endl;
1684  return SubJetiness_Numerator/SubJetiness_Denominator; }
1685  else return -2;
1686 }
1687 
1688 
1689 //______________________________________________________________________________________
1691 
1692  //WARNING!!! Only works for parent jets that are clustered with Anti-Kt! To change go to AliEmcalJetFinder.cxx and look at the Nsubjettiness() function
1693 
1694  //Algorithm==0 -> kt_axes;
1695  // Algorithm==1 -> ca_axes;
1696  //Algorithm==2 -> antikt_0p2_axes;
1697  //Algorithm==3 -> wta_kt_axes;
1698  //Algorithm==4 -> wta_ca_axes;
1699  //Algorithm==5 -> onepass_kt_axes;
1700  //Algorithm==6 -> onepass_ca_axes;
1701  //Algorithm==7 -> onepass_antikt_0p2_axes;
1702  //Algorithm==8 -> onepass_wta_kt_axes;
1703  //Algorithm==9 -> onepass_wta_ca_axes;
1704  //Algorithm==10 -> min_axes;
1705 
1706 
1707  //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.
1708 
1709 
1710  //Option==0 returns Nsubjettiness Value
1711  //Option==1 && N==2 returns opening angle between two subjet axes(Delta R?)
1712  //Option==2 && N==2 returns Delta R
1713  //Option==3 returns first splitting distance for soft dropped jet (CA)
1714  //Option==4 returns Symmetry measure (Zg) for soft dropped jet (CA)
1715  //Option==5 returns Pt of Subjet1
1716  //Option==6 returns Pt of Subjet2
1717  //Options==7 trutns deltaR of subjets...Is this different to before??
1718  //Option==8 Subjet1 Leading track Pt
1719  //Option==9 Subjet1 Leading track Pt
1720 
1721 
1722  Algorithm=fReclusteringAlgorithm;
1723  if (Jet->GetNumberOfTracks()>=N){
1724  if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==1) && (Algorithm==0) && (Beta==1.0) && (Option==0)){
1727  }
1728  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==0) && (Beta==1.0) && (Option==0)){
1731  }
1732  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==3) && (Algorithm==0) && (Beta==1.0) && (Option==0)){
1735  }
1736  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==0) && (Beta==1.0) && (Option==1)){
1739  }
1740  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==1) && (Algorithm==1) && (Beta==1.0) && (Option==0)){
1743  }
1744  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==1) && (Beta==1.0) && (Option==0)){
1747  }
1748  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==1) && (Beta==1.0) && (Option==1)){
1751  }
1752  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==1) && (Algorithm==2) && (Beta==1.0) && (Option==0)){
1755  }
1756  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==2) && (Beta==1.0) && (Option==0)){
1759  }
1760  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==2) && (Beta==1.0) && (Option==1)){
1763  }
1764  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==1) && (Algorithm==6) && (Beta==1.0) && (Option==0)){
1767  }
1768  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==6) && (Beta==1.0) && (Option==0)){
1771  }
1772  else if((fJetShapeSub==kDerivSub) && (JetContNb==0) && (N==2) && (Algorithm==6) && (Beta==1.0) && (Option==1)){
1775  }
1776  else{
1777  //Algorithm=fReclusteringAlgorithm;
1778  AliJetContainer *JetCont = GetJetContainer(JetContNb);
1779  AliEmcalJetFinder JetFinder("Nsubjettiness");
1780  JetFinder.SetJetMaxEta(0.9-fJetRadius);
1781  JetFinder.SetRadius(fJetRadius);
1782  JetFinder.SetJetAlgorithm(0); //0 for anti-kt 1 for kt //this is for the JET!!!!!!!!!! Not the SubJets
1783  JetFinder.SetRecombSheme(0);
1784  JetFinder.SetJetMinPt(Jet->Pt());
1786  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
1787  Double_t dVtx[3]={vert->GetX(),vert->GetY(),vert->GetZ()};
1788  return JetFinder.Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,fSubJetRadius,Beta,Option,0,Beta_SD,ZCut,fSoftDropOn);
1789  }
1790  else{
1791  Double_t dVtx[3]={1,1,1};
1792  if (!fNsubMeasure){
1793  return JetFinder.Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,fSubJetRadius,Beta,Option,0,Beta_SD,ZCut,fSoftDropOn);
1794  }
1795  else{
1796  if (Option==3) return JetFinder.Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,fSubJetRadius,Beta,Option,0,Beta_SD,ZCut,fSoftDropOn);
1797  else return JetFinder.Nsubjettiness(Jet,JetCont,dVtx,N,Algorithm,fSubJetRadius,Beta,Option,1);
1798  }
1799  }
1800  }
1801  }
1802  else return -2;
1803 }
1804 
1805 
1806 //________________________________________________________________________
1808  //
1809  // retrieve event objects
1810  //
1812  return kFALSE;
1813 
1814  return kTRUE;
1815 }
1816 
1817 
1818 //_______________________________________________________________________
1820 {
1821  // Called once at the end of the analysis.
1823 
1824  /*
1825  for (int i=1; i<=(XBinsJetPtSize)-1; i++){ //Rescales the fhJetPt Histograms according to the with of the variable bins and number of events
1826  fhJetPt->SetBinContent(i,((1.0/(XBinsJetPt[i]-XBinsJetPt[i-1]))*((fhJetPt->GetBinContent(i))*(1.0/(fhEventCounter->GetBinContent(1))))));
1827  fhSubJetPt->SetBinContent(i,((1.0/(XBinsJetPt[i]-XBinsJetPt[i-1]))*((fhSubJetPt->GetBinContent(i))*(1.0/(fhEventCounter->GetBinContent(5))))));
1828 
1829  //fhJetPt->SetBinContent(i,((1.0/(fhPt->GetBinWidth(i)))*((fhJetPt->GetBinContent(i))*(fhEventCounter->GetBinContent(1)))));
1830  // fhJetPt->SetBinContent(i,((1.0/((fhPt->GetLowEdge(i+1))-(fhJetPt->GetLowEdge(i))))*((fhPt->GetBinContent(i))*(fhEventCounter->GetBinContent(1)))));
1831  }
1832  for (int i=1; i<=(XBinsJetMassSize)-1; i++){ //Rescales the fhPt Histograms according to the with of the variable bins and number of events
1833  fhJetMass->SetBinContent(i,((1.0/(XBinsJetMass[i]-XBinsJetMass[i-1]))*((fhJetMass->GetBinContent(i))*(1.0/(fhEventCounter->GetBinContent(1))))));
1834  }
1835 
1836  fhJetPhi->Scale(Phi_Bins/((Phi_Up-Phi_Low)*((fhEventCounter->GetBinContent(1)))));
1837  fhJetEta->Scale(Eta_Bins/((Eta_Up-Eta_Low)*((fhEventCounter->GetBinContent(1)))));
1838  fhJetRadius->Scale(100/(fhEventCounter->GetBinContent(4))); //should this and JetAngularity be divided by Bin 1 or 4????
1839  fhNumberOfJetTracks->Scale(1.0/(fhEventCounter->GetBinContent(4)));
1840 
1841  fhJetCounter->Scale(1.0/(fhEventCounter->GetBinContent(1))); //is the first bin the correct one to look at?
1842 
1843 
1844  fhSubJetCounter->Scale(1.0/(fhEventCounter->GetBinContent(5)));
1845 
1846  */
1847  /*
1848 
1849  fhJetPt->Scale(1.0/(fhEventCounter->GetBinContent(1)));
1850  fhSubJetPt->Scale(1.0/(fhEventCounter->GetBinContent(5)));
1851  fhJetMass->Scale(1.0/(fhEventCounter->GetBinContent(1)));
1852  fhJetPhi->Scale(1.0/(fhEventCounter->GetBinContent(1)));
1853  fhJetEta->Scale(1.0/(fhEventCounter->GetBinContent(1)));
1854  fhJetRadius->Scale(1.0/(fhEventCounter->GetBinContent(4))); //should this and JetAngularity be divided by Bin 1 or 4????
1855  fhNumberOfJetTracks->Scale(1.0/(fhEventCounter->GetBinContent(4)));
1856  fhJetCounter->Scale(1.0/(fhEventCounter->GetBinContent(1))); //is the first bin the correct one to look at?
1857  fhSubJetCounter->Scale(1.0/(fhEventCounter->GetBinContent(5)));
1858  */
1859  }
1860 
1861 }
1862 
1864 std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> AliAnalysisTaskSubJetFraction::ModifyJet(AliEmcalJet* Jet, Int_t JetContNb, TString Modification){
1865  std::pair<fastjet::PseudoJet,fastjet::ClusterSequence *> Jet_ClusterSequence;
1866  std::vector<fastjet::PseudoJet> fInputVectors;
1867  AliJetContainer *fJetCont = GetJetContainer(JetContNb);
1868  AliParticleContainer *fTrackCont = fJetCont->GetParticleContainer();
1869  std::vector<fastjet::PseudoJet> Modified_Jet;
1870  fastjet::ClusterSequence *fClustSeqSA;
1871  Int_t NJetTracksTest = Jet->GetNumberOfTracks();
1872  if (fTrackCont) for (Int_t i=0; i<Jet->GetNumberOfTracks(); i++) {
1873  AliVParticle *fTrk = Jet->TrackAt(i, fTrackCont->GetArray());
1874  if (!fTrk) continue;
1875  fastjet::PseudoJet PseudoTracks(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
1876  PseudoTracks.set_user_index(Jet->TrackAt(i)+100); //100 is very arbitary....why is it here anyway?
1877  fInputVectors.push_back(PseudoTracks);
1878  }
1879  if (Modification=="Randomise") fInputVectors=RandomiseTracks(Jet,fInputVectors);
1880  if (Modification=="AddExtraProng_02_15") fInputVectors=AddExtraProng(fInputVectors,0.2,0.15);
1881  if (Modification=="AddExtraProng_02_30") fInputVectors=AddExtraProng(fInputVectors,0.2,0.30);
1882  if (Modification=="AddExtraProng_02_45") fInputVectors=AddExtraProng(fInputVectors,0.2,0.45);
1883  if (Modification=="AddExtraProng_02_60") fInputVectors=AddExtraProng(fInputVectors,0.2,0.60);
1884  if (Modification=="AddExtraProng_01_30") fInputVectors=AddExtraProng(fInputVectors,0.1,0.30);
1885  if (Modification=="AddExtraProng_01_45") fInputVectors=AddExtraProng(fInputVectors,0.1,0.45);
1886  if (Modification=="AddExtraProng_03_30") fInputVectors=AddExtraProng(fInputVectors,0.3,0.30);
1887  if (Modification=="AddExtraProng_03_45") fInputVectors=AddExtraProng(fInputVectors,0.3,0.45);
1888  fJetRadius=fJetRadius*2.0;
1889  try {
1890  fastjet::JetDefinition fJetDef(fastjet::antikt_algorithm, fJetRadius, static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
1891  fClustSeqSA=new fastjet::ClusterSequence(fInputVectors, fJetDef);
1892  Modified_Jet=fClustSeqSA->inclusive_jets(0);
1893  Jet_ClusterSequence.second=fClustSeqSA;
1894  } catch (fastjet::Error) {
1895  AliError(" [w] FJ Exception caught.");
1896  //return -1;
1897  }
1898  fJetRadius=fJetRadius/2.0;
1899  Jet_ClusterSequence.first=Modified_Jet[0];
1900  return Jet_ClusterSequence;
1901 }
1902 
1903 
1904 std::vector<fastjet::PseudoJet> AliAnalysisTaskSubJetFraction::RandomiseTracks(AliEmcalJet *Jet,std::vector<fastjet::PseudoJet> fInputVectors){
1905  TRandom3 fRandom1;
1906  fRandom1.SetSeed(0);
1907  Double_t Random_Phi=0.0;
1908  Double_t Random_Phi_Change=0.0;
1909  Double_t Random_Eta=0.0;
1910  Double_t Random_Eta_Change=0.0;
1911  Double_t Track_Pt=0.0;
1912  std::vector<fastjet::PseudoJet> Random_Track_Vector;
1913  Random_Track_Vector.clear();
1914  for (Int_t i=0; i< fInputVectors.size(); i++){
1915  Random_Phi_Change=fRandom1.Uniform(-1*fJetRadius,fJetRadius);
1916  Random_Phi=Jet->Phi()+Random_Phi_Change;
1917  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)));
1918  Random_Eta=Jet->Eta()+Random_Eta_Change;
1919  Track_Pt=fInputVectors[i].perp();
1920  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))));
1921  Random_Track.set_user_index(i);
1922  Random_Track_Vector.push_back(Random_Track);
1923  }
1924  fInputVectors.clear();
1925  return Random_Track_Vector;
1926 }
1927 
1928 std::vector<fastjet::PseudoJet> AliAnalysisTaskSubJetFraction::AddExtraProng(std::vector<fastjet::PseudoJet> fInputVectors, Double_t Distance, Double_t PtFrac){
1929  TRandom3 fRandom1;
1930  fRandom1.SetSeed(0);
1931  std::vector<fastjet::PseudoJet> Jet;
1932  Jet.clear();
1933  try {
1934  fastjet::JetDefinition fJetDef(fastjet::antikt_algorithm, fJetRadius*2, static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
1935  fastjet::ClusterSequence fClustSeqSA_Prong(fInputVectors, fJetDef);
1936  Jet= fClustSeqSA_Prong.inclusive_jets(0);
1937  } catch (fastjet::Error) {
1938  AliError(" [w] FJ Exception caught.");
1939  }
1940  Double_t Extra_Track_Phi_Change=fRandom1.Uniform(-1*Distance,Distance);
1941  Double_t Extra_Track_Phi=Jet[0].phi()+Extra_Track_Phi_Change;
1942  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)));
1943  Double_t Extra_Track_Eta=Jet[0].pseudorapidity()+Extra_Track_Eta_Change;
1944  Double_t Extra_Track_Pt=Jet[0].perp()*PtFrac;
1945  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))));
1946  ExtraTrack.set_user_index(150);
1947  fInputVectors.push_back(ExtraTrack);
1948  return fInputVectors;
1949 }
1951 
1952 //______________________________________________________________________________________
1953 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){
1954 
1955  //Only Works for kGenOnTheFly
1956 
1957  //WARNING!!! Only works for parent jets that are clustered with Anti-Kt! To change go to AliEmcalJetFinder.cxx and look at the Nsubjettiness() function
1958 
1959  //Algorithm==0 -> kt_axes;
1960  // Algorithm==1 -> ca_axes;
1961  //Algorithm==2 -> antikt_0p2_axes;
1962  //Algorithm==3 -> wta_kt_axes;
1963  //Algorithm==4 -> wta_ca_axes;
1964  //Algorithm==5 -> onepass_kt_axes;
1965  //Algorithm==6 -> onepass_ca_axes;
1966  //Algorithm==7 -> onepass_antikt_0p2_axes;
1967  //Algorithm==8 -> onepass_wta_kt_axes;
1968  //Algorithm==9 -> onepass_wta_ca_axes;
1969  //Algorithm==10 -> min_axes;
1970 
1971 
1972  //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.
1973 
1974 
1975  //Option==0 returns Nsubjettiness Value
1976  //Option==1 && N==2 returns opening angle between two subjet axes(Delta R?)
1977  //Option==2 && N==2 returns Delta R
1978  //Option==3 returns first splitting distance for soft dropped jet (CA)
1979  //Option==4 returns Symmetry measure (Zg) for soft dropped jet (CA)
1980  //Option==5 returns Pt of Subjet1
1981  //Option==6 returns Pt of Subjet2
1982  //Options==7 trutns deltaR of subjets...Is this different to before??
1983  //Option==8 Subjet1 Leading track Pt
1984  //Option==9 Subjet1 Leading track Pt
1985 
1986  fastjet::PseudoJet Jet=Jet_ClusterSequence.first;
1987  AliFJWrapper FJWrapper("FJWrapper", "FJWrapper");
1988  FJWrapper.SetR(fJetRadius);
1989  FJWrapper.SetMinJetPt(Jet.perp());
1990  FJWrapper.SetAlgorithm(fastjet::antikt_algorithm); //this is for the jet clustering not the subjet reclustering.
1991  FJWrapper.SetRecombScheme(static_cast<fastjet::RecombinationScheme>(0));
1992  if (Jet.constituents().size()>=N){
1993  Algorithm=fReclusteringAlgorithm;
1994  Double_t dVtx[3]={1,1,1};
1995  return FJWrapper.NSubjettinessDerivativeSub(N,Algorithm,fSubJetRadius,Beta,fJetRadius,Jet,Option,0,Beta_SD,ZCut,fSoftDropOn); //change this
1996  }
1997  else return -2;
1998 }
1999 
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
void Track()
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