TPC-GemGain/src/TPCGemGainCalb.cc

0001 #include "TPCGemGainCalb.h"
0002
0003 #include <fun4all/Fun4AllReturnCodes.h>
0004
0005 #include <phool/PHCompositeNode.h>
0006 #include <phool/getClass.h>
0007 #include <phool/phool.h>
0008 #include <TNtuple.h>
0009 #include <TFile.h>
0010
0011 /// Tracking includes
0012 #include <trackbase/TrkrDefs.h>
0013 #include <trackbase/TrkrClusterv2.h>
0014 #include <trackbase/TrkrClusterContainer.h>
0015 #include <trackbase/TrkrClusterHitAssocv2.h>
0016 #include <trackbase/TrkrHitTruthAssoc.h>
0017 #include <trackbase_historic/SvtxTrack_v2.h>
0018 #include <trackbase_historic/SvtxTrackMap.h>
0019 #include <trackbase_historic/SvtxVertexMap.h>
0020
0021
0022 using namespace std;
0023
0024 //____________________________________________________________________________..
0025 TPCGemGainCalb::TPCGemGainCalb(const std::string &name):
0026  SubsysReco(name)
0027 {
0028
0029   fout = new TFile("readback_ntuple.root","RECREATE");
0030
0031   ntp = new TNtuple("ntp", "ntp", "pt:x:y:z:dcaxy:dcaz:vtxid:nclus:qual");
0032
0033   //cout << "TPCGemGainCalb::TPCGemGainCalb(const std::string &name) Calling ctor" << endl;
0034 }
0035
0036 //____________________________________________________________________________..
0037 TPCGemGainCalb::~TPCGemGainCalb()
0038 {
0039
0040 }
0041
0042 //____________________________________________________________________________..
0043 int TPCGemGainCalb::InitRun(PHCompositeNode *topNode)
0044 {
0045
0046   int ret = GetNodes(topNode);
0047
0048   return ret;
0049 }
0050
0051 //____________________________________________________________________________..
0052 int TPCGemGainCalb::process_event(PHCompositeNode *topNode)
0053 {
0054   if (Verbosity() >= 1)
0055     cout << "TPCGemGainCalb::process_event(PHCompositeNode *topNode) Processing Event" << endl;
0056
0057   double vertex_x = 0;
0058   double vertex_y = 0;
0059   double vertex_z = 0;
0060
0061   // The vertex map
0062   std::cout << "Vertex map has " << _vertex_map->size() << " entries" << std::endl;
0063   for (auto vert_iter = _vertex_map->begin();
0064        vert_iter != _vertex_map->end();
0065        ++vert_iter)
0066     {
0067       auto vertex = vert_iter->second;
0068       vertex_x = vertex->get_x();
0069       vertex_y = vertex->get_y();
0070       vertex_z = vertex->get_z();
0071
0072       std::cout << " vertex (x,y,z) = " << vertex_x << " , " << vertex_y << " , " << vertex_z << std::endl;
0073     }
0074
0075
0076   // Loop over all SvtxTracks
0077
0078   int ntracks = 0;
0079   double mean_pt = 0;
0080   double mean_clusters = 0;
0081
0082   float pt, x, y, z, dcaxy, dcaz, vtxid, nclus, qual;
0083   std::cout << "Track map has " << _track_map->size() << " entries" << std::endl;
0084   for (auto phtrk_iter = _track_map->begin();
0085        phtrk_iter != _track_map->end();
0086        ++phtrk_iter)
0087     {
0088       _track = phtrk_iter->second;
0089
0090       pt = _track->get_pt();
0091       x = _track->get_x()  ;
0092       y = _track->get_y();
0093       z = _track->get_z();
0094       dcaxy = _track->get_dca3d_xy();
0095       dcaz = _track->get_dca3d_z();
0096       vtxid = (float) _track->get_vertex_id();
0097       nclus = (float) _track-> size_cluster_keys();
0098       qual = _track->get_chisq() / (float) _track->get_ndf();
0099
0100       if(Verbosity() > 0)
0101       std::cout
0102     << __LINE__
0103     << ": Processing itrack: " << phtrk_iter->first
0104     << ": nclus: " << nclus
0105     << ": pT " << pt
0106     << " : x,y,z " << x
0107     << ", " << y
0108     << ", " << z
0109     << ", dcaxy " << dcaxy
0110     << ", dcaz " << dcaz
0111     << ", vtxid " << vtxid
0112     << endl;
0113
0114       ntp->Fill(_track->get_pt(),
0115         _track->get_x(),
0116         _track->get_y(),
0117         _track->get_z(),
0118         _track->get_dca3d_xy(),
0119         _track->get_dca3d_z(),
0120         (float) _track->get_vertex_id(),
0121         (float) _track->size_cluster_keys(),
0122         qual
0123         );
0124
0125       ntracks ++;
0126       mean_pt += _track->get_pt();
0127       mean_clusters += _track->size_cluster_keys();
0128
0129       double cluster_avge_adc = 0.0;
0130       double cluster_avge_wt = 0.0;
0131
0132       // loop over associated clusters to get hits for track
0133       for (SvtxTrack::ConstClusterKeyIter iter = _track->begin_cluster_keys();
0134        iter != _track->end_cluster_keys();
0135        ++iter)
0136     {
0137       TrkrDefs::cluskey cluster_key = *iter;
0138       unsigned int layer = TrkrDefs::getLayer(cluster_key);
0139       unsigned int zelement = TrkrDefs::getZElement(cluster_key);
0140       unsigned int phielement = TrkrDefs::getPhiElement(cluster_key);
0141
0142       TrkrCluster *cluster =  _cluster_map->findCluster(cluster_key);
0143
0144       if(Verbosity() > 2)
0145         {
0146           std::cout << "   cluster " << cluster_key << " layer " << layer << " zelement " << zelement  << " phielement " << phielement << std::endl;
0147
0148           if(cluster)
0149         {
0150           double radius = sqrt( cluster->getX()*cluster->getX() + cluster->getY()*cluster->getY() );
0151           std::cout << "        cluster radius " << radius << " cluster adc " << cluster->getAdc() << std::endl;
0152           if(layer > 6)
0153             {
0154               cluster_avge_adc += (double) cluster->getAdc();
0155               cluster_avge_wt += 1.0;
0156             }
0157         }
0158           else
0159         std::cout << "Failed to find cluster with key :" << cluster_key << std::endl;
0160         }
0161
0162     }
0163       if(Verbosity() > 2)  std::cout << "   TPC cluster_avge_adc " << cluster_avge_adc/cluster_avge_wt << std::endl;
0164     }
0165
0166
0167   std::cout << "---- tracks with > 20 clusters " << ntracks << " mean nclusters " << mean_clusters /(double) ntracks << " mean pT " << mean_pt / (double) ntracks << std::endl;
0168
0169
0170   if (Verbosity() > 0)
0171     cout << "TPCGemGainCalb::process_event(PHCompositeNode *topNode) Leaving process_event" << endl;
0172
0173   return Fun4AllReturnCodes::EVENT_OK;
0174 }
0175
0176 //____________________________________________________________________________..
0177 int TPCGemGainCalb::EndRun(PHCompositeNode *topNode)
0178 {
0179
0180   return Fun4AllReturnCodes::EVENT_OK;
0181 }
0182
0183 //____________________________________________________________________________..
0184 int TPCGemGainCalb::End(PHCompositeNode *topNode)
0185 {
0186   fout->Write();
0187   //ntp->Write();
0188   //fout->Close();
0189
0190   return Fun4AllReturnCodes::EVENT_OK;
0191 }
0192
0193 //____________________________________________________________________________..
0194
0195 int  TPCGemGainCalb::GetNodes(PHCompositeNode* topNode)
0196 {
0197   //---------------------------------
0198   // Get Objects off of the Node Tree
0199   //---------------------------------
0200
0201   _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0202   if (!_cluster_map)
0203   {
0204     cerr << PHWHERE << " ERROR: Can't find node TRKR_CLUSTER" << endl;
0205     return Fun4AllReturnCodes::ABORTEVENT;
0206   }
0207
0208
0209   _vertex_map = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0210   if (!_vertex_map)
0211     {
0212       cerr << PHWHERE << " ERROR: Can't find SvtxVertexMap." << endl;
0213       return Fun4AllReturnCodes::ABORTEVENT;
0214     }
0215
0216
0217   _track_map = findNode::getClass<SvtxTrackMap>(topNode,  "SvtxTrackMap");
0218   if (!_track_map)
0219   {
0220     cerr << PHWHERE << " ERROR: Can't find SvtxTrackMap: " << endl;
0221     return Fun4AllReturnCodes::ABORTEVENT;
0222   }
0223
0224   return Fun4AllReturnCodes::EVENT_OK;
0225 }
0226