Changeset 2436 for trunk/MagicSoft/Mars

Timestamp:

10/27/03 08:58:48 (21 years ago)

Author:

wittek

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/macros/CT1Analysis.C (modified) (52 diffs)

trunk/MagicSoft/Mars/macros/CT1Analysis.C

@@ -210 +210 @@
       //-----------------------------------------------
 
-      TEnv env("macros/CT1env.rc");
-      Bool_t printEnv = kFALSE;
+      //TEnv env("macros/CT1env.rc");
+      //Bool_t printEnv = kFALSE;
 
     //************************************************************************
@@ -219 +219 @@
     //  - write root file for ON data (ON1.root);
 
-    Bool_t JobA_ON = kTRUE;  
-    Bool_t WHistON = kTRUE;   // write out histogram sigmabar vs. Theta ?
-    Bool_t WON1    = kTRUE;   // write out root file ON1.root ?
+    Bool_t JobA_ON = kFALSE;  
+    Bool_t WHistON = kFALSE;   // write out histogram sigmabar vs. Theta ?
+    Bool_t WON1    = kFALSE;   // write out root file ON1.root ?
 
 
@@ -234 +234 @@
 
     // Job B_RF_UP : read ON1.root (or MC1.root) file 
-    //  - depending on RLook : create (or read in) hadron and gamma matrix
-    //  - depending on RTree : create (or read in) trees
+    //  - if RTrainRF = TRUE : read in training matrices for hadrons and gammas
+    //  - if CTrainRF = TRUE : create  training matrices for hadrons and gammas
+    //  - if RTree    = TRUE : read in trees, otherwise create trees
     //  - calculate hadroness for method of RANDOM FOREST
-    //  - update the input files with the hadroness (ON1.root or MC1.root)
+    //  - update the input files with the hadroness (==>ON2.root or MC2.root)
 
     Bool_t JobB_RF_UP  = kFALSE;  
-    Bool_t RLookRF     = kFALSE;  // read in look-alike events
+    Bool_t RTrainRF    = kFALSE;  // read in matrices of training events
+    Bool_t CTrainRF    = kFALSE;  // create matrices of training events
     Bool_t RTree       = kFALSE;  // read in trees
     Bool_t WRF         = kFALSE;  // update input root file ?
@@ -247 +249 @@
 
 
-    // Job B_SC_UP : read ON1.root (or MC1.root) file 
+    // Job B_SC_UP : read ON2.root (or MC2.root) file 
     //  - depending on WParSC : create (or read in) supercuts parameter values
     //  - calculate hadroness for the SUPERCUTS
-    //  - update the input files with the hadroness (ON1.root or MC1.root)
+    //  - update the input files with the hadroness (==>ON3.root or MC3.root)
 
     Bool_t JobB_SC_UP  = kFALSE;
@@ -263 +265 @@
 
     // Job C: 
-    //  - read ON1.root and MC1.root files
+    //  - read ON3.root and MC3.root files
     //    which should have been updated to contain the hadronnesses  
     //    for the method of 
-    //              NEAREST NEIGHBORS  
+    //              RF
     //              SUPERCUTS and
-    //              RF
     //  - produce Neyman-Pearson plots
 
@@ -276 +277 @@
     // Job D :  
     //  - select g/h separation method XX
-    //  - read ON2 (or MC2) root file
+    //  - read ON3 (or MC3) root file
     //  - apply cuts in hadronness
     //  - make plots
@@ -296 +297 @@
 
     Bool_t JobE_XX  = kFALSE;  
-    Bool_t WEX      = kFALSE;  // write out root file  ?
+    Bool_t OEEst    = kFALSE;  // optimize energy estimator
+    Bool_t WEX      = kFALSE;  // update root file  ?
     Bool_t WRobert  = kFALSE;  // write out Robert's file  ?
 
@@ -492 +494 @@
     MEvtLoop evtloop;
     evtloop.SetParList(&pliston);
-    evtloop.ReadEnv(env, "", printEnv);
+    //evtloop.ReadEnv(env, "", printEnv);
     evtloop.SetProgressBar(&bar);
     //if (WON1)
@@ -822 +824 @@
     MEvtLoop evtloop;
     evtloop.SetParList(&plist);
-    evtloop.ReadEnv(env, "", printEnv);
+    //evtloop.ReadEnv(env, "", printEnv);
     evtloop.SetProgressBar(&bar);
     //if (WMC1)    
@@ -863 +865 @@
 
 
-    //  - create (or read in) the matrices of look-alike events for gammas 
+    //  - create (or read in) the matrices of training events for gammas 
     //    and hadrons
     //  - create (or read in) the trees
@@ -877 +879 @@
 
     gLog << "" << endl;
-    gLog << "Macro CT1Analysis : JobB_RF_UP, RLookRF, RTree, WRF = " 
+    gLog << "Macro CT1Analysis : JobB_RF_UP, RTrainRF, CTrainRF, RTree, WRF = "
          << (JobB_RF_UP ? "kTRUE" : "kFALSE")  << ",  " 
-         << (RLookRF    ? "kTRUE" : "kFALSE")  << ",  "
+         << (RTrainRF    ? "kTRUE" : "kFALSE")  << ",  "
+         << (CTrainRF    ? "kTRUE" : "kFALSE")  << ",  "
          << (RTree      ? "kTRUE" : "kFALSE")  << ",  "
          << (WRF        ? "kTRUE" : "kFALSE")  << endl;
@@ -890 +893 @@
     Int_t NdSize   =   1;
 
+    // cut in RF hadronness
+    TString hadRFName = "HadRF";
+    Float_t maxhadronness =  0.23;
+    Float_t maxalpha      =  20.0;
+    Float_t maxdist       =  10.0;
+
+    TString fHilName    = "MHillas"; 
+    TString fHilNameExt = "MHillasExt"; 
+    TString fHilNameSrc = "MHillasSrc"; 
+    TString fImgParName = "MNewImagePar"; 
+
 
     //--------------------------------------------
     // file to be updated (either ON or MC)
 
-    TString typeInput = "ON";
-    //TString typeInput = "MC";
+    //TString typeInput = "ON";
+    TString typeInput = "MC";
     gLog << "typeInput = " << typeInput << endl;
 
@@ -913 +927 @@
 
     //--------------------------------------------
-    // files to be read for generating the look-alike events
+    // files to be read for generating the matrices of training events
     // "hadrons" :
     TString filenameON = outPath;
@@ -932 +946 @@
     
     //--------------------------------------------
-    // files of look-alike events 
+    // files of training events 
 
     TString outNameGammas = outPath;
@@ -951 +965 @@
     matrixg.EnableGraphicalOutput();
 
+    matrixg.AddColumn("cos(MMcEvt.fTelescopeTheta)");
+    matrixg.AddColumn("MSigmabar.fSigmabar");
+    matrixg.AddColumn("log10(MHillas.fSize)");
+    matrixg.AddColumn("MHillasSrc.fDist");
+    matrixg.AddColumn("MHillas.fWidth");
+    matrixg.AddColumn("MHillas.fLength");
+    matrixg.AddColumn("log10(MHillas.fSize/(MHillas.fWidth*MHillas.fLength))");
+    matrixg.AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillasExt.fM3Long)");
+    matrixg.AddColumn("MNewImagePar.fConc");
+    matrixg.AddColumn("MNewImagePar.fLeakage1");
+
     MHMatrix matrixh("MatrixHadrons");
     matrixh.EnableGraphicalOutput();
 
+    matrixh.AddColumns(matrixg.GetColumns());
+
     //--------------------------------------------
     // file of trees of the random forest 
@@ -962 +989 @@
 
    //*************************************************************************
-   // read in matrices of look-alike events
-if (RLookRF)
+   // read in matrices of training events
+if (RTrainRF)
   {
     const char* mtxName = "MatrixGammas";
@@ -1005 +1032 @@
 
    //*************************************************************************
-   // create matrices of look-alike events
-if (!RLookRF)
+   // create matrices of training events
+if (CTrainRF)
   {
     gLog << "" << endl;
     gLog << "========================================================" << endl;
-    gLog << " Create matrices of look-alike events" << endl;
+    gLog << " Create matrices of training events" << endl;
     gLog << " Gammas :" << endl;
 
@@ -1085 +1112 @@
     MEvtLoop evtloopg;
     evtloopg.SetParList(&plistg);
-    evtloopg.ReadEnv(env, "", printEnv);
+    //evtloopg.ReadEnv(env, "", printEnv);
     evtloopg.SetProgressBar(&matrixbar);
 
@@ -1141 +1168 @@
     MEvtLoop evtlooph;
     evtlooph.SetParList(&plisth);
-    evtlooph.ReadEnv(env, "", printEnv);
+    //evtlooph.ReadEnv(env, "", printEnv);
     evtlooph.SetProgressBar(&matrixbar);
 
@@ -1155 +1182 @@
 
 
-    // write out look-alike events
+    // write out training events
 
     gLog << "" << endl;
     gLog << "========================================================" << endl;
-    gLog << "Write out look=alike events" << endl;
+    gLog << "Write out training events" << endl;
 
 
@@ -1171 +1198 @@
 
       gLog << "" << endl;
-      gLog << "Macro CT1Analysis : matrix of look-alike events for gammas written onto file "
+      gLog << "Macro CT1Analysis : matrix of training events for gammas written onto file "
            << outNameGammas << endl;
 
@@ -1183 +1210 @@
 
       gLog << "" << endl;
-      gLog << "Macro CT1Analysis : matrix of look-alike events for hadrons written onto file "
+      gLog << "Macro CT1Analysis : matrix of training events for hadrons written onto file "
            << outNameHadrons << endl;
 
   }
-   //**********   end of creating matrices of look-alike events   ***********
+   //**********   end of creating matrices of training events   ***********
 
 
@@ -1301 +1328 @@
 
 
-
     //-----------------------------------------------------------------
     // Update the input files with the RF hadronness
     //
-
  if (WRF)
   {
@@ -1328 +1353 @@
     read.DisableAutoScheme();
 
-    TString fHilName    = "MHillas"; 
-    TString fHilNameExt = "MHillasExt"; 
-    TString fHilNameSrc = "MHillasSrc"; 
-    TString fImgParName = "MNewImagePar"; 
 
     //.......................................................................
     // calculate hadronnes for method of RANDOM FOREST
 
-    TString hadRFName = "HadRF";
+
     MRanForestCalc rfcalc;
     rfcalc.SetHadronnessName(hadRFName);
@@ -1360 +1381 @@
 
     //-----------------------------------------------------------------
-    // geometry is needed in  MHHillas... classes 
-    MGeomCam *fGeom = 
-             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
-
-
-    Float_t maxhadronness =  0.40;
-    Float_t maxalpha      =  20.0;
-    Float_t maxdist       =  10.0;
+
 
     MFCT1SelFinal selfinalgh(fHilNameSrc);
@@ -1482 +1496 @@
     TH1  &alphaHist = absalpha->GetHist();
     alphaHist.SetXTitle("|alpha|  [\\circ]");
+    alphaHist.SetName("alpha-macro");
 
     Double_t alphasig = 13.1;
     Double_t alphamin = 30.0;
     Double_t alphamax = 90.0;
-    Int_t    degree   =    4;
+    Int_t    degree   =    2;
     Double_t significance = -99.0;
     Bool_t   drawpoly  = kTRUE;
@@ -1563 +1578 @@
 
     TString parSCfile = outPath;
-    parSCfile += "parSC_2209a";
+    parSCfile += "parSC_2310a";
 
     gLog << "parSCfile = " << parSCfile << endl;
@@ -1570 +1585 @@
     // file to be updated (either ON or MC)
 
-    TString typeInput = "ON";
-    //TString typeInput = "MC";
+    //TString typeInput = "ON";
+    TString typeInput = "MC";
     gLog << "typeInput = " << typeInput << endl;
 
@@ -1612 +1627 @@
 
     //--------------------------------------------
-    // files to contain the matrices (genertated from filenameTrain and
-    //                                                filenameTest)
+    // files to contain the matrices (generated from filenameTrain and
+    //                                               filenameTest)
     // 
     // for the training
@@ -1642 +1657 @@
     //--------------------------
     // create matrices and write them onto files 
-    if (!RMatrix  &&  (WOptimize || RTest)  )
+    if (!RMatrix)
     {
       MH3 &mh3 = *(new MH3("MHillas.fSize"));
@@ -1911 +1926 @@
       //MWriteRootFile write = new MWriteRootFile(outNameImage, "RECREATE");
 
-    /*
+    
       MWriteRootFile write(outNameImage, "RECREATE");
 
@@ -1927 +1942 @@
       write.AddContainer("HadRF",         "Events");
       write.AddContainer(hadSCName,       "Events");
-    */
+    
 
     //-----------------------------------------------------------------
@@ -2001 +2016 @@
     tliston.AddToList(&fillhadsc);
 
-    //tliston.AddToList(&write);
+    tliston.AddToList(&write);
     tliston.AddToList(&contfinalgh);
 
@@ -2292 +2307 @@
     // and definition of final selections
 
-    //TString XX("NN");
-    TString XX("SC");
-    //TString XX("RF");
+    //TString XX("SC");
+    TString XX("RF");
     TString fhadronnessName("Had");
     fhadronnessName += XX;
@@ -2300 +2314 @@
 
     // maximum values of the hadronness, |ALPHA| and DIST
-    Float_t maxhadronness   = 0.30;
+    Float_t maxhadronness   = 0.233;
     Float_t maxalpha        = 20.0;
     Float_t maxdist         = 10.0;
@@ -2355 +2369 @@
     contfinalgh.SetName("ContSelFinalgh");
 
-    MFillH fillhadnn("hadNN[MHHadronness]", "HadNN");
-    fillhadnn.SetName("HhadNN");
     MFillH fillhadsc("hadSC[MHHadronness]", "HadSC");
     fillhadsc.SetName("HhadSC");
@@ -2362 +2374 @@
     fillhadrf.SetName("HhadRF");
 
+    TString mh3name = "abs(Alpha)";
+    MBinning binsalphaabs("Binning"+mh3name);
+    binsalphaabs.SetEdges(50, -2.0, 98.0);
+
+    MH3 alphaabs("abs(MHillasSrc.fAlpha)");
+    alphaabs.SetName(mh3name);
+
+    TH1  &alphahist = alphaabs->GetHist();
+
+    MFillH alpha(&alphaabs);
+    alpha.SetName("FillAlphaAbs");
 
     MFillH hfill1("MHHillas",    fHilName);
@@ -2391 +2414 @@
     pliston.AddToList(&tliston);
     InitBinnings(&pliston);
-
+    pliston.AddToList(&binsalphaabs);
+    pliston.AddToList(&alphaabs);
 
     //*****************************
@@ -2400 +2424 @@
     tliston.AddToList(&contfinalgh);
 
-    tliston.AddToList(&fillhadnn);
     tliston.AddToList(&fillhadsc);
     tliston.AddToList(&fillhadrf);
 
+    tliston.AddToList(&alpha);
     tliston.AddToList(&hfill1);
     tliston.AddToList(&hfill2);
@@ -2434 +2458 @@
     //
 
-    pliston.FindObject("hadNN", "MHHadronness")->DrawClone();
     pliston.FindObject("hadRF", "MHHadronness")->DrawClone();
     pliston.FindObject("hadSC", "MHHadronness")->DrawClone();
@@ -2446 +2469 @@
 
     //-------------------------------------------
-    // fit alpha distribution to get the number of excess events
-    //
-
-    MHHillasSrc* hillasSrc = 
-      (MHHillasSrc*)(pliston.FindObject("MHHillasSrc"));
-    TH1F* alphaHist = (TH1F*)(hillasSrc->GetHistAlpha());
+
+    // fit alpha distribution to get the number of excess events and
+    // calculate significance of gamma signal in the alpha plot
   
-    MHOnSubtraction onsub;
-    onsub.Calc(alphaHist, kTRUE, 13.1);
+    MH3* absalpha = (MH3*)(pliston.FindObject(mh3name, "MH3"));
+    TH1  &alphaHist = absalpha->GetHist();
+    alphaHist.SetXTitle("|alpha|  [\\circ]");
+    alphaHist.SetName("alpha-JobD");
+
+    Double_t alphasig = 13.1;
+    Double_t alphamin = 30.0;
+    Double_t alphamax = 90.0;
+    Int_t    degree   =    2;
+    Double_t significance = -99.0;
+    Bool_t   drawpoly  = kTRUE;
+    Bool_t   fitgauss  = kTRUE;
+    Bool_t   print     = kTRUE;
+
+    MHFindSignificance findsig;
+    findsig.SetRebin(kTRUE);
+    findsig.SetReduceDegree(kFALSE);
+
+    findsig.FindSigma(&alphaHist, alphamin, alphamax, degree, 
+                        alphasig, drawpoly, fitgauss, print);
+    significance = findsig.GetSignificance();
+    Float_t alphasi = findsig.GetAlphasi();
+
+    gLog << "For file '" << filenameData << "' : " << endl;
+    gLog << "Significance of gamma signal after supercuts : "
+         << significance << " (for |alpha| < " << alphasi << " degrees)" 
+         << endl;
+
+    findsig.SigmaVsAlpha(&alphaHist, alphamin, alphamax, degree, print);
+
     //-------------------------------------------
 
@@ -2488 +2536 @@
 
     gLog << "" << endl;
-    gLog << "Macro CT1Analysis : JobE_XX, WEX = " 
+    gLog << "Macro CT1Analysis : JobE_XX, OEEst, WEX = " 
          << (JobE_XX ? "kTRUE" : "kFALSE")  << ",  " 
+         << (OEEst ?   "kTRUE" : "kFALSE")  << ",  " 
          << (WEX     ? "kTRUE" : "kFALSE")  << endl;
 
@@ -2507 +2556 @@
     // and definition of final selections
 
-    //TString XX("NN");
     TString XX("SC");
     //TString XX("RF");
@@ -2515 +2563 @@
 
     // maximum values of the hadronness, |ALPHA| and DIST
-    Float_t maxhadronness   = 0.40;
+    Float_t maxhadronness   = 0.4;
     Float_t maxalpha        = 20.0;
     Float_t maxdist         = 10.0;
@@ -2676 +2724 @@
     TString fImgParName = "MNewImagePar"; 
 
-  
-    //===========================================================
+
+ if (OEEst)
+ { 
+   //===========================================================
     //
     // Optimization of energy estimator
@@ -2691 +2741 @@
             binsE, binsTheta);
     gLog << "Macro CT1Analysis.C : returning from CT1EEst" << endl;
-
-
-  if (WEX)
-  {
+ }
+
+ if (WEX)
+ {
     //-----------------------------------------------------------
     //
@@ -2706 +2756 @@
     TFile enparam(energyParName);
     enparam.ls();
-
-    //MMcEnergyEst mcest("MMcEnergyEst"); 
-    //mcest.Read("MMcEnergyEst");
-    MMcEnergyEst &mcest = *((MMcEnergyEst*)gROOT->FindObject("MMcEnergyEst"));
-
+    MMcEnergyEst mcest("MMcEnergyEst"); 
+    mcest.Read("MMcEnergyEst");
+
+    //MMcEnergyEst &mcest = *((MMcEnergyEst*)gROOT->FindObject("MMcEnergyEst"));
     gLog << "Parameters of energy estimator were read in" << endl;
+
+
+    gLog << "Read in Migration matrix" << endl;   
+
+    MHMcEnergyMigration mighiston("MHMcEnergyMigration");
+    mighiston.Read("MHMcEnergyMigration");
+    //MHMcEnergyMigration &mighiston = 
+    //      *((MHMcEnergyMigration*)gROOT->FindObject("MHMcEnergyMigration"));
+
+    gLog << "Migration matrix was read in" << endl;
+
 
     TArrayD parA(mcest.GetNumCoeffA());
@@ -2720 +2780 @@
       parB[i] = mcest.GetCoeff( i+parA.GetSize() );
 
-    gLog << "Read in Migration matrix" << endl;   
-
-    //MHMcEnergyMigration mighiston("MHMcEnergyMigration");
-    //mighiston.Read("MHMcEnergyMigration");
-    MHMcEnergyMigration &mighiston = 
-          *((MHMcEnergyMigration*)gROOT->FindObject("MHMcEnergyMigration"));
-
-    gLog << "Migration matrix was read in" << endl;
-
     //*************************************************************************
     //
@@ -2758 +2809 @@
    
       //MWriteRootFile &write = *(new MWriteRootFile(filenameDataout));
+
 
       MWriteRootFile write(filenameDataout, "RECREATE");
@@ -2916 +2968 @@
 
     Int_t maxevents = -1;
-    //Int_t maxevents = 10000;
     if ( !evtloop.Eventloop(maxevents) )
         return;
@@ -2952 +3003 @@
     gLog << "before DeleteBinnings" << endl;
 
-    //DeleteBinnings(&pliston);
+    DeleteBinnings(&pliston);
 
     gLog << "before Robert's code" << endl;
@@ -3036 +3087 @@
   }
 
-    gLog << "Macro CT1Analysis : End of Job F_XX" << endl;
+    gLog << "Macro CT1Analysis : End of Job E_XX" << endl;
     gLog << "=======================================================" << endl;
  }