/* ======================================================================== *\ ! ! * ! * 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, 4/2004 ! ! Copyright: MAGIC Software Development, 2000-2004 ! ! \* ======================================================================== */ ////////////////////////////////////////////////////////////////////////////// // // MHCamEventRot // // Derotate MCamEvent before filling a 2D histogram. // // To be done: Set the sky position of the center of the display and // correct if the pointing position of the telescope is // different // ////////////////////////////////////////////////////////////////////////////// #include "MHCamEventRot.h" #include #include #include #include #include #include #include #include "MGeomCam.h" #include "MSrcPosCam.h" #include "MHillasSrc.h" #include "MTime.h" #include "MObservatory.h" #include "MPointingPos.h" #include "MAstroCatalog.h" #include "MAstroSky2Local.h" #include "MStatusDisplay.h" #include "MMath.h" #include "MBinning.h" #include "MParList.h" #include "MCamEvent.h" #include "MLog.h" #include "MLogManip.h" ClassImp(MHCamEventRot); using namespace std; // -------------------------------------------------------------------------- // // Default Constructor // MHCamEventRot::MHCamEventRot(const char *name, const char *title) : fTime(0), fPointPos(0), fObservatory(0), fType(0), fNameTime("MTime") { // // set the name and title of this object // fName = name ? name : "MHCamEventRot"; fTitle = title ? title : "Plot showing derotated MCamEvent data"; fHist.SetDirectory(NULL); fHist.SetName("Alpha"); fHist.SetTitle("2D-plot of MCamEvents (derotated)"); fHist.SetXTitle("x [\\circ]"); fHist.SetYTitle("y [\\circ]"); fHist.SetZTitle("Counts"); } // -------------------------------------------------------------------------- // // Set binnings (takes BinningCamEvent) and prepare filling of the // histogram. // // Also search for MTime, MObservatory and MPointingPos // Bool_t MHCamEventRot::SetupFill(const MParList *plist) { fGeom = (MGeomCam*)plist->FindObject("MGeomCam"); if (!fGeom) { *fLog << err << "MGeomCam not found... aborting." << endl; return kFALSE; } const MBinning *bins = (MBinning*)plist->FindObject("BinningCamEvent"); if (!bins) { const Float_t r = fGeom->GetMaxRadius()*fGeom->GetConvMm2Deg(); MBinning b; b.SetEdges(41, -r, r); SetBinning(&fHist, &b, &b); } else SetBinning(&fHist, bins, bins); fPointPos = (MPointingPos*)plist->FindObject(AddSerialNumber("MPointingPos")); if (!fPointPos) *fLog << warn << "MPointingPos not found... no derotation." << endl; fTime = (MTime*)plist->FindObject(AddSerialNumber(fNameTime)); if (!fTime) *fLog << warn << "MTime not found... no derotation." << endl; fObservatory = (MObservatory*)plist->FindObject(AddSerialNumber("MObservatory")); if (!fObservatory) *fLog << warn << "MObservatory not found... no derotation." << endl; // FIXME: Because the pointing position could change we must check // for the current pointing position and add a offset in the // Fill function! if (fPointPos) { fRa = fPointPos->GetRa(); fDec = fPointPos->GetDec(); } return kTRUE; } // -------------------------------------------------------------------------- // // Fill the histogram. For details see the code or the class description // Bool_t MHCamEventRot::Fill(const MParContainer *par, const Stat_t w) { const MCamEvent *evt = dynamic_cast(par); if (!evt) { *fLog << err << "MHCamEventRot::Fill: No container specified!" << endl; return kFALSE; } // Get max radius... // const Double_t maxr = fGeom->GetMaxRadius()*fGeom->GetConvMm2Deg(); // Get camera rotation angle Double_t rho = 0; if (fTime && fObservatory && fPointPos) rho = fPointPos->RotationAngle(*fObservatory, *fTime); //if (fPointPos) // rho = fPointPos->RotationAngle(*fObservatory); // Get number of bins and bin-centers const Int_t nx = fHist.GetNbinsX(); const Int_t ny = fHist.GetNbinsY(); Axis_t cx[nx]; Axis_t cy[ny]; fHist.GetXaxis()->GetCenter(cx); fHist.GetYaxis()->GetCenter(cy); for (int ix=0; ixmaxr) // continue; // rotate center of bin TVector2 v(cx[ix], cy[iy]); if (rho!=0) v=v.Rotate(rho); // FIXME: slow! Create Lookup table instead! const Int_t idx = fGeom->GetPixelIdxDeg(v); Double_t val; if (idx<0 || !evt->GetPixelContent(val, idx, *fGeom, fType)) continue; // Fill histogram fHist.Fill(cx[ix], cy[iy], val); } } return kTRUE; } // -------------------------------------------------------------------------- // // Get the MAstroCatalog corresponding to fRa, fDec. The limiting magnitude // is set to 9, while the fov is equal to the current fov of the false // source plot. // TObject *MHCamEventRot::GetCatalog() { const Double_t maxr = 0.98*TMath::Abs(fHist.GetBinCenter(1)); // Create catalog... MAstroCatalog stars; stars.SetLimMag(9); stars.SetGuiActive(kFALSE); stars.SetRadiusFOV(maxr); stars.SetRaDec(fRa*TMath::DegToRad()*15, fDec*TMath::DegToRad()); stars.ReadBSC("bsc5.dat"); TObject *o = (MAstroCatalog*)stars.Clone(); o->SetBit(kCanDelete); return o; } // -------------------------------------------------------------------------- // // Draw the histogram // void MHCamEventRot::Draw(Option_t *opt) { TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); pad->SetBorderMode(0); AppendPad(""); // draw the 2D histogram Sigmabar versus Theta fHist.Draw(opt); TObject *catalog = GetCatalog(); catalog->Draw("mirror same"); }