Changeset 2255 for trunk/MagicSoft/Mars/macros

Timestamp:

07/01/03 10:46:07 (21 years ago)

Author:

wittek

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/macros

Files:

• CT1Analysis.C (modified) (22 diffs)
• ONOFFCT1Analysis.C (modified) (6 diffs)
• fluxunfold.C (added)
• unfold.C (added)

trunk/MagicSoft/Mars/macros/CT1Analysis.C

@@ -232 +232 @@
     //  - write root file for ON data after final cuts (ON3.root))
 
-    Bool_t JobF_XX  = kFALSE;  
+    Bool_t JobF_XX  = kTRUE;  
     Bool_t WXX      = kFALSE;  // write out root file ON3.root ?
 
@@ -2175 +2175 @@
 
     // type of data to be analysed
-    //TString typeData = "ON";
+    TString typeData = "ON";
     //TString typeData = "OFF";
-    TString typeData = "MC";
+    //TString typeData = "MC";
     gLog << "typeData = " << typeData << endl;
 
@@ -2188 +2188 @@
 
     //TString XX("NN");
-    //TString XX("SC");
-    TString XX("RF");
+    TString XX("SC");
+    //TString XX("RF");
     TString fhadronnessName("Had");
     fhadronnessName += XX;
@@ -2195 +2195 @@
 
     // maximum values of the hadronness, |ALPHA| and DIST
-    Float_t maxhadronness   = 0.20;
+    Float_t maxhadronness   = 0.30;
     Float_t maxalpha        = 20.0;
     Float_t maxdist         = 10.0;
@@ -2295 +2295 @@
     tliston.AddToList(&contfinalgh);
 
-    tliston.AddToList(&contfinal);
-
     tliston.AddToList(&fillhadnn);
     tliston.AddToList(&fillhadsc);
@@ -2307 +2305 @@
     tliston.AddToList(&hfill5);
 
-
+    tliston.AddToList(&contfinal);
 
     //*****************************
@@ -2340 +2338 @@
     pliston.FindObject("MHNewImagePar")->DrawClone();
     pliston.FindObject("MHStarMap")->DrawClone();
+
+
+    //-------------------------------------------
+    // fit alpha distribution to get the number of excess events
+    //
+
+    MHHillasSrc* hillasSrc = 
+      (MHHillasSrc*)(pliston.FindObject("MHHillasSrc"));
+    TH1F* alphaHist = (TH1F*)(hillasSrc->GetHistAlpha());
+  
+    MHOnSubtraction onsub;
+    onsub.Calc(alphaHist, &pliston, kTRUE, 13.1);
+    //-------------------------------------------
+
 
     DeleteBinnings(&pliston);
@@ -2713 +2725 @@
 
     TString typeMC   = "MC";
-    TString ext      = "2.root";
-    TString extout   = "3.root";
+    TString ext      = "3.root";
+    TString extout   = "4.root";
 
     //------------------------------
@@ -2721 +2733 @@
 
     //TString XX("NN");
-    TString XX("SC");
-    //TString XX("RF");
+    //TString XX("SC");
+    TString XX("RF");
     TString fhadronnessName("Had");
     fhadronnessName += XX;
@@ -2735 +2747 @@
          << maxdist << endl;
 
+    //------------------------------
+    // name of MC file to be used for optimizing the energy estimator
+    TString filenameOpt(outPath);
+    filenameOpt += typeMC;
+    filenameOpt += ext; 
+    gLog << "filenameOpt = " << filenameOpt << endl;
+
+    //------------------------------
+    // name of file containing the parameters of the energy estimator
+    TString energyParName(outPath);
+    energyParName += "energyest_";
+    //energyParName += XX;
+    energyParName += ".root";
+    gLog << "energyParName = " << energyParName << endl;
 
     //------------------------------
@@ -2740 +2766 @@
     TString filenameArea(outPath);
     filenameArea += typeMC;
-    filenameArea += "_";
-    filenameArea += XX;
+    //filenameArea += "_";
+    //filenameArea += XX;
     filenameArea += ext; 
     gLog << "filenameArea = " << filenameArea << endl;
@@ -2749 +2775 @@
     TString collareaName(outPath);
     collareaName += "area_";
-    collareaName += XX;
+    //collareaName += XX;
     collareaName += ".root";
     gLog << "collareaName = " << collareaName << endl;
@@ -2757 +2783 @@
     TString filenameData(outPath);
     filenameData += typeData;
-    filenameData += "_";
-    filenameData += XX;
+    //filenameData += "_";
+    //filenameData += XX;
     filenameData += ext;
     gLog << "filenameData = " << filenameData << endl;
@@ -2766 +2792 @@
     TString filenameDataout(outPath);
     filenameDataout += typeData;
-    filenameDataout += "_";
-    filenameDataout += XX;
+    //filenameDataout += "_";
+    //filenameDataout += XX;
     filenameDataout += extout;
     gLog << "filenameDataout = " << filenameDataout << endl;
@@ -2790 +2816 @@
     MContinue conthadrons(&cuthadrons);
 
-    MHMcCT1CollectionArea* collarea = new MHMcCT1CollectionArea();
+    //MHMcCT1CollectionArea* collarea = new MHMcCT1CollectionArea();
+    //MHMcCT1CollectionArea* collarea;
 
     MFillH filler("MHMcCT1CollectionArea", "MMcEvt");
@@ -2800 +2827 @@
     parlist.AddToList(&tasklist);
     InitBinnings(&parlist);
-    parlist.AddToList(collarea);
+    //parlist.AddToList(collarea);
 
     //********************************
@@ -2827 +2854 @@
     // and display the histograms
     //
-    //MHMcCT1CollectionArea *collarea = parlist.FindObject("MHMcCT1CollectionArea");
+    MHMcCT1CollectionArea *collarea = 
+         (MHMcCT1CollectionArea*)parlist.FindObject("MHMcCT1CollectionArea");
     collarea->CalcEfficiency();
     collarea->DrawClone("lego");
@@ -2847 +2875 @@
 
 
-    //*************************************************************************
-    //
-    // Analyse the data
-    //
-
-    MTaskList tliston;
-    MParList pliston;
-
-    // geometry is needed in  MHHillas... classes 
-    MGeomCam *fGeom = 
-             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
-
-
     TString fHilName    = "MHillas"; 
     TString fHilNameExt = "MHillasExt"; 
     TString fHilNameSrc = "MHillasSrc"; 
     TString fImgParName = "MNewImagePar"; 
+
+
+    //===========================================================
+    //
+    // Optimization of energy estimator
+    //
+
+    TString inpath("");
+    TString outpath("");
+    Int_t howMany = 2000;
+    CT1EEst(inpath,   filenameOpt,   outpath, energyParName, 
+            fHilName, fHilNameSrc,   fhadronnessName,
+            howMany,  maxhadronness, maxalpha, maxdist);
+
+    //-----------------------------------------------------------
+    //
+    // Read in parameters of energy estimator ("MMcEnergyEst")
+    //                   and migration matrix ("MHMcEnergyMigration")
+    //
+    gLog << "================================================================"
+         << endl;
+    gLog << "Macro CT1Analysis.C : read parameters of energy estimator and migration matrix from file '"
+         << energyParName << "'" << endl;
+    TFile enparam(energyParName);
+    enparam.ls();
+
+    //MMcEnergyEst mcest("MMcEnergyEst"); 
+    //mcest.Read("MMcEnergyEst");
+    MMcEnergyEst &mcest = *((MMcEnergyEst*)gROOT->FindObject("MMcEnergyEst"));
+
+    gLog << "Parameters of energy estimator were read in" << endl;
+
+    TArrayD parA(5);
+    TArrayD parB(7);
+    for (Int_t i=0; i<parA.GetSize(); i++)
+      parA[i] = mcest.GetCoeff(i);
+    for (Int_t i=0; i<parB.GetSize(); i++)
+      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;
+
+    //*************************************************************************
+    //
+    // Analyse the data
+    //
+    gLog << "============================================================"
+         << endl;
+    gLog << "Analyse the data" << endl;
+
+    MTaskList tliston;
+    MParList pliston;
+
+    // geometry is needed in  MHHillas... classes 
+    MGeomCam *fGeom = 
+             (MGeomCam*)pliston->FindCreateObj("MGeomCamCT1", "MGeomCam");
+
 
     //-------------------------------------------
@@ -2914 +2992 @@
     fillhadrf.SetName("HhadRF");
 
+    //---------------------------
+    // calculate estimated energy
+
+    MEnergyEstParam eeston(fHilName);
+    eeston.Add(fHilNameSrc);
+
+    eeston.SetCoeffA(parA);
+    eeston.SetCoeffB(parB);
+    //---------------------------
+
 
     MFillH hfill1("MHHillas",    fHilName);
@@ -2955 +3043 @@
       tliston.AddToList(&write);
 
+    tliston.AddToList(&eeston);
+
     tliston.AddToList(&fillhadnn);
     tliston.AddToList(&fillhadsc);
@@ -2977 +3067 @@
     evtloop.SetProgressBar(&bar);
 
-    Int_t maxevents = -1;
-    //Int_t maxevents = 10000;
+    //Int_t maxevents = -1;
+    Int_t maxevents = 10000;
     if ( !evtloop.Eventloop(maxevents) )
         return;
@@ -3010 +3100 @@
 
 
+
+
+
+
+
+
+
+
+
+
+

trunk/MagicSoft/Mars/macros/ONOFFCT1Analysis.C

@@ -75 +75 @@
         gLog << "DeleteBinnings" << endl;
 
+        if (!plist)
+        {
+          gLog << "Deletebinnins : MParlist no longer existing" << endl;
+          return;
+        }
+
         TObject *bin;
 
@@ -113 +119 @@
         bin = plist->FindObject("BinningDiffsigma2");
         if (bin) delete bin;
+
+        gLog << "exit DeleteBinnings" << endl;
 }
 
@@ -134 +142 @@
 
       // path for input for Mars
-      TString inPath = "~magican/ct1test/wittek/marsoutput/optionA/";
+      TString inPath = "~magican/ct1test/wittek/marsoutput/optionC/";
 
       // path for output from Mars
-      TString outPath = "~magican/ct1test/wittek/marsoutput/optionA/";
+      TString outPath = "~magican/ct1test/wittek/marsoutput/optionC/";
 
       //-----------------------------------------------
@@ -153 +161 @@
     //  - write root file for ON (or OFF or MC) data (ON1.root, ...);
 
-    Bool_t JobA    = kFALSE;  
-    Bool_t WPad    = kFALSE;   // write out padding histograms ?
+    Bool_t JobA    = kTRUE;  
+    Bool_t WPad    = kTRUE;   // write out padding histograms ?
     Bool_t RPad    = kFALSE;   // read in padding histograms ?
     Bool_t Wout    = kFALSE;   // write out root file ON1.root 
@@ -203 +211 @@
     //  - make plots
 
-    Bool_t JobD  = kTRUE;  
+    Bool_t JobD  = kFALSE;  
 
 
@@ -2245 +2253 @@
 
     MHHillasSrc* hillasSrc = 
-      (MHHillasSrc*)(pliston->FindObject("MHHillasSrc"));
+      (MHHillasSrc*)(pliston.FindObject("MHHillasSrc"));
     TH1F* alphaHist = (TH1F*)(hillasSrc->GetHistAlpha());
   
     MHOnSubtraction onsub;
-    onsub.Calc(alphaHist, &pliston, kTRUE, 25);
+    onsub.Calc(alphaHist, &pliston, kTRUE, 13.1);
     //-------------------------------------------