hbcup/chiu/UPCMeson.cc

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