GenFitTest/AnaSvtxTracksForGenFit/AnaSvtxTracksForGenFit.cc

0001 ////////////////////////////////////////////////////////////////////////////////
0002 //
0003 // This module is desgined to grab svtx tracks and put truth and cluster
0004 // information into a TTree for GenFit testing
0005 //
0006 ////////////////////////////////////////////////////////////////////////////////
0007 //
0008 // Darren McGlinchey
0009 // 1 Apr 2016
0010 //
0011 ////////////////////////////////////////////////////////////////////////////////
0012
0013
0014 #include "AnaSvtxTracksForGenFit.h"
0015
0016 #include <phool/phool.h>
0017 #include <phool/getClass.h>
0018 #include <g4main/PHG4HitContainer.h>
0019 #include <g4main/PHG4TruthInfoContainer.h>
0020 #include <g4main/PHG4Particle.h>
0021 #include <g4main/PHG4Hit.h>
0022 #include <g4main/PHG4VtxPoint.h>
0023 #include <fun4all/PHTFileServer.h>
0024 #include <fun4all/Fun4AllServer.h>
0025
0026 #include <g4hough/SvtxVertexMap.h>
0027 #include <g4hough/SvtxVertex.h>
0028 #include <g4hough/SvtxTrackMap.h>
0029 #include <g4hough/SvtxTrack.h>
0030 #include <g4hough/SvtxClusterMap.h>
0031 #include <g4hough/SvtxCluster.h>
0032 #include <g4hough/SvtxHitMap.h>
0033 #include <g4hough/SvtxHit.h>
0034
0035 #include <g4eval/SvtxEvalStack.h>
0036 #include <g4eval/SvtxTrackEval.h>
0037 #include <g4eval/SvtxClusterEval.h>
0038 #include <g4eval/SvtxTruthEval.h>
0039 #include <g4eval/SvtxVertexEval.h>
0040 #include <g4eval/SvtxHitEval.h>
0041
0042 #include <TTree.h>
0043
0044 #include <iostream>
0045
0046 using namespace std;
0047
0048 //----------------------------------------------------------------------------//
0049 //-- Constructor:
0050 //--  simple initialization
0051 //----------------------------------------------------------------------------//
0052 AnaSvtxTracksForGenFit::AnaSvtxTracksForGenFit(const string &name):
0053   SubsysReco( name ),
0054   _flags( NONE ),
0055   _tracks( NULL ),
0056   _svtxevalstack( NULL )
0057 {
0058   //initialize
0059   _event = 0;
0060   _outfile = "AnaSvtxTracksForGenFit.root";
0061 }
0062
0063 //----------------------------------------------------------------------------//
0064 //-- Init():
0065 //--   Intialize all histograms, trees, and ntuples
0066 //----------------------------------------------------------------------------//
0067 int AnaSvtxTracksForGenFit::Init(PHCompositeNode *topNode)
0068 {
0069   cout << PHWHERE << " Openning file " << _outfile << endl;
0070   PHTFileServer::get().open( _outfile, "RECREATE");
0071
0072
0073   // create TTree
0074   _tracks = new TTree("tracks", "Svtx Tracks");
0075   _tracks->Branch("event", &event, "event/I");
0076   _tracks->Branch("gtrackID", &gtrackID, "gtrackID/I");
0077   _tracks->Branch("gflavor", &gflavor, "gflavor/I");
0078   _tracks->Branch("gpx", &gpx, "gpx/F");
0079   _tracks->Branch("gpy", &gpy, "gpy/F");
0080   _tracks->Branch("gpz", &gpz, "gpz/F");
0081   _tracks->Branch("gvx", &gvx, "gvx/F");
0082   _tracks->Branch("gvy", &gvy, "gvy/F");
0083   _tracks->Branch("gvz", &gvz, "gvz/F");
0084   _tracks->Branch("trackID", &trackID, "trackID/I");
0085   _tracks->Branch("charge", &charge, "charge/I");
0086   _tracks->Branch("nhits", &nhits, "nhits/I");
0087   _tracks->Branch("px", &px, "px/F");
0088   _tracks->Branch("py", &py, "py/F");
0089   _tracks->Branch("pz", &pz, "pz/F");
0090   _tracks->Branch("dca2d", &dca2d, "dca2d/F");
0091   _tracks->Branch("clusterID", &clusterID, "clusterID[nhits]/I");
0092   _tracks->Branch("layer", &layer, "layer[nhits]/I");
0093   _tracks->Branch("x", &x, "x[nhits]/F");
0094   _tracks->Branch("y", &y, "y[nhits]/F");
0095   _tracks->Branch("z", &z, "z[nhits]/F");
0096   _tracks->Branch("size_dphi", &size_dphi, "size_dphi[nhits]/F");
0097   _tracks->Branch("size_dz", &size_dz, "size_dz[nhits]/F");
0098
0099
0100   return 0;
0101 }
0102
0103 //----------------------------------------------------------------------------//
0104 //-- process_event():
0105 //--   Call user instructions for every event.
0106 //--   This function contains the analysis structure.
0107 //----------------------------------------------------------------------------//
0108 int AnaSvtxTracksForGenFit::process_event(PHCompositeNode *topNode)
0109 {
0110   _event++;
0111   if (_event % 1000 == 0)
0112     cout << PHWHERE << "Events processed: " << _event << endl;
0113
0114   GetNodes(topNode);
0115
0116   if (!_svtxevalstack) {
0117     _svtxevalstack = new SvtxEvalStack(topNode);
0118     _svtxevalstack->set_strict(false);
0119     _svtxevalstack->set_verbosity(verbosity + 1);
0120   } else {
0121     _svtxevalstack->next_event(topNode);
0122   }
0123
0124   fill_tree(topNode);
0125
0126   return 0;
0127 }
0128
0129 //----------------------------------------------------------------------------//
0130 //-- End():
0131 //--   End method, wrap everything up
0132 //----------------------------------------------------------------------------//
0133 int AnaSvtxTracksForGenFit::End(PHCompositeNode *topNode)
0134 {
0135
0136   PHTFileServer::get().cd( _outfile );
0137
0138   _tracks->Write();
0139
0140   if (_svtxevalstack) delete _svtxevalstack;
0141
0142   return 0;
0143 }
0144
0145
0146 //----------------------------------------------------------------------------//
0147 //-- fill_tree():
0148 //--   Fill the trees with truth, track fit, and cluster information
0149 //----------------------------------------------------------------------------//
0150 void AnaSvtxTracksForGenFit::fill_tree(PHCompositeNode *topNode)
0151 {
0152   // Make sure to reset all the TTree variables before trying to set them.
0153   reset_variables();
0154
0155   // get evaluators
0156   SvtxTrackEval*     trackeval = _svtxevalstack->get_track_eval();
0157   SvtxTruthEval*     trutheval = _svtxevalstack->get_truth_eval();
0158
0159   if (_truth_container)
0160   {
0161
0162     PHG4TruthInfoContainer::ConstRange range =
0163       _truth_container->GetPrimaryParticleRange();
0164
0165     for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
0166          iter != range.second;
0167          ++iter)
0168     {
0169
0170       PHG4Particle* g4particle = iter->second;
0171
0172       gtrackID = g4particle->get_track_id();
0173       gflavor  = g4particle->get_pid();
0174
0175       PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
0176       gvx      = vtx->get_x();
0177       gvy      = vtx->get_y();
0178       gvz      = vtx->get_z();
0179
0180             gpx          = g4particle->get_px();
0181             gpy          = g4particle->get_py();
0182             gpz          = g4particle->get_pz();
0183
0184       SvtxTrack* track = trackeval->best_track_from(g4particle);
0185       if (track)
0186       {
0187         trackID   = track->get_id();
0188         charge    = track->get_charge();
0189         nhits     = track->size_clusters();
0190         px        = track->get_px();
0191         py        = track->get_py();
0192         pz        = track->get_pz();
0193         dca2d     = track->get_dca2d();
0194
0195
0196         int iclus = 0;
0197         for (SvtxTrack::ConstClusterIter iter = track->begin_clusters();
0198              iter != track->end_clusters();
0199              ++iter)
0200         {
0201           unsigned int cluster_id = *iter;
0202           SvtxCluster* cluster = _clustermap->get(cluster_id);
0203           unsigned int l = cluster->get_layer();
0204
0205           clusterID[iclus] = (int)cluster_id;
0206           layer[iclus] = (int)l;
0207           x[iclus] = cluster->get_x();
0208           y[iclus] = cluster->get_y();
0209           z[iclus] = cluster->get_z();
0210           size_dphi[iclus] = cluster->get_phi_size();
0211           size_dz[iclus] = cluster->get_z_size();
0212
0213           ++iclus;
0214         }
0215
0216       } // if(track)
0217     } // for( iter)
0218   } //if (_truth_container)
0219
0220   _tracks->Fill();
0221   return;
0222
0223 }
0224
0225 //----------------------------------------------------------------------------//
0226 //-- reset_variables():
0227 //--   Reset all the tree variables to their default values.
0228 //--   Needs to be called at the start of every event
0229 //----------------------------------------------------------------------------//
0230 void AnaSvtxTracksForGenFit::reset_variables()
0231 {
0232   event = -9999;
0233
0234   //-- truth
0235   gtrackID = -9999;
0236   gflavor = -9999;
0237   gpx = -9999;
0238   gpy = -9999;
0239   gpz = -9999;
0240   gvx = -9999;
0241   gvy = -9999;
0242   gvz = -9999;
0243
0244   //-- reco
0245   trackID = -9999;
0246   charge = -9999;
0247   nhits = -9999;
0248   px = -9999;
0249   py = -9999;
0250   pz = -9999;
0251   dca2d = -9999;
0252
0253   //-- clusters
0254   for (int i = 0; i < 7; i++)
0255   {
0256     clusterID[i] = -9999;
0257     layer[i] = -9999;
0258     x[i] = -9999;
0259     y[i] = -9999;
0260     z[i] = -9999;
0261     size_dphi[i] = -9999;
0262     size_dz[i] = -9999;
0263   }
0264
0265 }
0266
0267 //----------------------------------------------------------------------------//
0268 //-- GetNodes():
0269 //--   Get all the all the required nodes off the node tree
0270 //----------------------------------------------------------------------------//
0271 void AnaSvtxTracksForGenFit::GetNodes(PHCompositeNode * topNode)
0272 {
0273   //DST objects
0274   //Truth container
0275   _truth_container = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0276   if (!_truth_container && _event < 2)
0277   {
0278     cout << PHWHERE
0279          << " PHG4TruthInfoContainer node not found on node tree"
0280          << endl;
0281   }
0282
0283   //Svtx Clusters
0284   _clustermap = findNode::getClass<SvtxClusterMap>(topNode, "SvtxClusterMap");
0285   if (!_clustermap && _event < 2)
0286   {
0287     cout << PHWHERE
0288          << " SvtxClusterMap node not found on node tree"
0289          << endl;
0290   }
0291
0292 }
0293
0294
0295
0296