AliPhysics  d497afb (d497afb)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskJetPP.cxx
Go to the documentation of this file.
1 #ifndef ALIANALYSISTASKSE_H
2 
3 #include <Riostream.h>
4 #include <TROOT.h>
5 #include <TFile.h>
6 #include <TChain.h>
7 #include <TTree.h>
8 #include <TKey.h>
9 #include <TProfile.h>
10 #include <TProfile2D.h>
11 #include <TH1.h>
12 #include <TH1F.h>
13 #include <TH2F.h>
14 #include <TH1D.h>
15 #include <TH1I.h>
16 #include <TArrayF.h>
17 #include <TArrayD.h>
18 #include <THnSparse.h>
19 #include <TCanvas.h>
20 #include <TList.h>
21 #include <TClonesArray.h>
22 #include <TObject.h>
23 #include <TMath.h>
24 #include <TSystem.h>
25 #include <TInterpreter.h>
26 #include "AliAnalysisTask.h"
27 #include "AliCentrality.h"
28 #include "AliStack.h"
29 #include "AliESDEvent.h"
30 #include "AliEmcalList.h"
31 #include "AliESDInputHandler.h"
32 #include "AliAODEvent.h"
33 #include "AliAODHandler.h"
34 #include "AliAnalysisManager.h"
35 #include "AliAnalysisTaskSE.h"
37 #include "AliParticleContainer.h"
38 #include "AliInputEventHandler.h"
39 #include "math.h"
40 #endif
41 
42 #include <time.h>
43 #include <TRandom3.h>
44 #include "AliGenEventHeader.h"
45 #include "AliGenPythiaEventHeader.h"
46 #include "AliGenHijingEventHeader.h"
47 #include "AliAODMCHeader.h"
48 #include "AliMCEvent.h"
49 #include "AliLog.h"
50 #include <AliEmcalJet.h>
51 #include <AliPicoTrack.h>
52 #include "AliVEventHandler.h"
53 #include "AliVParticle.h"
54 #include "AliAODMCParticle.h"
55 #include "AliAODTrack.h"
56 #include "AliAnalysisUtils.h"
57 #include "AliRhoParameter.h"
58 #include "TVector3.h"
59 #include "AliVVertex.h"
60 
61 #include <stdio.h>
62 #include <stdlib.h>
63 #include <vector>
64 
65 #include "AliJetContainer.h"
66 #include "AliAnalysisTaskEmcal.h"
68 #include "AliAnalysisTaskJetPP.h"
69 #include "AliHeader.h"
70 #include "AliRunLoader.h"
71 #include "AliVVZERO.h"
72 #include "AliVZDC.h"
73 #include "AliExternalTrackParam.h"
74 using namespace std;
75 
76 // Analyis of jets in pp collisions
77 // Author Peter Pribeli (12.Sep. 2017)
78 
80 ClassImp(AliAnalysisTaskJetPP)
82 
86  fUseDefaultVertexCut(1), fUsePileUpCut(1), fIsMC(0),
87  fSignalJetRadius(0.4),
88 fSignalJetEtaWindow(0.9 - fSignalJetRadius), fTrackEtaWindow(0.9), fMinTrackPt(0.150), fMinJetArea(0.0),
89 fCentralityType("V0A"), fHelperClass(0), fInitializedLocal(0),
90 fZVertexCut(10.0), fhJetPt(0x0), fhCuts(0x0), fhTrackPt(0x0),fhJetConstituentPt(0x0),fhJetEtaPt(0x0),fhAktJetEtaPhi(0x0),fhKtJetEtaPhi(0x0),fhJetAreaPt(0x0),fhJetPhiPt(0x0),fhZVertex(0x0),fhYVertex(0x0),fhXVertex(0x0),fhJetPtRho(0x0),fhJetPtConeRho(0x0),fhJetPtCMSRho(0x0),fhRho(0x0),fhConeRho(0x0),fhCMSRho(0x0),fhTrackEtaPt(0x0),fhGenTrackEtaPt(0x0),fhTrackPhiPt(0x0),fhTrackEtaPhi(0x0),fhKTJetPt(0x0),fhZVertexBC(0x0), fhRemx(0x0), fhPrimGenTrkPt(0x0), fhGenJetPt(0x0), fhRecTrkPt(0x0), fhFakeTrkPt(0x0),fhMult(0x0),fhTrackPhiCG(0x0),fhTrackPhiTPCG(0x0),fhAtimesRhoMedian(0x0),fhAtimesRhoCone(0x0),fhAtimesRhoCMS(0x0)
91 {
92  //Arrays initiation
93  for(int i=0;i<2;i++){
94  fhInvPtQVsPhi[i] = NULL;
95  fhInvPtQVsEta[i] = NULL;
96  fhInvPtQVsPhiASide[i] = NULL;
97  fhInvPtQVsPhiCSide[i] = NULL;
98  fhSigmaPtOverPtVsPt[i] = NULL;
99  }
100 }
101 
104 AliAnalysisTaskEmcalJet(name,kTRUE),
105  fUseDefaultVertexCut(1), fUsePileUpCut(1), fIsMC(0),
106  fSignalJetRadius(0.4),
107 fSignalJetEtaWindow(0.9 - fSignalJetRadius), fTrackEtaWindow(0.9), fMinTrackPt(0.150), fMinJetArea(0.0),
108 fCentralityType("V0A"), fHelperClass(0), fInitializedLocal(0),
109 fZVertexCut(10.0),fhJetPt(0x0),fhCuts(0x0), fhTrackPt(0x0),fhJetConstituentPt(0x0),fhJetEtaPt(0x0),fhAktJetEtaPhi(0x0),fhKtJetEtaPhi(0x0),fhJetAreaPt(0x0),fhJetPhiPt(0x0),fhZVertex(0x0),fhYVertex(0x0),fhXVertex(0x0),fhJetPtRho(0x0),fhJetPtConeRho(0x0),fhJetPtCMSRho(0x0),fhRho(0x0),fhConeRho(0x0),fhCMSRho(0x0),fhTrackEtaPt(0x0),fhTrackPhiPt(0x0),fhTrackEtaPhi(0x0),fhGenTrackEtaPt(0x0),fhKTJetPt(0x0),fhZVertexBC(0x0), fhRemx(0x0), fhPrimGenTrkPt(0x0), fhGenJetPt(0x0), fhRecTrkPt(0x0), fhFakeTrkPt(0x0),fhMult(0x0),fhTrackPhiCG(0x0),fhTrackPhiTPCG(0x0),fhAtimesRhoMedian(0x0),fhAtimesRhoCone(0x0),fhAtimesRhoCMS(0x0)
110 
111 {
112  //Arrays initiation
113  for(int i=0;i<2;i++){
114  fhInvPtQVsPhi[i] = NULL;
115  fhInvPtQVsEta[i] = NULL;
116  fhInvPtQVsPhiASide[i] = NULL;
117  fhInvPtQVsPhiCSide[i] = NULL;
118  fhSigmaPtOverPtVsPt[i] = NULL;
119  }
120 
121 
122  DefineOutput(1, AliEmcalList::Class());
123 }
124 
125 //________________________________________________________________________
127  //EVENT SELECTION RECONSTRUCTED DATA
128 
129  if(!event) return kFALSE;
130 
131  //___________________________________________________
132  //TEST PILE UP
133  if(fUsePileUpCut){
134  if(!fHelperClass || fHelperClass->IsPileUpEvent(event)){
135  return kFALSE;
136  }
137  if(fHelperClass && !fHelperClass->IsPileUpEvent(event)) fhCuts->Fill(1.5);//events that passed the pileup cut
138  }
139 
140  //___________________________________________________
141  //BEFORE VERTEX CUT
142  fhZVertexBC->Fill(event->GetPrimaryVertex()->GetZ());
144  if(!fHelperClass || !fHelperClass->IsVertexSelected2013pA(event)){
145  return kFALSE;
146  }
147  if(fHelperClass && fHelperClass->IsVertexSelected2013pA(event)) fhCuts->Fill(2.5);//events that passed the vertex cut
148  }else{
149 
150  if(TMath::Abs(event->GetPrimaryVertex()->GetZ()) > fZVertexCut){
151  return kFALSE;
152  }
153  else fhCuts->Fill(2.5);//events that passed the vertex cut
154  }
155  fhZVertex->Fill(event->GetPrimaryVertex()->GetZ());
156  fhXVertex->Fill(event->GetPrimaryVertex()->GetX());
157  fhYVertex->Fill(event->GetPrimaryVertex()->GetY());
158 
159  //___________________________________________________
160  //AFTER VERTEX CUT
161  return kTRUE;
162 }
163 
164 //________________________________________________________________________
166  // CHECK THE TRACK PT AND ETA RANGE
167  if(!track) return kFALSE;
168 
169  if(isprimary) {
170  if(!track->Charge()) return kFALSE;
171  if(!(static_cast<AliAODMCParticle*>(track))->IsPhysicalPrimary()) return kFALSE;
172  }
173 
174  if(TMath::Abs(track->Eta()) <= fTrackEtaWindow){ //APPLY TRACK ETA CUT
175  if(track->Pt() >= fMinTrackPt){ //APPLY TRACK MINIMUM PT CUT
176  return kTRUE;
177  }
178  }
179  return kFALSE;
180 }
181 
182 //________________________________________________________________________
184  //select jets in acceptance
185  if(!jet) return kFALSE;
186  if(TMath::Abs(jet->Eta()) <= fSignalJetEtaWindow){
187  if(jet->Area() >= fMinJetArea){
188  if(suppressGhost){
189  if(jet->Pt() >= fMinTrackPt) return kTRUE;
190  }else{
191  return kTRUE;
192  }
193  }
194  }
195  return kFALSE;
196 }
197 
200  // Initialization of jet containers done in AliAnalysisTaskEmcalJet::ExecOnce()
201  //Read arrays of jets and tracks
202  fInitializedLocal = kTRUE;
203 
204  // Initialize helper class (for vertex selection & pile up correction)
205  fHelperClass = new AliAnalysisUtils();
206  fHelperClass->SetCutOnZVertexSPD(kFALSE); // kFALSE: no cut; kTRUE: |zvtx-SPD - zvtx-TPC|<0.5cm
207 
208  return;
209 }
210 
214  if(!InputEvent()){
215  AliError("??? Event pointer == 0 ???");
216  return kFALSE;
217  }
218  fhCuts->Fill(0.5);
219  //Execute only once: Get tracks, jets from arrays if not already given
221 
222  //_________________________________________________________________
223  // FILL EVENT STATISTICS
224 
225  //Select events (vertex, pile-up,...)
226  if(!IsEventInAcceptance(InputEvent())) return kFALSE; //post data is in UserExec
227 
228  fhCuts->Fill(3.5); // to co prezilo
229 
230 
231  // JET+TRACK CONTAINERS
232  AliJetContainer *jetContRecAKT = NULL; //AKTjet container from reconstruced tracks
233  AliJetContainer *jetContRecKT = NULL; //KTjet container from reconstruced tracks
234  AliJetContainer *jetContGenAKT = NULL; //AKT jet container from generated MC tracks
235  AliJetContainer *jetContGenKT = NULL; //KT jet container from generated MC tracks
236 
237  AliEmcalJet *jetRec = NULL;
238  AliEmcalJet *jetGen = NULL;//MC
239 
240  AliParticleContainer *trkContRec = NULL; //track array of real reconstructed tracks
241  AliParticleContainer *parContGen = NULL; //track array of generated MC tracks
242 
243  AliVParticle *constTrackRec = NULL; //rec jet constituent
244  AliVParticle *constTrackGen = NULL; //gen jet constituent
245  //_________________________________________________________
246  //READ JET TRACK CONTAINERS
247  jetContRecAKT = GetJetContainer(kContainerOne); //AKT jet container
248  jetContRecKT = GetJetContainer(kContainerTwo); //KT jet container
249  if(fIsMC){
250  jetContGenAKT = GetJetContainer(kContainerThree); //MC AKT jet container
251  jetContGenKT = GetJetContainer(kContainerFour); //MC KT jet container
252  }
253 
254  trkContRec = GetParticleContainer(kContainerOne); //reconstructed hybrid tracks
255  if(fIsMC) parContGen = GetParticleContainer(kContainerTwo); //MC particles
256 
257  Int_t fMultCounter;
258  Double_t xyz[50];
259  Double_t pxpypz[50];
260  Double_t cv[21];
261  Int_t itrkq;
262 
263  AliAODTrack *constTrackRecAod=NULL;
264 
265 // Lop over reconstructed trcks
266  if(trkContRec){
267  fMultCounter=0;
268  for (auto trackIterator : trkContRec->accepted_momentum()){
269  constTrackRecAod = (AliAODTrack*)trackIterator.second;
270  if(!constTrackRecAod->IsHybridGlobalConstrainedGlobal()) continue;
271  if(!constTrackRecAod) continue;
272  if(!IsTrackInAcceptance(constTrackRecAod)) continue;
273  Double_t phi = Convert(constTrackRecAod->Phi());
274  fhTrackPt->Fill(constTrackRecAod->Pt());
275  fhTrackEtaPt->Fill(constTrackRecAod->Eta(),constTrackRecAod->Pt());
276  fhTrackPhiPt->Fill(phi,constTrackRecAod->Pt());
277  fhTrackEtaPhi->Fill(constTrackRecAod->Eta(),phi);
278 
279  if(constTrackRecAod->IsGlobalConstrained()){
280  fhTrackPhiCG->Fill(constTrackRecAod->Pt(), phi); //global constrained
281  }else{
282  fhTrackPhiTPCG->Fill(constTrackRecAod->Pt(), phi); //complementary
283  }
284  itrkq = (constTrackRecAod->Charge()<0) ? 0 : 1;
285 
286  fhInvPtQVsPhi[itrkq]->Fill(phi, 1.0/constTrackRecAod->Pt());
287  fhInvPtQVsEta[itrkq]->Fill(constTrackRecAod->Eta(), 1.0/constTrackRecAod->Pt());
288 
289  if(constTrackRecAod->Eta()>0){
290  fhInvPtQVsPhiASide[itrkq]->Fill(phi, 1.0/constTrackRecAod->Pt());
291  }else{
292  fhInvPtQVsPhiCSide[itrkq]->Fill(phi, 1.0/constTrackRecAod->Pt());
293  }
294  memset(cv, 0, sizeof(Double_t) * 21); //cleanup arrays
295  memset(pxpypz, 0, sizeof(Double_t) * 50);
296  memset(xyz, 0, sizeof(Double_t) * 50);
297  constTrackRecAod->GetXYZ(xyz);
298  constTrackRecAod->GetPxPyPz(pxpypz);
299  constTrackRecAod->GetCovarianceXYZPxPyPz(cv);
300  AliExternalTrackParam par(xyz, pxpypz, cv, constTrackRecAod->Charge());
301  fhSigmaPtOverPtVsPt[itrkq]->Fill(constTrackRecAod->Pt(), TMath::Abs(sqrt(par.GetSigma1Pt2())/par.GetSigned1Pt()));
302  fMultCounter++;
303  }
304  fhMult->Fill(fMultCounter);
305  }
306  //_________________________________________________
307  // Background
308  Double_t rho = EstimateBgKT(jetContRecKT);
309  fhRho->Fill(rho);
310 
311  Double_t conerho = EstimateLocalBg(jetContRecAKT,trkContRec);
312  fhConeRho->Fill(conerho);
313 
314  Double_t cmsrho = EstimateBgKTCMS(jetContRecKT);
315  fhCMSRho->Fill(cmsrho);
316  //_________________________________________________
317  // Loop over AKT jets
318  if(jetContRecAKT){
319  jetContRecAKT->ResetCurrentID();
320  for(auto jetIterator : jetContRecAKT->accepted_momentum()) {//loop over reconstructed jets
321  jetRec=(AliEmcalJet*)jetIterator.second;
322  if(!jetRec) continue;
323  if(!IsSignalJetInAcceptance(jetRec)) continue; //check jet eta and pt
324  Double_t phi = Convert(jetRec->Phi());
325  fhJetPt->Fill(jetRec->Pt()); //fill AKT jet pT
326  fhJetPtRho->Fill(jetRec->Pt()-rho*jetRec->Area()); //fill AKT jet pT without bkg
327  fhJetPtConeRho->Fill(jetRec->Pt()-conerho*jetRec->Area()); //fill AKT jet pT without local bkg
328  fhJetPtCMSRho->Fill(jetRec->Pt()-cmsrho*jetRec->Area()); //fill AKT jet pT without CMS bkg
329  fhJetEtaPt->Fill(jetRec->Eta(),jetRec->Pt());
330  fhAktJetEtaPhi->Fill(jetRec->Eta(),phi);
331  fhJetPhiPt->Fill(phi,jetRec->Pt());
332  fhJetAreaPt->Fill(jetRec->Area(),jetRec->Pt());
333  fhAtimesRhoMedian->Fill(rho*jetRec->Area());
334  fhAtimesRhoCone->Fill(conerho*jetRec->Area());
335  fhAtimesRhoCMS->Fill(cmsrho*jetRec->Area());
336 
337  for(Int_t iq=0; iq < jetRec->GetNumberOfTracks(); iq++){
338  constTrackRec = static_cast<AliVParticle*> (jetRec->TrackAt(iq, trkContRec->GetArray())); //matched rec and emb tracks
339  if(!constTrackRec) continue;
340  fhJetConstituentPt->Fill(constTrackRec->Pt());
341 
342  }
343  }
344  }
345  // Loop over KT jets
346  if(jetContRecKT){
347  for(auto jetIterator : jetContRecKT->accepted_momentum()) {
348  jetRec=(AliEmcalJet*)jetIterator.second;
349  if(!jetRec) continue;
350  if(!IsSignalJetInAcceptance(jetRec)) continue;
351 
352  fhKTJetPt->Fill(jetRec->Pt());
353  }
354  }
355 
356  // Single particle efficiency and contamination
357  //primary particles spectrum
358  if(parContGen){
359  for(auto trackIterator : parContGen->accepted_momentum()) {
360  constTrackGen=(AliVParticle*)trackIterator.second;
361  if(!constTrackGen) continue;
362  if(IsTrackInAcceptance(constTrackGen, kTRUE)){
363  fhPrimGenTrkPt->Fill(constTrackGen->Pt()); //pT spectrum of generator level particles eta-pt, phi-pt, eta-phi histogram
364  fhGenTrackEtaPt->Fill(constTrackGen->Eta(),constTrackGen->Pt()); //pT spectrum of generator level particles eta-pt, phi-pt, eta-phi histogram
365  }
366  }
367  }
368 
369 
370  //reconstructed primary + secondary particles
371  Bool_t bRecPrim = kFALSE; //tags the reconstructed primary particles
372  if(trkContRec && parContGen){
373  for(auto trackIterator : trkContRec->accepted_momentum()) {
374  constTrackRecAod=(AliAODTrack*)trackIterator.second;
375  if(!constTrackRecAod->IsHybridGlobalConstrainedGlobal()) continue;
376  if(!constTrackRecAod) continue;
377  if(!IsTrackInAcceptance(constTrackRecAod, kFALSE)) continue; //reconstructed level tracks
378  bRecPrim = kFALSE; //reset matching flag
379 
380  parContGen->ResetCurrentID();
381  for(auto trackIterator : parContGen->accepted_momentum()) {
382  constTrackGen=(AliAODTrack*)trackIterator.second;
383  if(!constTrackGen) continue;
384  if(!IsTrackInAcceptance(constTrackGen, kTRUE)) continue; //gen level physical primary
385  if(TMath::Abs(constTrackRecAod->GetLabel()) == TMath::Abs(constTrackGen->GetLabel())){
386  fhRecTrkPt->Fill(constTrackGen->Pt());
387  bRecPrim = kTRUE; //matched
388  break;
389  }
390  }
391  if(!bRecPrim){
392  fhFakeTrkPt->Fill(constTrackRecAod->Pt());//Pt of fake tracks
393  }
394  }
395  }
396  // Response matrix
397  Double_t ptGenCorr; //GEN jet pt corrected for rho
398  Double_t ptRecCorr; //REC jet pt corrected for rho
399 
400  //Response matrix normalization - spectrum of all generator level jets in acceptance
401  if(jetContGenAKT){
402  for(auto jetIterator : jetContGenAKT->accepted_momentum()) {
403  jetGen=(AliEmcalJet*)jetIterator.second;
404  if(!jetGen) continue;
405  if(!IsSignalJetInAcceptance(jetGen,kTRUE)) continue; //cuts on eta, pT ,are
406  fhGenJetPt->Fill(jetGen->Pt()); //Pt spectrum of MC jets
407  }
408  }
409  //Find closest gen level+rec level jets
410  if(jetContRecAKT){
411  for(auto jetIterator : jetContRecAKT->accepted_momentum()) {
412  jetRec=(AliEmcalJet*)jetIterator.second;
413  if(!jetRec) continue;
414  if(!IsSignalJetInAcceptance(jetRec,kTRUE)) continue; //cuts on eta, pT ,area
415 
416  jetGen = 0x0;
417  jetGen = jetRec->ClosestJet();
418 
419  if(!jetGen){ //did not find matching generator level jet
420  continue;
421  }
422  fhRemx->Fill(jetRec->Pt(),jetGen->Pt());//response matrix
423  }
424  }
425 
426 
427  return kTRUE;
428 }
429 
430 //________________________________________________________________________
432  //Treminate
433  PostData(1, fOutput);
434 
435  // Mandatory
436  fOutput = dynamic_cast<AliEmcalList*> (GetOutputData(1)); // '1' refers to the output slot
437  if(!fOutput) {
438  printf("ERROR: Output list not available\n");
439  return;
440  }
441 }
442 
443 //________________________________________________________________________
445  // Destructor. Clean-up the output list, but not the histograms that are put inside
446  if(fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
447  delete fOutput;
448  }
449  delete fHelperClass;
450 
451 }
452 
456  Bool_t oldStatus = TH1::AddDirectoryStatus();
457  TH1::AddDirectory(kFALSE);
458  TString name;
459 
460  //Histograms
461  fhJetPt = new TH1D("fhJetPt","Pt spectrum of the AKT jets",320,-20.0,300.0);
462  fhJetPtRho = new TH1D("fhJetPtRho","Pt spectrum of the AKT jets - kt bg",320,-20.0,300.0);
463  fhJetPtConeRho = new TH1D("fhJetPtConeRho","p_{T} spectrum of the a-k_{T} jets - local kt bg",320,-20.0,300.0);
464  fhJetPtCMSRho = new TH1D("fhJetPtCMSRho","p_{T} spectrum of the a-k_{T} jets - CMS bg",320,-20.0,300.0);
465  fhAtimesRhoMedian = new TH1D("fhAtimesRhoMedian","a-k_{T} background times area",320,-20.0,300.0);
466  fhAtimesRhoCone = new TH1D("fhAtimesRhoCone","cone background times area",320,-20.0,300.0);
467  fhAtimesRhoCMS = new TH1D("fhAtimesRhoCMS","CMS a-k_{T} background times area",320,-20.0,300.0);
468  fhJetConstituentPt = new TH1D("fhJetConstituentPt","Pt spectrum of the a-k_{T} jet constituents",300,0.0,300.0);
469  fhTrackPt = new TH1D("fhTrackPt","p_{T} spectrum of tracks",300,0.0,300.0);
470  fhJetAreaPt = new TH2D("fhJetAreaPt","Jet area vs pt",50,0.,2.,300,0.0,300.0);
471  fhJetEtaPt = new TH2D("fhJetEtaPt","#eta-p_{T} distribution of jets",50,-fSignalJetEtaWindow,fSignalJetEtaWindow,300,0.0,300.0);
472  fhAktJetEtaPhi = new TH2D("fhAktJetEtaPhi","Eta-Phi distribution of a-k_{T} jets",50,-fSignalJetEtaWindow,fSignalJetEtaWindow,50,-TMath::Pi(),TMath::Pi());
473  fhKtJetEtaPhi = new TH2D("fhKtJetEtaPhi","#eta-#phi distribution of k_{T} jets",50,-fSignalJetEtaWindow,fSignalJetEtaWindow,50,-TMath::Pi(),TMath::Pi());
474  fhTrackEtaPhi = new TH2D("fhTrackEtaPhi","#eta-#phi distribution of tracks",50,-fTrackEtaWindow,fTrackEtaWindow,50,-TMath::Pi(),TMath::Pi());
475  fhTrackEtaPt = new TH2D("fhTrackEtaPt","#eta-p_{T} distribution of tracks",50,-fTrackEtaWindow,fTrackEtaWindow,300,0.0,300.0);
476  fhGenTrackEtaPt = new TH2D("fhGenTrackEtaPt","#eta-p_{T} distribution of tracks",50,-fTrackEtaWindow,fTrackEtaWindow,300,0.0,300.0);
477  fhJetPhiPt = new TH2D("fhJetPhiPt","#phi-p_{T} distribution of jets",50,-TMath::Pi(),TMath::Pi(),300,0.0,300.0);
478  fhTrackPhiPt = new TH2D("fhTrackPhiPt","#phi-p_{T} distribution of tracks",50,-TMath::Pi(),TMath::Pi(),300,0.0,300.0);
479  fhZVertex = new TH1D("fhZVertex","Z vertex",300,-15.0,15.0);
480  fhZVertexBC = new TH1D("fhZVertexBC","Z vertex before the cut",200,-50.0,50.0);
481  fhXVertex = new TH1D("fhXVertex","X vertex",300,-15.0,15.0);
482  fhYVertex = new TH1D("fhYVertex","Y vertex",300,-15.0,15.0);
483  fhCuts = new TH1D("fhCuts","Cuts statistics",4,0,4);
484  fhRho = new TH1D("fhRho","k_{T} jet background",500,0,50.0);
485  fhConeRho = new TH1D("fhConeRho","Local a-k_{T} jet background",500,0,50.0);
486  fhCMSRho = new TH1D("fhCMSRho","CMS-style a-k_{T} jet background",500,0,50.0);
487  fhKTJetPt = new TH1D("fhKTJetPt","KT jets p_{T} spectrum",300,0.0,300.0);
488  fhPrimGenTrkPt = new TH1D("fhPrimGenTrkPt","MC track p_{T}",300,0.0,300.0);
489  fhGenJetPt = new TH1D("fhGenJetPt","MC jet p_{T}",300,0.0,300.0);
490  fhRecTrkPt = new TH1D("fhRecTrkPt","Correctly rec. track p_{T}",300,0.0,300.0);
491  fhFakeTrkPt = new TH1D("fhFakeTrkPt","Fake track p_{T}",300,0.0,300.0);
492  fhRemx = new TH2D("fhRemx","Response matrix",300,0.0,300.0,300,0.0,300.0);
493  fhMult = new TH1D("fhMult","Reconstructed track multiplicity",300,0,300.0);
494 
495  fhTrackPhiCG = new TH2D("fhTrackPhiCG","Global tracks azimuth vs p_{T} ",300,0,300,50,-TMath::Pi(),TMath::Pi());
496  fhTrackPhiTPCG = new TH2D("fhTrackPhiTPCG","Complementary tracks azimuth vs p_{T)",300,0,300.0,50,-TMath::Pi(),TMath::Pi());
497  for(int i=0;i<2;i++){
498  name = Form("fhSigmaPtOverPtVsPt%i",i);
499  fhSigmaPtOverPtVsPt[i] = new TH2D(name.Data(),"#sigma p_{T}/p_{T} vs p_{T}",100,0.0,100.0,300,0.0,3.0);
500  fOutput->Add(fhSigmaPtOverPtVsPt[i]);
501  name = Form("fhInvPtQVsPhi%i",i);
502  fhInvPtQVsPhi[i] = new TH2D(name.Data(),"Negative track azimuth vs 1/pt",50,-TMath::Pi(),TMath::Pi(),300,0.0,0.25);
503  fOutput->Add(fhInvPtQVsPhi[i]);
504  name = Form("fhInvPtQVsEta%i",i);
505  fhInvPtQVsEta[i] = new TH2D(name.Data(),"Negative track azimuth vs 1/pt",50,-fTrackEtaWindow,fTrackEtaWindow,300,0.0,0.25);
506  fOutput->Add(fhInvPtQVsEta[i]);
507  name = Form("fhInvPtQVsPhiASide%i",i);
508  fhInvPtQVsPhiASide[i] = new TH2D(name.Data(),"Negative phi vs 1/pt A-eta",50,-TMath::Pi(),TMath::Pi(),300,0.0,0.25);
509  fOutput->Add(fhInvPtQVsPhiASide[i]);
510  name = Form("fhInvPtQVsPhiCSide%i",i);
511  fhInvPtQVsPhiCSide[i] = new TH2D(name.Data(),"Negative phi vs 1/pt C-eta",50,-TMath::Pi(),TMath::Pi(),300,0.0,0.25);
512  fOutput->Add(fhInvPtQVsPhiCSide[i]);
513  }
514 
515  fOutput->Add(fhJetPt);
516  fOutput->Add(fhJetPtRho);
517  fOutput->Add(fhJetPtConeRho);
518  fOutput->Add(fhJetPtCMSRho);
520  fOutput->Add(fhAtimesRhoCone);
521  fOutput->Add(fhAtimesRhoCMS);
522  fOutput->Add(fhJetAreaPt);
523 
525  fOutput->Add(fhTrackPt);
526  fOutput->Add(fhJetEtaPt);
527  fOutput->Add(fhTrackEtaPt);
528  fOutput->Add(fhJetPhiPt);
529 
530  fOutput->Add(fhTrackPhiPt);
531  fOutput->Add(fhZVertex);
532  fOutput->Add(fhXVertex);
533  fOutput->Add(fhYVertex);
534  fOutput->Add(fhCuts);
535 
536  fOutput->Add(fhRho);
537  fOutput->Add(fhConeRho);
538  fOutput->Add(fhCMSRho);
539  fOutput->Add(fhKTJetPt);
540  fOutput->Add(fhZVertexBC);
541  fOutput->Add(fhPrimGenTrkPt);
542 
543  fOutput->Add(fhGenJetPt);
544  fOutput->Add(fhRecTrkPt);
545  fOutput->Add(fhFakeTrkPt);
546  fOutput->Add(fhRemx);
547  fOutput->Add(fhMult);
548  fOutput->Add(fhGenTrackEtaPt);
549 
550  fOutput->Add(fhTrackPhiCG);
551  fOutput->Add(fhTrackPhiTPCG);
552  fOutput->Add(fhAktJetEtaPhi);
553  fOutput->Add(fhKtJetEtaPhi);
554  fOutput->Add(fhTrackEtaPhi);
555  // =========== Switch on Sumw2 for all histos ===========
556  for(Int_t i=0; i<fOutput->GetEntries(); i++){
557  TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
558  if(h1){
559  h1->Sumw2();
560  continue;
561  }
562  THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
563  if(hn){
564  hn->Sumw2();
565  }
566  }
567  TH1::AddDirectory(oldStatus);
568 
569 
570  PostData(1, fOutput);
571 }
572 //________________________________________________________________________
574  //
575  // retrieve event objects
576  //
578 
579  return kTRUE;
580 }
583 {
584 
585  return kTRUE;
586 }
587 
592  Double_t rhoCone = 0.0;
593 
594  if(!jetCont) return rhoCone;
595 
596  AliEmcalJet* jet = NULL;
597  static Double_t rhovec[999];
598  Int_t nJetAcc = 0;
599  Double_t jetpt;
600 
601  for(auto jetIterator : jetCont->accepted_momentum()) {
602  jet=(AliEmcalJet*)jetIterator.second;
603  if(!jet) continue;
604  if(!IsSignalJetInAcceptance(jet,kFALSE)) continue; //check jet eta and pt
605 
606  jetpt = jet->Pt();
607  if(jetpt <0.005) jetpt = 0.; //set pt of ghost jets identical to zero
608  rhovec[nJetAcc] = jetpt/jet->Area();
609  nJetAcc++;
610  }
611 
612  if(nJetAcc>0){
613  rhoCone = TMath::Median(nJetAcc, rhovec); //take the median of pTjet/Ajet
614  }
615 
616  return rhoCone;
617 }
623  Double_t rhoCone = 0.0;
624 
625  if(!jetCont) return rhoCone;
626 
627  AliEmcalJet* jet = NULL;
628  AliVParticle* track = NULL;
629  Double_t trackpt=0;
630  Double_t tracketa=0;
631  Double_t trackphi=0;
632  Double_t coner = 0.4;//the radius of the perpendicular cone
633  Double_t coner2 = coner*coner;
634  Double_t highjetpt=-1;
635  Double_t highjeteta=0;
636  Double_t highjetphi=0;
637  Double_t perpphi1=0;
638  Double_t perpphi2=0;
639  Double_t bgpt=0;
640 
641  Double_t X=0;
642  Double_t Y=0;
643  Double_t Z=0;
644 
645  for(auto jetIterator : jetCont->accepted_momentum()) {
646  jet=(AliEmcalJet*)jetIterator.second;
647  if(!jet) continue;
648  if(!IsSignalJetInAcceptance(jet)) continue; //check jet eta and pt
649  if(jet->Pt()>highjetpt){
650  highjetpt = jet->Pt();
651  highjeteta = jet->Eta();
652  highjetphi = Convert(jet->Phi());
653  }
654  }
655  if(highjetpt<0) return 0.;
656  perpphi1=Convert(highjetphi+TMath::Pi()/2);//the two perpendicular angles
657  perpphi2=Convert(highjetphi-TMath::Pi()/2);
658 
659  for(auto trackIterator : trkCont->accepted_momentum()) {
660  track=(AliVParticle*)trackIterator.second;
661  if(!((AliAODTrack*)track)->IsHybridGlobalConstrainedGlobal()) continue;
662  if(!track) continue;
663  if(!IsTrackInAcceptance(track)) continue;
664 
665  trackpt = track->Pt();
666  tracketa = track->Eta();
667  trackphi = Convert(track->Phi());
668  X=tracketa-highjeteta;
669  Y=Convert(trackphi-perpphi1);
670  Z=Convert(trackphi-perpphi2);
671 
672  if((X*X+Y*Y) < coner2 || (X*X+Z*Z) < coner2){
673  bgpt=bgpt+trackpt;
674  }
675  }
676  rhoCone = bgpt/(2*TMath::Pi()*coner2); //calculate the pT density
677 
678  return rhoCone;
679 }
684  Double_t rhoCMS = 0.0;
685 
686  if(!jetCont) return rhoCMS;
687 
688  AliEmcalJet* jet = NULL;
689  static Double_t rhovec[999];
690  Int_t nJetAcc = 0;
691  Double_t jetpt = 0.0;
692  Double_t physJetArea = 0.0;
693  Double_t allJetArea = 0.0;
694 
695  for(auto jetIterator : jetCont->accepted_momentum()) {
696  jet=(AliEmcalJet*)jetIterator.second;
697  if(!jet) continue;
698  if(!IsSignalJetInAcceptance(jet,kFALSE)) continue; //check jet eta and pt
699 
700  if(jet->Pt() > 0.01){//Select only physical jets
701  physJetArea += jet->Area(); //Area of physical jets
702  jetpt = jet->Pt();
703  rhovec[nJetAcc] = jetpt/jet->Area();
704  nJetAcc++;
705  }
706  allJetArea += jet->Area();//Area of all jets
707  }
708 
709  if(nJetAcc>0 && allJetArea!=0){
710  rhoCMS = TMath::Median(nJetAcc, rhovec)*physJetArea/allJetArea; //Calculate the CMS background
711  }
712 
713  return rhoCMS;
714 }
718 Double_t AliAnalysisTaskJetPP::Convert(Double_t input){//Converts an angle to the range (-Pi,Pi)
719  Double_t phi;
720  phi = TMath::ATan2(TMath::Sin(input),TMath::Cos(input));
721  return phi;
722 }
void ExecOnceLocal()
Perform steps needed to initialize the analysis.
Double_t EstimateBgKT(AliJetContainer *jetCont)
Bool_t IsSignalJetInAcceptance(AliEmcalJet *jet, Bool_t suppressGhost=1)
Double_t Area() const
Definition: AliEmcalJet.h:130
double Double_t
Definition: External.C:58
TH1D * fhXVertex
! X vertex
TH2D * fhRemx
! Response matrix
Bool_t IsTrackInAcceptance(AliVParticle *track, Bool_t isprimary=0)
TH2D * fhTrackEtaPhi
! Track eta-phi distribution
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:327
TH2D * fhGenTrackEtaPt
! Generated track eta-pt distribution
AliJetContainer * GetJetContainer(Int_t i=0) const
Double_t Eta() const
Definition: AliEmcalJet.h:121
TH1D * fhZVertex
! Z vertex
TH1D * fhJetPtConeRho
! pt spectrum of AKT jets without local bgk
Double_t Phi() const
Definition: AliEmcalJet.h:117
TH2D * fhInvPtQVsEta[2]
! Inverse pt vs eta
Bool_t Run()
Run the nalysis.
TH2D * fhTrackPhiCG
! Global tracks
TH2D * fhTrackPhiTPCG
! complementary tracks
Double_t fMinJetArea
Min jet area to be accepted.
TH1D * fhRecTrkPt
! Pt spectrum of correctly reconstructed tracks
TH1D * fhPrimGenTrkPt
! Pt spectrum of MC particle
TH2D * fhTrackEtaPt
! Track eta-pt distribution
TH1D * fhMult
! Multiplicity
Declaration of the class AliAnalysisTaskJetPP.
TH2D * fhInvPtQVsPhiCSide[2]
! Inverse pt vs pi in close-side rapidity
TH1D * fhTrackPt
! Track pt
Double_t EstimateLocalBg(AliJetContainer *jetCont, AliParticleContainer *trkCont)
TH1D * fhJetConstituentPt
! Jet constituent pt
Double_t fSignalJetEtaWindow
+- window in eta for signal jets
Container for particles within the EMCAL framework.
TH1D * fhCMSRho
! CMS KT bgk
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:160
Pure virtual base class for AliAnalysisTaskJetPP.
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:139
TH1D * fhJetPt
! pt spectrum of AKT jets
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
TH2D * fhAktJetEtaPhi
! Jet et-phi distribution
Double_t Convert(Double_t input)
TH1D * fhAtimesRhoCone
! cone jet background times area
const AliJetIterableMomentumContainer accepted_momentum() const
int Int_t
Definition: External.C:63
Double_t fTrackEtaWindow
gc +- window in eta for tracks
TH2D * fhTrackPhiPt
! Track phi-pt distribution
TH2D * fhJetAreaPt
! Jet Area-pt distribution
Double_t fZVertexCut
zvertex cut
TH2D * fhKtJetEtaPhi
! Jet et-phi distribution
TH1D * fhJetPtRho
! pt spectrum of AKT jets without kt bgk
Definition: External.C:228
void UserCreateOutputObjects()
Overloads base class method. Creates output objects.
Definition: External.C:212
TH2D * fhInvPtQVsPhiASide[2]
! Inverse pt vs phi in away-side rapidity
TH1D * fhConeRho
! local KT bgk
TH1D * fhAtimesRhoCMS
! CMS jet background times area
TH2D * fhJetPhiPt
! Jet phi-pt distribution
TH1D * fhYVertex
! Y vertex
Bool_t fInitializedLocal
! gc trigger if tracks/jets are loaded initiates calling ExecOnce
TH1D * fhKTJetPt
! pt spectrum of KT jets
TH2D * fhSigmaPtOverPtVsPt[2]
! Pt vs sigmaPt/Pt
Double_t Pt() const
Definition: AliEmcalJet.h:109
Enhanced TList-derived class that implements correct merging for pt_hard binned production.
Definition: AliEmcalList.h:25
TH1D * fhCuts
! histogram for pilup/vertex cuts
TH2D * fhInvPtQVsPhi[2]
! Inverse pt vs phi
AliEmcalList * fOutput
!output list
AliAnalysisUtils * fHelperClass
! gc Vertex selection helper
Bool_t RetrieveEventObjects()
Retrieve common objects from event.
Double_t EstimateBgKTCMS(AliJetContainer *jetCont)
Double_t fMinTrackPt
gc Min track pt to be accepted
const AliParticleIterableMomentumContainer accepted_momentum() const
TH1D * fhJetPtCMSRho
! pt spectrum of AKT jets without CMS bgk
TH1D * fhAtimesRhoMedian
! kt jet background times area
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
TH1D * fhZVertexBC
! Z vertex before cut
const char Option_t
Definition: External.C:48
void UserCreateOutputObjects()
Main initialization function on the worker.
Bool_t IsEventInAcceptance(AliVEvent *event)
bool Bool_t
Definition: External.C:53
TH1D * fhGenJetPt
! Pt spectrum of MC jets
Bool_t fUseDefaultVertexCut
trigger if automatic vertex cut from helper class should be done
Bool_t fUsePileUpCut
trigger if pileup cut should be done
Container for jet within the EMCAL jet framework.
Definition: External.C:196
AliAnalysisTaskJetPP()
Default constructor.
TH1D * fhFakeTrkPt
! Pt spectrum of fake tracks
TH2D * fhJetEtaPt
! Jet eta-pt distribution