src/erhic/Pythia6EventBuilder.cxx

0001 /**
0002  \file
0003  Implementation of class erhic::Pythia6EventBuilder.
0004
0005  \author    Thomas Burton
0006  \date      2012-01-17
0007  \copyright 2012 Brookhaven National Lab
0008  */
0009
0010 #include "eicsmear/erhic/Pythia6ParticleBuilder.h"
0011
0012 #include <cmath>
0013 #include <iostream>
0014 #include <memory>
0015 #include <string>
0016
0017 #include <TBranch.h>
0018 #include <TClass.h>
0019 #include <TMCParticle.h>
0020 #include <TObjArray.h>
0021 #include <TProcessID.h>
0022 #include <TPythia6.h>
0023 #include <TTree.h>
0024
0025 #include "eicsmear/erhic/EventMCFilterABC.h"
0026 #include "eicsmear/erhic/EventPythia.h"
0027 #include "eicsmear/erhic/Kinematics.h"
0028 #include "eicsmear/erhic/ParticleMC.h"
0029 #include "eicsmear/erhic/Pythia6EventBuilder.h"
0030
0031 namespace erhic {
0032
0033 Pythia6EventBuilder::Pythia6EventBuilder(EventMCFilterABC* filter)
0034 : mNGenerated(0)
0035 , mNTrials(0)
0036 , mFilter(filter)
0037 , mEvent(NULL) {
0038 }
0039
0040 Pythia6EventBuilder::~Pythia6EventBuilder() {
0041   delete mFilter;
0042   mFilter = NULL;
0043 }
0044
0045 EventPythia* Pythia6EventBuilder::Create() {
0046   std::unique_ptr<EventPythia> event(BuildEvent());
0047   if (mFilter) {
0048     while (!mFilter->Accept(*event)) {
0049       event.reset(BuildEvent());
0050     }  // while
0051   }  // if
0052   return event.release();
0053 }
0054
0055 EventPythia* Pythia6EventBuilder::BuildEvent() {
0056   // Save current object count
0057   int objectNumber = TProcessID::GetObjectCount();
0058   // Generate a new PYTHIA event
0059   TPythia6* pythia = TPythia6::Instance();
0060   pythia->GenerateEvent();
0061   // Import all particles (not just final-state)
0062   TObjArray* particles = pythia->ImportParticles("All");
0063   // Construct the EventPythia object from the current
0064   // state of TPythia6.
0065   std::unique_ptr<EventPythia> event(new EventPythia);
0066   // Extract the event-wise quantities:
0067   event->SetNucleon(pythia->GetMSTI(12));
0068   event->SetTargetParton(pythia->GetMSTI(16));
0069   event->SetBeamParton(pythia->GetMSTI(15));
0070   event->SetGenEvent(1);
0071   event->SetTargetPartonX(pythia->GetPARI(34));
0072   event->SetBeamPartonX(pythia->GetPARI(33));
0073   event->SetBeamPartonTheta(pythia->GetPARI(53));
0074   event->SetLeptonPhi(pythia->GetVINT(313));
0075   event->SetHardS(pythia->GetPARI(14));
0076   event->SetHardT(pythia->GetPARI(15));
0077   event->SetHardU(pythia->GetPARI(16));
0078   event->SetHardQ2(pythia->GetPARI(22));
0079   event->SetHardPt2(pythia->GetPARI(18));
0080   event->SetPhotonFlux(pythia->GetVINT(319));
0081   event->SetProcess(pythia->GetMSTI(1));
0082   // We need the beam energies to compute the true x, W2 and nu.
0083   // The beam lepton should be the first particle
0084   // and the beam hadron the second particle.
0085   const double eLepton =
0086   static_cast<TMCParticle*>(particles->At(0))->GetEnergy();
0087   const double eHadron =
0088   static_cast<TMCParticle*>(particles->At(1))->GetEnergy();
0089   const double mHadron =
0090   static_cast<TMCParticle*>(particles->At(1))->GetMass();
0091   // x, W2 and nu are calculated from y, Q2 and the beam energies.
0092   // y and Q2 come from PYTHIA.
0093   // Use (approximate expression) Q2 = sxy, where s = 4.E_e.E_h
0094   double y = pythia->GetVINT(309);
0095   double Q2 = pythia->GetVINT(307);
0096   double x = Q2 / y / 4. / eLepton / eHadron;
0097   double W2 = pow(mHadron, 2.) + Q2 * (1. / x - 1.);
0098   double nu = (W2 + Q2 - pow(mHadron, 2.)) / 2. / mHadron;
0099   event->SetTrueY(y);
0100   event->SetTrueQ2(Q2);
0101   event->SetTrueX(x);
0102   event->SetTrueW2(W2);
0103   event->SetTrueNu(nu);
0104   // Set the number of event generation trials taken since the last call
0105   event->SetN(pythia->GetMSTI(5));
0106   event->SetGenEvent(pythia->GetPyint5()->NGEN[2][0] - mNTrials);
0107   mNTrials = pythia->GetPyint5()->NGEN[2][0];
0108   // Now populate the particle list.
0109   Pythia6ParticleBuilder builder;
0110   for (int i(0); i < particles->GetEntries(); ++i) {
0111     TMCParticle* p =
0112     static_cast<TMCParticle*>(particles->At(i));
0113     std::unique_ptr<ParticleMC> particle = builder.Create(*p);
0114     particle->SetIndex(i + 1);
0115     particle->SetEvent(event.get());
0116     event->AddLast(particle.get());
0117   }  // for
0118   // Compute derived event kinematics
0119   DisKinematics* nm = LeptonKinematicsComputer(*event).Calculate();
0120   DisKinematics* jb = JacquetBlondelComputer(*event).Calculate();
0121   DisKinematics* da = DoubleAngleComputer(*event).Calculate();
0122   if (nm) {
0123     event->SetLeptonKinematics(*nm);
0124   }  // if
0125   if (jb) {
0126     event->SetJacquetBlondelKinematics(*jb);
0127   }  // if
0128   if (da) {
0129     event->SetDoubleAngleKinematics(*da);
0130   }  // if
0131   // We also have to set the remaining variables not taken care of
0132   // by the general DIS event kinematic computations.
0133   // Find the beams, exchange boson, scattered lepton.
0134   BeamParticles beams;
0135   if (ParticleIdentifier::IdentifyBeams(*event, beams)) {
0136     const TLorentzVector h = beams.BeamHadron();
0137     TLorentzVector l = beams.BeamLepton();
0138     TLorentzVector s = beams.ScatteredLepton();
0139     TVector3 boost = -h.BoostVector();
0140     l.Boost(boost);
0141     s.Boost(boost);
0142     event->SetELeptonInNuclearFrame(l.E());
0143     event->SetEScatteredInNuclearFrame(s.E());
0144   }  // if
0145   for (unsigned i(0); i < event->GetNTracks(); ++i) {
0146     event->GetTrack(i)->ComputeEventDependentQuantities(*event);
0147   }  // for
0148   // Restore Object count
0149   // See example in $ROOTSYS/test/Event.cxx
0150   // To save space in the table keeping track of all referenced objects
0151   // we assume that our events do not address each other. We reset the
0152   // object count to what it was at the beginning of the event.
0153   TProcessID::SetObjectCount(objectNumber);
0154   return event.release();
0155 }
0156
0157 std::string Pythia6EventBuilder::EventName() const {
0158   return EventPythia::Class()->GetName();
0159 }
0160
0161 TBranch* Pythia6EventBuilder::Branch(TTree& tree, const std::string& name) {
0162   EventPythia* event(NULL);
0163   TBranch* branch =
0164   tree.Branch(name.c_str(), EventName().c_str(), &event, 32000, 99);
0165   tree.ResetBranchAddress(branch);
0166   return branch;
0167 }
0168
0169 void Pythia6EventBuilder::Fill(TBranch& branch) {
0170   if (mEvent) {
0171     branch.ResetAddress();
0172     delete mEvent;
0173     mEvent = NULL;
0174   }  // if
0175   mEvent = Create();
0176   branch.SetAddress(&mEvent);
0177 }
0178
0179 }  // namespace erhic