diff --git a/event.cc b/event.cc index 804d559..3c702ae 100644 --- a/event.cc +++ b/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) diff --git a/likelihood.cc b/likelihood.cc index 4d66e31..4f7d145 100644 --- a/likelihood.cc +++ b/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 diff --git a/likelihood.h b/likelihood.h index 7420b36..f54d321 100644 --- a/likelihood.h +++ b/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(); } diff --git a/myFit.cc b/myFit.cc index ecce00c..1f27b01 100644 --- a/myFit.cc +++ b/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++) diff --git a/predict.cc b/predict.cc index 2e0f593..8d84d73 100644 --- a/predict.cc +++ b/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,13 +51,17 @@ 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 lowerBounds; vector upperBounds; vector > allValues; + vector NMCperBin; for (Long_t i = 0; i < treeResult->GetEntries(); i++) { // This assumes the bins are sorted. @@ -85,6 +69,7 @@ main() lowerBounds.push_back(massLow); upperBounds.push_back(massHigh); allValues.push_back(vector(values,values+info->getNvars())); + NMCperBin.push_back(NMCevents); } const char *fnMC = "/data/zup/diefenbach/trees/EtaPr3Pi.root"; @@ -139,6 +124,18 @@ 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); @@ -146,10 +143,11 @@ main() 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);