diff --git a/core/BremstrPhoton.cpp b/core/BremstrPhoton.cpp index eebeb38..817a77f 100644 --- a/core/BremstrPhoton.cpp +++ b/core/BremstrPhoton.cpp @@ -39,7 +39,7 @@ namespace elSpectro{ //Also add the random phi angle... // BoostToParentWithRandPhi(parent,_gamma); products[0]->SetP4( _photon ); - + // products[1]->SetP4(parent - _photon); //electron return _weight; } diff --git a/core/CMakeLists.txt b/core/CMakeLists.txt index 49e00b4..f7bad02 100644 --- a/core/CMakeLists.txt +++ b/core/CMakeLists.txt @@ -20,6 +20,7 @@ ROOT_GENERATE_DICTIONARY(G__${ELSPECTRO} DecayModelQ2W.h DecayModelst0.h DecayModelst.h + DecayModelDnpee.h JpacModelst.h GenericModelst.h DecayVectors.h @@ -73,6 +74,7 @@ add_library(${ELSPECTRO} SHARED DecayModelst.cpp DecayModelst0.cpp JpacModelst.cpp + DecayModelDnpee.cpp GenericModelst.cpp DecayVectors.cpp DontDecay.cpp diff --git a/core/DecayModel.cpp b/core/DecayModel.cpp index f243aa0..7b9ed8e 100644 --- a/core/DecayModel.cpp +++ b/core/DecayModel.cpp @@ -85,8 +85,9 @@ namespace elSpectro{ } double DecayModel::PhaseSpaceWeightSq(double W){ + //std::cout<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. @@ -108,10 +109,12 @@ namespace elSpectro{ //Note in case there are additional unstable particle we //must subtract off their minimum masses - p->DetermineDynamicMass(-1,TCM-_unstableReservedMass[iu++]); + // std::cout<DetermineDynamicMass(-1,TCM-_unstableReservedMass[iu++]); TCM-=p->Mass(); + if(TCM<0){ - std::cout<<"DecayModel::PhaseSpaceWeightSq "<Pdg()<<" "<Mass()<<" "< +#include +#include +#include +#include +#include +#include + +namespace elSpectro{ + /////////////////////////////////////////////////////// + ///constructor includes subseqent decay of Ngamma* system + ///fix decay products + DecayModelDnpee::DecayModelDnpee() : + //_proton(new Particle{2212}),_neutron(new Particle{2112}),_ele(new Particle{11}),_pos(new Particle{-11}), + // DecayModel( {_proton,_neutron,_ele,_pos},{} ) + DecayModel( {},{2212,2112,11,-11} ) + { + _name={"DecayModelDnpee"}; + + if(Products().size()!=4){ + Fatal("DecayModelDnpee","Can only have four decay particles"); + } + + auto products = StableProducts(); + for(auto* prod:products){ + if(prod->Pdg()==2212) _proton=prod; + if(prod->Pdg()==2112) _neutron=prod; + if(prod->Pdg()==11) _ele=prod; + if(prod->Pdg()==-11) _pos=prod; + } + } + + void DecayModelDnpee::SetParent(DecayingParticle* pa){ + DecayModel::SetParent(pa); + //now can make cascade of decays + _phaseSpace.SetParentAndProducts(Parent(),{_proton,_neutron,_ele,_pos},{}); + // std::cout<<"DecayModelDnpee::SetParent ResetProducts "<<_phaseSpace.Products().size()<<" "<<_phaseSpace.Product(0)->Pdg()<<" "<<_phaseSpace.Product(0)->MinimumMassPossible()<<" "<<_phaseSpace.Product(1)->Pdg()< (info); //I need Reaction info + if(_prodInfo==nullptr){ + _isElProd=kFALSE; + _prodInfo= dynamic_cast (info); + } + _photon = _prodInfo->_photon; + _target = _prodInfo->_target; + _ebeam = _prodInfo->_ebeam; + + // double maxW = ( *(_prodInfo->_target) + *(_prodInfo->_ebeam) ).M(); + _Wmax = _prodInfo->_Wmax; + + _max = FindMaxOfIntensity()*1.08; //add 5% for Q2,meson mass effects etc. + + std::cout<<"DecayModelDnpee::PostInit max value "<<_max<P4().M(); + _s=_W*_W; + + // _dvolume = _dt; //NEED TO DEFINE PHASE SPACE RANGE AT CURRENT VALUE OF S. DEPENDS ON VARIABLES etc + _dvolume=(_W-_Wmin)*(_W-_Wmin);//THSIS IS NOT RIGHT, JUST A PROXY!!! + + double weight = DifferentialXSect() * _dvolume ;//must multiply by ranges for correct sampling + //std::cout<<"********************DecayModelDnpee::Intensity() "<<_W<<" weight = "<_sWeight; + // std::cout<<"********************DecayModelDnpee::Intensity() "<<_W<<" "<_sWeight<_Wmax; + _Wmin = Parent()->MinimumMassPossible(); + double Wrange=_Wmax-_Wmin; + std::cout<<"Search "<<_Wmin<<" to "<<_Wmax<SetMaxFunctionCalls(1000000); // for Minuit/Minuit2 + minimum->SetMaxIterations(1000); // for GSL + minimum->SetTolerance(0.0001); + minimum->SetPrintLevel(0); + + // create function wrapper for minimizer + // a IMultiGenFunction type + ROOT::Math::Functor wrapf(Fmax,1); //currently only dependent on W + double step[1] = {Wrange/100}; + double variable[1] = {_Wmin+Wrange/2}; + + minimum->SetFunction(wrapf); + + // Set the free variables to be minimized ! + minimum->SetVariable(0,"W",variable[0], step[0]); + + // do the minimization + minimum->Minimize(); + const double *xs = minimum->X(); + + auto minVal = -minimum->MinValue(); + auto minW= xs[0]; + std::cout << "Maximum : Probabiltiy Dist at ( W=" << minW << "): "<< minVal << std::endl; + return minVal; + } +} diff --git a/core/DecayModelDnpee.h b/core/DecayModelDnpee.h new file mode 100644 index 0000000..f104b8c --- /dev/null +++ b/core/DecayModelDnpee.h @@ -0,0 +1,109 @@ +////////////////////////////////////////////////////////////// +/// +///Class: DecayModelDnpee +///Description: +/// Control behaviour of npe+e- final states on deuteron +/// +/// +/// Note derived classes should include a constructor to initialise +/// DecayModelDnpee( particle_ptrs , const std::vector pdgs ); +#pragma once + +#include "DecayModel.h" +#include "PhaseSpaceDecay.h" +#include "FunctionsForElectronScattering.h" +#include "DecayingParticle.h" +#include + +namespace elSpectro{ + + + class DecayModelDnpee : public DecayModel { + + public: + + DecayModelDnpee(); + + // + double Intensity() const override; //this should perhaps be final so derived classes cannot overwrite... + + void PostInit(ReactionInfo* info) override; + + bool RegenerateOnFail() const noexcept override {return true;}; + bool HasAngularDistribution() override{return false; } + + const Particle* GetNeutron() const noexcept{return _neutron; } + const Particle* GetProton() const noexcept{return _proton; } + const Particle* GetElectron() const noexcept{return _ele; } + const Particle* GetPositron() const noexcept{return _pos; } + + void SetParent(DecayingParticle* pa) override ; + double PhaseSpaceWeightSq(double W) override{ + return _phaseSpace.PhaseSpaceWeightSq(W); + } + + double PhaseSpaceFactor() const noexcept { + // auto fluxPhaseSpace = p1*_W;//eqn 47.28b https://pdg.lbl.gov/2019/reviews/rpp2019-rev-kinematics.pdf + return PhaseSpaceNorm()/_s/PgammaCMsq(); + + } + double PgammaCMsq() const noexcept{ + // std::cout<<"PgammaCMsq() "<M())<<" "<<_W<<" "<<_target->M()<M()==0) return kine::PDK2(_W,0,_target->M()); + auto pgammaCM= PgammaCM(); + return pgammaCM* pgammaCM; //for dt phase space factor + } + + double PgammaCM()const noexcept{ + //in case no photon 4-vector yet + if(_photon->M()==0) return kine::PDK(_W,0,_target->M()); + //else PDK does not qork for virtual photons + auto cmBoost=Parent()->P4().BoostToCM(); + auto p1cm=boost(*_photon,cmBoost); + + // std::cout<<"PgammaCMsq M "<<_photon->M()<<" PLAB "<<_photon->P()<<" PCM "<M(),_target->M())<<" or "<M2(),_target->M())<<" or "<M())< nb + + virtual double DifferentialXSect() const{//dont let others call this as need _s, _W and _t set + //Note if your derived model already gives differential cross section + //you will need to divide by PhaseSpaceFactor to get MatrixElementSquared from it + //std::cout<<" DifferentialXSect() "<<_W<<" "<Epsilon()+_photonPol->Delta())*MatrixElementsSquared_L()); + } + + Particle* _neutron={nullptr}; + Particle* _proton={nullptr}; + Particle* _ele={nullptr}; + Particle* _pos={nullptr}; + + LorentzVector* _photon={nullptr}; + LorentzVector* _target={nullptr};//{0,0,0,escat::M_pr()}; + const LorentzVector* _ebeam={nullptr};//{0,0,0,escat::M_pr()}; + + PhaseSpaceDecay _phaseSpace; + ReactionPhotoProd* _prodInfo={nullptr}; + + double _Wmax={0}; + double _Wmin={0}; + mutable double _max={0}; + mutable double _dvolume={1}; + mutable double _s={0}; + mutable double _t={0}; + mutable double _W={0}; + + bool _isElProd={true}; + + ClassDefOverride(elSpectro::DecayModelDnpee,1); //class DecayModelst + + }; + +} + diff --git a/core/DecayModelQ2W.cpp b/core/DecayModelQ2W.cpp index 734c2ab..37a8234 100644 --- a/core/DecayModelQ2W.cpp +++ b/core/DecayModelQ2W.cpp @@ -136,7 +136,7 @@ namespace elSpectro{ //Q2 dependence of cross section weight*=Q2H1Rho(); - //std::cout<<" Q things "<_sWeight<<" "<_sWeight<<" "< #include "TFile.h" @@ -62,9 +63,12 @@ namespace elSpectro{ std::cout<<"DecayModelW::PostInit with W threshold = "<(GetGammaN()->Model()) == nullptr){ + _Wrealphoto_Dist.reset( new DistConst(1) ); + } + else{ + FindExcitationSpectra(); + } } //////////////////////////////////////////////////////// @@ -80,17 +84,9 @@ namespace elSpectro{ if(W < _threshold ) return 0.; - // auto epsilon = escat::virtualPhotonPolarisation(p4beam,p4tar,p4scat); - //auto delta = 2*escat::M2_el()/getQ2()*(1-epsilon); - - // _photonPol.SetEpsilon(epsilon); - // _photonPol.SetDelta(delta); - //Get envelope weight from integrated cross section double weight=_Wrealphoto_Dist->GetWeightFor( W ); - - - + _prodInfo->_sWeight=weight; //might be used in s and t diff --git a/core/DecayModelW.h b/core/DecayModelW.h index ed79d1a..47056f2 100644 --- a/core/DecayModelW.h +++ b/core/DecayModelW.h @@ -64,7 +64,8 @@ namespace elSpectro{ } void FindExcitationSpectra(); - DistTH1* GetApproxWDist() const {return _Wrealphoto_Dist.get();} + //DistTH1* GetApproxWDist() const {return _Wrealphoto_Dist.get();} + Distribution* GetApproxWDist() const {return _Wrealphoto_Dist.get();} double dsigma() const override { return dynamic_cast(_gstarNuc)->Model()->dsigma();} @@ -93,7 +94,8 @@ namespace elSpectro{ ReactionPhotoProd* _prodInfo={nullptr}; TH1D _hWPhaseSpace; - std::unique_ptr _Wrealphoto_Dist; + // std::unique_ptr _Wrealphoto_Dist; + std::unique_ptr _Wrealphoto_Dist; ClassDefOverride(elSpectro::DecayModelW,1); //class DecayModelW diff --git a/core/DecayModelst.cpp b/core/DecayModelst.cpp index bd54559..9f3eac4 100644 --- a/core/DecayModelst.cpp +++ b/core/DecayModelst.cpp @@ -17,6 +17,9 @@ namespace elSpectro{ { _name={"DecayModelst"}; + if(Products().size()!=2){ + Fatal("DecayModelst","Can only have two decay particles"); + } //need to find meson and baryon if(TDatabasePDG::Instance()->GetParticle(Products()[0]->Pdg())->ParticleClass()==TString("Baryon") ){ _baryon=Products()[0]; @@ -30,7 +33,7 @@ namespace elSpectro{ } ///////////////////////////////////////////////////////////////// void DecayModelst::PostInit(ReactionInfo* info){ - std::cout<<"DecayModelst::PostInit "<(_meson) ){ if( dynamic_cast(_meson)->Model()->CanUseSDME() ){ auto sdmeModel=dynamic_cast(dynamic_cast(_meson)->Model()); @@ -39,10 +42,8 @@ namespace elSpectro{ //_sdmeBaryon = _baryon->InitSDME(3,9); } } - std::cout<<"DecayModelst::PostInit "< (info); //I need Reaction info @@ -109,7 +110,7 @@ namespace elSpectro{ //Correct for W weighting which has already been applied weight/=_prodInfo->_sWeight; - // std::cout<<" s weight "<<_prodInfo->_sWeight<<" weight "<_sWeight<<" weight "<1){ std::cout<<" s weight "<<_prodInfo->_sWeight<<" Q2 "<<-_photon->M2()<<" 2Mmu "<<2*_target->M()*_photon->E() <<" W "<<_W<<" t "<<_t<<" new weight "<_sWeight<<" meson "<<_meson->Mass()<_sWeight<<" max "<<_max<<" val "<< weight*_prodInfo->_sWeight*_max<_target) + *(_prodInfo->_ebeam) ).M(); _Wmax = _prodInfo->_Wmax; - std::cout<<" DecayModelst::FindMaxOfIntensity() Wmin = "<M()<<" PLAB "<<_photon->P()<<" PCM "<M(),_target->M())<<" or "<M2(),_target->M())<<" or "<M())<Products(); + // std::cout<<" DecayingParticle::PostInit "<PostInit(info); + if(_decayer)_decayer->PostInit(info); + if(Pdg()!=0&&Pdg()!=-2211){//for real particles Double_t productMasses = 0.0; for(auto* prod: products){ @@ -62,16 +65,17 @@ namespace elSpectro{ } auto mmp=MaximumMassPossible(); if(mmpPostInit(info); - if(_decayer)_decayer->PostInit(info); + - //std::cout<<"DecayingParticle::PostInit pdg "<Products().size()<<" "<<" "<<_decay->Products()[0]->Pdg()<<" "<<_decay->Products()[1]->Pdg()<Products().size()<<" "<<" "<<_decay->Products()[0]->Pdg()<<" "<<_decay->UnstableProducts().size()<<" "<<_decay->StableProducts().size()<<" "<HasAngularDistribution()==false)samplingWeight=1; //Model has no angular distribution - + //samplingWeight ==0 => not physical (below threshold) if(samplingWeight==0) return DecayStatus::ReGenerate; //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; @@ -129,24 +134,19 @@ namespace elSpectro{ GenerateVertexPosition(); - // std::cout<<"+++++++++++++++++DecayingParticle "<UnstableProducts(); for(auto* prod: unproducts){ DecayStatus prodStatus=DecayStatus::ReGenerate; + // std::cout<<"+++++DecayingParticle start checking unstable "<StableProducts().size()<<" "<<_decay->Products().size()<<" "<GenerateProducts()) != DecayStatus::Decayed){ + if(prodStatus==DecayStatus::ReGenerate) return DecayStatus::ReGenerate; } } - /* - //stable products need a vertex too - auto& stproducts=_decay->StableProducts(); - // std::cout<<"DecayingPArticle size "<GetName()<SetVertex(_decayVertexID,_decayVertex); - } - */ - + return DecayStatus::Decayed; } diff --git a/core/DistFlatMass.h b/core/DistFlatMass.h index 12b51a0..568351d 100644 --- a/core/DistFlatMass.h +++ b/core/DistFlatMass.h @@ -40,20 +40,23 @@ namespace elSpectro{ double MinValue() const noexcept final{return 1;} double GetMinX() const noexcept final{return 0;} - double GetMaxX() const noexcept final{return 1;} + double GetMaxX() const noexcept final{return _maxX;} double GetValueFor(double valX,double valY=0) override {return 1.;} + void SetMaxX(double val){_maxX=val;} + protected : void SetIndex(uint index){_index=index;} uint Index() const noexcept {return _index;} + private: //no one should use default constructor DistFlatMass()=default; DistFlatMassMaster* _master={nullptr}; uint _index = {0}; - + double _maxX=0; ClassDef(elSpectro::DistFlatMass,1); //class Distribution diff --git a/core/DistVirtPhotFlux_xy.cpp b/core/DistVirtPhotFlux_xy.cpp index 2c81b4b..d8b6420 100644 --- a/core/DistVirtPhotFlux_xy.cpp +++ b/core/DistVirtPhotFlux_xy.cpp @@ -63,15 +63,16 @@ namespace elSpectro{ else _lnymax=TMath::Log(1.); - double xmin =XMin(ymin); - // FindMaxVal(); - + //Seacrh for lowest possible x... _maxPossiblexRange=1; ymin = TMath::Exp(_lnymin); double ymax = TMath::Exp(_lnymax); + double xmin =XMin(ymin); + + for(int i=0;i<1E5;i++){ // double avail_xmin = XMin(static_cast(i*(ymax-ymin) + ymin)/1E5); //Need xmin with no experiment based limits on Q2 or theta @@ -80,7 +81,9 @@ namespace elSpectro{ if(avail_xmin<_maxPossiblexRange) _maxPossiblexRange=avail_xmin; } - + + FindMaxVal(); //For sampling + std::cout<<"DistVirtPhotFlux_xy::SetWThreshold new Wmin "<Uniform(TMath::Log(_maxPossiblexRange),TMath::Log(1)); else lnx = gRandom->Uniform(TMath::Log(1E-50),TMath::Log(1)); - - //check if we are within allowed x-range + + + //check if we are within allowed x-range //if not return and throw another y value auto currx=TMath::Exp(lnx); if(currx>avail_xmax){ lnx=1;return; } @@ -172,7 +177,7 @@ namespace elSpectro{ }; getRandomXY();//intial sample - + // std::cout<<"DistVirtPhotFlux_xy::max "<<_max_val<Uniform()*_max_val > (_val=escat::flux_dlnxdlny(_ebeam,lnx,lny)) ) { if(_val>_max_val){ diff --git a/core/DistVirtPhotFlux_xy.h b/core/DistVirtPhotFlux_xy.h index 1dfec74..a1cf3b8 100644 --- a/core/DistVirtPhotFlux_xy.h +++ b/core/DistVirtPhotFlux_xy.h @@ -144,7 +144,8 @@ namespace elSpectro{ //evaluate photon flux as function of lnx and lny inline double DistVirtPhotFlux_xy::Eval(const double *x) const{ - double lnx=x[0]; + + double lnx=x[0]; double lny=x[1]; if(lny<_lnymin || lny>_lnymax) return 0; @@ -161,9 +162,9 @@ namespace elSpectro{ if(currx>avail_xmax){ return 0; } if(currxGetInteracting4Vector())); _beamNucl.SetP4(*(_targetptr->GetInteracting4Vector())); - _massIon=_beamNucl.PdgMass(); + _massIon=_beamNucl.PdgMass(); /* //not sure why this was there //For decaying diff --git a/core/FunctionsForJpac.h b/core/FunctionsForJpac.h index b979e48..b5220a7 100644 --- a/core/FunctionsForJpac.h +++ b/core/FunctionsForJpac.h @@ -1,9 +1,9 @@ -#include "reaction_kinematics.hpp" -#include "regge_trajectory.hpp" -#include "amplitudes/pseudoscalar_exchange.hpp" -#include "amplitudes/pomeron_exchange.hpp" -#include "amplitudes/amplitude_sum.hpp" +#include "core/reaction_kinematics.hpp" +#include "core/regge_trajectory.hpp" +#include "core/pseudoscalar_exchange.hpp" +#include "core/pomeron_exchange.hpp" +#include "core/amplitude_sum.hpp" #include #include diff --git a/core/Interface.h b/core/Interface.h index 18c4f8d..a939444 100644 --- a/core/Interface.h +++ b/core/Interface.h @@ -132,7 +132,7 @@ namespace elSpectro{ return generator().Reaction(); } - inline PhotoProduction* photoprod(CollidingParticle* el,CollidingParticle* pr,DecayModelW *totalXsec=nullptr){ + inline PhotoProduction* photoprod(CollidingParticle* el,CollidingParticle* pr,DecayModel *totalXsec=nullptr){ if(totalXsec!=nullptr){ model(totalXsec); //register ownership of model with manager diff --git a/core/JpacModelst.h b/core/JpacModelst.h index 7083922..f13b741 100644 --- a/core/JpacModelst.h +++ b/core/JpacModelst.h @@ -14,7 +14,7 @@ #include "DecayModelst.h" #include "SDME.h" #include "FunctionsForElectronScattering.h" -#include "amplitudes/amplitude.hpp" +#include "core/amplitude.hpp" namespace elSpectro{ @@ -34,7 +34,7 @@ namespace elSpectro{ double MatrixElementsSquared_T() const override { - _amp->_kinematics->set_mX( GetMeson()->Mass() ); + _amp->_kinematics->set_meson_mass( GetMeson()->Mass() ); // _amp->_kinematics->set_Q2( get_Q2() ); //std::cout<<"me "<Mass()<<" Q2 "<< get_Q2()<<" t "<_kinematics->Wth()<<" VAL "<<_amp->probability_distribution(get_s(),get_t())/4<_kinematics->Wth()) return 0; diff --git a/core/MassPhaseSpace.h b/core/MassPhaseSpace.h index 40de52a..dc39477 100644 --- a/core/MassPhaseSpace.h +++ b/core/MassPhaseSpace.h @@ -32,6 +32,7 @@ namespace elSpectro{ double PhaseSpaceWeight(double parentM) const noexcept{ double result=TMath::Sqrt(_model->PhaseSpaceWeightSq(parentM)); _weightCalcN++; + return result; } @@ -39,7 +40,7 @@ namespace elSpectro{ //sample all masses according to overall decay phase space if(_model==nullptr) return; if(_model!=amodel) return; //only 1 model controls phasespace - + // std::cout<<"MassPhaseSpace Find"<Uniform()*max*_suppressPhaseSpace ) { - + // if(wee>max){ // _sampledMax=wee; - // std::cout<<" weight > max "< max "< gRandom->Uniform()*max ) ? true : false; diff --git a/core/Particle.h b/core/Particle.h index 3b91b41..570193f 100644 --- a/core/Particle.h +++ b/core/Particle.h @@ -134,11 +134,12 @@ namespace elSpectro{ if(_massLocked==true) return; //someone else in charge... _dynamicMass=-1; auto minposs = MinimumMassPossible(); + while(_dynamicMassGetMaxX():xmax; if(minRange>maxRange){//unphysical - std::cout<<"Warning Particle::DetermineDynamicMass min "<GetCurrentWeight(); - // std::cout<<_pdg<<" DetermineDynamicMass( "<Boost(vboost);}); diff --git a/core/PhaseSpaceDecay.cpp b/core/PhaseSpaceDecay.cpp index a92fde3..f340dd7 100644 --- a/core/PhaseSpaceDecay.cpp +++ b/core/PhaseSpaceDecay.cpp @@ -14,6 +14,11 @@ namespace elSpectro{ } + void PhaseSpaceDecay::SetParentAndProducts(DecayingParticle* pa, const particle_ptrs stable, const decaying_ptrs unstable){ + DecayModel::SetParent(pa); + nBodyDecayer(pa,stable,unstable); + } + void PhaseSpaceDecay::SetParent(DecayingParticle* pa){ DecayModel::SetParent(pa); //now can make cascade of decays @@ -22,7 +27,25 @@ namespace elSpectro{ } void PhaseSpaceDecay::PostInit(ReactionInfo* info){ + std::cout<<"PhaseSpaceDecay::PostInit " <MaximumMassPossible()<<" pdg "<Pdg()<MaximumMassPossible(); + if(maxmass==0) maxmass = info->_Wmax; //if parent has no max just use reaction W + auto flatmass=dynamic_cast(Manager::Instance().Particles().GetMassDist(Parent()->Pdg())); + if(flatmass) flatmass->SetMaxX(maxmass); + + for(auto& product:Products()){ + auto flatmass=dynamic_cast(Manager::Instance().Particles().GetMassDist(product->Pdg())); + SumAllProducts(); + std::cout<<" PhaseSpaceDecay::PostInit product " <Pdg()<<" flat mass dist ? "<(Manager::Instance().Particles().GetMassDist(product->Pdg()))<<" "<SetMaxX(maxmass); + } + + if(_massMaster!=nullptr)_massMaster->SetMaxX(maxmass); + std::cout<<"PhaseSpaceDecay::PostInit " <MassDistribution()==nullptr){//production proces does not have mass distribution if(dynamic_cast(Parent())==nullptr&&Parent()->Pdg()!=-2211){ std::cerr<<"PhaseSpaceDecay::PostInit parent needs a mass distribution for pdg = "<Pdg(); @@ -33,6 +56,7 @@ namespace elSpectro{ // exit(0); } } + } void PhaseSpaceDecay::nBodyDecayer(DecayingParticle* parent, const particle_ptrs stable, const decaying_ptrs unstable ) //take copies of particle vectors @@ -58,10 +82,10 @@ namespace elSpectro{ //Register our X-1 particle type //We will construct this particle at the end of the rest auto pdgXNm1=particleMan.RegisterNewPdgParticle(0,new DistFlatMassMaster(parent,ps)); - std::cout<<"Particle* nBodyDecayer pdg "<< pdgXNm1 <(particleMan.GetMassDist(pdgXNm1)); - std::cout<<"Particle* nBodyDecayer pdg "<< massMaster <(particleMan.GetMassDist(pdgXNm1)); + std::cout<<"Particle* nBodyDecayer pdg "<< _massMaster <Pdg()<<" and "<Pdg()<Pdg()<<" and "<Pdg()<Pdg()<<" "<Pdg()< pdgs): - SDMEDecay{ps,pdgs} + VectorSDMEDecay::VectorSDMEDecay( particle_ptrs ps, const std::vector pdgs, int decayType): + SDMEDecay{ps,pdgs},_type{decayType} { _name={"VectorSDMEDecay"}; + if(pdgs[0]==11||pdgs[0]==13||pdgs[1]==11||pdgs[1]==13){ + if(_type!=1) + Fatal("VectorSDMEDecay"," using Spin 0 decay equations, but with leptons! If you want to continue to do this , change decaytype to 100, but better to change decay type to 1"); + else if(_type==100) //lets use of spin 0 equations for leptons, Should not really be used + _type = 0; + } } @@ -41,47 +47,49 @@ namespace elSpectro{ std::array W={0,0,0,0,0,0,0,0}; // eqn (31) Schilling,Seyboth and Wolf + factor const_3by4pi()* - - W[0] = ( - 0.5 * (1 -_rho->Re(0,0,0) ) - + 0.5 * (3*_rho->Re(0,0,0) - 1) *cosSqTh - - TMath::Sqrt2() * _rho->Re(0,1,0)*sin2Th*cosPh - - _rho->Re(0,1,-1)*sinSqTh*cos2Ph); - - W[1]=( - _rho->Re(1,1,1) * sinSqTh - + _rho->Re(1,0,0) * cosSqTh - - TMath::Sqrt2() * _rho->Re(1,1,0)*sin2Th*cosPh - - _rho->Re(1,1,-1)*sinSqTh*cos2Ph - ); - W[2]= ( - TMath::Sqrt2() * _rho->Im(2,1,0)*sin2Th*sinPh - + _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 - ); - + if(_type==1){//decay to lepton e.g. e+ e- + W[0]= W0lepto(cosSqTh,sin2Th,cosPh,sinSqTh,cos2Ph)*3./8/TMath::Pi(); + W[1]= W1lepto(sinSqTh,cosSqTh,sin2Th,cosPh,cos2Ph)*3./8/TMath::Pi(); + W[2]= W2lepto(sin2Th, sinPh,sinSqTh, sin2Ph)*3./8/TMath::Pi(); + W[3]= W3lepto(sin2Th,sinPh,sinSqTh,sin2Ph)*3./8/TMath::Pi(); + } + else{ //decay to spin0 e.g. pi+ pi- + W[0]= W0spin0(cosSqTh,sin2Th,cosPh,sinSqTh,cos2Ph)*3./4/TMath::Pi(); + W[1]= W1spin0(sinSqTh,cosSqTh,sin2Th,cosPh,cos2Ph)*3./4/TMath::Pi(); + W[2]= W2spin0(sin2Th, sinPh,sinSqTh, sin2Ph)*3./4/TMath::Pi(); + W[3]= W3spin0(sin2Th,sinPh,sinSqTh,sin2Ph)*3./4/TMath::Pi(); + + } //+ other elctroproduced see eqn(83-85) Schilling and Wolf //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] ) ; - } + } // std::cout<<(*_photonPol)[1]<<" "<Phi())<1 should check this + result/=1.1; //max seems to be slightly>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)<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)<Im(2,1,1)<<" "<<_rho->Im(2,1,-1)<Re(0,0,0); + // _countRho011+=_rho->Re(0,1,1); + // _countRho01m1+=_rho->Re(0,1,-1); + // _countRho100+=_rho->Re(1,0,0); + // _countRho111+=_rho->Re(1,1,1); + // _countRho11m1+=_rho->Re(1,1,-1); + // _countRho210+=_rho->Im(2,1,0); + // _countRho211+=_rho->Im(2,1,1); + // _countRho21m1+=_rho->Im(2,1,-1); + // ++_countN; + // std::cout<<"Mean SDMEs "<<" rho000 "<<_countRho000/_countN<<" rho011 "<<_countRho011/_countN<<" rho01m1 "<<_countRho01m1/_countN<<" rho100 "<<_countRho100/_countN<<" rho111 "<<_countRho111/_countN<<" rho11m1 "<<_countRho11m1/_countN<<" rho210 "<<_countRho210/_countN<<" rho211 "<<_countRho211/_countN<<" rho21m1 "<<_countRho21m1/_countN<<(*_photonPol)[1]<<" "<<(*_photonPol)[2]<Phi()); //return result = 0.5 + 0.5*TMath::Cos(2*_photonPol->Phi()); //std::cout<<"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()<<" phi " <<_photonPol->Phi()*TMath::RadToDeg()<Re(0,0,0)<=1) <<" "<< (_rho->Re(0,0,0)>=0) < pdgs ); + + //decay type = 0 for decay to spin 0; =1 for decay to leptons + VectorSDMEDecay( particle_ptrs , const std::vector pdgs, int decayType ); // Each model must define its intensity double Intensity() const final; @@ -29,10 +31,83 @@ namespace elSpectro{ short Spin() const final{ return 1;} + double W0spin0(double cosSqTh,double sin2Th,double cosPh,double sinSqTh,double cos2Ph) const noexcept; + double W1spin0(double sinSqTh,double cosSqTh,double sin2Th, double cosPh, double cos2Ph) const noexcept; + double W2spin0(double sin2Th, double sinPh,double sinSqTh, double sin2Ph) const noexcept; + double W3spin0(double sin2Th, double sinPh,double sinSqTh, double sin2Ph) const noexcept; + + + double W0lepto(double cosSqTh,double sin2Th,double cosPh,double sinSqTh,double cos2Ph) const noexcept; + double W1lepto(double sinSqTh,double cosSqTh,double sin2Th, double cosPh, double cos2Ph) const noexcept; + double W2lepto(double sin2Th, double sinPh,double sinSqTh, double sin2Ph) const noexcept; + double W3lepto(double sin2Th, double sinPh,double sinSqTh, double sin2Ph) const noexcept; + private: - + int _type=0; + // mutable double _countRho000=0; + // mutable double _countRho011=0; + // mutable double _countRho01m1=0; + // mutable double _countRho100=0; + // mutable double _countRho111=0; + // mutable double _countRho11m1=0; + // mutable double _countRho210=0; + // mutable double _countRho211=0; + // mutable double _countRho21m1=0; + // mutable double _countN=0; + ClassDefOverride(elSpectro::VectorSDMEDecay,1); //class VectorSDMEDecay };//class VectorSDMEDecay + inline double VectorSDMEDecay::W0spin0(double cosSqTh,double sin2Th,double cosPh,double sinSqTh,double cos2Ph) const noexcept{ + + return 0.5 * (1 -_rho->Re(0,0,0) ) + + 0.5 * (3*_rho->Re(0,0,0) - 1) *cosSqTh + - TMath::Sqrt2() * _rho->Re(0,1,0)*sin2Th*cosPh + - _rho->Re(0,1,-1)*sinSqTh*cos2Ph; + } + + inline double VectorSDMEDecay::W1spin0(double sinSqTh,double cosSqTh,double sin2Th, double cosPh, double cos2Ph) const noexcept{ + + return _rho->Re(1,1,1) * sinSqTh + + _rho->Re(1,0,0) * cosSqTh + - TMath::Sqrt2() * _rho->Re(1,1,0)*sin2Th*cosPh + - _rho->Re(1,1,-1)*sinSqTh*cos2Ph; + } + + inline double VectorSDMEDecay::W2spin0(double sin2Th, double sinPh,double sinSqTh, double sin2Ph) const noexcept{ + + return TMath::Sqrt2() * _rho->Im(2,1,0)*sin2Th*sinPh + + _rho->Im(2,1,-1)*sinSqTh*sin2Ph; + } + inline double VectorSDMEDecay::W3spin0(double sin2Th, double sinPh,double sinSqTh, double sin2Ph) const noexcept{ + + return TMath::Sqrt2() * _rho->Im(3,1,0)*sin2Th*sinPh + + _rho->Im(3,1,-1)*sinSqTh*sin2Ph; + } + + inline double VectorSDMEDecay::W0lepto(double cosSqTh,double sin2Th,double cosPh,double sinSqTh,double cos2Ph) const noexcept{ + + return 0.5 * (1 + _rho->Re(0,0,0) ) + - 0.5 * (3*_rho->Re(0,0,0) - 1) *cosSqTh + + TMath::Sqrt2() * _rho->Re(0,1,0)*sin2Th*cosPh + + _rho->Re(0,1,-1)*sinSqTh*cos2Ph; + } + inline double VectorSDMEDecay::W1lepto(double sinSqTh,double cosSqTh,double sin2Th, double cosPh, double cos2Ph) const noexcept{ + + return _rho->Re(1,1,1) * (1+cosSqTh) + + _rho->Re(1,0,0) * sinSqTh + + TMath::Sqrt2() * _rho->Re(1,1,0)*sin2Th*cosPh + + _rho->Re(1,1,-1)*sinSqTh*cos2Ph; + } + inline double VectorSDMEDecay::W2lepto(double sin2Th, double sinPh,double sinSqTh, double sin2Ph) const noexcept{ + + return -TMath::Sqrt2() * _rho->Im(2,1,0)*sin2Th*sinPh + - _rho->Im(2,1,-1)*sinSqTh*sin2Ph; + } + inline double VectorSDMEDecay::W3lepto(double sin2Th, double sinPh,double sinSqTh, double sin2Ph) const noexcept{ + + return -TMath::Sqrt2() * _rho->Im(3,1,0)*sin2Th*sinPh + - _rho->Im(3,1,-1)*sinSqTh*sin2Ph; + } }//namespace elSpectro diff --git a/core/amplitude_blend.cpp b/core/amplitude_blend.cpp index 8ee6971..5368633 100644 --- a/core/amplitude_blend.cpp +++ b/core/amplitude_blend.cpp @@ -12,26 +12,26 @@ std::complex jpacPhoto::amplitude_blend::helicity_amplitude(std::array helicities, double s, double t) { - int index = find_helicity(helicities, _kinematics->_jp[0], _kinematics->_mB); + int index = find_helicity(helicities, _kinematics->get_meson_JP()[0], _kinematics->get_baryon_JP()[0], _kinematics->is_photon()); if(s<_low_max){ - _amp_low->check_cache(s, t); - return _amp_low->_cached_helicity_amplitude[index]; - } + _amp_low->update_cache(s, t); + return _amp_low->get_cached_helicity_amplitude(index); + } else if(s>_high_min){ - _amp_high->check_cache(s, t); - return _amp_high->_cached_helicity_amplitude[index]; + _amp_high->update_cache(s, t); + return _amp_high->get_cached_helicity_amplitude(index); } else{ - _amp_low->check_cache(s, t); - _amp_high->check_cache(s, t); + _amp_low->update_cache(s, t); + _amp_high->update_cache(s, t); double blend_range= _high_min - _low_max; double fraction_low = 1 - (s-_low_max)/blend_range; - auto low_contrib = _amp_low->_cached_helicity_amplitude[index] * fraction_low; - auto high_contrib = _amp_high->_cached_helicity_amplitude[index] * (1 - fraction_low); + auto low_contrib = _amp_low->get_cached_helicity_amplitude(index) * fraction_low; + auto high_contrib = _amp_high->get_cached_helicity_amplitude(index) * (1 - fraction_low); return low_contrib + high_contrib ; } diff --git a/core/amplitude_blend.hpp b/core/amplitude_blend.hpp index e060f55..7bbcd2a 100644 --- a/core/amplitude_blend.hpp +++ b/core/amplitude_blend.hpp @@ -7,7 +7,7 @@ // --------------------------------------------------------------------------- #pragma once -#include "amplitudes/amplitude.hpp" +#include "core/amplitude.hpp" // --------------------------------------------------------------------------- // The amplitude_blend class can take a vector of the above amplitude objects @@ -41,16 +41,24 @@ namespace jpacPhoto }; //assume low amp allowedJP, - inline std::vector> allowedJP() - { - return _amp_low->allowedJP(); - }; + // inline std::vector> allowedJP() + // { + // //return _amp_low->allowedJP(); + // return _amp_low->allowed_meson_JP(); + // }; + + //take low amp allowed JP + std::vector> allowed_meson_JP() override{ return _amp_low->allowed_meson_JP();} + std::vector> allowed_baryon_JP() override{ + return _amp_low->allowed_baryon_JP(); + } + helicity_channel helicity_CM_frame() override {return _amp_low->helicity_CM_frame();} // TODO: Add a set_params which timesi in one vector and allocates approriaten number of // params to each sub amplitude // Evaluate the sum for given set of helicites, energy, and cos - std::complex helicity_amplitude(std::array helicities, double s, double t); + std::complex helicity_amplitude(std::array helicities, double s, double t) override; }; } diff --git a/examples/CLAS_JPAC_Jpsi.C b/examples/CLAS_JPAC_Jpsi.C new file mode 100644 index 0000000..b77489a --- /dev/null +++ b/examples/CLAS_JPAC_Jpsi.C @@ -0,0 +1,221 @@ +#include "Interface.h" +#include "LorentzVector.h" +#include "DistTF1.h" +#include "FunctionsForElectronScattering.h" +#include "LundWriter.h" +#include "TwoBody_stu.h" + +#include "reaction_kinematics.hpp" +#include "regge_trajectory.hpp" +#include "core/vector_exchange.hpp" +#include "core/pseudoscalar_exchange.hpp" +#include "core/pomeron_exchange.hpp" +#include "core/amplitude_sum.hpp" +#include "core/baryon_resonance.hpp" + +#include +#include +#include +#include + +jpacPhoto::amplitude_sum Amplitude(); + +// --------------------------------------------------------------------------- +// Diagnostic histograms +// --------------------------------------------------------------------------- + +double minMass = 0.2; +double maxMass = 4; +TH1F hQ2("Q2","Q2",1000,0,5); +TH2F hWQ2("WQ2","WQ2",100,4,5,100,0,5); +TH2F hWt("Wt","Wt",100,3.8,4.6,100,0,9); +TH1D heE("eE","eE",1000,0,4); +TH1D heTh("eTh","eTh",1000,0,18); +TH1D hYTh("YTh","YTh",1000,0,180); +TH1F hW("W","W",1000,4,5); +TH1F ht("t","t",1000,0,10); +TH1F hgE("gE","gE",1000,0,20); +TH1F hgTh("gTh","gTh",1000,0,180); +TH1F hXM("MesonM","; X to #pi#pi#pi Mass (GeV)",1000,minMass,maxMass); +TH1F hP1("P1","Pi1",100,0,10); +TH1F hP2("P2","Pi2",100,0,10); +TH1F hPm("Pm","Pi-",100,0,10); +TH1F hVectorCosTh("CosThGJ","cos(#theta_{GJ})",100,-1,1); +TH1F hVectorPhi("PhiGJ","#phi_{GJ}",100,-180,180); +TH1F hScatPhi("Phi_Y","#phi_{Y}",100,-180,180); +TH1F hScatPhidiff("Phi_diff","#phi_{Y}",100,-180*2,180*2); +TH1F hScatPhisum("Phi_sum","#phi_{Y}",100,-180*2,180*2); +TH2F hVectorvScatPhi("PhiGJPhi_Y","#phi_{GJ} v #phi_{Y}",50,-180,180,50,-180,180); +TH2F hVectorPhvVectorTh("ThPhiGJ","#phi_{GJ} v #theta_{GJ}",50,-1,1,50,-180,180); +TH2F hVectorPhvVectorTh2("ThPhiGJ2","#phi_{GJ} v #theta_{GJ}",50,-1,1,50,-180,180); + +void CLAS_JPAC_Jpsi(double ebeamE = 10.6, double nLumi=1E35, int nDays = 80) { + + Double_t crossingAngle=0; //0mrad + //define e- beam, pdg =11 + auto elBeam = initial(11,ebeamE); + auto elBeamP4=elBeam->GetInteracting4Vector(); + elBeam->SetAngleThetaPhi(0,0); + + //define pr beam, pdg =2212 + auto prBeam = initial(2212,0); + auto prBeamP4=prBeam->GetInteracting4Vector(); + + + //OR JPAC amplitude version + // not last argument ==1 =>use lepton decay distributions 0=> use decay to spin 0 + auto jpsi=particle(443,model(new VectorSDMEDecay({},{11,-11},1))); + auto jpac_amp = Amplitude(); + auto pGammaStarDecay = new JpacModelst{&jpac_amp, {jpsi},{2212} }; + auto photoprod = new DecayModelQ2W{0,pGammaStarDecay,new TwoBody_stu{1, .0, 1}}; + + //combine beam, target and reaction products + auto production=eic( elBeam,prBeam,photoprod ); + + + auto Jel = static_cast(jpsi)->Model()->Products()[1]; + auto Jpo = static_cast(jpsi)->Model()->Products()[0]; + + auto proton = pGammaStarDecay->GetBaryon(); + auto electron = dynamic_cast(production->Model())->GetScatteredElectron(); + // --------------------------------------------------------------------------- + // Initialize LUND + // --------------------------------------------------------------------------- + + writer(new EICSimpleWriter{Form("out_clas/ep_to_pJpsi_%d.dat",(int)(ebeamE))}); + + //initilase the generator, may take some time for making distribution tables + initGenerator(); + + // --------------------------------------------------------------------------- + //Set number of events via experimental luminosity and beamtime + // --------------------------------------------------------------------------- + production->SetCombinedBranchingFraction(0.06); //Just Jpsi->e+e- + generator().SetNEvents_via_LuminosityTime(nLumi,nDays*24*60*60); + //or can just do generator().SetNEvents(1E6); + auto fastIntegral=production->IntegrateCrossSectionFast(); + std::cout<<" check fast cross section "<Start("e"); + + while(finishedGenerator()==false){ + nextEvent(); + countGenEvent(); + if(generator().GetNDone()%10000==0) std::cout<<"event number "<P4(); + double W = (photon+ *prBeamP4).M(); + + //fill diagnostic histograms + double Q2 = -photon.M2(); + hWQ2.Fill(W,Q2); + hQ2.Fill(Q2); + hW.Fill(W); + hgE.Fill(photon.E()); + + auto elec = electron->P4(); + heTh.Fill(elec.Theta() *TMath::RadToDeg()); + heE.Fill(elec.E()); + + auto Xrec=Jel->P4()+Jpo->P4(); + hXM.Fill(Xrec.M()); + + hP1.Fill(Jel->P4().P()); + hPm.Fill(Jpo->P4().P()); + + double t = -1*(proton->P4()-*prBeamP4).M2();// + kine::t0(W,0,prbeam.M(),Xrec.M(),prbeam.M()); + ht.Fill(t); + hWt.Fill(W,t); + + //SDME related angles + MomentumVector elScatAngles; + auto gStarN=(photon+*prBeamP4); + kine::electroCMDecay(&gStarN,elBeamP4,&photon,&Xrec,&elScatAngles); + hScatPhi.Fill(elScatAngles.Phi()*TMath::RadToDeg()); + MomentumVector vectorAngles; + kine::mesonDecayGJ(&photon,&Xrec,&(proton->P4()),&Jel->P4(),&vectorAngles); + hVectorCosTh.Fill(TMath::Cos(vectorAngles.Theta())); + hVectorPhi.Fill(vectorAngles.Phi()*TMath::RadToDeg()); + hVectorvScatPhi.Fill(elScatAngles.Phi()*TMath::RadToDeg(),vectorAngles.Phi()*TMath::RadToDeg()); + hVectorPhvVectorTh.Fill(TMath::Cos(vectorAngles.Theta()),vectorAngles.Phi()*TMath::RadToDeg()); + if(W>4.4&&W<4.5) + hVectorPhvVectorTh2.Fill(TMath::Cos(vectorAngles.Theta()),vectorAngles.Phi()*TMath::RadToDeg()); + + + hScatPhidiff.Fill((vectorAngles.Phi()-elScatAngles.Phi())*TMath::RadToDeg()); + hScatPhisum.Fill((vectorAngles.Phi()+elScatAngles.Phi())*TMath::RadToDeg()); + + } + gBenchmark->Stop("e"); + gBenchmark->Print("e"); + + generator().Summary(); + // internally stored histograms for total ep cross section + TH1D *hWdist = (TH1D*)gDirectory->FindObject("Wdistlow"); + TH1D *hGenWdist = (TH1D*)gDirectory->FindObject("genWdist"); + + TFile *fout = TFile::Open(Form("out_clas/ep_to_pJpsi_%d.root",(int)ebeamE), "recreate"); + // total ep cross section inputs + if(hWdist)hWdist->Write(); + if(hGenWdist)hGenWdist->Write(); + + // generated event distributions + hWQ2.Write(); + hWt.Write(); + hQ2.Write(); + hW.Write(); + ht.Write(); + hgE.Write(); + heTh.Write(); + heE.Write(); + hXM.Write(); + fout->Close(); + + + +} + +jpacPhoto::amplitude_sum Amplitude(){ + using namespace jpacPhoto; + + reaction_kinematics * ptr = new reaction_kinematics(3.0969160);//Jpsi mass + ptr->set_meson_JP(1, -1); + + // --------------------------------------------------------------------------- + // S - CHANNEL + + // Two different pentaquarks + // masses and widths from 2015 LHCb paper [2] + auto P_c4450 =new baryon_resonance(ptr, 3, -1, 4.45, 0.040, "P_{c}(4450)"); + P_c4450->set_params({0.01, .7071}); // 2% branching fraction and equal photocouplings + + auto P_c4380 = new baryon_resonance(ptr, 5, +1, 4.38, 0.205, "P_{c}(4380)"); + P_c4380->set_params({0.01, .7071}); // 2% branching fraction and equal photocouplings + + // --------------------------------------------------------------------------- + // T - CHANNEL + + // Set up pomeron trajectory + // Best fit values from [1] + auto alpha =new linear_trajectory(+1, 0.941, 0.364, "pomeron"); + + // Create amplitude with kinematics and trajectory + auto background = new pomeron_exchange(ptr, alpha, false, "Background"); + + // normalization and t-slope + background->set_params({0.379, 0.12}); + + // --------------------------------------------------------------------------- + // SUM + // --------------------------------------------------------------------------- + // Incoherent sum of the s and t channels + amplitude_sum sum5q(ptr, {background}, "5q Sum"); + // amplitude_sum sum5q(ptr, {background, P_c4450}, "5q Sum"); + // amplitude_sum sum10q(ptr, {background, P_c4450, P_c4380}, "10q Sum"); + + return sum5q; +} diff --git a/examples/EIC_JPAC_X3872.C b/examples/EIC_JPAC_X3872.C index 2fd41c8..ad6d353 100644 --- a/examples/EIC_JPAC_X3872.C +++ b/examples/EIC_JPAC_X3872.C @@ -1,11 +1,11 @@ //Just need jpacPhoto headers #include "reaction_kinematics.hpp" #include "regge_trajectory.hpp" -#include "amplitudes/vector_exchange.hpp" -#include "amplitudes/pseudoscalar_exchange.hpp" -#include "amplitudes/pomeron_exchange.hpp" -#include "amplitudes/amplitude_sum.hpp" -#include "amplitudes/baryon_resonance.hpp" +#include "core/vector_exchange.hpp" +#include "core/pseudoscalar_exchange.hpp" +#include "core/pomeron_exchange.hpp" +#include "core/amplitude_sum.hpp" +#include "core/baryon_resonance.hpp" //Amplitude based on $JPACPHOTO/executables/XYZ_Plots/X3872_high.cpp @@ -13,16 +13,24 @@ //To set luminosity and days change last 2 arguments //e.g. for luminosoty 10^33 and 25 days, e- energy 100 and p energy 100 //with high energy paramterisation : -// 'EIC_JPAC_X3872.C("high",100,100,1E33,25)' +// elspectro 'EIC_JPAC_X3872.C(10,100,1E33,25)' // To just run a fixed number of events leave last // argument 0 and nLumi=number of events -// 'EIC_JPAC_X3872.C("high",100,100,1E4)' +// elspectro 'EIC_JPAC_X3872.C(10,100,1E4)' -void EIC_JPAC_X3872(string ampPar="high",double ebeamE = 5, double pbeamE = 41, double nLumi=100, int nDays = 0) { +void EIC_JPAC_X3872(double ebeamE = 5, double pbeamE = 41, double nLumi=100, int nDays = 0) { - LorentzVector elbeam(0,0,-1*ebeamE,escat::E_el(ebeamE)); - LorentzVector prbeam(0,0,pbeamE,escat::E_pr(pbeamE)); + Double_t crossingAngle=0; //30mrad + //define e- beam, pdg =11 + auto elBeam = initial(11,ebeamE); + auto elBeamP4=elBeam->GetInteracting4Vector(); + elBeam->SetAngleThetaPhi(TMath::Pi()-crossingAngle,0); + //define pr beam, pdg =2212 + auto prBeam = initial(2212,pbeamE); + auto prBeamP4=prBeam->GetInteracting4Vector(); + prBeam->SetAngleThetaPhi(crossingAngle,0); + // --------------------------------------------------------------------------- // AMPLITUDES // --------------------------------------------------------------------------- @@ -34,7 +42,7 @@ void EIC_JPAC_X3872(string ampPar="high",double ebeamE = 5, double pbeamE = 41, // X(3872) auto kX = reaction_kinematics{M_X3872}; - kX.set_JP(1, 1); + kX.set_meson_JP(1, 1); // Nucleon couplings and cutoffs double gV_omega = 16., gT_omega = 0.; @@ -51,27 +59,55 @@ void EIC_JPAC_X3872(string ampPar="high",double ebeamE = 5, double pbeamE = 41, // Linear trajectory for the rho auto alpha = linear_trajectory{-1, 0.5, 0.9, "#rho - #omega"}; + // --------------------------------------------------------------------------- + // High or Low - Energy Amplitudes + // --------------------------------------------------------------------------- + ////////////////// + // X(3872) + vector_exchange X_omega_high{&kX, &alpha, "#omega"}; + X_omega_high.set_params({gX_omega, gV_omega, gT_omega}); + X_omega_high.set_formfactor(true, LamOmega); + vector_exchange X_omega_low{&kX, M_OMEGA, "#omega"}; + X_omega_low.set_params({gX_omega, gV_omega, gT_omega}); + X_omega_low.set_formfactor(true, LamOmega); + + vector_exchange X_rho_high{&kX, &alpha, "#rho"}; + X_rho_high.set_params({gX_rho, gV_rho, gT_rho}); + X_rho_high.set_formfactor(true, LamRho); + vector_exchange X_rho_low{&kX, M_RHO, "#rho"}; + X_rho_low.set_params({gX_rho, gV_rho, gT_rho}); + X_rho_low.set_formfactor(true, LamRho); + + std::vector X_exchanges_low = {&X_omega_low, &X_rho_low}; + amplitude_sum jpac_amp_low(&kX, X_exchanges_low, "#it{X}(3872) low"); + std::vector X_exchanges_high = {&X_omega_high, &X_rho_high}; + amplitude_sum jpac_amp_high(&kX, X_exchanges_high, "#it{X}(3872) high"); + + double low_s = 7*7.;//GeV + double high_s = 20*20.;//GeV + amplitude_blend jpac_amp{&kX,&jpac_amp_low,low_s,&jpac_amp_high,high_s, "#it{X} blend"}; + // --------------------------------------------------------------------------- // High or Low - Energy Amplitudes // --------------------------------------------------------------------------- ////////////////// // X(3872) - vector_exchange *X_omega{nullptr}; - if(ampPar=="high")X_omega= new vector_exchange(&kX, &alpha, "#omega"); - else if(ampPar=="low") X_omega=new vector_exchange(&kX, M_OMEGA, "#omega"); - else {cerr<<"invalid amplitude parameterisation "<set_params({gX_omega, gV_omega, gT_omega}); - X_omega->set_formfactor(true, LamOmega); - - vector_exchange *X_rho{nullptr}; - if(ampPar=="high")X_rho= new vector_exchange(&kX, &alpha, "#rho"); - else if(ampPar=="low") X_rho=new vector_exchange(&kX, M_RHO, "#rho"); - else {cerr<<"invalid amplitude parameterisation "<set_params({gX_rho, gV_rho, gT_rho}); - X_rho->set_formfactor(true, LamRho); - - std::vector X_exchanges = {X_omega, X_rho}; - amplitude_sum jpac_amp(&kX, X_exchanges, "#it{X}(3872)"); + // vector_exchange *X_omega{nullptr}; + // if(ampPar=="high")X_omega= new vector_exchange(&kX, &alpha, "#omega"); + // else if(ampPar=="low") X_omega=new vector_exchange(&kX, M_OMEGA, "#omega"); + // else {cerr<<"invalid amplitude parameterisation "<set_params({gX_omega, gV_omega, gT_omega}); + // X_omega->set_formfactor(true, LamOmega); + + // vector_exchange *X_rho{nullptr}; + // if(ampPar=="high")X_rho= new vector_exchange(&kX, &alpha, "#rho"); + // else if(ampPar=="low") X_rho=new vector_exchange(&kX, M_RHO, "#rho"); + // else {cerr<<"invalid amplitude parameterisation "<set_params({gX_rho, gV_rho, gT_rho}); + // X_rho->set_formfactor(true, LamRho); + + // std::vector X_exchanges = {X_omega, X_rho}; + // amplitude_sum jpac_amp(&kX, X_exchanges, "#it{X}(3872)"); // --------------------------------------------------------------------------- @@ -81,7 +117,7 @@ void EIC_JPAC_X3872(string ampPar="high",double ebeamE = 5, double pbeamE = 41, //create a X decaying to J/psi pi+pi- auto jpsi=particle(443,model(new PhaseSpaceDecay({},{11,-11}))); //rho - mass_distribution(113,new DistTF1{TF1("hhRho","TMath::BreitWigner(x,0.775,0.151)",0.2,0.7)}); + mass_distribution(113,new DistTF1{TF1("hhRho","TMath::BreitWigner(x,0.775,0.151)",0.25,1)}); auto rho=particle(113,model(new PhaseSpaceDecay({},{211,-211}))); //x mass_distribution(9995,new DistTF1{TF1("hh","TMath::BreitWigner(x,3.872,0.001)",3.85,3.89)}); @@ -92,8 +128,8 @@ void EIC_JPAC_X3872(string ampPar="high",double ebeamE = 5, double pbeamE = 41, auto pGammaStarDecay = JpacModelst{&jpac_amp, {x},{2212} }; //photo-nucleon system //Decay g*p state, provide s channel and t-channel "shapes" //Note the amplitude will provide the actual t-distribution, this approximation speeds up sampling - //TwoBody_stu{0., 1.0, 2.5} => 0% s-schannel, 100% t channel with slope 2.5 - auto photoprod = DecayModelQ2W{0,&pGammaStarDecay,new TwoBody_stu{0., 1.0, 2.5}}; + //TwoBody_stu{0., 1.0, 2.5} => 0% s-schannel, 100% t channel with slope 2. + auto photoprod = DecayModelQ2W{0,&pGammaStarDecay,new TwoBody_stu{0., 2.0, 1.0}}; //combine beam, target and reaction products auto production=eic( ebeamE, pbeamE, &photoprod ); @@ -101,7 +137,7 @@ void EIC_JPAC_X3872(string ampPar="high",double ebeamE = 5, double pbeamE = 41, // --------------------------------------------------------------------------- // Initialize HepMC3 // --------------------------------------------------------------------------- - writer(new HepMC3Writer{Form("out/jpac_x3872_%s_%d_%d.txt",ampPar.data(),(int)ebeamE,(int)pbeamE)}); + writer(new HepMC3Writer{Form("out/jpac_x3872_%d_%d.txt",(int)ebeamE,(int)pbeamE)}); // --------------------------------------------------------------------------- @@ -139,7 +175,5 @@ void EIC_JPAC_X3872(string ampPar="high",double ebeamE = 5, double pbeamE = 41, generator().Summary(); - delete X_rho; - delete X_omega; } diff --git a/examples/EIC_JPAC_Y.C b/examples/EIC_JPAC_Y.C new file mode 100644 index 0000000..5cb023d --- /dev/null +++ b/examples/EIC_JPAC_Y.C @@ -0,0 +1,238 @@ +//Just need jpacPhoto headers +#include "reaction_kinematics.hpp" +#include "regge_trajectory.hpp" +#include "core/vector_exchange.hpp" +#include "core/pseudoscalar_exchange.hpp" +#include "core/pomeron_exchange.hpp" +#include "core/amplitude_sum.hpp" +#include "core/baryon_resonance.hpp" + +#include "FunctionsForGenvector.h" + +// --------------------------------------------------------------------------- +// Diagnostic histograms +// --------------------------------------------------------------------------- + +TH1F hQ2("Q2","Q2",1000,0,100); +TH1D hgPhi("gPhi","gPhi",100,-180,180); +TH1D heE("eE","eE",1000,0,50); +TH2D heThPhi("eThPhi","eThPhi",50,-180,180,500,178,180); +TH1D heTh("eTh","eTh",1000,0,180); +TH1D hPTh("NTh","NTh",1000,0,180); +TH2D hPThPhi("PThPhi","PThPhi",50,-180,180,50,0,10); +TH1D hPPhi("NPh","NPh",1000,-180,180); +TH1D hePhi("ePh","ePh",90,-180,180); +TH1F hW("W","W",1000,0,100); +TH1F ht("t","t",1000,0,10); +TH1F hgE("gE","gE",1000,0,20); +TH1F hMesonM("MesonM","; J/#psi#pi Mass (GeV)",1000,3.5,5); +TH1F h2PiM("M2Pi",";#pi+#pi- Mass (GeV)",1000,0,2); +TH1F hJpsiM("JpsiM","; e+e- Mass (GeV)",1000,2.9,3.3); +TH2F hElePVsEta("ElePVsEta","; #eta; p (GeV)",200,-5,5,500,0,50); +TH2F hPosPVsEta("PosPVsEta","; #eta; p (GeV)",200,-5,5,500,0,50); +TH2F hPionPVsEta("PionPVsEta","; #eta; p (GeV)",200,-5,5,500,0,50); +TH2F hRecoilPVsEta("RecoilPVsEta","; #eta; p (GeV)",200,0,10,1000,0,275); +TH2F hRecoilThetaVsP("RecoilThetaVsP","; p (GeV); #theta (mrad)",1000,0,275,200,0,200); +TH1F hRecoilPt("RecoilPt","; p_{T} (GeV)",200,0,5.0); + + +void EIC_JPAC_Y(double ebeamE = 10, double pbeamE = 100, double nLumi=1E33, int nDays = 100) { + + Double_t crossingAngle=0.0; //0mrad + //define e- beam, pdg =11 + auto elBeam = initial(11,ebeamE); + auto elBeamP4=elBeam->GetInteracting4Vector(); + elBeam->SetAngleThetaPhi(TMath::Pi()-crossingAngle,0); + + //define pr beam, pdg =2212 + auto prBeam = initial(2212,pbeamE); + auto prBeamP4=prBeam->GetInteracting4Vector(); + prBeam->SetAngleThetaPhi(crossingAngle,0); + + // --------------------------------------------------------------------------- + // AMPLITUDES + // --------------------------------------------------------------------------- + + // --------------------------------------------------------------------------- + // Preliminaries + // --------------------------------------------------------------------------- + + + + // Y(4260) + auto kY=reaction_kinematics(M_Y4260); + kY.set_meson_JP(1, -1); + // double R_Y = 1.55;//changed in JpacPhoto 11/10/2022 + double R_Y = 0.84; + + + // Same but high-energy + linear_trajectory alpha_traj_high{1, 1.15, 0.11, "HE"}; + double b_HE = 1.01; + double A_HE = 0.16; + pomeron_exchange Y_Amp_high{&kY, &alpha_traj_high, true, "#it{Y}(4260) high"}; + Y_Amp_high.set_params({A_HE * R_Y, b_HE}); + + // Low - energy trajectory and couplings + linear_trajectory alpha_traj_low{1, 0.94, 0.36, "LE"}; + double b_LE = 0.12; + double A_LE = 0.38; + pomeron_exchange Y_Amp_low{&kY, &alpha_traj_low, false, "#it{Y}(4260) low"}; + Y_Amp_low.set_params({A_LE * R_Y, b_LE}); + + double low_s = 7*7.;//GeV + double high_s = 20*20.;//GeV + amplitude_blend jpac_amp{&kY,&Y_Amp_low,low_s,&Y_Amp_high,high_s, "#it{X} blend"}; + + + // ------------------------------------------------------------------- + // elSpectro + // -------------------------------------------------------------------- + //create a V decaying to J/psi pi+pi- + auto jpsi=particle(443,model(new PhaseSpaceDecay({},{11,-11}))); + + //flat 2pi mass distribution + mass_distribution(9996,new DistTF1{TF1("hhsigma","1",0.25,4)}); + auto sigma=particle(9996,model(new PhaseSpaceDecay({},{211,-211}))); + + mass_distribution(9995,new DistTF1{TF1("hh","TMath::BreitWigner(x,4.22,0.05)",3.5,5)}); + auto Y=particle(9995,model(new PhaseSpaceDecay{{jpsi,sigma},{}})); + Y->SetPdgMass(M_Y4260); + + //create eic electroproduction of X + proton + auto pGammaStarDecay = JpacModelst{&jpac_amp, {Y},{2212} }; //photo-nucleon system + //Decay g*p state, provide s channel and t-channel "shapes" + //Note the amplitude will provide the actual t-distribution, this approximation speeds up sampling + //TwoBody_stu{0., 1.0, 2.5} => 0% s-schannel, 100% t channel with slope 2.5 + auto photoprod = DecayModelQ2W{0,&pGammaStarDecay,new TwoBody_stu{0., 1.0, 3}}; + //combine beam, target and reaction products + //auto production=eic( ebeamE, pbeamE, &photoprod ); + auto production=eic( elBeam,prBeam,&photoprod); + + // --------------------------------------------------------------------------- + // Initialize HepMC3 + // --------------------------------------------------------------------------- + writer(new HepMC3Writer{Form("outCollision/jpac_Y_%d_%d.txt",(int)ebeamE,(int)pbeamE)}); + // exit(0); + + // --------------------------------------------------------------------------- + //initilase the generator, may take some time for making distribution tables + // --------------------------------------------------------------------------- + initGenerator(); + + // --------------------------------------------------------------------------- + //Set number of events via experimental luminosity and beamtime + // --------------------------------------------------------------------------- + production->SetCombinedBranchingFraction(0.06*0.01); //Jpsi->e+e- and 1% Y + generator().SetNEvents_via_LuminosityTime(nLumi,24*60*60*nDays); + //generator().SetNEvents(10000); + // auto fastIntegral=production->IntegrateCrossSectionFast(); + //std::cout<<" check fast cross section "<(sigma)->Model()->Product(0); + auto pionm = static_cast(sigma)->Model()->Product(1); + auto posJ = static_cast(jpsi)->Model()->Products()[0]; + auto eleJ = static_cast(jpsi)->Model()->Products()[1]; + + auto electron = photoprod.GetScatteredElectron(); + auto neutron = photoprod.GetDecayBaryon(); + + // --------------------------------------------------------------------------- + // Generate events + // --------------------------------------------------------------------------- + + gBenchmark->Start("generator");//timer + + while(finishedGenerator()==false){ + nextEvent(); + countGenEvent(); + if(generator().GetNDone()%1000==0) std::cout<<"event number "<P4(); + double Q2 = -photon.M2(); + double W = (photon+*prBeamP4).M(); + double t = -1*(neutron->P4()- *prBeamP4).M2(); + hQ2.Fill(Q2); + hW.Fill(W); + ht.Fill(t); + hgE.Fill(photon.E()); + + auto forScPhi= electron->P4(); + hgPhi.Fill(forScPhi.Phi()*TMath::RadToDeg()); + + + auto elec = electron->P4(); + heThPhi.Fill(forScPhi.Phi() *TMath::RadToDeg(),forScPhi.Theta() *TMath::RadToDeg()); + heTh.Fill(elec.Theta() *TMath::RadToDeg()); + hePhi.Fill(elec.Phi() *TMath::RadToDeg()); + heE.Fill(elec.E()); + + auto jpsiP4 = eleJ->P4() + posJ->P4(); + hJpsiM.Fill(jpsiP4.M()); + auto sigmaP4 = pionp->P4() + pionm->P4(); + auto jpsi_pion=sigmaP4+jpsiP4; + hMesonM.Fill(jpsi_pion.M()); + h2PiM.Fill(sigmaP4.M()); + + hElePVsEta.Fill(eleJ->P4().Eta(), eleJ->P4().P()); + hPosPVsEta.Fill(posJ->P4().Eta(), posJ->P4().P()); + hPionPVsEta.Fill(pionp->P4().Eta(), pionp->P4().P()); + hRecoilPVsEta.Fill(neutron->P4().Eta(), neutron->P4().P()); + hRecoilThetaVsP.Fill(neutron->P4().P(), neutron->P4().Theta()*1000.); + hRecoilPt.Fill(ROOT::Math::VectorUtil::Perp(prBeamP4->Vect(),neutron->P4().Vect())); + + hPTh.Fill(neutron->P4().Theta()*TMath::RadToDeg()); + hPThPhi.Fill(neutron->P4().Phi()*TMath::RadToDeg(),neutron->P4().Theta()*TMath::RadToDeg()); + hPPhi.Fill(neutron->P4().Phi()*TMath::RadToDeg()); + + + } + + gBenchmark->Stop("generator"); + gBenchmark->Print("generator"); + + // --------------------------------------------------------------------------- + // Report generator statistics, can be used for optimising + // --------------------------------------------------------------------------- + + generator().Summary(); + + + // --------------------------------------------------------------------------- + // Write diagnostic histograms + // --------------------------------------------------------------------------- + TH1D *hWdist = (TH1D*)gDirectory->FindObject("Wdist"); + TH1D *hGenWdist = (TH1D*)gDirectory->FindObject("genWdist"); + + TFile *fout = TFile::Open(Form("out/jpac_y_%d_%d.txt",(int)ebeamE,(int)pbeamE), "recreate"); + // generated event distributions + hQ2.Write(); + hW.Write(); + ht.Write(); + hgE.Write(); + hgPhi.Write(); + heTh.Write(); + hPTh.Write(); + hPPhi.Write(); + hePhi.Write(); + heE.Write(); + hJpsiM.Write(); + hMesonM.Write(); + h2PiM.Write(); + hElePVsEta.Write(); + hPosPVsEta.Write(); + hPionPVsEta.Write(); + hRecoilPVsEta.Write(); + hRecoilThetaVsP.Write(); + hRecoilPt.Write(); + hPThPhi.Write(); + heThPhi.Write(); + hWdist->Write(); + fout->Close(); + +} diff --git a/examples/EIC_JPAC_nZc.C b/examples/EIC_JPAC_nZc.C index 6d29d3f..8bf4dcd 100644 --- a/examples/EIC_JPAC_nZc.C +++ b/examples/EIC_JPAC_nZc.C @@ -1,42 +1,83 @@ //Just need jpacPhoto headers #include "reaction_kinematics.hpp" #include "regge_trajectory.hpp" -#include "amplitudes/vector_exchange.hpp" -#include "amplitudes/pseudoscalar_exchange.hpp" -#include "amplitudes/pomeron_exchange.hpp" -#include "amplitudes/amplitude_sum.hpp" -#include "amplitudes/baryon_resonance.hpp" +#include "core/vector_exchange.hpp" +#include "core/pseudoscalar_exchange.hpp" +#include "core/pomeron_exchange.hpp" +#include "core/amplitude_sum.hpp" +#include "core/baryon_resonance.hpp" +#include "FunctionsForGenvector.h" + +// --------------------------------------------------------------------------- +// Diagnostic histograms +// --------------------------------------------------------------------------- + +TH1F hQ2("Q2","Q2",1000,0,100); +TH1D hgPhi("gPhi","gPhi",100,-180,180); +TH1D heE("eE","eE",1000,0,50); +TH2D heThPhi("eThPhi","eThPhi",50,-180,180,500,178,180); +TH1D heTh("eTh","eTh",1000,0,180); +TH1D hPTh("NTh","NTh",1000,0,180); +TH2D hPThPhi("PThPhi","PThPhi",50,-180,180,50,0,90); +TH2D hJeThPhi("JeThPhi","JEThPhi",50,-180,180,50,0,90); +TH1D hPPhi("NPh","NPh",1000,-180,180); +TH1D hePhi("ePh","ePh",90,-180,180); +TH1D hZPhi("ZPh","ZPh",90,-180,180); +TH1D hJPhi("JPh","JPh",90,-180,180); +TH1D hJePhi("JePh","JePh",90,-180,180); +TH1F hW("W","W",1000,0,50); +TH1F ht("t","t",1000,0,10); +TH1F hgE("gE","gE",1000,0,20); +TH1F hMesonM("MesonM","; J/#psi#pi Mass (GeV)",1000,3.5,4.3); +TH1F hJpsiM("JpsiM","; e+e- Mass (GeV)",1000,2.9,3.3); +TH2F hElePVsEta("ElePVsEta","; #eta; p (GeV)",200,-5,5,500,0,50); +TH2F hPosPVsEta("PosPVsEta","; #eta; p (GeV)",200,-5,5,500,0,50); +TH2F hPionPVsEta("PionPVsEta","; #eta; p (GeV)",200,-5,5,500,0,50); +TH2F hRecoilPVsEta("RecoilPVsEta","; #eta; p (GeV)",200,0,10,1000,0,275); +TH2F hRecoilThetaVsP("RecoilThetaVsP","; p (GeV); #theta (mrad)",1000,0,275,200,0,200); +TH1F hRecoilPt("RecoilPt","; p_{T} (GeV)",200,0,5.0); + +TH1F hScatPhi("Phi_Y","#phi_{Y}",100,-180,180); +TH1F hVectorCosTh("CosThGJ","cos(#theta_{GJ})",100,-1,1); +TH1F hVectorPhi("PhiGJ","#phi_{GJ}",100,-180,180); +TH2F hVectorvScatPhi("PhiGJPhi_Y","#phi_{GJ} v #phi_{Y}",50,-180,180,50,-180,180); //Amplitude based on $JPACPHOTO/executables/XYZ_Plots/Z_low.cpp and Z_high.cpp //To set luminosity and days change last 2 arguments //e.g. for luminosoty 10^33 and 25 days, e- energy 100 and p energy 100 //with high energy paramterisation : -// 'EIC_JPAC_X3872.C("high",100,100,1E33,25)' +//elspectro --i 'EIC_JPAC_nZc.C(10,100,1E33,25)' // To just run a fixed number of events leave last // argument 0 and nLumi=number of events -// 'EIC_JPAC_X3872.C("high",100,100,1E4)' +// elspectro --i 'EIC_JPAC_nZc.C(10,100,1E4)' -void EIC_JPAC_nZc(string ampPar="high",double ebeamE = 5, double pbeamE = 41, double nLumi=100, int nDays = 0) { +void EIC_JPAC_nZc(double ebeamE = 5, double pbeamE = 100, double nLumi=1E34, int nDays = 50) { - LorentzVector elbeam(0,0,-1*ebeamE,escat::E_el(ebeamE)); - LorentzVector prbeam(0,0,pbeamE,escat::E_pr(pbeamE)); + Double_t crossingAngle=0; + //define e- beam, pdg =11 + auto elBeam = initial(11,ebeamE); + auto elBeamP4=elBeam->GetInteracting4Vector(); + elBeam->SetAngleThetaPhi(TMath::Pi()-crossingAngle,0); + //define pr beam, pdg =2212 + auto prBeam = initial(2212,pbeamE); + auto prBeamP4=prBeam->GetInteracting4Vector(); + prBeam->SetAngleThetaPhi(crossingAngle,0); + // --------------------------------------------------------------------------- // AMPLITUDES // --------------------------------------------------------------------------- - // --------------------------------------------------------------------------- - // Preliminaries - // --------------------------------------------------------------------------- - double g_NN = sqrt(4. * M_PI * 13.81); // Nucleon coupling same for all + double g_NN = sqrt(2.) * sqrt(4. * M_PI * 13.81); // Nucleon coupling same for all double LamPi = .9; // 900 MeV cutoff for formfactor double bPi = 1. / (LamPi * LamPi); // Zc(3900) - auto kZc = reaction_kinematics{M_ZC3900}; - kZc.set_JP(1, 1); + double mZc = 3.8884; // GeV + auto kZc = reaction_kinematics{mZc}; + kZc.set_meson_JP(1, +1); double gc_Psi = 1.91; // psi coupling before VMD scaling double gc_Gamma = E * F_JPSI * gc_Psi / M_JPSI; @@ -48,21 +89,24 @@ void EIC_JPAC_nZc(string ampPar="high",double ebeamE = 5, double pbeamE = 41, do double alpha_0 = - alpha_prime * M2_PION; auto alpha = linear_trajectory{signature, alpha_0, alpha_prime}; - pseudoscalar_exchange* ampZc{nullptr}; // --------------------------------------------------------------------------- // low => Fixed-spin amplitudes // --------------------------------------------------------------------------- - if(ampPar=="low") ampZc = new pseudoscalar_exchange{&kZc, M_PION, "#it{Z_{c}} (3900)^{+}"}; + pseudoscalar_exchange ampZcLow{&kZc, M_PION, "#it{Z_{c}} (3900)^{+}"}; + ampZcLow.set_params(Zc_couplings); + ampZcLow.set_formfactor(true, bPi); // --------------------------------------------------------------------------- // high => Reggeized amplitudes // --------------------------------------------------------------------------- - else if(ampPar=="high") ampZc = new pseudoscalar_exchange(&kZc, &alpha, "#it{Z_{c}}(3900)^{+}"); - else {cerr<<"invalid amplitude parameterisation "<set_params(Zc_couplings); - ampZc->set_formfactor(true, bPi); + pseudoscalar_exchange ampZcHigh(&kZc, &alpha, "#it{Z_{c}}(3900)^{+}"); + ampZcHigh.set_params(Zc_couplings); + ampZcHigh.set_formfactor(true, bPi); + double low_s = 15*15.;//GeV + double high_s = 20*20.;//GeV + amplitude_blend ampZc{&kZc,&ZcLow,low_s,&ZcHigh,high_s, "#it{Z_{c}}(3900)^{+} blend"}; + //auto ampZc = ampZcHigh; // --------------------------------------------------------------------------- // elSpectro // --------------------------------------------------------------------------- @@ -78,21 +122,20 @@ void EIC_JPAC_nZc(string ampPar="high",double ebeamE = 5, double pbeamE = 41, do Z->SetPdgMass(M_ZC3900); //create eic electroproduction of X + proton - auto pGammaStarDecay = JpacModelst{ampZc, {Z},{2112} }; //photo-nucleon system + auto pGammaStarDecay = JpacModelst{&Zc, {Z},{2112} }; //photo-nucleon system //auto pGammaStarDecay = DecayModelst{ {Z},{2212} }; //photo-nucleon system //Decay g*p state, provide s channel and t-channel "shapes" //Note the amplitude will provide the actual t-distribution, this approximation speeds up sampling //TwoBody_stu{0., 1.0, 2.5} => 0% s-schannel, 100% t channel with slope 2.5 - auto photoprod = DecayModelQ2W{0,&pGammaStarDecay,new TwoBody_stu{0., 1.0, 2.5}}; + auto photoprod = DecayModelQ2W{0,&pGammaStarDecay,new TwoBody_stu{0., 1.0, 1.5}}; //combine beam, target and reaction products - auto production=eic( ebeamE, pbeamE, &photoprod ); - + auto production=eic( elBeam,prBeam,&photoprod ); // --------------------------------------------------------------------------- // Initialize HepMC3 // --------------------------------------------------------------------------- - writer(new HepMC3Writer{Form("out/jpac_Zc3900_%s_%d_%d.txt",ampPar.data(),(int)ebeamE,(int)pbeamE)}); - + writer(new EICSimpleWriter{Form("outCollision/jpac_Zc3900_%d_%d.txt",(int)ebeamE,(int)pbeamE)}); + // exit(0); // --------------------------------------------------------------------------- //initilase the generator, may take some time for making distribution tables @@ -102,12 +145,23 @@ void EIC_JPAC_nZc(string ampPar="high",double ebeamE = 5, double pbeamE = 41, do // --------------------------------------------------------------------------- //Set number of events via experimental luminosity and beamtime // --------------------------------------------------------------------------- - production->SetCombinedBranchingFraction(0.06); //Just Jpsi->e+e- + production->SetCombinedBranchingFraction(0.06*0.1); //Just Jpsi->e+e- generator().SetNEvents_via_LuminosityTime(nLumi,24*60*60*nDays); - - //or can just do generator().SetNEvents(1E6); + // generator().SetNEvents(10000); auto fastIntegral=production->IntegrateCrossSectionFast(); std::cout<<" check fast cross section "<(Z)->Model()->Products()[1]; + auto posJ = static_cast(jpsi)->Model()->Products()[0]; + auto eleJ = static_cast(jpsi)->Model()->Products()[1]; + + auto electron = photoprod.GetScatteredElectron(); + auto neutron = photoprod.GetDecayBaryon(); + ROOT::Math::RotationZYX rotateToZaxis(-elBeamP4->Phi(),-elBeamP4->Theta(),0); // --------------------------------------------------------------------------- // Generate events // --------------------------------------------------------------------------- @@ -118,6 +172,68 @@ void EIC_JPAC_nZc(string ampPar="high",double ebeamE = 5, double pbeamE = 41, do nextEvent(); countGenEvent(); if(generator().GetNDone()%1000==0) std::cout<<"event number "<P4(); + double Q2 = -photon.M2(); + double W = (photon+*prBeamP4).M(); + double t = -1*(neutron->P4()- *prBeamP4).M2(); + hQ2.Fill(Q2); + hW.Fill(W); + ht.Fill(t); + hgE.Fill(photon.E()); + + auto forScPhi= electron->P4(); + forScPhi=rotateToZaxis*forScPhi; + hgPhi.Fill(forScPhi.Phi()*TMath::RadToDeg()); + + LorentzVector totalIn = *elBeamP4 + *prBeamP4; + LorentzVector totalOut = electron->P4() + neutron->P4() + eleJ->P4() + posJ->P4() + pionp->P4(); + // cout<<"CHECK ENERGY CONSERVATION "<P4(); + // heThPhi.Fill(elec.Phi() *TMath::RadToDeg(),elec.Theta() *TMath::RadToDeg()); + heThPhi.Fill(forScPhi.Phi() *TMath::RadToDeg(),forScPhi.Theta() *TMath::RadToDeg()); + heTh.Fill(elec.Theta() *TMath::RadToDeg()); + hePhi.Fill(elec.Phi() *TMath::RadToDeg()); + heE.Fill(elec.E()); + + auto jpsiP4 = eleJ->P4() + posJ->P4(); + hJpsiM.Fill(jpsiP4.M()); + auto jpsi_pion=pionp->P4()+jpsiP4; + hMesonM.Fill(jpsi_pion.M()); + //ROOT::Math::RotationZYX _rotateToZaxis(poto); + //SDME related angles + MomentumVector elScatAngles; + auto gStarN=(photon+*prBeamP4); + kine::electroCMDecay(&gStarN,elBeamP4,&photon,&jpsi_pion,&elScatAngles); + hScatPhi.Fill(elScatAngles.Phi()*TMath::RadToDeg()); + MomentumVector vectorAngles; + // std::cout<<"***************************************** in script th = "<P4(),&vectorAngles); + hVectorCosTh.Fill(TMath::Cos(vectorAngles.Theta())); + hVectorPhi.Fill(vectorAngles.Phi()*TMath::RadToDeg()); + hVectorvScatPhi.Fill(elScatAngles.Phi()*TMath::RadToDeg(),vectorAngles.Phi()*TMath::RadToDeg()); + + hElePVsEta.Fill(eleJ->P4().Eta(), eleJ->P4().P()); + hPosPVsEta.Fill(posJ->P4().Eta(), posJ->P4().P()); + hPionPVsEta.Fill(pionp->P4().Eta(), pionp->P4().P()); + hRecoilPVsEta.Fill(neutron->P4().Eta(), neutron->P4().P()); + hRecoilThetaVsP.Fill(neutron->P4().P(), neutron->P4().Theta()*1000.); + hRecoilPt.Fill(ROOT::Math::VectorUtil::Perp(prBeamP4->Vect(),neutron->P4().Vect())); + + hJPhi.Fill((jpsiP4.Phi())*TMath::RadToDeg()); + hJePhi.Fill((eleJ->P4().Phi())*TMath::RadToDeg()); + hZPhi.Fill((jpsi_pion.Phi())*TMath::RadToDeg()); + hJeThPhi.Fill(eleJ->P4().Phi()*TMath::RadToDeg(),eleJ->P4().Theta()*TMath::RadToDeg()); + + hPTh.Fill(neutron->P4().Theta()*TMath::RadToDeg()); + hPThPhi.Fill(neutron->P4().Phi()*TMath::RadToDeg(),neutron->P4().Theta()*TMath::RadToDeg()); + hPPhi.Fill(neutron->P4().Phi()*TMath::RadToDeg()); + + } gBenchmark->Stop("generator"); @@ -129,6 +245,36 @@ void EIC_JPAC_nZc(string ampPar="high",double ebeamE = 5, double pbeamE = 41, do generator().Summary(); - delete ampZc; - + // delete ampZc; + + // --------------------------------------------------------------------------- + // Write diagnostic histograms + // --------------------------------------------------------------------------- + TH1D *hWdist = (TH1D*)gDirectory->FindObject("Wdist"); + + TFile *fout = TFile::Open(Form("outCollision/jpac_Zc3900_%d_%d.root",(int)ebeamE,(int)pbeamE), "recreate"); + // generated event distributions + hQ2.Write(); + hW.Write(); + ht.Write(); + hgE.Write(); + hgPhi.Write(); + heTh.Write(); + hPTh.Write(); + hPPhi.Write(); + hePhi.Write(); + heE.Write(); + hJpsiM.Write(); + hMesonM.Write(); + hElePVsEta.Write(); + hPosPVsEta.Write(); + hPionPVsEta.Write(); + hRecoilPVsEta.Write(); + hRecoilThetaVsP.Write(); + hRecoilPt.Write(); + hPThPhi.Write(); + heThPhi.Write(); + hWdist->Write(); + fout->Close(); + } diff --git a/examples/PhotoD_npDiLepton.C b/examples/PhotoD_npDiLepton.C new file mode 100644 index 0000000..5f13a31 --- /dev/null +++ b/examples/PhotoD_npDiLepton.C @@ -0,0 +1,129 @@ +#include "Interface.h" +#include "LorentzVector.h" +#include "DistTF1.h" +#include "FunctionsForElectronScattering.h" +#include "LundWriter.h" +#include "TwoBody_stu.h" + + +#include +#include +#include +#include + +// --------------------------------------------------------------------------- +// Diagnostic histograms +// --------------------------------------------------------------------------- + +double minMass = 0.2; +double maxMass = 3; +TH1F hW("W","W",1000,1,5); +TH1F ht("t","t",1000,0,10); +TH1F hgE("gE","gE",1000,0,20); +TH1F hpP("pP","proton momentum",100,0,10); +TH1F hpTheta("pTh","proton #theta",100,0,180); +TH1F hnP("nP","neutron momentum",100,0,10); +TH1F hnTheta("nTh","neutron #theta",100,0,180); +TH1F heP("eP","e- momentum",100,0,10); +TH1F heTheta("eTh","pe #theta",100,0,180); +TH1F hpoP("poP","proton momentum",100,0,10); +TH1F hpoTheta("poTh","proton #theta",100,0,180); + +void PhotoD_npDiLepton(double ebeamE=12,int nEvents = 10000) { + + // using namespace elSpectro; + // elSpectro::Manager::Instance(); + + auto bremPhoton = initial(22,0,11, + model(new Bremsstrahlung()), + new BremstrPhoton(ebeamE,0.1*ebeamE,ebeamE*0.999)); + //get 4-momentum of brems. photon + auto photonP4=bremPhoton->GetInteracting4Vector(); + + //deuteron target at rest + auto dTarget= initial(45,0); + auto dP4=dTarget->GetInteracting4Vector(); + + auto npee = dynamic_cast(model(new DecayModelDnpee{})); + auto production = photoprod( bremPhoton,dTarget,(new DecayModelW{0, npee}) ); + + + auto proton = npee->GetProton(); + auto neutron = npee->GetNeutron(); + auto ele = npee->GetElectron(); + auto pos = npee->GetPositron(); + + + + // --------------------------------------------------------------------------- + // Initialize LUND with 10000 events per file + // --------------------------------------------------------------------------- + + // writer(new LundWriter{Form("out_photo/gd_to_npDiLept_%d.dat",(int)ebeamE),10000}); + + //initilase the generator, may take some time for making distribution tables + initGenerator(); + + // --------------------------------------------------------------------------- + // Generate events + // --------------------------------------------------------------------------- + + gBenchmark->Start("e"); + + for(int i=0;iPdg()<<" "<P4().P()<<" "<Mass()<P4()+neutron->P4()+ele->P4()+pos->P4()<E()); + + hpP.Fill(proton->P4().P()); + hpTheta.Fill(proton->P4().Theta()*TMath::RadToDeg()); + hnP.Fill(neutron->P4().P()); + hnTheta.Fill(neutron->P4().Theta()*TMath::RadToDeg()); + heP.Fill(ele->P4().P()); + heTheta.Fill(ele->P4().Theta()*TMath::RadToDeg()); + hpoP.Fill(pos->P4().P()); + hpoTheta.Fill(pos->P4().Theta()*TMath::RadToDeg()); + + + } + gBenchmark->Stop("e"); + gBenchmark->Print("e"); + + generator().Summary(); + // internally stored histograms for total ep cross section + TH1D *hWdist = (TH1D*)gDirectory->FindObject("Wdist"); + TH1D *hGenWdist = (TH1D*)gDirectory->FindObject("genWdist"); + + TFile *fout = TFile::Open(Form("out_Dnpee/eD_to_npee_%d.root",(int)ebeamE), "recreate"); + // total ep cross section inputs + if(hWdist)hWdist->Write(); + if(hGenWdist)hGenWdist->Write(); + + // generated event distributions + hW.Write(); + hgE.Write(); + hpP.Write(); + hpTheta.Write(); + hnP.Write(); + hnTheta.Write(); + heP.Write(); + heTheta.Write(); + hpoP.Write(); + hpoTheta.Write(); + fout->Close(); + + + +} diff --git a/jpacPhoto b/jpacPhoto index 8ea3993..756e2ab 160000 --- a/jpacPhoto +++ b/jpacPhoto @@ -1 +1 @@ -Subproject commit 8ea3993d145fd0130cdb95a08a7dc2657f0ea8ac +Subproject commit 756e2ab5f90e590b8130e306898f9c2dee6afb2a