DimuonFitter_module.cc
Go to the documentation of this file.
1 ///
2 /// \file DimuonFitter_module.cc
3 /// \brief Run the break-point fitter
4 /// \author messier@indiana.edu
5 /// \version $Id:$
6 ///
7 #include <algorithm>
8 //
9 #include "TH2F.h"
10 #include "TNtuple.h"
11 //
16 //
17 #include "Geometry/Geometry.h"
19 #include "RecoBase/Vertex.h"
20 #include "RecoBase/Cluster.h"
21 #include "RecoBase/Prong.h"
22 #include "RecoBase/Track.h"
23 #include "RecoBase/FitSum.h"
24 //
25 #include "Utilities/AssociationUtil.h"
26 //
33 
34 namespace bpfit { class DimuonFitter; }
35 
37 {
38 public:
40 
41 public:
42  // Member data
43  int fRun; ///< Slice identifiers (run, subrun, event, slice)
44  int fSubrun; ///< are useful to have as class scope for logging
45  int fEvent; ///< debugging information in each of the methods
46  int fSlice; ///< used during the reconstruction.
47  int fVerbosity; ///< 0-100, how much text output to generate
48  TNtuple* fPass1Nt; ///< pass one ntuple
49 public:
50  //......................................................................
51 
53  {
54  fVerbosity = p.get<int>("Verbosity");
55  }
56 
57  //......................................................................
58 
59  explicit DimuonFitter(fhicl::ParameterSet const& p) :
60  fSliceToken(consumes<std::vector<rb::Cluster>>(p.get<std::string>("SlicerLabel"))),
61  fRun(0),
62  fSubrun(0),
63  fEvent(0),
64  fSlice(0)
65  {
66  this->produces< std::vector<rb::Vertex> >();
67  this->produces< std::vector<rb::Track> >();
68  this->reconfigure(p);
69  }
70 
71  //......................................................................
72 
73  virtual ~DimuonFitter() { }
74 
75  void beginJob(){
77  fPass1Nt = f->make<TNtuple>("p1",
78  "pass 1 ntuple",
79  "run:event:slice:"
80  "np:nhitx:nhity:"
81  "np5:"
82  "nc5l:nc5r:"
83  "nc5t:nc5b");
84 
85  }
86  //......................................................................
87 
88  void FindSlices(const art::Event& evt,
89  art::Handle< std::vector<rb::Cluster> >& h,
91  {
92  unsigned int i;
93  evt.getByToken(fSliceToken, h);
94  for (i=0; i<h->size(); ++i) {
96  }
97  if (fVerbosity>10) {
98  std::cout << __FILE__ << ":" << __LINE__
99  << " Found " << s.size() << " slices."
100  << std::endl;
101  }
102  }
103 
104  //......................................................................
105  //
106  // Make an analysis of the variables corresponding to a single slice.
107  //
108  // Return "true" if this is a slice we want to keep in the analysis
109  // "false" if we thinkg this slice is ~100% likely to be background.
110  //
111  // This selection is reduction pass 1 for the analysis
112  //
113  //
114  bool AnaSlice(const rb::Cluster* s) {
115  unsigned int i;
116  //
117  // Define the cell boundaries at the edges of the detector
118  //
119  int icl = 5; // 5 cells on the left side
120  int icr = 90; // 5 cells on the right side
121  int icb = 5; // 5 cells on the bottom face
122  int ict = 90; // 5 cells on the top face
123  int ipf = 5; // 5 planes on the front face
124  int nc5l = 0; // number hits in the boundary on the left
125  int nc5r = 0; // "" on the right
126  int nc5b = 0; // "" on the bottom
127  int nc5t = 0; // "" on the top
128  int npl = 0; // Total number of planes in cluster
129  int np5 = 0; // Number of hits in the first 5 planes
130  std::vector<int> hitp(1000);
131  for (i=0; i<s->NCell(); ++i) {
132  int p = s->Cell(i)->Plane();
133  int c = s->Cell(i)->Cell();
134  //
135  // Count the number of hit planes
136  //
137  if (hitp[s->Cell(i)->Plane()]==0) ++npl;
138  hitp[s->Cell(i)->Plane()] = 1;
139  //
140  // Count cell hits near walls. Don't double count hits in the
141  // corners (hence the "else")
142  //
143  if (p<ipf) {
144  ++np5;
145  }
146  else {
147  if (s->Cell(i)->View()==geo::kX) {
148  if (c<icl) ++nc5r;
149  if (c>icr) ++nc5l;
150  }
151  else if (s->Cell(i)->View()==geo::kY) {
152  if (c<icb) ++nc5b;
153  if (c>ict) ++nc5t;
154  }
155  }
156  }
157  float x[16];
158  x[0] = fRun;
159  x[1] = fEvent;
160  x[2] = fSlice;
161  x[3] = npl;
162  x[4] = s->NXCell();
163  x[5] = s->NYCell();
164  x[6] = np5;
165  x[7] = nc5l;
166  x[8] = nc5r;
167  x[9] = nc5b;
168  x[10] = nc5t;
169  fPass1Nt->Fill(x);
170 
171  return true;
172  }
173  //......................................................................
174 
175  double FindVertexZ(const rb::Cluster* c) {
176  if (fVerbosity>20) {
177  std::cout << __FILE__ << ":" << __LINE__
178  << " Start 'FindVertexZ'. "
179  << " Cluster has " << c->NCell() << " cells" << std::endl;
180  }
181  unsigned int i, j;
182  unsigned int ip, ic;
183  //
184  // The way we are going to find a candidate vertex z location is
185  // to identify the first plane which is followed by a reasonably
186  // continuous set of activity. So, we will loop over all the hits,
187  // record their plane locations. Then we will look for the first
188  // plane in the list where if I step n positions forward in the
189  // list the plane number doesn't change more than m. While we are
190  // making this list, keep track of the lowest z position we've
191  // seen. In cases where we don't have enough hits to apply the
192  // algorithm above, we'll just use that lowest z value.
193  //
194  std::vector<unsigned int> p(c->NCell());
195  unsigned int pvtx = 999999;
196  for (i=0; i<c->NCell(); ++i) {
197  ip = c->Cell(i)->Plane();
198  if (ip<pvtx) pvtx = ip;
199  p[i] = ip;
200  }
201  if (fVerbosity>20) {
202  std::cout << __FILE__ << ":" << __LINE__
203  << " Set default vertex plane to pvtx="
204  << pvtx << std::endl;
205  }
206  std::sort(p.begin(), p.end());
207 
208  const unsigned int N = 2; // How far to skip forward in list
209  const unsigned int M = N+2; // Plane numbers should be this close or closer
210  for (i=0, j=N; j+1<c->NCell(); ++i,++j) {
211  if ( (p[j]-p[i])<=M ) {
212  pvtx = p[i];
213  if (fVerbosity>20) {
214  std::cout << __FILE__ ":" << __LINE__
215  << " Plane " << pvtx << " satisfies "
216  << " N= " << N
217  << " M= " << M
218  << " pi = " << p[i] << " condition."
219  << " pj = " << p[j] << std::endl;
220  }
221  break;
222  }
223  }
224  //
225  // Now we have to match the plane number to a z location.
226  //
228  for (i=0; i<c->NCell(); ++i) {
229  ip = c->Cell(i)->Plane();
230  if (ip==pvtx) {
231  ic = c->Cell(i)->Cell();
232  double xyz[3];
233  geom->Plane(ip)->Cell(ic)->GetCenter(xyz);
234  if (fVerbosity>10) {
235  std::cout << __FILE__ << ":" << __LINE__
236  << " Found vertex z to z="
237  << xyz[2] << " cm." << std::endl;
238  }
239  return xyz[2];
240  }
241  }
242  //
243  // Shouldn't ever get here.
244  //
245  abort();
246  return 0.;
247  }
248 
249  //......................................................................
250 
251  void FitView(const rb::Cluster& clust, geo::View_t view, double zvtx)
252  {
253  if (fVerbosity>20) {
254  std::cout << __FILE__ << ":" << __LINE__
255  << " Starting FitView for view " << view
256  << std::endl;
257  }
258  unsigned int i;
259  const double zero_weight = 1e-9;
260  //
261  // My initial weighting scheme depends on the location of the hits
262  // inside the slice so first I need to work out the bounding box
263  // that contains the cluster
264  //
265  unsigned int plo = 999999;
266  unsigned int phi = 0;
267  unsigned int clo = 999999;
268  unsigned int chi = 0;
269  for (i=0; i<clust.NCell(); ++i) {
270  const rb::CellHit* h = clust.Cell(i).get();
271  unsigned int ip = h->Plane();
272  unsigned int ic = h->Cell();
273  if (ip<plo) plo = ip;
274  if (ip>phi) phi = ip;
275  if (ic<clo) clo = ic;
276  if (ic>chi) chi = ic;
277  }
278  if (fVerbosity>20) {
279  std::cout << __FILE__ << ":" << __LINE__
280  << " Found slice extent "
281  << " [" << plo << "," << phi << "] "
282  << " [" << clo << "," << chi << "]"
283  << std::endl;
284  }
285 
286  //
287  // Initialize the data for the fit. Need x,z locations, error bars
288  // on x and weights for each of the hitss
289  //
290  std::vector<double> x1; // Transverse hit locations (cm)
291  std::vector<double> z1; // Longitudinal hit locations (cm)
292  std::vector<double> s1; // Transverse uncertainty hit locations (cm)
293  std::vector<double> w1; // Weight of the hit to be associated with track 1
294 
295  std::vector<double> x2; // Transverse hit locations (cm)
296  std::vector<double> z2; // Longitudinal hit locations (cm)
297  std::vector<double> s2; // Transverse uncertainty hit locations (cm)
298  std::vector<double> w2; // Weight of the hit to be associated with track 1
299 
300  double zlo = 999999;
301  double zhi = -999999;
302  double rt12 = sqrt(12.0);
303  for (i=0; i<clust.NCell(); ++i) {
304  const rb::CellHit* h = clust.Cell(i).get();
305  unsigned int ip = h->Plane();
306  unsigned int ic = h->Cell();
307 
309  const geo::PlaneGeo* p = geom->Plane(ip);
310  const geo::CellGeo* c = p->Cell(ic);
311  geo::View_t v = p->View();
312  //
313  // If the hit doesn't match the view we've been asked to fit,
314  // skip
315  //
316  if (v != view) continue;
317 
318  double xyz[3];
319  c->GetCenter(xyz, clust.W(h));
320  if (xyz[2]<zlo) zlo = xyz[2];
321  if (xyz[2]>zhi) zhi = xyz[2];
322  double sigxy = c->HalfW()/rt12;
323 
324 
325  if (v==geo::kX) { x1.push_back(xyz[0]); x2.push_back(xyz[0]); }
326  else if (v==geo::kY) { x1.push_back(xyz[1]); x2.push_back(xyz[1]); }
327  else { abort(); }
328  z1.push_back(xyz[2]);
329  z2.push_back(xyz[2]);
330  s1.push_back(sigxy);
331  s2.push_back(sigxy);
332  //
333  // We have to break the symmetry in the weights somehow or the
334  // two tracks will converge to the same solution. I'm splitting
335  // according to cell number to that the "top" track generally
336  // will be track one and the "bottom" track generally will be
337  // track 2
338  //
339  double wp = (float)(phi-ip)/(float)(phi-plo);
340  double wc1 = (float)(chi-ic)/(float)(chi-clo);
341  double wc2 = (float)(ic-clo)/(float)(chi-clo);
342  wc1 = (wp*wp)*(wc1*wc1);
343  wc2 = (wp*wp)*(wc2*wc2);
344  wc1 = wc1/((wc1+wc2)+zero_weight);
345  wc2 = 1.0-wc1;
346  wc1 += zero_weight;
347  wc2 += zero_weight;
348  w1.push_back(wc1);
349  w2.push_back(wc2);
350  }
351  if (fVerbosity>30) {
352  for (i=0; i<z1.size(); ++i) {
353  std::cout << __FILE__ << ":" << __LINE__ << "\t"
354  << i << "\t"
355  << z1[i] << "\t" << z2[i] << "\t"
356  << x1[i] << "\t" << x2[i] << "\t"
357  << s1[i] << "\t" << s2[i] << "\t"
358  << w1[i] << "\t" << w2[i] << std::endl;
359  }
360  }
361  if (fVerbosity>20) {
362  std::cout << __FILE__ << ":" << __LINE__
363  << " Initialized fits."
364  << std::endl;
365  }
366 
367  //
368  // Construct the scattering surfaces. Do something simple to get
369  // started: A straight line down the middle of the
370  // detector. Should be good enough to estimate the total radiation
371  // lengths the tracks encounter.
372  //
373  int pdg_muon = 13;
374  double ds = 50.0;
375  double dx = 1.0;
376  double dt = 0.3;
377  double p1[3] = {3,3,zlo};
378  double p2[3] = {3,3,zhi};
379  std::cout << __FILE__ << ":" << __LINE__ << std::endl;
380  std::vector<double> s(z1);
382  bpfit::ScatteringSurfaces scatter(pdg_muon, ds, dx, dt);
383  sort(s.begin(), s.end());
384  scatter.MakeSurfaces(*geom, p1, p2, s);
385  if (fVerbosity>30) {
386  std::cout << __FILE__ << ":" << __LINE__ << " Made "
387  << scatter.N() << " scattering surfaces between z="
388  << zlo << ", " << zhi << std::endl;
389  }
390  //
391  // Now perform the fits
392  //
393  double lambda = 1/25.;
394  for (unsigned int itr=0; itr<5; ++itr) {
395  if (fVerbosity>30) {
396  std::cout << __FILE__ << ":" << __LINE__
397  << " Fitting iteration " << itr
398  << std::endl;
399  }
400  bpfit::Lutz lutz1(z1.size(), scatter.N());
401  bpfit::Lutz lutz2(z2.size(), scatter.N());
402  for (i=0; i<z1.size(); ++i) {
403  lutz1.SetMeasurement(i, z1[i], x1[i], s1[i]/w1[i]);
404  }
405  for (i=0; i<z2.size(); ++i) {
406  lutz2.SetMeasurement(i, z2[i], x2[i], s2[i]/w2[i]);
407  }
408  for (i=0; i<scatter.N(); ++i) {
409  lutz1.SetScatteringPlane(i, scatter.S(i), scatter.SigAlpha(i));
410  lutz2.SetScatteringPlane(i, scatter.S(i), scatter.SigAlpha(i));
411  }
412 
413  double a1, b1, chi1;
414  double a2, b2, chi2;
415  lutz1.Fit(&a1, &b1, 0, &chi1);
416  lutz2.Fit(&a2, &b2, 0, &chi2);
417  if (fVerbosity>20) {
418  std::cout << __FILE__ << ":" << __LINE__
419  << " Finished fits in view " << view << "."
420  << " chi1 = " << chi1
421  << " chi2 = " << chi2
422  << std::endl;
423  }
424  //
425  // Update the weights for the next iteration
426  //
427  lambda = (itr+1)*(itr+1)/25.0;
428  for (i=0; i<z1.size(); ++i) {
429  double ex1 = exp(-lutz1.Chi2XIi(i));
430  double ex2 = exp(-lutz2.Chi2XIi(i));
431  double ex3 = exp(-25.0*lambda);
432  w1[i] = ex1/(ex3+ex1+ex2);
433  w2[i] = ex2/(ex3+ex1+ex2);
434  if (w1[i]<zero_weight) w1[i] = zero_weight;
435  if (w2[i]<zero_weight) w2[i] = zero_weight;
436  }
437  }
438  }
439 
440  //......................................................................
441 
443  {
444  //
445  // Allocations for the track products and their associations
446  //
447  std::unique_ptr< std::vector<rb::Vertex> >
448  vtx(new std::vector<rb::Vertex>);
449 
450  fRun = evt.run();
451  fSubrun = evt.subRun();
452  fEvent = evt.id().event();
453  fSlice = 0;
454  //
455  // Extract the slices, vertices, and tracks we need as inputs
456  //
459  this->FindSlices(evt, sliceh, slicep);
460 
461  unsigned int islice;
462  for (islice=0; islice<slicep.size(); ++islice) {
463  //
464  // Get the index of this slices and a pointer to it for
465  // convenient access. Skip the noice slice, no physics there for
466  // us.
467  //
468  fSlice = islice;
469  const art::Ptr<rb::Cluster> slice = slicep[islice];
470  if (slice->IsNoise()) continue;
471  //
472  // Analyze the slice and see if it passes selections for analysis
473  //
474  this->AnaSlice(slice.get());
475 
476  //
477  // The track fits require a candidate z vertex location
478  //
479  double zvtx = this->FindVertexZ(slice.get());
480  vtx->push_back(rb::Vertex(0,0,zvtx,0));
481 
482  // Then we want to fit the x hits to two tracks and the y hits to
483  // two tracks. Likely this means creating a Lutz object, feeding in
484  // the hits from the x or y view. Will need to settle on some
485  // locations for scattering planes. Probably the hit associations
486  // can be done with some sort of annealing...
487  //
488  this->FitView(*slice.get(), geo::kX, zvtx);
489  this->FitView(*slice.get(), geo::kY, zvtx);
490 
491  // Now we'll need to match the four tracks across the two views.
492  // this->ViewMatch(...);
493 
494  //
495  // Annouce that we ran to the log file...
496  //
497  std::cout << __FILE__ << ":" << __LINE__
498  << " Slice #: " << islice << std::endl;
499  } // Loop on slices
500 
501  //
502  // To keep ART quiet for now. Eventually should produce real track
503  // fits
504  //
505  std::unique_ptr<std::vector<rb::Track> >
506  tracks(new std::vector<rb::Track>);
507  evt.put(std::move(vtx));
508  evt.put(std::move(tracks));
509  } // produce() method
510 };
511 
512 //......................................................................
513 
515 
516 ////////////////////////////////////////////////////////////////////////
int fRun
Slice identifiers (run, subrun, event, slice)
DimuonFitter(fhicl::ParameterSet const &p)
A 3D position and time representing an interaction vertex.
Definition: Vertex.h:15
SubRunNumber_t subRun() const
Definition: Event.h:72
int fVerbosity
0-100, how much text output to generate
virtual double W(const rb::CellHit *chit) const
Estimate the unmeasured coordinate of chit.
Definition: Cluster.cxx:121
void FindSlices(const art::Event &evt, art::Handle< std::vector< rb::Cluster > > &h, art::PtrVector< rb::Cluster > &s)
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
TH1F * a2
Definition: f2_nu.C:545
void SetMeasurement(unsigned int i, double z, double x, double sigx)
Add measurements.
Definition: Lutz.cxx:38
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
unsigned short Plane() const
Definition: CellHit.h:39
Construct scattering surfaces for Lutz.
void reconfigure(fhicl::ParameterSet const &p)
Float_t x1[n_points_granero]
Definition: compare.C:5
geo::View_t View() const
Definition: CellHit.h:41
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
Vertical planes which measure X.
Definition: PlaneGeo.h:28
double HalfW() const
Definition: CellGeo.cxx:191
A collection of associated CellHits.
Definition: Cluster.h:47
"Break-point" track fitter
TString ip
Definition: loadincs.C:5
void produce(art::Event &evt)
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
double chi2()
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
unsigned int N() const
The number of scattering planes constructed.
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
void FitView(const rb::Cluster &clust, geo::View_t view, double zvtx)
unsigned short Cell() const
Definition: CellHit.h:40
const XML_Char * s
Definition: expat.h:262
void MakeSurfaces(const geo::GeometryBase &geo, const double *p0, const double *p1, const std::vector< double > &s)
Build scattering surfaces.
Definition: Cand.cxx:23
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
double dx[NP][NC]
Break-point track fitter.
Definition: Lutz.h:20
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
TH1F * a1
Definition: f2_nu.C:476
T get(std::string const &key) const
Definition: ParameterSet.h:231
TNtuple * fPass1Nt
pass one ntuple
bool AnaSlice(const rb::Cluster *s)
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
Geometry information for a single readout plane.
Definition: PlaneGeo.h:36
int evt
Contain 3D hit information by extrapolating ortho view coordinate from surrounding points in the orth...
const double j
Definition: BetheBloch.cxx:29
Vertex location in position and time.
Perform a "2 point" Hough transform on a collection of hits.
A very simple service to remember what detector we&#39;re working in.
size_type size() const
Definition: PtrVector.h:308
Construct a set of scattering surfaces for Lutz.
OStream cout
Definition: OStream.cxx:6
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
int fSubrun
are useful to have as class scope for logging
EventNumber_t event() const
Definition: EventID.h:116
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
T * make(ARGS...args) const
T const * get() const
Definition: Ptr.h:321
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
int fEvent
debugging information in each of the methods
Construct a track-based (x,y,z) coordinate system.
void geom(int which=0)
Definition: geom.C:163
int fSlice
used during the reconstruction.
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
const art::ProductToken< std::vector< rb::Cluster > > fSliceToken
bool getByToken(ProductToken< PROD > const &token, Handle< PROD > &result) const
Definition: DataViewImpl.h:387
Float_t e
Definition: plot.C:35
RunNumber_t run() const
Definition: Event.h:77
Encapsulate the cell geometry.
Definition: CellGeo.h:25
ProductToken< T > consumes(InputTag const &)
EventID id() const
Definition: Event.h:56
Encapsulate the geometry of one entire detector (near, far, ndos)
double SigAlpha(unsigned int i) const
The predicted RMS scattering angle at plane i.
double FindVertexZ(const rb::Cluster *c)
double S(unsigned int i) const
The location of the ith scattering plane.