Index: /trunk/MagicSoft/Mars/macros/flux.C =================================================================== --- /trunk/MagicSoft/Mars/macros/flux.C (revision 1293) +++ /trunk/MagicSoft/Mars/macros/flux.C (revision 1294) @@ -16,6 +16,5 @@ ! ! -! Author(s): Thomas Bretz 1/2002 -! Author(s): Wolfgang Wittek 1/2002 +! Author(s): Wolfgang Wittek 4/2002 ! ! Copyright: MAGIC Software Development, 2000-2002 @@ -23,16 +22,48 @@ ! \* ======================================================================== */ - +////////////////////////////////////////////////////////////////////////////// +// // +// 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() { - // - // first we have to create our empty lists + //-------------------------------------------------------------- + // empty lists of parameter containers and tasks // MParList parlist; MTaskList tasklist; - parlist.AddToList(&tasklist); - - // + //-------------------------------------------------------------- + // Geometry information of the MAGIC camera + // + MGeomCamMagic geomcam; + + //-------------------------------------------------------------- // Define the two source positions // @@ -41,24 +72,46 @@ 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); - // - // 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(7, -2.5, 32.5); - - MBinning binstime("BinningTime"); - binstime.SetEdges(5, 0, 50); - - MBinning binsdifftime("BinningTimeDiff"); - binsdifftime.SetEdges(50, 0, 0.1); + parlist.AddToList(&binswidth); + parlist.AddToList(&binslength); + parlist.AddToList(&binsdist); parlist.AddToList(&binsalpha); @@ -68,14 +121,31 @@ parlist.AddToList(&binsdifftime); - // - // Setup our tasks - // - MReadMarsFile reader("Events", "~/data/Gamma*.root"); + //-------------------------------------------------------------- + // 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"); - MGeomCamMagic geomcam; - parlist.AddToList(&geomcam); - + //..................................... + // generation of event time MMcTimeGenerate rand; + + //..................................... + // collection area + MMcCollectionAreaCalc acalc; + MMcPedestalCopy pcopy; MMcPedestalNSBAdd pnsb; @@ -83,51 +153,92 @@ 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"); - MMcCollectionAreaCalc acalc; + //..................................... + // filters + MFTriggerLvl1 lvl1; const Float_t alpha0 = 15; // [deg] - MFAlpha fsrc ("HillasSrc", '>', alpha0); MFAlpha fasrc("HillasAntiSrc", '>', alpha0); - MFillH fillsp ("FluxSrcTime [MHAlphaEnergyTime]", "HillasSrc"); - MFillH fillasp ("FluxASrcTime [MHAlphaEnergyTime]", "HillasAntiSrc"); + //..................................... + // 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"); - - 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(&fasrc); fillasptheta.SetFilter(&fsrc); - MFillH fillontime ("EffOnTime [MHTimeDiffTime]", "MMcEvt"); - MFillH fillontheta("EffOnTheta [MHTimeDiffTheta]", "MMcEvt"); - - // - // Setup Task list + //..................................... + // 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(&fillontime); - tasklist.AddToList(&fillontheta); + tasklist.AddToList(&acalc); + tasklist.AddToList(&pcopy); tasklist.AddToList(&pnsb); @@ -139,22 +250,36 @@ 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(&fillsp); + + tasklist.AddToList(&fillsptime); + tasklist.AddToList(&fillasptime); + + tasklist.AddToList(&fillsptheta); 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 + //---------------------------------------------------------------------- + // Event loop // MEvtLoop magic; @@ -162,102 +287,177 @@ // - // Start to loop over all events + // loop over all events // if (!magic.Eventloop()) return; - - tasklist.PrintStatistics(); + //---------------------------------------------------------------------- + + + 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"); - 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"); - + //---------------------------------------------------------------------- + // 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"); /* - effontime.DrawClone(); - effontheta.DrawClone(); + dtimetime->DrawClone(); + dtimetheta->DrawClone(); */ - MHEffOnTimeTime ontime; - MHEffOnTimeTheta ontheta; - ontime.SetupFill(&parlist); - ontheta.SetupFill(&parlist); - - ontime.Calc(effontime.GetHist()); - ontheta.Calc(effontheta.GetHist()); + + //.................................... + // 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; + */ /* - ontime.DrawClone(); - ontheta.DrawClone(); - */ - - MHAlphaEnergyTime &fluxsp = *(MHAlphaEnergyTime*)parlist.FindObject("FluxSrcTime"); - MHAlphaEnergyTime &fluxasp = *(MHAlphaEnergyTime*)parlist.FindObject("FluxASrcTime"); - MHAlphaEnergyTheta &fluxsptheta = *(MHAlphaEnergyTheta*)parlist.FindObject("FluxSrcTheta"); - MHAlphaEnergyTheta &fluxasptheta = *(MHAlphaEnergyTheta*)parlist.FindObject("FluxASrcTheta"); + 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); + /* - fluxsp.DrawClone(); - fluxasp.DrawClone(); - fluxsptheta.DrawClone(); - fluxasptheta.DrawClone(); - */ - - MHAlphaEnergyTime 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)"); - - /* - TCanvas *c = new TCanvas("Unfold", "To be unfolded", 350, 500); - c->Divide(1, 2); - c->cd(1); - projecttime.DrawCopy(); - c->cd(2); - projecttheta.DrawCopy(); - */ + 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 (ontime.GetHist()->GetBinContent(i)==0) + if (effontime.GetHist()->GetBinContent(i)==0) continue; - TH1D &hist = *projecttime.ProjectionX("Number of Gammas vs. EnergyEst for a fixed time", i, i); + TH1D &hist = *evttime.ProjectionX("Number of Gammas vs. EnergyEst for a fixed time", i, i); /* UNFOLDING */ //hist->Divide(collareatime); - hist.Scale(1./ontime.GetHist()->GetBinContent(i)); + hist.Scale(1./effontime.GetHist()->GetBinContent(i)); for (int j=1; j<=binse.GetNumBins(); j++) @@ -270,9 +470,11 @@ sprintf(n, "Canv%d", j); c= new TCanvas(n, "Title"); + /* hist.DrawCopy(); + */ } - delete &projecttime; - delete &projecttheta; + delete &evttime; + delete &evttheta; return; @@ -284,2 +486,5 @@ TH1D *collareatheta = carea.GetHist(); // FIXME! } +//=========================================================================== + + Index: /trunk/MagicSoft/Mars/mhist/MHEffOnTimeTime.h =================================================================== --- /trunk/MagicSoft/Mars/mhist/MHEffOnTimeTime.h (revision 1293) +++ /trunk/MagicSoft/Mars/mhist/MHEffOnTimeTime.h (revision 1294) @@ -16,5 +16,8 @@ { private: - TH1D fHist; + TH1D fHEffOn; + TH1D fHChi2; + TH1D fHLambda; + TH1D fHN0del; public: @@ -24,6 +27,6 @@ virtual Bool_t Fill(const MParContainer *par); - const TH1D *GetHist() { return &fHist; } - const TH1D *GetHist() const { return &fHist; } + const TH1D *GetHist() { return &fHEffOn; } + const TH1D *GetHist() const { return &fHEffOn; } void Calc(TH2D *hist); @@ -32,5 +35,5 @@ TObject *DrawClone(Option_t *option="") const; - ClassDef(MHEffOnTimeTime, 1) //Histogram to store a 3-Dim histogram in alpha, Energy and time + ClassDef(MHEffOnTimeTime, 1) //1D-plot of Delta t vs. time };