source: trunk/FACT++/src/calcsourcemc.cc@ 19878

Last change on this file since 19878 was 19432, checked in by tbretz, 6 years ago
Added the calcsource for MC data:
File size: 14.8 KB
Line 
1#include "Database.h"
2
3#include "pal.h"
4#include "nova.h"
5#include "tools.h"
6#include "Time.h"
7#include "Configuration.h"
8
9#include <TROOT.h>
10#include <TVector3.h>
11#include <TRotation.h>
12
13using namespace std;
14
15// ------------------------------------------------------------------------
16
17void SetupConfiguration(Configuration &conf)
18{
19 po::options_description control("Calcsource options");
20 control.add_options()
21 ("uri,u", var<string>()
22#if BOOST_VERSION >= 104200
23 ->required()
24#endif
25 , "Database link as in\n\tuser:password@server[:port]/database[?compress=0|1].")
26 //("source-key", var<uint32_t>(5), "")
27 //("source-name", var<string>(""), "")
28 ("file", var<uint32_t>(uint32_t(0)),"FileId (YYMMDDXXX), defined the events to be processed (if omitted (=0), a list of possible numbers is printed)")
29 ("table.runinfo", var<string>("RunInfoMC"), "Name of the table where the run info is stored")
30 ("table.events", var<string>("EventsMC"), "Name of the table where the events are stored")
31 ("table.position", var<string>("PositionMC"), "Name of the table where the calculated posiiton will be stored")
32 ("engine", var<string>(""), "Database engine to be used when a new table is created")
33 ("row-format", var<string>(""), "Defines the ROW_FORMAT keyword for table creation")
34 ("ignore-errors", po_switch(), "Adds the IGNORE keyword to the INSERT query (turns errors into warnings, ignores rows with errors)")
35 ("force", po_switch(), "Force processing even if there is no database connection")
36 ("drop", po_switch(), "Drop the table (implies create)")
37 ("create", po_switch(), "Create the table if it does not yet exist")
38 ("update", po_switch(), "Uses the ON DUPLICATE KEY to update already existing extries")
39 ("focal-dist", var<double>(4889.), "Focal distance of the camera in millimeter")
40 ;
41
42 po::options_description debug("Debug options");
43 debug.add_options()
44 ("no-insert", po_switch(), "Does not insert or update any data to any table")
45 ("dry-run", po_switch(), "Skip any query which changes the databse (might result in consecutive failures)")
46 ("print-meta", po_switch(), "Print meta-queries (DROP, CREATE, DELETE, SELECT)")
47 ("print-insert", po_switch(), "Print the INSERT/UPDATE queries")
48 ("verbose,v", var<uint16_t>(1), "Verbosity (0: quiet, 1: default, 2: more, 3, ...)")
49 ;
50
51 po::positional_options_description p;
52 p.add("file", 1); // The 1st positional options (n=1)
53
54 conf.AddOptions(control);
55 conf.AddOptions(debug);
56 conf.SetArgumentPositions(p);
57}
58
59void PrintUsage()
60{
61 cout <<
62 "calcsourcemc - Fill a MC table with MC source positions\n"
63 "\n"
64 "This tool is to calculate the source position in the camera and fill "
65 "it into a database table for all events which come from a single run. "
66 "The run is idetified by its FileId."
67 "\n\n"
68 "Then for each event with the given FileId, the telescope pointing "
69 "position (`Zd`, `Az`) and the primary particle direction "
70 "(`ParticleTheta`,`ParticlePhi`) and its (`EvtNumber`) are requested "
71 "from the table given by --table.events. "
72 "Based on this information and the focal distance (--focal-distance) the "
73 "`X` and `Y` coordinate of the source in the camera plane is calculated. "
74 "The result can then be filled into a database."
75 "\n\n"
76 "The table to be filled or updated should contain the following columns:\n"
77 " - FileId INT UNSIGNED NOT NULL\n"
78 " - EvtNumber INT UNSIGNED NOT NULL\n"
79 " - X FLOAT NOT NULL\n"
80 " - Y FLOAT NOT NULL\n"
81 "\n"
82 "If the table does not exist, it can be created automatically specifying "
83 "the --crate option. If a table already exists and it should be dropped "
84 "before creation, the --drop option can be used:\n"
85 "\n"
86 " calcsource 170909009 --create\n"
87 "\n"
88 "Process the data from the run 009 from 09/09/2017. If no table exists "
89 "a table with the name given by --table.position is created.\n"
90 "\n"
91 " calcsource 170909009 --create --drop\n"
92 "\n"
93 "Same as before, but if a table with the name given by --table.position "
94 "exists, it is dropped before.\n"
95 "\n"
96 "For each event, a new row is inserted. If existing rows should be updated, "
97 "use:\n"
98 "\n"
99 " calcsource 170909009 --updated\n"
100 "\n"
101 "The --create option is compatible with that. The --drop option is ignored.\n"
102 "\n"
103 "To avoid failure in case an entry does already exist, you can add the IGNORE "
104 "keyword to the INSERT query by --ignore-errors, which essentially ignores "
105 "all errors and turns them into warnings which are printed after the query "
106 "succeeded.\n"
107 "\n"
108 "\n"
109 "For debugging purposes several print options and options to avoid irreversible "
110 "changes to the database exist."
111 "\n\n"
112 "Usage: calcsource YYMMDDXXX [-u URI] [options]\n"
113 "\n"
114 ;
115 cout << endl;
116}
117
118/*
119class MRotation : public TRotation
120{
121public:
122 MRotation() : TRotation(1, 0, 0, 0, -1, 0, 0, 0, 1)
123 {
124 }
125};*/
126
127
128int main(int argc, const char* argv[])
129{
130 Time start;
131
132 gROOT->SetBatch();
133
134 Configuration conf(argv[0]);
135 conf.SetPrintUsage(PrintUsage);
136 SetupConfiguration(conf);
137
138 if (!conf.DoParse(argc, argv))
139 return 127;
140
141 // ----------------------------- Evaluate options --------------------------
142 const string uri = conf.Get<string>("uri");
143 const string tab_runinfo = conf.Get<string>("table.runinfo");
144 const string tab_events = conf.Get<string>("table.events");
145 const string tab_position = conf.Get<string>("table.position");
146
147 const string engine = conf.Get<string>("engine");
148 const string row_format = conf.Get<string>("row-format");
149 const bool ignore_errors = conf.Get<bool>("ignore-errors");
150
151 const uint32_t file = conf.Get<uint32_t>("file");
152
153 const double focal_dist = conf.Get<double>("focal-dist");
154
155 const bool print_meta = conf.Get<bool>("print-meta");
156 const bool print_insert = conf.Get<bool>("print-insert");
157
158 const bool force = conf.Get<bool>("force");
159 const bool drop = conf.Get<bool>("drop");
160 const bool create = conf.Get<bool>("create") || drop;
161 const bool update = conf.Get<bool>("update");
162 const bool noinsert = conf.Get<bool>("no-insert");
163 const bool dry_run = conf.Get<bool>("dry-run");
164 const uint16_t verbose = conf.Get<uint16_t>("verbose");
165
166 // -------------------------------------------------------------------------
167 // Checking for database connection
168
169 Database connection(uri); // Keep alive while fetching rows
170
171 try
172 {
173 if (!force)
174 connection.connected();
175 }
176 catch (const exception &e)
177 {
178 cerr << "SQL connection failed: " << e.what() << endl;
179 return 1;
180 }
181
182 // -------------------------------------------------------------------------
183 // list file-ids in the events table
184
185 if (file==0)
186 {
187 const string query =
188 "SELECT FileId FROM `"+tab_runinfo+"`";
189
190 if (print_meta)
191 cout << query << endl;
192
193 const mysqlpp::StoreQueryResult res =
194 connection.query(query).store();
195
196 for (size_t i=0; i<res.num_rows(); i++)
197 cout << "calcsourcemc " << res[i][0] << '\n';
198 cout << endl;
199
200 return 0;
201 }
202
203 // -------------------------------------------------------------------------
204 // create INSERT/UPDATE query (calculate positions)
205
206 if (verbose>0)
207 {
208 cout << "\n-------------------------- Evaluating file ------------------------";
209 cout << "\nRequesting data from table `" << tab_events << "`" << endl;
210 }
211
212 const string query =
213 "SELECT ParticleTheta, ParticlePhi, Zd, Az, EvtNumber"
214 " FROM `"+tab_events+"`"
215 " WHERE FileId="+to_string(file);
216
217 if (print_meta)
218 cout << query << endl;
219
220 const mysqlpp::UseQueryResult res1 =
221 connection.query(query).use();
222
223 ostringstream ins;
224 ins << setprecision(16);
225
226 size_t count = 0;
227 while (auto row=res1.fetch_row())
228 {
229 count++;
230
231 const double part_zd = row[0];
232 const double part_az = row[1];
233 const double tel_zd = row[2]*M_PI/180;
234 const double tel_az = row[3]*M_PI/180;
235 const uint32_t event = row[4];
236
237 /*
238 // ============================ Mars ================================
239
240 MVector3 pos0, pos;
241 pos0.SetZdAz(fPointPos->GetZdRad(), fPointPos->GetAzRad());
242 pos.SetZdAz( fMcEvt->GetParticleTheta(), fMcEvt->GetParticlePhi());
243
244 srcpos = MAstro::GetDistOnPlain(pos0, pos, -fGeom->GetCameraDist()*1000);
245 */
246
247 // ================================= Nova ===============================
248
249 TVector3 pos;
250 TVector3 pos0;
251 pos.SetMagThetaPhi( 1, part_zd, part_az);
252 pos0.SetMagThetaPhi(1, tel_zd, tel_az);
253
254 pos.RotateZ(-pos0.Phi());
255 pos.RotateY(-pos0.Theta());
256 pos.RotateZ(-M_PI/2); // exchange x and y
257 pos *= -focal_dist/pos.Z();
258
259 TVector2 v = pos.XYvector();
260
261 ins << "( " << file << ", " << event << ", " << v.X() << ", " << v.Y() << " ),\n";
262 }
263
264 if (connection.errnum())
265 {
266 cerr << "SQL error fetching row: " << connection.error() << endl;
267 return 4;
268 }
269
270 if (verbose>0)
271 cout << "Processed " << count << " events.\n" << endl;
272
273 if (count==0)
274 {
275 if (verbose>0)
276 cout << "Total execution time: " << Time().UnixTime()-start.UnixTime() << "s\n" << endl;
277 return 0;
278 }
279
280 // -------------------------------------------------------------------------
281 // drop table if requested
282
283 if (drop)
284 {
285 try
286 {
287 if (verbose>0)
288 cout << "Dropping table `" << tab_position << "`" << endl;
289
290 if (!dry_run)
291 connection.query("DROP TABLE `"+tab_position+"`").execute();
292
293 if (verbose>0)
294 cout << "Table `" << tab_position << "` dropped.\n" << endl;
295 }
296 catch (const exception &e)
297 {
298 cerr << "DROP TABLE `" << tab_position << "`\n\n";
299 cerr << "SQL query failed: " << e.what() << endl;
300 return 5;
301 }
302 }
303
304 // -------------------------------------------------------------------------
305 // crate table if requested
306
307 if (create)
308 {
309 if (verbose>0)
310 cout << "Creating table `" << tab_position << "`" << endl;
311
312 string query2 =
313 "CREATE TABLE IF NOT EXISTS `"+tab_position+"`\n"
314 "(\n"
315 " FileId INT UNSIGNED NOT NULL,\n"
316 " EvtNumber INT UNSIGNED NOT NULL,\n"
317 " X FLOAT NOT NULL,\n"
318 " Y FLOAT NOT NULL,\n"
319 " PRIMARY KEY (FileId, EvtNumber)\n"
320 ")\n"
321 "DEFAULT CHARSET=latin1 COLLATE=latin1_general_ci\n";
322 if (!engine.empty())
323 query2 += "ENGINE="+engine+"\n";
324 if (!row_format.empty())
325 query2 += "ROW_FORMAT="+row_format+"\n";
326 query2 += "COMMENT='created by "+conf.GetName()+"'\n";
327
328 try
329 {
330 if (!dry_run)
331 connection.query(query2).execute();
332 }
333 catch (const exception &e)
334 {
335 cerr << query2 << "\n\n";
336 cerr << "SQL query failed: " << e.what() << endl;
337 return 6;
338 }
339
340 if (print_meta)
341 cout << query2 << endl;
342
343 if (verbose>0)
344 cout << "Table `" << tab_position << "` created.\n" << endl;
345 }
346
347 // -------------------------------------------------------------------------
348 // delete old entries from table
349
350 if (!drop)
351 {
352 if (verbose>0)
353 cout << "Deleting old entries from table `" << tab_position << "`" << endl;
354
355 const string query2 =
356 "DELETE FROM `"+tab_position+"` WHERE FileId="+to_string(file);
357
358 try
359 {
360 if (!dry_run)
361 {
362 const mysqlpp::SimpleResult res =
363 connection.query(query2).execute();
364
365 if (verbose>0)
366 cout << res.rows() << " row(s) affected.\n" << endl;
367 }
368 }
369 catch (const exception &e)
370 {
371 cerr << query2 << "\n\n";
372 cerr << "SQL query failed: " << e.what() << endl;
373 return 7;
374 }
375
376 if (print_meta)
377 cout << query2 << endl;
378 }
379
380 // -------------------------------------------------------------------------
381 // insert data into table
382
383 if (verbose>0)
384 cout << "Inserting data into table " << tab_position << "." << endl;
385
386 string query2 = "INSERT ";
387 if (ignore_errors)
388 query2 += "IGNORE ";
389 query2 += "`"+tab_position+"` (FileId, EvtNumber, X, Y) VALUES\n"+
390 ins.str().substr(0, ins.str().size()-2)+
391 "\n";
392 if (update)
393 query2 += "ON DUPLICATE KEY UPDATE X=VALUES(X), Y=VALUES(Y)\n";
394
395 try
396 {
397 if (!noinsert)
398 {
399 const mysqlpp::SimpleResult res =
400 connection.query(query2).execute();
401
402 if (verbose>0)
403 cout << res.info() << '\n' << endl;
404 }
405 }
406 catch (const exception &e)
407 {
408 cerr << query2 << "\n\n";
409 cerr << "SQL query (" << query2.length() << " bytes) failed:\n" << e.what() << endl;
410 return 8;
411 }
412
413 if (print_insert)
414 cout << query2 << endl;
415
416 if (verbose>0)
417 {
418 const auto sec = Time().UnixTime()-start.UnixTime();
419 cout << "Total execution time: " << sec << "s ";
420 cout << "(" << Tools::Fractional(sec/count) << "s/row)\n";
421
422 try
423 {
424 const auto resw =
425 connection.query("SHOW WARNINGS").store();
426
427 for (size_t i=0; i<resw.num_rows(); i++)
428 {
429 const mysqlpp::Row &roww = resw[i];
430
431 cout << roww["Level"] << '[' << roww["Code"] << "]: ";
432 cout << roww["Message"] << '\n';
433 }
434
435 cout << endl;
436 }
437 catch (const exception &e)
438 {
439 cerr << "\nSHOW WARNINGS\n\n";
440 cerr << "SQL query failed:\n" << e.what() << '\n' <<endl;
441 return 8;
442 }
443 }
444
445 return 0;
446}
Note: See TracBrowser for help on using the repository browser.