event_display/src/InttAna_MC_emc.cc

0001 #include "InttAna.h"
0002 #include <math.h>
0003 /// Cluster/Calorimeter includes
0004
0005 /// Fun4All includes
0006 #include <fun4all/Fun4AllHistoManager.h>
0007 #include <fun4all/Fun4AllReturnCodes.h>
0008 #include <phool/PHCompositeNode.h>
0009 #include <phool/getClass.h>
0010
0011 // #include <trackbase/TrkrClusterv4.h>
0012 // #include <trackbase/TrkrClusterContainerv3.h>
0013 #include <trackbase/TrkrCluster.h>
0014 #include <trackbase/TrkrClusterContainer.h>
0015 #include <trackbase/TrkrHit.h>
0016 #include <trackbase/TrkrHitSet.h>
0017 #include <trackbase/TrkrHitSetContainer.h>
0018 #include <trackbase/ActsGeometry.h>
0019 #include <trackbase/InttDefs.h>
0020
0021 #include <trackbase/InttEventInfo.h>
0022
0023 #include <globalvertex/GlobalVertex.h>
0024 #include <globalvertex/GlobalVertexMap.h>
0025
0026 // #include <globalvertex/MbdVertex.h>
0027 // #include <globalvertex/MbdVertexMap.h>
0028 #include <mbd/MbdOut.h>
0029 #include <ffarawobjects/Gl1RawHit.h>
0030
0031 #include <ffarawobjects/InttRawHit.h>
0032 #include <ffarawobjects/InttRawHitContainer.h>
0033
0034 #include "inttxyvertexfinder/InttVertex3D.h"
0035 #include "inttxyvertexfinder/InttVertex3DMap.h"
0036
0037 #include <globalvertex/SvtxVertex.h>
0038 #include <globalvertex/SvtxVertexMap.h>
0039 #include <trackbase_historic/SvtxTrack.h>
0040 #include <trackbase_historic/SvtxTrackMap.h>
0041
0042 #include <calobase/RawCluster.h>
0043 #include <calobase/RawClusterContainer.h>
0044
0045 /// ROOT includes
0046 #include <TFile.h>
0047 #include <TH1.h>
0048 #include <TH2.h>
0049 #include <TNtuple.h>
0050 #include <TTree.h>
0051
0052 /// C++ includes
0053 #include <cassert>
0054 #include <cmath>
0055 #include <sstream>
0056 #include <string>
0057
0058 #include "InttEvent.h"
0059 #include "InttRawData.h"
0060 // truth
0061 //  To get vertex
0062 #include <globalvertex/GlobalVertexMap.h>
0063 #include <globalvertex/GlobalVertex.h>
0064
0065 // TODO:
0066 /// HEPMC truth includes
0067 #pragma GCC diagnostic push
0068 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
0069 #include <HepMC/GenEvent.h>
0070 #include <HepMC/GenVertex.h>
0071 #pragma GCC diagnostic pop
0072
0073 #include <phhepmc/PHHepMCGenEvent.h>
0074 #include <phhepmc/PHHepMCGenEventMap.h>
0075
0076 #include <g4main/PHG4InEvent.h>
0077 #include <g4main/PHG4VtxPoint.h> // for PHG4Particle
0078 #include <g4main/PHG4Particle.h> // for PHG4Particle
0079 // truth end
0080
0081 using namespace std;
0082
0083 /**
0084  * This class demonstrates the construction and use of an analysis module
0085  * within the sPHENIX Fun4All framework. It is intended to show how to
0086  * obtain physics objects from the analysis tree, and then write them out
0087  * to a ROOT tree and file for further offline analysis.
0088  */
0089
0090 /**
0091  * Constructor of module
0092  */
0093 InttAna::InttAna(const std::string &name, const std::string &filename)
0094     : SubsysReco(name), _rawModule(nullptr), fname_(filename), anafile_(nullptr), h_zvtxseed_(nullptr)
0095 {
0096   xbeam_ = ybeam_ = 0.0;
0097 }
0098
0099 /**
0100  * Destructor of module
0101  */
0102 InttAna::~InttAna()
0103 {
0104 }
0105
0106 /**
0107  * Initialize the module and prepare looping over events
0108  */
0109 int InttAna::Init(PHCompositeNode * /*topNode*/)
0110 {
0111   if (Verbosity() > 5)
0112   {
0113     std::cout << "Beginning Init in InttAna" << std::endl;
0114   }
0115
0116   anafile_ = new TFile(fname_.c_str(), "recreate");
0117   h_dca2d_zero = new TH1F("h_dca2d_zero", "DCA2D to 0", 500, -10, 10);
0118   h2_dca2d_zero = new TH2F("h2_dca2d_zero", "DCA2D to 0 vs nclus", 500, -10, 10, 100, 0, 10000);
0119   h2_dca2d_len = new TH2F("h2_dca2d_len", "DCA2D to 0 vs len", 500, -10, 10, 500, -10, 10);
0120   h_zvtx = new TH1F("h_zvtx", "Zvertex", 400, -40, 40);
0121   h_ntp_clus = new TNtuple("ntp_clus", "Cluster Ntuple", "nclus:nclus2:bco_full:evt:size:adc:x:y:z:lay:lad:sen:lx:ly:phisize:zsize:zv;x_vtxsim;y_vtxsim;z_vtxsim");
0122   h_ntp_emcclus = new TNtuple("ntp_emcclus", "EMC Cluster Ntuple", "nemc:evt_emc:x_emc:y_emc:z_emc:e:ecore:size_emc");
0123
0124   h_ntp_cluspair = new TNtuple("ntp_cluspair", "Cluster Pair Ntuple", "nclus:nclus2:bco_full:evt:ang1:ang2:z1:z2:dca2d:len:unorm:l1:l2:vx:vy:vz:zvtx");
0125   h_ntp_emccluspair = new TNtuple("ntp_emccluspair", "INTT and EMC Cluster Pair Ntuple", "nclus:evt:ang1:ang2:z1:z2:dca2d:len:vx:vy:vz:eang:ez:ecore:esize:the1:the2:ethe");
0126   h_ntp_evt = new TNtuple("ntp_evt", "Event Ntuple", "nclus:nclus2:bco_full:evt:zv:zvs:zvm:zvc:bz:bqn:bqs:bfemclk:xvtx:yvtx:zvtx:nclusin:nclusout:nclszv:ntrkzv:chi2ndfzv:widthzv:ngroupzv:goodzv:peakratiozv:xvsim:yvsim:zvsim:xvsvtx:yvsvtx:zvsvtx:nmvtx0:nmvtx1:nmvtx2:ntrksvtx:nemc:nemc1");
0127
0128   h_zvtxseed_ = new TH1F("h_zvtxseed", "Zvertex Seed histogram", 200, -50, 50);
0129
0130   h_t_evt_bco = new TTree("t_evt_bco", "Event Tree for BCO");
0131   h_t_evt_bco->Branch("evt_gl1", &(m_evtbcoinfo.evt_gl1), "evt_gl1/I");
0132   h_t_evt_bco->Branch("evt_intt", &(m_evtbcoinfo.evt_intt), "evt_intt/i");
0133   h_t_evt_bco->Branch("evt_mbd", &(m_evtbcoinfo.evt_mbd), "evt_mbd/i");
0134   h_t_evt_bco->Branch("bco_gl1", &(m_evtbcoinfo.bco_gl1), "bco_gl1/l");
0135   h_t_evt_bco->Branch("bco_intt", &(m_evtbcoinfo.bco_intt), "bco_intt/l");
0136   h_t_evt_bco->Branch("bco_mbd", &(m_evtbcoinfo.bco_mbd), "bco_mbd/l");
0137   h_t_evt_bco->Branch("pevt_gl1", &(m_evtbcoinfo_prev.evt_gl1), "pevt_gl1/I");
0138   h_t_evt_bco->Branch("pevt_intt", &(m_evtbcoinfo_prev.evt_intt), "pevt_intt/i");
0139   h_t_evt_bco->Branch("pevt_mbd", &(m_evtbcoinfo_prev.evt_mbd), "pevt_mbd/i");
0140   h_t_evt_bco->Branch("pbco_gl1", &(m_evtbcoinfo_prev.bco_gl1), "pbco_gl1/l");
0141   h_t_evt_bco->Branch("pbco_intt", &(m_evtbcoinfo_prev.bco_intt), "pbco_intt/l");
0142   h_t_evt_bco->Branch("pbco_mbd", &(m_evtbcoinfo_prev.bco_mbd), "pbco_mbd/l");
0143
0144   m_hepmctree = new TTree("hepmctree", "A tree with hepmc truth particles");
0145   m_hepmctree->Branch("m_evt", &m_evt, "m_evt/I");
0146   m_hepmctree->Branch("m_xvtx", &m_xvtx, "m_xvtx/D");
0147   m_hepmctree->Branch("m_yvtx", &m_yvtx, "m_yvtx/D");
0148   m_hepmctree->Branch("m_zvtx", &m_zvtx, "m_zvtx/D");
0149   m_hepmctree->Branch("m_partid1", &m_partid1, "m_partid1/I");
0150   m_hepmctree->Branch("m_partid2", &m_partid2, "m_partid2/I");
0151   m_hepmctree->Branch("m_x1", &m_x1, "m_x1/D");
0152   m_hepmctree->Branch("m_x2", &m_x2, "m_x2/D");
0153   m_hepmctree->Branch("m_mpi", &m_mpi, "m_mpi/I");
0154   m_hepmctree->Branch("m_process_id", &m_process_id, "m_process_id/I");
0155   m_hepmctree->Branch("m_truthenergy", &m_truthenergy, "m_truthenergy/D");
0156   m_hepmctree->Branch("m_trutheta", &m_trutheta, "m_trutheta/D");
0157   m_hepmctree->Branch("m_truththeta", &m_truththeta, "m_truththeta/D");
0158   m_hepmctree->Branch("m_truthphi", &m_truthphi, "m_truthphi/D");
0159   m_hepmctree->Branch("m_status", &m_status, "m_status/I");
0160   m_hepmctree->Branch("m_truthpx", &m_truthpx, "m_truthpx/D");
0161   m_hepmctree->Branch("m_truthpy", &m_truthpy, "m_truthpy/D");
0162   m_hepmctree->Branch("m_truthpz", &m_truthpz, "m_truthpz/D");
0163   m_hepmctree->Branch("m_truthpt", &m_truthpt, "m_truthpt/D");
0164   m_hepmctree->Branch("m_numparticlesinevent", &m_numparticlesinevent, "m_numparticlesinevent/I");
0165   m_hepmctree->Branch("m_truthpid", &m_truthpid, "m_truthpid/I");
0166
0167   return 0;
0168 }
0169
0170 int InttAna::InitRun(PHCompositeNode * /*topNode*/)
0171 {
0172   cout << "InttAna::InitRun beamcenter " << xbeam_ << " " << ybeam_ << endl;
0173
0174   return 0;
0175 }
0176
0177 /**
0178  * Main workhorse function where each event is looped over and
0179  * data from each event is collected from the node tree for analysis
0180  */
0181 int InttAna::process_event(PHCompositeNode *topNode)
0182 {
0183   static int ievt = 0;
0184   cout << "InttEvt::process evt : " << ievt++ << endl;
0185
0186   if (Verbosity() > 5)
0187   {
0188     std::cout << "Beginning process_event in AnaTutorial" << std::endl;
0189   }
0190
0191   // readRawHit(topNode);
0192
0193   ActsGeometry *m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0194   if (!m_tGeometry)
0195   {
0196     std::cout << PHWHERE << "No ActsGeometry on node tree. Bailing." << std::endl;
0197     return Fun4AllReturnCodes::ABORTEVENT;
0198   }
0199   // TODO:
0200   PHHepMCGenEventMap *hepmceventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0201   if (!hepmceventmap)
0202   {
0203     cout << PHWHERE << "hepmceventmap node is missing." << endl;
0204     // return Fun4AllReturnCodes::ABORTEVENT;
0205   }
0206
0207   // TODO:
0208   PHG4InEvent *phg4inevent = findNode::getClass<PHG4InEvent>(topNode, "PHG4INEVENT");
0209   if (!phg4inevent)
0210   {
0211     cout << PHWHERE << "PHG4INEVENT node is missing." << endl;
0212     // return Fun4AllReturnCodes::ABORTEVENT;
0213   }
0214
0215   TrkrClusterContainer *m_clusterMap =
0216       findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0217   if (!m_clusterMap)
0218   {
0219     cout << PHWHERE << "TrkrClusterContainer node is missing." << endl;
0220     return Fun4AllReturnCodes::ABORTEVENT;
0221   }
0222
0223   GlobalVertexMap *vertexs =
0224       findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0225   if (!vertexs)
0226   {
0227     cout << PHWHERE << "GlobalVertexMap node is missing." << endl;
0228     // return Fun4AllReturnCodes::ABORTEVENT;
0229   }
0230
0231   InttEventInfo *inttevthead = findNode::getClass<InttEventInfo>(topNode, "INTTEVENTHEADER");
0232
0233   InttVertex3DMap *inttvertexmap = findNode::getClass<InttVertex3DMap>(topNode, "InttVertexMap");
0234
0235   double vtx_sim[3]{-9999, -9999, -9999};
0236
0237   std::map<GlobalVertex::VTXTYPE, std::string> sourcemap{
0238       {GlobalVertex::VTXTYPE::UNDEFINED, "UNDEFINED"},
0239       {GlobalVertex::VTXTYPE::TRUTH, "TRUTH"},
0240       {GlobalVertex::VTXTYPE::SMEARED, "SMEARED"},
0241       {GlobalVertex::VTXTYPE::MBD, "MBD"},
0242       {GlobalVertex::VTXTYPE::SVTX, "SVTX"},
0243       {GlobalVertex::VTXTYPE::SVTX_MBD, "SVTX_MBD"}};
0244   //   enum VTXTYPE
0245   std::string s_vtxid = "NULL";
0246
0247   GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0248   if (vertexmap)
0249   {
0250     if (!vertexmap->empty())
0251     {
0252
0253       for (auto vertex : *vertexmap)
0254       {
0255         std::cout << "map entry" << std::endl;
0256         std::map<GlobalVertex::VTXTYPE, std::string>::iterator itr_src;
0257         itr_src = sourcemap.find((GlobalVertex::VTXTYPE)vertex.first);
0258         if (itr_src != sourcemap.end())
0259         {
0260           s_vtxid = itr_src->second;
0261         }
0262
0263         if (vertex.first == GlobalVertex::VTXTYPE::TRUTH)
0264         {
0265           vtx_sim[0] = vertex.second->get_x();
0266           vtx_sim[1] = vertex.second->get_y();
0267           vtx_sim[2] = vertex.second->get_z();
0268         }
0269
0270         std::cout << std::endl
0271                   << "A " << vertex.second->get_x()
0272                   << " " << vertex.second->get_y()
0273                   << " " << vertex.second->get_z()
0274                   << " " << vertex.second->get_id()
0275                   << " " << s_vtxid << std::endl;
0276       }
0277     }
0278   }
0279
0280   /*
0281     MbdVertexMap *mbdvertexmap = findNode::getClass<MbdVertexMap>(topNode, "MbdVertexMap");
0282     MbdVertex* mbdvtx = nullptr;
0283     if (mbdvertexmap){
0284       cout<<"MbdVertexMap"<<endl;
0285       if (!mbdvertexmap->empty()){
0286
0287         for(auto mbdvertex : *mbdvertexmap)
0288         {
0289           std::cout<<"mbdmap entry"<<std::endl;
0290
0291           std::cout<<std::endl<<"Mbdz"
0292                               <<" " <<mbdvertex.second->get_z()
0293                               <<" " <<mbdvertex.second->get_id()<<endl;
0294           mbdvtx = mbdvertex;
0295         }
0296       }
0297       else {
0298         cout<<"MbdVertexMap empty"<<endl;
0299       }
0300     }
0301   */
0302
0303   SvtxVertexMap *svtxvertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0304   SvtxVertex *svtxvertex = nullptr;
0305   if (svtxvertexmap)
0306   {
0307     cout << "svtxvertex: size : " << svtxvertexmap->size() << endl;
0308     for (SvtxVertexMap::ConstIter iter = svtxvertexmap->begin(); iter != svtxvertexmap->end(); ++iter)
0309     {
0310       svtxvertex = iter->second;
0311       std::cout << std::endl
0312                 << "SvtxVertex"
0313                 << " " << svtxvertex->get_x()
0314                 << " " << svtxvertex->get_y()
0315                 << " " << svtxvertex->get_z()
0316                 << " " << svtxvertex->get_id() << endl;
0317     }
0318   }
0319   else
0320   {
0321     cout << "SvtxVertexMap is not available" << endl;
0322   }
0323
0324   MbdOut *mbdout = findNode::getClass<MbdOut>(topNode, "MbdOut");
0325   if (!mbdout)
0326   {
0327     cout << "No MbdOut" << endl;
0328   }
0329   double mbdqs = (mbdout != nullptr) ? mbdout->get_q(0) : -9999;
0330   double mbdqn = (mbdout != nullptr) ? mbdout->get_q(1) : -9999;
0331   //  double mbdt0 = (mbdout!=nullptr) ?  mbdout->get_t0() : -9999;
0332   double mbdz = (mbdout != nullptr) ? mbdout->get_zvtx() : -9999;
0333
0334   Gl1RawHit *gl1raw = findNode::getClass<Gl1RawHit>(topNode, "GL1RAWHIT");
0335   if (!gl1raw)
0336   {
0337     cout << "No Gl1Raw" << endl;
0338   }
0339
0340   InttRawHitContainer *inttrawmap = findNode::getClass<InttRawHitContainer>(topNode, "INTTRAWHIT");
0341   InttRawHit *inttraw = (inttrawmap != nullptr && inttrawmap->get_nhits() > 0) ? inttrawmap->get_hit(0) : nullptr;
0342   /////////////
0343   //  InttEvent from RawModule (not standard way)
0344   InttEvent *inttEvt = nullptr;
0345   if (_rawModule)
0346   {
0347     inttEvt = _rawModule->getInttEvent();
0348   }
0349   uint64_t bco = 0;
0350   int evtseq = -1;
0351   if (inttEvt != NULL)
0352   {
0353     bco = inttEvt->bco;
0354     evtseq = inttEvt->evtSeq;
0355   }
0356
0357   if (hepmceventmap != nullptr || phg4inevent != nullptr)
0358   {
0359     xvtx_sim = vertexs->find(GlobalVertex::TRUTH)->second->get_x();
0360     yvtx_sim = vertexs->find(GlobalVertex::TRUTH)->second->get_y();
0361     zvtx_sim = vertexs->find(GlobalVertex::TRUTH)->second->get_z();
0362
0363     // TODO:
0364     if (hepmceventmap != nullptr)
0365       getHEPMCTruth(topNode);
0366     else if (phg4inevent != nullptr)
0367       getPHG4Particle(topNode);
0368   }
0369
0370   if (inttevthead)
0371   {
0372     bco = inttevthead->get_bco_full();
0373     // evtseq = inttEvt->evtSeq;
0374   }
0375   cout << "bco 0x" << hex << bco << dec << endl;
0376
0377   if (gl1raw)
0378   {
0379     bco = gl1raw->get_bco();
0380     evtseq = gl1raw->getEvtSequence();
0381   }
0382
0383   double vertex[10][3]{};
0384
0385   InttVertex3D *zvtxobj = NULL;
0386   if (inttvertexmap)
0387   {
0388     int ivtx = 0;
0389     cout << "vertex map size : " << inttvertexmap->size() << endl;
0390     InttVertex3DMap::ConstIter biter_beg = inttvertexmap->begin();
0391     InttVertex3DMap::ConstIter biter_end = inttvertexmap->end();
0392     // cout<<"vertex map size : after size "<<endl;
0393     // inttvertexmap->identify();
0394
0395     for (InttVertex3DMap::ConstIter biter = biter_beg; biter != biter_end; ++biter)
0396     {
0397       InttVertex3D *vtxobj = biter->second;
0398       if (vtxobj)
0399       {
0400         cout << "vtxobj" << ivtx << endl;
0401         vtxobj->identify();
0402         vertex[ivtx][0] = vtxobj->get_x();
0403         vertex[ivtx][1] = vtxobj->get_y();
0404         vertex[ivtx][2] = vtxobj->get_z();
0405         if (ivtx == 2)
0406         {
0407           zvtxobj = vtxobj;
0408         }
0409       }
0410       else
0411       {
0412         cout << "no vtx object : " << ivtx << endl;
0413       }
0414       ivtx++;
0415       if (ivtx >= 10)
0416       {
0417         cout << "n_vertex exceeded : " << ivtx << endl;
0418         break;
0419       }
0420     }
0421   }
0422
0423   ////
0424   // fill evtbco_info
0425   m_evtbcoinfo.clear();
0426   m_evtbcoinfo.evt_gl1 = (gl1raw != nullptr) ? gl1raw->getEvtSequence() : -1;
0427   m_evtbcoinfo.evt_intt = (inttraw != nullptr) ? inttraw->get_event_counter() : -1;
0428   m_evtbcoinfo.evt_mbd = (mbdout != nullptr) ? mbdout->get_evt() : -1;
0429   m_evtbcoinfo.bco_gl1 = (gl1raw != nullptr) ? gl1raw->get_bco() : 0;
0430   m_evtbcoinfo.bco_intt = (inttevthead != nullptr) ? inttevthead->get_bco_full() : 0;
0431   m_evtbcoinfo.bco_mbd = (mbdout != nullptr) ? mbdout->get_femclock() : 0;
0432
0433   cout << "event : " << m_evtbcoinfo.evt_gl1 << " "
0434        << m_evtbcoinfo.evt_intt << " "
0435        << m_evtbcoinfo.evt_mbd << " :  "
0436        << hex << m_evtbcoinfo.bco_gl1 << dec << " "
0437        << hex << m_evtbcoinfo.bco_intt << dec << " "
0438        << hex << m_evtbcoinfo.bco_mbd << dec << " : "
0439        << hex << m_evtbcoinfo.bco_gl1 - m_evtbcoinfo_prev.bco_gl1 << dec << " "
0440        << hex << m_evtbcoinfo.bco_intt - m_evtbcoinfo_prev.bco_intt << dec << " "
0441        << hex << m_evtbcoinfo.bco_mbd - m_evtbcoinfo_prev.bco_mbd << dec << endl;
0442
0443   /////////////
0444   int nclusadd = 0, nclusadd2 = 0;
0445   int nclus_inner = 0, nclus_outer = 0;
0446   for (unsigned int inttlayer = 0; inttlayer < 4; inttlayer++)
0447   {
0448     for (const auto &hitsetkey : m_clusterMap->getHitSetKeys(TrkrDefs::TrkrId::inttId, inttlayer + 3))
0449     {
0450       auto range = m_clusterMap->getClusters(hitsetkey);
0451
0452       for (auto clusIter = range.first; clusIter != range.second; ++clusIter)
0453       {
0454         nclusadd++;
0455         const auto cluster = clusIter->second;
0456         if (cluster->getAdc() > 40)
0457           nclusadd2++;
0458
0459         if (inttlayer < 2)
0460           nclus_inner++;
0461         else
0462           nclus_outer++;
0463       }
0464     }
0465   }
0466   /////////////
0467   // MVTX cluster
0468   std::set<TrkrDefs::TrkrId> detectors;
0469   detectors.insert(TrkrDefs::TrkrId::mvtxId);
0470   int nclusmvtx[3] = {0, 0, 0};
0471   for (const auto &det : detectors)
0472   {
0473     for (const auto &layer : {0, 1, 2})
0474     {
0475       for (const auto &hitsetkey : m_clusterMap->getHitSetKeys(det, layer))
0476       {
0477         auto range = m_clusterMap->getClusters(hitsetkey);
0478         for (auto citer = range.first; citer != range.second; ++citer)
0479         {
0480           // const auto ckey = citer->first;
0481           // const auto cluster = citer->second;
0482           // const auto global = m_tGeometry->getGlobalPosition(ckey, cluster);
0483           nclusmvtx[layer]++;
0484         }
0485       }
0486     }
0487   }
0488
0489   /////////////
0490   // SvtxTrack
0491   SvtxTrackMap *svtxtrackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0492
0493   /////////////
0494   // EMC Cluster
0495   RawClusterContainer *clusterContainer = findNode::getClass<RawClusterContainer>(topNode, "CLUSTERINFO_POS_COR_CEMC");
0496   // RawClusterContainer *clusterContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_POS_COR_CEMC");
0497   // RawClusterContainer *clusterContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_CEMC");
0498   int nemc = -9999, nemc1 = -9999;
0499   if (!clusterContainer)
0500   {
0501     std::cout << PHWHERE << "InttAna::process_event - Fatal Error - CLUSTER_POS_COR_CEMC node is missing. " << std::endl;
0502   }
0503   else
0504   {
0505     float ntpval_emccls[20];
0506     // This is how we iterate over items in the container.
0507     RawClusterContainer::ConstRange clusterEnd = clusterContainer->getClusters();
0508     RawClusterContainer::ConstIterator clusterIter;
0509
0510     nemc = clusterContainer->size();
0511     nemc1 = 0;
0512
0513     for (clusterIter = clusterEnd.first; clusterIter != clusterEnd.second; clusterIter++)
0514     {
0515       RawCluster *recoCluster = clusterIter->second;
0516
0517       ntpval_emccls[0] = nemc;
0518       ntpval_emccls[1] = evtseq;
0519       ntpval_emccls[2] = recoCluster->get_x();
0520       ntpval_emccls[3] = recoCluster->get_y();
0521       ntpval_emccls[4] = recoCluster->get_z();
0522       ntpval_emccls[5] = recoCluster->get_energy();
0523       ntpval_emccls[6] = recoCluster->get_ecore();
0524       ntpval_emccls[7] = recoCluster->getNTowers();
0525
0526       h_ntp_emcclus->Fill(ntpval_emccls);
0527
0528       if (recoCluster->get_ecore() > 0.15)
0529         nemc1++;
0530     }
0531   }
0532
0533   /////////////
0534   static int evtCount = 0;
0535
0536   cout << "evtCount : " << evtCount << " " << evtseq << " " << hex << bco << dec << endl;
0537
0538   struct ClustInfo
0539   {
0540     int layer;
0541     Acts::Vector3 pos;
0542   };
0543
0544   float ntpval[20];
0545   int nCluster = 0;
0546   bool exceedNwrite = false;
0547   // std::vector<Acts::Vector3> clusters[2]; // inner=0; outer=1
0548   std::vector<ClustInfo> clusters[2]; // inner=0; outer=1
0549   for (unsigned int inttlayer = 0; inttlayer < 4; inttlayer++)
0550   {
0551     int layer = (inttlayer < 2 ? 0 : 1);
0552     for (const auto &hitsetkey : m_clusterMap->getHitSetKeys(TrkrDefs::TrkrId::inttId, inttlayer + 3))
0553     {
0554       auto range = m_clusterMap->getClusters(hitsetkey);
0555
0556       for (auto clusIter = range.first; clusIter != range.second; ++clusIter)
0557       {
0558         const auto cluskey = clusIter->first;
0559         const auto cluster = clusIter->second;
0560
0561         int ladder_z = InttDefs::getLadderZId(cluskey);
0562         int ladder_phi = InttDefs::getLadderPhiId(cluskey);
0563         int size = cluster->getSize();
0564
0565         // cout    <<cluster->getPosition(0)<<" "
0566         //         <<cluster->getPosition(1)<<" : "
0567         //         <<cluster->getAdc()<<" "<<size<<" "<<inttlayer<<" "<<ladder_z<<" "<<ladder_phi<<endl;
0568
0569         const auto globalPos = m_tGeometry->getGlobalPosition(cluskey, cluster);
0570
0571         if (nCluster < 5)
0572         {
0573           cout << "xyz : " << globalPos.x() << " " << globalPos.y() << " " << globalPos.z() << " :  "
0574                << cluster->getPosition(0) << " "
0575                << cluster->getPosition(1) << " : "
0576                << cluster->getAdc() << " " << size << " " << inttlayer << " " << ladder_z << " " << ladder_phi << endl;
0577         }
0578         else
0579         {
0580           if (!exceedNwrite)
0581           {
0582             cout << " exceed : ncluster limit.  no more cluster xyz printed" << endl;
0583             exceedNwrite = true;
0584           }
0585         }
0586
0587         ClustInfo info;
0588         info.layer = inttlayer;
0589         info.pos = globalPos;
0590
0591         // clusters[layer].push_back(globalPos);
0592         clusters[layer].push_back(info);
0593         nCluster++;
0594
0595         ntpval[0] = nclusadd;          // bco_full
0596         ntpval[1] = nclusadd2;         // bco_full
0597         ntpval[2] = bco;               // bco_full
0598         ntpval[3] = evtseq;            // evtCount;
0599         ntpval[4] = size;              // size
0600         ntpval[5] = cluster->getAdc(); // size
0601         ntpval[6] = globalPos.x();
0602         ntpval[7] = globalPos.y();
0603         ntpval[8] = globalPos.z();
0604         ntpval[9] = inttlayer;
0605         ntpval[10] = ladder_phi;
0606         ntpval[11] = ladder_z;
0607         ntpval[12] = cluster->getPosition(0);
0608         ntpval[13] = cluster->getPosition(1);
0609         ntpval[14] = cluster->getPhiSize();
0610         ntpval[15] = cluster->getZSize();
0611         ntpval[16] = vertex[2][2]; // zvtx;
0612         ntpval[17] = xvtx_sim;
0613         ntpval[18] = yvtx_sim;
0614         ntpval[19] = zvtx_sim;
0615         h_ntp_clus->Fill(ntpval);
0616         // = new TNtuple("ntp_clus", "Cluster Ntuple", "bco_full:evt:size:adc:x:y:z:lay:lad:sen");
0617       }
0618     }
0619   }
0620   cout << "nCluster : " << nCluster << endl;
0621
0622   double zvtx = -9999.;
0623
0624   struct cluspair
0625   {
0626     float p1_ang;
0627     float p2_ang;
0628     float p1_r;
0629     float p2_r;
0630     float p1_z;
0631     float p2_z;
0632     float dca_p0;
0633     float len_p0;
0634   };
0635
0636   vector<cluspair> vcluspair;
0637
0638   if (nCluster < 300)
0639   // if(nCluster<500)
0640   {
0641     float ntpval2[20];
0642     Acts::Vector3 beamspot = Acts::Vector3(xbeam_, ybeam_, 0);
0643     vector<double> vz_array;
0644     for (auto c1 = clusters[0].begin(); c1 != clusters[0].end(); ++c1) // inner
0645     {
0646       for (auto c2 = clusters[1].begin(); c2 != clusters[1].end(); ++c2) // outer
0647       {
0648         Acts::Vector3 p1 = c1->pos - beamspot;
0649         Acts::Vector3 p2 = c2->pos - beamspot;
0650         Acts::Vector3 u = p2 - p1;
0651         double unorm = sqrt(u.x() * u.x() + u.y() * u.y());
0652         if (unorm < 0.00001)
0653           continue;
0654
0655         // skip bad compbination
0656         double p1_ang = atan2(p1.y(), p1.x());
0657         double p2_ang = atan2(p2.y(), p2.x());
0658         double d_ang = p2_ang - p1_ang;
0659
0660         // unit vector in 2D
0661         double ux = u.x() / unorm;
0662         double uy = u.y() / unorm;
0663         double uz = u.z() / unorm; // same normalization factor(2D) is multiplied
0664
0665         // Acts::Vector3 p0   = beamspot - p1;
0666         Acts::Vector3 p0 = Acts::Vector3(0, 0, 0) - p1; // p1, p2 are already shifted by beam center
0667
0668         double dca_p0 = p0.x() * uy - p0.y() * ux; // cross product of p0 x u
0669         double len_p0 = p0.x() * ux + p0.y() * uy; // dot product of p0 . u
0670
0671         // beam center in X-Y plane
0672         double vx = len_p0 * ux + p1.x();
0673         double vy = len_p0 * uy + p1.y();
0674
0675         double vz = len_p0 * uz + p1.z();
0676
0677         double p1_r = sqrt(p1.y() * p1.y() + p1.x() * p1.x());
0678         double p2_r = sqrt(p2.y() * p2.y() + p2.x() * p2.x());
0679
0680         // if(unorm>4.5||d_ang<-0.005*3.1415||0.005*3.1415<d_ang)
0681         // if(unorm>4.5||d_ang<-0.25*3.1415||0.25*3.1415<d_ang)
0682         if (fabs(d_ang) < 0.05 && fabs(dca_p0) < 0.5)
0683         {
0684           h_zvtxseed_->Fill(vz);
0685           vz_array.push_back(vz);
0686
0687           cluspair clspair;
0688           clspair.p1_ang = p1_ang;
0689           clspair.p2_ang = p2_ang;
0690           clspair.p1_r = p1_r;
0691           clspair.p2_r = p2_r;
0692           clspair.p1_z = p1.z();
0693           clspair.p2_z = p2.z();
0694           clspair.dca_p0 = dca_p0;
0695           clspair.len_p0 = len_p0;
0696           vcluspair.push_back(clspair);
0697         }
0698
0699         h_dca2d_zero->Fill(dca_p0);
0700         h2_dca2d_zero->Fill(dca_p0, nCluster);
0701         h2_dca2d_len->Fill(dca_p0, len_p0);
0702         // cout<<"dca : "<<dca_p0<<endl;
0703         //
0704         ntpval2[0] = nclusadd;
0705         ntpval2[1] = nclusadd2;
0706         ntpval2[2] = 0;
0707         ntpval2[3] = evtCount;
0708         ntpval2[4] = p1_ang;
0709         ntpval2[5] = p2_ang;
0710         ntpval2[6] = p1.z();
0711         ntpval2[7] = p2.z();
0712         ntpval2[8] = dca_p0;
0713         ntpval2[9] = len_p0;
0714         ntpval2[10] = unorm;
0715         ntpval2[11] = c1->layer;
0716         ntpval2[12] = c2->layer;
0717         ntpval2[13] = vx;
0718         ntpval2[14] = vy;
0719         ntpval2[15] = vz;
0720         ntpval2[16] = vertex[2][2]; // zvtx;
0721
0722         if (evtCount < 10000)
0723           h_ntp_cluspair->Fill(ntpval2);
0724         // h_ntp_cluspair = new TNtuple("ntp_cluspair", "Cluster Pair Ntuple", "nclus:bco_full:evt:ang1:ang2:dca2d:len:unorm");
0725       }
0726     }
0727
0728     // calculate trancated mean of DCA~Z histogram as Z-vertex position
0729
0730     if (vz_array.size() > 3)
0731     {
0732       double zbin = h_zvtxseed_->GetMaximumBin();
0733       double zcenter = h_zvtxseed_->GetBinCenter(zbin);
0734       double zmean = h_zvtxseed_->GetMean();
0735       double zrms = h_zvtxseed_->GetRMS();
0736       if (zrms < 20)
0737         zrms = 20;
0738
0739       double zmax = zcenter + zrms; // 1 sigma
0740       double zmin = zcenter - zrms; // 1 sigma
0741
0742       double zsum = 0.;
0743       int zcount = 0;
0744       for (auto iz = vz_array.begin(); iz != vz_array.end(); ++iz)
0745       {
0746         double vz = (*iz);
0747         if (zmin < vz && vz < zmax)
0748         {
0749           zsum += vz;
0750           zcount++;
0751         }
0752       }
0753       if (zcount > 0)
0754         zvtx = zsum / zcount;
0755
0756       cout << "ZVTX: " << zvtx << " " << zcenter << " " << zmean << " " << zrms << " " << zbin << endl;
0757     }
0758
0759     h_zvtx->Fill(zvtx);
0760
0761     float ntpval_emcpair[20];
0762     if (clusterContainer != nullptr)
0763     {
0764       RawClusterContainer::ConstRange clusterEnd = clusterContainer->getClusters();
0765       RawClusterContainer::ConstIterator clusterIter;
0766
0767       /////////////
0768       //  association to EMCAL
0769       for (vector<cluspair>::iterator itrcp = vcluspair.begin(); itrcp != vcluspair.end(); ++itrcp)
0770       {
0771         cluspair &clspair = *itrcp;
0772
0773         for (clusterIter = clusterEnd.first; clusterIter != clusterEnd.second; clusterIter++)
0774         {
0775           RawCluster *recoCluster = clusterIter->second;
0776
0777           float the1 = atan2(clspair.p1_r, (clspair.p1_z - vertex[2][2]));
0778           float the2 = atan2(clspair.p2_r, (clspair.p2_z - vertex[2][2]));
0779
0780           Acts::Vector3 emcpos = Acts::Vector3(recoCluster->get_x(), recoCluster->get_y(), recoCluster->get_z());
0781           Acts::Vector3 emcposc = emcpos - beamspot;
0782           float ephi = atan2(emcposc.y(), emcposc.x());
0783           float er = sqrt(emcposc.x() * emcposc.x() + emcposc.y() * emcposc.y());
0784           float ethe = atan2(er, emcposc.z() - vertex[2][2]);
0785
0786           ntpval_emcpair[0] = nclusadd;
0787           ntpval_emcpair[1] = evtseq;
0788           ntpval_emcpair[2] = clspair.p1_ang; // INTT
0789           ntpval_emcpair[3] = clspair.p2_ang;
0790           ntpval_emcpair[4] = clspair.p1_z;
0791           ntpval_emcpair[5] = clspair.p2_z;
0792           ntpval_emcpair[6] = clspair.dca_p0;
0793           ntpval_emcpair[7] = clspair.len_p0;
0794           ntpval_emcpair[8] = vertex[1][0];  // x-vertex
0795           ntpval_emcpair[9] = vertex[1][1];  // y-vertex
0796           ntpval_emcpair[10] = vertex[2][2]; // y-vertex
0797           ntpval_emcpair[11] = ephi;
0798           ntpval_emcpair[12] = recoCluster->get_z(); // EMCal
0799           ntpval_emcpair[13] = recoCluster->get_ecore();
0800           ntpval_emcpair[14] = recoCluster->getNTowers();
0801           ntpval_emcpair[15] = the1;
0802           ntpval_emcpair[16] = the2;
0803           ntpval_emcpair[17] = ethe;
0804
0805           h_ntp_emccluspair->Fill(ntpval_emcpair);
0806         }
0807       }
0808     }
0809   }
0810
0811   float ntpval3[40];
0812   ntpval3[0] = nclusadd;  // bco_full
0813   ntpval3[1] = nclusadd2; // bco_full
0814   ntpval3[2] = bco;
0815   ntpval3[3] = evtseq;
0816   ntpval3[4] = vertex[2][2]; // zvtx;
0817   ntpval3[5] = zvtx;         // zrms;
0818   ntpval3[6] = 0;            // zmean;
0819   ntpval3[7] = 0;            // zcenter;
0820   ntpval3[8] = mbdz;
0821   ntpval3[9] = mbdqn;         // mbdt.bqn;
0822   ntpval3[10] = mbdqs;        // mbdt.bqs;
0823   ntpval3[11] = 0;            // mbdt.femclk;
0824   ntpval3[12] = vertex[1][0]; // x-vertex
0825   ntpval3[13] = vertex[1][1]; // y-vertex
0826   ntpval3[14] = vertex[2][2]; // y-vertex
0827   ntpval3[15] = nclus_inner;  // x-vertex
0828   ntpval3[16] = nclus_outer;  // y-vertex
0829
0830   ntpval3[17] = zvtxobj->get_nclus();
0831   ntpval3[18] = zvtxobj->get_ntracklet();
0832   ntpval3[19] = zvtxobj->get_chi2ndf();
0833   ntpval3[20] = zvtxobj->get_width();
0834   ntpval3[21] = zvtxobj->get_ngroup();
0835   ntpval3[22] = (zvtxobj->get_good()) ? 1 : 0;
0836   ntpval3[23] = zvtxobj->get_peakratio();
0837
0838   ntpval3[24] = vtx_sim[0]; // sim x-vertex
0839   ntpval3[25] = vtx_sim[1]; // sim y-vertex
0840   ntpval3[26] = vtx_sim[2]; // sim z-vertex
0841
0842   ntpval3[27] = (svtxvertex != nullptr) ? svtxvertex->get_x() : -9999; // svtx vertex x
0843   ntpval3[28] = (svtxvertex != nullptr) ? svtxvertex->get_y() : -9999; // svtx vertex y
0844   ntpval3[29] = (svtxvertex != nullptr) ? svtxvertex->get_z() : -9999; // svtx vertex z
0845
0846   ntpval3[30] = nclusmvtx[0];
0847   ntpval3[31] = nclusmvtx[1];
0848   ntpval3[32] = nclusmvtx[2];
0849   ntpval3[33] = (svtxtrackmap != nullptr) ? svtxtrackmap->size() : -9999;
0850   ntpval3[34] = nemc;
0851   ntpval3[35] = nemc1;
0852
0853   h_ntp_evt->Fill(ntpval3);
0854   // h_ntp_evt = new TNtuple("ntp_evt", "Event Ntuple",
0855   //"nclus:nclus2:bco_full:evt:zv:
0856   // zvs:zvm:zvc:bz:bqn:
0857   // bqs:bfemclk:xvtx:yvtx:zvtx:
0858   // nclusin:nclusout:nclszv:ntrkzv:chi2ndfzv:
0859   // widthzv:ngroupzv:goodzv:peakratiozv:xvsim:
0860   // yvsim:zvsim");
0861   //
0862
0863   h_t_evt_bco->Fill();
0864
0865   m_evtbcoinfo_prev.copy(m_evtbcoinfo);
0866
0867   evtCount++;
0868
0869   return Fun4AllReturnCodes::EVENT_OK;
0870 }
0871
0872 /**
0873  * End the module and finish any data collection. Clean up any remaining
0874  * loose ends
0875  */
0876 int InttAna::End(PHCompositeNode * /*topNode*/)
0877 {
0878   if (Verbosity() > 1)
0879   {
0880     std::cout << "Ending InttAna analysis package" << std::endl;
0881   }
0882
0883   m_hepmctree->Write();
0884
0885   if (anafile_ != nullptr)
0886   {
0887     anafile_->Write();
0888     anafile_->Close();
0889   }
0890
0891   return 0;
0892 }
0893
0894 void InttAna::readRawHit(PHCompositeNode *topNode)
0895 {
0896   TrkrHitSetContainer *m_hits = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0897   if (!m_hits)
0898   {
0899     cout << PHWHERE << "ERROR: Can't find node TRKR_HITSET" << endl;
0900     return;
0901   }
0902
0903   // uint8_t getLadderZId(TrkrDefs::hitsetkey key);
0904   // uint8_t getLadderPhiId(TrkrDefs::hitsetkey key);
0905   // int getTimeBucketId(TrkrDefs::hitsetkey key);
0906   //
0907   // uint16_t getCol(TrkrDefs::hitkey key); // z-strip = offline chip
0908   // uint16_t getRow(TrkrDefs::hitkey key); // y-strip = offline channel
0909
0910   // loop over the InttHitSet objects
0911   TrkrHitSetContainer::ConstRange hitsetrange =
0912       m_hits->getHitSets(TrkrDefs::TrkrId::inttId);
0913   for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first;
0914        hitsetitr != hitsetrange.second;
0915        ++hitsetitr)
0916   {
0917     // Each hitset contains only hits that are clusterizable - i.e. belong to a single sensor
0918     TrkrHitSet *hitset = hitsetitr->second;
0919
0920     if (Verbosity() > 1)
0921       cout << "InttClusterizer found hitsetkey " << hitsetitr->first << endl;
0922     if (Verbosity() > 2)
0923       hitset->identify();
0924
0925     // we have a single hitset, get the info that identifies the sensor
0926     int layer = TrkrDefs::getLayer(hitsetitr->first);
0927     int ladder_z_index = InttDefs::getLadderZId(hitsetitr->first);
0928     int ladder_phi_index = InttDefs::getLadderPhiId(hitsetitr->first);
0929     int ladder_bco_index = InttDefs::getTimeBucketId(hitsetitr->first);
0930
0931     cout << "layer " << layer << " " << ladder_z_index << " " << ladder_phi_index << " " << ladder_bco_index << endl;
0932
0933     //// we will need the geometry object for this layer to get the global position
0934     // CylinderGeomIntt* geom = dynamic_cast<CylinderGeomIntt*>(geom_container->GetLayerGeom(layer));
0935     // float pitch = geom->get_strip_y_spacing();
0936     // float length = geom->get_strip_z_spacing();
0937
0938     // fill a vector of hits to make things easier - gets every hit in the hitset
0939     std::vector<std::pair<TrkrDefs::hitkey, TrkrHit *>> hitvec;
0940     TrkrHitSet::ConstRange hitrangei = hitset->getHits();
0941     for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
0942          hitr != hitrangei.second;
0943          ++hitr)
0944     {
0945       hitvec.push_back(make_pair(hitr->first, hitr->second));
0946       int chip = InttDefs::getCol(hitr->first);
0947       int chan = InttDefs::getRow(hitr->first);
0948       int adc = (hitr->second)->getAdc();
0949       cout << "     hit : " << chip << " " << chan << " " << adc << endl;
0950     }
0951     cout << "hitvec.size(): " << hitvec.size() << endl;
0952   }
0953   void InttAna::getHEPMCTruth(PHCompositeNode *topNode)
0954   {
0955     cout << "getHEPMCTruth!!!" << endl;
0956     /// Get the node from the node tree
0957     PHHepMCGenEventMap *hepmceventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0958
0959     /// If the node was not properly put on the tree, return
0960     if (!hepmceventmap)
0961     {
0962       std::cout << PHWHERE
0963                 << "HEPMC event map node is missing, can't collected HEPMC truth particles"
0964                 << std::endl;
0965       return;
0966     }
0967
0968     /// Could have some print statements for debugging with verbosity
0969     if (Verbosity() > 1)
0970     {
0971       std::cout << "Getting HEPMC truth particles " << std::endl;
0972     }
0973
0974     /// You can iterate over the number of events in a hepmc event
0975     /// for pile up events where you have multiple hard scatterings per bunch crossing
0976     for (PHHepMCGenEventMap::ConstIter eventIter = hepmceventmap->begin();
0977          eventIter != hepmceventmap->end();
0978          ++eventIter)
0979     {
0980       /// Get the event
0981       PHHepMCGenEvent *hepmcevent = eventIter->second;
0982
0983       if (hepmcevent)
0984       {
0985         /// Get the event characteristics, inherited from HepMC classes
0986         HepMC::GenEvent *truthevent = hepmcevent->getEvent();
0987         if (!truthevent)
0988         {
0989           std::cout << PHWHERE
0990                     << "no evt pointer under phhepmvgeneventmap found "
0991                     << std::endl;
0992           return;
0993         }
0994
0995         /// Get the parton info
0996         HepMC::PdfInfo *pdfinfo = truthevent->pdf_info();
0997
0998         /// Get the parton info as determined from HEPMC
0999         m_partid1 = pdfinfo->id1();
1000         m_partid2 = pdfinfo->id2();
1001         m_x1 = pdfinfo->x1();
1002         m_x2 = pdfinfo->x2();
1003
1004         /// Are there multiple partonic intercations in a p+p event
1005         m_mpi = truthevent->mpi();
1006
1007         /// Get the PYTHIA signal process id identifying the 2-to-2 hard process
1008         m_process_id = truthevent->signal_process_id();
1009
1010         if (Verbosity() > 2)
1011         {
1012           std::cout << " Iterating over an event" << std::endl;
1013         }
1014
1015         /// Loop over all the truth particles and get their information
1016         for (HepMC::GenEvent::particle_const_iterator iter = truthevent->particles_begin();
1017              iter != truthevent->particles_end();
1018              ++iter)
1019         {
1020           /// Get each pythia particle characteristics
1021           m_truthenergy = (*iter)->momentum().e();
1022           m_truthpid = (*iter)->pdg_id();
1023
1024           // TODO: truth information
1025           m_xvtx = xvtx_sim;
1026           m_yvtx = yvtx_sim;
1027           m_zvtx = zvtx_sim;
1028           m_trutheta = (*iter)->momentum().pseudoRapidity();
1029           m_truththeta = 2 * atan(exp(-m_trutheta));
1030           // h_eta->Fill(m_trutheta);
1031           m_truthphi = (*iter)->momentum().phi();
1032           m_truthpx = (*iter)->momentum().px();
1033           m_truthpy = (*iter)->momentum().py();
1034           m_truthpz = (*iter)->momentum().pz();
1035
1036           m_truthpt = sqrt(m_truthpx * m_truthpx + m_truthpy * m_truthpy);
1037           m_status = (*iter)->status();
1038
1039           // m_truthphi = atan2(m_truthpy,m_truthpx);
1040
1041           // cout << "status: " << (*iter)->status() << " " << m_truthpid << endl;
1042
1043           h_phi->Fill(m_truthphi);
1044           h_theta->Fill(m_truththeta);
1045
1046           /*std::cout << "m_trutheta = " << m_trutheta << "  m_truthphi = " << m_truthphi << std::endl;*/
1047
1048           /// Fill the truth tree
1049           m_hepmctree->Fill();
1050           m_numparticlesinevent++;
1051         }
1052
1053         // for (HepMC::GenEventVertexRange::vertex_const_iterator viter = truthevent->vertices_begin();
1054         //      viter != truthevent->vertices_end();
1055         //      ++viter)
1056         // {
1057         //   /// Get each pythia particle characteristics
1058         //   m_vertex = (*viter)->vertex_range();
1059         //   cout<<m_vertex<<endl;
1060
1061         //   /// Fill the truth tree
1062         //   m_hepmctree->Fill();
1063         // }
1064       }
1065       // cout << "truth event = " << m_evt << endl;
1066       m_evt++;
1067     }
1068   }
1069
1070   void InttAna::getPHG4Particle(PHCompositeNode * topNode)
1071   {
1072     cout << "getPHG4Particle!!!" << endl;
1073     /// Get the node from the node tree
1074     PHG4InEvent *phg4inevent = findNode::getClass<PHG4InEvent>(topNode, "PHG4INEVENT");
1075
1076     /// If the node was not properly put on the tree, return
1077     if (!phg4inevent)
1078     {
1079       std::cout << PHWHERE
1080                 << "PHG4InEvent node is missing, can't collected PHG4 truth particles"
1081                 << std::endl;
1082       return;
1083     }
1084
1085     /// Could have some print statements for debugging with verbosity
1086     if (Verbosity() > 1)
1087     {
1088       std::cout << "Getting PHG4InEvent truth particles " << std::endl;
1089     }
1090
1091     /// You can iterate over the number of events in a hepmc event
1092     /// for pile up events where you have multiple hard scatterings per bunch crossing
1093     std::map<int, PHG4VtxPoint *>::const_iterator vtxiter;
1094     std::multimap<int, PHG4Particle *>::const_iterator particle_iter;
1095     std::pair<std::map<int, PHG4VtxPoint *>::const_iterator, std::map<int, PHG4VtxPoint *>::const_iterator> vtxbegin_end = phg4inevent->GetVertices();
1096
1097     for (vtxiter = vtxbegin_end.first; vtxiter != vtxbegin_end.second; ++vtxiter)
1098     {
1099       //*fout << "vtx number: " << vtxiter->first << std::endl;
1100       //(*vtxiter->second).identify(*fout);
1101       std::pair<std::multimap<int, PHG4Particle *>::const_iterator,
1102                 std::multimap<int, PHG4Particle *>::const_iterator>
1103           particlebegin_end = phg4inevent->GetParticles(vtxiter->first);
1104
1105       PHG4VtxPoint *vtx = vtxiter->second;
1106       m_xvtx = vtx->get_x();
1107       m_yvtx = vtx->get_y();
1108       m_zvtx = vtx->get_z();
1109       // virtual double get_t() const { return std::numeric_limits<double>::quiet_NaN(); }
1110
1111       for (particle_iter = particlebegin_end.first; particle_iter != particlebegin_end.second; ++particle_iter)
1112       {
1113         //      (particle_iter->second)->identify(*fout);
1114         PHG4Particle *part = particle_iter->second;
1115
1116         m_truthenergy = part->get_e();
1117         m_truthpid = part->get_pid();
1118
1119         m_truthpx = part->get_px();
1120         m_truthpy = part->get_py();
1121         m_truthpz = part->get_pz();
1122         m_truthphi = atan2(m_truthpy, m_truthpx);
1123         m_truthpt = sqrt(m_truthpx * m_truthpx + m_truthpy * m_truthpy);
1124
1125         m_truththeta = atan2(m_truthpt, m_truthpz);
1126         m_trutheta = -log(tan(0.5 * m_truththeta));
1127         m_status = 1;
1128
1129         // part->get_pid() const override { return fpid; }
1130         // part->get_name() const override { return fname; }
1131         // part->get_e() const override { return fe; }
1132         // part->get_track_id();
1133         // part->get_vtx_id() const override { return vtxid; }
1134         // part->get_parent_id() const override { return parentid; }
1135         // int get_primary_id() const override { return primaryid; }
1136
1137         /// Fill the truth tree
1138         h_phi->Fill(m_truthphi);
1139         h_theta->Fill(m_truththeta);
1140
1141         m_hepmctree->Fill();
1142         m_numparticlesinevent++;
1143       }
1144     }
1145     cout << "truth event = " << m_evt << endl;
1146     m_evt++;
1147
1148     /*
1149       for (PHHepMCGenEventMap::ConstIter eventIter = hepmceventmap->begin();
1150            eventIter != hepmceventmap->end();
1151            ++eventIter)
1152       {
1153         /// Get the event
1154         PHHepMCGenEvent *hepmcevent = eventIter->second;
1155
1156         if (hepmcevent)
1157         {
1158           /// Get the parton info
1159           HepMC::PdfInfo *pdfinfo = truthevent->pdf_info();
1160
1161           /// Get the parton info as determined from HEPMC
1162           m_partid1 = pdfinfo->id1();
1163           m_partid2 = pdfinfo->id2();
1164           m_x1 = pdfinfo->x1();
1165           m_x2 = pdfinfo->x2();
1166
1167           /// Are there multiple partonic intercations in a p+p event
1168           m_mpi = truthevent->mpi();
1169
1170           /// Get the PYTHIA signal process id identifying the 2-to-2 hard process
1171           m_process_id = truthevent->signal_process_id();
1172         }
1173         // cout << "truth event = " << m_evt << endl;
1174         m_evt++;
1175       }
1176     */
1177   }
1178 }