Index: trunk/MagicSoft/Mars/manalysis/MMatrixLoop.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MMatrixLoop.cc (revision 6922) +++ trunk/MagicSoft/Mars/manalysis/MMatrixLoop.cc (revision 6924) @@ -46,5 +46,5 @@ // through. // -MMatrixLoop::MMatrixLoop(MHMatrix *mat, const char *name, const char *title) : fMatrix(mat) +MMatrixLoop::MMatrixLoop(MHMatrix *mat, const char *name, const char *title) : fMatrix(mat), fOperationMode(kDefault) { fName = name ? name : gsDefName.Data(); @@ -75,5 +75,5 @@ Int_t MMatrixLoop::PreProcess(MParList *plist) { - fNumRow = 0; + fNumRow = fOperationMode==kOdd ? 1 : 0; return fMatrix ? kTRUE : kFALSE; @@ -87,4 +87,6 @@ Int_t MMatrixLoop::Process() { - return fMatrix->SetNumRow(fNumRow++); + const Int_t rc = fMatrix->SetNumRow(fNumRow); + fNumRow += fOperationMode==kDefault ? 1 : 2; + return rc; } Index: trunk/MagicSoft/Mars/manalysis/MMatrixLoop.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MMatrixLoop.h (revision 6922) +++ trunk/MagicSoft/Mars/manalysis/MMatrixLoop.h (revision 6924) @@ -10,4 +10,10 @@ class MMatrixLoop : public MRead { +public: + enum OperationMode_t { + kDefault, + kEven, + kOdd + }; private: // MMatrixLoop @@ -17,4 +23,6 @@ MHMatrix *fMatrix; Int_t fNumRow; //! Number of dimensions of histogram + + Byte_t fOperationMode; // MRead @@ -33,4 +41,6 @@ MMatrixLoop(MHMatrix *mat, const char *name=NULL, const char *title=NULL); + void SetOperationMode(OperationMode_t mode) { fOperationMode = mode; } + ClassDef(MMatrixLoop, 0) // Task 'reading' events from a MHMatrix }; Index: trunk/MagicSoft/Mars/mhflux/MHEnergyEst.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHEnergyEst.cc (revision 6922) +++ trunk/MagicSoft/Mars/mhflux/MHEnergyEst.cc (revision 6924) @@ -78,14 +78,14 @@ fHResolution.SetDirectory(NULL); fHResolution.SetName("EnergyRes"); - fHResolution.SetTitle("Histogram in \\frac{\\Delta E}{E_{mc}} vs E_{est} and E_{mc}"); + fHResolution.SetTitle("Histogram in \\frac{\\Delta lg(E)}{lg(E_{mc})} vs E_{est} and E_{mc}"); fHResolution.SetXTitle("E_{est} [GeV]"); fHResolution.SetYTitle("E_{mc} [GeV]"); - fHResolution.SetZTitle("\\frac{E_{est} - E_{mc}}{E_{mc}}"); + fHResolution.SetZTitle("\\frac{lg(E_{est}) - lg(E_{mc})}{lg(E_{mc})}"); fHImpact.SetDirectory(NULL); fHImpact.SetName("Impact"); - fHImpact.SetTitle("\\frac{\\Delta E}{E} vs Impact parameter"); + fHImpact.SetTitle("\\frac{\\Delta lg(E)}{lg(E)} vs Impact parameter"); fHImpact.SetXTitle("Impact parameter [m]"); - fHImpact.SetYTitle("\\frac{E_{est} - E_{mc}}{E_{mc}}"); + fHImpact.SetYTitle("\\frac{lg(E_{est}) - lg(E_{mc})}{lg(E_{mc})}"); fHEnergy.Sumw2(); @@ -97,5 +97,5 @@ binst.SetEdgesCos(50, 0, 60); binsi.SetEdges(10, 0, 400); - binsr.SetEdges(10, 0, 1); + binsr.SetEdges(50, -1.3, 1.3); SetBinning(&fHEnergy, &binse, &binse, &binst); @@ -164,11 +164,12 @@ const Double_t imp = fMatrix ? GetVal(1) : fMcEvt->GetImpact()/100; const Double_t theta = fMatrix ? GetVal(2) : fMcEvt->GetTelescopeTheta()*TMath::RadToDeg(); - const Double_t res = (eest-etru)/etru; + const Double_t res = (log10(eest)-log10(etru));///log10(etru); + const Double_t resE = (eest-etru)/etru; fHEnergy.Fill(eest, etru, theta, w); - fHResolution.Fill(eest, etru, TMath::Abs(res), w); - fHImpact.Fill(imp, TMath::Abs(res), w); - - fChisq += TMath::Abs(res);//*res; + fHResolution.Fill(eest, etru, resE, w); + fHImpact.Fill(imp, resE, w); + + fChisq += res*res; fBias += res; @@ -181,10 +182,9 @@ fBias /= GetNumExecutions(); - Double_t res = fChisq; //TMath::Sqrt(fChisq - fBias*fBias); - - fResult->SetVal(TMath::IsNaN(res)?0:res);/// GetNumExecutions()); - - *fLog << all << "Mean Energy Resoltuion: " << Form("%.1f%%", fResult->GetVal()*100) << endl; - *fLog << all << "Energy Bias at: " << Form("%.1f%%", fBias*100) << endl; + Double_t sigma = TMath::Sqrt(fChisq);//+TMath::Abs(fBias); + fResult->SetVal(sigma); + + *fLog << all << "Mean log10(Energy) Resoltion: " << Form("%.1f%%", TMath::Sqrt(fChisq-fBias*fBias)*100) << endl; + *fLog << all << "Mean log10(Energy) Bias: " << Form("%.1f%%", fBias*100) << endl; return kTRUE; @@ -239,5 +239,6 @@ } - //pad->GetPad(1)->GetPad(2)->cd(2); + pad->GetPad(1)->GetPad(2)->cd(2); + /*=*/fHResolution.ProjectionZ("Resolution"); ///*h =*/ fHImpact.ProjectionX("Impact", -1, 9999, "e"); } @@ -294,4 +295,6 @@ gPad->SetBorderMode(0); gPad->SetLogx(); + gPad->SetGridx(); + gPad->SetGridy(); //gROOT->GetListOfCleanups()->Add(gPad); // WHY? @@ -301,4 +304,5 @@ h2->SetBit(kCanDelete); h2->SetFillColor(kBlue); + h2->SetLineColor(kRed); TH1D *h1 = h2->ProfileX(Form("Prof%s", h.GetName()), -1, 9999, "s"); @@ -310,7 +314,6 @@ h1->SetStats(kFALSE); - - //h1->Draw("E3"); - h2->Draw(); + h2->Draw(""); + h1->Draw("E3same"); h1->Draw("Chistsame"); @@ -346,4 +349,6 @@ gPad->SetBorderMode(0); + gPad->SetGridx(); + gPad->SetGridy(); gPad->SetLogx(); h = (TH1D*)fHEnergy.Project3D("ey"); @@ -376,5 +381,20 @@ pad3->Divide(2, 1, 1e-10, 1e-10); pad3->cd(1); - gPad->SetBorderMode(0);/* + gPad->SetBorderMode(0); + gPad->SetGridx(); + gPad->SetGridy(); + h = fHEnergy.Project3D("ez"); + h->SetTitle("Zenith Angle Distribution"); + h->SetBit(TH1::kNoStats); + h->SetDirectory(NULL); + h->SetXTitle("\\Theta [\\circ]"); + h->SetBit(kCanDelete); + h->Draw(); + + pad3->cd(2); + gPad->SetBorderMode(0); + gPad->SetGridx(); + gPad->SetGridy(); + /* h = fHImpact.ProjectionX("Impact", -1, 9999, "e"); h->SetBit(TH1::kNoStats); @@ -384,14 +404,14 @@ h->SetBit(kCanDelete); h->Draw();*/ - h = fHEnergy.Project3D("ez"); - h->SetTitle("Distribution of Theta"); - h->SetBit(TH1::kNoStats); + // ---------------------------------------- + h = fHResolution.ProjectionZ("Resolution"); + //h->SetXTitle("\\frac{lg(E_{est}) - lg(E_{mc})}{lg(E_{mc})}"); + h->SetYTitle("Counts"); + h->SetTitle("Distribution of \\Delta lg(E) / lg(E_{mc})"); h->SetDirectory(NULL); - h->SetXTitle("\\Theta [\\circ]"); h->SetBit(kCanDelete); - h->Draw(); - - pad3->cd(2); - gPad->SetBorderMode(0); + h->GetXaxis()->SetRangeUser(-1.3, 1.3); + h->Draw(""); + // ---------------------------------------- pad->cd(2); @@ -413,16 +433,16 @@ h = MakePlot(fHResolution, "zy"); h->SetXTitle("E_{mc} [GeV]"); - h->SetYTitle("\\frac{E_{est} - E_{mc}}{E_{mc}}"); - h->SetTitle("Energy resolution \\Delta E / E vs Monte Carlo energy E_{mc}"); - h->SetMinimum(0); - h->SetMaximum(1); + h->SetYTitle("\\frac{lg(E_{est}) - lg(E_{mc})}{lg(E_{mc})}"); + h->SetTitle("Energy resolution \\Delta lg(E) / lg(E) vs Monte Carlo energy E_{mc}"); + h->SetMinimum(-1.3); + h->SetMaximum(1.3); pad2->cd(3); h = MakePlot(fHResolution, "zx"); h->SetXTitle("E_{est} [GeV]"); - h->SetYTitle("\\frac{E_{est} - E_{mc}}{E_{mc}}"); - h->SetTitle("Energy Resolution \\Delta E / E vs estimated Energy E_{est}"); - h->SetMinimum(0); - h->SetMaximum(1); + h->SetYTitle("\\frac{lg(E_{est}) - lg(E_{mc})}{lg(E_{mc})}"); + h->SetTitle("Energy Resolution \\Delta lg(E) / lg(E) vs estimated Energy E_{est}"); + h->SetMinimum(-1.3); + h->SetMaximum(1.3); } Index: trunk/MagicSoft/Mars/mjobs/MJOptimize.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJOptimize.cc (revision 6922) +++ trunk/MagicSoft/Mars/mjobs/MJOptimize.cc (revision 6924) @@ -228,5 +228,5 @@ } -MJOptimize::MJOptimize() : fDebug(-1), fNumEvents(0), fType(kSimplex), fNumMaxCalls(0), fTolerance(0) +MJOptimize::MJOptimize() : fDebug(-1), fNumEvents(0), fType(kSimplex), fNumMaxCalls(0), fTolerance(0), fTestTrain(kFALSE) { fRules.SetOwner(); @@ -642,4 +642,6 @@ } + MMatrixLoop *loop = dynamic_cast(parlist.FindTask("MRead")); + TString txt("Starting "); switch (fType) @@ -666,4 +668,9 @@ *fLog << inf << "Number of Parameters: " << fParameters.GetSize() << endl; + // In case the reader is the matrix loop and testrain is enabled + // switch on even mode... + if (loop && fTestTrain) + loop->SetOperationMode(MMatrixLoop::kEven); + if (!Optimize(evtloop)) return kFALSE; @@ -672,7 +679,17 @@ fEvtLoop->SetDisplay(fDisplay); - return Fcn(fParameters); + if (!Fcn(fParameters)) + return kFALSE; + + // In case the reader is the matrix loop and testrain is enabled + // switch on odd mode... + if (!loop || !fTestTrain) + return kTRUE; + + loop->SetOperationMode(MMatrixLoop::kOdd); + // Done already in Fcn // list.SetVariables(fParameters); + return Fcn(fParameters); } Index: trunk/MagicSoft/Mars/mjobs/MJOptimize.h =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJOptimize.h (revision 6922) +++ trunk/MagicSoft/Mars/mjobs/MJOptimize.h (revision 6924) @@ -84,4 +84,5 @@ UInt_t fNumMaxCalls; Float_t fTolerance; + Bool_t fTestTrain; Bool_t Optimize(MEvtLoop &evtloop); @@ -108,4 +109,5 @@ void SetNumMaxCalls(UInt_t num=0) { fNumMaxCalls=num; } void SetTolerance(Float_t tol=0) { fTolerance=tol; } + void EnableTestTrain(Bool_t b=kTRUE) { fTestTrain=b; } // Parameter access