Changeset 5332

Timestamp:

11/01/04 22:31:44 (20 years ago)

Author:

mase

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (2 diffs)
• mmuon/MMuonCalibPar.cc (modified) (4 diffs)
• mmuon/MMuonCalibParCalc.cc (modified) (11 diffs)
• mmuon/MMuonCalibParCalc.h (modified) (2 diffs)
• mmuon/Makefile (modified) (2 diffs)
• mmuon/MuonLinkDef.h (modified) (1 diff)

trunk/MagicSoft/Mars/Changelog

@@ -23 +23 @@
  2004/11/01: Markus Meyer and Keiichi Mase
 
-   * mmuon/MMuonCalibParCalc.[h,cc],
+   * mmuon/MMuonCalibParCalc.[h,cc]
      - changed according to Thomas Bretz's suggestion.
        The main part of the calculation is moved from MMuonCalibPar to 
@@ -31 +31 @@
        that we should get the Arc Width with uncleaned images.
 
+   * mmuon/MMuonCalibPar.[h,cc]
+     - changed according to Thomas Bretz's suggestion as written above.
+
    * mmuon/MMuonSearchPar.[h,cc]
      - changed according to Thomas Bretz's suggestion.
        The pointer of the funciotn is changed to the reference.
-
-   * mmuon/MMuonCalibPar.[h,cc],
-     - changed according to Thomas Bretz's suggestion as written above.
-
 
 

trunk/MagicSoft/Mars/mmuon/MMuonCalibPar.cc

@@ -46 +46 @@
 #include <fstream>
 
-#include <TH1.h>
-#include <TF1.h>
-#include <TMinuit.h>
-
 #include "MLog.h"
 #include "MLogManip.h"
-#include "MGeomCam.h"
-#include "MGeomPix.h"
 #include "MCerPhotEvt.h"
 #include "MCerPhotPix.h"
@@ -91 +85 @@
     fHistWidth->UseCurrentStyle();
 
-    fEnableImpactCalc = kFALSE; // By default the calculation of impact parameter is skipped.
-    fDisablePreCuts = kFALSE;   // By default the pre cuts will be applied.
-    fUseCleanForWidth = kFALSE; // By default all the pixels will be used for the histogram of arc width.
-
     fMargin = 60.;  // in mm
     fArcPhiThres  = 100.;
@@ -102 +92 @@
     fArcPhiHistEndVal   = 180.;  // deg.
     fArcWidthBinNum = 28;
-    fArcWidthHistStartVal = 0.3;   // deg.
-    fArcWidthHistEndVal   = 1.7;   // deg.
+    fArcWidthHistStartVal = 0.3; // deg.
+    fArcWidthHistEndVal   = 1.7; // deg.
 }
 
@@ -127 +117 @@
   fUseUnmap    = kFALSE;
   fPeakPhi     =  0.;
+
+  fHistPhi->Reset();
+  fHistWidth->Reset();
 }
-
-// --------------------------------------------------------------------------
-// 
-//  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. You can use the image after the cleaning by using the function of 
-// UseCleanForWidth(). See the manual of MMuonCalibParCalc.
-// 
-void MMuonCalibPar::FillHist
-( const MGeomCam &geom, const MCerPhotEvt &evt, 
-  const MMuonSearchPar &musearch)
-{
-  // preparation for a histgram
-  MBinning binsphi;
-  binsphi.SetEdges(fArcPhiBinNum, fArcPhiHistStartVal, fArcPhiHistEndVal);
-  binsphi.Apply(*fHistPhi);
-
-  MBinning binswid;
-  binswid.SetEdges(fArcWidthBinNum, fArcWidthHistStartVal, fArcWidthHistEndVal);
-  binswid.Apply(*fHistWidth);
-
- const Int_t entries = evt.GetNumPixels();
-
-  // the position of the center of a muon ring
-  const Float_t cenx = musearch.GetCenterX();
-  const Float_t ceny = musearch.GetCenterY();
-  
-  for (Int_t i=0; i<entries; i++ )
-    {
-      MCerPhotPix &pix = (evt)[i];
-      
-      const MGeomPix &gpix = (geom)[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 < musearch.GetRadius() + fMargin && dist > musearch.GetRadius() - fMargin)
-        {
-          // 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())
-            {
-              fUseUnmap = kTRUE;
-              continue;
-            }
-          fHistPhi->Fill(ang*kRad2Deg, pix.GetNumPhotons());
-          fMuonSize += pix.GetNumPhotons();
-        }
-
-      if(pix.GetPixId()>397)
-        continue;  // use only the inner pixles
-      
-      if(fUseCleanForWidth)
-        {
-          if(!pix.IsPixelUsed())
-            continue;
-        }
-
-      fHistWidth->Fill(dist*geom.GetConvMm2Deg(), pix.GetNumPhotons());
-    }
-
-
-  // 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<fArcPhiBinNum+1; i++)
-      {
-        Float_t rougherr  = TMath::Sqrt(TMath::Abs(fHistPhi->GetBinContent(i))*ADC2PhEl)/ADC2PhEl*Ffactor;
-        {
-          fHistPhi->SetBinError(i, rougherr);
-        }
-      }
-    for(Int_t i=0; i<fArcWidthBinNum+1; i++)
-      {
-        Float_t rougherr = TMath::Sqrt(TMath::Abs(fHistWidth->GetBinContent(i))*ADC2PhEl)/ADC2PhEl*Ffactor;
-        {
-          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 MMuonCalibPar::CalcPhi
-(const MGeomCam &geom, const MCerPhotEvt &evt,
- const MMuonSearchPar &musearch)
-{
-  Float_t convbin2val = (fArcPhiHistEndVal-fArcPhiHistStartVal)/
-    (Float_t)fArcPhiBinNum;
-
-    // adjust the peak to 0
-    Float_t maxval = 0.;
-    Int_t   maxbin = 0;
-    maxval = fHistPhi->GetMaximum();
-    maxbin = fHistPhi->GetMaximumBin();
-    fPeakPhi = 180.-(Float_t)(maxbin-1)*convbin2val; 
-    TArrayD tmp;
-    tmp.Set(fArcPhiBinNum+1);
-    for(Int_t i=1; i<fArcPhiBinNum+1; i++)
-      {
-        tmp[i] = fHistPhi->GetBinContent(i);
-      }
-    for(Int_t i=1; i<fArcPhiBinNum+1; i++)
-      {
-        Int_t id;
-        id = i + (maxbin-(Int_t)((Float_t)fArcPhiBinNum/2.)-1);
-        if(id>fArcPhiBinNum)
-          {
-            id-=(fArcPhiBinNum);
-          }
-        if(id<=0)
-          {
-            id+=(fArcPhiBinNum);
-          }
-        fHistPhi->SetBinContent(i,tmp[id]);
-      }
-    maxbin = (Int_t)((Float_t)fArcPhiBinNum/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<fArcPhiBinNum+1; i++)
-    {
-      Float_t content = fHistPhi->GetBinContent(i);
-      Float_t content_pre = fHistPhi->GetBinContent(i-1);
-      
-      if(content > fArcPhiThres && content_pre < fArcPhiThres)
-        {
-          startfitval = (Float_t)(i-1)*convbin2val+fArcPhiHistStartVal;
-        }
-      if(i==maxbin)
-        IsInMaxim = kTRUE;
-      
-      if(content < fArcPhiThres && IsInMaxim == kTRUE)
-        {
-          endfitval = (Float_t)(i-1)*convbin2val+fArcPhiHistStartVal;
-          break;
-        }
-      endfitval = fArcPhiHistEndVal;
-    }
-  
-  effbinnum = (Int_t)((endfitval-startfitval)/convbin2val);
-  
-  fArcPhi = effbinnum*convbin2val;
-  fArcLength = fArcPhi*3.1415926/180.*musearch.GetRadius()*geom.GetConvMm2Deg();
-  
-  if(fEnableImpactCalc)
-      CalcImpact(geom, musearch, 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 MMuonCalibPar::CalcImpact
-( const MGeomCam &geom, const MMuonSearchPar &musearch, 
-  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 = musearch.GetRadius()*geom.GetConvMm2Deg()*3.1415926/180.;
-  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.);
-  
-  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 = musearch.GetRadius()*geom.GetConvMm2Deg()*3.1415926/180.;
-  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.);
-  
-  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)
-    {
-      fEstImpact = -1.;
-    }
-  if(chi2 > chi1)
-    {
-      fEstImpact = estip1;
-      fChiArcPhi = chi1;
-    }
-  else
-    {
-      fEstImpact = estip2;
-      fChiArcPhi = 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 MMuonCalibPar::CalcWidth
-(const MGeomCam &geom, const MCerPhotEvt &evt,
- const MMuonSearchPar &musearch)
-{
-  Float_t convbin2val = (fArcWidthHistEndVal - fArcWidthHistStartVal)
-    /(Float_t)fArcWidthBinNum;
-
-    // determination of fitting region
-    Int_t maxbin = fHistWidth->GetMaximumBin();
-    Float_t startfitval = 0.;
-    Float_t endfitval   = 0.;
-    Bool_t   IsInMaxim = kFALSE;
-    Int_t    effbinnum = 0;
-    for(Int_t i=1; i<fArcWidthBinNum+1; i++)
-      {
-        Float_t content = fHistWidth->GetBinContent(i);
-        Float_t content_pre = fHistWidth->GetBinContent(i-1);
-
-        if(content > fArcWidthThres)
-          effbinnum++;
-
-        if(content > fArcWidthThres && content_pre < fArcWidthThres)
-          {
-            startfitval = (Float_t)(i-4)*convbin2val + fArcWidthHistStartVal;
-            if(startfitval<0.) startfitval = 0.;
-          }
-        if(i==maxbin)
-          IsInMaxim = kTRUE;
-
-        if(content < fArcWidthThres && IsInMaxim == kTRUE)
-          {
-            endfitval = (Float_t)(i+2)*convbin2val + fArcWidthHistStartVal;
-            if(endfitval>180.) endfitval = 180.;
-            break;
-          }
-        endfitval = fArcWidthHistEndVal;
-      }
-    effbinnum = (Int_t)((endfitval-startfitval)/convbin2val);
-
-    TF1 f1("f1","gaus",startfitval,endfitval);
-
-    fHistWidth->Fit("f1","QR","",startfitval,endfitval);
-    
-    if(effbinnum>3)
-      fChiArcWidth = 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 MMuonCalibPar::Calc
-(const MGeomCam &geom, const MCerPhotEvt &evt, 
- MMuonSearchPar &musearch, const Float_t *cuts)
-{
-  // sanity check
-  if(evt.GetNumPixels() < 3)
-    return kCONTINUE;
-
-  // If an event does not seem to be like muon, the calculation will be skipped.
-  if(musearch.IsNoMuon()) 
-    return kCONTINUE; 
-
-  // Pre Cuts 1
-  if(!fDisablePreCuts)
-    {
-      if(musearch.GetRadius() < cuts[0] || musearch.GetRadius() > cuts[1])
-        {
-          musearch.SetNoMuon();
-          return kCONTINUE;
-        }
-      if(musearch.GetDeviation() > cuts[2])
-        {
-          musearch.SetNoMuon();
-          return kCONTINUE;
-        }
-    }
-
-  Reset();
-
-  FillHist(geom,evt,musearch);
-
-  CalcPhi(geom,evt,musearch);
-
-  // Pre Cuts 2
-  if(!fDisablePreCuts)
-    {
-      if(fMuonSize < cuts[3] || fArcPhi < cuts[4])
-        {
-          musearch.SetNoMuon();
-          return kCONTINUE;
-        }
-    }
-
-  fArcWidth = CalcWidth(geom,evt,musearch); 
-  
-  SetReadyToSave();
-
-  return kTRUE;
-} 
 
 void MMuonCalibPar::Print(Option_t *) const

trunk/MagicSoft/Mars/mmuon/MMuonCalibParCalc.cc

@@ -30 +30 @@
 // Task to calculate the muon parameters
 //
-//  This class allows you to get more muon information especially usefule for
+//  This class allows you to get more muon information especially useful for
 // the calibration of our telescope. This class store the information into the
-// container of MMuonCalibPar. The actual calculation is done in the class of 
-// MMuonCalibPar. Please see the manual.
+// container of MMuonCalibPar. 
 //
 //  In order to make this class work, we need the information of the arc
@@ -53 +52 @@
 // the followings;
 //
+//   mucalibcalc.SetMargin(60.);
+//
+//  You can retrieve the histogram (TH1F) using the function of GetHistPhi() 
+// (also GetHistWid()). Therefore, you can draw the histogram such as
+//
 //   MParList  plist;
 //   MMuonCalibPar muparcalib;
 //   plist.AddToList(&muparcalib);.
-//   muparcalib.SetMargin(60.);
-//
-//  You can retrieve the histogram (TH1F) using the function of GetHistPhi() 
-// (also GetHistWid()). Therefore, you can draw the histogram such as
-//
 //   muparcalib.GetHistPhi().Draw();.
 //
@@ -73 +72 @@
 // namely;
 //
-//   muparcalib.EnableImpactCalc();.
+//   mucalibcalc.EnableImpactCalc();.
 //
 //  In addition, for the faster comutation, pre cuts to select the candidates
@@ -86 +85 @@
 // function of DisablePreCuts(), namely;
 //  
-//   muparcalib.DisablePreCuts();.
+//   mucalibcalc.DisablePreCuts();.
 //
 // 
@@ -97 +96 @@
 // done correctly.
 // 
-//  Because of this rough treatment of the error, it may be better to use
-// the image after the cleaning for the evaluation of the Arc Width. There
-// is a possibility to use the image after the cleaning for the Arc Width.
-// You can use the function of UseCleanForWid(). However, this treatment
-// should be temporal and if the problem of the error is fixed, this function
-// will be removed.
-//
 //
 //  Input Containers:
@@ -119 +111 @@
 #include <fstream>
 
+#include <TH1.h>
+#include <TF1.h>
+#include <TMinuit.h>
+
 #include "MParList.h"
 
 #include "MGeomCam.h"
+#include "MGeomPix.h"
 #include "MSrcPosCam.h"
 #include "MCerPhotEvt.h"
@@ -128 +125 @@
 #include "MLog.h"
 #include "MLogManip.h"
+#include "MBinning.h"
 
 using namespace std;
@@ -141 +139 @@
 //
 MMuonCalibParCalc::MMuonCalibParCalc(const char *name, const char *title)
-  : fCerPhotName("MCerPhotEvt")
+  : fNameCerPhot("MCerPhotEvt")
 {
     fName  = name  ? name  : gsDefName.Data();
@@ -151 +149 @@
     fPreCuts[3]     = 2000.;
     fPreCuts[4]     = 150.;
+
+    fMargin = 60.;
+    fArcPhiThres = 100.;
+    fArcWidthThres = 100.;
+
+    fEnableImpactCalc = kFALSE; // By default the calculation of impact parameter is skipped.
+    fDisablePreCuts =   kFALSE; // By default the pre cuts will be applied.
 }
 
@@ -157 +162 @@
 Int_t MMuonCalibParCalc::PreProcess(MParList *pList)
 {
-    fCerPhotEvt = (MCerPhotEvt*)pList->FindObject(AddSerialNumber(fCerPhotName), "MCerPhotEvt");
-    if (!fCerPhotEvt)
-    {
-        *fLog << err << fCerPhotName << " [MCerPhotEvt] not found... aborting." << endl;
-        return kFALSE;
-    }
-
-    fGeomCam = (MGeomCam*)pList->FindObject("MGeomCam");
-    if (!fGeomCam)
-    {
-        *fLog << err << "MGeomCam not found... aborting." << endl;
-        return kFALSE;
-    }
-
-    fMuonCalibPar = (MMuonCalibPar*)pList->FindCreateObj("MMuonCalibPar");
-    if (!fMuonCalibPar)
-        return kFALSE;
-
-    fMuonSearchPar = (MMuonSearchPar*)pList->FindCreateObj("MMuonSearchPar");
-    if (!fMuonSearchPar)
-        return kFALSE;
-
-    return kTRUE;
-}
+  fCerPhotEvt = (MCerPhotEvt*)pList->FindObject(AddSerialNumber(fNameCerPhot), "MCerPhotEvt");
+  if (!fCerPhotEvt)
+    {
+      *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl;
+      return kFALSE;
+    }
+  
+  fGeomCam = (MGeomCam*)pList->FindObject("MGeomCam");
+  if (!fGeomCam)
+    {
+      *fLog << dbginf << "MGeomCam (Camera Geometry) missing in Parameter List... aborting." << endl;
+      return kFALSE;
+    }
+  
+  fMuonCalibPar = (MMuonCalibPar*)pList->FindCreateObj("MMuonCalibPar", "MMuonCalibPar");
+  if (!fMuonCalibPar)
+    {
+      *fLog << dbginf << "MMuonCalibPar missing in Parameter List... aborting." << endl;
+      return kFALSE;
+    }
+  
+  fMuonSearchPar = (MMuonSearchPar*)pList->FindCreateObj("MMuonSearchPar", "MMuonSearchPar");
+  if (!fMuonSearchPar)
+    {
+      *fLog << dbginf << "MMuonSearchPar missing in Parameter List... aborting." << endl;
+      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 = fMuonCalibPar->fArcPhiBinNum;
+  Int_t binnumwid = fMuonCalibPar->fArcWidthBinNum;
+
+  // preparation for a histgram
+  MBinning binsphi;
+  binsphi.SetEdges(binnumphi, 
+                   fMuonCalibPar->fArcPhiHistStartVal, 
+                   fMuonCalibPar->fArcPhiHistEndVal);
+  binsphi.Apply(*(fMuonCalibPar->fHistPhi));
+
+  MBinning binswid;
+  binswid.SetEdges(binnumwid, 
+                   fMuonCalibPar->fArcWidthHistStartVal, 
+                   fMuonCalibPar->fArcWidthHistEndVal);
+  binswid.Apply(*(fMuonCalibPar->fHistWidth));
+
+  const Int_t entries = (*fCerPhotEvt).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++ )
+    {
+      MCerPhotPix &pix = (*fCerPhotEvt)[i];
+      
+      const MGeomPix &gpix = (*fGeomCam)[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() + fMuonCalibPar->GetMargin() 
+         && dist > (*fMuonSearchPar).GetRadius() - fMuonCalibPar->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();
+        }
+
+      if(pix.GetPixId()>397)
+        continue;  // use only the inner pixles
+      
+      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 = fMuonCalibPar->GetArcPhiThres();
+  Float_t startval = fMuonCalibPar->fArcPhiHistStartVal;
+  Float_t endval = fMuonCalibPar->fArcPhiHistEndVal;
+  Int_t   binnum = fMuonCalibPar->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 = fMuonCalibPar->fArcWidthHistStartVal;
+  Float_t endval = fMuonCalibPar->fArcWidthHistEndVal;
+  Int_t   binnum = fMuonCalibPar->fArcWidthBinNum;
+  Float_t thres = fMuonCalibPar->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);
+
+    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(const Float_t *cuts)
+{
+  // sanity check
+  if((*fCerPhotEvt).GetNumPixels() < 3)
+    return kCONTINUE;
+
+  // If an event does not seem to be like muon, the calculation will be skipped.
+  if((*fMuonSearchPar).IsNoMuon()) 
+    return kCONTINUE; 
+
+  // Pre Cuts 1
+  if(!fDisablePreCuts)
+    {
+      if((*fMuonSearchPar).GetRadius() < cuts[0] || (*fMuonSearchPar).GetRadius() > cuts[1])
+        {
+          (*fMuonSearchPar).SetNoMuon();
+          return kCONTINUE;
+        }
+      if((*fMuonSearchPar).GetDeviation() > cuts[2])
+        {
+          (*fMuonSearchPar).SetNoMuon();
+          return kCONTINUE;
+        }
+    }
+
+  // initialization
+  (*fMuonCalibPar).Reset();
+
+  // Fill signals to histograms
+  FillHist();
+
+  // Calculation of Arc Phi etc...
+  CalcPhi();
+
+  // Pre Cuts 2
+  if(!fDisablePreCuts)
+    {
+      if(fMuonCalibPar->GetMuonSize() < cuts[3] 
+         || fMuonCalibPar->GetArcPhi() < cuts[4])
+        {
+          (*fMuonSearchPar).SetNoMuon();
+          return kCONTINUE;
+        }
+    }
+
+  // Calculation of Arc Width etc...
+  fMuonCalibPar->SetArcWidth(CalcWidth()); 
+  
+  return kTRUE;
+} 
+
 
 // -------------------------------------------------------------------------
@@ -186 +569 @@
 Int_t MMuonCalibParCalc::Process()
 {
-
-  if(!fMuonCalibPar->Calc(*fGeomCam, *fCerPhotEvt, *fMuonSearchPar, fPreCuts))
+  fMuonCalibPar->SetMargin(fMargin);
+  fMuonCalibPar->SetArcPhiThres(fArcPhiThres);
+  fMuonCalibPar->SetArcWidthThres(fArcWidthThres);
+
+  if(!Calc(fPreCuts))
     return kCONTINUE;
 
   return kTRUE;
-}
-
-void MMuonCalibParCalc::SetMargin(Float_t margin)
-{
-  fMuonCalibPar->SetMargin(margin);
-}
-
-void MMuonCalibParCalc::EnableImpactCalc()
-{
-  fMuonCalibPar->EnableImpactCalc();
-}
-
-void MMuonCalibParCalc::DisablePreCuts()
-{
-  fMuonCalibPar->DisablePreCuts();
 }
 

trunk/MagicSoft/Mars/mmuon/MMuonCalibParCalc.h

@@ -20 +20 @@
     MMuonSearchPar *fMuonSearchPar; 
 
+    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
+    Bool_t fEnableImpactCalc; // If true, the impact calculation will be done, which consumes a lot of time.
+    Bool_t fDisablePreCuts;   // If true, the pre cuts to select muons for the calibration will be disabled.
+
     Float_t fPreCuts[5];  // The values for pre cuts.
 
@@ -25 +31 @@
     Int_t Process();
 
-    TString fCerPhotName;
+    TString fNameCerPhot;
 
 public:
     MMuonCalibParCalc(const char *name=NULL, const char *title=NULL);
 
-    void SetMargin(Float_t margin);
-    void EnableImpactCalc();
-    void DisablePreCuts();
+    void SetMargin(Float_t margin)       { fMargin = margin; }
+    void SetArcPhiThres(Float_t thres)   { fArcPhiThres = thres; }
+    void SetArcWidthThres(Float_t thres) { fArcWidthThres = thres; }
+    void EnableImpactCalc()              { fEnableImpactCalc = kTRUE; }
+    void DisablePreCuts()                { fDisablePreCuts = kTRUE; }
     void SetPreCuts(Float_t radcutlow, Float_t radcuthigh, Float_t devcuthigh,
                     Float_t musizecutlow, Float_t arcphicutlow);
 
-    void SetNameCerPhotEvt(const char *name) { fCerPhotName = name; }
+    void SetNameCerPhotEvt(const char *name) { fNameCerPhot = name; }
+
+    void FillHist();
+    void CalcPhi();
+    void CalcImpact(Int_t effbinnum, Float_t startfitval, Float_t endfitval);
+    Float_t CalcWidth();
+    Int_t   Calc(const Float_t *cuts);
 
     ClassDef(MMuonCalibParCalc, 0) // task to calculate muon parameters

trunk/MagicSoft/Mars/mmuon/Makefile

@@ -26 +26 @@
 
 SRCFILES = MMuonSearchPar.cc \
-           MMuonCalibPar.cc 
+           MMuonSearchParCalc.cc \
+           MMuonCalibPar.cc \
+           MMuonCalibParCalc.cc 
 
 ############################################################
@@ -37 +39 @@
 
 mrproper:       clean rmbak
+
+
+
+
+
+
+

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+;
 
 #endif