Index: trunk/MagicSoft/Mars/macros/flux.C
===================================================================
--- trunk/MagicSoft/Mars/macros/flux.C	(revision 1283)
+++ trunk/MagicSoft/Mars/macros/flux.C	(revision 1294)
@@ -16,6 +16,5 @@
 !
 !
-!   Author(s): Thomas Bretz  1/2002 <mailto:tbretz@uni-sw.gwdg.de>
-!   Author(s): Wolfgang Wittek 1/2002 <mailto:wittek@mppmu.mpg.de>
+!   Author(s): Wolfgang Wittek 4/2002 <mailto:wittek@mppmu.mpg.de>
 !
 !   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     (<Theta>,  time)
+    //                       (<Theta>, 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!
 }
+//===========================================================================
+
+
