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

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