All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AnalysisBaseDrawer.cxx
Go to the documentation of this file.
1 /// \file AnalysisBaseDrawer.cxx
2 /// \brief Class to aid in the rendering of AnalysisBase objects
3 /// \author msoderbe@syr.edu
4 
5 #include "larcore/Geometry/Geometry.h"
6 #include "lardataobj/AnalysisBase/Calorimetry.h"
7 #include "lardataobj/AnalysisBase/ParticleID.h"
8 #include "lardataobj/RecoBase/Track.h"
9 #include "lareventdisplay/EventDisplay/AnalysisBaseDrawer.h"
10 #include "lareventdisplay/EventDisplay/AnalysisDrawingOptions.h"
11 #include "lareventdisplay/EventDisplay/RecoDrawingOptions.h"
12 #include "lareventdisplay/EventDisplay/eventdisplay.h"
13 #include "nuevdb/EventDisplayBase/View2D.h"
14 
15 #include "art/Framework/Principal/Event.h"
16 #include "art/Framework/Principal/Handle.h"
17 #include "art/Framework/Services/Registry/ServiceHandle.h"
18 #include "canvas/Persistency/Common/FindMany.h"
19 #include "canvas/Persistency/Common/Ptr.h"
20 
21 #include "TLatex.h"
22 #include "TLine.h"
23 #include "TPolyMarker.h"
24 
25 #include <cmath>
26 
27 namespace evd{
28 
29  //......................................................................
30  AnalysisBaseDrawer::AnalysisBaseDrawer()
31  {
32 
33  }
34 
35  //......................................................................
36  AnalysisBaseDrawer::~AnalysisBaseDrawer()
37  {
38 
39  }
40 
41  //......................................................................
42  void AnalysisBaseDrawer::DrawDeDx(const art::Event& evt,
43  evdb::View2D* view)
44  {
45  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
46  art::ServiceHandle<evd::AnalysisDrawingOptions const> anaOpt;
47  art::ServiceHandle<geo::Geometry const> geom;
48 
49  for(size_t imod = 0; imod < recoOpt->fTrackLabels.size(); ++imod) {
50 
51 
52  //Get Track collection
53  art::InputTag which = recoOpt->fTrackLabels[imod];
54  art::Handle<std::vector<recob::Track> > trackListHandle;
55  evt.getByLabel(which,trackListHandle);
56  std::vector<art::Ptr<recob::Track> > tracklist;
57  art::fill_ptr_vector(tracklist, trackListHandle);
58 
59  //Loop over Calorimetry collections
60  for(size_t cmod = 0; cmod < anaOpt->fCalorimetryLabels.size(); ++cmod) {
61  std::string const callabel = anaOpt->fCalorimetryLabels[cmod];
62  //Association between Tracks and Calorimetry
63  art::FindMany<anab::Calorimetry> fmcal(trackListHandle, evt, callabel);
64  if (!fmcal.isValid()) continue;
65  //Loop over PID collections
66  for(size_t pmod = 0; pmod < anaOpt->fParticleIDLabels.size(); ++pmod) {
67  std::string const pidlabel = anaOpt->fParticleIDLabels[pmod];
68  //Association between Tracks and PID
69  art::FindMany<anab::ParticleID> fmpid(trackListHandle, evt, pidlabel);
70  if (!fmpid.isValid()) continue;
71 
72  //Loop over Tracks
73  int ntracks = 0;
74  for(size_t trkIter = 0; trkIter<tracklist.size(); ++trkIter){
75  if (anaOpt->fTrackID >=0 and tracklist[trkIter]->ID() != anaOpt->fTrackID) continue;
76  ++ntracks;
77  int color = tracklist[trkIter].key()%evd::kNCOLS;
78  std::vector<const anab::Calorimetry*> calos = fmcal.at(trkIter);
79  std::vector<const anab::ParticleID*> pids = fmpid.at(trkIter);
80  if (!calos.size()) continue;
81  if (calos.size()!=pids.size()) continue;
82  size_t bestplane = 0;
83  size_t calopl = 0;
84  //size_t pidpl = 0;
85  size_t nmaxhits = 0;
86  for (size_t icalo = 0; icalo < calos.size(); ++icalo){
87  if (calos[icalo]->dEdx().size() > nmaxhits){
88  nmaxhits = calos[icalo]->dEdx().size();
89  bestplane = calos[icalo]->PlaneID().Plane;
90  }
91  }
92  if (anaOpt->fCaloPlane>=0 and anaOpt->fCaloPlane<int(geom->Nplanes())){
93  for (size_t icalo = 0; icalo < calos.size(); ++icalo){
94  if (int(calos[icalo]->PlaneID().Plane)==anaOpt->fCaloPlane&&
95  calos[icalo]->dEdx().size())
96  bestplane = calos[icalo]->PlaneID().Plane;
97  }
98  }
99 
100  for (size_t icalo = 0; icalo < calos.size(); ++icalo){
101  if (calos[icalo]->PlaneID().Plane==bestplane){
102  calopl = icalo;
103  }
104  }
105  /*
106  for (size_t ipid = 0; ipid < pids.size(); ++ipid){
107  if (pids[ipid]->PlaneID().Plane==bestplane){
108  pidpl = ipid;
109  }
110  }
111  */
112 
113  TPolyMarker& pm = view->AddPolyMarker(calos[calopl]->dEdx().size(),evd::kColor[color],8,0.8);
114  for(size_t h = 0; h<calos[calopl]->dEdx().size();++h){
115  double xvalue = calos[calopl]->ResidualRange().at(h);
116  double yvalue = calos[calopl]->dEdx().at(h);
117  pm.SetPoint(h,xvalue,yvalue);
118 
119  double error = yvalue*(0.04231 + 0.0001783*(yvalue*yvalue));
120  TLine& l = view->AddLine(xvalue,yvalue-error,xvalue,yvalue+error);
121  l.SetLineColor(evd::kColor[color]);
122  }
123 
124  char trackinfo[80];
125  char pida[80];
126  char proton[80];
127  //char kaon[80];
128  char pion[80];
129  //char muon[80];
130  sprintf(trackinfo,"Track #%d: K.E. = %.1f MeV , Range = %.1f cm",
131  int(tracklist[trkIter].key()),
132  calos[calopl]->KineticEnergy(),
133  calos[calopl]->Range());
134  /*
135  sprintf(proton,"Proton Chi2 = %.1f, Kaon Chi2 = %.1f",
136  pids[pidpl]->Chi2Proton(),
137  pids[pidpl]->Chi2Kaon());
138 // sprintf(kaon,"Kaon Chi2 = %.1f",
139 // pids[pidpl]->Chi2Kaon());
140  sprintf(pion,"Pion Chi2 = %.1f, Muon Chi2 = %.1f",
141  pids[pidpl]->Chi2Pion(),
142  pids[pidpl]->Chi2Muon());
143 // sprintf(muon,"Muon Chi2 = %.1f",
144 // pids[pidpl]->Chi2Muon());
145  sprintf(pida,"Plane %d, PIDA = %.1f, NHits = %d",
146  calos[calopl]->PlaneID().Plane,
147  pids[pidpl]->PIDA(),
148  int(calos[calopl]->dEdx().size()));
149 
150  */
151  double offset = (ntracks - 1)*10.0;
152  TLatex& track_tex = view->AddLatex(13.0, (46.0) - offset,trackinfo);
153  TLatex& pida_tex = view->AddLatex(13.0, (46.0-2.5) - offset,pida);
154  TLatex& proton_tex = view->AddLatex(13.0, (46.0-5.0) - offset,proton);
155  //TLatex& kaon_tex = view->AddLatex(13.0, (46.0-4.0) - offset,kaon);
156  TLatex& pion_tex = view->AddLatex(13.0, (46.0-7.5) - offset,pion);
157  //TLatex& muon_tex = view->AddLatex(13.0, (46.0-8.0) - offset,muon);
158  track_tex.SetTextColor(evd::kColor[color]);
159  proton_tex.SetTextColor(evd::kColor[color]);
160  //kaon_tex.SetTextColor(evd::kColor[color]);
161  pion_tex.SetTextColor(evd::kColor[color]);
162  //muon_tex.SetTextColor(evd::kColor[color]);
163  pida_tex.SetTextColor(evd::kColor[color]);
164  track_tex.SetTextSize(0.05);
165  proton_tex.SetTextSize(0.05);
166  //kaon_tex.SetTextSize(0.05);
167  pion_tex.SetTextSize(0.05);
168  //muon_tex.SetTextSize(0.05);
169  pida_tex.SetTextSize(0.05);
170  }
171  }
172  }
173  }
174  }
175 
176  //......................................................................
177  void AnalysisBaseDrawer::DrawKineticEnergy(const art::Event& evt,
178  evdb::View2D* view)
179  {
180  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
181  art::ServiceHandle<evd::AnalysisDrawingOptions const> anaOpt;
182  art::ServiceHandle<geo::Geometry const> geom;
183  //add some legend-like labels with appropriate grayscale
184  char proton[80];
185  char kaon[80];
186  char pion[80];
187  char muon[80];
188  sprintf(proton,"proton");
189  sprintf(kaon,"kaon");
190  sprintf(pion,"pion");
191  sprintf(muon,"muon");
192  TLatex& proton_tex = view->AddLatex(2.0, 180.0,proton);
193  TLatex& kaon_tex = view->AddLatex(2.0, 165.0,kaon);
194  TLatex& pion_tex = view->AddLatex(2.0, 150.0,pion);
195  TLatex& muon_tex = view->AddLatex(2.0, 135.0,muon);
196  proton_tex.SetTextColor(kBlack);
197  kaon_tex.SetTextColor(kGray+2);
198  pion_tex.SetTextColor(kGray+1);
199  muon_tex.SetTextColor(kGray);
200  proton_tex.SetTextSize(0.075);
201  kaon_tex.SetTextSize(0.075);
202  pion_tex.SetTextSize(0.075);
203  muon_tex.SetTextSize(0.075);
204 
205  //now get the actual data
206  for(size_t imod = 0; imod < recoOpt->fTrackLabels.size(); ++imod) {
207  //Get Track collection
208  art::InputTag which = recoOpt->fTrackLabels[imod];
209  art::Handle<std::vector<recob::Track> > trackListHandle;
210  evt.getByLabel(which,trackListHandle);
211  std::vector<art::Ptr<recob::Track> > tracklist;
212  art::fill_ptr_vector(tracklist, trackListHandle);
213 
214  //Loop over Calorimetry collections
215  for(size_t cmod = 0; cmod < anaOpt->fCalorimetryLabels.size(); ++cmod) {
216  std::string const callabel = anaOpt->fCalorimetryLabels[cmod];
217  //Association between Tracks and Calorimetry
218  art::FindMany<anab::Calorimetry> fmcal(trackListHandle, evt, callabel);
219  if (!fmcal.isValid()) continue;
220 
221  //Loop over PID collections
222  for(size_t pmod = 0; pmod < anaOpt->fParticleIDLabels.size(); ++pmod) {
223  std::string const pidlabel = anaOpt->fParticleIDLabels[pmod];
224  //Association between Tracks and PID
225  art::FindMany<anab::ParticleID> fmpid(trackListHandle, evt, pidlabel);
226  if (!fmpid.isValid()) continue;
227 
228  //Loop over Tracks
229  for(size_t trkIter = 0; trkIter<tracklist.size(); ++trkIter){
230  if (anaOpt->fTrackID >=0 and tracklist[trkIter]->ID() != anaOpt->fTrackID) continue;
231  int color = tracklist[trkIter].key()%evd::kNCOLS;
232 
233  std::vector<const anab::Calorimetry*> calos = fmcal.at(trkIter);
234  if (!calos.size()) continue;
235  size_t bestplane = 0;
236  size_t nmaxhits = 0;
237  for (size_t icalo = 0; icalo < calos.size(); ++icalo){
238  if (calos[icalo]->dEdx().size() > nmaxhits){
239  nmaxhits = calos[icalo]->dEdx().size();
240  bestplane = icalo;
241  }
242  }
243  if (anaOpt->fCaloPlane>=0 and anaOpt->fCaloPlane<int(geom->Nplanes())){
244  for (size_t i = 0; i<geom->Nplanes(); ++i){
245  if (int(calos[i]->PlaneID().Plane)==anaOpt->fCaloPlane)
246  bestplane = i;
247  }
248  }
249 
250  double xvalue = calos[bestplane]->Range();
251  double yvalue = calos[bestplane]->KineticEnergy();
252  view->AddMarker(xvalue,yvalue,evd::kColor[color],8,0.8);
253  if(yvalue>0.0){
254  double error = yvalue*(0.6064/std::sqrt(yvalue));
255  TLine& l = view->AddLine(xvalue,yvalue-error,xvalue,yvalue+error);
256  l.SetLineColor(evd::kColor[color]);
257  }
258 
259  }
260  }
261  }
262  }
263 
264  }
265 
266  //......................................................................
267  void AnalysisBaseDrawer::CalorShower(const art::Event& evt,
268  evdb::View2D* view)
269  {
270  art::ServiceHandle<evd::RecoDrawingOptions const> recoOpt;
271  art::ServiceHandle<evd::AnalysisDrawingOptions const> anaOpt;
272 
273  for(size_t imod = 0; imod < recoOpt->fShowerLabels.size(); ++imod) {
274 
275  //Get Track collection
276  art::InputTag which = recoOpt->fShowerLabels[imod];
277  art::Handle<std::vector<anab::Calorimetry> > caloListHandle;
278  evt.getByLabel(which,caloListHandle);
279  std::vector<art::Ptr<anab::Calorimetry> > calolist;
280  art::fill_ptr_vector(calolist, caloListHandle);
281 
282 
283  //Loop over PID collections
284  for(size_t pmod = 0; pmod < anaOpt->fParticleIDLabels.size(); ++pmod) {
285  std::string const pidlabel = anaOpt->fParticleIDLabels[pmod];
286  //Association between Tracks and PID
287 
288  //Loop over Tracks
289  for(size_t shwIter = 0; shwIter<calolist.size(); ++shwIter){
290  int color = kRed;
291 
292  TPolyMarker& pm = view->AddPolyMarker((*calolist.at(shwIter)).dEdx().size(),color,8,0.8);
293  for(size_t h = 0; h<(*calolist.at(shwIter)).dEdx().size();++h){
294  pm.SetPoint(h,(*calolist.at(shwIter)).ResidualRange().at(h),(*calolist.at(shwIter)).dEdx().at(h));
295  }
296 
297 
298 
299  }
300  }
301  }
302 
303  char mip[80];
304  char mip2[80];
305 
306  sprintf(mip,"1 MIP");
307  sprintf(mip2,"2 MIP");
308  double offset = 0;
309 
310  double MIP = 2.12; // This is one mip in LAr, taken from uboone docdb #414
311  TLine & Line1Mip = view->AddLine(0, MIP, 100, MIP);
312  TLine & Line2Mip = view->AddLine(0, 2*MIP, 100, 2*MIP);
313 
314  TLatex& mip_tex = view->AddLatex(40.0, (23.0-20.0) - offset,mip);
315  TLatex& mip2_tex = view->AddLatex(40.0, (23.0-18.0) - offset,mip2);
316 
317  mip_tex.SetTextColor(kGray+3);
318  mip2_tex.SetTextColor(kGray+2);
319  mip_tex.SetTextSize(0.02);
320  mip2_tex.SetTextSize(0.02);
321 
322  Line1Mip.SetLineStyle(kDashed);
323  Line1Mip.SetLineColor(kGray+3);
324  Line2Mip.SetLineStyle(kDashed);
325  Line2Mip.SetLineColor(kGray+2);
326  }
327 
328 }// namespace
329 ////////////////////////////////////////////////////////////////////////
evd::AnalysisBaseDrawer::DrawDeDx
void DrawDeDx(const art::Event &evt, evdb::View2D *view)
Definition: AnalysisBaseDrawer.cxx:42
AnalysisBaseDrawer.h
Class to aid in the rendering of AnalysisBase objects.
offset
BEGIN_PROLOG TPC Trig offset(g4 rise time) ProjectToHeight
Definition: CORSIKAGen.fcl:7
Plane
BEGIN_PROLOG true icarus_rawdigitfilter FilterTools FilterPlane1 Plane
Definition: rawdigitfilter_icarus.fcl:8
test_msg.error
error
Definition: larcorealg/test/GeoAlgo/test_msg.py:32
icarus::ns::util::size
std::size_t size(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:561
evd::kNCOLS
static const int kNCOLS
Definition: eventdisplay.h:10
geo::fhicl::PlaneID
IDparameter< geo::PlaneID > PlaneID
Member type of validated geo::PlaneID parameter.
Definition: geo_types_fhicl.h:375
h
while getopts h
Definition: getStatistics.sh:13
and
return match has_match and(match.match_pdg==11 or match.match_pdg==-11)
tca::dEdx
float dEdx(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const TCSlice &slc, TP3D &tp3d)
Definition: PFPUtils.cxx:2687
evd::kColor
static const int kColor[kNCOLS]
Definition: eventdisplay.h:11
Track.h
Provides recob::Track data product.
eventdisplay.h
Place to keep constants for event display.
muon
process_name can override from command line with o or output muon
Definition: runPID.fcl:28
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
evd::AnalysisBaseDrawer::CalorShower
void CalorShower(const art::Event &evt, evdb::View2D *view)
Definition: AnalysisBaseDrawer.cxx:267
evd::AnalysisBaseDrawer::DrawKineticEnergy
void DrawKineticEnergy(const art::Event &evt, evdb::View2D *view)
Definition: AnalysisBaseDrawer.cxx:177
Geometry.h
art framework interface to geometry description
Generated by   doxygen 1.8.5