Index: /trunk/MagicSoft/Mars/mtemp/mifae/Changelog =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/Changelog (revision 3946) +++ /trunk/MagicSoft/Mars/mtemp/mifae/Changelog (revision 3947) @@ -18,14 +18,25 @@ -*-*- END OF LINE -*-*- + + 2004/05/03: Javier Rico + * library,programs,macros,*: + - new directories added + - directory macros contains macros (.CC) and scripts (.sh) + - directory programs contains root compilables and associated + [input] files + - directory library contains routines to be linked with libmars.so 2004/04/30: Javier López * alpha_plot.C - - macro to display hillas disttributions with corrected miss pointing + - macro to display hillas disttributions with corrected + miss pointing * srcPosRun.sh - - Script to produce the file with the positions of UMa51 using psffit executable + - 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 + - Tool to correct the position of the source for on and off data + from a file * psffit.C & psffit.cc @@ -33,5 +44,6 @@ * MPSFFit.[h,cc] & MPSFFitCalc.[h,cc] - - Container and taks to do calculation of the psf and position of stars. + - Container and taks to do calculation of the psf and position of + stars. 2004/04/28: Javier Rico Index: unk/MagicSoft/Mars/mtemp/mifae/IFAEIncl.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/IFAEIncl.h (revision 3946) +++ (revision ) @@ -1,3 +1,0 @@ -#ifndef __CINT__ - -#endif // __CINT__ Index: unk/MagicSoft/Mars/mtemp/mifae/IFAELinkDef.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/IFAELinkDef.h (revision 3946) +++ (revision ) @@ -1,12 +1,0 @@ -#ifdef __CINT__ - -#pragma link off all globals; -#pragma link off all classes; -#pragma link off all functions; - -#pragma link C++ class MDCA+; -#pragma link C++ class MPSFFit+; -#pragma link C++ class MPSFFitCalc+; -#pragma link C++ class MSrcPosFromFile+; - -#endif Index: unk/MagicSoft/Mars/mtemp/mifae/MDCA.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/MDCA.cc (revision 3946) +++ (revision ) @@ -1,261 +1,0 @@ - -///////////////////////////////////////////////////////////////////////////// -// // -// MDCA.cc // -// Container to store Hillas parameters and DCA stuff // -// // -// Author(s): S.C. Commichau, L.S. Stark, 7/2003 // -// // -///////////////////////////////////////////////////////////////////////////// - -#include "MDCA.h" - -using namespace std; - -ClassImp(MDCA); - -// Default constructor -MDCA::MDCA(const char *name, const char *title) -{ - fName = name ? name : "MDCA"; - fTitle = title ? title : "Storage container for Hillas parameters and DCA of one event"; - - Reset(); - - fEllipse = new TEllipse; - fRefCircle = new TEllipse; - - fLineL = new TLine; - fLineW = new TLine; - fLineX = new TLine; - fLineY = new TLine; - fLineDCA = new TLine; - fLineMean = new TLine; - -} - -// Destructor: Deletes ellipse and lines if they do exist - -MDCA::~MDCA() -{ - Clear(); -} - -// Initialize parameters with default values - -void MDCA::Reset() -{ - - fLength = -1; - fWidth = -1; - fDelta0 = 0; - fMeanX = 0; - fMeanY = 0; - fDelta1 = 0; - fDCA = -1000; - fX1W = 0; - fY1W = 0; - fX2W = 0; - fY2W = 0; - fX1L = 0; - fY1L = 0; - fX2L = 0; - fY2L = 0; - fXDCA = 0; - fYDCA = 0; - -} - - -// Print parameters to *fLog - -void MDCA::Print(Option_t *) const -{ - Double_t atg = atan2(fMeanY, fMeanX)*kRad2Deg; - - if (atg<0) - atg += 180; - - *fLog << all; - *fLog << "Basic Image Parameters (" << GetName() << ")" << endl; - *fLog << " - Length [mm] = " << fLength << endl; - *fLog << " - Width [mm] = " << fWidth << endl; - *fLog << " - Delta0 [deg] = " << fDelta0*kRad2Deg << endl; - *fLog << " - Meanx [mm] = " << fMeanX << endl; - *fLog << " - Meany [mm] = " << fMeanY << endl; - *fLog << " - atg(y/x) [deg] = " << atg << endl; - *fLog << " - DCA [mm] = " << fDCA << endl; - *fLog << " - Delta1 [deg] = " << fDelta1*kRad2Deg << endl; - - -} - -void MDCA::Paint(Option_t *opt) -{ - Clear(); - - if (fLength<=0 || fWidth<=0) //fLength<0 || fWidth<0 doesn't look nice... - return; //We get a circle! - - // Length line - fLineL = new TLine(fX1L, fY1L, fX2L, fY2L); - fLineL->SetLineWidth(2); - fLineL->SetLineColor(2); - fLineL->Draw(); - - // Width line - fLineW = new TLine(fX1W, fY1W, fX2W, fY2W); - fLineW->SetLineWidth(2); - fLineW->SetLineColor(2); - fLineW->Draw(); - - // Coordinate system - fLineX = new TLine(-600,fYRef,600,fYRef); - fLineY = new TLine(fXRef,-600,fXRef,600); - fLineX->SetLineWidth(1); - fLineX->SetLineColor(108); - fLineY->SetLineWidth(1); - fLineY->SetLineColor(108); - fLineX->Draw(); - fLineY->Draw(); - - // DCA line - fLineDCA = new TLine(fXRef,fYRef,fXDCA+fXRef,fYDCA+fYRef); - fLineDCA->SetLineWidth(2); - fLineDCA->SetLineColor(2); - fLineDCA->Draw(); - - // COG line - fLineMean = new TLine(fXRef,fYRef,fMeanX,fMeanY); - fLineMean->SetLineWidth(2); - fLineMean->SetLineColor(2); - fLineMean->Draw(); - - // Reference point marker - fRefCircle = new TEllipse(fXRef, fYRef, 10, 10, 0, 360, 0); - fRefCircle->SetLineColor(108); - fRefCircle->SetFillColor(108); - fRefCircle->Draw(); - - // Hillas ellipse - fEllipse = new TEllipse(fMeanX, fMeanY, fLength, fWidth, 0, 360, fDelta0*kRad2Deg+180); - fEllipse->SetLineWidth(2); - fEllipse->SetLineColor(2); - fEllipse->Draw(); - -} - - -// If an ellipse and lines exist they will be deleted - -void MDCA::Clear(Option_t *) -{ - if (!fEllipse && !fRefCircle && !fLineL && !fLineW && !fLineX && !fLineY && !fLineDCA && !fLineMean) - return; - - delete fEllipse; - delete fRefCircle; - delete fLineL; - delete fLineW; - delete fLineX; - delete fLineY; - delete fLineDCA; - delete fLineMean; - - fLineL = NULL; - fLineX = NULL; - fLineY = NULL; - fLineW = NULL; - fLineDCA = NULL; - fLineMean = NULL; - fEllipse = NULL; - fRefCircle = NULL; -} - - -Int_t MDCA::Calc(const MGeomCam &geom, const MCerPhotEvt &evt, const MHillas &hil) -{ - // Get basic Hillas parameters from MHillas - fDelta0 = hil.GetDelta(); - fMeanX = hil.GetMeanX(); - fMeanY = hil.GetMeanY(); - fLength = hil.GetLength(); - fWidth = hil.GetWidth(); - - // The Length Line - rotation and shift - fX1L = - (fLength+OffsetL)*cos(fDelta0) + fMeanX; // [mm] - fY1L = - (fLength+OffsetL)*sin(fDelta0) + fMeanY; // [mm] - fX2L = (fLength+OffsetL)*cos(fDelta0) + fMeanX; // [mm] - fY2L = (fLength+OffsetL)*sin(fDelta0) + fMeanY; // [mm] - - // The Width Line - rotation and shift - fX1W = (fWidth+OffsetW)*sin(fDelta0) + fMeanX; // [mm] - fY1W = - (fWidth+OffsetW)*cos(fDelta0) + fMeanY; // [mm] - fX2W = - (fWidth+OffsetW)*sin(fDelta0) + fMeanX; // [mm] - fY2W = (fWidth+OffsetW)*cos(fDelta0) + fMeanY; // [mm] - - // Vector of orientation of the shower axis - fr1 = fX2L-fX1L; - fr2 = fY2L-fY1L; - - // Determine parameters to calculate coordinates of the DCA vector - flambda = (fr1*(fXRef-fMeanX) + fr2*(fYRef-fMeanY))/(fr1*fr1 + fr2*fr2); - fmu = (fMeanY-fYRef)/fr1 + flambda*fr2/fr1; - - // Components of the DCA vector - fXDCA = -fmu*fr2; - fYDCA = fmu*fr1; - - // Components of vector going from intersection point of the DCA vector - // with the shower axis to the COG - fd1 = fMeanX + fmu*fr2 - fXRef; - fd2 = fMeanY - fmu*fr1 - fYRef; - - // Calculate DCA value - fDCA = sqrt(fXDCA*fXDCA + fYDCA*fYDCA); - - // Calculate angle of the shower axis with respect to the x-axis - fDelta1 = acos(fd1/sqrt(fd1*fd1 + fd2*fd2)); //botch - - // Calculate angle of the shower axis with respect to the y-axis - //fDelta1 = acos(fd2/sqrt(fd1*fd1 + fd2*fd2)); - - // Determine the correct sign of the DCA (cross product of DCA vector and - // vector going from the intersection point of the DCA vector with the shower axis - // to the COG) - if((fmu*(-fr2*(fMeanY-fYRef)-fr1*(fMeanX-fXRef)))<0) - fDCA = -fDCA; - - gRandom->Rannor(gx,gy); - gx = fabs(gx); - - // This is nice but does not remove the systematics in the profile plot... - //if(((1-0.6*gx)*(180-kRad2Deg*fDelta1)>120) || ((1-0.6*gx)*kRad2Deg*fDelta1>120)) - // fDCA = -1; - - // Enlarge the interval of Delta correctly... - if((fMeanY-fYRef-fmu*fr1)<0) - fDelta1 = TwoPi-fDelta1; - - // Enlarge the interval of Delta correctly... (Delta with respect to the y-axis) - // if(-(fMeanX-fXRef+fmu*fr2)<0) - // fDelta1 = TwoPi-fDelta1; - - // This has to be improved... - if(fr1 == 0 || fr2 == 0 || (fr1*fr1+fr2*fr2) == 0 || sqrt(fd1*fd1+fd2*fd2) == 0) - fDCA = -fDCA; - - SetReadyToSave(); - - return 0; -} - -void MDCA::SetRefPoint(const Float_t fXRef0, const Float_t fYRef0) -{ - fXRef = fXRef0; - fYRef = fYRef0; -} - - - - Index: unk/MagicSoft/Mars/mtemp/mifae/MDCA.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/MDCA.h (revision 3946) +++ (revision ) @@ -1,118 +1,0 @@ -#ifndef MARS_MDCA -#define MARS_MDCA - -///////////////////////////////////////////////////////////////////////////// -// // -// MDCA.cc // -// Container to store the DCA stuff // -// // -// Author(s): S.C. Commichau, L.S. Stark, 7/2003 // -// // -// // -///////////////////////////////////////////////////////////////////////////// - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include "MLog.h" -#include "MLogManip.h" -#include "MGeomCam.h" -#include "MGeomPix.h" -#include "MCerPhotPix.h" -#include "MCerPhotEvt.h" -#include "MHillas.h" - -// Add offsets to the Width and Length to draw nice lines -#define OffsetW 20.0 -#define OffsetL 440.0 - -#define kmm2Deg 0.00317460317 -#define TwoPi 6.28318530717958623 - -class TArrayF; -class TEllipse; -class MGeoCam; -class MCerPhotEvt; -class MHillas; - -class MDCA : public MParContainer -{ - private: - - Float_t fLength; // [mm] major axis of ellipse - Float_t fWidth; // [mm] minor axis of ellipse - Float_t fDelta0; // [rad] angle of major axis with x-axis (-Pi/2 .. Pi/2) - Float_t fDelta1; // [rad] my angle of major axis with x-axis (0 .. 2*Pi) - Float_t fMeanX; // [mm] x-coordinate of center of ellipse - Float_t fMeanY; // [mm] y-coordinate of center of ellipse - Float_t fX1W; // [mm] x-coordinate of 1st point of Lline - Float_t fY1W; // [mm] y-coordinate of 1st point of Lline - Float_t fX2W; // [mm] x-coordinate of 2nd point of Lline - Float_t fY2W; // [mm] y-coordinate of 2nd point of Lline - Float_t fX1L; // [mm] x-coordinate of 1st point of Wline - Float_t fY1L; // [mm] y-coordinate of 1st point of Wline - Float_t fX2L; // [mm] x-coordinate of 2nd point of Wline - Float_t fY2L; // [mm] y-coordinate of 2nd point of Wline - Double_t fDCA; // [mm] Distance of Cloasest Approach - Float_t fXDCA; // [mm] x-coordinate of 2nd point of the DCA-line - Float_t fYDCA; // [mm] y-coordinate of 2nd point of the DCA-line - Float_t fXRef; // [mm] x-coordinate of reference point - Float_t fYRef; // [mm] y-coordinate of reference point - Float_t fmu; - Float_t flambda; - Float_t fr1, fr2; // [mm] Coordinates of the orientation vector of the shower axis - Float_t fd1, fd2; // [mm] Coordinates of the DCA vector - Float_t gx, gy; - - TEllipse *fEllipse; // Graphical object to draw the ellipse, ROOT! - TEllipse *fRefCircle; // To draw reference point - - TLine *fLineL; // Shower axis - TLine *fLineW; // Line perpendicular to the shower axis - TLine *fLineX; // x-axis of the coordinate system - TLine *fLineY; // y-axis of the coordinate system - TLine *fLineDCA; // DCA line - TLine *fLineMean; // Line to COG of the shower - - public: - MDCA(const char *name=NULL, const char *title=NULL); - - ~MDCA(); - - void Reset(); - - Int_t Calc(const MGeomCam &geom, const MCerPhotEvt &pix, const MHillas &hil); - - void Print(Option_t *opt=NULL) const; - - void Paint(Option_t *opt=NULL); - - void Clear(Option_t *opt=NULL); - - void SetRefPoint(const Float_t fXRef0, const Float_t fYRef0); - - Float_t GetDCA() const { return fDCA; } - - Float_t GetDelta() const { return fDelta1; } - - ClassDef(MDCA, 1) - -}; - -#endif - - - - - Index: unk/MagicSoft/Mars/mtemp/mifae/MPSFFit.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/MPSFFit.cc (revision 3946) +++ (revision ) @@ -1,86 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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: unk/MagicSoft/Mars/mtemp/mifae/MPSFFit.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/MPSFFit.h (revision 3946) +++ (revision ) @@ -1,45 +1,0 @@ -#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: unk/MagicSoft/Mars/mtemp/mifae/MPSFFitCalc.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/MPSFFitCalc.cc (revision 3946) +++ (revision ) @@ -1,569 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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: unk/MagicSoft/Mars/mtemp/mifae/MPSFFitCalc.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/MPSFFitCalc.h (revision 3946) +++ (revision ) @@ -1,93 +1,0 @@ -///////////////////////////////////////////////////////////////////////////// -// // -// 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: unk/MagicSoft/Mars/mtemp/mifae/MSrcPosFromFile.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/MSrcPosFromFile.cc (revision 3946) +++ (revision ) @@ -1,248 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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: unk/MagicSoft/Mars/mtemp/mifae/MSrcPosFromFile.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/MSrcPosFromFile.h (revision 3946) +++ (revision ) @@ -1,55 +1,0 @@ -#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, 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: unk/MagicSoft/Mars/mtemp/mifae/Makefile =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/Makefile (revision 3946) +++ (revision ) @@ -1,80 +1,0 @@ -################################################################## -# -# makefile -# -# for the MARS IFAE library -# -################################################################## -# @maintitle - -# @code - -# -# please change all system depend values in the -# config.mk.${OSTYPE} file -# -# -include ../../Makefile.conf.$(OSTYPE) -include ../../Makefile.conf.general - -# -CINT = IFAE -# - -PROGRAMS = makeHillas psffit -#SOLIB = ../../libmars.so -#LIB = mifae.a - -#------------------------------------------------------------------------------ - -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../../mastro - - -#.SUFFIXES: .c .cc .h .o - -SRCFILES = \ - MDCA.cc \ - MPSFFit.cc \ - MPSFFitCalc.cc \ - MSrcPosFromFile.cc -# MSrcRotate.cc - -#SRCS = $(SRCFILES) -#HEADERS = $(SRCFILES:.cc=.h) -#OBJS = $(SRCFILES:.cc=.o) - -############################################################ -#all: $(LIB) $(PROGRAMS) -all: $(OBJS) - @echo " Done. " - @echo " " - -include ../../Makefile.rules - - -$(PROGRAMS): $(SOLIB) $(LIB) $(PROGRAMS:=.o) - @echo " Linking $@ ..." - $(CXX) $(CXXFLAGS) $(ROOTGLIBS) $(SOLIB) ${OBJS} -L../../lib -lmifae $@.o $(MARS_LIB) -o $@ - @echo " Removing $@.o ..." - rm $@.o - -mrproper: rmbin clean rmbak - -# @endcode Index: unk/MagicSoft/Mars/mtemp/mifae/alpha_plot.C =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/alpha_plot.C (revision 3946) +++ (revision ) @@ -1,834 +1,0 @@ -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 +! +! +\* ======================================================================== */ +#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/library/MPSFFit.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/library/MPSFFit.h (revision 3947) +++ /trunk/MagicSoft/Mars/mtemp/mifae/library/MPSFFit.h (revision 3947) @@ -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/library/MPSFFitCalc.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/library/MPSFFitCalc.cc (revision 3947) +++ /trunk/MagicSoft/Mars/mtemp/mifae/library/MPSFFitCalc.cc (revision 3947) @@ -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/library/MPSFFitCalc.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/library/MPSFFitCalc.h (revision 3947) +++ /trunk/MagicSoft/Mars/mtemp/mifae/library/MPSFFitCalc.h (revision 3947) @@ -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/library/MSrcPosFromFile.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcPosFromFile.cc (revision 3947) +++ /trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcPosFromFile.cc (revision 3947) @@ -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/library/MSrcPosFromFile.h =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcPosFromFile.h (revision 3947) +++ /trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcPosFromFile.h (revision 3947) @@ -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, 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/library/MSrcRotate.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcRotate.cc (revision 3947) +++ /trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcRotate.cc (revision 3947) @@ -0,0 +1,202 @@ +/* ======================================================================== *\ +! +! * +! * 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 Rico 04/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MSrcRotate +// +// Task to rotate the position of the source as a function of Azimuth and +// Zenith angles +// +// Input Containers: +// MSrcPosCam +// +// Output Containers: +// MSrcPosCam +// MDCA +// +////////////////////////////////////////////////////////////////////////////// + +#include +#include + +#include "MParList.h" +#include "MRawRunHeader.h" +#include "MRawEvtHeader.h" +#include "MSrcRotate.h" +#include "MSrcPosCam.h" +#include "MDCA.h" + +#include "MObservatory.h" +#include "MAstroSky2Local.h" +#include "MAstro.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MSrcRotate); + +using namespace std; + +static const TString gsDefName = "MSrcRotate"; +static const TString gsDefTitle = "De-rotate position of the source"; + +// ------------------------------------------------------------------------- +// +// Default constructor. The first (second) argument is the name of a container +// containing the source position in the camera plain, MScrPosCam (MDCA). +// The default is "MSrcPosCam" ("MDCA"). +// +MSrcRotate::MSrcRotate(const char* srcPos, const char* dca, const char *name, const char *title) + : fSrcPos(NULL), fDCA(NULL), fRA(0), fDEC(0), fRunNumber(0) +{ + fName = name ? name : gsDefName.Data(); + fTitle = title ? title : gsDefTitle.Data(); + + fSrcPosName = srcPos; + fDCAName = dca; +} + +// ------------------------------------------------------------------------- +// +// Look for/create the needed containers +// Check whether RA and DEC have been initialized +// +Int_t MSrcRotate::PreProcess(MParList *pList) +{ + // look for/create MSrcPosCam + fSrcPos = (MSrcPosCam*)pList->FindObject(AddSerialNumber(fSrcPosName), "MSrcPosCam"); + if (!fSrcPos) + { + *fLog << warn << AddSerialNumber(fSrcPosName) << " [MSrcPosCam] not found... creating default container." << endl; + fSrcPos = (MSrcPosCam*)pList->FindCreateObj("MSrcPosCam", AddSerialNumber(fSrcPosName)); + if(!fSrcPos) + return kFALSE; + } + fIniSrcPos = *fSrcPos; + + // look for/create MDCA + fDCA = (MDCA*)pList->FindObject(AddSerialNumber(fDCAName), "MDCA"); + if (!fDCA) + { + *fLog << warn << AddSerialNumber(fDCAName) << " [MDCA] not found... creating default container." << endl; + fDCA = (MDCA*)pList->FindCreateObj("MDCA", AddSerialNumber(fDCAName)); + if(!fDCA) + return kFALSE; + } + + // look for MRawRunHeader + fRunHeader = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); + if (!fRunHeader) + { + *fLog << err << "MSrcRotate::PreProcess Error: MRawRunHeader not found... aborting" << endl; + return kFALSE; + } + + // look for MRawEvtHeader + fEvtHeader = (MRawEvtHeader*)pList->FindObject("MRawEvtHeader"); + if (!fEvtHeader) + { + *fLog << err << "MSrcRotate::PreProcess Error: MRawEvtHeader not found... aborting." << endl; + return kFALSE; + } + + // look for/create the MObservatory + fObservatory = (MObservatory*)pList->FindCreateObj("MObservatory"); + if(!fObservatory) + { + *fLog << err << "MSrcRotate::PreProcess Error: MObservatory not found... aborting." << endl; + return kFALSE; + } + + // check RA and DEC + if(fRA==0. || fDEC ==0) + *fLog << warn << "MSrcRotate::PreProcess Warning: RA=" << fRA << "or DEC=" << fDEC << " could be not initialized" << endl; + + return kTRUE; +} + +// ------------------------------------------------------------------------- +// +// Perform the rotation of the source position +// +// FIXME: for the time being, this is computed by assuming constant event rate +// +Int_t MSrcRotate::Process() +{ + if(fRunHeader->GetRunNumber()!=fRunNumber) + { + fRunNumber=fRunHeader->GetRunNumber(); + + // save the number of events, initial and final times + fNEvts = fRunHeader->GetNumEvents(); + fFirstEvt = fEvtHeader->GetDAQEvtNumber(); + fIniTime = fRunHeader->GetRunStart(); + fFinTime = fRunHeader->GetRunEnd(); + fDeltaT = (fFinTime.GetMjd()-fIniTime.GetMjd())/fNEvts; + + if((ULong_t)round(fIniTime.GetMjd())!=(ULong_t)round(fFinTime.GetMjd())) + { + *fLog << err << "MSrcRotate::Process Error: Inial and final MJDs are different ("<GetRunNumber() << endl; + cout << "Number of events: " << fNEvts << endl; + cout << "First event: " << fFirstEvt << endl; + cout << "Initial MJD date: " << fIniTime.GetMjd() << endl; + cout << "Final MJD date: " << fFinTime.GetMjd() << endl; + cout << "Delta t: " << fDeltaT << endl; +#endif + } + + + // Compute the event time + // FIXME: for the time being, this is computed by assuming constant event rate + MTime eventTime; + Double_t newMJD = fIniTime.GetMjd() + (fFinTime.GetMjd()-fIniTime.GetMjd())*(fEvtHeader->GetDAQEvtNumber()-fFirstEvt)/fNEvts; + eventTime.SetMjd(newMJD); + + // de-rotate the source position + Double_t rotationAngle = fObservatory->RotationAngle(fRA,fDEC,eventTime); + + Float_t c = TMath::Cos(rotationAngle); + Float_t s = TMath::Sin(rotationAngle); + // perform a rotation of -rotationAngle to move the source back to the "initial" position + Float_t newX = c*fIniSrcPos.GetX()+s*fIniSrcPos.GetY(); + Float_t newY = -s*fIniSrcPos.GetX()+c*fIniSrcPos.GetY(); + +#ifdef DEBUG + Double_t rotationAngleComp = fObservatory->RotationAngle(0.1256,2.63); + cout << "Event " << fEvtHeader->GetDAQEvtNumber() << endl; + cout << "newMJD=" << newMJD <<", rotationAngle=" << rotationAngle <<", rotationAngleComp=" << rotationAngleComp << ", oldX="< +#include "MTime.h" +#include "MSrcPosCam.h" + +class MDCA; +class MObservatory; +class MRawEvtHeader; +class MRawRunHeader; + +class MSrcRotate : public MTask +{ +private: + MSrcPosCam* fSrcPos; //! Pointer to the source position + MDCA* fDCA; //! Pointer to the MDCA object + MRawEvtHeader* fEvtHeader; //! Pointer to the event header + MRawRunHeader* fRunHeader; //! Pointer to the run header + MObservatory* fObservatory; //! Pointer to the MObservatory + + TString fSrcPosName; + TString fDCAName; + + MSrcPosCam fIniSrcPos; // Copy of the initial source position + Double_t fRA; // [rad] Right ascention + Double_t fDEC; // [rad] Declination + UInt_t fNEvts; // Number of events in file + UInt_t fFirstEvt; // Event index for first event in file + MTime fIniTime; // Run initial time + MTime fFinTime; // Run final time + Double_t fDeltaT; // DeltaT between two events + UInt_t fRunNumber; // Current run number + + Int_t PreProcess(MParList *plist); + Int_t Process(); + +public: + MSrcRotate(const char* src="MSrcPosCam", const char* dca="MDCA", + const char* name=NULL, const char* title=NULL); + + void SetRAandDEC(Double_t ra, Double_t dec) {fRA=ra;fDEC=dec;}; + + ClassDef(MSrcRotate, 0) // task to rotate the position of the source as a function of Azimuth and Zenith angles +}; + +#endif + Index: /trunk/MagicSoft/Mars/mtemp/mifae/library/Makefile =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/library/Makefile (revision 3947) +++ /trunk/MagicSoft/Mars/mtemp/mifae/library/Makefile (revision 3947) @@ -0,0 +1,64 @@ +################################################################## +# +# makefile +# +# for the MARS IFAE library +# +################################################################## +# @maintitle + +# @code + +# +# please change all system depend values in the +# config.mk.${OSTYPE} file +# +# +include ../../../Makefile.conf.$(OSTYPE) +include ../../../Makefile.conf.general + +# +CINT = IFAE +# + +SOLIB = ../../../libmars.so + +#------------------------------------------------------------------------------ + +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../../../mastro + +SRCFILES = \ + MDCA.cc \ + MPSFFit.cc \ + MPSFFitCalc.cc \ + MSrcPosFromFile.cc +# MSrcRotate.cc + + +############################################################ + +all: $(OBJS) + +include ../../../Makefile.rules + +clean: rmcint rmobjs rmcore rmlib +mrproper: clean rmbak + +# @endcode Index: /trunk/MagicSoft/Mars/mtemp/mifae/macros/alpha_plot.C =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/macros/alpha_plot.C (revision 3947) +++ /trunk/MagicSoft/Mars/mtemp/mifae/macros/alpha_plot.C (revision 3947) @@ -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/macros/srcPosRun.sh =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/macros/srcPosRun.sh (revision 3947) +++ /trunk/MagicSoft/Mars/mtemp/mifae/macros/srcPosRun.sh (revision 3947) @@ -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 Index: unk/MagicSoft/Mars/mtemp/mifae/makeHillas.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/makeHillas.cc (revision 3946) +++ (revision ) @@ -1,417 +1,0 @@ -/*********************************/ -/* Compute the hillas parameters */ -/*********************************/ - -#include "TString.h" -#include "TArrayS.h" - -#include "MParList.h" -#include "MTaskList.h" -#include "MPedestalCam.h" -#include "MBadPixelsCam.h" -#include "MReadMarsFile.h" -#include "MGeomApply.h" -#include "MPedCalcPedRun.h" -#include "MEvtLoop.h" -#include "MGeomCamMagic.h" -#include "MExtractedSignalCam.h" -#include "MCalibrationChargeCam.h" -#include "MHCalibrationChargeCam.h" -#include "MHCalibrationRelTimeCam.h" -#include "MExtractFixedWindow.h" -#include "MCalibrationChargeCalc.h" -#include "MFCosmics.h" -#include "MContinue.h" -#include "MFillH.h" -#include "MLog.h" -#include "MCerPhotEvt.h" -#include "MPedPhotCam.h" -#include "MCalibrate.h" -#include "MPedPhotCalc.h" -#include "MHillas.h" -#include "MRawRunHeader.h" -#include "MSrcPosCam.h" -#include "MBlindPixelCalc.h" -#include "MImgCleanStd.h" -#include "MHillasSrcCalc.h" -#include "MHillasCalc.h" -#include "MWriteRootFile.h" -#include "MProgressBar.h" -#include "MArgs.h" -#include "MRunIter.h" -#include "MJPedestal.h" -#include "MJCalibration.h" - -#include -#include -#include - -using namespace std; - -Bool_t readDatacards(TString& filename); -void makeHillas(); - -// initial and final time slices to be used in signal extraction -const Byte_t hifirst = 1; -const Byte_t hilast = 14; -const Byte_t lofirst = 4; -const Byte_t lolast = 13; - -// declaration of variables read from datacards -TString outname; -TString idirname; -MRunIter caliter; -MRunIter pediter; -MRunIter datiter; -ULong_t nmaxevents=999999999; -Short_t calflag=1; -Float_t lcore = 3.0; -Float_t ltail = 1.5; -Int_t nfiles = 0; - -const TString defaultcard="input.datacard"; - -/*************************************************************/ -// makeHillas usage output -static void Usage() -{ - gLog <" << endl << endl; - gLog << " : datacards file name (dafault input.datacards)" << endl; - gLog << " -?/-h: This help" << endl << endl; -} - -/*************************************************************/ -// main program -int main(int argc, char **argv) -{ - // evaluate arguments - MArgs arg(argc, argv); - if (arg.HasOnly("-?") || arg.HasOnly("-h") || arg.HasOnly("--help")) - { - Usage(); - return -1; - } - - // get name of input datacard file - TString datacard = arg.GetArgumentStr(0); - if(!datacard.Length()) - datacard = defaultcard; - - // read the datacards - if(!readDatacards(datacard)) - { - cout << "Error reading datacards. Stoping" << endl; - return -1; - } - - // make the hillas file - makeHillas(); -} - -/*************************************************************/ -/*************************************************************/ -/*************************************************************/ -void makeHillas() -{ - - // set the signal extractor and calibration mode - MExtractFixedWindow extractor; - extractor.SetRange(hifirst,hilast,lofirst,lolast); - - MCalibrate::CalibrationMode_t calMode=MCalibrate::kDefault; - if(calflag==0) - calMode=MCalibrate::kNone; - if(calflag==-1) - calMode=MCalibrate::kDummy; - MCalibrate calibrate(calMode); - - - // general containers - MBadPixelsCam badcam; - MGeomCamMagic geomcam; - MGeomApply geomapl; - - // If you want to exclude pixels from the beginning, read - // an ascii-file with the corr. pixel numbers (see MBadPixelsCam) - // - // badcam.AsciiRead("badpixels.dat"); - - - /************************************/ - /* FIRST LOOP: PEDESTAL COMPUTATION */ - /************************************/ - - MJPedestal pedloop; - pedloop.SetInput(&pediter); - pedloop.SetBadPixels(badcam); - - cout << "*************************" << endl; - cout << "** COMPUTING PEDESTALS **" << endl; - cout << "*************************" << endl; - - if (!pedloop.Process()) - return; - - /*****************************/ - /* SECOND LOOP: CALIBRATION */ - /*****************************/ - MCalibrationQECam qecam; - MJCalibration calloop; - calloop.SetInput(&caliter); - calloop.SetExtractor(&extractor); - // - // Set the corr. cams: - // - calloop.SetQECam(qecam); - calloop.SetBadPixels(pedloop.GetBadPixels()); - - // - // Apply rel. time calibration: - // calloop.SetRelTimeCalibration(); - - // Use as arrival time extractor MArrivalTimeCalc2: - // calloop.SetArrivalTimeLevel(2); - - cout << "***************************" << endl; - cout << "** COMPUTING CALIBRATION **" << endl; - cout << "***************************" << endl; - - if (!calloop.Process(pedloop.GetPedestalCam())) - return; - - - /************************************************************************/ - /* THIRD LOOP: PEDESTAL CALIBRATION INTO PHOTONS */ - /************************************************************************/ - MParList plist3; - MTaskList tlist3; - plist3.AddToList(&tlist3); - - // containers - MCerPhotEvt nphot; - MPedPhotCam nphotrms; - - plist3.AddToList(&geomcam); - plist3.AddToList(&pedloop.GetPedestalCam()); - plist3.AddToList(&calloop.GetCalibrationCam()); - plist3.AddToList(&calloop.GetQECam()); - plist3.AddToList(&calloop.GetRelTimeCam()); - plist3.AddToList(&calloop.GetBadPixels()); - plist3.AddToList(&nphot); - plist3.AddToList(&nphotrms); - - //tasks - MReadMarsFile read3("Events"); - read3.DisableAutoScheme(); - static_cast(read3).AddFiles(pediter); - - MPedPhotCalc photrmscalc; - - tlist3.AddToList(&read3); - tlist3.AddToList(&geomapl); - tlist3.AddToList(&extractor); - tlist3.AddToList(&calibrate); - tlist3.AddToList(&photrmscalc); - - // Create and setup the eventloop - MEvtLoop evtloop3; - evtloop3.SetParList(&plist3); - if (!evtloop3.Eventloop()) - return; - - tlist3.PrintStatistics(); - - /************************************************************************/ - /* FOURTH LOOP: DATA CALIBRATION INTO PHOTONS */ - /************************************************************************/ - - MParList plist4; - MTaskList tlist4; - plist4.AddToList(&tlist4); - - // containers - MHillas hillas; - MSrcPosCam srcposcam; - MRawRunHeader runhead; - - plist4.AddToList(&geomcam); - plist4.AddToList(&pedloop.GetPedestalCam()); - plist4.AddToList(&calloop.GetCalibrationCam()); - plist4.AddToList(&calloop.GetQECam()); - plist4.AddToList(&calloop.GetRelTimeCam()); - plist4.AddToList(&calloop.GetBadPixels()); - plist4.AddToList(&nphot); - plist4.AddToList(&nphotrms); - plist4.AddToList(&srcposcam); - plist4.AddToList(&hillas); - plist4.AddToList(&runhead); - - //tasks - MReadMarsFile read4("Events"); - read4.DisableAutoScheme(); - static_cast(read4).AddFiles(datiter); - - MImgCleanStd clean(lcore,ltail); - MHillasCalc hcalc; - MHillasSrcCalc csrc1; - - MWriteRootFile write(outname,"RECREATE"); - - write.AddContainer("MHillas" , "Parameters"); - write.AddContainer("MHillasSrc" , "Parameters"); - write.AddContainer("MHillasExt" , "Parameters"); - write.AddContainer("MNewImagePar" , "Parameters"); - write.AddContainer("MRawEvtHeader" , "Parameters"); - write.AddContainer("MRawRunHeader" , "Parameters"); - write.AddContainer("MConcentration" , "Parameters"); - write.AddContainer("MSrcPosCam", "Parameters"); - - tlist4.AddToList(&read4); - tlist4.AddToList(&geomapl); - tlist4.AddToList(&extractor); - tlist4.AddToList(&calibrate); - tlist4.AddToList(&clean); - tlist4.AddToList(&hcalc); - tlist4.AddToList(&csrc1); - tlist4.AddToList(&write); - - // Create and setup the eventloop - MEvtLoop datloop; - datloop.SetParList(&plist4); - - cout << "*************************************************************" << endl; - cout << "*** COMPUTING DATA USING EXTRACTED SIGNAL (IN PHOTONS) ***" << endl; - cout << "*************************************************************" << endl; - - if (!datloop.Eventloop(nmaxevents)) - return; - - tlist4.PrintStatistics(); -} - -/******************************************************************************/ -Bool_t readDatacards(TString& filename) -{ - ifstream ifun(filename.Data()); - if(!ifun) - { - cout << "File " << filename << " not found" << endl; - return kFALSE; - } - - TString word; - - while(ifun >> word) - { - // skip comments - if(word[0]=='/' && word[1]=='/') - { - while(ifun.get()!='\n'); // skip line - continue; - } - - // number of events - if(strcmp(word.Data(),"NEVENTS")==0) - ifun >> nmaxevents; - - - // input file directory - if(strcmp(word.Data(),"IDIR")==0) - { - if(idirname.Length()) - cout << "readDataCards Warning: overriding input directory file name" << endl; - ifun >> idirname; - } - - // pedestal runs - if(strcmp(word.Data(),"PRUNS")==0) - { - if(pediter.GetNumRuns()) - cout << "readDataCards Warning: adding pedestal runs to the existing list" << endl; - ifun >> word; - pediter.AddRuns(word.Data(),idirname.Data()); - } - - // calibration runs - if(strcmp(word.Data(),"CRUNS")==0) - { - if(caliter.GetNumRuns()) - cout << "readDataCards Warning: adding calibration runs to the existing list" << endl; - ifun >> word; - caliter.AddRuns(word.Data(),idirname.Data()); - } - - // data runs - if(strcmp(word.Data(),"DRUNS")==0) - { - if(datiter.GetNumRuns()) - cout << "readDataCards Warning: adding data runs to the existing list" << endl; - ifun >> word; - datiter.AddRuns(word.Data(),idirname.Data()); - } - - // output file name - if(strcmp(word.Data(),"OUTFILE")==0) - { - if(outname.Length()) - cout << "readDataCards Warning: overriding output file name" << endl; - ifun >> outname; - } - - // calibration flag - if(strcmp(word.Data(),"CALFLAG")==0) - ifun >> calflag; - - // cleaning level - if(strcmp(word.Data(),"CLEANLEVEL")==0) - { - ifun >> lcore; - ifun >> ltail; - } - } - - pediter.Reset(); - caliter.Reset(); - datiter.Reset(); - TString pfile; - - // Dump read values - cout << "************************************************" << endl; - cout << "* Datacards read from file " << filename << endl; - cout << "************************************************" << endl; - cout << "Pedestal file (s): " << endl; - while(!(pfile=pediter.Next()).IsNull()) - cout << pfile << endl; - cout << "Calibration file (s): " << endl; - while(!(pfile=caliter.Next()).IsNull()) - cout << pfile << endl; - cout << "Data file (s): " << endl; - while(!(pfile=datiter.Next()).IsNull()) - cout << pfile << endl; - cout << "Maximum number of events: " << nmaxevents << endl; - cout << "Output file name: " << outname << endl; - cout << "Calibration flag: " << calflag << endl; - cout << "Cleaning level: ("< +#include +#include + +using namespace std; + +Bool_t readDatacards(TString& filename); +void makeHillas(); + +// initial and final time slices to be used in signal extraction +const Byte_t hifirst = 1; +const Byte_t hilast = 14; +const Byte_t lofirst = 4; +const Byte_t lolast = 13; + +// declaration of variables read from datacards +TString outname; +TString idirname; +MRunIter caliter; +MRunIter pediter; +MRunIter datiter; +ULong_t nmaxevents=999999999; +Short_t calflag=1; +Float_t lcore = 3.0; +Float_t ltail = 1.5; +Int_t nfiles = 0; + +const TString defaultcard="input.datacard"; +static void Usage() +{ + gLog <" << endl << endl; + gLog << " : datacards file name (dafault input.datacards)" << endl; + gLog << " -?/-h: This help" << endl << endl; +} + +int main(int argc, char **argv) +{ + // evaluate arguments + MArgs arg(argc, argv); + if (arg.HasOnly("-?") || arg.HasOnly("-h") || arg.HasOnly("--help")) + { + Usage(); + return -1; + } + + TString datacard = arg.GetArgumentStr(0); + if(!datacard.Length()) + datacard = defaultcard; + + if(!readDatacards(datacard)) + { + cout << "Error reading datacards. Stoping" << endl; + return -1; + } + makeHillas(); +} + +/*************************************************************/ +void makeHillas() +{ +#if 0 + MStatusDisplay *display = new MStatusDisplay; + display->SetUpdateTime(3000); + display->Resize(850,700); + + gStyle->SetOptStat(1111); + gStyle->SetOptFit(); +#endif + + /************************************/ + /* FIRST LOOP: PEDESTAL COMPUTATION */ + /************************************/ + + MParList plist1; + MTaskList tlist1; + plist1.AddToList(&tlist1); + + // containers + MPedestalCam pedcam; + MBadPixelsCam badcam; + // badcam.AsciiRead("badpixels.dat"); + + plist1.AddToList(&pedcam); + plist1.AddToList(&badcam); + + //tasks + MReadMarsFile read("Events"); + static_cast(read).AddFiles(pediter); + read.DisableAutoScheme(); + + MGeomApply geomapl; + MPedCalcPedRun pedcalc; + + tlist1.AddToList(&read); + tlist1.AddToList(&geomapl); + tlist1.AddToList(&pedcalc); + + // Create and setup the eventloop + MEvtLoop pedloop; + pedloop.SetParList(&plist1); + //pedloop.SetDisplay(display); + + cout << "*************************" << endl; + cout << "** COMPUTING PEDESTALS **" << endl; + cout << "*************************" << endl; + + if (!pedloop.Eventloop()) + return; + + tlist1.PrintStatistics(); + + /*****************************/ + /* SECOND LOOP: CALIBRATION */ + /*****************************/ + + MParList plist2; + MTaskList tlist2; + plist2.AddToList(&tlist2); + plist2.AddToList(&pedcam); + plist2.AddToList(&badcam); + + // new containers + MGeomCamMagic geomcam; + MExtractedSignalCam sigcam; + MCalibrationChargeCam calcam; + MHCalibrationChargeCam histcharge; + MHCalibrationRelTimeCam histtime; + + + plist2.AddToList(&geomcam); + plist2.AddToList(&sigcam); + plist2.AddToList(&calcam); + plist2.AddToList(&histcharge); + + //tasks + MReadMarsFile read2("Events"); + static_cast(read2).AddFiles(caliter); + read2.DisableAutoScheme(); + + MExtractSignal sigcalc; + sigcalc.SetRange(hifirst,hilast,lofirst,lolast); + MCalibrationChargeCalc calcalc; + MFCosmics cosmics; + MContinue cont(&cosmics); + + + calcalc.SetPulserColor(MCalibrationCam::kCT1); + MFillH fillcam ("MHCalibrationChargeCam" , "MExtractedSignalCam"); + + tlist2.AddToList(&read2); + tlist2.AddToList(&geomapl); + tlist2.AddToList(&sigcalc); + tlist2.AddToList(&cont); + tlist2.AddToList(&fillcam); + tlist2.AddToList(&calcalc); + + // Create and setup the eventloop + MEvtLoop calloop; + calloop.SetParList(&plist2); + //calloop.SetDisplay(display); + + cout << "***************************" << endl; + cout << "** COMPUTING CALIBRATION **" << endl; + cout << "***************************" << endl; + + if (!calloop.Eventloop()) + return; + + tlist2.PrintStatistics(); + + MLog gauglog; + gauglog.SetOutputFile(Form("%s%s",calcam.GetName(),".txt"),1); + calcam.SetLogStream(&gauglog); + calcam.Print(); + calcam.SetLogStream(&gLog); +#if 0 + // Create histograms to display + MHCamera disp1 (geomcam, "Cal;Charge", "Fitted Mean Charges"); + MHCamera disp2 (geomcam, "Cal;SigmaCharge", "Sigma of Fitted Charges"); + MHCamera disp3 (geomcam, "Cal;FitProb", "Probability of Fit"); + MHCamera disp4 (geomcam, "Cal;RSigma", "Reduced Sigmas"); + MHCamera disp5 (geomcam, "Cal;RSigma/Charge", "Reduced Sigma per Charge"); + MHCamera disp6 (geomcam, "Cal;FFactorPh", "Nr. of Photo-electrons (F-Factor Method)"); + MHCamera disp7 (geomcam, "Cal;FFactorConv", "Conversion Factor to photons (F-Factor Method)"); + MHCamera disp8 (geomcam, "Cal;FFactorFFactor", "Total F-Factor (F-Factor Method)"); + MHCamera disp9 (geomcam, "Cal;BlindPixPh", "Photon flux inside plexiglass (Blind Pixel Method)"); + MHCamera disp10 (geomcam, "Cal;BlindPixConv", "Conversion Factor to photons (Blind Pixel Method)"); + MHCamera disp11 (geomcam, "Cal;BlindPixFFactor","Total F-Factor (Blind Pixel Method)"); + MHCamera disp12 (geomcam, "Cal;PINDiodePh", "Photon flux outside plexiglass (PIN Diode Method)"); + MHCamera disp13 (geomcam, "Cal;PINDiodeConv", "Conversion Factor tp photons (PIN Diode Method)"); + MHCamera disp14 (geomcam, "Cal;PINDiodeFFactor","Total F-Factor (PIN Diode Method)"); + MHCamera disp15 (geomcam, "Cal;Excluded", "Pixels previously excluded"); + MHCamera disp16 (geomcam, "Cal;NotFitted", "Pixels that could not be fitted"); + MHCamera disp17 (geomcam, "Cal;NotFitValid", "Pixels with not valid fit results"); + MHCamera disp18 (geomcam, "Cal;HiGainOscillating", "Oscillating Pixels HI Gain"); + MHCamera disp19 (geomcam, "Cal;LoGainOscillating", "Oscillating Pixels LO Gain"); + MHCamera disp20 (geomcam, "Cal;HiGainPickup", "Number Pickup events Hi Gain"); + MHCamera disp21 (geomcam, "Cal;LoGainPickup", "Number Pickup events Lo Gain"); + MHCamera disp22 (geomcam, "Cal;Saturation", "Pixels with saturated Hi Gain"); + MHCamera disp23 (geomcam, "Cal;FFactorValid", "Pixels with valid F-Factor calibration"); + MHCamera disp24 (geomcam, "Cal;BlindPixelValid", "Pixels with valid BlindPixel calibration"); + MHCamera disp25 (geomcam, "Cal;PINdiodeFFactorValid", "Pixels with valid PINDiode calibration"); + + MHCamera disp26 (geomcam, "Cal;Ped", "Pedestals"); + MHCamera disp27 (geomcam, "Cal;PedRms", "Pedestal RMS"); + + MHCamera disp28 (geomcam, "time;Time", "Rel. Arrival Times"); + MHCamera disp29 (geomcam, "time;SigmaTime", "Sigma of Rel. Arrival Times"); + MHCamera disp30 (geomcam, "time;TimeProb", "Probability of Time Fit"); + MHCamera disp31 (geomcam, "time;NotFitValid", "Pixels with not valid fit results"); + MHCamera disp32 (geomcam, "time;Oscillating", "Oscillating Pixels"); + + MHCamera disp33 (geomcam, "Cal;AbsTimeMean", "Abs. Arrival Times"); + MHCamera disp34 (geomcam, "Cal;AbsTimeRms", "RMS of Arrival Times"); + + // Fitted charge means and sigmas + disp1.SetCamContent(calcam, 0); + disp1.SetCamError( calcam, 1); + disp2.SetCamContent(calcam, 2); + disp2.SetCamError( calcam, 3); + + // Fit probabilities + disp3.SetCamContent(calcam, 4); + + // Reduced Sigmas and reduced sigmas per charge + disp4.SetCamContent(calcam, 5); + disp4.SetCamError( calcam, 6); + disp5.SetCamContent(calcam, 7); + disp5.SetCamError( calcam, 8); + + // F-Factor Method + disp6.SetCamContent(calcam, 9); + disp6.SetCamError( calcam, 10); + disp7.SetCamContent(calcam, 11); + disp7.SetCamError( calcam, 12); + disp8.SetCamContent(calcam, 13); + disp8.SetCamError( calcam, 14); + + // Blind Pixel Method + disp9.SetCamContent(calcam, 15); + disp9.SetCamError( calcam, 16); + disp10.SetCamContent(calcam,17); + disp10.SetCamError( calcam,18); + disp11.SetCamContent(calcam,19); + disp11.SetCamError( calcam,20); + + // PIN Diode Method + disp12.SetCamContent(calcam,21); + disp12.SetCamError( calcam,22); + disp13.SetCamContent(calcam,23); + disp13.SetCamError( calcam,24); + disp14.SetCamContent(calcam,25); + disp14.SetCamError( calcam,26); + + // Pixels with defects + disp15.SetCamContent(calcam,27); + disp16.SetCamContent(calcam,28); + disp17.SetCamContent(badcam,9); + disp18.SetCamContent(badcam,16); + disp19.SetCamContent(badcam,15); + disp20.SetCamContent(calcam,29); + disp21.SetCamContent(calcam,30); + + // Lo Gain calibration + disp22.SetCamContent(calcam,31); + + // Valid flags + disp23.SetCamContent(calcam,32); + disp24.SetCamContent(calcam,33); + disp25.SetCamContent(calcam,34); + + // Pedestals + disp26.SetCamContent(calcam,35); + disp26.SetCamError( calcam,36); + disp27.SetCamContent(calcam,37); + disp27.SetCamError( calcam,38); + + // Relative Times + disp28.SetCamContent(histtime,0); + disp28.SetCamError( histtime,1); + disp29.SetCamContent(histtime,2); + disp29.SetCamError( histtime,3); + disp30.SetCamContent(histtime,4); + disp31.SetCamContent(histtime,5); + disp32.SetCamContent(histtime,6); + + // Absolute Times + disp33.SetCamContent(calcam,39); + disp33.SetCamError( calcam,40); + disp34.SetCamContent(calcam,41); + + disp1.SetYTitle("Charge [FADC units]"); + disp2.SetYTitle("\\sigma_{Charge} [FADC units]"); + disp3.SetYTitle("P_{Charge} [1]"); + + disp4.SetYTitle("\\sqrt{\\sigma^{2}_{Charge} - RMS^{2}_{Ped}} [FADC Counts]"); + disp5.SetYTitle("Reduced Sigma / Mean Charge [1]"); + + disp6.SetYTitle("Nr. Photo-electrons [1]"); + disp7.SetYTitle("Conversion Factor [Ph/FADC Count]"); + disp8.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1] "); + + disp9.SetYTitle("Photon flux [ph/mm^2]"); + disp10.SetYTitle("Conversion Factor [Phot/FADC Count]"); + disp11.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]"); + + disp12.SetYTitle("Photon flux [ph/mm^2]"); + disp13.SetYTitle("Conversion Factor [Phot/FADC Count]"); + disp14.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]"); + + disp15.SetYTitle("[1]"); + disp16.SetYTitle("[1]"); + disp17.SetYTitle("[1]"); + disp18.SetYTitle("[1]"); + disp19.SetYTitle("[1]"); + disp20.SetYTitle("[1]"); + disp21.SetYTitle("[1]"); + disp22.SetYTitle("[1]"); + disp23.SetYTitle("[1]"); + disp24.SetYTitle("[1]"); + disp25.SetYTitle("[1]"); + + disp26.SetYTitle("Ped [FADC Counts ]"); + disp27.SetYTitle("RMS_{Ped} [FADC Counts ]"); + + disp28.SetYTitle("Time Offset [ns]"); + disp29.SetYTitle("Timing resolution [ns]"); + disp30.SetYTitle("P_{Time} [1]"); + + disp31.SetYTitle("[1]"); + disp32.SetYTitle("[1]"); + + disp33.SetYTitle("Mean Abs. Time [FADC slice]"); + disp34.SetYTitle("RMS Abs. Time [FADC slices]"); + + gStyle->SetOptStat(1111); + gStyle->SetOptFit(); + + // Charges + TCanvas &c1 = display->AddTab("Fit.Charge"); + c1.Divide(2, 3); + + CamDraw(c1, disp1,calcam,1, 2 , 2); + CamDraw(c1, disp2,calcam,2, 2 , 2); + + // Fit Probability + TCanvas &c2 = display->AddTab("Fit.Prob"); + c2.Divide(1,3); + + CamDraw(c2, disp3,calcam,1, 1 , 4); + + // Reduced Sigmas + TCanvas &c3 = display->AddTab("Red.Sigma"); + c3.Divide(2,3); + + CamDraw(c3, disp4,calcam,1, 2 , 2); + CamDraw(c3, disp5,calcam,2, 2 , 2); + + // F-Factor Method + TCanvas &c4 = display->AddTab("F-Factor"); + c4.Divide(3,3); + + CamDraw(c4, disp6,calcam,1, 3 , 2); + CamDraw(c4, disp7,calcam,2, 3 , 2); + CamDraw(c4, disp8,calcam,3, 3 , 2); + + // Blind Pixel Method + TCanvas &c5 = display->AddTab("BlindPix"); + c5.Divide(3, 3); + + CamDraw(c5, disp9,calcam,1, 3 , 9); + CamDraw(c5, disp10,calcam,2, 3 , 2); + CamDraw(c5, disp11,calcam,3, 3 , 2); + + // PIN Diode Method + TCanvas &c6 = display->AddTab("PINDiode"); + c6.Divide(3,3); + + CamDraw(c6, disp12,calcam,1, 3 , 9); + CamDraw(c6, disp13,calcam,2, 3 , 2); + CamDraw(c6, disp14,calcam,3, 3 , 2); + + // Defects + TCanvas &c7 = display->AddTab("Defects"); + c7.Divide(4,2); + + CamDraw(c7, disp15,calcam,1,4, 0); + CamDraw(c7, disp16,calcam,2,4, 0); + CamDraw(c7, disp20,calcam,3,4, 0); + CamDraw(c7, disp21,calcam,4,4, 0); + + // BadCam + TCanvas &c8 = display->AddTab("Defects"); + c8.Divide(3,2); + + CamDraw(c8, disp17,badcam,1,3, 0); + CamDraw(c8, disp18,badcam,2,3, 0); + CamDraw(c8, disp19,badcam,3,3, 0); + + // Valid flags + TCanvas &c9 = display->AddTab("Validity"); + c9.Divide(4,2); + + CamDraw(c9, disp22,calcam,1,4,0); + CamDraw(c9, disp23,calcam,2,4,0); + CamDraw(c9, disp24,calcam,3,4,0); + CamDraw(c9, disp25,calcam,4,4,0); + + // Pedestals + TCanvas &c10 = display->AddTab("Pedestals"); + c10.Divide(2,3); + + CamDraw(c10,disp26,calcam,1,2,1); + CamDraw(c10,disp27,calcam,2,2,2); + + // Rel. Times + TCanvas &c11 = display->AddTab("Fitted Rel. Times"); + c11.Divide(3,3); + + CamDraw(c11,disp28,calcam,1,3,2); + CamDraw(c11,disp29,calcam,2,3,2); + CamDraw(c11,disp30,calcam,3,3,4); + + // Time Defects + TCanvas &c12 = display->AddTab("Time Def."); + c12.Divide(2,2); + + CamDraw(c12, disp31,calcam,1,2, 0); + CamDraw(c12, disp32,calcam,2,2, 0); + + // Abs. Times + TCanvas &c13 = display->AddTab("Abs. Times"); + c13.Divide(2,3); + + CamDraw(c13,disp33,calcam,1,2,2); + CamDraw(c13,disp34,calcam,2,2,2); +#endif + /************************************************************************/ + /* THIRD LOOP: PEDESTAL CALIBRATION INTO PHOTONS */ + /************************************************************************/ + MParList plist3; + MTaskList tlist3; + plist3.AddToList(&tlist3); + + // containers + MCerPhotEvt nphot; + MPedPhotCam nphotrms; + + plist3.AddToList(&geomcam); + plist3.AddToList(&pedcam); + plist3.AddToList(&calcam); + // plist3.AddToList(&badcam); + plist3.AddToList(&sigcam); + plist3.AddToList(&nphot); + plist3.AddToList(&nphotrms); + + + MCalibrate::CalibrationMode_t calMode=MCalibrate::kDefault; + if(calflag==0) + calMode=MCalibrate::kNone; + + + //tasks + MReadMarsFile read3("Events"); + static_cast(read3).AddFiles(pediter); + read3.DisableAutoScheme(); + + MExtractSignal extsig; + extsig.SetRange(hifirst,hilast,lofirst,lolast); + MCalibrate photcalc(calMode); + MPedPhotCalc photrmscalc; + + tlist3.AddToList(&read3); + tlist3.AddToList(&geomapl); + tlist3.AddToList(&extsig); + tlist3.AddToList(&photcalc); + tlist3.AddToList(&photrmscalc); + + // Create and setup the eventloop + MEvtLoop evtloop3; + evtloop3.SetParList(&plist3); + if (!evtloop3.Eventloop()) + return; + + tlist3.PrintStatistics(); + + /************************************************************************/ + /* FOURTH LOOP: DATA CALIBRATION INTO PHOTONS */ + /************************************************************************/ + + MParList plist4; + MTaskList tlist4; + plist4.AddToList(&tlist4); + + // containers + MHillas hillas; + MSrcPosCam source; + MRawRunHeader runhead; + + plist4.AddToList(&geomcam); + plist4.AddToList(&pedcam); + plist4.AddToList(&calcam); + // plist4.AddToList(&badcam); + plist4.AddToList(&nphot); + plist4.AddToList(&nphotrms); + plist4.AddToList(&source); + plist4.AddToList(&hillas); + plist4.AddToList(&runhead); + + //tasks + MReadMarsFile read4("Events"); + static_cast(read4).AddFiles(datiter); + read4.DisableAutoScheme(); + + // set bad pixels + MBlindPixelCalc blind; + MBlindPixelCalc blind2; + const Short_t x[16] = {330,395,329,396,389, + 323,388,322,384,385, + 386,387,321,320,319, + 394}; + const TArrayS bp(16,(Short_t*)x); + blind.SetPixelIndices(bp); + blind2.SetPixelIndices(bp); + + MImgCleanStd clean(lcore,ltail); + MHillasCalc hcalc; + MHillasSrcCalc csrc1; + + MWriteRootFile write(outname,"RECREATE"); + + // write.AddContainer("MGeomCam" , "RunHeaders"); + // write.AddContainer("MRawRunHeader" , "RunHeaders"); + // write.AddContainer("MSrcPosCam" , "RunHeaders"); + // write.AddContainer("MCalibrationChargeCam" , "RunHeaders"); + // write.AddContainer("MPedPhotCam" , "RunHeaders"); // Attention, was in Events - Tree!! + // write.AddContainer("MPedestalCam" , "RunHeaders"); + // write.AddContainer("MHCalibrationRelTimeCam","RunHeaders"); + + // write.AddContainer("MCerPhotEvt" , "Events"); + // write.AddContainer("MRawEvtHeader" , "Events"); + // write.AddContainer("MBadPixelsCam" , "Events"); + // write.AddContainer("MPedPhotCam" , "Events"); + + write.AddContainer("MHillas" , "Parameters"); + write.AddContainer("MHillasSrc" , "Parameters"); + write.AddContainer("MHillasExt" , "Parameters"); + write.AddContainer("MNewImagePar" , "Parameters"); + write.AddContainer("MRawEvtHeader" , "Parameters"); + write.AddContainer("MRawRunHeader" , "Parameters"); + write.AddContainer("MConcentration" , "Parameters"); + + tlist4.AddToList(&read4); + tlist4.AddToList(&geomapl); + tlist4.AddToList(&extsig); + tlist4.AddToList(&photcalc); + //tlist4.AddToList(&blind); + tlist4.AddToList(&clean); + //tlist4.AddToList(&blind2); + tlist4.AddToList(&hcalc); + // tlist4.AddToList(&srcposcalc); + tlist4.AddToList(&csrc1); + tlist4.AddToList(&write); + + // Create and setup the eventloop + MEvtLoop datloop; + datloop.SetParList(&plist4); + // MProgressBar bar; + // datloop.SetProgressBar(&bar); + + cout << "*************************************************************" << endl; + cout << "*** COMPUTING DATA USING EXTRACTED SIGNAL (IN PHOTONS) ***" << endl; + cout << "*************************************************************" << endl; + + if (!datloop.Eventloop(nmaxevents)) + return; + + tlist4.PrintStatistics(); + nphotrms.Print(); + +} + +#if 0 +void CamDraw(TCanvas &c, MHCamera &cam, MCamEvent &evt, Int_t i, Int_t j, Int_t fit) +{ + + c.cd(i); + gPad->SetBorderMode(0); + MHCamera *obj1=(MHCamera*)cam.DrawCopy("hist"); + // obj1->AddNotify(evt); + + c.cd(i+j); + gPad->SetBorderMode(0); + obj1->Draw(); + ((MHCamera*)obj1)->SetPrettyPalette(); + + if (fit != 0) + { + c.cd(i+2*j); + gPad->SetBorderMode(0); + TH1D *obj2 = (TH1D*)obj1->Projection(obj1.GetName()); + +// obj2->Sumw2(); + obj2->Draw(); + obj2->SetBit(kCanDelete); + + const Double_t min = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst()); + const Double_t max = obj2->GetBinCenter(obj2->GetXaxis()->GetLast()); + const Double_t integ = obj2->Integral("width")/2.5066283; + const Double_t mean = obj2->GetMean(); + const Double_t rms = obj2->GetRMS(); + const Double_t width = max-min; + + if (rms == 0. || width == 0. ) + return; + + switch (fit) + { + case 1: + TF1 *sgaus = new TF1("sgaus","gaus(0)",min,max); + sgaus->SetBit(kCanDelete); + sgaus->SetParNames("Area","#mu","#sigma"); + sgaus->SetParameters(integ/rms,mean,rms); + sgaus->SetParLimits(0,0.,integ); + sgaus->SetParLimits(1,min,max); + sgaus->SetParLimits(2,0,width/1.5); + obj2->Fit("sgaus","QLR"); + obj2->GetFunction("sgaus")->SetLineColor(kYellow); + break; + + case 2: + TString dgausform = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"; + dgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])"; + TF1 *dgaus = new TF1("dgaus",dgausform.Data(),min,max); + dgaus->SetBit(kCanDelete); + dgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}","A_{2}","#mu_{2}","#sigma_{2}"); + dgaus->SetParameters(integ,(min+mean)/2.,width/4., + integ/width/2.,(max+mean)/2.,width/4.); + // The left-sided Gauss + dgaus->SetParLimits(0,integ-1.5,integ+1.5); + dgaus->SetParLimits(1,min+(width/10.),mean); + dgaus->SetParLimits(2,0,width/2.); + // The right-sided Gauss + dgaus->SetParLimits(3,0,integ); + dgaus->SetParLimits(4,mean,max-(width/10.)); + dgaus->SetParLimits(5,0,width/2.); + obj2->Fit("dgaus","QLRM"); + obj2->GetFunction("dgaus")->SetLineColor(kYellow); + break; + + case 3: + TString tgausform = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"; + tgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])"; + tgausform += "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])"; + TF1 *tgaus = new TF1("tgaus",tgausform.Data(),min,max); + tgaus->SetBit(kCanDelete); + tgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}", + "A_{2}","#mu_{2}","#sigma_{2}", + "A_{3}","#mu_{3}","#sigma_{3}"); + tgaus->SetParameters(integ,(min+mean)/2,width/4., + integ/width/3.,(max+mean)/2.,width/4., + integ/width/3.,mean,width/2.); + // The left-sided Gauss + tgaus->SetParLimits(0,integ-1.5,integ+1.5); + tgaus->SetParLimits(1,min+(width/10.),mean); + tgaus->SetParLimits(2,width/15.,width/2.); + // The right-sided Gauss + tgaus->SetParLimits(3,0.,integ); + tgaus->SetParLimits(4,mean,max-(width/10.)); + tgaus->SetParLimits(5,width/15.,width/2.); + // The Gauss describing the outliers + tgaus->SetParLimits(6,0.,integ); + tgaus->SetParLimits(7,min,max); + tgaus->SetParLimits(8,width/4.,width/1.5); + obj2->Fit("tgaus","QLRM"); + obj2->GetFunction("tgaus")->SetLineColor(kYellow); + break; + case 4: + obj2->Fit("pol0","Q"); + obj2->GetFunction("pol0")->SetLineColor(kYellow); + break; + case 9: + break; + default: + obj2->Fit("gaus","Q"); + obj2->GetFunction("gaus")->SetLineColor(kYellow); + break; + } + + TArrayI s0(3); + s0[0] = 6; + s0[1] = 1; + s0[2] = 2; + + TArrayI s1(3); + s1[0] = 3; + s1[1] = 4; + s1[2] = 5; + + TArrayI inner(1); + inner[0] = 0; + + TArrayI outer(1); + outer[0] = 1; + + // Just to get the right (maximum) binning + TH1D *half[4]; + half[0] = obj1->ProjectionS(s0, inner, "Sector 6-1-2 Inner"); + half[1] = obj1->ProjectionS(s1, inner, "Sector 3-4-5 Inner"); + half[2] = obj1->ProjectionS(s0, outer, "Sector 6-1-2 Outer"); + half[3] = obj1->ProjectionS(s1, outer, "Sector 3-4-5 Outer"); + + for (int i=0; i<4; i++) + { + half[i]->SetLineColor(kRed+i); + half[i]->SetDirectory(0); + half[i]->SetBit(kCanDelete); + half[i]->Draw("same"); + } + + gPad->Modified(); + gPad->Update(); + + } +} +#endif + +Bool_t readDatacards(TString& filename) +{ + ifstream ifun(filename.Data()); + if(!ifun) + { + cout << "File " << filename << " not found" << endl; + return kFALSE; + } + + TString word; + + while(ifun >> word) + { + // skip comments + if(word[0]=='/' && word[1]=='/') + { + while(ifun.get()!='\n'); // skip line + continue; + } + + // number of events + if(strcmp(word.Data(),"NEVENTS")==0) + ifun >> nmaxevents; + + + // input file directory + if(strcmp(word.Data(),"IDIR")==0) + { + if(idirname.Length()) + cout << "readDataCards Warning: overriding input directory file name" << endl; + ifun >> idirname; + } + + // pedestal runs + if(strcmp(word.Data(),"PRUNS")==0) + { + if(pediter.GetNumRuns()) + cout << "readDataCards Warning: adding pedestal runs to the existing list" << endl; + ifun >> word; + pediter.AddRuns(word.Data(),idirname.Data()); + } + + // calibration runs + if(strcmp(word.Data(),"CRUNS")==0) + { + if(caliter.GetNumRuns()) + cout << "readDataCards Warning: adding calibration runs to the existing list" << endl; + ifun >> word; + caliter.AddRuns(word.Data(),idirname.Data()); + } + + // data runs + if(strcmp(word.Data(),"DRUNS")==0) + { + if(datiter.GetNumRuns()) + cout << "readDataCards Warning: adding data runs to the existing list" << endl; + ifun >> word; + datiter.AddRuns(word.Data(),idirname.Data()); + } + + // output file name + if(strcmp(word.Data(),"OUTFILE")==0) + { + if(outname.Length()) + cout << "readDataCards Warning: overriding output file name" << endl; + ifun >> outname; + } + + // calibration flag + if(strcmp(word.Data(),"CALFLAG")==0) + ifun >> calflag; + + // cleaning level + if(strcmp(word.Data(),"CLEANLEVEL")==0) + { + ifun >> lcore; + ifun >> ltail; + } + } + + pediter.Reset(); + caliter.Reset(); + datiter.Reset(); + TString pfile; + + // Dump read values + cout << "************************************************" << endl; + cout << "* Datacards read from file " << filename << endl; + cout << "************************************************" << endl; + cout << "Pedestal file (s): " << endl; + while(!(pfile=pediter.Next()).IsNull()) + cout << pfile << endl; + cout << "Calibration file (s): " << endl; + while(!(pfile=caliter.Next()).IsNull()) + cout << pfile << endl; + cout << "Data file (s): " << endl; + while(!(pfile=datiter.Next()).IsNull()) + cout << pfile << endl; + cout << "Maximum number of events: " << nmaxevents << endl; + cout << "Output file name: " << outname << endl; + cout << "Calibration flag: " << calflag << endl; + cout << "Cleaning level: ("< +! +! 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" + +using namespace std; + +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: unk/MagicSoft/Mars/mtemp/mifae/psffit.C =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/psffit.C (revision 3946) +++ (revision ) @@ -1,148 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 -! -! -\* ======================================================================== */ - - -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: unk/MagicSoft/Mars/mtemp/mifae/psffit.cc =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/psffit.cc (revision 3946) +++ (revision ) @@ -1,150 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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" - -using namespace std; - -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: unk/MagicSoft/Mars/mtemp/mifae/srcPosRun.sh =================================================================== --- /trunk/MagicSoft/Mars/mtemp/mifae/srcPosRun.sh (revision 3946) +++ (revision ) @@ -1,28 +1,0 @@ -#!/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