1 | #include <iostream> |
2 | #include <fstream> |
3 | #include <sstream> |
4 | #include <complex> |
5 | #include <string> |
6 | #include <vector> |
7 | #include <utility> |
8 | #include <map> |
9 | |
10 | #include "TSystem.h" |
11 | |
12 | #include "AMPTOOLS_DATAIO/ROOTDataReader.h" |
13 | #include "AMPTOOLS_DATAIO/ROOTDataReaderBootstrap.h" |
14 | #include "AMPTOOLS_DATAIO/ROOTDataReaderWithTCut.h" |
15 | #include "AMPTOOLS_DATAIO/ROOTDataReaderTEM.h" |
16 | #include "AMPTOOLS_AMPS/TwoPSAngles.h" |
17 | #include "AMPTOOLS_AMPS/TwoPSHelicity.h" |
18 | #include "AMPTOOLS_AMPS/TwoPiAngles.h" |
19 | #include "AMPTOOLS_AMPS/TwoPiAngles_amp.h" |
20 | #include "AMPTOOLS_AMPS/TwoPiWt_primakoff.h" |
21 | #include "AMPTOOLS_AMPS/TwoPiWt_sigma.h" |
22 | #include "AMPTOOLS_AMPS/TwoPiW_brokenetas.h" |
23 | #include "AMPTOOLS_AMPS/TwoPitdist.h" |
24 | #include "AMPTOOLS_AMPS/TwoPiNC_tdist.h" |
25 | #include "AMPTOOLS_AMPS/TwoPiEtas_tdist.h" |
26 | #include "AMPTOOLS_AMPS/TwoPiAngles_primakoff.h" |
27 | #include "AMPTOOLS_AMPS/ThreePiAngles.h" |
28 | #include "AMPTOOLS_AMPS/ThreePiAnglesSchilling.h" |
29 | #include "AMPTOOLS_AMPS/TwoPiAnglesRadiative.h" |
30 | #include "AMPTOOLS_AMPS/Zlm.h" |
31 | #include "AMPTOOLS_AMPS/BreitWigner.h" |
32 | #include "AMPTOOLS_AMPS/BreitWigner3body.h" |
33 | #include "AMPTOOLS_AMPS/b1piAngAmp.h" |
34 | #include "AMPTOOLS_AMPS/omegapiAngAmp.h" |
35 | #include "AMPTOOLS_AMPS/Uniform.h" |
36 | #include "AMPTOOLS_AMPS/polCoef.h" |
37 | #include "AMPTOOLS_AMPS/dblRegge.h" |
38 | #include "AMPTOOLS_AMPS/dblReggeMod.h" |
39 | #include "AMPTOOLS_AMPS/omegapi_amplitude.h" |
40 | #include "AMPTOOLS_AMPS/Vec_ps_refl.h" |
41 | #include "AMPTOOLS_AMPS/Piecewise.h" |
42 | |
43 | #include "MinuitInterface/MinuitMinimizationManager.h" |
44 | #include "IUAmpTools/AmpToolsInterface.h" |
45 | #include "IUAmpTools/FitResults.h" |
46 | #include "IUAmpTools/ConfigFileParser.h" |
47 | #include "IUAmpTools/ConfigurationInfo.h" |
48 | |
49 | using std::complex; |
50 | using namespace std; |
51 | |
52 | double runSingleFit(ConfigurationInfo* cfgInfo, bool useMinos, int maxIter, string seedfile) { |
53 | AmpToolsInterface ati( cfgInfo ); |
54 | |
55 | cout << "LIKELIHOOD BEFORE MINIMIZATION: " << ati.likelihood() << endl; |
56 | |
57 | MinuitMinimizationManager* fitManager = ati.minuitMinimizationManager(); |
58 | fitManager->setMaxIterations(maxIter); |
59 | |
60 | if( useMinos ){ |
61 | |
62 | fitManager->minosMinimization(); |
63 | } |
64 | else{ |
65 | |
66 | fitManager->migradMinimization(); |
67 | } |
68 | |
69 | bool fitFailed = |
70 | ( fitManager->status() != 0 && fitManager->eMatrixStatus() != 3 ); |
71 | |
72 | if( fitFailed ){ |
73 | cout << "ERROR: fit failed use results with caution..." << endl; |
74 | return 1e6; |
75 | } |
76 | |
77 | cout << "LIKELIHOOD AFTER MINIMIZATION: " << ati.likelihood() << endl; |
78 | |
79 | ati.finalizeFit(); |
80 | |
81 | if( seedfile.size() != 0 && !fitFailed ){ |
82 | ati.fitResults()->writeSeed( seedfile ); |
83 | } |
84 | |
85 | return ati.likelihood(); |
86 | } |
87 | |
88 | void runRndFits(ConfigurationInfo* cfgInfo, bool useMinos, int maxIter, string seedfile, int numRnd, double maxFraction) { |
89 | AmpToolsInterface ati( cfgInfo ); |
90 | string fitName = cfgInfo->fitName(); |
91 | |
92 | cout << "LIKELIHOOD BEFORE MINIMIZATION: " << ati.likelihood() << endl; |
93 | |
94 | MinuitMinimizationManager* fitManager = ati.minuitMinimizationManager(); |
95 | fitManager->setMaxIterations(maxIter); |
96 | |
97 | vector< vector<string> > parRangeKeywords = cfgInfo->userKeywordArguments("parRange"); |
98 | |
99 | |
100 | double minLL = 0; |
101 | int minFitTag = -1; |
102 | |
103 | for(int i=0; i<numRnd; i++) { |
104 | cout << endl << "###############################" << endl; |
105 | cout << "FIT " << i << " OF " << numRnd << endl; |
106 | cout << endl << "###############################" << endl; |
107 | |
108 | |
109 | ati.randomizeProductionPars(maxFraction); |
110 | for(size_t ipar=0; ipar<parRangeKeywords.size(); ipar++) { |
111 | ati.randomizeParameter(parRangeKeywords[ipar][0], atof(parRangeKeywords[ipar][1].c_str()), atof(parRangeKeywords[ipar][2].c_str())); |
112 | } |
113 | |
114 | if(useMinos) |
115 | fitManager->minosMinimization(); |
116 | else |
117 | fitManager->migradMinimization(); |
118 | |
119 | bool fitFailed = (fitManager->status() != 0 && fitManager->eMatrixStatus() != 3); |
120 | |
121 | if( fitFailed ) |
122 | cout << "ERROR: fit failed use results with caution..." << endl; |
123 | |
124 | cout << "LIKELIHOOD AFTER MINIMIZATION: " << ati.likelihood() << endl; |
125 | |
126 | ati.finalizeFit(to_string(i)); |
127 | |
128 | if( seedfile.size() != 0 && !fitFailed ){ |
129 | string seedfile_rand = seedfile + Form("_%d.txt", i); |
130 | ati.fitResults()->writeSeed( seedfile_rand ); |
131 | } |
132 | |
133 | |
134 | if( !fitFailed && ati.likelihood() < minLL ) { |
135 | minLL = ati.likelihood(); |
136 | minFitTag = i; |
137 | } |
138 | } |
139 | |
140 | |
141 | if(minFitTag < 0) cout << "ALL FITS FAILED!" << endl; |
142 | else { |
143 | cout << "MINIMUM LIKELIHOOD FROM " << minFitTag << " of " << numRnd << " RANDOM PRODUCTION PARS = " << minLL << endl; |
144 | gSystem->Exec(Form("cp %s_%d.fit %s.fit", fitName.data(), minFitTag, fitName.data())); |
145 | if( seedfile.size() != 0 ) |
146 | gSystem->Exec(Form("cp %s_%d.txt %s.txt", seedfile.data(), minFitTag, seedfile.data())); |
147 | } |
148 | } |
149 | |
150 | void runParScan(ConfigurationInfo* cfgInfo, bool useMinos, int maxIter, string seedfile, string parScan) { |
151 | double minVal, maxVal, stepSize; |
152 | int steps; |
| 7 | | 'steps' declared without an initial value | |
|
153 | |
154 | vector< vector<string> > parScanKeywords = cfgInfo->userKeywordArguments("parScan"); |
155 | |
156 | if(parScanKeywords.size()==0) { |
| |
157 | cout << "No parScan keyword found in configuration file. Set up at least one parameter for scanning! Aborting." << endl; |
158 | return; |
159 | } else { |
160 | for(size_t ipar=0; ipar<parScanKeywords.size(); ipar++) { |
| 9 | | Loop condition is true. Entering loop body | |
|
| 11 | | Loop condition is false. Execution continues on line 172 | |
|
161 | if(parScanKeywords[ipar][0]==parScan) { |
| |
162 | minVal = atof(parScanKeywords[ipar][1].c_str()); |
163 | maxVal = atof(parScanKeywords[ipar][2].c_str()); |
164 | stepSize = atof(parScanKeywords[ipar][3].c_str()); |
165 | steps = trunc((maxVal-minVal)/stepSize)+1; |
166 | break; |
167 | } else |
168 | cout << "Skipping configuration to scan " << parScanKeywords[ipar][0] << "since scanning of " << parScan << " was requested..." << endl; |
169 | } |
170 | } |
171 | |
172 | AmpToolsInterface ati( cfgInfo ); |
173 | |
174 | string fitName = cfgInfo->fitName(); |
175 | cout << "LIKELIHOOD BEFORE MINIMIZATION: " << ati.likelihood() << endl; |
176 | |
177 | ParameterManager* parMgr = ati.parameterManager(); |
178 | MinuitMinimizationManager* fitManager = ati.minuitMinimizationManager(); |
179 | fitManager->setMaxIterations(maxIter); |
180 | |
181 | |
182 | for(int i=0; i<steps; i++) { |
| 12 | | The right operand of '<' is a garbage value |
|
183 | cout << endl << "###############################" << endl; |
184 | cout << "FIT " << i << " OF " << steps << endl; |
185 | cout << endl << "###############################" << endl; |
186 | |
187 | |
188 | vector<ParameterInfo*> parInfoVec = cfgInfo->parameterList(); |
189 | |
190 | auto parItr = parInfoVec.begin(); |
191 | for( ; parItr != parInfoVec.end(); ++parItr ) { |
192 | if( (**parItr).parName() == parScan ) break; |
193 | } |
194 | |
195 | if( parItr == parInfoVec.end() ){ |
196 | cout << "ERROR: request to scan nonexistent parameter: " << parScan << endl; |
197 | return; |
198 | } |
199 | |
200 | |
201 | double value = minVal + i*stepSize; |
202 | parMgr->setAmpParameter( parScan, value ); |
203 | |
204 | cfgInfo->setFitName(fitName + "_scan"); |
205 | |
206 | if(useMinos) |
207 | fitManager->minosMinimization(); |
208 | else |
209 | fitManager->migradMinimization(); |
210 | |
211 | bool fitFailed = (fitManager->status() != 0 && fitManager->eMatrixStatus() != 3); |
212 | |
213 | if( fitFailed ) |
214 | cout << "ERROR: fit failed use results with caution..." << endl; |
215 | |
216 | cout << "LIKELIHOOD AFTER MINIMIZATION: " << ati.likelihood() << endl; |
217 | |
218 | ati.finalizeFit(to_string(i)); |
219 | |
220 | if( seedfile.size() != 0 && !fitFailed ){ |
221 | string seedfile_scan = seedfile + Form("_scan_%d.txt", i); |
222 | ati.fitResults()->writeSeed( seedfile_scan ); |
223 | } |
224 | } |
225 | } |
226 | |
227 | int main( int argc, char* argv[] ){ |
228 | |
229 | |
230 | |
231 | bool useMinos = false; |
232 | |
233 | string configfile; |
234 | string seedfile; |
235 | string scanPar; |
236 | int numRnd = 0; |
237 | int maxIter = 10000; |
238 | |
239 | |
240 | |
241 | for (int i = 1; i < argc; i++){ |
| 1 | Assuming 'i' is >= 'argc' | |
|
| 2 | | Loop condition is false. Execution continues on line 272 | |
|
242 | |
243 | string arg(argv[i]); |
244 | |
245 | if (arg == "-c"){ |
246 | if ((i+1 == argc) || (argv[i+1][0] == '-')) arg = "-h"; |
247 | else configfile = argv[++i]; } |
248 | if (arg == "-s"){ |
249 | if ((i+1 == argc) || (argv[i+1][0] == '-')) arg = "-h"; |
250 | else seedfile = argv[++i]; } |
251 | if (arg == "-r"){ |
252 | if ((i+1 == argc) || (argv[i+1][0] == '-')) arg = "-h"; |
253 | else numRnd = atoi(argv[++i]); } |
254 | if (arg == "-m"){ |
255 | if ((i+1 == argc) || (argv[i+1][0] == '-')) arg = "-h"; |
256 | else maxIter = atoi(argv[++i]); } |
257 | if (arg == "-n") useMinos = true; |
258 | if (arg == "-p"){ |
259 | if ((i+1 == argc) || (argv[i+1][0] == '-')) arg = "-h"; |
260 | else scanPar = argv[++i]; } |
261 | if (arg == "-h"){ |
262 | cout << endl << " Usage for: " << argv[0] << endl << endl; |
263 | cout << " -n \t\t\t\t\t use MINOS instead of MIGRAD" << endl; |
264 | cout << " -c <file>\t\t\t\t config file" << endl; |
265 | cout << " -s <output file>\t\t\t for seeding next fit based on this fit (optional)" << endl; |
266 | cout << " -r <int>\t\t\t Perform <int> fits each seeded with random parameters" << endl; |
267 | cout << " -p <parameter> \t\t\t\t Perform a scan of given parameter. Stepsize, min, max are to be set in cfg file" << endl; |
268 | cout << " -m <int>\t\t\t Maximum number of fit iterations" << endl; |
269 | exit(1);} |
270 | } |
271 | |
272 | if (configfile.size() == 0){ |
| |
273 | cout << "No config file specified" << endl; |
274 | exit(1); |
275 | } |
276 | |
277 | ConfigFileParser parser(configfile); |
278 | ConfigurationInfo* cfgInfo = parser.getConfigurationInfo(); |
279 | cfgInfo->display(); |
280 | |
281 | AmpToolsInterface::registerAmplitude( BreitWigner() ); |
282 | AmpToolsInterface::registerAmplitude( BreitWigner3body() ); |
283 | AmpToolsInterface::registerAmplitude( TwoPSAngles() ); |
284 | AmpToolsInterface::registerAmplitude( TwoPSHelicity() ); |
285 | AmpToolsInterface::registerAmplitude( TwoPiAngles() ); |
286 | AmpToolsInterface::registerAmplitude( TwoPiAngles_amp() ); |
287 | AmpToolsInterface::registerAmplitude( TwoPiAngles_primakoff() ); |
288 | AmpToolsInterface::registerAmplitude( TwoPiWt_primakoff() ); |
289 | AmpToolsInterface::registerAmplitude( TwoPiWt_sigma() ); |
290 | AmpToolsInterface::registerAmplitude( TwoPitdist() ); |
291 | AmpToolsInterface::registerAmplitude( ThreePiAngles() ); |
292 | AmpToolsInterface::registerAmplitude( ThreePiAnglesSchilling() ); |
293 | AmpToolsInterface::registerAmplitude( TwoPiAnglesRadiative() ); |
294 | AmpToolsInterface::registerAmplitude( Zlm() ); |
295 | AmpToolsInterface::registerAmplitude( b1piAngAmp() ); |
296 | AmpToolsInterface::registerAmplitude( omegapiAngAmp() ); |
297 | AmpToolsInterface::registerAmplitude( polCoef() ); |
298 | AmpToolsInterface::registerAmplitude( Uniform() ); |
299 | AmpToolsInterface::registerAmplitude( dblRegge() ); |
300 | AmpToolsInterface::registerAmplitude( omegapi_amplitude() ); |
301 | AmpToolsInterface::registerAmplitude( Vec_ps_refl() ); |
302 | AmpToolsInterface::registerAmplitude( Piecewise() ); |
303 | |
304 | AmpToolsInterface::registerDataReader( ROOTDataReader() ); |
305 | AmpToolsInterface::registerDataReader( ROOTDataReaderBootstrap() ); |
306 | AmpToolsInterface::registerDataReader( ROOTDataReaderWithTCut() ); |
307 | AmpToolsInterface::registerDataReader( ROOTDataReaderTEM() ); |
308 | |
309 | if(numRnd==0){ |
| |
310 | if(scanPar=="") |
| |
311 | runSingleFit(cfgInfo, useMinos, maxIter, seedfile); |
312 | else |
313 | runParScan(cfgInfo, useMinos, maxIter, seedfile, scanPar); |
| |
314 | } else { |
315 | runRndFits(cfgInfo, useMinos, maxIter, seedfile, numRnd, 0.5); |
316 | } |
317 | |
318 | return 0; |
319 | } |
320 | |
321 | |