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