AliPhysics  857879d (857879d)
 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 // ANALYSIS OF HIGH PT HADRON TRIGGER ASSOCIATED SPECTRUM OF RECOIL JETS IN P+PB
77 // Author Filip Krizek (7.Oct. 2015)
78 
80 //________________________________________________________________________________________
81 
84  fUseDefaultVertexCut(1), fUsePileUpCut(1), fIsMC(0),
85  fSignalJetRadius(0.4),
86 fSignalJetEtaWindow(0.9 - fSignalJetRadius), fTrackEtaWindow(0.9), fMinTrackPt(0.150), fMinJetArea(0.0),
87 fCentralityType("V0A"), fHelperClass(0), fInitializedLocal(0),
88 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),fhRho(0x0),fhConeRho(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),fhAtimesRho(0x0)
89 {
90  //default constructor
91 
92  //Initiate arrays= here in a cycle
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 
102 //________________________________________________________________________
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),fhRho(0x0),fhConeRho(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),fhAtimesRho(0x0)
110 
111 {
112 //Constructor
113  //Initiate arrays= here in a cycle
114  for(int i=0;i<2;i++){
115  fhInvPtQVsPhi[i] = NULL;
116  fhInvPtQVsEta[i] = NULL;
117  fhInvPtQVsPhiASide[i] = NULL;
118  fhInvPtQVsPhiCSide[i] = NULL;
119  fhSigmaPtOverPtVsPt[i] = NULL;
120  }
121 
122 
123  DefineOutput(1, AliEmcalList::Class());
124 }
125 
126 //________________________________________________________________________
128  //EVENT SELECTION RECONSTRUCTED DATA
129 
130  if(!event) return kFALSE;
131 
132  //___________________________________________________
133  //TEST PILE UP
134  if(fUsePileUpCut){
135  if(!fHelperClass || fHelperClass->IsPileUpEvent(event)){
136  return kFALSE;
137  }
138  if(fHelperClass && !fHelperClass->IsPileUpEvent(event)) fhCuts->Fill(1.5);//events that passed the pileup cut
139  }
140 
141  //___________________________________________________
142  //BEFORE VERTEX CUT
143 
144  fhZVertexBC->Fill(event->GetPrimaryVertex()->GetZ());
146  if(!fHelperClass || !fHelperClass->IsVertexSelected2013pA(event)){
147  return kFALSE;
148  }
149  if(fHelperClass && fHelperClass->IsVertexSelected2013pA(event)) fhCuts->Fill(2.5);//events that passed the vertex cut
150  }else{
151 
152  if(TMath::Abs(event->GetPrimaryVertex()->GetZ()) > fZVertexCut){ //DO VERTEX Z CUT
153  return kFALSE;
154  }
155  else fhCuts->Fill(2.5);//events that passed the vertex cut
156  }
157  fhZVertex->Fill(event->GetPrimaryVertex()->GetZ());
158  fhXVertex->Fill(event->GetPrimaryVertex()->GetX());
159  fhYVertex->Fill(event->GetPrimaryVertex()->GetY());
160 
161  //vertex z po cutu
162  //___________________________________________________
163  //AFTER VERTEX CUT
164  return kTRUE;
165 }
166 
167 //________________________________________________________________________
169  // CHECK THE TRACK PT AND ETA RANGE
170  if(!track) return kFALSE;
171 
172  if(isprimary) {
173  if(!track->Charge()) return kFALSE;
174  if(!(static_cast<AliAODMCParticle*>(track))->IsPhysicalPrimary()) return kFALSE;
175  }
176 
177  if(TMath::Abs(track->Eta()) <= fTrackEtaWindow){ //APPLY TRACK ETA CUT
178  if(track->Pt() >= fMinTrackPt){ //APPLY TRACK MINIMUM PT CUT
179  return kTRUE;
180  }
181  }
182  return kFALSE;
183 }
184 
185 //________________________________________________________________________
187  //select jets in acceptance
188  if(!jet) return kFALSE;
189  if(TMath::Abs(jet->Eta()) <= fSignalJetEtaWindow){
190  if(jet->Area() >= fMinJetArea){
191  if(suppressGhost){
192  if(jet->Pt() >= fMinTrackPt) return kTRUE;
193  }else{
194  return kTRUE;
195  }
196  }
197  }
198  return kFALSE;
199 }
200 
201 //________________________________________________________________________
203  // Initialization of jet containers done in AliAnalysisTaskEmcalJet::ExecOnce()
204  //Read arrays of jets and tracks
205  fInitializedLocal = kTRUE;
206 
207  // Initialize helper class (for vertex selection & pile up correction)
208  fHelperClass = new AliAnalysisUtils();
209  fHelperClass->SetCutOnZVertexSPD(kFALSE); // kFALSE: no cut; kTRUE: |zvtx-SPD - zvtx-TPC|<0.5cm
210 
211  return;
212 }
213 
214 
215 //________________________________________________________________________
217  // executed in each event
218  //called in AliAnalysisTaskEmcal::UserExec(Option_t *)
219  // Analyze the event and Fill histograms
220 
221  if(!InputEvent()){
222  AliError("??? Event pointer == 0 ???");
223  return kFALSE;
224  }
225 
226  fhCuts->Fill(0.5);
227  //Execute only once: Get tracks, jets from arrays if not already given
229 
230  //_________________________________________________________________
231  // FILL EVENT STATISTICS
232 
233  //Select events (vertex, pile-up,...)
234  if(!IsEventInAcceptance(InputEvent())) return kFALSE; //post data is in UserExec
235 
236  fhCuts->Fill(3.5);
237 
238  // END EVENT SELECTION
239 
240 
241  // JET+TRACK CONTAINERS
242  AliJetContainer *jetContRecAKT = NULL; //AKTjet container from reconstruced tracks
243  AliJetContainer *jetContRecKT = NULL; //KTjet container from reconstruced tracks
244  AliJetContainer *jetContGenAKT = NULL; //AKT jet container from generated MC tracks
245  AliJetContainer *jetContGenKT = NULL; //KT jet container from generated MC tracks
246 
247  AliEmcalJet *jetRec = NULL;
248  AliEmcalJet *jetGen = NULL;//MC
249 
250  AliParticleContainer *trkContRec = NULL; //track array of real reconstructed tracks
251  AliParticleContainer *parContGen = NULL; //track array of generated MC tracks
252 
253  AliVParticle *constTrackRec = NULL; //rec jet constituent
254  AliVParticle *constTrackGen = NULL; //gen jet constituent
255  //_________________________________________________________
256  //READ JET TRACK CONTAINERS
257  jetContRecAKT = GetJetContainer(kContainerOne); //AKT jet container
258  jetContRecKT = GetJetContainer(kContainerTwo); //KT jet container
259  if(fIsMC){
260  jetContGenAKT = GetJetContainer(kContainerThree); //MC AKT jet container
261  jetContGenKT = GetJetContainer(kContainerFour); //MC KT jet container
262  }
263 
264  trkContRec = GetParticleContainer(kContainerOne); //reconstructed hybrid tracks
265  if(fIsMC) parContGen = GetParticleContainer(kContainerTwo); //MC particles
266 
267  //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
268  //THIS IS LOOP OVER RECONSTRUCTED TRACKS
269  Int_t fMultCounter;
270  Double_t xyz[50];
271  Double_t pxpypz[50];
272  Double_t cv[21];
273  Int_t itrkq;
274 
275  AliAODTrack *constTrackRecAod=NULL;
276 
277  if(trkContRec){
278  fMultCounter=0;
279  trkContRec->ResetCurrentID();
280  while((constTrackRecAod = (AliAODTrack*) trkContRec->GetNextAcceptParticle())){
281  if(!constTrackRecAod->IsHybridGlobalConstrainedGlobal()) continue;
282  if(!constTrackRecAod) continue;
283  if(!IsTrackInAcceptance(constTrackRecAod)) continue;
284  Double_t phi = Convert(constTrackRecAod->Phi());
285  fhTrackPt->Fill(constTrackRecAod->Pt());
286  fhTrackEtaPt->Fill(constTrackRecAod->Eta(),constTrackRecAod->Pt());
287  fhTrackPhiPt->Fill(phi,constTrackRecAod->Pt());
288  fhTrackEtaPhi->Fill(constTrackRecAod->Eta(),phi);
289 
290  if(constTrackRecAod->IsGlobalConstrained()){
291  fhTrackPhiCG->Fill(constTrackRecAod->Pt(), phi); //global constrained
292  }else{
293  fhTrackPhiTPCG->Fill(constTrackRecAod->Pt(), phi); //complementary
294  }
295  itrkq = (constTrackRecAod->Charge()<0) ? 0 : 1;
296 
297  fhInvPtQVsPhi[itrkq]->Fill(phi, 1.0/constTrackRecAod->Pt());
298  fhInvPtQVsEta[itrkq]->Fill(constTrackRecAod->Eta(), 1.0/constTrackRecAod->Pt());
299 
300  if(constTrackRecAod->Eta()>0){
301  fhInvPtQVsPhiASide[itrkq]->Fill(phi, 1.0/constTrackRecAod->Pt());
302  }else{
303  fhInvPtQVsPhiCSide[itrkq]->Fill(phi, 1.0/constTrackRecAod->Pt());
304  }
305  //get sigma pT / pT
306  //Taken from AliEMCalTriggerExtraCuts::CalculateTPCTrackLength
307  memset(cv, 0, sizeof(Double_t) * 21); //cleanup arrays
308  memset(pxpypz, 0, sizeof(Double_t) * 50);
309  memset(xyz, 0, sizeof(Double_t) * 50);
310  constTrackRecAod->GetXYZ(xyz);
311  constTrackRecAod->GetPxPyPz(pxpypz);
312  constTrackRecAod->GetCovarianceXYZPxPyPz(cv);
313  AliExternalTrackParam par(xyz, pxpypz, cv, constTrackRecAod->Charge());
314  fhSigmaPtOverPtVsPt[itrkq]->Fill(constTrackRecAod->Pt(), TMath::Abs(sqrt(par.GetSigma1Pt2())/par.GetSigned1Pt()));
315  fMultCounter++;
316  }
317  fhMult->Fill(fMultCounter);
318  }
319  //_________________________________________________
320  //YOU CAN CALCULTE ENERGY DENSITY OF THE UNDERLYING EVENT (framework has for that also separate task, here i do it in my function)
321 
322  Double_t rho = EstimateBgKT(jetContRecKT);
323  fhRho->Fill(rho);
324 
325  Double_t conerho = EstimateLocalBg(jetContRecAKT,trkContRec);
326  fhConeRho->Fill(conerho);
327  //_________________________________________________
328  //LOOP OVER AKT JETS
329 
330  if(jetContRecAKT){
331  jetContRecAKT->ResetCurrentID();
332  while((jetRec = jetContRecAKT->GetNextAcceptJet())) { //loop over reconstructed jets
333  if(!jetRec) continue;
334  if(!IsSignalJetInAcceptance(jetRec)) continue; //check jet eta and pt
335  Double_t phi = Convert(jetRec->Phi());
336  fhJetPt->Fill(jetRec->Pt()); //fill AKT jet pT
337  fhJetPtRho->Fill(jetRec->Pt()-rho*jetRec->Area()); //fill AKT jet pT without bkg
338  fhJetPtConeRho->Fill(jetRec->Pt()-conerho*jetRec->Area()); //fill AKT jet pT without local bkg
339  fhJetEtaPt->Fill(jetRec->Eta(),jetRec->Pt());
340  fhAktJetEtaPhi->Fill(jetRec->Eta(),phi);
341  fhJetPhiPt->Fill(phi,jetRec->Pt());
342  fhJetAreaPt->Fill(jetRec->Area(),jetRec->Pt());
343  fhAtimesRho->Fill(rho*jetRec->Area());
344  //loop over jet constituents
345  for(Int_t iq=0; iq < jetRec->GetNumberOfTracks(); iq++){
346  constTrackRec = static_cast<AliVParticle*> (jetRec->TrackAt(iq, trkContRec->GetArray())); //matched rec and emb tracks
347  if(!constTrackRec) continue;
348  fhJetConstituentPt->Fill(constTrackRec->Pt());
349 
350  }
351  }
352  }
353 
355  if(jetContRecKT){
356  jetContRecKT->ResetCurrentID();
357  while((jetRec = jetContRecKT->GetNextAcceptJet())) {
358  if(!jetRec) continue;
359  if(!IsSignalJetInAcceptance(jetRec)) continue;
360 
361  fhKTJetPt->Fill(jetRec->Pt());
362  }
363  }
364 
365  //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
366  //SINGLE PARTICLE EFFICIENCY AND CONTAMINATION
367  //primary particles spectrum
368  if(parContGen){
369  parContGen->ResetCurrentID();
370  while((constTrackGen = (AliVParticle*)(parContGen->GetNextAcceptParticle()))){
371  if(!constTrackGen) continue;
372  if(IsTrackInAcceptance(constTrackGen, kTRUE)){
373  fhPrimGenTrkPt->Fill(constTrackGen->Pt()); //pT spectrum of generator level particles eta-pt, phi-pt, eta-phi histogram
374  fhGenTrackEtaPt->Fill(constTrackGen->Eta(),constTrackGen->Pt()); //pT spectrum of generator level particles eta-pt, phi-pt, eta-phi histogram
375  }
376  }
377  }
378 
379 
380  //reconstructed primary + secondary particles
381  Bool_t bRecPrim = kFALSE; //tags the reconstructed primary particles
382  if(trkContRec && parContGen){
383  trkContRec->ResetCurrentID();
384  while((constTrackRecAod =(AliAODTrack*) (trkContRec->GetNextAcceptParticle()))){
385  if(!constTrackRecAod->IsHybridGlobalConstrainedGlobal()) continue;
386  if(!constTrackRecAod) continue;
387  if(!IsTrackInAcceptance(constTrackRecAod, kFALSE)) continue; //reconstructed level tracks
388  bRecPrim = kFALSE; //reset matching flag
389 
390  parContGen->ResetCurrentID();
391  while((constTrackGen = (AliVParticle*)(parContGen->GetNextAcceptParticle()))){
392  if(!constTrackGen) continue;
393  if(!IsTrackInAcceptance(constTrackGen, kTRUE)) continue; //gen level physical primary
394  if(TMath::Abs(constTrackRecAod->GetLabel()) == TMath::Abs(constTrackGen->GetLabel())){
395  fhRecTrkPt->Fill(constTrackGen->Pt());
396  bRecPrim = kTRUE; //matched
397  break;
398  }
399  }
400  if(!bRecPrim){
401  fhFakeTrkPt->Fill(constTrackRecAod->Pt());//Pt of fake tracks
402  }
403  }
404  }
405  //+++++++++++++++++++++++++++++++++++++++++++++++++++++++
406  //2) FILL JET RESPONSE MATRIX
407  //+++++++++++++++++++++++++++++++++++++++++++++++++++++++
408  Double_t ptGenCorr; //GEN jet pt corrected for rho
409  Double_t ptRecCorr; //REC jet pt corrected for rho
410 
411  //Response matrix normalization - spectrum of all generator level jets in acceptance
412  if(jetContGenAKT){
413  jetContGenAKT->ResetCurrentID();
414  while((jetGen = jetContGenAKT->GetNextAcceptJet())){
415  if(!jetGen) continue;
416  if(!IsSignalJetInAcceptance(jetGen,kTRUE)) continue; //cuts on eta, pT ,are
417  fhGenJetPt->Fill(jetGen->Pt()); //Pt spectrum of MC jets
418  }
419  }
420  //Find closest gen level+rec level jets
421  if(jetContRecAKT){
422  jetContRecAKT->ResetCurrentID();
423  while((jetRec = jetContRecAKT->GetNextAcceptJet())) {
424  if(!jetRec) continue;
425  if(!IsSignalJetInAcceptance(jetRec,kTRUE)) continue; //cuts on eta, pT ,area
426 
427  jetGen = 0x0;
428  jetGen = jetRec->ClosestJet(); //z taggeru
429 
430  if(!jetGen){ //did not find matching generator level jet
431 
432  continue;
433  }
434  fhRemx->Fill(jetRec->Pt(),jetGen->Pt());//response matrix
435  }
436  }
437 
438 
439  return kTRUE;
440 }
441 
442 //________________________________________________________________________
444  //Treminate
445  PostData(1, fOutput);
446 
447  // Mandatory
448  fOutput = dynamic_cast<AliEmcalList*> (GetOutputData(1)); // '1' refers to the output slot
449  if(!fOutput) {
450  printf("ERROR: Output list not available\n");
451  return;
452  }
453 }
454 
455 //________________________________________________________________________
457  // Destructor. Clean-up the output list, but not the histograms that are put inside
458  // (the list is owner and will clean-up these histograms). Protect in PROOF case.
459  if(fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
460  delete fOutput;
461  }
462  delete fHelperClass;
463 
464 }
465 
466 //________________________________________________________________________
468  // called once to create user defined output objects like histograms, plots etc.
469  // and to put it on the output list.
470  // Note: Saving to file with e.g. OpenFile(0) is must be before creating other objects.
471  //fOutput TList defined in the mother class
473  Bool_t oldStatus = TH1::AddDirectoryStatus();
474  TH1::AddDirectory(kFALSE);
475  TString name;
476 
477  //Define your histograms here
478  fhJetPt = new TH1D("fhJetPt","Pt spectrum of the AKT jets",320,-20.0,300.0);
479  fhJetPtRho = new TH1D("fhJetPtRho","Pt spectrum of the AKT jets - bg",320,-20.0,300.0);
480  fhJetPtConeRho = new TH1D("fhJetPtConeRho","p_{T} spectrum of the a-k_{T} jets - bg",320,-20.0,300.0);
481  fhAtimesRho = new TH1D("fhAtimesRho","a-k_{T} background times area",320,-20.0,300.0);
482  fhJetConstituentPt = new TH1D("fhJetConstituentPt","Pt spectrum of the a-k_{T} jet constituents",300,0.0,300.0);
483  fhTrackPt = new TH1D("fhTrackPt","p_{T} spectrum of tracks",300,0.0,300.0);
484  fhJetAreaPt = new TH2D("fhJetAreaPt","Jet area vs pt",50,0.,2.,300,0.0,300.0);
485  fhJetEtaPt = new TH2D("fhJetEtaPt","#eta-p_{T} distribution of jets",50,-fSignalJetEtaWindow,fSignalJetEtaWindow,300,0.0,300.0);
486  fhAktJetEtaPhi = new TH2D("fhAktJetEtaPhi","Eta-Phi distribution of a-k_{T} jets",50,-fSignalJetEtaWindow,fSignalJetEtaWindow,50,-TMath::Pi(),TMath::Pi());
487  fhKtJetEtaPhi = new TH2D("fhKtJetEtaPhi","#eta-#phi distribution of k_{T} jets",50,-fSignalJetEtaWindow,fSignalJetEtaWindow,50,-TMath::Pi(),TMath::Pi());
488  fhTrackEtaPhi = new TH2D("fhTrackEtaPhi","#eta-#phi distribution of tracks",50,-fTrackEtaWindow,fTrackEtaWindow,50,-TMath::Pi(),TMath::Pi());
489  fhTrackEtaPt = new TH2D("fhTrackEtaPt","#eta-p_{T} distribution of tracks",50,-fTrackEtaWindow,fTrackEtaWindow,300,0.0,300.0);
490  fhGenTrackEtaPt = new TH2D("fhGenTrackEtaPt","#eta-p_{T} distribution of tracks",50,-fTrackEtaWindow,fTrackEtaWindow,300,0.0,300.0);
491  fhJetPhiPt = new TH2D("fhJetPhiPt","#phi-p_{T} distribution of jets",50,-TMath::Pi(),TMath::Pi(),300,0.0,300.0);
492  fhTrackPhiPt = new TH2D("fhTrackPhiPt","#phi-p_{T} distribution of tracks",50,-TMath::Pi(),TMath::Pi(),300,0.0,300.0);
493  fhZVertex = new TH1D("fhZVertex","Z vertex",300,-15.0,15.0);
494  fhZVertexBC = new TH1D("fhZVertexBC","Z vertex before the cut",200,-50.0,50.0);
495  fhXVertex = new TH1D("fhXVertex","X vertex",300,-15.0,15.0);
496  fhYVertex = new TH1D("fhYVertex","Y vertex",300,-15.0,15.0);
497  fhCuts = new TH1D("fhCuts","Cuts statistics",4,0,4);
498  fhRho = new TH1D("fhRho","k_{T} jet background",500,0,50.0);//histo rho*A
499  fhConeRho = new TH1D("fhConeRho","Local a-k_{T} jet background",500,0,50.0);
500  fhKTJetPt = new TH1D("fhKTJetPt","KT jets p_{T} spectrum",300,0.0,300.0);
501  fhPrimGenTrkPt = new TH1D("fhPrimGenTrkPt","MC track p_{T}",300,0.0,300.0);
502  fhGenJetPt = new TH1D("fhGenJetPt","MC jet p_{T}",300,0.0,300.0);
503  fhRecTrkPt = new TH1D("fhRecTrkPt","Correctly rec. track p_{T}",300,0.0,300.0);
504  fhFakeTrkPt = new TH1D("fhFakeTrkPt","Fake track p_{T}",300,0.0,300.0);
505  fhRemx = new TH2D("fhRemx","Response matrix",300,0.0,300.0,300,0.0,300.0);
506  fhMult = new TH1D("fhMult","Reconstructed track multiplicity",300,0,300.0);
507 
508  fhTrackPhiCG = new TH2D("fhTrackPhiCG","Global tracks azimuth vs p_{T} ",300,0,300,50,-TMath::Pi(),TMath::Pi());
509  fhTrackPhiTPCG = new TH2D("fhTrackPhiTPCG","Complementary tracks azimuth vs p_{T)",300,0,300.0,50,-TMath::Pi(),TMath::Pi());
510  for(int i=0;i<2;i++){
511  name = Form("fhSigmaPtOverPtVsPt%i",i);
512  fhSigmaPtOverPtVsPt[i] = new TH2D(name.Data(),"#sigma p_{T}/p_{T} vs p_{T}",100,0.0,100.0,300,0.0,3.0);
513  fOutput->Add((TH2D*)fhSigmaPtOverPtVsPt[i]);
514  name = Form("fhInvPtQVsPhi%i",i);
515  fhInvPtQVsPhi[i] = new TH2D(name.Data(),"Negative track azimuth vs 1/pt",50,-TMath::Pi(),TMath::Pi(),300,0.0,0.25);
516  fOutput->Add((TH2D*)fhInvPtQVsPhi[i]);
517  name = Form("fhInvPtQVsEta%i",i);
518  fhInvPtQVsEta[i] = new TH2D(name.Data(),"Negative track azimuth vs 1/pt",50,-fTrackEtaWindow,fTrackEtaWindow,300,0.0,0.25);
519  fOutput->Add((TH2D*)fhInvPtQVsEta[i]);
520  name = Form("fhInvPtQVsPhiASide%i",i);
521  fhInvPtQVsPhiASide[i] = new TH2D(name.Data(),"Negative phi vs 1/pt A-eta",50,-TMath::Pi(),TMath::Pi(),300,0.0,0.25);
522  fOutput->Add((TH2D*)fhInvPtQVsPhiASide[i]);
523  name = Form("fhInvPtQVsPhiCSide%i",i);
524  fhInvPtQVsPhiCSide[i] = new TH2D(name.Data(),"Negative phi vs 1/pt C-eta",50,-TMath::Pi(),TMath::Pi(),300,0.0,0.25);
525  fOutput->Add((TH2D*)fhInvPtQVsPhiCSide[i]);
526  }
527 
528  fOutput->Add((TH1D*)fhJetPt);
529  fOutput->Add((TH1D*)fhJetPtRho);
530  fOutput->Add((TH1D*)fhJetPtConeRho);
531  fOutput->Add((TH1D*)fhAtimesRho);
532  fOutput->Add((TH2D*)fhJetAreaPt);
533 
535  fOutput->Add((TH1D*)fhTrackPt);
536  fOutput->Add((TH2D*)fhJetEtaPt);
537  fOutput->Add((TH2D*)fhTrackEtaPt);
538  fOutput->Add((TH2D*)fhJetPhiPt);
539 
540  fOutput->Add((TH2D*)fhTrackPhiPt);
541  fOutput->Add((TH1D*)fhZVertex);
542  fOutput->Add((TH1D*)fhXVertex);
543  fOutput->Add((TH1D*)fhYVertex);
544  fOutput->Add((TH1D*)fhCuts);
545 
546  fOutput->Add((TH1D*)fhRho);
547  fOutput->Add((TH1D*)fhConeRho);
548  fOutput->Add((TH1D*)fhKTJetPt);
549  fOutput->Add((TH1D*)fhZVertexBC);
550  fOutput->Add((TH1D*)fhPrimGenTrkPt);
551 
552  fOutput->Add((TH1D*)fhGenJetPt);
553  fOutput->Add((TH1D*)fhRecTrkPt);
554  fOutput->Add((TH1D*)fhFakeTrkPt);
555  fOutput->Add((TH2D*)fhRemx);
556  fOutput->Add((TH1D*)fhMult);
557  fOutput->Add((TH2D*)fhGenTrackEtaPt);
558 
559  fOutput->Add((TH2D*)fhTrackPhiCG);
560  fOutput->Add((TH2D*)fhTrackPhiTPCG);
561  fOutput->Add((TH2D*)fhAktJetEtaPhi);
562  fOutput->Add((TH2D*)fhKtJetEtaPhi);
563  fOutput->Add((TH2D*)fhTrackEtaPhi);
564  // =========== Switch on Sumw2 for all histos ===========
565  for(Int_t i=0; i<fOutput->GetEntries(); i++){
566  TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
567  if(h1){
568  h1->Sumw2();
569  continue;
570  }
571  THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
572  if(hn){
573  hn->Sumw2();
574  }
575  }
576  TH1::AddDirectory(oldStatus);
577 
578 
579  PostData(1, fOutput);
580 }
581 //________________________________________________________________________
583  //
584  // retrieve event objects
585  //
587 
588  return kTRUE;
589 }
590 //________________________________________________________________________
592 {
593  // Run analysis code here, if needed. It will be executed before FillHistograms().
594 
595  return kTRUE;
596 }
597 
598 //________________________________________________________________________
599 
601  //Estimate rho by the standard area based approach (from KT jet median)
602  Double_t rhoCone = 0.0;
603 
604  if(!jetCont) return rhoCone;
605 
606  AliEmcalJet* jet = NULL;
607  static Double_t rhovec[999];
608  Int_t nJetAcc = 0;
609  Double_t jetpt;
610 
611  jetCont->ResetCurrentID();
612  while((jet = jetCont->GetNextAcceptJet())){ //loop over KT jets
613  if(!jet) continue;
614  if(!IsSignalJetInAcceptance(jet,kFALSE)) continue; //check jet eta and pt
615 
616  jetpt = jet->Pt();
617  if(jetpt <0.005) jetpt = 0.; //set pt of ghost jets identical to zero
618  rhovec[nJetAcc] = jetpt/jet->Area();
619  nJetAcc++;
620  }
621 
622  if(nJetAcc>0){
623  rhoCone = TMath::Median(nJetAcc, rhovec); //take the median of pTjet/Ajet
624  }
625 
626  return rhoCone;
627 }
628 //________________________________________________________________________
630  //Estimate background in a cone perpendicular to the leading jet
631  Double_t rhoCone = 0.0;
632 
633  if(!jetCont) return rhoCone;
634 
635  AliEmcalJet* jet = NULL;
636  AliVParticle* track = NULL;
637  Double_t trackpt=0;
638  Double_t tracketa=0;
639  Double_t trackphi=0;
640  Double_t coner = 0.4;//<<<<<<<<<<<<<<<<<the radius of the perpendicular cone
641  Double_t coner2 = coner*coner;
642  Double_t highjetpt=-1;
643  Double_t highjeteta=0;
644  Double_t highjetphi=0;
645  Double_t perpphi1=0;
646  Double_t perpphi2=0;
647  Double_t bgpt=0;
648 
649  Double_t X=0;
650  Double_t Y=0;
651  Double_t Z=0;
652 
653  jetCont->ResetCurrentID();
654  while((jet = jetCont->GetNextAcceptJet())){ //loop over KT jets
655  if(!jet) continue;
656  if(!IsSignalJetInAcceptance(jet)) continue; //check jet eta and pt
657  if(jet->Pt()>highjetpt){
658  highjetpt = jet->Pt();
659  highjeteta = jet->Eta();
660  highjetphi = Convert(jet->Phi());
661  }
662  }
663  if(highjetpt<0) return 0.;
664  perpphi1=Convert(highjetphi+TMath::Pi()/2);//the two perpendicular angles
665  perpphi2=Convert(highjetphi-TMath::Pi()/2);
666 
667  trkCont->ResetCurrentID();
668  while((track = trkCont->GetNextAcceptParticle())){
669  if(!((AliAODTrack*)track)->IsHybridGlobalConstrainedGlobal()) continue;
670  if(!track) continue;
671  if(!IsTrackInAcceptance(track)) continue;
672 
673  trackpt = track->Pt();
674  tracketa = track->Eta();
675  trackphi = track->Phi()-TMath::Pi();
676  X=tracketa-highjeteta;
677  Y=Convert(trackphi-perpphi1);
678  Z=Convert(trackphi-perpphi2);
679 
680  if(X*X+Y*Y<coner2 || X*X+Y*Y<coner2){
681  bgpt=bgpt+trackpt;
682  }
683  }
684  rhoCone = bgpt/(2*TMath::Pi()*coner2); //take the median of pTjet/Ajet
685 
686  return rhoCone;
687 }
688 Double_t AliAnalysisTaskJetPP::Convert(Double_t input){//Converts an angle to the range (-Pi,Pi)
689  return TMath::ATan2(TMath::Sin(input),TMath::Cos(input));
690 }
691 
Double_t EstimateBgKT(AliJetContainer *jetCont)
Bool_t IsSignalJetInAcceptance(AliEmcalJet *jet, Bool_t suppressGhost=1)
Double_t Area() const
Definition: AliEmcalJet.h:123
virtual AliVParticle * GetNextAcceptParticle()
double Double_t
Definition: External.C:58
Bool_t FillHistograms()
Function filling histograms.
TH1D * fhXVertex
Z vertex before cut.
TH2D * fhRemx
Track eta-phi distribution.
Bool_t IsTrackInAcceptance(AliVParticle *track, Bool_t isprimary=0)
TH2D * fhTrackEtaPhi
Track phi-pt distribution.
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:222
TH2D * fhGenTrackEtaPt
Track eta-pt distribution.
AliJetContainer * GetJetContainer(Int_t i=0) const
ClassImp(AliAnalysisTaskJetPP) AliAnalysisTaskJetPP
Double_t Eta() const
Definition: AliEmcalJet.h:114
TH1D * fhZVertex
Response matrix.
TH1D * fhJetPtConeRho
pt spectrum of AKT jets without kt bgk
Double_t Phi() const
Definition: AliEmcalJet.h:110
TH2D * fhInvPtQVsEta[2]
Inverse pt vs phi.
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.
TH2D * fhTrackPhiCG
Multiplicity.
TH2D * fhTrackPhiTPCG
Global tracks.
TH1D * fhRecTrkPt
Pt spectrum of MC jets.
TH1D * fhPrimGenTrkPt
Y vertex.
TH2D * fhTrackEtaPt
Jet et-phi distribution.
TH1D * fhMult
Pt spectrum of fake tracks.
TH2D * fhInvPtQVsPhiCSide[2]
Inverse pt vs phi in away-side rapidity.
TH1D * fhTrackPt
Jet constituent pt.
Double_t EstimateLocalBg(AliJetContainer *jetCont, AliParticleContainer *trkCont)
TH1D * fhJetConstituentPt
Histogram for pilup/vertex cuts.
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
TH1D * fhJetPt
gc trigger if tracks/jets are loaded initiates calling ExecOnce
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
TH2D * fhAktJetEtaPhi
Jet eta-pt distribution.
Double_t Convert(Double_t input)
int Int_t
Definition: External.C:63
TH2D * fhTrackPhiPt
Jet phi-pt distribution.
TH2D * fhJetAreaPt
local KT bgk
TH2D * fhKtJetEtaPhi
Jet et-phi distribution.
TH1D * fhJetPtRho
ptspectrum of KT jets
Definition: External.C:228
Definition: External.C:212
TH2D * fhInvPtQVsPhiASide[2]
Inverse pt vs eta.
TH2D * fhJetPhiPt
Generated track eta-pt distribution.
AliEmcalJet * GetNextAcceptJet()
TH1D * fhYVertex
X vertex.
Bool_t fInitializedLocal
gc Vertex selection helper
TH1D * fhKTJetPt
pt spectrum of AKT jets
TH2D * fhSigmaPtOverPtVsPt[2]
Inverse pt vs pi in close-side rapidity.
Double_t Pt() const
Definition: AliEmcalJet.h:102
Enhanced TList-derived class that implements correct merging for pt_hard binned production.
Definition: AliEmcalList.h:25
TH1D * fhCuts
Jet background times area.
TH2D * fhInvPtQVsPhi[2]
complementary tracks
AliEmcalList * fOutput
!output list
AliAnalysisUtils * fHelperClass
Bool_t RetrieveEventObjects()
Retrieve common objects from event.
TH1D * fhAtimesRho
pt spectrum of AKT jets without local bgk
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
TH1D * fhZVertexBC
Z vertex.
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 particle.
Container for jet within the EMCAL jet framework.
Definition: External.C:196
TH1D * fhFakeTrkPt
Pt spectrum of correctly reconstructed tracks.
TH2D * fhJetEtaPt
Jet Area-pt distribution.