KalmanVertexFit/src/SimpleVertexTree.cc

0001 #include "RecoVertex/KalmanVertexFit/interface/SimpleVertexTree.h"
0002
0003 #include "TROOT.h"
0004 #include "TTree.h"
0005 #include "TH1F.h"
0006 #include <iostream>
0007
0008 using namespace std;
0009
0010 SimpleVertexTree::SimpleVertexTree(const char *filterName, const MagneticField *magField) : theFitterName(filterName) {
0011   vertexTree = new TTree(filterName, "Vertex fit results");
0012   //   trackTest
0013   //     = SimpleConfigurable<bool> (false, "SimpleVertexTree:trackTest").value();
0014   //   if (trackTest) {
0015   //     maxTrack = SimpleConfigurable<int> (100, "SimpleVertexTree:maximumTracksToStore").value();
0016   //   } else
0017   maxTrack = 0;
0018   result = new VertexFitterResult(maxTrack, magField);
0019
0020   vertexTree->Branch("vertex", (void *)result->vertexPresent(), "vertex/I");
0021   vertexTree->Branch("simPos", (void *)result->simVertexPos(), "X/F:Y/F:Z/F");
0022   vertexTree->Branch("recPos", (void *)result->recVertexPos(), "X/F:Y/F:Z/F");
0023   vertexTree->Branch("recErr", (void *)result->recVertexErr(), "X/F:Y/F:Z/F");
0024   vertexTree->Branch("nbrTrk", (void *)result->trackInformation(), "Sim/I:Rec/I:Shared/I");
0025   vertexTree->Branch("chiTot", (void *)result->chi2Information(), "chiTot/F");
0026   vertexTree->Branch("ndf", (void *)(result->chi2Information() + 1), "ndf/F");
0027   vertexTree->Branch("chiProb", (void *)(result->chi2Information() + 2), "chiProb/F");
0028   vertexTree->Branch("time", (void *)result->time(), "time/F");
0029
0030   parameterNames[0] = new TString("ptinv");
0031   parameterNames[1] = new TString("theta");
0032   parameterNames[2] = new TString("phi");
0033   parameterNames[3] = new TString("timp");
0034   parameterNames[4] = new TString("limp");
0035
0036   vertexTree->Branch("simTracks", (void *)result->numberSimTracks(), "simTracks/I");
0037   vertexTree->Branch("simTrack_recIndex", (void *)result->simTrack_recIndex(), "simTrack_recIndex[simTracks]/I");
0038   defineTrackBranch("sim", "Par", &VertexFitterResult::simParameters, "simTracks");
0039
0040   vertexTree->Branch("recTracks", (void *)result->numberRecTracks(), "recTracks/I");
0041   vertexTree->Branch("recTrack_simIndex", (void *)result->recTrack_simIndex(), "recTrack_simIndex[recTracks]/I");
0042   vertexTree->Branch("recTrack_weight", (void *)result->recTrackWeight(), "recTrack_weight[recTracks]/F");
0043   defineTrackBranch("rec", "Par", &VertexFitterResult::recParameters, "recTracks");
0044   defineTrackBranch("ref", "Par", &VertexFitterResult::refParameters, "recTracks");
0045   defineTrackBranch("rec", "Err", &VertexFitterResult::recErrors, "recTracks");
0046   defineTrackBranch("ref", "Err", &VertexFitterResult::refErrors, "recTracks");
0047
0048   numberOfVertices = 0;
0049 }
0050
0051 void SimpleVertexTree::defineTrackBranch(const TString &prefix,
0052                                          const TString &type,
0053                                          const float *(VertexFitterResult::*pfunc)(const int) const,
0054                                          const TString &index) {
0055   TString branchName, branchVariables;
0056   for (int i = 0; i < 5; i++) {
0057     branchName = prefix + type + '_' + *parameterNames[i];
0058     branchVariables = branchName + '[' + index + "]/F";
0059     vertexTree->Branch(branchName, (void *)(result->*pfunc)(i), branchVariables);
0060   }
0061 }
0062
0063 SimpleVertexTree::~SimpleVertexTree() {
0064   std::cout << std::endl << "End of SimpleVertexTree for " << theFitterName << std::endl;
0065   std::cout << std::endl << "Number of vertices fit: " << numberOfVertices << std::endl;
0066
0067   //
0068   // save current root directory
0069   //
0070   TDirectory *rootDir = gDirectory;
0071   //
0072   // close files
0073   //
0074   vertexTree->GetDirectory()->cd();
0075   vertexTree->Write();
0076   if (numberOfVertices > 0) {
0077     TH1F *resX = new TH1F(theFitterName + "_ResX", "Residual x coordinate: " + theFitterName, 100, -0.03, 0.03);
0078     TH1F *resY = new TH1F(theFitterName + "_ResY", "Residual y coordinate: " + theFitterName, 100, -0.03, 0.03);
0079     TH1F *resZ = new TH1F(theFitterName + "_ResZ", "Residual z coordinate: " + theFitterName, 100, -0.03, 0.03);
0080     TH1F *pullX = new TH1F(theFitterName + "_PullX", "Pull x coordinate: " + theFitterName, 100, -10., 10.);
0081     TH1F *pullY = new TH1F(theFitterName + "_PullY", "Pull y coordinate: " + theFitterName, 100, -10., 10.);
0082     TH1F *pullZ = new TH1F(theFitterName + "_PullZ", "Pull z coordinate: " + theFitterName, 100, -10., 10.);
0083     TH1F *chiNorm = new TH1F(theFitterName + "_ChiNorm", "Normalized chi-square: " + theFitterName, 100, 0., 10.);
0084     TH1F *chiProb = new TH1F(theFitterName + "_ChiProb", "Chi-square probability: " + theFitterName, 100, 0., 1.);
0085     vertexTree->Project(theFitterName + "_ResX", "(simPos.X-recPos.X)");
0086     vertexTree->Project(theFitterName + "_ResY", "(simPos.Y-recPos.Y)");
0087     vertexTree->Project(theFitterName + "_ResZ", "(simPos.Z-recPos.Z)");
0088     vertexTree->Project(theFitterName + "_PullX", "(simPos.X-recPos.X)/recErr.X");
0089     vertexTree->Project(theFitterName + "_PullY", "(simPos.Y-recPos.Y)/recErr.Y");
0090     vertexTree->Project(theFitterName + "_PullZ", "(simPos.Z-recPos.Z)/recErr.Z");
0091     vertexTree->Project(theFitterName + "_ChiNorm", "chiTot/ndf");
0092     vertexTree->Project(theFitterName + "_ChiProb", "chiProb");
0093     std::cout << "Mean of Residual distribution X: " << resX->GetMean() << std::endl;
0094     std::cout << "Mean of Residual distribution Y: " << resY->GetMean() << std::endl;
0095     std::cout << "Mean of Residual distribution Z: " << resZ->GetMean() << std::endl;
0096     std::cout << "RMS of Residual distribution X:  " << resX->GetRMS() << std::endl;
0097     std::cout << "RMS of Residual distribution Y:  " << resY->GetRMS() << std::endl;
0098     std::cout << "RMS of Residual distribution Z:  " << resZ->GetRMS() << std::endl;
0099     std::cout << "Mean of Pull distribution X: " << pullX->GetMean() << std::endl;
0100     std::cout << "Mean of Pull distribution Y: " << pullY->GetMean() << std::endl;
0101     std::cout << "Mean of Pull distribution Z: " << pullZ->GetMean() << std::endl;
0102     std::cout << "RMS of Pull distribution X:  " << pullX->GetRMS() << std::endl;
0103     std::cout << "RMS of Pull distribution Y:  " << pullY->GetRMS() << std::endl;
0104     std::cout << "RMS of Pull distribution Z:  " << pullZ->GetRMS() << std::endl;
0105     std::cout << "Average chi-square probability: " << chiProb->GetMean() << std::endl;
0106     std::cout << "Average normalized chi-square : " << chiNorm->GetMean() << std::endl;
0107     resX->Write();
0108     resY->Write();
0109     resZ->Write();
0110     pullX->Write();
0111     pullY->Write();
0112     pullZ->Write();
0113     chiNorm->Write();
0114     chiProb->Write();
0115   }
0116   delete vertexTree;
0117   //
0118   // restore directory
0119   //
0120   rootDir->cd();
0121   std::cout << std::endl;
0122 }
0123
0124 void SimpleVertexTree::fill(const TransientVertex &recv,
0125                             const TrackingVertex *simv,
0126                             reco::RecoToSimCollection *recSimColl,
0127                             const float &time) {
0128   result->fill(recv, simv, recSimColl, time);
0129   fill();
0130 }
0131
0132 void SimpleVertexTree::fill(const TransientVertex &recv, const float &time) {
0133   result->fill(recv, nullptr, nullptr, time);
0134   fill();
0135 }
0136
0137 void SimpleVertexTree::fill(const TrackingVertex *simv) {
0138   result->fill(TransientVertex(), simv, nullptr);
0139   fill();
0140 }
0141
0142 // void SimpleVertexTree::fill(const RecVertex & recVertex,
0143 //          const std::vector < RecTrack > & recTrackV,
0144 //          const TkSimVertex * simv, const float &time)
0145 // {
0146 //   result->fill(recVertex, recTrackV, simv, time);
0147 //   fill();
0148 // }
0149 //
0150 // void SimpleVertexTree::fill(const std::vector < RecTrack > & recTrackV,
0151 //          const TkSimVertex * simv, const float &time)
0152 // {
0153 //   result->fill(RecVertex(), recTrackV, simv, time);
0154 //   fill();
0155 // }
0156
0157 void SimpleVertexTree::fill() {
0158   ++numberOfVertices;
0159   static std::atomic<int> nFill{0};
0160   TDirectory *rootDir = gDirectory;
0161   //
0162   // fill entry
0163   //
0164   vertexTree->GetDirectory()->cd();
0165   vertexTree->Fill();
0166   if ((++nFill) % 1000 == 0)
0167     vertexTree->AutoSave();
0168   //
0169   // restore directory
0170   //
0171   rootDir->cd();
0172   result->reset();
0173 }