disambigcheck_module.cc
Go to the documentation of this file.
1 #ifndef disambigcheck_h
2 #define disambigcheck_h
3 ////////////////////////////////////////////////////////////////////////
4 //
5 // Disambigcheck class
6 //
7 // tjyang@fnal.gov
8 //
9 // Based on Tom Junk's idea of 3-view matching
10 // HitFinder for 35t, input is undisambiguated hits, out is disambiguated hits
11 // Start from code written by talion@gmail.com
12 //
13 //
14 ////////////////////////////////////////////////////////////////////////
15 
16 extern "C" {
17 #include <sys/types.h>
18 #include <sys/stat.h>
19 }
20 #include <stdint.h>
21 #include <string>
22 
23 // Framework includes
27 #include "art_root_io/TFileService.h"
31 #include "fhiclcpp/ParameterSet.h"
32 
33 // LArSoft Includes
42 #include "DisambigAlg35t.h"
46 
47 // ROOT Includes
48 #include "TH1D.h"
49 #include "TH2D.h"
50 #include "TDecompSVD.h"
51 #include "TMath.h"
52 #include "TF1.h"
53 #include "TTree.h"
54 #include "TVectorD.h"
55 #include "TVector2.h"
56 #include "TVector3.h"
57 #include "TCanvas.h"
58 
59 namespace disambigcheck{
61  {
62 
63  public:
64 
65  explicit disambigcheck(fhicl::ParameterSet const& pset);
66 
67  void analyze(const art::Event& evt);
68  void beginJob();
69  void beginRun(const art::Run& run);
70  void endJob();
71  void reconfigure(fhicl::ParameterSet const& pset);
72 
73  virtual ~disambigcheck();
74 
75  private:
76 
80 
83 
84  TH1D* fCorrect; TH1D* fMissed; TH1D* fIncorrect;
85  TH2D* fXCorrect; TH2D* fXMissed; TH2D* fXIncorrect;
86  TH2D* fYCorrect; TH2D* fYMissed; TH2D* fZIncorrect;
87  TH2D* fZCorrect; TH2D* fZMissed; TH2D* fYIncorrect;
88  TH2D* fThetaCorrect; TH2D* fThetaMissed; TH2D* fThetaIncorrect;
89  TH2D* fPhiCorrect; TH2D* fPhiMissed; TH2D* fPhiIncorrect;
90 
91  //TCanvas* cXCorrect;
92 
93  protected:
94 
95  };
96 
97  //-------------------------------------------------
98  //-------------------------------------------------
100  : EDAnalyzer (pset)
101  {
102  this->reconfigure(pset);
103  // produces< std::vector<recob::Hit> >();
104  }
105 
106  disambigcheck::~disambigcheck() {}
107  //-------------------------------------------------
108  //-------------------------------------------------
110  {
111 
112  fChanHitDisambig = pset.get< std::string >("ChanHitDisambig");
113  fChanHitCheater = pset.get< std::string >("ChanHitCheater");
114 
115  }
116 
117  //-------------------------------------------------
118  //-------------------------------------------------
120  {
122 
123  fCorrect = tfs->make<TH1D>("correct" ,"Correct disambiguated" , 100, 0, 1.01);
124  fMissed = tfs->make<TH1D>("missed" ,"Missed hits" , 100, 0, 1.01);
125  fIncorrect = tfs->make<TH1D>("incorrect","Incorrect disambiguated", 100, 0, 1.01);
126 
127  fXCorrect = tfs->make<TH2D>("Xcorrect" ,"Correct disambiguated" , 40, -1.01, 1.01, 40, 0, 1.01);
128  fXMissed = tfs->make<TH2D>("Xmissed" ,"Missed hits" , 40, -1.01, 1.01, 40, 0, 1.01);
129  fXIncorrect = tfs->make<TH2D>("Xincorrect","Incorrect disambiguated", 40, -1.01, 1.01, 40, 0, 1.01);
130  fYCorrect = tfs->make<TH2D>("Ycorrect" ,"Correct disambiguated" , 40, -1.01, 1.01, 40, 0, 1.01);
131  fYMissed = tfs->make<TH2D>("Ymissed" ,"Missed hits" , 40, -1.01, 1.01, 40, 0, 1.01);
132  fYIncorrect = tfs->make<TH2D>("Yincorrect","Incorrect disambiguated", 40, -1.01, 1.01, 40, 0, 1.01);
133  fZCorrect = tfs->make<TH2D>("Zcorrect" ,"Correct disambiguated" , 40, -1.01, 1.01, 40, 0, 1.01);
134  fZMissed = tfs->make<TH2D>("Zmissed" ,"Missed hits" , 40, -1.01, 1.01, 40, 0, 1.01);
135  fZIncorrect = tfs->make<TH2D>("Zincorrect","Incorrect disambiguated", 40, -1.01, 1.01, 40, 0, 1.01);
136 
137  fThetaCorrect = tfs->make<TH2D>("Thetacorrect" ,"Correct disambiguated" , 20, 0 , 3.15, 40, 0, 1.01);
138  fThetaMissed = tfs->make<TH2D>("Thetamissed" ,"Missed hits" , 20, 0 , 3.15, 40, 0, 1.01);
139  fThetaIncorrect = tfs->make<TH2D>("Thetaincorrect","Incorrect disambiguated", 20, 0 , 3.15, 40, 0, 1.01);
140  fPhiCorrect = tfs->make<TH2D>("Phicorrect" ,"Correct disambiguated" , 40, -3.15, 3.15, 40, 0, 1.01);
141  fPhiMissed = tfs->make<TH2D>("Phimissed" ,"Missed hits" , 40, -3.15, 3.15, 40, 0, 1.01);
142  fPhiIncorrect = tfs->make<TH2D>("Phiincorrect" ,"Incorrect disambiguated", 40, -3.15, 3.15, 40, 0, 1.01);
143 
144  //cXCorrect = tfs->make<TCanvas>("XCorrect", "", 600, 500);
145 }
146 
147  //-------------------------------------------------
148  //-------------------------------------------------
149  void disambigcheck::beginRun(const art::Run& run){}
150  void disambigcheck::endJob()
151  {
152  //cXCorrect -> cd(0);
153  //fXCorrect -> Draw("colz");
154  //cXCorrect -> Write();
155  }
156 
157 
158  //-------------------------------------------------
160  {
161 
162  std::unique_ptr<std::vector<recob::Hit> > hcol(new std::vector<recob::Hit>);
163  std::vector< art::Ptr<recob::Hit> > ChHitsDisambig;
164  auto ChannelHitsDisambig = evt.getHandle< std::vector<recob::Hit> >(fChanHitDisambig);
165  if (ChannelHitsDisambig)
166  art::fill_ptr_vector(ChHitsDisambig, ChannelHitsDisambig);
167 
168  // Make unambiguous collection hits
169  std::vector< art::Ptr<recob::Hit> > ChHitsCheater;
170  auto ChannelHitsCheater = evt.getHandle< std::vector<recob::Hit> >(fChanHitCheater);
171  if (ChannelHitsCheater)
172  art::fill_ptr_vector(ChHitsCheater, ChannelHitsCheater);
173 
174  int correcthits = 0;
175  int incorrecthits = 0;
176  int missinghits = 0;
177  // Map of MCParticle Track IDs, Vector 1 ( matched, incorrect, missed ), Pair 1 ( Pair 2 (Theta, Phi) , Vector 2 direction cosines ( x, y, z ) )
178  std::map< int, std::pair< TVector3, std::pair < std::pair< double,double >, TVector3 > > > ParticleMap;
179  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
180  if (ChannelHitsDisambig.isValid()) {
181  for( size_t h = 0; h < ChHitsCheater.size(); h++ ){
182  if(ChHitsCheater[h]->View() == geo::kZ ) continue;
183  uint32_t cheatchannel = ChHitsCheater[h]->Channel();
184  double cheatpeaktime = ChHitsCheater[h]->PeakTime();
185  uint32_t cheatwire = ChHitsCheater[h]->WireID().Wire;
186 
187  std::vector<sim::IDE> ides;
188  std::vector<sim::TrackIDE> TrackIDs = bt_serv->HitToTrackIDEs(clockData, ChHitsCheater[h]);
189  double MaxE = 0;
190  int TrackID;
191  for(size_t e = 0; e < TrackIDs.size(); ++e){
192  if ( TrackIDs[e].energy > MaxE ) {
193  TrackID = TrackIDs[e].trackID;
194  MaxE = TrackIDs[e].energy;
195  }
196  }
197  // Now have trackID, so get PdG code and T0 etc.
198  const simb::MCParticle *particle = pi_serv->TrackIdToParticle_P(TrackID);
199  if ( ParticleMap.count(TrackID) == 0) {
200  //std::cout << "Looking at key " << TrackID << " for the first time so setting the angles." << std::endl;
201  TVector3 MC_NormMom;
202  MC_NormMom[0] = ( particle->Px() / particle->P() );
203  MC_NormMom[1] = ( particle->Py() / particle->P() );
204  MC_NormMom[2] = ( particle->Pz() / particle->P() );
205  double MC_Theta = MC_NormMom.Theta();
206  double MC_Phi = MC_NormMom.Phi();
207  ParticleMap[TrackID].second.second[0] = MC_NormMom[0];
208  ParticleMap[TrackID].second.second[1] = MC_NormMom[1];
209  ParticleMap[TrackID].second.second[2] = MC_NormMom[2];
210  ParticleMap[TrackID].second.first.first = MC_Theta;
211  ParticleMap[TrackID].second.first.second = MC_Phi;
212  } // If first instance of this trackID fill quantities
213  bool found = false;
214  for( size_t w = 0; w < ChHitsDisambig.size(); w++ ) {
215  uint32_t disambigchannel = ChHitsDisambig[w]->Channel();
216  double disambigpeaktime = ChHitsDisambig[w]->PeakTime();
217  uint32_t disambigwire = ChHitsDisambig[w]->WireID().Wire;
218 
219  if((cheatchannel == disambigchannel)&&(TMath::Abs(cheatpeaktime - disambigpeaktime)< 0.1)) {
220 
221  found=true;
222  if(cheatwire==disambigwire){
223  correcthits++;
224  ++ParticleMap[TrackID].first[0];
225  }
226  else{
227  incorrecthits++;
228  ++ParticleMap[TrackID].first[1];
229  }
230  }
231  } // Loop through disambiguated hits
232  if(!found){
233  missinghits++;
234  ++ParticleMap[TrackID].first[2];
235  }
236  //std::cout << "Particle Map Angles...." << ParticleMap[TrackID].first[0] << " " << ParticleMap[TrackID].first[1] << " " << ParticleMap[TrackID].first[2] <<"\n"
237  // << "And it now has " << ParticleMap[TrackID].second[0] << " correct, " << ParticleMap[TrackID].second[1] << " incorrect and " << ParticleMap[TrackID].second[2] << " missed hits." << std::endl;
238  } // size of cheated hits
239 
240  int totalhits = (correcthits + incorrecthits + missinghits);
241 
242  //double disambighitsFraction = ddisambighits / dtotalhits;
243  if ( totalhits != 0 ) {
244  fCorrect ->Fill( (double)correcthits / (double)totalhits );
245  fMissed ->Fill( (double)missinghits / (double)totalhits );
246  fIncorrect->Fill( (double)incorrecthits / (double)totalhits );
247  if ( (double)correcthits / (double)totalhits == 0 ) std::cout << "WHY IS THIS EVENT NOT MATCHING ANY HITS!!!????" << std::endl;
248  } // totalhits != 0
249  std::cout << "\nTotal hits " << totalhits << ", correct hits " << correcthits << ", incorrect hits " << incorrecthits << ", missing hits " << missinghits << std::endl;
250 
251  for (std::map< int, std::pair< TVector3, std::pair < std::pair< double,double >, TVector3 > > >::iterator ii = ParticleMap.begin(); ii!=ParticleMap.end(); ++ii){
252  int MapTotalHits = ii->second.first[0] + ii->second.first[1] + ii->second.first[2];
253  std::cout << "Total hits for this TrackID " << ii->first << " is; " << MapTotalHits << ", Correct " << ii->second.first[0] << ", Incorrect " << ii->second.first[1] << ", Missed " << ii->second.first[2] << std::endl;
254  if ( MapTotalHits != 0 ) {
255  fXCorrect -> Fill ( ii->second.second.second[0], (double)ii->second.first[0] / (double)MapTotalHits );
256  fXMissed -> Fill ( ii->second.second.second[0], (double)ii->second.first[1] / (double)MapTotalHits );
257  fXIncorrect -> Fill ( ii->second.second.second[0], (double)ii->second.first[2] / (double)MapTotalHits );
258  fYCorrect -> Fill ( ii->second.second.second[1], (double)ii->second.first[0] / (double)MapTotalHits );
259  fYMissed -> Fill ( ii->second.second.second[1], (double)ii->second.first[1] / (double)MapTotalHits );
260  fYIncorrect -> Fill ( ii->second.second.second[1], (double)ii->second.first[2] / (double)MapTotalHits );
261  fZCorrect -> Fill ( ii->second.second.second[2], (double)ii->second.first[0] / (double)MapTotalHits );
262  fZMissed -> Fill ( ii->second.second.second[2], (double)ii->second.first[1] / (double)MapTotalHits );
263  fZIncorrect -> Fill ( ii->second.second.second[2], (double)ii->second.first[2] / (double)MapTotalHits );
264 
265  fThetaCorrect -> Fill ( ii->second.second.first.first , (double)ii->second.first[0] / (double)MapTotalHits );
266  fThetaMissed -> Fill ( ii->second.second.first.first , (double)ii->second.first[1] / (double)MapTotalHits );
267  fThetaIncorrect -> Fill ( ii->second.second.first.first , (double)ii->second.first[2] / (double)MapTotalHits );
268  fPhiCorrect -> Fill ( ii->second.second.first.second, (double)ii->second.first[0] / (double)MapTotalHits );
269  fPhiMissed -> Fill ( ii->second.second.first.second, (double)ii->second.first[1] / (double)MapTotalHits );
270  fPhiIncorrect -> Fill ( ii->second.second.first.second, (double)ii->second.first[2] / (double)MapTotalHits );
271 
272  //Some cout statements to have a look at......
273  /*
274  if ( (double)ii->second.first[0] / (double)MapTotalHits < 0.8 ) {
275  std::cout << "\n\nThis track has efficiency " << (double)ii->second.first[0] / (double)MapTotalHits*100
276  << "\nTheta is " << ii->second.second.first.first << ", and phi is " << ii->second.second.first.second
277  << std::endl;
278  if ( ii->second.second.first.first < 1.7 && ii->second.second.first.first > 1.45 ) {
279  std::cout << "THIS TRACK HAS LOWER EFFICIENCY AT THETA IS ROUGHLY PI/2." << std::endl;
280  }
281  if ( ii->second.second.first.second < -1.45 && ii->second.second.first.second > -1.7 ) {
282  std::cout << "THIS TRACK HAS LOWER EFFICIENCY AT PHI IS ROUGHLY -PI/2." << std::endl;
283  }
284  std::cout <<"\n\n" << std::endl;
285  }//*/
286  } // MapTotalHits != 0
287  } // Loop through Map
288  } // If hits valid
289  ParticleMap.clear();
290  }
291 
293 
294 }
295 
296 
297 
298 #endif
def analyze(root, level, gtrees, gbranches, doprint)
Definition: rootstat.py:69
intermediate_table::iterator iterator
double Py(const int i=0) const
Definition: MCParticle.h:231
Encapsulate the construction of a single cyostat.
AdcChannelData::View View
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
std::string string
Definition: nybbler.cc:12
std::vector< TrackID > TrackIDs
double Px(const int i=0) const
Definition: MCParticle.h:230
Planes which measure Z direction.
Definition: geo_types.h:132
art framework interface to geometry description
Definition: Run.h:17
art::ServiceHandle< cheat::BackTrackerService > bt_serv
const double e
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
void beginJob()
Definition: Breakpoints.cc:14
virtual void reconfigure(fhicl::ParameterSet const &pset)
double P(const int i=0) const
Definition: MCParticle.h:234
T get(std::string const &key) const
Definition: ParameterSet.h:271
art::ServiceHandle< geo::Geometry > fGeom
Declaration of signal hit object.
Encapsulate the construction of a single detector plane.
double Pz(const int i=0) const
Definition: MCParticle.h:232
TCEvent evt
Definition: DataStructs.cxx:7
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
art::ServiceHandle< cheat::ParticleInventoryService > pi_serv
QTextStream & endl(QTextStream &s)
Encapsulate the construction of a single detector plane.