Index: trunk/MagicSoft/Mars/manalysis/MCT1FindSupercuts.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCT1FindSupercuts.cc (revision 2350) +++ trunk/MagicSoft/Mars/manalysis/MCT1FindSupercuts.cc (revision 2357) @@ -57,5 +57,6 @@ #include "MFEventSelector2.h" #include "MFillH.h" -#include "MGeomCamCT1Daniel.h" +//#include "MGeomCamCT1Daniel.h" +#include "MGeomCamCT1.h" #include "MH3.h" #include "MHCT1Supercuts.h" @@ -95,4 +96,6 @@ Double_t *par, Int_t iflag) { + //cout << "entry fcnSupercuts" << endl; + //------------------------------------------------------------- // save pointer to the MINUIT object for optimizing the supercuts @@ -106,4 +109,5 @@ MParList *plistfcn = (MParList*)evtloopfcn->GetParList(); + //MTaskList *tasklistfcn = (MTaskList*)plistfcn->FindObject("MTaskList"); MCT1Supercuts *super = (MCT1Supercuts*)plistfcn->FindObject("MCT1Supercuts"); @@ -130,4 +134,6 @@ // for testing //TArrayD checkparameters = super->GetParameters(); + //gLog << "fcnsupercuts : fNpari, fNparx =" << fNpari << ", " + // << fNparx << endl; //gLog << "fcnsupercuts : i, par, checkparameters =" << endl; //for (Int_t i=0; iEventloop(); + //tasklistfcn->PrintStatistics(0, kTRUE); + MH3* alpha = (MH3*)plistfcn->FindObject("AlphaFcn", "MH3"); if (!alpha) @@ -147,5 +155,5 @@ TH1 &alphaHist = alpha->GetHist(); - //alphaHist.SetXTitle("|\\alpha| [\\circ]"); + alphaHist.SetName("alpha-fcnSupercuts"); //------------------------------------------- @@ -161,5 +169,5 @@ const Double_t alphamin = 30.0; const Double_t alphamax = 90.0; - const Int_t degree = 4; + const Int_t degree = 2; Bool_t drawpoly; @@ -240,7 +248,8 @@ //--------------------------- // camera geometry is needed for conversion mm ==> degree - fCam = new MGeomCamCT1Daniel; - - // matrices to contain the the training/test samples + //fCam = new MGeomCamCT1Daniel; + fCam = new MGeomCamCT1; + + // matrices to contain the training/test samples fMatrixTrain = new MHMatrix("MatrixTrain"); fMatrixTest = new MHMatrix("MatrixTest"); @@ -277,5 +286,6 @@ // Bool_t MCT1FindSupercuts::DefineTrainMatrix(const TString &nametrain, - const Int_t howmanytrain, MH3 &hreftrain) + const Int_t howmanytrain, MH3 &hreftrain, + const TString &filetrain) { if (nametrain.IsNull() || howmanytrain <= 0) @@ -343,4 +353,18 @@ *fLog << "=============================================" << endl; + //-------------------------- + // write out training matrix + + if (filetrain != "") + { + TFile filetr(filetrain, "RECREATE"); + fMatrixTrain->Write(); + filetr.Close(); + + *fLog << "MCT1FindSupercuts::DefineTrainMatrix; Training matrix was written onto file '" + << filetrain << "'" << endl; + } + + return kTRUE; } @@ -356,5 +380,6 @@ // Bool_t MCT1FindSupercuts::DefineTestMatrix(const TString &nametest, - const Int_t howmanytest, MH3 &hreftest) + const Int_t howmanytest, MH3 &hreftest, + const TString &filetest) { if (nametest.IsNull() || howmanytest<=0) @@ -422,4 +447,18 @@ *fLog << "=============================================" << endl; + //-------------------------- + // write out test matrix + + if (filetest != "") + { + TFile filete(filetest, "RECREATE", ""); + fMatrixTest->Write(); + filete.Close(); + + *fLog << "MCT1FindSupercuts::DefineTestMatrix; Test matrix was written onto file '" + << filetest << "'" << endl; + } + + return kTRUE; } @@ -434,5 +473,6 @@ const TString &name, const Int_t howmanytrain, MH3 &hreftrain, - const Int_t howmanytest, MH3 &hreftest) + const Int_t howmanytest, MH3 &hreftest, + const TString &filetrain, const TString &filetest) { *fLog << "=============================================" << endl; @@ -511,5 +551,6 @@ evtloop.SetProgressBar(&bar); - if (!evtloop.Eventloop()) + Int_t maxev = -1; + if (!evtloop.Eventloop(maxev)) return kFALSE; @@ -520,5 +561,5 @@ if (TMath::Abs(generatedtrain-howmanytrain) > TMath::Sqrt(9.*howmanytrain)) { - *fLog << "MCT1FindSupercuts::DefineTrainMatrix; no.of generated events (" + *fLog << "MCT1FindSupercuts::DefineTrainTestMatrix; no.of generated events (" << generatedtrain << ") is incompatible with the no.of requested events (" @@ -530,5 +571,5 @@ if (TMath::Abs(generatedtest-howmanytest) > TMath::Sqrt(9.*howmanytest)) { - *fLog << "MCT1FindSupercuts::DefineTestMatrix; no.of generated events (" + *fLog << "MCT1FindSupercuts::DefineTrainTestMatrix; no.of generated events (" << generatedtest << ") is incompatible with the no.of requested events (" @@ -541,4 +582,30 @@ + //-------------------------- + // write out training matrix + + if (filetrain != "") + { + TFile filetr(filetrain, "RECREATE"); + fMatrixTrain->Write(); + filetr.Close(); + + *fLog << "MCT1FindSupercuts::DefineTrainTestMatrix; Training matrix was written onto file '" + << filetrain << "'" << endl; + } + + //-------------------------- + // write out test matrix + + if (filetest != "") + { + TFile filete(filetest, "RECREATE", ""); + fMatrixTest->Write(); + filete.Close(); + + *fLog << "MCT1FindSupercuts::DefineTrainTestMatrix; Test matrix was written onto file '" + << filetest << "'" << endl; + } + return kTRUE; } @@ -548,6 +615,4 @@ // Read training and test matrices from files // -// -// $$$$$$$$$$$$$$ this does not work !!! ??? $$$$$$$$$$$$$$$$$$$$$$ // @@ -563,4 +628,5 @@ *fLog << "MCT1FindSupercuts::ReadMatrix; Training matrix was read in from file '" << filetrain << "'" << endl; + filetr.Close(); @@ -574,45 +640,9 @@ *fLog << "MCT1FindSupercuts::ReadMatrix; Test matrix was read in from file '" << filetest << "'" << endl; - + filete.Close(); return kTRUE; } -// -------------------------------------------------------------------------- -// -// Write training and test matrices onto files -// -// -// $$$$$$$$$$$$$$ this does not work !!! ??? $$$$$$$$$$$$$$$$$$$$$$ -// -// -Bool_t MCT1FindSupercuts::WriteMatrix(const TString &filetrain, const TString &filetest) -{ - //-------------------------- - // write out training matrix - - TFile filetr(filetrain, "RECREATE"); - - *fLog << "nach TFile" << endl; - - fMatrixTrain->Write(); - - *fLog << "MCT1FindSupercuts::WriteMatrix; Training matrix was written onto file '" - << filetrain << "'" << endl; - - - //-------------------------- - // write out test matrix - - TFile filete(filetest, "RECREATE"); - fMatrixTest->Print("SizeCols"); - fMatrixTest->Write(); - - *fLog << "MCT1FindSupercuts::WriteMatrix; Test matrix was written onto file '" - << filetest << "'" << endl; - - - return kTRUE; -} //------------------------------------------------------------------------ @@ -645,13 +675,18 @@ // put the supercuts hadronness // -// - for the minimization, the starting values of the parameters are taken as -// MCT1Supercuts::GetParameters(fVinit) +// - for the minimization, the starting values of the parameters are taken +// - from file parSCinit (if it is != "") +// - or from the arrays params and/or steps // //---------------------------------------------------------------------- -Bool_t MCT1FindSupercuts::FindParams() +Bool_t MCT1FindSupercuts::FindParams(TString parSCinit, + TArrayD ¶ms, TArrayD &steps) { // Setup the event loop which will be executed in the // fcnSupercuts function of MINUIT // + // parSCinit is the name of the file containing the initial values + // of the parameters; + // if parSCinit = "" 'params' and 'steps' are taken as initial values *fLog << "=============================================" << endl; @@ -685,8 +720,43 @@ MMatrixLoop loop(fMatrixTrain); + //-------------------------------- // create container for the supercut parameters - // and set them to their default values + // and set them to their initial values MCT1Supercuts super; + // take initial values from file parSCinit + if (parSCinit != "") + { + TFile inparam(parSCinit); + super.Read("MCT1Supercuts"); + inparam.Close(); + *fLog << "MCT1FindSupercuts::FindParams; initial values of parameters are taken from file " + << parSCinit << endl; + } + + // take initial values from 'params' and/or 'steps' + else if (params.GetSize() != 0 || steps.GetSize() != 0 ) + { + if (params.GetSize() != 0) + { + *fLog << "MCT1FindSupercuts::FindParams; initial values of parameters are taken from 'params'" + << endl; + super.SetParameters(params); + } + if (steps.GetSize() != 0) + { + *fLog << "MCT1FindSupercuts::FindParams; initial step sizes are taken from 'steps'" + << endl; + super.SetStepsizes(steps); + } + } + else + { + *fLog << "MCT1FindSupercuts::FindParams; initial values and step sizes are taken from the MCT1Supercits constructor" + << endl; + } + + + //-------------------------------- // calculate supercuts hadronness fCalcHadTrain->SetHadronnessName(fHadronnessName); @@ -761,6 +831,7 @@ // - // get initial values of parameters from MCT1Supercuts + // get initial values of parameters fVinit = super.GetParameters(); + fStep = super.GetStepsizes(); TString name[fVinit.GetSize()]; @@ -776,16 +847,45 @@ name[i] = "p"; name[i] += i+1; - fStep[i] = TMath::Abs(fVinit[i]/10.0); + //fStep[i] = TMath::Abs(fVinit[i]/10.0); fLimlo[i] = -100.0; fLimup[i] = 100.0; fFix[i] = 0; - - // vary only first 48 parameters - if (i >= 48) - { - fStep[i] = 0.0; - fFix[i] = 1; - } - } + } + + // these parameters make no sense, fix them at 0.0 + fVinit[33] = 0.0; + fStep[33] = 0.0; + fFix[33] = 1; + + fVinit[36] = 0.0; + fStep[36] = 0.0; + fFix[36] = 1; + + fVinit[39] = 0.0; + fStep[39] = 0.0; + fFix[39] = 1; + + fVinit[41] = 0.0; + fStep[41] = 0.0; + fFix[41] = 1; + + fVinit[44] = 0.0; + fStep[44] = 0.0; + fFix[44] = 1; + + fVinit[47] = 0.0; + fStep[47] = 0.0; + fFix[47] = 1; + + // vary only first 48 parameters + //for (UInt_t i=0; i= 48) + // { + // fStep[i] = 0.0; + // fFix[i] = 1; + // } + //} + // ------------------------------------------- @@ -910,4 +1010,8 @@ MH3 alphabefore("MatrixTest[7]"); alphabefore.SetName(mh3NameB); + + TH1 &alphahistb = alphabefore.GetHist(); + alphahistb.SetName("alphaBefore-TestParams"); + MFillH fillalphabefore(&alphabefore); fillalphabefore.SetName("FillAlphaBefore"); @@ -924,4 +1028,8 @@ MH3 alphaafter("MatrixTest[7]"); alphaafter.SetName(mh3NameA); + + TH1 &alphahista = alphaafter.GetHist(); + alphahista.SetName("alphaAfter-TestParams"); + MFillH fillalphaafter(&alphaafter); fillalphaafter.SetName("FillAlphaAfter"); @@ -978,4 +1086,5 @@ TH1 &alphaHist2 = alphaAfter->GetHist(); alphaHist2.SetXTitle("|\\alpha| [\\circ]"); + alphaHist2.SetName("alpha-TestParams)"); TCanvas *c = new TCanvas("AlphaAfterSC", "AlphaTest", 600, 300); @@ -999,5 +1108,5 @@ const Double_t alphamin = 30.0; const Double_t alphamax = 90.0; - const Int_t degree = 4; + const Int_t degree = 2; const Bool_t drawpoly = kTRUE; const Bool_t fitgauss = kFALSE; @@ -1025,2 +1134,3 @@ return kTRUE; } + Index: trunk/MagicSoft/Mars/manalysis/MCT1FindSupercuts.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCT1FindSupercuts.h (revision 2350) +++ trunk/MagicSoft/Mars/manalysis/MCT1FindSupercuts.h (revision 2357) @@ -89,17 +89,20 @@ Bool_t DefineTrainMatrix(const TString &name, const Int_t howmany, - MH3 &href); + MH3 &href, const TString &filetrain); Bool_t DefineTestMatrix(const TString &name, const Int_t howmany, - MH3 &href); + MH3 &href, const TString &filetest); Bool_t DefineTrainTestMatrix(const TString &name, - const Int_t howmanytrain, MH3 &hreftrain, - const Int_t howmanytest, MH3 &hreftest); + const Int_t howmanytrain, MH3 &hreftrain, + const Int_t howmanytest, MH3 &hreftest, + const TString &filetrain, const TString &filetest); + + MHMatrix *GetMatrixTrain() { return fMatrixTrain; } + MHMatrix *GetMatrixTest() { return fMatrixTest; } Bool_t ReadMatrix( const TString &filetrain, const TString &filetest); - Bool_t WriteMatrix(const TString &filetrain, const TString &filetest); - Bool_t FindParams(); + Bool_t FindParams(TString parSCinit, TArrayD ¶ms, TArrayD &steps); Bool_t TestParams(); Index: trunk/MagicSoft/Mars/manalysis/MCT1Supercuts.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCT1Supercuts.cc (revision 2350) +++ trunk/MagicSoft/Mars/manalysis/MCT1Supercuts.cc (revision 2357) @@ -50,10 +50,12 @@ // MCT1Supercuts::MCT1Supercuts(const char *name, const char *title) - : fParameters(72), - fLengthUp(fParameters.GetArray()), fLengthLo(fParameters.GetArray()+8), + : fParameters(104), fStepsizes(104), + fLengthUp(fParameters.GetArray()), fLengthLo(fParameters.GetArray()+8), fWidthUp(fParameters.GetArray()+16), fWidthLo(fParameters.GetArray()+24), - fDistUp(fParameters.GetArray()+32), fDistLo(fParameters.GetArray()+40), - fAsymUp(fParameters.GetArray()+48), fAsymLo(fParameters.GetArray()+56), - fAlphaUp(fParameters.GetArray()+64) + fDistUp(fParameters.GetArray()+32), fDistLo(fParameters.GetArray()+40), + fAsymUp(fParameters.GetArray()+48), fAsymLo(fParameters.GetArray()+56), + fConcUp(fParameters.GetArray()+64), fConcLo(fParameters.GetArray()+72), + fLeakage1Up(fParameters.GetArray()+80), fLeakage1Lo(fParameters.GetArray()+88), + fAlphaUp(fParameters.GetArray()+96) { fName = name ? name : "MCT1Supercuts"; @@ -71,4 +73,7 @@ void MCT1Supercuts::InitParameters() { + //--------------------------------------------------- + // these are Daniel's original values for Mkn 421 + fLengthUp[0] = 0.315585; fLengthUp[1] = 0.001455; @@ -80,4 +85,13 @@ fLengthUp[7] = -0.013463; + fLengthLo[0] = 0.151530; + fLengthLo[1] = 0.028323; + fLengthLo[2] = 0.510707; + fLengthLo[3] = 0.053089; + fLengthLo[4] = 0.013708; + fLengthLo[5] = 2.357993; + fLengthLo[6] = 0.000080; + fLengthLo[7] = -0.007157; + fWidthUp[0] = 0.145412; fWidthUp[1] = -0.001771; @@ -89,4 +103,13 @@ fWidthUp[7] = -0.016703; + fWidthLo[0] = 0.089187; + fWidthLo[1] = -0.006430; + fWidthLo[2] = 0.074442; + fWidthLo[3] = 0.003738; + fWidthLo[4] = -0.004256; + fWidthLo[5] = -0.014101; + fWidthLo[6] = 0.006126; + fWidthLo[7] = -0.002849; + fDistUp[0] = 1.787943; fDistUp[1] = 0; @@ -98,22 +121,4 @@ fDistUp[7] = 0; - fLengthLo[0] = 0.151530; - fLengthLo[1] = 0.028323; - fLengthLo[2] = 0.510707; - fLengthLo[3] = 0.053089; - fLengthLo[4] = 0.013708; - fLengthLo[5] = 2.357993; - fLengthLo[6] = 0.000080; - fLengthLo[7] = -0.007157; - - fWidthLo[0] = 0.089187; - fWidthLo[1] = -0.006430; - fWidthLo[2] = 0.074442; - fWidthLo[3] = 0.003738; - fWidthLo[4] = -0.004256; - fWidthLo[5] = -0.014101; - fWidthLo[6] = 0.006126; - fWidthLo[7] = -0.002849; - fDistLo[0] = 0.589406; fDistLo[1] = 0; @@ -124,22 +129,61 @@ fDistLo[6] = -0.001750; fDistLo[7] = 0; - - fAsymUp[0] = 0.061267; - fAsymUp[1] = 0.014462; - fAsymUp[2] = 0.014327; - fAsymUp[3] = 0.014540; - fAsymUp[4] = 0.013391; - fAsymUp[5] = 0.012319; - fAsymUp[6] = 0.010444; - fAsymUp[7] = 0.008328; - - fAsymLo[0] = -0.012055; - fAsymLo[1] = 0.009157; - fAsymLo[2] = 0.005441; - fAsymLo[3] = 0.000399; - fAsymLo[4] = 0.001433; - fAsymLo[5] = -0.002050; - fAsymLo[6] = -0.000104; - fAsymLo[7] = -0.001188; + + + // dummy values + + fAsymUp[0] = 1.e10; + fAsymUp[1] = 0.0; + fAsymUp[2] = 0.0; + fAsymUp[3] = 0.0; + fAsymUp[4] = 0.0; + fAsymUp[5] = 0.0; + fAsymUp[6] = 0.0; + fAsymUp[7] = 0.0; + + fAsymLo[0] = -1.e10; + fAsymLo[1] = 0.0; + fAsymLo[2] = 0.0; + fAsymLo[3] = 0.0; + fAsymLo[4] = 0.0; + fAsymLo[5] = 0.0; + fAsymLo[6] = 0.0; + fAsymLo[7] = 0.0; + + fConcUp[0] = 1.e10; + fConcUp[1] = 0.0; + fConcUp[2] = 0.0; + fConcUp[3] = 0.0; + fConcUp[4] = 0.0; + fConcUp[5] = 0.0; + fConcUp[6] = 0.0; + fConcUp[7] = 0.0; + + fConcLo[0] = -1.e10; + fConcLo[1] = 0.0; + fConcLo[2] = 0.0; + fConcLo[3] = 0.0; + fConcLo[4] = 0.0; + fConcLo[5] = 0.0; + fConcLo[6] = 0.0; + fConcLo[7] = 0.0; + + fLeakage1Up[0] = 1.e10; + fLeakage1Up[1] = 0.0; + fLeakage1Up[2] = 0.0; + fLeakage1Up[3] = 0.0; + fLeakage1Up[4] = 0.0; + fLeakage1Up[5] = 0.0; + fLeakage1Up[6] = 0.0; + fLeakage1Up[7] = 0.0; + + fLeakage1Lo[0] = -1.e10; + fLeakage1Lo[1] = 0.0; + fLeakage1Lo[2] = 0.0; + fLeakage1Lo[3] = 0.0; + fLeakage1Lo[4] = 0.0; + fLeakage1Lo[5] = 0.0; + fLeakage1Lo[6] = 0.0; + fLeakage1Lo[7] = 0.0; fAlphaUp[0] = 13.123440; @@ -151,4 +195,139 @@ fAlphaUp[6] = 0; fAlphaUp[7] = 0; + + //--------------------------------------------------- + // fStepsizes + // if == 0.0 the parameter will be fixed in the minimization + // != 0.0 initial step sizes for the parameters + + // LengthUp + fStepsizes[0] = 0.03; + fStepsizes[1] = 0.0002; + fStepsizes[2] = 0.02; + fStepsizes[3] = 0.0006; + fStepsizes[4] = 0.0002; + fStepsizes[5] = 0.002; + fStepsizes[6] = 0.0008; + fStepsizes[7] = 0.002; + + // LengthLo + fStepsizes[8] = 0.02; + fStepsizes[9] = 0.003; + fStepsizes[10] = 0.05; + fStepsizes[11] = 0.006; + fStepsizes[12] = 0.002; + fStepsizes[13] = 0.3; + fStepsizes[14] = 0.0001; + fStepsizes[15] = 0.0008; + + // WidthUp + fStepsizes[16] = 0.02; + fStepsizes[17] = 0.0002; + fStepsizes[18] = 0.006; + fStepsizes[19] = 0.003; + fStepsizes[20] = 0.002; + fStepsizes[21] = 0.006; + fStepsizes[22] = 0.002; + fStepsizes[23] = 0.002; + + // WidthLo + fStepsizes[24] = 0.009; + fStepsizes[25] = 0.0007; + fStepsizes[26] = 0.008; + fStepsizes[27] = 0.0004; + fStepsizes[28] = 0.0005; + fStepsizes[29] = 0.002; + fStepsizes[30] = 0.0007; + fStepsizes[31] = 0.003; + + // DistUp + fStepsizes[32] = 0.2; + fStepsizes[33] = 0.0; + fStepsizes[34] = 0.3; + fStepsizes[35] = 0.02; + fStepsizes[36] = 0.0; + fStepsizes[37] = 0.03; + fStepsizes[38] = 0.02; + fStepsizes[39] = 0.0; + + // DistLo + fStepsizes[40] = 0.06; + fStepsizes[41] = 0.0; + fStepsizes[42] = 0.009; + fStepsizes[43] = 0.0008; + fStepsizes[44] = 0.0; + fStepsizes[45] = 0.005; + fStepsizes[46] = 0.0002; + fStepsizes[47] = 0.0; + + // AsymUp + fStepsizes[48] = 0.0; + fStepsizes[49] = 0.0; + fStepsizes[50] = 0.0; + fStepsizes[51] = 0.0; + fStepsizes[52] = 0.0; + fStepsizes[53] = 0.0; + fStepsizes[54] = 0.0; + fStepsizes[55] = 0.0; + + // AsymLo + fStepsizes[56] = 0.0; + fStepsizes[57] = 0.0; + fStepsizes[58] = 0.0; + fStepsizes[59] = 0.0; + fStepsizes[60] = 0.0; + fStepsizes[61] = 0.0; + fStepsizes[62] = 0.0; + fStepsizes[63] = 0.0; + + // ConcUp + fStepsizes[64] = 0.0; + fStepsizes[65] = 0.0; + fStepsizes[66] = 0.0; + fStepsizes[67] = 0.0; + fStepsizes[68] = 0.0; + fStepsizes[69] = 0.0; + fStepsizes[70] = 0.0; + fStepsizes[71] = 0.0; + + // ConcLo + fStepsizes[72] = 0.0; + fStepsizes[73] = 0.0; + fStepsizes[74] = 0.0; + fStepsizes[75] = 0.0; + fStepsizes[76] = 0.0; + fStepsizes[77] = 0.0; + fStepsizes[78] = 0.0; + fStepsizes[79] = 0.0; + + // Leakage1Up + fStepsizes[80] = 0.0; + fStepsizes[81] = 0.0; + fStepsizes[82] = 0.0; + fStepsizes[83] = 0.0; + fStepsizes[84] = 0.0; + fStepsizes[85] = 0.0; + fStepsizes[86] = 0.0; + fStepsizes[87] = 0.0; + + // Leakage1Lo + fStepsizes[88] = 0.0; + fStepsizes[89] = 0.0; + fStepsizes[90] = 0.0; + fStepsizes[91] = 0.0; + fStepsizes[92] = 0.0; + fStepsizes[93] = 0.0; + fStepsizes[94] = 0.0; + fStepsizes[95] = 0.0; + + // AlphaUp + fStepsizes[96] = 0.0; + fStepsizes[97] = 0.0; + fStepsizes[98] = 0.0; + fStepsizes[99] = 0.0; + fStepsizes[100] = 0.0; + fStepsizes[101] = 0.0; + fStepsizes[102] = 0.0; + fStepsizes[103] = 0.0; } @@ -173,13 +352,37 @@ } - - - - - - - - - - - +// -------------------------------------------------------------------------- +// +// Set the step sizes from the array 'd' +// +// +Bool_t MCT1Supercuts::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/manalysis/MCT1Supercuts.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCT1Supercuts.h (revision 2350) +++ trunk/MagicSoft/Mars/manalysis/MCT1Supercuts.h (revision 2357) @@ -14,4 +14,5 @@ private: TArrayD fParameters; // supercut parameters + TArrayD fStepsizes; // step sizes of supercut parameters Double_t *fLengthUp; //! @@ -23,5 +24,12 @@ Double_t *fAsymUp; //! Double_t *fAsymLo; //! + + Double_t *fConcUp; //! + Double_t *fConcLo; //! + Double_t *fLeakage1Up; //! + Double_t *fLeakage1Lo; //! + Double_t *fAlphaUp; //! + public: @@ -31,5 +39,8 @@ Bool_t SetParameters(const TArrayD &d); + Bool_t SetStepsizes(const TArrayD &d); + const TArrayD &GetParameters() const { return fParameters; } + const TArrayD &GetStepsizes() const { return fStepsizes; } const Double_t *GetLengthUp() const { return fLengthUp; } @@ -41,4 +52,11 @@ const Double_t *GetAsymUp() const { return fAsymUp; } const Double_t *GetAsymLo() const { return fAsymLo; } + + const Double_t *GetConcUp() const { return fConcUp; } + const Double_t *GetConcLo() const { return fConcLo; } + + const Double_t *GetLeakage1Up() const { return fLeakage1Up; } + const Double_t *GetLeakage1Lo() const { return fLeakage1Lo; } + const Double_t *GetAlphaUp() const { return fAlphaUp; } @@ -47,2 +65,5 @@ #endif + + + Index: trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.cc (revision 2350) +++ trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.cc (revision 2357) @@ -44,4 +44,5 @@ #include "MHillasExt.h" #include "MHillasSrc.h" +#include "MNewImagePar.h" #include "MMcEvt.hxx" #include "MCerPhotEvt.h" @@ -65,7 +66,9 @@ MCT1SupercutsCalc::MCT1SupercutsCalc(const char *hilname, - const char *hilsrcname, const char *name, const char *title) + const char *hilsrcname, + const char *name, const char *title) : fHadronnessName("MHadronness"), fHilName(hilname), fHilSrcName(hilsrcname), - fSuperName("MCT1Supercuts") + fHilExtName("MHillasExt"), fNewParName("MNewImagePar"), + fSuperName("MCT1Supercuts") { fName = name ? name : "MCT1SupercutsCalc"; @@ -101,5 +104,4 @@ return kFALSE; } - if (fMatrix) @@ -114,8 +116,22 @@ } + fHilExt = (MHillasExt*)pList->FindObject(fHilExtName, "MHillasExt"); + if (!fHilExt) + { + *fLog << err << fHilExtName << " [MHillasExt] not found... aborting." << endl; + return kFALSE; + } + fHilSrc = (MHillasSrc*)pList->FindObject(fHilSrcName, "MHillasSrc"); if (!fHilSrc) { *fLog << err << fHilSrcName << " [MHillasSrc] not found... aborting." << endl; + return kFALSE; + } + + fNewPar = (MNewImagePar*)pList->FindObject(fNewParName, "MNewImagePar"); + if (!fNewPar) + { + *fLog << err << fNewParName << " [MNewImagePar] not found... aborting." << endl; return kFALSE; } @@ -172,5 +188,13 @@ Double_t MCT1SupercutsCalc::GetVal(Int_t i) const { - return (*fMatrix)[fMap[i]]; + + Double_t val = (*fMatrix)[fMap[i]]; + + + //*fLog << "MCT1SupercutsCalc::GetVal; i, fMatrix, fMap, val = " + // << i << ", " << fMatrix << ", " << fMap[i] << ", " << val << endl; + + + return val; } @@ -187,5 +211,5 @@ { if (fMatrix) - return; + return; fMatrix = mat; @@ -199,10 +223,12 @@ fMap[6] = fMatrix->AddColumn("MHillasSrc.fDist"); fMap[7] = fMatrix->AddColumn("fabs(MHillasSrc.fAlpha)"); + fMap[8] = fMatrix->AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillasExt.fM3Long)"); + fMap[9] = fMatrix->AddColumn("MNewImagePar.fConc"); + fMap[10]= fMatrix->AddColumn("MNewImagePar.fLeakage1"); } // --------------------------------------------------------------------------- // -// Evaluate dynamical supercuts for CT1 Mkn421 2001 (Daniel Kranich) -// optimized for mkn 421 2001 data +// Evaluate dynamical supercuts // // set hadronness to 0.25 if cuts are fullfilled @@ -222,4 +248,12 @@ const Double_t dist0 = fMatrix ? GetVal(6) : fHilSrc->GetDist(); + Double_t help; + if (!fMatrix) + help = TMath::Sign(fHilExt->GetM3Long(), + fHilSrc->GetCosDeltaAlpha()); + const Double_t asym0 = fMatrix ? GetVal(8) : help; + const Double_t conc = fMatrix ? GetVal(9) : fNewPar->GetConc(); + const Double_t leakage = fMatrix ? GetVal(10): fNewPar->GetLeakage1(); + const Double_t newdist = dist0 * fMm2Deg; @@ -236,20 +270,67 @@ const Double_t length = length0 * fMm2Deg; const Double_t width = width0 * fMm2Deg; - - if (newdist < 1.05 && + const Double_t asym = asym0 * fMm2Deg; + + /* + *fLog << "newdist, length, width, asym, dist, conc, leakage = " + << newdist << ", " << length << ", " << width << ", " + << asym << ", " << dist << ", " << conc << ", " << leakage + << endl; + + *fLog << "upper cuts in newdist, length, width, asym, dist, conc, leakage = " + << CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) << ", " + + << CtsMCut (fSuper->GetLengthUp(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetLengthLo(), dmls, dmcza, dmls2, dd2) << ", " + + << CtsMCut (fSuper->GetWidthUp(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetWidthLo(), dmls, dmcza, dmls2, dd2) << ", " + + << CtsMCut (fSuper->GetAsymUp(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetAsymLo(), dmls, dmcza, dmls2, dd2) << ", " + + << CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) << ", " + + << CtsMCut (fSuper->GetConcUp(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetConcLo(), dmls, dmcza, dmls2, dd2) << ", " + + << CtsMCut (fSuper->GetLeakage1Up(), dmls, dmcza, dmls2, dd2) << ", " + << CtsMCut (fSuper->GetLeakage1Lo(), dmls, dmcza, dmls2, dd2) << ", " + << endl; + */ + + + if ( + //dist < 1.05 && + //newdist < 1.05 && + newdist < CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2) && newdist > CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) && - dist < 1.05 && + length < CtsMCut (fSuper->GetLengthUp(), dmls, dmcza, dmls2, dd2) && length > CtsMCut (fSuper->GetLengthLo(), dmls, dmcza, dmls2, dd2) && + width < CtsMCut (fSuper->GetWidthUp(), dmls, dmcza, dmls2, dd2) && width > CtsMCut (fSuper->GetWidthLo(), dmls, dmcza, dmls2, dd2) && - //asym < CtsMCut (fSuper->GetAsymUp(), dmls, dmcza, dmls2, dd2) && - //asym > CtsMCut (fSuper->GetAsymLo(), dmls, dmcza, dmls2, dd2) && + + asym < CtsMCut (fSuper->GetAsymUp(), dmls, dmcza, dmls2, dd2) && + asym > CtsMCut (fSuper->GetAsymLo(), dmls, dmcza, dmls2, dd2) && + dist < CtsMCut (fSuper->GetDistUp(), dmls, dmcza, dmls2, dd2) && - dist > CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) ) - fHadronness->SetHadronness(0.25); + dist > CtsMCut (fSuper->GetDistLo(), dmls, dmcza, dmls2, dd2) && + + conc < CtsMCut (fSuper->GetConcUp(), dmls, dmcza, dmls2, dd2) && + conc > CtsMCut (fSuper->GetConcLo(), dmls, dmcza, dmls2, dd2) && + + leakage < CtsMCut (fSuper->GetLeakage1Up(),dmls, dmcza, dmls2, dd2) && + leakage > CtsMCut (fSuper->GetLeakage1Lo(),dmls, dmcza, dmls2, dd2) ) + + fHadronness->SetHadronness(0.25); else - fHadronness->SetHadronness(0.75); + fHadronness->SetHadronness(0.75); + + //*fLog << "SChadroness = " << fHadronness->GetHadronness() << endl; fHadronness->SetReadyToSave(); @@ -257,2 +338,11 @@ return kTRUE; } + + + + + + + + + Index: trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.h (revision 2350) +++ trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.h (revision 2357) @@ -13,4 +13,6 @@ class MHillas; class MHillasSrc; +class MHillasExt; +class MNewImagePar; class MMcEvt; class MCerPhotEvt; @@ -25,4 +27,6 @@ MHillas *fHil; MHillasSrc *fHilSrc; + MHillasExt *fHilExt; + MNewImagePar *fNewPar; MMcEvt *fMcEvt; MHadronness *fHadronness; //! output container for hadronness @@ -32,9 +36,11 @@ TString fHilName; TString fHilSrcName; + TString fHilExtName; + TString fNewParName; TString fSuperName; // name of container for supercut parameters Double_t fMm2Deg; //! - Int_t fMap[8]; //! + Int_t fMap[11]; //! MHMatrix *fMatrix; //! @@ -63,2 +69,14 @@ #endif + + + + + + + + + + + + Index: trunk/MagicSoft/Mars/manalysis/MMinuitInterface.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MMinuitInterface.cc (revision 2350) +++ trunk/MagicSoft/Mars/manalysis/MMinuitInterface.cc (revision 2357) @@ -103,5 +103,5 @@ // 3 maximum output, showing progress of minimizations // - minuit.SetPrintLevel(1); + minuit.SetPrintLevel(-1); //.............................................. @@ -175,7 +175,8 @@ //.............................................. + // This doesn't seem to have any effect // Set maximum number of iterations (default = 500) //Int_t maxiter = 100000; - minuit.SetMaxIterations(10); + //minuit.SetMaxIterations(maxiter); //.............................................. @@ -211,9 +212,13 @@ if (nulloutput) fLog->SetNullOutput(kTRUE); - Double_t tmp = 0; - minuit.mnexcm("SIMPLEX", &tmp, 0, fErrMinimize); + Int_t maxcalls = 3000; + Double_t tolerance = 0.1; + Double_t tmp[2]; + tmp[0] = maxcalls; + tmp[1] = tolerance; + minuit.mnexcm("SIMPLEX", &tmp[0], 2, fErrMinimize); if (nulloutput) fLog->SetNullOutput(kFALSE); - //*fLog << "return from SIMPLEX" << endl; + *fLog << "return from SIMPLEX" << endl; } @@ -269,5 +274,5 @@ // 2 values, errors, step sizes, internal values // 3 values, errors, step sizes, 1st derivatives - // 4 values, paraboloc errors, MINOS errors + // 4 values, parabolic errors, MINOS errors //Int_t nkode = 4; @@ -292,2 +297,11 @@ return kTRUE; } + + + + + + + + +