Skip to content

Commit

Permalink
Merge pull request #10 from dglazier/deuteron
Browse files Browse the repository at this point in the history
Deuteron
  • Loading branch information
dglazier authored May 31, 2021
2 parents 33a8943 + b9cf710 commit 42fb9f0
Show file tree
Hide file tree
Showing 55 changed files with 1,391 additions and 532 deletions.
24 changes: 12 additions & 12 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -1,13 +1,17 @@
cmake_minimum_required(VERSION 2.9)
project(elSpectro)


#set(CMAKE_INSTALL_BINDIR ${CMAKE_CURRENT_SOURCE_DIR}/bin)
#set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_INSTALL_BINDIR})
#set(CMAKE_INSTALL_LIBDIR ${CMAKE_CURRENT_SOURCE_DIR}/lib)
#set(CMAKE_INSTALL_INCLUDEDIR ${CMAKE_CURRENT_SOURCE_DIR}/include)
#if -DCMAKE_INSTALL_PREFIX=/install/here use that
IF(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT)
SET(CMAKE_INSTALL_PREFIX ${CMAKE_CURRENT_SOURCE_DIR} CACHE PATH "install to source" FORCE)
ENDIF(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT)

#in case elspectro part of larger package
IF(NOT CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT)
set(CMAKE_INSTALL_BINDIR ${CMAKE_INSTALL_PREFIX}/bin)
set(CMAKE_INSTALL_LIBDIR ${CMAKE_INSTALL_PREFIX}/lib)
ENDIF(NOT CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT)

#if no default install dirs, install at source
if(NOT DEFINED CMAKE_INSTALL_BINDIR)
set(CMAKE_INSTALL_BINDIR ${CMAKE_CURRENT_SOURCE_DIR}/bin)
endif(NOT DEFINED CMAKE_INSTALL_BINDIR)
Expand All @@ -18,15 +22,11 @@ endif(NOT DEFINED CMAKE_INSTALL_LIBDIR)

set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_INSTALL_BINDIR})

message(INFO ${CMAKE_INSTALL_BINDIR} )

set(CMAKE_CXX_FLAGS "-fPIC -O3")
set(CMAKE_CXX_STANDARD 11)
set(CMAKE_CXX_STANDARD_REQUIRED YES)

#include($ENV{HEPMC3}/share/HepMC3/cmake/HepMC3Config.cmake)
#find_package(HepMC3 REQUIRED)

find_package(ROOT REQUIRED MathMore RooFit GenVector EG Minuit2)
find_package(ROOT REQUIRED MathMore RooFit GenVector EG)
list(APPEND CMAKE_PREFIX_PATH $ENV{ROOTSYS})
#---Define useful ROOT functions and macros (e.g. ROOT_GENERATE_DICTIONARY)
include(${ROOT_USE_FILE})
Expand Down
22 changes: 13 additions & 9 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,8 @@ Event Generator framework for incorporating Spectroscopy into electro / photopro

ROOT 6

with MathMore RooFit GenVector EG

We require jpacPhot but this is included as a submodule. You may link to your own version of you prefer.

# Installation
Expand All @@ -19,6 +21,10 @@ We require jpacPhot but this is included as a submodule. You may link to your ow
setenv ELSPECTRO /path/to/elSpectro (or $PWD)
setenv PATH ${PATH}:${ELSPECTRO}/bin

Note, if you want to use jpacPhoto with jpacBox you need to add the C++ boost library location to your path

setenv PATH ${PATH}:/where/is/boost

# Build with cmake

mkdir build; cd build; cmake ../
Expand All @@ -30,28 +36,26 @@ We require jpacPhot but this is included as a submodule. You may link to your ow

cd examples

Note the --i option means retain interactive root session, if not included elsepctro will exit once script is complete.

### Examples comparing to simple weighting of TGenPhaseSpace

1) Decay a rho meson to 2 pi

elspectro ComparePhaseSpaceto2.C
elspectro --i ComparePhaseSpaceto2.C

2) Decay g+p -> rho(pi+,pi-) p

elspectro ComparePhaseSpaceto3Rho.C
elspectro --i ComparePhaseSpaceto3Rho.C

3) Decay g+p -> X(rho(pi+,pi-) , phi(K+,K-) ) p

elspectro ComparePhaseSpaceto5RhoPhi.C
elspectro --i ComparePhaseSpaceto5RhoPhi.C

### Examples of ElectroProduction
4) Decay g+p -> X(rho(pi+,pi-) , 4pi ) p

1) EIC e + p -> e' rho(pi+,pi-) p
elspectro --i ComparePhaseSpaceto4PiRho.C

elspectro RunRhoProton.C
hW.Draw()
hQ2.Draw()
hRhoM.Draw()

### Examples of ElectroProduction of Jpac amplitudes

Expand Down
15 changes: 13 additions & 2 deletions core/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -3,12 +3,17 @@ set(DIRNAME core)
set(BASE $ENV{ELSPECTRO})

include_directories($ENV{JPACPHOTO}/include)
include_directories($ENV{JPACPHOTO}/include/amplitudes)

ROOT_GENERATE_DICTIONARY(G__${ELSPECTRO}
Particle.h
DecayingParticle.h
CollidingParticle.h
DecayModel.h
PhaseSpaceDecay.h
NuclearBreakup.h
SDMEDecay.h
TensorSDMEDecay.h
VectorSDMEDecay.h
DecayModelQ2W.h
DecayModelst.h
Expand All @@ -19,6 +24,7 @@ ROOT_GENERATE_DICTIONARY(G__${ELSPECTRO}
TwoBodyFlat.h
TwoBody_stu.h
ScatteredElectron_xy.h
QuasiFreeNucleon.h
ProductionProcess.h
ElectronScattering.h
Distribution.h
Expand All @@ -31,7 +37,7 @@ ROOT_GENERATE_DICTIONARY(G__${ELSPECTRO}
SDME.h
PhotonPolarisationVector.h
DistVirtPhotFlux_xy.h
CurrentEventInfo.h
# CurrentEventInfo.h
ParticleManager.h
DecayManager.h
MassPhaseSpace.h
Expand All @@ -47,8 +53,12 @@ ROOT_GENERATE_DICTIONARY(G__${ELSPECTRO}
add_library(${ELSPECTRO} SHARED
Particle.cpp
DecayingParticle.cpp
CollidingParticle.cpp
DecayModel.cpp
PhaseSpaceDecay.cpp
NuclearBreakup.cpp
SDMEDecay.cpp
TensorSDMEDecay.cpp
VectorSDMEDecay.cpp
DecayModelQ2W.cpp
DecayModelst.cpp
Expand All @@ -59,6 +69,7 @@ add_library(${ELSPECTRO} SHARED
TwoBodyFlat.cpp
TwoBody_stu.cpp
ScatteredElectron_xy.cpp
QuasiFreeNucleon.cpp
ProductionProcess.cpp
ElectronScattering.cpp
Distribution.cpp
Expand All @@ -69,7 +80,7 @@ add_library(${ELSPECTRO} SHARED
DistVirtPhotFlux_xy.cpp
SDME.cpp
PhotonPolarisationVector.cpp
CurrentEventInfo.cpp
# CurrentEventInfo.cpp
ParticleManager.cpp
DecayManager.cpp
MassPhaseSpace.cpp
Expand Down
95 changes: 95 additions & 0 deletions core/CollidingParticle.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,95 @@
#include "CollidingParticle.h"
#include "FunctionsForGenvector.h"
#include "Manager.h"

namespace elSpectro{

///////////////////////////////////////////////////////////
//cannot default construct at least need a 4-momentum
CollidingParticle::CollidingParticle(int pdg,Double_t momentum):
Particle(pdg){
//Set interacting particle pdg
_interactingPdg=pdg;
auto mass= PdgMass();
LorentzVector lv(0,0,momentum,TMath::Sqrt(momentum*momentum+mass*mass));
SetP4(lv);
_nominal=P4();
_interactingParticle=P4ptr();
}
/////////////////////////////////////////////////////////
//or a model to generate a 4-momentum
CollidingParticle::CollidingParticle(int pdg,Double_t momentum,int parentpdg,DecayModel* model,DecayVectors* decayer):
Particle(parentpdg),_model{model},_decayer{decayer}
{
//Set interacting particle pdg
_interactingPdg=pdg;
//set the parent LorentzVector
auto mass= PdgMass();
LorentzVector lv(0,0,momentum,TMath::Sqrt(momentum*momentum+mass*mass));
SetP4(lv);
_nominal=P4();

//Find the relevent particle pointer in the model
//this is the particle which is used in the production process
UInt_t position=0;
for(auto& p:_model->Products()){

if(p->Pdg()==pdg){

if(_interactingParticle!=nullptr){
std::cerr<<"CollidingParticle::CollidingParticle, multiple particles with pdg = "<<pdg<<std::endl; exit(0);
}
else{
_interactingParticle=p->P4ptr();
}
//Not a stable final state particle!
Manager::Instance().Particles().RemoveStable(p);
_model->SwapProducts(0,position); //make sure interacting particle is first in products vector, must be done after remove
}
else{ //not a beam particle but "spectator"
//remove from stable particle list
//as these will be boosted to lab
//whereas colliding particles and their
//products are already in lab frame

Manager::Instance().Particles().MoveStableToLab(p);

}
position++;
}
//move interacting particle to start of products vector

//We need a "nominal" 4-momentum for our interacting particle
//to do this we boost it from rest into lab frame of parent
_decayer->BoostToParent(P4(),(*_interactingParticle));
}
/////////////////////////////////////////////////////////
/// rotate beam angles
void CollidingParticle::SetAngleThetaPhi(Double_t th,Double_t phi){
_dirTheta=th;_dirPhi=phi;
auto p4=_nominal; //copy 4-vector
//and rotate it
genvector::LorentzRotateY(p4,th);
genvector::LorentzRotateZ(p4,phi);
//Set the rotated vector
SetP4(p4);
_nominal=p4;
}
/////////////////////////////////////////////////////////
void CollidingParticle::PostInit(ReactionInfo* info){
//decay vertex position

if(_model){
auto& products=_model->Products();
//same vertex as parent
for(auto* prod: products){
prod->SetVertex(VertexID(),VertexPosition());
}

_model->PostInit(info);
}
if(_decayer)_decayer->PostInit(info);

}

}
93 changes: 93 additions & 0 deletions core/CollidingParticle.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,93 @@
//////////////////////////////////////////////////////////////
///
///Class: CollidingParticle
///Description:
/// 1) Class for generating initial state LorentzVector
///


#pragma once

#include "Particle.h"
#include "DecayModel.h"
#include "DecayVectors.h"

namespace elSpectro{


class CollidingParticle : public Particle {

public:

CollidingParticle()=delete;
//cannot default construct at least need a 4-momentum
CollidingParticle(int pdg,Double_t momentum);
//or a model to generate a 4-momentum
CollidingParticle(int pdg,Double_t momentum,int parentpdg,DecayModel* model,DecayVectors* decayer);

DecayModel* Model()const {return _model;}

double Generate(){
if(_decayer.get()==nullptr) return 1.0;
return _decayer->Generate(P4(),_model->Products());
}

const DecayVectors* Decayer() const {return _decayer.get();}
void SetDecayer(DecayVectors* decayer){_decayer.reset(decayer);}

Double_t GenerateComponents(){
if(_model==nullptr )return 1.0;

Double_t weight=0;
while((weight)==0){ //sample until physical
weight=Generate();
weight*=_model->Intensity();
}
// can add beam divergence etc. here
// ...
return weight;
}
DecayType IsDecay() const noexcept final {return DecayType::Production;}

void PostInit(ReactionInfo* info) ;

const LorentzVector* GetInteracting4Vector() const {return _interactingParticle;}
Int_t GetInteractingPdg()const {return _interactingPdg;}

/*void SetVertexXYZT(double x,double y,double z,double t){
_productionVertex.SetXYZT(x,y,z,t);
}*/


void SetAngleThetaPhi(Double_t th,Double_t phi);
void SetHorSize(Double_t s){_sizeHor=s;}
void SetVerSize(Double_t s){_sizeVer=s;}
void SetHorDivergence(Double_t s){_divHor=s;}
void SetVerDivergence(Double_t s){_divVer=s;}


private:

DecayModel* _model={nullptr}; //not owner

std::unique_ptr<DecayVectors> _decayer={nullptr}; //owner

LorentzVector *_interactingParticle={nullptr}; //not owner;

LorentzVector _nominal; //nominal beam 4-momentum
// LorentzVector _productionVertex;
//int _prodVertexID={0};
int _interactingPdg={0};
Double_t _dirTheta={0};
Double_t _dirPhi={0};
Double_t _sizeHor={0};
Double_t _sizeVer={0};
Double_t _divHor={0};
Double_t _divVer={0};


ClassDefOverride(elSpectro::CollidingParticle,1); //class CollidingParticle

};//class CollidingParticle

}//namespace elSpectro
Loading

0 comments on commit 42fb9f0

Please sign in to comment.