AliPhysics  d7f9f85 (d7f9f85)
 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"
46 #include "AliMultSelection.h"
47 
48 //Globals
49 
50 
51 
55 Double_t Phi_Upper=2*(TMath::Pi());
56 Double_t Phi_Lower=(-1*(TMath::Pi()));
59 
60 
61 
62 
63 using std::cout;
64 using std::endl;
65 
67 
68 //________________________________________________________________________
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 //________________________________________________________________________
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  }
223  DefineOutput(1, TList::Class());
224  DefineOutput(2, TTree::Class());
225 }
226 
227 //________________________________________________________________________
229 {
230  // Destructor.
231 }
232 
233 //________________________________________________________________________
235 {
236  // Create user output.
237 
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] = "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  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);
281  fh2PtTriggerHadronJet= new TH2F("fh2PtTriggerHadronJet", "fh2PtTriggerHadronJet",1500,-0.5,149.5,1500,-0.5,149.5);
283  fhPhiTriggerHadronJet= new TH1F("fhPhiTriggerHadronJet", "fhPhiTriggerHadronJet",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
285  fhPhiTriggerHadronEventPlane= new TH1F("fhPhiTriggerHadronEventPlane", "fhPhiTriggerHadronEventPlane",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
287  fhPhiTriggerHadronEventPlaneTPC= new TH1F("fhPhiTriggerHadronEventPlaneTPC", "fhPhiTriggerHadronEventPlaneTPC",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
289 
290 
291  }
292 
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);
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);
315  fhDroppedBranches= new TH1F("fhDroppedBranches","Number of Softdropped branches",50,0,50);
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 );
321 
322 
323  }
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);
345  fhNumberOfJetTracks_Det= new TH1F("fhNumberOfJetTracks_Det", "Number of Tracks within a Jet Detector Level", 300, -0.5,299.5);
347  fhNumberOfJetTracks_True= new TH1F("fhNumberOfJetTracks_True", "Number of Tracks within a Jet Particle Level", 300, -0.5,299.5);
349 
350  fhJetCounter_Det= new TH1F("fhJetCounter_Det", "Jet Counter Detector Level", 150, -0.5, 149.5);
352  fhJetCounter_True= new TH1F("fhJetCounter_True", "Jet Counter Particle Level", 150, -0.5, 149.5);
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 //________________________________________________________________________
365 {
366  // Run analysis code here, if needed. It will be executed before FillHistograms().
367  return kTRUE;
368 }
369 
370 //________________________________________________________________________
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
396  if (TriggerHadronLabel==-99999) return 0; //Trigger Hadron Not Found
397  AliTrackContainer *PartCont =NULL;
398  AliParticleContainer *PartContMC = NULL;
399  if (fJetShapeSub==kConstSub){
401  else PartCont = GetTrackContainer(1);
402  }
403  else {
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 
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  fhJetRadius->Fill(TMath::Sqrt((Jet1->Area()/TMath::Pi()))); //Radius of Jets per event
468  fhNumberOfJetTracks->Fill(Jet1->GetNumberOfTracks());
469  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) fJetInfoVar[0]= Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
470  else fJetInfoVar[0]=Jet1->Pt();
471  fJetInfoVar[1]=0;
472  if(fDoSoftDrop) {
473  SoftDrop(Jet1,JetCont,fZCut,fBeta_SD,kFALSE);
474  //SoftDrop(Jet1,JetCont,fZCut,fBeta_SD,kTRUE);
475 
476  fJetInfoVar[3]=0;
477  fJetInfoVar[5]=0;
478  fJetInfoVar[7]=0;
479  fJetInfoVar[13]=0;
480  fJetInfoVar[15]=0;
481  fJetInfoVar[17]=0;
482  fJetInfoVar[19]=0;
483 
484  }
485  else{
486  fJetInfoVar[2]=0;
487  fJetInfoVar[3]=0;
488  fJetInfoVar[4]=0;
489  fJetInfoVar[5]=0;
490  fJetInfoVar[6]=0;
491  fJetInfoVar[7]=0;
492  fJetInfoVar[12]=0;
493  fJetInfoVar[13]=0;
494  fJetInfoVar[14]=0;
495  fJetInfoVar[15]=0;
496  fJetInfoVar[16]=0;
497  fJetInfoVar[17]=0;
498  fJetInfoVar[18]=0;
499  fJetInfoVar[19]=0;
500  }
501  fJetInfoVar[8]=PTD(Jet1,0);
502  fJetInfoVar[9]=0;
503  fJetInfoVar[10]=Angularity(Jet1,0);
504  fJetInfoVar[11]=0;
505  fTreeJetInfo->Fill();
506  }
507  }
508  fhJetCounter->Fill(JetCounter); //Number of Jets in Each Event
509  }
510  }
511 
513  if(fJetShapeType == kDetEmbPart){
514  AliEmcalJet *JetHybridS = NULL; //Subtracted hybrid Jet
515  AliEmcalJet *JetHybridUS = NULL; //Unsubtracted Hybrid Jet //For matching SubtractedHybrid->DetPythia this jet container is also Subtracted Hybrid
516  AliEmcalJet *JetPythDet = NULL; //Detector Level Pythia Jet
517  AliEmcalJet *JetPythTrue = NULL; //Particle Level Pyhtia Jet
518  AliJetContainer *JetContHybridS= GetJetContainer(0); //Jet Container for Subtracted Hybrid Jets
519  AliJetContainer *JetContHybridUS= GetJetContainer(1); //Jet Container for Unsubtracted Hybrid Jets
520  AliJetContainer *JetContPythDet= GetJetContainer(2); //Jet Container for Detector Level Pyhtia Jets
521  AliJetContainer *JetContPythTrue= GetJetContainer(3); //Jet Container for Particle Level Pythia Jets
522 
523 
524 
525  Bool_t JetsMatched = kFALSE;
526  Double_t JetPtThreshold;
527  JetContHybridS->ResetCurrentID();
528  JetContHybridUS->ResetCurrentID();
529  JetContPythDet->ResetCurrentID();
530  JetContPythTrue->ResetCurrentID();
531 
532  while((JetPythDet = JetContPythDet->GetNextAcceptJet())){
533  fhDetJetPt_Incl->Fill(JetPythDet->Pt()); //Fill histogram with all Detector level jets
534  }
535 
536 
537  while((JetHybridS = JetContHybridS->GetNextAcceptJet())){
538  if (fJetShapeSub==kConstSub) JetPtThreshold=JetHybridS->Pt();
539  else JetPtThreshold=JetHybridS->Pt()-(GetRhoVal(0)*JetHybridS->Area());
540  if ( (!JetHybridS) || (JetPtThreshold<fPtThreshold)) continue;
541  Int_t JetNumber=-1;
542  for(Int_t i = 0; i<JetContHybridUS->GetNJets(); i++) {
543  JetHybridUS = JetContHybridUS->GetJet(i);
544  if (!JetHybridUS) continue;
545 
546  if(JetHybridUS->GetLabel()==JetHybridS->GetLabel()) {
547  JetNumber=i;
548  }
549  }
550  if(JetNumber==-1) continue;
551  JetHybridUS=JetContHybridUS->GetJet(JetNumber);
552  if (JetContHybridUS->AliJetContainer::GetFractionSharedPt(JetHybridUS)<fSharedFractionPtMin && !fSubMatching) {
553  continue;
554  }
556  continue;
557  }
558  JetPythDet=JetHybridUS->ClosestJet();
559  if (!JetHybridUS) {
560  Printf("Unsubtracted embedded jet does not exist, returning");
561  continue;
562  }
563  if (!JetPythDet) continue;
564  UInt_t rejectionReason = 0;
565  if (!(JetContPythDet->AcceptJet(JetPythDet,rejectionReason))) continue;
566  fhDetJetPt_Matched->Fill(JetPythDet->Pt()); //Fill only matched detector level jets for tagging efficiency comparison
567  JetPythTrue=JetPythDet->ClosestJet();
568  if(!JetPythTrue) continue;
569  JetsMatched=kTRUE;
570 
571  if(fJetShapeSub==kConstSub) fJetInfoVar[0]=JetHybridS->Pt();
572  else fJetInfoVar[0]=JetHybridS->Pt()-(GetRhoVal(0)*JetHybridS->Area());
573  fJetInfoVar[1]=JetPythTrue->Pt();
574  if(fDoSoftDrop) {
575  SoftDrop(JetHybridS,JetContHybridS,fZCut,fBeta_SD,kFALSE);
576  SoftDrop(JetPythTrue,JetContPythTrue,fZCut,fBeta_SD,kTRUE);
577  }
578  else{
579  fJetInfoVar[2]=0;
580  fJetInfoVar[3]=0;
581  fJetInfoVar[4]=0;
582  fJetInfoVar[5]=0;
583  fJetInfoVar[6]=0;
584  fJetInfoVar[7]=0;
585  fJetInfoVar[12]=0;
586  fJetInfoVar[13]=0;
587  fJetInfoVar[14]=0;
588  fJetInfoVar[15]=0;
589  fJetInfoVar[16]=0;
590  fJetInfoVar[17]=0;
591  fJetInfoVar[18]=0;
592  fJetInfoVar[19]=0;
593  }
594  fJetInfoVar[8]=PTD(JetHybridS,0);
595  fJetInfoVar[9]=PTD(JetPythTrue,0);
596  fJetInfoVar[10]=Angularity(JetHybridS,0);
597  fJetInfoVar[11]=Angularity(JetPythTrue,0);
598  fTreeJetInfo->Fill();
599  }
600 
601  }
602 
604  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kTrueDet){ //Truth->Detector response
605  AliEmcalJet *JetDet = NULL; //Detector Level Jet
606  AliEmcalJet *JetTrue = NULL; //Particle Level Jet
607  AliJetContainer *JetContDet= GetJetContainer(0); //Jet Container for Detector Level Pythia
608  AliJetContainer *JetContTrue= GetJetContainer(1); //Jet Container for Particle Level Pythia
609  Double_t JetPhiDet=0;
610  Double_t JetPhiTrue=0;
611  Bool_t JetsMatched=kFALSE;
612  Double_t Pythia_Event_Weight=1;
613  Bool_t EventCounter=kFALSE;
614  Int_t JetCounter_Det=0,JetCounter_True=0;
615  JetContDet->ResetCurrentID();
616  JetContTrue->ResetCurrentID();
617  if(JetContDet) {
618  while((JetDet=JetContDet->GetNextAcceptJet())) {
619  if( (!JetDet) || ((JetDet->Pt())<fPtThreshold)) {
620  // fhEventCounter_1->Fill(3); //events where the jet had a problem
621  continue;
622  }
623  else {
624  /* if(fSemigoodCorrect){
625  Double_t HoleDistance=RelativePhi(Jet1->Phi(),fHolePos);
626  if(TMath::Abs(HoleDistance)<fHoleWidth) continue;
627  }*/
628  Float_t RecoilDeltaPhi = 0.;
629  if (fJetSelection == kRecoil){
630  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), JetDet->Phi());
631  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
632  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), JetDet->Pt());
633  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), JetDet->Phi()));
634  }
635  JetCounter_Det++;
636  fhJetPt_Det->Fill(JetDet->Pt());
637  JetPhiDet=JetDet->Phi();
638  if(JetPhiDet < -1*TMath::Pi()) JetPhiDet += (2*TMath::Pi());
639  else if (JetPhiDet > TMath::Pi()) JetPhiDet -= (2*TMath::Pi());
640  fhJetPhi_Det->Fill(JetPhiDet);
641  fhJetEta_Det->Fill(JetDet->Eta());
642  fhJetMass_Det->Fill(JetDet->M());
643  fhJetRadius_Det->Fill(TMath::Sqrt((JetDet->Area()/TMath::Pi())));
645  if((JetTrue = JetDet->ClosestJet())){
646  JetsMatched=kTRUE;
647  JetCounter_True++;
648  fhJetPt_True->Fill(JetTrue->Pt());
649  JetPhiTrue=JetTrue->Phi();
650  if(JetPhiTrue < -1*TMath::Pi()) JetPhiTrue += (2*TMath::Pi());
651  else if (JetPhiTrue > TMath::Pi()) JetPhiTrue -= (2*TMath::Pi());
652  fhJetPhi_True->Fill(JetPhiTrue);
653  fhJetEta_True->Fill(JetTrue->Eta());
654  fhJetMass_True->Fill(JetTrue->M());
655  fhJetRadius_True->Fill(TMath::Sqrt((JetTrue->Area()/TMath::Pi())));
657  fh2PtRatio->Fill(JetDet->Pt(),JetTrue->Pt(),Pythia_Event_Weight);
658  }
659  else continue;
660  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) fJetInfoVar[0]= JetDet->Pt()-(GetRhoVal(0)*JetDet->Area());
661  else fJetInfoVar[0]=JetDet->Pt();
662  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) fJetInfoVar[1]= JetTrue->Pt()-(GetRhoVal(0)*JetTrue->Area());
663  else fJetInfoVar[1]=JetTrue->Pt();
664  if(fDoSoftDrop) {
665  SoftDrop(JetDet,JetContDet,fZCut,fBeta_SD,kFALSE);
666  SoftDrop(JetTrue,JetContTrue,fZCut,fBeta_SD,kTRUE);
667  }
668  else{
669  fJetInfoVar[2]=0;
670  fJetInfoVar[3]=0;
671  fJetInfoVar[4]=0;
672  fJetInfoVar[5]=0;
673  fJetInfoVar[6]=0;
674  fJetInfoVar[7]=0;
675  fJetInfoVar[12]=0;
676  fJetInfoVar[13]=0;
677  fJetInfoVar[14]=0;
678  fJetInfoVar[15]=0;
679  fJetInfoVar[16]=0;
680  fJetInfoVar[17]=0;
681  fJetInfoVar[18]=0;
682  fJetInfoVar[19]=0;
683  }
684  fJetInfoVar[8]=PTD(JetDet,0);
685  fJetInfoVar[9]=PTD(JetTrue,0);
686  fJetInfoVar[10]=Angularity(JetDet,0);
687  fJetInfoVar[11]=Angularity(JetTrue,0);
688  fTreeJetInfo->Fill();
689 
690  JetsMatched=kFALSE;
691  }
692  }
693  fhJetCounter_Det->Fill(JetCounter_Det); //Number of Jets in Each Event Particle Level
694  fhJetCounter_True->Fill(JetCounter_True); //Number of Jets in Each Event Detector Level
695  }
696  }
697 
698 
701  AliEmcalJet *Jet1 = NULL; //Original Jet in the event
702  AliJetContainer *JetCont= GetJetContainer(0); //Jet Container for event
703  Int_t JetCounter=0; //Counts number of jets in event
704  Double_t JetPhi=0;
705  Bool_t EventCounter=kFALSE;
706  Double_t JetPt_ForThreshold=0;
707  AliEmcalJet *Jet2= NULL;
708  if(JetCont) {
709  JetCont->ResetCurrentID();
710  while((Jet1=JetCont->GetNextAcceptJet())) {
711  if(!Jet1) continue;
712  if (fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) JetPt_ForThreshold = Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
713  else JetPt_ForThreshold = Jet1->Pt();
714  if(JetPt_ForThreshold<fPtThreshold) {
715  continue;
716  }
717  else {
718  Float_t RecoilDeltaPhi = 0.;
719  if (fJetSelection == kRecoil){
720  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), Jet1->Phi());
721  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
722  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), Jet1->Pt());
723  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), Jet1->Phi()));
724  }
725  if (!EventCounter){
726  EventCounter=kTRUE;
727  }
728  JetCounter++;
729  fhJetPt->Fill(Jet1->Pt());
730  fhJetArea->Fill(Jet1->Area());
731  JetPhi=Jet1->Phi();
732  if(JetPhi < -1*TMath::Pi()) JetPhi += (2*TMath::Pi());
733  else if (JetPhi > TMath::Pi()) JetPhi -= (2*TMath::Pi());
734  fhJetPhi->Fill(JetPhi);
735  fhJetEta->Fill(Jet1->Eta());
736  fhJetMass->Fill(Jet1->M());
737  fhJetRadius->Fill(TMath::Sqrt((Jet1->Area()/TMath::Pi()))); //Radius of Jets per event
738  fhNumberOfJetTracks->Fill(Jet1->GetNumberOfTracks());
739  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) fJetInfoVar[0]= Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
740  else fJetInfoVar[0]=Jet1->Pt();
741  fJetInfoVar[1]=0;
742  if(fDoSoftDrop) {
743  SoftDrop(Jet1,JetCont,fZCut,fBeta_SD,kFALSE);
744  SoftDrop(Jet1,JetCont,fZCut,fBeta_SD,kTRUE);
745  }
746  else{
747  fJetInfoVar[2]=0;
748  fJetInfoVar[3]=0;
749  fJetInfoVar[4]=0;
750  fJetInfoVar[5]=0;
751  fJetInfoVar[6]=0;
752  fJetInfoVar[7]=0;
753  fJetInfoVar[12]=0;
754  fJetInfoVar[13]=0;
755  fJetInfoVar[14]=0;
756  fJetInfoVar[15]=0;
757  fJetInfoVar[16]=0;
758  fJetInfoVar[17]=0;
759  fJetInfoVar[18]=0;
760  fJetInfoVar[19]=0;
761  }
762  fJetInfoVar[8]=PTD(Jet1,0);
763  fJetInfoVar[9]=0;
764  fJetInfoVar[10]=Angularity(Jet1,0);
765  fJetInfoVar[11]=0;
766  fTreeJetInfo->Fill();
767  }
768  }
769  fhJetCounter->Fill(JetCounter); //Number of Jets in Each Event
770  }
771  }
772 
773 
774 
775  return kTRUE;
776 }
777 
778 //________________________________________________________________________
780 
781  if(Phi < -1*TMath::Pi()) Phi += (2*TMath::Pi());
782  else if (Phi > TMath::Pi()) Phi -= (2*TMath::Pi());
783  Double_t DeltaPhi=Phi-EventPlane;
784  if(DeltaPhi < -1*TMath::Pi()) DeltaPhi += (2*TMath::Pi());
785  else if (DeltaPhi > TMath::Pi()) DeltaPhi -= (2*TMath::Pi());
786  return DeltaPhi;
787 }
788 //________________________________________________________________________
790 
791  if(Phi1 < -1*TMath::Pi()) Phi1 += (2*TMath::Pi());
792  else if (Phi1 > TMath::Pi()) Phi1 -= (2*TMath::Pi());
793  if(Phi2 < -1*TMath::Pi()) Phi2 += (2*TMath::Pi());
794  else if (Phi2 > TMath::Pi()) Phi2 -= (2*TMath::Pi());
795  Double_t DeltaPhi=Phi2-Phi1;
796  if(DeltaPhi < -1*TMath::Pi()) DeltaPhi += (2*TMath::Pi());
797  else if (DeltaPhi > TMath::Pi()) DeltaPhi -= (2*TMath::Pi());
798  return DeltaPhi;
799 }
800 
801 
802 
803 
804 //--------------------------------------------------------------------------
806 
807  AliJetContainer *JetCont = GetJetContainer(JetContNb);
808  Double_t Angularity_Numerator=0;
809  Double_t Angularity_Denominator=0;
810  AliVParticle *Particle=0x0;
811  Double_t DeltaPhi=0;
812 
813 
814  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks
815  Particle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
816  if(!Particle) continue;
817  DeltaPhi=RelativePhi(Jet->Phi(),Particle->Phi());
818  Angularity_Numerator=Angularity_Numerator+(Particle->Pt()*TMath::Sqrt(((Particle->Eta()-Jet->Eta())*(Particle->Eta()-Jet->Eta()))+(DeltaPhi*DeltaPhi)));
819  Angularity_Denominator= Angularity_Denominator+Particle->Pt();
820  }
821  if(Angularity_Denominator!=0) return Angularity_Numerator/Angularity_Denominator;
822  else return -1;
823 
824 }
825 
826 
827 
828 //--------------------------------------------------------------------------
830 
831  AliJetContainer *JetCont = GetJetContainer(JetContNb);
832  Double_t PTD_Numerator=0; //Reset these values
833  Double_t PTD_Denominator=0;
834  AliVParticle *Particle=0x0;
835  Double_t DeltaPhi=0;
836  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks
837  Particle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
838  if(!Particle) continue;
839  PTD_Numerator=PTD_Numerator+(Particle->Pt()*Particle->Pt());
840  PTD_Denominator=PTD_Denominator+Particle->Pt();
841  }
842  if(PTD_Denominator!=0) return TMath::Sqrt(PTD_Numerator)/PTD_Denominator;
843  else return -1;
844 
845 }
846 
847 //_________________________________________________________________________
848  void AliAnalysisTaskRecoilJetYield::SoftDrop(AliEmcalJet *fJet,AliJetContainer *fJetCont, double zcut, double beta, Bool_t fTruthJet){
849  std::vector<fastjet::PseudoJet> fInputVectors;
850  Double_t JetInvMass=0, PseudJetInvMass=0, TrackMom = 0, TrackEnergy = 0;
851  AliParticleContainer *fTrackCont = fJetCont->GetParticleContainer();
852  Double_t JetEta=fJet->Eta(),JetPhi=fJet->Phi();
853  Double_t FJTrackEta[9999],FJTrackPhi[9999],FJTrackPt[9999],EmcalJetTrackEta[9999],EmcalJetTrackPhi[9999],EmcalJetTrackPt[9999];
854  UShort_t FJNTracks=0,EmcalJetNTracks=0;
855  if (fTrackCont) for (Int_t i=0; i<fJet->GetNumberOfTracks(); i++) {
856  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
857  JetInvMass += fTrk->M();
858  if (!fTrk) continue;
859  fastjet::PseudoJet PseudoTracks(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
860  TrackMom += TMath::Sqrt(TMath::Power(fTrk->Px(),2)+TMath::Power(fTrk->Py(),2)+TMath::Power(fTrk->Pz(),2));
861  TrackEnergy += fTrk->E();
862  PseudoTracks.set_user_index(fJet->TrackAt(i)+100);
863  PseudJetInvMass += PseudoTracks.m();
864  fInputVectors.push_back(PseudoTracks);
865  EmcalJetTrackEta[i]=fTrk->Eta();
866  EmcalJetTrackPhi[i]=fTrk->Phi();
867  EmcalJetTrackPt[i]=fTrk->Pt();
868  EmcalJetNTracks++;
869  }
870  fastjet::JetDefinition *fJetDef;
871  fastjet::ClusterSequence *fClustSeqSA;
872 
873 
874  fJetDef = new fastjet::JetDefinition(fastjet::antikt_algorithm, fJetRadius*2, static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
875 
876  try {
877  fClustSeqSA = new fastjet::ClusterSequence(fInputVectors, *fJetDef);
878  } catch (fastjet::Error) {
879  AliError(" [w] FJ Exception caught.");
880  //return -1;
881  }
882 
883  std::vector<fastjet::PseudoJet> fOutputJets;
884  fOutputJets.clear();
885  fOutputJets=fClustSeqSA->inclusive_jets(0);
886 
887  std::vector<fastjet::PseudoJet> jet_constituents = fOutputJets[0].constituents();
888  Float_t NSubjettinessResult[3], NSubBeta = 1, R0=0.4;
889  std::vector<fastjet::PseudoJet> Subjet_Axes;
890  fastjet::PseudoJet SubJet1_Axis,SubJet2_Axis;
891  Double_t DelR=-5;
892 
893 
894  for(Int_t j=1; j<3; j++){
895  if(jet_constituents.size() < j){
896  if(!fTruthJet){
897  if(j==1) fJetInfoVar[4]=-5;
898  else if(j==2) fJetInfoVar[6]=-5;
899  }
900  else {
901  if(j==1) fJetInfoVar[5]=-5;
902  else if(j==2) fJetInfoVar[7]=-5;
903  }
904  continue;
905  }
906  else{
907  fastjet::contrib::Nsubjettiness nSub(j,fastjet::contrib::KT_Axes(),fastjet::contrib::NormalizedMeasure(NSubBeta,R0));
908  NSubjettinessResult[j] = nSub.result(fOutputJets[0]);
909  if(j==2){ //find subjet axis to calculate delR
910  Subjet_Axes = nSub.currentAxes();
911  SubJet1_Axis = Subjet_Axes[0];
912  SubJet2_Axis = Subjet_Axes[1];
913 
914  Double_t SubJet1_Eta=SubJet1_Axis.pseudorapidity();
915  Double_t SubJet2_Eta=SubJet2_Axis.pseudorapidity();
916  Double_t SubJet1_Phi=SubJet1_Axis.phi();
917  if(SubJet1_Phi < -1*TMath::Pi()) SubJet1_Phi += (2*TMath::Pi());
918  else if (SubJet1_Phi > TMath::Pi()) SubJet1_Phi -= (2*TMath::Pi());
919  Double_t SubJet2_Phi=SubJet2_Axis.phi();
920  if(SubJet1_Phi < -1*TMath::Pi()) SubJet1_Phi += (2*TMath::Pi());
921  else if (SubJet1_Phi > TMath::Pi()) SubJet1_Phi -= (2*TMath::Pi());
922 
923  Double_t DeltaPhiSubJets,DeltaEtaSubJets;
924  DeltaPhiSubJets=SubJet1_Phi-SubJet2_Phi;
925  DeltaEtaSubJets=SubJet1_Eta-SubJet2_Eta;
926  if(DeltaPhiSubJets < -1*TMath::Pi()) DeltaPhiSubJets += (2*TMath::Pi());
927  else if (DeltaPhiSubJets > TMath::Pi()) DeltaPhiSubJets -= (2*TMath::Pi());
928 
929  DelR = TMath::Power(TMath::Power(DeltaPhiSubJets,2)+TMath::Power(DeltaEtaSubJets,2),0.5);
930  }
931  }
932  }
933  if(!fTruthJet){
934  fJetInfoVar[4]=NSubjettinessResult[1];
935  fJetInfoVar[6]=NSubjettinessResult[2];
936  // fJetInfoVar[12]=DelR;
937  }
938  else {
939  fJetInfoVar[5]=NSubjettinessResult[1];
940  fJetInfoVar[7]=NSubjettinessResult[2];
941  // fJetInfoVar[13]=DelR;
942  }
943 
944 
945 
946  fastjet::contrib::SoftDrop softdrop(beta, zcut);
947  //fastjet::contrib::SoftDrop softdrop_antikt(beta,zcut);
948  softdrop.set_verbose_structure(kTRUE);
949  //fastjet::JetDefinition jet_def_akt(fastjet::antikt_algorithm, 0.4);
950  // fastjet::contrib::Recluster *antiKT_Recluster(jet_def_akt);
951  fastjet::contrib::Recluster *recluster;
952  if(fReclusterAlgo == 2) recluster = new fastjet::contrib::Recluster(fastjet::kt_algorithm,1,true);
953  if(fReclusterAlgo == 1) recluster = new fastjet::contrib::Recluster(fastjet::antikt_algorithm,1,true);
954  if(fReclusterAlgo == 0) recluster = new fastjet::contrib::Recluster(fastjet::cambridge_algorithm,1,true);
955  softdrop.set_reclustering(true,recluster);
956  fastjet::PseudoJet finaljet = softdrop(fOutputJets[0]);
957  // fastjet::PseudoJet finaljet_antikt = softdrop_antikt(fOutputJets[0]);
958  //cout<< finaljet.structure_of<fastjet::contrib::SoftDrop>().symmetry()<<endl;
959  //cout<< finaljet_antikt.structure_of<fastjet::contrib::SoftDrop>().symmetry()<<endl;
960 
961 
962  AliEmcalJet* jet = new AliEmcalJet(finaljet.perp(), finaljet.eta(), finaljet.phi(), finaljet.m());
963  std::vector<fastjet::PseudoJet> fSDTracks=finaljet.constituents();
964  Double_t FastjetTrackDelR,EmcalTrackDelR;
965  for(Int_t i=0;i<fJet->GetNumberOfConstituents();i++){
966  if(i<=finaljet.constituents().size()){
967  FastjetTrackDelR = TMath::Sqrt(TMath::Power(fSDTracks[i].eta()-JetEta,2)+TMath::Power(fSDTracks[i].phi()-JetPhi,2));
968  FJTrackEta[i]=fSDTracks[i].eta();
969  FJTrackPhi[i]=fSDTracks[i].phi();
970  FJTrackPt[i]=fSDTracks[i].perp();
971  FJNTracks++;
972  }
973  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
974  EmcalTrackDelR = TMath::Sqrt(TMath::Power(fTrk->Eta()-JetEta,2)+TMath::Power(fTrk->Phi()-JetPhi,2));
975  }
976  Int_t NDroppedTracks = fJet->GetNumberOfTracks()-finaljet.constituents().size();
977  Int_t nConstituents(fClustSeqSA->constituents(finaljet).size());
978  jet->SetNumberOfTracks(nConstituents);
979  Double_t SymParam, Mu, DeltaR;
980  SymParam=(finaljet.structure_of<fastjet::contrib::SoftDrop>().symmetry());
981  Mu=(finaljet.structure_of<fastjet::contrib::SoftDrop>().mu());
982  DeltaR=(finaljet.structure_of<fastjet::contrib::SoftDrop>().delta_R());
983  //fhGroomedPtvJetPt->Fill(finaljet.perp(),fJet->Pt());
984  //fhDroppedBranches->Fill(finaljet.structure_of<fastjet::contrib::SoftDrop>().dropped_count());
985  if(!fTruthJet) fJetInfoVar[2]=SymParam;
986  else fJetInfoVar[3]=SymParam;
987  if(!fTruthJet) fJetInfoVar[12] = DeltaR;
988  else fJetInfoVar[13] = DeltaR;
989  if(!fTruthJet) fJetInfoVar[14]=finaljet.structure_of<fastjet::contrib::SoftDrop>().dropped_count();
990  else fJetInfoVar[15]=finaljet.structure_of<fastjet::contrib::SoftDrop>().dropped_count();
991  if(!(fJetShapeSub==kConstSub)){
992  if(!fTruthJet) fJetInfoVar[16]=(finaljet.perp()-(GetRhoVal(0)*fJet->Area()));
993  else fJetInfoVar[17]=finaljet.perp();
994  }
995  else{
996  if(!fTruthJet) fJetInfoVar[16]=(finaljet.perp());
997  else fJetInfoVar[17]=(finaljet.perp());
998  }
999  if(!fTruthJet) fJetInfoVar[18]=(finaljet.m());
1000  else fJetInfoVar[19]=(finaljet.m());
1001 
1002  return;
1003 
1004 
1005 }
1006 
1007 //--------------------------------------------------------------------------
1009 
1010  AliTrackContainer *PartCont = NULL;
1011  AliParticleContainer *PartContMC = NULL;
1012 
1013 
1014  if (fJetShapeSub==kConstSub){
1016  else PartCont = GetTrackContainer(1);
1017  }
1018  else{
1020  else PartCont = GetTrackContainer(0);
1021  }
1022 
1023  TClonesArray *TracksArray = NULL;
1024  TClonesArray *TracksArrayMC = NULL;
1025  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kGenOnTheFly) TracksArrayMC = PartContMC->GetArray();
1026  else TracksArray = PartCont->GetArray();
1027 
1029  if(!PartContMC || !TracksArrayMC) return -99999;
1030  }
1031  else {
1032  if(!PartCont || !TracksArray) return -99999;
1033  }
1034 
1035  AliAODTrack *Track = 0x0;
1036  Int_t Trigger_Index[100];
1037  for (Int_t i=0; i<100; i++) Trigger_Index[i] = 0;
1038  Int_t Trigger_Counter = 0;
1039  Int_t NTracks=0;
1040  if (fJetShapeType == AliAnalysisTaskRecoilJetYield::kGenOnTheFly) NTracks = TracksArrayMC->GetEntriesFast();
1041  else NTracks = TracksArray->GetEntriesFast();
1042  for(Int_t i=0; i < NTracks; i++){
1044  if((Track = static_cast<AliAODTrack*>(PartContMC->GetAcceptParticle(i)))){
1045  if (!Track) continue;
1046  if(TMath::Abs(Track->Eta())>0.9) continue;
1047  if (Track->Pt()<0.15) continue;
1048  if ((Track->Pt() >= PtMin) && (Track->Pt()< PtMax)) {
1049  Trigger_Index[Trigger_Counter] = i;
1050  Trigger_Counter++;
1051  }
1052  }
1053  }
1054  else{
1055  if((Track = static_cast<AliAODTrack*>(PartCont->GetAcceptTrack(i)))){
1056  if (!Track) continue;
1057  if(TMath::Abs(Track->Eta())>0.9) continue;
1058  if (Track->Pt()<0.15) continue;
1059  if ((Track->Pt() >= PtMin) && (Track->Pt()< PtMax)) {
1060  Trigger_Index[Trigger_Counter] = i;
1061  Trigger_Counter++;
1062  }
1063  }
1064  }
1065  }
1066  if (Trigger_Counter == 0) return -99999;
1067  Int_t RandomNumber = 0, Index = 0 ;
1068  TRandom3* Random = new TRandom3(0);
1069  RandomNumber = Random->Integer(Trigger_Counter);
1070  Index = Trigger_Index[RandomNumber];
1071  return Index;
1072 }
1073 
1074 
1075 //________________________________________________________________________
1077 {
1078  AliEmcalJet *jet2 = jet1->ClosestJet();
1079  if (!jet2) return -1;
1080 
1081  Double_t jetPt2 = jet2->Pt();
1082  if (jetPt2 <= 0) return -1;
1083 
1084  Int_t bgeom = kTRUE;
1085  if (!cont2) bgeom = kFALSE;
1086  Double_t sumPt = 0.;
1087  AliVParticle *vpf = 0x0;
1088  Int_t iFound = 0;
1089  for (Int_t icc = 0; icc < jet2->GetNumberOfTracks(); icc++) {
1090  Int_t idx = (Int_t)jet2->TrackAt(icc);
1091  //get particle
1092  AliVParticle *p2 = 0x0;
1093  if (bgeom) p2 = static_cast<AliVParticle*>(jet2->TrackAt(icc, cont2->GetArray()));
1094  iFound = 0;
1095  for (Int_t icf = 0; icf < jet1->GetNumberOfTracks(); icf++) {
1096  if (!bgeom && idx == jet1->TrackAt(icf) && iFound == 0 ) {
1097  iFound = 1;
1098  vpf = jet1->Track(icf);
1099  if (vpf) sumPt += vpf->Pt();
1100  continue;
1101  }
1102  if (bgeom){
1103  vpf = jet1->Track(icf);
1104  if (!vpf) continue;
1105  if (!SameParticle(vpf, p2, 1.e-4)) continue; //not the same particle
1106  sumPt += vpf->Pt();
1107  }
1108  }
1109  }
1110 
1111  Double_t fraction = sumPt / jetPt2;
1112 
1113  return fraction;
1114 }
1115 
1116 //________________________________________________________________________
1117 Bool_t AliAnalysisTaskRecoilJetYield::SameParticle(const AliVParticle* part1, const AliVParticle* part2, Double_t dist)
1118 {
1119  if(!part1) return kFALSE;
1120  if(!part2) return kFALSE;
1121  Double_t dPhi = TMath::Abs(part1->Phi() - part2->Phi());
1122  dPhi = TVector2::Phi_mpi_pi(dPhi);
1123  if (dPhi > dist) return kFALSE;
1124  return kTRUE;
1125 }
1126 
1127 //________________________________________________________________________
1129  //
1130  // retrieve event objects
1131  //
1133  return kFALSE;
1134 
1135  return kTRUE;
1136 }
1137 
1138 
1139 //_______________________________________________________________________
1141 {
1142  // Called once at the end of the analysis.
1143 
1144 
1145 }
Bool_t RetrieveEventObjects()
Retrieve common objects from event.
Double_t Area() const
Definition: AliEmcalJet.h:123
double Double_t
Definition: External.C:58
Double_t PTD(AliEmcalJet *Jet, Int_t JetContNb)
static Double_t DeltaPhi(Double_t phia, Double_t phib, Double_t rMin=-TMath::Pi()/2, Double_t rMax=3 *TMath::Pi()/2)
Definition: External.C:236
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:222
AliJetContainer * GetJetContainer(Int_t i=0) const
Bool_t FillHistograms()
Function filling histograms.
Double_t Eta() const
Definition: AliEmcalJet.h:114
Double_t Phi_Upper
Double_t Phi() const
Definition: AliEmcalJet.h:110
Container with name, TClonesArray and cuts for particles.
Double_t fEPV0
!event plane V0
Int_t GetLabel() const
Definition: AliEmcalJet.h:117
AliVParticle * Track(Int_t idx) const
virtual Bool_t AcceptJet(Int_t i, UInt_t &rejectionReason) const
Double_t Eta_Lower
Float_t Percentile_1
Double_t Eta_Upper
void SoftDrop(AliEmcalJet *fJet, AliJetContainer *fJetCont, double zcut, double beta, Bool_t fTruthJet)
Double_t Phi_Lower
UShort_t GetNumberOfConstituents() const
Definition: AliEmcalJet.h:133
Container for particles within the EMCAL framework.
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:153
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:132
TString kData
Declare data MC or deltaAOD.
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
AliParticleContainer * GetParticleContainer() const
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
Int_t SelectTriggerHadron(Float_t PtMin, Float_t PtMax)
Double_t GetFractionSharedPt_SubMatching(const AliEmcalJet *jet, AliParticleContainer *cont2=0x0) const
Int_t GetNJets() const
virtual AliVParticle * GetAcceptParticle(Int_t i=-1) const
Double_t fCent
!event centrality
const Int_t nBranch
Double_t Angularity(AliEmcalJet *Jet, Int_t JetContNb)
AliEmcalJet * GetNextAcceptJet()
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
Double_t RelativePhi(Double_t Phi1, Double_t Phi2)
static Bool_t SameParticle(const AliVParticle *part1, const AliVParticle *part2, Double_t dist=1.e-4)
Double_t Pt() const
Definition: AliEmcalJet.h:102
Double_t GetRhoVal(Int_t i=0) const
AliEmcalList * fOutput
!output list
AliTrackContainer * GetTrackContainer(Int_t i=0) const
Double_t DeltaR(const AliVParticle *part1, const AliVParticle *part2)
void SetMakeGeneralHistograms(Bool_t g)
virtual AliVTrack * GetAcceptTrack(Int_t i=-1) const
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
void SetNumberOfTracks(Int_t n)
Definition: AliEmcalJet.h:195
unsigned short UShort_t
Definition: External.C:28
const char Option_t
Definition: External.C:48
void UserCreateOutputObjects()
Main initialization function on the worker.
bool Bool_t
Definition: External.C:53
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.
Double_t RelativePhiEventPlane(Double_t EventPlane, Double_t Phi)
Double_t M() const
Definition: AliEmcalJet.h:113
Container for jet within the EMCAL jet framework.
AliEmcalJet * GetJet(Int_t i) const