analysis/SvtxSimPerformanceCheckReco/SvtxSimPerformanceCheckReco.C

0001
0002 #include "SvtxSimPerformanceCheckReco.h"
0003
0004 #include <phool/getClass.h>
0005 #include <fun4all/Fun4AllServer.h>
0006
0007 #include <g4main/PHG4TruthInfoContainer.h>
0008 #include <g4main/PHG4Particle.h>
0009 #include <g4main/PHG4VtxPoint.h>
0010
0011 #include <g4hough/SvtxVertexMap.h>
0012 #include <g4hough/SvtxVertex.h>
0013 #include <g4hough/SvtxTrackMap.h>
0014 #include <g4hough/SvtxTrack.h>
0015
0016 #include <g4eval/SvtxEvalStack.h>
0017 #include <g4eval/SvtxTrackEval.h>
0018 #include <g4eval/SvtxVertexEval.h>
0019 #include <g4eval/SvtxTruthEval.h>
0020
0021 #include <TH1D.h>
0022 #include <TH2D.h>
0023
0024 #include <iostream>
0025
0026 using namespace std;
0027
0028 SvtxSimPerformanceCheckReco::SvtxSimPerformanceCheckReco(const string &name)
0029   : SubsysReco(name),
0030     _event(0),
0031     _nlayers(7),
0032     _inner_layer_mask(0x7),
0033     _truept_dptoverpt(NULL),
0034     _truept_dca(NULL),
0035     _truept_particles_leaving7Hits(NULL),
0036     _truept_particles_recoWithExactHits(NULL),
0037     _truept_particles_recoWithin1Hit(NULL),
0038     _truept_particles_recoWithin2Hits(NULL),
0039     _truept_particles_recoWithExactInnerHits(NULL),
0040     _truept_particles_recoWithin1InnerHit(NULL),
0041     _truept_particles_recoWithin2InnerHits(NULL),
0042     _truept_particles_recoWithin3Percent(NULL),
0043     _truept_particles_recoWithin4Percent(NULL),
0044     _truept_particles_recoWithin5Percent(NULL),
0045     _recopt_tracks_all(NULL),
0046     _recopt_tracks_recoWithExactHits(NULL),
0047     _recopt_tracks_recoWithin1Hit(NULL),
0048     _recopt_tracks_recoWithin2Hits(NULL),
0049     _recopt_tracks_recoWithExactInnerHits(NULL),
0050     _recopt_tracks_recoWithin1InnerHit(NULL),
0051     _recopt_tracks_recoWithin2InnerHits(NULL),
0052     _recopt_tracks_recoWithin3Percent(NULL),
0053     _recopt_tracks_recoWithin4Percent(NULL),
0054     _recopt_tracks_recoWithin5Percent(NULL),
0055     _recopt_quality(NULL),
0056     _dx_vertex(NULL),
0057     _dy_vertex(NULL),
0058     _dz_vertex(NULL),
0059     _truept_quality_particles_recoWithin4Percent(NULL),
0060     _recopt_quality_tracks_all(NULL),
0061     _recopt_quality_tracks_recoWithin4Percent(NULL) {
0062 }
0063
0064 int SvtxSimPerformanceCheckReco::Init(PHCompositeNode *topNode) {
0065
0066   // register histograms
0067   Fun4AllServer *se = Fun4AllServer::instance();
0068
0069   _truept_dptoverpt = new TH2D("truept_dptoverpt",
0070                    "truept_dptoverpt",
0071                    40,0.0,40.0,
0072                    200,-0.5,0.5);
0073
0074   _truept_dca = new TH2D("truept_dca",
0075              "truept_dca",
0076              20,0.0,10.0,
0077              200,-0.1,0.1);
0078
0079   _truept_particles_leaving7Hits = new TH1D("truept_particles_leaving7Hits",
0080                         "truept_particles_leaving7Hits",
0081                         20,0.0,10.0);
0082
0083   _truept_particles_recoWithExactHits = new TH1D("truept_particles_recoWithExactHits",
0084                          "truept_particles_recoWithExactHits",
0085                          20,0.0,10.0);
0086
0087   _truept_particles_recoWithin1Hit = new TH1D("truept_particles_recoWithin1Hit",
0088                           "truept_particles_recoWithin1Hit",
0089                           20,0.0,10.0);
0090
0091   _truept_particles_recoWithin2Hits = new TH1D("truept_particles_recoWithin2Hits",
0092                            "truept_particles_recoWithin2Hits",
0093                            20,0.0,10.0);
0094
0095   _truept_particles_recoWithExactInnerHits = new TH1D("truept_particles_recoWithExactInnerHits",
0096                               "truept_particles_recoWithExactInnerHits",
0097                          20,0.0,10.0);
0098
0099   _truept_particles_recoWithin1InnerHit = new TH1D("truept_particles_recoWithin1InnerHit",
0100                            "truept_particles_recoWithin1InnerHit",
0101                            20,0.0,10.0);
0102
0103   _truept_particles_recoWithin2InnerHits = new TH1D("truept_particles_recoWithin2InnerHits",
0104                             "truept_particles_recoWithin2InnerHits",
0105                             20,0.0,10.0);
0106
0107
0108   _truept_particles_recoWithin3Percent = new TH1D("truept_particles_recoWithin3Percent",
0109                           "truept_particles_recoWithin3Percent",
0110                           20,0.0,10.0);
0111
0112   _truept_particles_recoWithin4Percent = new TH1D("truept_particles_recoWithin4Percent",
0113                           "truept_particles_recoWithin4Percent",
0114                           20,0.0,10.0);
0115
0116   _truept_particles_recoWithin5Percent = new TH1D("truept_particles_recoWithin5Percent",
0117                           "truept_particles_recoWithin5Percent",
0118                           20,0.0,10.0);
0119
0120   _recopt_tracks_all = new TH1D("recopt_tracks_all",
0121                 "recopt_tracks_all",
0122                 20,0.0,10.0);
0123
0124   _recopt_tracks_recoWithExactHits = new TH1D("recopt_tracks_recoWithExactHits",
0125                           "recopt_tracks_recoWithExactHits",
0126                           20,0.0,10.0);
0127
0128   _recopt_tracks_recoWithin1Hit = new TH1D("recopt_tracks_recoWithin1Hit",
0129                        "recopt_tracks_recoWithin1Hit",
0130                        20,0.0,10.0);
0131
0132   _recopt_tracks_recoWithin2Hits = new TH1D("recopt_tracks_recoWithin2Hits",
0133                         "recopt_tracks_recoWithin2Hits",
0134                         20,0.0,10.0);
0135
0136   _recopt_tracks_recoWithExactInnerHits = new TH1D("recopt_tracks_recoWithExactInnerHits",
0137                           "recopt_tracks_recoWithExactInnerHits",
0138                           20,0.0,10.0);
0139
0140   _recopt_tracks_recoWithin1InnerHit = new TH1D("recopt_tracks_recoWithin1InnerHit",
0141                        "recopt_tracks_recoWithin1InnerHit",
0142                        20,0.0,10.0);
0143
0144   _recopt_tracks_recoWithin2InnerHits = new TH1D("recopt_tracks_recoWithin2InnerHits",
0145                         "recopt_tracks_recoWithin2InnerHits",
0146                         20,0.0,10.0);
0147
0148   _recopt_tracks_recoWithin3Percent = new TH1D("recopt_tracks_recoWithin3Percent",
0149                         "recopt_tracks_recoWithin3Percent",
0150                         20,0.0,10.0);
0151
0152   _recopt_tracks_recoWithin4Percent = new TH1D("recopt_tracks_recoWithin4Percent",
0153                            "recopt_tracks_recoWithin4Percent",
0154                            20,0.0,10.0);
0155
0156   _recopt_tracks_recoWithin5Percent = new TH1D("recopt_tracks_recoWithin5Percent",
0157                            "recopt_tracks_recoWithin5Percent",
0158                            20,0.0,10.0);
0159
0160   _recopt_quality = new TH2D("recopt_quality",
0161                  "recopt_quality",
0162                  20,0.0,10.0,
0163                  100,0.0,5.0);
0164
0165   _dx_vertex = new TH1D("dx_vertex",
0166             "dx_vertex",
0167             200,-0.01,0.01);
0168
0169   _dy_vertex = new TH1D("dy_vertex",
0170             "dy_vertex",
0171             200,-0.01,0.01);
0172
0173   _dz_vertex = new TH1D("dz_vertex",
0174             "dz_vertex",
0175             200,-0.01,0.01);
0176
0177   _truept_quality_particles_recoWithin4Percent = new TH2D("truept_quality_particles_recoWithin4Percent",
0178                               "truept_quality_particles_recoWithin4Percent",
0179                               20,0.0,10.0,
0180                               100,0.0,5.0);
0181
0182   _recopt_quality_tracks_all = new TH2D("recopt_quality_tracks_all",
0183                     "recopt_quality_tracks_all",
0184                     20,0.0,10.0,
0185                     100,0.0,5.0);
0186
0187   _recopt_quality_tracks_recoWithin4Percent = new TH2D("recopt_quality_tracks_recoWithin4Percent",
0188                                "recopt_quality_tracks_recoWithin4Percent",
0189                                20,0.0,10.0,
0190                                100,0.0,5.0);
0191
0192
0193   se->registerHisto(_truept_dptoverpt);
0194
0195   se->registerHisto(_truept_dca);
0196
0197   se->registerHisto(_truept_particles_leaving7Hits);
0198
0199   se->registerHisto(_truept_particles_recoWithExactHits);
0200   se->registerHisto(_truept_particles_recoWithin1Hit);
0201   se->registerHisto(_truept_particles_recoWithin2Hits);
0202
0203   se->registerHisto(_truept_particles_recoWithExactInnerHits);
0204   se->registerHisto(_truept_particles_recoWithin1InnerHit);
0205   se->registerHisto(_truept_particles_recoWithin2InnerHits);
0206
0207   se->registerHisto(_truept_particles_recoWithin3Percent);
0208   se->registerHisto(_truept_particles_recoWithin4Percent);
0209   se->registerHisto(_truept_particles_recoWithin5Percent);
0210
0211   se->registerHisto(_recopt_tracks_all);
0212
0213   se->registerHisto(_recopt_tracks_recoWithExactHits);
0214   se->registerHisto(_recopt_tracks_recoWithin1Hit);
0215   se->registerHisto(_recopt_tracks_recoWithin2Hits);
0216
0217   se->registerHisto(_recopt_tracks_recoWithExactInnerHits);
0218   se->registerHisto(_recopt_tracks_recoWithin1InnerHit);
0219   se->registerHisto(_recopt_tracks_recoWithin2InnerHits);
0220
0221   se->registerHisto(_recopt_tracks_recoWithin3Percent);
0222   se->registerHisto(_recopt_tracks_recoWithin4Percent);
0223   se->registerHisto(_recopt_tracks_recoWithin5Percent);
0224
0225   se->registerHisto(_recopt_quality);
0226
0227   se->registerHisto(_dx_vertex);
0228   se->registerHisto(_dy_vertex);
0229   se->registerHisto(_dz_vertex);
0230
0231   se->registerHisto(_truept_quality_particles_recoWithin4Percent);
0232   se->registerHisto(_recopt_quality_tracks_all);
0233   se->registerHisto(_recopt_quality_tracks_recoWithin4Percent);
0234
0235   return 0;
0236 }
0237
0238 int SvtxSimPerformanceCheckReco::process_event(PHCompositeNode *topNode) {
0239
0240   ++_event;
0241
0242   // need things off of the DST...
0243   PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode,"G4TruthInfo");
0244   if (!truthinfo) {
0245     cerr << PHWHERE << " ERROR: Can't find G4TruthInfo" << endl;
0246     exit(-1);
0247   }
0248
0249   SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode,"SvtxTrackMap");
0250   if (!trackmap) {
0251     cerr << PHWHERE << " ERROR: Can't find SvtxTrackMap" << endl;
0252     exit(-1);
0253   }
0254
0255   SvtxVertexMap* vertexmap = findNode::getClass<SvtxVertexMap>(topNode,"SvtxVertexMap");
0256   if (!vertexmap) {
0257     cerr << PHWHERE << " ERROR: Can't find SvtxVertexMap" << endl;
0258     exit(-1);
0259   }
0260
0261   // create SVTX eval stack
0262   SvtxEvalStack svtxevalstack(topNode);
0263
0264   SvtxVertexEval*   vertexeval = svtxevalstack.get_vertex_eval();
0265   SvtxTrackEval*     trackeval = svtxevalstack.get_track_eval();
0266   SvtxTruthEval*     trutheval = svtxevalstack.get_truth_eval();
0267
0268   // loop over all truth particles
0269   PHG4TruthInfoContainer::Range range = truthinfo->GetPrimaryParticleRange();
0270   for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
0271        iter != range.second;
0272        ++iter) {
0273     PHG4Particle* g4particle = iter->second;
0274     if (trutheval->get_embed(g4particle) <= 0) continue;
0275
0276     std::set<PHG4Hit*> g4hits = trutheval->all_truth_hits(g4particle);
0277     float ng4hits = g4hits.size();
0278
0279     float truept = sqrt(pow(g4particle->get_px(),2)+pow(g4particle->get_py(),2));
0280
0281     // examine truth particles that leave 7 detector hits
0282     if (ng4hits == _nlayers) {
0283       _truept_particles_leaving7Hits->Fill(truept);
0284
0285       SvtxTrack* track = trackeval->best_track_from(g4particle);
0286
0287       if (!track) {continue;}
0288
0289       unsigned int nfromtruth = trackeval->get_nclusters_contribution(track,g4particle);
0290       float recopt = track->get_pt();
0291
0292       unsigned int ndiff = abs((int)nfromtruth-(int)_nlayers);
0293       if (ndiff <= 2) {
0294     _truept_particles_recoWithin2Hits->Fill(truept);
0295       }
0296       if (ndiff <= 1) {
0297     _truept_particles_recoWithin1Hit->Fill(truept);
0298       }
0299       if (ndiff == 0) {
0300     _truept_particles_recoWithExactHits->Fill(truept);
0301       }
0302
0303       unsigned int layersfromtruth = trackeval->get_nclusters_contribution_by_layer(track,g4particle);
0304       unsigned int innerhits = (layersfromtruth & _inner_layer_mask);
0305       unsigned int ninnerhitsfromtruth = 0;
0306       unsigned int ninnerlayers = 0;
0307       for (unsigned int i = 0; i < 32; ++i) {
0308     ninnerhitsfromtruth += (0x1 & (innerhits >> i));
0309     ninnerlayers += (0x1 & (_inner_layer_mask >> i));
0310       }
0311
0312       ndiff = abs((int)ninnerhitsfromtruth-(int)ninnerlayers);
0313       if (ndiff <= 2) {
0314     _truept_particles_recoWithin2InnerHits->Fill(truept);
0315       }
0316       if (ndiff <= 1) {
0317     _truept_particles_recoWithin1InnerHit->Fill(truept);
0318       }
0319       if (ndiff == 0) {
0320     _truept_particles_recoWithExactInnerHits->Fill(truept);
0321       }
0322
0323       float diff = fabs(recopt-truept)/truept;
0324       if (diff < 0.05) {
0325     _truept_particles_recoWithin5Percent->Fill(truept);
0326       }
0327       if (diff < 0.04) {
0328     _truept_particles_recoWithin4Percent->Fill(truept);
0329     _truept_quality_particles_recoWithin4Percent->Fill(truept,track->get_quality());
0330       }
0331       if (diff < 0.03) {
0332     _truept_particles_recoWithin3Percent->Fill(truept);
0333       }
0334     }
0335   }
0336
0337   // loop over all reco particles
0338   for (SvtxTrackMap::Iter iter = trackmap->begin();
0339        iter != trackmap->end();
0340        ++iter) {
0341
0342     SvtxTrack*    track      = iter->second;
0343     float recopt = track->get_pt();
0344
0345     PHG4Particle* g4particle = trackeval->max_truth_particle_by_nclusters(track);
0346     float truept = sqrt(pow(g4particle->get_px(),2)+pow(g4particle->get_py(),2));
0347
0348     if (trutheval->get_embed(g4particle) > 0) {
0349       // embedded results (quality or preformance measures)
0350       _truept_dptoverpt->Fill(truept,(recopt-truept)/truept);
0351       _truept_dca->Fill(truept,track->get_dca2d());
0352
0353       _recopt_quality->Fill(recopt,track->get_quality());
0354     } else {
0355       // non-embedded results (purity measures)
0356       _recopt_tracks_all->Fill(recopt);
0357       _recopt_quality_tracks_all->Fill(recopt,track->get_quality());
0358
0359       unsigned int nfromtruth = trackeval->get_nclusters_contribution(track,g4particle);
0360
0361       unsigned int ndiff = abs((int)nfromtruth-(int)_nlayers);
0362       if (ndiff <= 2) {
0363     _recopt_tracks_recoWithin2Hits->Fill(recopt);
0364       }
0365       if (ndiff <= 1) {
0366     _recopt_tracks_recoWithin1Hit->Fill(recopt);
0367       }
0368       if (ndiff == 0) {
0369     _recopt_tracks_recoWithExactHits->Fill(recopt);
0370       }
0371
0372       unsigned int layersfromtruth = trackeval->get_nclusters_contribution_by_layer(track,g4particle);
0373       unsigned int innerhits = (layersfromtruth & _inner_layer_mask);
0374       unsigned int ninnerhitsfromtruth = 0;
0375       unsigned int ninnerlayers = 0;
0376       for (unsigned int i = 0; i < 32; ++i) {
0377     ninnerhitsfromtruth += (0x1 & (innerhits >> i));
0378     ninnerlayers += (0x1 & (_inner_layer_mask >> i));
0379       }
0380
0381       ndiff = abs((int)ninnerhitsfromtruth-(int)ninnerlayers);
0382       if (ndiff <= 2) {
0383     _recopt_tracks_recoWithin2InnerHits->Fill(recopt);
0384       }
0385       if (ndiff <= 1) {
0386     _recopt_tracks_recoWithin1InnerHit->Fill(recopt);
0387       }
0388       if (ndiff == 0) {
0389     _recopt_tracks_recoWithExactInnerHits->Fill(recopt);
0390       }
0391
0392       float diff = fabs(recopt-truept)/truept;
0393       if (diff < 0.05) {
0394     _recopt_tracks_recoWithin5Percent->Fill(recopt);
0395       }
0396       if (diff < 0.04) {
0397     _recopt_tracks_recoWithin4Percent->Fill(recopt);
0398     _recopt_quality_tracks_recoWithin4Percent->Fill(recopt,track->get_quality());
0399       }
0400       if (diff < 0.03) {
0401     _recopt_tracks_recoWithin3Percent->Fill(recopt);
0402       }
0403     }
0404   }
0405
0406   // get leading vertex
0407   SvtxVertex* maxvertex = NULL;
0408   unsigned int maxtracks = 0;
0409   for (SvtxVertexMap::Iter iter = vertexmap->begin();
0410        iter != vertexmap->end();
0411        ++iter) {
0412     SvtxVertex* vertex = iter->second;
0413
0414     if (vertex->size_tracks() > maxtracks) {
0415       maxvertex = vertex;
0416       maxtracks = vertex->size_tracks();
0417     }
0418   }
0419
0420   if (maxvertex) {
0421     PHG4VtxPoint* point = vertexeval->max_truth_point_by_ntracks(maxvertex);
0422
0423     if (point) {
0424       _dx_vertex->Fill(maxvertex->get_x() - point->get_x());
0425       _dy_vertex->Fill(maxvertex->get_y() - point->get_y());
0426       _dz_vertex->Fill(maxvertex->get_z() - point->get_z());
0427     }
0428   }
0429
0430   return 0;
0431 }
0432
0433 int SvtxSimPerformanceCheckReco::End(PHCompositeNode *topNode) {
0434  return 0;
0435 }