add SuppressPhaseSpace option to speed up tests

dglazier · Jun 29, 2021 · 1195a5e · 1195a5e
1 parent 19a711e
commit 1195a5e
Show file tree

Hide file tree

Showing 3 changed files with 7 additions and 74 deletions.
diff --git a/core/DecayModelst.cpp b/core/DecayModelst.cpp
@@ -136,60 +136,14 @@ namespace elSpectro{
 	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];
-	if( _W < Wmin ) return 0.;
-	if( _W < 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(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;
-
-	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);
@@ -198,30 +152,7 @@ namespace elSpectro{
       double gridW=0;
       double gridt=0;
       double WtVals[2];
-      // 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;
-
+
       ROOT::Math::Minimizer* minimum =
 	ROOT::Math::Factory::CreateMinimizer("Genetic", "");
       //	ROOT::Math::Factory::CreateMinimizer("Minuit2", "Combined");
@@ -297,7 +228,7 @@ namespace elSpectro{
 	minVal=gridMin*1.05;
       }
 
-      return -minVal;//*dt ;//and convert to dt
+      return -minVal;
   }
   /*
   void DecayModelst::HistIntegratedXSection(TH1D& hist){

diff --git a/core/Manager.h b/core/Manager.h
@@ -83,6 +83,7 @@ namespace elSpectro{
 
 
      void SetModelForMassPhaseSpace(DecayModel* amodel){_massPhaseSpace.SetModel(amodel);}
+    void SuppressPhaseSpace(double val){_massPhaseSpace.SuppressPhaseSpace(val);}
      void  FindMassPhaseSpace(double parentM,const  DecayModel* amodel) {
        _massPhaseSpace.Find(parentM,amodel);
      }

diff --git a/core/MassPhaseSpace.h b/core/MassPhaseSpace.h
@@ -23,7 +23,7 @@ namespace elSpectro{
     void Print(){
       std::cout<<"MassPhaseSpace number calcs= "<<_weightCalcN<<" number of successes = "<<_successN<<" ratio  ="<< double(_successN)/_weightCalcN <<std::endl;
     }
-
+    void SuppressPhaseSpace(double val){_suppressPhaseSpace=val;}
   private:
 
     friend Manager; //only Manager can construct and use a MassPhaseSpace
@@ -51,7 +51,7 @@ namespace elSpectro{
       //as W dependence accounted for elsewere
       //PhaseSpaceWeight will try alternative masses
       double wee=0;
-      while( (wee=PhaseSpaceWeight(parentM)) < gRandom->Uniform()*max )
+      while( (wee=PhaseSpaceWeight(parentM)) < gRandom->Uniform()*max*_suppressPhaseSpace )
 	{
 	  // if(wee>max){
 	  // _sampledMax=wee;
@@ -95,6 +95,7 @@ namespace elSpectro{
     double _sampledMax = 0;
     mutable long _weightCalcN=0;
     mutable long _successN=0;
+    double _suppressPhaseSpace=1;
 
     ClassDef(elSpectro::MassPhaseSpace,1); //class MassPhaseSpace
   };