From c80363e2b27070f88cf063faa6d108776d1361c5 Mon Sep 17 00:00:00 2001 From: dglazier Date: Mon, 8 Mar 2021 16:07:17 +0000 Subject: [PATCH 01/21] give integrator name dependent on beam energy --- core/ElectronScattering.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/core/ElectronScattering.cpp b/core/ElectronScattering.cpp index e84d6ea..3c18df9 100644 --- a/core/ElectronScattering.cpp +++ b/core/ElectronScattering.cpp @@ -388,12 +388,12 @@ namespace elSpectro{ auto wrapPdf=ROOT::Math::Functor( flnXlnYt , 3); - auto pdf = RooFunctorPdfBinding("ElScatterIntegral", "ElScatterIntegral", wrapPdf, RooArgList(xvar,yvar,cthvar)); + auto pdf = RooFunctorPdfBinding(Form("ElScatterIntegral%lf",Eel), "ElScatterIntegral", wrapPdf, RooArgList(xvar,yvar,cthvar)); auto roovars= RooArgSet(xvar,yvar,cthvar); gBenchmark->Start("RooFitIntegral"); - - auto RFintegral=pdf.getNorm(roovars); + + auto RFintegral=pdf.getNorm(roovars); gBenchmark->Stop("RooFitIntegral"); gBenchmark->Print("RooFitIntegral"); From c8f1246f29fae335eab7a557e05f888e489b5cfe Mon Sep 17 00:00:00 2001 From: dglazier Date: Tue, 9 Mar 2021 14:26:34 +0000 Subject: [PATCH 02/21] protect from below threshold amplitude requests --- core/JpacModelst.h | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/core/JpacModelst.h b/core/JpacModelst.h index fbe849c..36ea57e 100644 --- a/core/JpacModelst.h +++ b/core/JpacModelst.h @@ -36,8 +36,9 @@ namespace elSpectro{ double MatrixElementsSquared_T() const override { _amp->_kinematics->set_mX( GetMeson()->Mass() ); // _amp->_kinematics->set_Q2( get_Q2() ); - // std::cout<<"me "<Mass()<<" Q2 "<< get_Q2()<<" t "<probability_distribution(get_s(),get_t())/4;// Average over initial state helicites; + //std::cout<<"me "<Mass()<<" Q2 "<< get_Q2()<<" t "<_kinematics->Wth()<_kinematics->Wth()) return 0; + return _amp->probability_distribution(get_s(),get_t())/4;// Average over initial state helicites; } void CalcMesonSDMEs() const override ; From 240e4d4d25cbec76f122819e2c9be11f21758052 Mon Sep 17 00:00:00 2001 From: dglazier Date: Tue, 20 Apr 2021 13:07:54 +0100 Subject: [PATCH 03/21] give install prefix, remove c++11 requirement, will be taken from ROOT install anyway --- CMakeLists.txt | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index fb3477a..60d33aa 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,12 +2,14 @@ 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) #in case elspectro part of larger package + +if(DEFINED CMAKE_INSTALL_PREFIX) + set(CMAKE_INSTALL_BINDIR ${CMAKE_INSTALL_PREFIX}/bin) + set(CMAKE_INSTALL_LIBDIR ${CMAKE_INSTALL_PREFIX}/lib) +endif(DEFINED CMAKE_INSTALL_PREFIX) + if(NOT DEFINED CMAKE_INSTALL_BINDIR) set(CMAKE_INSTALL_BINDIR ${CMAKE_CURRENT_SOURCE_DIR}/bin) endif(NOT DEFINED CMAKE_INSTALL_BINDIR) @@ -20,13 +22,8 @@ set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${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}) From e67b70dbef99b286902277a4fbbc52724b07fe9e Mon Sep 17 00:00:00 2001 From: dglazier Date: Tue, 20 Apr 2021 13:45:15 +0100 Subject: [PATCH 04/21] give install prefix, remove c++11 requirement, will be taken from ROOT install anyway --- CMakeLists.txt | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 60d33aa..a9662e8 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,15 +1,15 @@ cmake_minimum_required(VERSION 2.9) -project(elSpectro) - - +SET(CMAKE_INSTALL_PREFIX "DEFAULT") #do not /usr/local as default -#in case elspectro part of larger package +project(elSpectro) -if(DEFINED CMAKE_INSTALL_PREFIX) +#if -DCMAKE_INSTALL_PREFIX=/install/here use that +if(NOT CMAKE_INSTALL_PREFIX STREQUAL "DEFAULT") set(CMAKE_INSTALL_BINDIR ${CMAKE_INSTALL_PREFIX}/bin) set(CMAKE_INSTALL_LIBDIR ${CMAKE_INSTALL_PREFIX}/lib) -endif(DEFINED CMAKE_INSTALL_PREFIX) +endif(NOT CMAKE_INSTALL_PREFIX STREQUAL "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) @@ -20,6 +20,7 @@ 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") From 0ad971c96286e3204c3cd5d1fb9aaac05aa86284 Mon Sep 17 00:00:00 2001 From: dglazier Date: Tue, 20 Apr 2021 16:07:47 +0100 Subject: [PATCH 05/21] various edits, including renaming ElectronScattering integral, interactve mode for elspectro, ... --- core/DecayModel.cpp | 24 ++++++++++++++---------- core/DecayModel.h | 1 + core/DecayModelQ2W.cpp | 3 ++- core/DecayModelst.cpp | 2 +- core/DecayingParticle.cpp | 4 ++-- core/ElectronScattering.cpp | 27 +++++++++++++++++++-------- core/Manager.h | 4 +++- core/MassPhaseSpace.h | 2 +- core/ParticleManager.h | 3 +-- core/PhaseSpaceDecay.cpp | 2 +- core/ScatteredElectron_xy.cpp | 2 ++ core/VectorSDMEDecay.cpp | 2 +- core/src/elspectro.cpp | 14 ++++++++++---- 13 files changed, 58 insertions(+), 32 deletions(-) diff --git a/core/DecayModel.cpp b/core/DecayModel.cpp index c940c01..bee7873 100644 --- a/core/DecayModel.cpp +++ b/core/DecayModel.cpp @@ -52,7 +52,7 @@ namespace elSpectro{ p->PostInit(info); prodmasses.push_back(p->MinimumMassPossible()); } - + for(uint i=0;i<_unstables.size();++i){ //for each product //calculate mass reserved for subsequent unstable products float reserve=0; @@ -63,6 +63,7 @@ namespace elSpectro{ } } + void DecayModel::GetStableMasses( std::vector& masses) const{ @@ -85,9 +86,9 @@ namespace elSpectro{ } double DecayModel::PhaseSpaceWeightSq(double W){ //if(Parent()->Pdg()==-2211)std::cout<Mass(); } - + uint iu=0; //synch reserved mass vector for(auto* p:_unstables){ //This should be the only call to DetermineDynamicMass in the code. //Unless somewhere else uses LockMass in which case //this call to DetermineDynamicMass will not change its value - + //Note in case there are additional unstable particle we //must subtract off their minimum masses p->DetermineDynamicMass(-1,TCM-_unstableReservedMass[iu++]); TCM-=p->Mass(); - if(TCM<0){ std::cout<<"DecayModel::PhaseSpaceWeightSq "<Pdg()<<" "<Mass()<<" "<Pdg()<<" "<Mass()<<" "<Pdg()<<" "<Mass(),_products[1]->Mass()); - + + result *= kine::PDK2(W,_products[0]->Mass(),_products[1]->Mass()); + for(auto* p:_unstables) result*=p->PhaseSpaceWeightSq(); diff --git a/core/DecayModel.h b/core/DecayModel.h index 7c40334..af5bb01 100644 --- a/core/DecayModel.h +++ b/core/DecayModel.h @@ -49,6 +49,7 @@ namespace elSpectro{ virtual double Intensity() const=0; const particle_ptrs& Products() const{ return _products;} + const Particle* Product(UInt_t i) const{ return _products[i];} const decaying_ptrs& UnstableProducts() const{ return _unstables;} const particle_ptrs& StableProducts() const{ return _stables;} diff --git a/core/DecayModelQ2W.cpp b/core/DecayModelQ2W.cpp index 85c9443..17a4161 100644 --- a/core/DecayModelQ2W.cpp +++ b/core/DecayModelQ2W.cpp @@ -88,7 +88,8 @@ namespace elSpectro{ //////////////////////////////////////////////////////// double DecayModelQ2W::Intensity() const{ - // std::cout<<"DecayModelQ2W::Intensity "<P4().E() <P4().E() <GetBinCenter(ih); - + double val_at_t = F(tval)*(tmin-tmax); if(val_at_t>max_at_W) max_at_W=val_at_t; diff --git a/core/DecayingParticle.cpp b/core/DecayingParticle.cpp index 5f74b76..fa6e902 100644 --- a/core/DecayingParticle.cpp +++ b/core/DecayingParticle.cpp @@ -53,8 +53,8 @@ namespace elSpectro{ if(_decay)_decay->PostInit(info); if(_decayer)_decayer->PostInit(info); - std::cout<<"DecayingParticle::PostInit pdg "<(mutableDecayer()); - auto xvar = RooRealVar("lnxIntegral","lnxIntegral",TMath::Exp(photonFlux->Dist().GetMaxLnX()),TMath::Exp(photonFlux->Dist().GetMinLnX()),TMath::Exp(photonFlux->Dist().GetMaxLnX()),""); - auto yvar = RooRealVar("lnyIntegral","lnyIntegral",TMath::Exp(photonFlux->Dist().GetMaxLnY()),TMath::Exp(photonFlux->Dist().GetMinLnY()),TMath::Exp(photonFlux->Dist().GetMaxLnY()),""); + auto xvar = RooRealVar(Form("xIntegral%lf_%lf",photonFlux->Dist().GetMaxLnX(),photonFlux->Dist().GetMaxLnX()),"xIntegral",TMath::Exp(photonFlux->Dist().GetMaxLnX()),TMath::Exp(photonFlux->Dist().GetMinLnX()),TMath::Exp(photonFlux->Dist().GetMaxLnX()),""); + auto yvar = RooRealVar(Form("yIntegral%lf_%lf",photonFlux->Dist().GetMaxLnY(),photonFlux->Dist().GetMaxLnY()),"yIntegral",TMath::Exp(photonFlux->Dist().GetMaxLnY()),TMath::Exp(photonFlux->Dist().GetMinLnY()),TMath::Exp(photonFlux->Dist().GetMaxLnY()),""); auto cthvar = RooRealVar("CosThIntegral","CosThIntegral",-0.,-0.99999,0.99999,""); + xvar.Print(); + yvar.Print(); + cthvar.Print(); auto gStarModel =dynamic_cast(_gStarN->Model()); auto Q2WModel =dynamic_cast(Model()); @@ -324,13 +328,17 @@ namespace elSpectro{ photonFlux->Dist().SetWThresholdVal(gStarModel->GetMeson()->PdgMass()+gStarModel->GetBaryon()->PdgMass()); auto Eel=_nuclRestElec.E(); - auto flnXlnYt = [this,&photonFlux,&gStarModel,&Q2WModel,&Eel](const double *x) + + double_t threshW= gStarModel->GetMeson()->PdgMass()+gStarModel->GetBaryon()->PdgMass(); + std::cout<<"Eel "<Dist().Eval(x); if(TMath::IsNaN(val)) return 0.; if(val==0) return 0.; + //std::cout<<"a fXYcosth "<GenerateGivenXandY(P4(),Model()->Products(),x[0],x[1]); @@ -339,8 +347,10 @@ namespace elSpectro{ //auto W = gStarModel->Parent()->P4().M();//TMath::Sqrt(escat::M2_pr()+2*Eel*escat::M_pr()*x[1]-escat::Q2_xy( Eel,x[0],x[1])); //get value of dsigma/ds/dcosth cross section at x,y,costh val*=gStarModel->dsigma_costh(x[2]); + // std::cout<<"b fXYcosth "<Parent()->P4().M()<<" thesh"<Q2H1Rho(); return val; }; @@ -386,17 +396,18 @@ namespace elSpectro{ gBenchmark->Print("RooFitIntegralW"); */ - auto wrapPdf=ROOT::Math::Functor( flnXlnYt , 3); - + auto wrapPdf=ROOT::Math::Functor( fXYcosth , 3); + //add Eel so it recaluclate normalisation, otherwise RooFit returns cache + //auto pdf = RooFunctorPdfBinding("ElScatterIntegral", "ElScatterIntegral", wrapPdf, RooArgList(xvar,yvar,cthvar)); auto pdf = RooFunctorPdfBinding(Form("ElScatterIntegral%lf",Eel), "ElScatterIntegral", wrapPdf, RooArgList(xvar,yvar,cthvar)); auto roovars= RooArgSet(xvar,yvar,cthvar); gBenchmark->Start("RooFitIntegral"); - auto RFintegral=pdf.getNorm(roovars); + auto RFintegral=pdf.getNorm(roovars); gBenchmark->Stop("RooFitIntegral"); gBenchmark->Print("RooFitIntegral"); - + std::cout<<" ElectronScattering::IntegrateCrossSection() "<Dist().Integral()<<" nb"<IntegrateCrossSection(); _nEventsToGen=n_or_lum*1E-33*beamtime*_integralXSection*Reaction()->BranchingFraction();//1E-33(cm2tonb) std::cout<<"Manager::SetNEvents_via_LuminosityTime , going to generate "<<_nEventsToGen<<" events"<Uniform()*max ){ //reject this combintation - if(wee==0) std::cout<<"ps "<_sampledMax){ _sampledMax = wee; diff --git a/core/ParticleManager.h b/core/ParticleManager.h index 5cf4279..348497c 100644 --- a/core/ParticleManager.h +++ b/core/ParticleManager.h @@ -47,9 +47,8 @@ namespace elSpectro{ } int RegisterNewPdgParticle(double nominalMass,Distribution* dist=nullptr){ - std::cout<<"RegisterNewPdgParticle "<MassDistribution()==nullptr){//production proces does not have mass distribution - if(dynamic_cast(Parent())==nullptr){ + if(dynamic_cast(Parent())==nullptr&&Parent()->Pdg()!=-2211){ std::cerr<<"PhaseSpaceDecay::PostInit parent needs a mass distribution for pdg = "<Pdg(); std::cerr<<"\n you need to use \n >> mass_distribution(PDG,new DistTF1{TF1(\"massDist\",\"1\",MINMASS,MAXMASS)});"; std::cerr<<" \n where PDG (9995-9999) is the pdg number you assigned the decaying particle, and MINMAMSS and MAXMASS is the mass limits it will be allowed to have, for pure phase space this must be at least the kinematically allowed range"; diff --git a/core/ScatteredElectron_xy.cpp b/core/ScatteredElectron_xy.cpp index c21da26..4d3c8cc 100644 --- a/core/ScatteredElectron_xy.cpp +++ b/core/ScatteredElectron_xy.cpp @@ -132,7 +132,9 @@ namespace elSpectro{ _scattered.SetXYZT(x_sc,y_sc,z_sc,Esc);//scattered electron //Must make sure scattered e- is in the same frame as the parent //still in rest system of nucl, just need rotation + // std::cout<Pdg()==11){ products[0]->SetP4(_scattered); products[1]->SetP4(parent -_scattered); diff --git a/core/VectorSDMEDecay.cpp b/core/VectorSDMEDecay.cpp index fbf4d3f..81c537b 100644 --- a/core/VectorSDMEDecay.cpp +++ b/core/VectorSDMEDecay.cpp @@ -97,7 +97,7 @@ namespace elSpectro{ } std::cout<<_rho->Re(0,0,0)<<" 010 "<<_rho->Re(0,1,0)<<" 01-1 "<<_rho->Re(0,1,-1)<<" 111 "<<_rho->Re(1,1,1)<<" 100 "<<_rho->Re(1,0,0)<<" 110 "<<_rho->Re(1,1,0)<<" 11-1 "<<_rho->Re(1,1,-1)<<" 210 "<<_rho->Re(2,1,0)<<" 21-1 "<<_rho->Re(2,1,-1)<<" "<Epsilon()<<" delta " <<_photonPol->Delta()*TMath::RadToDeg()<<" phi " <<_photonPol->Phi()<Re(0,0,0)<=1) << (_rho->Re(0,0,0)>=0) <Re(0,0,0)<=1) <<" "<< (_rho->Re(0,0,0)>=0) <Re(0,1,-1))<=0.5*(1-_rho->Re(0,0,0)) )<Re(0,1,0)*_rho->Re(0,1,0) <=0.25*_rho->Re(0,0,0)*(2-_rho->Re(0,0,0)-_rho->Re(0,1,-1)) )<Re(1,0,0))<=_rho->Re(0,0,0)) <ProcessLine(".x $ELSPECTRO/core/src/Load.C"); - app->ProcessLine(Form(".x %s",macroName.Data())); - - app->Terminate(0); + + if(isInteractive == true) + app->Run(); + else{ + app->ProcessLine(Form(".x %s",macroName.Data())); + app->Terminate(0); + } return 0; From 71d2f182a86c94df9c020c583facae9730788673 Mon Sep 17 00:00:00 2001 From: dglazier Date: Tue, 20 Apr 2021 16:45:54 +0100 Subject: [PATCH 06/21] fix examples --- examples/ComparePhaseSpaceto2.C | 17 ++++++++--------- examples/ComparePhaseSpaceto3Rho.C | 8 ++++---- examples/ComparePhaseSpaceto4PiRho.C | 21 +++++++++------------ examples/ComparePhaseSpaceto5RhoPhi.C | 14 +++++++------- 4 files changed, 28 insertions(+), 32 deletions(-) diff --git a/examples/ComparePhaseSpaceto2.C b/examples/ComparePhaseSpaceto2.C index ce8de6e..3bfcd62 100644 --- a/examples/ComparePhaseSpaceto2.C +++ b/examples/ComparePhaseSpaceto2.C @@ -5,7 +5,7 @@ void ComparePhaseSpaceto2() { Double_t masses[2] = {0.139, 0.139} ; TGenPhaseSpace event; - TLorentzVector W(0,0,0.5,sqrt(0.78*0.78+0.5*0.5)); + TLorentzVector W(0,0,0.5,sqrt(0.775*0.775+0.5*0.5)); event.SetDecay(W, 2, masses); @@ -13,7 +13,7 @@ void ComparePhaseSpaceto2() { TH1F* h2p= new TH1F("pip","#pi+ momentum",1000,0,1); gBenchmark->Start("tgenphasespace"); - for (Int_t n=0;n<1E6;n++) { + for (Int_t n=0;n<1E7;n++) { Double_t weight = event.Generate(); @@ -35,23 +35,22 @@ void ComparePhaseSpaceto2() { using namespace elSpectro; elSpectro::Manager::Instance(); - //auto elSpectroPhaseSpace = model(new PhaseSpaceDecay{ {},{211,-211} }); - auto rho=dynamic_cast( particle(113,model(new PhaseSpaceDecay{{},{211,-211}}))); + auto rho=dynamic_cast( particle(113,model(new PhaseSpaceDecay{{},{211,-211}}))); generator().SetModelForMassPhaseSpace(rho->Model()); rho->SetXYZ(0,0,0.5); - auto pionp = rho->Model()->Products()[0]; - auto pionm = rho->Model()->Products()[1]; + auto pionp = rho->Model()->Product(0); + auto pionm = rho->Model()->Product(1); - TH1F* hel=new TH1F("RhoM","#rho Mass",1000,0.2,1.2); + TH1F* hel=new TH1F("RhoMelS","#rho Mass",1000,0.2,1.2); hel->SetLineColor(2); - TH1F* help= new TH1F("pip","#pi+ momentum",1000,0,1); + TH1F* help= new TH1F("pipelS","#pi+ momentum",1000,0,1); help->SetLineColor(2); gBenchmark->Start("elspectro"); - for (Int_t n=0;n<1E6;n++) { + for (Int_t n=0;n<1E7;n++) { rho->GenerateProducts(); auto rho_pions=pionp->P4()+pionm->P4(); hel->Fill(rho_pions.M()); diff --git a/examples/ComparePhaseSpaceto3Rho.C b/examples/ComparePhaseSpaceto3Rho.C index 4bde64b..fd9a4a4 100644 --- a/examples/ComparePhaseSpaceto3Rho.C +++ b/examples/ComparePhaseSpaceto3Rho.C @@ -77,13 +77,13 @@ void ComparePhaseSpaceto3Rho() { prho->SetXYZT(W.X(),W.Y(),W.Z(),W.T()); - auto pionp = rho->Model()->Products()[0]; - auto pionm = rho->Model()->Products()[1]; - auto proton = prho->Model()->Products()[1]; + auto pionp = rho->Model()->Product(0); + auto pionm = rho->Model()->Product(1); + auto proton = prho->Model()->Product(1); gBenchmark->Start("elspectro"); - + prho->PostInit(nullptr); for (Int_t n=0;nGenerateProducts(); auto p4rho =pionp->P4()+pionm->P4(); diff --git a/examples/ComparePhaseSpaceto4PiRho.C b/examples/ComparePhaseSpaceto4PiRho.C index 90179c2..7e6ae6a 100644 --- a/examples/ComparePhaseSpaceto4PiRho.C +++ b/examples/ComparePhaseSpaceto4PiRho.C @@ -16,7 +16,7 @@ void ComparePhaseSpaceto4PiRho() { TH2F *h2 = new TH2F("h2","h2", 50,1.1,20.8, 50,1.1,20.8); TH2F *h2weight = new TH2F("h2weight","h2weight",100,0,10.2,100,0,0.5); TH1F* h1 = new TH1F("h1","h1", 100,0,5); - TH1F* hMeson = new TH1F("hMeson","Meson Mass", 1000,0,10); + TH1F* hMeson = new TH1F("hMeson","Meson Mass", 100,0,10); TH1F* h1un = new TH1F("h1un","h1un", 100,0,5); h1un->SetLineColor(3); @@ -26,7 +26,7 @@ void ComparePhaseSpaceto4PiRho() { auto bwmax=bw.GetMaximum(); gBenchmark->Start("tgenphasespace"); Long64_t counter=0; - Long64_t Nevents=5000; + Long64_t Nevents=500; for (Int_t n=0;;n++) { @@ -73,7 +73,7 @@ void ComparePhaseSpaceto4PiRho() { hel2->SetMarkerColor(2); TH1F* hel1 = new TH1F("hel1","h1", 100,0,5); hel1->SetLineColor(2); - TH1F* hMeson2 = new TH1F("hMeson2","Meson Mass", 1000,0,10); + TH1F* hMeson2 = new TH1F("hMeson2","Meson Mass", 100,0,10); hMeson2->SetLineColor(2); using namespace elSpectro; @@ -99,17 +99,14 @@ void ComparePhaseSpaceto4PiRho() { prho->SetXYZT(W.X(),W.Y(),W.Z(),W.T()); - //auto pionp = prho->Model()->Products()[3]; - //auto pionm = prho->Model()->Products()[4]; - auto pionp = rho->Model()->Products()[0]; - auto pionm = rho->Model()->Products()[1]; + auto pionp = rho->Model()->Product(0); + auto pionm = rho->Model()->Product(1); Particle* proton =nullptr; //Particle* pionp =nullptr; //Particle* pionm =nullptr; - - for(auto& pp:particles().StableParticles()){ + for(auto& pp:particles().StableParticles()){ if(pp->Pdg()==2212){ proton=pp; } @@ -120,9 +117,9 @@ void ComparePhaseSpaceto4PiRho() { pionm=pp; }*/ } - - - gBenchmark->Start("elspectro"); + + prho->PostInit(nullptr); + gBenchmark->Start("elspectro"); double scale = 100; diff --git a/examples/ComparePhaseSpaceto5RhoPhi.C b/examples/ComparePhaseSpaceto5RhoPhi.C index 35a20de..69f4cac 100644 --- a/examples/ComparePhaseSpaceto5RhoPhi.C +++ b/examples/ComparePhaseSpaceto5RhoPhi.C @@ -15,7 +15,7 @@ void ComparePhaseSpaceto5RhoPhi() { TH1F* h1rho = new TH1F("hrho","hrho", 100,0,1.5); TH1F* h1phi = new TH1F("hphi","hphi", 100,0.9,1.5); - TH1F* h1X = new TH1F("hX","hrho", 100,0,50); + TH1F* h1X = new TH1F("hX","hX", 100,0,5); auto bwRho = TF1("hh","TMath::BreitWigner(x,0.78,0.1)",0.3,3.5); @@ -76,7 +76,7 @@ void ComparePhaseSpaceto5RhoPhi() { h1rho_el->SetLineColor(2); TH1F* h1phi_el = new TH1F("hphi_el","hphi", 100,0.9,1.5); h1phi_el->SetLineColor(2); - TH1F* h1X_el = new TH1F("hX_el","hX", 100,0,50); + TH1F* h1X_el = new TH1F("hX_el","hX", 100,0,5); h1X_el->SetLineColor(2); mass_distribution(113,new DistTF1{TF1("Mrho","TMath::BreitWigner(x,0.78,0.1)",0.2,3.5)}); @@ -93,15 +93,15 @@ void ComparePhaseSpaceto5RhoPhi() { pX->SetXYZT(W.X(),W.Y(),W.Z(),W.T()); - // pX->Print(); - auto pionp = rho->Model()->Products()[0]; - auto pionm = rho->Model()->Products()[1]; - auto Kp = phi->Model()->Products()[0]; - auto Km = phi->Model()->Products()[1]; + auto pionp = rho->Model()->Product(0); + auto pionm = rho->Model()->Product(1); + auto Kp = phi->Model()->Product(0); + auto Km = phi->Model()->Product(1); auto proton = pX->Model()->Products()[1]; + pX->PostInit(nullptr); gBenchmark->Start("elspectro"); int extra=1000; From ef03ebf5ca6df387c188c307b2d9d1e36d22cea0 Mon Sep 17 00:00:00 2001 From: dglazier Date: Thu, 22 Apr 2021 10:22:48 +0100 Subject: [PATCH 07/21] fix crossing angle phi rotation and add counter to ElectronScattering Integration to rpevent RooFit using cache --- core/CMakeLists.txt | 1 + core/DecayVectors.h | 2 +- core/ElectronScattering.cpp | 29 ++++++++++++++++++++--------- core/ElectronScattering.h | 11 ++++++++--- core/ScatteredElectron_xy.cpp | 24 +++++++++++++++++------- core/ScatteredElectron_xy.h | 15 ++++++++------- 6 files changed, 55 insertions(+), 27 deletions(-) diff --git a/core/CMakeLists.txt b/core/CMakeLists.txt index 37da615..bedfc57 100644 --- a/core/CMakeLists.txt +++ b/core/CMakeLists.txt @@ -3,6 +3,7 @@ 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 diff --git a/core/DecayVectors.h b/core/DecayVectors.h index f664233..d571bfa 100644 --- a/core/DecayVectors.h +++ b/core/DecayVectors.h @@ -97,7 +97,7 @@ namespace elSpectro{ //Apply random phi angle now z-axis is in correct direction _rotateAroundZaxis.SetAngle(RandomPhi()); - + child=_rotateToZaxis * child; child=_rotateAroundZaxis * child; diff --git a/core/ElectronScattering.cpp b/core/ElectronScattering.cpp index a680868..7834856 100644 --- a/core/ElectronScattering.cpp +++ b/core/ElectronScattering.cpp @@ -14,7 +14,7 @@ #include namespace elSpectro{ - + int ElectronScattering::NintegralsElectronScattering=0; // ElectronScattering::ElectronScattering(double ep,double ionp, // double eangle,double ionangle, // DecayModel* model, int ionpdg): @@ -343,15 +343,16 @@ namespace elSpectro{ //calculate scatered electron at x and y photonFlux->GenerateGivenXandY(P4(),Model()->Products(),x[0],x[1]); //calculate virtual photon - Model()->Intensity(); + Q2WModel->Intensity(); + //std::cout<<"ElectronScattering integrate g*p "<Parent()->P4()<Parent()->P4().M()<Parent()->P4().M();//TMath::Sqrt(escat::M2_pr()+2*Eel*escat::M_pr()*x[1]-escat::Q2_xy( Eel,x[0],x[1])); - //get value of dsigma/ds/dcosth cross section at x,y,costh + //get value of dsigma(s)/dcosth cross section at x,y,costh val*=gStarModel->dsigma_costh(x[2]); - // std::cout<<"b fXYcosth "<Parent()->P4().M()<<" thesh"<Parent()->P4().M()<<" thesh"<getQ2()<Q2H1Rho(); + val*=Q2WModel->Q2H1Rho(); return val; }; @@ -399,17 +400,27 @@ namespace elSpectro{ auto wrapPdf=ROOT::Math::Functor( fXYcosth , 3); //add Eel so it recaluclate normalisation, otherwise RooFit returns cache //auto pdf = RooFunctorPdfBinding("ElScatterIntegral", "ElScatterIntegral", wrapPdf, RooArgList(xvar,yvar,cthvar)); - auto pdf = RooFunctorPdfBinding(Form("ElScatterIntegral%lf",Eel), "ElScatterIntegral", wrapPdf, RooArgList(xvar,yvar,cthvar)); - auto roovars= RooArgSet(xvar,yvar,cthvar); + TString pdfname(Form("ElScatterIntegral%d",NintegralsElectronScattering)); + auto pdf = new RooFunctorPdfBinding(pdfname, "ElScatterIntegral", wrapPdf, RooArgList(xvar,yvar,cthvar)); + if(_cacheIntegrals==0) NintegralsElectronScattering++;//work around RooFit agressive caching! + pdf->Print(); + std::cout<<_cacheIntegrals<<" "<setDirtyInhibit(true); + + gBenchmark->Start("RooFitIntegral"); - auto RFintegral=pdf.getNorm(roovars); + auto RFintegral=pdf->getNorm(roovars); + delete pdf; + //auto RFintegral=integralRoo->getValV();//(roovars); gBenchmark->Stop("RooFitIntegral"); gBenchmark->Print("RooFitIntegral"); std::cout<<" ElectronScattering::IntegrateCrossSection() "<Dist().Integral()<<" nb"<Pdg()==11){ products[0]->SetP4(_scattered); products[1]->SetP4(parent -_scattered); diff --git a/core/ScatteredElectron_xy.h b/core/ScatteredElectron_xy.h index 220aabf..0a967ea 100644 --- a/core/ScatteredElectron_xy.h +++ b/core/ScatteredElectron_xy.h @@ -9,7 +9,7 @@ #include "DecayVectors.h" #include "DistVirtPhotFlux_xy.h" #include //for boosts etc. -#include +#include #include #include @@ -40,7 +40,7 @@ namespace elSpectro{ protected: - void RotateZaxisToCMDirection(const LorentzVector& parent); + void RotateZaxisToCMDirection(const LorentzVector& parent, LorentzVector& child); double CompleteGivenXandY(const LorentzVector& parent, const particle_ptrs& products,double xx, double yy); private: @@ -51,9 +51,9 @@ namespace elSpectro{ LorentzVector _scattered; LorentzVector _gamma_ion; //residual gamma* + ion system - + LorentzVector _parent_in_elFrame; - ROOT::Math::RotationY _rotateToZaxis; + ROOT::Math::RotationZYX _sc_rotateToZaxis; LorentzVector _cachedParent; // TH1D *hist, *histW; @@ -68,12 +68,13 @@ namespace elSpectro{ }; - inline void ScatteredElectron_xy::RotateZaxisToCMDirection(const LorentzVector& parent){ + inline void ScatteredElectron_xy::RotateZaxisToCMDirection(const LorentzVector& parent, LorentzVector& child){ if(_cachedParent!=parent){ //SetAngle is expensive (sin,cos calls) only call if necessary _cachedParent = parent; - _rotateToZaxis.SetAngle(_cachedParent.Theta()); + //_rotateToZaxis.SetAngle(_cachedParent.Theta()); + _sc_rotateToZaxis.SetComponents(-_cachedParent.Phi(),-_cachedParent.Theta(),0); } - _scattered=_rotateToZaxis*_scattered; + child=_sc_rotateToZaxis*child; } } From 071da82bcc2930f06818a289ea9a17c6865dc821 Mon Sep 17 00:00:00 2001 From: dglazier Date: Thu, 22 Apr 2021 18:32:54 +0100 Subject: [PATCH 08/21] revert to logy and x for integral calculation, swithc to genetic minimiser to find maximum --- core/DecayModelst.cpp | 69 +++++++++++++++++----- core/DecayModelst.h | 29 +++++----- core/DistVirtPhotFlux_xy.h | 7 ++- core/ElectronScattering.cpp | 112 +++++++++++++++--------------------- 4 files changed, 120 insertions(+), 97 deletions(-) diff --git a/core/DecayModelst.cpp b/core/DecayModelst.cpp index ac61d18..1b08903 100644 --- a/core/DecayModelst.cpp +++ b/core/DecayModelst.cpp @@ -49,7 +49,7 @@ namespace elSpectro{ double maxW = ( *(_prodInfo->_target) + *(_prodInfo->_ebeam) ).M(); - _max = FindMaxOfIntensity()*1.05; //add 20% for Q2,meson mass effects etc. + _max = FindMaxOfIntensity()*1.05; //add 5% for Q2,meson mass effects etc. std::cout<<"DecayModelst::PostInit max value "<<_max<<" "<<_meson<<" "<<_meson->Pdg()<<" "<<_sdmeMeson< kine::t0(_W,M1,M2,M3,M4) ) return 0.; - if( x[1]< kine::tmax(_W,M1,M2,M3,M4) ) return 0.; - _t=x[1]; - double val = DifferentialXSect()* (kine::t0(_W,M1,M2,M3,M4)-kine::tmax(_W,M1,M2,M3,M4)); + // if( x[1] > kine::t0(_W,M1,M2,M3,M4) ) return 0.; + // if( x[1]< kine::tmax(_W,M1,M2,M3,M4) ) return 0.; + + auto myt0=kine::t0(_W,M1,M2,M3,M4); + auto mytmax=kine::tmax(_W,M1,M2,M3,M4); + + auto currt=x[1]; + + /*if( TMath::IsNaN(val) ) { + if(TMath::Abs(x[1]-myt0) myt0){ + _t=myt0; + factor*=(1-1*(x[1]-myt0)); //lower result + } + else if( x[1]< mytmax ){ + _t=mytmax; + factor*=(1-1*(mytmax-x[1])); //lower result + } + else _t=x[1]; + double val = DifferentialXSect()* (myt0-mytmax)*factor; + // std::cout<SetMaxFunctionCalls(1000000); // for Minuit/Minuit2 - minimum->SetMaxIterations(10000); // for GSL + minimum->SetMaxIterations(1000); // for GSL minimum->SetTolerance(0.0001); minimum->SetPrintLevel(1); @@ -179,10 +214,11 @@ namespace elSpectro{ ROOT::Math::Functor wrapf(Fmax,2); //variable = W, variable 1 = t - double step[2] = {0.1,0.001}; + double step[2] = {Wrange/500,0.1}; // starting point double variable[2] = { gridW,gridt}; + //double variable[2] = { 10,-0.1}; minimum->SetFunction(wrapf); @@ -191,7 +227,7 @@ namespace elSpectro{ minimum->SetVariable(1,"t",variable[1], step[1]); minimum->SetVariableLimits(0,Wmin,_Wmax); minimum->SetVariableLimits(1,tmax,0); - + // do the minimization minimum->Minimize(); const double *xs = minimum->X(); @@ -217,13 +253,20 @@ namespace elSpectro{ //std::cout << "Minimum Mass Maximum : Probabiltiy Dist at ( W=" << minW << " , t = " << mint << "): "<< -minimum->MinValue() << std::endl; - if(minminValPdgMass()>M3) dynamic_cast(_meson)->TakePdgMass(); } + if(gridMin nb double PhaseSpaceFactor() const noexcept { - /* auto fluxPhaseSpace = p1*_W;//eqn 47.28b https://pdg.lbl.gov/2019/reviews/rpp2019-rev-kinematics.pdf - auto ans = 1./fluxPhaseSpace - * kine::PhaseSpaceFactorDt(_W,p1,_meson->Mass(),_baryon->Mass()) - * kine::PDK(_W,_meson->Mass(),_baryon->Mass())/_W - * PhaseSpaceNorm();//nbarn it - */ //Note above full calculation simplifies to + /* auto fluxPhaseSpace = p1*_W;//eqn 47.28b https://pdg.lbl.gov/2019/reviews/rpp2019-rev-kinematics.pdf + auto ans = 1./fluxPhaseSpace + * kine::PhaseSpaceFactorDt(_W,p1,_meson->Mass(),_baryon->Mass()) + * kine::PDK(_W,_meson->Mass(),_baryon->Mass())/_W + * PhaseSpaceNorm();//nbarn it + */ //Note above full calculation simplifies to //return PhaseSpaceNorm()/_s/kine::PDK2(_W,_photon->M(),_target->M()); //Please note kine::PDK2(_W,_photon->M(),_target->M()) does not give //correct momentum return PhaseSpaceNorm()/_s/PgammaCMsq(); //this would not be the case if the differential was dcosth rather than t } - - constexpr double PhaseSpaceNormCosTh() const {return 1./(2.56819E-6)/32/TMath::Pi();}// Convert from GeV^-2 -> nb double PhaseSpaceFactorCosTh() const noexcept { return PhaseSpaceNormCosTh()* kine::PDK(_W,_meson->Mass(),_baryon->Mass())/_s/PgammaCM(); } + protected: + + virtual double MatrixElementsSquared_L() const {return 0; } + virtual double MatrixElementsSquared_T() const {return 1; } //just real photon by default + + constexpr double PhaseSpaceNorm() const {return 1./(2.56819E-6)/64/TMath::Pi();}// Convert from GeV^-2 -> nb + + constexpr double PhaseSpaceNormCosTh() const {return 1./(2.56819E-6)/32/TMath::Pi();}// Convert from GeV^-2 -> nb + SDME* const GetMesonSDMEs() const {return _sdmeMeson;} SDME* const GetBaryonSDMEs() const {return _sdmeBaryon;} diff --git a/core/DistVirtPhotFlux_xy.h b/core/DistVirtPhotFlux_xy.h index 097f9d3..0293321 100644 --- a/core/DistVirtPhotFlux_xy.h +++ b/core/DistVirtPhotFlux_xy.h @@ -119,7 +119,7 @@ namespace elSpectro{ }; - /* //evaluate photon flux as function of lnx and lny + //evaluate photon flux as function of lnx and lny inline double DistVirtPhotFlux_xy::Eval(const double *x) const{ double lnx=x[0]; @@ -141,7 +141,8 @@ namespace elSpectro{ if(currx(mutableDecayer()); - auto xvar = RooRealVar(Form("xIntegral%lf_%lf",photonFlux->Dist().GetMaxLnX(),photonFlux->Dist().GetMaxLnX()),"xIntegral",TMath::Exp(photonFlux->Dist().GetMaxLnX()),TMath::Exp(photonFlux->Dist().GetMinLnX()),TMath::Exp(photonFlux->Dist().GetMaxLnX()),""); - auto yvar = RooRealVar(Form("yIntegral%lf_%lf",photonFlux->Dist().GetMaxLnY(),photonFlux->Dist().GetMaxLnY()),"yIntegral",TMath::Exp(photonFlux->Dist().GetMaxLnY()),TMath::Exp(photonFlux->Dist().GetMinLnY()),TMath::Exp(photonFlux->Dist().GetMaxLnY()),""); - auto cthvar = RooRealVar("CosThIntegral","CosThIntegral",-0.,-0.99999,0.99999,""); + // auto xvar = RooRealVar(Form("xIntegral%lf_%lf",photonFlux->Dist().GetMaxLnX(),photonFlux->Dist().GetMaxLnX()),"xIntegral",TMath::Exp(photonFlux->Dist().GetMinLnX()),TMath::Exp(photonFlux->Dist().GetMinLnX()),TMath::Exp(photonFlux->Dist().GetMaxLnX()),""); + //auto yvar = RooRealVar(Form("yIntegral%lf_%lf",photonFlux->Dist().GetMaxLnY(),photonFlux->Dist().GetMaxLnY()),"yIntegral",TMath::Exp(photonFlux->Dist().GetMinLnY()),TMath::Exp(photonFlux->Dist().GetMinLnY()),TMath::Exp(photonFlux->Dist().GetMaxLnY()),""); + + auto xvar = RooRealVar(Form("xIntegral%lf_%lf",photonFlux->Dist().GetMaxLnX(),photonFlux->Dist().GetMaxLnX()),"xIntegral",(photonFlux->Dist().GetMinLnX()),(photonFlux->Dist().GetMinLnX()),(photonFlux->Dist().GetMaxLnX()),""); + auto yvar = RooRealVar(Form("yIntegral%lf_%lf",photonFlux->Dist().GetMaxLnY(),photonFlux->Dist().GetMaxLnY()),"yIntegral",(photonFlux->Dist().GetMinLnY()),(photonFlux->Dist().GetMinLnY()),(photonFlux->Dist().GetMaxLnY()),""); + + + + // auto cthvar = RooRealVar("CosThIntegral","CosThIntegral",-0.,-0.99999,0.99999,""); + auto cthvar = RooRealVar("CosThIntegral","CosThIntegral",0.99,-0.99999,0.99999,""); xvar.Print(); yvar.Print(); cthvar.Print(); @@ -330,100 +337,71 @@ namespace elSpectro{ auto Eel=_nuclRestElec.E(); double_t threshW= gStarModel->GetMeson()->PdgMass()+gStarModel->GetBaryon()->PdgMass(); - std::cout<<"Eel "<Dist().Eval(x); if(TMath::IsNaN(val)) return 0.; if(val==0) return 0.; - //std::cout<<"a fXYcosth "<GenerateGivenXandY(P4(),Model()->Products(),x[0],x[1]); + // photonFlux->GenerateGivenXandY(P4(),Model()->Products(),x[0],x[1]); + photonFlux->GenerateGivenXandY(P4(),Model()->Products(),TMath::Exp(x[0]),TMath::Exp(x[1])); //calculate virtual photon Q2WModel->Intensity(); - //std::cout<<"ElectronScattering integrate g*p "<Parent()->P4()<Parent()->P4().M()<Parent()->P4().M();//TMath::Sqrt(escat::M2_pr()+2*Eel*escat::M_pr()*x[1]-escat::Q2_xy( Eel,x[0],x[1])); - //get value of dsigma(s)/dcosth cross section at x,y,costh - val*=gStarModel->dsigma_costh(x[2]); - //std::cout<<"b fXYcosth "<Parent()->P4().M()<<" thesh"<getQ2()<getQ2()<<" "<getW()<get_Q2()<<" "<get_W()<dsigma_costh(x[2]); + val*=dsigma_costh; + //std::cout<<"val "<< val<<" "<dsigma_costh(0.99)<Q2H1Rho(); return val; }; - /* - auto fWcosTh = [this,&photonFlux,&gStarModel,&Q2WModel,Eel](const double *x) - { - auto currx=x[2]; - auto W=x[0]; - auto y = (W*W - escat::M2_pr())/2/Eel/escat::M_pr(); - if(TMath::IsNaN(y)) return 0.; - Double_t xx[2]={currx,(y)}; - auto val = photonFlux->Dist().Eval(xx)*W/(Eel*escat::M_pr()); - //auto val =y; - if(val==0) return 0.; - //auto val = 1.0; - //calculate scatered electron at x and y - photonFlux->GenerateGivenXandY(P4(),Model()->Products(),(currx),y); - //photonFlux->GenerateGivenXandY(P4(),Model()->Products(),currx,y); - //calculate virtual photon - Model()->Intensity(); - //get value of dsigma/ds/dcosth cross section at x,y,costh - val*=gStarModel->dsigma_costhW(x[1],W); - //additional (not real photo) Q2dependence of cross section - if(TMath::IsNaN(val)) return 0.; - if(val<0) return 0.; - // val*=Q2WModel->Q2H1Rho(); - return val; - }; - - auto Wvar = RooRealVar("WIntegral","WIntegral",_Wmin,_Wmin, collision.M(),""); - - - auto wrapPdfW=ROOT::Math::Functor( fWcosTh , 3); - - auto pdfW = RooFunctorPdfBinding("ElScatterIntegral", "ElScatterIntegral", wrapPdfW, RooArgList(Wvar,cthvar,xvar)); - auto roovarsW= RooArgSet(Wvar,cthvar,xvar); - - gBenchmark->Start("RooFitIntegralW"); - - auto RFintegralW=pdfW.getNorm(roovarsW); - gBenchmark->Stop("RooFitIntegralW"); - gBenchmark->Print("RooFitIntegralW"); + double funcvars[3]; + funcvars[0]=0;funcvars[1]=0.1;TMath::Exp(photonFlux->Dist().GetMinLnY()); + for(Int_t i=0;i<100;i++){ + funcvars[0]=funcvars[0]+0.01; + funcvars[2]=-1; + std::cout<<" cosTh -1 "<Print(); - std::cout<<_cacheIntegrals<<" "<Print(); auto roovars= RooArgSet(xvar,yvar,cthvar); - //auto integralRoo=pdf.getNormIntegral(roovars); - //integralRoo->setDirtyInhibit(true); - + gBenchmark->Start("RooFitIntegral"); - auto RFintegral=pdf->getNorm(roovars); - delete pdf; - //auto RFintegral=integralRoo->getValV();//(roovars); + auto RFintegral=pdf.getNorm(roovars); + gBenchmark->Stop("RooFitIntegral"); gBenchmark->Print("RooFitIntegral"); std::cout<<" ElectronScattering::IntegrateCrossSection() "<Dist().Integral()<<" nb"<SetMaxFunctionCalls(1000000); // for Minuit/Minuit2 minimum->SetMaxIterations(1000); // for GSL minimum->SetTolerance(0.0001); - minimum->SetPrintLevel(1); + minimum->SetPrintLevel(0); // create function wrapper for minimizer // a IMultiGenFunction type ROOT::Math::Functor wrapf(Fmax,2); //variable = W, variable 1 = t - double step[2] = {Wrange/500,0.1}; + double step[2] = {Wrange/500,0.001}; // starting point double variable[2] = { gridW,gridt}; @@ -236,7 +241,7 @@ namespace elSpectro{ auto minW= xs[0]; auto mint= xs[1]; - std::cout << "Maximum : Probabiltiy Dist at ( W=" << minW << " , t = " << mint << "): "<< -minimum->MinValue() << std::endl; + std::cout << "Maximum : Probabiltiy Dist at ( W=" << minW << " , t = " << mint << "): "<< -minimum->MinValue() << " note t0 "<(_meson)){ //meson @@ -251,7 +256,7 @@ namespace elSpectro{ minW= xs[0]; mint= xs[1]; - //std::cout << "Minimum Mass Maximum : Probabiltiy Dist at ( W=" << minW << " , t = " << mint << "): "<< -minimum->MinValue() << std::endl; + std::cout << "Minimum Mass Maximum : Probabiltiy Dist at ( W=" << minW << " , t = " << mint << "): "<< -minimum->MinValue()<< " note t0 "< Date: Mon, 26 Apr 2021 14:21:08 +0100 Subject: [PATCH 10/21] add SetNEvents_via_LuminosityTimeFast --- core/Manager.h | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/core/Manager.h b/core/Manager.h index 4e623ec..0caa663 100644 --- a/core/Manager.h +++ b/core/Manager.h @@ -55,12 +55,21 @@ namespace elSpectro{ SetNEvents(n_or_lum); return; } - std::cout<<"SetNEvents_via_LuminosityTime "<IntegrateCrossSection(); _nEventsToGen=n_or_lum*1E-33*beamtime*_integralXSection*Reaction()->BranchingFraction();//1E-33(cm2tonb) std::cout<<"Manager::SetNEvents_via_LuminosityTime , going to generate "<<_nEventsToGen<<" events"<IntegrateCrossSectionFast(); + _nEventsToGen=n_or_lum*1E-33*beamtime*_integralXSection*Reaction()->BranchingFraction();//1E-33(cm2tonb) + std::cout<<"Manager::SetNEvents_via_LuminosityTimeFast , going to generate "<<_nEventsToGen<<" events"< Date: Mon, 26 Apr 2021 14:22:01 +0100 Subject: [PATCH 11/21] tidy and switch photon flux calculations to log x and y in Dist --- core/DecayModelst.cpp | 5 ++--- core/DistVirtPhotFlux_xy.cpp | 8 ++++---- core/ElectronScattering.cpp | 27 +++------------------------ 3 files changed, 9 insertions(+), 31 deletions(-) diff --git a/core/DecayModelst.cpp b/core/DecayModelst.cpp index a8af24c..4f5dfe3 100644 --- a/core/DecayModelst.cpp +++ b/core/DecayModelst.cpp @@ -258,8 +258,8 @@ namespace elSpectro{ std::cout << "Minimum Mass Maximum : Probabiltiy Dist at ( W=" << minW << " , t = " << mint << "): "<< -minimum->MinValue()<< " note t0 "<minVal){ + //std::cout<<"minmin "<<-minminVal<<" < "<<-minVal<(_gStarN->Model()); auto Q2WModel =dynamic_cast(Model()); @@ -345,39 +343,20 @@ namespace elSpectro{ auto val = photonFlux->Dist().Eval(x); if(TMath::IsNaN(val)) return 0.; if(val==0) return 0.; - // std::cout<<"a fXYcosth "<GenerateGivenXandY(P4(),Model()->Products(),x[0],x[1]); photonFlux->GenerateGivenXandY(P4(),Model()->Products(),TMath::Exp(x[0]),TMath::Exp(x[1])); //calculate virtual photon Q2WModel->Intensity(); - //std::cout<<"Q2WModel "<getQ2()<<" "<getW()<get_Q2()<<" "<get_W()<dsigma_costh(x[2]); val*=dsigma_costh; - //std::cout<<"val "<< val<<" "<dsigma_costh(0.99)<Q2H1Rho(); return val; }; - /* - double funcvars[3]; - funcvars[0]=0;funcvars[1]=0.1;TMath::Exp(photonFlux->Dist().GetMinLnY()); - for(Int_t i=0;i<100;i++){ - funcvars[0]=funcvars[0]+0.01; - funcvars[2]=-1; - std::cout<<" cosTh -1 "< Date: Mon, 26 Apr 2021 14:23:53 +0100 Subject: [PATCH 12/21] readme --- CMakeLists.txt | 10 ++++++---- README.md | 22 +++++++++++++--------- 2 files changed, 19 insertions(+), 13 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index a9662e8..28c0699 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,13 +1,15 @@ cmake_minimum_required(VERSION 2.9) -SET(CMAKE_INSTALL_PREFIX "DEFAULT") #do not /usr/local as default - project(elSpectro) #if -DCMAKE_INSTALL_PREFIX=/install/here use that -if(NOT CMAKE_INSTALL_PREFIX STREQUAL "DEFAULT") +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) + +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 STREQUAL "DEFAULT") +ENDIF(NOT CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT) #if no default install dirs, install at source if(NOT DEFINED CMAKE_INSTALL_BINDIR) diff --git a/README.md b/README.md index 962fdfa..b721e40 100644 --- a/README.md +++ b/README.md @@ -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 @@ -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 ../ @@ -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 From 291ca5868a62ce11bdabaa9279f004cab4fd160c Mon Sep 17 00:00:00 2001 From: dglazier Date: Tue, 27 Apr 2021 13:28:39 +0100 Subject: [PATCH 13/21] add Quasifree nucleon generator --- core/CMakeLists.txt | 6 ++-- core/ElSpectroLinkDef.h | 1 + core/QuasiFreeNucleon.cpp | 56 +++++++++++++++++++++++++++++++++++++ core/QuasiFreeNucleon.h | 59 +++++++++++++++++++++++++++++++++++++++ 4 files changed, 120 insertions(+), 2 deletions(-) create mode 100644 core/QuasiFreeNucleon.cpp create mode 100644 core/QuasiFreeNucleon.h diff --git a/core/CMakeLists.txt b/core/CMakeLists.txt index bedfc57..0012353 100644 --- a/core/CMakeLists.txt +++ b/core/CMakeLists.txt @@ -20,6 +20,7 @@ ROOT_GENERATE_DICTIONARY(G__${ELSPECTRO} TwoBodyFlat.h TwoBody_stu.h ScatteredElectron_xy.h + QuasiFreeNucleon.h ProductionProcess.h ElectronScattering.h Distribution.h @@ -32,7 +33,7 @@ ROOT_GENERATE_DICTIONARY(G__${ELSPECTRO} SDME.h PhotonPolarisationVector.h DistVirtPhotFlux_xy.h - CurrentEventInfo.h +# CurrentEventInfo.h ParticleManager.h DecayManager.h MassPhaseSpace.h @@ -60,6 +61,7 @@ add_library(${ELSPECTRO} SHARED TwoBodyFlat.cpp TwoBody_stu.cpp ScatteredElectron_xy.cpp + QuasiFreeNucleon.cpp ProductionProcess.cpp ElectronScattering.cpp Distribution.cpp @@ -70,7 +72,7 @@ add_library(${ELSPECTRO} SHARED DistVirtPhotFlux_xy.cpp SDME.cpp PhotonPolarisationVector.cpp - CurrentEventInfo.cpp +# CurrentEventInfo.cpp ParticleManager.cpp DecayManager.cpp MassPhaseSpace.cpp diff --git a/core/ElSpectroLinkDef.h b/core/ElSpectroLinkDef.h index e9e7aea..b5d9748 100644 --- a/core/ElSpectroLinkDef.h +++ b/core/ElSpectroLinkDef.h @@ -26,6 +26,7 @@ #pragma link C++ class elSpectro::TwoBodyFlat+; #pragma link C++ class elSpectro::TwoBody_stu+; #pragma link C++ class elSpectro::ScatteredElectron_xy+; +#pragma link C++ class elSpectro::QuasiFreeNucleon+; #pragma link C++ class elSpectro::ProductionProcess+; #pragma link C++ class elSpectro::ElectronScattering+; diff --git a/core/QuasiFreeNucleon.cpp b/core/QuasiFreeNucleon.cpp new file mode 100644 index 0000000..40d78b6 --- /dev/null +++ b/core/QuasiFreeNucleon.cpp @@ -0,0 +1,56 @@ +#include "QuasiFreeNucleon.h" +#include + +namespace elSpectro{ + + + //////////////////////////////////////////////////////////////////// + ///Caclulate two body decay from masses and random costh and phi + ///Return a weight that gives phase-space distribution + double QuasiFreeNucleon::Generate(const LorentzVector& parent, const particle_ptrs& products) { + _weight=1; + //Get Fermi momentum sampled from distribution given in input file + auto ptar = _fermiDist->SampleSingle(); + //Get random flat cos theta + auto costhtar = RandomCosTh(); + auto sinthtar = TMath::Sqrt( 1-costhtar*costhtar ); + //Get random flat phi + // done in boost + // phtar = RandomPhi(); + //calculate momentum components + auto pxtar = ptar * sinthtar; + auto pytar = 0; + auto pztar = ptar * costhtar; + + // Force spectator on mass shell + //get its mass from input data and PDG database + auto smass2=products[1]->M2(); + auto Espec = TMath::Sqrt(ptar*ptar + smass2); + //set spectator 4 momentum (opposite p to quasi target) + _spectator.SetXYZT(-pxtar,-pytar,-pztar, Espec); + + //calculate quasi target energy from energy conservation + //This will give offshell mass + auto Etar = parent.M() - Espec; + _nucleon.SetXYZT(pxtar,pytar,pztar, Etar); + + + //Must make sure scattered products are in the same frame as the parent + //theta=0 => moving along parent direction + //i.e. boost vector should only have z component + //Please note this needs checked for correct rotation + //Also add the random phi angle... + // std::cout<<"spectator "<<_spectator<<" "<<_spectator.M()<SetP4(_spectator); + products[0]->SetP4( _nucleon= parent - _spectator ); + //std::cout<<"nucleon boosted"<<_nucleon<<" "<<_nucleon.M()<nucleon + spectator +#pragma once + +#include "DecayVectors.h" +#include "Distribution.h" +#include "DistTF1.h" +#include "LorentzVector.h" +#include //for gRandom +#include //for boosts etc. + +namespace elSpectro{ + + using ROOT::Math::VectorUtil::boost; + + class QuasiFreeNucleon : public DecayVectors { + + public: + + //default deuteron Fermi distribution, can overwrite when create object + QuasiFreeNucleon(Distribution* dist=new DistTF1{TF1("deuteronFermiDist","(x*(0.26**2-0.0456**2)/(x**2+0.0456**2)/(x**2+0.26**2))**2",0.,1)}):_fermiDist{dist}{}; + + virtual ~QuasiFreeNucleon()=default; + QuasiFreeNucleon(const QuasiFreeNucleon& other); //need the virtual destructor...so rule of 5 + QuasiFreeNucleon(QuasiFreeNucleon&&)=default; + QuasiFreeNucleon& operator=(const QuasiFreeNucleon& other); + QuasiFreeNucleon& operator=(QuasiFreeNucleon&& other) = default; + + + double Generate(const LorentzVector& parent, + const particle_ptrs& products) final; + + double RandomCosTh() const noexcept{ + return gRandom->Uniform(-1,1); + } + + + protected : + + + private: + + LorentzVector _nucleon; + LorentzVector _spectator; + + std::unique_ptr _fermiDist; + + ClassDef(elSpectro::QuasiFreeNucleon,1); //class DecayVectors + + + }; + + +} + From 19ec22f7431ace6eb6d54a9c558a3f94b94dbcfb Mon Sep 17 00:00:00 2001 From: dglazier Date: Tue, 27 Apr 2021 17:31:05 +0100 Subject: [PATCH 14/21] introduce colliding particle class to allow sampling of initial state --- core/CMakeLists.txt | 4 ++ core/CollidingParticle.cpp | 30 ++++++++++++ core/CollidingParticle.h | 90 +++++++++++++++++++++++++++++++++++ core/DecayVectors.h | 3 -- core/ElSpectroLinkDef.h | 2 + core/ElectronScattering.cpp | 60 ++--------------------- core/Interface.h | 27 +++++++++-- core/NuclearBreakup.cpp | 15 ++++++ core/NuclearBreakup.h | 43 +++++++++++++++++ core/ParticleManager.cpp | 10 +++- core/ParticleManager.h | 2 +- core/ScatteredElectron_xy.cpp | 6 +-- core/TwoBodyFlat.h | 10 ++-- core/TwoBody_stu.h | 2 +- 14 files changed, 227 insertions(+), 77 deletions(-) create mode 100644 core/CollidingParticle.cpp create mode 100644 core/CollidingParticle.h create mode 100644 core/NuclearBreakup.cpp create mode 100644 core/NuclearBreakup.h diff --git a/core/CMakeLists.txt b/core/CMakeLists.txt index 0012353..c01d957 100644 --- a/core/CMakeLists.txt +++ b/core/CMakeLists.txt @@ -8,8 +8,10 @@ include_directories($ENV{JPACPHOTO}/include/amplitudes) ROOT_GENERATE_DICTIONARY(G__${ELSPECTRO} Particle.h DecayingParticle.h + CollidingParticle.h DecayModel.h PhaseSpaceDecay.h + NuclearBreakup.h VectorSDMEDecay.h DecayModelQ2W.h DecayModelst.h @@ -49,8 +51,10 @@ ROOT_GENERATE_DICTIONARY(G__${ELSPECTRO} add_library(${ELSPECTRO} SHARED Particle.cpp DecayingParticle.cpp + CollidingParticle.cpp DecayModel.cpp PhaseSpaceDecay.cpp + NuclearBreakup.cpp VectorSDMEDecay.cpp DecayModelQ2W.cpp DecayModelst.cpp diff --git a/core/CollidingParticle.cpp b/core/CollidingParticle.cpp new file mode 100644 index 0000000..733ea21 --- /dev/null +++ b/core/CollidingParticle.cpp @@ -0,0 +1,30 @@ +#include "CollidingParticle.h" + +namespace elSpectro{ + + /////////////////////////////////////////////////////////// + //cannot default construct at least need a 4-momentum + CollidingParticle::CollidingParticle(int pdg,LorentzVector lv): + Particle(pdg){ + SetP4(lv); + _interactingParticle=P4ptr(); + } + ///////////////////////////////////////////////////////// + //or a model to generate a 4-momentum + CollidingParticle::CollidingParticle(int pdg,int parentpdg,DecayModel* model,DecayVectors* decayer): + Particle(parentpdg),_model{model},_decayer{decayer} + { + //Find the relevent particle pointer in the model + for(auto& p:_model->Products()){ + if(p->Pdg()==pdg){ + if(_interactingParticle!=nullptr){ + std::cerr<<"CollidingParticle::CollidingParticle, multiple particles with pdg = "<P4ptr(); + + } + } + } + +} diff --git a/core/CollidingParticle.h b/core/CollidingParticle.h new file mode 100644 index 0000000..6b672b2 --- /dev/null +++ b/core/CollidingParticle.h @@ -0,0 +1,90 @@ +////////////////////////////////////////////////////////////// +/// +///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,LorentzVector lv); + //or a model to generate a 4-momentum + CollidingParticle(int pdg,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(){ + Double_t weight=Generate(); + + if(_model )weight*=_model->Intensity(); + + // can add beam divergence etc. here + // ... + return weight; + } + + void PostInit(ReactionInfo* info){ + + } + + LorentzVector* GetInteracting4Vector() const {return _interactingParticle;} + + void SetVertexXYZT(double x,double y,double z,double t){ + _productionVertex.SetXYZT(x,y,z,t); + } + + + void SetAngleTheta(Double_t angle){_dirTheta=angle;} + void SetAnglePhi(Double_t angle){_dirPhi=angle;} + 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 _decayer={nullptr}; //owner + + LorentzVector *_interactingParticle={nullptr}; //not owner; + + LorentzVector _productionVertex; + int _prodVertexID={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 diff --git a/core/DecayVectors.h b/core/DecayVectors.h index d571bfa..5e77c31 100644 --- a/core/DecayVectors.h +++ b/core/DecayVectors.h @@ -38,7 +38,6 @@ namespace elSpectro{ virtual double dsigma() const {return 1;} - void ForIntegrate(bool integ){_forIntegral=integ;} virtual double Probability() const {return 1;} protected: @@ -107,8 +106,6 @@ namespace elSpectro{ protected: - bool _forIntegral=false; - private: LorentzVector _cachedParent; diff --git a/core/ElSpectroLinkDef.h b/core/ElSpectroLinkDef.h index b5d9748..6421e2c 100644 --- a/core/ElSpectroLinkDef.h +++ b/core/ElSpectroLinkDef.h @@ -12,6 +12,7 @@ #pragma link C++ class elSpectro::Particle+; #pragma link C++ class elSpectro::DecayingParticle+; +#pragma link C++ class elSpectro::CollidingParticle+; #pragma link C++ class elSpectro::DecayModel+; #pragma link C++ class elSpectro::DecayModelQ2W+; @@ -20,6 +21,7 @@ #pragma link C++ class elSpectro::GenericModelst+; #pragma link C++ class elSpectro::DecayGammaN_Test+; #pragma link C++ class elSpectro::PhaseSpaceDecay+; +#pragma link C++ class elSpectro::NuclearBreakup+; #pragma link C++ class elSpectro::VectorSDMEDecay+; #pragma link C++ class elSpectro::DecayVectors+; diff --git a/core/ElectronScattering.cpp b/core/ElectronScattering.cpp index 8598f12..fe9862c 100644 --- a/core/ElectronScattering.cpp +++ b/core/ElectronScattering.cpp @@ -15,60 +15,7 @@ namespace elSpectro{ int ElectronScattering::NintegralsElectronScattering=0; - // ElectronScattering::ElectronScattering(double ep,double ionp, - // double eangle,double ionangle, - // DecayModel* model, int ionpdg): - // _pElectron{ep}, - // _pIon{ionp}, - // _angleElectron{eangle}, - // _angleIon{ionangle}, - // _pdgIon{ionpdg}, - // ProductionProcess{model} - // { - // SetBeamCondtion(); - // } - // ////////////////////////////////////////////////////////////////// - // ElectronScattering::ElectronScattering(double ep,double ionp, - // DecayModel* model, int ionpdg): - // _pElectron{ep}, - // _pIon{ionp}, - // _angleElectron{TMath::Pi()}, - // _angleIon{0}, - // _pdgIon{ionpdg}, - // ProductionProcess{model} - // { - // SetBeamCondtion(); - // } - ///////////////////////////////////////////////////////////////////// - /* ElectronScattering::ElectronScattering(double ep,double ionp, DecayVectors* decayer, DecayModel* model, int ionpdg): - _pElectron{ep}, - _pIon{ionp}, - _angleElectron{TMath::Pi()}, - _angleIon{0}, - _pdgIon{ionpdg}, - _beamElec{11}, - _beamNucl{ionpdg}, - ProductionProcess{0,decayer,model} - { - - SetBeamCondtion(); - } - ///////////////////////////////////////////////////////////////////// - ElectronScattering::ElectronScattering(double ep,double ionp, - double eangle,double ionangle, - DecayVectors* decayer, DecayModel* model, int ionpdg): - _pElectron{ep}, - _pIon{ionp}, - _angleElectron{eangle}, - _angleIon{ionangle}, - _pdgIon{ionpdg}, - _beamElec{11}, - _beamNucl{ionpdg}, - ProductionProcess{0,decayer,model} - { - - SetBeamCondtion(); - }*/ + ///////////////////////////////////////////////////////////////////// ElectronScattering::ElectronScattering(double ep,double ionp, DecayModel* model, int ionpdg): _pElectron{ep}, @@ -165,8 +112,6 @@ namespace elSpectro{ auto thresh=Q2WModel->getThreshold(); if(minMass //#include "ParticleManager.h" @@ -39,14 +40,32 @@ namespace elSpectro{ return particles().Take(new Particle{pdg}); } ////////////////////////////////////////////////////////////// - inline Particle* particle(int pdg,DecayModel* const model){ - return particles().Take(new DecayingParticle{pdg,model}); + inline DecayingParticle* particle(int pdg,DecayModel* const model){ + return dynamic_cast(particles().Take(new DecayingParticle{pdg,model})); } ////////////////////////////////////////////////////////////// - inline Particle* particle(int pdg,double mass,DecayModel* const model){ + inline DecayingParticle* particle(int pdg,double mass,DecayModel* const model){ auto p = particles().Take(new DecayingParticle{pdg,model}); p->SetPdgMass(mass); - return p; + return dynamic_cast(p); + } + ////////////////////////////////////////////////////////////// + inline CollidingParticle* initial(int pdg,int parentpdg,DecayModel* model,DecayVectors* decayer){ + return dynamic_cast(particles().Take(new CollidingParticle{pdg,parentpdg,model,decayer})); + } + ////////////////////////////////////////////////////////////// + inline CollidingParticle* initial(int pdg,LorentzVector lv){ + return dynamic_cast(particles().Take(new CollidingParticle{pdg,lv})); + } + ////////////////////////////////////////////////////////////// + inline CollidingParticle* initial(int pdg,Double_t momentum){ + auto mass= particles().GetMassFor(pdg); + LorentzVector lv(0,0,momentum,TMath::Sqrt(momentum*momentum+mass*mass)); + return initial(pdg,lv); + } + ///////////////////////////////////////////////////////////// + inline void add_particle_species(int pdg,double mass){ + particles().AddToPdgTable(pdg,mass); } ////////////////////////////////////////////////////////////// inline void mass_distribution(int pdg,Distribution *dist){ diff --git a/core/NuclearBreakup.cpp b/core/NuclearBreakup.cpp new file mode 100644 index 0000000..7501703 --- /dev/null +++ b/core/NuclearBreakup.cpp @@ -0,0 +1,15 @@ +#include "NuclearBreakup.h" + +namespace elSpectro{ + + /////////////////////////////////////////////////// + ///For nuclei breakup into 2 parts, 1 of which is used + /// as initial particle in collision + NuclearBreakup::NuclearBreakup(int pdg1, int pdg2 ): + DecayModel{{},{pdg1,pdg2}} + { + + + } + +} diff --git a/core/NuclearBreakup.h b/core/NuclearBreakup.h new file mode 100644 index 0000000..10fc347 --- /dev/null +++ b/core/NuclearBreakup.h @@ -0,0 +1,43 @@ +////////////////////////////////////////////////////////////// +/// +///Class: NuclearBreakup +///Description: +/// Break nucleus up into 2 pieces one of which used in collision + +#pragma once + +#include "DecayModel.h" + +namespace elSpectro{ + + + class NuclearBreakup : public DecayModel { + + public: + + NuclearBreakup()=default; + //only declaring default constructor + //so other 5 constructors also defaulted(rule of 5) + //constructor to decay into particles + //NuclearBreakup( particle_ptrs , const std::vector pdgs ); + NuclearBreakup( int pdg1, int pdg2 ); + + // Each model must define its intensity + // Phase space intensity is handled by MassPhaseSpace + double Intensity() const final{ + return 1.; + } + bool RegenerateOnFail() const noexcept final {return false;} + ///void SetParent(DecayingParticle* pa); + //void PostInit(ReactionInfo* info) override; + + private: + + + ClassDefOverride(elSpectro::NuclearBreakup,1); //class NuclearBreakup + + };//class NuclearBreakup + + + +}//namespace elSpectro diff --git a/core/ParticleManager.cpp b/core/ParticleManager.cpp index 36f2a2a..0f676f2 100644 --- a/core/ParticleManager.cpp +++ b/core/ParticleManager.cpp @@ -30,7 +30,8 @@ namespace elSpectro{ 0, kFALSE, 0, 0, "virtual", 9999); - pdgDB->AddParticle("deuteron","deuteron", 1.875612, kTRUE,0, 1, "Baryon", 45); + pdgDB->AddParticle("deuteron","deuteron", 1.875612, kTRUE,0, 1, "Baryon", 45); //Jlab CLAS numbering + pdgDB->AddParticle("deuteron","deuteron", 1.875612, kTRUE,0, 1, "Baryon", 1000010020); //PDG code numbering } Particle* ParticleManager::Take(Particle* p){ @@ -63,5 +64,12 @@ namespace elSpectro{ mass, kFALSE, 0, 0, "virtual", pdg); } + Double_t ParticleManager::GetMassFor(int pdg){ + TDatabasePDG *pdgDB = TDatabasePDG::Instance(); + if(pdgDB->GetParticle(pdg)==nullptr){ + std::cerr<<"ParticleManager::GetMassFor no particle with pdg code "<GetParticle(pdg)->Mass(); + } } diff --git a/core/ParticleManager.h b/core/ParticleManager.h index 348497c..db2336f 100644 --- a/core/ParticleManager.h +++ b/core/ParticleManager.h @@ -57,7 +57,7 @@ namespace elSpectro{ } void AddToPdgTable(int pdg,double mass); - + Double_t GetMassFor(int pdg); /* void RegisterMassSquaredDistribution(int pdg, Distribution* dist){ */ /* _mass2Dist[pdg]=dist_uptr{dist}; */ /* } */ diff --git a/core/ScatteredElectron_xy.cpp b/core/ScatteredElectron_xy.cpp index f01d7be..bfca365 100644 --- a/core/ScatteredElectron_xy.cpp +++ b/core/ScatteredElectron_xy.cpp @@ -109,11 +109,11 @@ namespace elSpectro{ double Ee = escat::E_el(_parent_in_elFrame.P()); //parent in rest frame of ion, momentum=e momentum - double Mion= _parent_in_elFrame.T()-Ee; // energy of parent = Mion + E(e-) + // double Mion= _parent_in_elFrame.T()-Ee; // energy of parent = Mion + E(e-) double Egamma = Ee * yy; - double W = sqrt( Mion*(Mion + 2*Egamma ) - escat::Q2_xy( Ee,xx,yy)); + // double W = sqrt( Mion*(Mion + 2*Egamma ) - escat::Q2_xy( Ee,xx,yy)); double Esc = Ee - Egamma; // if(escat::Q2_xy( Ee,xx,yy)<0.3&&escat::Q2_xy( Ee,xx,yy)>0.01) histy.Fill(yy); @@ -122,7 +122,7 @@ namespace elSpectro{ double costh = escat::CosTh_xy(Ee,xx,yy); costh = costh>1 ? 1 : costh; //protect <=1 - histyCosTh.Fill(costh,yy); + // histyCosTh.Fill(costh,yy); auto sinth=TMath::Sqrt(1-costh*costh); diff --git a/core/TwoBodyFlat.h b/core/TwoBodyFlat.h index 4ca2b30..c761310 100644 --- a/core/TwoBodyFlat.h +++ b/core/TwoBodyFlat.h @@ -30,21 +30,17 @@ namespace elSpectro{ double Generate(const LorentzVector& parent, const particle_ptrs& products) final; - virtual double MyRandomCosTh() const noexcept{ return gRandom->Uniform(-1,1); } virtual double RandomCosTh() const noexcept{ - if(_forIntegral==true) - //integration requires flat distribution - return gRandom->Uniform(-1,1); - return MyRandomCosTh(); + return gRandom->Uniform(-1,1); } + double Probability() const{return 1./4/TMath::Pi();} protected : double W() const noexcept{return _W;} - // void RotateZaxisToCMDirection(const LorentzVector& parent); - + private: LorentzVector _a; diff --git a/core/TwoBody_stu.h b/core/TwoBody_stu.h index f8f8303..5d4ae96 100644 --- a/core/TwoBody_stu.h +++ b/core/TwoBody_stu.h @@ -19,7 +19,7 @@ namespace elSpectro{ TwoBody_stu(double s,double t,double t_slope,double u=0,double u_slope=0); - double MyRandomCosTh() const noexcept final{ + double RandomCosTh() const noexcept final{ _weight=1; auto randChannel = gRandom->Uniform(); // std::cout<<" TwoBody_stu RandomCosTh() "< Date: Wed, 5 May 2021 10:14:45 +0100 Subject: [PATCH 15/21] Finish first implemenation of quasifree nucleons --- core/CollidingParticle.cpp | 76 ++++++++++++- core/CollidingParticle.h | 41 +++---- core/DecayModelQ2W.cpp | 47 ++++---- core/DecayModelQ2W.h | 3 +- core/DecayModelst.cpp | 19 ++-- core/DecayVectors.h | 14 +-- core/DecayingParticle.cpp | 7 +- core/DistTH1.h | 1 + core/DistVirtPhotFlux_xy.cpp | 20 ++-- core/DistVirtPhotFlux_xy.h | 7 +- core/ElectronScattering.cpp | 196 +++++++++++++--------------------- core/ElectronScattering.h | 8 +- core/Interface.h | 25 +++-- core/LundWriter.cpp | 10 +- core/LundWriter.h | 3 +- core/Particle.h | 1 + core/ParticleManager.cpp | 9 +- core/ParticleManager.h | 24 ++++- core/ProductionProcess.cpp | 57 +++++++++- core/ProductionProcess.h | 11 +- core/QuasiFreeNucleon.cpp | 87 ++++++++++++++- core/QuasiFreeNucleon.h | 6 +- core/ReactionInfo.h | 5 +- core/ScatteredElectron_xy.cpp | 9 +- core/ScatteredElectron_xy.h | 9 +- 25 files changed, 465 insertions(+), 230 deletions(-) diff --git a/core/CollidingParticle.cpp b/core/CollidingParticle.cpp index 733ea21..778e409 100644 --- a/core/CollidingParticle.cpp +++ b/core/CollidingParticle.cpp @@ -1,30 +1,96 @@ #include "CollidingParticle.h" +#include "FunctionsForGenvector.h" +#include "Manager.h" namespace elSpectro{ /////////////////////////////////////////////////////////// //cannot default construct at least need a 4-momentum - CollidingParticle::CollidingParticle(int pdg,LorentzVector lv): + 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,int parentpdg,DecayModel* model,DecayVectors* decayer): + CollidingParticle::CollidingParticle(int pdg,Double_t momentum,int parentpdg,DecayModel* model,DecayVectors* decayer): Particle(parentpdg),_model{model},_decayer{decayer} { - //Find the relevent particle pointer in the model - for(auto& p:_model->Products()){ + //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 = "<P4ptr(); + _model->SwapProducts(0,position); //make sure interacting particle is first in products vector + } + //Not a stable final state particle! + Manager::Instance().Particles().RemoveStable(p); + } + 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); + } } diff --git a/core/CollidingParticle.h b/core/CollidingParticle.h index 6b672b2..10ca03d 100644 --- a/core/CollidingParticle.h +++ b/core/CollidingParticle.h @@ -21,14 +21,14 @@ namespace elSpectro{ CollidingParticle()=delete; //cannot default construct at least need a 4-momentum - CollidingParticle(int pdg,LorentzVector lv); + CollidingParticle(int pdg,Double_t momentum); //or a model to generate a 4-momentum - CollidingParticle(int pdg,int parentpdg,DecayModel* model,DecayVectors* decayer); + 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; + if(_decayer.get()==nullptr) return 1.0; return _decayer->Generate(P4(),_model->Products()); } @@ -36,29 +36,31 @@ namespace elSpectro{ void SetDecayer(DecayVectors* decayer){_decayer.reset(decayer);} Double_t GenerateComponents(){ - Double_t weight=Generate(); + if(_model==nullptr )return 1.0; - if(_model )weight*=_model->Intensity(); - + 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){ - - } + void PostInit(ReactionInfo* info) ; - LorentzVector* GetInteracting4Vector() const {return _interactingParticle;} + const LorentzVector* GetInteracting4Vector() const {return _interactingParticle;} + Int_t GetInteractingPdg()const {return _interactingPdg;} - void SetVertexXYZT(double x,double y,double z,double t){ + /*void SetVertexXYZT(double x,double y,double z,double t){ _productionVertex.SetXYZT(x,y,z,t); - } + }*/ - - void SetAngleTheta(Double_t angle){_dirTheta=angle;} - void SetAnglePhi(Double_t angle){_dirPhi=angle;} - void SetHorSize(Double_t s){_sizeHor=s;} + + 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;} @@ -72,9 +74,10 @@ namespace elSpectro{ LorentzVector *_interactingParticle={nullptr}; //not owner; - LorentzVector _productionVertex; - int _prodVertexID={0}; - + 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}; diff --git a/core/DecayModelQ2W.cpp b/core/DecayModelQ2W.cpp index 17a4161..5b44444 100644 --- a/core/DecayModelQ2W.cpp +++ b/core/DecayModelQ2W.cpp @@ -69,7 +69,11 @@ namespace elSpectro{ _prodInfo->_photoN=_gstarNuc->P4ptr(); _prodInfo->_photon=&_gamma; _prodInfo->_photonPol=&_photonPol; - + if(_prodInfo->_Wmax==0){ + //For backward compatability, should be done with + //colliding particles now and set in ProductionProcess::PostInit(). + _prodInfo->_Wmax=( *(_prodInfo->_target) + *(_prodInfo->_ebeam) ).M(); + } auto gNprods=dynamic_cast(_gstarNuc)->Model()->Products(); // std::cout<<"DecayModelQ2W::PostInit "<<_prodInfo<<" "<Pdg()<<" "<Pdg()<_meson=gNprods[0]->P4ptr(); - std::cout<<"DecayModelQ2W::PostInit "<P4().E() < "<<_threshold<<" "<P4().E() <_sWeight=weight; //might be used in s and t - //now add Q2 depedence to get weighted here - - //Q2 dependence of phase space needed to effectively multiply st max value - //Here the Q2 dependence would multiply weight value and max ,thus cancelling - //i.e do not do weight*=PhaseSpaceFactorToQ2eq0(W,p4tar.M() ); - //_prodInfo->_sWeight*=PhaseSpaceFactorToQ2eq0(W,p4tar.M() ); - - + if(weight>1){ auto cmBoost=_gstarNuc->P4().BoostToCM(); auto p1cm=boost(_gamma,cmBoost); std::cout<<" Q2 DEPENDENECE "<_sWeight*=PhaseSpaceFactorToQ2eq0(W,p4tar.M() ); + //Q2 dependence of cross section weight*=Q2H1Rho(); - + return weight; - + } void DecayModelQ2W::FindExcitationSpectra(){ @@ -154,8 +158,9 @@ namespace elSpectro{ //for all values of meson mass, to do this take maximum from values at //meson threshold and pdg mass values //note the point is phase space and therefore xsection changes with mass - double maxW = ( *(_prodInfo->_target) + *(_prodInfo->_ebeam) ).M(); - + //double maxW = ( *(_prodInfo->_target) + *(_prodInfo->_ebeam) ).M(); + double maxW = _prodInfo->_Wmax; + std::cout<<"DecayModelQ2W::PostInit generating max cross section @W, may take some time... "<(_gstarNuc)->Model()->Products(); @@ -163,8 +168,8 @@ namespace elSpectro{ auto meson=gNprods[0]; DecayModelst* mesonBaryon = nullptr; - TH1D histlow("Wdistlow","Wdistlow",200,_threshold,maxW); - TH1D histpeak("Wdisthigh","Wdisthigh",200,_threshold,maxW); + TH1D histlow("Wdistlow","Wdistlow",400,_threshold,maxW); + TH1D histpeak("Wdisthigh","Wdisthigh",400,_threshold,maxW); double minMesonMass=-1; if( ( mesonBaryon=dynamic_cast(GetGammaN()->Model())) != nullptr){ //check for low mass meson limits @@ -204,7 +209,7 @@ namespace elSpectro{ } - TH1D HistFromLargestBinContents(const TH1D& h1,const TH1D& h2){ + TH1D HistFromLargestBinContents(const TH1D& h1,const TH1D& h2){ auto hist= TH1D{h1}; auto maxVal= h1.GetMaximum(); double max_so_far=0.; @@ -218,7 +223,11 @@ namespace elSpectro{ //We want envelope to contain this so once we get to the max //just stay there. Flux is low at high W so no big effect //on efficiency from this - if(valMass();} + //double getW() const noexcept{ return GetGammaN()->Mass();} + double getW() const noexcept{ return _gstarNuc->Mass();} double getThreshold() const noexcept{return _threshold;} void setThreshold(double val) noexcept{ if(val<_threshold) return; diff --git a/core/DecayModelst.cpp b/core/DecayModelst.cpp index 4f5dfe3..b6dafe1 100644 --- a/core/DecayModelst.cpp +++ b/core/DecayModelst.cpp @@ -112,12 +112,12 @@ namespace elSpectro{ auto M2 = _target->M(); auto M3 = _meson->Mass(); //should be pdg value here auto M4 = _baryon->Mass(); - // auto Wmin = M3+M4; auto Wmin = Parent()->MinimumMassPossible(); - - _Wmax = ( *(_prodInfo->_target) + *(_prodInfo->_ebeam) ).M(); - - std::cout<<" DecayModelst::FindMaxOfIntensity() "<_target) + *(_prodInfo->_ebeam) ).M(); + _Wmax = _prodInfo->_Wmax; + + std::cout<<" DecayModelst::FindMaxOfIntensity() Wmin = "<SetVariable(1,"t",variable[1], step[1]); - minimum->SetVariableLimits(0,Wmin,_Wmax); - minimum->SetVariableLimits(1,tmax,0); + // minimum->SetVariableLimits(0,Wmin,_Wmax); + // minimum->SetVariableLimits(1,tmax,0); // do the minimization minimum->Minimize(); @@ -258,8 +258,8 @@ namespace elSpectro{ std::cout << "Minimum Mass Maximum : Probabiltiy Dist at ( W=" << minW << " , t = " << mint << "): "<< -minimum->MinValue()<< " note t0 "<minVal){ - //std::cout<<"minmin "<<-minminVal<<" < "<<-minVal<Uniform(-TMath::Pi(),TMath::Pi()); } - + */ + protected: + mutable double _weight={1}; + virtual double RandomPhi() const noexcept { return gRandom->Uniform(-TMath::Pi(),TMath::Pi()); } + + private: LorentzVector _cachedParent; diff --git a/core/DecayingParticle.cpp b/core/DecayingParticle.cpp index fa6e902..c12c23d 100644 --- a/core/DecayingParticle.cpp +++ b/core/DecayingParticle.cpp @@ -37,7 +37,8 @@ namespace elSpectro{ //decay vertex position auto& products=_decay->Products(); - if(IsDecay()==DecayType::Detached||IsDecay()==DecayType::Production){ + // if(IsDecay()==DecayType::Detached||IsDecay()==DecayType::Production){ + if(IsDecay()==DecayType::Detached){ //create a new detached vertex _decayVertexID=Manager::Instance().AddVertex(&_decayVertex); for(auto* prod: products){ @@ -67,7 +68,7 @@ namespace elSpectro{ double _maxWeight=1; - // std::cout<<"DecayingParticle::GenerateProducts "<Products().size()<<" "<<" "<<_decay->Products()[0]->Pdg()<<" "<<_decay->Products()[1]->Pdg()<Products().size()<<" "<<" "<<_decay->Products()[0]->Pdg()<<" "<<_decay->Products()[1]->Pdg()<HasAngularDistribution()==false)samplingWeight=1; //Model has no angular distribution //samplingWeight ==0 => not physical (below threshold) @@ -83,6 +85,7 @@ namespace elSpectro{ //evaluate the model intensity for the product vectors double weight = 1; if(Model()!=nullptr) weight = Model()->Intensity(); + //std::cout<<"DecayingParticle::GenerateProducts "<GetName()<<" "<Products()[0]->Mass()<<" "<Products()[1]->Mass()<<" difference in weights "<(&_th1))->Interpolate(valX);} const TH1& GetTH1() const noexcept {return _th1;} + void Draw(const TString& opt) { _th1.Draw(opt);} private: //no one should use default constructor diff --git a/core/DistVirtPhotFlux_xy.cpp b/core/DistVirtPhotFlux_xy.cpp index 9a98f72..0a84458 100644 --- a/core/DistVirtPhotFlux_xy.cpp +++ b/core/DistVirtPhotFlux_xy.cpp @@ -49,7 +49,7 @@ namespace elSpectro{ double ymin = rmin/(2*_mTar*_ebeam); if(ymin==0) ymin=1E-50; //now find the y limits with this threshold - // std::cout<<"ymin "<ymin) @@ -64,8 +64,9 @@ namespace elSpectro{ _lnymax=TMath::Log(1.); double xmin =XMin(ymin); - - //random search for the maximum value + + + //random search for the maximum value _max_val=0; for(int i=0;i<1E5;i++){ @@ -84,7 +85,10 @@ namespace elSpectro{ if(avail_xmin<_maxPossiblexRange) _maxPossiblexRange=avail_xmin; } - + + std::cout<<"DistVirtPhotFlux_xy::SetWThreshold new Wmin "<GetInteractingPdg()} + { + //For convenience keep our own + //pointers to electron and target + _electronptr =electron; + _targetptr = target; + + SetNominalBeamCondtion(); + } + + ///////////////////////////////////////////////////////////////////// + void ElectronScattering::SetNominalBeamCondtion(){ + + + _beamElec.SetP4(*(_electronptr->GetInteracting4Vector())); + + _beamNucl.SetP4(*(_targetptr->GetInteracting4Vector())); + _massIon=_beamNucl.PdgMass(); + + /* //not sure why this was there + //For decaying + SetXYZT(_beamElec.P4().X(),_beamElec.P4().Y(), + _beamElec.P4().Z(),_beamElec.P4().T()); + */ + + //For nominal beam conditions set nucleon rest frame vectors + //in case info is needed in PostInit stage + //Boost into ion rest frame + auto prBoost=_beamNucl.P4().BoostToCM(); + _nuclRestNucl=LorentzVector(0,0,0,_beamNucl.Mass()); + _nuclRestElec= boost(_beamElec.P4(),prBoost); + //set inital lab particles + //this can be written to output file + std::cout<<"ElectronScattering vertex "<<_electronptr->VertexPosition()<Dist().SetWThreshold(minMass); _Wmin=minMass; + if(dynamic_cast(&(decayer->Dist())))_Wmin=dynamic_cast((&decayer->Dist()))->GetWMin();//can be effected by angle limits etc } - + std::cout<<"ElectronScattering::InitGen() final minimum W "<<_Wmin<SetMinMass(minMass); + _gStarN->SetMinMass(_Wmin); if(auto Q2WModel=dynamic_cast(Model())){ - Q2WModel->setThreshold(minMass); - } + Q2WModel->setThreshold(_Wmin); + } generator().SetModelForMassPhaseSpace(_gStarN->Model()); } - DecayingParticle::PostInit(dynamic_cast(&_reactionInfo)); + ProductionProcess::PostInit(dynamic_cast(&_reactionInfo)); //now scattered electron should be set //_scattered= _reactionInfo._scattered; @@ -240,7 +291,16 @@ namespace elSpectro{ LorentzVector ElectronScattering::MakeCollision(){ - + //Generate collision 4-momentum + if(_electronptr!=nullptr){ + _electronptr->GenerateComponents(); + _beamElec.SetP4(*(_electronptr->GetInteracting4Vector())); + } + if(_targetptr!=nullptr){ + _targetptr->GenerateComponents(); + _beamNucl.SetP4(*(_targetptr->GetInteracting4Vector())); + } + //First, Eventually want to sample from beam divergence distributions LorentzVector collision = _beamElec.P4() + _beamNucl.P4(); //Boost into ion rest frame @@ -250,6 +310,7 @@ namespace elSpectro{ _nuclRestElec= boost(_beamElec.P4(),prBoost); //set decay parent for e -> e'g* SetXYZT(collision.X(),collision.Y(),collision.Z(),collision.T()); + // std::cout<<"ElectronScattering::MakeCollision() "< e'g* - SetXYZT(collision.X(),collision.Y(),collision.Z(),collision.T()); - */ auto collision=MakeCollision(); //proceed through decay chain while(DecayingParticle::GenerateProducts()!=DecayStatus::Decayed){ - //std::cout<<"not decayed "<<_nsamples<Start("LoopIntegral"); - -// //auto eleDist=static_cast(&decayer->Dist()); -// // static_cast(& _gStarN->mutableDecayer()->Dist())->ForIntegrate(true); -// for(long iW=0;iW<_NIntegrateW;++iW){ -// LorentzVector collision = _beamElec.P4() + _beamNucl.P4(); - -// //Boost into ion rest frame -// auto prBoost=_beamNucl.P4().BoostToCM(); -// collision=boost(collision,prBoost); -// _nuclRestNucl=LorentzVector(0,0,0,_beamNucl.Mass()); -// _nuclRestElec= boost(_beamElec.P4(),prBoost); -// //set decay parent for e -> e'g* -// SetXYZT(collision.X(),collision.Y(),collision.Z(),collision.T()); - -// Decay(); //sample e' and calculate dsigma - -// auto pass=Model()->Intensity();//calculate photon -// if(pass==0) continue; //intensity 0 => does not contribute - -// double intet=0; -// Npass++; -// //dynamic_cast(_gStarN->Model())->kin_tFromWCosTh() -// _gStarN->Decay(); //sample t -// _gStarN->Model()->Intensity(); //calculate dsigma - - -// // _gStarN->Model()->DifferentialXSect(); - -// //correct for electron sampling distribution /Decayer()->Probabilty() -// //and differential cross section sampling /_gStarN->Decayer()->Probability(); - -// integral+=dsigma() -// /Decayer()->Probability() -// /_gStarN->Decayer()->Probability(); - -// // std::cout<<_gStarN->Model()->GetName()<<" "<Probability()<<" "<<_gStarN->Decayer()->Probability()<Probability()<<" "<<_gStarN->Decayer()->Probability()<Decay(); //random t/cosTh -// _gStarN->Model()->Intensity(); - - - -// //histt.Fill(static_cast(_gStarN->Model())->get_t()); - -// intet+=dsigma(); -// }*/ -// //integral+=intet; -// //TH1D histCross("XSection","XSection",1,_gStarN->Mass()-0.001,_gStarN->Mass()+0.001); -// //double Wmid=gRandom->Uniform(5.5,6.5); -// //TH1D histCross("XSection","XSection",1,Wmid-0.001,Wmid+0.001); -// // static_cast(_gStarN->Model())->HistIntegratedXSection(histCross); -// //std::cout<<"Integral for t "<<_gStarN->Mass()<<" "<M2()<<"phase Q2 "<(Model())->PhaseSpaceFactorToQ2eq0(_gStarN->Mass(),_beamNucl.Mass())<(&decayer->Dist())->ForIntegrate(false); -// double Wmax = (_nuclRestNucl + _nuclRestElec).M(); -// std::cout<<"Integral "<< _Wmin<<" - "<Stop("LoopFitIntegral"); -// gBenchmark->Print("LoopFitIntegral"); - -// //exit(0); -// /*TH1D histCross("XSection","XSection",10,Wmin,Wmax); -// std::cout<< "hist integral for "<<_gStarN->Model()->GetName()<(_gStarN->Model())->HistIntegratedXSection(histCross); -// TFile * crossFile=new TFile("XSection.root","recreate"); -// histCross.Write(); -// delete crossFile; - -// exit(0); -// */ -// */ diff --git a/core/ElectronScattering.h b/core/ElectronScattering.h index d26aac0..b11360a 100644 --- a/core/ElectronScattering.h +++ b/core/ElectronScattering.h @@ -39,6 +39,7 @@ namespace elSpectro{ double anglee,double anglep,DecayModel* model=new PhaseSpaceDecay{{},{11,-2211}},int ionpdg=2212); ElectronScattering(double ep,double ionp,DecayModel* model=new PhaseSpaceDecay{{},{11,-2211}},int ionpdg=2212); + ElectronScattering(CollidingParticle *electron,CollidingParticle* target, DecayModel* model); DecayStatus GenerateProducts( ) override; @@ -83,7 +84,8 @@ namespace elSpectro{ void SetBeamCondtion(); - + void SetNominalBeamCondtion(); + Particle _beamElec; Particle _beamNucl; LorentzVector _nuclRestElec; @@ -116,7 +118,9 @@ namespace elSpectro{ short _cacheIntegrals={0}; DecayingParticle* _gStarN={nullptr}; - + CollidingParticle* _electronptr={nullptr}; + CollidingParticle* _targetptr={nullptr}; + ClassDefOverride(elSpectro::ElectronScattering,1); //class ElectronScattering }; diff --git a/core/Interface.h b/core/Interface.h index 1f42337..65a6b8d 100644 --- a/core/Interface.h +++ b/core/Interface.h @@ -50,20 +50,14 @@ namespace elSpectro{ return dynamic_cast(p); } ////////////////////////////////////////////////////////////// - inline CollidingParticle* initial(int pdg,int parentpdg,DecayModel* model,DecayVectors* decayer){ - return dynamic_cast(particles().Take(new CollidingParticle{pdg,parentpdg,model,decayer})); - } - ////////////////////////////////////////////////////////////// - inline CollidingParticle* initial(int pdg,LorentzVector lv){ - return dynamic_cast(particles().Take(new CollidingParticle{pdg,lv})); + inline CollidingParticle* initial(int pdg,Double_t momentum,int parentpdg,DecayModel* model,DecayVectors* decayer){ + return dynamic_cast(particles().Take(new CollidingParticle{pdg,momentum,parentpdg,model,decayer})); } ////////////////////////////////////////////////////////////// inline CollidingParticle* initial(int pdg,Double_t momentum){ - auto mass= particles().GetMassFor(pdg); - LorentzVector lv(0,0,momentum,TMath::Sqrt(momentum*momentum+mass*mass)); - return initial(pdg,lv); + return dynamic_cast(particles().Take(new CollidingParticle{pdg,momentum})); } - ///////////////////////////////////////////////////////////// +///////////////////////////////////////////////////////////// inline void add_particle_species(int pdg,double mass){ particles().AddToPdgTable(pdg,mass); } @@ -124,5 +118,16 @@ namespace elSpectro{ return generator().Reaction(); } + inline ProductionProcess* mesonex(CollidingParticle* el,CollidingParticle* pr,DecayModelQ2W *totalXsec=nullptr){ + + if(totalXsec!=nullptr){ + model(totalXsec); //register ownership of model with manager + generator().Reaction(new ElectronScattering(el,pr, totalXsec )); + } + else + generator().Reaction(new ElectronScattering(el,pr,totalXsec)); + + return generator().Reaction(); + } }//namespace elSpectro diff --git a/core/LundWriter.cpp b/core/LundWriter.cpp index ba8aa72..d14f93a 100644 --- a/core/LundWriter.cpp +++ b/core/LundWriter.cpp @@ -25,7 +25,6 @@ namespace elSpectro{ void LundWriter::Init(){ Writer::Init(); - //need to find e- and baryon if(TDatabasePDG::Instance()->GetParticle(_initialParticles[0]->Pdg())->ParticleClass()==TString("Baryon") ){ _inTarget=_initialParticles[0]; @@ -58,19 +57,20 @@ namespace elSpectro{ ///////////////////////////////////////////////////////////// //write all the info required for this event void LundWriter::FillAnEvent(){ - - + ////fill _stream StreamEventInfo(); _id=1;//reset particle ID counter - + //final particles int final_status=1; for(const auto* p:_finalParticles){ + StreamParticle(p,final_status); - } + + } _nEvent++; diff --git a/core/LundWriter.h b/core/LundWriter.h index d6761fe..b00889c 100644 --- a/core/LundWriter.h +++ b/core/LundWriter.h @@ -45,7 +45,7 @@ namespace elSpectro{ // # of Particles, # of Target Nucleons, # of Target Protons, // Pol. of Target, Pol. of Electron, // BeamType, BeamEnergy,Target ID, ProcessID, Weight - + _stream<< "\t "<<_finalParticles.size()<<" "<<1<<" "<<1 <<" "<<0.<<" "<<0. <<" "<<_beamPdg<<" "<<_inBeam->P4().E()<<" "<<_targetPdg<<" "<< _inTarget->P4().E() <<" "<<0.<<"\n"; @@ -53,7 +53,6 @@ namespace elSpectro{ void StreamParticle(const Particle* p,int status){ auto p4=p->P4(); auto ver=p->VertexPosition(); - //second entry 0. == lifetime Could add to Particle.h _stream<<_id++<<" "<<0.<<" "<Pdg()<<" "<<0<<" "<<0<<" " diff --git a/core/Particle.h b/core/Particle.h index 2611a83..6ee8fb0 100644 --- a/core/Particle.h +++ b/core/Particle.h @@ -95,6 +95,7 @@ namespace elSpectro{ void SetVertex(int vertexID,const LorentzVector* v){ _vertexID=vertexID; _vertex=v; + std::cout<<"PArticle SetVertex "<<_pdg<<" "<<_vertexID<<" "<SetMassDist(_massDist.at(pdg).get()); auto dp=dynamic_cast( p ); + auto cp=dynamic_cast( p ); - //distinguish between stable and unstable particles + //distinguish between colliding, stable and unstable particles if(dp!=nullptr){ _unstables.push_back(dp); } + else if(cp!=nullptr){ + _initials.push_back(cp); + } else{ _stables.push_back(p); } @@ -60,9 +64,10 @@ namespace elSpectro{ void ParticleManager::AddToPdgTable(int pdg,double mass){ TDatabasePDG *pdgDB = TDatabasePDG::Instance(); + //Note make Baryons in case need to write as beam (see LundWriter) pdgDB->AddParticle(Form("Particle%d",pdg),"resonance", mass, kFALSE, - 0, 0, "virtual", pdg); + 0, 0, "Baryon", pdg); } Double_t ParticleManager::GetMassFor(int pdg){ TDatabasePDG *pdgDB = TDatabasePDG::Instance(); diff --git a/core/ParticleManager.h b/core/ParticleManager.h index db2336f..082293a 100644 --- a/core/ParticleManager.h +++ b/core/ParticleManager.h @@ -10,6 +10,7 @@ #pragma once #include "DecayingParticle.h" +#include "CollidingParticle.h" #include "Distribution.h" #include #include @@ -23,6 +24,8 @@ namespace elSpectro{ using particle_constptrs= std::vector ; using decaying_ptrs= std::vector ; using decaying_constptrs= std::vector ; + using colliding_ptrs= std::vector ; + using colliding_constptrs= std::vector ; using ROOT::Math::VectorUtil::boost; @@ -68,8 +71,25 @@ namespace elSpectro{ std::for_each(_stables.begin(),_stables.end(),[&vboost](Particle* p){p->Boost(vboost);}); } + void MoveStableToLab(Particle* particle){ + //So the particle does not get boosted + //but is written out to final state + _stables.erase(std::remove(_stables.begin(),_stables.end(),particle),_stables.end()); + _stableslab.push_back(particle); + } + void RemoveStable(Particle* particle){ + //So the particle does not get boosted + //or written out to final state + _stables.erase(std::remove(_stables.begin(),_stables.end(),particle),_stables.end()); + } + const decaying_ptrs UnstableParticles()const {return _unstables;} - const particle_ptrs StableParticles()const {return _stables;} + const particle_ptrs StableParticles()const { + particle_ptrs allstables; + allstables.insert(std::end(allstables), std::begin(_stables), std::end(_stables)); + allstables.insert(std::end(allstables), std::begin(_stableslab), std::end(_stableslab)); + return allstables; + } void BoostToFrame(const BetaVector& vboost,const LorentzVector& parent){ @@ -101,7 +121,9 @@ namespace elSpectro{ std::vector _particles; particle_ptrs _stables; //products which are stable (to be detected) + particle_ptrs _stableslab; //products which are stable and in lab frame decaying_ptrs _unstables; //products which decay + colliding_ptrs _initials; //intial interactin particles //distribitions for generating dynamic particle masses std::map _massDist; diff --git a/core/ProductionProcess.cpp b/core/ProductionProcess.cpp index 5c1bbfe..ad066b7 100644 --- a/core/ProductionProcess.cpp +++ b/core/ProductionProcess.cpp @@ -1,16 +1,69 @@ #include "ProductionProcess.h" +#include "Manager.h" namespace elSpectro{ ////////////////////////////////////////////////////////////////// - ProductionProcess::ProductionProcess(DecayModel* model): + ProductionProcess::ProductionProcess(CollidingParticle* p1,CollidingParticle* p2,DecayModel* model): + DecayingParticle{0,nullptr,model}, //will set decayer in Init of derived class + _in1{p1}, + _in2{p2} + { + Init(); + } + + ////////////////////////////////////////////////////////////////// + ProductionProcess::ProductionProcess(DecayModel* model): DecayingParticle{model} { - + Init(); } + ////////////////////////////////////////////////////////////////// ProductionProcess::ProductionProcess(int pdg,DecayVectors* decayer, DecayModel* model): DecayingParticle{pdg,decayer,model}{ + Init(); + + } + void ProductionProcess::Init() { + auto decayVertexID=Manager::Instance().AddVertex(&DecayVertexPosition()); + //Set my vertex position to my "decay vertex" + SetVertex(decayVertexID,&DecayVertexPosition()); + //decay vertex position for colliding + if(_in1){ + _in1->SetVertex(VertexID(),VertexPosition()); + } + if(_in2){ + _in2->SetVertex(VertexID(),VertexPosition()); + } + + + //decay vertex position for products + auto& products=Model()->Products(); + for(auto* prod: products){ + std::cout<<" ProductionProcess::Init() "<Pdg()<<" "<SetVertex(VertexID(),VertexPosition()); + } + } + void ProductionProcess::PostInit(ReactionInfo* info) { + + + if(_in1){ + _in1->PostInit(info); + } + if(_in2){ + _in2->PostInit(info); + } + if(_in1&&_in2){ + //Note for now take maximum from sum of both incident particles + //this ignores that one might be a quasifree nucleon + //but kinematically at high intiial nucleon momentum + //we can approach this maximum.... + info->_Wmax=(_in1->P4()+_in2->P4()).M(); + std::cout<<"ProductionProcess::PostInit maximum W = "<< info->_Wmax <M2(); + auto smass2=products[1]->M2(); //Note interacting particle is [0] auto Espec = TMath::Sqrt(ptar*ptar + smass2); //set spectator 4 momentum (opposite p to quasi target) _spectator.SetXYZT(-pxtar,-pytar,-pztar, Espec); @@ -40,17 +40,96 @@ namespace elSpectro{ //i.e. boost vector should only have z component //Please note this needs checked for correct rotation //Also add the random phi angle... - // std::cout<<"spectator "<<_spectator<<" "<<_spectator.M()<SetP4(_spectator); products[0]->SetP4( _nucleon= parent - _spectator ); - //std::cout<<"nucleon boosted"<<_nucleon<<" "<<_nucleon.M()< //for gRandom #include //for boosts etc. @@ -53,7 +54,10 @@ namespace elSpectro{ }; - + + namespace QuasiFree{ + DistTH1* CDBonnMomentum(const Int_t Nbins=1000,const Double_t pmin=0,const Double_t pmax=1); + } } diff --git a/core/ReactionInfo.h b/core/ReactionInfo.h index 1672062..55b696b 100644 --- a/core/ReactionInfo.h +++ b/core/ReactionInfo.h @@ -18,7 +18,9 @@ namespace elSpectro{ virtual ~ReactionInfo()=default; - + double _Wmax=0; + double _Wmin=0; + }; class ReactionPhotoProd : public ReactionInfo { @@ -36,7 +38,6 @@ namespace elSpectro{ PhotonPolarisationVector* _photonPol={nullptr}; mutable double _sWeight = {1}; //s=W^2 excitation function weight - }; class ReactionElectroProd : public ReactionPhotoProd { diff --git a/core/ScatteredElectron_xy.cpp b/core/ScatteredElectron_xy.cpp index bfca365..e145938 100644 --- a/core/ScatteredElectron_xy.cpp +++ b/core/ScatteredElectron_xy.cpp @@ -19,7 +19,7 @@ namespace elSpectro{ ///Sample a and y which are bound [0,1] so independent of eenrgy scale ///Return a weight that gives phase-space distribution double ScatteredElectron_xy::Generate(const LorentzVector& parent, const particle_ptrs& products) { - + double xx,yy; std::tie(xx,yy) = _random_xy.SamplePair(); CompleteGivenXandY(parent, products, xx, yy); @@ -142,9 +142,9 @@ namespace elSpectro{ //still in rest system of nucl, just need rotation //std::cout<<"1 ScatteredElectron_xy::CompleteGivenXandY"<<_parent_in_elFrame<<" "<<_scattered<<" "<<_parent_in_elFrame -_scattered<<" W "<<(_parent_in_elFrame -_scattered).M()<Pdg()==11){ products[0]->SetP4(_scattered); products[1]->SetP4(parent -_scattered); @@ -153,6 +153,7 @@ namespace elSpectro{ products[1]->SetP4(_scattered); products[0]->SetP4(parent -_scattered); } + return 1.; } } diff --git a/core/ScatteredElectron_xy.h b/core/ScatteredElectron_xy.h index 0a967ea..21539ab 100644 --- a/core/ScatteredElectron_xy.h +++ b/core/ScatteredElectron_xy.h @@ -54,7 +54,8 @@ namespace elSpectro{ LorentzVector _parent_in_elFrame; ROOT::Math::RotationZYX _sc_rotateToZaxis; - LorentzVector _cachedParent; + ROOT::Math::RotationZ _sc_rotateAroundZaxis;//save memory allocation + LorentzVector _cachedParent; // TH1D *hist, *histW; TH1D histy={"ydist","ydist",1000,0,1}; @@ -72,9 +73,13 @@ namespace elSpectro{ if(_cachedParent!=parent){ //SetAngle is expensive (sin,cos calls) only call if necessary _cachedParent = parent; //_rotateToZaxis.SetAngle(_cachedParent.Theta()); - _sc_rotateToZaxis.SetComponents(-_cachedParent.Phi(),-_cachedParent.Theta(),0); + _sc_rotateToZaxis.SetComponents(_cachedParent.Phi(),_cachedParent.Theta(),0); + _sc_rotateAroundZaxis.SetAngle(_cachedParent.Phi()); + // _sc_rotateToZaxis.SetComponents(0,-_cachedParent.Theta(),0); } child=_sc_rotateToZaxis*child; + child =_sc_rotateAroundZaxis*child; + } } From 814b9a931f68ca1db9745008f05d6126776cb399 Mon Sep 17 00:00:00 2001 From: dglazier Date: Wed, 5 May 2021 11:18:09 +0100 Subject: [PATCH 16/21] remove setVertex statement --- core/Particle.h | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/core/Particle.h b/core/Particle.h index 6ee8fb0..f4e372a 100644 --- a/core/Particle.h +++ b/core/Particle.h @@ -95,8 +95,7 @@ namespace elSpectro{ void SetVertex(int vertexID,const LorentzVector* v){ _vertexID=vertexID; _vertex=v; - std::cout<<"PArticle SetVertex "<<_pdg<<" "<<_vertexID<<" "< Date: Wed, 5 May 2021 15:18:08 +0100 Subject: [PATCH 17/21] mv pointer swap after remove --- core/CollidingParticle.cpp | 6 +++--- core/ParticleManager.cpp | 2 +- core/ParticleManager.h | 7 +++++-- 3 files changed, 9 insertions(+), 6 deletions(-) diff --git a/core/CollidingParticle.cpp b/core/CollidingParticle.cpp index 778e409..36ba205 100644 --- a/core/CollidingParticle.cpp +++ b/core/CollidingParticle.cpp @@ -34,7 +34,7 @@ namespace elSpectro{ UInt_t position=0; for(auto& p:_model->Products()){ - + std::cout<<"CollidingParticle prod "<Pdg()<Pdg()==pdg){ if(_interactingParticle!=nullptr){ @@ -42,11 +42,11 @@ namespace elSpectro{ } else{ _interactingParticle=p->P4ptr(); - _model->SwapProducts(0,position); //make sure interacting particle is first in products vector } //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 diff --git a/core/ParticleManager.cpp b/core/ParticleManager.cpp index 35f7ea3..35c2d33 100644 --- a/core/ParticleManager.cpp +++ b/core/ParticleManager.cpp @@ -58,7 +58,7 @@ namespace elSpectro{ _stables.push_back(p); } - + std::cout<<"Take Partice "<Pdg()<<" "<<_stables.size()<Pdg()<<" "<<_stables.size()<Pdg()<<" "<<_stables.size()< Date: Wed, 5 May 2021 15:19:59 +0100 Subject: [PATCH 18/21] mv pointer swap after remove --- core/CollidingParticle.cpp | 1 - core/ParticleManager.cpp | 1 - core/ParticleManager.h | 9 +++------ 3 files changed, 3 insertions(+), 8 deletions(-) diff --git a/core/CollidingParticle.cpp b/core/CollidingParticle.cpp index 36ba205..87cb86b 100644 --- a/core/CollidingParticle.cpp +++ b/core/CollidingParticle.cpp @@ -34,7 +34,6 @@ namespace elSpectro{ UInt_t position=0; for(auto& p:_model->Products()){ - std::cout<<"CollidingParticle prod "<Pdg()<Pdg()==pdg){ if(_interactingParticle!=nullptr){ diff --git a/core/ParticleManager.cpp b/core/ParticleManager.cpp index 35c2d33..5bae37f 100644 --- a/core/ParticleManager.cpp +++ b/core/ParticleManager.cpp @@ -58,7 +58,6 @@ namespace elSpectro{ _stables.push_back(p); } - std::cout<<"Take Partice "<Pdg()<<" "<<_stables.size()<Pdg()<<" "<<_stables.size()<Pdg()<<" "<<_stables.size()< Date: Wed, 12 May 2021 09:55:52 +0100 Subject: [PATCH 19/21] debug VectorSDME and finish TensorSDME --- core/CMakeLists.txt | 4 + core/DecayModelQ2W.cpp | 4 +- core/DecayModelst.cpp | 8 +- core/DecayingParticle.cpp | 3 +- core/ElSpectroLinkDef.h | 2 + core/JpacModelst.cpp | 52 +++++++++--- core/PhotonPolarisationVector.h | 10 +-- core/SDME.h | 2 +- core/SDMEDecay.cpp | 33 ++++++++ core/SDMEDecay.h | 56 +++++++++++++ core/TensorSDMEDecay.cpp | 138 ++++++++++++++++++++++++++++++++ core/TensorSDMEDecay.h | 38 +++++++++ core/VectorSDMEDecay.cpp | 63 +++++++-------- core/VectorSDMEDecay.h | 25 ++---- 14 files changed, 365 insertions(+), 73 deletions(-) create mode 100644 core/SDMEDecay.cpp create mode 100644 core/SDMEDecay.h create mode 100644 core/TensorSDMEDecay.cpp create mode 100644 core/TensorSDMEDecay.h diff --git a/core/CMakeLists.txt b/core/CMakeLists.txt index c01d957..8becf00 100644 --- a/core/CMakeLists.txt +++ b/core/CMakeLists.txt @@ -12,6 +12,8 @@ ROOT_GENERATE_DICTIONARY(G__${ELSPECTRO} DecayModel.h PhaseSpaceDecay.h NuclearBreakup.h + SDMEDecay.h + TensorSDMEDecay.h VectorSDMEDecay.h DecayModelQ2W.h DecayModelst.h @@ -55,6 +57,8 @@ add_library(${ELSPECTRO} SHARED DecayModel.cpp PhaseSpaceDecay.cpp NuclearBreakup.cpp + SDMEDecay.cpp + TensorSDMEDecay.cpp VectorSDMEDecay.cpp DecayModelQ2W.cpp DecayModelst.cpp diff --git a/core/DecayModelQ2W.cpp b/core/DecayModelQ2W.cpp index 5b44444..929a568 100644 --- a/core/DecayModelQ2W.cpp +++ b/core/DecayModelQ2W.cpp @@ -115,8 +115,9 @@ namespace elSpectro{ auto delta = 2*escat::M2_el()/getQ2()*(1-epsilon); _photonPol.SetEpsilon(epsilon); + //_photonPol.SetEpsilon(1); _photonPol.SetDelta(delta); - + //Get envelope weight from integrated cross section double weight=_Wrealphoto_Dist->GetWeightFor( W ); @@ -225,7 +226,6 @@ namespace elSpectro{ //on efficiency from this if(val #include @@ -32,7 +33,8 @@ namespace elSpectro{ if( dynamic_cast(_meson) ){ if( dynamic_cast(_meson)->Model()->CanUseSDME() ){ - _sdmeMeson = _meson->InitSDME(1,4); + auto sdmeModel=dynamic_cast(dynamic_cast(_meson)->Model()); + _sdmeMeson = _meson->InitSDME(sdmeModel->Spin(),4); //4=>photoproduction //could have electroproduced baryon spin 3/2 //_sdmeBaryon = _baryon->InitSDME(3,9); } @@ -424,9 +426,9 @@ namespace elSpectro{ _W=hist.GetXaxis()->GetBinCenter(ih); if( _W < Wmin ) hist.SetBinContent(ih, 0); - if( TMath::IsNaN(kine::tmax(_W,M1,M2,M3,M4)) ) + else if( TMath::IsNaN(kine::tmax(_W,M1,M2,M3,M4)) ) hist.SetBinContent(ih, 0); - if( TMath::IsNaN(kine::t0(_W,M1,M2,M3,M4)) ) + else if( TMath::IsNaN(kine::t0(_W,M1,M2,M3,M4)) ) hist.SetBinContent(ih, 0); else{ double max_at_W=0; diff --git a/core/DecayingParticle.cpp b/core/DecayingParticle.cpp index c12c23d..936a267 100644 --- a/core/DecayingParticle.cpp +++ b/core/DecayingParticle.cpp @@ -85,7 +85,7 @@ namespace elSpectro{ //evaluate the model intensity for the product vectors double weight = 1; if(Model()!=nullptr) weight = Model()->Intensity(); - //std::cout<<"DecayingParticle::GenerateProducts "<GetName()<<" "<Products()[0]->Mass()<<" "<Products()[1]->Mass()<<" difference in weights "< gRandom->Uniform()*_maxWeight ; if (decayed == false && (Model()->RegenerateOnFail()==false) ) return DecayStatus::TryAnother; diff --git a/core/ElSpectroLinkDef.h b/core/ElSpectroLinkDef.h index 6421e2c..313bdf4 100644 --- a/core/ElSpectroLinkDef.h +++ b/core/ElSpectroLinkDef.h @@ -22,6 +22,8 @@ #pragma link C++ class elSpectro::DecayGammaN_Test+; #pragma link C++ class elSpectro::PhaseSpaceDecay+; #pragma link C++ class elSpectro::NuclearBreakup+; +#pragma link C++ class elSpectro::SDMEDecay+; +#pragma link C++ class elSpectro::TensorSDMEDecay+; #pragma link C++ class elSpectro::VectorSDMEDecay+; #pragma link C++ class elSpectro::DecayVectors+; diff --git a/core/JpacModelst.cpp b/core/JpacModelst.cpp index 9689399..423709e 100644 --- a/core/JpacModelst.cpp +++ b/core/JpacModelst.cpp @@ -33,15 +33,49 @@ namespace elSpectro{ auto *sdme=GetMesonSDMEs(); if(sdme){ //note this is vector formalism - sdme->SetElement(0,0,0,(_amp->SDME(0, 0, 0, get_s(), get_t()))); - sdme->SetElement(0,1,0,(_amp->SDME(0, 1, 0, get_s(), get_t()))); - sdme->SetElement(0,1,-1,(_amp->SDME(0, 1, -1, get_s(), get_t()))); - sdme->SetElement(1,1,1,(_amp->SDME(1, 1, 1, get_s(), get_t()))); - sdme->SetElement(1,0,0,(_amp->SDME(1, 0, 0, get_s(), get_t()))); - sdme->SetElement(1,1,0,(_amp->SDME(1, 1, 0, get_s(), get_t()))); - sdme->SetElement(1,1,-1,(_amp->SDME(1, 1, -1, get_s(), get_t()))); - sdme->SetElement(2,1,0,(_amp->SDME(2, 1, 0, get_s(), get_t()))); - sdme->SetElement(2,1,-1,(_amp->SDME(2, 1, -1, get_s(), get_t()))); + if(sdme->Spin()==1){ + sdme->SetElement(0,0,0,(_amp->SDME(0, 0, 0, get_s(), get_t()))); + sdme->SetElement(0,1,0,(_amp->SDME(0, 1, 0, get_s(), get_t()))); + sdme->SetElement(0,1,-1,(_amp->SDME(0, 1, -1, get_s(), get_t()))); + sdme->SetElement(1,1,1,(_amp->SDME(1, 1, 1, get_s(), get_t()))); + sdme->SetElement(1,0,0,(_amp->SDME(1, 0, 0, get_s(), get_t()))); + sdme->SetElement(1,1,0,(_amp->SDME(1, 1, 0, get_s(), get_t()))); + sdme->SetElement(1,1,-1,(_amp->SDME(1, 1, -1, get_s(), get_t()))); + sdme->SetElement(2,1,0,(_amp->SDME(2, 1, 0, get_s(), get_t()))); + sdme->SetElement(2,1,-1,(_amp->SDME(2, 1, -1, get_s(), get_t()))); + } + else if(sdme->Spin()==2){ + sdme->SetElement(0,0,0,(_amp->SDME(0, 0, 0, get_s(), get_t()))); + sdme->SetElement(0,1,0,(_amp->SDME(0, 1, 0, get_s(), get_t()))); + sdme->SetElement(0,1,-1,(_amp->SDME(0, 1, -1, get_s(), get_t()))); + sdme->SetElement(0,1,1,(_amp->SDME(0, 1, 1, get_s(), get_t()))); + sdme->SetElement(0,2,-1,(_amp->SDME(0, 2, -1, get_s(), get_t()))); + sdme->SetElement(0,2,-2,0.25); + //sdme->SetElement(0,2,-2,(_amp->SDME(0, 2, -2, get_s(), get_t()))); + sdme->SetElement(0,2,0,(_amp->SDME(0, 2, 0, get_s(), get_t()))); + sdme->SetElement(0,2,1,(_amp->SDME(0, 2, 1, get_s(), get_t()))); + sdme->SetElement(0,2,2,(_amp->SDME(0, 2, 2, get_s(), get_t()))); + + sdme->SetElement(1,0,0,(_amp->SDME(0, 0, 0, get_s(), get_t()))); + sdme->SetElement(1,1,0,(_amp->SDME(0, 1, 0, get_s(), get_t()))); + sdme->SetElement(1,1,-1,(_amp->SDME(0, 1, -1, get_s(), get_t()))); + sdme->SetElement(1,1,1,(_amp->SDME(0, 1, 1, get_s(), get_t()))); + sdme->SetElement(1,2,-1,(_amp->SDME(0, 2, -1, get_s(), get_t()))); + sdme->SetElement(1,2,-2,(_amp->SDME(0, 2, -2, get_s(), get_t()))); + sdme->SetElement(1,2,0,(_amp->SDME(0, 2, 0, get_s(), get_t()))); + sdme->SetElement(1,2,1,(_amp->SDME(0, 2, 1, get_s(), get_t()))); + sdme->SetElement(1,2,2,(_amp->SDME(0, 2, 2, get_s(), get_t()))); + + sdme->SetElement(2,1,0,(_amp->SDME(0, 1, 0, get_s(), get_t()))); + sdme->SetElement(2,1,-1,(_amp->SDME(0, 1, -1, get_s(), get_t()))); + sdme->SetElement(2,2,-1,(_amp->SDME(0, 2, -1, get_s(), get_t()))); + sdme->SetElement(2,2,-2,(_amp->SDME(0, 2, -2, get_s(), get_t()))); + sdme->SetElement(2,2,0,(_amp->SDME(0, 2, 0, get_s(), get_t()))); + sdme->SetElement(2,2,1,(_amp->SDME(0, 2, 1, get_s(), get_t()))); + + + + } } diff --git a/core/PhotonPolarisationVector.h b/core/PhotonPolarisationVector.h index cd0aa7a..2d3065a 100644 --- a/core/PhotonPolarisationVector.h +++ b/core/PhotonPolarisationVector.h @@ -27,7 +27,7 @@ namespace elSpectro{ void SetEpsilonDeltaPhi(double eps,double delta,double phi){ SetEpsilonPhi(eps,phi);//[1] and [2] auto sinphi=TMath::Sin(phi); - auto cosphi=TMath::Sqrt(1- sinphi*sinphi); + auto cosphi=TMath::Cos(phi); _elements[4] = eps+delta; auto factor= TMath::Sqrt(2*eps*(1+eps+2*delta)); @@ -36,10 +36,10 @@ namespace elSpectro{ } void SetEpsilonPhi(double eps,double phi){ - auto sin2phi=TMath::Sin(2*phi); - auto cos2phi=TMath::Sqrt(1- sin2phi*sin2phi); - _elements[1] = - eps*cos2phi; - _elements[2] = - eps*sin2phi; + //auto sin2phi=TMath::Sin(2*phi); + //auto cos2phi=TMath::Sqrt(1- sin2phi*sin2phi); + _elements[1] = - eps*TMath::Cos(2*phi); + _elements[2] = - eps*TMath::Sin(2*phi); } diff --git a/core/SDME.h b/core/SDME.h index 19d426b..51a318c 100644 --- a/core/SDME.h +++ b/core/SDME.h @@ -43,7 +43,7 @@ namespace elSpectro{ _elements[alpha][x][_J+y] = val; } - + uint Spin()const {return _J;} private: uint _J={0}; diff --git a/core/SDMEDecay.cpp b/core/SDMEDecay.cpp new file mode 100644 index 0000000..c7980cf --- /dev/null +++ b/core/SDMEDecay.cpp @@ -0,0 +1,33 @@ +#include "SDMEDecay.h" +#include +#include + +namespace elSpectro{ + + SDMEDecay::SDMEDecay( particle_ptrs ps, const std::vector pdgs): + DecayModel{ps,pdgs} + { + + _name={"SDMEDecay"}; + + } + + /////////////////////////////////////////////////////////////// + void SDMEDecay::PostInit(ReactionInfo* info){ + DecayModel::PostInit(info); + auto prodInfo= dynamic_cast (info); //I need Reaction info + + _rho=Parent()->GetSDME(); //get pointer to SDME values + + if(prodInfo->_meson!=Parent()->P4ptr() ){ + std::cerr<<"SDMEDecay::PostInit this is not the photoproduced meson! I have pdg # "<Pdg()<<" while the photoproduced meson was "<_meson; + _baryon = prodInfo->_baryon; + _photon = prodInfo->_photon; + _photonPol = prodInfo->_photonPol; + _child1 = Products()[0]->P4ptr(); + } + +} diff --git a/core/SDMEDecay.h b/core/SDMEDecay.h new file mode 100644 index 0000000..59e1b0d --- /dev/null +++ b/core/SDMEDecay.h @@ -0,0 +1,56 @@ +////////////////////////////////////////////////////////////// +/// +///Class: SDMEDecay +///Description: +/// Calculate intensity based on vector SDME +/// See Schilling and Wolf formalism + +#pragma once + +#include "DecayModel.h" +#include "DecayingParticle.h" +#include "SDME.h" +#include "PhotonPolarisationVector.h" + +namespace elSpectro{ + + + class SDMEDecay : public DecayModel { + + public: + + SDMEDecay()=default; + //only declaring default constructor + //so other 5 constructors also defaulted(rule of 5) + //constructor to decay into particles + SDMEDecay( particle_ptrs , const std::vector pdgs ); + + // Each model must define its intensity + double Intensity() const override =0; + + //Keep trying with new DecayVectors until pass + bool RegenerateOnFail() const noexcept final {return true;} + + void PostInit(ReactionInfo* info) override; + + bool CanUseSDME()const noexcept final{return true;} + virtual short Spin() const { return 0;} + + + protected: + + const SDME* _rho ={nullptr}; + + LorentzVector* _photon={nullptr}; + LorentzVector* _meson={nullptr}; + LorentzVector* _baryon={nullptr}; + LorentzVector* _child1={nullptr}; + PhotonPolarisationVector* _photonPol={nullptr}; + + private: + + ClassDefOverride(elSpectro::SDMEDecay,1); //class SDMEDecay + + };//class SDMEDecay + +}//namespace elSpectro diff --git a/core/TensorSDMEDecay.cpp b/core/TensorSDMEDecay.cpp new file mode 100644 index 0000000..5aea39b --- /dev/null +++ b/core/TensorSDMEDecay.cpp @@ -0,0 +1,138 @@ +#include "TensorSDMEDecay.h" +#include +#include + +namespace elSpectro{ + + TensorSDMEDecay::TensorSDMEDecay( particle_ptrs ps, const std::vector pdgs): + SDMEDecay{ps,pdgs} + { + + _name={"TensorSDMEDecay"}; + + } + + + ////////////////////////////////////////////////////////////// + double TensorSDMEDecay::Intensity() const{ + // std::cout<<" TensorSDMEDecay::Intensity() "<<_photon<<" "<<_meson<<" "<<_baryon<<" "<<_child1<Val(1,1,1)<<" "< W={0,0,0,0,0,0,0,0}; + // Eqn E11 in https://arxiv.org/pdf/2005.01617.pdf "Exclusive tensor meson photoproduction, V. Mathieu" + + + W[0] = ( + 1./16*_rho->Re(0,0,0)*(1+3*cos2Th)*(1+3*cos2Th) + - 0.75*_rho->Re(0,1,-1)*sinSq2Th*cos2Ph + -TMath::Sqrt(3./8)* _rho->Re(0,1,0)*sin2Th*cosPh*(1+3*cos2Th) + +0.75*_rho->Re(0,1,1)*sinSq2Th + +3*_rho->Re(0,2,-1)*cosTh*sinCubeTh*cos3Ph + +3./4*_rho->Re(0,2,-2)*sinSqTh*sinSqTh*cos4Ph + +TMath::Sqrt(3./8)*_rho->Re(0,2,0)*sinSqTh*cos2Ph*(1+3*cos2Th) + -3*_rho->Re(0,2,1)*cosTh*sinCubeTh*cosPh + +0.75*_rho->Re(0,2,2)*sinSqTh*sinSqTh + ); + + W[1]=( + 1./16*_rho->Re(1,0,0)*(1+3*cos2Th)*(1+3*cos2Th) + - 0.75*_rho->Re(1,1,-1)*sinSq2Th*cos2Ph + -TMath::Sqrt(3./8)* _rho->Re(1,1,0)*sin2Th*cosPh*(1+3*cos2Th) + +0.75*_rho->Re(1,1,1)*sinSq2Th + +3*_rho->Re(1,2,-1)*cosTh*sinCubeTh*cos3Ph + +0.75*_rho->Re(1,2,-2)*sinSqTh*sinSqTh*cos4Ph + +TMath::Sqrt(3./8)*_rho->Re(1,2,0)*sinSqTh*cos2Ph*(1+3*cos2Th) + -3*_rho->Re(1,2,1)*cosTh*sinCubeTh*cosPh + +0.75*_rho->Re(1,2,2)*sinSqTh*sinSqTh + ); + //Assume 3/4i * _rho_(2,1,-1) = 3/4 *Im{rho2_1-1} as Re=0; + //i.e. i*_rho->Im(2,1,-1) cancels i in /4i + W[2]= ( + + TMath::Sqrt(3./8)* _rho->Im(2,1,0)*sin2Th*sinPh*(1+3*cos2Th) + + 0.75*_rho->Im(2,1,-1)*sinSq2Th*sin2Ph + - TMath::Sqrt(3./8)*_rho->Im(2,2,0)*sinSqTh*sin2Ph*(1+3*cos2Th) + +3*_rho->Im(2,2,1)*cosTh*sinCubeTh*sinPh + -3*_rho->Im(2,2,-1)*cosTh*sinCubeTh*sin3Ph + -0.75*_rho->Im(2,2,-2)*sinSqTh*sinSqTh*sin4Ph + ); + //Assume same form as W[2], not given in paper + W[3]= ( + + TMath::Sqrt(3./8)* _rho->Im(3,1,0)*sin2Th*sinPh*(1+3*cos2Th) + + 0.75*_rho->Im(3,1,-1)*sinSq2Th*sin2Ph + - TMath::Sqrt(3./8)*_rho->Im(3,2,0)*sinSqTh*sin2Ph*(1+3*cos2Th) + +3*_rho->Im(3,2,1)*cosTh*sinCubeTh*sinPh + -3*_rho->Im(3,2,-1)*cosTh*sinCubeTh*sin3Ph + -0.75*_rho->Im(3,2,-2)*sinSqTh*sinSqTh*sin4Ph + + ); + + //+ other elctroproduced see eqn(83-85) Schilling and Wolf for Vector + + //equation (82) Schilling and Wolf, cf eqn (29) Schilling Seyboth Wolf + double result=0; + _photonPol->Calc(); //epsilon,delta and phi should all now be set + for(uint alpha=0; alpha < 8; alpha++ ){ + result += ( W[alpha] * (*_photonPol)[alpha] ) ; + } + result/=1.05; //make sure < 1 + + if(result>1.0 || result<0){ + std::cout<<"TensorSDMEDecay invalid result "<< result<Re(0,0,0)<<" 010 "<<_rho->Re(0,1,0)<<" 01-1 "<<_rho->Re(0,1,-1)<<" 111 "<<_rho->Re(1,1,1)<<" 100 "<<_rho->Re(1,0,0)<<" 110 "<<_rho->Re(1,1,0)<<" 11-1 "<<_rho->Re(1,1,-1)<<" 210 "<<_rho->Re(2,1,0)<<" 21-1 "<<_rho->Re(2,1,-1)<<" "<Epsilon()<<" delta " <<_photonPol->Delta()*TMath::RadToDeg()<<" phi " <<_photonPol->Phi()<Re(0,0,0)<=1) <<" "<< (_rho->Re(0,0,0)>=0) <Re(0,1,-1))<=0.5*(1-_rho->Re(0,0,0)) )<Re(0,1,0)*_rho->Re(0,1,0) <=0.25*_rho->Re(0,0,0)*(2-_rho->Re(0,0,0)-_rho->Re(0,1,-1)) )<Re(1,0,0))<=_rho->Re(0,0,0)) <Re(1,1,1))<=0.5*(1-_rho->Re(0,0,0)) )<Re(1,1,-1))<=0.5*(1-_rho->Re(0,0,0))) <Re(1,1,0))<=TMath::Sqrt( 0.5*_rho->Re(0,0,0)*(1-_rho->Re(0,0,0)) )) < pdgs ); + + // Each model must define its intensity + double Intensity() const final; + + //void PostInit(ReactionInfo* info) final; + + short Spin() const final{ return 2;} + + private: + + ClassDefOverride(elSpectro::TensorSDMEDecay,1); //class TensorSDMEDecay + + };//class TensorSDMEDecay + +}//namespace elSpectro diff --git a/core/VectorSDMEDecay.cpp b/core/VectorSDMEDecay.cpp index 81c537b..306bdce 100644 --- a/core/VectorSDMEDecay.cpp +++ b/core/VectorSDMEDecay.cpp @@ -5,30 +5,14 @@ namespace elSpectro{ VectorSDMEDecay::VectorSDMEDecay( particle_ptrs ps, const std::vector pdgs): - DecayModel{ps,pdgs} + SDMEDecay{ps,pdgs} { _name={"VectorSDMEDecay"}; } - /////////////////////////////////////////////////////////////// - void VectorSDMEDecay::PostInit(ReactionInfo* info){ - DecayModel::PostInit(info); - auto prodInfo= dynamic_cast (info); //I need Reaction info - _rho=Parent()->GetSDME(); //get pointer to SDME values - - if(prodInfo->_meson!=Parent()->P4ptr() ){ - std::cerr<<"VectorSDMEDecay::PostInit this is not the photoproduced meson! I have pdg # "<Pdg()<<" while the photoproduced meson was "<_meson; - _baryon = prodInfo->_baryon; - _photon = prodInfo->_photon; - _photonPol = prodInfo->_photonPol; - _child1 = Products()[0]->P4ptr(); - } ////////////////////////////////////////////////////////////// double VectorSDMEDecay::Intensity() const{ // std::cout<<" VectorSDMEDecay::Intensity() "<<_photon<<" "<<_meson<<" "<<_baryon<<" "<<_child1<Im(2,1,0)*sin2Th*sinPh - - _rho->Im(2,1,-1)*sinSqTh*sin2Ph + + _rho->Im(2,1,-1)*sinSqTh*sin2Ph ); W[3]= ( TMath::Sqrt2() * _rho->Im(3,1,0)*sin2Th*sinPh - - _rho->Im(3,1,-1)*sinSqTh*sin2Ph + + _rho->Im(3,1,-1)*sinSqTh*sin2Ph ); //+ other elctroproduced see eqn(83-85) Schilling and Wolf - //equation (82) Schilling and Wolf, cf eqn (29) Schilling Seyboth Wolf - double result=0; + double result=0.; _photonPol->Calc(); //epsilon,delta and phi should all now be set for(uint alpha=0; alpha < 8; alpha++ ){ result += ( W[alpha] * (*_photonPol)[alpha] ) ; - } + } + + // std::cout<<(*_photonPol)[1]<<" "<Phi())<1 should check this + //std::cout<<" "<< W[0]<<" "<<(*_photonPol)[0]<<" "<Phi()<<" "<<0.5 + 0.5*TMath::Cos(2*_photonPol->Phi())<Re(0,1,1)<<" "<<_rho->Re(0,0,0)<<" "<<_rho->Re(0,1,-1)<Re(1,1,1)<<" "<<_rho->Re(1,0,0)<<" "<<_rho->Re(1,1,-1)<Phi()); + //return result = 0.5 + 0.5*TMath::Cos(2*_photonPol->Phi()); + //std::cout<<"result "<1.0 || result<0){ std::cout<<"VectorSDMEDecay invalid result "<< result<Re(0,0,0)<<" 010 "<<_rho->Re(0,1,0)<<" 01-1 "<<_rho->Re(0,1,-1)<<" 111 "<<_rho->Re(1,1,1)<<" 100 "<<_rho->Re(1,0,0)<<" 110 "<<_rho->Re(1,1,0)<<" 11-1 "<<_rho->Re(1,1,-1)<<" 210 "<<_rho->Re(2,1,0)<<" 21-1 "<<_rho->Re(2,1,-1)<<" "<Epsilon()<<" delta " <<_photonPol->Delta()*TMath::RadToDeg()<<" phi " <<_photonPol->Phi()<Epsilon()<<" delta " <<_photonPol->Delta()<<" phi " <<_photonPol->Phi()*TMath::RadToDeg()<Re(0,0,0)<=1) <<" "<< (_rho->Re(0,0,0)>=0) <Re(0,1,-1))<=0.5*(1-_rho->Re(0,0,0)) )<Re(0,1,0)*_rho->Re(0,1,0) <=0.25*_rho->Re(0,0,0)*(2-_rho->Re(0,0,0)-_rho->Re(0,1,-1)) )< Date: Mon, 17 May 2021 09:09:28 +0100 Subject: [PATCH 20/21] Tensor OK --- core/JpacModelst.cpp | 4 +--- core/TensorSDMEDecay.cpp | 18 ++---------------- 2 files changed, 3 insertions(+), 19 deletions(-) diff --git a/core/JpacModelst.cpp b/core/JpacModelst.cpp index 423709e..aae7a60 100644 --- a/core/JpacModelst.cpp +++ b/core/JpacModelst.cpp @@ -49,9 +49,7 @@ namespace elSpectro{ sdme->SetElement(0,1,0,(_amp->SDME(0, 1, 0, get_s(), get_t()))); sdme->SetElement(0,1,-1,(_amp->SDME(0, 1, -1, get_s(), get_t()))); sdme->SetElement(0,1,1,(_amp->SDME(0, 1, 1, get_s(), get_t()))); - sdme->SetElement(0,2,-1,(_amp->SDME(0, 2, -1, get_s(), get_t()))); - sdme->SetElement(0,2,-2,0.25); - //sdme->SetElement(0,2,-2,(_amp->SDME(0, 2, -2, get_s(), get_t()))); + sdme->SetElement(0,2,-2,(_amp->SDME(0, 2, -2, get_s(), get_t()))); sdme->SetElement(0,2,0,(_amp->SDME(0, 2, 0, get_s(), get_t()))); sdme->SetElement(0,2,1,(_amp->SDME(0, 2, 1, get_s(), get_t()))); sdme->SetElement(0,2,2,(_amp->SDME(0, 2, 2, get_s(), get_t()))); diff --git a/core/TensorSDMEDecay.cpp b/core/TensorSDMEDecay.cpp index 5aea39b..c0fa598 100644 --- a/core/TensorSDMEDecay.cpp +++ b/core/TensorSDMEDecay.cpp @@ -16,12 +16,11 @@ namespace elSpectro{ ////////////////////////////////////////////////////////////// double TensorSDMEDecay::Intensity() const{ // std::cout<<" TensorSDMEDecay::Intensity() "<<_photon<<" "<<_meson<<" "<<_baryon<<" "<<_child1<Val(1,1,1)<<" "< W={0,0,0,0,0,0,0,0}; // Eqn E11 in https://arxiv.org/pdf/2005.01617.pdf "Exclusive tensor meson photoproduction, V. Mathieu" From b9cf71087cd2297e1d90c8f3c2439d28fd4e7aee Mon Sep 17 00:00:00 2001 From: dglazier Date: Mon, 31 May 2021 10:11:35 +0100 Subject: [PATCH 21/21] fix _maxPossiblexRange to be independent of given theta and Q2 ranges for correct sampling behaviour --- core/DistVirtPhotFlux_xy.cpp | 7 +++++-- core/Interface.h | 10 ++++++++++ jpacPhoto | 2 +- 3 files changed, 16 insertions(+), 3 deletions(-) diff --git a/core/DistVirtPhotFlux_xy.cpp b/core/DistVirtPhotFlux_xy.cpp index 0a84458..7b1cc4c 100644 --- a/core/DistVirtPhotFlux_xy.cpp +++ b/core/DistVirtPhotFlux_xy.cpp @@ -81,8 +81,11 @@ namespace elSpectro{ ymin = TMath::Exp(_lnymin); double ymax = TMath::Exp(_lnymax); for(int i=0;i<1E5;i++){ - double avail_xmin = XMin(static_cast(i*(ymax-ymin) + ymin)/1E5); - if(avail_xmin<_maxPossiblexRange) + // double avail_xmin = XMin(static_cast(i*(ymax-ymin) + ymin)/1E5); + //Need xmin with no experiment based limits on Q2 or theta + double yformin=static_cast(i*(ymax-ymin) + ymin)/1E5; + double avail_xmin =escat::M2_el()*yformin/(2*_mTar*_ebeam)/(1-yformin); + if(avail_xmin<_maxPossiblexRange) _maxPossiblexRange=avail_xmin; } diff --git a/core/Interface.h b/core/Interface.h index 65a6b8d..de9b558 100644 --- a/core/Interface.h +++ b/core/Interface.h @@ -129,5 +129,15 @@ namespace elSpectro{ return generator().Reaction(); } + inline ProductionProcess* eic(CollidingParticle* el,CollidingParticle* pr,DecayModelQ2W *totalXsec=nullptr){ + + if(totalXsec!=nullptr){ + generator().Reaction(new ElectronScattering(el,pr, totalXsec )); + } + else + generator().Reaction(new ElectronScattering(el,pr,totalXsec)); + + return generator().Reaction(); + } }//namespace elSpectro diff --git a/jpacPhoto b/jpacPhoto index c67402b..8ea3993 160000 --- a/jpacPhoto +++ b/jpacPhoto @@ -1 +1 @@ -Subproject commit c67402be5965456b6c34117324aafda54da33001 +Subproject commit 8ea3993d145fd0130cdb95a08a7dc2657f0ea8ac