Index: /trunk/MagicSoft/Mars/msim/MHPhotonEvent.cc =================================================================== --- /trunk/MagicSoft/Mars/msim/MHPhotonEvent.cc (revision 9232) +++ /trunk/MagicSoft/Mars/msim/MHPhotonEvent.cc (revision 9232) @@ -0,0 +1,362 @@ +/* ======================================================================== *\ +! +! * +! * 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, 1/2009 +! +! Copyright: Software Development, 2000-2009 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MHPhotonEvent +// +// This is a histogram class to visualize the contents of a MPhotonEvent. +// +// Histograms for the distribution in x/y (z is just ignored), the +// arrival direction, a profile of the arrival time vs position and +// a spectrum is filles and displayed. +// +// There are several types of histograms (in the sense of binnings) +// for several purposes: +// +// Type 1: +// The maximum in x and y is determined from MCorsikaRunHeader +// (not yet implemented. Fixed to 25000) +// +// Type 2: +// The maximum in x and y is determined from MReflector->GetMaxR(); +// +// Type 3: +// The maximum in x and y is determined from MCamera->GetMaxR(); +// +// Type 4: +// As type 3 but divided by 10. +// +// Type 5: +// The maximum in x and y is determined from MGeomCam->GetMaxR(); +// +// Type 6: +// As type 5 but divided by 10. +// +// The binning is optimized using MH::FindGoodLimits. The number of bins +// in 100 in the default case and 50 for type 3-6. +// +// Fill expects a MPhotonEvent (the second argumnet in MFillH). +// +///////////////////////////////////////////////////////////////////////////// +#include "MHPhotonEvent.h" + +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MBinning.h" +#include "MParList.h" + +#include "MGeomCam.h" + +#include "MPhotonEvent.h" +#include "MPhotonData.h" + +#include "MCorsikaRunHeader.h" + +#include "../msimreflector/MSimReflector.h" // MCamera, MReflector + +ClassImp(MHPhotonEvent); + +using namespace std; + +void MHPhotonEvent::Init(const char *name, const char *title) +{ + fName = name ? name : "MHPhotonEvent"; + fTitle = title ? title : "Histogram to display the information of MPhotonEvents"; + + fHistXY.SetName("Position"); + fHistXY.SetTitle("Histogram of position x/y"); + fHistXY.SetXTitle("X [cm]"); + fHistXY.SetYTitle("Y [cm]"); + fHistXY.SetZTitle("Counts"); + fHistXY.SetDirectory(NULL); + fHistXY.Sumw2(); + + fHistUV.SetName("Direction"); + fHistUV.SetTitle("Histogram of arrival direction CosU/CosV"); + fHistUV.SetXTitle("CosU"); + fHistUV.SetYTitle("CosV"); + fHistUV.SetZTitle("Counts"); + fHistUV.SetDirectory(NULL); + fHistUV.Sumw2(); + + fHistT.SetName("Time"); + fHistT.SetTitle("Time profile in x/y"); + fHistT.SetXTitle("X [cm]"); + fHistT.SetYTitle("Y [cm]"); + fHistT.SetZTitle("T [ns]"); + fHistT.SetDirectory(NULL); + + fHistWL.SetName("Spectrum"); + fHistWL.SetTitle("Wavelength distribution"); + fHistWL.SetXTitle("\\lamba [nm]"); + fHistWL.SetYTitle("Counts"); + fHistWL.SetDirectory(NULL); + + // FIXME: Get this information from the corsika run-header + MBinning(70, 275, 625).Apply(fHistWL); +} + +// -------------------------------------------------------------------------- +// +// Default constructor. Creates and initializes the histograms. +// +MHPhotonEvent::MHPhotonEvent(Double_t max, const char *name, const char *title) + : fHistT("", "", 1, 0, 1, 1, 0, 1), fType(-1), fPermanentReset(kFALSE) +{ + // pre-initialization of the profile is necessary to get fZmin and fZmax set + + Init(name, title); + + MBinning binsd, binsa; + binsd.SetEdges(50, -max, max); + binsa.SetEdges(50, -1, 1); + + SetBinning(&fHistXY, &binsd, &binsd); + SetBinning(&fHistUV, &binsa, &binsa); + SetBinning(&fHistT, &binsd, &binsd); +} + +// -------------------------------------------------------------------------- +// +// Creates and initializes the histograms. +// +MHPhotonEvent::MHPhotonEvent(Int_t type, const char *name, const char *title) + : fHistT("", "", 1, 0, 1, 1, 0, 1), fType(type), fPermanentReset(kFALSE) +{ + // pre-initialization of the profile is necessary to get fZmin and fZmax set + + Init(name, title); + + MBinning binsd, bins; + bins.SetEdges(50, -1, 1); + + SetBinning(&fHistUV, &bins, &bins); +} + +// -------------------------------------------------------------------------- +// +// Search for MRflEvtData, MRflEvtHeader and MGeomCam +// +Bool_t MHPhotonEvent::SetupFill(const MParList *pList) +{ + Double_t xmax = -1; + Int_t num = 100; + + switch (fType) + { + case -1: + return kTRUE; + // case0: Take a value defined by the user + case 1: + { + MCorsikaRunHeader *h = (MCorsikaRunHeader*)pList->FindObject("MCorsikaRunHeader"); + if (!h) + { + *fLog << err << "MCorsikaRunHeader not found... aborting." << endl; + return kFALSE; + } + xmax = 25000;//h->GetImpactMax()*1.25; // Estimate scattering in the atmosphere + break; + } + case 2: + { + MReflector *r = (MReflector*)pList->FindObject("MReflector"); + if (!r) + { + *fLog << err << "MReflector not found... aborting." << endl; + return kFALSE; + } + xmax = r->GetMaxR(); + break; + } + case 3: + case 4: + { + MCamera *c = (MCamera*)pList->FindObject("MCamera"); + if (!c) + { + *fLog << err << "MCamera not found... aborting." << endl; + return kFALSE; + } + xmax = fType==3 ? c->GetMaxR() : c->GetMaxR()/10; + num = 50; + + break; + } + case 5: + case 6: + { + MGeomCam *c = (MGeomCam*)pList->FindObject("MGeomCam"); + if (!c) + { + *fLog << err << "MGeomCam not found... aborting." << endl; + return kFALSE; + } + xmax = fType==5 ? c->GetMaxRadius() : c->GetMaxRadius()/10; + num = 50; + + break; + } + default: + *fLog << err << "No valid binning (Type=" << fType << ") given... aborting." << endl; + return kFALSE; + } + + Double_t xmin = -xmax; + + MH::FindGoodLimits(num, num, xmin, xmax, kFALSE); + + MBinning binsd, binsa, binsz; + binsd.SetEdges(num, xmin, xmax); + + SetBinning(&fHistXY, &binsd, &binsd); + SetBinning(&fHistT, &binsd, &binsd); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Fill contents of MPhotonEvent into histograms +// +Int_t MHPhotonEvent::Fill(const MParContainer *par, const Stat_t weight) +{ + const MPhotonEvent *evt = dynamic_cast(par); + if (!evt) + { + *fLog << err << dbginf << "No MPhotonEvent found..." << endl; + return kERROR; + } + + // Check if we want to use this class as a single event display + if (fPermanentReset && evt->GetNumPhotons()>0) + { + fHistXY.Reset(); + fHistUV.Reset(); + fHistT.Reset(); + } + + // Get number of photons + const Int_t num = evt->GetNumPhotons(); + + // Set minimum to maximum possible value + Double_t min = FLT_MAX; + + // Loop over all photons and determine the time of the first photon + // FIXME: Should we get it from some statistics container? + for (Int_t idx=0; idxGetDir3(); + //dir *= -1./dir.Z(); + + fHistXY.Fill(x, y); + fHistUV.Fill(u, v); + fHistT.Fill(x, y, t); + fHistWL.Fill(w); + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +/* +Bool_t MHPhotonEvent::Finalize() +{ + const TAxis &axey = *fHistRad.GetYaxis(); + for (int y=1; y<=fHistRad.GetNbinsY(); y++) + { + const Double_t lo = axey.GetBinLowEdge(y); + const Double_t hi = axey.GetBinLowEdge(y+1); + + const Double_t A = (hi*hi-lo*lo)*TMath::Pi()*TMath::Pi(); + + for (int x=1; x<=fHistRad.GetNbinsX(); x++) + fHistRad.SetBinContent(x, y, fHistRad.GetBinContent(x, y)/A); + } + return kTRUE; +} +*/ + +// -------------------------------------------------------------------------- +// +// Make sure that the TProfile2D doesn't fix it's displayed minimum at 0. +// Set the pretty-palette. +// +void MHPhotonEvent::Paint(Option_t *opt) +{ + SetPalette("pretty"); + + fHistT.SetMinimum(); + fHistT.SetMinimum(fHistT.GetMinimum(0)); +} + +// -------------------------------------------------------------------------- +// +void MHPhotonEvent::Draw(Option_t *o) +{ + TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); + pad->SetBorderMode(0); + + AppendPad(); + + pad->Divide(3,2); + + pad->cd(1); + gPad->SetBorderMode(0); + gPad->SetGrid(); + fHistXY.Draw("colz"); + + pad->cd(2); + gPad->SetBorderMode(0); + gPad->SetGrid(); + fHistT.Draw("colz"); + + pad->cd(4); + gPad->SetBorderMode(0); + gPad->SetGrid(); + fHistUV.Draw("colz"); + + pad->cd(5); + gPad->SetBorderMode(0); + gPad->SetGrid(); + fHistWL.Draw(); +} Index: /trunk/MagicSoft/Mars/msim/MHPhotonEvent.h =================================================================== --- /trunk/MagicSoft/Mars/msim/MHPhotonEvent.h (revision 9232) +++ /trunk/MagicSoft/Mars/msim/MHPhotonEvent.h (revision 9232) @@ -0,0 +1,45 @@ +#ifndef MARS_MHPhotonEvent +#define MARS_MHPhotonEvent + +#ifndef MARS_MH +#include "MH.h" +#endif + +#ifndef ROOT_TH2 +#include +#endif + +#ifndef ROOT_TProfile2D +#include +#endif + +class MPhotonEvent; + +class MHPhotonEvent : public MH +{ +private: + TH2D fHistXY; + TH2D fHistUV; + TProfile2D fHistT; + TH1D fHistWL; + + Int_t fType; + Bool_t fPermanentReset; + + void Init(const char *name, const char *title); + + Bool_t SetupFill(const MParList *pList); + Int_t Fill(const MParContainer *par, const Stat_t weight=1); + //Bool_t Finalize(); + +public: + MHPhotonEvent(Double_t max=0, const char *name=0, const char *title=0); + MHPhotonEvent(Int_t type, const char *name=0, const char *title=0); + + void Draw(Option_t *o=""); + void Paint(Option_t *o=""); + + ClassDef(MHPhotonEvent, 1) // Histogram to display the information of MPhotonEvents +}; + +#endif Index: /trunk/MagicSoft/Mars/msim/MPhotonData.cc =================================================================== --- /trunk/MagicSoft/Mars/msim/MPhotonData.cc (revision 9232) +++ /trunk/MagicSoft/Mars/msim/MPhotonData.cc (revision 9232) @@ -0,0 +1,261 @@ +/* ======================================================================== *\ +! +! * +! * 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): Qi Zhe, 06/2007 +! +! Copyright: Software Development, 2000-2009 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MPhotonData +// +// Storage container to store Corsika events +// +// Version 1: +// ---------- +// - First implementation +// +///////////////////////////////////////////////////////////////////////////// +#include "MPhotonData.h" + +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MPhotonData); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MPhotonData::MPhotonData(/*const char *name, const char *title*/) + : fPosX(0), fPosY(0), fCosU(0), fCosV(0), fTime(0), fWavelength(0), + fNumPhotons(1), fProductionHeight(0), fPrimary(MMcEvtBasic::kUNDEFINED), + fTag(-1), fWeight(1) +{ + // fName = name ? name : "MPhotonData"; + // fTitle = title ? title : "Corsika Event Data Information"; +} + +/* +MPhotonData::MPhotonData(const MPhotonData &ph) +: fPosX(ph.fPosX), fPosY(ph.fPosY), fCosU(ph.fCosU), fCosV(ph.fCosV), +fTime(ph.fTime), fWavelength(ph.fWavelength), fNumPhotons(ph.fNumPhotons), +fProductionHeight(ph.fProductionHeight), fPrimary(ph.fPrimary), +fTag(ph.fTag), fWeight(ph.fWeight) +{ +} +*/ + +// -------------------------------------------------------------------------- +// +// Copy function. Copy all data members into obj. +// +void MPhotonData::Copy(TObject &obj) const +{ + MPhotonData &d = static_cast(obj); + + d.fNumPhotons = fNumPhotons; + d.fPosX = fPosX; + d.fPosY = fPosY; + d.fCosU = fCosU; + d.fCosV = fCosV; + d.fWavelength = fWavelength; + d.fPrimary = fPrimary; + d.fTime = fTime; + d.fTag = fTag; + d.fWeight = fWeight; + d.fProductionHeight = fProductionHeight; + + TObject::Copy(obj); +} + +// -------------------------------------------------------------------------- +// +// Set the data member according to the 8 floats read from a reflector-file. +// This function MUST reset all data-members, no matter whether these are +// contained in the input stream. +// +Int_t MPhotonData::FillRfl(Float_t f[8]) +{ + // Check coordinate system!!!! + fWavelength = TMath::Nint(f[0]); + fPosX = f[1]; // [cm] + fPosY = f[2]; // [cm] + fCosU = f[3]; // cos to x + fCosV = f[4]; // cos to y + fTime = f[5]; // [ns] + fProductionHeight = f[6]; + + // f[7]: Camera inclination angle + + fPrimary = MMcEvtBasic::kUNDEFINED; + fNumPhotons = 1; + fTag = -1; + fWeight = 1; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Set the data member according to the 7 floats read from a corsika-file. +// This function MUST reset all data-members, no matter whether these are +// contained in the input stream. +// +// Currently we exchange x and y and set y=-y to convert Corsikas coordinate +// system intpo our own. +// +Int_t MPhotonData::FillCorsika(Float_t f[7]) +{ + const UInt_t n = TMath::Nint(f[0]); + if (n==0) + return kCONTINUE; + + // This seems to be special to mmcs + fWavelength = n%1000; + fPrimary = MMcEvtBasic::ParticleId_t(n/100000); + fNumPhotons = (n/1000)%100; // Force this to be 1! + + if (fNumPhotons!=1) + { + // FIXME: Could be done in MPhotonEvent::ReadCorsikaEvent + gLog << err << "ERROR - MPhotonData::FillCorsika: fNummPhotons not 1, but " << fNumPhotons << endl; + gLog << " This is not yet supported." << endl; + return kERROR; + } + + // x=north, y=west + //fPosX = f[1]; // [cm] + //fPosY = f[2]; // [cm] + //fCosU = f[3]; // cos to x + //fCosV = f[4]; // cos to y + // x=east, y=north + fPosX = f[2]; // [cm] + fPosY = -f[1]; // [cm] + + fCosU = -f[4]; // cos to x + fCosV = f[3]; // cos to y + + fTime = f[5]; // [ns] + + fProductionHeight = f[6]; // [cm] + + // Now reset all data members which are not in the stream + fTag = -1; + fWeight = 1; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Read seven floats from the stream and call FillCorsika for them. +// +Int_t MPhotonData::ReadCorsikaEvt(istream &fin) +{ + Float_t f[7]; + fin.read((char*)&f, 7*4); + + const Int_t rc = FillCorsika(f); + + return rc==kTRUE ? !fin.eof() : rc; +} + +// -------------------------------------------------------------------------- +// +// Read eight floats from the stream and call FillRfl for them. +// +Int_t MPhotonData::ReadRflEvt(istream &fin) +{ + Float_t f[8]; + fin.read((char*)&f, 8*4); + + const Int_t rc = FillRfl(f); + + return rc==kTRUE ? !fin.eof() : rc; +} + +// -------------------------------------------------------------------------- +// +// Print contents. The tag and Weight are only printed if they are different +// from the default. +// +void MPhotonData::Print(Option_t *) const +{ + gLog << inf << endl; + gLog << "Num Photons: " << fNumPhotons << " from " << MMcEvtBasic::GetParticleName(fPrimary) << endl; + gLog << "Wavelength: " << fWavelength << "nm" << endl; + gLog << "Pos X/Y Cos U/V: " << fPosX << "/" << fPosY << " " << fCosU << "/" << fCosV << endl; + gLog << "Time/Prod.Height: " << fTime << "ns/" << fProductionHeight << "cm" << endl; + if (fTag>=0) + gLog << "Tag: " << fTag << endl; + if (fWeight!=1) + gLog << "Weight: " << fWeight << endl; +} + +/* +#include +#include + +// -------------------------------------------------------------------------- +// +// Fill radial photon distance scaled by scale into 1D histogram. +// +void MPhotonData::FillRad(TH1 &hist, Float_t scale) const +{ + hist.Fill(TMath::Hypot(fPosX, fPosY)*scale); +} + +// -------------------------------------------------------------------------- +// +// Fill radial photon distance scaled by scale versus x into 2D histogram. +// +void MPhotonData::FillRad(TH2 &hist, Double_t x, Float_t scale) const +{ + hist.Fill(x, TMath::Hypot(fPosX, fPosY)*scale); +} + +// -------------------------------------------------------------------------- +// +// Fill photon position (x,y) scaled by scale into 2D histogram. +// (Note north in the camera plane is -y, and east -x) +// +void MPhotonData::Fill(TH2 &hist, Float_t scale) const +{ + hist.Fill(fPosY*scale, fPosX*scale); +} + +// -------------------------------------------------------------------------- +// +// Fill photon position (x,y) scaled by scale versus z into 3D histogram. +// (Note north in the camera plane is -y, and east -x) +// +void MPhotonData::Fill(TH3 &hist, Double_t z, Float_t scale) const +{ + hist.Fill(fPosY*scale, fPosX*scale); +} + +*/ Index: /trunk/MagicSoft/Mars/msim/MPhotonData.h =================================================================== --- /trunk/MagicSoft/Mars/msim/MPhotonData.h (revision 9232) +++ /trunk/MagicSoft/Mars/msim/MPhotonData.h (revision 9232) @@ -0,0 +1,153 @@ +#ifndef MARS_MPhotonData +#define MARS_MPhotonData + +#ifndef MARS_MMcEvtBasic +#include "MMcEvtBasic.h" +#endif + +#ifndef ROOT_TVector2 +#include +#endif + +#ifndef ROOT_TVector3 +#include +#endif + +#ifndef ROOT_TQuaternion +#include +#endif + +// gcc 3.2 +//class ifstream; +#include +/* +class TH1; +class TH2; +class TH3; +*/ +class MCorsikaRunHeader; + +class MPhotonData : public TObject //MParContainer +{ +private: + Float_t fPosX; // [cm] X (north) at observation level + Float_t fPosY; // [cm] Y (west) at observation level + + Float_t fCosU; // [cos x] U direction cosine to x-axis + Float_t fCosV; // [cos y] V direction cosine to y-axis + + Float_t fTime; // [ns] Time since first interaction or entrance into atmosphere + // 17M + UShort_t fWavelength; // [nm] Wavelength + // 19M + UInt_t fNumPhotons; // Number of cherenkov photons ins bunch + Float_t fProductionHeight; // [cm] Height of bunch production + MMcEvtBasic::ParticleId_t fPrimary; // Type of emitting particle + // 22M + // gzip + // 25M + // raw + // 32M + + Int_t fTag; //! A tag for external use + Float_t fWeight; //! A weight for external use + +protected: + virtual Int_t FillCorsika(Float_t f[7]); + virtual Int_t FillRfl(Float_t f[8]); + +public: + MPhotonData(/*const char *name=NULL, const char *title=NULL*/); + //MPhotonData(const MPhotonData &ph); + + // Getter 1D + Float_t GetPosX() const { return fPosX; } + Float_t GetPosY() const { return fPosY; } + + Float_t GetCosU() const { return fCosU; } + Float_t GetCosV() const { return fCosV; } + Double_t GetCosW() const { return TMath::Sqrt(1 - fCosU*fCosU - fCosV*fCosV); } // cos(Theta) + Double_t GetTheta() const { return TMath::ACos(GetCosW()); } + + Double_t GetTime() const { return fTime; } + + // Getter 2D + TVector2 GetPos2() const { return TVector2(fPosX, fPosY); } + TVector2 GetDir2() const { return TVector2(fCosU, fCosV);} +// TVector2 GetDir2() const { return TVector2(fCosU, fCosV)/(1-TMath::Hypot(fCosU, fCosV)); } + + // Getter 3D + TVector3 GetPos3() const { return TVector3(fPosX, fPosY, 0); } + TVector3 GetDir3() const { return TVector3(fCosU, fCosV, -GetCosW()); } + + // Getter 4D // cm // ns + TQuaternion GetPosQ() const { return TQuaternion(GetPos3(), fTime); } + // FIXME: v in air! // m/s cm? / ns + TQuaternion GetDirQ() const { return TQuaternion(GetDir3(), 1./(TMath::C()*100/1e9)); } + + // Getter Others + UShort_t GetWavelength() const { return fWavelength; } + MMcEvtBasic::ParticleId_t GetPrimary() const { return fPrimary; } + + //virtual Float_t GetWeight() const { return 1; } + virtual Float_t GetWeight() const { return fWeight; } + void SetWeight(Float_t w=1) { fWeight=w; } + + // Setter + void SetPosition(Float_t x, Float_t y) { fPosX=x; fPosY=y; } + void SetDirection(Float_t u, Float_t v) { fCosU=u; fCosV=v; } + + void SetPosition(const TVector2 &p) { fPosX=p.X(); fPosY=p.Y(); } + void SetDirection(const TVector2 &d) { fCosU=d.X(); fCosV=d.Y(); } + + void SetPosition(const TQuaternion &p) { fPosX=p.fVectorPart.X(); fPosY=p.fVectorPart.Y(); fTime=p.fRealPart; } + void SetDirection(const TQuaternion &d) { fCosU=d.fVectorPart.X(); fCosV=d.fVectorPart.Y(); } + + void SetPrimary(MMcEvtBasic::ParticleId_t p) { fPrimary=p; } + void SetWavelength(UShort_t wl) { fWavelength=wl; } + + void AddTime(Double_t dt) { fTime += dt; } + void SetTime(Double_t t) { fTime = t; } + + void Copy(TObject &obj) const; + + void SetTag(Int_t tag) { fTag=tag; } + Int_t GetTag() const { return fTag; } + + //void Clear(Option_t * = NULL); + void Print(Option_t * = NULL) const; + //void Draw (Option_t * = NULL); + Bool_t IsSortable() const { return kTRUE; } + Int_t Compare(const TObject *obj) const + { + const MPhotonData &d = *static_cast(obj); + if (fTimed.fTime) + return 1; + return 0; + } +/* + void FillRad(TH1 &hist, Float_t scale=1) const; + void FillRad(TH2 &hist, Double_t x, Float_t scale=1) const; + void Fill(TH2 &hist, Float_t scale=1) const; + void Fill(TH3 &hist, Double_t z, Float_t scale=1) const; + */ + Int_t ReadCorsikaEvt(istream &fin); + Int_t ReadRflEvt(istream &fin); + + ClassDef(MPhotonData, 1) //Container to store a cherenkov photon bunch from a CORSUKA file +}; +/* +class MPhotonDataWeighted : public MPhotonData +{ +private: + Float_t fWeight; // A weight for external use + +public: + Float_t GetWeight() const { return fWeight; } + + ClassDef(MPhotonData, 1) +}; +*/ +#endif Index: /trunk/MagicSoft/Mars/msim/MPhotonEvent.cc =================================================================== --- /trunk/MagicSoft/Mars/msim/MPhotonEvent.cc (revision 9232) +++ /trunk/MagicSoft/Mars/msim/MPhotonEvent.cc (revision 9232) @@ -0,0 +1,424 @@ +/* ======================================================================== *\ +! +! * +! * 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): Qi Zhe, 06/2007 +! +! Copyright: Software Development, 2000-2009 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MPhotonEvent +// +// Storage container to store photon collections +// +// The class is designed to be extremely fast which is important taking into +// account the extremely high number of photons. This has some impacts on +// its handling. +// +// The list has to be kept consistent, i.e. without holes. +// +// There are two ways to achieve this: +// +// a) Use RemoveAt to remove an entry somewhere +// b) Compress() the TClonesArray afterwards +// +// Compress is not the fastes, so there is an easier way. +// +// a) When you loop over the list and want to remove an entry copy all +// following entry backward in the list, so that the hole will +// be created at its end. +// b) Call Shrink(n) with n the number of valid entries in your list. +// +// To loop over the TClonesArray you can use a TIter which has some +// unnecessary overhead and therefore is slower than necessary. +// +// Since the list is kept consistent you can use a simple loop saving +// a lot of CPU time taking into account the high number of calls to +// TObjArrayIter::Next which you would create. +// +// Here is an example (how to remove every second entry) +// +// Int_t cnt = 0; +// +// const Int_t num = event->GetNumPhotons(); +// for (Int_t idx=0; idxShrink(cnt); +// +// ---------------------------------- or ------------------------------- +// +// TClonesArray &arr = MPhotonEvent->GetArray(); +// +// Int_t cnt = 0; +// +// const Int_t num = arr.GetEntriesFast(); +// for (Int_t idx=0; idx(arr.UncheckedAt(idx)); +// +// *static_cast(arr.UncheckedAt(cnt++)) = *dat; +// } +// +// MPhotonEvent->Shrink(cnt); +// +// +// Version 1: +// ---------- +// - First implementation +// +///////////////////////////////////////////////////////////////////////////// +#include "MPhotonEvent.h" + +#include +#include + +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MPhotonData.h" + +ClassImp(MPhotonEvent); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default constructor. It initializes all arrays with zero size. +// +MPhotonEvent::MPhotonEvent(const char *name, const char *title) + : fData("MPhotonData", 1) +{ + fName = name ? name : "MPhotonEvent"; + fTitle = title ? title : "Corsika Event Data Information"; + + fData.SetBit(TClonesArray::kForgetBits); + fData.BypassStreamer(kFALSE); +} + +const char *MPhotonEvent::GetClassName() const +{ + return static_cast(fData.GetClass())->GetName(); +} + +// -------------------------------------------------------------------------- +// +// If n is smaller than the current allocated array size a reference to +// the n-th entry is returned, otherwise an entry at n is created +// calling the default constructor. Note, that there is no range check +// but it is not recommended to call this function with +// n>fData.GetSize() +// +MPhotonData &MPhotonEvent::Add(Int_t n) +{ + // Do not modify this. It is optimized for execution + // speed and flexibility! + TObject *o = n(*o); +} + +// -------------------------------------------------------------------------- +// +// Add a new photon (MPhtonData) at the end of the array. +// In this case the default constructor of MPhotonData is called. +// +// A reference to the new object is returned. +// +MPhotonData &MPhotonEvent::Add() +{ + return Add(GetNumPhotons()); +} + +// -------------------------------------------------------------------------- +// +// Get the i-th photon from the array. Not, for speed reasons there is no +// range check so you are responsible that you do not excess the number +// of photons (GetNumPhotons) +// +MPhotonData &MPhotonEvent::operator[](UInt_t idx) +{ + return *static_cast(fData.UncheckedAt(idx)); +} + +// -------------------------------------------------------------------------- +// +// Get the i-th photon from the array. Not, for speed reasons there is no +// range check so you are responsible that you do not excess the number +// of photons (GetNumPhotons) +// +const MPhotonData &MPhotonEvent::operator[](UInt_t idx) const +{ + return *static_cast(fData.UncheckedAt(idx)); +} + +Int_t MPhotonEvent::ReadCorsikaEvt(istream &fin) +{ + Int_t n = 0; + + while (1) + { + // Check the first four bytes + char c[4]; + fin.read(c, 4); + + // End of stream + if (!fin) + return kFALSE; + + // Check if we found the end of the event + if (!memcmp(c, "EVTE", 4)) + break; + + // The first for byte contained data already --> go back + fin.seekg(-4, ios::cur); + + // Do not modify this. It is optimized for execution + // speed and flexibility! + MPhotonData &ph = Add(n); + // It checks how many entries the lookup table has. If it has enough + // entries and the entry was already allocated, we can re-use it, + // otherwise we have to allocate it. + + // Now we read a single cherenkov bunch. Note that for speed reason we have not + // called the constructor if the event was already constructed (virtual table + // set), consequently we must make sure that ReadCorsikaEvent does reset + // all data mebers no matter whether they are read or not. + const Int_t rc = ph.ReadCorsikaEvt(fin); + + // Evaluate result from reading event + switch (rc) + { + case kCONTINUE: continue; // No data in this bunch... skip it. + case kFALSE: return kFALSE; // End of stream + case kERROR: return kERROR; // Error occured + } + + // FIXME: If fNumPhotons!=1 add the photon more than once + + // Now increase the number of entries which are kept, + // i.e. keep this photon(s) + n++; + } + + Shrink(n); + + SetReadyToSave(); + + //*fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl; + return kTRUE; +} + +Int_t MPhotonEvent::ReadRflEvt(std::istream &fin) +{ + Int_t n = 0; + + while (1) + { + // Check the first four bytes + char c[13]; + fin.read(c, 13); + + // End of stream + if (!fin) + return kFALSE; + + // Check if we found the end of the event + if (!memcmp(c, "\nEND---EVENT\n", 13)) + break; + + // The first for byte contained data already --> go back + fin.seekg(-13, ios::cur); + + // Do not modify this. It is optimized for execution + // speed and flexibility! + //TObject *o = n(o)->ReadRflEvt(fin); + const Int_t rc = Add(n).ReadRflEvt(fin); + + // Evaluate result from reading event + switch (rc) + { + case kCONTINUE: continue; // No data in this bunch... skip it. + case kFALSE: return kFALSE; // End of stream + case kERROR: return kERROR; // Error occured + } + + // Now increase the number of entries which are kept, + // i.e. keep this photon(s) + n++; + } + + Shrink(n); + + SetReadyToSave(); + + //*fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl; + return kTRUE; +} + +void MPhotonEvent::Print(Option_t *) const +{ + fData.Print(); +} + +/* +// -------------------------------------------------------------------------- +// +// Fills all photon distances scaled with scale (default=1) into a TH1 +// +void MPhotonEvent::FillRad(TH1 &hist, Float_t scale) const +{ + MPhotonData *ph=NULL; + + TIter Next(&fData); + while ((ph=(MPhotonData*)Next())) + ph->FillRad(hist, scale); +} + +// -------------------------------------------------------------------------- +// +// Fills all photon distances scaled with scale (default=1) versus x +// into a TH2 +// +void MPhotonEvent::FillRad(TH2 &hist, Double_t x, Float_t scale) const +{ + MPhotonData *ph=NULL; + + TIter Next(&fData); + while ((ph=(MPhotonData*)Next())) + ph->FillRad(hist, x, scale); +} + +// -------------------------------------------------------------------------- +// +// Fills all photons (x,y) scaled with scale (default=1) into a TH2. +// +void MPhotonEvent::Fill(TH2 &hist, Float_t scale) const +{ + MPhotonData *ph=NULL; + + TIter Next(&fData); + while ((ph=(MPhotonData*)Next())) + ph->Fill(hist, scale); +} + +// -------------------------------------------------------------------------- +// +// Fills all photons (x,y) scaled with scale (default=1) versus z into a TH3 +// +void MPhotonEvent::Fill(TH3 &hist, Double_t z, Float_t scale) const +{ + MPhotonData *ph=NULL; + + TIter Next(&fData); + while ((ph=(MPhotonData*)Next())) + ph->Fill(hist, z, scale); +} + +void MPhotonEvent::DrawPixelContent(Int_t num) const +{ +} + +#include "MHexagon.h" +#include "MGeomCam.h" +#include + +// ------------------------------------------------------------------------ +// +// Fill a reflector event. Sums all pixels in each pixel as the +// pixel contents. +// +// WARNING: Due to the estimation in DistanceToPrimitive and the +// calculation in pixels instead of x, y this is only a +// rough estimate. +// +Bool_t MPhotonEvent::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const +{ + // + // sum the photons content in each pixel + // + val = 0; + + MHexagon hex(cam[idx]); + + MPhotonData *ph=NULL; + + TIter Next(&fData); + + switch (type) + { + + case 1: // time + while ((ph=(MPhotonData*)Next())) + if (hex.DistanceToPrimitive(ph->GetPosY(), ph->GetPosX())<=0) + val += ph->GetTime(); + return val/fData.GetEntriesFast()>0; + + default: // signal + while ((ph=(MPhotonData*)Next())) + if (hex.DistanceToPrimitive(ph->GetPosY()*10, ph->GetPosX()*10)<=0) + val += cam.GetPixRatio(idx); + return val>0; + } + + return kFALSE; + +} +*/ + +// ------------------------------------------------------------------------ +// +// You can call Draw() to add the photons to the current pad. +// The photons are painted each tim ethe pad is updated. +// Make sure that you use the right (world) coordinate system, +// like created, eg. by the MHCamera histogram. +// +void MPhotonEvent::Paint(Option_t *) +{ + MPhotonData *ph=NULL; + + TMarker m; + m.SetMarkerStyle(kFullDotMedium); // Gtypes.h + + TIter Next(&fData); + while ((ph=(MPhotonData*)Next())) + { + m.SetX(ph->GetPosY()*10); // north + m.SetY(ph->GetPosX()*10); // east + m.Paint(); + } +} Index: /trunk/MagicSoft/Mars/msim/MPhotonEvent.h =================================================================== --- /trunk/MagicSoft/Mars/msim/MPhotonEvent.h (revision 9232) +++ /trunk/MagicSoft/Mars/msim/MPhotonEvent.h (revision 9232) @@ -0,0 +1,173 @@ +#ifndef MARS_MPhotonEvent +#define MARS_MPhotonEvent + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif +/* +#ifndef MARS_MCamEvent +#include "MCamEvent.h" +#endif +*/ +#ifndef ROOT_TClonesArray +#include +#endif + +// gcc 3.2 +//class ifstream; +#include + +/* +class TH1; +class TH2; +class TH3; +*/ + +class MPhotonData; +class MCorsikaRunHeader; + +class MyClonesArray : public TClonesArray +{ +public: + TObject **FirstRef() { return fCont; } + + void FastRemove(Int_t idx1, Int_t idx2) + { + // Remove object at index idx. + + //if (!BoundsOk("RemoveAt", idx1)) return 0; + //if (!BoundsOk("RemoveAt", idx2)) return 0; + + Long_t dtoronly = TObject::GetDtorOnly(); + + idx1 -= fLowerBound; + idx2 -= fLowerBound; + + for (TObject **obj=fCont+idx1; obj<=fCont+idx2; obj++) + { + if (!*obj) + continue; + + if ((*obj)->TestBit(kNotDeleted)) { + // Tell custom operator delete() not to delete space when + // object fCont[i] is deleted. Only destructors are called + // for this object. + TObject::SetDtorOnly(*obj); + delete *obj; + } + + *obj = 0; + // recalculate array size + } + TObject::SetDtorOnly((void*)dtoronly); + + if (idx1<=fLast && fLast<=idx2) + { + do { + fLast--; + } while (fLast >= 0 && fCont[fLast] == 0); + } + + Changed(); + } +}; + +class MPhotonEvent : public MParContainer//, public MCamEvent +{ +private: + TClonesArray fData; + +public: + MPhotonEvent(const char *name=NULL, const char *title=NULL); + + //void Clear(Option_t * = NULL); + void Print(Option_t * = NULL) const; + + const char *GetClassName() const; + + TClonesArray &GetArray() { return fData; } + const TClonesArray &GetArray() const { return fData; } + void Remove(TObject *o) { fData.Remove(o); } + + MPhotonData &Add(Int_t n); + MPhotonData &Add(); + + MPhotonData &operator[](UInt_t idx); + const MPhotonData &operator[](UInt_t idx) const; + + Int_t Shrink(Int_t n) + { + // The number of object written by the streamer is defined by + // GetEntriesFast(), i.e. this number must be shrinked to + // the real array size. We use RemoveAt instead of ExpandCreate + // because RemoveAt doesn't free memory. Thus in the next + // iteration it can be reused and doesn't need to be reallocated. + // Do not change this. It is optimized for execution speed + // for (int i=fData.GetSize()-1; i>=n; i--) + // fData.RemoveAt(i); + static_cast(fData).FastRemove(n, fData.GetSize()-1); + + return fData.GetEntriesFast(); + } +/* + void Expand(Int_t n) + { + for (int i=fData.GetSize(); i=fTimeFirst; } + + ClassDef(MPhotonStatistics, 1) +}; +#endif Index: /trunk/MagicSoft/Mars/msim/MSimAbsorption.cc =================================================================== --- /trunk/MagicSoft/Mars/msim/MSimAbsorption.cc (revision 9232) +++ /trunk/MagicSoft/Mars/msim/MSimAbsorption.cc (revision 9232) @@ -0,0 +1,335 @@ +/* ======================================================================== *\ +! +! * +! * 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 1/2009 +! +! Copyright: Software Development, 2000-2009 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MSimAbsorption +// +// Task to calculate wavelength or incident angle dependent absorption +// +// Input Containers: +// MPhotonEvent +// MCorsikaEvtHeader +// MCorsikaRunHeader +// +// Output Containers: +// MPhotonEvent +// +////////////////////////////////////////////////////////////////////////////// +#include "MSimAbsorption.h" + +#include + +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" + +#include "MBinning.h" +#include "MArrayD.h" + +#include "MSpline3.h" + +#include "MCorsikaEvtHeader.h" +#include "MCorsikaRunHeader.h" +#include "MPhotonEvent.h" +#include "MPhotonData.h" + +ClassImp(MSimAbsorption); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +MSimAbsorption::MSimAbsorption(const char* name, const char *title) + : fEvt(0), fHeader(0), fSpline(0), fUseTheta(kFALSE) +{ + fName = name ? name : "MSimAbsorption"; + fTitle = title ? title : "Task to calculate wavelength dependent absorption"; +} + +// -------------------------------------------------------------------------- +// +// Calls Clear() +// +MSimAbsorption::~MSimAbsorption() +{ + Clear(); +} + +// -------------------------------------------------------------------------- +// +// Delete fSpline and set to 0. +// +void MSimAbsorption::Clear(Option_t *) +{ + if (fSpline) + delete fSpline; + fSpline=0; +} + +// -------------------------------------------------------------------------- +// +// Read a TGraph from a file and return a newly allocated MSpline3. +// +MSpline3 *MSimAbsorption::ReadSpline(const char *fname) +{ + *fLog << inf << "Reading spline from " << fname << endl; + + TGraph g(fname); + if (g.GetN()==0) + { + *fLog << err << "ERROR - No data points from " << fname << "." << endl; + return 0; + } + + // option: b1/e1 b2/e2 (first second derivative?) + // option: valbeg/valend (first second derivative?) + + return new MSpline3(g); +} + +// -------------------------------------------------------------------------- +// +// Initializes a spline from min to max with n points with 1. +// +void MSimAbsorption::InitUnity(UInt_t n, Float_t min, Float_t max) +{ + // Delete the existing spline + Clear(); + + // We need at least two points (the edges) + if (n<2) + return; + + // Initialize an array with the x-values + const MBinning bins(n-1, min, max); + + // Initialize an array with all one + MArrayD y(n); + y.Reset(1); + + // Set the spline + fSpline = new MSpline3(bins.GetEdges(), y.GetArray(), n); +} + +// -------------------------------------------------------------------------- +// +// Takes the existing fSpline and multiplies it with the given spline. +// As x-points the values from fSpline are used. +// +void MSimAbsorption::Multiply(const TSpline3 &spline) +{ + if (!fSpline) + { + Clear(); + // WARNING WARNING WARNING: This is a very dangerous cast! + fSpline = (MSpline3*)spline.Clone(); + return; + } + + // Initialize a TGraph with the number of knots from a TSpline + TGraph g(fSpline->GetNp()); + + // Loop over all knots + for (int i=0; iGetNp(); i++) + { + // Get th i-th knot + Double_t x, y; + fSpline->GetKnot(i, x, y); + + // Multiply y by the value from the given spline + y *= spline.Eval(x); + + // push the point "back" + g.SetPoint(i, x, y); + } + + // Clear the spline and initialize a new one from the updated TGraph + Clear(); + fSpline = new MSpline3(g); +} + +// -------------------------------------------------------------------------- +// +// Initializes a TSpline3 with n knots x and y. Call Multiply for it. +// +void MSimAbsorption::Multiply(UInt_t n, const Double_t *x, const Double_t *y) +{ + const TSpline3 spline("Spline", const_cast(x), const_cast(y), n); + Multiply(spline); +} + +// -------------------------------------------------------------------------- +// +// Read a Spline from a file using ReadSpline. +// Call Multiply for it. +// +void MSimAbsorption::Multiply(const char *fname) +{ + const MSpline3 *spline = ReadSpline(fname); + if (!spline) + return; + + fFileName = ""; + + Multiply(*spline); + + delete spline; +} + +// -------------------------------------------------------------------------- +// +// Read a spline from a file and set fSpline accfordingly. +// +Bool_t MSimAbsorption::ReadFile(const char *fname) +{ + if (fname) + fFileName = fname; + + MSpline3 *spline = ReadSpline(fFileName); + if (!spline) + return kFALSE; + + + // option: b1/e1 b2/e2 (first second derivative?) + // option: valbeg/valend (first second derivative?) + + Clear(); + fSpline = spline; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Search for the needed parameter containers. Read spline from file +// calling ReadFile(); +// +Int_t MSimAbsorption::PreProcess(MParList *pList) +{ + fEvt = (MPhotonEvent*)pList->FindObject("MPhotonEvent"); + if (!fEvt) + { + *fLog << err << "MPhotonEvent not found... aborting." << endl; + return kFALSE; + } + + fHeader = (MCorsikaEvtHeader*)pList->FindObject("MCorsikaEvtHeader"); + if (!fHeader) + { + *fLog << err << "MCorsikaEvtHeader not found... aborting." << endl; + return kFALSE; + } + + return ReadFile(); +} + +// -------------------------------------------------------------------------- +// +Bool_t MSimAbsorption::ReInit(MParList *pList) +{ + MCorsikaRunHeader *h = (MCorsikaRunHeader*)pList->FindObject("MCorsikaRunHeader"); + if (!h) + { + *fLog << err << "MCorsikaRunHeader not found... aborting." << endl; + return kFALSE; + } + + if (h->GetWavelengthMin()GetXmin()) + *fLog << warn << "WARNING - Lower bound of wavelength bandwidth exceeds lower bound of spline." << endl; + + if (h->GetWavelengthMax()>fSpline->GetXmax()) + *fLog << warn << "WARNING - Upper bound of wavelength bandwidth exceeds upper bound of spline." << endl; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Throw all events out of the MPhotonEvent which don't survive. +// Depending on fUseTheta either the wavelength or the incident angle +// is used. +// +Int_t MSimAbsorption::Process() +{ + // Get the number of photons in the list + const Int_t num = fEvt->GetNumPhotons(); + + // Counter for number of total and final events + Int_t cnt = 0; + for (Int_t i=0; iEval(wl); + + // Get a random value between 0 and 1 to determine whether the photn will survive + // gRandom->Rndm() = [0;1[ + if (gRandom->Rndm()>=eff) + continue; + + // Copy the surviving events bakc in the list + (*fEvt)[cnt++] = ph; + } + + // Now we shrink the array to the number of new entries. + fEvt->Shrink(cnt); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// FileName: reflectivity.txt +// UseTheta: No +// +Int_t MSimAbsorption::ReadEnv(const TEnv &env, TString prefix, Bool_t print) +{ + Bool_t rc = kFALSE; + if (IsEnvDefined(env, prefix, "FileName", print)) + { + rc = kTRUE; + SetFileName(GetEnvValue(env, prefix, "FileName", fFileName)); + } + + if (IsEnvDefined(env, prefix, "UseTheta", print)) + { + rc = kTRUE; + SetUseTheta(GetEnvValue(env, prefix, "UseTheta", fUseTheta)); + } + + return rc; +} Index: /trunk/MagicSoft/Mars/msim/MSimAbsorption.h =================================================================== --- /trunk/MagicSoft/Mars/msim/MSimAbsorption.h (revision 9232) +++ /trunk/MagicSoft/Mars/msim/MSimAbsorption.h (revision 9232) @@ -0,0 +1,59 @@ +#ifndef MARS_MSimAbsorption +#define MARS_MSimAbsorption + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +class TSpline3; + +class MParList; +class MSpline3; +class MPhotonEvent; +class MCorsikaEvtHeader; + +class MSimAbsorption : public MTask +{ +private: + MPhotonEvent *fEvt; //! Event stroing the photons + MCorsikaEvtHeader *fHeader; //! Header storing event information + + MSpline3 *fSpline; //! Spline to interpolate wavelength or incident angle + + TString fFileName; // FileName of the file containing the curve + Bool_t fUseTheta; // Switches between using wavelength or incident angle + + // MParContainer + Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print=kFALSE); + + // MTask + Int_t PreProcess(MParList *pList); + Bool_t ReInit(MParList *pList); + Int_t Process(); + + // MSimAbsorption + MSpline3 *ReadSpline(const char *fname); + +public: + MSimAbsorption(const char *name=NULL, const char *title=NULL); + ~MSimAbsorption(); + + // MSimAbsorption + Bool_t ReadFile(const char *fname=0); + void SetFileName(const char *fname) { fFileName=fname; } + + void SetUseTheta(Bool_t b=kTRUE) { fUseTheta = b; } + + void InitUnity(UInt_t n, Float_t min, Float_t max); + + void Multiply(const char *fname); + void Multiply(const TSpline3 &spline); + void Multiply(UInt_t n, const Double_t *x, const Double_t *y); + + // TObject + void Clear(Option_t *o=""); + + ClassDef(MSimAbsorption, 0) // Task to calculate wavelength or incident angle dependent absorption +}; + +#endif Index: /trunk/MagicSoft/Mars/msim/MSimPointingPos.cc =================================================================== --- /trunk/MagicSoft/Mars/msim/MSimPointingPos.cc (revision 9232) +++ /trunk/MagicSoft/Mars/msim/MSimPointingPos.cc (revision 9232) @@ -0,0 +1,135 @@ +/* ======================================================================== *\ +! +! * +! * 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 1/2009 +! +! Copyright: MAGIC Software Development, 2000-2009 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MSimPointingPos +// +// This task is meant to simulate the pointing position (mirror orientation). +// This depends on the direction from which the shower is coming but also +// on the user request (e.g. Wobble mode). +// +// WARNING: Currently the telescope orientation is just fixed to the +// direction in the run-header (if a view cone was given) or +// the direction in the evt-header (if no view cone given) +// +// Input Containers: +// MCorsikaRunHeader +// MCorsikaEvtHeader +// +// Output Containers: +// MPointingPos +// +////////////////////////////////////////////////////////////////////////////// +#include "MSimPointingPos.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" + +#include "MCorsikaEvtHeader.h" +#include "MCorsikaRunHeader.h" + +#include "MPointingPos.h" + +ClassImp(MSimPointingPos); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +MSimPointingPos::MSimPointingPos(const char* name, const char *title) + : fRunHeader(0), fEvtHeader(0), fPointing(0) +{ + fName = name ? name : "MSimPointingPos"; + fTitle = title ? title : "Task to simulate the pointing position (mirror orientation)"; +} + + +// -------------------------------------------------------------------------- +// +// Search for all necessary containers +// +Int_t MSimPointingPos::PreProcess(MParList *pList) +{ + fPointing = (MPointingPos*)pList->FindCreateObj("MPointingPos"); + if (!fPointing) + return kFALSE; + + fRunHeader = (MCorsikaRunHeader*)pList->FindObject("MCorsikaRunHeader"); + if (!fRunHeader) + { + *fLog << err << "MCorsikaRunHeader not found... aborting." << endl; + return kFALSE; + } + + fEvtHeader = (MCorsikaEvtHeader*)pList->FindObject("MCorsikaEvtHeader"); + if (!fEvtHeader) + { + *fLog << err << "MCorsikaEvtHeader not found... aborting." << endl; + return kFALSE; + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +Int_t MSimPointingPos::Process() +{ + // If a view cone is given use the fixed telescope orientation + const Bool_t fixed = fRunHeader->HasViewCone(); + + // Local sky coordinates (direction of telescope axis) + const Double_t zd = fixed ? fRunHeader->GetZdMin() : fEvtHeader->GetZd()*TMath::RadToDeg(); // x==north + const Double_t az = fixed ? fRunHeader->GetAzMin() : fEvtHeader->GetAz()*TMath::RadToDeg(); // y==west + + // Setup the pointing position + fPointing->SetLocalPosition(zd, az); + + return kTRUE; +} + +/* +Int_t MSimPointingPos::ReadEnv(const TEnv &env, TString prefix, Bool_t print) +{ + Bool_t rc = kFALSE; + if (IsEnvDefined(env, prefix, "FileName", print)) + { + rc = kTRUE; + SetFileName(GetEnvValue(env, prefix, "FileName", fFileName)); + } + + if (IsEnvDefined(env, prefix, "UseTheta", print)) + { + rc = kTRUE; + SetUseTheta(GetEnvValue(env, prefix, "UseTheta", fUseTheta)); + } + + return rc; +} +*/ Index: /trunk/MagicSoft/Mars/msim/MSimPointingPos.h =================================================================== --- /trunk/MagicSoft/Mars/msim/MSimPointingPos.h (revision 9232) +++ /trunk/MagicSoft/Mars/msim/MSimPointingPos.h (revision 9232) @@ -0,0 +1,39 @@ +#ifndef MARS_MSimPointingPos +#define MARS_MSimPointingPos + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +class MParList; +class MCorsikaEvtHeader; +class MCorsikaRunHeader; +class MPointingPos; + +class MSimPointingPos : public MTask +{ +private: + MCorsikaRunHeader *fRunHeader; //! Header storing event information + MCorsikaEvtHeader *fEvtHeader; //! Header storing event information + MPointingPos *fPointing; //! Output storing telescope poiting position + + // MParContainer + //Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print=kFALSE); + + // MTask + Int_t PreProcess(MParList *pList); + Int_t Process(); + + +public: + MSimPointingPos(const char *name=NULL, const char *title=NULL); + + // MSimPointingPos + + // TObject + + ClassDef(MSimPointingPos, 0) // Task to simulate the pointing position (mirror orientation) +}; + +#endif + Index: /trunk/MagicSoft/Mars/msim/Makefile =================================================================== --- /trunk/MagicSoft/Mars/msim/Makefile (revision 9232) +++ /trunk/MagicSoft/Mars/msim/Makefile (revision 9232) @@ -0,0 +1,36 @@ +################################################################## +# +# makefile +# +# for the MARS software +# +################################################################## +include ../Makefile.conf.$(OSTYPE) +include ../Makefile.conf.general + +#------------------------------------------------------------------------------ + +# +# Handling name of the Root Dictionary Files +# +CINT = Sim + +# +# connect the include files defined in the config.mk file +# +INCLUDES = -I. -I../mbase -I../mmc -I../mgeom -I../mgui -I../mcorsika \ + -I../mpointing -I../mreflector -I../mhbase + +SRCFILES = MPhotonData.cc \ + MPhotonEvent.cc \ + MHPhotonEvent.cc \ + MSimAbsorption.cc \ + MSimPointingPos.cc + +############################################################ + +all: $(OBJS) + +include ../Makefile.rules + +mrproper: clean rmbak Index: /trunk/MagicSoft/Mars/msim/SimIncl.h =================================================================== --- /trunk/MagicSoft/Mars/msim/SimIncl.h (revision 9232) +++ /trunk/MagicSoft/Mars/msim/SimIncl.h (revision 9232) @@ -0,0 +1,3 @@ +#ifndef __CINT__ + +#endif // __CINT__ Index: /trunk/MagicSoft/Mars/msim/SimLinkDef.h =================================================================== --- /trunk/MagicSoft/Mars/msim/SimLinkDef.h (revision 9232) +++ /trunk/MagicSoft/Mars/msim/SimLinkDef.h (revision 9232) @@ -0,0 +1,17 @@ +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class MPhotonEvent+; +#pragma link C++ class MPhotonData+; + +#pragma link C++ class MHPhotonEvent+; + +#pragma link C++ class MPhotonStatistics+; + +#pragma link C++ class MSimPointingPos+; +#pragma link C++ class MSimAbsorption+; + +#endif