SplitTracks_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: SplitTracks
3 // Module Type: producer
4 // File: SplitTracks_module.cc
5 //
6 // Generated at Tue Apr 23 11:01:33 2013 by Denver Whittington using artmod
7 // from art v1_02_06.
8 ////////////////////////////////////////////////////////////////////////
9 
10 // ART includes
17 
18 // C++ includes
19 #include <vector>
20 #include <string>
21 #include <cstdlib>
22 #include <iostream>
23 
24 // NOvASoft includes
26 #include "Geometry/Geometry.h"
27 #include "GeometryObjects/Geo.h"
29 #include "RecoBase/Track.h"
32 #include "Utilities/AssociationUtil.h"
33 
34 // ROOT includes
35 #include "TTree.h"
36 #include "TH1.h"
37 #include "TH2.h"
38 
39 namespace align {
40  class SplitTracks;
41 }
42 
44 public:
45  explicit SplitTracks(fhicl::ParameterSet const &pset);
46  virtual ~SplitTracks();
47 
48  void analyze(art::Event const &evt) override;
49  //void produce(art::Event &evt);
50  void beginJob() override;
51  //void endJob() override;
52  void ClearVectors();
53 
54  double DemingRegFit( const std::vector<double>& x,
55  const std::vector<double>& y,
56  const std::vector<double>& w,
57  double& sigXY, double& sigZ,
58  double& x1, double& y1,
59  double& x2, double& y2,
60  double& sigy1, double& sigy2 );
61 
62  TVector3 GeoEstimatePos( art::Ptr<rb::CellHit> cellHit, double W );
63  void GetCellEndpoints( art::Ptr<rb::CellHit> cell, double *cellLo, double *cellHi );
64 
65 private:
66 
69  unsigned int fMinPlanesInSubCluster;
70 
73  TTree *fGeoTree;
74  TTree *fOutTree;
75 
76  // Geometry details
77  std::vector<double> blockCenterX;
78  std::vector<double> blockCenterY;
79  std::vector<double> blockCenterZ;
80 
81  // Track details
82  int trk_ID;
85  std::vector<int> subtrk_Block;
86  std::vector<double> subtrk_StartX, subtrk_StartY, subtrk_StartZ;
88  std::vector<double> subtrk_StopX, subtrk_StopY, subtrk_StopZ;
90 
91  // Monitoring histograms
92  TH1D *xzViewTrkRes;
93  TH1D *yzViewTrkRes;
110 
114 
118 
119 };
120 
121 // Constructor
123  : EDAnalyzer(pset)
124  , fTrackInLabel(pset.get<std::string>("TrackInLabel"))
125  , fMinBlocksInTrack(pset.get<int>("MinBlocksInTrack"))
126  , fMinPlanesInSubCluster(pset.get<unsigned int>("MinPlanesInSubCluster"))
127  , fTrkAlg(0)
128  , fTrackAlgEnum(pset.get< int >("TrackingAlgEnumeration", 0))
129 {
130 
131  switch(fTrackAlgEnum){
132  case 0:
133  fTrkAlg = new trk::CosmicTrackAlg(pset.get<fhicl::ParameterSet>("TrackAlgPSet"));
134  break;
135  case 1:
136  fTrkAlg = new trk::WindowTrackingAlg(pset.get<fhicl::ParameterSet>("TrackAlgPSet"));
137  break;
138  default:
139  fTrkAlg = new trk::CosmicTrackAlg(pset.get<fhicl::ParameterSet>("TrackAlgPSet"));
140  break;
141  }
142 }
143 
144 // Destructor
146 {
147  if (fTrkAlg) delete fTrkAlg;
148 }
149 
150 // Begin job
153 
154  // Output TTree
155  fGeoTree = tfs->make<TTree>("GeometryInfo","Geometry Info");
156  fOutTree = tfs->make<TTree>("SplitTracks","Split Tracks");
157 
158  // Geometry details
159  fGeoTree->Branch("blockCenterX", &blockCenterX);
160  fGeoTree->Branch("blockCenterY", &blockCenterY);
161  fGeoTree->Branch("blockCenterZ", &blockCenterZ);
162  // Get Geo details...
163  fGeoTree->Fill();
164 
165  // Track details
166  fOutTree->Branch("trk_ID", &trk_ID);
167  fOutTree->Branch("trk_StartX", &trk_StartX);
168  fOutTree->Branch("trk_StartY", &trk_StartY);
169  fOutTree->Branch("trk_StartZ", &trk_StartZ);
170  fOutTree->Branch("trk_StopX", &trk_StopX);
171  fOutTree->Branch("trk_StopY", &trk_StopY);
172  fOutTree->Branch("trk_StopZ", &trk_StopZ);
173 
174  fOutTree->Branch("subtrk_Block", &subtrk_Block);
175  fOutTree->Branch("subtrk_StartX", &subtrk_StartX);
176  fOutTree->Branch("subtrk_StartY", &subtrk_StartY);
177  fOutTree->Branch("subtrk_StartZ", &subtrk_StartZ);
178  fOutTree->Branch("subtrk_ErrStartX", &subtrk_ErrStartX);
179  fOutTree->Branch("subtrk_ErrStartY", &subtrk_ErrStartY);
180  fOutTree->Branch("subtrk_ErrStartZ", &subtrk_ErrStartZ);
181  fOutTree->Branch("subtrk_StopX", &subtrk_StopX);
182  fOutTree->Branch("subtrk_StopY", &subtrk_StopY);
183  fOutTree->Branch("subtrk_StopZ", &subtrk_StopZ);
184  fOutTree->Branch("subtrk_ErrStopX", &subtrk_ErrStopX);
185  fOutTree->Branch("subtrk_ErrStopY", &subtrk_ErrStopY);
186  fOutTree->Branch("subtrk_ErrStopZ", &subtrk_ErrStopZ);
187 
188  subtrkXYstart = tfs->make<TH2D>("subtrkXYstart","Downstream Subtrack Start (XY);X-view Cell;Y-view Cell",100,-775,775,100,-775,775);
189  subtrkXZstart = tfs->make<TH2D>("subtrkXZstart","Downstream Subtrack Start (XZ);Plane;X-view Cell",600,0,6000,1550,-775,775);
190  subtrkYZstart = tfs->make<TH2D>("subtrkYZstart","Downstream Subtrack Start (YZ);Plane;Y-view Cell",600,0,6000,1550,-775,775);
191 
192  xzViewTrkRes = tfs->make<TH1D>("xzViewTrkRes","Track Residual (X-View);Residual",100,-10,10);
193  yzViewTrkRes = tfs->make<TH1D>("yzViewTrkRes","Track Residual (Y-View);Residual",100,-10,10);
194  xzViewTrkResVsX = tfs->make<TH2D>("xzViewTrkResVsX","Track Residual (X-View) vs Cell Position (X);Cell Number;Residual",32*12,0,32*12,100,-10,10);
195  xzViewTrkResVsZ = tfs->make<TH2D>("xzViewTrkResVsZ","Track Residual (X-View) vs Cell Position (Z);Plane Number;Residual",32*28,0,32*28,100,-10,10);
196  yzViewTrkResVsY = tfs->make<TH2D>("yzViewTrkResVsY","Track Residual (Y-View) vs Cell Position (Y);Cell Number;Residual",32*12,0,32*12,100,-10,10);
197  yzViewTrkResVsZ = tfs->make<TH2D>("yzViewTrkResVsZ","Track Residual (Y-View) vs Cell Position (Z);Plane Number;Residual",32*28,0,32*28,100,-10,10);
198 
199  xzViewSubtrkRes = tfs->make<TH1D>("xzViewSubtrkRes","SubTrack Residual (X-View);Residual",100,-10,10);
200  yzViewSubtrkRes = tfs->make<TH1D>("yzViewSubtrkRes","SubTrack Residual (Y-View);Residual",100,-10,10);
201  xzViewSubtrkResVsX = tfs->make<TH2D>("xzViewSubtrkResVsX","SubTrack Residual (X-View) vs Cell Position (X);Cell Number;Residual",32*12,0,32*12,100,-10,10);
202  xzViewSubtrkResVsZ = tfs->make<TH2D>("xzViewSubtrkResVsZ","SubTrack Residual (X-View) vs Cell Position (Z);Plane Number;Residual",32*28,0,32*28,100,-10,10);
203  yzViewSubtrkResVsY = tfs->make<TH2D>("yzViewSubtrkResVsY","SubTrack Residual (Y-View) vs Cell Position (Y);Cell Number;Residual",32*12,0,32*12,100,-10,10);
204  yzViewSubtrkResVsZ = tfs->make<TH2D>("yzViewSubtrkResVsZ","SubTrack Residual (Y-View) vs Cell Position (Z);Plane Number;Residual",32*28,0,32*28,100,-10,10);
205 
206  xzViewRefitRes = tfs->make<TH1D>("xzViewRefitRes","Refitted SubTrack Residual (X-View);Residual",100,-10,10);
207  yzViewRefitRes = tfs->make<TH1D>("yzViewRefitRes","Refitted SubTrack Residual (Y-View);Residual",100,-10,10);
208  xzViewRefitResVsX = tfs->make<TH2D>("xzViewRefitResVsX","Refitted SubTrack Residual (X-View) vs Cell Position (X);Cell Number;Residual",32*12,0,32*12,100,-10,10);
209  xzViewRefitResVsZ = tfs->make<TH2D>("xzViewRefitResVsZ","Refitted SubTrack Residual (X-View) vs Cell Position (Z);Plane Number;Residual",32*28,0,32*28,100,-10,10);
210  yzViewRefitResVsY = tfs->make<TH2D>("yzViewRefitResVsY","Refitted SubTrack Residual (Y-View) vs Cell Position (Y);Cell Number;Residual",32*12,0,32*12,100,-10,10);
211  yzViewRefitResVsZ = tfs->make<TH2D>("yzViewRefitResVsZ","Refitted SubTrack Residual (Y-View) vs Cell Position (Z);Plane Number;Residual",32*28,0,32*28,100,-10,10);
212 
213  trackDCAtoCell = tfs->make<TH1D>("trackDCAtoCell","Track DCA to Cell",100,0,10);
214  subTrackDCAtoCell = tfs->make<TH1D>("subTrackDCAtoCell","SubTrack DCA to Cell",100,0,10);
215  refitDCAtoCell = tfs->make<TH1D>("refitDCAtoCell","Refitted SubTrack DCA to Cell",100,0,10);
216 
217  return;
218 }
219 
221  trk_ID = 0;
222  trk_StartX = 0; trk_StartY = 0; trk_StartZ = 0;
223  trk_StopX = 0; trk_StopY = 0; trk_StopZ = 0;
224  subtrk_Block.clear();
225  subtrk_StartX.clear();
226  subtrk_StartY.clear();
227  subtrk_StartZ.clear();
228  subtrk_ErrStartX.clear();
229  subtrk_ErrStartY.clear();
230  subtrk_ErrStartZ.clear();
231  subtrk_StopX.clear();
232  subtrk_StopY.clear();
233  subtrk_StopZ.clear();
234  subtrk_ErrStopX.clear();
235  subtrk_ErrStopY.clear();
236  subtrk_ErrStopZ.clear();
237  return;
238 }
239 
240 // Produce
242 //void align::SplitTracks::produce(art::Event &evt)
243 {
245  const int detId = geoSvc->DetId(); // What detector am I in?
247 
248  // Get the tracks that were made by CosmicTrack
249  // and split along alignment boundaries
250  /* (blocks for now) */
252  evt.getByLabel(fTrackInLabel,CosmicTracks);
253  if ( CosmicTracks->empty() ) {
254  std::cout << "No tracks found in " << fTrackInLabel << " for this event!" << std::endl;
255  return;
256  }
257 
258  // Loop through CosmicTracks
259  std::unique_ptr< std::vector<rb::Track> > subTrackColl(new std::vector<rb::Track>);
260  for ( size_t trkItr = 0; trkItr < CosmicTracks->size(); ++trkItr ) {
261  art::Ptr<rb::Track> trackIn(CosmicTracks,trkItr);
262 
263  if ( !trackIn->Is3D() ) continue;
264  this->ClearVectors();
265 
266  /* Track Requirements:
267  * Crosses a plane boundary
268  * Contains at least 3 hits in XZ and YZ in each plane
269  */
270 
271  // Crosses a plane boundary?
272  daqchannelmap::lchan logchan = cmapSvc->Map()->encodeLChan(detId, trackIn->MaxPlane(), 1);
273  daqchannelmap::lchan logchan1 = cmapSvc->Map()->encodeLChan(detId, trackIn->MinPlane(), 1);
274  int maxBlock = cmapSvc->Map()->computeBlock(logchan);
275  int minBlock = cmapSvc->Map()->computeBlock(logchan1);
276  if ( maxBlock < minBlock ) {
277  std::cout << "MaxBlock is smaller than MinBlock for this track!" << std::endl;
278  std::swap(maxBlock,minBlock);
279  // This shouldn't happen, so throw a warning if it does...
280  }
281  if ( (maxBlock - minBlock + 1) < fMinBlocksInTrack ) continue;
282 
283  // Does the track go through the top and bottom of the detector bigbox?
284  /*
285  double start[3] = { trackIn->Start().X(), trackIn->Start().Y(), trackIn->Start().Z() };
286  double dir[3] = { trackIn->Dir().X(), trackIn->Dir().Y(), trackIn->Dir().Z() };
287  double boxMin[3];
288  double boxMax[3];
289  geoSvc->DetectorBigBox(&boxMin[0],&boxMax[0],&boxMin[1],&boxMax[1],&boxMin[2],&boxMax[2]);
290  bool thruTop = geo::IntersectsBox( start, dir, boxMin[0], boxMax[0], boxMax[1], boxMax[1], boxMin[2], boxMax[2] );
291  bool thruBot = geo::IntersectsBox( start, dir, boxMin[0], boxMax[0], boxMin[1], boxMin[1], boxMin[2], boxMax[2] );
292  */
293  //std::cout << "Min: " << trackIn->MinCell( geo::kY ) << " Max: " << trackIn->MaxCell( geo::kY ) << std::endl;
294  bool thruTop = ( trackIn->MaxCell( geo::kY ) == 32*12-1 );
295  bool thruBot = ( trackIn->MinCell( geo::kY ) == 0 );
296  if (!(thruTop&&thruBot) ) continue;
297 
298  // Record the parameters for this track
299  trk_ID = trackIn->ID();
300  trk_StartX = trackIn->Start().X();
301  trk_StartY = trackIn->Start().Y();
302  trk_StartZ = trackIn->Start().Z();
303  trk_StopX = trackIn->Stop().X();
304  trk_StopY = trackIn->Stop().Y();
305  trk_StopZ = trackIn->Stop().Z();
306 
307  // Group CellHits into new clusters by block
308  std::vector<rb::Cluster> subClusters( maxBlock-minBlock+1, rb::Cluster() );
309  // Loop through hits in track and sort into subClusters by block
310  for ( size_t cellItr = 0; cellItr < trackIn->NCell(); ++cellItr ) {
311  const art::Ptr<rb::CellHit>& cell = trackIn->Cell(cellItr);
312  daqchannelmap::lchan logchan2 = cmapSvc->Map()->encodeLChan(detId, cell->Plane(),1);
313  int block = cmapSvc->Map()->computeBlock(logchan2);
314  subClusters.at(block-minBlock).Add(cell); // WEIGHTS???
315  /*
316  double W = trackIn->W(cell.get());
317  double cellPos[3];
318  geoSvc->Plane(cell->Plane())->Cell(cell->Cell())->GetCenter(cellPos,W);
319  */
320  double cellLo[3], cellHi[3];
321  this->GetCellEndpoints( cell, cellLo, cellHi );
322  double dOnTrack, dInCell;
323  TVector3 pcaOnTrack, pcaInCell;
324  double dca = geo::ClosestApproach( trackIn->Start(), trackIn->Stop(), TVector3(cellLo), TVector3(cellHi), &dOnTrack, &dInCell, &pcaOnTrack, &pcaInCell );
325  trackDCAtoCell->Fill(dca);
326  double cellPos[3] = { pcaInCell.X(), pcaInCell.Y(), pcaInCell.Z() };
327  double trkX(0), trkY(0);
328  trackIn->InterpolateXY(cellPos[2],trkX,trkY);
329  if ( cell->View() == geo::kX ) {
330  double res = trkX-cellPos[0];
331  xzViewTrkRes->Fill(res);
332  xzViewTrkResVsX->Fill(cell->Cell(),res);
333  xzViewTrkResVsZ->Fill(cell->Plane(),res);
334  }
335  if ( cell->View() == geo::kY ) {
336  double res = trkY-cellPos[1];
337  yzViewTrkRes->Fill(res);
338  yzViewTrkResVsY->Fill(cell->Cell(),res);
339  yzViewTrkResVsZ->Fill(cell->Plane(),res);
340  }
341  }
342 
343  /*
344  // Report!
345  std::cout << "Hey! This track had " << trackIn->NCell() << " cell hits." << std::endl;
346  std::cout << " They were distributed across " << subClusters.size() << " blocks, with" << std::endl;
347  for ( size_t i = 0; i < subClusters.size(); ++i ) {
348  std::cout << " " << subClusters[i].NCell() << " in block " << i << std::endl;
349  }
350  */
351 
352  // Loop through our subClusters and fit them to get end-points and errors
353  for ( size_t clusItr = 0; clusItr < subClusters.size(); ++clusItr ) {
354 
355  // Make a track so we can refine its hit positions with the extrapolated W positions
356  std::vector<rb::Track> subTracks = fTrkAlg->MakeTrack(&subClusters.at(clusItr));
357 
358  // only use subtracks if the vector is of size 1
359  if(subTracks.size() < 1 || subTracks.size() > 1) continue;
360 
361  rb::Track subTrack = subTracks.front();
362 
363  // Does this cluster contain enough Planes?
364  unsigned int nPlanes = subTrack.ExtentPlane();
365  if ( nPlanes < fMinPlanesInSubCluster ) continue;
366  if ( subTrack.NXCell() < 2 ) continue;
367  if ( subTrack.NYCell() < 2 ) continue;
368  subTrackColl->push_back(subTrack);
369 
370  // Refit the two views to get the uncertainties
371  std::vector<double> xzViewX, xzViewZ, yzViewY, yzViewZ, xzViewW, yzViewW;
372  double zMin = 10000000;
373  double zMax = 0.;
374  for ( size_t i = 0; i < subTrack.NCell(); ++i ) {
375  art::Ptr<rb::CellHit> cellHit = subTrack.Cell(i);
376  /*
377  double W = subTrack.W(cellHit.get());
378  double pos[3];
379  geoSvc->Plane(cellHit->Plane())->Cell(cellHit->Cell())->GetCenter(pos,W);
380  */
381  double cellLo[3], cellHi[3];
382  this->GetCellEndpoints( cellHit, cellLo, cellHi );
383  double dOnTrack, dInCell;
384  TVector3 pcaOnTrack, pcaInCell;
385  double dca = geo::ClosestApproach( subTrack.Start(), subTrack.Stop(), TVector3(cellLo), TVector3(cellHi), &dOnTrack, &dInCell, &pcaOnTrack, &pcaInCell );
386  subTrackDCAtoCell->Fill(dca);
387  double pos[3] = { pcaInCell.X(), pcaInCell.Y(), pcaInCell.Z() };
388  //double localXCheck[3]; geoSvc->Plane(cellHit->Plane())->Cell(cellHit->Cell())->WorldToLocal(pos,localXCheck);
389  //std::cout << localXCheck[0] << " " << localXCheck[1] << " " << localXCheck[2] << std::endl;
390  double subtrkX(0), subtrkY(0);
391  subTrack.InterpolateXY(pos[2],subtrkX,subtrkY);
392  if ( cellHit->View() == geo::kX ) {
393  double res = subtrkX-pos[0];
394  xzViewSubtrkRes->Fill(res);
395  xzViewSubtrkResVsX->Fill(cellHit->Cell(),res);
396  xzViewSubtrkResVsZ->Fill(cellHit->Plane(),res);
397  xzViewX.push_back(pos[0]);
398  xzViewZ.push_back(pos[2]);
399  xzViewW.push_back(1.);
400  }
401  if ( cellHit->View() == geo::kY ) {
402  double res = subtrkY-pos[1];
403  yzViewSubtrkRes->Fill(res);
404  yzViewSubtrkResVsY->Fill(cellHit->Cell(),res);
405  yzViewSubtrkResVsZ->Fill(cellHit->Plane(),res);
406  yzViewY.push_back(pos[1]);
407  yzViewZ.push_back(pos[2]);
408  yzViewW.push_back(1.);
409  }
410  if ( pos[2] < zMin ) zMin = pos[2];
411  if ( pos[2] > zMax ) zMax = pos[2];
412  }
413 
414  // Fits
415  double sigX = 3./sqrt(12);
416  double sigY = 3./sqrt(12);
417  double sigZ = 5./sqrt(12);
418  double xStart, yStart, zStart;
419  double xStartErr, yStartErr, zStartErr;
420  double xStop, yStop, zStop;
421  double xStopErr, yStopErr, zStopErr;
422  zStart = zMin; zStartErr = 0.;//sigZ;
423  zStop = zMax; zStopErr = 0.;//sigZ;
424  this->DemingRegFit( xzViewZ, xzViewX, xzViewW, sigZ, sigX, zStart, xStart, zStop, xStop, xStartErr, xStopErr );
425  this->DemingRegFit( yzViewZ, yzViewY, yzViewW, sigZ, sigY, zStart, yStart, zStop, yStop, yStartErr, yStopErr );
426 
427  // Residuals from refits
428  for ( size_t i = 0; i < subTrack.NXCell(); ++i ) {
429  art::Ptr<rb::CellHit> cellHit = subTrack.XCell(i);
430  /*
431  double W = subTrack.W(cellHit.get());
432  double pos[3];
433  geoSvc->Plane(cellHit->Plane())->Cell(cellHit->Cell())->GetCenter(pos,W);
434  */
435  double cellLo[3], cellHi[3];
436  this->GetCellEndpoints( cellHit, cellLo, cellHi );
437  double dOnTrack, dInCell;
438  TVector3 pcaOnTrack, pcaInCell;
439  double dca = geo::ClosestApproach( TVector3(xStart,yStart,zStart), TVector3(xStop,yStop,zStop), TVector3(cellLo), TVector3(cellHi), &dOnTrack, &dInCell, &pcaOnTrack, &pcaInCell );
440  refitDCAtoCell->Fill(dca);
441  double pos[3] = { pcaInCell.X(), pcaInCell.Y(), pcaInCell.Z() };
442  double slopeX = (xStop-xStart)/(zStop-zStart);
443  double extrapX = slopeX*(pos[2]-zStart)+xStart;
444  double res = extrapX-pos[0];
445  xzViewRefitRes->Fill(res);
446  xzViewRefitResVsX->Fill(cellHit->Cell(),res);
447  xzViewRefitResVsZ->Fill(cellHit->Plane(),res);
448  }
449  for ( size_t i = 0; i < subTrack.NYCell(); ++i ) {
450  art::Ptr<rb::CellHit> cellHit = subTrack.YCell(i);
451  /*
452  double W = subTrack.W(cellHit.get());
453  double pos[3];
454  geoSvc->Plane(cellHit->Plane())->Cell(cellHit->Cell())->GetCenter(pos,W);
455  */
456  double cellLo[3], cellHi[3];
457  this->GetCellEndpoints( cellHit, cellLo, cellHi );
458  double dOnTrack, dInCell;
459  TVector3 pcaOnTrack, pcaInCell;
460  double dca = geo::ClosestApproach( TVector3(xStart,yStart,zStart), TVector3(xStop,yStop,zStop), TVector3(cellLo), TVector3(cellHi), &dOnTrack, &dInCell, &pcaOnTrack, &pcaInCell );
461  refitDCAtoCell->Fill(dca);
462  double pos[3] = { pcaInCell.X(), pcaInCell.Y(), pcaInCell.Z() };
463  double slopeY = (yStop-yStart)/(zStop-zStart);
464  double extrapY = slopeY*(pos[2]-zStart)+yStart;
465  double res = extrapY-pos[1];
466  yzViewRefitRes->Fill(res);
467  yzViewRefitResVsY->Fill(cellHit->Cell(),res);
468  yzViewRefitResVsZ->Fill(cellHit->Plane(),res);
469  }
470  subtrkXYstart->Fill(xStart,yStart);
471  subtrkXZstart->Fill(zStart,xStart);
472  subtrkYZstart->Fill(zStart,yStart);
473 
474  //if ( xStartErr > 1000 || xStopErr > 1000 ) continue;
475 
476  // Record the parameters for this subtrack
477  daqchannelmap::lchan logchan3 = cmapSvc->Map()->encodeLChan(detId, subTrack.MinPlane(),1);
478  subtrk_Block.push_back( cmapSvc->Map()->computeBlock(logchan3));
479  subtrk_StartX.push_back( xStart );
480  subtrk_StartY.push_back( yStart );
481  subtrk_StartZ.push_back( zStart );
482  subtrk_ErrStartX.push_back( xStartErr );
483  subtrk_ErrStartY.push_back( yStartErr );
484  subtrk_ErrStartZ.push_back( zStartErr );
485  subtrk_StopX.push_back( xStop );
486  subtrk_StopY.push_back( yStop );
487  subtrk_StopZ.push_back( zStop );
488  subtrk_ErrStopX.push_back( xStopErr );
489  subtrk_ErrStopY.push_back( yStopErr );
490  subtrk_ErrStopZ.push_back( zStopErr );
491 
492  }
493 
494  //std::cout << "This track had " << TracksOut.size() << " subtracks." << std::endl;
495 
496  // Save properties of this track and each of its subtracks to ROOT TTree
497  fOutTree->Fill();
498 
499  } // End loop through CosmicTracks
500  //evt.put(std::move(subTrackColl));
501 
502  return;
503 }
504 
505 //......................................................................
506 ///
507 /// \brief Find the best-fit line to a collection of points in 2-D
508 /// using Deming regression
509 ///
510 /// This is similar to LinFitMinDPerp, but uses errors on hits
511 /// http://en.wikipedia.org/wiki/Deming_regression
512 ///
513 /// @param x - input vector of x coordinates
514 /// @param y - input vector of y coordinates
515 /// @param w - input vector of weights for the points
516 /// @param sigX - uncertainty on cell hit position along x-axis
517 /// @param sigY - uncertainty on cell hit position along y-axis
518 /// @param x1 - position of one end of segment (given)
519 /// @param x2 - position of other end of segment (given)
520 /// @param y1 - coordinate of line extrapolated to x1
521 /// @param y2 - coordinate of line extrapolated to x2
522 /// @param sigy1 - uncertainty on extrapolated coordinate y1
523 /// @param sigy2 - uncertainty on extrapolated coordinate y2
524 /// @returns The chi^2 value defined by chi^2 = sum_i[w_i d^2_i]
525 ///
526 double align::SplitTracks::DemingRegFit(const std::vector<double>& x,
527  const std::vector<double>& y,
528  const std::vector<double>& w,
529  double& sigX, double& sigY,
530  double& x1, double& y1,
531  double& x2, double& y2,
532  double& sigy1, double& sigy2 )
533 {
534  double delta = sigY / sigX;
535 
536  // Before going ahead, make sure we have sensible arrays
537  assert(x.size()==y.size());
538  assert(x.size()==w.size());
539  assert(x.size()>=2);
540  unsigned int i;
541 
542  // Accumulate the sums needed to compute the fit line
543  double xAve = 0.0;
544  double yAve = 0.0;
545  int nHit = 0;
546  for (i=0; i<w.size(); ++i) {
547  if (w[i] == 0) continue;
548  nHit++;
549  xAve += x[i];
550  yAve += y[i];
551  }
552  xAve /= double(nHit);
553  yAve /= double(nHit);
554  double Sxy = 0.0;
555  double Syy = 0.0;
556  double Sxx = 0.0;
557  for (i=0; i<w.size(); ++i) {
558  if (w[i] == 0) continue;
559  Sxx += (x[i] - xAve)*(x[i] - xAve);
560  Syy += (y[i] - yAve)*(y[i] - yAve);
561  Sxy += (x[i] - xAve)*(y[i] - yAve);
562  }
563  Sxx /= double(nHit-1);
564  Syy /= double(nHit-1);
565  Sxy /= double(nHit-1);
566 
567  // The best fit through the points
568  double M1;
569  //double B1;
570  double xydiff = Syy - delta*Sxx;
571  double radic = sqrt( (xydiff*xydiff) + 4*delta*Sxy*Sxy );
572  double numer = xydiff + radic;
573  double denom = 2.*Sxy;
574  M1 = ( numer ) / ( denom );
575  //B1 = yAve - M1*xAve;
576 
577  y1 = M1*(x1-xAve) + yAve;
578  y2 = M1*(x2-xAve) + yAve;
579 
580  // Compute errors
581  double sig2y1(0), sig2y2(0);
582  for (i=0; i<w.size(); ++i) {
583  if (w[i] == 0) continue;
584  double dnumdx = 0.;
585  double dnumdy = 0.;
586  dnumdx += delta*(2/double(nHit-1))*(x[i]-xAve);
587  dnumdx += 1./(2*radic)*( 2*xydiff*( -delta*2/double(nHit-1)*(x[i]-xAve) )
588  + 8*delta*Sxy*1./double(nHit-1)*(y[i]-yAve) );
589  dnumdy += (2/double(nHit-1))*(y[i]-yAve);
590  dnumdx += 1./(2*radic)*( 2*xydiff*( 2/double(nHit-1)*(y[i]-yAve) )
591  + 8*Sxy*1./double(nHit-1)*(x[i]-xAve) );
592  double dM1dx = denom != 0. ? ( dnumdx - M1*2/double(nHit-1)*(y[i]-yAve) ) / ( denom ) : 0.;
593  double dM1dy = denom != 0. ? ( dnumdy - M1*2/double(nHit-1)*(x[i]-xAve) ) / ( denom ) : 0.;
594  /*
595  double dB1dx = -( M1/double(nHit) + xAve*dM1dx );
596  double dB1dy = 1/double(nHit) - xAve*dM1dy;
597  double dx1dx = 0;
598  double dx2dx = 0;
599  double dy1dx = dM1dx*x1 + M1*dx1dx + dB1dx;
600  double dy1dy = dM1dy*x1 + dB1dy;
601  double dy2dx = dM1dx*x2 + M1*dx2dx + dB1dx;
602  double dy2dy = dM1dy*x2 + dB1dy;
603  */
604  double dy1dx = dM1dx*(x1-xAve) + M1;
605  double dy1dy = dM1dy*(x1-xAve) + 1./double(nHit);
606  double dy2dx = dM1dx*(x2-xAve) + M1;
607  double dy2dy = dM1dy*(x2-xAve) + 1./double(nHit);
608  sig2y1 += dy1dx*dy1dx*sigX*sigX + dy1dy*dy1dy*sigY*sigY;
609  sig2y2 += dy2dx*dy2dx*sigX*sigX + dy2dy*dy2dy*sigY*sigY;
610  }
611  sigy1 = sqrt(sig2y1);
612  sigy2 = sqrt(sig2y2);
613 
614  //if ( sigy1 > 1000 || sigy2 > 1000 ) {
615  // std::cout << "Anomalous solution! ( " << sigy1 << " " << sigy2 << " )" << std::endl;
616  // for ( i = 0; i < w.size(); ++i ) {
617  // if ( w[i] == 0 ) continue;
618  // std::cout << "(" << y[i] << "," << x[i] << ") ";
619  // }
620  // std::cout << "Endpoints: (" << y1 << "," << x1 << ") (" << y2 << "," << x2 << ")" <<std::endl;
621  // std::cout << "Slope " << M1 << " Intercept " << B1 << std::endl;
622  //}
623 
624  // Compute the chi^2 function for return (Doesn't use errors!)
625  double chi2 = 0.;
626 
627  return chi2;
628 }
629 
632  size_t plane = cellHit->Plane();
633  size_t cell = cellHit->Cell();
634  const geo::CellGeo* theCell = geoSvc->Plane(plane)->Cell(cell);
635  double localHitPos[3] = { 0., 0., XY };
636  double worldHitPos[3];
637  theCell->LocalToWorld(localHitPos,worldHitPos);
638  return TVector3(worldHitPos);
639 }
640 
641 void align::SplitTracks::GetCellEndpoints( art::Ptr<rb::CellHit> cellHit, double *cellLo, double *cellHi ) {
643  size_t plane = cellHit->Plane();
644  size_t cell = cellHit->Cell();
645  const geo::CellGeo* theCell = geoSvc->Plane(plane)->Cell(cell);
646  double localLo[3] = { 0., 0., -theCell->HalfL() };
647  double localHi[3] = { 0., 0., theCell->HalfL() };
648  theCell->LocalToWorld(localLo,cellLo);
649  theCell->LocalToWorld(localHi,cellHi);
650  return;
651 }
652 
std::vector< double > subtrk_ErrStopY
void beginJob() override
An algorithm to perform cosmic ray track fitting.
Definition: TrackAlg.h:25
double HalfL() const
Definition: CellGeo.cxx:198
Simple analyzer module to print to text file information required for c++ based block alignment study...
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
std::vector< double > subtrk_StartZ
art::Ptr< rb::CellHit > XCell(unsigned int xIdx) const
Get the ith cell in the x-view.
Definition: Cluster.cxx:157
Float_t y1[n_points_granero]
Definition: compare.C:5
double delta
Definition: runWimpSim.h:98
void LocalToWorld(const double *local, double *world) const
Definition: CellGeo.cxx:80
unsigned short Plane() const
Definition: CellHit.h:39
std::vector< double > subtrk_ErrStartZ
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
double DemingRegFit(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &w, double &sigXY, double &sigZ, double &x1, double &y1, double &x2, double &y2, double &sigy1, double &sigy2)
Find the best-fit line to a collection of points in 2-D using Deming regression.
std::vector< double > blockCenterX
T sqrt(T number)
Definition: d0nt_math.hpp:156
Vertical planes which measure X.
Definition: PlaneGeo.h:28
std::vector< double > subtrk_ErrStopZ
A collection of associated CellHits.
Definition: Cluster.h:47
const daqchannelmap::DAQChannelMap * Map() const
Definition: CMap.h:57
unsigned int MaxCell(geo::View_t view) const
Definition: Cluster.cxx:518
trk::TrackAlg * fTrkAlg
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
Track finder for cosmic rays.
double chi2()
lchan encodeLChan(int detId, plane_t plane, cell_t cell) const
An algorithm to perform cosmic ray track fitting.
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
unsigned short Cell() const
Definition: CellHit.h:40
SplitTracks(fhicl::ParameterSet const &pset)
block
print "ROW IS " print row
Definition: elec2geo.py:31
std::vector< double > subtrk_ErrStartX
art::Ptr< rb::CellHit > YCell(unsigned int yIdx) const
Get the ith cell in the y-view.
Definition: Cluster.cxx:165
T get(std::string const &key) const
Definition: ParameterSet.h:231
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
int evt
std::vector< int > subtrk_Block
std::vector< double > subtrk_StopX
Collect Geo headers and supply basic geometry functions.
void cellPos()
Definition: cellShifts.C:26
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
An algorithm to perform cosmic ray track fitting.
std::vector< double > subtrk_StopY
OStream cout
Definition: OStream.cxx:6
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
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
std::vector< double > subtrk_ErrStartY
T * make(ARGS...args) const
double ClosestApproach(const double point[], const double intercept[], const double slopes[], double closest[])
Find the distance of closest approach between point and line.
Definition: Geo.cxx:222
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
void analyze(art::Event const &evt) override
std::vector< double > blockCenterZ
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
unsigned int MinCell(geo::View_t view) const
Definition: Cluster.cxx:472
std::vector< double > subtrk_StartY
unsigned int ExtentPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.h:250
std::vector< double > subtrk_StopZ
assert(nhit_max >=nhit_nbins)
std::vector< double > subtrk_ErrStopX
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
const int ID() const
Definition: Cluster.h:75
virtual void InterpolateXY(double z, double &x, double &y) const
Definition: Track.cxx:325
std::vector< double > blockCenterY
virtual std::vector< rb::Track > MakeTrack(const rb::Cluster *slice)=0
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
virtual block_t computeBlock(lchan logicalchan) const =0
Which block is this lchan in?
unsigned int fMinPlanesInSubCluster
Encapsulate the cell geometry.
Definition: CellGeo.h:25
Track finder for cosmic rays.
Float_t w
Definition: plot.C:20
#define W(x)
std::vector< double > subtrk_StartX
Encapsulate the geometry of one entire detector (near, far, ndos)
virtual bool Is3D() const
Definition: Prong.h:71
void GetCellEndpoints(art::Ptr< rb::CellHit > cell, double *cellLo, double *cellHi)
TVector3 GeoEstimatePos(art::Ptr< rb::CellHit > cellHit, double W)