Index: trunk/MagicSoft/Mars/mtemp/mifae/Changelog =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/Changelog (revision 5879) +++ trunk/MagicSoft/Mars/mtemp/mifae/Changelog (revision 5904) @@ -18,4 +18,20 @@ -*-*- END OF LINE -*-*- + + 2005/01/20 Eva Domingo + * library/MDispParameters.[cc,h] + - Added substituting previous MDisp.[cc,h] classes. + Old MDisp::Calc function moved to MDispCalc task. + * library/MDispCalc.[cc,h] + - Now includes the DISP parameterization, defined in + MDispCalc::Calc. MDispCalc::InitMapping now defines + training and test matrices with only the columns of + variables needed in the DISP expression. + * library/MHDisp.[cc,h] + - Changed the Chi2 name by MinPar (to clarify it accounts for + the parameter asked to be minimized in the Disp optimization). + Added MHDisp::InitMapping, defines matrices with only the columns + of variables needed for computing the minimization parameter + and for filling the DISP histograms. 2005/01/18 Javier Rico & Markus Gaug Index: trunk/MagicSoft/Mars/mtemp/mifae/library/IFAELinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/IFAELinkDef.h (revision 5879) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/IFAELinkDef.h (revision 5904) @@ -27,5 +27,5 @@ #pragma link C++ class MTopologyCalc+; #pragma link C++ class MImageParDisp+; -#pragma link C++ class MDisp+; +#pragma link C++ class MDispParameters+; #pragma link C++ class MDispCalc+; #pragma link C++ class MHDisp+; Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MDisp.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MDisp.cc (revision 5879) +++ (revision ) @@ -1,203 +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 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// // -// MDisp // -// // -// this is the container for the Disp parameterization // -// // -///////////////////////////////////////////////////////////////////////////// -#include "MDisp.h" - -#include -#include - -#include - -#include "MImageParDisp.h" - -#include "MParList.h" - -#include "MLog.h" -#include "MLogManip.h" - - -ClassImp(MDisp); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// constructor -// -MDisp::MDisp(const char *name, const char *title) - : fLogSize0(3.), fLength0(0.1), fWidth0(0.05), fConc0(0.35), - fLeakage10(0.05),fLeakage20(0.1) -{ - fName = name ? name : "MDisp"; - fTitle = title ? title : "Container for the Disp parameters"; - - // set arrays size to the number of parameters used - fParameters.Set(5); - fStepsizes.Set(5); - - // set Disp parameters to their default values - InitParameters(); -} - - -// -------------------------------------------------------------------------- -// -// set default values for the Disp parameters -// -void MDisp::InitParameters() -{ - //--------------------------------------------------- - // these are the default values - - fParameters[0] = 1.0; - fParameters[1] = 0.6; - fParameters[2] = -0.8; - fParameters[3] = -0.8; - fParameters[4] = -1.2; - - - //--------------------------------------------------- - // 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.006; - fStepsizes[2] = 0.008; - fStepsizes[3] = 0.008; - fStepsizes[4] = 0.012; -} - - -// -------------------------------------------------------------------------- -// -// Set the parameter values from the array 'd' -// -// -Bool_t MDisp::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 MDisp::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; -} - - -// -------------------------------------------------------------------------- -// -// Set the Disp parameterization and Calculate Disp -// -// -Double_t MDisp::Calc(TArrayD &imagevar) -{ - // get parameters - Int_t numpar = fParameters.GetSize(); - TArrayD p(numpar); - for (Int_t i=0; i -#endif - -class MDisp : public MParContainer -{ -private: - - TArrayD fParameters; // Disp parameters - TArrayD fStepsizes; // step sizes of Disp parameters - - Double_t fLogSize0; // Variables introduced in the Disp expression - Double_t fLength0; // to shift the minimization around the mean - Double_t fWidth0; // values of the variables, so it makes easier - Double_t fConc0; // to MINUIT to converge. Default values are set - Double_t fLeakage10; // in the constructor (to their standard mean - Double_t fLeakage20; // distribution values) but can be changed with - // the set functions. - -public: - - MDisp(const char *name=NULL, const char *title=NULL); - - void InitParameters(); - - Bool_t SetParameters(const TArrayD &d); - Bool_t SetStepsizes(const TArrayD &d); - - const TArrayD &GetParameters() const { return fParameters; } - const TArrayD &GetStepsizes() const { return fStepsizes; } - - void SetLogSize0(Double_t newmeanval) { fLogSize0 = newmeanval; } - void SetWidth0(Double_t newmeanval) { fWidth0 = newmeanval; } - void SetLength0(Double_t newmeanval) { fLength0 = newmeanval; } - void SetConc0(Double_t newmeanval) { fConc0 = newmeanval; } - void SetLeakage10(Double_t newmeanval) { fLeakage10 = newmeanval; } - void SetLeakage20(Double_t newmeanval) { fLeakage20 = newmeanval; } - - Double_t Calc(TArrayD &imagevar); - - ClassDef(MDisp, 1) // Container for the Disp parameters -}; - -#endif - - - Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MDispCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MDispCalc.cc (revision 5879) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MDispCalc.cc (revision 5904) @@ -28,8 +28,8 @@ // 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 // +// This task calculates Disp with a given parameterization // +// (parameters are stored in the MDispParameters container) // +// Also the matrix of variables to be used in the Disp // +// parameterization is defined // // // // // @@ -48,5 +48,5 @@ #include "MPointingPos.h" #include "MImageParDisp.h" -#include "MDisp.h" +#include "MDispParameters.h" #include "MHMatrix.h" @@ -62,11 +62,20 @@ using namespace std; +static const Int_t nPars=5; // number of parameters used in the Disp expression + +enum dispVar_t {kSize,kWidth,kLength,kConc,kLeakage1,kLeakage2,kTotVar}; // enum variables for the + // matrix columns mapping + // -------------------------------------------------------------------------- // // constructor // -MDispCalc::MDispCalc(const char *imagepardispname, const char *dispname, +MDispCalc::MDispCalc(const char *imagepardispname, + const char *dispparametersname, const char *name, const char *title) - : fImageParDispName(imagepardispname), fDispName(dispname) + : fImageParDispName(imagepardispname), + fDispParametersName(dispparametersname), + fLogSize0(3.), fLength0(0.1), fWidth0(0.05), + fConc0(0.35), fLeakage10(0.05), fLeakage20(0.1) { fName = name ? name : "MDispCalc"; @@ -76,5 +85,5 @@ // initialize mapping array dimension to the number of columns of fMatrix - fMap.Set(18); + fMap.Set(kTotVar); } @@ -106,10 +115,12 @@ } - fDisp = (MDisp*)pList->FindObject(fDispName, "MDisp"); - if (!fDisp) - { - *fLog << err << fDispName << " [MDisp] not found... aborting." << endl; - return kFALSE; - } + fDispParameters = (MDispParameters*)pList->FindObject(fDispParametersName, "MDispParameters"); + if (!fDispParameters) + { + *fLog << err << fDispParametersName << " [MDispParameters] not found... aborting." << endl; + return kFALSE; + } + fDispParameters->GetParameters().Set(nPars); + fDispParameters->GetStepsizes().Set(nPars); @@ -160,8 +171,6 @@ if (!fPointing) { - *fLog << err << "MPointingPos not found... " - << "MMcEvt is going to be used to get Theta and Phi." - << endl; - // return kFALSE; + *fLog << err << "MPointingPos not found... aborting." << endl; + return kFALSE; } @@ -175,9 +184,9 @@ // // 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 +// given containers for the Disp parameterization. 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 +// MMatrixLoop), MDispCalc::Calc will take the values from the matrix // instead of the containers. // @@ -189,47 +198,11 @@ 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; + fMap[kSize] = fMatrix->AddColumn("MHillas.fSize"); + fMap[kWidth] = fMatrix->AddColumn("MHillas.fWidth"); + fMap[kLength] = fMatrix->AddColumn("MHillas.fLength"); + + fMap[kConc] = fMatrix->AddColumn("MNewImagePar.fConc"); + fMap[kLeakage1] = fMatrix->AddColumn("MNewImagePar.fLeakage1"); + fMap[kLeakage2] = fMatrix->AddColumn("MNewImagePar.fLeakage2"); } @@ -258,10 +231,10 @@ { // 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(); + const Double_t size = fMatrix ? GetVal(kSize) : fHil->GetSize(); + const Double_t width0 = fMatrix ? GetVal(kWidth) : fHil->GetWidth(); + const Double_t length0 = fMatrix ? GetVal(kLength) : fHil->GetLength(); + const Double_t conc = fMatrix ? GetVal(kConc) : fNewPar->GetConc(); + const Double_t leakage1 = fMatrix ? GetVal(kLeakage1) : fNewPar->GetLeakage1(); + const Double_t leakage2 = fMatrix ? GetVal(kLeakage2) : fNewPar->GetLeakage2(); // convert to deg @@ -269,16 +242,15 @@ 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; + // create an array to pass the variables to MDispCalc::Calc + TArrayD imagevars(kTotVar); + imagevars[kSize] = log10(size); + imagevars[kWidth] = width; + imagevars[kLength] = length; + imagevars[kConc] = conc; + imagevars[kLeakage1] = leakage1; + imagevars[kLeakage2] = leakage2; // compute Disp - Double_t disp = fDisp->Calc(imagevars); + Double_t disp = Calc(imagevars); // save value into MImageParDisp container @@ -288,4 +260,65 @@ return kTRUE; } + + +// -------------------------------------------------------------------------- +// +// Set the Disp parameterization and Calculate Disp +// +// +Double_t MDispCalc::Calc(TArrayD &imagevar) +{ + // get parameters + const TArrayD& p = fDispParameters->GetParameters(); + + // get image variables to be used in the Disp parameterization + Double_t logsize = imagevar[0]; + Double_t width = imagevar[1]; + Double_t length = imagevar[2]; + // Double_t conc = imagevar[3]; + // Double_t leakage1 = imagevar[4]; + // Double_t leakage2 = imagevar[5]; + + // Disp parameterization to be optimized + // Note: fLogSize0, fLength0... variables are introduced to uncorrelate as much + // as possible the parameters in the Disp expression, with the purpose of + // helping the minimization algorithm to converge. They are set approx. + // to their distribution mean value in the MDisp constructor, but can be + // changed using the corresponding set function. + // + + // Double_t disp = p[0] + p[1]*(logsize-fLogSize0) + // + (p[2] + p[3]*(logsize-fLogSize0))*width/length; + + Double_t disp = p[0] + p[1]*(logsize-fLogSize0) + p[4]*(length-fLength0) + + (p[2] + p[3]*(logsize-fLogSize0))*width/length; + + // Double_t disp = p[0] + p[1]*(logsize-fLogSize0) + p[4]*(length-fLength0) + // + (p[2] + p[3]*(logsize-fLogSize0))*length/width; + + // Double_t disp = p[0] + p[1]*(logsize-fLogSize0) + p[4]*(length-fLength0) + // + (p[2] + p[3]*(logsize-fLogSize0) + p[5]*(length-fLength0))*width/length; + + // Double_t disp = p[0] + p[1]*(logsize-fLogSize0) + p[4]*(width-fWidth0) + // + (p[2] + p[3]*(logsize-fLogSize0))*width/length; // + p[5]*(width-fWidth0))*width/length; + + // Double_t disp = p[0] + p[1]*(logsize-fLogSize0) + p[4]*(conc-fConc0) + // + (p[2] + p[3]*(logsize-fLogSize0))*width/length; // + p[5]*(conc-fConc0))*width/length; + + // Double_t disp = ( p[0] + p[1]*(logsize-fLogSize0) + // + p[2]*pow(logsize-fLogSize0,2) + // + p[3]*pow(logsize-fLogSize0,3) + // + p[4]*pow(logsize-fLogSize0,4) ) + // *( 1/(1+width/length) ); + + // Double_t disp = ( p[0] + p[1]*(logsize-fLogSize0) + // + p[2]*pow(logsize-fLogSize0,2) + // + p[3]*pow(logsize-fLogSize0,3) + // + p[4]*pow(logsize-fLogSize0,4) ) + // *(1-width/length); + + return disp; +} + //=========================================================================== Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MDispCalc.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MDispCalc.h (revision 5879) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MDispCalc.h (revision 5904) @@ -21,5 +21,5 @@ class MPointingPos; class MImageParDisp; -class MDisp; +class MDispParameters; class MHMatrix; class MParList; @@ -29,29 +29,35 @@ 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 + MSrcPosCam *fSrcPos; + MHillas *fHil; + MHillasExt *fHilExt; + MNewImagePar *fNewPar; + MMcEvt *fMcEvt; + MPointingPos *fPointing; + + MImageParDisp *fImageParDisp; //! output container for Disp value + MDispParameters *fDispParameters; //! container for Disp parameters TString fImageParDispName; // name of container to store Disp - TString fDispName; // name of container for Disp parameters + TString fDispParametersName; // name of container for Disp parameters + + Double_t fLogSize0; // Variables introduced in the Disp expression + Double_t fLength0; // to shift the minimization around the mean + Double_t fWidth0; // values of the variables, so it makes easier + Double_t fConc0; // to MINUIT to converge. Default values are set + Double_t fLeakage10; // in the constructor (to their standard mean + Double_t fLeakage20; // distribution values) but can be changed with + // the corresponding set function. 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) + // variables wanted to parameterize Disp, + // 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 + TArrayI fMap; // array storing the matrix mapping column + // numbers corresponding to each variable + // added to fMatrix @@ -59,6 +65,6 @@ Int_t Process(); - Double_t GetVal(Int_t i) const; // get value of variable entered in the matrix - // at column fMap[i] + Double_t GetVal(Int_t i) const; // get value of variable entered in the + // matrix at column fMap[i] @@ -70,5 +76,14 @@ void InitMapping(MHMatrix *mat); // define the matrix of variables - MHMatrix* GetMatrixMap(TArrayI &map); // get matrix and its mapping array + // needed for the Disp parameterization + + Double_t Calc(TArrayD &imagevar); // calculate Disp with a given expression + + void SetLogSize0(Double_t newmeanval) { fLogSize0 = newmeanval; } + void SetWidth0(Double_t newmeanval) { fWidth0 = newmeanval; } + void SetLength0(Double_t newmeanval) { fLength0 = newmeanval; } + void SetConc0(Double_t newmeanval) { fConc0 = newmeanval; } + void SetLeakage10(Double_t newmeanval) { fLeakage10 = newmeanval; } + void SetLeakage20(Double_t newmeanval) { fLeakage20 = newmeanval; } ClassDef(MDispCalc, 0) // Task to evaluate Disp Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MDispParameters.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MDispParameters.cc (revision 5904) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MDispParameters.cc (revision 5904) @@ -0,0 +1,132 @@ +/* ======================================================================== *\ +! +! * +! * 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 Disp parameters to their default values + InitParameters(); +} + + +// -------------------------------------------------------------------------- +// +// set default values for the Disp parameters +// +void MDispParameters::InitParameters() +{ + //--------------------------------------------------- + // these are the default values + + fParameters[0] = 1.0; + fParameters[1] = 0.6; + fParameters[2] = -0.8; + fParameters[3] = -0.8; + fParameters[4] = -1.2; + + + //--------------------------------------------------- + // 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.006; + fStepsizes[2] = 0.008; + fStepsizes[3] = 0.008; + fStepsizes[4] = 0.012; +} + + +// -------------------------------------------------------------------------- +// +// 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 5904) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MDispParameters.h (revision 5904) @@ -0,0 +1,40 @@ +#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/MFindDisp.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MFindDisp.cc (revision 5879) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MFindDisp.cc (revision 5904) @@ -48,5 +48,5 @@ #include "MBinning.h" #include "MContinue.h" -#include "MDisp.h" +#include "MDispParameters.h" #include "MDispCalc.h" #include "MDataElement.h" @@ -72,4 +72,5 @@ #include "MMatrixLoop.h" #include "MMinuitInterface.h" +#include "MParameters.h" #include "MParList.h" #include "MProgressBar.h" @@ -95,5 +96,5 @@ // // - the parameters to be varied in the minimization are the parameters -// appearing in the parametrization of Disp (MDisp::Calc()) +// appearing in the parametrization of Disp (MDispCalc::Calc()) // //------------------------------------------------------------------------ @@ -104,5 +105,5 @@ cout << "entry fcnDisp" << endl; - // save pointer to the MINUIT oject (for optimizing Disp) for the case + // 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; @@ -116,8 +117,8 @@ // get needed Disp containers from the existing parList - MDisp *disp = (MDisp*)plistfcn->FindObject("MDisp"); - if (!disp) - { - gLog << "fcnDisp : MDisp object '" << "MDisp" + MDispParameters *dispparams = (MDispParameters*)plistfcn->FindObject("MDispParameters"); + if (!dispparams) + { + gLog << "fcnDisp : MDispParameters object '" << "MDispParameters" << "' not found... aborting" << endl; return; @@ -132,6 +133,8 @@ } + MParameterD *minpar = (MParameterD*)plistfcn->FindObject("MinimizationParameter"); + // - // transfer current parameter values to MDisp + // transfer current parameter values to MDispParameters // // Attention : npar is the number of variable parameters @@ -142,8 +145,8 @@ gMinuit->mnstat(fMin, fEdm, fErrdef, fNpari, fNparx, fIstat); - disp->SetParameters(TArrayD(fNparx, par)); + dispparams->SetParameters(TArrayD(fNparx, par)); // checking that parameters have been properly set - TArrayD checkparameters = disp->GetParameters(); + TArrayD checkparameters = dispparams->GetParameters(); gLog << "fcnDisp : fNpari, fNparx =" << fNpari << ", " << fNparx << endl; @@ -155,6 +158,6 @@ // - // loop over the events in the training matrix to compute the Chi^2 - // for the current parameter values + // loop over the events in the training matrices to compute the parameter + // to minimize for the current Disp parameter values evtloopfcn->Eventloop(); @@ -162,6 +165,14 @@ //------------------------------------------- - // get the Chi^2 value for the current values of the parameters - f = hdisp->GetChi2(); + // 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; + } } @@ -184,19 +195,28 @@ // (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 + 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"); fDispCalcTest = new MDispCalc("DispTest"); + // 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 - fDispCalcTrain->InitMapping(fMatrixTrain); - fDispCalcTest->InitMapping(fMatrixTest); + // Train matrix + fDispCalcTrain->InitMapping(fMatrixTrainCalc); + fHDispTrain->InitMapping(fMatrixTrainHists); + // Test matrix + fDispCalcTest->InitMapping(fMatrixTestCalc); + fHDispTest->InitMapping(fMatrixTestHists); } + // -------------------------------------------------------------------------- @@ -207,8 +227,12 @@ { delete fCam; - delete fMatrixTrain; - delete fMatrixTest; + delete fMatrixTrainCalc; + delete fMatrixTrainHists; + delete fMatrixTestCalc; + delete fMatrixTestHists; delete fDispCalcTrain; delete fDispCalcTest; + delete fHDispTrain; + delete fHDispTest; } @@ -216,5 +240,6 @@ // -------------------------------------------------------------------------- // -// Define the matrix 'fMatrixTrain' for the TRAINING sample +// Define the matrices 'fMatrixTrainCalc' and 'fMatrixTrainHists' +// for the TRAINING sample // // alltogether 'howmanytrain' events are read from file 'nametrain'; @@ -223,6 +248,6 @@ // Bool_t MFindDisp::DefineTrainMatrix( - const TString &nametrain, MH3 &hreftrain, - const Int_t howmanytrain, const TString &filetrain, + const TString &nametrain, MH3 &hreftrain, + const Int_t howmanytrain, const TString &filetrain, Int_t iflag) { @@ -231,5 +256,5 @@ *fLog << "=============================================" << endl; - *fLog << "Fill TRAINING Matrix from file '" << nametrain << endl; + *fLog << "Fill TRAINING Matrices from file '" << nametrain << endl; *fLog << " select " << howmanytrain << " events " << endl; if (!hreftrain.GetHist().GetEntries()==0) @@ -266,7 +291,11 @@ seltrain.SetName("selectTrain"); - MFillH filltrain(fMatrixTrain); - filltrain.SetFilter(&seltrain); - filltrain.SetName("fillMatrixTrain"); + MFillH filltraincalc(fMatrixTrainCalc); + filltraincalc.SetFilter(&seltrain); + filltraincalc.SetName("fillMatrixTrainCalc"); + + MFillH filltrainhists(fMatrixTrainHists); + filltrainhists.SetFilter(&seltrain); + filltrainhists.SetName("fillMatrixTrainHists"); //****************************** @@ -275,5 +304,6 @@ plist.AddToList(&tlist); plist.AddToList(fCam); - plist.AddToList(fMatrixTrain); + plist.AddToList(fMatrixTrainCalc); + plist.AddToList(fMatrixTrainHists); //****************************** @@ -284,5 +314,6 @@ tlist.AddToList(&contdisp); tlist.AddToList(&seltrain); - tlist.AddToList(&filltrain); + tlist.AddToList(&filltraincalc); + tlist.AddToList(&filltrainhists); //****************************** @@ -302,33 +333,46 @@ - // print the filled Training Matrix - fMatrixTrain->Print("SizeCols"); + // print the filled Training Matrices + fMatrixTrainCalc->Print("SizeCols"); + fMatrixTrainHists->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)) + Int_t howmanygeneratedcalc = fMatrixTrainCalc->GetM().GetNrows(); + if (TMath::Abs(howmanygeneratedcalc-howmanytrain) > TMath::Sqrt(9.*howmanytrain)) { *fLog << "MFindDisp::DefineTrainMatrix; no.of generated events (" - << howmanygenerated + << howmanygeneratedcalc << ") is incompatible with the no.of requested events (" << howmanytrain << ")" << endl; } - *fLog << "TRAINING Matrix was filled" << 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 matrix + // write out training matrices if (filetrain != "") { TFile filetr(filetrain, "RECREATE"); - fMatrixTrain->Write(); + fMatrixTrainCalc->Write(); + fMatrixTrainHists->Write(); filetr.Close(); - *fLog << "MFindDisp::DefineTrainMatrix; Training matrix was written onto file '" + *fLog << "MFindDisp::DefineTrainMatrix; Training matrices were written onto file '" << filetrain << "'" << endl; } + if (display != NULL) @@ -341,5 +385,6 @@ // -------------------------------------------------------------------------- // -// Define the matrix 'fMatrixTest' for the TEST sample +// Define the matrices 'fMatrixTestCalc' and 'fMatrixTestHists' +// for the TEST sample // // alltogether 'howmanytest' events are read from file 'nametest' @@ -356,5 +401,5 @@ *fLog << "=============================================" << endl; - *fLog << "Fill TEST Matrix from file '" << nametest << endl; + *fLog << "Fill TEST Matrices from file '" << nametest << endl; *fLog << " select " << howmanytest << " events " << endl; if (!hreftest.GetHist().GetEntries()==0) @@ -390,7 +435,11 @@ seltest.SetName("selectTest"); - MFillH filltest(fMatrixTest); - filltest.SetFilter(&seltest); - filltest.SetName("fillMatrixTest"); + MFillH filltestcalc(fMatrixTestCalc); + filltestcalc.SetFilter(&seltest); + filltestcalc.SetName("fillMatrixTestCalc"); + + MFillH filltesthists(fMatrixTestHists); + filltesthists.SetFilter(&seltest); + filltesthists.SetName("fillMatrixTestHists"); //****************************** @@ -399,5 +448,6 @@ plist.AddToList(&tlist); plist.AddToList(fCam); - plist.AddToList(fMatrixTest); + plist.AddToList(fMatrixTestCalc); + plist.AddToList(fMatrixTestHists); //****************************** @@ -408,5 +458,6 @@ tlist.AddToList(&contdisp); tlist.AddToList(&seltest); - tlist.AddToList(&filltest); + tlist.AddToList(&filltestcalc); + tlist.AddToList(&filltesthists); //****************************** @@ -426,33 +477,45 @@ - // print the filled Test Matrix - fMatrixTest->Print("SizeCols"); + // 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 howmanygenerated = fMatrixTest->GetM().GetNrows(); - if (TMath::Abs(howmanygenerated-howmanytest) > TMath::Sqrt(9.*howmanytest)) + const Int_t howmanygeneratedcalc = fMatrixTestCalc->GetM().GetNrows(); + if (TMath::Abs(howmanygeneratedcalc-howmanytest) > TMath::Sqrt(9.*howmanytest)) { *fLog << "MFindDisp::DefineTestMatrix; no.of generated events (" - << howmanygenerated + << howmanygeneratedcalc << ") is incompatible with the no.of requested events (" << howmanytest << ")" << endl; } - *fLog << "TEST Matrix was filled" << 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 matrix + // write out test matrices if (filetest != "") { TFile filete(filetest, "RECREATE"); - fMatrixTest->Write(); + fMatrixTestCalc->Write(); + fMatrixTestHists->Write(); filete.Close(); - *fLog << "MFindDisp::DefineTestMatrix; Test matrix was written onto file '" + *fLog << "MFindDisp::DefineTestMatrix; Test matrices were written onto file '" << filetest << "'" << endl; } + if (display != NULL) @@ -477,5 +540,5 @@ { *fLog << "=============================================" << endl; - *fLog << "Fill TRAINING and TEST Matrix from file '" << name << endl; + *fLog << "Fill TRAINING and TEST Matrices from file '" << name << endl; *fLog << " select " << howmanytrain << " training and " << howmanytest << " test events " << endl; @@ -521,7 +584,11 @@ split.SetSelectionRatio(prob); - MFillH filltrain(fMatrixTrain); - filltrain.SetFilter(&split); - filltrain.SetName("fillMatrixTrain"); + 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 @@ -530,6 +597,9 @@ conttrain.SetName("ContTrain"); - MFillH filltest(fMatrixTest); - filltest.SetName("fillMatrixTest"); + MFillH filltestcalc(fMatrixTestCalc); + filltestcalc.SetName("fillMatrixTestCalc"); + + MFillH filltesthists(fMatrixTestHists); + filltesthists.SetName("fillMatrixTestHists"); @@ -539,6 +609,8 @@ plist.AddToList(&tlist); plist.AddToList(fCam); - plist.AddToList(fMatrixTrain); - plist.AddToList(fMatrixTest); + plist.AddToList(fMatrixTrainCalc); + plist.AddToList(fMatrixTrainHists); + plist.AddToList(fMatrixTestCalc); + plist.AddToList(fMatrixTestHists); //****************************** @@ -549,7 +621,9 @@ tlist.AddToList(&contdisp); tlist.AddToList(&cont); - tlist.AddToList(&filltrain); + tlist.AddToList(&filltraincalc); + tlist.AddToList(&filltrainhists); tlist.AddToList(&conttrain); - tlist.AddToList(&filltest); + tlist.AddToList(&filltestcalc); + tlist.AddToList(&filltesthists); //****************************** @@ -570,59 +644,79 @@ - // print the filled Train Matrix - fMatrixTrain->Print("SizeCols"); + // 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 generatedtrain = fMatrixTrain->GetM().GetNrows(); - if (TMath::Abs(generatedtrain-howmanytrain) > TMath::Sqrt(9.*howmanytrain)) + const Int_t generatedtraincalc = fMatrixTrainCalc->GetM().GetNrows(); + if (TMath::Abs(generatedtraincalc-howmanytrain) > TMath::Sqrt(9.*howmanytrain)) { *fLog << "MFindDisp::DefineTrainTestMatrix; no.of generated events (" - << generatedtrain + << generatedtraincalc << ") is incompatible with the no.of requested events (" << howmanytrain << ")" << endl; } - - - // print the filled Test Matrix - fMatrixTest->Print("SizeCols"); + 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 generatedtest = fMatrixTest->GetM().GetNrows(); - if (TMath::Abs(generatedtest-howmanytest) > TMath::Sqrt(9.*howmanytest)) + const Int_t generatedtestcalc = fMatrixTestCalc->GetM().GetNrows(); + if (TMath::Abs(generatedtestcalc-howmanytest) > TMath::Sqrt(9.*howmanytest)) { *fLog << "MFindDisp::DefineTrainTestMatrix; no.of generated events (" - << generatedtest + << generatedtestcalc << ") is incompatible with the no.of requested events (" << howmanytest << ")" << endl; } - - - *fLog << "TRAINING and TEST Matrix were filled" << 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 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 + // write out test matrices if (filetest != "") { TFile filete(filetest, "RECREATE"); - fMatrixTest->Write(); + fMatrixTestCalc->Write(); + fMatrixTestHists->Write(); filete.Close(); - *fLog << "MFindDisp::DefineTrainTestMatrix; Test matrix was written onto file '" + *fLog << "MFindDisp::DefineTrainTestMatrix; Test matrices were written onto file '" << filetest << "'" << endl; } @@ -643,11 +737,13 @@ { //-------------------------- - // read in training matrix + // read in training matrices TFile filetr(filetrain); - fMatrixTrain->Read("MatrixTrain"); - fMatrixTrain->Print("SizeCols"); - - *fLog << "MFindDisp::ReadMatrix; Training matrix was read in from file '" + 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(); @@ -655,11 +751,13 @@ //-------------------------- - // read in test matrix + // read in test matrices TFile filete(filetest); - fMatrixTest->Read("MatrixTest"); - fMatrixTest->Print("SizeCols"); - - *fLog << "MFindDisp::ReadMatrix; Test matrix was read in from file '" + 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(); @@ -686,12 +784,12 @@ // // - call TMinuit to do the minimization : -// the fcnDisp function calculates the Chi^2 for the current -// parameter values; +// 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() ==> 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() +// 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) @@ -703,5 +801,5 @@ // - from the file parDispInit (if it is != "") // - or from the arrays params and/or steps -// - or from the Disp constructor +// - or from the DispParameters constructor // //---------------------------------------------------------------------- @@ -721,14 +819,14 @@ - if (fMatrixTrain == NULL) - { - *fLog << "MFindDisp::FindParams; training matrix is not defined... aborting" + if (fMatrixTrainCalc == NULL || fMatrixTrainHists == NULL) + { + *fLog << "MFindDisp::FindParams; training matrices are not defined... aborting" << endl; return kFALSE; } - if (fMatrixTrain->GetM().GetNrows() <= 0) - { - *fLog << "MFindDisp::FindParams; training matrix has no entries" + if (fMatrixTrainCalc->GetM().GetNrows() <= 0 || fMatrixTrainHists->GetM().GetNrows() <= 0) + { + *fLog << "MFindDisp::FindParams; training matrices have no entries" << endl; return kFALSE; @@ -740,10 +838,11 @@ // loop over rows of matrix - MMatrixLoop loop(fMatrixTrain); + MMatrixLoop loopcalc(fMatrixTrainCalc); + MMatrixLoop loophists(fMatrixTrainHists); //-------------------------------- // create container for the Disp parameters // and set them to their initial values - MDisp disp; + MDispParameters dispparams; // take initial values from file parDispInit @@ -751,5 +850,5 @@ { TFile inparam(parDispInit); - disp.Read("MDisp"); + dispparams.Read("MDispParameters"); inparam.Close(); *fLog << "MFindDisp::FindParams; initial values of parameters are taken from file " @@ -764,5 +863,5 @@ *fLog << "MFindDisp::FindParams; initial values of parameters are taken from 'params'" << endl; - disp.SetParameters(params); + dispparams.SetParameters(params); } if (steps.GetSize() != 0) @@ -770,5 +869,5 @@ *fLog << "MFindDisp::FindParams; initial step sizes are taken from 'steps'" << endl; - disp.SetStepsizes(steps); + dispparams.SetStepsizes(steps); } } @@ -776,18 +875,8 @@ { *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 + << "from the MDispParameters constructor" << endl; + } + + // fill the plots for Disp and sum up the Minimization parameter contributions MFillH filldispplots("MHDisp", ""); @@ -796,13 +885,15 @@ parlistfcn.AddToList(&tasklistfcn); - parlistfcn.AddToList(&disp); - parlistfcn.AddToList(&hdisp); + parlistfcn.AddToList(&dispparams); + parlistfcn.AddToList(fHDispTrain); parlistfcn.AddToList(fCam); - parlistfcn.AddToList(fMatrixTrain); + parlistfcn.AddToList(fMatrixTrainCalc); + parlistfcn.AddToList(fMatrixTrainHists); //****************************** // entries in MTaskList - tasklistfcn.AddToList(&loop); + tasklistfcn.AddToList(&loopcalc); + tasklistfcn.AddToList(&loophists); tasklistfcn.AddToList(fDispCalcTrain); tasklistfcn.AddToList(&filldispplots); @@ -834,6 +925,6 @@ // get initial values of parameters - fVinit = disp.GetParameters(); - fStep = disp.GetStepsizes(); + fVinit = dispparams.GetParameters(); + fStep = dispparams.GetStepsizes(); TString name[fVinit.GetSize()]; @@ -901,5 +992,5 @@ TFile outparam(fFilenameParam, "RECREATE"); - disp.Write(); + dispparams.Write(); outparam.Close(); @@ -907,5 +998,5 @@ << fFilenameParam << "' :" << endl; - const TArrayD &check = disp.GetParameters(); + const TArrayD &check = dispparams.GetParameters(); for (Int_t i=0; iDraw(); *fLog << "End of Disp optimization part" << endl; @@ -932,14 +1023,14 @@ Bool_t MFindDisp::TestParams() { - if (fMatrixTest == NULL) - { - *fLog << "MFindDisp::TestParams; test matrix is not defined... aborting" + if (fMatrixTestCalc == NULL || fMatrixTestHists == NULL) + { + *fLog << "MFindDisp::TestParams; test matrices are not defined... aborting" << endl; return kFALSE; } - if (fMatrixTest->GetM().GetNrows() <= 0) - { - *fLog << "MFindDisp::TestParams; test matrix has no entries" + if (fMatrixTestCalc->GetM().GetNrows() <= 0 || fMatrixTestHists->GetM().GetNrows() <= 0) + { + *fLog << "MFindDisp::TestParams; test matrices have no entries" << endl; return kFALSE; @@ -955,6 +1046,6 @@ TFile inparam(fFilenameParam); - MDisp dispin; - dispin.Read("MDisp"); + MDispParameters dispin; + dispin.Read("MDispParameters"); inparam.Close(); @@ -973,22 +1064,22 @@ MTaskList tasklist2; - MDisp disp; - disp.SetParameters(dispPar); - - MMatrixLoop loop(fMatrixTest); - - // get the matrix pointer and mapping from MDispCalc object to set the same in MHDisp object - TArrayI map; - MHMatrix *tmpmatrix = NULL; - tmpmatrix = fDispCalcTest->GetMatrixMap(map); + MDispParameters dispparams; + 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 Chi^2 contributions + // fill the plots for Disp and sum up the Minimization parameter contributions MHDisp hdisp1("DispTest"); hdisp1.SetName("MHDispCorr"); hdisp1.SetSelectedPos(1); - hdisp1.SetMatrixMap(tmpmatrix,map); + hdisp1.SetMatrixMap(imatrix,imap); MFillH filldispplots1("MHDispCorr[MHDisp]", ""); @@ -996,5 +1087,5 @@ hdisp2.SetName("MHDispWrong"); hdisp2.SetSelectedPos(2); - hdisp2.SetMatrixMap(tmpmatrix,map); + hdisp2.SetMatrixMap(imatrix,imap); MFillH filldispplots2("MHDispWrong[MHDisp]", ""); @@ -1002,5 +1093,5 @@ hdisp3.SetName("MHDispM3Long"); hdisp3.SetSelectedPos(3); - hdisp3.SetMatrixMap(tmpmatrix,map); + hdisp3.SetMatrixMap(imatrix,imap); MFillH filldispplots3("MHDispM3Long[MHDisp]", ""); @@ -1008,5 +1099,5 @@ hdisp4.SetName("MHDispAsym"); hdisp4.SetSelectedPos(4); - hdisp4.SetMatrixMap(tmpmatrix,map); + hdisp4.SetMatrixMap(imatrix,imap); MFillH filldispplots4("MHDispAsym[MHDisp]", ""); @@ -1016,5 +1107,5 @@ parlist2.AddToList(&tasklist2); - parlist2.AddToList(&disp); + parlist2.AddToList(&dispparams); parlist2.AddToList(&hdisp1); parlist2.AddToList(&hdisp2); @@ -1022,10 +1113,12 @@ parlist2.AddToList(&hdisp4); parlist2.AddToList(fCam); - parlist2.AddToList(fMatrixTest); + parlist2.AddToList(fMatrixTestCalc); + parlist2.AddToList(fMatrixTestHists); //****************************** // entries in MTaskList - tasklist2.AddToList(&loop); + tasklist2.AddToList(&loopcalc); + tasklist2.AddToList(&loophists); tasklist2.AddToList(fDispCalcTest); tasklist2.AddToList(&filldispplots1); @@ -1064,10 +1157,2 @@ return kTRUE; } - - - - - - - - Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MFindDisp.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MFindDisp.h (revision 5879) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MFindDisp.h (revision 5904) @@ -16,4 +16,5 @@ class MHMatrix; class MDispCalc; +class MHDisp; class MFDisp; class MFilter; @@ -34,9 +35,14 @@ Int_t fHowManyTest; // number of events for testing - MHMatrix *fMatrixTrain; // training matrix - MHMatrix *fMatrixTest; // testing matrix + 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 @@ -93,6 +99,8 @@ Int_t iflag=0); - MHMatrix *GetMatrixTrain() { return fMatrixTrain; } - MHMatrix *GetMatrixTest() { return fMatrixTest; } + 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); Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MHDisp.cc =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MHDisp.cc (revision 5879) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MHDisp.cc (revision 5904) @@ -29,5 +29,5 @@ // // // container holding the histograms for Disp // -// also computes the Chi^2 of the Disp optimization // +// also computes the minimization parameter of the Disp optimization // // // ///////////////////////////////////////////////////////////////////////////// @@ -54,4 +54,5 @@ #include "MHMatrix.h" +#include "MParameters.h" #include "MParList.h" @@ -63,4 +64,8 @@ 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 // -------------------------------------------------------------------------- @@ -75,7 +80,10 @@ fTitle = title ? title : "Histograms for Disp"; - fSelectedPos = 1; // default MHDisp flag (selects Correct Disp srcpos) + 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); //-------------------------------------------------- @@ -110,10 +118,10 @@ fSkymapXY->SetYTitle("Y Disp source position [deg]"); - fHistChi2 = new TH1F("fHistChi2" , + fHistMinPar = new TH1F("fHistMinPar" , "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"); + fHistMinPar->SetDirectory(NULL); + fHistMinPar->UseCurrentStyle(); + fHistMinPar->SetXTitle("Minimization parameter = d^2 Disp to real srcpos [deg^2]"); + fHistMinPar->SetYTitle("# events"); fHistDuDv = new TH2F("fHistDuDv", @@ -125,17 +133,17 @@ 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]"); + 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", @@ -170,5 +178,5 @@ "Fraction events srcpos well assign vs. Size", 50, 2., 4., 0., 1.); fEvCorrAssign->SetDirectory(NULL); - fEvCorrAssign->UseCurrentStyle(); + fEvCorrAssign->SetStats(0); fEvCorrAssign->SetXTitle("log10(Size) [#phot]"); fEvCorrAssign->SetYTitle("Fraction of events with srcpos well assign"); @@ -186,8 +194,8 @@ delete fHistcosZA; delete fSkymapXY; - delete fHistChi2; + delete fHistMinPar; delete fHistDuDv; - delete fHistChi2Energy; - delete fHistChi2Size; + delete fHistMinParEnergy; + delete fHistMinParSize; delete fHistDuEnergy; delete fHistDuSize; @@ -204,9 +212,16 @@ Bool_t MHDisp::SetupFill(const MParList *pList) { - // reset all histograms and Chi^2 computing variables + // reset all histograms and Minimization parameter computing variables // before each new eventloop fNumEv = 0; - fSumChi2 = 0.; - fChi2 = 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(); @@ -214,8 +229,8 @@ fHistcosZA->Reset(); fSkymapXY->Reset(); - fHistChi2->Reset(); + fHistMinPar->Reset(); fHistDuDv->Reset(); - fHistChi2Energy->Reset(); - fHistChi2Size->Reset(); + fHistMinParEnergy->Reset(); + fHistMinParSize->Reset(); fHistDuEnergy->Reset(); fHistDuSize->Reset(); @@ -250,4 +265,6 @@ // if fMatrix exists, skip preprocessing and just read events from matrix; // if doesn't exist, read variables values from containers, so look for them + cout << "MHDisp::SetupFill: fMatrix = " << fMatrix << endl; + if (fMatrix) return kTRUE; @@ -291,8 +308,6 @@ if (!fPointing) { - *fLog << err << "MPointingPos not found... " - << "MMcEvt is going to be used to get Theta and Phi." - << endl; - // return kFALSE; + *fLog << err << "MPointingPos not found... aborting." << endl; + return kFALSE; } @@ -306,5 +321,5 @@ // // Set which selection algorithm for the Disp estimated source position -// we want to follow when filling the histograms: +// 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) @@ -323,5 +338,4 @@ // // Sets the Matrix and the array of mapped values for each Matrix column -// (see MDispCalc) // void MHDisp::SetMatrixMap(MHMatrix *matrix, TArrayI &map) @@ -334,4 +348,52 @@ // -------------------------------------------------------------------------- // +// 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 // @@ -349,5 +411,6 @@ // -------------------------------------------------------------------------- // -// Fill the histograms and sum up the Chi^2 outcome of each processed event +// 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) @@ -356,48 +419,30 @@ 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; - } + energy = fMatrix ? GetVal(kEnergy) : fMcEvt->GetEnergy(); + impact = fMatrix ? GetVal(kImpact) : fMcEvt->GetImpact(); + xmax = fMatrix ? GetVal(kLongitmax) : fMcEvt->GetLongitmax(); } 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(); + 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 width = width0 * fMm2Deg; - // Double_t length = length0 * fMm2Deg; Double_t xsrcpos = xsrcpos0 * fMm2Deg; Double_t ysrcpos = ysrcpos0 * fMm2Deg; @@ -453,5 +498,5 @@ // Distance between estimated Disp and real source position Double_t d2 = (xdisp-xsrcpos)*(xdisp-xsrcpos) + (ydisp-ysrcpos)*(ydisp-ysrcpos); - fHistChi2->Fill(d2); + fHistMinPar->Fill(d2); // Longitudinal and transversal distances between Disp and real source positon @@ -465,7 +510,7 @@ fHistcosZA->Fill(cos(theta/kRad2Deg)); - // to check the size dependence of the optimization - fHistChi2Energy->Fill(log10(energy),d2); - fHistChi2Size->Fill(log10(size),d2); + // to check the size and energy dependence of the optimization + fHistMinParEnergy->Fill(log10(energy),d2); + fHistMinParSize->Fill(log10(size),d2); fHistDuEnergy->Fill(log10(energy),Du); fHistDuSize->Fill(log10(size),Du); @@ -473,7 +518,7 @@ fHistDvSize->Fill(log10(size),Dv); - // variables for Chi^2 calculation (= d^2) + // variables for the Minimization parameter calculation (= d^2) fNumEv += 1; - fSumChi2 += d2; + fSumMinPar += d2; return kTRUE; @@ -483,15 +528,15 @@ // -------------------------------------------------------------------------- // -// Calculates the final Chi^2 of the Disp optimization +// Calculates the final Minimization parameter of the Disp optimization // Bool_t MHDisp::Finalize() { - fChi2 = fSumChi2/fNumEv; + fMinPar->SetVal(fSumMinPar/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; + *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(); @@ -542,27 +587,29 @@ 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]"); + 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); - 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]"); + 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); - profChi2Size->SetTitleOffset(1.2,"Y"); - profChi2Size->DrawClone(opt); - profChi2Size->SetBit(kCanDelete); + profMinParSize->SetTitleOffset(1.2,"Y"); + profMinParSize->SetStats(0); + profMinParSize->DrawClone(opt); + profMinParSize->SetBit(kCanDelete); gPad->Modified(); @@ -604,4 +651,5 @@ gPad->SetBorderMode(0); profDuEnergy->SetTitleOffset(1.2,"Y"); + profDuEnergy->SetStats(0); profDuEnergy->DrawClone(opt); profDuEnergy->SetBit(kCanDelete); @@ -614,4 +662,5 @@ gPad->SetBorderMode(0); profDuSize->SetTitleOffset(1.2,"Y"); + profDuSize->SetStats(0); profDuSize->DrawClone(opt); profDuSize->SetBit(kCanDelete); @@ -634,4 +683,5 @@ gPad->SetBorderMode(0); profDvEnergy->SetTitleOffset(1.2,"Y"); + profDvEnergy->SetStats(0); profDvEnergy->DrawClone(opt); profDvEnergy->SetBit(kCanDelete); @@ -644,4 +694,5 @@ gPad->SetBorderMode(0); profDvSize->SetTitleOffset(1.2,"Y"); + profDvSize->SetStats(0); profDvSize->DrawClone(opt); profDvSize->SetBit(kCanDelete); Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MHDisp.h =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/MHDisp.h (revision 5879) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/MHDisp.h (revision 5904) @@ -21,4 +21,5 @@ class MPointingPos; class MHMatrix; +class MParameterD; class MHDisp : public MH @@ -36,6 +37,6 @@ TString fImageParDispName; - Int_t fSelectedPos; // flag to select which Disp estimated position - // we want to fill the histograms (see Set function) + 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 @@ -43,10 +44,10 @@ 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) + 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 *fHistChi2Energy; // Chi^2 (= d^2) vs. Energy - TH2F *fHistChi2Size; // Chi^2 (= d^2) vs. Size + 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 @@ -58,13 +59,13 @@ 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 + 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 fSumChi2; // current sum of Chi^2 - Double_t fChi2; // final Chi^2 of the Disp optimization + Double_t fSumMinPar; // current sum of the minimization parameter + MParameterD *fMinPar; // final minimization parameters of the Disp optimization public: @@ -81,6 +82,8 @@ void SetMatrixMap(MHMatrix *matrix, TArrayI &map); + void GetMatrixMap(MHMatrix* &matrix, TArrayI &map); // get matrix and its mapping array - Double_t GetChi2() { return fChi2; } + void InitMapping(MHMatrix *mat); // define the matrix of variables + // needed for the Disp optimization TH1F *GetHistEnergy() { return fHistEnergy; } @@ -88,8 +91,8 @@ TH1F *GetHistcosZA() { return fHistcosZA; } TH2F *GetSkymapXY() { return fSkymapXY; } - TH1F *GetHistChi2() { return fHistChi2; } + TH1F *GetHistMinPar() { return fHistMinPar; } TH2F *GetHistDuDv() { return fHistDuDv; } - TH2F *GetHistChi2Energy() { return fHistChi2Energy; } - TH2F *GetHistChi2Size() { return fHistChi2Size; } + TH2F *GetHistMinParEnergy() { return fHistMinParEnergy; } + TH2F *GetHistMinParSize() { return fHistMinParSize; } TH2F *GetHistDuEnergy() { return fHistDuEnergy; } TH2F *GetHistDuSize() { return fHistDuSize; } @@ -100,5 +103,6 @@ void Draw(Option_t *opt=""); - ClassDef(MHDisp, 1) // Container holding the Histograms for Disp + ClassDef(MHDisp, 1) // Container holding the Histograms for Disp and + // the parameter to minimize in the Disp optimization }; Index: trunk/MagicSoft/Mars/mtemp/mifae/library/Makefile =================================================================== --- trunk/MagicSoft/Mars/mtemp/mifae/library/Makefile (revision 5879) +++ trunk/MagicSoft/Mars/mtemp/mifae/library/Makefile (revision 5904) @@ -74,5 +74,5 @@ MTopologyCalc.cc \ MImageParDisp.cc \ - MDisp.cc \ + MDispParameters.cc \ MDispCalc.cc \ MHDisp.cc \