Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 3208) +++ /trunk/MagicSoft/Mars/Changelog (revision 3209) @@ -4,4 +4,15 @@ -*-*- END OF LINE -*-*- + + + + 2004/02/17: Wolfgang Wittek + + * manalysis/MSourcPosfromStarPos.[h,cc] + - new class; calculates the position of the source in the camera from + the position of a star in the camera and puts it into + MSrcPosCam + + 2004/02/16: Abelardo Moralejo Index: /trunk/MagicSoft/Mars/macros/ONOFFAnalysis.C =================================================================== --- /trunk/MagicSoft/Mars/macros/ONOFFAnalysis.C (revision 3208) +++ /trunk/MagicSoft/Mars/macros/ONOFFAnalysis.C (revision 3209) @@ -172,5 +172,5 @@ //----------------------------------------------- //const char *offfile = "~magican/ct1test/wittek/offdata.preproc"; - const char *offfile = "/data/MAGIC/rootdata/2003_11_29/20031128_032*_D_OffCrab1_E.root"; + const char *offfile = "20040126_OffCrab_"; //const char *onfile = "~magican/ct1test/wittek/mkn421_on.preproc"; @@ -203,5 +203,5 @@ // - write root file for ON (or OFF or MC) data (ON1.root, ...); - Bool_t JobA = kFALSE; + Bool_t JobA = kTRUE; Bool_t GPad = kFALSE; // generate padding histograms? Bool_t WPad = kFALSE; // write out padding histograms ? @@ -233,5 +233,5 @@ // - update the input files with the hadroness (==>ON3.root or MC3.root) - Bool_t JobB_SC_UP = kTRUE; + Bool_t JobB_SC_UP = kFALSE; Bool_t CMatrix = kFALSE; // create training and test matrices Bool_t RMatrix = kTRUE; // read training and test matrices from file @@ -239,5 +239,5 @@ // and write supercuts parameter values // onto the file parSCfile - Bool_t RTest = kTRUE; // test the supercuts using the test matrix + Bool_t RTest = kFALSE; // test the supercuts using the test matrix Bool_t WSC = kFALSE; // update input root file ? @@ -313,10 +313,20 @@ // for generating the histograms to be used in the padding + TString fileON = inPath; fileON += onfile; fileON += "CalibratedEvts"; fileON += ".root"; - TString fileOFF = offfile; + + TString fileOFF = inPath; + fileOFF += offfile; + fileOFF += "CalibratedEvts"; + fileOFF += ".root"; + TString fileMC = mcfile; + fileMC += mcfile; + fileMC += "CalibratedEvts"; + fileMC += ".root"; + gLog << "fileON, fileOFF, fileMC = " << fileON << ", " << fileOFF << ", " << fileMC << endl; @@ -329,6 +339,6 @@ //-------------------------------------------------- // type of data to be padded - TString typeInput = "ON"; - //TString typeInput = "OFF"; + //TString typeInput = "ON"; + TString typeInput = "OFF"; //TString typeInput = "MC"; gLog << "typeInput = " << typeInput << endl; @@ -345,6 +355,13 @@ // name of output root file + if (typeInput == "ON") + TString file(onfile); + else if (typeInput == "OFF") + TString file(offfile); + else if (typeInput == "MC") + TString file(mcfile); + TString outNameImage = outPath; - outNameImage += onfile; + outNameImage += file; outNameImage += "Hillas"; outNameImage += typeInput; @@ -1765,6 +1782,6 @@ TString parSCinit = outPath; - parSCinit += "parSC_060204a"; - //parSCinit = ""; + //parSCinit += "parSC_060204"; + parSCinit = ""; if (parSCinit != "") @@ -1785,5 +1802,5 @@ TString parSCfile = outPath; - parSCfile += "parSC_060204b"; + parSCfile += "parSC_130204a"; gLog << "parSCfile = " << parSCfile << endl; @@ -1946,4 +1963,5 @@ if (parSCinit == "") { + /* Double_t vparams[104] = { // LengthUp @@ -2027,6 +2045,90 @@ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }; - - + */ + + + Double_t vparams[104] = { + // LengthUp + 0.315585, 0.001455, 0.203198, 0.005532, -0.001670, -0.020362, + 0.007388, -0.013463, + // LengthLo + 0.151530, 0.028323, 0.510707, 0.053089, 0.013708, 2.357993, + 0.000080, -0.007157, + // WidthUp + 0.145412, -0.001771, 0.054462, 0.022280, -0.009893, 0.056353, + 0.020711, -0.016703, + // WidthLo + 0.089187, -0.006430, 0.074442, 0.003738, -0.004256, -0.014101, + 0.006126, -0.002849, + // DistUp + 1.787943, 0.0, 2.942310, 0.199815, 0.0, 0.249909, + 0.189697, 0.0, + // DistLo + 0.589406, 0.0, -0.083964,-0.007975, 0.0, 0.045374, + -0.001750, 0.0, + // AsymUp + 1.e10, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // AsymLo + -1.e10, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // ConcUp + 1.e10, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // ConcLo + -1.e10, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // Leakage1Up + 1.e10, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // Leakage1Lo + -1.e10, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // AlphaUp + 13.12344, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0 }; + + Double_t vsteps[104] = { + // LengthUp + 0.03, 0.0002, 0.02, 0.0006, 0.0002, 0.002, + 0.0008, 0.002, + // LengthLo + 0.02, 0.003, 0.05, 0.006, 0.002, 0.3, + 0.0001, 0.0008, + // WidthUp + 0.02, 0.0002, 0.006, 0.003, 0.002, 0.006, + 0.002, 0.002, + // WidthLo + 0.009, 0.0007, 0.008, 0.0004, 0.0005, 0.002, + 0.0007, 0.003, + // DistUp + 0.2, 0.0, 0.3, 0.02, 0.0, 0.03, + 0.02, 0.0 + // DistLo + 0.06, 0.0, 0.009, 0.0008, 0.0, 0.005, + 0.0002, 0.0 + // AsymUp + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // AsymLo + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // ConcUp + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // ConcLo + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // Leakage1Up + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // Leakage1Lo + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, + // AlphaUp + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0 }; + + params.Set(104, vparams); steps.Set (104, vsteps ); Index: /trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h =================================================================== --- /trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 3208) +++ /trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 3209) @@ -80,4 +80,5 @@ #pragma link C++ class MArrivalTimeCalc2+; +#pragma link C++ class MSourcePosfromStarPos+; #pragma link C++ class MMcCalibrationUpdate+; @@ -86,2 +87,4 @@ #endif + + Index: /trunk/MagicSoft/Mars/manalysis/MSourcePosfromStarPos.cc =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MSourcePosfromStarPos.cc (revision 3209) +++ /trunk/MagicSoft/Mars/manalysis/MSourcePosfromStarPos.cc (revision 3209) @@ -0,0 +1,448 @@ +/* ======================================================================== *\ +! +! * +! * 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 02/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MSourcePosfromStarPos +// +// This is a task which +// - calculates the position of the source in the camera +// from the position of a known star in the camera +// - and puts the source position into the container MSrcPosCam +// +// Input : +// ASCII file containing for each run +// - the run number +// - the direction (theta, phi) the telescope is pointing to in [deg] +// - the position (xStar, yStar) of a known star in the camera in [mm] +// - the error (dxStar, dyStar) of this position in [mm] +// +// Output Containers : +// MSrcPosCam +// +///////////////////////////////////////////////////////////////////////////// +#include +#include + +#include + +#include "MSourcePosfromStarPos.h" + +#include "MParList.h" +#include "MRawRunHeader.h" +#include "MGeomCam.h" +#include "MSrcPosCam.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MSourcePosfromStarPos); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MSourcePosfromStarPos::MSourcePosfromStarPos( + const char *name, const char *title) + : fIn(NULL) +{ + fName = name ? name : "MSourcePosfromStarPos"; + fTitle = title ? title : "Calculate source position from star position"; + + fFileNames = new TList; + fFileNames->SetOwner(); +} + +// -------------------------------------------------------------------------- +// +// Delete the filename list and the input stream if one exists. +// +MSourcePosfromStarPos::~MSourcePosfromStarPos() +{ + delete fFileNames; + if (fIn) + delete fIn; +} + + +// -------------------------------------------------------------------------- +// +// Set the sky coordinates of the source and of the star +// +void MSourcePosfromStarPos::SetSourceAndStarPosition( + Double_t decSource, Double_t raSource, + Double_t decStar, Double_t raStar) +{ + fDecSource = decSource; + fRaSource = raSource; + + fDecStar = decStar; + fRaStar = raStar; + + *fLog << all << "MSourcePosfromStarPos::SetSourceAndStarPosition; fDecSource, fRaSource, fDecStar, fRaStar were set to : " + << fDecSource << ", " << fRaSource << ", " + << fDecStar << ", " << fRaStar << endl; +} + +// -------------------------------------------------------------------------- +// +// +Int_t MSourcePosfromStarPos::PreProcess(MParList *pList) +{ + MGeomCam *geom = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam")); + if (!geom) + { + *fLog << err << "MSourcePosfromStarPos : MGeomCam (Camera Geometry) missing in Parameter List... aborting." << endl; + return kFALSE; + } + fMm2Deg = geom->GetConvMm2Deg(); + + + fRun = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); + if (!fRun) + { + *fLog << err << "MRawRunHeader not found... aborting." << endl; + return kFALSE; + } + + + fSrcPos = (MSrcPosCam*)pList->FindObject(AddSerialNumber("MSrcPosCam")); + if (!fSrcPos) + { + *fLog << err << "MSourcePosfromStarPos : MSrcPosCam not found... aborting" << endl; + return kFALSE; + } + + //--------------------------------------------------------------------- + // read all files and call ReadData() to read and store the information + // + + fRuns = 0; + fSize = 0; + + while(1) + { + if (!OpenNextFile()) break; + + if (fIn->eof()) + if (!OpenNextFile()) break; + + // FIXME! Set InputStreamID + + ReadData(); + } + + FixSize(); + //------------------------------------------------------------- + + return kTRUE; +} + +//========================================================================= +// +// SourcefromStar +// +// this routine calculates the position of a source (for example Crab) in the camera +// from the position of a star (for example ZetaTauri) in the camera. The latter +// position may have been determined by analysing the DC currents in the different +// pixels. +// +// Input : thetaTel, phiTel the direction the telescope is pointing to, +// in local coordinates +// f the distance between camera and reflector +// decStar, raStar the position of the star in sky coordinates +// decSource, raSource the position of the source in sky coordinates +// xStar, yStar the position of the star in the camera +// dxStar, dyStar error of the position of the star in the camera +// +// Output : xSource, ySource the calculated position of the source in the camera +// dxSource, dySource error of the calculated position of the source in +// the camera +// +// Useful formulas can be found in TDAS 00-11 and TDAS 01-05 +// + +void MSourcePosfromStarPos::SourcefromStar(Double_t &f, + Double_t &decStar, Double_t &raStar, + Double_t &decSource, Double_t &raSource, + Double_t &thetaTel, Double_t &phiTel, + Double_t &xStar, Double_t &yStar, + Double_t &dxStar, Double_t &dyStar, + Double_t &xSource, Double_t &ySource, + Double_t &dxSource, Double_t &dySource) +{ + // the units are assumed to be radians for theta, phi, dec and ra + // and mm for f, x and y + + // calculate coordinates of star and source in system B (see TDAS 00-11, eqs. (2)) + Double_t xB = cos(decStar) * cos(raStar); + Double_t yB = cos(decStar) * sin(raStar); + Double_t zB = -sin(decStar); + + Double_t xB0 = cos(decSource) * cos(raSource); + Double_t yB0 = cos(decSource) * sin(raSource); + Double_t zB0 = -sin(decSource); + + // calculate rotation angle alpha of sky image in camera + // (see TDAS 00-11, eqs. (18) and (20)) + // a1 = cos(Lat), a3 = -sin(Lat), where Lat is the geographical latitude of La Palma + Double_t a1 = 0.876627; + Double_t a3 = -0.481171; + + Double_t denom = 1./ sqrt( sin(thetaTel)*sin(phiTel) * sin(thetaTel)*sin(phiTel) + + ( a1*cos(thetaTel)+a3*sin(thetaTel)*cos(phiTel) ) * + ( a1*cos(thetaTel)+a3*sin(thetaTel)*cos(phiTel) ) ); + Double_t cosal = - (a3 * sin(thetaTel) + a1 * cos(thetaTel) * cos(phiTel)) * denom; + Double_t sinal = a1 * sin(phiTel) * denom; + + + // calculate coordinates of star in a system with the basis vectors e1, e2, e3 + // where e1 is in the direction (r0 x a) + // e2 is in the direction (e1 x r0) + // and e3 is in the direction -r0; + // r0 is the direction the telescope is pointing to + // and a is the earth rotation axis (pointing to the celestial north pole) + // + Double_t x = (-xB*yB0 + xB0*yB) / sqrt( xB0*xB0 + yB0*yB0 ); + Double_t y = ( xB*xB0*zB0 + yB*yB0*zB0 - zB*(xB0*xB0 + yB0*yB0) ) + / sqrt( xB0*xB0 + yB0*yB0 ); + Double_t z = -(xB*xB0 + yB*yB0 + zB*zB0); + + // calculate coordinates of star in camera + Double_t xtilde = -f/z * (cosal*x - sinal*y); + Double_t ytilde = -f/z * (sinal*x + cosal*y); + + // calculate coordinates of source in camera + // note : in real camera signs are inverted (therefore s = -1.0) + Double_t s = -1.0; + xSource = s * (s*xStar - xtilde); + ySource = s * (s*yStar - ytilde); + + *fLog << "thetaTel, phiTel, cosal, sinal = " << thetaTel << ", " + << phiTel << ", " << cosal << ", " << sinal << endl; +} + + +// -------------------------------------------------------------------------- +// +// Get the source position and put it into MSrcPosCam +// +// +Bool_t MSourcePosfromStarPos::ReInit(MParList *pList) +{ + if (fDecStar == 0.0 || fRaStar == 0.0 || + fDecSource == 0.0 || fRaSource == 0.0) + { + *fLog << err << "MSourcePosfromStarPos::ReInit; sky coordinates of star and source are not defined ... aborting" + << endl; + return kFALSE; + } + + // f is the distance of the camera from the reflector center in [mm] + Double_t f = kRad2Deg / fMm2Deg; + + for (Int_t i=0; iGetRunNumber(); + if (run == fRunNr[i]) + { + MSourcePosfromStarPos::SourcefromStar( f, + fDecStar, fRaStar, fDecSource, fRaSource, + fThetaTel[i], fPhiTel[i], + fxStar[i], fyStar[i], + fdxStar[i], fdyStar[i], + fxSource, fySource, + fdxSource, fdySource); + + fSrcPos->SetXY(fxSource, fySource); + + *fLog << all << "MSourcePosfromStarPos::ReInit; f = " << f << endl; + *fLog << all << "MSourcePosfromStarPos::ReInit; fRunNr, fxSource, fySource = " + << fRunNr[i] << ", " << fxSource << ", " << fySource + << endl; + + fSrcPos->SetReadyToSave(); + } + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// +Int_t MSourcePosfromStarPos::Process() +{ + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// +Int_t MSourcePosfromStarPos::PostProcess() +{ + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// read the data from the ASCII file and store them +// +void MSourcePosfromStarPos::FixSize() +{ + fSize = fRuns; + + fRunNr.Set(fSize); + + fThetaTel.Set(fSize); + fPhiTel.Set(fSize); + fdThetaTel.Set(fSize); + fdPhiTel.Set(fSize); + + fxStar.Set(fSize); + fyStar.Set(fSize); + fdxStar.Set(fSize); + fdyStar.Set(fSize); +} + +// -------------------------------------------------------------------------- +// +// read the data from the ASCII file and store them +// +void MSourcePosfromStarPos::ReadData() +{ + Float_t val; + + if (fRuns >= fSize) + { + fSize += 100; + + fRunNr.Set(fSize); + + fThetaTel.Set(fSize); + fPhiTel.Set(fSize); + fdThetaTel.Set(fSize); + fdPhiTel.Set(fSize); + + fxStar.Set(fSize); + fyStar.Set(fSize); + fdxStar.Set(fSize); + fdyStar.Set(fSize); + } + + fRuns += 1; + + *fIn >> val; + fRunNr.AddAt(val, fRuns-1); + + *fIn >> val; + fThetaTel.AddAt(val, fRuns-1); + *fIn >> val; + fPhiTel.AddAt(val, fRuns-1); + + *fIn >> val; + fdThetaTel.AddAt(val, fRuns-1); + *fIn >> val; + fdPhiTel.AddAt(val, fRuns-1); + + *fIn >> val; + fxStar.AddAt(val, fRuns-1); + *fIn >> val; + fyStar.AddAt(val, fRuns-1); + + *fIn >> val; + fdxStar.AddAt(val, fRuns-1); + *fIn >> val; + fdyStar.AddAt(val, fRuns-1); +} + +// -------------------------------------------------------------------------- +// +// Add this file as the last entry in the chain +// +Int_t MSourcePosfromStarPos::AddFile(const char *txt, Int_t) +{ + TNamed *name = new TNamed(txt, ""); + fFileNames->AddLast(name); + return 1; +} + +// -------------------------------------------------------------------------- +// +// This opens the next file in the list and deletes its name from the list. +// +Bool_t MSourcePosfromStarPos::OpenNextFile() +{ + // + // open the input stream and check if it is really open (file exists?) + // + if (fIn) + delete fIn; + fIn = NULL; + + // + // Check for the existence of a next file to read + // + TNamed *file = (TNamed*)fFileNames->First(); + if (!file) + return kFALSE; + + // + // open the file which is the first one in the chain + // + const char *name = file->GetName(); + + const char *expname = gSystem->ExpandPathName(name); + fIn = new ifstream(expname); + delete [] expname; + + const Bool_t noexist = !(*fIn); + + if (noexist) + *fLog << dbginf << "Cannot open file '" << name << "'" << endl; + else + *fLog << "Open file: '" << name << "'" << endl; + + // + // Remove this file from the list of pending files + // + fFileNames->Remove(file); + + return !noexist; +} + +// -------------------------------------------------------------------------- + + + + + + + + + + Index: /trunk/MagicSoft/Mars/manalysis/MSourcePosfromStarPos.h =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MSourcePosfromStarPos.h (revision 3209) +++ /trunk/MagicSoft/Mars/manalysis/MSourcePosfromStarPos.h (revision 3209) @@ -0,0 +1,99 @@ +#ifndef MARS_MSourcePosfromStarPos +#define MARS_MSourcePosfromStarPos + +///////////////////////////////////////////////////////////////////////////// +// // +// MSourcePosfromStarPos // +// // +// Task to calculate the position of the source in the camera from +// the position of a known star in the camera +// // +///////////////////////////////////////////////////////////////////////////// + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#ifndef ROOT_TArrayI +#include +#endif + +#ifndef ROOT_TArrayD +#include +#endif + +class TList; +class MRawRunHeader; +class MGeomCam; +class MSrcPosCam; + + +class MSourcePosfromStarPos : public MTask +{ + private: + const MRawRunHeader *fRun; //! + const MGeomCam *fGeomCam; //! Camera Geometry used to calculate Hillas + MSrcPosCam *fSrcPos; //! + + ifstream *fIn; // input file + TList *fFileNames; // array which contains the \0-terminated file names + + Float_t fMm2Deg; + + Int_t fRuns; // current number of entries in TArray + Int_t fSize; // final number of entries in TArray + + Double_t fDecStar; + Double_t fRaStar; + + Double_t fDecSource; + Double_t fRaSource; + + Double_t fxSource; + Double_t fySource; + Double_t fdxSource; + Double_t fdySource; + + TArrayI fRunNr; + TArrayD fThetaTel; + TArrayD fPhiTel; + TArrayD fdThetaTel; + TArrayD fdPhiTel; + TArrayD fxStar; + TArrayD fyStar; + TArrayD fdxStar; + TArrayD fdyStar; + + Int_t AddFile(const char *fname, Int_t dummy=-1); + Bool_t OpenNextFile(); + void ReadData(); + void FixSize(); + + Int_t PreProcess(MParList *pList); + Bool_t ReInit(MParList *pList); + Int_t Process(); + Int_t PostProcess(); + +public: + MSourcePosfromStarPos(const char *name=NULL, const char *title=NULL); + ~MSourcePosfromStarPos(); + + void SetSourceAndStarPosition(Double_t decSource, Double_t raSource, + Double_t decStar, Double_t raStar); + + void SourcefromStar(Double_t &, Double_t &, Double_t &, + Double_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &, + Double_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t & ); + + ClassDef(MSourcePosfromStarPos, 0) // Task to calculate the source position from a star position +}; + +#endif + + + + + + + + Index: /trunk/MagicSoft/Mars/manalysis/Makefile =================================================================== --- /trunk/MagicSoft/Mars/manalysis/Makefile (revision 3208) +++ /trunk/MagicSoft/Mars/manalysis/Makefile (revision 3209) @@ -88,4 +88,5 @@ MArrivalTimeCalc.cc \ MArrivalTimeCalc2.cc \ + MSourcePosfromStarPos.cc \ MMcCalibrationUpdate.cc