25 #include "Utilities/AssociationUtil.h" 34 namespace bpfit {
class DimuonFitter; }
54 fVerbosity = p.
get<
int>(
"Verbosity");
66 this->produces< std::vector<rb::Vertex> >();
67 this->produces< std::vector<rb::Track> >();
77 fPass1Nt = f->
make<TNtuple>(
"p1",
94 for (i=0; i<
h->size(); ++
i) {
99 <<
" Found " << s.
size() <<
" slices." 130 std::vector<int> hitp(1000);
131 for (i=0; i<s->
NCell(); ++
i) {
178 <<
" Start 'FindVertexZ'. " 194 std::vector<unsigned int>
p(c->
NCell());
195 unsigned int pvtx = 999999;
196 for (i=0; i<c->
NCell(); ++
i) {
198 if (ip<pvtx) pvtx =
ip;
203 <<
" Set default vertex plane to pvtx=" 206 std::sort(
p.begin(),
p.end());
208 const unsigned int N = 2;
209 const unsigned int M = N+2;
210 for (i=0, j=N; j+1<c->
NCell(); ++
i,++
j) {
211 if ( (
p[j]-
p[i])<=M ) {
215 <<
" Plane " << pvtx <<
" satisfies " 218 <<
" pi = " <<
p[
i] <<
" condition." 228 for (i=0; i<c->
NCell(); ++
i) {
236 <<
" Found vertex z to z=" 255 <<
" Starting FitView for view " << view
259 const double zero_weight = 1
e-9;
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) {
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;
280 <<
" Found slice extent " 281 <<
" [" << plo <<
"," << phi <<
"] " 282 <<
" [" << clo <<
"," << chi <<
"]" 290 std::vector<double>
x1;
291 std::vector<double> z1;
292 std::vector<double> s1;
293 std::vector<double> w1;
295 std::vector<double>
x2;
296 std::vector<double> z2;
297 std::vector<double> s2;
298 std::vector<double> w2;
301 double zhi = -999999;
302 double rt12 =
sqrt(12.0);
303 for (i=0; i<clust.
NCell(); ++
i) {
306 unsigned int ic = h->
Cell();
316 if (v != view)
continue;
320 if (xyz[2]<zlo) zlo = xyz[2];
321 if (xyz[2]>zhi) zhi = xyz[2];
322 double sigxy = c->
HalfW()/rt12;
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]); }
328 z1.push_back(xyz[2]);
329 z2.push_back(xyz[2]);
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);
352 for (i=0; i<z1.size(); ++
i) {
353 std::cout << __FILE__ <<
":" << __LINE__ <<
"\t" 355 << z1[
i] <<
"\t" << z2[
i] <<
"\t" 356 << x1[
i] <<
"\t" << x2[
i] <<
"\t" 357 << s1[
i] <<
"\t" << s2[
i] <<
"\t" 363 <<
" Initialized fits." 377 double p1[3] = {3,3,zlo};
378 double p2[3] = {3,3,zhi};
380 std::vector<double>
s(z1);
383 sort(s.begin(), s.end());
386 std::cout << __FILE__ <<
":" << __LINE__ <<
" Made " 387 << scatter.
N() <<
" scattering surfaces between z=" 393 double lambda = 1/25.;
394 for (
unsigned int itr=0; itr<5; ++itr) {
397 <<
" Fitting iteration " << itr
402 for (i=0; i<z1.size(); ++
i) {
405 for (i=0; i<z2.size(); ++
i) {
406 lutz2.SetMeasurement(i, z2[i], x2[i], s2[i]/w2[i]);
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));
415 lutz1.Fit(&a1, &b1, 0, &chi1);
416 lutz2.Fit(&a2, &b2, 0, &chi2);
419 <<
" Finished fits in view " << view <<
"." 420 <<
" chi1 = " << chi1
421 <<
" chi2 = " << chi2
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;
447 std::unique_ptr< std::vector<rb::Vertex> >
448 vtx(
new std::vector<rb::Vertex>);
462 for (islice=0; islice<slicep.
size(); ++islice) {
470 if (slice->
IsNoise())
continue;
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));
int fRun
Slice identifiers (run, subrun, event, slice)
DimuonFitter(fhicl::ParameterSet const &p)
A 3D position and time representing an interaction vertex.
SubRunNumber_t subRun() const
int fVerbosity
0-100, how much text output to generate
virtual double W(const rb::CellHit *chit) const
Estimate the unmeasured coordinate of chit.
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.
void GetCenter(double *xyz, double localz=0.0) const
void SetMeasurement(unsigned int i, double z, double x, double sigx)
Add measurements.
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
unsigned short Plane() const
Construct scattering surfaces for Lutz.
void reconfigure(fhicl::ParameterSet const &p)
Float_t x1[n_points_granero]
const CellGeo * Cell(int icell) const
Vertical planes which measure X.
A collection of associated CellHits.
"Break-point" track fitter
void produce(art::Event &evt)
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
Horizontal planes which measure Y.
unsigned int N() const
The number of scattering planes constructed.
ProductID put(std::unique_ptr< PROD > &&product)
void FitView(const rb::Cluster &clust, geo::View_t view, double zvtx)
unsigned short Cell() const
void MakeSurfaces(const geo::GeometryBase &geo, const double *p0, const double *p1, const std::vector< double > &s)
Build scattering surfaces.
View_t View() const
Which coordinate does this plane measure.
Break-point track fitter.
void push_back(Ptr< U > const &p)
T get(std::string const &key) const
TNtuple * fPass1Nt
pass one ntuple
bool AnaSlice(const rb::Cluster *s)
fvar< T > exp(const fvar< T > &x)
Geometry information for a single readout plane.
Contain 3D hit information by extrapolating ortho view coordinate from surrounding points in the orth...
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're working in.
Construct a set of scattering surfaces for Lutz.
unsigned int NYCell() const
Number of cells in the y-view.
int fSubrun
are useful to have as class scope for logging
EventNumber_t event() const
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
A rawdata::RawDigit with channel information decoded.
T * make(ARGS...args) const
unsigned int NXCell() const
Number of cells in the x-view.
int fEvent
debugging information in each of the methods
Construct a track-based (x,y,z) coordinate system.
int fSlice
used during the reconstruction.
bool IsNoise() const
Is the noise flag set?
const art::ProductToken< std::vector< rb::Cluster > > fSliceToken
bool getByToken(ProductToken< PROD > const &token, Handle< PROD > &result) const
Encapsulate the cell geometry.
ProductToken< T > consumes(InputTag const &)
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.