Skip to content

Commit

Permalink
tidy up
Browse files Browse the repository at this point in the history
  • Loading branch information
@dglazier
dglazier committed Jun 11, 2021
1 parent b90202b commit 8a77129
Show file tree
Hide file tree
Showing 4 changed files with 3 additions and 126 deletions.
19 changes: 1 addition & 18 deletions core/DecayModelQ2W.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,15 +18,6 @@ namespace elSpectro{
}
///////////////////////////////////////////////////////
///constructor includes subseqent decay of Ngamma* system
/* DecayModelQ2W::DecayModelQ2W( double thresh, DecayModel* gNmodel) :
_threshold{thresh},
DecayModel{{ new DecayingParticle{-2211,primary} },{11}}
{
_name={"DecayModelQ2W_with_primary_decay"};
Init();
}*/
///////////////////////////////////////////////////////
///constructor includes subseqent decay of Ngamma* system
DecayModelQ2W::DecayModelQ2W( double thresh,
DecayModel* gNmodel,DecayVectors* gNdecayer) :
_threshold{thresh},
Expand Down Expand Up @@ -115,7 +106,6 @@ 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
Expand Down Expand Up @@ -193,14 +183,7 @@ namespace elSpectro{
std::cout<<"DecayModelQ2W::FindExcitationSpectra() result "<<hist.GetMaximum()<<" "<<hist.GetBinCenter(hist.GetMaximumBin())<<" "<<hist.GetNbinsX()<<std::endl;
hist.SetName("Wdist");

/*
TFile ff("debug.root","recreate");
hist.Write();
histpeak.Write();
histlow.Write();
ff.Close();
*/


_Wrealphoto_Dist.reset( new DistTH1(hist) );
}
else{
Expand Down
18 changes: 0 additions & 18 deletions core/DecayModelst.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -125,25 +125,14 @@ namespace elSpectro{
{
_s = x[0]*x[0];
_W=x[0];
//if( _W < Wmin ) return 1E2;
//if( _W < M3+M4 ) return 1E2;
if( _W < Wmin ) return 0.;
if( _W < M3+M4 ) return 0.;
// 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)<TMath::Abs(x[1]-myt0) )
curr_t = 0;
else x[1] = mytmax- myt0;
}*/

// if( TMath::IsNaN(x[1]) ) return 1E2;
if( TMath::IsNaN(x[1]) ) return 0.;

//make smooth function at boundary, but must be lower than max
Expand All @@ -152,23 +141,16 @@ namespace elSpectro{
if( x[1] > myt0){
_t=myt0;
factor*=(1-1*(x[1]-myt0)); //lower result
//return 0.;
// std::cout<<"too low "<<x[1]<<" "<<_t<<" "<<myt0<<std::endl;
}
else if( x[1]< mytmax ){
_t=mytmax;
factor*=(1-1*(mytmax-x[1])); //lower result
// std::cout<<"too high "<<x[1]<<" "<<_t<<" "<<mytmax<<std::endl;
//return 0.;
}
else _t=x[1];
if(factor<0) factor = 0;

double val = DifferentialXSect()* (myt0-mytmax)*factor;
//std::cout<<factor<<" "<<val<<" "<<(myt0-mytmax)<<"t "<<_t<<" "<<myt0<<" "<<mytmax<<" W "<<_W<<" x1 "<<x[1]<<" "<<DifferentialXSect()<<std::endl;
//if(val==0)return 1E2;
if( TMath::IsNaN(val) ) return 0.;
//if( TMath::IsNaN(val) ) return 1E2;
return -(val); //using a minimiser!
};

Expand Down
89 changes: 2 additions & 87 deletions core/DistVirtPhotFlux_xy.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ namespace elSpectro{
_ebeam(eb),
_mTar(mion)
{
// SetWThreshold(Wmin);

}

void DistVirtPhotFlux_xy::SetWThreshold(double Wmin){
Expand Down Expand Up @@ -94,11 +94,6 @@ namespace elSpectro{

std::cout<<"DistVirtPhotFlux_xy max "<<_max_val<<" within y limits "<<TMath::Exp(_lnymin)<<" "<<TMath::Exp(_lnymax)<<" minimum possible x "<<_maxPossiblexRange<<std::endl;
std::cout<<" Other limits : "<<std::endl;

// if(_maxPossiblexRange==0)
// _maxPossiblexRange=200;
//else
// _maxPossiblexRange=-TMath::Log(_maxPossiblexRange);

if(_requestQ2min!=0)std::cout<<"\tQ2min "<<_requestQ2min<<std::endl;
if(_requestQ2max!=0)std::cout<<"\tQ2max "<<_requestQ2max<<std::endl;
Expand All @@ -117,8 +112,6 @@ namespace elSpectro{
//Finally, Integrate over photon flux
auto xvar = RooRealVar("x","x",-(_lnxmax-_lnxmin)/2,_lnxmin,_lnxmax,"");
auto yvar = RooRealVar("y","y",-(_lnymax-_lnymin)/2,_lnymin,_lnymax,"");
// auto xvar = RooRealVar("x","x",TMath::Exp((_lnxmax-_lnxmin)/2),TMath::Exp(_lnxmin),TMath::Exp(_lnxmax),"");
//auto yvar = RooRealVar("y","y",TMath::Exp((_lnymax-_lnymin)/2),TMath::Exp(_lnymin),TMath::Exp(_lnymax),"");
xvar.Print();
yvar.Print();

Expand Down Expand Up @@ -179,7 +172,7 @@ namespace elSpectro{
(_val=escat::flux_dlnxdlny(_ebeam,lnx,lny)) ) {
if(_val>_max_val){
_max_val=_val;
std::cout<<" MAX REACHED "<<_val<<" "<<_max_val<<std::endl;
std::cout<<"DistVirtPhotFlux_xy::FindWithAcceptReject() MAX REACHED "<<_val<<" "<<_max_val<<std::endl;
exit(0);
}
getRandomXY(); //sample another pair
Expand All @@ -198,83 +191,5 @@ namespace elSpectro{

}

/*
void DistVirtPhotFlux_xy::FindFlat(){
_val=0;
_xy=std::make_pair(-1,-1); //unphysical should never be used
double y = gRandom->Uniform( TMath::Exp(_lnymin),TMath::Exp(_lnymax) );
//double y = gRandom->Uniform( 1E-8,1 );
//if(y<0.1) return;
//calculate the fraction of x-space available
//now calculate x limits
//y = r/2ME and x = Q2/r
double avail_xmax = XMax(y);
double avail_xmin = XMin(y);
//if(avail_xmin>=1){ return; }
//_maxPossiblexRange=0;
double x=-1;
if(_maxPossiblexRange)
//for efficiency we need not sample below lowest possible x value
x = gRandom->Uniform(_maxPossiblexRange,1);
else
x = gRandom->Uniform(0,1);
//check if we are within allowed x-range
if( x<avail_xmax && x>avail_xmin )
_val=escat::flux_dxdy(_ebeam,x,y);
//return x and y values
_xy=std::make_pair(x,y);
return;
}
*/
/*
void DistVirtPhotFlux_xy::FindFlat(){
TF1 sampleX("sampleX","TMath::Exp(-x[0])",0,1);
sampleX.SetRange(0,1);
sampleX.SetNpx(1E4);
TF1 sampleY("sampleY","TMath::Exp(-x[0])",0,1);
sampleY.SetRange(TMath::Exp(_lnymin),TMath::Exp(_lnymax));
sampleY.SetNpx(1E4);
_val=0;
_xy=std::make_pair(-1,-1); //unphysical should never be used
double y = sampleY.GetRandom();
double sampleYVal=sampleY.Eval(y)/sampleYIntegral;
//calculate the fraction of x-space available
//now calculate x limits
//y = r/2ME and x = Q2/r
double avail_xmax = XMax(y);
double avail_xmin = XMin(y);
if(avail_xmin>=1){ return; }
double x=-1;
while(x<_maxPossiblexRange)
x = sampleX.GetRandom();
double sampleXVal=sampleX.Eval(x)/sampleXIntegral; //value of PDF
//check if we are within allowed x-range
if( x<avail_xmax && x>avail_xmin )
_val=escat::flux_dxdy(_ebeam,x,y)/sampleYVal/sampleXVal;
//return x and y values
_xy=std::make_pair(x,y);
return;
}

*/
}//namespace
3 changes: 0 additions & 3 deletions core/DistVirtPhotFlux_xy.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,15 +19,13 @@ namespace elSpectro{

public :

//DistVirtPhotFlux_xy(double ebeam,float xmin,float xmax,float ymin,float ymax);
DistVirtPhotFlux_xy(double eb, double mion, double Wmin);

double SampleSingle() noexcept final {
return 0;
}

dist_pair SamplePair() noexcept final {
// _forIntegral==false ? FindWithAcceptReject() : FindFlat();
FindWithAcceptReject();
return _xy;
}
Expand All @@ -48,7 +46,6 @@ namespace elSpectro{
double GetWMin() const noexcept {return TMath::Sqrt(_Wthresh2);}

void FindWithAcceptReject();
// void FindFlat();

void SetElecE(double ee){_ebeam=ee;}
void SetM(double m){_mTar=m;}
Expand Down

0 comments on commit 8a77129

Please sign in to comment.