Tracking/FastTrackingEval/FastTrackingEval.cc

0001 /*!
0002  *  \file       FastTrackingEval.cc
0003  *  \brief      Evaluation module for PHG4TrackFastSim output
0004  *  \details    input: PHG4TruthInfoContainer, SvtxTrackMap with SvtxTrack_FastSim inside
0005  *  \author     Haiwang Yu <yuhw@nmsu.edu>
0006  */
0007
0008 #include "FastTrackingEval.h"
0009
0010 #include <phool/phool.h>
0011 #include <phool/getClass.h>
0012 #include <fun4all/Fun4AllReturnCodes.h>
0013 #include <g4main/PHG4HitContainer.h>
0014 #include <g4main/PHG4TruthInfoContainer.h>
0015 #include <g4main/PHG4Particle.h>
0016 #include <g4main/PHG4Hit.h>
0017 #include <g4main/PHG4VtxPoint.h>
0018 #include <fun4all/PHTFileServer.h>
0019 #include <fun4all/Fun4AllServer.h>
0020
0021 /*
0022 #include <g4hough/SvtxVertexMap.h>
0023 #include <g4hough/SvtxVertex.h>
0024 #include <g4hough/SvtxTrackMap.h>
0025 #include <g4hough/SvtxTrack.h>
0026 #include <g4hough/SvtxTrack_FastSim.h>
0027 #include <g4hough/SvtxClusterMap.h>
0028 #include <g4hough/SvtxCluster.h>
0029 #include <g4hough/SvtxHitMap.h>
0030 #include <g4hough/SvtxHit.h>
0031 */
0032
0033 #include <trackbase_historic/SvtxVertexMap.h>
0034 #include <trackbase_historic/SvtxVertex.h>
0035 #include <trackbase_historic/SvtxTrackMap.h>
0036 #include <trackbase_historic/SvtxTrack.h>
0037 #include <trackbase_historic/SvtxTrack_FastSim.h>
0038 //#include <trackbase_historic/SvtxClusterMap.h>
0039 //#include <trackbase_historic/SvtxCluster.h>
0040 //#include <trackbase_historic/SvtxHitMap.h>
0041 //#include <trackbase_historic/SvtxHit.h>
0042
0043 #include <g4eval/SvtxEvalStack.h>
0044 #include <g4eval/SvtxTrackEval.h>
0045 #include <g4eval/SvtxClusterEval.h>
0046 #include <g4eval/SvtxTruthEval.h>
0047 #include <g4eval/SvtxVertexEval.h>
0048 #include <g4eval/SvtxHitEval.h>
0049
0050 #include <TTree.h>
0051 #include <TH2D.h>
0052 #include <TVector3.h>
0053
0054 #include <math.h>
0055
0056 #include <iostream>
0057
0058 #define LogDebug(exp)       std::cout<<"DEBUG: "    <<__FILE__<<": "<<__LINE__<<": "<< exp <<"\n"
0059 #define LogError(exp)       std::cout<<"ERROR: "    <<__FILE__<<": "<<__LINE__<<": "<< exp <<"\n"
0060 #define LogWarning(exp) std::cout<<"WARNING: "  <<__FILE__<<": "<<__LINE__<<": "<< exp <<"\n"
0061
0062 using namespace std;
0063
0064 //----------------------------------------------------------------------------//
0065 //-- Constructor:
0066 //--  simple initialization
0067 //----------------------------------------------------------------------------//
0068 FastTrackingEval::FastTrackingEval(const string &name, const string &filename, const string &trackmapname) :
0069   SubsysReco(name),
0070   _outfile_name(filename),
0071   _trackmapname(trackmapname),
0072   _event(0),
0073   _flags(NONE),
0074   _eval_tree_tracks( NULL)
0075 {
0076 }
0077
0078 //----------------------------------------------------------------------------//
0079 //-- Init():
0080 //--   Intialize all histograms, trees, and ntuples
0081 //----------------------------------------------------------------------------//
0082 int FastTrackingEval::Init(PHCompositeNode *topNode) {
0083     cout << PHWHERE << " Openning file " << _outfile_name << endl;
0084     PHTFileServer::get().open(_outfile_name, "RECREATE");
0085
0086     // create TTree
0087     _eval_tree_tracks = new TTree("tracks", "FastSim Eval => tracks");
0088     _eval_tree_tracks->Branch("event", &event, "event/I");
0089     _eval_tree_tracks->Branch("gtrackID", &gtrackID, "gtrackID/I");
0090     _eval_tree_tracks->Branch("gflavor", &gflavor, "gflavor/I");
0091     _eval_tree_tracks->Branch("gpx", &gpx, "gpx/F");
0092     _eval_tree_tracks->Branch("gpy", &gpy, "gpy/F");
0093     _eval_tree_tracks->Branch("gpz", &gpz, "gpz/F");
0094     _eval_tree_tracks->Branch("gpt", &gpt, "gpt/F");
0095     _eval_tree_tracks->Branch("gp", &gp, "gp/F");
0096     _eval_tree_tracks->Branch("gtheta", &gtheta, "gtheta/F");
0097     _eval_tree_tracks->Branch("geta", &geta, "geta/F");
0098     _eval_tree_tracks->Branch("gphi", &gphi, "gphi/F");
0099     _eval_tree_tracks->Branch("gvx", &gvx, "gvx/F");
0100     _eval_tree_tracks->Branch("gvy", &gvy, "gvy/F");
0101     _eval_tree_tracks->Branch("gvz", &gvz, "gvz/F");
0102     _eval_tree_tracks->Branch("trackID", &trackID, "trackID/I");
0103     _eval_tree_tracks->Branch("charge", &charge, "charge/I");
0104     _eval_tree_tracks->Branch("nhits", &nhits, "nhits/I");
0105     _eval_tree_tracks->Branch("px", &px, "px/F");
0106     _eval_tree_tracks->Branch("py", &py, "py/F");
0107     _eval_tree_tracks->Branch("pz", &pz, "pz/F");
0108     _eval_tree_tracks->Branch("pt", &pt, "pt/F");
0109     _eval_tree_tracks->Branch("p", &p, "p/F");
0110     _eval_tree_tracks->Branch("theta", &theta, "theta/F");
0111     _eval_tree_tracks->Branch("eta", &eta, "eta/F");
0112     _eval_tree_tracks->Branch("phi", &phi, "phi/F");
0113     _eval_tree_tracks->Branch("dca2d", &dca2d, "dca2d/F");
0114
0115     _h2d_Delta_mom_vs_truth_eta = new TH2D("_h2d_Delta_mom_vs_truth_eta",
0116             "#frac{#Delta p}{truth p} vs. truth #eta", 54, -4.5, +4.5, 1000, -1,
0117             1);
0118
0119     _h2d_Delta_mom_vs_truth_mom = new TH2D("_h2d_Delta_mom_vs_truth_mom",
0120             "#frac{#Delta p}{truth p} vs. truth p", 41, -0.5, 40.5, 1000, -1,
0121             1);
0122
0123     return Fun4AllReturnCodes::EVENT_OK;
0124 }
0125
0126 //----------------------------------------------------------------------------//
0127 //-- process_event():
0128 //--   Call user instructions for every event.
0129 //--   This function contains the analysis structure.
0130 //----------------------------------------------------------------------------//
0131 int FastTrackingEval::process_event(PHCompositeNode *topNode) {
0132     _event++;
0133     if (Verbosity() >= 2 and _event % 1000 == 0)
0134         cout << PHWHERE << "Events processed: " << _event << endl;
0135
0136     //std::cout << "Opening nodes" << std::endl;
0137     GetNodes(topNode);
0138
0139     //std::cout << "Filling trees" << std::endl;
0140     fill_tree(topNode);
0141     //std::cout << "DONE" << std::endl;
0142
0143     return Fun4AllReturnCodes::EVENT_OK;
0144 }
0145
0146 //----------------------------------------------------------------------------//
0147 //-- End():
0148 //--   End method, wrap everything up
0149 //----------------------------------------------------------------------------//
0150 int FastTrackingEval::End(PHCompositeNode *topNode) {
0151
0152     PHTFileServer::get().cd(_outfile_name);
0153
0154     _eval_tree_tracks->Write();
0155
0156     _h2d_Delta_mom_vs_truth_eta->Write();
0157     _h2d_Delta_mom_vs_truth_mom->Write();
0158
0159     //PHTFileServer::get().close();
0160
0161     return Fun4AllReturnCodes::EVENT_OK;
0162 }
0163
0164 //----------------------------------------------------------------------------//
0165 //-- fill_tree():
0166 //--   Fill the trees with truth, track fit, and cluster information
0167 //----------------------------------------------------------------------------//
0168 void FastTrackingEval::fill_tree(PHCompositeNode *topNode) {
0169
0170     // Make sure to reset all the TTree variables before trying to set them.
0171     reset_variables();
0172     //std::cout << "A1" << std::endl;
0173     event = _event;
0174
0175     if (!_truth_container) {
0176         LogError("_truth_container not found!");
0177         return;
0178     }
0179
0180     if (!_trackmap) {
0181         LogError("_trackmap not found!");
0182         return;
0183     }
0184
0185     PHG4TruthInfoContainer::ConstRange range =
0186             _truth_container->GetPrimaryParticleRange();
0187     //std::cout << "A2" << std::endl;
0188     for (PHG4TruthInfoContainer::ConstIterator truth_itr = range.first;
0189          truth_itr != range.second; ++truth_itr) {
0190
0191         PHG4Particle* g4particle = truth_itr->second;
0192         if(!g4particle) {
0193             LogDebug("");
0194             continue;
0195         }
0196         //std::cout << "B1" << std::endl;
0197
0198         SvtxTrack_FastSim* track = NULL;
0199
0200         //std::cout << "TRACKmap size " << _trackmap->size() << std::endl;
0201         for (SvtxTrackMap::ConstIter track_itr = _trackmap->begin();
0202              track_itr != _trackmap->end();
0203              track_itr++) {
0204           //std::cout << "TRACK * " << track_itr->first << std::endl;
0205           SvtxTrack_FastSim* temp = dynamic_cast<SvtxTrack_FastSim*>(track_itr->second);
0206         if(!temp) {
0207           std::cout << "ERROR CASTING PARTICLE!" << std::endl;
0208           continue;
0209         }
0210         //std::cout << " PARTICLE!" << std::endl;
0211
0212             if ((temp->get_truth_track_id() - g4particle->get_track_id()) == 0) {
0213                 track = temp;
0214             }
0215         }
0216
0217         //std::cout << "B2" << std::endl;
0218         gtrackID = g4particle->get_track_id();
0219         gflavor = g4particle->get_pid();
0220
0221         gpx = g4particle->get_px();
0222         gpy = g4particle->get_py();
0223         gpz = g4particle->get_pz();
0224         gpt = sqrt(gpx*gpx+gpy*gpy);
0225         gp = sqrt(gpx*gpx+gpy*gpy+gpz*gpz);
0226         gtheta = atan2(gpt,gpz);
0227         geta=-1.*log(tan(gtheta/2.));
0228         gphi = atan(gpy/gpx);
0229
0230         if (track) {
0231           //std::cout << "C1" << std::endl;
0232             trackID = track->get_id();
0233             charge = track->get_charge();
0234             nhits = track->size_clusters();
0235
0236             px = track->get_px();
0237             py = track->get_py();
0238             pz = track->get_pz();
0239             pt = sqrt(px*px+py*py);
0240             p = sqrt(px*px+py*py+pz*pz);
0241             theta = atan2(pt,pz);
0242             eta=-1.*log(tan(theta/2.));
0243             phi = atan(py/px);
0244             dca2d = track->get_dca2d();
0245
0246             TVector3 truth_mom(gpx,gpy,gpz);
0247             TVector3 reco_mom(px, py, pz);
0248             //std::cout << "C2" << std::endl;
0249
0250             _h2d_Delta_mom_vs_truth_mom->Fill(truth_mom.Mag(), (reco_mom.Mag()-truth_mom.Mag())/truth_mom.Mag());
0251             _h2d_Delta_mom_vs_truth_eta->Fill(truth_mom.Eta(), (reco_mom.Mag()-truth_mom.Mag())/truth_mom.Mag());
0252         }
0253         //std::cout << "B3" << std::endl;
0254
0255         _eval_tree_tracks->Fill();
0256
0257     }
0258     //std::cout << "A3" << std::endl;
0259
0260     return;
0261
0262 }
0263
0264 //----------------------------------------------------------------------------//
0265 //-- reset_variables():
0266 //--   Reset all the tree variables to their default values.
0267 //--   Needs to be called at the start of every event
0268 //----------------------------------------------------------------------------//
0269 void FastTrackingEval::reset_variables() {
0270     event = -9999;
0271
0272     //-- truth
0273     gtrackID = -9999;
0274     gflavor = -9999;
0275     gpx = -9999;
0276     gpy = -9999;
0277     gpz = -9999;
0278     gpt = -9999;
0279     gp  = -9999;
0280     gtheta = -9999;
0281     geta = -9999;
0282     gphi = -9999;
0283     gvx = -9999;
0284     gvy = -9999;
0285     gvz = -9999;
0286
0287     //-- reco
0288     trackID = -9999;
0289     charge = -9999;
0290     nhits = -9999;
0291     px = -9999;
0292     py = -9999;
0293     pz = -9999;
0294     pt = -9999;
0295     p  = -9999;
0296     theta = -9999;
0297     eta = -9999;
0298     phi = -9999;
0299     dca2d = -9999;
0300 }
0301
0302 //----------------------------------------------------------------------------//
0303 //-- GetNodes():
0304 //--   Get all the all the required nodes off the node tree
0305 //----------------------------------------------------------------------------//
0306 int FastTrackingEval::GetNodes(PHCompositeNode * topNode) {
0307     //DST objects
0308     //Truth container
0309     _truth_container = findNode::getClass<PHG4TruthInfoContainer>(topNode,
0310             "G4TruthInfo");
0311     if (!_truth_container && _event < 2) {
0312         cout << PHWHERE << " PHG4TruthInfoContainer node not found on node tree"
0313                 << endl;
0314         return Fun4AllReturnCodes::ABORTEVENT;
0315     }
0316
0317     _trackmap = findNode::getClass<SvtxTrackMap>(topNode,
0318             _trackmapname);
0319     //std::cout << _trackmapname.c_str() << std::endl;
0320     if (!_trackmap && _event < 2) {
0321         cout << PHWHERE << "SvtxTrackMap node with name "
0322              << _trackmapname
0323              <<" not found on node tree"
0324              << endl;
0325         return Fun4AllReturnCodes::ABORTEVENT;
0326     }
0327
0328     return Fun4AllReturnCodes::EVENT_OK;
0329 }
0330