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

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