Ignore:
Timestamp:
04/27/05 16:46:24 (20 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/mmuon
Files:
2 added
14 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Mars/mmuon/MHMuonPar.cc

    r6973 r6979  
    1616!
    1717!
    18 !   Author(s): Wolfgang Wittek, 03/2003 <mailto:wittek@mppmu.mpg.de>
    19 !   Author(s): Thomas Bretz, 04/2003 <mailto:tbretz@astro.uni-wuerzburg.de>
    2018!   Author(s): Markus Meyer, 02/2005 <mailto:meyer@astro.uni-wuerzburg.de>
    21 !
    22 !   Copyright: MAGIC Software Development, 2000-2004
     19!   Author(s): Thomas Bretz, 04/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
     20!
     21!   Copyright: MAGIC Software Development, 2000-2005
    2322!
    2423!
     
    2827//
    2928// MHMuonPar
     29//
    3030// This class is a histogram class for displaying muonparameters.
    3131// The folowing histgrams will be plotted:
     
    3737//
    3838// Inputcontainer:
    39 // MMuonSearchPar
    40 // MMuonCalibPar
     39//   - MGeomCam
     40//   - MMuonSearchPar
     41//   - MMuonCalibPar
    4142//
    4243////////////////////////////////////////////////////////////////////////////
    4344#include "MHMuonPar.h"
    44 
    45 #include <math.h>
    4645
    4746#include <TH1.h>
     
    4948#include <TCanvas.h>
    5049#include <TProfile.h>
     50
    5151#include "MLog.h"
    5252#include "MLogManip.h"
     
    5656#include "MParList.h"
    5757
    58 #include "MHillas.h"
    59 #include "MHMuonPar.h"
    6058#include "MMuonSearchPar.h"
    6159#include "MMuonCalibPar.h"
     
    6967// Setup histograms
    7068//
    71 MHMuonPar::MHMuonPar(const char *name, const char *title):fGeomCam(NULL), fMuonSearchPar(NULL),
    72     fMuonCalibPar(NULL)
     69MHMuonPar::MHMuonPar(const char *name, const char *title) :
     70    fMuonSearchPar(NULL), fMuonCalibPar(NULL)
    7371{
    7472    fName  = name  ? name  : "MHMuonPar";
     
    7674
    7775    fHistRadius.SetName("Radius");
    78     fHistRadius.SetTitle("Radius");
    79     fHistRadius.SetXTitle("Radius[deg]");
     76    fHistRadius.SetTitle("Distribution of Radius'");
     77    fHistRadius.SetXTitle("R [\\circ]");
    8078    fHistRadius.SetYTitle("Counts");
    8179    fHistRadius.SetDirectory(NULL);
     
    8482
    8583    fHistArcWidth.SetName("ArcWidth");
    86     fHistArcWidth.SetTitle("ArcWidth");
    87     fHistArcWidth.SetXTitle("ArcWidth[deg]");
     84    fHistArcWidth.SetTitle("Distribution of ArcWidth");
     85    fHistArcWidth.SetXTitle("W [\\circ]");
    8886    fHistArcWidth.SetYTitle("Counts");
    8987    fHistArcWidth.SetDirectory(NULL);
     
    9189    fHistArcWidth.SetFillStyle(4000);
    9290
    93     fHistBroad.SetName("Ringbroadening");
    94     fHistBroad.SetTitle("Ringbroadening");
    95     fHistBroad.SetXTitle("Radius[deg]");
    96     fHistBroad.SetYTitle("ArcWidth/Radius");
     91    fHistBroad.SetName("RingBroadening");
     92    fHistBroad.SetTitle("Profile of ArcWidth/Radius versus Radius");
     93    fHistBroad.SetXTitle("R [\\circ]");
     94    fHistBroad.SetYTitle("W/R [1]");
    9795    fHistBroad.SetDirectory(NULL);
    9896    fHistBroad.UseCurrentStyle();
     
    10098
    10199    fHistSize.SetName("SizeVsRadius");
    102     fHistSize.SetXTitle("Radius");
    103     fHistSize.SetYTitle("MuonSize");
     100    fHistSize.SetTitle("Profile of Muon Size vs. Radius");
     101    fHistSize.SetXTitle("Rs [[\circ]");
     102    fHistSize.SetYTitle("S [phe]");
    104103    fHistSize.SetDirectory(NULL);
    105104    fHistSize.UseCurrentStyle();
     
    119118    bins.SetEdges(20, 0.5, 1.5);
    120119    bins.Apply(fHistSize);
    121 
    122120}
    123121
     
    129127Bool_t MHMuonPar::SetupFill(const MParList *plist)
    130128{
    131     fGeomCam = (MGeomCam*)plist->FindObject("MGeomCam");
    132     if (!fGeomCam)
     129    MGeomCam *geom = (MGeomCam*)plist->FindObject("MGeomCam");
     130    if (!geom)
    133131    {
    134132        *fLog << warn << "MGeomCam not found... abort." << endl;
    135133        return kFALSE;
    136134    }
     135    fMm2Deg = geom->GetConvMm2Deg();
     136
    137137    fMuonSearchPar = (MMuonSearchPar*)plist->FindObject("MMuonSearchPar");
    138138    if (!fMuonSearchPar)
     
    148148    }
    149149
    150     fMm2Deg = fGeomCam->GetConvMm2Deg();
    151 
    152     ApplyBinning(*plist, "Radius", &fHistRadius);
    153 
    154     ApplyBinning(*plist, "ArcWidth",  &fHistArcWidth);
    155 
    156     ApplyBinning(*plist, "Ringbroadening",  &fHistBroad);
    157 
    158     ApplyBinning(*plist, "SizeVsRadius",  &fHistSize);
     150    ApplyBinning(*plist, "Radius",          &fHistRadius);
     151    ApplyBinning(*plist, "ArcWidth",        &fHistArcWidth);
     152    ApplyBinning(*plist, "RingBroadening",  &fHistBroad);
     153    ApplyBinning(*plist, "SizeVsRadius",    &fHistSize);
    159154
    160155    return kTRUE;
     
    211206    fHistBroad.Draw();
    212207}
    213 
     208/*
    214209TH1 *MHMuonPar::GetHistByName(const TString name)
    215210{
     
    220215    return NULL;
    221216}
    222 
     217*/
  • trunk/MagicSoft/Mars/mmuon/MHMuonPar.h

    r6973 r6979  
    1212#endif
    1313
    14 class MHillas;
    1514class MMuonSearchPar;
    1615class MMuonCalibPar;
     
    2019{
    2120private:
    22     TH1F fHistRadius;  //
     21    TH1F     fHistRadius;    //
     22    TH1F     fHistArcWidth;  //
    2323
    24     TH1F fHistArcWidth;   //
    25 
    26     TProfile fHistBroad;  // Area of used pixels
    27 
     24    TProfile fHistBroad;     // Area of used pixels
    2825    TProfile fHistSize;      // [ratio] concentration ratio: sum of the two highest pixels / fSize
    2926
    30     MGeomCam *fGeomCam; //! Camera geometry for plots (for the moment this is a feature for a loop only!)
    31 
    32     MMuonSearchPar *fMuonSearchPar;//!
    33 
    34     MMuonCalibPar *fMuonCalibPar;//!
     27    MMuonSearchPar *fMuonSearchPar; //!
     28    MMuonCalibPar  *fMuonCalibPar;  //!
    3529
    3630    Float_t fMm2Deg;
     
    4236    Bool_t Fill(const MParContainer *par, const Stat_t w=1);
    4337
    44     TH1 *GetHistByName(const TString name);
     38    //TH1 *GetHistByName(const TString name);
    4539
    46     TH1F &GetHistRadius()  { return fHistRadius; }
    47 
    48     TH1F &GetHistArcWidth()   { return fHistArcWidth; }
    49 
    50     TProfile &GetHistBroad()  { return fHistBroad; }
    51 
    52     TProfile &GetHistSize()      { return fHistSize; }
     40    const TH1F&     GetHistRadius() const    { return fHistRadius; }
     41    const TH1F&     GetHistArcWidth() const  { return fHistArcWidth; }
     42    const TProfile& GetHistBroad() const     { return fHistBroad; }
     43    const TProfile& GetHistSize() const      { return fHistSize; }
    5344
    5445    void Draw(Option_t *opt="");
  • trunk/MagicSoft/Mars/mmuon/MMuonCalibPar.cc

    r6973 r6979  
    1717!
    1818!   Author(s): Keiichi Mase 10/2004 <mailto:mase@mppmu.mpg.de>
    19 !             Markus Meyer 10/2004 <mailto:meyer@astro.uni-wuerzburg.de>
     19!   Author(s): Markus Meyer 10/2004 <mailto:meyer@astro.uni-wuerzburg.de>
    2020!
    21 !   Copyright: MAGIC Software Development, 2000-2004
     21!   Copyright: MAGIC Software Development, 2000-2005
    2222!
    2323!
     
    3030// Storage Container for muon
    3131//
    32 //  This class holds some information for a calibraion using muons. Muons
     32//  This class holds some information for a calibration using muons. Muons
    3333// are identified by using the class of the MMuonSearchParCalc. You can fill
    3434// these information by using the MMuonCalibParCalc. See also these class
    3535// manuals.
    3636//
    37 // 
    38 //  Input Containers:
    39 //   [MGeomCam]
    40 //   [MSignalCam]
    41 //   [MMuonSearchPar]
    42 //
    4337/////////////////////////////////////////////////////////////////////////////
    4438#include "MMuonCalibPar.h"
    4539
    46 #include <fstream>
    47 
    4840#include "MLog.h"
    4941#include "MLogManip.h"
    50 #include "MSignalCam.h"
    51 #include "MSignalPix.h"
    52 #include "MMuonSearchPar.h"
    53 #include "MBinning.h"
    5442
    5543using namespace std;
     
    6654    fTitle = title ? title : "Muon calibration parameters";
    6755
    68     fHistPhi   = new TH1F;
    69     fHistWidth = new TH1F;
    70 
    71     fHistPhi->SetName("HistPhi");
    72     fHistPhi->SetTitle("HistPhi");
    73     fHistPhi->SetXTitle("phi [deg.]");
    74     fHistPhi->SetYTitle("sum of ADC");
    75     fHistPhi->SetDirectory(NULL);
    76     fHistPhi->SetFillStyle(4000);
    77     fHistPhi->UseCurrentStyle();
    78 
    79     fHistWidth->SetName("HistWidth");
    80     fHistWidth->SetTitle("HistWidth");
    81     fHistWidth->SetXTitle("distance from the ring center [deg.]");
    82     fHistWidth->SetYTitle("sum of ADC");
    83     fHistWidth->SetDirectory(NULL);
    84     fHistWidth->SetFillStyle(4000);
    85     fHistWidth->UseCurrentStyle();
    86 
    87 }
    88 
    89 // --------------------------------------------------------------------------
    90 //
    91 MMuonCalibPar::~MMuonCalibPar()
    92 {
    93   delete fHistPhi;
    94   delete fHistWidth;
     56    Reset();
    9557}
    9658
     
    9961void MMuonCalibPar::Reset()
    10062{
    101   fArcLength   = -1.;
    102   fArcPhi      =  0.;
    103   fArcWidth    = -1.;
    104   fChiArcPhi   = -1.;
    105   fChiArcWidth = -1.;
    106   fMuonSize    =  0.;
    107   fEstImpact   = -1.;
    108   fUseUnmap    = kFALSE;
    109   fPeakPhi     =  0.;
    110 
    111   fHistPhi->Reset();
    112   fHistWidth->Reset();
     63    fArcLength   = -1.;
     64    fArcPhi      =  0.;
     65    fArcWidth    = -1.;
     66    fChiArcPhi   = -1.;
     67    fChiArcWidth = -1.;
     68    fMuonSize    =  0.;
     69    fEstImpact   = -1.;
     70    fPeakPhi     =  0.;
    11371}
    11472
     
    11674{
    11775    *fLog << all;
    118     *fLog << "Muon Parameters (" << GetName() << ")"    << endl;
    119     *fLog << " - Arc Length    [deg.]  = " << fArcLength     << endl;
    120     *fLog << " - Arc Phi       [deg.]  = " << fArcPhi     << endl;
    121     *fLog << " - Arc Width     [deg.]  = " << fArcWidth     << endl;
    122     *fLog << " - Chi Arc Phi   [x2/ndf]= " << fChiArcPhi  << endl;
    123     *fLog << " - Chi Arc Width [x2/ndf]= " << fChiArcWidth  << endl;
    124     *fLog << " - Est. I. P.    [m]     = " << fEstImpact  << endl;
    125     *fLog << " - Size of muon          = " << fMuonSize   << endl;
    126     *fLog << " - Peak Phi      [deg.]  = " << fPeakPhi    << endl;
    127     *fLog << " - UseUnmap              = " << fUseUnmap   << endl;
     76    *fLog << "Muon Parameters (" << GetName() << ")"       << endl;
     77    *fLog << " - Arc Length    [deg]   = " << fArcLength   << endl;
     78    *fLog << " - Arc Phi       [deg]   = " << fArcPhi      << endl;
     79    *fLog << " - Arc Width     [deg]   = " << fArcWidth    << endl;
     80    *fLog << " - Chi Arc Phi   [x2/ndf]= " << fChiArcPhi   << endl;
     81    *fLog << " - Chi Arc Width [x2/ndf]= " << fChiArcWidth << endl;
     82    *fLog << " - Est. I. P.    [m]     = " << fEstImpact   << endl;
     83    *fLog << " - Size of muon  [phe]   = " << fMuonSize    << endl;
     84    *fLog << " - Peak Phi      [deg]   = " << fPeakPhi     << endl;
    12885}
    12986
  • trunk/MagicSoft/Mars/mmuon/MMuonCalibPar.h

    r6973 r6979  
    1616{
    1717private:
    18   Float_t fArcLength;     // An arc length of a muon along the arc [deg.]
    19   Float_t fArcPhi;        // A opening angle of a muon arc [deg.]
    20   Float_t fArcWidth;      // A width of a muon [deg.] (1 sigma of gaussian fit)
    21   Float_t fChiArcPhi;     // A chisquare value of the cosine fit for arc phi
    22   Float_t fChiArcWidth;   // A chisquare value of the cosine fit for arc wid
    23   Float_t fMuonSize;      // A SIZE of muon which is defined as a SIZE around the estimated circle
    24   Float_t fEstImpact;     // An estimated impact parameter from the photon distribution along the arc image
    25   Bool_t  fUseUnmap;      // This is a flag to know the Unmapped pixels are used. Refer to the class of MImgCleanStd
    26   Float_t fPeakPhi;       // The angle which indicates the peak position in the estimated circle
     18    Float_t fArcLength;     // An arc length of a muon along the arc [deg.]
     19    Float_t fArcPhi;        // A opening angle of a muon arc [deg.]
     20    Float_t fArcWidth;      // A width of a muon [deg.] (1 sigma of gaussian fit)
     21    Float_t fChiArcPhi;     // A chisquare value of the cosine fit for arc phi
     22    Float_t fChiArcWidth;   // A chisquare value of the cosine fit for arc wid
     23    Float_t fMuonSize;      // A SIZE of muon which is defined as a SIZE around the estimated circle
     24    Float_t fEstImpact;     // An estimated impact parameter from the photon distribution along the arc image
     25    //Bool_t  fUseUnmap;      // This is a flag to know the Unmapped pixels are used. Refer to the class of MImgCleanStd
     26    Float_t fPeakPhi;       // The angle which indicates the peak position in the estimated circle
    2727
    2828public:
    29   MMuonCalibPar(const char *name=NULL, const char *title=NULL);
    30   ~MMuonCalibPar();
    31  
    32   TH1F *fHistPhi;   // Histogram of photon distribution along the arc.
    33   TH1F *fHistWidth; // Histogram of radial photon distribution of the arc.
    34  
    35   void Reset();
    36  
    37   Float_t GetArcLength()      const { return fArcLength; }
    38   Float_t GetArcPhi()         const { return fArcPhi; }
    39   Float_t GetArcWidth()       const { return fArcWidth; }
    40   Float_t GetChiArcPhi()      const { return fChiArcPhi; }
    41   Float_t GetChiArcWidth()    const { return fChiArcWidth; }
    42   Float_t GetMuonSize()       const { return fMuonSize; }
    43   Float_t GetEstImpact()      const { return fEstImpact; }
    44   Bool_t  IsUseUnmap()        const { return fUseUnmap; }
    45   Float_t GetPeakPhi()        const { return fPeakPhi; }
    46   TH1F    *GetHistPhi()       { return fHistPhi; }
    47   TH1F    *GetHistWidth()     { return fHistWidth; }
    48  
    49   void    SetArcLength(Float_t len)       { fArcLength = len; }
    50   void    SetArcPhi(Float_t phi)          { fArcPhi = phi; }
    51   void    SetArcWidth(Float_t wid)        { fArcWidth = wid; }
    52   void    SetChiArcPhi(Float_t chi)       { fChiArcPhi = chi; }
    53   void    SetChiArcWidth(Float_t chi)     { fChiArcWidth = chi; }
    54   void    SetMuonSize(Float_t size)       { fMuonSize = size; }
    55   void    SetEstImpact(Float_t impact)    { fEstImpact = impact; }
    56   void    SetUseUnmap()                   { fUseUnmap = kTRUE; }
    57   void    SetPeakPhi(Float_t phi)         { fPeakPhi = phi; }
    58  
    59   void    Print(Option_t *opt=NULL) const;
     29    MMuonCalibPar(const char *name=NULL, const char *title=NULL);
     30    //~MMuonCalibPar();
    6031
    61   ClassDef(MMuonCalibPar, 1) // Container to hold muon calibration parameters
     32    void Reset();
     33
     34    Float_t GetArcLength()      const { return fArcLength; }
     35    Float_t GetArcPhi()         const { return fArcPhi; }
     36    Float_t GetArcWidth()       const { return fArcWidth; }
     37    Float_t GetChiArcPhi()      const { return fChiArcPhi; }
     38    Float_t GetChiArcWidth()    const { return fChiArcWidth; }
     39    Float_t GetMuonSize()       const { return fMuonSize; }
     40    Float_t GetEstImpact()      const { return fEstImpact; }
     41    //Bool_t  IsUseUnmap()        const { return fUseUnmap; }
     42    Float_t GetPeakPhi()        const { return fPeakPhi; }
     43    //  TH1F    *GetHistPhi()       { return fHistPhi; }
     44    //  TH1F    *GetHistWidth()     { return fHistWidth; }
     45
     46    void    SetArcLength(Float_t len)       { fArcLength = len; }
     47    void    SetArcPhi(Float_t phi)          { fArcPhi = phi; }
     48    void    SetArcWidth(Float_t wid)        { fArcWidth = wid; }
     49    void    SetChiArcPhi(Float_t chi)       { fChiArcPhi = chi; }
     50    void    SetChiArcWidth(Float_t chi)     { fChiArcWidth = chi; }
     51    void    SetMuonSize(Float_t size)       { fMuonSize = size; }
     52    void    SetEstImpact(Float_t impact)    { fEstImpact = impact; }
     53    //void    SetUseUnmap()                   { fUseUnmap = kTRUE; }
     54    void    SetPeakPhi(Float_t phi)         { fPeakPhi = phi; }
     55
     56    void    Print(Option_t *opt=NULL) const;
     57
     58    ClassDef(MMuonCalibPar, 1) // Container to hold muon calibration parameters
    6259};
    6360   
  • trunk/MagicSoft/Mars/mmuon/MMuonCalibParCalc.cc

    r6973 r6979  
    6767//
    6868//
    69 // For the faster computation, by default, the calculation of impact
    70 // parameter is suppressed. If you want to calculate the impact parameter
    71 // from the muon image, you can use the function of EnableImpactCalc(),
    72 // namely;
    73 //
    74 //   mucalibcalc.EnableImpactCalc();.
     69// // For the faster computation, by default, the calculation of impact
     70// // parameter is suppressed. If you want to calculate the impact parameter
     71// // from the muon image, you can use the function of EnableImpactCalc(),
     72// // namely;
     73// //
     74// //   mucalibcalc.EnableImpactCalc();.
    7575//
    7676//  In addition, for the faster comutation, pre cuts to select the candidates
     
    9898//
    9999//  Input Containers:
    100 //   [MGeomCam]
    101 //   [MSignalCam]
    102 //   [MMuonSearchPar]
     100//   MGeomCam
     101//   MSignalCam
     102//   MMuonSearchPar
    103103//
    104104//  Output Containers:
    105 //   [MMuonCalibPar]
     105//   MMuonCalibPar
    106106//
    107107//////////////////////////////////////////////////////////////////////////////
    108108
    109109#include "MMuonCalibParCalc.h"
    110 
    111 #include <fstream>
    112110
    113111#include <TH1.h>
     
    115113#include <TMinuit.h>
    116114
     115#include "MLog.h"
     116#include "MLogManip.h"
     117
    117118#include "MParList.h"
    118119
    119120#include "MGeomCam.h"
    120121#include "MGeomPix.h"
    121 #include "MSrcPosCam.h"
     122
    122123#include "MSignalCam.h"
    123 #include "MMuonSearchPar.h"
     124
    124125#include "MMuonCalibPar.h"
    125126#include "MMuonSetup.h"
    126 #include "MLog.h"
    127 #include "MLogManip.h"
    128 #include "MBinning.h"
     127#include "MMuonSearchPar.h"
     128#include "MHSingleMuon.h"
    129129
    130130using namespace std;
     
    140140//
    141141MMuonCalibParCalc::MMuonCalibParCalc(const char *name, const char *title)
     142//    : fEnableImpactCalc(kFALSE)
    142143{
    143144    fName  = name  ? name  : gsDefName.Data();
    144145    fTitle = title ? title : gsDefTitle.Data();
    145 
    146     fEnableImpactCalc = kFALSE; // By default the calculation of impact parameter is skipped.
    147146}
    148147
     
    151150Int_t MMuonCalibParCalc::PreProcess(MParList *pList)
    152151{
    153   fSignalCam = (MSignalCam*)pList->FindObject("MSignalCam");
    154   if (!fSignalCam)
     152    fGeomCam = (MGeomCam*)pList->FindObject("MGeomCam");
     153    if (!fGeomCam)
    155154    {
    156       *fLog << dbginf << "MSignalCam not found... aborting." << endl;
    157       return kFALSE;
     155        *fLog << err << "MGeomCam not found... abort." << endl;
     156        return kFALSE;
    158157    }
    159  
    160   fGeomCam = (MGeomCam*)pList->FindObject("MGeomCam");
    161   if (!fGeomCam)
     158
     159    fHist = (MHSingleMuon*)pList->FindObject("MHSingleMuon");
     160    if (!fHist)
    162161    {
    163       *fLog << dbginf << "MGeomCam (Camera Geometry) missing in Parameter List... aborting." << endl;
    164       return kFALSE;
     162        *fLog << err << "MHSingleMuon not found... abort." << endl;
     163        return kFALSE;
    165164    }
    166  
    167   fMuonCalibPar = (MMuonCalibPar*)pList->FindCreateObj("MMuonCalibPar", "MMuonCalibPar");
    168   if (!fMuonCalibPar)
    169       return kFALSE;
    170  
    171   fMuonSearchPar = (MMuonSearchPar*)pList->FindCreateObj("MMuonSearchPar", "MMuonSearchPar");
    172   if (!fMuonSearchPar)
    173       return kFALSE;
    174 
    175   fMuonSetup = (MMuonSetup*)pList->FindCreateObj("MMuonSetup", "MMuonSetup");
    176   if (!fMuonSetup)
    177       return kFALSE;
    178  
    179   return kTRUE;
     165
     166    fMuonSetup = (MMuonSetup*)pList->FindObject("MMuonSetup");
     167    if (!fMuonSetup)
     168    {
     169        *fLog << err << "MMuonSetup not found... abort." << endl;
     170        return kFALSE;
     171    }
     172
     173    fMuonCalibPar = (MMuonCalibPar*)pList->FindCreateObj("MMuonCalibPar");
     174    if (!fMuonCalibPar)
     175        return kFALSE;
     176
     177    fMuonSearchPar = (MMuonSearchPar*)pList->FindCreateObj("MMuonSearchPar");
     178    if (!fMuonSearchPar)
     179        return kFALSE;
     180
     181    return kTRUE;
    180182}
    181 
    182 // --------------------------------------------------------------------------
    183 //
    184 //  This function fill the histograms in order to get muon parameters.
    185 // For the evaluation of the Arc Width, we use only the signals in the inner
    186 // part.
    187 //
    188 void MMuonCalibParCalc::FillHist()
    189 {
    190   Float_t MuonSize = 0.;
    191 
    192   Int_t binnumphi = fMuonSetup->fArcPhiBinNum;
    193   Int_t binnumwid = fMuonSetup->fArcWidthBinNum;
    194 
    195   // preparation for a histgram
    196   MBinning binsphi;
    197   binsphi.SetEdges(binnumphi,
    198                    fMuonSetup->fArcPhiHistStartVal,
    199                    fMuonSetup->fArcPhiHistEndVal);
    200   binsphi.Apply(*(fMuonCalibPar->fHistPhi));
    201 
    202   MBinning binswid;
    203   binswid.SetEdges(binnumwid,
    204                    fMuonSetup->fArcWidthHistStartVal,
    205                    fMuonSetup->fArcWidthHistEndVal);
    206   binswid.Apply(*(fMuonCalibPar->fHistWidth));
    207 
    208   const Int_t entries = (*fSignalCam).GetNumPixels();
    209 
    210   // the position of the center of a muon ring
    211   const Float_t cenx = (*fMuonSearchPar).GetCenterX();
    212   const Float_t ceny = (*fMuonSearchPar).GetCenterY();
    213  
    214   for (Int_t i=0; i<entries; i++ )
    215     {
    216       MSignalPix &pix = (*fSignalCam)[i];
    217      
    218       const MGeomPix &gpix = (*fGeomCam)[i/*pix.GetPixId()*/];
    219      
    220       const Float_t dx = gpix.GetX() - cenx;
    221       const Float_t dy = gpix.GetY() - ceny;
    222      
    223      const Float_t dist = TMath::Sqrt(dx*dx+dy*dy);
    224      
    225       Float_t ang = TMath::ACos(dx/dist);
    226       if(dy>0)
    227         ang *= -1.0;
    228      
    229       // if the signal is not near the estimated circle, it is ignored.
    230       if(dist < (*fMuonSearchPar).GetRadius() + fMuonSetup->GetMargin()
    231          && dist > (*fMuonSearchPar).GetRadius() - fMuonSetup->GetMargin())
    232         {
    233           // check whether ummapped pixel is used or not.
    234           // if it is so, ingnore the pixel information since the pixels totally deteriorate the muon information.
    235           if(pix.IsPixelUnmapped())
    236             {
    237               fMuonCalibPar->SetUseUnmap();
    238               continue;
    239             }
    240           fMuonCalibPar->fHistPhi->Fill(ang*kRad2Deg, pix.GetNumPhotons());
    241           MuonSize += pix.GetNumPhotons();
    242         }
    243 
    244       // use only the inner pixles. This is geometry dependent. This has to
    245       // be fixed!
    246       if(i>397)
    247         continue; 
    248      
    249       fMuonCalibPar->fHistWidth
    250         ->Fill(dist*(*fGeomCam).GetConvMm2Deg(), pix.GetNumPhotons());
    251     }
    252 
    253   fMuonCalibPar->SetMuonSize(MuonSize);
    254 
    255   // error estimation (temporaly)
    256   //  The error is estimated from the signal. In order to do so, we have to
    257   // once convert the signal from ADC to photo-electron. Then we can get
    258   // the fluctuation such as F-factor*sqrt(phe).
    259   //  Up to now, the error of pedestal is not taken into accout. This is not
    260   // of course correct. We will include this soon.
    261     Double_t ADC2PhEl = 0.14;
    262     Double_t Ffactor = 1.4;
    263     for(Int_t i=0; i<binnumphi+1; i++)
    264       {
    265         Float_t rougherr  = TMath::Sqrt(TMath::Abs(fMuonCalibPar->fHistPhi->GetBinContent(i))*ADC2PhEl)/ADC2PhEl*Ffactor;
    266         {
    267           fMuonCalibPar->fHistPhi->SetBinError(i, rougherr);
    268         }
    269       }
    270     for(Int_t i=0; i<binnumwid+1; i++)
    271       {
    272         Float_t rougherr = TMath::Sqrt(TMath::Abs(fMuonCalibPar->fHistWidth->GetBinContent(i))*ADC2PhEl)/ADC2PhEl*Ffactor;
    273         {
    274           fMuonCalibPar->fHistWidth->SetBinError(i, rougherr);
    275         }
    276       }
    277 }
    278 
    279 // --------------------------------------------------------------------------
    280 //
    281 //  Photon distribution along the estimated circle is fitted with theoritical
    282 // function in order to get some more information such as Arc Phi and Arc
    283 // Length.
    284 //
    285 void MMuonCalibParCalc::CalcPhi()
    286 {
    287   Float_t thres = fMuonSetup->GetArcPhiThres();
    288   Float_t startval = fMuonSetup->fArcPhiHistStartVal;
    289   Float_t endval = fMuonSetup->fArcPhiHistEndVal;
    290   Int_t   binnum = fMuonSetup->fArcPhiBinNum;
    291 
    292   Float_t convbin2val = (endval-startval)/(Float_t)binnum;
    293 
    294     // adjust the peak to 0
    295     Float_t maxval = 0.;
    296     Int_t   maxbin = 0;
    297     maxval = fMuonCalibPar->fHistPhi->GetMaximum();
    298     maxbin = fMuonCalibPar->fHistPhi->GetMaximumBin();
    299     fMuonCalibPar->SetPeakPhi(180.-(Float_t)(maxbin-1.)*convbin2val);
    300     TArrayD tmp;
    301     tmp.Set(binnum+1);
    302     for(Int_t i=1; i<binnum+1; i++)
    303       {
    304         tmp[i] = fMuonCalibPar->fHistPhi->GetBinContent(i);
    305       }
    306     for(Int_t i=1; i<binnum+1; i++)
    307       {
    308         Int_t id;
    309         id = i + (maxbin-(Int_t)((Float_t)binnum/2.)-1);
    310         if(id>binnum)
    311           {
    312             id-=(binnum);
    313           }
    314         if(id<=0)
    315           {
    316             id+=(binnum);
    317           }
    318         fMuonCalibPar->fHistPhi->SetBinContent(i,tmp[id]);
    319       }
    320     maxbin = (Int_t)((Float_t)binnum/2.)+1;
    321 
    322   // Determination of fitting region
    323   // The threshold is fixed with 100 [photons or ADC] in a bin. Therefore,
    324   // if you change the bin number, YOU HAVE TO CHANGE THIS VALUE!!!
    325   Float_t startfitval = 0.;
    326   Float_t endfitval   = 0.;
    327   Bool_t  IsInMaxim = kFALSE;
    328   Int_t effbinnum = 0;
    329   for(Int_t i=1; i<binnum+1; i++)
    330     {
    331       Float_t content = fMuonCalibPar->fHistPhi->GetBinContent(i);
    332       Float_t content_pre = fMuonCalibPar->fHistPhi->GetBinContent(i-1);
    333      
    334       if(content > thres && content_pre < thres)
    335         {
    336           startfitval = (Float_t)(i-1)*convbin2val+startval;
    337         }
    338       if(i==maxbin)
    339         IsInMaxim = kTRUE;
    340      
    341       if(content < thres && IsInMaxim == kTRUE)
    342         {
    343           endfitval = (Float_t)(i-1)*convbin2val+startval;
    344           break;
    345         }
    346       endfitval = endval;
    347     }
    348  
    349   effbinnum = (Int_t)((endfitval-startfitval)/convbin2val);
    350 
    351   fMuonCalibPar->SetArcPhi(endfitval-startfitval);
    352 
    353   fMuonCalibPar->SetArcLength( fMuonCalibPar->GetArcPhi()*TMath::DegToRad()*(*fMuonSearchPar).GetRadius()*(*fGeomCam).GetConvMm2Deg());
    354  
    355   if(fEnableImpactCalc)
    356       CalcImpact(effbinnum, startfitval, endfitval);
    357 }
    358 
    359 // --------------------------------------------------------------------------
    360 //
    361 //  An impact parameter is calculated by fitting the histogram of photon
    362 // distribution along the circle with a theoritical model.
    363 // (See G. Vacanti et. al., Astroparticle Physics 2, 1994, 1-11.
    364 // The function (6) is used here.)
    365 //
    366 //  By default this calculation is suppressed because this calculation is
    367 // very time consuming. If you want to calculate an impact parameter,
    368 // you can call the function of EnableImpactCalc().
    369 //
    370 void MMuonCalibParCalc::CalcImpact
    371 (Int_t effbinnum, Float_t startfitval, Float_t endfitval)
    372 {
    373   // Fit the distribution with Vacanti function. The function is different
    374   // for the impact parameter of inside or outside of our reflector.
    375   // Then two different functions are applied to the photon distribution,
    376   // and the one which give us smaller chisquare value is taken as a
    377   // proper one.
    378   Double_t val1,err1,val2,err2;
    379   // impact parameter inside mirror radius (8.5m)
    380   TString func1;
    381   Float_t tmpval = (*fMuonSearchPar).GetRadius()*(*fGeomCam).GetConvMm2Deg()*TMath::DegToRad();
    382   tmpval = sin(2.*tmpval)*8.5;
    383   func1 += "[0]*";
    384   func1 += tmpval;
    385   func1 += "*(sqrt(1.-([1]/8.5)**2*sin((x-[2])*3.1415926/180.)**2)+([1]/8.5)*cos((x-[2])*3.1415926/180.))";
    386  
    387   TF1 f1("f1",func1,startfitval,endfitval);
    388   f1.SetParameters(2000,3,0);
    389   f1.SetParLimits(1,0,8.5);
    390   f1.SetParLimits(2,-180.,180.);
    391  
    392   fMuonCalibPar->fHistPhi->Fit("f1","RQEM");
    393  
    394   Float_t chi1 = -1;
    395   Float_t chi2 = -1;
    396   if(effbinnum>3)
    397     chi1 = f1.GetChisquare()/((Float_t)(effbinnum-3));
    398  
    399   gMinuit->GetParameter(1,val1,err1);  // get the estimated IP
    400   Float_t estip1 = val1;
    401  
    402   // impact parameter beyond mirror area (8.5m)
    403   TString func2;
    404   Float_t tmpval2 = (*fMuonSearchPar).GetRadius()*(*fGeomCam).GetConvMm2Deg()*TMath::DegToRad();
    405   tmpval2 = sin(2.*tmpval2)*8.5*2.;
    406   func2 += "[0]*";
    407   func2 += tmpval2;
    408   func2 += "*sqrt(1.-(([1]/8.5)*sin((x-[2])*3.1415926/180.))**2)";
    409  
    410   TF1 f2("f2",func2,startfitval,endfitval);
    411   f2.SetParameters(2000,20,0);
    412   f2.SetParLimits(1,8.5,300.);
    413   f2.SetParLimits(2,-180.,180.);
    414  
    415   fMuonCalibPar->fHistPhi->Fit("f2","RQEM+");
    416  
    417   if(effbinnum>3)
    418     chi2 = f2.GetChisquare()/((Float_t)(effbinnum-3));
    419  
    420   gMinuit->GetParameter(1,val2,err2);  // get the estimated IP
    421   Float_t estip2 = val2;
    422  
    423   if(effbinnum<=3)
    424     {
    425       fMuonCalibPar->SetEstImpact(-1.);
    426     }
    427   if(chi2 > chi1)
    428     {
    429       fMuonCalibPar->SetEstImpact(estip1);
    430       fMuonCalibPar->SetChiArcPhi(chi1);
    431     }
    432   else
    433     {
    434       fMuonCalibPar->SetEstImpact(estip2);
    435       fMuonCalibPar->SetChiArcPhi(chi2);
    436     }
    437 }
    438 
    439 // --------------------------------------------------------------------------
    440 //
    441 //  Photon distribution of distance from the center of estimated ring is
    442 // fitted in order to get some more information such as ARC WIDTH which
    443 // can represent to the PSF of our reflector.
    444 //
    445 Float_t MMuonCalibParCalc::CalcWidth()
    446 {
    447   Float_t startval = fMuonSetup->fArcWidthHistStartVal;
    448   Float_t endval = fMuonSetup->fArcWidthHistEndVal;
    449   Int_t   binnum = fMuonSetup->fArcWidthBinNum;
    450   Float_t thres = fMuonSetup->GetArcWidthThres();
    451 
    452   Float_t convbin2val = (endval - startval)
    453     /(Float_t)binnum;
    454 
    455     // determination of fitting region
    456     Int_t maxbin = fMuonCalibPar->fHistWidth->GetMaximumBin();
    457     Float_t startfitval = 0.;
    458     Float_t endfitval   = 0.;
    459     Bool_t   IsInMaxim = kFALSE;
    460     Int_t    effbinnum = 0;
    461     for(Int_t i=1; i<binnum+1; i++)
    462       {
    463         Float_t content = fMuonCalibPar->fHistWidth->GetBinContent(i);
    464         Float_t content_pre = fMuonCalibPar->fHistWidth->GetBinContent(i-1);
    465 
    466         if(content > thres)
    467           effbinnum++;
    468 
    469         if(content > thres && content_pre < thres)
    470           {
    471             startfitval = (Float_t)(i-4)*convbin2val + startval;
    472             if(startfitval<0.) startfitval = 0.;
    473           }
    474         if(i==maxbin)
    475           IsInMaxim = kTRUE;
    476 
    477         if(content < thres && IsInMaxim == kTRUE)
    478           {
    479             endfitval = (Float_t)(i+2)*convbin2val + startval;
    480             if(endfitval>180.) endfitval = 180.;
    481             break;
    482           }
    483         endfitval = endval;
    484       }
    485     effbinnum = (Int_t)((endfitval-startfitval)/convbin2val);
    486 
    487     TF1 f1("f1","gaus",startfitval,endfitval);
    488 
    489     // options : N  do not store the function, do not draw
    490     //           I  use integral of function in bin rather than value at bin center
    491     //           R  use the range specified in the function range
    492     //           Q  quiet mode
    493     fMuonCalibPar->fHistWidth->Fit("f1","QR","",startfitval,endfitval);
    494    
    495     if(effbinnum>3)
    496       fMuonCalibPar->SetChiArcWidth(f1.GetChisquare()/((Float_t)(effbinnum-3)));
    497 
    498     Double_t val,err;
    499     gMinuit->GetParameter(2,val,err); // get the sigma value
    500 
    501     return val;
    502 }
    503 
    504 // --------------------------------------------------------------------------
    505 //
    506 //  Calculation of muon parameters
    507 //
    508 Int_t MMuonCalibParCalc::Calc()
    509 {
    510   // initialization
    511   (*fMuonCalibPar).Reset();
    512 
    513   // Fill signals to histograms
    514   FillHist();
    515 
    516   // Calculation of Arc Phi etc...
    517   CalcPhi();
    518 
    519   // Calculation of Arc Width etc...
    520   fMuonCalibPar->SetArcWidth(CalcWidth());
    521 
    522   if(fMuonCalibPar->GetArcPhi()>160 && fMuonSearchPar->GetRadius()>170 &&
    523      fMuonSearchPar->GetRadius()<400 && fMuonSearchPar->GetDeviation()<50)
    524       fMuonCalibPar->SetReadyToSave();
    525  
    526   return kTRUE;
    527 }
    528183
    529184// -------------------------------------------------------------------------
     
    531186Int_t MMuonCalibParCalc::Process()
    532187{
    533 
    534   Calc();
    535 
    536   return kTRUE;
     188    // Calculation of Arc Phi etc...
     189    const Float_t thresphi   = fMuonSetup->GetThresholdArcPhi()  /fGeomCam->GetConvMm2Deg();
     190    const Float_t threswidth = fMuonSetup->GetThresholdArcWidth()/fGeomCam->GetConvMm2Deg();
     191
     192    Double_t peakphi, arcphi;
     193    Double_t width, chi;
     194
     195    if (!fHist->CalcPhi(thresphi, peakphi, arcphi))
     196        return kTRUE;
     197
     198    if (!fHist->CalcWidth(threswidth, width, chi))
     199        return kTRUE;
     200
     201    // Get Muon Size
     202    fMuonCalibPar->SetMuonSize(fHist->GetHistPhi().Integral());
     203
     204    // Calculate Arc Length
     205    fMuonCalibPar->SetPeakPhi(peakphi);
     206    fMuonCalibPar->SetArcPhi(arcphi);
     207
     208    const Float_t conv = TMath::RadToDeg()*fGeomCam->GetConvMm2Deg();
     209    fMuonCalibPar->SetArcLength(fMuonCalibPar->GetArcPhi()*fMuonSearchPar->GetRadius()*conv);
     210
     211    // Calculation of Arc Width etc...
     212    fMuonCalibPar->SetChiArcWidth(chi);
     213    fMuonCalibPar->SetArcWidth(width);
     214
     215    // Check if this is a 'Write-Out' candidate
     216    if (fMuonCalibPar->GetArcPhi()>160    && fMuonSearchPar->GetRadius()<400 &&
     217        fMuonSearchPar->GetDeviation()<50 && fMuonSearchPar->GetRadius()>170)
     218        fMuonCalibPar->SetReadyToSave();
     219
     220    return kTRUE;
    537221}
  • trunk/MagicSoft/Mars/mmuon/MMuonCalibParCalc.h

    r6973 r6979  
    88class MMuonSearchPar;
    99class MMuonCalibPar;
    10 class MSrcPosCam;
    1110class MGeomCam;
    12 class MSignalCam;
    1311class MMuonSetup;
     12class MHSingleMuon;
    1413
    1514class MMuonCalibParCalc : public MTask
    1615{
    1716private:
    18     MGeomCam       *fGeomCam;
    19     MSignalCam     *fSignalCam;
    20     MMuonCalibPar  *fMuonCalibPar; 
    21     MMuonSearchPar *fMuonSearchPar;
    22     MMuonSetup     *fMuonSetup;
     17    MGeomCam       *fGeomCam;         //!
     18    MMuonCalibPar  *fMuonCalibPar;    //!
     19    MMuonSearchPar *fMuonSearchPar;   //!
     20    MMuonSetup     *fMuonSetup;       //!
     21    MHSingleMuon   *fHist;            //!
    2322
    24     Bool_t fEnableImpactCalc; // If true, the impact calculation will be done, which consumes a lot of time.
     23    //Bool_t fEnableImpactCalc; // If true, the impact calculation will be done, which consumes a lot of time.
    2524
    26     Int_t PreProcess(MParList *plist);
    27     Int_t Process();
     25    Int_t   PreProcess(MParList *plist);
     26    Int_t   Process();
     27
     28    void    FillHist();
     29    void    CalcPhi();
     30    void    CalcImpact(Int_t effbinnum, Float_t startfitval, Float_t endfitval);
     31    Float_t CalcWidth();
    2832
    2933public:
    3034    MMuonCalibParCalc(const char *name=NULL, const char *title=NULL);
    3135
    32     void EnableImpactCalc()              { fEnableImpactCalc = kTRUE; }
    33 
    34     void FillHist();
    35     void CalcPhi();
    36     void CalcImpact(Int_t effbinnum, Float_t startfitval, Float_t endfitval);
    37     Float_t CalcWidth();
    38     Int_t   Calc();
     36    //void EnableImpactCalc(Bool_t b=kTRUE) { fEnableImpactCalc = b; }
    3937
    4038    ClassDef(MMuonCalibParCalc, 0) // task to calculate muon parameters
  • trunk/MagicSoft/Mars/mmuon/MMuonSearchPar.cc

    r6973 r6979  
    1717!
    1818!   Author(s): Keiichi Mase 10/2004 <mailto:mase@mppmu.mpg.de>
    19 !             Markus Meyer 10/2004 <mailto:meyer@astro.uni-wuerzburg.de>
    20 !
    21 !   Copyright: MAGIC Software Development, 2000-2004
     19!   Author(s): Markus Meyer 10/2004 <mailto:meyer@astro.uni-wuerzburg.de>
     20!
     21!   Copyright: MAGIC Software Development, 2000-2005
    2222!
    2323!
     
    4040// MMuonCalibPar. The information will be available by using the task of
    4141// MMuonCalibParCalc.
    42 //
    4342//
    4443// --- How to search muons ---
     
    4746//
    4847// 1. A temporal center position of a circle is determined by using
    49 //  the Hillas parameters. Assumed that the center position will be on the
    50 //  line which is perpendicular to the longitudinal image axis and the
    51 //  distance from the gravity center of the image to the center position of
    52 //  a ring is approximately 1 deg. (corresponding to the Cherenkov angle.).
    53 //  Therefore, we will have two candidates of the center positions.
     48//    the Hillas parameters. Assumed that the center position will be on the
     49//    line which is perpendicular to the longitudinal image axis and the
     50//    distance from the gravity center of the image to the center position of
     51//    a ring is approximately 1 deg. (corresponding to the Cherenkov angle.).
     52//    Therefore, we will have two candidates of the center positions.
     53//
    5454// 2. Find the ring radius which gives the minimum RMS between the camera
    55 //  images and the estimated circle.
     55//    images and the estimated circle.
     56//
    5657// 3. Select one temporal position which gives smaller RMS as a true temporal
    57 //  center position.
     58//    center position.
     59//
    5860// 4. Changing the center position of a circle on the camera plane from the
    59 // determined temporal center position, find the position which gives the
    60 // minimum RMS of the fit.
    61 //
    62 //
    63 // --- Remark ---
    64 //  This method to search for muons is not fully optimized yet. However,
    65 // it is good idea to use the temporal estimated center position from
    66 // the Hillas parameters in order to reduce time to search. In addition,
    67 // This method is faster than the MINUIT.
     61//    determined temporal center position, find the position which gives the
     62//    minimum RMS of the fit.
    6863//
    6964//
    7065//  Input Containers:
    71 //   [MGeomCam]
    72 //   [MHillas]
    73 //   [MSignalCam]
     66//   MGeomCam
     67//   MHillas
     68//   MSignalCam
    7469//
    7570/////////////////////////////////////////////////////////////////////////////
    7671#include "MMuonSearchPar.h"
    7772
    78 #include <fstream>
     73#include <TMinuit.h>
    7974
    8075#include "MLog.h"
    8176#include "MLogManip.h"
     77
    8278#include "MHillas.h"
     79
    8380#include "MGeomCam.h"
    8481#include "MGeomPix.h"
     82
    8583#include "MSignalPix.h"
    8684#include "MSignalCam.h"
     
    9694MMuonSearchPar::MMuonSearchPar(const char *name, const char *title)
    9795{
    98   fName  = name  ? name  : "MMuonSearchPar";
    99   fTitle = title ? title : "Muon search parameters";
     96    fName  = name  ? name  : "MMuonSearchPar";
     97    fTitle = title ? title : "Muon search parameters";
    10098}
    10199
     
    104102void MMuonSearchPar::Reset()
    105103{
    106   fRadius  = -1.;
    107   fDeviation  = -1.;
    108   fCenterX = 0.;
    109   fCenterY = 0.;
    110 }
    111 
    112 // --------------------------------------------------------------------------
    113 //
    114 //  Get the tempolary center of a ring from the Hillas parameters.
    115 //  Two candidates of the position is returened.
    116 //
    117 void MMuonSearchPar::CalcTempCenter(const MHillas &hillas,
    118       Float_t &xtmp1, Float_t &ytmp1, Float_t &xtmp2, Float_t &ytmp2)
    119 {
    120   Float_t a,dx,dy;
    121   Float_t tmp_r = 300.;  // assume that the temporal cherenkov angle is 1 deg. (300 mm).
    122 
    123   a = TMath::Tan(hillas.GetDelta());
    124 
    125   dx = a/TMath::Sqrt(tmp_r+a*a)/3.;
    126   dy = -tmp_r/TMath::Sqrt(1+a*a)/3.;
    127 
    128   xtmp1 = hillas.GetMeanX() + dx;
    129   ytmp1 = hillas.GetMeanY() + dy;
    130   xtmp2 = hillas.GetMeanX() - dx;
    131   ytmp2 = hillas.GetMeanY() - dy;
     104    fRadius    = -1;
     105    fDeviation = -1;
     106    fCenterX   =  0;
     107    fCenterY   =  0;
     108}
     109
     110// --------------------------------------------------------------------------
     111//
     112// This is a wrapper function to have direct access to the data members
     113// in the function calculating the minimization value.
     114//
     115void MMuonSearchPar::fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
     116{
     117    const MMuonSearchPar *optim = (MMuonSearchPar*)gMinuit->GetObjectFit();
     118    f = optim->Fcn(par);
    132119}
    133120
     
    137124//  and its RMS for the input position.
    138125//
    139 Bool_t MMuonSearchPar::CalcRadius(const MGeomCam &geom, const MSignalCam &evt,
    140       Float_t x, Float_t y, Float_t &r, Float_t &sigma)
    141 {
    142   Float_t mean_r=0., dev_r=0., sums=0., tmp=0.;
    143 
    144   const Int_t entries = evt.GetNumPixels();
    145 
    146   for (Int_t i=0; i<entries; i++ ){
    147     const MSignalPix &pix = evt[i];
    148 
    149     if (!pix.IsPixelUsed())
    150       continue;
    151 
    152     const MGeomPix &gpix = geom[i/*pix.GetPixId()*/];
    153 
    154     tmp=TMath::Sqrt((gpix.GetX()-x)*(gpix.GetX()-x)
    155                     +(gpix.GetY()-y)*(gpix.GetY()-y));
    156 
    157     mean_r += pix.GetNumPhotons()*tmp;
    158     dev_r  += pix.GetNumPhotons()*tmp*tmp;
    159     sums   += pix.GetNumPhotons();
    160   }
    161 
    162   if(sums<1.E-10)
    163       return kFALSE;
    164 
    165   r = mean_r/sums;
    166 
    167   if(dev_r/sums-(r)*(r)<1.E-10)
    168       return kFALSE;
    169 
    170   sigma = TMath::Sqrt(dev_r/sums-(r)*(r));
    171 
    172   return kTRUE;
    173 }
     126Double_t MMuonSearchPar::Fcn(Double_t *par) const
     127{
     128    const Int_t entries = fSignal.GetSize();
     129
     130    Double_t meanr=0;
     131    Double_t devr =0;
     132    Double_t sums =0;
     133
     134    // It seems that the loop is easy enough for a compiler optimization.
     135    // Using pointer arithmetics doesn't improve the speed of the fit.
     136    for (Int_t i=0; i<entries; i++ )
     137    {
     138        Double_t tmp = TMath::Hypot(fX[i]-par[0], fY[i]-par[1]);
     139
     140        sums  += fSignal[i];
     141        meanr += fSignal[i] * tmp;
     142        devr  += fSignal[i] * tmp*tmp;
     143    }
     144
     145    par[2] = meanr/sums;
     146    par[3] = devr/sums - par[2]*par[2];
     147
     148    return par[3];
     149}
    174150
    175151// --------------------------------------------------------------------------
     
    178154// RMS of the fit, changing the center position of the circle.
    179155//
    180 void MMuonSearchPar::CalcMinimumDeviation
    181 ( const MGeomCam &geom, const MSignalCam &evt, Float_t x, Float_t y,
    182  Float_t xcog, Float_t ycog, Float_t sigma, Float_t &opt_rad,
    183  Float_t &new_sigma, Float_t &newx, Float_t &newy )
    184 {
    185   Float_t delta = 3.;  // 3 mm (1/10 of an inner pixel size) Step to move.
    186   Float_t rad_tmp,sig_tmp;
    187   Float_t r2;
    188 
    189   while(1)
    190   {
    191       r2=(xcog-x)*(xcog-x)+(ycog-y)*(ycog-y);
    192       // Exit if the new estimated radius is above 2 deg. (600 mm).
    193       if(r2 > 360000.)
    194       {
    195           new_sigma=sigma;
    196           opt_rad=rad_tmp;
    197           newx=x;
    198           newy=y;
    199           break;
    200       }
    201       if(CalcRadius(geom,evt,x,y+delta,rad_tmp,sig_tmp))
    202       {
    203           if(sig_tmp<sigma)
    204           {
    205               sigma=sig_tmp;
    206               opt_rad=rad_tmp;
    207               y=y+delta;
    208               continue;
    209           }
    210       }
    211       if(CalcRadius(geom,evt,x-delta,y,rad_tmp,sig_tmp))
    212       {
    213           if(sig_tmp<sigma)
    214           {
    215               sigma=sig_tmp;
    216               opt_rad=rad_tmp;
    217               x=x-delta;
    218               continue;
    219           }
    220       }
    221       if(CalcRadius(geom,evt,x+delta,y,rad_tmp,sig_tmp))
    222       {
    223           if(sig_tmp<sigma)
    224           {
    225               sigma=sig_tmp;
    226               opt_rad=rad_tmp;
    227               x=x+delta;
    228               continue;
    229           }
    230       }
    231       if(CalcRadius(geom,evt,x,y-delta,rad_tmp,sig_tmp))
    232       {
    233           if(sig_tmp<sigma)
    234           {
    235               sigma=sig_tmp;
    236               opt_rad=rad_tmp;
    237               y=y-delta;
    238               continue;
    239           }
    240       }
    241       new_sigma=sigma;
    242       newx=x;
    243       newy=y;
    244       break;
     156void MMuonSearchPar::CalcMinimumDeviation(const MGeomCam &geom,
     157                                          const MSignalCam &evt,
     158                                          Double_t &x, Double_t &y,
     159                                          Double_t &sigma, Double_t &radius)
     160{
     161    // ------- Make a temporaray copy of the signal ---------
     162    const Int_t n = geom.GetNumPixels();
     163
     164    fSignal.Set(n);
     165    fX.Set(n);
     166    fY.Set(n);
     167
     168    Int_t p=0;
     169    for (int i=0; i<n; i++)
     170    {
     171        const MSignalPix &pix = evt[i];
     172        if (pix.IsPixelUsed())
     173        {
     174            fSignal[p] = pix.GetNumPhotons();
     175
     176            fX[p] = geom[i].GetX();
     177            fY[p] = geom[i].GetY();
     178            p++;
     179        }
    245180    }
     181    fSignal.Set(p);
     182
     183
     184    // ----------------- Setup and call minuit -------------------
     185    const Float_t  delta = 30.;  // 3 mm (1/10 of an inner pixel size) Step to move.
     186    const Double_t r     = geom.GetMaxRadius()*2;
     187
     188    // Save gMinuit
     189    TMinuit *minsave = gMinuit;
     190
     191    // Initialize TMinuit with 4 parameters
     192    TMinuit minuit;
     193    minuit.SetPrintLevel(-1);     // Switch off printing
     194    minuit.Command("set nowarn"); // Switch off warning
     195    // Define Parameters
     196    minuit.DefineParameter(0, "x",     x, delta, -r, r);
     197    minuit.DefineParameter(1, "y",     y, delta, -r, r);
     198    minuit.DefineParameter(2, "r",     0, 1,      0, 0);
     199    minuit.DefineParameter(3, "sigma", 0, 1,      0, 0);
     200    // Fix return parameters
     201    minuit.FixParameter(2);
     202    minuit.FixParameter(3);
     203    // Setup Minuit for 'this' and use fit function fcn
     204    minuit.SetObjectFit(this);
     205    minuit.SetFCN(fcn);
     206
     207    // Perform Simplex minimization
     208    Int_t err;
     209    Double_t tmp[2] = { 0, 0 };
     210    minuit.mnexcm("simplex", tmp, 2, err);
     211
     212    // Get resulting parameters
     213    Double_t error;
     214    minuit.GetParameter(0, x,      error);
     215    minuit.GetParameter(1, y,      error);
     216    minuit.GetParameter(2, radius, error);
     217    minuit.GetParameter(3, sigma,  error);
     218
     219    sigma = sigma>0 ? TMath::Sqrt(sigma) : 0;
     220
     221    gMinuit = minsave;
    246222}
    247223
     
    250226//  Calculation of muon parameters
    251227//
    252 void MMuonSearchPar::Calc
    253 (const MGeomCam &geom, const MSignalCam &evt, const MHillas &hillas)
    254 {
    255   Reset();
    256  
    257   Float_t xtmp1,xtmp2,ytmp1,ytmp2;
    258   Float_t rad,dev,rad2,dev2;
    259   Float_t opt_rad,new_sigma,newx,newy;
    260    
    261   // gets temporaly center
    262   CalcTempCenter(hillas,xtmp1,ytmp1,xtmp2,ytmp2);
    263  
    264   // determine which position will be the true position. Here mainly
    265   // the curvature of a muon arc is relied on.
    266   CalcRadius(geom, evt, xtmp1,ytmp1,rad,dev);
    267   CalcRadius(geom, evt, xtmp2,ytmp2,rad2,dev2);
    268   if(dev2<dev){
    269     xtmp1=xtmp2; ytmp1=ytmp2; dev=dev2; rad=rad2;
    270   }
    271  
    272   // find the best fit.
    273   CalcMinimumDeviation(geom, evt, xtmp1,ytmp1,hillas.GetMeanX(),
    274                        hillas.GetMeanY(), dev, opt_rad, new_sigma,
    275                        newx, newy);
    276 
    277   fRadius = opt_rad;
    278   fDeviation = new_sigma;
    279  
    280   fCenterX = newx;
    281   fCenterY = newy;
    282  
    283   //SetReadyToSave();
     228void MMuonSearchPar::Calc(const MGeomCam &geom, const MSignalCam &evt,
     229                          const MHillas &hillas)
     230{
     231    Double_t x = hillas.GetMeanX();
     232    Double_t y = hillas.GetMeanY();
     233
     234    // -------------------------------------------------
     235    // Keiichi suggested trying to precalculate the Muon
     236    // center a bit better, but it neither improves the
     237    // fit result nor the speed
     238    //
     239    // const Float_t tmpr = 300.;  // assume that the temporal cherenkov angle is 1 deg. (300 mm).
     240    //
     241    // const Double_t a = TMath::Tan(hillas.GetDelta());
     242    //
     243    // const Double_t dx =     a/TMath::Sqrt(tmpr+a*a)/3.;
     244    // const Double_t dy = -tmpr/TMath::Sqrt(1+a*a)/3.;
     245    //
     246    // Double_t par1[] = { x+dx, y+dy, 0, 0 };
     247    // Double_t par2[] = { x-dx, y-dy, 0, 0 };
     248    //
     249    // const Double_t dev1 = MMuonSearchPar::Fcn(par1);
     250    // const Double_t dev2 = MMuonSearchPar::Fcn(par2);
     251    //
     252    // if (dev1<dev2)
     253    // {
     254    //     x += dx;
     255    //     y += dy;
     256    // }
     257    // else
     258    // {
     259    //     x -= dx;
     260    //     y -= dy;
     261    // }
     262    // -------------------------------------------------
     263
     264    Double_t sigma, rad;
     265
     266    // find the best fit.
     267    CalcMinimumDeviation(geom, evt, x, y, sigma, rad);
     268
     269    fCenterX   = x;
     270    fCenterY   = y;
     271    fRadius    = rad;
     272    fDeviation = sigma;
     273
     274    //SetReadyToSave();
    284275}
    285276
     
    298289    *fLog << all;
    299290    *fLog << "Muon Parameters (" << GetName() << ")" << endl;
    300     *fLog << " - Est. Radius   [deg.]  = " << fRadius*geom.GetConvMm2Deg()   << endl;
    301     *fLog << " - Deviation     [deg.]  = " << fDeviation*geom.GetConvMm2Deg()   << endl;
    302     *fLog << " - Center Pos. X [deg.]  = " << fCenterX*geom.GetConvMm2Deg()  << endl;
    303     *fLog << " - Center Pos. Y [deg.]  = " << fCenterY*geom.GetConvMm2Deg()  << endl;
    304 }
    305 
    306 
    307 
     291    *fLog << " - Est. Radius   [deg] = " << fRadius*geom.GetConvMm2Deg()   << endl;
     292    *fLog << " - Deviation     [deg] = " << fDeviation*geom.GetConvMm2Deg()   << endl;
     293    *fLog << " - Center Pos. X [deg] = " << fCenterX*geom.GetConvMm2Deg()  << endl;
     294    *fLog << " - Center Pos. Y [deg] = " << fCenterY*geom.GetConvMm2Deg()  << endl;
     295}
  • trunk/MagicSoft/Mars/mmuon/MMuonSearchPar.h

    r6973 r6979  
    44#ifndef MARS_MParContainer
    55#include "MParContainer.h"
     6#endif
     7
     8#ifndef MARS_MArrayF
     9#include "MArrayF.h"
    610#endif
    711
     
    1822    Float_t fCenterY;   // An estimated center position in Y of the muon ring [mm]
    1923
     24    MArrayF fSignal;    //! Temporary storage for signal
     25    MArrayF fX;         //! Temporary storage for pixels X-position
     26    MArrayF fY;         //! Temporary storage for pixels Y-position
     27
     28    static void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
     29    Double_t Fcn(Double_t *par) const;
     30
    2031public:
    2132    MMuonSearchPar(const char *name=NULL, const char *title=NULL);
     
    2839    Float_t GetCenterY()   const { return fCenterY; }
    2940
    30     void   CalcTempCenter(const MHillas &hillas, Float_t &xtmp1,
    31                           Float_t &ytmp1, Float_t &xtmp2, Float_t &ytmp2);
    32     Bool_t CalcRadius(const MGeomCam &geom, const MSignalCam &evt, Float_t x,
    33                       Float_t y, Float_t &r, Float_t &sigma);
    3441    void   CalcMinimumDeviation(const MGeomCam &geom, const MSignalCam &evt,
    35                              Float_t x, Float_t y, Float_t xcog,
    36                              Float_t ycog, Float_t sigma, Float_t &opt_rad,
    37                              Float_t &new_sigma, Float_t &newx, Float_t &newy);
     42                                Double_t &x, Double_t &y, Double_t &sigma, Double_t &rad);
     43
    3844    void   Calc(const MGeomCam &geom, const MSignalCam &evt,
    3945                const MHillas &hillas);
  • trunk/MagicSoft/Mars/mmuon/MMuonSearchParCalc.cc

    r6973 r6979  
    1717!
    1818!   Author(s): Keiichi Mase 09/2004 <mailto:mase@mppmu.mpg.de>
    19 !             Markus Meyer 09/2004 <mailto:meyer@astro.uni-wuerzburg.de>
     19!   Author(s): Markus Meyer 09/2004 <mailto:meyer@astro.uni-wuerzburg.de>
    2020!
    21 !   Copyright: MAGIC Software Development, 2000-2004
     21!   Copyright: MAGIC Software Development, 2000-2005
    2222!
    2323!
     
    4141//
    4242//  Input Containers:
    43 //   [MGeomCam]
    44 //   [MHillas]
    45 //   [MSignalCam]
     43//   MGeomCam
     44//   MHillas
     45//   MSignalCam
    4646//
    4747//  Output Containers:
    48 //   [MMuonSearchPar]
     48//   MMuonSearchPar
    4949//
    5050//////////////////////////////////////////////////////////////////////////////
    5151#include "MMuonSearchParCalc.h"
    5252
    53 #include <fstream>
     53#include "MParList.h"
    5454
    55 #include "MParList.h"
     55#include "MLog.h"
     56#include "MLogManip.h"
     57
    5658#include "MGeomCam.h"
    5759#include "MSignalCam.h"
    5860#include "MMuonSearchPar.h"
    59 #include "MLog.h"
    60 #include "MLogManip.h"
    6161
    6262using namespace std;
     
    7171// Default constructor.
    7272//
    73 MMuonSearchParCalc::MMuonSearchParCalc
    74 (const char *mupar, const char *name, const char *title)
     73MMuonSearchParCalc::MMuonSearchParCalc(const char *name, const char *title)
    7574  : fHillas(NULL), fMuonPar(NULL)
    7675{
    7776    fName  = name  ? name  : gsDefName.Data();
    7877    fTitle = title ? title : gsDefTitle.Data();
    79 
    80     fMuparName  = mupar;
    81     fHillasInput = "MHillas";
    8278}
    8379
     
    8682Int_t MMuonSearchParCalc::PreProcess(MParList *pList)
    8783{
    88     fHillas = (MHillas*)pList->FindObject(fHillasInput, "MHillas");
     84    fHillas = (MHillas*)pList->FindObject("MHillas");
    8985    if (!fHillas)
    9086    {
    91         *fLog << err << fHillasInput << " [MHillas] not found... aborting." << endl;
     87        *fLog << err << "MHillas not found... aborting." << endl;
    9288        return kFALSE;
    9389    }
     
    107103    }
    108104
    109     fMuonPar = (MMuonSearchPar*)pList->FindCreateObj("MMuonSearchPar", fMuparName);
     105    fMuonPar = (MMuonSearchPar*)pList->FindCreateObj("MMuonSearchPar");
    110106    if (!fMuonPar)
    111107        return kFALSE;
     
    118114Int_t MMuonSearchParCalc::Process()
    119115{
    120 
    121   fMuonPar->Calc(*fGeomCam, *fSignalCam, *fHillas);
    122 
    123   return kTRUE;
     116    fMuonPar->Calc(*fGeomCam, *fSignalCam, *fHillas);
     117    return kTRUE;
    124118}
    125 
  • trunk/MagicSoft/Mars/mmuon/MMuonSearchParCalc.h

    r6973 r6979  
    1414{
    1515private:
    16     MGeomCam    *fGeomCam;
    17     MSignalCam  *fSignalCam;
    18     MHillas     *fHillas;      //! Pointer to the source independent hillas parameters
    19     MMuonSearchPar *fMuonPar;  //! Pointer to the output container for the new image parameters
    20 
    21     TString fMuparName;
    22     TString fHillasInput;
     16    MGeomCam       *fGeomCam;     //!
     17    MSignalCam     *fSignalCam;   //!
     18    MHillas        *fHillas;      //! Pointer to the source independent hillas parameters
     19    MMuonSearchPar *fMuonPar;     //! Pointer to the output container for the new image parameters
    2320
    2421    Int_t PreProcess(MParList *plist);
     
    2623
    2724public:
    28     MMuonSearchParCalc(const char *mupar="MMuonSearchPar",
    29                        const char *name=NULL, const char *title=NULL);
    30 
    31     void SetInput(TString hilname) { fHillasInput = hilname; }
     25    MMuonSearchParCalc(const char *name=NULL, const char *title=NULL);
    3226
    3327    ClassDef(MMuonSearchParCalc, 0) // task to calculate muon parameters
  • trunk/MagicSoft/Mars/mmuon/MMuonSetup.cc

    r6973 r6979  
    1616!
    1717!
    18 !   Author(s): Keiichi Mase 10/2004 <mailto:mase@mppmu.mpg.de>
    19 !              Markus Meyer 10/2004 <mailto:meyer@astro.uni-wuerzburg.de>
     18!   Author(s): Markus Meyer 04/2005 <mailto:meyer@astro.uni-wuerzburg.de>
    2019!
    21 !   Copyright: MAGIC Software Development, 2000-2004
     20!   Copyright: MAGIC Software Development, 2000-2005
    2221!
    2322!
     
    2827// MMuonSetup
    2928//
    30 // Storage Container for parameter for the muon analysis
    31 //
    32 //
    33 //  Input Containers:
    34 //   [MSignalCam]
     29// Storage Container for setup parameter for the muon analysis
    3530//
    3631/////////////////////////////////////////////////////////////////////////////
     
    5045//
    5146MMuonSetup::MMuonSetup(const char *name, const char *title)
     47    : fMargin(0.2), fThresholdArcPhi(0.1), fThresholdArcWidth(0.1)
    5248{
    5349    fName  = name  ? name  : "MMuonSetup";
    54     fTitle = title ? title : "Muon calibration parameters";
    55 
    56     fMargin = 60.;  // in mm
    57     fArcPhiThres  = 30.;
    58     fArcWidthThres  = 30.;
    59     fArcPhiBinNum = 20;
    60     fArcPhiHistStartVal = -180.; // deg.
    61     fArcPhiHistEndVal   = 180.;  // deg.
    62     fArcWidthBinNum = 28;
    63     fArcWidthHistStartVal = 0.3; // deg.
    64     fArcWidthHistEndVal   = 1.7; // deg.
     50    fTitle = title ? title : "Muon calibration setup parameters";
    6551}
    6652
    6753// --------------------------------------------------------------------------
    6854//
    69 MMuonSetup::~MMuonSetup()
     55// Read the setup from a TEnv, eg:
     56//   MMuonSetup.Margin:            0.2
     57//   MMuonSetup.ThresholdArcPhi:   0.1
     58//   MMuonSetup.ThresholdArcWidth: 0.1
     59//
     60Int_t MMuonSetup::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
    7061{
     62    Bool_t rc = kFALSE;
     63    if (IsEnvDefined(env, prefix, "Margin", print))
     64    {
     65        rc = kTRUE;
     66        fMargin = GetEnvValue(env, prefix, "Margin", fMargin);
     67    }
     68    if (IsEnvDefined(env, prefix, "ThresholdArcPhi", print))
     69    {
     70        rc = kTRUE;
     71        fThresholdArcPhi = GetEnvValue(env, prefix, "ThresholdArcPhi", fThresholdArcPhi);
     72    }
     73    if (IsEnvDefined(env, prefix, "ThresholdArcWidth", print))
     74    {
     75        rc = kTRUE;
     76        fThresholdArcWidth = GetEnvValue(env, prefix, "ThresholdArcWidth", fThresholdArcWidth);
     77    }
     78
     79    return rc;
    7180}
  • trunk/MagicSoft/Mars/mmuon/MMuonSetup.h

    r6973 r6979  
    99{
    1010private:
    11   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
    12   Float_t fArcPhiThres;   // The threshold value to define arc phi
    13   Float_t fArcWidthThres; // The threshold value to define arc width
     11    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
     12    Float_t fThresholdArcPhi;   // [deg] The threshold value to define arc phi
     13    Float_t fThresholdArcWidth; // [deg] The threshold value to define arc width
    1414
     15    Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print);
    1516
    1617public:
    17   MMuonSetup(const char *name=NULL, const char *title=NULL);
    18   ~MMuonSetup();
    19  
    20   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.
    21   Int_t   fArcWidthBinNum;       // The bin number for the histogram of arc wid
    22   Float_t fArcPhiHistStartVal;   // The starting value for the histogram of arc phi
    23   Float_t fArcPhiHistEndVal;     // The end value for the histogram of arc phi
    24   Float_t fArcWidthHistStartVal; // The starting value for the histogram of arc width
    25   Float_t fArcWidthHistEndVal;   // The end value for the histogram of arc width
     18    MMuonSetup(const char *name=NULL, const char *title=NULL);
    2619
    27   Float_t GetMargin()         const { return fMargin; }
    28   Float_t GetArcPhiThres()    const { return fArcPhiThres; }
    29   Float_t GetArcWidthThres()  const { return fArcWidthThres; }
    30   Float_t GetArcPhiBinNum()   const { return fArcPhiBinNum; }
    31   Float_t GetArcWidthBinNum() const { return fArcWidthBinNum; }
    32  
    33   void    SetMargin(Float_t margin)       { fMargin = margin; }
    34   void    SetArcPhiThres(Float_t thres)   { fArcPhiThres = thres; }
    35   void    SetArcWidthThres(Float_t thres) { fArcWidthThres = thres; }
    36   void    SetArcPhiBinNum(Int_t num)      { fArcPhiBinNum = num; }
    37   void    SetArcWidthBinNum(Int_t num)    { fArcWidthBinNum = num; }
    38  
    39   ClassDef(MMuonSetup, 1) // Container to hold muon parameters for the whole sample
     20    // Getter
     21    Float_t GetMargin() const { return fMargin; }
     22    Float_t GetThresholdArcPhi() const { return fThresholdArcPhi; }
     23    Float_t GetThresholdArcWidth() const { return fThresholdArcWidth; }
     24
     25    // Setter
     26    void SetMargin(Float_t margin)           { fMargin = margin; }
     27    void SetThresholdArcPhi(Float_t thres)   { fThresholdArcPhi = thres; }
     28    void SetThresholdArcWidth(Float_t thres) { fThresholdArcWidth = thres; }
     29
     30    ClassDef(MMuonSetup, 1) // Container to hold setup for muon analysis
    4031};
    4132   
  • trunk/MagicSoft/Mars/mmuon/Makefile

    r6973 r6979  
    3030           MMuonCalibParCalc.cc \
    3131           MHMuonPar.cc \
     32           MHSingleMuon.cc \
    3233           MMuonSetup.cc
    3334
  • trunk/MagicSoft/Mars/mmuon/MuonLinkDef.h

    r6973 r6979  
    66
    77#pragma link C++ class MMuonSearchPar+;
    8 #pragma link C++ class MMuonSearchParCalc+;
    98#pragma link C++ class MMuonCalibPar+;
    10 #pragma link C++ class MMuonCalibParCalc+;
    11 #pragma link C++ class MHMuonPar+;
    129#pragma link C++ class MMuonSetup+;
    1310
     11#pragma link C++ class MMuonSearchParCalc+;
     12#pragma link C++ class MMuonCalibParCalc+;
     13
     14#pragma link C++ class MHMuonPar+;
     15#pragma link C++ class MHSingleMuon+;
     16
    1417#endif
    15 
    16 
    17 
    18 
    19 
    20 
    21 
Note: See TracChangeset for help on using the changeset viewer.