diff --git a/examples/B3piKMNoAdler.dat b/examples/B3piKMNoAdler.dat deleted file mode 100644 index da79d74..0000000 --- a/examples/B3piKMNoAdler.dat +++ /dev/null @@ -1,45 +0,0 @@ -# -# Copyright 2016 University of Warwick -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# -# Laura++ package authors: -# John Back -# Paul Harrison -# Thomas Latham -# -# File for pi-pi K-matrix S-wave amplitude. -# Coeffs from [hep-ex]0804.2089, table I, page 10 (RHS). -# First, define the indices of the phase-space channels defined by -# the LauKMatrixPropagator::KMatrixChannel enum -Channels 1 2 3 4 5 -# Next, define the bare pole masses and coupling constants g^(alpha)_j. -# First int = alpha (pole index number), 2nd number = pole mass (GeV/c^2), while -# the other numbers are constants for each channel j = 1 to N for the given pole -Pole 1 0.65100 0.22889 -0.55377 0.0 -0.39899 -0.34639 -Pole 2 1.20360 0.94128 0.55095 0.0 0.39065 0.31503 -Pole 3 1.55817 0.36856 0.23888 0.55639 0.18340 0.18681 -Pole 4 1.21000 0.33650 0.40907 0.85679 0.19906 -0.00984 -Pole 5 1.82206 0.18171 -0.17558 -0.79658 -0.00355 0.22358 -# Next, define the scattering f_{ij} constants. Here, the first integer is -# the scattering channel index, while the other numbers = f_{ij}, j = 1 to N -Scatt 1 0.23399 0.15044 -0.20545 0.32825 0.35412 -# We usually assume symmetry: f_{ji} = f_{ij}. Otherwise, uncomment this line -#ScattSymmetry 0 -# Finally, define the "Adler-zero" constants -# mSq0 s0Scatt s0Prod sA sA0 -mSq0 1.0 -s0Scatt -3.92637 -s0Prod 0.0 -sA 0.0 -sA0 0.0 diff --git a/examples/B3piKMNoAdler.json b/examples/B3piKMNoAdler.json new file mode 100644 index 0000000..8efdd3b --- /dev/null +++ b/examples/B3piKMNoAdler.json @@ -0,0 +1,59 @@ +{ + "comment" : "File for pi-pi K-matrix S-wave amplitude. Coeffs from [hep-ex]0804.2089, table I, page 10 (RHS).", + "channels" : [ "PiPi", "KK", "FourPi", "EtaEta", "EtaEtaP" ], + "poles" : [ + { + "index" : 1, + "mass" : 0.65100, + "couplings" : [ 0.22889, -0.55377, 0.0, -0.39899, -0.34639 ] + }, + { + "index" : 2, + "mass" : 1.20360, + "couplings" : [ 0.94128, 0.55095, 0.0, 0.39065, 0.31503 ] + }, + { + "index" : 3, + "mass" : 1.55817, + "couplings" : [ 0.36856, 0.23888, 0.55639, 0.18340, 0.18681 ] + }, + { + "index" : 4, + "mass" : 1.21000, + "couplings" : [ 0.33650, 0.40907, 0.85679, 0.19906, -0.00984 ] + }, + { + "index" : 5, + "mass" : 1.82206, + "couplings" : [ 0.18171, -0.17558, -0.79658, -0.00355, 0.22358 ] + } + ], + "scattering" : [ + { + "index" : 1, + "couplings" : [ 0.23399, 0.15044, -0.20545, 0.32825, 0.35412 ] + } + ], + "parameters" : [ + { + "name" : "mSq0", + "value" : 1.0 + }, + { + "name" : "s0Scatt", + "value" : -3.92637 + }, + { + "name" : "s0Prod", + "value" : 0.0 + }, + { + "name" : "sA", + "value" : 0.0 + }, + { + "name" : "sA0", + "value" : 0.0 + } + ] +} diff --git a/examples/B3piKMPoles.dat b/examples/B3piKMPoles.dat deleted file mode 100644 index 03288e0..0000000 --- a/examples/B3piKMPoles.dat +++ /dev/null @@ -1,45 +0,0 @@ -# -# Copyright 2016 University of Warwick -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# -# Laura++ package authors: -# John Back -# Paul Harrison -# Thomas Latham -# -# File for pi-pi K-matrix S-wave amplitude. -# Coeffs from [hep-ex]0804.2089, table I, page 10 (RHS). -# First, define the indices of the phase-space channels defined by -# the LauKMatrixPropagator::KMatrixChannel enum -Channels 1 2 3 4 5 -# Next, define the bare pole masses and coupling constants g^(alpha)_j. -# First int = alpha (pole index number), 2nd number = pole mass (GeV/c^2), while -# the other numbers are constants for each channel j = 1 to N for the given pole -Pole 1 0.65100 0.22889 -0.55377 0.0 -0.39899 -0.34639 -Pole 2 1.20360 0.94128 0.55095 0.0 0.39065 0.31503 -Pole 3 1.55817 0.36856 0.23888 0.55639 0.18340 0.18681 -Pole 4 1.21000 0.33650 0.40907 0.85679 0.19906 -0.00984 -Pole 5 1.82206 0.18171 -0.17558 -0.79658 -0.00355 0.22358 -# Next, define the scattering f_{ij} constants. Here, the first integer is -# the scattering channel index, while the other numbers = f_{ij}, j = 1 to N -Scatt 1 0.0 0.0 0.0 0.0 0.0 -# We usually assume symmetry: f_{ji} = f_{ij}. Otherwise, uncomment this line -#ScattSymmetry 0 -# Finally, define the "Adler-zero" constants -# mSq0 s0Scatt s0Prod sA sA0 -mSq0 1.0 -s0Scatt -3.92637 -s0Prod 0.0 -sA 0.0 -sA0 0.0 diff --git a/examples/B3piKMPoles.json b/examples/B3piKMPoles.json new file mode 100644 index 0000000..7cc7104 --- /dev/null +++ b/examples/B3piKMPoles.json @@ -0,0 +1,59 @@ +{ + "comment" : "File for pi-pi K-matrix S-wave amplitude. Coeffs from [hep-ex]0804.2089, table I, page 10 (RHS).", + "channels" : [ "PiPi", "KK", "FourPi", "EtaEta", "EtaEtaP" ], + "poles" : [ + { + "index" : 1, + "mass" : 0.65100, + "couplings" : [ 0.22889, -0.55377, 0.0, -0.39899, -0.34639 ] + }, + { + "index" : 2, + "mass" : 1.20360, + "couplings" : [ 0.94128, 0.55095, 0.0, 0.39065, 0.31503 ] + }, + { + "index" : 3, + "mass" : 1.55817, + "couplings" : [ 0.36856, 0.23888, 0.55639, 0.18340, 0.18681 ] + }, + { + "index" : 4, + "mass" : 1.21000, + "couplings" : [ 0.33650, 0.40907, 0.85679, 0.19906, -0.00984 ] + }, + { + "index" : 5, + "mass" : 1.82206, + "couplings" : [ 0.18171, -0.17558, -0.79658, -0.00355, 0.22358 ] + } + ], + "scattering" : [ + { + "index" : 1, + "couplings" : [ 0.0, 0.0, 0.0, 0.0, 0.0 ] + } + ], + "parameters" : [ + { + "name" : "mSq0", + "value" : 1.0 + }, + { + "name" : "s0Scatt", + "value" : -3.92637 + }, + { + "name" : "s0Prod", + "value" : 0.0 + }, + { + "name" : "sA", + "value" : 0.0 + }, + { + "name" : "sA0", + "value" : 0.0 + } + ] +} diff --git a/examples/B3piKMSVPs.dat b/examples/B3piKMSVPs.dat deleted file mode 100644 index c624758..0000000 --- a/examples/B3piKMSVPs.dat +++ /dev/null @@ -1,45 +0,0 @@ -# -# Copyright 2016 University of Warwick -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# -# Laura++ package authors: -# John Back -# Paul Harrison -# Thomas Latham -# -# File for pi-pi K-matrix S-wave amplitude. -# Coeffs from [hep-ex]0804.2089, table I, page 10 (RHS). -# First, define the indices of the phase-space channels defined by -# the LauKMatrixPropagator::KMatrixChannel enum -Channels 1 2 3 4 5 -# Next, define the bare pole masses and coupling constants g^(alpha)_j. -# First int = alpha (pole index number), 2nd number = pole mass (GeV/c^2), while -# the other numbers are constants for each channel j = 1 to N for the given pole -#Pole 1 0.65100 0.0 0.0 0.0 0.0 0.0 -#Pole 2 1.20360 0.0 0.0 0.0 0.0 0.0 -#Pole 3 1.55817 0.0 0.0 0.0 0.0 0.0 -#Pole 4 1.21000 0.0 0.0 0.0 0.0 0.0 -#Pole 5 1.82206 0.0 0.0 0.0 0.0 0.0 -# Next, define the scattering f_{ij} constants. Here, the first integer is -# the scattering channel index, while the other numbers = f_{ij}, j = 1 to N -Scatt 1 0.23399 0.15044 -0.20545 0.32825 0.35412 -# We usually assume symmetry: f_{ji} = f_{ij}. Otherwise, uncomment this line -#ScattSymmetry 0 -# Finally, define the "Adler-zero" constants -# mSq0 s0Scatt s0Prod sA sA0 -mSq0 1.0 -s0Scatt -3.92637 -s0Prod 0.0 -sA 0.0 -sA0 0.0 diff --git a/examples/B3piKMSVPs.json b/examples/B3piKMSVPs.json new file mode 100644 index 0000000..05fd8f3 --- /dev/null +++ b/examples/B3piKMSVPs.json @@ -0,0 +1,34 @@ +{ + "comment" : "File for pi-pi K-matrix S-wave amplitude. Coeffs from [hep-ex]0804.2089, table I, page 10 (RHS).", + "channels" : [ "PiPi", "KK", "FourPi", "EtaEta", "EtaEtaP" ], + "poles" : [ + ], + "scattering" : [ + { + "index" : 1, + "couplings" : [ 0.23399, 0.15044, -0.20545, 0.32825, 0.35412 ] + } + ], + "parameters" : [ + { + "name" : "mSq0", + "value" : 1.0 + }, + { + "name" : "s0Scatt", + "value" : -3.92637 + }, + { + "name" : "s0Prod", + "value" : 0.0 + }, + { + "name" : "sA", + "value" : 0.0 + }, + { + "name" : "sA0", + "value" : 0.0 + } + ] +} diff --git a/examples/B3piKMatrixCoeff.dat b/examples/B3piKMatrixCoeff.dat deleted file mode 100644 index 0b27f82..0000000 --- a/examples/B3piKMatrixCoeff.dat +++ /dev/null @@ -1,45 +0,0 @@ -# -# Copyright 2016 University of Warwick -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# -# Laura++ package authors: -# John Back -# Paul Harrison -# Thomas Latham -# -# File for pi-pi K-matrix S-wave amplitude. -# Coeffs from [hep-ex]0804.2089, table I, page 10 (RHS). -# First, define the indices of the phase-space channels defined by -# the LauKMatrixPropagator::KMatrixChannel enum -Channels 1 2 3 4 5 -# Next, define the bare pole masses and coupling constants g^(alpha)_j. -# First int = alpha (pole index number), 2nd number = pole mass (GeV/c^2), while -# the other numbers are constants for each channel j = 1 to N for the given pole -Pole 1 0.65100 0.22889 -0.55377 0.0 -0.39899 -0.34639 -Pole 2 1.20360 0.94128 0.55095 0.0 0.39065 0.31503 -Pole 3 1.55817 0.36856 0.23888 0.55639 0.18340 0.18681 -Pole 4 1.21000 0.33650 0.40907 0.85679 0.19906 -0.00984 -Pole 5 1.82206 0.18171 -0.17558 -0.79658 -0.00355 0.22358 -# Next, define the scattering f_{ij} constants. Here, the first integer is -# the scattering channel index, while the other numbers = f_{ij}, j = 1 to N -Scatt 1 0.23399 0.15044 -0.20545 0.32825 0.35412 -# We usually assume symmetry: f_{ji} = f_{ij}. Otherwise, uncomment this line -#ScattSymmetry 0 -# Finally, define the "Adler-zero" constants -# mSq0 s0Scatt s0Prod sA sA0 -mSq0 1.0 -s0Scatt -3.92637 -s0Prod -0.07 -sA 1.0 -sA0 -0.15 diff --git a/examples/B3piKMatrixCoeff.json b/examples/B3piKMatrixCoeff.json new file mode 100644 index 0000000..51240e4 --- /dev/null +++ b/examples/B3piKMatrixCoeff.json @@ -0,0 +1,59 @@ +{ + "comment" : "File for pi-pi K-matrix S-wave amplitude. Coeffs from [hep-ex]0804.2089, table I, page 10 (RHS).", + "channels" : [ "PiPi", "KK", "FourPi", "EtaEta", "EtaEtaP" ], + "poles" : [ + { + "index" : 1, + "mass" : 0.65100, + "couplings" : [ 0.22889, -0.55377, 0.0, -0.39899, -0.34639 ] + }, + { + "index" : 2, + "mass" : 1.20360, + "couplings" : [ 0.94128, 0.55095, 0.0, 0.39065, 0.31503 ] + }, + { + "index" : 3, + "mass" : 1.55817, + "couplings" : [ 0.36856, 0.23888, 0.55639, 0.18340, 0.18681 ] + }, + { + "index" : 4, + "mass" : 1.21000, + "couplings" : [ 0.33650, 0.40907, 0.85679, 0.19906, -0.00984 ] + }, + { + "index" : 5, + "mass" : 1.82206, + "couplings" : [ 0.18171, -0.17558, -0.79658, -0.00355, 0.22358 ] + } + ], + "scattering" : [ + { + "index" : 1, + "couplings" : [ 0.23399, 0.15044, -0.20545, 0.32825, 0.35412 ] + } + ], + "parameters" : [ + { + "name" : "mSq0", + "value" : 1.0 + }, + { + "name" : "s0Scatt", + "value" : -3.92637 + }, + { + "name" : "s0Prod", + "value" : -0.07 + }, + { + "name" : "sA", + "value" : 1.0 + }, + { + "name" : "sA0", + "value" : -0.15 + } + ] +} diff --git a/examples/B3piKMatrixMassProj.cc b/examples/B3piKMatrixMassProj.cc.in similarity index 97% rename from examples/B3piKMatrixMassProj.cc rename to examples/B3piKMatrixMassProj.cc.in index 2e2e9ce..dcfa63a 100644 --- a/examples/B3piKMatrixMassProj.cc +++ b/examples/B3piKMatrixMassProj.cc.in @@ -1,433 +1,434 @@ /* Copyright 2016 University of Warwick Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0 Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License. */ /* Laura++ package authors: John Back Paul Harrison Thomas Latham */ // Plot the mass projections of K-matrix amplitude terms generated over the Dalitz plot so that // we can see what shapes the poles and SVP components have. Also show the rho for comparison #include using std::cout; using std::endl; #include "LauCartesianCPCoeffSet.hh" #include "LauDaughters.hh" #include "LauEffModel.hh" #include "LauIsobarDynamics.hh" #include "LauKMatrixPropagator.hh" #include "LauParameter.hh" #include "LauResonanceMaker.hh" #include "LauSimpleFitModel.hh" #include "LauVetoes.hh" #include "TGraph.h" #include "TCanvas.h" #include "TStyle.h" #include "TROOT.h" #include "TSystem.h" #include "TString.h" #include "TH1.h" #include "TH2.h" #include "TTree.h" #include "TFile.h" #include "THStack.h" #include "TLegend.h" #include #include void createHistos(const std::string& rootFileName); void plotHistos(const std::string& rootFileName); std::vector genMassHistos(Int_t index, Int_t NSignal); Bool_t useProdAdler = kFALSE; int main(const int argc, const char** argv) { int histFlag(0), prodFlag(0); if (argc > 1) { histFlag = std::atoi(argv[1]); } if (argc > 2) { prodFlag = std::atoi(argv[2]); if (prodFlag != 0) {useProdAdler = kTRUE;} } std::string rootFileName("KMatrixMassHistos.root"); if (histFlag) {createHistos(rootFileName);} plotHistos(rootFileName); } void plotHistos(const std::string& rootFileName) { cout<<"Running plotHistos for "<SetStyle("Plain"); gStyle->SetOptStat(0); theCanvas->Clear(); theCanvas->UseCurrentStyle(); TFile* rootFile = TFile::Open(rootFileName.c_str(), "read"); rootFile->cd(); // Plot rho with K-matrix poles TLegend* poleLegend = new TLegend(0.7, 0.5, 0.90, 0.90, ""); poleLegend->SetFillColor(0); poleLegend->SetTextSize(0.03); THStack* poleStack = new THStack(); TH1* rhoHist = dynamic_cast(rootFile->Get("Rho_m13")); rhoHist->SetLineColor(kBlack); rhoHist->SetLineWidth(2); rhoHist->Scale(1.0/rhoHist->Integral()); poleStack->Add(rhoHist); poleLegend->AddEntry(rhoHist, "#rho(770)", "l"); Int_t colours[5] = {kRed, kViolet, kBlue, kGreen+2, kOrange+2}; Int_t i(0), iN(5); for (i = 0; i < iN; i++) { Int_t j = i + 1; TString poleName("KMPole"); poleName += j; poleName += "_m13"; TH1* poleHist = dynamic_cast(rootFile->Get(poleName.Data())); poleHist->SetLineColor(colours[i]); poleHist->SetLineWidth(2); // Normalise histogram to unit area poleHist->Scale(1.0/poleHist->Integral()); poleStack->Add(poleHist); TString poleLabel("Pole"); poleLabel += j; poleLegend->AddEntry(poleHist, poleLabel, "l"); } poleStack->Draw("nostackc"); TAxis* poleXAxis = poleStack->GetXaxis(); poleXAxis->SetTitle("m(#pi^{+}#pi^{-}) (GeV/c^{2})"); poleXAxis->CenterTitle(kTRUE); TAxis* poleYAxis = poleStack->GetYaxis(); poleYAxis->SetTitle("Normalised intensity"); poleYAxis->CenterTitle(kTRUE); poleYAxis->SetTitleOffset(1.25); poleLegend->Draw(); theCanvas->Update(); theCanvas->Print("B3piKMatrixProjPoles.png"); theCanvas->Print("B3piKMatrixProjPoles.eps"); // Plot NR with K-matrix SVPs TLegend* SVPLegend = new TLegend(0.7, 0.5, 0.90, 0.90, ""); SVPLegend->SetFillColor(0); SVPLegend->SetTextSize(0.03); THStack* SVPStack = new THStack(); SVPStack->Add(rhoHist); SVPLegend->AddEntry(rhoHist, "rho(770)", "l"); for (i = 0; i < iN; i++) { Int_t j = i + 1; TString SVPName("KMSVP"); SVPName += j; SVPName += "_m13"; TH1* SVPHist = dynamic_cast(rootFile->Get(SVPName.Data())); SVPHist->SetLineColor(colours[i]); SVPHist->SetLineWidth(2); SVPHist->Scale(1.0/SVPHist->Integral()); SVPStack->Add(SVPHist); TString SVPLabel("SVP"); SVPLabel += j; SVPLegend->AddEntry(SVPHist, SVPLabel, "l"); } SVPStack->Draw("nostackc"); TAxis* SVPXAxis = SVPStack->GetXaxis(); SVPXAxis->SetTitle("m(#pi^{+}#pi^{-}) (GeV/c^{2})"); SVPXAxis->CenterTitle(kTRUE); TAxis* SVPYAxis = SVPStack->GetYaxis(); SVPYAxis->SetTitle("Normalised intensity"); SVPYAxis->CenterTitle(kTRUE); SVPYAxis->SetTitleOffset(1.25); SVPLegend->Draw(); theCanvas->Update(); theCanvas->Print("B3piKMatrixProjSVPs.png"); theCanvas->Print("B3piKMatrixProjSVPs.eps"); delete SVPStack; delete SVPLegend; delete poleStack; delete poleLegend; rootFile->Close(); } void createHistos(const std::string& rootFileName) { cout<<"rootFileName = "< > histMap; // Create mass/DP histograms for the different amplitude terms for (i = 0; i < iN; i++) { std::vector histVect = genMassHistos(i, NSignal); histMap[i] = histVect; } TFile* rootFile = TFile::Open(rootFileName.c_str(), "recreate"); rootFile->cd(); for (i = 0; i < iN; i++) { cout<<"Writing histograms for i = "< histVect = histMap[i]; std::vector::iterator iter; for (iter = histVect.begin(); iter != histVect.end(); ++iter) { TH1* theHist = *iter; theHist->Write(); } } rootFile->Close(); } std::vector genMassHistos(Int_t index, Int_t NSignal) { // Vector to store output mass/DP histograms std::vector histVect; // Define Dalitz plot: kinematics and resonances Bool_t squareDP = kTRUE; LauDaughters* theDaughters = new LauDaughters("B+", "pi+", "pi+", "pi-", squareDP); // Apply some vetoes to the DP LauVetoes* vetoes = new LauVetoes(); // Define the efficiency model LauEffModel* theEffModel = new LauEffModel(theDaughters, vetoes); // Create the isobar model LauResonanceMaker& resMaker = LauResonanceMaker::get(); resMaker.setDefaultBWRadius(LauBlattWeisskopfFactor::Category::Parent, 4.0); resMaker.setDefaultBWRadius(LauBlattWeisskopfFactor::Category::Light, 4.0); resMaker.fixBWRadius(LauBlattWeisskopfFactor::Category::Parent, kTRUE); resMaker.fixBWRadius(LauBlattWeisskopfFactor::Category::Light, kTRUE); // Define the "S-wave" K-matrix propagator - Int_t nChannels = 5; - Int_t nPoles = 5; - Int_t resPairInt = 1; - Int_t KMatrixIndex = 1; // for S-wave + const UInt_t nChannels = 5; + const UInt_t nPoles = 5; + const UInt_t resPairInt = 1; + const UInt_t KMatrixIndex = 1; // for S-wave // Add one pole or SVP amplitude for the K-matrix and plot its mass projections // so that we can see the amplitude/intensity shape of the different terms Double_t aSqMaxValue = 0.3; Double_t nSigEvents = NSignal*1.0; Bool_t fixNSigEvents = kTRUE; Int_t nExpt = 1; Int_t firstExpt = 0; LauIsobarDynamics* theSigModel = new LauIsobarDynamics(theDaughters, theEffModel); theSigModel->setASqMaxValue(aSqMaxValue); TString theName("KMPole"); if (index >= 5) {theName = "KMSVP";} TString parName(""); if (index == 11) { // Simple Non-resonant cout<<"Non-resonant"<addResonance("NonReson", 0, LauAbsResonance::ResonanceModel::FlatNR); } else if (index == 10) { // Rho resonance cout<<"Using rho(770)"<addResonance("rho0(770)", 1, LauAbsResonance::ResonanceModel::GS); } else { // K-matrix amplitude - theSigModel->defineKMatrixPropagator("KMSWave", "B3piKMatrixCoeff.dat", resPairInt, + const TString jsonFileName { "@CMAKE_INSTALL_PREFIX@/@CMAKE_INSTALL_DATADIR@/@CMAKE_PROJECT_NAME@/B3piKMatrixCoeff.json" }; + theSigModel->defineKMatrixPropagator("KMSWave", jsonFileName, resPairInt, nChannels, nPoles, KMatrixIndex); cout<<"Implementing K matrix term number "<addKMatrixProdSVP(theName.Data(), "KMSWave", j, useProdAdler); } } // Create the fit model LauSimpleFitModel* fitModel = new LauSimpleFitModel(theSigModel); std::vector> coeffset; Bool_t fixPar = kTRUE; Bool_t fixCP = kTRUE; Bool_t doTwoStageFit = kFALSE; coeffset.push_back(std::make_unique(parName.Data(), 1.0, 0.0, 0.0, 0.0, fixPar, fixPar, fixCP, fixCP, doTwoStageFit, doTwoStageFit)); for ( auto& coeff : coeffset ) { fitModel->setAmpCoeffSet( std::move(coeff) ); } // Set the number of signal events and the number of experiments LauParameter* signalEvents = new LauParameter("signalEvents", nSigEvents, -2.0*nSigEvents, 2.0*nSigEvents, fixNSigEvents); fitModel->setNSigEvents(signalEvents); fitModel->setNExpts(nExpt, firstExpt); // Execute the generation/fit TString treeName(""); TString dummyName(""); TString tableFileName(""); TString dataFileName("Gen"); dataFileName += theName; dataFileName += ".root"; fitModel->run("gen", dataFileName, treeName, dummyName, tableFileName); // Create DP and mass projection plots so that we can see the K matrix // amplitude term shapes // Histogram of the invariant mass for pi+ pi- pairs: m13 & m23 TString m13Name(theName); m13Name += "_m13"; TH1D* m13Hist = new TH1D(m13Name.Data(), "", 100, 0.0, 5.5); m13Hist->SetXTitle("m(#pi^{+}#pi^{-})"); m13Hist->SetDirectory(0); histVect.push_back(m13Hist); // Histogram of the wrong sign invariant mass pi+ pi+: m12 TString m12Name(theName); m12Name += "_m12"; TH1D* m12Hist = new TH1D(m12Name.Data(), "", 100, 0.0, 5.5); m12Hist->SetXTitle("m(#pi^{+}#pi^{+})"); m12Hist->SetDirectory(0); histVect.push_back(m12Hist); // Dalitz plot 2d histogram TString DPName(theName); DPName += "_DP"; TH2D* DPHist = new TH2D(DPName.Data(), "", 100, 0.0, 14.0, 100, 0.0, 28.0); //TH2D* DPHist = new TH2D("DPHist", "", 100, 0.0, 28.0, 100, 0.0, 28.0); DPHist->SetXTitle("m^{2}(#pi^{+}#pi^{-})_{min}"); DPHist->SetYTitle("m^{2}(#pi^{+}#pi^{-})_{max}"); DPHist->SetDirectory(0); histVect.push_back(DPHist); // Simulation results file TFile* genFile = TFile::Open(dataFileName.Data(), "read"); genFile->cd(); TTree* genTree = dynamic_cast(genFile->Get("genResults")); genTree->SetBranchStatus("*", 0); genTree->SetBranchStatus("m12", 1); genTree->SetBranchStatus("m23", 1); genTree->SetBranchStatus("m13", 1); Double_t m12, m23, m13; genTree->SetBranchAddress("m12", &m12); genTree->SetBranchAddress("m23", &m23); genTree->SetBranchAddress("m13", &m13); Int_t nEntries = genTree->GetEntries(); cout<<"nEntries = "<GetEntry(k); Double_t mMin = m13; Double_t mMax = m23; if (m23 < m13) { mMin = m23; mMax = m13; } m13Hist->Fill(mMin); m13Hist->Fill(mMax); m12Hist->Fill(m12); DPHist->Fill(mMin*mMin, mMax*mMax); } // Cleanup genFile->Close(); delete signalEvents; delete fitModel; delete theSigModel; delete theEffModel; delete vetoes; delete theDaughters; // Return histograms return histVect; } diff --git a/examples/B3piKMatrixPlots.cc b/examples/B3piKMatrixPlots.cc.in similarity index 97% rename from examples/B3piKMatrixPlots.cc rename to examples/B3piKMatrixPlots.cc.in index cec74a3..f4d2ca2 100644 --- a/examples/B3piKMatrixPlots.cc +++ b/examples/B3piKMatrixPlots.cc.in @@ -1,1466 +1,1466 @@ /* Copyright 2016 University of Warwick Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0 Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License. */ /* Laura++ package authors: John Back Paul Harrison Thomas Latham */ // Code to plot amplitudes for the pi pi K matrix used for B -> 3pi: // K-matrix itself, propagator terms and production amplitudes #include using std::cout; using std::endl; #include "LauKMatrixPropagator.hh" #include "LauComplex.hh" #include "LauDaughters.hh" #include "LauKMatrixProdPole.hh" #include "LauKMatrixProdSVP.hh" #include "TGraph.h" #include "TCanvas.h" #include "TStyle.h" #include "TROOT.h" #include "TSystem.h" #include "TAxis.h" #include "TLine.h" #include "TMath.h" #include "TMatrixD.h" #include "TMultiGraph.h" #include "TString.h" #include "TList.h" #include "TPad.h" #include "TLegend.h" #include "TText.h" #include #include #include Double_t radToDeg = 180.0/TMath::Pi(); Bool_t doEPS = kFALSE; Bool_t doProdAdler = kFALSE; void doKScattPlots() { // Create plots of K_ij (with i >= j, since its symmetric), showing // only the pole and SVP contributions separately, along with the full // K matrix without and with the Adler zero suppression term for low s. // 1st page: K_1j (5), 2nd page: K_2j (4), 3rd page: K_3j, K_4j and K_55 (6) Int_t resPairInt(1); Int_t nChannels(5); Int_t KMatrixIndex(1); // pi pi S-wave Int_t nPoles(5); // Mass range for plots Double_t mMin(0.28); Double_t mMax(2.0); Int_t nSteps = 2000; Double_t dm = (mMax - mMin)/(nSteps*1.0); // Propagator with only the pole terms - TString KFile1("B3piKMPoles.dat"); + const TString KFile1 { "@CMAKE_INSTALL_PREFIX@/@CMAKE_INSTALL_DATADIR@/@CMAKE_PROJECT_NAME@/B3piKMPoles.json" }; LauKMatrixPropagator* propPoles = new LauKMatrixPropagator("propPoles", KFile1, resPairInt, nChannels, nPoles, KMatrixIndex); // Propagator with only the SVP terms - TString KFile2("B3piKMSVPs.dat"); + const TString KFile2 { "@CMAKE_INSTALL_PREFIX@/@CMAKE_INSTALL_DATADIR@/@CMAKE_PROJECT_NAME@/B3piKMSVPs.json" }; LauKMatrixPropagator* propSVPs = new LauKMatrixPropagator("propSVPs", KFile2, resPairInt, nChannels, 0, KMatrixIndex); // Propagator with poles and SVP terms, but no Adler zero scaling factor - TString KFile3("B3piKMNoAdler.dat"); + const TString KFile3 { "@CMAKE_INSTALL_PREFIX@/@CMAKE_INSTALL_DATADIR@/@CMAKE_PROJECT_NAME@/B3piKMNoAdler.json" }; LauKMatrixPropagator* propNoAdler = new LauKMatrixPropagator("propNoAdler", KFile3, resPairInt, nChannels, nPoles, KMatrixIndex); // Full K-matrix parameterisation with Adler zero - TString KFile4("B3piKMatrixCoeff.dat"); + const TString KFile4 { "@CMAKE_INSTALL_PREFIX@/@CMAKE_INSTALL_DATADIR@/@CMAKE_PROJECT_NAME@/B3piKMatrixCoeff.json" }; LauKMatrixPropagator* propFull = new LauKMatrixPropagator("propFull", KFile4, resPairInt, nChannels, nPoles, KMatrixIndex); // For each K_ij, there are 4 graphs. Therefore, need 15 K_ij multigraphs std::vector multiGraphVect; // Legend for the graphs (same for all multigraphs) TLegend* theLegend = new TLegend(0.2, 0.2, 0.8, 0.7); theLegend->SetHeader("Symmetric K matrix components"); Int_t iColours[4] = {kRed, kBlue, kMagenta, kBlack}; Double_t yLimit(5.0); Bool_t initLegend(kFALSE); Int_t iChannel(0), jChannel(0); for (iChannel = 0; iChannel < nChannels; iChannel++) { for (jChannel = iChannel; jChannel < nChannels; jChannel++) { TString label("Graphs_"); label += iChannel; label += "_"; label += jChannel; TString multiName("KMatrix"); multiName += label; TMultiGraph* multiGraph = new TMultiGraph(); multiGraph->SetName(multiName.Data()); // To limit the graphs abscissas getting to large (pole "singularities") multiGraph->SetMaximum(yLimit); multiGraph->SetMinimum(-yLimit); // Add to the storage vector multiGraphVect.push_back(multiGraph); // Graph with poles only TGraph* poleGraph = new TGraph(nSteps); TString poleName("Pole"); poleName += label; poleGraph->SetName(poleName.Data()); poleGraph->SetMarkerColor(iColours[0]); poleGraph->SetLineColor(iColours[0]); poleGraph->SetLineWidth(1); multiGraph->Add(poleGraph, "lp"); // Graph with SVPs only TGraph* SVPGraph = new TGraph(nSteps); TString SVPName("SVP"); SVPName += label; SVPGraph->SetName(SVPName.Data()); SVPGraph->SetMarkerColor(iColours[1]); SVPGraph->SetLineColor(iColours[1]); SVPGraph->SetLineWidth(1); multiGraph->Add(SVPGraph, "lp"); // Graph with poles and SVP, but no Adler zero factor TGraph* noAdlerGraph = new TGraph(nSteps); TString noAdlerName("NoAdler"); noAdlerName += label; noAdlerGraph->SetName(noAdlerName.Data()); noAdlerGraph->SetMarkerColor(iColours[2]); noAdlerGraph->SetLineColor(iColours[2]); noAdlerGraph->SetLineWidth(1); multiGraph->Add(noAdlerGraph, "lp"); // Full K matrix graphs TGraph* fullGraph = new TGraph(nSteps); TString fullName("Full"); fullName += label; fullGraph->SetName(fullName.Data()); fullGraph->SetMarkerColor(iColours[3]); fullGraph->SetLineColor(iColours[3]); fullGraph->SetLineWidth(1); multiGraph->Add(fullGraph, "lp"); // Initialise the graph legend if (!initLegend) { theLegend->AddEntry(poleGraph, "Poles only", "l"); theLegend->AddEntry(SVPGraph, "SVPs only (f_{1j} #neq 0)", "l"); theLegend->AddEntry(noAdlerGraph, "Poles & SVPs, no Adler zero", "l"); theLegend->AddEntry(fullGraph, "Poles & SVPs with Adler zero", "l"); initLegend = kTRUE; } } } for (Int_t iM = 0; iM < nSteps; iM++) { Double_t m = dm*iM + mMin; Double_t s = m*m; // Find the K matrices for each pole/SVP/Adler condition propPoles->updatePropagator(s); propSVPs->updatePropagator(s); propNoAdler->updatePropagator(s); propFull->updatePropagator(s); TMatrixD KMPoles = propPoles->getKMatrix(); TMatrixD KMSVPs = propSVPs->getKMatrix(); TMatrixD KMNoAdler = propNoAdler->getKMatrix(); TMatrixD KMFull = propFull->getKMatrix(); Int_t vectIndex(-1); for (iChannel = 0; iChannel < nChannels; iChannel++) { for (jChannel = iChannel; jChannel < nChannels; jChannel++) { vectIndex++; TMultiGraph* theGraphs = multiGraphVect[vectIndex]; TList* graphList = theGraphs->GetListOfGraphs(); TGraph* poleGraph = dynamic_cast(graphList->At(0)); poleGraph->SetPoint(iM, m, KMPoles[iChannel][jChannel]); TGraph* SVPGraph = dynamic_cast(graphList->At(1)); SVPGraph->SetPoint(iM, m, KMSVPs[iChannel][jChannel]); TGraph* noAdlerGraph = dynamic_cast(graphList->At(2)); noAdlerGraph->SetPoint(iM, m, KMNoAdler[iChannel][jChannel]); TGraph* fullGraph = dynamic_cast(graphList->At(3)); fullGraph->SetPoint(iM, m, KMFull[iChannel][jChannel]); } } } // Create plots TCanvas* theCanvas = new TCanvas("theCanvas", "", 1400, 1000); gROOT->SetStyle("Plain"); gStyle->SetOptStat(0); theCanvas->Clear(); theCanvas->UseCurrentStyle(); // K_1j graphs, j = 1 to 5 theCanvas->Divide(2,3); for (Int_t iL = 0; iL < 5; iL++) { theCanvas->cd(iL+1); TMultiGraph* graphs = multiGraphVect[iL]; graphs->Draw("a"); // Update the axes, which can only be done after drawing the graphs TAxis* xAxis = graphs->GetXaxis(); xAxis->SetTitle("m(#pi^{+}#pi^{-}) (GeV/c^{2})"); xAxis->SetLabelSize(0.05); xAxis->SetTitleSize(0.05); xAxis->CenterTitle(kTRUE); xAxis->SetTitleOffset(0.9); TString yLabel("K(1,"); yLabel += (iL+1); yLabel += ")"; TAxis* yAxis = graphs->GetYaxis(); yAxis->SetTitleOffset(0.9); yAxis->SetTitle(yLabel.Data()); yAxis->SetLabelSize(0.05); yAxis->SetTitleSize(0.05); yAxis->CenterTitle(kTRUE); // Update the canvas with the new axes settings theCanvas->Update(); } // Print the graph legend theCanvas->cd(6); theLegend->Draw(); theCanvas->Print("K_1jPlots.png"); if (doEPS) {theCanvas->Print("K_1jPlots.eps");} // K_2j graphs, j = 2 to 5 theCanvas->Clear(); theCanvas->Divide(2,2); for (Int_t iL = 0; iL < 4; iL++) { theCanvas->cd(iL+1); TMultiGraph* graphs = multiGraphVect[iL+5]; graphs->Draw("a"); // Update the axes, which can only be done after drawing the graphs TAxis* xAxis = graphs->GetXaxis(); xAxis->SetTitle("m(#pi^{+}#pi^{-}) (GeV/c^{2})"); xAxis->SetLabelSize(0.05); xAxis->SetTitleSize(0.05); xAxis->CenterTitle(kTRUE); xAxis->SetTitleOffset(0.9); TString yLabel("K(2,"); yLabel += (iL+2); yLabel += ")"; TAxis* yAxis = graphs->GetYaxis(); yAxis->SetTitleOffset(0.9); yAxis->SetTitle(yLabel.Data()); yAxis->SetLabelSize(0.05); yAxis->SetTitleSize(0.05); yAxis->CenterTitle(kTRUE); // Update the canvas with the new axes settings theCanvas->Update(); } theCanvas->Print("K_2jPlots.png"); if (doEPS) {theCanvas->Print("K_2jPlots.eps");} // K_3j (j = 3,4,5), K_4j (j = 4,5) and K_55 theCanvas->Clear(); theCanvas->Divide(2,3); for (Int_t iL = 0; iL < 6; iL++) { theCanvas->cd(iL+1); TMultiGraph* graphs = multiGraphVect[iL+9]; graphs->Draw("a"); // Update the axes, which can only be done after drawing the graphs TAxis* xAxis = graphs->GetXaxis(); xAxis->SetTitle("m(#pi^{+}#pi^{-}) (GeV/c^{2})"); xAxis->SetLabelSize(0.05); xAxis->SetTitleSize(0.05); xAxis->CenterTitle(kTRUE); xAxis->SetTitleOffset(0.9); TString yLabel("K("); if (iL == 0) { yLabel += "3,3)"; } else if (iL == 1) { yLabel += "3,4)"; } else if (iL == 2) { yLabel += "3,5)"; } else if (iL == 3) { yLabel += "4,4)"; } else if (iL == 4) { yLabel += "4,5)"; } else { yLabel += "5,5)"; } TAxis* yAxis = graphs->GetYaxis(); yAxis->SetTitleOffset(0.9); yAxis->SetTitle(yLabel.Data()); yAxis->SetLabelSize(0.05); yAxis->SetTitleSize(0.05); yAxis->CenterTitle(kTRUE); // Update the canvas with the new axes settings theCanvas->Update(); } theCanvas->Print("K_345jPlots.png"); if (doEPS) {theCanvas->Print("K_345jPlots.eps");} // Cleanup delete propPoles; delete propSVPs; delete propNoAdler; delete propFull; std::vector::iterator multiIter; for (multiIter = multiGraphVect.begin(); multiIter != multiGraphVect.end(); ++multiIter) { TMultiGraph* multiGraph = *multiIter; delete multiGraph; } delete theLegend; delete theCanvas; } void doPropagatorPlots() { // Create plots of the s variation of the elements of the // propagator (I - i K rho)^-1 Int_t resPairInt(1); Int_t nChannels(5); Int_t KMatrixIndex(1); // pi pi S-wave Int_t nPoles(5); // Mass range for plots Double_t mMin(0.28); Double_t mMax(2.0); Int_t nSteps = 2000; Double_t dm = (mMax - mMin)/(nSteps*1.0); // Find the point numbers for m = mMin, 0.5, 0.75, 1, 1.25, 1.5, 1.75 and 2 GeV std::vector points; std::vector colours, styles; std::vector massPoints; massPoints.push_back(mMin); massPoints.push_back(0.5); massPoints.push_back(0.75); massPoints.push_back(1.0); massPoints.push_back(1.25); massPoints.push_back(1.5); massPoints.push_back(1.75); massPoints.push_back(1.99999); Int_t NMain = static_cast(massPoints.size()); for (Int_t iP = 0; iP < NMain; iP++) { Int_t pointInt = static_cast((massPoints[iP] - mMin)/dm); points.push_back(pointInt); // Red or Blue Int_t iColour = kRed; if (iP > 3) {iColour = kBlue;} Int_t iStyle = kFullCircle; Int_t iPMod = iP%4; if (iPMod == 1) { iStyle = kOpenCircle; } else if (iPMod == 2) { iStyle = kFullSquare; } else if (iPMod == 3) { iStyle = kOpenSquare; } colours.push_back(iColour); styles.push_back(iStyle); } // Full K-matrix parameterisation - TString KFile("B3piKMatrixCoeff.dat"); + const TString KFile { "@CMAKE_INSTALL_PREFIX@/@CMAKE_INSTALL_DATADIR@/@CMAKE_PROJECT_NAME@/B3piKMatrixCoeff.json" }; LauKMatrixPropagator* propFull = new LauKMatrixPropagator("propFull", KFile, resPairInt, nChannels, nPoles, KMatrixIndex); // Maps of // Argand diagram graphs (full plot, then others showing only selected points) std::map > argMap; // Real, imaginary and magnitude amplitude multigraphs std::map > ampMap; // TLegend for argand graphs TLegend* argLegend = new TLegend(0.1, 0.1, 0.7, 0.995); argLegend->SetHeader("m(#pi^{+}#pi^{-}) points (GeV/c^{2})"); argLegend->SetTextSize(0.075); Bool_t initArgLegend(kFALSE); // Legend for the amplitude graphs (same for all multigraphs) TLegend* ampLegend = new TLegend(0.2, 0.2, 0.8, 0.7); ampLegend->SetHeader("Propagator (I - iK#rho)^{-1} components"); Bool_t initAmpLegend(kFALSE); Int_t iChannel(0), jChannel(0); for (iChannel = 0; iChannel < nChannels; iChannel++) { std::vector argGraphVect; std::vector ampGraphVect; for (jChannel = 0; jChannel < nChannels; jChannel++) { // Common numbering scheme TString numbers(""); numbers += iChannel; numbers += "_"; numbers += jChannel; // Argand plot showing phase motion TMultiGraph* argGraphs = new TMultiGraph(); TString argMultiLabel("PropMultiArg_"); argMultiLabel += numbers; argGraphs->SetName(argMultiLabel.Data()); argGraphVect.push_back(argGraphs); TGraph* argGraph = new TGraph(nSteps); TString argLabel("PropArg_"); argLabel += numbers; argGraph->SetName(argLabel.Data()); argGraph->SetLineWidth(1); argGraphs->Add(argGraph, "lp"); // Argand graphs showing only one highlighted point for (Int_t iG = 0; iG < NMain; iG++) { TGraph* pointGraph = new TGraph(1); TString pointLabel("PropPoint_"); pointLabel += iG; pointLabel += "_"; pointLabel += numbers; pointGraph->SetName(pointLabel.Data()); pointGraph->SetMarkerColor(colours[iG]); pointGraph->SetMarkerSize(1.1); pointGraph->SetMarkerStyle(styles[iG]); argGraphs->Add(pointGraph, "p"); } if (!initArgLegend) { TList* argList = argGraphs->GetListOfGraphs(); for (Int_t iG = 0; iG < NMain; iG++) { TGraph* pointGraph = dynamic_cast(argList->At(iG+1)); argLegend->AddEntry(pointGraph, Form("%.2f",massPoints[iG]), "p"); initArgLegend = kTRUE; } } // Real, imaginary and magnitude amplitude TMultiGraph* ampGraphs = new TMultiGraph(); TString multiLabel("PropMultiAmp_"); multiLabel += numbers; ampGraphs->SetName(multiLabel.Data()); ampGraphVect.push_back(ampGraphs); TGraph* realGraph = new TGraph(nSteps); TString realLabel("PropReal_"); realLabel += numbers; realGraph->SetName(realLabel.Data()); realGraph->SetLineColor(kRed); realGraph->SetMarkerColor(kRed); realGraph->SetLineWidth(1); ampGraphs->Add(realGraph, "lp"); TGraph* imagGraph = new TGraph(nSteps); TString imagLabel("PropImag_"); imagLabel += numbers; imagGraph->SetName(imagLabel.Data()); imagGraph->SetLineColor(kBlue); imagGraph->SetMarkerColor(kBlue); imagGraph->SetLineWidth(1); ampGraphs->Add(imagGraph, "lp"); TGraph* magGraph = new TGraph(nSteps); TString magLabel("PropMag_"); magLabel += numbers; magGraph->SetName(magLabel.Data()); magGraph->SetLineColor(kBlack); magGraph->SetMarkerColor(kBlack); magGraph->SetLineWidth(1); ampGraphs->Add(magGraph, "lp"); if (!initAmpLegend) { ampLegend->AddEntry(realGraph, "Real part", "l"); ampLegend->AddEntry(imagGraph, "Imag part", "l"); ampLegend->AddEntry(magGraph, "Magnitude", "l"); initAmpLegend = kTRUE; } } argMap[iChannel] = argGraphVect; ampMap[iChannel] = ampGraphVect; } // Fill the graphs for (Int_t iM = 0; iM < nSteps; iM++) { Double_t m = dm*iM + mMin; Double_t s = m*m; propFull->updatePropagator(s); TMatrixD realPropMatrix = propFull->getRealPropMatrix(); TMatrixD negImagPropMatrix = propFull->getNegImagPropMatrix(); for (iChannel = 0; iChannel < nChannels; iChannel++) { std::vector argGraphVect = argMap[iChannel]; std::vector ampGraphVect = ampMap[iChannel]; for (jChannel = 0; jChannel < nChannels; jChannel++) { // Get the real and imaginary parts of the propagator term Double_t realProp = realPropMatrix[iChannel][jChannel]; Double_t imagProp = -negImagPropMatrix[iChannel][jChannel]; Double_t magProp = sqrt(realProp*realProp + imagProp*imagProp); // Argand diagram: full plot and only selected points TMultiGraph* argGraphs = argGraphVect[jChannel]; TList* argList = argGraphs->GetListOfGraphs(); TGraph* argGraph = dynamic_cast(argList->At(0)); argGraph->SetPoint(iM, realProp, imagProp); // Graphs showing specific points for (Int_t iP = 0; iP < NMain; iP++) { if (points[iP] == iM) { TGraph* pointGraph = dynamic_cast(argList->At(iP+1)); pointGraph->SetPoint(0, realProp, imagProp); } } // Amplitude multigraphs TMultiGraph* ampGraphs = ampGraphVect[jChannel]; TList* ampList = ampGraphs->GetListOfGraphs(); TGraph* realGraph = dynamic_cast(ampList->At(0)); realGraph->SetPoint(iM, m, realProp); TGraph* imagGraph = dynamic_cast(ampList->At(1)); imagGraph->SetPoint(iM, m, imagProp); TGraph* magGraph = dynamic_cast(ampList->At(2)); magGraph->SetPoint(iM, m, magProp); } } } // Create plots TCanvas* theCanvas = new TCanvas("theCanvas", "", 1000, 1000); gROOT->SetStyle("Plain"); gStyle->SetOptStat(0); for (iChannel = 0; iChannel < nChannels; iChannel++) { std::vector argGraphVect = argMap[iChannel]; std::vector ampGraphVect = ampMap[iChannel]; theCanvas->Clear(); theCanvas->UseCurrentStyle(); theCanvas->Divide(2,3); // Argand plots for (jChannel = 0; jChannel < nChannels; jChannel++) { theCanvas->cd(jChannel+1); TMultiGraph* argGraphs = argGraphVect[jChannel]; argGraphs->Draw("a"); TString title("Propagator("); title += iChannel+1; title += ","; title += jChannel+1; title += ")"; argGraphs->SetTitle(title.Data()); TAxis* xAxis = argGraphs->GetXaxis(); xAxis->SetTitle("Real"); xAxis->SetLabelSize(0.05); xAxis->SetTitleSize(0.05); xAxis->CenterTitle(kTRUE); xAxis->SetTitleOffset(0.9); TAxis* yAxis = argGraphs->GetYaxis(); yAxis->SetTitle("Imag"); yAxis->SetLabelSize(0.05); yAxis->SetTitleSize(0.05); yAxis->CenterTitle(kTRUE); theCanvas->Update(); } theCanvas->cd(6); argLegend->Draw(); TString argPlotName("PropArg_"); argPlotName += iChannel; argPlotName += ".png"; theCanvas->Print(argPlotName.Data()); if (doEPS) { TString argEpsPlotName(argPlotName); argEpsPlotName.ReplaceAll(".png", ".eps"); theCanvas->Print(argEpsPlotName.Data()); } // Amplitude plots theCanvas->Clear(); theCanvas->UseCurrentStyle(); theCanvas->Divide(2,3); for (jChannel = 0; jChannel < nChannels; jChannel++) { theCanvas->cd(jChannel+1); TMultiGraph* ampGraphs = ampGraphVect[jChannel]; ampGraphs->Draw("a"); TAxis* xAxis = ampGraphs->GetXaxis(); xAxis->SetTitle("m(#pi^{+}#pi^{-}) (GeV/c^{2})"); xAxis->SetLabelSize(0.05); xAxis->SetTitleSize(0.05); xAxis->CenterTitle(kTRUE); xAxis->SetTitleOffset(0.9); TAxis* yAxis = ampGraphs->GetYaxis(); TString yLabel("Propagator("); yLabel += iChannel+1; yLabel += ","; yLabel += jChannel+1; yLabel += ")"; yAxis->SetTitleOffset(0.9); yAxis->SetTitle(yLabel.Data()); yAxis->SetLabelSize(0.05); yAxis->SetTitleSize(0.05); yAxis->CenterTitle(kTRUE); theCanvas->Update(); } theCanvas->cd(6); ampLegend->Draw(); TString ampPlotName("PropAmp_"); ampPlotName += iChannel; ampPlotName += ".png"; theCanvas->Print(ampPlotName.Data()); if (doEPS) { TString ampEpsPlotName(ampPlotName); ampEpsPlotName.ReplaceAll(".png", ".eps"); theCanvas->Print(ampEpsPlotName.Data()); } } // Cleanup delete propFull; std::map >::iterator argIter; for (argIter = argMap.begin(); argIter != argMap.end(); ++argIter) { std::vector argGraphs = argIter->second; std::vector::iterator theIter; for (theIter = argGraphs.begin(); theIter != argGraphs.end(); ++theIter) { TMultiGraph* aGraph = *theIter; delete aGraph; } argGraphs.clear(); } std::map >::iterator ampIter; for (ampIter = ampMap.begin(); ampIter != ampMap.end(); ++ampIter) { std::vector ampGraphs = ampIter->second; std::vector::iterator theIter; for (theIter = ampGraphs.begin(); theIter != ampGraphs.end(); ++theIter) { TMultiGraph* aGraph = *theIter; delete aGraph; } ampGraphs.clear(); } delete ampLegend; delete theCanvas; } void doProdPolePlots(Bool_t useProdAdler) { // Create plots of the pole production amplitudes using the first row of // the propagator (for pi pi S-wave). For each pole we sum over all channels Int_t resPairInt(1); Int_t nChannels(5); Int_t KMatrixIndex(1); // pi pi S-wave Int_t nPoles(5); // Mass range for plots Double_t mMin(0.28); Double_t mMax(2.0); Int_t nSteps = 2000; Double_t dm = (mMax - mMin)/(nSteps*1.0); // Find the point numbers for m = mMin, 0.5, 0.75, 1, 1.25, 1.5, 1.75 and 2 GeV std::vector points; std::vector colours, styles; std::vector massPoints; massPoints.push_back(mMin); massPoints.push_back(0.5); massPoints.push_back(0.75); massPoints.push_back(1.0); massPoints.push_back(1.25); massPoints.push_back(1.5); massPoints.push_back(1.75); massPoints.push_back(1.99999); Int_t NMain = static_cast(massPoints.size()); for (Int_t iP = 0; iP < NMain; iP++) { Int_t pointInt = static_cast((massPoints[iP] - mMin)/dm); points.push_back(pointInt); // Red or Blue Int_t iColour = kRed; if (iP > 3) {iColour = kBlue;} Int_t iStyle = kFullCircle; Int_t iPMod = iP%4; if (iPMod == 1) { iStyle = kOpenCircle; } else if (iPMod == 2) { iStyle = kFullSquare; } else if (iPMod == 3) { iStyle = kOpenSquare; } colours.push_back(iColour); styles.push_back(iStyle); } Double_t mPi = 0.13957018; Double_t mPiSq = mPi*mPi; Double_t mB = 5.279; Double_t mBSq = mB*mB; Double_t mSqTerm = 2.0*((mBSq/3.0) + mPiSq); // Full K-matrix parameterisation - TString KFile("B3piKMatrixCoeff.dat"); + const TString KFile { "@CMAKE_INSTALL_PREFIX@/@CMAKE_INSTALL_DATADIR@/@CMAKE_PROJECT_NAME@/B3piKMatrixCoeff.json" }; LauKMatrixPropagator* propFull = new LauKMatrixPropagator("propFull", KFile, resPairInt, nChannels, nPoles, KMatrixIndex); // Define kinematics pointer Bool_t squareDP = kTRUE; LauDaughters* theDaughters = new LauDaughters("B+", "pi+", "pi+", "pi-", squareDP); // Kinematics LauKinematics* kinematics = theDaughters->getKinematics(); std::vector poleVect; Int_t iPole(0); for (iPole = 0; iPole < nPoles; iPole++) { TString poleName("Pole"); poleName += iPole; LauKMatrixProdPole* thePole = new LauKMatrixProdPole(poleName, iPole+1, resPairInt, propFull, theDaughters, useProdAdler); poleVect.push_back(thePole); } // Argand diagram graphs (full plot, then others showing only selected points) std::vector argGraphVect; // Real, imaginary and magnitude amplitude multigraphs std::vector ampGraphVect; // TLegend for argand graphs TLegend* argLegend = new TLegend(0.1, 0.1, 0.7, 0.995); argLegend->SetHeader("m(#pi^{+}#pi^{-}) points (GeV/c^{2})"); argLegend->SetTextSize(0.075); Bool_t initArgLegend(kFALSE); // Legend for the graphs (same for all multigraphs) TLegend* ampLegend = new TLegend(0.2, 0.2, 0.8, 0.7); ampLegend->SetHeader("Production pole components"); Bool_t initLegend(kFALSE); for (iPole = 0; iPole < nPoles; iPole++) { // Argand plot showing phase motion TMultiGraph* argGraphs = new TMultiGraph(); TString argMultiLabel("PoleMultiArg_"); argMultiLabel += iPole; argGraphs->SetName(argMultiLabel.Data()); argGraphVect.push_back(argGraphs); TGraph* argGraph = new TGraph(nSteps); TString argLabel("PoleArg_"); argLabel += iPole; argGraph->SetName(argLabel.Data()); argGraph->SetLineWidth(1); argGraphs->Add(argGraph, "lp"); // Argand graphs showing only one highlighted point for (Int_t iG = 0; iG < NMain; iG++) { TGraph* pointGraph = new TGraph(1); TString pointLabel("PolePoint_"); pointLabel += iG; pointGraph->SetName(pointLabel.Data()); pointGraph->SetMarkerColor(colours[iG]); pointGraph->SetMarkerSize(1.1); pointGraph->SetMarkerStyle(styles[iG]); argGraphs->Add(pointGraph, "p"); } if (!initArgLegend) { TList* argList = argGraphs->GetListOfGraphs(); for (Int_t iG = 0; iG < NMain; iG++) { TGraph* pointGraph = dynamic_cast(argList->At(iG+1)); argLegend->AddEntry(pointGraph, Form("%.2f",massPoints[iG]), "p"); initArgLegend = kTRUE; } } // Real, imaginary and magnitude amplitude TMultiGraph* ampGraphs = new TMultiGraph(); TString multiLabel("PoleMulti_"); multiLabel += iPole; ampGraphs->SetName(multiLabel.Data()); ampGraphVect.push_back(ampGraphs); TGraph* realGraph = new TGraph(nSteps); TString realLabel("PoleReal_"); realLabel += iPole; realGraph->SetName(realLabel.Data()); realGraph->SetLineColor(kRed); realGraph->SetMarkerColor(kRed); realGraph->SetLineWidth(1); ampGraphs->Add(realGraph, "lp"); TGraph* imagGraph = new TGraph(nSteps); TString imagLabel("PoleImag_"); imagLabel += iPole; imagGraph->SetName(imagLabel.Data()); imagGraph->SetLineColor(kBlue); imagGraph->SetMarkerColor(kBlue); imagGraph->SetLineWidth(1); ampGraphs->Add(imagGraph, "lp"); TGraph* magGraph = new TGraph(nSteps); TString magLabel("PoleMag_"); magLabel += iPole; magGraph->SetName(magLabel.Data()); magGraph->SetLineColor(kBlack); magGraph->SetMarkerColor(kBlack); magGraph->SetLineWidth(1); ampGraphs->Add(magGraph, "lp"); if (!initLegend) { ampLegend->AddEntry(realGraph, "Real part", "l"); ampLegend->AddEntry(imagGraph, "Imag part", "l"); ampLegend->AddEntry(magGraph, "Magnitude", "l"); initLegend = kTRUE; } } for (Int_t iM = 0; iM < nSteps; iM++) { Double_t m = dm*iM + mMin; Double_t s = m*m; // Update kinematics. The K-matrix only really needs s = m23Sq, i.e. m13Sq is "ignored". // Here, take the diagonal where m12Sq = 0.5*(mB^2 + 3mpi^2), the DP centre and use // m13Sq = mBSq + 3mPiSq - m12Sq - m23Sq Double_t m13Sq = mSqTerm - s; kinematics->updateKinematics(m13Sq, s); for (iPole = 0; iPole < nPoles; iPole++) { // Get the amplitude for the ith pole LauComplex amp = poleVect[iPole]->amplitude(kinematics); Double_t reAmp = amp.re(); Double_t imAmp = amp.im(); Double_t magAmp = sqrt(reAmp*reAmp + imAmp*imAmp); // Argand diagrams: full plot and only selected points TMultiGraph* argGraphs = argGraphVect[iPole]; TList* argList = argGraphs->GetListOfGraphs(); TGraph* argGraph = dynamic_cast(argList->At(0)); argGraph->SetPoint(iM, reAmp, imAmp); // Graphs showing specific points for (Int_t iP = 0; iP < NMain; iP++) { if (points[iP] == iM) { TGraph* pointGraph = dynamic_cast(argList->At(iP+1)); pointGraph->SetPoint(0, reAmp, imAmp); } } // Amplitude multigraphs TMultiGraph* ampGraphs = ampGraphVect[iPole]; TList* graphList = ampGraphs->GetListOfGraphs(); TGraph* realGraph = dynamic_cast(graphList->At(0)); realGraph->SetPoint(iM, m, reAmp); TGraph* imagGraph = dynamic_cast(graphList->At(1)); imagGraph->SetPoint(iM, m, imAmp); TGraph* magGraph = dynamic_cast(graphList->At(2)); magGraph->SetPoint(iM, m, magAmp); } } // Create plots TCanvas* theCanvas = new TCanvas("theCanvas", "", 1000, 1000); gROOT->SetStyle("Plain"); gStyle->SetOptStat(0); theCanvas->Clear(); theCanvas->UseCurrentStyle(); theCanvas->Divide(2,3); // Argand plots for (iPole = 0; iPole < nPoles; iPole++) { theCanvas->cd(iPole+1); TMultiGraph* argGraphs = argGraphVect[iPole]; argGraphs->Draw("a"); TString title("Production Pole "); title += iPole+1; argGraphs->SetTitle(title.Data()); TAxis* xAxis = argGraphs->GetXaxis(); xAxis->SetTitle("Real"); xAxis->SetLabelSize(0.05); xAxis->SetTitleSize(0.05); xAxis->CenterTitle(kTRUE); xAxis->SetTitleOffset(0.9); TAxis* yAxis = argGraphs->GetYaxis(); yAxis->SetTitle("Imag"); yAxis->SetLabelSize(0.05); yAxis->SetTitleSize(0.05); yAxis->CenterTitle(kTRUE); theCanvas->Update(); } theCanvas->cd(6); argLegend->Draw(); theCanvas->Print("ProdPoleArgand.png"); if (doEPS) {theCanvas->Print("ProdPoleArgand.eps");} // Amplitude plots theCanvas->Clear(); theCanvas->UseCurrentStyle(); theCanvas->Divide(2,3); for (iPole = 0; iPole < nPoles; iPole++) { theCanvas->cd(iPole+1); TMultiGraph* ampGraphs = ampGraphVect[iPole]; ampGraphs->Draw("a"); TAxis* xAxis = ampGraphs->GetXaxis(); xAxis->SetTitle("m(#pi^{+}#pi^{-}) (GeV/c^{2})"); xAxis->SetLabelSize(0.05); xAxis->SetTitleSize(0.05); xAxis->CenterTitle(kTRUE); xAxis->SetTitleOffset(0.9); TAxis* yAxis = ampGraphs->GetYaxis(); TString yLabel("Production Pole "); yLabel += iPole+1; yAxis->SetTitleOffset(0.9); yAxis->SetTitle(yLabel.Data()); yAxis->SetLabelSize(0.05); yAxis->SetTitleSize(0.05); yAxis->CenterTitle(kTRUE); theCanvas->Update(); } // Legend theCanvas->cd(6); ampLegend->Draw(); theCanvas->Print("ProdPoleAmp.png"); if (doEPS) {theCanvas->Print("ProdPoleAmp.eps");} // Cleanup std::vector::iterator argIter; for (argIter = argGraphVect.begin(); argIter != argGraphVect.end(); ++argIter) { TMultiGraph* theGraph = *argIter; delete theGraph; } std::vector::iterator ampIter; for (ampIter = ampGraphVect.begin(); ampIter != ampGraphVect.end(); ++ampIter) { TMultiGraph* theGraph = *ampIter; delete theGraph; } std::vector::iterator poleIter; for (poleIter = poleVect.begin(); poleIter != poleVect.end(); ++poleIter) { LauKMatrixProdPole* thePole = *poleIter; delete thePole; } delete theDaughters; delete propFull; delete ampLegend; delete theCanvas; } void doProdSVPPlots(Bool_t useProdAdler) { // Create plots of the SVP production amplitudes using the first row // of the propagator (for pipi). Each SVP is a pipi -> channel mode Int_t resPairInt(1); Int_t nChannels(5); Int_t KMatrixIndex(1); // pi pi S-wave Int_t nSVPs(5); // Mass range for plots Double_t mMin(0.28); Double_t mMax(2.0); Int_t nSteps = 2000; Double_t dm = (mMax - mMin)/(nSteps*1.0); // Find the point numbers for m = mMin, 0.5, 0.75, 1, 1.25, 1.5, 1.75 and 2 GeV std::vector points; std::vector colours, styles; std::vector massPoints; massPoints.push_back(mMin); massPoints.push_back(0.5); massPoints.push_back(0.75); massPoints.push_back(1.0); massPoints.push_back(1.25); massPoints.push_back(1.5); massPoints.push_back(1.75); massPoints.push_back(1.99999); Int_t NMain = static_cast(massPoints.size()); for (Int_t iP = 0; iP < NMain; iP++) { Int_t pointInt = static_cast((massPoints[iP] - mMin)/dm); points.push_back(pointInt); // Red or Blue Int_t iColour = kRed; if (iP > 3) {iColour = kBlue;} Int_t iStyle = kFullCircle; Int_t iPMod = iP%4; if (iPMod == 1) { iStyle = kOpenCircle; } else if (iPMod == 2) { iStyle = kFullSquare; } else if (iPMod == 3) { iStyle = kOpenSquare; } colours.push_back(iColour); styles.push_back(iStyle); } Double_t mPi = 0.13957018; Double_t mPiSq = mPi*mPi; Double_t mB = 5.279; Double_t mBSq = mB*mB; Double_t mSqTerm = 2.0*((mBSq/3.0) + mPiSq); // Full K-matrix parameterisation - TString KFile("B3piKMatrixCoeff.dat"); + const TString KFile { "@CMAKE_INSTALL_PREFIX@/@CMAKE_INSTALL_DATADIR@/@CMAKE_PROJECT_NAME@/B3piKMatrixCoeff.json" }; LauKMatrixPropagator* propFull = new LauKMatrixPropagator("propFull", KFile, resPairInt, nChannels, nSVPs, KMatrixIndex); // Define kinematics pointer Bool_t squareDP = kTRUE; LauDaughters* theDaughters = new LauDaughters("B+", "pi+", "pi+", "pi-", squareDP); // Kinematics LauKinematics* kinematics = theDaughters->getKinematics(); std::vector SVPVect; Int_t iSVP(0); for (iSVP = 0; iSVP < nSVPs; iSVP++) { TString SVPName("SVP"); SVPName += iSVP; LauKMatrixProdSVP* theSVP = new LauKMatrixProdSVP(SVPName, iSVP+1, resPairInt, propFull, theDaughters, useProdAdler); SVPVect.push_back(theSVP); } // Argand diagram graphs (full plot, then others showing only selected points) std::vector argGraphVect; // Real, imaginary and magnitude amplitude multigraphs std::vector ampGraphVect; // TLegend for argand graphs TLegend* argLegend = new TLegend(0.1, 0.1, 0.7, 0.995); argLegend->SetHeader("m(#pi^{+}#pi^{-}) points (GeV/c^{2})"); argLegend->SetTextSize(0.075); Bool_t initArgLegend(kFALSE); // Legend for the graphs (same for all multigraphs) TLegend* ampLegend = new TLegend(0.2, 0.2, 0.8, 0.7); ampLegend->SetHeader("Production SVP components"); Bool_t initLegend(kFALSE); for (iSVP = 0; iSVP < nSVPs; iSVP++) { // Argand plot showing phase motion TMultiGraph* argGraphs = new TMultiGraph(); TString argMultiLabel("SVPMultiArg_"); argMultiLabel += iSVP; argGraphs->SetName(argMultiLabel.Data()); argGraphVect.push_back(argGraphs); TGraph* argGraph = new TGraph(nSteps); TString argLabel("SVPArg_"); argLabel += iSVP; argGraph->SetName(argLabel.Data()); argGraph->SetLineWidth(1); argGraphs->Add(argGraph, "lp"); // Argand graphs showing only one highlighted point for (Int_t iG = 0; iG < NMain; iG++) { TGraph* pointGraph = new TGraph(1); TString pointLabel("SVPPoint_"); pointLabel += iG; pointGraph->SetName(pointLabel.Data()); pointGraph->SetMarkerColor(colours[iG]); pointGraph->SetMarkerSize(1.1); pointGraph->SetMarkerStyle(styles[iG]); argGraphs->Add(pointGraph, "p"); } if (!initArgLegend) { TList* argList = argGraphs->GetListOfGraphs(); for (Int_t iG = 0; iG < NMain; iG++) { TGraph* pointGraph = dynamic_cast(argList->At(iG+1)); argLegend->AddEntry(pointGraph, Form("%.2f",massPoints[iG]), "p"); initArgLegend = kTRUE; } } // Real, imaginary and magnitude amplitude TMultiGraph* ampGraphs = new TMultiGraph(); TString multiLabel("SVPMulti_"); multiLabel += iSVP; ampGraphs->SetName(multiLabel.Data()); ampGraphVect.push_back(ampGraphs); TGraph* realGraph = new TGraph(nSteps); TString realLabel("SVPReal_"); realLabel += iSVP; realGraph->SetName(realLabel.Data()); realGraph->SetLineColor(kRed); realGraph->SetMarkerColor(kRed); realGraph->SetLineWidth(1); ampGraphs->Add(realGraph, "lp"); TGraph* imagGraph = new TGraph(nSteps); TString imagLabel("SVPImag_"); imagLabel += iSVP; imagGraph->SetName(imagLabel.Data()); imagGraph->SetLineColor(kBlue); imagGraph->SetMarkerColor(kBlue); imagGraph->SetLineWidth(1); ampGraphs->Add(imagGraph, "lp"); TGraph* magGraph = new TGraph(nSteps); TString magLabel("SVPMag_"); magLabel += iSVP; magGraph->SetName(magLabel.Data()); magGraph->SetLineColor(kBlack); magGraph->SetMarkerColor(kBlack); magGraph->SetLineWidth(1); ampGraphs->Add(magGraph, "lp"); if (!initLegend) { ampLegend->AddEntry(realGraph, "Real part", "l"); ampLegend->AddEntry(imagGraph, "Imag part", "l"); ampLegend->AddEntry(magGraph, "Magnitude", "l"); initLegend = kTRUE; } } for (Int_t iM = 0; iM < nSteps; iM++) { Double_t m = dm*iM + mMin; Double_t s = m*m; // Update kinematics. The K-matrix only really needs s = m23Sq, i.e. m13Sq is "ignored". // Here, take the diagonal where m12Sq = 0.5*(mB^2 + 3mpi^2), the DP centre and use // m13Sq = mBSq + 3mPiSq - m12Sq - m23Sq Double_t m13Sq = mSqTerm - s; kinematics->updateKinematics(m13Sq, s); for (iSVP = 0; iSVP < nSVPs; iSVP++) { // Get the amplitude for the ith SVP LauComplex amp = SVPVect[iSVP]->amplitude(kinematics); Double_t reAmp = amp.re(); Double_t imAmp = amp.im(); Double_t magAmp = sqrt(reAmp*reAmp + imAmp*imAmp); // Argand diagrams: full plot and only selected points TMultiGraph* argGraphs = argGraphVect[iSVP]; TList* argList = argGraphs->GetListOfGraphs(); TGraph* argGraph = dynamic_cast(argList->At(0)); argGraph->SetPoint(iM, reAmp, imAmp); // Graphs showing specific points for (Int_t iP = 0; iP < NMain; iP++) { if (points[iP] == iM) { TGraph* pointGraph = dynamic_cast(argList->At(iP+1)); pointGraph->SetPoint(0, reAmp, imAmp); } } // Amplitude multigraphs TMultiGraph* ampGraphs = ampGraphVect[iSVP]; TList* graphList = ampGraphs->GetListOfGraphs(); TGraph* realGraph = dynamic_cast(graphList->At(0)); realGraph->SetPoint(iM, m, reAmp); TGraph* imagGraph = dynamic_cast(graphList->At(1)); imagGraph->SetPoint(iM, m, imAmp); TGraph* magGraph = dynamic_cast(graphList->At(2)); magGraph->SetPoint(iM, m, magAmp); } } // Create plots TCanvas* theCanvas = new TCanvas("theCanvas", "", 1000, 1000); gROOT->SetStyle("Plain"); gStyle->SetOptStat(0); theCanvas->Clear(); theCanvas->UseCurrentStyle(); theCanvas->Divide(2,3); // Argand plots for (iSVP = 0; iSVP < nSVPs; iSVP++) { theCanvas->cd(iSVP+1); TMultiGraph* argGraphs = argGraphVect[iSVP]; argGraphs->Draw("a"); TString title("Production SVP "); title += iSVP+1; argGraphs->SetTitle(title.Data()); TAxis* xAxis = argGraphs->GetXaxis(); xAxis->SetTitle("Real"); xAxis->SetLabelSize(0.05); xAxis->SetTitleSize(0.05); xAxis->CenterTitle(kTRUE); xAxis->SetTitleOffset(0.9); TAxis* yAxis = argGraphs->GetYaxis(); yAxis->SetTitle("Imag"); yAxis->SetLabelSize(0.05); yAxis->SetTitleSize(0.05); yAxis->CenterTitle(kTRUE); theCanvas->Update(); } theCanvas->cd(6); argLegend->Draw(); theCanvas->Print("ProdSVPArgand.png"); if (doEPS) {theCanvas->Print("ProdSVPArgand.eps");} // Amplitude plots theCanvas->Clear(); theCanvas->UseCurrentStyle(); theCanvas->Divide(2,3); for (iSVP = 0; iSVP < nSVPs; iSVP++) { theCanvas->cd(iSVP+1); TMultiGraph* ampGraphs = ampGraphVect[iSVP]; ampGraphs->Draw("a"); TAxis* xAxis = ampGraphs->GetXaxis(); xAxis->SetTitle("m(#pi^{+}#pi^{-}) (GeV/c^{2})"); xAxis->SetLabelSize(0.05); xAxis->SetTitleSize(0.05); xAxis->CenterTitle(kTRUE); xAxis->SetTitleOffset(0.9); TAxis* yAxis = ampGraphs->GetYaxis(); TString yLabel("Production SVP "); yLabel += iSVP+1; yAxis->SetTitleOffset(0.9); yAxis->SetTitle(yLabel.Data()); yAxis->SetLabelSize(0.05); yAxis->SetTitleSize(0.05); yAxis->CenterTitle(kTRUE); theCanvas->Update(); } // Legend theCanvas->cd(6); ampLegend->Draw(); theCanvas->Print("ProdSVPAmp.png"); if (doEPS) {theCanvas->Print("ProdSVPAmp.eps");} // Cleanup std::vector::iterator argIter; for (argIter = argGraphVect.begin(); argIter != argGraphVect.end(); ++argIter) { TMultiGraph* theGraph = *argIter; delete theGraph; } std::vector::iterator ampIter; for (ampIter = ampGraphVect.begin(); ampIter != ampGraphVect.end(); ++ampIter) { TMultiGraph* theGraph = *ampIter; delete theGraph; } std::vector::iterator SVPIter; for (SVPIter = SVPVect.begin(); SVPIter != SVPVect.end(); ++SVPIter) { LauKMatrixProdSVP* theSVP = *SVPIter; delete theSVP; } delete theDaughters; delete propFull; delete ampLegend; delete theCanvas; } int main(const int argc, const char** argv) { Int_t prodAdlerInt(0), epsInt(0); if (argc > 1) { prodAdlerInt = std::atoi(argv[1]); if (prodAdlerInt != 0) {doProdAdler = kTRUE;} } if (argc > 2) { epsInt = std::atoi(argv[2]); if (epsInt != 0) {doEPS = kTRUE;} } // Scattering K-matrix components doKScattPlots(); // First row of the propagator (pi pi channel) doPropagatorPlots(); // SVP and SVP production amplitudes (pi pi channel) doProdPolePlots(doProdAdler); doProdSVPPlots(doProdAdler); return 0; } diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 6ed679c..df151c1 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -1,83 +1,89 @@ add_subdirectory(ProgOpts) # List of files (without extensions) that need pre-processing before compiling list(APPEND EXAMPLE_SOURCES_CONFIGURE + B3piKMatrixMassProj + B3piKMatrixPlots KMatrixDto3pi + KMatrixExample + PlotKMatrixTAmp ) # List of files (without extensions) that do not need pre-processing before compiling list(APPEND EXAMPLE_SOURCES - B3piKMatrixMassProj - B3piKMatrixPlots CalcChiSq FlatPhaseSpace FlatPhaseSpace-CustomMasses FlatSqDalitz FlatSqDalitz-CustomMasses GenFit3K GenFit3KS GenFit3pi GenFitBelleCPKpipi GenFitDpipi GenFitDs2KKpi GenFitEFKLLM GenFitIncoherent_Bs2KSKpi GenFitKpipi GenFitNoDP GenFitNoDPMultiDim GenFitTimeDep GenFitTimeDep_Bs2KSKpi GenFitTimeDep_Bd2Dpipi - KMatrixExample MergeDataFiles mixedSampleTest - PlotKMatrixTAmp PlotResults point2PointTestSample QuasiFlatSqDalitz QuasiFlatSqDalitz-CustomMasses ResultsExtractor SimFitCoordinator SimFitTask SimFitTaskRooFit ) if(NOT LAURA_BUILD_ROOFIT_TASK) list(REMOVE_ITEM EXAMPLE_SOURCES SimFitTaskRooFit) endif() # Configure (if appropriate), build and install examples foreach( _example ${EXAMPLE_SOURCES_CONFIGURE}) configure_file(${_example}.cc.in ${_example}.cc @ONLY) add_executable(${_example} ${CMAKE_CURRENT_BINARY_DIR}/${_example}.cc) target_link_libraries(${_example} PRIVATE Laura++ ProgramOptions) install(TARGETS ${_example} RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) endforeach() foreach( _example ${EXAMPLE_SOURCES}) add_executable(${_example} ${_example}.cc) target_link_libraries(${_example} PRIVATE Laura++ ProgramOptions) install(TARGETS ${_example} RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) endforeach() # Also install the python script version of GenFit3pi configure_file(GenFit3pi.py.in GenFit3pi.py @ONLY) install(PROGRAMS ${CMAKE_CURRENT_BINARY_DIR}/GenFit3pi.py DESTINATION ${CMAKE_INSTALL_BINDIR}) # Also install the files to configure the models, coeffs, etc. # Some need pre-processing list(APPEND EXAMPLE_CONF_CONFIGURE KMatrixDto3pi.json ) foreach( _json ${EXAMPLE_CONF_CONFIGURE}) configure_file(${_json}.in ${_json} @ONLY) endforeach() list(APPEND EXAMPLE_CONF + B3piKMNoAdler.json + B3piKMPoles.json + B3piKMSVPs.json + B3piKMatrixCoeff.json + FOCUSD3pi.json GenFit3pi.json - FOCUSD3pi.dat + KMatrixKPiCoeff.json + KMatrixPiPiCoeff.json ) list(TRANSFORM EXAMPLE_CONF PREPEND ${CMAKE_CURRENT_SOURCE_DIR}/) list(TRANSFORM EXAMPLE_CONF_CONFIGURE PREPEND ${CMAKE_CURRENT_BINARY_DIR}/) list(APPEND EXAMPLE_CONF ${EXAMPLE_CONF_CONFIGURE}) install(FILES ${EXAMPLE_CONF} DESTINATION ${CMAKE_INSTALL_DATADIR}/${CMAKE_PROJECT_NAME}) diff --git a/examples/DToKspipiKMatrixCoeff.dat b/examples/DToKspipiKMatrixCoeff.dat deleted file mode 100644 index f975d62..0000000 --- a/examples/DToKspipiKMatrixCoeff.dat +++ /dev/null @@ -1,45 +0,0 @@ -# -# Copyright 2015 University of Warwick -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# -# Laura++ package authors: -# John Back -# Paul Harrison -# Thomas Latham -# -# File for pi-pi K-matrix S-wave amplitude. -# Coeffs from [hep-ex]0804.2089, table I, page 10 (RHS). -# First, define the indices of the phase-space channels defined by -# the LauKMatrixPropagator::KMatrixChannel enum -Channels 1 2 3 4 5 -# Next, define the bare pole masses and coupling constants g^(alpha)_j. -# First int = alpha (pole index number), 2nd number = pole mass (GeV/c^2), while -# the other numbers are constants for each channel j = 1 to N for the given pole -Pole 1 0.65100 0.22889 -0.55377 0.0 -0.39899 -0.34639 -Pole 2 1.20360 0.94128 0.55095 0.0 0.39065 0.31503 -Pole 3 1.55817 0.36856 0.23888 0.55639 0.18340 0.18681 -Pole 4 1.21000 0.33650 0.40907 0.85679 0.19906 -0.00984 -Pole 5 1.82206 0.18171 -0.17558 -0.79658 -0.00355 0.22358 -# Next, define the scattering f_{ij} constants. Here, the first integer is -# the scattering channel index, while the other numbers = f_{ij}, j = 1 to N -Scatt 1 0.23399 0.15044 -0.20545 0.32825 0.35412 -# We usually assume symmetry: f_{ji} = f_{ij}. Otherwise, uncomment this line -#ScattSymmetry 0 -# Finally, define the "Adler-zero" constants -# mSq0 s0Scatt s0Prod sA sA0 -mSq0 1.0 -s0Scatt -3.92637 -s0Prod -0.07 -sA 1.0 -sA0 -0.15 diff --git a/examples/FOCUSD3pi.dat b/examples/FOCUSD3pi.dat deleted file mode 100644 index d6584dc..0000000 --- a/examples/FOCUSD3pi.dat +++ /dev/null @@ -1,45 +0,0 @@ -# -# Copyright 2008 University of Warwick -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# -# Laura++ package authors: -# John Back -# Paul Harrison -# Thomas Latham -# -# File for pi-pi K-matrix S-wave amplitude. -# Coeffs from FOCUS D->3pi paper (hep-ex/0312040) -# First, define the indices of the phase-space channels defined by -# the LauKMatrixPropagator::KMatrixChannel enum -Channels 1 2 3 4 5 -# Next, define the bare pole masses and coupling constants g^(alpha)_j. -# First int = alpha (pole index number), 2nd number = pole mass (GeV/c^2), while -# the other numbers are constants for each channel j = 1 to N for the given pole -Pole 1 0.65100 0.24844 -0.52523 0.00000 -0.38878 -0.36397 -Pole 2 1.20720 0.91779 0.55427 0.00000 0.38705 0.29448 -Pole 3 1.56122 0.37024 0.23591 0.62605 0.18409 0.18923 -Pole 4 1.21257 0.34501 0.39642 0.97644 0.19746 0.00357 -Pole 5 1.81746 0.15770 -0.17915 -0.90100 -0.00931 0.20689 -# Next, define the scattering f_{ij} constants. Here, the first integer is -# the scattering channel index, while the other numbers = f_{ij}, j = 1 to N -Scatt 1 0.26681 0.16583 -0.19840 0.32808 0.31193 -# We usually assume symmetry: f_{ji} = f_{ij}. Otherwise, uncomment this line -#ScattSymmetry 0 -# Finally, define the "Adler-zero" constants -# mSq0 s0Scatt s0Prod sA sA0 -mSq0 1.0 -s0Scatt -3.30564 -s0Prod -1.0 -sA 1.0 -sA0 -0.2 diff --git a/examples/FOCUSD3pi.json b/examples/FOCUSD3pi.json new file mode 100644 index 0000000..a601a89 --- /dev/null +++ b/examples/FOCUSD3pi.json @@ -0,0 +1,59 @@ +{ + "comment" : "File for pi-pi K-matrix S-wave amplitude. Coeffs from FOCUS D->3pi paper (hep-ex/0312040).", + "channels" : [ "PiPi", "KK", "FourPi", "EtaEta", "EtaEtaP" ], + "poles" : [ + { + "index" : 1, + "mass" : 0.65100, + "couplings" : [ 0.24844, -0.52523, 0.0, -0.38878, -0.36397 ] + }, + { + "index" : 2, + "mass" : 1.20720, + "couplings" : [ 0.91779, 0.55427, 0.0, 0.38705, 0.29448 ] + }, + { + "index" : 3, + "mass" : 1.56122, + "couplings" : [ 0.37024, 0.23591, 0.62605, 0.18409, 0.18923 ] + }, + { + "index" : 4, + "mass" : 1.21257, + "couplings" : [ 0.34501, 0.39642, 0.97644, 0.19746, 0.00357 ] + }, + { + "index" : 5, + "mass" : 1.81746, + "couplings" : [ 0.15770, -0.17915, -0.90100, -0.00931, 0.20689 ] + } + ], + "scattering" : [ + { + "index" : 1, + "couplings" : [ 0.26681, 0.16583, -0.19840, 0.32808, 0.31193 ] + } + ], + "parameters" : [ + { + "name" : "mSq0", + "value" : 1.0 + }, + { + "name" : "s0Scatt", + "value" : -3.30564 + }, + { + "name" : "s0Prod", + "value" : -1.0 + }, + { + "name" : "sA", + "value" : 1.0 + }, + { + "name" : "sA0", + "value" : -0.2 + } + ] +} diff --git a/examples/KMatrixDto3pi.json.in b/examples/KMatrixDto3pi.json.in index bda8b80..3bae255 100644 --- a/examples/KMatrixDto3pi.json.in +++ b/examples/KMatrixDto3pi.json.in @@ -1,246 +1,246 @@ { "coeffs": { "coeffs": [ { "DeltaX": 0.0, "DeltaXBlind": false, "DeltaXFixed": true, "DeltaXSecondStage": false, "DeltaY": 0.0, "DeltaYBlind": false, "DeltaYFixed": true, "DeltaYSecondStage": false, "X": 1.0, "XBlind": false, "XFixed": true, "XSecondStage": false, "Y": 0.0, "YBlind": false, "YFixed": true, "YSecondStage": false, "clone": false, "name": "rho0(770)", "type": "CartesianCP" }, { "DeltaX": 0.0, "DeltaXBlind": false, "DeltaXFixed": true, "DeltaXSecondStage": false, "DeltaY": 0.0, "DeltaYBlind": false, "DeltaYFixed": true, "DeltaYSecondStage": false, "X": 0.45, "XBlind": false, "XFixed": false, "XSecondStage": false, "Y": 0.35, "YBlind": false, "YFixed": false, "YSecondStage": false, "clone": false, "name": "f_2(1270)", "type": "CartesianCP" }, { "DeltaX": 0.0, "DeltaXBlind": false, "DeltaXFixed": true, "DeltaXSecondStage": false, "DeltaY": 0.0, "DeltaYBlind": false, "DeltaYFixed": true, "DeltaYSecondStage": false, "X": 0.85, "XBlind": false, "XFixed": false, "XSecondStage": false, "Y": -0.55, "YBlind": false, "YFixed": false, "YSecondStage": false, "clone": false, "name": "KMPole1", "type": "CartesianCP" }, { "DeltaX": 0.0, "DeltaXBlind": false, "DeltaXFixed": true, "DeltaXSecondStage": false, "DeltaY": 0.0, "DeltaYBlind": false, "DeltaYFixed": true, "DeltaYSecondStage": false, "X": -1.1, "XBlind": false, "XFixed": false, "XSecondStage": false, "Y": -0.3, "YBlind": false, "YFixed": false, "YSecondStage": false, "clone": false, "name": "KMPole2", "type": "CartesianCP" }, { "DeltaX": 0.0, "DeltaXBlind": false, "DeltaXFixed": true, "DeltaXSecondStage": false, "DeltaY": 0.0, "DeltaYBlind": false, "DeltaYFixed": true, "DeltaYSecondStage": false, "X": -0.84, "XBlind": false, "XFixed": false, "XSecondStage": false, "Y": 0.35, "YBlind": false, "YFixed": false, "YSecondStage": false, "clone": false, "name": "KMPole3", "type": "CartesianCP" }, { "DeltaX": 0.0, "DeltaXBlind": false, "DeltaXFixed": true, "DeltaXSecondStage": false, "DeltaY": 0.0, "DeltaYBlind": false, "DeltaYFixed": true, "DeltaYSecondStage": false, "X": 0.23, "XBlind": false, "XFixed": false, "XSecondStage": false, "Y": 0.352, "YBlind": false, "YFixed": false, "YSecondStage": false, "clone": false, "name": "KMSVP1", "type": "CartesianCP" }, { "DeltaX": 0.0, "DeltaXBlind": false, "DeltaXFixed": true, "DeltaXSecondStage": false, "DeltaY": 0.0, "DeltaYBlind": false, "DeltaYFixed": true, "DeltaYSecondStage": false, "X": 0.33, "XBlind": false, "XFixed": false, "XSecondStage": false, "Y": 0.42, "YBlind": false, "YFixed": false, "YSecondStage": false, "clone": false, "name": "KMSVP2", "type": "CartesianCP" } ], "nCoeffs": 7 }, "model": { "commonSettings": { "setBWBachelorRestFrame": "ResonanceFrame", "setBWType": "BWPrimeBarrier", "setDefaultBWRadius": [ { "category": "Parent", "fix": true, "value": 4.0 }, { "category": "Light", "fix": true, "value": 5.3 } ], "setSpinFormalism": "Zemach_P" }, "kmatrix": [ { "poles": [ { "poleIndex": 1, "poleName": "KMPole1" }, { "poleIndex": 2, "poleName": "KMPole2" }, { "poleIndex": 3, "poleName": "KMPole3" } ], "propagator": { "nChannels": 5, "nPoles": 5, - "paramFileName": "@CMAKE_INSTALL_PREFIX@/@CMAKE_INSTALL_DATADIR@/@CMAKE_PROJECT_NAME@/FOCUSD3pi.dat", + "paramFileName": "@CMAKE_INSTALL_PREFIX@/@CMAKE_INSTALL_DATADIR@/@CMAKE_PROJECT_NAME@/FOCUSD3pi.json", "propName": "KMSWave", "resPairAmpInt": 1 }, "svps": [ { "channelIndex": 1, "svpName": "KMSVP1" }, { "channelIndex": 2, "svpName": "KMSVP2" } ] } ], "resonances": [ { "changeResonance": { "mass": { "fix": true, "value": 0.77526 }, "spin": { "value": 1 }, "width": { "fix": true, "value": 0.1478 } }, "resName": "rho0(770)", "resPairAmpInt": 1, "resType": "RelBW" }, { "changeResonance": { "mass": { "fix": true, "value": 1.2755 }, "spin": { "value": 2 }, "width": { "fix": true, "value": 0.1867 } }, "resName": "f_2(1270)", "resPairAmpInt": 1, "resType": "RelBW" } ] } } diff --git a/examples/KMatrixExample.cc b/examples/KMatrixExample.cc.in similarity index 95% rename from examples/KMatrixExample.cc rename to examples/KMatrixExample.cc.in index ef5031b..65eee56 100644 --- a/examples/KMatrixExample.cc +++ b/examples/KMatrixExample.cc.in @@ -1,219 +1,220 @@ /* Copyright 2008 University of Warwick Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0 Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License. */ /* Laura++ package authors: John Back Paul Harrison Thomas Latham */ #include using std::cout; using std::cerr; using std::endl; #include #include #include "LauDaughters.hh" #include "LauVetoes.hh" #include "LauEffModel.hh" #include "LauIsobarDynamics.hh" #include "LauSimpleFitModel.hh" #include "LauAbsCoeffSet.hh" #include "LauRealImagCoeffSet.hh" #include "TString.h" void usage( std::ostream& out, const TString& progName ) { out<<"Usage:\n"; out< [nExpt = 1] [firstExpt = 0]"< [nExpt = 1] [firstExpt = 0] if ( argc < 2 ) { usage( std::cerr, argv[0] ); return EXIT_FAILURE; } TString command = argv[1]; command.ToLower(); Int_t iFit(0); Int_t nExpt(1); Int_t firstExpt(0); if ( command == "gen" ) { if ( argc > 2 ) { nExpt = atoi( argv[2] ); if ( argc > 3 ) { firstExpt = atoi( argv[3] ); } } } else if ( command == "fit" ) { if ( argc < 3 ) { usage( std::cerr, argv[0] ); return EXIT_FAILURE; } iFit = atoi( argv[2] ); if ( argc > 3 ) { nExpt = atoi( argv[3] ); if ( argc > 4 ) { firstExpt = atoi( argv[4] ); } } } else { usage( std::cerr, argv[0] ); return EXIT_FAILURE; } // If you want to use square DP histograms for efficiency, // backgrounds or you just want the square DP co-ordinates // stored in the toy MC ntuple then set this to kTRUE Bool_t squareDP = kFALSE; // This defines the DP => decay is B0 -> pi+ pi- K_S0 // Particle 1 = pi+ // Particle 2 = pi- // Particle 3 = K_S0 // The DP is defined in terms of m13Sq and m23Sq LauDaughters* daughters = new LauDaughters("B0", "pi+", "pi-", "K_S0", squareDP); LauVetoes* vetoes = new LauVetoes(); LauEffModel* effModel = new LauEffModel(daughters, vetoes); LauIsobarDynamics* sigModel = new LauIsobarDynamics(daughters, effModel); // Add relativistic BWs for the rho and K* resonances sigModel->addResonance("rho0(770)", 3, LauAbsResonance::ResonanceModel::RelBW); sigModel->addResonance("K*+(892)", 2, LauAbsResonance::ResonanceModel::RelBW); sigModel->addResonance("K*+_2(1430)", 2, LauAbsResonance::ResonanceModel::RelBW); // Define the pi+ pi- "S-wave" K-matrix propagator // The defineKMatrixPropagator requires four arguments // name of propagator, parameter file, mass pair index (1,2,3) for "s" variable, number of channels, number of // mass poles and the row index defining the K-matrix state (e.g. index = 1 for the first row -> S-wave). - Int_t resPairInt = 3; - Int_t nChannels = 5; - Int_t nPoles = 5; - Int_t KMatrixIndex = 1; - sigModel->defineKMatrixPropagator("KM_pipi", "LauKMatrixCoeff.dat", resPairInt, nChannels, nPoles, KMatrixIndex); + const UInt_t resPairInt = 3; + const UInt_t nChannels = 5; + const UInt_t nPoles = 5; + const UInt_t KMatrixIndex = 1; + const TString jsonFileName { "@CMAKE_INSTALL_PREFIX@/@CMAKE_INSTALL_DATADIR@/@CMAKE_PROJECT_NAME@/KMatrixPiPiCoeff.json" }; + sigModel->defineKMatrixPropagator("KM_pipi", jsonFileName, resPairInt, nChannels, nPoles, KMatrixIndex); // Now define the pole and "slowly-varying part" (SVP) amplitudes for this K-matrix. // Each part will have their own complex coefficients which can (should) be fitted. // The addKMatrixProd functions need: name of pole/SVP, name of propagator and the index // number of the amplitude (from 1 to nPoles or nChannels) sigModel->addKMatrixProdPole("KM_pipi_Pole1", "KM_pipi", 1); sigModel->addKMatrixProdPole("KM_pipi_Pole2", "KM_pipi", 2); sigModel->addKMatrixProdPole("KM_pipi_Pole3", "KM_pipi", 3); sigModel->addKMatrixProdPole("KM_pipi_Pole4", "KM_pipi", 4); sigModel->addKMatrixProdPole("KM_pipi_Pole5", "KM_pipi", 5); sigModel->addKMatrixProdSVP("KM_pipi_SVP1", "KM_pipi", 1); sigModel->addKMatrixProdSVP("KM_pipi_SVP2", "KM_pipi", 2); sigModel->addKMatrixProdSVP("KM_pipi_SVP3", "KM_pipi", 3); sigModel->setASqMaxValue(3.1); LauSimpleFitModel* fitModel = new LauSimpleFitModel(sigModel); std::vector> coeffset; // Define the coefficients for the P-wave rho and K* resonances coeffset.push_back(std::make_unique("rho0(770)", 1.00, 0.00, kTRUE, kTRUE)); coeffset.push_back(std::make_unique("K*+(892)", 0.97, -1.10, kFALSE, kFALSE)); coeffset.push_back(std::make_unique("K*+_2(1430)", 0.38, -0.54, kFALSE, kFALSE)); // Define the coefficients for the K-matrix amplitudes defined above ("beta" and "f" production coeffs) coeffset.push_back(std::make_unique("KM_pipi_Pole1", -0.18, -0.91, kFALSE, kFALSE)); coeffset.push_back(std::make_unique("KM_pipi_Pole2", -1.02, -0.38, kFALSE, kFALSE)); coeffset.push_back(std::make_unique("KM_pipi_Pole3", -0.37, 0.50, kFALSE, kFALSE)); coeffset.push_back(std::make_unique("KM_pipi_Pole4", -0.01, 0.92, kFALSE, kFALSE)); coeffset.push_back(std::make_unique("KM_pipi_Pole5", 0.01, -0.05, kFALSE, kFALSE)); coeffset.push_back(std::make_unique("KM_pipi_SVP1", 0.22, 0.76, kFALSE, kFALSE)); coeffset.push_back(std::make_unique("KM_pipi_SVP2", 0.89, -0.31, kFALSE, kFALSE)); coeffset.push_back(std::make_unique("KM_pipi_SVP3", -0.97, 0.93, kFALSE, kFALSE)); for ( auto& coeff : coeffset ) { fitModel->setAmpCoeffSet( std::move(coeff) ); } Double_t nSigEvents = 5000; Bool_t fixSigEvents = kTRUE; LauParameter* nSig = new LauParameter("signalEvents",nSigEvents,-1.0*nSigEvents,2.0*nSigEvents,fixSigEvents); fitModel->setNSigEvents(nSig); fitModel->setNExpts(nExpt, firstExpt); // Configure various fit options // Switch on/off calculation of asymmetric errors. fitModel->useAsymmFitErrors(kFALSE); // Randomise initial fit values for the signal mode fitModel->useRandomInitFitPars(kTRUE); const Bool_t haveBkgnds = ( fitModel->nBkgndClasses() > 0 ); // Switch on/off Poissonian smearing of total number of events fitModel->doPoissonSmearing(haveBkgnds); // Switch on/off Extended ML Fit option fitModel->doEMLFit(haveBkgnds); // Generate toy from the fitted parameters //TString fitToyFileName("fitToyMC_KMatrixExample_"); //fitToyFileName += iFit; //fitToyFileName += ".root"; //fitModel->compareFitData(100, fitToyFileName); // Write out per-event likelihoods and sWeights //TString splotFileName("splot_KMatrixExample_"); //splotFileName += iFit; //splotFileName += ".root"; //fitModel->writeSPlotData(splotFileName, "splot", kFALSE); // Set the names of the files to read/write TString dataFile("gen-KMatrixExample.root"); TString treeName("genResults"); TString rootFileName(""); TString tableFileName(""); if (command == "fit") { rootFileName = "fitKMatrixExample_"; rootFileName += iFit; rootFileName += "_expt_"; rootFileName += firstExpt; rootFileName += "-"; rootFileName += (firstExpt+nExpt-1); rootFileName += ".root"; tableFileName = "fitKMatrixExampleResults_"; tableFileName += iFit; } else { rootFileName = "dummy.root"; tableFileName = "genKMatrixExampleResults"; } // Execute the generation/fit fitModel->run( command, dataFile, treeName, rootFileName, tableFileName ); return EXIT_SUCCESS; } diff --git a/examples/KMatrixKPiCoeff.json b/examples/KMatrixKPiCoeff.json new file mode 100644 index 0000000..2ea1b99 --- /dev/null +++ b/examples/KMatrixKPiCoeff.json @@ -0,0 +1,49 @@ +{ + "comment" : "File for K-pi K-matrix S-wave amplitude. Coeffs are purely guesses.", + "channels" : [ "KPi", "KEtaP", "KThreePi" ], + "poles" : [ + { + "index" : 1, + "mass" : 1.235, + "couplings" : [ 1.739, 1.0, 0.112 ] + }, + { + "index" : 2, + "mass" : 1.810, + "couplings" : [ 0.454, 0.8, 0.391 ] + }, + { + "index" : 3, + "mass" : 1.947, + "couplings" : [ 0.493, 0.3, -0.368 ] + } + ], + "scattering" : [ + { + "index" : 1, + "couplings" : [ 0.880, 0.430, 0.170 ] + } + ], + "parameters" : [ + { + "name" : "mSq0", + "value" : 1.0 + }, + { + "name" : "s0Scatt", + "value" : -3.92637 + }, + { + "name" : "s0Prod", + "value" : -3.0 + }, + { + "name" : "sA", + "value" : 1.0 + }, + { + "name" : "sA0", + "value" : -0.15 + } + ] +} diff --git a/examples/KMatrixPiPiCoeff.json b/examples/KMatrixPiPiCoeff.json new file mode 100644 index 0000000..51240e4 --- /dev/null +++ b/examples/KMatrixPiPiCoeff.json @@ -0,0 +1,59 @@ +{ + "comment" : "File for pi-pi K-matrix S-wave amplitude. Coeffs from [hep-ex]0804.2089, table I, page 10 (RHS).", + "channels" : [ "PiPi", "KK", "FourPi", "EtaEta", "EtaEtaP" ], + "poles" : [ + { + "index" : 1, + "mass" : 0.65100, + "couplings" : [ 0.22889, -0.55377, 0.0, -0.39899, -0.34639 ] + }, + { + "index" : 2, + "mass" : 1.20360, + "couplings" : [ 0.94128, 0.55095, 0.0, 0.39065, 0.31503 ] + }, + { + "index" : 3, + "mass" : 1.55817, + "couplings" : [ 0.36856, 0.23888, 0.55639, 0.18340, 0.18681 ] + }, + { + "index" : 4, + "mass" : 1.21000, + "couplings" : [ 0.33650, 0.40907, 0.85679, 0.19906, -0.00984 ] + }, + { + "index" : 5, + "mass" : 1.82206, + "couplings" : [ 0.18171, -0.17558, -0.79658, -0.00355, 0.22358 ] + } + ], + "scattering" : [ + { + "index" : 1, + "couplings" : [ 0.23399, 0.15044, -0.20545, 0.32825, 0.35412 ] + } + ], + "parameters" : [ + { + "name" : "mSq0", + "value" : 1.0 + }, + { + "name" : "s0Scatt", + "value" : -3.92637 + }, + { + "name" : "s0Prod", + "value" : -0.07 + }, + { + "name" : "sA", + "value" : 1.0 + }, + { + "name" : "sA0", + "value" : -0.15 + } + ] +} diff --git a/examples/LauKMatrixCoeff.dat b/examples/LauKMatrixCoeff.dat deleted file mode 100644 index 178f981..0000000 --- a/examples/LauKMatrixCoeff.dat +++ /dev/null @@ -1,43 +0,0 @@ -# -# Copyright 2008 University of Warwick -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# -# Laura++ package authors: -# John Back -# Paul Harrison -# Thomas Latham -# -# File for pi-pi K-matrix S-wave amplitude. -# Coeffs from [hep-ex]0804.2089, table I, page 10 (RHS). -# First, define the indices of the phase-space channels defined by -# the LauKMatrixPropagator::KMatrixChannel enum -Channels 1 2 3 4 5 -# Next, define the bare pole masses and coupling constants g^(alpha)_j. -# First int = alpha (pole index number), 2nd number = pole mass (GeV/c^2), while -# the other numbers are constants for each channel j = 1 to N for the given pole -Pole 1 0.65100 0.22889 -0.55377 0.0 -0.39899 -0.34639 -Pole 2 1.20360 0.94128 0.55095 0.0 0.39065 0.31503 -Pole 3 1.55817 0.36856 0.23888 0.55639 0.18340 0.18681 -Pole 4 1.21000 0.33650 0.40907 0.85679 0.19906 -0.00984 -Pole 5 1.82206 0.18171 -0.17558 -0.79658 -0.00355 0.22358 -# Next, define the scattering f_{ij} constants. Here, the first integer is -# the scattering channel index, while the other numbers = f_{ij}, j = 1 to N -Scatt 1 0.23399 0.15044 -0.20545 0.32825 0.35412 -# Finally, define the "Adler-zero" constants -# mSq0 s0Scatt s0Prod sA sA0 -mSq0 1.0 -s0Scatt -3.92637 -s0Prod -0.07 -sA 1.0 -sA0 -0.15 diff --git a/examples/LauKMatrixCoeff2.dat b/examples/LauKMatrixCoeff2.dat deleted file mode 100644 index 5fb0a9b..0000000 --- a/examples/LauKMatrixCoeff2.dat +++ /dev/null @@ -1,40 +0,0 @@ -# -# Copyright 2008 University of Warwick -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# -# Laura++ package authors: -# John Back -# Paul Harrison -# Thomas Latham -# -# File for K-pi K-matrix S-wave amplitude. Values are guesses - need updating. -# First, define the indices of the phase-space channels defined by -# the LauKMatrixPropagator::KMatrixChannel enum -Channels 6 7 8 -# Next, define the bare pole masses and coupling constants g^(alpha)_j. -# First int = alpha (pole index number), 2nd number = pole mass (GeV/c^2), while -# the other numbers are constants for each channel j = 1 to N for the given pole -Pole 1 1.235 1.739 1.0 0.112 -Pole 2 1.810 0.454 0.8 0.391 -Pole 3 1.947 0.493 0.3 -0.368 -# Next, define the scattering f_{ij} constants. Here, the first integer is -# the scattering channel index, while the other numbers = f_{ij}, j = 1 to N -Scatt 1 0.880 0.430 0.170 -# Finally, define the "Alder-zero constants" -# mSq0 s0Scatt s0Prod sA sA0 -mSq0 1.0 -s0Scatt -3.92637 -s0Prod -3.0 -sA 1.0 -sA0 -0.15 diff --git a/examples/PlotKMatrixTAmp.cc b/examples/PlotKMatrixTAmp.cc.in similarity index 96% rename from examples/PlotKMatrixTAmp.cc rename to examples/PlotKMatrixTAmp.cc.in index 2128ab3..41e1861 100644 --- a/examples/PlotKMatrixTAmp.cc +++ b/examples/PlotKMatrixTAmp.cc.in @@ -1,229 +1,230 @@ /* Copyright 2014 University of Warwick Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0 Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License. */ /* Laura++ package authors: John Back Paul Harrison Thomas Latham */ // Code to plot the transition amplitude T for the K-matrix pi pi S-wave #include using std::cout; using std::endl; #include "LauKMatrixPropagator.hh" #include "LauComplex.hh" #include "TGraph.h" #include "TCanvas.h" #include "TStyle.h" #include "TROOT.h" #include "TSystem.h" #include "TAxis.h" #include "TLine.h" #include "TMath.h" #include "TLatex.h" #include "TArrow.h" int main(const int /*argc*/, const char** /*argv*/) { // Define the "S-wave" K-matrix propagator - Int_t nChannels = 5; - Int_t nPoles = 5; - Int_t resPairInt = 1; - Int_t KMatrixIndex = 1; // for S-wave - LauKMatrixPropagator* propagator = new LauKMatrixPropagator("KMSWave", "DToKspipiKMatrixCoeff.dat", + const UInt_t nChannels = 5; + const UInt_t nPoles = 5; + const UInt_t resPairInt = 1; + const UInt_t KMatrixIndex = 1; + const TString jsonFileName { "@CMAKE_INSTALL_PREFIX@/@CMAKE_INSTALL_DATADIR@/@CMAKE_PROJECT_NAME@/KMatrixPiPiCoeff.json" }; + LauKMatrixPropagator* propagator = new LauKMatrixPropagator("KMSWave", jsonFileName, resPairInt, nChannels, nPoles, KMatrixIndex); // Find the transition amplitude T for the pi pi mode for // invariant mass squared s between (2*mpi)^2 and 4.0 GeV^2 Double_t mMin(0.28); Double_t mMax(2.0); Int_t nS = 2000; Double_t dm = (mMax - mMin)/(nS*1.0); // Argand plot TGraph* TxyGraph = new TGraph(nS); // Intensity TGraph* TSqGraph = new TGraph(nS); // Phase TGraph* TDeltaGraph = new TGraph(nS); // Inelasticity eta TGraph* TEtaGraph = new TGraph(nS); Double_t radToDeg = 180.0/TMath::Pi(); // Keep track of the real value of T to know when it // changes sign from -ve to +ve => 180 phase shift Double_t oldReT(0.0); int n180Shift(0); for (Int_t i = 0; i < nS; i++) { Double_t m = dm*i + mMin; Double_t s = m*m; LauComplex TAmp = propagator->getTransitionAmp(s, KMatrixIndex); Double_t realT = TAmp.re(); Double_t imagT = TAmp.im(); // Argand diagram showing the transition amplitude phase motion TxyGraph->SetPoint(i, realT, imagT); // Transition amplitude squared TSqGraph->SetPoint(i, m, TAmp.abs2()); // T - E = 0.5*i, where E = inelasticity vector, pointing to T from (0,i/2) Double_t ReE = realT; Double_t ImE = imagT - 0.5; LauComplex EVect(ReE, ImE); // Inelasticity Double_t EMag = EVect.abs(); Double_t eta = 2.0*EMag; TEtaGraph->SetPoint(i, m, eta); // Find the phase shift angle, delta Double_t delta = 0.5*TMath::ATan2(ImE, fabs(ReE))*radToDeg + 45.0; if (realT < 0.0) {delta = 180.0 - delta;} // Have we gone through 180 deg (or 2*delta through 360 deg)? if (oldReT < 0.0 && realT > 0.0) {n180Shift += 1;} // Add 180 if required delta += 180.0*n180Shift; TDeltaGraph->SetPoint(i, m, delta); // Keep track of the real value of T to see if it changes sign oldReT = realT; } TCanvas* theCanvas = new TCanvas("theCanvas", "", 1400, 1000); gROOT->SetStyle("Plain"); gStyle->SetOptStat(0); theCanvas->Clear(); theCanvas->UseCurrentStyle(); theCanvas->Divide(2,2); // Argand plot theCanvas->cd(1); TxyGraph->SetTitle("Transition amplitude T#; S #equiv I + 2iT"); TxyGraph->GetXaxis()->SetTitle("Re(T)"); TxyGraph->GetYaxis()->SetTitle("Im(T)"); TxyGraph->GetXaxis()->CenterTitle(kTRUE); TxyGraph->GetYaxis()->CenterTitle(kTRUE); TxyGraph->SetLineWidth(1); TxyGraph->Draw("apl"); // Overlay the elastic boundaries TLine lineA1(-0.5, 0.0, -0.5, 1.0); lineA1.SetLineStyle(kDotted); lineA1.SetLineWidth(1); lineA1.Draw(); TLine lineA2(0.5, 0.0, 0.5, 1.0); lineA2.SetLineStyle(kDotted); lineA2.SetLineWidth(1); lineA2.Draw(); TLine lineA3(-0.5, 1.0, 0.5, 1.0); lineA3.SetLineStyle(kDotted); lineA3.SetLineWidth(1); lineA3.Draw(); // Intensity theCanvas->cd(2); TSqGraph->SetTitle("Transition amplitude intensity"); TSqGraph->GetXaxis()->SetTitle("m(#pi^{+}#pi^{-}) (GeV/c^{2})"); TSqGraph->GetYaxis()->SetTitle("|T|^{2}"); TSqGraph->SetLineWidth(1); TSqGraph->GetXaxis()->CenterTitle(kTRUE); TSqGraph->GetYaxis()->CenterTitle(kTRUE); TSqGraph->Draw("apl"); // Overlay the unitarity limit TAxis* xAxis = TSqGraph->GetXaxis(); TLine lineB1(xAxis->GetXmin(), 1.0, xAxis->GetXmax(), 1.0); lineB1.SetLineStyle(kDotted); lineB1.SetLineWidth(1); lineB1.Draw(); // Also overlay "interesting regions" TLatex text; text.SetTextSize(0.045); // sigma text.DrawLatex(0.65, 0.6, "f_{0}(500)"); text.DrawLatex(0.68, 0.525, "or #sigma"); // f0(980) text.SetTextSize(0.045); text.DrawLatex(0.9, 1.02, "f_{0}(980)"); TArrow arrow1(0.98, 1.0, 0.98, 0.8, 0.01); arrow1.SetLineWidth(1); arrow1.Draw(); // f0(1370) text.DrawLatex(1.25, 0.91, "f_{0}(1370)"); TArrow arrow2(1.37, 0.9, 1.37, 0.7, 0.01); arrow2.SetLineWidth(1); arrow2.Draw(); // f0(1500) text.DrawLatex(1.45, 0.76, "f_{0}(1500)"); TArrow arrow3(1.50, 0.72, 1.50, 0.52, 0.01); arrow3.SetLineWidth(1); arrow3.Draw(); // f0(1710) text.DrawLatex(1.60, 0.61, "f_{0}(1710)"); TArrow arrow4(1.71, 0.6, 1.71, 0.45, 0.01); arrow4.SetLineWidth(1); arrow4.Draw(); // Phase theCanvas->cd(3); TDeltaGraph->SetTitle("Phase shift"); TDeltaGraph->GetXaxis()->SetTitle("m(#pi^{+}#pi^{-}) (GeV/c^{2})"); TDeltaGraph->GetYaxis()->SetTitle("#delta (degrees)"); TDeltaGraph->GetXaxis()->CenterTitle(kTRUE); TDeltaGraph->GetYaxis()->CenterTitle(kTRUE); TDeltaGraph->SetLineWidth(1); TDeltaGraph->Draw("apl"); // Inelasticity theCanvas->cd(4); TEtaGraph->SetTitle("Inelasticity, #eta #equiv |2T - iI| = |S|"); TEtaGraph->GetXaxis()->SetTitle("m(#pi^{+}#pi^{-}) (GeV/c^{2})"); TEtaGraph->GetYaxis()->SetTitle("#eta"); TEtaGraph->GetXaxis()->CenterTitle(kTRUE); TEtaGraph->GetYaxis()->CenterTitle(kTRUE); TEtaGraph->SetLineWidth(1); TEtaGraph->Draw("apl"); // Overlay the elastic limit xAxis = TEtaGraph->GetXaxis(); TLine lineD1(xAxis->GetXmin(), 1.0, xAxis->GetXmax(), 1.0); lineD1.SetLineStyle(kDotted); lineD1.SetLineWidth(1); lineD1.Draw(); theCanvas->Print("KMatrixTAmpPlots.png"); theCanvas->Print("KMatrixTAmpPlots.eps"); } diff --git a/inc/LauKMatrixPropagator.hh b/inc/LauKMatrixPropagator.hh index 96ada5e..24fe7e7 100644 --- a/inc/LauKMatrixPropagator.hh +++ b/inc/LauKMatrixPropagator.hh @@ -1,657 +1,661 @@ /* Copyright 2008 University of Warwick Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0 Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License. */ /* Laura++ package authors: John Back Paul Harrison Thomas Latham */ /*! \file LauKMatrixPropagator.hh \brief File containing declaration of LauKMatrixPropagator class. */ /*! \class LauKMatrixPropagator \brief Class for defining a K-matrix propagator. Class used to define a K-matrix propagator. See the following papers for info: hep-ph/0204328, hep-ex/0312040, [hep-ex]0804.2089 and hep-ph/9705401. */ #ifndef LAU_KMATRIX_PROPAGATOR #define LAU_KMATRIX_PROPAGATOR #include "LauDatabasePDG.hh" #include "LauResonanceMaker.hh" #include "LauResonanceInfo.hh" #include "TMatrixD.h" #include "TString.h" #include #include #include class LauParameter; class LauKinematics; class LauComplex; class LauKMatrixPropagator { public: //! Constructor /*! \param [in] name name of the propagator \param [in] paramFileName the parameter file name \param [in] resPairAmpInt the number of the daughter not produced by the resonance \param [in] nChannels the number of channels \param [in] nPoles the number of poles \param [in] rowIndex this specifies which row of the propagator should be used when summing over the amplitude channels */ - LauKMatrixPropagator( const TString& name, const TString& paramFileName, - const Int_t resPairAmpInt, const Int_t nChannels, const Int_t nPoles, - const Int_t rowIndex = 1 ); + LauKMatrixPropagator( const TString& name, const TString& paramFileName, + const UInt_t resPairAmpInt, const UInt_t nChannels, const UInt_t nPoles, + const UInt_t rowIndex = 1 ); //! Destructor virtual ~LauKMatrixPropagator(); //! Calculate the K-matrix propagator for the given s value /*! \param [in] s the invariant mass squared */ void updatePropagator(const Double_t s); //! Read an input file to set parameters /*! \param [in] inputFile name of the input file */ void setParameters(const TString& inputFile); //! Set flag to ignore Blatt-Weisskopf-like barrier factor void ignoreBWBarrierFactor() {includeBWBarrierFactor_=kFALSE;} //! Get the scattering K matrix /*! \return the real, symmetric scattering K matrix */ TMatrixD getKMatrix() const {return ScattKMatrix_;} //! Get the real part of the propagator full matrix /*! \return the real part of the propagator full matrix */ TMatrixD getRealPropMatrix() const {return realProp_;} //! Get the negative imaginary part of the full propagator matrix /*! \return the negative imaginary part of the full propagator matrix */ TMatrixD getNegImagPropMatrix() const {return negImagProp_;} //! Get the real part of the term of the propagator /*! \param [in] channelIndex the channel number \return the real part of the propagator term */ - Double_t getRealPropTerm(const Int_t channelIndex) const; + Double_t getRealPropTerm(const UInt_t channelIndex) const; //! Get the imaginary part of the term of the propagator /*! \param [in] channelIndex the channel number \return the imaginiary part of the propagator term */ - Double_t getImagPropTerm(const Int_t channelIndex) const; + Double_t getImagPropTerm(const UInt_t channelIndex) const; //! Get the 1/(m_pole^2 -s) terms for the scattering and production K-matrix formulae /*! \param [in] poleIndex the number of the pole required \return the value of 1/(m_pole^2 -s) */ - Double_t getPoleDenomTerm(const Int_t poleIndex) const; + Double_t getPoleDenomTerm(const UInt_t poleIndex) const; //! Get spin of K-matrix /*! \param [in] iChannel the index of the channel \return the value of the orbital angular momentum, L_, for this channel */ - Int_t getL(const Int_t iChannel) const {return L_[iChannel];} + UInt_t getL(const UInt_t iChannel) const {return L_[iChannel];} //! Get index of final channel /*! \return the index of the channel into which the scattering happens */ - Int_t getIndex() const {return index_;}; + UInt_t getIndex() const {return index_;}; //! Get pole mass parameters, set according to the input file /*! \param [in] poleIndex number of the required pole \return the parameter of the pole mass */ - LauParameter& getPoleMassSqParameter(const Int_t poleIndex); + LauParameter& getPoleMassSqParameter(const UInt_t poleIndex); //! Get coupling constants that were loaded from the input file /*! \param [in] poleIndex number of the required pole \param [in] channelIndex number of the required channel \return the value of the coupling constant */ - Double_t getCouplingConstant(const Int_t poleIndex, const Int_t channelIndex) const; + Double_t getCouplingConstant(const UInt_t poleIndex, const UInt_t channelIndex) const; //! Get coupling parameters, set according to the input file /*! \param [in] poleIndex number of the required pole \param [in] channelIndex number of the required channel \return the parameter of the coupling constant */ - LauParameter& getCouplingParameter(const Int_t poleIndex, const Int_t channelIndex); + LauParameter& getCouplingParameter(const UInt_t poleIndex, const UInt_t channelIndex); //! Get scattering constants that were loaded from the input file /*! \param [in] channel1Index number of the first channel index \param [in] channel2Index number of the second channel index \return the value of the scattering constant */ - Double_t getScatteringConstant(const Int_t channel1Index, const Int_t channel2Index) const; + Double_t getScatteringConstant(const UInt_t channel1Index, const UInt_t channel2Index) const; //! Get scattering parameters, set according to the input file /*! \param [in] channel1Index number of the first channel index \param [in] channel2Index number of the second channel index \return the parameter of the scattering constant */ - LauParameter& getScatteringParameter(const Int_t channel1Index, const Int_t channel2Index); + LauParameter& getScatteringParameter(const UInt_t channel1Index, const UInt_t channel2Index); //! Get mSq0 production parameter /*! \return the mSq0 parameter */ LauParameter& getmSq0() {return mSq0_;} //! Get s0Scatt production parameter /*! \return the s0Scatt parameter */ LauParameter& gets0Scatt() {return s0Scatt_;} //! Get s0 production parameter /*! \return the s0Prod parameter */ LauParameter& gets0Prod() {return s0Prod_;} //! Get sA production parameter /*! \return the sA parameter */ LauParameter& getsA() {return sA_;} //! Get sA0 production parameter /*! \return the sA0 parameter */ LauParameter& getsA0() {return sA0_;} //! Get the "slowly-varying part" term of the amplitude /*! \return the svp term */ Double_t getProdSVPTerm() const {return prodSVP_;} //! Get the full complex propagator term for a given channel /*! \param [in] channelIndex the number of the required channel \return the complex propagator term */ - LauComplex getPropTerm(const Int_t channelIndex) const; + LauComplex getPropTerm(const UInt_t channelIndex) const; //! Get the DP axis identifier /*! \return the value to identify the DP axis in question */ - Int_t getResPairAmpInt() const {return resPairAmpInt_;} + UInt_t getResPairAmpInt() const {return resPairAmpInt_;} //! Get the number of channels /*! \return the number of channels */ - Int_t getNChannels() const {return nChannels_;} + UInt_t getNChannels() const {return nChannels_;} //! Get the number of poles /*! \return the number of poles */ - Int_t getNPoles() const {return nPoles_;} + UInt_t getNPoles() const {return nPoles_;} //! Get the propagator name /*! \return the name of the propagator */ TString getName() const {return name_;} //! Get the unitary transition amplitude for the given channel /*! \param [in] s The invariant mass squared \param [in] channel The index number of the channel process \return the complex amplitude T */ - LauComplex getTransitionAmp(const Double_t s, const Int_t channel); + LauComplex getTransitionAmp(const Double_t s, const UInt_t channel); //! Get the complex phase space term for the given channel and invariant mass squared /*! \param [in] s The invariant mass squared \param [in] channel The index number of the channel process \return the complex phase space term rho(channel, channel) */ - LauComplex getPhaseSpaceTerm(const Double_t s, const Int_t channel); + LauComplex getPhaseSpaceTerm(const Double_t s, const UInt_t channel); //! Get the Adler zero factor, which is set when updatePropagator is called /*! \return the Adler zero factor */ Double_t getAdlerZero() const {return adlerZeroFactor_;} //! Get the THat amplitude for the given s and channel number /*! \param [in] s The invariant mass squared \param [in] channel The index number of the channel process \return the complex THat amplitude */ - LauComplex getTHat(const Double_t s, const Int_t channel); + LauComplex getTHat(const Double_t s, const UInt_t channel); //! Create a JSON object containing the current settings nlohmann::json writeSettingsToJson() const; - protected: - // Integers to specify the allowed channels for the phase space calculations. - // Please keep Zero at the start and leave TotChannels at the end - // whenever more channels are added to this. - //! Integers to specify the allowed channels for the phase space calculations - enum class KMatrixChannels {Zero, PiPi, KK, FourPi, EtaEta, EtaEtaP, - KPi, KEtaP, KThreePi, D0K, Dstar0K, TotChannels}; + //! Enumeration to specify the allowed channels for the phase space calculations + enum class KMatrixChannels { + PiPi, + KK, + FourPi, + EtaEta, + EtaEtaP, + KPi, + KEtaP, + KThreePi, + D0K, + Dstar0K + }; + protected: //! Calculate the scattering K-matrix for the given value of s /*! \param [in] s the invariant mass squared */ void calcScattKMatrix(const Double_t s); //! Calculate the real and imaginary part of the phase space density diagonal matrix /*! \param [in] s the invariant mass squared */ void calcRhoMatrix(const Double_t s); //! Retrieve the complex phasespace density for a given channel /*! \param [in] s the invariant mass squared \param [in] phaseSpaceIndex the phasespace index of the channel \return the complex phasespace density */ LauComplex getRho(const Double_t s, const LauKMatrixPropagator::KMatrixChannels) const; //! Calculate the (real) gamma angular-momentum barrier matrix /*! \param [in] s the invariant mass squared */ void calcGammaMatrix(const Double_t s); //! Calculate the gamma angular-momentum barrier /*! \param [in] iCh the channel index \param [in] s the invariant mass squared \return the centrifugal barrier factor for L=0,1, or 2 */ - Double_t calcGamma(const Int_t iCh, const Double_t s) const; + Double_t calcGamma(const UInt_t iCh, const Double_t s) const; //! Calulate the term 1/(m_pole^2 - s) for the scattering and production K-matrix formulae /*! \param [in] s the invariant mass squared */ void calcPoleDenomVect(const Double_t s); //! Calculate the D0K+ phase space factor /*! \param [in] s the invariant mass squared \return the complex phase space factor */ LauComplex calcD0KRho(const Double_t s) const; //! Calculate the D*0K+ phase space factor /*! \param [in] s the invariant mass squared \return the complex phase space factor */ LauComplex calcDstar0KRho(const Double_t s) const; //! Calculate the pipi phase space factor /*! \param [in] s the invariant mass squared \return the complex phase space factor */ LauComplex calcPiPiRho(const Double_t s) const; //! Calculate the KK phase space factor /*! \param [in] s the invariant mass squared \return the complex phase space factor */ LauComplex calcKKRho(const Double_t s) const; //! Calculate the 4 pi phase space factor /*! \param [in] s the invariant mass squared \return the complex phase space factor */ LauComplex calcFourPiRho(const Double_t s) const; //! Calculate the eta-eta phase space factor /*! \param [in] s the invariant mass squared \return the complex phase space factor */ LauComplex calcEtaEtaRho(const Double_t s) const; //! Calculate the eta-eta' phase space factor /*! \param [in] s the invariant mass squared \return the complex phase space factor */ LauComplex calcEtaEtaPRho(const Double_t s) const; //! Calculate the Kpi phase space factor /*! \param [in] s the invariant mass squared \return the complex phase space factor */ LauComplex calcKPiRho(const Double_t s) const; //! Calculate the K-eta' phase space factor /*! \param [in] s the invariant mass squared \return the complex phase space factor */ LauComplex calcKEtaPRho(const Double_t s) const; //! Calculate the Kpipipi phase space factor /*! \param [in] s the invariant mass squared \return the complex phase space factor */ LauComplex calcKThreePiRho(const Double_t s) const; //! Calculate the "slow-varying part" /*! \param [in] s the invariant mass squared \param [in] s0 the invariant mass squared at the Adler zero \return the SVP term */ Double_t calcSVPTerm(const Double_t s, const Double_t s0) const; //! Update the scattering "slowly-varying part" /*! \param [in] s the invariant mass squared */ void updateScattSVPTerm(const Double_t s); //! Update the production "slowly-varying part" /*! \param [in] s the invariant mass squared */ void updateProdSVPTerm(const Double_t s); //! Calculate the multiplicative factor containing severa Adler zero constants /*! \param [in] s the invariant mass squared */ void updateAdlerZeroFactor(const Double_t s); - //! Check the phase space factors that need to be used - /*! - \param [in] phaseSpaceInt phase space types - \return true of false - */ - Bool_t checkPhaseSpaceType(const Int_t phaseSpaceInt) const; - - //! Get the unitary transition amplitude matrix for the given kinematics /*! \param [in] kinematics The pointer to the constant kinematics */ void getTMatrix(const LauKinematics* kinematics); //! Get the unitary transition amplitude matrix for the given kinematics /*! \param [in] s The invariant mass squared of the system */ void getTMatrix(const Double_t s); //! Get the square root of the phase space matrix void getSqrtRhoMatrix(); private: //! Copy constructor (not implemented) LauKMatrixPropagator(const LauKMatrixPropagator& rhs)=delete; //! Copy assignment operator (not implemented) LauKMatrixPropagator& operator=(const LauKMatrixPropagator& rhs)=delete; //! Initialise and set the dimensions for the internal matrices and parameter arrays void initialiseMatrices(); - //! Store the (phase space) channel indices from a line in the parameter file + //! Store the (phase space) channel indices from the parameter file /*! - \param [in] theLine Vector of strings corresponding to the line from the parameter file + \param [in] j "channels" element of the parameter file + \return whether the JSON could be successfully parsed */ - void storeChannels(const std::vector& theLine); + bool storeChannels(const nlohmann::json& j); - //! Store the pole mass and couplings from a line in the parameter file + //! Store the pole masses and couplings from the parameter file /*! - \param [in] theLine Vector of strings corresponding to the line from the parameter file + \param [in] j "poles" element of the parameter file + \return whether the JSON could be successfully parsed */ - void storePole(const std::vector& theLine); + bool storePoles(const nlohmann::json& j); - //! Store the scattering coefficients from a line in the parameter file + //! Store the scattering coefficients from the parameter file /*! - \param [in] theLine Vector of strings corresponding to the line from the parameter file + \param [in] j "scattering" element of the parameter file + \return whether the JSON could be successfully parsed */ - void storeScattering(const std::vector& theLine); + bool storeScattering(const nlohmann::json& j); //! Store the channels' characteristic radii from a line in the parameter file /*! - \param [in] theLine Vector of strings corresponding to the line from the parameter file + \param [in] j "radii" element of the parameter file + \return whether the JSON could be successfully parsed */ - void storeRadii(const std::vector& theLine); + bool storeRadii(const nlohmann::json& j); - //! Store the channels' angular momenta from a line in the parameter file + //! Store the channels' angular momenta from the parameter file /*! - \param [in] theLine Vector of strings corresponding to the line from the parameter file - \param [in] a Vector of integers corresponding to parameter in the propagator denominator + \param [in] j "angularmomentum" element of the parameter file + \param [in] a vector of integers corresponding to parameter in the propagator denominator + \return whether the JSON could be successfully parsed */ - void storeOrbitalAngularMomenta(const std::vector& theLine, std::vector& a); + bool storeOrbitalAngularMomenta(const nlohmann::json& j, std::vector& a); - //! Store the barrier-factor parameter from a line in the parameter file + //! Store the barrier-factor parameter from the parameter file /*! - \param [in] theLine Vector of strings corresponding to the line from the parameter file - \param [in] a Vector of integers corresponding to parameter in the propagator denominator + \param [in] j "barrierfactorparameter" element of the parameter file + \param [in] a vector of integers corresponding to parameter in the propagator denominator + \return whether the JSON could be successfully parsed */ - void storeBarrierFactorParameter(const std::vector& theLine, std::vector& a); - + bool storeBarrierFactorParameter(const nlohmann::json& j, std::vector& a); - //! Store miscelleanous parameters from a line in the parameter file + //! Store one of the miscelleanous parameters from the parameter file /*! - \param [in] keyword the name of the parameter to be set - \param [in] parString the string containing the value of the parameter + \param [in] j "parameters" element of the parameter file + \return whether the JSON could be successfully parsed */ - void storeParameter(const TString& keyword, const TString& parString); + bool storeParameter(const nlohmann::json& j); //! String to store the propagator name TString name_; //! Name of the input parameter file TString paramFileName_; //! Number to identify the DP axis in question - Int_t resPairAmpInt_; + UInt_t resPairAmpInt_; //! Row index - 1 - Int_t index_; + UInt_t index_; //! s value of the previous pole Double_t previousS_{0.0}; //! "slowly-varying part" for the scattering K-matrix Double_t scattSVP_{0.0}; //! "slowly-varying part" for the production K-matrix Double_t prodSVP_{0.0}; //! Real part of the propagator matrix TMatrixD realProp_; //! Imaginary part of the propagator matrix TMatrixD negImagProp_; //! Scattering K-matrix TMatrixD ScattKMatrix_; //! Real part of the phase space density diagonal matrix TMatrixD ReRhoMatrix_; //! Imaginary part of the phase space density diagonal matrix TMatrixD ImRhoMatrix_; //! Gamma angular-momentum barrier matrix TMatrixD GammaMatrix_; //! Identity matrix TMatrixD IMatrix_; //! Null matrix TMatrixD zeroMatrix_; //! Real part of the square root of the phase space density diagonal matrix TMatrixD ReSqrtRhoMatrix_; //! Imaginary part of the square root of the phase space density diagonal matrix TMatrixD ImSqrtRhoMatrix_; //! Real part of the unitary T matrix TMatrixD ReTMatrix_; //! Imaginary part of the unitary T matrix TMatrixD ImTMatrix_; //! Number of channels - Int_t nChannels_; + UInt_t nChannels_; //! Number of poles - Int_t nPoles_; + UInt_t nPoles_; //! Vector of orbital angular momenta - std::vector L_; + std::vector L_; //! Boolean to indicate whether storeBarrierFactorParameter has been called - Bool_t haveCalled_storeBarrierFactorParameter{kFALSE}; + Bool_t haveCalledStoreBarrierFactorParameter_{kFALSE}; //! Boolean to dictate whether to include Blatt-Weisskopf-like denominator in K-matrix centrifugal barrier factor Bool_t includeBWBarrierFactor_{kTRUE}; //! Vector of squared pole masses std::vector mSqPoles_; //! Array of coupling constants LauParArray gCouplings_; //! Array of scattering SVP values LauParArray fScattering_; //! Vector of characteristic radii std::vector radii_; //! Vector of ratio a/R^2 std::vector gamAInvRadSq_; //! Vector of phase space types std::vector phaseSpaceTypes_; //! Vector of squared masses std::vector mSumSq_; //! Vector of mass differences std::vector mDiffSq_; //! Vector of 1/(m_pole^2 - s) terms for scattering and production K-matrix formulae std::vector poleDenomVect_; //! Constant from input file LauParameter mSq0_; //! Constant from input file LauParameter s0Scatt_; //! Constant from input file LauParameter s0Prod_; //! Constant from input file LauParameter sA_; //! Constant from input file LauParameter sA0_; //! Defined as 0.5*sA*mPi*mPi Double_t sAConst_{0.0}; //! Charged pion mass const Double_t mPi_{LauDatabasePDG::particleMass("pi+")}; //! Charged pion mass squared const Double_t mPiSq_{mPi_*mPi_}; //! Charged kaon mass const Double_t mK_{LauDatabasePDG::particleMass("K+")}; //! Eta mass const Double_t mEta_{LauDatabasePDG::particleMass("eta")}; //! Eta-prime mass const Double_t mEtaPrime_{LauDatabasePDG::particleMass("eta'")}; //! Neutral D mass const Double_t mD0_{LauDatabasePDG::particleMass("D0")}; //! Neutral Dstar mass const Double_t mDstar0_{LauResonanceMaker::get().getResInfo("D*0")->getMass()->value()}; //! Defined as 4*mPi*mPi const Double_t m2piSq_{4.0*mPiSq_}; //! Defined as 4*mK*mK const Double_t m2KSq_{4.0*mK_*mK_}; //! Defined as 4*mEta*mEta const Double_t m2EtaSq_{4.0*mEta_*mEta_}; //! Defined as (mEta+mEta')^2 const Double_t mEtaEtaPSumSq_{(mEtaPrime_ + mEta_)*(mEtaPrime_ + mEta_)}; //! Defined as (mEta-mEta')^2 const Double_t mEtaEtaPDiffSq_{(mEtaPrime_ - mEta_)*(mEtaPrime_ - mEta_)}; //! Defined as (mK+mPi)^2 const Double_t mKpiSumSq_{(mK_ + mPi_)*(mK_ + mPi_)}; //! Defined as (mK-mPi)^2 const Double_t mKpiDiffSq_{(mK_ - mPi_)*(mK_ - mPi_)}; //! Defined as (mK+mEta')^2 const Double_t mKEtaPSumSq_{(mK_ + mEtaPrime_)*(mK_ + mEtaPrime_)}; //! Defined as (mK-mEta')^2 const Double_t mKEtaPDiffSq_{(mK_ - mEtaPrime_)*(mK_ - mEtaPrime_)}; //! Defined as (mK-3*mPi)^2 const Double_t mK3piDiffSq_{(mK_ - 3.0*mPi_)*(mK_ - 3.0*mPi_)}; //! Factor used to calculate the Kpipipi phase space term const Double_t k3piFactor_{TMath::Power((1.44 - mK3piDiffSq_)/1.44, -2.5)}; //! Factor used to calculate the pipipipi phase space term const Double_t fourPiFactor1_{16.0*mPiSq_}; //! Factor used to calculate the pipipipi phase space term const Double_t fourPiFactor2_{TMath::Sqrt(1.0 - fourPiFactor1_)}; //! Defined as (mD0+mK)^2 const Double_t mD0KSumSq_{(mD0_ + mK_)*(mD0_ + mK_)}; //! Defined as (mD0-mK)^2 const Double_t mD0KDiffSq_{(mD0_ - mK_)*(mD0_ - mK_)}; //! Defined as (mD*0+mK)^2 const Double_t mDstar0KSumSq_{(mDstar0_ + mK_)*(mDstar0_ + mK_)}; //! Defined as (mD*0-mK)^2 const Double_t mDstar0KDiffSq_{(mDstar0_ - mK_)*(mDstar0_ - mK_)}; //! Multiplicative factor containing various Adler zero constants Double_t adlerZeroFactor_{0.0}; //! Tracks if all params have been set Bool_t parametersSet_{kFALSE}; - //! Control the output of the functions - static constexpr Bool_t verbose_{kFALSE}; - //! Control if scattering constants are channel symmetric: f_ji = f_ij Bool_t scattSymmetry_{kFALSE}; + //! Control the output of the functions + static constexpr Bool_t verbose_{kFALSE}; + ClassDef(LauKMatrixPropagator,0) // K-matrix amplitude model }; #endif diff --git a/src/LauKMatrixPropagator.cc b/src/LauKMatrixPropagator.cc index 46b92cb..398a668 100644 --- a/src/LauKMatrixPropagator.cc +++ b/src/LauKMatrixPropagator.cc @@ -1,1447 +1,1474 @@ /* Copyright 2008 University of Warwick Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0 Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License. */ /* Laura++ package authors: John Back Paul Harrison Thomas Latham */ /*! \file LauKMatrixPropagator.cc \brief File containing implementation of LauKMatrixPropagator class. */ #include "LauKMatrixPropagator.hh" -#include "LauTextFileParser.hh" -#include "LauKinematics.hh" + #include "LauComplex.hh" +#include "LauJsonTools.hh" +#include "LauKinematics.hh" #include "TMath.h" #include "TSystem.h" #include #include #include #include #include -using std::cout; -using std::endl; -using std::cerr; - ClassImp(LauKMatrixPropagator) +//! \cond DOXYGEN_IGNORE +// map LauKMatrixPropagator::KMatrixChannels values to JSON as strings +NLOHMANN_JSON_SERIALIZE_ENUM( LauKMatrixPropagator::KMatrixChannels, { + {LauKMatrixPropagator::KMatrixChannels::PiPi, "PiPi"}, + {LauKMatrixPropagator::KMatrixChannels::KK, "KK"}, + {LauKMatrixPropagator::KMatrixChannels::FourPi, "FourPi"}, + {LauKMatrixPropagator::KMatrixChannels::EtaEta, "EtaEta"}, + {LauKMatrixPropagator::KMatrixChannels::EtaEtaP, "EtaEtaP"}, + {LauKMatrixPropagator::KMatrixChannels::KPi, "KPi"}, + {LauKMatrixPropagator::KMatrixChannels::KEtaP, "KEtaP"}, + {LauKMatrixPropagator::KMatrixChannels::KThreePi, "KThreePi"}, + {LauKMatrixPropagator::KMatrixChannels::D0K, "D0K"}, + {LauKMatrixPropagator::KMatrixChannels::Dstar0K, "Dstar0K"}, + }) + +std::ostream& operator<<( std::ostream& out, const LauKMatrixPropagator::KMatrixChannels channel ) +{ + switch ( channel ) { + case LauKMatrixPropagator::KMatrixChannels::PiPi : + out << "PiPi"; + break; + case LauKMatrixPropagator::KMatrixChannels::KK : + out << "KK"; + break; + case LauKMatrixPropagator::KMatrixChannels::FourPi : + out << "FourPi"; + break; + case LauKMatrixPropagator::KMatrixChannels::EtaEta : + out << "EtaEta"; + break; + case LauKMatrixPropagator::KMatrixChannels::EtaEtaP : + out << "EtaEtaP"; + break; + case LauKMatrixPropagator::KMatrixChannels::KPi : + out << "KPi"; + break; + case LauKMatrixPropagator::KMatrixChannels::KEtaP : + out << "KEtaP"; + break; + case LauKMatrixPropagator::KMatrixChannels::KThreePi : + out << "KThreePi"; + break; + case LauKMatrixPropagator::KMatrixChannels::D0K : + out << "D0K"; + break; + case LauKMatrixPropagator::KMatrixChannels::Dstar0K : + out << "Dstar0K"; + break; + } + + return out; +} +//! \endcond + LauKMatrixPropagator::LauKMatrixPropagator(const TString& name, const TString& paramFile, - Int_t resPairAmpInt, Int_t nChannels, - Int_t nPoles, Int_t rowIndex) : - name_(name), - paramFileName_(paramFile), - resPairAmpInt_(resPairAmpInt), - index_(rowIndex - 1), - nChannels_(nChannels), - nPoles_(nPoles) + const UInt_t resPairAmpInt, const UInt_t nChannels, + const UInt_t nPoles, const UInt_t rowIndex) : + name_{name}, + paramFileName_{paramFile}, + resPairAmpInt_{resPairAmpInt}, + index_{rowIndex - 1}, + nChannels_{nChannels}, + nPoles_{nPoles} { // Constructor // Check that the index is OK if (index_ < 0 || index_ >= nChannels_) { - std::cerr << "ERROR in LauKMatrixPropagator constructor. The rowIndex, which is set to " - << rowIndex << ", must be between 1 and the number of channels " - << nChannels_ << std::endl; + std::cerr << "ERROR in LauKMatrixPropagator constructor. The rowIndex, which is set to " << rowIndex << ", must be between 1 and the number of channels " << nChannels_ << std::endl; gSystem->Exit(EXIT_FAILURE); } this->setParameters(paramFile); } LauKMatrixPropagator::~LauKMatrixPropagator() { // Destructor realProp_.Clear(); negImagProp_.Clear(); ScattKMatrix_.Clear(); ReRhoMatrix_.Clear(); ImRhoMatrix_.Clear(); GammaMatrix_.Clear(); ReTMatrix_.Clear(); ImTMatrix_.Clear(); IMatrix_.Clear(); zeroMatrix_.Clear(); } -LauComplex LauKMatrixPropagator::getPropTerm(const Int_t channel) const +LauComplex LauKMatrixPropagator::getPropTerm(const UInt_t channel) const { // Get the (i,j) = (index_, channel) term of the propagator // matrix. This allows us not to return the full propagator matrix. Double_t realTerm = this->getRealPropTerm(channel); Double_t imagTerm = this->getImagPropTerm(channel); LauComplex propTerm(realTerm, imagTerm); return propTerm; } -Double_t LauKMatrixPropagator::getRealPropTerm(const Int_t channel) const +Double_t LauKMatrixPropagator::getRealPropTerm(const UInt_t channel) const { // Get the real part of the (i,j) = (index_, channel) term of the propagator // matrix. This allows us not to return the full propagator matrix. if (parametersSet_ == kFALSE) {return 0.0;} Double_t propTerm = realProp_[index_][channel]; return propTerm; } -Double_t LauKMatrixPropagator::getImagPropTerm(const Int_t channel) const +Double_t LauKMatrixPropagator::getImagPropTerm(const UInt_t channel) const { // Get the imaginary part of the (i,j) = (index_, channel) term of the propagator // matrix. This allows us not to return the full propagator matrix. if (parametersSet_ == kFALSE) {return 0.0;} Double_t imagTerm = -1.0*negImagProp_[index_][channel]; return imagTerm; } void LauKMatrixPropagator::updatePropagator(const Double_t s) { // Calculate the K-matrix propagator for the given s value. // The K-matrix amplitude is given by // T_i = sum_{ij} (I - iK*rho)^-1 * P_j, where P is the production K-matrix. // i = index for the state (e.g. S-wave index = 0). // Here, we only find the (I - iK*rho)^-1 matrix part. // Check if we have almost the same s value as before. If so, don't re-calculate // the propagator nor any of the pole mass summation terms. if (TMath::Abs(s - previousS_) < 1e-6*s) { - //cout<<"Already got propagator for s = "<calcPoleDenomVect(s); this->updateAdlerZeroFactor(s); // Calculate the scattering K-matrix (real and symmetric) this->calcScattKMatrix(s); // Calculate the phase space density matrix, which is diagonal, but can be complex // if the quantity s is below various threshold values (analytic continuation). this->calcRhoMatrix(s); // Calculate the angular momentum barrier matrix, which is real and diagonal this->calcGammaMatrix(s); // Calculate K*rho*(gamma^2) (real and imaginary parts, since rho can be complex) TMatrixD GammaMatrixSq = (GammaMatrix_*GammaMatrix_); TMatrixD K_realRhoGammaSq(ScattKMatrix_); K_realRhoGammaSq *= ReRhoMatrix_; K_realRhoGammaSq *= GammaMatrixSq; TMatrixD K_imagRhoGammaSq(ScattKMatrix_); K_imagRhoGammaSq *= ImRhoMatrix_; K_imagRhoGammaSq *= GammaMatrixSq; // A = I + K*Imag(rho)Gamma^2, B = -K*Real(Rho)Gamma^2 // Calculate C and D matrices such that (A + iB)*(C + iD) = I, // ie. C + iD = (I - i K*rhoGamma^2)^-1, separated into real and imaginary parts. // realProp C = (A + B A^-1 B)^-1, imagProp D = -A^-1 B C TMatrixD A(IMatrix_); A += K_imagRhoGammaSq; TMatrixD B(zeroMatrix_); B -= K_realRhoGammaSq; TMatrixD invA(TMatrixD::kInverted, A); TMatrixD invA_B(invA); invA_B *= B; TMatrixD B_invA_B(B); B_invA_B *= invA_B; TMatrixD invC(A); invC += B_invA_B; // Set the real part of the propagator matrix ("C") realProp_ = TMatrixD(TMatrixD::kInverted, invC); // Set the (negative) imaginary part of the propagator matrix ("-D") TMatrixD BC(B); BC *= realProp_; negImagProp_ = TMatrixD(invA); negImagProp_ *= BC; // Pre-multiply by the Gamma matrix: realProp_ = GammaMatrix_ * realProp_; negImagProp_ = GammaMatrix_ * negImagProp_; if(verbose_) { std::cout << "In LauKMatrixPropagator::updatePropagator(s). D[1-iKrhoD^2]^-1: " << std::endl; TString realOutput("Real part:"), imagOutput("Imag part:"); - for (int iChannel = 0; iChannel < nChannels_; iChannel++) + for (UInt_t iChannel = 0; iChannel < nChannels_; iChannel++) { - for (int jChannel = 0; jChannel < nChannels_; jChannel++) + for (UInt_t jChannel = 0; jChannel < nChannels_; jChannel++) { realOutput += Form("\t%.6f",realProp_[iChannel][jChannel]); imagOutput += Form("\t%.6f",-1*negImagProp_[iChannel][jChannel]); } realOutput += "\n "; imagOutput += "\n "; } std::cout << realOutput << std::endl; std::cout << imagOutput << std::endl; } // Also calculate the production SVP term, since this uses Adler-zero parameters // defined in the parameter file. this->updateProdSVPTerm(s); // Finally, keep track of the value of s we just used. previousS_ = s; } void LauKMatrixPropagator::setParameters(const TString& inputFile) { - // Read an input file that specifies the values of the coupling constants g_i for // the given number of poles and their (bare) masses. Also provided are the f_{ab} // slow-varying constants. The input file should also provide the Adler zero // constants s_0, s_A and s_A0. parametersSet_ = kFALSE; - cout<<"Initialising K-matrix propagator "<Exit(EXIT_FAILURE); + if ( j.is_null() ) { + std::cerr << "ERROR in LauKMatrixPropagator::setParameters : unable to retrieve JSON object from root element of file \"" << inputFile << "\"" << std::endl; + return; } - UInt_t iLine(0); - - for (iLine = 1; iLine <= nTotLines; iLine++) { - - // Get the line of strings - std::vector theLine = readFile.getLine(iLine); - - // There should always be at least two strings: a keyword and at least 1 value - if (theLine.size() < 2) {continue;} - - TString keyword(theLine[0].c_str()); - keyword.ToLower(); // Use lowercase - - if (!keyword.CompareTo("channels")) { - - // Channel indices for phase-space factors - this->storeChannels(theLine); - - } else if (!keyword.CompareTo("pole")) { - - // Pole terms - this->storePole(theLine); - - } else if (!keyword.CompareTo("scatt")) { - - // Scattering terms - this->storeScattering(theLine); - - } else if (!keyword.CompareTo("angularmomentum")) { - - // Orbital angular momentum state for each channel & set default a values if called before storeBarrierFactorParameter - this->storeOrbitalAngularMomenta(theLine, a); - - } else if (!keyword.CompareTo("barrierfactorparameter")) { + const std::vector mandatoryElements { + std::make_pair("channels", JsonType::Array), + std::make_pair("poles", JsonType::Array), + std::make_pair("scattering", JsonType::Array), + std::make_pair("parameters", JsonType::Array) + }; + if ( ! LauJsonTools::checkObjectElements( j, mandatoryElements ) ) { + std::cerr << "ERROR in LauKMatrixPropagator::setParameters : aborting processing due to mis-formatted elements" << std::endl; + return; + } - // Value of parameter "a" in denominator of centrifugal barrier factor, gamma - this->storeBarrierFactorParameter(theLine, a); + // Channel indices for phase-space factors + if ( ! this->storeChannels( j.at("channels") ) ) { + return; + } - } else if (!keyword.CompareTo("radii")) { + // Pole terms + if ( ! this->storePoles( j.at("poles") ) ) { + return; + } - // Values of characteristic radius - this->storeRadii(theLine); + // Scattering terms + if ( ! this->storeScattering( j.at("scattering") ) ) { + return; + } - } else { + // Orbital angular momentum state for each channel & set default a values if called before storeBarrierFactorParameter + if ( j.contains("angularmomentum") ) { + if ( ! this->storeOrbitalAngularMomenta( j.at("angularmomentum"), a ) ) { + return; + } + } - // Usually Adler-zero constants - TString parString(theLine[1].c_str()); - this->storeParameter(keyword, parString); + // Value of parameter "a" in denominator of centrifugal barrier factor, gamma + if ( j.contains("barrierfactorparameter") ) { + if ( ! this->storeBarrierFactorParameter( j.at("barrierfactorparameter"), a ) ) { + return; + } + } + // Values of characteristic radius + if ( j.contains("radii") ) { + if ( ! this->storeRadii( j.at("radii") ) ) { + return; } + } + // Other miscellaneous parameters, e.g. Adler-zero constants + for ( auto& parameter : j.at("parameters") ) { + if ( ! this->storeParameter( parameter ) ) { + return; + } } sAConst_ = 0.5*sA_.unblindValue()*mPiSq_; // Symmetrise scattering parameters if enabled - if (scattSymmetry_ == kTRUE) { - - for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) { - - for (Int_t jChannel = iChannel; jChannel < nChannels_; jChannel++) { - - LauParameter fPar = fScattering_[iChannel][jChannel]; - fScattering_[jChannel][iChannel] = LauParameter(fPar); - + if ( scattSymmetry_ ) { + for (UInt_t iChannel{0}; iChannel < nChannels_; ++iChannel) { + for (UInt_t jChannel{iChannel}; jChannel < nChannels_; ++jChannel) { + fScattering_[jChannel][iChannel] = fScattering_[iChannel][jChannel]; } } } // Now that radii and barrier-factor-denominator parameters have been set, cache the value of "a/(R*R)" - for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) { + for (UInt_t iChannel{0}; iChannel < nChannels_; ++iChannel) { gamAInvRadSq_[iChannel] = a[iChannel]/(radii_[iChannel]*radii_[iChannel]); } // All required parameters have been set parametersSet_ = kTRUE; - cout<<"Finished initialising K-matrix propagator "< >. // Set their sizes using the number of poles and channels defined in the constructor gCouplings_.clear(); gCouplings_.resize(nPoles_); - for (Int_t iPole = 0; iPole < nPoles_; iPole++) { + for (UInt_t iPole = 0; iPole < nPoles_; iPole++) { gCouplings_[iPole].resize(nChannels_); } fScattering_.clear(); fScattering_.resize(nChannels_); - for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) { + for (UInt_t iChannel = 0; iChannel < nChannels_; iChannel++) { fScattering_[iChannel].resize(nChannels_); } // Clear other vectors phaseSpaceTypes_.clear(); phaseSpaceTypes_.resize(nChannels_); mSqPoles_.clear(); mSqPoles_.resize(nPoles_); - haveCalled_storeBarrierFactorParameter = kFALSE; + haveCalledStoreBarrierFactorParameter_ = kFALSE; } -void LauKMatrixPropagator::storeChannels(const std::vector& theLine) +bool LauKMatrixPropagator::storeChannels(const nlohmann::json& j) { - // Get the list of channel indices to specify what phase space factors should be used // e.g. pipi, Kpi, eta eta', 4pi etc.. - // Check that the line has nChannels_+1 strings - Int_t nTypes = static_cast(theLine.size()) - 1; - if (nTypes != nChannels_) { - cerr<<"Error in LauKMatrixPropagator::storeChannels. The input file defines " - <checkPhaseSpaceType(phaseSpaceInt); - - if (checkChannel == kTRUE) { - cout<<"Adding phase space channel index "<(phaseSpaceInt); - } else { - cerr<<"Phase space channel index "<(LauKMatrixPropagator::KMatrixChannels::TotChannels)-1<(); + for ( UInt_t iChannel{0}; iChannel < nChannels_; ++iChannel) { + std::cout << "INFO in LauKMatrixPropagator::storeChannels : added phase space channel " << phaseSpaceTypes_[iChannel] << " with index " << iChannel+1 << " to propagator " << name_ << std::endl; } + return true; } -void LauKMatrixPropagator::storePole(const std::vector& theLine) +bool LauKMatrixPropagator::storePoles(const nlohmann::json& j) { + // For each pole, store the pole mass and its coupling constants for each channel - // Store the pole mass and its coupling constants for each channel. - // Each line will contain: Pole poleNumber poleMass poleCouplingsPerChannel - - // Check that the line has nChannels_ + 3 strings - Int_t nWords = static_cast(theLine.size()); - Int_t nExpect = nChannels_ + 3; - - if (nWords == nExpect) { - - Int_t poleIndex = std::atoi(theLine[1].c_str()) - 1; - if (poleIndex >= 0 && poleIndex < nPoles_) { - - Double_t poleMass = std::atof(theLine[2].c_str()); - Double_t poleMassSq = poleMass*poleMass; - LauParameter mPoleParam(Form("KM_%s_poleMassSq_%i",name_.Data(),poleIndex),poleMassSq); - mSqPoles_[poleIndex] = mPoleParam; - - cout<<"Added bare pole mass "< mandatoryElements { + std::make_pair("index", JsonType::Number_Unsigned), + std::make_pair("mass", JsonType::Number_Float), + std::make_pair("couplings", JsonType::Array) + }; + + for ( auto& elem : j ) { + if ( ! checkObjectElements( elem, mandatoryElements ) ) { + std::cerr << "ERROR in LauKMatrixPropagator::storePoles : aborting processing due to mis-formatted elements" << std::endl; + return false; + } + } - cout<<"Added coupling parameter g^{"<= nPoles_ ) { + std::cerr << "ERROR in LauKMatrixPropagator::storePoles : invalid pole index: " << poleIndex+1 << ", aborting processing" << std::endl; + return false; + } - } + const auto poleMass { LauJsonTools::getValue( poleDef, "mass" ) }; + const auto poleMassSq { poleMass * poleMass }; + const auto couplings { LauJsonTools::getValue>( poleDef, "couplings" ) }; + if ( couplings.size() != nChannels_ ) { + std::cerr << "ERROR in LauKMatrixPropagator::storePoles : invalid number of couplings ( " << couplings.size() << " ) provided for pole index " << poleIndex+1 << ", aborting processing" << std::endl; + return false; } - } else { + mSqPoles_[poleIndex] = LauParameter{ Form("KM_%s_poleMassSq_%i", name_.Data(), poleIndex+1), poleMassSq }; - cerr<<"Error in LauKMatrixPropagator::storePole. Expecting "< mandatoryElements { + std::make_pair("index", JsonType::Number_Unsigned), + std::make_pair("couplings", JsonType::Array) + }; + + for ( auto& elem : j ) { + if ( ! checkObjectElements( elem, mandatoryElements ) ) { + std::cerr << "ERROR in LauKMatrixPropagator::storeScattering : aborting processing due to mis-formatted elements" << std::endl; + return false; } + } - } else { + for ( auto& scattDef : j ) { + const auto scattIndex { LauJsonTools::getValue( scattDef, "index" ) - 1 }; + if ( scattIndex >= nChannels_ ) { + std::cerr << "ERROR in LauKMatrixPropagator::storeScattering : invalid channel index: " << scattIndex+1 << ", aborting processing" << std::endl; + return false; + } - cerr<<"Error in LauKMatrixPropagator::storeScattering. Expecting "<>( scattDef, "couplings" ) }; + if ( couplings.size() != nChannels_ ) { + std::cerr << "ERROR in LauKMatrixPropagator::storeScattering : invalid number of couplings ( " << couplings.size() << " ) provided for scattering index " << scattIndex+1 << ", aborting processing" << std::endl; + return false; + } + for ( UInt_t iChannel{0}; iChannel < nChannels_; ++iChannel) { + const Double_t couplingConst { couplings[iChannel] }; + fScattering_[scattIndex][iChannel] = LauParameter{ Form("KM_%s_fScatteringConst_%i_%i", name_.Data(), scattIndex+1, iChannel+1), couplingConst }; + std::cout << "INFO in LauKMatrixPropagator::storeScattering : Added coupling parameter f(" << scattIndex+1 << "," << iChannel+1 << ") = " << couplingConst << std::endl; + } } + return true; } -void LauKMatrixPropagator::storeOrbitalAngularMomenta(const std::vector& theLine, std::vector& a) +bool LauKMatrixPropagator::storeOrbitalAngularMomenta(const nlohmann::json& j, std::vector& a) { - // Store the orbital angular momentum for each channel - // Each line will contain: angularmomentum OrbitalAngularMomentumPerChannel - - // Check that the line has nChannels_ + 1 strings - Int_t nWords = static_cast(theLine.size()); - Int_t nExpect = nChannels_ + 1; - - if (nWords == nExpect) { - for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) { - - Int_t angularMomentum = std::atoi(theLine[iChannel+1].c_str()); - L_[iChannel] = angularMomentum; + if ( ! LauJsonTools::checkValueType( j, LauJsonTools::JsonType::Array ) ) { + std::cerr << "ERROR in LauKMatrixPropagator::storeOrbitalAngularMomenta : \"angularmomentum\" element is not an array, aborting processing" << std::endl; + return false; + } - cout<<"Defined K-matrix orbital angular momentum "<(); - } else { + for ( UInt_t iChannel{0}; iChannel < nChannels_; ++iChannel) { + std::cout << "INFO in LauKMatrixPropagator::storeOrbitalAngularMomenta : defined K-matrix orbital angular momentum " << L_[iChannel] << " for channel " << iChannel+1 << std::endl; + } - cerr<<"Error in LauKMatrixPropagator::storeOrbitalAngularMomenta. Expecting "<Exit(EXIT_FAILURE); - } + // Otherwise, set default value of spin-dependent centrifugal-barrier-factor parameter + for ( UInt_t iCh{0}; iCh < nChannels_; ++iCh) { + switch ( L_[iCh] ) { + case 0: + a[iCh] = 0; + break; + case 1: + a[iCh] = 1; + break; + case 2: + a[iCh] = 3; + break; + default: + std::cerr << "ERROR in LauKMatrixPropagator::storeOrbitalAngularMomenta : barrier factor and angular-momentum terms of K-matrix are only defined for S-, P-, or D-wave." << std::endl; + return false; } } + return true; } -void LauKMatrixPropagator::storeRadii(const std::vector& theLine) +bool LauKMatrixPropagator::storeRadii(const nlohmann::json& j) { + // Store the characteristic barrier radii (measured in GeV^{-1}) - // Store the characteristic radii (measured in GeV^{-1}) - // Each line will contain: Radii RadiusConstantsPerChannel - - // Check that the line has nChannels_ + 1 strings - Int_t nWords = static_cast(theLine.size()); - Int_t nExpect = nChannels_ + 1; - - if (nWords == nExpect) { - - for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) { - - Double_t radiusConst = std::atof(theLine[iChannel+1].c_str()); - radii_[iChannel] = radiusConst; - - cout<<"Added K-matrix radius "<(); + for ( UInt_t iChannel{0}; iChannel < nChannels_; ++iChannel) { + std::cout << "INFO in LauKMatrixPropagator::storeRadii : added K-matrix barrier radius " << radii_[iChannel] << " for channel " << iChannel+1 << std::endl; } + return true; } -void LauKMatrixPropagator::storeBarrierFactorParameter(const std::vector& theLine, std::vector& a) +bool LauKMatrixPropagator::storeBarrierFactorParameter(const nlohmann::json& j, std::vector& a) { - // Store the parameter of the barrier factor - // Each line will contain: barrierfactorparameter ParameterValuePerchannel - - // Check that the line has nChannels_ + 1 strings - Int_t nWords = static_cast(theLine.size()); - Int_t nExpect = nChannels_ + 1; - if (nWords == nExpect) { - - for (Int_t iChannel = 0; iChannel < nChannels_; iChannel++) { - - Double_t parameterValue = std::atof(theLine[iChannel+1].c_str()); - a[iChannel] = parameterValue; - - cout<<"Added K-matrix barrier factor parameter value "<>(); - } else { + for ( UInt_t iChannel{0}; iChannel < nChannels_; ++iChannel) { + std::cout << "INFO in LauKMatrixPropagator::storeBarrierFactorParameter : added K-matrix barrier factor parameter " << a[iChannel] << " for channel " << iChannel+1 << std::endl; + } - cerr<<"Error in LauKMatrixPropagator::storeBarrierFactorParameter. Expecting "< mandatoryElements { + std::make_pair("name", JsonType::String), + std::make_pair("value", JsonType::Any) + }; + if ( ! checkObjectElements( j, mandatoryElements ) ) { + std::cerr << "ERROR in LauKMatrixPropagator::storeParameter : aborting processing due to mis-formatted elements" << std::endl; + return false; + } + + auto keyword { getValue( j, "name" ) }; + keyword.ToLower(); - if (!keyword.CompareTo("msq0")) { + if ( keyword == "msq0" ) { - Double_t mSq0Value = std::atof(parString.Data()); - cout<<"Adler zero constant m0Sq = "<( j, "value" ) }; + mSq0_ = LauParameter{ Form("KM_%s_mSq0",name_.Data()), mSq0Value }; + std::cout << "INFO in LauKMatrixPropagator::storeParameter : Adler zero constant m0Sq = " << mSq0Value << std::endl; - } else if (!keyword.CompareTo("s0scatt")) { + } else if ( keyword == "s0scatt" ) { - Double_t s0ScattValue = std::atof(parString.Data()); - cout<<"Adler zero constant s0Scatt = "<( j, "value" ) }; + s0Scatt_ = LauParameter{ Form("KM_%s_s0Scatt",name_.Data()), s0ScattValue }; + std::cout << "INFO in LauKMatrixPropagator::storeParameter : Adler zero constant s0Scatt = " << s0ScattValue << std::endl; - } else if (!keyword.CompareTo("s0prod")) { + } else if ( keyword == "s0prod" ) { - Double_t s0ProdValue = std::atof(parString.Data()); - cout<<"Adler zero constant s0Prod = "<( j, "value" ) }; + s0Prod_ = LauParameter{ Form("KM_%s_s0Prod",name_.Data()), s0ProdValue }; + std::cout << "INFO in LauKMatrixPropagator::storeParameter : Adler zero constant s0Prod = " << s0ProdValue << std::endl; - } else if (!keyword.CompareTo("sa0")) { + } else if ( keyword == "sa0" ) { - Double_t sA0Value = std::atof(parString.Data()); - cout<<"Adler zero constant sA0 = "<( j, "value" ) }; + sA0_ = LauParameter{ Form("KM_%s_sA0",name_.Data()), sA0Value }; + std::cout << "INFO in LauKMatrixPropagator::storeParameter : Adler zero constant sA0 = " << sA0Value << std::endl; - } else if (!keyword.CompareTo("sa")) { + } else if ( keyword == "sa" ) { - Double_t sAValue = std::atof(parString.Data()); - cout<<"Adler zero constant sA = "<( j, "value" ) }; + sA_ = LauParameter{ Form("KM_%s_sA",name_.Data()), sAValue }; + std::cout << "INFO in LauKMatrixPropagator::storeParameter : Adler zero constant sA = " << sAValue << std::endl; - } else if (!keyword.CompareTo("scattsymmetry")) { + } else if ( keyword == "scattsymmetry" ) { - Int_t flag = std::atoi(parString.Data()); - if (flag == 0) { - cout<<"Turning off scattering parameter symmetry: f_ji = f_ij will not be assumed"<( j, "value" ); + if ( ! scattSymmetry_ ) { + std::cout << "INFO in LauKMatrixPropagator::storeParameter : Turning off scattering parameter symmetry: f_ji = f_ij will not be assumed" << std::endl; } } + return true; } void LauKMatrixPropagator::calcScattKMatrix(const Double_t s) { // Calculate the scattering K-matrix for the given value of s. // We need to obtain the complete matrix (not just one row/column) // to get the correct inverted (I - i K rho) terms for the propagator. - if (verbose_) {cout<<"Within calcScattKMatrix for s = "<updateScattSVPTerm(s); // Now loop over iChannel, jChannel to calculate Kij = Kji. for (iChannel = 0; iChannel < nChannels_; iChannel++) { // Scattering matrix is real and symmetric. Start j loop from i. for (jChannel = iChannel; jChannel < nChannels_; jChannel++) { Double_t Kij(0.0); // Calculate pole mass summation term for (iPole = 0; iPole < nPoles_; iPole++) { Double_t g_i = this->getCouplingConstant(iPole, iChannel); Double_t g_j = this->getCouplingConstant(iPole, jChannel); Kij += poleDenomVect_[iPole]*g_i*g_j; - if (verbose_) {cout<<"1: Kij for i = "<getScatteringConstant(iChannel, jChannel); Kij += fij*scattSVP_; Kij *= adlerZeroFactor_; - if (verbose_) {cout<<"2: Kij for i = "< 1.0e-6) {invPoleTerm = 1.0/poleTerm;} poleDenomVect_.push_back(invPoleTerm); } } -Double_t LauKMatrixPropagator::getPoleDenomTerm(const Int_t poleIndex) const +Double_t LauKMatrixPropagator::getPoleDenomTerm(const UInt_t poleIndex) const { if (parametersSet_ == kFALSE) {return 0.0;} Double_t poleDenom(0.0); poleDenom = poleDenomVect_[poleIndex]; return poleDenom; } -LauParameter& LauKMatrixPropagator::getPoleMassSqParameter(const Int_t poleIndex) +LauParameter& LauKMatrixPropagator::getPoleMassSqParameter(const UInt_t poleIndex) { if ( (parametersSet_ == kFALSE) || (poleIndex < 0 || poleIndex >= nPoles_) ) { std::cerr << "ERROR from LauKMatrixPropagator::getPoleMassSqParameter(). Invalid pole." << std::endl; gSystem->Exit(EXIT_FAILURE); } return mSqPoles_[poleIndex]; } -Double_t LauKMatrixPropagator::getCouplingConstant(const Int_t poleIndex, const Int_t channelIndex) const +Double_t LauKMatrixPropagator::getCouplingConstant(const UInt_t poleIndex, const UInt_t channelIndex) const { if (parametersSet_ == kFALSE) {return 0.0;} if (poleIndex < 0 || poleIndex >= nPoles_) {return 0.0;} if (channelIndex < 0 || channelIndex >= nChannels_) {return 0.0;} Double_t couplingConst = gCouplings_[poleIndex][channelIndex].unblindValue(); return couplingConst; } -LauParameter& LauKMatrixPropagator::getCouplingParameter(const Int_t poleIndex, const Int_t channelIndex) +LauParameter& LauKMatrixPropagator::getCouplingParameter(const UInt_t poleIndex, const UInt_t channelIndex) { if ( (parametersSet_ == kFALSE) || (poleIndex < 0 || poleIndex >= nPoles_) || (channelIndex < 0 || channelIndex >= nChannels_) ) { std::cerr << "ERROR from LauKMatrixPropagator::getCouplingParameter(). Invalid coupling." << std::endl; gSystem->Exit(EXIT_FAILURE); } //std::cout << "Minvalue + range for " << poleIndex << ", " << channelIndex << ": " << gCouplings_[poleIndex][channelIndex].minValue() << " => + " << gCouplings_[poleIndex][channelIndex].range() << // " and init value: " << gCouplings_[poleIndex][channelIndex].initValue() << std::endl; return gCouplings_[poleIndex][channelIndex]; } -Double_t LauKMatrixPropagator::getScatteringConstant(const Int_t channel1Index, const Int_t channel2Index) const +Double_t LauKMatrixPropagator::getScatteringConstant(const UInt_t channel1Index, const UInt_t channel2Index) const { if (parametersSet_ == kFALSE) {return 0.0;} if (channel1Index < 0 || channel1Index >= nChannels_) {return 0.0;} if (channel2Index < 0 || channel2Index >= nChannels_) {return 0.0;} Double_t scatteringConst = fScattering_[channel1Index][channel2Index].unblindValue(); return scatteringConst; } -LauParameter& LauKMatrixPropagator::getScatteringParameter(const Int_t channel1Index, const Int_t channel2Index) +LauParameter& LauKMatrixPropagator::getScatteringParameter(const UInt_t channel1Index, const UInt_t channel2Index) { if ( (parametersSet_ == kFALSE) || (channel1Index < 0 || channel1Index >= nChannels_) || (channel2Index < 0 || channel2Index >= nChannels_) ) { std::cerr << "ERROR from LauKMatrixPropagator::getScatteringParameter(). Invalid chanel index." << std::endl; gSystem->Exit(EXIT_FAILURE); } return fScattering_[channel1Index][channel2Index]; } Double_t LauKMatrixPropagator::calcSVPTerm(const Double_t s, const Double_t s0) const { if (parametersSet_ == kFALSE) {return 0.0;} // Calculate the "slowly-varying part" (SVP) Double_t result(0.0); Double_t deltaS = s - s0; if (TMath::Abs(deltaS) > 1.0e-6) { result = (mSq0_.unblindValue() - s0)/deltaS; } return result; } void LauKMatrixPropagator::updateScattSVPTerm(const Double_t s) { // Update the scattering "slowly-varying part" (SVP) Double_t s0Scatt = s0Scatt_.unblindValue(); scattSVP_ = this->calcSVPTerm(s, s0Scatt); } void LauKMatrixPropagator::updateProdSVPTerm(const Double_t s) { // Update the production "slowly-varying part" (SVP) Double_t s0Prod = s0Prod_.unblindValue(); prodSVP_ = this->calcSVPTerm(s, s0Prod); } void LauKMatrixPropagator::updateAdlerZeroFactor(const Double_t s) { // Calculate the multiplicative factor containing various Adler zero // constants. adlerZeroFactor_ = 0.0; Double_t sA0Val = sA0_.unblindValue(); Double_t deltaS = s - sA0Val; if (TMath::Abs(deltaS) > 1e-6) { adlerZeroFactor_ = (s - sAConst_)*(1.0 - sA0Val)/deltaS; } } void LauKMatrixPropagator::calcGammaMatrix(const Double_t s) { // Calculate the gamma angular momentum barrier matrix // for the given invariant mass squared quantity, s. // Initialise all entries to zero GammaMatrix_.Zero(); Double_t gamma(0.0); - for (Int_t iChannel (0); iChannel < nChannels_; ++iChannel) { + for (UInt_t iChannel (0); iChannel < nChannels_; ++iChannel) { if ( L_[iChannel] != 0 ) { gamma = this->calcGamma(iChannel,s); } else { gamma = 1.0; // S-wave } if (verbose_) { - cout<<"GammaMatrix("<Exit(EXIT_FAILURE); } return rho; } LauComplex LauKMatrixPropagator::calcD0KRho(const Double_t s) const { // Calculate the D0K+ phase space factor LauComplex rho(0.0, 0.0); if (TMath::Abs(s) < 1e-10) {return rho;} Double_t sqrtTerm1 = (-mD0KSumSq_/s) + 1.0; Double_t sqrtTerm2 = (-mD0KDiffSq_/s) + 1.0; Double_t sqrtTerm = sqrtTerm1*sqrtTerm2; if (sqrtTerm < 0.0) { rho.setImagPart( TMath::Sqrt(-sqrtTerm) ); } else { rho.setRealPart( TMath::Sqrt(sqrtTerm) ); } return rho; } LauComplex LauKMatrixPropagator::calcDstar0KRho(const Double_t s) const { // Calculate the Dstar0K+ phase space factor LauComplex rho(0.0, 0.0); if (TMath::Abs(s) < 1e-10) {return rho;} Double_t sqrtTerm1 = (-mDstar0KSumSq_/s) + 1.0; Double_t sqrtTerm2 = (-mDstar0KDiffSq_/s) + 1.0; Double_t sqrtTerm = sqrtTerm1*sqrtTerm2; if (sqrtTerm < 0.0) { rho.setImagPart( TMath::Sqrt(-sqrtTerm) ); } else { rho.setRealPart( TMath::Sqrt(sqrtTerm) ); } return rho; } LauComplex LauKMatrixPropagator::calcPiPiRho(const Double_t s) const { // Calculate the pipi phase space factor LauComplex rho(0.0, 0.0); if (TMath::Abs(s) < 1e-10) {return rho;} Double_t sqrtTerm = (-m2piSq_/s) + 1.0; if (sqrtTerm < 0.0) { rho.setImagPart( TMath::Sqrt(-sqrtTerm) ); } else { rho.setRealPart( TMath::Sqrt(sqrtTerm) ); } return rho; } LauComplex LauKMatrixPropagator::calcKKRho(const Double_t s) const { // Calculate the KK phase space factor LauComplex rho(0.0, 0.0); if (TMath::Abs(s) < 1e-10) {return rho;} Double_t sqrtTerm = (-m2KSq_/s) + 1.0; if (sqrtTerm < 0.0) { rho.setImagPart( TMath::Sqrt(-sqrtTerm) ); } else { rho.setRealPart( TMath::Sqrt(sqrtTerm) ); } return rho; } LauComplex LauKMatrixPropagator::calcFourPiRho(const Double_t s) const { // Calculate the 4pi phase space factor. This uses a 6th-order polynomial // parameterisation that approximates the multi-body phase space double integral // defined in Eq 4 of the A&S paper hep-ph/0204328. This form agrees with the // BaBar model (another 6th order polynomial from s^4 down to 1/s^2), but avoids the // exponential increase at small values of s (~< 0.1) arising from 1/s and 1/s^2. // Eq 4 is evaluated for each value of s by assuming incremental steps of 1e-3 GeV^2 // for s1 and s2, the invariant energy squared of each of the di-pion states, // with the integration limits of s1 = (2*mpi)^2 to (sqrt(s) - 2*mpi)^2 and // s2 = (2*mpi)^2 to (sqrt(s) - sqrt(s1))^2. The mass M of the rho is taken to be // 0.775 GeV and the energy-dependent width of the 4pi system // Gamma(s) = gamma_0*rho1^3(s), where rho1 = sqrt(1.0 - 4*mpiSq/s) and gamma_0 is // the "width" of the 4pi state at s = 1, which is taken to be 0.3 GeV // (~75% of the total width from PDG estimates of the f0(1370) -> 4pi state). // The normalisation term rho_0 is found by ensuring that the phase space integral // at s = 1 is equal to sqrt(1.0 - 16*mpiSq/s). Note that the exponent for this // factor in hep-ph/0204328 is wrong; it should be 0.5, i.e. sqrt, not n = 1 to 5. // Plotting the value of this double integral as a function of s can then be fitted // to a 6th-order polynomial (for s < 1), which is the result used below LauComplex rho(0.0, 0.0); if (TMath::Abs(s) < 1e-10) {return rho;} if (s <= 1.0) { Double_t rhoTerm = ((1.07885*s + 0.13655)*s - 0.29744)*s - 0.20840; rhoTerm = ((rhoTerm*s + 0.13851)*s - 0.01933)*s + 0.00051; // For some values of s (below 2*mpi), this term is a very small // negative number. Check for this and set the rho term to zero. if (rhoTerm < 0.0) {rhoTerm = 0.0;} rho.setRealPart( rhoTerm ); } else { rho.setRealPart( TMath::Sqrt(1.0 - (fourPiFactor1_/s)) ); } return rho; } LauComplex LauKMatrixPropagator::calcEtaEtaRho(const Double_t s) const { // Calculate the eta-eta phase space factor LauComplex rho(0.0, 0.0); if (TMath::Abs(s) < 1e-10) {return rho;} Double_t sqrtTerm = (-m2EtaSq_/s) + 1.0; if (sqrtTerm < 0.0) { rho.setImagPart( TMath::Sqrt(-sqrtTerm) ); } else { rho.setRealPart( TMath::Sqrt(sqrtTerm) ); } return rho; } LauComplex LauKMatrixPropagator::calcEtaEtaPRho(const Double_t s) const { // Calculate the eta-eta' phase space factor. Note that the // mass difference term m_eta - m_eta' is not included, // since this corresponds to a "t or u-channel crossing", // which means that we cannot simply analytically continue // this part of the phase space factor below threshold, which // we can do for s-channel contributions. This is actually an // unsolved problem, e.g. see Guo et al 1409.8652, and // Danilkin et al 1409.7708. Anisovich and Sarantsev in // hep-ph/0204328 "solve" this issue by setting the mass // difference term to unity, which is what we do here... LauComplex rho(0.0, 0.0); if (TMath::Abs(s) < 1e-10) {return rho;} Double_t sqrtTerm = (-mEtaEtaPSumSq_/s) + 1.0; if (sqrtTerm < 0.0) { rho.setImagPart( TMath::Sqrt(-sqrtTerm) ); } else { rho.setRealPart( TMath::Sqrt(sqrtTerm) ); } return rho; } LauComplex LauKMatrixPropagator::calcKPiRho(const Double_t s) const { // Calculate the K-pi phase space factor LauComplex rho(0.0, 0.0); if (TMath::Abs(s) < 1e-10) {return rho;} Double_t sqrtTerm1 = (-mKpiSumSq_/s) + 1.0; Double_t sqrtTerm2 = (-mKpiDiffSq_/s) + 1.0; Double_t sqrtTerm = sqrtTerm1*sqrtTerm2; if (sqrtTerm < 0.0) { rho.setImagPart( TMath::Sqrt(-sqrtTerm) ); } else { rho.setRealPart( TMath::Sqrt(sqrtTerm) ); } return rho; } LauComplex LauKMatrixPropagator::calcKEtaPRho(const Double_t s) const { // Calculate the K-eta' phase space factor LauComplex rho(0.0, 0.0); if (TMath::Abs(s) < 1e-10) {return rho;} Double_t sqrtTerm1 = (-mKEtaPSumSq_/s) + 1.0; Double_t sqrtTerm2 = (-mKEtaPDiffSq_/s) + 1.0; Double_t sqrtTerm = sqrtTerm1*sqrtTerm2; if (sqrtTerm < 0.0) { rho.setImagPart( TMath::Sqrt(-sqrtTerm) ); } else { rho.setRealPart( TMath::Sqrt(sqrtTerm) ); } return rho; } LauComplex LauKMatrixPropagator::calcKThreePiRho(const Double_t s) const { // Calculate the Kpipipi + multimeson phase space factor. // Use the simplest definition in hep-ph/9705401 (Eq 14), which is the form // used for the rest of that paper (thankfully, the amplitude does not depend // significantly on the form used for the K3pi phase space factor). LauComplex rho(0.0, 0.0); if (TMath::Abs(s) < 1e-10) {return rho;} if (s < 1.44) { Double_t powerTerm = (-mK3piDiffSq_/s) + 1.0; if (powerTerm < 0.0) { rho.setImagPart( k3piFactor_*TMath::Power(-powerTerm, 2.5) ); } else { rho.setRealPart( k3piFactor_*TMath::Power(powerTerm, 2.5) ); } } else { rho.setRealPart( 1.0 ); } return rho; } -Bool_t LauKMatrixPropagator::checkPhaseSpaceType(const Int_t phaseSpaceInt) const -{ - Bool_t passed(kFALSE); - - if (phaseSpaceInt >= 1 && phaseSpaceInt < static_cast(LauKMatrixPropagator::KMatrixChannels::TotChannels)) { - passed = kTRUE; - } - - return passed; -} - -LauComplex LauKMatrixPropagator::getTransitionAmp(const Double_t s, const Int_t channel) +LauComplex LauKMatrixPropagator::getTransitionAmp(const Double_t s, const UInt_t channel) { // Get the complex (unitary) transition amplitude T for the given channel LauComplex TAmp(0.0, 0.0); if (channel <= 0 || channel > nChannels_) {return TAmp;} this->getTMatrix(s); TAmp.setRealPart(ReTMatrix_[index_][channel-1]); TAmp.setImagPart(ImTMatrix_[index_][channel-1]); return TAmp; } -LauComplex LauKMatrixPropagator::getPhaseSpaceTerm(const Double_t s, const Int_t channel) +LauComplex LauKMatrixPropagator::getPhaseSpaceTerm(const Double_t s, const UInt_t channel) { // Get the complex (unitary) transition amplitude T for the given channel LauComplex rho(0.0, 0.0); if (channel <= 0 || channel > nChannels_) {return rho;} // If s has changed from the previous value, recalculate rho if (TMath::Abs(s - previousS_) > 1e-6*s) { this->calcRhoMatrix(s); } rho.setRealPart(ReRhoMatrix_[channel][channel-1]); rho.setImagPart(ImRhoMatrix_[channel][channel-1]); return rho; } void LauKMatrixPropagator::getTMatrix(const LauKinematics* kinematics) { // Find the unitary T matrix, where T = [sqrt(rho)]^{*} T_hat sqrt(rho), // and T_hat = (I - i K rho)^-1 * K is the Lorentz-invariant T matrix, // which has phase-space factors included (rho). This function is not // needed to calculate the K-matrix amplitudes, but allows us // to check the variation of T as a function of s (kinematics) if (!kinematics) {return;} // Get the invariant mass squared (s) from the kinematics object. // Use the resPairAmpInt to find which mass-squared combination to use. Double_t s(0.0); if (resPairAmpInt_ == 1) { s = kinematics->getm23Sq(); } else if (resPairAmpInt_ == 2) { s = kinematics->getm13Sq(); } else if (resPairAmpInt_ == 3) { s = kinematics->getm12Sq(); } this->getTMatrix(s); } void LauKMatrixPropagator::getTMatrix(const Double_t s) { // Find the unitary transition T matrix, where // T = [sqrt(rho)]^{*} T_hat sqrt(rho), and // T_hat = (I - i K rho)^-1 * K is the Lorentz-invariant T matrix, // which has phase-space factors included (rho). Note that the first // sqrt of the rho matrix is complex conjugated. // This function is not needed to calculate the K-matrix amplitudes, but // allows us to check the variation of T as a function of s (kinematics) // Initialse the real and imaginary parts of the T matrix to zero ReTMatrix_.Zero(); ImTMatrix_.Zero(); if (parametersSet_ == kFALSE) {return;} // Update K, rho and the propagator (I - i K rho)^-1 this->updatePropagator(s); // Find the real and imaginary T_hat matrices TMatrixD THatReal = realProp_*ScattKMatrix_; TMatrixD THatImag(zeroMatrix_); THatImag -= negImagProp_*ScattKMatrix_; // Find the square-root of the phase space matrix this->getSqrtRhoMatrix(); // Let sqrt(rho) = A + iB and T_hat = C + iD // => T = A(CA-DB) + B(DA+CB) + i[A(DA+CB) + B(DB-CA)] TMatrixD CA(THatReal); CA *= ReSqrtRhoMatrix_; TMatrixD DA(THatImag); DA *= ReSqrtRhoMatrix_; TMatrixD CB(THatReal); CB *= ImSqrtRhoMatrix_; TMatrixD DB(THatImag); DB *= ImSqrtRhoMatrix_; TMatrixD CAmDB(CA); CAmDB -= DB; TMatrixD DApCB(DA); DApCB += CB; TMatrixD DBmCA(DB); DBmCA -= CA; // Find the real and imaginary parts of the transition matrix T ReTMatrix_ = ReSqrtRhoMatrix_*CAmDB + ImSqrtRhoMatrix_*DApCB; ImTMatrix_ = ReSqrtRhoMatrix_*DApCB + ImSqrtRhoMatrix_*DBmCA; } void LauKMatrixPropagator::getSqrtRhoMatrix() { // Find the square root of the (current) phase space matrix so that // we can find T = [sqrt(rho)}^{*} T_hat sqrt(rho), where T_hat is the // Lorentz-invariant T matrix = (I - i K rho)^-1 * K; note that the first // sqrt of rho matrix is complex conjugated // If rho = rho_i + i rho_r = a + i b, then sqrt(rho) = c + i d, where // c = sqrt(0.5*(r+a)) and d = sqrt(0.5(r-a)), where r = sqrt(a^2 + b^2). // Since rho is diagonal, then the square root of rho will also be diagonal, // with its real and imaginary matrix elements equal to c and d, respectively // Initialise the real and imaginary parts of the square root of // the rho matrix to zero ReSqrtRhoMatrix_.Zero(); ImSqrtRhoMatrix_.Zero(); - for (Int_t iChannel (0); iChannel < nChannels_; ++iChannel) { + for (UInt_t iChannel (0); iChannel < nChannels_; ++iChannel) { Double_t realRho = ReRhoMatrix_[iChannel][iChannel]; Double_t imagRho = ImRhoMatrix_[iChannel][iChannel]; Double_t rhoMag = sqrt(realRho*realRho + imagRho*imagRho); Double_t rhoSum = rhoMag + realRho; Double_t rhoDiff = rhoMag - realRho; Double_t reSqrtRho(0.0), imSqrtRho(0.0); if (rhoSum > 0.0) {reSqrtRho = sqrt(0.5*rhoSum);} if (rhoDiff > 0.0) {imSqrtRho = sqrt(0.5*rhoDiff);} ReSqrtRhoMatrix_[iChannel][iChannel] = reSqrtRho; ImSqrtRhoMatrix_[iChannel][iChannel] = imSqrtRho; } } -LauComplex LauKMatrixPropagator::getTHat(const Double_t s, const Int_t channel) { +LauComplex LauKMatrixPropagator::getTHat(const Double_t s, const UInt_t channel) { LauComplex THat(0.0, 0.0); if (channel <= 0 || channel > nChannels_) {return THat;} this->updatePropagator(s); // Find the real and imaginary T_hat matrices TMatrixD THatReal = realProp_*ScattKMatrix_; TMatrixD THatImag(zeroMatrix_); THatImag -= negImagProp_*ScattKMatrix_; // Return the specific THat component THat.setRealPart(THatReal[index_][channel-1]); THat.setImagPart(THatImag[index_][channel-1]); return THat; } nlohmann::json LauKMatrixPropagator::writeSettingsToJson() const { using nlohmann::json; json j = json::object(); j[ "propName" ] = name_.Data(); j[ "paramFileName" ] = paramFileName_.Data(); j[ "resPairAmpInt" ] = resPairAmpInt_; j[ "nChannels" ] = nChannels_; j[ "nPoles" ] = nPoles_; if ( index_ > 0 ) { j[ "rowIndex" ] = index_+1; } if ( ! includeBWBarrierFactor_ ) { j[ "ignoreBWBarrierFactor" ] = true; } return j; }