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