Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:13:58 

0001 #include "InttAna.h"
0002 
0003 using namespace std;
0004 
0005 int InttAna::evtCount = 0;
0006 int InttAna::ievt = 0;
0007 
0008 InttAna::InttAna(const std::string &name, const std::string &filename)
0009     : SubsysReco(name), _rawModule(nullptr), fname_(filename), anafile_(nullptr), h_zvtxseed_(nullptr)
0010 {
0011   xbeam_ = ybeam_ = 0.0;
0012 
0013   // object
0014   zvtxobj_ = nullptr;
0015 
0016   // nodes
0017   inttrawmap = nullptr;
0018   m_tGeometry = nullptr;
0019   hepmceventmap = nullptr;
0020   phg4inevent = nullptr;
0021   m_clusterMap = nullptr;
0022   vertices = nullptr;
0023   inttevthead = nullptr;
0024   svtxvertexmap = nullptr;
0025   mbdout = nullptr;
0026   svtxvertex = nullptr;
0027   //  inttvertexmap = nullptr;
0028   intt_vertex_map = nullptr;
0029   svtxtrackmap = nullptr;
0030 
0031   // vector object
0032   vertex_pos_intt_ = new TVector3();
0033 }
0034 
0035 InttAna::~InttAna()
0036 {
0037 }
0038 
0039 int InttAna::Init(PHCompositeNode * /*topNode*/)
0040 {
0041   if (Verbosity() > 5)
0042   {
0043     std::cout << "Beginning Init in InttAna" << std::endl;
0044   }
0045 
0046   anafile_ = new TFile( fname_.c_str(), "recreate");
0047 
0048   this->InitHists();
0049   this->InitTrees();
0050   this->InitTuples();
0051 
0052   return 0;
0053 }
0054 
0055 void InttAna::InitHists()
0056 {
0057 
0058   h_dca2d_zero  = new TH1F("h_dca2d_zero", "DCA2D to 0", 500, -10, 10);
0059   h2_dca2d_zero = new TH2F("h2_dca2d_zero", "DCA2D to 0 vs nclus", 500, -10, 10, 100, 0, 10000);
0060   h2_dca2d_len  = new TH2F("h2_dca2d_len", "DCA2D to 0 vs len", 500, -10, 10, 500, -10, 10);
0061   h_zvtx    = new TH1F("h_zvtx", "Zvertex", 400, -40, 40);
0062   h_eta     = new TH1F("h_eta", "h_eta", 100, -6, 6);
0063   h_phi     = new TH1F("h_phi", "h_phi", 100, -6, 6);
0064   h_theta   = new TH1F("h_theta", "h_theta", 100, -4, 4);
0065 
0066   h_zvtxseed_ = new TH1F("h_zvtxseed", "Zvertex Seed histogram", 200, -50, 50);
0067 
0068 }
0069 
0070 void InttAna::InitTrees()
0071 {
0072 
0073   h_t_evt_bco = new TTree("t_evt_bco", "Event Tree for BCO");
0074   h_t_evt_bco->Branch("evt_gl1"     , &(m_evtbcoinfo.evt_gl1    ), "evt_gl1/I"      );
0075   h_t_evt_bco->Branch("evt_intt"    , &(m_evtbcoinfo.evt_intt   ), "evt_intt/i"     );
0076   h_t_evt_bco->Branch("evt_mbd"     , &(m_evtbcoinfo.evt_mbd    ), "evt_mbd/i"      );
0077   h_t_evt_bco->Branch("bco_gl1"     , &(m_evtbcoinfo.bco_gl1    ), "bco_gl1/l"      );
0078   h_t_evt_bco->Branch("bco_intt"    , &(m_evtbcoinfo.bco_intt   ), "bco_intt/l"     );
0079   h_t_evt_bco->Branch("bco_mbd"     , &(m_evtbcoinfo.bco_mbd    ), "bco_mbd/l"      );
0080   h_t_evt_bco->Branch("pevt_gl1"    , &(m_evtbcoinfo_prev.evt_gl1   ), "pevt_gl1/I"     );
0081   h_t_evt_bco->Branch("pevt_intt"   , &(m_evtbcoinfo_prev.evt_intt  ), "pevt_intt/i"    );
0082   h_t_evt_bco->Branch("pevt_mbd"    , &(m_evtbcoinfo_prev.evt_mbd   ), "pevt_mbd/i"     );
0083   h_t_evt_bco->Branch("pbco_gl1"    , &(m_evtbcoinfo_prev.bco_gl1   ), "pbco_gl1/l"     );
0084   h_t_evt_bco->Branch("pbco_intt"   , &(m_evtbcoinfo_prev.bco_intt  ), "pbco_intt/l"    );
0085   h_t_evt_bco->Branch("pbco_mbd"    , &(m_evtbcoinfo_prev.bco_mbd   ), "pbco_mbd/l"     );
0086 
0087   tree_event_ = new TTree( "tree_event", "Event-base TTree" );
0088   tree_event_->Branch( "nclus", &nclusadd_, "nclus/I" );
0089   tree_event_->Branch( "nclus_in", &nclus_inner_, "nclus_in/I" );
0090   tree_event_->Branch( "nclus_out", &nclus_outer_, "nclus_out/I" );
0091   tree_event_->Branch( "vertex_z_mbd", &vertex_z_mbd_, "vertex_z_mnd/D" );
0092   tree_event_->Branch( "vertex_intt", &vertex_pos_intt_ ); // TVector3  
0093   tree_event_->Branch( "bco_gl1"        , &(m_evtbcoinfo.bco_gl1    ), "bco_gl1/l"      );
0094   tree_event_->Branch( "bco_intt"   , &(m_evtbcoinfo.bco_intt   ), "bco_intt/l"     );
0095   
0096 
0097   m_hepmctree = new TTree("hepmctree", "A tree with hepmc truth particles");
0098   m_hepmctree->Branch("m_evt", &m_evt, "m_evt/I");
0099   m_hepmctree->Branch("m_xvtx", &m_xvtx, "m_xvtx/D");
0100   m_hepmctree->Branch("m_yvtx", &m_yvtx, "m_yvtx/D");
0101   m_hepmctree->Branch("m_zvtx", &m_zvtx, "m_zvtx/D");
0102   m_hepmctree->Branch("m_partid1", &m_partid1, "m_partid1/I");
0103   m_hepmctree->Branch("m_partid2", &m_partid2, "m_partid2/I");
0104   m_hepmctree->Branch("m_x1", &m_x1, "m_x1/D");
0105   m_hepmctree->Branch("m_x2", &m_x2, "m_x2/D");
0106   m_hepmctree->Branch("m_mpi", &m_mpi, "m_mpi/I");
0107   m_hepmctree->Branch("m_process_id", &m_process_id, "m_process_id/I");
0108   m_hepmctree->Branch("m_truthenergy", &m_truthenergy, "m_truthenergy/D");
0109   m_hepmctree->Branch("m_trutheta", &m_trutheta, "m_trutheta/D");
0110   m_hepmctree->Branch("m_truththeta", &m_truththeta, "m_truththeta/D");
0111   m_hepmctree->Branch("m_truthphi", &m_truthphi, "m_truthphi/D");
0112   m_hepmctree->Branch("m_status", &m_status, "m_status/I");
0113   m_hepmctree->Branch("m_truthpx", &m_truthpx, "m_truthpx/D");
0114   m_hepmctree->Branch("m_truthpy", &m_truthpy, "m_truthpy/D");
0115   m_hepmctree->Branch("m_truthpz", &m_truthpz, "m_truthpz/D");
0116   m_hepmctree->Branch("m_truthpt", &m_truthpt, "m_truthpt/D");
0117   m_hepmctree->Branch("m_numparticlesinevent", &m_numparticlesinevent, "m_numparticlesinevent/I");
0118   m_hepmctree->Branch("m_truthpid", &m_truthpid, "m_truthpid/I");
0119 
0120 }
0121 
0122 void InttAna::InitTuples()
0123 {
0124 
0125   h_ntp_clus = new TNtuple("ntp_clus",
0126                "Cluster Ntuple",
0127                (
0128                 string( "nclus:nclus2:bco_full:evt:size:adc:x:y:z:lay:" ) // 0-9
0129                 + "lad:sen:lx:ly:phisize:zsize:zv:x_vtx:y_vtx:z_vtx" // 10-19
0130                 ).c_str()
0131                );
0132   //  h_ntp_emcclus = new TNtuple("ntp_emcclus", "EMC Cluster Ntuple", "nemc:evt_emc:x_emc:y_emc:z_emc:e:ecore:size_emc");
0133 
0134   h_ntp_cluspair = new TNtuple("ntp_cluspair",
0135                    "Cluster Pair Ntuple",
0136                    "nclus:nclus2:bco_full:evt:ang1:ang2:z1:z2:dca2d:len:unorm:l1:l2:vx:vy:vz:zvtx" );
0137 
0138   //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");
0139 
0140   h_ntp_evt = new TNtuple("ntp_evt",
0141               "Event Ntuple",
0142               (
0143                string( "nclus:nclus2:bco_full:evt:zv:zvs:zvm:zvc:bz:bqn:" ) // 0-9
0144                + "bqs:bfemclk:xvtx:yvtx:zvtx:nclusin:nclusout:nclszv:ntrkzv:chi2ndfzv:" // 10-19
0145                + "widthzv:ngroupzv:goodzv:peakratiozv:xvsim:yvsim:zvsim:xvsvtx:yvsvtx:zvsvtx:" // 20-20
0146                + "nmvtx0:nmvtx1:nmvtx2:ntrksvtx:nemc:nemc1" //
0147                ).c_str()
0148               );
0149 
0150 }
0151 
0152 int InttAna::InitRun(PHCompositeNode * /*topNode*/)
0153 {
0154   if( Verbosity() > 1 )
0155     {
0156       cout << "InttAna::InitRun beamcenter " << xbeam_ << " " << ybeam_ << endl;
0157     }
0158   
0159   return 0;
0160 }
0161 
0162 int InttAna::GetNodes(PHCompositeNode *topNode)
0163 {
0164 
0165   m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0166   if (!m_tGeometry)
0167   {
0168     if( Verbosity() > 1 )
0169       {
0170     std::cerr << PHWHERE << "No ActsGeometry on node tree. Bailing." << std::endl;
0171       }
0172     return Fun4AllReturnCodes::ABORTEVENT;
0173   }
0174 
0175   inttrawmap = findNode::getClass<InttRawHitContainer>(topNode, "INTTRAWHIT");
0176   if( !inttrawmap )
0177     {
0178       if( Verbosity() > 1 )
0179     {
0180       cerr << PHWHERE << "INTTRAWHIT code is missing." << endl;
0181     }
0182     }
0183   
0184   // TODO:
0185   hepmceventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0186   if (!hepmceventmap)
0187   {
0188     if( Verbosity() > 3 )
0189       {
0190     cerr << PHWHERE << "hepmceventmap node is missing." << endl;
0191     // return Fun4AllReturnCodes::ABORTEVENT;
0192       }
0193   }
0194   
0195   // TODO:
0196   phg4inevent = findNode::getClass<PHG4InEvent>(topNode, "PHG4INEVENT");
0197   if (!phg4inevent)
0198     {
0199       if( Verbosity() > 3 )
0200     {
0201       cerr << PHWHERE << "PHG4INEVENT node is missing." << endl;
0202       // return Fun4AllReturnCodes::ABORTEVENT;
0203     }
0204     }
0205   
0206   m_clusterMap =
0207     findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0208   if (!m_clusterMap)
0209     {
0210       if( Verbosity() > 1 )
0211     {
0212       cerr << PHWHERE << "TrkrClusterContainer node is missing." << endl;
0213     }
0214       return Fun4AllReturnCodes::ABORTEVENT;
0215     }
0216   
0217   inttevthead = findNode::getClass<InttEventInfo>(topNode, "INTTEVENTHEADER");
0218   
0219   
0220   //gl1raw_ = findNode::getClass<Gl1RawHit>(topNode, "GL1RAWHIT");  
0221   gl1raw_ = findNode::getClass<Gl1Packetv2>(topNode, "GL1RAWHIT");  
0222   if( !gl1raw_ )
0223   {
0224     if( Verbosity() > 3 )
0225       {
0226     cerr << "No Gl1Raw" << endl;
0227       }
0228     
0229   }
0230 
0231   mbdout = findNode::getClass<MbdOut>(topNode, "MbdOut");
0232   if (!mbdout)
0233   {
0234     if( Verbosity() > 3 )
0235       {
0236     cerr << "MbdOut node is missing." << endl;
0237       }    
0238   }
0239 
0240   vertices =
0241       findNode::getClass<GlobalVertexMapv1>(topNode, "GlobalVertexMap");
0242   if (!vertices)
0243   {
0244     if( Verbosity() > 3 )
0245       {
0246     cerr << PHWHERE << "GlobalVertexMap node is missing." << endl;
0247       }
0248     // return Fun4AllReturnCodes::ABORTEVENT;
0249   }
0250   
0251 
0252   svtxvertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0253   if( !svtxvertexmap )
0254     {
0255       if( Verbosity() > 3 )
0256     {
0257       cerr << PHWHERE << "SvtxVertexMap node is missing." << endl;
0258     }
0259     }
0260 
0261   // inttvertexmap = findNode::getClass<InttVertex3DMap>(topNode, "InttVertexMap");
0262   // if( !inttvertexmap )
0263   //   {
0264   //     if( Verbosity() > 1 )
0265   //    {
0266   //      cerr << PHWHERE  << "InttVertex3DMap node is missing." << endl;
0267   //    }
0268       
0269   //     //return Fun4AllReturnCodes::ABORTEVENT;
0270   //   }
0271 
0272   intt_vertex_map = findNode::getClass<InttVertexMapv1>(topNode, "InttVertexMap");
0273   if( !intt_vertex_map )
0274     {
0275       if( Verbosity() > 1 )
0276     {
0277       cerr << PHWHERE  << "InttVertexMap node is missing." << endl;
0278     }
0279       
0280     }
0281 
0282   /////////////
0283   // SvtxTrack
0284   svtxtrackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0285 
0286   return Fun4AllReturnCodes::EVENT_OK;
0287 } // end of int InttAna::GetNodes(PHCompositeNode *topNode)
0288 
0289 int InttAna::process_event(PHCompositeNode *topNode)
0290 {
0291   
0292   if( Verbosity() > 1 )
0293     {
0294       cout << endl << "InttEvt::process evt : " << ievt++ << endl;
0295     }
0296   
0297   // readRawHit(topNode);
0298   this->GetNodes( topNode );    
0299   this->process_event_gl1( topNode );
0300   this->process_event_mbd( topNode );
0301   this->process_event_global_vertex( topNode );
0302   // this->process_event_svtx_vertex( topNode );
0303   
0304   this->process_event_intt_raw( topNode );
0305   this->process_event_intt_vertex( topNode ); // bco_ is overwritten, this should be before process_event_intt_cluster because vertex_pos_intt_ is needed in the function
0306   this->process_event_intt_cluster( topNode );
0307   this->process_event_intt_cluster_pair( topNode );
0308   // this->process_event_mvtx( topNode );
0309   // this->process_event_emcal( topNode );
0310   
0311   // fill evtbco_info
0312   m_evtbcoinfo.clear();
0313   m_evtbcoinfo.evt_gl1  = (gl1raw_  != nullptr)    ? gl1raw_->getEvtSequence()     : -1;
0314   m_evtbcoinfo.evt_intt = (inttraw_ != nullptr)    ? inttraw_->get_event_counter() : -1;
0315   m_evtbcoinfo.evt_mbd  = (mbdout   != nullptr)    ? mbdout->get_evt()             : -1;
0316   //m_evtbcoinfo.bco_gl1  = (gl1raw_  != nullptr)    ? gl1raw_->get_bco()            : 0;  
0317   m_evtbcoinfo.bco_gl1  = (gl1raw_  != nullptr)    ? gl1raw_->getBCO()            : 0;
0318 
0319   if( inttevthead != nullptr )
0320     m_evtbcoinfo.bco_intt = inttevthead->get_bco_full();
0321   else if( inttrawmap != nullptr )
0322     {
0323       if( inttrawmap->get_nhits() > 0 )
0324     m_evtbcoinfo.bco_intt = (inttrawmap->get_hit( 0 )->get_bco());
0325     }
0326   else
0327     m_evtbcoinfo.bco_intt =  0;
0328   
0329   m_evtbcoinfo.bco_mbd  = (mbdout   != nullptr)    ? mbdout->get_femclock()        : 0;
0330 
0331   if( Verbosity() > 3 )
0332     {
0333       cout << "event : "
0334        << m_evtbcoinfo.evt_gl1 << " "
0335        << m_evtbcoinfo.evt_intt << " "
0336        << m_evtbcoinfo.evt_mbd << " :  "
0337        << "INTT BCO " << setw( 20 ) << m_evtbcoinfo.bco_intt << "\t"
0338        << hex << m_evtbcoinfo.bco_gl1 << dec << " "
0339        << hex << m_evtbcoinfo.bco_intt << dec << " "
0340        << hex << m_evtbcoinfo.bco_mbd << dec << " : "
0341        << hex << m_evtbcoinfo.bco_gl1 - m_evtbcoinfo_prev.bco_gl1 << dec << " "
0342        << hex << m_evtbcoinfo.bco_intt - m_evtbcoinfo_prev.bco_intt << dec << " "
0343        << hex << m_evtbcoinfo.bco_mbd - m_evtbcoinfo_prev.bco_mbd << dec << endl;
0344 
0345       //cout << "evtCount : " << evtCount << " " << evtseq_ << " " << hex << bco_ << dec << endl;
0346     }
0347   
0348 
0349   this->process_event_fill( topNode );
0350   m_evtbcoinfo_prev.copy(m_evtbcoinfo);
0351 
0352   evtCount++;
0353 
0354   return Fun4AllReturnCodes::EVENT_OK;
0355 }
0356 
0357 int InttAna::ResetEvent(PHCompositeNode *topNode)
0358 {
0359 
0360   for( int i=0; i<10; i++ )
0361     for( int j=0; j<3; j++ )
0362       vertex_[i][j] = -9999;
0363 
0364   for( int i=0; i<3; i++ )
0365     {
0366       vtx_sim_[i] = -9999;
0367       nclusmvtx_[i] = 0;
0368     }
0369 
0370   // double 
0371   mbdqn_ = mbdqs_ = mbdz_ = 0;
0372   zvtx_
0373     = m_xvtx = m_yvtx = m_zvtx
0374     = m_x1 = m_x2
0375     = m_truthenergy = m_trutheta = m_truththeta = m_truthphi
0376     = m_truthpx = m_truthpy = m_truthpz = m_truthpt
0377     = m_truthp = m_vertex
0378     = vertex_z_mbd_
0379     = -9999.;
0380   
0381   // uint64
0382   bco_ = 0;
0383 
0384   // int
0385   evtseq_ = -1;
0386   nclusadd_ = nclusadd2_ = nclus_inner_ = nclus_outer_ = nemc_ = nemc1_
0387     = m_partid1 = m_partid2 = m_mpi = m_process_id = m_status
0388     = m_numparticlesinevent =  m_truthpid
0389     = 0;
0390 
0391   // vector
0392   for( auto& cluster : clusters_ )
0393     cluster.erase( cluster.begin(), cluster.end() );
0394 
0395   return Fun4AllReturnCodes::EVENT_OK;
0396 }
0397 
0398 int InttAna::End(PHCompositeNode * /*topNode*/)
0399 {
0400   if (Verbosity() > 1)
0401   {
0402     std::cout << "Ending InttAna analysis package" << std::endl;
0403   }
0404 
0405   //  m_hepmctree->Write();
0406 
0407   if (anafile_ != nullptr)
0408   {
0409     //anafile_->Write();
0410     anafile_->WriteTObject( h_dca2d_zero, h_dca2d_zero->GetName() );
0411     anafile_->WriteTObject( h2_dca2d_zero, h2_dca2d_zero->GetName() );
0412     anafile_->WriteTObject( h2_dca2d_len, h2_dca2d_len->GetName() );
0413     anafile_->WriteTObject( h_zvtx, h_zvtx->GetName() );
0414     anafile_->WriteTObject( h_eta, h_eta->GetName() );
0415     anafile_->WriteTObject( h_phi, h_phi->GetName() );
0416     anafile_->WriteTObject( h_theta, h_theta->GetName() );
0417 
0418     anafile_->WriteTObject( h_ntp_clus, h_ntp_clus->GetName() );
0419     //    anafile_->WriteTObject( h_ntp_emcclus, h_ntp_emcclus->GetName() );
0420     anafile_->WriteTObject( h_ntp_cluspair, h_ntp_cluspair->GetName() );
0421     //    anafile_->WriteTObject( h_ntp_emccluspair, h_ntp_emccluspair->GetName() );
0422     anafile_->WriteTObject( h_ntp_evt, h_ntp_evt->GetName() );
0423     anafile_->WriteTObject( h_zvtxseed_, h_zvtxseed_->GetName() );
0424     anafile_->WriteTObject( h_t_evt_bco, h_t_evt_bco->GetName() );
0425     anafile_->WriteTObject( m_hepmctree, m_hepmctree->GetName() );
0426     anafile_->WriteTObject( tree_event_, tree_event_->GetName() );
0427     anafile_->Close();
0428   }
0429 
0430   return 0;
0431 }
0432 
0433 int InttAna::process_event_gl1(PHCompositeNode *topNode )
0434 {
0435 
0436   if (gl1raw_)
0437   {
0438     //bco_ = gl1raw_->get_bco();
0439     bco_ = gl1raw_->getBCO();
0440     evtseq_ = gl1raw_->getEvtSequence();
0441   }
0442 
0443   return 0;
0444 }
0445 
0446 int InttAna::process_event_mbd(PHCompositeNode *topNode )
0447 {
0448   if( !mbdout )
0449     return 0;
0450 
0451   mbdqs_ = (mbdout != nullptr) ? mbdout->get_q(0) : -9999;
0452   mbdqn_ = (mbdout != nullptr) ? mbdout->get_q(1) : -9999;
0453   //  double mbdt0 = (mbdout!=nullptr) ?  mbdout->get_t0() : -9999;
0454   mbdz_ = (mbdout != nullptr) ? mbdout->get_zvtx() : -9999;
0455   
0456   if( mbdout != nullptr && mbdout->isValid() )
0457     {
0458       vtx_sim_[2] = mbdout->get_zvtx();
0459     }
0460   
0461   // MbdVertexMap *mbdvertexmap = findNode::getClass<MbdVertexMap>(topNode, "MbdVertexMap");
0462   // MbdVertex* mbdvtx = nullptr;
0463   // if (mbdvertexmap)
0464   //   {
0465   //     if( Verbosity() > 1 )
0466   //    {
0467   //      cout << "MbdVertexMap" << endl;
0468   //    }
0469       
0470   //   if (!mbdvertexmap->empty())
0471   //     {
0472     
0473   //    for(auto mbdvertex : *mbdvertexmap)
0474   //      {
0475   //        if( Verbosity() > 1 )
0476   //          {
0477   //        std::cout << "mbdmap entry" << std::endl;
0478   //        std::cout << std::endl << "Mbdz"
0479   //              << " "  << mbdvertex.second->get_z()
0480   //              << " "  << mbdvertex.second->get_id() << endl;
0481   //          }
0482         
0483   //        mbdvtx = mbdvertex.second;
0484   //      }
0485   //     }
0486   //   else if( Verbosity() > 1 )
0487   //     {
0488   //    cout << "MbdVertexMap empty" << endl;
0489   //     }
0490   //   }
0491   
0492   return 0;
0493 }
0494 
0495 
0496 int InttAna::process_event_global_vertex(PHCompositeNode *topNode )
0497 {
0498 
0499   if( !vertices )
0500     return 0;
0501 
0502   std::map < GlobalVertex::VTXTYPE, std::string > sourcemap
0503     {
0504       { GlobalVertex::VTXTYPE::UNDEFINED, "UNDEFINED" },
0505       { GlobalVertex::VTXTYPE::TRUTH,     "TRUTH"     },
0506       { GlobalVertex::VTXTYPE::SMEARED,   "SMEARED"   },
0507       { GlobalVertex::VTXTYPE::MBD,       "MBD"       },
0508       { GlobalVertex::VTXTYPE::SVTX,      "SVTX"      },
0509       { GlobalVertex::VTXTYPE::SVTX_MBD,  "SVTX_MBD"  }
0510     };
0511   
0512   //   enum VTXTYPE
0513   std::string s_vtxid = "NULL";
0514 
0515   if (vertices)
0516   {
0517 
0518     vertices->identify();
0519     
0520     if (!vertices->empty())
0521     {
0522 
0523       for (auto vertex : *vertices)
0524       {
0525         std::map<GlobalVertex::VTXTYPE, std::string>::iterator itr_src;
0526         itr_src = sourcemap.find((GlobalVertex::VTXTYPE)vertex.first);
0527         if (itr_src != sourcemap.end())
0528         {
0529           s_vtxid = itr_src->second;
0530         }
0531 
0532         if (vertex.first == GlobalVertex::VTXTYPE::TRUTH)
0533         {
0534           vtx_sim_[0] = vertex.second->get_x();
0535           vtx_sim_[1] = vertex.second->get_y();
0536           vtx_sim_[2] = vertex.second->get_z();
0537         }
0538     else if( vertex.first == GlobalVertex::VTXTYPE::MBD )
0539       {
0540         vertex_z_mbd_ = vertex.second->get_z();
0541 
0542       }
0543         if( Verbosity() > 2 )
0544       {
0545         std::cout << "GlobalVertex map entry" << std::endl;
0546         std::cout << std::endl
0547               << " " << vertex.second->get_x()
0548               << " " << vertex.second->get_y()
0549               << " " << vertex.second->get_z()
0550               << " " << vertex.second->get_id()
0551               << " " << s_vtxid << std::endl;
0552       } // end of if( Verbosity() > 2 )
0553     
0554       } // end of for (auto vertex : *vertices)
0555     } // end of if (!vertices->empty())
0556   } // if (vertices)
0557 
0558 
0559   if (hepmceventmap != nullptr || phg4inevent != nullptr)
0560   {
0561     xvtx_sim = vertices->find(GlobalVertex::TRUTH)->second->get_x();
0562     yvtx_sim = vertices->find(GlobalVertex::TRUTH)->second->get_y();
0563     zvtx_sim = vertices->find(GlobalVertex::TRUTH)->second->get_z();
0564     
0565     // TODO:
0566     if (hepmceventmap != nullptr)
0567       getHEPMCTruth(topNode);
0568     else if (phg4inevent != nullptr)
0569       getPHG4Particle(topNode);
0570   }
0571   
0572   
0573   return 0;
0574 }
0575 
0576 int InttAna::process_event_svtx_vertex(PHCompositeNode *topNode )
0577 {
0578 
0579   svtxvertex = nullptr;
0580   if( svtxvertexmap )
0581   {
0582     if( Verbosity() > 3 )
0583       {
0584     cout << "svtxvertex: size : " << svtxvertexmap->size() << endl;
0585       }
0586     
0587     for (SvtxVertexMap::ConstIter iter = svtxvertexmap->begin(); iter != svtxvertexmap->end(); ++iter)
0588       {
0589     svtxvertex = iter->second;
0590     if( Verbosity() > 2 )
0591       {
0592         std::cout << std::endl
0593               << "SvtxVertex"
0594               << " " << svtxvertex->get_x()
0595               << " " << svtxvertex->get_y()
0596               << " " << svtxvertex->get_z()
0597               << " " << svtxvertex->get_id() << endl;
0598       }
0599     
0600       }
0601   }
0602 
0603   return 0;
0604 }
0605 
0606 
0607 int InttAna::process_event_intt_raw(PHCompositeNode *topNode )
0608 {
0609 
0610   if( inttrawmap == nullptr )
0611     return 0;
0612   
0613   inttraw_ = (inttrawmap != nullptr && inttrawmap->get_nhits() > 0) ? inttrawmap->get_hit(0) : nullptr;
0614 
0615   //cout << "#INTT raw hit: " << inttrawmap->get_nhits() << endl;
0616   
0617   /////////////
0618   //  InttEvent from RawModule (not standard way)
0619   InttEvent *inttEvt = nullptr;
0620   if (_rawModule)
0621   {
0622     inttEvt = _rawModule->getInttEvent();
0623   }
0624   
0625   bco_ = 0;
0626   evtseq_ = -1;
0627   
0628   if (inttEvt != NULL)
0629   {
0630     bco_ = inttEvt->bco;
0631     evtseq_ = inttEvt->evtSeq;
0632   }
0633 
0634   if (inttevthead)
0635   {
0636     bco_ = inttevthead->get_bco_full();
0637     // evtseq = inttEvt->evtSeq;
0638   }
0639 
0640   return 0;
0641 }
0642 
0643 int InttAna::process_event_intt_cluster(PHCompositeNode *topNode )
0644 {
0645 
0646   ////////////////////////////////
0647   // single cluster analysis    //
0648   ////////////////////////////////
0649   float ntpval[20];
0650   int nCluster = 0;
0651   bool exceedNwrite = false;
0652   // std::vector<Acts::Vector3> clusters_[2]; // inner=0; outer=1
0653   
0654   for( unsigned int inttlayer = 0; inttlayer < 4; inttlayer++ )
0655   {
0656 
0657     int layer = (inttlayer < 2 ? 0 : 1);
0658     
0659     for( const auto &hitsetkey : m_clusterMap->getHitSetKeys(TrkrDefs::TrkrId::inttId, inttlayer + 3) )
0660     {
0661       
0662       auto range = m_clusterMap->getClusters(hitsetkey);
0663 
0664       for( auto clusIter = range.first; clusIter != range.second; ++clusIter )
0665       {
0666 
0667         nclusadd_++;
0668 
0669         const auto cluskey = clusIter->first;
0670         const auto cluster = clusIter->second;
0671 
0672     if (cluster->getAdc() > 40)
0673           nclusadd2_++;
0674     
0675         if (inttlayer < 2)
0676           nclus_inner_++;
0677         else
0678           nclus_outer_++;
0679 
0680         int ladder_z   = InttDefs::getLadderZId(cluskey);
0681         int ladder_phi = InttDefs::getLadderPhiId(cluskey);
0682         int size       = cluster->getSize();
0683 
0684         // cout    <<cluster->getPosition(0)<<" "
0685         //         <<cluster->getPosition(1)<<" : "
0686         //         <<cluster->getAdc()<<" "<<size<<" "<<inttlayer<<" "<<ladder_z<<" "<<ladder_phi<<endl;
0687 
0688         const auto globalPos = m_tGeometry->getGlobalPosition(cluskey, cluster);
0689 
0690         if (nCluster < 5)
0691         {
0692           if( Verbosity() > 2 )
0693         {
0694           cout << "cluster xyz : " << globalPos.x() << " " << globalPos.y() << " " << globalPos.z() << " :  "
0695            << cluster->getPosition(0) << " "
0696            << cluster->getPosition(1) << " : "
0697            << cluster->getAdc() << " " << size << " " << inttlayer << " " << ladder_z << " " << ladder_phi << endl;
0698         }
0699       
0700         }
0701         else
0702         {
0703           if (!exceedNwrite)
0704           {
0705             if( Verbosity() > 2 )
0706           {
0707         cout << " exceed : ncluster limit.  no more cluster xyz printed" << endl;
0708           }
0709         
0710             exceedNwrite = true;
0711           }
0712         }
0713 
0714         ClustInfo info;
0715         info.layer = inttlayer;
0716         info.pos   = globalPos;
0717 
0718         // clusters[layer].push_back(globalPos);
0719         clusters_[layer].push_back(info);
0720         nCluster++;
0721 
0722         ntpval[0] = nclusadd_;          // bco_full
0723         ntpval[1] = nclusadd2_;         // bco_full
0724         ntpval[2] = bco_;               // bco_full
0725         ntpval[3] = evtseq_;            // evtCount;
0726         ntpval[4] = size;              // size
0727         ntpval[5] = cluster->getAdc(); // size
0728         ntpval[6] = globalPos.x();
0729         ntpval[7] = globalPos.y();
0730         ntpval[8] = globalPos.z();
0731         ntpval[9] = inttlayer;
0732         ntpval[10] = ladder_phi;
0733         ntpval[11] = ladder_z;
0734         ntpval[12] = cluster->getPosition(0);
0735         ntpval[13] = cluster->getPosition(1);
0736         ntpval[14] = cluster->getPhiSize();
0737         ntpval[15] = cluster->getZSize();
0738         ntpval[16] = vertex_[2][2]; // zvtx;
0739         // ntpval[17] = xvtx_sim;
0740         // ntpval[18] = yvtx_sim;
0741         // ntpval[19] = zvtx_sim;
0742     ntpval[17] = vertex_pos_intt_->X();
0743     ntpval[18] = vertex_pos_intt_->Y();
0744     ntpval[19] = vertex_pos_intt_->Z();
0745 
0746     h_ntp_clus->Fill(ntpval);
0747         // = new TNtuple("ntp_clus", "Cluster Ntuple", "bco_full:evt:size:adc:x:y:z:lay:lad:sen");
0748       }
0749     }
0750   }
0751 
0752   //  cout << "#INTT cluster: " << nclusadd_ << endl;
0753 
0754   return 0;
0755 }
0756 
0757 int InttAna::process_event_intt_cluster_pair(PHCompositeNode *topNode )
0758 {
0759 
0760   // cluster pair analysis
0761   struct cluspair
0762   {
0763     float p1_ang;
0764     float p2_ang;
0765     float p1_r;
0766     float p2_r;
0767     float p1_z;
0768     float p2_z;
0769     float dca_p0;
0770     float len_p0;
0771   };
0772 
0773   vector<cluspair> vcluspair;
0774 
0775   int nCluster = nclusadd_;
0776   if (nCluster < 300)
0777   // if(nCluster<500)
0778   {
0779     // if( Verbosity() > 2 )
0780     //   {
0781     //  cout << "#cluster " << nCluster << " is less than 300." << endl;
0782     //   }
0783     
0784     float ntpval2[20];
0785     Acts::Vector3 beamspot = Acts::Vector3(xbeam_, ybeam_, 0);
0786     vector<double> vz_array;
0787     for (auto c1 = clusters_[0].begin(); c1 != clusters_[0].end(); ++c1) // inner
0788     {
0789 
0790       for (auto c2 = clusters_[1].begin(); c2 != clusters_[1].end(); ++c2) // outer
0791       {
0792         Acts::Vector3 p1 = c1->pos - beamspot;
0793         Acts::Vector3 p2 = c2->pos - beamspot;
0794         Acts::Vector3 u = p2 - p1;
0795         double unorm = sqrt(u.x() * u.x() + u.y() * u.y());
0796         if (unorm < 0.00001)
0797           continue;
0798 
0799         // skip bad compbination
0800         double p1_ang = atan2(p1.y(), p1.x());
0801         double p2_ang = atan2(p2.y(), p2.x());
0802         double d_ang = p2_ang - p1_ang;
0803 
0804         // unit vector in 2D
0805         double ux = u.x() / unorm;
0806         double uy = u.y() / unorm;
0807         double uz = u.z() / unorm; // same normalization factor(2D) is multiplied
0808 
0809         // Acts::Vector3 p0   = beamspot - p1;
0810         Acts::Vector3 p0 = Acts::Vector3(0, 0, 0) - p1; // p1, p2 are already shifted by beam center
0811 
0812         double dca_p0 = p0.x() * uy - p0.y() * ux; // cross product of p0 x u
0813         double len_p0 = p0.x() * ux + p0.y() * uy; // dot product of p0 . u
0814 
0815         // beam center in X-Y plane
0816         double vx = len_p0 * ux + p1.x();
0817         double vy = len_p0 * uy + p1.y();
0818 
0819         double vz = len_p0 * uz + p1.z();
0820 
0821         double p1_r = sqrt(p1.y() * p1.y() + p1.x() * p1.x());
0822         double p2_r = sqrt(p2.y() * p2.y() + p2.x() * p2.x());
0823 
0824         // if(unorm>4.5||d_ang<-0.005*3.1415||0.005*3.1415<d_ang)
0825         // if(unorm>4.5||d_ang<-0.25*3.1415||0.25*3.1415<d_ang)
0826         if (fabs(d_ang) < 0.05 && fabs(dca_p0) < 0.5)
0827         {
0828           h_zvtxseed_->Fill(vz);
0829           vz_array.push_back(vz);
0830 
0831           cluspair clspair;
0832           clspair.p1_ang = p1_ang;
0833           clspair.p2_ang = p2_ang;
0834           clspair.p1_r = p1_r;
0835           clspair.p2_r = p2_r;
0836           clspair.p1_z = p1.z();
0837           clspair.p2_z = p2.z();
0838           clspair.dca_p0 = dca_p0;
0839           clspair.len_p0 = len_p0;
0840           vcluspair.push_back(clspair);
0841         }
0842 
0843         h_dca2d_zero->Fill(dca_p0);
0844         h2_dca2d_zero->Fill(dca_p0, nCluster);
0845         h2_dca2d_len->Fill(dca_p0, len_p0);
0846         // cout<<"dca : "<<dca_p0<<endl;
0847         //
0848         ntpval2[0] = nclusadd_;
0849         ntpval2[1] = nclusadd2_;
0850         ntpval2[2] = 0;
0851         ntpval2[3] = evtCount;
0852         ntpval2[4] = p1_ang;
0853         ntpval2[5] = p2_ang;
0854         ntpval2[6] = p1.z();
0855         ntpval2[7] = p2.z();
0856         ntpval2[8] = dca_p0;
0857         ntpval2[9] = len_p0;
0858         ntpval2[10] = unorm;
0859         ntpval2[11] = c1->layer;
0860         ntpval2[12] = c2->layer;
0861         ntpval2[13] = vx;
0862         ntpval2[14] = vy;
0863         ntpval2[15] = vz;
0864         ntpval2[16] = vertex_[2][2]; // zvtx;
0865 
0866         if (evtCount < 10000)
0867           h_ntp_cluspair->Fill(ntpval2);
0868         // h_ntp_cluspair = new TNtuple("ntp_cluspair", "Cluster Pair Ntuple", "nclus:bco_full:evt:ang1:ang2:dca2d:len:unorm");
0869       } // end of loop over 2nd cluster
0870     } // end of loop over 1st cluster
0871 
0872     if( Verbosity() > 2 )
0873       {
0874     cout << "End of loop over cluster pairs" << endl;
0875       }
0876     
0877     // calculate trancated mean of DCA~Z histogram as Z-vertex position
0878 
0879     if (vz_array.size() > 3)
0880     {
0881       if( Verbosity() > 2 )
0882     {
0883       cout << "zvtx determination as vz_array size " << vz_array.size() << " is larger than 3" << endl;
0884     }
0885       
0886       double zbin = h_zvtxseed_->GetMaximumBin();
0887       double zcenter = h_zvtxseed_->GetBinCenter(zbin);
0888       double zmean = h_zvtxseed_->GetMean();
0889       double zrms = h_zvtxseed_->GetRMS();
0890       if (zrms < 20)
0891         zrms = 20;
0892 
0893       double zmax = zcenter + zrms; // 1 sigma
0894       double zmin = zcenter - zrms; // 1 sigma
0895 
0896       double zsum = 0.;
0897       int zcount = 0;
0898       for (auto iz = vz_array.begin(); iz != vz_array.end(); ++iz)
0899       {
0900         double vz = (*iz);
0901         if (zmin < vz && vz < zmax)
0902         {
0903           zsum += vz;
0904           zcount++;
0905         }
0906       }
0907       if (zcount > 0)
0908         zvtx_ = zsum / zcount;
0909 
0910       if( Verbosity() > 2 )
0911     {
0912       cout << "ZVTX: " << zvtx_ << " " << zcenter << " " << zmean << " " << zrms << " " << zbin << endl;
0913     }
0914       
0915     }
0916 
0917     h_zvtx->Fill(zvtx_);
0918 
0919     // float ntpval_emcpair[20];
0920     // if (clusterContainer != nullptr)
0921     // {
0922     //   RawClusterContainer::ConstRange clusterEnd = clusterContainer->getClusters();
0923     //   RawClusterContainer::ConstIterator clusterIter;
0924 
0925     //   /////////////
0926     //   //  association to EMCAL
0927     //   for (vector<cluspair>::iterator itrcp = vcluspair.begin(); itrcp != vcluspair.end(); ++itrcp)
0928     //   {
0929     //     cluspair &clspair = *itrcp;
0930 
0931     //     for (clusterIter = clusterEnd.first; clusterIter != clusterEnd.second; clusterIter++)
0932     //     {
0933     //       RawCluster *recoCluster = clusterIter->second;
0934 
0935     //       float the1 = atan2(clspair.p1_r, (clspair.p1_z - vertex_[2][2]));
0936     //       float the2 = atan2(clspair.p2_r, (clspair.p2_z - vertex_[2][2]));
0937 
0938     //       Acts::Vector3 emcpos = Acts::Vector3(recoCluster->get_x(), recoCluster->get_y(), recoCluster->get_z());
0939     //       Acts::Vector3 emcposc = emcpos - beamspot;
0940     //       float ephi = atan2(emcposc.y(), emcposc.x());
0941     //       float er = sqrt(emcposc.x() * emcposc.x() + emcposc.y() * emcposc.y());
0942     //       float ethe = atan2(er, emcposc.z() - vertex_[2][2]);
0943 
0944     //       ntpval_emcpair[0] = nclusadd;
0945     //       ntpval_emcpair[1] = evtseq;
0946     //       ntpval_emcpair[2] = clspair.p1_ang; // INTT
0947     //       ntpval_emcpair[3] = clspair.p2_ang;
0948     //       ntpval_emcpair[4] = clspair.p1_z;
0949     //       ntpval_emcpair[5] = clspair.p2_z;
0950     //       ntpval_emcpair[6] = clspair.dca_p0;
0951     //       ntpval_emcpair[7] = clspair.len_p0;
0952     //       ntpval_emcpair[8] = vertex_[1][0];  // x-vertex
0953     //       ntpval_emcpair[9] = vertex_[1][1];  // y-vertex
0954     //       ntpval_emcpair[10] = vertex_[2][2]; // y-vertex
0955     //       ntpval_emcpair[11] = ephi;
0956     //       ntpval_emcpair[12] = recoCluster->get_z(); // EMCal
0957     //       ntpval_emcpair[13] = recoCluster->get_ecore();
0958     //       ntpval_emcpair[14] = recoCluster->getNTowers();
0959     //       ntpval_emcpair[15] = the1;
0960     //       ntpval_emcpair[16] = the2;
0961     //       ntpval_emcpair[17] = ethe;
0962 
0963     //       h_ntp_emccluspair->Fill(ntpval_emcpair);
0964     //     }
0965     //   }
0966     // }
0967   }
0968   
0969   return 0;
0970 }
0971 
0972 int InttAna::process_event_intt_vertex(PHCompositeNode *topNode )
0973 {
0974 
0975   vertex_pos_intt_ = new TVector3();
0976 
0977   if( intt_vertex_map )
0978     {
0979       if( Verbosity() > 2 )
0980     {
0981       cout << "#vertex from InttVertexMapV1: " << intt_vertex_map->size() << endl;
0982       intt_vertex_map->identify();
0983     }
0984       
0985       InttVertexMap::ConstIter iter_begin = intt_vertex_map->begin();
0986       InttVertexMap::ConstIter iter_end = intt_vertex_map->end();
0987 
0988       int ivtx = 0;
0989       for( auto iter=iter_begin; iter!= iter_end; iter++, ivtx++ )
0990     {
0991       if( ivtx >= 10 )
0992         {         
0993           break;
0994         }
0995       
0996       InttVertex* intt_vertex = iter->second;
0997       if( Verbosity() > 2 )
0998       {
0999         cout << endl << "INTT vertex map\t"
1000          << ivtx << "\t"          
1001          << intt_vertex->isValid() << "\t"
1002            << "( "
1003          << setw(8) << setprecision(5) << intt_vertex->get_x() << ", "
1004          << setw(8) << setprecision(5) << intt_vertex->get_y() << ", "
1005          << setw(8) << setprecision(5) << intt_vertex->get_z()
1006          << " )"
1007          << endl;
1008       }
1009       
1010       
1011       if( intt_vertex )
1012         {
1013           //intt_vertex->identify();
1014           vertex_[ ivtx ][ 0 ] = intt_vertex->get_x();
1015           vertex_[ ivtx ][ 1 ] = intt_vertex->get_y();
1016           vertex_[ ivtx ][ 2 ] = intt_vertex->get_z();
1017           if( ivtx == 2 )
1018         {
1019           zvtxobj_ = intt_vertex;
1020           vertex_pos_intt_->SetXYZ( intt_vertex->get_x(), 
1021                         intt_vertex->get_y(),
1022                         vertex_[2][2] );
1023         }
1024         }
1025       
1026       // if( iter == iter_begin )
1027       //   {
1028       //     vertex_pos_intt_->SetXYZ( intt_vertex->get_x(), 
1029       //                intt_vertex->get_y(),
1030       //                intt_vertex->get_z() );
1031       //   }
1032     }
1033       
1034       
1035     }
1036   
1037   // //  zvtxobj_ = NULL;
1038   // if( inttvertexmap )
1039   //   {
1040   //     int ivtx = 0;
1041   //     if( Verbosity() > 2 )
1042   //    {
1043   //      cout << "vertex map size : " << inttvertexmap->size() << endl;
1044   //    }
1045       
1046   //     InttVertex3DMap::ConstIter biter_beg = inttvertexmap->begin();
1047   //     InttVertex3DMap::ConstIter biter_end = inttvertexmap->end();
1048   //     // cout<<"vertex map size : after size "<<endl;
1049   //     // inttvertexmap->identify();
1050       
1051   //     for (InttVertex3DMap::ConstIter biter = biter_beg; biter != biter_end; ++biter)
1052   //    {
1053   //      InttVertex3D *vtxobj = biter->second;
1054   //      if (vtxobj)
1055   //        {
1056   //          if( Verbosity() > 2 )
1057   //        {
1058   //          cout << "vtxobj" << ivtx << endl;
1059   //        }
1060           
1061   //          vtxobj->identify();
1062   //          vertex_[ivtx][0] = vtxobj->get_x();
1063   //          vertex_[ivtx][1] = vtxobj->get_y();
1064   //          vertex_[ivtx][2] = vtxobj->get_z();
1065   //          // if (ivtx == 2)
1066   //          //    {
1067   //          //      zvtxobj_ = vtxobj;
1068   //          //    }
1069   //        }
1070   //      else
1071   //        {
1072   //          if( Verbosity() > 2 )
1073   //        {
1074   //          cout << "no vtx object : " << ivtx << endl;
1075   //        }
1076   //        }
1077   //      ivtx++;
1078   //      if (ivtx >= 10)
1079   //        {
1080   //          if( Verbosity() > 2 )
1081   //        {
1082   //          cout << "n_vertex exceeded : " << ivtx << endl;
1083   //        }
1084           
1085   //          break;
1086   //        }
1087   //    }
1088   //   }
1089   // else
1090   //   {
1091   //     // zvtxobj_ = new InttVertex3D(); // impossible as constructor protected
1092   //   }
1093 
1094   return 0;
1095 }
1096 
1097 
1098 int InttAna::process_event_mvtx(PHCompositeNode *topNode )
1099 {
1100 
1101   /////////////
1102   // MVTX cluster
1103   std::set<TrkrDefs::TrkrId> detectors;
1104   detectors.insert(TrkrDefs::TrkrId::mvtxId);
1105 
1106   for (const auto &det : detectors)
1107   {
1108     for (const auto &layer : {0, 1, 2})
1109     {
1110       for (const auto &hitsetkey : m_clusterMap->getHitSetKeys(det, layer))
1111       {
1112         auto range = m_clusterMap->getClusters(hitsetkey);
1113         for (auto citer = range.first; citer != range.second; ++citer)
1114         {
1115           // const auto ckey = citer->first;
1116           // const auto cluster = citer->second;
1117           // const auto global = m_tGeometry->getGlobalPosition(ckey, cluster);
1118           nclusmvtx_[layer]++;
1119         }
1120       }
1121     }
1122   }
1123 
1124   return 0;
1125 }
1126 
1127 
1128 int InttAna::process_event_emcal(PHCompositeNode *topNode )
1129 {
1130 
1131   /////////////
1132   // EMC Cluster
1133   //RawClusterContainer *clusterContainer = findNode::getClass<RawClusterContainer>(topNode, "CLUSTERINFO_POS_COR_CEMC");
1134   // RawClusterContainer *clusterContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_POS_COR_CEMC");
1135   // RawClusterContainer *clusterContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_CEMC");
1136 
1137   // if (!clusterContainer)
1138   // {
1139   //   std::cout << PHWHERE << "InttAna::process_event - Fatal Error - CLUSTER_POS_COR_CEMC node is missing. " << std::endl;
1140   // }
1141   // else
1142   // {
1143   //   float ntpval_emccls[20];
1144   //   // This is how we iterate over items in the container.
1145   //   RawClusterContainer::ConstRange clusterEnd = clusterContainer->getClusters();
1146   //   RawClusterContainer::ConstIterator clusterIter;
1147 
1148   //   nemc = clusterContainer->size();
1149   //   nemc1 = 0;
1150 
1151   //   for (clusterIter = clusterEnd.first; clusterIter != clusterEnd.second; clusterIter++)
1152   //   {
1153   //     RawCluster *recoCluster = clusterIter->second;
1154 
1155   //     ntpval_emccls[0] = nemc;
1156   //     ntpval_emccls[1] = evtseq;
1157   //     ntpval_emccls[2] = recoCluster->get_x();
1158   //     ntpval_emccls[3] = recoCluster->get_y();
1159   //     ntpval_emccls[4] = recoCluster->get_z();
1160   //     ntpval_emccls[5] = recoCluster->get_energy();
1161   //     ntpval_emccls[6] = recoCluster->get_ecore();
1162   //     ntpval_emccls[7] = recoCluster->getNTowers();
1163 
1164   //     h_ntp_emcclus->Fill(ntpval_emccls);
1165 
1166   //     if (recoCluster->get_ecore() > 0.15)
1167   //       nemc1++;
1168   //   }
1169   // }
1170 
1171   return 0;
1172 }
1173 
1174 int InttAna::process_event_fill(PHCompositeNode *topNode )
1175 {
1176 
1177   if( Verbosity() > 2 )
1178     {
1179       cout << "value assignments for fill" << endl;
1180     }
1181   
1182   float ntpval3[40] = { -9999.9 };
1183   ntpval3[0] = nclusadd_;  // bco_full
1184   ntpval3[1] = nclusadd2_; // bco_full
1185   ntpval3[2] = bco_;
1186   ntpval3[3] = evtseq_;
1187   ntpval3[4] = vertex_[2][2]; // zvtx;
1188   ntpval3[5] = zvtx_;         // zrms;
1189   ntpval3[6] = 0;            // zmean;
1190   ntpval3[7] = 0;            // zcenter;
1191   ntpval3[8] = mbdz_;
1192   ntpval3[9] = mbdqn_;         // mbdt.bqn;
1193 
1194   ntpval3[10] = mbdqs_;        // mbdt.bqs;
1195   ntpval3[11] = 0;            // mbdt.femclk;
1196   ntpval3[12] = vertex_[1][0]; // x-vertex
1197   ntpval3[13] = vertex_[1][1]; // y-vertex
1198   ntpval3[14] = vertex_[2][2]; // y-vertex
1199 
1200 
1201   ntpval3[15] = nclus_inner_;  // 
1202   ntpval3[16] = nclus_outer_;  //
1203 
1204   if( zvtxobj_ != nullptr )
1205     {
1206       ntpval3[17] = zvtxobj_->get_nclus();
1207       ntpval3[18] = zvtxobj_->get_ntracklet();
1208       ntpval3[19] = zvtxobj_->get_chi2ndf();
1209       
1210       ntpval3[20] = zvtxobj_->get_width();
1211       ntpval3[21] = zvtxobj_->get_ngroup();
1212       ntpval3[22] = (zvtxobj_->get_good()) ? 1 : 0;
1213       ntpval3[23] = zvtxobj_->get_peakratio();
1214 
1215     }
1216   else
1217     {
1218 
1219       ntpval3[17] = ntpval3[18] = ntpval3[19]
1220     = ntpval3[20] = ntpval3[21] = ntpval3[22]
1221     = ntpval3[23]
1222     = 0.0;
1223 
1224     }
1225   
1226   ntpval3[24] = vtx_sim_[0]; // sim x-vertex
1227   ntpval3[25] = vtx_sim_[1]; // sim y-vertex
1228   ntpval3[26] = vtx_sim_[2]; // sim z-vertex
1229 
1230 
1231   if( svtxvertex != nullptr )
1232     {
1233 
1234       ntpval3[27] = svtxvertex->get_x();
1235       ntpval3[28] = svtxvertex->get_y();
1236       ntpval3[29] = svtxvertex->get_z();
1237 
1238     }
1239   else
1240     {
1241 
1242       ntpval3[27] = ntpval3[28] = ntpval3[29]
1243     = -9999;
1244     }
1245 
1246 
1247   ntpval3[30] = nclusmvtx_[0];
1248   ntpval3[31] = nclusmvtx_[1];
1249   ntpval3[32] = nclusmvtx_[2];
1250   ntpval3[33] = (svtxtrackmap != nullptr) ? svtxtrackmap->size() : -9999;
1251   // ntpval3[34] = nemc_;
1252   // ntpval3[35] = nemc1_;
1253   ntpval3[34] = 0;
1254   ntpval3[35] = 0;
1255 
1256   h_ntp_evt->Fill(ntpval3);
1257   h_t_evt_bco->Fill();
1258   tree_event_->Fill();
1259   
1260   return 0;
1261 }
1262   
1263 void InttAna::readRawHit(PHCompositeNode *topNode)
1264 {
1265   TrkrHitSetContainer *m_hits = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
1266   if (!m_hits)
1267   {
1268     if( Verbosity() > 1 )
1269       {
1270     cout << PHWHERE << "ERROR: Can't find node TRKR_HITSET" << endl;
1271       }
1272     
1273     return;
1274   }
1275 
1276   // uint8_t getLadderZId(TrkrDefs::hitsetkey key);
1277   // uint8_t getLadderPhiId(TrkrDefs::hitsetkey key);
1278   // int getTimeBucketId(TrkrDefs::hitsetkey key);
1279   //
1280   // uint16_t getCol(TrkrDefs::hitkey key); // z-strip = offline chip
1281   // uint16_t getRow(TrkrDefs::hitkey key); // y-strip = offline channel
1282 
1283   // loop over the InttHitSet objects
1284   TrkrHitSetContainer::ConstRange hitsetrange =
1285       m_hits->getHitSets(TrkrDefs::TrkrId::inttId);
1286   for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first;
1287        hitsetitr != hitsetrange.second;
1288        ++hitsetitr)
1289   {
1290     // Each hitset contains only hits that are clusterizable - i.e. belong to a single sensor
1291     TrkrHitSet *hitset = hitsetitr->second;
1292 
1293     if (Verbosity() > 1)
1294       cout << "InttClusterizer found hitsetkey " << hitsetitr->first << endl;
1295     if (Verbosity() > 2)
1296       hitset->identify();
1297 
1298     // we have a single hitset, get the info that identifies the sensor
1299     int layer = TrkrDefs::getLayer(hitsetitr->first);
1300     int ladder_z_index = InttDefs::getLadderZId(hitsetitr->first);
1301     int ladder_phi_index = InttDefs::getLadderPhiId(hitsetitr->first);
1302     int ladder_bco_index = InttDefs::getTimeBucketId(hitsetitr->first);
1303 
1304     if( Verbosity() > 2 )
1305       {
1306     cout << "layer " << layer << " " << ladder_z_index << " " << ladder_phi_index << " " << ladder_bco_index << endl;
1307       }
1308     
1309     //// we will need the geometry object for this layer to get the global position
1310     // CylinderGeomIntt* geom = dynamic_cast<CylinderGeomIntt*>(geom_container->GetLayerGeom(layer));
1311     // float pitch = geom->get_strip_y_spacing();
1312     // float length = geom->get_strip_z_spacing();
1313 
1314     // fill a vector of hits to make things easier - gets every hit in the hitset
1315     std::vector<std::pair<TrkrDefs::hitkey, TrkrHit *>> hitvec;
1316     TrkrHitSet::ConstRange hitrangei = hitset->getHits();
1317     for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
1318          hitr != hitrangei.second;
1319          ++hitr)
1320     {
1321       hitvec.push_back(make_pair(hitr->first, hitr->second));
1322       int chip = InttDefs::getCol(hitr->first);
1323       int chan = InttDefs::getRow(hitr->first);
1324       int adc = (hitr->second)->getAdc();
1325       if( Verbosity() > 2 )
1326     {
1327       cout << "     hit : " << chip << " " << chan << " " << adc << endl;
1328     }
1329       
1330     }
1331     if( Verbosity() > 2 )
1332       {
1333     cout << "hitvec.size(): " << hitvec.size() << endl;
1334       }
1335     
1336   }
1337 }
1338 
1339 void InttAna::getHEPMCTruth(PHCompositeNode *topNode)
1340 {
1341   if( Verbosity() > 2 )
1342     {
1343       cout << "getHEPMCTruth!!!" << endl;
1344     }
1345   
1346   /// Get the node from the node tree
1347   hepmceventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
1348 
1349   /// If the node was not properly put on the tree, return
1350   if (!hepmceventmap)
1351   {
1352     if( Verbosity() > 2 )
1353       {
1354     std::cout << PHWHERE
1355           << "HEPMC event map node is missing, can't collected HEPMC truth particles"
1356           << std::endl;
1357       }
1358     
1359     return;
1360   }
1361 
1362   /// Could have some print statements for debugging with verbosity
1363   if (Verbosity() > 1)
1364   {
1365     std::cout << "Getting HEPMC truth particles " << std::endl;
1366   }
1367 
1368   /// You can iterate over the number of events in a hepmc event
1369   /// for pile up events where you have multiple hard scatterings per bunch crossing
1370   for (PHHepMCGenEventMap::ConstIter eventIter = hepmceventmap->begin();
1371        eventIter != hepmceventmap->end();
1372        ++eventIter)
1373   {
1374     /// Get the event
1375     PHHepMCGenEvent *hepmcevent = eventIter->second;
1376 
1377     if (hepmcevent)
1378     {
1379       /// Get the event characteristics, inherited from HepMC classes
1380       HepMC::GenEvent *truthevent = hepmcevent->getEvent();
1381       if (!truthevent)
1382       {
1383         if( Verbosity() > 1 )
1384       {
1385         std::cout << PHWHERE
1386                   << "no evt pointer under phhepmvgeneventmap found "
1387                   << std::endl;
1388       }
1389     
1390         return;
1391       }
1392 
1393       /// Get the parton info
1394       HepMC::PdfInfo *pdfinfo = truthevent->pdf_info();
1395 
1396       /// Get the parton info as determined from HEPMC
1397       m_partid1 = pdfinfo->id1();
1398       m_partid2 = pdfinfo->id2();
1399       m_x1 = pdfinfo->x1();
1400       m_x2 = pdfinfo->x2();
1401 
1402       /// Are there multiple partonic intercations in a p+p event
1403       m_mpi = truthevent->mpi();
1404 
1405       /// Get the PYTHIA signal process id identifying the 2-to-2 hard process
1406       m_process_id = truthevent->signal_process_id();
1407 
1408       if (Verbosity() > 2)
1409       {
1410         std::cout << " Iterating over an event" << std::endl;
1411       }
1412 
1413       /// Loop over all the truth particles and get their information
1414       for (HepMC::GenEvent::particle_const_iterator iter = truthevent->particles_begin();
1415            iter != truthevent->particles_end();
1416            ++iter)
1417       {
1418         /// Get each pythia particle characteristics
1419         m_truthenergy = (*iter)->momentum().e();
1420         m_truthpid = (*iter)->pdg_id();
1421 
1422         // TODO: truth information
1423         m_xvtx = xvtx_sim;
1424         m_yvtx = yvtx_sim;
1425         m_zvtx = zvtx_sim;
1426         m_trutheta = (*iter)->momentum().pseudoRapidity();
1427         m_truththeta = 2 * atan(exp(-m_trutheta));
1428         // h_eta->Fill(m_trutheta);
1429         m_truthphi = (*iter)->momentum().phi();
1430         m_truthpx = (*iter)->momentum().px();
1431         m_truthpy = (*iter)->momentum().py();
1432         m_truthpz = (*iter)->momentum().pz();
1433 
1434         m_truthpt = sqrt(m_truthpx * m_truthpx + m_truthpy * m_truthpy);
1435         m_status = (*iter)->status();
1436 
1437         // m_truthphi = atan2(m_truthpy,m_truthpx);
1438 
1439         // cout << "status: " << (*iter)->status() << " " << m_truthpid << endl;
1440 
1441         h_phi->Fill(m_truthphi);
1442         h_theta->Fill(m_truththeta);
1443 
1444         /*std::cout << "m_trutheta = " << m_trutheta << "  m_truthphi = " << m_truthphi << std::endl;*/
1445 
1446         /// Fill the truth tree
1447         m_hepmctree->Fill();
1448         m_numparticlesinevent++;
1449       }
1450 
1451       // for (HepMC::GenEventVertexRange::vertex_const_iterator viter = truthevent->vertices_begin();
1452       //      viter != truthevent->vertices_end();
1453       //      ++viter)
1454       // {
1455       //   /// Get each pythia particle characteristics
1456       //   m_vertex = (*viter)->vertex_range();
1457       //   cout<<m_vertex<<endl;
1458 
1459       //   /// Fill the truth tree
1460       //   m_hepmctree->Fill();
1461       // }
1462     }
1463     // cout << "truth event = " << m_evt << endl;
1464     m_evt++;
1465   }
1466 }
1467 
1468 void InttAna::getPHG4Particle(PHCompositeNode * topNode)
1469 {
1470   if( Verbosity() > 1 )
1471     {
1472       cout << "getPHG4Particle!!!" << endl;
1473     }
1474   
1475   /// Get the node from the node tree
1476   phg4inevent = findNode::getClass<PHG4InEvent>(topNode, "PHG4INEVENT");
1477 
1478   /// If the node was not properly put on the tree, return
1479   if (!phg4inevent)
1480   {
1481     if( Verbosity() > 2 )
1482       {
1483     std::cout << PHWHERE
1484               << "PHG4InEvent node is missing, can't collected PHG4 truth particles"
1485               << std::endl;
1486       }
1487     
1488     return;
1489   }
1490 
1491   /// Could have some print statements for debugging with verbosity
1492   if (Verbosity() > 1)
1493   {
1494     std::cout << "Getting PHG4InEvent truth particles " << std::endl;
1495   }
1496 
1497   /// You can iterate over the number of events in a hepmc event
1498   /// for pile up events where you have multiple hard scatterings per bunch crossing
1499   std::map<int, PHG4VtxPoint *>::const_iterator vtxiter;
1500   std::multimap<int, PHG4Particle *>::const_iterator particle_iter;
1501   std::pair<std::map<int, PHG4VtxPoint *>::const_iterator, std::map<int, PHG4VtxPoint *>::const_iterator> vtxbegin_end = phg4inevent->GetVertices();
1502 
1503   for (vtxiter = vtxbegin_end.first; vtxiter != vtxbegin_end.second; ++vtxiter)
1504   {
1505     //*fout << "vtx number: " << vtxiter->first << std::endl;
1506     //(*vtxiter->second).identify(*fout);
1507     std::pair<std::multimap<int, PHG4Particle *>::const_iterator,
1508               std::multimap<int, PHG4Particle *>::const_iterator>
1509         particlebegin_end = phg4inevent->GetParticles(vtxiter->first);
1510 
1511     PHG4VtxPoint *vtx = vtxiter->second;
1512     m_xvtx = vtx->get_x();
1513     m_yvtx = vtx->get_y();
1514     m_zvtx = vtx->get_z();
1515     // virtual double get_t() const { return std::numeric_limits<double>::quiet_NaN(); }
1516 
1517     for (particle_iter = particlebegin_end.first; particle_iter != particlebegin_end.second; ++particle_iter)
1518     {
1519       //      (particle_iter->second)->identify(*fout);
1520       PHG4Particle *part = particle_iter->second;
1521 
1522       m_truthenergy = part->get_e();
1523       m_truthpid = part->get_pid();
1524 
1525       m_truthpx = part->get_px();
1526       m_truthpy = part->get_py();
1527       m_truthpz = part->get_pz();
1528       m_truthphi = atan2(m_truthpy, m_truthpx);
1529       m_truthpt = sqrt(m_truthpx * m_truthpx + m_truthpy * m_truthpy);
1530 
1531       m_truththeta = atan2(m_truthpt, m_truthpz);
1532       m_trutheta = -log(tan(0.5 * m_truththeta));
1533       m_status = 1;
1534 
1535       // part->get_pid() const override { return fpid; }
1536       // part->get_name() const override { return fname; }
1537       // part->get_e() const override { return fe; }
1538       // part->get_track_id();
1539       // part->get_vtx_id() const override { return vtxid; }
1540       // part->get_parent_id() const override { return parentid; }
1541       // int get_primary_id() const override { return primaryid; }
1542 
1543       /// Fill the truth tree
1544       h_phi->Fill(m_truthphi);
1545       h_theta->Fill(m_truththeta);
1546 
1547       m_hepmctree->Fill();
1548       m_numparticlesinevent++;
1549     }
1550   }
1551 
1552   if( Verbosity() > 2 )
1553     {
1554       cout << "truth event = " << m_evt << endl;
1555     }
1556   
1557   m_evt++;
1558 
1559   /*
1560     for (PHHepMCGenEventMap::ConstIter eventIter = hepmceventmap->begin();
1561          eventIter != hepmceventmap->end();
1562          ++eventIter)
1563     {
1564       /// Get the event
1565       PHHepMCGenEvent *hepmcevent = eventIter->second;
1566 
1567       if (hepmcevent)
1568       {
1569         /// Get the parton info
1570         HepMC::PdfInfo *pdfinfo = truthevent->pdf_info();
1571 
1572         /// Get the parton info as determined from HEPMC
1573         m_partid1 = pdfinfo->id1();
1574         m_partid2 = pdfinfo->id2();
1575         m_x1 = pdfinfo->x1();
1576         m_x2 = pdfinfo->x2();
1577 
1578         /// Are there multiple partonic intercations in a p+p event
1579         m_mpi = truthevent->mpi();
1580 
1581         /// Get the PYTHIA signal process id identifying the 2-to-2 hard process
1582         m_process_id = truthevent->signal_process_id();
1583       }
1584       // cout << "truth event = " << m_evt << endl;
1585       m_evt++;
1586     }
1587   */
1588 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!