Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MDispCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MDispCalc.cc (revision 5671) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MDispCalc.cc (revision 5671) @@ -0,0 +1,294 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// // +// MDispCalc // +// // +// This task calculates Disp // +// (parameters are taken from the container MDisp) // +// Also the matrix of variables to be used in the Disp optimization // +// is defined // +// // +// // +///////////////////////////////////////////////////////////////////////////// +#include "MDispCalc.h" + +#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 "MDisp.h" + +#include "MHMatrix.h" + +#include "MParList.h" + +#include "MLog.h" +#include "MLogManip.h" + + +ClassImp(MDispCalc); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// constructor +// +MDispCalc::MDispCalc(const char *imagepardispname, const char *dispname, + const char *name, const char *title) + : fImageParDispName(imagepardispname), fDispName(dispname) +{ + fName = name ? name : "MDispCalc"; + fTitle = title ? title : "Class to calculate Disp"; + + fMatrix = NULL; + + // initialize mapping array dimension to the number of columns of fMatrix + fMap.Set(18); +} + + +// -------------------------------------------------------------------------- +// +Int_t MDispCalc::PreProcess(MParList *pList) +{ + // 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->FindCreateObj("MImageParDisp", + fImageParDispName); + if (!fImageParDisp) + { + *fLog << err << fImageParDispName + << " [MImageParDisp] not found... aborting." << endl; + return kFALSE; + } + + fDisp = (MDisp*)pList->FindObject(fDispName, "MDisp"); + if (!fDisp) + { + *fLog << err << fDispName << " [MDisp] 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... " + << "MMcEvt is going to be used to get Theta and Phi." + << endl; + // return kFALSE; + } + + //------------------------------------------ + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// 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), MDispCalc::Process will take the values from the matrix +// instead of the containers. +// +void MDispCalc::InitMapping(MHMatrix *mat) +{ + if (fMatrix) + return; + + fMatrix = mat; + + fMap[0] = fMatrix->AddColumn("MSrcPosCam.fX"); + fMap[1] = fMatrix->AddColumn("MSrcPosCam.fY"); + fMap[2] = fMatrix->AddColumn("MHillas.fMeanX"); + fMap[3] = fMatrix->AddColumn("MHillas.fMeanY"); + fMap[4] = fMatrix->AddColumn("MHillas.fDelta"); + + fMap[5] = fMatrix->AddColumn("MHillas.fSize"); + fMap[6] = fMatrix->AddColumn("MHillas.fWidth"); + fMap[7] = fMatrix->AddColumn("MHillas.fLength"); + + fMap[8] = fMatrix->AddColumn("MNewImagePar.fConc"); + fMap[9] = fMatrix->AddColumn("MNewImagePar.fLeakage1"); + fMap[10] = fMatrix->AddColumn("MNewImagePar.fLeakage2"); + + fMap[11] = fMatrix->AddColumn("MHillasExt.fM3Long"); + fMap[12] = fMatrix->AddColumn("MHillasExt.fAsym"); + + fMap[13] = fMatrix->AddColumn("MMcEvt.fEnergy"); + fMap[14] = fMatrix->AddColumn("MMcEvt.fImpact"); + fMap[15] = fMatrix->AddColumn("MMcEvt.fLongitmax"); + + if (fPointing) + { + fMap[16] = fMatrix->AddColumn("MPointingPos.fZd"); + fMap[17] = fMatrix->AddColumn("MPointingPos.fAz"); + } + else + { + fMap[16] = fMatrix->AddColumn("MMcEvt.fTelescopeTheta*kRad2Deg"); + fMap[17] = fMatrix->AddColumn("MMcEvt.fTelescopePhi*kRad2Deg"); + } +} + + +// -------------------------------------------------------------------------- +// +// Returns the Matrix and the mapped value for each Matrix column +// +MHMatrix* MDispCalc::GetMatrixMap(TArrayI &map) +{ + map = fMap; + + return fMatrix; +} + + +// -------------------------------------------------------------------------- +// +// Returns a mapped value from the Matrix +// +Double_t MDispCalc::GetVal(Int_t i) const +{ + Double_t val = (*fMatrix)[fMap[i]]; + + // *fLog << "MDispCalc::GetVal; i, fMatrix, fMap, val = " + // << i << ", " << fMatrix << ", " << fMap[i] << ", " << val << endl; + + return val; +} + + +// --------------------------------------------------------------------------- +// +// Calculate Disp +// +// +Int_t MDispCalc::Process() +{ + // get variables needed to compute disp from the matrix or the containers + const Double_t size = fMatrix ? GetVal(5) : fHil->GetSize(); + const Double_t width0 = fMatrix ? GetVal(6) : fHil->GetWidth(); + const Double_t length0 = fMatrix ? GetVal(7) : fHil->GetLength(); + const Double_t conc = fMatrix ? GetVal(8) : fNewPar->GetConc(); + const Double_t leakage1 = fMatrix ? GetVal(9) : fNewPar->GetLeakage1(); + const Double_t leakage2 = fMatrix ? GetVal(10) : fNewPar->GetLeakage2(); + + // convert to deg + const Double_t width = width0 * fMm2Deg; + const Double_t length = length0 * fMm2Deg; + + // create an array to pass the variables to MDisp::Calc + Int_t numimagevars = 6; + TArrayD imagevars(numimagevars); + imagevars[0] = log10(size); + imagevars[1] = width; + imagevars[2] = length; + imagevars[3] = conc; + imagevars[4] = leakage1; + imagevars[5] = leakage2; + + // compute Disp + Double_t disp = fDisp->Calc(imagevars); + + // save value into MImageParDisp container + fImageParDisp->SetDisp(disp); + fImageParDisp->SetReadyToSave(); + + return kTRUE; +} +//=========================================================================== + + + + Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MDispCalc.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MDispCalc.h (revision 5671) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MDispCalc.h (revision 5671) @@ -0,0 +1,89 @@ +#ifndef MARS_MDispCalc +#define MARS_MDispCalc + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#ifndef ROOT_TArrayD +#include +#endif + +#ifndef ROOT_TArrayI +#include +#endif + +class MSrcPosCam; +class MHillas; +class MHillasExt; +class MNewImagePar; +class MMcEvt; +class MPointingPos; +class MImageParDisp; +class MDisp; +class MHMatrix; +class MParList; + +class MDispCalc : public MTask +{ +private: + + MSrcPosCam *fSrcPos; + MHillas *fHil; + MHillasExt *fHilExt; + MNewImagePar *fNewPar; + MMcEvt *fMcEvt; + MPointingPos *fPointing; // NOTE: take (theta,phi) from this container + // when it is available (MC calibrated files + // or data files) + + MImageParDisp *fImageParDisp; //! output container for Disp + MDisp *fDisp; // container for Disp parameters + + TString fImageParDispName; // name of container to store Disp + TString fDispName; // name of container for Disp parameters + + Double_t fMm2Deg; // conversion factor from mm to deg + + MHMatrix *fMatrix; // matrix defined to have as much columns as + // variables wanted to be accessible during + // Disp optimization, and as much rows as + // events used for the optimization + // (filled at MFindDisp) + + TArrayI fMap; // array to store the matrix mapping column + // numbers corresponding to each added variable + + + Int_t PreProcess(MParList *pList); + Int_t Process(); + + Double_t GetVal(Int_t i) const; // get value of variable entered in the matrix + // at column fMap[i] + + +public: + + MDispCalc(const char *imagepardispname = "MImageParDisp", + const char *dispname = "MDisp", + const char *name=NULL, const char *title=NULL); + + void InitMapping(MHMatrix *mat); // define the matrix of variables + MHMatrix* GetMatrixMap(TArrayI &map); // get matrix and its mapping array + + ClassDef(MDispCalc, 0) // Task to evaluate Disp +}; + +#endif + + + + + + + + + + + + Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MFDisp.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MFDisp.cc (revision 5671) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MFDisp.cc (revision 5671) @@ -0,0 +1,291 @@ +/* ======================================================================== *\ +! +! * +! * 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 "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; + } + + 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 = fNewImgPar->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 5671) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MFDisp.h (revision 5671) @@ -0,0 +1,88 @@ +#ifndef MARS_MFDisp +#define MARS_MFDisp + +#ifndef MARS_MFilter +#include "MFilter.h" +#endif + +class MHillas; +class MHillasSrc; +class MNewImagePar; + +class MFDisp : public MFilter +{ +private: + + MHillas *fHil; + MHillasSrc *fHilSrc; + 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 5671) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MFindDisp.cc (revision 5671) @@ -0,0 +1,1073 @@ +/* ======================================================================== *\ +! +! * +! * 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 "MDisp.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 "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 (MDisp::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 oject (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 + MDisp *disp = (MDisp*)plistfcn->FindObject("MDisp"); + if (!disp) + { + gLog << "fcnDisp : MDisp object '" << "MDisp" + << "' not found... aborting" << endl; + return; + } + + MHDisp *hdisp = (MHDisp*)plistfcn->FindObject("MHDisp"); + if (!hdisp) + { + gLog << "fcnDisp : MHDisp object '" << "MHDisp" + << "' not found... aborting" << endl; + return; + } + + // + // transfer current parameter values to MDisp + // + // 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); + + disp->SetParameters(TArrayD(fNparx, par)); + + // checking that parameters have been properly set + TArrayD checkparameters = disp->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 Chi^2 value for the current values of the parameters + f = hdisp->GetChi2(); +} + + +// -------------------------------------------------------------------------- +// +// 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; + *fLog << inf << "MFindDisp::MFindDisp; address of MFDisp object = " + << fDispFilter << endl; + + // matrices to contain the training/test samples + fMatrixTrain = new MHMatrix("MatrixTrain"); + fMatrixTest = new MHMatrix("MatrixTest"); + + // objects of MDispCalc to which these matrices are attached + fDispCalcTrain = new MDispCalc("DispTrain"); + fDispCalcTest = new MDispCalc("DispTest"); + + // define columns of matrices + fDispCalcTrain->InitMapping(fMatrixTrain); + fDispCalcTest->InitMapping(fMatrixTest); +} + +// -------------------------------------------------------------------------- +// +// Default destructor. +// +MFindDisp::~MFindDisp() +{ + delete fCam; + delete fMatrixTrain; + delete fMatrixTest; + delete fDispCalcTrain; + delete fDispCalcTest; +} + + +// -------------------------------------------------------------------------- +// +// Define the matrix 'fMatrixTrain' 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 Matrix 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 filltrain(fMatrixTrain); + filltrain.SetFilter(&seltrain); + filltrain.SetName("fillMatrixTrain"); + + //****************************** + // entries in MParList + + plist.AddToList(&tlist); + plist.AddToList(fCam); + plist.AddToList(fMatrixTrain); + + //****************************** + // entries in MTaskList + + tlist.AddToList(&read); + if (fDispFilter != NULL) + tlist.AddToList(&contdisp); + tlist.AddToList(&seltrain); + tlist.AddToList(&filltrain); + + //****************************** + + 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 Matrix + fMatrixTrain->Print("SizeCols"); + + // check that number of generated events is compatible with the resquested + Int_t howmanygenerated = fMatrixTrain->GetM().GetNrows(); + if (TMath::Abs(howmanygenerated-howmanytrain) > TMath::Sqrt(9.*howmanytrain)) + { + *fLog << "MFindDisp::DefineTrainMatrix; no.of generated events (" + << howmanygenerated + << ") is incompatible with the no.of requested events (" + << howmanytrain << ")" << endl; + } + + *fLog << "TRAINING Matrix was filled" << endl; + *fLog << "=============================================" << endl; + + + //-------------------------- + // write out training matrix + + if (filetrain != "") + { + TFile filetr(filetrain, "RECREATE"); + fMatrixTrain->Write(); + filetr.Close(); + + *fLog << "MFindDisp::DefineTrainMatrix; Training matrix was written onto file '" + << filetrain << "'" << endl; + } + + if (display != NULL) + delete display; + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Define the matrix 'fMatrixTest' 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 Matrix 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 filltest(fMatrixTest); + filltest.SetFilter(&seltest); + filltest.SetName("fillMatrixTest"); + + //****************************** + // entries in MParList + + plist.AddToList(&tlist); + plist.AddToList(fCam); + plist.AddToList(fMatrixTest); + + //****************************** + // entries in MTaskList + + tlist.AddToList(&read); + if (fDispFilter != NULL) + tlist.AddToList(&contdisp); + tlist.AddToList(&seltest); + tlist.AddToList(&filltest); + + //****************************** + + 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 Matrix + fMatrixTest->Print("SizeCols"); + + // check that number of generated events is compatible with the resquested + const Int_t howmanygenerated = fMatrixTest->GetM().GetNrows(); + if (TMath::Abs(howmanygenerated-howmanytest) > TMath::Sqrt(9.*howmanytest)) + { + *fLog << "MFindDisp::DefineTestMatrix; no.of generated events (" + << howmanygenerated + << ") is incompatible with the no.of requested events (" + << howmanytest << ")" << endl; + } + + *fLog << "TEST Matrix was filled" << endl; + *fLog << "=============================================" << endl; + + + //-------------------------- + // write out test matrix + + if (filetest != "") + { + TFile filete(filetest, "RECREATE"); + fMatrixTest->Write(); + filete.Close(); + + *fLog << "MFindDisp::DefineTestMatrix; Test matrix was 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 Matrix 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 filltrain(fMatrixTrain); + filltrain.SetFilter(&split); + filltrain.SetName("fillMatrixTrain"); + + // 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 filltest(fMatrixTest); + filltest.SetName("fillMatrixTest"); + + + //****************************** + // entries in MParList + + plist.AddToList(&tlist); + plist.AddToList(fCam); + plist.AddToList(fMatrixTrain); + plist.AddToList(fMatrixTest); + + //****************************** + // entries in MTaskList + + tlist.AddToList(&read); + if (fDispFilter != NULL) + tlist.AddToList(&contdisp); + tlist.AddToList(&cont); + tlist.AddToList(&filltrain); + tlist.AddToList(&conttrain); + tlist.AddToList(&filltest); + + //****************************** + + 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 Matrix + fMatrixTrain->Print("SizeCols"); + + // check that number of generated events is compatible with the resquested + const Int_t generatedtrain = fMatrixTrain->GetM().GetNrows(); + if (TMath::Abs(generatedtrain-howmanytrain) > TMath::Sqrt(9.*howmanytrain)) + { + *fLog << "MFindDisp::DefineTrainTestMatrix; no.of generated events (" + << generatedtrain + << ") is incompatible with the no.of requested events (" + << howmanytrain << ")" << endl; + } + + + // print the filled Test Matrix + fMatrixTest->Print("SizeCols"); + + // check that number of generated events is compatible with the resquested + const Int_t generatedtest = fMatrixTest->GetM().GetNrows(); + if (TMath::Abs(generatedtest-howmanytest) > TMath::Sqrt(9.*howmanytest)) + { + *fLog << "MFindDisp::DefineTrainTestMatrix; no.of generated events (" + << generatedtest + << ") is incompatible with the no.of requested events (" + << howmanytest << ")" << endl; + } + + + *fLog << "TRAINING and TEST Matrix were filled" << endl; + *fLog << "=============================================" << endl; + + + //-------------------------- + // write out training matrix + + if (filetrain != "") + { + TFile filetr(filetrain, "RECREATE"); + fMatrixTrain->Write(); + filetr.Close(); + + *fLog << "MFindDisp::DefineTrainTestMatrix; Training matrix was written onto file '" + << filetrain << "'" << endl; + } + + //-------------------------- + // write out test matrix + + if (filetest != "") + { + TFile filete(filetest, "RECREATE"); + fMatrixTest->Write(); + filete.Close(); + + *fLog << "MFindDisp::DefineTrainTestMatrix; Test matrix was 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 matrix + + TFile filetr(filetrain); + fMatrixTrain->Read("MatrixTrain"); + fMatrixTrain->Print("SizeCols"); + + *fLog << "MFindDisp::ReadMatrix; Training matrix was read in from file '" + << filetrain << "'" << endl; + filetr.Close(); + + + //-------------------------- + // read in test matrix + + TFile filete(filetest); + fMatrixTest->Read("MatrixTest"); + fMatrixTest->Print("SizeCols"); + + *fLog << "MFindDisp::ReadMatrix; Test matrix was 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 Chi^2 for the current +// parameter values; +// for this - Disp is calculated in the event loop by calling +// MDispCalc::Process() ==> MDisp::Calc() +// - the Chi^2 contributions are summed up in the event loop +// by calling MHDisp::Fill() +// - after the event loop the Chi^2 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 Disp 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 (fMatrixTrain == NULL) + { + *fLog << "MFindDisp::FindParams; training matrix is not defined... aborting" + << endl; + return kFALSE; + } + + if (fMatrixTrain->GetM().GetNrows() <= 0) + { + *fLog << "MFindDisp::FindParams; training matrix has no entries" + << endl; + return kFALSE; + } + + //--------------------------------------------------------- + MParList parlistfcn; + MTaskList tasklistfcn; + + // loop over rows of matrix + MMatrixLoop loop(fMatrixTrain); + + //-------------------------------- + // create container for the Disp parameters + // and set them to their initial values + MDisp disp; + + // take initial values from file parDispInit + if (parDispInit != "") + { + TFile inparam(parDispInit); + disp.Read("MDisp"); + 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; + disp.SetParameters(params); + } + if (steps.GetSize() != 0) + { + *fLog << "MFindDisp::FindParams; initial step sizes are taken from 'steps'" + << endl; + disp.SetStepsizes(steps); + } + } + else + { + *fLog << "MFindDisp::FindParams; initial values and step sizes are taken " + << "from the MDisp constructor" << endl; + } + + // get the matrix pointer and mapping from MDispCalc object to set the same in MHDisp object + TArrayI map; + MHMatrix *tmpmatrix = NULL; + tmpmatrix = fDispCalcTrain->GetMatrixMap(map); + + // 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("DispTrain"); + MHDisp hdisp("DispTrain"); + hdisp.SetMatrixMap(tmpmatrix,map); + // fill the plots for Disp and sum up the Chi^2 contributions + MFillH filldispplots("MHDisp", ""); + + //****************************** + // entries in MParList + + parlistfcn.AddToList(&tasklistfcn); + parlistfcn.AddToList(&disp); + parlistfcn.AddToList(&hdisp); + parlistfcn.AddToList(fCam); + parlistfcn.AddToList(fMatrixTrain); + + //****************************** + // entries in MTaskList + + tasklistfcn.AddToList(&loop); + 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 = disp.GetParameters(); + fStep = disp.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; iGetM().GetNrows() <= 0) + { + *fLog << "MFindDisp::TestParams; test matrix has 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); + MDisp dispin; + dispin.Read("MDisp"); + 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; iGetMatrixMap(map); + + // 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 Chi^2 contributions + MHDisp hdisp1("DispTest"); + hdisp1.SetName("MHDispCorr"); + hdisp1.SetSelectedPos(1); + hdisp1.SetMatrixMap(tmpmatrix,map); + MFillH filldispplots1("MHDispCorr[MHDisp]", ""); + + MHDisp hdisp2("DispTest"); + hdisp2.SetName("MHDispWrong"); + hdisp2.SetSelectedPos(2); + hdisp2.SetMatrixMap(tmpmatrix,map); + MFillH filldispplots2("MHDispWrong[MHDisp]", ""); + + MHDisp hdisp3("DispTest"); + hdisp3.SetName("MHDispM3Long"); + hdisp3.SetSelectedPos(3); + hdisp3.SetMatrixMap(tmpmatrix,map); + MFillH filldispplots3("MHDispM3Long[MHDisp]", ""); + + MHDisp hdisp4("DispTest"); + hdisp4.SetName("MHDispAsym"); + hdisp4.SetSelectedPos(4); + hdisp4.SetMatrixMap(tmpmatrix,map); + MFillH filldispplots4("MHDispAsym[MHDisp]", ""); + + + //****************************** + // entries in MParList + + parlist2.AddToList(&tasklist2); + parlist2.AddToList(&disp); + parlist2.AddToList(&hdisp1); + parlist2.AddToList(&hdisp2); + parlist2.AddToList(&hdisp3); + parlist2.AddToList(&hdisp4); + parlist2.AddToList(fCam); + parlist2.AddToList(fMatrixTest); + + //****************************** + // entries in MTaskList + + tasklist2.AddToList(&loop); + 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 5671) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MFindDisp.h (revision 5671) @@ -0,0 +1,126 @@ +#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 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 *fMatrixTrain; // training matrix + MHMatrix *fMatrixTest; // testing matrix + + MDispCalc *fDispCalcTrain; + MDispCalc *fDispCalcTest; + + 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 *GetMatrixTrain() { return fMatrixTrain; } + MHMatrix *GetMatrixTest() { return fMatrixTest; } + + 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 5671) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MHDisp.cc (revision 5671) @@ -0,0 +1,664 @@ +/* ======================================================================== *\ +! +! * +! * 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 Chi^2 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 "MParList.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MHDisp); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// 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 srcpos) + + fMatrix = NULL; + + //-------------------------------------------------- + // 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]"); + + fHistChi2 = new TH1F("fHistChi2" , + "Distance^2 between Disp and real srcpos", 100, 0., 2.); + fHistChi2->SetDirectory(NULL); + fHistChi2->UseCurrentStyle(); + fHistChi2->SetXTitle("Chi^2 = distance^2 Disp to real srcpos [deg^2]"); + fHistChi2->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]"); + + fHistChi2Energy = new TH2F("fHistChi2Energy", + "Chi^2 vs. Energy", 50, 1., 3., 100, 0., 2.); + fHistChi2Energy->SetDirectory(NULL); + fHistChi2Energy->UseCurrentStyle(); + fHistChi2Energy->SetXTitle("log10(Energy) [GeV]"); + fHistChi2Energy->SetYTitle("Chi^2 = distance^2 Disp to real srcpos [deg^2]"); + + fHistChi2Size = new TH2F("fHistChi2Size", + "Chi^2 vs. Size", 50, 2., 4., 100, 0., 2.); + fHistChi2Size->SetDirectory(NULL); + fHistChi2Size->UseCurrentStyle(); + fHistChi2Size->SetXTitle("log10(Size) [#phot]"); + fHistChi2Size->SetYTitle("Chi^2 = distance^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->UseCurrentStyle(); + 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 fHistChi2; + delete fHistDuDv; + delete fHistChi2Energy; + delete fHistChi2Size; + 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 Chi^2 computing variables + // before each new eventloop + fNumEv = 0; + fSumChi2 = 0.; + fChi2 = 0.; + + fHistEnergy->Reset(); + fHistSize->Reset(); + fHistcosZA->Reset(); + fSkymapXY->Reset(); + fHistChi2->Reset(); + fHistDuDv->Reset(); + fHistChi2Energy->Reset(); + fHistChi2Size->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... aborting." << endl; + return kFALSE; + } + + fPointing = (MPointingPos*)pList->FindObject("MPointingPos"); + if (!fPointing) + { + *fLog << err << "MPointingPos not found... " + << "MMcEvt is going to be used to get Theta and Phi." + << endl; + // return kFALSE; + } + + //------------------------------------------ + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Set which selection algorithm for the Disp estimated source position +// 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 +// (see MDispCalc) +// +void MHDisp::SetMatrixMap(MHMatrix *matrix, TArrayI &map) +{ + fMatrix = matrix; + fMap = map; +} + + +// -------------------------------------------------------------------------- +// +// 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 Chi^2 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.; + Double_t theta = 0.; + Double_t phi = 0.; + + if ( fMatrix || (!fMatrix && fMcEvt) ) + { + energy = fMatrix ? GetVal(13) : fMcEvt->GetEnergy(); + impact = fMatrix ? GetVal(14) : fMcEvt->GetImpact(); + xmax = fMatrix ? GetVal(15) : fMcEvt->GetLongitmax(); + + if (fPointing) + { + theta = fMatrix ? GetVal(16) : fPointing->GetZd(); + phi = fMatrix ? GetVal(17) : fPointing->GetAz(); + } + else + { + theta = fMatrix ? GetVal(16) : fMcEvt->GetTelescopeTheta()*kRad2Deg; + phi = fMatrix ? GetVal(17) : fMcEvt->GetTelescopePhi()*kRad2Deg; + } + } + else + *fLog << "No MMcEvt container available (no Energy,ImpactPar or Xmax)" << endl; + + Double_t xsrcpos0 = fMatrix ? GetVal(0) : fSrcPos->GetX(); + Double_t ysrcpos0 = fMatrix ? GetVal(1) : fSrcPos->GetY(); + Double_t xmean0 = fMatrix ? GetVal(2) : fHil->GetMeanX(); + Double_t ymean0 = fMatrix ? GetVal(3) : fHil->GetMeanY(); + Double_t delta = fMatrix ? GetVal(4) : fHil->GetDelta(); + + Double_t size = fMatrix ? GetVal(5) : fHil->GetSize(); + // Double_t width0 = fMatrix ? GetVal(6) : fHil->GetWidth(); + // Double_t length0 = fMatrix ? GetVal(7) : fHil->GetLength(); + + // Double_t conc = fMatrix ? GetVal(8) : fNewPar->GetConc(); + // Double_t leakage1 = fMatrix ? GetVal(9) : fNewPar->GetLeakage1(); + // Double_t leakage2 = fMatrix ? GetVal(10) : fNewPar->GetLeakage2(); + + Double_t m3long = fMatrix ? GetVal(11) : fHilExt->GetM3Long(); + Double_t asym = fMatrix ? GetVal(12) : fHilExt->GetAsym(); + + //------------------------------------------ + // convert to deg + // Double_t width = width0 * fMm2Deg; + // Double_t length = length0 * fMm2Deg; + 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); + fHistChi2->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 + fHistEnergy->Fill(log10(energy)); + fHistSize->Fill(log10(size)); + fHistcosZA->Fill(cos(theta/kRad2Deg)); + + // to check the size dependence of the optimization + fHistChi2Energy->Fill(log10(energy),d2); + fHistChi2Size->Fill(log10(size),d2); + fHistDuEnergy->Fill(log10(energy),Du); + fHistDuSize->Fill(log10(size),Du); + fHistDvEnergy->Fill(log10(energy),Dv); + fHistDvSize->Fill(log10(size),Dv); + + // variables for Chi^2 calculation (= d^2) + fNumEv += 1; + fSumChi2 += d2; + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Calculates the final Chi^2 of the Disp optimization +// +Bool_t MHDisp::Finalize() +{ + fChi2 = fSumChi2/fNumEv; + *fLog << endl; + *fLog << "MHDisp::Finalize: SumChi2, NumEv = " << fSumChi2 << ", " << fNumEv << endl; + *fLog << "MHDisp::Finalize: Final Chi2 = " << fChi2 << endl; + + fChi2 = fHistChi2->GetMean(); + *fLog << "MHDisp::Finalize: Final Chi2 = " << fChi2 << 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); + fHistChi2->SetTitleOffset(1.2,"Y"); + fHistChi2->DrawClone(opt); + fHistChi2->SetBit(kCanDelete); + gPad->Modified(); + + TProfile *profChi2Energy = fHistChi2Energy->ProfileX("Chi^2 vs. Energy",0,99999,"s"); + profChi2Energy->SetXTitle("log10(Energy) [GeV]"); + profChi2Energy->SetYTitle("Chi^2 = distance^2 Disp to real srcpos [deg^2]"); + pad->cd(5); + gPad->SetBorderMode(0); + profChi2Energy->SetTitleOffset(1.2,"Y"); + profChi2Energy->DrawClone(opt); + profChi2Energy->SetBit(kCanDelete); + gPad->Modified(); + + TProfile *profChi2Size = fHistChi2Size->ProfileX("Chi^2 vs. Size",0,99999,"s"); + profChi2Size->SetXTitle("log10(Size) [#phot]"); + profChi2Size->SetYTitle("Chi^2 = distance^2 Disp to real srcpos [deg^2]"); + pad->cd(6); + gPad->SetBorderMode(0); + profChi2Size->SetTitleOffset(1.2,"Y"); + profChi2Size->DrawClone(opt); + profChi2Size->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->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->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->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->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 5671) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MHDisp.h (revision 5671) @@ -0,0 +1,107 @@ +#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 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 Disp estimated position + // 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 *fHistChi2; // Histogram of the event Chi^2 (= 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 *fHistChi2Energy; // Chi^2 (= d^2) vs. Energy + TH2F *fHistChi2Size; // Chi^2 (= 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 Disp optimization + TArrayI fMap; // array storing the matrix mapping of columns defined + // in MDispCalc::InitMapping + + Double_t GetVal(Int_t i) const; + + Int_t fNumEv; // total number of events + Double_t fSumChi2; // current sum of Chi^2 + Double_t fChi2; // final Chi^2 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); + + Double_t GetChi2() { return fChi2; } + + TH1F *GetHistEnergy() { return fHistEnergy; } + TH1F *GetHistSize() { return fHistSize; } + TH1F *GetHistcosZA() { return fHistcosZA; } + TH2F *GetSkymapXY() { return fSkymapXY; } + TH1F *GetHistChi2() { return fHistChi2; } + TH2F *GetHistDuDv() { return fHistDuDv; } + TH2F *GetHistChi2Energy() { return fHistChi2Energy; } + TH2F *GetHistChi2Size() { return fHistChi2Size; } + 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 +}; + +#endif + + Index: trunk/MagicSoft/Mars/mtemp/mifae/macros/OptimizeDisp.C =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/macros/OptimizeDisp.C (revision 5671) +++ trunk/MagicSoft/Mars/mtemp/mifae/macros/OptimizeDisp.C (revision 5671) @@ -0,0 +1,560 @@ +//************************************************************************ +// +// Authors : Eva Domingo, 12/2004 +// Wolfgang Wittek, 12/2004 +// +// +// Macro for estimating the DISP parameter +// --------------------------------------- +// +// DISP is the distance between the center of gravity of the shower image +// and the estimated source position, assumed to lie on the major +// axis of the shower +// +// In order to get an estimate of DISP +// - one assumes a parametrization DISP = f(a,b,c,... ; s1,s2,s3,... ) +// where s1,s2,s3,... are measurable quantities like +// image parameters, (theta, phi), etc. +// and a,b,c, ... are free parameters +// +// - and determines the free parameters a,b,c ... from MC gamma data +// by minimizing a chi2 which measures the deviation of the estimated +// source position (using DISP) and the true source position +// +// The following classes are used : +// +// MDisp container holding the parameters p[0],p[1],p[2],... +// +// ::Calc member function calculating DISP for a given event; +// it is called by MDispCalc::Process() +// +// MDispCalc task calculating DISP with the parameters stored in MDisp +// +// MFindDisp class with several functions : +// +// ::DefineTrainMatrix \ member functions which generate the training +// ::DefineTestMatrix | and test samples (in the form of matrices) +// ::DefineTrainTestMatrix / from the MC gamma data +// +// ::FindParams member function steering the minimization +// (definition of the fcn function for Minuit, setting up the +// event loop to be executed in each minimization step, +// call to Minuit); +// for the minimization the training matrix is used +// +// ::TestParams member function testing the quality of the DISP estimate; +// for the test the test matrix is used +// +// MHDisp container for histograms which are useful for judging +// the quality of the DISP estimate +// Also task calculating the chi2 contribution of an event +// (in Fill()) +// and computing the final chi2 of the whole event sample +// (in Finalize()) +// +// MFDisp filter to select sample for the optimization +// +// +// The classes are stored in the CVS directory Mars/mtemp/mifae/library +// and Mars/mtemp/mifae/macros +// +// +//************************************************************************ + +void OptimizeDisp() +{ + //************************************************************************ + + Bool_t CMatrix = kFALSE; // Create training and test matrices + Bool_t WOptimize = kFALSE; // Do optimization using the training sample + // and write Disp parameter values + // onto the file parDispfile + // Optimize Disp with : + //TString typeOpt = "Data"; + TString typeOpt = "MC"; + Bool_t RTest = kFALSE; // Test the quality of the Disp estimate + // using the test matrix + Bool_t WDisp = kTRUE; // Make Disp plots for the data of type : + //TString typeInput = "ON"; + //TString typeInput = "OFF"; + TString typeInput = "MC"; + + //************************************************************************ + + gLog.SetNoColors(); + + if (gRandom) + delete gRandom; + gRandom = new TRandom3(0); + + //----------------------------------------------- + //names of files to be read for generating the training and test matrices + const char *filetrain = "..."; + const char *filetest = "..."; + + //----------------------------------------------- + // path for input for Mars + TString inPathON = + "/home/pcmagic14/wittek/CalibData/CrabLow42/2004_10_16/"; + TString inPathOFF = + "/home/pcmagic14/wittek/CalibData/CrabLow42/2004_10_17/"; + TString inPathMC = + "~domingo/MC/"; + + // path for output from Mars + TString outPath = "/mnt/home/pcmagic04/domingo/DispOptimization/"; + + //----------------------------------------------- + // names of files for which Disp plots should be made + const char *offfile = "*OffCrabLow42*"; + const char *onfile = "*CrabLowEn42*"; + const char *mcfile = "starmc_Gammas_30-20_Period21.27ReducedMC.Zbins0-12"; + //----------------------------------------------- + + + //--------------------------------------------------------------------- + gLog << "=====================================================" << endl; + gLog << "Macro OptimizeDisp : Start " << endl; + + gLog << "" << endl; + gLog << "Macro OptimizeDisp : CMatrix, WOptimize, RTest, WDisp = " + << (CMatrix ? "kTRUE" : "kFALSE") << ", " + << (WOptimize ? "kTRUE" : "kFALSE") << ", " + << (RTest ? "kTRUE" : "kFALSE") << ", " + << (WDisp ? "kTRUE" : "kFALSE") << endl; + + + //-------------------------------------------- + // files to contain the matrices (generated from filenameTrain and + // filenameTest) + // + // for the training + TString fileMatrixTrain = outPath; + fileMatrixTrain += "MatrixTrainDisp"; + fileMatrixTrain += ".root"; + gLog << "" << endl; + gLog << "Files containing the Training and Test matrices :" << endl; + gLog << " fileMatrixTrain = " << fileMatrixTrain << endl; + + // for testing + TString fileMatrixTest = outPath; + fileMatrixTest += "MatrixTestDisp"; + fileMatrixTest += ".root"; + gLog << " fileMatrixTest = " << fileMatrixTest << endl; + gLog << "" << endl; + + + //--------------- + // file to contain the optimum Disp parameter values + TString parDispFile = outPath; + // parDispFile += "parDispMC.root"; + parDispFile += "parDispMC_withcuts.root"; + + gLog << "" << endl; + gLog << "File containing the optimum Disp parameter values : parDispFile = " + << parDispFile << endl; + + + // Set the filter cuts to select a sample of events for the Disp optimization + // (islandsmin,islandsmax,usedpixelsmin,usedpixelsmax,corepixelsmin, + // corepixelsmax,sizemin,sizemax,leakage1min,leakage1max,leakage2min, + // leakage2max,lengthmin,widthmin); + MFDisp *fdisp = NULL; + fdisp = new MFDisp; + fdisp->SetCuts(0,2,7,600,0,600,0.,3000.,0.,0.,0.,0.,0.,0.); + fdisp->SetName("FilterSelector2"); + + + // Create the MFindDisp object and set the file to store optimium Disp parameters + MFindDisp finddisp(fdisp); + finddisp.SetFilenameParam(parDispFile); + + + + //====================================================================== + // Create matrices and write them onto files + //====================================================================== + + if (CMatrix) + { + gLog << "-----------------------------------------------------" << endl; + gLog << "Generate the training and test samples" << endl; + + //-------------------------------------------- + // files to be read for optimizing the Disp parameters + // + // for the training + TString filenameTrain = inPathMC; + filenameTrain += mcfile; + filenameTrain += ".root"; + Int_t howManyTrain = 30000; + gLog << "filenameTrain = " << filenameTrain << ", howManyTrain = " + << howManyTrain << endl; + + // for testing + TString filenameTest = inPathMC; + filenameTest += mcfile; + filenameTest += ".root"; + Int_t howManyTest = 30000; + gLog << "filenameTest = " << filenameTest << ", howManyTest = " + << howManyTest << endl; + + + //-------------------------------------------- + // For the events in the matrices define requested distribution + // in some quantities; this may be a 1-dim, 2-dim or 3-dim distribution + + /* + // Define the cos(theta) distribution + TString mname("costheta"); + MBinning bincost("Binning"+mname); + bincost.SetEdges(10, 0., 1.0); + + //MH3 mh3("cos((MPointingPos.fZd)/kRad2Deg)"); + MH3 mh3("cos(MMcEvt.fTelescopeTheta)"); + mh3.SetName(mname); + MH::SetBinning(&mh3.GetHist(), &bincost); + + for (Int_t i=1; i<=mh3.GetNbins(); i++) + mh3.GetHist().SetBinContent(i, 1.0); + */ + + // Define the log10(Etrue) distribution + TString mh3name("log10Etrue"); + MBinning binslogE("Binning"+mh3name); + binslogE.SetEdges(50, 1., 3.); + + MH3 mh3("log10(MMcEvt.fEnergy)"); + mh3.SetName(mh3name); + MH::SetBinning(&mh3.GetHist(), &binslogE); + + // Flat energy distribution + // for (Int_t i=1; i<=mh3.GetNbins(); i++) + // mh3.GetHist().SetBinContent(i, 1.0); + + // Power law energy distribution + // for (Int_t i=1; i<=mh3.GetNbins(); i++) + // { + // Double_t weight = pow((Double_t)i, -1.7); + // mh3.GetHist().SetBinContent(i, weight); + // } + + /* + // define the 'cos(Theta) vs. log10(Etrue)' distribution + TString mh3name("cosThetaVslog10Etrue"); + MBinning binslogE("Binning"+mh3name+"X"); + binslogE.SetEdges(18, 1.5, 2.2); + MBinning binscost("Binning"+mh3name+"Y"); + binscost.SetEdges(10, 0., 1.0); + + // MH3 mh3("log10(MMcEvt.fEnergy)", "cos((MPointingPos.fZd)/kRad2Deg)"); + MH3 mh3("log10(MMcEvt.fEnergy)", "cos(MMcEvt.fTelescopeTheta)"); + mh3.SetName(mh3name); + MH::SetBinning((TH2*)&mh3.GetHist(), &binslogE, &binscost); + + for (Int_t i=1; i<=((TH2*)mh3.GetHist()).GetNbinsX(); i++) + { + // Double_t weight = pow((Double_t)i, -1.7); + for (Int_t j=1; j<=((TH2*)mh3.GetHist()).GetNbinsY(); j++) + { + // mh3.GetHist().SetBinContent(i, j, weight); + mh3.GetHist().SetBinContent(i, j, 1.); + } + } + */ + + //-------------------------- + // Training and test samples are generated from the same input files + if (filenameTrain == filenameTest) + { + if ( !finddisp.DefineTrainTestMatrix( + filenameTrain, mh3, + howManyTrain, howManyTest, + fileMatrixTrain, fileMatrixTest, 0) ) + { + *fLog << "OptimizeDisp.C : DefineTrainTestMatrix failed" << endl; + return; + } + } + + //-------------------------- + // Training and test samples are generated from different input files + else + { + if ( !finddisp.DefineTrainMatrix(filenameTrain, mh3, + howManyTrain, fileMatrixTrain, 0) ) + { + *fLog << "OptimizeDisp.C : DefineTrainMatrix failed" << endl; + return; + } + + if ( !finddisp.DefineTestMatrix( filenameTest, mh3, + howManyTest, fileMatrixTest, 0) ) + { + *fLog << "OptimizeDisp.C : DefineTestMatrix failed" << endl; + return; + } + } + + gLog << "Generation of training and test samples finished" << endl; + gLog << "-----------------------------------------------------" << endl; + } + + + + //====================================================================== + // Optimize Disp parameters using the training sample + // + // the initial values of the parameters are taken + // - from the file parDispInit (if != "") + // - or from the arrays params and steps (if their sizes are != 0) + // - or from the MDisp constructor + //====================================================================== + + if (WOptimize) + { + gLog << "-----------------------------------------------------" << endl; + gLog << "Optimize the Disp parameters over a " + << typeOpt << " sample, using the training matrix" << endl; + + // Read training matrix from file + finddisp.ReadMatrix(fileMatrixTrain, fileMatrixTest); + + //-------------------------------------------- + // file which contains the initial values for the Disp parameters; + // if parDispInit ="" + // - the initial values are taken from the arrays params and steps + // if their sizes are different from 0 + // - otherwise they are taken from the MDisp constructor + + TString parDispInit = outPath; + //parDispInit += "parDispInit.root"; + parDispInit = ""; + + //-------------------------------------------- + + TArrayD params(0); + TArrayD steps(0); + + if (parDispInit == "") + { + Double_t vparams[5] = { + // p[0], p[1], p[2], p[3], p[4], p[5] + 1.0, 0.6, -0.8, -0.8, -1.2/*, 0.5*/}; + // 0.5, 0., 0.03, 0., 0./*, 0.5*/}; + // 0.8, 0., 0., 0., 0./*, 0.5*/}; + + Double_t vsteps[5] = { + // dp[0], dp[1], dp[2], dp[3], dp[4], dp[5] + 0.01, 0.005, 0.005, 0.005, 0.01/*, 0.001*/}; + // 0.01, 0.01, 0.01, 0.01, 0.01/*, 0.001*/}; + // 0.01, 0.01, 0.01, 0.01, 0.01/*, 0.001*/}; + + params.Set(5, vparams); + steps.Set (5, vsteps ); + } + + + Bool_t rf; + rf = finddisp.FindParams(parDispInit, params, steps); + + if (!rf) + { + gLog << "OptimizeDisp.C : optimization of the Disp parameters failed" << endl; + return; + } + + gLog << "Optimization of Disp parameters finished" << endl; + gLog << "-----------------------------------------------------" << endl; + } + + + + //====================================================================== + // Test the Disp parameters on the test sample + //====================================================================== + + if (RTest) + { + // Read test matrix from file + finddisp.ReadMatrix(fileMatrixTrain, fileMatrixTest); + + gLog << "-----------------------------------------------------" << endl; + gLog << "Test the Disp parameters using on the test matrix" << endl; + Bool_t rt = finddisp.TestParams(); + if (!rt) + { + gLog << "Test of the Disp parameters on the test matrix failed" + << endl; + } + gLog << "Test of the Disp parameters finished" << endl; + gLog << "-----------------------------------------------------" << endl; + } + + + + + //====================================================================== + // Make Disp plots + //====================================================================== + + if (WDisp) + { + gLog << "" << endl; + gLog << "-----------------------------------------------------" << endl; + gLog << "Make plots for Disp for data of the type " << typeInput << endl; + + //-------------------------------------------- + // type of data to be read (ON, OFF or MC) + if (typeInput == "ON") + TString file(onfile); + else if (typeInput == "OFF") + TString file(offfile); + else if (typeInput == "MC") + TString file(mcfile); + + // name of input root file + TString filenameData = inPathMC; + filenameData += mcfile; + filenameData += ".root"; + gLog << "Input file '" << filenameData << endl; + + //---------------------------------------------------- + // read in optimum Disp parameter values + + TFile inparam(parDispFile); + MDisp dispin; + dispin.Read("MDisp"); + inparam.Close(); + + gLog << "Disp parameter values were read in from file '" + << parDispFile << "'" << endl; + + TArrayD dispPar; + dispPar = dispin.GetParameters(); + + TArrayD dispStep; + dispStep = dispin.GetStepsizes(); + + gLog << "optimum parameter values for calculating Disp : " << endl; + for (Int_t i=0; i