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

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