eic-smear  1.0.3
A collection of ROOT classes for Monte Carlo events and a fast-smearing code simulating detector effects for the Electron-Ion Collider task force
Pythia6EventBuilder.cxx
Go to the documentation of this file.
1 
11 
12 #include <cmath>
13 #include <iostream>
14 #include <memory>
15 #include <string>
16 
17 #include <TBranch.h>
18 #include <TClass.h>
19 #include <TMCParticle.h>
20 #include <TObjArray.h>
21 #include <TProcessID.h>
22 #include <TPythia6.h>
23 #include <TTree.h>
24 
30 
31 namespace erhic {
32 
34 : mNGenerated(0)
35 , mNTrials(0)
36 , mFilter(filter)
37 , mEvent(NULL) {
38 }
39 
41  delete mFilter;
42  mFilter = NULL;
43 }
44 
46  std::auto_ptr<EventPythia> event(BuildEvent());
47  if (mFilter) {
48  while (!mFilter->Accept(*event)) {
49  event.reset(BuildEvent());
50  } // while
51  } // if
52  return event.release();
53 }
54 
56  // Save current object count
57  int objectNumber = TProcessID::GetObjectCount();
58  // Generate a new PYTHIA event
59  TPythia6* pythia = TPythia6::Instance();
60  pythia->GenerateEvent();
61  // Import all particles (not just final-state)
62  TObjArray* particles = pythia->ImportParticles("All");
63  // Construct the EventPythia object from the current
64  // state of TPythia6.
65  std::auto_ptr<EventPythia> event(new EventPythia);
66  // Extract the event-wise quantities:
67  event->SetNucleon(pythia->GetMSTI(12));
68  event->SetTargetParton(pythia->GetMSTI(16));
69  event->SetBeamParton(pythia->GetMSTI(15));
70  event->SetGenEvent(1);
71  event->SetTargetPartonX(pythia->GetPARI(34));
72  event->SetBeamPartonX(pythia->GetPARI(33));
73  event->SetBeamPartonTheta(pythia->GetPARI(53));
74  event->SetLeptonPhi(pythia->GetVINT(313));
75  event->SetHardS(pythia->GetPARI(14));
76  event->SetHardT(pythia->GetPARI(15));
77  event->SetHardU(pythia->GetPARI(16));
78  event->SetHardQ2(pythia->GetPARI(22));
79  event->SetHardPt2(pythia->GetPARI(18));
80  event->SetPhotonFlux(pythia->GetVINT(319));
81  event->SetProcess(pythia->GetMSTI(1));
82  // We need the beam energies to compute the true x, W2 and nu.
83  // The beam lepton should be the first particle
84  // and the beam hadron the second particle.
85  const double eLepton =
86  static_cast<TMCParticle*>(particles->At(0))->GetEnergy();
87  const double eHadron =
88  static_cast<TMCParticle*>(particles->At(1))->GetEnergy();
89  const double mHadron =
90  static_cast<TMCParticle*>(particles->At(1))->GetMass();
91  // x, W2 and nu are calculated from y, Q2 and the beam energies.
92  // y and Q2 come from PYTHIA.
93  // Use (approximate expression) Q2 = sxy, where s = 4.E_e.E_h
94  double y = pythia->GetVINT(309);
95  double Q2 = pythia->GetVINT(307);
96  double x = Q2 / y / 4. / eLepton / eHadron;
97  double W2 = pow(mHadron, 2.) + Q2 * (1. / x - 1.);
98  double nu = (W2 + Q2 - pow(mHadron, 2.)) / 2. / mHadron;
99  event->SetTrueY(y);
100  event->SetTrueQ2(Q2);
101  event->SetTrueX(x);
102  event->SetTrueW2(W2);
103  event->SetTrueNu(nu);
104  // Set the number of event generation trials taken since the last call
105  event->SetN(pythia->GetMSTI(5));
106  event->SetGenEvent(pythia->GetPyint5()->NGEN[2][0] - mNTrials);
107  mNTrials = pythia->GetPyint5()->NGEN[2][0];
108  // Now populate the particle list.
109  Pythia6ParticleBuilder builder;
110  for (int i(0); i < particles->GetEntries(); ++i) {
111  TMCParticle* p =
112  static_cast<TMCParticle*>(particles->At(i));
113  std::auto_ptr<ParticleMC> particle = builder.Create(*p);
114  particle->SetIndex(i + 1);
115  particle->SetEvent(event.get());
116  event->AddLast(particle.get());
117  } // for
118  // Compute derived event kinematics
122  if (nm) {
123  event->SetLeptonKinematics(*nm);
124  } // if
125  if (jb) {
126  event->SetJacquetBlondelKinematics(*jb);
127  } // if
128  if (da) {
129  event->SetDoubleAngleKinematics(*da);
130  } // if
131  // We also have to set the remaining variables not taken care of
132  // by the general DIS event kinematic computations.
133  // Find the beams, exchange boson, scattered lepton.
134  BeamParticles beams;
135  if (ParticleIdentifier::IdentifyBeams(*event, beams)) {
136  const TLorentzVector h = beams.BeamHadron();
137  TLorentzVector l = beams.BeamLepton();
138  TLorentzVector s = beams.ScatteredLepton();
139  TVector3 boost = -h.BoostVector();
140  l.Boost(boost);
141  s.Boost(boost);
142  event->SetELeptonInNuclearFrame(l.E());
143  event->SetEScatteredInNuclearFrame(s.E());
144  } // if
145  for (unsigned i(0); i < event->GetNTracks(); ++i) {
146  event->GetTrack(i)->ComputeEventDependentQuantities(*event);
147  } // for
148  // Restore Object count
149  // See example in $ROOTSYS/test/Event.cxx
150  // To save space in the table keeping track of all referenced objects
151  // we assume that our events do not address each other. We reset the
152  // object count to what it was at the beginning of the event.
153  TProcessID::SetObjectCount(objectNumber);
154  return event.release();
155 }
156 
157 std::string Pythia6EventBuilder::EventName() const {
158  return EventPythia::Class()->GetName();
159 }
160 
161 TBranch* Pythia6EventBuilder::Branch(TTree& tree, const std::string& name) {
162  EventPythia* event(NULL);
163  TBranch* branch =
164  tree.Branch(name.c_str(), EventName().c_str(), &event, 32000, 99);
165  tree.ResetBranchAddress(branch);
166  return branch;
167 }
168 
169 void Pythia6EventBuilder::Fill(TBranch& branch) {
170  if (mEvent) {
171  branch.ResetAddress();
172  delete mEvent;
173  mEvent = NULL;
174  } // if
175  mEvent = Create();
176  branch.SetAddress(&mEvent);
177 }
178 
179 } // namespace erhic
erhic::Pythia6EventBuilder::mFilter
EventMCFilterABC * mFilter
Definition: Pythia6EventBuilder.h:65
Pythia6EventBuilder.h
erhic
Definition: EventDis.cxx:14
erhic::DisKinematics
A collection of DIS kinematic variables.
Definition: Kinematics.h:31
EventPythia.h
tree
Definition: tree.py:1
BeamParticles::ScatteredLepton
const TLorentzVector & ScatteredLepton() const
Definition: BeamParticles.h:89
erhic::LeptonKinematicsComputer
Computes DIS event kinematics from the scattered lepton.
Definition: Kinematics.h:64
erhic::JacquetBlondelComputer
Computes DIS event kinematics from final-state hadrons using the Jacquet-Blondel method.
Definition: Kinematics.h:86
erhic::Pythia6EventBuilder::mEvent
EventPythia * mEvent
Definition: Pythia6EventBuilder.h:66
erhic::Pythia6EventBuilder::Pythia6EventBuilder
Pythia6EventBuilder(EventMCFilterABC *=NULL)
Constructor.
Definition: Pythia6EventBuilder.cxx:33
erhic::Pythia6EventBuilder::EventName
virtual std::string EventName() const
Returns a string with the full (including namespace) class name of the event type produced.
Definition: Pythia6EventBuilder.cxx:157
ParticleIdentifier::IdentifyBeams
static bool IdentifyBeams(const erhic::VirtualEvent &, BeamParticles &)
Identify the beams from an event and store their properties in a BeamParticles object.
Definition: ParticleIdentifier.cxx:139
erhic::Pythia6ParticleBuilder::Create
std::auto_ptr< ParticleMC > Create(const TMCParticle &) const
Generate a ParticleMC from a ROOT TMCParticle.
Definition: Pythia6ParticleBuilder.cxx:22
erhic::DoubleAngleComputer::Calculate
virtual DisKinematics * Calculate()
Definition: Kinematics.cxx:455
erhic::EventMCFilterABC
Abstract base class for a filter on Monte Carlo events.
Definition: EventMCFilterABC.h:22
EventMCFilterABC.h
erhic::EventPythia
Pythia 6 DIS event.
Definition: erhic/EventPythia.h:28
erhic::Pythia6EventBuilder::BuildEvent
EventPythia * BuildEvent()
Definition: Pythia6EventBuilder.cxx:55
erhic::Pythia6EventBuilder::Branch
virtual TBranch * Branch(TTree &, const std::string &)
Add a branch named "name" for the event type generated by this factory to a ROOT TTree.
Definition: Pythia6EventBuilder.cxx:161
erhic::Pythia6ParticleBuilder
Factory class for Monte Carlo particles.
Definition: Pythia6ParticleBuilder.h:24
BeamParticles::BeamLepton
const TLorentzVector & BeamLepton() const
Definition: BeamParticles.h:81
erhic::EventMCFilterABC::Accept
virtual bool Accept(const VirtualEvent &) const =0
Implement this method in a derived class such that it returns true when the event passes the filter's...
erhic::Pythia6EventBuilder::Create
virtual EventPythia * Create()
Generates an event from the current state of ROOT::TPythia6.
Definition: Pythia6EventBuilder.cxx:45
BeamParticles::BeamHadron
const TLorentzVector & BeamHadron() const
Definition: BeamParticles.h:77
erhic::Pythia6EventBuilder::Fill
virtual void Fill(TBranch &)
Calls Create() to generate an event and fills the provided branch with that event.
Definition: Pythia6EventBuilder.cxx:169
erhic::LeptonKinematicsComputer::Calculate
virtual DisKinematics * Calculate()
Definition: Kinematics.cxx:253
erhic::JacquetBlondelComputer::Calculate
virtual DisKinematics * Calculate()
Definition: Kinematics.cxx:332
Pythia6ParticleBuilder.h
erhic::Pythia6EventBuilder::~Pythia6EventBuilder
virtual ~Pythia6EventBuilder()
Destructor.
Definition: Pythia6EventBuilder.cxx:40
erhic::DoubleAngleComputer
Computes DIS event kinematics from a mixture of hadronic and lepton variables using the double-angle ...
Definition: Kinematics.h:117
ParticleMC.h
BeamParticles
Wrapper class around energy-momentum 4-vectors defining the incident and scattered beams and the exch...
Definition: BeamParticles.h:20
erhic::Pythia6EventBuilder::mNTrials
Long64_t mNTrials
Definition: Pythia6EventBuilder.h:64
Kinematics.h