Index: /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MFindStars.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MFindStars.cc (revision 5120) +++ /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MFindStars.cc (revision 5120) @@ -0,0 +1,1274 @@ +/* ======================================================================== *\ +! +! * +! * 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, 2/2004 +! Javier López , 4/2004 +! Wolfgang Wittek, 8/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MFindStars +// +///////////////////////////////////////////////////////////////////////////// +#include "MFindStars.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +#include "MObservatory.h" +#include "MAstroCamera.h" +#include "MMcConfigRunHeader.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MHCamera.h" + +#include "MGeomCam.h" +#include "MGeomPix.h" +#include "MCameraDC.h" +#include "MTime.h" +#include "MReportDrive.h" +#include "MStarCam.h" +#include "MStarPos.h" + +#include "MParList.h" +#include "MTaskList.h" + +ClassImp(MFindStars); +using namespace std; + +const Float_t sqrt2 = sqrt(2.); +const Float_t sqrt3 = sqrt(3.); +const UInt_t lastInnerPixel = 396; + + +//______________________________________________________________________________ +// +// The 2D uncorrelated gaussian function used to fit the spot of the star +// +static Double_t func(float x,float y,Double_t *par) +{ + Double_t value=par[0]*exp( -(x-par[1])*(x-par[1])/(2*par[2]*par[2]) ) + *exp( -(y-par[3])*(y-par[3])/(2*par[4]*par[4]) ); + return value; +} + +//______________________________________________________________________________ +// +// The 2D correlated gaussian function used to fit the spot of the star +// +static Double_t funcCG(float x,float y,Double_t *par) +{ + Double_t temp = 1.0-par[5]*par[5]; + // Double_t value = par[0] / (2.0*TMath::Pi()*par[2]*par[4]*sqrt(temp)) + Double_t value = par[0] + * exp( -0.5/temp * ( (x-par[1])*(x-par[1])/(par[2]*par[2]) + -2.0*par[5] * (x-par[1])*(y-par[3])/(par[2]*par[4]) + + (y-par[3])*(y-par[3])/(par[4]*par[4]) ) ); + return value; +} + +//______________________________________________________________________________ +// +// Function used by Minuit to do the fit +// +static void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) +{ + + MParList* plist = (MParList*)gMinuit->GetObjectFit(); + MTaskList* tlist = (MTaskList*)plist->FindObject("MTaskList"); + MFindStars* find = (MFindStars*)tlist->FindObject("MFindStars"); + MStarCam* stars = + (MStarCam*)plist->FindObject("MStarCam","MStarCam"); + MGeomCam* geom = (MGeomCam*)plist->FindObject("MGeomCam"); + + MHCamera& display = (MHCamera&)find->GetDisplay(); + + Float_t innerped = stars->GetInnerPedestalDC(); + Float_t innerrms = stars->GetInnerPedestalRMSDC(); + Float_t outerped = stars->GetOuterPedestalDC(); + Float_t outerrms = stars->GetOuterPedestalRMSDC(); + + UInt_t numPixels = geom->GetNumPixels(); + + +//calculate chisquare + Double_t chisq = 0; + Double_t delta; + Double_t x,y,z; + Double_t errorz=0; + + UInt_t usedPx=0; + for (UInt_t pixid=1; pixidlastInnerPixel?outerped:innerped); + errorz=(pixid>lastInnerPixel?outerrms:innerrms); + + + if (errorz > 0.0) + { + usedPx++; + + Double_t fu; + if (find->GetUseCorrelatedGauss()) + fu = funcCG(x,y,par); + else + fu = func(x,y,par); + + delta = (z-fu) / errorz; + chisq += delta*delta; + + if (iflag == 3) + { + gLog << "fcn : usedPx, pixid, content, pedestal,z, fu, errorz, delta = " + << usedPx << ", " << pixid << ", " + << display.GetBinContent(pixid+1) << ", " + << (pixid>lastInnerPixel?outerped:innerped) + << ", " << z << ", " << fu << ", " << errorz + << ", " << delta << endl; + } + } + else + cerr << " TMinuit::fcn errorz[" << pixid << "] " << errorz << endl; + } + } + f = chisq; + + find->SetChisquare(chisq); + find->SetDegreesofFreedom(usedPx); + + //gLog << "fcn : chisq, usedPx = " << chisq << ", " << usedPx << endl; +} + +//------------------------------------------------------------------------- +// +// Constructor +// + +MFindStars::MFindStars(const char *name, const char *title): + fGeomCam(NULL), fCurr(NULL), fTimeCurr(NULL), fDrive(NULL), fStars(NULL), + fNumVar(6) +{ + fName = name ? name : "MFindStars"; + fTitle = title ? title : "Tool to find stars from DC Currents"; + + // the correlated Gauss function + // is fitted by default + fNumVar = 6; + fUseCorrelatedGauss = kTRUE; + + fNumIntegratedEvents=0; + fMaxNumIntegratedEvents = 10; + fRingInterest = 125.; //[mm] ~ 0.4 deg + fDCTailCut = 4; + + fPixelsUsed.Set(577); + fPixelsUsed.Reset((Char_t)kTRUE); + + //Fitting(Minuit) initialitation + const Float_t pixelSize = 31.5; //[mm] + fMinuitPrintOutLevel = -1; + //fMinuitPrintOutLevel = 3; + + fVname = new TString[fNumVar]; + fVinit.Set(fNumVar); + fStep.Set(fNumVar); + fLimlo.Set(fNumVar); + fLimup.Set(fNumVar); + fFix.Set(fNumVar); + + fVname[0] = "max"; + fVinit[0] = 10.*fMaxNumIntegratedEvents; + fStep[0] = fVinit[0]/sqrt2; + fLimlo[0] = fMinDCForStars; + fLimup[0] = 30.*fMaxNumIntegratedEvents; + fFix[0] = 0; + + fVname[1] = "meanx"; + fVinit[1] = 0.; + fStep[1] = fVinit[1]/sqrt2; + fLimlo[1] = -600.; + fLimup[1] = 600.; + fFix[1] = 0; + + fVname[2] = "sigmax"; + fVinit[2] = pixelSize; + fStep[2] = fVinit[2]/sqrt2; + fLimlo[2] = pixelSize/(2*sqrt3); + fLimup[2] = 500.; + fFix[2] = 0; + + fVname[3] = "meany"; + fVinit[3] = 0.; + fStep[3] = fVinit[3]/sqrt2; + fLimlo[3] = -600.; + fLimup[3] = 600.; + fFix[3] = 0; + + fVname[4] = "sigmay"; + fVinit[4] = pixelSize; + fStep[4] = fVinit[4]/sqrt2; + fLimlo[4] = pixelSize/(2*sqrt3); + fLimup[4] = 500.; + fFix[4] = 0; + + if (fUseCorrelatedGauss) + { + fVname[5] = "xycorr"; + fVinit[5] = 0.0; + fStep[5] = 0.1; + fLimlo[5] = -1.0; + fLimup[5] = 1.0; + fFix[5] = 0; + } + + fObjectFit = NULL; + // fMethod = "SIMPLEX"; + fMethod = "MIGRAD"; + // fMethod = "MINIMIZE"; + fNulloutput = kFALSE; + + // Set output level + // fLog->SetOutputLevel(3); // No dbg messages + + fGeometryFile=""; + fBSCFile=""; +} + +//------------------------------------------------------------------------- +// +// Set 'fUseCorrelatedGauss' flag +// +// if 'fUseCorrelatedGauss' is TRUE a 2dim Gaussian with correlation +// will be fitted +// + +void MFindStars::SetUseCorrelatedGauss(Bool_t usecorrgauss) +{ + fUseCorrelatedGauss = usecorrgauss; + + if (usecorrgauss) + fNumVar = 6; + else + fNumVar = 5; +} + +//------------------------------------------------------------------------- +// +// PreProcess +// + +Int_t MFindStars::PreProcess(MParList *pList) +{ + fGeomCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam")); + + if (!fGeomCam) + { + *fLog << err << AddSerialNumber("MGeomCam") << " not found ... aborting" << endl; + return kFALSE; + } + + // Initialize camera display with the MGeomCam information + fDisplay.SetGeometry(*fGeomCam,"FindStarsDisplay","FindStarsDisplay"); + fDisplay.SetUsed(fPixelsUsed); + + fCurr = (MCameraDC*)pList->FindObject(AddSerialNumber("MCameraDC")); + + if (!fCurr) + { + *fLog << err << AddSerialNumber("MCameraDC") << " not found ... aborting" << endl; + return kFALSE; + } + + fTimeCurr = (MTime*)pList->FindObject(AddSerialNumber("MTimeCurrents")); + + if (!fTimeCurr) + { + *fLog << err << AddSerialNumber("MTimeCurrents") << " not found ... aborting" << endl; + return kFALSE; + } + + fDrive = (MReportDrive*)pList->FindObject(AddSerialNumber("MReportDrive")); + + if (!fDrive) + { + + *fLog << warn << AddSerialNumber("MReportDrive") << " not found ... ignored." << endl; + + } + else + { + MObservatory magic1; + + MMcConfigRunHeader *config=0; + MGeomCam *geom=0; + + TFile file(fGeometryFile); + TTree *tree = dynamic_cast(file.Get("RunHeaders")); + tree->SetBranchAddress("MMcConfigRunHeader.", &config); + if (tree->GetBranch("MGeomCam.")) tree->SetBranchAddress("MGeomCam.", &geom); + + TClonesArray mirrorArray = *config->GetMirrors(); + fAstro.SetMirrors(*config->GetMirrors()); + fAstro.SetGeom(*geom); + fAstro.ReadBSC(fBSCFile); + fAstro.SetObservatory(magic1); + + } + + + fStars = (MStarCam*)pList->FindCreateObj(AddSerialNumber("MStarCam"),"MStarCam"); + if (!fStars) + { + *fLog << err << AddSerialNumber("MStarCam") << " cannot be created ... aborting" << endl; + return kFALSE; + } + + fMinDCForStars = 1.*fMaxNumIntegratedEvents; //[uA] + + // Initialize the TMinuit object + + TMinuit *gMinuit = new TMinuit(fNumVar); //initialize TMinuit with a maximum of params + gMinuit->SetFCN(fcn); + + Double_t arglist[10]; + Int_t ierflg = 0; + + arglist[0] = 1; + gMinuit->mnexcm("SET ERR", arglist ,1,ierflg); + arglist[0] = fMinuitPrintOutLevel; + gMinuit->mnexcm("SET PRI", arglist ,1,ierflg); + + // suppress warnings + Int_t errWarn; + Double_t tmpwar = 0; + gMinuit->mnexcm("SET NOW", &tmpwar, 0, errWarn); + + + // Set MParList object pointer to allow mimuit to access internally + gMinuit->SetObjectFit(pList); + return kTRUE; +} + +//------------------------------------------------------------------------ +// +// Process : +// +// +// +// + +Int_t MFindStars::Process() +{ + UInt_t numPixels = fGeomCam->GetNumPixels(); + TArrayC origPixelsUsed; + origPixelsUsed.Set(numPixels); + + if (fNumIntegratedEvents >= fMaxNumIntegratedEvents) { + + //Fist delete the previous stars in the list + fStars->GetList()->Delete(); + + if (fDrive) { + + //If drive information is provided we take RaDec info + //from the drive and let the star list fill by the astrocamera. + + //FIXME: rwagner: Doesn't work as expected + //FIXME: rwagner: For the time being set manually. + //fAstro.SetRaDec(fDrive->GetRa(), fDrive->GetDec()); + + fAstro.SetTime(*fTimeCurr); + fAstro.SetGuiActive(); + fAstro.FillStarList(fStars->GetList()); + + cout << "Number of Stars added by astrocamera is " + <GetList()->GetSize() << endl; + + MStarPos* starpos; + TIter Next(fStars->GetList()); + while ((starpos=(MStarPos*)Next())) { + //starpos->SetCalcValues(40,40,starpos->GetXExp(), starpos->GetYExp(), + // fRingInterest/2, fRingInterest/2, + // 0., 0., 0., 0.); + //starpos->SetFitValues (40,40,starpos->GetXExp(), starpos->GetYExp(), + // fRingInterest/2, fRingInterest/2, 0., + // 0., 0., 0., 0., -1); + } + + for (UInt_t pix=1; pixGetInnerPedestalDC()+fDCTailCut*fStars->GetInnerPedestalRMSDC(); + Float_t outermin = fStars->GetOuterPedestalDC()+fDCTailCut*fStars->GetOuterPedestalRMSDC(); + fMinDCForStars = innermin>outermin?innermin:outermin; + if (fMinuitPrintOutLevel>=0) *fLog << dbg << "fMinDCForStars = " << fMinDCForStars << endl; + + // Find the star candidates searching the most brights pairs of pixels + Float_t maxPixelDC; + MGeomPix maxPixel; + + while(FindPixelWithMaxDC(maxPixelDC, maxPixel)) { + MStarPos *starpos = new MStarPos; + + starpos->SetExpValues(maxPixelDC,maxPixel.GetX(),maxPixel.GetY()); + + //starpos->SetCalcValues(maxPixelDC, maxPixelDC, + // maxPixel.GetX(), maxPixel.GetY(), + // fRingInterest/2, fRingInterest/2, + // 0., 0., 0., 0.); + + //starpos->SetFitValues(maxPixelDC, maxPixelDC, + // maxPixel.GetX(), maxPixel.GetY(), + // fRingInterest/2, fRingInterest/2, 0., + // 0., 0., 0., 0., -1); + + fStars->GetList()->Add(starpos); + ShadowStar(starpos); + } + fDisplay.SetUsed(origPixelsUsed); + } + } + + //Show the stars found + fStars->Print("namepossizechierr"); + + //loop to extract position of stars on the camera + if (fStars->GetList()->GetSize() == 0) { + *fLog << err << GetName() << "No stars candidates in the camera." << endl; + return kCONTINUE; + } + else + { + *fLog << inf << GetName() << " found " << fStars->GetList()->GetSize() + << " star candidates in the camera." << endl; + } + + for (UInt_t pix=1; pixGetList()); + MStarPos* star; + while ((star=(MStarPos*)Next())) { + FindStar(star); + ShadowStar(star); + } + + //Show the stars found: Here it is interesting to see what FindStars + //made out of the positions we gave to it. + if (fMinuitPrintOutLevel>=0) + fStars->Print("namepossizchierr"); + + //After finding stars reset all variables + fDisplay.Reset(); + fStars->GetDisplay().Reset(); //FIXME: Put this display just in the container + fDisplay.SetUsed(origPixelsUsed); + fNumIntegratedEvents=0; + } + + for (UInt_t pix=1; pix=0) *fLog << dbg << "MFindStars::SetBlindPixels fDisplay.IsUsed(" <GetNumPixels(); + Float_t ped; + Float_t rms; + + //TH1F **dchist = new TH1F*[2]; + TH1F *dchist[2]; + TString tit; + TString nam; + for (UInt_t i=0; i<2; i++) + { + nam = i; + tit = "dchist["; + tit += i; + tit += "]"; + dchist[i] = new TH1F(nam, tit, + 26,0.4*fMaxNumIntegratedEvents,3.*fMaxNumIntegratedEvents); + } + + for (UInt_t pix=1; pix<=lastInnerPixel; pix++) + dchist[0]->Fill(fDisplay.GetBinContent(pix+1)); + for (UInt_t pix=lastInnerPixel+1; pixFill(fDisplay.GetBinContent(pix+1)); + + // inner/outer pixels + for (UInt_t i=0; i<2; i++) + { + Float_t nummaxprobdc = dchist[i]->GetBinContent(dchist[i]->GetMaximumBin()); + Float_t maxprobdc = dchist[i]->GetBinCenter(dchist[i]->GetMaximumBin()); + UInt_t bin = dchist[i]->GetMaximumBin(); + do + { + bin++; + } + while(dchist[i]->GetBinContent(bin)/nummaxprobdc > 0.5); + Float_t halfmaxprobdc = dchist[i]->GetBinCenter(bin); + + if (fMinuitPrintOutLevel>=0) *fLog << dbg << " maxprobdc[high] " << maxprobdc << "[" << nummaxprobdc << "] "; + if (fMinuitPrintOutLevel>=0) *fLog << dbg << " halfmaxprobdc[high] " << halfmaxprobdc << "[" << dchist[i]->GetBinContent(bin) << "]" << endl; + + Float_t rmsguess = TMath::Abs(maxprobdc-halfmaxprobdc); + Float_t min = maxprobdc-3*rmsguess; + min = (min<0.?0.:min); + Float_t max = maxprobdc+3*rmsguess; + + if (fMinuitPrintOutLevel>=0) *fLog << dbg << " maxprobdc " << maxprobdc << " rmsguess " << rmsguess << endl; + + TF1 func("func","gaus",min,max); + func.SetParameters(nummaxprobdc, maxprobdc, rmsguess); + + dchist[i]->Fit("func","QR0"); + + UInt_t aproxnumdegrees = 6*(bin-dchist[i]->GetMaximumBin()); + Float_t chiq = func.GetChisquare(); + ped = func.GetParameter(1); + rms = func.GetParameter(2); + + if (fMinuitPrintOutLevel>=0) *fLog << dbg << " ped " << ped << " rms " << rms << " chiq/ndof " << chiq << "/" << aproxnumdegrees << endl; + + Int_t numsigmas = 5; + Axis_t minbin = ped-numsigmas*rms/dchist[i]->GetBinWidth(1); + minbin=minbin<1?1:minbin; + Axis_t maxbin = ped+numsigmas*rms/dchist[i]->GetBinWidth(1); + if (fMinuitPrintOutLevel>=0) *fLog << dbg << " Number of pixels with dc under " << numsigmas << " sigmas = " << dchist[i]->Integral((int)minbin,(int)maxbin) << endl; + + //Check results from the fit are consistent + if (ped < 0. || rms < 0.) + { + *fLog << dbg << "dchist[i]->GetEntries()" << dchist[i]->GetEntries(); +// TCanvas *c1 = new TCanvas("c2","c2",500,800); +// dchist[i]->Draw(); +// c1->Print("dchist.ps"); +// delete c1; +// exit(1); + } + else if (TMath::Abs(ped-maxprobdc) > rmsguess || rms > rmsguess) + { + *fLog << warn << GetName() << " Pedestal DC fit give non consistent results for " << (i==0?"Inner":"Outer") << "pixels." << endl; + *fLog << warn << " maxprobdc " << maxprobdc << " rmsguess " << rmsguess << endl; + *fLog << warn << " ped " << ped << " rms " << rms << " chiq/ndof " << chiq << "/" << aproxnumdegrees << endl; + ped = maxprobdc; + rms = rmsguess/1.175; // FWHM=2.35*rms + } + + if (i == 0) + { + fStars->SetInnerPedestalDC(ped); + fStars->SetInnerPedestalRMSDC(rms); + } + else + { + fStars->SetOuterPedestalDC(ped); + fStars->SetOuterPedestalRMSDC(rms); + } + } + + + + for (UInt_t i=0; i<2; i++) + { + delete dchist[i]; + } + //delete [] dchist; + + //================================================================= + + // reset gMinuit to the MINUIT object for optimizing the supercuts + gMinuit = savePointer; + //------------------------------------------- + + if (fStars->GetInnerPedestalDC() < 0. || fStars->GetInnerPedestalRMSDC() < 0. || fStars->GetOuterPedestalDC() < 0. || fStars->GetOuterPedestalRMSDC() < 0.) + return kFALSE; + + return kTRUE; +} + +Bool_t MFindStars::FindPixelWithMaxDC(Float_t &maxDC, MGeomPix &maxPix) +{ + UInt_t numPixels = fGeomCam->GetNumPixels(); + +// Find the two close pixels with the maximun dc + UInt_t maxPixIdx[2]; + + maxDC = 0; + + for (UInt_t pix=1; pix maxDC*2) + { + if(dc[0]>=dc[1]) + { + maxPixIdx[0] = pix; + maxPixIdx[1] = swneighbor; + maxDC = dc[0]; + } + else + { + maxPixIdx[1] = pix; + maxPixIdx[0] = swneighbor; + maxDC = dc[1]; + } + } + } + } + } + } + + if (maxDC == 0) + { + *fLog << warn << " No found pixels with maximum dc" << endl; + return kFALSE; + } + + maxPix = (*fGeomCam)[maxPixIdx[0]]; + + if (fMinuitPrintOutLevel>=0) *fLog << dbg << "Star candidate maxDC(" << setw(3) << maxDC << " uA) x position(" << setw(3) << maxPix.GetX() << " mm) x position(" << setw(3) << maxPix.GetY() << " mm) swnumber(" << maxPixIdx[0] << ")" << endl; + + return kTRUE; +} + + +//---------------------------------------------------------------------------- +// +// FindStar +// +// +// + +Bool_t MFindStars::FindStar(MStarPos* star) +{ + UInt_t numPixels = fGeomCam->GetNumPixels(); + Float_t innerped = fStars->GetInnerPedestalDC(); + Float_t outerped = fStars->GetOuterPedestalDC(); + + TArrayC origPixelsUsed; + origPixelsUsed.Set(numPixels); + + for (UInt_t pix=1; pixGetXExp(); + Float_t expY = star->GetYExp(); + + Float_t max =0; + UInt_t pixmax =0; + Float_t meanX =0; + Float_t meanY =0; + Float_t meanSqX =0; + Float_t meanSqY =0; + Float_t sumCharge =0; + Float_t sumSqCharge =0; + UInt_t usedInnerPx =0; + UInt_t usedOuterPx =0; + Float_t factor = 1.0; + + Float_t meanXold = 1.e10; + Float_t meanYold = 1.e10; + + const Float_t meanPresition = 3.; //[mm] + const UInt_t maxNumIterations = 10; + UInt_t numIterations = 0; + + //-------------------- start of iteration loop ----------------------- + for (UInt_t it=0; itmax) + { + max=charge; + pixmax=pix; + } + } + + // 3.) make it new center + + expX = (*fGeomCam)[pixmax].GetX(); + expY = (*fGeomCam)[pixmax].GetY(); + + for (UInt_t pix=1; pixlastInnerPixel?usedOuterPx++:usedInnerPx++; + + Float_t charge = fDisplay.GetBinContent(pix+1); + Float_t pixXpos = (*fGeomCam)[pix].GetX(); + Float_t pixYpos = (*fGeomCam)[pix].GetY(); + + if (charge>max) + { + max=charge; + pixmax=pix; + } + + meanX += charge*pixXpos; + meanY += charge*pixYpos; + meanSqX += charge*pixXpos*pixXpos; + meanSqY += charge*pixYpos*pixYpos; + sumCharge += charge; + sumSqCharge += charge*charge; + } + } + + if (fMinuitPrintOutLevel>=0) *fLog << dbg << " usedInnerPx " << usedInnerPx << " usedOuterPx " << usedOuterPx << endl; + + meanX /= sumCharge; + meanY /= sumCharge; + meanSqX /= sumCharge; + meanSqY /= sumCharge; + factor = sqrt( sumSqCharge/(sumCharge*sumCharge) ); + + // stop iteration when (meanX, meanY) becomes stable + if ( TMath::Abs(meanX-meanXold) < meanPresition && + TMath::Abs(meanY-meanYold) < meanPresition ) + break; + + meanXold = meanX; + meanYold = meanY; + numIterations++; + } + // this statement was commented out because the condition + // was never fulfilled : + //while(TMath::Abs(meanX-expX) > meanPresition || + // TMath::Abs(meanY-expY) > meanPresition); + //-------------------- end of iteration loop ----------------------- + + if (numIterations >= maxNumIterations) + { + *fLog << warn << GetName() << "Mean calculation not converged after " + << maxNumIterations << " iterations" << endl; + } + + + //Float_t rmsX = (meanSqX - meanX*meanX) - fRingInterest*fRingInterest/12; + //Float_t rmsY = (meanSqY - meanY*meanY) - fRingInterest*fRingInterest/12; + + Float_t rmsX = (meanSqX - meanX*meanX); + Float_t rmsY = (meanSqY - meanY*meanY); + + if ( rmsX > 0) + rmsX = TMath::Sqrt(rmsX); + else + { + *fLog << warn << " MFindStars::FindStar negative rmsX² " << rmsX << endl; + *fLog << warn << " meanSqX " << meanSqX << " meanX " << meanX << " fRingInterest " << fRingInterest << " sumCharge " << sumCharge << endl; + rmsX = 0.; + } + + if ( rmsY > 0) + rmsY = TMath::Sqrt(rmsY); + else + { + *fLog << warn << " MFindStars::FindStar negative rmsY² " << rmsY << endl; + *fLog << warn << " meanSqY " << meanSqY << " meanY " << meanY << " fRingInterest " << fRingInterest << " sumCharge " << sumCharge<< endl; + rmsY = 0.; + } + + Float_t deltameanX = factor * rmsX; + Float_t deltameanY = factor * rmsY; + + + // Substrack pedestal DC + sumCharge-= (usedInnerPx*innerped+usedOuterPx*outerped)/(usedInnerPx+usedOuterPx); + max-=pixmax>lastInnerPixel?outerped:innerped; + + + star->SetCalcValues( sumCharge, max, meanX, meanY, rmsX, rmsY, + 0., deltameanX*deltameanX, 0., deltameanY*deltameanY); + + if (rmsX <= 0. || rmsY <= 0.) + return kFALSE; + + +// fit the star spot using TMinuit + + + for (UInt_t pix=1; pix=0) *fLog << dbg << "[fit the star spot] fDisplay.IsUsed(" << pix << ") kTRUE" << endl; + + if (fMinuitPrintOutLevel>=0) *fLog << dbg << "fMinDCForStars " << fMinDCForStars << " pixmax>lastInnerPixel?outerped:innerped " << (pixmax>lastInnerPixel?outerped:innerped) << " fMaxNumIntegratedEvents " << fMaxNumIntegratedEvents << endl; + + //Initialize variables for fit + fVinit[0] = max; + fLimlo[0] = fMinDCForStars-(pixmax>lastInnerPixel?outerped:innerped); + fLimup[0] = 30*fMaxNumIntegratedEvents-(pixmax>lastInnerPixel?outerped:innerped); + fVinit[1] = meanX; + fVinit[2] = rmsX; + fVinit[3] = meanY; + fVinit[4] = rmsY; + //Init steps + for(Int_t i=0; i=0) *fLog << dbg << " fVinit[" << i << "] " << fVinit[i]; + if (fMinuitPrintOutLevel>=0) *fLog << dbg << " fStep[" << i << "] " << fStep[i]; + if (fMinuitPrintOutLevel>=0) *fLog << dbg << " fLimlo[" << i << "] " << fLimlo[i]; + if (fMinuitPrintOutLevel>=0) *fLog << dbg << " fLimup[" << i << "] " << fLimup[i] << endl; + } + + + // ------------------------------------------- + // call MINUIT + + Double_t arglist[10]; + Int_t ierflg = 0; + + for (Int_t i=0; imnparm(i, fVname[i], fVinit[i], fStep[i], fLimlo[i], fLimup[i], ierflg); + + TStopwatch clock; + clock.Start(); + +// Now ready for minimization step + arglist[0] = 500; + arglist[1] = 1.; + gMinuit->mnexcm(fMethod, arglist ,2,ierflg); + + // call fcn with iflag = 3 + //arglist[0] = 3; + //Int_t ierflg3; + //gMinuit->mnexcm("CALL", arglist, 1, ierflg3); + + clock.Stop(); + + if(fMinuitPrintOutLevel>=0) + { + if (fMinuitPrintOutLevel>=0) *fLog << dbg << "Time spent for the minimization in MINUIT : " << endl;; + clock.Print(); + } + + + //---------- for the uncorrelated Gauss fit (start) --------- + if (!fUseCorrelatedGauss) + { + Double_t integratedCharge; + Double_t maxFit, maxFitError; + Double_t meanXFit, meanXFitError; + Double_t sigmaX, sigmaXError; + Double_t meanYFit, meanYFitError; + Double_t sigmaY, sigmaYError; + Float_t chisquare = GetChisquare(); + Int_t degrees = GetDegreesofFreedom()-fNumVar; + + if (!ierflg) + { + gMinuit->GetParameter(0, maxFit, maxFitError); + gMinuit->GetParameter(1, meanXFit, meanXFitError); + gMinuit->GetParameter(2, sigmaX, sigmaXError); + gMinuit->GetParameter(3, meanYFit, meanYFitError); + gMinuit->GetParameter(4, sigmaY, sigmaYError); + + //FIXME: Do the integral properlly + integratedCharge = 0.; + } + else + { + maxFit = 0.; + meanXFit = 0.; + sigmaX = 0.; + meanYFit = 0.; + sigmaY = 0.; + integratedCharge = 0.; + + *fLog << err << "TMinuit::Call error " << ierflg << endl; + } + + // get error matrix + Double_t matrix[5][5]; + gMinuit->mnemat(&matrix[0][0],5); + + star->SetFitValues(integratedCharge,maxFit, meanXFit, meanYFit, + sigmaX, sigmaY, 0.0, + matrix[1][1], matrix[1][3], matrix[3][3], + chisquare, degrees); + } + //---------- for the uncorrelated Gauss fit (end) --------- + + //---------- for the correlated Gauss fit (start) --------- + if (fUseCorrelatedGauss) + { + TArrayD fValues(fNumVar); + TArrayD fErrors(fNumVar); + + const static Int_t fNdim = 6; + Double_t matrix[fNdim][fNdim]; + + Double_t fmin = 0.0; + Double_t fedm = 0.0; + Double_t errdef = 0.0; + Int_t npari = 0; + Int_t nparx = 0; + Int_t istat = 0; + gMinuit->mnstat(fmin, fedm, errdef, npari, nparx, istat); + if (fMinuitPrintOutLevel>=0) + { + *fLog << dbg + << "Status of Correlated Gauss fit to the DC currents : " << endl; + } + if (fMinuitPrintOutLevel>=0) + { + *fLog << dbg + << " fmin, fedm, errdef, npari, nparx, istat = " + << fmin << ", " << fedm << ", " << errdef << ", " << npari + << ", " << nparx << ", " << istat << endl; + } + + if (!ierflg) + { + for (Int_t k=0; kGetParameter( k, fValues[k], fErrors[k] ); + } + else + { + *fLog << err << "Correlated Gauss fit failed : ierflg, istat = " + << ierflg << ", " << istat << endl; + } + + Float_t charge = 0.0; + Float_t chi2 = fmin; + Int_t ndof = GetDegreesofFreedom()-fNumVar; + + //get error matrix + if (istat >= 1) + { + gMinuit->mnemat( &matrix[0][0], fNdim); + + // copy covariance matrix into a matrix which includes also the fixed + // parameters + TString name; + Double_t bnd1, bnd2, val, err; + Int_t jvarbl; + Int_t kvarbl; + for (Int_t j=0; jmnpout(j, name, val, err, bnd1, bnd2, jvarbl); + + for (Int_t k=0; kmnpout(k, name, val, err, bnd1, bnd2, kvarbl); + + matrix[j][k] = jvarbl==0 || kvarbl==0 ? 0 : matrix[jvarbl-1][kvarbl-1]; + } + } + } + else + { + // error matrix was not calculated, construct it + if (fMinuitPrintOutLevel>=0) + { + *fLog << warn + << "error matrix is not defined, construct it" << endl; + } + for (Int_t j=0; j=0) + { + *fLog << "Before calling SetFitValues : " << endl; + *fLog << "fValues = " << fValues[0] << ", " << fValues[1] << ", " + << fValues[2] << ", " << fValues[3] << ", " << fValues[4] + << ", " << fValues[5] << endl; + + *fLog << "fErrors = " << fErrors[0] << ", " << fErrors[1] << ", " + << fErrors[2] << ", " << fErrors[3] << ", " << fErrors[4] + << ", " << fErrors[5] << endl; + + *fLog << "error matrix = " << matrix[0][0] << " " << matrix[0][1] + << " " << matrix[0][2] << " " << matrix[0][3] + << " " << matrix[0][4] << " " << matrix[0][5] + << endl; + *fLog << " " << matrix[1][0] << " " << matrix[1][1] + << " " << matrix[1][2] << " " << matrix[1][3] + << " " << matrix[1][4] << " " << matrix[1][5] + << endl; + *fLog << " " << matrix[2][0] << " " << matrix[2][1] + << " " << matrix[2][2] << " " << matrix[2][3] + << " " << matrix[2][4] << " " << matrix[2][5] + << endl; + *fLog << " " << matrix[3][0] << " " << matrix[3][1] + << " " << matrix[3][2] << " " << matrix[3][3] + << " " << matrix[3][4] << " " << matrix[3][5] + << endl; + *fLog << " " << matrix[4][0] << " " << matrix[4][1] + << " " << matrix[4][2] << " " << matrix[4][3] + << " " << matrix[4][4] << " " << matrix[4][5] + << endl; + *fLog << " " << matrix[5][0] << " " << matrix[5][1] + << " " << matrix[5][2] << " " << matrix[5][3] + << " " << matrix[5][4] << " " << matrix[5][5] + << endl; + } + + star->SetFitValues(charge, fValues[0], fValues[1], fValues[3], + fValues[2], fValues[4], fValues[5], + matrix[1][1], matrix[1][3],matrix[3][3], + chi2, ndof); + } + + //---------- for the correlated Gauss fit (end) --------- + + + // reset the display to the starting values + fDisplay.SetUsed(origPixelsUsed); + + if (ierflg) + return kCONTINUE; + return kTRUE; +} + +Bool_t MFindStars::ShadowStar(MStarPos* star) +{ + UInt_t numPixels = fGeomCam->GetNumPixels(); + +// Define an area around the star which will be set unused. + UInt_t shadowPx=0; + for (UInt_t pix=1; pixGetMeanX(); + Float_t starYpos = star->GetMeanY(); + + Float_t starSize = 3*sqrt( star->GetSigmaX()*star->GetSigmaX() + + star->GetSigmaY()*star->GetSigmaY() ); + + Float_t dist = sqrt((pixXpos-starXpos)*(pixXpos-starXpos)+ + (pixYpos-starYpos)*(pixYpos-starYpos)); + + if (dist > starSize && fDisplay.IsUsed(pix)) + fPixelsUsed[pix]=(Char_t)kTRUE; + else + { + fPixelsUsed[pix]=(Char_t)kFALSE; + shadowPx++; + } + } + + if (fMinuitPrintOutLevel>=0) *fLog << dbg << " shadowPx " << shadowPx << endl; + + fDisplay.SetUsed(fPixelsUsed); + + return kTRUE; +} + + + + + + + + + + + + Index: /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MFindStars.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MFindStars.h (revision 5120) +++ /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MFindStars.h (revision 5120) @@ -0,0 +1,123 @@ +#ifndef MARS_MFindStars +#define MARS_MFindStars + +#ifndef ROOT_TArrayS +#include +#endif + +#ifndef ROOT_TArrayC +#include +#endif + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#ifndef MARS_MHCamera +#include "MHCamera.h" +#endif + +#ifndef MARS_MAstroCamera +#include "MAstroCamera.h" +#endif + +class MGeomCam; +class MGeomPix; +class MCameraDC; +class MTime; +class MReportDrive; +class MStarCam; +class MStarPos; + +class MFindStars : public MTask +{ + +private: + + MGeomCam *fGeomCam; + MCameraDC *fCurr; + MTime *fTimeCurr; + MReportDrive *fDrive; + MStarCam *fStars; + + MAstroCamera fAstro; + TArrayC fPixelsUsed; + MHCamera fDisplay; + + UInt_t fMaxNumIntegratedEvents; + UInt_t fNumIntegratedEvents; + + Float_t fRingInterest; //[mm] + Float_t fMinDCForStars; //[uA] + + Float_t fDCTailCut; + + //Fitting(Minuit) variables + Int_t fNumVar; + Float_t fTempChisquare; + Int_t fTempDegreesofFreedom; + Int_t fMinuitPrintOutLevel; + + Bool_t fUseCorrelatedGauss; + + TString *fVname; + TArrayD fVinit; + TArrayD fStep; + TArrayD fLimlo; + TArrayD fLimup; + TArrayI fFix; + TObject *fObjectFit; + TString fMethod; + Bool_t fNulloutput; + + Bool_t DCPedestalCalc(); + Bool_t FindPixelWithMaxDC(Float_t &maxDC, MGeomPix &maxPix); + Bool_t FindStar(MStarPos* star); + Bool_t ShadowStar(MStarPos* star); + + TString fGeometryFile; + TString fBSCFile; + + public: + + MFindStars(const char *name=NULL, const char *title=NULL); + + Int_t PreProcess(MParList *pList); + Int_t Process(); + Int_t PostProcess(); + + void SetUseCorrelatedGauss(Bool_t usecorrgauss = kTRUE); + Bool_t GetUseCorrelatedGauss() {return fUseCorrelatedGauss;} + + // setters + void SetNumIntegratedEvents(UInt_t max) {fMaxNumIntegratedEvents=max;} + void SetRingInterest(Float_t ring) {fRingInterest=ring;} + void SetBlindPixels(TArrayS blindpixels); + void SetMinuitPrintOutLevel(Int_t level) {fMinuitPrintOutLevel=level;} + void SetDCTailCut(Float_t cut) {fDCTailCut=cut;} + + void SetChisquare(Float_t chi) {fTempChisquare=chi;} + void SetDegreesofFreedom(Int_t free) {fTempDegreesofFreedom=free;} + + void SetGeometryFile(TString f) {fGeometryFile=f;} + void SetBSCFile(TString f) {fBSCFile=f;} + void SetRaDec(Double_t ra, Double_t dec) {fAstro.SetRaDec(ra,dec);} + void SetRaDec(TVector3 &v) { fAstro.SetRaDec(v); } + void SetLimMag(Double_t mag) { fAstro.SetLimMag(mag); } + void SetRadiusFOV(Double_t deg) { fAstro.SetRadiusFOV(deg); } + + + //Getters + MHCamera& GetDisplay() { return fDisplay; } + + Float_t GetChisquare() {return fTempChisquare;} + Int_t GetDegreesofFreedom() {return fTempDegreesofFreedom;} + UInt_t GetNumIntegratedEvents() {return fMaxNumIntegratedEvents;} + Float_t GetRingInterest() {return fRingInterest;} + + ClassDef(MFindStars, 0) // Tool to find stars from DC Currents +}; + +#endif + + Index: /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MHTelAxisFromStars.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MHTelAxisFromStars.cc (revision 5120) +++ /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MHTelAxisFromStars.cc (revision 5120) @@ -0,0 +1,638 @@ +/* ======================================================================== *\ +! +! * +! * 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): Wolfgang Wittek 07/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MHTelAxisFromStars +// +// This class contains histograms for MTelAxisFromStars +// +///////////////////////////////////////////////////////////////////////////// +#include "MHTelAxisFromStars.h" + +#include + +#include +#include +#include +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" + +#include "MGeomCam.h" +#include "MBinning.h" + +#include "MStarCam.h" +#include "MStarPos.h" +#include "MSkyCamTrans.h" +#include "MSrcPosCam.h" + + +ClassImp(MHTelAxisFromStars); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Setup the histograms +// +MHTelAxisFromStars::MHTelAxisFromStars(const char *name, const char *title) + : fMm2Deg(1), fUseMmScale(kFALSE) +{ + fName = name ? name : "MHTelAxisFromStars"; + fTitle = title ? title : "Plots for MTelAxisFromStars"; + + fNStars = new TH1D("NStars", "No. of stars", 10, -0.5, 9.5); + fNStars->SetDirectory(NULL); + fNStars->SetXTitle("Numder of stars"); + fNStars->SetYTitle("Counts"); + + fNdoF = new TH1D("NdoF", "No. of degrees of freedom", 20, -0.5, 19.5); + fNdoF->SetDirectory(NULL); + fNdoF->SetXTitle("Numder of degrees of freedom"); + fNdoF->SetYTitle("Counts"); + + fLog10Chi2 = new TH1D("log10Chi2", "log10Chi2", 60, -10, 5); + fLog10Chi2->SetDirectory(NULL); + fLog10Chi2->SetXTitle("log10(Chi2)"); + fLog10Chi2->SetYTitle("Counts"); + + fChi2Prob = new TH1D("Chi2-Prob", "Chi2 probability", 40, 0.0, 1.0); + fChi2Prob->SetDirectory(NULL); + fChi2Prob->SetXTitle("Chi2 probability"); + fChi2Prob->SetYTitle("Counts"); + + + fNumIter = new TH1D("NumIter", "Number of iterations", 50, -0.5, 49.5); + fNumIter->SetDirectory(NULL); + fNumIter->SetXTitle("Number of iterations"); + fNumIter->SetYTitle("Counts"); + + + fLambda = new TH1D("Lambda", "Scale factor lambda", 80, 0.90, 1.10); + fLambda->SetDirectory(NULL); + fLambda->SetXTitle("Scale factor lambda"); + fLambda->SetYTitle("Counts"); + + + fAlfa = new TH1D("Alfa", "Rotation angle alfa", 100, -2.5, 2.5); + fAlfa->SetDirectory(NULL); + fAlfa->SetXTitle("Rotation angle alfa [\\circ]"); + fAlfa->SetYTitle("Counts"); + + + fShift = new TH2D("Shift", "Sky-Cam transformnation : (x,y) shift", + 72, -534.0707, 534.0707, 72, -534.0707, 534.0707); + fShift->SetDirectory(NULL); + fShift->SetZTitle("Counts"); + fShift->SetXTitle("x-shift [\\circ]"); + fShift->SetYTitle("y-shift [\\circ]"); + + + fEstPos1 = new TH2D("EstPos1", "Estimated position1", + 72, -534.0707, 534.0707, 72, -534.0707, 534.0707); + fEstPos1->SetDirectory(NULL); + fEstPos1->SetZTitle("Counts"); + fEstPos1->SetXTitle("Estimated position1 x [\\circ]"); + fEstPos1->SetYTitle("Estimated position1 y [\\circ]"); + + fEstPos2 = new TH2D("EstPos2", "Estimated position2", + 72, -534.0707, 534.0707, 72, -534.0707, 534.0707); + fEstPos2->SetDirectory(NULL); + fEstPos2->SetZTitle("Counts"); + fEstPos2->SetXTitle("Estimated position2 x [\\circ]"); + fEstPos2->SetYTitle("Estimated position2 y [\\circ]"); + + fEstPos3 = new TH2D("EstPos3", "Estimated position3", + 72, -534.0707, 534.0707, 72, -534.0707, 534.0707); + fEstPos3->SetDirectory(NULL); + fEstPos3->SetZTitle("Counts"); + fEstPos3->SetXTitle("Estimated position3 x [\\circ]"); + fEstPos3->SetYTitle("Estimated position3 y [\\circ]"); + + + // default input type : results from Gauss fit + fInputType = 1; +} + +// -------------------------------------------------------------------------- +// +// Set the type of the input +// +// type = 0 calculated star positions (by averaging) +// type = 1 fitted star positions (by Gauss fit) +// +void MHTelAxisFromStars::SetInputType(Int_t type) +{ + *fLog << all << "MHTelAxisFromStars::SetInputType; input type is set equal to : " + << type ; + + if (type == 0) + *fLog << " (calculated star positions)" << endl; + else + *fLog << " (fitted star positions)" << endl; + + + fInputType = type; +} + + +// -------------------------------------------------------------------------- +// +// Delete the histograms +// +MHTelAxisFromStars::~MHTelAxisFromStars() +{ + //*fLog << "MHTelAxisFromStars::Destructor" << endl; + + delete fNStars; + delete fNdoF; + delete fLog10Chi2; + delete fChi2Prob; + delete fNumIter; + delete fLambda; + delete fAlfa; + + delete fShift; + delete fEstPos1; + delete fEstPos2; + delete fEstPos3; +} + +// -------------------------------------------------------------------------- +// +// Setup the Binning for the histograms. +// +// Find the pointers to the containers. +// +// Use this function if you want to set the conversion factor which +// is used to convert the mm-scale in the camera plane into the deg-scale. +// The conversion factor is part of the camera geometry. Please create a +// corresponding MGeomCam container. +// +Bool_t MHTelAxisFromStars::SetupFill(const MParList *plist) +{ + + fStarCam = (MStarCam*)plist->FindObject("MStarCam", "MStarCam"); + if (!fStarCam) + { + *fLog << err << "MHTelAxisFromStars::SetupFill; container 'MStarCam' not found... aborting." << endl; + return kFALSE; + } + + + fSourceCam = (MStarCam*)plist->FindObject("MSourceCam", "MStarCam"); + if (!fSourceCam) + { + *fLog << err << "MHTelAxisFromStars::SetupFill; container 'MSourceCam' not found... " << endl; + } + + + + fSrcPos = (MSrcPosCam*)plist->FindObject(AddSerialNumber("MSrcPosCam")); + if (!fSrcPos) + { + *fLog << err << "MHTelAxisFromStars::SetupFill; MSrcPosCam not found... aborting" << endl; + return kFALSE; + } + + fSkyCamTrans = (MSkyCamTrans*)plist->FindObject(AddSerialNumber("MSkyCamTrans")); + if (!fSkyCamTrans) + { + *fLog << err << "MHTelAxisFromStars::SetupFill; MSkyCamTrans not found... aborting" << endl; + return kFALSE; + } + + + //------------------------------------------------------------------ + const MGeomCam *geom = (MGeomCam*)plist->FindObject("MGeomCam"); + if (!geom) + { + *fLog << warn << GetDescriptor() + << ": No Camera Geometry available. Using mm-scale for histograms." + << endl; + SetMmScale(kTRUE); + } + else + { + fMm2Deg = geom->GetConvMm2Deg(); + SetMmScale(kFALSE); + } + + ApplyBinning(*plist, "NStars", fNStars); + ApplyBinning(*plist, "NdoF", fNdoF); + ApplyBinning(*plist, "Log10Chi2",fLog10Chi2); + ApplyBinning(*plist, "Chi2Prob", fChi2Prob); + ApplyBinning(*plist, "NumIter", fNumIter); + ApplyBinning(*plist, "Lambda", fLambda); + ApplyBinning(*plist, "Alfa", fAlfa); + + const MBinning *bins = (MBinning*)plist->FindObject("BinningCamera"); + if (!bins) + { + float r = geom ? geom->GetMaxRadius() : 600; + r *= 0.9; + if (!fUseMmScale) + r *= fMm2Deg; + + MBinning b; + b.SetEdges(61, -r, r); + SetBinning(fShift, &b, &b); + SetBinning(fEstPos1, &b, &b); + SetBinning(fEstPos2, &b, &b); + SetBinning(fEstPos3, &b, &b); + } + else + { + SetBinning(fShift, bins, bins); + SetBinning(fEstPos1, bins, bins); + SetBinning(fEstPos2, bins, bins); + SetBinning(fEstPos3, bins, bins); + } + + //------------------------------------------- + // copy names from MStarCam to the histograms + MStarPos* starpos; + Int_t istar = 0; + TIter Next(fSourceCam->GetList()); + + while ((starpos=(MStarPos*)Next())) { + fStarnames[istar] = starpos->GetName(); + //*fLog << "istar, star name = " << istar << ", " + // << fStarnames[istar] << endl; + istar++; + if (istar >= fNstarnames) break; + } + + if (fSourceCam) + { + MStarPos *starSource = 0; + TIter nextSource(fSourceCam->GetList()); + while ( (starSource = (MStarPos*)nextSource()) ) + { + if ( fNstarnames > 0 && starSource->GetName() == fStarnames[0] ) + fEstPos1->SetName(starSource->GetName()); + + else if( fNstarnames > 1 && starSource->GetName() == fStarnames[1] ) + fEstPos2->SetName(starSource->GetName()); + + else if( fNstarnames > 2 && starSource->GetName() == fStarnames[2] ) + fEstPos3->SetName(starSource->GetName()); + } + } + + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Use this function to setup your own conversion factor between degrees +// and millimeters. The conversion factor should be the one calculated in +// MGeomCam. Use this function with Caution: You could create wrong values +// by setting up your own scale factor. +// +void MHTelAxisFromStars::SetMm2Deg(Float_t mmdeg) +{ + if (mmdeg<0) + { + *fLog << warn << dbginf + << "Warning - Conversion factor < 0 - nonsense. Ignored." + << endl; + return; + } + + if (fMm2Deg>=0) + *fLog << warn << dbginf + << "Warning - Conversion factor already set. Overwriting" + << endl; + + fMm2Deg = mmdeg; +} + +// -------------------------------------------------------------------------- +// +// With this function you can convert the histogram ('on the fly') between +// degrees and millimeters. +// +void MHTelAxisFromStars::SetMmScale(Bool_t mmscale) +{ + if (fUseMmScale == mmscale) + return; + + if (fMm2Deg<0) + { + *fLog << warn << dbginf + << "Warning - Sorry, no conversion factor for conversion available." + << endl; + return; + } + + const Double_t scale = mmscale ? 1./fMm2Deg : fMm2Deg; + MH::ScaleAxis(fShift, scale, scale); + MH::ScaleAxis(fEstPos1, scale, scale); + MH::ScaleAxis(fEstPos2, scale, scale); + MH::ScaleAxis(fEstPos3, scale, scale); + + if (mmscale) + { + fShift->SetXTitle("x-shift [mm]"); + fShift->SetYTitle("y-shift [mm]"); + + fEstPos1->SetXTitle("Estimated position1 x [mm]"); + fEstPos1->SetYTitle("Estimated position1 y [mm]"); + + fEstPos2->SetXTitle("Estimated position2 x [mm]"); + fEstPos2->SetYTitle("Estimated position2 y [mm]"); + + fEstPos3->SetXTitle("Estimated position3 x [mm]"); + fEstPos3->SetYTitle("Estimated position3 y [mm]"); + } + else + { + fShift->SetXTitle("x-shift [\\circ]"); + fShift->SetYTitle("y-shift [\\circ]"); + + fEstPos1->SetXTitle("Estimated position1 x [\\circ]"); + fEstPos1->SetYTitle("Estimated position1 y [\\circ]"); + + fEstPos2->SetXTitle("Estimated position2 x [\\circ]"); + fEstPos2->SetYTitle("Estimated position2 y [\\circ]"); + + fEstPos3->SetXTitle("Estimated position3 x [\\circ]"); + fEstPos3->SetYTitle("Estimated position3 y [\\circ]"); + } + + fUseMmScale = mmscale; +} + +// -------------------------------------------------------------------------- +// +// Fill the histograms +// +// +Bool_t MHTelAxisFromStars::Fill(const MParContainer *par, const Stat_t w) +{ + Int_t Ndof = fSkyCamTrans->GetNdof(); + if (Ndof < 0) + return kTRUE; + + const Double_t scale = fUseMmScale ? 1 : fMm2Deg; + + + fNStars ->Fill(fSkyCamTrans->GetNStars(), w); + fNdoF ->Fill(fSkyCamTrans->GetNdof(), w); + if (fSkyCamTrans->GetChiSquare() > 0.0) + fLog10Chi2->Fill( log10(fSkyCamTrans->GetChiSquare() ), w); + fChi2Prob ->Fill(fSkyCamTrans->GetChiSquareProb(), w); + fNumIter ->Fill(fSkyCamTrans->GetNumIter(), w); + fLambda ->Fill(fSkyCamTrans->GetLambda(), w); + fAlfa ->Fill(fSkyCamTrans->GetAlfa(), w); + + + Double_t lowx; + Double_t higx; + Double_t lowy; + Double_t higy; + Double_t x; + Double_t y; + Int_t nbix; + Int_t nbiy; + + nbix = fShift->GetXaxis()->GetNbins(); + lowx = fShift->GetXaxis()->GetBinLowEdge(1); + higx = fShift->GetXaxis()->GetBinLowEdge(nbix+1); + + nbiy = fShift->GetYaxis()->GetNbins(); + lowy = fShift->GetYaxis()->GetBinLowEdge(1); + higy = fShift->GetYaxis()->GetBinLowEdge(nbiy+1); + + x = scale * (fSkyCamTrans->GetShiftD())[0]; + y = scale * (fSkyCamTrans->GetShiftD())[1]; + if (x>lowx && xlowy && yFill(x, y, w); + } + + if (fSourceCam) + { + MStarPos *starSource = 0; + TIter nextSource(fSourceCam->GetList()); + + while( (starSource = (MStarPos*)nextSource()) ) + { + if (fInputType == 1) + { + x = scale * starSource->GetMeanXFit(); + y = scale * starSource->GetMeanYFit(); + } + else + { + x = scale * starSource->GetMeanXCalc(); + y = scale * starSource->GetMeanYCalc(); + } + + if (x>lowx && xlowy && y 0 && starSource->GetName() == fStarnames[0] ) + fEstPos1->Fill(x, y, w); + + else if( fNstarnames > 1 && starSource->GetName() == fStarnames[1] ) + fEstPos2->Fill(x, y, w); + + else if( fNstarnames > 2 && starSource->GetName() == fStarnames[2] ) + fEstPos3->Fill(x, y, w); + } + } + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Finalize : +// +// it is called in the postprocessing +// reset pointers in order to allow cloning of the object +// +Bool_t MHTelAxisFromStars::Finalize() +{ + //*fLog << "MHTelAxisFromStars::Finalize; fSourceCam = " + // << fSourceCam << endl; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Creates a new canvas and draws the histograms. +// Be careful: The histograms belongs to this object and won't get deleted +// together with the canvas. +// + +void MHTelAxisFromStars::Draw(Option_t *opt) +{ + TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); + pad->SetBorderMode(0); + + AppendPad(""); + + //TCanvas *pad = new TCanvas("TelAxisFromStars", "TelAxis plots", 1200, 900); + //gROOT->SetSelectedPad(NULL); + + pad->Divide(4,3); + + pad->cd(1); + gPad->SetBorderMode(0); + fNStars->Draw(opt); + + pad->cd(2); + gPad->SetBorderMode(0); + fNdoF->Draw(opt); + + pad->cd(3); + gPad->SetBorderMode(0); + fLog10Chi2->Draw(opt); + + pad->cd(4); + gPad->SetBorderMode(0); + fChi2Prob->Draw(opt); + + pad->cd(5); + gPad->SetBorderMode(0); + fNumIter->Draw(opt); + + pad->cd(6); + gPad->SetBorderMode(0); + fLambda->Draw(opt); + + pad->cd(7); + gPad->SetBorderMode(0); + fAlfa->Draw(opt); + + pad->cd(8); + gPad->SetBorderMode(0); + //SetColors(); + //fShift->Draw("colz"); + fShift->Draw(""); + + + //----------------------------------------------- + // plot the expected positions of some sources + //*fLog << "fSourcsCam = " << fSourceCam << endl; + + if (fSourceCam) + { + *fLog << "MHTelAxisFromSrars::Draw; plotting" << endl; + + pad->cd(9); + gPad->SetBorderMode(0); + //SetColors(); + fEstPos1->Draw(""); + + pad->cd(10); + gPad->SetBorderMode(0); + //SetColors(); + fEstPos2->Draw(""); + + pad->cd(11); + gPad->SetBorderMode(0); + //SetColors(); + fEstPos3->Draw(""); + } + + pad->Modified(); + pad->Update(); +} + +//--------------------------------------------------------------------- +// + +TH1 *MHTelAxisFromStars::GetHistByName(const TString name) +{ + if (name.Contains("NStars", TString::kIgnoreCase)) + return fNStars; + if (name.Contains("NdoF", TString::kIgnoreCase)) + return fNdoF; + if (name.Contains("Log10Chi2", TString::kIgnoreCase)) + return fLog10Chi2; + if (name.Contains("Chi2Prob", TString::kIgnoreCase)) + return fChi2Prob; + if (name.Contains("NumIter", TString::kIgnoreCase)) + return fNumIter; + if (name.Contains("Lambda", TString::kIgnoreCase)) + return fLambda; + if (name.Contains("Alfa", TString::kIgnoreCase)) + return fAlfa; + + if (name.Contains("Shift", TString::kIgnoreCase)) + return fShift; + if (name.Contains("EstPos1", TString::kIgnoreCase)) + return fEstPos1; + if (name.Contains("EstPos2", TString::kIgnoreCase)) + return fEstPos2; + if (name.Contains("EstPos3", TString::kIgnoreCase)) + return fEstPos3; + + return NULL; +} + +// -------------------------------------------------------------------------- +// +// Setup an inversed deep blue sea palette. +// + +void MHTelAxisFromStars::SetColors() const +{ + gStyle->SetPalette(51, NULL); + Int_t c[50]; + for (int i=0; i<50; i++) + c[49-i] = gStyle->GetColorPalette(i); + gStyle->SetPalette(50, c); +} + + +//--------------------------------------------------------------------- +// + +void MHTelAxisFromStars::Paint(Option_t *opt) +{ + SetColors(); + MH::Paint(); +} + +//========================================================================== + + + + + + + + Index: /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MHTelAxisFromStars.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MHTelAxisFromStars.h (revision 5120) +++ /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MHTelAxisFromStars.h (revision 5120) @@ -0,0 +1,86 @@ +#ifndef MARS_MHTelAxisFromStars +#define MARS_MHTelAxisFromStars + +#ifndef MARS_MH +#include "MH.h" +#endif + +class TH1D; +class TH2D; +class MTelAxisFromStars; +class MSrcPosCam; +class MSkyCamTrans; +class MStarCam; + +class MHTelAxisFromStars : public MH +{ +private: + + TH1D *fNStars; //-> Number of stars + TH1D *fNdoF; //-> Number of degrees of freedom + TH1D *fLog10Chi2; //-> log10(Chi2) + TH1D *fChi2Prob; //-> Chi2 probability + TH1D *fNumIter; //-> Number of iterations + TH1D *fLambda; //-> Scale factor lambda + TH1D *fAlfa; //-> Rotation angle alfa + + TH2D *fShift; //-> Shift between Sky and Camera system + TH2D *fEstPos1; //-> Estimated position 1 + TH2D *fEstPos2; //-> Estimated position 2 + TH2D *fEstPos3; //-> Estimated position 3 + + MStarCam *fStarCam; //! + MStarCam *fSourceCam; + MSrcPosCam *fSrcPos; //! + MSkyCamTrans *fSkyCamTrans; //! + + Float_t fMm2Deg; //! + Bool_t fUseMmScale; //! + Int_t fInputType; //! + + const static Int_t fNstarnames = 3; //! + TString fStarnames[fNstarnames]; //! + + void SetColors() const; + void Paint(Option_t *opt=""); + +public: + MHTelAxisFromStars(const char *name=NULL, const char *title=NULL); + ~MHTelAxisFromStars(); + + void SetInputType(Int_t type=1); + void SetMmScale(Bool_t mmscale=kTRUE); + virtual void SetMm2Deg(Float_t mmdeg); + + Bool_t SetupFill(const MParList *pList); + Bool_t Fill(const MParContainer *par, const Stat_t w=1); + Bool_t Finalize(); + + TH1 *GetHistByName(const TString name); + + TH1D *GetHistNStars() { return fNStars; } + TH1D *GetHistNdoF() { return fNdoF; } + TH1D *GetHistLog10Chi2() { return fLog10Chi2; } + TH1D *GetHistChi2Prob() { return fChi2Prob; } + + TH1D *GetHistNumIter() { return fNumIter; } + TH1D *GetHistLambda() { return fLambda; } + + TH1D *GetHistAlfa() { return fAlfa; } + + TH2D *GetHistShift() { return fShift; } + TH2D *GetHistEstPos1() { return fEstPos1; } + TH2D *GetHistEstPos2() { return fEstPos2; } + TH2D *GetHistEstPos3() { return fEstPos3; } + + void Draw(Option_t *opt=NULL); + + ClassDef(MHTelAxisFromStars, 1) // Container which holds histograms for MTelAxisFromStars +}; + +#endif + + + + + Index: /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MSkyCamTrans.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MSkyCamTrans.cc (revision 5120) +++ /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MSkyCamTrans.cc (revision 5120) @@ -0,0 +1,118 @@ +/* ======================================================================== *\ +! +! * +! * 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): Wolfgang Wittek , 7/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +/////////////////////////////////////////////////////////////////////////// +// +// MSkyCamTrans +// +// container holding the parameters of the transformation from +// sky directions 'a' (projected onto the camera) +// to positions 'b' in the camera +// +// ( cos(fAlfa) -sin(fAlfa) ) +// b = fLambda * fA * a + fD fA = ( ) +// ^ ^ ^ ( sin(fAlfa) cos(fAlfa) ) +// | | | +// scale rotation shift +// factor matrix +// +// fNumIter number of iterations +// fNdof number of degrees of freedom +// fChiSquare chi-square value +// fChiSquareProb chi-square probability +// +// The units are assumed to be +// [degrees] for fAlfa +// [mm] for a, b, fD +// [1] for fLambda +// +// +// The container is filled by the task 'MTelAxisFromStars' +// +/////////////////////////////////////////////////////////////////////////// + + +#include "MSkyCamTrans.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MSkyCamTrans); + +using namespace std; + +// --------------------------------------------------------------------- +// +// +// +MSkyCamTrans::MSkyCamTrans(const char *name, const char *title) +{ + fName = name ? name : "MSkyCamTrans"; + fTitle = title ? title : "Sky-Cam transformation parameters"; +} + +// --------------------------------------------------------------------- +// +// +// +void MSkyCamTrans::SetParameters(Double_t &lambda, Double_t &alfa, + Double_t a[2][2], Double_t d[2], Double_t errd[2][2], + Int_t &nstars, Int_t &numiter, + Int_t &ndof, Double_t &chisquare, Double_t &chisquareprob) +{ + fLambda = lambda; + fAlfa = alfa; + + fA[0][0] = a[0][0]; + fA[0][1] = a[0][1]; + fA[1][0] = a[1][0]; + fA[1][1] = a[1][1]; + + fD[0] = d[0]; + fD[1] = d[1]; + + fErrD[0][0] = errd[0][0]; + fErrD[0][1] = errd[0][1]; + fErrD[1][0] = errd[1][0]; + fErrD[1][1] = errd[1][1]; + + fNStars = nstars; + fNumIter = numiter; + fNdof = ndof; + fChiSquare = chisquare; + fChiSquareProb = chisquareprob; +} +// --------------------------------------------------------------------- + + + + + + + + + + + + Index: /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MSkyCamTrans.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MSkyCamTrans.h (revision 5120) +++ /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MSkyCamTrans.h (revision 5120) @@ -0,0 +1,55 @@ +#ifndef MARS_MSkyCamTrans +#define MARS_MSkyCamTrans + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +class MSkyCamTrans : public MParContainer +{ +private: + + //Parameters of transformation from sky to camera + + Double_t fLambda; // scale factor lambda + Double_t fAlfa; // rotation angle alfa [degrees] + Double_t fA[2][2]; // rotation matrix A + Double_t fD[2]; // shift D [mm] + Double_t fErrD[2][2]; // error matrix of shift D [mm*mm] + + Int_t fNStars; + Int_t fNumIter; + Int_t fNdof; + Double_t fChiSquare; + Double_t fChiSquareProb; + +public: + + MSkyCamTrans(const char *name=NULL, const char *title=NULL); + + void SetParameters(Double_t &, Double_t &, + Double_t[2][2], Double_t[2], Double_t[2][2], + Int_t &, Int_t &, Int_t &, Double_t &, Double_t &); + + Int_t GetNStars() { return fNStars; } + Int_t GetNumIter() { return fNumIter; } + Int_t GetNdof() { return fNdof; } + Double_t GetChiSquare() { return fChiSquare; } + Double_t GetChiSquareProb() { return fChiSquareProb; } + Double_t GetLambda() { return fLambda; } + Double_t GetAlfa() { return fAlfa; } + + Double_t *GetRotationMatrix() { return &fA[0][0]; } + Double_t *GetShiftD() { return &fD[0]; } + Double_t *GetErrMatrixShiftD() { return &fErrD[0][0]; } + + ClassDef(MSkyCamTrans, 1) // Container holding the sky-camera transformation parameters +}; + +#endif + + + + + + Index: /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MSourceDirections.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MSourceDirections.cc (revision 5120) +++ /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MSourceDirections.cc (revision 5120) @@ -0,0 +1,190 @@ +/* ======================================================================== *\ +! +! * +! * 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 expressed +! * or implied warranty. +! * +! +! Author(s): Robert Wagner, 8/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MSourceDirections +// +///////////////////////////////////////////////////////////////////////////// +#include "MSourceDirections.h" + +#include +#include +#include +#include +#include +#include +#include +#include + + +#include "MObservatory.h" +#include "MAstroCamera.h" +#include "MMcConfigRunHeader.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MHCamera.h" + +#include "MGeomCam.h" +#include "MGeomPix.h" +#include "MCameraDC.h" +#include "MTime.h" +#include "MReportDrive.h" +#include "MStarCam.h" +#include "MStarPos.h" + +#include "MParList.h" +#include "MTaskList.h" + +ClassImp(MSourceDirections); +using namespace std; + +MSourceDirections::MSourceDirections(const char *name, const char *title): + fGeomCam(NULL), fTimeCurr(NULL), fDrive(NULL), fStars(NULL) +{ + fName = name ? name : "MSourceDirections"; + fTitle = title ? title : "Convert Ra-Dec source positions into camera coordinates"; + fGeometryFile=""; +} + +Int_t MSourceDirections::PreProcess(MParList *pList) +{ + fGeomCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam")); + if (!fGeomCam) { + *fLog << err << AddSerialNumber("MGeomCam") << " not found ... aborting" << endl; + return kFALSE; + } + + fTimeCurr = (MTime*)pList->FindObject(AddSerialNumber("MTimeCurrents")); + if (!fTimeCurr) { + *fLog << err << AddSerialNumber("MTimeCurrents") << " not found ... aborting" << endl; + return kFALSE; + } + + fDrive = (MReportDrive*)pList->FindObject(AddSerialNumber("MReportDrive")); + if (!fDrive) { + *fLog << warn << AddSerialNumber("MReportDrive") << " not found ... aborting" << endl; + return kFALSE; + } + + + fStars = (MStarCam*)pList->FindCreateObj(AddSerialNumber("MStarCam"),"MSourceCam"); + if (!fStars) { + *fLog << err << AddSerialNumber("MStarCam") << " with name '" + << "MSourceCam" << " cannot be created ... aborting" << endl; + return kFALSE; + } + + + + MObservatory magic1; + + MMcConfigRunHeader *config=0; + MGeomCam *geom=0; + + TFile file(fGeometryFile); + //TTree *tree = (TTree*)file.Get("RunHeaders"); + TTree *tree = dynamic_cast(file.Get("RunHeaders")); + tree->SetBranchAddress("MMcConfigRunHeader.", &config); + if (tree->GetBranch("MGeomCam.")) tree->SetBranchAddress("MGeomCam.", &geom); + tree->GetEntry(0); + + fAstro.SetMirrors(*config->GetMirrors()); + fAstro.SetGeom(*geom); + fAstro.SetObservatory(magic1); + + + // make the starlist available already in the preprocessing + fStars->GetList()->Delete(); + fAstro.SetTime(*fTimeCurr); + fAstro.SetGuiActive(); + fAstro.FillStarList(fStars->GetList()); + //*fLog << inf << "in preprocessing " << GetName() << " found " + // << fStars->GetList()->GetSize() + // << " directions inside the chosen FOV." << endl; + + + return kTRUE; +} + +Int_t MSourceDirections::AddDirection(Float_t ra, Float_t dec, Float_t mag, TString name) +{ + *fLog << "MSourceDirections::AddDirection; add the direction : ra, dec, mag, name = " + << ra << ", " << dec << ", " << mag << ", " << name << endl; + + Int_t rc = fAstro.AddObject(ra,dec,1,name); +// if (rc) { +// MStarPos *starpos = new MStarPos; +// starpos->SetName(name); +// fStars->GetList()->Add(starpos); +// } + return rc; +} + +Int_t MSourceDirections::Process() +{ + //First delete the previous directions in the list + fStars->GetList()->Delete(); + + fAstro.SetTime(*fTimeCurr); + fAstro.SetGuiActive(); + fAstro.FillStarList(fStars->GetList()); + + MStarPos* starpos; + TIter Next(fStars->GetList()); + while ((starpos=(MStarPos*)Next())) { + //starpos->SetCalcValues(40,40,starpos->GetXExp(),starpos->GetYExp(), + // 0.,0.,0., 0.,0.,0.); + //starpos->SetFitValues (40,40,starpos->GetXExp(),starpos->GetYExp(), + // 0.,0.,0., 0.,0.,0., 0., -1); + } + + if (fStars->GetList()->GetSize() == 0) { + *fLog << err << GetName() << "No directions inside the chosen FOV." << endl; + } + else { + *fLog << inf << GetName() << " found " << fStars->GetList()->GetSize() + << " directions inside the chosen FOV." << endl; + } + return kTRUE; +} + +Int_t MSourceDirections::PostProcess() +{ + return kTRUE; +} + + + + + + + + + + + + + Index: /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MSourceDirections.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MSourceDirections.h (revision 5120) +++ /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MSourceDirections.h (revision 5120) @@ -0,0 +1,49 @@ +#ifndef MARS_MSourceDirections +#define MARS_MSourceDirections + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#ifndef MARS_MAstroCamera +#include "MAstroCamera.h" +#endif + +class MTime; +class MReportDrive; +class MStarCam; +class MStarPos; + +class MSourceDirections : public MTask +{ + +private: + + MGeomCam *fGeomCam; + MTime *fTimeCurr; + MReportDrive *fDrive; + MStarCam *fStars; + + MAstroCamera fAstro; + + TString fGeometryFile; + + public: + + MSourceDirections(const char *name=NULL, const char *title=NULL); + + Int_t AddDirection(Float_t ra, Float_t dec, Float_t mag, TString name=""); + Int_t PreProcess(MParList *pList); + Int_t Process(); + Int_t PostProcess(); + + // setters + void SetGeometryFile(TString f) {fGeometryFile=f;} + void SetRaDec(Double_t ra, Double_t dec) {fAstro.SetRaDec(ra,dec);} + void SetRaDec(TVector3 &v) { fAstro.SetRaDec(v); } + void SetRadiusFOV(Double_t deg) { fAstro.SetRadiusFOV(deg); } + + ClassDef(MSourceDirections, 0) // Tool to translate RaDec Source Directions into camera coordinates +}; + +#endif Index: /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MStarLocalCam.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MStarLocalCam.cc (revision 5120) +++ /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MStarLocalCam.cc (revision 5120) @@ -0,0 +1,111 @@ +/* ======================================================================== *\ +! +! * +! * 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): Javier Lopez 04/2004 +! Author(s): Jordi Albert 04/2004 +! +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MStarLocalCam +// // +// // +///////////////////////////////////////////////////////////////////////////// + +#include "MStarLocalCam.h" + +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MStarLocalPos.h" + +using namespace std; + +// -------------------------------------------------------------------------- +// Default constructor. +// +// +MStarLocalCam::MStarLocalCam(const char *name, const char *title) +{ + fName = name ? name : "MStarLocalCam"; + fTitle = title ? title : ""; + + fStars = new TList; + + fInnerPedestalDC = 0.; + fOuterPedestalDC = 0.; + fInnerPedestalRMSDC = 0.; + fOuterPedestalRMSDC = 0.; + +} + +MStarLocalCam::~MStarLocalCam() +{ + fStars->SetOwner(); + fStars->Delete(); +} + +// -------------------------------------------------------------------------- +// +// Get i-th +// +MStarLocalPos &MStarLocalCam::operator[] (Int_t i) +{ + MStarLocalPos& star = *static_cast(fStars->At(i)); + return star; +} + +// -------------------------------------------------------------------------- +// +// Get i-th +// +const MStarLocalPos &MStarLocalCam::operator[] (Int_t i) const +{ + return *static_cast(fStars->At(i)); +} + +void MStarLocalCam::Paint(Option_t *o) +{ + TIter Next(fStars); + MStarLocalPos* star; + while ((star=(MStarLocalPos*)Next())) + star->Paint(o); +} + +void MStarLocalCam::Print(Option_t *o) const +{ + *fLog << inf << "DCs baseline:" << endl; + *fLog << inf << " Inner Pixels Mean pedestal \t" << setw(4) << fInnerPedestalDC << " uA and RMS " << setw(4) << fInnerPedestalRMSDC << " uA for DCs" << endl; + *fLog << inf << " Outer Pixels Mean pedestal \t" << setw(4) << fOuterPedestalDC << " uA and RMS " << setw(4) << fOuterPedestalRMSDC << " uA for DCs" << endl; + + TIter Next(fStars); + MStarLocalPos* star; + UInt_t starnum = 0; + while ((star=(MStarLocalPos*)Next())) + { + *fLog << inf << "Star[" << starnum << "] info:" << endl; + star->Print(o); + starnum++; + } +} Index: /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MStarLocalCam.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MStarLocalCam.h (revision 5120) +++ /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MStarLocalCam.h (revision 5120) @@ -0,0 +1,64 @@ +#ifndef MARS_MStarLocalCam +#define MARS_MStarLocalCam + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +#ifndef MARS_MHCamera +#include "MHCamera.h" +#endif + +class TList; + +class MGeomCam; +class MStarLocalPos; + +class MStarLocalCam : public MParContainer +{ +private: + + TList *fStars; //-> FIXME: Change TClonesArray away from a pointer? + + // BaseLine DCs info + Float_t fInnerPedestalDC; //[ua] + Float_t fOuterPedestalDC; //[ua] + + Float_t fInnerPedestalRMSDC; //[ua] + Float_t fOuterPedestalRMSDC; //[ua] + + MHCamera fDisplay; + + public: + + MStarLocalCam(const char *name=NULL, const char *title=NULL); + ~MStarLocalCam(); + + MStarLocalPos &operator[] (Int_t i); + const MStarLocalPos &operator[] (Int_t i) const; + + TList *GetList() const { return fStars; } + UInt_t GetNumStars() const { return fStars->GetSize(); } + + //Getters + + Float_t GetInnerPedestalDC() {return fInnerPedestalDC;} + Float_t GetOuterPedestalDC() {return fOuterPedestalDC;} + Float_t GetInnerPedestalRMSDC() { return fInnerPedestalRMSDC;} + Float_t GetOuterPedestalRMSDC() { return fOuterPedestalRMSDC;} + + MHCamera& GetDisplay() { return fDisplay; } + + //Setters + void SetInnerPedestalDC(Float_t ped) {fInnerPedestalDC = ped;} + void SetOuterPedestalDC(Float_t ped) {fOuterPedestalDC = ped;} + void SetInnerPedestalRMSDC(Float_t rms){fInnerPedestalRMSDC = rms;} + void SetOuterPedestalRMSDC(Float_t rms){fOuterPedestalRMSDC = rms;} + + void Paint(Option_t *o=NULL); + void Print(Option_t *o=NULL) const; + + ClassDef(MStarLocalCam, 1) // Storage Container for star positions in the camera +}; + +#endif Index: /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MStarLocalPos.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MStarLocalPos.cc (revision 5120) +++ /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MStarLocalPos.cc (revision 5120) @@ -0,0 +1,309 @@ +/* ======================================================================== *\ +! +! * +! * 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 expressed +! * or implied warranty. +! * +! +! +! Author(s): Javier López , 4/2004 +! Robert Wagner, 7/2004 +! Wolfgang Wittek, 8/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ +#include "MStarLocalPos.h" + +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MStarLocalPos); + +using namespace std; + +MStarLocalPos::MStarLocalPos(const char *name, const char *title) +{ + + fName = name ? name : "MStarLocalPos"; + fTitle = title ? title : ""; + + Reset(); +} + +void MStarLocalPos::Reset() +{ + + //Expected position on camera + fMagExp = 0.; + fXExp = 0.; + fYExp = 0.; + + //Info from calculation + fMagCalc = 0.; + fMaxCalc = 0.; + fMeanXCalc = 0.; + fMeanYCalc = 0.; + fSigmaMinorAxisCalc = 0.; + fSigmaMajorAxisCalc = 0.; + + //Info from uncorrelated Gauss fit + fMagFit = 0.; + fMaxFit = 0.; + fMeanXFit = 0.; + fMeanYFit = 0.; + fSigmaMinorAxisFit = 0.; + fSigmaMajorAxisFit = 0.; + + fChiSquare = 0.; + fNdof = 1; + + //Info from correlated Gauss fit + fMagCGFit = 0.; + fMaxCGFit = 0.; + fMeanXCGFit = 0.; + fMeanYCGFit = 0.; + fSigmaXCGFit = 0.; + fSigmaYCGFit = 0.; + fCorrXYCGFit = 0.; + fXXErrCGFit = 0.; + fXYErrCGFit = 0.; + fYYErrCGFit = 0.; + + fChiSquareCGFit = 0.; + fNdofCGFit = 1; + +} + +void MStarLocalPos::SetExpValues(Float_t mag, Float_t x, Float_t y) +{ + fMagExp = mag; + fXExp = x; + fYExp = y; +} + +void MStarLocalPos::SetIdealValues(Float_t mag, Float_t x, Float_t y) +{ + fMagIdeal = mag; + fXIdeal = x; + fYIdeal = y; +} + +void MStarLocalPos::SetCalcValues(Float_t mag, Float_t max, + Float_t x, Float_t y, Float_t sigmaMinorAxis, Float_t sigmaMajorAxis) +{ + fMagCalc = mag; + fMaxCalc = max; + fMeanXCalc = x; + fMeanYCalc = y; + fSigmaMinorAxisCalc = sigmaMinorAxis; + fSigmaMajorAxisCalc = sigmaMajorAxis; +} + +void MStarLocalPos::SetFitValues(Float_t mag, Float_t max, + Float_t x, Float_t y, Float_t sigmaMinorAxis, Float_t sigmaMajorAxis, + Float_t chiSquare, Int_t ndof) +{ + fMagFit = mag; + fMaxFit = max; + fMeanXFit = x; + fMeanYFit = y; + fSigmaMinorAxisFit = sigmaMinorAxis; + fSigmaMajorAxisFit = sigmaMajorAxis; + fChiSquare = chiSquare; + fNdof = ndof; +} + +void MStarLocalPos::SetFitValues(Float_t mag, Float_t max, + Float_t x, Float_t y, Float_t sigmaMinorAxis, Float_t sigmaMajorAxis, + Float_t chiSquare, Int_t ndof, + Float_t xx, Float_t xy, Float_t yy) +{ + SetFitValues(mag, max, x, y, sigmaMinorAxis, sigmaMajorAxis, chiSquare, ndof); + fXXErr = xx; + fYYErr = yy; + fXYErr = xy; +} + +void MStarLocalPos::SetCGFitValues( + Float_t mag, Float_t max, Float_t x, Float_t y, + Float_t sigmaX, Float_t sigmaY, Float_t correlation, + Float_t xx, Float_t xy, Float_t yy, + Float_t chiSquare, Int_t ndof) +{ + fMagCGFit = mag; + fMaxCGFit = max; + fMeanXCGFit = x; + fMeanYCGFit = y; + fSigmaXCGFit = sigmaX; + fSigmaYCGFit = sigmaY; + fCorrXYCGFit = correlation; + fXXErrCGFit = xx; + fXYErrCGFit = xy; + fYYErrCGFit = yy; + + fChiSquareCGFit = chiSquare; + fNdofCGFit = ndof; +} + + +// -------------------------------------------------------------------------- +// +// Paint the ellipse corresponding to the parameters +// +void MStarLocalPos::Paint(Option_t *opt) +{ + //Print a cross in the expected position + TMarker mexp(fXExp, fYExp, 29); + mexp.SetMarkerSize(3); + mexp.SetMarkerColor(94); + mexp.Paint(); + + if (fSigmaMinorAxisCalc>0. && fSigmaMajorAxisCalc>0.) + { + TEllipse ecalc(fMeanXCalc, fMeanYCalc, fSigmaMinorAxisCalc, fSigmaMajorAxisCalc, 0, 360, 0); + ecalc.SetLineWidth(3); + ecalc.SetLineColor(kBlue); + ecalc.Paint(); + } + + if (fSigmaMinorAxisFit>0. && fSigmaMajorAxisFit>0.) + { + TEllipse efit(fMeanXFit, fMeanYFit, fSigmaMinorAxisFit, fSigmaMajorAxisFit, 0, 360, 0); + efit.SetLineWidth(3); + efit.SetLineColor(kBlack); + efit.Paint(); + } + + if (fSigmaXCGFit>0. && fSigmaYCGFit>0.) + { + //Print a cross in the fitted position + //TMarker mCGFit(fMeanXCGFit, fMeanYCGFit, 3); + //mCGFit.SetMarkerSize(3); + //mCGFit.SetMarkerColor(1); + //mCGFit.Paint(); + + Double_t cxx = fSigmaXCGFit*fSigmaXCGFit; + Double_t cyy = fSigmaYCGFit*fSigmaYCGFit; + Double_t d = cyy - cxx; + Double_t cxy = fCorrXYCGFit * fSigmaXCGFit * fSigmaYCGFit; + Double_t tandel; + if (cxy != 0.0) + tandel = ( d + sqrt(d*d + 4.0*cxy*cxy) ) / (2.0*cxy); + else + tandel = 0.0; + + Double_t sindel = tandel / sqrt(1.0 + tandel*tandel); + Double_t delta = TMath::ASin(sindel); + + Double_t major = (cxx + 2.0*tandel*cxy + tandel*tandel*cyy) + / (1.0 + tandel*tandel); + + Double_t minor = (tandel*tandel*cxx - 2.0*tandel*cxy + cyy) + / (1.0 + tandel*tandel); + + TEllipse efit(fMeanXCGFit, fMeanYCGFit, sqrt(major), sqrt(minor), + 0, 360, delta*kRad2Deg); + efit.SetLineWidth(3); + efit.SetLineColor(kMagenta); + efit.Paint(); + } +} + +void MStarLocalPos::Print(Option_t *opt) const +{ + TString o = opt; + + if (o.Contains("name", TString::kIgnoreCase) || opt == NULL) + { + *fLog << inf << "Star Name: \"" << this->GetName() << "\"" << endl; + } + + if (o.Contains("mag", TString::kIgnoreCase) || opt == NULL) + { + + *fLog << inf << "Star magnitude:" << endl; + *fLog << inf << " Expected \t" << setw(4) << fMagExp << endl; + *fLog << inf << " Calcultated \t " << setw(4) << fMagCalc << endl; + *fLog << inf << " Fitted \t " << setw(4) << fMagFit << endl; + *fLog << inf << " CGFitted \t " << setw(4) << fMagCGFit << endl; + } + + if (o.Contains("max", TString::kIgnoreCase) || opt == NULL) + { + *fLog << inf << "Star Maximum:" << endl; + *fLog << inf << " Calcultated \t " << setw(4) << fMaxCalc << " uA" + << endl; + *fLog << inf << " Fitted \t " << setw(4) << fMaxFit << " uA" << endl; + *fLog << inf << " CGFitted \t " << setw(4) << fMaxCGFit << " uA" << endl; + } + + if (o.Contains("pos", TString::kIgnoreCase) || opt == NULL) + { + *fLog << inf << "Star position:" << endl; + *fLog << inf << " Expected \t X " << setw(4) << fXExp + << " mm \tY " << setw(4) << fYExp << " mm" << endl; + *fLog << inf << " Calcultated \t X " << setw(4) << fMeanXCalc + << " mm \tY " << setw(4) << fMeanYCalc << " mm" << endl; + *fLog << inf << " Fitted \t X " << setw(4) << fMeanXFit + << " mm \tY " << setw(4) << fMeanYFit << " mm" << endl; + *fLog << inf << " CGFitted \t X " << setw(4) << fMeanXCGFit + << " mm \tY " << setw(4) << fMeanYCGFit << " mm" << endl; + } + + if (o.Contains("siz", TString::kIgnoreCase) || opt == NULL) + { + *fLog << inf << "Star size:" << endl; + *fLog << inf << " Calcultated \t X " << setw(4) << fSigmaMinorAxisCalc + << " mm \tY " << setw(4) << fSigmaMajorAxisCalc << " mm" << endl; + *fLog << inf << " Fitted \t X " << setw(4) << fSigmaMinorAxisFit + << " mm \tY " << setw(4) << fSigmaMajorAxisFit << " mm" << endl; + *fLog << inf << " CGFitted \t X " << setw(4) << fSigmaXCGFit + << " mm \tY " << setw(4) << fSigmaYCGFit << " mm \t correlation" + << setw(4) << fCorrXYCGFit << endl; + } + + if (o.Contains("chi", TString::kIgnoreCase) || opt == NULL) + { + *fLog << inf << "Star Fit Quality:" << endl; + *fLog << inf << " ChiSquare/Ndof \t " << setw(3) << fChiSquare + << "/" << fNdof << endl; + + *fLog << inf << "Star CGFit Quality:" << endl; + *fLog << inf << " ChiSquareCGFit/NdofCGFit \t " << setw(3) + << fChiSquareCGFit << "/" << fNdofCGFit << endl; + } + + if (o.Contains("err", TString::kIgnoreCase) || opt == NULL) + { + *fLog << inf << "CGFit Error Matrix of (fMeanXCGFit,fMeanYCGFit) :" + << endl; + *fLog << inf << " xx,xy,yy \t " << setw(3) << fXXErrCGFit << ", " + << fXYErrCGFit << ", " << fYYErrCGFit << endl; + } + + +} +//-------------------------------------------------------------------------- + + + + + + + + + Index: /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MStarLocalPos.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MStarLocalPos.h (revision 5120) +++ /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MStarLocalPos.h (revision 5120) @@ -0,0 +1,142 @@ +#ifndef MARS_MStarLocalPos +#define MARS_MStarLocalPos + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +class MStarLocalPos : public MParContainer +{ +private: + + //Expected position on camera + Float_t fMagExp; + Float_t fXExp; //[mm] + Float_t fYExp; //[mm] + + //Ideal position on camera + Float_t fMagIdeal; + Float_t fXIdeal; //[mm] + Float_t fYIdeal; //[mm] + + //Info from calculation + Float_t fMagCalc; + Float_t fMaxCalc; //[uA] + Float_t fMeanXCalc; //[mm] + Float_t fMeanYCalc; //[mm] + Float_t fSigmaMinorAxisCalc; //[mm] + Float_t fSigmaMajorAxisCalc; //[mm] + + //Info from uncorrelated Gauss fit + Float_t fMagFit; + Float_t fMaxFit; //[uA] + Float_t fMeanXFit; //[mm] + Float_t fMeanYFit; //[mm] + Float_t fSigmaMinorAxisFit; //[mm] + Float_t fSigmaMajorAxisFit; //[mm] + Float_t fXXErr; + Float_t fXYErr; + Float_t fYYErr; + + Float_t fChiSquare; + Int_t fNdof; + + //Info from correlated Gauss fit + Float_t fMagCGFit; + Float_t fMaxCGFit; //[uA] + Float_t fMeanXCGFit; //[mm] + Float_t fMeanYCGFit; //[mm] + Float_t fSigmaXCGFit; //[mm] + Float_t fSigmaYCGFit; //[mm] + Float_t fCorrXYCGFit; // correlation coefficient + Float_t fXXErrCGFit; // error matrix of (fMeanXCGFit,fMeanYCGFit) + Float_t fXYErrCGFit; + Float_t fYYErrCGFit; + + Float_t fChiSquareCGFit; + Int_t fNdofCGFit; + + +public: + + MStarLocalPos(const char *name=NULL, const char *title=NULL); + //~MStarLocalPos(); + + Float_t GetMagExp() {return fMagExp;} + Float_t GetXExp() {return fXExp;} + Float_t GetYExp() {return fYExp;} + + Float_t GetMagIdeal() {return fMagIdeal;} + Float_t GetXIdeal() {return fXIdeal;} + Float_t GetYIdeal() {return fYIdeal;} + + Float_t GetMagCalc() {return fMagCalc;} + Float_t GetMaxCalc() {return fMaxCalc;} + Float_t GetMeanXCalc() {return fMeanXCalc;} + Float_t GetMeanYCalc() {return fMeanYCalc;} + Float_t GetSigmaMinorAxisCalc() {return fSigmaMinorAxisCalc;} + Float_t GetSigmaMajorAxisCalc() {return fSigmaMajorAxisCalc;} + + Float_t GetMagFit() {return fMagFit;} + Float_t GetMaxFit() {return fMaxFit;} + Float_t GetMeanXFit() {return fMeanXFit;} + Float_t GetMeanYFit() {return fMeanYFit;} + Float_t GetSigmaMinorAxisFit() {return fSigmaMinorAxisFit;} + Float_t GetSigmaMajorAxisFit() {return fSigmaMajorAxisFit;} + Float_t GetChiSquare() {return fChiSquare;} + UInt_t GetNdof() {return fNdof;} + Float_t GetChiSquareNdof() {return fChiSquare/fNdof;} + + Float_t GetMeanX() {return fMeanXFit!=0?fMeanXFit:fMeanXCalc;} + Float_t GetMeanY() {return fMeanXFit!=0?fMeanYFit:fMeanYCalc;} + Float_t GetSigmaMinorAxis() {return fSigmaMinorAxisFit!=0?fSigmaMinorAxisFit:fSigmaMinorAxisCalc;} + Float_t GetSigmaMajorAxis() {return fSigmaMajorAxisFit!=0?fSigmaMajorAxisFit:fSigmaMajorAxisCalc;} + + // getters for the correlated Gauss fit + Float_t GetMagCGFit() {return fMagCGFit;} + Float_t GetMaxCGFit() {return fMaxCGFit;} + Float_t GetMeanXCGFit() {return fMeanXCGFit;} + Float_t GetMeanYCGFit() {return fMeanYCGFit;} + Float_t GetSigmaXCGFit() {return fSigmaXCGFit;} + Float_t GetSigmaYCGFit() {return fSigmaYCGFit;} + Float_t GetCorrXYCGFit() {return fCorrXYCGFit;} + Float_t GetXXErrCGFit() {return fXXErrCGFit;} + Float_t GetXYErrCGFit() {return fXYErrCGFit;} + Float_t GetYYErrCGFit() {return fYYErrCGFit;} + Float_t GetChiSquareCGFit() {return fChiSquareCGFit;} + UInt_t GetNdofCGFit() {return fNdofCGFit;} + Float_t GetChiSquareNdofCGFit() {return fChiSquareCGFit/fNdofCGFit;} + + + void Reset(); + + void SetExpValues(Float_t mag, Float_t x, Float_t y); + + void SetIdealValues(Float_t mag, Float_t x, Float_t y); + + void SetCalcValues(Float_t mag, Float_t max, Float_t x, Float_t y, + Float_t sigmaMinorAxis, Float_t sigmaMajorAxis); + + void SetFitValues(Float_t mag, Float_t max, Float_t x, Float_t y, + Float_t sigmaMinorAxis, Float_t sigmaMajorAxis, + Float_t chi, Int_t ndof); + + void SetFitValues(Float_t mag, Float_t max, Float_t x, Float_t y, + Float_t sigmaX, Float_t sigmaY, + Float_t chi, Int_t ndof, + Float_t xx, Float_t xy, Float_t yy); + + void SetCGFitValues(Float_t mag, Float_t max, Float_t x, Float_t y, + Float_t sigmaX, Float_t sigmaY, Float_t correlation, + Float_t xx, Float_t xy, Float_t yy, + Float_t chi, Int_t ndof); + + void Paint(Option_t *opt=NULL); + void Print(Option_t *opt=NULL) const; + + ClassDef(MStarLocalPos, 1) // Container that holds the star information in the PMT camera +}; + +#endif + + Index: /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MTelAxisFromStars.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MTelAxisFromStars.cc (revision 5120) +++ /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MTelAxisFromStars.cc (revision 5120) @@ -0,0 +1,1268 @@ +/* ======================================================================== *\ +! +! * +! * 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): Wolfgang Wittek 07/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MTelAxisFromStars +// +// This task +// - determines the transformation from expected positions of stars +// in the camera to measured positions of these stars in the camera +// - applies this transformation to expected positions of other objects +// to obtain the estimated positions of these objects in the camera +// - puts the estimated positions into the relevant containers +// +// Input Containers : +// MStarCam[MStarCam], MStarCamSource[MStarCam] +// +// Output Containers : +// MSkyCamTrans, MSrcPosCam +// +///////////////////////////////////////////////////////////////////////////// +#include +#include +#include + +#include + +#include "MTelAxisFromStars.h" + +#include "MParList.h" + + +#include "MLog.h" +#include "MLogManip.h" + +#include "MReportDrive.h" +#include "MPointingPos.h" +#include "MSrcPosCam.h" + +#include "MStarCam.h" +#include "MStarPos.h" +#include "MSkyCamTrans.h" +#include "MStarCamTrans.h" + +ClassImp(MTelAxisFromStars); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Constructor +// +MTelAxisFromStars::MTelAxisFromStars(const char *name, const char *title) +{ + fName = name ? name : "MTelAxisFromStars"; + fTitle = title ? title : "Calculate source position from star positions"; + + // if scale factor fLambda should NOT be fixed set fFixdScaleFactor to + // -1.0; otherwise set it to the value requested + fFixedScaleFactor = 1.0; + + // if rotation angle fAlfa should NOT be fixed set fFixdRotationAngle to + // -1.0; otherwise set it to the requested value + fFixedRotationAngle = 0.0; + + // default type of input is : the result of the Gauss fit + // type 0 : result from the weighted average + // type 1 : result from the Gauss fit + fInputType = 1; + + // default value of fAberr + // the value 1.07 is valid if the expected position (with aberration) + // in the camera is calculated as the average of the positions obtained + // from the reflections at the individual mirrors + fAberr = 1.07; + +} + +// -------------------------------------------------------------------------- +// +// Destructor +// +MTelAxisFromStars::~MTelAxisFromStars() +{ + delete fStarCamTrans; +} + +// -------------------------------------------------------------------------- +// +// Set links to containers +// + +Int_t MTelAxisFromStars::PreProcess(MParList *pList) +{ + fDrive = (MReportDrive*)pList->FindObject(AddSerialNumber("MReportDrive")); + if (!fDrive) + { + *fLog << err << AddSerialNumber("MReportDrive") + << " not found... aborting." << endl; + return kFALSE; + } + + + fStarCam = (MStarCam*)pList->FindObject("MStarCam", "MStarCam"); + if (!fStarCam) + { + *fLog << err << "MStarCam not found... aborting." << endl; + return kFALSE; + } + + + fSourceCam = (MStarCam*)pList->FindObject("MSourceCam", "MStarCam"); + if (!fSourceCam) + { + *fLog << warn << "MSourceCam[MStarCam] not found... continue " << endl; + } + + + //----------------------------------------------------------------- + fSrcPos = (MSrcPosCam*)pList->FindCreateObj(AddSerialNumber("MSrcPosCam"), "MSrcPosCam"); + if (!fSrcPos) + return kFALSE; + + fPntPos = (MSrcPosCam*)pList->FindCreateObj(AddSerialNumber("MSrcPosCam"), "MPntPosCam"); + if (!fPntPos) + return kFALSE; + + fPointPos = (MPointingPos*)pList->FindCreateObj(AddSerialNumber("MPointingPos"), "MPointingPos"); + if (!fPointPos) + return kFALSE; + + fSourcePos = (MPointingPos*)pList->FindCreateObj(AddSerialNumber("MPointingPos"), "MSourcePos"); + if (!fSourcePos) + return kFALSE; + + fSkyCamTrans = (MSkyCamTrans*)pList->FindCreateObj(AddSerialNumber("MSkyCamTrans")); + if (!fSkyCamTrans) + return kFALSE; + + + //----------------------------------------------------------------- + // book an MStarCamTrans object + // this is needed when calling one of the member functions of MStarCamTrans + + MGeomCam *geom = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam")); + if (!geom) + { + *fLog << err << AddSerialNumber("MGeomCam") + << " not found... aborting." << endl; + return kFALSE; + } + + MObservatory *obs = (MObservatory*)pList->FindObject(AddSerialNumber("MObservatory")); + if (!obs) + { + *fLog << err << AddSerialNumber("MObservatory") + << " not found... aborting." << endl; + return kFALSE; + } + + //----------------------------------------------------------------- + fStarCamTrans = new MStarCamTrans(*geom, *obs); + + *fLog << all << "MTelAxisFromStars::Preprocess; the optical aberration factor is set equal to : " + << fAberr ; + + *fLog << all << "MTelAxisFromStars::Preprocess; input type is set equal to : " + << fInputType ; + if (fInputType == 0) + *fLog << " (calculated star positions)" << endl; + else + *fLog << " (fitted star positions)" << endl; + + *fLog << all << "MTelAxisFromStars::Preprocess; scale factor will be fixed at : " + << fFixedScaleFactor << endl; + + *fLog << all << "MTelAxisFromStars::Preprocess; rotation angle will be fixed at : " + << fFixedRotationAngle << endl; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Set optical aberration factor +// +// fAberr is the ratio between +// the distance from the camera center with optical aberration and +// the distance from the camera center with an ideal imaging +// +// fAberr = r_real/r_ideal +// +void MTelAxisFromStars::SetOpticalAberr(Double_t aberr) +{ + fAberr = aberr; +} + +// -------------------------------------------------------------------------- +// +// Set the type of the input +// +// type = 0 calculated star positions (by averaging) +// type = 1 fitted star positions (by Gauss fit) +// +void MTelAxisFromStars::SetInputType(Int_t type) +{ + fInputType = type; +} + +// -------------------------------------------------------------------------- +// +// Fix the scale factor fLambda +// +// +void MTelAxisFromStars::FixScaleFactorAt(Double_t lambda) +{ + fFixedScaleFactor = lambda; +} + + +// -------------------------------------------------------------------------- +// +// Fix rotation angle fAlfa +// +// +void MTelAxisFromStars::FixRotationAngleAt(Double_t alfa) +{ + fFixedRotationAngle = alfa; // [degrees] +} + + +// -------------------------------------------------------------------------- +// +// Process +// +// call FindSkyCamTrans to find the Sky-Camera transformation +// call TransSkyCam to transform some sky directions +// into the camera system +// put the estimated source position into MSrcPosCam +// + +Int_t MTelAxisFromStars::Process() +{ + //Int_t run = fRun->GetRunNumber(); + //*fLog << "MTelAxisFromStars::Process; run = " << run << endl; + + //-------------------------------------- + // Define the input for FindSkyCamTrans + // + + // get the expected (axy[0], axy[1]) and the measured positions + // (bxy[0], bxy[1]) of stars in the camera from MStarCam + Int_t fNumStars = fStarCam->GetNumStars(); + + if (fNumStars <= 0){ + *fLog << err << "No stars found!" << endl; + return kTRUE; + } + + TArrayD axy[2]; + axy[0].Set(fNumStars); + axy[1].Set(fNumStars); + + TArrayD bxy[2]; + bxy[0].Set(fNumStars); + bxy[1].Set(fNumStars); + + // error matrix of bxy + TArrayD exy[2][2]; + exy[0][0].Set(fNumStars); + exy[0][1].Set(fNumStars); + exy[1][0].Set(fNumStars); + exy[1][1].Set(fNumStars); + + // transformation parameters + Double_t fLambda; + Double_t fAlfa; + Double_t fA[2][2]; + Double_t fD[2]; + Double_t fErrD[2][2]; + Int_t fNumIter; + Int_t fNdof; + Double_t fChi2; + Double_t fChi2Prob; + + MStarPos *star = 0; + TIter next(fStarCam->GetList()); + Int_t ix = 0; + + // loop over all stars + while ( (star = (MStarPos*)next()) ) + { + axy[0][ix] = star->GetXExp(); + axy[1][ix] = star->GetYExp(); + + if (fInputType == 0) + { + // values from averaging + bxy[0][ix] = star->GetMeanXCalc(); + bxy[1][ix] = star->GetMeanYCalc(); + + // this is the error matrix for (MeanXCalc, MeanYCalc); + // this is the error matrix which should be used + exy[0][0][ix] = star->GetXXErrCalc(); + exy[0][1][ix] = star->GetXYErrCalc(); + exy[1][0][ix] = star->GetXYErrCalc(); + exy[1][1][ix] = star->GetYYErrCalc(); + + //exy[0][0][ix] = star->GetSigmaXCalc()*star->GetSigmaXCalc(); + //exy[0][1][ix] = 0.0; + //exy[1][0][ix] = 0.0; + //exy[1][1][ix] = star->GetSigmaYCalc()*star->GetSigmaYCalc(); + } + + else if (fInputType == 1) + { + // values from Gauss fit + bxy[0][ix] = star->GetMeanXFit(); + bxy[1][ix] = star->GetMeanYFit(); + + // this is the error matrix for (MeanXFit, MeanYFit); + // this is the error matrix which should be used + exy[0][0][ix] = star->GetXXErrFit(); + exy[0][1][ix] = star->GetXYErrFit(); + exy[1][0][ix] = star->GetXYErrFit(); + exy[1][1][ix] = star->GetYYErrFit(); + + // this is the error matrix constructed from SigmaXFit and SigmaYFit; + // it is used because the errors above are too small, at present + //exy[0][0][ix] = star->GetSigmaXFit() * star->GetSigmaXFit(); + //exy[0][1][ix] = star->GetCorrXYFit() * + // star->GetSigmaXFit() * star->GetSigmaYFit(); + //exy[1][0][ix] = exy[0][1][ix]; + //exy[1][1][ix] = star->GetSigmaYFit() * star->GetSigmaYFit(); + } + + else + { + *fLog << err << "MTelAxisFromStars::Process; type of input is not defined" + << endl; + return kFALSE; + } + + // don't include stars with undefined error + Double_t deter = exy[0][0][ix]*exy[1][1][ix] + - exy[0][1][ix]*exy[1][0][ix]; + + //*fLog << "ix ,deter, xx, xy, yy = " << ix << ": " + // << deter << ", " << exy[0][0][ix] << ", " + // << exy[0][1][ix] << ", " << exy[1][1][ix] << endl; + if (deter <= 0.0) + continue; + + //*fLog << "MTelAxisFromStars : " << endl; + //*fLog << " ix, XExp, YExp, XFit, YFit, SigmaX2, SigmaXY, SigmaY2 = " + // << ix << " : " + // << axy[0][ix] << ", " << axy[1][ix] << ", " + // << bxy[0][ix] << ", " << bxy[1][ix] << ", " + // << exy[0][0][ix] << ", " << exy[0][1][ix] << ", " + // << exy[1][1][ix] << endl; + + ix++; + } + + //-------------------------------------- + // Find the transformation from expected positions (axy[1], axy[2]) + // to measured positions (bxy[1], bxy[2]) in the camera + + Int_t fNStars = ix; + + if (ix < fNumStars) + { + // reset the sizes of the arrays + Int_t fNStars = ix; + axy[0].Set(fNStars); + axy[1].Set(fNStars); + + bxy[0].Set(fNStars); + bxy[1].Set(fNStars); + + exy[0][0].Set(fNStars); + exy[0][1].Set(fNStars); + exy[1][0].Set(fNStars); + exy[1][1].Set(fNStars); + } + + Bool_t fitOK; + if (fNStars < 1) + { + //*fLog << "MTelAxisFromStars::Process; no star for MTelAxisFromStars" + // << endl; + fitOK = kFALSE; + } + else + { + fitOK = FindSkyCamTrans(axy, bxy, exy, + fFixedRotationAngle, fFixedScaleFactor, fLambda, + fAlfa , fA, fD, fErrD, + fNumIter, fNdof, fChi2, fChi2Prob); + } + + if (!fitOK && fNStars >= 1) + { + *fLog << err + << "MTelAxisFromStars::Process; Fit to find transformation from star to camera system failed" + << endl; + + //if (fNStars > 0) + //{ + // *fLog << err + // << " fNumIter, fNdof, fChi2, fChi2Prob = " << fNumIter + // << ", " << fNdof << ", " << fChi2 << ", " << fChi2Prob << endl; + //} + + return kTRUE; + } + + + //-------------------------------------- + // Put the transformation parameters into the MSkyCamTrans container + + fSkyCamTrans->SetParameters(fLambda, fAlfa, fA, fD, fErrD, + fNumStars, fNumIter, fNdof, fChi2, fChi2Prob); + fSkyCamTrans->SetReadyToSave(); + + //-------------------------------------- + // Put the camera position (X, Y) + // obtained by transforming the camera center (0, 0) + // into MPntPosCam[MSrcPosCam] + + fPntPos->SetXY(fD[0], fD[1]); + fPntPos->SetReadyToSave(); + + + //-------------------------------------- + // Put the corrected pointing position into MPointingPos + // + SetPointingPosition(fStarCamTrans, fDrive, fSkyCamTrans, fPointPos); + + + //-------------------------------------- + // Put the estimated position of the source into SrcPosCam + // + // get the source direction from MReportDrive + // Note : this has to be changed for the wobble mode, where the source + // isn't in the center of the camera + Double_t decsource = fDrive->GetDec(); + Double_t rasource = fDrive->GetRa(); + // + Double_t radrive = fDrive->GetRa(); + Double_t hdrive = fDrive->GetHa(); + Double_t hsource = hdrive + radrive - rasource; + fSourcePos->SetSkyPosition(rasource, decsource, hsource); + + SetSourcePosition(fStarCamTrans, fPointPos, fSourcePos, fSrcPos); + + *fLog << "after calling SetSourcePosition : , X, Y ,fD = " + << fSrcPos->GetX() << ", " << fSrcPos->GetY() << ", " + << fD[0] << ", " << fD[1] << endl; + + //-------------------------------------- + // Apply the transformation to some expected positions (asxy[1], asxy[2]) + // to obtain estimated positions (bsxy[1], bsxy[2]) in the camera + // and their error matrices esxy[2][2] + + // get the expected positions (asxy[1], asxy[2]) from another MStarCam + // container (with the name "MSourceCam") + Int_t fNumStarsSource = 0; + + if (fSourceCam) + fNumStarsSource = fSourceCam->GetNumStars(); + + //*fLog << "MTelAxisFromStars::Process; fNumStarsSource = " + // << fNumStarsSource << endl; + + if (fNumStarsSource > 0) + { + TArrayD asxy[2]; + asxy[0].Set(fNumStarsSource); + asxy[1].Set(fNumStarsSource); + + TArrayD bsxy[2]; + bsxy[0].Set(fNumStarsSource); + bsxy[1].Set(fNumStarsSource); + + TArrayD esxy[2][2]; + esxy[0][0].Set(fNumStarsSource); + esxy[0][1].Set(fNumStarsSource); + esxy[1][0].Set(fNumStarsSource); + esxy[1][1].Set(fNumStarsSource); + + MStarPos *starSource = 0; + TIter nextSource(fSourceCam->GetList()); + ix = 0; + while ( (starSource = (MStarPos*)nextSource()) ) + { + asxy[0][ix] = starSource->GetXExp(); + asxy[1][ix] = starSource->GetYExp(); + + ix++; + } + + TransSkyCam(fLambda, fA, fD, fErrD, asxy, bsxy, esxy); + + // put the estimated positions into the MStarCam container + // with name "MSourceCam" + TIter setnextSource(fSourceCam->GetList()); + ix = 0; + while ( (starSource = (MStarPos*)setnextSource()) ) + { + Double_t corr = esxy[0][1][ix] / + TMath::Sqrt( esxy[0][0][ix] * esxy[1][1][ix] ); + if (fInputType == 1) + { + starSource->SetFitValues(100.0, 100.0, bsxy[0][ix], bsxy[1][ix], + TMath::Sqrt(esxy[0][0][ix]), TMath::Sqrt(esxy[1][1][ix]), corr, + esxy[0][0][ix], esxy[0][1][ix], esxy[1][1][ix], + fChi2, fNdof); + } + else + { + starSource->SetCalcValues(100.0, 100.0, bsxy[0][ix], bsxy[1][ix], + TMath::Sqrt(esxy[0][0][ix]), TMath::Sqrt(esxy[1][1][ix]), corr, + esxy[0][0][ix], esxy[0][1][ix], esxy[1][1][ix]); + } + + ix++; + } + + } + + //-------------------------------------- + + return kTRUE; +} + + + +//--------------------------------------------------------------------------- +// +// FindSkyCamTrans +// +// This routine determines the transformation +// +// ( cos(alfa) -sin(alfa) ) +// b = lambda * A * a + d A = ( ) +// ^ ^ ^ ( sin(alfa) cos(alfa) ) +// | | | +// scale rotation shift +// factor matrix +// +// from sky coordinates 'a' (projected onto the camera) to camera +// coordinates 'b', using the positions of known stars in the camera. +// The latter positions may have been determined by analysing the +// DC currents in the different pixels. +// +// Input : a[2] x and y coordinates of stars projected onto the camera; +// they were obtained from (RA, dec) of the stars and +// (ThetaTel, PhiTel) and the time of observation; +// these are the 'expected positions' of stars in the camera +// b[2] 'measured positions' of these stars in the camera; +// they may have been obtained from the DC currents +// e[2][2] error matrix of b[2] +// fixedrotationangle value [in degrees] at which rotation angle +// alfa should be fixed; -1 means don't fix +// fixedscalefactor value at which scale factor lambda +// should be fixed; -1 means don't fix +// +// Output : lambda, alfadeg, A[2][2], d[2] fit results; +// parameters describing the transformation +// from 'expected positions' to the 'measured +// positions' in the camera +// errd[2][2] error matrix of d[2] +// fNumIter number of iterations +// fNdoF number of degrees of freedom +// fChi2 chi-square value +// fChi2Prob chi-square probability +// +// The units are assumed to be +// [degrees] for alfadeg +// [mm] for a, b, d +// [1] for lambda + +Bool_t MTelAxisFromStars::FindSkyCamTrans( + TArrayD a[2], TArrayD b[2], TArrayD e[2][2], + Double_t &fixedrotationang, Double_t &fixedscalefac, Double_t &lambda, + Double_t &alfadeg, Double_t A[2][2], Double_t d[2], Double_t errd[2][2], + Int_t &fNumIter, Int_t &fNdof, Double_t &fChi2, Double_t &fChi2Prob) +{ + Int_t fNumStars = a[0].GetSize(); + + //*fLog << "MTelAxisFromStars::FindSkyCamTrans; expected and measured positions :" + // << endl; + for (Int_t ix=0; ix 1.1 ) + lambda = 1.1; + + dlambda = lambda - lambdaold; + lambdaold = lambda; + + //*fLog << "num, denom, lambda, dlambda = " << num + // << ", " << denom << ", " << lambda << ", " + // << dlambda << endl; + } + + } + //------- end of iteration ***************************************** + + alfadeg = alfa * TMath::RadToDeg(); + + // calculate error matrix of d[2] + errd[0][0] = sumcinv[0][0]; + errd[0][1] = sumcinv[0][1]; + errd[1][0] = sumcinv[1][0]; + errd[1][1] = sumcinv[1][1]; + + // evaluate quality of fit + + // calculate chi2 + for (Int_t ix=0; ixGetDec(); + Double_t hdrive = fdrive->GetHa(); + Double_t radrive = fdrive->GetRa(); + + // this is the estimated position (with optical aberration) in the camera + // corresponding to the direction in MReportDrive + Double_t Xpoint = (ftrans->GetShiftD())[0]; + Double_t Ypoint = (ftrans->GetShiftD())[1]; + + // get the sky direction corresponding to the position (0,0) in the camera + Double_t decpoint = 0.0; + Double_t hpoint = 0.0; + fstarcamtrans->CelCamToCel0(decdrive, hdrive, + Xpoint/fAberr, Ypoint/fAberr, decpoint, hpoint); + Double_t rapoint = radrive - hpoint + hdrive; + fpointpos->SetSkyPosition(rapoint, decpoint, hpoint); + + // get the local direction corresponding to the position (0,0) in the camera + Double_t thetadrive = fdrive->GetNominalZd(); + Double_t phidrive = fdrive->GetNominalAz(); + Double_t thetapoint = 0.0; + Double_t phipoint = 0.0; + fstarcamtrans->LocCamToLoc0(thetadrive, phidrive, + Xpoint/fAberr, Ypoint/fAberr, thetapoint, phipoint); + fpointpos->SetLocalPosition(thetapoint, phipoint); + fpointpos->SetReadyToSave(); + + //*fLog << "SetPointingPosition : decdrive, hdrive, radrive Xpoint, Ypoint = " + // << decdrive << ", " << hdrive << ", " << radrive << ", " + // << Xpoint << ", " << Ypoint << endl; + + //*fLog << "SetPointingPosition : thetadrive, phidrive, thetapoint, phipoint = " + // << thetadrive << ", " << phidrive << ", " << thetapoint << ", " + // << phipoint << endl; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// SetSourcePosition +// +// put the estimated position of the source in the camera into +// MSrcPosCam[MSrcPosCam] +// +// and the estimated local direction of the source into +// MSourcePos[MPointingPos] +// + +Bool_t MTelAxisFromStars::SetSourcePosition(MStarCamTrans *fstarcamtrans, + MPointingPos *fpointpos, MPointingPos *fsourcepos, MSrcPosCam *fsrcpos) +{ + // get the corrected pointing direction + // corresponding to the position (0,0) in the camera + Double_t decpoint = fpointpos->GetDec(); + Double_t hpoint = fpointpos->GetHa(); + + // get the sky direction of the source + Double_t decsource = fsourcepos->GetDec(); + Double_t hsource = fsourcepos->GetHa(); + + // get the estimated position (Xsource, Ysource) of the source in the camera; + // this is a position for an ideal imaging, without optical aberration + Double_t Xsource = 0.0; + Double_t Ysource = 0.0; + fstarcamtrans->Cel0CelToCam(decpoint, hpoint, + decsource, hsource, Xsource, Ysource); + fsrcpos->SetXY(Xsource*fAberr, Ysource*fAberr); + fsrcpos->SetReadyToSave(); + + // get the estimated local direction of the source + Double_t thetapoint = fpointpos->GetZd(); + Double_t phipoint = fpointpos->GetAz(); + Double_t thetasource = 0.0; + Double_t phisource = 0.0; + fstarcamtrans->Loc0CamToLoc(thetapoint, phipoint, + Xsource, Ysource, thetasource, phisource); + fsourcepos->SetLocalPosition(thetasource, phisource); + fsourcepos->SetReadyToSave(); + + //*fLog << "SetSourcePosition : decpoint, hpoint, decsource, hsource, Xsource, Ysource = " + // << decpoint << ", " << hpoint << ", " << decsource << ", " + // << hsource << ", " << Xsource << ", " << Ysource << endl; + //*fLog << "SetSourcePosition : thetapoint, phipoint, thetasource, phisource = " + // << thetapoint << ", " << phipoint << ", " << thetasource << ", " + // << phisource << endl; + + return kTRUE; +} + +// -------------------------------------------------------------------------- + + + + + + + + + + + + Index: /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MTelAxisFromStars.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MTelAxisFromStars.h (revision 5120) +++ /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MTelAxisFromStars.h (revision 5120) @@ -0,0 +1,88 @@ +#ifndef MARS_MTelAxisFromStars +#define MARS_MTelAxisFromStars + +///////////////////////////////////////////////////////////////////////////// +// // +// MTelAxisFromStars // +// // +// Task to calculate the position (in the camera) +// of certain sky directions (source position, tracking direction, ...) +// from the positions (in the camera) of known stars +// // +///////////////////////////////////////////////////////////////////////////// + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#ifndef ROOT_TArrayD +#include +#endif + +class MReportDrive; +class MPointingPos; +class MSrcPosCam; +class MStarCam; +class MSkyCamTrans; +class MStarCamTrans; + +class MTelAxisFromStars : public MTask +{ + private: + + const MStarCam *fStarCam; //! + const MStarCam *fSourceCam; //! + + MReportDrive *fDrive; //! + MPointingPos *fPointPos; //! + MPointingPos *fSourcePos; //! + MSrcPosCam *fSrcPos; //! + MSrcPosCam *fPntPos; //! + MSkyCamTrans *fSkyCamTrans; //! + MStarCamTrans *fStarCamTrans; //! + + Double_t fFixedScaleFactor; //! + Double_t fFixedRotationAngle; //! [degrees] + Double_t fAberr; //! optical aberration factor + + Int_t fInputType; //! type of input + + Int_t PreProcess(MParList *pList); + Int_t Process(); + + public: + + MTelAxisFromStars(const char *name=NULL, const char *title=NULL); + ~MTelAxisFromStars(); + + void FixScaleFactorAt(Double_t lambda = 1.0); + void FixRotationAngleAt(Double_t alfa = 0.0); // alfa in [degrees] + + void SetOpticalAberr(Double_t aberr); + void SetInputType(Int_t type = 1); + + Bool_t FindSkyCamTrans(TArrayD[2], TArrayD[2], TArrayD[2][2], + Double_t &, Double_t &, Double_t &, + Double_t &, Double_t[2][2] , Double_t[2], Double_t[2][2], + Int_t &, Int_t &, Double_t &, Double_t &); + + void TransSkyCam(Double_t &, Double_t[2][2], Double_t[2], Double_t[2][2], + TArrayD[2], TArrayD[2], TArrayD[2][2]); + + Bool_t SetPointingPosition(MStarCamTrans *fstarcamtrans, + MReportDrive *fdrive, MSkyCamTrans *ftrans, MPointingPos *fpointpos); + Bool_t SetSourcePosition(MStarCamTrans *fstarcamtrans, + MPointingPos *fpointpos, MPointingPos *fsourcepos, MSrcPosCam *fsrcpos); + + ClassDef(MTelAxisFromStars, 0) // Task to calculate the source position from star positions +}; + +#endif + + + + + + + + Index: /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MVPObject.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MVPObject.cc (revision 5120) +++ /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MVPObject.cc (revision 5120) @@ -0,0 +1,280 @@ +/* ======================================================================== *\ +! +! * +! * 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 +// "-", eg. "1959+650" +// +void MVPObject::SetObjectByName(const char *object) +{ + fObjectName=object; + fGotName=kTRUE; + + Int_t ra, dec; + sscanf(object, "%d%d", &ra, &dec); + + fRA = fgHrsToRad * (0.01*ra + (r%100)/60.); + fDec = fgDegToRad * 0.1 * dec; + + fGotRA=kTRUE; + fGotDec=kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Sets RA position of object. Position is to be provided in hours, minutes, +// seconds, and microseconds (if needed) +// +void MVPObject::SetRA(Int_t rh, Int_t rm, Int_t rs, Int_t ru) +{ + // Rect is a timelike value... + fRA = fgHrsToRad*((Double_t)rh + (Double_t)rm/60 + (Double_t)rs/(60*60) + (Double_t)ru/(36000)); + fBody = 10; + fGotRA = kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Sets Dec position of object. Position is to be provided in degrees, +// minutes, seconds, and microseconds (if needed) +// +void MVPObject::SetDec(Int_t dh, Int_t dm, Int_t ds, Int_t du) +{ + // Dec is an anglelike value + fDec = fgDegToRad*((Double_t)dh + (Double_t)dm/60 + (Double_t)ds/(60*60) + (Double_t)du/(36000)); + fBody = 10; + fGotDec = kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Alternatively to providing RA, Dec and Name of an object, you may provide +// a solar system object (which has no fixed RA, Dec, by the way!) with +// MVPObject::SetObject. +// - +// UInt_t body | Object Sun and Moon will be objects needed at most, +// 0 | Sun presumably. +// 1 | Mercury +// 2 | Venus +// 3 | Moon +// 4 | Mars +// 5 | Jupiter +// 6 | Saturn +// 7 | Uranus +// 8 | Neptune +// 9 | Pluto +// +Bool_t MVPObject::SetObject(UInt_t body) +{ + if (body > 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/mpisa/classes/mispointing/MVPObject.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MVPObject.h (revision 5120) +++ /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MVPObject.h (revision 5120) @@ -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/mpisa/classes/mispointing/MVPPlotter.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MVPPlotter.cc (revision 5120) +++ /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MVPPlotter.cc (revision 5120) @@ -0,0 +1,587 @@ +/* ======================================================================== *\ +! +! * +! * 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]; + memset(visibility, 0, 13*18*sizeof(Float_t)); + /* + 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 = TMath::Min(distance, fMars->GetDistance(fObject)); + distance = TMath::Min(distance, fJupiter->GetDistance(fObject)); + distance = TMath::Min(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)); + } + printf("\n"); + } + + int vistimestart=0; + int vistimeend=0; + 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) { + 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++; + + const Int_t startday2 = fTime->MJDStartOfYear(year); + const 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/30; + maxza += 80; + + } else { //colors=yellow to red + maxza -= 30; + maxza *= 11/60; + maxza += 89; + } + + 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=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 *= 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+=pow(10.,(6.15 - rho/40.)); /* eqn 21 */ + else if (fabs(rho) > 0.25) + 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/mpisa/classes/mispointing/MVPPlotter.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MVPPlotter.h (revision 5120) +++ /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MVPPlotter.h (revision 5120) @@ -0,0 +1,77 @@ +#ifndef MARS_MVPPlotter +#define MARS_MVPPlotter + +#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/mpisa/classes/mispointing/MVPTime.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MVPTime.cc (revision 5120) +++ /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MVPTime.cc (revision 5120) @@ -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/mpisa/classes/mispointing/MVPTime.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MVPTime.h (revision 5120) +++ /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MVPTime.h (revision 5120) @@ -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/mpisa/classes/mispointing/Makefile =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/Makefile (revision 5120) +++ /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/Makefile (revision 5120) @@ -0,0 +1,57 @@ +################################################################## +# +# subdirectory makefile +# +# for the MARS software +# +################################################################## +include ../../../../Makefile.conf.$(OSTYPE) +include ../../../../Makefile.conf.general + +#------------------------------------------------------------------------------ + +CINT = Pisa + +#------------------------------------------------------------------------------ + +INCLUDES = -I. \ + -I../../../../mbase \ + -I../../../../mjobs \ + -I../../../../mpedestal \ + -I../../../../mbadpixels \ + -I../../../../mfileio \ + -I../../../../mraw \ + -I../../../../manalysis \ + -I../../../../mgui \ + -I../../../../mgeom \ + -I../../../../msignal \ + -I../../../../mcalib \ + -I../../../../mfilter \ + -I../../../../mhbase \ + -I../../../../mimage \ + -I../../../../mpointing \ + -I../../../../mcamera \ + -I../../../../mhist \ + -I../../../../mmc \ + -I../../../../mreport \ + -I../../../../mastro \ + -I../../../../mstarcam + +SRCFILES = \ + ./MFindStars.cc \ + ./MTelAxisFromStars.cc \ + ./MHTelAxisFromStars.cc \ + ./MSkyCamTrans.cc \ + ./MSourceDirections.cc \ + ./MStarLocalCam.cc \ + ./MStarLocalPos.cc + + +############################################################ + +all: $(OBJS) + +include ../../../../Makefile.rules + +clean: rmcint rmobjs rmcore rmlib +mrproper: clean rmbak Index: /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MispointingLinkDef.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MispointingLinkDef.h (revision 5120) +++ /trunk/MagicSoft/Mars/mtemp/mpisa/classes/mispointing/MispointingLinkDef.h (revision 5120) @@ -0,0 +1,15 @@ +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class MFindStars+; +#pragma link C++ class MTelAxisFromStars+; +#pragma link C++ class MHTelAxisFromStars+; +#pragma link C++ class MSkyCamTrans+; +#pragma link C++ class MSourceDirections+; +#pragma link C++ class MStarLocalCam+; +#pragma link C++ class MStarLocalPos; + +#endif