| 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(7, -2.5, 32.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 | reader.EnableBranch("MMcEvt.fTheta");
|
|---|
| 75 |
|
|---|
| 76 | MGeomCamMagic geomcam;
|
|---|
| 77 | parlist.AddToList(&geomcam);
|
|---|
| 78 |
|
|---|
| 79 | MMcTimeGenerate rand;
|
|---|
| 80 | MMcPedestalCopy pcopy;
|
|---|
| 81 | MMcPedestalNSBAdd pnsb;
|
|---|
| 82 | MCerPhotCalc ncalc;
|
|---|
| 83 | MImgCleanStd clean;
|
|---|
| 84 | MBlindPixelCalc blind;
|
|---|
| 85 | MHillasCalc hcalc;
|
|---|
| 86 | MHillasSrcCalc hsrc1("Source", "HillasSrc");
|
|---|
| 87 | MHillasSrcCalc hsrc2("AntiSource", "HillasAntiSrc");
|
|---|
| 88 | MEnergyEstimate estim;
|
|---|
| 89 |
|
|---|
| 90 | MFillH hfill1h("MHHillas", "MHillas");
|
|---|
| 91 | MFillH hfill1m("MHStarMap", "MHillas");
|
|---|
| 92 | MFillH hfill2s("HSource [MHHillasSrc]", "HillasSrc");
|
|---|
| 93 | MFillH hfill2a("HAntiSource [MHHillasSrc]", "HillasAntiSrc");
|
|---|
| 94 |
|
|---|
| 95 | MMcCollectionAreaCalc acalc;
|
|---|
| 96 |
|
|---|
| 97 | const Float_t alpha0 = 15; // [deg]
|
|---|
| 98 |
|
|---|
| 99 | MFAlpha fsrc ("HillasSrc", '>', alpha0);
|
|---|
| 100 | MFAlpha fasrc("HillasAntiSrc", '>', alpha0);
|
|---|
| 101 |
|
|---|
| 102 | MFillH fillsp ("FluxSrcTime [MHAlphaEnergyTime]", "HillasSrc");
|
|---|
| 103 | MFillH fillasp ("FluxASrcTime [MHAlphaEnergyTime]", "HillasAntiSrc");
|
|---|
| 104 | MFillH fillsptheta ("FluxSrcTheta [MHAlphaEnergyTheta]", "HillasSrc");
|
|---|
| 105 | MFillH fillasptheta("FluxASrcTheta [MHAlphaEnergyTheta]", "HillasAntiSrc");
|
|---|
| 106 |
|
|---|
| 107 | MFTriggerLvl1 lvl1;
|
|---|
| 108 | MFillH fethetaall("AllTheta [MHEnergyTheta]", "MMcEvt");
|
|---|
| 109 | MFillH fethetasel("SelTheta [MHEnergyTheta]", "MMcEvt");
|
|---|
| 110 | MFillH fetimeall ("AllTime [MHEnergyTime]", "MMcEvt");
|
|---|
| 111 | MFillH fetimesel ("SelTime [MHEnergyTime]", "MMcEvt");
|
|---|
| 112 |
|
|---|
| 113 | fetimesel.SetFilter(&lvl1);
|
|---|
| 114 | fethetasel.SetFilter(&lvl1);
|
|---|
| 115 |
|
|---|
| 116 | fillsp.SetFilter(&fasrc);
|
|---|
| 117 | fillasp.SetFilter(&fsrc);
|
|---|
| 118 | fillsptheta.SetFilter(&fsrc);
|
|---|
| 119 | fillasptheta.SetFilter(&fasrc);
|
|---|
| 120 |
|
|---|
| 121 | MFillH fillontime ("EffOnTime [MHTimeDiffTime]", "MMcEvt");
|
|---|
| 122 | MFillH fillontheta("EffOnTheta [MHTimeDiffTheta]", "MMcEvt");
|
|---|
| 123 |
|
|---|
| 124 | //
|
|---|
| 125 | // Setup Task list
|
|---|
| 126 | //
|
|---|
| 127 | tasklist.AddToList(&reader);
|
|---|
| 128 | tasklist.AddToList(&rand);
|
|---|
| 129 | tasklist.AddToList(&fillontime);
|
|---|
| 130 | tasklist.AddToList(&fillontheta);
|
|---|
| 131 | tasklist.AddToList(&acalc);
|
|---|
| 132 | tasklist.AddToList(&pcopy);
|
|---|
| 133 | tasklist.AddToList(&pnsb);
|
|---|
| 134 | tasklist.AddToList(&ncalc);
|
|---|
| 135 | tasklist.AddToList(&clean);
|
|---|
| 136 | tasklist.AddToList(&blind);
|
|---|
| 137 | tasklist.AddToList(&hcalc);
|
|---|
| 138 | tasklist.AddToList(&hsrc1);
|
|---|
| 139 | tasklist.AddToList(&hsrc2);
|
|---|
| 140 | tasklist.AddToList(&estim);
|
|---|
| 141 | tasklist.AddToList(&hfill1h);
|
|---|
| 142 | tasklist.AddToList(&hfill1m);
|
|---|
| 143 | tasklist.AddToList(&hfill2s);
|
|---|
| 144 | tasklist.AddToList(&hfill2a);
|
|---|
| 145 | tasklist.AddToList(&fsrc);
|
|---|
| 146 | tasklist.AddToList(&fasrc);
|
|---|
| 147 | tasklist.AddToList(&fillsp);
|
|---|
| 148 | tasklist.AddToList(&fillasptheta);
|
|---|
| 149 | tasklist.AddToList(&fillasp);
|
|---|
| 150 | tasklist.AddToList(&fillsptheta);
|
|---|
| 151 | tasklist.AddToList(&lvl1);
|
|---|
| 152 | tasklist.AddToList(&fethetaall);
|
|---|
| 153 | tasklist.AddToList(&fethetasel);
|
|---|
| 154 | tasklist.AddToList(&fetimeall);
|
|---|
| 155 | tasklist.AddToList(&fetimesel);
|
|---|
| 156 |
|
|---|
| 157 | //
|
|---|
| 158 | // set up the loop for the processing
|
|---|
| 159 | //
|
|---|
| 160 | MEvtLoop magic;
|
|---|
| 161 | magic.SetParList(&parlist);
|
|---|
| 162 |
|
|---|
| 163 | //
|
|---|
| 164 | // Start to loop over all events
|
|---|
| 165 | //
|
|---|
| 166 | if (!magic.Eventloop())
|
|---|
| 167 | return;
|
|---|
| 168 |
|
|---|
| 169 | tasklist.PrintStatistics();
|
|---|
| 170 |
|
|---|
| 171 | /*
|
|---|
| 172 | parlist.FindObject("HSource")->DrawClone();;
|
|---|
| 173 | parlist.FindObject("MHHillas")->DrawClone();
|
|---|
| 174 | parlist.FindObject("MHStarMap")->DrawClone();
|
|---|
| 175 | parlist.FindObject("HAntiSource")->DrawClone();
|
|---|
| 176 | parlist.FindObject("MHMcCollectionArea")->DrawClone();
|
|---|
| 177 | */
|
|---|
| 178 |
|
|---|
| 179 | /*
|
|---|
| 180 | MHEnergyTime &alltime = *(MHEnergyTime*)parlist.FindObject("AllTime");
|
|---|
| 181 | MHEnergyTheta &alltheta = *(MHEnergyTheta*)parlist.FindObject("AllTheta");
|
|---|
| 182 | MHEnergyTime &seltime = *(MHEnergyTime*)parlist.FindObject("SelTime");
|
|---|
| 183 | MHEnergyTheta &seltheta = *(MHEnergyTheta*)parlist.FindObject("SelTheta");
|
|---|
| 184 |
|
|---|
| 185 | MHEnergyTime collareatime;
|
|---|
| 186 | MHEnergyTheta collareatheta;
|
|---|
| 187 | collareatime.Divide(&seltime, &alltime);
|
|---|
| 188 | collareatheta.Divide(&seltheta, &alltheta);
|
|---|
| 189 | */
|
|---|
| 190 |
|
|---|
| 191 | MHTimeDiffTime &effontime = *(MHTimeDiffTime*)parlist.FindObject("EffOnTime");
|
|---|
| 192 | MHTimeDiffTheta &effontheta = *(MHTimeDiffTheta*)parlist.FindObject("EffOnTheta");
|
|---|
| 193 |
|
|---|
| 194 | /*
|
|---|
| 195 | effontime.DrawClone();
|
|---|
| 196 | effontheta.DrawClone();
|
|---|
| 197 | */
|
|---|
| 198 |
|
|---|
| 199 | MHEffOnTimeTime ontime;
|
|---|
| 200 | MHEffOnTimeTheta ontheta;
|
|---|
| 201 | ontime.SetupFill(&parlist);
|
|---|
| 202 | ontheta.SetupFill(&parlist);
|
|---|
| 203 |
|
|---|
| 204 | ontime.Calc(effontime.GetHist());
|
|---|
| 205 | ontheta.Calc(effontheta.GetHist());
|
|---|
| 206 |
|
|---|
| 207 | /*
|
|---|
| 208 | ontime.DrawClone();
|
|---|
| 209 | ontheta.DrawClone();
|
|---|
| 210 | */
|
|---|
| 211 |
|
|---|
| 212 | MHAlphaEnergyTime &fluxsp = *(MHAlphaEnergyTime*)parlist.FindObject("FluxSrcTime");
|
|---|
| 213 | MHAlphaEnergyTime &fluxasp = *(MHAlphaEnergyTime*)parlist.FindObject("FluxASrcTime");
|
|---|
| 214 | MHAlphaEnergyTheta &fluxsptheta = *(MHAlphaEnergyTheta*)parlist.FindObject("FluxSrcTheta");
|
|---|
| 215 | MHAlphaEnergyTheta &fluxasptheta = *(MHAlphaEnergyTheta*)parlist.FindObject("FluxASrcTheta");
|
|---|
| 216 |
|
|---|
| 217 | /*
|
|---|
| 218 | fluxsp.DrawClone();
|
|---|
| 219 | fluxasp.DrawClone();
|
|---|
| 220 | fluxsptheta.DrawClone();
|
|---|
| 221 | fluxasptheta.DrawClone();
|
|---|
| 222 | */
|
|---|
| 223 |
|
|---|
| 224 | MHAlphaEnergyTime resulttime;
|
|---|
| 225 | MHAlphaEnergyTheta resulttheta;
|
|---|
| 226 | resulttime.Substract(&fluxsp, &fluxasp);
|
|---|
| 227 | resulttheta.Substract(&fluxsptheta, &fluxasptheta);
|
|---|
| 228 |
|
|---|
| 229 | /*
|
|---|
| 230 | resulttime.DrawClone();
|
|---|
| 231 | resulttheta.DrawClone();
|
|---|
| 232 | */
|
|---|
| 233 |
|
|---|
| 234 | TH2D &projecttime = *resulttime.GetAlphaProjection(-10, 10);
|
|---|
| 235 | TH2D &projecttheta = *resulttheta.GetAlphaProjection(-10, 10);
|
|---|
| 236 |
|
|---|
| 237 | projecttime.SetTitle("Number of Gammas vs. EnergyEst and Time (Alpha integrated between -10, 10deg)");
|
|---|
| 238 | projecttheta.SetTitle("Number of Gammas vs. EnergyEst and Theta (Alpha integrated between -10, 10deg)");
|
|---|
| 239 |
|
|---|
| 240 | /*
|
|---|
| 241 | TCanvas *c = new TCanvas("Unfold", "To be unfolded", 350, 500);
|
|---|
| 242 | c->Divide(1, 2);
|
|---|
| 243 | c->cd(1);
|
|---|
| 244 | projecttime.DrawCopy();
|
|---|
| 245 | c->cd(2);
|
|---|
| 246 | projecttheta.DrawCopy();
|
|---|
| 247 | */
|
|---|
| 248 |
|
|---|
| 249 | return;
|
|---|
| 250 |
|
|---|
| 251 | for (int i=1; i<=binstime.GetNumBins(); i++)
|
|---|
| 252 | {
|
|---|
| 253 | if (ontime.GetHist()->GetBinContent(i)==0)
|
|---|
| 254 | continue;
|
|---|
| 255 |
|
|---|
| 256 | TH1D &hist = *projecttime.ProjectionX("Number of Gammas vs. EnergyEst for a fixed time", i, i);
|
|---|
| 257 |
|
|---|
| 258 | /* UNFOLDING */
|
|---|
| 259 |
|
|---|
| 260 | //hist->Divide(collareatime);
|
|---|
| 261 | hist.Scale(1./ontime.GetHist()->GetBinContent(i));
|
|---|
| 262 |
|
|---|
| 263 | for (int j=1; j<=binse.GetNumBins(); j++)
|
|---|
| 264 | hist.SetBinContent(j, hist.GetBinContent(j)/hist.GetBinWidth(j));
|
|---|
| 265 |
|
|---|
| 266 | hist.SetName("Flux");
|
|---|
| 267 | hist.SetTitle("Flux[Gammas/s/m^2/GeV] vs. EnergyTrue for a fixed Time");
|
|---|
| 268 |
|
|---|
| 269 | char n[100];
|
|---|
| 270 | sprintf(n, "Canv%d", j);
|
|---|
| 271 | c= new TCanvas(n, "Title");
|
|---|
| 272 | hist.DrawCopy();
|
|---|
| 273 | }
|
|---|
| 274 |
|
|---|
| 275 | delete &projecttime;
|
|---|
| 276 | delete &projecttheta;
|
|---|
| 277 |
|
|---|
| 278 | return;
|
|---|
| 279 |
|
|---|
| 280 | // ------------------------------------------
|
|---|
| 281 |
|
|---|
| 282 | MHMcCollectionArea carea;
|
|---|
| 283 | TH1D *collareatime = carea.GetHist(); // FIXME!
|
|---|
| 284 | TH1D *collareatheta = carea.GetHist(); // FIXME!
|
|---|
| 285 | }
|
|---|
