AliPhysics  ec7afe5 (ec7afe5)
 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  fReclusterAlgo(0),
121  fSubMatching(kFALSE)
122 
123 {
124  for(Int_t i=0;i<nBranch;i++){
125  fJetInfoVar[i]=0;
126  }
127  SetMakeGeneralHistograms(kTRUE);
128  DefineOutput(1, TList::Class());
129  DefineOutput(2, TTree::Class());
130 }
131 
132 //________________________________________________________________________
134  AliAnalysisTaskEmcalJet(name, kTRUE),
135  fContainer(0),
136  fMinFractionShared(0),
137  fJetShapeType(kData),
138  fJetShapeSub(kNoSub),
139  fJetSelection(kInclusive),
140  fPtThreshold(-9999.),
141  fRMatching(0.2),
142  fPtMinTriggerHadron(20.),
143  fPtMaxTriggerHadron(50.),
144  fRecoilAngularWindow(0.6),
145  fSemigoodCorrect(0),
146  fHolePos(0),
147  fHoleWidth(0),
148  fCentSelectOn(kTRUE),
149  fCentMin(0),
150  fCentMax(10),
151  fSubJetAlgorithm(0),
152  fSubJetRadius(0.1),
153  fSubJetMinPt(1),
154  fJetRadius(0.4),
155  fRMatched(0.2),
156  fSharedFractionPtMin(0.5),
157  fDerivSubtrOrder(0),
158  fFullTree(kFALSE),
159  fBeta_SD(0),
160  fZCut(0.1),
161  fRho(1e-6),
162  fRhoParam(0),
163  fNsubMeasure(kFALSE),
164  fDoSoftDrop(kFALSE),
165  fhJetPt(0x0),
166  fhJetPhi(0x0),
167  fhJetEta(0x0),
168  fhJetMass(0x0),
169  fhJetRadius(0x0),
170  fhJetCounter(0x0),
171  fhNumberOfJetTracks(0x0),
172  fTreeJetInfo(0),
173  fhJetArea(0x0),
174  fhTrackPt(0x0),
175  fhGroomedPtvJetPt(0x0),
176  fhDroppedBranches(0x0),
177  fhPtTriggerHadron(0x0),
178  fh2PtTriggerHadronJet(0x0),
179  fhPhiTriggerHadronJet(0x0),
180  fhPhiTriggerHadronEventPlane(0x0),
181  fhPhiTriggerHadronEventPlaneTPC(0x0),
182  fhDetJetPt_Incl(0x0),
183  fhDetJetPt_Matched(0x0),
184  fReclusterAlgo(0),
185  fSubMatching(kFALSE)
186 
187 {
188  // Standard constructor.
189  for(Int_t i=0;i<nBranch;i++){
190  fJetInfoVar[i]=0;
191  }
193  DefineOutput(1, TList::Class());
194  DefineOutput(2, TTree::Class());
195 }
196 
197 //________________________________________________________________________
199 {
200  // Destructor.
201 }
202 
203 //________________________________________________________________________
205 {
206  // Create user output.
207 
209 
210  Bool_t oldStatus = TH1::AddDirectoryStatus();
211  TH1::AddDirectory(kFALSE);
212  TH1::AddDirectory(oldStatus);
213  const char* nameoutput = GetOutputSlot(2)->GetContainer()->GetName();
214  fTreeJetInfo = new TTree(nameoutput, nameoutput);
215 
216  if (!fFullTree){
217  const Int_t nVarMin = 18;
218  TString *fJetInfoVarNames = new TString [nVarMin];
219 
220  fJetInfoVarNames[0] = "Pt";
221  fJetInfoVarNames[1] = "Pt_Truth";
222  fJetInfoVarNames[2] = "SymParam";
223  fJetInfoVarNames[3] = "SymParam_Truth";
224  fJetInfoVarNames[4] = "Tau1";
225  fJetInfoVarNames[5] = "Tau1_Truth";
226  fJetInfoVarNames[6] = "Tau2";
227  fJetInfoVarNames[7] = "Tau2_Truth";
228  fJetInfoVarNames[8] = "PTD";
229  fJetInfoVarNames[9] = "PTD_Truth";
230  fJetInfoVarNames[10] = "Angularity";
231  fJetInfoVarNames[11] = "Angularity_Truth";
232  fJetInfoVarNames[12] = "DelR";
233  fJetInfoVarNames[13] = "DelR_Truth";
234  fJetInfoVarNames[14] = "N_Groomed_Branches";
235  fJetInfoVarNames[15] = "N_Groomed_Branches_Truth";
236  fJetInfoVarNames[16] = "Groomed_Jet_Pt";
237  fJetInfoVarNames[17] = "Groomed_Jet_Pt_Truth";
238 
239 
240 
241  for(Int_t ivar=0; ivar < nVarMin; ivar++){
242  cout<<"looping over variables"<<endl;
243  fTreeJetInfo->Branch(fJetInfoVarNames[ivar].Data(), &fJetInfoVar[ivar], Form("%s/D", fJetInfoVarNames[ivar].Data()));
244  }
245  }
246 
248 
249  fhJetPt= new TH1F("fhJetPt", "Jet Pt",150,-0.5,149.5 );
250  fOutput->Add(fhJetPt);
251  fhJetPhi= new TH1F("fhJetPhi", "Jet Phi",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
252  fOutput->Add(fhJetPhi);
253  fhJetEta= new TH1F("fhJetEta", "Jet Eta", Eta_Bins_1, Eta_Lower, Eta_Upper);
254  fOutput->Add(fhJetEta);
255  fhJetMass= new TH1F("fhJetMass", "Jet Mass", 4000,-0.5, 39.5);
256  fOutput->Add(fhJetMass);
257  fhJetRadius= new TH1F("fhJetRadius", "Jet Radius", 100, -0.05,0.995);
258  fOutput->Add(fhJetRadius);
259  fhNumberOfJetTracks= new TH1F("fhNumberOfJetTracks", "Number of Tracks within a Jet", 300, -0.5,299.5);
261  fhJetCounter= new TH1F("fhJetCounter", "Jet Counter", 150, -0.5, 149.5);
262  fOutput->Add(fhJetCounter);
263  fhJetArea= new TH1F("fhJetArea", "Jet Area", 400,-0.5, 1.95);
264  fOutput->Add(fhJetArea);
265  fhTrackPt= new TH1F("fhTrackPt", "Track Pt",600,-0.5,59.5);
266  fOutput->Add(fhTrackPt);
267  fhGroomedPtvJetPt= new TH2F("fhGroomedPtvJetPt","Groomed Jet p_{T} v Original Jet p_{T}",150,0,150,150,0,150);
269  fhDroppedBranches= new TH1F("fhDroppedBranches","Number of Softdropped branches",50,0,50);
271  fhDetJetPt_Incl= new TH1F("fhDetJetPt_Incl", "Jet Pt",200,-0.5,199.5 );
272  fOutput->Add(fhDetJetPt_Incl);
273  fhDetJetPt_Matched= new TH1F("fhDetJetPt_Matched", "Jet Pt",200,-0.5,199.5 );
275 
276  if(fJetSelection == kRecoil){
277  fhPtTriggerHadron= new TH1F("fhPtTriggerHadron", "fhPtTriggerHadron",1500,-0.5,149.5);
279  fh2PtTriggerHadronJet= new TH2F("fh2PtTriggerHadronJet", "fh2PtTriggerHadronJet",1500,-0.5,149.5,1500,-0.5,149.5);
281  fhPhiTriggerHadronJet= new TH1F("fhPhiTriggerHadronJet", "fhPhiTriggerHadronJet",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
283  fhPhiTriggerHadronEventPlane= new TH1F("fhPhiTriggerHadronEventPlane", "fhPhiTriggerHadronEventPlane",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
285  fhPhiTriggerHadronEventPlaneTPC= new TH1F("fhPhiTriggerHadronEventPlaneTPC", "fhPhiTriggerHadronEventPlaneTPC",360 , -1.5*(TMath::Pi()), 1.5*(TMath::Pi()));
287 
288 
289  }
290  }
291  PostData(1,fOutput);
292  PostData(2,fTreeJetInfo);
293  //cout<<"End of UserCreateOutputObjects"<<endl;
294  // delete [] fShapesVarNames;
295 }
296 
297 //________________________________________________________________________
299 {
300  // Run analysis code here, if needed. It will be executed before FillHistograms().
301  return kTRUE;
302 }
303 
304 //________________________________________________________________________
306 {
307  AliMultSelection *MultSelection = 0x0;
308  AliAODEvent *fEvent = dynamic_cast<AliAODEvent*>(InputEvent());
309  if (!fEvent) {
310  Error("UserExec","AOD not available");
311  return 0;
312  }
313  MultSelection = (AliMultSelection * ) fEvent->FindListObject("MultSelection");
314  if(!MultSelection) {
315  AliWarning("AliMultSelection object not found!");
316  }
317  else{
318  Percentile_1 = MultSelection->GetMultiplicityPercentile("V0M");
319  }
320  if (fCentSelectOn){
321  if ((fCent>fCentMax) || (fCent<fCentMin)) return 0;
322  }
323 
325  AliAODTrack *TriggerHadron = 0x0;
326  if (fJetSelection == kRecoil) {
327  //you have to set a flag and the limits of the pT interval for your trigger
329  if (TriggerHadronLabel==-99999) return 0; //Trigger Hadron Not Found
330  AliTrackContainer *PartCont =NULL;
331  if (fJetShapeSub==kConstSub) PartCont = GetTrackContainer(1);
332  else PartCont = GetTrackContainer(0);
333  TClonesArray *TrackArray = PartCont->GetArray();
334  TriggerHadron = static_cast<AliAODTrack*>(TrackArray->At(TriggerHadronLabel));
335  if (!TriggerHadron) return 0;//No trigger hadron with label found
336  if(fSemigoodCorrect){
337  Double_t HoleDistance=RelativePhi(TriggerHadron->Phi(),fHolePos);
338  if(TMath::Abs(HoleDistance)+fHoleWidth+fJetRadius>TMath::Pi()-fRecoilAngularWindow) return 0;
339  }
340  fhPtTriggerHadron->Fill(TriggerHadron->Pt()); //Needed for per trigger Normalisation
341  fhPhiTriggerHadronEventPlane->Fill(TMath::Abs(RelativePhiEventPlane(fEPV0,TriggerHadron->Phi()))); //fEPV0 is the event plane from AliAnalysisTaskEmcal
342  fhPhiTriggerHadronEventPlaneTPC->Fill(TMath::Abs(RelativePhiEventPlane(((AliVAODHeader*)(InputEvent()->GetHeader()))->GetEventplane(),TriggerHadron->Phi()))); //TPC event plane
343  }
344 
347  // cout<<"entering kData"<<endl;
348 
349  AliEmcalJet *Jet1 = NULL; //Original Jet in the event
350  AliJetContainer *JetCont= GetJetContainer(0); //Jet Container for event
351  Int_t JetCounter=0; //Counts number of jets in event
352  Double_t JetPhi=0;
353  Bool_t EventCounter=kFALSE;
354  Double_t JetPt_ForThreshold=0;
355  AliEmcalJet *Jet2= NULL;
356  if(JetCont) {
357  JetCont->ResetCurrentID();
358  while((Jet1=JetCont->GetNextAcceptJet())) {
359  if(!Jet1) continue;
360  if (fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) JetPt_ForThreshold = Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
361  else JetPt_ForThreshold = Jet1->Pt();
362  if(JetPt_ForThreshold<fPtThreshold) {
363  continue;
364  }
365  else {
366  Float_t RecoilDeltaPhi = 0.;
367  if (fJetSelection == kRecoil){
368  RecoilDeltaPhi = RelativePhi(TriggerHadron->Phi(), Jet1->Phi());
369  if (TMath::Abs(RecoilDeltaPhi) < (TMath::Pi() - fRecoilAngularWindow)) continue; //accept the jet only if it overlaps with the recoil phi area of the trigger
370  fh2PtTriggerHadronJet->Fill(TriggerHadron->Pt(), Jet1->Pt());
371  fhPhiTriggerHadronJet->Fill(RelativePhi(TriggerHadron->Phi(), Jet1->Phi()));
372  }
373  if (!EventCounter){
374  EventCounter=kTRUE;
375  }
376  JetCounter++;
377  fhJetPt->Fill(Jet1->Pt());
378  fhJetArea->Fill(Jet1->Area());
379  JetPhi=Jet1->Phi();
380  if(JetPhi < -1*TMath::Pi()) JetPhi += (2*TMath::Pi());
381  else if (JetPhi > TMath::Pi()) JetPhi -= (2*TMath::Pi());
382  fhJetPhi->Fill(JetPhi);
383  fhJetEta->Fill(Jet1->Eta());
384  fhJetMass->Fill(Jet1->M());
385  fhJetRadius->Fill(TMath::Sqrt((Jet1->Area()/TMath::Pi()))); //Radius of Jets per event
386  fhNumberOfJetTracks->Fill(Jet1->GetNumberOfTracks());
387  if(fJetShapeSub==kNoSub || fJetShapeSub==kDerivSub) fJetInfoVar[0]= Jet1->Pt()-(GetRhoVal(0)*Jet1->Area());
388  else fJetInfoVar[0]=Jet1->Pt();
389  fJetInfoVar[1]=0;
390  if(fDoSoftDrop) {
391  SoftDrop(Jet1,JetCont,fZCut,fBeta_SD,kFALSE);
392  SoftDrop(Jet1,JetCont,fZCut,fBeta_SD,kTRUE);
393  }
394  else{
395  fJetInfoVar[2]=0;
396  fJetInfoVar[3]=0;
397  fJetInfoVar[4]=0;
398  fJetInfoVar[5]=0;
399  fJetInfoVar[6]=0;
400  fJetInfoVar[7]=0;
401  fJetInfoVar[12]=0;
402  fJetInfoVar[13]=0;
403  fJetInfoVar[14]=0;
404  fJetInfoVar[15]=0;
405  fJetInfoVar[16]=0;
406  fJetInfoVar[17]=0;
407  }
408  fJetInfoVar[8]=PTD(Jet1,0);
409  fJetInfoVar[9]=0;
410  fJetInfoVar[10]=Angularity(Jet1,0);
411  fJetInfoVar[11]=0;
412  fTreeJetInfo->Fill();
413  }
414  }
415  fhJetCounter->Fill(JetCounter); //Number of Jets in Each Event
416  }
417  }
418 
420  if(fJetShapeType == kDetEmbPart){
421  AliEmcalJet *JetHybridS = NULL; //Subtracted hybrid Jet
422  AliEmcalJet *JetHybridUS = NULL; //Unsubtracted Hybrid Jet //For matching SubtractedHybrid->DetPythia this jet container is also Subtracted Hybrid
423  AliEmcalJet *JetPythDet = NULL; //Detector Level Pythia Jet
424  AliEmcalJet *JetPythTrue = NULL; //Particle Level Pyhtia Jet
425  AliJetContainer *JetContHybridS= GetJetContainer(0); //Jet Container for Subtracted Hybrid Jets
426  AliJetContainer *JetContHybridUS= GetJetContainer(1); //Jet Container for Unsubtracted Hybrid Jets
427  AliJetContainer *JetContPythDet= GetJetContainer(2); //Jet Container for Detector Level Pyhtia Jets
428  AliJetContainer *JetContPythTrue= GetJetContainer(3); //Jet Container for Particle Level Pythia Jets
429 
430 
431 
432  Bool_t JetsMatched = kFALSE;
433  Double_t JetPtThreshold;
434  JetContHybridS->ResetCurrentID();
435  JetContHybridUS->ResetCurrentID();
436  JetContPythDet->ResetCurrentID();
437  JetContPythTrue->ResetCurrentID();
438 
439  while((JetPythDet = JetContPythDet->GetNextAcceptJet())){
440  fhDetJetPt_Incl->Fill(JetPythDet->Pt()); //Fill histogram with all Detector level jets
441  }
442 
443 
444  while((JetHybridS = JetContHybridS->GetNextAcceptJet())){
445  if (fJetShapeSub==kConstSub) JetPtThreshold=JetHybridS->Pt();
446  else JetPtThreshold=JetHybridS->Pt()-(GetRhoVal(0)*JetHybridS->Area());
447  if ( (!JetHybridS) || (JetPtThreshold<fPtThreshold)) continue;
448  if (fJetShapeSub==kConstSub){
449  Int_t JetNumber=-1;
450  for(Int_t i = 0; i<JetContHybridUS->GetNJets(); i++) {
451  JetHybridUS = JetContHybridUS->GetJet(i);
452  if (!JetHybridUS) continue;
453 
454  if(JetHybridUS->GetLabel()==JetHybridS->GetLabel()) {
455  JetNumber=i;
456  }
457  }
458  if(JetNumber==-1) continue;
459  JetHybridUS=JetContHybridUS->GetJet(JetNumber);
460  //if(JetHybridUS) cout<<"Matched to jet i = "<< JetNumber<<endl;
461  if (JetContHybridUS->AliJetContainer::GetFractionSharedPt(JetHybridUS)<fSharedFractionPtMin && !fSubMatching) {
462  //cout<<"Fraction shared pt below cut = "<<JetContHybridUS->AliJetContainer::GetFractionSharedPt(JetHybridUS)<<endl;
463  continue;
464  }
466  continue;
467  }
468  JetPythDet=JetHybridUS->ClosestJet();
469  }
470  if(!(fJetShapeSub==kConstSub)){
471  if (JetContHybridS->AliJetContainer::GetFractionSharedPt(JetHybridS)<fSharedFractionPtMin) continue;
472  JetPythDet = JetHybridS->ClosestJet();
473  }
474  if (!JetHybridUS) {
475  Printf("Unsubtracted embedded jet does not exist, returning");
476  continue;
477  }
478  if (!JetPythDet) continue;
479  fhDetJetPt_Matched->Fill(JetPythDet->Pt()); //Fill only matched detector level jets for tagging efficiency comparison
480  JetPythTrue=JetPythDet->ClosestJet();
481  if(!JetPythTrue) continue;
482  JetsMatched=kTRUE;
483 
484  if(fJetShapeSub==kConstSub) fJetInfoVar[0]=JetHybridS->Pt();
485  else fJetInfoVar[0]=JetHybridS->Pt()-(GetRhoVal(0)*JetHybridS->Area());
486  fJetInfoVar[1]=JetPythTrue->Pt();
487  if(fDoSoftDrop) {
488  SoftDrop(JetHybridS,JetContHybridS,fZCut,fBeta_SD,kFALSE);
489  SoftDrop(JetPythTrue,JetContPythTrue,fZCut,fBeta_SD,kTRUE);
490  }
491  else{
492  fJetInfoVar[2]=0;
493  fJetInfoVar[3]=0;
494  fJetInfoVar[4]=0;
495  fJetInfoVar[5]=0;
496  fJetInfoVar[6]=0;
497  fJetInfoVar[7]=0;
498  fJetInfoVar[12]=0;
499  fJetInfoVar[13]=0;
500  fJetInfoVar[14]=0;
501  fJetInfoVar[15]=0;
502  fJetInfoVar[16]=0;
503  fJetInfoVar[17]=0;
504  }
505  fJetInfoVar[8]=PTD(JetHybridS,0);
506  fJetInfoVar[9]=PTD(JetPythTrue,0);
507  fJetInfoVar[10]=Angularity(JetHybridS,0);
508  fJetInfoVar[11]=Angularity(JetPythTrue,0);
509  fTreeJetInfo->Fill();
510  }
511 
512  }
513 
514 
515 
516  return kTRUE;
517 }
518 
519 //________________________________________________________________________
521 
522  if(Phi < -1*TMath::Pi()) Phi += (2*TMath::Pi());
523  else if (Phi > TMath::Pi()) Phi -= (2*TMath::Pi());
524  Double_t DeltaPhi=Phi-EventPlane;
525  if(DeltaPhi < -1*TMath::Pi()) DeltaPhi += (2*TMath::Pi());
526  else if (DeltaPhi > TMath::Pi()) DeltaPhi -= (2*TMath::Pi());
527  return DeltaPhi;
528 }
529 //________________________________________________________________________
531 
532  if(Phi1 < -1*TMath::Pi()) Phi1 += (2*TMath::Pi());
533  else if (Phi1 > TMath::Pi()) Phi1 -= (2*TMath::Pi());
534  if(Phi2 < -1*TMath::Pi()) Phi2 += (2*TMath::Pi());
535  else if (Phi2 > TMath::Pi()) Phi2 -= (2*TMath::Pi());
536  Double_t DeltaPhi=Phi2-Phi1;
537  if(DeltaPhi < -1*TMath::Pi()) DeltaPhi += (2*TMath::Pi());
538  else if (DeltaPhi > TMath::Pi()) DeltaPhi -= (2*TMath::Pi());
539  return DeltaPhi;
540 }
541 
542 
543 
544 
545 //--------------------------------------------------------------------------
547 
548  AliJetContainer *JetCont = GetJetContainer(JetContNb);
549  Double_t Angularity_Numerator=0;
550  Double_t Angularity_Denominator=0;
551  AliVParticle *Particle=0x0;
552  Double_t DeltaPhi=0;
553 
554 
555  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks
556  Particle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
557  if(!Particle) continue;
558  DeltaPhi=RelativePhi(Jet->Phi(),Particle->Phi());
559  Angularity_Numerator=Angularity_Numerator+(Particle->Pt()*TMath::Sqrt(((Particle->Eta()-Jet->Eta())*(Particle->Eta()-Jet->Eta()))+(DeltaPhi*DeltaPhi)));
560  Angularity_Denominator= Angularity_Denominator+Particle->Pt();
561  }
562  if(Angularity_Denominator!=0) return Angularity_Numerator/Angularity_Denominator;
563  else return -1;
564 
565 }
566 
567 
568 
569 //--------------------------------------------------------------------------
571 
572  AliJetContainer *JetCont = GetJetContainer(JetContNb);
573  Double_t PTD_Numerator=0; //Reset these values
574  Double_t PTD_Denominator=0;
575  AliVParticle *Particle=0x0;
576  Double_t DeltaPhi=0;
577  for (Int_t i=0; i< Jet->GetNumberOfTracks(); i++){ //loops through all tracks
578  Particle = static_cast<AliVParticle*>(Jet->TrackAt(i, JetCont->GetParticleContainer()->GetArray()));
579  if(!Particle) continue;
580  PTD_Numerator=PTD_Numerator+(Particle->Pt()*Particle->Pt());
581  PTD_Denominator=PTD_Denominator+Particle->Pt();
582  }
583  if(PTD_Denominator!=0) return TMath::Sqrt(PTD_Numerator)/PTD_Denominator;
584  else return -1;
585 
586 }
587 
588 //_________________________________________________________________________
589  void AliAnalysisTaskRecoilJetYield::SoftDrop(AliEmcalJet *fJet,AliJetContainer *fJetCont, double zcut, double beta, Bool_t fTruthJet){
590  std::vector<fastjet::PseudoJet> fInputVectors;
591  Double_t JetInvMass=0, PseudJetInvMass=0, TrackMom = 0, TrackEnergy = 0;
592  AliParticleContainer *fTrackCont = fJetCont->GetParticleContainer();
593  Double_t JetEta=fJet->Eta(),JetPhi=fJet->Phi();
594  Double_t FJTrackEta[9999],FJTrackPhi[9999],FJTrackPt[9999],EmcalJetTrackEta[9999],EmcalJetTrackPhi[9999],EmcalJetTrackPt[9999];
595  UShort_t FJNTracks=0,EmcalJetNTracks=0;
596  if (fTrackCont) for (Int_t i=0; i<fJet->GetNumberOfTracks(); i++) {
597  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
598  JetInvMass += fTrk->M();
599  if (!fTrk) continue;
600  fastjet::PseudoJet PseudoTracks(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
601  TrackMom += TMath::Sqrt(TMath::Power(fTrk->Px(),2)+TMath::Power(fTrk->Py(),2)+TMath::Power(fTrk->Pz(),2));
602  TrackEnergy += fTrk->E();
603  PseudoTracks.set_user_index(fJet->TrackAt(i)+100);
604  PseudJetInvMass += PseudoTracks.m();
605  fInputVectors.push_back(PseudoTracks);
606  EmcalJetTrackEta[i]=fTrk->Eta();
607  EmcalJetTrackPhi[i]=fTrk->Phi();
608  EmcalJetTrackPt[i]=fTrk->Pt();
609  EmcalJetNTracks++;
610  }
611  fastjet::JetDefinition *fJetDef;
612  fastjet::ClusterSequence *fClustSeqSA;
613 
614 
615 
616  fJetDef = new fastjet::JetDefinition(fastjet::antikt_algorithm, fJetRadius*2, static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
617 
618  try {
619  fClustSeqSA = new fastjet::ClusterSequence(fInputVectors, *fJetDef);
620  } catch (fastjet::Error) {
621  AliError(" [w] FJ Exception caught.");
622  //return -1;
623  }
624 
625  std::vector<fastjet::PseudoJet> fOutputJets;
626  fOutputJets.clear();
627  fOutputJets=fClustSeqSA->inclusive_jets(0);
628 
629 
630  std::vector<fastjet::PseudoJet> jet_constituents = fOutputJets[0].constituents();
631  Float_t NSubjettinessResult[3], NSubBeta = 1, R0=0.4;
632  std::vector<fastjet::PseudoJet> Subjet_Axes;
633  fastjet::PseudoJet SubJet1_Axis,SubJet2_Axis;
634  Double_t DelR=-5;
635 
636 
637  for(Int_t j=1; j<3; j++){
638  if(jet_constituents.size() < j){
639  if(!fTruthJet){
640  if(j==1) fJetInfoVar[4]=-5;
641  else if(j==2) fJetInfoVar[6]=-5;
642  }
643  else {
644  if(j==1) fJetInfoVar[5]=-5;
645  else if(j==2) fJetInfoVar[7]=-5;
646  }
647  continue;
648  }
649  else{
650  fastjet::contrib::Nsubjettiness nSub(j,fastjet::contrib::KT_Axes(),fastjet::contrib::NormalizedMeasure(NSubBeta,R0));
651  NSubjettinessResult[j] = nSub.result(fOutputJets[0]);
652  if(j==2){ //find subjet axis to calculate delR
653  Subjet_Axes = nSub.currentAxes();
654  SubJet1_Axis = Subjet_Axes[0];
655  SubJet2_Axis = Subjet_Axes[1];
656 
657  Double_t SubJet1_Eta=SubJet1_Axis.pseudorapidity();
658  Double_t SubJet2_Eta=SubJet2_Axis.pseudorapidity();
659  Double_t SubJet1_Phi=SubJet1_Axis.phi();
660  if(SubJet1_Phi < -1*TMath::Pi()) SubJet1_Phi += (2*TMath::Pi());
661  else if (SubJet1_Phi > TMath::Pi()) SubJet1_Phi -= (2*TMath::Pi());
662  Double_t SubJet2_Phi=SubJet2_Axis.phi();
663  if(SubJet1_Phi < -1*TMath::Pi()) SubJet1_Phi += (2*TMath::Pi());
664  else if (SubJet1_Phi > TMath::Pi()) SubJet1_Phi -= (2*TMath::Pi());
665 
666  Double_t DeltaPhiSubJets,DeltaEtaSubJets;
667  DeltaPhiSubJets=SubJet1_Phi-SubJet2_Phi;
668  DeltaEtaSubJets=SubJet1_Eta-SubJet2_Eta;
669  if(DeltaPhiSubJets < -1*TMath::Pi()) DeltaPhiSubJets += (2*TMath::Pi());
670  else if (DeltaPhiSubJets > TMath::Pi()) DeltaPhiSubJets -= (2*TMath::Pi());
671 
672  DelR = TMath::Power(TMath::Power(DeltaPhiSubJets,2)+TMath::Power(DeltaEtaSubJets,2),0.5);
673  }
674  }
675  }
676  if(!fTruthJet){
677  fJetInfoVar[4]=NSubjettinessResult[1];
678  fJetInfoVar[6]=NSubjettinessResult[2];
679  // fJetInfoVar[12]=DelR;
680  }
681  else {
682  fJetInfoVar[5]=NSubjettinessResult[1];
683  fJetInfoVar[7]=NSubjettinessResult[2];
684  // fJetInfoVar[13]=DelR;
685  }
686 
687 
688 
689  fastjet::contrib::SoftDrop softdrop(beta, zcut);
690  //fastjet::contrib::SoftDrop softdrop_antikt(beta,zcut);
691  softdrop.set_verbose_structure(kTRUE);
692  //fastjet::JetDefinition jet_def_akt(fastjet::antikt_algorithm, 0.4);
693  // fastjet::contrib::Recluster *antiKT_Recluster(jet_def_akt);
694  fastjet::contrib::Recluster *recluster;
695  if(fReclusterAlgo == 2) recluster = new fastjet::contrib::Recluster(fastjet::kt_algorithm,1,true);
696  if(fReclusterAlgo == 1) recluster = new fastjet::contrib::Recluster(fastjet::antikt_algorithm,1,true);
697  if(fReclusterAlgo == 0) recluster = new fastjet::contrib::Recluster(fastjet::cambridge_algorithm,1,true);
698  softdrop.set_reclustering(true,recluster);
699  fastjet::PseudoJet finaljet = softdrop(fOutputJets[0]);
700  // fastjet::PseudoJet finaljet_antikt = softdrop_antikt(fOutputJets[0]);
701  //cout<< finaljet.structure_of<fastjet::contrib::SoftDrop>().symmetry()<<endl;
702  //cout<< finaljet_antikt.structure_of<fastjet::contrib::SoftDrop>().symmetry()<<endl;
703 
704 
705  AliEmcalJet* jet = new AliEmcalJet(finaljet.perp(), finaljet.eta(), finaljet.phi(), finaljet.m());
706  std::vector<fastjet::PseudoJet> fSDTracks=finaljet.constituents();
707  Double_t FastjetTrackDelR,EmcalTrackDelR;
708  for(Int_t i=0;i<fJet->GetNumberOfConstituents();i++){
709  if(i<=finaljet.constituents().size()){
710  FastjetTrackDelR = TMath::Sqrt(TMath::Power(fSDTracks[i].eta()-JetEta,2)+TMath::Power(fSDTracks[i].phi()-JetPhi,2));
711  FJTrackEta[i]=fSDTracks[i].eta();
712  FJTrackPhi[i]=fSDTracks[i].phi();
713  FJTrackPt[i]=fSDTracks[i].perp();
714  FJNTracks++;
715  }
716  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
717  EmcalTrackDelR = TMath::Sqrt(TMath::Power(fTrk->Eta()-JetEta,2)+TMath::Power(fTrk->Phi()-JetPhi,2));
718  }
719  Int_t NDroppedTracks = fJet->GetNumberOfTracks()-finaljet.constituents().size();
720  Int_t nConstituents(fClustSeqSA->constituents(finaljet).size());
721  jet->SetNumberOfTracks(nConstituents);
722  Double_t SymParam, Mu, DeltaR;
723  SymParam=(finaljet.structure_of<fastjet::contrib::SoftDrop>().symmetry());
724  Mu=(finaljet.structure_of<fastjet::contrib::SoftDrop>().mu());
725  DeltaR=(finaljet.structure_of<fastjet::contrib::SoftDrop>().delta_R());
726  fhGroomedPtvJetPt->Fill(finaljet.perp(),fJet->Pt());
727  fhDroppedBranches->Fill(finaljet.structure_of<fastjet::contrib::SoftDrop>().dropped_count());
728  if(!fTruthJet) fJetInfoVar[2]=SymParam;
729  else fJetInfoVar[3]=SymParam;
730  if(!fTruthJet) fJetInfoVar[12] = DeltaR;
731  else fJetInfoVar[13] = DeltaR;
732  if(!fTruthJet) fJetInfoVar[14]=finaljet.structure_of<fastjet::contrib::SoftDrop>().dropped_count();
733  else fJetInfoVar[15]=finaljet.structure_of<fastjet::contrib::SoftDrop>().dropped_count();
734  if(!fTruthJet) fJetInfoVar[16]=finaljet.perp();
735  else fJetInfoVar[17]=finaljet.perp();
736 
737 
738  return;
739 
740 
741 }
742 
743 //________________________________________________________________________
745  AliTrackContainer *PartCont = NULL;
746  if (fJetShapeSub==kConstSub) PartCont = GetTrackContainer(1);
747  else PartCont = GetTrackContainer(0);
748  TClonesArray *TracksArray = PartCont->GetArray();
749  if(!PartCont || !TracksArray) return -99999;
750  AliAODTrack *Track = 0x0;
751  Int_t Trigger_Index[10000];
752  for (Int_t i=0; i<100; i++) Trigger_Index[i] = 0;
753  Int_t Trigger_Counter = 0;
754  for(Int_t i=0; i < TracksArray->GetEntriesFast(); i++){
755  if((Track = static_cast<AliAODTrack*>(PartCont->GetAcceptTrack(i)))){
756  if (!Track) continue;
757  fhTrackPt->Fill(Track->Pt());
758  if(TMath::Abs(Track->Eta())>0.9) continue;
759  if (Track->Pt()<0.15) continue;
760  if ((Track->Pt() >= PtMin) && (Track->Pt()< PtMax)) {
761  Trigger_Index[Trigger_Counter] = i;
762  Trigger_Counter++;
763  }
764  }
765  }
766  if (Trigger_Counter == 0) return -99999;
767  Int_t RandomNumber = 0, Index = 0 ;
768  TRandom3* Random = new TRandom3(0);
769  RandomNumber = Random->Integer(Trigger_Counter);
770  Index = Trigger_Index[RandomNumber];
771  return Index;
772 }
773 
774 
775 //________________________________________________________________________
777 {
778  AliEmcalJet *jet2 = jet1->ClosestJet();
779  if (!jet2) return -1;
780 
781  Double_t jetPt2 = jet2->Pt();
782  if (jetPt2 <= 0) return -1;
783 
784  Int_t bgeom = kTRUE;
785  if (!cont2) bgeom = kFALSE;
786  Double_t sumPt = 0.;
787  AliVParticle *vpf = 0x0;
788  Int_t iFound = 0;
789  for (Int_t icc = 0; icc < jet2->GetNumberOfTracks(); icc++) {
790  Int_t idx = (Int_t)jet2->TrackAt(icc);
791  //get particle
792  AliVParticle *p2 = 0x0;
793  if (bgeom) p2 = static_cast<AliVParticle*>(jet2->TrackAt(icc, cont2->GetArray()));
794  iFound = 0;
795  for (Int_t icf = 0; icf < jet1->GetNumberOfTracks(); icf++) {
796  if (!bgeom && idx == jet1->TrackAt(icf) && iFound == 0 ) {
797  iFound = 1;
798  vpf = jet1->Track(icf);
799  if (vpf) sumPt += vpf->Pt();
800  continue;
801  }
802  if (bgeom){
803  vpf = jet1->Track(icf);
804  if (!vpf) continue;
805  if (!SameParticle(vpf, p2, 1.e-4)) continue; //not the same particle
806  sumPt += vpf->Pt();
807  }
808  }
809  }
810 
811  Double_t fraction = sumPt / jetPt2;
812 
813  return fraction;
814 }
815 
816 //________________________________________________________________________
817 Bool_t AliAnalysisTaskRecoilJetYield::SameParticle(const AliVParticle* part1, const AliVParticle* part2, Double_t dist)
818 {
819  if(!part1) return kFALSE;
820  if(!part2) return kFALSE;
821  Double_t dPhi = TMath::Abs(part1->Phi() - part2->Phi());
822  dPhi = TVector2::Phi_mpi_pi(dPhi);
823  if (dPhi > dist) return kFALSE;
824  return kTRUE;
825 }
826 
827 //________________________________________________________________________
829  //
830  // retrieve event objects
831  //
833  return kFALSE;
834 
835  return kTRUE;
836 }
837 
838 
839 //_______________________________________________________________________
841 {
842  // Called once at the end of the analysis.
843 
844 
845 }
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:221
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
ClassImp(AliAnalysisTaskRecoilJetYield) AliAnalysisTaskRecoilJetYield
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
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() const
int Int_t
Definition: External.C:63
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
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:194
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