Tracking/TrackHitContributions/hit_matching.C

0001 #include <edm4hep/SimTrackerHitCollection.h>
0002 #include <edm4eic/TrackerHitCollection.h>
0003 #include <edm4eic/Measurement2DCollection.h>
0004 #include <edm4eic/TrackCollection.h>
0005 #include <edm4eic/MCRecoTrackerHitAssociationCollection.h>
0006 #include <podio/Frame.h>
0007 #include <podio/ROOTReader.h>
0008
0009 //SimHit collections
0010 std::vector<std::string> sim_coll_names{
0011                     "VertexBarrelHits","SiBarrelHits","TrackerEndcapHits",
0012                     "MPGDBarrelHits","BackwardMPGDEndcapHits","ForwardMPGDEndcapHits","OuterMPGDBarrelHits",
0013                     "TOFBarrelHits","TOFEndcapHits"
0014                     };
0015
0016 //TrackerHit RecHit (digitized) collections
0017 std::vector<std::string> rec_coll_names{
0018                     "SiBarrelVertexRecHits","SiBarrelTrackerRecHits","SiEndcapTrackerRecHits",
0019                     "MPGDBarrelRecHits","BackwardMPGDEndcapRecHits","ForwardMPGDEndcapRecHits","OuterMPGDBarrelRecHits",
0020                     "TOFBarrelRecHits","TOFEndcapRecHits"
0021                     };
0022
0023 //Collection IDs for RecHit collections.
0024 std::vector<unsigned int> rec_coll_ids;
0025
0026 //Simulation file to use
0027 //std::string input_file = "input/eicrecon_etof_test_local.root";
0028 //std::string input_file = "input/eicrecon_etof_test.root";
0029 std::string input_file = "input/eicrecon_out_2GeV.root";
0030
0031
0032 //------------------
0033 //Template to access 'sign' of radius in (x,y) plane
0034 template <typename T> int sgn(T val) {
0035     return (T(0) < val) - (val < T(0));
0036 }
0037
0038 //-----------------
0039 //Function to set Collection IDs
0040 vector<unsigned int> get_coll_ids(vector<string> coll_names){
0041
0042     vector<unsigned int> coll_ids;
0043
0044         TFile *f = new TFile(input_file.c_str());
0045         TTree *t = (TTree*) f->Get("podio_metadata");
0046
0047         TTreeReader tr(t);
0048         TTreeReaderArray<unsigned int> c_ids(tr,"events___idTable.m_collectionIDs");
0049         TTreeReaderArray<string> c_names(tr,"events___idTable.m_names");
0050
0051     tr.Next();
0052
0053     //Find the associated collection ID
0054     for(int iname = 0; iname < coll_names.size(); iname++) {
0055         for(int icol=0;icol<c_ids.GetSize();icol++){
0056             if( c_names[icol] == coll_names[iname] )
0057                 coll_ids.push_back(c_ids[icol]);
0058         }
0059     }
0060
0061     delete f;
0062
0063     return coll_ids;
0064 }
0065
0066 //-----------------
0067 // Main function
0068 void hit_matching(){
0069
0070     //Track collection to use: real-seeded or truth-seeded
0071     std::string trk_coll = "CentralCKFTracks";
0072
0073     //Print out information event-by-event
0074     bool print_evt_info = false;
0075
0076     //Using new Barrel MPDG digitization
0077     bool new_mpgd_digi = false;
0078
0079     //Plot x-y positions of last Forward TOF layer
0080     bool ftof_back_xy = false;
0081
0082     //Set collection Ids
0083     rec_coll_ids = get_coll_ids(rec_coll_names);
0084
0085     //Create map between RecHit collection IDs and names
0086         std::unordered_map<unsigned int, std::string> RecHitCollMap;
0087         for (int i = 0; i < rec_coll_ids.size(); i++) {
0088         RecHitCollMap[rec_coll_ids[i]] = rec_coll_names[i];
0089         }
0090
0091     //Print out Map information
0092     cout << endl << "Created Map between Collection IDs and Names:"<<endl;
0093     for (const auto& pair : RecHitCollMap) {
0094             cout << "Collection ID: "<< pair.first << " | Collection Name: " << pair.second << endl;
0095     }
0096
0097     //Create map between RecHit collection names and vector index
0098     std::unordered_map<std::string, int> index_map;
0099         for (size_t i = 0; i < rec_coll_names.size(); i++) {
0100             index_map[rec_coll_names[i]] = i;
0101         }
0102
0103     //Define Style
0104         gStyle->SetOptStat(0);
0105         gStyle->SetPadBorderMode(0);
0106         gStyle->SetFrameBorderMode(0);
0107         gStyle->SetFrameLineWidth(2);
0108         gStyle->SetLabelSize(0.035,"X");
0109     gStyle->SetLabelSize(0.035,"Y");
0110         //gStyle->SetLabelOffset(0.01,"X");
0111         //gStyle->SetLabelOffset(0.01,"Y");
0112         gStyle->SetTitleXSize(0.04);
0113         gStyle->SetTitleXOffset(0.9);
0114         gStyle->SetTitleYSize(0.04);
0115         gStyle->SetTitleYOffset(0.9);
0116
0117     //Define histograms
0118     TH1 *h1 = new TH1D("h1","Number of SimHits associated w/ primary particle",10,0,10);
0119     h1->SetLineWidth(3);h1->SetLineColor(kBlue);
0120     h1->GetXaxis()->SetTitle("Number of hits");h1->GetXaxis()->CenterTitle();
0121     h1->GetXaxis()->CenterLabels();
0122     h1->GetYaxis()->SetTitle("Number of tracks");h1->GetYaxis()->CenterTitle();
0123
0124     TH1 *h2 = new TH1D("h2","Fraction of associated SimHits that survive 'digitization'",25,0,1.1);
0125     h2->SetLineWidth(3);h2->SetLineColor(kBlue);
0126     h2->GetXaxis()->SetTitle("Fraction");h2->GetXaxis()->CenterTitle();
0127     h2->GetYaxis()->SetTitle("Number of tracks");h2->GetYaxis()->CenterTitle();
0128
0129     TH1 *h3 = new TH1D("h3","Fraction of track measurement hits that are associated w/ primary particle",25,0,1.1);
0130     h3->SetLineWidth(3);h3->SetLineColor(kBlue);
0131     h3->GetXaxis()->SetTitle("Fraction");h3->GetXaxis()->CenterTitle();
0132     h3->GetYaxis()->SetTitle("Number of tracks");h3->GetYaxis()->CenterTitle();
0133
0134     TH1 *h4 = new TH1D("h4","Fraction of associated digitized hits classified as measurement hits",25,0,1.1);
0135     h4->SetLineWidth(3);h4->SetLineColor(kBlue);
0136     h4->GetXaxis()->SetTitle("Fraction");h4->GetXaxis()->CenterTitle();
0137     h4->GetYaxis()->SetTitle("Number of tracks");h4->GetYaxis()->CenterTitle();
0138
0139     TH1 *h5 = new TH1D("h5","Fraction of track outlier hits that are associated w/ primary particle",25,0,1.1);
0140     h5->SetLineWidth(3);h5->SetLineColor(kBlue);
0141     h5->GetXaxis()->SetTitle("Fraction");h5->GetXaxis()->CenterTitle();
0142     h5->GetYaxis()->SetTitle("Number of tracks");h5->GetYaxis()->CenterTitle();
0143
0144     TH1 *h6 = new TH1D("h6","Fraction of associated digitized hits classified as outlier hits",25,0,1.1);
0145     h6->SetLineWidth(3);h6->SetLineColor(kBlue);
0146     h6->GetXaxis()->SetTitle("Fraction");h6->GetXaxis()->CenterTitle();
0147     h6->GetYaxis()->SetTitle("Number of tracks");h6->GetYaxis()->CenterTitle();
0148
0149     TH1 *h7 = new TH1D("h7","Fraction of associated digitized hits missing from track",25,0,1.1);
0150     h7->SetLineWidth(3);h7->SetLineColor(kBlue);
0151     h7->GetXaxis()->SetTitle("Fraction");h7->GetXaxis()->CenterTitle();
0152     h7->GetYaxis()->SetTitle("Number of tracks");h7->GetYaxis()->CenterTitle();
0153
0154     //Number of RecHits per collection
0155     //1) All associated RecHits for events w/ a reconstructed track
0156     //2) All associated measurement RecHits for events w/ a reconstructed track
0157     //3) All associated outlier RecHits for events w/ a reconstructed track
0158     int nBins = rec_coll_names.size();
0159
0160         TH1 *h8a = new TH1D("h8a","Events with a reconstructed track", nBins, 0, nBins); //Associated digitized hits
0161     h8a->SetLineWidth(3);h8a->SetLineColor(kBlue);
0162     h8a->GetYaxis()->SetTitle("Counts");h8a->GetYaxis()->CenterTitle();
0163
0164     TH1 *h8b = new TH1D("h8b","Events with a reconstructed track", nBins, 0, nBins); //Associated track measurement hits
0165         h8b->SetFillColor(kGreen+2);h8b->SetFillStyle(3004);h8b->SetLineWidth(0);
0166         h8b->GetYaxis()->SetTitle("Counts");h8b->GetYaxis()->CenterTitle();
0167
0168     TH1 *h8c = new TH1D("h8c","Events with a reconstructed track", nBins, 0, nBins); //Associated outlier hits
0169         h8c->SetFillColor(kRed);h8c->SetFillStyle(3005);h8c->SetLineWidth(0);
0170         h8c->GetYaxis()->SetTitle("Counts");h8c->GetYaxis()->CenterTitle();
0171
0172         //Assign bin labels
0173         for (size_t i = 0; i < rec_coll_names.size(); i++) {
0174             h8a->GetXaxis()->SetBinLabel(i + 1, rec_coll_names[i].c_str());
0175         h8b->GetXaxis()->SetBinLabel(i + 1, rec_coll_names[i].c_str());
0176         h8c->GetXaxis()->SetBinLabel(i + 1, rec_coll_names[i].c_str());
0177         }
0178
0179     //Number of RecHits per collection
0180         //1) All associated RecHits
0181         //2) All associated Measurement2D Hits
0182         TH1 *h9a = new TH1D("h9a","", nBins, 0, nBins); //Associated RecHits
0183         h9a->SetLineWidth(3);h9a->SetLineColor(kBlue);
0184         h9a->GetYaxis()->SetTitle("Counts");h9a->GetYaxis()->CenterTitle();
0185
0186         TH1 *h9b = new TH1D("h9b","", nBins, 0, nBins); //Associated Measurement2D
0187         h9b->SetFillColor(kMagenta+2);h9b->SetFillStyle(3004);h9b->SetLineWidth(0);
0188         h9b->GetYaxis()->SetTitle("Counts");h9b->GetYaxis()->CenterTitle();
0189
0190     //Assign bin labels
0191     for (size_t i = 0; i < rec_coll_names.size(); i++) {
0192                 h9a->GetXaxis()->SetBinLabel(i + 1, rec_coll_names[i].c_str());
0193         h9b->GetXaxis()->SetBinLabel(i + 1, rec_coll_names[i].c_str());
0194     }
0195
0196     //R vs. Z: Events with a reconstructed track
0197     TH2 *hh1a = new TH2D("hh1a","All associated digitized hits",1000,-1800,2000,1000,-800,800);
0198     hh1a->GetXaxis()->SetTitle("z [mm]");hh1a->GetXaxis()->CenterTitle();
0199     hh1a->GetYaxis()->SetTitle("signed r [mm]");hh1a->GetYaxis()->CenterTitle();
0200
0201     TH2 *hh1b = new TH2D("hh1b","Associated track measurement hits",1000,-1800,2000,1000,-800,800);
0202     hh1b->GetXaxis()->SetTitle("z [mm]");hh1b->GetXaxis()->CenterTitle();
0203     hh1b->GetYaxis()->SetTitle("signed r [mm]");hh1b->GetYaxis()->CenterTitle();
0204
0205     TH2 *hh1c = new TH2D("hh1c","Associated track outlier hits",1000,-1800,2000,1000,-800,800);
0206     hh1c->GetXaxis()->SetTitle("z [mm]");hh1c->GetXaxis()->CenterTitle();
0207     hh1c->GetYaxis()->SetTitle("signed r [mm]");hh1c->GetYaxis()->CenterTitle();
0208
0209     TH2 *hh1d = new TH2D("hh1d","Associated digitized hits missing from track",1000,-1800,2000,1000,-800,800);
0210     hh1d->GetXaxis()->SetTitle("z [mm]");hh1d->GetXaxis()->CenterTitle();
0211     hh1d->GetYaxis()->SetTitle("signed r [mm]");hh1d->GetYaxis()->CenterTitle();
0212
0213     //Forward TOF (last layer) y vs. x: Events with a reconstructed track
0214     TH2 *hh2a = new TH2D("hh2a","All associated digitized hits",1000,-600,600,1000,-600,600);
0215         hh2a->GetXaxis()->SetTitle("x [mm]");hh2a->GetXaxis()->CenterTitle();
0216         hh2a->GetYaxis()->SetTitle("y [mm]");hh2a->GetYaxis()->CenterTitle();
0217
0218         TH2 *hh2b = new TH2D("hh2b","Associated track measurement hits",1000,-600,600,1000,-600,600);
0219         hh2b->GetXaxis()->SetTitle("x [mm]");hh2b->GetXaxis()->CenterTitle();
0220         hh2b->GetYaxis()->SetTitle("y [mm]");hh2b->GetYaxis()->CenterTitle();
0221
0222         TH2 *hh2c = new TH2D("hh2c","Associated track outlier hits",1000,-600,600,1000,-600,600);
0223         hh2c->GetXaxis()->SetTitle("x [mm]");hh2c->GetXaxis()->CenterTitle();
0224         hh2c->GetYaxis()->SetTitle("y [mm]");hh2c->GetYaxis()->CenterTitle();
0225
0226         TH2 *hh2d = new TH2D("hh2d","Associated digitized hits missing from track",1000,-600,600,1000,-600,600);
0227         hh2d->GetXaxis()->SetTitle("x [mm]");hh2d->GetXaxis()->CenterTitle();
0228         hh2d->GetYaxis()->SetTitle("y [mm]");hh2d->GetYaxis()->CenterTitle();
0229
0230
0231     //Open file
0232     podio::ROOTReader r;
0233     r.openFile(input_file);
0234
0235     auto nevents = r.getEntries(podio::Category::Event);
0236     cout<<"---------------"<<endl;
0237     cout<<"Total number of events = "<<nevents<<"!"<<endl;
0238     cout<<"---------------"<<endl;
0239
0240     //Loop over events. Can I do this more easily?
0241     for( unsigned int ievent = 0; ievent < nevents; ievent++){
0242
0243         if((ievent%1000)==0) cout<<"Processed event "<<ievent<<endl;
0244
0245         //Get event
0246         auto f = podio::Frame(r.readNextEntry(podio::Category::Event));
0247         if(print_evt_info)
0248             cout<<"For event "<<ievent<<":"<<endl;
0249
0250         //RawHit / SimHit association collection
0251         auto& raw_hit_assocs = f.get<edm4eic::MCRecoTrackerHitAssociationCollection>("CentralTrackingRawHitAssociations");
0252
0253         //Total number of SimHits
0254         int tot_simhits(0);
0255         //Number of SimHits that are associated with primary particle
0256         int matched_simhits(0);
0257
0258         //Loop over Simhit collections
0259         if(print_evt_info)
0260             cout<<"Number of SimHit collections = "<<sim_coll_names.size()<<endl;
0261         for(int icoll = 0; icoll < sim_coll_names.size(); icoll++){
0262             //Simhits
0263             auto coll_name = sim_coll_names.at(icoll);
0264             auto& simhit_coll = f.get<edm4hep::SimTrackerHitCollection>(coll_name);
0265             auto num_simhits = simhit_coll.size();
0266
0267             if(print_evt_info){
0268                 printf("\nNumber of %s = %lu \n",coll_name.c_str(),num_simhits);
0269                 cout<<"Hit Collection ID | Hit Index | Hit CellID | Hit Quality | MC Index | MC PDG ID | MC Status"<<endl;
0270             }
0271             for(int ihit = 0; ihit<num_simhits; ihit++){
0272                 auto hit = simhit_coll.at(ihit);
0273                 auto hit_collid = hit.getObjectID().collectionID;
0274                 auto hit_index = hit.getObjectID().index;
0275                 auto hit_cellid = hit.getCellID();
0276                 auto quality = hit.getQuality();
0277
0278                 //MCParticle associated with simhit
0279                 auto hit_mcpart = hit.getMCParticle();
0280                 auto mc_index = hit_mcpart.getObjectID().index;
0281                 auto mc_pdg = hit_mcpart.getPDG();
0282                 auto mc_status = hit_mcpart.getGeneratorStatus();
0283
0284                 if(print_evt_info)
0285                     printf("%u | %d | %lu | %d | %d| %d | %d \n",hit_collid, hit_index, hit_cellid, quality, mc_index, mc_pdg, mc_status);
0286
0287                 if(mc_status==1 && quality==0) matched_simhits++;
0288                 tot_simhits++;
0289
0290                 //For new Barrel MPGDs digitization, a single SimHit usually is converted into 2 RecHits
0291                 //So, we add an additional count for those to keep the numbers consistent
0292                 if( new_mpgd_digi && (coll_name=="MPGDBarrelHits" || coll_name=="OuterMPGDBarrelHits") ){
0293                     if(mc_status==1 && quality==0) matched_simhits++;
0294                                     tot_simhits++;
0295                 }
0296             } //Loop over SimHits in a given collection
0297         } //Loop over SimHit collections
0298
0299         //-------------
0300         //Total number of RecHits associated to primary particle
0301         int tot_rechits(0);
0302
0303         //Loop over RecHit collections
0304         if(print_evt_info){
0305             cout<<endl<<endl;
0306             cout<<"Number of RecHit collections = "<<rec_coll_names.size()<<endl;
0307         }
0308         for(int icoll = 0; icoll < rec_coll_names.size(); icoll++){
0309             //RecHits
0310             auto coll_name = rec_coll_names.at(icoll);
0311             auto& rechit_coll = f.get<edm4eic::TrackerHitCollection>(coll_name);
0312             auto num_rechits = rechit_coll.size();
0313
0314             if(print_evt_info){
0315                 printf("\nNumber of %s = %lu \n",coll_name.c_str(),num_rechits);
0316                 cout<<"Hit Collection ID | Hit Index | Hit CellID"<<endl;
0317             }
0318             for(int ihit = 0; ihit<num_rechits; ihit++){
0319
0320                 auto hit = rechit_coll.at(ihit);
0321                 auto hit_collid = hit.getObjectID().collectionID;
0322                 auto hit_index = hit.getObjectID().index;
0323                 auto hit_cellid = hit.getCellID();
0324
0325                 //RecHit positions
0326                 auto hit_x = hit.getPosition().x;
0327                 auto hit_y = hit.getPosition().y;
0328                 auto hit_z = hit.getPosition().z;
0329                 double hit_radius = sgn<double>(hit_x) * sqrt(hit_x*hit_x + hit_y*hit_y);
0330
0331                 if(print_evt_info)
0332                     printf("%u | %d | %lu \n",hit_collid, hit_index, hit_cellid);
0333
0334                 //Find corresponding SimHit and check if it is associated w/ primary particle
0335                 auto raw_hit = hit.getRawHit();
0336
0337                 for (const auto raw_hit_assoc : raw_hit_assocs) {
0338                     if (raw_hit_assoc.getRawHit() == raw_hit) {
0339                         auto sim_hit = raw_hit_assoc.getSimHit();
0340                                     auto mc_particle = sim_hit.getMCParticle();
0341                         auto mc_status = sim_hit.getMCParticle().getGeneratorStatus();
0342                         auto quality = sim_hit.getQuality();
0343
0344                         if(mc_status==1 && quality==0){
0345                             tot_rechits++;
0346
0347                             //If there is at least 1 track, fill histogram of collection RecHits
0348                             auto& track_coll = f.get<edm4eic::TrackCollection>(trk_coll);
0349                                     auto num_tracks = track_coll.size();
0350
0351                             if(num_tracks > 0 && (index_map.find(coll_name) != index_map.end()) ){
0352                                 h8a->Fill( index_map[coll_name] );
0353                                 hh1a->Fill(hit_z,hit_radius);
0354
0355                                 if(coll_name=="TOFEndcapRecHits" && hit_z>1860.)
0356                                     hh2a->Fill(hit_x,hit_y);
0357                             }
0358
0359                             //Without track requirement
0360                             if( index_map.find(coll_name) != index_map.end() ){
0361                                                                 h9a->Fill( index_map[coll_name] );
0362                                                         }
0363
0364                         } //If hit associated to primary particle
0365                     }
0366                 } //Loop over RawHit/SimHit associations
0367             } //Loop over RecHits in a given collection
0368         } //Loop over RecHit collections
0369
0370         //-------------
0371         //Loop over Measurement2D hits (independent of track reconstruction)
0372         auto& meas2d_coll = f.get<edm4eic::Measurement2DCollection>("CentralTrackerMeasurements");
0373         for(int imeas = 0; imeas<meas2d_coll.size(); imeas++){
0374
0375             auto rec_hit = meas2d_coll.at(imeas).getHits().at(0); //One RecHit exists for a given Measurement2D
0376
0377             auto cellid = rec_hit.getCellID();
0378                         auto collid = rec_hit.getObjectID().collectionID;
0379
0380             //Get collection name for RecHit
0381                         std::string coll_name;
0382                         if (RecHitCollMap.find(collid) != RecHitCollMap.end())
0383                             coll_name = RecHitCollMap.at(collid);
0384                         else
0385                             coll_name = "UnknownCollection";
0386
0387             //Find corresponding SimHit and check if it is associated w/ primary particle
0388                         auto raw_hit = rec_hit.getRawHit();
0389
0390                         for (const auto raw_hit_assoc : raw_hit_assocs) {
0391                             if (raw_hit_assoc.getRawHit() == raw_hit) {
0392                                     auto sim_hit = raw_hit_assoc.getSimHit();
0393                                         auto mc_particle = sim_hit.getMCParticle();
0394                                         auto mc_status = sim_hit.getMCParticle().getGeneratorStatus();
0395                                         auto quality = sim_hit.getQuality();
0396
0397                                         if(mc_status==1 && quality==0){
0398                                             if( index_map.find(coll_name) != index_map.end() ){
0399                                                     h9b->Fill( index_map[coll_name] );
0400                                                 }
0401                                         } //If hit associated to primary particle
0402                             }
0403             } //Loop over RawHit/SimHit associations
0404         } //Loop over Measurement2D points
0405
0406         //-------------
0407         //RecHits that are measurement/outlier hits of a given track
0408         auto& track_coll = f.get<edm4eic::TrackCollection>(trk_coll);
0409         auto num_tracks = track_coll.size();
0410
0411         //Will want to adjust this for events with multiple tracks (not single particle)
0412         int num_meas(0); //number of measurements
0413         int num_out(0); //number of outliers
0414         int num_meas_corr(0); //number of correctly-identified measurements
0415         int num_out_wrong(0); //number of wrongly-identified outliers
0416
0417         //Print track and associated measurement RecHit info
0418         if(print_evt_info)
0419             cout<<endl<<endl<<"Number of tracks = "<<num_tracks<<endl;
0420
0421         for(int itrack = 0; itrack<num_tracks; itrack++){
0422
0423             auto track = track_coll.at(itrack); //Track
0424
0425             //-------------
0426             //Get measurements
0427             if(print_evt_info)
0428                 cout<<"Track Index | Measurement number | Measurement CellID | Collection | Should be measurement?"<<endl;
0429             auto meas2d = track.getMeasurements(); //Associated Measurement2D array
0430             for( int imeas = 0; imeas < meas2d.size(); imeas++){
0431
0432                 auto rec_hit = meas2d.at(imeas).getHits().at(0); //One RecHit exists for a given Measurement2D
0433
0434                 auto cellid = rec_hit.getCellID();
0435                 auto collid = rec_hit.getObjectID().collectionID;
0436
0437                 //RecHit positions
0438                             auto hit_x = rec_hit.getPosition().x;
0439                             auto hit_y = rec_hit.getPosition().y;
0440                             auto hit_z = rec_hit.getPosition().z;
0441                             double hit_radius = sgn<double>(hit_x) * sqrt(hit_x*hit_x + hit_y*hit_y);
0442
0443                 //Get collection name for RecHit
0444                 std::string coll_name;
0445                 if (RecHitCollMap.find(collid) != RecHitCollMap.end())
0446                     coll_name = RecHitCollMap.at(collid);
0447                 else
0448                     coll_name = "UnknownCollection";
0449
0450                 //Find associated SimHit and see if hit should be used in track fit
0451                 std::string should_be_meas = "N";
0452                 auto raw_hit = rec_hit.getRawHit();
0453
0454                 for (const auto raw_hit_assoc : raw_hit_assocs) {
0455                     if (raw_hit_assoc.getRawHit() == raw_hit) {
0456                         auto sim_hit = raw_hit_assoc.getSimHit();
0457                                     auto mc_particle = sim_hit.getMCParticle();
0458                         auto mc_status = sim_hit.getMCParticle().getGeneratorStatus();
0459                         auto quality = sim_hit.getQuality();
0460
0461                         if(mc_status==1 && quality==0){
0462                             should_be_meas = "Y";
0463                             num_meas_corr++;
0464
0465                             //Fill associated measurement hit collection histogram
0466                             if( (index_map.find(coll_name) != index_map.end()) ){
0467                                                                 h8b->Fill( index_map[coll_name] );
0468                                 hh1b->Fill(hit_z,hit_radius);
0469
0470                                 if(coll_name=="TOFEndcapRecHits" && hit_z>1860.)
0471                                                                         hh2b->Fill(hit_x,hit_y);
0472                                                         }
0473
0474                         } //If hit associated to primary particle
0475                     }
0476                 } //Loop over RawHit/SimHit associations
0477
0478                 if(print_evt_info)
0479                     printf("%d | %d | %lu | %s | %s \n",itrack, imeas, cellid, coll_name.c_str(), should_be_meas.c_str());
0480                 num_meas++;
0481             }
0482
0483             //-------------
0484             //Get outliers (currently stored in Trajectory)
0485             if(print_evt_info)
0486                 cout<<"Track Index | Outlier number | Measurement CellID | Collection | Should be measurement?"<<endl;
0487             auto traj = track.getTrajectory();
0488             auto out2d = traj.getOutliers_deprecated();
0489
0490             for( int iout = 0; iout < out2d.size(); iout++){
0491
0492                 auto out_hit = out2d.at(iout).getHits().at(0); //One RecHit exists for a given Measurement2D
0493
0494                 auto cellid = out_hit.getCellID();
0495                 auto collid = out_hit.getObjectID().collectionID;
0496
0497                 //RecHit positions
0498                                 auto hit_x = out_hit.getPosition().x;
0499                                 auto hit_y = out_hit.getPosition().y;
0500                                 auto hit_z = out_hit.getPosition().z;
0501                                 double hit_radius = sgn<double>(hit_x) * sqrt(hit_x*hit_x + hit_y*hit_y);
0502
0503                                 //Get collection name for RecHit
0504                                 std::string coll_name;
0505                                 if (RecHitCollMap.find(collid) != RecHitCollMap.end())
0506                                         coll_name = RecHitCollMap.at(collid);
0507                                 else
0508                                         coll_name = "UnknownCollection";
0509
0510                 //Find associated SimHit and see if hit should be used in track fit
0511                 std::string should_be_meas = "N";
0512                 auto raw_hit = out_hit.getRawHit();
0513
0514                 for (const auto raw_hit_assoc : raw_hit_assocs) {
0515                     if (raw_hit_assoc.getRawHit() == raw_hit) {
0516                         auto sim_hit = raw_hit_assoc.getSimHit();
0517                                     auto mc_particle = sim_hit.getMCParticle();
0518                         auto mc_status = sim_hit.getMCParticle().getGeneratorStatus();
0519                         auto quality = sim_hit.getQuality();
0520
0521                         if(mc_status==1 && quality==0){
0522                             should_be_meas = "Y";
0523                             num_out_wrong++;
0524
0525                             //Fill associated outlier hit collection histogram
0526                                                         if( (index_map.find(coll_name) != index_map.end()) ){
0527                                                                 h8c->Fill( index_map[coll_name] );
0528                                 hh1c->Fill(hit_z,hit_radius);
0529
0530                                 if(coll_name=="TOFEndcapRecHits" && hit_z>1860.)
0531                                                                         hh2c->Fill(hit_x,hit_y);
0532                                                         }
0533                         }
0534                     }
0535                 } //Loop over RawHit/SimHit associations
0536
0537                 if(print_evt_info)
0538                     printf("%d | %d | %lu | %s | %s \n",itrack, iout, cellid, coll_name.c_str(), should_be_meas.c_str());
0539                 num_out++;
0540             }
0541
0542         } //End loop over tracks
0543
0544         //Print some event statistics
0545         if(print_evt_info){
0546             cout<<endl<<"Event statistics:"<<endl;
0547             cout<<"SimHits (Barrel MPGD corrected) associated with primary particle = "<<matched_simhits<<" / "<<tot_simhits<<endl;
0548             cout<<"RecHits associated with primary particle = "<<tot_rechits<<endl;
0549             cout<<"Number of correctly-identified measurement hits = "<<num_meas_corr<<" / "<<num_meas<<endl;
0550             cout<<"Number of outlier hits that should be measurement hits = "<<num_out_wrong<<" / "<<num_out<<endl;
0551         }
0552
0553         //Total number of hits (meas. + out.) in track that are associated to a primary-particle hit
0554         auto track_assoc_hit_tot = num_meas_corr+num_out_wrong;
0555         //Total number of digitized hits (RecHits) associated with primary particle completely missing from track
0556         auto simhits_missing = matched_simhits - track_assoc_hit_tot;
0557         auto rechits_missing = simhits_missing - ( matched_simhits - tot_rechits ); //Don't count any hits lost during digitization
0558
0559         if(print_evt_info)
0560             cout<<"Number of RecHits associated with primary particle completely missing from track = "<<rechits_missing<<" / "<<tot_rechits<<endl;
0561
0562         //Add line between events
0563         if(print_evt_info)
0564             cout<<endl<<"---------------"<<endl;
0565
0566         //Fill histograms
0567         if(num_tracks>0){ //Require at least 1 track
0568             h1->Fill(matched_simhits);
0569             if(matched_simhits>0)
0570                 h2->Fill( ((float) tot_rechits) / matched_simhits );
0571             if(num_meas>0)
0572                 h3->Fill( ((float) num_meas_corr) / num_meas );
0573             if(tot_rechits>0)
0574                 h4->Fill( ((float) num_meas_corr) / tot_rechits );
0575             if(num_out>0)
0576                 h5->Fill( ((float) num_out_wrong) / num_out );
0577             if(tot_rechits>0)
0578                 h6->Fill( ((float) num_out_wrong) / tot_rechits );
0579             if(tot_rechits>0)
0580                 h7->Fill( ((float) rechits_missing) / tot_rechits );
0581         }
0582
0583     } //End loop over events
0584
0585     //Scale and Subtract histograms
0586     TH1* h8bs = (TH1D*)h8b->Clone("h8bs");
0587     h8bs->GetYaxis()->SetTitle("Fraction of associated digitized hits");
0588     h8bs->Divide(h8a);
0589
0590     TH1* h8cs = (TH1D*)h8c->Clone("h8cs");
0591     h8cs->GetYaxis()->SetTitle("Fraction of associated digitized hits");
0592         h8cs->Divide(h8a);
0593
0594     hh1d->Add(hh1a,hh1b,1,-1); //Subtract track measurement hits from associated hits
0595     hh1d->Add(hh1c,-1); //Also subtract outlier hits
0596
0597     hh2d->Add(hh2a,hh2b,1,-1); //Subtract track measurement hits from associated hits
0598         hh2d->Add(hh2c,-1); //Also subtract outlier hits
0599
0600     //Histogram stacking
0601     THStack *stack1 = new THStack("stack1", "Stacked Histograms");
0602         stack1->Add(h8b);
0603         stack1->Add(h8c);
0604
0605     THStack *stack2 = new THStack("stack2", "Events with a reconstructed track; ; Fraction of associated digitized hits");
0606         stack2->Add(h8bs);
0607         stack2->Add(h8cs);
0608
0609     //Make plots
0610     TCanvas *c1 = new TCanvas("c1");
0611     h1->Draw();
0612
0613     TLatex *tex1 = new TLatex();
0614     tex1->SetTextSize(0.03);
0615     //tex1->DrawLatexNDC(0.2,0.7,"Events with a reconstructed truth-seeded track"); //Uncomment if using truth-seeded tracks
0616     tex1->DrawLatexNDC(0.2,0.7,"Events with a reconstructed real-seeded track");
0617
0618     TCanvas *c2 = new TCanvas("c2");
0619     h2->Draw();
0620
0621     TCanvas *c3 = new TCanvas("c3");
0622     h3->Draw();
0623
0624     TCanvas *c4 = new TCanvas("c4");
0625     h4->Draw();
0626
0627     TCanvas *c5 = new TCanvas("c5");
0628     h5->Draw();
0629
0630     TCanvas *c6 = new TCanvas("c6");
0631     h6->Draw();
0632
0633     TCanvas *c7 = new TCanvas("c7");
0634     h7->Draw();
0635
0636     TCanvas *c8 = new TCanvas("c8");
0637         h8a->Draw();
0638         gPad->Update(); // Update the pad to set axis range correctly
0639     stack1->Draw("HIST same");
0640
0641     TLegend *legend = new TLegend(0.4, 0.75, 0.9, 0.9);
0642     legend->AddEntry(h8a,"Associated Digitized Hits (RecHits)","l");
0643     legend->AddEntry(h8b,"Associated Tracker Measurement Hits", "f");
0644     legend->AddEntry(h8c,"Associated Tracker Outlier Hits", "f");
0645     legend->SetBorderSize(0);
0646         legend->SetFillStyle(0);
0647         legend->SetTextSize(0.035);
0648     legend->Draw();
0649
0650     TCanvas *c8s = new TCanvas("c8s");
0651     c8s->SetGrid();
0652     stack2->Draw("HIST");
0653
0654     TCanvas *c9 = new TCanvas("c9");
0655     h9a->Draw();
0656     h9b->Draw("same");
0657
0658     TLegend *legend2 = new TLegend(0.4, 0.75, 0.9, 0.9);
0659         legend2->AddEntry(h9a,"Associated Digitized Hits (RecHits)","l");
0660         legend2->AddEntry(h9b,"Associated Measurement2D Hits", "f");
0661         legend2->SetBorderSize(0);
0662         legend2->SetFillStyle(0);
0663         legend2->SetTextSize(0.035);
0664         legend2->Draw();
0665
0666     TCanvas *c10a = new TCanvas("c10a");
0667     hh1a->Draw();
0668
0669     TCanvas *c10b = new TCanvas("c10b");
0670         hh1b->Draw();
0671
0672     TCanvas *c10c = new TCanvas("c10c");
0673         hh1c->Draw();
0674
0675     TCanvas *c10d = new TCanvas("c10d");
0676         hh1d->Draw();
0677
0678     TCanvas *c11a = new TCanvas("c11a");
0679         hh2a->Draw();
0680
0681     TCanvas *c11b = new TCanvas("c11b");
0682         hh2b->Draw();
0683
0684         TCanvas *c11c = new TCanvas("c11c");
0685         hh2c->Draw();
0686
0687         TCanvas *c11d = new TCanvas("c11d");
0688         hh2d->Draw();
0689
0690     //Print plots to file
0691     cout<<endl;
0692     c1->Print("plots/hit_matching.pdf[");
0693     c1->Print("plots/hit_matching.pdf");
0694     c2->Print("plots/hit_matching.pdf");
0695     c3->Print("plots/hit_matching.pdf");
0696     c4->Print("plots/hit_matching.pdf");
0697     c5->Print("plots/hit_matching.pdf");
0698     c6->Print("plots/hit_matching.pdf");
0699     c7->Print("plots/hit_matching.pdf");
0700     c8->Print("plots/hit_matching.pdf");
0701     c8s->Print("plots/hit_matching.pdf");
0702     c9->Print("plots/hit_matching.pdf");
0703     c10a->Print("plots/hit_matching.pdf");
0704     c10b->Print("plots/hit_matching.pdf");
0705     c10c->Print("plots/hit_matching.pdf");
0706     c10d->Print("plots/hit_matching.pdf");
0707     c11a->Print("plots/hit_matching.pdf");
0708     c11b->Print("plots/hit_matching.pdf");
0709     c11c->Print("plots/hit_matching.pdf");
0710     c11d->Print("plots/hit_matching.pdf");
0711     c11d->Print("plots/hit_matching.pdf]");
0712 }