Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MDispParameters.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MDispParameters.cc (revision 7136) +++ (revision ) @@ -1,129 +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): Eva Domingo, 12/2004 -! Wolfgang Wittek, 12/2004 -! -! Copyright: MAGIC Software Development, 2000-2005 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// // -// MDispParameters // -// // -// this is the container for the Disp parameters // -// // -///////////////////////////////////////////////////////////////////////////// -#include "MDispParameters.h" - -#include -#include - -#include - -#include "MImageParDisp.h" - -#include "MParList.h" - -#include "MLog.h" -#include "MLogManip.h" - - -ClassImp(MDispParameters); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// constructor -// -MDispParameters::MDispParameters(const char *name, const char *title) -{ - fName = name ? name : "MDispParameters"; - fTitle = title ? title : "Container for the Disp parameters"; -} - - -// -------------------------------------------------------------------------- -// -// set default values for the Disp parameters -// -void MDispParameters::InitParameters() -{ - //--------------------------------------------------- - // these are the default values - - fParameters[0] = 0.93227; - fParameters[1] = 0.187132; - fParameters[2] = -0.192808; - fParameters[3] = 0.0350143; - fParameters[4] = -0.0110078; - - - //--------------------------------------------------- - // fStepsizes - // if == 0.0 the parameter will be fixed in the minimization - // != 0.0 initial step sizes for the parameters - - fStepsizes[0] = 0.010; - fStepsizes[1] = 0.002; - fStepsizes[2] = 0.002; - fStepsizes[3] = 0.0004; - fStepsizes[4] = 0.0001; -} - - -// -------------------------------------------------------------------------- -// -// Set the parameter values from the array 'd' -// -// -Bool_t MDispParameters::SetParameters(const TArrayD &d) -{ - if (d.GetSize() != fParameters.GetSize()) - { - *fLog << err << "Sizes of d and of fParameters are different : " - << d.GetSize() << ", " << fParameters.GetSize() << endl; - return kFALSE; - } - - fParameters = d; - - return kTRUE; -} - - -// -------------------------------------------------------------------------- -// -// Set the step sizes from the array 'd' -// -// -Bool_t MDispParameters::SetStepsizes(const TArrayD &d) -{ - if (d.GetSize() != fStepsizes.GetSize()) - { - *fLog << err << "Sizes of d and of fStepsizes are different : " - << d.GetSize() << ", " << fStepsizes.GetSize() << endl; - return kFALSE; - } - - fStepsizes = d; - - return kTRUE; -} Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MDispParameters.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MDispParameters.h (revision 7136) +++ (revision ) @@ -1,40 +1,0 @@ -#ifndef MARS_MDispParameters -#define MARS_MDispParameters - -#ifndef MARS_MParContainer -#include "MParContainer.h" -#endif - -#ifndef ROOT_TArrayD -#include -#endif - -class MDispParameters : public MParContainer -{ -private: - - TArrayD fParameters; // Disp parameters - TArrayD fStepsizes; //! step sizes of Disp parameters - - -public: - - MDispParameters(const char *name=NULL, const char *title=NULL); - - void InitParameters(); - - Bool_t SetParameters(const TArrayD &d); - Bool_t SetStepsizes(const TArrayD &d); - - void SetVariables(const TArrayD &d) { SetParameters(d); } - - TArrayD &GetParameters() { return fParameters; } - TArrayD &GetStepsizes() { return fStepsizes; } - - ClassDef(MDispParameters, 1) // Container for the Disp parameters -}; - -#endif - - - Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MDisplay.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MDisplay.cc (revision 7136) +++ (revision ) @@ -1,176 +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 Rico 05/2004 -! -! Copyright: MAGIC Software Development, 2000-2004 -! -! -\* ======================================================================== */ - -////////////////////////////////////////////////////////////////////////////// -// -// MDisplay -// -// Class to display camera events (MCamEvent) -// You can set an event-by-event display with pause between two consecutive -// events. You can set an output PS file. -// -// Input containers (to be provided to the constructor): -// -// MCamEvent -// MGeomCam -// -// Output containers: -// -// -////////////////////////////////////////////////////////////////////////////// - -#include -#include - -#include "TCanvas.h" -#include "TPostScript.h" - -#include "MParList.h" -#include "MDisplay.h" -#include "MCamEvent.h" -#include "MGeomCam.h" -#include "MHCamera.h" -#include "MRawRunHeader.h" -#include "MRawEvtHeader.h" - -#include "MLog.h" -#include "MLogManip.h" - -ClassImp(MDisplay); - -using namespace std; - -static const TString gsDefName = "MDisplay"; -static const TString gsDefTitle = "Camera display task"; -static const TString gsDefPSFileName = "display.ps"; - -// ------------------------------------------------------------------------- -// -// Constructor. Need to provide the MCamEvent container to be displayed -// and camera geometry . Also the display type can be specified -// (see the MHCamera documentation for more details) -// -MDisplay::MDisplay(MCamEvent* event, MGeomCam* geom, Int_t type, const char* name, const char* title) - : fGeomCam(geom), fCamEvent(event), fCanvas(NULL), fPSFile(NULL), fDisplayType(type), fCreatePSFile(kFALSE), fPause(kTRUE) -{ - fName = name ? name : gsDefName.Data(); - fTitle = title ? title : gsDefTitle.Data(); - - fDisplay = new MHCamera(*geom); - fDisplay->SetPrettyPalette(); - - fPSFileName = gsDefPSFileName; -} -// ------------------------------------------------------------------------- -// -// Destructor -// -MDisplay::~MDisplay() -{ - delete fDisplay; - if(fCanvas) - delete fCanvas; - if(fPSFile) - delete fPSFile; -} - -// ------------------------------------------------------------------------- -// -// Create the canvas, eventually set the batch mode and open ps file -// -Int_t MDisplay::PreProcess(MParList* pList) -{ - fCanvas = new TCanvas("myCanvas","Event Display",600,600); - if(fCreatePSFile) - fPSFile = new TPostScript(fPSFileName,111); - if(!fPause) - fCanvas->SetBatch(); - fCanvas->cd(1); - fDisplay->Draw(); - - //look for the run and event headers - fEvtHeader = (MRawEvtHeader*)pList->FindObject("MRawEvtHeader"); - if (!fEvtHeader) - { - *fLog << err << "MRawEvtHeader not found... aborting." << endl; - return kFALSE; - } - - fRunHeader = (MRawRunHeader*)pList->FindObject(AddSerialNumber("MRawRunHeader")); - if (!fRunHeader) - { - *fLog << err << AddSerialNumber("MRawRunHeader") << " not found... aborting." << endl; - return kFALSE; - } - - return kTRUE; -} - -// ------------------------------------------------------------------------- -// -// Set the new containt of the camera event and update the display. -// Set new page if ps file is requested -// Pause execution if event-by-event display is chosen -// -Int_t MDisplay::Process() -{ - // new page in ps file - if(fPSFile) - fPSFile->NewPage(); - - // update the display contents - fDisplay->SetCamContent(*fCamEvent,fDisplayType); - fCanvas->GetPad(1)->Modified(); - fCanvas->GetPad(1)->Update(); - - *fLog << all << "Run: " << fRunHeader->GetRunNumber() << ", Event: " << fEvtHeader->GetDAQEvtNumber() << endl; - - // pause execution - if(fPause) - { - cout << "Type 'q' to exit, 'p' to print event into ps file, to go on: "; - TString input; - input =cin.get(); - - if (input=='q') - return kFALSE; - if(input=='p') - { - Char_t psfile[200]; - sprintf(psfile,"Run%06dEvent%07d.ps", fRunHeader->GetRunNumber(), fEvtHeader->GetDAQEvtNumber()); - fCanvas->Print(psfile); - } - } - - return kTRUE; -} - -// ------------------------------------------------------------------------- -// -// Close ps file if it was open -// -Int_t MDisplay::PostProcess() -{ - if(fPSFile) fPSFile->Close(); - return kTRUE; -} Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MDisplay.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MDisplay.h (revision 7136) +++ (revision ) @@ -1,57 +1,0 @@ -#ifndef MARS_MDisplay -#define MARS_MDisplay - -#ifndef MARS_MTask -#include "MTask.h" -#endif - -class MRawRunHeader; -class MRawEvtHeader; -class MHCamera; -class MCamEvent; -class MGeomCam; -class TCanvas; -class TPostScript; - -class MDisplay : public MTask -{ - private: - MRawEvtHeader* fEvtHeader; // event header - MRawRunHeader* fRunHeader; // run header - MHCamera* fDisplay; // pointer to the camera display - MGeomCam* fGeomCam; // pointer to the camera geometry - MCamEvent* fCamEvent; // pointer to camera event - TCanvas* fCanvas; // pointer to the canvas - TPostScript* fPSFile; // pointer to ps file - TString fPSFileName; // name for ps file - Int_t fDisplayType; // display type (see MHCamera) - Bool_t fCreatePSFile; // flag to produce a ps file with events - Bool_t fPause; // flag to pause execution between events - - - virtual Int_t PostProcess(); - - protected: - virtual Int_t PreProcess(MParList *plist); - virtual Int_t Process(); - - public: - MDisplay(MCamEvent* event, MGeomCam* geom, Int_t type=0, const char* name=NULL, const char* title=NULL); - virtual ~MDisplay(); - - virtual void Paint(Option_t* option) {}; - - MGeomCam* GetGeomCam() {return fGeomCam;} - Bool_t GetPauseMode() {return fPause;} - Bool_t GetCreatePSFile() {return fCreatePSFile;} - - void SetDisplayType(Int_t type) {fDisplayType=type;} - void SetPSFile(Bool_t set=kTRUE) {fCreatePSFile=set;} - void SetPSFileName(TString name) {fPSFileName=name;} - void SetPause(Bool_t set=kTRUE) {fPause=set;} - - ClassDef(MDisplay, 0) // Task to display camera containers -}; - -#endif - Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MFDisp.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MFDisp.cc (revision 7136) +++ (revision ) @@ -1,299 +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): Eva Domingo, 12/2004 -! Wolfgang Wittek, 12/2004 -! -! -! Copyright: MAGIC Software Development, 2000-2005 -! -! -\* ======================================================================== */ - -////////////////////////////////////////////////////////////////////////////// -// // -// MFDisp // -// // -// This is a class to set cuts to select an event sample // -// for the Disp optimization // -// // -// // -////////////////////////////////////////////////////////////////////////////// - -#include "MFDisp.h" - -#include "MGeomCam.h" -#include "MHillas.h" -#include "MHillasSrc.h" -#include "MImagePar.h" -#include "MNewImagePar.h" - -#include "MParList.h" - -#include "MLog.h" -#include "MLogManip.h" - - -ClassImp(MFDisp); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// Default constructor. -// -MFDisp::MFDisp(const char *name, const char *title) -{ - fName = name ? name : "MFDisp"; - fTitle = title ? title : "Cuts for Disp optimization"; - - // default values of cuts - SetCuts(0, 1000, 0, 1000, 0, 1000, - 0.0, 1.e10, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0); -} - - -// -------------------------------------------------------------------------- -// -// Set the values for the cuts -// -void MFDisp::SetCuts(Int_t islandsmin, Int_t islandsmax, - Int_t usedpixelsmin, Int_t usedpixelsmax, - Int_t corepixelsmin, Int_t corepixelsmax, - Float_t sizemin, Float_t sizemax, - Float_t leakage1min, Float_t leakage1max, - Float_t leakage2min, Float_t leakage2max, - Float_t lengthmin, Float_t widthmin) -{ - fIslandsMin = islandsmin; - fIslandsMax = islandsmax; - - fUsedPixelsMin = usedpixelsmin; - fUsedPixelsMax = usedpixelsmax; - - fCorePixelsMin = corepixelsmin; - fCorePixelsMax = corepixelsmax; - - fSizeMin = sizemin; - fSizeMax = sizemax; - - fLeakage1Min = leakage1min; - fLeakage1Max = leakage1max; - - fLeakage2Min = leakage2min; - fLeakage2Max = leakage2max; - - fLengthMin = lengthmin; - fWidthMin = widthmin; -} - - -// -------------------------------------------------------------------------- -// -Int_t MFDisp::PreProcess(MParList *pList) -{ - MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam"); - if (!cam) - { - *fLog << err << "MGeomCam (Camera Geometry) not found... aborting." << endl; - return kFALSE; - } - - fMm2Deg = cam->GetConvMm2Deg(); - - - fHil = (MHillas*)pList->FindObject("MHillas"); - if (!fHil) - { - *fLog << err << "MHillas not found... aborting." << endl; - return kFALSE; - } - - fHilSrc = (MHillasSrc*)pList->FindObject("MHillasSrc"); - if (!fHilSrc) - { - *fLog << err << "MHillasSrc not found... aborting." << endl; - return kFALSE; - } - - fImgPar = (MImagePar*)pList->FindObject("MImagePar"); - if (!fImgPar) - { - *fLog << err << "MImagePar not found... aborting." << endl; - return kFALSE; - } - - fNewImgPar = (MNewImagePar*)pList->FindObject("MNewImagePar"); - if (!fNewImgPar) - { - *fLog << err << "MNewImagePar not found... aborting." << endl; - return kFALSE; - } - - - memset(fCut, 0, sizeof(fCut)); - - - //-------------------- - *fLog << inf << "MFDisp cut values :" << endl; - - *fLog << inf << " fIslandsMin, fIslandsMax = " - << fIslandsMin << ", " << fIslandsMax << endl; - - *fLog << inf << " fUsedPixelsMin, fUsedPixelsMax = " - << fUsedPixelsMin << ", " << fUsedPixelsMax << endl; - - *fLog << inf << " fCorePixelsMin, fCorePixelsMax = " - << fCorePixelsMin << ", " << fCorePixelsMax << endl; - - *fLog << inf << " fSizeMin, fSizeMax = " - << fSizeMin << ", " << fSizeMax << endl; - - *fLog << inf << " fLeakage1Min, fLeakage1Max = " - << fLeakage1Min << ", " << fLeakage1Max << endl; - - *fLog << inf << " fLeakage2Min, fLeakage2Max = " - << fLeakage2Min << ", " << fLeakage2Max << endl; - - *fLog << inf << " fLengthMin, fWidthMin = " - << fLengthMin << ", " << fWidthMin << endl; - //-------------------- - - return kTRUE; -} - - -// -------------------------------------------------------------------------- -// -// Sum up the number of events passing each cut -// -Bool_t MFDisp::Set(Int_t rc) -{ - fResult = kTRUE; - fCut[rc]++; - return kTRUE; -} - - -// -------------------------------------------------------------------------- -// -// Evaluate cuts for Disp optimization -// if selections are fulfilled set fResult = kTRUE; -// -Int_t MFDisp::Process() -{ - const Double_t length = fHil->GetLength() * fMm2Deg; - const Double_t width = fHil->GetWidth() * fMm2Deg; - //const Double_t dist = fHilSrc->GetDist()* fMm2Deg; - //const Double_t delta = fHil->GetDelta() * kRad2Deg; - const Double_t size = fHil->GetSize(); - - const Double_t leakage1 = fNewImgPar->GetLeakage1(); - const Double_t leakage2 = fNewImgPar->GetLeakage2(); - - const Int_t numusedpixels = fNewImgPar->GetNumUsedPixels(); - const Int_t numcorepixels = fNewImgPar->GetNumCorePixels(); - const Int_t numislands = fImgPar->GetNumIslands(); - - fResult = kFALSE; - - if (numislandsfIslandsMax ) - return Set(1); - - if (numusedpixelsfUsedPixelsMax ) - return Set(2); - - if (numcorepixelsfCorePixelsMax ) - return Set(3); - - if (sizefSizeMax) - return Set(4); - - if (leakage1fLeakage1Max) - return Set(5); - - if (leakage2fLeakage2Max) - return Set(6); - - if (length " << fIslandsMax << endl; - - *fLog << " " << setw(7) << fCut[2] << " (" << setw(3); - *fLog << (int)(fCut[2]*100/GetNumExecutions()+0.5) ; - *fLog << "%) Evts skipped due to: No.of used pixels < " << fUsedPixelsMin ; - *fLog << " or > " << fUsedPixelsMax << endl; - - *fLog << " " << setw(7) << fCut[3] << " (" << setw(3); - *fLog << (int)(fCut[3]*100/GetNumExecutions()+0.5) ; - *fLog << "%) Evts skipped due to: No.of core pixels < " << fCorePixelsMin ; - *fLog << " or > " << fCorePixelsMax << endl; - - *fLog << " " << setw(7) << fCut[4] << " (" << setw(3) ; - *fLog << (int)(fCut[4]*100/GetNumExecutions()+0.5) ; - *fLog << "%) Evts skipped due to: SIZE < " << fSizeMin; - *fLog << " or > " << fSizeMax << endl; - - *fLog << " " << setw(7) << fCut[5] << " (" << setw(3) ; - *fLog << (int)(fCut[5]*100/GetNumExecutions()+0.5) ; - *fLog << "%) Evts skipped due to: Leakage1 < " << fLeakage1Min; - *fLog << " or > " << fLeakage1Max << endl; - - *fLog << " " << setw(7) << fCut[6] << " (" << setw(3) ; - *fLog << (int)(fCut[6]*100/GetNumExecutions()+0.5) ; - *fLog << "%) Evts skipped due to: Leakage2 < " << fLeakage2Min; - *fLog << " or > " << fLeakage2Max << endl; - - *fLog << " " << setw(7) << fCut[7] << " (" << setw(3) ; - *fLog << (int)(fCut[7]*100/GetNumExecutions()+0.5) ; - *fLog << "%) Evts skipped due to: LENGTH <= " << fLengthMin; - *fLog << " or WIDTH <= " << fWidthMin << endl; - - *fLog << " " << fCut[0] << " (" ; - *fLog << (int)(fCut[0]*100/GetNumExecutions()+0.5) ; - *fLog << "%) Evts survived the Disp cuts!" << endl; - *fLog << endl; - - return kTRUE; -} - - - - - Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MFDisp.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MFDisp.h (revision 7136) +++ (revision ) @@ -1,90 +1,0 @@ -#ifndef MARS_MFDisp -#define MARS_MFDisp - -#ifndef MARS_MFilter -#include "MFilter.h" -#endif - -class MHillas; -class MHillasSrc; -class MImagePar; -class MNewImagePar; - -class MFDisp : public MFilter -{ -private: - - MHillas *fHil; - MHillasSrc *fHilSrc; - MImagePar *fImgPar; - MNewImagePar *fNewImgPar; - - //-------------------------- - - // cutting variables to select the sample - // for the Disp optimization - - Int_t fIslandsMin; - Int_t fIslandsMax; - - Int_t fUsedPixelsMin; - Int_t fUsedPixelsMax; - - Int_t fCorePixelsMin; - Int_t fCorePixelsMax; - - Float_t fSizeMin; - Float_t fSizeMax; - - Float_t fLeakage1Min; - Float_t fLeakage1Max; - - Float_t fLeakage2Min; - Float_t fLeakage2Max; - - Float_t fLengthMin; - Float_t fWidthMin; - - //-------------------------- - - Double_t fMm2Deg; // conversion mm to degrees in camera - - Int_t fCut[8]; // array to save cuts statistics - - Bool_t fResult; - - Int_t PreProcess(MParList *pList); - Int_t Process(); - Int_t PostProcess(); - - Bool_t IsExpressionTrue() const { return fResult; } - - Bool_t Set(Int_t rc); - -public: - - MFDisp(const char *name=NULL, const char *title=NULL); - - void SetCuts(Int_t islandsmin, Int_t islandsmax, - Int_t usedpixelsmin, Int_t usedpixelsmax, - Int_t corepixelsmin, Int_t corepixelsmax, - Float_t sizemin, Float_t sizemax, - Float_t leakage1min, Float_t leakage1max, - Float_t leakage2min, Float_t leakage2max, - Float_t lengthmin, Float_t widthmin); - - ClassDef(MFDisp, 0) // Class to set cuts for Disp optimization -}; - -#endif - - - - - - - - - - - Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MFindDisp.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MFindDisp.cc (revision 7136) +++ (revision ) @@ -1,1158 +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): Eva Domingo, 12/2004 -! Wolfgang Wittek, 12/2004 -! -! Copyright: MAGIC Software Development, 2000-2005 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// // -// MFindDisp // -// // -// Class for otimizing the estimation of Disp // -// // -// // -// // -///////////////////////////////////////////////////////////////////////////// -#include "MFindDisp.h" - -#include // fabs - -#include -#include -#include -#include -#include -#include -#include -#include - -#include "MBinning.h" -#include "MContinue.h" -#include "MDispParameters.h" -#include "MDispCalc.h" -#include "MDataElement.h" -#include "MDataMember.h" - -#include "MEvtLoop.h" -#include "MFSelFinal.h" -#include "MF.h" -#include "MFDisp.h" -#include "MFEventSelector.h" -#include "MFEventSelector2.h" -#include "MFillH.h" -#include "MFEventSelector.h" -#include "MGeomCamMagic.h" -#include "MH3.h" -#include "MHDisp.h" -#include "MHFindSignificance.h" -#include "MHMatrix.h" -#include "MHOnSubtraction.h" - -#include "MLog.h" -#include "MLogManip.h" -#include "MMatrixLoop.h" -#include "MMinuitInterface.h" -#include "MParameters.h" -#include "MParList.h" -#include "MProgressBar.h" -#include "MReadMarsFile.h" -#include "MReadTree.h" -#include "MStatusDisplay.h" -#include "MTaskList.h" - - -ClassImp(MFindDisp); - -using namespace std; - - -//------------------------------------------------------------------------ -// fcnDisp -//------------------------------------------------------------------------ -// -// - calculates the quantity to be minimized (by TMinuit) -// -// - the quantity to be minimized is defined in -// MHDisp::Fill() and MHDisp::Finalize() -// -// - the parameters to be varied in the minimization are the parameters -// appearing in the parametrization of Disp (MDispCalc::Calc()) -// -//------------------------------------------------------------------------ - -static void fcnDisp(Int_t &npar, Double_t *gin, Double_t &f, - Double_t *par, Int_t iflag) -{ - cout << "entry fcnDisp" << endl; - - // save pointer to the MINUIT object (for optimizing Disp) for the case - // that it is overwritten by the pointer to another Minuit object - // TMinuit *savePointer = gMinuit; - //------------------------------------------------------------- - - - MEvtLoop *evtloopfcn = (MEvtLoop*)gMinuit->GetObjectFit(); - - MParList *plistfcn = (MParList*)evtloopfcn->GetParList(); - MTaskList *tasklistfcn = (MTaskList*)plistfcn->FindObject("MTaskList"); - - // get needed Disp containers from the existing parList - MDispParameters *dispparams = (MDispParameters*)plistfcn->FindObject("MDispParameters"); - if (!dispparams) - { - gLog << "fcnDisp : MDispParameters object '" << "MDispParameters" - << "' not found... aborting" << endl; - return; - } - - MHDisp *hdisp = (MHDisp*)plistfcn->FindObject("MHDisp"); - if (!hdisp) - { - gLog << "fcnDisp : MHDisp object '" << "MHDisp" - << "' not found... aborting" << endl; - return; - } - - MParameterD *minpar = (MParameterD*)plistfcn->FindObject("MinimizationParameter"); - - // - // transfer current parameter values to MDispParameters - // - // Attention : npar is the number of variable parameters - // not the total number of parameters - // - Double_t fMin, fEdm, fErrdef; - Int_t fNpari, fNparx, fIstat; - gMinuit->mnstat(fMin, fEdm, fErrdef, fNpari, fNparx, fIstat); - - dispparams->SetParameters(TArrayD(fNparx, par)); - - // checking that parameters have been properly set - TArrayD checkparameters = dispparams->GetParameters(); - gLog << "fcnDisp : fNpari, fNparx =" << fNpari << ", " - << fNparx << endl; - gLog << "fcnDisp : i, par, checkparameters =" << endl; - for (Int_t i=0; iEventloop(); - - tasklistfcn->PrintStatistics(0, kTRUE); - - //------------------------------------------- - // get the Minimization parameter value for the current values - // of the Disp parameters - if (minpar) - f = minpar->GetVal(); - else - { - gLog << "fcnDisp : MParameterD object '" << "MinimizationParameter" - << "' not found... aborting" << endl; - return; - } -} - - -// -------------------------------------------------------------------------- -// -// Default constructor. -// -MFindDisp::MFindDisp(MFDisp *fdisp, const char *name, const char *title) - : fHowManyTrain(10000), fHowManyTest(10000) -{ - fName = name ? name : "MFindDisp"; - fTitle = title ? title : "Optimizer of Disp"; - - //--------------------------- - // camera geometry is needed for conversion mm ==> degree - fCam = new MGeomCamMagic; - - // initialize MFDisp filter cuts to default (no cutting) values - // (done at MFDisp constructor) - fDispFilter = fdisp; - - // matrices to contain the training/test samples - fMatrixTrainCalc = new MHMatrix("MatrixTrainCalc"); - fMatrixTrainHists = new MHMatrix("MatrixTrainHists"); - fMatrixTestCalc = new MHMatrix("MatrixTestCalc"); - fMatrixTestHists = new MHMatrix("MatrixTestHists"); - - // objects of MDispCalc where the first part of the matrices mapping is defined - fDispCalcTrain = new MDispCalc("DispTrain","MDispParameters"); - fDispCalcTest = new MDispCalc("DispTest","MDispParametersTest"); - - // objects of MHDisp where the second part of the matrices mapping is defined - fHDispTrain = new MHDisp("DispTrain"); - fHDispTest = new MHDisp("DispTest"); - - // define columns of matrices - // Train matrix - fDispCalcTrain->InitMapping(fMatrixTrainCalc); - fHDispTrain->InitMapping(fMatrixTrainHists); - // Test matrix - fDispCalcTest->InitMapping(fMatrixTestCalc); - fHDispTest->InitMapping(fMatrixTestHists); -} - - -// -------------------------------------------------------------------------- -// -// Default destructor. -// -MFindDisp::~MFindDisp() -{ - delete fCam; - delete fMatrixTrainCalc; - delete fMatrixTrainHists; - delete fMatrixTestCalc; - delete fMatrixTestHists; - delete fDispCalcTrain; - delete fDispCalcTest; - delete fHDispTrain; - delete fHDispTest; -} - - -// -------------------------------------------------------------------------- -// -// Define the matrices 'fMatrixTrainCalc' and 'fMatrixTrainHists' -// for the TRAINING sample -// -// alltogether 'howmanytrain' events are read from file 'nametrain'; -// the events are selected according to a target distribution 'hreftrain' -// -// -Bool_t MFindDisp::DefineTrainMatrix( - const TString &nametrain, MH3 &hreftrain, - const Int_t howmanytrain, const TString &filetrain, - Int_t iflag) -{ - if (nametrain.IsNull() || howmanytrain <= 0) - return kFALSE; - - *fLog << "=============================================" << endl; - *fLog << "Fill TRAINING Matrices from file '" << nametrain << endl; - *fLog << " select " << howmanytrain << " events " << endl; - if (!hreftrain.GetHist().GetEntries()==0) - { - *fLog << " according to a distribution given by the MH3 object '" - << hreftrain.GetName() << "'" << endl; - } - else - { - *fLog << " randomly" << endl; - } - *fLog << "=============================================" << endl; - - - MParList plist; - MTaskList tlist; - - // initialize display to check the reference distribution - // for the event selection (just if iflag==1) - MStatusDisplay *display = NULL; - if (iflag) - display = new MStatusDisplay; - - MReadMarsFile read("Events", nametrain); - read.DisableAutoScheme(); - - // apply cuts for event selection - MContinue contdisp(fDispFilter); - contdisp.SetName("ContFilterSelector2"); - - // choose a reference distribution - MFEventSelector2 seltrain(hreftrain); - seltrain.SetNumMax(howmanytrain); - seltrain.SetName("selectTrain"); - - MFillH filltraincalc(fMatrixTrainCalc); - filltraincalc.SetFilter(&seltrain); - filltraincalc.SetName("fillMatrixTrainCalc"); - - MFillH filltrainhists(fMatrixTrainHists); - filltrainhists.SetFilter(&seltrain); - filltrainhists.SetName("fillMatrixTrainHists"); - - //****************************** - // entries in MParList - - plist.AddToList(&tlist); - plist.AddToList(fCam); - plist.AddToList(fMatrixTrainCalc); - plist.AddToList(fMatrixTrainHists); - - //****************************** - // entries in MTaskList - - tlist.AddToList(&read); - if (fDispFilter != NULL) - tlist.AddToList(&contdisp); - tlist.AddToList(&seltrain); - tlist.AddToList(&filltraincalc); - tlist.AddToList(&filltrainhists); - - //****************************** - - MProgressBar bar; - MEvtLoop evtloop; - if (display != NULL) - evtloop.SetDisplay(display); - evtloop.SetParList(&plist); - evtloop.SetName("EvtLoopMatrixTrain"); - evtloop.SetProgressBar(&bar); - - if (!evtloop.Eventloop()) - return kFALSE; - - tlist.PrintStatistics(0, kTRUE); - - - // print the filled Training Matrices - fMatrixTrainCalc->Print("SizeCols"); - fMatrixTrainHists->Print("SizeCols"); - - // check that number of generated events is compatible with the resquested - Int_t howmanygeneratedcalc = fMatrixTrainCalc->GetM().GetNrows(); - if (TMath::Abs(howmanygeneratedcalc-howmanytrain) > TMath::Sqrt(9.*howmanytrain)) - { - *fLog << "MFindDisp::DefineTrainMatrix; no.of generated events (" - << howmanygeneratedcalc - << ") is incompatible with the no.of requested events (" - << howmanytrain << ")" << endl; - } - - Int_t howmanygeneratedhists = fMatrixTrainHists->GetM().GetNrows(); - if (TMath::Abs(howmanygeneratedhists-howmanytrain) > TMath::Sqrt(9.*howmanytrain)) - { - *fLog << "MFindDisp::DefineTrainMatrix; no.of generated events (" - << howmanygeneratedhists - << ") is incompatible with the no.of requested events (" - << howmanytrain << ")" << endl; - } - - - *fLog << "TRAINING Matrices were filled" << endl; - *fLog << "=============================================" << endl; - - - //-------------------------- - // write out training matrices - - if (filetrain != "") - { - TFile filetr(filetrain, "RECREATE"); - fMatrixTrainCalc->Write(); - fMatrixTrainHists->Write(); - filetr.Close(); - - *fLog << "MFindDisp::DefineTrainMatrix; Training matrices were written onto file '" - << filetrain << "'" << endl; - } - - - if (display != NULL) - delete display; - - return kTRUE; -} - - -// -------------------------------------------------------------------------- -// -// Define the matrices 'fMatrixTestCalc' and 'fMatrixTestHists' -// for the TEST sample -// -// alltogether 'howmanytest' events are read from file 'nametest' -// the events are selected according to a target distribution 'hreftest' -// -// -Bool_t MFindDisp::DefineTestMatrix( - const TString &nametest, MH3 &hreftest, - const Int_t howmanytest, const TString &filetest, - Int_t iflag) -{ - if (nametest.IsNull() || howmanytest<=0) - return kFALSE; - - *fLog << "=============================================" << endl; - *fLog << "Fill TEST Matrices from file '" << nametest << endl; - *fLog << " select " << howmanytest << " events " << endl; - if (!hreftest.GetHist().GetEntries()==0) - { - *fLog << " according to a distribution given by the MH3 object '" - << hreftest.GetName() << "'" << endl; - } - else - { - *fLog << " randomly" << endl; - } - - - MParList plist; - MTaskList tlist; - - // initialize display to check the reference distribution - // for the event selection (just if iflag==1) - MStatusDisplay *display = NULL; - if (iflag) - display = new MStatusDisplay; - - MReadMarsFile read("Events", nametest); - read.DisableAutoScheme(); - - // apply cuts for event selection - MContinue contdisp(fDispFilter); - contdisp.SetName("ContFilterSelector2"); - - // choose a reference distribution - MFEventSelector2 seltest(hreftest); - seltest.SetNumMax(howmanytest); - seltest.SetName("selectTest"); - - MFillH filltestcalc(fMatrixTestCalc); - filltestcalc.SetFilter(&seltest); - filltestcalc.SetName("fillMatrixTestCalc"); - - MFillH filltesthists(fMatrixTestHists); - filltesthists.SetFilter(&seltest); - filltesthists.SetName("fillMatrixTestHists"); - - //****************************** - // entries in MParList - - plist.AddToList(&tlist); - plist.AddToList(fCam); - plist.AddToList(fMatrixTestCalc); - plist.AddToList(fMatrixTestHists); - - //****************************** - // entries in MTaskList - - tlist.AddToList(&read); - if (fDispFilter != NULL) - tlist.AddToList(&contdisp); - tlist.AddToList(&seltest); - tlist.AddToList(&filltestcalc); - tlist.AddToList(&filltesthists); - - //****************************** - - MProgressBar bar; - MEvtLoop evtloop; - if (display != NULL) - evtloop.SetDisplay(display); - evtloop.SetParList(&plist); - evtloop.SetName("EvtLoopMatrixTest"); - evtloop.SetProgressBar(&bar); - - if (!evtloop.Eventloop()) - return kFALSE; - - tlist.PrintStatistics(0, kTRUE); - - - // print the filled Test Matrices - fMatrixTestCalc->Print("SizeCols"); - fMatrixTestHists->Print("SizeCols"); - - // check that number of generated events is compatible with the resquested - const Int_t howmanygeneratedcalc = fMatrixTestCalc->GetM().GetNrows(); - if (TMath::Abs(howmanygeneratedcalc-howmanytest) > TMath::Sqrt(9.*howmanytest)) - { - *fLog << "MFindDisp::DefineTestMatrix; no.of generated events (" - << howmanygeneratedcalc - << ") is incompatible with the no.of requested events (" - << howmanytest << ")" << endl; - } - - const Int_t howmanygeneratedhists = fMatrixTestHists->GetM().GetNrows(); - if (TMath::Abs(howmanygeneratedhists-howmanytest) > TMath::Sqrt(9.*howmanytest)) - { - *fLog << "MFindDisp::DefineTestMatrix; no.of generated events (" - << howmanygeneratedhists - << ") is incompatible with the no.of requested events (" - << howmanytest << ")" << endl; - } - - *fLog << "TEST Matrices were filled" << endl; - *fLog << "=============================================" << endl; - - - //-------------------------- - // write out test matrices - - if (filetest != "") - { - TFile filete(filetest, "RECREATE"); - fMatrixTestCalc->Write(); - fMatrixTestHists->Write(); - filete.Close(); - - *fLog << "MFindDisp::DefineTestMatrix; Test matrices were written onto file '" - << filetest << "'" << endl; - } - - - if (display != NULL) - delete display; - - return kTRUE; -} - - -// -------------------------------------------------------------------------- -// -// Define the matrices for the TRAINING and TEST sample respectively -// -// - this is for the case that training and test samples are generated -// from the same input file (which has the name 'name') -// -Bool_t MFindDisp::DefineTrainTestMatrix( - const TString &name, MH3 &href, - const Int_t howmanytrain, const Int_t howmanytest, - const TString &filetrain, const TString &filetest, - Int_t iflag) -{ - *fLog << "=============================================" << endl; - *fLog << "Fill TRAINING and TEST Matrices from file '" << name << endl; - *fLog << " select " << howmanytrain - << " training and " << howmanytest << " test events " << endl; - if (!href.GetHist().GetEntries()==0) - { - *fLog << " according to a distribution given by the MH3 object '" - << href.GetName() << "'" << endl; - } - else - { - *fLog << " randomly" << endl; - } - - - MParList plist; - MTaskList tlist; - - // initialize display to check the reference distribution - // for the event selection (just if iflag==1) - MStatusDisplay *display = NULL; - if (iflag) - display = new MStatusDisplay; - - MReadMarsFile read("Events", name); - read.DisableAutoScheme(); - - // apply cuts for event selection - MContinue contdisp(fDispFilter); - contdisp.SetName("ContFilterSelector2"); - - // choose a reference distribution - MFEventSelector2 selector(href); - selector.SetNumMax(howmanytrain+howmanytest); - selector.SetName("selectTrainTest"); - selector.SetInverted(); - MContinue cont(&selector); - cont.SetName("ContTrainTest"); - - // choose randomly the events to fill the Training Matrix - Double_t prob = ( (Double_t) howmanytrain ) - / ( (Double_t)(howmanytrain+howmanytest) ); - MFEventSelector split; - split.SetSelectionRatio(prob); - - MFillH filltraincalc(fMatrixTrainCalc); - filltraincalc.SetFilter(&split); - filltraincalc.SetName("fillMatrixTrainCalc"); - - MFillH filltrainhists(fMatrixTrainHists); - filltrainhists.SetFilter(&split); - filltrainhists.SetName("fillMatrixTrainHists"); - - // consider this event as candidate for a test event - // only if event was not accepted as a training event - MContinue conttrain(&split); - conttrain.SetName("ContTrain"); - - MFillH filltestcalc(fMatrixTestCalc); - filltestcalc.SetName("fillMatrixTestCalc"); - - MFillH filltesthists(fMatrixTestHists); - filltesthists.SetName("fillMatrixTestHists"); - - - //****************************** - // entries in MParList - - plist.AddToList(&tlist); - plist.AddToList(fCam); - plist.AddToList(fMatrixTrainCalc); - plist.AddToList(fMatrixTrainHists); - plist.AddToList(fMatrixTestCalc); - plist.AddToList(fMatrixTestHists); - - //****************************** - // entries in MTaskList - - tlist.AddToList(&read); - if (fDispFilter != NULL) - tlist.AddToList(&contdisp); - tlist.AddToList(&cont); - tlist.AddToList(&filltraincalc); - tlist.AddToList(&filltrainhists); - tlist.AddToList(&conttrain); - tlist.AddToList(&filltestcalc); - tlist.AddToList(&filltesthists); - - //****************************** - - MProgressBar bar; - MEvtLoop evtloop; - if (display != NULL) - evtloop.SetDisplay(display); - evtloop.SetParList(&plist); - evtloop.SetName("EvtLoopMatrixTrainTest"); - evtloop.SetProgressBar(&bar); - - Int_t maxev = -1; - if (!evtloop.Eventloop(maxev)) - return kFALSE; - - tlist.PrintStatistics(0, kTRUE); - - - // print the filled Train Matrices - fMatrixTrainCalc->Print("SizeCols"); - fMatrixTrainHists->Print("SizeCols"); - - // check that number of generated events is compatible with the resquested - const Int_t generatedtraincalc = fMatrixTrainCalc->GetM().GetNrows(); - if (TMath::Abs(generatedtraincalc-howmanytrain) > TMath::Sqrt(9.*howmanytrain)) - { - *fLog << "MFindDisp::DefineTrainTestMatrix; no.of generated events (" - << generatedtraincalc - << ") is incompatible with the no.of requested events (" - << howmanytrain << ")" << endl; - } - const Int_t generatedtrainhists = fMatrixTrainHists->GetM().GetNrows(); - if (TMath::Abs(generatedtrainhists-howmanytrain) > TMath::Sqrt(9.*howmanytrain)) - { - *fLog << "MFindDisp::DefineTrainTestMatrix; no.of generated events (" - << generatedtrainhists - << ") is incompatible with the no.of requested events (" - << howmanytrain << ")" << endl; - } - - - // print the filled Test Matrices - fMatrixTestCalc->Print("SizeCols"); - fMatrixTestHists->Print("SizeCols"); - - // check that number of generated events is compatible with the resquested - const Int_t generatedtestcalc = fMatrixTestCalc->GetM().GetNrows(); - if (TMath::Abs(generatedtestcalc-howmanytest) > TMath::Sqrt(9.*howmanytest)) - { - *fLog << "MFindDisp::DefineTrainTestMatrix; no.of generated events (" - << generatedtestcalc - << ") is incompatible with the no.of requested events (" - << howmanytest << ")" << endl; - } - const Int_t generatedtesthists = fMatrixTestHists->GetM().GetNrows(); - if (TMath::Abs(generatedtesthists-howmanytest) > TMath::Sqrt(9.*howmanytest)) - { - *fLog << "MFindDisp::DefineTrainTestMatrix; no.of generated events (" - << generatedtesthists - << ") is incompatible with the no.of requested events (" - << howmanytest << ")" << endl; - } - - - *fLog << "TRAINING and TEST Matrices were filled" << endl; - *fLog << "=============================================" << endl; - - - //---------------------------- - // write out training matrices - - if (filetrain != "") - { - TFile filetr(filetrain, "RECREATE"); - fMatrixTrainCalc->Write(); - fMatrixTrainHists->Write(); - filetr.Close(); - - *fLog << "MFindDisp::DefineTrainTestMatrix; Training matrices were written onto file '" - << filetrain << "'" << endl; - } - - //-------------------------- - // write out test matrices - - if (filetest != "") - { - TFile filete(filetest, "RECREATE"); - fMatrixTestCalc->Write(); - fMatrixTestHists->Write(); - filete.Close(); - - *fLog << "MFindDisp::DefineTrainTestMatrix; Test matrices were written onto file '" - << filetest << "'" << endl; - } - - if (display != NULL) - delete display; - - return kTRUE; -} - - -// -------------------------------------------------------------------------- -// -// Read training and test matrices from files -// -// -Bool_t MFindDisp::ReadMatrix(const TString &filetrain, const TString &filetest) -{ - //-------------------------- - // read in training matrices - - TFile filetr(filetrain); - fMatrixTrainCalc->Read("MatrixTrainCalc"); - fMatrixTrainHists->Read("MatrixTrainHists"); - fMatrixTrainCalc->Print("SizeCols"); - fMatrixTrainHists->Print("SizeCols"); - - *fLog << "MFindDisp::ReadMatrix; Training matrices were read in from file '" - << filetrain << "'" << endl; - filetr.Close(); - - - //-------------------------- - // read in test matrices - - TFile filete(filetest); - fMatrixTestCalc->Read("MatrixTestCalc"); - fMatrixTestHists->Read("MatrixTestHists"); - fMatrixTestCalc->Print("SizeCols"); - fMatrixTestHists->Print("SizeCols"); - - *fLog << "MFindDisp::ReadMatrix; Test matrices were read in from file '" - << filetest << "'" << endl; - filete.Close(); - - return kTRUE; -} - - -//------------------------------------------------------------------------ -// -// Steering program for optimizing Disp -// ------------------------------------ -// -// the criterion for the 'optimum' is defined in -// MHDisp::Fill() and MHDisp::Finalize(); -// for example : smallest sum (over all events) of d^2, where d is the -// distance between the estimated source position -// (calculated using the current value of Disp) and -// the true source position -// -// The various steps are : -// -// - setup the event loop to be executed for each call to fcnDisp -// -// - call TMinuit to do the minimization : -// the fcnDisp function calculates the parameter to minimize -// for the current Disp parameter values; -// for this - Disp is calculated in the event loop by calling -// MDispCalc::Process() ==> MDispCalc::Calc() -// - the Minimization parameter contributions are summed up -// in the event loop by calling MHDisp::Fill() -// - after the event loop the final value of the Minimization -// parameter is calculated by calling MHDisp::Finalize() -// -// Needed as input : (to be set by the Set functions) -// -// - fFilenameParam name of file to which optimum values of the -// parameters are written -// -// - for the minimization, the starting values of the parameters are taken -// - from the file parDispInit (if it is != "") -// - or from the arrays params and/or steps -// - or from the DispParameters constructor -// -//---------------------------------------------------------------------- -Bool_t MFindDisp::FindParams(TString parDispInit, - TArrayD ¶ms, TArrayD &steps) -{ - // Setup the event loop which will be executed in the - // fcnDisp function of MINUIT - // - // parDispInit is the name of the file containing the initial values - // of the parameters; - // if parDispInit = "" 'params' and 'steps' are taken as initial values - // - - *fLog << "=============================================" << endl; - *fLog << "Setup event loop for fcnDisp" << endl; - - - if (fMatrixTrainCalc == NULL || fMatrixTrainHists == NULL) - { - *fLog << "MFindDisp::FindParams; training matrices are not defined... aborting" - << endl; - return kFALSE; - } - - if (fMatrixTrainCalc->GetM().GetNrows() <= 0 || fMatrixTrainHists->GetM().GetNrows() <= 0) - { - *fLog << "MFindDisp::FindParams; training matrices have no entries" - << endl; - return kFALSE; - } - - //--------------------------------------------------------- - MParList parlistfcn; - MTaskList tasklistfcn; - - // loop over rows of matrix - MMatrixLoop loopcalc(fMatrixTrainCalc); - MMatrixLoop loophists(fMatrixTrainHists); - - //-------------------------------- - // create container for the Disp parameters - // and set them to their initial values - MDispParameters *dispparams = fDispCalcTrain->GetDispParameters(); - - // take initial values from file parDispInit - if (parDispInit != "") - { - TFile inparam(parDispInit); - dispparams->Read("MDispParameters"); - inparam.Close(); - *fLog << "MFindDisp::FindParams; initial values of parameters are taken from file " - << parDispInit << endl; - } - - // take initial values from 'params' and/or 'steps' - else if (params.GetSize() != 0 || steps.GetSize() != 0 ) - { - if (params.GetSize() != 0) - { - *fLog << "MFindDisp::FindParams; initial values of parameters are taken from 'params'" - << endl; - dispparams->SetParameters(params); - } - if (steps.GetSize() != 0) - { - *fLog << "MFindDisp::FindParams; initial step sizes are taken from 'steps'" - << endl; - dispparams->SetStepsizes(steps); - } - } - else - { - *fLog << "MFindDisp::FindParams; initial values and step sizes are taken " - << "from the MDispParameters constructor" << endl; - } - - // fill the plots for Disp and sum up the Minimization parameter contributions - MFillH filldispplots("MHDisp", ""); - - //****************************** - // entries in MParList - - parlistfcn.AddToList(&tasklistfcn); - parlistfcn.AddToList(dispparams); - parlistfcn.AddToList(fHDispTrain); - parlistfcn.AddToList(fCam); - parlistfcn.AddToList(fMatrixTrainCalc); - parlistfcn.AddToList(fMatrixTrainHists); - - //****************************** - // entries in MTaskList - - tasklistfcn.AddToList(&loopcalc); - tasklistfcn.AddToList(&loophists); - tasklistfcn.AddToList(fDispCalcTrain); - tasklistfcn.AddToList(&filldispplots); - - //****************************** - - MEvtLoop evtloopfcn("EvtLoopFCN"); - evtloopfcn.SetParList(&parlistfcn); - *fLog << "Event loop for fcnDisp has been setup" << endl; - - // address of evtloopfcn is used in CallMinuit() - - - //----------------------------------------------------------------------- - // - //---------- Start of minimization part -------------------- - // - // Do the minimization with MINUIT - // - // Be careful: This is not thread safe - // - *fLog << "========================================================" << endl; - *fLog << "Start minimization for Disp" << endl; - - - // ------------------------------------------- - // prepare call to MINUIT - // - - // get initial values of parameters - fVinit = dispparams->GetParameters(); - fStep = dispparams->GetStepsizes(); - - TString name[fVinit.GetSize()]; - fStep.Set(fVinit.GetSize()); - fLimlo.Set(fVinit.GetSize()); - fLimup.Set(fVinit.GetSize()); - fFix.Set(fVinit.GetSize()); - - fNpar = fVinit.GetSize(); - - // define names, step sizes, lower and upper limits of the parameters - // fFix[] = 0; means the parameter is not fixed - for (UInt_t i=0; iWrite(); - outparam.Close(); - - *fLog << "Optimum parameter values for Disp were written onto file '" - << fFilenameParam << "' :" << endl; - - const TArrayD &check = dispparams->GetParameters(); - for (Int_t i=0; iDraw(); - - *fLog << "End of Disp optimization part" << endl; - *fLog << "======================================================" << endl; - - return kTRUE; -} - - -// ----------------------------------------------------------------------- -// -// Test the result of the Disp optimization on the test sample -// - -Bool_t MFindDisp::TestParams() -{ - if (fMatrixTestCalc == NULL || fMatrixTestHists == NULL) - { - *fLog << "MFindDisp::TestParams; test matrices are not defined... aborting" - << endl; - return kFALSE; - } - - if (fMatrixTestCalc->GetM().GetNrows() <= 0 || fMatrixTestHists->GetM().GetNrows() <= 0) - { - *fLog << "MFindDisp::TestParams; test matrices have no entries" - << endl; - return kFALSE; - } - - // ------------- TEST the Disp optimization ------------------- - // - *fLog << "Test the Disp optimization on the test sample" << endl; - - // ----------------------------------------------------------------- - // read optimum parameter values from file fFilenameParam - // into array 'dispPar' - - TFile inparam(fFilenameParam); - MDispParameters dispin; - dispin.Read("MDispParameters"); - inparam.Close(); - - *fLog << "Optimum parameter values for Disp were read from file '"; - *fLog << fFilenameParam << "' :" << endl; - - const TArrayD &dispPar = dispin.GetParameters(); - for (Int_t i=0; iGetDispParameters(); - dispparams->SetParameters(dispPar); - - MMatrixLoop loopcalc(fMatrixTestCalc); - MMatrixLoop loophists(fMatrixTestHists); - - MHMatrix *imatrix = NULL; - TArrayI imap; - fHDispTest->GetMatrixMap(imatrix,imap); - - // attention: argument of MHDisp is the name of MImageParDisp container, that should - // be the same than the name given to it when creating MDispCalc object at the MFindDisp - // constructor: fDispCalcTrain = new MDispCalc("DispTest"); - // fill the plots for Disp and sum up the Minimization parameter contributions - MHDisp hdisp1("DispTest"); - hdisp1.SetName("MHDispCorr"); - hdisp1.SetSelectedPos(1); - hdisp1.SetMatrixMap(imatrix,imap); - MFillH filldispplots1("MHDispCorr[MHDisp]", ""); - - MHDisp hdisp2("DispTest"); - hdisp2.SetName("MHDispWrong"); - hdisp2.SetSelectedPos(2); - hdisp2.SetMatrixMap(imatrix,imap); - MFillH filldispplots2("MHDispWrong[MHDisp]", ""); - - MHDisp hdisp3("DispTest"); - hdisp3.SetName("MHDispM3Long"); - hdisp3.SetSelectedPos(3); - hdisp3.SetMatrixMap(imatrix,imap); - MFillH filldispplots3("MHDispM3Long[MHDisp]", ""); - - MHDisp hdisp4("DispTest"); - hdisp4.SetName("MHDispAsym"); - hdisp4.SetSelectedPos(4); - hdisp4.SetMatrixMap(imatrix,imap); - MFillH filldispplots4("MHDispAsym[MHDisp]", ""); - - - //****************************** - // entries in MParList - - parlist2.AddToList(&tasklist2); - parlist2.AddToList(dispparams); - parlist2.AddToList(&hdisp1); - parlist2.AddToList(&hdisp2); - parlist2.AddToList(&hdisp3); - parlist2.AddToList(&hdisp4); - parlist2.AddToList(fCam); - parlist2.AddToList(fMatrixTestCalc); - parlist2.AddToList(fMatrixTestHists); - - //****************************** - // entries in MTaskList - - tasklist2.AddToList(&loopcalc); - tasklist2.AddToList(&loophists); - tasklist2.AddToList(fDispCalcTest); - tasklist2.AddToList(&filldispplots1); - tasklist2.AddToList(&filldispplots2); - tasklist2.AddToList(&filldispplots3); - tasklist2.AddToList(&filldispplots4); - - //****************************** - - MProgressBar bar2; - MEvtLoop evtloop2; - evtloop2.SetParList(&parlist2); - evtloop2.SetName("EvtLoopTestParams"); - evtloop2.SetProgressBar(&bar2); - Int_t maxevents2 = -1; - if (!evtloop2.Eventloop(maxevents2)) - return kFALSE; - - tasklist2.PrintStatistics(0, kTRUE); - - - //------------------------------------------- - // draw the plots - - hdisp1.Draw(); - hdisp2.Draw(); - hdisp3.Draw(); - hdisp4.Draw(); - - //------------------------------------------- - - - *fLog << "Test of Disp otimization finished" << endl; - *fLog << "======================================================" << endl; - - return kTRUE; -} Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MFindDisp.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MFindDisp.h (revision 7136) +++ (revision ) @@ -1,134 +1,0 @@ -#ifndef MARS_MFindDisp -#define MARS_MFindDisp - -#ifndef MARS_MParContainer -#include "MParContainer.h" -#endif - -#ifndef ROOT_TArrayD -#include -#endif - -#ifndef ROOT_TArrayI -#include -#endif - -class MHMatrix; -class MDispCalc; -class MHDisp; -class MFDisp; -class MFilter; -class MH3; -class MGeomCam; -class MEvtLoop; - -class MFindDisp : public MParContainer -{ -private: - - TString fFilenameParam; // filename to store optimum Disp parameters - - TString fFilenameTrain; - TString fFilenameTest; - - Int_t fHowManyTrain; // number of events for training - Int_t fHowManyTest; // number of events for testing - - MHMatrix *fMatrixTrainCalc; // training matrix with variables needed in MDispCalc - MHMatrix *fMatrixTrainHists; // training matrix with variables needed in MHDisp - MHMatrix *fMatrixTestCalc; // testing matrix with variables needed in MDispCalc - MHMatrix *fMatrixTestHists; // testing matrix with variables needed in MHDisp - - MDispCalc *fDispCalcTrain; - MDispCalc *fDispCalcTest; - - MHDisp *fHDispTrain; - MHDisp *fHDispTest; - - MFDisp *fDispFilter; // filter to select an events sample - - MGeomCam *fCam; - - MEvtLoop *fObjectFit; - - //-------------------------------------------- - // To comunicate with MINUIT - //-------------------------------------------- - // attention : dimensions must agree with those in - // MMinuitInterface::CallMinuit() - //char fParName [80][100]; - TArrayD fVinit; - TArrayD fStep; - TArrayD fLimlo; - TArrayD fLimup; - TArrayI fFix; - - UInt_t fNpar; - - TString fMethod; - - Double_t fMin, fEdm, fErrdef; - Int_t fNpari, fNparx, fIstat; - Int_t fErrMinimize; - //-------------------------------------------- - -public: - - MFindDisp(MFDisp *fdisp=NULL, const char *name=NULL, const char *title=NULL); - ~MFindDisp(); - - void SetFilenameParam(const TString &name) {fFilenameParam = name;} - - void SetFilenameTraining(const TString &name, const Int_t howmany) - {fFilenameTrain = name; fHowManyTrain = howmany; } - - void SetFilenameTest(const TString &name, const Int_t howmany) - {fFilenameTest = name; fHowManyTest = howmany; } - - Bool_t DefineTrainMatrix(const TString &name, MH3 &href, - const Int_t howmany, const TString &filetrain, - Int_t iflag=0); - - Bool_t DefineTestMatrix(const TString &name, MH3 &href, - const Int_t howmany, const TString &filetest, - Int_t iflag=0); - - Bool_t DefineTrainTestMatrix(const TString &name, MH3 &href, - const Int_t howmanytrain, const Int_t howmanytest, - const TString &filetrain, const TString &filetest, - Int_t iflag=0); - - MHMatrix *GetMatrixTrainCalc() { return fMatrixTrainCalc; } - MHMatrix *GetMatrixTrainHists() { return fMatrixTrainHists; } - MHMatrix *GetMatrixTestCalc() { return fMatrixTestCalc; } - MHMatrix *GetMatrixTestHists() { return fMatrixTestHists; } - - Bool_t ReadMatrix( const TString &filetrain, const TString &filetest); - - Bool_t FindParams(TString parSCinit, TArrayD ¶ms, TArrayD &steps); - - Bool_t TestParams(); - - ClassDef(MFindDisp, 1) // Class for optimizing Disp -}; - -#endif - - - - - - - - - - - - - - - - - - - Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MHDisp.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MHDisp.cc (revision 7136) +++ (revision ) @@ -1,717 +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): Eva Domingo, 12/2004 -! Wolfgang Wittek, 12/2004 -! -! Copyright: MAGIC Software Development, 2000-2005 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// // -// MHDisp // -// // -// container holding the histograms for Disp // -// also computes the minimization parameter of the Disp optimization // -// // -///////////////////////////////////////////////////////////////////////////// -#include "MHDisp.h" - -#include - -#include -#include -#include -#include -#include -#include -#include - -#include "MGeomCam.h" -#include "MSrcPosCam.h" -#include "MHillas.h" -#include "MHillasExt.h" -#include "MNewImagePar.h" -#include "MMcEvt.hxx" -#include "MPointingPos.h" -#include "MImageParDisp.h" - -#include "MHMatrix.h" -#include "MParameters.h" - -#include "MParList.h" - -#include "MLog.h" -#include "MLogManip.h" - -ClassImp(MHDisp); - -using namespace std; - -enum dispVar_t {kX,kY,kMeanX,kMeanY,kDelta,kSize,kM3Long,kAsym, - kEnergy,kImpact,kLongitmax,kZd,kAz,kTotVar}; // enum variables for the - // matrix columns mapping - -// -------------------------------------------------------------------------- -// -// Constructor -// -MHDisp::MHDisp(const char *imagepardispname, - const char *name, const char *title) - : fImageParDispName(imagepardispname) -{ - fName = name ? name : "MHDisp"; - fTitle = title ? title : "Histograms for Disp"; - - fSelectedPos = 1; // default MHDisp flag (selects Correct Disp source position solution) - - fMatrix = NULL; - - // initialize mapping array dimension to the number of columns of fMatrix - fMap.Set(kTotVar); - - //-------------------------------------------------- - // creating the Disp related histograms - - fHistEnergy = new TH1F("fHistEnergy", - "log10(Energy)", 50, 1., 3.); - fHistEnergy->SetDirectory(NULL); - fHistEnergy->UseCurrentStyle(); - fHistEnergy->SetXTitle("log10(Energy) [GeV]"); - fHistEnergy->SetYTitle("# events"); - - fHistSize = new TH1F("fHistSize", - "log10(Size)", 50, 2., 4.); - fHistSize->SetDirectory(NULL); - fHistSize->UseCurrentStyle(); - fHistSize->SetXTitle("log10(Size) [#phot]"); - fHistSize->SetYTitle("# events"); - - fHistcosZA = new TH1F("fHistcosZA", - "cos(Zenith Angle)", 10, 0., 1.); - fHistcosZA->SetDirectory(NULL); - fHistcosZA->UseCurrentStyle(); - fHistcosZA->SetXTitle("cos(Theta)"); - fHistcosZA->SetYTitle("# events"); - - fSkymapXY = new TH2F("fSkymapXY", - "Disp estimated source positions Skymap", 71, -2., 2., 71, -2., 2.); - fSkymapXY->SetDirectory(NULL); - fSkymapXY->UseCurrentStyle(); - fSkymapXY->SetXTitle("X Disp source position [deg]"); - fSkymapXY->SetYTitle("Y Disp source position [deg]"); - - fHistMinPar = new TH1F("fHistMinPar" , - "Distance^2 between Disp and real srcpos", 100, 0., 2.); - fHistMinPar->SetDirectory(NULL); - fHistMinPar->UseCurrentStyle(); - fHistMinPar->SetXTitle("Minimization parameter = d^2 Disp to real srcpos [deg^2]"); - fHistMinPar->SetYTitle("# events"); - - fHistDuDv = new TH2F("fHistDuDv", - "Du vs. Dv (distances between Disp and real srcpos)", - 100, -2., 2., 100, -2., 2.); - fHistDuDv->SetDirectory(NULL); - fHistDuDv->UseCurrentStyle(); - fHistDuDv->SetXTitle("Dv = transveral distance [deg]"); - fHistDuDv->SetYTitle("Du = longitudinal distance [deg]"); - - fHistMinParEnergy = new TH2F("fHistMinParEnergy", - "Minimization parameter vs. Energy", 50, 1., 3., 100, 0., 2.); - fHistMinParEnergy->SetDirectory(NULL); - fHistMinParEnergy->UseCurrentStyle(); - fHistMinParEnergy->SetXTitle("log10(Energy) [GeV]"); - fHistMinParEnergy->SetYTitle("Minimization parameter = d^2 Disp to real srcpos [deg^2]"); - - fHistMinParSize = new TH2F("fHistMinParSize", - "Minimization parameter vs. Size", 50, 2., 4., 100, 0., 2.); - fHistMinParSize->SetDirectory(NULL); - fHistMinParSize->UseCurrentStyle(); - fHistMinParSize->SetXTitle("log10(Size) [#phot]"); - fHistMinParSize->SetYTitle("Minimization parameter = d^2 Disp to real srcpos [deg^2]"); - - fHistDuEnergy = new TH2F("fHistDuEnergy", - "Du vs. Energy", 50, 1., 3., 100, -2., 2.); - fHistDuEnergy->SetDirectory(NULL); - fHistDuEnergy->UseCurrentStyle(); - fHistDuEnergy->SetXTitle("log10(Energy) [GeV]"); - fHistDuEnergy->SetYTitle("Du = longitudinal distance Disp to real srcpos [deg]"); - - fHistDuSize = new TH2F("fHistDuSize", - "Du vs. Size", 50, 2., 4., 100, -2., 2.); - fHistDuSize->SetDirectory(NULL); - fHistDuSize->UseCurrentStyle(); - fHistDuSize->SetXTitle("log10(Size) [#phot]"); - fHistDuSize->SetYTitle("Du = longitudinal distance Disp to real srcpos [deg]"); - - fHistDvEnergy = new TH2F("fHistDvEnergy", - "Dv vs. Energy", 50, 1., 3., 100, -2., 2.); - fHistDvEnergy->SetDirectory(NULL); - fHistDvEnergy->UseCurrentStyle(); - fHistDvEnergy->SetXTitle("log10(Energy) [GeV]"); - fHistDvEnergy->SetYTitle("Dv = transveral distance Disp to real srcpos [deg]"); - - fHistDvSize = new TH2F("fHistDvSize", - "Dv vs. Size", 50, 2., 4., 100, -2., 2.); - fHistDvSize->SetDirectory(NULL); - fHistDvSize->UseCurrentStyle(); - fHistDvSize->SetXTitle("log10(Size) [#phot]"); - fHistDvSize->SetYTitle("Dv = transveral distance Disp to real srcpos [deg]"); - - fEvCorrAssign = new TProfile("fEvCorrAssign", - "Fraction events srcpos well assign vs. Size", 50, 2., 4., 0., 1.); - fEvCorrAssign->SetDirectory(NULL); - fEvCorrAssign->SetStats(0); - fEvCorrAssign->SetXTitle("log10(Size) [#phot]"); - fEvCorrAssign->SetYTitle("Fraction of events with srcpos well assign"); -} - - -// -------------------------------------------------------------------------- -// -// Destructor -// -MHDisp::~MHDisp() -{ - delete fHistEnergy; - delete fHistSize; - delete fHistcosZA; - delete fSkymapXY; - delete fHistMinPar; - delete fHistDuDv; - delete fHistMinParEnergy; - delete fHistMinParSize; - delete fHistDuEnergy; - delete fHistDuSize; - delete fHistDvEnergy; - delete fHistDvSize; - delete fEvCorrAssign; -} - -// -------------------------------------------------------------------------- -// -// Set the pointers to the containers -// -// -Bool_t MHDisp::SetupFill(const MParList *pList) -{ - // reset all histograms and Minimization parameter computing variables - // before each new eventloop - fNumEv = 0; - fSumMinPar = 0.; - fMinPar = (MParameterD*)const_cast(pList)->FindCreateObj("MParameterD", "MinimizationParameter"); - if (!fMinPar) - { - *fLog << err << "MParameterD (MinimizationParameter) not found and could not be created... aborting." - << endl; - return kFALSE; - } - fMinPar->SetVal(0); - - fHistEnergy->Reset(); - fHistSize->Reset(); - fHistcosZA->Reset(); - fSkymapXY->Reset(); - fHistMinPar->Reset(); - fHistDuDv->Reset(); - fHistMinParEnergy->Reset(); - fHistMinParSize->Reset(); - fHistDuEnergy->Reset(); - fHistDuSize->Reset(); - fHistDvEnergy->Reset(); - fHistDvSize->Reset(); - fEvCorrAssign->Reset(); - - - // look for the defined camera geometry to get mm to deg conversion factor - MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam"); - if (!cam) - { - *fLog << err << "MGeomCam (Camera Geometry) not found... aborting." - << endl; - return kFALSE; - } - fMm2Deg = cam->GetConvMm2Deg(); - - - // look for Disp related containers - fImageParDisp = (MImageParDisp*)pList->FindObject(fImageParDispName, - "MImageParDisp"); - if (!fImageParDisp) - { - *fLog << err << fImageParDispName - << " [MImageParDisp] not found... aborting." << endl; - return kFALSE; - } - - - //----------------------------------------------------------- - // if fMatrix exists, skip preprocessing and just read events from matrix; - // if doesn't exist, read variables values from containers, so look for them - if (fMatrix) - return kTRUE; - - fSrcPos = (MSrcPosCam*)pList->FindObject("MSrcPosCam"); - if (!fSrcPos) - { - *fLog << err << "MSrcPosCam not found... aborting." << endl; - return kFALSE; - } - - fHil = (MHillas*)pList->FindObject("MHillas"); - if (!fHil) - { - *fLog << err << "MHillas not found... aborting." << endl; - return kFALSE; - } - - fHilExt = (MHillasExt*)pList->FindObject("MHillasExt"); - if (!fHilExt) - { - *fLog << err << "MHillasExt not found... aborting." << endl; - return kFALSE; - } - - fNewPar = (MNewImagePar*)pList->FindObject("MNewImagePar"); - if (!fNewPar) - { - *fLog << err << "MNewImagePar not found... aborting." << endl; - return kFALSE; - } - - fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); - if (!fMcEvt) - { - *fLog << err << "MMcEvt not found... This is not a MC file," - << " you are not trying to optimize Disp, just calculating it." - << endl; - // return kFALSE; - } - - fPointing = (MPointingPos*)pList->FindObject("MPointingPos"); - if (!fPointing) - { - *fLog << err << "MPointingPos not found... aborting." << endl; - return kFALSE; - } - - //------------------------------------------ - - return kTRUE; -} - - -// -------------------------------------------------------------------------- -// -// Set which selection algorithm for the Disp estimated source position -// solutions we want to follow when filling the histograms: -// * iflag == 1 --> Correct source position, the closest to the real srcpos -// (only applicable to MC or well known source real data) -// * iflag == 2 --> Wrong source position, the furthest to the real srcpos -// (same applicability than case 1) -// * iflag == 3 --> Source position assigned as M3Long sign indicates -// * iflag == 4 --> Source position assigned as Asym sign indicates -// -void MHDisp::SetSelectedPos(Int_t iflag) -{ - fSelectedPos = iflag; -} - - -// -------------------------------------------------------------------------- -// -// Sets the Matrix and the array of mapped values for each Matrix column -// -void MHDisp::SetMatrixMap(MHMatrix *matrix, TArrayI &map) -{ - fMatrix = matrix; - fMap = map; -} - - -// -------------------------------------------------------------------------- -// -// Returns the Matrix and the mapped value for each Matrix column -// -void MHDisp::GetMatrixMap(MHMatrix* &matrix, TArrayI &map) -{ - map = fMap; - matrix = fMatrix; -} - - -// -------------------------------------------------------------------------- -// -// You can use this function if you want to use a MHMatrix instead of the -// given containers for the Disp optimization. This function adds all -// necessary columns to the given matrix. Afterwards, you should fill -// the matrix with the corresponding data (eg from a file by using -// MHMatrix::Fill). Then, if you loop through the matrix (eg using -// MMatrixLoop), MHDisp::Fill will take the values from the matrix -// instead of the containers. -// -void MHDisp::InitMapping(MHMatrix *mat) -{ - if (fMatrix) - return; - - fMatrix = mat; - - fMap[kX] = fMatrix->AddColumn("MSrcPosCam.fX"); - fMap[kY] = fMatrix->AddColumn("MSrcPosCam.fY"); - fMap[kMeanX] = fMatrix->AddColumn("MHillas.fMeanX"); - fMap[kMeanY] = fMatrix->AddColumn("MHillas.fMeanY"); - fMap[kDelta] = fMatrix->AddColumn("MHillas.fDelta"); - - fMap[kSize] = fMatrix->AddColumn("MHillas.fSize"); - - fMap[kM3Long] = fMatrix->AddColumn("MHillasExt.fM3Long"); - fMap[kAsym] = fMatrix->AddColumn("MHillasExt.fAsym"); - - fMap[kEnergy] = fMatrix->AddColumn("MMcEvt.fEnergy"); - fMap[kImpact] = fMatrix->AddColumn("MMcEvt.fImpact"); - fMap[kLongitmax] = fMatrix->AddColumn("MMcEvt.fLongitmax"); - - fMap[kZd] = fMatrix->AddColumn("MPointingPos.fZd"); - fMap[kAz] = fMatrix->AddColumn("MPointingPos.fAz"); -} - - -// -------------------------------------------------------------------------- -// -// Returns a mapped value from the Matrix -// -Double_t MHDisp::GetVal(Int_t i) const -{ - Double_t val = (*fMatrix)[fMap[i]]; - - //*fLog << "MHDisp::GetVal; i, fMatrix, fMap, val = " - // << i << ", " << fMatrix << ", " << fMap[i] << ", " << val << endl; - - return val; -} - - -// -------------------------------------------------------------------------- -// -// Fill the histograms and sum up the Minimization paramter outcome -// of each processed event -// -Bool_t MHDisp::Fill(const MParContainer *par, const Stat_t w) -{ - Double_t energy = 0.; - Double_t impact = 0.; - Double_t xmax = 0.; - - if ( fMatrix || (!fMatrix && fMcEvt) ) - { - energy = fMatrix ? GetVal(kEnergy) : fMcEvt->GetEnergy(); - impact = fMatrix ? GetVal(kImpact) : fMcEvt->GetImpact(); - xmax = fMatrix ? GetVal(kLongitmax) : fMcEvt->GetLongitmax(); - } - - Double_t theta = fMatrix ? GetVal(kZd) : fPointing->GetZd(); - // Double_t phi = fMatrix ? GetVal(kAz) : fPointing->GetAz(); - - Double_t xsrcpos0 = fMatrix ? GetVal(kX) : fSrcPos->GetX(); - Double_t ysrcpos0 = fMatrix ? GetVal(kY) : fSrcPos->GetY(); - Double_t xmean0 = fMatrix ? GetVal(kMeanX) : fHil->GetMeanX(); - Double_t ymean0 = fMatrix ? GetVal(kMeanY) : fHil->GetMeanY(); - Double_t delta = fMatrix ? GetVal(kDelta) : fHil->GetDelta(); - - Double_t size = fMatrix ? GetVal(kSize) : fHil->GetSize(); - - Double_t m3long = fMatrix ? GetVal(kM3Long) : fHilExt->GetM3Long(); - Double_t asym = fMatrix ? GetVal(kAsym) : fHilExt->GetAsym(); - - //------------------------------------------ - // convert to deg - Double_t xsrcpos = xsrcpos0 * fMm2Deg; - Double_t ysrcpos = ysrcpos0 * fMm2Deg; - Double_t xmean = xmean0 * fMm2Deg; - Double_t ymean = ymean0 * fMm2Deg; - - // calculate cosinus of the angle between d and a vectors - Double_t a = (xmean-xsrcpos)*cos(delta) + (ymean-ysrcpos)*sin(delta); - Double_t b = sqrt( (xmean-xsrcpos)*(xmean-xsrcpos) + (ymean-ysrcpos)*(ymean-ysrcpos) ); - Double_t cosda = a/b; - - // sign of cosda - Int_t s0 = cosda>0 ? 1 : -1; - - // get the sign of M3Long and Asym variables - Int_t sm3 = m3long>0 ? 1 : -1; - Int_t sa = asym>0 ? 1 : -1; - - // set the sign "s" that will select one Disp source position for each - // shower, according to each of the possible algorithms for solution selection: - // SelectedPos = 1 means choose the right source position - // 2 wrong - // 3 the position according to M3Long - // 4 the position according to Asym - Int_t s = s0; - if (fSelectedPos == 1) - s = s0; - else if (fSelectedPos == 2) - s = -s0; - else if (fSelectedPos == 3) - s = sm3; - else if (fSelectedPos == 4) - s = sa; - else - *fLog << "Illegal value for Disp srcpos selection algorithm: " - << " fSelectedPos = " << fSelectedPos << endl; - - // count the number of events where source position has been correctly assigned - if (s == s0) - fEvCorrAssign->Fill(log10(size), 1.); - else - fEvCorrAssign->Fill(log10(size), 0.); - - // get estimated Disp value - Double_t disp = fImageParDisp->GetDisp(); - - //------------------------------------------ - // Disp estimated source position - Double_t xdisp = xmean - s*cos(delta)*disp; - Double_t ydisp = ymean - s*sin(delta)*disp; - fSkymapXY->Fill(xdisp,ydisp); - - // Distance between estimated Disp and real source position - Double_t d2 = (xdisp-xsrcpos)*(xdisp-xsrcpos) + (ydisp-ysrcpos)*(ydisp-ysrcpos); - fHistMinPar->Fill(d2); - - // Longitudinal and transversal distances between Disp and real source positon - Double_t Du = -s*( (xdisp-xsrcpos)*cos(delta) + (ydisp-ysrcpos)*sin(delta)); - Double_t Dv = -s*(-(xdisp-xsrcpos)*sin(delta) + (ydisp-ysrcpos)*cos(delta)); - fHistDuDv->Fill(Dv,Du); - - // Fill Energy, Size and ZA distributions - if (fMatrix || (!fMatrix && fMcEvt)) - fHistEnergy->Fill(log10(energy)); - fHistSize->Fill(log10(size)); - fHistcosZA->Fill(cos(theta/kRad2Deg)); - - // to check the size and energy dependence of the optimization - if (fMatrix || (!fMatrix && fMcEvt)) - { - fHistMinParEnergy->Fill(log10(energy),d2); - fHistDuEnergy->Fill(log10(energy),Du); - fHistDvEnergy->Fill(log10(energy),Dv); - } - fHistMinParSize->Fill(log10(size),d2); - fHistDuSize->Fill(log10(size),Du); - fHistDvSize->Fill(log10(size),Dv); - - // variables for the Minimization parameter calculation (= d^2) - fNumEv += 1; - fSumMinPar += d2; - - return kTRUE; -} - - -// -------------------------------------------------------------------------- -// -// Calculates the final Minimization parameter of the Disp optimization -// -Bool_t MHDisp::Finalize() -{ - fMinPar->SetVal(fSumMinPar/fNumEv); - *fLog << endl; - *fLog << "MHDisp::Finalize: SumMinPar, NumEv = " << fSumMinPar << ", " << fNumEv << endl; - *fLog << "MHDisp::Finalize: Final MinPar = " << fMinPar->GetVal() << endl; - - fMinPar->SetVal(fHistMinPar->GetMean()); - *fLog << "MHDisp::Finalize: Final MinPar = " << fMinPar->GetVal() << endl; - - SetReadyToSave(); - - return kTRUE; -} - - -// -------------------------------------------------------------------------- -// -// Creates a new canvas and draws the Disp related histograms into it. -// Be careful: The histograms belongs to this object and won't get deleted -// together with the canvas. -// -void MHDisp::Draw(Option_t *opt) -{ - gStyle->SetPalette(1); - - TString title = GetName(); - title += ": "; - title += GetTitle(); - - TCanvas *pad = new TCanvas(GetName(),title,0,0,900,1500); - pad->SetBorderMode(0); - pad->Divide(3,5); - - pad->cd(1); - gPad->SetBorderMode(0); - fHistcosZA->SetTitleOffset(1.2,"Y"); - fHistcosZA->DrawClone(opt); - fHistcosZA->SetBit(kCanDelete); - gPad->Modified(); - - pad->cd(2); - gPad->SetBorderMode(0); - fHistEnergy->SetTitleOffset(1.2,"Y"); - fHistEnergy->DrawClone(opt); - fHistEnergy->SetBit(kCanDelete); - gPad->Modified(); - - pad->cd(3); - gPad->SetBorderMode(0); - fHistSize->SetTitleOffset(1.2,"Y"); - fHistSize->DrawClone(opt); - fHistSize->SetBit(kCanDelete); - gPad->Modified(); - - pad->cd(4); - gPad->SetBorderMode(0); - fHistMinPar->SetTitleOffset(1.2,"Y"); - fHistMinPar->DrawClone(opt); - fHistMinPar->SetBit(kCanDelete); - gPad->Modified(); - - TProfile *profMinParEnergy = fHistMinParEnergy->ProfileX("Minimization parameter vs. Energy",0,99999,"s"); - profMinParEnergy->SetXTitle("log10(Energy) [GeV]"); - profMinParEnergy->SetYTitle("Minimization parameter = d^2 Disp to real srcpos [deg^2]"); - pad->cd(5); - gPad->SetBorderMode(0); - profMinParEnergy->SetTitleOffset(1.2,"Y"); - profMinParEnergy->SetStats(0); - profMinParEnergy->DrawClone(opt); - profMinParEnergy->SetBit(kCanDelete); - gPad->Modified(); - - TProfile *profMinParSize = fHistMinParSize->ProfileX("Minimization parameter vs. Size",0,99999,"s"); - profMinParSize->SetXTitle("log10(Size) [#phot]"); - profMinParSize->SetYTitle("Minimization parameter = d^2 Disp to real srcpos [deg^2]"); - pad->cd(6); - gPad->SetBorderMode(0); - profMinParSize->SetTitleOffset(1.2,"Y"); - profMinParSize->SetStats(0); - profMinParSize->DrawClone(opt); - profMinParSize->SetBit(kCanDelete); - gPad->Modified(); - - pad->cd(7); - gPad->SetBorderMode(0); - fSkymapXY->SetTitleOffset(1.2,"Y"); - fSkymapXY->DrawClone("COLZ"); - fSkymapXY->SetBit(kCanDelete); - gPad->Modified(); - - pad->cd(8); - gPad->SetBorderMode(0); - fEvCorrAssign->SetTitleOffset(1.2,"Y"); - fEvCorrAssign->DrawClone(opt); - fEvCorrAssign->SetBit(kCanDelete); - gPad->Modified(); - - pad->cd(9); - gPad->SetBorderMode(0); - fHistDuDv->SetTitleOffset(1.2,"Y"); - fHistDuDv->DrawClone("COLZ"); - fHistDuDv->SetBit(kCanDelete); - gPad->Modified(); - - TH1F *histDu = (TH1F*)fHistDuDv->ProjectionY("histDu"); - histDu->SetTitle("Longitudinal distance Du"); - histDu->SetXTitle("Du = longitudinal distance [deg]"); - pad->cd(10); - gPad->SetBorderMode(0); - histDu->SetTitleOffset(1.2,"Y"); - histDu->DrawClone(opt); - histDu->SetBit(kCanDelete); - gPad->Modified(); - - TProfile *profDuEnergy = fHistDuEnergy->ProfileX("Du vs. Energy",0,99999,"s"); - profDuEnergy->SetXTitle("log10(Energy) [GeV]"); - profDuEnergy->SetYTitle("Du = longitudinal distance [deg]"); - pad->cd(11); - gPad->SetBorderMode(0); - profDuEnergy->SetTitleOffset(1.2,"Y"); - profDuEnergy->SetStats(0); - profDuEnergy->DrawClone(opt); - profDuEnergy->SetBit(kCanDelete); - gPad->Modified(); - - TProfile *profDuSize = fHistDuSize->ProfileX("Du vs. Size",0,99999,"s"); - profDuSize->SetXTitle("log10(Size) [#phot]"); - profDuSize->SetYTitle("Du = longitudinal distance [deg]"); - pad->cd(12); - gPad->SetBorderMode(0); - profDuSize->SetTitleOffset(1.2,"Y"); - profDuSize->SetStats(0); - profDuSize->DrawClone(opt); - profDuSize->SetBit(kCanDelete); - gPad->Modified(); - - TH1F *histDv = (TH1F*)fHistDuDv->ProjectionX("histDv"); - histDv->SetTitle("Transversal distance Dv"); - histDv->SetXTitle("Dv = transversal distance [deg]"); - pad->cd(13); - gPad->SetBorderMode(0); - histDv->SetTitleOffset(1.2,"Y"); - histDv->DrawClone(opt); - histDv->SetBit(kCanDelete); - gPad->Modified(); - - TProfile *profDvEnergy = fHistDvEnergy->ProfileX("Dv vs. Energy",0,99999,"s"); - profDvEnergy->SetXTitle("log10(Energy) [GeV]"); - profDvEnergy->SetYTitle("Dv = transversal distance [deg]"); - pad->cd(14); - gPad->SetBorderMode(0); - profDvEnergy->SetTitleOffset(1.2,"Y"); - profDvEnergy->SetStats(0); - profDvEnergy->DrawClone(opt); - profDvEnergy->SetBit(kCanDelete); - gPad->Modified(); - - TProfile *profDvSize = fHistDvSize->ProfileX("Dv vs. Size",0,99999,"s"); - profDvSize->SetXTitle("log10(Size) [#phot]"); - profDvSize->SetYTitle("Dv = transversal distance [deg]"); - pad->cd(15); - gPad->SetBorderMode(0); - profDvSize->SetTitleOffset(1.2,"Y"); - profDvSize->SetStats(0); - profDvSize->DrawClone(opt); - profDvSize->SetBit(kCanDelete); - gPad->Modified(); -} - - - - - - - - - - - - - - - Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MHDisp.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MHDisp.h (revision 7136) +++ (revision ) @@ -1,111 +1,0 @@ -#ifndef MARS_MHDisp -#define MARS_MHDisp - -#ifndef MARS_MH -#include "MH.h" -#endif - -#ifndef ROOT_TArrayI -#include -#endif - -class TH1F; -class TH2F; -class TProfile; -class MImageParDisp; -class MSrcPosCam; -class MHillas; -class MHillasExt; -class MNewImagePar; -class MMcEvt; -class MPointingPos; -class MHMatrix; -class MParameterD; - -class MHDisp : public MH -{ -private: - - MImageParDisp *fImageParDisp; // container with the estimated Disp - MSrcPosCam *fSrcPos; - MHillas *fHil; - MHillasExt *fHilExt; - MNewImagePar *fNewPar; - MMcEvt *fMcEvt; - MPointingPos *fPointing; - - TString fImageParDispName; - - Int_t fSelectedPos; // flag to select which of the two Disp source position - // solutions we want to fill the histograms (see Set function) - - TH1F *fHistEnergy; // Energy distribution of events - TH1F *fHistSize; // Size distribution of events - TH1F *fHistcosZA; // cosinus Zenith angle distribution of events - TH2F *fSkymapXY; // 2D histogram for Disp estimated source positions - TH1F *fHistMinPar; // Histogram of the event Minimization Parameter (= d^2 = - // = distance^2 between Disp estimated and true source position) - TH2F *fHistDuDv; // Distribution of longitudinal (Du) and transversal - // (Dv) distances between Disp and true source position - TH2F *fHistMinParEnergy; // Minimization Parameter (= d^2) vs. Energy - TH2F *fHistMinParSize; // Minimization Parameter (= d^2) vs. Size - TH2F *fHistDuEnergy; // Du vs. Energy - TH2F *fHistDuSize; // Du vs. Size - TH2F *fHistDvEnergy; // Dv vs. Energy - TH2F *fHistDvSize; // Dv vs. Size - TProfile *fEvCorrAssign; // % events with source position well assign vs. Size - - - Double_t fMm2Deg; // conversion factor from mm to deg - - MHMatrix *fMatrix; // matrix storing variables needed for the Disp optimization - TArrayI fMap; // array storing the matrix mapping column numbers corresponding - // to each variable added to fMatrix - - Double_t GetVal(Int_t i) const; - - Int_t fNumEv; // total number of events - Double_t fSumMinPar; // current sum of the minimization parameter - MParameterD *fMinPar; // final minimization parameters of the Disp optimization - -public: - - MHDisp(const char *imagepardispname = "MImageParDisp", - const char *name=NULL, const char *title=NULL); - ~MHDisp(); - - Bool_t SetupFill(const MParList *plist); - Bool_t Fill(const MParContainer *par, const Stat_t w=1); - Bool_t Finalize(); - - void SetSelectedPos(Int_t iflag); - - void SetMatrixMap(MHMatrix *matrix, TArrayI &map); - void GetMatrixMap(MHMatrix* &matrix, TArrayI &map); // get matrix and its mapping array - - void InitMapping(MHMatrix *mat); // define the matrix of variables - // needed for the Disp optimization - - TH1F *GetHistEnergy() { return fHistEnergy; } - TH1F *GetHistSize() { return fHistSize; } - TH1F *GetHistcosZA() { return fHistcosZA; } - TH2F *GetSkymapXY() { return fSkymapXY; } - TH1F *GetHistMinPar() { return fHistMinPar; } - TH2F *GetHistDuDv() { return fHistDuDv; } - TH2F *GetHistMinParEnergy() { return fHistMinParEnergy; } - TH2F *GetHistMinParSize() { return fHistMinParSize; } - TH2F *GetHistDuEnergy() { return fHistDuEnergy; } - TH2F *GetHistDuSize() { return fHistDuSize; } - TH2F *GetHistDvEnergy() { return fHistDvEnergy; } - TH2F *GetHistDvSize() { return fHistDvSize; } - TProfile *GetEvCorrAssign() { return fEvCorrAssign; } - - void Draw(Option_t *opt=""); - - ClassDef(MHDisp, 1) // Container holding the Histograms for Disp and - // the parameter to minimize in the Disp optimization -}; - -#endif - - Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MImageParDisp.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MImageParDisp.cc (revision 7136) +++ (revision ) @@ -1,69 +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): Eva Domingo , 12/2004 -! Wolfgang Wittek, 12/2004 -! -! Copyright: MAGIC Software Development, 2000-2005 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// -// MImageParDisp -// -// Storage Container for estimated distance to source position (Disp) -// -///////////////////////////////////////////////////////////////////////////// -#include "MImageParDisp.h" - -#include "MLog.h" -#include "MLogManip.h" - - -ClassImp(MImageParDisp); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// Default constructor. -// -MImageParDisp::MImageParDisp(const char *name, const char *title) -{ - fName = name ? name : "MImageParDisp"; - fTitle = title ? title : "Estimated distance to source position (Disp)"; - - Reset(); -} - -// -------------------------------------------------------------------------- -// -void MImageParDisp::Reset() -{ - fDisp = -1; -} - -// -------------------------------------------------------------------------- -// -void MImageParDisp::Print(Option_t *) const -{ - *fLog << all; - *fLog << "Estimated distance to source position (" << GetName() << ")" << endl; - *fLog << " - Disp [deg] = " << fDisp << endl; -} Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MImageParDisp.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MImageParDisp.h (revision 7136) +++ (revision ) @@ -1,30 +1,0 @@ -#ifndef MARS_MImageParDisp -#define MARS_MImageParDisp - -#ifndef MARS_MParContainer -#include "MParContainer.h" -#endif - -class MImageParDisp : public MParContainer -{ -private: - - Float_t fDisp; // [deg] distance between the image center - // and the estimated source position, - // along the major axis of the image - -public: - - MImageParDisp(const char *name=NULL, const char *title=NULL); - - void Reset(); - - void SetDisp(Float_t disp) { fDisp = disp; } - Float_t GetDisp() const { return fDisp; } - - void Print(Option_t *opt=NULL) const; - - ClassDef(MImageParDisp, 1) // Container to hold estimated distance to source position (Disp) -}; - -#endif Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MLiveTime.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MLiveTime.cc (revision 7136) +++ (revision ) @@ -1,75 +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 , 7/2004 -! -! Copyright: MAGIC Software Development, 2000-2004 -! -! -\* ======================================================================== */ - -#include "MLiveTime.h" - -#include "MLog.h" -#include "MLogManip.h" - -ClassImp(MLiveTime); - -using namespace std; - -MLiveTime::MLiveTime(UInt_t numberbins, const char *name, const char *title) -{ - - fName = name ? name : "MLiveTime"; - fTitle = title ? title : "Container to hold the live time for a certain time bin"; - - Set(numberbins); -} - -void MLiveTime::Set(UInt_t numberbins) -{ - fNumberTimeBins = numberbins; - fLiveTimeBin.Set(numberbins); - fMeanRealTimeBin.Set(numberbins); - fWidthRealTimeBin.Set(numberbins); -} - -void MLiveTime::Print(const Option_t* o) const -{ - TString opt = o; - - Double_t totalLiveTime = 0; - - if (opt.Contains("last", TString::kIgnoreCase)) - { - *fLog << GetName() << ": Present real time bin " << setprecision(10) << fMeanRealTimeBin[fNumberTimeBins-1] << " +- " << setprecision(5) << fWidthRealTimeBin[fNumberTimeBins-1] << " [" << 2*fWidthRealTimeBin[fNumberTimeBins-1]*24*60*60 << " sec] MJD" << endl; - *fLog << GetName() << ": Present live time " << fLiveTimeBin[fNumberTimeBins-1] << " sec" << endl; - } - else if (opt.Contains("all", TString::kIgnoreCase)) - { - *fLog << GetName() << ": " << fNumberTimeBins << " time bins" << endl; - for (UInt_t bin = 0; bin -#endif - -class MLiveTime : public MParContainer -{ - private: - - UInt_t fNumberTimeBins; - TArrayD fLiveTimeBin; - TArrayD fMeanRealTimeBin; - TArrayD fWidthRealTimeBin; - - void Set(UInt_t numberbins); - - public: - - MLiveTime(UInt_t numberbins = 1, const char *name=NULL, const char *title=NULL); - //~MLiveTime; - - void AddBin() { Set(fNumberTimeBins+1); } - void AddTime(Double_t time) { fLiveTimeBin[fNumberTimeBins-1]+=time; } - void SetRealTime (Double_t mean, Double_t width) - { - fMeanRealTimeBin[fNumberTimeBins-1]=mean; - fWidthRealTimeBin[fNumberTimeBins-1]=width; - } - - UInt_t GetNumberTimeBins() { return fNumberTimeBins;} - Double_t GetLiveTime() { return fLiveTimeBin[fNumberTimeBins-1]; } - Double_t GetMeanRealTime() { return fMeanRealTimeBin[fNumberTimeBins-1]; } - Double_t GetWidthRealTime() { return fWidthRealTimeBin[fNumberTimeBins-1]; } - - TArrayD& GetLiveTimeTArray() { return fLiveTimeBin; } - TArrayD& GetMeanRealTimeTArray() { return fMeanRealTimeBin; } - TArrayD& GetWidthRealTimeTArray() { return fWidthRealTimeBin; } - - void Print(const Option_t*) const; - - ClassDef(MLiveTime, 1) // Storage for the live time extracted from real data - -}; - -#endif Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MLiveTimeCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MLiveTimeCalc.cc (revision 7136) +++ (revision ) @@ -1,360 +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 , 7/2004 -! -! Copyright: MAGIC Software Development, 2000-2004 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// -// MLiveTimeCalc -// -///////////////////////////////////////////////////////////////////////////// -#include "MLiveTimeCalc.h" - -#include "MRawRunHeader.h" -#include "MRawEvtHeader.h" -#include "MLiveTime.h" - -#include "MLog.h" -#include "MLogManip.h" - -#include "MParList.h" -#include "MTaskList.h" - -ClassImp(MLiveTimeCalc) -using namespace std; - -Bool_t Debug = kFALSE; -Bool_t PrintNewRun = kTRUE; - -MLiveTimeCalc::MLiveTimeCalc(const char *name, const char *title) : kSecTomSec(1e3), kDayToSec(24.*60.*60.), fRunHeader(NULL), fEvtHeader(NULL), fPresentEventTime(NULL), fLiveTime(NULL) -{ - - fName = name ? name : "MLiveTimeCalc"; - fTitle = title ? title : "Task to compute the live time from real data."; - - for (Int_t i=0; i<10; i++) - fNumberLostEvents[i] = 0; - - fFirstEventMjd = 0.; - fPresentEventMjd = 0.; - fLastEventMjd = 0.; - - fPresentEventNumber = 0; - fLastEventNumber = 0; - - fPresentEventRunNumber = 0; - fLastEventRunNumber = 0; - - fRealTimeBinSize = 0.; - - fPrintNextEvent = kFALSE; -} - -Int_t MLiveTimeCalc::PreProcess(MParList *pList) -{ - - fEvtHeader = (MRawEvtHeader*)pList->FindObject(AddSerialNumber("MRawEvtHeader")); - if (!fEvtHeader) - { - *fLog << err << AddSerialNumber("MRawEvtHeader") << " not found ... aborting" << endl; - return kFALSE; - } - - fRunHeader = (MRawRunHeader*)pList->FindObject(AddSerialNumber("MRawRunHeader")); - if (!fRunHeader) - { - *fLog << err << AddSerialNumber("MRawRunHeader") << " not found ... aborting" << endl; - return kFALSE; - } - - fPresentEventTime = (MTime*)pList->FindObject(AddSerialNumber("MTime")); - if (!fPresentEventTime) - { - *fLog << err << AddSerialNumber("MTime") << " not found ... aborting" << endl; - return kFALSE; - } - - - fLiveTime = (MLiveTime*)pList->FindCreateObj(AddSerialNumber("MLiveTime")); - if (!fLiveTime) - { - *fLog << err << AddSerialNumber("MLiveTime") << " cannot be created ... aborting" << endl; - return kFALSE; - } - - return kTRUE; -} - -Int_t MLiveTimeCalc::Process() -{ - - fPresentEventMjd = fPresentEventTime->GetMjd(); - fLastEventMjd = fLastEventTime.GetMjd(); - fPreaviusLastEventMjd = fPreaviusLastEventTime.GetMjd(); - - fPresentEventNumber = fEvtHeader->GetDAQEvtNumber(); - fPresentEventRunNumber = fRunHeader->GetRunNumber();; - - if (fPrintNextEvent && Debug) - { - *fLog << dbg << GetName() << ": Printing next event" << endl; - Print("all"); - fLiveTime->Print("last"); - *fLog << inf << endl; - fPrintNextEvent = kFALSE; - } - - if (fFirstEventMjd == 0) - { - fFirstEventMjd = fPresentEventMjd; - - fPreaviusLastEventTime = fLastEventTime; - fLastEventTime = *fPresentEventTime; - - fPreaviusLastEventRunNumber = fLastEventRunNumber; - fLastEventRunNumber = fPresentEventRunNumber; - - fPreaviusLastEventNumber = fLastEventNumber; - fLastEventNumber = fPresentEventNumber; - } - else - { - if (isTimeStampOk()) - { - if (fRealTimeBinSize > 0 && - (fPresentEventMjd - fFirstEventMjd)*kDayToSec > fRealTimeBinSize) - { - - Double_t width = (fLastEventMjd - fFirstEventMjd)/2; - Double_t mean = fFirstEventMjd + width; - - fLiveTime->SetRealTime(mean,width); - - if (Debug) - { - *fLog << inf << GetName() << ": New time bin" << endl; - Print("all"); - *fLog << GetName() << ": First event time " << setprecision(10) << fFirstEventMjd << setprecision(5) << endl; - fLiveTime->Print("last"); - } - - fLiveTime->AddBin(); - fFirstEventMjd = fPresentEventMjd; - } - - if (fPresentEventRunNumber!=fLastEventRunNumber) - { - if (fLastEventRunNumber != 0 && PrintNewRun) - { - *fLog << dbg << GetName() << ": New run" << endl; - Print("all"); - fLiveTime->Print("last"); - *fLog << inf << endl; - } - fLastEventRunNumber = fPresentEventRunNumber; - } - else - fLiveTime->AddTime((fPresentEventMjd-fLastEventMjd)*kDayToSec); - - - - fPreaviusLastEventTime = fLastEventTime; - fLastEventTime = *fPresentEventTime; - - fPreaviusLastEventRunNumber = fLastEventRunNumber; - fLastEventRunNumber = fPresentEventRunNumber; - - fPreaviusLastEventNumber = fLastEventNumber; - fLastEventNumber = fPresentEventNumber; - - } - } - - return kTRUE; -} - -Int_t MLiveTimeCalc::PostProcess() -{ - Double_t width = (fLastEventMjd - fFirstEventMjd)/2; - Double_t mean = fFirstEventMjd + width; - - fLiveTime->SetRealTime(mean,width); - - *fLog << dbg << endl; - *fLog << dbg << GetName() << ": PostProcess" << endl; - fLiveTime->Print("all"); - *fLog << inf << endl; - - *fLog << GetName() << " execution statistics:" << endl; - *fLog << dec << setfill(' '); - - ULong_t fTotalNumberLostEvents = 0; - for (Int_t i=0; i<6; i++) - fTotalNumberLostEvents += fNumberLostEvents[i]; - - *fLog << " " << setw(7) << fTotalNumberLostEvents << " (" << setw(3) ; - *fLog << (Int_t)(fTotalNumberLostEvents*100/GetNumExecutions()) ; - *fLog << "%) bad time stamp events" << endl; - - *fLog << "\t\t(" << setw(3) << (Int_t)(fNumberLostEvents[0]*100/fTotalNumberLostEvents) ; - *fLog << "%) time stamp == 0" << endl; - - *fLog << "\t\t(" << setw(3) << (Int_t)(fNumberLostEvents[1]*100/fTotalNumberLostEvents) ; - *fLog << "%) last time stamp == 0" << endl; - - *fLog << "\t\t(" << setw(3) << (Int_t)(fNumberLostEvents[2]*100/fTotalNumberLostEvents) ; - *fLog << "%) time stamp in the past" << endl; - - *fLog << "\t\t(" << setw(3) << (Int_t)(fNumberLostEvents[3]*100/fTotalNumberLostEvents) ; - *fLog << "%) time stamp == last one" << endl; - - *fLog << "\t\t(" << setw(3) << (Int_t)(fNumberLostEvents[4]*100/fTotalNumberLostEvents) ; - *fLog << "%) time stamp just with integer part" << endl; - - *fLog << "\t\t(" << setw(3) << (Int_t)(fNumberLostEvents[5]*100/fTotalNumberLostEvents) ; - *fLog << "%) run number < last one" << endl; - - return kTRUE; -} - -Bool_t MLiveTimeCalc::isTimeStampOk() -{ - - Bool_t result = kTRUE; - - - if (fPresentEventMjd == 0) - { - - if (Debug) - { - *fLog << err << GetName() << ": Present event time stamp equal to 0" << endl; - Print("all"); - } - - fNumberLostEvents[0]++; - result = kFALSE; - fPrintNextEvent = kTRUE; - } - else if (fLastEventMjd == 0) - { - - if (Debug) - { - *fLog << err << GetName() << ": Last event time stamp equal to 0" << endl; - Print("all"); - } - - fNumberLostEvents[1]++; - result = kFALSE; - fPrintNextEvent = kTRUE; - } - else if (fPresentEventMjd-fLastEventMjd < 0) - { - - if (Debug) - { - *fLog << err << GetName() << ": Present event with time stamp in the past" << endl; - Print("all"); - } - - fNumberLostEvents[2]++; - result = kFALSE; - fPrintNextEvent = kTRUE; - } - else if (fPresentEventMjd-fLastEventMjd == 0) - { - - if (Debug) - { - *fLog << err << GetName() << ": Present event time stamp equal than last event" << endl; - Print("all"); - } - - fNumberLostEvents[3]++; - result = kFALSE; - fPrintNextEvent = kTRUE; - } - else if (fPresentEventNumber- fLastEventNumber<= 0 && fPresentEventRunNumber == fLastEventRunNumber) - { - - if (Debug) - { - *fLog << warn << GetName() << ": Present event number equal or smaller than last event" << endl; - Print("all"); - } - - result = kTRUE; - fPrintNextEvent = kTRUE; - } - else if ((Int_t)fPresentEventMjd == fPresentEventMjd) - { - - if (Debug) - { - *fLog << err << GetName() << ": Present event time stamp idetical to midnight" << endl; - Print("all"); - } - - fNumberLostEvents[4]++; - result = kFALSE; - fPrintNextEvent = kTRUE; - } - else if ((fPresentEventMjd-fLastEventMjd)*kDayToSec > 1.) - { - - if (Debug) - { - *fLog << warn << GetName() << ": Time from last event bigger than 1 sec" << endl; - Print("all"); - } - - result = kTRUE; - fPrintNextEvent = kTRUE; - } - else if (fPresentEventRunNumber-fLastEventRunNumber < 0) - { - - if (Debug) - { - *fLog << warn << GetName() << ": Present run number previuos than last one" << endl; - Print("all"); - } - - fNumberLostEvents[5]++; - result = kTRUE; - fPrintNextEvent = kTRUE; - } - - return result; -} - -void MLiveTimeCalc::Print(const Option_t *o) const -{ - - *fLog << "Present event run number [" << fPresentEventRunNumber << "] event number [" << fPresentEventNumber << ']' << endl; - *fLog << "Last event run number [" << fLastEventRunNumber << "] event number [" << fLastEventNumber << ']' << endl; - *fLog << "Preavius last event run number [" << fPreaviusLastEventRunNumber << "] event number [" << fPreaviusLastEventNumber << ']' << endl; - *fLog << "Present, last and preavius to last event times:" << endl; - fPresentEventTime->Print(); - fLastEventTime.Print(); - fPreaviusLastEventTime.Print(); -} Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MLiveTimeCalc.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MLiveTimeCalc.h (revision 7136) +++ (revision ) @@ -1,71 +1,0 @@ -#ifndef MARS_MLiveTimeCalc -#define MARS_MLiveTimeCalc - -#ifndef MARS_MTask -#include "MTask.h" -#endif - -#ifndef MARS_MTime -#include "MTime.h" -#endif - -class MRawRunHeader; -class MRawEvtHeader; -class MLiveTime; - -class MLiveTimeCalc : public MTask -{ - public: - - const Double_t kSecTomSec; - const Double_t kDayToSec; - - private: - - MRawRunHeader* fRunHeader; - MRawEvtHeader* fEvtHeader; - - MTime* fPresentEventTime; - MTime fLastEventTime; - MTime fPreaviusLastEventTime; - - MLiveTime* fLiveTime; - - - UInt_t fNumberLostEvents[10]; - - Double_t fFirstEventMjd; - Double_t fPresentEventMjd; - Double_t fLastEventMjd; - Double_t fPreaviusLastEventMjd; - - Int_t fPresentEventNumber; - Int_t fLastEventNumber; - Double_t fPreaviusLastEventNumber; - - Int_t fPresentEventRunNumber; - Int_t fLastEventRunNumber; - Double_t fPreaviusLastEventRunNumber; - - Double_t fRealTimeBinSize; - - Bool_t fPrintNextEvent; - - public: - - MLiveTimeCalc(const char *name=NULL, const char *title=NULL); - //~MLiveTimeCalc; - - Int_t PreProcess(MParList *pList); - Int_t Process(); - Int_t PostProcess(); - - void SetRealTimeBinSize(Double_t time) { fRealTimeBinSize=time; } - - Bool_t isTimeStampOk(); - void Print(const Option_t *o) const; - - ClassDef(MLiveTimeCalc, 0) // Task to compute the live time from real data -}; - -#endif Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MTopologyCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MTopologyCalc.cc (revision 7136) +++ (revision ) @@ -1,133 +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): Josep Flix 09/2004 -! -! Copyright: MAGIC Software Development, 2004 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// -// MTopologyCalc -// -// Taks that evaluates the Topological Parameters of images after image cleaning -// -// fDistance = SUM_(i,j) D_ij , where i,jFindObject("MCerPhotEvt"); - if (!fCerPhotEvt) - { - *fLog << err << dbginf << "MCerPhotEvt not found... aborting." << endl; - return kFALSE; - } - - fGeomCam = (MGeomCam*)pList->FindObject("MGeomCam"); - if (!fGeomCam) - { - *fLog << err << dbginf << "MGeomCam not found... aborting." << endl; - return kFALSE; - } - - fTopology = (MTopology*)pList->FindCreateObj("MTopology"); - if (!fTopology) - return kFALSE; - - return kTRUE; - -} - -// ------------------------------------------------------------------------- -// -Int_t MTopologyCalc::Process() -{ - - Int_t rc = fTopology->Calc(*fGeomCam, *fCerPhotEvt); - - if (rc == 1) - fErrors++; - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Prints some statistics about the Topology calculation. The percentage -// is calculated with respect to the number of executions of this task. -// -Int_t MTopologyCalc::PostProcess() -{ - if (GetNumExecutions()==0) - return kTRUE; - - *fLog << inf << endl; - *fLog << GetDescriptor() << " execution statistics:" << endl; - *fLog << dec << setfill(' '); - *fLog << " " << fErrors << " (" << (int)(fErrors*100/GetNumExecutions()) - << "%) Evts skipped due to UsedPixels == -1" << endl; - *fLog << endl; - - /* - delete fHillas; - delete fHillasSrc; - delete fSrcPos; - delete fTopology; - */ - - return kTRUE; -} Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MTopologyCalc.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MTopologyCalc.h (revision 7136) +++ (revision ) @@ -1,34 +1,0 @@ -#ifndef MARS_MTopologyCalc -#define MARS_MTopologyCalc - -#ifndef MARS_MTask -#include "MTask.h" -#endif - -class MGeomCam; -class MCerPhotEvt; -class MTopology; - -class MTopologyCalc : public MTask -{ -private: - - const MGeomCam *fGeomCam; - const MCerPhotEvt *fCerPhotEvt; - - MTopology *fTopology; - - Int_t fErrors; - - Int_t PreProcess(MParList *plist); - Int_t Process(); - Int_t PostProcess(); - -public: - - MTopologyCalc(const char *name=NULL, const char *title=NULL); - - ClassDef(MTopologyCalc, 0) // task to calculate.... -}; - -#endif