| 1 | /* ======================================================================== *\
|
|---|
| 2 | !
|
|---|
| 3 | ! *
|
|---|
| 4 | ! * This file is part of MARS, the MAGIC Analysis and Reconstruction
|
|---|
| 5 | ! * Software. It is distributed to you in the hope that it can be a useful
|
|---|
| 6 | ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
|
|---|
| 7 | ! * It is distributed WITHOUT ANY WARRANTY.
|
|---|
| 8 | ! *
|
|---|
| 9 | ! * Permission to use, copy, modify and distribute this software and its
|
|---|
| 10 | ! * documentation for any purpose is hereby granted without fee,
|
|---|
| 11 | ! * provided that the above copyright notice appear in all copies and
|
|---|
| 12 | ! * that both that copyright notice and this permission notice appear
|
|---|
| 13 | ! * in supporting documentation. It is provided "as is" without express
|
|---|
| 14 | ! * or implied warranty.
|
|---|
| 15 | ! *
|
|---|
| 16 | !
|
|---|
| 17 | !
|
|---|
| 18 | ! Author(s): Thomas Bretz 1/2002 <mailto:tbretz@uni-sw.gwdg.de>
|
|---|
| 19 | ! Author(s): Wolfgang Wittek 1/2002 <mailto:wittek@mppmu.mpg.de>
|
|---|
| 20 | !
|
|---|
| 21 | ! Copyright: MAGIC Software Development, 2000-2002
|
|---|
| 22 | !
|
|---|
| 23 | !
|
|---|
| 24 | \* ======================================================================== */
|
|---|
| 25 |
|
|---|
| 26 | void flux()
|
|---|
| 27 | {
|
|---|
| 28 | //
|
|---|
| 29 | // first we have to create our empty lists
|
|---|
| 30 | //
|
|---|
| 31 | MParList parlist;
|
|---|
| 32 | MTaskList tasklist;
|
|---|
| 33 |
|
|---|
| 34 | parlist.AddToList(&tasklist);
|
|---|
| 35 |
|
|---|
| 36 | //
|
|---|
| 37 | // Define the two source positions
|
|---|
| 38 | //
|
|---|
| 39 | MSrcPosCam source("Source");
|
|---|
| 40 | MSrcPosCam antisrc("AntiSource");
|
|---|
| 41 | source.SetXY(0, 0);
|
|---|
| 42 | antisrc.SetXY(+240, 0);
|
|---|
| 43 | parlist.AddToList(&source);
|
|---|
| 44 | parlist.AddToList(&antisrc);
|
|---|
| 45 |
|
|---|
| 46 | //
|
|---|
| 47 | // Setup binning for the histograms
|
|---|
| 48 | //
|
|---|
| 49 | MBinning binsalpha("BinningAlpha");
|
|---|
| 50 | binsalpha.SetEdges(45, -90, 90);
|
|---|
| 51 |
|
|---|
| 52 | MBinning binse("BinningE");
|
|---|
| 53 | binse.SetEdgesLog(10, 10, 1e3);
|
|---|
| 54 |
|
|---|
| 55 | MBinning binstheta("BinningTheta");
|
|---|
| 56 | binstheta.SetEdges(5, 2.5, 27.5);
|
|---|
| 57 |
|
|---|
| 58 | MBinning binstime("BinningTime");
|
|---|
| 59 | binstime.SetEdges(5, 0, 50);
|
|---|
| 60 |
|
|---|
| 61 | MBinning binsdifftime("BinningTimeDiff");
|
|---|
| 62 | binsdifftime.SetEdges(50, 0, 0.1);
|
|---|
| 63 |
|
|---|
| 64 | parlist.AddToList(&binsalpha);
|
|---|
| 65 | parlist.AddToList(&binse);
|
|---|
| 66 | parlist.AddToList(&binstheta);
|
|---|
| 67 | parlist.AddToList(&binstime);
|
|---|
| 68 | parlist.AddToList(&binsdifftime);
|
|---|
| 69 |
|
|---|
| 70 | //
|
|---|
| 71 | // Setup our tasks
|
|---|
| 72 | //
|
|---|
| 73 | MReadMarsFile reader("Events", "~/data/Gamma*.root");
|
|---|
| 74 |
|
|---|
| 75 | MGeomCamMagic geomcam;
|
|---|
| 76 | parlist.AddToList(&geomcam);
|
|---|
| 77 |
|
|---|
| 78 | MMcTimeGenerate rand;
|
|---|
| 79 | MMcPedestalCopy pcopy;
|
|---|
| 80 | MMcPedestalNSBAdd pnsb;
|
|---|
| 81 | MCerPhotCalc ncalc;
|
|---|
| 82 | MImgCleanStd clean;
|
|---|
| 83 | MBlindPixelCalc blind;
|
|---|
| 84 | MHillasCalc hcalc;
|
|---|
| 85 | MHillasSrcCalc hsrc1("Source", "HillasSrc");
|
|---|
| 86 | MHillasSrcCalc hsrc2("AntiSource", "HillasAntiSrc");
|
|---|
| 87 | MEnergyEstimate estim;
|
|---|
| 88 |
|
|---|
| 89 | MFillH hfill1h("MHHillas", "MHillas");
|
|---|
| 90 | MFillH hfill1m("MHStarMap", "MHillas");
|
|---|
| 91 | MFillH hfill2s("HSource [MHHillasSrc]", "HillasSrc");
|
|---|
| 92 | MFillH hfill2a("HAntiSource [MHHillasSrc]", "HillasAntiSrc");
|
|---|
| 93 |
|
|---|
| 94 | MMcCollectionAreaCalc acalc;
|
|---|
| 95 |
|
|---|
| 96 | const Float_t alpha0 = 15; // [deg]
|
|---|
| 97 |
|
|---|
| 98 | MFAlpha fsrc ("HillasSrc", '>', alpha0);
|
|---|
| 99 | MFAlpha fasrc("HillasAntiSrc", '>', alpha0);
|
|---|
| 100 |
|
|---|
| 101 | MFillH fillsp ("FluxSrcTime [MHAlphaEnergyTime]", "HillasSrc");
|
|---|
| 102 | MFillH fillasp ("FluxASrcTime [MHAlphaEnergyTime]", "HillasAntiSrc");
|
|---|
| 103 | MFillH fillsptheta ("FluxSrcTheta [MHAlphaEnergyTheta]", "HillasSrc");
|
|---|
| 104 | MFillH fillasptheta("FluxASrcTheta [MHAlphaEnergyTheta]", "HillasAntiSrc");
|
|---|
| 105 |
|
|---|
| 106 | MFTriggerLvl1 lvl1;
|
|---|
| 107 | MFillH fethetaall("AllTheta [MHEnergyTheta]", "MMcEvt");
|
|---|
| 108 | MFillH fethetasel("SelTheta [MHEnergyTheta]", "MMcEvt");
|
|---|
| 109 | MFillH fetimeall ("AllTime [MHEnergyTime]", "MMcEvt");
|
|---|
| 110 | MFillH fetimesel ("SelTime [MHEnergyTime]", "MMcEvt");
|
|---|
| 111 |
|
|---|
| 112 | fetimesel.SetFilter(&lvl1);
|
|---|
| 113 | fethetasel.SetFilter(&lvl1);
|
|---|
| 114 |
|
|---|
| 115 | fillsp.SetFilter(&fasrc);
|
|---|
| 116 | fillasp.SetFilter(&fsrc);
|
|---|
| 117 | fillsptheta.SetFilter(&fsrc);
|
|---|
| 118 | fillasptheta.SetFilter(&fasrc);
|
|---|
| 119 |
|
|---|
| 120 | MFillH fillontime ("EffOnTime [MHTimeDiffTime]", "MMcEvt");
|
|---|
| 121 | MFillH fillontheta("EffOnTheta [MHTimeDiffTheta]", "MMcEvt");
|
|---|
| 122 |
|
|---|
| 123 | //
|
|---|
| 124 | // Setup Task list
|
|---|
| 125 | //
|
|---|
| 126 | tasklist.AddToList(&reader);
|
|---|
| 127 | tasklist.AddToList(&rand);
|
|---|
| 128 | tasklist.AddToList(&fillontime);
|
|---|
| 129 | tasklist.AddToList(&fillontheta);
|
|---|
| 130 | tasklist.AddToList(&acalc);
|
|---|
| 131 | tasklist.AddToList(&pcopy);
|
|---|
| 132 | tasklist.AddToList(&pnsb);
|
|---|
| 133 | tasklist.AddToList(&ncalc);
|
|---|
| 134 | tasklist.AddToList(&clean);
|
|---|
| 135 | tasklist.AddToList(&blind);
|
|---|
| 136 | tasklist.AddToList(&hcalc);
|
|---|
| 137 | tasklist.AddToList(&hsrc1);
|
|---|
| 138 | tasklist.AddToList(&hsrc2);
|
|---|
| 139 | tasklist.AddToList(&estim);
|
|---|
| 140 | tasklist.AddToList(&hfill1h);
|
|---|
| 141 | tasklist.AddToList(&hfill1m);
|
|---|
| 142 | tasklist.AddToList(&hfill2s);
|
|---|
| 143 | tasklist.AddToList(&hfill2a);
|
|---|
| 144 | tasklist.AddToList(&fsrc);
|
|---|
| 145 | tasklist.AddToList(&fasrc);
|
|---|
| 146 | tasklist.AddToList(&fillsp);
|
|---|
| 147 | tasklist.AddToList(&fillasptheta);
|
|---|
| 148 | tasklist.AddToList(&fillasp);
|
|---|
| 149 | tasklist.AddToList(&fillsptheta);
|
|---|
| 150 | tasklist.AddToList(&lvl1);
|
|---|
| 151 | tasklist.AddToList(&fethetaall);
|
|---|
| 152 | tasklist.AddToList(&fethetasel);
|
|---|
| 153 | tasklist.AddToList(&fetimeall);
|
|---|
| 154 | tasklist.AddToList(&fetimesel);
|
|---|
| 155 |
|
|---|
| 156 | //
|
|---|
| 157 | // set up the loop for the processing
|
|---|
| 158 | //
|
|---|
| 159 | MEvtLoop magic;
|
|---|
| 160 | magic.SetParList(&parlist);
|
|---|
| 161 |
|
|---|
| 162 | //
|
|---|
| 163 | // Start to loop over all events
|
|---|
| 164 | //
|
|---|
| 165 | if (!magic.Eventloop())
|
|---|
| 166 | return;
|
|---|
| 167 |
|
|---|
| 168 | tasklist.PrintStatistics();
|
|---|
| 169 |
|
|---|
| 170 | return;
|
|---|
| 171 |
|
|---|
| 172 | parlist.FindObject("MHMcCollectionArea")->DrawClone();
|
|---|
| 173 | parlist.FindObject("MHStarMap")->DrawClone();
|
|---|
| 174 |
|
|---|
| 175 | // ------------ Eff On Time -----------------
|
|---|
| 176 |
|
|---|
| 177 | MHEnergyTime &alltime = *(MHEnergyTime*)parlist.FindObject("AllTime");
|
|---|
| 178 | MHEnergyTheta &alltheta = *(MHEnergyTheta*)parlist.FindObject("AllTheta");
|
|---|
| 179 | MHEnergyTime &seltime = *(MHEnergyTime*)parlist.FindObject("SelTime");
|
|---|
| 180 | MHEnergyTheta &seltheta = *(MHEnergyTheta*)parlist.FindObject("SelTheta");
|
|---|
| 181 |
|
|---|
| 182 | MHEnergyTime collareatime;
|
|---|
| 183 | MHEnergyTheta collareatheta;
|
|---|
| 184 | collareatime.Divide(&seltime, &alltime);
|
|---|
| 185 | collareatheta.Divide(&seltheta, &alltheta);
|
|---|
| 186 |
|
|---|
| 187 | MHTimeDiffTime &effontime = *(MHTimeDiffTime*)parlist.FindObject("EffOnTime");
|
|---|
| 188 | MHTimeDiffTheta &effontheta = *(MHTimeDiffTheta*)parlist.FindObject("EffOnTheta");
|
|---|
| 189 |
|
|---|
| 190 | effontime.DrawClone();
|
|---|
| 191 | effontheta.DrawClone();
|
|---|
| 192 |
|
|---|
| 193 | MHEffOnTimeTime ontime;
|
|---|
| 194 | MHEffOnTimeTheta ontheta;
|
|---|
| 195 | ontime.SetupFill(&parlist);
|
|---|
| 196 | ontheta.SetupFill(&parlist);
|
|---|
| 197 |
|
|---|
| 198 | ontime.Calc(effontime.GetHist());
|
|---|
| 199 | ontheta.Calc(effontheta.GetHist());
|
|---|
| 200 |
|
|---|
| 201 | ontime.DrawClone();
|
|---|
| 202 | ontheta.DrawClone();
|
|---|
| 203 |
|
|---|
| 204 | parlist.FindObject("HSource")->DrawClone();;
|
|---|
| 205 | parlist.FindObject("HAntiSource")->DrawClone();
|
|---|
| 206 | parlist.FindObject("MHHillas")->DrawClone();
|
|---|
| 207 |
|
|---|
| 208 | MHAlphaEnergyTime &fluxsp = *(MHAlphaEnergyTime)parlist.FindObject("HillasSrc");
|
|---|
| 209 | MHAlphaEnergyTime &fluxasp = *(MHAlphaEnergyTime)parlist.FindObject("HillasAntiSrc");
|
|---|
| 210 | MHAlphaEnergyTheta &fluxsptheta = *(MHAlphaEnergyTime)parlist.FindObject("HillasSrc");
|
|---|
| 211 | MHAlphaEnergyTheta &fluxasptheta = *(MHAlphaEnergyTime)parlist.FindObject("HillasAntiSrc");
|
|---|
| 212 |
|
|---|
| 213 | fluxsp.DrawClone();
|
|---|
| 214 | fluxasp.DrawClone();
|
|---|
| 215 | fluxsptheta.DrawClone();
|
|---|
| 216 | fluxasptheta.DrawClone();
|
|---|
| 217 |
|
|---|
| 218 | MHAlphaEnergyTheta resulttime;
|
|---|
| 219 | MHAlphaEnergyTheta resulttheta;
|
|---|
| 220 | resulttime.Substract(&fluxsp, &fluxasp);
|
|---|
| 221 | resulttheta.Substract(&fluxsptheta, &fluxasptheta);
|
|---|
| 222 |
|
|---|
| 223 | resulttime.DrawClone();
|
|---|
| 224 | resulttheta.DrawClone();
|
|---|
| 225 |
|
|---|
| 226 | TH2D &projecttime = *resulttime.GetAlphaProjection(-10, 10);
|
|---|
| 227 | TH2D &projecttheta = *resulttheta.GetAlphaProjection(-10, 10);
|
|---|
| 228 |
|
|---|
| 229 | projecttime.SetTitle("Number of Gammas vs. EnergyEst and Time (Alpha integrated between -10, 10deg)");
|
|---|
| 230 | projecttheta.SetTitle("Number of Gammas vs. EnergyEst and Theta (Alpha integrated between -10, 10deg)");
|
|---|
| 231 |
|
|---|
| 232 | c = new TCanvas("To be unfolded");
|
|---|
| 233 | c->Divide(2,2);
|
|---|
| 234 | c->cd(1);
|
|---|
| 235 | projecttime.DrawCopy();
|
|---|
| 236 | c->cd(2);
|
|---|
| 237 | projecttheta.DrawCopy();
|
|---|
| 238 |
|
|---|
| 239 | for (int i=1; i<=binstime.GetNumBins(); i++)
|
|---|
| 240 | {
|
|---|
| 241 | if (ontime.GetHist()->GetBinContent(i)==0)
|
|---|
| 242 | continue;
|
|---|
| 243 |
|
|---|
| 244 | TH1D &hist = *projecttime.ProjectionX("Number of Gammas vs. EnergyEst for a fixed time", i, i);
|
|---|
| 245 |
|
|---|
| 246 | /* UNFOLDING */
|
|---|
| 247 |
|
|---|
| 248 | //hist->Divide(collareatime);
|
|---|
| 249 | hist.Scale(1./ontime.GetHist()->GetBinContent(i));
|
|---|
| 250 |
|
|---|
| 251 | for (int j=1; j<=binse.GetNumBins(); j++)
|
|---|
| 252 | hist.SetBinContent(j, hist.GetBinContent(j)/hist.GetBinWidth(j));
|
|---|
| 253 |
|
|---|
| 254 | hist.SetName("Flux");
|
|---|
| 255 | hist.SetTitle("Flux[Gammas/s/m^2/GeV] vs. EnergyTrue for a fixed Time");
|
|---|
| 256 |
|
|---|
| 257 | char n[100];
|
|---|
| 258 | sprintf(n, "Canv%d", j);
|
|---|
| 259 | c= new TCanvas(n, "Title");
|
|---|
| 260 | hist.DrawCopy();
|
|---|
| 261 | }
|
|---|
| 262 |
|
|---|
| 263 | delete &projecttime;
|
|---|
| 264 | delete &projecttheta;
|
|---|
| 265 |
|
|---|
| 266 | return;
|
|---|
| 267 |
|
|---|
| 268 | // ------------------------------------------
|
|---|
| 269 |
|
|---|
| 270 | MHMcCollectionArea carea;
|
|---|
| 271 | TH1D *collareatime = carea.GetHist(); // FIXME!
|
|---|
| 272 | TH1D *collareatheta = carea.GetHist(); // FIXME!
|
|---|
| 273 | }
|
|---|
