source: branches/MarsMoreSimulationTruth/manalysisct1/MCT1ReadPreProc.cc

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