packages/TrackingDiagnostics/BeamCrossingAnalysis.cc

0001 #include "BeamCrossingAnalysis.h"
0002
0003 #include <fun4all/Fun4AllReturnCodes.h>
0004 #include <phool/PHCompositeNode.h>
0005 #include <phool/getClass.h>
0006
0007 #include <trackbase/TrackVertexCrossingAssoc.h>
0008
0009 #include <utility>
0010
0011 #include <TFile.h>
0012 #include <TH1.h>
0013 #include <TNtuple.h>
0014
0015 BeamCrossingAnalysis::BeamCrossingAnalysis(const std::string& name)
0016   : SubsysReco(name)
0017 {
0018 }
0019
0020 int BeamCrossingAnalysis::InitRun(PHCompositeNode* topNode)
0021 {
0022   const char* cfilepath = filepath.c_str();
0023   fout = new TFile(cfilepath, "recreate");
0024   ntp_vertex = new TNtuple("ntp_vertex", "Vertex ntuple","crossing:event:vtxid:x:y:z:ntracks");
0025   ntp_track = new TNtuple("ntp_track", "Track ntuple", "crossing:event:trackid:quality:x:y:z:px:py:pz");
0026
0027   getNodes(topNode);
0028
0029   hcross = new TH1D("hcross", "Crossings", 1000, -150, 250);  // root histogram arguments: name,title,bins,minvalx,maxvalx
0030   hvertz = new TH1D("hvertz", "Vertex z", 1000, -20, 20);
0031   htrackz = new TH1D("htrackz", "Track z", 1000, -20, 20);
0032   hvertcross = new TH1D("hvertcross", "Vertex crossing", 1000, -150, 250);
0033   htrackcross = new TH1D("htrackcross", "Track crossing", 1000, -150, 250);
0034
0035   _event = 0;
0036
0037   return 0;
0038 }
0039
0040 int BeamCrossingAnalysis::process_event(PHCompositeNode* /**topNode*/)
0041 {
0042
0043   std::set<short int> crossing_set = m_track_vertex_crossing_map->getCrossings();
0044
0045   for (const auto& cross : crossing_set)
0046     {
0047       if(Verbosity() > 0)
0048     {
0049
0050       std::cout  << "Event " << _event << " Crossing " << cross << std::endl;
0051     }
0052
0053       // get all vertices for this crossing
0054       auto vtxit = m_track_vertex_crossing_map->getVertices(cross);
0055
0056       for (auto itr = vtxit.first; itr != vtxit.second; ++itr)
0057     {
0058       unsigned int vtxid = itr->second;
0059       const SvtxVertex* svtxVertex = m_vertexMap->get(vtxid);
0060       if(!svtxVertex) { continue; }
0061
0062       float vx=0;
0063       float vy=0;
0064       float vz=0;
0065           unsigned int ntracks = 0;
0066       vx = svtxVertex->get_x() ;
0067       vy = svtxVertex->get_y() ;
0068       vz = svtxVertex->get_z();
0069
0070       ntracks = svtxVertex->size_tracks();
0071
0072       // fill vertex ntuple
0073           float vertex_info[] = {(float) cross, (float) _event,  (float) vtxid, vx, vy, vz, (float) ntracks };
0074       ntp_vertex->Fill(vertex_info);
0075
0076       hvertz->Fill(vz);
0077       hvertcross->Fill(cross);
0078
0079       if(Verbosity() > 0)
0080         {
0081           std::cout << "    crossing  " << cross << " vtxid " << vtxid  << " ntracks " << ntracks
0082             << " x " << vx << " y " << vy << " z " << vz << " cm " <<  std::endl;
0083         }
0084     }
0085
0086        // get all tracks for this crossing
0087       unsigned int tcount = 0;
0088       auto trit = m_track_vertex_crossing_map->getTracks(cross);
0089       for (auto itr = trit.first; itr != trit.second; ++itr)
0090     {
0091       unsigned int  trid = itr->second;
0092
0093       const SvtxTrack *track = m_svtxTrackMap->get(trid);
0094           if (!track) { continue; }
0095
0096       float tx=0;
0097       float ty=0;
0098       float tz=0;
0099       float px=0;
0100       float py=0;
0101       float pz=0;
0102
0103
0104       tx = track->get_x() ;
0105       ty = track->get_y() ;
0106       tz = track->get_z() ;
0107
0108       px = track->get_px();
0109       py = track->get_py();
0110       pz = track->get_pz();
0111
0112       float quality = track->get_quality();
0113
0114       // fill track ntuple
0115           float track_info[] = {(float) cross, (float) _event, (float) trid, quality, tx, ty, tz, px, py, pz };
0116       ntp_track->Fill(track_info);
0117
0118       if(Verbosity() > 0)
0119         {
0120           std::cout << "    cross  " << cross << " evt " << _event << " trid " << trid
0121             << " qual " << quality << " x " << tx << " y " << ty << " z " << tz
0122             << " px " << px << " py " << py << " pz " << pz   << std::endl;
0123         }
0124
0125       // count tracks per crossing
0126       tcount++;
0127
0128       htrackz->Fill(tz);
0129       htrackcross->Fill(cross);
0130     }
0131
0132       if(Verbosity() > 0)
0133     {
0134       std::cout << "Crossing " << cross << " tcount " << tcount << std::endl;
0135     }
0136
0137       if(tcount > 0)
0138     {
0139       hcross->Fill(cross);
0140     }
0141
0142       // code here to identify single track crossings and crossings with no vertex (if different)
0143
0144
0145     }
0146
0147   _event++;
0148
0149   return 0;
0150 }
0151
0152
0153 int BeamCrossingAnalysis::getNodes(PHCompositeNode* topNode)
0154 {
0155   m_svtxTrackMap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0156   if (!m_svtxTrackMap)
0157     {
0158       std::cout << PHWHERE << "No SvtxTrackMap on node tree, exiting." << std::endl;
0159       return Fun4AllReturnCodes::ABORTEVENT;
0160     }
0161
0162   m_vertexMap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0163   if (!m_vertexMap)
0164     {
0165       std::cout << PHWHERE << "No vertexMap on node tree, exiting." << std::endl;
0166       return Fun4AllReturnCodes::ABORTEVENT;
0167     }
0168
0169   m_track_vertex_crossing_map = findNode::getClass<TrackVertexCrossingAssoc>(topNode,"TrackVertexCrossingAssocMap");
0170   if(!m_track_vertex_crossing_map)
0171     {
0172       std::cout << PHWHERE << "No TrackVertexCrossingAssocMap on node tree, exiting." << std::endl;
0173       return Fun4AllReturnCodes::ABORTEVENT;
0174     }
0175   return Fun4AllReturnCodes::EVENT_OK;
0176 }
0177
0178 int BeamCrossingAnalysis::End(PHCompositeNode* /**topNode*/)
0179 {
0180   fout->cd();
0181   ntp_vertex->Write();
0182   ntp_track->Write();
0183   hcross->Write();
0184   htrackz->Write();
0185   hvertz->Write();
0186   htrackcross->Write();
0187   hvertcross->Write();
0188   fout->Close();
0189
0190   return 0;
0191 }