AliPhysics  ff07904 (ff07904)
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), fASPDCvsTCut(65.), fBSPDCvsTCut(4.), 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), fASPDCvsTCut(65.), fBSPDCvsTCut(4.), 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  //BEFORE VERTEX CUT
133  fhZVertexBC->Fill(event->GetPrimaryVertex()->GetZ());
135  if(!fHelperClass || !fHelperClass->IsVertexSelected2013pA(event)){
136  return kFALSE;
137  }
138  if(fHelperClass && fHelperClass->IsVertexSelected2013pA(event)) fhCuts->Fill(1.5);//events that passed the vertex cut
139  }else{
140 
141  if(TMath::Abs(event->GetPrimaryVertex()->GetZ()) > fZVertexCut){
142  return kFALSE;
143  }
144  else fhCuts->Fill(1.5);//events that passed the vertex cut
145  }
146  fhZVertex->Fill(event->GetPrimaryVertex()->GetZ());
147  fhXVertex->Fill(event->GetPrimaryVertex()->GetX());
148  fhYVertex->Fill(event->GetPrimaryVertex()->GetY());
149 
150  //___________________________________________________
151  //AFTER VERTEX CUT
152  return kTRUE;
153 }
154 //PILEUP CUT
156  if(fUsePileUpCut){
157  Int_t nClustersLayer0 = event->GetNumberOfITSClusters(0);
158  Int_t nClustersLayer1 = event->GetNumberOfITSClusters(1);
159  Int_t nTracklets = event->GetMultiplicity()->GetNumberOfTracklets();
160  if (nClustersLayer0 + nClustersLayer1 > fASPDCvsTCut + nTracklets*fBSPDCvsTCut){
161  return kTRUE;
162  }
163  }
164  fhCuts->Fill(2.5);
165  return kFALSE;
166 }
167 
168 //________________________________________________________________________
170  // CHECK THE TRACK PT AND ETA RANGE
171  if(!track) return kFALSE;
172 
173  if(isprimary) {
174  if(!track->Charge()) return kFALSE;
175  if(!(static_cast<AliAODMCParticle*>(track))->IsPhysicalPrimary()) return kFALSE;
176  }
177 
178  if(TMath::Abs(track->Eta()) <= fTrackEtaWindow){ //APPLY TRACK ETA CUT
179  if(track->Pt() >= fMinTrackPt){ //APPLY TRACK MINIMUM PT CUT
180  return kTRUE;
181  }
182  }
183  return kFALSE;
184 }
185 
186 //________________________________________________________________________
188  //select jets in acceptance
189  if(!jet) return kFALSE;
190  if(TMath::Abs(jet->Eta()) <= fSignalJetEtaWindow){
191  if(jet->Area() >= fMinJetArea){
192  if(suppressGhost){
193  if(jet->Pt() >= fMinTrackPt) return kTRUE;
194  }else{
195  return kTRUE;
196  }
197  }
198  }
199  return kFALSE;
200 }
201 
204  // Initialization of jet containers done in AliAnalysisTaskEmcalJet::ExecOnce()
205  //Read arrays of jets and tracks
206  fInitializedLocal = kTRUE;
207 
208  // Initialize helper class (for vertex selection & pile up correction)
209  fHelperClass = new AliAnalysisUtils();
210  fHelperClass->SetCutOnZVertexSPD(kFALSE); // kFALSE: no cut; kTRUE: |zvtx-SPD - zvtx-TPC|<0.5cm
211 
212  return;
213 }
214 
218  if(!InputEvent()){
219  AliError("??? Event pointer == 0 ???");
220  return kFALSE;
221  }
222  fhCuts->Fill(0.5);
223  //Execute only once: Get tracks, jets from arrays if not already given
225 
226  //_________________________________________________________________
227  // FILL EVENT STATISTICS
228 
229  //Select events (vertex, pile-up,...)
230  if(!IsEventInAcceptance(InputEvent())) return kFALSE; //post data is in UserExec
231  if(IsSPDClusterVsTrackletBG(InputEvent())) return kFALSE;
232 
233 
234  // JET+TRACK CONTAINERS
235  AliJetContainer *jetContRecAKT = NULL; //AKTjet container from reconstruced tracks
236  AliJetContainer *jetContRecKT = NULL; //KTjet container from reconstruced tracks
237  AliJetContainer *jetContGenAKT = NULL; //AKT jet container from generated MC tracks
238  AliJetContainer *jetContGenKT = NULL; //KT jet container from generated MC tracks
239 
240  AliEmcalJet *jetRec = NULL;
241  AliEmcalJet *jetGen = NULL;//MC
242 
243  AliParticleContainer *trkContRec = NULL; //track array of real reconstructed tracks
244  AliParticleContainer *parContGen = NULL; //track array of generated MC tracks
245 
246  AliVParticle *constTrackRec = NULL; //rec jet constituent
247  AliVParticle *constTrackGen = NULL; //gen jet constituent
248  //_________________________________________________________
249  //READ JET TRACK CONTAINERS
250  jetContRecAKT = GetJetContainer(kContainerOne); //AKT jet container
251  jetContRecKT = GetJetContainer(kContainerTwo); //KT jet container
252  if(fIsMC){
253  jetContGenAKT = GetJetContainer(kContainerThree); //MC AKT jet container
254  jetContGenKT = GetJetContainer(kContainerFour); //MC KT jet container
255  }
256 
257  trkContRec = GetParticleContainer(kContainerOne); //reconstructed hybrid tracks
258  if(fIsMC) parContGen = GetParticleContainer(kContainerTwo); //MC particles
259 
260  Int_t fMultCounter;
261  Double_t xyz[50];
262  Double_t pxpypz[50];
263  Double_t cv[21];
264  Int_t itrkq;
265 
266  AliAODTrack *constTrackRecAod=NULL;
267 
268 // Lop over reconstructed trcks
269  if(trkContRec){
270  fMultCounter=0;
271  for (auto trackIterator : trkContRec->accepted_momentum()){
272  constTrackRecAod = (AliAODTrack*)trackIterator.second;
273  if(!constTrackRecAod->IsHybridGlobalConstrainedGlobal()) continue;
274  if(!constTrackRecAod) continue;
275  if(!IsTrackInAcceptance(constTrackRecAod)) continue;
276  Double_t phi = Convert(constTrackRecAod->Phi());
277  fhTrackPt->Fill(constTrackRecAod->Pt());
278  fhTrackEtaPt->Fill(constTrackRecAod->Eta(),constTrackRecAod->Pt());
279  fhTrackPhiPt->Fill(phi,constTrackRecAod->Pt());
280  fhTrackEtaPhi->Fill(constTrackRecAod->Eta(),phi);
281 
282  if(constTrackRecAod->IsGlobalConstrained()){
283  fhTrackPhiCG->Fill(constTrackRecAod->Pt(), phi); //global constrained
284  }else{
285  fhTrackPhiTPCG->Fill(constTrackRecAod->Pt(), phi); //complementary
286  }
287  itrkq = (constTrackRecAod->Charge()<0) ? 0 : 1;
288 
289  fhInvPtQVsPhi[itrkq]->Fill(phi, 1.0/constTrackRecAod->Pt());
290  fhInvPtQVsEta[itrkq]->Fill(constTrackRecAod->Eta(), 1.0/constTrackRecAod->Pt());
291 
292  if(constTrackRecAod->Eta()>0){
293  fhInvPtQVsPhiASide[itrkq]->Fill(phi, 1.0/constTrackRecAod->Pt());
294  }else{
295  fhInvPtQVsPhiCSide[itrkq]->Fill(phi, 1.0/constTrackRecAod->Pt());
296  }
297  memset(cv, 0, sizeof(Double_t) * 21); //cleanup arrays
298  memset(pxpypz, 0, sizeof(Double_t) * 50);
299  memset(xyz, 0, sizeof(Double_t) * 50);
300  constTrackRecAod->GetXYZ(xyz);
301  constTrackRecAod->GetPxPyPz(pxpypz);
302  constTrackRecAod->GetCovarianceXYZPxPyPz(cv);
303  AliExternalTrackParam par(xyz, pxpypz, cv, constTrackRecAod->Charge());
304  fhSigmaPtOverPtVsPt[itrkq]->Fill(constTrackRecAod->Pt(), TMath::Abs(sqrt(par.GetSigma1Pt2())/par.GetSigned1Pt()));
305  fMultCounter++;
306  }
307  fhMult->Fill(fMultCounter);
308  }
309  //_________________________________________________
310  // Background
311  Double_t rho = EstimateBgKT(jetContRecKT);
312  fhRho->Fill(rho);
313 
314  Double_t conerho = EstimateLocalBg(jetContRecAKT,trkContRec);
315  fhConeRho->Fill(conerho);
316 
317  Double_t cmsrho = EstimateBgKTCMS(jetContRecKT);
318  fhCMSRho->Fill(cmsrho);
319  //_________________________________________________
320  // Loop over AKT jets
321  if(jetContRecAKT){
322  jetContRecAKT->ResetCurrentID();
323  for(auto jetIterator : jetContRecAKT->accepted_momentum()) {//loop over reconstructed jets
324  jetRec=(AliEmcalJet*)jetIterator.second;
325  if(!jetRec) continue;
326  if(!IsSignalJetInAcceptance(jetRec)) continue; //check jet eta and pt
327  Double_t phi = Convert(jetRec->Phi());
328  fhJetPt->Fill(jetRec->Pt()); //fill AKT jet pT
329  fhJetPtRho->Fill(jetRec->Pt()-rho*jetRec->Area()); //fill AKT jet pT without bkg
330  fhJetPtConeRho->Fill(jetRec->Pt()-conerho*jetRec->Area()); //fill AKT jet pT without local bkg
331  fhJetPtCMSRho->Fill(jetRec->Pt()-cmsrho*jetRec->Area()); //fill AKT jet pT without CMS bkg
332  fhJetEtaPt->Fill(jetRec->Eta(),jetRec->Pt());
333  fhAktJetEtaPhi->Fill(jetRec->Eta(),phi);
334  fhJetPhiPt->Fill(phi,jetRec->Pt());
335  fhJetAreaPt->Fill(jetRec->Area(),jetRec->Pt());
336  fhAtimesRhoMedian->Fill(rho*jetRec->Area());
337  fhAtimesRhoCone->Fill(conerho*jetRec->Area());
338  fhAtimesRhoCMS->Fill(cmsrho*jetRec->Area());
339 
340  for(Int_t iq=0; iq < jetRec->GetNumberOfTracks(); iq++){
341  constTrackRec = static_cast<AliVParticle*> (jetRec->TrackAt(iq, trkContRec->GetArray())); //matched rec and emb tracks
342  if(!constTrackRec) continue;
343  fhJetConstituentPt->Fill(constTrackRec->Pt());
344 
345  }
346  }
347  }
348  // Loop over KT jets
349  if(jetContRecKT){
350  for(auto jetIterator : jetContRecKT->accepted_momentum()) {
351  jetRec=(AliEmcalJet*)jetIterator.second;
352  if(!jetRec) continue;
353  if(!IsSignalJetInAcceptance(jetRec)) continue;
354 
355  fhKTJetPt->Fill(jetRec->Pt());
356  }
357  }
358 
359  // Single particle efficiency and contamination
360  //primary particles spectrum
361  if(parContGen){
362  for(auto trackIterator : parContGen->accepted_momentum()) {
363  constTrackGen=(AliVParticle*)trackIterator.second;
364  if(!constTrackGen) continue;
365  if(IsTrackInAcceptance(constTrackGen, kTRUE)){
366  fhPrimGenTrkPt->Fill(constTrackGen->Pt()); //pT spectrum of generator level particles eta-pt, phi-pt, eta-phi histogram
367  fhGenTrackEtaPt->Fill(constTrackGen->Eta(),constTrackGen->Pt()); //pT spectrum of generator level particles eta-pt, phi-pt, eta-phi histogram
368  }
369  }
370  }
371 
372 
373  //reconstructed primary + secondary particles
374  Bool_t bRecPrim = kFALSE; //tags the reconstructed primary particles
375  if(trkContRec && parContGen){
376  for(auto trackIterator : trkContRec->accepted_momentum()) {
377  constTrackRecAod=(AliAODTrack*)trackIterator.second;
378  if(!constTrackRecAod->IsHybridGlobalConstrainedGlobal()) continue;
379  if(!constTrackRecAod) continue;
380  if(!IsTrackInAcceptance(constTrackRecAod, kFALSE)) continue; //reconstructed level tracks
381  bRecPrim = kFALSE; //reset matching flag
382 
383  parContGen->ResetCurrentID();
384  for(auto trackIterator : parContGen->accepted_momentum()) {
385  constTrackGen=(AliAODTrack*)trackIterator.second;
386  if(!constTrackGen) continue;
387  if(!IsTrackInAcceptance(constTrackGen, kTRUE)) continue; //gen level physical primary
388  if(TMath::Abs(constTrackRecAod->GetLabel()) == TMath::Abs(constTrackGen->GetLabel())){
389  fhRecTrkPt->Fill(constTrackGen->Pt());
390  bRecPrim = kTRUE; //matched
391  break;
392  }
393  }
394  if(!bRecPrim){
395  fhFakeTrkPt->Fill(constTrackRecAod->Pt());//Pt of fake tracks
396  }
397  }
398  }
399  // Response matrix
400  Double_t ptGenCorr; //GEN jet pt corrected for rho
401  Double_t ptRecCorr; //REC jet pt corrected for rho
402 
403  //Response matrix normalization - spectrum of all generator level jets in acceptance
404  if(jetContGenAKT){
405  for(auto jetIterator : jetContGenAKT->accepted_momentum()) {
406  jetGen=(AliEmcalJet*)jetIterator.second;
407  if(!jetGen) continue;
408  if(!IsSignalJetInAcceptance(jetGen,kTRUE)) continue; //cuts on eta, pT ,are
409  fhGenJetPt->Fill(jetGen->Pt()); //Pt spectrum of MC jets
410  }
411  }
412  //Find closest gen level+rec level jets
413  if(jetContRecAKT){
414  for(auto jetIterator : jetContRecAKT->accepted_momentum()) {
415  jetRec=(AliEmcalJet*)jetIterator.second;
416  if(!jetRec) continue;
417  if(!IsSignalJetInAcceptance(jetRec,kTRUE)) continue; //cuts on eta, pT ,area
418 
419  jetGen = 0x0;
420  jetGen = jetRec->ClosestJet();
421 
422  if(!jetGen){ //did not find matching generator level jet
423  continue;
424  }
425  fhRemx->Fill(jetRec->Pt(),jetGen->Pt());//response matrix
426  }
427  }
428 
429 
430  return kTRUE;
431 }
432 
433 //________________________________________________________________________
435  //Treminate
436  PostData(1, fOutput);
437 
438  // Mandatory
439  fOutput = dynamic_cast<AliEmcalList*> (GetOutputData(1)); // '1' refers to the output slot
440  if(!fOutput) {
441  printf("ERROR: Output list not available\n");
442  return;
443  }
444 }
445 
446 //________________________________________________________________________
448  // Destructor. Clean-up the output list, but not the histograms that are put inside
449  if(fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
450  delete fOutput;
451  }
452  delete fHelperClass;
453 
454 }
455 
459  Bool_t oldStatus = TH1::AddDirectoryStatus();
460  TH1::AddDirectory(kFALSE);
461  TString name;
462 
463  //Histograms
464  fhJetPt = new TH1D("fhJetPt","Pt spectrum of the AKT jets",320,-20.0,300.0);
465  fhJetPtRho = new TH1D("fhJetPtRho","Pt spectrum of the AKT jets - kt bg",320,-20.0,300.0);
466  fhJetPtConeRho = new TH1D("fhJetPtConeRho","p_{T} spectrum of the a-k_{T} jets - local kt bg",320,-20.0,300.0);
467  fhJetPtCMSRho = new TH1D("fhJetPtCMSRho","p_{T} spectrum of the a-k_{T} jets - CMS bg",320,-20.0,300.0);
468  fhAtimesRhoMedian = new TH1D("fhAtimesRhoMedian","a-k_{T} background times area",320,-20.0,300.0);
469  fhAtimesRhoCone = new TH1D("fhAtimesRhoCone","cone background times area",320,-20.0,300.0);
470  fhAtimesRhoCMS = new TH1D("fhAtimesRhoCMS","CMS a-k_{T} background times area",320,-20.0,300.0);
471  fhJetConstituentPt = new TH1D("fhJetConstituentPt","Pt spectrum of the a-k_{T} jet constituents",300,0.0,300.0);
472  fhTrackPt = new TH1D("fhTrackPt","p_{T} spectrum of tracks",300,0.0,300.0);
473  fhJetAreaPt = new TH2D("fhJetAreaPt","Jet area vs pt",50,0.,2.,300,0.0,300.0);
474  fhJetEtaPt = new TH2D("fhJetEtaPt","#eta-p_{T} distribution of jets",50,-fSignalJetEtaWindow,fSignalJetEtaWindow,300,0.0,300.0);
475  fhAktJetEtaPhi = new TH2D("fhAktJetEtaPhi","Eta-Phi distribution of a-k_{T} jets",50,-fSignalJetEtaWindow,fSignalJetEtaWindow,50,-TMath::Pi(),TMath::Pi());
476  fhKtJetEtaPhi = new TH2D("fhKtJetEtaPhi","#eta-#phi distribution of k_{T} jets",50,-fSignalJetEtaWindow,fSignalJetEtaWindow,50,-TMath::Pi(),TMath::Pi());
477  fhTrackEtaPhi = new TH2D("fhTrackEtaPhi","#eta-#phi distribution of tracks",50,-fTrackEtaWindow,fTrackEtaWindow,50,-TMath::Pi(),TMath::Pi());
478  fhTrackEtaPt = new TH2D("fhTrackEtaPt","#eta-p_{T} distribution of tracks",50,-fTrackEtaWindow,fTrackEtaWindow,300,0.0,300.0);
479  fhGenTrackEtaPt = new TH2D("fhGenTrackEtaPt","#eta-p_{T} distribution of tracks",50,-fTrackEtaWindow,fTrackEtaWindow,300,0.0,300.0);
480  fhJetPhiPt = new TH2D("fhJetPhiPt","#phi-p_{T} distribution of jets",50,-TMath::Pi(),TMath::Pi(),300,0.0,300.0);
481  fhTrackPhiPt = new TH2D("fhTrackPhiPt","#phi-p_{T} distribution of tracks",50,-TMath::Pi(),TMath::Pi(),300,0.0,300.0);
482  fhZVertex = new TH1D("fhZVertex","Z vertex",300,-15.0,15.0);
483  fhZVertexBC = new TH1D("fhZVertexBC","Z vertex before the cut",200,-50.0,50.0);
484  fhXVertex = new TH1D("fhXVertex","X vertex",300,-15.0,15.0);
485  fhYVertex = new TH1D("fhYVertex","Y vertex",300,-15.0,15.0);
486  fhCuts = new TH1D("fhCuts","Cuts statistics",3,0,3);
487  fhRho = new TH1D("fhRho","k_{T} jet background",500,0,50.0);
488  fhConeRho = new TH1D("fhConeRho","Local a-k_{T} jet background",500,0,50.0);
489  fhCMSRho = new TH1D("fhCMSRho","CMS-style a-k_{T} jet background",500,0,50.0);
490  fhKTJetPt = new TH1D("fhKTJetPt","KT jets p_{T} spectrum",300,0.0,300.0);
491  fhPrimGenTrkPt = new TH1D("fhPrimGenTrkPt","MC track p_{T}",300,0.0,300.0);
492  fhGenJetPt = new TH1D("fhGenJetPt","MC jet p_{T}",300,0.0,300.0);
493  fhRecTrkPt = new TH1D("fhRecTrkPt","Correctly rec. track p_{T}",300,0.0,300.0);
494  fhFakeTrkPt = new TH1D("fhFakeTrkPt","Fake track p_{T}",300,0.0,300.0);
495  fhRemx = new TH2D("fhRemx","Response matrix",300,0.0,300.0,300,0.0,300.0);
496  fhMult = new TH1D("fhMult","Reconstructed track multiplicity",300,0,300.0);
497 
498  fhTrackPhiCG = new TH2D("fhTrackPhiCG","Global tracks azimuth vs p_{T} ",300,0,300,50,-TMath::Pi(),TMath::Pi());
499  fhTrackPhiTPCG = new TH2D("fhTrackPhiTPCG","Complementary tracks azimuth vs p_{T)",300,0,300.0,50,-TMath::Pi(),TMath::Pi());
500  for(int i=0;i<2;i++){
501  name = Form("fhSigmaPtOverPtVsPt%i",i);
502  fhSigmaPtOverPtVsPt[i] = new TH2D(name.Data(),"#sigma p_{T}/p_{T} vs p_{T}",100,0.0,100.0,300,0.0,3.0);
503  fOutput->Add(fhSigmaPtOverPtVsPt[i]);
504  name = Form("fhInvPtQVsPhi%i",i);
505  fhInvPtQVsPhi[i] = new TH2D(name.Data(),"Negative track azimuth vs 1/pt",50,-TMath::Pi(),TMath::Pi(),300,0.0,0.25);
506  fOutput->Add(fhInvPtQVsPhi[i]);
507  name = Form("fhInvPtQVsEta%i",i);
508  fhInvPtQVsEta[i] = new TH2D(name.Data(),"Negative track azimuth vs 1/pt",50,-fTrackEtaWindow,fTrackEtaWindow,300,0.0,0.25);
509  fOutput->Add(fhInvPtQVsEta[i]);
510  name = Form("fhInvPtQVsPhiASide%i",i);
511  fhInvPtQVsPhiASide[i] = new TH2D(name.Data(),"Negative phi vs 1/pt A-eta",50,-TMath::Pi(),TMath::Pi(),300,0.0,0.25);
512  fOutput->Add(fhInvPtQVsPhiASide[i]);
513  name = Form("fhInvPtQVsPhiCSide%i",i);
514  fhInvPtQVsPhiCSide[i] = new TH2D(name.Data(),"Negative phi vs 1/pt C-eta",50,-TMath::Pi(),TMath::Pi(),300,0.0,0.25);
515  fOutput->Add(fhInvPtQVsPhiCSide[i]);
516  }
517 
518  fOutput->Add(fhJetPt);
519  fOutput->Add(fhJetPtRho);
520  fOutput->Add(fhJetPtConeRho);
521  fOutput->Add(fhJetPtCMSRho);
523  fOutput->Add(fhAtimesRhoCone);
524  fOutput->Add(fhAtimesRhoCMS);
525  fOutput->Add(fhJetAreaPt);
526 
528  fOutput->Add(fhTrackPt);
529  fOutput->Add(fhJetEtaPt);
530  fOutput->Add(fhTrackEtaPt);
531  fOutput->Add(fhJetPhiPt);
532 
533  fOutput->Add(fhTrackPhiPt);
534  fOutput->Add(fhZVertex);
535  fOutput->Add(fhXVertex);
536  fOutput->Add(fhYVertex);
537  fOutput->Add(fhCuts);
538 
539  fOutput->Add(fhRho);
540  fOutput->Add(fhConeRho);
541  fOutput->Add(fhCMSRho);
542  fOutput->Add(fhKTJetPt);
543  fOutput->Add(fhZVertexBC);
544  fOutput->Add(fhPrimGenTrkPt);
545 
546  fOutput->Add(fhGenJetPt);
547  fOutput->Add(fhRecTrkPt);
548  fOutput->Add(fhFakeTrkPt);
549  fOutput->Add(fhRemx);
550  fOutput->Add(fhMult);
551  fOutput->Add(fhGenTrackEtaPt);
552 
553  fOutput->Add(fhTrackPhiCG);
554  fOutput->Add(fhTrackPhiTPCG);
555  fOutput->Add(fhAktJetEtaPhi);
556  fOutput->Add(fhKtJetEtaPhi);
557  fOutput->Add(fhTrackEtaPhi);
558  // =========== Switch on Sumw2 for all histos ===========
559  for(Int_t i=0; i<fOutput->GetEntries(); i++){
560  TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
561  if(h1){
562  h1->Sumw2();
563  continue;
564  }
565  THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
566  if(hn){
567  hn->Sumw2();
568  }
569  }
570  TH1::AddDirectory(oldStatus);
571 
572 
573  PostData(1, fOutput);
574 }
575 //________________________________________________________________________
577  //
578  // retrieve event objects
579  //
581 
582  return kTRUE;
583 }
586 {
587 
588  return kTRUE;
589 }
590 
595  Double_t rhoCone = 0.0;
596 
597  if(!jetCont) return rhoCone;
598 
599  AliEmcalJet* jet = NULL;
600  static Double_t rhovec[999];
601  Int_t nJetAcc = 0;
602  Double_t jetpt;
603 
604  for(auto jetIterator : jetCont->accepted_momentum()) {
605  jet=(AliEmcalJet*)jetIterator.second;
606  if(!jet) continue;
607  if(!IsSignalJetInAcceptance(jet,kFALSE)) continue; //check jet eta and pt
608 
609  jetpt = jet->Pt();
610  if(jetpt <0.005) jetpt = 0.; //set pt of ghost jets identical to zero
611  rhovec[nJetAcc] = jetpt/jet->Area();
612  nJetAcc++;
613  }
614 
615  if(nJetAcc>0){
616  rhoCone = TMath::Median(nJetAcc, rhovec); //take the median of pTjet/Ajet
617  }
618 
619  return rhoCone;
620 }
626  Double_t rhoCone = 0.0;
627 
628  if(!jetCont) return rhoCone;
629 
630  AliEmcalJet* jet = NULL;
631  AliVParticle* track = NULL;
632  Double_t trackpt=0;
633  Double_t tracketa=0;
634  Double_t trackphi=0;
635  Double_t coner = 0.4;//the radius of the perpendicular cone
636  Double_t coner2 = coner*coner;
637  Double_t highjetpt=-1;
638  Double_t highjeteta=0;
639  Double_t highjetphi=0;
640  Double_t perpphi1=0;
641  Double_t perpphi2=0;
642  Double_t bgpt=0;
643 
644  Double_t X=0;
645  Double_t Y=0;
646  Double_t Z=0;
647 
648  for(auto jetIterator : jetCont->accepted_momentum()) {
649  jet=(AliEmcalJet*)jetIterator.second;
650  if(!jet) continue;
651  if(!IsSignalJetInAcceptance(jet)) continue; //check jet eta and pt
652  if(jet->Pt()>highjetpt){
653  highjetpt = jet->Pt();
654  highjeteta = jet->Eta();
655  highjetphi = Convert(jet->Phi());
656  }
657  }
658  if(highjetpt<0) return 0.;
659  perpphi1=Convert(highjetphi+TMath::Pi()/2);//the two perpendicular angles
660  perpphi2=Convert(highjetphi-TMath::Pi()/2);
661 
662  for(auto trackIterator : trkCont->accepted_momentum()) {
663  track=(AliVParticle*)trackIterator.second;
664  if(!((AliAODTrack*)track)->IsHybridGlobalConstrainedGlobal()) continue;
665  if(!track) continue;
666  if(!IsTrackInAcceptance(track)) continue;
667 
668  trackpt = track->Pt();
669  tracketa = track->Eta();
670  trackphi = Convert(track->Phi());
671  X=tracketa-highjeteta;
672  Y=Convert(trackphi-perpphi1);
673  Z=Convert(trackphi-perpphi2);
674 
675  if((X*X+Y*Y) < coner2 || (X*X+Z*Z) < coner2){
676  bgpt=bgpt+trackpt;
677  }
678  }
679  rhoCone = bgpt/(2*TMath::Pi()*coner2); //calculate the pT density
680 
681  return rhoCone;
682 }
687  Double_t rhoCMS = 0.0;
688 
689  if(!jetCont) return rhoCMS;
690 
691  AliEmcalJet* jet = NULL;
692  static Double_t rhovec[999];
693  Int_t nJetAcc = 0;
694  Double_t jetpt = 0.0;
695  Double_t physJetArea = 0.0;
696  Double_t allJetArea = 0.0;
697 
698  for(auto jetIterator : jetCont->accepted_momentum()) {
699  jet=(AliEmcalJet*)jetIterator.second;
700  if(!jet) continue;
701  if(!IsSignalJetInAcceptance(jet,kFALSE)) continue; //check jet eta and pt
702 
703  if(jet->Pt() > 0.01){//Select only physical jets
704  physJetArea += jet->Area(); //Area of physical jets
705  jetpt = jet->Pt();
706  rhovec[nJetAcc] = jetpt/jet->Area();
707  nJetAcc++;
708  }
709  allJetArea += jet->Area();//Area of all jets
710  }
711 
712  if(nJetAcc>0 && allJetArea!=0){
713  rhoCMS = TMath::Median(nJetAcc, rhovec)*physJetArea/allJetArea; //Calculate the CMS background
714  }
715 
716  return rhoCMS;
717 }
721 Double_t AliAnalysisTaskJetPP::Convert(Double_t input){//Converts an angle to the range (-Pi,Pi)
722  Double_t phi;
723  phi = TMath::ATan2(TMath::Sin(input),TMath::Cos(input));
724  return phi;
725 }
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
Bool_t IsSPDClusterVsTrackletBG(AliVEvent *event)
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
Float_t fBSPDCvsTCut
SPD cluster-vs-tracklet cut.
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
Float_t fASPDCvsTCut
SPD cluster-vs-tracklet cut.
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