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

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