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