Changeset 1886

Timestamp:

04/01/03 17:49:19 (22 years ago)

Author:

moralejo

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• macros/CT1EnergyEst.C (modified) (15 diffs)
• manalysis/MEnergyEstParam.cc (modified) (3 diffs)
• manalysis/MEnergyEstParam.h (modified) (2 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -1 +1 @@
                                                  -*-*- END OF LINE -*-*-
+
+ 2003/04/01: Abelardo Moralejo
+
+   * macros/CT1EnergyEst.C
+     - added argument (maximum dist parameter), changed (reduced) output 
+       histograms, added writing to (and reading from) a file the energy 
+       estimation parameters and the histograms. Added comments.
+
+   * manalysis/MEnergyEstParam.[h,cc]
+     - added member function GetCoeff. Changed comment.
+
 
  2003/03/31: Thomas Bretz

trunk/MagicSoft/Mars/macros/CT1EnergyEst.C

@@ -16 +16 @@
 !
 !
-!   Author(s): Thomas Bretz et al,  09/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
+!   Author(s): Thomas Bretz,  09/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
+!              Abelardo Moralejo, 03/2003 <mailto:moralejo@pd.infn.it>
 !
 !   Copyright: MAGIC Software Development, 2000-2002
@@ -23 +24 @@
 \* ======================================================================== */
 
+//
+// fcn is the function to be minimized using TMinuit::Migrad
+//
 void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
 {
   MEvtLoop *evtloop = (MEvtLoop*)gMinuit->GetObjectFit();
 
-  MTaskList *tlist  = (MTaskList*)evtloop->GetParList()->FindObject("MTaskList"); //GetTaskList();
+  MTaskList *tlist  = (MTaskList*)evtloop->GetParList()->FindObject("MTaskList"); // GetTaskList();
 
   MChisqEval      *eval = (MChisqEval*)     tlist->FindObject("MChisqEval");
@@ -41 +45 @@
 //------------------------------------------------------------------------
 
-void CT1EnergyEst(Float_t maxhadronness=1.)
+void CT1EnergyEst(Float_t maxhadronness=0.22, Float_t maxdist = 1.)
 {
-  //  Bool_t evalenergy=kFALSE;
+  //
+  // Upper cut in Maxdist taken from D. Kranich's Ph D.
+  // We have to convert maxdist from degrees to mm:
+  //
+
+  MGeomCamCT1 gct1;
+  maxdist /= gct1.GetConvMm2Deg();
+
   //
   // Fill events into a MHMatrix
@@ -50 +61 @@
   MHMatrix matrix;
 
-  Int_t col = matrix.AddColumn("MMcEvt.fEnergy");
+  //
+  // colenergy and colimpact are the indexes of the matrix columns containing 
+  // the MC energy and impact parameter.
+  //
+  Int_t colenergy = matrix.AddColumn("MMcEvt.fEnergy");
+  Int_t colimpact = matrix.AddColumn("MMcEvt.fImpact");
 
   MEnergyEstParam eest("Hillas");
@@ -56 +72 @@
   eest.InitMapping(&matrix);
 
+  //
+  // Here we read in the events for the training. Second argument of AddFile is
+  // the number of events. With 2000 it is rather fast, but only the lowest zenith
+  // angles will enter in the minimization.
+  //
   MReadTree read("Events");
-  read.AddFile("MC_ON2.root", 200);
+  read.AddFile("MC_ON2.root",2000);
   read.DisableAutoScheme();
 
+  //
+  // Filter to keep the most gamma-like events:
+  //
   TString hcut("MHadronness.fHadronness<");
   hcut += maxhadronness;
+  hcut += " && HillasSrc.fDist < ";
+  hcut += maxdist;
 
   MF filterhadrons(hcut);
 
-  if (!matrix.Fill(&parlist, &read))
+  //
+  // Fill the matrix used later in the minimization loop by Minuit:
+  //
+  if (!matrix.Fill(&parlist, &read, &filterhadrons))
     return;
 
@@ -76 +105 @@
   MMatrixLoop loop(&matrix);
 
-
   MChisqEval eval;
-  eval.SetY1(new MDataElement(&matrix, col));
+
+  eval.SetY1(new MDataElement(&matrix, colenergy));
   eval.SetY2(new MDataMember("MEnergyEst.fEnergy"));
+
   eval.SetOwner();
 
@@ -97 +127 @@
 
   // Ready for: minuit.mnexcm("SET ERR", arglist, 1, ierflg)
+
   if (minuit.SetErrorDef(1))
     {
@@ -104 +135 @@
 
   //
-  // Set initial values
-  //
-  TArrayD fA(5);
-  fA[0] =1;//4916.4;   */-2414.75;
-  fA[1] =1;//149.549;  */   1134.28;
-  fA[2] =1;//-558.209; */   132.932;
-  fA[3] =1;//0.270725; */  0.292845;
-  fA[4] =1;//107.001;  */   107.001;
-  
+  // Set initial values of the parameters (close enough to the final ones, taken
+  // from previous runs of the procedure). Parameter fA[4] is not used in the default
+  // energy estimation model (from D. Kranich).
+  //
+  TArrayD fA(5);  
   TArrayD fB(7);
-  fB[0] = 1;//-8234.12; */ -4282.25;
-  fB[1] = 1;// 23.2153; */    18.892;
-  fB[2] = 1;// 0.416372;*/  0.193373;
-  fB[3] = 1;// 332.42;  */  203.803;
-  fB[4] = 1;// -0.701764;*/ -0.534876;
-  fB[5] = 1;//-0.0131774;*/ -0.00789539;
-  fB[6] = 1;//-0.162687;*/   0.111913;
-  
+
+  fA[0] =  21006.2;
+  fA[1] = -43.2648;
+  fA[2] = -690.403;
+  fA[3] = -0.0428544;
+  fA[4] =  1.;
+  fB[0] = -3391.05;
+  fB[1] =  136.58;
+  fB[2] =  0.253807;
+  fB[3] =  254.363;
+  fB[4] =  61.0386;
+  fB[5] = -0.0190984;
+  fB[6] = -421695;
+
+  //
   // Set starting values and step sizes for parameters
+  //
   for (Int_t i=0; i<fA.GetSize(); i++)
     {
@@ -166 +201 @@
   minuit.SetObjectFit(&evtloop);
 
-   
   TStopwatch clock;
   clock.Start();
 
-  // Now ready for minimization step: minuit.mnexcm("MIGRAD", arglist, 1, ierflg)
+  // Now ready for minimization step:
+
   gLog.SetNullOutput(kTRUE);
   Bool_t rc = minuit.Migrad();
   gLog.SetNullOutput(kFALSE);
-
+  
   if (rc)
     {
@@ -199 +234 @@
         }
 
-      cout << i << ":  " << val << "  +-  " << er << endl;
+      cout << i << ":  " << val << endl; 
+        //      "  +-  " << er << endl;
     }
 
@@ -208 +244 @@
     minuit.mnstat(amin, edm, errdef, nvpar, nparx, icstat);
     minuit.mnprin(3, amin);
+    eest.Print();
   */
-  eest.Print();
-
-
-
-  //---------------------------------------------------------------
-
+
+  eest.StopMapping();
+
+
+  //
+  // Now write the parameters of the energy estimator to file:
+  //
+  TVector* EnergyParams = new TVector(12);
+  for (Int_t i=0; i<eest.GetNumCoeff(); i++)
+    EnergyParams(i) = eest.GetCoeff(i);
+  TFile out("energyest.root", "recreate");
+  EnergyParams->Write("EnergyParams");
+  out.Close();
+
+  //
+  // End of Part 1 (estimation of the parameters)
+  //
+
+
+  ////////////////////////////////////////////////////////////////////////////
+  ////////////////////////////////////////////////////////////////////////////
+  ////////////////////////////////////////////////////////////////////////////
+  ////////////////////////////////////////////////////////////////////////////
+
+
+  //
   // Part 2: Now test how the energy estimation method works.
   //
@@ -222 +279 @@
   //
 
-
   MTaskList tlist2;
-  parlist.Replace(&tlist2);
+  MParList  parlist2;  
+
+  parlist2.AddToList(&tlist2);
 
   //
@@ -231 +289 @@
   //
 
-
-  read->SetEventNum(0);
-
-  //
-  // Use this to change the binnign of the histograms to CT1-style
-  //
-  Bool_t usect1 = kTRUE;
-
-  MH3 mh3e("MMcEvt.fEnergy",     "(MEnergyEst.fEnergy/MMcEvt.fEnergy-1)*(MEnergyEst.fEnergy/MMcEvt.fEnergy-1)");
-  MH3 mh3i("MMcEvt.fImpact/100", "(MEnergyEst.fImpact/MMcEvt.fImpact-1)*(MEnergyEst.fImpact/MMcEvt.fImpact-1)");
-  MH3 mh3eo("MMcEvt.fEnergy",     "MEnergyEst.fEnergy/MMcEvt.fEnergy-1");
-  MH3 mh3io("MMcEvt.fImpact/100", "MEnergyEst.fImpact/MMcEvt.fImpact-1");
-
-  MH3 mh3e2("MEnergyEst.fEnergy",     "(MEnergyEst.fEnergy/MMcEvt.fEnergy-1)*(MEnergyEst.fEnergy/MMcEvt.fEnergy-1)");
-  MH3 mh3i2("MEnergyEst.fImpact/100", "(MEnergyEst.fImpact/MMcEvt.fImpact-1)*(MEnergyEst.fImpact/MMcEvt.fImpact-1)");
-  MH3 mh3eo2("MEnergyEst.fEnergy",     "MEnergyEst.fEnergy/MMcEvt.fEnergy-1");
-  MH3 mh3io2("MEnergyEst.fImpact/100", "MEnergyEst.fImpact/MMcEvt.fImpact-1");
-
-  MH3 mhe("MMcEvt.fEnergy",     "MEnergyEst.fEnergy");
-  MH3 mhi("MMcEvt.fImpact/100", "MEnergyEst.fImpact/100");
-
-  mh3e.SetName("HistEnergy");
-  mh3i.SetName("HistImpact");
-  mh3eo.SetName("HistEnergyOffset");
-  mh3io.SetName("HistImpactOffset");
-
-  mh3e2.SetName("HistEnergy");
-  mh3i2.SetName("HistImpact");
-  mh3eo2.SetName("HistEnergyOffset");
-  mh3io2.SetName("HistImpactOffset");
+  MReadTree read2("Events");
+  read2.AddFile("MC_ON2.root");
+  read2.DisableAutoScheme();
+
+  //
+  // Create some histograms to display the results:
+  //
+
+  MH3 mh3ed("log10(MMcEvt.fEnergy)", "MEnergyEst.fEnergy/MMcEvt.fEnergy-1");
+  MH3 mh3ed2("log10(MEnergyEst.fEnergy)", "MEnergyEst.fEnergy/MMcEvt.fEnergy-1");
+  MH3 mhe("log10(MMcEvt.fEnergy)", "log10(MEnergyEst.fEnergy)");
+
+  MFillH hfilled(&mh3ed);
+  MFillH hfilled2(&mh3ed2);
+  MFillH hfillee(&mhe);
 
   mhe.SetName("HistEE");
-  mhi.SetName("HistII");
-
-  MFillH hfille(&mh3e);
-  MFillH hfilli(&mh3i);
-  MFillH hfilleo(&mh3eo);
-  MFillH hfillio(&mh3io);
-
-  MFillH hfille2(&mh3e2);
-  MFillH hfilli2(&mh3i2);
-  MFillH hfilleo2(&mh3eo2);
-  MFillH hfillio2(&mh3io2);
-
-  MFillH hfillee(&mhe);
-  MFillH hfillii(&mhi);
-
-  MBinning binsex("BinningHistEnergyX");
-  MBinning binsey("BinningHistEnergyY");
-  MBinning binsix("BinningHistImpactX");
-  MBinning binsiy("BinningHistImpactY");
-  MBinning binseox("BinningHistEnergyOffsetX");
-  MBinning binseoy("BinningHistEnergyOffsetY");
-  MBinning binsiox("BinningHistImpactOffsetX");
-  MBinning binsioy("BinningHistImpactOffsetY");
+  mh3ed.SetName("HistEnergyDelta");
+  mh3ed2.SetName("HistEnergyDelta");
+
+  MBinning binsedx("BinningHistEnergyDeltaX");
+  MBinning binsedy("BinningHistEnergyDeltaY");
   MBinning binseex("BinningHistEEX");
-  MBinning binsiix("BinningHistIIX");
   MBinning binseey("BinningHistEEY");
-  MBinning binsiiy("BinningHistIIY");
-
-  binsex.SetEdgesLog(50, usect1 ? 300: 10, usect1 ? 50000 : 1e4);
-  binsey.SetEdges(50, 0, usect1 ? 0.8 : 1.75);
-  binseox.SetEdgesLog(50, usect1 ? 300 : 10, usect1 ? 50000 : 1e4);
-  binseoy.SetEdges(50, usect1 ? -0.75 : -1.75, usect1 ? 0.75 : 1.75);
-
-  binsix.SetEdges(50, 0, usect1 ? 275 : 300);
-  binsiy.SetEdges(50, 0, usect1 ? 0.2 : 1.75);
-  binsiox.SetEdges(50, 0, usect1 ? 275 : 300);
-  binsioy.SetEdges(50, usect1 ? -0.75 : -1.75, usect1 ? 0.75 : 1.75);
-
-  binseex.SetEdgesLog(50, usect1 ? 300 : 10, usect1 ? 50000 : 15e3);
-  binseey.SetEdgesLog(50, usect1 ? 300 : 1,  usect1 ? 50000 : 2e3);
-  binsiix.SetEdges(50, 0, usect1 ? 275 : 300);
-  binsiiy.SetEdges(50, 0, usect1 ? 275 : 150);
-
-  parlist.AddToList(&binsex);
-  parlist.AddToList(&binsey);
-  parlist.AddToList(&binsix);
-  parlist.AddToList(&binsiy);
-  parlist.AddToList(&binseox);
-  parlist.AddToList(&binseoy);
-  parlist.AddToList(&binsiox);
-  parlist.AddToList(&binsioy);
-  parlist.AddToList(&binseex);
-  parlist.AddToList(&binseey);
-  parlist.AddToList(&binsiix);
-  parlist.AddToList(&binsiiy);
-
-  eest.StopMapping();
-  eest.Add("HillasSrc");
+
+  binsedx.SetEdges(25, 2.5, 5.);
+  binsedy.SetEdges(30, -1., 1.);
+  binseex.SetEdges(25, 2.5, 5.);
+  binseey.SetEdges(30, 2., 5.);
+
+  parlist2.AddToList(&binsedx);
+  parlist2.AddToList(&binsedy);
+  parlist2.AddToList(&binseex);
+  parlist2.AddToList(&binseey);
+
+  //
+  // Create energy estimator task, add necessary containers and 
+  // initialize parameters read from file:
+  //
+
+  MEnergyEstParam eest2("Hillas");
+  eest2.Add("HillasSrc");
+
+  TFile enparam("energyest.root");
+  TArrayD parA(5);
+  TArrayD parB(7);
+
+  for (Int_t i=0; i<parA.GetSize(); i++)
+    parA[i] = EnergyParams(i);
+  for (Int_t i=0; i<parB.GetSize(); i++)
+    parB[i] = EnergyParams(i+parA.GetSize());
+
+  eest2.SetCoeffA(parA);
+  eest2.SetCoeffB(parB);
+
+  enparam.Close();
 
   //
   //  Setup tasklists
   //
-  tlist2.AddToList(&read);
-
+
+  tlist2.AddToList(&read2);
+
+  //
+  // Include filter on hadronness and maximum value of dist:
+  //
   TString hcut2("MHadronness.fHadronness>");
   hcut2 += maxhadronness;
+  hcut2 += "|| HillasSrc.fDist > ";
+  hcut2 += maxdist;
+
   MContinue cont(hcut2);
 
   tlist2.AddToList(&cont);
-
-  tlist2.AddToList(&eest);
+  tlist2.AddToList(&eest2);
+
   // tlist2.AddToList(new MPrint("MMcEvt"));
   // tlist2.AddToList(new MPrint("MEnergyEst"));
 
-  tlist2.AddToList(&hfille);
-  tlist2.AddToList(&hfilli);
-  tlist2.AddToList(&hfilleo);
-  tlist2.AddToList(&hfillio);
-
-  tlist2.AddToList(&hfille2);
-  tlist2.AddToList(&hfilli2);
-  tlist2.AddToList(&hfilleo2);
-  tlist2.AddToList(&hfillio2);
-
+  tlist2.AddToList(&hfilled);
+  tlist2.AddToList(&hfilled2);
   tlist2.AddToList(&hfillee);
-  tlist2.AddToList(&hfillii);
 
   //
@@ -357 +379 @@
   MEvtLoop evtloop2;
   evtloop2.SetProgressBar(&bar);
-  evtloop2.SetParList(&parlist);
+  evtloop2.SetParList(&parlist2);
 
   //
@@ -365 +387 @@
     return;
 
-  tlist2.PrintStatistics();
-
-  const TString text = "\\sqrt{<y>}=%.0f%%";
-
-  char txt[1000];
-
-  TCanvas *c=new TCanvas("Est1", "Estimates vs. E_{true}");
+  //
+  // Plot results:
+  //
+  mhe.GetHist()->GetXaxis()->SetTitle("log_{10} E_{MC} (GeV)");
+  mhe.GetHist()->GetYaxis()->SetTitle("log_{10} E_{EST} (GeV)");
+  mhe.GetHist()->GetXaxis()->SetTitleOffset(1.2);
+
+  TCanvas *c = new TCanvas("energy","CT1 Energy estimation");
   c->Divide(2,2);
   c->cd(1);
-  mh3i.DrawClone("PROFXnonew");
-  sprintf(txt, text.Data(), sqrt(mh3i.GetHist().GetMean(2))*100);
-  TLatex *t = new TLatex(180, 0.15, txt);
-  t->Draw();
+  energy_1->SetBottomMargin(0.12);
+  mhe.DrawClone("PROFXnonew");
+
+  mh3ed.GetHist()->GetXaxis()->SetTitle("log_{10} E_{MC} (GeV)");
+  mh3ed.GetHist()->GetYaxis()->SetTitle("\\frac{E_{EST} - E_{MC}}{E_{MC}}");
+  mh3ed.GetHist()->GetXaxis()->SetTitleOffset(1.2);
+  mh3ed.GetHist()->GetYaxis()->SetTitleOffset(1.5);
   c->cd(2);
-  mh3e.DrawClone("PROFXnonew");
-  sprintf(txt, text.Data(), sqrt(mh3e.GetHist().GetMean(2))*100);
-  t = new TLatex(3.5, 0.6, txt);
-  t->Draw();
+  energy_2->SetLeftMargin(0.15);
+  energy_2->SetBottomMargin(0.12);
+  mh3ed.DrawClone("PROFXnonew");
+
+  mh3ed2.GetHist()->GetXaxis()->SetTitle("log_{10} E_{EST} (GeV)");
+  mh3ed2.GetHist()->GetYaxis()->SetTitle("\\frac{E_{EST} - E_{MC}}{E_{MC}}");
+  mh3ed2.GetHist()->GetXaxis()->SetTitleOffset(1.2);
+  mh3ed2.GetHist()->GetYaxis()->SetTitleOffset(1.5);
   c->cd(3);
-  mh3io.DrawClone("PROFXnonew");
+  energy_3->SetLeftMargin(0.15);
+  energy_3->SetBottomMargin(0.12);
+  mh3ed2.DrawClone("PROFXnonew");
+
+  ((TH2*)mh3ed2.GetHist())->ProjectionY("deltae");
+
   c->cd(4);
-  mh3eo.DrawClone("PROFXnonew");
-
-  c=new TCanvas("Est2", "Estimates vs. E_{est}");
-  c->Divide(2,2);
-  c->cd(1);
-  mh3i2.DrawClone("PROFXnonew");
-  c->cd(2);
-  mh3e2.DrawClone("PROFXnonew");
-  c->cd(3);
-  mh3io2.DrawClone("PROFXnonew");
-  c->cd(4);
-  mh3eo2.DrawClone("PROFXnonew");
-
-  mhe.DrawClone("PROFX");
-  mhi.DrawClone("PROFX");
+  energy_4->SetLeftMargin(0.1);
+  energy_4->SetRightMargin(0.05);
+  energy_4->SetBottomMargin(0.12);
+  deltae.GetXaxis()->SetTitleOffset(1.2);
+  deltae.GetXaxis()->SetTitle("(E_{EST} - E_{MC}) / E_{MC}");
+  deltae.Draw("e");
+  
+  TFile out2("energyest.root", "update");
+  ((TH2*)mh3ed.GetHist())->Write("mh3ed");
+  ((TH2*)mh3ed2.GetHist())->Write("mh3ed2");
+  ((TH2*)mhe.GetHist())->Write("mhe");
+  deltae.Write();
+  out2.Close();
+
+  return;
+
 }

trunk/MagicSoft/Mars/manalysis/MEnergyEstParam.cc

@@ -58 +58 @@
     fA[3] =  0.006;
     fA[4] =   58.3;
+
     /*
-     fA[0] =   39.667402; // [cm]
-     fA[1] =  143.081912; // [cm]
-     fA[2] = -395.501677; // [cm]
-     fA[3] =    0.006193;
-     fA[4] =    0;
-     */
+      fA[0] =   39.667402; // [cm]
+      fA[1] =  143.081912; // [cm]
+      fA[2] = -395.501677; // [cm]
+      fA[3] =    0.006193;
+      fA[4] =    0;
+    */
 
     fB.Set(7);
@@ -74 +75 @@
     fB[5] =  4.13e-8;   // [GeV]
     fB[6] =     0.05;
+
     /*
-     fB[0] =   45.037867; // [GeV]
-     fB[1] =    0.087222; // [GeV]
-     fB[2] =   -0.208065; // [GeV/cm]
-     fB[3] =   78.164942; // [GeV]
-     fB[4] = -159.361283; // [GeV]
-     fB[5] =   -0.130455; // [GeV/cm]
-     fB[6] =    0.051139;
-     */
+      fB[0] =   45.037867; // [GeV]
+      fB[1] =    0.087222; // [GeV]
+      fB[2] =   -0.208065; // [GeV/cm]
+      fB[3] =   78.164942; // [GeV]
+      fB[4] = -159.361283; // [GeV]
+      fB[5] =   -0.130455; // [GeV/cm]
+      fB[6] =    0.051139;
+    */
 }
 
@@ -340 +342 @@
         const Double_t dist = fMatrix ? GetVal(col++) : ((MHillasSrc*)NextH())->GetDist();
 
-
         /* DANIEL - MARCOS */
         const Double_t ir = i0 * (i1 + fA[1]*dist); // [cm]
         Double_t er = e0 * (e1 + e2*ir);       // [GeV]
 
-
-        /* MY PARAM */
+        /* THOMAS BRETZ */
         // if (width==0) return kCONTINUE;
         // const Double_t ir = i0 * (i1 + dist*(fA[1]/width + fA[4]/log10(size))); // [cm]

trunk/MagicSoft/Mars/manalysis/MEnergyEstParam.h

@@ -40 +40 @@
     Double_t GetVal(Int_t i) const;
 
+
 public:
     MEnergyEstParam(const char *hil="MHillas", const char *name=NULL, const char *title=NULL);
@@ -58 +59 @@
     void SetCoeffB(const TArrayD &arr);
 
+    Double_t GetCoeff(Int_t i) { return i<fA.GetSize()? fA[i] : fB[i-fA.GetSize()]; }
+
     void Print(Option_t *o=NULL);