/* ======================================================================== *\ ! ! * ! * 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): Wolfgang Wittek 4/2002 ! ! Copyright: MAGIC Software Development, 2000-2002 ! ! \* ======================================================================== */ ////////////////////////////////////////////////////////////////////////////// // // // This macro calculates photon fluxes as a function of // // - energy and time // // - energy and Theta // // // // It is assumed that the data have been taken in the wobble mode. // // This means that there is only one set of data, from which both // // 'on' and 'off' data are constructed. // // // // Necessary input from MC : // // - migration matrix (E_est, E_true) as a functioin of Theta // // - effective collection areas as a function of energy and Theta // // // // // // The input from MC has to refer to the wobble mode too. // // // // The macro includes : // // - the calculation of Hillas parameters // // - the calculation of the effective on time // // - the unfolding of the distributions in the estimated energy E_est // // // // For testing purposes (as long as no real data, no energy estimator, // // no migration matrices and no effective collection areas are available) // // - event times are generated (according to a Poissonian with dead time)// // - a dummy energy estimator is provided // // - migration matrices are generated // // - dummy effective collection areas are used // // // ////////////////////////////////////////////////////////////////////////////// void flux() { //-------------------------------------------------------------- // empty lists of parameter containers and tasks // MParList parlist; MTaskList tasklist; //-------------------------------------------------------------- // Geometry information of the MAGIC camera // MGeomCamMagic geomcam; //-------------------------------------------------------------- // Define the two source positions // MSrcPosCam source("Source"); MSrcPosCam antisrc("AntiSource"); source.SetXY(0, 0); antisrc.SetXY(+240, 0); //-------------------------------------------------------------- // Setup the binning for the histograms // //............................................................... // These are NOT the binnings for the flux determination MBinning binswidth("BinningWidth"); binswidth.SetEdges(100, 0, 0.5); // 100 bins from 0 to 0.5 deg MBinning binslength("BinningLength"); binslength.SetEdges(100, 0, 1); // 100 bins from 0 to 1 deg MBinning binscamera("BinningCamera"); binscamera.SetEdges(90, -2.25, 2.25); // 90 bins from -2.25 to 2.25 deg MBinning binsdist("BinningDist"); binsdist.SetEdges(90, 0, 2.25); // 90 bins from 0 to 2.25 deg MBinning binsalpha("BinningAlpha"); binsalpha.SetEdges(100, -100, 100); MBinning binsasym("BinningAsym"); binsasym.SetEdges(100, -0.5, 0.5); // 100 bins from -0.5 to 0.5 deg MBinning binsasymn("BinningAsymn"); binsasymn.SetEdges(100, -0.005, 0.005); // 100 bins from -0.005 to 0.005 deg MBinning binsm3long("BinningM3Long"); binsm3long.SetEdges(100, -0.5, 0.5); // 100 bins from -0.5 to 0.5 deg MBinning binsm3trans("BinningM3Trans"); binsm3trans.SetEdges(100, -0.5, 0.5); // 100 bins from -0.5 to 0.5 deg //............................................................... // These ARE the binnings for the flux determination MBinning binsalphaflux("BinningAlphaFlux"); binsalphaflux.SetEdges(100, 0, 100); MBinning binse("BinningE"); binse.SetEdgesLog(10, 10, 1e3); MBinning binstheta("BinningTheta"); binstheta.SetEdges(7, -2.5, 32.5); MBinning binstime("BinningTime"); binstime.SetEdges(5, 0, 100); MBinning binsdifftime("BinningTimeDiff"); binsdifftime.SetEdges(50, 0, 0.1); //-------------------------------------------------------------- // Fill list of parameter containers // parlist.AddToList(&tasklist); parlist.AddToList(&geomcam); parlist.AddToList(&source); parlist.AddToList(&antisrc); parlist.AddToList(&binswidth); parlist.AddToList(&binslength); parlist.AddToList(&binscamera); parlist.AddToList(&binsdist); parlist.AddToList(&binsalpha); parlist.AddToList(&binsasym); parlist.AddToList(&binsasymn); parlist.AddToList(&binsm3long); parlist.AddToList(&binsm3trans); parlist.AddToList(&binsalphaflux); parlist.AddToList(&binse); parlist.AddToList(&binstheta); parlist.AddToList(&binstime); parlist.AddToList(&binsdifftime); //-------------------------------------------------------------- // Setup the tasks // //..................................... // input file // // next statement commented out // MReadMarsFile reader("Events", "~/data/Gamma*.root"); MReadMarsFile reader("Events"); //reader.AddFile("/hd31/rkb/Camera/mpi/Gamma*.root"); //reader.AddFile("/hd31/rkb/Camera/mpi/prot_N_*.root"); //reader.AddFile("/hd31/rkb/Camera/mpi/prot_NS_*.root"); //reader.AddFile("protondata/gamma_10_n*.root"); reader.AddFile("protondata/gamma_10_cn*.root"); reader.EnableBranch("MMcEvt.fTheta"); //..................................... // filters MFTriggerLvl1 lvl1; //..................................... // generation of event time MMcTimeGenerate rand; //..................................... MMcPedestalCopy pcopy; MMcPedestalNSBAdd pnsb; MCerPhotCalc ncalc; MImgCleanStd clean; MBlindPixelCalc blind; //..................................... // source independent image parameters MHillasCalc hcalc; MHillasExt hext; parlist.AddToList(&hext); //..................................... // source dependent image parameters MHillasSrcCalc hsrc1("Source", "HillasSrc"); MHillasSrcCalc hsrc2("AntiSource", "HillasAntiSrc"); //..................................... // energy estimation MEnergyEstimate estim; //..................................... // migration matrix (MC) MFillH eestetrue("MHMcEnergyMigration", "HillasSrc"); eestetrue.SetTitle("Task to determine the migration matrix for the energy"); //..................................... // plots of source independent parameters (length, width, delta, size, center) MFillH hfill1h("MHHillas", "MHillas"); hfill1h.SetTitle("Task to plot the source independent image parameters"); hfill1h.SetFilter(&lvl1); MFillH hfill1x("MHHillasExt", "MHillas"); hfill1x.SetTitle("Task to plot the extended image parameters"); hfill1x.SetFilter(&lvl1); //..................................... // plots of source dependent parameters (alpha, dist) MFillH hfill2s("HSource [MHHillasSrc]", "HillasSrc"); hfill2s.SetTitle("Task to plot the source dependent image parameters (Source)"); hfill2s.SetFilter(&lvl1); MFillH hfill2a("HAntiSource [MHHillasSrc]", "HillasAntiSrc"); hfill2a.SetTitle("Task to plot the source dependent image parameters (AntiSource)"); hfill2a.SetFilter(&lvl1); //..................................... // star map MFillH hfill1m("MHStarMap", "MHillas"); hfill1m.SetTitle("Task to plot the Star map"); //..................................... const Float_t alpha0 = 15; // [deg] MFAlpha fsrc ("HillasSrc", '>', alpha0); fsrc.SetTitle("Task to check alphaSrc > alpha0"); MFAlpha fasrc("HillasAntiSrc", '>', alpha0); fasrc.SetTitle("Task to check alphaAntiSrc > alpha0"); //..................................... // 1-D profile plots (, time) // (, Theta) MFillH fillthetabartime ("ThetabarTime [MHThetabarTime]", "MMcEvt"); fillthetabartime.SetTitle("Task to plot vs time"); // fillthetabartime.SetFilter(&lvl1); MFillH fillthetabartheta("ThetabarTheta [MHThetabarTheta]", "MMcEvt"); fillthetabartheta.SetTitle("Task to plot vs theta"); // fillthetabartheta.SetFilter(&lvl1); //..................................... // 2-D plots (Delta(time), time) // (Delta(time), Theta) MFillH filldtimetime ("EffOnTime [MHTimeDiffTime]", "MMcEvt"); filldtimetime.SetTitle("Task to plot Delta(time) vs time"); // filldtime.SetFilter(&lvl1); MFillH filldtimetheta("EffOnTheta [MHTimeDiffTheta]", "MMcEvt"); filldtimetheta.SetTitle("Task to plot Delta(time) vs theta"); // fillontheta.SetFilter(&lvl1); //..................................... // 3-D plots (alpha, E_est, time) MFillH fillsptime ("SrcTime [MHAlphaEnergyTime]", "HillasSrc"); fillsptime.SetTitle("Task for 3D plots (alpha, E_est, time) (Source)"); fillsptime.SetFilter(&fasrc); MFillH fillasptime ("ASrcTime [MHAlphaEnergyTime]", "HillasAntiSrc"); fillasptime.SetTitle("Task for 3D plots (alpha, E_est, time) (AntiSource)"); fillasptime.SetFilter(&fsrc); //..................................... // 3-D plots (alpha, E_est, Theta) MFillH fillsptheta ("SrcTheta [MHAlphaEnergyTheta]", "HillasSrc"); fillsptheta.SetTitle("Task for 3D plots (alpha, E_est, Theta) (Source)"); fillsptheta.SetFilter(&fasrc); MFillH fillasptheta("ASrcTheta [MHAlphaEnergyTheta]", "HillasAntiSrc"); fillasptheta.SetTitle("Task for 3D plots (alpha, E_est, Theta) (AntiSource)"); fillasptheta.SetFilter(&fsrc); //..................................... // MFillH fethetaall("AllTheta [MHEnergyTheta]", "MMcEvt"); // MFillH fethetasel("SelTheta [MHEnergyTheta]", "MMcEvt"); // fethetasel.SetFilter(&lvl1); // MFillH fetimeall ("AllTime [MHEnergyTime]", "MMcEvt"); // MFillH fetimesel ("SelTime [MHEnergyTime]", "MMcEvt"); // fetimesel.SetFilter(&lvl1); //--------------------------------------------------------------------- // Setup the task list // tasklist.AddToList(&reader); tasklist.AddToList(&lvl1); tasklist.AddToList(&rand); 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(&eestetrue); tasklist.AddToList(&hfill1h); tasklist.AddToList(&hfill1x); tasklist.AddToList(&hfill1m); tasklist.AddToList(&hfill2s); tasklist.AddToList(&hfill2a); tasklist.AddToList(&fillthetabartime); tasklist.AddToList(&fillthetabartheta); tasklist.AddToList(&filldtimetime); tasklist.AddToList(&filldtimetheta); // tasklist.AddToList(&fethetaall); // tasklist.AddToList(&fethetasel); // tasklist.AddToList(&fetimeall); // tasklist.AddToList(&fetimesel); tasklist.AddToList(&fsrc); tasklist.AddToList(&fasrc); tasklist.AddToList(&fillsptime); tasklist.AddToList(&fillasptime); tasklist.AddToList(&fillsptheta); tasklist.AddToList(&fillasptheta); //---------------------------------------------------------------------- // Event loop // MEvtLoop magic; magic.SetParList(&parlist); // // loop over all events // if (!magic.Eventloop()) return; //---------------------------------------------------------------------- tasklist.PrintStatistics(10); MHHillas *mhhillas = (MHHillas*)parlist.FindObject("MHHillas"); mhhillas->SetTitle("Source indep. image parameters"); mhhillas->DrawClone(); MHHillasExt *mhhillasext = (MHHillasExt*)parlist.FindObject("MHHillasExt"); mhhillasext->SetTitle("Extended image parameters"); mhhillasext->DrawClone(); MHHillasSrc *hsource = (MHHillasSrc*)parlist.FindObject("HSource"); hsource->SetTitle("Source dependent image parameters (Source)"); hsource->DrawClone(); MHHillasSrc *hantisource = (MHHillasSrc*)parlist.FindObject("HAntiSource"); hantisource->SetTitle("Source dependent image parameters (AntiSource)"); hantisource->DrawClone(); /* parlist.FindObject("MHStarMap")->DrawClone(); */ //---------------------------------------------------------------------- // average Theta versus time // and versus Theta // This is needed for selecting later the right collection area // (at the right Theta) for each bin in time or Theta // MHThetabarTime *bartime = (MHThetabarTime*)parlist.FindObject("ThetabarTime"); MHThetabarTheta *bartheta = (MHThetabarTheta*)parlist.FindObject("ThetabarTheta"); /* bartime->DrawClone(); bartheta->DrawClone(); */ //---------------------------------------------------------------------- // Effective on time //.................................... // get plots of time differences for different intervals in time // and plots of time differences for different intervals in Theta MHTimeDiffTime *dtimetime = (MHTimeDiffTime*)parlist.FindObject("EffOnTime"); MHTimeDiffTheta *dtimetheta = (MHTimeDiffTheta*)parlist.FindObject("EffOnTheta"); /* dtimetime->DrawClone(); dtimetheta->DrawClone(); */ //.................................... // calculate effective on time for different intervals in Time // and for different intervals in Theta MHEffOnTime effontime ("Time", "[s]"); MHEffOnTime effontheta("Theta", "[\\circ]"); effontime.SetupFill(&parlist); effontheta.SetupFill(&parlist); // Draw == 0 don't draw the individual distributions of time differences // != 0 draw them Bool_t draw=kFALSE; effontime.Calc (dtimetime->GetHist(), draw); effontheta.Calc(dtimetheta->GetHist(),draw); // plot effective on time versus Time // and versus Theta /* effontime.DrawClone(); effontheta.DrawClone(); */ //====================================================================== // // A : fully differential plots (Alpha, E-est, Var) // //---------------------------------------------------------------------- // 3D-plots of alpha, E_est and time // and of alpha, E_est and Theta // with alpha calculated relative to the source position // and relative to the anti-source position MHAlphaEnergyTime *evtsptime = (MHAlphaEnergyTime*)parlist.FindObject("SrcTime"); MHAlphaEnergyTime *evtasptime = (MHAlphaEnergyTime*)parlist.FindObject("ASrcTime"); MHAlphaEnergyTheta *evtsptheta = (MHAlphaEnergyTheta*)parlist.FindObject("SrcTheta"); MHAlphaEnergyTheta *evtasptheta = (MHAlphaEnergyTheta*)parlist.FindObject("ASrcTheta"); /* evtsptime->SetTitle("Source : 3D-plot Al-En-time, alpha\_{asp} .gt. alpha\_0"); evtsptime->DrawClone(); evtasptime->SetTitle("AntiSource : 3D-plot Al-En-time, alpha\_{sp} .gt. alpha\_0"); evtasptime->DrawClone(); evtsptheta->SetTitle("Source : 3D-plot Al-En-Theta, alpha\_{asp} .gt. alpha\_0"); evtsptheta->DrawClone(); evtasptheta->SetTitle("AntiSource : 3D-plot Al-En-Theta, alpha\_{sp} .gt. alpha\_0"); evtasptheta->DrawClone(); */ //---------------------------------------------------------------------- // Difference between source and antisource (= 'gamma' sample) // for 3D-plots of alpha, E_est and time // and of alpha, E_est and Theta // 3rd and 4th argument are name and title of the resulting histogram MHGamma gamma; TH3D *hsubtime = gamma.Subtract( evtsptime->GetHist(), evtasptime->GetHist(), "Al-En-time", "3D-plot of Alpha,E-est,time", draw); TH3D *hsubtheta = gamma.Subtract( evtsptheta->GetHist(), evtasptheta->GetHist(), "Al-En-time", "3D-plot of Alpha,E-est,Theta", draw); //---------------------------------------------------------------------- // get number of 'gammas' in the alpha range (lo, up) as a function of // E_est and time // or E_est and Theta Axis_t lo = 0; // [deg] Axis_t up = 10; // [deg] const TH2D &evttime = *(gamma.GetAlphaProjection(hsubtime, lo, up, draw)); const TH2D &evttheta = *(gamma.GetAlphaProjection(hsubtheta, lo, up, draw)); //---------------------------------------------------------------------- // plot migration matrix: E-true E-est Theta /* MHMcEnergyMigration *migrationmatrix = (MHMcEnergyMigration*)parlist.FindObject("MHMcEnergyMigration"); migrationmatrix->SetTitle("Migration matrix");; migrationmatrix->DrawClone(); */ //---------------------------------------------------------------------- // determine gamma flux from distributions of E-est : // - do unfolding to get the distributions in E-unfold // - divide by bin width in energy // effective ontime and // effective collection area // to get the absolute gamma flux //............................................. // get flux spectrum for different bins in Time // // use dummy histogram *aeff which eventually will contain // the effective collection area as a function of Etrue and Theta TH2D aeff; // arguments of MHFlux constructor are : // - 2D histogram containing the no.of gammas as a function of // E-est and Var // - name of variable Var ("Time" or "Theta") // - units of variable Var ( "[s]" or "[\\circ]") // Draw == kTRUE draw the no.of photons vs. E-est // for the individual bins of the variable Var draw=kTRUE; MHFlux fluxtime(evttime, draw, "Time", "[s]"); fluxtime.Unfold(draw); fluxtime.CalcFlux(effontime.GetHist(), bartime.GetHist(), &aeff, draw); fluxtime.DrawClone(); //.............................................. // get flux spectrum for different bins in Theta MHFlux fluxtheta(evttheta, draw, "Theta", "[\\circ]"); fluxtheta.Unfold(draw); fluxtheta.CalcFlux(effontheta.GetHist(), bartheta.GetHist(), &aeff, draw); fluxtheta.DrawClone(); return; //---------------------------------------------------------------------- //====================================================================== // // B : plots integrated over E-est // //---------------------------------------------------------------------- // integrate 3D-plot of (alpha, E_est and Time) over E-est // to get 2D-plot of (alpha, Time) Draw = kFALSE; TH2D *evttimesp = evtsptime.IntegrateEest ("AlphaSP vs. Time", Draw); TH2D *evttimeasp = evtasptime.IntegrateEest("AlphaASP vs. Time",Draw); //====================================================================== // // C : plots integrated over Time // //---------------------------------------------------------------------- // integrate 3D-plot of (alpha, E_est and Time) over the Time // to get 2D-plot of (alpha, E_est) Draw = kFALSE; TH2D *evteestsp = evtsptime.IntegrateTime ("AlphaSP vs. E-est", Draw); TH2D *evteestasp = evtasptime.IntegrateTime("AlphaASP vs. E-est",Draw); //====================================================================== // // D : plots integrated over E-est and Time // //---------------------------------------------------------------------- // integrate 3D-plot of (alpha, E_est and Time) over E-est and Time // to get 1D-plot of (alpha) Draw = kFALSE; TH1D *evtalphasp = evtsptime.IntegrateEestTime ("AlphaSP", Draw); TH1D *evtalphaasp = evtasptime.IntegrateEestTime("AlphaASP",Draw); //---------------------------------------------------------------------- // list of tasks for processing cout << " " << endl; cout << "List of tasks for processing :" << endl; cout << "==============================" << endl; cout << " " << endl; // was folgt geht nicht : // Error: Function Next() is not defined in current scope FILE:macros/flux.C LINE:636 //while ( (task=(MTask*)Next()) ) //{ // cout << tasklist.GetName() << " : " << tasklist.GetTitle() << endl; //} return; //====================================================================== } //===========================================================================