Changeset 2357 for trunk/MagicSoft/Mars/manalysis

Timestamp:

09/24/03 13:29:00 (21 years ago)

Author:

wittek

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/manalysis

Files:

• MCT1FindSupercuts.cc (modified) (29 diffs)
• MCT1FindSupercuts.h (modified) (1 diff)
• MCT1Supercuts.cc (modified) (8 diffs)
• MCT1Supercuts.h (modified) (5 diffs)
• MCT1SupercutsCalc.cc (modified) (10 diffs)
• MCT1SupercutsCalc.h (modified) (4 diffs)
• MMinuitInterface.cc (modified) (5 diffs)

trunk/MagicSoft/Mars/manalysis/MCT1FindSupercuts.cc

@@ -57 +57 @@
 #include "MFEventSelector2.h"
 #include "MFillH.h"
-#include "MGeomCamCT1Daniel.h"
+//#include "MGeomCamCT1Daniel.h"
+#include "MGeomCamCT1.h"
 #include "MH3.h"
 #include "MHCT1Supercuts.h"
@@ -95 +96 @@
                          Double_t *par, Int_t iflag)
 {
+    //cout <<  "entry fcnSupercuts" << endl;
+
     //-------------------------------------------------------------
     // save pointer to the MINUIT object for optimizing the supercuts
@@ -106 +109 @@
 
     MParList *plistfcn   = (MParList*)evtloopfcn->GetParList();
+    //MTaskList *tasklistfcn   = (MTaskList*)plistfcn->FindObject("MTaskList");
 
     MCT1Supercuts *super = (MCT1Supercuts*)plistfcn->FindObject("MCT1Supercuts");
@@ -130 +134 @@
     // for testing
     //TArrayD checkparameters = super->GetParameters();
+    //gLog << "fcnsupercuts : fNpari, fNparx =" << fNpari << ",  " 
+    //     << fNparx  << endl;
     //gLog << "fcnsupercuts : i, par, checkparameters =" << endl;
     //for (Int_t i=0; i<fNparx; i++)
@@ -142 +148 @@
     evtloopfcn->Eventloop();
 
+    //tasklistfcn->PrintStatistics(0, kTRUE);
+
     MH3* alpha = (MH3*)plistfcn->FindObject("AlphaFcn", "MH3");
     if (!alpha)
@@ -147 +155 @@
 
     TH1 &alphaHist = alpha->GetHist();
-    //alphaHist.SetXTitle("|\\alpha|  [\\circ]");
+    alphaHist.SetName("alpha-fcnSupercuts");
 
     //-------------------------------------------
@@ -161 +169 @@
     const Double_t alphamin = 30.0;
     const Double_t alphamax = 90.0;
-    const Int_t    degree   =    4;
+    const Int_t    degree   =    2;
 
     Bool_t drawpoly;
@@ -240 +248 @@
     //---------------------------
     // camera geometry is needed for conversion mm ==> degree
-    fCam = new MGeomCamCT1Daniel; 
-
-    // matrices to contain the the training/test samples
+    //fCam = new MGeomCamCT1Daniel; 
+    fCam = new MGeomCamCT1; 
+
+    // matrices to contain the training/test samples
     fMatrixTrain = new MHMatrix("MatrixTrain");
     fMatrixTest  = new MHMatrix("MatrixTest");
@@ -277 +286 @@
 //
 Bool_t MCT1FindSupercuts::DefineTrainMatrix(const TString &nametrain,
-                          const Int_t howmanytrain, MH3 &hreftrain)
+                          const Int_t howmanytrain, MH3 &hreftrain,
+                          const TString &filetrain)
 {
     if (nametrain.IsNull() || howmanytrain <= 0)
@@ -343 +353 @@
     *fLog << "=============================================" << endl;
 
+    //--------------------------
+    // write out training matrix
+
+    if (filetrain != "")
+    {
+      TFile filetr(filetrain, "RECREATE");
+      fMatrixTrain->Write();
+      filetr.Close();
+
+      *fLog << "MCT1FindSupercuts::DefineTrainMatrix; Training matrix was written onto file '"
+            << filetrain << "'" << endl;
+    }
+
+
     return kTRUE;
 }
@@ -356 +380 @@
 //
 Bool_t MCT1FindSupercuts::DefineTestMatrix(const TString &nametest,
-                          const Int_t howmanytest, MH3 &hreftest)
+                          const Int_t howmanytest, MH3 &hreftest,
+                          const TString &filetest)
 {
     if (nametest.IsNull() || howmanytest<=0)
@@ -422 +447 @@
     *fLog << "=============================================" << endl;
 
+    //--------------------------
+    // write out test matrix
+
+    if (filetest != "")
+    {
+      TFile filete(filetest, "RECREATE", "");
+      fMatrixTest->Write();
+      filete.Close();
+
+      *fLog << "MCT1FindSupercuts::DefineTestMatrix; Test matrix was written onto file '"
+            << filetest << "'" << endl;
+    }
+
+
   return kTRUE;
 }
@@ -434 +473 @@
                           const TString &name,
                           const Int_t howmanytrain, MH3 &hreftrain,
-                          const Int_t howmanytest,  MH3 &hreftest)
+                          const Int_t howmanytest,  MH3 &hreftest,
+                          const TString &filetrain, const TString &filetest)
 {
     *fLog << "=============================================" << endl;
@@ -511 +551 @@
     evtloop.SetProgressBar(&bar);
 
-    if (!evtloop.Eventloop())
+    Int_t maxev = -1;
+    if (!evtloop.Eventloop(maxev))
       return kFALSE;
 
@@ -520 +561 @@
     if (TMath::Abs(generatedtrain-howmanytrain) > TMath::Sqrt(9.*howmanytrain))
     {
-      *fLog << "MCT1FindSupercuts::DefineTrainMatrix; no.of generated events ("
+      *fLog << "MCT1FindSupercuts::DefineTrainTestMatrix; no.of generated events ("
             << generatedtrain 
             << ") is incompatible with the no.of requested events ("
@@ -530 +571 @@
     if (TMath::Abs(generatedtest-howmanytest) > TMath::Sqrt(9.*howmanytest))
     {
-      *fLog << "MCT1FindSupercuts::DefineTestMatrix; no.of generated events ("
+      *fLog << "MCT1FindSupercuts::DefineTrainTestMatrix; no.of generated events ("
             << generatedtest 
             << ") is incompatible with the no.of requested events ("
@@ -541 +582 @@
 
 
+    //--------------------------
+    // write out training matrix
+
+    if (filetrain != "")
+    {
+      TFile filetr(filetrain, "RECREATE");
+      fMatrixTrain->Write();
+      filetr.Close();
+
+      *fLog << "MCT1FindSupercuts::DefineTrainTestMatrix; Training matrix was written onto file '"
+            << filetrain << "'" << endl;
+    }
+
+    //--------------------------
+    // write out test matrix
+
+    if (filetest != "")
+    {
+      TFile filete(filetest, "RECREATE", "");
+      fMatrixTest->Write();
+      filete.Close();
+
+      *fLog << "MCT1FindSupercuts::DefineTrainTestMatrix; Test matrix was written onto file '"
+            << filetest << "'" << endl;
+    }
+
   return kTRUE;
 }
@@ -548 +615 @@
 // Read training and test matrices from files
 //
-//
-// $$$$$$$$$$$$$$ this does not work !!! ??? $$$$$$$$$$$$$$$$$$$$$$
 //
 
@@ -563 +628 @@
   *fLog << "MCT1FindSupercuts::ReadMatrix; Training matrix was read in from file '"
         << filetrain << "'" << endl;
+  filetr.Close();
 
 
@@ -574 +640 @@
   *fLog << "MCT1FindSupercuts::ReadMatrix; Test matrix was read in from file '"
         << filetest << "'" << endl;
-
+  filete.Close();
 
   return kTRUE;  
 }
 
-// --------------------------------------------------------------------------
-//
-// Write training and test matrices onto files
-//
-//
-// $$$$$$$$$$$$$$ this does not work !!! ??? $$$$$$$$$$$$$$$$$$$$$$
-//
-//
-Bool_t MCT1FindSupercuts::WriteMatrix(const TString &filetrain, const TString &filetest)
-{
-  //--------------------------
-  // write out training matrix
-
-  TFile filetr(filetrain, "RECREATE");
-
-  *fLog << "nach TFile" << endl;
-
-  fMatrixTrain->Write();
-
-  *fLog << "MCT1FindSupercuts::WriteMatrix; Training matrix was written onto file '"
-        << filetrain << "'" << endl;
-
-
-  //--------------------------
-  // write out test matrix
-
-  TFile filete(filetest, "RECREATE");
-  fMatrixTest->Print("SizeCols");
-  fMatrixTest->Write();
-
-  *fLog << "MCT1FindSupercuts::WriteMatrix; Test matrix was written onto file '"
-        << filetest << "'" << endl;
-
-
-  return kTRUE;  
-}
 
 //------------------------------------------------------------------------
@@ -645 +675 @@
 //                       put the supercuts hadronness
 //
-// - for the minimization, the starting values of the parameters are taken as 
-//   MCT1Supercuts::GetParameters(fVinit)
+// - for the minimization, the starting values of the parameters are taken  
+//     - from file parSCinit (if it is != "")
+//     - or from the arrays params and/or steps 
 //
 //----------------------------------------------------------------------
-Bool_t MCT1FindSupercuts::FindParams()
+Bool_t MCT1FindSupercuts::FindParams(TString parSCinit,
+                                     TArrayD &params, TArrayD &steps)
 {
     // Setup the event loop which will be executed in the 
     //                 fcnSupercuts function  of MINUIT
     //
+    // parSCinit is the name of the file containing the initial values 
+    // of the parameters; 
+    // if parSCinit = ""   'params' and 'steps' are taken as initial values
 
     *fLog << "=============================================" << endl;
@@ -685 +720 @@
     MMatrixLoop loop(fMatrixTrain);
 
+    //--------------------------------
     // create container for the supercut parameters
-    // and set them to their default values
+    // and set them to their initial values
     MCT1Supercuts super;
 
+    // take initial values from file parSCinit
+    if (parSCinit != "")
+    {
+      TFile inparam(parSCinit);
+      super.Read("MCT1Supercuts");
+      inparam.Close();
+      *fLog << "MCT1FindSupercuts::FindParams; initial values of parameters are taken from file "
+            << parSCinit << endl;
+    }
+
+    // take initial values from 'params' and/or 'steps'
+    else if (params.GetSize() != 0  || steps.GetSize()  != 0 )
+    {
+      if (params.GetSize()  != 0)
+      {
+        *fLog << "MCT1FindSupercuts::FindParams; initial values of parameters are taken from 'params'"
+              << endl;
+        super.SetParameters(params);
+      }
+      if (steps.GetSize()  != 0)
+      {
+        *fLog << "MCT1FindSupercuts::FindParams; initial step sizes are taken from 'steps'"
+              << endl;
+        super.SetStepsizes(steps);
+      }
+    }
+    else
+    {
+        *fLog << "MCT1FindSupercuts::FindParams; initial values and step sizes are taken from the MCT1Supercits constructor"
+              << endl;
+    }
+
+
+    //--------------------------------
     // calculate supercuts hadronness
     fCalcHadTrain->SetHadronnessName(fHadronnessName);
@@ -761 +831 @@
     //
 
-    // get initial values of parameters from MCT1Supercuts
+    // get initial values of parameters 
     fVinit = super.GetParameters();
+    fStep  = super.GetStepsizes();
 
     TString name[fVinit.GetSize()];
@@ -776 +847 @@
         name[i]   = "p";
         name[i]  += i+1;
-        fStep[i]  = TMath::Abs(fVinit[i]/10.0);
+        //fStep[i]  = TMath::Abs(fVinit[i]/10.0);
         fLimlo[i] = -100.0;
         fLimup[i] =  100.0;
         fFix[i]   =     0;
-
-        // vary only first 48 parameters
-        if (i >= 48)
-        {
-          fStep[i] = 0.0;
-          fFix[i]  =   1;
-        }
-    }
+    }
+
+    // these parameters make no sense, fix them at 0.0
+    fVinit[33] = 0.0;
+    fStep[33]  = 0.0;
+    fFix[33]   = 1;
+
+    fVinit[36] = 0.0;
+    fStep[36]  = 0.0;
+    fFix[36]   = 1;
+
+    fVinit[39] = 0.0;
+    fStep[39]  = 0.0;
+    fFix[39]   = 1;
+
+    fVinit[41] = 0.0;
+    fStep[41]  = 0.0;
+    fFix[41]   = 1;
+
+    fVinit[44] = 0.0;
+    fStep[44]  = 0.0;
+    fFix[44]   = 1;
+
+    fVinit[47] = 0.0;
+    fStep[47]  = 0.0;
+    fFix[47]   = 1;
+
+    // vary only first 48 parameters
+    //for (UInt_t i=0; i<fNpar; i++)
+    //{
+    //    if (i >= 48)
+    //  {
+    //      fStep[i] = 0.0;
+    //      fFix[i]  =   1;
+    //  }
+    //}
+ 
 
     // -------------------------------------------
@@ -910 +1010 @@
     MH3 alphabefore("MatrixTest[7]");
     alphabefore.SetName(mh3NameB);
+
+    TH1 &alphahistb = alphabefore.GetHist();
+    alphahistb.SetName("alphaBefore-TestParams");
+
     MFillH fillalphabefore(&alphabefore);
     fillalphabefore.SetName("FillAlphaBefore");
@@ -924 +1028 @@
     MH3 alphaafter("MatrixTest[7]");
     alphaafter.SetName(mh3NameA);
+
+    TH1 &alphahista = alphaafter.GetHist();
+    alphahista.SetName("alphaAfter-TestParams");
+
     MFillH fillalphaafter(&alphaafter);
     fillalphaafter.SetName("FillAlphaAfter");
@@ -978 +1086 @@
     TH1  &alphaHist2 = alphaAfter->GetHist();
     alphaHist2.SetXTitle("|\\alpha|  [\\circ]");
+    alphaHist2.SetName("alpha-TestParams)");
 
     TCanvas *c = new TCanvas("AlphaAfterSC", "AlphaTest", 600, 300);
@@ -999 +1108 @@
     const Double_t alphamin = 30.0;
     const Double_t alphamax = 90.0;
-    const Int_t    degree   =    4;
+    const Int_t    degree   =    2;
     const Bool_t   drawpoly  = kTRUE;
     const Bool_t   fitgauss  = kFALSE;
@@ -1025 +1134 @@
     return kTRUE;
 }
+

trunk/MagicSoft/Mars/manalysis/MCT1FindSupercuts.h

@@ -89 +89 @@
 
   Bool_t DefineTrainMatrix(const TString &name, const Int_t howmany,
-                           MH3 &href);
+                           MH3 &href, const TString &filetrain); 
 
   Bool_t DefineTestMatrix(const TString &name, const Int_t howmany,
-                          MH3 &href);
+                          MH3 &href,  const TString &filetest);
 
   Bool_t DefineTrainTestMatrix(const TString &name, 
-                               const Int_t howmanytrain, MH3 &hreftrain,
-                               const Int_t howmanytest,  MH3 &hreftest);
+                         const Int_t howmanytrain, MH3 &hreftrain,
+                         const Int_t howmanytest,  MH3 &hreftest,
+                         const TString &filetrain, const TString &filetest);
+
+  MHMatrix *GetMatrixTrain() { return fMatrixTrain; }
+  MHMatrix *GetMatrixTest()  { return fMatrixTest;  }
 
   Bool_t ReadMatrix( const TString &filetrain, const TString &filetest);
-  Bool_t WriteMatrix(const TString &filetrain, const TString &filetest);
 
-  Bool_t FindParams();
+  Bool_t FindParams(TString parSCinit, TArrayD &params, TArrayD &steps);
   Bool_t TestParams();
 

trunk/MagicSoft/Mars/manalysis/MCT1Supercuts.cc

@@ -50 +50 @@
 //
 MCT1Supercuts::MCT1Supercuts(const char *name, const char *title)
-    : fParameters(72),
-    fLengthUp(fParameters.GetArray()), fLengthLo(fParameters.GetArray()+8),
+  : fParameters(104), fStepsizes(104),
+    fLengthUp(fParameters.GetArray()),   fLengthLo(fParameters.GetArray()+8),
     fWidthUp(fParameters.GetArray()+16), fWidthLo(fParameters.GetArray()+24),
-    fDistUp(fParameters.GetArray()+32), fDistLo(fParameters.GetArray()+40),
-    fAsymUp(fParameters.GetArray()+48), fAsymLo(fParameters.GetArray()+56),
-    fAlphaUp(fParameters.GetArray()+64)
+    fDistUp(fParameters.GetArray()+32),  fDistLo(fParameters.GetArray()+40),
+    fAsymUp(fParameters.GetArray()+48),  fAsymLo(fParameters.GetArray()+56),
+    fConcUp(fParameters.GetArray()+64),  fConcLo(fParameters.GetArray()+72),
+    fLeakage1Up(fParameters.GetArray()+80), fLeakage1Lo(fParameters.GetArray()+88),
+    fAlphaUp(fParameters.GetArray()+96)
 {
     fName  = name  ? name  : "MCT1Supercuts";
@@ -71 +73 @@
 void MCT1Supercuts::InitParameters()
 {
+    //---------------------------------------------------
+    //  these are Daniel's original values for Mkn 421
+
     fLengthUp[0] =  0.315585;
     fLengthUp[1] =  0.001455;
@@ -80 +85 @@
     fLengthUp[7] = -0.013463;
 
+    fLengthLo[0] =  0.151530;
+    fLengthLo[1] =  0.028323;
+    fLengthLo[2] =  0.510707;
+    fLengthLo[3] =  0.053089;
+    fLengthLo[4] =  0.013708;
+    fLengthLo[5] =  2.357993;
+    fLengthLo[6] =  0.000080;
+    fLengthLo[7] = -0.007157;
+
     fWidthUp[0] =  0.145412;
     fWidthUp[1] = -0.001771;
@@ -89 +103 @@
     fWidthUp[7] = -0.016703;
 
+    fWidthLo[0] =  0.089187;
+    fWidthLo[1] = -0.006430;
+    fWidthLo[2] =  0.074442;
+    fWidthLo[3] =  0.003738;
+    fWidthLo[4] = -0.004256;
+    fWidthLo[5] = -0.014101;
+    fWidthLo[6] =  0.006126;
+    fWidthLo[7] = -0.002849;
+
     fDistUp[0] =  1.787943;
     fDistUp[1] =  0;
@@ -98 +121 @@
     fDistUp[7] =  0;
 
-    fLengthLo[0] =  0.151530;
-    fLengthLo[1] =  0.028323;
-    fLengthLo[2] =  0.510707;
-    fLengthLo[3] =  0.053089;
-    fLengthLo[4] =  0.013708;
-    fLengthLo[5] =  2.357993;
-    fLengthLo[6] =  0.000080;
-    fLengthLo[7] = -0.007157;
-
-    fWidthLo[0] =  0.089187;
-    fWidthLo[1] = -0.006430;
-    fWidthLo[2] =  0.074442;
-    fWidthLo[3] =  0.003738;
-    fWidthLo[4] = -0.004256;
-    fWidthLo[5] = -0.014101;
-    fWidthLo[6] =  0.006126;
-    fWidthLo[7] = -0.002849;
-
     fDistLo[0] =  0.589406;
     fDistLo[1] =  0;
@@ -124 +129 @@
     fDistLo[6] = -0.001750;
     fDistLo[7] =  0;
-
-    fAsymUp[0] =  0.061267;
-    fAsymUp[1] =  0.014462;
-    fAsymUp[2] =  0.014327;
-    fAsymUp[3] =  0.014540;
-    fAsymUp[4] =  0.013391;
-    fAsymUp[5] =  0.012319;
-    fAsymUp[6] =  0.010444;
-    fAsymUp[7] =  0.008328;
-
-    fAsymLo[0] = -0.012055;
-    fAsymLo[1] =  0.009157;
-    fAsymLo[2] =  0.005441;
-    fAsymLo[3] =  0.000399;
-    fAsymLo[4] =  0.001433;
-    fAsymLo[5] = -0.002050;
-    fAsymLo[6] = -0.000104;
-    fAsymLo[7] = -0.001188;
+    
+
+    // dummy values
+
+    fAsymUp[0] =  1.e10;
+    fAsymUp[1] =  0.0;
+    fAsymUp[2] =  0.0;
+    fAsymUp[3] =  0.0;
+    fAsymUp[4] =  0.0;
+    fAsymUp[5] =  0.0;
+    fAsymUp[6] =  0.0;
+    fAsymUp[7] =  0.0;
+
+    fAsymLo[0] = -1.e10;
+    fAsymLo[1] =  0.0;
+    fAsymLo[2] =  0.0;
+    fAsymLo[3] =  0.0;
+    fAsymLo[4] =  0.0;
+    fAsymLo[5] =  0.0;
+    fAsymLo[6] =  0.0;
+    fAsymLo[7] =  0.0;
+
+    fConcUp[0] =  1.e10;
+    fConcUp[1] =  0.0;
+    fConcUp[2] =  0.0;
+    fConcUp[3] =  0.0;
+    fConcUp[4] =  0.0;
+    fConcUp[5] =  0.0;
+    fConcUp[6] =  0.0;
+    fConcUp[7] =  0.0;
+
+    fConcLo[0] = -1.e10;
+    fConcLo[1] =  0.0;
+    fConcLo[2] =  0.0;
+    fConcLo[3] =  0.0;
+    fConcLo[4] =  0.0;
+    fConcLo[5] =  0.0;
+    fConcLo[6] =  0.0;
+    fConcLo[7] =  0.0;
+
+    fLeakage1Up[0] =  1.e10;
+    fLeakage1Up[1] =  0.0;
+    fLeakage1Up[2] =  0.0;
+    fLeakage1Up[3] =  0.0;
+    fLeakage1Up[4] =  0.0;
+    fLeakage1Up[5] =  0.0;
+    fLeakage1Up[6] =  0.0;
+    fLeakage1Up[7] =  0.0;
+
+    fLeakage1Lo[0] = -1.e10;
+    fLeakage1Lo[1] =  0.0;
+    fLeakage1Lo[2] =  0.0;
+    fLeakage1Lo[3] =  0.0;
+    fLeakage1Lo[4] =  0.0;
+    fLeakage1Lo[5] =  0.0;
+    fLeakage1Lo[6] =  0.0;
+    fLeakage1Lo[7] =  0.0;
 
     fAlphaUp[0] = 13.123440; 
@@ -151 +195 @@
     fAlphaUp[6] = 0;
     fAlphaUp[7] = 0;
+
+    //---------------------------------------------------
+    // fStepsizes 
+    // if == 0.0    the parameter will be fixed in the minimization
+    //    != 0.0    initial step sizes for the parameters
+
+    // LengthUp
+    fStepsizes[0] = 0.03;
+    fStepsizes[1] = 0.0002;
+    fStepsizes[2] = 0.02;
+    fStepsizes[3] = 0.0006;
+    fStepsizes[4] = 0.0002;
+    fStepsizes[5] = 0.002;
+    fStepsizes[6] = 0.0008;
+    fStepsizes[7] = 0.002;
+
+    // LengthLo
+    fStepsizes[8]  = 0.02;
+    fStepsizes[9]  = 0.003;
+    fStepsizes[10] = 0.05;
+    fStepsizes[11] = 0.006;
+    fStepsizes[12] = 0.002;
+    fStepsizes[13] = 0.3;
+    fStepsizes[14] = 0.0001;
+    fStepsizes[15] = 0.0008;
+
+    // WidthUp
+    fStepsizes[16] = 0.02;
+    fStepsizes[17] = 0.0002;
+    fStepsizes[18] = 0.006;
+    fStepsizes[19] = 0.003;
+    fStepsizes[20] = 0.002;
+    fStepsizes[21] = 0.006;
+    fStepsizes[22] = 0.002;
+    fStepsizes[23] = 0.002;
+
+    // WidthLo
+    fStepsizes[24] = 0.009;
+    fStepsizes[25] = 0.0007;
+    fStepsizes[26] = 0.008;
+    fStepsizes[27] = 0.0004;
+    fStepsizes[28] = 0.0005;
+    fStepsizes[29] = 0.002;
+    fStepsizes[30] = 0.0007;
+    fStepsizes[31] = 0.003;
+
+    // DistUp
+    fStepsizes[32] = 0.2;
+    fStepsizes[33] = 0.0;
+    fStepsizes[34] = 0.3;
+    fStepsizes[35] = 0.02;
+    fStepsizes[36] = 0.0;
+    fStepsizes[37] = 0.03;
+    fStepsizes[38] = 0.02;
+    fStepsizes[39] = 0.0;
+
+    // DistLo
+    fStepsizes[40] = 0.06;
+    fStepsizes[41] = 0.0;
+    fStepsizes[42] = 0.009;
+    fStepsizes[43] = 0.0008;
+    fStepsizes[44] = 0.0;
+    fStepsizes[45] = 0.005;
+    fStepsizes[46] = 0.0002;
+    fStepsizes[47] = 0.0;
+
+    // AsymUp
+    fStepsizes[48] = 0.0;
+    fStepsizes[49] = 0.0;
+    fStepsizes[50] = 0.0;
+    fStepsizes[51] = 0.0;
+    fStepsizes[52] = 0.0;
+    fStepsizes[53] = 0.0;
+    fStepsizes[54] = 0.0;
+    fStepsizes[55] = 0.0;
+
+    // AsymLo
+    fStepsizes[56] = 0.0;
+    fStepsizes[57] = 0.0;
+    fStepsizes[58] = 0.0;
+    fStepsizes[59] = 0.0;
+    fStepsizes[60] = 0.0;
+    fStepsizes[61] = 0.0;
+    fStepsizes[62] = 0.0;
+    fStepsizes[63] = 0.0;
+
+    // ConcUp
+    fStepsizes[64] = 0.0;
+    fStepsizes[65] = 0.0;
+    fStepsizes[66] = 0.0;
+    fStepsizes[67] = 0.0;
+    fStepsizes[68] = 0.0;
+    fStepsizes[69] = 0.0;
+    fStepsizes[70] = 0.0;
+    fStepsizes[71] = 0.0;
+
+    // ConcLo
+    fStepsizes[72] = 0.0;
+    fStepsizes[73] = 0.0;
+    fStepsizes[74] = 0.0;
+    fStepsizes[75] = 0.0;
+    fStepsizes[76] = 0.0;
+    fStepsizes[77] = 0.0;
+    fStepsizes[78] = 0.0;
+    fStepsizes[79] = 0.0;
+
+    // Leakage1Up
+    fStepsizes[80] = 0.0;
+    fStepsizes[81] = 0.0;
+    fStepsizes[82] = 0.0;
+    fStepsizes[83] = 0.0;
+    fStepsizes[84] = 0.0;
+    fStepsizes[85] = 0.0;
+    fStepsizes[86] = 0.0;
+    fStepsizes[87] = 0.0;
+
+    // Leakage1Lo
+    fStepsizes[88] = 0.0;
+    fStepsizes[89] = 0.0;
+    fStepsizes[90] = 0.0;
+    fStepsizes[91] = 0.0;
+    fStepsizes[92] = 0.0;
+    fStepsizes[93] = 0.0;
+    fStepsizes[94] = 0.0;
+    fStepsizes[95] = 0.0;
+
+    // AlphaUp
+    fStepsizes[96]  = 0.0;
+    fStepsizes[97]  = 0.0;
+    fStepsizes[98]  = 0.0;
+    fStepsizes[99]  = 0.0;
+    fStepsizes[100] = 0.0;
+    fStepsizes[101] = 0.0;
+    fStepsizes[102] = 0.0;
+    fStepsizes[103] = 0.0;
 }
 
@@ -173 +352 @@
 }
 
-
-
-
-
-
-
-
-
-
-
-
+// --------------------------------------------------------------------------
+//
+// Set the step sizes from the array 'd'
+//
+//
+Bool_t MCT1Supercuts::SetStepsizes(const TArrayD &d)
+{
+    if (d.GetSize() != fStepsizes.GetSize())
+    {
+        *fLog << err << "Sizes of d and of fStepsizes are different : "
+              << d.GetSize() << ",  " << fStepsizes.GetSize() << endl;
+        return kFALSE;
+    }
+
+    fStepsizes = d;
+
+    return kTRUE;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

trunk/MagicSoft/Mars/manalysis/MCT1Supercuts.h

@@ -14 +14 @@
 private:
     TArrayD fParameters; // supercut parameters
+    TArrayD fStepsizes;  // step sizes of supercut parameters
 
     Double_t *fLengthUp; //!
@@ -23 +24 @@
     Double_t *fAsymUp;   //!
     Double_t *fAsymLo;   //!
+
+    Double_t *fConcUp;   //!
+    Double_t *fConcLo;   //!
+    Double_t *fLeakage1Up;   //!
+    Double_t *fLeakage1Lo;   //!
+
     Double_t *fAlphaUp;  //!
+
 
 public:
@@ -31 +39 @@
 
     Bool_t SetParameters(const TArrayD &d);
+    Bool_t SetStepsizes(const TArrayD &d);
+
     const TArrayD &GetParameters() const { return fParameters; }
+    const TArrayD &GetStepsizes()  const { return fStepsizes;  }
 
     const Double_t *GetLengthUp() const { return fLengthUp; }
@@ -41 +52 @@
     const Double_t *GetAsymUp() const   { return fAsymUp; }
     const Double_t *GetAsymLo() const   { return fAsymLo; }
+
+    const Double_t *GetConcUp() const   { return fConcUp; }
+    const Double_t *GetConcLo() const   { return fConcLo; }
+
+    const Double_t *GetLeakage1Up() const   { return fLeakage1Up; }
+    const Double_t *GetLeakage1Lo() const   { return fLeakage1Lo; }
+
     const Double_t *GetAlphaUp() const  { return fAlphaUp; }
 
@@ -47 +65 @@
 
 #endif
+
+
+

trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.cc

@@ -44 +44 @@
 #include "MHillasExt.h"
 #include "MHillasSrc.h"
+#include "MNewImagePar.h"
 #include "MMcEvt.hxx"
 #include "MCerPhotEvt.h"
@@ -65 +66 @@
 
 MCT1SupercutsCalc::MCT1SupercutsCalc(const char *hilname, 
-                                     const char *hilsrcname, const char *name, const char *title)
+                                     const char *hilsrcname, 
+                                     const char *name, const char *title)
   : fHadronnessName("MHadronness"), fHilName(hilname), fHilSrcName(hilsrcname),
-    fSuperName("MCT1Supercuts")
+    fHilExtName("MHillasExt"), fNewParName("MNewImagePar"), 
+    fSuperName("MCT1Supercuts") 
 {
     fName  = name  ? name  : "MCT1SupercutsCalc";
@@ -101 +104 @@
         return kFALSE;
     }
-
 
     if (fMatrix)
@@ -114 +116 @@
     }
 
+    fHilExt = (MHillasExt*)pList->FindObject(fHilExtName, "MHillasExt");
+    if (!fHilExt)
+    {
+        *fLog << err << fHilExtName << " [MHillasExt] not found... aborting." << endl;
+        return kFALSE;
+    }
+
     fHilSrc = (MHillasSrc*)pList->FindObject(fHilSrcName, "MHillasSrc");
     if (!fHilSrc)
     {
         *fLog << err << fHilSrcName << " [MHillasSrc] not found... aborting." << endl;
+        return kFALSE;
+    }
+
+    fNewPar = (MNewImagePar*)pList->FindObject(fNewParName, "MNewImagePar");
+    if (!fNewPar)
+    {
+        *fLog << err << fNewParName << " [MNewImagePar] not found... aborting." << endl;
         return kFALSE;
     }
@@ -172 +188 @@
 Double_t MCT1SupercutsCalc::GetVal(Int_t i) const
 {
-    return (*fMatrix)[fMap[i]];
+
+    Double_t val = (*fMatrix)[fMap[i]];
+
+
+    //*fLog << "MCT1SupercutsCalc::GetVal; i, fMatrix, fMap, val = "
+    //    << i << ",  " << fMatrix << ",  " << fMap[i] << ",  " << val << endl;
+
+
+    return val;
 }
 
@@ -187 +211 @@
 {
     if (fMatrix)
-        return;
+      return;
 
     fMatrix = mat;
@@ -199 +223 @@
     fMap[6] = fMatrix->AddColumn("MHillasSrc.fDist");
     fMap[7] = fMatrix->AddColumn("fabs(MHillasSrc.fAlpha)");
+    fMap[8] = fMatrix->AddColumn("sgn(MHillasSrc.fCosDeltaAlpha)*(MHillasExt.fM3Long)");
+    fMap[9] = fMatrix->AddColumn("MNewImagePar.fConc");
+    fMap[10]= fMatrix->AddColumn("MNewImagePar.fLeakage1");
 }
 
 // ---------------------------------------------------------------------------
 //
-// Evaluate dynamical supercuts for CT1 Mkn421 2001 (Daniel Kranich)
-// optimized for mkn 421 2001 data
+// Evaluate dynamical supercuts 
 // 
 //          set hadronness to 0.25 if cuts are fullfilled
@@ -222 +248 @@
     const Double_t dist0   = fMatrix ? GetVal(6) : fHilSrc->GetDist();
 
+    Double_t help;
+    if (!fMatrix)
+      help = TMath::Sign(fHilExt->GetM3Long(), 
+                         fHilSrc->GetCosDeltaAlpha());
+    const Double_t asym0   = fMatrix ? GetVal(8) : help;
+    const Double_t conc    = fMatrix ? GetVal(9) : fNewPar->GetConc();
+    const Double_t leakage = fMatrix ? GetVal(10): fNewPar->GetLeakage1();
+
     const Double_t newdist = dist0 * fMm2Deg;
 
@@ -236 +270 @@
     const Double_t length  = length0 * fMm2Deg;
     const Double_t width   = width0  * fMm2Deg;
-
-    if (newdist < 1.05                                                     &&
+    const Double_t asym    = asym0   * fMm2Deg;
+
+    /*
+    *fLog << "newdist, length, width, asym, dist, conc, leakage = " 
+          << newdist << ",  " << length << ",  " << width << ",  "
+          << asym    << ",  " << dist   << ",  " << conc  << ",  " << leakage
+          << endl;
+  
+    *fLog << "upper cuts in newdist, length, width, asym, dist, conc, leakage = " 
+          << CtsMCut (fSuper->GetDistUp(),   dmls, dmcza, dmls2, dd2) << ",  "
+          << CtsMCut (fSuper->GetDistLo(),   dmls, dmcza, dmls2, dd2) << ",  "
+
+          << CtsMCut (fSuper->GetLengthUp(),   dmls, dmcza, dmls2, dd2) << ",  "
+          << CtsMCut (fSuper->GetLengthLo(),   dmls, dmcza, dmls2, dd2) << ",  "
+
+          << CtsMCut (fSuper->GetWidthUp(),   dmls, dmcza, dmls2, dd2) << ",  "
+          << CtsMCut (fSuper->GetWidthLo(),   dmls, dmcza, dmls2, dd2) << ",  "
+
+          << CtsMCut (fSuper->GetAsymUp(),   dmls, dmcza, dmls2, dd2) << ",  "
+          << CtsMCut (fSuper->GetAsymLo(),   dmls, dmcza, dmls2, dd2) << ",  "
+
+          << CtsMCut (fSuper->GetDistUp(),   dmls, dmcza, dmls2, dd2) << ",  "
+          << CtsMCut (fSuper->GetDistLo(),   dmls, dmcza, dmls2, dd2) << ",  "
+
+          << CtsMCut (fSuper->GetConcUp(),   dmls, dmcza, dmls2, dd2) << ",  "
+          << CtsMCut (fSuper->GetConcLo(),   dmls, dmcza, dmls2, dd2) << ",  "
+
+          << CtsMCut (fSuper->GetLeakage1Up(),   dmls, dmcza, dmls2, dd2) << ",  "
+          << CtsMCut (fSuper->GetLeakage1Lo(),   dmls, dmcza, dmls2, dd2) << ",  "
+          << endl;
+    */
+
+
+    if (
+        //dist    < 1.05                                                     &&
+        //newdist < 1.05                                                     &&
+
         newdist < CtsMCut (fSuper->GetDistUp(),   dmls, dmcza, dmls2, dd2) &&
         newdist > CtsMCut (fSuper->GetDistLo(),   dmls, dmcza, dmls2, dd2) &&
-        dist    < 1.05                                                     &&
+
         length  < CtsMCut (fSuper->GetLengthUp(), dmls, dmcza, dmls2, dd2) &&
         length  > CtsMCut (fSuper->GetLengthLo(), dmls, dmcza, dmls2, dd2) &&
+
         width   < CtsMCut (fSuper->GetWidthUp(),  dmls, dmcza, dmls2, dd2) &&
         width   > CtsMCut (fSuper->GetWidthLo(),  dmls, dmcza, dmls2, dd2) &&
-        //asym  < CtsMCut (fSuper->GetAsymUp(),   dmls, dmcza, dmls2, dd2) &&
-        //asym  > CtsMCut (fSuper->GetAsymLo(),   dmls, dmcza, dmls2, dd2) &&
+
+        asym    < CtsMCut (fSuper->GetAsymUp(),   dmls, dmcza, dmls2, dd2) &&
+        asym    > CtsMCut (fSuper->GetAsymLo(),   dmls, dmcza, dmls2, dd2) &&
+
         dist    < CtsMCut (fSuper->GetDistUp(),   dmls, dmcza, dmls2, dd2) &&
-        dist    > CtsMCut (fSuper->GetDistLo(),   dmls, dmcza, dmls2, dd2)  )
-        fHadronness->SetHadronness(0.25);
+        dist    > CtsMCut (fSuper->GetDistLo(),   dmls, dmcza, dmls2, dd2) &&
+
+        conc    < CtsMCut (fSuper->GetConcUp(),   dmls, dmcza, dmls2, dd2) &&
+        conc    > CtsMCut (fSuper->GetConcLo(),   dmls, dmcza, dmls2, dd2) &&
+
+        leakage < CtsMCut (fSuper->GetLeakage1Up(),dmls, dmcza, dmls2, dd2) &&
+        leakage > CtsMCut (fSuper->GetLeakage1Lo(),dmls, dmcza, dmls2, dd2)  ) 
+
+      fHadronness->SetHadronness(0.25);
     else
-        fHadronness->SetHadronness(0.75);
+      fHadronness->SetHadronness(0.75);
+
+    //*fLog << "SChadroness = " << fHadronness->GetHadronness() << endl;
 
     fHadronness->SetReadyToSave();
@@ -257 +338 @@
     return kTRUE;
 }
+
+
+
+
+
+
+
+
+

trunk/MagicSoft/Mars/manalysis/MCT1SupercutsCalc.h

@@ -13 +13 @@
 class MHillas;
 class MHillasSrc;
+class MHillasExt;
+class MNewImagePar;
 class MMcEvt;
 class MCerPhotEvt;
@@ -25 +27 @@
     MHillas       *fHil;
     MHillasSrc    *fHilSrc;
+    MHillasExt    *fHilExt;
+    MNewImagePar  *fNewPar;
     MMcEvt        *fMcEvt;
     MHadronness   *fHadronness; //! output container for hadronness
@@ -32 +36 @@
     TString  fHilName;
     TString  fHilSrcName;
+    TString  fHilExtName;
+    TString  fNewParName;
     TString  fSuperName;        // name of container for supercut parameters
 
     Double_t fMm2Deg;           //!
 
-    Int_t     fMap[8];          //!
+    Int_t     fMap[11];         //!
     MHMatrix *fMatrix;          //!
 
@@ -63 +69 @@
 
 #endif
+
+
+
+
+
+
+
+
+
+
+
+

trunk/MagicSoft/Mars/manalysis/MMinuitInterface.cc

@@ -103 +103 @@
     //  3   maximum output, showing progress of minimizations
     //
-    minuit.SetPrintLevel(1);
+    minuit.SetPrintLevel(-1);
 
     //..............................................
@@ -175 +175 @@
 
     //..............................................
+    // This doesn't seem to have any effect
     // Set maximum number of iterations (default = 500)
     //Int_t maxiter = 100000;
-    minuit.SetMaxIterations(10);
+    //minuit.SetMaxIterations(maxiter);
 
     //..............................................
@@ -211 +212 @@
         if (nulloutput)
             fLog->SetNullOutput(kTRUE);
-        Double_t tmp = 0;
-        minuit.mnexcm("SIMPLEX", &tmp, 0, fErrMinimize);
+        Int_t    maxcalls  = 3000;
+        Double_t tolerance = 0.1;
+        Double_t tmp[2];
+        tmp[0] = maxcalls;
+        tmp[1] = tolerance;
+        minuit.mnexcm("SIMPLEX", &tmp[0], 2, fErrMinimize);
         if (nulloutput)
             fLog->SetNullOutput(kFALSE);
-        //*fLog << "return from SIMPLEX" << endl;
+        *fLog << "return from SIMPLEX" << endl;
     }
 
@@ -269 +274 @@
     //            2    values, errors, step sizes, internal values
     //            3    values, errors, step sizes, 1st derivatives
-    //            4    values, paraboloc errors, MINOS errors
+    //            4    values, parabolic errors, MINOS errors
 
     //Int_t nkode = 4;
@@ -292 +297 @@
     return kTRUE;
 }
+
+
+
+
+
+
+
+
+