Ignore:
Timestamp:
04/01/03 17:49:19 (22 years ago)
Author:
moralejo
Message:
*** empty log message ***
File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/macros/CT1EnergyEst.C

    r1852 r1886  
    1616!
    1717!
    18 !   Author(s): Thomas Bretz et al,  09/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
     18!   Author(s): Thomas Bretz,  09/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
     19!              Abelardo Moralejo, 03/2003 <mailto:moralejo@pd.infn.it>
    1920!
    2021!   Copyright: MAGIC Software Development, 2000-2002
    … …  
    2324\* ======================================================================== */
    2425
     26//
     27// fcn is the function to be minimized using TMinuit::Migrad
     28//
    2529void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
    2630{
    2731  MEvtLoop *evtloop = (MEvtLoop*)gMinuit->GetObjectFit();
    2832
    29   MTaskList *tlist  = (MTaskList*)evtloop->GetParList()->FindObject("MTaskList"); //GetTaskList();
     33  MTaskList *tlist  = (MTaskList*)evtloop->GetParList()->FindObject("MTaskList"); // GetTaskList();
    3034
    3135  MChisqEval      *eval = (MChisqEval*)     tlist->FindObject("MChisqEval");
    … …  
    4145//------------------------------------------------------------------------
    4246
    43 void CT1EnergyEst(Float_t maxhadronness=1.)
     47void CT1EnergyEst(Float_t maxhadronness=0.22, Float_t maxdist = 1.)
    4448{
    45   //  Bool_t evalenergy=kFALSE;
     49  //
     50  // Upper cut in Maxdist taken from D. Kranich's Ph D.
     51  // We have to convert maxdist from degrees to mm:
     52  //
     53
     54  MGeomCamCT1 gct1;
     55  maxdist /= gct1.GetConvMm2Deg();
     56
    4657  //
    4758  // Fill events into a MHMatrix
    … …  
    5061  MHMatrix matrix;
    5162
    52   Int_t col = matrix.AddColumn("MMcEvt.fEnergy");
     63  //
     64  // colenergy and colimpact are the indexes of the matrix columns containing
     65  // the MC energy and impact parameter.
     66  //
     67  Int_t colenergy = matrix.AddColumn("MMcEvt.fEnergy");
     68  Int_t colimpact = matrix.AddColumn("MMcEvt.fImpact");
    5369
    5470  MEnergyEstParam eest("Hillas");
    … …  
    5672  eest.InitMapping(&matrix);
    5773
     74  //
     75  // Here we read in the events for the training. Second argument of AddFile is
     76  // the number of events. With 2000 it is rather fast, but only the lowest zenith
     77  // angles will enter in the minimization.
     78  //
    5879  MReadTree read("Events");
    59   read.AddFile("MC_ON2.root", 200);
     80  read.AddFile("MC_ON2.root",2000);
    6081  read.DisableAutoScheme();
    6182
     83  //
     84  // Filter to keep the most gamma-like events:
     85  //
    6286  TString hcut("MHadronness.fHadronness<");
    6387  hcut += maxhadronness;
     88  hcut += " && HillasSrc.fDist < ";
     89  hcut += maxdist;
    6490
    6591  MF filterhadrons(hcut);
    6692
    67   if (!matrix.Fill(&parlist, &read))
     93  //
     94  // Fill the matrix used later in the minimization loop by Minuit:
     95  //
     96  if (!matrix.Fill(&parlist, &read, &filterhadrons))
    6897    return;
    6998
    … …  
    76105  MMatrixLoop loop(&matrix);
    77106
    78 
    79107  MChisqEval eval;
    80   eval.SetY1(new MDataElement(&matrix, col));
     108
     109  eval.SetY1(new MDataElement(&matrix, colenergy));
    81110  eval.SetY2(new MDataMember("MEnergyEst.fEnergy"));
     111
    82112  eval.SetOwner();
    83113
    … …  
    97127
    98128  // Ready for: minuit.mnexcm("SET ERR", arglist, 1, ierflg)
     129
    99130  if (minuit.SetErrorDef(1))
    100131    {
    … …  
    104135
    105136  //
    106   // Set initial values
    107   //
    108   TArrayD fA(5);
    109   fA[0] =1;//4916.4;   */-2414.75;
    110   fA[1] =1;//149.549;  */   1134.28;
    111   fA[2] =1;//-558.209; */   132.932;
    112   fA[3] =1;//0.270725; */  0.292845;
    113   fA[4] =1;//107.001;  */   107.001;
    114  
     137  // Set initial values of the parameters (close enough to the final ones, taken
     138  // from previous runs of the procedure). Parameter fA[4] is not used in the default
     139  // energy estimation model (from D. Kranich).
     140  //
     141  TArrayD fA(5); 
    115142  TArrayD fB(7);
    116   fB[0] = 1;//-8234.12; */ -4282.25;
    117   fB[1] = 1;// 23.2153; */    18.892;
    118   fB[2] = 1;// 0.416372;*/  0.193373;
    119   fB[3] = 1;// 332.42;  */  203.803;
    120   fB[4] = 1;// -0.701764;*/ -0.534876;
    121   fB[5] = 1;//-0.0131774;*/ -0.00789539;
    122   fB[6] = 1;//-0.162687;*/   0.111913;
    123  
     143
     144  fA[0] =  21006.2;
     145  fA[1] = -43.2648;
     146  fA[2] = -690.403;
     147  fA[3] = -0.0428544;
     148  fA[4] =  1.;
     149  fB[0] = -3391.05;
     150  fB[1] =  136.58;
     151  fB[2] =  0.253807;
     152  fB[3] =  254.363;
     153  fB[4] =  61.0386;
     154  fB[5] = -0.0190984;
     155  fB[6] = -421695;
     156
     157  //
    124158  // Set starting values and step sizes for parameters
     159  //
    125160  for (Int_t i=0; i<fA.GetSize(); i++)
    126161    {
    … …  
    166201  minuit.SetObjectFit(&evtloop);
    167202
    168    
    169203  TStopwatch clock;
    170204  clock.Start();
    171205
    172   // Now ready for minimization step: minuit.mnexcm("MIGRAD", arglist, 1, ierflg)
     206  // Now ready for minimization step:
     207
    173208  gLog.SetNullOutput(kTRUE);
    174209  Bool_t rc = minuit.Migrad();
    175210  gLog.SetNullOutput(kFALSE);
    176 
     211 
    177212  if (rc)
    178213    {
    … …  
    199234        }
    200235
    201       cout << i << ":  " << val << "  +-  " << er << endl;
     236      cout << i << ":  " << val << endl;
     237        //      "  +-  " << er << endl;
    202238    }
    203239
    … …  
    208244    minuit.mnstat(amin, edm, errdef, nvpar, nparx, icstat);
    209245    minuit.mnprin(3, amin);
     246    eest.Print();
    210247  */
    211   eest.Print();
    212 
    213 
    214 
    215   //---------------------------------------------------------------
    216 
     248
     249  eest.StopMapping();
     250
     251
     252  //
     253  // Now write the parameters of the energy estimator to file:
     254  //
     255  TVector* EnergyParams = new TVector(12);
     256  for (Int_t i=0; i<eest.GetNumCoeff(); i++)
     257    EnergyParams(i) = eest.GetCoeff(i);
     258  TFile out("energyest.root", "recreate");
     259  EnergyParams->Write("EnergyParams");
     260  out.Close();
     261
     262  //
     263  // End of Part 1 (estimation of the parameters)
     264  //
     265
     266
     267  ////////////////////////////////////////////////////////////////////////////
     268  ////////////////////////////////////////////////////////////////////////////
     269  ////////////////////////////////////////////////////////////////////////////
     270  ////////////////////////////////////////////////////////////////////////////
     271
     272
     273  //
    217274  // Part 2: Now test how the energy estimation method works.
    218275  //
    … …  
    222279  //
    223280
    224 
    225281  MTaskList tlist2;
    226   parlist.Replace(&tlist2);
     282  MParList  parlist2; 
     283
     284  parlist2.AddToList(&tlist2);
    227285
    228286  //
    … …  
    231289  //
    232290
    233 
    234   read->SetEventNum(0);
    235 
    236   //
    237   // Use this to change the binnign of the histograms to CT1-style
    238   //
    239   Bool_t usect1 = kTRUE;
    240 
    241   MH3 mh3e("MMcEvt.fEnergy",     "(MEnergyEst.fEnergy/MMcEvt.fEnergy-1)*(MEnergyEst.fEnergy/MMcEvt.fEnergy-1)");
    242   MH3 mh3i("MMcEvt.fImpact/100", "(MEnergyEst.fImpact/MMcEvt.fImpact-1)*(MEnergyEst.fImpact/MMcEvt.fImpact-1)");
    243   MH3 mh3eo("MMcEvt.fEnergy",     "MEnergyEst.fEnergy/MMcEvt.fEnergy-1");
    244   MH3 mh3io("MMcEvt.fImpact/100", "MEnergyEst.fImpact/MMcEvt.fImpact-1");
    245 
    246   MH3 mh3e2("MEnergyEst.fEnergy",     "(MEnergyEst.fEnergy/MMcEvt.fEnergy-1)*(MEnergyEst.fEnergy/MMcEvt.fEnergy-1)");
    247   MH3 mh3i2("MEnergyEst.fImpact/100", "(MEnergyEst.fImpact/MMcEvt.fImpact-1)*(MEnergyEst.fImpact/MMcEvt.fImpact-1)");
    248   MH3 mh3eo2("MEnergyEst.fEnergy",     "MEnergyEst.fEnergy/MMcEvt.fEnergy-1");
    249   MH3 mh3io2("MEnergyEst.fImpact/100", "MEnergyEst.fImpact/MMcEvt.fImpact-1");
    250 
    251   MH3 mhe("MMcEvt.fEnergy",     "MEnergyEst.fEnergy");
    252   MH3 mhi("MMcEvt.fImpact/100", "MEnergyEst.fImpact/100");
    253 
    254   mh3e.SetName("HistEnergy");
    255   mh3i.SetName("HistImpact");
    256   mh3eo.SetName("HistEnergyOffset");
    257   mh3io.SetName("HistImpactOffset");
    258 
    259   mh3e2.SetName("HistEnergy");
    260   mh3i2.SetName("HistImpact");
    261   mh3eo2.SetName("HistEnergyOffset");
    262   mh3io2.SetName("HistImpactOffset");
     291  MReadTree read2("Events");
     292  read2.AddFile("MC_ON2.root");
     293  read2.DisableAutoScheme();
     294
     295  //
     296  // Create some histograms to display the results:
     297  //
     298
     299  MH3 mh3ed("log10(MMcEvt.fEnergy)", "MEnergyEst.fEnergy/MMcEvt.fEnergy-1");
     300  MH3 mh3ed2("log10(MEnergyEst.fEnergy)", "MEnergyEst.fEnergy/MMcEvt.fEnergy-1");
     301  MH3 mhe("log10(MMcEvt.fEnergy)", "log10(MEnergyEst.fEnergy)");
     302
     303  MFillH hfilled(&mh3ed);
     304  MFillH hfilled2(&mh3ed2);
     305  MFillH hfillee(&mhe);
    263306
    264307  mhe.SetName("HistEE");
    265   mhi.SetName("HistII");
    266 
    267   MFillH hfille(&mh3e);
    268   MFillH hfilli(&mh3i);
    269   MFillH hfilleo(&mh3eo);
    270   MFillH hfillio(&mh3io);
    271 
    272   MFillH hfille2(&mh3e2);
    273   MFillH hfilli2(&mh3i2);
    274   MFillH hfilleo2(&mh3eo2);
    275   MFillH hfillio2(&mh3io2);
    276 
    277   MFillH hfillee(&mhe);
    278   MFillH hfillii(&mhi);
    279 
    280   MBinning binsex("BinningHistEnergyX");
    281   MBinning binsey("BinningHistEnergyY");
    282   MBinning binsix("BinningHistImpactX");
    283   MBinning binsiy("BinningHistImpactY");
    284   MBinning binseox("BinningHistEnergyOffsetX");
    285   MBinning binseoy("BinningHistEnergyOffsetY");
    286   MBinning binsiox("BinningHistImpactOffsetX");
    287   MBinning binsioy("BinningHistImpactOffsetY");
     308  mh3ed.SetName("HistEnergyDelta");
     309  mh3ed2.SetName("HistEnergyDelta");
     310
     311  MBinning binsedx("BinningHistEnergyDeltaX");
     312  MBinning binsedy("BinningHistEnergyDeltaY");
    288313  MBinning binseex("BinningHistEEX");
    289   MBinning binsiix("BinningHistIIX");
    290314  MBinning binseey("BinningHistEEY");
    291   MBinning binsiiy("BinningHistIIY");
    292 
    293   binsex.SetEdgesLog(50, usect1 ? 300: 10, usect1 ? 50000 : 1e4);
    294   binsey.SetEdges(50, 0, usect1 ? 0.8 : 1.75);
    295   binseox.SetEdgesLog(50, usect1 ? 300 : 10, usect1 ? 50000 : 1e4);
    296   binseoy.SetEdges(50, usect1 ? -0.75 : -1.75, usect1 ? 0.75 : 1.75);
    297 
    298   binsix.SetEdges(50, 0, usect1 ? 275 : 300);
    299   binsiy.SetEdges(50, 0, usect1 ? 0.2 : 1.75);
    300   binsiox.SetEdges(50, 0, usect1 ? 275 : 300);
    301   binsioy.SetEdges(50, usect1 ? -0.75 : -1.75, usect1 ? 0.75 : 1.75);
    302 
    303   binseex.SetEdgesLog(50, usect1 ? 300 : 10, usect1 ? 50000 : 15e3);
    304   binseey.SetEdgesLog(50, usect1 ? 300 : 1,  usect1 ? 50000 : 2e3);
    305   binsiix.SetEdges(50, 0, usect1 ? 275 : 300);
    306   binsiiy.SetEdges(50, 0, usect1 ? 275 : 150);
    307 
    308   parlist.AddToList(&binsex);
    309   parlist.AddToList(&binsey);
    310   parlist.AddToList(&binsix);
    311   parlist.AddToList(&binsiy);
    312   parlist.AddToList(&binseox);
    313   parlist.AddToList(&binseoy);
    314   parlist.AddToList(&binsiox);
    315   parlist.AddToList(&binsioy);
    316   parlist.AddToList(&binseex);
    317   parlist.AddToList(&binseey);
    318   parlist.AddToList(&binsiix);
    319   parlist.AddToList(&binsiiy);
    320 
    321   eest.StopMapping();
    322   eest.Add("HillasSrc");
     315
     316  binsedx.SetEdges(25, 2.5, 5.);
     317  binsedy.SetEdges(30, -1., 1.);
     318  binseex.SetEdges(25, 2.5, 5.);
     319  binseey.SetEdges(30, 2., 5.);
     320
     321  parlist2.AddToList(&binsedx);
     322  parlist2.AddToList(&binsedy);
     323  parlist2.AddToList(&binseex);
     324  parlist2.AddToList(&binseey);
     325
     326  //
     327  // Create energy estimator task, add necessary containers and
     328  // initialize parameters read from file:
     329  //
     330
     331  MEnergyEstParam eest2("Hillas");
     332  eest2.Add("HillasSrc");
     333
     334  TFile enparam("energyest.root");
     335  TArrayD parA(5);
     336  TArrayD parB(7);
     337
     338  for (Int_t i=0; i<parA.GetSize(); i++)
     339    parA[i] = EnergyParams(i);
     340  for (Int_t i=0; i<parB.GetSize(); i++)
     341    parB[i] = EnergyParams(i+parA.GetSize());
     342
     343  eest2.SetCoeffA(parA);
     344  eest2.SetCoeffB(parB);
     345
     346  enparam.Close();
    323347
    324348  //
    325349  //  Setup tasklists
    326350  //
    327   tlist2.AddToList(&read);
    328 
     351
     352  tlist2.AddToList(&read2);
     353
     354  //
     355  // Include filter on hadronness and maximum value of dist:
     356  //
    329357  TString hcut2("MHadronness.fHadronness>");
    330358  hcut2 += maxhadronness;
     359  hcut2 += "|| HillasSrc.fDist > ";
     360  hcut2 += maxdist;
     361
    331362  MContinue cont(hcut2);
    332363
    333364  tlist2.AddToList(&cont);
    334 
    335   tlist2.AddToList(&eest);
     365  tlist2.AddToList(&eest2);
     366
    336367  // tlist2.AddToList(new MPrint("MMcEvt"));
    337368  // tlist2.AddToList(new MPrint("MEnergyEst"));
    338369
    339   tlist2.AddToList(&hfille);
    340   tlist2.AddToList(&hfilli);
    341   tlist2.AddToList(&hfilleo);
    342   tlist2.AddToList(&hfillio);
    343 
    344   tlist2.AddToList(&hfille2);
    345   tlist2.AddToList(&hfilli2);
    346   tlist2.AddToList(&hfilleo2);
    347   tlist2.AddToList(&hfillio2);
    348 
     370  tlist2.AddToList(&hfilled);
     371  tlist2.AddToList(&hfilled2);
    349372  tlist2.AddToList(&hfillee);
    350   tlist2.AddToList(&hfillii);
    351373
    352374  //
    … …  
    357379  MEvtLoop evtloop2;
    358380  evtloop2.SetProgressBar(&bar);
    359   evtloop2.SetParList(&parlist);
     381  evtloop2.SetParList(&parlist2);
    360382
    361383  //
    … …  
    365387    return;
    366388
    367   tlist2.PrintStatistics();
    368 
    369   const TString text = "\\sqrt{<y>}=%.0f%%";
    370 
    371   char txt[1000];
    372 
    373   TCanvas *c=new TCanvas("Est1", "Estimates vs. E_{true}");
     389  //
     390  // Plot results:
     391  //
     392  mhe.GetHist()->GetXaxis()->SetTitle("log_{10} E_{MC} (GeV)");
     393  mhe.GetHist()->GetYaxis()->SetTitle("log_{10} E_{EST} (GeV)");
     394  mhe.GetHist()->GetXaxis()->SetTitleOffset(1.2);
     395
     396  TCanvas *c = new TCanvas("energy","CT1 Energy estimation");
    374397  c->Divide(2,2);
    375398  c->cd(1);
    376   mh3i.DrawClone("PROFXnonew");
    377   sprintf(txt, text.Data(), sqrt(mh3i.GetHist().GetMean(2))*100);
    378   TLatex *t = new TLatex(180, 0.15, txt);
    379   t->Draw();
     399  energy_1->SetBottomMargin(0.12);
     400  mhe.DrawClone("PROFXnonew");
     401
     402  mh3ed.GetHist()->GetXaxis()->SetTitle("log_{10} E_{MC} (GeV)");
     403  mh3ed.GetHist()->GetYaxis()->SetTitle("\\frac{E_{EST} - E_{MC}}{E_{MC}}");
     404  mh3ed.GetHist()->GetXaxis()->SetTitleOffset(1.2);
     405  mh3ed.GetHist()->GetYaxis()->SetTitleOffset(1.5);
    380406  c->cd(2);
    381   mh3e.DrawClone("PROFXnonew");
    382   sprintf(txt, text.Data(), sqrt(mh3e.GetHist().GetMean(2))*100);
    383   t = new TLatex(3.5, 0.6, txt);
    384   t->Draw();
     407  energy_2->SetLeftMargin(0.15);
     408  energy_2->SetBottomMargin(0.12);
     409  mh3ed.DrawClone("PROFXnonew");
     410
     411  mh3ed2.GetHist()->GetXaxis()->SetTitle("log_{10} E_{EST} (GeV)");
     412  mh3ed2.GetHist()->GetYaxis()->SetTitle("\\frac{E_{EST} - E_{MC}}{E_{MC}}");
     413  mh3ed2.GetHist()->GetXaxis()->SetTitleOffset(1.2);
     414  mh3ed2.GetHist()->GetYaxis()->SetTitleOffset(1.5);
    385415  c->cd(3);
    386   mh3io.DrawClone("PROFXnonew");
     416  energy_3->SetLeftMargin(0.15);
     417  energy_3->SetBottomMargin(0.12);
     418  mh3ed2.DrawClone("PROFXnonew");
     419
     420  ((TH2*)mh3ed2.GetHist())->ProjectionY("deltae");
     421
    387422  c->cd(4);
    388   mh3eo.DrawClone("PROFXnonew");
    389 
    390   c=new TCanvas("Est2", "Estimates vs. E_{est}");
    391   c->Divide(2,2);
    392   c->cd(1);
    393   mh3i2.DrawClone("PROFXnonew");
    394   c->cd(2);
    395   mh3e2.DrawClone("PROFXnonew");
    396   c->cd(3);
    397   mh3io2.DrawClone("PROFXnonew");
    398   c->cd(4);
    399   mh3eo2.DrawClone("PROFXnonew");
    400 
    401   mhe.DrawClone("PROFX");
    402   mhi.DrawClone("PROFX");
     423  energy_4->SetLeftMargin(0.1);
     424  energy_4->SetRightMargin(0.05);
     425  energy_4->SetBottomMargin(0.12);
     426  deltae.GetXaxis()->SetTitleOffset(1.2);
     427  deltae.GetXaxis()->SetTitle("(E_{EST} - E_{MC}) / E_{MC}");
     428  deltae.Draw("e");
     429 
     430  TFile out2("energyest.root", "update");
     431  ((TH2*)mh3ed.GetHist())->Write("mh3ed");
     432  ((TH2*)mh3ed2.GetHist())->Write("mh3ed2");
     433  ((TH2*)mhe.GetHist())->Write("mhe");
     434  deltae.Write();
     435  out2.Close();
     436
     437  return;
     438
    403439}
Note: See TracChangeset for help on using the changeset viewer.