1 #ifndef disambigcheck_h 2 #define disambigcheck_h 17 #include <sys/types.h> 27 #include "art_root_io/TFileService.h" 50 #include "TDecompSVD.h" 84 TH1D* fCorrect; TH1D*
fMissed; TH1D* fIncorrect;
85 TH2D* fXCorrect; TH2D*
fXMissed; TH2D* fXIncorrect;
87 TH2D* fZCorrect; TH2D*
fZMissed; TH2D* fYIncorrect;
89 TH2D* fPhiCorrect; TH2D*
fPhiMissed; TH2D* fPhiIncorrect;
106 disambigcheck::~disambigcheck() {}
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);
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);
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);
150 void disambigcheck::endJob()
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)
169 std::vector< art::Ptr<recob::Hit> > ChHitsCheater;
170 auto ChannelHitsCheater = evt.
getHandle< std::vector<recob::Hit> >(fChanHitCheater);
171 if (ChannelHitsCheater)
175 int incorrecthits = 0;
178 std::map< int, std::pair< TVector3, std::pair < std::pair< double,double >, TVector3 > > > ParticleMap;
180 if (ChannelHitsDisambig.isValid()) {
181 for(
size_t h = 0;
h < ChHitsCheater.size();
h++ ){
183 uint32_t cheatchannel = ChHitsCheater[
h]->Channel();
184 double cheatpeaktime = ChHitsCheater[
h]->PeakTime();
185 uint32_t cheatwire = ChHitsCheater[
h]->WireID().Wire;
187 std::vector<sim::IDE> ides;
188 std::vector<sim::TrackIDE>
TrackIDs = bt_serv->HitToTrackIDEs(clockData, ChHitsCheater[
h]);
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;
199 if ( ParticleMap.count(TrackID) == 0) {
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;
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;
219 if((cheatchannel == disambigchannel)&&(TMath::Abs(cheatpeaktime - disambigpeaktime)< 0.1)) {
222 if(cheatwire==disambigwire){
224 ++ParticleMap[
TrackID].first[0];
228 ++ParticleMap[
TrackID].first[1];
234 ++ParticleMap[
TrackID].first[2];
240 int totalhits = (correcthits + incorrecthits + missinghits);
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;
249 std::cout <<
"\nTotal hits " << totalhits <<
", correct hits " << correcthits <<
", incorrect hits " << incorrecthits <<
", missing hits " << missinghits <<
std::endl;
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 );
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 );
def analyze(root, level, gtrees, gbranches, doprint)
double Py(const int i=0) const
Encapsulate the construction of a single cyostat.
AdcChannelData::View View
Handle< PROD > getHandle(SelectorBase const &) const
std::vector< TrackID > TrackIDs
double Px(const int i=0) const
Planes which measure Z direction.
art framework interface to geometry description
art::ServiceHandle< cheat::BackTrackerService > bt_serv
#define DEFINE_ART_MODULE(klass)
virtual void reconfigure(fhicl::ParameterSet const &pset)
double P(const int i=0) const
T get(std::string const &key) const
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
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
art::ServiceHandle< cheat::ParticleInventoryService > pi_serv
std::string fChanHitDisambig
QTextStream & endl(QTextStream &s)
Encapsulate the construction of a single detector plane.
std::string fChanHitCheater