TutProducer_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief TODO
3 /// \author bckhouse@hep.caltech.edu
4 ////////////////////////////////////////////////////////////////////////
5 
6 #include <string>
7 #include <vector>
8 
9 // Framework includes
16 
17 #include "Geometry/Geometry.h"
18 #include "GeometryObjects/Geo.h"
20 #include "MCCheater/BackTracker.h"
21 #include "RecoBase/Prong.h"
22 
23 // ROOT includes
24 #include "TVector3.h"
25 
26 
27 /// Tutorial modules
28 namespace tut
29 {
31  {
32  public:
33  explicit TutProducer(const fhicl::ParameterSet& pset);
34  ~TutProducer();
35 
36  void produce(art::Event& evt);
37 
38  void reconfigure(const fhicl::ParameterSet& pset);
39 
40  void beginJob();
41 
42  protected:
44 
45  rb::Prong FitPhoton(art::Handle<std::vector<rb::CellHit> > chits,
46  TVector3 vtx, int photId) const;
47  };
48 }
49 
50 
51 ////////////////////////////////////////////////////////////////////////
52 namespace tut
53 {
54  //.......................................................................
56  {
57  reconfigure(pset);
58 
59  produces<std::vector<rb::Prong> >();
60  }
61 
62  //......................................................................
64  {
65  }
66 
67  //......................................................................
69  {
70  fGeantLabel = pset.get<std::string>("GeantLabel");
71  fCellHitLabel = pset.get<std::string>("CellHitLabel");
72  }
73 
74  //......................................................................
76  {
77  }
78 
79  //......................................................................
81  {
82  // Declare a container for Prong objects to be stored in the art::event
83  std::unique_ptr<std::vector<rb::Prong> > prongcol(new std::vector<rb::Prong>);
84 
85  // This bit is all duplicated out of TutFilter. Need to find that pi0
86  // again. In a non-demo context we probably wouldn't bother with the filter
87  // at all.
88 
90  evt.getByLabel(fCellHitLabel, chits);
91 
92  // get the particle navigator from the BackTracker
94  const sim::ParticleNavigator& nav = bt->ParticleNavigator();
95 
96  // pizero should be a primary
97  for(int primIdx = 0; primIdx < nav.NumberOfPrimaries(); ++primIdx){
98  const sim::Particle* pizero = nav.Primary(primIdx);
99  // look for pi0's only
100  if(pizero->PdgCode() != 111) continue;
101 
102  if(pizero->NumberDaughters() != 2) continue;
103 
104  // Check the daughters are what they should be
105  const sim::Particle* phot0 = nav[pizero->Daughter(0)];
106  if(phot0->PdgCode() != 22) continue;
107  const sim::Particle* phot1 = nav[pizero->Daughter(1)];
108  if(phot1->PdgCode() != 22) continue;
109 
110  // New code starts here
111 
112  // True vertex
113  const TVector3 vtx = phot0->Position(0).Vect();
114 
115  // Stick photons in an array so we can loop over them
116  const sim::Particle* phots[2] = {phot0, phot1};
117 
118  for(int photIdx = 0; photIdx < 2; ++photIdx){
119  // This function is getting long, so do all the actual work in
120  // FitPhoton(). Put the result in the prongs collection
121  const int photId = phots[photIdx]->TrackId();
122  prongcol->push_back(FitPhoton(chits, vtx, photId));
123  } // end for photIdx
124  } // end for primIdx
125 
126  // Don't forget to save our prongs...
127  evt.put(std::move(prongcol));
128  }
129 
130  //......................................................................
131  rb::Prong TutProducer::FitPhoton(art::Handle<std::vector<rb::CellHit> > chits,
132  TVector3 vtx,
133  int photId) const
134  {
135  // Services we're going to need
138 
139  // We're going to store all the constituent hits here
140  art::PtrVector<rb::CellHit> prongHits;
141 
142  // Fit each view seperately and store the gradients in here
143  double dvdz[2];
144  for(geo::View_t view: {geo::kX, geo::kY}){
145 
146  // We're going to accumulate information about all the matching hits
147  // to do a straight line fit with.
148  std::vector<double> zs, vs, pes;
149 
150  // go through hits and cheat to find what hits are from what photon
151  for(unsigned int chitIdx = 0; chitIdx < chits->size(); ++chitIdx){
152  art::Ptr<rb::CellHit> chit(chits, chitIdx);
153 
154  if(bt->IsNoise(chit)) continue; // if hit is noise, ignore
155 
156  if(chit->View() != view) continue; // Wrong view
157 
158  // Did this hit make any major contribution to the prong?
159  bool any = false;
160 
161  // get track IDs of tracks that made light in this cell
162  const std::vector<cheat::TrackIDE> ides = bt->HitToTrackIDE(chit);
163 
164  // go through IDs
165  for(unsigned int idIdx = 0; idIdx < ides.size(); ++idIdx){
166  const int id = ides[idIdx].trackID;
167 
168  // Did this ID contribute to the right photon?
169  if(id == photId){
170  // good, so what fraction of the total energy deposited came
171  // from this particle?
172  const double frac = ides[idIdx].energyFrac;
173 
174  // If it's the majority contributor, keep it
175  if(frac > 0.5) any = true;
176 
177  // Now find out the physical position of the cell
178  double xyz[3];
179  geo::View_t junk; // Already know the view
180  geom->CellInfo(chit->Plane(), chit->Cell(), &junk, xyz);
181 
182  zs.push_back(xyz[2]);
183  // The view we actually measure with this cell
184  vs.push_back(xyz[view]);
185  // Don't know where the track goes exactly yet so can't calibrate
186  pes.push_back(frac*chit->PE());
187  } // end if
188  } // end for idIdx
189 
190  if(any) prongHits.push_back(chit);
191  } // end for chitIdx
192 
193  // Not enough hits to do a fit
194  if(vs.size() < 2){
195  dvdz[view] = 0;
196  continue;
197  }
198 
199  double z1, v1, z2, v2;
200  geo::LinFitMinDperp(zs, vs, pes, &z1, &v1, &z2, &v2);
201 
202  dvdz[view] = (v2-v1)/(z2-z1);
203  } // end for view
204 
205  // Combine two fits into one 3D vector
206  const TVector3 dir(dvdz[0], dvdz[1], 1);
207 
208  // For the sake of the demo, use the known true vertex
209  return rb::Prong(prongHits, vtx, dir);
210  }
211 
212 } // end namespace tut
213 
214 
215 ////////////////////////////////////////////////////////////////////////
216 namespace tut
217 {
219 }
back track the reconstruction to the simulation
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:218
int PdgCode() const
Definition: MCParticle.h:211
rb::Prong FitPhoton(art::Handle< std::vector< rb::CellHit > > chits, TVector3 vtx, int photId) const
std::vector< TrackIDE > HitToTrackIDE(const rb::CellHit &hit, bool useBirksE=false) const
Convenience function. HitsToTrackIDE but for a single hit.
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
Definition: BackTracker.h:744
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
unsigned short Plane() const
Definition: CellHit.h:39
geo::View_t View() const
Definition: CellHit.h:41
Vertical planes which measure X.
Definition: PlaneGeo.h:28
DEFINE_ART_MODULE(TestTMapFile)
int NumberDaughters() const
Definition: MCParticle.h:216
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
void CellInfo(unsigned int ip, unsigned int ic, View_t *view=0, double *pos=0, double *dpos=0) const
bool IsNoise(const art::Ptr< rb::CellHit > &hit) const
Is this hit not associated with any particles?
int TrackId() const
Definition: MCParticle.h:209
int Daughter(const int i) const
Definition: MCParticle.cxx:112
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
unsigned short Cell() const
Definition: CellHit.h:40
const sim::Particle * Primary(const int) const
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
double LinFitMinDperp(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &w, double *x1, double *y1, double *x2, double *y2)
Find the best-fit line to a collection of points in 2-D by minimizing the squared perpendicular dista...
Definition: Geo.cxx:449
T get(std::string const &key) const
Definition: ParameterSet.h:231
TODO.
float PE() const
Definition: CellHit.h:42
Collect Geo headers and supply basic geometry functions.
double frac(double x)
Fractional part.
std::string fCellHitLabel
void produce(art::Event &evt)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
A Cluster with defined start position and direction.
Definition: Prong.h:19
TDirectory * dir
Definition: macro.C:5
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
void geom(int which=0)
Definition: geom.C:163
TutProducer(const fhicl::ParameterSet &pset)
Encapsulate the geometry of one entire detector (near, far, ndos)
void reconfigure(const fhicl::ParameterSet &pset)