AliPhysics  19b3b9d (19b3b9d)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskRecoilJetYield.cxx
Go to the documentation of this file.
1 //
2 // My analysis task
3 //
4 // Author: Harry Andrews
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"
45 #include "AliAnalysisTaskRecoilJetYield.h"
46 #include "AliMultSelection.h"
47 
48 //Globals
49 
50 
51 
52 Double_t Eta_Upper=1.00;
53 Double_t Eta_Lower=-1.00;
54 Int_t Eta_Bins_1=100;
55 Double_t Phi_Upper=2*(TMath::Pi());
56 Double_t Phi_Lower=(-1*(TMath::Pi()));
57 Int_t Phi_Bins_1=100;
58 Float_t Percentile_1 = 300;
59 
60 
61 
62 
63 using std::cout;
64 using std::endl;
65 
66 ClassImp(AliAnalysisTaskRecoilJetYield)
67 
68 //________________________________________________________________________
69 AliAnalysisTaskRecoilJetYield::AliAnalysisTaskRecoilJetYield() :
70  AliAnalysisTaskEmcalJet("AliAnalysisTaskRecoilJetYield", kTRUE),
71  fContainer(0),
72  fMinFractionShared(0),
73  fJetShapeType(kData),
74  fJetShapeSub(kNoSub),
75  fJetSelection(kInclusive),
76  fPtThreshold(-9999.),
77  fRMatching(0.2),
78  fPtMinTriggerHadron(20.),
79  fPtMaxTriggerHadron(50.),
80  fRecoilAngularWindow(0.6),
81  fSemigoodCorrect(0),
82  fHolePos(0),
83  fHoleWidth(0),
84  fCentSelectOn(kTRUE),
85  fCentMin(0),
86  fCentMax(10),
87  fSubJetAlgorithm(0),
88  fSubJetRadius(0.1),
89  fSubJetMinPt(1),
90  fJetRadius(0.4),
91  fRMatched(0.2),
92  fSharedFractionPtMin(0.5),
93  fDerivSubtrOrder(0),
94  fFullTree(kFALSE),
95  fBeta_SD(0),
96  fZCut(0.1),
97  fRho(1e-6),
98  fRhoParam(0),
99  fNsubMeasure(kFALSE),
100  fDoSoftDrop(kFALSE),
101  fhJetPt(0x0),
102  fhJetPhi(0x0),
103  fhJetEta(0x0),
104  fhJetMass(0x0),
105  fhJetRadius(0x0),
106  fhJetCounter(0x0),
107  fhNumberOfJetTracks(0x0),
108  fTreeJetInfo(0),
109  fhJetArea(0x0),
110  fhTrackPt(0x0),
111  fhGroomedPtvJetPt(0x0),
112  fhDroppedBranches(0x0),
113  fhPtTriggerHadron(0x0),
114  fh2PtTriggerHadronJet(0x0),
115  fhPhiTriggerHadronJet(0x0),
116  fhPhiTriggerHadronEventPlane(0x0),
117  fhPhiTriggerHadronEventPlaneTPC(0x0),
118  fhDetJetPt_Incl(0x0),
119  fhDetJetPt_Matched(0x0),
120  fhJetPt_Det(0x0),
121  fhJetPt_True(0x0),
122  fhJetPhi_Det(0x0),
123  fhJetPhi_True(0x0),
124  fhJetEta_Det(0x0),
125  fhJetEta_True(0x0),
126  fhJetMass_Det(0x0),
127  fhJetMass_True(0x0),
128  fhJetRadius_Det(0x0),
129  fhJetRadius_True(0x0),
130  fhJetCounter_Det(0x0),
131  fhJetCounter_True(0x0),
132  fhNumberOfJetTracks_Det(0x0),
133  fhNumberOfJetTracks_True(0x0),
134  fh2PtRatio(0x0),
135  fReclusterAlgo(0),
136  fSubMatching(kFALSE)
137 
138 {
139  for(Int_t i=0;i<nBranch;i++){
140  fJetInfoVar[i]=0;
141  }
142  SetMakeGeneralHistograms(kTRUE);
143  DefineOutput(1, TList::Class());
144  DefineOutput(2, TTree::Class());
145 }
146 
147 //________________________________________________________________________
148 AliAnalysisTaskRecoilJetYield::AliAnalysisTaskRecoilJetYield(const char *name) :
149  AliAnalysisTaskEmcalJet(name, kTRUE),
150  fContainer(0),
151  fMinFractionShared(0),
152  fJetShapeType(kData),
153  fJetShapeSub(kNoSub),
154  fJetSelection(kInclusive),
155  fPtThreshold(-9999.),
156  fRMatching(0.2),
157  fPtMinTriggerHadron(20.),
158  fPtMaxTriggerHadron(50.),
159  fRecoilAngularWindow(0.6),
160  fSemigoodCorrect(0),
161  fHolePos(0),
162  fHoleWidth(0),
163  fCentSelectOn(kTRUE),
164  fCentMin(0),
165  fCentMax(10),
166  fSubJetAlgorithm(0),
167  fSubJetRadius(0.1),
168  fSubJetMinPt(1),
169  fJetRadius(0.4),
170  fRMatched(0.2),
171  fSharedFractionPtMin(0.5),
172  fDerivSubtrOrder(0),
173  fFullTree(kFALSE),
174  fBeta_SD(0),
175  fZCut(0.1),
176  fRho(1e-6),
177  fRhoParam(0),
178  fNsubMeasure(kFALSE),
179  fDoSoftDrop(kFALSE),
180  fhJetPt(0x0),
181  fhJetPhi(0x0),
182  fhJetEta(0x0),
183  fhJetMass(0x0),
184  fhJetRadius(0x0),
185  fhJetCounter(0x0),
186  fhNumberOfJetTracks(0x0),
187  fTreeJetInfo(0),
188  fhJetArea(0x0),
189  fhTrackPt(0x0),
190  fhGroomedPtvJetPt(0x0),
191  fhDroppedBranches(0x0),
192  fhPtTriggerHadron(0x0),
193  fh2PtTriggerHadronJet(0x0),
194  fhPhiTriggerHadronJet(0x0),
195  fhPhiTriggerHadronEventPlane(0x0),
196  fhPhiTriggerHadronEventPlaneTPC(0x0),
197  fhDetJetPt_Incl(0x0),
198  fhDetJetPt_Matched(0x0),
199  fhJetPt_Det(0x0),
200  fhJetPt_True(0x0),
201  fhJetPhi_Det(0x0),
202  fhJetPhi_True(0x0),
203  fhJetEta_Det(0x0),
204  fhJetEta_True(0x0),
205  fhJetMass_Det(0x0),
206  fhJetMass_True(0x0),
207  fhJetRadius_Det(0x0),
208  fhJetRadius_True(0x0),
209  fhJetCounter_Det(0x0),
210  fhJetCounter_True(0x0),
211  fhNumberOfJetTracks_Det(0x0),
212  fhNumberOfJetTracks_True(0x0),
213  fh2PtRatio(0x0),
214  fReclusterAlgo(0),
215  fSubMatching(kFALSE)
216 
217 {
218  // Standard constructor.
219  for(Int_t i=0;i<nBranch;i++){
220  fJetInfoVar[i]=0;
221  }
222  SetMakeGeneralHistograms(kTRUE);
223  DefineOutput(1, TList::Class());
224  DefineOutput(2, TTree::Class());
225 }
226 
227 //________________________________________________________________________
228 AliAnalysisTaskRecoilJetYield::~AliAnalysisTaskRecoilJetYield()
229 {
230  // Destructor.
231 }
232 
233 //________________________________________________________________________
234  void AliAnalysisTaskRecoilJetYield::UserCreateOutputObjects()
235 {
236  // Create user output.
237 
238  AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
239 
240  Bool_t oldStatus = TH1::AddDirectoryStatus();
241  TH1::AddDirectory(kFALSE);
242  TH1::AddDirectory(oldStatus);
243  const char* nameoutput = GetOutputSlot(2)->GetContainer()->GetName();
244  fTreeJetInfo = new TTree(nameoutput, nameoutput);
245 
246  if (!fFullTree){
247  const Int_t nVarMin = 18;
248  TString *fJetInfoVarNames = new TString [nVarMin];
249 
250  fJetInfoVarNames[0] = "Pt";
251  fJetInfoVarNames[1] = "Pt_Truth";
252  fJetInfoVarNames[2] = "SymParam";
253  fJetInfoVarNames[3] = "SymParam_Truth";
254  fJetInfoVarNames[4] = "Tau1";
255  fJetInfoVarNames[5] = "Tau1_Truth";
256  fJetInfoVarNames[6] = "Tau2";
257  fJetInfoVarNames[7] = "Tau2_Truth";
258  fJetInfoVarNames[8] = "PTD";
259  fJetInfoVarNames[9] = "PTD_Truth";
260  fJetInfoVarNames[10] = "Angularity";
261  fJetInfoVarNames[11] = "Angularity_Truth";
262  fJetInfoVarNames[12] = "DelR";
263  fJetInfoVarNames[13] = "DelR_Truth";
264  fJetInfoVarNames[14] = "N_Groomed_Branches";
265  fJetInfoVarNames[15] = "N_Groomed_Branches_Truth";
266  fJetInfoVarNames[16] = "Groomed_Jet_Pt";
267  fJetInfoVarNames[17] = "Groomed_Jet_Pt_Truth";
268 
269 
270 
271  for(Int_t ivar=0; ivar < nVarMin; ivar++){
272  cout<<"looping over variables"<<endl;
273  fTreeJetInfo->Branch(fJetInfoVarNames[ivar].Data(), &fJetInfoVar[ivar], Form("%s/D", fJetInfoVarNames[ivar].Data()));
274  }
275  }
276  if(fJetSelection == kRecoil){
277  fhPtTriggerHadron= new TH1F("fhPtTriggerHadron", "fhPtTriggerHadron",1500,-0.5,149.5);
278  fOutput->Add(fhPtTriggerHadron);
279  fh2PtTriggerHadronJet= new TH2F("fh2PtTriggerHadronJet", "fh2PtTriggerHadronJet",1500,-0.5,149.5,1500,-0.5,149.5);
280  fOutput->Add(fh2PtTriggerHadronJet);
281  fhPhiTriggerHadronJet= new TH1F("fhPhiTriggerHadronJet", "fhPhiTriggerHadronJet",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
282  fOutput->Add(fhPhiTriggerHadronJet);
283  fhPhiTriggerHadronEventPlane= new TH1F("fhPhiTriggerHadronEventPlane", "fhPhiTriggerHadronEventPlane",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
284  fOutput->Add(fhPhiTriggerHadronEventPlane);
285  fhPhiTriggerHadronEventPlaneTPC= new TH1F("fhPhiTriggerHadronEventPlaneTPC", "fhPhiTriggerHadronEventPlaneTPC",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
286  fOutput->Add(fhPhiTriggerHadronEventPlaneTPC);
287 
288 
289  }
290 
291  if (fJetShapeType==AliAnalysisTaskRecoilJetYield::kData || fJetShapeType==AliAnalysisTaskRecoilJetYield::kDetEmbPart || fJetShapeType==AliAnalysisTaskRecoilJetYield::kGenOnTheFly){
292 
293  fhJetPt= new TH1F("fhJetPt", "Jet Pt",150,-0.5,149.5 );
294  fOutput->Add(fhJetPt);
295  fhJetPhi= new TH1F("fhJetPhi", "Jet Phi",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
296  fOutput->Add(fhJetPhi);
297  fhJetEta= new TH1F("fhJetEta", "Jet Eta", Eta_Bins_1, Eta_Lower, Eta_Upper);
298  fOutput->Add(fhJetEta);
299  fhJetMass= new TH1F("fhJetMass", "Jet Mass", 4000,-0.5, 39.5);
300  fOutput->Add(fhJetMass);
301  fhJetRadius= new TH1F("fhJetRadius", "Jet Radius", 100, -0.05,0.995);
302  fOutput->Add(fhJetRadius);
303  fhNumberOfJetTracks= new TH1F("fhNumberOfJetTracks", "Number of Tracks within a Jet", 300, -0.5,299.5);
304  fOutput->Add(fhNumberOfJetTracks);
305  fhJetCounter= new TH1F("fhJetCounter", "Jet Counter", 150, -0.5, 149.5);
306  fOutput->Add(fhJetCounter);
307  fhJetArea= new TH1F("fhJetArea", "Jet Area", 400,-0.5, 1.95);
308  fOutput->Add(fhJetArea);
309  fhTrackPt= new TH1F("fhTrackPt", "Track Pt",600,-0.5,59.5);
310  fOutput->Add(fhTrackPt);
311  fhGroomedPtvJetPt= new TH2F("fhGroomedPtvJetPt","Groomed Jet p_{T} v Original Jet p_{T}",150,0,150,150,0,150);
312  fOutput->Add(fhGroomedPtvJetPt);
313  fhDroppedBranches= new TH1F("fhDroppedBranches","Number of Softdropped branches",50,0,50);
314  fOutput->Add(fhDroppedBranches);
315  fhDetJetPt_Incl= new TH1F("fhDetJetPt_Incl", "Jet Pt",200,-0.5,199.5 );
316  fOutput->Add(fhDetJetPt_Incl);
317  fhDetJetPt_Matched= new TH1F("fhDetJetPt_Matched", "Jet Pt",200,-0.5,199.5 );
318  fOutput->Add(fhDetJetPt_Matched);
319 
320 
321  }
322  if(fJetShapeType==AliAnalysisTaskRecoilJetYield::kTrueDet){
323  fhJetPt_Det= new TH1F("fhJetPt_Det", "Jet Pt Detector Level",1500,-0.5,149.5 );
324  fOutput->Add(fhJetPt_Det);
325  fhJetPt_True= new TH1F("fhJetPt_True", "Jet Pt Particle Level",1500,-0.5,149.5 );
326  fOutput->Add(fhJetPt_True);
327  fhJetPhi_Det= new TH1F("fhJetPhi_Det", "Jet Phi Detector Level",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
328  fOutput->Add(fhJetPhi_Det);
329  fhJetPhi_True= new TH1F("fhJetPhi_True", "Jet Phi Particle Level",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
330  fOutput->Add(fhJetPhi_True);
331  fhJetEta_Det= new TH1F("fhJetEta_Det", "Jet Eta Detector Level", Eta_Bins_1, Eta_Lower, Eta_Upper);
332  fOutput->Add(fhJetEta_Det);
333  fhJetEta_True= new TH1F("fhJetEta_True", "Jet Eta Particle Level", Eta_Bins_1, Eta_Lower, Eta_Upper);
334  fOutput->Add(fhJetEta_True);
335  fhJetMass_Det= new TH1F("fhJetMass_Det", "Jet Mass Detector Level", 4000,-0.5, 39.5);
336  fOutput->Add(fhJetMass_Det);
337  fhJetMass_True= new TH1F("fhJetMass_True", "Jet Mass Particle Level", 4000,-0.5, 39.5);
338  fOutput->Add(fhJetMass_True);
339  fhJetRadius_Det= new TH1F("fhJetRadius_Det", "Jet Radius Detector Level", 100, -0.05,0.995);
340  fOutput->Add(fhJetRadius_Det);
341  fhJetRadius_True= new TH1F("fhJetRadius_True", "Jet Radius Particle Level", 100, -0.05,0.995);
342  fOutput->Add(fhJetRadius_True);
343  fhNumberOfJetTracks_Det= new TH1F("fhNumberOfJetTracks_Det", "Number of Tracks within a Jet Detector Level", 300, -0.5,299.5);
344  fOutput->Add(fhNumberOfJetTracks_Det);
345  fhNumberOfJetTracks_True= new TH1F("fhNumberOfJetTracks_True", "Number of Tracks within a Jet Particle Level", 300, -0.5,299.5);
346  fOutput->Add(fhNumberOfJetTracks_True);
347 
348  fhJetCounter_Det= new TH1F("fhJetCounter_Det", "Jet Counter Detector Level", 150, -0.5, 149.5);
349  fOutput->Add(fhJetCounter_Det);
350  fhJetCounter_True= new TH1F("fhJetCounter_True", "Jet Counter Particle Level", 150, -0.5, 149.5);
351  fOutput->Add(fhJetCounter_True);
352  fh2PtRatio= new TH2F("fhPtRatio", "MC pt for detector level vs particle level jets",1500,-0.5,149.5,1500,-0.5,149.5);
353  fOutput->Add(fh2PtRatio);
354  }
355  PostData(1,fOutput);
356  PostData(2,fTreeJetInfo);
357  //cout<<"End of UserCreateOutputObjects"<<endl;
358  // delete [] fShapesVarNames;
359 }
360 
361 //________________________________________________________________________
362 Bool_t AliAnalysisTaskRecoilJetYield::Run()
363 {
364  // Run analysis code here, if needed. It will be executed before FillHistograms().
365  return kTRUE;
366 }
367 
368 //________________________________________________________________________
369 Bool_t AliAnalysisTaskRecoilJetYield::FillHistograms()
370 {
371  //AliMultSelection *MultSelection = 0x0;
372  AliAODEvent *fEvent = dynamic_cast<AliAODEvent*>(InputEvent());
373  if (!fEvent) {
374  Error("UserExec","AOD not available");
375  return 0;
376  }
377  /*MultSelection = (AliMultSelection * ) fEvent->FindListObject("MultSelection");
378  if(!MultSelection) {
379  AliWarning("AliMultSelection object not found!");
380  }
381  else{
382  Percentile_1 = MultSelection->GetMultiplicityPercentile("V0M");
383  }
384  */
385  if (fCentSelectOn){
386  if ((fCent>fCentMax) || (fCent<fCentMin)) return 0;
387  }
388 
390  AliAODTrack *TriggerHadron = 0x0;
391  if (fJetSelection == kRecoil) {
392  //you have to set a flag and the limits of the pT interval for your trigger
393  Int_t TriggerHadronLabel = SelectTriggerHadron(fPtMinTriggerHadron, fPtMaxTriggerHadron);
394  if (TriggerHadronLabel==-99999) return 0; //Trigger Hadron Not Found
395  AliTrackContainer *PartCont =NULL;
396  AliParticleContainer *PartContMC = NULL;
397  if (fJetShapeSub==kConstSub){
398  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kGenOnTheFly) PartContMC = GetParticleContainer(1);
399  else PartCont = GetTrackContainer(1);
400  }
401  else {
402  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kGenOnTheFly) PartContMC = GetParticleContainer(0);
403  else PartCont = GetTrackContainer(0);
404  }
405  TClonesArray *TrackArray = NULL;
406  TClonesArray *TrackArrayMC = NULL;
407  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kGenOnTheFly) TrackArrayMC = PartContMC->GetArray();
408  else TrackArray = PartCont->GetArray();
409  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kGenOnTheFly) TriggerHadron = static_cast<AliAODTrack*>(TrackArrayMC->At(TriggerHadronLabel));
410  else TriggerHadron = static_cast<AliAODTrack*>(TrackArray->At(TriggerHadronLabel));
411  if (!TriggerHadron) return 0;//No trigger hadron with label found
412  if(fSemigoodCorrect){
413  Double_t HoleDistance=RelativePhi(TriggerHadron->Phi(),fHolePos);
414  if(TMath::Abs(HoleDistance)+fHoleWidth+fJetRadius>TMath::Pi()-fRecoilAngularWindow) return 0;
415  }
416  fhPtTriggerHadron->Fill(TriggerHadron->Pt()); //Needed for per trigger Normalisation
417  if (fJetShapeType != AliAnalysisTaskRecoilJetYield::kGenOnTheFly) fhPhiTriggerHadronEventPlane->Fill(TMath::Abs(RelativePhiEventPlane(fEPV0,TriggerHadron->Phi()))); //fEPV0 is the event plane from AliAnalysisTaskEmcal
418  if (fJetShapeType != AliAnalysisTaskRecoilJetYield::kGenOnTheFly) fhPhiTriggerHadronEventPlaneTPC->Fill(TMath::Abs(RelativePhiEventPlane(((AliVAODHeader*)(InputEvent()->GetHeader()))->GetEventplane(),TriggerHadron->Phi()))); //TPC event plane
419  }
420 
421 
422 
423 
424 
426  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kData){
427  // cout<<"entering kData"<<endl;
428 
429  AliEmcalJet *Jet1 = NULL; //Original Jet in the event
430  AliJetContainer *JetCont= GetJetContainer(0); //Jet Container for event
431  Int_t JetCounter=0; //Counts number of jets in event
432  Double_t JetPhi=0;
433  Bool_t EventCounter=kFALSE;
434  Double_t JetPt_ForThreshold=0;
435  AliEmcalJet *Jet2= NULL;
436  if(JetCont) {
437  JetCont->ResetCurrentID();
438  while((Jet1=JetCont->GetNextAcceptJet())) {
439  if(!Jet1) continue;
440  if (fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) JetPt_ForThreshold = Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
441  else JetPt_ForThreshold = Jet1->Pt();
442  if(JetPt_ForThreshold<fPtThreshold) {
443  continue;
444  }
445  else {
446  Float_t RecoilDeltaPhi = 0.;
447  if (fJetSelection == kRecoil){
448  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), Jet1->Phi());
449  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
450  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), Jet1->Pt());
451  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), Jet1->Phi()));
452  }
453  if (!EventCounter){
454  EventCounter=kTRUE;
455  }
456  JetCounter++;
457  fhJetPt->Fill(Jet1->Pt());
458  fhJetArea->Fill(Jet1->Area());
459  JetPhi=Jet1->Phi();
460  if(JetPhi < -1*TMath::Pi()) JetPhi += (2*TMath::Pi());
461  else if (JetPhi > TMath::Pi()) JetPhi -= (2*TMath::Pi());
462  fhJetPhi->Fill(JetPhi);
463  fhJetEta->Fill(Jet1->Eta());
464  fhJetMass->Fill(Jet1->M());
465  fhJetRadius->Fill(TMath::Sqrt((Jet1->Area()/TMath::Pi()))); //Radius of Jets per event
466  fhNumberOfJetTracks->Fill(Jet1->GetNumberOfTracks());
467  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) fJetInfoVar[0]= Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
468  else fJetInfoVar[0]=Jet1->Pt();
469  fJetInfoVar[1]=0;
470  if(fDoSoftDrop) {
471  SoftDrop(Jet1,JetCont,fZCut,fBeta_SD,kFALSE);
472  //SoftDrop(Jet1,JetCont,fZCut,fBeta_SD,kTRUE);
473 
474  fJetInfoVar[3]=0;
475  fJetInfoVar[5]=0;
476  fJetInfoVar[7]=0;
477  fJetInfoVar[13]=0;
478  fJetInfoVar[15]=0;
479  fJetInfoVar[17]=0;
480 
481  }
482  else{
483  fJetInfoVar[2]=0;
484  fJetInfoVar[3]=0;
485  fJetInfoVar[4]=0;
486  fJetInfoVar[5]=0;
487  fJetInfoVar[6]=0;
488  fJetInfoVar[7]=0;
489  fJetInfoVar[12]=0;
490  fJetInfoVar[13]=0;
491  fJetInfoVar[14]=0;
492  fJetInfoVar[15]=0;
493  fJetInfoVar[16]=0;
494  fJetInfoVar[17]=0;
495  }
496  fJetInfoVar[8]=PTD(Jet1,0);
497  fJetInfoVar[9]=0;
498  fJetInfoVar[10]=Angularity(Jet1,0);
499  fJetInfoVar[11]=0;
500  fTreeJetInfo->Fill();
501  }
502  }
503  fhJetCounter->Fill(JetCounter); //Number of Jets in Each Event
504  }
505  }
506 
508  if(fJetShapeType == kDetEmbPart){
509  AliEmcalJet *JetHybridS = NULL; //Subtracted hybrid Jet
510  AliEmcalJet *JetHybridUS = NULL; //Unsubtracted Hybrid Jet //For matching SubtractedHybrid->DetPythia this jet container is also Subtracted Hybrid
511  AliEmcalJet *JetPythDet = NULL; //Detector Level Pythia Jet
512  AliEmcalJet *JetPythTrue = NULL; //Particle Level Pyhtia Jet
513  AliJetContainer *JetContHybridS= GetJetContainer(0); //Jet Container for Subtracted Hybrid Jets
514  AliJetContainer *JetContHybridUS= GetJetContainer(1); //Jet Container for Unsubtracted Hybrid Jets
515  AliJetContainer *JetContPythDet= GetJetContainer(2); //Jet Container for Detector Level Pyhtia Jets
516  AliJetContainer *JetContPythTrue= GetJetContainer(3); //Jet Container for Particle Level Pythia Jets
517 
518 
519 
520  Bool_t JetsMatched = kFALSE;
521  Double_t JetPtThreshold;
522  JetContHybridS->ResetCurrentID();
523  JetContHybridUS->ResetCurrentID();
524  JetContPythDet->ResetCurrentID();
525  JetContPythTrue->ResetCurrentID();
526 
527  while((JetPythDet = JetContPythDet->GetNextAcceptJet())){
528  fhDetJetPt_Incl->Fill(JetPythDet->Pt()); //Fill histogram with all Detector level jets
529  }
530 
531 
532  while((JetHybridS = JetContHybridS->GetNextAcceptJet())){
533  if (fJetShapeSub==kConstSub) JetPtThreshold=JetHybridS->Pt();
534  else JetPtThreshold=JetHybridS->Pt()-(GetRhoVal(0)*JetHybridS->Area());
535  if ( (!JetHybridS) || (JetPtThreshold<fPtThreshold)) continue;
536  Int_t JetNumber=-1;
537  for(Int_t i = 0; i<JetContHybridUS->GetNJets(); i++) {
538  JetHybridUS = JetContHybridUS->GetJet(i);
539  if (!JetHybridUS) continue;
540 
541  if(JetHybridUS->GetLabel()==JetHybridS->GetLabel()) {
542  JetNumber=i;
543  }
544  }
545  if(JetNumber==-1) continue;
546  JetHybridUS=JetContHybridUS->GetJet(JetNumber);
547  if (JetContHybridUS->AliJetContainer::GetFractionSharedPt(JetHybridUS)<fSharedFractionPtMin && !fSubMatching) {
548  continue;
549  }
550  else if(GetFractionSharedPt_SubMatching(JetHybridUS,JetContHybridUS->GetParticleContainer())<fSharedFractionPtMin && fSubMatching){
551  continue;
552  }
553  JetPythDet=JetHybridUS->ClosestJet();
554  if (!JetHybridUS) {
555  Printf("Unsubtracted embedded jet does not exist, returning");
556  continue;
557  }
558  if (!JetPythDet) continue;
559  UInt_t rejectionReason = 0;
560  if (!(JetContPythDet->AcceptJet(JetPythDet,rejectionReason))) continue;
561  fhDetJetPt_Matched->Fill(JetPythDet->Pt()); //Fill only matched detector level jets for tagging efficiency comparison
562  JetPythTrue=JetPythDet->ClosestJet();
563  if(!JetPythTrue) continue;
564  JetsMatched=kTRUE;
565 
566  if(fJetShapeSub==kConstSub) fJetInfoVar[0]=JetHybridS->Pt();
567  else fJetInfoVar[0]=JetHybridS->Pt()-(GetRhoVal(0)*JetHybridS->Area());
568  fJetInfoVar[1]=JetPythTrue->Pt();
569  if(fDoSoftDrop) {
570  SoftDrop(JetHybridS,JetContHybridS,fZCut,fBeta_SD,kFALSE);
571  SoftDrop(JetPythTrue,JetContPythTrue,fZCut,fBeta_SD,kTRUE);
572  }
573  else{
574  fJetInfoVar[2]=0;
575  fJetInfoVar[3]=0;
576  fJetInfoVar[4]=0;
577  fJetInfoVar[5]=0;
578  fJetInfoVar[6]=0;
579  fJetInfoVar[7]=0;
580  fJetInfoVar[12]=0;
581  fJetInfoVar[13]=0;
582  fJetInfoVar[14]=0;
583  fJetInfoVar[15]=0;
584  fJetInfoVar[16]=0;
585  fJetInfoVar[17]=0;
586  }
587  fJetInfoVar[8]=PTD(JetHybridS,0);
588  fJetInfoVar[9]=PTD(JetPythTrue,0);
589  fJetInfoVar[10]=Angularity(JetHybridS,0);
590  fJetInfoVar[11]=Angularity(JetPythTrue,0);
591  fTreeJetInfo->Fill();
592  }
593 
594  }
595 
597  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kTrueDet){ //Truth->Detector response
598  AliEmcalJet *JetDet = NULL; //Detector Level Jet
599  AliEmcalJet *JetTrue = NULL; //Particle Level Jet
600  AliJetContainer *JetContDet= GetJetContainer(0); //Jet Container for Detector Level Pythia
601  AliJetContainer *JetContTrue= GetJetContainer(1); //Jet Container for Particle Level Pythia
602  Double_t JetPhiDet=0;
603  Double_t JetPhiTrue=0;
604  Bool_t JetsMatched=kFALSE;
605  Double_t Pythia_Event_Weight=1;
606  Bool_t EventCounter=kFALSE;
607  Int_t JetCounter_Det=0,JetCounter_True=0;
608  JetContDet->ResetCurrentID();
609  JetContTrue->ResetCurrentID();
610  if(JetContDet) {
611  while((JetDet=JetContDet->GetNextAcceptJet())) {
612  if( (!JetDet) || ((JetDet->Pt())<fPtThreshold)) {
613  // fhEventCounter_1->Fill(3); //events where the jet had a problem
614  continue;
615  }
616  else {
617  /* if(fSemigoodCorrect){
618  Double_t HoleDistance=RelativePhi(Jet1->Phi(),fHolePos);
619  if(TMath::Abs(HoleDistance)<fHoleWidth) continue;
620  }*/
621  Float_t RecoilDeltaPhi = 0.;
622  if (fJetSelection == kRecoil){
623  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), JetDet->Phi());
624  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
625  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), JetDet->Pt());
626  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), JetDet->Phi()));
627  }
628  JetCounter_Det++;
629  fhJetPt_Det->Fill(JetDet->Pt());
630  JetPhiDet=JetDet->Phi();
631  if(JetPhiDet < -1*TMath::Pi()) JetPhiDet += (2*TMath::Pi());
632  else if (JetPhiDet > TMath::Pi()) JetPhiDet -= (2*TMath::Pi());
633  fhJetPhi_Det->Fill(JetPhiDet);
634  fhJetEta_Det->Fill(JetDet->Eta());
635  fhJetMass_Det->Fill(JetDet->M());
636  fhJetRadius_Det->Fill(TMath::Sqrt((JetDet->Area()/TMath::Pi())));
637  fhNumberOfJetTracks_Det->Fill(JetDet->GetNumberOfTracks());
638  if((JetTrue = JetDet->ClosestJet())){
639  JetsMatched=kTRUE;
640  JetCounter_True++;
641  fhJetPt_True->Fill(JetTrue->Pt());
642  JetPhiTrue=JetTrue->Phi();
643  if(JetPhiTrue < -1*TMath::Pi()) JetPhiTrue += (2*TMath::Pi());
644  else if (JetPhiTrue > TMath::Pi()) JetPhiTrue -= (2*TMath::Pi());
645  fhJetPhi_True->Fill(JetPhiTrue);
646  fhJetEta_True->Fill(JetTrue->Eta());
647  fhJetMass_True->Fill(JetTrue->M());
648  fhJetRadius_True->Fill(TMath::Sqrt((JetTrue->Area()/TMath::Pi())));
649  fhNumberOfJetTracks_True->Fill(JetTrue->GetNumberOfTracks());
650  fh2PtRatio->Fill(JetDet->Pt(),JetTrue->Pt(),Pythia_Event_Weight);
651  }
652  else continue;
653  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) fJetInfoVar[0]= JetDet->Pt()-(GetRhoVal(0)*JetDet->Area());
654  else fJetInfoVar[0]=JetDet->Pt();
655  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) fJetInfoVar[1]= JetTrue->Pt()-(GetRhoVal(0)*JetTrue->Area());
656  else fJetInfoVar[1]=JetTrue->Pt();
657  if(fDoSoftDrop) {
658  SoftDrop(JetDet,JetContDet,fZCut,fBeta_SD,kFALSE);
659  SoftDrop(JetTrue,JetContTrue,fZCut,fBeta_SD,kTRUE);
660  }
661  else{
662  fJetInfoVar[2]=0;
663  fJetInfoVar[3]=0;
664  fJetInfoVar[4]=0;
665  fJetInfoVar[5]=0;
666  fJetInfoVar[6]=0;
667  fJetInfoVar[7]=0;
668  fJetInfoVar[12]=0;
669  fJetInfoVar[13]=0;
670  fJetInfoVar[14]=0;
671  fJetInfoVar[15]=0;
672  fJetInfoVar[16]=0;
673  fJetInfoVar[17]=0;
674  }
675  fJetInfoVar[8]=PTD(JetDet,0);
676  fJetInfoVar[9]=PTD(JetTrue,0);
677  fJetInfoVar[10]=Angularity(JetDet,0);
678  fJetInfoVar[11]=Angularity(JetTrue,0);
679  fTreeJetInfo->Fill();
680 
681  JetsMatched=kFALSE;
682  }
683  }
684  fhJetCounter_Det->Fill(JetCounter_Det); //Number of Jets in Each Event Particle Level
685  fhJetCounter_True->Fill(JetCounter_True); //Number of Jets in Each Event Detector Level
686  }
687  }
688 
689 
691  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kGenOnTheFly){
692  AliEmcalJet *Jet1 = NULL; //Original Jet in the event
693  AliJetContainer *JetCont= GetJetContainer(0); //Jet Container for event
694  Int_t JetCounter=0; //Counts number of jets in event
695  Double_t JetPhi=0;
696  Bool_t EventCounter=kFALSE;
697  Double_t JetPt_ForThreshold=0;
698  AliEmcalJet *Jet2= NULL;
699  if(JetCont) {
700  JetCont->ResetCurrentID();
701  while((Jet1=JetCont->GetNextAcceptJet())) {
702  if(!Jet1) continue;
703  if (fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) JetPt_ForThreshold = Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
704  else JetPt_ForThreshold = Jet1->Pt();
705  if(JetPt_ForThreshold<fPtThreshold) {
706  continue;
707  }
708  else {
709  Float_t RecoilDeltaPhi = 0.;
710  if (fJetSelection == kRecoil){
711  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), Jet1->Phi());
712  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
713  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), Jet1->Pt());
714  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), Jet1->Phi()));
715  }
716  if (!EventCounter){
717  EventCounter=kTRUE;
718  }
719  JetCounter++;
720  fhJetPt->Fill(Jet1->Pt());
721  fhJetArea->Fill(Jet1->Area());
722  JetPhi=Jet1->Phi();
723  if(JetPhi < -1*TMath::Pi()) JetPhi += (2*TMath::Pi());
724  else if (JetPhi > TMath::Pi()) JetPhi -= (2*TMath::Pi());
725  fhJetPhi->Fill(JetPhi);
726  fhJetEta->Fill(Jet1->Eta());
727  fhJetMass->Fill(Jet1->M());
728  fhJetRadius->Fill(TMath::Sqrt((Jet1->Area()/TMath::Pi()))); //Radius of Jets per event
729  fhNumberOfJetTracks->Fill(Jet1->GetNumberOfTracks());
730  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) fJetInfoVar[0]= Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
731  else fJetInfoVar[0]=Jet1->Pt();
732  fJetInfoVar[1]=0;
733  if(fDoSoftDrop) {
734  SoftDrop(Jet1,JetCont,fZCut,fBeta_SD,kFALSE);
735  SoftDrop(Jet1,JetCont,fZCut,fBeta_SD,kTRUE);
736  }
737  else{
738  fJetInfoVar[2]=0;
739  fJetInfoVar[3]=0;
740  fJetInfoVar[4]=0;
741  fJetInfoVar[5]=0;
742  fJetInfoVar[6]=0;
743  fJetInfoVar[7]=0;
744  fJetInfoVar[12]=0;
745  fJetInfoVar[13]=0;
746  fJetInfoVar[14]=0;
747  fJetInfoVar[15]=0;
748  fJetInfoVar[16]=0;
749  fJetInfoVar[17]=0;
750  }
751  fJetInfoVar[8]=PTD(Jet1,0);
752  fJetInfoVar[9]=0;
753  fJetInfoVar[10]=Angularity(Jet1,0);
754  fJetInfoVar[11]=0;
755  fTreeJetInfo->Fill();
756  }
757  }
758  fhJetCounter->Fill(JetCounter); //Number of Jets in Each Event
759  }
760  }
761 
762 
763 
764  return kTRUE;
765 }
766 
767 //________________________________________________________________________
768 Double_t AliAnalysisTaskRecoilJetYield::RelativePhiEventPlane(Double_t EventPlane, Double_t Phi){
769 
770  if(Phi < -1*TMath::Pi()) Phi += (2*TMath::Pi());
771  else if (Phi > TMath::Pi()) Phi -= (2*TMath::Pi());
772  Double_t DeltaPhi=Phi-EventPlane;
773  if(DeltaPhi < -1*TMath::Pi()) DeltaPhi += (2*TMath::Pi());
774  else if (DeltaPhi > TMath::Pi()) DeltaPhi -= (2*TMath::Pi());
775  return DeltaPhi;
776 }
777 //________________________________________________________________________
778 Double_t AliAnalysisTaskRecoilJetYield::RelativePhi(Double_t Phi1, Double_t Phi2){
779 
780  if(Phi1 < -1*TMath::Pi()) Phi1 += (2*TMath::Pi());
781  else if (Phi1 > TMath::Pi()) Phi1 -= (2*TMath::Pi());
782  if(Phi2 < -1*TMath::Pi()) Phi2 += (2*TMath::Pi());
783  else if (Phi2 > TMath::Pi()) Phi2 -= (2*TMath::Pi());
784  Double_t DeltaPhi=Phi2-Phi1;
785  if(DeltaPhi < -1*TMath::Pi()) DeltaPhi += (2*TMath::Pi());
786  else if (DeltaPhi > TMath::Pi()) DeltaPhi -= (2*TMath::Pi());
787  return DeltaPhi;
788 }
789 
790 
791 
792 
793 //--------------------------------------------------------------------------
794 Double_t AliAnalysisTaskRecoilJetYield::Angularity(AliEmcalJet *Jet, Int_t JetContNb){
795 
796  AliJetContainer *JetCont = GetJetContainer(JetContNb);
797  Double_t Angularity_Numerator=0;
798  Double_t Angularity_Denominator=0;
799  AliVParticle *Particle=0x0;
800  Double_t DeltaPhi=0;
801 
802 
803  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks
804  Particle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
805  if(!Particle) continue;
806  DeltaPhi=RelativePhi(Jet->Phi(),Particle->Phi());
807  Angularity_Numerator=Angularity_Numerator+(Particle->Pt()*TMath::Sqrt(((Particle->Eta()-Jet->Eta())*(Particle->Eta()-Jet->Eta()))+(DeltaPhi*DeltaPhi)));
808  Angularity_Denominator= Angularity_Denominator+Particle->Pt();
809  }
810  if(Angularity_Denominator!=0) return Angularity_Numerator/Angularity_Denominator;
811  else return -1;
812 
813 }
814 
815 
816 
817 //--------------------------------------------------------------------------
818 Double_t AliAnalysisTaskRecoilJetYield::PTD(AliEmcalJet *Jet, Int_t JetContNb){
819 
820  AliJetContainer *JetCont = GetJetContainer(JetContNb);
821  Double_t PTD_Numerator=0; //Reset these values
822  Double_t PTD_Denominator=0;
823  AliVParticle *Particle=0x0;
824  Double_t DeltaPhi=0;
825  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks
826  Particle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
827  if(!Particle) continue;
828  PTD_Numerator=PTD_Numerator+(Particle->Pt()*Particle->Pt());
829  PTD_Denominator=PTD_Denominator+Particle->Pt();
830  }
831  if(PTD_Denominator!=0) return TMath::Sqrt(PTD_Numerator)/PTD_Denominator;
832  else return -1;
833 
834 }
835 
836 //_________________________________________________________________________
837  void AliAnalysisTaskRecoilJetYield::SoftDrop(AliEmcalJet *fJet,AliJetContainer *fJetCont, double zcut, double beta, Bool_t fTruthJet){
838  std::vector<fastjet::PseudoJet> fInputVectors;
839  Double_t JetInvMass=0, PseudJetInvMass=0, TrackMom = 0, TrackEnergy = 0;
840  AliParticleContainer *fTrackCont = fJetCont->GetParticleContainer();
841  Double_t JetEta=fJet->Eta(),JetPhi=fJet->Phi();
842  Double_t FJTrackEta[9999],FJTrackPhi[9999],FJTrackPt[9999],EmcalJetTrackEta[9999],EmcalJetTrackPhi[9999],EmcalJetTrackPt[9999];
843  UShort_t FJNTracks=0,EmcalJetNTracks=0;
844  if (fTrackCont) for (Int_t i=0; i<fJet->GetNumberOfTracks(); i++) {
845  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
846  JetInvMass += fTrk->M();
847  if (!fTrk) continue;
848  fastjet::PseudoJet PseudoTracks(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
849  TrackMom += TMath::Sqrt(TMath::Power(fTrk->Px(),2)+TMath::Power(fTrk->Py(),2)+TMath::Power(fTrk->Pz(),2));
850  TrackEnergy += fTrk->E();
851  PseudoTracks.set_user_index(fJet->TrackAt(i)+100);
852  PseudJetInvMass += PseudoTracks.m();
853  fInputVectors.push_back(PseudoTracks);
854  EmcalJetTrackEta[i]=fTrk->Eta();
855  EmcalJetTrackPhi[i]=fTrk->Phi();
856  EmcalJetTrackPt[i]=fTrk->Pt();
857  EmcalJetNTracks++;
858  }
859  fastjet::JetDefinition *fJetDef;
860  fastjet::ClusterSequence *fClustSeqSA;
861 
862 
863  fJetDef = new fastjet::JetDefinition(fastjet::antikt_algorithm, fJetRadius*2, static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
864 
865  try {
866  fClustSeqSA = new fastjet::ClusterSequence(fInputVectors, *fJetDef);
867  } catch (fastjet::Error) {
868  AliError(" [w] FJ Exception caught.");
869  //return -1;
870  }
871 
872  std::vector<fastjet::PseudoJet> fOutputJets;
873  fOutputJets.clear();
874  fOutputJets=fClustSeqSA->inclusive_jets(0);
875 
876  std::vector<fastjet::PseudoJet> jet_constituents = fOutputJets[0].constituents();
877  Float_t NSubjettinessResult[3], NSubBeta = 1, R0=0.4;
878  std::vector<fastjet::PseudoJet> Subjet_Axes;
879  fastjet::PseudoJet SubJet1_Axis,SubJet2_Axis;
880  Double_t DelR=-5;
881 
882 
883  for(Int_t j=1; j<3; j++){
884  if(jet_constituents.size() < j){
885  if(!fTruthJet){
886  if(j==1) fJetInfoVar[4]=-5;
887  else if(j==2) fJetInfoVar[6]=-5;
888  }
889  else {
890  if(j==1) fJetInfoVar[5]=-5;
891  else if(j==2) fJetInfoVar[7]=-5;
892  }
893  continue;
894  }
895  else{
896  fastjet::contrib::Nsubjettiness nSub(j,fastjet::contrib::KT_Axes(),fastjet::contrib::NormalizedMeasure(NSubBeta,R0));
897  NSubjettinessResult[j] = nSub.result(fOutputJets[0]);
898  if(j==2){ //find subjet axis to calculate delR
899  Subjet_Axes = nSub.currentAxes();
900  SubJet1_Axis = Subjet_Axes[0];
901  SubJet2_Axis = Subjet_Axes[1];
902 
903  Double_t SubJet1_Eta=SubJet1_Axis.pseudorapidity();
904  Double_t SubJet2_Eta=SubJet2_Axis.pseudorapidity();
905  Double_t SubJet1_Phi=SubJet1_Axis.phi();
906  if(SubJet1_Phi < -1*TMath::Pi()) SubJet1_Phi += (2*TMath::Pi());
907  else if (SubJet1_Phi > TMath::Pi()) SubJet1_Phi -= (2*TMath::Pi());
908  Double_t SubJet2_Phi=SubJet2_Axis.phi();
909  if(SubJet1_Phi < -1*TMath::Pi()) SubJet1_Phi += (2*TMath::Pi());
910  else if (SubJet1_Phi > TMath::Pi()) SubJet1_Phi -= (2*TMath::Pi());
911 
912  Double_t DeltaPhiSubJets,DeltaEtaSubJets;
913  DeltaPhiSubJets=SubJet1_Phi-SubJet2_Phi;
914  DeltaEtaSubJets=SubJet1_Eta-SubJet2_Eta;
915  if(DeltaPhiSubJets < -1*TMath::Pi()) DeltaPhiSubJets += (2*TMath::Pi());
916  else if (DeltaPhiSubJets > TMath::Pi()) DeltaPhiSubJets -= (2*TMath::Pi());
917 
918  DelR = TMath::Power(TMath::Power(DeltaPhiSubJets,2)+TMath::Power(DeltaEtaSubJets,2),0.5);
919  }
920  }
921  }
922  if(!fTruthJet){
923  fJetInfoVar[4]=NSubjettinessResult[1];
924  fJetInfoVar[6]=NSubjettinessResult[2];
925  // fJetInfoVar[12]=DelR;
926  }
927  else {
928  fJetInfoVar[5]=NSubjettinessResult[1];
929  fJetInfoVar[7]=NSubjettinessResult[2];
930  // fJetInfoVar[13]=DelR;
931  }
932 
933 
934 
935  fastjet::contrib::SoftDrop softdrop(beta, zcut);
936  //fastjet::contrib::SoftDrop softdrop_antikt(beta,zcut);
937  softdrop.set_verbose_structure(kTRUE);
938  //fastjet::JetDefinition jet_def_akt(fastjet::antikt_algorithm, 0.4);
939  // fastjet::contrib::Recluster *antiKT_Recluster(jet_def_akt);
940  fastjet::contrib::Recluster *recluster;
941  if(fReclusterAlgo == 2) recluster = new fastjet::contrib::Recluster(fastjet::kt_algorithm,1,true);
942  if(fReclusterAlgo == 1) recluster = new fastjet::contrib::Recluster(fastjet::antikt_algorithm,1,true);
943  if(fReclusterAlgo == 0) recluster = new fastjet::contrib::Recluster(fastjet::cambridge_algorithm,1,true);
944  softdrop.set_reclustering(true,recluster);
945  fastjet::PseudoJet finaljet = softdrop(fOutputJets[0]);
946  // fastjet::PseudoJet finaljet_antikt = softdrop_antikt(fOutputJets[0]);
947  //cout<< finaljet.structure_of<fastjet::contrib::SoftDrop>().symmetry()<<endl;
948  //cout<< finaljet_antikt.structure_of<fastjet::contrib::SoftDrop>().symmetry()<<endl;
949 
950 
951  AliEmcalJet* jet = new AliEmcalJet(finaljet.perp(), finaljet.eta(), finaljet.phi(), finaljet.m());
952  std::vector<fastjet::PseudoJet> fSDTracks=finaljet.constituents();
953  Double_t FastjetTrackDelR,EmcalTrackDelR;
954  for(Int_t i=0;i<fJet->GetNumberOfConstituents();i++){
955  if(i<=finaljet.constituents().size()){
956  FastjetTrackDelR = TMath::Sqrt(TMath::Power(fSDTracks[i].eta()-JetEta,2)+TMath::Power(fSDTracks[i].phi()-JetPhi,2));
957  FJTrackEta[i]=fSDTracks[i].eta();
958  FJTrackPhi[i]=fSDTracks[i].phi();
959  FJTrackPt[i]=fSDTracks[i].perp();
960  FJNTracks++;
961  }
962  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
963  EmcalTrackDelR = TMath::Sqrt(TMath::Power(fTrk->Eta()-JetEta,2)+TMath::Power(fTrk->Phi()-JetPhi,2));
964  }
965  Int_t NDroppedTracks = fJet->GetNumberOfTracks()-finaljet.constituents().size();
966  Int_t nConstituents(fClustSeqSA->constituents(finaljet).size());
967  jet->SetNumberOfTracks(nConstituents);
968  Double_t SymParam, Mu, DeltaR;
969  SymParam=(finaljet.structure_of<fastjet::contrib::SoftDrop>().symmetry());
970  Mu=(finaljet.structure_of<fastjet::contrib::SoftDrop>().mu());
971  DeltaR=(finaljet.structure_of<fastjet::contrib::SoftDrop>().delta_R());
972  //fhGroomedPtvJetPt->Fill(finaljet.perp(),fJet->Pt());
973  //fhDroppedBranches->Fill(finaljet.structure_of<fastjet::contrib::SoftDrop>().dropped_count());
974  if(!fTruthJet) fJetInfoVar[2]=SymParam;
975  else fJetInfoVar[3]=SymParam;
976  if(!fTruthJet) fJetInfoVar[12] = DeltaR;
977  else fJetInfoVar[13] = DeltaR;
978  if(!fTruthJet) fJetInfoVar[14]=finaljet.structure_of<fastjet::contrib::SoftDrop>().dropped_count();
979  else fJetInfoVar[15]=finaljet.structure_of<fastjet::contrib::SoftDrop>().dropped_count();
980  if(!(fJetShapeSub==kConstSub)){
981  if(!fTruthJet) fJetInfoVar[16]=(finaljet.perp()-(GetRhoVal(0)*fJet->Area()));
982  else fJetInfoVar[17]=finaljet.perp();
983  }
984  else{
985  if(!fTruthJet) fJetInfoVar[16]=(finaljet.perp());
986  else fJetInfoVar[17]=(finaljet.perp());
987  }
988 
989  return;
990 
991 
992 }
993 
994 //--------------------------------------------------------------------------
995 Int_t AliAnalysisTaskRecoilJetYield::SelectTriggerHadron(Float_t PtMin, Float_t PtMax){
996 
997  AliTrackContainer *PartCont = NULL;
998  AliParticleContainer *PartContMC = NULL;
999 
1000 
1001  if (fJetShapeSub==kConstSub){
1002  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kGenOnTheFly) PartContMC = GetParticleContainer(1);
1003  else PartCont = GetTrackContainer(1);
1004  }
1005  else{
1006  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kGenOnTheFly) PartContMC = GetParticleContainer(0);
1007  else PartCont = GetTrackContainer(0);
1008  }
1009 
1010  TClonesArray *TracksArray = NULL;
1011  TClonesArray *TracksArrayMC = NULL;
1012  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kGenOnTheFly) TracksArrayMC = PartContMC->GetArray();
1013  else TracksArray = PartCont->GetArray();
1014 
1015  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kGenOnTheFly){
1016  if(!PartContMC || !TracksArrayMC) return -99999;
1017  }
1018  else {
1019  if(!PartCont || !TracksArray) return -99999;
1020  }
1021 
1022  AliAODTrack *Track = 0x0;
1023  Int_t Trigger_Index[100];
1024  for (Int_t i=0; i<100; i++) Trigger_Index[i] = 0;
1025  Int_t Trigger_Counter = 0;
1026  Int_t NTracks=0;
1027  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kGenOnTheFly) NTracks = TracksArrayMC->GetEntriesFast();
1028  else NTracks = TracksArray->GetEntriesFast();
1029  for(Int_t i=0; i < NTracks; i++){
1030  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kGenOnTheFly){
1031  if((Track = static_cast<AliAODTrack*>(PartContMC->GetAcceptParticle(i)))){
1032  if (!Track) continue;
1033  if(TMath::Abs(Track->Eta())>0.9) continue;
1034  if (Track->Pt()<0.15) continue;
1035  if ((Track->Pt() >= PtMin) && (Track->Pt()< PtMax)) {
1036  Trigger_Index[Trigger_Counter] = i;
1037  Trigger_Counter++;
1038  }
1039  }
1040  }
1041  else{
1042  if((Track = static_cast<AliAODTrack*>(PartCont->GetAcceptTrack(i)))){
1043  if (!Track) continue;
1044  if(TMath::Abs(Track->Eta())>0.9) continue;
1045  if (Track->Pt()<0.15) continue;
1046  if ((Track->Pt() >= PtMin) && (Track->Pt()< PtMax)) {
1047  Trigger_Index[Trigger_Counter] = i;
1048  Trigger_Counter++;
1049  }
1050  }
1051  }
1052  }
1053  if (Trigger_Counter == 0) return -99999;
1054  Int_t RandomNumber = 0, Index = 0 ;
1055  TRandom3* Random = new TRandom3(0);
1056  RandomNumber = Random->Integer(Trigger_Counter);
1057  Index = Trigger_Index[RandomNumber];
1058  return Index;
1059 }
1060 
1061 
1062 //________________________________________________________________________
1063 Double_t AliAnalysisTaskRecoilJetYield::GetFractionSharedPt_SubMatching(const AliEmcalJet *jet1, AliParticleContainer *cont2) const
1064 {
1065  AliEmcalJet *jet2 = jet1->ClosestJet();
1066  if (!jet2) return -1;
1067 
1068  Double_t jetPt2 = jet2->Pt();
1069  if (jetPt2 <= 0) return -1;
1070 
1071  Int_t bgeom = kTRUE;
1072  if (!cont2) bgeom = kFALSE;
1073  Double_t sumPt = 0.;
1074  AliVParticle *vpf = 0x0;
1075  Int_t iFound = 0;
1076  for (Int_t icc = 0; icc < jet2->GetNumberOfTracks(); icc++) {
1077  Int_t idx = (Int_t)jet2->TrackAt(icc);
1078  //get particle
1079  AliVParticle *p2 = 0x0;
1080  if (bgeom) p2 = static_cast<AliVParticle*>(jet2->TrackAt(icc, cont2->GetArray()));
1081  iFound = 0;
1082  for (Int_t icf = 0; icf < jet1->GetNumberOfTracks(); icf++) {
1083  if (!bgeom && idx == jet1->TrackAt(icf) && iFound == 0 ) {
1084  iFound = 1;
1085  vpf = jet1->Track(icf);
1086  if (vpf) sumPt += vpf->Pt();
1087  continue;
1088  }
1089  if (bgeom){
1090  vpf = jet1->Track(icf);
1091  if (!vpf) continue;
1092  if (!SameParticle(vpf, p2, 1.e-4)) continue; //not the same particle
1093  sumPt += vpf->Pt();
1094  }
1095  }
1096  }
1097 
1098  Double_t fraction = sumPt / jetPt2;
1099 
1100  return fraction;
1101 }
1102 
1103 //________________________________________________________________________
1104 Bool_t AliAnalysisTaskRecoilJetYield::SameParticle(const AliVParticle* part1, const AliVParticle* part2, Double_t dist)
1105 {
1106  if(!part1) return kFALSE;
1107  if(!part2) return kFALSE;
1108  Double_t dPhi = TMath::Abs(part1->Phi() - part2->Phi());
1109  dPhi = TVector2::Phi_mpi_pi(dPhi);
1110  if (dPhi > dist) return kFALSE;
1111  return kTRUE;
1112 }
1113 
1114 //________________________________________________________________________
1115 Bool_t AliAnalysisTaskRecoilJetYield::RetrieveEventObjects() {
1116  //
1117  // retrieve event objects
1118  //
1119  if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
1120  return kFALSE;
1121 
1122  return kTRUE;
1123 }
1124 
1125 
1126 //_______________________________________________________________________
1127 void AliAnalysisTaskRecoilJetYield::Terminate(Option_t *)
1128 {
1129  // Called once at the end of the analysis.
1130 
1131 
1132 }
AliAnalysisTaskRecoilJetYield::RetrieveEventObjects
Bool_t RetrieveEventObjects()
Retrieve common objects from event.
Definition: AliAnalysisTaskRecoilJetYield.cxx:1115
AliEmcalJet::Area
Double_t Area() const
Definition: AliEmcalJet.h:123
AliAnalysisTaskRecoilJetYield.h
AliPicoTrack.h
Double_t
double Double_t
Definition: External.C:58
AliEmcalJetFinder.h
AliAnalysisTaskRecoilJetYield::PTD
Double_t PTD(AliEmcalJet *Jet, Int_t JetContNb)
Definition: AliAnalysisTaskRecoilJetYield.cxx:818
AliAnalysisTaskEmcal::DeltaPhi
static Double_t DeltaPhi(Double_t phia, Double_t phib, Double_t rMin=-TMath::Pi()/2, Double_t rMax=3 *TMath::Pi()/2)
Definition: AliAnalysisTaskEmcal.h:932
TH2F
Definition: External.C:236
AliEmcalJet::ClosestJet
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:222
AliAnalysisTaskEmcalJet::GetJetContainer
AliJetContainer * GetJetContainer(Int_t i=0) const
Definition: AliAnalysisTaskEmcalJet.cxx:361
AliAnalysisTaskRecoilJetYield::FillHistograms
Bool_t FillHistograms()
Function filling histograms.
Definition: AliAnalysisTaskRecoilJetYield.cxx:369
AliEmcalJet::Eta
Double_t Eta() const
Definition: AliEmcalJet.h:114
Phi_Upper
Double_t Phi_Upper
Definition: AliAnalysisTaskRecoilJetYield.cxx:55
AliEmcalJet::Phi
Double_t Phi() const
Definition: AliEmcalJet.h:110
ClassImp
ClassImp(AliAnalysisTaskRecoilJetYield) AliAnalysisTaskRecoilJetYield
Definition: AliAnalysisTaskRecoilJetYield.cxx:66
AliTrackContainer
Container with name, TClonesArray and cuts for particles.
Definition: AliTrackContainer.h:29
AliAnalysisTaskEmcal::fEPV0
Double_t fEPV0
!event plane V0
Definition: AliAnalysisTaskEmcal.h:876
AliEmcalJet::GetLabel
Int_t GetLabel() const
Definition: AliEmcalJet.h:117
AliEmcalParticle.h
AliEmcalJet::Track
AliVParticle * Track(Int_t idx) const
Definition: AliEmcalJet.cxx:683
AliJetContainer::AcceptJet
virtual Bool_t AcceptJet(Int_t i, UInt_t &rejectionReason) const
Definition: AliJetContainer.cxx:483
Eta_Lower
Double_t Eta_Lower
Definition: AliAnalysisTaskRecoilJetYield.cxx:53
Percentile_1
Float_t Percentile_1
Definition: AliAnalysisTaskRecoilJetYield.cxx:58
Eta_Upper
Double_t Eta_Upper
Definition: AliAnalysisTaskRecoilJetYield.cxx:52
AliAnalysisTaskRecoilJetYield::SoftDrop
void SoftDrop(AliEmcalJet *fJet, AliJetContainer *fJetCont, double zcut, double beta, Bool_t fTruthJet)
Definition: AliAnalysisTaskRecoilJetYield.cxx:837
Phi_Bins_1
Int_t Phi_Bins_1
Definition: AliAnalysisTaskRecoilJetYield.cxx:57
Phi_Lower
Double_t Phi_Lower
Definition: AliAnalysisTaskRecoilJetYield.cxx:56
AliEmcalJet::GetNumberOfConstituents
UShort_t GetNumberOfConstituents() const
Definition: AliEmcalJet.h:133
AliParticleContainer
Container for particles within the EMCAL framework.
Definition: AliParticleContainer.h:31
AliEmcalJet::TrackAt
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:153
AliEmcalJet::GetNumberOfTracks
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:132
kData
TString kData
Declare data MC or deltaAOD.
Definition: AddTaskCaloTrackCorr.C:22
AliAnalysisTaskEmcal::GetParticleContainer
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
Definition: AliAnalysisTaskEmcal.cxx:1581
AliJetContainer::GetParticleContainer
AliParticleContainer * GetParticleContainer() const
Definition: AliJetContainer.h:180
Int_t
int Int_t
Definition: External.C:63
UInt_t
unsigned int UInt_t
Definition: External.C:33
Float_t
float Float_t
Definition: External.C:68
AliAnalysisTaskRecoilJetYield::SelectTriggerHadron
Int_t SelectTriggerHadron(Float_t PtMin, Float_t PtMax)
Definition: AliAnalysisTaskRecoilJetYield.cxx:995
AliAnalysisTaskRecoilJetYield::GetFractionSharedPt_SubMatching
Double_t GetFractionSharedPt_SubMatching(const AliEmcalJet *jet, AliParticleContainer *cont2=0x0) const
Definition: AliAnalysisTaskRecoilJetYield.cxx:1063
AliParticleContainer.h
AliJetContainer::GetNJets
Int_t GetNJets() const
Definition: AliJetContainer.h:153
AliJetContainer.h
AliParticleContainer::GetAcceptParticle
virtual AliVParticle * GetAcceptParticle(Int_t i=-1) const
Definition: AliParticleContainer.cxx:113
AliAnalysisTaskEmcal::fCent
Double_t fCent
!event centrality
Definition: AliAnalysisTaskEmcal.h:874
nBranch
const Int_t nBranch
Definition: AliAnalysisTaskRecoilJetYield.h:22
AliAnalysisTaskRecoilJetYield::Angularity
Double_t Angularity(AliEmcalJet *Jet, Int_t JetContNb)
Definition: AliAnalysisTaskRecoilJetYield.cxx:794
AliRhoParameter.h
AliJetContainer::GetNextAcceptJet
AliEmcalJet * GetNextAcceptJet()
Definition: AliJetContainer.cxx:299
Data
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)
Definition: charmCutsOptimization.C:371
AliAnalysisTaskRecoilJetYield::RelativePhi
Double_t RelativePhi(Double_t Phi1, Double_t Phi2)
Definition: AliAnalysisTaskRecoilJetYield.cxx:778
AliAnalysisTaskRecoilJetYield::SameParticle
static Bool_t SameParticle(const AliVParticle *part1, const AliVParticle *part2, Double_t dist=1.e-4)
Definition: AliAnalysisTaskRecoilJetYield.cxx:1104
AliEmcalJet::Pt
Double_t Pt() const
Definition: AliEmcalJet.h:102
AliAnalysisTaskEmcalJet::GetRhoVal
Double_t GetRhoVal(Int_t i=0) const
Definition: AliAnalysisTaskEmcalJet.cxx:610
AliAnalysisTaskEmcal::fOutput
AliEmcalList * fOutput
!output list
Definition: AliAnalysisTaskEmcal.h:896
AliAnalysisTaskEmcal::GetTrackContainer
AliTrackContainer * GetTrackContainer(Int_t i=0) const
Definition: AliAnalysisTaskEmcal.h:223
NS_AliAnalysisTaskEmcalJetCDF::DeltaR
Double_t DeltaR(const AliVParticle *part1, const AliVParticle *part2)
Definition: AliAnalysisTaskEmcalJetCDF.h:144
AliAnalysisTaskEmcal::SetMakeGeneralHistograms
void SetMakeGeneralHistograms(Bool_t g)
Definition: AliAnalysisTaskEmcal.h:256
AliTrackContainer::GetAcceptTrack
virtual AliVTrack * GetAcceptTrack(Int_t i=-1) const
Definition: AliTrackContainer.cxx:208
AliAnalysisTaskEmcalJet
Base task in the EMCAL jet framework.
Definition: AliAnalysisTaskEmcalJet.h:30
AliEmcalJet
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
AliEmcalJet::SetNumberOfTracks
void SetNumberOfTracks(Int_t n)
Definition: AliEmcalJet.h:195
Eta_Bins_1
Int_t Eta_Bins_1
Definition: AliAnalysisTaskRecoilJetYield.cxx:54
UShort_t
unsigned short UShort_t
Definition: External.C:28
Option_t
const char Option_t
Definition: External.C:48
AliAnalysisTaskEmcal::UserCreateOutputObjects
void UserCreateOutputObjects()
Main initialization function on the worker.
Definition: AliAnalysisTaskEmcal.cxx:331
Bool_t
bool Bool_t
Definition: External.C:53
AliAnalysisTaskRecoilJetYield::Run
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.
Definition: AliAnalysisTaskRecoilJetYield.cxx:362
AliAnalysisTaskRecoilJetYield::RelativePhiEventPlane
Double_t RelativePhiEventPlane(Double_t EventPlane, Double_t Phi)
Definition: AliAnalysisTaskRecoilJetYield.cxx:768
AliEmcalJet::M
Double_t M() const
Definition: AliEmcalJet.h:113
AliEmcalJet.h
AliJetContainer
Container for jet within the EMCAL jet framework.
Definition: AliJetContainer.h:37
AliJetContainer::GetJet
AliEmcalJet * GetJet(Int_t i) const
Definition: AliJetContainer.cxx:270