1 | #include <iostream> |
2 | #include <fstream> |
3 | #include <string> |
4 | #include <vector> |
5 | |
6 | #include "TClass.h" |
7 | #include "TApplication.h" |
8 | #include "TGClient.h" |
9 | #include "TROOT.h" |
10 | #include "TH1.h" |
11 | #include "TStyle.h" |
12 | #include "TClass.h" |
13 | #include "TFile.h" |
14 | |
15 | #include "IUAmpTools/AmpToolsInterface.h" |
16 | #include "IUAmpTools/FitResults.h" |
17 | |
18 | #include "AmpPlotter/PlotterMainWindow.h" |
19 | #include "AmpPlotter/PlotFactory.h" |
20 | |
21 | #include "AMPTOOLS_DATAIO/TwoPiPlotGenerator.h" |
22 | #include "AMPTOOLS_DATAIO/ROOTDataReader.h" |
23 | #include "AMPTOOLS_AMPS/TwoPiAngles.h" |
24 | #include "AMPTOOLS_AMPS/TwoPiAngles_amp.h" |
25 | #include "AMPTOOLS_AMPS/BreitWigner.h" |
26 | |
27 | typedef TwoPiPlotGenerator PlotGen; |
28 | |
29 | void atiSetup(){ |
30 | |
31 | AmpToolsInterface::registerAmplitude( TwoPiAngles() ); |
32 | AmpToolsInterface::registerAmplitude( TwoPiAngles_amp() ); |
33 | AmpToolsInterface::registerAmplitude( BreitWigner() ); |
34 | AmpToolsInterface::registerDataReader( ROOTDataReader() ); |
35 | } |
36 | |
37 | using namespace std; |
38 | |
39 | int main( int argc, char* argv[] ){ |
40 | |
41 | |
42 | |
43 | |
44 | |
45 | |
46 | cout << endl << " *** Viewing Results Using AmpPlotter and writing root histograms *** " << endl << endl; |
47 | |
48 | if (argc < 2){ |
49 | cout << "Usage:" << endl << endl; |
50 | cout << "\ttwopi_plotter <results file name> -o <output file name>" << endl << endl; |
51 | return 0; |
52 | } |
53 | |
54 | bool showGui = false; |
55 | string outName = "twopi_plot.root"; |
56 | string resultsName(argv[1]); |
57 | for (int i = 2; i < argc; i++){ |
58 | |
59 | string arg(argv[i]); |
60 | |
61 | if (arg == "-g"){ |
62 | showGui = true; |
63 | } |
64 | if (arg == "-o"){ |
65 | outName = argv[++i]; |
66 | } |
67 | if (arg == "-h"){ |
68 | cout << endl << " Usage for: " << argv[0] << endl << endl; |
69 | cout << "\t -o <file>\t output file path" << endl; |
70 | cout << "\t -g <file>\t show GUI" << endl; |
71 | exit(1); |
72 | } |
73 | } |
74 | |
75 | |
76 | |
77 | |
78 | |
79 | |
80 | cout << "Fit results file name = " << resultsName << endl; |
81 | cout << "Output file name = " << outName << endl << endl; |
82 | |
83 | |
84 | |
85 | |
86 | |
87 | FitResults results( resultsName ); |
88 | if( !results.valid() ){ |
89 | |
90 | cout << "Invalid fit results in file: " << resultsName << endl; |
91 | exit( 1 ); |
92 | } |
93 | |
94 | |
95 | |
96 | |
97 | |
98 | atiSetup(); |
99 | PlotGen plotGen( results ); |
100 | |
101 | |
102 | |
103 | |
104 | |
105 | TFile* plotfile = new TFile( outName.c_str(), "recreate"); |
106 | TH1::AddDirectory(kFALSE); |
107 | |
108 | string reactionName = results.reactionList()[0]; |
109 | plotGen.enableReaction( reactionName ); |
110 | vector<string> sums = plotGen.uniqueSums(); |
111 | |
112 | |
113 | |
114 | for (unsigned int isum = 0; isum <= sums.size(); isum++){ |
115 | |
116 | |
117 | for (unsigned int i = 0; i < sums.size(); i++){ |
118 | plotGen.enableSum(i); |
119 | } |
120 | |
121 | |
122 | if (isum < sums.size()){ |
123 | for (unsigned int i = 0; i < sums.size(); i++){ |
124 | if (i != isum) plotGen.disableSum(i); |
125 | } |
126 | } |
127 | |
128 | |
129 | |
130 | for (unsigned int iplot = 0; iplot < PlotGenerator::kNumTypes; iplot++){ |
131 | if (isum < sums.size() && iplot == PlotGenerator::kData) continue; |
132 | |
133 | |
134 | for (unsigned int ivar = 0; ivar < TwoPiPlotGenerator::kNumHists; ivar++){ |
135 | |
136 | |
137 | string histname = ""; |
138 | if (ivar == TwoPiPlotGenerator::k2PiMass) histname += "M2pi"; |
139 | else if (ivar == TwoPiPlotGenerator::kPiPCosTheta) histname += "cosTheta"; |
140 | else if (ivar == TwoPiPlotGenerator::kPhi) histname += "Phi"; |
141 | else if (ivar == TwoPiPlotGenerator::kphi) histname += "phi"; |
142 | else if (ivar == TwoPiPlotGenerator::kPsi) histname += "psi"; |
143 | else if (ivar == TwoPiPlotGenerator::kt) histname += "t"; |
144 | else continue; |
145 | |
146 | if (iplot == PlotGenerator::kData) histname += "dat"; |
147 | if (iplot == PlotGenerator::kAccMC) histname += "acc"; |
148 | if (iplot == PlotGenerator::kGenMC) histname += "gen"; |
149 | |
150 | if (isum < sums.size()){ |
151 | |
152 | |
153 | |
154 | |
155 | string sumName = sums[isum]; |
156 | histname += "_"; |
157 | histname += sumName; |
158 | } |
159 | |
160 | Histogram* hist = plotGen.projection(ivar, reactionName, iplot); |
161 | TH1* thist = hist->toRoot(); |
162 | thist->SetName(histname.c_str()); |
163 | plotfile->cd(); |
164 | thist->Write(); |
165 | |
166 | } |
167 | } |
168 | } |
169 | |
170 | plotfile->Close(); |
171 | |
172 | |
173 | |
174 | |
175 | |
176 | |
177 | vector< string > pars; |
178 | pars.push_back("Pi+Pi-::helplusN+::g1VM1_re"); |
179 | pars.push_back("Pi+Pi-::helplusN+::g1VM0_re"); |
180 | pars.push_back("Pi+Pi-::helplusN+::g1VM0_im"); |
181 | |
182 | pars.push_back("Pi+Pi-::helplusN+::g1VM-1_re"); |
183 | pars.push_back("Pi+Pi-::helplusN+::g1VM-1_im"); |
184 | pars.push_back("Pi+Pi-::helplusN+::g-1VM1_re"); |
185 | pars.push_back("Pi+Pi-::helplusN+::g-1VM1_im"); |
186 | |
187 | pars.push_back("Pi+Pi-::helplusN+::g-1VM0_re"); |
188 | pars.push_back("Pi+Pi-::helplusN+::g-1VM0_im"); |
189 | pars.push_back("Pi+Pi-::helplusN+::g-1VM-1_re"); |
190 | pars.push_back("Pi+Pi-::helplusN+::g-1VM-1_im"); |
191 | |
192 | |
193 | ofstream outfile; |
194 | outfile.open( "twopi_fitPars.txt" ); |
195 | |
196 | for(unsigned int i = 0; i<pars.size(); i++) { |
197 | double parValue = results.parValue( pars[i] ); |
198 | double parError = results.parError( pars[i] ); |
199 | outfile << parValue << "\t" << parError << "\t"; |
200 | } |
201 | |
202 | |
203 | |
204 | |
205 | vector< vector< double > > covMatrix; |
206 | covMatrix = results.errorMatrix(); |
207 | |
208 | double SigmaN = results.parValue(pars[3]) + results.parValue(pars[6]); |
209 | double SigmaN_err = covMatrix[5][5] + covMatrix[8][8] + 2*covMatrix[5][8]; |
210 | |
211 | double SigmaD = 0.5*(1 - results.parValue(pars[0])) + results.parValue(pars[2]); |
212 | double SigmaD_err = 0.5*0.5*covMatrix[2][2] + covMatrix[4][4] - 2*0.5*covMatrix[2][4]; |
213 | |
214 | double Sigma = SigmaN/SigmaD; |
215 | 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 |
216 | |
217 | double P = 2*results.parValue(pars[6]) - results.parValue(pars[4]); |
218 | double P_err = sqrt(2*2*covMatrix[8][8] + covMatrix[6][6] - 2*2*covMatrix[6][8]); |
219 | |
220 | Sigma = Sigma_err = P = P_err = 0; |
221 | outfile << Sigma << "\t" << Sigma_err << "\t"; |
222 | outfile << P << "\t" << P_err << "\t"; |
223 | |
224 | outfile << endl; |
225 | |
226 | |
227 | |
228 | |
229 | |
230 | if(showGui) { |
231 | |
232 | cout << ">> Plot generator ready, starting GUI..." << endl; |
233 | |
234 | int dummy_argc = 0; |
235 | char* dummy_argv[] = {}; |
236 | TApplication app( "app", &dummy_argc, dummy_argv ); |
237 | |
238 | gStyle->SetFillColor(10); |
239 | gStyle->SetCanvasColor(10); |
240 | gStyle->SetPadColor(10); |
241 | gStyle->SetFillStyle(1001); |
242 | gStyle->SetPalette(1); |
243 | gStyle->SetFrameFillColor(10); |
244 | gStyle->SetFrameFillStyle(1001); |
245 | |
246 | PlotFactory factory( plotGen ); |
247 | PlotterMainWindow mainFrame( gClient(TGClient::Instance())->GetRoot(), factory ); |
248 | |
249 | app.Run(); |
250 | } |
251 | |
252 | return 0; |
253 | |
254 | } |
255 | |