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

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