analysis/SimpleTrackingAnalysis/SimpleTrackingAnalysis.C

0001 #include "SimpleTrackingAnalysis.h"
0002
0003 #include <phool/getClass.h>
0004 #include <fun4all/Fun4AllServer.h>
0005 #include <fun4all/Fun4AllReturnCodes.h>
0006
0007 #include <phool/PHCompositeNode.h>
0008 //#include <g4hough/PHG4HoughTransform.h> // CAUSES ERRORS DUE TO VertexFinder.h FAILING TO FIND <Eigen/LU>
0009
0010
0011 #include <g4main/PHG4TruthInfoContainer.h>
0012 #include <g4main/PHG4Particle.h>
0013 #include <g4main/PHG4VtxPoint.h>
0014
0015 #include <g4hough/SvtxVertexMap.h>
0016 #include <g4hough/SvtxVertex.h>
0017 #include <g4hough/SvtxTrackMap.h>
0018 #include <g4hough/SvtxTrack.h>
0019
0020 #include <g4eval/SvtxEvalStack.h>
0021 #include <g4eval/SvtxTrackEval.h>
0022 #include <g4eval/SvtxVertexEval.h>
0023 #include <g4eval/SvtxTruthEval.h>
0024
0025 // --- common to all calorimeters
0026 #include <calobase/RawTowerGeomContainer.h>
0027 #include <calobase/RawTowerContainer.h>
0028 #include <calobase/RawTower.h>
0029 #include <calobase/RawCluster.h>
0030 #include <calobase/RawClusterContainer.h>
0031 #include <g4eval/CaloEvalStack.h>
0032 #include <g4eval/CaloRawClusterEval.h>
0033 #include <g4eval/CaloRawTowerEval.h>
0034
0035 #include <TH1D.h>
0036 #include <TH2D.h>
0037 #include <TProfile2D.h>
0038
0039 #include <iostream>
0040 //#include <cassert>
0041
0042 using namespace std;
0043
0044
0045
0046 SimpleTrackingAnalysis::SimpleTrackingAnalysis(const string &name) : SubsysReco(name)
0047 {
0048   cout << "SimpleTrackingAnalysis class constructor called " << endl;
0049   nevents = 0;
0050   nerrors = 0;
0051   nwarnings = 0;
0052   nlayers = 7;
0053   verbosity = 0;
0054   magneticfield = 1.4; // default is 1.5 but Mike's tracking stuff is using 1.4 for now...
0055   //docalocuts = false;
0056 }
0057
0058
0059
0060 int SimpleTrackingAnalysis::Init(PHCompositeNode *topNode)
0061 {
0062
0063   // --------------------------------------------------------------------------------
0064   // --- This is the class initialization method
0065   // --- Get a pointer to the Fun4All server to register histograms
0066   // --- Histograms are declared in class header file
0067   // --- Histograms are instantiated here and then registered with the Fun4All server
0068   // --------------------------------------------------------------------------------
0069
0070   cout << "SimpleTrackingAnalysis::Init called" << endl;
0071
0072   // --- get instance of
0073   Fun4AllServer *se = Fun4AllServer::instance();
0074
0075
0076
0077   // --- special stuff...
0078
0079
0080
0081
0082   // --- additional tracking histograms for studying quality
0083
0084   _recopt_quality = new TH2D("recopt_quality", "", 20,0.0,10.0, 100,0.0,5.0);
0085   _recopt_quality_tracks_all = new TH2D("recopt_quality_tracks_all", "", 20,0.0,10.0, 100,0.0,5.0);
0086   _recopt_quality_tracks_recoWithin4Percent = new TH2D("recopt_quality_tracks_recoWithin4Percent", "", 20,0.0,10.0, 100,0.0,5.0);
0087   _truept_quality_particles_recoWithin4Percent = new TH2D("truept_quality_particles_recoWithin4Percent", "", 20,0.0,10.0, 100,0.0,5.0);
0088
0089   se->registerHisto(_recopt_quality);
0090   se->registerHisto(_recopt_quality_tracks_all);
0091   se->registerHisto(_recopt_quality_tracks_recoWithin4Percent);
0092   se->registerHisto(_truept_quality_particles_recoWithin4Percent);
0093
0094
0095
0096   // --- histograms over true pt, used for finding efficiencies
0097
0098   _truept_dca = new TH2D("truept_dca", "", 20,0.0,10.0, 200,-0.1,0.1);
0099   _truept_dptoverpt = new TH2D("truept_dptoverpt", "", 40,0.0,40.0, 200,-0.5,0.5);
0100   _truept_particles_leavingAllHits = new TH1D("truept_particles_leavingAllHits", "", 20,0.0,10.0);
0101   _truept_particles_recoWithExactHits = new TH1D("truept_particles_recoWithExactHits", "", 20,0.0,10.0);
0102   _truept_particles_recoWithin1Hit = new TH1D("truept_particles_recoWithin1Hit", "", 20,0.0,10.0);
0103   _truept_particles_recoWithin2Hits = new TH1D("truept_particles_recoWithin2Hits", "", 20,0.0,10.0);
0104   _truept_particles_recoWithin3Percent = new TH1D("truept_particles_recoWithin3Percent", "", 20,0.0,10.0);
0105   _truept_particles_recoWithin4Percent = new TH1D("truept_particles_recoWithin4Percent", "", 20,0.0,10.0);
0106   _truept_particles_recoWithin5Percent = new TH1D("truept_particles_recoWithin5Percent", "", 20,0.0,10.0);
0107
0108   se->registerHisto(_truept_dca);
0109   se->registerHisto(_truept_dptoverpt);
0110   se->registerHisto(_truept_particles_leavingAllHits);
0111   se->registerHisto(_truept_particles_recoWithExactHits);
0112   se->registerHisto(_truept_particles_recoWithin1Hit);
0113   se->registerHisto(_truept_particles_recoWithin2Hits);
0114   se->registerHisto(_truept_particles_recoWithin3Percent);
0115   se->registerHisto(_truept_particles_recoWithin4Percent);
0116   se->registerHisto(_truept_particles_recoWithin5Percent);
0117
0118
0119
0120   // --- (mostly) the same set of histograms over reconstructed pt, used for studying purity
0121
0122   _recopt_tracks_all = new TH1D("recopt_tracks_all", "", 20,0.0,10.0);
0123   _recopt_tracks_recoWithExactHits = new TH1D("recopt_tracks_recoWithExactHits", "", 20,0.0,10.0);
0124   _recopt_tracks_recoWithin1Hit = new TH1D("recopt_tracks_recoWithin1Hit", "", 20,0.0,10.0);
0125   _recopt_tracks_recoWithin2Hits = new TH1D("recopt_tracks_recoWithin2Hits", "", 20,0.0,10.0);
0126   _recopt_tracks_recoWithin3Percent = new TH1D("recopt_tracks_recoWithin3Percent", "", 20,0.0,10.0);
0127   _recopt_tracks_recoWithin4Percent = new TH1D("recopt_tracks_recoWithin4Percent", "", 20,0.0,10.0);
0128   _recopt_tracks_recoWithin5Percent = new TH1D("recopt_tracks_recoWithin5Percent", "", 20,0.0,10.0);
0129
0130   se->registerHisto(_recopt_tracks_all);
0131   se->registerHisto(_recopt_tracks_recoWithExactHits);
0132   se->registerHisto(_recopt_tracks_recoWithin1Hit);
0133   se->registerHisto(_recopt_tracks_recoWithin2Hits);
0134   se->registerHisto(_recopt_tracks_recoWithin3Percent);
0135   se->registerHisto(_recopt_tracks_recoWithin4Percent);
0136   se->registerHisto(_recopt_tracks_recoWithin5Percent);
0137
0138
0139
0140   // --- purity study with calorimeter cuts
0141
0142   th2d_recopt_tracks_withcalocuts_all = new TH2D(Form("th2d_recopt_tracks_withcalocuts_all"), "", 80,0.0,40.0, 20,0.0,2.0);
0143   th2d_recopt_tracks_withcalocuts_recoWithExactHits = new TH2D(Form("th2d_recopt_tracks_withcalocuts_recoWithExactHits"), "", 80,0.0,40.0, 20,0.0,2.0);
0144   th2d_recopt_tracks_withcalocuts_recoWithin1Hit = new TH2D(Form("th2d_recopt_tracks_withcalocuts_recoWithin1Hit"), "", 80,0.0,40.0, 20,0.0,2.0);
0145   th2d_recopt_tracks_withcalocuts_recoWithin2Hits = new TH2D(Form("th2d_recopt_tracks_withcalocuts_recoWithin2Hits"), "", 80,0.0,40.0, 20,0.0,2.0);
0146   th2d_recopt_tracks_withcalocuts_recoWithin3Percent = new TH2D(Form("th2d_recopt_tracks_withcalocuts_recoWithin3Percent"), "", 80,0.0,40.0, 20,0.0,2.0);
0147   th2d_recopt_tracks_withcalocuts_recoWithin4Percent = new TH2D(Form("th2d_recopt_tracks_withcalocuts_recoWithin4Percent"), "", 80,0.0,40.0, 20,0.0,2.0);
0148   th2d_recopt_tracks_withcalocuts_recoWithin5Percent = new TH2D(Form("th2d_recopt_tracks_withcalocuts_recoWithin5Percent"), "", 80,0.0,40.0, 20,0.0,2.0);
0149   th2d_recopt_tracks_withcalocuts_recoWithin1Sigma = new TH2D(Form("th2d_recopt_tracks_withcalocuts_recoWithin1Sigma"), "", 80,0.0,40.0, 20,0.0,2.0);
0150   th2d_recopt_tracks_withcalocuts_recoWithin2Sigma = new TH2D(Form("th2d_recopt_tracks_withcalocuts_recoWithin2Sigma"), "", 80,0.0,40.0, 20,0.0,2.0);
0151   th2d_recopt_tracks_withcalocuts_recoWithin3Sigma = new TH2D(Form("th2d_recopt_tracks_withcalocuts_recoWithin3Sigma"), "", 80,0.0,40.0, 20,0.0,2.0);
0152
0153   se->registerHisto(th2d_recopt_tracks_withcalocuts_all);
0154   se->registerHisto(th2d_recopt_tracks_withcalocuts_recoWithExactHits);
0155   se->registerHisto(th2d_recopt_tracks_withcalocuts_recoWithin1Hit);
0156   se->registerHisto(th2d_recopt_tracks_withcalocuts_recoWithin2Hits);
0157   se->registerHisto(th2d_recopt_tracks_withcalocuts_recoWithin3Percent);
0158   se->registerHisto(th2d_recopt_tracks_withcalocuts_recoWithin4Percent);
0159   se->registerHisto(th2d_recopt_tracks_withcalocuts_recoWithin5Percent);
0160   se->registerHisto(th2d_recopt_tracks_withcalocuts_recoWithin1Sigma);
0161   se->registerHisto(th2d_recopt_tracks_withcalocuts_recoWithin2Sigma);
0162   se->registerHisto(th2d_recopt_tracks_withcalocuts_recoWithin3Sigma);
0163
0164
0165
0166   // --- efficiency study with calorimter cuts
0167
0168   th2d_truept_particles_withcalocuts_leavingAllHits = new TH2D(Form("th2d_truept_particles_withcalocuts_leavingAllHits"), "", 80,0.0,40.0, 20,0.0,2.0);
0169   th2d_truept_particles_withcalocuts_recoWithExactHits = new TH2D(Form("th2d_truept_particles_withcalocuts_recoWithExactHits"), "", 80,0.0,40.0, 20,0.0,2.0);
0170   th2d_truept_particles_withcalocuts_recoWithin1Hit = new TH2D(Form("th2d_truept_particles_withcalocuts_recoWithin1Hit"), "", 80,0.0,40.0, 20,0.0,2.0);
0171   th2d_truept_particles_withcalocuts_recoWithin2Hits = new TH2D(Form("th2d_truept_particles_withcalocuts_recoWithin2Hits"), "", 80,0.0,40.0, 20,0.0,2.0);
0172   th2d_truept_particles_withcalocuts_recoWithin3Percent = new TH2D(Form("th2d_truept_particles_withcalocuts_recoWithin3Percent"), "", 80,0.0,40.0, 20,0.0,2.0);
0173   th2d_truept_particles_withcalocuts_recoWithin4Percent = new TH2D(Form("th2d_truept_particles_withcalocuts_recoWithin4Percent"), "", 80,0.0,40.0, 20,0.0,2.0);
0174   th2d_truept_particles_withcalocuts_recoWithin5Percent = new TH2D(Form("th2d_truept_particles_withcalocuts_recoWithin5Percent"), "", 80,0.0,40.0, 20,0.0,2.0);
0175   th2d_truept_particles_withcalocuts_recoWithin1Sigma = new TH2D(Form("th2d_truept_particles_withcalocuts_recoWithin1Sigma"), "", 80,0.0,40.0, 20,0.0,2.0);
0176   th2d_truept_particles_withcalocuts_recoWithin2Sigma = new TH2D(Form("th2d_truept_particles_withcalocuts_recoWithin2Sigma"), "", 80,0.0,40.0, 20,0.0,2.0);
0177   th2d_truept_particles_withcalocuts_recoWithin3Sigma = new TH2D(Form("th2d_truept_particles_withcalocuts_recoWithin3Sigma"), "", 80,0.0,40.0, 20,0.0,2.0);
0178
0179   se->registerHisto(th2d_truept_particles_withcalocuts_leavingAllHits);
0180   se->registerHisto(th2d_truept_particles_withcalocuts_recoWithExactHits);
0181   se->registerHisto(th2d_truept_particles_withcalocuts_recoWithin1Hit);
0182   se->registerHisto(th2d_truept_particles_withcalocuts_recoWithin2Hits);
0183   se->registerHisto(th2d_truept_particles_withcalocuts_recoWithin3Percent);
0184   se->registerHisto(th2d_truept_particles_withcalocuts_recoWithin4Percent);
0185   se->registerHisto(th2d_truept_particles_withcalocuts_recoWithin5Percent);
0186   se->registerHisto(th2d_truept_particles_withcalocuts_recoWithin1Sigma);
0187   se->registerHisto(th2d_truept_particles_withcalocuts_recoWithin2Sigma);
0188   se->registerHisto(th2d_truept_particles_withcalocuts_recoWithin3Sigma);
0189
0190
0191
0192   // --- new additional (eventual replacement?) histograms for purity study
0193
0194   th2d_reco_calo_nhits8 = new TH2D("th2d_reco_calo_nhits8", "", 80,0.0,40.0, 20,0.0,2.0);
0195   th2d_reco_calo_nhits7 = new TH2D("th2d_reco_calo_nhits7", "", 80,0.0,40.0, 20,0.0,2.0);
0196   th2d_reco_calo_nhits6 = new TH2D("th2d_reco_calo_nhits6", "", 80,0.0,40.0, 20,0.0,2.0);
0197   th2d_reco_calo_nhits5 = new TH2D("th2d_reco_calo_nhits5", "", 80,0.0,40.0, 20,0.0,2.0);
0198   th2d_reco_calo_nhits4 = new TH2D("th2d_reco_calo_nhits4", "", 80,0.0,40.0, 20,0.0,2.0);
0199   th2d_reco_calo_nhits3 = new TH2D("th2d_reco_calo_nhits3", "", 80,0.0,40.0, 20,0.0,2.0);
0200   th2d_reco_calo_nhits2 = new TH2D("th2d_reco_calo_nhits2", "", 80,0.0,40.0, 20,0.0,2.0);
0201   th2d_reco_calo_nhits1 = new TH2D("th2d_reco_calo_nhits1", "", 80,0.0,40.0, 20,0.0,2.0);
0202
0203   th2d_reco_calo_pt1sigma = new TH2D("th2d_reco_calo_pt1sigma", "", 80,0.0,40.0, 20,0.0,2.0);
0204   th2d_reco_calo_pt2sigma = new TH2D("th2d_reco_calo_pt2sigma", "", 80,0.0,40.0, 20,0.0,2.0);
0205   th2d_reco_calo_pt3sigma = new TH2D("th2d_reco_calo_pt3sigma", "", 80,0.0,40.0, 20,0.0,2.0);
0206   th2d_reco_calo_pt4sigma = new TH2D("th2d_reco_calo_pt4sigma", "", 80,0.0,40.0, 20,0.0,2.0);
0207   th2d_reco_calo_pt5sigma = new TH2D("th2d_reco_calo_pt5sigma", "", 80,0.0,40.0, 20,0.0,2.0);
0208   th2d_reco_calo_pt6sigma = new TH2D("th2d_reco_calo_pt6sigma", "", 80,0.0,40.0, 20,0.0,2.0);
0209
0210
0211
0212   // --- new additional (eventual replacement?) histograms for purity study
0213
0214   th2d_true_calo_nhits8 = new TH2D("th2d_true_calo_nhits8", "", 80,0.0,40.0, 20,0.0,2.0);
0215   th2d_true_calo_nhits7 = new TH2D("th2d_true_calo_nhits7", "", 80,0.0,40.0, 20,0.0,2.0);
0216   th2d_true_calo_nhits6 = new TH2D("th2d_true_calo_nhits6", "", 80,0.0,40.0, 20,0.0,2.0);
0217   th2d_true_calo_nhits5 = new TH2D("th2d_true_calo_nhits5", "", 80,0.0,40.0, 20,0.0,2.0);
0218   th2d_true_calo_nhits4 = new TH2D("th2d_true_calo_nhits4", "", 80,0.0,40.0, 20,0.0,2.0);
0219   th2d_true_calo_nhits3 = new TH2D("th2d_true_calo_nhits3", "", 80,0.0,40.0, 20,0.0,2.0);
0220   th2d_true_calo_nhits2 = new TH2D("th2d_true_calo_nhits2", "", 80,0.0,40.0, 20,0.0,2.0);
0221   th2d_true_calo_nhits1 = new TH2D("th2d_true_calo_nhits1", "", 80,0.0,40.0, 20,0.0,2.0);
0222
0223   th2d_true_calo_pt1sigma = new TH2D("th2d_true_calo_pt1sigma", "", 80,0.0,40.0, 20,0.0,2.0);
0224   th2d_true_calo_pt2sigma = new TH2D("th2d_true_calo_pt2sigma", "", 80,0.0,40.0, 20,0.0,2.0);
0225   th2d_true_calo_pt3sigma = new TH2D("th2d_true_calo_pt3sigma", "", 80,0.0,40.0, 20,0.0,2.0);
0226   th2d_true_calo_pt4sigma = new TH2D("th2d_true_calo_pt4sigma", "", 80,0.0,40.0, 20,0.0,2.0);
0227   th2d_true_calo_pt5sigma = new TH2D("th2d_true_calo_pt5sigma", "", 80,0.0,40.0, 20,0.0,2.0);
0228   th2d_true_calo_pt6sigma = new TH2D("th2d_true_calo_pt6sigma", "", 80,0.0,40.0, 20,0.0,2.0);
0229
0230
0231
0232   // --- vertex residual histograms
0233
0234   _dx_vertex = new TH1D("dx_vertex", "dx_vertex", 200,-0.03,0.03);
0235   _dy_vertex = new TH1D("dy_vertex", "dy_vertex", 200,-0.03,0.03);
0236   _dz_vertex = new TH1D("dz_vertex", "dz_vertex", 200,-0.03,0.03);
0237   se->registerHisto(_dy_vertex);
0238   se->registerHisto(_dx_vertex);
0239   se->registerHisto(_dz_vertex);
0240
0241   hmult = new TH1D("hmult","",5000,-0.5,4999.5);
0242   hmult_vertex = new TH1D("hmult_vertex","",5000,-0.5,4999.5);
0243   se->registerHisto(hmult);
0244   se->registerHisto(hmult_vertex);
0245
0246
0247   return 0;
0248
0249 }
0250
0251
0252
0253 int SimpleTrackingAnalysis::process_event(PHCompositeNode *topNode)
0254 {
0255
0256   // --- This is the class process_event method
0257   // --- This is where the bulk of the analysis is done
0258   // --- Here we get the various data nodes we need to do the analysis
0259   // --- Then we use variables (accessed through class methods) to perform calculations
0260
0261   if ( verbosity > -1 )
0262     {
0263       cout << endl;
0264       cout << "------------------------------------------------------------------------------------" << endl;
0265       cout << "Now processing event number " << nevents << endl; // would be good to add verbosity switch
0266     }
0267
0268   ++nevents; // You may as youtself, why ++nevents (pre-increment) rather
0269   // than nevents++ (post-increment)?  The short answer is performance.
0270   // For simple types it probably doesn't matter, but it can really help
0271   // for complex types (like the iterators below).
0272
0273
0274   // --- Truth level information
0275   PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode,"G4TruthInfo");
0276   if ( !truthinfo )
0277     {
0278       cerr << PHWHERE << " ERROR: Can't find G4TruthInfo" << endl;
0279       exit(-1);
0280     }
0281
0282   // --- SvtxTrackMap
0283   SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode,"SvtxTrackMap");
0284   if ( !trackmap )
0285     {
0286       cerr << PHWHERE << " ERROR: Can't find SvtxTrackMap" << endl;
0287       exit(-1);
0288     }
0289
0290   // --- SvtxVertexMap
0291   SvtxVertexMap* vertexmap = findNode::getClass<SvtxVertexMap>(topNode,"SvtxVertexMap");
0292   if ( !vertexmap )
0293     {
0294       cerr << PHWHERE << " ERROR: Can't find SvtxVertexMap" << endl;
0295       exit(-1);
0296     }
0297
0298
0299
0300   // --- Create SVTX eval stack
0301   SvtxEvalStack svtxevalstack(topNode);
0302   // --- Get evaluator objects from the eval stack
0303   SvtxVertexEval*   vertexeval = svtxevalstack.get_vertex_eval();
0304   SvtxTrackEval*     trackeval = svtxevalstack.get_track_eval();
0305   SvtxTruthEval*     trutheval = svtxevalstack.get_truth_eval();
0306
0307
0308
0309   if ( verbosity > 0 ) cout << "Now going to loop over truth partcles..." << endl; // need verbosity switch
0310
0311   // --- Loop over all truth particles
0312   PHG4TruthInfoContainer::Range range = truthinfo->GetPrimaryParticleRange();
0313   for ( PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter )
0314     {
0315       PHG4Particle* g4particle = iter->second; // You may ask yourself, why second?
0316       // In C++ the iterator is a map, which has two members
0317       // first is the key (analogous the index of an arry),
0318       // second is the value (analogous to the value stored for the array index)
0319       int particleID = g4particle->get_pid();
0320
0321       if ( trutheval->get_embed(g4particle) <= 0 && fabs(particleID) == 11 && verbosity > 0 )
0322         {
0323           cout << "NON EMBEDDED ELECTRON!!!  WHEE!!! " << particleID << " " << iter->first << endl;
0324         }
0325
0326       if ( trutheval->get_embed(g4particle) <= 0 ) continue; // only look at embedded particles // no good for hits files
0327       bool iselectron = fabs(particleID) == 11;
0328       bool ispion = fabs(particleID) == 211;
0329       if ( verbosity > 0 ) cout << "embedded particle ID is " << particleID << " ispion " << ispion << " iselectron " << iselectron << " " << iter->first << endl;
0330
0331       set<PHG4Hit*> g4hits = trutheval->all_truth_hits(g4particle);
0332       float ng4hits = g4hits.size();
0333
0334       float truept = sqrt(pow(g4particle->get_px(),2)+pow(g4particle->get_py(),2));
0335       float true_energy = g4particle->get_e();
0336
0337       // --- Get the reconsructed SvtxTrack based on the best candidate from the truth info
0338       SvtxTrack* track = trackeval->best_track_from(g4particle);
0339       if (!track) continue;
0340       float recopt = track->get_pt();
0341       float recop = track->get_p();
0342
0343       if ( verbosity > 0 )
0344     {
0345       cout << "truept is " << truept << endl;
0346       cout << "recopt is " << recopt << endl;
0347       cout << "true energy is " << true_energy << endl;
0348     }
0349
0350       // --- energy variables directly from the track object
0351       float emc_energy_track = track->get_cal_energy_3x3(SvtxTrack::CEMC);
0352       float hci_energy_track = track->get_cal_energy_3x3(SvtxTrack::HCALIN);
0353       float hco_energy_track = track->get_cal_energy_3x3(SvtxTrack::HCALOUT);
0354
0355       if ( verbosity > 0 )
0356     {
0357       cout << "emc_energy_track is " << emc_energy_track << endl;
0358       cout << "hci_energy_track is " << hci_energy_track << endl;
0359       cout << "hco_energy_track is " << hco_energy_track << endl;
0360     }
0361
0362       // -------------------------------------------------------------------------------------
0363       // --- IMPORTANT NOTE: according to Jin, dphi and deta will not work correctly in HIJING
0364
0365       float emc_dphi_track = track->get_cal_dphi(SvtxTrack::CEMC);
0366       float hci_dphi_track = track->get_cal_dphi(SvtxTrack::HCALIN);
0367       float hco_dphi_track = track->get_cal_dphi(SvtxTrack::HCALOUT);
0368
0369       float emc_deta_track = track->get_cal_deta(SvtxTrack::CEMC);
0370       float hci_deta_track = track->get_cal_deta(SvtxTrack::HCALIN);
0371       float hco_deta_track = track->get_cal_deta(SvtxTrack::HCALOUT);
0372
0373       float assoc_dphi = 0.1; // adjust as needed, consider class set method
0374       float assoc_deta = 0.1; // adjust as needed, consider class set method
0375       bool good_emc_assoc = fabs(emc_dphi_track) < assoc_dphi && fabs(emc_deta_track) < assoc_deta;
0376       bool good_hci_assoc = fabs(hci_dphi_track) < assoc_dphi && fabs(hci_deta_track) < assoc_deta;
0377       bool good_hco_assoc = fabs(hco_dphi_track) < assoc_dphi && fabs(hco_deta_track) < assoc_deta;
0378
0379       // ------------------------------------------------------------------------------------------
0380
0381
0382       float hct_energy_track = 0;
0383       if ( hci_energy_track >= 0 ) hct_energy_track += hci_energy_track;
0384       if ( hco_energy_track >= 0 ) hct_energy_track += hco_energy_track;
0385
0386       float total_energy_dumb = 0;
0387       if ( emc_energy_track >= 0 ) total_energy_dumb += emc_energy_track;
0388       if ( hci_energy_track >= 0 ) total_energy_dumb += hci_energy_track;
0389       if ( hco_energy_track >= 0 ) total_energy_dumb += hco_energy_track;
0390
0391       float total_energy_smart = 0;
0392       if ( good_emc_assoc ) total_energy_smart += emc_energy_track;
0393       if ( good_hci_assoc ) total_energy_smart += hci_energy_track;
0394       if ( good_hco_assoc ) total_energy_smart += hco_energy_track;
0395
0396
0397
0398       // ----------------------------------------------------------------------
0399       // ----------------------------------------------------------------------
0400       // ----------------------------------------------------------------------
0401
0402       //cout << "starting the main part of the truth analysis" << endl;
0403
0404       // examine truth particles that leave all (7 or 8 depending on design) detector hits
0405       if ( ng4hits == nlayers )
0406     {
0407       _truept_particles_leavingAllHits->Fill(truept);
0408
0409       unsigned int nfromtruth = trackeval->get_nclusters_contribution(track,g4particle);
0410
0411       unsigned int ndiff = abs((int)nfromtruth-(int)nlayers);
0412       if ( ndiff <= 2 ) _truept_particles_recoWithin2Hits->Fill(truept);
0413       if ( ndiff <= 1 ) _truept_particles_recoWithin1Hit->Fill(truept);
0414       if ( ndiff == 0 ) _truept_particles_recoWithExactHits->Fill(truept);
0415
0416       float diff = fabs(recopt-truept)/truept;
0417       if ( diff < 0.05 ) _truept_particles_recoWithin5Percent->Fill(truept);
0418       if ( diff < 0.04 )
0419         {
0420           _truept_particles_recoWithin4Percent->Fill(truept);
0421           _truept_quality_particles_recoWithin4Percent->Fill(truept,track->get_quality());
0422         }
0423       if ( diff < 0.03 ) _truept_particles_recoWithin3Percent->Fill(truept);
0424
0425           double good_energy = total_energy_dumb - 3.14;
0426
0427
0428           double eoverp = good_energy/recop;
0429           double sigmapt = 0.011 + 0.0008*recopt;
0430           th2d_truept_particles_withcalocuts_leavingAllHits->Fill(truept,eoverp);
0431           if ( ndiff <= 2 ) th2d_truept_particles_withcalocuts_recoWithin2Hits->Fill(truept,eoverp);
0432           if ( ndiff <= 1 ) th2d_truept_particles_withcalocuts_recoWithin1Hit->Fill(truept,eoverp);
0433           if ( ndiff == 0 ) th2d_truept_particles_withcalocuts_recoWithExactHits->Fill(truept,eoverp);
0434           if ( diff < 0.05 ) th2d_truept_particles_withcalocuts_recoWithin5Percent->Fill(truept,eoverp);
0435           if ( diff < 0.04 ) th2d_truept_particles_withcalocuts_recoWithin4Percent->Fill(truept,eoverp);
0436           if ( diff < 0.03 ) th2d_truept_particles_withcalocuts_recoWithin3Percent->Fill(truept,eoverp);
0437           if ( diff < 1.0*sigmapt ) th2d_truept_particles_withcalocuts_recoWithin1Sigma->Fill(recopt,eoverp);
0438           if ( diff < 2.0*sigmapt ) th2d_truept_particles_withcalocuts_recoWithin2Sigma->Fill(recopt,eoverp);
0439           if ( diff < 3.0*sigmapt ) th2d_truept_particles_withcalocuts_recoWithin3Sigma->Fill(recopt,eoverp);
0440
0441
0442     } // end of requirement of ng4hits == nlayers
0443
0444     } // end of loop over truth particles
0445
0446
0447
0448   // loop over all reco particles
0449   int ntracks = 0;
0450   for ( SvtxTrackMap::Iter iter = trackmap->begin(); iter != trackmap->end(); ++iter )
0451     {
0452
0453       // --- Get the StxTrack object (from the iterator)
0454       SvtxTrack* track = iter->second;
0455       float recopt = track->get_pt();
0456       float recop = track->get_p();
0457
0458       // --- Get the truth particle from the evaluator
0459       PHG4Particle* g4particle = trackeval->max_truth_particle_by_nclusters(track);
0460       float truept = sqrt(pow(g4particle->get_px(),2)+pow(g4particle->get_py(),2));
0461       int particleID = g4particle->get_pid();
0462       if ( verbosity > 5 ) cout << "particle ID is " << particleID << endl;
0463       bool iselectron = fabs(particleID) == 11;
0464       bool ispion = fabs(particleID) == 211;
0465
0466       // ---------------------
0467       // --- calorimeter stuff
0468       // ---------------------
0469
0470       // --- get the energy values directly from the track
0471       float emc_energy_track = track->get_cal_energy_3x3(SvtxTrack::CEMC);
0472       float hci_energy_track = track->get_cal_energy_3x3(SvtxTrack::HCALIN);
0473       float hco_energy_track = track->get_cal_energy_3x3(SvtxTrack::HCALOUT);
0474
0475       float total_energy = 0;
0476       if ( emc_energy_track > 0 ) total_energy += emc_energy_track;
0477       if ( hci_energy_track > 0 ) total_energy += hci_energy_track;
0478       if ( hco_energy_track > 0 ) total_energy += hco_energy_track;
0479
0480       if ( verbosity > 2 ) cout << "total calo energy is " << total_energy << endl;
0481
0482       if (trutheval->get_embed(g4particle) > 0)
0483     {
0484       // embedded results (quality or performance measures)
0485       _truept_dptoverpt->Fill(truept,(recopt-truept)/truept);
0486       _truept_dca->Fill(truept,track->get_dca2d());
0487       _recopt_quality->Fill(recopt,track->get_quality());
0488           if ( verbosity > 0 ) cout << "embedded particle ID is " << particleID << " ispion " << ispion << " iselectron " << iselectron << endl;
0489           // ---
0490     } // end if (embedded results)
0491       else
0492     {
0493           // electron and pion (hadron) id
0494
0495       // non-embedded results (purity measures)
0496       _recopt_tracks_all->Fill(recopt);
0497       _recopt_quality_tracks_all->Fill(recopt,track->get_quality());
0498
0499       unsigned int nfromtruth = trackeval->get_nclusters_contribution(track,g4particle);
0500
0501       unsigned int ndiff = abs((int)nfromtruth-(int)nlayers);
0502       if ( ndiff <= 2 ) _recopt_tracks_recoWithin2Hits->Fill(recopt);
0503       if ( ndiff <= 1 ) _recopt_tracks_recoWithin1Hit->Fill(recopt);
0504       if ( ndiff == 0 ) _recopt_tracks_recoWithExactHits->Fill(recopt);
0505
0506       float diff = fabs(recopt-truept)/truept;
0507       if ( diff < 0.05 ) _recopt_tracks_recoWithin5Percent->Fill(recopt);
0508       if ( diff < 0.04 )
0509         {
0510           _recopt_tracks_recoWithin4Percent->Fill(recopt);
0511           _recopt_quality_tracks_recoWithin4Percent->Fill(recopt,track->get_quality());
0512         }
0513       if ( diff < 0.03 ) _recopt_tracks_recoWithin3Percent->Fill(recopt);
0514
0515
0516       // --------------------------------------
0517       // --- same but now with calorimeter cuts
0518       // --------------------------------------
0519
0520           double good_energy = total_energy - 3.14;
0521
0522           double eoverp = good_energy/recop;
0523           double sigmapt = 0.011 + 0.0008*recopt;
0524           th2d_recopt_tracks_withcalocuts_all->Fill(recopt,eoverp);
0525           if ( ndiff <= 2 ) th2d_recopt_tracks_withcalocuts_recoWithin2Hits->Fill(recopt,eoverp);
0526           if ( ndiff <= 1 ) th2d_recopt_tracks_withcalocuts_recoWithin1Hit->Fill(recopt,eoverp);
0527           if ( ndiff == 0 ) th2d_recopt_tracks_withcalocuts_recoWithExactHits->Fill(recopt,eoverp);
0528           if ( diff < 0.05 ) th2d_recopt_tracks_withcalocuts_recoWithin5Percent->Fill(recopt,eoverp);
0529           if ( diff < 0.04 ) th2d_recopt_tracks_withcalocuts_recoWithin4Percent->Fill(recopt,eoverp);
0530           if ( diff < 0.03 ) th2d_recopt_tracks_withcalocuts_recoWithin3Percent->Fill(recopt,eoverp);
0531           if ( diff < 1.0*sigmapt ) th2d_recopt_tracks_withcalocuts_recoWithin1Sigma->Fill(recopt,eoverp);
0532           if ( diff < 2.0*sigmapt ) th2d_recopt_tracks_withcalocuts_recoWithin2Sigma->Fill(recopt,eoverp);
0533           if ( diff < 3.0*sigmapt ) th2d_recopt_tracks_withcalocuts_recoWithin3Sigma->Fill(recopt,eoverp);
0534
0535       // --- done with reco tracks
0536
0537     } // else (non-embedded results)
0538
0539       ++ntracks;
0540     } // loop over reco tracks
0541
0542
0543   hmult->Fill(ntracks);
0544
0545   // --- Get the leading vertex
0546   SvtxVertex* maxvertex = NULL;
0547   unsigned int maxtracks = 0;
0548   for ( SvtxVertexMap::Iter iter = vertexmap->begin(); iter != vertexmap->end(); ++iter )
0549     {
0550       SvtxVertex* vertex = iter->second;
0551       if ( vertex->size_tracks() > maxtracks )
0552     {
0553       maxvertex = vertex;
0554       maxtracks = vertex->size_tracks();
0555     }
0556     }
0557   if ( !maxvertex )
0558     {
0559       cerr << PHWHERE << " ERROR: cannot get reconstructed vertex (event number " << nevents << ")" << endl;
0560       ++nerrors;
0561       return Fun4AllReturnCodes::DISCARDEVENT;
0562     }
0563
0564   // --- Get the coordinates for the vertex from the evaluator
0565   PHG4VtxPoint* point = vertexeval->max_truth_point_by_ntracks(maxvertex);
0566   if ( !point )
0567     {
0568       cerr << PHWHERE << " ERROR: cannot get truth vertex (event number " << nevents << ")" << endl;
0569       ++nerrors;
0570       return Fun4AllReturnCodes::DISCARDEVENT;
0571     }
0572   _dx_vertex->Fill(maxvertex->get_x() - point->get_x());
0573   _dy_vertex->Fill(maxvertex->get_y() - point->get_y());
0574   _dz_vertex->Fill(maxvertex->get_z() - point->get_z());
0575
0576   hmult_vertex->Fill(ntracks);
0577
0578
0579
0580   return Fun4AllReturnCodes::EVENT_OK;
0581
0582 }
0583
0584
0585
0586 int SimpleTrackingAnalysis::End(PHCompositeNode *topNode)
0587 {
0588   cout << "End called, " << nevents << " events processed with " << nerrors << " errors and " << nwarnings << " warnings." << endl;
0589   return 0;
0590 }
0591