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 | |
40 | using namespace std; |
41 | |
42 | |
43 | vector<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 | |
60 | typedef TwoZPiPlotGenerator PlotGen; |
61 | |
62 | void 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 | |
76 | using namespace std; |
77 | |
78 | int main( int argc, char* argv[] ){ |
79 | |
80 | |
81 | |
82 | |
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 | |
117 | |
118 | |
119 | cout << "Fit results file name = " << resultsName << endl; |
120 | cout << "Output file name = " << outName << endl << endl; |
121 | |
122 | |
123 | |
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 | |
135 | |
136 | |
137 | atiSetup(); |
138 | PlotGen plotGen( results ); |
139 | |
140 | |
141 | |
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 | |
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 | |
163 | |
164 | |
165 | string reactionName = results.reactionList()[0]; |
166 | plotGen.enableReaction( reactionName ); |
167 | vector<string> sums = plotGen.uniqueSums(); |
168 | |
169 | |
170 | |
171 | |
172 | plotGen.enableSum(0); |
173 | plotGen.enableSum(1); |
174 | |
175 | |
176 | |
177 | |
178 | |
179 | |
180 | |
181 | for (unsigned int iamp = 0; iamp < amplist.size(); iamp++){ |
182 | |
183 | |
184 | for (unsigned int i = 0; i < amplist.size(); i++){ |
185 | plotGen.enableAmp(i); |
186 | } |
187 | } |
188 | |
189 | |
190 | |
191 | for (unsigned int iplot = 0; iplot < PlotGenerator::kNumTypes; iplot++){ |
192 | |
193 | |
194 | for (unsigned int ivar = 0; ivar < TwoZPiPlotGenerator::kNumHists; ivar++){ |
195 | |
196 | |
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 | |
224 | for (unsigned int iamp = 0; iamp < amplist.size(); iamp++){ |
225 | |
226 | |
227 | for (unsigned int i = 0; i < amplist.size(); i++){ |
228 | plotGen.enableAmp(i); |
229 | } |
230 | |
231 | |
232 | for (unsigned int i = 0; i < amplist.size(); i++){ |
233 | if (i != iamp) plotGen.disableAmp(i); |
234 | } |
235 | |
236 | |
237 | |
238 | for (unsigned int iplot = 0; iplot < PlotGenerator::kNumTypes; iplot++){ |
239 | |
240 | |
241 | for (unsigned int ivar = 0; ivar < TwoZPiPlotGenerator::kNumHists; ivar++){ |
242 | |
243 | |
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 | |
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 | |
279 | |
280 | |
281 | |
282 | |
283 | |
284 | vector< string > pars; |
285 | |
286 | |
287 | |
288 | |
289 | |
290 | |
291 | |
292 | |
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 | |
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 | |
319 | |
320 | |
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); |
| Value stored to 'Sigma_err' during its initialization is never read |
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]); |
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 | |
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 | |
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 | |