Alignment_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief Apply mis-alignment with TranslatePlane() method and then
3 /// Alignment code calculates residuals for tracks for each cell.
4 /// Populates a TH2 using D2M method that can be used to calculate
5 /// alignment correction inter-cell.
6 /// \author gsdavies@iastate.edu
7 ////////////////////////////////////////////////////////////////////////
8 
9 // C/C++ includes
10 #include <cmath>
11 #include <iostream>
12 #include <map>
13 #include <vector>
14 
15 // NOvASoft includes
16 #include "Geometry/Geometry.h"
17 #include "GeometryObjects/Geo.h"
20 #include "RecoBase/CellHit.h"
21 #include "RecoBase/HitMap.h"
22 #include "RecoBase/Track.h"
23 #include "Utilities/func/MathUtil.h"
24 
25 // ROOT includes
26 #include "TH1.h"
27 #include "TH2.h"
28 
29 // ART includes
37 
38 
39 /// Relative alignment of detector cells and blocks
40 namespace align {
41 
42  class Alignment : public art::EDAnalyzer {
43  public:
44  explicit Alignment(fhicl::ParameterSet const& pset);
45  virtual ~Alignment();
46 
47  void analyze(art::Event const& evt);
48  void reconfigure(const fhicl::ParameterSet& pset);
49 
50  void beginJob();
51  void endJob();
52 
53  private:
54 
56  double x,
57  double y,
58  double z); ///< yes/no is "track" entering within a buffer window
59 
60  std::string fCellHitsListLabel; ///< Where to find CellHit
61  std::string fTrackLabel; ///< Where to find Tracks
62 
63  unsigned int fMuonNhit; ///< Required # hits on muon track
64  unsigned int fMuonNhitX; ///< # required in x view
65  unsigned int fMuonNhitY; ///< # required in y view
66  double fDetEdgeBuffer; ///< Detector edge buffer zone for entering cosmic ray track decision [cm]
67 
68  //......................................................................
69  struct Val_t{
70  TH1 *ResidV, *ResidZ;
71  TH2 *HCellPos;
72  };
73  Val_t GetChannelHists(int plane, int cell);
74 
75  /// This is where \ref GetChannelHist stores its histograms
76  std::map<geo::OfflineChan, Val_t> fChannelMap;
77  };
78 }
79 
80 ////////////////////////////////////////////////////////////////////////
81 namespace align{
82 
84  EDAnalyzer(pset),
85  fCellHitsListLabel(pset.get< std::string >("CellHitsListLabel")),
86  fTrackLabel (pset.get< std::string >("TrackLabel") ),
87  fMuonNhit (pset.get< unsigned int >("MuonNhit") ),
88  fMuonNhitX (pset.get< unsigned int >("MuonNhitX") ),
89  fMuonNhitY (pset.get< unsigned int >("MuonNhitY") ),
90  fDetEdgeBuffer (pset.get< double >("DetEdgeBuffer") )
91  {
92  reconfigure(pset);
93  }
94 
95  //......................................................................
97  {
98  }
99 
100  //......................................................................
102  {
103  }
104 
105  //......................................................................
107  {
109 
110  // Get the tracks that were made by CosmicTrack
111  // and test that there are tracks in the event
113  evt.getByLabel(fTrackLabel,trks);
114  if (trks->size()==0) {
115  mf::LogWarning("AlignWarn") << "There are " << trks->size()
116  << " Tracks in this event, skip";
117  return;
118  }
119  LOG_DEBUG("Alignment") << "There are " << trks->size() << " Tracks in this event";
120 
121 
122  const unsigned int trkMax = trks->size();
123 
124  for(unsigned int trkIdx = 0; trkIdx < trkMax; ++trkIdx){
125  art::Ptr<rb::Track> track(trks, trkIdx);
126 
127  TVector3 start(track->Start());
128  TVector3 end(track->Stop());
129  if(end.Y() > start.Y()) std::swap(start, end);
130 
131  // Is this a good track?? i.e. entering and extends across at least 8 planes
132  if (this->IsGoodTrack(track,start.X(),start.Y(),start.Z()) == false) return;
133  assert(track->ExtentPlane() >= 8);
134 
135  const unsigned int hitMax = track->NCell(); // maximum number of CellHits on track
136  const unsigned int hitX = track->NXCell(); // maximum number of CellHits on track in XZ
137  const unsigned int hitY = track->NYCell(); // maximum number of CellHits on track in YZ
138 
139  // Do NOT select tracks if max number of hits or XZ- YZ-hits less than configurable values
140  if ( (hitMax < fMuonNhit) || (hitX < fMuonNhitX) || (hitY < fMuonNhitY) ) return;
141  assert(hitMax >= fMuonNhit);
142 
143  rb::HitMap hmap;
144  // Loop through both views and make a list of every CellHit in this track
145  for(int view = geo::kX; view <= geo::kY; ++view){
146  // C++ is a bit strict about the conversions. Save typing.
147  const geo::View_t geoview = geo::View_t(view);
148  // Make a list of every hit in this track
149  const unsigned int hitMaxView = track->NCell(geoview);
150  for(unsigned int hitIdx = 0; hitIdx < hitMaxView; ++hitIdx){
151  art::Ptr<rb::CellHit> hit = track->Cell(geoview, hitIdx);
152  hmap.Add(hit);
153  } // end for hitIdx
154  } // view
155 
156  const unsigned int pMax = geom->NPlanes();
157 
158  for( unsigned int pIdx = 0; pIdx < pMax; pIdx++){ // loop through planes
159 
160  geo::PlaneGeo* p = new geo::PlaneGeo(*geom->Plane(pIdx));
161  const unsigned int cMax = p->Ncells();
162 
163  // This command is what actually applies the Translation/Rotation.
164  // In this case it applies a 50cm translation in the +z-axis
165  // For now this is the default. Will re-work code so that this is a
166  // fcl parameter input.
167  p->TranslatePlane(0,0,50);
168 
169 
170  for( unsigned int cIdx = 0; cIdx < cMax; cIdx++){ // loop through cells
171  if(hmap.CellExists(pIdx,cIdx)) {
172  double xyz[3];
173  geo::CellGeo* c = new geo::CellGeo(*p->Cell(cIdx));
174 
175  c->GetCenter(xyz, 0.0);
176 
177  const geo::View_t view = geom->Plane(pIdx)->View();
178 
179  double residV[2];
180  double residZ[2];
181  TVector3 closest;
182  double vz,zz;
183 
184  // New format (correct) whereby residual is V(Z0) - V0 so that V0' = V0 + Residual[V]
185  if(view == geo::kX){
186  TVector3 start_vec(start.X(),0.0,start.Z());
187  TVector3 end_vec(end.X(),0.0, end.Z());
188  TVector3 xyz_vec(xyz[0],0.0,xyz[2]);
189  geo::ClosestApproach(xyz_vec,start_vec,(end_vec - start_vec),closest);
190  vz = closest[0];
191  zz = closest[2];
192  residV[view] = vz - xyz[0];
193  residZ[view] = zz - xyz[2];
194  }
195  if(view == geo::kY){
196  TVector3 start_vec(0.0,start.Y(),start.Z());
197  TVector3 end_vec(0.0,end.Y(),end.Z());
198  TVector3 xyz_vec(0.0,xyz[1],xyz[2]);
199  geo::ClosestApproach(xyz_vec,start_vec,(end_vec - start_vec),closest);
200  vz = closest[1];
201  zz = closest[2];
202  residV[view] = vz - xyz[1];
203  residZ[view] = zz - xyz[2];
204  }
205  Val_t hists = GetChannelHists(pIdx, cIdx);
206  hists.ResidV->Fill(residV[view]);
207  hists.ResidZ->Fill(residZ[view]);
208  hists.HCellPos->SetEntries(1);
209  //loop over 2D bins of a histo
210  double bs,bt,d,moda;
211  for(int ibin = 1; ibin <= hists.HCellPos->GetNbinsX(); ibin++){
212  for(int jbin = 1; jbin <= hists.HCellPos->GetNbinsY(); jbin++){
213  bs = hists.HCellPos->GetXaxis()->GetBinCenter(ibin);
214  bt = hists.HCellPos->GetYaxis()->GetBinCenter(jbin);
215 
216  moda = util::pythag(residV[view],residZ[view]);
217  d = moda - ((residV[view]*bs + residZ[view]*bt)/moda);
218  int bin = hists.HCellPos->GetBin(ibin,jbin);
219  hists.HCellPos->AddBinContent(bin,util::sqr(d));
220  } // jbin
221  } // ibin
222  } // if cell exists
223  } // cIdx
224  } // pIdx
225  } // end for trkIdx
226 
227 
228  return;
229  }
230 
231  //......................................................................
233  {
234  }
235  //......................................................................
237  {
238  }
239 
240  //......................................................................
242  {
243  geo::OfflineChan key(plane, cell);
244  std::map<geo::OfflineChan, Val_t>::iterator it = fChannelMap.find(key);
245  if(it != fChannelMap.end()) return it->second;
247  TH1* residv =
248  tfs->make<TH1I>(TString::Format("residualmcV_%03d_%03d", plane, cell),
249  ";Residual in V (cm);count", 1000, -50, 50);
250  TH1* residz =
251  tfs->make<TH1I>(TString::Format("residualmcZ_%03d_%03d", plane, cell),
252  ";Residual in Z (cm);count", 1000, -50, 50);
253  TH2* hCellPos =
254  tfs->make<TH2D>(TString::Format("hCellPos_%03d_%03d", plane, cell),
255  ";Position in Z (cm); Position in Y (cm)", 40, -10, 10, 40, -10, 10);
256 
257  Val_t v;
258  v.ResidV = residv;
259  v.ResidZ = residz;
260  v.HCellPos = hCellPos;
261  fChannelMap[key] = v;
262  return v;
263  }
264 
265  //.....................................................................
266  bool Alignment::IsGoodTrack(art::Ptr<rb::Track> t, double x, double y, double z)
267  {
269  if (t->ExtentPlane() >= 8) return true;
270  else return false;
271  if (y > geom->DetHalfHeight()-fDetEdgeBuffer) return true;
272  if (fabs(x) > geom->DetHalfWidth()-fDetEdgeBuffer) return true;
273  if (z < fDetEdgeBuffer) return true;
274  if (z > geom->DetLength()-fDetEdgeBuffer) return true;
275  return false;
276  }
277 } //namespace align
278 
279 ////////////////////////////////////////////////////////////////////////
281 ////////////////////////////////////////////////////////////////////////
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
Simple analyzer module to print to text file information required for c++ based block alignment study...
std::string fCellHitsListLabel
Where to find CellHit.
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
set< int >::iterator it
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
TH2D * hCellPos[pmax][cmax]
Definition: cellShifts.C:21
Alignment(fhicl::ParameterSet const &pset)
std::map< geo::OfflineChan, Val_t > fChannelMap
This is where GetChannelHist stores its histograms.
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
Vertical planes which measure X.
Definition: PlaneGeo.h:28
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
double DetLength() const
unsigned int fMuonNhitX
required in x view
Double_t zz
Definition: plot.C:277
const PlaneGeo * Plane(unsigned int i) const
TString hists[nhists]
Definition: bdt_com.C:3
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
Provides efficient lookup of CellHits by plane and cell number.
Definition: HitMap.h:22
void analyze(art::Event const &evt)
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
std::string fTrackLabel
Where to find Tracks.
double fDetEdgeBuffer
Detector edge buffer zone for entering cosmic ray track decision [cm].
void reconfigure(const fhicl::ParameterSet &pset)
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
Geometry information for a single readout plane.
Definition: PlaneGeo.h:36
Collect Geo headers and supply basic geometry functions.
Float_t d
Definition: plot.C:236
bool IsGoodTrack(art::Ptr< rb::Track > t, double x, double y, double z)
yes/no is "track" entering within a buffer window
double DetHalfHeight() const
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
float bin[41]
Definition: plottest35.C:14
z
Definition: test.py:28
void TranslatePlane(double dx, double dy=0, double dz=0, double phi=0, double theta=0, double psi=0)
Definition: PlaneGeo.cxx:179
Val_t GetChannelHists(int plane, int cell)
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
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
double DetHalfWidth() 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
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
Definition: event.h:1
A (plane, cell) pair.
Definition: OfflineChan.h:17
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
void geom(int which=0)
Definition: geom.C:163
unsigned int ExtentPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.h:250
assert(nhit_max >=nhit_nbins)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
bool CellExists(unsigned int planeIdx, unsigned int cellIdx) const
Does the map contain any cell at this position?
Definition: HitMap.cxx:81
unsigned int NPlanes() const
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
Encapsulate the cell geometry.
Definition: CellGeo.h:25
unsigned int fMuonNhitY
required in y view
Float_t track
Definition: plot.C:35
Simple object representing a (plane, cell) pair.
Encapsulate the geometry of one entire detector (near, far, ndos)
void Add(const art::Ptr< rb::CellHit > &chit)
Definition: HitMap.cxx:37
unsigned int fMuonNhit
Required # hits on muon track.