/* ======================================================================== *\ ! ! * ! * 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 // // - a kind of effective collection areas are calculated // // // ////////////////////////////////////////////////////////////////////////////// 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 binning for the histograms // MBinning binswidth("BinningWidth"); binswidth.SetEdges(100, 0, 1); // 100 bins from 0 to 1 deg MBinning binslength("BinningLength"); binslength.SetEdges(100, 0, 1); // 100 bins from 0 to 1 deg MBinning binsdist("BinningDist"); binsdist.SetEdges(100, 0, 2); // 100 bins from 0 to 2 deg MBinning binsalpha("BinningAlpha"); binsalpha.SetEdges(45, -90, 90); 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(&binsdist); parlist.AddToList(&binsalpha); 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", "/hd31/rkb/Camera/mpi/Gamma*.root"); //reader.AddFile("/hd31/rkb/Camera/mpi/prot_N_*.root"); // next 2 statements added MReadMarsFile reader("Events"); reader.AddFile("protondata/gamma_10_cn*.root"); reader.EnableBranch("MMcEvt.fTheta"); //..................................... // generation of event time MMcTimeGenerate rand; //..................................... // collection area MMcCollectionAreaCalc acalc; MMcPedestalCopy pcopy; MMcPedestalNSBAdd pnsb; MCerPhotCalc ncalc; MImgCleanStd clean; MBlindPixelCalc blind; //..................................... // source independent image parameters MHillasCalc hcalc; //..................................... // source dependent image parameters MHillasSrcCalc hsrc1("Source", "HillasSrc"); MHillasSrcCalc hsrc2("AntiSource", "HillasAntiSrc"); //..................................... // energy estimation MEnergyEstimate estim; //..................................... // migration matrix (MC) MFillH eestetrue("MHMcEnergyMigration", "HillasSrc"); //..................................... // plots of width and length MFillH hfill1h("MHHillas", "MHillas"); //..................................... // star map MFillH hfill1m("MHStarMap", "MHillas"); //..................................... // plots of alpha and dist MFillH hfill2s("HSource [MHHillasSrc]", "HillasSrc"); MFillH hfill2a("HAntiSource [MHHillasSrc]", "HillasAntiSrc"); //..................................... // filters MFTriggerLvl1 lvl1; const Float_t alpha0 = 15; // [deg] MFAlpha fsrc ("HillasSrc", '>', alpha0); MFAlpha fasrc("HillasAntiSrc", '>', alpha0); //..................................... // 1-D profile plots (, time) // (, Theta) MFillH fillthetabartime ("ThetabarTime [MHThetabarTime]", "MMcEvt"); // fillthetabartime.SetFilter(&lvl1); MFillH fillthetabartheta("ThetabarTheta [MHThetabarTheta]", "MMcEvt"); // fillthetabartheta.SetFilter(&lvl1); //..................................... // 2-D plots (Delta(time), time) // (Delta(time), Theta) MFillH filldtimetime ("EffOnTime [MHTimeDiffTime]", "MMcEvt"); // filldtime.SetFilter(&lvl1); MFillH filldtimetheta("EffOnTheta [MHTimeDiffTheta]", "MMcEvt"); // fillontheta.SetFilter(&lvl1); //..................................... // 3-D plots (alpha, E_est, time) MFillH fillsptime ("FluxSrcTime [MHAlphaEnergyTime]", "HillasSrc"); fillsptime.SetFilter(&fasrc); MFillH fillasptime ("FluxASrcTime [MHAlphaEnergyTime]", "HillasAntiSrc"); fillasptime.SetFilter(&fsrc); //..................................... // 3-D plots (alpha, E_est, Theta) MFillH fillsptheta ("FluxSrcTheta [MHAlphaEnergyTheta]", "HillasSrc"); fillsptheta.SetFilter(&fasrc); MFillH fillasptheta("FluxASrcTheta [MHAlphaEnergyTheta]", "HillasAntiSrc"); 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(&rand); 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(&eestetrue); tasklist.AddToList(&hfill1h); tasklist.AddToList(&hfill1m); tasklist.AddToList(&hfill2s); tasklist.AddToList(&hfill2a); tasklist.AddToList(&fillthetabartime); tasklist.AddToList(&fillthetabartheta); tasklist.AddToList(&filldtimetime); tasklist.AddToList(&filldtimetheta); // tasklist.AddToList(&lvl1); // 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); /* parlist.FindObject("HSource")->DrawClone();; parlist.FindObject("MHHillas")->DrawClone(); */ /* parlist.FindObject("MHStarMap")->DrawClone(); */ /* parlist.FindObject("HAntiSource")->DrawClone(); */ /* parlist.FindObject("MHMcCollectionArea")->DrawClone(); */ //---------------------------------------------------------------------- // average Theta versus time // and versus Theta // MHThetabarTime *bartime = (MHThetabarTime*)parlist.FindObject("ThetabarTime"); MHThetabarTheta *bartheta = (MHThetabarTheta*)parlist.FindObject("ThetabarTheta"); /* bartime->DrawClone(); bartheta->DrawClone(); */ //---------------------------------------------------------------------- // this is something like the collection area // as a function of time // and as a function of Theta /* MHEnergyTime &alltime = *(MHEnergyTime*)parlist.FindObject("AllTime"); MHEnergyTime &seltime = *(MHEnergyTime*)parlist.FindObject("SelTime"); MHEnergyTime collareatime; collareatime.Divide(&seltime, &alltime); MHEnergyTheta &alltheta = *(MHEnergyTheta*)parlist.FindObject("AllTheta"); MHEnergyTheta &seltheta = *(MHEnergyTheta*)parlist.FindObject("SelTheta"); MHEnergyTheta collareatheta; collareatheta.Divide(&seltheta, &alltheta); */ //---------------------------------------------------------------------- // 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 MHEffOnTimeTime effontime; MHEffOnTimeTheta effontheta; effontime.SetupFill(&parlist); effontheta.SetupFill(&parlist); effontime.Calc(dtimetime->GetHist()); effontheta.Calc(dtimetheta->GetHist()); // plot effective on time versus time // and versus Theta effontime.DrawClone(); effontheta.DrawClone(); //---------------------------------------------------------------------- // 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("FluxSrcTime"); MHAlphaEnergyTime *evtasptime = (MHAlphaEnergyTime*)parlist.FindObject("FluxASrcTime"); MHAlphaEnergyTheta *evtsptheta = (MHAlphaEnergyTheta*)parlist.FindObject("FluxSrcTheta"); MHAlphaEnergyTheta *evtasptheta = (MHAlphaEnergyTheta*)parlist.FindObject("FluxASrcTheta"); /* this test shows that the addresses 'evtsptime' and '&evtspobj' are identical MHAlphaEnergyTime &evtspobj = *(MHAlphaEnergyTime*)parlist.FindObject("FluxSrcTime"); cout << "evtsptime = " << evtsptime << endl; cout << "&evtspobj = " << &evtspobj << endl; */ /* evtsptime->SetTitle("3D-plot Al-En-time, alpha\_{asp} .gt. alpha\_0"); evtsptime->DrawClone(); evtasptime->SetTitle("3D-plot Al-En-time, alpha\_{sp} .gt. alpha\_0"); evtasptime->DrawClone(); evtsptheta->SetTitle("3D-plot Al-En-Theta, alpha\_{asp} .gt. alpha\_0"); evtsptheta->DrawClone(); evtasptheta->SetTitle("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 MHAlphaEnergyTime gammatime("Al-En-time","3D-plot Al-En-time"); MHAlphaEnergyTheta gammatheta("Al-En-Theta","3D-plot Al-En-Theta"); gammatime.SetupFill(&parlist); gammatime.Subtract(evtsptime, evtasptime); gammatheta.SetupFill(&parlist); gammatheta.Subtract(evtsptheta, evtasptheta); /* gammatime.SetTitle("3D-plot Al-En-time for \'gamma\' sample"); gammatime.DrawClone(); gammatheta.SetTitle("3D-plot Al-En-Theta for \'gamma\' sample"); gammatheta.DrawClone(); */ //---------------------------------------------------------------------- // 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 = -10; Axis_t up = 10; // TH2D &evttime = *gammatime.GetAlphaProjection (lo, up); // TH2D &evttheta = *gammatheta.GetAlphaProjection(lo, up); TH2D *evttime = gammatime.DrawAlphaProjection (lo, up); TH2D *evttheta = gammatheta.DrawAlphaProjection(lo, up); //---------------------------------------------------------------------- // plot migration matrix: E-true E-est Theta parlist.FindObject("MHMcEnergyMigration")->DrawClone(); return; //---------------------------------------------------------------------- for (int i=1; i<=binstime.GetNumBins(); i++) { if (effontime.GetHist()->GetBinContent(i)==0) continue; TH1D &hist = *evttime.ProjectionX("Number of Gammas vs. EnergyEst for a fixed time", i, i); /* UNFOLDING */ //hist->Divide(collareatime); hist.Scale(1./effontime.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 &evttime; delete &evttheta; return; // ------------------------------------------ MHMcCollectionArea carea; TH1D *collareatime = carea.GetHist(); // FIXME! TH1D *collareatheta = carea.GetHist(); // FIXME! } //===========================================================================