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

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