source: trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.cc@ 1736

Last change on this file since 1736 was 1736, checked in by rwagner, 22 years ago
*** empty log message ***
File size: 32.1 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Thomas Bretz 11/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2003
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MCT1ReadPreProc
28//
29// Reads a output file of the CT1 preproc.
30//
31// Implements usage of a selector (see MRead) Use such a filter to skip
32// events before reading! But never use a filter which needs read data
33// as input...
34//
35// Input Containers:
36// -/-
37//
38// Output Containers:
39// MCerPhotEvt the data container for all data.
40// MPedestalCam ct1 pedestals
41// MMcEvt monte carlo data container for MC files
42// MMcTrig mc data container for trigger information
43// MSrcPosCam source position in the camera
44// MBlindPixels Array holding blind pixels
45//
46/////////////////////////////////////////////////////////////////////////////
47#include "MCT1ReadPreProc.h"
48
49#include <fstream.h>
50
51#include <TList.h>
52#include <TSystem.h>
53
54#define LINUX
55#define HISTO void
56#define HBOOK_FILE int
57#include "defines.h"
58#include "structures.h"
59
60#include "MLog.h"
61#include "MLogManip.h"
62
63#include "MTime.h"
64#include "MFilter.h"
65
66#include "MParList.h"
67#include "MCerPhotEvt.h"
68
69#include "MPedestalPix.h"
70#include "MPedestalCam.h"
71
72#include "MGeomCam.h"
73#include "MSrcPosCam.h"
74#include "MBlindPixels.h"
75
76#include "MRawRunHeader.h"
77#include "MTaskList.h"
78
79#include "MMcEvt.hxx"
80#include "MMcTrig.hxx"
81
82ClassImp(MCT1ReadPreProc);
83
84// --------------------------------------------------------------------------
85//
86// Default constructor. Creates an array which stores the file names of
87// the files which should be read. If a filename is given it is added
88// to the list.
89//
90MCT1ReadPreProc::MCT1ReadPreProc(const char *fname, const char *name,
91 const char *title) : fIn(NULL), fEntries(0)
92{
93 fName = name ? name : "MCT1ReadPreProc";
94 fTitle = title ? title : "Reads a CT1 preproc data file";
95
96 //
97 // remember file name for opening the file in the preprocessor
98 //
99 fFileNames = new TList;
100 fFileNames->SetOwner();
101
102 if (fname)
103 AddFile(fname);
104}
105
106// --------------------------------------------------------------------------
107//
108// Delete the filename list and the input stream if one exists.
109//
110MCT1ReadPreProc::~MCT1ReadPreProc()
111{
112 delete fFileNames;
113 if (fIn)
114 delete fIn;
115}
116
117// --------------------------------------------------------------------------
118//
119// Add this file as the last entry in the chain
120//
121void MCT1ReadPreProc::AddFile(const char *txt)
122{
123 const char *name = gSystem->ExpandPathName(txt);
124 TString fname(name);
125 delete [] name;
126
127 if (!CheckHeader(fname))
128 {
129 *fLog << warn << "WARNING - Problem reading header... ignored." << endl;
130 return;
131 }
132
133 const Int_t n = GetNumEvents(fname);
134 if (n==0)
135 {
136 *fLog << warn << "WARNING - File contains no data... ignored." << endl;
137 return;
138 }
139
140 fEntries += n;
141
142 *fLog << inf << "File " << txt << " contains " << n << " events (Total=" << fEntries << ")" << endl;
143
144 fFileNames->AddLast(new TNamed(txt, ""));
145}
146
147// --------------------------------------------------------------------------
148//
149// Print data from the header to the screen and analyse the header data,
150// means store and copy the needed data into Mars structures and
151// data members
152//
153void MCT1ReadPreProc::ProcessRunHeader(const struct outputpars &outpars)
154{
155 if (outpars.inumpixels != iMAXNUMPIX)
156 *fLog << warn << "WARNING! File doesn't contain " << iMAXNUMPIX << " Pixels... maybe corrupt." << endl;
157
158 fNumEventsInRun = 0;
159
160 //
161 // ------------------- Output some stuff -----------------------
162 //
163
164 // int itelescope; // number of the CT which took the data
165 *fLog << inf << "Telescope: CT" << outpars.itelescope;
166
167 // float flongitude_deg; // longitude (counted positive towards West) of CT position */
168 // float flatitude_deg; // latitude (counted positive towards North) of CT position */
169 *fLog << " located @ Longitude=" << outpars.flongitude_deg;
170 *fLog << "deg Latitude=" << outpars.flatitude_deg << "deg" << endl;
171
172 // int irunnum; // run number (from parameters file)
173 *fLog << inf << "Run number" << outpars.irunnum;
174
175 fRawRunHeader->SetRunNumber(outpars.irunnum);
176 fRawRunHeader->SetReadyToSave();
177
178 // enum onoroff {NEITHER_ON_NOR_OFF, OFF_SOURCE, ON_SOURCE} eruntype; // runtype
179 *fLog << "Run: #" << outpars.irunnum << " (";
180 switch (outpars.eruntype)
181 {
182 case NEITHER_ON_NOR_OFF: *fLog << "unknown"; break;
183 case OFF_SOURCE: *fLog << "off-source"; break;
184 case ON_SOURCE: *fLog << "on-source"; break;
185 default: *fLog << (int)outpars.eruntype; break;
186 }
187 *fLog << ", ";
188 switch (outpars.etrackmode)
189 {
190 case NORMAL: *fLog << "normal tracking"; break;
191 case REVERSE: *fLog << "reverse tracking"; break;
192 case DUNNO: *fLog << "unknown tracking"; break;
193 default: *fLog << (int)outpars.etrackmode; break;
194 }
195 *fLog << ")" << endl;
196
197 //double dsourcera_hours; // right ascension of observed source in hours
198 //double dsourcedec_deg; // declination of observed source in degrees
199 *fLog << "Source: RA=" << outpars.dsourcera_hours << "h DEC=";
200 *fLog << outpars.dsourcedec_deg << "deg" << endl;
201
202 //int inummuonpixels; // number of pixels in the muon shield
203 //int inumcointdcs; // number of coincidence tdcs recorded in the runfile
204 //float fpixdiameter_deg; // smallest pixel diameter (degrees) (from parameters file) */
205
206 // enum axes {RA, DEC, ALT, AZ} ese1_is; // name of the axis to which shaft encoder 1 is attached (implies the type of mount)
207 *fLog << "Shaftencoder 1 @ ";
208 switch (outpars.ese1_is)
209 {
210 case RA: *fLog << "RA"; break;
211 case DEC: *fLog << "DEC"; break;
212 case ALT: *fLog << "ALT"; break;
213 case AZ: *fLog << "AZ"; break;
214 default: *fLog << (int)outpars.ese1_is; break;
215 }
216 *fLog << endl;
217
218 // int isezeropos[2]; // zero position of shaftencoders 1 and 2 (from parameters file)
219 *fLog << "SE Zero: SE(1)=" << outpars.isezeropos[0] << " ";
220 *fLog << "SE(2)=" << outpars.isezeropos[1] << endl;
221
222 // int iaz_rev_track_corr; // correction for the azimuth shaft encoder (ALT/AZ mount only) in reverse tracking mode
223 // int ialt_rev_track_corr; // correction for the altitude shaft encoder (ALT/AZ mount only) in reverse tracking mode
224 *fLog << "Reverse tracking corrections: SE(az)=" << outpars.iaz_rev_track_corr;
225 *fLog << " SE(alt)=" << outpars.ialt_rev_track_corr << endl;
226
227 // float fbendingcorr; // bending correction factor (ALT/AZ mount only)
228 // float fextinction; // atmospheric extinction (typically taken from the Carlsberg Meridian Circle data)
229 *fLog << "Bending: Correction factor=" << outpars.fbendingcorr << " ";
230 *fLog << "Extinction=" << outpars.fextinction << endl;
231
232 // Boolean bdontusepix[iMAXNUMPIX]; // bdontusepix is set true if the pixel should not be used in image analysis, otherwise it is true;
233 fBlinds->Clear();
234 *fLog << "Don't use pixels: ";
235 for (int i=0; i<iMAXNUMPIX; i++)
236 if (outpars.bdontusepix[i])
237 {
238 *fLog << i << " ";
239 fBlinds->SetPixelBlind(i);
240 }
241 *fLog << endl;
242
243 *fLog << "Exclude: ";
244
245 // Boolean bexcludepix[iMAXNUMPIX];
246 for (int i=0; i<iMAXNUMPIX; i++)
247 if (outpars.bexcludepix[i])
248 {
249 *fLog << i << " ";
250 fBlinds->SetPixelBlind(i);
251 }
252 *fLog << endl;
253
254 fBlinds->SetReadyToSave();
255
256 /* bexcludepix[] is set TRUE (== exclude from pedestal, Laser
257 * calibration and the further analysis) when the Mean value
258 * of a pixel inside a pedestal Run is larger than 50 or ( || )
259 * if the pedestal RMS value of this pixel is larger than 5.0
260 * This is reflected in the (new for versions >= 0.4)
261 * variable "pixonoff" in the ntuple written by preproc:
262 * preproc.nt.hbook
263 *
264 * When the pixel is excluded by the user it will get a -2 otherwise
265 * pixonoff = 0.
266 * Additive to this a -1 is added when preproc excludes the pixel
267 * for a given Run. So the actual value tells you whether you caught
268 * it already by hand or not.
269 *
270 * A plot of pixonoff may also be handy to tell you the status of your
271 * ADC equipment. */
272
273 // float fphotoel_per_adccnt[iMAXNUMPIX]; // conversion factors for the pixel signals */
274 /*
275 float padc = outpars.fphotoel_per_adccnt[0];
276 *fLog << "Phe/ADC (pixel 0): " << padc << endl;
277 for (int i=0; i<iMAXNUMPIX; i++)
278 *fLog << outpars.fphotoel_per_adccnt[i] << " ";
279 *fLog << endl;
280 */
281 /*
282 --- USEFULL? NEEDED? ---
283 int irubminusutc_usecs; // difference between rubidium clock and UTC in microseconds
284 int isum_thresh_phot; // threshold for the total sum of photoelectrons filter
285 int i2out_thresh_phot; // threshold for the two-pixels-out-of-all software
286 int imuoncut_thresh_adccnt[iMAXNUMMUONPIX]; // thresholds for the muon cut
287 Boolean bmuon_suppression; // "Dolby" noise reduction flag
288 float ftolerated_pointerror_deg; // maximum tolerated pointing error in the position
289 */
290
291 // float fxpointcorr_deg; // pointing correction (to be added along the camera x axis) e.g. from point run */
292 // float fypointcorr_deg; // pointing correction (to be added along the camera y axis) e.g. from point run */
293 *fLog << "Pointing correction: dx=" << outpars.fxpointcorr_deg << "deg ";
294 *fLog << "dy=" << outpars.fypointcorr_deg << "deg" << endl;
295
296 // FIXME? Is x-y echanged between Mars CT1 geometry and CT1 definition?
297 fSrcPos->SetXY(-outpars.fypointcorr_deg/fGeom->GetConvMm2Deg(),
298 -outpars.fxpointcorr_deg/fGeom->GetConvMm2Deg());
299 fSrcPos->SetReadyToSave();
300
301 /*
302 --- USEFULL? NEEDED? ---
303 float fcamera_align_angle_deg; // the angle between the camera y-axis and the meridian when a culminating object is observed (defined counter-clockwise looking at the sky)
304 int iratecalc_numevents_odd; // number of events used in the rate calculation (must be odd)
305 int inumpedfile; // number of the pedestal file used
306 int inumpedrun; // number of the pedestal run used in the file (starting at 0)
307 int inumcalfile; // number of the calibration file used
308 int inumlaserrun; // number of the laserrun used in the file (starting at 0)
309 int inumtellogfile; // number of the TelLog file to be used
310 int inumtellogrun; // number of the tellog entry (Runnumber) used from the log file
311 int imcparticle; // CORSIKA-coded Monte Carlo particle type.
312 */
313
314 // ----- preprocessing results -----
315
316 // int istart_mjdate_day; // MJD of run start (first event) */
317 // int iend_mjdate_day; // MJD of run end (last event) */
318 // int irunduration_secs; // difference between start and end time (secs) */
319 *fLog << "Run Time: From " << outpars.istart_mjdate_day << " to ";
320 *fLog << outpars.iend_mjdate_day << " (MJD), Duration=";
321 *fLog << outpars.irunduration_secs/3600 << "h";
322 *fLog << (outpars.irunduration_secs/60)%60 << "m";
323 *fLog << outpars.irunduration_secs%60 << "s" << endl;
324
325 /*
326 --- USEFULL? NEEDED? ---
327 int iproc_mjdate; // MJD of data processing (i.e. creation of this file)
328 */
329
330 // int iproc_evts; // number of events processed */
331 *fLog << "Number of processed events: " << outpars.iproc_evts << endl;
332
333 // --- USEFULL? NEEDED? ---
334 // double dactual_sourcera_hours; // for off runs: the false source (that should have been) observed */
335
336 // float frms_pedsig_phot[iMAXNUMPIX]; // standard deviation of the calibrated signals from the pedestal run */
337 fPedest->InitSize(iMAXNUMPIX);
338 for (Int_t i=0; i<iMAXNUMPIX; i++)
339 (*fPedest)[i].SetMeanRms(outpars.frms_pedsig_phot[i]);
340
341 fPedest->SetReadyToSave();
342
343 // Used to communicate the mean over all pixels
344 // pedestal RMS into the Runs NTuple, as it might
345 // be used for e.g. quality cuts.
346 // float fpedrms_mean;
347 *fLog << "Pedestal RMS: " << outpars.fpedrms_mean << endl;
348
349 // The average current over the active pixels
350 // for this run. This is done separately for
351 // ON and OF runs.
352 //float fcurrent_mean;
353
354 // enum eERRORTOLERANCE {CAUTIOUS=0, GOODPHYSICS, TANK} eerrortolerance;
355 /* 0 == "cautious", exits on any reason (but tells in
356 * the .err file,
357 * 1 == "goodphysics", exits when physics could be affected
358 * by the error,
359 * 2 == "tank", never exits except on coredumps and when
360 * all files have been processed. Do not use such files for
361 * physics analysis!
362 *
363 * NOTE: the capital letter words are the enums, the small letter
364 * words must be used inside the parameter file. */
365
366 // enum eMCTRIGGERFLAG {ALL=0, FLAG, NOFLAG} emctriggerflag;
367 /* all: all events which survive the filter are written to the
368 * events NTuple.
369 * flag: When Dorota's triggerflag is set to 1 for a particular
370 * event, it shall be written to the output. All others shall
371 * just be disregarded. (Default)
372 * noflag: Opposite of 'flag': only events with triggerflag = 0 shall
373 * be treated further. */
374
375 *fLog << "File is a ";
376 *fLog << (outpars.bmontecarlo ? "Monte Carlo" : "Real Data");
377 *fLog << " file." << endl;
378
379 fIsMcFile = outpars.bmontecarlo==TRUE;
380
381}
382
383// --------------------------------------------------------------------------
384//
385// Read CT1 PreProc File Header:
386//
387Int_t MCT1ReadPreProc::ReadRunHeader()
388{
389 char cheadertitle[iHEADERTITLELENGTH];
390 fIn->read(cheadertitle, iHEADERTITLELENGTH);
391
392 TString s = cheadertitle;
393 TString m = cTITLE_TEMPLATE;
394
395 if (!s.BeginsWith(m(0, m.First('%'))))
396 return kFALSE;
397
398 *fLog << cheadertitle << flush;
399
400 // cTITLE_TEMPLATE "PREPROC V%f/S%f CT %d RUN %d %d PROCMJD %d\n"
401 struct outputpars outpars;
402
403 int dummy;
404
405 Float_t fpreprocversion, structversion;
406 sscanf(cheadertitle, cTITLE_TEMPLATE,
407 &fpreprocversion, &structversion,
408 &outpars.itelescope, &outpars.irunnum,
409 &dummy/*&outpars.eruntype*/, &outpars.iproc_mjdate);
410
411 if (fpreprocversion<0.6)
412 {
413 *fLog << err << "Sorry, only files from PreProc V0.6 and newer are supported." << endl;
414 return kFALSE;
415 }
416
417 //
418 // This is a stupid way of getting rid of numerical uncertanties when
419 // comparing floating point numbers (Argh...)
420 //
421 TString s1 = Form("%.2f", structversion);
422 TString s2 = Form("%.2f", STRUCT_VERSION);
423
424 if (s1 != s2)
425 {
426 *fLog << warn << "WARNING: Version of C-structures of file (V";
427 *fLog << s1 << ") not identical with current structures (V";
428 *fLog << s2 << ")" << endl;
429 }
430
431 fIn->read((Byte_t*)&outpars, sizeof(struct outputpars));
432
433 ProcessRunHeader(outpars);
434
435 //rwagner: ReInit whenever new run commences
436 // rc==-1 means: ReInit didn't work out
437 if (!fTaskList->ReInit(fParList)) return -1;
438
439 return kTRUE;
440}
441
442Int_t MCT1ReadPreProc::ReadRunFooter()
443{
444 char cheadertitle[iHEADERTITLELENGTH];
445 fIn->read(cheadertitle, iHEADERTITLELENGTH);
446 /*
447 sscanf(cheadertitle, cEND_EVENTS_TEMPLATE,
448 &filterres.ifilter_passed_evts);
449 */
450
451 TString s = cheadertitle;
452 TString m = cEND_EVENTS_TEMPLATE;
453 Int_t p = m.First('%');
454
455 if (!s.BeginsWith(m(0,p)))
456 {
457 fIn->seekg(-iHEADERTITLELENGTH, ios::cur);
458 return 0;
459 }
460
461 *fLog << inf << cheadertitle << flush;
462
463 struct filterresults filterres;
464 fIn->read((Byte_t*)&filterres, sizeof(struct filterresults));
465 /*
466 int imax_alt_arcs; // maximum altitude reached during the run
467 int iaz_at_max_alt_arcs; // azimuth at the time the max. alt. was reached
468 int itimeaverage_alt_arcs; // altitude averaged over the runtime
469 int icoord_inconsist_evts; // number of events with time-coordinate inconsistency in this run
470 int ifilter_passed_evts; // number of events which passed the filter
471 int imuon_fail_evts; // number of events rejected as muons (other filters passed)
472 int i2out_fail_evts; // number of events which failed in the two out of all pixels software trigger
473 int imuon_and_2out_fail_evts; // number of events failing in both muon and 2out filter
474 int isum_fail_evts; // number of events which failed the sum-of-all-calibrated ADC counts filter
475 int isum_and_muon_fail_evts; // number of events which failed in both the sum and the muon filter
476 int isum_and_2out_fail_evts; // number of events which failed in both the sum and the 2out filter
477 int iall_filters_fail_evts; // number of events which failed in all filters
478 float favg_event_rate_hz; // average rate before filtering
479 float fstddev_event_rate_hz; // standard deviation of the rate before filtering
480 */
481
482 if (fNumEventsInRun!=(UInt_t)filterres.ifilter_passed_evts)
483 {
484 *fLog << err << "ERROR! Number of events in run (" << (UInt_t)filterres.ifilter_passed_evts;
485 *fLog << ") doesn't match number of read events (";
486 *fLog << fNumEventsInRun << ")" << endl;
487 *fLog << " File corrupted." << endl;
488 return -1;
489 }
490
491 fNumFilterEvts += fNumEventsInRun;
492 fNumRuns++;
493
494 *fLog << inf << "Read " << fNumEventsInRun << " events from run (Total=";
495 *fLog << fNumFilterEvts << "/" << fEntries << " [";
496 *fLog << 100*fNumFilterEvts/fEntries << "%], Runs=" << fNumRuns << ")";
497 *fLog << endl;
498
499 return 1;
500}
501
502// --------------------------------------------------------------------------
503//
504// This opens the next file in the list and deletes its name from the list.
505//
506Bool_t MCT1ReadPreProc::OpenNextFile()
507{
508 //
509 // open the input stream and check if it is really open (file exists?)
510 //
511 if (fIn)
512 delete fIn;
513 fIn = NULL;
514
515 //
516 // Check for the existance of a next file to read
517 //
518 TNamed *file = (TNamed*)fFileNames->First();
519 if (!file)
520 return kFALSE;
521
522 //
523 // open the file which is the first one in the chain
524 //
525 const TString name = file->GetName();
526
527 const char *expname = gSystem->ExpandPathName(name);
528 const TString fname(expname);
529 delete [] expname;
530
531 //
532 // Remove this file from the list of pending files
533 //
534 fFileNames->Remove(file);
535
536 *fLog << inf << "Open file: '" << name << "'" << endl;
537
538 if (!CheckHeader(fname))
539 return kFALSE;
540
541 fIn = new ifstream(fname);
542
543 *fLog << inf << "-----------------------------------------------------------------------" << endl;
544
545 switch (ReadRunHeader())
546 {
547 case kFALSE:
548 *fLog << warn << "Unable to read first run header... skipping file." << endl;
549 return kFALSE;
550 break;
551 case -1:
552 *fLog << warn << "ReInit of Tasklist didn't succeed." << endl;
553 return kFALSE;
554 break;
555 default:
556 return kTRUE;
557 break;
558 }
559
560 return kTRUE;
561}
562
563Bool_t MCT1ReadPreProc::CheckHeader(const TString fname) const
564{
565 ifstream fin(fname);
566 if (!fin)
567 {
568 *fLog << dbginf << err << "ERROR - Cannot open file '" << fname << "'" << endl;
569 return kFALSE;
570 }
571
572 char cheadertitle[iHEADERTITLELENGTH];
573 fin.read(cheadertitle, iHEADERTITLELENGTH);
574
575 Float_t fpreprocversion, structversion;
576 Int_t dummyi;
577
578 sscanf(cheadertitle, cTITLE_TEMPLATE,
579 &fpreprocversion, &structversion,
580 &dummyi, &dummyi, &dummyi, &dummyi);
581
582 if (fpreprocversion < 0.6)
583 {
584 *fLog << dbginf << err << "ERROR - You must use PreProc V0.6 or higher." << endl;
585 return kFALSE;
586 }
587
588 if (STRUCT_VERSION > structversion)
589 {
590 *fLog << warn << "WARNING: Version of C-structures of file (V";
591 *fLog << structversion << ") newer than current structures (V";
592 *fLog << STRUCT_VERSION << ")" << endl;
593 }
594
595 *fLog << "Current structures: " << STRUCT_VERSION << " ";
596 *fLog << "Structures in file: " << structversion << " ";
597 *fLog << "Used preproc version: " << fpreprocversion << endl;
598
599 return kTRUE;
600}
601
602
603Int_t MCT1ReadPreProc::GetNumEvents(const TString fname) const
604{
605 *fLog << inf << "Scanning file " << fname << " for size" << flush;
606
607 ifstream fin(fname);
608 if (!fin)
609 {
610 *fLog << dbginf << err << "ERROR - Opening file." << endl;
611 return 0;
612 }
613
614 const TString m(cEND_EVENTS_TEMPLATE);
615 const Int_t p = m.First('%');
616 const TString test = m(0, p);
617
618 Int_t nevts = 0;
619 Int_t nruns = 0;
620
621 while (!fin.eof() && fin.peek()!=EOF)
622 {
623 fin.seekg(iHEADERTITLELENGTH, ios::cur);
624 fin.seekg(sizeof(struct outputpars), ios::cur);
625
626 while (1)
627 {
628 if (fin.peek()==cEND_EVENTS_TEMPLATE[0])
629 {
630 char cheadertitle[iHEADERTITLELENGTH];
631 fin.read(cheadertitle, iHEADERTITLELENGTH);
632
633 const TString s = cheadertitle;
634 if (s.BeginsWith(test))
635 {
636 fin.seekg(sizeof(struct filterresults), ios::cur);
637 nruns++;
638 break;
639 }
640
641 fin.seekg(-iHEADERTITLELENGTH, ios::cur);
642 }
643
644 fin.seekg(sizeof(struct eventrecord), ios::cur);
645 if (fin.eof())
646 break;
647
648 nevts++;
649 }
650 *fLog << "." << flush;
651 }
652
653 *fLog << "done." << endl;
654 *fLog << "Found " << nevts << " events in " << nruns << " runs." << endl;
655
656 return nevts;
657}
658
659// --------------------------------------------------------------------------
660//
661// Open the first file in the list. Check for the output containers or create
662// them if they don't exist.
663//
664// Initialize the size of the MPedestalCam container to 127 pixels (CT1 camera)
665//
666Bool_t MCT1ReadPreProc::PreProcess(MParList *pList)
667{
668
669 fParList = pList;
670
671 //
672 // look for the MCerPhotEvt class in the plist
673 //
674 fNphot = (MCerPhotEvt*)pList->FindCreateObj("MCerPhotEvt");
675 if (!fNphot)
676 return kFALSE;
677
678 //
679 // look for the pedestal class in the plist
680 //
681 fPedest = (MPedestalCam*)pList->FindCreateObj("MPedestalCam");
682 if (!fPedest)
683 return kFALSE;
684
685 //
686 // look for the time class in the plist
687 //
688 fTime = (MTime*)pList->FindCreateObj("MTime");
689 if (!fTime)
690 return kFALSE;
691
692 //
693 // look for the pedestal class in the plist
694 //
695 fBlinds = (MBlindPixels*)pList->FindCreateObj("MBlindPixels");
696 if (!fBlinds)
697 return kFALSE;
698
699 //
700 // look for the source position in the camera
701 //
702 fSrcPos = (MSrcPosCam*)pList->FindCreateObj("MSrcPosCam");
703 if (!fSrcPos)
704 return kFALSE;
705
706 //
707 // look for the camera geometry
708 //
709 fGeom = (MGeomCam*)pList->FindCreateObj("MGeomCamCT1", "MGeomCam");
710 if (!fGeom)
711 return kFALSE;
712
713 //
714 // look for the mc event class
715 //
716 fMcEvt = (MMcEvt*)pList->FindCreateObj("MMcEvt");
717 if (!fMcEvt)
718 return kFALSE;
719
720 //
721 // look for the mc trigger class
722 //
723 fMcTrig = (MMcTrig*)pList->FindCreateObj("MMcTrig");
724 if (!fMcTrig)
725 return kFALSE;
726
727 //
728 // look for the raw run header class
729 //
730 fRawRunHeader = (MRawRunHeader*)pList->FindCreateObj("MRawRunHeader");
731 if (!fRawRunHeader)
732 return kFALSE;
733
734 //
735 // look for the task list
736 //
737 fTaskList = (MTaskList*)pList->FindCreateObj("MTaskList");
738 if (!fTaskList)
739 return kFALSE;
740
741
742 fNumFilterEvts = 0;
743 fNumEvents = 0;
744 fNumRuns = 0;
745
746 return GetSelector() ? GetSelector()->CallPreProcess(pList) : kTRUE;
747}
748
749// --------------------------------------------------------------------------
750//
751// Analyse the event data, means store and copy the needed data into
752// Mars structures and data members
753//
754void MCT1ReadPreProc::ProcessEvent(const struct eventrecord &event)
755{
756 // int isecs_since_midday; // seconds passed since midday before sunset (JD of run start)
757 // int isecfrac_200ns; // fractional part of isecs_since_midday
758 fTime->SetTime(event.isecfrac_200ns, event.isecs_since_midday);
759 fTime->SetReadyToSave();
760
761 /*
762 --- USEFULL? NEEDED? ---
763 short snot_ok_flags; // the bits in these two bytes are flags for additional information on the event: Everything OK =: all Bits = 0
764
765 // for ALT-AZ mount telescopes: rotation angle of the field of
766 // view; this angle is defined mathematically positive looking
767 // towards the sky as the angle between the hour circle through
768 // the object being tracked and the line through pixel 1 and 2
769 int ifieldrot_arcs;
770
771 // event rate in milli Hertz before filtering calculated by
772 // iratecalc_numevents_odd/(time[i+iratecalc_numevents_odd/2] -
773 // time[i-iratecalc_numevents_odd/2])
774 // For the first and the last iratecalc_numevents_odd/2
775 // events the rate is assumed to be constant
776 unsigned short srate_millihz;
777
778 // This is the angle between the observation of this event and the
779 // culmination point. It is going to be written into the events NTuple.
780 float fhourangle;
781 */
782
783 //
784 // read in the number of cerenkov photons and add the 'new' pixel
785 // too the list with it's id, number of photons and error
786 //
787 fNphot->InitSize(iMAXNUMPIX);
788
789 // number of photoelectrons measured in each pixel only the
790 // actual number of pixels (outputpars.inumpixels) is written out
791 // short spixsig_10thphot[iMAXNUMPIX];
792 for (Int_t i=0; i<iMAXNUMPIX; i++)
793 {
794 if (event.spixsig_10thphot[i]==0)
795 continue;
796
797 fNphot->AddPixel(i, 0.1*event.spixsig_10thphot[i],
798 (*fPedest)[i].GetMeanRms());
799 }
800 fNphot->SetReadyToSave();
801
802 // int ipreproc_alt_arcs; // "should be" alt according to preproc (arcseconds)
803 // int ipreproc_az_arcs; // "should be" az according to preproc (arcseconds)
804
805 fMcEvt->Fill(0, /*fEvtNum*/
806 fIsMcFile ? event.imcparticle : 0, /*corsika particle type*/
807 fIsMcFile ? event.fmcenergy_tev*1000 : 0,
808 0, /* fThi0 */
809 0, /* fFirTar */
810 0, /* fzFirInt */
811 // 0, /* fThet*/
812 // rwagner: The following should be theta, right? Check with
813 // altitude fill some lines down...
814 0, // altitude (arcseconds)
815 0, /* fPhii */
816 0, /* fCorD */
817 0, /* fCorX */
818 0, /* fCorY */
819 fIsMcFile ? event.imcimpact_m*100 : 0,
820 TMath::Pi()/180*event.iaz_arcs/3600, // azimuth (arcseconds)
821 TMath::Pi()*(0.5-1./180*event.ialt_arcs/3600), // altitude (arcseconds)
822 0, /* fTFirst */
823 0, /* fTLast */
824 0, /* fL_Nmax */
825 0, /* fL_t0 */
826 0, /* fL_tmax */
827 0, /* fL_a */
828 0, /* fL_b */
829 0, /* fL_c */
830 0, /* fL_chi2 */
831 0, /* uiPin */
832 0, /* uiPat */
833 0, /* uiPre */
834 0, /* uiPco */
835 0, /* uiPelS */
836 fIsMcFile ? event.fmcsize_phel : 0, /* uiPelC, Simulated SIZE */
837 0, /* elec */
838 0, /* muon */
839 0 /* other */
840 );
841
842 fMcTrig->SetFirstLevel(event.imctriggerflag); // MC data from Dorota get a triggerflag: 1 means triggered, 0 not. */
843
844 fMcTrig->SetReadyToSave();
845 fMcEvt->SetReadyToSave();
846}
847
848// --------------------------------------------------------------------------
849//
850// Because of the file format of the preproc output we have to check at any
851// event where in the file stream we are...
852//
853Bool_t MCT1ReadPreProc::CheckFilePosition()
854{
855 //
856 // Means: If no file is open (first call) try to open the first file
857 //
858 if (!fIn)
859 return kFALSE;
860
861 //
862 // Because we can have 0-event runs in the file we loop as often
863 // as we don't find a new footer-header combination.
864 //
865 while (1)
866 {
867 //
868 // If the first character isn't the first of the footer it must be
869 // an event
870 //
871 if (fIn->peek()!=cEND_EVENTS_TEMPLATE[0])
872 return kTRUE;
873
874 //
875 // Try reading the footer. If this isn't successful...
876 // must be an event
877 //
878 switch (ReadRunFooter())
879 {
880 case -1:
881 return kFALSE;
882 case 0:
883 return kTRUE;
884 }
885
886 *fLog << inf << "Footer found." << endl;
887
888 const char c = fIn->peek();
889
890 //
891 // No after reading the footer check if we reached the end of the file
892 //
893 if (fIn->eof() || c==EOF)
894 {
895 *fLog << "End of file." << endl;
896 return kFALSE;
897 }
898
899 //
900 // If the eof isn't reached a new header must follow. Check for it.
901 //
902 if (c!=cTITLE_TEMPLATE[0])
903 {
904 *fLog << inf << "Error finding new run header in file (possible EOF)... skipping rest of file." << endl;
905 return kFALSE;
906 }
907
908 *fLog << "-----------------------------------------------------------------------" << endl;
909
910 switch (ReadRunHeader())
911 {
912 case kTRUE:
913 return kTRUE;
914 break;
915 case -1:
916 *fLog << warn << "ReInit of Tasklist didn't succeed." << endl;
917 return kFALSE;
918 break;
919 default:
920 return kTRUE;
921 }
922 }
923}
924
925// --------------------------------------------------------------------------
926//
927// Check for the event number and depending on this number decide if
928// pedestals or event data has to be read.
929//
930// If the end of the file is reached try to open the next in the list. If
931// there is now next file stop the eventloop.
932//
933Bool_t MCT1ReadPreProc::Process()
934{
935 //
936 // Check where in the file we are. If neither a new event, nor a
937 // footer/header combination is detected go to the next file.
938 //
939 if (!CheckFilePosition())
940 if (!OpenNextFile())
941 return kFALSE;
942
943 //
944 // Check for a selector. If one is given and returns kFALSE
945 // skip this event.
946 //
947 if (GetSelector())
948 {
949 //
950 // Make sure selector is processed
951 //
952 if (!GetSelector()->CallProcess())
953 {
954 *fLog << err << dbginf << "Processing Selector failed." << endl;
955 return kFALSE;
956 }
957
958 //
959 // Skip Event
960 //
961 if (!GetSelector()->IsExpressionTrue())
962 {
963 fIn->seekg(sizeof(struct eventrecord), ios::cur);
964
965 fNumEvents++;
966 fNumEventsInRun++;
967
968 return kCONTINUE;
969 }
970 }
971
972 // event data to be read from the file
973 struct eventrecord event;
974
975 // read the eventrecord from the file
976 fIn->read((Byte_t*)&event, sizeof(struct eventrecord));
977
978 ProcessEvent(event);
979
980 fNumEvents++;
981 fNumEventsInRun++;
982
983 return kTRUE;
984}
985
986Bool_t MCT1ReadPreProc::PostProcess()
987{
988 *fLog << all;
989 *fLog << "Number events passed the filter: " << fNumFilterEvts << endl;
990 *fLog << "Number of Events read from file: " << fNumEvents << endl;
991 *fLog << "Number of Runs read from file: " << fNumRuns << endl;
992 *fLog << "Number of events detected first: " << fEntries << endl;
993
994 if (fNumEvents!=fNumFilterEvts)
995 {
996 *fLog << warn << "WARNING! Number of events in file doesn't match number of read events..." << endl;
997 *fLog << " File might be corrupt." << endl;
998 }
999
1000 return GetSelector() ? GetSelector()->CallPostProcess() : kTRUE;
1001}
Note: See TracBrowser for help on using the repository browser.