From c63cd3b07d74b3076b9775a60b0f81a44d0c187b Mon Sep 17 00:00:00 2001 From: Cameron Bravo Date: Mon, 26 Mar 2018 15:59:16 -0700 Subject: [PATCH] Initial Commit of BaryoGEN, forked from sphaleronMC --- .gitignore | 5 + LICENSE | 439 +++++++++++++++++++++++++++++++ Makefile | 9 + README.md | 19 ++ include/BaryoGEN.h | 2 + include/LHEWriter.h | 22 ++ include/LinkDef.h | 3 + include/configBuilder.h | 22 ++ include/decay.h | 123 +++++++++ include/mcTree.h | 172 +++++++++++++ include/particles.h | 73 ++++++ include/sphalerons.h | 46 ++++ macros/mcTree.C | 43 ++++ src/BaryoGEN.cc | 557 ++++++++++++++++++++++++++++++++++++++++ src/LHEWriter.cc | 70 +++++ src/configBuilder.cc | 181 +++++++++++++ src/particles.cc | 76 ++++++ 17 files changed, 1862 insertions(+) create mode 100644 .gitignore create mode 100644 LICENSE create mode 100644 Makefile create mode 100644 README.md create mode 100644 include/BaryoGEN.h create mode 100644 include/LHEWriter.h create mode 100644 include/LinkDef.h create mode 100644 include/configBuilder.h create mode 100644 include/decay.h create mode 100644 include/mcTree.h create mode 100644 include/particles.h create mode 100644 include/sphalerons.h create mode 100644 macros/mcTree.C create mode 100644 src/BaryoGEN.cc create mode 100644 src/LHEWriter.cc create mode 100644 src/configBuilder.cc create mode 100644 src/particles.cc diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..a8beba3 --- /dev/null +++ b/.gitignore @@ -0,0 +1,5 @@ +MyDict.cxx +MyDict_rdict.pcm +genTree +*.root +*.lhe diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..3f3cb87 --- /dev/null +++ b/LICENSE @@ -0,0 +1,439 @@ +Attribution-NonCommercial-ShareAlike 4.0 International + +======================================================================= + +Creative Commons Corporation ("Creative Commons") is not a law firm and +does not provide legal services or legal advice. Distribution of +Creative Commons public licenses does not create a lawyer-client or +other relationship. Creative Commons makes its licenses and related +information available on an "as-is" basis. Creative Commons gives no +warranties regarding its licenses, any material licensed under their +terms and conditions, or any related information. Creative Commons +disclaims all liability for damages resulting from their use to the +fullest extent possible. + +Using Creative Commons Public Licenses + +Creative Commons public licenses provide a standard set of terms and +conditions that creators and other rights holders may use to share +original works of authorship and other material subject to copyright +and certain other rights specified in the public license below. The +following considerations are for informational purposes only, are not +exhaustive, and do not form part of our licenses. + + Considerations for licensors: Our public licenses are + intended for use by those authorized to give the public + permission to use material in ways otherwise restricted by + copyright and certain other rights. Our licenses are + irrevocable. Licensors should read and understand the terms + and conditions of the license they choose before applying it. + Licensors should also secure all rights necessary before + applying our licenses so that the public can reuse the + material as expected. Licensors should clearly mark any + material not subject to the license. This includes other CC- + licensed material, or material used under an exception or + limitation to copyright. More considerations for licensors: + wiki.creativecommons.org/Considerations_for_licensors + + Considerations for the public: By using one of our public + licenses, a licensor grants the public permission to use the + licensed material under specified terms and conditions. If + the licensor's permission is not necessary for any reason--for + example, because of any applicable exception or limitation to + copyright--then that use is not regulated by the license. Our + licenses grant only permissions under copyright and certain + other rights that a licensor has authority to grant. Use of + the licensed material may still be restricted for other + reasons, including because others have copyright or other + rights in the material. A licensor may make special requests, + such as asking that all changes be marked or described. + Although not required by our licenses, you are encouraged to + respect those requests where reasonable. More_considerations + for the public: + wiki.creativecommons.org/Considerations_for_licensees + +======================================================================= + +Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International +Public License + +By exercising the Licensed Rights (defined below), You accept and agree +to be bound by the terms and conditions of this Creative Commons +Attribution-NonCommercial-ShareAlike 4.0 International Public License +("Public License"). To the extent this Public License may be +interpreted as a contract, You are granted the Licensed Rights in +consideration of Your acceptance of these terms and conditions, and the +Licensor grants You such rights in consideration of benefits the +Licensor receives from making the Licensed Material available under +these terms and conditions. + + +Section 1 -- Definitions. + + a. Adapted Material means material subject to Copyright and Similar + Rights that is derived from or based upon the Licensed Material + and in which the Licensed Material is translated, altered, + arranged, transformed, or otherwise modified in a manner requiring + permission under the Copyright and Similar Rights held by the + Licensor. For purposes of this Public License, where the Licensed + Material is a musical work, performance, or sound recording, + Adapted Material is always produced where the Licensed Material is + synched in timed relation with a moving image. + + b. Adapter's License means the license You apply to Your Copyright + and Similar Rights in Your contributions to Adapted Material in + accordance with the terms and conditions of this Public License. + + c. BY-NC-SA Compatible License means a license listed at + creativecommons.org/compatiblelicenses, approved by Creative + Commons as essentially the equivalent of this Public License. + + d. Copyright and Similar Rights means copyright and/or similar rights + closely related to copyright including, without limitation, + performance, broadcast, sound recording, and Sui Generis Database + Rights, without regard to how the rights are labeled or + categorized. For purposes of this Public License, the rights + specified in Section 2(b)(1)-(2) are not Copyright and Similar + Rights. + + e. Effective Technological Measures means those measures that, in the + absence of proper authority, may not be circumvented under laws + fulfilling obligations under Article 11 of the WIPO Copyright + Treaty adopted on December 20, 1996, and/or similar international + agreements. + + f. Exceptions and Limitations means fair use, fair dealing, and/or + any other exception or limitation to Copyright and Similar Rights + that applies to Your use of the Licensed Material. + + g. License Elements means the license attributes listed in the name + of a Creative Commons Public License. The License Elements of this + Public License are Attribution, NonCommercial, and ShareAlike. + + h. Licensed Material means the artistic or literary work, database, + or other material to which the Licensor applied this Public + License. + + i. Licensed Rights means the rights granted to You subject to the + terms and conditions of this Public License, which are limited to + all Copyright and Similar Rights that apply to Your use of the + Licensed Material and that the Licensor has authority to license. + + j. Licensor means the individual(s) or entity(ies) granting rights + under this Public License. + + k. NonCommercial means not primarily intended for or directed towards + commercial advantage or monetary compensation. For purposes of + this Public License, the exchange of the Licensed Material for + other material subject to Copyright and Similar Rights by digital + file-sharing or similar means is NonCommercial provided there is + no payment of monetary compensation in connection with the + exchange. + + l. Share means to provide material to the public by any means or + process that requires permission under the Licensed Rights, such + as reproduction, public display, public performance, distribution, + dissemination, communication, or importation, and to make material + available to the public including in ways that members of the + public may access the material from a place and at a time + individually chosen by them. + + m. Sui Generis Database Rights means rights other than copyright + resulting from Directive 96/9/EC of the European Parliament and of + the Council of 11 March 1996 on the legal protection of databases, + as amended and/or succeeded, as well as other essentially + equivalent rights anywhere in the world. + + n. You means the individual or entity exercising the Licensed Rights + under this Public License. Your has a corresponding meaning. + + +Section 2 -- Scope. + + a. License grant. + + 1. Subject to the terms and conditions of this Public License, + the Licensor hereby grants You a worldwide, royalty-free, + non-sublicensable, non-exclusive, irrevocable license to + exercise the Licensed Rights in the Licensed Material to: + + a. reproduce and Share the Licensed Material, in whole or + in part, for NonCommercial purposes only; and + + b. produce, reproduce, and Share Adapted Material for + NonCommercial purposes only. + + 2. Exceptions and Limitations. For the avoidance of doubt, where + Exceptions and Limitations apply to Your use, this Public + License does not apply, and You do not need to comply with + its terms and conditions. + + 3. Term. The term of this Public License is specified in Section + 6(a). + + 4. Media and formats; technical modifications allowed. The + Licensor authorizes You to exercise the Licensed Rights in + all media and formats whether now known or hereafter created, + and to make technical modifications necessary to do so. The + Licensor waives and/or agrees not to assert any right or + authority to forbid You from making technical modifications + necessary to exercise the Licensed Rights, including + technical modifications necessary to circumvent Effective + Technological Measures. For purposes of this Public License, + simply making modifications authorized by this Section 2(a) + (4) never produces Adapted Material. + + 5. Downstream recipients. + + a. Offer from the Licensor -- Licensed Material. Every + recipient of the Licensed Material automatically + receives an offer from the Licensor to exercise the + Licensed Rights under the terms and conditions of this + Public License. + + b. Additional offer from the Licensor -- Adapted Material. + Every recipient of Adapted Material from You + automatically receives an offer from the Licensor to + exercise the Licensed Rights in the Adapted Material + under the conditions of the Adapter's License You apply. + + c. No downstream restrictions. You may not offer or impose + any additional or different terms or conditions on, or + apply any Effective Technological Measures to, the + Licensed Material if doing so restricts exercise of the + Licensed Rights by any recipient of the Licensed + Material. + + 6. No endorsement. Nothing in this Public License constitutes or + may be construed as permission to assert or imply that You + are, or that Your use of the Licensed Material is, connected + with, or sponsored, endorsed, or granted official status by, + the Licensor or others designated to receive attribution as + provided in Section 3(a)(1)(A)(i). + + b. Other rights. + + 1. Moral rights, such as the right of integrity, are not + licensed under this Public License, nor are publicity, + privacy, and/or other similar personality rights; however, to + the extent possible, the Licensor waives and/or agrees not to + assert any such rights held by the Licensor to the limited + extent necessary to allow You to exercise the Licensed + Rights, but not otherwise. + + 2. Patent and trademark rights are not licensed under this + Public License. + + 3. To the extent possible, the Licensor waives any right to + collect royalties from You for the exercise of the Licensed + Rights, whether directly or through a collecting society + under any voluntary or waivable statutory or compulsory + licensing scheme. In all other cases the Licensor expressly + reserves any right to collect such royalties, including when + the Licensed Material is used other than for NonCommercial + purposes. + + +Section 3 -- License Conditions. + +Your exercise of the Licensed Rights is expressly made subject to the +following conditions. + + a. Attribution. + + 1. If You Share the Licensed Material (including in modified + form), You must: + + a. retain the following if it is supplied by the Licensor + with the Licensed Material: + + i. identification of the creator(s) of the Licensed + Material and any others designated to receive + attribution, in any reasonable manner requested by + the Licensor (including by pseudonym if + designated); + + ii. a copyright notice; + + iii. a notice that refers to this Public License; + + iv. a notice that refers to the disclaimer of + warranties; + + v. a URI or hyperlink to the Licensed Material to the + extent reasonably practicable; + + b. indicate if You modified the Licensed Material and + retain an indication of any previous modifications; and + + c. indicate the Licensed Material is licensed under this + Public License, and include the text of, or the URI or + hyperlink to, this Public License. + + 2. You may satisfy the conditions in Section 3(a)(1) in any + reasonable manner based on the medium, means, and context in + which You Share the Licensed Material. For example, it may be + reasonable to satisfy the conditions by providing a URI or + hyperlink to a resource that includes the required + information. + 3. If requested by the Licensor, You must remove any of the + information required by Section 3(a)(1)(A) to the extent + reasonably practicable. + + b. ShareAlike. + + In addition to the conditions in Section 3(a), if You Share + Adapted Material You produce, the following conditions also apply. + + 1. The Adapter's License You apply must be a Creative Commons + license with the same License Elements, this version or + later, or a BY-NC-SA Compatible License. + + 2. You must include the text of, or the URI or hyperlink to, the + Adapter's License You apply. You may satisfy this condition + in any reasonable manner based on the medium, means, and + context in which You Share Adapted Material. + + 3. You may not offer or impose any additional or different terms + or conditions on, or apply any Effective Technological + Measures to, Adapted Material that restrict exercise of the + rights granted under the Adapter's License You apply. + + +Section 4 -- Sui Generis Database Rights. + +Where the Licensed Rights include Sui Generis Database Rights that +apply to Your use of the Licensed Material: + + a. for the avoidance of doubt, Section 2(a)(1) grants You the right + to extract, reuse, reproduce, and Share all or a substantial + portion of the contents of the database for NonCommercial purposes + only; + + b. if You include all or a substantial portion of the database + contents in a database in which You have Sui Generis Database + Rights, then the database in which You have Sui Generis Database + Rights (but not its individual contents) is Adapted Material, + including for purposes of Section 3(b); and + + c. You must comply with the conditions in Section 3(a) if You Share + all or a substantial portion of the contents of the database. + +For the avoidance of doubt, this Section 4 supplements and does not +replace Your obligations under this Public License where the Licensed +Rights include other Copyright and Similar Rights. + + +Section 5 -- Disclaimer of Warranties and Limitation of Liability. + + a. UNLESS OTHERWISE SEPARATELY UNDERTAKEN BY THE LICENSOR, TO THE + EXTENT POSSIBLE, THE LICENSOR OFFERS THE LICENSED MATERIAL AS-IS + AND AS-AVAILABLE, AND MAKES NO REPRESENTATIONS OR WARRANTIES OF + ANY KIND CONCERNING THE LICENSED MATERIAL, WHETHER EXPRESS, + IMPLIED, STATUTORY, OR OTHER. THIS INCLUDES, WITHOUT LIMITATION, + WARRANTIES OF TITLE, MERCHANTABILITY, FITNESS FOR A PARTICULAR + PURPOSE, NON-INFRINGEMENT, ABSENCE OF LATENT OR OTHER DEFECTS, + ACCURACY, OR THE PRESENCE OR ABSENCE OF ERRORS, WHETHER OR NOT + KNOWN OR DISCOVERABLE. WHERE DISCLAIMERS OF WARRANTIES ARE NOT + ALLOWED IN FULL OR IN PART, THIS DISCLAIMER MAY NOT APPLY TO YOU. + + b. TO THE EXTENT POSSIBLE, IN NO EVENT WILL THE LICENSOR BE LIABLE + TO YOU ON ANY LEGAL THEORY (INCLUDING, WITHOUT LIMITATION, + NEGLIGENCE) OR OTHERWISE FOR ANY DIRECT, SPECIAL, INDIRECT, + INCIDENTAL, CONSEQUENTIAL, PUNITIVE, EXEMPLARY, OR OTHER LOSSES, + COSTS, EXPENSES, OR DAMAGES ARISING OUT OF THIS PUBLIC LICENSE OR + USE OF THE LICENSED MATERIAL, EVEN IF THE LICENSOR HAS BEEN + ADVISED OF THE POSSIBILITY OF SUCH LOSSES, COSTS, EXPENSES, OR + DAMAGES. WHERE A LIMITATION OF LIABILITY IS NOT ALLOWED IN FULL OR + IN PART, THIS LIMITATION MAY NOT APPLY TO YOU. + + c. The disclaimer of warranties and limitation of liability provided + above shall be interpreted in a manner that, to the extent + possible, most closely approximates an absolute disclaimer and + waiver of all liability. + + +Section 6 -- Term and Termination. + + a. This Public License applies for the term of the Copyright and + Similar Rights licensed here. However, if You fail to comply with + this Public License, then Your rights under this Public License + terminate automatically. + + b. Where Your right to use the Licensed Material has terminated under + Section 6(a), it reinstates: + + 1. automatically as of the date the violation is cured, provided + it is cured within 30 days of Your discovery of the + violation; or + + 2. upon express reinstatement by the Licensor. + + For the avoidance of doubt, this Section 6(b) does not affect any + right the Licensor may have to seek remedies for Your violations + of this Public License. + + c. For the avoidance of doubt, the Licensor may also offer the + Licensed Material under separate terms or conditions or stop + distributing the Licensed Material at any time; however, doing so + will not terminate this Public License. + + d. Sections 1, 5, 6, 7, and 8 survive termination of this Public + License. + + +Section 7 -- Other Terms and Conditions. + + a. The Licensor shall not be bound by any additional or different + terms or conditions communicated by You unless expressly agreed. + + b. Any arrangements, understandings, or agreements regarding the + Licensed Material not stated herein are separate from and + independent of the terms and conditions of this Public License. + + +Section 8 -- Interpretation. + + a. For the avoidance of doubt, this Public License does not, and + shall not be interpreted to, reduce, limit, restrict, or impose + conditions on any use of the Licensed Material that could lawfully + be made without permission under this Public License. + + b. To the extent possible, if any provision of this Public License is + deemed unenforceable, it shall be automatically reformed to the + minimum extent necessary to make it enforceable. If the provision + cannot be reformed, it shall be severed from this Public License + without affecting the enforceability of the remaining terms and + conditions. + + c. No term or condition of this Public License will be waived and no + failure to comply consented to unless expressly agreed to by the + Licensor. + + d. Nothing in this Public License constitutes or may be interpreted + as a limitation upon, or waiver of, any privileges and immunities + that apply to the Licensor or You, including from the legal + processes of any jurisdiction or authority. + +======================================================================= + +Creative Commons is not a party to its public +licenses. Notwithstanding, Creative Commons may elect to apply one of +its public licenses to material it publishes and in those instances +will be considered the “Licensor.” The text of the Creative Commons +public licenses is dedicated to the public domain under the CC0 Public +Domain Dedication. Except for the limited purpose of indicating that +material is shared under a Creative Commons public license or as +otherwise permitted by the Creative Commons policies published at +creativecommons.org/policies, Creative Commons does not authorize the +use of the trademark "Creative Commons" or any other trademark or logo +of Creative Commons without its prior written consent including, +without limitation, in connection with any unauthorized modifications +to any of its public licenses or any other arrangements, +understandings, or agreements concerning use of licensed material. For +the avoidance of doubt, this paragraph does not form part of the +public licenses. + +Creative Commons may be contacted at creativecommons.org. + + diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..4740b79 --- /dev/null +++ b/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` diff --git a/README.md b/README.md new file mode 100644 index 0000000..e96857c --- /dev/null +++ b/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. diff --git a/include/BaryoGEN.h b/include/BaryoGEN.h new file mode 100644 index 0000000..2e74a3c --- /dev/null +++ b/include/BaryoGEN.h @@ -0,0 +1,2 @@ +#pragma link C++ class std::vector < TLorentzVector >+; +#pragma link C++ class std::vector < std::vector >+; diff --git a/include/LHEWriter.h b/include/LHEWriter.h new file mode 100644 index 0000000..0ca49ea --- /dev/null +++ b/include/LHEWriter.h @@ -0,0 +1,22 @@ +#include +#include +#include + +#include "particles.h" +#include "TLorentzVector.h" + +using namespace std; + +class LHEWriter +{ + private: + ofstream oF; + + public: + LHEWriter(string fName); + ~LHEWriter(); + + int writeEvent(vector outParts, double Q); + int close(); + +}; diff --git a/include/LinkDef.h b/include/LinkDef.h new file mode 100644 index 0000000..08578ca --- /dev/null +++ b/include/LinkDef.h @@ -0,0 +1,3 @@ +#pragma link C++ class std::vector +; +#pragma link C++ class std::vector < std::vector >+; + diff --git a/include/configBuilder.h b/include/configBuilder.h new file mode 100644 index 0000000..1f5e400 --- /dev/null +++ b/include/configBuilder.h @@ -0,0 +1,22 @@ +#include "particles.h" + +#include "TRandom3.h" + +#include + +using namespace std; + +class configBuilder +{ + private: + vector parts; + TRandom3 rand; + + public: + configBuilder(); + ~configBuilder(); + + vector build(int iQ1, int iQ2, float pNCS, bool bCancel); + +}; + diff --git a/include/decay.h b/include/decay.h new file mode 100644 index 0000000..4ea30e0 --- /dev/null +++ b/include/decay.h @@ -0,0 +1,123 @@ +#include +#include + +#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; + + +}; diff --git a/include/mcTree.h b/include/mcTree.h new file mode 100644 index 0000000..80e4864 --- /dev/null +++ b/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 +#include +#include + +// 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 *daughtersPt; + vector *daughtersEta; + vector *daughtersPhi; + vector *daughtersE; + vector *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 diff --git a/include/particles.h b/include/particles.h new file mode 100644 index 0000000..e2cfbd5 --- /dev/null +++ b/include/particles.h @@ -0,0 +1,73 @@ +#ifndef _PARTICLEH_ +#define _PARTICLEH_ + +#include +#include "TLorentzVector.h" + +#define ELE_MASS 0.000511 +#define ENU_MASS 0.0 +#define MU_MASS 0.10566 +#define MNU_MASS 0.0 +#define TAU_MASS 1.7768 +#define TNU_MASS 0.0 +#define TQ_MASS 173.34 +#define BQ_MASS 4.18 +#define CQ_MASS 1.275 +#define SQ_MASS 0.095 +#define UQ_MASS 0.0023 +#define DQ_MASS 0.0048 + +#define ELE_PID 11 +#define ENU_PID 12 +#define MU_PID 13 +#define MNU_PID 14 +#define TAU_PID 15 +#define TNU_PID 16 +#define TQ_PID 6 +#define BQ_PID 5 +#define CQ_PID 4 +#define SQ_PID 3 +#define UQ_PID 2 +#define DQ_PID 1 + +#define ELE_Q3 -3 +#define ENU_Q3 0 +#define MU_Q3 -3 +#define MNU_Q3 0 +#define TAU_Q3 -3 +#define TNU_Q3 0 +#define TQ_Q3 2 +#define BQ_Q3 -1 +#define CQ_Q3 2 +#define SQ_Q3 -1 +#define UQ_Q3 2 +#define DQ_Q3 -1 + +struct particle +{ + int pid; + float mass; + int color; + int q3; + TLorentzVector p4v; + int m1; + int m2; + int ic; +}; + +using namespace std; + +class particleBase +{ + private: + vector partLUT; + + public: + particleBase(); + ~particleBase(); + + particle getParticle(int pid); + +}; + +#endif diff --git a/include/sphalerons.h b/include/sphalerons.h new file mode 100644 index 0000000..03e2f6f --- /dev/null +++ b/include/sphalerons.h @@ -0,0 +1,46 @@ +#include "particles.h" + +double sphal[2][12]; + +int fillSphal() +{ + + sphal[0][0] = ELE_MASS; + sphal[1][0] = ELE_PID; + + sphal[0][1] = MU_MASS; + sphal[1][1] = MU_PID; + + sphal[0][2] = TAU_MASS; + sphal[1][2] = TAU_PID; + + sphal[0][3] = TQ_MASS; + sphal[1][3] = TQ_PID; + + sphal[0][4] = TQ_MASS; + sphal[1][4] = TQ_PID; + + sphal[0][5] = BQ_MASS; + sphal[1][5] = BQ_PID; + + sphal[0][6] = CQ_MASS; + sphal[1][6] = CQ_PID; + + sphal[0][7] = CQ_MASS; + sphal[1][7] = CQ_PID; + + sphal[0][8] = SQ_MASS; + sphal[1][8] = SQ_PID; + + sphal[0][9] = UQ_MASS; + sphal[1][9] = UQ_PID; + + sphal[0][10] = UQ_MASS; + sphal[1][10] = UQ_PID; + + sphal[0][11] = DQ_MASS; + sphal[1][11] = DQ_PID; + + return 0; +} + diff --git a/macros/mcTree.C b/macros/mcTree.C new file mode 100644 index 0000000..da9c57e --- /dev/null +++ b/macros/mcTree.C @@ -0,0 +1,43 @@ +#define mcTree_cxx +#include "mcTree.h" +#include +#include +#include + +void mcTree::Loop() +{ +// In a ROOT session, you can do: +// root> .L mcTree.C +// root> mcTree t +// root> t.GetEntry(12); // Fill t data members with entry number 12 +// root> t.Show(); // Show values of entry 12 +// root> t.Show(16); // Read and show values of entry 16 +// root> t.Loop(); // Loop on all entries +// + +// This is the loop skeleton where: +// jentry is the global entry number in the chain +// ientry is the entry number in the current Tree +// Note that the argument to GetEntry must be: +// jentry for TChain::GetEntry +// ientry for TTree::GetEntry and TBranch::GetEntry +// +// To read only selected branches, Insert statements like: +// METHOD1: +// fChain->SetBranchStatus("*",0); // disable all branches +// fChain->SetBranchStatus("branchname",1); // activate branchname +// METHOD2: replace line +// fChain->GetEntry(jentry); //read all branches +//by b_branchname->GetEntry(ientry); //read only this branch + if (fChain == 0) return; + + Long64_t nentries = fChain->GetEntriesFast(); + + Long64_t nbytes = 0, nb = 0; + for (Long64_t jentry=0; jentryGetEntry(jentry); nbytes += nb; + // if (Cut(ientry) < 0) continue; + } +} diff --git a/src/BaryoGEN.cc b/src/BaryoGEN.cc new file mode 100644 index 0000000..68e5a44 --- /dev/null +++ b/src/BaryoGEN.cc @@ -0,0 +1,557 @@ +#include +#include +#include +#include +#include + +#include "LHAPDF/LHAPDF.h" + +#include "TLorentzVector.h" +#include "TGenPhaseSpace.h" +#include "TRandom3.h" +#include "TTree.h" +#include "TFile.h" +#include "TF1.h" +#include "TH1.h" +#include "TH2.h" +#include "TMath.h" +#include "TString.h" + +#include "../include/configBuilder.h" +#include "../include/BaryoGEN.h" +#include "../include/LHEWriter.h" + +#define ARGS 7 +#define MPW 5.0e-7 + +using namespace LHAPDF; +using namespace std; + +int main(int argc, char* argv[]) +{ + + //read arguements and do conversions and checks on them + if(argc != ARGS+1) + { + cout << "genTree sqrt(s) threshold maxweight Nevents pNCS bCancel Filename" << endl; + return -99; + } + + float SQRTS = atof(argv[1]); + float thr = atof(argv[2]); + float MCW = atof(argv[3]); + if(MCW < 0.0) + { + cout << "maxweight must be greater than zero" << endl; + return -99; + } + int Nevt = atoi(argv[4]); + float pNCS = atof(argv[5]); + if(pNCS > 1 || pNCS < 0) + { + cout << "pNCS must be in [0,1]" << endl; + return -99; + } + int bCancel = atoi(argv[6]); + if(bCancel < 0 || bCancel > 1) + { + cout << "bCancel must be 0 or 1" << endl; + return -99; + } + bool bCan = bCancel; + string ofName = string(argv[7]); + + //Init structures needed during generatation + double maxwt = 0; + double maxMCtot = 0.0; + + int NF = 0; + TRandom3 rand; + rand.SetSeed(0); + particleBase *partBase = new particleBase(); + configBuilder confBuild; + vector daughtersPt; + vector daughtersEta; + vector daughtersPhi; + vector daughtersE; + vector daughtersID; + TLorentzVector mom(0.0,0.0,0.0,thr); + TLorentzVector u1(0.0,0.0,0.0,0.0); + TLorentzVector u2(0.0,0.0,0.0,0.0); + double masses[15]; + + double weight = 0.0; + + TGenPhaseSpace gen; + double x1 = 0.0; + double x2 = 0.0; + int iq1 = 0; + int iq2 = 0; + double Q2 = 0.0; + double momM = 0.0; + double pz = 0.0; + + //Init output files + TFile *myF = new TFile((ofName+".root").c_str(),"RECREATE","Holds daughters from sphaleron decay"); + LHEWriter lheF(ofName); + + double minx = thr*thr/(SQRTS*SQRTS); + + //Init parton distribution functions + const PDF* LHApdf = mkPDF("NNPDF30_lo_as_0118",0); + cout << LHApdf->xfxQ2(2, 0.5, SQRTS*SQRTS) << endl; + + //Initialize histograms for debugging + TH1D *x1_h = new TH1D("x1_h","x1 inclusive",1000,0.0,1.0); + TH1D *mcTot_h = new TH1D("mcTot_h","Monte Carlo Probabilities",100,0.0,MCW); + TH1D *sumInterQ3_h = new TH1D("sumInterQ3_h","Intermediate particle charges",21,-10.5,10.5); + TH1D *sphM_h = new TH1D("sphM_h","Sphaleron Transition Energy;Invariant Mass [TeV];Events / 100 GeV",int(SQRTS/100.0),0.0,SQRTS/1000.0); + TH1D *st_h = new TH1D("st_h","S_{T};S_{T} [TeV];Events / 100 GeV", int(SQRTS/100.0), 0.0, SQRTS/1000.0); + TH1D *sphPz_h = new TH1D("sphPz_h","Sphaleron p_{z};p_{z} [GeV];Events / 100 GeV",80,-4000.0,4000.0); + TH1D *outID_h = new TH1D("outID_h","Outgoing PDG IDs;PDG ID;Entries",33,-16.5,16.5); + TH1D *p1x_h = new TH1D("p1x_h","Parton 1 Momentum Fraction;x_{1};Events / 0.01",60,0.4,1.0); + TH1D *p1id_h = new TH1D("p1id_h","Parton 1 PDG ID;PDG ID;Events",13,-6.5,6.5); + + TH2D *inQid_h = new TH2D("inQid_h","Colliding Parton Species;Parton 2 PDG ID;Parton 1 PDG ID",11,-5.5,5.5,11,-5.5,5.5); + TH2D *frac2D_h = new TH2D("frac2D_h","Parton Momentum Fractions;x_{2};x_{1}",60,0.4,1.0,60,0.4,1.0); + + vector pt_hv; + vector eta_hv; + vector phi_hv; + for(int iq = -6; iq < 7; iq++) + { + string nameBuf; + string titBuf; + titBuf = Form("Quark %i p_{T};p_{T} [GeV];Entries / 20 GeV",iq); + if(iq < 0) nameBuf = Form("pt_iqm%i_h",int(fabs(iq))); + else nameBuf = Form("pt_iqp%i_h",int(fabs(iq))); + TH1D hBufpt(nameBuf.c_str(),titBuf.c_str(),250,0.0,5000.0); + pt_hv.push_back(hBufpt); + + titBuf = Form("Quark %i #eta;#eta;Entries / 0.1",iq); + if(iq < 0) nameBuf = Form("eta_iqm%i_h",int(fabs(iq))); + else nameBuf = Form("eta_iqp%i_h",int(fabs(iq))); + TH1D hBufeta(nameBuf.c_str(),titBuf.c_str(),100,-5,5); + eta_hv.push_back(hBufeta); + + titBuf = Form("Quark %i #phi;#phi;Entries / 0.1 #pi",iq); + if(iq < 0) nameBuf = Form("phi_iqm%i_h",int(fabs(iq))); + else nameBuf = Form("phi_iqp%i_h",int(fabs(iq))); + TH1D hBufphi(nameBuf.c_str(),titBuf.c_str(),20,-1.0*TMath::Pi(),TMath::Pi()); + phi_hv.push_back(hBufphi); + } + + //Initialize TTree to save event info + TTree *myT = new TTree("mcTree","mcTree"); + myT->Branch("daughtersPt",&daughtersPt); + myT->Branch("daughtersEta",&daughtersEta); + myT->Branch("daughtersPhi",&daughtersPhi); + myT->Branch("daughtersE",&daughtersE); + myT->Branch("daughtersID",&daughtersID); + myT->Branch("weight",&weight); + myT->Branch("pz",&pz); + myT->Branch("x1",&x1); + myT->Branch("x2",&x2); + myT->Branch("iq1",&iq1); + myT->Branch("iq2",&iq2); + myT->Branch("Q2",&Q2); + myT->Branch("momM",&momM); + + int pdfN = 0; + + while(NF < Nevt)//Main event generation loop + { + //Make sure everything is reset from the last event + daughtersPt.clear(); + daughtersEta.clear(); + daughtersPhi.clear(); + daughtersE.clear(); + daughtersID.clear(); + + Q2 = 0.0; + double mcP = 1.1; + bool mcPass = false; + //Collide protons + while(!mcPass) + { + //Choose x1 and x2 in proper range + Q2 = 0.0; + while(Q2 < thr*thr) + { + x1 = (1.0-minx)*rand.Uniform()+minx; + x2 = (1.0-minx)*rand.Uniform()+minx; + Q2 = x1*x2*SQRTS*SQRTS; + } + mcP = MCW*rand.Uniform(); + double mcTot = 0.0; + for(int i1 = -5; i1 < 6; i1++) + { + if(i1 == 0) continue; + iq1 = i1; + double x1p = LHApdf->xfxQ2(iq1, x1, Q2)/x1; + for(int i2 = -5; i2 < 6; i2++) + { + if(i2 == 0) continue; + iq2 = i2; + //double x2p = pdf->parton(iq2,x2,SQRTS); + double x2p = LHApdf->xfxQ2(iq2, x2, Q2)/x2; + mcTot += x1p*x2p; + if(mcTot > mcP) {mcPass = true; break;} + } + if(mcPass) break; + } + pdfN++; + x1_h->Fill(x1); + mcTot_h->Fill(mcTot); + if(maxMCtot < mcTot) {maxMCtot = mcTot; cout << "Max MC Total: " << maxMCtot << endl;} + } + + //Build incoming particles, sphaleron, and prepare to decay + vector inParts; + particle partBuf = partBase->getParticle(iq1); + if(iq1 != partBuf.pid) cout << "iq1 = " << iq1 << " != partBuf.pid = " << partBuf.pid << endl; + partBuf.color = 501; + partBuf.p4v.SetXYZM(0.0,0.0,x1*6500,partBuf.mass); + inParts.push_back(partBuf); //Push first incoming quark + + partBuf = partBase->getParticle(iq2); + partBuf.color = 502; + partBuf.p4v.SetXYZM(0.0,0.0,x2*(-6500),partBuf.mass); + inParts.push_back(partBuf); //Push second incoming quark + + int colNow = 503;//This is used to keep track of the color line numbers already used + + //Build the mother particle from incoming partons + particle mother; + mother.p4v = inParts[0].p4v + inParts[1].p4v; + mother.mass = mother.p4v.M(); + + //Generate vector of outgoing fermionic configuration + vector confBuf = confBuild.build(iq1,iq2,pNCS,bCan); + int Nline = 0; + for(int i = 0; i < confBuf.size(); i++) + { + masses[i] = confBuf[i].mass; + } + + u1.SetXYZM(0.0,0.0,x1*6500,UQ_MASS); + u2.SetXYZM(0.0,0.0,x2*(-6500),UQ_MASS); + mom = u1 + u2; + momM = mom.M(); + pz = mom.Pz(); + + //"Decay" mother to the generated configuration + weight = 1.0; + while(weight > MPW*rand.Uniform()) + { + gen.SetDecay(mom, confBuf.size(), masses); + + weight = gen.Generate(); + if(weight > MPW) cout << "The code needs to recompiled with a higher MPW" << endl; + } + + //Extract Decay 4-vectors and assign colors to non-spectator quarks + for(int ii = 0; ii < confBuf.size(); ii++) + { + TLorentzVector prod = *gen.GetDecay(ii); + confBuf[ii].p4v = prod; + confBuf[ii].m1 = 0; + confBuf[ii].m2 = 0; + if(confBuf[ii].color != 501 && confBuf[ii].color != 502 && fabs(confBuf[ii].pid) < 7) confBuf[ii].color = colNow++; + if(fabs(confBuf[ii].pid) > 7) confBuf[ii].color = 0; + int kinI = confBuf[ii].pid + 6; + if(confBuf[ii].color != 501 && confBuf[ii].color != 502 && fabs(confBuf[ii].pid) < 7) + { + pt_hv[kinI].Fill(confBuf[ii].p4v.Pt()); + eta_hv[kinI].Fill(confBuf[ii].p4v.Eta()); + phi_hv[kinI].Fill(confBuf[ii].p4v.Phi()); + } + outID_h->Fill(confBuf[ii].pid); + } + + vector g1q; + vector g2q; + vector g3q; + vector leps; + vector specParts; + + vector qIg1; + vector qIg2; + vector qIg3; + int Nlep = 0; + vector lepI; + //Seperate leptons and quarks + for(int i = 0; i < confBuf.size(); i++) + { + if(confBuf[i].color == 501 || confBuf[i].color == 502) + { + confBuf[i].m1 = 1; + confBuf[i].m2 = 2; + specParts.push_back(confBuf[i]); + continue; + } + if(fabs(confBuf[i].pid) > 7) //lepton + { + leps.push_back(confBuf[i]); + lepI.push_back(i); + Nlep++; + } + else //quark + { + if(fabs(confBuf[i].pid) > 4) { qIg3.push_back(i); g3q.push_back(confBuf[i]); } + else if(fabs(confBuf[i].pid) > 2) { qIg2.push_back(i); g2q.push_back(confBuf[i]); } + else { qIg1.push_back(i); g1q.push_back(confBuf[i]); } + } + } + int interCol = 0; + int icolS = 0; + if(specParts.size() == 1) + { + if(specParts[0].color == 501) + { + interCol = 502; + if(iq2 > 0) icolS = 1; + else icolS = -1; + } + if(specParts[0].color == 502) + { + interCol = 501; + if(iq1 > 0) icolS = 1; + else icolS = -1; + } + } + + vector fileParts; + fileParts.push_back(inParts[0]); + fileParts.push_back(inParts[1]); + + //Calculate color flow case + int specN = specParts.size(); + vector anParts; + for(int i = 0; i < inParts.size(); i++) + { + bool mF = false; + for(int ii = 0; ii < specParts.size(); ii++) if(inParts[i].color == specParts[ii].color) mF = true; + if(!mF) anParts.push_back(inParts[i]); + } + bool dcg3 = (specParts.size() == 0 && g3q.size() == 2); + bool dcg2 = (specParts.size() == 0 && g2q.size() == 2); + bool dcg1 = (specParts.size() == 0 && g1q.size() == 2); + + //Push third Gen quarks and fake susy mediator onto output stack + TLorentzVector interP(0.0,0.0,0.0,0.0); + int interQ3 = 0; + for(int ii = 0; ii < g3q.size(); ii++) + { + if(g3q[ii].color == 501 || g3q[ii].color == 502) continue; + interP = interP + g3q[ii].p4v; + interQ3 += g3q[ii].q3; + } + sumInterQ3_h->Fill(interQ3); + particle interBuf; + interBuf.mass = interP.M(); + interBuf.color = 0; + if(g3q.size() == 2) + { + interBuf.color = interCol; + interBuf.ic = icolS; + } + interBuf.p4v = interP; + interBuf.q3 = interQ3; + interBuf.m1 = 1; + interBuf.m2 = 2; + if(interQ3 == 0) interBuf.pid = 1000022; + if(fabs(interQ3) == 3) interBuf.pid = interQ3*1006213/fabs(interQ3); + if(fabs(interQ3) == 6) interBuf.pid = interQ3*1006223/fabs(interQ3); + if(g3q.size() == 2) for(int i = 0; i < anParts.size(); i++) if(fabs(anParts[i].pid) == 5 || fabs(anParts[i].pid) == 6) + { + interBuf.pid = anParts[i].pid; + interBuf.color = anParts[i].color; + } + bool useInter = (g3q.size() > 1); + if(useInter) fileParts.push_back(interBuf); + int interI = fileParts.size(); + for(int ii = 0; ii < g3q.size(); ii++) + { + g3q[ii].m1 = interI; + g3q[ii].m2 = interI; + if(!useInter) + { + g3q[ii].m1 = 1; + g3q[ii].m2 = 2; + g3q[ii].ic = icolS; + } + fileParts.push_back(g3q[ii]); + daughtersPt.push_back(g3q[ii].p4v.Pt()); + daughtersEta.push_back(g3q[ii].p4v.Eta()); + daughtersPhi.push_back(g3q[ii].p4v.Phi()); + daughtersE.push_back(g3q[ii].p4v.E()); + daughtersID.push_back(g3q[ii].pid); + } + + //Push second Gen quarks onto output stack + interP.SetPxPyPzE(0.0,0.0,0.0,0.0); + interQ3 = 0; + for(int ii = 0; ii < g2q.size(); ii++) + { + if(g2q[ii].color == 501 || g2q[ii].color == 502) continue; + interP = interP + g2q[ii].p4v; + interQ3 += g2q[ii].q3; + } + sumInterQ3_h->Fill(interQ3); + interBuf.mass = interP.M(); + interBuf.color = 0; + if(g2q.size() == 2) + { + interBuf.color = interCol; + interBuf.ic = icolS; + } + interBuf.p4v = interP; + interBuf.q3 = interQ3; + if(interQ3 == 0) interBuf.pid = 1000022; + if(fabs(interQ3) == 3) interBuf.pid = interQ3*1006213/fabs(interQ3); + if(fabs(interQ3) == 6) interBuf.pid = interQ3*1006223/fabs(interQ3); + if(dcg2 && dcg1) for(int i = 0; i < anParts.size(); i++) if(fabs(anParts[i].pid) == 3 || fabs(anParts[i].pid) == 4) + { + interBuf.pid = anParts[i].pid; + interBuf.color = anParts[i].color; + } + useInter = (g2q.size() == 3 || (dcg2 && dcg1)); + if(useInter) fileParts.push_back(interBuf); + interI = fileParts.size(); + for(int ii = 0; ii < g2q.size(); ii++) + { + g2q[ii].m1 = interI; + g2q[ii].m2 = interI; + if(!useInter) + { + g2q[ii].m1 = 1; + g2q[ii].m2 = 2; + g2q[ii].ic = icolS; + } + fileParts.push_back(g2q[ii]); + daughtersPt.push_back(g2q[ii].p4v.Pt()); + daughtersEta.push_back(g2q[ii].p4v.Eta()); + daughtersPhi.push_back(g2q[ii].p4v.Phi()); + daughtersE.push_back(g2q[ii].p4v.E()); + daughtersID.push_back(g2q[ii].pid); + } + + //Push first Gen quarks onto output stack + interP.SetPxPyPzE(0.0,0.0,0.0,0.0); + interQ3 = 0; + for(int ii = 0; ii < g1q.size(); ii++) + { + if(g1q[ii].color == 501 || g1q[ii].color == 502) continue; + interP = interP + g1q[ii].p4v; + interQ3 += g1q[ii].q3; + } + sumInterQ3_h->Fill(interQ3); + interBuf.mass = interP.M(); + interBuf.color = 0; + if(g1q.size() == 2) + { + interBuf.color = interCol; + interBuf.ic = icolS; + } + interBuf.p4v = interP; + interBuf.q3 = interQ3; + if(interQ3 == 0) interBuf.pid = 1000022; + if(fabs(interQ3) == 3) interBuf.pid = interQ3*1006213/fabs(interQ3); + if(fabs(interQ3) == 6) interBuf.pid = interQ3*1006223/fabs(interQ3); + useInter = (g1q.size() == 3); + if(useInter) fileParts.push_back(interBuf); + interI = fileParts.size(); + for(int ii = 0; ii < g1q.size(); ii++) + { + g1q[ii].m1 = interI; + g1q[ii].m2 = interI; + if(!useInter) + { + g1q[ii].m1 = 1; + g1q[ii].m2 = 2; + } + fileParts.push_back(g1q[ii]); + daughtersPt.push_back(g1q[ii].p4v.Pt()); + daughtersEta.push_back(g1q[ii].p4v.Eta()); + daughtersPhi.push_back(g1q[ii].p4v.Phi()); + daughtersE.push_back(g1q[ii].p4v.E()); + daughtersID.push_back(g1q[ii].pid); + } + + //Push leptons onto output stack + for(int i = 0; i < leps.size(); i++) + { + leps[i].m1 = 1; + leps[i].m2 = 2; + fileParts.push_back(leps[i]); + daughtersPt.push_back(leps[i].p4v.Pt()); + daughtersEta.push_back(leps[i].p4v.Eta()); + daughtersPhi.push_back(leps[i].p4v.Phi()); + daughtersE.push_back(leps[i].p4v.E()); + daughtersID.push_back(leps[i].pid); + } + + //Push spectator Quarks if there are any + for(int i = 0; i < specParts.size(); i++) + { + fileParts.push_back(specParts[i]); + daughtersPt.push_back(specParts[i].p4v.Pt()); + daughtersEta.push_back(specParts[i].p4v.Eta()); + daughtersPhi.push_back(specParts[i].p4v.Phi()); + daughtersE.push_back(specParts[i].p4v.E()); + daughtersID.push_back(specParts[i].pid); + } + + if(weight > maxwt) + { + maxwt = weight; + cout << "MaxWt: " << maxwt << endl; + } + + inQid_h->Fill(iq2,iq1); + frac2D_h->Fill(x2,x1); + sphM_h->Fill(momM/1000.0); + sphPz_h->Fill(pz); + p1x_h->Fill(x1); + p1id_h->Fill(iq1); + + lheF.writeEvent(fileParts,sqrt(Q2)); + myT->Fill(); + NF++; + if(NF%(Nevt/10) == 0) cout << "Produced Event " << NF << " pdfN : " << pdfN << endl; + //cout << "Produced Event " << NF << " pdfN : " << pdfN << endl; + } + + cout << "Max Weight: " << maxwt << endl; + + lheF.close(); + myF->cd(); + x1_h->Write(); + mcTot_h->Write(); + inQid_h->Write(); + frac2D_h->Write(); + p1x_h->Write(); + p1id_h->Write(); + outID_h->Write(); + sphM_h->Write(); + st_h->Write(); + sphPz_h->Write(); + for(int i = 0; i < pt_hv.size(); i++) + { + pt_hv[i].Write(); + eta_hv[i].Write(); + phi_hv[i].Write(); + } + myT->Write(); + myF->Close(); + + + return 0; +}; + + + + + + diff --git a/src/LHEWriter.cc b/src/LHEWriter.cc new file mode 100644 index 0000000..1cbb8e7 --- /dev/null +++ b/src/LHEWriter.cc @@ -0,0 +1,70 @@ +#include "../include/LHEWriter.h" +#include +#include +using namespace std; + +LHEWriter::LHEWriter(string fName) +{ + oF.open((fName+".lhe").c_str()); + oF << std::scientific; + oF << "" << endl; + oF << "
" << endl; + oF << "
" << endl; + oF << "" << endl; + oF << "\t2212 2212 0.65000000000e+04 0.65000000000e+04 292200 292200 292200 292200 3 1" << endl;//beam conditions information + oF << "\t7.3 1.0 1.0 7000" << endl;//sphaleron process information + oF << "" << endl; +} + +LHEWriter::~LHEWriter() +{ +} + +int LHEWriter::writeEvent(vector outParts, double Q) +{ + //cout << "size check: " << decayK.size() << "\t" << decayPIDs.size() << "\t" << decayColz.size() << endl; + oF << "" << endl; + oF << "\t" << outParts.size() << " 7000 1 "<< Q << " 7.3e-03 0.118" << endl; + for(int i = 0; i < int(outParts.size()); i++) + { + oF << "\t" << outParts[i].pid; + + if(i < 2) oF << "\t-1\t"; + else if(fabs(outParts[i].pid) > 1000 || ( i < outParts.size() - 2 && (outParts[i].color == 501 || outParts[i].color == 502) ) ) oF << "\t2\t"; + else oF << "\t1\t"; + oF << outParts[i].m1 << "\t" << outParts[i].m2 << "\t"; + + if(fabs(outParts[i].pid) < 1000) + { + if(outParts[i].pid > 0) oF << outParts[i].color << "\t0\t"; + else oF << "0\t" << outParts[i].color << "\t"; + } + else + { + if(outParts[i].ic > 0) + { + oF << outParts[i].color << "\t0\t"; + } + else + { + oF << "0\t" << outParts[i].color << "\t"; + } + } + + oF << outParts[i].p4v.Px() << "\t"; + oF << outParts[i].p4v.Py() << "\t"; + oF << outParts[i].p4v.Pz() << "\t"; + oF << outParts[i].p4v.E() << "\t"; + oF << outParts[i].p4v.M() << "\t"; + oF << "0\t9" << endl; + } + oF << "" << endl; + return 1; +} + +int LHEWriter::close() +{ + oF << "" << endl; + oF.close(); + return 1; +} diff --git a/src/configBuilder.cc b/src/configBuilder.cc new file mode 100644 index 0000000..ce13aec --- /dev/null +++ b/src/configBuilder.cc @@ -0,0 +1,181 @@ +#include "../include/configBuilder.h" +#include + +using namespace std; + +configBuilder::configBuilder() +{ + //Seed random number generator + rand.SetSeed(0); + + //Build Fermion Doublets, the order of them is SUPER important + particle partBuf; + partBuf.mass = ELE_MASS; + partBuf.pid = ELE_PID; + partBuf.color = 0; + partBuf.q3 = ELE_Q3; + parts.push_back(partBuf); // 0 + + partBuf.mass = 0; + partBuf.pid = ENU_PID; + partBuf.color = 0; + partBuf.q3 = ENU_Q3; + parts.push_back(partBuf); // 1 + + partBuf.mass = DQ_MASS; + partBuf.pid = DQ_PID; + partBuf.color = 1; + partBuf.q3 = DQ_Q3; + parts.push_back(partBuf); // 2 + + partBuf.mass = UQ_MASS; + partBuf.pid = UQ_PID; + partBuf.color = 1; + partBuf.q3 = UQ_Q3; + parts.push_back(partBuf); // 3 + + partBuf.mass = MU_MASS; + partBuf.pid = MU_PID; + partBuf.color = 0; + partBuf.q3 = MU_Q3; + parts.push_back(partBuf); // 4 + + partBuf.mass = 0; + partBuf.pid = MNU_PID; + partBuf.color = 0; + partBuf.q3 = MNU_Q3; + parts.push_back(partBuf); // 5 + + partBuf.mass = SQ_MASS; + partBuf.pid = SQ_PID; + partBuf.color = 1; + partBuf.q3 = SQ_Q3; + parts.push_back(partBuf); // 6 + + partBuf.mass = CQ_MASS; + partBuf.pid = CQ_PID; + partBuf.color = 1; + partBuf.q3 = CQ_Q3; + parts.push_back(partBuf); // 7 + + partBuf.mass = TAU_MASS; + partBuf.pid = TAU_PID; + partBuf.color = 0; + partBuf.q3 = TAU_Q3; + parts.push_back(partBuf); // 8 + + partBuf.mass = 0; + partBuf.pid = TNU_PID; + partBuf.color = 0; + partBuf.q3 = TNU_Q3; + parts.push_back(partBuf); // 9 + + partBuf.mass = BQ_MASS; + partBuf.pid = BQ_PID; + partBuf.color = 1; + partBuf.q3 = BQ_Q3; + parts.push_back(partBuf); // 10 + + partBuf.mass = TQ_MASS; + partBuf.pid = TQ_PID; + partBuf.color = 1; + partBuf.q3 = TQ_Q3; + parts.push_back(partBuf); // 11 +} + +configBuilder::~configBuilder() +{ +} + +vector configBuilder::build(int iQ1, int iQ2, float pNCS, bool bCancel) +{ + vector conf; + vector doublets; + int dCS = 2*int(rand.Uniform() < pNCS) - 1;//Choose matter (1) or antimatter (-1) + for(int i = 0; i < 12; i++) + { + doublets.push_back(i); + } + while(doublets.size() > 0) + { + int j = rand.Integer(doublets.size()); //Choose a random doublet from the list of the remaining doublets + int doub = doublets[j]; + //Calculate index of particle in partBuf, this part conserves charge while also giving the correct number of each doublet + int partI = (doub/4)*4 + 2*int(bool(doub%4)) + (12-doublets.size())%2;//This is the magic part of this whole thing + particle part = parts[partI]; + part.pid = dCS*part.pid; + part.q3 = dCS*part.q3; + conf.push_back(part); + doublets.erase(doublets.begin()+j); + } + + bool spec1 = true; + bool spec2 = true; + + if(bCancel) + { + for(int i = 0; i < 12; i++) + { + if(iQ1 + conf[i].pid == 0 && spec1) + { + conf.erase(conf.begin()+i); + spec1 = false; + } + else if(iQ2 + conf[i].pid == 0 && spec2) + { + conf.erase(conf.begin()+i); + spec2 = false; + } + } + } + if(spec1) + { + //Look up spectator quark particle info and push on configuration + particle specBuf; + for(int i = 0; i < parts.size(); i++) + { + if(fabs(iQ1) == parts[i].pid) + { + specBuf = parts[i]; + specBuf.pid = iQ1; + specBuf.color = 501; + } + } + conf.push_back(specBuf); + } + + if(spec2) + { + //Look up spectator quark particle info and push on configuration + particle specBuf; + for(int i = 0; i < parts.size(); i++) + { + if(fabs(iQ2) == parts[i].pid) + { + specBuf = parts[i]; + specBuf.pid = iQ2; + specBuf.color = 502; + } + } + conf.push_back(specBuf); + } + + /*conf.clear(); + for(int i = 0; i < parts.size(); i++) + { + if(i==1 || i==5 || i==9) continue; + particle partBuf; + partBuf.mass = parts[i].mass; + partBuf.pid = dCS*parts[i].pid; + partBuf.color = parts[i].color; + conf.push_back(partBuf); + if(i == 2 || i == 6 || i == 10) conf.push_back(partBuf); + }*/ + return conf; +} + + + + + + diff --git a/src/particles.cc b/src/particles.cc new file mode 100644 index 0000000..828448a --- /dev/null +++ b/src/particles.cc @@ -0,0 +1,76 @@ +//This class serves as a LUT of the particles that can come from sphaleron transitions +#include "../include/particles.h" +#include +#include + +particleBase::particleBase() +{ + partLUT.clear(); + particle partBuf; + //Put in a blank first so the entry in partLUT is equal to fabs(PID) + partBuf.pid = -1; partBuf.mass = -1.0; partBuf.color = -99; partBuf.q3 = -4; partBuf.p4v = TLorentzVector(0.0,0.0,0.0,0.0); + partLUT.push_back(partBuf); //0 + //1 Down Quark, the p4v will always stay zero at this stage, as the color will also be established later + //Also I only store the matter version of particles and use getParticle to make the correction if antimatter is requested (negative pid) + partBuf.pid = DQ_PID; partBuf.mass = DQ_MASS; partBuf.q3 = DQ_Q3; + partLUT.push_back(partBuf); + //2 Up Quark + partBuf.pid = UQ_PID; partBuf.mass = UQ_MASS; partBuf.q3 = UQ_Q3; + partLUT.push_back(partBuf); + //3 Strange Quark + partBuf.pid = SQ_PID; partBuf.mass = SQ_MASS; partBuf.q3 = SQ_Q3; + partLUT.push_back(partBuf); + //4 Charm Quark + partBuf.pid = CQ_PID; partBuf.mass = CQ_MASS; partBuf.q3 = CQ_Q3; + partLUT.push_back(partBuf); + //5 Bottom Quark + partBuf.pid = BQ_PID; partBuf.mass = BQ_MASS; partBuf.q3 = BQ_Q3; + partLUT.push_back(partBuf); + //6 Top Quark + partBuf.pid = TQ_PID; partBuf.mass = TQ_MASS; partBuf.q3 = TQ_Q3; + partLUT.push_back(partBuf); + //Put in some more blanks to get to ELE_PID + partBuf.pid = -1; partBuf.mass = -1.0; partBuf.q3 = -4; + partLUT.push_back(partBuf); //7 + partLUT.push_back(partBuf); //8 + partLUT.push_back(partBuf); //9 + partLUT.push_back(partBuf); //10 + //11 Electron + partBuf.pid = ELE_PID; partBuf.mass = ELE_MASS; partBuf.q3 = ELE_Q3; + partLUT.push_back(partBuf); + //12 Electron Neutrino + partBuf.pid = ENU_PID; partBuf.mass = ENU_MASS; partBuf.q3 = ENU_Q3; + partLUT.push_back(partBuf); + //13 Muon + partBuf.pid = MU_PID; partBuf.mass = MU_MASS; partBuf.q3 = MU_Q3; + partLUT.push_back(partBuf); + //14 Muon Neutrino + partBuf.pid = MNU_PID; partBuf.mass = MNU_MASS; partBuf.q3 = MNU_Q3; + partLUT.push_back(partBuf); + //15 Tau + partBuf.pid = TAU_PID; partBuf.mass = TAU_MASS; partBuf.q3 = TAU_Q3; + partLUT.push_back(partBuf); + //16 Tau Neutrino + partBuf.pid = TNU_PID; partBuf.mass = TNU_MASS; partBuf.q3 = TNU_Q3; + partLUT.push_back(partBuf); + + +} + +particleBase::~particleBase() +{ +} + +particle particleBase::getParticle(int pid) +{ + int ii = fabs(pid); + int ss = pid/ii; + if(!(ss == 1 || ss == -1)) cout << "Something odd happened in particleBase::getParticle(int pid)" << endl; + particle partBuf = partLUT[ii]; + partBuf.pid = pid; partBuf.q3 = ss*partBuf.q3; + partBuf.m1 = 0; + partBuf.m2 = 0; + partBuf.ic = 0; + return partBuf; + +}