/* ======================================================================== *\ ! ! * ! * 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, 2008 ! ! \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// // // MCorsikaEvtHeader // ///////////////////////////////////////////////////////////////////////////// #include "MCorsikaEvtHeader.h" #include #include #include "MLog.h" #include "MLogManip.h" #include "MMcEvt.hxx" ClassImp(MCorsikaEvtHeader); using namespace std; // -------------------------------------------------------------------------- // // Default constructor. Create the array to store the data. // MCorsikaEvtHeader::MCorsikaEvtHeader(const char *name, const char *title) { fName = name ? name : "MCorsikaEvtHeader"; fTitle = title ? title : "Raw Event Header Information"; } // -------------------------------------------------------------------------- // // This member function prints all Data of one Event to *fLog. // void MCorsikaEvtHeader::Print(Option_t *o) const { *fLog << all; *fLog << "Event Number: " << dec << fEvtNumber << endl; *fLog << "Particle ID: " << MMcEvt::GetParticleName(fParticleID) << endl; *fLog << "Energy: " << fTotalEnergy << "GeV" << endl; *fLog << "Starting Altitude: " << fStartAltitude << "g/cm²" << endl; *fLog << "Number of 1st Target: " << fFirstTargetNum << endl, *fLog << "Height of 1st Interaction: " << fFirstInteractionHeight/100. << "m" << endl; *fLog << "Momentum X/Y/Z (GeV/c): " << fMomentumX << "/" << fMomentumY << "/" << fMomentumZ << endl; *fLog << "Zenith/Azimuth Angle: " << fZd*TMath::RadToDeg() << "°/" << fAz*TMath::RadToDeg() << "°" << endl; *fLog << "Impact X/Y: " << fX/100. << "m/" << fY/100. << "m (r=" << TMath::Hypot(fX, fY)/100. << "m)" << endl; *fLog << "Weighted Num Photons: " << fWeightedNumPhotons << endl; *fLog << endl; } // -------------------------------------------------------------------------- // // read the EVTH information from the input stream // return FALSE if there is no header anymore, else TRUE // Int_t MCorsikaEvtHeader::ReadEvt(std::istream &fin) { char evth[4]; fin.read(evth, 4); if (memcmp(evth, "EVTH", 4)) { fin.seekg(-4, ios::cur); return kFALSE; } Float_t f[273]; fin.read((char*)&f, 273*4); fEvtNumber = TMath::Nint(f[0]); fParticleID = TMath::Nint(f[1]); fTotalEnergy = f[2]; fStartAltitude = f[3]; fFirstTargetNum = f[4]; fFirstInteractionHeight = f[5]; // Pointing opposite particle direction // (x=north, y=west, z=upwards) //fMomentumX = f[6]; //fMomentumY = f[7]; //fMomentumZ = f[8]; // Pointing along particle direction // (x=east, y=north, z=upwards) fMomentumX = -f[7]; fMomentumY = f[6]; fMomentumZ = -f[8]; fZd = f[9]; fAz = TMath::Pi()-f[10]; const Int_t n = TMath::Nint(f[96]); if (n!=1) { *fLog << err << "ERROR - Currently only one impact parameter per event is supported." << endl; return kFALSE; } //fX = f[97]; //fY = f[117]; fX = f[117]; fY = -f[97]; /* if (fEvtNumber==1) { header.fZdMin = f[79]; header.fZdMax = f[80]; header.fAzMin = 180-f[81]; header.fAzMax = 180-f[82]; header.fViewConeInnerAngle = f[151]; header.fViewConeOuterAngle = f[152]; header.Print(); } */ fin.seekg(1088-273*4, ios::cur); return !fin.eof(); } // this member function is for reading the event end block Bool_t MCorsikaEvtHeader::ReadEvtEnd(std::istream &fin) { //fin.seekg(-1088,ios::cur); Float_t f[2]; fin.read((char*)&f, 2*4); const UInt_t evtnum = TMath::Nint(f[0]); if (evtnum!=fEvtNumber) { *fLog << err << "ERROR - Mismatch in stream: Event number in EVTE ("; *fLog << evtnum << ") doesn't match EVTH (" << fEvtNumber << ")." << endl; return kERROR; } fWeightedNumPhotons = f[1]; fin.seekg(1080,ios::cur); return !fin.eof(); }