Index: /trunk/MagicSoft/Mars/mtemp/MApplyPadding.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/MApplyPadding.h (revision 1678) +++ /trunk/MagicSoft/Mars/mtemp/MApplyPadding.h (revision 1678) @@ -0,0 +1,55 @@ +#ifndef MARS_MApplyPadding +#define MARS_MApplyPadding + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#include "TRandom3.h" +#include "TH1.h" +#include "TH2.h" + +class MGeomCam; +class MCerPhotEvt; +class MPedestalCam; +class MMcEvt; +class MPedestalCam; +class MSigmabar; +class MParList; + +class MApplyPadding : public MTask +{ +private: + const MGeomCam *fCam; + MCerPhotEvt *fEvt; + MSigmabar *fSigmabar; + TRandom3 *fRnd; + Int_t fRunType; + Int_t fGroup; + char *fDatabaseFilename; + TH1D *fHSigmabarMax; + MMcEvt *fMcEvt; + MPedestalCam *fPed; + TH2D *fTest; + Bool_t fUseHistogram; + Double_t fFixedSigmabar; + +public: + MApplyPadding(const char *name=NULL, const char *title=NULL); + ~MApplyPadding(); + + Bool_t PreProcess(MParList *pList); + Bool_t Process(); + Bool_t PostProcess(); + + void SetRunType(Int_t runtype) { fRunType = runtype; } + void SetGroup(Int_t group) { fGroup = group; } + void SetDatabaseFile(char *filename) { fDatabaseFilename = filename; } + void SetTargetLevel(Double_t sigmabar) { fFixedSigmabar = sigmabar; fUseHistogram=kFALSE; } + Bool_t SetDefiningHistogram(TH1D *histo); + + ClassDef(MApplyPadding, 1) // task for applying padding +}; + +#endif + Index: /trunk/MagicSoft/Mars/mtemp/MHSigmaPixel.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/MHSigmaPixel.cc (revision 1678) +++ /trunk/MagicSoft/Mars/mtemp/MHSigmaPixel.cc (revision 1678) @@ -0,0 +1,149 @@ +/* ======================================================================== *\ +! +! * +! * 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): Robert Wagner 10/2002 +! Copyright: MAGIC Software Development, 2000-2002 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// // +// MHSigmaPixel // +// // +// 2D-Histogram pedestal sigma vs pixel number // +// // +////////////////////////////////////////////////////////////////////////////// + +#include "MHSigmaPixel.h" + +#include + +#include + +#include "MPedestalCam.h" +#include "MPedestalPix.h" +#include "MSigmabar.h" + +#include "MBinning.h" +#include "MParList.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MHSigmaPixel); + +// -------------------------------------------------------------------------- +// +// Default Constructor. It sets name and title of the histogram. +// +MHSigmaPixel::MHSigmaPixel(const char *name, const char *title) +{ + // + // set the name and title of this object + // + fName = name ? name : "MHSigmaPixel"; + fTitle = title ? title : "2-D histogram in sigma and pixel"; + + fHist.SetDirectory(NULL); + + fHist.SetTitle("Sigma vs pixel #"); + fHist.SetXTitle("pixel #"); + fHist.SetYTitle("\\sigma"); + fHist.SetZTitle("N"); +} + +// -------------------------------------------------------------------------- +// +// Set binnings and prepare filling of the histogram. The binning for the +// pixel axis is derived automagically from the MPedestalCam container +// expected to be found in the MParList *plist. +// +Bool_t MHSigmaPixel::SetupFill(const MParList *plist) +{ + fPedestalCam = (MPedestalCam*)plist->FindObject("MPedestalCam"); + if (!fPedestalCam) + { + *fLog << err << dbginf << "MHSigmaPixel: MPedestalCam not found... aborting." << endl; + return kFALSE; + } + + MBinning* binssigma = (MBinning*)plist->FindObject("BinningSigma"); + MBinning* binspixel = new MBinning(); + binspixel->SetEdges(fPedestalCam->GetSize(), -0.5, -0.5+fPedestalCam->GetSize()); + + if (!binssigma) + { + *fLog << err << dbginf << "MHSigmaPixel: BinningSigma not found... aborting." << endl; + return kFALSE; + } + + SetBinning(&fHist, binspixel, binssigma); + + fHist.Sumw2(); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Fill the histogram +// +Bool_t MHSigmaPixel::Fill(const MParContainer *par) +{ + MPedestalCam &ped = *(MPedestalCam*)par; + for (Int_t i=0;i<(ped.GetSize());i++) + { + const MPedestalPix pix = ped[i]; + fHist.Fill(i, pix.GetSigma()); + } + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Draw the histogram +// +void MHSigmaPixel::Draw(Option_t *opt) +{ + //gStyle->SetOptStat(1000); + if (!gPad) + MakeDefCanvas("SigmaPixel", fTitle); + + fHist.Draw(opt); + + gPad->Modified(); + gPad->Update(); +} + +// -------------------------------------------------------------------------- +// +// Draw copies of the histogram +// +TObject *MHSigmaPixel::DrawClone(Option_t *opt) const +{ + //gStyle->SetOptStat(1000); + TCanvas &c = *MakeDefCanvas("SigmaPixel", fTitle); + + ((TH2&)fHist).DrawCopy(opt); + + c.Modified(); + c.Update(); + + return &c; +} + Index: /trunk/MagicSoft/Mars/mtemp/MHSigmaPixel.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/MHSigmaPixel.h (revision 1678) +++ /trunk/MagicSoft/Mars/mtemp/MHSigmaPixel.h (revision 1678) @@ -0,0 +1,46 @@ +#ifndef MARS_MHSigmaPixel +#define MARS_MHSigmaPixel + +#ifndef MARS_MH +#include "MH.h" +#endif + +#ifndef ROOT_TH2 +#include "TH2.h" +#endif + +class MMcEvt; +class MPedestalCam; +class TH2D; +class MParList; + + +class MHSigmaPixel : public MH +{ +private: + MPedestalCam *fPedestalCam; + TH2D fHist; + +public: + MHSigmaPixel(const char *name=NULL, const char *title=NULL); + + virtual Bool_t SetupFill(const MParList *pList); + virtual Bool_t Fill(const MParContainer *par); + + const TH2D *GetHist() { return &fHist; } + const TH2D *GetHist() const { return &fHist; } + + TH1 *GetHistByName(const TString name) { return &fHist; } + + void Draw(Option_t *option=""); + TObject *DrawClone(Option_t *option="") const; + + ClassDef(MHSigmaPixel, 1) //2D-histogram in Sigma and Pixel number +}; + +#endif + + + + + Index: /trunk/MagicSoft/Mars/mtemp/MHSigmabarTheta.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/MHSigmabarTheta.cc (revision 1678) +++ /trunk/MagicSoft/Mars/mtemp/MHSigmabarTheta.cc (revision 1678) @@ -0,0 +1,144 @@ +/* ======================================================================== *\ +! +! * +! * 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): Robert Wagner 10/2002 +! Copyright: MAGIC Software Development, 2000-2002 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// // +// MHSigmabarTheta // +// // +// 2D-Histogram in Sigmabar and Theta // +// // +////////////////////////////////////////////////////////////////////////////// + +#include "MHSigmabarTheta.h" + +#include + +#include + +#include "MMcEvt.hxx" +#include "MSigmabar.h" + +#include "MBinning.h" +#include "MParList.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MHSigmabarTheta); + +// -------------------------------------------------------------------------- +// +// Default Constructor. It sets name and title of the histogram. +// +MHSigmabarTheta::MHSigmabarTheta(const char *name, const char *title) +{ + // + // set the name and title of this object + // + fName = name ? name : "MHSigmabarTheta"; + fTitle = title ? title : "3-D histogram in sigmabar and theta"; + + fHist.SetDirectory(NULL); + + fHist.SetTitle("3D-plot of sigmabar and theta"); + fHist.SetXTitle("\\theta [\\circ]"); + fHist.SetYTitle("\\overline{\\sigma}"); + fHist.SetZTitle("N"); +} + +// -------------------------------------------------------------------------- +// +// Set binnings and prepare filling of the histogram +// +Bool_t MHSigmabarTheta::SetupFill(const MParList *plist) +{ + fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt"); + if (!fMcEvt) + { + *fLog << err << dbginf << "MHSigmabarTheta : MMcEvt not found... aborting." << endl; + return kFALSE; + } + + fSigmabar = (MSigmabar*)plist->FindObject("MSigmabar"); + if (!fSigmabar) + { + *fLog << err << dbginf << "MHSigmabarTheta : MSigmabar not found... aborting." << endl; + return kFALSE; + } + + MBinning* binssigmabar = (MBinning*)plist->FindObject("BinningSigmabar"); + MBinning* binstheta = (MBinning*)plist->FindObject("BinningTheta"); + if (!binssigmabar || !binstheta) + { + *fLog << err << dbginf << "MHSigmabarTheta : At least one MBinning not found... aborting." << endl; + return kFALSE; + } + + SetBinning(&fHist, binstheta, binssigmabar); + + fHist.Sumw2(); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Fill the histogram +// +Bool_t MHSigmabarTheta::Fill(const MParContainer *par) +{ + fHist.Fill(fMcEvt->GetTheta()*kRad2Deg, fSigmabar->GetSigmabar()); + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Draw the histogram +// +void MHSigmabarTheta::Draw(Option_t *opt) +{ + if (!gPad) + MakeDefCanvas("SigmabarTheta", fTitle); + + fHist.Draw(opt); + + gPad->Modified(); + gPad->Update(); +} + +// -------------------------------------------------------------------------- +// +// Draw copies of the histogram +// +TObject *MHSigmabarTheta::DrawClone(Option_t *opt) const +{ + TCanvas &c = *MakeDefCanvas("SigmabarTheta", fTitle); + + ((TH2&)fHist).DrawCopy(opt); + + c.Modified(); + c.Update(); + + return &c; +} + Index: /trunk/MagicSoft/Mars/mtemp/MHSigmabarTheta.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/MHSigmabarTheta.h (revision 1678) +++ /trunk/MagicSoft/Mars/mtemp/MHSigmabarTheta.h (revision 1678) @@ -0,0 +1,48 @@ +#ifndef MARS_MHSigmabarTheta +#define MARS_MHSigmabarTheta + +#ifndef MARS_MH +#include "MH.h" +#endif + +#ifndef ROOT_TH2 +#include "TH2.h" +#endif + +class MMcEvt; +class MSigmabar; +class TH2D; +class MParList; + + +class MHSigmabarTheta : public MH +{ +private: + MMcEvt *fMcEvt; + MSigmabar *fSigmabar; + + TH2D fHist; + +public: + MHSigmabarTheta(const char *name=NULL, const char *title=NULL); + + virtual Bool_t SetupFill(const MParList *pList); + virtual Bool_t Fill(const MParContainer *par); + + const TH2D *GetHist() { return &fHist; } + const TH2D *GetHist() const { return &fHist; } + + TH1 *GetHistByName(const TString name) { return &fHist; } + + void Draw(Option_t *option=""); + TObject *DrawClone(Option_t *option="") const; + + ClassDef(MHSigmabarTheta, 1) //3D-histogram in alpha, Energy and theta +}; + +#endif + + + + + Index: /trunk/MagicSoft/Mars/mtemp/MObservatoryLocation.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/MObservatoryLocation.cc (revision 1678) +++ /trunk/MagicSoft/Mars/mtemp/MObservatoryLocation.cc (revision 1678) @@ -0,0 +1,69 @@ +/* ======================================================================== *\ +! +! * +! * 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): Robert Wagner 10/2002 +! +! Copyright: MAGIC Software Development, 2002 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MObservatoryLocation // +// // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MObservatoryLocation.h" + +#include + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MObservatoryLocation); + +MObservatoryLocation::MObservatoryLocation(const char *name, const char *title) +{ + fName = name ? name : "MObservatoryLocation"; + fTitle = title ? title : "Storage container for coordinates of an observatory"; + // TH1F fHorizon=new TH1F; + fgDegToRad=2*TMath::Pi()/360; + fLatitude = 28.7594 * fgDegToRad; // rad; 28 45 34 + fLongitude = 17.8761 * fgDegToRad; // rad; 17 52 34; + // slalib uses + for WEST !!! + fElevation = 2300; // m + fObsName = "Observatorio del Roque de los Muchachos"; +} + +//Double_t GetHorizon(Double_t phi); +//TF1 SetHorizonLine() { return fHorizon; } +//TF1 GetHorizonLine() { return fHorizon; } + +MObservatoryLocation::~MObservatoryLocation() +{ + // do nothing special. +} + +void MObservatoryLocation::Print(Option_t *) const +{ + *fLog << all; + *fLog << fObsName << endl; + *fLog << "Latitude " << (fLatitude > 0 ? (fLatitude/fgDegToRad) : -(fLatitude/fgDegToRad)) << " deg " << (fLatitude > 0 ? "W" : "E") << endl; + *fLog << "Longitude " << (fLongitude > 0 ? (fLongitude/fgDegToRad) : -(fLongitude/fgDegToRad)) <<" deg " << (fLongitude < 0 ? "N" : "S") << endl; + *fLog << "Elevation " << fElevation << "m" << endl; +} Index: /trunk/MagicSoft/Mars/mtemp/MObservatoryLocation.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/MObservatoryLocation.h (revision 1678) +++ /trunk/MagicSoft/Mars/mtemp/MObservatoryLocation.h (revision 1678) @@ -0,0 +1,44 @@ +#ifndef MARS_MObservatoryLocation +#define MARS_MObservatoryLocation + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +class MObservatoryLocation : public MParContainer +{ +private: + char* fObsName; + Double_t fLatitude; + Double_t fLongitude; + Double_t fElevation; + Double_t fgDegToRad; + // TH1F fHorizon; + +public: + MObservatoryLocation(const char *name=NULL, const char *title=NULL); + ~MObservatoryLocation(); + + inline void SetLatitude(Double_t latitude) { fLatitude = latitude; } + inline void SetLongitude(Double_t longitude) { fLongitude = longitude; } + inline void SetElevation(Double_t elevation) { fElevation = elevation; } + inline void SetObservatoryName(char* name) { fObsName = name; } + + void MObservatoryLocation::Print(Option_t *) const; + + inline Double_t GetLatitude() { return fLatitude/fgDegToRad; } + inline Double_t GetLongitude() { return fLongitude/fgDegToRad; } + inline Double_t GetElevation() { return fElevation; } + Double_t GetLatitudeRad() { return fLatitude; } + Double_t GetLongitudeRad() { return fLongitude; } + char* GetObservatoryName() { return fObsName; } + // Double_t GetHorizon(Double_t phi); + // void SetHorizonLine(TF1 hor) { fHorizon = hor; } + // TH1F GetHorizonLine() { return fHorizon; } + + ClassDef(MObservatoryLocation, 1) + +}; + +#endif + Index: /trunk/MagicSoft/Mars/mtemp/MSigmabar.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/MSigmabar.cc (revision 1678) +++ /trunk/MagicSoft/Mars/mtemp/MSigmabar.cc (revision 1678) @@ -0,0 +1,150 @@ +/* ======================================================================== *\ +! +! * +! * 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): Robert Wagner 10/2002 +! +! Copyright: MAGIC Software Development, 2002 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MSigmabar // +// // +// This is the storage container to hold information about the mean sigma // +// (aka Sigmabar) of all pedestals // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MSigmabar.h" + +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MGeomCam.h" +#include "MPedestalCam.h" +#include "MGeomPix.h" +#include "MPedestalPix.h" + +ClassImp(MSigmabar); + +MSigmabar::MSigmabar(const char *name, const char *title) +{ + fName = name ? name : "MSigmabar"; + fTitle = title ? title : "Storage container for Sigmabar"; + + fSigmabar = 0.0; + fSigmabarInner = 0.0; + fSigmabarOuter = 0.0; + fRatioA = 0.0; + + fCalcPixNum=kTRUE; +} + +MSigmabar::~MSigmabar() +{ + // do nothing special. +} + + +// -------------------------------------------------------------------------- +// +// Actual calculation of sigmabar. This is done for each of the six sectors +// separately due to their possibly different HV behavior. Also inner and +// outer pixels are treated separately +// +// Preliminary! Works for CT1 test, for real MAGIC crosschecks still have +// to be done. Also implementation details will be updated, like +// determination of sector to which a respective pixel belongs +// +Bool_t MSigmabar::Calc(const MGeomCam &geom, const MPedestalCam &ped) +{ + const UInt_t npix = ped.GetSize(); + Float_t innerSquaredSum[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; + Float_t outerSquaredSum[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; + Int_t innerPixels[6] = {0,0,0,0,0,0}; + Int_t outerPixels[6] = {0,0,0,0,0,0}; + Int_t currentSector; + Float_t angle; + + for (UInt_t i=0; i 10/2002 +! +! Copyright: MAGIC Software Development, 2000-2002 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MSigmabarCalc // +// // +// This task calculates Sigmabar using the MSigmabar container and stores // +// its extremal values together with the corresponding Theta values // +// in MSigmabarParam. For the time being, Theta is taken from a Monte // +// Carlo container. This is preliminary and to be changed ASAP. // + // // +///////////////////////////////////////////////////////////////////////////// +#include "MSigmabarCalc.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" +#include "MGeomCam.h" +#include "MPedestalCam.h" +#include "MSigmabar.h" +#include "MSigmabarParam.h" +#include "MMcEvt.hxx" + +ClassImp(MSigmabarCalc); + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MSigmabarCalc::MSigmabarCalc(const char *name, const char *title) +{ + fName = name ? name : "MSigmabarCalc"; + fTitle = title ? title : "Task to calculate Sigmabar"; + + Reset(); +} + +MSigmabarCalc::~MSigmabarCalc() +{ + // nothing special yet. +} + +// -------------------------------------------------------------------------- +// +// check if necessary containers exists in the Parameter list already. +// if not create one and add them to the list +// +Bool_t MSigmabarCalc::PreProcess(MParList *pList) +{ + fCam = (MGeomCam*)pList->FindObject("MGeomCam"); + if (!fCam) + { + *fLog << dbginf << "MGeomCam not found (no geometry information available)... aborting." << endl; + return kFALSE; + } + + fPed = (MPedestalCam*)pList->FindObject("MPedestalCam"); + if (!fPed) + { + *fLog << dbginf << "MPedestalCam not found... aborting." << endl; + return kFALSE; + } + + fSig = (MSigmabar*)pList->FindCreateObj("MSigmabar"); + if (!fSig) + { + *fLog << dbginf << "MSigmabar not found... aborting." << endl; + return kFALSE; + } + + fSigParam = (MSigmabarParam*)pList->FindCreateObj("MSigmabarParam"); + if (!fSigParam) + { + *fLog << dbginf << "MSigmabarParam not found... aborting." << endl; + return kFALSE; + } + + fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); + // This is needed for determining min/max Theta + if (!fMcEvt) + { + *fLog << err << dbginf << "MHSigmabarTheta : MMcEvt not found... aborting." << endl; + return kFALSE; + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Calculating a new Sigmabar is not necessary on event basis as long as +// we deal with CT1 data. Therefore, the real calculation is done in +// the ReInit function. Process just takes care for finding the extremal +// values. Preliminary. +// +Bool_t MSigmabarCalc::Process() +{ + Bool_t rc = fSig->Calc(*fCam, *fPed); + fSigmabarMax = TMath::Max(fSig->GetSigmabar(), fSigmabarMax); + fSigmabarMin = TMath::Min(fSig->GetSigmabar(), fSigmabarMin); + + if (fMcEvt->GetTheta()*kRad2Deg < 90) + fThetaMax = TMath::Max(fMcEvt->GetTheta()*kRad2Deg, fThetaMax); + fThetaMin = TMath::Min(fMcEvt->GetTheta()*kRad2Deg, fThetaMin); + + return rc; +} + +// -------------------------------------------------------------------------- +// +// Calculates Sigmabar (for CT1 only needed at Reinit, i.e. when reading +// a new file) +// +Bool_t MSigmabarCalc::ReInit(MParList *pList) +{ + + fSigParam->SetParams(1, fSigmabarMin, fSigmabarMax, fThetaMin, fThetaMax); + Reset(); + + return kTRUE; +} + +void MSigmabarCalc::Reset() +{ + fThetaMin = 200; //there must be a function which gives me the hightest + fThetaMax = 0; // value allowed for a certain type! + fSigmabarMin = 200; + fSigmabarMax = 0; +} + Index: /trunk/MagicSoft/Mars/mtemp/MSigmabarCalc.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/MSigmabarCalc.h (revision 1678) +++ /trunk/MagicSoft/Mars/mtemp/MSigmabarCalc.h (revision 1678) @@ -0,0 +1,54 @@ +#ifndef MARS_MSigmabarCalc +#define MARS_MSigmabarCalc + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#ifndef MARS_MGeomCam +#include "MGeomCam.h" +#endif + +#ifndef MARS_MMcEvt +#include "MMcEvt.hxx" +#endif + +#ifndef MARS_MPedestalCam +#include "MPedestalCam.h" +#endif + +#ifndef MARS_MSigmabar +#include "MSigmabar.h" +#endif + +#ifndef MARS_MSigmabarParam +#include "MSigmabarParam.h" +#endif + +class MSigmabarCalc : public MTask +{ +private: + const MGeomCam *fCam; + const MPedestalCam *fPed; + MSigmabar *fSig; + Double_t fSigmabarMin; // Parametrization + Double_t fSigmabarMax; + Double_t fThetaMin; + Double_t fThetaMax; + MSigmabarParam *fSigParam; + MMcEvt *fMcEvt; + void MSigmabarCalc::Reset(); + +public: + MSigmabarCalc(const char *name=NULL, const char *title=NULL); + ~MSigmabarCalc(); + + Bool_t PreProcess(MParList *pList); + Bool_t ReInit(MParList *pList); + Bool_t Process(); + + ClassDef(MSigmabarCalc, 2) // task for calculating sigmabar +}; + +#endif + Index: /trunk/MagicSoft/Mars/mtemp/MSigmabarParam.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/MSigmabarParam.cc (revision 1678) +++ /trunk/MagicSoft/Mars/mtemp/MSigmabarParam.cc (revision 1678) @@ -0,0 +1,80 @@ +/* ======================================================================== *\ +! +! * +! * 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): Robert Wagner 10/2002 +! +! Copyright: MAGIC Software Development, 2002 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MSigmabarParam // +// // +// Storage Container for parameters characterizing a distribution of // +// events in the Sigmabar-Theta plane // +// // +// For CT1 tests, we just store Sigmabar_max, Sigmabar_min, Theta_max, // +// Theta_min. Later MJD and perhaps more than two points on the // +// distribution might follow. // +// // +// This implementation is still PRELIMINARY and requires some workarounds // +// put in SPECIFICALLY FOR THE CT1 TESTS, since a database to access is // +// missing. It is not the FINAL MAGIC VERSION. // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MSigmabarParam.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MSigmabarParam); + +MSigmabarParam::MSigmabarParam(const char *name, const char *title) : fSigmabarMin(0), fSigmabarMax(0), fMjdMin(0), fMjdMax(0), fThetaMin(0), fThetaMax(0), fRunNumber(0) +{ + fName = name ? name : "MSigmabarParam"; + fTitle = title ? title : "Storage container for characterizing a distribution of events in the Sigmabar-Theta plane"; +} + +MSigmabarParam::~MSigmabarParam() +{ + // nothing special yet +} + +void MSigmabarParam::SetRunNumber(Int_t r) +{ + fRunNumber = r; +} + +void MSigmabarParam::SetParams(Int_t r, Double_t si, Double_t sx, Double_t ti, Double_t tx, Double_t mi, Double_t mx) +{ + fSigmabarMin = si; + fSigmabarMax = sx; + fThetaMin = ti; + fThetaMax = tx; + fMjdMin = mi; + fMjdMax = mx; + //fRunNumber = r; +} + +void MSigmabarParam::Print() +{ + *fLog << endl << "Run " << fRunNumber << " | " + << "Sigmabar Min, Max: " << fSigmabarMin << " " << fSigmabarMax + << "| Theta Min, Max: " << fThetaMin << " " << fThetaMax << endl; +} Index: /trunk/MagicSoft/Mars/mtemp/MSigmabarParam.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/MSigmabarParam.h (revision 1678) +++ /trunk/MagicSoft/Mars/mtemp/MSigmabarParam.h (revision 1678) @@ -0,0 +1,42 @@ +#ifndef MARS_MSigmabarParam +#define MARS_MSigmabarParam + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +class MSigmabarParam : public MParContainer +{ +private: + Double_t fSigmabarMin; + Double_t fSigmabarMax; + Double_t fMjdMin; + Double_t fMjdMax; + Double_t fThetaMin; + Double_t fThetaMax; + Int_t fRunNumber; + +public: + + MSigmabarParam(const char *name=NULL, const char *title=NULL); + ~MSigmabarParam(); + + // void Print(Option_t *) const; + + void SetRunNumber(Int_t r); + + void SetParams(Int_t r, Double_t si, Double_t sx, Double_t ti, Double_t tx, Double_t mi=0, Double_t mx=0); + + Double_t GetSigmabarMin() { return fSigmabarMin; } + Double_t GetSigmabarMax() { return fSigmabarMax; } + Double_t GetThetaMin() { return fThetaMin; } + Double_t GetThetaMax() { return fThetaMax; } + + void Print(); + + ClassDef(MSigmabarParam, 1) // Storage container for characterizing a distribution of events in the Sigmabar-Theta plane + +}; + +#endif + Index: /trunk/MagicSoft/Mars/mtemp/MVPObject.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/MVPObject.cc (revision 1678) +++ /trunk/MagicSoft/Mars/mtemp/MVPObject.cc (revision 1678) @@ -0,0 +1,312 @@ +/* ======================================================================== *\ +! +! * +! * 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): Robert Wagner 10/2002 +! +! Copyright: MAGIC Software Development, 2000-2002 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MVPObject // +// // +// Class used by the visibility plotter to convert RA/Dec to Alt/Az // +// // +// This class represents an object and is used with the Visibility // +// macro. It must be provided with its RA/Dec coordinates and an // +// object name (cf. MVPObject::SetRA, MVPObject::SetDec, and // +// MVPObject::SetName). Alternatively, you can require the MVPObject // +// to be a solar system object like the Sun, Mars or the Moon // +// (cf. MVPObject::SetObject). // +// // +// MVPObject is ready to be used in a Mars Eventloop. You must provide // +// an Observatory Location as well as a time at which the position // +// of the MVPObject is to be calculated. MVPObject::PreProcess // +// checks the existence of the required containers and also makes sure // +// all necessary setters have been called. MVPObject::Process // +// then calculates the Alt/Az position of the object, as well as the // +// Zenith angle and the object diameter (Solar system objects). // +// // +// The astronomical algorithms used are taken from SLALIB 2.4-8. // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MVPObject.h" + +#include + +#include "MLog.h" +#include "MLogManip.h" +#include "MParList.h" + +#include "../../slalib/slalib.h" + +ClassImp(MVPObject); + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MVPObject::MVPObject(const char *name, const char *title) : fDiameter(0), fCalcEc(kFALSE), fUT1(52000), fBody(10), fGotRA(kFALSE), fGotDec(kFALSE), fGotName(kFALSE) +{ + fName = name ? name : "MVPObject"; + fTitle = title ? title : "Task to calculate Alt, Az of a given object"; + + fgDegToRad=2*TMath::Pi()/360; + fgHrsToRad=2*TMath::Pi()/24; +} + +MVPObject::~MVPObject() +{ + //Destructor: nothing special yet. +} + +// -------------------------------------------------------------------------- +// +// Check if necessary containers exist in the parameter list already. +// We need an ObservatoryLocation and a MVPTime object. +// +Bool_t MVPObject::PreProcess(MParList *pList) +{ + fObservatory = (MObservatoryLocation*)pList->FindObject("MObservatoryLocation"); + if (!fObservatory) + { + *fLog << dbginf << "MObservatoryLocation not found... aborting." << endl; + return kFALSE; + } + + fTime = (MVPTime*)pList->FindObject("MVPTime"); + if (!fTime) + { + *fLog << dbginf << "MVPTime not found... aborting." << endl; + return kFALSE; + } + + if (!fGotRA || !fGotDec || !fGotName) + { + *fLog << dbginf << "Object information is not complete." << endl; + return kFALSE; + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Sets coordinates from object name. Instead of providing RA, Dec and Name +// of an object, you may also just provide the object name in the from +// HHMMsDDT, where RA is given in hours and minutes and Declination is +// given by degrees DD and tenths of degrees T. "s" may be "+" or +// "-" +// +void MVPObject::SetObjectByName(char* object) +{ + fObjectName=object; + fGotName=kTRUE; + +// cout<<"OBJ:"< 9) + { + *fLog << dbginf << "No solar system object associated with value " << body <<"! Ignoring request." << endl; + return kFALSE; + } + else // We are working on a solar system body. + { + switch (body) + { + case 1: fObjectName="Mercury"; break; + case 2: fObjectName="Venus"; break; + case 3: fObjectName="Moon"; break; + case 4: fObjectName="Mars"; break; + case 5: fObjectName="Jupiter"; break; + case 6: fObjectName="Saturn"; break; + case 7: fObjectName="Uranus"; break; + case 8: fObjectName="Neptune"; break; + case 9: fObjectName="Pluto"; break; + default: fObjectName="Sun"; + } + } + + fBody = body; + fGotRA = fGotDec = fGotName = kTRUE; + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Given RA, Dec or a solar system object as well as an observatory +// location and a MVPTime, MVPObject::Process() calculates +// Alt, Az, ZA and (in the case of solar system objects) the apparent +// object diameter +// +Bool_t MVPObject::Process() +{ + Double_t diameter = 0.0; + + if (fBody < 10) // We are working on a solar system body. + { + slaRdplan(fTime->GetMJD(), fBody, fObservatory->GetLongitudeRad(), fObservatory->GetLatitudeRad(), &fRA, &fDec, &diameter); + } + + if (fCalcEc) slaEqecl(fRA, fDec, fTime->GetMJD(), &fEcLong, &fEcLat); + + Float_t azimuth; + Float_t elevation; + + Float_t hourAngle = (Float_t)UT1ToGMST(fTime->GetMJD()) - fRA; + + // cout << "ha: " << hourAngle << " ra: " << fRA << " dec " << fDec <GetLatitudeRad(), &azimuth, &elevation); + + fZA = slaZd(hourAngle, fDec, fObservatory->GetLatitudeRad()); + fAlt = (Double_t)elevation; + fAz = (Double_t)azimuth; + fDiameter = diameter; + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Returns distance of given object to this object in degrees +// +Double_t MVPObject::GetDistance(MVPObject* object) +{ + return slaSep(fRA, fDec, object->GetRARad(), object->GetDecRad())/fgDegToRad; +} + +// -------------------------------------------------------------------------- +// +// Returns distance of given object to this object in radians +// +Double_t MVPObject::GetDistanceRad(MVPObject* object) +{ + return slaSep(fRA, fDec, object->GetRARad(), object->GetDecRad()); +} + + +// -------------------------------------------------------------------------- +// +// Converts UT1 (given as MJD) to Greenwich mean star time in radians +// +Double_t MVPObject::UT1ToGMST(Double_t ut1) +{ + return slaGmst(ut1); +} + +void MVPObject::Print(Option_t *) const +{ + *fLog << all; + *fLog << "Position of "<< fObjectName << + ": Dec " << fDec/fgDegToRad << " deg, " << + "RA " << fRA/fgHrsToRad << " hrs" << endl; + if (fCalcEc) *fLog << "Ecliptic Long: " << fEcLong/fgDegToRad << " deg, " << + "Ecliptic Lat: " << fEcLat/fgDegToRad << " deg, " << endl; +} + + Index: /trunk/MagicSoft/Mars/mtemp/MVPObject.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/MVPObject.h (revision 1678) +++ /trunk/MagicSoft/Mars/mtemp/MVPObject.h (revision 1678) @@ -0,0 +1,100 @@ +#ifndef MARS_MVPObject +#define MARS_MVPObject + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#ifndef MARS_MObservatoryLocation +#include "MObservatoryLocation.h" +#endif + +#ifndef MARS_MVPTime +#include "MVPTime.h" +#endif + +class MVPObject : public MTask +{ +private: + MObservatoryLocation *fObservatory; + + Double_t fRA; + Double_t fDec; + Double_t fAlt; + Double_t fAz; + Double_t fZA; + Double_t fDiameter; + Double_t fEcLat; + Double_t fEcLong; + Bool_t fCalcEc; + + Double_t fUT1; + MVPTime *fTime; + + char* fObjectName; + UInt_t fBody; + + Double_t UT1ToGMST(Double_t ut1); + + Double_t fgDegToRad; + Double_t fgHrsToRad; + + Bool_t fGotRA; + Bool_t fGotDec; + Bool_t fGotName; + +public: + + MVPObject(const char *name=NULL, const char *title=NULL); + ~MVPObject(); + + Bool_t PreProcess(MParList *pList); + Bool_t Process(); + + void SetObservatory(MObservatoryLocation *observatory); + + void SetRA(Double_t rad) { fRA = rad; fBody = 10; fGotRA = kTRUE; } + void SetDec(Double_t rad) { fDec = rad; fBody = 10; fGotDec = kTRUE; } + + void SetRA(Int_t rh, Int_t rm, Int_t rs, Int_t ru = 0); + void SetDec(Int_t dh, Int_t dm, Int_t ds, Int_t du = 0); + + void SetCalcEc(Bool_t calcEc) { fCalcEc = calcEc; } + + void SetObjectName(char* name) { fObjectName = name; fGotName = kTRUE; } + void SetObjectByName(char* object); + + Bool_t SetObject(UInt_t body); + + Double_t GetRA() { return fRA/fgHrsToRad; } + Double_t GetDec() { return fDec/fgDegToRad; } + + Double_t GetRARad() { return fRA; } + Double_t GetDecRad() { return fDec; } + + Double_t GetZA() { return fZA; } + Double_t GetAltitude() { return fAlt; } + Double_t GetAzimut() { return fAz; } + Double_t GetDiameter() { return fDiameter; } + + Double_t GetEcLat() { return fEcLat; } + Double_t GetEcLong() { return fEcLong; } + + Double_t GetZADeg() { return fZA/fgDegToRad; } + Double_t GetAltitudeDeg() { return fAlt/fgDegToRad; } + Double_t GetAzimutDeg() { return fAz/fgDegToRad; } + + Double_t GetDistance(MVPObject* object); + Double_t GetDistanceRad(MVPObject* object); + + char* GetObjectName() { return fObjectName; } + + void Print(Option_t *) const; + + Double_t MJDStartOfYear(UInt_t year); + + ClassDef(MVPObject, 1) +}; + +#endif + Index: /trunk/MagicSoft/Mars/mtemp/MVPPlotter.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/MVPPlotter.cc (revision 1678) +++ /trunk/MagicSoft/Mars/mtemp/MVPPlotter.cc (revision 1678) @@ -0,0 +1,579 @@ +/* ======================================================================== *\ +! +! * +! * 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): Robert Wagner 12/2002 +! +! Copyright: MAGIC Software Development, 2000-2002 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MVPPlotter // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MVPPlotter.h" + +#include + +#include "TCanvas.h" +#include "TH1.h" +#include "TH2.h" +#include "TPaveText.h" +#include "TPaveLabel.h" +#include "TGraph.h" +#include "TString.h" +#include "TStyle.h" + +#include "MParList.h" +#include "MVPTime.h" + +ClassImp(MVPPlotter); + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MVPPlotter::MVPPlotter(const char *name, const char *title) : fUseSun(kTRUE), fUseMoon(kTRUE), fUsePlanets(kTRUE), fAstronomicalDarkness(-18.0) +{ + // fName = name ? name : "MVPPlotter"; + // fTitle = title ? title : "Generates visibility histograms and information"; + + fgSecPerDay = 86400; + fgMJD010170 = 40586; // 01-01-70 is JD 40586 + fgDegToRad = 2*TMath::Pi()/360; + +} + +MVPPlotter::~MVPPlotter() +// -------------------------------------------------------------------------- +// +// Destructor. Deletes objects if allocated. +// +{ + if (fUseSun) delete fSun; + if (fUseMoon) delete fMoon; + if (fUsePlanets) { + delete fVenus; + delete fMars; + delete fJupiter; + delete fSaturn; + } +} + +void MVPPlotter::SetupObjects() +// -------------------------------------------------------------------------- +// +// Create instances for Sun, Moon, Planets if requested. +// +{ + + fTime = new MVPTime(); + fPlist.AddToList(fTime); + fPlist.AddToList(fObs); + + if (fObject==NULL) { + cout << "You didn't specify an object!" << endl; + } + + fObject->PreProcess(&fPlist); + cout << "Processing object " << fObject->GetObjectName() << endl; + + if (fUseSun) { + fSun = new MVPObject(); + fSun->SetObject(0); // Sun + fSun->PreProcess(&fPlist); + fSun->SetCalcEc(kTRUE); + } + + if (fUseMoon) { + fMoon = new MVPObject(); + fMoon->SetObject(3); // Moon + fMoon->PreProcess(&fPlist); + fMoon->SetCalcEc(kTRUE); + } + + if (fUsePlanets) { + fVenus = new MVPObject(); + fVenus->SetObject(2); + fVenus->PreProcess(&fPlist); + fVenus->SetCalcEc(kTRUE); + + fMars = new MVPObject(); + fMars->SetObject(4); + fMars->PreProcess(&fPlist); + fMars->SetCalcEc(kTRUE); + + fJupiter = new MVPObject(); + fJupiter->SetObject(5); + fJupiter->PreProcess(&fPlist); + fJupiter->SetCalcEc(kTRUE); + + fSaturn = new MVPObject(); + fSaturn->SetObject(6); + fSaturn->PreProcess(&fPlist); + fSaturn->SetCalcEc(kTRUE); + } +} + +// -------------------------------------------------------------------------- +// +// Plots for a single object and a whole year are generated. +// +// Currently we do the following: +// - Create a plot MJD vs UTC for one year +// - Create a plot maxObjHeight vs Week for one year +// - Create visibility tables for one year (ZA, hours) +// +Bool_t MVPPlotter::CalcYear(Int_t year, UInt_t daySlices=96) +{ + SetupObjects(); + + UInt_t startday = (UInt_t)fTime->MJDStartOfYear(year); + UInt_t stopday = (UInt_t)fTime->MJDStartOfYear(year+1)-1; + + cout << "Processing period MJD "<< startday << " to MJD "<< stopday << endl; + + ULong_t startdayROOT = fgSecPerDay * (startday-fgMJD010170); + ULong_t stopdayROOT = fgSecPerDay * (stopday-fgMJD010170); + + // 2D Plot ZA vs MJD and UTC + fMjdUtcYear = new TH2D("fMjdUtcYear", "Visibility of object", + stopday-startday+1,startdayROOT,stopdayROOT, + daySlices+1,-450,fgSecPerDay+450); + + // Observability hours per day MJD + fMjdObsHours = new TH1D("fMjdObsHours", "Observation hours per day", + stopday-startday+1,startdayROOT,stopdayROOT); + + if (fUseMoon) { + // 2D Plot ZA of moon vs MJD and UTC + fMjdUtcYearMoon =new TH2D("fMjdUtcYearMoon", "Moon ZA", + stopday-startday+1,startdayROOT,stopdayROOT, + daySlices+1,-450,fgSecPerDay+450); + + // Moon phase vs MJD + fMjdMoonPhase =new TH1D("fMjdMoonPhase", "Moon phase", + stopday-startday+1,startdayROOT,stopdayROOT); + // Moon distance of object vs MJD + fMjdMoonDistance =new TH1D("fMjdMoonDistance", "Moon distance", + stopday-startday+1,startdayROOT,stopdayROOT); + // Moon intensity at object vs MJD + fMjdMoonIntensity =new TH1D("fMjdMoonIntensity", "Moon intensity at locus of object", + stopday-startday+1,startdayROOT,stopdayROOT); + + } + + if (fUsePlanets) { + // Distance of closest planet vs MJD + fMjdPlanetDistance =new TH1D("fMjdPlanetDistance", "PlanetDistance", + stopday-startday+1,startdayROOT,stopdayROOT); + } + + + // Array which holds total visibility time for whole year [0] and + // each month [1]..[12] + Float_t visibility[13][18]; + for (int m=0;m<13;m++) + for (int z=0;z<18;z++) + visibility[m][z]=0; + + int fday, ftime; + Double_t phase=0; + Double_t obsHours; + + for (UInt_t day=startday; day=(daySlices/2)) { + fday=day; + ftime=i-(daySlices/2); + } else { + fday=day-1; + ftime=i+(daySlices/2); + } + + // Objects access fTime via parameter list... + fTime->SetMJD(day,(Double_t)i/daySlices); + + if (fUseSun) fSun->Process(); + if (fUseMoon) fMoon->Process(); + + if (fUseSun && fUseMoon) { + + // Calculate moon phase... + phase = fSun->GetEcLong() - fMoon->GetEcLong(); + phase = TMath::Pi()-(acos(cos(phase)*cos(fMoon->GetEcLat()))); + phase = phase*180/TMath::Pi(); + todaysPhase+=phase; + + } + + // If sun is not up (or we should not use sun information...) + if (fSun->GetAltitudeDeg() < fAstronomicalDarkness) { + // Calculate Position of object: + fObject->Process(); + + // Calculate moon brightness at locus of object + // now this is gonna be fun... + + /* Evaluates predicted LUNAR part of sky brightness, in + * V magnitudes per square arcsecond, + */ + + moonIntensity = LunSkyBright(phase/180, fObject->GetDistance(fMoon), + fMoon->GetAltitudeDeg(), fObject->GetAltitudeDeg()); + fMjdMoonIntensity->Fill(fgSecPerDay*(fday-fgMJD010170),moonIntensity); + + + + // If moon is not up (or we should not use moon information...) + if ((!fUseMoon)||(fMoon->GetAltitudeDeg()<=0)||(moonIntensity<60)) { + // Fill MJD-UTC histogram + fMjdUtcYear->Fill(fgSecPerDay*(fday-fgMJD010170),fgSecPerDay*ftime/daySlices,fObject->GetAltitudeDeg()); + } + + // Sum up visibility time (only if moon not visible or moon + // info shouldn't be used at all) + if ((!fUseMoon)||(fMoon->GetAltitudeDeg()<=0)||(moonIntensity<60)) { + // Calculate + for (int z=0;z<18;z++) { + if (fObject->GetAltitudeDeg()>(z*5)) { + visibility[0][z]+=(Double_t)(60*24/daySlices); + visibility[(Int_t)fTime->GetMonth()][z]+=(Double_t)(60*24/daySlices); + } + } //for + + if ((fObject->GetAltitudeDeg())>40) { + fMjdObsHours->Fill(fgSecPerDay*(fday-fgMJD010170),(Double_t)(24/(Double_t)daySlices)); + } + + } + + + } //fi sun + + // Optional: Fill histo with moon-up times... + // if (fMoon->GetAltitudeDeg() >0) + // fMjdUtcYearMoon->Fill(fgSecPerDay*(day-fgMJD010170),fgSecPerDay*i/daySlices,phase); + // fMoon->Process(); + // Double_t phase; + // phase = fSun->GetEcLong() - moon->GetEcLong(); + // phase = TMath::Pi()-(acos(cos(phase)*cos(moon->GetEcLat()))); + // cout << "Phase: " << phase*180/TMath::Pi() << endl; + + } //for daySlices + + // Distance of Venus to object + + if (fUsePlanets) { + fVenus->Process(); + fMars->Process(); + fJupiter->Process(); + fSaturn->Process(); + + Double_t distance = fVenus->GetDistance(fObject); + distance = distance < fMars->GetDistance(fObject) ? distance : fMars->GetDistance(fObject); + distance = distance < fJupiter->GetDistance(fObject) ? distance : fJupiter->GetDistance(fObject); + distance = distance < fSaturn->GetDistance(fObject) ? distance : fSaturn->GetDistance(fObject); + + fMjdPlanetDistance->Fill(fgSecPerDay*(fday-fgMJD010170),distance); + } + + fMjdMoonPhase->Fill(fgSecPerDay*(fday-fgMJD010170),todaysPhase/i); + fMjdMoonDistance->Fill(fgSecPerDay*(fday-fgMJD010170),fMoon->GetDistance(fObject)); + + + } //for days + + + // Here we print the tables with visibilities... + cout << "Visibility time [hours]: " << endl; + + for (int z=1;z<17;z++) { + if (visibility[0][z]==0) break; + printf("Alt>%2d|%6d|",z*5,(Int_t)(visibility[0][z]/60)); + for (int m=1;m<13;m++) { + printf("%5d ",(Int_t)(visibility[m][z]/60)); + } + cout<20)&&((visibility[n][8])/60 <= 20)) { + vistimestart=m; // we're on the rising slope + } + } + + for (int m=1;m<13;m++) { + int n = (m==1 ? 12 : m-1); + if (((visibility[m][8]/60)<20)&&((visibility[n][8]/60)>=20)) { + vistimeend=n; // we're on the falling slope + } + } + + cout << "Visibility (Alt>40) during months: " << vistimestart << "-" << vistimeend << " approx " << visibility[0][8]/60 << "hrs" <Reset(); + + TCanvas *cMjdUtcYear = new TCanvas("cMjdUtcYear", "Object Visibility MjdUtcYear", 1100,500); + cMjdUtcYear->cd(0); + + gStyle->SetPalette(1); + // gStyle->SetOptStat(0); + gStyle->SetOptFit(0); + gStyle->SetFillStyle(0); + gStyle->SetFillColor(10); + gStyle->SetCanvasColor(10); + gStyle->SetDrawBorder(0); + gStyle->SetPadColor(10); + gStyle->SetPadBorderSize(0); + gStyle->SetPadLeftMargin(0.12); + gStyle->SetTitleYOffset(1.2); + gStyle->SetTitleXOffset(1.2); + + gStyle->SetPadLeftMargin(0.01); + gStyle->SetPadRightMargin(0.1); + gStyle->SetPadTopMargin(0.03); + gStyle->SetPadBottomMargin(0.12); + + fMjdUtcYear->SetTitle(fObject->GetObjectName()); + fMjdUtcYear->GetXaxis()->SetTimeDisplay(1); + fMjdUtcYear->GetXaxis()->SetTimeFormat("%b %y"); + fMjdUtcYear->GetYaxis()->SetTimeDisplay(1); + fMjdUtcYear->GetYaxis()->SetTimeFormat("%Hh"); + gStyle->SetTimeOffset(43200); + fMjdUtcYear->GetYaxis()->SetLabelOffset(0.01); + fMjdUtcYear->SetMinimum(0); + fMjdUtcYear->SetMaximum(90); + + fMjdUtcYear->GetYaxis()->SetTitle("Hour"); + fMjdUtcYear->GetXaxis()->SetTitle("Day of year"); + fMjdUtcYear->GetZaxis()->SetTitle("Altitude"); + fMjdUtcYear->Draw("SURF2BB9Z"); + gPad->SetPhi(0); + gPad->SetTheta(90); + gPad->Modified(); + + TPaveText *pt = new TPaveText(-0.74,-0.35,-0.58,0.52,""); + char ptLine[80]; + pt->AddText("Visibility time [hrs]:"); + pt->SetTextAlign(13); + pt->SetFillStyle(0); + pt->SetBorderSize(0); + pt->SetTextFont(42); + for (int j=1;j<17;j++) { + sprintf (ptLine, "Alt>%i: %i", j*5, (Int_t)visibility[0][j]/60); + pt->AddText(ptLine); + if (visibility[0][j]==0) break; + } + pt->Draw(); + + + if (fUseMoon) { + TCanvas *cMjdMoon = new TCanvas("cMjdMoon", "the Moon phase", 600,600); + gStyle->SetOptTitle(1); + cMjdMoon->Divide(1,3); + cMjdMoon->cd(1); + fMjdMoonPhase->Draw(); + cMjdMoon->cd(2); + fMjdMoonDistance->Draw(); + cMjdMoon->cd(3); + fMjdMoonIntensity->Draw(); + } + + TCanvas *cObsHours = new TCanvas("cObsHours", "ObsHours per day", 600,400); + fMjdObsHours->Draw(); + + + + if (fUsePlanets) { + TCanvas *cMjdPlanets = new TCanvas("cMjdPlanets", "the distance to planets", 600,200); + cMjdPlanets->cd(0); + fMjdPlanetDistance->Draw(); + } + + + // next histogram + Float_t objectHeight[55]; + Float_t simpleCounter[55]; + Int_t weekCounter=0; + + simpleCounter[weekCounter]=0; + objectHeight[weekCounter]=0; + weekCounter++; + + Int_t startday2 = fTime->MJDStartOfYear(year); + Int_t stopday2 = fTime->MJDStartOfYear(year+1)-1; + + for (int day=startday2; day=48) { + fday=day; + ftime=i-48; + } else { + fday=day-1; + ftime=i+48; + } + + fTime->SetMJD(day,(Double_t)i/daySlices); + fObject->Process(); + fSun->Process(); + + if (fSun->GetAltitudeDeg() < -25.0) { + if (objectHeight[weekCounter]<(fObject->GetAltitudeDeg())) + objectHeight[weekCounter]=fObject->GetAltitudeDeg(); + } + + } //i + weekCounter++; + } //day + simpleCounter[weekCounter]=weekCounter-2; + objectHeight[weekCounter]=0; + weekCounter++; + + TString months[12] = {"Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"}; + + TCanvas *c2 = new TCanvas("c2", "Object Visibility",600,100); + + // gStyle->SetOptTitle(0); +// gStyle->SetPadLeftMargin(0.000001); +// gStyle->SetPadRightMargin(0.000001); +// gStyle->SetPadTopMargin(0.001); +// gStyle->SetPadBottomMargin(0.000001); + gPad->SetGrid(); + + c2->SetGrid(); + TGraph *tg=new TGraph(weekCounter,simpleCounter,objectHeight); + + tg->SetMinimum(1); + tg->SetMaximum(90); + + Double_t maxza = abs(fObject->GetDec()-fObs->GetLatitude()); + + if (maxza > 90) maxza=90; + if (maxza < 0) maxza=00; + + cout << "MaxZA: "; + cout << maxza << endl; + + if (maxza < 30) { //colors=green to yellow + maxza=9*maxza/30; + maxza=80+maxza; + + } else { //colors=yellow to red + maxza-=30; + maxza=11*maxza/60; + maxza=89+maxza; + } + + tg->SetFillColor((Int_t)maxza); + tg->SetLineColor((Int_t)maxza); + tg->SetLineWidth((Int_t)maxza); + tg->Draw("AF"); + + tg->GetXaxis()->SetLimits(0,52); + tg->GetXaxis()->SetTickLength(0.1); + tg->GetXaxis()->SetNdivisions(0); + tg->GetYaxis()->SetNdivisions(202); + + TPaveLabel* l= new TPaveLabel(2,46,35,88, fObject->GetObjectName()); + l->SetBorderSize(0); + l->SetFillStyle(0); + l->SetTextAlign(12); + l->Draw(); + + if ((vistimestart<13)&&(vistimeend>0)) { + TString label; + label=months[vistimestart-1]+"-"+months[vistimeend-1]; + TPaveLabel* l2= new TPaveLabel(35,46,50,88, label); + l2->SetBorderSize(0); + l2->SetFillStyle(0); + l2->SetTextAlign(32); + l2->Draw(); + + } + + c2->Modified(); + c2->Update(); + +return kTRUE; +} + +Double_t MVPPlotter::LunSkyBright(Double_t moon_phase,Double_t rho,Double_t altmoon,Double_t alt) +{ +/* Evaluates predicted LUNAR part of sky brightness, in + * V magnitudes per square arcsecond, following K. Krisciunas + * and B. E. Schaeffer (1991) PASP 103, 1033. + * + * moon_phase = Phase of the Moon, between 0. (no moon) and 1. (full moon), + * rho (deg) = separation of moon and object, + * altmoon (deg) = altitude of moon above horizon, + * alt (deg) = altitude of object above horizon + */ + + double kzen=1.; + + double rho_rad = rho*fgDegToRad; + double alpha = 180.*(1. - moon_phase); + double Zmoon = (90. - altmoon)*fgDegToRad; + double Z = (90. - alt)*fgDegToRad; + + double istar = -0.4*(3.84 + 0.026*fabs(alpha) + 4.0e-9*pow(alpha,4.)); /*eqn 20*/ + istar = pow(10.,istar); + if(fabs(alpha) < 7.) /* crude accounting for opposition effect */ + istar = istar * (1.35 - 0.05 * fabs(istar)); + /* 35 per cent brighter at full, effect tapering linearly to + zero at 7 degrees away from full. mentioned peripherally in + Krisciunas and Scheafer, p. 1035. */ + double fofrho = 229087. * (1.06 + cos(rho_rad)*cos(rho_rad)); + if(fabs(rho) > 10.) + fofrho=fofrho+pow(10.,(6.15 - rho/40.)); /* eqn 21 */ + else if (fabs(rho) > 0.25) + fofrho= fofrho+ 6.2e7 / (rho*rho); /* eqn 19 */ + else fofrho = fofrho+9.9e8; /*for 1/4 degree -- radius of moon! */ + double Xzm = sqrt(1.0 - 0.96*sin(Zmoon)*sin(Zmoon)); + if(Xzm != 0.) Xzm = 1./Xzm; + else Xzm = 10000.; + double Xo = sqrt(1.0 - 0.96*sin(Z)*sin(Z)); + if(Xo != 0.) Xo = 1./Xo; + else Xo = 10000.; + double Bmoon = fofrho * istar * pow(10.,(-0.4*kzen*Xzm)) + * (1. - pow(10.,(-0.4*kzen*Xo))); /* nanoLamberts */ + // cout << " Bmoon=" << Bmoon; + if(Bmoon > 0.001) + return(Bmoon); + else return(99999.); +} + + Index: /trunk/MagicSoft/Mars/mtemp/MVPPlotter.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/MVPPlotter.h (revision 1678) +++ /trunk/MagicSoft/Mars/mtemp/MVPPlotter.h (revision 1678) @@ -0,0 +1,85 @@ +#ifndef MARS_MVPPlotter +#define MARS_MVPPlotter + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#ifndef MARS_MParList +#include "MParList.h" +#endif + +#ifndef MARS_MTaskList +#include "MTaskList.h" +#endif + +#ifndef MARS_MVPObject +#include "MVPObject.h" +#endif + +#ifndef MARS_MVPTime +#include "MVPTime.h" +#endif + +#ifndef MARS_MObservatoryLocation +#include "MObservatoryLocation.h" +#endif + +class TH1D; +class TH2D; +class TCanvas; + +class MVPPlotter +{ +private: + + MVPTime *fTime; + MObservatoryLocation *fObs; + + MParList fPlist; + MTaskList fTlist; + + Bool_t fUseSun; + Bool_t fUseMoon; + Bool_t fUsePlanets; + + Double_t fAstronomicalDarkness; + + Int_t fgSecPerDay; + Int_t fgMJD010170; // 01-01-70 is JD 40586 + Double_t fgDegToRad; + + TH2D *fMjdUtcYear; + TH2D *fMjdUtcYearMoon; + TH1D *fMjdMoonPhase; + TH1D *fMjdMoonDistance; + TH1D *fMjdPlanetDistance; + TH1D *fMjdMoonIntensity; + TH1D *fMjdObsHours; + + MVPObject *fObject; + MVPObject *fSun; + MVPObject *fMoon; + MVPObject *fVenus; + MVPObject *fMars; + MVPObject *fJupiter; + MVPObject *fSaturn; + + void SetupObjects(); + +public: + MVPPlotter(const char *name=NULL, const char *title=NULL); + ~MVPPlotter(); + + void SetObject(MVPObject *o) { fObject=o; } + void SetObservatory(MObservatoryLocation *l) { fObs=l; } + void SetAstronomicalDarkness(Double_t d) { fAstronomicalDarkness=d; } + + Double_t LunSkyBright(Double_t moon_phase,Double_t rho,Double_t altmoon,Double_t alt); + Bool_t CalcYear(Int_t year, UInt_t daySlices); + + ClassDef(MVPPlotter, 2) // Visibility Plotter: The Plotter Routine itself +}; + +#endif + Index: /trunk/MagicSoft/Mars/mtemp/MVPTime.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/MVPTime.cc (revision 1678) +++ /trunk/MagicSoft/Mars/mtemp/MVPTime.cc (revision 1678) @@ -0,0 +1,125 @@ +/* ======================================================================== *\ +! +! * +! * 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): Robert Wagner 11/2002 +! +! Copyright: MAGIC Software Development, 2002 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MVPTime // +// // +// Time Container used by the visibility plotter storing times in MJD // +// // +// This is a simple time container for the Visibility Plotter. It takes // +// "normal" or MJD-formatted time, converts it to MJD and can return MJD. // +// // +// The astronomical algorithms used are taken from SLALIB 2.4-8. // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MVPTime.h" + +#include + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MVPTime); + +MVPTime::MVPTime(const char *name, const char *title) +{ + fName = name ? name : "MVPTime"; + fTitle = title ? title : "Storage container for time"; +} + +MVPTime::~MVPTime() +{ + // do nothing special. +} + +void MVPTime::SetTime(UInt_t yr, UInt_t mo, UInt_t dy, UInt_t hr, UInt_t mi, UInt_t se) +{ + fUT1 = CalendarToUT1(yr, mo, dy, hr, mi, se); + fYear = yr; + fMonth = mo; + fDay = dy; +} + +void MVPTime::SetMJD(Double_t mjd, Double_t fracMjd) +{ + fUT1 = mjd + fracMjd; + Double_t fd; + Int_t status; + slaDjcl (mjd, &fYear, &fMonth, &fDay, &fd, &status); +} + +// -------------------------------------------------------------------------- +// +// +Double_t MVPTime::CalendarToUT1(UInt_t yr, UInt_t mo, UInt_t dy, UInt_t hr, UInt_t mi, UInt_t se) +{ + int status; + Double_t mjd; + slaCldj(yr, mo, dy, &mjd, &status); + switch (status) + { + case 1: + *fLog << "CalendarToMJD Warn: bad year" << endl; + break; + case 2: + *fLog << "CalendarToMJD Warn: bad month" << endl; + break; + case 3: + *fLog << "CalendarToMJD Warn: bad day" << endl; + break; + } + + mjd += ((Double_t)hr/24) + ((Double_t)mi/(24*60)) + ((Double_t)se/(24*60*60)); + return mjd; +} + +// -------------------------------------------------------------------------- +// +// +Double_t MVPTime::MJDStartOfYear(UInt_t year) +{ + int status; + Double_t mjd; + slaCldj(year, 1, 1, &mjd, &status); + switch (status) + { + case 1: + *fLog << "CalendarToMJD Warn: bad year" << endl; + break; + case 2: + *fLog << "CalendarToMJD Warn: bad month" << endl; + break; + case 3: + *fLog << "CalendarToMJD Warn: bad day" << endl; + break; + } + return mjd; +} + +void MVPTime::Print(Option_t *) const +{ + *fLog << all; + *fLog << "Time (MJD) is: " << fUT1 << endl; +} Index: /trunk/MagicSoft/Mars/mtemp/MVPTime.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/MVPTime.h (revision 1678) +++ /trunk/MagicSoft/Mars/mtemp/MVPTime.h (revision 1678) @@ -0,0 +1,37 @@ +#ifndef MARS_MVPTime +#define MARS_MVPTime + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +class MVPTime : public MParContainer +{ +private: + Double_t fUT1; + Int_t fYear; + Int_t fMonth; + Int_t fDay; + Double_t MVPTime::CalendarToUT1(UInt_t yr, UInt_t mo, UInt_t dy, UInt_t hr, UInt_t mi, UInt_t se); + +public: + MVPTime(const char *name=NULL, const char *title=NULL); + ~MVPTime(); + + void SetTime(UInt_t yr, UInt_t mo, UInt_t dy, UInt_t hr, UInt_t mi, UInt_t se); + void SetMJD(Double_t mjd, Double_t fracMjd = 0); + inline Double_t GetMJD() { return fUT1; } + + void Print(Option_t *) const; + + Double_t MJDStartOfYear(UInt_t year); + Double_t GetYear() { return fYear; } + Double_t GetMonth() { return fMonth; } + Double_t GetDay() { return fDay; } + + ClassDef(MVPTime, 1) + +}; + +#endif + Index: /trunk/MagicSoft/Mars/mtemp/Makefile =================================================================== --- /trunk/MagicSoft/Mars/mtemp/Makefile (revision 1678) +++ /trunk/MagicSoft/Mars/mtemp/Makefile (revision 1678) @@ -0,0 +1,57 @@ +################################################################## +# +# makefile +# +# for the MARS software +# +################################################################## +include ../Makefile.conf.$(OSTYPE) +include ../Makefile.conf.general + +# +# Handling name of the Root Dictionary Files +# +CINT = Temp + +# +# Library name to creatre +# +LIB = mtemp.a + +# +# connect the include files defined in the config.mk file +# +INCLUDES = -I. -I../mbase -I../mraw -I../manalysis -I../mmc \ + -I../mgui -I../mgeom -I../mdata -I../mhist -I../../slalib + +.SUFFIXES: .c .cc .cxx .h .hxx .o + +SRCFILES = \ + MApplyPadding.cc \ + MHSigmaPixel.cc \ + MHSigmabarTheta.cc \ + MObservatoryLocation.cc \ + MSigmabar.cc \ + MSigmabarCalc.cc \ + MSigmabarParam.cc \ + MVPObject.cc \ + MVPPlotter.cc \ + MVPTime.cc + +SRCS = $(SRCFILES) +HEADERS = $(SRCFILES:.cc=.h) +OBJS = $(SRCFILES:.cc=.o) + +############################################################ + +all: $(LIB) + +include ../Makefile.rules + +clean: rmcint rmobjs rmcore rmlib + +mrproper: clean rmbak + +# @endcode + + Index: /trunk/MagicSoft/Mars/mtemp/TempIncl.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/TempIncl.h (revision 1678) +++ /trunk/MagicSoft/Mars/mtemp/TempIncl.h (revision 1678) @@ -0,0 +1,5 @@ +#ifndef __CINT__ + +#include + +#endif // __CINT__ Index: /trunk/MagicSoft/Mars/mtemp/TempLinkDef.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/TempLinkDef.h (revision 1678) +++ /trunk/MagicSoft/Mars/mtemp/TempLinkDef.h (revision 1678) @@ -0,0 +1,19 @@ +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class MApplyPadding+; + +#pragma link C++ class MHSigmaPixel+; +#pragma link C++ class MHSigmabarTheta+; + +#pragma link C++ class MSigmabar+; +#pragma link C++ class MSigmabarCalc+; +#pragma link C++ class MSigmabarParam+; + +#pragma link C++ class MObservatoryLocation+; +#pragma link C++ class MVPObject+; +#pragma link C++ class MVPPlotter+; +#pragma link C++ class MVPTime+;