pythiaeRHIC/erhic/factory.cpp

0001 /**
0002  \file factory.cpp
0003
0004  \author Thomas Burton
0005  \date 10/9/12
0006  \copyright 2012 BNL.
0007  */
0008
0009 #include "factory.h"
0010
0011 #include <TBranch.h>
0012 #include <TLorentzVector.h>
0013 #include <TTree.h>
0014 #include <TVector3.h>
0015
0016 #include <eicsmear/erhic/EventPythia.h>
0017 #include <eicsmear/erhic/Kinematics.h>
0018
0019 #include "pythia_commons.h"
0020 #include "pythia_erhic.h"
0021
0022 // =====================================================================
0023 // =====================================================================
0024 Factory::Factory()
0025 : mEvent(new erhic::EventPythia()) {
0026 }
0027
0028 // =====================================================================
0029 // =====================================================================
0030 Factory::~Factory() {
0031    if(mEvent) {
0032       delete mEvent;
0033       mEvent = NULL;
0034    } // if
0035 }
0036
0037 // =====================================================================
0038 // Populate the stored event with the current contents of
0039 // PYTHIA and return a pointer to the event.
0040 // Do not delete this pointer!
0041 // =====================================================================
0042 erhic::EventPythia* Factory::Create() {
0043    mEvent->Clear("");
0044    // See the PYTHIA manual for more details about the meaning
0045    // of the various msti(), pari(), vint() variables.
0046    mEvent->SetGenEvent(__pythia6_MOD_genevent);
0047    mEvent->SetProcess(msti(1));
0048    mEvent->SetNucleon(msti(12));
0049    mEvent->SetTargetParton(msti(16));
0050    mEvent->SetBeamParton(msti(15));
0051    mEvent->SetTargetPartonX(pari(34));
0052    mEvent->SetBeamPartonX(pari(33));
0053    mEvent->SetBeamPartonTheta(pari(53));
0054    mEvent->SetLeptonPhi(vint(313));
0055    mEvent->SetF1(__pythia6_MOD_f1);
0056    mEvent->SetF2(__pythia6_MOD_f2);
0057    mEvent->SetSigmaRad(__pythia6_MOD_sigma_rad);
0058    mEvent->SetHardS(pari(14));
0059    mEvent->SetHardT(pari(15));
0060    mEvent->SetHardU(pari(16));
0061    mEvent->SetHardQ2(pari(22));
0062    mEvent->SetHardPt2(pari(18));
0063    mEvent->SetSigRadCor(__pythia6_MOD_sigradcor);
0064    mEvent->SetEBrems(__pythia6_MOD_ebrems);
0065    mEvent->SetPhotonFlux(vint(319));
0066    mEvent->SetTrueY(vint(309));
0067    mEvent->SetTrueQ2(vint(307));
0068    mEvent->SetTrueX(__pythia6_MOD_truex);
0069    mEvent->SetTrueW2(__pythia6_MOD_truew2);
0070    mEvent->SetTrueNu(__pythia6_MOD_truenu);
0071    mEvent->SetR(__pythia6_MOD_r);
0072    static std::stringstream stream;
0073    const int ntracks = pyjets_.n;
0074    // Loop with indices from [1, N] as per Fortran convention.
0075    // Accessor functions k(), p(), v() take care of correct indexing
0076    // of the equivalent C++ arrays.
0077    for(int i(1); i <= ntracks; ++i) {
0078       erhic::ParticleMC* track = new erhic::ParticleMC;
0079       track->SetIndex(i);
0080       track->SetStatus(k(i, 1));
0081       track->SetId(k(i, 2));
0082       track->SetParentIndex(k(i, 3));
0083       track->SetChild1Index(k(i, 4));
0084       track->SetChildNIndex(k(i, 5));
0085       track->Set4Vector(TLorentzVector(p(i, 1), p(i, 2),
0086                                        p(i, 3), p(i, 4)));
0087       track->SetM(p(i, 5));
0088       track->SetVertex(TVector3(v(i, 1), v(i, 2), v(i, 3)));
0089       track->SetEvent(mEvent);
0090       mEvent->AddLast(track);
0091    } // for
0092    // If a radiative photon was generated append it to the track record.
0093    if(__pythia6_MOD_ebrems > 0.) {
0094       erhic::ParticleMC* track = new erhic::ParticleMC;
0095       track->SetIndex(mEvent->GetNTracks() + 1);
0096       track->SetStatus(55);
0097       track->SetId(22);
0098       track->SetParentIndex(1);
0099       track->SetChild1Index(0);
0100       track->SetChildNIndex(0);
0101       track->Set4Vector(TLorentzVector(__pythia6_MOD_radgamp[0],
0102                                        __pythia6_MOD_radgamp[1],
0103                                        __pythia6_MOD_radgamp[2],
0104                                        __pythia6_MOD_radgame));
0105       track->SetVertex(TVector3(0., 0., 0.));
0106       track->SetEvent(mEvent);
0107       mEvent->AddLast(track);
0108    } // if
0109    // After adding all tracks, compute event kinematics.
0110    // Compute using all three methods: electron, Jacquet-Blondel
0111    // and double angle.
0112    erhic::DisKinematics* nm =
0113    erhic::LeptonKinematicsComputer(*mEvent).Calculate();
0114    erhic::DisKinematics* jb =
0115    erhic::JacquetBlondelComputer(*mEvent).Calculate();
0116    erhic::DisKinematics* da =
0117    erhic::DoubleAngleComputer(*mEvent).Calculate();
0118    if(nm) {
0119       mEvent->SetLeptonKinematics(*nm);
0120       delete nm;
0121    } // if
0122    if(jb) {
0123       mEvent->SetJacquetBlondelKinematics(*jb);
0124       delete jb;
0125    } // if
0126    if(da) {
0127       mEvent->SetDoubleAngleKinematics(*da);
0128       delete da;
0129    } // if
0130    // Set particle properties that depend on other particles.
0131    // This has to come last as some of these properties depend
0132    // on the computed event kinematics.
0133    for(unsigned i(0); i < mEvent->GetNTracks(); ++i) {
0134       mEvent->GetTrack(i)->ComputeEventDependentQuantities(*mEvent);
0135    } // for
0136    return mEvent;
0137 }
0138
0139 // =====================================================================
0140 // Returns a string with the full (including namespace) class name
0141 // of the event type produced.
0142 // This is important for use with ROOT TTree to ensure the correct
0143 // event type in branches.
0144 // =====================================================================
0145 std::string Factory::EventName() const {
0146    return mEvent->Class()->GetName();
0147 }
0148
0149 // =====================================================================
0150 // Add a branch named "name" for the event type generated
0151 // by this factory to a ROOT TTree.
0152 // Returns a pointer to the branch, or NULL in the case of an error.
0153 // =====================================================================
0154 TBranch* Factory::Branch(TTree& tree, const std::string& name) {
0155    return tree.Branch(name.c_str(), &mEvent);
0156 }