Changeset 6979 for trunk/MagicSoft

Timestamp:

04/27/05 16:46:24 (20 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• Makefile (modified) (1 diff)
• NEWS (modified) (1 diff)
• callisto.cc (modified) (8 diffs)
• ganymed.cc (modified) (1 diff)
• mastro/MAstroCamera.cc (modified) (4 diffs)
• mastro/MAstroCamera.h (modified) (1 diff)
• mbase/MStatusArray.h (modified) (1 diff)
• mcalib/MCalibColorSet.cc (modified) (3 diffs)
• mjobs/MJCalibTest.cc (modified) (6 diffs)
• mjobs/MJCalibTest.h (modified) (1 diff)
• mjobs/MJCalibrateSignal.cc (modified) (6 diffs)
• mjobs/MJCalibrateSignal.h (modified) (2 diffs)
• mjobs/MJCalibration.cc (modified) (11 diffs)
• mjobs/MJCalibration.h (modified) (2 diffs)
• mjobs/MJCut.cc (modified) (1 diff)
• mjobs/MJCut.h (modified) (1 diff)
• mjobs/MJPedestal.cc (modified) (7 diffs)
• mjobs/MJPedestal.h (modified) (1 diff)
• mmuon/MHMuonPar.cc (modified) (15 diffs)
• mmuon/MHMuonPar.h (modified) (3 diffs)
• mmuon/MHSingleMuon.cc (added)
• mmuon/MHSingleMuon.h (added)
• mmuon/MMuonCalibPar.cc (modified) (5 diffs)
• mmuon/MMuonCalibPar.h (modified) (1 diff)
• mmuon/MMuonCalibParCalc.cc (modified) (6 diffs)
• mmuon/MMuonCalibParCalc.h (modified) (1 diff)
• mmuon/MMuonSearchPar.cc (modified) (9 diffs)
• mmuon/MMuonSearchPar.h (modified) (3 diffs)
• mmuon/MMuonSearchParCalc.cc (modified) (6 diffs)
• mmuon/MMuonSearchParCalc.h (modified) (2 diffs)
• mmuon/MMuonSetup.cc (modified) (3 diffs)
• mmuon/MMuonSetup.h (modified) (1 diff)
• mmuon/Makefile (modified) (1 diff)
• mmuon/MuonLinkDef.h (modified) (1 diff)
• star.rc (modified) (1 diff)

trunk/MagicSoft/Mars/Changelog

@@ -21 +21 @@
 
                                                  -*-*- END OF LINE -*-*-
+
+ 2005/04/27 Thomas Bretz
+
+   * Makefile:
+     - added mmuon
+     - remobed mstarcam
+
+   * callisto.cc, ganymed.cc, star.cc:
+     - renamed ProcessFile to Process
+
+   * star.rc:
+     - added some muon parameters
+
+   * mastro/MAstroCamera.[h,cc]:
+     - temporarily removed interface to MStarPos
+
+   * mbase/MStatusArray.h:
+     - added default constructor
+
+   * mcalib/MCalibColorSet.cc:
+     - added runs 39942, 39944, 44834, 39941, 39943 and 44833
+       (undocumented change from the BCN cvs)
+
+   * mjobs/MJCalib.[h,cc], mjobs/MJCalibTest.[h,cc],
+     mjobs/MJCalibrateSignal.[h,cc], mjobs/MJCalibration.[h,cc],
+     mjobs/MJCut.[h,cc], mjobs/MJPedestal.[h,cc]:
+     - removed support for MRunIter (use the setter of MSequence
+       instead) -- this makes the code a lot easier to maintain
+     - removed support for autodetection if the output already exists --
+       this makes the code a lot easier to maintain
+     - renamed ProcessFile to Process - which was missleading
+
+   * mmuon/MHMuonPar.[h,cc]:
+     - changes to axis labels etc.
+
+   * mmuon/MMuonCalibPar.[h,cc]:
+     - removed the histograms and all obsolete variables
+     - removed obsolete SetUseUnmap (this cannot happen 
+       by definition of Unmap)
+
+   * mmuon/MMuonCalibParCalc.[h,cc]:
+     - moved the code for calculation the parameters to new class
+       MHSingleMuon
+
+   * mmuon/MMuonSearchPar.[h,cc]:
+     - replaced arbitrary fir by minuit (faster and more accurate)
+     - removed precalculation of muon center - makes fit worse
+
+   * mmuon/MMuonSearchParCalc.[h,cc]:
+     - fixes to comments
+     - fixes to includes
+
+   * mmuon/MMuonSetup.[h,cc]:
+     - binnings removed (replaces by MBinning)
+
+   * mmuon/Makefile, mmuon/MuonLinkDef.h:
+     - added MHSingleMuon
+
+
 
  2005/04/25 Thomas Bretz

trunk/MagicSoft/Mars/Makefile

@@ -68 +68 @@
           mjobs \
           mtools \
-          mstarcam
+          mmuon
 
 #LIBRARIES = $(SUBDIRS:%=lib/lib%.a)

trunk/MagicSoft/Mars/NEWS

@@ -25 +25 @@
      excess events versus size. The energy estimation is done in
      MJSpectrum (sponde)
+
+   - added support for the runs 39942, 39944, 44834, 39941, 39943, 44833
+     in the calibration (MCalibColorSet)
+
+   - MJCalibration.MHCalibrationChargeCam.ProbLimit has been set to
+     1e-18 in callisto_Dec04_Jan05.rc
+
+   - ProcessFile has been renamed to Process in all job classes, because
+     ProcessFile is missleading
+
+   - support for MRunIter has been removed from the job classes (use
+     the setter functions of MSeqeunce instead)
+
+   - added muon support to star
 
 

trunk/MagicSoft/Mars/callisto.cc

@@ -360 +360 @@
         // job1.SetPathIn(kInpathC);   // not yet needed
 
-        if (!job1.ProcessFile())
+        if (!job1.Process())
         {
             gLog << err << "Calculation of pedestal failed." << endl << endl;
@@ -396 +396 @@
         // Please check the Changelog of 2005/04/20 about further deatils of the next comment
         //if (job1.GetExtractor().InheritsFrom("MExtractTimeAndCharge"))
-        if (!job2.ProcessFile())
+        if (!job2.Process())
         {
             gLog << err << "Calculation of pedestal resolution failed." << endl << endl;
@@ -426 +426 @@
         job3.SetExtractorCam(job2.GetPedestalCam());
 
-        if (!job3.ProcessFile(job1.GetPedestalCam()))
+        if (!job3.Process(job1.GetPedestalCam()))
         {
             gLog << err << "Calculation of calibration failed." << endl << endl;
@@ -454 +454 @@
             job4.SetDataType(kDataType);
 
-            if (!job4.ProcessFile(job1.GetPedestalCam()))
+            if (!job4.Process(job1.GetPedestalCam()))
             {
                 gLog << err << "Calibration of calibration failed." << endl << endl;
@@ -494 +494 @@
         job1.SetUseHists(kMoon);
 
-        if (!job1.ProcessFile())
+        if (!job1.Process())
         {
             gLog << err << "Calculation of fundamental pedestal failed." << endl << endl;
@@ -532 +532 @@
         // Please check the Changelog of 2005/04/20 about further deatils of the next comment
         //if (job1.GetExtractor().InheritsFrom("MExtractTimeAndCharge"))
-        if (!job2.ProcessFile())
+        if (!job2.Process())
         {
             gLog << err << "Calculation of pedestal from extrtactor (random) failed." << endl << endl;
@@ -568 +568 @@
         job3.SetUseHists(kMoon);
 
-        if (!job3.ProcessFile())
+        if (!job3.Process())
         {
             gLog << err << "Calculation of pedestal from extractor failed." << endl << endl;
@@ -597 +597 @@
 
         // Where to search for calibration files
-        if (!job4.ProcessFile(job1.GetPedestalCam(), job2.GetPedestalCam(), job3.GetPedestalCam()))
+        if (!job4.Process(job1.GetPedestalCam(), job2.GetPedestalCam(), job3.GetPedestalCam()))
             return 2;
 

trunk/MagicSoft/Mars/ganymed.cc

@@ -250 +250 @@
     //    job.EnableStorageOfResult();
 
-    if (!job.ProcessFile(seq))
+    if (!job.Process(seq))
     {
         gLog << err << "Calculation of cuts failed." << endl << endl;

trunk/MagicSoft/Mars/mastro/MAstroCamera.cc

@@ -93 +93 @@
                       // HOW TO GET RID OF IT? Move MHCamera to mgeom?
 
-#include "MStarPos.h"
+//#include "MStarPos.h"
 
 ClassImp(MAstroCamera);
@@ -429 +429 @@
 // Currently, the mean spot when averaging over all mirrors is used.
 //
+/*
 void MAstroCamera::FillStarList(TList* list)
 {
@@ -464 +465 @@
       }
       mean *= 1./num;
+
       MStarPos *starpos = new MStarPos;
       starpos->SetExpValues(mag,mean(0),mean(1));
@@ -477 +479 @@
     }
 }
+*/
 
 // ------------------------------------------------------------------------

trunk/MagicSoft/Mars/mastro/MAstroCamera.h

@@ -35 +35 @@
     void SetGeom(const MGeomCam &cam);
 
-    void FillStarList(TList *list);
+    //void FillStarList(TList *list);
 
     void GetDiffZdAz(Double_t x, Double_t y, Double_t &dzd, Double_t &daz);

trunk/MagicSoft/Mars/mbase/MStatusArray.h

@@ -20 +20 @@
 
 public:
+    MStatusArray() : TObjArray() { }
     TObject *DisplayIn(Option_t *o=0) const;         // *MENU*
     void     DisplayIn(MStatusDisplay &d, const char *tab=0) const;

trunk/MagicSoft/Mars/mcalib/MCalibColorSet.cc

@@ -203 +203 @@
     case 34814:
     case 35415:
+    case 39942:
+    case 39944:
     case 44768:
     case 44976:
@@ -261 +263 @@
     case 26568:
     case 26924:
+    case 44834:
     case 45051:
     case 45084:
@@ -277 +280 @@
       break;
       
+    case 39941:
+    case 39943:
+    case 44833:
     case 45086:
     case 45088:

trunk/MagicSoft/Mars/mjobs/MJCalibTest.cc

@@ -95 +95 @@
 // Default constructor. 
 //
-// Sets fUseCosmicsFilter to kTRUE, fRuns to 0, fExtractor to NULL, fTimeExtractor to NULL
+// Sets fUseCosmicsFilter to kTRUE, fExtractor to NULL, fTimeExtractor to NULL
 // fDisplay to kNormalDisplay
 //
@@ -210 +210 @@
 const char* MJCalibTest::GetOutputFile() const
 {
-    const TString name(GetOutputFileName());
-    if (name.IsNull())
-        return "";
-
-    return Form("%s/%s", fPathOut.Data(), name.Data());
-}
-
-
-const char* MJCalibTest::GetOutputFileName() const
-{
-
-  if (fSequence.IsValid())
-    return Form("calib%08d.root", fSequence.GetSequence());
-  
-  if (!fRuns)
-    return "";
-  
-  return Form("%s-F1.root", (const char*)fRuns->GetRunsAsFileName());
-  
+    return Form("%s/calib%08d.root", fPathOut.Data(), fSequence.GetSequence());
 }
 
@@ -279 +261 @@
 }
 
-Bool_t MJCalibTest::ProcessFile(MPedestalCam &pedcam)
+Bool_t MJCalibTest::Process(MPedestalCam &pedcam)
 {
-  
-  if (!fSequence.IsValid())
-    {
-      if (!fRuns)
-        {
-            *fLog << err << "ERROR - Sequence invalid and no runs chosen!" << endl;
-            return kFALSE;
-        }
-      
-      if (fRuns->GetNumRuns() != fRuns->GetNumEntries())
-        {
-          *fLog << err << "Number of files found doesn't match number of runs... abort."
-                << fRuns->GetNumRuns() << " vs. " << fRuns->GetNumEntries() << endl;
-          return kFALSE;
-        }
-      *fLog << "Calibrate data from ";
-      *fLog << "Runs " << fRuns->GetRunsAsString() << endl;
-      *fLog << endl;
-    }
-  
+    if (!fSequence.IsValid())
+    {
+          *fLog << err << "ERROR - Sequence invalid..." << endl;
+          return kFALSE;
+    }
+
+    *fLog << inf;
+    fLog->Separator(GetDescriptor());
+    *fLog << "Calculate calibration test from Sequence #";
+    *fLog << fSequence.GetSequence() << endl << endl;
+
   CheckEnv();
   
@@ -307 +280 @@
   
   *fLog << "Calibrate Calibration data from ";
-  if (fSequence.IsValid())
-    *fLog << "Sequence #" << fSequence.GetSequence() << endl;
-  else
-    *fLog << "Runs " << fRuns->GetRunsAsString() << endl;
+  *fLog << "Sequence #" << fSequence.GetSequence() << endl;
   *fLog << endl;
   
@@ -415 +385 @@
 
   if (IsUseRawData())
-    rawread.AddFiles(fSequence.IsValid() ? iter : *fRuns);
+    rawread.AddFiles(iter);
   else
-    static_cast<MRead&>(read).AddFiles(fSequence.IsValid() ? iter : *fRuns);
+    static_cast<MRead&>(read).AddFiles(iter);
 
   // Check for interleaved events
@@ -475 +445 @@
   MCalibrationTestCalc  testcalc;
 
-  if (!fSequence.IsValid())
-    {
-      testcalc.SetOutputPath(fPathOut);
-      testcalc.SetOutputFile(Form("%s-TestCalibStat.txt",(const char*)fRuns->GetRunsAsFileName()));
-    }
-  
   MHCamEvent evt0(0,"Signal", "Un-Calibrated Signal;;S [FADC cnts]" );
   MHCamEvent evt1(0,"CalSig", "Cal. and Interp. Sig. by Pixel Size Ratio;;S [phe]");

trunk/MagicSoft/Mars/mjobs/MJCalibTest.h

@@ -74 +74 @@
   void SetNormalDisplay()    { fDisplayType = kNormalDisplay;    }
   
-  Bool_t ProcessFile(MPedestalCam &pedcam);
+  Bool_t Process(MPedestalCam &pedcam);
   
   ClassDef(MJCalibTest, 0) // Tool to calibrate and test the calibration run itself

trunk/MagicSoft/Mars/mjobs/MJCalibrateSignal.cc

@@ -149 +149 @@
 }
 
-const char* MJCalibrateSignal::GetInputFile() const
-{
-    const TString name(GetInputFileName());
-    if (name.IsNull())
-        return "";
-
-    return Form("%s/%s", fPathIn.Data(), name.Data());
-}
-
-const char* MJCalibrateSignal::GetInputFileName() const
-{
-
-  if (fSequence.IsValid())
-    return Form("calib%08d.root", fSequence.GetSequence());
-  
-  //  if (!fCruns)
-    return "";
-  
-  //  return Form("%s-F1.root", (const char*)fCruns->GetRunsAsFileName());
-}
-
 Bool_t MJCalibrateSignal::WriteResult(TObjArray &cont)
 {
@@ -181 +160 @@
 Bool_t MJCalibrateSignal::ReadCalibration(TObjArray &l, MBadPixelsCam &cam, MExtractor* &ext2, MExtractor* &ext3, TString &geom) const
 {
-
-    const TString fname = GetInputFile();
+    TString fname = Form("%s/calib%08d.root", fPathIn.Data(), fSequence.GetSequence());
 
     *fLog << inf << "Reading from file: " << fname << endl;
@@ -274 +252 @@
 }
 
-Bool_t MJCalibrateSignal::ProcessFile(MPedestalCam &pedcamab, MPedestalCam &pedcambias, 
-                                      MPedestalCam &pedcamextr)
+Bool_t MJCalibrateSignal::Process(MPedestalCam &pedcamab, MPedestalCam &pedcambias,
+                                  MPedestalCam &pedcamextr)
 {
     if (!fSequence.IsValid())
     {
-        if (!fRuns)
-        {
-          *fLog << err << "ERROR - Sequence invalid and no runs chosen!" << endl;
+          *fLog << err << "ERROR - Sequence invalid..." << endl;
           return kFALSE;
-        }
-
-        if (fRuns->GetNumRuns() != fRuns->GetNumEntries())
-        {
-            *fLog << err << "Number of files found doesn't match number of runs... abort."
-                << fRuns->GetNumRuns() << " vs. " << fRuns->GetNumEntries() << endl;
-            return kFALSE;
-        }
-        *fLog << "Calibrate data from ";
-        *fLog << "Runs " << fRuns->GetRunsAsString() << endl;
-        *fLog << endl;
-    }
+    }
+
+    *fLog << inf;
+    fLog->Separator(GetDescriptor());
+    *fLog << "Calculate calibrated data from Sequence #";
+    *fLog << fSequence.GetSequence() << endl << endl;
+
 
     //if (!CheckEnv())
@@ -303 +274 @@
     // --------------------------------------------------------------------------------
 
-    *fLog << inf;
-    fLog->Separator(GetDescriptor());
-    *fLog << "Calculate calibrated data from runs ";
-    *fLog << fSequence.GetName() << endl;
-    *fLog << endl;
-
-    // --------------------------------------------------------------------------------
-
     MDirIter iter;
     if (fSequence.IsValid())
-      {
+    {
         const Int_t n0 = fSequence.SetupDatRuns(iter, fPathData, IsUseRawData());
         const Int_t n1 = fSequence.GetNumDatRuns();
         if (n0==0)
-          {
+        {
             *fLog << err << "ERROR - No input files of sequence found!" << endl;
             return kFALSE;
-          }
+        }
         if (n0!=n1)
-          {
+        {
             *fLog << err << "ERROR - Number of files found (" << n0 << ") doesn't match number of files in sequence (" << n1 << ")" << endl;
             if (fLog->GetDebugLevel()>4)
@@ -330 +293 @@
             }
             return kFALSE;
-          }
-      }
-    
+        }
+    }
+
     // Read File
     MCalibrationIntensityChargeCam      ichcam;
@@ -464 +427 @@
     case kIsUseRootData: read = &readreal; break;
     }
-    read->AddFiles(fSequence.IsValid() ? iter : *fRuns);
+    read->AddFiles(iter);
 
     const TString fname(Form("%s{s/_D_/_Y_}{s/.raw$/.root}", fPathOut.Data()));

trunk/MagicSoft/Mars/mjobs/MJCalibrateSignal.h

@@ -34 +34 @@
     Bool_t ReadExtractorCosmics(MExtractor* &ext1) const;
 
-    const char*  GetInputFileName() const;
-    const char*  GetInputFile() const;
+    const char* GetInputFile() const;
     
 public:
@@ -42 +41 @@
     ~MJCalibrateSignal();
 
-    Bool_t ProcessFile(MPedestalCam &camab, MPedestalCam &cam1, MPedestalCam &cam2);
+    Bool_t Process(MPedestalCam &camab, MPedestalCam &cam1, MPedestalCam &cam2);
 
     void SetInterlaced     ( const Bool_t b=kTRUE )  { fIsInterlaced      = b; }

trunk/MagicSoft/Mars/mjobs/MJCalibration.cc

@@ -187 +187 @@
 // Default constructor. 
 //
-// - Sets fRuns to 0, fExtractor to NULL, fTimeExtractor to NULL, fColor to kNONE, 
+// - fExtractor to NULL, fTimeExtractor to NULL, fColor to kNONE,
 //   fDisplay to kNormalDisplay, kRelTimes to kFALSE, kataCheck to kFALSE, kDebug to kFALSE,
 //   kIntensity to kFALSE
@@ -315 +315 @@
     TString title = fDisplay->GetTitle();
     title += "--  Calibration ";
-    title += fSequence.IsValid() ? Form("calib%08d", fSequence.GetSequence()) : (const char*)fRuns->GetRunsAsString();
+    title += Form("calib%08d", fSequence.GetSequence());
     title += "  --";
     fDisplay->SetTitle(title);
@@ -1377 +1377 @@
 // Retrieve the output file written by WriteResult()
 // 
-const char* MJCalibration::GetOutputFile() const
-{
-    const TString name(GetOutputFileName());
-    if (name.IsNull())
-        return "";
-
-    return Form("%s/%s", fPathOut.Data(), name.Data());
-}
-
 const char* MJCalibration::GetOutputFileName() const
 {
-
-  if (fSequence.IsValid())
-    return Form("calib%08d.root", fSequence.GetSequence());
-  if (!fRuns)
-    return "";
-  
-  return Form("%s-F1.root", (const char*)fRuns->GetRunsAsFileName());
+    return Form("%s/calib%08d.root", fPathOut.Data(), fSequence.GetSequence());
 }
 
@@ -1597 +1582 @@
 }
 
-// --------------------------------------------------------------------------
-//
-// Call the ProcessFile(MPedestalCam)
-// 
-Bool_t MJCalibration::Process(MPedestalCam &pedcam)
-{
-    if (!ReadCalibrationCam())
-        return ProcessFile(pedcam);
-
-    return kTRUE;
-}
-
 void MJCalibration::InitBlindPixel(MExtractBlindPixel &blindext,
                                    MHCalibrationChargeBlindCam &blindcam)
 {
 
-  Int_t run = fSequence.IsValid() ? fSequence.GetLastRun() : fRuns->GetRuns()[fRuns->GetNumRuns()-1];
+  Int_t run = fSequence.GetLastRun();
   
   //
@@ -1684 +1657 @@
 // Execute the task list and the eventloop:
 //
-// - Check if there are fRuns, otherwise return
 // - Check the colour of the files in fRuns (FindColor()), otherwise return
 // - Check for consistency between run numbers and number of files
@@ -1716 +1688 @@
 // - WriteResult()
 //
-Bool_t MJCalibration::ProcessFile(MPedestalCam &pedcam)
+Bool_t MJCalibration::Process(MPedestalCam &pedcam)
 {
     if (!fSequence.IsValid())
     {
-        if (!fRuns)
-        {
-            *fLog << err << "No Runs choosen... abort." << endl;
-            return kFALSE;
-        }
-
-        if (fRuns->GetNumRuns() != fRuns->GetNumEntries())
-        {
-            *fLog << err << "Number of files found doesn't match number of runs... abort."
-                << fRuns->GetNumRuns() << " vs. " << fRuns->GetNumEntries() << endl;
-            return kFALSE;
-        }
-    }
-
-    // --------------------------------------------------------------------------------
+          *fLog << err << "ERROR - Sequence invalid..." << endl;
+          return kFALSE;
+    }
 
     *fLog << inf;
     fLog->Separator(GetDescriptor());
-
-    *fLog << "Calculate MCalibrationCam from ";
-    if (fSequence.IsValid())
-        *fLog << "Sequence #" << fSequence.GetSequence() << endl;
-    else
-        *fLog << "Runs " << fRuns->GetRunsAsString() << endl;
-    *fLog << endl;
+    *fLog << "Calculate calibration constants from Sequence #";
+    *fLog << fSequence.GetSequence() << endl << endl;
 
     // --------------------------------------------------------------------------------
     
     if (!CheckEnv())
-      return kFALSE;
+        return kFALSE;
 
     // --------------------------------------------------------------------------------
@@ -1760 +1715 @@
 
     MDirIter iter;
-    if (fSequence.IsValid())
-    {
-        const Int_t n0 = fSequence.SetupCalRuns(iter, fPathData, IsUseRawData());
-        const Int_t n1 = fSequence.GetNumCalRuns();
-        if (n0==0)
+    const Int_t n0 = fSequence.SetupCalRuns(iter, fPathData, IsUseRawData());
+    const Int_t n1 = fSequence.GetNumCalRuns();
+    if (n0==0)
+    {
+        *fLog << err << "ERROR - No input files of sequence found!" << endl;
+        return kFALSE;
+    }
+    if (n0!=n1)
+    {
+        *fLog << err << "ERROR - Number of files found ("
+            << n0 << ") doesn't match number of files in sequence ("
+            << n1 << ")" << endl;
+        if (fLog->GetDebugLevel()>4)
         {
-            *fLog << err << "ERROR - No input files of sequence found!" << endl;
-            return kFALSE;
+            *fLog << dbg << "Files which are searched:" << endl;
+            iter.Print();
         }
-        if (n0!=n1)
-        {
-            *fLog << err << "ERROR - Number of files found (" 
-                  << n0 << ") doesn't match number of files in sequence (" 
-                  << n1 << ")" << endl;
-            if (fLog->GetDebugLevel()>4)
-            {
-                *fLog << dbg << "Files which are searched:" << endl;
-                iter.Print();
-            }
-            return kFALSE;
-        }
+        return kFALSE;
     }
 
@@ -1834 +1786 @@
     if (IsUseRawData())
     {
-        rawread.AddFiles(fSequence.IsValid() ? iter : *fRuns);
+        rawread.AddFiles(iter);
         tlist.AddToList(&rawread);
     }
@@ -1840 +1792 @@
     {
         read.DisableAutoScheme();
-        read.AddFiles(fSequence.IsValid() ? iter : *fRuns);
+        read.AddFiles(iter);
         tlist.AddToList(&read);
     }
@@ -1871 +1823 @@
     calcalc.SetOutputFile("");
     timecalc.SetOutputFile("");
-
-    if (!fSequence.IsValid())
-    {
-        calcalc.SetOutputPath(fPathOut);
-        calcalc.SetOutputFile(Form("%s-ChargeCalibStat.txt",(const char*)fRuns->GetRunsAsFileName()));
-        timecalc.SetOutputPath(fPathOut);
-        timecalc.SetOutputFile(Form("%s-ChargeCalibStat.txt",(const char*)fRuns->GetRunsAsFileName()));
-    }
 
     if (IsDebug())
@@ -2051 +1995 @@
     
     return kTRUE;
-}
-
-// --------------------------------------------------------------------------
-//
-// Read the following containers from GetOutputFile()
-// - MCalibrationChargeCam
-// - MCalibrationQECam
-// - MBadPixelsCam
-//
-Bool_t MJCalibration::ReadCalibrationCam()
-{
-  
-  if (IsNoStorage())
-    return kFALSE;
-
-  const TString fname = GetOutputFile();
-  
-  if (gSystem->AccessPathName(fname, kFileExists))
-  {
-      *fLog << err << "Input file " << fname << " doesn't exist." << endl;
-      return kFALSE;
-  }
-
-  *fLog << inf << "Reading from file: " << fname << endl;
-
-  TFile file(fname, "READ");
-  if (fCalibrationCam.Read()<=0)
-  {
-      *fLog << err << "Unable to read MCalibrationChargeCam from " << fname << endl;
-      return kFALSE;
-  }
-
-  if (fQECam.Read()<=0)
-  {
-      *fLog << err << "Unable to read MCalibrationQECam from " << fname << endl;
-      return kFALSE;
-  }
-
-
-  if (file.FindKey("MCalibrationRelTimeCam"))
-      if (fRelTimeCam.Read()<=0)
-      {
-          *fLog << err << "Unable to read MCalibrationRelTimeCam from " << fname << endl;
-          return kFALSE;
-      }
-
-  if (file.FindKey("MBadPixelsCam"))
-  {
-      MBadPixelsCam bad;
-      if (bad.Read()<=0)
-      {
-          *fLog << err << "Unable to read MBadPixelsCam from " << fname << endl;
-          return kFALSE;
-      }
-      fBadPixels.Merge(bad);
-  }
-
-  if (fDisplay /*&& !fDisplay->GetCanvas("Pedestals")*/) // FIXME!
-      fDisplay->Read();
-
-  return kTRUE;
 }
 

trunk/MagicSoft/Mars/mjobs/MJCalibration.h

@@ -161 +161 @@
   MJCalibration(const char *name=NULL, const char *title=NULL);
   
-  const char* GetOutputFile() const;
-  
   MCalibrationIntensityChargeCam  &GetIntensCalibrationCam() { return fIntensCalibCam;   }
   MCalibrationIntensityRelTimeCam &GetIntensRelTimeCam()     { return fIntensRelTimeCam; }
@@ -199 +197 @@
 
   // Precessing
-  Bool_t ReadCalibrationCam();
-  Bool_t ProcessFile(MPedestalCam &pedcam);
   Bool_t Process(MPedestalCam &pedcam);
 

trunk/MagicSoft/Mars/mjobs/MJCut.cc

@@ -354 +354 @@
 }
 
-Bool_t MJCut::ProcessFile(const MDataSet &set)
+Bool_t MJCut::Process(const MDataSet &set)
 {
     if (!set.IsValid())

trunk/MagicSoft/Mars/mjobs/MJCut.h

@@ -43 +43 @@
     ~MJCut();
 
-    Bool_t ProcessFile(const MDataSet &set);
+    Bool_t Process(const MDataSet &set);
 
     void EnableStorageOfSummary(Bool_t b=kTRUE)  { fStoreSummary = b; } // See SetNameSummary

trunk/MagicSoft/Mars/mjobs/MJPedestal.cc

@@ -148 +148 @@
 }
 
-const char* MJPedestal::GetOutputFile() const
-{
-    const TString name(GetOutputFileName());
-    if (name.IsNull())
-        return "";
-
-    return Form("%s/%s", fPathOut.Data(), name.Data());
-}
-
 const char* MJPedestal::GetOutputFileName() const
 {
-
-  if (fSequence.IsValid())
     return Form("pedest%08d.root", fSequence.GetSequence());
-  
-  if (!fRuns)
-    return "";
-  
-  return Form("%s-F0.root", (const char*)fRuns->GetRunsAsFileName());
-}
-
-//---------------------------------------------------------------------------------
-//
-// Try to read an existing MPedestalCam from a previously created output file. 
-// If found, also an MBadPixelsCam and the corresponding display is read.
-//
-// In case of Storage type "kNoStorage" or if the file is not found or the 
-// MPedestalCam cannot be read, return kFALSE, otherwise kTRUE.
-//
-Bool_t MJPedestal::ReadPedestalCam()
-{
-    const TString fname = GetOutputFile();
-
-    const Bool_t fileexist = !gSystem->AccessPathName(fname, kFileExists);
-    if (!fileexist)
-      {
-        *fLog << inf << "Input file " << fname << " not found, will create pedestals" << endl;
-        return kFALSE;
-      }
-
-    *fLog << inf << "Reading pedestals from file: " << fname << endl;
-
-    TFile file(fname, "READ");
-    if (fPedestalCamIn.Read()<=0)
-    {
-        *fLog << err << "Unable to read incoming MPedestalCam from " << fname << endl;
-        return kFALSE;
-    }
-
-    if (fPedestalCamOut.Read()<=0)
-    {
-        *fLog << err << "Unable to read outgoing MPedestalCam from " << fname << endl;
-        return kFALSE;
-    }
-
-    if (file.FindKey("MBadPixelsCam"))
-    {
-        MBadPixelsCam bad;
-        if (bad.Read()<=0)
-        {
-            *fLog << err << "Unable to read MBadPixelsCam from " << fname << endl;
-            return kFALSE;
-        }
-        fBadPixels.Merge(bad);
-    }
-
-    if (fDisplay && !fDisplay->GetCanvas("Pedestals"))
-        fDisplay->Read();
-
-    return kTRUE;
 }
 
@@ -243 +176 @@
 }
 
-Bool_t MJPedestal::WriteExtractor() const
-{
-
-  const TString name  = Form("pedy%08d.root",fSequence.GetSequence());
-  const TString fname = Form("%s/%s",fPathIn.Data(),name.Data());
-
-  *fLog << inf << "Updating extractor in file: " << fname << endl;
-
-  TFile file(fname, fOverwrite?"RECREATE":"NEW");
-  if (!file.IsOpen())
-    {
-      *fLog << err << dbginf << "ERROR - Could not open file " << fname << endl;
-      *fLog << err << dbginf << "Maybe, file " << fname << " already exists, call callisto with option -f then." << endl;
-      return kFALSE;
-    }
-
-  TObjArray cont;
-  cont.Add(fExtractor);
-
-  return WriteContainer(cont);
-}
-
 //---------------------------------------------------------------------------------
 //
@@ -280 +191 @@
     TString title = fDisplay->GetTitle();
     title += "--  Pedestal ";
-    if (fSequence.IsValid())
-        title += fSequence.GetName();
-    else
-        if (fRuns)  // FIXME: What to do if an environmentfile was used?
-            title += fRuns->GetRunsAsString();
+    title += fSequence.GetName();
     title += "  --";
     fDisplay->SetTitle(title);
@@ -902 +809 @@
     cont.Add(&fBadPixels);
 
-    return WriteContainer(cont, GetOutputFileName(),fOverwrite?"RECREATE":"NEW");
+    return WriteContainer(cont, GetOutputFileName(), fOverwrite?"RECREATE":"NEW");
 }
 
 Bool_t MJPedestal::Process()
 {
-    if (!ReadPedestalCam())
-        return ProcessFile();
-
-    return kTRUE;
-}
-
-Bool_t MJPedestal::ProcessFile()
-{
     if (!fSequence.IsValid())
     {
-        if (!fRuns)
-        {
-            *fLog << err << "Neither AddRuns nor SetSequence nor SetEnv was called... abort." << endl;
-            return kFALSE;
-        }
-        if (fRuns && fRuns->GetNumRuns() != fRuns->GetNumEntries())
-        {
-            *fLog << err << "Number of files found doesn't match number of runs... abort."
-                << fRuns->GetNumRuns() << " vs. " << fRuns->GetNumEntries() << endl;
-            return kFALSE;
-        }
-    }
+          *fLog << err << "ERROR - Sequence invalid..." << endl;
+          return kFALSE;
+    }
+
+    *fLog << inf;
+    fLog->Separator(GetDescriptor());
+    *fLog << "Calculate pedestal from Sequence #";
+    *fLog << fSequence.GetSequence() << endl << endl;
 
     if (!CheckEnv())
@@ -940 +835 @@
     fLog->Separator(GetDescriptor());
     *fLog << "Calculate MPedestalCam from " << type << "-runs ";
-    if (fSequence.IsValid())
-        *fLog << fSequence.GetName() << endl;
-    else
-        if (fRuns)
-            *fLog << fRuns->GetRunsAsString() << endl;
-        else
-            *fLog << "in Resource File." << endl;
+    *fLog << fSequence.GetName() << endl;
     *fLog << endl;
 
@@ -987 +876 @@
     if (IsUseRawData())
     {
-        if (fRuns || fSequence.IsValid())
-            rawread.AddFiles(fSequence.IsValid() ? iter : *fRuns);
+        rawread.AddFiles(iter);
         tlist.AddToList(&rawread);
     }
@@ -994 +882 @@
     {
         read.DisableAutoScheme();
-        if (fRuns || fSequence.IsValid())
-          read.AddFiles(fSequence.IsValid() ? iter : *fRuns);
+        read.AddFiles(iter);
         tlist.AddToList(&read);
     }

trunk/MagicSoft/Mars/mjobs/MJPedestal.h

@@ -99 +99 @@
     const MBadPixelsCam &GetBadPixels() const { return fBadPixels;   }
 
-    const char*  GetOutputFile() const;
-
     //    const MHPedestalCam &GetPedestalHist() const { return fPedestalHist;  }
 
     const Bool_t IsUseData() const { return fExtractType == kUseData; }
 
-    Bool_t Process    ();
-    Bool_t ProcessFile();
+    Bool_t Process();
 
     void SetBadPixels(const MBadPixelsCam &bad) { bad.Copy(fBadPixels); }

trunk/MagicSoft/Mars/mmuon/MHMuonPar.cc

@@ -16 +16 @@
 !
 !
-!   Author(s): Wolfgang Wittek, 03/2003 <mailto:wittek@mppmu.mpg.de>
-!   Author(s): Thomas Bretz, 04/2003 <mailto:tbretz@astro.uni-wuerzburg.de>
 !   Author(s): Markus Meyer, 02/2005 <mailto:meyer@astro.uni-wuerzburg.de>
-!
-!   Copyright: MAGIC Software Development, 2000-2004
+!   Author(s): Thomas Bretz, 04/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
+!
+!   Copyright: MAGIC Software Development, 2000-2005
 !
 !
@@ -28 +27 @@
 //
 // MHMuonPar
+//
 // This class is a histogram class for displaying muonparameters.
 // The folowing histgrams will be plotted:
@@ -37 +37 @@
 //
 // Inputcontainer:
-// MMuonSearchPar
-// MMuonCalibPar
+//   - MGeomCam
+//   - MMuonSearchPar
+//   - MMuonCalibPar
 //
 ////////////////////////////////////////////////////////////////////////////
 #include "MHMuonPar.h"
-
-#include <math.h>
 
 #include <TH1.h>
@@ -49 +48 @@
 #include <TCanvas.h>
 #include <TProfile.h>
+
 #include "MLog.h"
 #include "MLogManip.h"
@@ -56 +56 @@
 #include "MParList.h"
 
-#include "MHillas.h"
-#include "MHMuonPar.h"
 #include "MMuonSearchPar.h"
 #include "MMuonCalibPar.h"
@@ -69 +67 @@
 // Setup histograms 
 //
-MHMuonPar::MHMuonPar(const char *name, const char *title):fGeomCam(NULL), fMuonSearchPar(NULL),
-    fMuonCalibPar(NULL)
+MHMuonPar::MHMuonPar(const char *name, const char *title) :
+    fMuonSearchPar(NULL), fMuonCalibPar(NULL)
 {
     fName  = name  ? name  : "MHMuonPar";
@@ -76 +74 @@
 
     fHistRadius.SetName("Radius");
-    fHistRadius.SetTitle("Radius");
-    fHistRadius.SetXTitle("Radius[deg]");
+    fHistRadius.SetTitle("Distribution of Radius'");
+    fHistRadius.SetXTitle("R [\\circ]");
     fHistRadius.SetYTitle("Counts");
     fHistRadius.SetDirectory(NULL);
@@ -84 +82 @@
 
     fHistArcWidth.SetName("ArcWidth");
-    fHistArcWidth.SetTitle("ArcWidth");
-    fHistArcWidth.SetXTitle("ArcWidth[deg]");
+    fHistArcWidth.SetTitle("Distribution of ArcWidth");
+    fHistArcWidth.SetXTitle("W [\\circ]");
     fHistArcWidth.SetYTitle("Counts");
     fHistArcWidth.SetDirectory(NULL);
@@ -91 +89 @@
     fHistArcWidth.SetFillStyle(4000);
 
-    fHistBroad.SetName("Ringbroadening");
-    fHistBroad.SetTitle("Ringbroadening");
-    fHistBroad.SetXTitle("Radius[deg]");
-    fHistBroad.SetYTitle("ArcWidth/Radius");
+    fHistBroad.SetName("RingBroadening");
+    fHistBroad.SetTitle("Profile of ArcWidth/Radius versus Radius");
+    fHistBroad.SetXTitle("R [\\circ]");
+    fHistBroad.SetYTitle("W/R [1]");
     fHistBroad.SetDirectory(NULL);
     fHistBroad.UseCurrentStyle();
@@ -100 +98 @@
 
     fHistSize.SetName("SizeVsRadius");
-    fHistSize.SetXTitle("Radius");
-    fHistSize.SetYTitle("MuonSize");
+    fHistSize.SetTitle("Profile of Muon Size vs. Radius");
+    fHistSize.SetXTitle("Rs [[\circ]");
+    fHistSize.SetYTitle("S [phe]");
     fHistSize.SetDirectory(NULL);
     fHistSize.UseCurrentStyle();
@@ -119 +118 @@
     bins.SetEdges(20, 0.5, 1.5);
     bins.Apply(fHistSize);
-
 }
 
@@ -129 +127 @@
 Bool_t MHMuonPar::SetupFill(const MParList *plist)
 {
-    fGeomCam = (MGeomCam*)plist->FindObject("MGeomCam");
-    if (!fGeomCam)
+    MGeomCam *geom = (MGeomCam*)plist->FindObject("MGeomCam");
+    if (!geom)
     {
         *fLog << warn << "MGeomCam not found... abort." << endl;
         return kFALSE;
     }
+    fMm2Deg = geom->GetConvMm2Deg();
+
     fMuonSearchPar = (MMuonSearchPar*)plist->FindObject("MMuonSearchPar");
     if (!fMuonSearchPar)
@@ -148 +148 @@
     }
 
-    fMm2Deg = fGeomCam->GetConvMm2Deg();
-
-    ApplyBinning(*plist, "Radius", &fHistRadius);
-
-    ApplyBinning(*plist, "ArcWidth",  &fHistArcWidth);
-
-    ApplyBinning(*plist, "Ringbroadening",  &fHistBroad);
-
-    ApplyBinning(*plist, "SizeVsRadius",  &fHistSize);
+    ApplyBinning(*plist, "Radius",          &fHistRadius);
+    ApplyBinning(*plist, "ArcWidth",        &fHistArcWidth);
+    ApplyBinning(*plist, "RingBroadening",  &fHistBroad);
+    ApplyBinning(*plist, "SizeVsRadius",    &fHistSize);
 
     return kTRUE;
@@ -211 +206 @@
     fHistBroad.Draw();
 }
-
+/*
 TH1 *MHMuonPar::GetHistByName(const TString name)
 {
@@ -220 +215 @@
     return NULL;
 }
-
+*/

trunk/MagicSoft/Mars/mmuon/MHMuonPar.h

@@ -12 +12 @@
 #endif
 
-class MHillas;
 class MMuonSearchPar;
 class MMuonCalibPar;
@@ -20 +19 @@
 {
 private:
-    TH1F fHistRadius;  //
+    TH1F     fHistRadius;    //
+    TH1F     fHistArcWidth;  //
 
-    TH1F fHistArcWidth;   //
-
-    TProfile fHistBroad;  // Area of used pixels
-
+    TProfile fHistBroad;     // Area of used pixels
     TProfile fHistSize;      // [ratio] concentration ratio: sum of the two highest pixels / fSize
 
-    MGeomCam *fGeomCam; //! Camera geometry for plots (for the moment this is a feature for a loop only!)
-
-    MMuonSearchPar *fMuonSearchPar;//!
-
-    MMuonCalibPar *fMuonCalibPar;//!
+    MMuonSearchPar *fMuonSearchPar; //!
+    MMuonCalibPar  *fMuonCalibPar;  //!
 
     Float_t fMm2Deg;
@@ -42 +36 @@
     Bool_t Fill(const MParContainer *par, const Stat_t w=1);
 
-    TH1 *GetHistByName(const TString name);
+    //TH1 *GetHistByName(const TString name);
 
-    TH1F &GetHistRadius()  { return fHistRadius; }
-
-    TH1F &GetHistArcWidth()   { return fHistArcWidth; }
-
-    TProfile &GetHistBroad()  { return fHistBroad; }
-
-    TProfile &GetHistSize()      { return fHistSize; }
+    const TH1F&     GetHistRadius() const    { return fHistRadius; }
+    const TH1F&     GetHistArcWidth() const  { return fHistArcWidth; }
+    const TProfile& GetHistBroad() const     { return fHistBroad; }
+    const TProfile& GetHistSize() const      { return fHistSize; }
 
     void Draw(Option_t *opt="");

trunk/MagicSoft/Mars/mmuon/MMuonCalibPar.cc

@@ -17 +17 @@
 !
 !   Author(s): Keiichi Mase 10/2004 <mailto:mase@mppmu.mpg.de>
-!              Markus Meyer 10/2004 <mailto:meyer@astro.uni-wuerzburg.de>
+!   Author(s): Markus Meyer 10/2004 <mailto:meyer@astro.uni-wuerzburg.de>
 !
-!   Copyright: MAGIC Software Development, 2000-2004
+!   Copyright: MAGIC Software Development, 2000-2005
 !
 !
@@ -30 +30 @@
 // Storage Container for muon
 //
-//  This class holds some information for a calibraion using muons. Muons 
+//  This class holds some information for a calibration using muons. Muons
 // are identified by using the class of the MMuonSearchParCalc. You can fill 
 // these information by using the MMuonCalibParCalc. See also these class
 // manuals.
 //
-//  
-//  Input Containers:
-//   [MGeomCam]
-//   [MSignalCam]
-//   [MMuonSearchPar]
-//
 /////////////////////////////////////////////////////////////////////////////
 #include "MMuonCalibPar.h"
 
-#include <fstream>
-
 #include "MLog.h"
 #include "MLogManip.h"
-#include "MSignalCam.h"
-#include "MSignalPix.h"
-#include "MMuonSearchPar.h"
-#include "MBinning.h"
 
 using namespace std;
@@ -66 +54 @@
     fTitle = title ? title : "Muon calibration parameters";
 
-    fHistPhi   = new TH1F;
-    fHistWidth = new TH1F;
-
-    fHistPhi->SetName("HistPhi");
-    fHistPhi->SetTitle("HistPhi");
-    fHistPhi->SetXTitle("phi [deg.]");
-    fHistPhi->SetYTitle("sum of ADC");
-    fHistPhi->SetDirectory(NULL);
-    fHistPhi->SetFillStyle(4000);
-    fHistPhi->UseCurrentStyle();
-
-    fHistWidth->SetName("HistWidth");
-    fHistWidth->SetTitle("HistWidth");
-    fHistWidth->SetXTitle("distance from the ring center [deg.]");
-    fHistWidth->SetYTitle("sum of ADC");
-    fHistWidth->SetDirectory(NULL);
-    fHistWidth->SetFillStyle(4000);
-    fHistWidth->UseCurrentStyle();
-
-}
-
-// --------------------------------------------------------------------------
-//
-MMuonCalibPar::~MMuonCalibPar()
-{
-  delete fHistPhi;
-  delete fHistWidth;
+    Reset();
 }
 
@@ -99 +61 @@
 void MMuonCalibPar::Reset()
 {
-  fArcLength   = -1.;
-  fArcPhi      =  0.;
-  fArcWidth    = -1.;
-  fChiArcPhi   = -1.;
-  fChiArcWidth = -1.;
-  fMuonSize    =  0.;
-  fEstImpact   = -1.;
-  fUseUnmap    = kFALSE;
-  fPeakPhi     =  0.;
-
-  fHistPhi->Reset();
-  fHistWidth->Reset();
+    fArcLength   = -1.;
+    fArcPhi      =  0.;
+    fArcWidth    = -1.;
+    fChiArcPhi   = -1.;
+    fChiArcWidth = -1.;
+    fMuonSize    =  0.;
+    fEstImpact   = -1.;
+    fPeakPhi     =  0.;
 }
 
@@ -116 +74 @@
 {
     *fLog << all;
-    *fLog << "Muon Parameters (" << GetName() << ")"    << endl;
-    *fLog << " - Arc Length    [deg.]  = " << fArcLength     << endl;
-    *fLog << " - Arc Phi       [deg.]  = " << fArcPhi     << endl;
-    *fLog << " - Arc Width     [deg.]  = " << fArcWidth     << endl;
-    *fLog << " - Chi Arc Phi   [x2/ndf]= " << fChiArcPhi  << endl;
-    *fLog << " - Chi Arc Width [x2/ndf]= " << fChiArcWidth  << endl;
-    *fLog << " - Est. I. P.    [m]     = " << fEstImpact  << endl;
-    *fLog << " - Size of muon          = " << fMuonSize   << endl;
-    *fLog << " - Peak Phi      [deg.]  = " << fPeakPhi    << endl;
-    *fLog << " - UseUnmap              = " << fUseUnmap   << endl;
+    *fLog << "Muon Parameters (" << GetName() << ")"       << endl;
+    *fLog << " - Arc Length    [deg]   = " << fArcLength   << endl;
+    *fLog << " - Arc Phi       [deg]   = " << fArcPhi      << endl;
+    *fLog << " - Arc Width     [deg]   = " << fArcWidth    << endl;
+    *fLog << " - Chi Arc Phi   [x2/ndf]= " << fChiArcPhi   << endl;
+    *fLog << " - Chi Arc Width [x2/ndf]= " << fChiArcWidth << endl;
+    *fLog << " - Est. I. P.    [m]     = " << fEstImpact   << endl;
+    *fLog << " - Size of muon  [phe]   = " << fMuonSize    << endl;
+    *fLog << " - Peak Phi      [deg]   = " << fPeakPhi     << endl;
 }
 

trunk/MagicSoft/Mars/mmuon/MMuonCalibPar.h

@@ -16 +16 @@
 {
 private:
-  Float_t fArcLength;     // An arc length of a muon along the arc [deg.]
-  Float_t fArcPhi;        // A opening angle of a muon arc [deg.]
-  Float_t fArcWidth;      // A width of a muon [deg.] (1 sigma of gaussian fit)
-  Float_t fChiArcPhi;     // A chisquare value of the cosine fit for arc phi
-  Float_t fChiArcWidth;   // A chisquare value of the cosine fit for arc wid
-  Float_t fMuonSize;      // A SIZE of muon which is defined as a SIZE around the estimated circle
-  Float_t fEstImpact;     // An estimated impact parameter from the photon distribution along the arc image
-  Bool_t  fUseUnmap;      // This is a flag to know the Unmapped pixels are used. Refer to the class of MImgCleanStd
-  Float_t fPeakPhi;       // The angle which indicates the peak position in the estimated circle
+    Float_t fArcLength;     // An arc length of a muon along the arc [deg.]
+    Float_t fArcPhi;        // A opening angle of a muon arc [deg.]
+    Float_t fArcWidth;      // A width of a muon [deg.] (1 sigma of gaussian fit)
+    Float_t fChiArcPhi;     // A chisquare value of the cosine fit for arc phi
+    Float_t fChiArcWidth;   // A chisquare value of the cosine fit for arc wid
+    Float_t fMuonSize;      // A SIZE of muon which is defined as a SIZE around the estimated circle
+    Float_t fEstImpact;     // An estimated impact parameter from the photon distribution along the arc image
+    //Bool_t  fUseUnmap;      // This is a flag to know the Unmapped pixels are used. Refer to the class of MImgCleanStd
+    Float_t fPeakPhi;       // The angle which indicates the peak position in the estimated circle
 
 public:
-  MMuonCalibPar(const char *name=NULL, const char *title=NULL);
-  ~MMuonCalibPar();
-  
-  TH1F *fHistPhi;   // Histogram of photon distribution along the arc.
-  TH1F *fHistWidth; // Histogram of radial photon distribution of the arc.
-  
-  void Reset();
-  
-  Float_t GetArcLength()      const { return fArcLength; }
-  Float_t GetArcPhi()         const { return fArcPhi; }
-  Float_t GetArcWidth()       const { return fArcWidth; }
-  Float_t GetChiArcPhi()      const { return fChiArcPhi; }
-  Float_t GetChiArcWidth()    const { return fChiArcWidth; }
-  Float_t GetMuonSize()       const { return fMuonSize; }
-  Float_t GetEstImpact()      const { return fEstImpact; }
-  Bool_t  IsUseUnmap()        const { return fUseUnmap; }
-  Float_t GetPeakPhi()        const { return fPeakPhi; }
-  TH1F    *GetHistPhi()       { return fHistPhi; }
-  TH1F    *GetHistWidth()     { return fHistWidth; }
-  
-  void    SetArcLength(Float_t len)       { fArcLength = len; }
-  void    SetArcPhi(Float_t phi)          { fArcPhi = phi; }
-  void    SetArcWidth(Float_t wid)        { fArcWidth = wid; }
-  void    SetChiArcPhi(Float_t chi)       { fChiArcPhi = chi; }
-  void    SetChiArcWidth(Float_t chi)     { fChiArcWidth = chi; }
-  void    SetMuonSize(Float_t size)       { fMuonSize = size; }
-  void    SetEstImpact(Float_t impact)    { fEstImpact = impact; }
-  void    SetUseUnmap()                   { fUseUnmap = kTRUE; }
-  void    SetPeakPhi(Float_t phi)         { fPeakPhi = phi; }
-  
-  void    Print(Option_t *opt=NULL) const;
+    MMuonCalibPar(const char *name=NULL, const char *title=NULL);
+    //~MMuonCalibPar();
 
-  ClassDef(MMuonCalibPar, 1) // Container to hold muon calibration parameters
+    void Reset();
+
+    Float_t GetArcLength()      const { return fArcLength; }
+    Float_t GetArcPhi()         const { return fArcPhi; }
+    Float_t GetArcWidth()       const { return fArcWidth; }
+    Float_t GetChiArcPhi()      const { return fChiArcPhi; }
+    Float_t GetChiArcWidth()    const { return fChiArcWidth; }
+    Float_t GetMuonSize()       const { return fMuonSize; }
+    Float_t GetEstImpact()      const { return fEstImpact; }
+    //Bool_t  IsUseUnmap()        const { return fUseUnmap; }
+    Float_t GetPeakPhi()        const { return fPeakPhi; }
+    //  TH1F    *GetHistPhi()       { return fHistPhi; }
+    //  TH1F    *GetHistWidth()     { return fHistWidth; }
+
+    void    SetArcLength(Float_t len)       { fArcLength = len; }
+    void    SetArcPhi(Float_t phi)          { fArcPhi = phi; }
+    void    SetArcWidth(Float_t wid)        { fArcWidth = wid; }
+    void    SetChiArcPhi(Float_t chi)       { fChiArcPhi = chi; }
+    void    SetChiArcWidth(Float_t chi)     { fChiArcWidth = chi; }
+    void    SetMuonSize(Float_t size)       { fMuonSize = size; }
+    void    SetEstImpact(Float_t impact)    { fEstImpact = impact; }
+    //void    SetUseUnmap()                   { fUseUnmap = kTRUE; }
+    void    SetPeakPhi(Float_t phi)         { fPeakPhi = phi; }
+
+    void    Print(Option_t *opt=NULL) const;
+
+    ClassDef(MMuonCalibPar, 1) // Container to hold muon calibration parameters
 };
     

trunk/MagicSoft/Mars/mmuon/MMuonCalibParCalc.cc

@@ -67 +67 @@
 //
 // 
-//  For the faster computation, by default, the calculation of impact
-// parameter is suppressed. If you want to calculate the impact parameter
-// from the muon image, you can use the function of EnableImpactCalc(),
-// namely;
-//
-//   mucalibcalc.EnableImpactCalc();.
+// // For the faster computation, by default, the calculation of impact
+// // parameter is suppressed. If you want to calculate the impact parameter
+// // from the muon image, you can use the function of EnableImpactCalc(),
+// // namely;
+// // 
+// //   mucalibcalc.EnableImpactCalc();.
 //
 //  In addition, for the faster comutation, pre cuts to select the candidates
@@ -98 +98 @@
 //
 //  Input Containers:
-//   [MGeomCam]
-//   [MSignalCam]
-//   [MMuonSearchPar]
+//   MGeomCam
+//   MSignalCam
+//   MMuonSearchPar
 //
 //  Output Containers:
-//   [MMuonCalibPar]
+//   MMuonCalibPar
 //
 //////////////////////////////////////////////////////////////////////////////
 
 #include "MMuonCalibParCalc.h"
-
-#include <fstream>
 
 #include <TH1.h>
@@ -115 +113 @@
 #include <TMinuit.h>
 
+#include "MLog.h"
+#include "MLogManip.h"
+
 #include "MParList.h"
 
 #include "MGeomCam.h"
 #include "MGeomPix.h"
-#include "MSrcPosCam.h"
+
 #include "MSignalCam.h"
-#include "MMuonSearchPar.h"
+
 #include "MMuonCalibPar.h"
 #include "MMuonSetup.h"
-#include "MLog.h"
-#include "MLogManip.h"
-#include "MBinning.h"
+#include "MMuonSearchPar.h"
+#include "MHSingleMuon.h"
 
 using namespace std;
@@ -140 +140 @@
 //
 MMuonCalibParCalc::MMuonCalibParCalc(const char *name, const char *title)
+//    : fEnableImpactCalc(kFALSE)
 {
     fName  = name  ? name  : gsDefName.Data();
     fTitle = title ? title : gsDefTitle.Data();
-
-    fEnableImpactCalc = kFALSE; // By default the calculation of impact parameter is skipped.
 }
 
@@ -151 +150 @@
 Int_t MMuonCalibParCalc::PreProcess(MParList *pList)
 {
-  fSignalCam = (MSignalCam*)pList->FindObject("MSignalCam");
-  if (!fSignalCam)
+    fGeomCam = (MGeomCam*)pList->FindObject("MGeomCam");
+    if (!fGeomCam)
     {
-      *fLog << dbginf << "MSignalCam not found... aborting." << endl;
-      return kFALSE;
+        *fLog << err << "MGeomCam not found... abort." << endl;
+        return kFALSE;
     }
-  
-  fGeomCam = (MGeomCam*)pList->FindObject("MGeomCam");
-  if (!fGeomCam)
+
+    fHist = (MHSingleMuon*)pList->FindObject("MHSingleMuon");
+    if (!fHist)
     {
-      *fLog << dbginf << "MGeomCam (Camera Geometry) missing in Parameter List... aborting." << endl;
-      return kFALSE;
+        *fLog << err << "MHSingleMuon not found... abort." << endl;
+        return kFALSE;
     }
-  
-  fMuonCalibPar = (MMuonCalibPar*)pList->FindCreateObj("MMuonCalibPar", "MMuonCalibPar");
-  if (!fMuonCalibPar)
-      return kFALSE;
-  
-  fMuonSearchPar = (MMuonSearchPar*)pList->FindCreateObj("MMuonSearchPar", "MMuonSearchPar");
-  if (!fMuonSearchPar)
-      return kFALSE;
-
-  fMuonSetup = (MMuonSetup*)pList->FindCreateObj("MMuonSetup", "MMuonSetup");
-  if (!fMuonSetup)
-      return kFALSE;
-  
-  return kTRUE;
+
+    fMuonSetup = (MMuonSetup*)pList->FindObject("MMuonSetup");
+    if (!fMuonSetup)
+    {
+        *fLog << err << "MMuonSetup not found... abort." << endl;
+        return kFALSE;
+    }
+
+    fMuonCalibPar = (MMuonCalibPar*)pList->FindCreateObj("MMuonCalibPar");
+    if (!fMuonCalibPar)
+        return kFALSE;
+
+    fMuonSearchPar = (MMuonSearchPar*)pList->FindCreateObj("MMuonSearchPar");
+    if (!fMuonSearchPar)
+        return kFALSE;
+
+    return kTRUE;
 }
-
-// --------------------------------------------------------------------------
-// 
-//  This function fill the histograms in order to get muon parameters.
-// For the evaluation of the Arc Width, we use only the signals in the inner 
-// part. 
-// 
-void MMuonCalibParCalc::FillHist()
-{
-  Float_t MuonSize = 0.;
-
-  Int_t binnumphi = fMuonSetup->fArcPhiBinNum;
-  Int_t binnumwid = fMuonSetup->fArcWidthBinNum;
-
-  // preparation for a histgram
-  MBinning binsphi;
-  binsphi.SetEdges(binnumphi, 
-                   fMuonSetup->fArcPhiHistStartVal,
-                   fMuonSetup->fArcPhiHistEndVal);
-  binsphi.Apply(*(fMuonCalibPar->fHistPhi));
-
-  MBinning binswid;
-  binswid.SetEdges(binnumwid, 
-                   fMuonSetup->fArcWidthHistStartVal,
-                   fMuonSetup->fArcWidthHistEndVal);
-  binswid.Apply(*(fMuonCalibPar->fHistWidth));
-
-  const Int_t entries = (*fSignalCam).GetNumPixels();
-
-  // the position of the center of a muon ring
-  const Float_t cenx = (*fMuonSearchPar).GetCenterX();
-  const Float_t ceny = (*fMuonSearchPar).GetCenterY();
-  
-  for (Int_t i=0; i<entries; i++ )
-    {
-      MSignalPix &pix = (*fSignalCam)[i];
-      
-      const MGeomPix &gpix = (*fGeomCam)[i/*pix.GetPixId()*/];
-      
-      const Float_t dx = gpix.GetX() - cenx;
-      const Float_t dy = gpix.GetY() - ceny;
-      
-     const Float_t dist = TMath::Sqrt(dx*dx+dy*dy);
-      
-      Float_t ang = TMath::ACos(dx/dist);
-      if(dy>0)
-        ang *= -1.0;
-      
-      // if the signal is not near the estimated circle, it is ignored.
-      if(dist < (*fMuonSearchPar).GetRadius() + fMuonSetup->GetMargin()
-         && dist > (*fMuonSearchPar).GetRadius() - fMuonSetup->GetMargin())
-        {
-          // check whether ummapped pixel is used or not.
-          // if it is so, ingnore the pixel information since the pixels totally deteriorate the muon information.
-          if(pix.IsPixelUnmapped())
-            {
-              fMuonCalibPar->SetUseUnmap();
-              continue;
-            }
-          fMuonCalibPar->fHistPhi->Fill(ang*kRad2Deg, pix.GetNumPhotons());
-          MuonSize += pix.GetNumPhotons();
-        }
-
-      // use only the inner pixles. This is geometry dependent. This has to
-      // be fixed!
-      if(i>397)
-        continue;  
-      
-      fMuonCalibPar->fHistWidth
-        ->Fill(dist*(*fGeomCam).GetConvMm2Deg(), pix.GetNumPhotons());
-    }
-
-  fMuonCalibPar->SetMuonSize(MuonSize);
-
-  // error estimation (temporaly)
-  //  The error is estimated from the signal. In order to do so, we have to
-  // once convert the signal from ADC to photo-electron. Then we can get
-  // the fluctuation such as F-factor*sqrt(phe). 
-  //  Up to now, the error of pedestal is not taken into accout. This is not 
-  // of course correct. We will include this soon.
-    Double_t ADC2PhEl = 0.14;
-    Double_t Ffactor = 1.4;
-    for(Int_t i=0; i<binnumphi+1; i++)
-      {
-        Float_t rougherr  = TMath::Sqrt(TMath::Abs(fMuonCalibPar->fHistPhi->GetBinContent(i))*ADC2PhEl)/ADC2PhEl*Ffactor;
-        {
-          fMuonCalibPar->fHistPhi->SetBinError(i, rougherr);
-        }
-      }
-    for(Int_t i=0; i<binnumwid+1; i++)
-      {
-        Float_t rougherr = TMath::Sqrt(TMath::Abs(fMuonCalibPar->fHistWidth->GetBinContent(i))*ADC2PhEl)/ADC2PhEl*Ffactor;
-        {
-          fMuonCalibPar->fHistWidth->SetBinError(i, rougherr);
-        }
-      }
-}
-
-// --------------------------------------------------------------------------
-//
-//  Photon distribution along the estimated circle is fitted with theoritical
-// function in order to get some more information such as Arc Phi and Arc 
-// Length.
-//
-void MMuonCalibParCalc::CalcPhi()
-{
-  Float_t thres = fMuonSetup->GetArcPhiThres();
-  Float_t startval = fMuonSetup->fArcPhiHistStartVal;
-  Float_t endval = fMuonSetup->fArcPhiHistEndVal;
-  Int_t   binnum = fMuonSetup->fArcPhiBinNum;
-
-  Float_t convbin2val = (endval-startval)/(Float_t)binnum;
-
-    // adjust the peak to 0
-    Float_t maxval = 0.;
-    Int_t   maxbin = 0;
-    maxval = fMuonCalibPar->fHistPhi->GetMaximum();
-    maxbin = fMuonCalibPar->fHistPhi->GetMaximumBin();
-    fMuonCalibPar->SetPeakPhi(180.-(Float_t)(maxbin-1.)*convbin2val); 
-    TArrayD tmp;
-    tmp.Set(binnum+1);
-    for(Int_t i=1; i<binnum+1; i++)
-      {
-        tmp[i] = fMuonCalibPar->fHistPhi->GetBinContent(i);
-      }
-    for(Int_t i=1; i<binnum+1; i++)
-      {
-        Int_t id;
-        id = i + (maxbin-(Int_t)((Float_t)binnum/2.)-1);
-        if(id>binnum)
-          {
-            id-=(binnum);
-          }
-        if(id<=0)
-          {
-            id+=(binnum);
-          }
-        fMuonCalibPar->fHistPhi->SetBinContent(i,tmp[id]);
-      }
-    maxbin = (Int_t)((Float_t)binnum/2.)+1;
-
-  // Determination of fitting region
-  // The threshold is fixed with 100 [photons or ADC] in a bin. Therefore,
-  // if you change the bin number, YOU HAVE TO CHANGE THIS VALUE!!!
-  Float_t startfitval = 0.;
-  Float_t endfitval   = 0.;
-  Bool_t  IsInMaxim = kFALSE;
-  Int_t effbinnum = 0;
-  for(Int_t i=1; i<binnum+1; i++)
-    {
-      Float_t content = fMuonCalibPar->fHistPhi->GetBinContent(i);
-      Float_t content_pre = fMuonCalibPar->fHistPhi->GetBinContent(i-1);
-      
-      if(content > thres && content_pre < thres)
-        {
-          startfitval = (Float_t)(i-1)*convbin2val+startval;
-        }
-      if(i==maxbin)
-        IsInMaxim = kTRUE;
-      
-      if(content < thres && IsInMaxim == kTRUE)
-        {
-          endfitval = (Float_t)(i-1)*convbin2val+startval;
-          break;
-        }
-      endfitval = endval;
-    }
-  
-  effbinnum = (Int_t)((endfitval-startfitval)/convbin2val);
-
-  fMuonCalibPar->SetArcPhi(endfitval-startfitval);
-
-  fMuonCalibPar->SetArcLength( fMuonCalibPar->GetArcPhi()*TMath::DegToRad()*(*fMuonSearchPar).GetRadius()*(*fGeomCam).GetConvMm2Deg());
-  
-  if(fEnableImpactCalc)
-      CalcImpact(effbinnum, startfitval, endfitval);
-}
-
-// --------------------------------------------------------------------------
-//
-//  An impact parameter is calculated by fitting the histogram of photon
-// distribution along the circle with a theoritical model. 
-// (See G. Vacanti et. al., Astroparticle Physics 2, 1994, 1-11. 
-// The function (6) is used here.) 
-//
-//  By default this calculation is suppressed because this calculation is 
-// very time consuming. If you want to calculate an impact parameter,
-// you can call the function of EnableImpactCalc().
-//
-void MMuonCalibParCalc::CalcImpact
-(Int_t effbinnum, Float_t startfitval, Float_t endfitval)
-{
-  // Fit the distribution with Vacanti function. The function is different
-  // for the impact parameter of inside or outside of our reflector. 
-  // Then two different functions are applied to the photon distribution,
-  // and the one which give us smaller chisquare value is taken as a 
-  // proper one.
-  Double_t val1,err1,val2,err2;
-  // impact parameter inside mirror radius (8.5m)
-  TString func1;
-  Float_t tmpval = (*fMuonSearchPar).GetRadius()*(*fGeomCam).GetConvMm2Deg()*TMath::DegToRad();
-  tmpval = sin(2.*tmpval)*8.5;
-  func1 += "[0]*";
-  func1 += tmpval;
-  func1 += "*(sqrt(1.-([1]/8.5)**2*sin((x-[2])*3.1415926/180.)**2)+([1]/8.5)*cos((x-[2])*3.1415926/180.))";
-  
-  TF1 f1("f1",func1,startfitval,endfitval);
-  f1.SetParameters(2000,3,0);
-  f1.SetParLimits(1,0,8.5);
-  f1.SetParLimits(2,-180.,180.);
-  
-  fMuonCalibPar->fHistPhi->Fit("f1","RQEM");
-  
-  Float_t chi1 = -1;
-  Float_t chi2 = -1;
-  if(effbinnum>3)
-    chi1 = f1.GetChisquare()/((Float_t)(effbinnum-3));
-  
-  gMinuit->GetParameter(1,val1,err1);  // get the estimated IP
-  Float_t estip1 = val1;
-  
-  // impact parameter beyond mirror area (8.5m)
-  TString func2;
-  Float_t tmpval2 = (*fMuonSearchPar).GetRadius()*(*fGeomCam).GetConvMm2Deg()*TMath::DegToRad();
-  tmpval2 = sin(2.*tmpval2)*8.5*2.;
-  func2 += "[0]*";
-  func2 += tmpval2;
-  func2 += "*sqrt(1.-(([1]/8.5)*sin((x-[2])*3.1415926/180.))**2)";
-  
-  TF1 f2("f2",func2,startfitval,endfitval);
-  f2.SetParameters(2000,20,0);
-  f2.SetParLimits(1,8.5,300.);
-  f2.SetParLimits(2,-180.,180.);
-  
-  fMuonCalibPar->fHistPhi->Fit("f2","RQEM+");
-  
-  if(effbinnum>3)
-    chi2 = f2.GetChisquare()/((Float_t)(effbinnum-3));
-  
-  gMinuit->GetParameter(1,val2,err2);  // get the estimated IP
-  Float_t estip2 = val2;
-  
-  if(effbinnum<=3)
-    {
-      fMuonCalibPar->SetEstImpact(-1.);
-    }
-  if(chi2 > chi1)
-    {
-      fMuonCalibPar->SetEstImpact(estip1);
-      fMuonCalibPar->SetChiArcPhi(chi1);
-    }
-  else
-    {
-      fMuonCalibPar->SetEstImpact(estip2);
-      fMuonCalibPar->SetChiArcPhi(chi2);
-    }
-}
-
-// --------------------------------------------------------------------------
-//
-//  Photon distribution of distance from the center of estimated ring is 
-// fitted in order to get some more information such as ARC WIDTH which 
-// can represent to the PSF of our reflector.
-//
-Float_t MMuonCalibParCalc::CalcWidth()
-{
-  Float_t startval = fMuonSetup->fArcWidthHistStartVal;
-  Float_t endval = fMuonSetup->fArcWidthHistEndVal;
-  Int_t   binnum = fMuonSetup->fArcWidthBinNum;
-  Float_t thres = fMuonSetup->GetArcWidthThres();
-
-  Float_t convbin2val = (endval - startval)
-    /(Float_t)binnum;
-
-    // determination of fitting region
-    Int_t maxbin = fMuonCalibPar->fHistWidth->GetMaximumBin();
-    Float_t startfitval = 0.;
-    Float_t endfitval   = 0.;
-    Bool_t   IsInMaxim = kFALSE;
-    Int_t    effbinnum = 0;
-    for(Int_t i=1; i<binnum+1; i++)
-      {
-        Float_t content = fMuonCalibPar->fHistWidth->GetBinContent(i);
-        Float_t content_pre = fMuonCalibPar->fHistWidth->GetBinContent(i-1);
-
-        if(content > thres)
-          effbinnum++;
-
-        if(content > thres && content_pre < thres)
-          {
-            startfitval = (Float_t)(i-4)*convbin2val + startval;
-            if(startfitval<0.) startfitval = 0.;
-          }
-        if(i==maxbin)
-          IsInMaxim = kTRUE;
-
-        if(content < thres && IsInMaxim == kTRUE)
-          {
-            endfitval = (Float_t)(i+2)*convbin2val + startval;
-            if(endfitval>180.) endfitval = 180.;
-            break;
-          }
-        endfitval = endval;
-      }
-    effbinnum = (Int_t)((endfitval-startfitval)/convbin2val);
-
-    TF1 f1("f1","gaus",startfitval,endfitval);
-
-    // options : N  do not store the function, do not draw
-    //           I  use integral of function in bin rather than value at bin center
-    //           R  use the range specified in the function range
-    //           Q  quiet mode
-    fMuonCalibPar->fHistWidth->Fit("f1","QR","",startfitval,endfitval);
-    
-    if(effbinnum>3)
-      fMuonCalibPar->SetChiArcWidth(f1.GetChisquare()/((Float_t)(effbinnum-3)));
-
-    Double_t val,err;
-    gMinuit->GetParameter(2,val,err); // get the sigma value
-
-    return val;
-}
-
-// --------------------------------------------------------------------------
-//
-//  Calculation of muon parameters
-//
-Int_t MMuonCalibParCalc::Calc()
-{
-  // initialization
-  (*fMuonCalibPar).Reset();
-
-  // Fill signals to histograms
-  FillHist();
-
-  // Calculation of Arc Phi etc...
-  CalcPhi();
-
-  // Calculation of Arc Width etc...
-  fMuonCalibPar->SetArcWidth(CalcWidth());
-
-  if(fMuonCalibPar->GetArcPhi()>160 && fMuonSearchPar->GetRadius()>170 &&
-     fMuonSearchPar->GetRadius()<400 && fMuonSearchPar->GetDeviation()<50)
-      fMuonCalibPar->SetReadyToSave();
-  
-  return kTRUE;
-} 
 
 // -------------------------------------------------------------------------
@@ -531 +186 @@
 Int_t MMuonCalibParCalc::Process()
 {
-
-  Calc();
-
-  return kTRUE;
+    // Calculation of Arc Phi etc...
+    const Float_t thresphi   = fMuonSetup->GetThresholdArcPhi()  /fGeomCam->GetConvMm2Deg();
+    const Float_t threswidth = fMuonSetup->GetThresholdArcWidth()/fGeomCam->GetConvMm2Deg();
+
+    Double_t peakphi, arcphi;
+    Double_t width, chi;
+
+    if (!fHist->CalcPhi(thresphi, peakphi, arcphi))
+        return kTRUE;
+
+    if (!fHist->CalcWidth(threswidth, width, chi))
+        return kTRUE;
+
+    // Get Muon Size
+    fMuonCalibPar->SetMuonSize(fHist->GetHistPhi().Integral());
+
+    // Calculate Arc Length
+    fMuonCalibPar->SetPeakPhi(peakphi);
+    fMuonCalibPar->SetArcPhi(arcphi);
+
+    const Float_t conv = TMath::RadToDeg()*fGeomCam->GetConvMm2Deg();
+    fMuonCalibPar->SetArcLength(fMuonCalibPar->GetArcPhi()*fMuonSearchPar->GetRadius()*conv);
+
+    // Calculation of Arc Width etc...
+    fMuonCalibPar->SetChiArcWidth(chi);
+    fMuonCalibPar->SetArcWidth(width);
+
+    // Check if this is a 'Write-Out' candidate
+    if (fMuonCalibPar->GetArcPhi()>160    && fMuonSearchPar->GetRadius()<400 &&
+        fMuonSearchPar->GetDeviation()<50 && fMuonSearchPar->GetRadius()>170)
+        fMuonCalibPar->SetReadyToSave();
+
+    return kTRUE;
 }

trunk/MagicSoft/Mars/mmuon/MMuonCalibParCalc.h

@@ -8 +8 @@
 class MMuonSearchPar;
 class MMuonCalibPar;
-class MSrcPosCam;
 class MGeomCam;
-class MSignalCam;
 class MMuonSetup;
+class MHSingleMuon;
 
 class MMuonCalibParCalc : public MTask
 {
 private:
-    MGeomCam       *fGeomCam;
-    MSignalCam     *fSignalCam;
-    MMuonCalibPar  *fMuonCalibPar;  
-    MMuonSearchPar *fMuonSearchPar; 
-    MMuonSetup     *fMuonSetup;
+    MGeomCam       *fGeomCam;         //!
+    MMuonCalibPar  *fMuonCalibPar;    //!
+    MMuonSearchPar *fMuonSearchPar;   //!
+    MMuonSetup     *fMuonSetup;       //!
+    MHSingleMuon   *fHist;            //!
 
-    Bool_t fEnableImpactCalc; // If true, the impact calculation will be done, which consumes a lot of time.
+    //Bool_t fEnableImpactCalc; // If true, the impact calculation will be done, which consumes a lot of time.
 
-    Int_t PreProcess(MParList *plist);
-    Int_t Process();
+    Int_t   PreProcess(MParList *plist);
+    Int_t   Process();
+
+    void    FillHist();
+    void    CalcPhi();
+    void    CalcImpact(Int_t effbinnum, Float_t startfitval, Float_t endfitval);
+    Float_t CalcWidth();
 
 public:
     MMuonCalibParCalc(const char *name=NULL, const char *title=NULL);
 
-    void EnableImpactCalc()              { fEnableImpactCalc = kTRUE; }
-
-    void FillHist();
-    void CalcPhi();
-    void CalcImpact(Int_t effbinnum, Float_t startfitval, Float_t endfitval);
-    Float_t CalcWidth();
-    Int_t   Calc();
+    //void EnableImpactCalc(Bool_t b=kTRUE) { fEnableImpactCalc = b; }
 
     ClassDef(MMuonCalibParCalc, 0) // task to calculate muon parameters

trunk/MagicSoft/Mars/mmuon/MMuonSearchPar.cc

@@ -17 +17 @@
 !
 !   Author(s): Keiichi Mase 10/2004 <mailto:mase@mppmu.mpg.de>
-!              Markus Meyer 10/2004 <mailto:meyer@astro.uni-wuerzburg.de>
-!
-!   Copyright: MAGIC Software Development, 2000-2004
+!   Author(s): Markus Meyer 10/2004 <mailto:meyer@astro.uni-wuerzburg.de>
+!
+!   Copyright: MAGIC Software Development, 2000-2005
 !
 !
@@ -40 +40 @@
 // MMuonCalibPar. The information will be available by using the task of 
 // MMuonCalibParCalc.
-// 
 // 
 // --- How to search muons --- 
@@ -47 +46 @@
 // 
 // 1. A temporal center position of a circle is determined by using 
-//  the Hillas parameters. Assumed that the center position will be on the
-//  line which is perpendicular to the longitudinal image axis and the 
-//  distance from the gravity center of the image to the center position of
-//  a ring is approximately 1 deg. (corresponding to the Cherenkov angle.).
-//  Therefore, we will have two candidates of the center positions.
+//    the Hillas parameters. Assumed that the center position will be on the
+//    line which is perpendicular to the longitudinal image axis and the
+//    distance from the gravity center of the image to the center position of
+//    a ring is approximately 1 deg. (corresponding to the Cherenkov angle.).
+//    Therefore, we will have two candidates of the center positions.
+//
 // 2. Find the ring radius which gives the minimum RMS between the camera
-//  images and the estimated circle.
+//    images and the estimated circle.
+//
 // 3. Select one temporal position which gives smaller RMS as a true temporal
-//  center position.
+//    center position.
+//
 // 4. Changing the center position of a circle on the camera plane from the 
-// determined temporal center position, find the position which gives the 
-// minimum RMS of the fit.
-//
-//
-// --- Remark ---
-//  This method to search for muons is not fully optimized yet. However,
-// it is good idea to use the temporal estimated center position from
-// the Hillas parameters in order to reduce time to search. In addition,
-// This method is faster than the MINUIT.
+//    determined temporal center position, find the position which gives the
+//    minimum RMS of the fit.
 //
 //
 //  Input Containers:
-//   [MGeomCam]
-//   [MHillas]
-//   [MSignalCam]
+//   MGeomCam
+//   MHillas
+//   MSignalCam
 //
 /////////////////////////////////////////////////////////////////////////////
 #include "MMuonSearchPar.h"
 
-#include <fstream>
+#include <TMinuit.h>
 
 #include "MLog.h"
 #include "MLogManip.h"
+
 #include "MHillas.h"
+
 #include "MGeomCam.h"
 #include "MGeomPix.h"
+
 #include "MSignalPix.h"
 #include "MSignalCam.h"
@@ -96 +94 @@
 MMuonSearchPar::MMuonSearchPar(const char *name, const char *title)
 {
-  fName  = name  ? name  : "MMuonSearchPar";
-  fTitle = title ? title : "Muon search parameters";
+    fName  = name  ? name  : "MMuonSearchPar";
+    fTitle = title ? title : "Muon search parameters";
 }
 
@@ -104 +102 @@
 void MMuonSearchPar::Reset()
 {
-  fRadius  = -1.;
-  fDeviation  = -1.;
-  fCenterX = 0.;
-  fCenterY = 0.;
-}
-
-// --------------------------------------------------------------------------
-//
-//  Get the tempolary center of a ring from the Hillas parameters.
-//  Two candidates of the position is returened.
-//
-void MMuonSearchPar::CalcTempCenter(const MHillas &hillas, 
-      Float_t &xtmp1, Float_t &ytmp1, Float_t &xtmp2, Float_t &ytmp2)
-{
-  Float_t a,dx,dy;
-  Float_t tmp_r = 300.;  // assume that the temporal cherenkov angle is 1 deg. (300 mm).
-
-  a = TMath::Tan(hillas.GetDelta());
-
-  dx = a/TMath::Sqrt(tmp_r+a*a)/3.;
-  dy = -tmp_r/TMath::Sqrt(1+a*a)/3.;
-
-  xtmp1 = hillas.GetMeanX() + dx;
-  ytmp1 = hillas.GetMeanY() + dy;
-  xtmp2 = hillas.GetMeanX() - dx;
-  ytmp2 = hillas.GetMeanY() - dy;
+    fRadius    = -1;
+    fDeviation = -1;
+    fCenterX   =  0;
+    fCenterY   =  0;
+}
+
+// --------------------------------------------------------------------------
+//
+// This is a wrapper function to have direct access to the data members
+// in the function calculating the minimization value.
+//
+void MMuonSearchPar::fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
+{
+    const MMuonSearchPar *optim = (MMuonSearchPar*)gMinuit->GetObjectFit();
+    f = optim->Fcn(par);
 }
 
@@ -137 +124 @@
 //  and its RMS for the input position.
 //
-Bool_t MMuonSearchPar::CalcRadius(const MGeomCam &geom, const MSignalCam &evt,
-      Float_t x, Float_t y, Float_t &r, Float_t &sigma)
-{
-  Float_t mean_r=0., dev_r=0., sums=0., tmp=0.;
-
-  const Int_t entries = evt.GetNumPixels();
-
-  for (Int_t i=0; i<entries; i++ ){
-    const MSignalPix &pix = evt[i];
-
-    if (!pix.IsPixelUsed())
-      continue;
-
-    const MGeomPix &gpix = geom[i/*pix.GetPixId()*/];
-
-    tmp=TMath::Sqrt((gpix.GetX()-x)*(gpix.GetX()-x)
-                    +(gpix.GetY()-y)*(gpix.GetY()-y));
-
-    mean_r += pix.GetNumPhotons()*tmp;
-    dev_r  += pix.GetNumPhotons()*tmp*tmp;
-    sums   += pix.GetNumPhotons();
-  }
-
-  if(sums<1.E-10)
-      return kFALSE;
-
-  r = mean_r/sums;
-
-  if(dev_r/sums-(r)*(r)<1.E-10)
-      return kFALSE;
-
-  sigma = TMath::Sqrt(dev_r/sums-(r)*(r)); 
-
-  return kTRUE;
-} 
+Double_t MMuonSearchPar::Fcn(Double_t *par) const
+{
+    const Int_t entries = fSignal.GetSize();
+
+    Double_t meanr=0;
+    Double_t devr =0;
+    Double_t sums =0;
+
+    // It seems that the loop is easy enough for a compiler optimization.
+    // Using pointer arithmetics doesn't improve the speed of the fit.
+    for (Int_t i=0; i<entries; i++ )
+    {
+        Double_t tmp = TMath::Hypot(fX[i]-par[0], fY[i]-par[1]);
+
+        sums  += fSignal[i];
+        meanr += fSignal[i] * tmp;
+        devr  += fSignal[i] * tmp*tmp;
+    }
+
+    par[2] = meanr/sums;
+    par[3] = devr/sums - par[2]*par[2];
+
+    return par[3];
+}
 
 // --------------------------------------------------------------------------
@@ -178 +154 @@
 // RMS of the fit, changing the center position of the circle.
 //
-void MMuonSearchPar::CalcMinimumDeviation
-( const MGeomCam &geom, const MSignalCam &evt, Float_t x, Float_t y,
- Float_t xcog, Float_t ycog, Float_t sigma, Float_t &opt_rad, 
- Float_t &new_sigma, Float_t &newx, Float_t &newy )
-{
-  Float_t delta = 3.;  // 3 mm (1/10 of an inner pixel size) Step to move.
-  Float_t rad_tmp,sig_tmp;
-  Float_t r2;
-
-  while(1) 
-  {
-      r2=(xcog-x)*(xcog-x)+(ycog-y)*(ycog-y);
-      // Exit if the new estimated radius is above 2 deg. (600 mm).
-      if(r2 > 360000.)
-      { 
-          new_sigma=sigma;
-          opt_rad=rad_tmp;
-          newx=x; 
-          newy=y;
-          break;
-      }
-      if(CalcRadius(geom,evt,x,y+delta,rad_tmp,sig_tmp))
-      {
-          if(sig_tmp<sigma)
-          {
-              sigma=sig_tmp; 
-              opt_rad=rad_tmp;
-              y=y+delta;
-              continue;
-          }
-      }
-      if(CalcRadius(geom,evt,x-delta,y,rad_tmp,sig_tmp))
-      {
-          if(sig_tmp<sigma)
-          {
-              sigma=sig_tmp; 
-              opt_rad=rad_tmp;
-              x=x-delta;
-              continue;
-          }
-      }
-      if(CalcRadius(geom,evt,x+delta,y,rad_tmp,sig_tmp))
-      {
-          if(sig_tmp<sigma)
-          {
-              sigma=sig_tmp; 
-              opt_rad=rad_tmp;
-              x=x+delta;
-              continue;
-          }
-      }
-      if(CalcRadius(geom,evt,x,y-delta,rad_tmp,sig_tmp))
-      {
-          if(sig_tmp<sigma)
-          {
-              sigma=sig_tmp; 
-              opt_rad=rad_tmp;
-              y=y-delta;
-              continue;
-          }
-      }
-      new_sigma=sigma;
-      newx=x; 
-      newy=y;
-      break;
+void MMuonSearchPar::CalcMinimumDeviation(const MGeomCam &geom,
+                                          const MSignalCam &evt,
+                                          Double_t &x, Double_t &y,
+                                          Double_t &sigma, Double_t &radius)
+{
+    // ------- Make a temporaray copy of the signal ---------
+    const Int_t n = geom.GetNumPixels();
+
+    fSignal.Set(n);
+    fX.Set(n);
+    fY.Set(n);
+
+    Int_t p=0;
+    for (int i=0; i<n; i++)
+    {
+        const MSignalPix &pix = evt[i];
+        if (pix.IsPixelUsed())
+        {
+            fSignal[p] = pix.GetNumPhotons();
+
+            fX[p] = geom[i].GetX();
+            fY[p] = geom[i].GetY();
+            p++;
+        }
     }
+    fSignal.Set(p);
+
+
+    // ----------------- Setup and call minuit -------------------
+    const Float_t  delta = 30.;  // 3 mm (1/10 of an inner pixel size) Step to move.
+    const Double_t r     = geom.GetMaxRadius()*2;
+
+    // Save gMinuit
+    TMinuit *minsave = gMinuit;
+
+    // Initialize TMinuit with 4 parameters
+    TMinuit minuit;
+    minuit.SetPrintLevel(-1);     // Switch off printing
+    minuit.Command("set nowarn"); // Switch off warning
+    // Define Parameters
+    minuit.DefineParameter(0, "x",     x, delta, -r, r);
+    minuit.DefineParameter(1, "y",     y, delta, -r, r);
+    minuit.DefineParameter(2, "r",     0, 1,      0, 0);
+    minuit.DefineParameter(3, "sigma", 0, 1,      0, 0);
+    // Fix return parameters
+    minuit.FixParameter(2);
+    minuit.FixParameter(3);
+    // Setup Minuit for 'this' and use fit function fcn
+    minuit.SetObjectFit(this);
+    minuit.SetFCN(fcn);
+
+    // Perform Simplex minimization
+    Int_t err;
+    Double_t tmp[2] = { 0, 0 };
+    minuit.mnexcm("simplex", tmp, 2, err);
+
+    // Get resulting parameters
+    Double_t error;
+    minuit.GetParameter(0, x,      error);
+    minuit.GetParameter(1, y,      error);
+    minuit.GetParameter(2, radius, error);
+    minuit.GetParameter(3, sigma,  error);
+
+    sigma = sigma>0 ? TMath::Sqrt(sigma) : 0;
+
+    gMinuit = minsave;
 }
 
@@ -250 +226 @@
 //  Calculation of muon parameters
 //
-void MMuonSearchPar::Calc
-(const MGeomCam &geom, const MSignalCam &evt, const MHillas &hillas)
-{
-  Reset();
-  
-  Float_t xtmp1,xtmp2,ytmp1,ytmp2;
-  Float_t rad,dev,rad2,dev2;
-  Float_t opt_rad,new_sigma,newx,newy;
-    
-  // gets temporaly center
-  CalcTempCenter(hillas,xtmp1,ytmp1,xtmp2,ytmp2);
-  
-  // determine which position will be the true position. Here mainly
-  // the curvature of a muon arc is relied on.
-  CalcRadius(geom, evt, xtmp1,ytmp1,rad,dev);
-  CalcRadius(geom, evt, xtmp2,ytmp2,rad2,dev2);
-  if(dev2<dev){
-    xtmp1=xtmp2; ytmp1=ytmp2; dev=dev2; rad=rad2;
-  }
-  
-  // find the best fit.
-  CalcMinimumDeviation(geom, evt, xtmp1,ytmp1,hillas.GetMeanX(),
-                       hillas.GetMeanY(), dev, opt_rad, new_sigma, 
-                       newx, newy);
-
-  fRadius = opt_rad;
-  fDeviation = new_sigma;
-  
-  fCenterX = newx;
-  fCenterY = newy;
-  
-  //SetReadyToSave();
+void MMuonSearchPar::Calc(const MGeomCam &geom, const MSignalCam &evt,
+                          const MHillas &hillas)
+{
+    Double_t x = hillas.GetMeanX();
+    Double_t y = hillas.GetMeanY();
+
+    // -------------------------------------------------
+    // Keiichi suggested trying to precalculate the Muon
+    // center a bit better, but it neither improves the
+    // fit result nor the speed
+    //
+    // const Float_t tmpr = 300.;  // assume that the temporal cherenkov angle is 1 deg. (300 mm).
+    //
+    // const Double_t a = TMath::Tan(hillas.GetDelta());
+    //
+    // const Double_t dx =     a/TMath::Sqrt(tmpr+a*a)/3.;
+    // const Double_t dy = -tmpr/TMath::Sqrt(1+a*a)/3.;
+    //
+    // Double_t par1[] = { x+dx, y+dy, 0, 0 };
+    // Double_t par2[] = { x-dx, y-dy, 0, 0 };
+    //
+    // const Double_t dev1 = MMuonSearchPar::Fcn(par1);
+    // const Double_t dev2 = MMuonSearchPar::Fcn(par2);
+    //
+    // if (dev1<dev2)
+    // {
+    //     x += dx;
+    //     y += dy;
+    // }
+    // else
+    // {
+    //     x -= dx;
+    //     y -= dy;
+    // }
+    // -------------------------------------------------
+
+    Double_t sigma, rad;
+
+    // find the best fit.
+    CalcMinimumDeviation(geom, evt, x, y, sigma, rad);
+
+    fCenterX   = x;
+    fCenterY   = y;
+    fRadius    = rad;
+    fDeviation = sigma;
+
+    //SetReadyToSave();
 } 
 
@@ -298 +289 @@
     *fLog << all;
     *fLog << "Muon Parameters (" << GetName() << ")" << endl;
-    *fLog << " - Est. Radius   [deg.]  = " << fRadius*geom.GetConvMm2Deg()   << endl;
-    *fLog << " - Deviation     [deg.]  = " << fDeviation*geom.GetConvMm2Deg()   << endl;
-    *fLog << " - Center Pos. X [deg.]  = " << fCenterX*geom.GetConvMm2Deg()  << endl;
-    *fLog << " - Center Pos. Y [deg.]  = " << fCenterY*geom.GetConvMm2Deg()  << endl;
-}
-
-
-
+    *fLog << " - Est. Radius   [deg] = " << fRadius*geom.GetConvMm2Deg()   << endl;
+    *fLog << " - Deviation     [deg] = " << fDeviation*geom.GetConvMm2Deg()   << endl;
+    *fLog << " - Center Pos. X [deg] = " << fCenterX*geom.GetConvMm2Deg()  << endl;
+    *fLog << " - Center Pos. Y [deg] = " << fCenterY*geom.GetConvMm2Deg()  << endl;
+}

trunk/MagicSoft/Mars/mmuon/MMuonSearchPar.h

@@ -4 +4 @@
 #ifndef MARS_MParContainer
 #include "MParContainer.h"
+#endif
+
+#ifndef MARS_MArrayF
+#include "MArrayF.h"
 #endif
 
@@ -18 +22 @@
     Float_t fCenterY;   // An estimated center position in Y of the muon ring [mm]
 
+    MArrayF fSignal;    //! Temporary storage for signal
+    MArrayF fX;         //! Temporary storage for pixels X-position
+    MArrayF fY;         //! Temporary storage for pixels Y-position
+
+    static void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
+    Double_t Fcn(Double_t *par) const;
+
 public:
     MMuonSearchPar(const char *name=NULL, const char *title=NULL);
@@ -28 +39 @@
     Float_t GetCenterY()   const { return fCenterY; }
 
-    void   CalcTempCenter(const MHillas &hillas, Float_t &xtmp1, 
-                          Float_t &ytmp1, Float_t &xtmp2, Float_t &ytmp2);
-    Bool_t CalcRadius(const MGeomCam &geom, const MSignalCam &evt, Float_t x,
-                      Float_t y, Float_t &r, Float_t &sigma);
     void   CalcMinimumDeviation(const MGeomCam &geom, const MSignalCam &evt,
-                             Float_t x, Float_t y, Float_t xcog, 
-                             Float_t ycog, Float_t sigma, Float_t &opt_rad, 
-                             Float_t &new_sigma, Float_t &newx, Float_t &newy);
+                                Double_t &x, Double_t &y, Double_t &sigma, Double_t &rad);
+
     void   Calc(const MGeomCam &geom, const MSignalCam &evt,
                 const MHillas &hillas);

trunk/MagicSoft/Mars/mmuon/MMuonSearchParCalc.cc

@@ -17 +17 @@
 !
 !   Author(s): Keiichi Mase 09/2004 <mailto:mase@mppmu.mpg.de>
-!              Markus Meyer 09/2004 <mailto:meyer@astro.uni-wuerzburg.de>
+!   Author(s): Markus Meyer 09/2004 <mailto:meyer@astro.uni-wuerzburg.de>
 !
-!   Copyright: MAGIC Software Development, 2000-2004
+!   Copyright: MAGIC Software Development, 2000-2005
 !
 !
@@ -41 +41 @@
 //
 //  Input Containers:
-//   [MGeomCam]
-//   [MHillas]
-//   [MSignalCam]
+//   MGeomCam
+//   MHillas
+//   MSignalCam
 //
 //  Output Containers:
-//   [MMuonSearchPar]
+//   MMuonSearchPar
 //
 //////////////////////////////////////////////////////////////////////////////
 #include "MMuonSearchParCalc.h"
 
-#include <fstream>
+#include "MParList.h"
 
-#include "MParList.h"
+#include "MLog.h"
+#include "MLogManip.h"
+
 #include "MGeomCam.h"
 #include "MSignalCam.h"
 #include "MMuonSearchPar.h"
-#include "MLog.h"
-#include "MLogManip.h"
 
 using namespace std;
@@ -71 +71 @@
 // Default constructor. 
 //
-MMuonSearchParCalc::MMuonSearchParCalc
-(const char *mupar, const char *name, const char *title)
+MMuonSearchParCalc::MMuonSearchParCalc(const char *name, const char *title)
   : fHillas(NULL), fMuonPar(NULL)
 {
     fName  = name  ? name  : gsDefName.Data();
     fTitle = title ? title : gsDefTitle.Data();
-
-    fMuparName  = mupar;
-    fHillasInput = "MHillas";
 }
 
@@ -86 +82 @@
 Int_t MMuonSearchParCalc::PreProcess(MParList *pList)
 {
-    fHillas = (MHillas*)pList->FindObject(fHillasInput, "MHillas");
+    fHillas = (MHillas*)pList->FindObject("MHillas");
     if (!fHillas)
     {
-        *fLog << err << fHillasInput << " [MHillas] not found... aborting." << endl;
+        *fLog << err << "MHillas not found... aborting." << endl;
         return kFALSE;
     }
@@ -107 +103 @@
     }
 
-    fMuonPar = (MMuonSearchPar*)pList->FindCreateObj("MMuonSearchPar", fMuparName);
+    fMuonPar = (MMuonSearchPar*)pList->FindCreateObj("MMuonSearchPar");
     if (!fMuonPar)
         return kFALSE;
@@ -118 +114 @@
 Int_t MMuonSearchParCalc::Process()
 {
-
-  fMuonPar->Calc(*fGeomCam, *fSignalCam, *fHillas);
-
-  return kTRUE;
+    fMuonPar->Calc(*fGeomCam, *fSignalCam, *fHillas);
+    return kTRUE;
 }
-

trunk/MagicSoft/Mars/mmuon/MMuonSearchParCalc.h

@@ -14 +14 @@
 {
 private:
-    MGeomCam    *fGeomCam;
-    MSignalCam  *fSignalCam;
-    MHillas     *fHillas;      //! Pointer to the source independent hillas parameters
-    MMuonSearchPar *fMuonPar;  //! Pointer to the output container for the new image parameters
-
-    TString fMuparName;
-    TString fHillasInput;
+    MGeomCam       *fGeomCam;     //!
+    MSignalCam     *fSignalCam;   //!
+    MHillas        *fHillas;      //! Pointer to the source independent hillas parameters
+    MMuonSearchPar *fMuonPar;     //! Pointer to the output container for the new image parameters
 
     Int_t PreProcess(MParList *plist);
@@ -26 +23 @@
 
 public:
-    MMuonSearchParCalc(const char *mupar="MMuonSearchPar",
-                       const char *name=NULL, const char *title=NULL);
-
-    void SetInput(TString hilname) { fHillasInput = hilname; }
+    MMuonSearchParCalc(const char *name=NULL, const char *title=NULL);
 
     ClassDef(MMuonSearchParCalc, 0) // task to calculate muon parameters

trunk/MagicSoft/Mars/mmuon/MMuonSetup.cc

@@ -16 +16 @@
 !
 !
-!   Author(s): Keiichi Mase 10/2004 <mailto:mase@mppmu.mpg.de>
-!              Markus Meyer 10/2004 <mailto:meyer@astro.uni-wuerzburg.de>
+!   Author(s): Markus Meyer 04/2005 <mailto:meyer@astro.uni-wuerzburg.de>
 !
-!   Copyright: MAGIC Software Development, 2000-2004
+!   Copyright: MAGIC Software Development, 2000-2005
 !
 !
@@ -28 +27 @@
 // MMuonSetup
 //
-// Storage Container for parameter for the muon analysis
-//
-//
-//  Input Containers:
-//   [MSignalCam]
+// Storage Container for setup parameter for the muon analysis
 //
 /////////////////////////////////////////////////////////////////////////////
@@ -50 +45 @@
 //
 MMuonSetup::MMuonSetup(const char *name, const char *title)
+    : fMargin(0.2), fThresholdArcPhi(0.1), fThresholdArcWidth(0.1)
 {
     fName  = name  ? name  : "MMuonSetup";
-    fTitle = title ? title : "Muon calibration parameters";
-
-    fMargin = 60.;  // in mm
-    fArcPhiThres  = 30.;
-    fArcWidthThres  = 30.;
-    fArcPhiBinNum = 20;
-    fArcPhiHistStartVal = -180.; // deg.
-    fArcPhiHistEndVal   = 180.;  // deg.
-    fArcWidthBinNum = 28;
-    fArcWidthHistStartVal = 0.3; // deg.
-    fArcWidthHistEndVal   = 1.7; // deg.
+    fTitle = title ? title : "Muon calibration setup parameters";
 }
 
 // --------------------------------------------------------------------------
 //
-MMuonSetup::~MMuonSetup()
+// Read the setup from a TEnv, eg:
+//   MMuonSetup.Margin:            0.2
+//   MMuonSetup.ThresholdArcPhi:   0.1
+//   MMuonSetup.ThresholdArcWidth: 0.1
+//
+Int_t MMuonSetup::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
 {
+    Bool_t rc = kFALSE;
+    if (IsEnvDefined(env, prefix, "Margin", print))
+    {
+        rc = kTRUE;
+        fMargin = GetEnvValue(env, prefix, "Margin", fMargin);
+    }
+    if (IsEnvDefined(env, prefix, "ThresholdArcPhi", print))
+    {
+        rc = kTRUE;
+        fThresholdArcPhi = GetEnvValue(env, prefix, "ThresholdArcPhi", fThresholdArcPhi);
+    }
+    if (IsEnvDefined(env, prefix, "ThresholdArcWidth", print))
+    {
+        rc = kTRUE;
+        fThresholdArcWidth = GetEnvValue(env, prefix, "ThresholdArcWidth", fThresholdArcWidth);
+    }
+
+    return rc;
 }

trunk/MagicSoft/Mars/mmuon/MMuonSetup.h

@@ -9 +9 @@
 {
 private:
-  Float_t fMargin;        // margin to evaluate muons [mm]. The defaut value is 60 mm, corresponding to 0.2 deg. This value can be changed by using the function of SetMargin
-  Float_t fArcPhiThres;   // The threshold value to define arc phi
-  Float_t fArcWidthThres; // The threshold value to define arc width
+    Float_t fMargin;            // [deg] margin to evaluate muons. The defaut value is 60 mm, corresponding to 0.2 deg. This value can be changed by using the function of SetMargin
+    Float_t fThresholdArcPhi;   // [deg] The threshold value to define arc phi
+    Float_t fThresholdArcWidth; // [deg] The threshold value to define arc width
 
+    Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print);
 
 public:
-  MMuonSetup(const char *name=NULL, const char *title=NULL);
-  ~MMuonSetup();
-  
-  Int_t   fArcPhiBinNum;         // The bin number for the histogram of arc phi. You may change this value. However, if you change this, YOU ALSO HAVE TO CHANGE THE THRESHOLD VALUE TO GET ARC LENGTH.
-  Int_t   fArcWidthBinNum;       // The bin number for the histogram of arc wid
-  Float_t fArcPhiHistStartVal;   // The starting value for the histogram of arc phi
-  Float_t fArcPhiHistEndVal;     // The end value for the histogram of arc phi
-  Float_t fArcWidthHistStartVal; // The starting value for the histogram of arc width
-  Float_t fArcWidthHistEndVal;   // The end value for the histogram of arc width
+    MMuonSetup(const char *name=NULL, const char *title=NULL);
 
-  Float_t GetMargin()         const { return fMargin; }
-  Float_t GetArcPhiThres()    const { return fArcPhiThres; }
-  Float_t GetArcWidthThres()  const { return fArcWidthThres; }
-  Float_t GetArcPhiBinNum()   const { return fArcPhiBinNum; }
-  Float_t GetArcWidthBinNum() const { return fArcWidthBinNum; }
-  
-  void    SetMargin(Float_t margin)       { fMargin = margin; }
-  void    SetArcPhiThres(Float_t thres)   { fArcPhiThres = thres; }
-  void    SetArcWidthThres(Float_t thres) { fArcWidthThres = thres; }
-  void    SetArcPhiBinNum(Int_t num)      { fArcPhiBinNum = num; }
-  void    SetArcWidthBinNum(Int_t num)    { fArcWidthBinNum = num; }
-  
-  ClassDef(MMuonSetup, 1) // Container to hold muon parameters for the whole sample
+    // Getter
+    Float_t GetMargin() const { return fMargin; }
+    Float_t GetThresholdArcPhi() const { return fThresholdArcPhi; }
+    Float_t GetThresholdArcWidth() const { return fThresholdArcWidth; }
+
+    // Setter
+    void SetMargin(Float_t margin)           { fMargin = margin; }
+    void SetThresholdArcPhi(Float_t thres)   { fThresholdArcPhi = thres; }
+    void SetThresholdArcWidth(Float_t thres) { fThresholdArcWidth = thres; }
+
+    ClassDef(MMuonSetup, 1) // Container to hold setup for muon analysis
 };
     

trunk/MagicSoft/Mars/mmuon/Makefile

@@ -30 +30 @@
            MMuonCalibParCalc.cc \
            MHMuonPar.cc \
+           MHSingleMuon.cc \
            MMuonSetup.cc
 

trunk/MagicSoft/Mars/mmuon/MuonLinkDef.h

@@ -6 +6 @@
 
 #pragma link C++ class MMuonSearchPar+;
-#pragma link C++ class MMuonSearchParCalc+;
 #pragma link C++ class MMuonCalibPar+;
-#pragma link C++ class MMuonCalibParCalc+;
-#pragma link C++ class MHMuonPar+;
 #pragma link C++ class MMuonSetup+;
 
+#pragma link C++ class MMuonSearchParCalc+;
+#pragma link C++ class MMuonCalibParCalc+;
+
+#pragma link C++ class MHMuonPar+;
+#pragma link C++ class MHSingleMuon+;
+
 #endif
-
-
-
-
-
-
-

trunk/MagicSoft/Mars/star.rc

@@ -46 +46 @@
 
 # -------------------------------------------------------------------------
+# Setup the muon analysis here
+# -------------------------------------------------------------------------
+# MMuonSetup.Margin:            0.2
+# MMuonSetup.ThresholdArcPhi:   0.1
+# MMuonSetup.ThresholdArcWidth: 0.1
+# BinningRadius.Raw:            20  0.5  1.5
+# BinningArcWidth.Raw:          60  0.0  0.3
+# BinningRingBroadening.Raw:    20  0.5  1.5
+# BinningSizeVsArcRadius.Raw:   20  0.5  1.5
+# BinningsArcPhi.Raw:           21 -180  180
+# BinningsMuonWidth.Raw:        28  0.3  1.7
+
+# -------------------------------------------------------------------------
 # -------------------------------------------------------------------------
 MJStar.MImgCleanStd.CleanLevel1: 4.5