Index: /trunk/MagicSoft/Mars/macros/flux.C =================================================================== --- /trunk/MagicSoft/Mars/macros/flux.C (revision 1416) +++ /trunk/MagicSoft/Mars/macros/flux.C (revision 1417) @@ -49,5 +49,5 @@ // - a dummy energy estimator is provided // // - migration matrices are generated // -// - a kind of effective collection areas are calculated // +// - dummy effective collection areas are used // // // ////////////////////////////////////////////////////////////////////////////// @@ -74,17 +74,42 @@ //-------------------------------------------------------------- - // Setup binning for the histograms - // + // Setup the binning for the histograms + // + + //............................................................... + // These are NOT the binnings for the flux determination + MBinning binswidth("BinningWidth"); - binswidth.SetEdges(100, 0, 1); // 100 bins from 0 to 1 deg + 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(100, 0, 2); // 100 bins from 0 to 2 deg + binsdist.SetEdges(90, 0, 2.25); // 90 bins from 0 to 2.25 deg MBinning binsalpha("BinningAlpha"); - binsalpha.SetEdges(45, -90, 90); + 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"); @@ -113,7 +138,16 @@ 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); @@ -131,12 +165,17 @@ // 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("/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; //..................................... @@ -145,7 +184,4 @@ //..................................... - // collection area - MMcCollectionAreaCalc acalc; - MMcPedestalCopy pcopy; MMcPedestalNSBAdd pnsb; @@ -157,4 +193,6 @@ // source independent image parameters MHillasCalc hcalc; + MHillasExt hext; + parlist.AddToList(&hext); //..................................... @@ -170,26 +208,37 @@ // migration matrix (MC) MFillH eestetrue("MHMcEnergyMigration", "HillasSrc"); - - - //..................................... - // plots of width and length + 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"); - - //..................................... - // plots of alpha and dist - MFillH hfill2s("HSource [MHHillasSrc]", "HillasSrc"); - MFillH hfill2a("HAntiSource [MHHillasSrc]", "HillasAntiSrc"); - - //..................................... - // filters - MFTriggerLvl1 lvl1; - + 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"); //..................................... @@ -197,6 +246,9 @@ // (, 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); @@ -205,20 +257,29 @@ // (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 ("FluxSrcTime [MHAlphaEnergyTime]", "HillasSrc"); + MFillH fillsptime ("SrcTime [MHAlphaEnergyTime]", "HillasSrc"); + fillsptime.SetTitle("Task for 3D plots (alpha, E_est, time) (Source)"); fillsptime.SetFilter(&fasrc); - MFillH fillasptime ("FluxASrcTime [MHAlphaEnergyTime]", "HillasAntiSrc"); + + 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 ("FluxSrcTheta [MHAlphaEnergyTheta]", "HillasSrc"); + MFillH fillsptheta ("SrcTheta [MHAlphaEnergyTheta]", "HillasSrc"); + fillsptheta.SetTitle("Task for 3D plots (alpha, E_est, Theta) (Source)"); fillsptheta.SetFilter(&fasrc); - MFillH fillasptheta("FluxASrcTheta [MHAlphaEnergyTheta]", "HillasAntiSrc"); + + MFillH fillasptheta("ASrcTheta [MHAlphaEnergyTheta]", "HillasAntiSrc"); + fillasptheta.SetTitle("Task for 3D plots (alpha, E_est, Theta) (AntiSource)"); fillasptheta.SetFilter(&fsrc); @@ -237,7 +298,7 @@ // tasklist.AddToList(&reader); + + tasklist.AddToList(&lvl1); tasklist.AddToList(&rand); - - tasklist.AddToList(&acalc); tasklist.AddToList(&pcopy); @@ -246,5 +307,7 @@ tasklist.AddToList(&clean); tasklist.AddToList(&blind); + tasklist.AddToList(&hcalc); + tasklist.AddToList(&hsrc1); tasklist.AddToList(&hsrc2); @@ -253,4 +316,6 @@ tasklist.AddToList(&hfill1h); + tasklist.AddToList(&hfill1x); + tasklist.AddToList(&hfill1m); tasklist.AddToList(&hfill2s); @@ -263,5 +328,5 @@ tasklist.AddToList(&filldtimetheta); - // tasklist.AddToList(&lvl1); + // tasklist.AddToList(&fethetaall); @@ -296,16 +361,27 @@ tasklist.PrintStatistics(10); - /* - parlist.FindObject("HSource")->DrawClone();; - parlist.FindObject("MHHillas")->DrawClone(); - */ - /* - parlist.FindObject("MHStarMap")->DrawClone(); - */ - /* - parlist.FindObject("HAntiSource")->DrawClone(); - */ - /* - parlist.FindObject("MHMcCollectionArea")->DrawClone(); + + 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(); */ @@ -313,30 +389,15 @@ // 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(); */ - - //---------------------------------------------------------------------- - // 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 @@ -348,4 +409,5 @@ MHTimeDiffTime *dtimetime = (MHTimeDiffTime*)parlist.FindObject("EffOnTime"); MHTimeDiffTheta *dtimetheta = (MHTimeDiffTheta*)parlist.FindObject("EffOnTheta"); + /* dtimetime->DrawClone(); @@ -355,22 +417,42 @@ //.................................... - // calculate effective on time for different intervals in time + // calculate effective on time for different intervals in Time // and for different intervals in Theta - MHEffOnTimeTime effontime; - MHEffOnTimeTheta effontheta; + + char *varnametime; + varnametime = "Time"; + char *unittime; + unittime = "[ s]"; + MHEffOnTime effontime(varnametime, unittime); + + char *varnametheta; + varnametheta = "Theta"; + char *unittheta; + unittheta = "[ \\circ]"; + MHEffOnTime effontheta(varnametheta, unittheta); + effontime.SetupFill(&parlist); effontheta.SetupFill(&parlist); - effontime.Calc(dtimetime->GetHist()); - effontheta.Calc(dtimetheta->GetHist()); - - // plot effective on time versus time + + // 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 @@ -379,50 +461,39 @@ // 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; + 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(); */ - - /* - 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(); - */ + // 3rd and 4th argument are name and title of the resulting histogram + + MHGamma gamma; + + Bool_t Draw = kFALSE; + 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); + //---------------------------------------------------------------------- @@ -431,58 +502,134 @@ // 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); - + Axis_t lo = 0; // [deg] + Axis_t up = 10; // [deg] + Bool_t Draw = kFALSE; + 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 - parlist.FindObject("MHMcEnergyMigration")->DrawClone(); + + /* + 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 + Bool_t 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 + Bool_t Draw=kTRUE; + 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; - //---------------------------------------------------------------------- - 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! + //====================================================================== } //===========================================================================