/* ======================================================================== *\ ! ! * ! * 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"); }