event_display/src/InttAna_MC.cc

0001 //for MC
0002 #include "InttAna.h"
0003
0004 #include <math.h>
0005
0006 /// Cluster/Calorimeter includes
0007
0008 /// Fun4All includes
0009 #include <fun4all/Fun4AllHistoManager.h>
0010 #include <fun4all/Fun4AllReturnCodes.h>
0011 #include <phool/PHCompositeNode.h>
0012 #include <phool/getClass.h>
0013
0014 // #include <trackbase/TrkrClusterv4.h>
0015 // #include <trackbase/TrkrClusterContainerv3.h>
0016 #include <trackbase/TrkrCluster.h>
0017 #include <trackbase/TrkrClusterContainer.h>
0018 #include <trackbase/TrkrHit.h>
0019 #include <trackbase/TrkrHitSet.h>
0020 #include <trackbase/TrkrHitSetContainer.h>
0021 #include <trackbase/ActsGeometry.h>
0022 #include <trackbase/InttDefs.h>
0023
0024 /// ROOT includes
0025 #include <TFile.h>
0026 #include <TH1.h>
0027 #include <TH2.h>
0028 #include <TNtuple.h>
0029 #include <TTree.h>
0030
0031 /// C++ includes
0032 #include <cassert>
0033 #include <cmath>
0034 #include <sstream>
0035 #include <string>
0036
0037 #include "InttEvent.h"
0038 #include "InttRawData.h"
0039
0040 // To get vertex
0041 #include <globalvertex/GlobalVertexMap.h>
0042 #include <globalvertex/GlobalVertex.h>
0043
0044 // TODO:
0045 /// HEPMC truth includes
0046 #pragma GCC diagnostic push
0047 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
0048 #include <HepMC/GenEvent.h>
0049 #include <HepMC/GenVertex.h>
0050 #pragma GCC diagnostic pop
0051
0052 #include <phhepmc/PHHepMCGenEvent.h>
0053 #include <phhepmc/PHHepMCGenEventMap.h>
0054
0055 #include <g4main/PHG4InEvent.h>
0056 #include <g4main/PHG4VtxPoint.h> // for PHG4Particle
0057 #include <g4main/PHG4Particle.h> // for PHG4Particle
0058
0059 using namespace std;
0060
0061 /**
0062  * This class demonstrates the construction and use of an analysis module
0063  * within the sPHENIX Fun4All framework. It is intended to show how to
0064  * obtain physics objects from the analysis tree, and then write them out
0065  * to a ROOT tree and file for further offline analysis.
0066  */
0067
0068 /**
0069  * Constructor of module
0070  */
0071 InttAna::InttAna(const std::string &name, const std::string &filename)
0072     : SubsysReco(name), _rawModule(nullptr), fname_(filename), anafile_(nullptr), h_zvtxseed_(nullptr)
0073 {
0074   xbeam_ = ybeam_ = 0.0;
0075 }
0076
0077 /**
0078  * Destructor of module
0079  */
0080 InttAna::~InttAna()
0081 {
0082 }
0083
0084 /**
0085  * Initialize the module and prepare looping over events
0086  */
0087 int InttAna::Init(PHCompositeNode * /*topNode*/)
0088 {
0089   if (Verbosity() > 5)
0090   {
0091     std::cout << "Beginning Init in InttAna" << std::endl;
0092   }
0093
0094   anafile_ = new TFile(fname_.c_str(), "recreate");
0095   h_dca2d_zero = new TH1F("h_dca2d_zero", "DCA2D to 0", 500, -10, 10);
0096   h2_dca2d_zero = new TH2F("h2_dca2d_zero", "DCA2D to 0 vs nclus", 500, -10, 10, 100, 0, 10000);
0097   h2_dca2d_len = new TH2F("h2_dca2d_len", "DCA2D to 0 vs len", 500, -10, 10, 500, -10, 10);
0098   h_zvtx = new TH1F("h_zvtx", "Zvertex", 400, -40, 40);
0099   h_eta = new TH1F("h_eta", "h_eta", 100, -6, 6);
0100   h_phi = new TH1F("h_phi", "h_phi", 100, -6, 6);
0101   h_theta = new TH1F("h_theta", "h_theta", 100, -4, 4);
0102
0103   h_ntp_clus = new TNtuple("ntp_clus", "Cluster Ntuple", "nclus:nclus2:bco_full:evt:size:adc:x:y:z:lay:lad:sen:x_vtx:y_vtx:z_vtx");
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");
0105
0106   h_zvtxseed_ = new TH1F("h_zvtxseed", "Zvertex Seed histogram", 200, -50, 50);
0107
0108   // TODO:
0109   m_hepmctree = new TTree("hepmctree", "A tree with hepmc truth particles");
0110   m_hepmctree->Branch("m_evt", &m_evt, "m_evt/I");
0111   m_hepmctree->Branch("m_xvtx", &m_xvtx, "m_xvtx/D");
0112   m_hepmctree->Branch("m_yvtx", &m_yvtx, "m_yvtx/D");
0113   m_hepmctree->Branch("m_zvtx", &m_zvtx, "m_zvtx/D");
0114   m_hepmctree->Branch("m_partid1", &m_partid1, "m_partid1/I");
0115   m_hepmctree->Branch("m_partid2", &m_partid2, "m_partid2/I");
0116   m_hepmctree->Branch("m_x1", &m_x1, "m_x1/D");
0117   m_hepmctree->Branch("m_x2", &m_x2, "m_x2/D");
0118   m_hepmctree->Branch("m_mpi", &m_mpi, "m_mpi/I");
0119   m_hepmctree->Branch("m_process_id", &m_process_id, "m_process_id/I");
0120   m_hepmctree->Branch("m_truthenergy", &m_truthenergy, "m_truthenergy/D");
0121   m_hepmctree->Branch("m_trutheta", &m_trutheta, "m_trutheta/D");
0122   m_hepmctree->Branch("m_truththeta", &m_truththeta, "m_truththeta/D");
0123   m_hepmctree->Branch("m_truthphi", &m_truthphi, "m_truthphi/D");
0124   m_hepmctree->Branch("m_status", &m_status, "m_status/I");
0125   m_hepmctree->Branch("m_truthpx", &m_truthpx, "m_truthpx/D");
0126   m_hepmctree->Branch("m_truthpy", &m_truthpy, "m_truthpy/D");
0127   m_hepmctree->Branch("m_truthpz", &m_truthpz, "m_truthpz/D");
0128   m_hepmctree->Branch("m_truthpt", &m_truthpt, "m_truthpt/D");
0129   m_hepmctree->Branch("m_numparticlesinevent", &m_numparticlesinevent, "m_numparticlesinevent/I");
0130   m_hepmctree->Branch("m_truthpid", &m_truthpid, "m_truthpid/I");
0131
0132   return 0;
0133 }
0134
0135 int InttAna::InitRun(PHCompositeNode * /*topNode*/)
0136 {
0137   cout << "InttAna::InitRun beamcenter " << xbeam_ << " " << ybeam_ << endl;
0138
0139   return 0;
0140 }
0141
0142 /**
0143  * Main workhorse function where each event is looped over and
0144  * data from each event is collected from the node tree for analysis
0145  */
0146 int InttAna::process_event(PHCompositeNode *topNode)
0147 {
0148   static int ievt = 0;
0149   cout << "InttEvt::process evt : " << ievt++ << endl;
0150
0151   if (Verbosity() > 5)
0152   {
0153     std::cout << "Beginning process_event in AnaTutorial" << std::endl;
0154   }
0155   // TODO:
0156   // getHEPMCTruth(topNode);
0157
0158   // readRawHit(topNode);
0159
0160   ActsGeometry *m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0161   if (!m_tGeometry)
0162   {
0163     std::cout << PHWHERE << "No ActsGeometry on node tree. Bailing." << std::endl;
0164     return Fun4AllReturnCodes::ABORTEVENT;
0165   }
0166   // TODO:
0167   PHHepMCGenEventMap *hepmceventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0168   if (!hepmceventmap)
0169   {
0170     cout << PHWHERE << "hepmceventmap node is missing." << endl;
0171     // return Fun4AllReturnCodes::ABORTEVENT;
0172   }
0173
0174   // TODO:
0175   PHG4InEvent *phg4inevent = findNode::getClass<PHG4InEvent>(topNode, "PHG4INEVENT");
0176   if (!phg4inevent)
0177   {
0178     cout << PHWHERE << "PHG4INEVENT node is missing." << endl;
0179     // return Fun4AllReturnCodes::ABORTEVENT;
0180   }
0181
0182   TrkrClusterContainer *m_clusterMap =
0183       findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0184   if (!m_clusterMap)
0185   {
0186     cout << PHWHERE << "TrkrClusterContainer node is missing." << endl;
0187     return Fun4AllReturnCodes::ABORTEVENT;
0188   }
0189
0190   GlobalVertexMap *vertexs =
0191       findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0192   if (!vertexs)
0193   {
0194     cout << PHWHERE << "GlobalVertexMap node is missing." << endl;
0195     // return Fun4AllReturnCodes::ABORTEVENT;
0196   }
0197
0198   /////////////
0199   //  InttEvent from RawModule (not standard way)
0200   InttEvent *inttEvt = nullptr;
0201   if (_rawModule)
0202   {
0203     inttEvt = _rawModule->getInttEvent();
0204   }
0205   uint64_t bco = 0;
0206   //-- int evtseq = -1;
0207   if (inttEvt != NULL)
0208   {
0209     bco = inttEvt->bco;
0210     //-- evtseq = inttEvt->evtSeq;
0211   }
0212
0213   // xvtx_sim = 0; //vertexs->find(GlobalVertex::TRUTH)->second->get_x();
0214   // yvtx_sim = 0; //vertexs->find(GlobalVertex::TRUTH)->second->get_y();
0215   // zvtx_sim = 0; //vertexs->find(GlobalVertex::TRUTH)->second->get_z();
0216   xvtx_sim = vertexs->find(GlobalVertex::TRUTH)->second->get_x();
0217   yvtx_sim = vertexs->find(GlobalVertex::TRUTH)->second->get_y();
0218   zvtx_sim = vertexs->find(GlobalVertex::TRUTH)->second->get_z();
0219
0220   // TODO:
0221   if (hepmceventmap != nullptr)
0222     getHEPMCTruth(topNode);
0223   else if (phg4inevent != nullptr)
0224     getPHG4Particle(topNode);
0225   else
0226   {
0227     cout << "no inputinfo available. hepmctree is empty" << endl;
0228   }
0229
0230   /////////////
0231   int nclusadd = 0, nclusadd2 = 0;
0232   for (unsigned int inttlayer = 0; inttlayer < 4; inttlayer++)
0233   {
0234     for (const auto &hitsetkey : m_clusterMap->getHitSetKeys(TrkrDefs::TrkrId::inttId, inttlayer + 3))
0235     {
0236       auto range = m_clusterMap->getClusters(hitsetkey);
0237
0238       for (auto clusIter = range.first; clusIter != range.second; ++clusIter)
0239       {
0240         nclusadd++;
0241         const auto cluster = clusIter->second;
0242         if (cluster->getAdc() > 40)
0243           nclusadd2++;
0244       }
0245     }
0246   }
0247   /////////////
0248
0249   static int evtCount = 0;
0250
0251   // cout << "evtCount : " << evtCount << " " << evtseq << " " << hex << bco << dec << endl;
0252   cout << "nclus : " << nclusadd << " " << nclusadd2 << endl;
0253
0254   struct ClustInfo
0255   {
0256     int layer;
0257     Acts::Vector3 pos;
0258   };
0259
0260   float ntpval[20];
0261   int nCluster = 0;
0262   bool exceedNwrite = false;
0263   // std::vector<Acts::Vector3> clusters[2]; // inner=0; outer=1
0264   std::vector<ClustInfo> clusters[2]; // inner=0; outer=1
0265   for (unsigned int inttlayer = 0; inttlayer < 4; inttlayer++)
0266   {
0267     int layer = (inttlayer < 2 ? 0 : 1);
0268     for (const auto &hitsetkey : m_clusterMap->getHitSetKeys(TrkrDefs::TrkrId::inttId, inttlayer + 3))
0269     {
0270       auto range = m_clusterMap->getClusters(hitsetkey);
0271
0272       for (auto clusIter = range.first; clusIter != range.second; ++clusIter)
0273       {
0274         const auto cluskey = clusIter->first;
0275         const auto cluster = clusIter->second;
0276         const auto globalPos = m_tGeometry->getGlobalPosition(cluskey, cluster);
0277         int ladder_z = InttDefs::getLadderZId(cluskey);
0278         int ladder_phi = InttDefs::getLadderPhiId(cluskey);
0279         int size = cluster->getSize();
0280
0281         if (nCluster < 5)
0282         {
0283           cout << "xyz : " << globalPos.x() << " " << globalPos.y() << " " << globalPos.z() << " :  "
0284                << cluster->getAdc() << " " << size << " " << inttlayer << " " << ladder_z << " " << ladder_phi << endl;
0285         }
0286         else
0287         {
0288           if (!exceedNwrite)
0289           {
0290             cout << " exceed : ncluster limit.  no more cluster xyz printed" << endl;
0291             exceedNwrite = true;
0292           }
0293         }
0294
0295         ClustInfo info;
0296         info.layer = inttlayer;
0297         info.pos = globalPos;
0298
0299         // clusters[layer].push_back(globalPos);
0300         clusters[layer].push_back(info);
0301         nCluster++;
0302
0303         ntpval[0] = nclusadd;  // bco_full
0304         ntpval[1] = nclusadd2; // bco_full
0305         ntpval[2] = bco;       // bco_full
0306         // ntpval[3] = evtseq;            // evtCount;
0307         ntpval[3] = evtCount;          // evtCount;
0308         ntpval[4] = size;              // size
0309         ntpval[5] = cluster->getAdc(); // size
0310         ntpval[6] = globalPos.x();
0311         ntpval[7] = globalPos.y();
0312         ntpval[8] = globalPos.z();
0313         ntpval[9] = inttlayer;
0314         ntpval[10] = ladder_phi;
0315         ntpval[11] = ladder_z;
0316         ntpval[12] = xvtx_sim;
0317         ntpval[13] = yvtx_sim;
0318         ntpval[14] = zvtx_sim;
0319         h_ntp_clus->Fill(ntpval);
0320
0321         // cout << "track event = " << evtCount << endl;
0322
0323         // = new TNtuple("ntp_clus", "Cluster Ntuple", "bco_full:evt:size:adc:x:y:z:lay:lad:sen");
0324       }
0325     }
0326   }
0327   cout << "nCluster : " << nCluster << "-----" << endl;
0328
0329   /////////////////////////////////
0330   /////////////////////////////////
0331
0332   /////////////////////////////////
0333   /////////////////////////////////
0334
0335   h_zvtxseed_->Reset();
0336
0337   if (nCluster < 300)
0338   // if(nCluster<500)
0339   {
0340     float ntpval2[20];
0341     Acts::Vector3 beamspot = Acts::Vector3(xbeam_, ybeam_, 0);
0342     vector<double> vz_array;
0343     for (auto c1 = clusters[0].begin(); c1 != clusters[0].end(); ++c1) // inner
0344     {
0345       for (auto c2 = clusters[1].begin(); c2 != clusters[1].end(); ++c2) // outer
0346       {
0347         Acts::Vector3 p1 = c1->pos - beamspot;
0348         Acts::Vector3 p2 = c2->pos - beamspot;
0349         Acts::Vector3 u = p2 - p1;
0350         double unorm = sqrt(u.x() * u.x() + u.y() * u.y());
0351         if (unorm < 0.00001)
0352           continue;
0353
0354         // skip bad compbination
0355         double p1_ang = atan2(p1.y(), p1.x());
0356         double p2_ang = atan2(p2.y(), p2.x());
0357         double d_ang = p2_ang - p1_ang;
0358
0359         // unit vector in 2D
0360         double ux = u.x() / unorm;
0361         double uy = u.y() / unorm;
0362         double uz = u.z() / unorm; // same normalization factor(2D) is multiplied
0363
0364         Acts::Vector3 p0 = beamspot - p1;
0365
0366         double dca_p0 = p0.x() * uy - p0.y() * ux; // cross product of p0 x u
0367         double len_p0 = p0.x() * ux + p0.y() * uy; // dot product of p0 . u
0368
0369         // beam center in X-Y plane
0370         double vx = len_p0 * ux + p1.x();
0371         double vy = len_p0 * uy + p1.y();
0372
0373         double vz = len_p0 * uz + p1.z();
0374
0375         // if(unorm>4.5||d_ang<-0.005*3.1415||0.005*3.1415<d_ang)
0376         // if(unorm>4.5||d_ang<-0.25*3.1415||0.25*3.1415<d_ang)
0377         if (fabs(d_ang) < 0.1 && fabs(dca_p0) < 0.5)
0378         {
0379           h_zvtxseed_->Fill(vz);
0380           vz_array.push_back(vz);
0381         }
0382
0383         h_dca2d_zero->Fill(dca_p0);
0384         h2_dca2d_zero->Fill(dca_p0, nCluster);
0385         h2_dca2d_len->Fill(dca_p0, len_p0);
0386         // cout<<"dca : "<<dca_p0<<endl;
0387         //
0388         ntpval2[0] = nclusadd;
0389         ntpval2[1] = nclusadd2;
0390         ntpval2[2] = 0;
0391         ntpval2[3] = evtCount;
0392         ntpval2[4] = p1_ang;
0393         ntpval2[5] = p2_ang;
0394         ntpval2[6] = p1.z();
0395         ntpval2[7] = p2.z();
0396         ntpval2[8] = dca_p0;
0397         ntpval2[9] = len_p0;
0398         ntpval2[10] = unorm;
0399         ntpval2[11] = c1->layer;
0400         ntpval2[12] = c2->layer;
0401         ntpval2[13] = vx;
0402         ntpval2[14] = vy;
0403         ntpval2[15] = vz;
0404         h_ntp_cluspair->Fill(ntpval2);
0405         // h_ntp_cluspair = new TNtuple("ntp_cluspair", "Cluster Pair Ntuple", "nclus:bco_full:evt:ang1:ang2:dca2d:len:unorm");
0406       }
0407     }
0408
0409     // calculate trancated mean of DCA~Z histogram as Z-vertex position
0410     double zvtx = -9999.;
0411
0412     if (vz_array.size() > 3)
0413     {
0414       double zbin = h_zvtxseed_->GetMaximumBin();
0415       double zcenter = h_zvtxseed_->GetBinCenter(zbin);
0416       //-- double zmean = h_zvtxseed_->GetMean();
0417       double zrms = h_zvtxseed_->GetRMS();
0418       if (zrms < 20)
0419         zrms = 20;
0420
0421       double zmax = zcenter + zrms; // 1 sigma
0422       double zmin = zcenter - zrms; // 1 sigma
0423
0424       double zsum = 0.;
0425       int zcount = 0;
0426       for (auto iz = vz_array.begin(); iz != vz_array.end(); ++iz)
0427       {
0428         double vz = (*iz);
0429         if (zmin < vz && vz < zmax)
0430         {
0431           zsum += vz;
0432           zcount++;
0433         }
0434       }
0435       if (zcount > 0)
0436         zvtx = zsum / zcount;
0437
0438       // cout << "ZVTX: " << zvtx << " " << zcenter << " " << zmean << " " << zrms << " " << zbin << endl;
0439     }
0440
0441     h_zvtx->Fill(zvtx);
0442   }
0443
0444   evtCount++;
0445
0446   return Fun4AllReturnCodes::EVENT_OK;
0447 }
0448
0449 /**
0450  * End the module and finish any data collection. Clean up any remaining
0451  * loose ends
0452  */
0453 int InttAna::End(PHCompositeNode * /*topNode*/)
0454 {
0455   if (Verbosity() > 1)
0456   {
0457     std::cout << "Ending InttAna analysis package" << std::endl;
0458   }
0459   anafile_->cd();
0460   // TODO:
0461   m_hepmctree->Write();
0462   if (anafile_ != nullptr)
0463   {
0464     anafile_->Write();
0465     anafile_->Close();
0466   }
0467
0468   return 0;
0469 }
0470 // TODO:
0471
0472 void InttAna::readRawHit(PHCompositeNode *topNode)
0473 {
0474   TrkrHitSetContainer *m_hits = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0475   if (!m_hits)
0476   {
0477     cout << PHWHERE << "ERROR: Can't find node TRKR_HITSET" << endl;
0478     return;
0479   }
0480
0481   // uint8_t getLadderZId(TrkrDefs::hitsetkey key);
0482   // uint8_t getLadderPhiId(TrkrDefs::hitsetkey key);
0483   // int getTimeBucketId(TrkrDefs::hitsetkey key);
0484   //
0485   // uint16_t getCol(TrkrDefs::hitkey key); // z-strip = offline chip
0486   // uint16_t getRow(TrkrDefs::hitkey key); // y-strip = offline channel
0487
0488   // loop over the InttHitSet objects
0489   TrkrHitSetContainer::ConstRange hitsetrange =
0490       m_hits->getHitSets(TrkrDefs::TrkrId::inttId);
0491   for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first;
0492        hitsetitr != hitsetrange.second;
0493        ++hitsetitr)
0494   {
0495     // Each hitset contains only hits that are clusterizable - i.e. belong to a single sensor
0496     TrkrHitSet *hitset = hitsetitr->second;
0497
0498     if (Verbosity() > 1)
0499       cout << "InttClusterizer found hitsetkey " << hitsetitr->first << endl;
0500     if (Verbosity() > 2)
0501       hitset->identify();
0502
0503     // we have a single hitset, get the info that identifies the sensor
0504     // int layer = TrkrDefs::getLayer(hitsetitr->first);
0505     // int ladder_z_index = InttDefs::getLadderZId(hitsetitr->first);
0506     // int ladder_phi_index = InttDefs::getLadderPhiId(hitsetitr->first);
0507     // int ladder_bco_index = InttDefs::getTimeBucketId(hitsetitr->first);
0508
0509     // cout << "layer " << layer << " " << ladder_z_index << " " << ladder_phi_index << " " << ladder_bco_index << endl;
0510
0511     //// we will need the geometry object for this layer to get the global position
0512     // CylinderGeomIntt* geom = dynamic_cast<CylinderGeomIntt*>(geom_container->GetLayerGeom(layer));
0513     // float pitch = geom->get_strip_y_spacing();
0514     // float length = geom->get_strip_z_spacing();
0515
0516     // fill a vector of hits to make things easier - gets every hit in the hitset
0517     std::vector<std::pair<TrkrDefs::hitkey, TrkrHit *>> hitvec;
0518     TrkrHitSet::ConstRange hitrangei = hitset->getHits();
0519     for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
0520          hitr != hitrangei.second;
0521          ++hitr)
0522     {
0523       hitvec.push_back(make_pair(hitr->first, hitr->second));
0524       // int chip = InttDefs::getCol(hitr->first);
0525       // int chan = InttDefs::getRow(hitr->first);
0526       // int adc = (hitr->second)->getAdc();
0527       // cout << "     hit : " << chip << " " << chan << " " << adc << endl;
0528     }
0529     // cout << "hitvec.size(): " << hitvec.size() << endl;
0530   }
0531 }
0532
0533 void InttAna::getHEPMCTruth(PHCompositeNode *topNode)
0534 {
0535   cout << "getHEPMCTruth!!!" << endl;
0536   /// Get the node from the node tree
0537   PHHepMCGenEventMap *hepmceventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0538
0539   /// If the node was not properly put on the tree, return
0540   if (!hepmceventmap)
0541   {
0542     std::cout << PHWHERE
0543               << "HEPMC event map node is missing, can't collected HEPMC truth particles"
0544               << std::endl;
0545     return;
0546   }
0547
0548   /// Could have some print statements for debugging with verbosity
0549   if (Verbosity() > 1)
0550   {
0551     std::cout << "Getting HEPMC truth particles " << std::endl;
0552   }
0553
0554   /// You can iterate over the number of events in a hepmc event
0555   /// for pile up events where you have multiple hard scatterings per bunch crossing
0556   for (PHHepMCGenEventMap::ConstIter eventIter = hepmceventmap->begin();
0557        eventIter != hepmceventmap->end();
0558        ++eventIter)
0559   {
0560     /// Get the event
0561     PHHepMCGenEvent *hepmcevent = eventIter->second;
0562
0563     if (hepmcevent)
0564     {
0565       /// Get the event characteristics, inherited from HepMC classes
0566       HepMC::GenEvent *truthevent = hepmcevent->getEvent();
0567       if (!truthevent)
0568       {
0569         std::cout << PHWHERE
0570                   << "no evt pointer under phhepmvgeneventmap found "
0571                   << std::endl;
0572         return;
0573       }
0574
0575       /// Get the parton info
0576       HepMC::PdfInfo *pdfinfo = truthevent->pdf_info();
0577
0578       /// Get the parton info as determined from HEPMC
0579       m_partid1 = pdfinfo->id1();
0580       m_partid2 = pdfinfo->id2();
0581       m_x1 = pdfinfo->x1();
0582       m_x2 = pdfinfo->x2();
0583
0584       /// Are there multiple partonic intercations in a p+p event
0585       m_mpi = truthevent->mpi();
0586
0587       /// Get the PYTHIA signal process id identifying the 2-to-2 hard process
0588       m_process_id = truthevent->signal_process_id();
0589
0590       if (Verbosity() > 2)
0591       {
0592         std::cout << " Iterating over an event" << std::endl;
0593       }
0594
0595       /// Loop over all the truth particles and get their information
0596       for (HepMC::GenEvent::particle_const_iterator iter = truthevent->particles_begin();
0597            iter != truthevent->particles_end();
0598            ++iter)
0599       {
0600         /// Get each pythia particle characteristics
0601         m_truthenergy = (*iter)->momentum().e();
0602         m_truthpid = (*iter)->pdg_id();
0603
0604         // TODO: truth information
0605         m_xvtx = xvtx_sim;
0606         m_yvtx = yvtx_sim;
0607         m_zvtx = zvtx_sim;
0608         m_trutheta = (*iter)->momentum().pseudoRapidity();
0609         m_truththeta = 2 * atan(exp(-m_trutheta));
0610         // h_eta->Fill(m_trutheta);
0611         m_truthphi = (*iter)->momentum().phi();
0612         m_truthpx = (*iter)->momentum().px();
0613         m_truthpy = (*iter)->momentum().py();
0614         m_truthpz = (*iter)->momentum().pz();
0615
0616         m_truthpt = sqrt(m_truthpx * m_truthpx + m_truthpy * m_truthpy);
0617         m_status = (*iter)->status();
0618
0619         // m_truthphi = atan2(m_truthpy,m_truthpx);
0620
0621         // cout << "status: " << (*iter)->status() << " " << m_truthpid << endl;
0622
0623         h_phi->Fill(m_truthphi);
0624         h_theta->Fill(m_truththeta);
0625
0626         /*std::cout << "m_trutheta = " << m_trutheta << "  m_truthphi = " << m_truthphi << std::endl;*/
0627
0628         /// Fill the truth tree
0629         m_hepmctree->Fill();
0630         m_numparticlesinevent++;
0631       }
0632
0633       // for (HepMC::GenEventVertexRange::vertex_const_iterator viter = truthevent->vertices_begin();
0634       //      viter != truthevent->vertices_end();
0635       //      ++viter)
0636       // {
0637       //   /// Get each pythia particle characteristics
0638       //   m_vertex = (*viter)->vertex_range();
0639       //   cout<<m_vertex<<endl;
0640
0641       //   /// Fill the truth tree
0642       //   m_hepmctree->Fill();
0643       // }
0644     }
0645     // cout << "truth event = " << m_evt << endl;
0646     m_evt++;
0647   }
0648 }
0649
0650 void InttAna::getPHG4Particle(PHCompositeNode *topNode)
0651 {
0652   cout << "getPHG4Particle!!!" << endl;
0653   /// Get the node from the node tree
0654   PHG4InEvent *phg4inevent = findNode::getClass<PHG4InEvent>(topNode, "PHG4INEVENT");
0655
0656   /// If the node was not properly put on the tree, return
0657   if (!phg4inevent)
0658   {
0659     std::cout << PHWHERE
0660               << "PHG4InEvent node is missing, can't collected PHG4 truth particles"
0661               << std::endl;
0662     return;
0663   }
0664
0665   /// Could have some print statements for debugging with verbosity
0666   if (Verbosity() > 1)
0667   {
0668     std::cout << "Getting PHG4InEvent truth particles " << std::endl;
0669   }
0670
0671   /// You can iterate over the number of events in a hepmc event
0672   /// for pile up events where you have multiple hard scatterings per bunch crossing
0673   std::map<int, PHG4VtxPoint *>::const_iterator vtxiter;
0674   std::multimap<int, PHG4Particle *>::const_iterator particle_iter;
0675   std::pair<std::map<int, PHG4VtxPoint *>::const_iterator, std::map<int, PHG4VtxPoint *>::const_iterator> vtxbegin_end = phg4inevent->GetVertices();
0676
0677   for (vtxiter = vtxbegin_end.first; vtxiter != vtxbegin_end.second; ++vtxiter)
0678   {
0679     //*fout << "vtx number: " << vtxiter->first << std::endl;
0680     //(*vtxiter->second).identify(*fout);
0681     std::pair<std::multimap<int, PHG4Particle *>::const_iterator,
0682               std::multimap<int, PHG4Particle *>::const_iterator>
0683         particlebegin_end = phg4inevent->GetParticles(vtxiter->first);
0684
0685     PHG4VtxPoint *vtx = vtxiter->second;
0686     m_xvtx = vtx->get_x();
0687     m_yvtx = vtx->get_y();
0688     m_zvtx = vtx->get_z();
0689     // virtual double get_t() const { return std::numeric_limits<double>::quiet_NaN(); }
0690
0691     for (particle_iter = particlebegin_end.first; particle_iter != particlebegin_end.second; ++particle_iter)
0692     {
0693       //      (particle_iter->second)->identify(*fout);
0694       PHG4Particle *part = particle_iter->second;
0695
0696       m_truthenergy = part->get_e();
0697       m_truthpid = part->get_pid();
0698
0699       m_truthpx = part->get_px();
0700       m_truthpy = part->get_py();
0701       m_truthpz = part->get_pz();
0702       m_truthphi = atan2(m_truthpy, m_truthpx);
0703       m_truthpt = sqrt(m_truthpx * m_truthpx + m_truthpy * m_truthpy);
0704
0705       m_truththeta = atan2(m_truthpt, m_truthpz);
0706       m_trutheta = -log(tan(0.5 * m_truththeta));
0707       m_status = 1;
0708
0709       // part->get_pid() const override { return fpid; }
0710       // part->get_name() const override { return fname; }
0711       // part->get_e() const override { return fe; }
0712       // part->get_track_id();
0713       // part->get_vtx_id() const override { return vtxid; }
0714       // part->get_parent_id() const override { return parentid; }
0715       // int get_primary_id() const override { return primaryid; }
0716
0717       /// Fill the truth tree
0718       h_phi->Fill(m_truthphi);
0719       h_theta->Fill(m_truththeta);
0720
0721       m_hepmctree->Fill();
0722       m_numparticlesinevent++;
0723     }
0724   }
0725   cout << "truth event = " << m_evt << endl;
0726   m_evt++;
0727
0728   /*
0729     for (PHHepMCGenEventMap::ConstIter eventIter = hepmceventmap->begin();
0730          eventIter != hepmceventmap->end();
0731          ++eventIter)
0732     {
0733       /// Get the event
0734       PHHepMCGenEvent *hepmcevent = eventIter->second;
0735
0736       if (hepmcevent)
0737       {
0738         /// Get the parton info
0739         HepMC::PdfInfo *pdfinfo = truthevent->pdf_info();
0740
0741         /// Get the parton info as determined from HEPMC
0742         m_partid1 = pdfinfo->id1();
0743         m_partid2 = pdfinfo->id2();
0744         m_x1 = pdfinfo->x1();
0745         m_x2 = pdfinfo->x2();
0746
0747         /// Are there multiple partonic intercations in a p+p event
0748         m_mpi = truthevent->mpi();
0749
0750         /// Get the PYTHIA signal process id identifying the 2-to-2 hard process
0751         m_process_id = truthevent->signal_process_id();
0752       }
0753       // cout << "truth event = " << m_evt << endl;
0754       m_evt++;
0755     }
0756   */
0757 }