/* ======================================================================== *\ ! ! * ! * This file is part of MARS, the MAGIC Analysis and Reconstruction ! * Software. It is distributed to you in the hope that it can be a useful ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. ! * It is distributed WITHOUT ANY WARRANTY. ! * ! * Permission to use, copy, modify and distribute this software and its ! * documentation for any purpose is hereby granted without fee, ! * provided that the above copyright notice appear in all copies and ! * that both that copyright notice and this permission notice appear ! * in supporting documentation. It is provided "as is" without express ! * or implied warranty. ! * ! ! ! Author(s): Thomas Bretz 12/2000 ! Author(s): Harald Kornmayer 1/2001 (harald@mppmu.mpg.de) ! ! Copyright: MAGIC Software Development, 2000-2001 ! ! \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// // // MCT1ReadPreProc // // Reads a output file of the CT1 preproc. // // Input Containers: // -/- // // Output Containers: // MCerPhotEvt the data container for all data. // MPedestalCam ct1 pedestals // MMcEvt monte carlo data container for MC files // MMcTrig mc data container for trigger information // MSrcPosCam source position in the camera // MBlindPixels Array holding blind pixels // ///////////////////////////////////////////////////////////////////////////// #include "MCT1ReadPreProc.h" #include #include #include #define LINUX #define HISTO void #define HBOOK_FILE int #include "defines.h" #include "structures.h" #include "MLog.h" #include "MLogManip.h" #include "MParList.h" #include "MCerPhotEvt.h" #include "MPedestalPix.h" #include "MPedestalCam.h" #include "MGeomCam.h" #include "MSrcPosCam.h" #include "MBlindPixels.h" #include "MMcEvt.hxx" #include "MMcTrig.hxx" ClassImp(MCT1ReadPreProc); // -------------------------------------------------------------------------- // // Default constructor. Creates an array which stores the file names of // the files which should be read. If a filename is given it is added // to the list. // MCT1ReadPreProc::MCT1ReadPreProc(const char *fname, const char *name, const char *title) : fIn(NULL), fEntries(0) { fName = name ? name : "MRead"; fTitle = title ? title : "Task to loop over events in CT1 ascii file"; // // remember file name for opening the file in the preprocessor // fFileNames = new TList; fFileNames->SetOwner(); if (fname) AddFile(fname); } // -------------------------------------------------------------------------- // // Delete the filename list and the input stream if one exists. // MCT1ReadPreProc::~MCT1ReadPreProc() { delete fFileNames; if (fIn) delete fIn; } // -------------------------------------------------------------------------- // // Add this file as the last entry in the chain // void MCT1ReadPreProc::AddFile(const char *txt) { ifstream *store = fIn; fIn = new ifstream(gSystem->ExpandPathName(txt)); if (!(*fIn)) { *fLog << warn << "Cannot open file '" << txt << "'... ignored." << endl; fIn = store; return; } fEntries += GetNumEvents(); delete fIn; fIn = store; fFileNames->AddLast(new TNamed(txt, "")); } // -------------------------------------------------------------------------- // // Print data from the header to the screen and analyse the header data, // means store and copy the needed data into Mars structures and // data members // void MCT1ReadPreProc::ProcessRunHeader(const struct outputpars &outpars) { if (outpars.inumpixels != iMAXNUMPIX) *fLog << warn << "WARNING! File doesn't contain " << iMAXNUMPIX << " Pixels... maybe corrupt." << endl; fNumEventsInRun = 0; // // ------------------- Output some stuff ----------------------- // // int itelescope; // number of the CT which took the data *fLog << inf << "Telescope: CT" << outpars.itelescope; // float flongitude_deg; // longitude (counted positive towards West) of CT position */ // float flatitude_deg; // latitude (counted positive towards North) of CT position */ *fLog << " located @ Longitude=" << outpars.flongitude_deg; *fLog << "deg Latitude=" << outpars.flatitude_deg << "deg" << endl; // int irunnum; // run number (from parameters file) // enum onoroff {NEITHER_ON_NOR_OFF, OFF_SOURCE, ON_SOURCE} eruntype; // runtype *fLog << "Run: #" << outpars.irunnum << " ("; switch (outpars.eruntype) { case NEITHER_ON_NOR_OFF: *fLog << "unknown"; break; case OFF_SOURCE: *fLog << "off-source"; break; case ON_SOURCE: *fLog << "on-source"; break; default: *fLog << (int)outpars.eruntype; break; } *fLog << ", "; switch (outpars.etrackmode) { case NORMAL: *fLog << "normal tracking"; break; case REVERSE: *fLog << "reverse tracking"; break; case DUNNO: *fLog << "unknown tracking"; break; default: *fLog << (int)outpars.etrackmode; break; } *fLog << ")" << endl; //double dsourcera_hours; // right ascension of observed source in hours //double dsourcedec_deg; // declination of observed source in degrees *fLog << "Source: RA=" << outpars.dsourcera_hours << "h DEC="; *fLog << outpars.dsourcedec_deg << "deg" << endl; //int inummuonpixels; // number of pixels in the muon shield //int inumcointdcs; // number of coincidence tdcs recorded in the runfile //float fpixdiameter_deg; // smallest pixel diameter (degrees) (from parameters file) */ // enum axes {RA, DEC, ALT, AZ} ese1_is; // name of the axis to which shaft encoder 1 is attached (implies the type of mount) *fLog << "Shaftencoder 1 @ "; switch (outpars.ese1_is) { case RA: *fLog << "RA"; break; case DEC: *fLog << "DEC"; break; case ALT: *fLog << "ALT"; break; case AZ: *fLog << "AZ"; break; default: *fLog << (int)outpars.ese1_is; break; } *fLog << endl; // int isezeropos[2]; // zero position of shaftencoders 1 and 2 (from parameters file) *fLog << "SE Zero: SE(1)=" << outpars.isezeropos[0] << " "; *fLog << "SE(2)=" << outpars.isezeropos[1] << endl; // int iaz_rev_track_corr; // correction for the azimuth shaft encoder (ALT/AZ mount only) in reverse tracking mode // int ialt_rev_track_corr; // correction for the altitude shaft encoder (ALT/AZ mount only) in reverse tracking mode *fLog << "Reverse tracking corrections: SE(az)=" << outpars.iaz_rev_track_corr; *fLog << " SE(alt)=" << outpars.ialt_rev_track_corr << endl; // float fbendingcorr; // bending correction factor (ALT/AZ mount only) // float fextinction; // atmospheric extinction (typically taken from the Carlsberg Meridian Circle data) *fLog << "Bending: Correction factor=" << outpars.fbendingcorr << " "; *fLog << "Extinction=" << outpars.fextinction << endl; // Boolean bdontusepix[iMAXNUMPIX]; // bdontusepix is set true if the pixel should not be used in image analysis, otherwise it is true; fBlinds->Clear(); *fLog << "Don't use pixels: "; for (int i=0; iSetPixelBlind(i); } *fLog << endl; *fLog << "Exclude: "; // Boolean bexcludepix[iMAXNUMPIX]; for (int i=0; iSetPixelBlind(i); } *fLog << endl; /* bexcludepix[] is set TRUE (== exclude from pedestal, Laser * calibration and the further analysis) when the Mean value * of a pixel inside a pedestal Run is larger than 50 or ( || ) * if the pedestal RMS value of this pixel is larger than 5.0 * This is reflected in the (new for versions >= 0.4) * variable "pixonoff" in the ntuple written by preproc: * preproc.nt.hbook * * When the pixel is excluded by the user it will get a -2 otherwise * pixonoff = 0. * Additive to this a -1 is added when preproc excludes the pixel * for a given Run. So the actual value tells you whether you caught * it already by hand or not. * * A plot of pixonoff may also be handy to tell you the status of your * ADC equipment. */ // float fphotoel_per_adccnt[iMAXNUMPIX]; // conversion factors for the pixel signals */ /* float padc = outpars.fphotoel_per_adccnt[0]; *fLog << "Phe/ADC (pixel 0): " << padc << endl; for (int i=0; iSetXY(-outpars.fypointcorr_deg/fGeom->GetConvMm2Deg(), -outpars.fxpointcorr_deg/fGeom->GetConvMm2Deg()); fSrcPos->SetReadyToSave(); /* --- USEFULL? NEEDED? --- 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) int iratecalc_numevents_odd; // number of events used in the rate calculation (must be odd) int inumpedfile; // number of the pedestal file used int inumpedrun; // number of the pedestal run used in the file (starting at 0) int inumcalfile; // number of the calibration file used int inumlaserrun; // number of the laserrun used in the file (starting at 0) int inumtellogfile; // number of the TelLog file to be used int inumtellogrun; // number of the tellog entry (Runnumber) used from the log file int imcparticle; // CORSIKA-coded Monte Carlo particle type. */ // ----- preprocessing results ----- // int istart_mjdate_day; // MJD of run start (first event) */ // int iend_mjdate_day; // MJD of run end (last event) */ // int irunduration_secs; // difference between start and end time (secs) */ *fLog << "Run Time: From " << outpars.istart_mjdate_day << " to "; *fLog << outpars.iend_mjdate_day << " (MJD), Duration="; *fLog << outpars.irunduration_secs/3600 << "h"; *fLog << (outpars.irunduration_secs/60)%60 << "m"; *fLog << outpars.irunduration_secs%60 << "s" << endl; /* --- USEFULL? NEEDED? --- int iproc_mjdate; // MJD of data processing (i.e. creation of this file) */ // int iproc_evts; // number of events processed */ *fLog << "Number of processed events: " << outpars.iproc_evts << endl; // --- USEFULL? NEEDED? --- // double dactual_sourcera_hours; // for off runs: the false source (that should have been) observed */ // float frms_pedsig_phot[iMAXNUMPIX]; // standard deviation of the calibrated signals from the pedestal run */ fPedest->InitSize(iMAXNUMPIX); for (Int_t i=0; iSetReadyToSave(); } // -------------------------------------------------------------------------- // // Read CT1 PreProc File Header: // Bool_t MCT1ReadPreProc::ReadRunHeader() { char cheadertitle[iHEADERTITLELENGTH]; fIn->read(cheadertitle, iHEADERTITLELENGTH); TString s = cheadertitle; TString m = cTITLE_TEMPLATE; if (!s.BeginsWith(m(0, m.First('%')))) return kFALSE; *fLog << cheadertitle << flush; // cTITLE_TEMPLATE "PREPROC V%f/S%f CT %d RUN %d %d PROCMJD %d\n" struct outputpars outpars; int dummy; Float_t fpreprocversion, structversion; sscanf(cheadertitle, cTITLE_TEMPLATE, &fpreprocversion, &structversion, &outpars.itelescope, &outpars.irunnum, &dummy/*&outpars.eruntype*/, &outpars.iproc_mjdate); if (STRUCT_VERSION != structversion) { *fLog << warn << "WARNING: Version of C-structures of file (V"; *fLog << structversion << ") not identical with current structures (V"; *fLog << STRUCT_VERSION << ")" << endl; } fIn->read((Byte_t*)&outpars, sizeof(struct outputpars)); ProcessRunHeader(outpars); return kTRUE; } Bool_t MCT1ReadPreProc::ReadRunFooter() { char cheadertitle[iHEADERTITLELENGTH]; fIn->read(cheadertitle, iHEADERTITLELENGTH); /* sscanf(cheadertitle, cEND_EVENTS_TEMPLATE, &filterres.ifilter_passed_evts); */ TString s = cheadertitle; TString m = cEND_EVENTS_TEMPLATE; Int_t p = m.First('%'); if (!s.BeginsWith(m(0,p))) { fIn->seekg(-iHEADERTITLELENGTH, ios::cur); return kFALSE; } *fLog << inf << cheadertitle << flush; struct filterresults filterres; fIn->read((Byte_t*)&filterres, sizeof(struct filterresults)); /* int imax_alt_arcs; // maximum altitude reached during the run int iaz_at_max_alt_arcs; // azimuth at the time the max. alt. was reached int itimeaverage_alt_arcs; // altitude averaged over the runtime int icoord_inconsist_evts; // number of events with time-coordinate inconsistency in this run int ifilter_passed_evts; // number of events which passed the filter int imuon_fail_evts; // number of events rejected as muons (other filters passed) int i2out_fail_evts; // number of events which failed in the two out of all pixels software trigger int imuon_and_2out_fail_evts; // number of events failing in both muon and 2out filter int isum_fail_evts; // number of events which failed the sum-of-all-calibrated ADC counts filter int isum_and_muon_fail_evts; // number of events which failed in both the sum and the muon filter int isum_and_2out_fail_evts; // number of events which failed in both the sum and the 2out filter int iall_filters_fail_evts; // number of events which failed in all filters float favg_event_rate_hz; // average rate before filtering float fstddev_event_rate_hz; // standard deviation of the rate before filtering */ fNumFilterEvts += filterres.ifilter_passed_evts; fNumRuns++; *fLog << inf << "Read " << fNumEventsInRun << " events from run." << endl; if (fNumEventsInRun!=(UInt_t)filterres.ifilter_passed_evts) { *fLog << warn << "WARNING! Number of events in run doesn't match number of read events." << endl; *fLog << " File might be corrupt." << endl; } return kTRUE; } // -------------------------------------------------------------------------- // // This opens the next file in the list and deletes its name from the list. // Bool_t MCT1ReadPreProc::OpenNextFile() { // // open the input stream and check if it is really open (file exists?) // if (fIn) delete fIn; fIn = NULL; // // Check for the existance of a next file to read // TNamed *file = (TNamed*)fFileNames->First(); if (!file) return kFALSE; // // open the file which is the first one in the chain // const char *name = file->GetName(); fIn = new ifstream(gSystem->ExpandPathName(name)); if (!(*fIn)) { *fLog << dbginf << "Cannot open file '" << name << "'" << endl; return kFALSE; } *fLog << "Open file: '" << name << "'" << endl; // // Remove this file from the list of pending files // fFileNames->Remove(file); *fLog << inf << "-----------------------------------------------------------------------" << endl; *fLog << "File contains " << GetNumEvents() << " events." << endl; // WORKAROUND for not working seekg(0) in GetNumEvents fIn->close(); fIn->open(gSystem->ExpandPathName(name)); Bool_t rc = ReadRunHeader(); if (!rc) *fLog << warn << "Unable to read first run header... skipping file." << endl; return rc; } Int_t MCT1ReadPreProc::GetNumEvents() { *fLog << inf << "Scanning file for size" << flush; const TString m(cEND_EVENTS_TEMPLATE); const Int_t p = m.First('%'); const TString test = m(0, p); Int_t nevts = 0; Int_t nruns = 0; while (!fIn->eof()) { fIn->seekg(iHEADERTITLELENGTH, ios::cur); fIn->seekg(sizeof(struct outputpars), ios::cur); while (1) { if (fIn->peek()==cEND_EVENTS_TEMPLATE[0]) { char cheadertitle[iHEADERTITLELENGTH]; fIn->read(cheadertitle, iHEADERTITLELENGTH); const TString s = cheadertitle; if (s.BeginsWith(test)) { fIn->seekg(sizeof(struct filterresults), ios::cur); nruns++; break; } fIn->seekg(-iHEADERTITLELENGTH, ios::cur); } fIn->seekg(sizeof(struct eventrecord), ios::cur); if (fIn->eof()) break; nevts++; } *fLog << "." << flush; } *fLog << "done." << endl; *fLog << "Found " << nevts << " events in " << nruns << " runs." << endl; return nevts; } // -------------------------------------------------------------------------- // // Open the first file in the list. Check for the output containers or create // them if they don't exist. // // Initialize the size of the MPedestalCam container to 127 pixels (CT1 camera) // Bool_t MCT1ReadPreProc::PreProcess(MParList *pList) { // // look for the MCerPhotEvt class in the plist // fNphot = (MCerPhotEvt*)pList->FindCreateObj("MCerPhotEvt"); if (!fNphot) return kFALSE; // // look for the pedestal class in the plist // fPedest = (MPedestalCam*)pList->FindCreateObj("MPedestalCam"); if (!fPedest) return kFALSE; // // look for the pedestal class in the plist // fBlinds = (MBlindPixels*)pList->FindCreateObj("MBlindPixels"); if (!fBlinds) return kFALSE; // // look for the source position in the camera // fSrcPos = (MSrcPosCam*)pList->FindCreateObj("MSrcPosCam", "Source"); if (!fSrcPos) return kFALSE; // // look for the camera geometry // fGeom = (MGeomCam*)pList->FindCreateObj("MGeomCamCT1", "MGeomCam"); if (!fGeom) return kFALSE; // // look for the mc event class // fMcEvt = (MMcEvt*)pList->FindCreateObj("MMcEvt"); if (!fMcEvt) return kFALSE; // // look for the mc trigger class // fMcTrig = (MMcTrig*)pList->FindCreateObj("MMcTrig"); if (!fMcTrig) return kFALSE; fNumFilterEvts = 0; fNumRuns = 0; return kTRUE; } // -------------------------------------------------------------------------- // // Analyse the event data, means store and copy the needed data into // Mars structures and data members // void MCT1ReadPreProc::ProcessEvent(const struct eventrecord &event) { /* --- USEFULL? NEEDED? --- int isecs_since_midday; // seconds passed since midday before sunset (JD of run start) int isecfrac_200ns; // fractional part of isecs_since_midday short snot_ok_flags; // the bits in these two bytes are flags for additional information on the event: Everything OK =: all Bits = 0 // for ALT-AZ mount telescopes: rotation angle of the field of // view; this angle is defined mathematically positive looking // towards the sky as the angle between the hour circle through // the object being tracked and the line through pixel 1 and 2 int ifieldrot_arcs; // event rate in milli Hertz before filtering calculated by // iratecalc_numevents_odd/(time[i+iratecalc_numevents_odd/2] - // time[i-iratecalc_numevents_odd/2]) // For the first and the last iratecalc_numevents_odd/2 // events the rate is assumed to be constant unsigned short srate_millihz; // This is the angle between the observation of this event and the // culmination point. It is going to be written into the events NTuple. float fhourangle; */ // // read in the number of cerenkov photons and add the 'new' pixel // too the list with it's id, number of photons and error // fNphot->InitSize(iMAXNUMPIX); // number of photoelectrons measured in each pixel only the // actual number of pixels (outputpars.inumpixels) is written out // short spixsig_10thphot[iMAXNUMPIX]; for (Int_t i=0; iAddPixel(i, 0.1*event.spixsig_10thphot[i], (*fPedest)[i].GetMeanRms()); } fNphot->SetReadyToSave(); // int ipreproc_alt_arcs; // "should be" alt according to preproc (arcseconds) // int ipreproc_az_arcs; // "should be" az according to preproc (arcseconds) fMcEvt->SetTheta(TMath::Pi()*(0.5-1./180*event.ialt_arcs/3600)); // altitude (arcseconds) fMcEvt->SetPhi(TMath::Pi()/180*event.iaz_arcs/3600); // azimuth (arcseconds) fMcEvt->SetReadyToSave(); if (!fIsMcFile) return; fMcEvt->SetPartId(event.imcparticle); // corsika particle type fMcEvt->SetEnergy(event.fmcenergy_tev*1000); // simulated energy fMcEvt->SetImpact(event.imcimpact_m*100); // simulated impact fMcTrig->SetFirstLevel(event.imctriggerflag); // MC data from Dorota get a triggerflag: 1 means triggered, 0 not. */ fMcTrig->SetReadyToSave(); //float fmcsize_phel; // Simulated SIZE } // -------------------------------------------------------------------------- // // Because of the file format of the preproc output we have to check at any // event where in the file stream we are... // Bool_t MCT1ReadPreProc::CheckFilePosition() { // // Means: If no file is open (first call) try to open the first file // if (!fIn) return kFALSE; // // Because we can have 0-event runs in the file we loop as often // as we don't find a new footer-header combination. // while (1) { // // If the first character isn't the first of the footer it must be // an event // if (fIn->peek()!=cEND_EVENTS_TEMPLATE[0]) return kTRUE; // // Try reading the footer. If it isn't succefull... // must be an event // if (!ReadRunFooter()) return kTRUE; *fLog << inf << "Footer found." << endl; // // No after reading the footer check if we reached the end of the file // if (fIn->eof()) { *fLog << "End of file." << endl; return kFALSE; } // // If the eof isn't reached a new header must follow. Check for it. // if (fIn->peek()!=cTITLE_TEMPLATE[0]) { *fLog << inf << "Error finding new run header in file (possible EOF)... skipping rest of file." << endl; return kFALSE; } *fLog << "-----------------------------------------------------------------------" << endl; if (!ReadRunHeader()) return kTRUE; } } // -------------------------------------------------------------------------- // // Check for the event number and depending on this number decide if // pedestals or event data has to be read. // // If the end of the file is reached try to open the next in the list. If // there is now next file stop the eventloop. // Bool_t MCT1ReadPreProc::Process() { // // Check where in the file we are. If neither a new event, nor a // footer/header combination is detected go to the next file. // if (!CheckFilePosition()) if (!OpenNextFile()) return kFALSE; // event data to be read from the file struct eventrecord event; // read the eventrecord from the file fIn->read((Byte_t*)&event, sizeof(struct eventrecord)); ProcessEvent(event); fNumEvents++; fNumEventsInRun++; return kTRUE; } Bool_t MCT1ReadPreProc::PostProcess() { *fLog << all; *fLog << "Number events passed the filter: " << fNumFilterEvts << endl; *fLog << "Number of Events read from file: " << fNumEvents << endl; *fLog << "Number of Runs read from file: " << fNumRuns << endl; if (fNumEvents!=fNumFilterEvts) { *fLog << warn << "WARNING! Number of events in file doesn't match number of read events..." << endl; *fLog << " File might be corrupt." << endl; } return kTRUE; }