source: trunk/MagicSoft/Mars/sinope.cc@ 8761

Last change on this file since 8761 was 8106, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 16.0 KB
Line 
1/* ======================================================================== *\
2! $Name: not supported by cvs2svn $:$Id: sinope.cc,v 1.12 2006-10-17 17:15:58 tbretz Exp $
3! --------------------------------------------------------------------------
4!
5! *
6! * This file is part of MARS, the MAGIC Analysis and Reconstruction
7! * Software. It is distributed to you in the hope that it can be a useful
8! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
9! * It is distributed WITHOUT ANY WARRANTY.
10! *
11! * Permission to use, copy, modify and distribute this software and its
12! * documentation for any purpose is hereby granted without fee,
13! * provided that the above copyright notice appear in all copies and
14! * that both that copyright notice and this permission notice appear
15! * in supporting documentation. It is provided "as is" without express
16! * or implied warranty.
17! *
18!
19!
20! Author(s): Thomas Bretz
21! Author(s): Daniela Doener
22!
23! Copyright: MAGIC Software Development, 2000-2006
24!
25!
26\* ======================================================================== */
27#include <errno.h>
28#include <fstream>
29
30#include <TROOT.h>
31#include <TApplication.h>
32#include <TObjectTable.h>
33
34#include <TH1.h>
35#include <TLine.h>
36#include <TString.h>
37#include <TVirtualPad.h>
38
39#include "MRawRunHeader.h"
40#include "MRawEvtData.h"
41#include "MRawEvtPixelIter.h"
42#include "MRawFileRead.h"
43#include "MArrayD.h"
44#include "MFTriggerPattern.h"
45#include "MTriggerPatternDecode.h"
46#include "MContinue.h"
47#include "MTaskInteractive.h"
48#include "MLog.h"
49#include "MLogManip.h"
50#include "MArgs.h"
51#include "MSequence.h"
52#include "MStatusDisplay.h"
53#include "MParList.h"
54#include "MTaskList.h"
55#include "MEvtLoop.h"
56#include "MRawFileRead.h"
57#include "MDirIter.h"
58
59using namespace std;
60
61MRawEvtData data;
62MRawRunHeader header;
63
64MArrayD artime;
65MArrayD height(256);
66Int_t entries=0;
67Int_t events =0;
68
69Int_t PreProcess(MParList *plist)
70{
71 artime.Set(header.GetNumSamplesHiGain() + header.GetNumSamplesLoGain());
72 return kTRUE;
73}
74
75Int_t Process()
76{
77 events++;
78
79 // This is the workaround to put hi- and lo-gains together
80 const Int_t nhigain = header.GetNumSamplesHiGain();
81 const Int_t nlogain = header.GetNumSamplesLoGain();
82
83 const Int_t n = nhigain+nlogain;
84
85 // Real Process
86 MRawEvtPixelIter pixel(&data);
87
88 Byte_t slices[n];
89
90 while (pixel.Next())
91 {
92 // This is the fast workaround to put hi- and lo-gains together
93 memcpy(slices, pixel.GetHiGainSamples(), nhigain);
94 memcpy(slices+nhigain, pixel.GetLoGainSamples(), nlogain);
95
96 Byte_t *max = slices;
97 Byte_t *min = slices;
98
99 for (Byte_t *b=slices+1; b<slices+n; b++)
100 {
101 if (*b>=*max)
102 max = b;
103 if (*b<=*min)
104 min = b;
105 }
106
107 const Int_t smax = max-slices;
108 const Int_t diff = *max-*min;
109
110 if (diff<50 || diff>235) // no-signal
111 continue;
112
113 height[diff]++;
114 artime[smax]++;
115 entries++;
116 }
117 return kTRUE;
118}
119
120Int_t GetFWHM(const TH1D &h, Int_t &min, Int_t &max)
121{
122 const Double_t hmax = h.GetMaximum()/2;
123 const Int_t bin = h.GetMaximumBin();
124
125 for (min=bin; min>1; min--)
126 if (h.GetBinContent(min)<hmax)
127 break;
128
129 for (max=bin; max<h.GetNbinsX(); max++)
130 if (h.GetBinContent(max)<hmax)
131 break;
132
133 return max-min;
134}
135
136TString kOutpath="";
137TString kOutfile="";
138MStatusDisplay *d=0;
139
140Int_t PostProcess()
141{
142 if (entries==0)
143 {
144 gLog << warn << "No entries processed..." << endl;
145
146 ofstream fout(Form("%stxt", kOutpath.Data()));
147 if (!fout)
148 {
149 gLog << err << "Cannot open file: " << strerror(errno) << endl;
150 return kERROR;
151 }
152
153 fout << "Events: " << events << endl;
154 fout << endl;
155
156 return kTRUE;
157 }
158
159 TH1D h1("Arrival", "Arrival Time distribution for signals", artime.GetSize(), -0.5, artime.GetSize()-0.5);
160 TH1D h2("Height", "Pulse height distribution", height.GetSize(), -0.5, height.GetSize()-0.5);
161 h1.SetXTitle("Arrival Time [slice]");
162 h2.SetXTitle("Pulse Height [cts]");
163 h1.SetDirectory(0);
164 h2.SetDirectory(0);
165 h1.SetEntries(entries);
166 h2.SetEntries(entries);
167 memcpy(h1.GetArray()+1, artime.GetArray(), artime.GetSize()*sizeof(Double_t));
168 memcpy(h2.GetArray()+1, height.GetArray(), height.GetSize()*sizeof(Double_t));
169
170 TLine l;
171 l.SetLineColor(kGreen);
172 l.SetLineWidth(2);
173
174 Int_t min, max;
175
176 //MStatusDisplay *d = new MStatusDisplay;
177 d->AddTab("Time");
178 h1.DrawCopy();
179 l.DrawLine(h1.GetMaximumBin()-1, 0, h1.GetMaximumBin()-1, h1.GetMaximum());
180 const Int_t fwhm1 = GetFWHM(h1, min, max);
181 const Bool_t asym1 = TMath::Abs((h1.GetMaximumBin()-min)-(max-h1.GetMaximumBin()))>fwhm1/2;;
182 l.DrawLine(min-1, h1.GetMaximum()/2, max-1, h1.GetMaximum()/2);
183 gPad->Update();
184
185 d->AddTab("Pulse");
186 h2.DrawCopy();
187 l.DrawLine(h2.GetMaximumBin()-1, 0, h2.GetMaximumBin()-1, h2.GetMaximum());
188 const Int_t fwhm2 = GetFWHM(h2, min, max);
189 const Bool_t asym2 = TMath::Abs((h2.GetMaximumBin()-min)-(max-h2.GetMaximumBin()))>fwhm2/2;;
190 l.DrawLine(min-1, h2.GetMaximum()/2, max-1, h2.GetMaximum()/2);
191 gPad->Update();
192
193 d->SaveAsRoot(Form("%sroot", kOutpath.Data()));
194
195 ofstream fout(Form("%stxt", kOutpath.Data()));
196 if (!fout)
197 {
198 gLog << err << "Cannot open file: " << strerror(errno) << endl;
199 return kERROR;
200 }
201
202 fout << "Events: " << events << endl;
203 fout << "HasSignal: " << (fwhm1>10?"no":"yes") << endl;
204 fout << "HasPedestal: " << (fwhm1<20?"no":"yes") << endl;
205 fout << endl;
206 fout << "PositionSignal: " << h1.GetMaximumBin()-1 << endl;
207 fout << "PositionFWHM: " << fwhm1 << endl;
208 fout << "PositionAsym: " << (asym1?"yes":"no") << endl;
209 fout << endl;
210 fout << "HeightSignal: " << h2.GetMaximumBin()-1 << endl;
211 fout << "HeightFWHM: " << fwhm2 << endl;
212 fout << "HeightAsym: " << (asym2?"yes":"no") << endl;
213 fout << endl;
214
215 return kTRUE;
216}
217
218static void StartUpMessage()
219{
220 gLog << all << endl;
221
222 // 1 2 3 4 5
223 // 12345678901234567890123456789012345678901234567890
224 gLog << "==================================================" << endl;
225 gLog << " Sinope - MARS V" << MARSVER << endl;
226 gLog << " MARS -- SImple Non Online Pulse Evaluation" << endl;
227 gLog << " Compiled with ROOT v" << ROOT_RELEASE << " on <" << __DATE__ << ">" << endl;
228 gLog << "==================================================" << endl;
229 gLog << endl;
230}
231
232static void Usage()
233{
234 // 1 2 3 4 5 6 7 8
235 // 12345678901234567890123456789012345678901234567890123456789012345678901234567890
236 gLog << all << endl;
237 gLog << "Sorry the usage is:" << endl;
238 gLog << " sinope [options] --run={number} --date={yyyy-mm-dd}" << endl << endl;
239 gLog << " Arguments:" << endl;
240 gLog << " --run={number}: Run number of run to process" << endl;
241 gLog << " --date={yy-mm-dd}: Night the run belongs to" << endl << endl;
242 gLog << " Root Options:" << endl;
243 gLog << " -b Batch mode (no graphical output to screen)" << endl<<endl;
244 gLog << " Options:" << endl;
245 gLog.Usage();
246 // gLog << " --debug-env=0 Disable debugging setting resources <default>" << endl;
247 // gLog << " --debug-env[=1] Display untouched resources after program execution" << endl;
248 // gLog << " --debug-env=2 Display untouched resources after eventloop setup" << endl;
249 // gLog << " --debug-env=3 Debug setting resources from resource file" << endl;
250 gLog << " --debug-mem Debug memory usage" << endl << endl;
251 gLog << endl;
252 gLog << " -q Quit when job is finished" << endl;
253 gLog << " -f Force overwrite of existing files" << endl;
254 gLog << " -dat Run sinope only for data events in the run" << endl;
255 gLog << " -cal Run sinope only for calibration events in the run" << endl;
256 gLog << " --ind=path Path where to search for the data file" << endl;
257 gLog << " [default=standard path in datacenter]" << endl;
258 gLog << " --out=path Path to write the all results to [def=local path]" << endl;
259 gLog << " --outf=filename Filename of the outputfiles [def=sinope{runnumber}.*]" << endl;
260 gLog << " --num={number} Number of events to process (default=1000)" << endl;
261 gLog << " --print-seq Print Sequence information" << endl;
262 gLog << " --print-files Print Files taken from Sequence" << endl;
263 gLog << " --print-found Print Files found from Sequence" << endl;
264 // gLog << " --config=callisto.rc Resource file [default=callisto.rc]" << endl;
265 gLog << endl;
266 gLog << " --version, -V Show startup message with version number" << endl;
267 gLog << " -?, -h, --help This help" << endl << endl;
268 gLog << "Background:" << endl;
269 gLog << " Sinope is Jupiter's sixteenth moon. Sinope is 28km in diameter and" << endl;
270 gLog << " and orbits 23,700,000km from Jupiter. Sinope has a mass of 8e16kg." << endl;
271 gLog << " It orbits Jupiter in 758days and is in a retrograde orbit (orbiting" << endl;
272 gLog << " opposite to the direction of Jupiter). Very little is known about" << endl;
273 gLog << " Sinope. Sinope was discovered by S.Nicholson in 1914." << endl << endl;
274 gLog << "Example:" << endl;
275 gLog << " sinope -f --date=2004-05-06 --run=32456" << endl;
276 gLog << endl;
277}
278
279static void PrintFiles(const MSequence &seq, const TString &kInpathD, Bool_t ball)
280{
281 const char *prep = ball ? "Found" : "Scheduled";
282
283 MDirIter Next;
284 seq.SetupAllRuns(Next, kInpathD, kTRUE);
285
286 gLog << all;
287 gLog.Separator(Form("%s Files", prep));
288 Next.Print(ball?"all":"");
289 gLog << endl;
290}
291
292int main(int argc, char **argv)
293{
294 if (!MARS::CheckRootVer())
295 return 0xff;
296
297 MLog::RedirectErrorHandler(MLog::kColor);
298
299 //
300 // Evaluate arguments
301 //
302 MArgs arg(argc, argv, kTRUE);
303 gLog.Setup(arg);
304
305 StartUpMessage();
306
307 if (arg.HasOnly("-V") || arg.HasOnly("--version"))
308 return 0;
309
310 if (arg.HasOnly("-?") || arg.HasOnly("-h") || arg.HasOnly("--help"))
311 {
312 Usage();
313 return 2;
314 }
315
316 //const TString kConfig = arg.GetStringAndRemove("--config=", "callisto.rc");
317
318 const Bool_t kPrintSeq = arg.HasOnlyAndRemove("--print-seq");
319 const Bool_t kPrintFiles = arg.HasOnlyAndRemove("--print-files");
320 const Bool_t kPrintFound = arg.HasOnlyAndRemove("--print-found");
321 const Bool_t kDebugMem = arg.HasOnlyAndRemove("--debug-mem");
322 //Int_t kDebugEnv = arg.HasOnlyAndRemove("--debug-env") ? 1 : 0;
323 //kDebugEnv = arg.GetIntAndRemove("--debug-env=", kDebugEnv);
324
325 const Bool_t kQuit = arg.HasOnlyAndRemove("-q");
326 const Bool_t kBatch = arg.HasOnlyAndRemove("-b");
327 const Bool_t kOverwrite = arg.HasOnlyAndRemove("-f");
328 const Bool_t kData = arg.HasOnlyAndRemove("-dat");
329 const Bool_t kCal = arg.HasOnlyAndRemove("-cal");
330 //const Bool_t kForceExec = arg.HasOnlyAndRemove("-ff");
331
332 const TString kInpathD = arg.GetStringAndRemove("--ind=", "");
333 /*const TString*/ kOutpath = arg.GetStringAndRemove("--out=", "");
334 /*const TString*/ kOutfile = arg.GetStringAndRemove("--outf=", "");
335
336 const Int_t kNumEvents = arg.GetIntAndRemove("--num=", header.GetNumEvents());
337 const Int_t kRunNumber = arg.GetIntAndRemove("--run=", -1);
338 const TString kDate = arg.GetStringAndRemove("--date=", "");
339
340 if (arg.GetNumOptions()>0)
341 {
342 gLog << warn << "WARNING - Unknown commandline options..." << endl;
343 arg.Print("options");
344 gLog << endl;
345 Usage();
346 return 2;
347 }
348
349 if (kRunNumber<0)
350 {
351 gLog << warn << "ERROR - No '--run=' option given... required." << endl;
352 gLog << endl;
353 Usage();
354 return 2;
355 }
356 if (kDate.IsNull())
357 {
358 gLog << warn << "ERROR - No '--date=' option given... required." << endl;
359 gLog << endl;
360 Usage();
361 return 2;
362 }
363
364 //
365 // check for the right usage of the program
366 //
367 if (arg.GetNumArguments()>0)
368 {
369 Usage();
370 return 2;
371 }
372
373 if (kDebugMem)
374 TObject::SetObjectStat(kTRUE);
375
376 //
377 // Setup sequence and check its validity
378 //
379 MSequence seq;
380 seq.SetNight(kDate);
381 seq.AddRuns(kRunNumber);
382 if (kPrintSeq)
383 {
384 gLog << all;
385 gLog.Separator();
386 seq.Print();
387 gLog << endl;
388 }
389 if (!seq.IsValid())
390 {
391 gLog << err << "Sequence invalid!" << endl << endl;
392 return 2;
393 }
394
395 //
396 // Process print options
397 //
398 if (kPrintFiles)
399 PrintFiles(seq, kInpathD, kFALSE);
400 if (kPrintFound)
401 PrintFiles(seq, kInpathD, kTRUE);
402
403 if (kOutpath.IsNull())
404 {
405 kOutpath = seq.GetStandardPath();
406 kOutpath += "rawfiles/";
407 }
408 if (!kOutpath.EndsWith("/"))
409 kOutpath += "/";
410 if (kOutfile.IsNull())
411 kOutpath += Form("sinope%08d.", kRunNumber);
412 else
413 kOutpath += kOutfile+".";
414
415 if (!kOverwrite)
416 {
417 TString file = Form("%sroot", kOutpath.Data());
418 if (!gSystem->AccessPathName(file, kFileExists))
419 {
420 gLog << err << "Sorry, file '" << file << "' exists... use -f option.." << endl;
421 return 2;
422 }
423 file = Form("%stxt", kOutpath.Data());
424 if (!gSystem->AccessPathName(file, kFileExists))
425 {
426 gLog << err << "Sorry, file '" << file << "' exists... use -f option.." << endl;
427 return 2;
428 }
429 }
430
431 //
432 // Initialize root
433 //
434 MArray::Class()->IgnoreTObjectStreamer();
435 MParContainer::Class()->IgnoreTObjectStreamer();
436
437 TApplication app("sinope", &argc, argv);
438 if (!gROOT->IsBatch() && !gClient || gROOT->IsBatch() && !kBatch)
439 {
440 gLog << err << "Bombing... maybe your DISPLAY variable is not set correctly!" << endl;
441 return 1;
442 }
443
444 // ----------------------------------------------------------
445
446 /*MStatusDisplay **/d = new MStatusDisplay;
447
448 // From now on each 'Exit' means: Terminate the application
449 d->SetBit(MStatusDisplay::kExitLoopOnExit);
450 d->SetTitle(Form("Sinope #%d", kRunNumber));
451
452 MDirIter iter;
453 seq.SetupAllRuns(iter, 0, kTRUE);
454
455 MRawFileRead read;
456 read.AddFiles(iter);
457
458 MTaskInteractive task;
459
460 MTriggerPatternDecode decode;
461 MFTriggerPattern ftp;
462 if (kCal || kData)
463 {
464 ftp.SetDefault(kTRUE);
465 ftp.DenyCalibration();
466 if (!kCal)
467 {
468 ftp.DenyPedestal();
469 ftp.SetInverted();
470 }
471 }
472 MContinue conttp(&ftp, "ContTrigPattern");
473
474 task.SetPreProcess(PreProcess);
475 task.SetProcess(Process);
476 task.SetPostProcess(PostProcess);
477
478 MTaskList tlist;
479 tlist.AddToList(&read);
480 if (kCal || kData)
481 {
482 tlist.AddToList(&decode);
483 tlist.AddToList(&conttp);
484 }
485 tlist.AddToList(&task);
486
487 MParList plist;
488 plist.AddToList(&tlist);
489
490 plist.AddToList(&data);
491 plist.AddToList(&header);
492
493 MEvtLoop evtloop;
494 evtloop.SetParList(&plist);
495 evtloop.SetDisplay(d);
496
497 if (!evtloop.Eventloop(kNumEvents))
498 return 1;
499
500 if (!evtloop.GetDisplay())
501 {
502 gLog << warn << "Display closed by user... execution aborted." << endl << endl;
503 return 0;
504 }
505
506 if (kBatch || kQuit)
507 delete d;
508 else
509 {
510 // From now on each 'Close' means: Terminate the application
511 d->SetBit(MStatusDisplay::kExitLoopOnClose);
512
513 // Wait until the user decides to exit the application
514 app.Run(kFALSE);
515 }
516
517 if (TObject::GetObjectStat())
518 {
519 TObject::SetObjectStat(kFALSE);
520 gObjectTable->Print();
521 }
522
523 return 0;
524}
Note: See TracBrowser for help on using the repository browser.