event_display/src/InttAna_emc.cc

0001 #include "InttAna.h"
0002
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
0061 using namespace std;
0062
0063 /**
0064  * This class demonstrates the construction and use of an analysis module
0065  * within the sPHENIX Fun4All framework. It is intended to show how to
0066  * obtain physics objects from the analysis tree, and then write them out
0067  * to a ROOT tree and file for further offline analysis.
0068  */
0069
0070 /**
0071  * Constructor of module
0072  */
0073 InttAna::InttAna(const std::string &name, const std::string &filename)
0074     : SubsysReco(name), _rawModule(nullptr), fname_(filename), anafile_(nullptr), h_zvtxseed_(nullptr)
0075 {
0076   xbeam_ = ybeam_ = 0.0;
0077 }
0078
0079 /**
0080  * Destructor of module
0081  */
0082 InttAna::~InttAna()
0083 {
0084 }
0085
0086 /**
0087  * Initialize the module and prepare looping over events
0088  */
0089 int InttAna::Init(PHCompositeNode * /*topNode*/)
0090 {
0091   if (Verbosity() > 5)
0092   {
0093     std::cout << "Beginning Init in InttAna" << std::endl;
0094   }
0095
0096   anafile_ = new TFile(fname_.c_str(), "recreate");
0097   h_dca2d_zero = new TH1F("h_dca2d_zero", "DCA2D to 0", 500, -10, 10);
0098   h2_dca2d_zero = new TH2F("h2_dca2d_zero", "DCA2D to 0 vs nclus", 500, -10, 10, 100, 0, 10000);
0099   h2_dca2d_len = new TH2F("h2_dca2d_len", "DCA2D to 0 vs len", 500, -10, 10, 500, -10, 10);
0100   h_zvtx = new TH1F("h_zvtx", "Zvertex", 400, -40, 40);
0101   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");
0102   h_ntp_emcclus = new TNtuple("ntp_emcclus", "EMC Cluster Ntuple", "nemc:evt_emc:x_emc:y_emc:z_emc:e:ecore:size_emc");
0103
0104   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");
0105   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");
0106   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");
0107
0108   h_zvtxseed_ = new TH1F("h_zvtxseed", "Zvertex Seed histogram", 200, -50, 50);
0109
0110   h_t_evt_bco = new TTree("t_evt_bco", "Event Tree for BCO");
0111   h_t_evt_bco->Branch("evt_gl1", &(m_evtbcoinfo.evt_gl1), "evt_gl1/I");
0112   h_t_evt_bco->Branch("evt_intt", &(m_evtbcoinfo.evt_intt), "evt_intt/i");
0113   h_t_evt_bco->Branch("evt_mbd", &(m_evtbcoinfo.evt_mbd), "evt_mbd/i");
0114   h_t_evt_bco->Branch("bco_gl1", &(m_evtbcoinfo.bco_gl1), "bco_gl1/l");
0115   h_t_evt_bco->Branch("bco_intt", &(m_evtbcoinfo.bco_intt), "bco_intt/l");
0116   h_t_evt_bco->Branch("bco_mbd", &(m_evtbcoinfo.bco_mbd), "bco_mbd/l");
0117   h_t_evt_bco->Branch("pevt_gl1", &(m_evtbcoinfo_prev.evt_gl1), "pevt_gl1/I");
0118   h_t_evt_bco->Branch("pevt_intt", &(m_evtbcoinfo_prev.evt_intt), "pevt_intt/i");
0119   h_t_evt_bco->Branch("pevt_mbd", &(m_evtbcoinfo_prev.evt_mbd), "pevt_mbd/i");
0120   h_t_evt_bco->Branch("pbco_gl1", &(m_evtbcoinfo_prev.bco_gl1), "pbco_gl1/l");
0121   h_t_evt_bco->Branch("pbco_intt", &(m_evtbcoinfo_prev.bco_intt), "pbco_intt/l");
0122   h_t_evt_bco->Branch("pbco_mbd", &(m_evtbcoinfo_prev.bco_mbd), "pbco_mbd/l");
0123
0124   return 0;
0125 }
0126
0127 int InttAna::InitRun(PHCompositeNode * /*topNode*/)
0128 {
0129   cout << "InttAna::InitRun beamcenter " << xbeam_ << " " << ybeam_ << endl;
0130
0131   return 0;
0132 }
0133
0134 /**
0135  * Main workhorse function where each event is looped over and
0136  * data from each event is collected from the node tree for analysis
0137  */
0138 int InttAna::process_event(PHCompositeNode *topNode)
0139 {
0140   static int ievt = 0;
0141   cout << "InttEvt::process evt : " << ievt++ << endl;
0142
0143   if (Verbosity() > 5)
0144   {
0145     std::cout << "Beginning process_event in AnaTutorial" << std::endl;
0146   }
0147
0148   // readRawHit(topNode);
0149
0150   ActsGeometry *m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0151   if (!m_tGeometry)
0152   {
0153     std::cout << PHWHERE << "No ActsGeometry on node tree. Bailing." << std::endl;
0154     return Fun4AllReturnCodes::ABORTEVENT;
0155   }
0156
0157   TrkrClusterContainer *m_clusterMap =
0158       findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0159   if (!m_clusterMap)
0160   {
0161     cout << PHWHERE << "TrkrClusterContainer node is missing." << endl;
0162     return Fun4AllReturnCodes::ABORTEVENT;
0163   }
0164
0165   InttEventInfo *inttevthead = findNode::getClass<InttEventInfo>(topNode, "INTTEVENTHEADER");
0166
0167   InttVertex3DMap *inttvertexmap = findNode::getClass<InttVertex3DMap>(topNode, "InttVertexMap");
0168
0169   double vtx_sim[3]{-9999, -9999, -9999};
0170
0171   std::map<GlobalVertex::VTXTYPE, std::string> sourcemap{
0172       {GlobalVertex::VTXTYPE::UNDEFINED, "UNDEFINED"},
0173       {GlobalVertex::VTXTYPE::TRUTH, "TRUTH"},
0174       {GlobalVertex::VTXTYPE::SMEARED, "SMEARED"},
0175       {GlobalVertex::VTXTYPE::MBD, "MBD"},
0176       {GlobalVertex::VTXTYPE::SVTX, "SVTX"},
0177       {GlobalVertex::VTXTYPE::SVTX_MBD, "SVTX_MBD"}};
0178   //   enum VTXTYPE
0179   std::string s_vtxid = "NULL";
0180
0181   GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0182   if (vertexmap)
0183   {
0184     if (!vertexmap->empty())
0185     {
0186
0187       for (auto vertex : *vertexmap)
0188       {
0189         std::cout << "map entry" << std::endl;
0190         std::map<GlobalVertex::VTXTYPE, std::string>::iterator itr_src;
0191         itr_src = sourcemap.find((GlobalVertex::VTXTYPE)vertex.first);
0192         if (itr_src != sourcemap.end())
0193         {
0194           s_vtxid = itr_src->second;
0195         }
0196
0197         if (vertex.first == GlobalVertex::VTXTYPE::TRUTH)
0198         {
0199           vtx_sim[0] = vertex.second->get_x();
0200           vtx_sim[1] = vertex.second->get_y();
0201           vtx_sim[2] = vertex.second->get_z();
0202         }
0203
0204         std::cout << std::endl
0205                   << "A " << vertex.second->get_x()
0206                   << " " << vertex.second->get_y()
0207                   << " " << vertex.second->get_z()
0208                   << " " << vertex.second->get_id()
0209                   << " " << s_vtxid << std::endl;
0210       }
0211     }
0212   }
0213
0214   /*
0215     MbdVertexMap *mbdvertexmap = findNode::getClass<MbdVertexMap>(topNode, "MbdVertexMap");
0216     MbdVertex* mbdvtx = nullptr;
0217     if (mbdvertexmap){
0218       cout<<"MbdVertexMap"<<endl;
0219       if (!mbdvertexmap->empty()){
0220
0221         for(auto mbdvertex : *mbdvertexmap)
0222         {
0223           std::cout<<"mbdmap entry"<<std::endl;
0224
0225           std::cout<<std::endl<<"Mbdz"
0226                               <<" " <<mbdvertex.second->get_z()
0227                               <<" " <<mbdvertex.second->get_id()<<endl;
0228           mbdvtx = mbdvertex;
0229         }
0230       }
0231       else {
0232         cout<<"MbdVertexMap empty"<<endl;
0233       }
0234     }
0235   */
0236
0237   SvtxVertexMap *svtxvertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0238   SvtxVertex *svtxvertex = nullptr;
0239   if (svtxvertexmap)
0240   {
0241     cout << "svtxvertex: size : " << svtxvertexmap->size() << endl;
0242     for (SvtxVertexMap::ConstIter iter = svtxvertexmap->begin(); iter != svtxvertexmap->end(); ++iter)
0243     {
0244       svtxvertex = iter->second;
0245       std::cout << std::endl
0246                 << "SvtxVertex"
0247                 << " " << svtxvertex->get_x()
0248                 << " " << svtxvertex->get_y()
0249                 << " " << svtxvertex->get_z()
0250                 << " " << svtxvertex->get_id() << endl;
0251     }
0252   }
0253   else
0254   {
0255     cout << "SvtxVertexMap is not available" << endl;
0256   }
0257
0258   MbdOut *mbdout = findNode::getClass<MbdOut>(topNode, "MbdOut");
0259   if (!mbdout)
0260   {
0261     cout << "No MbdOut" << endl;
0262   }
0263   double mbdqs = (mbdout != nullptr) ? mbdout->get_q(0) : -9999;
0264   double mbdqn = (mbdout != nullptr) ? mbdout->get_q(1) : -9999;
0265   //  double mbdt0 = (mbdout!=nullptr) ?  mbdout->get_t0() : -9999;
0266   double mbdz = (mbdout != nullptr) ? mbdout->get_zvtx() : -9999;
0267
0268   Gl1RawHit *gl1raw = findNode::getClass<Gl1RawHit>(topNode, "GL1RAWHIT");
0269   if (!gl1raw)
0270   {
0271     cout << "No Gl1Raw" << endl;
0272   }
0273
0274   InttRawHitContainer *inttrawmap = findNode::getClass<InttRawHitContainer>(topNode, "INTTRAWHIT");
0275   InttRawHit *inttraw = (inttrawmap != nullptr && inttrawmap->get_nhits() > 0) ? inttrawmap->get_hit(0) : nullptr;
0276   /////////////
0277   //  InttEvent from RawModule (not standard way)
0278   InttEvent *inttEvt = nullptr;
0279   if (_rawModule)
0280   {
0281     inttEvt = _rawModule->getInttEvent();
0282   }
0283   uint64_t bco = 0;
0284   int evtseq = -1;
0285   if (inttEvt != NULL)
0286   {
0287     bco = inttEvt->bco;
0288     evtseq = inttEvt->evtSeq;
0289   }
0290   if (inttevthead)
0291   {
0292     bco = inttevthead->get_bco_full();
0293     // evtseq = inttEvt->evtSeq;
0294   }
0295   cout << "bco 0x" << hex << bco << dec << endl;
0296
0297   if (gl1raw)
0298   {
0299     bco = gl1raw->get_bco();
0300     evtseq = gl1raw->getEvtSequence();
0301   }
0302
0303   double vertex[10][3]{};
0304
0305   InttVertex3D *zvtxobj = NULL;
0306   if (inttvertexmap)
0307   {
0308     int ivtx = 0;
0309     cout << "vertex map size : " << inttvertexmap->size() << endl;
0310     InttVertex3DMap::ConstIter biter_beg = inttvertexmap->begin();
0311     InttVertex3DMap::ConstIter biter_end = inttvertexmap->end();
0312     // cout<<"vertex map size : after size "<<endl;
0313     // inttvertexmap->identify();
0314
0315     for (InttVertex3DMap::ConstIter biter = biter_beg; biter != biter_end; ++biter)
0316     {
0317       InttVertex3D *vtxobj = biter->second;
0318       if (vtxobj)
0319       {
0320         cout << "vtxobj" << ivtx << endl;
0321         vtxobj->identify();
0322         vertex[ivtx][0] = vtxobj->get_x();
0323         vertex[ivtx][1] = vtxobj->get_y();
0324         vertex[ivtx][2] = vtxobj->get_z();
0325         if (ivtx == 2)
0326         {
0327           zvtxobj = vtxobj;
0328         }
0329       }
0330       else
0331       {
0332         cout << "no vtx object : " << ivtx << endl;
0333       }
0334       ivtx++;
0335       if (ivtx >= 10)
0336       {
0337         cout << "n_vertex exceeded : " << ivtx << endl;
0338         break;
0339       }
0340     }
0341   }
0342
0343   ////
0344   // fill evtbco_info
0345   m_evtbcoinfo.clear();
0346   m_evtbcoinfo.evt_gl1 = (gl1raw != nullptr) ? gl1raw->getEvtSequence() : -1;
0347   m_evtbcoinfo.evt_intt = (inttraw != nullptr) ? inttraw->get_event_counter() : -1;
0348   m_evtbcoinfo.evt_mbd = (mbdout != nullptr) ? mbdout->get_evt() : -1;
0349   m_evtbcoinfo.bco_gl1 = (gl1raw != nullptr) ? gl1raw->get_bco() : 0;
0350   m_evtbcoinfo.bco_intt = (inttevthead != nullptr) ? inttevthead->get_bco_full() : 0;
0351   m_evtbcoinfo.bco_mbd = (mbdout != nullptr) ? mbdout->get_femclock() : 0;
0352
0353   cout << "event : " << m_evtbcoinfo.evt_gl1 << " "
0354        << m_evtbcoinfo.evt_intt << " "
0355        << m_evtbcoinfo.evt_mbd << " :  "
0356        << hex << m_evtbcoinfo.bco_gl1 << dec << " "
0357        << hex << m_evtbcoinfo.bco_intt << dec << " "
0358        << hex << m_evtbcoinfo.bco_mbd << dec << " : "
0359        << hex << m_evtbcoinfo.bco_gl1 - m_evtbcoinfo_prev.bco_gl1 << dec << " "
0360        << hex << m_evtbcoinfo.bco_intt - m_evtbcoinfo_prev.bco_intt << dec << " "
0361        << hex << m_evtbcoinfo.bco_mbd - m_evtbcoinfo_prev.bco_mbd << dec << endl;
0362
0363   /////////////
0364   int nclusadd = 0, nclusadd2 = 0;
0365   int nclus_inner = 0, nclus_outer = 0;
0366   for (unsigned int inttlayer = 0; inttlayer < 4; inttlayer++)
0367   {
0368     for (const auto &hitsetkey : m_clusterMap->getHitSetKeys(TrkrDefs::TrkrId::inttId, inttlayer + 3))
0369     {
0370       auto range = m_clusterMap->getClusters(hitsetkey);
0371
0372       for (auto clusIter = range.first; clusIter != range.second; ++clusIter)
0373       {
0374         nclusadd++;
0375         const auto cluster = clusIter->second;
0376         if (cluster->getAdc() > 40)
0377           nclusadd2++;
0378
0379         if (inttlayer < 2)
0380           nclus_inner++;
0381         else
0382           nclus_outer++;
0383       }
0384     }
0385   }
0386   /////////////
0387   // MVTX cluster
0388   std::set<TrkrDefs::TrkrId> detectors;
0389   detectors.insert(TrkrDefs::TrkrId::mvtxId);
0390   int nclusmvtx[3] = {0, 0, 0};
0391   for (const auto &det : detectors)
0392   {
0393     for (const auto &layer : {0, 1, 2})
0394     {
0395       for (const auto &hitsetkey : m_clusterMap->getHitSetKeys(det, layer))
0396       {
0397         auto range = m_clusterMap->getClusters(hitsetkey);
0398         for (auto citer = range.first; citer != range.second; ++citer)
0399         {
0400           // const auto ckey = citer->first;
0401           // const auto cluster = citer->second;
0402           // const auto global = m_tGeometry->getGlobalPosition(ckey, cluster);
0403           nclusmvtx[layer]++;
0404         }
0405       }
0406     }
0407   }
0408
0409   /////////////
0410   // SvtxTrack
0411   SvtxTrackMap *svtxtrackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0412
0413   /////////////
0414   // EMC Cluster
0415   RawClusterContainer *clusterContainer = findNode::getClass<RawClusterContainer>(topNode, "CLUSTERINFO_POS_COR_CEMC");
0416   // RawClusterContainer *clusterContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_POS_COR_CEMC");
0417   // RawClusterContainer *clusterContainer = findNode::getClass<RawClusterContainer>(topNode,"CLUSTER_CEMC");
0418   int nemc = -9999, nemc1 = -9999;
0419   if (!clusterContainer)
0420   {
0421     std::cout << PHWHERE << "InttAna::process_event - Fatal Error - CLUSTER_POS_COR_CEMC node is missing. " << std::endl;
0422   }
0423   else
0424   {
0425     float ntpval_emccls[20];
0426     // This is how we iterate over items in the container.
0427     RawClusterContainer::ConstRange clusterEnd = clusterContainer->getClusters();
0428     RawClusterContainer::ConstIterator clusterIter;
0429
0430     nemc = clusterContainer->size();
0431     nemc1 = 0;
0432
0433     for (clusterIter = clusterEnd.first; clusterIter != clusterEnd.second; clusterIter++)
0434     {
0435       RawCluster *recoCluster = clusterIter->second;
0436
0437       ntpval_emccls[0] = nemc;
0438       ntpval_emccls[1] = evtseq;
0439       ntpval_emccls[2] = recoCluster->get_x();
0440       ntpval_emccls[3] = recoCluster->get_y();
0441       ntpval_emccls[4] = recoCluster->get_z();
0442       ntpval_emccls[5] = recoCluster->get_energy();
0443       ntpval_emccls[6] = recoCluster->get_ecore();
0444       ntpval_emccls[7] = recoCluster->getNTowers();
0445
0446       h_ntp_emcclus->Fill(ntpval_emccls);
0447
0448       if (recoCluster->get_ecore() > 0.15)
0449         nemc1++;
0450     }
0451   }
0452
0453   /////////////
0454   static int evtCount = 0;
0455
0456   cout << "evtCount : " << evtCount << " " << evtseq << " " << hex << bco << dec << endl;
0457
0458   struct ClustInfo
0459   {
0460     int layer;
0461     Acts::Vector3 pos;
0462   };
0463
0464   float ntpval[20];
0465   int nCluster = 0;
0466   bool exceedNwrite = false;
0467   // std::vector<Acts::Vector3> clusters[2]; // inner=0; outer=1
0468   std::vector<ClustInfo> clusters[2]; // inner=0; outer=1
0469   for (unsigned int inttlayer = 0; inttlayer < 4; inttlayer++)
0470   {
0471     int layer = (inttlayer < 2 ? 0 : 1);
0472     for (const auto &hitsetkey : m_clusterMap->getHitSetKeys(TrkrDefs::TrkrId::inttId, inttlayer + 3))
0473     {
0474       auto range = m_clusterMap->getClusters(hitsetkey);
0475
0476       for (auto clusIter = range.first; clusIter != range.second; ++clusIter)
0477       {
0478         const auto cluskey = clusIter->first;
0479         const auto cluster = clusIter->second;
0480
0481         int ladder_z = InttDefs::getLadderZId(cluskey);
0482         int ladder_phi = InttDefs::getLadderPhiId(cluskey);
0483         int size = cluster->getSize();
0484
0485         // cout    <<cluster->getPosition(0)<<" "
0486         //         <<cluster->getPosition(1)<<" : "
0487         //         <<cluster->getAdc()<<" "<<size<<" "<<inttlayer<<" "<<ladder_z<<" "<<ladder_phi<<endl;
0488
0489         const auto globalPos = m_tGeometry->getGlobalPosition(cluskey, cluster);
0490
0491         if (nCluster < 5)
0492         {
0493           cout << "xyz : " << globalPos.x() << " " << globalPos.y() << " " << globalPos.z() << " :  "
0494                << cluster->getPosition(0) << " "
0495                << cluster->getPosition(1) << " : "
0496                << cluster->getAdc() << " " << size << " " << inttlayer << " " << ladder_z << " " << ladder_phi << endl;
0497         }
0498         else
0499         {
0500           if (!exceedNwrite)
0501           {
0502             cout << " exceed : ncluster limit.  no more cluster xyz printed" << endl;
0503             exceedNwrite = true;
0504           }
0505         }
0506
0507         ClustInfo info;
0508         info.layer = inttlayer;
0509         info.pos = globalPos;
0510
0511         // clusters[layer].push_back(globalPos);
0512         clusters[layer].push_back(info);
0513         nCluster++;
0514
0515         ntpval[0] = nclusadd;          // bco_full
0516         ntpval[1] = nclusadd2;         // bco_full
0517         ntpval[2] = bco;               // bco_full
0518         ntpval[3] = evtseq;            // evtCount;
0519         ntpval[4] = size;              // size
0520         ntpval[5] = cluster->getAdc(); // size
0521         ntpval[6] = globalPos.x();
0522         ntpval[7] = globalPos.y();
0523         ntpval[8] = globalPos.z();
0524         ntpval[9] = inttlayer;
0525         ntpval[10] = ladder_phi;
0526         ntpval[11] = ladder_z;
0527         ntpval[12] = cluster->getPosition(0);
0528         ntpval[13] = cluster->getPosition(1);
0529         ntpval[14] = cluster->getPhiSize();
0530         ntpval[15] = cluster->getZSize();
0531         ntpval[16] = vertex[2][2]; // zvtx;
0532         h_ntp_clus->Fill(ntpval);
0533         // = new TNtuple("ntp_clus", "Cluster Ntuple", "bco_full:evt:size:adc:x:y:z:lay:lad:sen");
0534       }
0535     }
0536   }
0537   cout << "nCluster : " << nCluster << endl;
0538
0539   double zvtx = -9999.;
0540
0541   struct cluspair
0542   {
0543     float p1_ang;
0544     float p2_ang;
0545     float p1_r;
0546     float p2_r;
0547     float p1_z;
0548     float p2_z;
0549     float dca_p0;
0550     float len_p0;
0551   };
0552
0553   vector<cluspair> vcluspair;
0554
0555   if (nCluster < 300)
0556   // if(nCluster<500)
0557   {
0558     float ntpval2[20];
0559     Acts::Vector3 beamspot = Acts::Vector3(xbeam_, ybeam_, 0);
0560     vector<double> vz_array;
0561     for (auto c1 = clusters[0].begin(); c1 != clusters[0].end(); ++c1) // inner
0562     {
0563       for (auto c2 = clusters[1].begin(); c2 != clusters[1].end(); ++c2) // outer
0564       {
0565         Acts::Vector3 p1 = c1->pos - beamspot;
0566         Acts::Vector3 p2 = c2->pos - beamspot;
0567         Acts::Vector3 u = p2 - p1;
0568         double unorm = sqrt(u.x() * u.x() + u.y() * u.y());
0569         if (unorm < 0.00001)
0570           continue;
0571
0572         // skip bad compbination
0573         double p1_ang = atan2(p1.y(), p1.x());
0574         double p2_ang = atan2(p2.y(), p2.x());
0575         double d_ang = p2_ang - p1_ang;
0576
0577         // unit vector in 2D
0578         double ux = u.x() / unorm;
0579         double uy = u.y() / unorm;
0580         double uz = u.z() / unorm; // same normalization factor(2D) is multiplied
0581
0582         // Acts::Vector3 p0   = beamspot - p1;
0583         Acts::Vector3 p0 = Acts::Vector3(0, 0, 0) - p1; // p1, p2 are already shifted by beam center
0584
0585         double dca_p0 = p0.x() * uy - p0.y() * ux; // cross product of p0 x u
0586         double len_p0 = p0.x() * ux + p0.y() * uy; // dot product of p0 . u
0587
0588         // beam center in X-Y plane
0589         double vx = len_p0 * ux + p1.x();
0590         double vy = len_p0 * uy + p1.y();
0591
0592         double vz = len_p0 * uz + p1.z();
0593
0594         double p1_r = sqrt(p1.y() * p1.y() + p1.x() * p1.x());
0595         double p2_r = sqrt(p2.y() * p2.y() + p2.x() * p2.x());
0596
0597         // if(unorm>4.5||d_ang<-0.005*3.1415||0.005*3.1415<d_ang)
0598         // if(unorm>4.5||d_ang<-0.25*3.1415||0.25*3.1415<d_ang)
0599         if (fabs(d_ang) < 0.05 && fabs(dca_p0) < 0.5)
0600         {
0601           h_zvtxseed_->Fill(vz);
0602           vz_array.push_back(vz);
0603
0604           cluspair clspair;
0605           clspair.p1_ang = p1_ang;
0606           clspair.p2_ang = p2_ang;
0607           clspair.p1_r = p1_r;
0608           clspair.p2_r = p2_r;
0609           clspair.p1_z = p1.z();
0610           clspair.p2_z = p2.z();
0611           clspair.dca_p0 = dca_p0;
0612           clspair.len_p0 = len_p0;
0613           vcluspair.push_back(clspair);
0614         }
0615
0616         h_dca2d_zero->Fill(dca_p0);
0617         h2_dca2d_zero->Fill(dca_p0, nCluster);
0618         h2_dca2d_len->Fill(dca_p0, len_p0);
0619         // cout<<"dca : "<<dca_p0<<endl;
0620         //
0621         ntpval2[0] = nclusadd;
0622         ntpval2[1] = nclusadd2;
0623         ntpval2[2] = 0;
0624         ntpval2[3] = evtCount;
0625         ntpval2[4] = p1_ang;
0626         ntpval2[5] = p2_ang;
0627         ntpval2[6] = p1.z();
0628         ntpval2[7] = p2.z();
0629         ntpval2[8] = dca_p0;
0630         ntpval2[9] = len_p0;
0631         ntpval2[10] = unorm;
0632         ntpval2[11] = c1->layer;
0633         ntpval2[12] = c2->layer;
0634         ntpval2[13] = vx;
0635         ntpval2[14] = vy;
0636         ntpval2[15] = vz;
0637         ntpval2[16] = vertex[2][2]; // zvtx;
0638
0639         if (evtCount < 10000)
0640           h_ntp_cluspair->Fill(ntpval2);
0641         // h_ntp_cluspair = new TNtuple("ntp_cluspair", "Cluster Pair Ntuple", "nclus:bco_full:evt:ang1:ang2:dca2d:len:unorm");
0642       }
0643     }
0644
0645     // calculate trancated mean of DCA~Z histogram as Z-vertex position
0646
0647     if (vz_array.size() > 3)
0648     {
0649       double zbin = h_zvtxseed_->GetMaximumBin();
0650       double zcenter = h_zvtxseed_->GetBinCenter(zbin);
0651       double zmean = h_zvtxseed_->GetMean();
0652       double zrms = h_zvtxseed_->GetRMS();
0653       if (zrms < 20)
0654         zrms = 20;
0655
0656       double zmax = zcenter + zrms; // 1 sigma
0657       double zmin = zcenter - zrms; // 1 sigma
0658
0659       double zsum = 0.;
0660       int zcount = 0;
0661       for (auto iz = vz_array.begin(); iz != vz_array.end(); ++iz)
0662       {
0663         double vz = (*iz);
0664         if (zmin < vz && vz < zmax)
0665         {
0666           zsum += vz;
0667           zcount++;
0668         }
0669       }
0670       if (zcount > 0)
0671         zvtx = zsum / zcount;
0672
0673       cout << "ZVTX: " << zvtx << " " << zcenter << " " << zmean << " " << zrms << " " << zbin << endl;
0674     }
0675
0676     h_zvtx->Fill(zvtx);
0677
0678     float ntpval_emcpair[20];
0679     if (clusterContainer != nullptr)
0680     {
0681       RawClusterContainer::ConstRange clusterEnd = clusterContainer->getClusters();
0682       RawClusterContainer::ConstIterator clusterIter;
0683
0684       /////////////
0685       //  association to EMCAL
0686       for (vector<cluspair>::iterator itrcp = vcluspair.begin(); itrcp != vcluspair.end(); ++itrcp)
0687       {
0688         cluspair &clspair = *itrcp;
0689
0690         for (clusterIter = clusterEnd.first; clusterIter != clusterEnd.second; clusterIter++)
0691         {
0692           RawCluster *recoCluster = clusterIter->second;
0693
0694           float the1 = atan2(clspair.p1_r, (clspair.p1_z - vertex[2][2]));
0695           float the2 = atan2(clspair.p2_r, (clspair.p2_z - vertex[2][2]));
0696
0697           Acts::Vector3 emcpos = Acts::Vector3(recoCluster->get_x(), recoCluster->get_y(), recoCluster->get_z());
0698           Acts::Vector3 emcposc = emcpos - beamspot;
0699           float ephi = atan2(emcposc.y(), emcposc.x());
0700           float er = sqrt(emcposc.x() * emcposc.x() + emcposc.y() * emcposc.y());
0701           float ethe = atan2(er, emcposc.z() - vertex[2][2]);
0702
0703           ntpval_emcpair[0] = nclusadd;
0704           ntpval_emcpair[1] = evtseq;
0705           ntpval_emcpair[2] = clspair.p1_ang; // INTT
0706           ntpval_emcpair[3] = clspair.p2_ang;
0707           ntpval_emcpair[4] = clspair.p1_z;
0708           ntpval_emcpair[5] = clspair.p2_z;
0709           ntpval_emcpair[6] = clspair.dca_p0;
0710           ntpval_emcpair[7] = clspair.len_p0;
0711           ntpval_emcpair[8] = vertex[1][0];  // x-vertex
0712           ntpval_emcpair[9] = vertex[1][1];  // y-vertex
0713           ntpval_emcpair[10] = vertex[2][2]; // y-vertex
0714           ntpval_emcpair[11] = ephi;
0715           ntpval_emcpair[12] = recoCluster->get_z(); // EMCal
0716           ntpval_emcpair[13] = recoCluster->get_ecore();
0717           ntpval_emcpair[14] = recoCluster->getNTowers();
0718           ntpval_emcpair[15] = the1;
0719           ntpval_emcpair[16] = the2;
0720           ntpval_emcpair[17] = ethe;
0721
0722           h_ntp_emccluspair->Fill(ntpval_emcpair);
0723         }
0724       }
0725     }
0726   }
0727
0728   float ntpval3[40];
0729   ntpval3[0] = nclusadd;  // bco_full
0730   ntpval3[1] = nclusadd2; // bco_full
0731   ntpval3[2] = bco;
0732   ntpval3[3] = evtseq;
0733   ntpval3[4] = vertex[2][2]; // zvtx;
0734   ntpval3[5] = zvtx;         // zrms;
0735   ntpval3[6] = 0;            // zmean;
0736   ntpval3[7] = 0;            // zcenter;
0737   ntpval3[8] = mbdz;
0738   ntpval3[9] = mbdqn;         // mbdt.bqn;
0739   ntpval3[10] = mbdqs;        // mbdt.bqs;
0740   ntpval3[11] = 0;            // mbdt.femclk;
0741   ntpval3[12] = vertex[1][0]; // x-vertex
0742   ntpval3[13] = vertex[1][1]; // y-vertex
0743   ntpval3[14] = vertex[2][2]; // y-vertex
0744   ntpval3[15] = nclus_inner;  // x-vertex
0745   ntpval3[16] = nclus_outer;  // y-vertex
0746
0747   ntpval3[17] = zvtxobj->get_nclus();
0748   ntpval3[18] = zvtxobj->get_ntracklet();
0749   ntpval3[19] = zvtxobj->get_chi2ndf();
0750   ntpval3[20] = zvtxobj->get_width();
0751   ntpval3[21] = zvtxobj->get_ngroup();
0752   ntpval3[22] = (zvtxobj->get_good()) ? 1 : 0;
0753   ntpval3[23] = zvtxobj->get_peakratio();
0754
0755   ntpval3[24] = vtx_sim[0]; // sim x-vertex
0756   ntpval3[25] = vtx_sim[1]; // sim y-vertex
0757   ntpval3[26] = vtx_sim[2]; // sim z-vertex
0758
0759   ntpval3[27] = (svtxvertex != nullptr) ? svtxvertex->get_x() : -9999; // svtx vertex x
0760   ntpval3[28] = (svtxvertex != nullptr) ? svtxvertex->get_y() : -9999; // svtx vertex y
0761   ntpval3[29] = (svtxvertex != nullptr) ? svtxvertex->get_z() : -9999; // svtx vertex z
0762
0763   ntpval3[30] = nclusmvtx[0];
0764   ntpval3[31] = nclusmvtx[1];
0765   ntpval3[32] = nclusmvtx[2];
0766   ntpval3[33] = (svtxtrackmap != nullptr) ? svtxtrackmap->size() : -9999;
0767   ntpval3[34] = nemc;
0768   ntpval3[35] = nemc1;
0769
0770   h_ntp_evt->Fill(ntpval3);
0771   // h_ntp_evt = new TNtuple("ntp_evt", "Event Ntuple",
0772   //"nclus:nclus2:bco_full:evt:zv:
0773   // zvs:zvm:zvc:bz:bqn:
0774   // bqs:bfemclk:xvtx:yvtx:zvtx:
0775   // nclusin:nclusout:nclszv:ntrkzv:chi2ndfzv:
0776   // widthzv:ngroupzv:goodzv:peakratiozv:xvsim:
0777   // yvsim:zvsim");
0778   //
0779
0780   h_t_evt_bco->Fill();
0781
0782   m_evtbcoinfo_prev.copy(m_evtbcoinfo);
0783
0784   evtCount++;
0785
0786   return Fun4AllReturnCodes::EVENT_OK;
0787 }
0788
0789 /**
0790  * End the module and finish any data collection. Clean up any remaining
0791  * loose ends
0792  */
0793 int InttAna::End(PHCompositeNode * /*topNode*/)
0794 {
0795   if (Verbosity() > 1)
0796   {
0797     std::cout << "Ending InttAna analysis package" << std::endl;
0798   }
0799
0800   if (anafile_ != nullptr)
0801   {
0802     anafile_->Write();
0803     anafile_->Close();
0804   }
0805
0806   return 0;
0807 }
0808
0809 void InttAna::readRawHit(PHCompositeNode *topNode)
0810 {
0811   TrkrHitSetContainer *m_hits = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0812   if (!m_hits)
0813   {
0814     cout << PHWHERE << "ERROR: Can't find node TRKR_HITSET" << endl;
0815     return;
0816   }
0817
0818   // uint8_t getLadderZId(TrkrDefs::hitsetkey key);
0819   // uint8_t getLadderPhiId(TrkrDefs::hitsetkey key);
0820   // int getTimeBucketId(TrkrDefs::hitsetkey key);
0821   //
0822   // uint16_t getCol(TrkrDefs::hitkey key); // z-strip = offline chip
0823   // uint16_t getRow(TrkrDefs::hitkey key); // y-strip = offline channel
0824
0825   // loop over the InttHitSet objects
0826   TrkrHitSetContainer::ConstRange hitsetrange =
0827       m_hits->getHitSets(TrkrDefs::TrkrId::inttId);
0828   for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first;
0829        hitsetitr != hitsetrange.second;
0830        ++hitsetitr)
0831   {
0832     // Each hitset contains only hits that are clusterizable - i.e. belong to a single sensor
0833     TrkrHitSet *hitset = hitsetitr->second;
0834
0835     if (Verbosity() > 1)
0836       cout << "InttClusterizer found hitsetkey " << hitsetitr->first << endl;
0837     if (Verbosity() > 2)
0838       hitset->identify();
0839
0840     // we have a single hitset, get the info that identifies the sensor
0841     int layer = TrkrDefs::getLayer(hitsetitr->first);
0842     int ladder_z_index = InttDefs::getLadderZId(hitsetitr->first);
0843     int ladder_phi_index = InttDefs::getLadderPhiId(hitsetitr->first);
0844     int ladder_bco_index = InttDefs::getTimeBucketId(hitsetitr->first);
0845
0846     cout << "layer " << layer << " " << ladder_z_index << " " << ladder_phi_index << " " << ladder_bco_index << endl;
0847
0848     //// we will need the geometry object for this layer to get the global position
0849     // CylinderGeomIntt* geom = dynamic_cast<CylinderGeomIntt*>(geom_container->GetLayerGeom(layer));
0850     // float pitch = geom->get_strip_y_spacing();
0851     // float length = geom->get_strip_z_spacing();
0852
0853     // fill a vector of hits to make things easier - gets every hit in the hitset
0854     std::vector<std::pair<TrkrDefs::hitkey, TrkrHit *>> hitvec;
0855     TrkrHitSet::ConstRange hitrangei = hitset->getHits();
0856     for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
0857          hitr != hitrangei.second;
0858          ++hitr)
0859     {
0860       hitvec.push_back(make_pair(hitr->first, hitr->second));
0861       int chip = InttDefs::getCol(hitr->first);
0862       int chan = InttDefs::getRow(hitr->first);
0863       int adc = (hitr->second)->getAdc();
0864       cout << "     hit : " << chip << " " << chan << " " << adc << endl;
0865     }
0866     cout << "hitvec.size(): " << hitvec.size() << endl;
0867   }
0868 }