Index: trunk/MagicSoft/Mars/macros/flux.C =================================================================== --- trunk/MagicSoft/Mars/macros/flux.C (revision 1211) +++ trunk/MagicSoft/Mars/macros/flux.C (revision 1211) @@ -0,0 +1,273 @@ +/* ======================================================================== *\ +! +! * +! * 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 1/2002 +! Author(s): Wolfgang Wittek 1/2002 +! +! Copyright: MAGIC Software Development, 2000-2002 +! +! +\* ======================================================================== */ + +void flux() +{ + // + // first we have to create our empty lists + // + MParList parlist; + MTaskList tasklist; + + parlist.AddToList(&tasklist); + + // + // Define the two source positions + // + MSrcPosCam source("Source"); + MSrcPosCam antisrc("AntiSource"); + source.SetXY(0, 0); + antisrc.SetXY(+240, 0); + parlist.AddToList(&source); + parlist.AddToList(&antisrc); + + // + // Setup binning for the histograms + // + MBinning binsalpha("BinningAlpha"); + binsalpha.SetEdges(45, -90, 90); + + MBinning binse("BinningE"); + binse.SetEdgesLog(10, 10, 1e3); + + MBinning binstheta("BinningTheta"); + binstheta.SetEdges(5, 2.5, 27.5); + + MBinning binstime("BinningTime"); + binstime.SetEdges(5, 0, 50); + + MBinning binsdifftime("BinningTimeDiff"); + binsdifftime.SetEdges(50, 0, 0.1); + + parlist.AddToList(&binsalpha); + parlist.AddToList(&binse); + parlist.AddToList(&binstheta); + parlist.AddToList(&binstime); + parlist.AddToList(&binsdifftime); + + // + // Setup our tasks + // + MReadMarsFile reader("Events", "~/data/Gamma*.root"); + + MGeomCamMagic geomcam; + parlist.AddToList(&geomcam); + + MMcTimeGenerate rand; + MMcPedestalCopy pcopy; + MMcPedestalNSBAdd pnsb; + MCerPhotCalc ncalc; + MImgCleanStd clean; + MBlindPixelCalc blind; + MHillasCalc hcalc; + MHillasSrcCalc hsrc1("Source", "HillasSrc"); + MHillasSrcCalc hsrc2("AntiSource", "HillasAntiSrc"); + MEnergyEstimate estim; + + MFillH hfill1h("MHHillas", "MHillas"); + MFillH hfill1m("MHStarMap", "MHillas"); + MFillH hfill2s("HSource [MHHillasSrc]", "HillasSrc"); + MFillH hfill2a("HAntiSource [MHHillasSrc]", "HillasAntiSrc"); + + MMcCollectionAreaCalc acalc; + + const Float_t alpha0 = 15; // [deg] + + MFAlpha fsrc ("HillasSrc", '>', alpha0); + MFAlpha fasrc("HillasAntiSrc", '>', alpha0); + + MFillH fillsp ("FluxSrcTime [MHAlphaEnergyTime]", "HillasSrc"); + MFillH fillasp ("FluxASrcTime [MHAlphaEnergyTime]", "HillasAntiSrc"); + MFillH fillsptheta ("FluxSrcTheta [MHAlphaEnergyTheta]", "HillasSrc"); + MFillH fillasptheta("FluxASrcTheta [MHAlphaEnergyTheta]", "HillasAntiSrc"); + + MFTriggerLvl1 lvl1; + MFillH fethetaall("AllTheta [MHEnergyTheta]", "MMcEvt"); + MFillH fethetasel("SelTheta [MHEnergyTheta]", "MMcEvt"); + MFillH fetimeall ("AllTime [MHEnergyTime]", "MMcEvt"); + MFillH fetimesel ("SelTime [MHEnergyTime]", "MMcEvt"); + + fetimesel.SetFilter(&lvl1); + fethetasel.SetFilter(&lvl1); + + fillsp.SetFilter(&fasrc); + fillasp.SetFilter(&fsrc); + fillsptheta.SetFilter(&fsrc); + fillasptheta.SetFilter(&fasrc); + + MFillH fillontime ("EffOnTime [MHTimeDiffTime]", "MMcEvt"); + MFillH fillontheta("EffOnTheta [MHTimeDiffTheta]", "MMcEvt"); + + // + // Setup Task list + // + tasklist.AddToList(&reader); + tasklist.AddToList(&rand); + tasklist.AddToList(&fillontime); + tasklist.AddToList(&fillontheta); + tasklist.AddToList(&acalc); + tasklist.AddToList(&pcopy); + tasklist.AddToList(&pnsb); + tasklist.AddToList(&ncalc); + tasklist.AddToList(&clean); + tasklist.AddToList(&blind); + tasklist.AddToList(&hcalc); + tasklist.AddToList(&hsrc1); + tasklist.AddToList(&hsrc2); + tasklist.AddToList(&estim); + tasklist.AddToList(&hfill1h); + tasklist.AddToList(&hfill1m); + tasklist.AddToList(&hfill2s); + tasklist.AddToList(&hfill2a); + tasklist.AddToList(&fsrc); + tasklist.AddToList(&fasrc); + tasklist.AddToList(&fillsp); + tasklist.AddToList(&fillasptheta); + tasklist.AddToList(&fillasp); + tasklist.AddToList(&fillsptheta); + tasklist.AddToList(&lvl1); + tasklist.AddToList(&fethetaall); + tasklist.AddToList(&fethetasel); + tasklist.AddToList(&fetimeall); + tasklist.AddToList(&fetimesel); + + // + // set up the loop for the processing + // + MEvtLoop magic; + magic.SetParList(&parlist); + + // + // Start to loop over all events + // + if (!magic.Eventloop()) + return; + + tasklist.PrintStatistics(); + + return; + + parlist.FindObject("MHMcCollectionArea")->DrawClone(); + parlist.FindObject("MHStarMap")->DrawClone(); + + // ------------ Eff On Time ----------------- + + MHEnergyTime &alltime = *(MHEnergyTime*)parlist.FindObject("AllTime"); + MHEnergyTheta &alltheta = *(MHEnergyTheta*)parlist.FindObject("AllTheta"); + MHEnergyTime &seltime = *(MHEnergyTime*)parlist.FindObject("SelTime"); + MHEnergyTheta &seltheta = *(MHEnergyTheta*)parlist.FindObject("SelTheta"); + + MHEnergyTime collareatime; + MHEnergyTheta collareatheta; + collareatime.Divide(&seltime, &alltime); + collareatheta.Divide(&seltheta, &alltheta); + + MHTimeDiffTime &effontime = *(MHTimeDiffTime*)parlist.FindObject("EffOnTime"); + MHTimeDiffTheta &effontheta = *(MHTimeDiffTheta*)parlist.FindObject("EffOnTheta"); + + effontime.DrawClone(); + effontheta.DrawClone(); + + MHEffOnTimeTime ontime; + MHEffOnTimeTheta ontheta; + ontime.SetupFill(&parlist); + ontheta.SetupFill(&parlist); + + ontime.Calc(effontime.GetHist()); + ontheta.Calc(effontheta.GetHist()); + + ontime.DrawClone(); + ontheta.DrawClone(); + + parlist.FindObject("HSource")->DrawClone();; + parlist.FindObject("HAntiSource")->DrawClone(); + parlist.FindObject("MHHillas")->DrawClone(); + + MHAlphaEnergyTime &fluxsp = *(MHAlphaEnergyTime)parlist.FindObject("HillasSrc"); + MHAlphaEnergyTime &fluxasp = *(MHAlphaEnergyTime)parlist.FindObject("HillasAntiSrc"); + MHAlphaEnergyTheta &fluxsptheta = *(MHAlphaEnergyTime)parlist.FindObject("HillasSrc"); + MHAlphaEnergyTheta &fluxasptheta = *(MHAlphaEnergyTime)parlist.FindObject("HillasAntiSrc"); + + fluxsp.DrawClone(); + fluxasp.DrawClone(); + fluxsptheta.DrawClone(); + fluxasptheta.DrawClone(); + + MHAlphaEnergyTheta resulttime; + MHAlphaEnergyTheta resulttheta; + resulttime.Substract(&fluxsp, &fluxasp); + resulttheta.Substract(&fluxsptheta, &fluxasptheta); + + resulttime.DrawClone(); + resulttheta.DrawClone(); + + TH2D &projecttime = *resulttime.GetAlphaProjection(-10, 10); + TH2D &projecttheta = *resulttheta.GetAlphaProjection(-10, 10); + + projecttime.SetTitle("Number of Gammas vs. EnergyEst and Time (Alpha integrated between -10, 10deg)"); + projecttheta.SetTitle("Number of Gammas vs. EnergyEst and Theta (Alpha integrated between -10, 10deg)"); + + c = new TCanvas("To be unfolded"); + c->Divide(2,2); + c->cd(1); + projecttime.DrawCopy(); + c->cd(2); + projecttheta.DrawCopy(); + + for (int i=1; i<=binstime.GetNumBins(); i++) + { + if (ontime.GetHist()->GetBinContent(i)==0) + continue; + + TH1D &hist = *projecttime.ProjectionX("Number of Gammas vs. EnergyEst for a fixed time", i, i); + + /* UNFOLDING */ + + //hist->Divide(collareatime); + hist.Scale(1./ontime.GetHist()->GetBinContent(i)); + + for (int j=1; j<=binse.GetNumBins(); j++) + hist.SetBinContent(j, hist.GetBinContent(j)/hist.GetBinWidth(j)); + + hist.SetName("Flux"); + hist.SetTitle("Flux[Gammas/s/m^2/GeV] vs. EnergyTrue for a fixed Time"); + + char n[100]; + sprintf(n, "Canv%d", j); + c= new TCanvas(n, "Title"); + hist.DrawCopy(); + } + + delete &projecttime; + delete &projecttheta; + + return; + + // ------------------------------------------ + + MHMcCollectionArea carea; + TH1D *collareatime = carea.GetHist(); // FIXME! + TH1D *collareatheta = carea.GetHist(); // FIXME! +}