AliPhysics  9c66e61 (9c66e61)
AliAnalysisTaskTrackDensity.cxx
Go to the documentation of this file.
1 /*
2  * AliAnalysisTaskTrackDensity.cxx
3  *
4  * Created on: Mar 2, 2016
5  * Author: markus
6  */
7 #include <algorithm>
8 #include <THashList.h>
9 #include <THistManager.h>
10 #include <TMath.h>
11 
12 #include "AliAODMCParticle.h"
13 #include "AliEmcalJet.h"
14 #include "AliEmcalTrackSelection.h"
15 #include "AliParticleContainer.h"
16 #include "AliVEvent.h"
17 
19 
21 
22 namespace EMCalTriggerPtAnalysis {
23 
24 AliAnalysisTaskTrackDensity::AliAnalysisTaskTrackDensity() :
26  fHistos(NULL),
27  fTrackSelection(NULL),
28  fMCJetContainerName(""),
29  fMCParticleContainerName(""),
30  fJetRadii(),
31  fJetPtBins(),
32  fPtMinSteps(),
33  fParticlePtBinning()
34 {
35 }
36 
38  AliAnalysisTaskEmcalJet(name, kTRUE),
39  fHistos(NULL),
40  fTrackSelection(NULL),
43  fJetRadii(),
44  fJetPtBins(),
45  fPtMinSteps(),
47 {
49  double defaultJetRadii[] = {0, 0.01, 0.02, 0.05, 0.1, 0.2, 0.4},
50  defaultPtMinSteps[] = {0.5, 1, 2, 5, 10, 20},
51  defaultJetPtBins[] = {20, 40, 60, 80, 100, 150, 200, 1000},
52  defaultParticlePtBinning[] = {0., 0.5, 1., 1.5, 2., 3., 4., 5., 7.5, 10., 15., 20., 30., 40., 60., 80., 100.};
53  fJetRadii.Set(sizeof(defaultJetRadii)/sizeof(double), defaultJetRadii);
54  fJetPtBins.Set(sizeof(defaultJetPtBins)/sizeof(double), defaultJetPtBins);
55  fPtMinSteps.Set(sizeof(defaultPtMinSteps)/sizeof(double), defaultPtMinSteps);
56  fParticlePtBinning.Set(sizeof(defaultParticlePtBinning)/sizeof(double), defaultParticlePtBinning);
57 }
58 
61 }
62 
65  fHistos = new THistManager("trackdensityhistos");
67 
68  TArrayD linearBinning(102);
69  int counter = 0;
70  for(double val = -0.5; val <= 100.5; val += 1.){
71  linearBinning[counter++] = val;
72  }
73 
74  fHistos->CreateTH1("ParticlePt","p_{t} of particles; p_{t} (GeV/c); Yield", 100, 0., 100);
75  fHistos->CreateTH1("ParticlePtAcc","p_{t} of particles; p_{t} (GeV/c); Yield", 100, 0., 100);
76  fHistos->CreateTH2("ParticlePtJet","p_{t} of particles; p_{t, jet} (GeV/c); p_{t, particle} (GeV/c)", 200, 0., 400., 100, 0., 100);
77  fHistos->CreateTH2("ParticlePtAccJet","p_{t} of particles; p_{t, jet} (GeV/c); p_{t, particle} (GeV/c)", 200, 0., 400., 100, 0., 100);
78  fHistos->CreateTH2("JetMulitplicity","Number of contributors in jet; p_{t, jet} (GeV/c); Number of contributors", 200, 0., 400, 102, -0.5, 100.5);
79  for(int irad = 0; irad < fJetRadii.GetSize()-1; irad++){
80  for(int ptstep = 0; ptstep <fPtMinSteps.GetSize(); ptstep++){
81  fHistos->CreateTH2(Form("trackDensityJet_r%d_%d_minpt%d",
82  static_cast<int>(fJetRadii[irad] * 100.),
83  static_cast<int>(fJetRadii[irad+1] * 100.),
84  static_cast<int>(fPtMinSteps[ptstep] * 10.)),
85  Form("Density of tracks with p_{t} > %f GeV/c in r [%.2f, %.2f]; p_{t, jet} (GeV/c); Number of tracks",
86  fPtMinSteps[ptstep],
87  fJetRadii[irad],
88  fJetRadii[irad+1]),
89  200, 0., 200.,
90  102, -0.5, 100.5);
91  fHistos->CreateTH2(Form("trackDensityJetRec_r%d_%d_minpt%d",
92  static_cast<int>(fJetRadii[irad] * 100.),
93  static_cast<int>(fJetRadii[irad+1] * 100.),
94  static_cast<int>(fPtMinSteps[ptstep] * 10.)),
95  Form("Density of reconstructed tracks with p_{t} > %f GeV/c in r [%.2f, %.2f]; p_{t, jet} (GeV/c); Number of tracks",
96  fPtMinSteps[ptstep],
97  fJetRadii[irad],
98  fJetRadii[irad+1]),
99  200, 0., 200.,
100  102, -0.5, 100.5);
101  }
102  for(int jetptbin = 0 ; jetptbin < fJetPtBins.GetSize()-1; jetptbin++){
103  fHistos->CreateTH2(Form("trackDensityParticle_r%d_%d_jetpt%d_%d",
104  static_cast<int>(fJetRadii[irad] * 100.),
105  static_cast<int>(fJetRadii[irad+1] * 100.),
106  static_cast<int>(fJetPtBins[jetptbin]),
107  static_cast<int>(fJetPtBins[jetptbin+1])),
108  Form("Density of tracks in jet with p_{t} [%.1f, %.1f] in r[%.2f,%2f]",
109  fJetPtBins[jetptbin],
110  fJetPtBins[jetptbin+1],
111  fJetRadii[irad],
112  fJetRadii[irad+1]),
113  fParticlePtBinning, linearBinning);
114  fHistos->CreateTH2(Form("trackDensityParticleRec_r%d_%d_jetpt%d_%d",
115  static_cast<int>(fJetRadii[irad] * 100.),
116  static_cast<int>(fJetRadii[irad+1] * 100.),
117  static_cast<int>(fJetPtBins[jetptbin]),
118  static_cast<int>(fJetPtBins[jetptbin+1])),
119  Form("Density of reconstructed tracks in jet with p_{t} [%.1f, %.1f] in r[%.2f,%2f]",
120  fJetPtBins[jetptbin],
121  fJetPtBins[jetptbin+1],
122  fJetRadii[irad],
123  fJetRadii[irad+1]),
124  fParticlePtBinning, linearBinning);
125  }
126  }
127 
128  for(TIter histiter = TIter(fHistos->GetListOfHistograms()).Begin(); histiter != TIter::End(); ++histiter){
129  fOutput->Add(*histiter);
130  }
131  PostData(1, fOutput);
132 }
133 
135  std::vector<int> recoAcceptLabels;
136  GetAcceptLabels(*(InputEvent()), recoAcceptLabels);
137  AliEmcalJet *myjet = NULL;
138  AliAODMCParticle *jetparticle = NULL;
139  AliJetContainer *mcjetcontainer = GetJetContainer(fMCJetContainerName.Data());
141  if(!mcparticleContainer) printf("Error getting particle container with name %s\n", fMCParticleContainerName.Data());
142  mcparticleContainer->ResetCurrentID(-1);
143  while((jetparticle = static_cast<AliAODMCParticle *>(mcparticleContainer->GetNextAcceptParticle()))){
144  if(TMath::Abs(jetparticle->Eta()) > 0.8) continue;
145  if(!jetparticle->Charge()) continue;
146  if(!jetparticle->IsPhysicalPrimary()) continue;
147  fHistos->FillTH1("ParticlePt", jetparticle->Pt());
148  if(std::binary_search(recoAcceptLabels.begin(), recoAcceptLabels.end(), TMath::Abs(jetparticle->GetLabel())))
149  fHistos->FillTH1("ParticlePtAcc", jetparticle->Pt());
150  }
151  double jetptmin, jetptmax;
152  mcjetcontainer->ResetCurrentID(-1);
153  while((myjet = mcjetcontainer->GetNextAcceptJet())){
154  fHistos->FillTH2("JetMulitplicity", myjet->Pt(), myjet->GetNumberOfConstituents());
155  for(int iconst = 0; iconst < myjet->GetNumberOfTracks(); iconst++){
156  jetparticle = static_cast<AliAODMCParticle *>(myjet->TrackAt(iconst, mcparticleContainer->GetArray()));
157  if(TMath::Abs(jetparticle->Eta()) > 0.8) continue;
158  if(!jetparticle->Charge()) continue;
159  if(!jetparticle->IsPhysicalPrimary()) continue;
160  fHistos->FillTH2("ParticlePtJet", myjet->Pt(), jetparticle->Pt());
161  if(std::binary_search(recoAcceptLabels.begin(), recoAcceptLabels.end(), TMath::Abs(jetparticle->GetLabel())))
162  fHistos->FillTH2("ParticlePtAccJet", myjet->Pt(), jetparticle->Pt());
163  }
164  for(int irad = 0 ; irad < fJetRadii.GetSize()-1; irad++){
165  for(int ptstep = 0; ptstep < fPtMinSteps.GetSize(); ptstep++){
166  fHistos->FillTH2(Form("trackDensityJet_r%d_%d_minpt%d", static_cast<int>(fJetRadii[irad] * 100.), static_cast<int>(fJetRadii[irad+1] * 100.), static_cast<int>(fPtMinSteps[ptstep] * 10.)),
167  TMath::Abs(myjet->Pt()), GetParticleMultiplicity(*myjet, *mcparticleContainer, fPtMinSteps[ptstep], 10000., fJetRadii[irad], fJetRadii[irad+1]));
168  fHistos->FillTH2(Form("trackDensityJetRec_r%d_%d_minpt%d", static_cast<int>(fJetRadii[irad] * 100.), static_cast<int>(fJetRadii[irad+1] * 100.), static_cast<int>(fPtMinSteps[ptstep] * 10.)),
169  TMath::Abs(myjet->Pt()), GetParticleMultiplicity(*myjet, *mcparticleContainer, fPtMinSteps[ptstep], 10000., fJetRadii[irad], fJetRadii[irad+1], &recoAcceptLabels));
170  }
171  FindJetPtBin(myjet, jetptmin, jetptmax);
172  if(jetptmin > 0 && jetptmax > 0){
173  for(int ptstep = 0; ptstep < fParticlePtBinning.GetSize()-1; ptstep++){
174  double mean = (fParticlePtBinning[ptstep] + fParticlePtBinning[ptstep+1])/2.;
175  fHistos->FillTH2(Form("trackDensityParticle_r%d_%d_jetpt%d_%d", static_cast<int>(fJetRadii[irad] * 100.),static_cast<int>(fJetRadii[irad+1] * 100.),
176  static_cast<int>(jetptmin), static_cast<int>(jetptmax)), mean, GetParticleMultiplicity(*myjet, *mcparticleContainer, fParticlePtBinning[ptstep], fParticlePtBinning[+1], fJetRadii[irad], fJetRadii[irad+1]));
177  fHistos->FillTH2(Form("trackDensityParticleRec_r%d_%d_jetpt%d_%d", static_cast<int>(fJetRadii[irad] * 100.),static_cast<int>(fJetRadii[irad+1] * 100.),
178  static_cast<int>(jetptmin), static_cast<int>(jetptmax)), mean, GetParticleMultiplicity(*myjet, *mcparticleContainer, fParticlePtBinning[ptstep], fParticlePtBinning[+1], fJetRadii[irad], fJetRadii[irad+1], &recoAcceptLabels));
179  }
180  }
181  }
182  }
183 
184  return true;
185 }
186 
187 int AliAnalysisTaskTrackDensity::GetParticleMultiplicity(const AliEmcalJet &jet, const AliParticleContainer &partcont, double ptmin, double ptmax, double rmin, double rmax, const std::vector<int> *labels) const {
188  AliDebug(1, Form("Next jet: %s\n", jet.toString().Data()));
189  TLorentzVector jetaxis(jet.Px(), jet.Py(), jet.Pz(), jet.E());
190  int nselected = 0;
191  const AliAODMCParticle *jetparticle(NULL);
192  for(int ipart = 0; ipart < jet.GetNumberOfTracks(); ipart++){
193  jetparticle = static_cast<AliAODMCParticle *>(jet.TrackAt(ipart, partcont.GetArray()));
194  if(TMath::Abs(jetparticle->Eta()) > 0.8) continue;
195  if(!TMath::Abs(jetparticle->Charge())) continue; // select charged particles only
196  if(labels){
197  if(!std::binary_search(labels->begin(), labels->end(), TMath::Abs(jetparticle->GetLabel()))) continue;
198  }
199  double partpt = TMath::Abs(jetparticle->Pt());
200  if(partpt >= ptmin && partpt < ptmax){
201  TLorentzVector partvector(jetparticle->Px(), jetparticle->Py(), jetparticle->Pz(), jetparticle->E());
202  double r = TMath::Abs(jetaxis.DeltaR(partvector));
203  if(r >= rmin && r < rmax) nselected++;
204  }
205  }
206  return nselected;
207 }
208 
209 void AliAnalysisTaskTrackDensity::FindJetPtBin(const AliEmcalJet *const jet, double &ptmin, double &ptmax) const {
210  ptmin = ptmax = -1;
211  double jetpt = TMath::Abs(jet->Pt());
212  for(int ptstep = 0; ptstep < fJetPtBins.GetSize() - 1; ptstep++){
213  if(jetpt >= fJetPtBins[ptstep] && jetpt < fJetPtBins[ptstep+1]){
214  ptmin = fJetPtBins[ptstep];
215  ptmax = fJetPtBins[ptstep+1];
216  break;
217  }
218  }
219 }
220 
221 void AliAnalysisTaskTrackDensity::GetAcceptLabels(const AliVEvent &event, std::vector<int> &labels) const {
222  AliVTrack * track = NULL;
223  for(int itrk = 0; itrk < event.GetNumberOfTracks(); ++itrk){
224  if(fTrackSelection->IsTrackAccepted((track = static_cast<AliVTrack *>(event.GetTrack(itrk))))) labels.push_back(TMath::Abs(track->GetLabel()));
225  }
226  std::sort(labels.begin(), labels.end());
227 }
228 
229 } /* namespace EMCalTriggerPtAnalysis */
AliJetContainer * GetJetContainer(Int_t i=0) const
Double_t Py() const
Definition: AliEmcalJet.h:107
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
Fill a 2D histogram within the container.
TString fMCParticleContainerName
Name of the MC jet container.
Double_t E() const
Definition: AliEmcalJet.h:119
void FindJetPtBin(const AliEmcalJet *const jet, double &ptmin, double &ptmax) const
UShort_t GetNumberOfConstituents() const
Definition: AliEmcalJet.h:140
Container for particles within the EMCAL framework.
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:160
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:139
Double_t Px() const
Definition: AliEmcalJet.h:106
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
int GetParticleMultiplicity(const AliEmcalJet &jet, const AliParticleContainer &partcont, double ptmin, double ptmax, double rmin, double rmax, const std::vector< int > *labels=NULL) const
TH2 * CreateTH2(const char *name, const char *title, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax, Option_t *opt="")
Create a new TH2 within the container.
THashList * GetListOfHistograms() const
Get the list of histograms.
Definition: THistManager.h:671
const Double_t ptmax
TString toString() const
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
Create a new TH1 within the container.
void ReleaseOwner()
Definition: THistManager.h:254
const Double_t ptmin
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
Fill a 1D histogram within the container.
AliEmcalJet * GetNextAcceptJet()
Double_t Pt() const
Definition: AliEmcalJet.h:109
virtual PWG::EMCAL::AliEmcalTrackSelResultPtr IsTrackAccepted(AliVTrack *const trk)=0
Interface for track selection code.
void GetAcceptLabels(const AliVEvent &event, std::vector< int > &ref) const
AliEmcalList * fOutput
!output list
Analysis of high- tracks in triggered events.
void SetMakeGeneralHistograms(Bool_t g)
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
Container class for histograms.
Definition: THistManager.h:99
TArrayD fJetRadii
Name of the MC particle container.
Double_t Pz() const
Definition: AliEmcalJet.h:108
virtual bool Run()
Run function. This is the core function of the analysis and contains the user code. Therefore users have to implement this function.
void UserCreateOutputObjects()
Main initialization function on the worker.
Container for jet within the EMCAL jet framework.