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

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