AliPhysics  f05a842 (f05a842)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
37 AliAnalysisTaskTrackDensity::AliAnalysisTaskTrackDensity(const char *name) :
38  AliAnalysisTaskEmcalJet(name, kTRUE),
39  fHistos(NULL),
40  fTrackSelection(NULL),
41  fMCJetContainerName(""),
42  fMCParticleContainerName(""),
43  fJetRadii(),
44  fJetPtBins(),
45  fPtMinSteps(),
46  fParticlePtBinning()
47 {
48  SetMakeGeneralHistograms(true);
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 
59 AliAnalysisTaskTrackDensity::~AliAnalysisTaskTrackDensity() {
60  if(fTrackSelection) delete fTrackSelection;
61 }
62 
63 void AliAnalysisTaskTrackDensity::UserCreateOutputObjects() {
65  fHistos = new THistManager("trackdensityhistos");
66  fHistos->ReleaseOwner();
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 
134 bool AliAnalysisTaskTrackDensity::Run(){
135  std::vector<int> recoAcceptLabels;
136  GetAcceptLabels(*(InputEvent()), recoAcceptLabels);
137  AliEmcalJet *myjet = NULL;
138  AliAODMCParticle *jetparticle = NULL;
139  AliJetContainer *mcjetcontainer = GetJetContainer(fMCJetContainerName.Data());
140  AliParticleContainer *mcparticleContainer = GetParticleContainer(fMCParticleContainerName.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 */
ClassImp(EMCalTriggerPtAnalysis::AliAnalysisTaskTrackDensity) namespace EMCalTriggerPtAnalysis
Double_t Py() const
Definition: AliEmcalJet.h:91
Double_t E() const
Definition: AliEmcalJet.h:103
UShort_t GetNumberOfConstituents() const
Definition: AliEmcalJet.h:124
Container for particles within the EMCAL framework.
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:123
Double_t Px() const
Definition: AliEmcalJet.h:90
const Double_t ptmax
TString toString() const
const Double_t ptmin
AliEmcalJet * GetNextAcceptJet()
Double_t Pt() const
Definition: AliEmcalJet.h:93
Short_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:144
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
Container class for histograms for the high- charged particle analysis.
Definition: THistManager.h:43
Double_t Pz() const
Definition: AliEmcalJet.h:92
Container for jet within the EMCAL jet framework.