Skip to content

Commit

Permalink
Change max corss section to use only Genetic algoithm, i.e. no grid s…
Browse files Browse the repository at this point in the history 
…earch and switch search variable to costh due to its fixed range
  • Loading branch information
@dglazier
dglazier committed Jun 25, 2021
1 parent 6c3b6e2 commit 19a711e
Show file tree
Hide file tree
Showing 6 changed files with 86 additions and 44 deletions.
3 changes: 1 addition & 2 deletions core/DecayModelQ2W.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -99,8 +99,7 @@ namespace elSpectro{

_gamma = p4beam-p4scat;//can now use getQ2
// std::cout<<"DecayModelQ2W "<<Parent()->Pdg()<<" "<<Parent()->P4().Vect().Unit()<<" "<<_gamma.Vect().Unit()<<" "<<(p4scat + _gstarNuc->P4()).Vect().Unit()<<"initial "<<(p4beam+p4tar).Vect().Unit()<<std::endl;

//calculate photon polarisation
//calculate photon polarisation
auto epsilon = escat::virtualPhotonPolarisation(p4beam,p4tar,p4scat);
auto delta = 2*escat::M2_el()/getQ2()*(1-epsilon);

Expand Down
107 changes: 74 additions & 33 deletions core/DecayModelst.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -72,7 +72,8 @@ namespace elSpectro{
_dsigma=0;
//check above threshold for meson and baryon masses
if( _W < (_meson->P4().M()+_baryon->P4().M()) ) return 0;

//std::cout<<"DecayModelst "<<Parent()->Pdg()<<" "<<_meson->P4().M()<<" "<<_baryon->P4().M()<<std::endl;

//now we can define production/polarisation plane
MomentumVector decayAngles;
kine::electroCMDecay(&Parent()->P4(),_ebeam,_photon,_meson->P4ptr(),&decayAngles);
Expand Down Expand Up @@ -118,9 +119,45 @@ namespace elSpectro{
// _Wmax = ( *(_prodInfo->_target) + *(_prodInfo->_ebeam) ).M();
_Wmax = _prodInfo->_Wmax;

std::cout<<" DecayModelst::FindMaxOfIntensity() Wmin = "<<Wmin<<" Wmax = "<<_Wmax<<" meson mass = "<<M3<<std::endl;
std::cout<<" DecayModelst::FindMaxOfIntensity() Wmin = "<<Wmin<<" Wmax = "<<_Wmax<<" meson mass = "<<M3<<" baryon mass = "<<M4<<" target mass = "<<M2<<std::endl;

auto Fmax = [&M1,&M2,&M3,&M4,&Wmin,this](const double *x)
{
_s = x[0]*x[0];
_W=x[0];
if( _W < Wmin ) return 0.;
if( _W < M3+M4 ) return 0.;
if( _W > _Wmax ) return 0.;

auto currt=kine::tFromcosthW(x[1],_W,M1,M2,M3,M4);
auto myt0=kine::t0(_W,M1,M2,M3,M4);
auto mytmax=kine::tmax(_W,M1,M2,M3,M4);
if(currt>myt0) return 0.;
if(currt<mytmax) return 0.;
if( TMath::IsNaN(x[1]) ) return 0.;

//make smooth function at boundary, but must be lower than max
//this helps with the minimisation
/*Double_t factor=1;
if( x[1] > 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];
if(factor<0) factor = 0;
*/
_t=currt;
auto dt=4* TMath::Sqrt(PgammaCMsq()) * kine::PDK(_W,M3,M4 );

double val = DifferentialXSect()*dt;
if( TMath::IsNaN(val) ) return 0.;
return -(val); //using a minimiser!
};
/* auto Fmax = [&M1,&M2,&M3,&M4,&Wmin,this](const double *x)
{
_s = x[0]*x[0];
_W=x[0];
Expand Down Expand Up @@ -148,11 +185,11 @@ namespace elSpectro{
else _t=x[1];
if(factor<0) factor = 0;
double val = DifferentialXSect()* (myt0-mytmax)*factor;
double val = DifferentialXSect()*(myt0-mytmax)*factor;
if( TMath::IsNaN(val) ) return 0.;
return -(val); //using a minimiser!
};

*/
//First perform grid search for intital values
double Wrange=_Wmax-Wmin;
double tmax=kine::tmax(_Wmax,M1,M2,M3,M4);
Expand All @@ -161,26 +198,27 @@ namespace elSpectro{
double gridW=0;
double gridt=0;
double WtVals[2];
int Npoints=50;
for(int iW=1;iW<Npoints;iW++){
WtVals[0]=Wmin+iW*Wrange/Npoints - Wrange/Npoints/2;//mid point
double tming=kine::t0(WtVals[0],M1,M2,M3,M4);
double tmaxg=kine::tmax(WtVals[0],M1,M2,M3,M4);
double trange= tming-tmaxg;

for(int it=0;it<Npoints;it++){
WtVals[1]=tming-it*trange/Npoints;
auto val = Fmax(WtVals);
// if(it%10==0)std::cout<<WtVals[0]<<" "<<WtVals[1]<<" val "<<TMath::Exp(-val)<<" "<<PhaseSpaceFactor()<<" s "<<_s<<" pgam "<<PgammaCMsq()<<" at Q2 0 = "<<kine::PDK2(_W,0,_target->M())<<std::endl;
if(val<gridMin) {
gridMin=val;
gridW=WtVals[0];
gridt=WtVals[1];
}
}
}

std::cout<<"FindMaxOfProbabilityDistribution grid search max= "<<-gridMin<<" at W = "<<gridW<<" and t = "<<gridt<<" note t0 "<<kine::t0(gridW,M1,M2,M3,M4)<<std::endl;
// int Npoints=500;
// for(int iW=1;iW<Npoints;iW++){
// WtVals[0]=Wmin+iW*Wrange/Npoints - Wrange/Npoints/2;//mid point
// double tming=kine::t0(WtVals[0],M1,M2,M3,M4);
// double tmaxg=kine::tmax(WtVals[0],M1,M2,M3,M4);
// double trange= tming-tmaxg;

// for(int it=0;it<Npoints;it++){
// WtVals[1]=tming-it*trange/Npoints;
// auto val = Fmax(WtVals);
// //if(it%10==0)
// std::cout<<WtVals[0]<<" "<<WtVals[1]<<" val "<<(-val)<<" "<<PhaseSpaceFactor()<<" s "<<_s<<" pgam "<<PgammaCMsq()<<" at Q2 0 = "<<kine::PDK2(_W,0,_target->M())<<std::endl;
// if(val<gridMin) {
// gridMin=val;
// gridW=WtVals[0];
// gridt=WtVals[1];
// }
// }
// }

// std::cout<<"FindMaxOfProbabilityDistribution grid search max= "<<-gridMin<<" at W = "<<gridW<<" and t = "<<gridt<<" note t0 "<<kine::t0(gridW,M1,M2,M3,M4)<<std::endl;
//gridW+=Wrange/100;
//gridt+=tmax/1000;

Expand All @@ -202,11 +240,12 @@ namespace elSpectro{
ROOT::Math::Functor wrapf(Fmax,2);

//variable = W, variable 1 = t
double step[2] = {Wrange/500,0.001};
// double step[2] = {Wrange/500,0.001};
double step[2] = {Wrange/100,2./100};
// starting point

double variable[2] = { gridW,gridt};
//double variable[2] = { 10,-0.1};
// double variable[2] = { gridW,gridt};
double variable[2] = {Wmin+Wrange/2,0};

minimum->SetFunction(wrapf);

Expand All @@ -222,7 +261,7 @@ namespace elSpectro{

auto minVal = minimum->MinValue();
auto minW= xs[0];
auto mint= xs[1];
auto mint= kine::tFromcosthW(xs[1],minW,M1,M2,M3,M4);//xs[1];

std::cout << "Maximum : Probabiltiy Dist at ( W=" << minW << " , t = " << mint << "): "<< -minimum->MinValue() << " note t0 "<<kine::t0(minW,M1,M2,M3,M4)<< std::endl;

Expand All @@ -236,17 +275,18 @@ namespace elSpectro{
const double *xs = minimum->X();

auto minminVal = minimum->MinValue();
minW= xs[0];
mint= xs[1];

std::cout << "Minimum Mass Maximum : Probabiltiy Dist at ( W=" << minW << " , t = " << mint << "): "<< -minimum->MinValue()<< " note t0 "<<kine::t0(minW,M1,M2,M3,M4) << std::endl;

if(minminVal<minVal){
std::cout<<"minmin "<<-minminVal<<" < "<<-minVal<<std::endl;
std::cout << "Minimum Mass Maximum : Probabiltiy Dist at ( W=" << minW << " , t = " << mint << "): "<< -minimum->MinValue()<< " note t0 "<<kine::t0(minW,M1,M2,M3,M4) << std::endl;
std::cout<<"minmin "<<-minminVal<<" < "<<-minVal<<std::endl;
minVal=minminVal;
minW= xs[0];
mint= xs[1];
}
//back to PDg mass if exists
if(_meson->PdgMass()>M3)
M3=_meson->PdgMass();
dynamic_cast<DecayingParticle*>(_meson)->TakePdgMass();

}
Expand All @@ -256,7 +296,8 @@ namespace elSpectro{
std::cout<<"gridMin "<<gridMin<<" "<<minVal<<std::endl;
minVal=gridMin*1.05;
}
return -minVal ;

return -minVal;//*dt ;//and convert to dt
}
/*
void DecayModelst::HistIntegratedXSection(TH1D& hist){
Expand Down
5 changes: 3 additions & 2 deletions core/DecayModelst.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -163,8 +163,9 @@ namespace elSpectro{
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
return _dsigma=PhaseSpaceFactor() * ( MatrixElementsSquared_T() +
(_photonPol->Epsilon()+_photonPol->Delta())*MatrixElementsSquared_L()); //eqn from Seyboth and Wolf
// std::cout<<" DifferentialXSect() "<<PhaseSpaceFactor()<<" "<<" "<<PgammaCMsq()<<std::endl;
return _dsigma=PhaseSpaceFactor() *
( MatrixElementsSquared_T() + (_photonPol->Epsilon()+_photonPol->Delta())*MatrixElementsSquared_L()); //eqn from Seyboth and Wolf
}

SDME* _sdmeMeson={nullptr};
Expand Down
11 changes: 6 additions & 5 deletions core/DecayVectors.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,8 +45,9 @@ namespace elSpectro{


virtual void BoostToParentWithRandPhi(const LorentzVector& parent, LorentzVector& child){
if(parent.P()==0) return; //no boost to be done, or direction

if(parent.P()==0){ return;} //no boost to be done, or direction
//std::cout<<"BoostToParentW "<<child<<" "<<child.M()<<" "<<parent.M()<<std::endl;

//Need axis to rotate around to align
//z-axis with with parent vector
//This is just parent X z-axis
Expand All @@ -58,14 +59,14 @@ namespace elSpectro{
ROOT::Math::AxisAngle rotAroundParent(parent.Vect().Unit(),RandomPhi());
child= rotAroundParent* child;

// auto check2 = LorentzVector(0,0,1,1);
//auto check2 = LorentzVector(0,0,1,1);
//check2=rot1*check2;
// std::cout<<" BoostToParentWithRandPhi "<<check2.Vect().Unit()<<" "<<parent.Vect().Unit()<<" axis "<<axis1<<" rot "<<rot1<<std::endl;
//std::cout<<" BoostToParentWithRandPhi "<<check2.Vect().Unit()<<" "<<parent.Vect().Unit()<<" axis "<<axis1<<" rot "<<rot1<<std::endl;

//boost from parent rest frame
auto boostFromParent=-parent.BoostToCM();
child=boost(child,boostFromParent);//ROOT::Math::VectorUtil::boost;

//std::cout<<"Done BoostToParentW "<<child<<" "<<child.M()<<std::endl;
}

protected:
Expand Down
2 changes: 1 addition & 1 deletion core/JpacModelst.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,7 @@ namespace elSpectro{
double MatrixElementsSquared_T() const override {
_amp->_kinematics->set_mX( GetMeson()->Mass() );
// _amp->_kinematics->set_Q2( get_Q2() );
//std::cout<<"me "<<GetMeson()->Mass()<<" Q2 "<< get_Q2()<<" t "<<get_t()<<" s "<<get_s()<<" W "<<get_W()<<" jpac "<<_amp->_kinematics->Wth()<<std::endl;
//std::cout<<"me "<<GetMeson()->Mass()<<" Q2 "<< get_Q2()<<" t "<<get_t()<<" s "<<get_s()<<" W "<<get_W()<<" jpac "<<_amp->_kinematics->Wth()<<" VAL "<<_amp->probability_distribution(get_s(),get_t())/4<<std::endl;
if(get_W()<_amp->_kinematics->Wth()) return 0;
return _amp->probability_distribution(get_s(),get_t())/4;// Average over initial state helicites;
}
Expand Down
2 changes: 1 addition & 1 deletion core/TwoBodyFlat.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,7 @@ namespace elSpectro{
BoostToParentWithRandPhi(parent,_a);
products[0]->SetP4(_a);
products[1]->SetP4( parent - _a );

// std::cout<<"TwoBody "<<products[0]->P4()<<" "<<products[0]->P4().M()<<" "<<products[1]->P4()<<" "<<products[1]->P4().M()<<std::endl;
return _weight;
}

Expand Down

0 comments on commit 19a711e

Please sign in to comment.