From 8a7712905c9047fc0b0444ffee77c614d2b4866b Mon Sep 17 00:00:00 2001
From: dglazier <Derek.Glazier@glasgow.ac.uk>
Date: Fri, 11 Jun 2021 15:00:19 +0100
Subject: [PATCH] tidy up

---
 core/DecayModelQ2W.cpp       | 19 +-------
 core/DecayModelst.cpp        | 18 --------
 core/DistVirtPhotFlux_xy.cpp | 89 +-----------------------------------
 core/DistVirtPhotFlux_xy.h   |  3 --
 4 files changed, 3 insertions(+), 126 deletions(-)

diff --git a/core/DecayModelQ2W.cpp b/core/DecayModelQ2W.cpp
index 929a568..5761132 100644
--- a/core/DecayModelQ2W.cpp
+++ b/core/DecayModelQ2W.cpp
@@ -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},
@@ -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
@@ -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{
diff --git a/core/DecayModelst.cpp b/core/DecayModelst.cpp
index 69d7986..f35f0aa 100644
--- a/core/DecayModelst.cpp
+++ b/core/DecayModelst.cpp
@@ -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
@@ -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!
       };
     
diff --git a/core/DistVirtPhotFlux_xy.cpp b/core/DistVirtPhotFlux_xy.cpp
index 08d36a1..6b02eaa 100644
--- a/core/DistVirtPhotFlux_xy.cpp
+++ b/core/DistVirtPhotFlux_xy.cpp
@@ -16,7 +16,7 @@ namespace elSpectro{
     _ebeam(eb),
     _mTar(mion)
    {
-     // SetWThreshold(Wmin);
+
    }
   
   void DistVirtPhotFlux_xy::SetWThreshold(double Wmin){
@@ -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;
@@ -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();
     
@@ -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
@@ -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
diff --git a/core/DistVirtPhotFlux_xy.h b/core/DistVirtPhotFlux_xy.h
index 04ac2e4..9566b40 100644
--- a/core/DistVirtPhotFlux_xy.h
+++ b/core/DistVirtPhotFlux_xy.h
@@ -19,7 +19,6 @@ 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 {
@@ -27,7 +26,6 @@ namespace elSpectro{
     }
     
     dist_pair SamplePair()   noexcept final {
-      //    _forIntegral==false ? FindWithAcceptReject() : FindFlat();
       FindWithAcceptReject();
       return _xy;
     }
@@ -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;}