Initial Commit of BaryoGEN, forked from sphaleronMC

.gitignore

@@ -0,0 +1,5 @@
+MyDict.cxx
+MyDict_rdict.pcm
+genTree
+*.root
+*.lhe

Makefile

@@ -0,0 +1,9 @@
+CXXFLAGS=$(shell root-config --cflags)
+GLIBS=$(shell root-config --glibs)
+
+SRCS:=$(wildcard src/*.cc)
+
+all: BaryoGEN
+
+BaryoGEN: $(SRCS)
+	g++ $(CXXFLAGS) $^ -o $@ -I$(HOME)/local/include/ $(GLIBS) `lhapdf-config --cflags --ldflags`

README.md

@@ -0,0 +1,19 @@
+# BaryoGEN
+Compiles with gcc version >= 5.3.0
+Requires ROOT version >= 6.06/01
+Requires LHAPDF == 6.1.6
+
+mstwpdf is from the mstw2008 package, available publically at https://mstwpdf.hepforge.org/
+    A. D. Martin, W. J. Stirling, R. S. Thorne and G. Watt,
+    "Parton distributions for the LHC",
+    Eur. Phys. J. C63 (2009) 189-285
+    [arXiv:0901.0002 [hep-ph]].
+
+The LHAPDF code which can be found at:
+    "LHAPDF6: parton density access in the LHC precision era"
+    Eur.Phys.J. C75 (2015) 3, 132
+    http://arxiv.org/abs/1412.7420
+
+The pdf grid string is hard coded, and can be changed before compiling.
+End user must make sure they have downloaded the desired pdf grid and
+installed it into their installation of LHAPDF.

include/BaryoGEN.h

@@ -0,0 +1,2 @@
+#pragma link C++ class std::vector < TLorentzVector >+;
+#pragma link C++ class std::vector < std::vector<TLorentzVector> >+;

include/LHEWriter.h

@@ -0,0 +1,22 @@
+#include <fstream>
+#include <vector>
+#include <string>
+
+#include "particles.h"
+#include "TLorentzVector.h"
+
+using namespace std;
+
+class LHEWriter
+{
+    private:
+        ofstream oF;
+
+    public:
+        LHEWriter(string fName);
+        ~LHEWriter();
+
+        int writeEvent(vector<particle> outParts, double Q);
+        int close();
+
+};

include/LinkDef.h

@@ -0,0 +1,3 @@
+#pragma link C++ class std::vector <TLorentzVector> +;
+#pragma link C++ class std::vector < std::vector<TLorentzVector> >+;
+

include/configBuilder.h

@@ -0,0 +1,22 @@
+#include "particles.h"
+
+#include "TRandom3.h"
+
+#include <vector>
+
+using namespace std;
+
+class configBuilder
+{
+    private:
+        vector<particle> parts;
+        TRandom3 rand;
+
+    public:
+        configBuilder();
+        ~configBuilder();
+
+        vector<particle> build(int iQ1, int iQ2, float pNCS, bool bCancel);
+
+};
+

include/decay.h

@@ -0,0 +1,123 @@
+#include <iostream>
+#include <cmath>
+
+#include "TRandom3.h"
+#include "TLorentzVector.h"
+#include "TVector3.h"
+
+using namespace std;
+
+struct biP 
+{
+    TLorentzVector p1;
+    TLorentzVector p2;
+};
+
+struct triP 
+{
+    TLorentzVector p1;
+    TLorentzVector p2;
+    TLorentzVector p3;
+};
+
+biP decayTwo(TLorentzVector mother,double m1, double m2)
+{
+    TRandom3 rand;
+    rand.SetSeed();
+
+    //Get boost vector
+    TVector3 boost;
+    boost = mother.BoostVector();
+    double mass = mother.M();
+
+    if(m1 + m2 > mass) 
+    {
+        cout << "m1: " << m1 << endl;
+        cout << "m2: " << m2 << endl;
+        cout << "mass: " << mass << endl;
+        throw std::invalid_argument( "Daughter mass sum must be less than Mother mass" );
+    }
+
+    //Decay in random direction
+    double th1 = acos(2.0*rand.Uniform()-1.0);
+    double phi1 = 2.0*M_PI*rand.Uniform() - M_PI;
+
+    //Solve for everything else
+    double E1 = (mass*mass - m2*m2 + m1*m1)/(2.0*mass);
+    double E2 = mass - E1;
+    double mp1 = sqrt(E1*E1 - m1*m1);
+
+    TVector3 p1(mp1*cos(phi1)*sin(th1),mp1*sin(phi1)*sin(th1),mp1*cos(th1));
+    TVector3 p2 = -1.0*p1;
+
+    biP output;
+    output.p1.SetPxPyPzE(p1.Px(),p1.Py(),p1.Pz(),E1);
+    output.p2.SetPxPyPzE(p2.Px(),p2.Py(),p2.Pz(),E2);
+
+    if(abs(output.p1.Eta()) > 100 )
+    {
+        cout << "mass: " << mass << endl;
+        cout << "th1: " << th1 << endl;
+        cout << "E1: " << E1 << endl;
+        cout << "phi1: " << phi1 << endl;
+        cout << "mp1: " << mp1 << endl;
+    }
+
+    //Boost into lab frame
+    output.p1.Boost(boost);
+    output.p2.Boost(boost);
+
+    return output;
+
+
+};
+
+
+triP decayThree(TLorentzVector mother,double m1, double m2,double m3)
+{
+    TRandom3 rand;
+    rand.SetSeed();
+
+    //Get boost vector
+    TVector3 boost;
+    boost = mother.BoostVector();
+    double mass = mother.M();
+
+    if(m1 + m2 + m3 > mass) throw std::invalid_argument( "Daughter mass sum must be less than Mother mass" );
+
+    //Decay in random direction
+    double th1 = acos(2.0*rand.Uniform()-1.0);
+    double phi1 = 2.0*M_PI*rand.Uniform() - M_PI;
+
+    //Solve for everything else
+    double E1 = (mass*mass - m2*m2 + m1*m1)/(2.0*mass);
+    double E2 = mass - E1;
+    double E3 = mass - E1;
+    double mp1 = sqrt(E1*E1 - m1*m1);
+
+    TVector3 p1(mp1*cos(phi1)*sin(th1),mp1*sin(phi1)*sin(th1),mp1*cos(th1));
+    TVector3 p2 = -1.0*p1;
+    TVector3 p3 = -2.0*p1;
+
+    triP output;
+    output.p1.SetPxPyPzE(p1.Px(),p1.Py(),p1.Pz(),E1);
+    output.p2.SetPxPyPzE(p2.Px(),p2.Py(),p2.Pz(),E2);
+    output.p3.SetPxPyPzE(p3.Px(),p3.Py(),p3.Pz(),E3);
+
+    if(abs(output.p1.Eta()) > 100 )
+    {
+        cout << "mass: " << mass << endl;
+        cout << "th1: " << th1 << endl;
+        cout << "E1: " << E1 << endl;
+        cout << "phi1: " << phi1 << endl;
+        cout << "mp1: " << mp1 << endl;
+    }
+
+    //Boost into lab frame
+    output.p1.Boost(boost);
+    output.p2.Boost(boost);
+
+    return output;
+
+
+};

include/mcTree.h

@@ -0,0 +1,172 @@
+//////////////////////////////////////////////////////////
+// This class has been automatically generated on
+// Mon Jul 31 14:53:13 2017 by ROOT version 6.06/01
+// from TTree mcTree/mcTree
+// found on file: test/testRun.root
+//////////////////////////////////////////////////////////
+
+#ifndef mcTree_h
+#define mcTree_h
+
+#include <TROOT.h>
+#include <TChain.h>
+#include <TFile.h>
+
+// Header file for the classes stored in the TTree if any.
+#include "vector"
+
+class mcTree {
+public :
+   TTree          *fChain;   //!pointer to the analyzed TTree or TChain
+   Int_t           fCurrent; //!current Tree number in a TChain
+
+// Fixed size dimensions of array or collections stored in the TTree if any.
+
+   // Declaration of leaf types
+   vector<double>  *daughtersPt;
+   vector<double>  *daughtersEta;
+   vector<double>  *daughtersPhi;
+   vector<double>  *daughtersE;
+   vector<double>  *daughtersID;
+   Double_t        weight;
+   Double_t        pz;
+   Double_t        x1;
+   Double_t        x2;
+   Int_t           iq1;
+   Int_t           iq2;
+   Double_t        Q2;
+   Double_t        momM;
+
+   // List of branches
+   TBranch        *b_daughtersPt;   //!
+   TBranch        *b_daughtersEta;   //!
+   TBranch        *b_daughtersPhi;   //!
+   TBranch        *b_daughtersE;   //!
+   TBranch        *b_daughtersID;   //!
+   TBranch        *b_weight;   //!
+   TBranch        *b_pz;   //!
+   TBranch        *b_x1;   //!
+   TBranch        *b_x2;   //!
+   TBranch        *b_iq1;   //!
+   TBranch        *b_iq2;   //!
+   TBranch        *b_Q2;   //!
+   TBranch        *b_momM;   //!
+
+   mcTree(TTree *tree=0);
+   virtual ~mcTree();
+   virtual Int_t    Cut(Long64_t entry);
+   virtual Int_t    GetEntry(Long64_t entry);
+   virtual Long64_t LoadTree(Long64_t entry);
+   virtual void     Init(TTree *tree);
+   virtual void     Loop();
+   virtual Bool_t   Notify();
+   virtual void     Show(Long64_t entry = -1);
+};
+
+#endif
+
+#ifdef mcTree_cxx
+mcTree::mcTree(TTree *tree) : fChain(0) 
+{
+// if parameter tree is not specified (or zero), connect the file
+// used to generate this class and read the Tree.
+   if (tree == 0) {
+      TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("test/testRun.root");
+      if (!f || !f->IsOpen()) {
+         f = new TFile("test/testRun.root");
+      }
+      f->GetObject("mcTree",tree);
+
+   }
+   Init(tree);
+}
+
+mcTree::~mcTree()
+{
+   if (!fChain) return;
+   delete fChain->GetCurrentFile();
+}
+
+Int_t mcTree::GetEntry(Long64_t entry)
+{
+// Read contents of entry.
+   if (!fChain) return 0;
+   return fChain->GetEntry(entry);
+}
+Long64_t mcTree::LoadTree(Long64_t entry)
+{
+// Set the environment to read one entry
+   if (!fChain) return -5;
+   Long64_t centry = fChain->LoadTree(entry);
+   if (centry < 0) return centry;
+   if (fChain->GetTreeNumber() != fCurrent) {
+      fCurrent = fChain->GetTreeNumber();
+      Notify();
+   }
+   return centry;
+}
+
+void mcTree::Init(TTree *tree)
+{
+   // The Init() function is called when the selector needs to initialize
+   // a new tree or chain. Typically here the branch addresses and branch
+   // pointers of the tree will be set.
+   // It is normally not necessary to make changes to the generated
+   // code, but the routine can be extended by the user if needed.
+   // Init() will be called many times when running on PROOF
+   // (once per file to be processed).
+
+   // Set object pointer
+   daughtersPt = 0;
+   daughtersEta = 0;
+   daughtersPhi = 0;
+   daughtersE = 0;
+   daughtersID = 0;
+   // Set branch addresses and branch pointers
+   if (!tree) return;
+   fChain = tree;
+   fCurrent = -1;
+   fChain->SetMakeClass(1);
+
+   fChain->SetBranchAddress("daughtersPt", &daughtersPt, &b_daughtersPt);
+   fChain->SetBranchAddress("daughtersEta", &daughtersEta, &b_daughtersEta);
+   fChain->SetBranchAddress("daughtersPhi", &daughtersPhi, &b_daughtersPhi);
+   fChain->SetBranchAddress("daughtersE", &daughtersE, &b_daughtersE);
+   fChain->SetBranchAddress("daughtersID", &daughtersID, &b_daughtersID);
+   fChain->SetBranchAddress("weight", &weight, &b_weight);
+   fChain->SetBranchAddress("pz", &pz, &b_pz);
+   fChain->SetBranchAddress("x1", &x1, &b_x1);
+   fChain->SetBranchAddress("x2", &x2, &b_x2);
+   fChain->SetBranchAddress("iq1", &iq1, &b_iq1);
+   fChain->SetBranchAddress("iq2", &iq2, &b_iq2);
+   fChain->SetBranchAddress("Q2", &Q2, &b_Q2);
+   fChain->SetBranchAddress("momM", &momM, &b_momM);
+   Notify();
+}
+
+Bool_t mcTree::Notify()
+{
+   // The Notify() function is called when a new file is opened. This
+   // can be either for a new TTree in a TChain or when when a new TTree
+   // is started when using PROOF. It is normally not necessary to make changes
+   // to the generated code, but the routine can be extended by the
+   // user if needed. The return value is currently not used.
+
+   return kTRUE;
+}
+
+void mcTree::Show(Long64_t entry)
+{
+// Print contents of entry.
+// If entry is not specified, print current entry
+   if (!fChain) return;
+   fChain->Show(entry);
+}
+Int_t mcTree::Cut(Long64_t entry)
+{
+// This function may be called from Loop.
+// returns  1 if entry is accepted.
+// returns -1 otherwise.
+   return 1;
+}
+#endif // #ifdef mcTree_cxx