/* ======================================================================== *\ ! ! * ! * 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 Qi Zhe, 06/2007 ! Copyright: Software Development, 2000-2009 ! ! \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// // // MCorsikaRunHeader // // Root storage container for the RUN HEADER information // // Class Version 2: // ---------------- // + UInt_t fParticleID // + Float_t fImpactMax // //////////////////////////////////////////////////////////////////////////// #include "MCorsikaRunHeader.h" #include #include #include "MLog.h" #include "MLogManip.h" #include "MMcEvt.hxx" ClassImp(MCorsikaRunHeader); using namespace std; // -------------------------------------------------------------------------- // // Default constructor. Creates array which stores the pixel assignment. // // MCorsikaRunHeader::MCorsikaRunHeader(const char *name, const char *title) : fNumObsLevel(0), fZdMin(0), fZdMax(-1), fAzMin(0), fAzMax(0), fViewConeInnerAngle(0), fViewConeOuterAngle(-1) { fName = name ? name : "MCorsikaRunHeader"; fTitle = title ? title : "Raw Run Header Information"; } // -------------------------------------------------------------------------- // // Read in one run header from the binary file // Bool_t MCorsikaRunHeader::ReadEvt(istream& fin) { char runh[4]; fin.read(runh, 4); if (memcmp(runh, "RUNH", 4)) { *fLog << err << "ERROR - Wrong identifier: RUNH expected." << endl; return kFALSE; } Float_t f[68]; fin.read((char*)f, 68); fRunNumber = TMath::Nint(f[0]); fNumEvents = 0; fRunStart.SetCorsikaTime(f[1]); fProgramVersion = f[2]; //FIXME: INT??? fNumObsLevel = TMath::Nint(f[3]); if (fNumObsLevel!=1) { *fLog << err << "ERROR - Currently only one observation level is allowed." << endl; return kFALSE; } memset(fObsLevel, 0, 10*4); memcpy(fObsLevel, f+4, fNumObsLevel*4); fSlopeSpectrum = f[14]; fEnergyMin = f[15]; fEnergyMax = f[16]; fin.seekg(1020, ios::cur); // skip the remaining data of this block // -------------------- Read first event header ------------------- // FIXME: Add sanity checks! // f[76] Cherenkov flag: // bit(1) : CERENKOV option compiled in // bit(2) : IACT option compiled in // bit(3) : CEFFIC option compiled in // bit(4) : ATMEXT option compiled in // bit(5) : ATMEXT option used with refraction enabled // bit(6) : VOLUMEDET option compiled in // bit(7) : CURVED option compiled in // bit(9) : SLATN option compiled in // 11-21 : table number for externam athmosphere (but<1024) // // f[78] Curved athmosphere? (0=flat, 1=curved) // f[84] cherenkov bunch size // f[93] flag for additinal muon information of particle output file // f[145] Muon multiple scattering flag // f[156] altitude of horizontal shower axis char evth[4]; fin.read(evth, 4); if (memcmp(evth, "EVTH", 4)) { *fLog << err << "ERROR - Wrong identifier: EVTH expected." << endl; return kFALSE; } Float_t g[273]; fin.read((char*)&g, 273*4); if (fin.eof()) return kFALSE; fin.seekg(-274*4, ios::cur); const Int_t n = TMath::Nint(g[96]); // Number i of uses of each cherenkov event if (n!=1) { *fLog << err << "ERROR - Currently only one impact parameter per event is supported." << endl; return kFALSE; } fParticleID = TMath::Nint(g[1]); fImpactMax = g[86]; fZdMin = g[79]; // lower edge of theta in ° fZdMax = g[80]; // upper edge of theta in ° fAzMin = 180-g[81]; // lower edge of phi in ° fAzMax = 180-g[82]; // upper edge of phi in ° // FIXME: Correct for direction of magnetic field! fWavelengthMin = g[95]; // Cherenkov bandwidth lower end in nm fWavelengthMax = g[96]; // Cherenkov bandwidth upper end in nm fViewConeInnerAngle = g[151]; // inner angle of view cone (°) fViewConeOuterAngle = g[152]; // outer angle of view cone (°) return kTRUE; } Bool_t MCorsikaRunHeader::ReadEvtEnd(istream& fin) { char runh[4]; fin.read(runh, 4); if (memcmp(runh, "RUNE", 4)) { *fLog << err << "ERROR - Wrong identifier: RUNE expected." << endl; return kFALSE; } Float_t f[2]; fin.read((char*)f, 2*4); const UInt_t runnum = TMath::Nint(f[0]); if (runnum!=fRunNumber) { *fLog << err << "ERROR - Mismatch in stream: Run number in RUNE ("; *fLog << runnum << ") doesn't match RUNH (" << fRunNumber << ")." << endl; return kFALSE; } fNumEvents = TMath::Nint(f[1]); fin.seekg(270*4, ios::cur); // skip the remaining data of this block return kTRUE; } Bool_t MCorsikaRunHeader::SeekEvtEnd(istream &fin) { // Search subblockwise backward (Block: 5733*4 = 21*273*4) for (int i=1; i<22; i++) { fin.seekg(-i*273*4, ios::end); char runh[4]; fin.read(runh, 4); if (!memcmp(runh, "RUNE", 4)) { fin.seekg(-4, ios::cur); return kTRUE; } } return kFALSE; } // -------------------------------------------------------------------------- // // print run header information on *fLog. The option 'header' supresses // the pixel index translation table. // void MCorsikaRunHeader::Print(Option_t *t) const { *fLog << all << endl; *fLog << "Run Number: " << fRunNumber << " (" << fRunStart.GetStringFmt("%d.%m.%Y") << ", V" << fProgramVersion << ")" << endl; *fLog << "Particle ID: " << MMcEvt::GetParticleName(fParticleID) << endl; if (fNumEvents>0) *fLog << "Num Events: " << fNumEvents << endl; *fLog << "Obs Level: "; for (Byte_t i=0; i0) *fLog << "ViewCone: " << fViewConeInnerAngle << "deg-" << fViewConeOuterAngle << "deg" << endl; if (fZdMax>=0) { *fLog << "Zd/Az: " << fZdMin << "deg"; if (fZdMin==fZdMax) *fLog << " (fixed)"; else *fLog << "-" << fZdMax << "deg"; *fLog << " / " << fAzMin << "deg"; if (fAzMin==fAzMax) *fLog << " (fixed)"; else *fLog << "-" << fAzMax << "deg"; *fLog << endl; } }