/* ======================================================================== *\ ! ! * ! * 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/2008 ! ! Copyright: Software Development, 2000-2008 ! ! \* ======================================================================== */ ////////////////////////////////////////////////////////////////////////////// // // MCorsikaRead // // Input Containers: // -/- // // Output Containers: // MCorsikaRunHeader // MCorsikaEvtHeader // MPhotonEvent // ////////////////////////////////////////////////////////////////////////////// #include "MCorsikaRead.h" #include #include #include #include "MLog.h" #include "MLogManip.h" #include "MParList.h" #include "MStatusDisplay.h" #include "MCorsikaFormat.h" #include "MCorsikaRunHeader.h" #include "MCorsikaEvtHeader.h" #include "MPhotonEvent.h" ClassImp(MCorsikaRead); using namespace std; /* ----------- please don't delete and don't care about (Thomas) ------------ #define kBUFSZ 64 //1024*1024*64 #include class bifstream : public istream, public streambuf { private: char fBuffer[kBUFSZ]; //! FILE *fd; int sync() { memset(fBuffer, 0, kBUFSZ); return 0; } int underflow() { int sz=fread(fBuffer, kBUFSZ, 1, fd); //int sz=fread(fBuffer, 1, kBUFSZ, fd); setg(fBuffer, fBuffer, fBuffer+kBUFSZ); return sz==1 ? *(unsigned char*)fBuffer : EOF;//EOF; //return sz==kBUFSZ ? *(unsigned char*)fBuffer : EOF;//EOF; } public: bifstream(const char *name) : istream(this) { fd = fopen(name, "rb"); setbuf(fBuffer, kBUFSZ); } }; */ // -------------------------------------------------------------------------- // // Default constructor. It tries to open the given file. // MCorsikaRead::MCorsikaRead(const char *fname, const char *name, const char *title) : fRunHeader(0), fEvtHeader(0), fEvent(0), fTelescopeIdx(-1), fForceMode(kFALSE), fFileNames(0), fNumFile(0), fNumEvents(0), fNumTotalEvents(0), /*fIn(0),*/ fInFormat(0), fParList(0), fNumTelescopes(1) { fName = name ? name : "MRead"; fTitle = title ? title : "Read task to read DAQ binary files"; fFileNames = new TList; fFileNames->SetOwner(); if (fname!=NULL) AddFile(fname); } // -------------------------------------------------------------------------- // // Destructor. Delete input stream. // MCorsikaRead::~MCorsikaRead() { delete fFileNames; // if (fIn) // delete fIn; if (fInFormat) delete fInFormat; } /* Byte_t MCorsikaRead::IsFileValid(const char *name) { MZlib fin(name); if (!fin) return 0; Byte_t c[4]; fin.read((char*)c, 4); if (!fin) return 0; if (c[0]!=0xc0) return 0; if (c[1]==0xc0) return 1; if (c[1]==0xc1) return 2; return 0; } */ // -------------------------------------------------------------------------- // // Add a new file to a list of files to be processed, Returns the number // of files added. (We can enhance this with a existance chack and // wildcard support) // Int_t MCorsikaRead::AddFile(const char *fname, Int_t entries) { TNamed *name = new TNamed(fname, ""); fFileNames->AddLast(name); return 1; } Bool_t MCorsikaRead::ReadEvtEnd() { if (!fInFormat->SeekEvtEnd()) { *fLog << (fForceMode?warn:err) << "Error: RUNE section not found in file." << endl; if (!fForceMode) return kFALSE; } if (!fRunHeader->ReadEvtEnd(fInFormat)) { *fLog << (fForceMode?warn:err) << "Error: Reading RUNE section failed." << endl; if (!fForceMode) return kFALSE; } return kTRUE; } // -------------------------------------------------------------------------- // // This opens the next file in the list and deletes its name from the list. // Int_t MCorsikaRead::OpenNextFile(Bool_t print) { // // open the input stream and check if it is really open (file exists?) // /* if (fIn) delete fIn; fIn = NULL; */ if (fInFormat) delete fInFormat; fInFormat = NULL; // // Check for the existance of a next file to read // TObject *file = fFileNames->At(fNumFile); 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); fInFormat = MCorsikaFormat::CorsikaFormatFactory(expname); delete [] expname; if (fInFormat == NULL) return kERROR; *fLog << inf << "Open file: '" << name << "'" << endl; if (fDisplay) { // Show the number of the last event after // which we now open a new file TString txt = GetFileName(); txt += " @ "; txt += GetNumExecutions()-1; fDisplay->SetStatusLine2(txt); } fNumFile++; // // Read RUN HEADER (see specification) from input stream // if (!fRunHeader->ReadEvt(fInFormat)) return kERROR; // if (!fEvtHeader->ReadRunHeader(*fIn, *fRunHeader)) // return kERROR; fNumTelescopes = 1; if (fInFormat->IsEventioFormat()) { fInFormat->StorePos(); fInFormat->Rewind(); if (!fInFormat->SeekNextBlock(0, 1201)) { *fLog << err << "ERROR - Object 1201 not found in file." << endl; return kERROR; } fInFormat->Read(&fNumTelescopes); if (fTelescopeIdx>=fNumTelescopes) { *fLog << err << "ERROR - Requested telescope index " << fTelescopeIdx << " exceeds number of telescopes " << fNumTelescopes << " in file." << endl; return kERROR; } fTelescopeX.Set(fNumTelescopes); fTelescopeY.Set(fNumTelescopes); fTelescopeZ.Set(fNumTelescopes); fTelescopeR.Set(fNumTelescopes); fInFormat->Read(fTelescopeX.GetArray(), 4*fNumTelescopes); fInFormat->Read(fTelescopeY.GetArray(), 4*fNumTelescopes); fInFormat->Read(fTelescopeZ.GetArray(), 4*fNumTelescopes); fInFormat->Read(fTelescopeR.GetArray(), 4*fNumTelescopes); *fLog << all; *fLog << "Number of telescopes: " << fNumTelescopes << " (using "; if (fTelescopeIdx>=0) *fLog << "telescope " << fTelescopeIdx; else *fLog << "all telescopes"; *fLog << ")" << endl; const Int_t lo = fTelescopeIdx<0 ? 0 : fTelescopeIdx; const Int_t up = fTelescopeIdx<0 ? fNumTelescopes : fTelescopeIdx+1; for (int i=lo; iResetPos(); } fInFormat->StorePos(); if (!ReadEvtEnd()) return kERROR; fInFormat->ResetPos(); fNumEvents += fRunHeader->GetNumEvents(); fRunHeader->SetReadyToSave(); // // Print Run Header // We print it after the first event was read because // we still miss information which is stored in the event header?!? if (print) fRunHeader->Print(); if (!fParList) return kTRUE; // // Search for MTaskList // MTask *tlist = (MTask*)fParList->FindObject("MTaskList"); if (!tlist) { *fLog << err << dbginf << "MTaskList not found... abort." << endl; return kERROR; } // // A new file has been opened and new headers have been read. // --> ReInit tasklist // return tlist->ReInit(fParList) ? kTRUE : kERROR; } // -------------------------------------------------------------------------- // // Return file name of current file. // TString MCorsikaRead::GetFullFileName() const { const TObject *file = fFileNames->At(fNumFile-1); return file ? file->GetName() : ""; } // -------------------------------------------------------------------------- // // Restart with the first file // Bool_t MCorsikaRead::Rewind() { fNumFile=0; fNumEvents=0; return OpenNextFile()==kTRUE; } Bool_t MCorsikaRead::CalcNumTotalEvents() { fNumTotalEvents = 0; Bool_t rc = kTRUE; while (1) { switch (OpenNextFile(kFALSE)) { case kFALSE: break; case kERROR: rc = kFALSE; break; case kTRUE: if (!ReadEvtEnd()) { rc = kFALSE; break; } fNumTotalEvents += fRunHeader->GetNumEvents()*fRunHeader->GetNumReuse()* (fTelescopeIdx<0 ? fNumTelescopes : 1); continue; } break; } /* if (fIn) delete fIn; fIn = NULL; */ return rc; } // -------------------------------------------------------------------------- // // The PreProcess of this task checks for the following containers in the // list: // MCorsikaRunHeader if not found it is created // MCorsikaEvtHeader if not found it is created // MCorsikaEvtData if not found it is created // MCorsikaCrateArray if not found it is created // MCorsikaEvtTime if not found it is created (MTime) // // If all containers are found or created the run header is read from the // binary file and printed. If the Magic-Number (file identification) // doesn't match we stop the eventloop. // // Now the EvtHeader and EvtData containers are initialized. // Int_t MCorsikaRead::PreProcess(MParList *pList) { // // open the input stream // first of all check if opening the file in the constructor was // successfull // fParList = 0; /* if (!OpenStream()) return kFALSE; */ // // check if all necessary containers exist in the Parameter list. // if not create one and add them to the list // fRunHeader = (MCorsikaRunHeader*)pList->FindCreateObj("MCorsikaRunHeader"); if (!fRunHeader) return kFALSE; fEvtHeader = (MCorsikaEvtHeader*)pList->FindCreateObj("MCorsikaEvtHeader"); if (!fEvtHeader) return kFALSE; fEvent = (MPhotonEvent*)pList->FindCreateObj("MPhotonEvent"); if (!fEvent) return kFALSE; *fLog << inf << "Calculating number of total events..." << flush; if (!CalcNumTotalEvents()) return kFALSE; *fLog << inf << " " << fNumTotalEvents << " found." << endl; fNumFile=0; fNumEvents=0; fParList = pList; return kTRUE; } // -------------------------------------------------------------------------- // // Read a single event from the stream // Int_t MCorsikaRead::ReadEvent() { if (fInFormat->IsEventioFormat()) { if (fNumTelescope==0) { const Int_t rc1 = fEvtHeader->ReadEvt(fInFormat); if (rc1!=kTRUE) return rc1; if (fEvtHeader->GetNumReuse()==fRunHeader->GetNumReuse()-1) fEvtHeader->ResetNumReuse(); fEvtHeader->IncNumReuse(); cout << "===== Array " << fEvtHeader->GetNumReuse() << " =====" << endl; } const Int_t rc2 = fEvent->ReadCorsikaEvt(fInFormat, fTelescopeIdx); if (rc2!=kTRUE) return rc2; fEvtHeader->InitXY(); const Double_t x = fTelescopeX[fNumTelescope]; const Double_t y = fTelescopeY[fNumTelescope]; fEvtHeader->AddXY(y, -x); fEvent->AddXY(fEvtHeader->GetX(), fEvtHeader->GetY()); fEvent->SimWavelength(fRunHeader->GetWavelengthMin(), fRunHeader->GetWavelengthMax()); fNumTelescope++; fNumTelescope%=fNumTelescopes; return kTRUE; } else { // Store the position to re-read a single event several times if (fEvtHeader->GetNumReuse()GetNumReuse()-1) fInFormat->ResetPos(); else { fInFormat->StorePos(); fEvtHeader->ResetNumReuse(); } // Read the event header const Int_t rc1 = fEvtHeader->ReadEvt(fInFormat); if (rc1!=kTRUE) return rc1; fEvtHeader->IncNumReuse(); fEvtHeader->InitXY(); // In the case of corsika raw data (MMCS only) we have too loop over one event // several times to read all impact parameters (reusage of events) // In case of the eventio format we have to decide, the data of which telescope // we want to extract // Read the photons corresponding to the i-th core location // EventIO: Number of telescope // Raw: Number of re-use const Int_t rc2 = fEvent->ReadCorsikaEvt(fInFormat, fEvtHeader->GetNumReuse()+1); if (rc2!=kTRUE) return rc2; // read event end return fEvtHeader->ReadEvtEnd(fInFormat); } } // -------------------------------------------------------------------------- // // The Process reads one event from the binary file: // - The event header is read // - the run header is read // - all crate information is read // - the raw data information of one event is read // Int_t MCorsikaRead::Process() { while (1) { if (fInFormat) { // Read a single event from file const Int_t rc = ReadEvent(); // kFALSE means: end of file (try next one) if (rc!=kFALSE) return rc; fInFormat->UnreadLastHeader(); if (!fRunHeader->ReadEvtEnd(fInFormat)) if (!fForceMode) return kERROR; } // // If an event could not be read from file try to open new file // const Int_t rc = OpenNextFile(); if (rc!=kTRUE) return rc; } return kTRUE; } // -------------------------------------------------------------------------- // // Close the file. Check whether the number of read events differs from // the number the file should containe (MCorsikaRunHeader). Prints a warning // if it doesn't match. // Int_t MCorsikaRead::PostProcess() { const UInt_t n = fNumEvents*fRunHeader->GetNumReuse()*(fTelescopeIdx<0?fNumTelescopes:1); // // Sanity check for the number of events // if (n==GetNumExecutions()-1 || GetNumExecutions()==0) return kTRUE; *fLog << warn << dec; *fLog << "Warning - number of read events (" << GetNumExecutions()-1; *fLog << ") doesn't match number in run header(s) ("; *fLog << fRunHeader->GetNumReuse() << "*" << fNumEvents << "=" << n << ")." << endl; return kTRUE; } // -------------------------------------------------------------------------- // // Telescope: 1 // // In the case of eventio-format select which telescope to extract from the // file. -1 will extract all telescopes // Int_t MCorsikaRead::ReadEnv(const TEnv &env, TString prefix, Bool_t print) { Bool_t rc = kFALSE; if (IsEnvDefined(env, prefix, "Telescope", print)) { rc = kTRUE; fTelescopeIdx = GetEnvValue(env, prefix, "Telescope", fTelescopeIdx); } return rc; }