Fix normalization.

event.cc

@@ -38,9 +38,9 @@ event::decayAmplitude(int reflectivity, int l, int m) const
     spherical = ROOT::Math::sph_legendre(l, m, this->theta);
 
   // This absorbs the factor 2 from e^i \pm e^-i = 2 {i sin, cos}
-  double factor = .5; //sqrt(nrdevents);
-  if (m != 0)
-    factor *= sqrt(2.);    
+  double factor = sqrt(2); //.5; //sqrt(nrdevents);
+  if (m == 0)
+    factor = 1;
 
   double result;
   if (reflectivity == +1)

likelihood.cc

@@ -124,7 +124,7 @@ likelihood::MCweight(int reflectivity, const wave& w1, const wave& w2) const
     = flatMC ? flatMCevents : binnedMCevents[currentBin];
 
   size_t nEv = pMCevents.size();
-#pragma omp parallel for reduction (+:sum)
+  //#pragma omp parallel for reduction (+:sum)
   for (size_t i = 0; i < nEv; i++)
     {
       // Forming this sum as REAL sum and with no conjugate, because
@@ -137,14 +137,14 @@ likelihood::MCweight(int reflectivity, const wave& w1, const wave& w2) const
     }
 
   /*
-    cout << "calculated MCweight " << sum / pMCevents.size()
+  cout << "calculated MCweight " << 16*atan(1) * sum / nEv
     << " from " << pMCevents.size() << " MC events for "
     << "(l1,m1,l2,m2) = "
     << "(" << w1.l << "," << w1.m << "," << w2.l << "," << w2.m << ")"
     << endl;
   */
 
-  return weights[id] = sum / pMCevents.size();
+  return weights[id] = 4*M_PI * sum / nEv;
 }
 
 #if 0

likelihood.h

@@ -79,6 +79,9 @@ class likelihood : public ROOT::Minuit2::FCNBase {
   size_t
   eventsInBin() const { return binnedRDevents[currentBin].size(); }
 
+  size_t
+  MCeventsInBin() const { return (size_t)(binnedMCevents[currentBin].size()/* / binnedEtaAcc[currentBin]*/); }
+
   void
   clearWeights() { weights.clear(); }
 

myFit.cc

@@ -123,11 +123,13 @@ class combinedLikelihood : public ROOT::Minuit2::FCNBase {
 	// RD likelihood calculation as follows.
 	result += (myLs[i]->eventsInBin()*log(BR) + myLs[i]->calc_rd_likelihood(x)
 		   - BR*myLs[i]->calc_mc_likelihood(x));
+	//cout << "nevents = " << myLs[i]->eventsInBin() << " likelihood = " << myLs[i]->calc_rd_likelihood(x) << " - " << myLs[i]->calc_mc_likelihood(x) << endl;
       }
     return -result;
   }
 
   size_t getNChannels() const { return myLs.size(); }
+  size_t MCeventsInBin() const { size_t sum = 0; for (size_t i = 0; i < myLs.size(); i++) sum += myLs[i]->MCeventsInBin(); return sum; }
 };
 
 
@@ -432,12 +434,14 @@ myFit()
   TTree *outTree = new TTree("tFitResults", "fit results tree");
   double *values = new double[lastIdx];
   TMatrixDSym *covMat = new TMatrixDSym(lastIdx);
+  UInt_t NMCevents;
   char branchDesc[99];
   snprintf(branchDesc, 99, "values[%zd]/D", lastIdx);
   outTree->Branch("massLow", &massLow, "massLow/D");
   outTree->Branch("massHigh", &massHigh, "massHigh/D");
   outTree->Branch("values", values, branchDesc);
   outTree->Branch("covMat", "TMatrixDSym", &covMat);
+  outTree->Branch("NMCevents", &NMCevents, "NMCevents/i");
   outTree->GetUserInfo()->Add(new fitInfo(modelName, ws, startingValues, nBins, threshold, binWidth, lastIdx));
 
   /*  TTree* predictTree = new TTree("predictTree","Weighted MC");
@@ -734,7 +738,10 @@ myFit()
 	  MCevents.reserve(NFLATMCEVENTS);
 	  for (int iMC = 0; iMC < NFLATMCEVENTS; iMC++)
 	    {
-	      event e(acos(gRandom->Uniform(-1,1)), gRandom->Uniform(-M_PI, M_PI));
+	      double costh = gRandom->Uniform(-1,1);
+	      double phi = gRandom->Uniform(-M_PI, M_PI);
+	      //cout << costh << " " << acos(costh) << " " << phi << endl;
+	      event e(acos(costh), phi);
 	      MCevents.push_back(e);
 	    }
 	}
@@ -906,6 +913,7 @@ myFit()
 		    }
 	      }
 
+	  NMCevents = myL.MCeventsInBin();
 	  outTree->Fill();
 
 	  for (int j = 0; j < minuit->GetNumberTotalParameters(); j++)

predict.cc

@@ -8,32 +8,12 @@ using namespace std;
 
 #include "TFile.h"
 #include "TTree.h"
-#include "TRandom1.h"
+#include "TLorentzVector.h"
 
 #include "gHist.h"
 #include "fitInfo.h"
 
 
-double
-decayAmplitude(int reflectivity, int l, int m, double costh, double phi)
-{
-  double spherical;
-  spherical = ROOT::Math::sph_legendre(l, m, acos(costh));
-  double factor = .5; //sqrt(nrdevents);
-  if (m != 0)
-    factor *= sqrt(2.);
-
-  double result;
-  if (reflectivity == +1)
-    result = factor*spherical*sin(m*phi);
-  else
-    result = factor*spherical*cos(m*phi);
-
-  //lookupTable[id] = result;
-  return result; // * this->tPrime*exp(-7.*this->tPrime);
-}
-
-
 int
 main()
 {
@@ -71,20 +51,25 @@ main()
 
   double *values = new double[info->getNvars()];
   double massLow, massHigh;
+  UInt_t NMCevents;
   treeResult->SetBranchAddress("massLow", &massLow);
   treeResult->SetBranchAddress("massHigh", &massHigh);
   treeResult->SetBranchAddress("values", values);
+  treeResult->SetBranchAddress("NMCevents", &NMCevents);
+
 
   vector<double> lowerBounds;
   vector<double> upperBounds;
   vector<vector<double> > allValues;
+  vector<UInt_t> NMCperBin;
   for (Long_t i = 0; i < treeResult->GetEntries(); i++)
     {
       // This assumes the bins are sorted.
       treeResult->GetEntry(i);
       lowerBounds.push_back(massLow);
       upperBounds.push_back(massHigh);
       allValues.push_back(vector<double>(values,values+info->getNvars()));
+      NMCperBin.push_back(NMCevents);
     }
 
   const char *fnMC = "/data/zup/diefenbach/trees/EtaPr3Pi.root";
@@ -139,17 +124,30 @@ main()
   float pg2Y; t->SetBranchAddress("pg2Y", &pg2Y);
   float pg2Z; t->SetBranchAddress("pg2Z", &pg2Z);
   float pg2E; t->SetBranchAddress("pg2E", &pg2E);
+  float pPim3X; t->SetBranchAddress("pPim3X", &pPim3X);
+  float pPim3Y; t->SetBranchAddress("pPim3Y", &pPim3Y);
+  float pPim3Z; t->SetBranchAddress("pPim3Z", &pPim3Z);
+  float pPim3E; t->SetBranchAddress("pPim3E", &pPim3E);
+  float pPip3X; t->SetBranchAddress("pPip3X", &pPip3X);
+  float pPip3Y; t->SetBranchAddress("pPip3Y", &pPip3Y);
+  float pPip3Z; t->SetBranchAddress("pPip3Z", &pPip3Z);
+  float pPip3E; t->SetBranchAddress("pPip3E", &pPip3E);
+  float pEta3X; t->SetBranchAddress("pEta3X", &pEta3X);
+  float pEta3Y; t->SetBranchAddress("pEta3Y", &pEta3Y);
+  float pEta3Z; t->SetBranchAddress("pEta3Z", &pEta3Z);
+  float pEta3E; t->SetBranchAddress("pEta3E", &pEta3E);
   float pPimX; t->SetBranchAddress("pPimX", &pPimX);
   float pPimY; t->SetBranchAddress("pPimY", &pPimY);
   float pPimZ; t->SetBranchAddress("pPimZ", &pPimZ);
   float pPimE; t->SetBranchAddress("pPimE", &pPimE);
 
   TFile *fOutput = TFile::Open("predict.root", "RECREATE");
 
-  for (Long_t i = 0; i < t->GetEntries(); i++)
+  Long_t nTotal = t->GetEntries();
+  for (Long_t i = 0; i < nTotal; i++)
     {
-      if ((i+1)*10LL/t->GetEntries() > i*10/t->GetEntries())
-	cout << (i+1)*100/t->GetEntries() << "% " << flush;
+      if ((i+1)*10LL/nTotal > i*10/nTotal)
+	cout << (i+1)*100/nTotal << "% " << flush;
 
       t->GetEntry(i);
       if (!acc)
@@ -181,24 +179,37 @@ main()
 	  continue;
 	}
 
-      double weight = info->getWaveSet().getEventWeight(allValues[idx], e);
+      TLorentzVector lvEta3(pEta3X, pEta3Y, pEta3Z, pEta3E);
+      TLorentzVector lvPim3(pPim3X, pPim3Y, pPim3Z, pPim3E);
+      TLorentzVector lvPip3(pPip3X, pPip3Y, pPip3Z, pPip3E);
+
+      double weight = info->getWaveSet()[0].getEventWeight(allValues[idx], e) / NMCperBin[idx];
       gHist.Fill("hMvsCosth", "m vs cos th",
 		 info->getNbins(), info->getLower(), info->getUpper(), 100, -1, 1,
 		 mX, costh, weight);
-      gHist.Fill("hpvz", "pvz", 500, -100, 0,
-		 pvz, weight);
-      gHist.Fill("htPr", "t'", 1000, 0, 1,
-		 tPr, weight);
-      gHist.Fill("hXYECAL", "XY ECAL", 200, -.04, .04, 200, -0.04, 0.04,
-		 pg1X / pg1Z, pg1Y / pg1Z, weight);
-      gHist.Fill("hXYECAL", "XY ECAL", 200, -.04, .04, 200, -0.04, 0.04,
-		 pg2X / pg2Z, pg2Y / pg2Z, weight);
-      gHist.Fill("hEvsAngle", "E vs angle", 200, 0, 200, 100, 0, 0.06,
-		 pg1E, hypot(pg1X, pg1Y) / pg1Z, weight);
-      gHist.Fill("hEvsAngle", "E vs angle", 200, 0, 200, 100, 0, 0.06,
-		 pg2E, hypot(pg2X, pg2Y) / pg2Z, weight);
-      if (mX < 2)
+      //if (1.9 < mX && mX < 2.1)
 	{
+	  gHist.Fill("hpvz", "pvz", 500, -100, 0,
+		     pvz, weight);
+	  gHist.Fill("htPr", "t'", 1000, 0, 1,
+		     tPr, weight);
+	  gHist.Fill("hXYECAL", "XY ECAL", 200, -.04, .04, 200, -0.04, 0.04,
+		     pg1X / pg1Z, pg1Y / pg1Z, weight);
+	  gHist.Fill("hXYECAL", "XY ECAL", 200, -.04, .04, 200, -0.04, 0.04,
+		     pg2X / pg2Z, pg2Y / pg2Z, weight);
+	  gHist.Fill("hEvsAngle", "E vs angle", 200, 0, 200, 100, 0, 0.06,
+		     pg1E, hypot(pg1X, pg1Y) / pg1Z, weight);
+	  gHist.Fill("hEvsAngle", "E vs angle", 200, 0, 200, 100, 0, 0.06,
+		     pg2E, hypot(pg2X, pg2Y) / pg2Z, weight);
+
+	  gHist.Fill("hDalitz", "Dalitz Plot;pi+eta;pi-eta", 50, 0.45, 0.7, 50, 0.45, 0.7,
+		     (lvEta3 + lvPip3).M2(), (lvEta3 + lvPim3).M2());
+	  gHist.Fill("hDalitzFlipped", "Dalitz Plot;pi+eta;pi-pi+", 50, 0.45, 0.75, 50, 0.06, 0.18,
+		     (lvEta3 + lvPip3).M2(), (lvPip3 + lvPim3).M2());
+
+	  //}
+	  //      if (mX < 2)
+	  //{
 	  gHist.Fill("hpPimE", "energy of non-eta' #pi-",
 		     500, 0, 200,
 		     pPimE, weight);