AliPhysics  a17849b (a17849b)
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 = 20;
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] = "Mass";
255  fJetInfoVarNames[5] = "Mass_Truth";
256  fJetInfoVarNames[6] = "SLSubJetMass";
257  fJetInfoVarNames[7] = "SLSubJetMass_Truth";
258  fJetInfoVarNames[8] = "LeadingParton";
259  fJetInfoVarNames[9] = "LeadingParton_Truth";
260  fJetInfoVarNames[10] = "SLSubJetPt";
261  fJetInfoVarNames[11] = "SLSubJetPt_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  fJetInfoVarNames[18] = "Groomed_Mass";
269  fJetInfoVarNames[19] = "Groomed_Mass_Truth";
270 
271 
272 
273  for(Int_t ivar=0; ivar < nVarMin; ivar++){
274  cout<<"looping over variables"<<endl;
275  fTreeJetInfo->Branch(fJetInfoVarNames[ivar].Data(), &fJetInfoVar[ivar], Form("%s/D", fJetInfoVarNames[ivar].Data()));
276  }
277  }
278  if(fJetSelection == kRecoil){
279  fhPtTriggerHadron= new TH1F("fhPtTriggerHadron", "fhPtTriggerHadron",1500,-0.5,149.5);
280  fOutput->Add(fhPtTriggerHadron);
281  fh2PtTriggerHadronJet= new TH2F("fh2PtTriggerHadronJet", "fh2PtTriggerHadronJet",1500,-0.5,149.5,1500,-0.5,149.5);
282  fOutput->Add(fh2PtTriggerHadronJet);
283  fhPhiTriggerHadronJet= new TH1F("fhPhiTriggerHadronJet", "fhPhiTriggerHadronJet",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
284  fOutput->Add(fhPhiTriggerHadronJet);
285  fhPhiTriggerHadronEventPlane= new TH1F("fhPhiTriggerHadronEventPlane", "fhPhiTriggerHadronEventPlane",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
286  fOutput->Add(fhPhiTriggerHadronEventPlane);
287  fhPhiTriggerHadronEventPlaneTPC= new TH1F("fhPhiTriggerHadronEventPlaneTPC", "fhPhiTriggerHadronEventPlaneTPC",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
288  fOutput->Add(fhPhiTriggerHadronEventPlaneTPC);
289 
290 
291  }
292 
293  if (fJetShapeType==AliAnalysisTaskRecoilJetYield::kData || fJetShapeType==AliAnalysisTaskRecoilJetYield::kDetEmbPart || fJetShapeType==AliAnalysisTaskRecoilJetYield::kGenOnTheFly){
294 
295  fhJetPt= new TH1F("fhJetPt", "Jet Pt",150,-0.5,149.5 );
296  fOutput->Add(fhJetPt);
297  fhJetPhi= new TH1F("fhJetPhi", "Jet Phi",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
298  fOutput->Add(fhJetPhi);
299  fhJetEta= new TH1F("fhJetEta", "Jet Eta", Eta_Bins_1, Eta_Lower, Eta_Upper);
300  fOutput->Add(fhJetEta);
301  fhJetMass= new TH1F("fhJetMass", "Jet Mass", 4000,-0.5, 39.5);
302  fOutput->Add(fhJetMass);
303  fhJetRadius= new TH1F("fhJetRadius", "Jet Radius", 100, -0.05,0.995);
304  fOutput->Add(fhJetRadius);
305  fhNumberOfJetTracks= new TH1F("fhNumberOfJetTracks", "Number of Tracks within a Jet", 300, -0.5,299.5);
306  fOutput->Add(fhNumberOfJetTracks);
307  fhJetCounter= new TH1F("fhJetCounter", "Jet Counter", 150, -0.5, 149.5);
308  fOutput->Add(fhJetCounter);
309  fhJetArea= new TH1F("fhJetArea", "Jet Area", 400,-0.5, 1.95);
310  fOutput->Add(fhJetArea);
311  fhTrackPt= new TH1F("fhTrackPt", "Track Pt",600,-0.5,59.5);
312  fOutput->Add(fhTrackPt);
313  fhGroomedPtvJetPt= new TH2F("fhGroomedPtvJetPt","Groomed Jet p_{T} v Original Jet p_{T}",150,0,150,150,0,150);
314  fOutput->Add(fhGroomedPtvJetPt);
315  fhDroppedBranches= new TH1F("fhDroppedBranches","Number of Softdropped branches",50,0,50);
316  fOutput->Add(fhDroppedBranches);
317  fhDetJetPt_Incl= new TH1F("fhDetJetPt_Incl", "Jet Pt",200,-0.5,199.5 );
318  fOutput->Add(fhDetJetPt_Incl);
319  fhDetJetPt_Matched= new TH1F("fhDetJetPt_Matched", "Jet Pt",200,-0.5,199.5 );
320  fOutput->Add(fhDetJetPt_Matched);
321 
322 
323  }
324  if(fJetShapeType==AliAnalysisTaskRecoilJetYield::kTrueDet){
325  fhJetPt_Det= new TH1F("fhJetPt_Det", "Jet Pt Detector Level",1500,-0.5,149.5 );
326  fOutput->Add(fhJetPt_Det);
327  fhJetPt_True= new TH1F("fhJetPt_True", "Jet Pt Particle Level",1500,-0.5,149.5 );
328  fOutput->Add(fhJetPt_True);
329  fhJetPhi_Det= new TH1F("fhJetPhi_Det", "Jet Phi Detector Level",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
330  fOutput->Add(fhJetPhi_Det);
331  fhJetPhi_True= new TH1F("fhJetPhi_True", "Jet Phi Particle Level",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
332  fOutput->Add(fhJetPhi_True);
333  fhJetEta_Det= new TH1F("fhJetEta_Det", "Jet Eta Detector Level", Eta_Bins_1, Eta_Lower, Eta_Upper);
334  fOutput->Add(fhJetEta_Det);
335  fhJetEta_True= new TH1F("fhJetEta_True", "Jet Eta Particle Level", Eta_Bins_1, Eta_Lower, Eta_Upper);
336  fOutput->Add(fhJetEta_True);
337  fhJetMass_Det= new TH1F("fhJetMass_Det", "Jet Mass Detector Level", 4000,-0.5, 39.5);
338  fOutput->Add(fhJetMass_Det);
339  fhJetMass_True= new TH1F("fhJetMass_True", "Jet Mass Particle Level", 4000,-0.5, 39.5);
340  fOutput->Add(fhJetMass_True);
341  fhJetRadius_Det= new TH1F("fhJetRadius_Det", "Jet Radius Detector Level", 100, -0.05,0.995);
342  fOutput->Add(fhJetRadius_Det);
343  fhJetRadius_True= new TH1F("fhJetRadius_True", "Jet Radius Particle Level", 100, -0.05,0.995);
344  fOutput->Add(fhJetRadius_True);
345  fhNumberOfJetTracks_Det= new TH1F("fhNumberOfJetTracks_Det", "Number of Tracks within a Jet Detector Level", 300, -0.5,299.5);
346  fOutput->Add(fhNumberOfJetTracks_Det);
347  fhNumberOfJetTracks_True= new TH1F("fhNumberOfJetTracks_True", "Number of Tracks within a Jet Particle Level", 300, -0.5,299.5);
348  fOutput->Add(fhNumberOfJetTracks_True);
349 
350  fhJetCounter_Det= new TH1F("fhJetCounter_Det", "Jet Counter Detector Level", 150, -0.5, 149.5);
351  fOutput->Add(fhJetCounter_Det);
352  fhJetCounter_True= new TH1F("fhJetCounter_True", "Jet Counter Particle Level", 150, -0.5, 149.5);
353  fOutput->Add(fhJetCounter_True);
354  fh2PtRatio= new TH2F("fhPtRatio", "MC pt for detector level vs particle level jets",1500,-0.5,149.5,1500,-0.5,149.5);
355  fOutput->Add(fh2PtRatio);
356  }
357  PostData(1,fOutput);
358  PostData(2,fTreeJetInfo);
359  //cout<<"End of UserCreateOutputObjects"<<endl;
360  // delete [] fShapesVarNames;
361 }
362 
363 //________________________________________________________________________
364 Bool_t AliAnalysisTaskRecoilJetYield::Run()
365 {
366  // Run analysis code here, if needed. It will be executed before FillHistograms().
367  return kTRUE;
368 }
369 
370 //________________________________________________________________________
371 Bool_t AliAnalysisTaskRecoilJetYield::FillHistograms()
372 {
373  //AliMultSelection *MultSelection = 0x0;
374  AliAODEvent *fEvent = dynamic_cast<AliAODEvent*>(InputEvent());
375  if (!fEvent) {
376  Error("UserExec","AOD not available");
377  return 0;
378  }
379  /*MultSelection = (AliMultSelection * ) fEvent->FindListObject("MultSelection");
380  if(!MultSelection) {
381  AliWarning("AliMultSelection object not found!");
382  }
383  else{
384  Percentile_1 = MultSelection->GetMultiplicityPercentile("V0M");
385  }
386  */
387  if (fCentSelectOn){
388  if ((fCent>fCentMax) || (fCent<fCentMin)) return 0;
389  }
390 
392  AliAODTrack *TriggerHadron = 0x0;
393  if (fJetSelection == kRecoil) {
394  //you have to set a flag and the limits of the pT interval for your trigger
395  Int_t TriggerHadronLabel = SelectTriggerHadron(fPtMinTriggerHadron, fPtMaxTriggerHadron);
396  if (TriggerHadronLabel==-99999) return 0; //Trigger Hadron Not Found
397  AliTrackContainer *PartCont =NULL;
398  AliParticleContainer *PartContMC = NULL;
399  if (fJetShapeSub==kConstSub){
400  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kGenOnTheFly) PartContMC = GetParticleContainer(1);
401  else PartCont = GetTrackContainer(1);
402  }
403  else {
404  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kGenOnTheFly) PartContMC = GetParticleContainer(0);
405  else PartCont = GetTrackContainer(0);
406  }
407  TClonesArray *TrackArray = NULL;
408  TClonesArray *TrackArrayMC = NULL;
409  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kGenOnTheFly) TrackArrayMC = PartContMC->GetArray();
410  else TrackArray = PartCont->GetArray();
411  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kGenOnTheFly) TriggerHadron = static_cast<AliAODTrack*>(TrackArrayMC->At(TriggerHadronLabel));
412  else TriggerHadron = static_cast<AliAODTrack*>(TrackArray->At(TriggerHadronLabel));
413  if (!TriggerHadron) return 0;//No trigger hadron with label found
414  if(fSemigoodCorrect){
415  Double_t HoleDistance=RelativePhi(TriggerHadron->Phi(),fHolePos);
416  if(TMath::Abs(HoleDistance)+fHoleWidth+fJetRadius>TMath::Pi()-fRecoilAngularWindow) return 0;
417  }
418  fhPtTriggerHadron->Fill(TriggerHadron->Pt()); //Needed for per trigger Normalisation
419  if (fJetShapeType != AliAnalysisTaskRecoilJetYield::kGenOnTheFly) fhPhiTriggerHadronEventPlane->Fill(TMath::Abs(RelativePhiEventPlane(fEPV0,TriggerHadron->Phi()))); //fEPV0 is the event plane from AliAnalysisTaskEmcal
420  if (fJetShapeType != AliAnalysisTaskRecoilJetYield::kGenOnTheFly) fhPhiTriggerHadronEventPlaneTPC->Fill(TMath::Abs(RelativePhiEventPlane(((AliVAODHeader*)(InputEvent()->GetHeader()))->GetEventplane(),TriggerHadron->Phi()))); //TPC event plane
421  }
422 
423 
424 
425 
426 
428  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kData){
429  // cout<<"entering kData"<<endl;
430 
431  AliEmcalJet *Jet1 = NULL; //Original Jet in the event
432  AliJetContainer *JetCont= GetJetContainer(0); //Jet Container for event
433  Int_t JetCounter=0; //Counts number of jets in event
434  Double_t JetPhi=0;
435  Bool_t EventCounter=kFALSE;
436  Double_t JetPt_ForThreshold=0;
437  AliEmcalJet *Jet2= NULL;
438  if(JetCont) {
439  JetCont->ResetCurrentID();
440  while((Jet1=JetCont->GetNextAcceptJet())) {
441  if(!Jet1) continue;
442  if (fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) JetPt_ForThreshold = Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
443  else JetPt_ForThreshold = Jet1->Pt();
444  if(JetPt_ForThreshold<fPtThreshold) {
445  continue;
446  }
447  else {
448  Float_t RecoilDeltaPhi = 0.;
449  if (fJetSelection == kRecoil){
450  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), Jet1->Phi());
451  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
452  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), Jet1->Pt());
453  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), Jet1->Phi()));
454  }
455  if (!EventCounter){
456  EventCounter=kTRUE;
457  }
458  JetCounter++;
459  fhJetPt->Fill(Jet1->Pt());
460  fhJetArea->Fill(Jet1->Area());
461  JetPhi=Jet1->Phi();
462  if(JetPhi < -1*TMath::Pi()) JetPhi += (2*TMath::Pi());
463  else if (JetPhi > TMath::Pi()) JetPhi -= (2*TMath::Pi());
464  fhJetPhi->Fill(JetPhi);
465  fhJetEta->Fill(Jet1->Eta());
466  fhJetMass->Fill(Jet1->M());
467  fJetInfoVar[4]=Jet1->M();
468  fhJetRadius->Fill(TMath::Sqrt((Jet1->Area()/TMath::Pi()))); //Radius of Jets per event
469  fhNumberOfJetTracks->Fill(Jet1->GetNumberOfTracks());
470  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) fJetInfoVar[0]= Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
471  else fJetInfoVar[0]=Jet1->Pt();
472  fJetInfoVar[1]=0;
473  if(fDoSoftDrop) {
474  SoftDrop(Jet1,JetCont,fZCut,fBeta_SD,kFALSE);
475  //SoftDrop(Jet1,JetCont,fZCut,fBeta_SD,kTRUE);
476 
477  fJetInfoVar[3]=0;
478  fJetInfoVar[5]=0;
479  fJetInfoVar[7]=0;
480  fJetInfoVar[9]=0;
481  fJetInfoVar[11]=0;
482  fJetInfoVar[13]=0;
483  fJetInfoVar[15]=0;
484  fJetInfoVar[17]=0;
485  fJetInfoVar[19]=0;
486 
487  }
488  else{
489  fJetInfoVar[2]=0;
490  fJetInfoVar[3]=0;
491  fJetInfoVar[4]=0;
492  fJetInfoVar[5]=0;
493  fJetInfoVar[6]=0;
494  fJetInfoVar[7]=0;
495  fJetInfoVar[8]=0;
496  fJetInfoVar[9]=0;
497  fJetInfoVar[10]=0;
498  fJetInfoVar[11]=0;
499  fJetInfoVar[12]=0;
500  fJetInfoVar[13]=0;
501  fJetInfoVar[14]=0;
502  fJetInfoVar[15]=0;
503  fJetInfoVar[16]=0;
504  fJetInfoVar[17]=0;
505  fJetInfoVar[18]=0;
506  fJetInfoVar[19]=0;
507  }
508  fTreeJetInfo->Fill();
509  }
510  }
511  fhJetCounter->Fill(JetCounter); //Number of Jets in Each Event
512  }
513  }
514 
516  if(fJetShapeType == kDetEmbPart){
517  AliEmcalJet *JetHybridS = NULL; //Subtracted hybrid Jet
518  AliEmcalJet *JetHybridUS = NULL; //Unsubtracted Hybrid Jet //For matching SubtractedHybrid->DetPythia this jet container is also Subtracted Hybrid
519  AliEmcalJet *JetPythDet = NULL; //Detector Level Pythia Jet
520  AliEmcalJet *JetPythTrue = NULL; //Particle Level Pyhtia Jet
521  AliJetContainer *JetContHybridS= GetJetContainer(0); //Jet Container for Subtracted Hybrid Jets
522  AliJetContainer *JetContHybridUS= GetJetContainer(1); //Jet Container for Unsubtracted Hybrid Jets
523  AliJetContainer *JetContPythDet= GetJetContainer(2); //Jet Container for Detector Level Pyhtia Jets
524  AliJetContainer *JetContPythTrue= GetJetContainer(3); //Jet Container for Particle Level Pythia Jets
525 
526 
527 
528  Bool_t JetsMatched = kFALSE;
529  Double_t JetPtThreshold;
530  JetContHybridS->ResetCurrentID();
531  JetContHybridUS->ResetCurrentID();
532  JetContPythDet->ResetCurrentID();
533  JetContPythTrue->ResetCurrentID();
534 
535  while((JetPythDet = JetContPythDet->GetNextAcceptJet())){
536  fhDetJetPt_Incl->Fill(JetPythDet->Pt()); //Fill histogram with all Detector level jets
537  }
538 
539 
540  while((JetHybridS = JetContHybridS->GetNextAcceptJet())){
541  if (fJetShapeSub==kConstSub) JetPtThreshold=JetHybridS->Pt();
542  else JetPtThreshold=JetHybridS->Pt()-(GetRhoVal(0)*JetHybridS->Area());
543  if ( (!JetHybridS) || (JetPtThreshold<fPtThreshold)) continue;
544  Int_t JetNumber=-1;
545  for(Int_t i = 0; i<JetContHybridUS->GetNJets(); i++) {
546  JetHybridUS = JetContHybridUS->GetJet(i);
547  if (!JetHybridUS) continue;
548 
549  if(JetHybridUS->GetLabel()==JetHybridS->GetLabel()) {
550  JetNumber=i;
551  }
552  }
553  if(JetNumber==-1) continue;
554  JetHybridUS=JetContHybridUS->GetJet(JetNumber);
555  if (JetContHybridUS->AliJetContainer::GetFractionSharedPt(JetHybridUS)<fSharedFractionPtMin && !fSubMatching) {
556  continue;
557  }
558  else if(GetFractionSharedPt_SubMatching(JetHybridUS,JetContHybridUS->GetParticleContainer())<fSharedFractionPtMin && fSubMatching){
559  continue;
560  }
561  JetPythDet=JetHybridUS->ClosestJet();
562  if (!JetHybridUS) {
563  Printf("Unsubtracted embedded jet does not exist, returning");
564  continue;
565  }
566  if (!JetPythDet) continue;
567  UInt_t rejectionReason = 0;
568  if (!(JetContPythDet->AcceptJet(JetPythDet,rejectionReason))) continue;
569  fhDetJetPt_Matched->Fill(JetPythDet->Pt()); //Fill only matched detector level jets for tagging efficiency comparison
570  JetPythTrue=JetPythDet->ClosestJet();
571  if(!JetPythTrue) continue;
572  JetsMatched=kTRUE;
573 
574  if(fJetShapeSub==kConstSub) fJetInfoVar[0]=JetHybridS->Pt();
575  else fJetInfoVar[0]=JetHybridS->Pt()-(GetRhoVal(0)*JetHybridS->Area());
576  fJetInfoVar[1]=JetPythTrue->Pt();
577  if(fDoSoftDrop) {
578  SoftDrop(JetHybridS,JetContHybridS,fZCut,fBeta_SD,kFALSE);
579  SoftDrop(JetPythTrue,JetContPythTrue,fZCut,fBeta_SD,kTRUE);
580  }
581  else{
582  fJetInfoVar[2]=0;
583  fJetInfoVar[3]=0;
584  fJetInfoVar[4]=0;
585  fJetInfoVar[5]=0;
586  fJetInfoVar[6]=0;
587  fJetInfoVar[7]=0;
588  fJetInfoVar[8]=0;
589  fJetInfoVar[9]=0;
590  fJetInfoVar[10]=0;
591  fJetInfoVar[11]=0;
592  fJetInfoVar[12]=0;
593  fJetInfoVar[13]=0;
594  fJetInfoVar[14]=0;
595  fJetInfoVar[15]=0;
596  fJetInfoVar[16]=0;
597  fJetInfoVar[17]=0;
598  fJetInfoVar[18]=0;
599  fJetInfoVar[19]=0;
600  }
601  fJetInfoVar[4]=JetHybridS->M();
602  fJetInfoVar[5]=JetPythTrue->M();
603  fTreeJetInfo->Fill();
604  }
605 
606  }
607 
609  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kTrueDet){ //Truth->Detector response
610  AliEmcalJet *JetDet = NULL; //Detector Level Jet
611  AliEmcalJet *JetTrue = NULL; //Particle Level Jet
612  AliJetContainer *JetContDet= GetJetContainer(0); //Jet Container for Detector Level Pythia
613  AliJetContainer *JetContTrue= GetJetContainer(1); //Jet Container for Particle Level Pythia
614  Double_t JetPhiDet=0;
615  Double_t JetPhiTrue=0;
616  Bool_t JetsMatched=kFALSE;
617  Double_t Pythia_Event_Weight=1;
618  Bool_t EventCounter=kFALSE;
619  Int_t JetCounter_Det=0,JetCounter_True=0;
620  JetContDet->ResetCurrentID();
621  JetContTrue->ResetCurrentID();
622  if(JetContDet) {
623  while((JetDet=JetContDet->GetNextAcceptJet())) {
624  if( (!JetDet) || ((JetDet->Pt())<fPtThreshold)) {
625  // fhEventCounter_1->Fill(3); //events where the jet had a problem
626  continue;
627  }
628  else {
629  /* if(fSemigoodCorrect){
630  Double_t HoleDistance=RelativePhi(Jet1->Phi(),fHolePos);
631  if(TMath::Abs(HoleDistance)<fHoleWidth) continue;
632  }*/
633  Float_t RecoilDeltaPhi = 0.;
634  if (fJetSelection == kRecoil){
635  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), JetDet->Phi());
636  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
637  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), JetDet->Pt());
638  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), JetDet->Phi()));
639  }
640  JetCounter_Det++;
641  fhJetPt_Det->Fill(JetDet->Pt());
642  JetPhiDet=JetDet->Phi();
643  if(JetPhiDet < -1*TMath::Pi()) JetPhiDet += (2*TMath::Pi());
644  else if (JetPhiDet > TMath::Pi()) JetPhiDet -= (2*TMath::Pi());
645  fhJetPhi_Det->Fill(JetPhiDet);
646  fhJetEta_Det->Fill(JetDet->Eta());
647  fhJetMass_Det->Fill(JetDet->M());
648  fhJetRadius_Det->Fill(TMath::Sqrt((JetDet->Area()/TMath::Pi())));
649  fhNumberOfJetTracks_Det->Fill(JetDet->GetNumberOfTracks());
650  if((JetTrue = JetDet->ClosestJet())){
651  JetsMatched=kTRUE;
652  JetCounter_True++;
653  fhJetPt_True->Fill(JetTrue->Pt());
654  JetPhiTrue=JetTrue->Phi();
655  if(JetPhiTrue < -1*TMath::Pi()) JetPhiTrue += (2*TMath::Pi());
656  else if (JetPhiTrue > TMath::Pi()) JetPhiTrue -= (2*TMath::Pi());
657  fhJetPhi_True->Fill(JetPhiTrue);
658  fhJetEta_True->Fill(JetTrue->Eta());
659  fhJetMass_True->Fill(JetTrue->M());
660  fhJetRadius_True->Fill(TMath::Sqrt((JetTrue->Area()/TMath::Pi())));
661  fhNumberOfJetTracks_True->Fill(JetTrue->GetNumberOfTracks());
662  fh2PtRatio->Fill(JetDet->Pt(),JetTrue->Pt(),Pythia_Event_Weight);
663  }
664  else continue;
665  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) fJetInfoVar[0]= JetDet->Pt()-(GetRhoVal(0)*JetDet->Area());
666  else fJetInfoVar[0]=JetDet->Pt();
667  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) fJetInfoVar[1]= JetTrue->Pt()-(GetRhoVal(0)*JetTrue->Area());
668  else fJetInfoVar[1]=JetTrue->Pt();
669  if(fDoSoftDrop) {
670  SoftDrop(JetDet,JetContDet,fZCut,fBeta_SD,kFALSE);
671  SoftDrop(JetTrue,JetContTrue,fZCut,fBeta_SD,kTRUE);
672  }
673  else{
674  fJetInfoVar[2]=0;
675  fJetInfoVar[3]=0;
676  fJetInfoVar[4]=0;
677  fJetInfoVar[5]=0;
678  fJetInfoVar[6]=0;
679  fJetInfoVar[7]=0;
680  fJetInfoVar[8]=0;
681  fJetInfoVar[9]=0;
682  fJetInfoVar[10]=0;
683  fJetInfoVar[11]=0;
684  fJetInfoVar[12]=0;
685  fJetInfoVar[13]=0;
686  fJetInfoVar[14]=0;
687  fJetInfoVar[15]=0;
688  fJetInfoVar[16]=0;
689  fJetInfoVar[17]=0;
690  fJetInfoVar[18]=0;
691  fJetInfoVar[19]=0;
692  }
693  fJetInfoVar[4]=JetDet->M();
694  fJetInfoVar[5]=JetTrue->M();
695  fTreeJetInfo->Fill();
696 
697  JetsMatched=kFALSE;
698  }
699  }
700  fhJetCounter_Det->Fill(JetCounter_Det); //Number of Jets in Each Event Particle Level
701  fhJetCounter_True->Fill(JetCounter_True); //Number of Jets in Each Event Detector Level
702  }
703  }
704 
705 
707  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kGenOnTheFly){
708  AliEmcalJet *Jet1 = NULL; //Original Jet in the event
709  AliJetContainer *JetCont= GetJetContainer(0); //Jet Container for event
710  Int_t JetCounter=0; //Counts number of jets in event
711  Double_t JetPhi=0;
712  Bool_t EventCounter=kFALSE;
713  Double_t JetPt_ForThreshold=0;
714  AliEmcalJet *Jet2= NULL;
715  if(JetCont) {
716  JetCont->ResetCurrentID();
717  while((Jet1=JetCont->GetNextAcceptJet())) {
718  if(!Jet1) continue;
719  if (fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) JetPt_ForThreshold = Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
720  else JetPt_ForThreshold = Jet1->Pt();
721  if(JetPt_ForThreshold<fPtThreshold) {
722  continue;
723  }
724  else {
725  Float_t RecoilDeltaPhi = 0.;
726  if (fJetSelection == kRecoil){
727  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), Jet1->Phi());
728  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
729  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), Jet1->Pt());
730  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), Jet1->Phi()));
731  }
732  if (!EventCounter){
733  EventCounter=kTRUE;
734  }
735  JetCounter++;
736  fhJetPt->Fill(Jet1->Pt());
737  fhJetArea->Fill(Jet1->Area());
738  JetPhi=Jet1->Phi();
739  if(JetPhi < -1*TMath::Pi()) JetPhi += (2*TMath::Pi());
740  else if (JetPhi > TMath::Pi()) JetPhi -= (2*TMath::Pi());
741  fhJetPhi->Fill(JetPhi);
742  fhJetEta->Fill(Jet1->Eta());
743  fhJetMass->Fill(Jet1->M());
744  fhJetRadius->Fill(TMath::Sqrt((Jet1->Area()/TMath::Pi()))); //Radius of Jets per event
745  fhNumberOfJetTracks->Fill(Jet1->GetNumberOfTracks());
746  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) fJetInfoVar[0]= Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
747  else fJetInfoVar[0]=Jet1->Pt();
748  fJetInfoVar[1]=0;
749  if(fDoSoftDrop) {
750  SoftDrop(Jet1,JetCont,fZCut,fBeta_SD,kFALSE);
751  SoftDrop(Jet1,JetCont,fZCut,fBeta_SD,kTRUE);
752  }
753  else{
754  fJetInfoVar[2]=0;
755  fJetInfoVar[3]=0;
756  fJetInfoVar[4]=0;
757  fJetInfoVar[5]=0;
758  fJetInfoVar[6]=0;
759  fJetInfoVar[7]=0;
760  fJetInfoVar[8]=0;
761  fJetInfoVar[9]=0;
762  fJetInfoVar[10]=0;
763  fJetInfoVar[11]=0;
764  fJetInfoVar[12]=0;
765  fJetInfoVar[13]=0;
766  fJetInfoVar[14]=0;
767  fJetInfoVar[15]=0;
768  fJetInfoVar[16]=0;
769  fJetInfoVar[17]=0;
770  fJetInfoVar[18]=0;
771  fJetInfoVar[19]=0;
772  }
773  fJetInfoVar[4]=Jet1->M();
774  fJetInfoVar[5]=0;
775  fTreeJetInfo->Fill();
776  }
777  }
778  fhJetCounter->Fill(JetCounter); //Number of Jets in Each Event
779  }
780  }
781 
782 
783 
784  return kTRUE;
785 }
786 
787 //________________________________________________________________________
788 Double_t AliAnalysisTaskRecoilJetYield::RelativePhiEventPlane(Double_t EventPlane, Double_t Phi){
789 
790  if(Phi < -1*TMath::Pi()) Phi += (2*TMath::Pi());
791  else if (Phi > TMath::Pi()) Phi -= (2*TMath::Pi());
792  Double_t DeltaPhi=Phi-EventPlane;
793  if(DeltaPhi < -1*TMath::Pi()) DeltaPhi += (2*TMath::Pi());
794  else if (DeltaPhi > TMath::Pi()) DeltaPhi -= (2*TMath::Pi());
795  return DeltaPhi;
796 }
797 //________________________________________________________________________
798 Double_t AliAnalysisTaskRecoilJetYield::RelativePhi(Double_t Phi1, Double_t Phi2){
799 
800  if(Phi1 < -1*TMath::Pi()) Phi1 += (2*TMath::Pi());
801  else if (Phi1 > TMath::Pi()) Phi1 -= (2*TMath::Pi());
802  if(Phi2 < -1*TMath::Pi()) Phi2 += (2*TMath::Pi());
803  else if (Phi2 > TMath::Pi()) Phi2 -= (2*TMath::Pi());
804  Double_t DeltaPhi=Phi2-Phi1;
805  if(DeltaPhi < -1*TMath::Pi()) DeltaPhi += (2*TMath::Pi());
806  else if (DeltaPhi > TMath::Pi()) DeltaPhi -= (2*TMath::Pi());
807  return DeltaPhi;
808 }
809 
810 
811 
812 
813 //--------------------------------------------------------------------------
814 Double_t AliAnalysisTaskRecoilJetYield::Angularity(AliEmcalJet *Jet, Int_t JetContNb){
815 
816  AliJetContainer *JetCont = GetJetContainer(JetContNb);
817  Double_t Angularity_Numerator=0;
818  Double_t Angularity_Denominator=0;
819  AliVParticle *Particle=0x0;
820  Double_t DeltaPhi=0;
821 
822 
823  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks
824  Particle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
825  if(!Particle) continue;
826  DeltaPhi=RelativePhi(Jet->Phi(),Particle->Phi());
827  Angularity_Numerator=Angularity_Numerator+(Particle->Pt()*TMath::Sqrt(((Particle->Eta()-Jet->Eta())*(Particle->Eta()-Jet->Eta()))+(DeltaPhi*DeltaPhi)));
828  Angularity_Denominator= Angularity_Denominator+Particle->Pt();
829  }
830  if(Angularity_Denominator!=0) return Angularity_Numerator/Angularity_Denominator;
831  else return -1;
832 
833 }
834 
835 
836 
837 //--------------------------------------------------------------------------
838 Double_t AliAnalysisTaskRecoilJetYield::PTD(AliEmcalJet *Jet, Int_t JetContNb){
839 
840  AliJetContainer *JetCont = GetJetContainer(JetContNb);
841  Double_t PTD_Numerator=0; //Reset these values
842  Double_t PTD_Denominator=0;
843  AliVParticle *Particle=0x0;
844  Double_t DeltaPhi=0;
845  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks
846  Particle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
847  if(!Particle) continue;
848  PTD_Numerator=PTD_Numerator+(Particle->Pt()*Particle->Pt());
849  PTD_Denominator=PTD_Denominator+Particle->Pt();
850  }
851  if(PTD_Denominator!=0) return TMath::Sqrt(PTD_Numerator)/PTD_Denominator;
852  else return -1;
853 
854 }
855 
856 //_________________________________________________________________________
857  void AliAnalysisTaskRecoilJetYield::SoftDrop(AliEmcalJet *fJet,AliJetContainer *fJetCont, double zcut, double beta, Bool_t fTruthJet){
858  std::vector<fastjet::PseudoJet> fInputVectors;
859  Double_t JetInvMass=0, PseudJetInvMass=0, TrackMom = 0, TrackEnergy = 0;
860  AliParticleContainer *fTrackCont = fJetCont->GetParticleContainer();
861  Double_t JetEta=fJet->Eta(),JetPhi=fJet->Phi();
862  Double_t FJTrackEta[9999],FJTrackPhi[9999],FJTrackPt[9999],EmcalJetTrackEta[9999],EmcalJetTrackPhi[9999],EmcalJetTrackPt[9999];
863  UShort_t FJNTracks=0,EmcalJetNTracks=0;
864  if (fTrackCont) for (Int_t i=0; i<fJet->GetNumberOfTracks(); i++) {
865  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
866  JetInvMass += fTrk->M();
867  if (!fTrk) continue;
868  fastjet::PseudoJet PseudoTracks(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
869  TrackMom += TMath::Sqrt(TMath::Power(fTrk->Px(),2)+TMath::Power(fTrk->Py(),2)+TMath::Power(fTrk->Pz(),2));
870  TrackEnergy += fTrk->E();
871  PseudoTracks.set_user_index(fJet->TrackAt(i)+100);
872  PseudJetInvMass += PseudoTracks.m();
873  fInputVectors.push_back(PseudoTracks);
874  EmcalJetTrackEta[i]=fTrk->Eta();
875  EmcalJetTrackPhi[i]=fTrk->Phi();
876  EmcalJetTrackPt[i]=fTrk->Pt();
877  EmcalJetNTracks++;
878  }
879  fastjet::JetDefinition *fJetDef;
880  fastjet::ClusterSequence *fClustSeqSA;
881 
882 
883  fJetDef = new fastjet::JetDefinition(fastjet::antikt_algorithm, fJetRadius*2, static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
884 
885  try {
886  fClustSeqSA = new fastjet::ClusterSequence(fInputVectors, *fJetDef);
887  } catch (fastjet::Error) {
888  AliError(" [w] FJ Exception caught.");
889  //return -1;
890  }
891 
892  std::vector<fastjet::PseudoJet> fOutputJets;
893  fOutputJets.clear();
894  fOutputJets=fClustSeqSA->inclusive_jets(0);
895 
896  std::vector<fastjet::PseudoJet> jet_constituents = fOutputJets[0].constituents();
897  Float_t NSubjettinessResult[3], NSubBeta = 1, R0=0.4;
898  std::vector<fastjet::PseudoJet> Subjet_Axes;
899  fastjet::PseudoJet SubJet1_Axis,SubJet2_Axis;
900  Double_t DelR=-5;
901 
902 
903  for(Int_t j=1; j<3; j++){
904  if(jet_constituents.size() < j){
905  if(!fTruthJet){
906  //if(j==1) fJetInfoVar[4]=-5;
907  // if(j==2) fJetInfoVar[6]=-5;
908  }
909  else {
910  //if(j==1) fJetInfoVar[5]=-5;
911  // if(j==2) fJetInfoVar[7]=-5;
912  }
913  continue;
914  }
915  else{
916  fastjet::contrib::Nsubjettiness nSub(j,fastjet::contrib::KT_Axes(),fastjet::contrib::NormalizedMeasure(NSubBeta,R0));
917  NSubjettinessResult[j] = nSub.result(fOutputJets[0]);
918  if(j==2){ //find subjet axis to calculate delR
919  Subjet_Axes = nSub.currentAxes();
920  SubJet1_Axis = Subjet_Axes[0];
921  SubJet2_Axis = Subjet_Axes[1];
922 
923  Double_t SubJet1_Eta=SubJet1_Axis.pseudorapidity();
924  Double_t SubJet2_Eta=SubJet2_Axis.pseudorapidity();
925  Double_t SubJet1_Phi=SubJet1_Axis.phi();
926  if(SubJet1_Phi < -1*TMath::Pi()) SubJet1_Phi += (2*TMath::Pi());
927  else if (SubJet1_Phi > TMath::Pi()) SubJet1_Phi -= (2*TMath::Pi());
928  Double_t SubJet2_Phi=SubJet2_Axis.phi();
929  if(SubJet1_Phi < -1*TMath::Pi()) SubJet1_Phi += (2*TMath::Pi());
930  else if (SubJet1_Phi > TMath::Pi()) SubJet1_Phi -= (2*TMath::Pi());
931 
932  Double_t DeltaPhiSubJets,DeltaEtaSubJets;
933  DeltaPhiSubJets=SubJet1_Phi-SubJet2_Phi;
934  DeltaEtaSubJets=SubJet1_Eta-SubJet2_Eta;
935  if(DeltaPhiSubJets < -1*TMath::Pi()) DeltaPhiSubJets += (2*TMath::Pi());
936  else if (DeltaPhiSubJets > TMath::Pi()) DeltaPhiSubJets -= (2*TMath::Pi());
937 
938  DelR = TMath::Power(TMath::Power(DeltaPhiSubJets,2)+TMath::Power(DeltaEtaSubJets,2),0.5);
939  }
940  }
941  }
942  if(!fTruthJet){
943  //fJetInfoVar[4]=NSubjettinessResult[1];
944  //fJetInfoVar[6]=NSubjettinessResult[2];
945  // fJetInfoVar[12]=DelR;
946  }
947  else {
948  //fJetInfoVar[5]=NSubjettinessResult[1];
949  //fJetInfoVar[7]=NSubjettinessResult[2];
950  // fJetInfoVar[13]=DelR;
951  }
952 
953 
954 
955  fastjet::contrib::SoftDrop softdrop(beta, zcut);
956  //fastjet::contrib::SoftDrop softdrop_antikt(beta,zcut);
957  softdrop.set_verbose_structure(kTRUE);
958  //fastjet::JetDefinition jet_def_akt(fastjet::antikt_algorithm, 0.4);
959  // fastjet::contrib::Recluster *antiKT_Recluster(jet_def_akt);
960  fastjet::contrib::Recluster *recluster;
961  if(fReclusterAlgo == 2) recluster = new fastjet::contrib::Recluster(fastjet::kt_algorithm,1,true);
962  if(fReclusterAlgo == 1) recluster = new fastjet::contrib::Recluster(fastjet::antikt_algorithm,1,true);
963  if(fReclusterAlgo == 0) recluster = new fastjet::contrib::Recluster(fastjet::cambridge_algorithm,1,true);
964  softdrop.set_reclustering(true,recluster);
965  fastjet::PseudoJet finaljet = softdrop(fOutputJets[0]);
966  // fastjet::PseudoJet finaljet_antikt = softdrop_antikt(fOutputJets[0]);
967  //cout<< finaljet.structure_of<fastjet::contrib::SoftDrop>().symmetry()<<endl;
968  //cout<< finaljet_antikt.structure_of<fastjet::contrib::SoftDrop>().symmetry()<<endl;
969 
970 
971  AliEmcalJet* jet = new AliEmcalJet(finaljet.perp(), finaljet.eta(), finaljet.phi(), finaljet.m());
972  std::vector<fastjet::PseudoJet> fSDTracks=finaljet.constituents();
973  Double_t FastjetTrackDelR,EmcalTrackDelR;
974  for(Int_t i=0;i<fJet->GetNumberOfConstituents();i++){
975  if(i<=finaljet.constituents().size()){
976  FastjetTrackDelR = TMath::Sqrt(TMath::Power(fSDTracks[i].eta()-JetEta,2)+TMath::Power(fSDTracks[i].phi()-JetPhi,2));
977  FJTrackEta[i]=fSDTracks[i].eta();
978  FJTrackPhi[i]=fSDTracks[i].phi();
979  FJTrackPt[i]=fSDTracks[i].perp();
980  FJNTracks++;
981  }
982  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
983  EmcalTrackDelR = TMath::Sqrt(TMath::Power(fTrk->Eta()-JetEta,2)+TMath::Power(fTrk->Phi()-JetPhi,2));
984  }
985  Int_t NDroppedTracks = fJet->GetNumberOfTracks()-finaljet.constituents().size();
986  Int_t nConstituents(fClustSeqSA->constituents(finaljet).size());
987  jet->SetNumberOfTracks(nConstituents);
988  Double_t SymParam, Mu, DeltaR;
989  SymParam=(finaljet.structure_of<fastjet::contrib::SoftDrop>().symmetry());
990  Mu=(finaljet.structure_of<fastjet::contrib::SoftDrop>().mu());
991  DeltaR=(finaljet.structure_of<fastjet::contrib::SoftDrop>().delta_R());
992  //fhGroomedPtvJetPt->Fill(finaljet.perp(),fJet->Pt());
993  //fhDroppedBranches->Fill(finaljet.structure_of<fastjet::contrib::SoftDrop>().dropped_count());
994  std::vector<fastjet::PseudoJet> subjets;
995  if ( finaljet.has_pieces() ) {
996  subjets = finaljet.pieces();
997  fastjet::PseudoJet subjet1 = subjets[0];
998  fastjet::PseudoJet subjet2 = subjets[1];
999  if(!fTruthJet){
1000  if(subjets[0].perp() > subjets[1].perp()){
1001  fJetInfoVar[6]=subjets[1].m();
1002  fJetInfoVar[8]=LeadingTrackPt(subjets[1]);
1003  fJetInfoVar[10]=subjets[1].perp();
1004  }
1005  else {
1006  fJetInfoVar[6]= subjets[0].m();
1007  fJetInfoVar[8]=LeadingTrackPt(subjets[0]);
1008  fJetInfoVar[10]= subjets[0].perp();
1009  }
1010  }
1011  else {
1012  if(subjets[0].perp() > subjets[1].perp()){
1013  fJetInfoVar[7]=subjets[1].m();
1014  fJetInfoVar[9]=LeadingTrackPt(subjets[1]);
1015  fJetInfoVar[11]=subjets[1].perp();
1016  }
1017  else{
1018  fJetInfoVar[7]= subjets[0].m();
1019  fJetInfoVar[9]=LeadingTrackPt(subjets[0]);
1020  fJetInfoVar[11]=subjets[0].perp();
1021  }
1022  }
1023 
1024  }
1025  else {
1026  if(!fTruthJet){
1027  fJetInfoVar[6]=0;
1028  fJetInfoVar[10]=-1;
1029  }
1030  else {
1031  fJetInfoVar[7]=0;
1032  fJetInfoVar[11]=-1;
1033  }
1034  }
1035 
1036  if(!fTruthJet) fJetInfoVar[2]=SymParam;
1037  else fJetInfoVar[3]=SymParam;
1038  if(!fTruthJet) fJetInfoVar[12] = DeltaR;
1039  else fJetInfoVar[13] = DeltaR;
1040  if(!fTruthJet) fJetInfoVar[14]=finaljet.structure_of<fastjet::contrib::SoftDrop>().dropped_count();
1041  else fJetInfoVar[15]=finaljet.structure_of<fastjet::contrib::SoftDrop>().dropped_count();
1042  if(!(fJetShapeSub==kConstSub)){
1043  if(!fTruthJet) fJetInfoVar[16]=(finaljet.perp()-(GetRhoVal(0)*fJet->Area()));
1044  else fJetInfoVar[17]=finaljet.perp();
1045  }
1046  else{
1047  if(!fTruthJet) fJetInfoVar[16]=(finaljet.perp());
1048  else fJetInfoVar[17]=(finaljet.perp());
1049  }
1050  if(!fTruthJet) fJetInfoVar[18]=(finaljet.m());
1051  else fJetInfoVar[19]=(finaljet.m());
1052 
1053  return;
1054 
1055 
1056 }
1057 
1058 //--------------------------------------------------------------------------
1059 Double_t AliAnalysisTaskRecoilJetYield::LeadingTrackPt(fastjet::PseudoJet jet){
1060  std::vector< fastjet::PseudoJet > constituents = jet.constituents();
1061  fastjet::PseudoJet leadingtrack = constituents[0];
1062  Double_t track_pt;
1063  for(size_t i=0; i<constituents.size(); i++){
1064  track_pt=constituents[i].perp();
1065  if (track_pt > leadingtrack.perp()){
1066  leadingtrack = constituents[i];
1067  }
1068  }
1069  return leadingtrack.perp();
1070 }
1071 
1072 
1073 //--------------------------------------------------------------------------
1074 Int_t AliAnalysisTaskRecoilJetYield::SelectTriggerHadron(Float_t PtMin, Float_t PtMax){
1075 
1076  AliTrackContainer *PartCont = NULL;
1077  AliParticleContainer *PartContMC = NULL;
1078 
1079 
1080  if (fJetShapeSub==kConstSub){
1081  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kGenOnTheFly) PartContMC = GetParticleContainer(1);
1082  else PartCont = GetTrackContainer(1);
1083  }
1084  else{
1085  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kGenOnTheFly) PartContMC = GetParticleContainer(0);
1086  else PartCont = GetTrackContainer(0);
1087  }
1088 
1089  TClonesArray *TracksArray = NULL;
1090  TClonesArray *TracksArrayMC = NULL;
1091  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kGenOnTheFly) TracksArrayMC = PartContMC->GetArray();
1092  else TracksArray = PartCont->GetArray();
1093 
1094  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kGenOnTheFly){
1095  if(!PartContMC || !TracksArrayMC) return -99999;
1096  }
1097  else {
1098  if(!PartCont || !TracksArray) return -99999;
1099  }
1100 
1101  AliAODTrack *Track = 0x0;
1102  Int_t Trigger_Index[100];
1103  for (Int_t i=0; i<100; i++) Trigger_Index[i] = 0;
1104  Int_t Trigger_Counter = 0;
1105  Int_t NTracks=0;
1106  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kGenOnTheFly) NTracks = TracksArrayMC->GetEntriesFast();
1107  else NTracks = TracksArray->GetEntriesFast();
1108  for(Int_t i=0; i < NTracks; i++){
1109  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kGenOnTheFly){
1110  if((Track = static_cast<AliAODTrack*>(PartContMC->GetAcceptParticle(i)))){
1111  if (!Track) continue;
1112  if(TMath::Abs(Track->Eta())>0.9) continue;
1113  if (Track->Pt()<0.15) continue;
1114  if ((Track->Pt() >= PtMin) && (Track->Pt()< PtMax)) {
1115  Trigger_Index[Trigger_Counter] = i;
1116  Trigger_Counter++;
1117  }
1118  }
1119  }
1120  else{
1121  if((Track = static_cast<AliAODTrack*>(PartCont->GetAcceptTrack(i)))){
1122  if (!Track) continue;
1123  if(TMath::Abs(Track->Eta())>0.9) continue;
1124  if (Track->Pt()<0.15) continue;
1125  if ((Track->Pt() >= PtMin) && (Track->Pt()< PtMax)) {
1126  Trigger_Index[Trigger_Counter] = i;
1127  Trigger_Counter++;
1128  }
1129  }
1130  }
1131  }
1132  if (Trigger_Counter == 0) return -99999;
1133  Int_t RandomNumber = 0, Index = 0 ;
1134  TRandom3* Random = new TRandom3(0);
1135  RandomNumber = Random->Integer(Trigger_Counter);
1136  Index = Trigger_Index[RandomNumber];
1137  return Index;
1138 }
1139 
1140 
1141 //________________________________________________________________________
1142 Double_t AliAnalysisTaskRecoilJetYield::GetFractionSharedPt_SubMatching(const AliEmcalJet *jet1, AliParticleContainer *cont2) const
1143 {
1144  AliEmcalJet *jet2 = jet1->ClosestJet();
1145  if (!jet2) return -1;
1146 
1147  Double_t jetPt2 = jet2->Pt();
1148  if (jetPt2 <= 0) return -1;
1149 
1150  Int_t bgeom = kTRUE;
1151  if (!cont2) bgeom = kFALSE;
1152  Double_t sumPt = 0.;
1153  AliVParticle *vpf = 0x0;
1154  Int_t iFound = 0;
1155  for (Int_t icc = 0; icc < jet2->GetNumberOfTracks(); icc++) {
1156  Int_t idx = (Int_t)jet2->TrackAt(icc);
1157  //get particle
1158  AliVParticle *p2 = 0x0;
1159  if (bgeom) p2 = static_cast<AliVParticle*>(jet2->TrackAt(icc, cont2->GetArray()));
1160  iFound = 0;
1161  for (Int_t icf = 0; icf < jet1->GetNumberOfTracks(); icf++) {
1162  if (!bgeom && idx == jet1->TrackAt(icf) && iFound == 0 ) {
1163  iFound = 1;
1164  vpf = jet1->Track(icf);
1165  if (vpf) sumPt += vpf->Pt();
1166  continue;
1167  }
1168  if (bgeom){
1169  vpf = jet1->Track(icf);
1170  if (!vpf) continue;
1171  if (!SameParticle(vpf, p2, 1.e-4)) continue; //not the same particle
1172  sumPt += vpf->Pt();
1173  }
1174  }
1175  }
1176 
1177  Double_t fraction = sumPt / jetPt2;
1178 
1179  return fraction;
1180 }
1181 
1182 //________________________________________________________________________
1183 Bool_t AliAnalysisTaskRecoilJetYield::SameParticle(const AliVParticle* part1, const AliVParticle* part2, Double_t dist)
1184 {
1185  if(!part1) return kFALSE;
1186  if(!part2) return kFALSE;
1187  Double_t dPhi = TMath::Abs(part1->Phi() - part2->Phi());
1188  dPhi = TVector2::Phi_mpi_pi(dPhi);
1189  if (dPhi > dist) return kFALSE;
1190  return kTRUE;
1191 }
1192 
1193 //________________________________________________________________________
1194 Bool_t AliAnalysisTaskRecoilJetYield::RetrieveEventObjects() {
1195  //
1196  // retrieve event objects
1197  //
1198  if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
1199  return kFALSE;
1200 
1201  return kTRUE;
1202 }
1203 
1204 
1205 //_______________________________________________________________________
1206 void AliAnalysisTaskRecoilJetYield::Terminate(Option_t *)
1207 {
1208  // Called once at the end of the analysis.
1209 
1210 
1211 }
AliAnalysisTaskRecoilJetYield::RetrieveEventObjects
Bool_t RetrieveEventObjects()
Retrieve common objects from event.
Definition: AliAnalysisTaskRecoilJetYield.cxx:1194
AliEmcalJet::Area
Double_t Area() const
Definition: AliEmcalJet.h:130
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:838
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:945
TH2F
Definition: External.C:236
AliEmcalJet::ClosestJet
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:327
AliAnalysisTaskEmcalJet::GetJetContainer
AliJetContainer * GetJetContainer(Int_t i=0) const
Definition: AliAnalysisTaskEmcalJet.cxx:377
AliAnalysisTaskRecoilJetYield::FillHistograms
Bool_t FillHistograms()
Function filling histograms.
Definition: AliAnalysisTaskRecoilJetYield.cxx:371
AliEmcalJet::Eta
Double_t Eta() const
Definition: AliEmcalJet.h:121
Phi_Upper
Double_t Phi_Upper
Definition: AliAnalysisTaskRecoilJetYield.cxx:55
AliEmcalJet::Phi
Double_t Phi() const
Definition: AliEmcalJet.h:117
AliTrackContainer
Container with name, TClonesArray and cuts for particles.
Definition: AliTrackContainer.h:54
AliAnalysisTaskEmcal::fEPV0
Double_t fEPV0
!event plane V0
Definition: AliAnalysisTaskEmcal.h:888
AliEmcalJet::GetLabel
Int_t GetLabel() const
Definition: AliEmcalJet.h:124
AliEmcalParticle.h
AliEmcalJet::Track
AliVParticle * Track(Int_t idx) const
Definition: AliEmcalJet.cxx:690
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:857
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:140
AliParticleContainer
Container for particles within the EMCAL framework.
Definition: AliParticleContainer.h:31
AliEmcalJet::TrackAt
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:160
AliEmcalJet::GetNumberOfTracks
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:139
kData
TString kData
Declare data MC or deltaAOD.
Definition: AddTaskCaloTrackCorr.C:63
AliAnalysisTaskEmcal::GetParticleContainer
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
Definition: AliAnalysisTaskEmcal.cxx:1623
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:1074
AliAnalysisTaskRecoilJetYield::GetFractionSharedPt_SubMatching
Double_t GetFractionSharedPt_SubMatching(const AliEmcalJet *jet, AliParticleContainer *cont2=0x0) const
Definition: AliAnalysisTaskRecoilJetYield.cxx:1142
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:886
nBranch
const Int_t nBranch
Definition: AliAnalysisTaskRecoilJetYield.h:24
AliAnalysisTaskRecoilJetYield::Angularity
Double_t Angularity(AliEmcalJet *Jet, Int_t JetContNb)
Definition: AliAnalysisTaskRecoilJetYield.cxx:814
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:798
AliAnalysisTaskRecoilJetYield::SameParticle
static Bool_t SameParticle(const AliVParticle *part1, const AliVParticle *part2, Double_t dist=1.e-4)
Definition: AliAnalysisTaskRecoilJetYield.cxx:1183
AliEmcalJet::Pt
Double_t Pt() const
Definition: AliEmcalJet.h:109
AliAnalysisTaskEmcalJet::GetRhoVal
Double_t GetRhoVal(Int_t i=0) const
Definition: AliAnalysisTaskEmcalJet.cxx:626
AliAnalysisTaskEmcal::fOutput
AliEmcalList * fOutput
!output list
Definition: AliAnalysisTaskEmcal.h:909
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:252
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:51
AliEmcalJet::SetNumberOfTracks
void SetNumberOfTracks(Int_t n)
Definition: AliEmcalJet.h:266
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:336
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:364
AliAnalysisTaskRecoilJetYield::RelativePhiEventPlane
Double_t RelativePhiEventPlane(Double_t EventPlane, Double_t Phi)
Definition: AliAnalysisTaskRecoilJetYield.cxx:788
AliAnalysisTaskRecoilJetYield::LeadingTrackPt
Double_t LeadingTrackPt(fastjet::PseudoJet jet)
Definition: AliAnalysisTaskRecoilJetYield.cxx:1059
AliEmcalJet::M
Double_t M() const
Definition: AliEmcalJet.h:120
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