Bug Summary

File:programs/AmplitudeAnalysis/twopi_plotter_primakoff/twopi_plotter_primakoff.cc
Location:line 334, column 10
Description:Value stored to 'P_err' during its initialization is never read

Annotated Source Code

1#include <iostream>
2#include <fstream>
3#include <string>
4#include <vector>
5#include <complex>
6#include <algorithm>
7#include <cassert>
8
9#include "TClass.h"
10#include "TApplication.h"
11#include "TGClient.h"
12#include "TROOT.h"
13#include "TH1.h"
14#include "TStyle.h"
15#include "TClass.h"
16#include "TFile.h"
17
18#include "IUAmpTools/PlotGenerator.h"
19#include "IUAmpTools/AmplitudeManager.h"
20#include "IUAmpTools/NormIntInterface.h"
21#include "IUAmpTools/AmpToolsInterface.h"
22#include "IUAmpTools/FitResults.h"
23
24#include "AmpPlotter/PlotterMainWindow.h"
25#include "AmpPlotter/PlotFactory.h"
26
27#include "AMPTOOLS_DATAIO/TwoZPiPlotGenerator.h"
28#include "AMPTOOLS_DATAIO/ROOTDataReader.h"
29#include "AMPTOOLS_AMPS/TwoPiAngles.h"
30#include "AMPTOOLS_AMPS/TwoPiWt_primakoff.h"
31#include "AMPTOOLS_AMPS/TwoPiWt_sigma.h"
32#include "AMPTOOLS_AMPS/TwoPiW_brokenetas.h"
33#include "AMPTOOLS_AMPS/TwoPitdist.h"
34#include "AMPTOOLS_AMPS/TwoPiNC_tdist.h"
35#include "AMPTOOLS_AMPS/TwoPiEtas_tdist.h"
36#include "AMPTOOLS_AMPS/TwoPiAngles_primakoff.h"
37#include "AMPTOOLS_AMPS/BreitWigner.h"
38
39
40using namespace std;
41
42// for string delimiter
43vector<string> stringSplit (string s, string delimiter) {
44 size_t pos_start = 0, pos_end, delim_len = delimiter.length();
45 string token;
46 vector<string> res;
47
48 while ((pos_end = s.find (delimiter, pos_start)) != string::npos) {
49 token = s.substr (pos_start, pos_end - pos_start);
50 pos_start = pos_end + delim_len;
51 res.push_back (token);
52 }
53
54 res.push_back (s.substr (pos_start));
55 return res;
56}
57
58
59
60typedef TwoZPiPlotGenerator PlotGen;
61
62void atiSetup(){
63
64 AmpToolsInterface::registerAmplitude( TwoPiAngles() );
65 AmpToolsInterface::registerAmplitude( TwoPiAngles_primakoff() );
66 AmpToolsInterface::registerAmplitude( TwoPiWt_primakoff() );
67 AmpToolsInterface::registerAmplitude( TwoPiWt_sigma() );
68 AmpToolsInterface::registerAmplitude( TwoPiW_brokenetas() );
69 AmpToolsInterface::registerAmplitude( TwoPitdist() );
70 AmpToolsInterface::registerAmplitude( TwoPiNC_tdist() );
71 AmpToolsInterface::registerAmplitude( TwoPiEtas_tdist() );
72 AmpToolsInterface::registerAmplitude( BreitWigner() );
73 AmpToolsInterface::registerDataReader( ROOTDataReader() );
74}
75
76using namespace std;
77
78int main( int argc, char* argv[] ){
79
80
81 // ************************
82 // usage
83 // ************************
84
85 cout << endl << " *** Viewing Results Using AmpPlotter and writing root histograms *** " << endl << endl;
86
87 if (argc < 2){
88 cout << "Usage:" << endl << endl;
89 cout << "\ttwopi_plotter <results file name> -o <output file name>" << endl << endl;
90 return 0;
91 }
92
93 bool showGui = false;
94 string outName = "twopi_plot.root";
95 string resultsName(argv[1]);
96 for (int i = 2; i < argc; i++){
97
98 string arg(argv[i]);
99
100 if (arg == "-g"){
101 showGui = true;
102 }
103 if (arg == "-o"){
104 outName = argv[++i];
105 }
106 if (arg == "-h"){
107 cout << endl << " Usage for: " << argv[0] << endl << endl;
108 cout << "\t -o <file>\t output file path" << endl;
109 cout << "\t -g <file>\t show GUI" << endl;
110 exit(1);
111 }
112 }
113
114
115 // ************************
116 // parse the command line parameters
117 // ************************
118
119 cout << "Fit results file name = " << resultsName << endl;
120 cout << "Output file name = " << outName << endl << endl;
121
122 // ************************
123 // load the results and display the configuration info
124 // ************************
125
126 FitResults results( resultsName );
127 if( !results.valid() ){
128
129 cout << "Invalid fit results in file: " << resultsName << endl;
130 exit( 1 );
131 }
132
133 // ************************
134 // set up the plot generator
135 // ************************
136
137 atiSetup();
138 PlotGen plotGen( results );
139
140 // ************************
141 // set up an output ROOT file to store histograms
142 // ************************
143
144 TFile* plotfile = new TFile( outName.c_str(), "recreate");
145 TH1::AddDirectory(kFALSE);
146
147 vector <string> paramlist;
148 paramlist = results.ampList("Primakoff");
149
150 vector <string> amplist;
151
152 for (auto parsum : paramlist) {
153 // cout << " sum segment=" << parsum << endl;
154 vector<string> parbreak = stringSplit (parsum, "::");
155
156 if (parbreak[1] == "Aplus" || parbreak[1] == "IAplus" || parbreak[1] == "EAplus" ) {
157 amplist.push_back(parbreak[2]);
158 cout << " amp =" << parbreak[2] << endl;
159 }
160 }
161
162 // for (auto amp : amplist) cout << " amp=" << amp << endl;
163
164
165 string reactionName = results.reactionList()[0];
166 plotGen.enableReaction( reactionName );
167 vector<string> sums = plotGen.uniqueSums();
168
169 //for (auto isum : sums) cout << " isum=" << isum << endl;
170
171 // Enable both Aplus and Aminus sum
172 plotGen.enableSum(0);
173 plotGen.enableSum(1);
174 // for (unsigned int isum = 0; isum < sums.size(); isum++){
175 // plotGen.disableSum(isum);
176 // }
177
178
179
180 // loop over sum, with all amplitudes turned on
181 for (unsigned int iamp = 0; iamp < amplist.size(); iamp++){
182
183 // turn on all amplist by default
184 for (unsigned int i = 0; i < amplist.size(); i++){
185 plotGen.enableAmp(i);
186 }
187 }
188
189
190 // loop over sum, accMC, and genMC and kBkgnd
191 for (unsigned int iplot = 0; iplot < PlotGenerator::kNumTypes; iplot++){
192
193 // loop over different variables
194 for (unsigned int ivar = 0; ivar < TwoZPiPlotGenerator::kNumHists; ivar++){
195
196 // set unique histogram name for each plot (could put in directories...)
197 string histname = "";
198 if (ivar == TwoZPiPlotGenerator::k2PiMass) histname += "M2pi";
199 else if (ivar == TwoZPiPlotGenerator::kPiPCosTheta) histname += "cosTheta";
200 else if (ivar == TwoZPiPlotGenerator::kPhi) histname += "Phi";
201 else if (ivar == TwoZPiPlotGenerator::kphi) histname += "phi";
202 else if (ivar == TwoZPiPlotGenerator::kPsi) histname += "psi";
203 else if (ivar == TwoZPiPlotGenerator::kt) histname += "t";
204 else if (ivar == TwoZPiPlotGenerator::ktheta_scat) histname += "theta_scat";
205 else continue;
206
207 if (iplot == PlotGenerator::kData) histname += "dat";
208 if (iplot == PlotGenerator::kAccMC) histname += "acc";
209 if (iplot == PlotGenerator::kGenMC) histname += "gen";
210 if (iplot == PlotGenerator::kBkgnd) histname += "bkgnd";
211
212 Histogram* hist = plotGen.projection(ivar, reactionName, iplot);
213 TH1* thist = hist->toRoot();
214 thist->SetName(histname.c_str());
215 plotfile->cd();
216 thist->Write();
217
218 }
219 }
220
221
222
223 // loop over sum, once per each amplitude
224 for (unsigned int iamp = 0; iamp < amplist.size(); iamp++){
225
226 // turn on all amplist by default
227 for (unsigned int i = 0; i < amplist.size(); i++){
228 plotGen.enableAmp(i);
229 }
230
231 // for individual contributions turn off all amplist but the one of interest
232 for (unsigned int i = 0; i < amplist.size(); i++){
233 if (i != iamp) plotGen.disableAmp(i);
234 }
235
236
237 // loop over data, accMC, and genMC and kBkgnd
238 for (unsigned int iplot = 0; iplot < PlotGenerator::kNumTypes; iplot++){
239
240 // loop over different variables
241 for (unsigned int ivar = 0; ivar < TwoZPiPlotGenerator::kNumHists; ivar++){
242
243 // set unique histogram name for each plot (could put in directories...)
244 string histname = "";
245 if (ivar == TwoZPiPlotGenerator::k2PiMass) histname += "M2pi";
246 else if (ivar == TwoZPiPlotGenerator::kPiPCosTheta) histname += "cosTheta";
247 else if (ivar == TwoZPiPlotGenerator::kPhi) histname += "Phi";
248 else if (ivar == TwoZPiPlotGenerator::kphi) histname += "phi";
249 else if (ivar == TwoZPiPlotGenerator::kPsi) histname += "psi";
250 else if (ivar == TwoZPiPlotGenerator::kt) histname += "t";
251 else if (ivar == TwoZPiPlotGenerator::ktheta_scat) histname += "theta_scat";
252 else continue;
253
254 if (iplot == PlotGenerator::kData) histname += "dat";
255 if (iplot == PlotGenerator::kAccMC) histname += "acc";
256 if (iplot == PlotGenerator::kGenMC) histname += "gen";
257 if (iplot == PlotGenerator::kBkgnd) histname += "bkgnd";
258
259 // get name of amp for naming histogram
260 string ampName = amplist[iamp];
261 histname += "_";
262 histname += ampName;
263
264 Histogram* hist = plotGen.projection(ivar, reactionName, iplot);
265 TH1* thist = hist->toRoot();
266 thist->SetName(histname.c_str());
267 plotfile->cd();
268 thist->Write();
269
270 }
271 }
272 }
273
274
275 plotfile->Close();
276
277 // ************************
278 // retrieve amplitudes for output
279 // ************************
280
281 // get parameter list
282
283 // parameters to check
284 vector< string > pars;
285 /* pars.push_back("Primakoff::Aplus::g1V00_re");
286 pars.push_back("Primakoff::Aplus::g1V00_im");
287 pars.push_back("Primakoff::Aplus::g1V11_re");
288 pars.push_back("Primakoff::Aplus::g1V11_im");
289 pars.push_back("Primakoff::Aplus::g1V10_re");
290 pars.push_back("Primakoff::Aplus::g1V10_im");
291 pars.push_back("Primakoff::Aplus::g1V1-1_re");
292 pars.push_back("Primakoff::Aplus::g1V1-1_im");*/
293
294 vector <string> parlist;
295 parlist = results.ampList("Primakoff");
296 for(unsigned int j=0; j<parlist.size(); j++) {
297 cout << " j=" << j << " parlist[j]=" << parlist[j] << " " << results.realProdParName(parlist[j]) << " " << results.imagProdParName(parlist[j]) << endl;
298 if ((parlist[j].find("Aplus") != string::npos) ) {
299 pars.push_back(results.realProdParName(parlist[j]));
300 pars.push_back(results.imagProdParName(parlist[j]));
301 }
302 }
303
304 // file for writing parameters (later switch to putting in ROOT file)
305 ofstream outfile;
306 outfile.open( "twopi_fitPars.txt" );
307 cout << "Opened Output File twopi_fitPars.txt" << " pars.size=" << pars.size() << endl;
308
309 for(unsigned int i = 0; i<pars.size(); i++) {
310 double parValue = results.parValue( pars[i] );
311 double parError = results.parError( pars[i] );
312 int ifindg = pars[i].find("g");
313 outfile << pars[i].substr(ifindg) << "\t" << parValue << "\t" << parError << "\t";
314 }
315
316
317
318 // Note: For twopi_primakoff_plotter: The following computations are nonsense for amplitudes
319
320 // covariance matrix
321 vector< vector< double > > covMatrix;
322 covMatrix = results.errorMatrix();
323
324 double SigmaN = results.parValue(pars[0]) + results.parValue(pars[0]);
325 double SigmaN_err = covMatrix[0][0] + covMatrix[0][0] + 2*covMatrix[0][0];
326
327 double SigmaD = 0.5*(1 - results.parValue(pars[0])) + results.parValue(pars[0]);
328 double SigmaD_err = 0.5*0.5*covMatrix[0][0] + covMatrix[0][0] - 2*0.5*covMatrix[0][0];
329
330 double Sigma = SigmaN/SigmaD;
331 double Sigma_err = fabs(Sigma) * sqrt(SigmaN_err/SigmaN/SigmaN + SigmaD_err/SigmaD/SigmaD);
332
333 double P = 2*results.parValue(pars[0]) - results.parValue(pars[0]);
334 double P_err = sqrt(2*2*covMatrix[0][0] + covMatrix[0][0] - 2*2*covMatrix[0][0]);
Value stored to 'P_err' during its initialization is never read
335
336 Sigma = Sigma_err = P = P_err = 0;
337 outfile << "Sigma" << "\t" << Sigma << "\t" << Sigma_err << "\t";
338 outfile << "P" << "\t" << P << "\t" << P_err << "\t";
339
340 // output covariance matrix. Output only half since A+ and A- are constrained to be the same.
341 for (unsigned int j=0; j< covMatrix.size()/2; j++) {
342 outfile << endl;
343 for (unsigned int jj=0; jj< covMatrix.size()/2; jj++) {
344 outfile.width(20);
345 outfile << covMatrix[j][jj];
346 }
347 }
348
349
350 outfile << endl;
351
352 // ************************
353 // start the GUI
354 // ************************
355
356 if(showGui) {
357
358 cout << ">> Plot generator ready, starting GUI..." << endl;
359
360 int dummy_argc = 0;
361 char* dummy_argv[] = {};
362 TApplication app( "app", &dummy_argc, dummy_argv );
363
364 gStyle->SetFillColor(10);
365 gStyle->SetCanvasColor(10);
366 gStyle->SetPadColor(10);
367 gStyle->SetFillStyle(1001);
368 gStyle->SetPalette(1);
369 gStyle->SetFrameFillColor(10);
370 gStyle->SetFrameFillStyle(1001);
371
372 PlotFactory factory( plotGen );
373 PlotterMainWindow mainFrame( gClient(TGClient::Instance())->GetRoot(), factory );
374
375 app.Run();
376 }
377
378 return 0;
379
380}
381