Index: trunk/MagicSoft/Mars/mtemp/mifae/Changelog =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/Changelog (revision 3899) +++ trunk/MagicSoft/Mars/mtemp/mifae/Changelog (revision 3900) @@ -19,4 +19,20 @@ -*-*- END OF LINE -*-*- + 2004/04/30: Javier López + * alpha_plot.C + - macro to display hillas disttributions with corrected miss pointing + + * srcPosRun.sh + - Script to produce the file with the positions of UMa51 using psffit executable + + * MSrcPosFromFile.[h,cc] + - Tool to correct the position of the source for on and off data from a file + + * psffit.C & psffit.cc + - Macro and executable to calculate the mean position of a star. + + * MPSFFit.[h,cc] & MPSFFitCalc.[h,cc] + - Container and taks to do calculation of the psf and position of stars. + 2004/04/28: Javier Rico * makeHillas.cc, makehillas.datacard Index: trunk/MagicSoft/Mars/mtemp/mifae/MPSFFit.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/MPSFFit.cc (revision 3900) +++ trunk/MagicSoft/Mars/mtemp/mifae/MPSFFit.cc (revision 3900) @@ -0,0 +1,86 @@ +/* ======================================================================== *\ +! +! * +! * 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 López , 4/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ +#include "MPSFFit.h" + +#include + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MPSFFit); + +using namespace std; + +MPSFFit::MPSFFit(const char *name, const char *title) +{ + + fName = name ? name : "MPSFFit"; + fTitle = title ? title : "Container that holds the values of the PSF fit."; + + fMaximun = -666.; + fMeanMinorAxis = -666.; + fMeanMajorAxis = -666.; + fSigmaMinorAxis = -666.; + fSigmaMajorAxis = -666.; + fChisquare = -666.; +} + +void MPSFFit::Print(Option_t *opt) const +{ + + *fLog << all << GetDescriptor() << " " << GetTitle() << endl; + + *fLog << all << " Maximun \t " << fMaximun << " uA" << endl; + *fLog << all << " Mean Minor Axis \t " << fMeanMinorAxis << " mm \t"; + *fLog << all << " Sigma Minor Axis \t " << fSigmaMinorAxis << " mm" << endl; + *fLog << all << " Mean Major Axis \t " << fMeanMajorAxis << " mm \t"; + *fLog << all << " Sigma Major Axis \t " << fSigmaMajorAxis << " mm" << endl; + *fLog << all << " Chi Square \t " << fChisquare; + *fLog << all << endl; + +} + +// -------------------------------------------------------------------------- +// +// Paint the ellipse corresponding to the parameters +// +void MPSFFit::Paint(Option_t *opt) +{ + if (fSigmaMinorAxis<0 || fSigmaMajorAxis<0) + return; + + TEllipse e(fMeanMinorAxis, fMeanMajorAxis, fSigmaMinorAxis, fSigmaMajorAxis, 0, 360, 0); + e.SetLineWidth(2); + e.Paint(); +} + +void MPSFFit::Reset() +{ + SetMaximun(0.0); + SetMeanMinorAxis(0.0); + SetMeanMajorAxis(0.0); + SetSigmaMinorAxis(0.0); + SetSigmaMajorAxis(0.0); + SetChisquare(0.0); +} Index: trunk/MagicSoft/Mars/mtemp/mifae/MPSFFit.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/MPSFFit.h (revision 3900) +++ trunk/MagicSoft/Mars/mtemp/mifae/MPSFFit.h (revision 3900) @@ -0,0 +1,45 @@ +#ifndef MARS_MPSFFit +#define MARS_MPSFFit + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +class MPSFFit : public MParContainer +{ +private: + + Float_t fMaximun; + Float_t fMeanMinorAxis; + Float_t fMeanMajorAxis; + Float_t fSigmaMinorAxis; + Float_t fSigmaMajorAxis; + Float_t fChisquare; + +public: + MPSFFit(const char *name=NULL, const char *title=NULL); + //~MPSFFit(); + + void SetMaximun(Float_t Maximun_) {fMaximun=Maximun_;} + void SetMeanMinorAxis(Float_t MeanMinorAxis_) {fMeanMinorAxis=MeanMinorAxis_;} + void SetMeanMajorAxis(Float_t MeanMajorAxis_) {fMeanMajorAxis=MeanMajorAxis_;} + void SetSigmaMinorAxis(Float_t SigmaMinorAxis_) {fSigmaMinorAxis=SigmaMinorAxis_;} + void SetSigmaMajorAxis(Float_t SigmaMajorAxis_) {fSigmaMajorAxis=SigmaMajorAxis_;} + void SetChisquare(Float_t Chisquare_) {fChisquare=Chisquare_;} + + Float_t GetMaximun() {return fMaximun;} + Float_t GetMeanMinorAxis() {return fMeanMinorAxis;} + Float_t GetMeanMajorAxis() {return fMeanMajorAxis;} + Float_t GetSigmaMinorAxis() {return fSigmaMinorAxis;} + Float_t GetSigmaMajorAxis() {return fSigmaMajorAxis;} + Float_t GetChisquare() {return fChisquare;} + + void Reset(); + + void Paint(Option_t *opt=NULL); + void Print(Option_t *opt=NULL) const; + + ClassDef(MPSFFit, 1) // Container that holds the values of the PSF fit. +}; + +#endif Index: trunk/MagicSoft/Mars/mtemp/mifae/MPSFFitCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/MPSFFitCalc.cc (revision 3900) +++ trunk/MagicSoft/Mars/mtemp/mifae/MPSFFitCalc.cc (revision 3900) @@ -0,0 +1,569 @@ +/* ======================================================================== *\ +! +! * +! * 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 López , 4/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ +#include "MPSFFitCalc.h" + +#include + + +#include +#include +#include +#include +#include +#include +#include + +#include "MGeomPix.h" +#include "MGeomCam.h" +#include "MMinuitInterface.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" + +ClassImp(MPSFFitCalc); + +using namespace std; + + +const Int_t numVar = 5; +const Float_t pixelSize = 31.5; //[mm] +const Float_t sqrt2 = sqrt(2.); +const Float_t sqrt3 = sqrt(3.); +const Bool_t usePrintOut = kFALSE; +const Int_t minuitPrintOut = -1; + +UInt_t numPixels; +Bool_t isPixelUsed[577]; +UInt_t numPixelsUsed; +Float_t pixelPosX[577]; +Float_t pixelPosY[577]; +Float_t pixelValue[577]; +Float_t ChiSquare; + +MPSFFitCalc::MPSFFitCalc(ImgCleanMode_t imgmode, const char *name, const char *title) +{ + fName = name ? name : "MPSFFitCalc"; + fTitle = title ? title : "Task that fits the PSF using the dc readout of the camera."; + + fImgCleanMode = imgmode; + fNumRings = 6; + fPedLevel = 3.0; + + +// Initialization of the camera dc mask 'fIsPixelused[]' + numPixels = 577; + for (UInt_t pixid=0; pixid0.) + { + x = pixelPosX[pixid]; + y = pixelPosY[pixid]; + z = pixelValue[pixid]; + + if (errorz > 0.0) + { + delta = (z-func(x,y,par))/errorz; + chisq += delta*delta; + } + else + cerr << "ERROR fcnPSF: This should never happen errorz[" << pixid << "] " << errorz << endl; + } + } + f = chisq; + ChiSquare = chisq; + +} + + +void MPSFFitCalc::InitFitVariables() +{ + + for (UInt_t pixid=1; pixid<577; pixid++) + pixelValue[pixid]=(Float_t)((*fCamera)[pixid]); + + numPixelsUsed = 0; + Float_t totalDC = 0.0; + + fVinit[0] = fMaxDC; + + if(usePrintOut) + *fLog << dbg << "Pixels used in the fit \t"; + + for (UInt_t idx=0; idxFill(pixelValue[pixidx]); + + fPedestalDC = fPedestalDCHist->GetBinCenter(fPedestalDCHist->GetMaximumBin()); + + for (UInt_t pixid=1; pixid<577; pixid++) + pixelValue[pixid]-=fPedestalDC; + + if(usePrintOut) + { + *fLog << dbg << "numPixelsUsed \t" << numPixelsUsed << endl; + *fLog << dbg << "fPedestalDC \t" << fPedestalDC << endl; + *fLog << dbg << "Maximun DC Init. value \t" << fVinit[0] << endl; + *fLog << dbg << "Mean Minor Axis Init. value \t" << fVinit[1] << " mm\t"; + *fLog << dbg << "Sigma Minor Axis Init. value \t" << fVinit[2] << endl; + *fLog << dbg << "Mean Major Axis Init. value \t" << fVinit[3] << " mm\t"; + *fLog << dbg << "Sigma Major Axis Init. value \t" << fVinit[4] << endl; + } +} + +void MPSFFitCalc::RingImgClean() +{ + + Bool_t isPixelUsedTmp[577]; + + fMaxDC=0; + UInt_t maxDCpix[2]; + + +// Look for the two neighbor pixels with the maximun signal and set all pixels as unused + Float_t dc[2]; + Float_t dcsum; + +// Find the two close pixels with the maximun dc + for(UInt_t pixidx=0; pixidx fMaxDC*2) + { + maxDCpix[0] = pixidx; + maxDCpix[1] = swneighbor; + fMaxDC = dcsum/2; + } + } + } + } + } + +// Those variables are: +// 1. an array of 2 dimensions +// 1.1 the first dimension store the ring around the 'maxiun signal pixel' +// 1.2 the sw numbers of the pixels in this ring +// 2. an array with the number of pixels in each ring + UInt_t isPixelUesdInRing[fNumRings][577]; + UInt_t numPixelsUsedInRing[fNumRings]; + +// Store the 'maximun signal pixel in the [0] ring and set it as used + + for (Int_t core=0; core<2; core++) + { + isPixelUesdInRing[0][0] = maxDCpix[core]; + numPixelsUsedInRing[0] = 1; + isPixelUsedTmp[isPixelUesdInRing[0][0]] = kTRUE; + + UInt_t count; + +// Loop over the neighbors of the previus ring and store the sw numbers in the 2D array to be +// use in the next loop/ring + for (UShort_t ring=0; ringFindObject(AddSerialNumber("MGeomCam")); + + if (!fGeomCam) + { + *fLog << err << AddSerialNumber("MGeomCam") << " not found ... aborting" << endl; + return kFALSE; + } + +// Initialize positions of the pixels to be used in the minuit minimizations + + for (UInt_t pixid=1; pixid<577; pixid++) + { + MGeomPix &pix=(*fGeomCam)[pixid]; + pixelPosX[pixid] = pix.GetX(); + pixelPosY[pixid] = pix.GetY(); + } + + fCamera = (MCameraDC*)pList->FindObject(AddSerialNumber("MCameraDC")); + + if (!fCamera) + { + *fLog << err << AddSerialNumber("MCameraDC") << " not found ... aborting" << endl; + return kFALSE; + } + + fFit = (MPSFFit*)pList->FindCreateObj(AddSerialNumber("MPSFFit")); + if (!fFit) + { + *fLog << err << AddSerialNumber("MPSFFit") << " cannot be created ... aborting" << endl; + return kFALSE; + } + + // Minuit initialization + + TMinuit *gMinuit = new TMinuit(6); //initialize TMinuit with a maximum of 5 params + gMinuit->SetFCN(fcnPSF); + + Double_t arglist[10]; + Int_t ierflg = 0; + + arglist[0] = 1; + gMinuit->mnexcm("SET ERR", arglist ,1,ierflg); + arglist[0] = minuitPrintOut; + gMinuit->mnexcm("SET PRI", arglist ,1,ierflg); + + return kTRUE; + +} + +Int_t MPSFFitCalc::Process() +{ + + // ------------------------------------------ + // image cleaning + + switch(fImgCleanMode) + { + case kRing: + { + RingImgClean(); + break; + } + case kRms: + { + RmsImgClean(fPedLevel); + break; + } + case kMask: + { + MaskImgClean(fBlindPixels); + break; + } + case kCombined: + { + MaskImgClean(fBlindPixels); + RingImgClean(); + MaskImgClean(fBlindPixels); + break; + } + default: + { + *fLog << err << "Image Cleaning mode " << fImgCleanMode << " not defined" << endl; + return kFALSE; + } + } + + InitFitVariables(); + + // ------------------------------------------- + // 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); + + clock.Stop(); + + if(usePrintOut) + { + *fLog << dbg << "Time spent for the minimization in MINUIT : " << endl;; + clock.Print(); + } + + if (ierflg) + { + *fLog << err << "MMinuitInterface::CallMinuit error " << ierflg << endl; + return kCONTINUE; + } + + fNumTotalMeanFits++; + + Double_t parm,parmerr; + + gMinuit->GetParameter(0,parm,parmerr); + fFit->SetMaximun(parm); + fTotalMeanFit.SetMaximun(fTotalMeanFit.GetMaximun()+parm); + + gMinuit->GetParameter(1,parm,parmerr); + fFit->SetMeanMinorAxis(parm); + fTotalMeanFit.SetMeanMinorAxis(fTotalMeanFit.GetMeanMinorAxis()+parm); + + gMinuit->GetParameter(2,parm,parmerr); + fFit->SetSigmaMinorAxis(parm); + fTotalMeanFit.SetSigmaMinorAxis(fTotalMeanFit.GetSigmaMinorAxis()+parm); + + gMinuit->GetParameter(3,parm,parmerr); + fFit->SetMeanMajorAxis(parm); + fTotalMeanFit.SetMeanMajorAxis(fTotalMeanFit.GetMeanMajorAxis()+parm); + + gMinuit->GetParameter(4,parm,parmerr); + fFit->SetSigmaMajorAxis(parm); + fTotalMeanFit.SetSigmaMajorAxis(fTotalMeanFit.GetSigmaMajorAxis()+parm); + + fFit->SetChisquare(ChiSquare/(numPixelsUsed-numVar)); + fTotalMeanFit.SetChisquare(fTotalMeanFit.GetChisquare()+ChiSquare/(numPixelsUsed-numVar)); + + if(usePrintOut) + { + fFit->Print(); + fTotalMeanFit.Print(); + } + + return kTRUE; +} + +Int_t MPSFFitCalc::PostProcess() +{ + + fTotalMeanFit.SetMaximun(fTotalMeanFit.GetMaximun()/fNumTotalMeanFits); + fTotalMeanFit.SetMeanMinorAxis(fTotalMeanFit.GetMeanMinorAxis()/fNumTotalMeanFits); + fTotalMeanFit.SetSigmaMinorAxis(fTotalMeanFit.GetSigmaMinorAxis()/fNumTotalMeanFits); + fTotalMeanFit.SetMeanMajorAxis(fTotalMeanFit.GetMeanMajorAxis()/fNumTotalMeanFits); + fTotalMeanFit.SetSigmaMajorAxis(fTotalMeanFit.GetSigmaMajorAxis()/fNumTotalMeanFits); + fTotalMeanFit.SetChisquare(fTotalMeanFit.GetChisquare()/fNumTotalMeanFits); + + fFit->SetMaximun(fTotalMeanFit.GetMaximun()); + fFit->SetMeanMinorAxis(fTotalMeanFit.GetMeanMinorAxis()); + fFit->SetSigmaMinorAxis(fTotalMeanFit.GetSigmaMinorAxis()); + fFit->SetMeanMajorAxis(fTotalMeanFit.GetMeanMajorAxis()); + fFit->SetSigmaMajorAxis(fTotalMeanFit.GetSigmaMajorAxis()); + fFit->SetChisquare(fTotalMeanFit.GetChisquare()); + + if(usePrintOut) + fTotalMeanFit.Print(); + + return kTRUE; +} Index: trunk/MagicSoft/Mars/mtemp/mifae/MPSFFitCalc.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/MPSFFitCalc.h (revision 3900) +++ trunk/MagicSoft/Mars/mtemp/mifae/MPSFFitCalc.h (revision 3900) @@ -0,0 +1,93 @@ +///////////////////////////////////////////////////////////////////////////// +// // +// MPSFFitCalc // +// // +// // +///////////////////////////////////////////////////////////////////////////// +#ifndef MARS_MPSFFitCalc +#define MARS_MPSFFitCalc + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#ifndef ROOT_TArrayD +#include +#endif +#ifndef ROOT_TArrayI +#include +#endif +#ifndef ROOT_TArrayS +#include +#endif + +class TSring; +class TH1D; + +class MGeomCam; +#ifndef MARS_MCameraDC +#include "MCameraDC.h" +#endif +#ifndef MARS_MPSFFit +#include "MPSFFit.h" +#endif + +class MPSFFitCalc : public MTask +{ + +private: + + MGeomCam *fGeomCam; + MCameraDC *fCamera; + MPSFFit *fFit; + + MPSFFit fTotalMeanFit; + UInt_t fNumTotalMeanFits; + + TH1D *fPedestalDCHist; + Float_t fPedestalDC; + + TString *fVname; + TArrayD fVinit; + TArrayD fStep; + TArrayD fLimlo; + TArrayD fLimup; + TArrayI fFix; + TObject *fObjectFit; + TString fMethod; + Bool_t fNulloutput; + + Float_t fMaxDC; + UShort_t fImgCleanMode; + UShort_t fNumRings; + Float_t fPedLevel; + TArrayS fBlindPixels; + + void RingImgClean(); + void RmsImgClean(Float_t pedestal); + void MaskImgClean(TArrayS blindpixels); + + void InitFitVariables(); + + Int_t PreProcess(MParList *pList); + Int_t Process(); + Int_t PostProcess(); + +public: + + enum ImgCleanMode_t{kNone=0,kRing,kRms,kMask,kCombined}; + static const ImgCleanMode_t kDefault = kRing; + + MPSFFitCalc(ImgCleanMode_t imgmode=kDefault, const char *name=NULL, const char *title=NULL); + ~MPSFFitCalc(); + + void SetImgCleanMode(ImgCleanMode_t imgmode=kDefault){ fImgCleanMode=imgmode;}; + void SetNumRings(UShort_t numRings) {fNumRings=numRings;} + void SetPedLevel(Float_t pedestal) {fPedLevel=pedestal;} + void SetBlindPixels(TArrayS blindpixels) {fBlindPixels=blindpixels;} + + + ClassDef(MPSFFitCalc, 0) // Task that fits the PSF using the dc readout of the camera. +}; + +#endif Index: trunk/MagicSoft/Mars/mtemp/mifae/MSrcPosFromFile.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/MSrcPosFromFile.cc (revision 3900) +++ trunk/MagicSoft/Mars/mtemp/mifae/MSrcPosFromFile.cc (revision 3900) @@ -0,0 +1,248 @@ +/* ======================================================================== *\ +! +! * +! * 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 López 04/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MSrcPosFromFile +// +// Task to calculate the position of the source as a function of run number +// +// Output Containers: +// MSrcPosCam +// +////////////////////////////////////////////////////////////////////////////// + +#include + +#include + +#include "MParList.h" + +#include "MSrcPosFromFile.h" + +#include "MRawRunHeader.h" +#include "MSrcPosCam.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MSrcPosFromFile); + +using namespace std; + +static const TString gsDefName = "MSrcPosFromFile"; +static const TString gsDefTitle = "Calculate position of the (off axis) source"; + +// ------------------------------------------------------------------------- +// +// Default constructor +// +MSrcPosFromFile::MSrcPosFromFile(TString cardpath, OnOffMode_t mode, const char *name, const char *title) + : fRawRunHeader(NULL), fSrcPos(NULL), fMode(mode), fSourcePositionFilePath(cardpath) +{ + fName = name ? name : gsDefName.Data(); + fTitle = title ? title : gsDefTitle.Data(); + + fLastRun = 0; + +// Count the number of runs in the card with the source poistions + ReadSourcePositionsFile(kCount); + + fRunList = new Int_t[fNumRuns]; + fRunSrcPos = new MSrcPosCam[fNumRuns]; + fRunMap = new TExMap(fNumRuns); + + Float_t cameraSize = 600; //[mm] + Float_t binPrecision = 3; //[mm] ~ 0.01 deg + + UInt_t nbins = (UInt_t)(cameraSize*2/binPrecision); + + fHistSrcPos = new TH2F("HistSrcPos","",nbins,-cameraSize,cameraSize,nbins,-cameraSize,cameraSize); + +// Read card with the source poistions + ReadSourcePositionsFile(kRead); +} + +MSrcPosFromFile::~MSrcPosFromFile() +{ + delete [] fRunList; + delete [] fRunSrcPos; + delete fRunMap; + delete fHistSrcPos; +} + +// ------------------------------------------------------------------------- +// +Int_t MSrcPosFromFile::PreProcess(MParList *pList) +{ + + fRawRunHeader = (MRawRunHeader*)pList->FindObject(AddSerialNumber("MRawRunHeader")); + if (!fRawRunHeader) + { + *fLog << err << AddSerialNumber("MRawRunHeader") << " not found ... aborting" << endl; + return kFALSE; + } + + fSrcPos = (MSrcPosCam*)pList->FindCreateObj("MSrcPosCam"); + if (!fSrcPos) + return kFALSE; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// +Bool_t MSrcPosFromFile::ReInit(MParList *pList) +{ + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// +Int_t MSrcPosFromFile::Process() +{ + + switch(fMode) + { + case kOn: + { + const UInt_t run = fRawRunHeader->GetRunNumber(); + + MSrcPosCam* srcpos = (MSrcPosCam*)fRunMap->GetValue(run); + + Float_t x; + Float_t y; + + if (srcpos) + { + x = srcpos->GetX(); + y = srcpos->GetY(); + + + if (srcpos && run != fLastRun) + { + fSrcPos->SetXY(x,y); + + *fLog << inf << "Source position for run " << run; + *fLog << inf << "\tX\t" << setprecision(2) << x; + *fLog << inf << "\tY\t" << setprecision(2) << y << endl; + } + + fLastRun = run; + } + + x = fSrcPos->GetX(); + y = fSrcPos->GetY(); + + fHistSrcPos->Fill(x,y); + break; + } + case kOff: + { + Axis_t x; + Axis_t y; + + fHistSrcPos->GetRandom2(x,y); + fSrcPos->SetXY(x,y); + + break; + } + default: + *fLog << err << "Wrond mode " << fMode << endl; + return kFALSE; + } + + return kTRUE; +} + + +Int_t MSrcPosFromFile::PostProcess() +{ + + *fLog << dbg << endl; + *fLog << dbg << "fHistSrcPos->GetEntries() " << fHistSrcPos->GetEntries() << endl; + *fLog << dbg << "fHistSrcPos->ProjectionX()->GetMean() " << fHistSrcPos->ProjectionX()->GetMean() << endl; + *fLog << dbg << "fHistSrcPos->ProjectionX()->GetRMS() " << fHistSrcPos->ProjectionX()->GetRMS() << endl; + *fLog << dbg << "fHistSrcPos->ProjectionY()->GetMean() " << fHistSrcPos->ProjectionY()->GetMean() << endl; + *fLog << dbg << "fHistSrcPos->ProjectionY()->GetRMS() " << fHistSrcPos->ProjectionY()->GetRMS() << endl; + + return kTRUE; +} + + +Int_t MSrcPosFromFile::ReadSourcePositionsFile(UShort_t readmode) +{ + + ifstream fin(fSourcePositionFilePath); + if(!fin) + { + *fLog << err << "MSrcPosFromFile::ReadSourcePositionsFile. Cannot open file " << fSourcePositionFilePath << endl; + return kFALSE; + } + + UInt_t run; + Float_t x,y; + + fNumRuns=0; + + *fLog << dbg << "MSrcPosFromFile::ReadSourcePositionsFile(" << readmode << ')' << endl; + while(1) + { + fin >> run >> x >> y; + if(fin.eof()) + break; + + switch(readmode) + { + case kCount: + { + + *fLog << dbg << "Source position for run " << run; + *fLog << dbg << "\tX\t" << x << " mm"; + *fLog << dbg << "\tY\t" << y << " mm" << endl; + + fNumRuns++; + break; + } + case kRead: + { + fRunList[fNumRuns] = run; + fRunSrcPos[fNumRuns].SetXY(x,y); + fRunMap->Add(fRunList[fNumRuns],(Long_t)&(fRunSrcPos[fNumRuns])); + fNumRuns++; + break; + } + default: + *fLog << err << "Read mode " << readmode << " node defined" << endl; + return kFALSE; + } + } + + fin.close(); + + + return kTRUE; +} Index: trunk/MagicSoft/Mars/mtemp/mifae/MSrcPosFromFile.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/MSrcPosFromFile.h (revision 3900) +++ trunk/MagicSoft/Mars/mtemp/mifae/MSrcPosFromFile.h (revision 3900) @@ -0,0 +1,55 @@ +#ifndef MARS_MSrcPosFromFile +#define MARS_MSrcPosFromFile + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#ifndef ROOT_TMap +#include +#endif + +class TH2F; + +class MRawRunHeader; +class MSrcPosCam; + +class MSrcPosFromFile : public MTask +{ + private: + + MRawRunHeader *fRawRunHeader; + MSrcPosCam *fSrcPos; // Pointer to the source position + + Int_t fNumRuns; + UInt_t fLastRun; + + Int_t *fRunList; + MSrcPosCam *fRunSrcPos; + TExMap *fRunMap; // list of run numbers positions + + TH2F *fHistSrcPos; + UShort_t fMode; + + TString fSourcePositionFilePath; + Int_t ReadSourcePositionsFile(UShort_t readmode); + + Int_t PreProcess(MParList *plist); + Bool_t ReInit(MParList *plist); + Int_t Process(); + Int_t PostProcess(); + +public: + + enum ReadMode_t {kCount=0,kRead}; + enum OnOffMode_t {kOn=0,kOff}; + + MSrcPosFromFile(TString cardpath=NULL, OnOffMode_t mode=kOn, const char *name=NULL, const char *title=NULL); + ~MSrcPosFromFile(); + + void SetMode(OnOffMode_t mode) {fMode=mode;} + + ClassDef(MSrcPosFromFile, 0) // task to calculate the position of the source as a function of the run number +}; + +#endif Index: trunk/MagicSoft/Mars/mtemp/mifae/alpha_plot.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/alpha_plot.C (revision 3900) +++ trunk/MagicSoft/Mars/mtemp/mifae/alpha_plot.C (revision 3900) @@ -0,0 +1,834 @@ +Double_t ChiSquareNDof(TH1D *h1, TH1D *h2) +{ + Double_t chiq = 0.; + Double_t chi; + Double_t error; + Int_t nbinsnozero = 0; + + Int_t nbins = h1->GetNbinsX(); + if (nbins != h2->GetNbinsX() || nbins == 0) + return -1; + + for (UInt_t bin=1; bin<=nbins; bin++) + { + error = sqrt(h1->GetBinError(bin)*h1->GetBinError(bin) + + h2->GetBinError(bin)*h2->GetBinError(bin)); + if (error != 0) + { + chi = (h1->GetBinContent(bin)-h2->GetBinContent(bin))/error; + chiq += chi*chi; + nbinsnozero++; + } + } + + return chiq/nbinsnozero; +} + +//void alpha_plot(TString f_on_name = "../HillasFiles/20040319_Mrk421_30_15.KDummy*.root", +void alpha_plot(TString f_on_name = "../HillasFiles/20040319_OffMrk421_30_15.KDummy*.root", + TString f_off_name = "../HillasFiles/20040319_OffMrk421_30_15.KDummy*.root", + TString f_src_name = "../HillasFiles/20040319_OffMrk421.fake.pos") +/*void alpha_plot(TString f_on_name = "../HillasFiles/mrk421OnMisp3015*.root", + TString f_off_name = "../HillasFiles/mrk421OffMisp3015*.root", + TString f_src_name = "../HillasFiles/20040215_Mrk421.pos") +*/ +{ + + const UInt_t numEntries = 100000000; + + //cuts + Float_t sizemin = 500.; //[ADC] + Float_t sizemax = 10000000000.; //[ADC] + Float_t widthmin = 0.; + Float_t widthmax = 0.8; + Float_t lengthmin = 0.; + Float_t lengthmax = 0.8; + Float_t distmin = 0.; + Float_t distmax = 2.; + Float_t alphamin = 0.; + Float_t alphamax = 90.; + + //integration + Float_t sigexccmin = 0.; + Float_t sigexccmax = 30.; + Float_t bkgnormmin = 0.; + Float_t bkgnormmax = 90.; + + gStyle->SetOptStat(111111); + gStyle->SetOptFit(); + + // + // Make a loop only for the ON data: + // + + MParList plist_on; + MTaskList tlist_on; + plist_on.AddToList(&tlist_on); + + // ON containers + MGeomCamMagic geomcam; + MSrcPosCam source_on; + MHillas hillas; + MHillasSrc hillasscr; + + plist_on.AddToList(&geomcam); + plist_on.AddToList(&source_on); + plist_on.AddToList(&hillas); + plist_on.AddToList(&hillasscr); + + //create some 1-dim histo to test only for the ON distribution of dist, width , length, size... + MH3 hDist_on("MHillasSrc.fDist/315."); + hDist_on.SetName("Dist_on"); + plist_on.AddToList(&hDist_on); + MBinning binsDist_on("BinningDist_on"); + Int_t nbins_Dist = 20; + Float_t min_Dist = 0.; + Float_t max_Dist = distmax*1.2; + binsDist_on.SetEdges(nbins_Dist, min_Dist, max_Dist); + plist_on.AddToList(&binsDist_on); + + MH3 hWidth_on("MHillas.fWidth/315."); + hWidth_on.SetName("Width_on"); + plist_on.AddToList(&hWidth_on); + MBinning binsWidth_on("BinningWidth_on"); + Int_t nbins_Width = 20; + Float_t min_Width = 0.; + Float_t max_Width = widthmax*1.2; + binsWidth_on.SetEdges(nbins_Width, min_Width, max_Width); + plist_on.AddToList(&binsWidth_on); + + MH3 hLength_on("MHillas.fLength/315."); + hLength_on.SetName("Length_on"); + plist_on.AddToList(&hLength_on); + MBinning binsLength_on("BinningLength_on"); + Int_t nbins_Length = 20; + Float_t min_Length = 0.; + Float_t max_Length = lengthmax*1.2; + binsLength_on.SetEdges(nbins_Length, min_Length, max_Length); + plist_on.AddToList(&binsLength_on); + + MH3 hSize_on("log10(MHillas.fSize)"); + hSize_on.SetName("Size_on"); + plist_on.AddToList(&hSize_on); + MBinning binsSize_on("BinningSize_on"); + Int_t nbins_Size = 60; + Float_t min_Size = log10(sizemin)*0.8; + Float_t max_Size = log10(1000000)*1.2; + binsSize_on.SetEdges(nbins_Size, min_Size, max_Size); + plist_on.AddToList(&binsSize_on); + + //create a histo to fill the alpha values: one alpha plot form 0 to +90 deg in abs value + MH3 hAlpha_on_abs("abs(MHillasSrc.fAlpha)"); + hAlpha_on_abs.SetName("Alpha_on_abs"); + plist_on.AddToList(&hAlpha_on_abs); + MBinning binsAlpha_on_abs("BinningAlpha_on_abs"); + Int_t nbins_abs = 9; + Float_t minalpha_abs = 0.; + Float_t maxalpha_abs =90.; + binsAlpha_on_abs.SetEdges(nbins_abs, minalpha_abs, maxalpha_abs); + plist_on.AddToList(&binsAlpha_on_abs); + + //create a histo to fill the alpha values: one alpha plot form -90 to +90 deg. + MH3 hAlpha_on("MHillasSrc.fAlpha"); + hAlpha_on.SetName("Alpha_on"); + plist_on.AddToList(&hAlpha_on); + MBinning binsAlpha_on("BinningAlpha_on"); + Int_t nbins = nbins_abs*2; + Float_t minalpha = -90.; + Float_t maxalpha = 90.; + binsAlpha_on.SetEdges(nbins, minalpha, maxalpha); + plist_on.AddToList(&binsAlpha_on); + + + MH3 hSrcPos_on("MSrcPosCam.fX","MSrcPosCam.fY"); + hSrcPos_on.SetName("SrcPos_on"); + plist_on.AddToList(&hSrcPos_on); + MBinning binsSrcPos_onX("BinningSrcPos_onX"); + MBinning binsSrcPos_onY("BinningSrcPos_onY"); + Int_t nbins_srcpos = 400; + Float_t minsrcpos = -600.; + Float_t maxsrcpos = 600.; + binsSrcPos_onX.SetEdges(nbins_srcpos, minsrcpos, maxsrcpos); + binsSrcPos_onY.SetEdges(nbins_srcpos, minsrcpos, maxsrcpos); + plist_on.AddToList(&binsSrcPos_onX); + plist_on.AddToList(&binsSrcPos_onY); + + MH3 hDAQEvtNumber_on("MRawEvtHeader.fDAQEvtNumber"); + hDAQEvtNumber_on.SetName("DAQEvtNumber_on"); + plist_on.AddToList(&hDAQEvtNumber_on); + MBinning binsDAQEvtNumber_onX("BinningDAQEvtNumber_onX"); + Int_t nbins_evtnum = 1000; + Float_t minevtnum = 0.; + Float_t maxevtnum = 1000.; + binsDAQEvtNumber_onX.SetEdges(nbins_evtnum,minevtnum,maxevtnum); + plist_on.AddToList(&binsDAQEvtNumber_onX); + + // + //tasks + // + + MReadTree read_on("Parameters", f_on_name); + read_on.DisableAutoScheme(); + + + //cuts + TString sizestr = "MHillas.fSize < "; + sizestr += sizemin; + sizestr += " || "; + sizestr += "MHillas.fSize > "; + sizestr += sizemax; + MF sizefilter(sizestr); + + TString widthstr = "{MHillas.fWidth/315.} < "; + widthstr += widthmin; + widthstr += " || "; + widthstr += "{MHillas.fWidth/315.} > "; + widthstr += widthmax; + MF widthfilter(widthstr); + + TString lengthstr = "{MHillas.fLength/315.} < "; + lengthstr += lengthmin; + lengthstr += " || "; + lengthstr += "{MHillas.fLength/315.} > "; + lengthstr += lengthmax; + MF lengthfilter(lengthstr); + + TString diststr = "{MHillasSrc.fDist/315.} < "; + diststr += distmin; + diststr += " || "; + diststr += "{MHillasSrc.fDist/315.} > "; + diststr += distmax; + MF distfilter(diststr); + + TString alphastr = "{abs(MHillasSrc.fAlpha)} < "; + alphastr += alphamin; + alphastr += " || "; + alphastr += "{abs(MHillasSrc.fAlpha)} > "; + alphastr += alphamax; + MF alphafilter(alphastr); + + MF evenfilter("{MRawEvtHeader.fDAQEvtNumber%3}<0.5"); + MF oddfilter("{MRawEvtHeader.fDAQEvtNumber%3}>0.5"); + + MContinue cont_size(&sizefilter); + MContinue cont_width(&widthfilter); + MContinue cont_length(&lengthfilter); + MContinue cont_dist(&distfilter); + MContinue cont_alpha(&alphafilter); + MContinue cont_even(&evenfilter); + MContinue cont_odd(&oddfilter); + + MSrcPosFromFile srccalc(f_src_name); + + MHillasSrcCalc csrc_on; + + // fill all histograms + MFillH falpha_on_abs(&hAlpha_on_abs); + MFillH falpha_on(&hAlpha_on); + MFillH fdist_on(&hDist_on); + MFillH fwidth_on(&hWidth_on); + MFillH flength_on(&hLength_on); + MFillH fsize_on(&hSize_on); + MFillH fsrcpos_on(&hSrcPos_on); + MFillH fevtnum_on(&hDAQEvtNumber_on); + + + // prints + MPrint pevent("MRawEvtHeader"); + MPrint phillas("MHillas"); + MPrint phillassrc("MHillasSrc"); + MPrint psrcpos("MSrcPosCam"); + + //tasklist + tlist_on.AddToList(&read_on); + tlist_on.AddToList(&srccalc); + tlist_on.AddToList(&csrc_on); + tlist_on.AddToList(&fsrcpos_on); + tlist_on.AddToList(&cont_odd); + tlist_on.AddToList(&cont_size); + tlist_on.AddToList(&cont_width); + tlist_on.AddToList(&cont_length); + tlist_on.AddToList(&cont_dist); + tlist_on.AddToList(&cont_alpha); + tlist_on.AddToList(&falpha_on_abs); + tlist_on.AddToList(&falpha_on); + tlist_on.AddToList(&fdist_on); + tlist_on.AddToList(&fwidth_on); + tlist_on.AddToList(&flength_on); + tlist_on.AddToList(&fsize_on); + tlist_on.AddToList(&fevtnum_on); + + // Create and setup the eventloop + MEvtLoop loop_on; + loop_on.SetParList(&plist_on); + //loop_on.SetDisplay(display); + + MProgressBar bar; + loop_on.SetProgressBar(&bar); + + if (!loop_on.Eventloop(numEntries)) + return; + + tlist_on.PrintStatistics(); + + // + // Make a loop only for the OFF data: + // + + MParList plist_off; + MTaskList tlist_off; + plist_off.AddToList(&tlist_off); + + MSrcPosCam source_off; + + plist_off.AddToList(&geomcam); + plist_off.AddToList(&source_off); + plist_off.AddToList(&hillas); + plist_off.AddToList(&hillasscr); + + //create some 1-dim histo to test only for the OFF distribution of dist, width , length, size... + MH3 hDist_off("MHillasSrc.fDist/315."); + hDist_off.SetName("Dist_off"); + plist_off.AddToList(&hDist_off); + MBinning binsDist_off("BinningDist_off"); + binsDist_off.SetEdges(nbins_Dist, min_Dist, max_Dist); + plist_off.AddToList(&binsDist_off); + + MH3 hWidth_off("MHillas.fWidth/315."); + hWidth_off.SetName("Width_off"); + plist_off.AddToList(&hWidth_off); + MBinning binsWidth_off("BinningWidth_off"); + binsWidth_off.SetEdges(nbins_Width, min_Width, max_Width); + plist_off.AddToList(&binsWidth_off); + + MH3 hLength_off("MHillas.fLength/315."); + hLength_off.SetName("Length_off"); + plist_off.AddToList(&hLength_off); + MBinning binsLength_off("BinningLength_off"); + binsLength_off.SetEdges(nbins_Length, min_Length, max_Length); + plist_off.AddToList(&binsLength_off); + + MH3 hSize_off("log10(MHillas.fSize)"); + hSize_off.SetName("Size_off"); + plist_off.AddToList(&hSize_off); + MBinning binsSize_off("BinningSize_off"); + binsSize_off.SetEdges(nbins_Size, min_Size, max_Size); + plist_off.AddToList(&binsSize_off); + + //create a histo to fill the alpha values: from 0 to 90 deg -> abs value + MH3 hAlpha_off_abs("abs(MHillasSrc.fAlpha)"); + hAlpha_off_abs.SetName("Alpha_off_abs"); + plist_off.AddToList(&hAlpha_off_abs); + MBinning binsAlpha_off_abs("BinningAlpha_off_abs"); + binsAlpha_off_abs.SetEdges(nbins_abs, minalpha_abs, maxalpha_abs); + plist_off.AddToList(&binsAlpha_off_abs); + + //create a histo to fill the alpha values: from -90 to 90 deg + MH3 hAlpha_off("MHillasSrc.fAlpha"); + hAlpha_off.SetName("Alpha_off"); + plist_off.AddToList(&hAlpha_off); + MBinning binsAlpha_off("BinningAlpha_off"); + binsAlpha_off.SetEdges(nbins, minalpha, maxalpha); + plist_off.AddToList(&binsAlpha_off); + + + MH3 hSrcPos_off("MSrcPosCam.fX","MSrcPosCam.fY"); + hSrcPos_off.SetName("SrcPos_off"); + plist_off.AddToList(&hSrcPos_off); + MBinning binsSrcPos_offX("BinningSrcPos_offX"); + MBinning binsSrcPos_offY("BinningSrcPos_offY"); + binsSrcPos_offX.SetEdges(nbins_srcpos, minsrcpos, maxsrcpos); + binsSrcPos_offY.SetEdges(nbins_srcpos, minsrcpos, maxsrcpos); + plist_off.AddToList(&binsSrcPos_offX); + plist_off.AddToList(&binsSrcPos_offY); + + MH3 hDAQEvtNumber_off("MRawEvtHeader.fDAQEvtNumber"); + hDAQEvtNumber_off.SetName("DAQEvtNumber_off"); + plist_off.AddToList(&hDAQEvtNumber_off); + MBinning binsDAQEvtNumber_offX("BinningDAQEvtNumber_offX"); + Int_t nbins_evtnum = 100; + Float_t minevtnum = 0.; + Float_t maxevtnum = 100.; + binsDAQEvtNumber_offX.SetEdges(nbins_evtnum,minevtnum,maxevtnum); + plist_off.AddToList(&binsDAQEvtNumber_offX); + + //tasks + MReadTree read_off("Parameters", f_off_name); + read_off.DisableAutoScheme(); + + srccalc.SetMode(MSrcPosFromFile::kOff); + + MHillasSrcCalc csrc_off; + + // fill all histograms + MFillH falpha_off_abs(&hAlpha_off_abs); + MFillH falpha_off(&hAlpha_off); + MFillH fdist_off(&hDist_off); + MFillH fwidth_off(&hWidth_off); + MFillH flength_off(&hLength_off); + MFillH fsize_off(&hSize_off); + MFillH fsrcpos_off(&hSrcPos_off); + MFillH fevtnum_off(&hDAQEvtNumber_off); + + //tasklist + tlist_off.AddToList(&read_off); + tlist_off.AddToList(&srccalc); + tlist_off.AddToList(&csrc_off); + tlist_off.AddToList(&fsrcpos_off); + tlist_off.AddToList(&cont_even); + tlist_off.AddToList(&cont_size); + tlist_off.AddToList(&cont_width); + tlist_off.AddToList(&cont_length); + tlist_off.AddToList(&cont_dist); + tlist_off.AddToList(&cont_alpha); + tlist_off.AddToList(&falpha_off_abs); + tlist_off.AddToList(&falpha_off); + tlist_off.AddToList(&fdist_off); + tlist_off.AddToList(&fwidth_off); + tlist_off.AddToList(&flength_off); + tlist_off.AddToList(&fsize_off); + tlist_off.AddToList(&fevtnum_off); + + // Create and setup the eventloop + MEvtLoop loop_off; + loop_off.SetParList(&plist_off); + //loop_off.SetDisplay(display); + + MProgressBar bar_off; + loop_off.SetProgressBar(&bar_off); + + if (!loop_off.Eventloop(numEntries)) + return; + + tlist_off.PrintStatistics(); + + // ############################################################################ + // look for the histograms + // ############################################################################ + + TH1F *hist_size_on = (TH1F*)hSize_on.GetHist(); + TH1F *hist_size_off = (TH1F*)hSize_off.GetHist(); + + TH1F *hist_dist_on = (TH1F*)hDist_on.GetHist(); + TH1F *hist_dist_off = (TH1F*)hDist_off.GetHist(); + + TH1F *hist_width_on = (TH1F*)hWidth_on.GetHist(); + TH1F *hist_width_off = (TH1F*)hWidth_off.GetHist(); + + TH1F *hist_length_on = (TH1F*)hLength_on.GetHist(); + TH1F *hist_length_off = (TH1F*)hLength_off.GetHist(); + + TH1F *hist_on_abs = (TH1F*)hAlpha_on_abs.GetHist(); + TH1F *hist_off_abs = (TH1F*)hAlpha_off_abs.GetHist(); + + TH1F *hist_on = (TH1F*)hAlpha_on.GetHist(); + TH1F *hist_off = (TH1F*)hAlpha_off.GetHist(); + + + // ############################################################################ + // Calculate significance and excess: + // ############################################################################ + + Double_t norm_on_abs = (Double_t) hist_on_abs->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90+1); + Double_t exces_on_abs = (Double_t) hist_on_abs->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90+1); + Double_t norm_off_abs = (Double_t) hist_off_abs->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90+1); + Double_t exces_off_abs = (Double_t) hist_off_abs->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90+1); + Double_t norm = norm_on_abs/norm_off_abs; + + char text_tit_alpha[256]; + sprintf(text_tit_alpha, " Alpha Plot On and Off "); + hist_off_abs->SetTitle(text_tit_alpha); + hist_on_abs->SetTitle(text_tit_alpha); + + Double_t excess = exces_on_abs - exces_off_abs*norm; + Double_t sign = excess / sqrt( exces_on_abs + norm*norm*exces_off_abs ); + Double_t int_off = (Double_t) hist_off_abs->Integral(1, 18); + int hist_on_entries = (int) hist_on_abs->GetEntries(); + int hist_off_entries = (int) hist_off_abs->GetEntries(); + + cout << "---> Normalization F factor =\t" << norm < Excess =\t\t\t" << excess < Significancia =\t\t" << sign < entries on =\t\t" << hist_on_entries < entries off =\t\t" << hist_off_entries < integral off =\t\t" << int_off < from now on, all histos are plotted + MStatusDisplay *display = new MStatusDisplay; + display->SetUpdateTime(3000); + display->Resize(850,700); + + // ############################################################################ + // Draw SIZE + // ############################################################################ + display->AddTab("SIZE"); + + gPad->cd(); + + gPad->SetLogy(); + hist_size_on->Sumw2(); + hist_size_off->Sumw2(); + hist_size_off->Scale(norm); + hist_size_on->SetLineColor(kBlack); + hist_size_on->SetMarkerStyle(21); + hist_size_on->SetMarkerSize(0.7); + hist_size_on->SetMarkerColor(kBlack); + hist_size_off->SetFillColor(46); + hist_size_off->SetLineColor(46); + hist_size_off->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines + hist_size_off->SetMinimum(0.1); + hist_size_on->SetMinimum(0.1); + hist_size_on->SetTitle("SIZE distribution"); + hist_size_off->SetTitle("SIZE distribution"); + + hist_size_on->DrawCopy("E1P"); + + // move stat box to make them all visible + gPad->Update(); + TPaveStats* pavs_on_size = (TPaveStats*) hist_size_on->GetListOfFunctions()->FindObject("stats"); + if(pavs_on_size){ + shiftx = pavs_on_size->GetX2NDC() - pavs_on_size->GetX1NDC(); + pavs_on_size->SetX1NDC(pavs_on_size->GetX1NDC() - shiftx); + pavs_on_size->SetX2NDC(pavs_on_size->GetX2NDC() - shiftx); + } + gPad->Modified(); + gPad->Update(); + + hist_size_off->DrawCopy("HISTSAME"); + hist_size_off->DrawCopy("ESAME"); + + gPad->Modified(); + gPad->Update(); + + Double_t chisize = ChiSquareNDof((TH1D*)hist_size_on,(TH1D*)hist_size_off); + + Double_t x_label_pos = log10(1000000)*0.7; + Double_t y_label_pos = log10((hist_size_on->GetBinContent(hist_size_on->GetMaximumBin()))/2.); + Double_t textsize = 0.03; + + char text_size[256]; + sprintf(text_size,"ChiSquare/NDof = %4.2f",chisize); + + TLatex *tsize = new TLatex(x_label_pos, y_label_pos, text_size); + tsize->SetTextSize(textsize); +// tsize->Draw(); + + gPad->Modified(); + gPad->Update(); + + // ############################################################################ + // DrawCopy DIST + // ############################################################################ + display->AddTab("DIST"); + + gPad->cd(); + + hist_dist_on->Sumw2(); + hist_dist_off->Sumw2(); + hist_dist_off->Scale(norm); + hist_dist_on->SetLineColor(kBlack); + hist_dist_on->SetMarkerStyle(21); + hist_dist_on->SetMarkerSize(0.7); + hist_dist_on->SetMarkerColor(kBlack); + hist_dist_off->SetFillColor(46); + hist_dist_off->SetLineColor(46); + hist_dist_off->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines + hist_dist_off->SetMinimum(0.); + hist_dist_on->SetTitle("DIST distribution"); + hist_dist_off->SetTitle("DIST distribution"); + + hist_dist_on->DrawCopy("E1P"); + + // move stat box to make them all visible + gPad->Update(); + TPaveStats* pavs_on_dist = (TPaveStats*) hist_dist_on->GetListOfFunctions()->FindObject("stats"); + if(pavs_on_dist){ + shiftx = pavs_on_dist->GetX2NDC() - pavs_on_dist->GetX1NDC(); + pavs_on_dist->SetX1NDC(pavs_on_dist->GetX1NDC() - shiftx); + pavs_on_dist->SetX2NDC(pavs_on_dist->GetX2NDC() - shiftx); + } + gPad->Modified(); + gPad->Update(); + + hist_dist_off->DrawCopy("HISTSAME"); + hist_dist_off->DrawCopy("ESAME"); + hist_dist_on->DrawCopy("E1PSAME"); + + Double_t chidist = ChiSquareNDof((TH1D*)hist_dist_on,(TH1D*)hist_dist_off); + + x_label_pos = distmax*0.7; + y_label_pos = hist_dist_on->GetBinContent(hist_dist_on->GetMaximumBin())/2.; + + char text_dist[256]; + sprintf(text_size,"ChiSquare/NDof = %4.2f",chidist); + + TLatex *tdist = new TLatex(x_label_pos, y_label_pos, text_dist); + tdist->SetTextSize(textsize); +// tdist->Draw(); + + gPad->Modified(); + gPad->Update(); + + // ############################################################################ + // DrawCopy WIDTH + // ############################################################################ + display->AddTab("WIDTH"); + + gPad->cd(); + + hist_width_off->Sumw2(); + hist_width_off->Scale(norm); + hist_width_on->SetLineColor(kBlack); + hist_width_on->SetMarkerStyle(21); + hist_width_on->SetMarkerSize(0.7); + hist_width_on->SetMarkerColor(kBlack); + hist_width_off->SetFillColor(46); + hist_width_off->SetLineColor(46); + hist_width_off->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines + hist_width_off->SetMinimum(0.); + hist_width_on->SetTitle("WIDTH distribution"); + hist_width_off->SetTitle("WIDTH distribution"); + + hist_width_on->DrawCopy("E1P"); + + // move stat box to make them all visible + gPad->Update(); + TPaveStats* pavs_on_width = (TPaveStats*) hist_width_on->GetListOfFunctions()->FindObject("stats"); + if(pavs_on_width){ + shiftx = pavs_on_width->GetX2NDC() - pavs_on_width->GetX1NDC(); + pavs_on_width->SetX1NDC(pavs_on_width->GetX1NDC() - shiftx); + pavs_on_width->SetX2NDC(pavs_on_width->GetX2NDC() - shiftx); + } + gPad->Modified(); + gPad->Update(); + + hist_width_off->DrawCopy("HISTSAME"); + hist_width_off->DrawCopy("ESAME"); + hist_width_on->DrawCopy("E1PSAME"); + + Double_t chiwidth = ChiSquareNDof((TH1D*)hist_width_on,(TH1D*)hist_width_off); + + x_label_pos = widthmax*0.7; + y_label_pos = hist_width_on->GetBinContent(hist_width_on->GetMaximumBin())/2.; + + char text_width[256]; + sprintf(text_size,"ChiSquare/NDof = %4.2f",chiwidth); + + TLatex *twidth = new TLatex(x_label_pos, y_label_pos, text_width); + twidth->SetTextSize(textsize); +// twidth->Draw(); + + gPad->Modified(); + gPad->Update(); + + // ############################################################################ + // DrawCopy LENGTH + // ############################################################################ + display->AddTab("LENGTH"); + + gPad->cd(); + + hist_length_on->Sumw2(); + hist_length_off->Sumw2(); + hist_length_off->Scale(norm); + hist_length_on->SetLineColor(kBlack); + hist_length_on->SetMarkerStyle(21); + hist_length_on->SetMarkerSize(0.7); + hist_length_on->SetMarkerColor(kBlack); + hist_length_off->SetFillColor(46); + hist_length_off->SetLineColor(46); + hist_length_off->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines + hist_length_off->SetMinimum(0.); + hist_length_on->SetTitle("LENGTH distribution"); + hist_length_off->SetTitle("LENGTH distribution"); + + hist_length_on->DrawCopy("E1P"); + + // move stat box to make them all visible + gPad->Update(); + TPaveStats* pavs_on_length = (TPaveStats*) hist_length_on->GetListOfFunctions()->FindObject("stats"); + if(pavs_on_length){ + shiftx = pavs_on_length->GetX2NDC() - pavs_on_length->GetX1NDC(); + pavs_on_length->SetX1NDC(pavs_on_length->GetX1NDC() - shiftx); + pavs_on_length->SetX2NDC(pavs_on_length->GetX2NDC() - shiftx); + } + gPad->Modified(); + gPad->Update(); + + hist_length_off->DrawCopy("HISTSAME"); + hist_length_off->DrawCopy("ESAME"); + hist_length_on->DrawCopy("E1PSAME"); + + Double_t chilength = ChiSquareNDof((TH1D*)hist_length_on,(TH1D*)hist_length_off); + + x_label_pos = lengthmax*0.7; + y_label_pos = hist_length_on->GetBinContent(hist_length_on->GetMaximumBin())/2.; + + char text_length[256]; + sprintf(text_size,"ChiSquare/NDof = %4.2f",chilength); + + TLatex *tlength = new TLatex(x_label_pos, y_label_pos, text_length); + tlength->SetTextSize(textsize); +// tlength->Draw(); + + gPad->Modified(); + gPad->Update(); + + // ############################################################################ + // DrawCopy normalized ALPHA plot + // ############################################################################ + display->AddTab("ALPHA"); + + gPad->cd(); + + hist_on_abs->Sumw2(); + hist_off_abs->SetStats(0); + hist_off_abs->Sumw2(); + hist_off_abs->Scale(norm); + hist_on_abs->SetStats(0); //-> Do NOT show the legend with statistics + hist_on_abs->SetLineColor(kBlack); + hist_on_abs->SetMarkerStyle(21); + //hist_on_abs->SetMarkerSize(); + hist_on_abs->SetMarkerColor(kBlack); + hist_on_abs->SetMarkerSize(0.7); + hist_off_abs->SetFillColor(46); + hist_off_abs->SetLineColor(46); + hist_off_abs->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines + hist_off_abs->SetMinimum(0.); + hist_on_abs->SetTitle("Alpha plot"); + hist_off_abs->SetTitle("Alpha plot"); + + + hist_on_abs->DrawCopy("E1P"); + hist_off_abs->DrawCopy("HISTSAME"); + hist_off_abs->DrawCopy("ESAME"); + hist_on_abs->DrawCopy("E1PSAME"); + + + //draw the LEGEND with excess and significance values in the alpha plot: + char text_Fnorm[256], text_excess[256], text_sign[256]; + char text_entries_on[256], text_entries_off[256], text_integral_off[256]; + int hist_on_entries = (int) hist_on_abs->GetEntries(); + int hist_off_entries = (int) hist_off_abs->GetEntries(); + sprintf(text_Fnorm, " F norm = %.3f", norm); + sprintf(text_excess, " Excess = %.3f", excess); + sprintf(text_sign, " Significance = %.3f", sign); + sprintf(text_entries_on, " Entries ON = %d", hist_on_entries); + sprintf(text_entries_off, " Entries OFF = %d", hist_off_entries); + sprintf(text_integral_off," Integral OFF = %d", int_off); + + x_label_pos = alphamax*0.7; + y_label_pos = (hist_on_abs->GetBinContent(hist_on_abs->GetMaximumBin()))/1.6; //2.; + Double_t y_label_step = y_label_pos / 8.; + + TLatex *t0 = new TLatex(x_label_pos, y_label_pos - y_label_step*0, text_Fnorm); + t0->SetTextSize(textsize); + t0->Draw(); + TLatex *t1 = new TLatex(x_label_pos, y_label_pos - y_label_step*1, text_excess); + t1->SetTextSize(textsize); + t1->Draw(); + TLatex *t2 = new TLatex(x_label_pos, y_label_pos - y_label_step*2, text_sign); + t2->SetTextSize(textsize); + t2->Draw(); + TLatex *t3 = new TLatex(x_label_pos, y_label_pos - y_label_step*3, text_entries_on); + t3->SetTextSize(textsize); + t3->Draw(); + TLatex *t4 = new TLatex(x_label_pos, y_label_pos - y_label_step*4, text_entries_off); + t4->SetTextSize(textsize); + t4->Draw(); + TLatex *t5 = new TLatex(x_label_pos, y_label_pos - y_label_step*5, text_integral_off); + t5->SetTextSize(textsize); + t5->Draw(); + + + Double_t chialpha = ChiSquareNDof((TH1D*)hist_on_abs,(TH1D*)hist_off_abs); + + y_label_pos = (hist_on_abs->GetBinContent(hist_on_abs->GetMaximumBin()))/2.; + + char text_alpha[256]; + sprintf(text_size,"ChiSquare/NDof = %4.2f",chialpha); + + TLatex *talpha = new TLatex(x_label_pos, y_label_pos, text_alpha); + talpha->SetTextSize(textsize); +// talpha->Draw(); + + gPad->Modified(); + gPad->Update(); + + // ############################################################################ + // DrawCopy normalized alpha histos for alpha form -90 to 90 deg. + // ############################################################################ + display->AddTab("ALPHA +-90"); + + gPad->cd(); + + hist_on->Sumw2(); + hist_off->SetStats(0); + hist_off->Sumw2(); + hist_off->Scale(norm); + hist_off->SetFillColor(46); + hist_off->SetLineColor(46); + hist_off->SetFillStyle(3004); //(1001)-> To set the pad NOT transparent and solid; (3004)-> pattern lines + hist_off->SetMinimum(0.); + hist_on->SetStats(0); //-> Do NOT show the legend with statistics + hist_on->SetLineColor(kBlack); + hist_on->SetMarkerStyle(21); + hist_on->SetMarkerSize(0.7); + hist_on->SetMarkerColor(kBlack); + hist_on->SetTitle("Alpha plot form -90 to 90 deg"); + hist_off->SetTitle("Alpha plot form -90 to 90 deg"); + + hist_on->DrawCopy("E1P"); + hist_off->DrawCopy("HISTSAME"); + hist_off->DrawCopy("ESAME"); + hist_on->DrawCopy("E1PSAME"); + + Double_t chialpha90 = ChiSquareNDof((TH1D*)hist_on,(TH1D*)hist_off); + + x_label_pos = alphamax*0.5; + y_label_pos = hist_on->GetBinContent(hist_on->GetMaximumBin())/2.; + + char text_alpha90[256]; + sprintf(text_alpha90,"ChiSquare/NDof = %4.2f",chialpha90); + + TLatex *talpha90 = new TLatex(x_label_pos, y_label_pos, text_alpha90); + talpha90->SetTextSize(textsize); +// talpha90->Draw(); + + gPad->Update(); + gPad->Modified(); + + + cout << "---> ChiSquare/NDof [Size] =\t\t" << chisize << endl; + cout << "---> ChiSquare/NDof [Dist] =\t\t" << chidist << endl; + cout << "---> ChiSquare/NDof [Width] =\t\t" << chiwidth << endl; + cout << "---> ChiSquare/NDof [Length] =\t\t" << chilength << endl; + cout << "---> ChiSquare/NDof [Abs(Alpha)] =\t" << chialpha << endl; + cout << "---> ChiSquare/NDof [Alpha] =\t\t" << chialpha90 << endl; + + + display->AddTab("SRCPOS ON"); + TH2F *hist_srcpos_on = (TH2F*)hSrcPos_on.GetHist(); + hist_srcpos_on->DrawCopy("BOX"); + + display->AddTab("SRCPOS OFF"); + TH2F *hist_srcpos_off = (TH2F*)hSrcPos_off.GetHist(); + hist_srcpos_off->DrawCopy("BOX"); + + display->AddTab("EVTNUM ON"); + TH1F *hist_evtnum_on = (TH1F*)hDAQEvtNumber_on.GetHist(); + hist_evtnum_on->DrawCopy(); + + display->AddTab("EVTNUM OFF"); + TH1F *hist_evtnum_off = (TH1F*)hDAQEvtNumber_off.GetHist(); + hist_evtnum_off->DrawCopy(); + + cout << "Done!!" < +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + + +Bool_t HandleInput() +{ + TTimer timer("gSystem->ProcessEvents();", 50, kFALSE); + while (1) + { + // + // While reading the input process gui events asynchronously + // + timer.TurnOn(); + TString input = Getline("Type 'q' to exit, to go on: "); + timer.TurnOff(); + + if (input=="q\n") + return kFALSE; + + if (input=="\n") + return kTRUE; + }; + + return kFALSE; +} + + +void psffit(const TString filename="dc_2004_02_15_01_51_22_17286_Mrk421.root", const TString directory="/nfs/magic/CaCodata/2004_02_15/", const UInt_t numEvents = 100000) +{ + + // + // Create a empty Parameter List and an empty Task List + // The tasklist is identified in the eventloop by its name + // + MParList plist; + + MTaskList tlist; + plist.AddToList(&tlist); + + + MGeomCamMagic geomcam; + MCameraDC dccam; + MPSFFit psffit; + + plist.AddToList(&geomcam); + plist.AddToList(&dccam); + plist.AddToList(&psffit); + + // + // Now setup the tasks and tasklist: + // --------------------------------- + // + + // Reads the trees of the root file and the analysed branches + MReadReports read; + read.AddTree("Currents"); + read.AddFile(directory+filename); // after the reading of the trees!!! + read.AddToBranchList("MReportCurrents.*"); + + MGeomApply geomapl; + + const Int_t numrings = 3; + const Int_t numblind = 23; + const Short_t x[numblind] = { 8, 27, 224, 279, 339, + 507, 508, 509, 510, 511, 512, 513, 514, + 543, + 559, 560, 561, 562, 563, 564, 565, 566, 567}; + const TArrayS blindpixels(numblind,(Short_t*)x); + MPSFFitCalc psfcalc; + //psfcalc.SetImgCleanMode(MPSFFitCalc::kRing); + psfcalc.SetImgCleanMode(MPSFFitCalc::kCombined); + psfcalc.SetNumRings(numrings); + psfcalc.SetBlindPixels(blindpixels); + + + tlist.AddToList(&geomapl); + tlist.AddToList(&read); + tlist.AddToList(&psfcalc, "Currents"); + + // + // Create and setup the eventloop + // + MEvtLoop evtloop; + evtloop.SetParList(&plist); + + // + // Execute your analysis + // + + if (numEvents > 0) + { + if (!evtloop.Eventloop(numEvents)) + return; + } + else + { + if (!evtloop.PreProcess()) + return; + + MHCamera display(geomcam); + display.SetPrettyPalette(); + display.Draw(); + gPad->cd(1); + psffit.Draw(); + + while (tlist.Process()) + { + display.SetCamContent(dccam); + gPad->Modified(); + gPad->Update(); + // Remove the comments if you want to go through the file + // event-by-event: + if (!HandleInput()) + break; + } + + evtloop.PostProcess(); + } + + tlist.PrintStatistics(); + + psffit.Print(); + cout << "RUN " << psffit.GetMeanMinorAxis() << ' ' << psffit.GetSigmaMinorAxis() << ' ' << psffit.GetMeanMajorAxis() << ' ' << psffit.GetSigmaMajorAxis() << ' ' << psffit.GetChisquare() << endl; + +} + + Index: trunk/MagicSoft/Mars/mtemp/mifae/psffit.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/psffit.cc (revision 3900) +++ trunk/MagicSoft/Mars/mtemp/mifae/psffit.cc (revision 3900) @@ -0,0 +1,147 @@ +/* ======================================================================== *\ +! +! * +! * 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 López, 04/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +#include + +#include +#include + +#include "MParList.h" +#include "MTaskList.h" +#include "MGeomCamMagic.h" +#include "MCameraDC.h" +#include "MPSFFit.h" +#include "MReadReports.h" +#include "MGeomApply.h" +#include "MPSFFitCalc.h" +#include "MEvtLoop.h" + +int main(int argc, char *argv[]) +{ + if(argc!=2 && argc!=3) + { + printf("\n usage: %s filename numberEnvents \n\n",argv[0]); + return 1; + } + + UInt_t len=strlen(argv[1]); + + if(argv[1][len] == '/') + { + argv[1][len]='\0'; + } + + TString filename; + filename=argv[1]; + + UInt_t numEvents = 1000000; + if(argc == 3) + { + len=strlen(argv[2]); + + if(argv[2][len] == '/') + { + argv[2][len]='\0'; + } + + numEvents=atoi(argv[2]); + } + + // + // Create a empty Parameter List and an empty Task List + // The tasklist is identified in the eventloop by its name + // + MParList plist; + + MTaskList tlist; + plist.AddToList(&tlist); + + + MGeomCamMagic geomcam; + MCameraDC dccam; + MPSFFit psffit; + + plist.AddToList(&geomcam); + plist.AddToList(&dccam); + plist.AddToList(&psffit); + + // + // Now setup the tasks and tasklist: + // --------------------------------- + // + + // Reads the trees of the root file and the analysed branches + MReadReports read; + read.AddTree("Currents"); + read.AddFile(filename); // after the reading of the trees!!! + read.AddToBranchList("MReportCurrents.*"); + + MGeomApply geomapl; + + const Int_t numrings = 3; + const Int_t numblind = 23; + const Short_t x[numblind] = { 8, 27, 224, 279, 339, + 507, 508, 509, 510, 511, 512, 513, 514, + 543, + 559, 560, 561, 562, 563, 564, 565, 566, 567}; + + // 2004_02_15 + /* const Int_t numblind = 28; + const Short_t x[numblind] = { 8, 224, 279, 339, + 433, 434, 435, 436, 437, 438, 439, + 475, 476, 477, 478, 479, 480, 481, 482, + 523, 524, 525, 526, 527, 528, 529, 530, 531}; + */ + const TArrayS blindpixels(numblind,(Short_t*)x); + MPSFFitCalc psfcalc; + //psfcalc.SetImgCleanMode(MPSFFitCalc::kRing); + psfcalc.SetImgCleanMode(MPSFFitCalc::kCombined); + psfcalc.SetNumRings(numrings); + psfcalc.SetBlindPixels(blindpixels); + + + tlist.AddToList(&geomapl); + tlist.AddToList(&read); + tlist.AddToList(&psfcalc, "Currents"); + + // + // Create and setup the eventloop + // + MEvtLoop evtloop; + evtloop.SetParList(&plist); + + // + // Execute your analysis + // + + if (!evtloop.Eventloop(numEvents)) + return kFALSE; + + // tlist.PrintStatistics(); + + cout << "RUN " << psffit.GetMeanMinorAxis() << ' ' << psffit.GetSigmaMinorAxis() << ' ' << psffit.GetMeanMajorAxis() << ' ' << psffit.GetSigmaMajorAxis() << ' ' << psffit.GetChisquare() << endl; + +} + + Index: trunk/MagicSoft/Mars/mtemp/mifae/srcPosRun.sh =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/srcPosRun.sh (revision 3900) +++ trunk/MagicSoft/Mars/mtemp/mifae/srcPosRun.sh (revision 3900) @@ -0,0 +1,28 @@ +#!/bin/bash + +DCFILESPATH=/nfs/magic/CaCodata/2004_03_19 +SOURCE=Mrk421 +BIN=./psffit + +DCFILES=$(ls ${DCFILESPATH}/*${SOURCE}*.root) +#echo DCFILES ${DCFILES} + +echo RUN SRCX SRCY +for DCFILE in ${DCFILES} +do + DCFILE=${DCFILE//${DCFILESPATH}/.} + #echo DCFILE ${DCFILE} + RUN=${DCFILE:25:5} + #echo RUN ${RUN} + if [ RUN != '00000' ]; then + SRCPOS=$(./${BIN} ${DCFILESPATH}/${DCFILE} | tail -1) + #echo SRCPOS ${SRCPOS} + SRCX=$(echo ${SRCPOS} | gawk '{print $2}') + SIGX=$(echo ${SRCPOS} | gawk '{print $3}') + SRCY=$(echo ${SRCPOS} | gawk '{print $4}') + SIGY=$(echo ${SRCPOS} | gawk '{print $5}') + CHI2=$(echo ${SRCPOS} | gawk '{print $6}') + #echo ${RUN} ${SRCX} ${SIGX} ${SRCY} ${SIGY} ${CHI2} + echo ${RUN} ${SRCX} ${SRCY} + fi +done