Public Types | Public Member Functions | Static Public Member Functions | Public Attributes | Protected Member Functions | List of all members
bpfit::DimuonFitter Class Reference
Inheritance diagram for bpfit::DimuonFitter:
art::EDProducer art::ProducerBase art::Consumer art::EngineCreator art::ProductRegistryHelper

Public Types

using ModuleType = EDProducer
 
using WorkerType = WorkerT< EDProducer >
 
template<typename UserConfig , typename KeysToIgnore = void>
using Table = ProducerBase::Table< UserConfig, KeysToIgnore >
 

Public Member Functions

void reconfigure (fhicl::ParameterSet const &p)
 
 DimuonFitter (fhicl::ParameterSet const &p)
 
virtual ~DimuonFitter ()
 
void beginJob ()
 
void FindSlices (const art::Event &evt, art::Handle< std::vector< rb::Cluster > > &h, art::PtrVector< rb::Cluster > &s)
 
bool AnaSlice (const rb::Cluster *s)
 
double FindVertexZ (const rb::Cluster *c)
 
void FitView (const rb::Cluster &clust, geo::View_t view, double zvtx)
 
void produce (art::Event &evt)
 
template<typename PROD , BranchType B = InEvent>
ProductID getProductID (std::string const &instanceName={}) const
 
template<typename PROD , BranchType B>
ProductID getProductID (ModuleDescription const &moduleDescription, std::string const &instanceName) const
 
bool modifiesEvent () const
 
template<typename T , BranchType = InEvent>
ProductToken< Tconsumes (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ProductToken< Tconsumes (InputTag const &it)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ViewToken< TconsumesView (InputTag const &it)
 
template<typename T , BranchType = InEvent>
ProductToken< TmayConsume (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ProductToken< TmayConsume (InputTag const &it)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ViewToken< TmayConsumeView (InputTag const &it)
 
base_engine_tcreateEngine (seed_t seed)
 
base_engine_tcreateEngine (seed_t seed, std::string const &kind_of_engine_to_make)
 
base_engine_tcreateEngine (seed_t seed, std::string const &kind_of_engine_to_make, label_t const &engine_label)
 
seed_t get_seed_value (fhicl::ParameterSet const &pset, char const key[]="seed", seed_t const implicit_seed=-1)
 

Static Public Member Functions

static cet::exempt_ptr< Consumernon_module_context ()
 

Public Attributes

const art::ProductToken< std::vector< rb::Cluster > > fSliceToken
 
int fRun
 Slice identifiers (run, subrun, event, slice) More...
 
int fSubrun
 are useful to have as class scope for logging More...
 
int fEvent
 debugging information in each of the methods More...
 
int fSlice
 used during the reconstruction. More...
 
int fVerbosity
 0-100, how much text output to generate More...
 
TNtuple * fPass1Nt
 pass one ntuple More...
 

Protected Member Functions

CurrentProcessingContext const * currentContext () const
 
void validateConsumedProduct (BranchType const bt, ProductInfo const &pi)
 
void prepareForJob (fhicl::ParameterSet const &pset)
 
void showMissingConsumes () const
 

Detailed Description

Definition at line 36 of file DimuonFitter_module.cc.

Member Typedef Documentation

using art::EDProducer::ModuleType = EDProducer
inherited

Definition at line 34 of file EDProducer.h.

template<typename UserConfig , typename KeysToIgnore = void>
using art::EDProducer::Table = ProducerBase::Table<UserConfig, KeysToIgnore>
inherited

Definition at line 43 of file EDProducer.h.

using art::EDProducer::WorkerType = WorkerT<EDProducer>
inherited

Definition at line 35 of file EDProducer.h.

Constructor & Destructor Documentation

bpfit::DimuonFitter::DimuonFitter ( fhicl::ParameterSet const &  p)
inlineexplicit

Definition at line 59 of file DimuonFitter_module.cc.

References reconfigure().

59  :
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  }
int fRun
Slice identifiers (run, subrun, event, slice)
void reconfigure(fhicl::ParameterSet const &p)
const char * p
Definition: xmltok.h:285
int fSubrun
are useful to have as class scope for logging
int fEvent
debugging information in each of the methods
int fSlice
used during the reconstruction.
const art::ProductToken< std::vector< rb::Cluster > > fSliceToken
ProductToken< T > consumes(InputTag const &)
enum BeamMode string
virtual bpfit::DimuonFitter::~DimuonFitter ( )
inlinevirtual

Definition at line 73 of file DimuonFitter_module.cc.

73 { }

Member Function Documentation

bool bpfit::DimuonFitter::AnaSlice ( const rb::Cluster s)
inline

Definition at line 114 of file DimuonFitter_module.cc.

References plot_validation_datamc::c, rb::CellHit::Cell(), rb::Cluster::Cell(), fEvent, fRun, fSlice, MECModelEnuComparisons::i, geo::kX, geo::kY, rb::Cluster::NCell(), rb::Cluster::NXCell(), rb::Cluster::NYCell(), rb::CellHit::Plane(), rb::CellHit::View(), and submit_syst::x.

Referenced by produce().

114  {
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  }
int fRun
Slice identifiers (run, subrun, event, slice)
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
unsigned short Plane() const
Definition: CellHit.h:39
geo::View_t View() const
Definition: CellHit.h:41
const char * p
Definition: xmltok.h:285
Vertical planes which measure X.
Definition: PlaneGeo.h:28
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
unsigned short Cell() const
Definition: CellHit.h:40
TNtuple * fPass1Nt
pass one ntuple
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
int fEvent
debugging information in each of the methods
int fSlice
used during the reconstruction.
void bpfit::DimuonFitter::beginJob ( )
inlinevirtual

Reimplemented from art::EDProducer.

Definition at line 75 of file DimuonFitter_module.cc.

References MakeMiniprodValidationCuts::f, and art::TFileDirectory::make().

75  {
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  }
TNtuple * fPass1Nt
pass one ntuple
T * make(ARGS...args) const
template<typename T , BranchType = InEvent>
ProductToken<T> art::Consumer::consumes ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ProductToken<T> art::Consumer::consumes ( InputTag const &  it)
inherited

Definition at line 146 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

147 {
148  if (!moduleContext_)
149  return ProductToken<T>::invalid();
150 
151  consumables_[BT].emplace_back(ConsumableType::Product,
152  TypeID{typeid(T)},
153  it.label(),
154  it.instance(),
155  it.process());
156  return ProductToken<T>{it};
157 }
set< int >::iterator it
static ProductToken< T > invalid()
Definition: ProductToken.h:47
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename T , art::BranchType BT>
void art::Consumer::consumesMany ( )
inherited

Definition at line 161 of file Consumer.h.

References T.

162 {
163  if (!moduleContext_)
164  return;
165 
166  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
167 }
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::Consumer::consumesView ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ViewToken<T> art::Consumer::consumesView ( InputTag const &  it)
inherited

Definition at line 171 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

172 {
173  if (!moduleContext_)
174  return ViewToken<T>::invalid();
175 
176  consumables_[BT].emplace_back(ConsumableType::ViewElement,
177  TypeID{typeid(T)},
178  it.label(),
179  it.instance(),
180  it.process());
181  return ViewToken<T>{it};
182 }
set< int >::iterator it
static ViewToken< Element > invalid()
Definition: ProductToken.h:75
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed)
inherited
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed,
std::string const &  kind_of_engine_to_make 
)
inherited
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed,
std::string const &  kind_of_engine_to_make,
label_t const &  engine_label 
)
inherited
CurrentProcessingContext const* art::EDProducer::currentContext ( ) const
protectedinherited
void bpfit::DimuonFitter::FindSlices ( const art::Event evt,
art::Handle< std::vector< rb::Cluster > > &  h,
art::PtrVector< rb::Cluster > &  s 
)
inline

Definition at line 88 of file DimuonFitter_module.cc.

References om::cout, allTimeWatchdog::endl, art::DataViewImpl::getByToken(), make_syst_table_plots::h, MECModelEnuComparisons::i, art::PtrVector< T >::push_back(), and art::PtrVector< T >::size().

Referenced by produce().

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  }
int fVerbosity
0-100, how much text output to generate
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
size_type size() const
Definition: PtrVector.h:308
OStream cout
Definition: OStream.cxx:6
const art::ProductToken< std::vector< rb::Cluster > > fSliceToken
bool getByToken(ProductToken< PROD > const &token, Handle< PROD > &result) const
Definition: DataViewImpl.h:387
double bpfit::DimuonFitter::FindVertexZ ( const rb::Cluster c)
inline

Definition at line 175 of file DimuonFitter_module.cc.

References rb::CellHit::Cell(), geo::PlaneGeo::Cell(), rb::Cluster::Cell(), om::cout, allTimeWatchdog::endl, geom(), geo::CellGeo::GetCenter(), MECModelEnuComparisons::i, ip, calib::j, rb::Cluster::NCell(), rb::CellHit::Plane(), and geo::GeometryBase::Plane().

Referenced by produce().

175  {
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  }
int fVerbosity
0-100, how much text output to generate
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
unsigned short Plane() const
Definition: CellHit.h:39
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
TString ip
Definition: loadincs.C:5
const PlaneGeo * Plane(unsigned int i) const
unsigned short Cell() const
Definition: CellHit.h:40
const double j
Definition: BetheBloch.cxx:29
OStream cout
Definition: OStream.cxx:6
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
void geom(int which=0)
Definition: geom.C:163
void bpfit::DimuonFitter::FitView ( const rb::Cluster clust,
geo::View_t  view,
double  zvtx 
)
inline

Definition at line 251 of file DimuonFitter_module.cc.

References a1, a2, plot_validation_datamc::c, rb::CellHit::Cell(), geo::PlaneGeo::Cell(), rb::Cluster::Cell(), chi2(), om::cout, dx, e, allTimeWatchdog::endl, stan::math::exp(), check_time_usage::float, geom(), art::Ptr< T >::get(), geo::CellGeo::GetCenter(), make_syst_table_plots::h, geo::CellGeo::HalfW(), MECModelEnuComparisons::i, ip, geo::kX, geo::kY, bpfit::ScatteringSurfaces::MakeSurfaces(), bpfit::ScatteringSurfaces::N(), rb::Cluster::NCell(), plot_validation_datamc::p1, plot_validation_datamc::p2, rb::CellHit::Plane(), geo::GeometryBase::Plane(), bpfit::ScatteringSurfaces::S(), bpfit::Lutz::SetMeasurement(), bpfit::ScatteringSurfaces::SigAlpha(), std::sqrt(), registry_explorer::v, geo::PlaneGeo::View(), rb::Cluster::W(), x1, and submit_syst::x2.

Referenced by produce().

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  }
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
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.
Float_t x1[n_points_granero]
Definition: compare.C:5
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
TString ip
Definition: loadincs.C:5
const PlaneGeo * Plane(unsigned int i) const
double chi2()
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
unsigned short Cell() const
Definition: CellHit.h:40
const XML_Char * s
Definition: expat.h:262
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
TH1F * a1
Definition: f2_nu.C:476
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
Geometry information for a single readout plane.
Definition: PlaneGeo.h:36
A very simple service to remember what detector we&#39;re working in.
OStream cout
Definition: OStream.cxx:6
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
T const * get() const
Definition: Ptr.h:321
void geom(int which=0)
Definition: geom.C:163
Float_t e
Definition: plot.C:35
Encapsulate the cell geometry.
Definition: CellGeo.h:25
seed_t art::EngineCreator::get_seed_value ( fhicl::ParameterSet const &  pset,
char const  key[] = "seed",
seed_t const  implicit_seed = -1 
)
inherited
template<typename PROD , BranchType B>
ProductID art::EDProducer::getProductID ( std::string const &  instanceName = {}) const
inlineinherited

Definition at line 123 of file EDProducer.h.

References art::EDProducer::moduleDescription_.

Referenced by skim::NueSkimmer::CopyMichelSlice(), and skim::NueSkimmer::CopyMichelTrack().

124  {
125  return ProducerBase::getProductID<PROD, B>(moduleDescription_,
126  instanceName);
127  }
ModuleDescription moduleDescription_
Definition: EDProducer.h:115
template<typename PROD , BranchType B>
ProductID art::ProducerBase::getProductID ( ModuleDescription const &  moduleDescription,
std::string const &  instanceName 
) const
inherited

Definition at line 56 of file ProducerBase.h.

References art::ModuleDescription::moduleLabel().

Referenced by art::ProducerBase::modifiesEvent().

58  {
59  auto const& pd =
60  get_ProductDescription<PROD>(B, md.moduleLabel(), instanceName);
61  return pd.productID();
62  }
template<typename T , BranchType = InEvent>
ProductToken<T> art::Consumer::mayConsume ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ProductToken<T> art::Consumer::mayConsume ( InputTag const &  it)
inherited

Definition at line 189 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

190 {
191  if (!moduleContext_)
192  return ProductToken<T>::invalid();
193 
194  consumables_[BT].emplace_back(ConsumableType::Product,
195  TypeID{typeid(T)},
196  it.label(),
197  it.instance(),
198  it.process());
199  return ProductToken<T>{it};
200 }
set< int >::iterator it
static ProductToken< T > invalid()
Definition: ProductToken.h:47
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename T , art::BranchType BT>
void art::Consumer::mayConsumeMany ( )
inherited

Definition at line 204 of file Consumer.h.

References T.

205 {
206  if (!moduleContext_)
207  return;
208 
209  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
210 }
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::Consumer::mayConsumeView ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ViewToken<T> art::Consumer::mayConsumeView ( InputTag const &  it)
inherited

Definition at line 214 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

215 {
216  if (!moduleContext_)
217  return ViewToken<T>::invalid();
218 
219  consumables_[BT].emplace_back(ConsumableType::ViewElement,
220  TypeID{typeid(T)},
221  it.label(),
222  it.instance(),
223  it.process());
224  return ViewToken<T>{it};
225 }
set< int >::iterator it
static ViewToken< Element > invalid()
Definition: ProductToken.h:75
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
bool art::ProducerBase::modifiesEvent ( ) const
inlineinherited

Definition at line 40 of file ProducerBase.h.

References art::ProducerBase::getProductID(), and string.

41  {
42  return true;
43  }
static cet::exempt_ptr<Consumer> art::Consumer::non_module_context ( )
staticinherited
void art::Consumer::prepareForJob ( fhicl::ParameterSet const &  pset)
protectedinherited
void bpfit::DimuonFitter::produce ( art::Event evt)
inlinevirtual

Implements art::EDProducer.

Definition at line 442 of file DimuonFitter_module.cc.

References AnaSlice(), om::cout, DEFINE_ART_MODULE(), allTimeWatchdog::endl, art::EventID::event(), FindSlices(), FindVertexZ(), FitView(), art::Ptr< T >::get(), art::Event::id(), rb::Cluster::IsNoise(), geo::kX, geo::kY, art::Event::put(), art::Event::run(), art::PtrVector< T >::size(), and art::Event::subRun().

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
int fRun
Slice identifiers (run, subrun, event, slice)
A 3D position and time representing an interaction vertex.
Definition: Vertex.h:15
SubRunNumber_t subRun() const
Definition: Event.h:72
void FindSlices(const art::Event &evt, art::Handle< std::vector< rb::Cluster > > &h, art::PtrVector< rb::Cluster > &s)
Vertical planes which measure X.
Definition: PlaneGeo.h:28
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
void FitView(const rb::Cluster &clust, geo::View_t view, double zvtx)
bool AnaSlice(const rb::Cluster *s)
size_type size() const
Definition: PtrVector.h:308
OStream cout
Definition: OStream.cxx:6
int fSubrun
are useful to have as class scope for logging
EventNumber_t event() const
Definition: EventID.h:116
T const * get() const
Definition: Ptr.h:321
int fEvent
debugging information in each of the methods
int fSlice
used during the reconstruction.
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
RunNumber_t run() const
Definition: Event.h:77
EventID id() const
Definition: Event.h:56
double FindVertexZ(const rb::Cluster *c)
void bpfit::DimuonFitter::reconfigure ( fhicl::ParameterSet const &  p)
inline

Definition at line 52 of file DimuonFitter_module.cc.

References fhicl::ParameterSet::get().

Referenced by DimuonFitter().

53  {
54  fVerbosity = p.get<int>("Verbosity");
55  }
int fVerbosity
0-100, how much text output to generate
const char * p
Definition: xmltok.h:285
void art::Consumer::showMissingConsumes ( ) const
protectedinherited

Referenced by art::RootOutput::endJob().

void art::Consumer::validateConsumedProduct ( BranchType const  bt,
ProductInfo const &  pi 
)
protectedinherited

Member Data Documentation

int bpfit::DimuonFitter::fEvent

debugging information in each of the methods

Definition at line 45 of file DimuonFitter_module.cc.

Referenced by AnaSlice().

TNtuple* bpfit::DimuonFitter::fPass1Nt

pass one ntuple

Definition at line 48 of file DimuonFitter_module.cc.

int bpfit::DimuonFitter::fRun

Slice identifiers (run, subrun, event, slice)

Definition at line 43 of file DimuonFitter_module.cc.

Referenced by AnaSlice().

int bpfit::DimuonFitter::fSlice

used during the reconstruction.

Definition at line 46 of file DimuonFitter_module.cc.

Referenced by AnaSlice().

const art::ProductToken<std::vector<rb::Cluster> > bpfit::DimuonFitter::fSliceToken

Definition at line 39 of file DimuonFitter_module.cc.

int bpfit::DimuonFitter::fSubrun

are useful to have as class scope for logging

Definition at line 44 of file DimuonFitter_module.cc.

int bpfit::DimuonFitter::fVerbosity

0-100, how much text output to generate

Definition at line 47 of file DimuonFitter_module.cc.


The documentation for this class was generated from the following file: