inlcuded pipi fit model

fitModel.cc

@@ -193,7 +193,7 @@ fitModelKKbarSD::evaluateAt(double mass_, const vector<double>& x)
   double m2 = mK;
   double phaseSpace = breakupMomentum(mass*mass, m1, m2);
 
-  double mred = mass-2*mK;
+  //double mred = mass-2*mK;
 
   //par = &x[13];
 
@@ -256,6 +256,106 @@ fitModelKKbarSD::valueForWave(const char* name) const
   exit(1);
 }
 
+class fitModelpipiSD : public fitModel {
+public:
+  fitModelpipiSD() : mass(0), Swave(0), SwaveBG(0), Dwave(0), DwaveBG(0) {};
+  void evaluateAt(double mass_, const std::vector<double>& x);
+  std::complex<double> valueForWave(const char* name) const;
+private:
+  double mass;
+  std::complex<double> Swave;
+  std::complex<double> SwaveBG;
+  std::complex<double> Dwave;
+  std::complex<double> DwaveBG;
+  std::complex<double> phaseS;
+
+  std::complex<double> f2;
+
+  std::complex<double> f0980;
+
+  std::complex<double> f01370;
+  std::complex<double> f01500;
+  std::complex<double> f01710;
+};
+
+void
+fitModelpipiSD::evaluateAt(double mass_, const vector<double>& x)
+{
+  mass = mass_;
+
+  double m1 = mPi;
+  double m2 = mPi;
+  double phaseSpace = breakupMomentum(mass*mass, m1, m2);
+
+  //double mred = mass-2*mPi;
+
+  //par = &x[13];
+
+  //cout << x[0] << endl;
+
+  DwaveBG = complex<double>(x[4],x[5]) * pow(phaseSpace, 2.5)/mass * exp(-1*x[6]*phaseSpace - x[7]*phaseSpace*phaseSpace);
+
+  f2 = complex<double>(x[0],x[1]) * BW(mass*mass, m1, m2, 2, x[2], x[3]);
+
+  Dwave = f2 + DwaveBG;
+
+  SwaveBG = x[8] * pow(phaseSpace, 0.5)/mass * exp(-1*x[9]*phaseSpace - x[10]*phaseSpace*phaseSpace);
+
+  //SwaveBG = x[8]*BW(mass*mass, m1, m2, 0, x[9], x[10]);
+
+  f0980 = complex<double>(x[11], x[12])*BW(mass*mass, m1, m2, 0, x[13], x[14]);
+  //f01370 = complex<double>(x[15], x[16])*BW(mass*mass, m1, m2, 0, x[17], x[18]);
+
+  //f01500 = complex<double>(x[19], x[20])*BW(mass*mass, m1, m2, 0, x[21], x[22]);
+  //f01710 = complex<double>(x[23], x[24])*BW(mass*mass, m1, m2, 0, x[25], x[26]);
+
+  Swave = SwaveBG + f0980;// + f01370 + f01500 + f01710;
+  //+ SwaveBG;
+
+  phaseS = 1;
+  if (abs(Swave) != 0)
+    phaseS = Swave / abs(Swave);
+  Swave /= phaseS;
+
+  Dwave /= phaseS;
+  //  f01500 /= phaseS;
+  //  f01370 /= phaseS;
+  //  f01710 /= phaseS;
+
+  //SwaveBG = complex<double>(x[21],x[22])*(x[23] + mass*x[24] + mass*mass*x[25]);
+}
+
+
+complex<double>
+fitModelpipiSD::valueForWave(const char* name) const
+{
+  if (!strcmp(name, "S0"))
+    return Swave;
+  else if (!strcmp(name, "S0BG"))
+    return SwaveBG;
+  else if (!strcmp(name, "D0"))
+    return Dwave;
+  else if (!strcmp(name, "D0BG"))
+    return DwaveBG;
+  else if (!strcmp(name, "f2"))
+    return f2;
+  else if (!strcmp(name, "f0980"))
+    return f0980;
+  else if (!strcmp(name, "f01370"))
+    return f01370;
+  else if (!strcmp(name, "f01500"))
+    return f01500;
+  else if (!strcmp(name, "f01710"))
+    return f01710;
+  /*else if (!strcmp(name, "f2prime"))
+    return f2prime;
+    else if (!strcmp(name, "f01370"))
+    return f01370;*/
+
+  cerr << "unknown wave '" << name << "' requested" << endl;
+  exit(1);
+}
+
 fitModel*
 fitModel::getFitModelForName(const string& name)
 {
@@ -265,6 +365,11 @@ fitModel::getFitModelForName(const string& name)
   if (!strcmp(name.c_str(), "etaPpi"))
     return new fitModelEtaPpi();
 
+  if (!strcmp(name.c_str(), ""))
+    return new fitModelpipiSD();
+  if (!strcmp(name.c_str(), "pipiSD"))
+    return new fitModelpipiSD();
+
   if (!strcmp(name.c_str(), "KKbarSD"))
     return new fitModelKKbarSD();