eic-smear/src/functions.cxx

0001 /**
0002  \file
0003  Global function implementations.
0004
0005  \author    Thomas Burton
0006  \date      2011-07-07
0007  \copyright 2011 Brookhaven National Lab
0008  */
0009
0010 #include "eicsmear/functions.h"
0011
0012 #include <cmath>
0013 #include <fstream>
0014 #include <set>
0015 #include <sstream>
0016 #include <string>
0017
0018 #include <TLorentzVector.h>
0019 #include <TMath.h>
0020 #include <TParticlePDG.h>
0021 #include <TVector2.h>
0022 #include <TVector3.h>
0023
0024 #include "eicsmear/erhic/EventMC.h"
0025 #include "eicsmear/erhic/ParticleMC.h"
0026
0027 /**
0028  Returns the first non-blank character in a line.
0029  Returns \0 if there are no non-blank characters in the line.
0030  */
0031 char getFirstNonBlank(const std::string& line) {
0032   char first('\0');  // NULL character by default
0033                       // Find the index first non-blank character.
0034   size_t index = line.find_first_not_of(" \t");
0035
0036   if (index != std::string::npos) {
0037     first = line.at(index);
0038   }  // if
0039
0040   return first;
0041 }
0042
0043 /**
0044  Calculate the hadron azimuthal angle around the virtual photon direction with
0045  respect to the lepton scattering plane in the proton rest frame.
0046  We use the HERMES convention, returning an angle in the range [0,2pi].
0047  The vectors passed as arguments should already be boosted to the proton rest
0048  frame. Incident and scattered leptons, incident protons and virtual photons
0049  all return -999.
0050  */
0051 double
0052 computeHermesPhiH(const TLorentzVector& hadronInPrf,
0053                   const TLorentzVector& leptonInPrf,
0054                   const TLorentzVector& photonInPrf) {
0055   const TVector3& q = photonInPrf.Vect();
0056   const TVector3& k = leptonInPrf.Vect();
0057   const TVector3& Ph = hadronInPrf.Vect();
0058   const TVector3& qCrossK = q.Cross(k);
0059   const TVector3& qCrossPh = q.Cross(Ph);
0060   // Used to determine spin direction:
0061   double qCrossKDotPh = qCrossK.Dot(Ph);
0062   if (::fabs(qCrossKDotPh) < 1.0e-5) return -999.0;
0063
0064   // First get phi in the range [0,pi]:
0065   double phi = TMath::ACos(qCrossK.Dot(qCrossPh)
0066                            / qCrossK.Mag()
0067                            / qCrossPh.Mag());
0068
0069   // Handle spin direction to get result in range [-pi,pi] by
0070   // scaling by [-1,+1], depending on hemisphere.
0071   phi *= qCrossKDotPh / ::fabs(qCrossKDotPh);
0072
0073   return TVector2::Phi_0_2pi(phi);
0074 }
0075
0076 //
0077 // class EventToDot.
0078 //
0079
0080 struct Pair {
0081   int a;
0082   int b;
0083   Pair(int x, int y) : a(x), b(y) { }
0084   // Less-than operator required for entry in sorted container.
0085   bool operator<(const Pair& other) const {
0086     if (other.a == a) return b < other.b;
0087     return a < other.a;
0088   }
0089 };
0090
0091
0092 void EventToDot::Generate(const erhic::EventMC& event,
0093                           const std::string& name) const {
0094   std::ofstream file(name.c_str());
0095   std::ostringstream oss;
0096
0097   file << "digraph G {" << std::endl;
0098   oss << "   label = \"Event " << event.GetN() << "\"';";
0099   file << oss.str() << std::endl;
0100
0101   std::set<Pair> pairs;
0102
0103   // Keep track of which particles we use in the diagram
0104   // so we can set node attributes at the end.
0105   // Use a set to avoid worrying about adding duplicates.
0106   std::set<const erhic::ParticleMC*> used;
0107
0108   // Loop over all tracks in the event and accumulate indices giving
0109   // parent/child relationships.
0110   // We look from parent->child and child->parent because they
0111   // aren't always fully indexed both ways.
0112   for (unsigned i(0); i < event.GetNTracks(); ++i) {
0113     // Check parent of particle
0114     const erhic::ParticleMC* parent = event.GetTrack(i)->GetParent();
0115     if (parent) {
0116       pairs.insert(Pair(parent->GetIndex(), event.GetTrack(i)->GetIndex()));
0117       used.insert(parent);
0118       used.insert(event.GetTrack(i));
0119     }  // if
0120     // Check children of particle
0121     for (unsigned j(0); j < event.GetTrack(i)->GetNChildren(); ++j) {
0122       pairs.insert(Pair(event.GetTrack(i)->GetIndex(),
0123                         event.GetTrack(i)->GetChild(j)->GetIndex()));
0124       used.insert(event.GetTrack(i));
0125       used.insert(event.GetTrack(i)->GetChild(j));
0126     }  // for
0127   }  // for
0128   // Insert relationships between particles.
0129   for (std::set<Pair>::const_iterator i = pairs.begin();
0130       i != pairs.end();
0131       ++i) {
0132     const erhic::ParticleMC* a = event.GetTrack(i->a - 1);
0133     const erhic::ParticleMC* b = event.GetTrack(i->b - 1);
0134     oss.str("");
0135     oss << "   " << a->GetIndex() << " -> " <<
0136     b->GetIndex();
0137     file << oss.str() << std::endl;
0138   }  // for
0139   file << "#   Node attributes" << std::endl;
0140   // Apply labels, formatting to used particles.
0141   for (std::set<const erhic::ParticleMC*>::const_iterator i = used.begin();
0142       i != used.end();
0143       ++i) {
0144     std::string shape("ellipse");
0145     if ((*i)->GetStatus() == 1) {
0146       shape = "box";
0147     }  // if
0148     if ((*i)->GetIndex() < 3) {
0149       shape = "diamond";
0150     }  // if
0151     oss.str("");
0152     oss << "   " << (*i)->GetIndex() << " [label=\"" << (*i)->GetIndex() << " ";
0153     if ((*i)->Id().Info()) {
0154       oss << (*i)->Id().Info()->GetName();
0155     } else {
0156       oss << "unknown PDG " << (*i)->Id().Code();
0157     }  // if
0158     oss << "\", shape=" << shape << "];";
0159     file << oss.str() << std::endl;
0160   }  // for
0161   file << "}";
0162 }