Index: /trunk/MagicSoft/Mars/macros/estct1.C =================================================================== --- /trunk/MagicSoft/Mars/macros/estct1.C (revision 1843) +++ /trunk/MagicSoft/Mars/macros/estct1.C (revision 1843) @@ -0,0 +1,243 @@ +/* ======================================================================== *\ +! +! * +! * 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 et al, 09/2002 +! +! Copyright: MAGIC Software Development, 2000-2002 +! +! +\* ======================================================================== */ + + +void estct1() +{ + // + // This is a demonstration program which calculates the Hillas + // parameter out of a Magic root file (raw data file). + // + + // + // Create a empty Parameter List and an empty Task List + // The tasklist is identified in the eventloop by its name + // + MParList plist; + + + MTaskList tlist; + plist.AddToList(&tlist); + + // + // Now setup the tasks and tasklist: + // --------------------------------- + // + MReadTree read("Events", "~/ct1/MC_ON2.root"); + + //MReadMarsFile read("Events"); + read.DisableAutoScheme(); + /* + read.AddFile("star.root"); + read.AddFile("star2.root"); + */ + //read.AddFile("~/ct1/MC_ON2.root"); + + MEnergyEstParam eest("Hillas"); + eest.Add("HillasSrc"); + + // + // Use this to change the binnign of the histograms to CT1-style + // + Bool_t usect1 = kTRUE; + + // + // Here you set up the coefficients of the parametrization + // (MEnergyEstParam) + // + TArrayD fA(5); + fA[0] = -3907.74; //4916.4; //-2414.75; + fA[1] = 1162.3; //149.549; // 1134.28; + fA[2] = 199.351;//-558.209; // 132.932; + fA[3] = 0.403192;//0.270725; //0.292845; + fA[4] = 121.921;//107.001; // 107.001; + + TArrayD fB(7); + fB[0] = 677.821;//-8234.12; //-4282.25; + fB[1] = 11.3302;//23.2153; // 18.892; + fB[2] = -0.0211177;//0.416372; //0.193373; + fB[3] = -23.0217;//332.42; //203.803; + fB[4] = -0.785242;//-0.701764; //-0.534876; + fB[5] = 5.6413e-5;//-0.0131774; //-0.00789539; + fB[6] = -0.146595;//-0.162687; // 0.111913; + + eest.SetCoeffA(fA); + eest.SetCoeffB(fB); + + 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"); + + 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"); + 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); + + plist.AddToList(&binsex); + plist.AddToList(&binsey); + plist.AddToList(&binsix); + plist.AddToList(&binsiy); + plist.AddToList(&binseox); + plist.AddToList(&binseoy); + plist.AddToList(&binsiox); + plist.AddToList(&binsioy); + plist.AddToList(&binseex); + plist.AddToList(&binseey); + plist.AddToList(&binsiix); + plist.AddToList(&binsiiy); + + // + // Setup tasklists + // + tlist.AddToList(&read); + + tlist.AddToList(&eest); + tlist.AddToList(&hfille); + tlist.AddToList(&hfilli); + tlist.AddToList(&hfilleo); + tlist.AddToList(&hfillio); + + tlist.AddToList(&hfille2); + tlist.AddToList(&hfilli2); + tlist.AddToList(&hfilleo2); + tlist.AddToList(&hfillio2); + + tlist.AddToList(&hfillee); + tlist.AddToList(&hfillii); + + /* + MPrint p("MEnergyEst"); + tlist2.AddToList(&p); + */ + + // + // Create and setup the eventloop + // + MProgressBar bar; + + MEvtLoop evtloop; + evtloop.SetProgressBar(&bar); + evtloop.SetParList(&plist); + + // + // Execute your analysis + // + if (!evtloop.Eventloop()) + return; + + tlist.PrintStatistics(); + + const TString text = "\\sqrt{}=%.0f%%"; + + char txt[1000]; + + TCanvas *c=new TCanvas("Est1", "Estimates vs. E_{true}"); + 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(); + 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(); + c->cd(3); + mh3io.DrawClone("PROFXnonew"); + 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"); +}