AliPhysics  608b256 (608b256)
AliAnalysisTaskParticleInJet.cxx
Go to the documentation of this file.
1 /*
2  * AliAnalysisTaskParticleInJet.cxx
3  *
4  * Created on: Feb 17, 2016
5  * Author: markus
6  */
7 #include <map>
8 #include <string>
9 
10 #include <TArrayD.h>
11 #include <THashList.h>
12 #include <THistManager.h>
13 #include <TMath.h>
14 #include <TString.h>
15 
16 #include "AliAODMCParticle.h"
17 #include "AliEmcalJet.h"
18 #include "AliEmcalTrackSelection.h"
19 #include "AliInputEventHandler.h"
20 #include "AliJetContainer.h"
21 #include "AliMCEvent.h"
22 #include "AliMCParticle.h"
23 #include "AliParticleContainer.h"
24 #include "AliVEvent.h"
25 #include "AliVParticle.h"
26 #include "AliVTrack.h"
27 
29 
31 
34 fHistMgr(NULL),
35 fTrackSelection(NULL),
36 fParticleContainerNameRec(""),
37 fParticleContainerNameMC(""),
38 fJetContainerNameRec(""),
39 fJetContainerNameMC("")
40 {
41 }
42 
45 fHistMgr(NULL),
46 fTrackSelection(NULL),
47 fParticleContainerNameRec(""),
48 fParticleContainerNameMC(""),
49 fJetContainerNameRec(""),
50 fJetContainerNameMC("")
51 {
52 }
53 
56 }
57 
59  fHistMgr = new THistManager("histos");
60 
61  TArrayD ptbinning, jetptbinning;
62  CreatePtBinning(ptbinning);
63  CreateLinearBinning(jetptbinning, 200, 0., 200);
64 
65  std::map<std::string, std::string> histmap1D, histmap2D;
66  histmap1D.insert(std::pair<std::string, std::string>("hMCall","MC true particles in full acceptance"));
67  histmap1D.insert(std::pair<std::string, std::string>("hMCcont","Accepted true particles in particle containers"));
68  histmap1D.insert(std::pair<std::string, std::string>("hTracksMB","All accepted tracks in the EMCAL acceptance in MB events"));
69  histmap1D.insert(std::pair<std::string, std::string>("hTracksEG1","All accepted tracks in the EMCAL acceptance in EG1 events"));
70  histmap1D.insert(std::pair<std::string, std::string>("hTracksEG2","All accepted tracks in the EMCAL acceptance in EG2 events"));
71  histmap1D.insert(std::pair<std::string, std::string>("hTracksEJ1","All accepted tracks in the EMCAL acceptance in EJ1 events"));
72  histmap1D.insert(std::pair<std::string, std::string>("hTracksEJ2","All accepted tracks in the EMCAL acceptance in EJ2 events"));
73  histmap2D.insert(std::pair<std::string, std::string>("hMCjetTrack","Particles in jets with jet pt"));
74  histmap2D.insert(std::pair<std::string, std::string>("hRecjetTrackMB","Particles in jets with jet pt in MB events"));
75  histmap2D.insert(std::pair<std::string, std::string>("hRecjetTrackEG1","Particles in jets with jet pt in EG1 events"));
76  histmap2D.insert(std::pair<std::string, std::string>("hRecjetTrackEG2","Particles in jets with jet pt in EG2 events"));
77  histmap2D.insert(std::pair<std::string, std::string>("hRecjetTrackEJ1","Particles in jets with jet pt in EJ1 events"));
78  histmap2D.insert(std::pair<std::string, std::string>("hRecjetTrackEJ2","Particles in jets with jet pt in EJ2 events"));
79 
80  for(std::map<std::string, std::string>::iterator it = histmap1D.begin(); it != histmap1D.end(); ++it){
81  fHistMgr->CreateTH1(it->first.c_str(), it->second.c_str(), ptbinning);
82  }
83  for(std::map<std::string, std::string>::iterator it = histmap2D.begin(); it != histmap2D.end(); ++it){
84  fHistMgr->CreateTH2(it->first.c_str(), it->second.c_str(), jetptbinning, ptbinning);
85  }
86 
87  fHistMgr->GetListOfHistograms()->SetOwner(kFALSE);
88  for(TIter histiter = TIter(fHistMgr->GetListOfHistograms()).Begin(); histiter != TIter::End(); ++histiter)
89  fHistosQA->Add(*histiter);
90 }
91 
93  TString triggerstring = fInputEvent->GetFiredTriggerClasses();
94  Bool_t isMinBias = fInputHandler->IsEventSelected() & AliVEvent::kINT7,
95  isEG1 = triggerstring.Contains("EG1"),
96  isEG2 = triggerstring.Contains("EG2"),
97  isEJ1 = triggerstring.Contains("EJ1"),
98  isEJ2 = triggerstring.Contains("EJ2");
99  if(!(isMinBias || isEG1 || isEG2 || isEJ1 || isEJ2)) return kFALSE;
100 
101  if(MCEvent()){
102  AliVParticle *mctrack = 0;
103  for(Int_t ipart = 0; ipart < MCEvent()->GetNumberOfTracks(); ipart++){
104  mctrack = MCEvent()->GetTrack(ipart);
105  if(!IsPhysicalPrimary(mctrack)) continue;
106  if(!AcceptParticle(mctrack)) continue;
107  fHistMgr->FillTH1("hMCall", mctrack->Pt());
108  }
109  std::vector<const AliVParticle *> particles = GetSelectedParticles(GetParticleContainer(fParticleContainerNameMC.Data()));
110  for(std::vector<const AliVParticle *>::iterator it = particles.begin(); it != particles.end(); ++it) {
111  const Double_t &pt = (*it)->Pt();
112  fHistMgr->FillTH1("hMCcont", pt);
113  }
114 
115  // Look at the MC Jet container
116  GetJetContainer(fJetContainerNameMC.Data())->ResetCurrentID();
117  AliEmcalJet *mcjet = GetJetContainer(fJetContainerNameMC.Data())->GetNextAcceptJet();
118  do{
119  for(int itrk = 0; itrk < mcjet->GetNumberOfTracks(); itrk++){
120  mctrack = mcjet->TrackAt(itrk, GetParticleContainer(0)->GetArray());
121  if(!AcceptParticle(mctrack)) continue;
122  fHistMgr->FillTH2("hMCjetTrack", mcjet->Pt(), mctrack->Pt());
123  }
124  } while ((mcjet = GetJetContainer(fJetContainerNameMC.Data())->GetNextAcceptJet()));
125  }
126 
127  // Loop over particles, select tracks, fill histograms of tracks in EMCAL
128  std::vector<const AliVParticle *> tracks = GetSelectedParticles(GetParticleContainer(fParticleContainerNameRec.Data()));
129  for(std::vector<const AliVParticle *>::iterator it = tracks.begin(); it != tracks.end(); ++it) {
130  const Double_t &pt = (*it)->Pt();
131  if(isMinBias) fHistMgr->FillTH1("hTracksMB", pt);
132  if(isEJ1) fHistMgr->FillTH1("hTracksEJ1", pt);
133  if(isEJ2) fHistMgr->FillTH1("hTracksEJ2", pt);
134  if(isEG1) fHistMgr->FillTH1("hTracksEG1", pt);
135  if(isEG2) fHistMgr->FillTH1("hTracksEG2", pt);
136  }
137 
138  // Look at the MC Jet container
139  AliVParticle *jettrack = NULL;
140  GetJetContainer(fJetContainerNameRec.Data())->ResetCurrentID();
141  AliEmcalJet *recjet = GetJetContainer(fJetContainerNameRec.Data())->GetNextAcceptJet();
142  do{
143  for(int itrk = 0; itrk < recjet->GetNumberOfTracks(); itrk++){
144  jettrack = recjet->TrackAt(itrk, GetParticleContainer(1)->GetArray());
145  if(!AcceptParticle(jettrack)) continue;
146  if(!AcceptTrack(dynamic_cast<AliVTrack *>(jettrack))) continue;
147  if(isMinBias) fHistMgr->FillTH2("hRecjetTrackMB", recjet->Pt(), jettrack->Pt());
148  if(isEJ1) fHistMgr->FillTH2("hRecjetTrackEJ1", recjet->Pt(), jettrack->Pt());
149  if(isEJ2) fHistMgr->FillTH2("hRecjetTrackEJ2", recjet->Pt(), jettrack->Pt());
150  if(isEG1) fHistMgr->FillTH2("hRecjetTrackEG1", recjet->Pt(), jettrack->Pt());
151  if(isEG2) fHistMgr->FillTH2("hRecjetTrackEG2", recjet->Pt(), jettrack->Pt());
152  }
153  } while ((recjet = GetJetContainer(fJetContainerNameRec.Data())->GetNextAcceptJet()));
154 
155  return true;
156 }
157 
158 std::vector<const AliVParticle *> AliAnalysisTaskParticleInJet::GetSelectedParticles(AliParticleContainer *const cont) const {
159  std::vector<const AliVParticle *> result;
160  cont->ResetCurrentID();
161  AliVParticle * test = cont->GetNextAcceptParticle();
162  AliVTrack *track = NULL;
163  do{
164  if(!AcceptParticle(test)) continue;
165  if((track = dynamic_cast<AliVTrack *>(test))){
166  // apply extra track cuts
167  if(!AcceptTrack(track)) continue;
168  }
169  result.push_back(test);
170  } while ((test = cont->GetNextAcceptParticle()));
171 
172  return result;
173 }
174 
175 Bool_t AliAnalysisTaskParticleInJet::AcceptParticle(const AliVParticle * const part) const{
176  if(part->Phi() < 1.4 || part->Phi() > 3.1) return false;
177  if(TMath::Abs(part->Eta()) > 0.5) return false;
178  if(!part->Charge()) return false;
179  return true;
180 }
181 
182 Bool_t AliAnalysisTaskParticleInJet::AcceptTrack(AliVTrack * const track) const{
183  if(!fTrackSelection->IsTrackAccepted(track)) return false;
184  return true;
185 }
186 
187 Bool_t AliAnalysisTaskParticleInJet::IsPhysicalPrimary(const AliVParticle *const part) const {
188  const AliMCParticle *mcpart = dynamic_cast<const AliMCParticle *>(part);
189  if(mcpart){
190  return MCEvent()->IsPhysicalPrimary(part->GetLabel());
191  } else {
192  const AliAODMCParticle *aodpart = dynamic_cast<const AliAODMCParticle *>(part);
193  if(aodpart) return aodpart->IsPhysicalPrimary();
194  }
195  return false;
196 }
197 
199  std::vector<double> mybinning;
200  std::map<double,double> definitions;
201  definitions.insert(std::pair<double, double>(1, 0.05));
202  definitions.insert(std::pair<double, double>(2, 0.1));
203  definitions.insert(std::pair<double, double>(4, 0.2));
204  definitions.insert(std::pair<double, double>(7, 0.5));
205  definitions.insert(std::pair<double, double>(16, 1));
206  definitions.insert(std::pair<double, double>(36, 2));
207  definitions.insert(std::pair<double, double>(40, 4));
208  definitions.insert(std::pair<double, double>(50, 5));
209  definitions.insert(std::pair<double, double>(100, 10));
210  definitions.insert(std::pair<double, double>(200, 20));
211  double currentval = 0.;
212  mybinning.push_back(currentval);
213  for(std::map<double,double>::iterator id = definitions.begin(); id != definitions.end(); ++id){
214  double limit = id->first, binwidth = id->second;
215  while(currentval < limit){
216  currentval += binwidth;
217  mybinning.push_back(currentval);
218  }
219  }
220  binning.Set(mybinning.size());
221  int ib = 0;
222  for(std::vector<double>::iterator it = mybinning.begin(); it != mybinning.end(); ++it)
223  binning[ib++] = *it;
224 }
225 
226 void AliAnalysisTaskParticleInJet::CreateLinearBinning(TArrayD& binning, int nbins, double min, double max) const {
227  double binwidth = (max-min)/static_cast<double>(nbins);
228  binning.Set(nbins+1);
229  binning[0] = min;
230  double currentlimit = min + binwidth;
231  for(int ibin = 0; ibin < nbins; ibin++){
232  binning[ibin+1] = currentlimit;
233  currentlimit += binwidth;
234  }
235 }
Bool_t IsPhysicalPrimary(const AliVParticle *const part) const
virtual AliVParticle * GetNextAcceptParticle()
double Double_t
Definition: External.C:58
AliJetContainer * GetJetContainer(Int_t i=0) const
std::vector< const AliVParticle * > GetSelectedParticles(AliParticleContainer *const cont) const
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
Fill a 2D histogram within the container.
Bool_t AcceptTrack(AliVTrack *const track) const
Bool_t AcceptParticle(const AliVParticle *const part) const
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
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
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.
int Int_t
Definition: External.C:63
THashList * GetListOfHistograms() const
Get the list of histograms.
Definition: THistManager.h:671
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.
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 CreatePtBinning(TArrayD &binning) const
AliEmcalTrackSelection * fTrackSelection
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
Fill a 1D histogram within the container.
Double_t Pt() const
Definition: AliEmcalJet.h:109
void CreateLinearBinning(TArrayD &binning, int nbins, double min, double max) const
virtual PWG::EMCAL::AliEmcalTrackSelResultPtr IsTrackAccepted(AliVTrack *const trk)=0
Interface for track selection code.
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
void test(int runnumber=195345)
const Int_t nbins
bool Bool_t
Definition: External.C:53
void ResetCurrentID(Int_t i=-1)
Reset the iterator to a given index.