Removed a redundant check in msr2data

src/external/libGapIntegrals/TGapIntegrals.cpp

@@ -74,7 +74,7 @@ TGapDWave::TGapDWave() {
 }
 
 TGapAnSWave::TGapAnSWave() {
-  TAnSWaveGapIntegralDivonne *gapint = new TAnSWaveGapIntegralDivonne();
+  TAnSWaveGapIntegralSuave *gapint = new TAnSWaveGapIntegralSuave();
   fGapIntegral = gapint;
   gapint = 0;
   delete gapint;
@@ -230,16 +230,17 @@ double TGapDWave::operator()(double t, const vector<double> &par) const {
 
     double ds;
     vector<double> intPar; // parameters for the integral, T & Delta(T)
-    intPar.push_back(t);
+    intPar.push_back(2.0*0.08617384436*t); // 2 kB T, kB in meV/K
     intPar.push_back(par[1]*tanh(1.82*pow(1.018*(par[0]/t-1.0),0.51)));
-
+    intPar.push_back(4.0*(t+intPar[1])); // upper limit of energy-integration: cutoff energy
+    intPar.push_back(TMath::PiOver2()); // upper limit of phi-integration
 
 //    double xl[] = {0.0, 0.0}; // lower bound E, phi
-//    double xu[] = {2.0*(t+intPar[1]), 0.5*PI}; // upper bound E, phi
+//    double xu[] = {4.0*(t+intPar[1]), 0.5*PI}; // upper bound E, phi
 
     fGapIntegral->SetParameters(intPar);
 //    ds = 1.0+4.0/PI*fGapIntegral->IntegrateFunc(2, xl, xu);
-    ds = 1.0+4.0/PI*fGapIntegral->IntegrateFunc();
+    ds = 1.0-intPar[2]/intPar[0]*fGapIntegral->IntegrateFunc();
 
     intPar.clear();
 
@@ -300,17 +301,18 @@ double TGapAnSWave::operator()(double t, const vector<double> &par) const {
 
     double ds;
     vector<double> intPar; // parameters for the integral, T & Delta(T)
-    intPar.push_back(t);
+    intPar.push_back(2.0*0.08617384436*t); // 2 kB T, kB in meV/K
     intPar.push_back(par[1]*tanh(1.82*pow(1.018*(par[0]/t-1.0),0.51)));
     intPar.push_back(par[2]);
-
+    intPar.push_back(4.0*(t+(1.0+par[2])*intPar[1])); // upper limit of energy-integration: cutoff energy
+    intPar.push_back(TMath::PiOver2()); // upper limit of phi-integration
 
 //    double xl[] = {0.0, 0.0}; // lower bound E, phi
 //    double xu[] = {4.0*(t+intPar[1]), 0.5*PI}; // upper bound E, phi
 
     fGapIntegral->SetParameters(intPar);
 //    ds = 1.0+4.0/PI*fGapIntegral->IntegrateFunc(2, xl, xu);
-    ds = 1.0+4.0/PI*fGapIntegral->IntegrateFunc();
+    ds = 1.0-intPar[3]/intPar[0]*fGapIntegral->IntegrateFunc();
 
     intPar.clear();
 

src/external/libGapIntegrals/TGapIntegrals.h

@@ -89,7 +89,7 @@ public:
   double operator()(double, const vector<double>&) const;
 
 private:
-  TAnSWaveGapIntegralDivonne *fGapIntegral;
+  TAnSWaveGapIntegralSuave *fGapIntegral;
   mutable vector<double> fTemp;
   mutable vector<double>::const_iterator fTempIter;
   mutable vector<double> fIntegralValues;

src/external/libGapIntegrals/TIntegrator.cpp

@@ -59,13 +59,10 @@ double TDWaveGapIntegralCuhre::IntegrateFunc()
 }
 
 void TDWaveGapIntegralCuhre::Integrand(const int *ndim, const double x[],
-                      const int *ncomp, double f[]) // x = {E, phi}, fPar = {T, Delta(T)}
+                      const int *ncomp, double f[]) // x = {E, phi}, fPar = {twokBT, Delta(T), Ec, phic}
 {
-  double twokt(2.0*0.08617384436*fPar[0]); // kB in meV/K
+  double deltasq(TMath::Power(fPar[1]*TMath::Cos(2.0*x[1]*fPar[3]),2.0));
-  double Ec(4.0*(fPar[0]+fPar[1])); // upper limit of energy-integration: cutoff energy
+  f[0] = 1.0/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[2]*fPar[2]+deltasq)/fPar[0]),2.0);
-  double phic(TMath::PiOver2()); // upper limit of phi-integration
-  double deltasq(TMath::Power(fPar[1]*TMath::Cos(2.0*x[1]*phic),2.0));
-  f[0] = -phic*Ec/(2.0*twokt*TMath::CosH(TMath::Sqrt(x[0]*x[0]*Ec*Ec+deltasq)/twokt)*TMath::CosH(TMath::Sqrt(x[0]*x[0]*Ec*Ec+deltasq)/twokt));
   return;
 }
 
@@ -78,7 +75,7 @@ double TAnSWaveGapIntegralCuhre::IntegrateFunc()
   const double EPSABS (1e-6);
   const unsigned int VERBOSE (0);
   const unsigned int LAST (4);
-  const unsigned int MINEVAL (0);
+  const unsigned int MINEVAL (10000);
   const unsigned int MAXEVAL (1000000);
 
   const unsigned int KEY (13);
@@ -95,13 +92,10 @@ double TAnSWaveGapIntegralCuhre::IntegrateFunc()
 }
 
 void TAnSWaveGapIntegralCuhre::Integrand(const int *ndim, const double x[],
-                      const int *ncomp, double f[]) // x = {E, phi}, fPar = {T, Delta(T),a}
+                      const int *ncomp, double f[]) // x = {E, phi}, fPar = {twokBT, Delta(T),a, Ec, phic}
 {
-  double twokt(2.0*0.08617384436*fPar[0]); // kB in meV/K
+  double deltasq(TMath::Power(fPar[1]*(1.0+fPar[2]*TMath::Cos(4.0*x[1]*fPar[4])),2.0));
-  double Ec(4.0*(fPar[0]+(1.0+fPar[2])*fPar[1])); // upper limit of energy-integration: cutoff energy
+  f[0] = 1.0/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[3]*fPar[3]+deltasq)/fPar[0]),2.0);
-  double phic(TMath::PiOver2()); // upper limit of phi-integration
-  double deltasq(TMath::Power(fPar[1]*(1.0+fPar[2]*TMath::Cos(4.0*x[1]*phic)),2.0));
-  f[0] = -phic*Ec/(2.0*twokt*TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*Ec*Ec+deltasq)/twokt),2.0));
   return;
 }
 
@@ -113,7 +107,7 @@ double TAnSWaveGapIntegralDivonne::IntegrateFunc()
   const double EPSREL (1e-4);
   const double EPSABS (1e-6);
   const unsigned int VERBOSE (0);
-  const unsigned int MINEVAL (0);
+  const unsigned int MINEVAL (10000);
   const unsigned int MAXEVAL (1000000);
   const unsigned int KEY1 (47);
   const unsigned int KEY2 (1);
@@ -139,12 +133,43 @@ double TAnSWaveGapIntegralDivonne::IntegrateFunc()
 }
 
 void TAnSWaveGapIntegralDivonne::Integrand(const int *ndim, const double x[],
-                      const int *ncomp, double f[]) // x = {E, phi}, fPar = {T, Delta(T),a}
+                      const int *ncomp, double f[]) // x = {E, phi}, fPar = {twokBT, Delta(T),a, Ec, phic}
 {
-  double twokt(2.0*0.08617384436*fPar[0]); // kB in meV/K
+  double deltasq(TMath::Power(fPar[1]*(1.0+fPar[2]*TMath::Cos(4.0*x[1]*fPar[4])),2.0));
-  double Ec(4.0*(fPar[0]+(1.0+fPar[2])*fPar[1])); // upper limit of energy-integration: cutoff energy
+  f[0] = 1.0/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[3]*fPar[3]+deltasq)/fPar[0]),2.0);
-  double phic(TMath::PiOver2()); // upper limit of phi-integration
+  return;
-  double deltasq(TMath::Power(fPar[1]*(1.0+fPar[2]*TMath::Cos(4.0*x[1]*phic)),2.0));
+}
-  f[0] = -phic*Ec/(2.0*twokt*TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*Ec*Ec+deltasq)/twokt),2.0));
+
+std::vector<double> TAnSWaveGapIntegralSuave::fPar;
+
+double TAnSWaveGapIntegralSuave::IntegrateFunc()
+{
+  const unsigned int NCOMP(1);
+  const double EPSREL (1e-4);
+  const double EPSABS (1e-6);
+  const unsigned int VERBOSE (0);
+  const unsigned int LAST (4);
+  const unsigned int MINEVAL (10000);
+  const unsigned int MAXEVAL (1000000);
+
+  const unsigned int NNEW (1000);
+  const double FLATNESS (25.);
+
+  int nregions, neval, fail;
+  double integral[NCOMP], error[NCOMP], prob[NCOMP];
+
+  Suave(fNDim, NCOMP, Integrand,
+    EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
+    NNEW, FLATNESS,
+    &nregions, &neval, &fail, integral, error, prob);
+
+  return integral[0];
+}
+
+void TAnSWaveGapIntegralSuave::Integrand(const int *ndim, const double x[],
+                      const int *ncomp, double f[]) // x = {E, phi}, fPar = {twokBT, Delta(T),a, Ec, phic}
+{
+  double deltasq(TMath::Power(fPar[1]*(1.0+fPar[2]*TMath::Cos(4.0*x[1]*fPar[4])),2.0));
+  f[0] = 1.0/TMath::Power(TMath::CosH(TMath::Sqrt(x[0]*x[0]*fPar[3]*fPar[3]+deltasq)/fPar[0]),2.0);
   return;
 }

src/external/libGapIntegrals/TIntegrator.h

@@ -170,6 +170,20 @@ class TAnSWaveGapIntegralDivonne {
 
 };
 
+class TAnSWaveGapIntegralSuave {
+  public:
+    TAnSWaveGapIntegralSuave() : fNDim(2) {}
+    ~TAnSWaveGapIntegralSuave() { fPar.clear(); }
+    void SetParameters(const std::vector<double> &par) { fPar=par; }
+    static void Integrand(const int*, const double[], const int*, double[]);
+    double IntegrateFunc();
+
+  protected:
+    static vector<double> fPar;
+    unsigned int fNDim;
+
+};
+
 // To be integrated: x*y dx dy
 
 class T2DTest : public TMCIntegrator {

src/msr2data.cpp

@@ -178,8 +178,6 @@ int msr2data_doFitting(vector<string> &arg, bool &chainfit)
       temp = -1;
       chainfit = false;
       iter = arg.erase(iter);
-      if (iter == arg.end())
-        break;
     }
     else if (!iter->substr(0,4).compare("fit-")) {
       if (temp) {
@@ -194,8 +192,6 @@ int msr2data_doFitting(vector<string> &arg, bool &chainfit)
       iss.str(s);
       iss >> temp;
       iter = arg.erase(iter);
-      if (iter == arg.end())
-        break;
     } else {
       iter++;
     }