diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1DQ2_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1DQ2_nu.cxx index 2764225..ee52182 100644 --- a/src/T2K/T2K_CC1pip_CH_XSec_1DQ2_nu.cxx +++ b/src/T2K/T2K_CC1pip_CH_XSec_1DQ2_nu.cxx @@ -1,204 +1,201 @@ #include #include "T2K_SignalDef.h" #include "T2K_CC1pip_CH_XSec_1DQ2_nu.h" // The constructor T2K_CC1pip_CH_XSec_1DQ2_nu::T2K_CC1pip_CH_XSec_1DQ2_nu(std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile){ EnuMin = 0.; EnuMax = 100.; fIsDiag = false; // Here we can give either MB (kMB), extended MB (keMB) or Delta (kDelta) if (type.find("eMB") != std::string::npos) { fT2KSampleType = keMB; fName = "T2K_CC1pip_CH_XSec_1DQ2eMB_nu"; fPlotTitles = "; Q^{2}_{eMB} (GeV^{2}); d#sigma/dQ^{2}_{eMB} (cm^{2}/GeV^{2}/nucleon)"; } else if (type.find("MB") != std::string::npos) { fT2KSampleType = kMB; fName = "T2K_CC1pip_CH_XSec_1DQ2MB_nu"; fPlotTitles = "; Q^{2}_{MB} (GeV^{2}); d#sigma/dQ^{2}_{MB} (cm^{2}/GeV^{2}/nucleon)"; } else if (type.find("Delta") != std::string::npos) { fT2KSampleType = kDelta; fName = "T2K_CC1pip_CH_XSec_1DQ2delta_nu"; fPlotTitles = "; Q^{2}_{#Delta} (GeV^{2}); d#sigma/dQ^{2}_{#Delta} (cm^{2}/GeV^{2}/nucleon)"; } else { LOG(SAM) << "Found no specified type, using MiniBooNE E_nu/Q2 definition" << std::endl; fT2KSampleType = kMB; fName = "T2K_CC1pip_CH_XSec_1DQ2MB_nu"; fPlotTitles = "; Q^{2}_{MB} (GeV^{2}); d#sigma/dQ^{2}_{MB} (cm^{2}/GeV^{2}/nucleon)"; } Measurement1D::SetupMeasurement(inputfile, type, rw, fakeDataFile); //type = keMB; //type = kDelta; if (fT2KSampleType == kMB) { this->SetDataValues(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/Q2_MB.root"); this->SetCovarMatrix(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/Q2_MB.root"); } else if (fT2KSampleType == keMB) { this->SetDataValues(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/Q2_extendedMB.root"); this->SetCovarMatrix(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/Q2_extendedMB.root"); } else if (fT2KSampleType == kDelta) { this->SetDataValues(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/Q2_Delta.root"); this->SetCovarMatrix(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/Q2_Delta.root"); } else { ERR(FTL) << "No data type set for " << fName << std::endl; ERR(FTL) << __FILE__ << ":" << __LINE__ << std::endl; exit(-1); } + SetShapeCovar(); this->SetupDefaultHist(); - + this->fScaleFactor = (GetEventHistogram()->Integral("width")*1E-38)/double(fNEvents)/TotalIntegratedFlux("width"); }; // Override this for now // Should really have Measurement1D do this properly though void T2K_CC1pip_CH_XSec_1DQ2_nu::SetDataValues(std::string fileLocation) { LOG(DEB) << "Reading: " << this->fName << "\nData: " << fileLocation.c_str() << std::endl; TFile *dataFile = new TFile(fileLocation.c_str()); //truly great .root file! // Don't want the last bin of dataCopy TH1D *dataCopy = (TH1D*)(dataFile->Get("hResult_sliced_0_1"))->Clone(); const int nPoints = dataCopy->GetNbinsX()-6; LOG(DEB) << nPoints << std::endl; double *binEdges = new double[nPoints+1]; for (int i = 0; i < nPoints+1; i++) { binEdges[i] = dataCopy->GetBinLowEdge(i+1); } for (int i = 0; i < nPoints+1; i++) { LOG(DEB) << "binEdges[" << i << "] = " << binEdges[i] << std::endl; } fDataHist = new TH1D((fName+"_data").c_str(), (fName+"_data"+fPlotTitles).c_str(), nPoints, binEdges); for (int i = 0; i < fDataHist->GetNbinsX(); i++) { fDataHist->SetBinContent(i+1, dataCopy->GetBinContent(i+1)*1E-38); fDataHist->SetBinError(i+1, dataCopy->GetBinError(i+1)*1E-38); LOG(DEB) << fDataHist->GetBinLowEdge(i+1) << " " << fDataHist->GetBinContent(i+1) << " " << fDataHist->GetBinError(i+1) << std::endl; } fDataHist->SetDirectory(0); //should disassociate fDataHist with dataFile fDataHist->SetNameTitle((fName+"_data").c_str(), (fName+"_MC"+fPlotTitles).c_str()); dataFile->Close(); }; // Override this for now // Should really have Measurement1D do this properly though void T2K_CC1pip_CH_XSec_1DQ2_nu::SetCovarMatrix(std::string fileLocation) { LOG(DEB) << "Covariance: " << fileLocation.c_str() << std::endl; TFile *dataFile = new TFile(fileLocation.c_str()); //truly great .root file! TH2D *covarMatrix = (TH2D*)(dataFile->Get("TMatrixDBase;1"))->Clone(); int nBinsX = covarMatrix->GetXaxis()->GetNbins(); int nBinsY = covarMatrix->GetYaxis()->GetNbins(); LOG(DEB) << nBinsX << std::endl; LOG(DEB) << fDataHist->GetNbinsX() << std::endl; if ((nBinsX != nBinsY)) ERR(WRN) << "covariance matrix not square!" << std::endl; - this->covar = new TMatrixDSym(nBinsX-7); this->fFullCovar = new TMatrixDSym(nBinsX-7); // First two entries are BS // Last entry is BS for (int i = 0; i < nBinsX-7; i++) { for (int j = 0; j < nBinsY-7; j++) { - (*this->covar)(i, j) = covarMatrix->GetBinContent(i+3, j+3); //adds syst+stat covariances (*this->fFullCovar)(i, j) = covarMatrix->GetBinContent(i+3, j+3); //adds syst+stat covariances - LOG(DEB) << "covar(" << i << ", " << j << ") = " << (*this->covar)(i,j) << std::endl; + LOG(DEB) << "fFullCovar(" << i << ", " << j << ") = " << (*this->fFullCovar)(i,j) << std::endl; } } //should now have set covariance, I hope - TDecompChol tempMat = TDecompChol(*this->covar); - this->covar = new TMatrixDSym(nBinsX, tempMat.Invert().GetMatrixArray(), ""); - // *this->covar *= 1E-38*1E-38; - + this->fDecomp = StatUtils::GetDecomp(this->fFullCovar); + this->covar = StatUtils::GetInvert(this->fFullCovar); return; }; void T2K_CC1pip_CH_XSec_1DQ2_nu::FillEventVariables(FitEvent *event) { if (event->NumFSParticle(13) == 0 || event->NumFSParticle(211) == 0) return; TLorentzVector Pnu = event->GetNeutrinoIn()->fP; TLorentzVector Ppip = event->GetHMFSParticle(211)->fP; TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; double q2 = -999; switch(fT2KSampleType) { // First int refers to how we reconstruct Enu // 0 uses true neutrino energy (not used here but common for other analyses when they unfold to true Enu from reconstructed Enu) // 1 uses "extended MiniBooNE" method // 2 uses "MiniBooNE reconstructed" method // 3 uses Delta resonance mass for reconstruction // // The last bool refers to if pion directional information was used // // Use MiniBooNE reconstructed Enu; uses Michel tag so no pion direction information case kMB: q2 = FitUtils::Q2CC1piprec(Pnu, Pmu, Ppip, 2, false); break; // Use Extended MiniBooNE reconstructed Enu // Needs pion information to reconstruct so bool is true (did not include Michel e tag) case keMB: q2 = FitUtils::Q2CC1piprec(Pnu, Pmu, Ppip, 1, true); break; // Use Delta resonance reconstructed Enu // Uses Michel electron so don't have pion directional information case kDelta: q2 = FitUtils::Q2CC1piprec(Pnu, Pmu, Ppip, 3, false); break; } fXVar = q2; return; }; //******************************************************************** bool T2K_CC1pip_CH_XSec_1DQ2_nu::isSignal(FitEvent *event) { //******************************************************************** // Warning: The CH analysis has different signal definition to the H2O analysis! // Often to do with the Michel tag switch(fT2KSampleType) { // Using MiniBooNE formula for Enu reconstruction on the Q2 variable // Does have Michel e tag, set bool to true! case kMB: return SignalDef::isCC1pip_T2K_CH(event, EnuMin, EnuMax, true); break; // Using extended MiniBooNE formula for Enu reconstruction on the Q2 variable // Does not have Michel e tag because we need directional information to reconstruct Q2 case keMB: return SignalDef::isCC1pip_T2K_CH(event, EnuMin, EnuMax, false); break; // Using Delta resonance for Enu reconstruction on the Q2 variable // Does have Michel e tag, bool to true case kDelta: return SignalDef::isCC1pip_T2K_CH(event, EnuMin, EnuMax, true); break; } // Default to return false return false; } diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1DWrec_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1DWrec_nu.cxx index 0961a84..c31b6aa 100644 --- a/src/T2K/T2K_CC1pip_CH_XSec_1DWrec_nu.cxx +++ b/src/T2K/T2K_CC1pip_CH_XSec_1DWrec_nu.cxx @@ -1,118 +1,116 @@ #include #include "T2K_SignalDef.h" #include "T2K_CC1pip_CH_XSec_1DWrec_nu.h" // The constructor T2K_CC1pip_CH_XSec_1DWrec_nu::T2K_CC1pip_CH_XSec_1DWrec_nu(std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile){ fName = "T2K_CC1pip_CH_XSec_1DWrec_nu"; fPlotTitles = "; W_{rec} (GeV/c); d#sigma/dW_{rec} (cm^{2}/(GeV/c^{2})/nucleon)"; EnuMin = 0.; EnuMax = 100.; fIsDiag = false; Measurement1D::SetupMeasurement(inputfile, type, rw, fakeDataFile); this->SetDataValues(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/W.root"); this->SetCovarMatrix(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/W.root"); - + SetShapeCovar(); this->SetupDefaultHist(); this->fScaleFactor = (GetEventHistogram()->Integral("width")*1E-38)/double(fNEvents)/TotalIntegratedFlux("width"); }; // Override this for now // Should really have Measurement1D do this properly though void T2K_CC1pip_CH_XSec_1DWrec_nu::SetDataValues(std::string fileLocation) { LOG(DEB) << "Reading: " << this->fName << "\nData: " << fileLocation.c_str() << std::endl; TFile *dataFile = new TFile(fileLocation.c_str()); //truly great .root file! // Don't want the first and last bin of dataCopy TH1D *dataCopy = (TH1D*)(dataFile->Get("hResult_sliced_0_1"))->Clone(); LOG(DEB) << dataCopy->GetNbinsX() << std::endl; const int dataPoints = dataCopy->GetNbinsX()-2; double *binEdges = new double[dataPoints+1]; // Want to skip the first bin here for (int i = 0; i < dataPoints+1; i++) { binEdges[i] = dataCopy->GetBinLowEdge(i+2); } for (int i = 0; i < dataPoints+1; i++) { LOG(DEB) << "binEdges[" << i << "] = " << binEdges[i] << std::endl; } fDataHist = new TH1D((fName+"_data").c_str(), (fName+"_data"+fPlotTitles).c_str(), dataPoints, binEdges); for (int i = 0; i < fDataHist->GetNbinsX(); i++) { fDataHist->SetBinContent(i+1, dataCopy->GetBinContent(i+2)*1E-38); fDataHist->SetBinError(i+1, dataCopy->GetBinError(i+2)*1E-38); LOG(DEB) << fDataHist->GetBinLowEdge(i+1) << " " << fDataHist->GetBinContent(i+1) << " " << fDataHist->GetBinError(i+1) << std::endl; } fDataHist->SetDirectory(0); //should disassociate fDataHist with dataFile //fDataHist->SetNameTitle((fName+"_data").c_str(), (fName+"_MC"+fPlotTitles).c_str()); dataFile->Close(); }; // Override this for now // Should really have Measurement1D do this properly though void T2K_CC1pip_CH_XSec_1DWrec_nu::SetCovarMatrix(std::string fileLocation) { LOG(DEB) << "Covariance: " << fileLocation.c_str() << std::endl; TFile *dataFile = new TFile(fileLocation.c_str()); //truly great .root file! TH2D *covarMatrix = (TH2D*)(dataFile->Get("TMatrixDBase;1"))->Clone(); const int nBinsX = covarMatrix->GetXaxis()->GetNbins(); const int nBinsY = covarMatrix->GetYaxis()->GetNbins(); if ((nBinsX != nBinsY)) ERR(WRN) << "covariance matrix not square!" << std::endl; LOG(DEB) << nBinsX << std::endl; LOG(DEB) << fDataHist->GetNbinsX() << std::endl; - this->covar = new TMatrixDSym(nBinsX-3); this->fFullCovar = new TMatrixDSym(nBinsX-3); for (int i = 0; i < nBinsX-3; i++) { for (int j = 0; j < nBinsY-3; j++) { - (*this->covar)(i, j) = covarMatrix->GetBinContent(i+4, j+4); //adds syst+stat covariances - LOG(DEB) << "covar(" << i << ", " << j << ") = " << (*this->covar)(i,j) << std::endl; + (*this->fFullCovar)(i, j) = covarMatrix->GetBinContent(i+4, j+4); //adds syst+stat covariances + LOG(DEB) << "fFullCovar(" << i << ", " << j << ") = " << (*this->fFullCovar)(i,j) << std::endl; } } //should now have set covariance, I hope - TDecompChol tempMat = TDecompChol(*this->covar); - this->covar = new TMatrixDSym(nBinsX, tempMat.Invert().GetMatrixArray(), ""); - // *this->covar *= 1E-38*1E-38; + this->fDecomp = StatUtils::GetDecomp(this->fFullCovar); + this->covar = StatUtils::GetInvert(this->fFullCovar); return; }; void T2K_CC1pip_CH_XSec_1DWrec_nu::FillEventVariables(FitEvent *event) { if (event->NumFSParticle(13) == 0 || event->NumFSParticle(211) == 0) return; TLorentzVector Pnu = event->GetNeutrinoIn()->fP; TLorentzVector Ppip = event->GetHMFSParticle(211)->fP; TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; double Wrec = FitUtils::WrecCC1pip_T2K_MB(Pnu, Pmu, Ppip); fXVar = Wrec; return; }; //******************************************************************** bool T2K_CC1pip_CH_XSec_1DWrec_nu::isSignal(FitEvent *event) { //******************************************************************** // This sample includes the Michel e tag so do not have to cut into the pion phase space return SignalDef::isCC1pip_T2K_CH(event, EnuMin, EnuMax, true); } diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1Dpmu_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1Dpmu_nu.cxx index b56e2b7..702a16d 100644 --- a/src/T2K/T2K_CC1pip_CH_XSec_1Dpmu_nu.cxx +++ b/src/T2K/T2K_CC1pip_CH_XSec_1Dpmu_nu.cxx @@ -1,156 +1,154 @@ #include "T2K_CC1pip_CH_XSec_1Dpmu_nu.h" #include //******************************************************************** T2K_CC1pip_CH_XSec_1Dpmu_nu::T2K_CC1pip_CH_XSec_1Dpmu_nu(nuiskey samplekey) { //******************************************************************** // Sample overview --------------------------------------------------- std::string descrip = "T2K_CC1pip_CH_XSec_1Dpmu_nu sample. \n" \ "Target: CH \n" \ "Flux: T2k Forward Horn Current nue + nuebar \n" \ "Signal: Any event with 1 electron, any nucleons, and no other FS particles \n"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetDescription(descrip); fSettings.SetXTitle("p_{#mu} (GeV/c)"); fSettings.SetYTitle("d#sigma/dp_{#mu} (cm^{2}/(GeV/c)/nucleon)"); fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG"); fSettings.SetEnuRange(0.0, 100.0); fSettings.DefineAllowedTargets("C,H"); // CCQELike plot information fSettings.SetTitle("T2K_CC1pip_CH_XSec_1Dpmu_nu"); fSettings.SetDataInput( GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/Pmu.root" ); fSettings.SetCovarInput( GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/Pmu.root" ); fSettings.DefineAllowedSpecies("numu"); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = (GetEventHistogram()->Integral("width")*1E-38)/double(fNEvents)/TotalIntegratedFlux("width"); // Plot Setup ------------------------------------------------------- SetDataValues( fSettings.GetDataInput() ); SetCovarMatrix( fSettings.GetCovarInput() ); + SetShapeCovar(); + // Final setup --------------------------------------------------- FinaliseMeasurement(); }; // Override this for now // Should really have Measurement1D do this properly though void T2K_CC1pip_CH_XSec_1Dpmu_nu::SetDataValues(std::string fileLocation) { LOG(DEB) << "Reading: " << this->fName << "\nData: " << fileLocation.c_str() << std::endl; TFile *dataFile = new TFile(fileLocation.c_str()); //truly great .root file! // Don't want the last bin of dataCopy TH1D *dataCopy = (TH1D*)(dataFile->Get("hResult_sliced_0_1"))->Clone(); LOG(DEB) << "nbins = " << dataCopy->GetNbinsX() << std::endl; double *binEdges = new double[dataCopy->GetNbinsX()-1]; for (int i = 0; i < dataCopy->GetNbinsX()-1; i++) { binEdges[i] = dataCopy->GetBinLowEdge(i+1); } binEdges[dataCopy->GetNbinsX()-1] = dataCopy->GetBinLowEdge(dataCopy->GetNbinsX()); for (int i = 0; i < dataCopy->GetNbinsX()+5; i++) { LOG(DEB) << "binEdges[" << i << "] = " << binEdges[i] << std::endl; } fDataHist = new TH1D((fName+"_data").c_str(), (fName+"_data"+fPlotTitles).c_str(), dataCopy->GetNbinsX()-1, binEdges); for (int i = 0; i < fDataHist->GetNbinsX(); i++) { fDataHist->SetBinContent(i+1, dataCopy->GetBinContent(i+1)*1E-38); fDataHist->SetBinError(i+1, dataCopy->GetBinError(i+1)*1E-38); LOG(DEB) << fDataHist->GetBinLowEdge(i+1) << " " << fDataHist->GetBinContent(i+1) << std::endl; } fDataHist->SetDirectory(0); //should disassociate fDataHist with dataFile fDataHist->SetNameTitle((fName+"_data").c_str(), (fName+"_MC"+fPlotTitles).c_str()); dataFile->Close(); }; // Override this for now // Should really have Measurement1D do this properly though void T2K_CC1pip_CH_XSec_1Dpmu_nu::SetCovarMatrix(std::string fileLocation) { LOG(DEB) << "Covariance: " << fileLocation.c_str() << std::endl; TFile *dataFile = new TFile(fileLocation.c_str()); //truly great .root file! TH2D *covarMatrix = (TH2D*)(dataFile->Get("TMatrixDBase;1"))->Clone(); int nBinsX = covarMatrix->GetXaxis()->GetNbins(); int nBinsY = covarMatrix->GetYaxis()->GetNbins(); if ((nBinsX != nBinsY)) ERR(WRN) << "covariance matrix not square!" << std::endl; - this->covar = new TMatrixDSym(nBinsX-2); this->fFullCovar = new TMatrixDSym(nBinsX-2); // First two entries are BS // Last entry is BS for (int i = 1; i < nBinsX-1; i++) { for (int j = 1; j < nBinsY-1; j++) { - (*this->covar)(i-1, j-1) = covarMatrix->GetBinContent(i+1, j+1); //adds syst+stat covariances (*this->fFullCovar)(i-1, j-1) = covarMatrix->GetBinContent(i+1, j+1); //adds syst+stat covariances - LOG(DEB) << "covar(" << i-1 << ", " << j-1 << ") = " << (*this->covar)(i-1,j-1) << std::endl; + LOG(DEB) << "fFullCovar(" << i-1 << ", " << j-1 << ") = " << (*this->fFullCovar)(i-1,j-1) << std::endl; } } //should now have set covariance, I hope - TDecompChol tempMat = TDecompChol(*this->covar); - this->covar = new TMatrixDSym(nBinsX, tempMat.Invert().GetMatrixArray(), ""); - // *this->covar *= 1E-38*1E-38; - + this->fDecomp = StatUtils::GetDecomp(this->fFullCovar); + this->covar = StatUtils::GetInvert(this->fFullCovar); return; }; void T2K_CC1pip_CH_XSec_1Dpmu_nu::FillEventVariables(FitEvent *event) { if (event->NumFSParticle(13) == 0 || event->NumFSParticle(211) == 0) return; TLorentzVector Pnu = event->GetNeutrinoIn()->fP; TLorentzVector Ppip = event->GetHMFSParticle(211)->fP; TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; double pmu = FitUtils::p(Pmu); fXVar = pmu; return; }; //******************************************************************** bool T2K_CC1pip_CH_XSec_1Dpmu_nu::isSignal(FitEvent *event) { //******************************************************************** // Warning: The CH analysis has different signal definition to the H2O analysis! // // If Michel e- is used for pion PID we don't have directional info on pion; set the bool to true // The bool is if we use Michel e- or not // Also, for events binned in muon momentum/angle there's no cut on the pion kinematics // // Additionally, the 2D distribution uses 0.8 > cos th mu > 0 and no pion phase space reduction applied // Also no muon momentum reduction applied // // This uses a slightly custom signal definition where a cut is only placed on the muon angle, not the momentum if (!SignalDef::isCC1pi(event, 14, 211, EnuMin, EnuMax)) return false; TLorentzVector Pnu = event->GetHMISParticle(14)->fP; TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; double cos_th_mu = cos(FitUtils::th(Pnu,Pmu)); if (cos_th_mu >= 0.2) return true; return false; }; diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1Dppi_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1Dppi_nu.cxx index 423d9c1..f787b76 100644 --- a/src/T2K/T2K_CC1pip_CH_XSec_1Dppi_nu.cxx +++ b/src/T2K/T2K_CC1pip_CH_XSec_1Dppi_nu.cxx @@ -1,179 +1,177 @@ #include #include "T2K_SignalDef.h" #include "T2K_CC1pip_CH_XSec_1Dppi_nu.h" //******************************************************************** T2K_CC1pip_CH_XSec_1Dppi_nu::T2K_CC1pip_CH_XSec_1Dppi_nu(nuiskey samplekey) { //******************************************************************** // Sample overview --------------------------------------------------- std::string descrip = "T2K_CC1pip_CH_XSec_1Dppi_nu sample. \n" \ "Target: CH \n" \ "Flux: T2k Forward Horn Current nue + nuebar \n" \ "Signal: Any event with 1 electron, any nucleons, and no other FS particles \n"; // Setup common settings fSettings = LoadSampleSettings(samplekey); fSettings.SetTitle("T2K_CC1pip_CH_XSec_1Dppi_nu"); fSettings.SetDescription(descrip); fSettings.SetXTitle("p_{#pi} (GeV/c)"); fSettings.SetYTitle("d#sigma/dW_{rec} (cm^{2}/(GeV/c)/nucleon)"); fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/DIAG"); fSettings.SetEnuRange(0.0, 100.0); fSettings.DefineAllowedTargets("C,H"); fSettings.DefineAllowedSpecies("numu"); FinaliseSampleSettings(); // Scaling Setup --------------------------------------------------- // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38) / double(fNEvents) / TotalIntegratedFlux("width"); // Plot Setup ------------------------------------------------------- if (fSettings.GetS("type").find("Michel") != std::string::npos) { useMichel = true; fName += "_Michel"; SetDataValues(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/Ppi.root"); SetCovarMatrix(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/Ppi.root"); } else { useMichel = false; - fName += "_kin"; + // fName += "_kin"; SetDataValues(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/Ppi_noME.root"); SetCovarMatrix(GeneralUtils::GetTopLevelDir() + "/data/T2K/CC1pip/CH/Ppi_noME.root"); } + + SetShapeCovar(); // Final setup --------------------------------------------------- FinaliseMeasurement(); }; // Override this for now // Should really have Measurement1D do this properly though void T2K_CC1pip_CH_XSec_1Dppi_nu::SetDataValues(std::string fileLocation) { LOG(DEB) << "Reading: " << this->fName << "\nData: " << fileLocation.c_str() << std::endl; TFile *dataFile = new TFile(fileLocation.c_str()); //truly great .root file! // Don't want the last bin of dataCopy TH1D *dataCopy = (TH1D*)(dataFile->Get("hResult_sliced_0_1"))->Clone(); LOG(DEB) << dataCopy->GetNbinsX() << std::endl; double *binEdges = new double[dataCopy->GetNbinsX() - 1]; for (int i = 0; i < dataCopy->GetNbinsX() - 1; i++) { binEdges[i] = dataCopy->GetBinLowEdge(i + 1); } binEdges[dataCopy->GetNbinsX() - 1] = dataCopy->GetBinLowEdge(dataCopy->GetNbinsX()); fDataHist = new TH1D((fName + "_data").c_str(), (fName + "_data" + fPlotTitles).c_str(), dataCopy->GetNbinsX() - 2, binEdges); for (int i = 0; i < fDataHist->GetNbinsX(); i++) { fDataHist->SetBinContent(i + 1, dataCopy->GetBinContent(i + 1) * 1E-38); fDataHist->SetBinError(i + 1, dataCopy->GetBinError(i + 1) * 1E-38); LOG(DEB) << fDataHist->GetBinLowEdge(i + 1) << " " << fDataHist->GetBinContent(i + 1) << std::endl; } fDataHist->SetDirectory(0); //should disassociate fDataHist with dataFile fDataHist->SetNameTitle((fName + "_data").c_str(), (fName + "_MC" + fPlotTitles).c_str()); dataFile->Close(); }; // Override this for now // Should really have Measurement1D do this properly though void T2K_CC1pip_CH_XSec_1Dppi_nu::SetCovarMatrix(std::string fileLocation) { LOG(DEB) << "Covariance: " << fileLocation.c_str() << std::endl; TFile *dataFile = new TFile(fileLocation.c_str()); //truly great .root file! TH2D *covarMatrix = (TH2D*)(dataFile->Get("TMatrixDBase;1"))->Clone(); int nBinsX = covarMatrix->GetXaxis()->GetNbins(); int nBinsY = covarMatrix->GetYaxis()->GetNbins(); if ((nBinsX != nBinsY)) ERR(WRN) << "covariance matrix not square!" << std::endl; - this->covar = new TMatrixDSym(nBinsX - 2); - this->fFullCovar = new TMatrixDSym(nBinsX - 2); + this->fFullCovar = new TMatrixDSym(nBinsX - 3); // First two entries are BS // Last entry is BS - for (int i = 1; i < nBinsX - 1; i++) { - for (int j = 1; j < nBinsY - 1; j++) { + for (int i = 2; i < nBinsX-1; i++) { + for (int j = 2; j < nBinsY-1; j++) { LOG(DEB) << "(" << i << ", " << j << ") = " << covarMatrix->GetBinContent(i + 1, j + 1) << std::endl; - (*this->covar)(i - 1, j - 1) = covarMatrix->GetBinContent(i, j); //adds syst+stat covariances - (*this->fFullCovar)(i - 1, j - 1) = covarMatrix->GetBinContent(i, j); //adds syst+stat covariances + (*this->fFullCovar)(i - 2, j - 2) = covarMatrix->GetBinContent(i, j); //adds syst+stat covariances } } //should now have set covariance, I hope - TDecompChol tempMat = TDecompChol(*this->covar); - this->covar = new TMatrixDSym(nBinsX, tempMat.Invert().GetMatrixArray(), ""); - // *this->covar *= 1E-38*1E-38; - + this->fDecomp = StatUtils::GetDecomp(this->fFullCovar); + this->covar = StatUtils::GetInvert(this->fFullCovar); return; }; void T2K_CC1pip_CH_XSec_1Dppi_nu::FillEventVariables(FitEvent *event) { if (event->NumFSParticle(13) == 0 || event->NumFSParticle(211) == 0) return; TLorentzVector Pnu = event->GetNeutrinoIn()->fP; TLorentzVector Ppip = event->GetHMFSParticle(211)->fP; TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; double ppip = FitUtils::p(Ppip); fXVar = ppip; return; }; //******************************************************************** bool T2K_CC1pip_CH_XSec_1Dppi_nu::isSignal(FitEvent *event) { //******************************************************************** // This distribution uses a somewhat different signal definition so might as well implement it separately here // If we use Michel tag sample we don't cut into the pion phase space, only the muon phase space // The last bool refers to if we have Michel e or not if (useMichel) { return SignalDef::isCC1pip_T2K_CH(event, EnuMin, EnuMax, true); } else { // Custom signal definition if we aren't using Michel tag; cut on muon and cut only on pion angle // does the event pass the muon cut? bool muonPass = SignalDef::isCC1pip_T2K_CH(event, EnuMin, EnuMax, true); // If the event doesn't pass the muon cut return false if (!muonPass) { return false; } // does the event pass the pion angle cut? // we already know there's just one muon in the event if it passes muonPass so don't need to make an event loop rejection // Need the neutrino four-vector to get the angle between pion and neutrino TLorentzVector Pnu = event->PartInfo(0)->fP; TLorentzVector Ppip; for (unsigned int j = 2; j < event->Npart(); j++) { if (!((event->PartInfo(j))->fIsAlive) && (event->PartInfo(j))->fNEUTStatusCode != 0) continue; //move on if NOT ALIVE and NOT NORMAL int PID = (event->PartInfo(j))->fPID; if (PID == 211) { Ppip = event->PartInfo(j)->fP; // Once the pion is found we can break break; } } double cos_th_pi = cos(FitUtils::th(Pnu, Ppip)); // Now check the angle of the pion if (cos_th_pi <= 0.2) { return false; } else { return true; } } // Unnecessary default to false return false; } diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1Dq3_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1Dq3_nu.cxx index 28f6e7b..6a34455 100644 --- a/src/T2K/T2K_CC1pip_CH_XSec_1Dq3_nu.cxx +++ b/src/T2K/T2K_CC1pip_CH_XSec_1Dq3_nu.cxx @@ -1,113 +1,109 @@ #include #include "T2K_SignalDef.h" #include "T2K_CC1pip_CH_XSec_1Dq3_nu.h" // The constructor T2K_CC1pip_CH_XSec_1Dq3_nu::T2K_CC1pip_CH_XSec_1Dq3_nu(std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile){ fName = "T2K_CC1pip_CH_XSec_1Dq3_nu"; fPlotTitles = "; q_{3} (GeV/c); d#sigma/dq_{3} (cm^{2}/(GeV/c)/nucleon)"; EnuMin = 0.; EnuMax = 100.; fIsDiag = false; Measurement1D::SetupMeasurement(inputfile, type, rw, fakeDataFile); this->SetDataValues(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/Q3.root"); this->SetCovarMatrix(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/Q3.root"); - + SetShapeCovar(); this->SetupDefaultHist(); this->fScaleFactor = (GetEventHistogram()->Integral("width")*1E-38)/double(fNEvents)/TotalIntegratedFlux("width"); }; // Override this for now // Should really have Measurement1D do this properly though void T2K_CC1pip_CH_XSec_1Dq3_nu::SetDataValues(std::string fileLocation) { LOG(DEB) << "Reading: " << this->fName << "\nData: " << fileLocation.c_str() << std::endl; TFile *dataFile = new TFile(fileLocation.c_str()); //truly great .root file! // Don't want the last bin of dataCopy TH1D *dataCopy = (TH1D*)(dataFile->Get("hResult_sliced_0_1"))->Clone(); - double *binEdges = new double[dataCopy->GetNbinsX()-1]; + double *binEdges = new double[dataCopy->GetNbinsX()-2]; LOG(DEB) << dataCopy->GetNbinsX() << std::endl; - for (int i = 1; i < dataCopy->GetNbinsX(); i++) { + for (int i = 1; i < dataCopy->GetNbinsX()-1; i++) { binEdges[i-1] = dataCopy->GetBinLowEdge(i+1); LOG(DEB) << i-1 << " " << binEdges[i-1] << " from binLowEdge " << i+1 << std::endl; } fDataHist = new TH1D((fName+"_data").c_str(), (fName+"_data"+fPlotTitles).c_str(), dataCopy->GetNbinsX()-2, binEdges); for (int i = 0; i < fDataHist->GetNbinsX(); i++) { fDataHist->SetBinContent(i+1, dataCopy->GetBinContent(i+2)*1E-38); fDataHist->SetBinError(i+1, dataCopy->GetBinError(i+2)*1E-38); LOG(DEB) << fDataHist->GetBinLowEdge(i+1) << " " << fDataHist->GetBinContent(i+1) << " " << fDataHist->GetBinError(i+1) << std::endl; } fDataHist->SetDirectory(0); //should disassociate fDataHist with dataFile fDataHist->SetNameTitle((fName+"_data").c_str(), (fName+"_MC"+fPlotTitles).c_str()); dataFile->Close(); }; // Override this for now // Should really have Measurement1D do this properly though void T2K_CC1pip_CH_XSec_1Dq3_nu::SetCovarMatrix(std::string fileLocation) { LOG(DEB) << "Covariance: " << fileLocation.c_str() << std::endl; TFile *dataFile = new TFile(fileLocation.c_str()); //truly great .root file! TH2D *covarMatrix = (TH2D*)(dataFile->Get("TMatrixDBase;1"))->Clone(); int nBinsX = covarMatrix->GetXaxis()->GetNbins(); int nBinsY = covarMatrix->GetYaxis()->GetNbins(); if ((nBinsX != nBinsY)) ERR(WRN) << "covariance matrix not square!" << std::endl; - this->covar = new TMatrixDSym(nBinsX-2); - this->fFullCovar = new TMatrixDSym(nBinsX-2); + this->fFullCovar = new TMatrixDSym(nBinsX-3); // First two entries are BS // Last entry is BS for (int i = 2; i < nBinsX-1; i++) { for (int j = 2; j < nBinsY-1; j++) { - (*this->covar)(i-2, j-2) = covarMatrix->GetBinContent(i+1, j+1); //adds syst+stat covariances (*this->fFullCovar)(i-2, j-2) = covarMatrix->GetBinContent(i+1, j+1); //adds syst+stat covariances - LOG(DEB) << "covar(" << i-2 << ", " << j-2 << ") = " << (*this->covar)(i-2,j-2) << std::endl; + LOG(DEB) << "fFullCovar(" << i-2 << ", " << j-2 << ") = " << (*this->fFullCovar)(i-2,j-2) << std::endl; } } //should now have set covariance, I hope - TDecompChol tempMat = TDecompChol(*this->covar); - this->covar = new TMatrixDSym(nBinsX, tempMat.Invert().GetMatrixArray(), ""); - // *this->covar *= 1E-38*1E-38; - + this->fDecomp = StatUtils::GetDecomp(this->fFullCovar); + this->covar = StatUtils::GetInvert(this->fFullCovar); return; }; void T2K_CC1pip_CH_XSec_1Dq3_nu::FillEventVariables(FitEvent *event) { if (event->NumFSParticle(13) == 0 || event->NumFSParticle(211) == 0) return; TLorentzVector Pnu = event->GetNeutrinoIn()->fP; TLorentzVector Ppip = event->GetHMFSParticle(211)->fP; TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; double q3 = FitUtils::q3_CC1pip_T2K(Pnu, Pmu, Ppip); fXVar = q3; return; }; //******************************************************************** bool T2K_CC1pip_CH_XSec_1Dq3_nu::isSignal(FitEvent *event) { //******************************************************************** // Has a Michel e sample in so no phase space cut on pion required return SignalDef::isCC1pip_T2K_CH(event, EnuMin, EnuMax, true); } diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1Dthmupi_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1Dthmupi_nu.cxx index b191482..da6cc55 100644 --- a/src/T2K/T2K_CC1pip_CH_XSec_1Dthmupi_nu.cxx +++ b/src/T2K/T2K_CC1pip_CH_XSec_1Dthmupi_nu.cxx @@ -1,113 +1,109 @@ #include #include "T2K_SignalDef.h" #include "T2K_CC1pip_CH_XSec_1Dthmupi_nu.h" // The constructor T2K_CC1pip_CH_XSec_1Dthmupi_nu::T2K_CC1pip_CH_XSec_1Dthmupi_nu(std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile){ fName = "T2K_CC1pip_CH_XSec_1Dthmupi_nu"; fPlotTitles = "; #theta_{#pi,#mu} (radians); d#sigma/d#theta_{#pi} (cm^{2}/nucleon)"; EnuMin = 0.; EnuMax = 100.; fIsDiag = false; Measurement1D::SetupMeasurement(inputfile, type, rw, fakeDataFile); this->SetDataValues(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/Thetapimu.root"); this->SetCovarMatrix(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/Thetapimu.root"); - + SetShapeCovar(); this->SetupDefaultHist(); this->fScaleFactor = (GetEventHistogram()->Integral("width")*1E-38)/double(fNEvents)/TotalIntegratedFlux("width"); }; // Override this for now // Should really have Measurement1D do this properly though void T2K_CC1pip_CH_XSec_1Dthmupi_nu::SetDataValues(std::string fileLocation) { LOG(DEB) << "Reading: " << this->fName << "\nData: " << fileLocation.c_str() << std::endl; TFile *dataFile = new TFile(fileLocation.c_str()); //truly great .root file! // Don't want the first and last bin of dataCopy TH1D *dataCopy = (TH1D*)(dataFile->Get("hResult_sliced_0_1"))->Clone(); LOG(DEB) << "dataCopy->GetNbinsX() = " << dataCopy->GetNbinsX() << std::endl; double *binEdges = new double[dataCopy->GetNbinsX()]; for (int i = 0; i < dataCopy->GetNbinsX(); i++) { binEdges[i] = dataCopy->GetBinLowEdge(i+1); } binEdges[dataCopy->GetNbinsX()] = dataCopy->GetBinLowEdge(dataCopy->GetNbinsX()+1); for (int i = 0; i < dataCopy->GetNbinsX()+5; i++) { LOG(DEB) << "binEdges[" << i << "] = " << binEdges[i] << std::endl; } fDataHist = new TH1D((fName+"_data").c_str(), (fName+"_data"+fPlotTitles).c_str(), dataCopy->GetNbinsX(), binEdges); for (int i = 0; i < fDataHist->GetNbinsX()+1; i++) { fDataHist->SetBinContent(i+1, dataCopy->GetBinContent(i+1)*1E-38); fDataHist->SetBinError(i+1, dataCopy->GetBinError(i+1)*1E-38); LOG(DEB) << fDataHist->GetBinLowEdge(i+1) << " " << fDataHist->GetBinContent(i+1) << " " << fDataHist->GetBinError(i+1) << std::endl; } fDataHist->SetDirectory(0); //should disassociate fDataHist with dataFile dataFile->Close(); }; // Override this for now // Should really have Measurement1D do this properly though void T2K_CC1pip_CH_XSec_1Dthmupi_nu::SetCovarMatrix(std::string fileLocation) { LOG(DEB) << "Covariance: " << fileLocation.c_str() << std::endl; TFile *dataFile = new TFile(fileLocation.c_str()); //truly great .root file! TH2D *covarMatrix = (TH2D*)(dataFile->Get("TMatrixDBase;1"))->Clone(); int nBinsX = covarMatrix->GetXaxis()->GetNbins(); int nBinsY = covarMatrix->GetYaxis()->GetNbins(); if ((nBinsX != nBinsY)) ERR(WRN) << "covariance matrix not square!" << std::endl; - this->covar = new TMatrixDSym(nBinsX-1); this->fFullCovar = new TMatrixDSym(nBinsX-1); for (int i = 1; i < nBinsX; i++) { for (int j = 1; j < nBinsY; j++) { - (*this->covar)(i-1, j-1) = covarMatrix->GetBinContent(i, j); //adds syst+stat covariances (*this->fFullCovar)(i-1, j-1) = covarMatrix->GetBinContent(i, j); //adds syst+stat covariances - LOG(DEB) << "covar(" << i-1 << ", " << j-1 << ") = " << (*this->covar)(i-1,j-1) << std::endl; + LOG(DEB) << "fFullCovar(" << i-1 << ", " << j-1 << ") = " << (*this->fFullCovar)(i-1,j-1) << std::endl; } } //should now have set covariance, I hope - TDecompChol tempMat = TDecompChol(*this->covar); - this->covar = new TMatrixDSym(nBinsX, tempMat.Invert().GetMatrixArray(), ""); - // *this->covar *= 1E-38*1E-38; - + this->fDecomp = StatUtils::GetDecomp(this->fFullCovar); + this->covar = StatUtils::GetInvert(this->fFullCovar); dataFile->Close(); }; void T2K_CC1pip_CH_XSec_1Dthmupi_nu::FillEventVariables(FitEvent *event) { if (event->NumFSParticle(13) == 0 || event->NumFSParticle(211) == 0) return; TLorentzVector Pnu = event->GetNeutrinoIn()->fP; TLorentzVector Ppip = event->GetHMFSParticle(211)->fP; TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; double thmupi = FitUtils::th(Pmu, Ppip); fXVar = thmupi; return; }; //******************************************************************** bool T2K_CC1pip_CH_XSec_1Dthmupi_nu::isSignal(FitEvent *event) { //******************************************************************** // This sample requires directional information on the pion so can't use Michel tag sample return SignalDef::isCC1pip_T2K_CH(event, EnuMin, EnuMax, false); } diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1Dthpi_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1Dthpi_nu.cxx index 4324fe6..b8be6c5 100644 --- a/src/T2K/T2K_CC1pip_CH_XSec_1Dthpi_nu.cxx +++ b/src/T2K/T2K_CC1pip_CH_XSec_1Dthpi_nu.cxx @@ -1,114 +1,110 @@ #include #include "T2K_SignalDef.h" #include "T2K_CC1pip_CH_XSec_1Dthpi_nu.h" // The constructor T2K_CC1pip_CH_XSec_1Dthpi_nu::T2K_CC1pip_CH_XSec_1Dthpi_nu(std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile){ fName = "T2K_CC1pip_CH_XSec_1Dthpi_nu"; fPlotTitles = "; #theta_{#pi} (radians); d#sigma/d#theta_{#pi} (cm^{2}/nucleon)"; EnuMin = 0.; EnuMax = 100.; fIsDiag = false; Measurement1D::SetupMeasurement(inputfile, type, rw, fakeDataFile); this->SetDataValues(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/Thetapi.root"); this->SetCovarMatrix(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/Thetapi.root"); - + SetShapeCovar(); this->SetupDefaultHist(); this->fScaleFactor = (GetEventHistogram()->Integral("width")*1E-38)/double(fNEvents)/TotalIntegratedFlux("width"); }; // Override this for now // Should really have Measurement1D do this properly though void T2K_CC1pip_CH_XSec_1Dthpi_nu::SetDataValues(std::string fileLocation) { LOG(DEB) << "Reading: " << this->fName << "\nData: " << fileLocation.c_str() << std::endl; TFile *dataFile = new TFile(fileLocation.c_str()); //truly great .root file! // Don't want the first and last bin of dataCopy TH1D *dataCopy = (TH1D*)(dataFile->Get("hResult_sliced_0_1"))->Clone(); LOG(DEB) << "dataCopy->GetNbinsX() = " << dataCopy->GetNbinsX() << std::endl; double *binEdges = new double[dataCopy->GetNbinsX()-4]; for (int i = 0; i < dataCopy->GetNbinsX()-4; i++) { binEdges[i] = dataCopy->GetBinLowEdge(i+1); } binEdges[dataCopy->GetNbinsX()-4] = dataCopy->GetBinLowEdge(dataCopy->GetNbinsX()-3); for (int i = 0; i < dataCopy->GetNbinsX(); i++) { LOG(DEB) << "binEdges[" << i << "] = " << binEdges[i] << std::endl; } fDataHist = new TH1D((fName+"_data").c_str(), (fName+"_data"+fPlotTitles).c_str(), dataCopy->GetNbinsX()-4, binEdges); for (int i = 0; i < fDataHist->GetNbinsX(); i++) { fDataHist->SetBinContent(i+1, dataCopy->GetBinContent(i+1)*1E-38); fDataHist->SetBinError(i+1, dataCopy->GetBinError(i+1)*1E-38); LOG(DEB) << fDataHist->GetBinLowEdge(i+1) << " " << fDataHist->GetBinContent(i+1) << " " << fDataHist->GetBinError(i+1) << std::endl; } fDataHist->SetDirectory(0); //should disassociate fDataHist with dataFile dataFile->Close(); }; // Override this for now // Should really have Measurement1D do this properly though void T2K_CC1pip_CH_XSec_1Dthpi_nu::SetCovarMatrix(std::string fileLocation) { LOG(DEB) << "Covariance: " << fileLocation.c_str() << std::endl; TFile *dataFile = new TFile(fileLocation.c_str()); //truly great .root file! TH2D *covarMatrix = (TH2D*)(dataFile->Get("TMatrixDBase;1"))->Clone(); int nBinsX = covarMatrix->GetXaxis()->GetNbins(); int nBinsY = covarMatrix->GetYaxis()->GetNbins(); if ((nBinsX != nBinsY)) ERR(WRN) << "covariance matrix not square!" << std::endl; - this->covar = new TMatrixDSym(nBinsX-5); this->fFullCovar = new TMatrixDSym(nBinsX-5); for (int i = 2; i < nBinsX-3; i++) { for (int j = 2; j < nBinsY-3; j++) { - (*this->covar)(i-2, j-2) = covarMatrix->GetBinContent(i, j); //adds syst+stat covariances (*this->fFullCovar)(i-2, j-2) = covarMatrix->GetBinContent(i, j); //adds syst+stat covariances - LOG(DEB) << "covar(" << i-2 << ", " << j-2 << ") = " << (*this->covar)(i-2,j-2) << std::endl; + LOG(DEB) << "fFullCovar(" << i-2 << ", " << j-2 << ") = " << (*this->fFullCovar)(i-2,j-2) << std::endl; } } //should now have set covariance, I hope - TDecompChol tempMat = TDecompChol(*this->covar); - this->covar = new TMatrixDSym(nBinsX, tempMat.Invert().GetMatrixArray(), ""); - // *this->covar *= 1E-38*1E-38; - + this->fDecomp = StatUtils::GetDecomp(this->fFullCovar); + this->covar = StatUtils::GetInvert(this->fFullCovar); dataFile->Close(); }; void T2K_CC1pip_CH_XSec_1Dthpi_nu::FillEventVariables(FitEvent *event) { if (event->NumFSParticle(13) == 0 || event->NumFSParticle(211) == 0) return; TLorentzVector Pnu = event->GetNeutrinoIn()->fP; TLorentzVector Ppip = event->GetHMFSParticle(211)->fP; TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; double thpi = FitUtils::th(Pnu, Ppip); fXVar = thpi; return; }; //******************************************************************** bool T2K_CC1pip_CH_XSec_1Dthpi_nu::isSignal(FitEvent *event) { //******************************************************************** // This sample uses directional info on the pion so Michel e tag sample can not be included // i.e. we need reduce the pion variable phase space return SignalDef::isCC1pip_T2K_CH(event, EnuMin, EnuMax, false); } diff --git a/src/T2K/T2K_CC1pip_CH_XSec_1Dthq3pi_nu.cxx b/src/T2K/T2K_CC1pip_CH_XSec_1Dthq3pi_nu.cxx index e3a1cbb..1426e80 100644 --- a/src/T2K/T2K_CC1pip_CH_XSec_1Dthq3pi_nu.cxx +++ b/src/T2K/T2K_CC1pip_CH_XSec_1Dthq3pi_nu.cxx @@ -1,119 +1,115 @@ #include #include "T2K_SignalDef.h" #include "T2K_CC1pip_CH_XSec_1Dthq3pi_nu.h" // The constructor T2K_CC1pip_CH_XSec_1Dthq3pi_nu::T2K_CC1pip_CH_XSec_1Dthq3pi_nu(std::string inputfile, FitWeight *rw, std::string type, std::string fakeDataFile){ fName = "T2K_CC1pip_CH_XSec_1Dthq3pi_nu"; fPlotTitles = "; #theta_{q_{3},#pi} (radians); d#sigma/d#theta_{q_{3},#pi} (cm^{2}/(radian)/nucleon)"; EnuMin = 0.; EnuMax = 100.; fIsDiag = false; Measurement1D::SetupMeasurement(inputfile, type, rw, fakeDataFile); this->SetDataValues(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/ThetaQ3Pi.root"); this->SetCovarMatrix(GeneralUtils::GetTopLevelDir()+"/data/T2K/CC1pip/CH/ThetaQ3Pi.root"); - + SetShapeCovar(); this->SetupDefaultHist(); this->fScaleFactor = (GetEventHistogram()->Integral("width")*1E-38)/double(fNEvents)/TotalIntegratedFlux("width"); }; // Override this for now // Should really have Measurement1D do this properly though void T2K_CC1pip_CH_XSec_1Dthq3pi_nu::SetDataValues(std::string fileLocation) { LOG(DEB) << "Reading: " << this->fName << "\nData: " << fileLocation.c_str() << std::endl; TFile *dataFile = new TFile(fileLocation.c_str()); //truly great .root file! // Don't want the last bin of dataCopy TH1D *dataCopy = (TH1D*)(dataFile->Get("hResult_sliced_0_1"))->Clone(); double *binEdges = new double[dataCopy->GetNbinsX()-1]; LOG(DEB) << dataCopy->GetNbinsX() << std::endl; for (int i = 1; i < dataCopy->GetNbinsX()+1; i++) { binEdges[i-1] = dataCopy->GetBinLowEdge(i); LOG(DEB) << i-1 << " " << binEdges[i-1] << std::endl; } for (int i = 0; i < dataCopy->GetNbinsX()+5; i++) { LOG(DEB) << "binEdges[" << i << "] = " << binEdges[i] << std::endl; } fDataHist = new TH1D((fName+"_data").c_str(), (fName+"_data"+fPlotTitles).c_str(), dataCopy->GetNbinsX()-1, binEdges); for (int i = 0; i < fDataHist->GetNbinsX(); i++) { fDataHist->SetBinContent(i+1, dataCopy->GetBinContent(i+1)*1E-38); fDataHist->SetBinError(i+1, dataCopy->GetBinError(i+1)*1E-38); LOG(DEB) << fDataHist->GetBinLowEdge(i+1) << " " << fDataHist->GetBinContent(i+1) << " " << fDataHist->GetBinError(i+1) << std::endl; } fDataHist->SetDirectory(0); //should disassociate fDataHist with dataFile fDataHist->SetNameTitle((fName+"_data").c_str(), (fName+"_MC"+fPlotTitles).c_str()); dataFile->Close(); }; // Override this for now // Should really have Measurement1D do this properly though void T2K_CC1pip_CH_XSec_1Dthq3pi_nu::SetCovarMatrix(std::string fileLocation) { LOG(DEB) << "Covariance: " << fileLocation.c_str() << std::endl; TFile *dataFile = new TFile(fileLocation.c_str()); //truly great .root file! TH2D *covarMatrix = (TH2D*)(dataFile->Get("TMatrixDBase;1"))->Clone(); int nBinsX = covarMatrix->GetXaxis()->GetNbins(); int nBinsY = covarMatrix->GetYaxis()->GetNbins(); if ((nBinsX != nBinsY)) ERR(WRN) << "covariance matrix not square!" << std::endl; // First bin is underflow, last bin is overflow - this->covar = new TMatrixDSym(nBinsX-2); this->fFullCovar = new TMatrixDSym(nBinsX-2); // First two entries are BS // Last entry is BS for (int i = 1; i < nBinsX-1; i++) { for (int j = 1; j < nBinsY-1; j++) { - (*this->covar)(i-1, j-1) = covarMatrix->GetBinContent(i+1, j+1); //adds syst+stat covariances (*this->fFullCovar)(i-1, j-1) = covarMatrix->GetBinContent(i+1, j+1); //adds syst+stat covariances - LOG(DEB) << "covar(" << i-1 << ", " << j-1 << ") = " << (*this->covar)(i-1,j-1) << std::endl; + LOG(DEB) << "fFullCovar(" << i-1 << ", " << j-1 << ") = " << (*this->fFullCovar)(i-1,j-1) << std::endl; } } //should now have set covariance, I hope - TDecompChol tempMat = TDecompChol(*this->covar); - this->covar = new TMatrixDSym(nBinsX, tempMat.Invert().GetMatrixArray(), ""); - // *this->covar *= 1E-38*1E-38; - + this->fDecomp = StatUtils::GetDecomp(this->fFullCovar); + this->covar = StatUtils::GetInvert(this->fFullCovar); return; }; void T2K_CC1pip_CH_XSec_1Dthq3pi_nu::FillEventVariables(FitEvent *event) { if (event->NumFSParticle(13) == 0 || event->NumFSParticle(211) == 0) return; TLorentzVector Pnu = event->GetNeutrinoIn()->fP; TLorentzVector Ppip = event->GetHMFSParticle(211)->fP; TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; double th_q3_pi = FitUtils::thq3pi_CC1pip_T2K(Pnu, Pmu, Ppip); fXVar = th_q3_pi; return; }; //******************************************************************** bool T2K_CC1pip_CH_XSec_1Dthq3pi_nu::isSignal(FitEvent *event) { //******************************************************************** // This sample requires pion directional information so can not include Michel tag sample // i.e. will need to cut the pion phase space return SignalDef::isCC1pip_T2K_CH(event, EnuMin, EnuMax, false); }