Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 2106) +++ trunk/MagicSoft/Mars/Changelog (revision 2107) @@ -1,3 +1,17 @@ -*-*- END OF LINE -*-*- + + 2003/05/09: Abelardo Moralejo + + * mhistmc/MHMcEnergyMigration.[h,cc] + - Added histograms, changed Draw() to display them. Still + provisional, many changes in the whole part of the energy + estimator are needed. + + * macros/CT1EEst.C, CT1EnergyEst.C + - Removed + + * macros/CT1EgyEst.C + - Added example on how to use the energy estimation for CT1. + Very provisional! 2003/05/09: Wolfgang Wittek Index: trunk/MagicSoft/Mars/macros/CT1EEst.C =================================================================== --- trunk/MagicSoft/Mars/macros/CT1EEst.C (revision 2106) +++ (revision ) @@ -1,618 +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): Thomas Bretz, 09/2002 -! Abelardo Moralejo, 03/2003 -! Wolfgang Wittek, 04/2003 -! -! Copyright: MAGIC Software Development, 2000-2003 -! -! -\* ======================================================================== */ -// -//------------------------------------------------------------------------ -// -// fcn calculates the function to be minimized (using TMinuit::Migrad) -// -void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) -{ - MEvtLoop *evtloop = (MEvtLoop*)gMinuit->GetObjectFit(); - - MTaskList *tlist = (MTaskList*)evtloop->GetParList()->FindObject("MTaskList"); // GetTaskList(); - - MChisqEval *eval = (MChisqEval*) tlist->FindObject("MChisqEval"); - MEnergyEstParam *eest = (MEnergyEstParam*)tlist->FindObject("MEnergyEstParam"); - - eest->SetCoeff(TArrayD(eest->GetNumCoeff(), par)); - - evtloop->Eventloop(); - - f = eval->GetChisq(); -} -//------------------------------------------------------------------------ - -void CT1EgyEst() -{ - TString inPath = "~magican/ct1test/wittek/marsoutput/optionA/"; - TString fileNameIn = "MC2.root"; - - TString outPath = "~magican/ct1test/wittek/marsoutput/optionA/"; - TString energyParName = "energyest.root"; - - TString hilName = "MHillas"; - TString hilSrcName = "MHillasSrc"; - - TString hadronnessName = "MHadronness"; - - Int_t howMany = 2200; - - Float_t maxhadronness = 0.4; - Float_t maxalpha = 20.0; - - CT1EEst(inPath, fileNameIn, - outPath, energyParName, - hilName, hilSrcName, hadronnessName, - howMany, maxhadronness, maxalpha); -} - - - -//------------------------------------------------------------------------ -// -// Arguments of CT1EEst : -// ------------------------ -// -// inPath, fileNameIn path and name of input file (MC gamma events) -// outPath, energyParName path and name of file containing the parameters -// of the energy estimator -// hilName, hilSrcName names of "MHillas" and "MHillasSrc" containers -// hadronnessName name of container holding the hadronness -// howMany no.of events to be read from input file -// maxhadronness maximum value of hadronness to be accepted -// maxalpha maximum value of |ALPHA| to be accepted -// -void CT1EEst(TString inPath, TString fileNameIn, - TString outPath, TString energyParName, - TString hilName, TString hilSrcName, TString hadronnessName, - Int_t howMany, - Float_t maxhadronness, Float_t maxalpha) -{ - cout << "================================================================" - << endl; - cout << "Start of energy estimation part" << endl; - cout << "" << endl; - cout << "CT1EEst input values : " << endl; - cout << "---------------------- " << endl; - cout << "inPath, fileNameIn = " - << inPath << ", " << fileNameIn << endl; - cout << "outPath, energyParName = " - << outPath << ", " << energyParName << endl; - cout << "hilName, hilSrcName, hadronnessName = " << hilName << ", " - << hilSrcName << ", " << hadronnessName << endl; - cout << "howMany, maxhadronness, maxalpha = " << howMany << ", " - << maxhadronness << ", " << maxalpha << endl; - cout << "" << endl; - - - // convert from [deg] to [mm] - // - MGeomCamCT1 gct1; - //maxdist /= gct1.GetConvMm2Deg(); - - - //========================================================================== - // Start of Part 1 (determination of the parameters of the energy estimator) - // - - //----------------------------------------------------------------------- - // Fill events into a MHMatrix, - // to be used later in the minimization by MINUIT - // - - MParList parlist; - - TString inputfile(outPath); - inputfile += fileNameIn; - - MFEventSelector selector; - selector.SetNumSelectEvts(howMany); - - cout << "Read events from file '" << inputfile << "'" << endl; - MReadTree read("Events"); - read.AddFile(inputfile); - read.DisableAutoScheme(); - read.SetSelector(&selector); - - MFCT1SelFinal filterhadrons; - filterhadrons.SetCuts(maxhadronness, maxalpha); - filterhadrons.SetHadronnessName(hadronnessName); - filterhadrons.SetInverted(); - - cout << "Define columns of matrix" << endl; - MHMatrix matrix; - Int_t colenergy = matrix.AddColumn("MMcEvt.fEnergy"); - Int_t colimpact = matrix.AddColumn("MMcEvt.fImpact"); - - if (colenergy < 0 || colimpact < 0) - { - cout << "colenergy, colimpact = " << colenergy << ", " - << colimpact << endl; - return; - } - - MEnergyEstParam eest(hilName); - eest.Add(hilSrcName); - eest.InitMapping(&matrix); - - cout << "--------------------------------------" << endl; - cout << "Fill events into the matrix" << endl; - if ( !matrix.Fill(&parlist, &read, &filterhadrons) ) - return; - cout << "Matrix was filled with " << matrix.GetNumRows() - << " events" << endl; - - //----------------------------------------------------------------------- - // - // Setup the eventloop which will be executed in the fcn of MINUIT - // - cout << "--------------------------------------" << endl; - cout << "Setup event loop to be used in 'fcn'" << endl; - - MParList parlist; - MTaskList tasklist; - - MMatrixLoop loop(&matrix); - - - MChisqEval eval; - eval.SetY1(new MDataElement(&matrix, colenergy)); - eval.SetY2(new MDataMember("MEnergyEst.fEnergy")); - eval.SetOwner(); - - //******************************** - // Entries in MParList - - parlist.AddToList(&tasklist); - - //******************************** - // Entries in MTaskList - - tasklist.AddToList(&loop); - tasklist.AddToList(&eest); - tasklist.AddToList(&eval); - - //******************************** - - cout << "Event loop was setup" << endl; - MEvtLoop evtloop; - evtloop.SetParList(&parlist); - - - - //.......................................................................... - - //---------- Start of minimization part -------------------- - // - // Do the minimization with MINUIT - // - // Be careful: This is not thread safe - // - cout << "--------------------------------------" << endl; - cout << "Start minimization" << endl; - - TMinuit minuit(12); - minuit.SetPrintLevel(-1); - minuit.SetFCN(fcn); - - // Ready for: minuit.mnexcm("SET ERR", arglist, 1, ierflg) - - if (minuit.SetErrorDef(1)) - { - cout << "SetErrorDef failed." << endl; - return; - } - - // - // Set initial values of the parameters (close enough to the final ones, taken - // from previous runs of the procedure). Parameter fA[4] is not used in the default - // energy estimation model (from D. Kranich). - // - TArrayD fA(5); - TArrayD fB(7); - - fA[0] = 21006.2; - fA[1] = -43.2648; - fA[2] = -690.403; - fA[3] = -0.0428544; - fA[4] = 1.; - fB[0] = -3391.05; - fB[1] = 136.58; - fB[2] = 0.253807; - fB[3] = 254.363; - fB[4] = 61.0386; - fB[5] = -0.0190984; - fB[6] = -421695; - - // - // Set starting values and step sizes for parameters - // - for (Int_t i=0; iWrite("EnergyParams"); - outfile.Close(); - - cout << "--------------------------------------" << endl; - cout << "Write parameters of energy estimator onto file '" - << paramout << endl; - cout << "--------------------------------------" << endl; - // - // End of Part 1 (determination of the parameters of the energy estimator) - //========================================================================= - - - //////////////////////////////////////////////////////////////////////////// - //////////////////////////////////////////////////////////////////////////// - //////////////////////////////////////////////////////////////////////////// - //////////////////////////////////////////////////////////////////////////// - - - //========================================================================== - // Start of Part 2 (test quality of energy estimation) - // - // - - //-------------------------------------------- - // Read the parameters of the energy estimator - // - - TString paramin(outPath); - paramin += energyParName; - - cout << "--------------------------------------" << endl; - cout << "Read parameters of energy estimator from file '" - << paramin << endl; - TFile enparam(paramin); - - //======================================================================= - // Setup the event loop - // - cout << "--------------------------------------" << endl; - cout << "Setup event loop for checking quality of energy estimation" << endl; - - - MTaskList tlist2; - MParList parlist2; - - //----------------------------------------------- - // Read events - // - - TString inputfile(inPath); - inputfile += fileNameIn; - - cout << "Read events from file '" << inputfile << "'" << endl; - MReadTree read2("Events"); - read2.AddFile(inputfile); - read2.DisableAutoScheme(); - - - //----------------------------------------------- - // Select events - // - cout << "Select events with hadronness < " << maxhadronness - << " and |alpha| < " << maxalpha << endl; - MFCT1SelFinal hcut2; - hcut2.SetCuts(maxhadronness, maxalpha); - hcut2.SetHadronnessName(hadronnessName); - - MContinue cont(&hcut2); - - - //----------------------------------------------- - // Create some histograms to display the results - // - - MH3 mh3ed ("log10(MMcEvt.fEnergy)", "MEnergyEst.fEnergy/MMcEvt.fEnergy-1.0"); - MH3 mh3ed2("log10(MEnergyEst.fEnergy)", "MEnergyEst.fEnergy/MMcEvt.fEnergy-1.0"); - MH3 mhe ("log10(MMcEvt.fEnergy)", "log10(MEnergyEst.fEnergy)"); - - mhe.SetName("HistEE"); - mh3ed.SetName("HistEnergyDelta"); - mh3ed2.SetName("HistEnergyDelta"); - - MBinning binsedx("BinningHistEnergyDeltaX"); - MBinning binsedy("BinningHistEnergyDeltaY"); - MBinning binseex("BinningHistEEX"); - MBinning binseey("BinningHistEEY"); - - binsedx.SetEdges(35, 2.0, 5.5); - binsedy.SetEdges(40,-1.5, 2.5); - - binseex.SetEdges(35, 2.0, 5.5); - binseey.SetEdges(35, 2.0, 5.5); - - MFillH hfilled(&mh3ed); - MFillH hfilled2(&mh3ed2); - MFillH hfillee(&mhe); - - - //----------------------------------------------- - // Create energy estimator task, add necessary containers and - // initialize parameters read from file: - // - - MEnergyEstParam eest2(hilName); - eest2.Add(hilSrcName); - - TArrayD parA(5); - TArrayD parB(7); - - for (Int_t i=0; iGetXaxis()->SetTitle("log_{10} E_{MC} (GeV)"); - mhe.GetHist()->GetYaxis()->SetTitle("log_{10} E_{EST} (GeV)"); - mhe.GetHist()->GetXaxis()->SetTitleOffset(1.2); - - TCanvas *c = new TCanvas("energy","CT1 Energy estimation"); - c->Divide(2,2); - c->cd(1); - energy_1->SetBottomMargin(0.12); - mhe.DrawClone("PROFXnonew"); - - mh3ed.GetHist()->GetXaxis()->SetTitle("log_{10} E_{MC} (GeV)"); - mh3ed.GetHist()->GetYaxis()->SetTitle("\\frac{E_{EST} - E_{MC}}{E_{MC}}"); - mh3ed.GetHist()->GetXaxis()->SetTitleOffset(1.2); - mh3ed.GetHist()->GetYaxis()->SetTitleOffset(1.5); - c->cd(2); - energy_2->SetLeftMargin(0.15); - energy_2->SetBottomMargin(0.12); - mh3ed.DrawClone("PROFXnonew"); - - mh3ed2.GetHist()->GetXaxis()->SetTitle("log_{10} E_{EST} (GeV)"); - mh3ed2.GetHist()->GetYaxis()->SetTitle("\\frac{E_{EST} - E_{MC}}{E_{MC}}"); - mh3ed2.GetHist()->GetXaxis()->SetTitleOffset(1.2); - mh3ed2.GetHist()->GetYaxis()->SetTitleOffset(1.5); - c->cd(3); - energy_3->SetLeftMargin(0.15); - energy_3->SetBottomMargin(0.12); - mh3ed2.DrawClone("PROFXnonew"); - - ((TH2*)mh3ed2.GetHist())->ProjectionY("deltae"); - - c->cd(4); - energy_4->SetLeftMargin(0.1); - energy_4->SetRightMargin(0.05); - energy_4->SetBottomMargin(0.12); - deltae.GetXaxis()->SetTitleOffset(1.2); - deltae.GetXaxis()->SetTitle("(E_{EST} - E_{MC}) / E_{MC}"); - deltae.Draw("e"); - - //------------------------------------------- - // - // include the histograms - // into the root file containing the parameters of the energy estimator - // - - - TString paramout(outPath); - paramout += energyParName; - - TFile out2(paramout, "UPDATE"); - ((TH2*)mh3ed.GetHist())->Write("mh3ed"); - ((TH2*)mh3ed2.GetHist())->Write("mh3ed2"); - ((TH2*)mhe.GetHist())->Write("mhe"); - deltae.Write(); - out2.Close(); - - cout << "Quality histograms were added onto the file '" << paramout << endl; - - cout << "" << endl; - cout << "End of energy estimation part" << endl; - cout << "================================================================" - << endl; - // - // End of Part 2 (test quality of energy estimation) - //========================================================================== - - - return; - -} -//============================================================================ - - - - - - - - - Index: trunk/MagicSoft/Mars/macros/CT1EgyEst.C =================================================================== --- trunk/MagicSoft/Mars/macros/CT1EgyEst.C (revision 2107) +++ trunk/MagicSoft/Mars/macros/CT1EgyEst.C (revision 2107) @@ -0,0 +1,301 @@ +/* ======================================================================== *\ +! +! * +! * 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): Thomas Bretz, 09/2002 +! Abelardo Moralejo, 03/2003 +! Wolfgang Wittek, 04/2003 +! +! Copyright: MAGIC Software Development, 2000-2003 +! +! +\* ======================================================================== */ + +#include "MParList.h" +#include "MTaskList.h" +#include "MGeomCamCT1.h" +#include "MFEventSelector.h" +#include "MReadTree.h" +#include "MFCT1SelFinal.h" +#include "MHMatrix.h" +#include "MEnergyEstParam.h" +#include "MMatrixLoop.h" +#include "MChisqEval.h" +#include "MLog.h" +#include "MLogManip.h" + +void CT1EgyEst() +{ + // TString inPath = "~magican/ct1test/wittek/marsoutput/optionA/"; + TString inPath = "./"; + TString fileNameIn = "MC1.root"; + + // TString outPath = "~magican/ct1test/wittek/marsoutput/optionA/"; + TString outPath = "./"; + TString energyParName = "energyest.root"; + + TString hilName = "MHillas"; + TString hilSrcName = "MHillasSrc"; + + // TString hadronnessName = "MHadronness"; + TString hadronnessName = "HadNN"; + + // Int_t howMany = 10000; + Int_t howMany = 2000; + + Float_t maxhadronness = 0.3; + Float_t maxalpha = 20.0; + Float_t maxdist = 1.; + + CT1EEst(inPath, fileNameIn, + outPath, energyParName, + hilName, hilSrcName, hadronnessName, + howMany, maxhadronness, maxalpha, maxdist); +} + +//------------------------------------------------------------------------ +// +// Arguments of CT1EEst : +// ------------------------ +// +// inPath, fileNameIn path and name of input file (MC gamma events) +// outPath, energyParName path and name of file containing the parameters +// of the energy estimator +// hilName, hilSrcName names of "MHillas" and "MHillasSrc" containers +// hadronnessName name of container holding the hadronness +// howMany no.of events to be read from input file +// maxhadronness maximum value of hadronness to be accepted +// maxalpha maximum value of |ALPHA| to be accepted +// maxdist maximum value of DIST (degrees) to be accepted +// + +void CT1EEst(TString inPath, TString fileNameIn, + TString outPath, TString energyParName, + TString hilName, TString hilSrcName, TString hadronnessName, + Int_t howMany, + Float_t maxhadronness, Float_t maxalpha, Float_t maxdist) +{ + cout << "================================================================" + << endl; + cout << "Start of energy estimation part" << endl; + cout << "" << endl; + cout << "CT1EEst input values : " << endl; + cout << "---------------------- " << endl; + cout << "inPath, fileNameIn = " + << inPath << ", " << fileNameIn << endl; + cout << "outPath, energyParName = " + << outPath << ", " << energyParName << endl; + cout << "hilName, hilSrcName, hadronnessName = " << hilName << ", " + << hilSrcName << ", " << hadronnessName << endl; + cout << "howMany, maxhadronness, maxalpha, maxdist = " << howMany << ", " + << maxhadronness << ", " << maxalpha << ", " << maxdist << endl; + cout << "" << endl; + + + //========================================================================== + // Start of Part 1 (determination of the parameters of the energy estimator) + // + + //----------------------------------------------------------------------- + // Fill events into a MHMatrix, + // to be used later in the minimization by MINUIT + // + + MMcEnergyEst Optimize; + + TString inputfile(outPath); + inputfile += fileNameIn; + Optimize.SetInFile(inputfile); + + TString paramout(outPath); + paramout += energyParName; + Optimize.SetOutFile(paramout); + + MFCT1SelFinal filterhadrons; + filterhadrons.SetHadronnessName(hadronnessName); + filterhadrons.SetCuts(maxhadronness, maxalpha, maxdist); + filterhadrons.SetInverted(); + Optimize.SetEventFilter(&filterhadrons); + + Optimize.SetNevents(howMany); + + // + // Find the optimal parameters for the energy estimator: + // + Optimize.FindParams(); + + cout << "--------------------------------------" << endl; + cout << "Write parameters of energy estimator onto file '" + << paramout << endl; + cout << "--------------------------------------" << endl; + + Optimize.Print(); + + TFile out("energyest.root","recreate"); + Optimize.SetReadyToSave(); + Optimize.Write(); + out.Close(); + + // + // END of Part 1 + //========================================================================== + + + + //========================================================================== + // Start of Part 2 (test quality of energy estimation) + // + // + + //-------------------------------------------- + // Read the parameters of the energy estimator + // + + TString paramin(outPath); + paramin += energyParName; + + TFile enparam(paramin); + MMcEnergyEst mcest; + mcest.Read("MMcEnergyEst"); + enparam.Close(); + + cout << "--------------------------------------" << endl; + cout << "Read parameters of energy estimator from file '" + << paramin << endl; + + TArrayD parA(5); + TArrayD parB(7); + + for (Int_t i=0; i -! Abelardo Moralejo, 03/2003 -! -! Copyright: MAGIC Software Development, 2000-2002 -! -! -\* ======================================================================== */ - -// -// fcn is the function to be minimized using TMinuit::Migrad -// -void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) -{ - MEvtLoop *evtloop = (MEvtLoop*)gMinuit->GetObjectFit(); - - MTaskList *tlist = (MTaskList*)evtloop->GetParList()->FindObject("MTaskList"); // GetTaskList(); - - MChisqEval *eval = (MChisqEval*) tlist->FindObject("MChisqEval"); - MEnergyEstParam *eest = (MEnergyEstParam*)tlist->FindObject("MEnergyEstParam"); - - eest->SetCoeff(TArrayD(eest->GetNumCoeff(), par)); - - evtloop->Eventloop(); - - f = eval->GetChisq(); -} - -//------------------------------------------------------------------------ - -void CT1EnergyEst(Float_t maxhadronness=0.22, Float_t maxdist = 1.) -{ - // - // Upper cut in Maxdist taken from D. Kranich's Ph D. - // We have to convert maxdist from degrees to mm: - // - - MGeomCamCT1 gct1; - maxdist /= gct1.GetConvMm2Deg(); - - // - // Fill events into a MHMatrix - // - MParList parlist; - MHMatrix matrix; - - // - // colenergy and colimpact are the indexes of the matrix columns containing - // the MC energy and impact parameter. - // - Int_t colenergy = matrix.AddColumn("MMcEvt.fEnergy"); - Int_t colimpact = matrix.AddColumn("MMcEvt.fImpact"); - - MEnergyEstParam eest("Hillas"); - eest.Add("HillasSrc"); - eest.InitMapping(&matrix); - - // - // Here we read in the events for the training. Second argument of AddFile is - // the number of events. With 2000 it is rather fast, but only the lowest zenith - // angles will enter in the minimization. - // - MReadTree read("Events"); - read.AddFile("MC_ON2.root",2000); - read.DisableAutoScheme(); - - // - // Filter to keep the most gamma-like events: - // - TString hcut("MHadronness.fHadronness<"); - hcut += maxhadronness; - hcut += " && HillasSrc.fDist < "; - hcut += maxdist; - - MF filterhadrons(hcut); - - // - // Fill the matrix used later in the minimization loop by Minuit: - // - if (!matrix.Fill(&parlist, &read, &filterhadrons)) - return; - - // - // Setup the tasklist used to evaluate the needed chisq - // - MTaskList tasklist; - parlist.AddToList(&tasklist); - - MMatrixLoop loop(&matrix); - - MChisqEval eval; - - eval.SetY1(new MDataElement(&matrix, colenergy)); - eval.SetY2(new MDataMember("MEnergyEst.fEnergy")); - - eval.SetOwner(); - - tasklist.AddToList(&loop); - tasklist.AddToList(&eest); - tasklist.AddToList(&eval); - - MEvtLoop evtloop; - evtloop.SetParList(&parlist); - - // - // Be careful: This is not thread safe - // - TMinuit minuit(12); - minuit.SetPrintLevel(-1); - minuit.SetFCN(fcn); - - // Ready for: minuit.mnexcm("SET ERR", arglist, 1, ierflg) - - if (minuit.SetErrorDef(1)) - { - cout << "SetErrorDef failed." << endl; - return; - } - - // - // Set initial values of the parameters (close enough to the final ones, taken - // from previous runs of the procedure). Parameter fA[4] is not used in the default - // energy estimation model (from D. Kranich). - // - TArrayD fA(5); - TArrayD fB(7); - - fA[0] = 21006.2; - fA[1] = -43.2648; - fA[2] = -690.403; - fA[3] = -0.0428544; - fA[4] = 1.; - fB[0] = -3391.05; - fB[1] = 136.58; - fB[2] = 0.253807; - fB[3] = 254.363; - fB[4] = 61.0386; - fB[5] = -0.0190984; - fB[6] = -421695; - - // - // Set starting values and step sizes for parameters - // - for (Int_t i=0; iWrite("EnergyParams"); - out.Close(); - - // - // End of Part 1 (estimation of the parameters) - // - - - //////////////////////////////////////////////////////////////////////////// - //////////////////////////////////////////////////////////////////////////// - //////////////////////////////////////////////////////////////////////////// - //////////////////////////////////////////////////////////////////////////// - - - // - // Part 2: Now test how the energy estimation method works. - // - // - // Create a empty Parameter List and an empty Task List - // The tasklist is identified in the eventloop by its name - // - - MTaskList tlist2; - MParList parlist2; - - parlist2.AddToList(&tlist2); - - // - // Now setup the tasks and tasklist: - // --------------------------------- - // - - MReadTree read2("Events"); - read2.AddFile("MC_ON2.root"); - read2.DisableAutoScheme(); - - // - // Create some histograms to display the results: - // - - MH3 mh3ed("log10(MMcEvt.fEnergy)", "MEnergyEst.fEnergy/MMcEvt.fEnergy-1"); - MH3 mh3ed2("log10(MEnergyEst.fEnergy)", "MEnergyEst.fEnergy/MMcEvt.fEnergy-1"); - MH3 mhe("log10(MMcEvt.fEnergy)", "log10(MEnergyEst.fEnergy)"); - - MFillH hfilled(&mh3ed); - MFillH hfilled2(&mh3ed2); - MFillH hfillee(&mhe); - - mhe.SetName("HistEE"); - mh3ed.SetName("HistEnergyDelta"); - mh3ed2.SetName("HistEnergyDelta"); - - MBinning binsedx("BinningHistEnergyDeltaX"); - MBinning binsedy("BinningHistEnergyDeltaY"); - MBinning binseex("BinningHistEEX"); - MBinning binseey("BinningHistEEY"); - - binsedx.SetEdges(25, 2.5, 5.); - binsedy.SetEdges(30, -1., 1.); - binseex.SetEdges(25, 2.5, 5.); - binseey.SetEdges(30, 2., 5.); - - parlist2.AddToList(&binsedx); - parlist2.AddToList(&binsedy); - parlist2.AddToList(&binseex); - parlist2.AddToList(&binseey); - - // - // Create energy estimator task, add necessary containers and - // initialize parameters read from file: - // - - MEnergyEstParam eest2("Hillas"); - eest2.Add("HillasSrc"); - - TFile enparam("energyest.root"); - TArrayD parA(5); - TArrayD parB(7); - - for (Int_t i=0; i"); - hcut2 += maxhadronness; - hcut2 += "|| HillasSrc.fDist > "; - hcut2 += maxdist; - - MContinue cont(hcut2); - - tlist2.AddToList(&cont); - tlist2.AddToList(&eest2); - - // tlist2.AddToList(new MPrint("MMcEvt")); - // tlist2.AddToList(new MPrint("MEnergyEst")); - - tlist2.AddToList(&hfilled); - tlist2.AddToList(&hfilled2); - tlist2.AddToList(&hfillee); - - // - // Create and setup the eventloop - // - MProgressBar bar; - - MEvtLoop evtloop2; - evtloop2.SetProgressBar(&bar); - evtloop2.SetParList(&parlist2); - - // - // Execute your analysis - // - if (!evtloop2.Eventloop()) - return; - - // - // Plot results: - // - mhe.GetHist()->GetXaxis()->SetTitle("log_{10} E_{MC} (GeV)"); - mhe.GetHist()->GetYaxis()->SetTitle("log_{10} E_{EST} (GeV)"); - mhe.GetHist()->GetXaxis()->SetTitleOffset(1.2); - - TCanvas *c = new TCanvas("energy","CT1 Energy estimation"); - c->Divide(2,2); - c->cd(1); - energy_1->SetBottomMargin(0.12); - mhe.DrawClone("PROFXnonew"); - - mh3ed.GetHist()->GetXaxis()->SetTitle("log_{10} E_{MC} (GeV)"); - mh3ed.GetHist()->GetYaxis()->SetTitle("\\frac{E_{EST} - E_{MC}}{E_{MC}}"); - mh3ed.GetHist()->GetXaxis()->SetTitleOffset(1.2); - mh3ed.GetHist()->GetYaxis()->SetTitleOffset(1.5); - c->cd(2); - energy_2->SetLeftMargin(0.15); - energy_2->SetBottomMargin(0.12); - mh3ed.DrawClone("PROFXnonew"); - - mh3ed2.GetHist()->GetXaxis()->SetTitle("log_{10} E_{EST} (GeV)"); - mh3ed2.GetHist()->GetYaxis()->SetTitle("\\frac{E_{EST} - E_{MC}}{E_{MC}}"); - mh3ed2.GetHist()->GetXaxis()->SetTitleOffset(1.2); - mh3ed2.GetHist()->GetYaxis()->SetTitleOffset(1.5); - c->cd(3); - energy_3->SetLeftMargin(0.15); - energy_3->SetBottomMargin(0.12); - mh3ed2.DrawClone("PROFXnonew"); - - ((TH2*)mh3ed2.GetHist())->ProjectionY("deltae"); - - c->cd(4); - energy_4->SetLeftMargin(0.1); - energy_4->SetRightMargin(0.05); - energy_4->SetBottomMargin(0.12); - deltae.GetXaxis()->SetTitleOffset(1.2); - deltae.GetXaxis()->SetTitle("(E_{EST} - E_{MC}) / E_{MC}"); - deltae.Draw("e"); - - TFile out2("energyest.root", "update"); - ((TH2*)mh3ed.GetHist())->Write("mh3ed"); - ((TH2*)mh3ed2.GetHist())->Write("mh3ed2"); - ((TH2*)mhe.GetHist())->Write("mhe"); - deltae.Write(); - out2.Close(); - - return; - -} Index: trunk/MagicSoft/Mars/mhistmc/MHMcEnergyMigration.cc =================================================================== --- trunk/MagicSoft/Mars/mhistmc/MHMcEnergyMigration.cc (revision 2106) +++ trunk/MagicSoft/Mars/mhistmc/MHMcEnergyMigration.cc (revision 2107) @@ -17,4 +17,5 @@ ! ! Author(s): Wolfgang Wittek 4/2002 +! Abelardo Moralejo 5/2003 ! ! Copyright: MAGIC Software Development, 2000-2002 @@ -45,4 +46,5 @@ #include "MLog.h" #include "MLogManip.h" +#include ClassImp(MHMcEnergyMigration); @@ -54,5 +56,5 @@ // MHMcEnergyMigration::MHMcEnergyMigration(const char *name, const char *title) - : fHist() + : fHist(), fHist2(), fHistImp() { // @@ -63,9 +65,20 @@ fHist.SetDirectory(NULL); - - fHist.SetTitle("3D-plot E-true E-est Theta"); - fHist.SetXTitle("E_{true} [GeV]"); - fHist.SetYTitle("E_{est} [GeV]"); + fHist.SetTitle("3D-plot E_{TRUE} E_{EST} \\Theta"); + fHist.SetXTitle("E_{TRUE} [GeV]"); + fHist.SetYTitle("E_{EST} [GeV]"); fHist.SetZTitle("\\Theta [\\circ]"); + + fHist2.SetDirectory(NULL); + fHist2.SetTitle("3D-plot \\Delta E / E vs E_{TRUE} E_{EST}"); + fHist2.SetXTitle("E_{TRUE} [GeV]"); + fHist2.SetYTitle("E_{EST} [GeV]"); + fHist2.SetZTitle("\\frac{E_{EST} - E_{TRUE}}{E_{TRUE}}"); + + fHistImp.SetDirectory(NULL); + fHistImp.SetTitle("\\Delta E / E vs Impact parameter"); + fHistImp.SetXTitle("Impact parameter (m)"); + fHistImp.SetYTitle("\\frac{E_{EST} - E_{TRUE}}{E_{TRUE}}"); + } @@ -99,6 +112,14 @@ SetBinning(&fHist, binsenergy, binsenergy, binstheta); - fHist.Sumw2(); + + const MBinning* binsde = (MBinning*)plist->FindObject("BinningDE"); + const MBinning* binsimpact = (MBinning*)plist->FindObject("BinningImpact"); + SetBinning(&fHistImp, binsimpact, binsde); + + fHistImp.Sumw2(); + + SetBinning(&fHist2, binsenergy, binsenergy, binsde); + fHist2.Sumw2(); return kTRUE; @@ -124,8 +145,9 @@ gPad->SetLogx(); h = fHist.Project3D("ex_pro"); - h->SetTitle("Distribution of E_{true}"); - h->SetXTitle("E_{true} [GeV]"); + h->SetTitle("Distribution of E_{TRUE}"); + h->SetXTitle("E_{TRUE} [GeV]"); h->SetBit(kCanDelete); h->Draw(opt); + h->SetDirectory(NULL); pad->cd(2); @@ -133,23 +155,133 @@ gPad->SetLogx(); h = fHist.Project3D("ey_pro"); - h->SetTitle("Distribution of E_{est}"); + h->SetTitle("Distribution of E_{EST}"); h->SetXTitle("E_{est} [GeV]"); h->SetBit(kCanDelete); h->Draw(opt); + Double_t minEest = h->GetXaxis()->GetXmin(); + h->SetDirectory(NULL); pad->cd(3); gPad->SetBorderMode(0); - h = fHist.Project3D("ez_pro"); + h = fHist.Project3D("z_pro"); h->SetTitle("Distribution of \\Theta"); h->SetXTitle("\\Theta [\\circ]"); h->SetBit(kCanDelete); + h->SetLineWidth(2); h->Draw(opt); + h->SetDirectory(NULL); pad->cd(4); gPad->SetBorderMode(0); - fHist.Draw(opt); + TH1D* hpx; + hpx = fHistImp.ProjectionX("_px", 1, fHistImp.GetNbinsY(),"e"); + hpx->SetTitle("Distribution of Impact parameter"); + hpx->SetXTitle("Impact parameter (m)"); + hpx->SetBit(kCanDelete); + hpx->Draw(opt); + hpx->SetDirectory(NULL); + fHistImp.SetDirectory(NULL); pad->Modified(); pad->Update(); + + TVirtualPad *pad2 = MakeDefCanvas((TString)GetName()+"2"); + pad2->SetBorderMode(0); + + AppendPad(""); + + pad2->Divide(2,2); + + TH2D *h2; + + pad2->cd(1); + gPad->SetBorderMode(0); + gPad->SetLogx(); + gPad->SetLogy(); + h2 = (TH2D*) fHist.Project3D("xy"); + + TProfile* h2pfx; + h2pfx = h2->ProfileX("_pfx", 1, h2->GetNbinsY(),"S"); + h2pfx->SetXTitle("E_{TRUE} (GeV)"); + h2pfx->SetYTitle("E_{EST} (GeV)"); + h2pfx->SetTitle("E_{EST} vs E_{TRUE}"); + h2pfx->SetBit(kCanDelete); + h2pfx->SetFillColor(41); + h2pfx->SetMinimum(minEest); + h2pfx->SetStats(kFALSE); + h2pfx->Draw("E3"); + + h2->SetBit(kCanDelete); + h2->SetFillColor(1); + h2->Draw("box,same"); + h2->SetDirectory(NULL); + + h2pfx->SetLineColor(2); + h2pfx->SetLineWidth(2); + h2pfx->Draw("L,histo,same"); + h2pfx->SetDirectory(NULL); + + pad2->cd(2); + gPad->SetBorderMode(0); + gPad->SetLogx(); + gPad->SetLeftMargin(0.15); + h2 = (TH2D*) fHist2.Project3D("zx"); + h2->SetBit(kCanDelete); + h2pfx = h2->ProfileX("_pfx", 1, h2->GetNbinsY(),"S"); + h2pfx->SetTitle("\\Delta E / E vs E_{TRUE}"); + h2pfx->SetXTitle("E_{TRUE} (GeV)"); + h2pfx->SetYTitle("\\frac{E_{EST} - E_{TRUE}}{E_{TRUE}}"); + h2pfx->SetBit(kCanDelete); + h2pfx->SetFillColor(41); + h2pfx->GetYaxis()->SetTitleOffset(1.4); + h2pfx->SetStats(kFALSE); + h2pfx->Draw("E3"); + h2->SetFillColor(1); + h2->Draw("same,box"); + h2->SetDirectory(NULL); + h2pfx->SetLineColor(2); + h2pfx->SetLineWidth(2); + h2pfx->Draw("L,histo,same"); + h2pfx->SetDirectory(NULL); + + pad2->cd(3); + gPad->SetBorderMode(0); + gPad->SetLogx(); + gPad->SetLeftMargin(0.15); + h2 = (TH2D*) fHist2.Project3D("zy"); + h2->SetBit(kCanDelete); + h2pfx = h2->ProfileX("_pfx", 1, h2->GetNbinsY(),"S"); + h2pfx->SetTitle("\\Delta E / E vs E_{EST}"); + h2pfx->SetXTitle("E_{TRUE} (GeV)"); + h2pfx->SetYTitle("\\frac{E_{EST} - E_{TRUE}}{E_{TRUE}}"); + h2pfx->SetBit(kCanDelete); + h2pfx->SetFillColor(41); + h2pfx->GetYaxis()->SetTitleOffset(1.4); + h2pfx->SetStats(kFALSE); + h2pfx->SetMinimum(-1.); + h2pfx->SetMaximum(1.); + h2pfx->Draw("E3"); + + h2->SetFillColor(1); + h2->Draw("same,box"); + h2->SetDirectory(NULL); + h2pfx->SetLineColor(2); + h2pfx->SetLineWidth(2); + h2pfx->Draw("L,histo,same"); + h2pfx->SetDirectory(NULL); + + pad2->cd(4); + gPad->SetBorderMode(0); + h = fHist2.ProjectionZ("_pz",1,fHist2.GetNbinsX(),1,fHist2.GetNbinsY(),"e"); + h->SetBit(kCanDelete); + h->Draw(); + h->SetDirectory(NULL); + + pad2->Modified(); + pad2->Update(); + + fHist.SetDirectory(NULL); + fHist2.SetDirectory(NULL); + } @@ -160,9 +292,11 @@ Bool_t MHMcEnergyMigration::Fill(const MParContainer *par, const Stat_t w) { - // MHillasSrc &hil = *(MHillasSrc*)par; - // fHist.Fill(fMcEvt->GetEnergy(), hil.GetSize()); - // get E-true from fMcEvt and E-est from fEnergy - fHist.Fill(fMcEvt->GetEnergy(), fEnergy->GetEnergy(), fMcEvt->GetTheta()*kRad2Deg); + + fHist.Fill(fMcEvt->GetEnergy(), fEnergy->GetEnergy(), fMcEvt->GetTelescopeTheta()*kRad2Deg, w); + + fHist2.Fill(fMcEvt->GetEnergy(), fEnergy->GetEnergy(), (fEnergy->GetEnergy()-fMcEvt->GetEnergy())/fMcEvt->GetEnergy(), w); + + fHistImp.Fill(fMcEvt->GetImpact()/100., (fEnergy->GetEnergy()-fMcEvt->GetEnergy())/fMcEvt->GetEnergy(), w); return kTRUE; Index: trunk/MagicSoft/Mars/mhistmc/MHMcEnergyMigration.h =================================================================== --- trunk/MagicSoft/Mars/mhistmc/MHMcEnergyMigration.h (revision 2106) +++ trunk/MagicSoft/Mars/mhistmc/MHMcEnergyMigration.h (revision 2107) @@ -8,4 +8,8 @@ #ifndef ROOT_TH3 #include +#endif + +#ifndef ROOT_TH2 +#include #endif @@ -21,4 +25,6 @@ TH3D fHist; + TH3D fHist2; + TH2D fHistImp; public: @@ -30,6 +36,8 @@ const TH3D *GetHist() { return &fHist; } const TH3D *GetHist() const { return &fHist; } + TH1 *GetHistByName(const TString name) { return &fHist; } - TH1 *GetHistByName(const TString name) { return &fHist; } + const TH3D *GetHist2() { return &fHist2; } + const TH2D *GetHistImp() { return &fHistImp; } void Draw(Option_t *option=""); @@ -40,5 +48,2 @@ #endif - - -