hbcup/src/AnaUPC.cc

0001 #include "AnaUPC.h"
0002
0003 /// Cluster/Calorimeter includes
0004 /*
0005 #include <calobase/RawCluster.h>
0006 #include <calobase/RawClusterContainer.h>
0007 #include <calobase/RawClusterUtility.h>
0008 #include <calobase/RawTower.h>
0009 #include <calobase/RawTowerContainer.h>
0010 #include <calobase/RawTowerGeom.h>
0011 #include <calobase/RawTowerGeomContainer.h>
0012 #include <calotrigger/CaloTriggerInfo.h>
0013 */
0014
0015 /// Tracking includes
0016 #include <globalvertex/GlobalVertex.h>
0017 #include <globalvertex/GlobalVertexMap.h>
0018 #include <trackbase_historic/SvtxTrack.h>
0019 #include <trackbase_historic/SvtxTrackMap.h>
0020 #include <globalvertex/SvtxVertex.h>
0021 #include <globalvertex/SvtxVertexMap.h>
0022
0023 /// Truth evaluation includes
0024 //#include <g4eval/JetEvalStack.h>
0025 #include <g4eval/SvtxEvalStack.h>
0026
0027 /// HEPMC truth includes
0028 #pragma GCC diagnostic push
0029 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
0030 #include <HepMC/GenEvent.h>
0031 #include <HepMC/GenVertex.h>
0032 #pragma GCC diagnostic pop
0033
0034 #include <phhepmc/PHHepMCGenEvent.h>
0035 #include <phhepmc/PHHepMCGenEventMap.h>
0036
0037 /// Fun4All includes
0038 #include <fun4all/Fun4AllReturnCodes.h>
0039 #include <g4main/PHG4Hit.h>
0040 #include <g4main/PHG4Particle.h>
0041 #include <g4main/PHG4TruthInfoContainer.h>
0042 #include <phool/PHCompositeNode.h>
0043 #include <phool/getClass.h>
0044 #include <ffaobjects/EventHeader.h>
0045
0046 /// ROOT includes
0047 #include <TFile.h>
0048 #include <TH1.h>
0049 #include <TH2.h>
0050 #include <TNtuple.h>
0051 #include <TTree.h>
0052 #include <Math/Vector4D.h>
0053 #include <Math/VectorUtil.h>
0054 #include <TDatabasePDG.h>
0055
0056
0057 /// C++ includes
0058 #include <cassert>
0059 #include <cmath>
0060 #include <sstream>
0061 #include <string>
0062
0063 /**
0064  * Constructor of module
0065  */
0066 AnaUPC::AnaUPC(const std::string &name, const std::string &filename)
0067   : SubsysReco(name)
0068   , m_outfilename(filename)
0069   //, m_mincluspt(0.25)
0070   , m_analyzeTracks(true)
0071   //, m_analyzeClusters(true)
0072   , m_analyzeTruth(true)
0073 {
0074   /// Initialize variables and trees so we don't accidentally access
0075   /// memory that was never allocated
0076   initializeVariables();
0077 }
0078
0079 /**
0080  * Destructor of module
0081  */
0082 AnaUPC::~AnaUPC()
0083 {
0084   std::cout << PHWHERE << "In ~AnaUPC()" << std::endl;
0085   /*
0086   if ( m_hepmctree!=nullptr )  delete m_hepmctree;
0087   if ( m_tracktree!=nullptr )  delete m_tracktree;
0088   if ( m_truthtree!=nullptr )  delete m_truthtree;
0089   if ( m_pairtree!=nullptr )   delete m_pairtree;
0090   if ( m_globaltree!=nullptr ) delete m_globaltree;
0091   */
0092 }
0093
0094 /**
0095  * Initialize the module and prepare looping over events
0096  */
0097 int AnaUPC::Init(PHCompositeNode * /*topNode*/)
0098 {
0099   if (Verbosity() > 5)
0100   {
0101     std::cout << "Beginning Init in AnaUPC" << std::endl;
0102   }
0103
0104   m_outfile = new TFile(m_outfilename.c_str(), "RECREATE");
0105   initializeTrees();
0106
0107   h_phi[0] = new TH1F("h_phi", "#phi [rad]", 60, -M_PI, M_PI);
0108   h2_eta_phi[0] = new TH2F("h2_phi_eta", ";#eta;#phi [rad]", 24, -5.0, 5.0, 60, -M_PI, M_PI);
0109   h_mass[0] = new TH1F("h_mass", "mass [GeV]", 1200, 0, 6);
0110   h_pt[0] = new TH1F("h_pt", "p_{T}", 200, 0, 2);
0111   h_y[0] = new TH1F("h_y", "y", 24, -1.2, 1.2);
0112   h_eta[0] = new TH1F("h_eta", "#eta", 24, -5.0, 5.0);
0113
0114   // like-sign pairs
0115   h_phi[1] = new TH1F("h_phi_ls", "#phi [rad]", 60, -M_PI, M_PI);
0116   h2_eta_phi[1] = new TH2F("h2_phi_eta_ls", ";#eta;#phi [rad]", 24, -5.0, 5.0, 60, -M_PI, M_PI);
0117   h_mass[1] = new TH1F("h_mass_ls", "mass [GeV]", 1200, 0, 6);
0118   h_pt[1] = new TH1F("h_pt_ls", "p_{T}", 200, 0, 2);
0119   h_y[1] = new TH1F("h_y_ls", "y", 24, -1.2, 1.2);
0120   h_eta[1] = new TH1F("h_eta_ls", "#eta", 24, -5.0, 5.0);
0121
0122   h_trig = new TH1F("h_trig", "trig", 16, 0, 16);
0123   h_ntracks = new TH1F("h_ntracks", "num tracks", 2000, 0, 2000);
0124   h2_ntrksvsb = new TH2F("h2_ntrksvsb", "num tracks vs b", 220, 0, 22, 2001, -0.5, 2000.5);
0125   h2_ntrksvsb->SetXTitle("b [fm]");
0126   h2_ntrksvsb->SetYTitle("N_{TRKS}");
0127
0128   h_b_mb = new TH1F("h_b_mb", "b, MB events", 200, 0, 20);
0129   h_npart_mb = new TH1F("h_npart_mb", "npart, MB events", 401, -0.5, 400.5);
0130   h_ncoll_mb = new TH1F("h_ncoll_mb", "ncoll, MB events", 1301, -0.5, 1300.5);
0131   h_b = new TH1F("h_b", "b", 200, 0, 20);
0132   h_npart = new TH1F("h_npart", "npart", 401, -0.5, 400.5);
0133   h_ncoll = new TH1F("h_ncoll", "ncoll", 1301, -0.5, 1300.5);
0134
0135   return 0;
0136 }
0137
0138 /**
0139  * Main workhorse function where each event is looped over and
0140  * data from each event is collected from the node tree for analysis
0141  */
0142 int AnaUPC::GetNodes(PHCompositeNode *topNode)
0143 {
0144   /// EventHeader node
0145   _evthdr = findNode::getClass<EventHeader>(topNode, "EventHeader");
0146
0147   /// SVTX tracks node
0148   _trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0149   if (!_trackmap)
0150   {
0151     std::cout << PHWHERE << "SvtxTrackMap node is missing, can't collect tracks" << std::endl;
0152     return Fun4AllReturnCodes::ABORTEVENT;
0153   }
0154
0155   /// G4 truth particle node
0156   _truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0157   if (!_truthinfo)
0158   {
0159     static int ctr = 0;
0160     if ( ctr<4 )
0161     {
0162       std::cout << PHWHERE << "PHG4TruthInfoContainer node is missing, can't collect G4 truth particles" << std::endl;
0163       ctr++;
0164     }
0165   }
0166
0167   /// HEPMC info
0168   _genevent_map = findNode::getClass<PHHepMCGenEventMap>(topNode,"PHHepMCGenEventMap");
0169   if (!_genevent_map)
0170   {
0171     static int ctr = 0;
0172     if ( ctr<4 )
0173     {
0174       std::cout << PHWHERE << "PHHepMCGenEventMap node is missing, can't collect HEPMC info" << std::endl;
0175       ctr++;
0176     }
0177   }
0178
0179   return Fun4AllReturnCodes::EVENT_OK;
0180 }
0181
0182 /**
0183  * Main workhorse function where each event is looped over and
0184  * data from each event is collected from the node tree for analysis
0185  */
0186 int AnaUPC::process_event(PHCompositeNode *topNode)
0187 {
0188   if ( m_evt%1000 == 1 )
0189   {
0190     std::cout << "m_evt " << m_evt << std::endl;
0191   }
0192
0193   Reset(topNode);
0194
0195   if (Verbosity() > 5)
0196   {
0197     std::cout << "Beginning process_event in AnaUPC, " << m_evt << std::endl;
0198   }
0199
0200   h_trig->Fill( 0 );  // event counter
0201
0202   /// Get all the data nodes
0203   int status = GetNodes(topNode);
0204   if ( status != Fun4AllReturnCodes::EVENT_OK )
0205   {
0206     return status;
0207   }
0208
0209   /// Get the run and eventnumber
0210   if ( _evthdr )
0211   {
0212     m_run = _evthdr->get_RunNumber();
0213     m_evt = _evthdr->get_EvtSequence();
0214   }
0215
0216   if ( m_ntracks!=0 || m_ntrk_sphenix!= 0 )
0217   {
0218     std::cout << PHWHERE << " ERROR, evt m_ntracks m_ntrk_sphenix = " << m_evt << "\t"
0219       << m_ntracks << "\t" << m_ntrk_sphenix << std::endl;
0220   }
0221
0222   /// Get global info
0223
0224   if ( _genevent_map )
0225   {
0226     if (Verbosity()>5)
0227     {
0228       std::cout << PHWHERE << "processing global info from sim" << std::endl;
0229     }
0230     PHHepMCGenEvent *genevent = (_genevent_map->begin())->second;
0231     HepMC::GenEvent *hepmc_event{nullptr};
0232     HepMC::HeavyIon *hepmc_hi{nullptr};
0233     if (genevent)
0234     {
0235       hepmc_event = genevent->getEvent();
0236       hepmc_hi = hepmc_event->heavy_ion();
0237       if ( hepmc_hi )
0238       {
0239         m_npart_targ =  hepmc_hi->Npart_targ();
0240         m_npart_proj =  hepmc_hi->Npart_proj();
0241         m_npart = m_npart_targ + m_npart_proj;
0242         m_ncoll =  hepmc_hi->Ncoll();
0243         m_ncoll_hard =  hepmc_hi->Ncoll_hard();
0244         //std::cout << "ncoll " << m_ncoll << "\t" << m_ncoll_hard << std::endl;
0245         m_bimpact =  hepmc_hi->impact_parameter();
0246         //std::cout << "b ntracks " << m_bimpact << "\t" << m_ntracks << std::endl;
0247
0248         //m_globaltree->Fill();
0249
0250         h_b_mb->Fill( m_bimpact );
0251         h_npart_mb->Fill( m_npart );
0252         h_ncoll_mb->Fill( m_ncoll );
0253       }
0254     }
0255   }
0256
0257   /// Get the tracks
0258   if (m_analyzeTracks)
0259   {
0260     status = getTracks(topNode);
0261     if ( status != Fun4AllReturnCodes::EVENT_OK )
0262     {
0263       return status;
0264     }
0265   }
0266
0267   /// Get the truth track information
0268   if (m_analyzeTruth && _truthinfo )
0269   {
0270     getPHG4Truth();
0271   }
0272
0273   // Fill Global info for events that pass cuts
0274   if ( _genevent_map )
0275   {
0276     h2_ntrksvsb->Fill( m_bimpact, m_ntracks );
0277
0278     if ( m_ntracks<=3 )
0279     {
0280       h_b->Fill( m_bimpact );
0281       h_npart->Fill( m_npart );
0282       h_ncoll->Fill( m_ncoll );
0283     }
0284
0285     m_globaltree->Fill();
0286   }
0287
0288   /// Get calorimeter information
0289   /*
0290      if (m_analyzeClusters)
0291      {
0292      getEMCalClusters(topNode);
0293      }
0294      */
0295
0296   return Fun4AllReturnCodes::EVENT_OK;
0297 }
0298
0299 /**
0300  * End the module and finish any data collection. Clean up any remaining
0301  * loose ends
0302  */
0303 int AnaUPC::End(PHCompositeNode * /*topNode*/)
0304 {
0305   if (Verbosity() > 5)
0306   {
0307     std::cout << "Ending AnaUPC analysis package" << std::endl;
0308   }
0309
0310   /// Write and close the outfile
0311   m_outfile->Write();
0312   m_outfile->Close();
0313
0314   if (Verbosity() > 1)
0315   {
0316     std::cout << "Finished AnaUPC analysis package" << std::endl;
0317   }
0318
0319   return 0;
0320 }
0321
0322 /**
0323  * This method gets all of the HEPMC truth particles from the node tree
0324  * and stores them in a ROOT TTree. The HEPMC truth particles are what,
0325  * for example, directly comes out of PYTHIA and thus gives you all of
0326  * the associated parton information
0327  */
0328 void AnaUPC::getHEPMCTruth()
0329 {
0330   /// Could have some print statements for debugging with verbosity
0331   if (Verbosity() > 1)
0332   {
0333     std::cout << "Getting HEPMC truth particles " << std::endl;
0334   }
0335
0336   /// You can iterate over the number of events in a hepmc event
0337   /// for pile up events where you have multiple hard scatterings per bunch crossing
0338   for (PHHepMCGenEventMap::ConstIter eventIter = _genevent_map->begin(); eventIter != _genevent_map->end(); ++eventIter)
0339   {
0340     /// Get the event
0341     PHHepMCGenEvent *hepmcevent = eventIter->second;
0342
0343     if (hepmcevent)
0344     {
0345       /// Get the event characteristics, inherited from HepMC classes
0346       HepMC::GenEvent *truthevent = hepmcevent->getEvent();
0347       if (!truthevent)
0348       {
0349         std::cout << PHWHERE << "no evt pointer under phhepmvgeneventmap found " << std::endl;
0350         return;
0351       }
0352
0353       /// Get the parton info
0354       HepMC::PdfInfo *pdfinfo = truthevent->pdf_info();
0355
0356       /// Get the parton info as determined from HEPMC
0357       m_partid1 = pdfinfo->id1();
0358       m_partid2 = pdfinfo->id2();
0359       m_x1 = pdfinfo->x1();
0360       m_x2 = pdfinfo->x2();
0361
0362       /// Are there multiple partonic intercations in a p+p event
0363       m_mpi = truthevent->mpi();
0364
0365       /// Get the PYTHIA signal process id identifying the 2-to-2 hard process
0366       m_process_id = truthevent->signal_process_id();
0367
0368       if (Verbosity() > 2)
0369       {
0370         std::cout << " Iterating over an event" << std::endl;
0371       }
0372       /// Loop over all the truth particles and get their information
0373       for (HepMC::GenEvent::particle_const_iterator iter = truthevent->particles_begin(); iter != truthevent->particles_end(); ++iter)
0374       {
0375         /// Get each pythia particle characteristics
0376         m_truthenergy = (*iter)->momentum().e();
0377         m_truthpid = (*iter)->pdg_id();
0378
0379         m_trutheta = (*iter)->momentum().pseudoRapidity();
0380         m_truthphi = (*iter)->momentum().phi();
0381         m_truthpx = (*iter)->momentum().px();
0382         m_truthpy = (*iter)->momentum().py();
0383         m_truthpz = (*iter)->momentum().pz();
0384         m_truthpt = sqrt(m_truthpx * m_truthpx + m_truthpy * m_truthpy);
0385
0386         /// Fill the truth tree
0387         m_hepmctree->Fill();
0388         m_numparticlesinevent++;
0389       }
0390     }
0391   }
0392 }
0393
0394 /**
0395  * This function collects the truth PHG4 stable particles that
0396  * are produced from PYTHIA, or some other event generator. These
0397  * are the stable particles, e.g. there are not any (for example)
0398  * partons here.
0399  */
0400 void AnaUPC::getPHG4Truth()
0401 {
0402   /// Get the primary particle range
0403   PHG4TruthInfoContainer::Range range = _truthinfo->GetPrimaryParticleRange();
0404
0405   ROOT::Math::XYZTVector v1;
0406
0407   /// Loop over the G4 truth (stable) particles
0408   for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter)
0409   {
0410     /// Get this truth particle
0411     const PHG4Particle *truth = iter->second;
0412
0413     /// Get this particles momentum, etc.
0414     m_truthpx = truth->get_px();
0415     m_truthpy = truth->get_py();
0416     m_truthpz = truth->get_pz();
0417     m_truthp = sqrt(m_truthpx * m_truthpx + m_truthpy * m_truthpy + m_truthpz * m_truthpz);
0418     m_truthenergy = truth->get_e();
0419
0420     v1.SetPxPyPzE( m_truthpx, m_truthpy, m_truthpz, m_truthenergy );
0421     m_truthpt = v1.Pt();
0422     m_truthphi = v1.Phi();
0423     m_trutheta = v1.Eta();
0424
0425     /// Check for nans
0426     if (!std::isfinite(m_trutheta))
0427     {
0428       m_trutheta = -99;
0429     }
0430     m_truthpid = truth->get_pid();
0431
0432     // Get the singleton instance of the PDG database
0433     TDatabasePDG* pdgDB = TDatabasePDG::Instance();
0434
0435     // Get particle information by PDG code
0436     TParticlePDG* particle = pdgDB->GetParticle(m_truthpid);
0437
0438     if (particle)
0439     {
0440       m_truthcharge = particle->Charge();
0441       if ( m_truthcharge != 0 )
0442       {
0443         m_ntrk_mc++;  // found charged track in sphenix accept
0444       }
0445       if ( (fabs(m_trutheta) < 1.1) && (m_truthpt>0.4) && (m_truthcharge!=0) )
0446       {
0447         m_ntrk_sphenix++;
0448         m_truthtree->Fill();
0449       }
0450     }
0451     else
0452     {
0453       std::cout << "Particle not found!" << std::endl;
0454     }
0455
0456     //std::cout << "pid ch\t" << m_truthpid << "\t" << m_truthcharge << std::endl;
0457
0458   }
0459 }
0460
0461 /**
0462  * This method gets the tracks as reconstructed from the tracker.
0463  */
0464 int AnaUPC::getTracks(PHCompositeNode *topNode)
0465 {
0466   // make a cut on low ntracks
0467   m_ntracks = _trackmap->size();
0468   h_ntracks->Fill( m_ntracks );
0469   if (Verbosity() > 1)
0470   {
0471     std::cout << "ntracks " << m_ntracks << std::endl;
0472   }
0473
0474   if ( m_ntracks > 3 || m_ntracks < 2 )
0475   //if ( m_ntracks != 2 )
0476   {
0477     return Fun4AllReturnCodes::DISCARDEVENT;
0478   }
0479
0480   /// EvalStack for truth track matching
0481   if ( !m_svtxEvalStack && _truthinfo )
0482   {
0483     std::cout << "getting svtx eval stack" << std::endl;
0484     m_svtxEvalStack = new SvtxEvalStack(topNode);
0485     m_svtxEvalStack->set_verbosity(Verbosity());
0486   }
0487
0488   SvtxTrackEval *trackeval = nullptr;
0489   if ( m_svtxEvalStack )
0490   {
0491     m_svtxEvalStack->next_event(topNode);
0492
0493     // Get the track evaluator (only available if eval is run during sim production)
0494     trackeval = m_svtxEvalStack->get_track_eval();
0495   }
0496
0497   if (Verbosity() > 1)
0498   {
0499     std::cout << "Get the SVTX tracks " << m_ntracks << std::endl;
0500   }
0501
0502   for (auto &iter : *_trackmap)
0503   {
0504     SvtxTrack *track = iter.second;
0505
0506     /// Get the reconstructed track info
0507     m_tr_px = track->get_px();
0508     m_tr_py = track->get_py();
0509     m_tr_pz = track->get_pz();
0510     m_tr_p = sqrt(m_tr_px * m_tr_px + m_tr_py * m_tr_py + m_tr_pz * m_tr_pz);
0511
0512     m_tr_pt = sqrt(m_tr_px * m_tr_px + m_tr_py * m_tr_py);
0513
0514     // Make some cuts on the track to clean up sample
0515     if (m_tr_pt < 0.5)
0516     {
0517       continue;
0518     }
0519     m_tr_phi = track->get_phi();
0520     m_tr_eta = track->get_eta();
0521
0522     m_charge = track->get_charge();
0523     m_chisq = track->get_chisq();
0524     m_ndf = track->get_ndf();
0525     m_dca = track->get_dca();
0526     m_tr_x = track->get_x();
0527     m_tr_y = track->get_y();
0528     m_tr_z = track->get_z();
0529
0530     /// Get truth track info that matches this reconstructed track
0531     PHG4Particle *truthtrack = nullptr;
0532     if ( trackeval != nullptr )
0533     {
0534       truthtrack = trackeval->max_truth_particle_by_nclusters(track);
0535     }
0536     if ( truthtrack != nullptr )
0537     {
0538       m_truth_is_primary = _truthinfo->is_primary(truthtrack);
0539
0540       m_truthtrackpx = truthtrack->get_px();
0541       m_truthtrackpy = truthtrack->get_py();
0542       m_truthtrackpz = truthtrack->get_pz();
0543       m_truthtrackp = std::sqrt(m_truthtrackpx * m_truthtrackpx + m_truthtrackpy * m_truthtrackpy + m_truthtrackpz * m_truthtrackpz);
0544
0545       m_truthtracke = truthtrack->get_e();
0546
0547       m_truthtrackpt = sqrt(m_truthtrackpx * m_truthtrackpx + m_truthtrackpy * m_truthtrackpy);
0548       m_truthtrackphi = atan(m_truthtrackpy / m_truthtrackpx);
0549       m_truthtracketa = atanh(m_truthtrackpz / m_truthtrackp);
0550       m_truthtrackpid = truthtrack->get_pid();
0551     }
0552     else
0553     {
0554       // Seems that we often miss the truth track?
0555       //std::cout << "Missing truth track" << std::endl;
0556       m_truth_is_primary = -9999;
0557
0558       m_truthtrackpx = 0.;
0559       m_truthtrackpy = 0.;
0560       m_truthtrackpz = 0.;
0561       m_truthtrackp = 0.;
0562
0563       m_truthtracke = 0.;
0564
0565       m_truthtrackpt = 0.;
0566       m_truthtrackphi = 0.;
0567       m_truthtracketa = 0.;
0568       m_truthtrackpid = 0;
0569     }
0570
0571     m_tracktree->Fill();
0572   }
0573
0574   ROOT::Math::XYZTVector v1, v2;
0575
0576   // make pairs
0577   for (auto iter1 = _trackmap->begin(); iter1 != _trackmap->end(); iter1++)
0578   {
0579     for (auto iter2 = iter1; iter2 != _trackmap->end(); iter2++)
0580     {
0581       if ( iter2 == iter1 ) continue;
0582
0583       //SvtxTrack *track2 = iter2.second;
0584       //std::cout << "XXX " << iter1->first << "\t" << iter2->first << std::endl;
0585       SvtxTrack *track1 = iter1->second;
0586       SvtxTrack *track2 = iter2->second;
0587
0588       // same sign or opposite
0589       m_pq1 = track1->get_charge();
0590       m_pq2 = track2->get_charge();
0591       //std::cout << "charge " << m_pq1 << "\t" << m_pq2 << std::endl;
0592       int type = 0;
0593       if ( m_pq1*m_pq2 > 0 )
0594       {
0595         type = 1;
0596       }
0597
0598       double px1 = track1->get_px();
0599       double py1 = track1->get_py();
0600       double pz1 = track1->get_pz();
0601       double e1 = sqrt( E_MASS*E_MASS + px1*px1 + py1*py1 + pz1*pz1 );
0602       v1.SetPxPyPzE( px1, py1, pz1, e1 );
0603
0604       double px2 = track2->get_px();
0605       double py2 = track2->get_py();
0606       double pz2 = track2->get_pz();
0607       double e2 = sqrt( E_MASS*E_MASS + px2*px2 + py2*py2 + pz2*pz2 );
0608       v2.SetPxPyPzE( px2, py2, pz2, e2 );
0609
0610       //TLorentzVector sum = v1 + v2;
0611       ROOT::Math::XYZTVector sum = v1 + v2;
0612       m_pm = sum.M();
0613       m_ppt = sum.Pt();
0614       m_pphi = sum.Phi();
0615       m_py = sum.Rapidity();
0616       m_peta = sum.Eta();
0617       //m_pdphi = ROOT::Math::VectorUtil::DeltaPhi(v1,v2);
0618       m_pdphi = ROOT::Math::VectorUtil::DeltaPhi(v1,v2);
0619       m_ppt1 = v1.Pt();
0620       m_ppz1 = v1.Pz();
0621       m_pphi1 = v1.Phi();
0622       m_peta1 = v1.Eta();
0623       m_ppt2 = v2.Pt();
0624       m_ppz2 = v2.Pz();
0625       m_pphi2 = v2.Phi();
0626       m_peta2 = v2.Eta();
0627
0628       h_mass[type]->Fill( m_pm );
0629       h_pt[type]->Fill( m_ppt );
0630       h_y[type]->Fill( m_py );
0631       h_eta[type]->Fill( m_peta );
0632       h2_eta_phi[type]->Fill( m_peta, m_pphi );
0633       h_phi[type]->Fill( m_pphi );
0634
0635       m_pairtree->Fill();
0636     }
0637   }
0638
0639   return Fun4AllReturnCodes::EVENT_OK;
0640 }
0641
0642
0643 /**
0644  * This method gets clusters from the EMCal and stores them in a tree. It
0645  * also demonstrates how to get trigger emulator information. Clusters from
0646  * other containers can be obtained in a similar way (e.g. clusters from
0647  * the IHCal, etc.)
0648  */
0649 /*
0650    void AnaUPC::getEMCalClusters(PHCompositeNode *topNode)
0651    {
0652 /// Get the raw cluster container
0653 /// Note: other cluster containers exist as well. Check out the node tree when
0654 /// you run a simulation
0655 RawClusterContainer *clusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_CEMC");
0656
0657 if (!clusters)
0658 {
0659 std::cout << PHWHERE
0660 << "EMCal cluster node is missing, can't collect EMCal clusters"
0661 << std::endl;
0662 return;
0663 }
0664
0665 /// Get the global vertex to determine the appropriate pseudorapidity of the clusters
0666 GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0667 if (!vertexmap)
0668 {
0669 std::cout << "AnaUPC::getEmcalClusters - Fatal Error - GlobalVertexMap node is missing. Please turn on the do_global flag in the main macro in order to reconstruct the global vertex." << std::endl;
0670 assert(vertexmap);  // force quit
0671
0672 return;
0673 }
0674
0675 if (vertexmap->empty())
0676 {
0677 std::cout << "AnaUPC::getEmcalClusters - Fatal Error - GlobalVertexMap node is empty. Please turn on the do_global flag in the main macro in order to reconstruct the global vertex." << std::endl;
0678 return;
0679 }
0680
0681 /// just take a bbc vertex
0682 GlobalVertex *vtx = nullptr;
0683 for(auto iter = vertexmap->begin(); iter!= vertexmap->end(); ++iter)
0684 {
0685 GlobalVertex* vertex = iter->second;
0686 if(vertex->find_vtxids(GlobalVertex::MBD) != vertex->end_vtxids())
0687 {
0688 vtx = vertex;
0689 }
0690 }
0691 if (vtx == nullptr)
0692 {
0693 return;
0694 }
0695
0696 /// Trigger emulator
0697 CaloTriggerInfo *trigger = findNode::getClass<CaloTriggerInfo>(topNode, "CaloTriggerInfo");
0698
0699 /// Can obtain some trigger information if desired
0700 if (trigger)
0701 {
0702 m_E_4x4 = trigger->get_best_EMCal_4x4_E();
0703 }
0704 RawClusterContainer::ConstRange begin_end = clusters->getClusters();
0705 RawClusterContainer::ConstIterator clusIter;
0706
0707 /// Loop over the EMCal clusters
0708 for (clusIter = begin_end.first; clusIter != begin_end.second; ++clusIter)
0709 {
0710 /// Get this cluster
0711 const RawCluster *cluster = clusIter->second;
0712
0713 /// Get cluster characteristics
0714 /// This helper class determines the photon characteristics
0715 /// depending on the vertex position
0716 /// This is important for e.g. eta determination and E_T determination
0717 CLHEP::Hep3Vector vertex(vtx->get_x(), vtx->get_y(), vtx->get_z());
0718 CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetECoreVec(*cluster, vertex);
0719 m_clusenergy = E_vec_cluster.mag();
0720 m_cluseta = E_vec_cluster.pseudoRapidity();
0721 m_clustheta = E_vec_cluster.getTheta();
0722 m_cluspt = E_vec_cluster.perp();
0723 m_clusphi = E_vec_cluster.getPhi();
0724
0725 if (m_cluspt < m_mincluspt)
0726 {
0727   continue;
0728 }
0729
0730 m_cluspx = m_cluspt * cos(m_clusphi);
0731 m_cluspy = m_cluspt * sin(m_clusphi);
0732 m_cluspz = sqrt(m_clusenergy * m_clusenergy - m_cluspx * m_cluspx - m_cluspy * m_cluspy);
0733
0734 // fill the cluster tree with all emcal clusters
0735 m_clustertree->Fill();
0736 }
0737 }
0738 */
0739
0740 /**
0741  * This function puts all of the tree branch assignments in one place so as to not
0742  * clutter up the AnaUPC::Init function.
0743  */
0744 void AnaUPC::initializeTrees()
0745 {
0746   m_tracktree = new TTree("tracktree", "A tree with svtx tracks");
0747   m_tracktree->Branch("m_tr_px", &m_tr_px, "m_tr_px/D");
0748   m_tracktree->Branch("m_tr_py", &m_tr_py, "m_tr_py/D");
0749   m_tracktree->Branch("m_tr_pz", &m_tr_pz, "m_tr_pz/D");
0750   m_tracktree->Branch("m_tr_p", &m_tr_p, "m_tr_p/D");
0751   m_tracktree->Branch("m_tr_pt", &m_tr_pt, "m_tr_pt/D");
0752   m_tracktree->Branch("m_tr_phi", &m_tr_phi, "m_tr_phi/D");
0753   m_tracktree->Branch("m_tr_eta", &m_tr_eta, "m_tr_eta/D");
0754   m_tracktree->Branch("m_charge", &m_charge, "m_charge/I");
0755   m_tracktree->Branch("m_chisq", &m_chisq, "m_chisq/D");
0756   m_tracktree->Branch("m_ndf", &m_ndf, "m_ndf/I");
0757   m_tracktree->Branch("m_dca", &m_dca, "m_dca/D");
0758   m_tracktree->Branch("m_tr_x", &m_tr_x, "m_tr_x/D");
0759   m_tracktree->Branch("m_tr_y", &m_tr_y, "m_tr_y/D");
0760   m_tracktree->Branch("m_tr_z", &m_tr_z, "m_tr_z/D");
0761   m_tracktree->Branch("m_truth_is_primary", &m_truth_is_primary, "m_truth_is_primary/I");
0762   m_tracktree->Branch("m_truthtrackpx", &m_truthtrackpx, "m_truthtrackpx/D");
0763   m_tracktree->Branch("m_truthtrackpy", &m_truthtrackpy, "m_truthtrackpy/D");
0764   m_tracktree->Branch("m_truthtrackpz", &m_truthtrackpz, "m_truthtrackpz/D");
0765   m_tracktree->Branch("m_truthtrackp", &m_truthtrackp, "m_truthtrackp/D");
0766   m_tracktree->Branch("m_truthtracke", &m_truthtracke, "m_truthtracke/D");
0767   m_tracktree->Branch("m_truthtrackpt", &m_truthtrackpt, "m_truthtrackpt/D");
0768   m_tracktree->Branch("m_truthtrackphi", &m_truthtrackphi, "m_truthtrackphi/D");
0769   m_tracktree->Branch("m_truthtracketa", &m_truthtracketa, "m_truthtracketa/D");
0770   m_tracktree->Branch("m_truthtrackpid", &m_truthtrackpid, "m_truthtrackpid/I");
0771
0772   m_globaltree = new TTree("globaltree", "Global Info");
0773   m_globaltree->Branch("run", &m_run, "run/I");
0774   m_globaltree->Branch("evt", &m_evt, "evt/I");
0775   m_globaltree->Branch("npart", &m_npart, "npart/I");
0776   m_globaltree->Branch("ncoll", &m_ncoll, "ncoll/I");
0777   m_globaltree->Branch("b", &m_bimpact, "b/F");
0778   m_globaltree->Branch("totntrks", &m_ntracks, "tottrks/I");
0779   m_globaltree->Branch("sphntrks", &m_ntrk_sphenix, "sphntrks/I");
0780   m_globaltree->Branch("hijntrks", &m_ntrk_mc, "mcntrks/I");
0781
0782   m_truthtree = new TTree("truthg4tree", "A tree with truth g4 particles");
0783   m_truthtree->Branch("evt", &m_evt, "evt/I");
0784   m_truthtree->Branch("m_e", &m_truthenergy, "m_truthe/D");
0785   m_truthtree->Branch("m_p", &m_truthp, "m_truthp/D");
0786   m_truthtree->Branch("m_px", &m_truthpx, "m_truthpx/D");
0787   m_truthtree->Branch("m_py", &m_truthpy, "m_truthpy/D");
0788   m_truthtree->Branch("m_pz", &m_truthpz, "m_truthpz/D");
0789   m_truthtree->Branch("m_pt", &m_truthpt, "m_truthpt/D");
0790   m_truthtree->Branch("m_phi", &m_truthphi, "m_truthphi/D");
0791   m_truthtree->Branch("m_eta", &m_trutheta, "m_trutheta/D");
0792   m_truthtree->Branch("m_pid", &m_truthpid, "m_truthpid/I");
0793   m_truthtree->Branch("m_q", &m_truthcharge, "m_truthcharge/I");
0794
0795   m_pairtree = new TTree("pairs", "opp sign pairs");
0796   m_pairtree->Branch("prun", &m_run, "prun/I");
0797   m_pairtree->Branch("pevt", &m_evt, "pevt/I");
0798   m_pairtree->Branch("pm", &m_pm, "pm/D");
0799   m_pairtree->Branch("ppt", &m_ppt, "ppt/D");
0800   m_pairtree->Branch("pphi", &m_pphi, "pphi/D");
0801   m_pairtree->Branch("py", &m_py, "py/D");
0802   m_pairtree->Branch("peta", &m_peta, "peta/D");
0803   m_pairtree->Branch("pdphi", &m_pdphi, "pdphi/D");
0804   m_pairtree->Branch("ppt1", &m_ppt1, "ppt1/D");
0805   m_pairtree->Branch("ppz1", &m_ppz1, "ppz1/D");
0806   m_pairtree->Branch("pphi1", &m_pphi1, "pphi1/D");
0807   m_pairtree->Branch("peta1", &m_peta1, "peta1/D");
0808   m_pairtree->Branch("ppt2", &m_ppt2, "ppt2/D");
0809   m_pairtree->Branch("ppz2", &m_ppz2, "ppz2/D");
0810   m_pairtree->Branch("pphi2", &m_pphi2, "pphi2/D");
0811   m_pairtree->Branch("peta2", &m_peta2, "peta2/D");
0812   m_pairtree->Branch("pq1", &m_pq1, "pq1/S");
0813   m_pairtree->Branch("pq2", &m_pq2, "pq2/S");
0814
0815   /*
0816      m_clustertree = new TTree("clustertree", "A tree with emcal clusters");
0817      m_clustertree->Branch("m_clusenergy", &m_clusenergy, "m_clusenergy/D");
0818      m_clustertree->Branch("m_cluseta", &m_cluseta, "m_cluseta/D");
0819      m_clustertree->Branch("m_clustheta", &m_clustheta, "m_clustheta/D");
0820      m_clustertree->Branch("m_cluspt", &m_cluspt, "m_cluspt/D");
0821      m_clustertree->Branch("m_clusphi", &m_clusphi, "m_clusphi/D");
0822      m_clustertree->Branch("m_cluspx", &m_cluspx, "m_cluspx/D");
0823      m_clustertree->Branch("m_cluspy", &m_cluspy, "m_cluspy/D");
0824      m_clustertree->Branch("m_cluspz", &m_cluspz, "m_cluspz/D");
0825      m_clustertree->Branch("m_E_4x4", &m_E_4x4, "m_E_4x4/D");
0826      */
0827 }
0828
0829 /**
0830  * This function initializes all of the member variables in this class so that there
0831  * are no variables that might not be set before e.g. writing them to the output
0832  * trees.
0833  */
0834 void AnaUPC::initializeVariables()
0835 {
0836   m_partid1 = -99;
0837   m_partid2 = -99;
0838   m_x1 = -99;
0839   m_x2 = -99;
0840   m_mpi = -99;
0841   m_process_id = -99;
0842   m_truthenergy = -99;
0843   m_trutheta = -99;
0844   m_truthphi = -99;
0845   m_truthp = -99;
0846   m_truthpx = -99;
0847   m_truthpy = -99;
0848   m_truthpz = -99;
0849   m_truthpt = -99;
0850   m_numparticlesinevent = -99;
0851   m_truthpid = -99;
0852
0853   m_tr_px = -99;
0854   m_tr_py = -99;
0855   m_tr_pz = -99;
0856   m_tr_p = -99;
0857   m_tr_pt = -99;
0858   m_tr_phi = -99;
0859   m_tr_eta = -99;
0860   m_charge = -99;
0861   m_chisq = -99;
0862   m_ndf = -99;
0863   m_dca = -99;
0864   m_tr_x = -99;
0865   m_tr_y = -99;
0866   m_tr_z = -99;
0867   m_truth_is_primary = -99;
0868   m_truthtrackpx = -99;
0869   m_truthtrackpy = -99;
0870   m_truthtrackpz = -99;
0871   m_truthtrackp = -99;
0872   m_truthtracke = -99;
0873   m_truthtrackpt = -99;
0874   m_truthtrackphi = -99;
0875   m_truthtracketa = -99;
0876   m_truthtrackpid = -99;
0877
0878   m_ntracks = 0;
0879   m_ntrk_sphenix = 0;
0880   m_ntrk_mc = 0;
0881
0882   /*
0883      m_recojetpt = -99;
0884      m_recojetid = -99;
0885      m_recojetpx = -99;
0886      m_recojetpy = -99;
0887      m_recojetpz = -99;
0888      m_recojetphi = -99;
0889      m_recojetp = -99;
0890      m_recojetenergy = -99;
0891      m_recojeteta = -99;
0892      m_truthjetid = -99;
0893      m_truthjetp = -99;
0894      m_truthjetphi = -99;
0895      m_truthjeteta = -99;
0896      m_truthjetpt = -99;
0897      m_truthjetenergy = -99;
0898      m_truthjetpx = -99;
0899      m_truthjetpy = -99;
0900      m_truthjetpz = -99;
0901      m_dR = -99;
0902      */
0903 }
0904
0905 int AnaUPC::Reset(PHCompositeNode * /*topNode*/)
0906 {
0907   //std::cout << "In Reset()" << std::endl;
0908   initializeVariables();
0909   return 0;
0910 }
0911