Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 1885) +++ /trunk/MagicSoft/Mars/Changelog (revision 1886) @@ -1,3 +1,14 @@ -*-*- END OF LINE -*-*- + + 2003/04/01: Abelardo Moralejo + + * macros/CT1EnergyEst.C + - added argument (maximum dist parameter), changed (reduced) output + histograms, added writing to (and reading from) a file the energy + estimation parameters and the histograms. Added comments. + + * manalysis/MEnergyEstParam.[h,cc] + - added member function GetCoeff. Changed comment. + 2003/03/31: Thomas Bretz Index: /trunk/MagicSoft/Mars/macros/CT1EnergyEst.C =================================================================== --- /trunk/MagicSoft/Mars/macros/CT1EnergyEst.C (revision 1885) +++ /trunk/MagicSoft/Mars/macros/CT1EnergyEst.C (revision 1886) @@ -16,5 +16,6 @@ ! ! -! Author(s): Thomas Bretz et al, 09/2002 +! Author(s): Thomas Bretz, 09/2002 +! Abelardo Moralejo, 03/2003 ! ! Copyright: MAGIC Software Development, 2000-2002 @@ -23,9 +24,12 @@ \* ======================================================================== */ +// +// 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(); + MTaskList *tlist = (MTaskList*)evtloop->GetParList()->FindObject("MTaskList"); // GetTaskList(); MChisqEval *eval = (MChisqEval*) tlist->FindObject("MChisqEval"); @@ -41,7 +45,14 @@ //------------------------------------------------------------------------ -void CT1EnergyEst(Float_t maxhadronness=1.) +void CT1EnergyEst(Float_t maxhadronness=0.22, Float_t maxdist = 1.) { - // Bool_t evalenergy=kFALSE; + // + // 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 @@ -50,5 +61,10 @@ MHMatrix matrix; - Int_t col = matrix.AddColumn("MMcEvt.fEnergy"); + // + // 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"); @@ -56,14 +72,27 @@ 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", 200); + 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); - if (!matrix.Fill(&parlist, &read)) + // + // Fill the matrix used later in the minimization loop by Minuit: + // + if (!matrix.Fill(&parlist, &read, &filterhadrons)) return; @@ -76,8 +105,9 @@ MMatrixLoop loop(&matrix); - MChisqEval eval; - eval.SetY1(new MDataElement(&matrix, col)); + + eval.SetY1(new MDataElement(&matrix, colenergy)); eval.SetY2(new MDataMember("MEnergyEst.fEnergy")); + eval.SetOwner(); @@ -97,4 +127,5 @@ // Ready for: minuit.mnexcm("SET ERR", arglist, 1, ierflg) + if (minuit.SetErrorDef(1)) { @@ -104,23 +135,27 @@ // - // Set initial values - // - TArrayD fA(5); - fA[0] =1;//4916.4; */-2414.75; - fA[1] =1;//149.549; */ 1134.28; - fA[2] =1;//-558.209; */ 132.932; - fA[3] =1;//0.270725; */ 0.292845; - fA[4] =1;//107.001; */ 107.001; - + // 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); - fB[0] = 1;//-8234.12; */ -4282.25; - fB[1] = 1;// 23.2153; */ 18.892; - fB[2] = 1;// 0.416372;*/ 0.193373; - fB[3] = 1;// 332.42; */ 203.803; - fB[4] = 1;// -0.701764;*/ -0.534876; - fB[5] = 1;//-0.0131774;*/ -0.00789539; - fB[6] = 1;//-0.162687;*/ 0.111913; - + + 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. // @@ -222,7 +279,8 @@ // - MTaskList tlist2; - parlist.Replace(&tlist2); + MParList parlist2; + + parlist2.AddToList(&tlist2); // @@ -231,122 +289,86 @@ // - - read->SetEventNum(0); - - // - // Use this to change the binnign of the histograms to CT1-style - // - Bool_t usect1 = kTRUE; - - MH3 mh3e("MMcEvt.fEnergy", "(MEnergyEst.fEnergy/MMcEvt.fEnergy-1)*(MEnergyEst.fEnergy/MMcEvt.fEnergy-1)"); - MH3 mh3i("MMcEvt.fImpact/100", "(MEnergyEst.fImpact/MMcEvt.fImpact-1)*(MEnergyEst.fImpact/MMcEvt.fImpact-1)"); - MH3 mh3eo("MMcEvt.fEnergy", "MEnergyEst.fEnergy/MMcEvt.fEnergy-1"); - MH3 mh3io("MMcEvt.fImpact/100", "MEnergyEst.fImpact/MMcEvt.fImpact-1"); - - MH3 mh3e2("MEnergyEst.fEnergy", "(MEnergyEst.fEnergy/MMcEvt.fEnergy-1)*(MEnergyEst.fEnergy/MMcEvt.fEnergy-1)"); - MH3 mh3i2("MEnergyEst.fImpact/100", "(MEnergyEst.fImpact/MMcEvt.fImpact-1)*(MEnergyEst.fImpact/MMcEvt.fImpact-1)"); - MH3 mh3eo2("MEnergyEst.fEnergy", "MEnergyEst.fEnergy/MMcEvt.fEnergy-1"); - MH3 mh3io2("MEnergyEst.fImpact/100", "MEnergyEst.fImpact/MMcEvt.fImpact-1"); - - MH3 mhe("MMcEvt.fEnergy", "MEnergyEst.fEnergy"); - MH3 mhi("MMcEvt.fImpact/100", "MEnergyEst.fImpact/100"); - - mh3e.SetName("HistEnergy"); - mh3i.SetName("HistImpact"); - mh3eo.SetName("HistEnergyOffset"); - mh3io.SetName("HistImpactOffset"); - - mh3e2.SetName("HistEnergy"); - mh3i2.SetName("HistImpact"); - mh3eo2.SetName("HistEnergyOffset"); - mh3io2.SetName("HistImpactOffset"); + 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"); - mhi.SetName("HistII"); - - MFillH hfille(&mh3e); - MFillH hfilli(&mh3i); - MFillH hfilleo(&mh3eo); - MFillH hfillio(&mh3io); - - MFillH hfille2(&mh3e2); - MFillH hfilli2(&mh3i2); - MFillH hfilleo2(&mh3eo2); - MFillH hfillio2(&mh3io2); - - MFillH hfillee(&mhe); - MFillH hfillii(&mhi); - - MBinning binsex("BinningHistEnergyX"); - MBinning binsey("BinningHistEnergyY"); - MBinning binsix("BinningHistImpactX"); - MBinning binsiy("BinningHistImpactY"); - MBinning binseox("BinningHistEnergyOffsetX"); - MBinning binseoy("BinningHistEnergyOffsetY"); - MBinning binsiox("BinningHistImpactOffsetX"); - MBinning binsioy("BinningHistImpactOffsetY"); + mh3ed.SetName("HistEnergyDelta"); + mh3ed2.SetName("HistEnergyDelta"); + + MBinning binsedx("BinningHistEnergyDeltaX"); + MBinning binsedy("BinningHistEnergyDeltaY"); MBinning binseex("BinningHistEEX"); - MBinning binsiix("BinningHistIIX"); MBinning binseey("BinningHistEEY"); - MBinning binsiiy("BinningHistIIY"); - - binsex.SetEdgesLog(50, usect1 ? 300: 10, usect1 ? 50000 : 1e4); - binsey.SetEdges(50, 0, usect1 ? 0.8 : 1.75); - binseox.SetEdgesLog(50, usect1 ? 300 : 10, usect1 ? 50000 : 1e4); - binseoy.SetEdges(50, usect1 ? -0.75 : -1.75, usect1 ? 0.75 : 1.75); - - binsix.SetEdges(50, 0, usect1 ? 275 : 300); - binsiy.SetEdges(50, 0, usect1 ? 0.2 : 1.75); - binsiox.SetEdges(50, 0, usect1 ? 275 : 300); - binsioy.SetEdges(50, usect1 ? -0.75 : -1.75, usect1 ? 0.75 : 1.75); - - binseex.SetEdgesLog(50, usect1 ? 300 : 10, usect1 ? 50000 : 15e3); - binseey.SetEdgesLog(50, usect1 ? 300 : 1, usect1 ? 50000 : 2e3); - binsiix.SetEdges(50, 0, usect1 ? 275 : 300); - binsiiy.SetEdges(50, 0, usect1 ? 275 : 150); - - parlist.AddToList(&binsex); - parlist.AddToList(&binsey); - parlist.AddToList(&binsix); - parlist.AddToList(&binsiy); - parlist.AddToList(&binseox); - parlist.AddToList(&binseoy); - parlist.AddToList(&binsiox); - parlist.AddToList(&binsioy); - parlist.AddToList(&binseex); - parlist.AddToList(&binseey); - parlist.AddToList(&binsiix); - parlist.AddToList(&binsiiy); - - eest.StopMapping(); - eest.Add("HillasSrc"); + + 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(&eest); + tlist2.AddToList(&eest2); + // tlist2.AddToList(new MPrint("MMcEvt")); // tlist2.AddToList(new MPrint("MEnergyEst")); - tlist2.AddToList(&hfille); - tlist2.AddToList(&hfilli); - tlist2.AddToList(&hfilleo); - tlist2.AddToList(&hfillio); - - tlist2.AddToList(&hfille2); - tlist2.AddToList(&hfilli2); - tlist2.AddToList(&hfilleo2); - tlist2.AddToList(&hfillio2); - + tlist2.AddToList(&hfilled); + tlist2.AddToList(&hfilled2); tlist2.AddToList(&hfillee); - tlist2.AddToList(&hfillii); // @@ -357,5 +379,5 @@ MEvtLoop evtloop2; evtloop2.SetProgressBar(&bar); - evtloop2.SetParList(&parlist); + evtloop2.SetParList(&parlist2); // @@ -365,39 +387,53 @@ return; - tlist2.PrintStatistics(); - - const TString text = "\\sqrt{}=%.0f%%"; - - char txt[1000]; - - TCanvas *c=new TCanvas("Est1", "Estimates vs. E_{true}"); + // + // 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); - mh3i.DrawClone("PROFXnonew"); - sprintf(txt, text.Data(), sqrt(mh3i.GetHist().GetMean(2))*100); - TLatex *t = new TLatex(180, 0.15, txt); - t->Draw(); + 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); - mh3e.DrawClone("PROFXnonew"); - sprintf(txt, text.Data(), sqrt(mh3e.GetHist().GetMean(2))*100); - t = new TLatex(3.5, 0.6, txt); - t->Draw(); + 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); - mh3io.DrawClone("PROFXnonew"); + energy_3->SetLeftMargin(0.15); + energy_3->SetBottomMargin(0.12); + mh3ed2.DrawClone("PROFXnonew"); + + ((TH2*)mh3ed2.GetHist())->ProjectionY("deltae"); + c->cd(4); - mh3eo.DrawClone("PROFXnonew"); - - c=new TCanvas("Est2", "Estimates vs. E_{est}"); - c->Divide(2,2); - c->cd(1); - mh3i2.DrawClone("PROFXnonew"); - c->cd(2); - mh3e2.DrawClone("PROFXnonew"); - c->cd(3); - mh3io2.DrawClone("PROFXnonew"); - c->cd(4); - mh3eo2.DrawClone("PROFXnonew"); - - mhe.DrawClone("PROFX"); - mhi.DrawClone("PROFX"); + 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/manalysis/MEnergyEstParam.cc =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MEnergyEstParam.cc (revision 1885) +++ /trunk/MagicSoft/Mars/manalysis/MEnergyEstParam.cc (revision 1886) @@ -58,11 +58,12 @@ fA[3] = 0.006; fA[4] = 58.3; + /* - fA[0] = 39.667402; // [cm] - fA[1] = 143.081912; // [cm] - fA[2] = -395.501677; // [cm] - fA[3] = 0.006193; - fA[4] = 0; - */ + fA[0] = 39.667402; // [cm] + fA[1] = 143.081912; // [cm] + fA[2] = -395.501677; // [cm] + fA[3] = 0.006193; + fA[4] = 0; + */ fB.Set(7); @@ -74,13 +75,14 @@ fB[5] = 4.13e-8; // [GeV] fB[6] = 0.05; + /* - fB[0] = 45.037867; // [GeV] - fB[1] = 0.087222; // [GeV] - fB[2] = -0.208065; // [GeV/cm] - fB[3] = 78.164942; // [GeV] - fB[4] = -159.361283; // [GeV] - fB[5] = -0.130455; // [GeV/cm] - fB[6] = 0.051139; - */ + fB[0] = 45.037867; // [GeV] + fB[1] = 0.087222; // [GeV] + fB[2] = -0.208065; // [GeV/cm] + fB[3] = 78.164942; // [GeV] + fB[4] = -159.361283; // [GeV] + fB[5] = -0.130455; // [GeV/cm] + fB[6] = 0.051139; + */ } @@ -340,11 +342,9 @@ const Double_t dist = fMatrix ? GetVal(col++) : ((MHillasSrc*)NextH())->GetDist(); - /* DANIEL - MARCOS */ const Double_t ir = i0 * (i1 + fA[1]*dist); // [cm] Double_t er = e0 * (e1 + e2*ir); // [GeV] - - /* MY PARAM */ + /* THOMAS BRETZ */ // if (width==0) return kCONTINUE; // const Double_t ir = i0 * (i1 + dist*(fA[1]/width + fA[4]/log10(size))); // [cm] Index: /trunk/MagicSoft/Mars/manalysis/MEnergyEstParam.h =================================================================== --- /trunk/MagicSoft/Mars/manalysis/MEnergyEstParam.h (revision 1885) +++ /trunk/MagicSoft/Mars/manalysis/MEnergyEstParam.h (revision 1886) @@ -40,4 +40,5 @@ Double_t GetVal(Int_t i) const; + public: MEnergyEstParam(const char *hil="MHillas", const char *name=NULL, const char *title=NULL); @@ -58,4 +59,6 @@ void SetCoeffB(const TArrayD &arr); + Double_t GetCoeff(Int_t i) { return i