Index: trunk/MagicSoft/Simulation/Detector/include-MC/MCEventHeader_2.hxx =================================================================== --- trunk/MagicSoft/Simulation/Detector/include-MC/MCEventHeader_2.hxx (revision 1562) +++ trunk/MagicSoft/Simulation/Detector/include-MC/MCEventHeader_2.hxx (revision 1562) @@ -0,0 +1,338 @@ +///////////////////////////////////////////////////////////////// +// +// MCEventHeader_2 +// +// Created: Sat Oct 19 12:11:43 CEST 2002 +// Author: Oscar Blanch Bigas +// Purpose: Base class for reflector 0.6 EventHeader +// Notes: +// +///////////////////////////////////////////////////////////////// + +// @T \newpage + +// @section Source code of {\tt MCEventHeader_2.hxx} + +/* @text +This section shows the include file {\tt MCEventHeader_2.hxx} +@endtext */ + +#ifndef MCEventHeader_2_Class +#define MCEventHeader_2_Class + +// @subsection Include files + +// @code +#ifdef __ROOT__ +#include "TROOT.h" +#include "TObject.h" +#else // not __ROOT__ +#include "Rtypes.h" +#endif + +#include +#include +#include +#include +#include + +#include "COREventHeader.hxx" + +#define SIZE_OF_MCEVENTHEADER_2 178 /* floats */ + +// @endcode + +// @subsection Class {\em MCEventHeader_2}: Definition + +// @code +class MCEventHeader_2 { + +protected: + Float_t EvtNumber; + Float_t PrimaryID; + Float_t Etotal; + Float_t Thick0; + Float_t FirstTarget; + Float_t zFirstInt; + Float_t p[3]; + Float_t Theta; + Float_t Phi; + + Float_t NumRndSeq; + Float_t RndData[10][3]; + + Float_t RunNumber; + Float_t DateRun; + Float_t VersionPGM; + + Float_t NumObsLev; + Float_t HeightLev[10]; + + Float_t SlopeSpec; + Float_t ELowLim; + Float_t EUppLim; + + Float_t Ecutoffh; + Float_t Ecutoffm; + Float_t Ecutoffe; + Float_t Ecutoffg; + + Float_t NFLAIN; + Float_t NFLDIF; + Float_t NFLPI0; + Float_t NFLPIF; + Float_t NFLCHE; + Float_t NFRAGM; + + Float_t Bx; + Float_t By; + + Float_t EGS4yn; + Float_t NKGyn; + Float_t GHEISHAyn; + Float_t VENUSyn; + Float_t CERENKOVyn; + Float_t NEUTRINOyn; + Float_t HORIZONTyn; + Float_t COMPUTER; + + Float_t ThetaMin; + Float_t ThetaMax; + Float_t PhiMin; + Float_t PhiMax; + + Float_t CBunchSize; + Float_t CDetInX,CDetInY; + Float_t CSpacInX,CSpacInY; + Float_t CLenInX,CLenInY; + Float_t COutput; + + Float_t AngleNorthX; + Float_t MuonInfo; + + Float_t StepLength; + Float_t CWaveLower; + Float_t CWaveUpper; + Float_t Multipl; + Float_t CorePos[2][20]; + Float_t SIBYLL[2]; + Float_t QGSJET[2]; + Float_t DPMJET[2]; + Float_t VENUS_cross; + Float_t mu_mult_scat; + Float_t NKG_range; + Float_t EFRCTHN[2]; + Float_t WMAX[2]; + Float_t rthin_rmax; + + Float_t viewcone_angles[2]; + /* (degrees) Inner and outer angles in + * Corsika's VIEWCONE option. This is + * only possible with Corsika>6 versions. In + * that case, PhiMin=PhiMax and + * ThetaMin=ThetaMax (also in this header) + * indicate the axis of this cone. + */ + + /* ^^^ Up to here, the info from the CORSIKA event header. */ + + /* Telescope orientation: */ + Float_t telescopePhi; /* rad */ + Float_t telescopeTheta; /* rad */ + + /* Time of first and last photon: */ + Float_t TimeFirst; + Float_t TimeLast; + + /* 6 parameters and chi2 of the NKG fit to the longitudinal + * particle distribution (see CORSIKA manual for explanation): + */ + Float_t longi_Nmax; + Float_t longi_t0; + Float_t longi_tmax; + Float_t longi_a; + Float_t longi_b; + Float_t longi_c; + Float_t longi_chi2; + + Float_t CORSIKAPhs; // Original photons written by Corsika + Float_t AtmAbsPhs; // Photons absorved by the atmosphere + Float_t MirrAbsPhs; // Photons absorved by the mirrors + Float_t OutOfMirrPhs; // Photons outside the mirror + Float_t BlackSpotPhs; // Photons lost in the "black spot" + Float_t OutOfChamPhs; // Photons outside the chamber + Float_t CPhotons; // Photons reaching the chamber + + Float_t dummy[10]; /* not used */ + +public: + MCEventHeader_2(void) {} // default constructor + + virtual ~MCEventHeader_2(void) {} // default destructor + + // reads EventHeader from binary input stream + Int_t read ( ifstream &is ) { + is.read ( (char *)this, mysize() ); + return is.gcount(); + } + + // writes EventHeader to binary output stream + Int_t write ( ofstream &os ) { + os.write ( (char *)this, mysize() ); + return 0; + } + + // get information about the EventHeader + + // get the primary type (GEANT code) + inline Float_t get_primary ( void ) { return ( PrimaryID ); } + + // get the total primary energy + inline Float_t get_energy ( void ) { return ( Etotal ); } + + // get the initial zenith angle + inline Float_t get_theta ( void ) { return ( Theta ); } + + // get the initial phi angle + inline Float_t get_phi ( void ) { return ( Phi ); } + + // get the telescope zenith angle + inline Float_t get_theta_CT ( void ) { return ( telescopeTheta ); } + + // get the telescope phi angle + inline Float_t get_phi_CT ( void ) { return ( telescopePhi ); } + + // get the spectrum slope + inline Float_t get_slope ( void ) { return ( SlopeSpec ); } + + // get height of first interaction (in cm) + inline Float_t get_height ( void ) { return ( zFirstInt ); } + + // get the energy range of this run + inline void get_energy_range ( Float_t *elow, Float_t *eup ) { + *elow = ELowLim; + *eup = EUppLim; + } + + // get the theta range of this run + inline void get_theta_range ( Float_t *thlow, Float_t *thup ) { + *thlow = ThetaMin; + *thup = ThetaMax; + } + + // get the energy range of this run + inline void get_phi_range ( Float_t *plow, Float_t *pup ) { + *plow = PhiMin; + *pup = PhiMax; + } + + // get the core position + inline Float_t get_core ( Float_t *x, Float_t *y, Int_t ncore = 0 ) { + *x = CorePos[0][ncore]; + *y = CorePos[1][ncore]; + return ( (Float_t) sqrt ( (*x)*(*x) + (*y)*(*y) ) ); + } + + // get the core position + inline Float_t get_core ( Int_t ncore = 0 ) { + return ( (Float_t) sqrt ( ((CorePos[0][ncore])*(CorePos[0][ncore]))+ + ((CorePos[1][ncore])*(CorePos[1][ncore])) ) ); + } + + // get the core position in X + inline Float_t get_coreX ( Int_t ncore = 0 ) { + return ( (Float_t) CorePos[0][ncore] ); + } + + // get the core position in Y + inline Float_t get_coreY ( Int_t ncore = 0 ) { + return ( (Float_t) CorePos[1][ncore] ); + } + + //get photons at different stages + inline Float_t get_CORSIKA (){ return CORSIKAPhs;} + + inline Float_t get_AtmAbs (){ return AtmAbsPhs;} + + inline Float_t get_MirrAbs (){ return MirrAbsPhs;} + + inline Float_t get_OutOfMirr (){ return OutOfMirrPhs;} + + inline Float_t get_BlackSpot (){ return BlackSpotPhs;} + + inline Float_t get_OutOfCham (){ return OutOfChamPhs;} + + inline Float_t get_CPhotons (){ return CPhotons;} + + inline Float_t get_NumObsLev (){return NumObsLev;} + + inline Float_t get_HeightLev (Int_t i) {return HeightLev[i];} + + inline Float_t get_VersionPGM () {return VersionPGM;} + + inline Float_t get_RunNumber () {return RunNumber;} + + inline Float_t get_DateRun () {return DateRun;} + + inline Float_t get_NumRndSeq () {return NumRndSeq;} + + inline Float_t get_CWaveLower() {return CWaveLower;} + + inline Float_t get_CWaveUpper() {return CWaveUpper;} + + // transport from COREventHeader to MCEventHeader_2 + void transport ( COREventHeader *e ); + + // write extreme times + inline void put_times ( Float_t t1, Float_t t2 ) { + TimeFirst = t1; + TimeLast = t2; + } + + // get extreme times + inline Float_t get_times ( Float_t *t1, Float_t *t2 ) { + *t1 = TimeFirst; + *t2 = TimeLast; + return ( TimeLast - TimeFirst ); + } + + inline int mysize(void) + { return ( sizeof( float ) * SIZE_OF_MCEVENTHEADER_2 ); } + + // get deviations of the CT from the original shower direction + inline Float_t get_deviations ( Float_t *t1, Float_t *t2 ) { + + float ct1,st1,cp1,sp1; + float ct2,st2,cp2,sp2; + float x1,y1,z1; + float x2,y2,z2; + + ct1 = cos(Theta); + st1 = sin(Theta); + cp1 = cos(Phi); + sp1 = sin(Phi); + + ct2 = cos(telescopeTheta); + st2 = sin(telescopeTheta); + cp2 = cos(telescopePhi); + sp2 = sin(telescopePhi); + + x1 = st1*cp1; y1 = st1*sp1; z1 = ct1; + x2 = st2*cp2; y2 = st2*sp2; z2 = ct2; + + *t1 = (telescopeTheta-Theta); + *t2 = (telescopePhi-Phi); + + return ( acos(x1*x2 + y1*y2 + z1*z2) ); + } + + inline void print ( void ) { + float *ptr = (float *)this; + for(int i=0; i