Index: trunk/MagicSoft/Mars/mtemp/mmpi/MPointingPosCalcWA.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mmpi/MPointingPosCalcWA.cc (revision 5490) +++ trunk/MagicSoft/Mars/mtemp/mmpi/MPointingPosCalcWA.cc (revision 5490) @@ -0,0 +1,181 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Thomas Bretz, 11/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MPointingPosCalc +// +// Currently: +// +// * MC files: Copy the simulated telescope position (Telescope Theta, +// Telescope Phi in MMcEvt) to MPointingPosition +// +// * Real Data: Copy the nominal poiting position (Nominal Zd, Nominal Az +// in MReportDrive) to MPointingPosition +// +// Future: Interpolate the pointing position for each event between two +// consecutive drive reports. +// +// Input Container: +// MRawRunHeader +// [MMcEvt, MReportDrive] +// +// Output Container: +// MPointingPosition +// +///////////////////////////////////////////////////////////////////////////// +#include "MPointingPosCalcWA.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" +#include "MVector3.h" + +#include "MPointingPos.h" +#include "MAstro.h" +#include "MObservatory.h" +#include "MAstroSky2Local.h" +#include "MRawRunHeader.h" +#include "MReportDrive.h" +#include "MMcEvt.hxx" + +ClassImp(MPointingPosCalcWA); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Search for MRawRunHeader. Get the run type from there. Depending on +// the run type search either for MMcEvt or MReportDrive. +// +Bool_t MPointingPosCalcWA::ReInit(MParList *plist) +{ + fObs = (MObservatory*)plist->FindCreateObj("MObservatory"); + if (!fObs) + { + *fLog << err << "MObservatory not found... aborting." << endl; + return kFALSE; + } + + fRunHeader = (MRawRunHeader*)plist->FindObject("MRawRunHeader"); + if (!fRunHeader) + { + *fLog << err << "MRawRunHeader not found... aborting." << endl; + return kFALSE; + } + + fRunType = fRunHeader->GetRunType(); + + switch (fRunType) + { + case MRawRunHeader::kRTData: + fReport = (MReportDrive*)plist->FindObject("MReportDrive"); + if (!fReport) + { + *fLog << err << "MReportDrive not found... will set Az and Zd to a fixed value." << endl; + return kTRUE; + } + return kTRUE; + + case MRawRunHeader::kRTMonteCarlo: + fMcEvt = (MMcEvt*)plist->FindObject("MMcEvt"); + if (!fMcEvt) + { + *fLog << err << "MMcEvt not found... aborting." << endl; + return kFALSE; + } + return kTRUE; + + case MRawRunHeader::kRTPedestal: + *fLog << err << "Cannot work in a pedestal Run!... aborting." << endl; + return kFALSE; + + case MRawRunHeader::kRTCalibration: + *fLog << err << "Cannot work in a calibration Run!... aborting." << endl; + return kFALSE; + + default: + *fLog << err << "Run Type " << fRunType << " unknown!... aborting." << endl; + return kFALSE; + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Search for 'MPointingPos'. Create if not found. +// +Int_t MPointingPosCalcWA::PreProcess(MParList *plist) +{ + fPosition = (MPointingPos*)plist->FindCreateObj("MPointingPos"); + return fPosition ? kTRUE : kFALSE; +} + +// -------------------------------------------------------------------------- +// +// See class description. +// +Int_t MPointingPosCalcWA::Process() +{ + + const char plus = '+'; + Char_t Rsgn; + UShort_t Rhour, Rmin, Rsec; + Double_t ra; + + switch (fRunType) + { + case MRawRunHeader::kRTData: + if(!fReport) + { + MTime start = fRunHeader->GetRunStart(); + + MVector3 v; + v.SetRaDec(fRa,fDec); + v *= MAstroSky2Local(start,*fObs); + Double_t ZA = v.Theta()*180./acos(-1.); + Double_t Az = v.Phi()*180./acos(-1.); +// cout << " Za = " << ZA << " v.Theta() " << v.Theta() << " Az " << Az << " v.Phi " << v.Phi() << endl; + fPosition->SetLocalPosition(ZA, Az); + MAstro::Rad2Hms(fRa, Rsgn, Rhour, Rmin, Rsec); + if (Rsgn==plus) ra = Rhour + Rmin/60. + Rsec/3600.; + else ra = - (Rhour + Rmin/60. + Rsec/3600.); + fPosition->SetSkyPosition(ra, TMath::RadToDeg()*fDec); + return kTRUE; + } + else + { + fPosition->SetLocalPosition(fReport->GetNominalZd(), fReport->GetNominalAz()); + fPosition->SetSkyPosition(fReport->GetRa(), fReport->GetDec()); + return kTRUE; + } + + case MRawRunHeader::kRTMonteCarlo: + fPosition->SetLocalPosition(fMcEvt->GetTelescopeTheta()*TMath::RadToDeg(), fMcEvt->GetTelescopePhi()*TMath::RadToDeg()); + return kTRUE; + } + return kTRUE; +} Index: trunk/MagicSoft/Mars/mtemp/mmpi/MPointingPosCalcWA.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mmpi/MPointingPosCalcWA.h (revision 5490) +++ trunk/MagicSoft/Mars/mtemp/mmpi/MPointingPosCalcWA.h (revision 5490) @@ -0,0 +1,44 @@ +#ifndef MARS_MPointingPositionCalcWA +#define MARS_MPointingPositionCalcWA + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +class MMcEvt; +class MReportDrive; +class MPointingPos; +class MRawRunHeader; +class MObservatory; + +class MPointingPosCalcWA : public MTask +{ +private: + MMcEvt *fMcEvt; //! MMcEvt to get simulated poiting position from + MReportDrive *fReport; //! MReportDrive to get real poiting position from + MPointingPos *fPosition; //! Output container to store pointing position + MObservatory *fObs; + MRawRunHeader *fRunHeader; + Double_t fRa; + Double_t fDec; + + + UShort_t fRunType; //! Run Type to decide where to get pointing position from + + Bool_t ReInit(MParList *plist); + Int_t PreProcess(MParList *plist); + Int_t Process(); + +public: + MPointingPosCalcWA() + { + fName = "MPointingPosCalcWA"; + fTitle = "Task calculating the pointing position"; + } + + void SetRaDec(Double_t ra, Double_t dec) {fRa = ra; fDec = dec;} + + ClassDef(MPointingPosCalcWA, 0) //Task calculating the pointing position +}; + +#endif Index: trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.cc (revision 5489) +++ trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.cc (revision 5490) @@ -123,4 +123,5 @@ #include #include +#include #include #include @@ -281,6 +282,10 @@ } + kSaveAlphaPlots=kTRUE; + kSaveSkyPlots=kTRUE; + kSaveNexPlot=kTRUE; fAlphaHName = "alphaplot.root"; fSkyHName = "skyplot.root"; + fNexHName = "Nexcess.gif"; } @@ -352,4 +357,8 @@ fMinYGrid = ymin; fMaxYGrid = ymax; + + *fLog << endl << inf << " SetSkyPlot: fMinXGrid, fMaxXGrid, fMinYGrid, fMaxYGrid, fBinStepGrid: " << endl; + *fLog << inf << " " << fMinXGrid << ", " << fMaxXGrid << ", " << fMinYGrid << ", " + << fMaxYGrid<< ", " << fBinStepGrid << endl; } @@ -473,6 +482,4 @@ // for the current pointing position and add a offset in the // Fill function! - kSaveAlphaPlots=kTRUE; - kSaveSkyPlots=kTRUE; // prepare skyplot @@ -480,4 +487,7 @@ fNumStepsY = (int) ((fMaxYGrid - fMinYGrid) / fBinStepGrid + 1.5); fNumalphahist = (int) (fNumStepsX * fNumStepsY + 0.5); + + *fLog << inf << "SetSkyPlot: fNumStepsX, fNumStepsY, fNumalphahist: " + << fNumStepsX << ", " << fNumStepsY << ", " << fNumalphahist << endl; // fHistSignif.SetName("SPSignif2ndOrder"); @@ -519,4 +529,5 @@ fHistAlpha[i].SetDirectory(NULL); } +//cout << " fHistAlpha[10].GetBinContent(5) " << fHistAlpha[10].GetBinContent(5) << endl; delete temp; @@ -626,9 +637,10 @@ */ +// make the center of the plot different from the center of the camera /* x_0 = fPntPosCam->GetX()*fMm2Deg; y_0 = fPntPosCam->GetY()*fMm2Deg; */ - x_0 = 0.; + x_0 = 0.; y_0 = 0.; @@ -641,5 +653,7 @@ for (Int_t gridx = 0; gridx < fNumStepsX; gridx++) { + xsource = fMinXGrid + gridx * fBinStepGrid; + /* derotation : rotate into camera coordinates */ /* formel: (x_cam) (cos(gam) -sin(gam)) (xtilde) (x_0) @@ -649,5 +663,9 @@ xsourcam = cosgam * xsource - singam * ysource + x_0; ysourcam = singam * xsource + cosgam * ysource + y_0; + + /* end derotatiom */ +// xsourcam = xsource; +// ysourcam = ysource; /* calculate ALPHA und DIST according to the source position */ @@ -698,6 +716,4 @@ const Int_t onbinsalpha = TMath::Nint(fAlphaONMax/fHistAlphaBinWidth); -// fit with gaus - fHistSignifGaus.Fit("gaus"); // fit function for the background @@ -713,4 +729,9 @@ alpha_iterator = fHistAlpha.begin(); + fHistNexcess.Reset(); + fHistOn.Reset(); + fHistSignif.Reset(); + fHistSignifGaus.Reset(); + while( alpha_iterator != fHistAlpha.end() ) { // find the bin in the 2dim histogram @@ -751,4 +772,8 @@ } +// fit with gaus + fHistSignifGaus.Fit("gaus","N"); + + //temp /* @@ -770,4 +795,7 @@ // TString stri2 = "skyplots.root"; if(kSaveSkyPlots==kTRUE) SaveSkyPlots(fSkyHName); + +// save nex plot: + if(kSaveNexPlot==kTRUE) SaveNexPlot(fNexHName); fHistAlpha.clear(); @@ -800,4 +828,40 @@ } +void MSkyPlot::SaveNexPlot(TString nexplotname) +{ + gStyle->SetCanvasBorderMode(0); + gStyle->SetCanvasBorderSize(0); + gStyle->SetCanvasColor(10); +// gStyle->SetPadBorderMode(0); +// gStyle->SetPadBorderSize(0); + gStyle->SetPadColor(10); + gStyle->SetOptFit(0); + gStyle->SetOptStat(0); + gStyle->SetStatColor(10); + gStyle->SetPalette(1); + gStyle->SetPadRightMargin(0.16); + gStyle->SetPadLeftMargin(0.13); + + TH2D tmp = fHistNexcess; + tmp.SetMaximum(470); + tmp.SetMinimum(-90); + TCanvas can("SkyPlot","SkyPlot", 0, 0, 800, 400); + can.Divide(2,1); + can.cd(1); + tmp.GetYaxis()->SetTitleOffset(1.3); + tmp.Draw("colz"); + TPaveLabel myname(fMinXGrid-0.5,fMinYGrid-0.4,fMinXGrid+0.3,fMinYGrid-0.2,"by D. Mazin"); + myname.Draw(); + + can.cd(2); + fHistNexcess.SetMaximum(470); + fHistNexcess.SetMinimum(-40); + fHistNexcess.GetXaxis()->SetTitleOffset(1.3); + fHistNexcess.GetYaxis()->SetTitleOffset(1.6); + gPad->SetTheta(20); + fHistNexcess.Draw("lego2"); + can.Print(nexplotname); +// can.Print("test.root"); +} void MSkyPlot::SaveSkyPlots(TString skyplotfilename) Index: trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.h (revision 5489) +++ trunk/MagicSoft/Mars/mtemp/mmpi/MSkyPlot.h (revision 5490) @@ -80,4 +80,5 @@ Bool_t kSaveAlphaPlots; Bool_t kSaveSkyPlots; + Bool_t kSaveNexPlot; Float_t fMinXGrid; // [deg] , left edge of the skyplot @@ -105,4 +106,5 @@ TString fAlphaHName; // name for histogram with alpha plots TString fSkyHName; // name for histogram with sky plots + TString fNexHName; // name for canvas with Nex plot Int_t DistancetoPrimitive(Int_t px, Int_t py); @@ -128,18 +130,21 @@ void SetSizeMin(Float_t size) { fSizeMin = size; } // Absolute minimum Size void SetSizeMax(Float_t size) { fSizeMax = size; } // Absolute maximum Size -// void SetMinLD(Float_t ratio) { fMinLD = ratio; } // Minimum ratio between length/dist -// void SetMaxLD(Float_t ratio) { fMaxLD = ratio; } // Maximum ratio between length/dist void ReadCuts(const TString parSCinit); void SetSkyPlot(Float_t xmin, Float_t xmax, Float_t ymin, Float_t ymax, Float_t step); Double_t CalcLimit(Double_t *a, Double_t ls, Double_t ls2, Double_t dd2); - //void SaveSkyPlots(const TString skyplotfilename="skyplots.root"); + + void SetOutputSkyName(TString outname2) { fSkyHName = outname2; } + void SaveSkyPlots(TString stri); + void SetNotSaveSkyPlots() { kSaveSkyPlots = kFALSE; } + void SetOutputAlphaName(TString outname1) { fAlphaHName = outname1; } - void SaveSkyPlots(TString stri); + void SaveAlphaPlots(const TString stri2); + void SetNotSaveAlphaPlots() { kSaveAlphaPlots = kFALSE; } + void SetOutputNexName(TString outname3) { fNexHName = outname3; } + void SaveNexPlot(const TString stri3); + void SetNotSaveNexPlot() { kSaveNexPlot = kFALSE; } - //void SaveAlphaPlots(const TString alphaplotfilename="alphaplots.root"); - void SetOutputSkyName(TString outname2) { fSkyHName = outname2; } - void SaveAlphaPlots(const TString stri2); void SetAlphaCut(Float_t alpha);