Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 4459) +++ trunk/MagicSoft/Mars/Changelog (revision 4460) @@ -87,5 +87,5 @@ - new - * manalysis/MCT1*: + * manalysis/MCT1*, mfileio/structures.h, mfileio/MCT1Read*.[h,cc]: - moved to new directory manalysisct1 @@ -94,4 +94,7 @@ - removed MPad because long time ago changed discuseed were not made (it doesn't compile anymore, because MBlindPixels is missing now) + + * mfileio/Makefile, mfileio/FileIOLinkDef.h: + - updated Index: trunk/MagicSoft/Mars/manalysisct1/AnalysisCT1LinkDef.h =================================================================== --- trunk/MagicSoft/Mars/manalysisct1/AnalysisCT1LinkDef.h (revision 4459) +++ trunk/MagicSoft/Mars/manalysisct1/AnalysisCT1LinkDef.h (revision 4460) @@ -11,4 +11,6 @@ #pragma link C++ class MCT1SupercutsCalc+; #pragma link C++ class MCT1FindSupercuts+; +#pragma link C++ class MCT1ReadPreProc+; +#pragma link C++ class MCT1ReadPreAscii+; #endif Index: trunk/MagicSoft/Mars/manalysisct1/MCT1ReadAscii.cc =================================================================== --- trunk/MagicSoft/Mars/manalysisct1/MCT1ReadAscii.cc (revision 4460) +++ trunk/MagicSoft/Mars/manalysisct1/MCT1ReadAscii.cc (revision 4460) @@ -0,0 +1,304 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MCT1ReadAscii +// +// Reads a ascii file with CT1 data. The file description and some example +// files can be found on the Magic homepage. +// +// Input Containers: +// -/- +// +// Output Containers: +// MCerPhotEvt +// +///////////////////////////////////////////////////////////////////////////// +#include "MCT1ReadAscii.h" + +#include +#include + +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" +#include "MCerPhotEvt.h" + +#include "MPedPhotPix.h" +#include "MPedPhotCam.h" + +ClassImp(MCT1ReadAscii); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// 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. +// +MCT1ReadAscii::MCT1ReadAscii(const char *fname, + const char *name, + const char *title) + : fIn(NULL) +{ + fName = name ? name : "MCT1ReadAscii"; + 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. +// +MCT1ReadAscii::~MCT1ReadAscii() +{ + delete fFileNames; + if (fIn) + delete fIn; +} + +// -------------------------------------------------------------------------- +// +// Add this file as the last entry in the chain +// +Int_t MCT1ReadAscii::AddFile(const char *txt, Int_t) +{ + TNamed *name = new TNamed(txt, ""); + fFileNames->AddLast(name); + return 1; +} + +// -------------------------------------------------------------------------- +// +// This opens the next file in the list and deletes its name from the list. +// +Bool_t MCT1ReadAscii::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(); + + const char *expname = gSystem->ExpandPathName(name); + fIn = new ifstream(expname); + + const Bool_t noexist = !(*fIn); + if (noexist) + { + *fLog << err << "Cannot open file " << expname << ": "; + *fLog << strerror(errno) << endl; + } + else + *fLog << inf << "Open file: '" << name << "'" << endl; + + delete [] expname; + + // + // Remove this file from the list of pending files + // + fFileNames->Remove(file); + + return !noexist; +} + +// -------------------------------------------------------------------------- +// +// 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 MPedPhotCam container to 127 pixels (CT1 camera) +// +Int_t MCT1ReadAscii::PreProcess(MParList *pList) +{ + // + // Preprocessing + // + + // + // Try to open at least one (the first) file + // + if (!OpenNextFile()) + return kFALSE; + + // + // 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 = (MPedPhotCam*)pList->FindCreateObj("MPedPhotCam"); + if (!fPedest) + return kFALSE; + + MGeomCam *geom = (MGeomCam*)pList->FindObject("MGeomCam"); + if (!geom) + { + *fLog << err << "MGeomCam not found... abort." << endl; + return kFALSE; + } + + fPedest->Init(*geom); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Read apedestal entry (line) from the file +// +void MCT1ReadAscii::ReadPedestals() +{ + // + // skip the next 4 values + // + Float_t val; + + *fIn >> val; + *fIn >> val; + *fIn >> val; + *fIn >> val; + + // + // read in the next 127 numbers as the pedestals + // + for (Int_t i = 0; i<127; i++) + { + *fIn >> val; + + if (val>0) + (*fPedest)[i].Set(0, val); + } +} + +// -------------------------------------------------------------------------- +// +// Read a data entry (line) from the file +// +void MCT1ReadAscii::ReadData() +{ + // + // five unsused numbers + // + Int_t val; + + *fIn >> val; // ener + *fIn >> val; // zenang + *fIn >> val; // sec1 + *fIn >> val; // sec2 + + // + // read in the number of cerenkov photons and add the 'new' pixel + // too the list with it's id, number of photons and error + // + for (Int_t i = 0; i<127; i++ ) + { + Float_t nphot; + + *fIn >> nphot; + + if (nphot > 0.0) + fNphot->AddPixel(i, nphot, (*fPedest)[i].GetRms()); + } + fNphot->FixSize(); +} + +// -------------------------------------------------------------------------- +// +// 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. +// +Int_t MCT1ReadAscii::Process() +{ + // + // FIXME. This function should switch between reading pedestals and + // reading event data by the 'switch entry'. + // After reading it should set the InputStreamID correctly. + // (should use MPedPhotCam ) + // + + // + // read in the event nr + // + Int_t evtnr; + *fIn >> evtnr; + + // + // check if we are done. Try to open the next file in chain. + // If it was possible start reading. If not break the event loop + // + if (fIn->eof()) + return OpenNextFile() ? kCONTINUE : kFALSE; + + // + // if the first number is negativ this is a pedestal line: + // read in pedestals + // + // FIXME! Set InputStreamID + + if (evtnr < 0) + { + ReadPedestals(); + return kCONTINUE; + } + + ReadData(); + + return kTRUE; +} Index: trunk/MagicSoft/Mars/manalysisct1/MCT1ReadAscii.h =================================================================== --- trunk/MagicSoft/Mars/manalysisct1/MCT1ReadAscii.h (revision 4460) +++ trunk/MagicSoft/Mars/manalysisct1/MCT1ReadAscii.h (revision 4460) @@ -0,0 +1,41 @@ +#ifndef MARS_MCT1ReadAscii +#define MARS_MCT1ReadAscii + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +class TList; +class MCerPhotEvt; +class MPedPhotCam; + +class MCT1ReadAscii : public MTask +{ +private: + ifstream *fIn; // the inputfile + MCerPhotEvt *fNphot; // the data container for all data. + MPedPhotCam *fPedest; // CT1 pedestals + TList *fFileNames; // Array which stores the \0-terminated filenames + + Bool_t OpenNextFile(); + + void ReadPedestals(); + void ReadData(); + + Int_t PreProcess(MParList *pList); + Int_t Process(); + +public: + MCT1ReadAscii(const char *filename=NULL, + const char *name=NULL, + const char *title=NULL); + + ~MCT1ReadAscii(); + + Int_t AddFile(const char *fname, Int_t dummy=-1); + + ClassDef(MCT1ReadAscii, 0) // Reads the CT1 data file +}; + +#endif + Index: trunk/MagicSoft/Mars/manalysisct1/MCT1ReadPreProc.cc =================================================================== --- trunk/MagicSoft/Mars/manalysisct1/MCT1ReadPreProc.cc (revision 4460) +++ trunk/MagicSoft/Mars/manalysisct1/MCT1ReadPreProc.cc (revision 4460) @@ -0,0 +1,1171 @@ +/* ======================================================================== *\ +! +! * +! * 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, 11/2002 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MCT1ReadPreProc +// +// Reads a output file of the CT1 preproc. +// +// Implements usage of a selector (see MRead) Use such a filter to skip +// events before reading! But never use a filter which needs read data +// as input... +// +// Input Containers: +// -/- +// +// Output Containers: +// MCerPhotEvt the data container for all data. +// MPedPhotCam CT1 pedestals +// MMcEvt monte carlo data container for MC files +// MMcTrig mc data container for trigger information +// MSrcPosCam source position in the camera +// MBadPixelsCam Array holding blind pixels +// +///////////////////////////////////////////////////////////////////////////// +#include "MCT1ReadPreProc.h" + +#include +#include + +#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 "MTime.h" +#include "MFilter.h" + +#include "MParList.h" +#include "MCerPhotEvt.h" + +#include "MPedPhotPix.h" +#include "MPedPhotCam.h" + +#include "MGeomCam.h" +#include "MSrcPosCam.h" +#include "MBadPixelsCam.h" +#include "MBadPixelsPix.h" + +#include "MRawRunHeader.h" +#include "MTaskList.h" + +#include "MMcEvt.hxx" +#include "MMcTrig.hxx" +#include "MBinning.h" + +#include "MParameters.h" + +ClassImp(MCT1ReadPreProc); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// 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 : "Reads a CT1 preproc data 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 +// +Int_t MCT1ReadPreProc::AddFile(const char *txt, int) +{ + const char *name = gSystem->ExpandPathName(txt); + + TString fname(name); + delete [] name; + + if (!CheckHeader(fname)) + { + *fLog << warn << "WARNING - Problem reading header... ignored." << endl; + return 0; + } + + const Int_t n = GetNumEvents(fname); + if (n==0) + { + *fLog << warn << "WARNING - File contains no data... ignored." << endl; + return 0; + } + + fEntries += n; + + *fLog << inf << "File " << txt << " contains " << n << " events (Total=" << fEntries << ")" << endl; + + fFileNames->AddLast(new TNamed(txt, "")); + return 1; +} + +// -------------------------------------------------------------------------- +// +// 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) + fRawRunHeader->SetRunNumber(outpars.irunnum); + fRawRunHeader->SetReadyToSave(); + + // 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; iSetReadyToSave(); + + /* 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; + fRawRunHeader->SetRunTime(outpars.istart_mjdate_day, outpars.iend_mjdate_day); + fRawRunHeader->SetReadyToSave(); + + /* + --- 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->Init(*fGeom); + + fPedRMS.Set(iMAXNUMPIX); + + *fLog << "PedestalRMS : "; + for (Int_t i=0; iSetReadyToSave(); + + // Used to communicate the mean over all pixels + // pedestal RMS into the Runs NTuple, as it might + // be used for e.g. quality cuts. + // float fpedrms_mean; + *fLog << "Pedestal RMS: " << outpars.fpedrms_mean << endl; + + // The average current over the active pixels + // for this run. This is done separately for + // ON and OF runs. + //float fcurrent_mean; + + // enum eERRORTOLERANCE {CAUTIOUS=0, GOODPHYSICS, TANK} eerrortolerance; + /* 0 == "cautious", exits on any reason (but tells in + * the .err file, + * 1 == "goodphysics", exits when physics could be affected + * by the error, + * 2 == "tank", never exits except on coredumps and when + * all files have been processed. Do not use such files for + * physics analysis! + * + * NOTE: the capital letter words are the enums, the small letter + * words must be used inside the parameter file. */ + + // enum eMCTRIGGERFLAG {ALL=0, FLAG, NOFLAG} emctriggerflag; + /* all: all events which survive the filter are written to the + * events NTuple. + * flag: When Dorota's triggerflag is set to 1 for a particular + * event, it shall be written to the output. All others shall + * just be disregarded. (Default) + * noflag: Opposite of 'flag': only events with triggerflag = 0 shall + * be treated further. */ + + *fLog << "Particle Id #" << outpars.imcparticle << endl; + *fLog << "Right Ascension: " << outpars.dsourcera_hours << "h" << endl; + *fLog << "Declination: " << outpars.dsourcedec_deg << "deg" << endl; + + // Next statement commented out because bmontecarlo was set wrongly + //fIsMcFile = outpars.bmontecarlo==TRUE; + fIsMcFile = (outpars.dsourcera_hours==0 && outpars.dsourcedec_deg==0 && + outpars.imcparticle != 0); + + if (fIsMcFile != (outpars.bmontecarlo==TRUE)) + { + *fLog << "File tells you that it is a "; + *fLog << (outpars.bmontecarlo ? "Monte Carlo" : "Real Data"); + *fLog << " file." << endl; + } + + *fLog << "File detected as a "; + *fLog << (fIsMcFile ? "Monte Carlo" : "Real Data"); + *fLog << " file." << endl; + *fLog << " " << endl; +} + +// -------------------------------------------------------------------------- +// +// Read CT1 PreProc File Header: +// +Int_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 (fpreprocversion<0.6) + { + *fLog << err << "Sorry, only files from PreProc V0.6 and newer are supported." << endl; + return kFALSE; + } + + // + // This is a stupid way of getting rid of numerical uncertanties when + // comparing floating point numbers (Argh...) + // + TString s1 = Form("%.2f", structversion); + TString s2 = Form("%.2f", STRUCT_VERSION); + + if (s1 != s2) + { + *fLog << warn << "WARNING: Version of C-structures of file (V"; + *fLog << s1 << ") not identical with current structures (V"; + *fLog << s2 << ")" << endl; + } + + fIn->read((char*)&outpars, sizeof(struct outputpars)); + + ProcessRunHeader(outpars); + + //rwagner: ReInit whenever new run commences + // rc==-1 means: ReInit didn't work out + + MTaskList *tlist = (MTaskList*)fParList->FindCreateObj("MTaskList"); + if (!tlist) + return -1; + + if (!tlist->ReInit(fParList)) + return -1; + + return kTRUE; +} + +Int_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 0; + } + + *fLog << inf << cheadertitle << flush; + + struct filterresults filterres; + fIn->read((char*)&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 + */ + + if (fNumEventsInRun!=(UInt_t)filterres.ifilter_passed_evts) + { + *fLog << err << "ERROR! Number of events in run (" << (UInt_t)filterres.ifilter_passed_evts; + *fLog << ") doesn't match number of read events ("; + *fLog << fNumEventsInRun << ")" << endl; + *fLog << " File corrupted." << endl; + return -1; + } + + fNumFilterEvts += fNumEventsInRun; + fNumRuns++; + + *fLog << inf << "Read " << fNumEventsInRun << " events from run (Total="; + *fLog << fNumFilterEvts << "/" << fEntries << " ["; + *fLog << 100*fNumFilterEvts/fEntries << "%], Runs=" << fNumRuns << ")"; + *fLog << endl; + + return 1; +} + +// -------------------------------------------------------------------------- +// +// 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 existence of a next file to read + // + if (fNumFile >= (UInt_t)fFileNames->GetSize()) + return kFALSE; + + TNamed *file = (TNamed*)fFileNames->At(fNumFile); + + //TNamed *file = (TNamed*)fFileNames->GetFirst(); + //if (!file) + // return kFALSE; + + // + // open the file which is the first one in the chain + // + const TString name = file->GetName(); + + const char *expname = gSystem->ExpandPathName(name); + const TString fname(expname); + delete [] expname; + + // + // Remove this file from the list of pending files + // + //fFileNames->Remove(file); + + *fLog << inf << "Open file: '" << name << "'" << endl; + + if (!CheckHeader(fname)) + return kFALSE; + + fNumFile++; + + fIn = new ifstream(fname); + if (!*fIn) + { + *fLog << err << "Cannot open file " << fname << ": "; + *fLog << strerror(errno) << endl; + return kFALSE; + } + + *fLog << inf << "-----------------------------------------------------------------------" << endl; + + switch (ReadRunHeader()) + { + case kFALSE: + *fLog << warn << "Unable to read first run header... skipping file." << endl; + return kFALSE; + case -1: + *fLog << warn << "ReInit of Tasklist didn't succeed." << endl; + return kFALSE; + default: + *fLog << "After opening next file: Number of Events #" << fNumEventsInRun << endl; + return kTRUE; + } +} + +Bool_t MCT1ReadPreProc::CheckHeader(const TString fname) const +{ + ifstream fin(fname); + if (!fin) + { + *fLog << dbginf << err << "ERROR - Cannot open file '" << fname << "'" << endl; + return kFALSE; + } + + char cheadertitle[iHEADERTITLELENGTH]; + fin.read(cheadertitle, iHEADERTITLELENGTH); + + Float_t fpreprocversion, structversion; + Int_t dummyi; + + sscanf(cheadertitle, cTITLE_TEMPLATE, + &fpreprocversion, &structversion, + &dummyi, &dummyi, &dummyi, &dummyi); + + if (fpreprocversion < 0.6) + { + *fLog << dbginf << err << "ERROR - You must use PreProc V0.6 or higher." << endl; + return kFALSE; + } + + if (STRUCT_VERSION > structversion) + { + *fLog << warn << "WARNING: Version of C-structures of file (V"; + *fLog << structversion << ") newer than current structures (V"; + *fLog << STRUCT_VERSION << ")" << endl; + } + + *fLog << "Current structures: " << STRUCT_VERSION << " "; + *fLog << "Structures in file: " << structversion << " "; + *fLog << "Used preproc version: " << fpreprocversion << endl; + + return kTRUE; +} + + +Int_t MCT1ReadPreProc::GetNumEvents(const TString fname) const +{ + *fLog << inf << "Scanning file " << fname << " for size" << flush; + + ifstream fin(fname); + if (!fin) + { + *fLog << dbginf << err << "ERROR - Opening file." << endl; + return 0; + } + + 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.peek()!=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; +} + +Bool_t MCT1ReadPreProc::Rewind() +{ + fNumFilterEvts = 0; + fNumEvents = 0; + fNumRuns = 0; + fNumFile = 0; + if (fIn) + delete fIn; + + fIn=NULL; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// 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 MPedPhotCam container to 127 pixels (CT1 camera) +// +Int_t MCT1ReadPreProc::PreProcess(MParList *pList) +{ + fParList = pList; + + // + // look for the HourAngle container in the plist + // + fHourAngle = (MParameterD*)pList->FindCreateObj("MParameterD", "HourAngle"); + if (!fHourAngle) + return kFALSE; + fHourAngle->SetTitle("Store the CT1 hour angle [deg]"); + + // + // look for the ThetaOrig container in the plist + // + fThetaOrig = (MParameterD*)pList->FindCreateObj("MParameterD", "ThetaOrig"); + if (!fThetaOrig) + return kFALSE; + fThetaOrig->SetTitle("Store the original CT1 zenith angle [rad]"); + + // + // 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 = (MPedPhotCam*)pList->FindCreateObj("MPedPhotCam"); + if (!fPedest) + return kFALSE; + + // + // look for the time class in the plist + // + fTime = (MTime*)pList->FindCreateObj("MTime"); + if (!fTime) + return kFALSE; + + // + // look for the pedestal class in the plist + // + fBlinds = (MBadPixelsCam*)pList->FindCreateObj("MBadPixelsCam"); + if (!fBlinds) + return kFALSE; + + // + // look for the source position in the camera + // + fSrcPos = (MSrcPosCam*)pList->FindCreateObj("MSrcPosCam"); + 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; + + // + // look for the raw run header class + // + fRawRunHeader = (MRawRunHeader*)pList->FindCreateObj("MRawRunHeader"); + if (!fRawRunHeader) + return kFALSE; + + fBinningT = (MBinning*)pList->FindObject("BinningTheta"); + + Rewind(); + + fPedest->Init(*fGeom); + + return GetSelector() ? GetSelector()->CallPreProcess(pList) : kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Smear Theta uniformly in a bin of Theta +// theta [rad] +// SmearTheta [rad] +// +Double_t MCT1ReadPreProc::SmearTheta(Double_t theta) +{ + if (!fBinningT) + return theta; + + const Int_t bin = fBinningT->FindLoEdge(theta * 180/TMath::Pi()); + if (bin<0) + return theta; + + // smear Theta within the Theta bin + const Double_t low = fBinningT->GetEdges()[bin]; + const Double_t up = fBinningT->GetEdges()[bin+1]; + + // "up-": Do not allow the upper edge + return (up - gRandom->Uniform() * (up-low)) * TMath::Pi()/180; +} + +// -------------------------------------------------------------------------- +// +// Discretize Theta according to the binning in Theta +// theta [rad] +// DiscreteTheta [rad] (= bin center) +// +Double_t MCT1ReadPreProc::DiscreteTheta(Double_t theta) +{ + if (!fBinningT) + return theta; + + const Int_t bin = fBinningT->FindLoEdge(theta * 180/TMath::Pi()); + if (bin<0) + return theta; + + const Double_t low = fBinningT->GetEdges()[bin]; + const Double_t up = fBinningT->GetEdges()[bin+1]; + + return 0.5*(up+low) * TMath::Pi()/180; +} + +// -------------------------------------------------------------------------- +// +// Analyse the event data, means store and copy the needed data into +// Mars structures and data members +// +Bool_t MCT1ReadPreProc::ProcessEvent(const struct eventrecord &event) +{ + /* + if (fRawRunHeader->GetRunNumber() == 1) + { + *fLog << "eventrecord" << endl; + *fLog << "isecs_since_midday = " << event.isecs_since_midday << endl; + *fLog << "isecfrac_200ns = " << event.isecfrac_200ns << endl; + *fLog << "snot_ok_flags = " << event.snot_ok_flags << endl; + *fLog << "ialt_arcs = " << event.ialt_arcs << endl; + *fLog << "iaz_arcs = " << event.iaz_arcs << endl; + *fLog << "ipreproc_alt_arcs = " << event.ipreproc_alt_arcs << endl; + *fLog << "ipreproc_az_arcs = " << event.ipreproc_az_arcs << endl; + *fLog << "ifieldrot_arcs = " << event.ifieldrot_arcs << endl; + + *fLog << "srate_millihz = " << event.srate_millihz << endl; + *fLog << "fhourangle = " << event.fhourangle << endl; + *fLog << "fmcenergy_tev = " << event.fmcenergy_tev << endl; + *fLog << "fmcsize_phel = " << event.fmcsize_phel << endl; + *fLog << "imcimpact_m = " << event.imcimpact_m << endl; + *fLog << "imcparticle = " << event.imcparticle << endl; + *fLog << "imctriggerflag = " << event.imctriggerflag << endl; + } + */ + + // reset blind pixels for this event + fBlinds->Clear(); + for (int i=0; iSetCT1Time((UInt_t)fRawRunHeader->GetRunStart().GetMjd(), event.isecfrac_200ns, event.isecs_since_midday); + fTime->SetReadyToSave(); + + /* + --- USEFULL? NEEDED? --- + 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 + // 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].GetRms()); + } + fNphot->FixSize(); + fNphot->SetReadyToSave(); + + // int ipreproc_alt_arcs; // "should be" alt according to preproc (arcseconds) + // int ipreproc_az_arcs; // "should be" az according to preproc (arcseconds) + + // smear Theta in its Theta bin + const Double_t theta = TMath::Pi()*(0.5-1./180*event.ialt_arcs/3600); + fThetaOrig->SetVal(theta); + + // store hour angle + fHourAngle->SetVal(event.fhourangle); + + fMcEvt->Fill(event.isecs_since_midday, //0, /*fEvtNum*/ + fIsMcFile ? event.imcparticle : 0, /*corsika particle type*/ + fIsMcFile ? event.fmcenergy_tev*1000 : 0, + 0, /* fThi0 */ + 0, /* fFirTar */ + 0, /* fzFirInt */ + // 0, /* fThet*/ + // rwagner: The following should be theta, right? Check with + // altitude fill some lines down... + 0, // altitude (arcseconds) + 0, /* fPhii */ + 0, /* fCorD */ + 0, /* fCorX */ + 0, /* fCorY */ + fIsMcFile ? event.imcimpact_m*100 : 0, + TMath::Pi()/180*event.iaz_arcs/3600, // azimuth (arcseconds) + // fIsMcFile ? SmearTheta(theta) : theta, + DiscreteTheta(theta), + 0, /* fTFirst */ + 0, /* fTLast */ + 0, /* fL_Nmax */ + 0, /* fL_t0 */ + 0, /* fL_tmax */ + 0, /* fL_a */ + 0, /* fL_b */ + 0, /* fL_c */ + 0, /* fL_chi2 */ + 0, /* uiPin */ + 0, /* uiPat */ + 0, /* uiPre */ + 0, /* uiPco */ + 0, /* uiPelS */ + (int)(fIsMcFile ? event.fmcsize_phel : 0), /* uiPelC, Simulated SIZE */ + 0, /* elec */ + 0, /* muon */ + 0 /* other */ + ); + + fMcTrig->SetFirstLevel(event.imctriggerflag); // MC data from Dorota get a triggerflag: 1 means triggered, 0 not. */ + + fMcTrig->SetReadyToSave(); + fMcEvt->SetReadyToSave(); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// 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 this isn't successful... + // must be an event + // + switch (ReadRunFooter()) + { + case -1: + return kFALSE; + case 0: + return kTRUE; + } + + *fLog << inf << "Footer found." << endl; + + const char c = fIn->peek(); + + // + // No after reading the footer check if we reached the end of the file + // + if (fIn->eof() || c==EOF) + { + *fLog << "End of file." << endl; + return kFALSE; + } + + // + // If the eof isn't reached a new header must follow. Check for it. + // + if (c!=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() < 0) + { + *fLog << warn << "ReInit of Tasklist didn't succeed." << endl; + return kFALSE; + } + } +} + +// -------------------------------------------------------------------------- +// +// 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. +// +Int_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; + + // + // Check for a selector. If one is given and returns kFALSE + // skip this event. + // + if (GetSelector()) + { + // + // Make sure selector is processed + // + if (!GetSelector()->CallProcess()) + { + *fLog << err << dbginf << "Processing Selector failed." << endl; + return kFALSE; + } + + // + // Skip Event + // + if (!GetSelector()->IsConditionTrue()) + { + fIn->seekg(sizeof(struct eventrecord), ios::cur); + + fNumEvents++; + fNumEventsInRun++; + + return kCONTINUE; + } + } + + // event data to be read from the file + struct eventrecord event; + + + + // read the eventrecord from the file + fIn->read((char*)&event, sizeof(struct eventrecord)); + + switch (ProcessEvent(event)) + { + case kFALSE: + return kFALSE; + case kCONTINUE: + return kCONTINUE; + } + + fNumEvents++; + fNumEventsInRun++; + + return kTRUE; +} + +Int_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; + *fLog << "Number of events detected first: " << fEntries << 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; + } + + delete fIn; + fIn = NULL; + + return GetSelector() ? GetSelector()->CallPostProcess() : kTRUE; +} Index: trunk/MagicSoft/Mars/manalysisct1/MCT1ReadPreProc.h =================================================================== --- trunk/MagicSoft/Mars/manalysisct1/MCT1ReadPreProc.h (revision 4460) +++ trunk/MagicSoft/Mars/manalysisct1/MCT1ReadPreProc.h (revision 4460) @@ -0,0 +1,104 @@ +#ifndef MARS_MCT1ReadPreProc +#define MARS_MCT1ReadPreProc + +#ifndef ROOT_TArrayF +#include +#endif + +#ifndef ROOT_TArrayI +#include +#endif + +#ifndef MARS_MRead +#include "MRead.h" +#endif + +class TList; +class MTime; +class MMcEvt; +class MMcTrig; +class MGeomCam; +class MSrcPosCam; +class MCerPhotEvt; +class MPedPhotCam; +class MBadPixelsCam; +class MRawRunHeader; +class MTaskList; +class MParList; +class MParameterD; +class MBinning; + +struct outputpars; +struct eventrecord; + +class MCT1ReadPreProc : public MRead +{ +private: + ifstream *fIn; // the inputfile + TList *fFileNames; // Array which stores the \0-terminated filenames + + MGeomCam *fGeom; // camera geometry + MCerPhotEvt *fNphot; // the data container for all data. + MPedPhotCam *fPedest; // ct1 pedestals + MTime *fTime; // event time + MMcEvt *fMcEvt; // monte carlo data container for MC files + MMcTrig *fMcTrig; // mc data container for trigger information + MSrcPosCam *fSrcPos; // source position in the camera + MBadPixelsCam *fBlinds; // Array holding blind pixels + MRawRunHeader *fRawRunHeader; // raw run header + MParList *fParList; // parameter list + MParameterD *fHourAngle; // hour angle [deg] + MParameterD *fThetaOrig; // original zenith angle [rad] + MBinning *fBinningT; // Theta binning for the smearing + + Bool_t fIsMcFile; // Flag whether current run is a MC run + + UInt_t fNumFile; + UInt_t fNumEvents; // number of events counted in all runs in all files + UInt_t fNumEventsInRun; // number of events in the counted in th ecurrent run + UInt_t fNumRuns; // number of processed runs of all files + UInt_t fEntries; // entries of all files succesfully added + UInt_t fNumFilterEvts; // number of events mentioned in the runs footers + + TArrayF fPedRMS; + TArrayI fBlnd; + + Bool_t OpenNextFile(); + + Int_t GetNumEvents(const TString name) const; + Bool_t CheckHeader(const TString fname) const; + + void ReadPedestals(); + Int_t ReadRunHeader(); + Int_t ReadRunFooter(); + Bool_t CheckFilePosition(); + void ProcessRunHeader(const struct outputpars &outpars); + Bool_t ProcessEvent(const struct eventrecord &event); + + Double_t SmearTheta(Double_t theta); + Double_t DiscreteTheta(Double_t theta); + + Int_t PreProcess(MParList *pList); + Int_t Process(); + Int_t PostProcess(); + + Bool_t Rewind(); + +public: + MCT1ReadPreProc(const char *filename=NULL, + const char *name=NULL, + const char *title=NULL); + + ~MCT1ReadPreProc(); + + Int_t AddFile(const char *fname, int i=0); + + UInt_t GetEntries() { return fEntries; } + + ClassDef(MCT1ReadPreProc, 0) // Reads the CT1 preproc data file +}; + +#endif + + + Index: trunk/MagicSoft/Mars/manalysisct1/Makefile =================================================================== --- trunk/MagicSoft/Mars/manalysisct1/Makefile (revision 4459) +++ trunk/MagicSoft/Mars/manalysisct1/Makefile (revision 4460) @@ -27,5 +27,7 @@ MCT1PointingCorrCalc.cc \ MCT1Supercuts.cc \ - MCT1SupercutsCalc + MCT1SupercutsCalc \ + MCT1ReadPreProc.cc \ + MCT1ReadAscii.cc ############################################################ Index: trunk/MagicSoft/Mars/manalysisct1/structures.h =================================================================== --- trunk/MagicSoft/Mars/manalysisct1/structures.h (revision 4460) +++ trunk/MagicSoft/Mars/manalysisct1/structures.h (revision 4460) @@ -0,0 +1,882 @@ +/******************************************************************/ +/* GGG AAA M M M M AAA SSS */ +/* G G A A MM MM MM MM A A S S */ +/* G A A M M M M M M A A SS */ +/* G GG AAAAA M M M M AAAAA SS */ +/* G G A A M M M M A A S S */ +/* GGG A A M M M M A A SSS */ +/* Gamma Astronomer's Munich Madrid Analysis Suite */ +/* An Atmospheric Cherenkov Telescope Data Analysis Software */ +/* MPI f"ur Physik, M"unchen & Universidad Complutense, Madrid */ +/******************************************************************/ + +/******************************************************************/ +/* This File belongs to the programs */ +/* */ +/* P R E P R O C E S S O R and I M A G E R */ +/* */ +/* Purpose: provide the structure definitions common to both */ +/* programs, especially the output structures of the */ +/* preprocessor */ +/* */ +/******************************************************************/ + +/******************************************************************/ +/* @(#) File name structures.h + @(#) latest SID 0.5 + @(#) Date of latest delta 26/6/00 + Author: D. Petry + changed by M. Kestel */ +/******************************************************************/ + +#ifndef _structures_h_ +#define _structures_h_ "__structures_h__" + +#define STRUCT_VERSION 0.6 + +#define iHEADERTITLELENGTH 60 +/* length of the string written at the beginning of the */ +/* output for each runfile */ +#define cTITLE_TEMPLATE "PREPROC V%f/S%f CT %d RUN %d %d PROCMJD %d\n" + +#define cEND_EVENTS_TEMPLATE "EVENTS AFTER FILTER = %d\n" + +#define iMAXNUMPIX 127 +/* maximum value for the number of camera pixels (excl. muon shield) */ +#define iMAXNUMMUONPIX 5 +/* maximum value for the number of pixels in the muon shield */ +#define iMAXNUMCOINTDCS 5 +/* maximum value for the number of coincidence tdcs */ + + +#define iBADMACTIME_BIT 1 +/* snot_ok_flags is ORed with this number to set the bit which + * indicates that the mactime measurement didn't work */ +#define iBADRUBTIME_BIT 2 +/* snot_ok_flags is ORed with this number to set the bit which + * indicates that the Rubidium time measurement didn't work */ +#define iBADPOINTING_BIT 4 +/* snot_ok_flags is ORed with this number to set the bit which + * indicates that the shaft encoder values were inconsistent with + * time OR source coordinates */ +#define iMUON_BIT 8 + /* snot_ok_flags is ORed with this number to set the bit which + * indicates that the event didn't pass the muon filter + * (for test puposes) */ + +#define i2OUTFAIL_BIT 16 +/* snot_ok_flags is ORed with this number to set the bit which + * indicates that the event didn't pass the "two out of all" + * software trigger */ +#define iSUMFAIL_BIT 32 +/* snot_ok_flags is ORed with this number to set the bit which + * indicates that the event didn't pass the "sum of all + * photoelectrons" filter */ + +#define iPARTICLETYPE_FACTOR 1024 +/* for MC data: the particle type multiplied with this number + * is ORed with snot_ok_flags */ + +enum onoroff {NEITHER_ON_NOR_OFF, OFF_SOURCE, ON_SOURCE}; +enum trackmodes {NORMAL, REVERSE, DUNNO}; +enum axes {RA, DEC, ALT, AZ}; +enum eERRORTOLERANCE {CAUTIOUS=0, GOODPHYSICS, TANK}; +enum eMCTRIGGERFLAG {ALL=0, FLAG, NOFLAG}; + +struct outputpars { + /* preprocessing parameters that are passed to imager. + * Note: Not all parameters from the parameter file are passed to + * imager, only those that are repeated in this struct outputpars + * here !! That's why these two structs must be the same for both + * preproc and imager. */ + + int itelescope; + /* number of the CT which took the data (from parameters file */ + /* crosschecked with run file name) */ + float flongitude_deg; + /* longitude (counted positive towards West) of CT position */ + float flatitude_deg; + /* latitude (counted positive towards North) of CT position */ + int irunnum; + /* run number (from parameters file cross checked with run */ + /* file name) */ + enum onoroff eruntype; +#ifdef LINUX + int dummy; +#endif + /* indicates if the */ + /* run is on- or off source (from runfile cross checked with */ + /* file name) */ + double dsourcera_hours; + /* right ascension of observed source in hours */ + double dsourcedec_deg; + /* declination of observed source in degrees */ + int inumpixels; + /* number of pixels in the camera (from parameters file cross */ + /* checked with run file, calibration file and pedestal file) */ + int inummuonpixels; + /* number of pixels in the muon shield (from parameters file */ + /* cross checked with run file, calibration file and pedestal */ + /* file) */ + int inumcointdcs; + /* number of coincidence tdcs recorded in the runfile (from */ + /* parameters file cross checked */ + /* with run file) */ + float fpixdiameter_deg; + /* smallest pixel diameter (degrees) (from parameters file) */ + enum axes ese1_is; + /* name of the axis to which shaft encoder 1 is attached */ + /* (implies the type of mount) */ + int isezeropos[2]; + /* zero position of shaftencoders 1 and 2 (from parameters file) */ + 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 */ + float fbendingcorr; + /* bending correction factor (ALT/AZ mount only) */ + float fextinction; + /* atmospheric extinction (typically taken from the Carlsberg */ + /* Meridian Circle data) */ + Boolean bdontusepix[iMAXNUMPIX]; + /* bdontusepix is set true if the pixel should due to whatever + * reason not be used in image analysis, otherwise it is true; + * + * bdontusepix is set by the user inside preproc.par and is + * untouched by preproc (version >= 0.4). */ + Boolean bexcludepix[iMAXNUMPIX]; + /* 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 */ + int irubminusutc_usecs; + /* difference between rubidium clock and UTC in microseconds */ + int isum_thresh_phot; + /* threshold for the total sum of photoelectrons filter */ + /* from the parameters file */ + int i2out_thresh_phot; + /* threshold for the two-pixels-out-of-all software */ + /* trigger from parameters file */ + int imuoncut_thresh_adccnt[iMAXNUMMUONPIX]; + /* thresholds for the muon cut */ + Boolean bmuon_suppression; + /* "Dolby" noise reduction flag */ + float ftolerated_pointerror_deg; + /* maximum tolerated pointing error in the position + * check in iprocessdata */ + float fxpointcorr_deg; + /* pointing correction (to be added along the camera + * x axis) e.g. from point run */ + float fypointcorr_deg; + /* pointing correction (to be added along the camera + * y axis) e.g. from point run */ + 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. + * This is given once per run. */ + + /* 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) */ + int iproc_mjdate; + /* MJD of data processing (i.e. creation of this file) */ + enum trackmodes etrackmode; + /* tracking mode (for ALT/AZ CTs) */ + int iproc_evts; +#ifdef LINUX + int dummy2; +#endif + /* number of events processed */ + 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 */ + float fpedrms_mean; /* Used to communicate the mean over all pixels + * pedestal RMS into the Runs NTuple, as it might + * be used for e.g. quality cuts. */ + float fcurrent_mean; /* The average current over the active pixels + * for this run. This is done separately for + * ON and OF runs. */ + + enum eERRORTOLERANCE eerrortolerance; + /* 0 == "cautious", exits on any reason (but tells in + * the .err file, + * 1 == "goodphysics", exits when physics could be affected + * by the error, + * 2 == "tank", never exits except on coredumps and when + * all files have been processed. Do not use such files for + * physics analysis! + * + * NOTE: the capital letter words are the enums, the small letter + * words must be used inside the parameter file. */ + enum eMCTRIGGERFLAG emctriggerflag; + /* all: all events which survive the filter are written to the + * events NTuple. + * flag: When Dorota's triggerflag is set to 1 for a particular + * event, it shall be written to the output. All others shall + * just be disregarded. (Default) + * noflag: Opposite of 'flag': only events with triggerflag = 0 shall + * be treated further. */ + Boolean bmontecarlo; + /* if TRUE we have a monte carlo dataset before us, if FALSE + * (default value) it is just normal data. */ +}; + +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 */ + +}; + +struct eventrecord { + 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 */ + int ialt_arcs; + /* altitude (arcseconds) */ + int iaz_arcs; + /* azimuth (arcseconds) */ + int ipreproc_alt_arcs; + /* "should be" alt according to preproc (arcseconds) */ + int ipreproc_az_arcs; + /* "should be" az according to preproc (arcseconds) */ + int ifieldrot_arcs; + /* 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 */ + unsigned short srate_millihz; + /* 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 */ + float fhourangle; /* 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 fmcenergy_tev; /* Simulated Energy.... dropping through to + * the Events NTuple. */ + float fmcsize_phel; /* Simulated SIZE.... dropping through to + * the Events NTuple. */ + int imcimpact_m; + /* MC data contain the impact parameter, which is given here in + * meters. */ + int imcparticle; + /* MC data know which particle they are.... all in CORSIKA standard. */ + int imctriggerflag; + /* MC data from Dorota get a triggerflag: 1 means triggered, 0 not. */ + short spixsig_10thphot[iMAXNUMPIX]; + /* number of photoelectrons measured in each pixel only the + * actual number of pixels (outputpars.inumpixels) is written out */ +}; + +struct camera { /* camera parameters for imaging */ + int inumpixels; + int inumrings; + double dpixdiameter_deg; + double dxc[iMAXNUMPIX]; +/* Pixel coordinates in camera coordinate system (x points from + * pixel 1 to 2). */ + double dyc[iMAXNUMPIX]; + /* The numbering of the pixels in these arrays starts at 0! */ + double dxpointcorr_deg; + /* correction of the pixel coordinates; to be added to dxc[] + * to get correct value */ + double dypointcorr_deg; + /* correction of the pixel coordinates; to be added to dxc[] + * to get correct value */ +}; + + +/* two structures for better file handling */ +struct inputfile { + char cname[iMAXFILENAMELENGTH]; + /* filename (including path) */ + FILE *pointer; + /* filepointer */ + char ccannotopentext[161]; + /* Error text printed when file cannot be opened */ + int icannotopencode; + /* Error code for the exit statement when file could not be opened */ + int iearlyeofcode; + /* Error code for the exit statement for unexpected EOF */ + int ierroratclosecode; + /* Error code for the exit statement for error while trying to close */ +}; + +struct outputfile { + char cname[iMAXFILENAMELENGTH]; /* filename (including path) */ + FILE *pointer; /* filepointer */ + char ccannotopentext[161]; + /* Error text printed when file cannot be opened */ + int icannotopencode; + /* Error code for the exit statement when file could not be opened */ + int icannotwritecode; + /* Error code for the exit statement for failed fprintf */ + int ierroratclosecode; + /* Error code for the exit statement for error while trying to close */ +}; + +struct P2ALLHISTOS { + + HISTO *p2adchist; + HISTO *p2pedonhist; + HISTO *p2pedofhist; + HISTO *p2pedonmeanhist; + HISTO *p2pedofmeanhist; + HISTO *p2pedonjitthist; + HISTO *p2pedofjitthist; + HISTO *p2calonhist; + HISTO *p2calofhist; + HISTO *p2calonmeanhist; + HISTO *p2calofmeanhist; + HISTO *p2calonjitthist; + HISTO *p2calofjitthist; + HISTO *p2phehist; + HISTO *p2zenonhist; + HISTO *p2zenofhist; + HISTO *p2corrpedhist; + HISTO *p2pedfwhmhist; +}; + +struct CALNT_DAT { + /* the structure for an entry in the calibration ntuple */ + float CTrunnum; + float runtype; + float ALT_deg; /* the ALT taken from the first event of the run */ + float extinct; /* extinction */ + float RtiSecBF; /* runtime in seconds before filter */ + float pixnum; /* number of the pixel for which the following + * data is given */ + float pixonoff; /* Quantity telling whether this pixel was switched + * off by the user (-= 2.) or by preproc (-= 1.). + * If ON, its value = 0. */ + float ped_adcc; /* the standard deviation of the pedestal as + * calculated by TIJARAFE */ + float pedrmsac; /* the standard deviation of the pedestal in + * ADC counts as calculated by TIJARAFE */ + float pedrmsPE; /* the calibrated ped RMS in photoelectrons + * which is also given to the imager */ + float las_adcc; /* the mean signal for this pixel from the + * laser run */ + float lasrmsac; /* RMS of the laser events after cleaning from cosmics */ + float conv_fac; /* conversion factor = /Var(Q) * F^2 */ + float F; /* The F from the line before */ + float rel_gain; /* the relative gain of the pixel */ + float numtrig; /* number of events in which the pixel was + * above trigger threshold */ + float num3sig; /* number of events in which the pixel was + * 3 sigma above its pedestal */ + float HV1; /* Adjustment of HV1 for this run, if available. */ + float HV2; /* Adjustment of HV2 for this run, if available. */ + float ThrCurr; + float AvgCurr; /* the sum of all currents during this run (with the + * particular runtype in question) from the TelLog + * file. The currents are averaged over the 2-minute + * intervals but summed over all pixels with currents + * higher than the threshold current value. This sum + * is then stored in the runs ntuple. + * Same for the scaler values --> */ + float AvgScal; +}; + +struct inputpars { + Boolean bkeywordgiven[iNUMKEYWORDS+1]; /* if a valid keyword + * is given in the + * parameters file + * the corresponding + * element in this + * array is set to true */ + int itelescope; /* number of the CT which took the data */ + float flongitude_deg; /* longitude (counted positive + * towards West) of CT position */ + float flatitude_deg; /* latitude (counted positive + * towards North) of CT position */ + int ifirstrunnum; /* first run number to which these + * parameters apply*/ + int ilastrunnum; /* last run number to which these + * parameters apply, */ + /* i.e. this run will be processed + * with these parameters */ + int inumpixels; /* number of pixels in the camera */ + 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) */ + enum axes ese1_is; /* name of the axis to which shaft + * encoder 1 is attached + * (implies the type of mount) (the + * type axes is declared in structures.h) */ + float fdegrees_per_step[2]; /* angular resolution of shaft + * encoders 1 and 2 */ + int isezeropos[2]; /* zero position of shaftencoders 1 and + * 2 from parameters file */ + 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 */ + float fbendingcorr; /* bending correction factor + * (ALT/AZ mount only) */ + Boolean bdontusepix[iMAXNUMPIX]; /* bdontusepix is set true + * if the pixel should due + * to whatever reason not be + * used in image analysis, + * otherwise it is false; + * this is a copy of the + * input from the parameters file */ + Boolean bdontusepix_in_trig[iMAXNUMPIX]; /* is set true if the + * pixel should due + * to whatever reason not + * be used in the two out + * of all trigger, otherwise + * it is false; this is a + * copy of the input from + * the parameters file */ + float fphotoel_per_adccnt[iMAXNUMPIX]; /* conversion factors for + * the pixel signals */ + float fextinction; /* atmospheric extinction (typically + * taken from the Carlsberg Meridian + * Circle data) */ + int irubminusutc_usecs; /* difference between rubidium clock + * and UTC in microseconds */ + int isum_thresh_phot; /* threshold for the total sum of + * photoelectrons filter from the + * parameters file */ + int i2out_thresh_phot; /* threshold for the + * two-pixels-out-of-all software + * trigger from parameters file */ + int imuoncut_thresh_adccnt[iMAXNUMMUONPIX]; /* thresholds for + * the muon cut */ + Boolean bmuon_suppression; /* if true, the events which do + * not pass the muon cut are not + * written to the output */ + float ftolerated_pointerror_deg; /* maximum tolerated pointing + * error in the position check + * in iprocessdata */ + float fxpointcorr_deg; /* pointing correction (to be added + * along the camera x axis) e.g. + * from point run */ + float fypointcorr_deg; /* pointing correction (to be added + * along the camera y axis) e.g. + * from point run */ + 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) */ + enum pedsearchdirs {MATCH=0, BEFORE, AFTER, PINGPONG, + PEAKFWHM, NONE} epedsearchdir ; + /* MATCH = only same number as run file (if ipedendofsearch is 0) or + exactly ipedendofsearch; + BEFORE = search backwards until pedendofsearch; + AFTER = search forward until pedendofsearch, + PINGPONG = Try to fibnd the matching partner, try run#-1, run#+1, + run#-2, run#+2 etc. up to the third argument, the partnerfindrange + NONE = use none of both. This is only useful for Monte Carlo + Analysis, for everything else it is highly unrecommended.*/ + int ipedendofsearch; /* pedestal file number until which to + * search (see epedsearchdir) */ + int ipedsearchdepth; /* depth, until which preproc should + * try to find TelPed* partnerfiles */ + enum pedsearchdirs ecalsearchdir; + /* MATCH = only same number as run file (if ipedendofsearch is 0) or + exactly ipedendofsearch; + BEFORE = search backwards until pedendofsearch; + AFTER = search forward until pedendofsearch, + PINGPONG = Try to fibnd the matching partner, try run#-1, run#+1, + run#-2, run#+2 etc. up to the third argument, the partnerfindrange + NONE = use none of both. This is only useful for Monte Carlo + Analysis, for everything else it is highly unrecommended.*/ + int icalendofsearch; /* calibration file number until which + * to search (see ecalsearchdir)*/ + int icalsearchdepth; /* depth, until which preproc should + * try to find TelCal* partnerfiles */ + double dsourcera_hours; /* right ascension of observed source + * in hours */ + double dsourcedec_deg; /* declination of observed source in + * degrees */ + float fpadlevel_phot_pix[iMAXNUMPIX]; /* if > 0., software + * padding to this + * level is applied + * for the given pixel */ + float fseed1; /* seed for random generators */ + float fseed2; + + int iswap_partner_of_pix[iMAXNUMPIX]; /* if >= 0, pixel i is + * swapped with + * iswap_partner_of_pix[i] + * in ireadrunfile, -pedfile, + * -calfile (internal + * numbering starting at 0)*/ + + Boolean bnewhistograms; /* set to true if new histograms shall be booked */ + Boolean bhisto_on; /* set to true if histogramming is on */ + + enum eERRORTOLERANCE eerrortolerance; + /* 0 == "cautious", exits on any reason (but tells in + * the .err file, + * 1 == "normal", exits when physics could be affected + * by the error, + * 2 == "tank", never exits except on coredumps and when + * all files have been processed. Do not use such files for + * physics analysis! + * + * NOTE: the capital letter words are the enums, the small letter + * words are to be used inside the parameter file. */ + Boolean bmontecarlo; + /* if TRUE we have a monte carlo dataset before us, if FALSE + * (default value) it is just normal data. */ + enum eMCTRIGGERFLAG emctriggerflag; + /* all: all events which survive the filter are written to the + * events NTuple. + * flag: When Dorota's triggerflag is set to 1 for a particular + * event, it shall be written to the output. All others shall + * just be disregarded. (Default) + * noflag: Opposite of 'flag': only events with triggerflag = 0 shall + * be treated further. */ + char mcfile[FILENAME_MAX]; + float fmean_nsb_pe_pix[iMAXNUMPIX]; + Boolean bfitmypedrms; + /* If TRUE, the pedestal MEAN and SIGMA=sqrt(MEAN) is fitted + * from pickup corrected pedestal data and not calculated from + * the pickup corrected ADC values, as it was done up to + * ~November 2000. */ +}; + +struct calibrationdata { /* Calibration (Laser Run) Data to be read + * from the Te?Cal files */ + int iid; /* no useful information, is long in the Macintosh format */ + int itelescope; /* number of the CT which took the data */ + int irunnum; /* run number from the file name */ + enum onoroff eruntype; /* indicates if the run is on- or off + * (taken from filename) */ + int inumruns; /* number of laser runs recorded in this file, is + * long in the Macintosh format */ + char cname[iMAXFILENAMELENGTH]; /* copy of the filename used when + * the data was read */ + double dmean_adc[iMAXNUMLASERRUNS][iMAXNUMPIX]; + /* means of pixel signals over all shots calculated by TIJARAFE / + * Roberto: in processdata.c this number is recalculated. */ + double dmean_lasoutput[iMAXNUMLASERRUNS][iMAXNUMLASEREVENTS]; + /* This new (April 26 2000) quantity monitors the laser mean + * output recorded in events, which are not considered cosmics. */ + double dmean_alladc[iMAXNUMPIX]; + /* means of pixel signals over all triggers calculated by + * preproc. */ + double dmean_laserjitter[iMAXNUMPIX]; + /* Here we calc. the deviation from the actual lasergain in a + * single laser shot to the laser mean as averaged over the camera */ + double drms_adc[iMAXNUMPEDRUNS][iMAXNUMPIX]; + /* the standard deviations of the pixel signals */ + int inumpixels; /* number of pixels in the camera, short in + * Macintosh format */ + int inumevents[iMAXNUMLASERRUNS]; /* number laser shots, + * short in Macintosh format */ + int ipixsig_adc[iMAXNUMLASERRUNS][iMAXNUMLASEREVENTS][iMAXNUMPIX]; + /* the signals obtained from each pixel for each laser shot and + * each laser run, unsigned short in the Macintosh format */ +}; + +struct pedestaldata { /* Pedestal Run Data to be read from the + * Te?Ped files */ + int iid; /* no useful information, is long in the Macintosh + * format */ + int itelescope; /* number of the CT which took the data */ + int irunnum; /* run number from the file name */ + enum onoroff eruntype; /* indicates if the run is on- or of + * (taken from filename) */ + int inumruns; /* number of pedestal events recorded in this file, + * is long in the Macintosh format */ + char cname[iMAXFILENAMELENGTH]; /* copy of the filename used + * when the data was read */ + double dmean_adc[iMAXNUMPEDRUNS][iMAXNUMPIX]; + /* means of pixel signals over all triggers calculated by + * TIJARAFE using a sigma cut of 2 */ + double dmean_alladc[iMAXNUMPIX]; + /* means of pixel signals over all triggers calculated by + * preproc. */ + double dmean_allPE[iMAXNUMPIX]; + /* means of pixel signals over all triggers calculated by + * preproc. */ + double dmean_pedlevel[iMAXNUMPEDRUNS][iMAXNUMPEDEVENTS]; + /* This new (preproc_0.4) quantity monitors the pedestal mean + * output recorded in events, which are not considered cosmics. + * The calc.ing is done inside processdata.c */ + double dmean_pedlevelPE[iMAXNUMPEDRUNS][iMAXNUMPEDEVENTS]; + /* This new (preproc_0.4) quantity monitors the pedestal mean + * output recorded in events, which are not considered cosmics. + * The calc.ing is done inside processdata.c */ + double dmean_pedoffset[iMAXNUMPIX]; + /* Mean Offset from the ADC spectrum of a pixel to the mean of the + * pedestal events. */ + double dmean_pedoffsetPE[iMAXNUMPIX]; + /* Mean Offset from the ADC spectrum of a pixel to the mean of the + * pedestal events. */ + double drms_adc[iMAXNUMPEDRUNS][iMAXNUMPIX]; + /* the standard deviations of the pixel signals used in the + * sigma cut mentioned above */ + double drms_alladc[iMAXNUMPIX]; + /* The RMS value of the as-is ADC spectrum for all events inside a + * Pedestal file. */ + double drms_allPE[iMAXNUMPIX]; + /* The RMS value of the as-is ADC spectrum for all events inside a + * Pedestal file. */ + double drms_nopickup[iMAXNUMPIX]; + /* the standard deviations of the pixel signals' deviation from + * the events' mean (in ADC counts, of course...) --- + * This is the same as the RMS of the pickup corrected pedestals! */ + double drms_nopickupPE[iMAXNUMPIX]; + /* the standard deviations of the pixel signals' deviation from + * the events' mean (in ADC counts, of course...) --- + * This is the same as the RMS of the pickup corrected pedestals! */ + int inumpixels; /* number of pixels in the camera, short in + * Macintosh format */ + int inumevents[iMAXNUMPEDRUNS]; /* number of random triggers, + * short in Macintosh format */ + int ipixsig_adc[iMAXNUMPEDRUNS][iMAXNUMPEDEVENTS][iMAXNUMPIX]; + /* the signals obtained from each pixel for each random trigger + * and each pedestal event, short in the Macintosh format */ +}; + +struct runfiledata { /* Data read from Te?Run files */ + int iformatid; /* Format ID of the data, long in Macintosh format */ + int iversid; /* Version ID of the data, long in Macintosh format */ + int irevid; /* Revision ID..., long in Macintosh format */ + int itelescope; /* number of the CT which took the data */ + int irunnum; /* run number from the file name */ + int irunnumx; /* number read from the file, long in Macintosh format */ + int inumevents; /* number of events recorded in this run, long in + * Macintosh format */ + char cname[iMAXFILENAMELENGTH]; /* copy of the filename used when + * the data was read */ + int ieventheaderlength; /* number of bytes in the event header, + * long in Macintosh format */ + int ieventheaderoffset; /* unknown meaning, long in Macintosh + * format */ + int isignalslength; /* number of bytes used by the pixel signals, + * long in Macintosh format */ + int isignalsoffset; /* unknown meaning, long in Macintosh format */ + int itdccountslength; /* number of bytes used by the TDC signals, + * long in Macintosh format */ + int itdccountsoffset; /* unknown meaning, long in Macintosh format */ + int iobservtime_mins; /* (planned) observational time (minutes), + * short in Macintosh format */ + enum onoroff eruntype; /* indicates if the run is on- or off-source; + * in the Macintosh format, this is a short + * where 1 means off and 2 means on */ + int ileapsecs; /* Leap seconds, introduced into UTC in order to + * keep track of some very small differences + * accumulated through the years. */ + int ihv1; + int ihv2; /* these two numbers represent the HV adjustment made + * by the shiftpersonnel. The last three items are + * only read for Roberto files, but not for Tijarafe + * files, as e.g. the HVs are not written out correctly. */ + char cdummy[40]; /* unknown meaning */ + int imactime_secs[iMAXNUMRUNEVENTS]; + /* Macintosh time (time_t mydatum) in seconds since 1970, this is + * put together from two Macintosh shorts and converted to the time + * format required by SYSTEM */ + unsigned int uigpstime_secs[iMAXNUMRUNEVENTS]; + /* GPS time (time_t mydatum) in seconds since 1970, this is + * put together from two Macintosh shorts and converted to the time + * format required by SYSTEM + * This time is new from Roberto and is coming from the PC-GPS card. + * There is another flag available in the runheader, which gives an + * overall flag for the run, whether or not the GPS time had been + * working. It is, however, not a good idea to use this flag, as + * the workability of this card can also be dedued from the GPS time + * itself. This allows to check for failures on a event-by-event + * basis, which seems more reasonable than using the global flag. */ + unsigned int uigpstime_usecs[iMAXNUMRUNEVENTS]; + /* GPS microseconds -- available from the card, that's why it is + * put to the data, here it is. + * This is put together from two Macintosh shorts and converted + * to the time format required by SYSTEM */ + unsigned int uirubsecs_secs[iMAXNUMRUNEVENTS]; + /* the seconds of the Rubidium clock time this is put together + * from two Macintosh shorts and converted to the time format + * required by SYSTEM */ + unsigned int uirubsecfrac_200ns[iMAXNUMRUNEVENTS]; + /* the fractional part of the Rubidium clock time second, this is + * put together from two Macintosh shorts, unit is 200 ns */ + int isepos[2][iMAXNUMRUNEVENTS]; + /* positions of shaft encoders 1 and 2, short in Macintosh format */ + float fhourangle; /* This is the angle between the observation of this + * event and the culmination point. It is going to + * be written into the events NTuple. */ + int ieventnumber[iMAXNUMRUNEVENTS]; +/* number of the event, short in Macintosh format */ + int inumpixels; /* number of pixels in the camera, short in + * Macintosh format */ + int inummuonpixels; /* number of pixels in the muon shield; + * the sum of inumpixels and inummuonpixels is + * part of the event record but should be the + * same for all events, so it is not put into + * an array */ + int inumcointdcs; /* number of coincidence TDCs for which counts + * were recorded, short in Macintosh format + * this value is part of the event record but + * should be the same for all events, so it is + * not put into an array */ + int ipixsig_adc[iMAXNUMRUNEVENTS][iMAXNUMPIX]; + /* signals from the camera photo multipliers (ADC counts) */ + int imuonpixsig_adc[iMAXNUMRUNEVENTS][iMAXNUMMUONPIX]; + /* signals from the muon shield photo multipliers (ADC counts) */ + int itdcsig_cnts[iMAXNUMRUNEVENTS][iMAXNUMCOINTDCS]; + /* counts from the coincidence TDCs if itdccountslength > 0 */ +}; + +struct TELLOGDATA { + int irunnum; /* The runnumber which we are interested in. + * Only data for this runnumber are written into + * this structure and then copied to the calnt NTuple. + * + * This runnumber is generally the closest match + * inside the TelLog file. */ + float pixel_timeaveraged_current[127]; + float pixel_timeaveraged_scaler[127]; +}; + +/************************************************ + * structures in output file * + * (Konopelko files format) * + * (from structures.h file 4.4) * + ************************************************/ +struct mcbankheader { /* Once in every file */ + int iparticle_type; + /* Primary particle type: 0 = gamma; +-1= e+-; + * +-2 = mu+-; +-3 = p+-, n; +-4 = pi+- */ + int inum_cts; + /* Number of simulated telescopes */ + int inum_nfl_shifts; + /* Number of NFL shifts for hadrons (per event) */ + int itrigger_thresh_phot; + /* Only images with two pixels above this threshold are + * in the MC data file */ + int imin_signal; + /* "Tail cut" = smallest signal saved */ + int inum_sim_showers; + /* Total number of showers simulated for this file */ + float fatmothick_gcm2; + /* Atmosphere above the telescope [g/cm2] */ + float fcos_zenangle; + /* Cosine of the zenith angle */ + float fnfl_shift_radius_deg; + /* Maximum angular radius within which a shifted + * hadron shower is sampled [degrees] */ + float fenergybin_bound_tev[iNUMENERGYBINS+1]; + /* Boundaries for the 14 energy bins [TeV] + * (newline after the first 8 numbers) */ + float fzenithanglebin_bound[iNUMZENANGLEBINS+1]; + /* Boundaries for the 11 zenith angle bins [degrees] */ + int inum_show_in_bin[iNUMENERGYBINS]; + /* Number of simulated (or saved) showers + * per energy bin (newline after the first 8 numbers) */ + float fmaxdist_impact_m[iNUMIMPACTBINS]; + /* Maximum distance of the impact point to the + * central telescope [m] for each energy bin + * (newline after the first 8 numbers) */ + int inumpixels_for_ct[iMAXNUMCTS]; + /* Number of pixels in the camera of each simulated CT */ + float fpixwidth_deg_ct[iMAXNUMCTS]; + /* Pixel width [degrees] for each CT */ + float fmirrsize_m2_ct[iMAXNUMCTS]; + /* Mirror area [m^2] for each CT */ + float fmean_nsb_phot_ct[iMAXNUMCTS]; + /* Mean signal caused by the NSB in each pixel for each CT + * [photoelectrons]. This is the simulation of the NSB, + * not the electronic noise */ + float fphes_per_photon_ct[iMAXNUMCTS]; + /* Conversion factor photoelectron per photon */ + float frelative_x_ct[iMAXNUMCTS]; + /* x position relative to the central CT for each CT */ + float frelative_y_ct[iMAXNUMCTS]; + /* y position relative to the central CT for each CT */ +}; +#endif + + + + + + + Index: trunk/MagicSoft/Mars/mfileio/FileIOLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mfileio/FileIOLinkDef.h (revision 4459) +++ trunk/MagicSoft/Mars/mfileio/FileIOLinkDef.h (revision 4460) @@ -18,6 +18,3 @@ #pragma link C++ class MWriteRootFile+; -#pragma link C++ class MCT1ReadAscii+; -#pragma link C++ class MCT1ReadPreProc+; - #endif Index: trunk/MagicSoft/Mars/mfileio/MCT1ReadAscii.cc =================================================================== --- trunk/MagicSoft/Mars/mfileio/MCT1ReadAscii.cc (revision 4459) +++ (revision ) @@ -1,297 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 -! -! Copyright: MAGIC Software Development, 2000-2003 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// // -// MCT1ReadAscii // -// // -// Reads a ascii file with CT1 data. The file description and some example // -// files can be found on the Magic homepage. // -// // -// Input Containers: // -// -/- // -// // -// Output Containers: // -// MCerPhotEvt // -// // -///////////////////////////////////////////////////////////////////////////// -#include "MCT1ReadAscii.h" - -#include -#include - -#include -#include - -#include "MLog.h" -#include "MLogManip.h" - -#include "MParList.h" -#include "MCerPhotEvt.h" - -#include "MPedPhotPix.h" -#include "MPedPhotCam.h" - -ClassImp(MCT1ReadAscii); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// 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. -// -MCT1ReadAscii::MCT1ReadAscii(const char *fname, - const char *name, - const char *title) - : fIn(NULL) -{ - fName = name ? name : "MCT1ReadAscii"; - 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. -// -MCT1ReadAscii::~MCT1ReadAscii() -{ - delete fFileNames; - if (fIn) - delete fIn; -} - -// -------------------------------------------------------------------------- -// -// Add this file as the last entry in the chain -// -Int_t MCT1ReadAscii::AddFile(const char *txt, Int_t) -{ - TNamed *name = new TNamed(txt, ""); - fFileNames->AddLast(name); - return 1; -} - -// -------------------------------------------------------------------------- -// -// This opens the next file in the list and deletes its name from the list. -// -Bool_t MCT1ReadAscii::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(); - - const char *expname = gSystem->ExpandPathName(name); - fIn = new ifstream(expname); - - const Bool_t noexist = !(*fIn); - if (noexist) - { - *fLog << err << "Cannot open file " << expname << ": "; - *fLog << strerror(errno) << endl; - } - else - *fLog << inf << "Open file: '" << name << "'" << endl; - - delete [] expname; - - // - // Remove this file from the list of pending files - // - fFileNames->Remove(file); - - return !noexist; -} - -// -------------------------------------------------------------------------- -// -// 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 MPedPhotCam container to 127 pixels (CT1 camera) -// -Int_t MCT1ReadAscii::PreProcess(MParList *pList) -{ - // - // Preprocessing - // - - // - // Try to open at least one (the first) file - // - if (!OpenNextFile()) - return kFALSE; - - // - // 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 = (MPedPhotCam*)pList->FindCreateObj("MPedPhotCam"); - if (!fPedest) - return kFALSE; - - fPedest->InitSize(127); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Read apedestal entry (line) from the file -// -void MCT1ReadAscii::ReadPedestals() -{ - // - // skip the next 4 values - // - Float_t val; - - *fIn >> val; - *fIn >> val; - *fIn >> val; - *fIn >> val; - - // - // read in the next 127 numbers as the pedestals - // - for (Int_t i = 0; i<127; i++) - { - *fIn >> val; - - if (val>0) - (*fPedest)[i].Set(0, val); - } -} - -// -------------------------------------------------------------------------- -// -// Read a data entry (line) from the file -// -void MCT1ReadAscii::ReadData() -{ - // - // five unsused numbers - // - Int_t val; - - *fIn >> val; // ener - *fIn >> val; // zenang - *fIn >> val; // sec1 - *fIn >> val; // sec2 - - // - // read in the number of cerenkov photons and add the 'new' pixel - // too the list with it's id, number of photons and error - // - for (Int_t i = 0; i<127; i++ ) - { - Float_t nphot; - - *fIn >> nphot; - - if (nphot > 0.0) - fNphot->AddPixel(i, nphot, (*fPedest)[i].GetRms()); - } - fNphot->FixSize(); -} - -// -------------------------------------------------------------------------- -// -// 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. -// -Int_t MCT1ReadAscii::Process() -{ - // - // FIXME. This function should switch between reading pedestals and - // reading event data by the 'switch entry'. - // After reading it should set the InputStreamID correctly. - // (should use MPedPhotCam ) - // - - // - // read in the event nr - // - Int_t evtnr; - *fIn >> evtnr; - - // - // check if we are done. Try to open the next file in chain. - // If it was possible start reading. If not break the event loop - // - if (fIn->eof()) - return OpenNextFile() ? kCONTINUE : kFALSE; - - // - // if the first number is negativ this is a pedestal line: - // read in pedestals - // - // FIXME! Set InputStreamID - - if (evtnr < 0) - { - ReadPedestals(); - return kCONTINUE; - } - - ReadData(); - - return kTRUE; -} Index: trunk/MagicSoft/Mars/mfileio/MCT1ReadAscii.h =================================================================== --- trunk/MagicSoft/Mars/mfileio/MCT1ReadAscii.h (revision 4459) +++ (revision ) @@ -1,41 +1,0 @@ -#ifndef MARS_MCT1ReadAscii -#define MARS_MCT1ReadAscii - -#ifndef MARS_MTask -#include "MTask.h" -#endif - -class TList; -class MCerPhotEvt; -class MPedPhotCam; - -class MCT1ReadAscii : public MTask -{ -private: - ifstream *fIn; // the inputfile - MCerPhotEvt *fNphot; // the data container for all data. - MPedPhotCam *fPedest; // CT1 pedestals - TList *fFileNames; // Array which stores the \0-terminated filenames - - Bool_t OpenNextFile(); - - void ReadPedestals(); - void ReadData(); - - Int_t PreProcess(MParList *pList); - Int_t Process(); - -public: - MCT1ReadAscii(const char *filename=NULL, - const char *name=NULL, - const char *title=NULL); - - ~MCT1ReadAscii(); - - Int_t AddFile(const char *fname, Int_t dummy=-1); - - ClassDef(MCT1ReadAscii, 0) // Reads the CT1 data file -}; - -#endif - Index: trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.cc =================================================================== --- trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.cc (revision 4459) +++ (revision ) @@ -1,1171 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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, 11/2002 -! -! Copyright: MAGIC Software Development, 2000-2004 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// -// MCT1ReadPreProc -// -// Reads a output file of the CT1 preproc. -// -// Implements usage of a selector (see MRead) Use such a filter to skip -// events before reading! But never use a filter which needs read data -// as input... -// -// Input Containers: -// -/- -// -// Output Containers: -// MCerPhotEvt the data container for all data. -// MPedPhotCam CT1 pedestals -// MMcEvt monte carlo data container for MC files -// MMcTrig mc data container for trigger information -// MSrcPosCam source position in the camera -// MBadPixelsCam Array holding blind pixels -// -///////////////////////////////////////////////////////////////////////////// -#include "MCT1ReadPreProc.h" - -#include -#include - -#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 "MTime.h" -#include "MFilter.h" - -#include "MParList.h" -#include "MCerPhotEvt.h" - -#include "MPedPhotPix.h" -#include "MPedPhotCam.h" - -#include "MGeomCam.h" -#include "MSrcPosCam.h" -#include "MBadPixelsCam.h" -#include "MBadPixelsPix.h" - -#include "MRawRunHeader.h" -#include "MTaskList.h" - -#include "MMcEvt.hxx" -#include "MMcTrig.hxx" -#include "MBinning.h" - -#include "MParameters.h" - -ClassImp(MCT1ReadPreProc); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// 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 : "Reads a CT1 preproc data 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 -// -Int_t MCT1ReadPreProc::AddFile(const char *txt, int) -{ - const char *name = gSystem->ExpandPathName(txt); - - TString fname(name); - delete [] name; - - if (!CheckHeader(fname)) - { - *fLog << warn << "WARNING - Problem reading header... ignored." << endl; - return 0; - } - - const Int_t n = GetNumEvents(fname); - if (n==0) - { - *fLog << warn << "WARNING - File contains no data... ignored." << endl; - return 0; - } - - fEntries += n; - - *fLog << inf << "File " << txt << " contains " << n << " events (Total=" << fEntries << ")" << endl; - - fFileNames->AddLast(new TNamed(txt, "")); - return 1; -} - -// -------------------------------------------------------------------------- -// -// 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) - fRawRunHeader->SetRunNumber(outpars.irunnum); - fRawRunHeader->SetReadyToSave(); - - // 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; iSetReadyToSave(); - - /* 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; - fRawRunHeader->SetRunTime(outpars.istart_mjdate_day, outpars.iend_mjdate_day); - fRawRunHeader->SetReadyToSave(); - - /* - --- 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->Init(*fGeom); - - fPedRMS.Set(iMAXNUMPIX); - - *fLog << "PedestalRMS : "; - for (Int_t i=0; iSetReadyToSave(); - - // Used to communicate the mean over all pixels - // pedestal RMS into the Runs NTuple, as it might - // be used for e.g. quality cuts. - // float fpedrms_mean; - *fLog << "Pedestal RMS: " << outpars.fpedrms_mean << endl; - - // The average current over the active pixels - // for this run. This is done separately for - // ON and OF runs. - //float fcurrent_mean; - - // enum eERRORTOLERANCE {CAUTIOUS=0, GOODPHYSICS, TANK} eerrortolerance; - /* 0 == "cautious", exits on any reason (but tells in - * the .err file, - * 1 == "goodphysics", exits when physics could be affected - * by the error, - * 2 == "tank", never exits except on coredumps and when - * all files have been processed. Do not use such files for - * physics analysis! - * - * NOTE: the capital letter words are the enums, the small letter - * words must be used inside the parameter file. */ - - // enum eMCTRIGGERFLAG {ALL=0, FLAG, NOFLAG} emctriggerflag; - /* all: all events which survive the filter are written to the - * events NTuple. - * flag: When Dorota's triggerflag is set to 1 for a particular - * event, it shall be written to the output. All others shall - * just be disregarded. (Default) - * noflag: Opposite of 'flag': only events with triggerflag = 0 shall - * be treated further. */ - - *fLog << "Particle Id #" << outpars.imcparticle << endl; - *fLog << "Right Ascension: " << outpars.dsourcera_hours << "h" << endl; - *fLog << "Declination: " << outpars.dsourcedec_deg << "deg" << endl; - - // Next statement commented out because bmontecarlo was set wrongly - //fIsMcFile = outpars.bmontecarlo==TRUE; - fIsMcFile = (outpars.dsourcera_hours==0 && outpars.dsourcedec_deg==0 && - outpars.imcparticle != 0); - - if (fIsMcFile != (outpars.bmontecarlo==TRUE)) - { - *fLog << "File tells you that it is a "; - *fLog << (outpars.bmontecarlo ? "Monte Carlo" : "Real Data"); - *fLog << " file." << endl; - } - - *fLog << "File detected as a "; - *fLog << (fIsMcFile ? "Monte Carlo" : "Real Data"); - *fLog << " file." << endl; - *fLog << " " << endl; -} - -// -------------------------------------------------------------------------- -// -// Read CT1 PreProc File Header: -// -Int_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 (fpreprocversion<0.6) - { - *fLog << err << "Sorry, only files from PreProc V0.6 and newer are supported." << endl; - return kFALSE; - } - - // - // This is a stupid way of getting rid of numerical uncertanties when - // comparing floating point numbers (Argh...) - // - TString s1 = Form("%.2f", structversion); - TString s2 = Form("%.2f", STRUCT_VERSION); - - if (s1 != s2) - { - *fLog << warn << "WARNING: Version of C-structures of file (V"; - *fLog << s1 << ") not identical with current structures (V"; - *fLog << s2 << ")" << endl; - } - - fIn->read((char*)&outpars, sizeof(struct outputpars)); - - ProcessRunHeader(outpars); - - //rwagner: ReInit whenever new run commences - // rc==-1 means: ReInit didn't work out - - MTaskList *tlist = (MTaskList*)fParList->FindCreateObj("MTaskList"); - if (!tlist) - return -1; - - if (!tlist->ReInit(fParList)) - return -1; - - return kTRUE; -} - -Int_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 0; - } - - *fLog << inf << cheadertitle << flush; - - struct filterresults filterres; - fIn->read((char*)&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 - */ - - if (fNumEventsInRun!=(UInt_t)filterres.ifilter_passed_evts) - { - *fLog << err << "ERROR! Number of events in run (" << (UInt_t)filterres.ifilter_passed_evts; - *fLog << ") doesn't match number of read events ("; - *fLog << fNumEventsInRun << ")" << endl; - *fLog << " File corrupted." << endl; - return -1; - } - - fNumFilterEvts += fNumEventsInRun; - fNumRuns++; - - *fLog << inf << "Read " << fNumEventsInRun << " events from run (Total="; - *fLog << fNumFilterEvts << "/" << fEntries << " ["; - *fLog << 100*fNumFilterEvts/fEntries << "%], Runs=" << fNumRuns << ")"; - *fLog << endl; - - return 1; -} - -// -------------------------------------------------------------------------- -// -// 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 existence of a next file to read - // - if (fNumFile >= (UInt_t)fFileNames->GetSize()) - return kFALSE; - - TNamed *file = (TNamed*)fFileNames->At(fNumFile); - - //TNamed *file = (TNamed*)fFileNames->GetFirst(); - //if (!file) - // return kFALSE; - - // - // open the file which is the first one in the chain - // - const TString name = file->GetName(); - - const char *expname = gSystem->ExpandPathName(name); - const TString fname(expname); - delete [] expname; - - // - // Remove this file from the list of pending files - // - //fFileNames->Remove(file); - - *fLog << inf << "Open file: '" << name << "'" << endl; - - if (!CheckHeader(fname)) - return kFALSE; - - fNumFile++; - - fIn = new ifstream(fname); - if (!*fIn) - { - *fLog << err << "Cannot open file " << fname << ": "; - *fLog << strerror(errno) << endl; - return kFALSE; - } - - *fLog << inf << "-----------------------------------------------------------------------" << endl; - - switch (ReadRunHeader()) - { - case kFALSE: - *fLog << warn << "Unable to read first run header... skipping file." << endl; - return kFALSE; - case -1: - *fLog << warn << "ReInit of Tasklist didn't succeed." << endl; - return kFALSE; - default: - *fLog << "After opening next file: Number of Events #" << fNumEventsInRun << endl; - return kTRUE; - } -} - -Bool_t MCT1ReadPreProc::CheckHeader(const TString fname) const -{ - ifstream fin(fname); - if (!fin) - { - *fLog << dbginf << err << "ERROR - Cannot open file '" << fname << "'" << endl; - return kFALSE; - } - - char cheadertitle[iHEADERTITLELENGTH]; - fin.read(cheadertitle, iHEADERTITLELENGTH); - - Float_t fpreprocversion, structversion; - Int_t dummyi; - - sscanf(cheadertitle, cTITLE_TEMPLATE, - &fpreprocversion, &structversion, - &dummyi, &dummyi, &dummyi, &dummyi); - - if (fpreprocversion < 0.6) - { - *fLog << dbginf << err << "ERROR - You must use PreProc V0.6 or higher." << endl; - return kFALSE; - } - - if (STRUCT_VERSION > structversion) - { - *fLog << warn << "WARNING: Version of C-structures of file (V"; - *fLog << structversion << ") newer than current structures (V"; - *fLog << STRUCT_VERSION << ")" << endl; - } - - *fLog << "Current structures: " << STRUCT_VERSION << " "; - *fLog << "Structures in file: " << structversion << " "; - *fLog << "Used preproc version: " << fpreprocversion << endl; - - return kTRUE; -} - - -Int_t MCT1ReadPreProc::GetNumEvents(const TString fname) const -{ - *fLog << inf << "Scanning file " << fname << " for size" << flush; - - ifstream fin(fname); - if (!fin) - { - *fLog << dbginf << err << "ERROR - Opening file." << endl; - return 0; - } - - 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.peek()!=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; -} - -Bool_t MCT1ReadPreProc::Rewind() -{ - fNumFilterEvts = 0; - fNumEvents = 0; - fNumRuns = 0; - fNumFile = 0; - if (fIn) - delete fIn; - - fIn=NULL; - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// 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 MPedPhotCam container to 127 pixels (CT1 camera) -// -Int_t MCT1ReadPreProc::PreProcess(MParList *pList) -{ - fParList = pList; - - // - // look for the HourAngle container in the plist - // - fHourAngle = (MParameterD*)pList->FindCreateObj("MParameterD", "HourAngle"); - if (!fHourAngle) - return kFALSE; - fHourAngle->SetTitle("Store the CT1 hour angle [deg]"); - - // - // look for the ThetaOrig container in the plist - // - fThetaOrig = (MParameterD*)pList->FindCreateObj("MParameterD", "ThetaOrig"); - if (!fThetaOrig) - return kFALSE; - fThetaOrig->SetTitle("Store the original CT1 zenith angle [rad]"); - - // - // 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 = (MPedPhotCam*)pList->FindCreateObj("MPedPhotCam"); - if (!fPedest) - return kFALSE; - - // - // look for the time class in the plist - // - fTime = (MTime*)pList->FindCreateObj("MTime"); - if (!fTime) - return kFALSE; - - // - // look for the pedestal class in the plist - // - fBlinds = (MBadPixelsCam*)pList->FindCreateObj("MBadPixelsCam"); - if (!fBlinds) - return kFALSE; - - // - // look for the source position in the camera - // - fSrcPos = (MSrcPosCam*)pList->FindCreateObj("MSrcPosCam"); - 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; - - // - // look for the raw run header class - // - fRawRunHeader = (MRawRunHeader*)pList->FindCreateObj("MRawRunHeader"); - if (!fRawRunHeader) - return kFALSE; - - fBinningT = (MBinning*)pList->FindObject("BinningTheta"); - - Rewind(); - - fPedest->Init(*fGeom); - - return GetSelector() ? GetSelector()->CallPreProcess(pList) : kTRUE; -} - - -// -------------------------------------------------------------------------- -// -// Smear Theta uniformly in a bin of Theta -// theta [rad] -// SmearTheta [rad] -// -Double_t MCT1ReadPreProc::SmearTheta(Double_t theta) -{ - if (!fBinningT) - return theta; - - const Int_t bin = fBinningT->FindLoEdge(theta * 180/TMath::Pi()); - if (bin<0) - return theta; - - // smear Theta within the Theta bin - const Double_t low = fBinningT->GetEdges()[bin]; - const Double_t up = fBinningT->GetEdges()[bin+1]; - - // "up-": Do not allow the upper edge - return (up - gRandom->Uniform() * (up-low)) * TMath::Pi()/180; -} - -// -------------------------------------------------------------------------- -// -// Discretize Theta according to the binning in Theta -// theta [rad] -// DiscreteTheta [rad] (= bin center) -// -Double_t MCT1ReadPreProc::DiscreteTheta(Double_t theta) -{ - if (!fBinningT) - return theta; - - const Int_t bin = fBinningT->FindLoEdge(theta * 180/TMath::Pi()); - if (bin<0) - return theta; - - const Double_t low = fBinningT->GetEdges()[bin]; - const Double_t up = fBinningT->GetEdges()[bin+1]; - - return 0.5*(up+low) * TMath::Pi()/180; -} - -// -------------------------------------------------------------------------- -// -// Analyse the event data, means store and copy the needed data into -// Mars structures and data members -// -Bool_t MCT1ReadPreProc::ProcessEvent(const struct eventrecord &event) -{ - /* - if (fRawRunHeader->GetRunNumber() == 1) - { - *fLog << "eventrecord" << endl; - *fLog << "isecs_since_midday = " << event.isecs_since_midday << endl; - *fLog << "isecfrac_200ns = " << event.isecfrac_200ns << endl; - *fLog << "snot_ok_flags = " << event.snot_ok_flags << endl; - *fLog << "ialt_arcs = " << event.ialt_arcs << endl; - *fLog << "iaz_arcs = " << event.iaz_arcs << endl; - *fLog << "ipreproc_alt_arcs = " << event.ipreproc_alt_arcs << endl; - *fLog << "ipreproc_az_arcs = " << event.ipreproc_az_arcs << endl; - *fLog << "ifieldrot_arcs = " << event.ifieldrot_arcs << endl; - - *fLog << "srate_millihz = " << event.srate_millihz << endl; - *fLog << "fhourangle = " << event.fhourangle << endl; - *fLog << "fmcenergy_tev = " << event.fmcenergy_tev << endl; - *fLog << "fmcsize_phel = " << event.fmcsize_phel << endl; - *fLog << "imcimpact_m = " << event.imcimpact_m << endl; - *fLog << "imcparticle = " << event.imcparticle << endl; - *fLog << "imctriggerflag = " << event.imctriggerflag << endl; - } - */ - - // reset blind pixels for this event - fBlinds->Clear(); - for (int i=0; iSetCT1Time((UInt_t)fRawRunHeader->GetRunStart().GetMjd(), event.isecfrac_200ns, event.isecs_since_midday); - fTime->SetReadyToSave(); - - /* - --- USEFULL? NEEDED? --- - 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 - // 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].GetRms()); - } - fNphot->FixSize(); - fNphot->SetReadyToSave(); - - // int ipreproc_alt_arcs; // "should be" alt according to preproc (arcseconds) - // int ipreproc_az_arcs; // "should be" az according to preproc (arcseconds) - - // smear Theta in its Theta bin - const Double_t theta = TMath::Pi()*(0.5-1./180*event.ialt_arcs/3600); - fThetaOrig->SetVal(theta); - - // store hour angle - fHourAngle->SetVal(event.fhourangle); - - fMcEvt->Fill(event.isecs_since_midday, //0, /*fEvtNum*/ - fIsMcFile ? event.imcparticle : 0, /*corsika particle type*/ - fIsMcFile ? event.fmcenergy_tev*1000 : 0, - 0, /* fThi0 */ - 0, /* fFirTar */ - 0, /* fzFirInt */ - // 0, /* fThet*/ - // rwagner: The following should be theta, right? Check with - // altitude fill some lines down... - 0, // altitude (arcseconds) - 0, /* fPhii */ - 0, /* fCorD */ - 0, /* fCorX */ - 0, /* fCorY */ - fIsMcFile ? event.imcimpact_m*100 : 0, - TMath::Pi()/180*event.iaz_arcs/3600, // azimuth (arcseconds) - // fIsMcFile ? SmearTheta(theta) : theta, - DiscreteTheta(theta), - 0, /* fTFirst */ - 0, /* fTLast */ - 0, /* fL_Nmax */ - 0, /* fL_t0 */ - 0, /* fL_tmax */ - 0, /* fL_a */ - 0, /* fL_b */ - 0, /* fL_c */ - 0, /* fL_chi2 */ - 0, /* uiPin */ - 0, /* uiPat */ - 0, /* uiPre */ - 0, /* uiPco */ - 0, /* uiPelS */ - (int)(fIsMcFile ? event.fmcsize_phel : 0), /* uiPelC, Simulated SIZE */ - 0, /* elec */ - 0, /* muon */ - 0 /* other */ - ); - - fMcTrig->SetFirstLevel(event.imctriggerflag); // MC data from Dorota get a triggerflag: 1 means triggered, 0 not. */ - - fMcTrig->SetReadyToSave(); - fMcEvt->SetReadyToSave(); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// 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 this isn't successful... - // must be an event - // - switch (ReadRunFooter()) - { - case -1: - return kFALSE; - case 0: - return kTRUE; - } - - *fLog << inf << "Footer found." << endl; - - const char c = fIn->peek(); - - // - // No after reading the footer check if we reached the end of the file - // - if (fIn->eof() || c==EOF) - { - *fLog << "End of file." << endl; - return kFALSE; - } - - // - // If the eof isn't reached a new header must follow. Check for it. - // - if (c!=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() < 0) - { - *fLog << warn << "ReInit of Tasklist didn't succeed." << endl; - return kFALSE; - } - } -} - -// -------------------------------------------------------------------------- -// -// 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. -// -Int_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; - - // - // Check for a selector. If one is given and returns kFALSE - // skip this event. - // - if (GetSelector()) - { - // - // Make sure selector is processed - // - if (!GetSelector()->CallProcess()) - { - *fLog << err << dbginf << "Processing Selector failed." << endl; - return kFALSE; - } - - // - // Skip Event - // - if (!GetSelector()->IsConditionTrue()) - { - fIn->seekg(sizeof(struct eventrecord), ios::cur); - - fNumEvents++; - fNumEventsInRun++; - - return kCONTINUE; - } - } - - // event data to be read from the file - struct eventrecord event; - - - - // read the eventrecord from the file - fIn->read((char*)&event, sizeof(struct eventrecord)); - - switch (ProcessEvent(event)) - { - case kFALSE: - return kFALSE; - case kCONTINUE: - return kCONTINUE; - } - - fNumEvents++; - fNumEventsInRun++; - - return kTRUE; -} - -Int_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; - *fLog << "Number of events detected first: " << fEntries << 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; - } - - delete fIn; - fIn = NULL; - - return GetSelector() ? GetSelector()->CallPostProcess() : kTRUE; -} Index: trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.h =================================================================== --- trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.h (revision 4459) +++ (revision ) @@ -1,104 +1,0 @@ -#ifndef MARS_MCT1ReadPreProc -#define MARS_MCT1ReadPreProc - -#ifndef ROOT_TArrayF -#include -#endif - -#ifndef ROOT_TArrayI -#include -#endif - -#ifndef MARS_MRead -#include "MRead.h" -#endif - -class TList; -class MTime; -class MMcEvt; -class MMcTrig; -class MGeomCam; -class MSrcPosCam; -class MCerPhotEvt; -class MPedPhotCam; -class MBadPixelsCam; -class MRawRunHeader; -class MTaskList; -class MParList; -class MParameterD; -class MBinning; - -struct outputpars; -struct eventrecord; - -class MCT1ReadPreProc : public MRead -{ -private: - ifstream *fIn; // the inputfile - TList *fFileNames; // Array which stores the \0-terminated filenames - - MGeomCam *fGeom; // camera geometry - MCerPhotEvt *fNphot; // the data container for all data. - MPedPhotCam *fPedest; // ct1 pedestals - MTime *fTime; // event time - MMcEvt *fMcEvt; // monte carlo data container for MC files - MMcTrig *fMcTrig; // mc data container for trigger information - MSrcPosCam *fSrcPos; // source position in the camera - MBadPixelsCam *fBlinds; // Array holding blind pixels - MRawRunHeader *fRawRunHeader; // raw run header - MParList *fParList; // parameter list - MParameterD *fHourAngle; // hour angle [deg] - MParameterD *fThetaOrig; // original zenith angle [rad] - MBinning *fBinningT; // Theta binning for the smearing - - Bool_t fIsMcFile; // Flag whether current run is a MC run - - UInt_t fNumFile; - UInt_t fNumEvents; // number of events counted in all runs in all files - UInt_t fNumEventsInRun; // number of events in the counted in th ecurrent run - UInt_t fNumRuns; // number of processed runs of all files - UInt_t fEntries; // entries of all files succesfully added - UInt_t fNumFilterEvts; // number of events mentioned in the runs footers - - TArrayF fPedRMS; - TArrayI fBlnd; - - Bool_t OpenNextFile(); - - Int_t GetNumEvents(const TString name) const; - Bool_t CheckHeader(const TString fname) const; - - void ReadPedestals(); - Int_t ReadRunHeader(); - Int_t ReadRunFooter(); - Bool_t CheckFilePosition(); - void ProcessRunHeader(const struct outputpars &outpars); - Bool_t ProcessEvent(const struct eventrecord &event); - - Double_t SmearTheta(Double_t theta); - Double_t DiscreteTheta(Double_t theta); - - Int_t PreProcess(MParList *pList); - Int_t Process(); - Int_t PostProcess(); - - Bool_t Rewind(); - -public: - MCT1ReadPreProc(const char *filename=NULL, - const char *name=NULL, - const char *title=NULL); - - ~MCT1ReadPreProc(); - - Int_t AddFile(const char *fname, int i=0); - - UInt_t GetEntries() { return fEntries; } - - ClassDef(MCT1ReadPreProc, 0) // Reads the CT1 preproc data file -}; - -#endif - - - Index: trunk/MagicSoft/Mars/mfileio/Makefile =================================================================== --- trunk/MagicSoft/Mars/mfileio/Makefile (revision 4459) +++ trunk/MagicSoft/Mars/mfileio/Makefile (revision 4460) @@ -12,6 +12,5 @@ INCLUDES = -I. -I../mbase -I../mraw -I../mmc -I../mreflector -I../mgui \ - -I../mdata -I../manalysis -I../mgeom -I../mhbase \ - -I../mpointing -I../mpedestal + -I../mdata -I../mbadpixels #mraw: MRawRunHeader (MReadMaraFile) #mmc: MMcRunHeader (MReadMarsFile) @@ -19,9 +18,4 @@ #mgui: MCamEvent (MReadRflFile) #mdata: MData* (MWriteAsciiFile) -#manalysis: MCerPhotEvt (MCT1ReadAscii) -#mgeom: MGeomCam (MCT1ReadPreProc) -#mhbase: MBinning (MCT1ReadPreProc) -#mpointing: MPointingPos (MCT1ReadPreProc) -#mpedestal: MPedestalCam (MCT1ReadAscii) CINT = FileIO @@ -35,7 +29,5 @@ MWriteFile.cc \ MWriteAsciiFile.cc \ - MWriteRootFile.cc \ - MCT1ReadAscii.cc \ - MCT1ReadPreProc.cc + MWriteRootFile.cc ############################################################ Index: trunk/MagicSoft/Mars/mfileio/structures.h =================================================================== --- trunk/MagicSoft/Mars/mfileio/structures.h (revision 4459) +++ (revision ) @@ -1,882 +1,0 @@ -/******************************************************************/ -/* GGG AAA M M M M AAA SSS */ -/* G G A A MM MM MM MM A A S S */ -/* G A A M M M M M M A A SS */ -/* G GG AAAAA M M M M AAAAA SS */ -/* G G A A M M M M A A S S */ -/* GGG A A M M M M A A SSS */ -/* Gamma Astronomer's Munich Madrid Analysis Suite */ -/* An Atmospheric Cherenkov Telescope Data Analysis Software */ -/* MPI f"ur Physik, M"unchen & Universidad Complutense, Madrid */ -/******************************************************************/ - -/******************************************************************/ -/* This File belongs to the programs */ -/* */ -/* P R E P R O C E S S O R and I M A G E R */ -/* */ -/* Purpose: provide the structure definitions common to both */ -/* programs, especially the output structures of the */ -/* preprocessor */ -/* */ -/******************************************************************/ - -/******************************************************************/ -/* @(#) File name structures.h - @(#) latest SID 0.5 - @(#) Date of latest delta 26/6/00 - Author: D. Petry - changed by M. Kestel */ -/******************************************************************/ - -#ifndef _structures_h_ -#define _structures_h_ "__structures_h__" - -#define STRUCT_VERSION 0.6 - -#define iHEADERTITLELENGTH 60 -/* length of the string written at the beginning of the */ -/* output for each runfile */ -#define cTITLE_TEMPLATE "PREPROC V%f/S%f CT %d RUN %d %d PROCMJD %d\n" - -#define cEND_EVENTS_TEMPLATE "EVENTS AFTER FILTER = %d\n" - -#define iMAXNUMPIX 127 -/* maximum value for the number of camera pixels (excl. muon shield) */ -#define iMAXNUMMUONPIX 5 -/* maximum value for the number of pixels in the muon shield */ -#define iMAXNUMCOINTDCS 5 -/* maximum value for the number of coincidence tdcs */ - - -#define iBADMACTIME_BIT 1 -/* snot_ok_flags is ORed with this number to set the bit which - * indicates that the mactime measurement didn't work */ -#define iBADRUBTIME_BIT 2 -/* snot_ok_flags is ORed with this number to set the bit which - * indicates that the Rubidium time measurement didn't work */ -#define iBADPOINTING_BIT 4 -/* snot_ok_flags is ORed with this number to set the bit which - * indicates that the shaft encoder values were inconsistent with - * time OR source coordinates */ -#define iMUON_BIT 8 - /* snot_ok_flags is ORed with this number to set the bit which - * indicates that the event didn't pass the muon filter - * (for test puposes) */ - -#define i2OUTFAIL_BIT 16 -/* snot_ok_flags is ORed with this number to set the bit which - * indicates that the event didn't pass the "two out of all" - * software trigger */ -#define iSUMFAIL_BIT 32 -/* snot_ok_flags is ORed with this number to set the bit which - * indicates that the event didn't pass the "sum of all - * photoelectrons" filter */ - -#define iPARTICLETYPE_FACTOR 1024 -/* for MC data: the particle type multiplied with this number - * is ORed with snot_ok_flags */ - -enum onoroff {NEITHER_ON_NOR_OFF, OFF_SOURCE, ON_SOURCE}; -enum trackmodes {NORMAL, REVERSE, DUNNO}; -enum axes {RA, DEC, ALT, AZ}; -enum eERRORTOLERANCE {CAUTIOUS=0, GOODPHYSICS, TANK}; -enum eMCTRIGGERFLAG {ALL=0, FLAG, NOFLAG}; - -struct outputpars { - /* preprocessing parameters that are passed to imager. - * Note: Not all parameters from the parameter file are passed to - * imager, only those that are repeated in this struct outputpars - * here !! That's why these two structs must be the same for both - * preproc and imager. */ - - int itelescope; - /* number of the CT which took the data (from parameters file */ - /* crosschecked with run file name) */ - float flongitude_deg; - /* longitude (counted positive towards West) of CT position */ - float flatitude_deg; - /* latitude (counted positive towards North) of CT position */ - int irunnum; - /* run number (from parameters file cross checked with run */ - /* file name) */ - enum onoroff eruntype; -#ifdef LINUX - int dummy; -#endif - /* indicates if the */ - /* run is on- or off source (from runfile cross checked with */ - /* file name) */ - double dsourcera_hours; - /* right ascension of observed source in hours */ - double dsourcedec_deg; - /* declination of observed source in degrees */ - int inumpixels; - /* number of pixels in the camera (from parameters file cross */ - /* checked with run file, calibration file and pedestal file) */ - int inummuonpixels; - /* number of pixels in the muon shield (from parameters file */ - /* cross checked with run file, calibration file and pedestal */ - /* file) */ - int inumcointdcs; - /* number of coincidence tdcs recorded in the runfile (from */ - /* parameters file cross checked */ - /* with run file) */ - float fpixdiameter_deg; - /* smallest pixel diameter (degrees) (from parameters file) */ - enum axes ese1_is; - /* name of the axis to which shaft encoder 1 is attached */ - /* (implies the type of mount) */ - int isezeropos[2]; - /* zero position of shaftencoders 1 and 2 (from parameters file) */ - 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 */ - float fbendingcorr; - /* bending correction factor (ALT/AZ mount only) */ - float fextinction; - /* atmospheric extinction (typically taken from the Carlsberg */ - /* Meridian Circle data) */ - Boolean bdontusepix[iMAXNUMPIX]; - /* bdontusepix is set true if the pixel should due to whatever - * reason not be used in image analysis, otherwise it is true; - * - * bdontusepix is set by the user inside preproc.par and is - * untouched by preproc (version >= 0.4). */ - Boolean bexcludepix[iMAXNUMPIX]; - /* 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 */ - int irubminusutc_usecs; - /* difference between rubidium clock and UTC in microseconds */ - int isum_thresh_phot; - /* threshold for the total sum of photoelectrons filter */ - /* from the parameters file */ - int i2out_thresh_phot; - /* threshold for the two-pixels-out-of-all software */ - /* trigger from parameters file */ - int imuoncut_thresh_adccnt[iMAXNUMMUONPIX]; - /* thresholds for the muon cut */ - Boolean bmuon_suppression; - /* "Dolby" noise reduction flag */ - float ftolerated_pointerror_deg; - /* maximum tolerated pointing error in the position - * check in iprocessdata */ - float fxpointcorr_deg; - /* pointing correction (to be added along the camera - * x axis) e.g. from point run */ - float fypointcorr_deg; - /* pointing correction (to be added along the camera - * y axis) e.g. from point run */ - 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. - * This is given once per run. */ - - /* 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) */ - int iproc_mjdate; - /* MJD of data processing (i.e. creation of this file) */ - enum trackmodes etrackmode; - /* tracking mode (for ALT/AZ CTs) */ - int iproc_evts; -#ifdef LINUX - int dummy2; -#endif - /* number of events processed */ - 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 */ - float fpedrms_mean; /* Used to communicate the mean over all pixels - * pedestal RMS into the Runs NTuple, as it might - * be used for e.g. quality cuts. */ - float fcurrent_mean; /* The average current over the active pixels - * for this run. This is done separately for - * ON and OF runs. */ - - enum eERRORTOLERANCE eerrortolerance; - /* 0 == "cautious", exits on any reason (but tells in - * the .err file, - * 1 == "goodphysics", exits when physics could be affected - * by the error, - * 2 == "tank", never exits except on coredumps and when - * all files have been processed. Do not use such files for - * physics analysis! - * - * NOTE: the capital letter words are the enums, the small letter - * words must be used inside the parameter file. */ - enum eMCTRIGGERFLAG emctriggerflag; - /* all: all events which survive the filter are written to the - * events NTuple. - * flag: When Dorota's triggerflag is set to 1 for a particular - * event, it shall be written to the output. All others shall - * just be disregarded. (Default) - * noflag: Opposite of 'flag': only events with triggerflag = 0 shall - * be treated further. */ - Boolean bmontecarlo; - /* if TRUE we have a monte carlo dataset before us, if FALSE - * (default value) it is just normal data. */ -}; - -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 */ - -}; - -struct eventrecord { - 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 */ - int ialt_arcs; - /* altitude (arcseconds) */ - int iaz_arcs; - /* azimuth (arcseconds) */ - int ipreproc_alt_arcs; - /* "should be" alt according to preproc (arcseconds) */ - int ipreproc_az_arcs; - /* "should be" az according to preproc (arcseconds) */ - int ifieldrot_arcs; - /* 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 */ - unsigned short srate_millihz; - /* 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 */ - float fhourangle; /* 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 fmcenergy_tev; /* Simulated Energy.... dropping through to - * the Events NTuple. */ - float fmcsize_phel; /* Simulated SIZE.... dropping through to - * the Events NTuple. */ - int imcimpact_m; - /* MC data contain the impact parameter, which is given here in - * meters. */ - int imcparticle; - /* MC data know which particle they are.... all in CORSIKA standard. */ - int imctriggerflag; - /* MC data from Dorota get a triggerflag: 1 means triggered, 0 not. */ - short spixsig_10thphot[iMAXNUMPIX]; - /* number of photoelectrons measured in each pixel only the - * actual number of pixels (outputpars.inumpixels) is written out */ -}; - -struct camera { /* camera parameters for imaging */ - int inumpixels; - int inumrings; - double dpixdiameter_deg; - double dxc[iMAXNUMPIX]; -/* Pixel coordinates in camera coordinate system (x points from - * pixel 1 to 2). */ - double dyc[iMAXNUMPIX]; - /* The numbering of the pixels in these arrays starts at 0! */ - double dxpointcorr_deg; - /* correction of the pixel coordinates; to be added to dxc[] - * to get correct value */ - double dypointcorr_deg; - /* correction of the pixel coordinates; to be added to dxc[] - * to get correct value */ -}; - - -/* two structures for better file handling */ -struct inputfile { - char cname[iMAXFILENAMELENGTH]; - /* filename (including path) */ - FILE *pointer; - /* filepointer */ - char ccannotopentext[161]; - /* Error text printed when file cannot be opened */ - int icannotopencode; - /* Error code for the exit statement when file could not be opened */ - int iearlyeofcode; - /* Error code for the exit statement for unexpected EOF */ - int ierroratclosecode; - /* Error code for the exit statement for error while trying to close */ -}; - -struct outputfile { - char cname[iMAXFILENAMELENGTH]; /* filename (including path) */ - FILE *pointer; /* filepointer */ - char ccannotopentext[161]; - /* Error text printed when file cannot be opened */ - int icannotopencode; - /* Error code for the exit statement when file could not be opened */ - int icannotwritecode; - /* Error code for the exit statement for failed fprintf */ - int ierroratclosecode; - /* Error code for the exit statement for error while trying to close */ -}; - -struct P2ALLHISTOS { - - HISTO *p2adchist; - HISTO *p2pedonhist; - HISTO *p2pedofhist; - HISTO *p2pedonmeanhist; - HISTO *p2pedofmeanhist; - HISTO *p2pedonjitthist; - HISTO *p2pedofjitthist; - HISTO *p2calonhist; - HISTO *p2calofhist; - HISTO *p2calonmeanhist; - HISTO *p2calofmeanhist; - HISTO *p2calonjitthist; - HISTO *p2calofjitthist; - HISTO *p2phehist; - HISTO *p2zenonhist; - HISTO *p2zenofhist; - HISTO *p2corrpedhist; - HISTO *p2pedfwhmhist; -}; - -struct CALNT_DAT { - /* the structure for an entry in the calibration ntuple */ - float CTrunnum; - float runtype; - float ALT_deg; /* the ALT taken from the first event of the run */ - float extinct; /* extinction */ - float RtiSecBF; /* runtime in seconds before filter */ - float pixnum; /* number of the pixel for which the following - * data is given */ - float pixonoff; /* Quantity telling whether this pixel was switched - * off by the user (-= 2.) or by preproc (-= 1.). - * If ON, its value = 0. */ - float ped_adcc; /* the standard deviation of the pedestal as - * calculated by TIJARAFE */ - float pedrmsac; /* the standard deviation of the pedestal in - * ADC counts as calculated by TIJARAFE */ - float pedrmsPE; /* the calibrated ped RMS in photoelectrons - * which is also given to the imager */ - float las_adcc; /* the mean signal for this pixel from the - * laser run */ - float lasrmsac; /* RMS of the laser events after cleaning from cosmics */ - float conv_fac; /* conversion factor = /Var(Q) * F^2 */ - float F; /* The F from the line before */ - float rel_gain; /* the relative gain of the pixel */ - float numtrig; /* number of events in which the pixel was - * above trigger threshold */ - float num3sig; /* number of events in which the pixel was - * 3 sigma above its pedestal */ - float HV1; /* Adjustment of HV1 for this run, if available. */ - float HV2; /* Adjustment of HV2 for this run, if available. */ - float ThrCurr; - float AvgCurr; /* the sum of all currents during this run (with the - * particular runtype in question) from the TelLog - * file. The currents are averaged over the 2-minute - * intervals but summed over all pixels with currents - * higher than the threshold current value. This sum - * is then stored in the runs ntuple. - * Same for the scaler values --> */ - float AvgScal; -}; - -struct inputpars { - Boolean bkeywordgiven[iNUMKEYWORDS+1]; /* if a valid keyword - * is given in the - * parameters file - * the corresponding - * element in this - * array is set to true */ - int itelescope; /* number of the CT which took the data */ - float flongitude_deg; /* longitude (counted positive - * towards West) of CT position */ - float flatitude_deg; /* latitude (counted positive - * towards North) of CT position */ - int ifirstrunnum; /* first run number to which these - * parameters apply*/ - int ilastrunnum; /* last run number to which these - * parameters apply, */ - /* i.e. this run will be processed - * with these parameters */ - int inumpixels; /* number of pixels in the camera */ - 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) */ - enum axes ese1_is; /* name of the axis to which shaft - * encoder 1 is attached - * (implies the type of mount) (the - * type axes is declared in structures.h) */ - float fdegrees_per_step[2]; /* angular resolution of shaft - * encoders 1 and 2 */ - int isezeropos[2]; /* zero position of shaftencoders 1 and - * 2 from parameters file */ - 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 */ - float fbendingcorr; /* bending correction factor - * (ALT/AZ mount only) */ - Boolean bdontusepix[iMAXNUMPIX]; /* bdontusepix is set true - * if the pixel should due - * to whatever reason not be - * used in image analysis, - * otherwise it is false; - * this is a copy of the - * input from the parameters file */ - Boolean bdontusepix_in_trig[iMAXNUMPIX]; /* is set true if the - * pixel should due - * to whatever reason not - * be used in the two out - * of all trigger, otherwise - * it is false; this is a - * copy of the input from - * the parameters file */ - float fphotoel_per_adccnt[iMAXNUMPIX]; /* conversion factors for - * the pixel signals */ - float fextinction; /* atmospheric extinction (typically - * taken from the Carlsberg Meridian - * Circle data) */ - int irubminusutc_usecs; /* difference between rubidium clock - * and UTC in microseconds */ - int isum_thresh_phot; /* threshold for the total sum of - * photoelectrons filter from the - * parameters file */ - int i2out_thresh_phot; /* threshold for the - * two-pixels-out-of-all software - * trigger from parameters file */ - int imuoncut_thresh_adccnt[iMAXNUMMUONPIX]; /* thresholds for - * the muon cut */ - Boolean bmuon_suppression; /* if true, the events which do - * not pass the muon cut are not - * written to the output */ - float ftolerated_pointerror_deg; /* maximum tolerated pointing - * error in the position check - * in iprocessdata */ - float fxpointcorr_deg; /* pointing correction (to be added - * along the camera x axis) e.g. - * from point run */ - float fypointcorr_deg; /* pointing correction (to be added - * along the camera y axis) e.g. - * from point run */ - 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) */ - enum pedsearchdirs {MATCH=0, BEFORE, AFTER, PINGPONG, - PEAKFWHM, NONE} epedsearchdir ; - /* MATCH = only same number as run file (if ipedendofsearch is 0) or - exactly ipedendofsearch; - BEFORE = search backwards until pedendofsearch; - AFTER = search forward until pedendofsearch, - PINGPONG = Try to fibnd the matching partner, try run#-1, run#+1, - run#-2, run#+2 etc. up to the third argument, the partnerfindrange - NONE = use none of both. This is only useful for Monte Carlo - Analysis, for everything else it is highly unrecommended.*/ - int ipedendofsearch; /* pedestal file number until which to - * search (see epedsearchdir) */ - int ipedsearchdepth; /* depth, until which preproc should - * try to find TelPed* partnerfiles */ - enum pedsearchdirs ecalsearchdir; - /* MATCH = only same number as run file (if ipedendofsearch is 0) or - exactly ipedendofsearch; - BEFORE = search backwards until pedendofsearch; - AFTER = search forward until pedendofsearch, - PINGPONG = Try to fibnd the matching partner, try run#-1, run#+1, - run#-2, run#+2 etc. up to the third argument, the partnerfindrange - NONE = use none of both. This is only useful for Monte Carlo - Analysis, for everything else it is highly unrecommended.*/ - int icalendofsearch; /* calibration file number until which - * to search (see ecalsearchdir)*/ - int icalsearchdepth; /* depth, until which preproc should - * try to find TelCal* partnerfiles */ - double dsourcera_hours; /* right ascension of observed source - * in hours */ - double dsourcedec_deg; /* declination of observed source in - * degrees */ - float fpadlevel_phot_pix[iMAXNUMPIX]; /* if > 0., software - * padding to this - * level is applied - * for the given pixel */ - float fseed1; /* seed for random generators */ - float fseed2; - - int iswap_partner_of_pix[iMAXNUMPIX]; /* if >= 0, pixel i is - * swapped with - * iswap_partner_of_pix[i] - * in ireadrunfile, -pedfile, - * -calfile (internal - * numbering starting at 0)*/ - - Boolean bnewhistograms; /* set to true if new histograms shall be booked */ - Boolean bhisto_on; /* set to true if histogramming is on */ - - enum eERRORTOLERANCE eerrortolerance; - /* 0 == "cautious", exits on any reason (but tells in - * the .err file, - * 1 == "normal", exits when physics could be affected - * by the error, - * 2 == "tank", never exits except on coredumps and when - * all files have been processed. Do not use such files for - * physics analysis! - * - * NOTE: the capital letter words are the enums, the small letter - * words are to be used inside the parameter file. */ - Boolean bmontecarlo; - /* if TRUE we have a monte carlo dataset before us, if FALSE - * (default value) it is just normal data. */ - enum eMCTRIGGERFLAG emctriggerflag; - /* all: all events which survive the filter are written to the - * events NTuple. - * flag: When Dorota's triggerflag is set to 1 for a particular - * event, it shall be written to the output. All others shall - * just be disregarded. (Default) - * noflag: Opposite of 'flag': only events with triggerflag = 0 shall - * be treated further. */ - char mcfile[FILENAME_MAX]; - float fmean_nsb_pe_pix[iMAXNUMPIX]; - Boolean bfitmypedrms; - /* If TRUE, the pedestal MEAN and SIGMA=sqrt(MEAN) is fitted - * from pickup corrected pedestal data and not calculated from - * the pickup corrected ADC values, as it was done up to - * ~November 2000. */ -}; - -struct calibrationdata { /* Calibration (Laser Run) Data to be read - * from the Te?Cal files */ - int iid; /* no useful information, is long in the Macintosh format */ - int itelescope; /* number of the CT which took the data */ - int irunnum; /* run number from the file name */ - enum onoroff eruntype; /* indicates if the run is on- or off - * (taken from filename) */ - int inumruns; /* number of laser runs recorded in this file, is - * long in the Macintosh format */ - char cname[iMAXFILENAMELENGTH]; /* copy of the filename used when - * the data was read */ - double dmean_adc[iMAXNUMLASERRUNS][iMAXNUMPIX]; - /* means of pixel signals over all shots calculated by TIJARAFE / - * Roberto: in processdata.c this number is recalculated. */ - double dmean_lasoutput[iMAXNUMLASERRUNS][iMAXNUMLASEREVENTS]; - /* This new (April 26 2000) quantity monitors the laser mean - * output recorded in events, which are not considered cosmics. */ - double dmean_alladc[iMAXNUMPIX]; - /* means of pixel signals over all triggers calculated by - * preproc. */ - double dmean_laserjitter[iMAXNUMPIX]; - /* Here we calc. the deviation from the actual lasergain in a - * single laser shot to the laser mean as averaged over the camera */ - double drms_adc[iMAXNUMPEDRUNS][iMAXNUMPIX]; - /* the standard deviations of the pixel signals */ - int inumpixels; /* number of pixels in the camera, short in - * Macintosh format */ - int inumevents[iMAXNUMLASERRUNS]; /* number laser shots, - * short in Macintosh format */ - int ipixsig_adc[iMAXNUMLASERRUNS][iMAXNUMLASEREVENTS][iMAXNUMPIX]; - /* the signals obtained from each pixel for each laser shot and - * each laser run, unsigned short in the Macintosh format */ -}; - -struct pedestaldata { /* Pedestal Run Data to be read from the - * Te?Ped files */ - int iid; /* no useful information, is long in the Macintosh - * format */ - int itelescope; /* number of the CT which took the data */ - int irunnum; /* run number from the file name */ - enum onoroff eruntype; /* indicates if the run is on- or of - * (taken from filename) */ - int inumruns; /* number of pedestal events recorded in this file, - * is long in the Macintosh format */ - char cname[iMAXFILENAMELENGTH]; /* copy of the filename used - * when the data was read */ - double dmean_adc[iMAXNUMPEDRUNS][iMAXNUMPIX]; - /* means of pixel signals over all triggers calculated by - * TIJARAFE using a sigma cut of 2 */ - double dmean_alladc[iMAXNUMPIX]; - /* means of pixel signals over all triggers calculated by - * preproc. */ - double dmean_allPE[iMAXNUMPIX]; - /* means of pixel signals over all triggers calculated by - * preproc. */ - double dmean_pedlevel[iMAXNUMPEDRUNS][iMAXNUMPEDEVENTS]; - /* This new (preproc_0.4) quantity monitors the pedestal mean - * output recorded in events, which are not considered cosmics. - * The calc.ing is done inside processdata.c */ - double dmean_pedlevelPE[iMAXNUMPEDRUNS][iMAXNUMPEDEVENTS]; - /* This new (preproc_0.4) quantity monitors the pedestal mean - * output recorded in events, which are not considered cosmics. - * The calc.ing is done inside processdata.c */ - double dmean_pedoffset[iMAXNUMPIX]; - /* Mean Offset from the ADC spectrum of a pixel to the mean of the - * pedestal events. */ - double dmean_pedoffsetPE[iMAXNUMPIX]; - /* Mean Offset from the ADC spectrum of a pixel to the mean of the - * pedestal events. */ - double drms_adc[iMAXNUMPEDRUNS][iMAXNUMPIX]; - /* the standard deviations of the pixel signals used in the - * sigma cut mentioned above */ - double drms_alladc[iMAXNUMPIX]; - /* The RMS value of the as-is ADC spectrum for all events inside a - * Pedestal file. */ - double drms_allPE[iMAXNUMPIX]; - /* The RMS value of the as-is ADC spectrum for all events inside a - * Pedestal file. */ - double drms_nopickup[iMAXNUMPIX]; - /* the standard deviations of the pixel signals' deviation from - * the events' mean (in ADC counts, of course...) --- - * This is the same as the RMS of the pickup corrected pedestals! */ - double drms_nopickupPE[iMAXNUMPIX]; - /* the standard deviations of the pixel signals' deviation from - * the events' mean (in ADC counts, of course...) --- - * This is the same as the RMS of the pickup corrected pedestals! */ - int inumpixels; /* number of pixels in the camera, short in - * Macintosh format */ - int inumevents[iMAXNUMPEDRUNS]; /* number of random triggers, - * short in Macintosh format */ - int ipixsig_adc[iMAXNUMPEDRUNS][iMAXNUMPEDEVENTS][iMAXNUMPIX]; - /* the signals obtained from each pixel for each random trigger - * and each pedestal event, short in the Macintosh format */ -}; - -struct runfiledata { /* Data read from Te?Run files */ - int iformatid; /* Format ID of the data, long in Macintosh format */ - int iversid; /* Version ID of the data, long in Macintosh format */ - int irevid; /* Revision ID..., long in Macintosh format */ - int itelescope; /* number of the CT which took the data */ - int irunnum; /* run number from the file name */ - int irunnumx; /* number read from the file, long in Macintosh format */ - int inumevents; /* number of events recorded in this run, long in - * Macintosh format */ - char cname[iMAXFILENAMELENGTH]; /* copy of the filename used when - * the data was read */ - int ieventheaderlength; /* number of bytes in the event header, - * long in Macintosh format */ - int ieventheaderoffset; /* unknown meaning, long in Macintosh - * format */ - int isignalslength; /* number of bytes used by the pixel signals, - * long in Macintosh format */ - int isignalsoffset; /* unknown meaning, long in Macintosh format */ - int itdccountslength; /* number of bytes used by the TDC signals, - * long in Macintosh format */ - int itdccountsoffset; /* unknown meaning, long in Macintosh format */ - int iobservtime_mins; /* (planned) observational time (minutes), - * short in Macintosh format */ - enum onoroff eruntype; /* indicates if the run is on- or off-source; - * in the Macintosh format, this is a short - * where 1 means off and 2 means on */ - int ileapsecs; /* Leap seconds, introduced into UTC in order to - * keep track of some very small differences - * accumulated through the years. */ - int ihv1; - int ihv2; /* these two numbers represent the HV adjustment made - * by the shiftpersonnel. The last three items are - * only read for Roberto files, but not for Tijarafe - * files, as e.g. the HVs are not written out correctly. */ - char cdummy[40]; /* unknown meaning */ - int imactime_secs[iMAXNUMRUNEVENTS]; - /* Macintosh time (time_t mydatum) in seconds since 1970, this is - * put together from two Macintosh shorts and converted to the time - * format required by SYSTEM */ - unsigned int uigpstime_secs[iMAXNUMRUNEVENTS]; - /* GPS time (time_t mydatum) in seconds since 1970, this is - * put together from two Macintosh shorts and converted to the time - * format required by SYSTEM - * This time is new from Roberto and is coming from the PC-GPS card. - * There is another flag available in the runheader, which gives an - * overall flag for the run, whether or not the GPS time had been - * working. It is, however, not a good idea to use this flag, as - * the workability of this card can also be dedued from the GPS time - * itself. This allows to check for failures on a event-by-event - * basis, which seems more reasonable than using the global flag. */ - unsigned int uigpstime_usecs[iMAXNUMRUNEVENTS]; - /* GPS microseconds -- available from the card, that's why it is - * put to the data, here it is. - * This is put together from two Macintosh shorts and converted - * to the time format required by SYSTEM */ - unsigned int uirubsecs_secs[iMAXNUMRUNEVENTS]; - /* the seconds of the Rubidium clock time this is put together - * from two Macintosh shorts and converted to the time format - * required by SYSTEM */ - unsigned int uirubsecfrac_200ns[iMAXNUMRUNEVENTS]; - /* the fractional part of the Rubidium clock time second, this is - * put together from two Macintosh shorts, unit is 200 ns */ - int isepos[2][iMAXNUMRUNEVENTS]; - /* positions of shaft encoders 1 and 2, short in Macintosh format */ - float fhourangle; /* This is the angle between the observation of this - * event and the culmination point. It is going to - * be written into the events NTuple. */ - int ieventnumber[iMAXNUMRUNEVENTS]; -/* number of the event, short in Macintosh format */ - int inumpixels; /* number of pixels in the camera, short in - * Macintosh format */ - int inummuonpixels; /* number of pixels in the muon shield; - * the sum of inumpixels and inummuonpixels is - * part of the event record but should be the - * same for all events, so it is not put into - * an array */ - int inumcointdcs; /* number of coincidence TDCs for which counts - * were recorded, short in Macintosh format - * this value is part of the event record but - * should be the same for all events, so it is - * not put into an array */ - int ipixsig_adc[iMAXNUMRUNEVENTS][iMAXNUMPIX]; - /* signals from the camera photo multipliers (ADC counts) */ - int imuonpixsig_adc[iMAXNUMRUNEVENTS][iMAXNUMMUONPIX]; - /* signals from the muon shield photo multipliers (ADC counts) */ - int itdcsig_cnts[iMAXNUMRUNEVENTS][iMAXNUMCOINTDCS]; - /* counts from the coincidence TDCs if itdccountslength > 0 */ -}; - -struct TELLOGDATA { - int irunnum; /* The runnumber which we are interested in. - * Only data for this runnumber are written into - * this structure and then copied to the calnt NTuple. - * - * This runnumber is generally the closest match - * inside the TelLog file. */ - float pixel_timeaveraged_current[127]; - float pixel_timeaveraged_scaler[127]; -}; - -/************************************************ - * structures in output file * - * (Konopelko files format) * - * (from structures.h file 4.4) * - ************************************************/ -struct mcbankheader { /* Once in every file */ - int iparticle_type; - /* Primary particle type: 0 = gamma; +-1= e+-; - * +-2 = mu+-; +-3 = p+-, n; +-4 = pi+- */ - int inum_cts; - /* Number of simulated telescopes */ - int inum_nfl_shifts; - /* Number of NFL shifts for hadrons (per event) */ - int itrigger_thresh_phot; - /* Only images with two pixels above this threshold are - * in the MC data file */ - int imin_signal; - /* "Tail cut" = smallest signal saved */ - int inum_sim_showers; - /* Total number of showers simulated for this file */ - float fatmothick_gcm2; - /* Atmosphere above the telescope [g/cm2] */ - float fcos_zenangle; - /* Cosine of the zenith angle */ - float fnfl_shift_radius_deg; - /* Maximum angular radius within which a shifted - * hadron shower is sampled [degrees] */ - float fenergybin_bound_tev[iNUMENERGYBINS+1]; - /* Boundaries for the 14 energy bins [TeV] - * (newline after the first 8 numbers) */ - float fzenithanglebin_bound[iNUMZENANGLEBINS+1]; - /* Boundaries for the 11 zenith angle bins [degrees] */ - int inum_show_in_bin[iNUMENERGYBINS]; - /* Number of simulated (or saved) showers - * per energy bin (newline after the first 8 numbers) */ - float fmaxdist_impact_m[iNUMIMPACTBINS]; - /* Maximum distance of the impact point to the - * central telescope [m] for each energy bin - * (newline after the first 8 numbers) */ - int inumpixels_for_ct[iMAXNUMCTS]; - /* Number of pixels in the camera of each simulated CT */ - float fpixwidth_deg_ct[iMAXNUMCTS]; - /* Pixel width [degrees] for each CT */ - float fmirrsize_m2_ct[iMAXNUMCTS]; - /* Mirror area [m^2] for each CT */ - float fmean_nsb_phot_ct[iMAXNUMCTS]; - /* Mean signal caused by the NSB in each pixel for each CT - * [photoelectrons]. This is the simulation of the NSB, - * not the electronic noise */ - float fphes_per_photon_ct[iMAXNUMCTS]; - /* Conversion factor photoelectron per photon */ - float frelative_x_ct[iMAXNUMCTS]; - /* x position relative to the central CT for each CT */ - float frelative_y_ct[iMAXNUMCTS]; - /* y position relative to the central CT for each CT */ -}; -#endif - - - - - - - Index: trunk/MagicSoft/Mars/mmain/MOnlineDisplay.cc =================================================================== --- trunk/MagicSoft/Mars/mmain/MOnlineDisplay.cc (revision 4459) +++ trunk/MagicSoft/Mars/mmain/MOnlineDisplay.cc (revision 4460) @@ -69,5 +69,4 @@ #include "MHillasCalc.h" // MHillasCalc #include "MHillasSrcCalc.h" // MHillasSrcCalc -#include "MBlindPixelCalc.h" // MBlindPixelCalc #include "MFillH.h" // MFillH #include "MGTask.h" // MGTask Index: trunk/MagicSoft/Mars/status.cc =================================================================== --- trunk/MagicSoft/Mars/status.cc (revision 4459) +++ trunk/MagicSoft/Mars/status.cc (revision 4460) @@ -49,5 +49,5 @@ #include "MHillasCalc.h" #include "MHillasSrcCalc.h" -#include "MCT1SupercutsCalc.h" +//#include "MCT1SupercutsCalc.h" #include "MHCamEvent.h" #include "MFillH.h"