Changeset 5332 for trunk/MagicSoft/Mars


Ignore:
Timestamp:
11/01/04 22:31:44 (20 years ago)
Author:
mase
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/Changelog

    r5331 r5332  
    2323 2004/11/01: Markus Meyer and Keiichi Mase
    2424
    25    * mmuon/MMuonCalibParCalc.[h,cc],
     25   * mmuon/MMuonCalibParCalc.[h,cc]
    2626     - changed according to Thomas Bretz's suggestion.
    2727       The main part of the calculation is moved from MMuonCalibPar to
    … …  
    3131       that we should get the Arc Width with uncleaned images.
    3232
     33   * mmuon/MMuonCalibPar.[h,cc]
     34     - changed according to Thomas Bretz's suggestion as written above.
     35
    3336   * mmuon/MMuonSearchPar.[h,cc]
    3437     - changed according to Thomas Bretz's suggestion.
    3538       The pointer of the funciotn is changed to the reference.
    36 
    37    * mmuon/MMuonCalibPar.[h,cc],
    38      - changed according to Thomas Bretz's suggestion as written above.
    39 
    4039
    4140

  • trunk/MagicSoft/Mars/mmuon/MMuonCalibPar.cc

    r5210 r5332  
    4646#include <fstream>
    4747
    48 #include <TH1.h>
    49 #include <TF1.h>
    50 #include <TMinuit.h>
    51 
    5248#include "MLog.h"
    5349#include "MLogManip.h"
    54 #include "MGeomCam.h"
    55 #include "MGeomPix.h"
    5650#include "MCerPhotEvt.h"
    5751#include "MCerPhotPix.h"
    … …  
    9185    fHistWidth->UseCurrentStyle();
    9286
    93     fEnableImpactCalc = kFALSE; // By default the calculation of impact parameter is skipped.
    94     fDisablePreCuts = kFALSE;   // By default the pre cuts will be applied.
    95     fUseCleanForWidth = kFALSE; // By default all the pixels will be used for the histogram of arc width.
    96 
    9787    fMargin = 60.;  // in mm
    9888    fArcPhiThres  = 100.;
    … …  
    10292    fArcPhiHistEndVal   = 180.;  // deg.
    10393    fArcWidthBinNum = 28;
    104     fArcWidthHistStartVal = 0.3;   // deg.
    105     fArcWidthHistEndVal   = 1.7;   // deg.
     94    fArcWidthHistStartVal = 0.3; // deg.
     95    fArcWidthHistEndVal   = 1.7; // deg.
    10696}
    10797
    … …  
    127117  fUseUnmap    = kFALSE;
    128118  fPeakPhi     =  0.;
     119
     120  fHistPhi->Reset();
     121  fHistWidth->Reset();
    129122}
    130 
    131 // --------------------------------------------------------------------------
    132 //
    133 //  This function fill the histograms in order to get muon parameters.
    134 // For the evaluation of the Arc Width, we use only the signals in the inner
    135 // part. You can use the image after the cleaning by using the function of
    136 // UseCleanForWidth(). See the manual of MMuonCalibParCalc.
    137 //
    138 void MMuonCalibPar::FillHist
    139 ( const MGeomCam &geom, const MCerPhotEvt &evt,
    140   const MMuonSearchPar &musearch)
    141 {
    142   // preparation for a histgram
    143   MBinning binsphi;
    144   binsphi.SetEdges(fArcPhiBinNum, fArcPhiHistStartVal, fArcPhiHistEndVal);
    145   binsphi.Apply(*fHistPhi);
    146 
    147   MBinning binswid;
    148   binswid.SetEdges(fArcWidthBinNum, fArcWidthHistStartVal, fArcWidthHistEndVal);
    149   binswid.Apply(*fHistWidth);
    150 
    151  const Int_t entries = evt.GetNumPixels();
    152 
    153   // the position of the center of a muon ring
    154   const Float_t cenx = musearch.GetCenterX();
    155   const Float_t ceny = musearch.GetCenterY();
    156  
    157   for (Int_t i=0; i<entries; i++ )
    158     {
    159       MCerPhotPix &pix = (evt)[i];
    160      
    161       const MGeomPix &gpix = (geom)[pix.GetPixId()];
    162      
    163       const Float_t dx = gpix.GetX() - cenx;
    164       const Float_t dy = gpix.GetY() - ceny;
    165      
    166      const Float_t dist = TMath::Sqrt(dx*dx+dy*dy);
    167      
    168       Float_t ang = TMath::ACos(dx/dist);
    169       if(dy>0)
    170         ang *= -1.0;
    171      
    172       // if the signal is not near the estimated circle, it is ignored.
    173       if(dist < musearch.GetRadius() + fMargin && dist > musearch.GetRadius() - fMargin)
    174         {
    175           // check whether ummapped pixel is used or not.
    176           // if it is so, ingnore the pixel information since the pixels totally deteriorate the muon information.
    177           if(pix.IsPixelUnmapped())
    178             {
    179               fUseUnmap = kTRUE;
    180               continue;
    181             }
    182           fHistPhi->Fill(ang*kRad2Deg, pix.GetNumPhotons());
    183           fMuonSize += pix.GetNumPhotons();
    184         }
    185 
    186       if(pix.GetPixId()>397)
    187         continue;  // use only the inner pixles
    188      
    189       if(fUseCleanForWidth)
    190         {
    191           if(!pix.IsPixelUsed())
    192             continue;
    193         }
    194 
    195       fHistWidth->Fill(dist*geom.GetConvMm2Deg(), pix.GetNumPhotons());
    196     }
    197 
    198 
    199   // error estimation (temporaly)
    200   //  The error is estimated from the signal. In order to do so, we have to
    201   // once convert the signal from ADC to photo-electron. Then we can get
    202   // the fluctuation such as F-factor*sqrt(phe).
    203   //  Up to now, the error of pedestal is not taken into accout. This is not
    204   // of course correct. We will include this soon.
    205     Double_t ADC2PhEl = 0.14;
    206     Double_t Ffactor = 1.4;
    207     for(Int_t i=0; i<fArcPhiBinNum+1; i++)
    208       {
    209         Float_t rougherr  = TMath::Sqrt(TMath::Abs(fHistPhi->GetBinContent(i))*ADC2PhEl)/ADC2PhEl*Ffactor;
    210         {
    211           fHistPhi->SetBinError(i, rougherr);
    212         }
    213       }
    214     for(Int_t i=0; i<fArcWidthBinNum+1; i++)
    215       {
    216         Float_t rougherr = TMath::Sqrt(TMath::Abs(fHistWidth->GetBinContent(i))*ADC2PhEl)/ADC2PhEl*Ffactor;
    217         {
    218           fHistWidth->SetBinError(i, rougherr);
    219         }
    220       }
    221 }
    222 
    223 // --------------------------------------------------------------------------
    224 //
    225 //  Photon distribution along the estimated circle is fitted with theoritical
    226 // function in order to get some more information such as Arc Phi and Arc Length.
    227 //
    228 void MMuonCalibPar::CalcPhi
    229 (const MGeomCam &geom, const MCerPhotEvt &evt,
    230  const MMuonSearchPar &musearch)
    231 {
    232   Float_t convbin2val = (fArcPhiHistEndVal-fArcPhiHistStartVal)/
    233     (Float_t)fArcPhiBinNum;
    234 
    235     // adjust the peak to 0
    236     Float_t maxval = 0.;
    237     Int_t   maxbin = 0;
    238     maxval = fHistPhi->GetMaximum();
    239     maxbin = fHistPhi->GetMaximumBin();
    240     fPeakPhi = 180.-(Float_t)(maxbin-1)*convbin2val;
    241     TArrayD tmp;
    242     tmp.Set(fArcPhiBinNum+1);
    243     for(Int_t i=1; i<fArcPhiBinNum+1; i++)
    244       {
    245         tmp[i] = fHistPhi->GetBinContent(i);
    246       }
    247     for(Int_t i=1; i<fArcPhiBinNum+1; i++)
    248       {
    249         Int_t id;
    250         id = i + (maxbin-(Int_t)((Float_t)fArcPhiBinNum/2.)-1);
    251         if(id>fArcPhiBinNum)
    252           {
    253             id-=(fArcPhiBinNum);
    254           }
    255         if(id<=0)
    256           {
    257             id+=(fArcPhiBinNum);
    258           }
    259         fHistPhi->SetBinContent(i,tmp[id]);
    260       }
    261     maxbin = (Int_t)((Float_t)fArcPhiBinNum/2.)+1;
    262 
    263   // Determination of fitting region
    264   // The threshold is fixed with 100 [photons or ADC] in a bin. Therefore,
    265   // if you change the bin number, YOU HAVE TO CHANGE THIS VALUE!!!
    266   Float_t startfitval = 0.;
    267   Float_t endfitval   = 0.;
    268   Bool_t  IsInMaxim = kFALSE;
    269   Int_t   effbinnum = 0;
    270   for(Int_t i=1; i<fArcPhiBinNum+1; i++)
    271     {
    272       Float_t content = fHistPhi->GetBinContent(i);
    273       Float_t content_pre = fHistPhi->GetBinContent(i-1);
    274      
    275       if(content > fArcPhiThres && content_pre < fArcPhiThres)
    276         {
    277           startfitval = (Float_t)(i-1)*convbin2val+fArcPhiHistStartVal;
    278         }
    279       if(i==maxbin)
    280         IsInMaxim = kTRUE;
    281      
    282       if(content < fArcPhiThres && IsInMaxim == kTRUE)
    283         {
    284           endfitval = (Float_t)(i-1)*convbin2val+fArcPhiHistStartVal;
    285           break;
    286         }
    287       endfitval = fArcPhiHistEndVal;
    288     }
    289  
    290   effbinnum = (Int_t)((endfitval-startfitval)/convbin2val);
    291  
    292   fArcPhi = effbinnum*convbin2val;
    293   fArcLength = fArcPhi*3.1415926/180.*musearch.GetRadius()*geom.GetConvMm2Deg();
    294  
    295   if(fEnableImpactCalc)
    296       CalcImpact(geom, musearch, effbinnum, startfitval, endfitval);
    297 }
    298 
    299 // --------------------------------------------------------------------------
    300 //
    301 //  An impact parameter is calculated by fitting the histogram of photon
    302 // distribution along the circle with a theoritical model.
    303 // (See G. Vacanti et. al., Astroparticle Physics 2, 1994, 1-11.
    304 // The function (6) is used here.)
    305 //
    306 //  By default this calculation is suppressed because this calculation is
    307 // very time consuming. If you want to calculate an impact parameter,
    308 // you can call the function of EnableImpactCalc().
    309 //
    310 void MMuonCalibPar::CalcImpact
    311 ( const MGeomCam &geom, const MMuonSearchPar &musearch,
    312   Int_t effbinnum, Float_t startfitval, Float_t endfitval)
    313 {
    314   // Fit the distribution with Vacanti function. The function is different
    315   // for the impact parameter of inside or outside of our reflector.
    316   // Then two different functions are applied to the photon distribution,
    317   // and the one which give us smaller chisquare value is taken as a
    318   // proper one.
    319   Double_t val1,err1,val2,err2;
    320   // impact parameter inside mirror radius (8.5m)
    321   TString func1;
    322   Float_t tmpval = musearch.GetRadius()*geom.GetConvMm2Deg()*3.1415926/180.;
    323   tmpval = sin(2.*tmpval)*8.5;
    324   func1 += "[0]*";
    325   func1 += tmpval;
    326   func1 += "*(sqrt(1.-([1]/8.5)**2*sin((x-[2])*3.1415926/180.)**2)+([1]/8.5)*cos((x-[2])*3.1415926/180.))";
    327  
    328   TF1 f1("f1",func1,startfitval,endfitval);
    329   f1.SetParameters(2000,3,0);
    330   f1.SetParLimits(1,0,8.5);
    331   f1.SetParLimits(2,-180.,180.);
    332  
    333   fHistPhi->Fit("f1","RQEM");
    334  
    335   Float_t chi1 = -1;
    336   Float_t chi2 = -1;
    337   if(effbinnum>3)
    338     chi1 = f1.GetChisquare()/((Float_t)(effbinnum-3));
    339  
    340   gMinuit->GetParameter(1,val1,err1);  // get the estimated IP
    341   Float_t estip1 = val1;
    342  
    343   // impact parameter beyond mirror area (8.5m)
    344   TString func2;
    345   Float_t tmpval2 = musearch.GetRadius()*geom.GetConvMm2Deg()*3.1415926/180.;
    346   tmpval2 = sin(2.*tmpval2)*8.5*2.;
    347   func2 += "[0]*";
    348   func2 += tmpval2;
    349   func2 += "*sqrt(1.-(([1]/8.5)*sin((x-[2])*3.1415926/180.))**2)";
    350  
    351   TF1 f2("f2",func2,startfitval,endfitval);
    352   f2.SetParameters(2000,20,0);
    353   f2.SetParLimits(1,8.5,300.);
    354   f2.SetParLimits(2,-180.,180.);
    355  
    356   fHistPhi->Fit("f2","RQEM+");
    357  
    358   if(effbinnum>3)
    359     chi2 = f2.GetChisquare()/((Float_t)(effbinnum-3));
    360  
    361   gMinuit->GetParameter(1,val2,err2);  // get the estimated IP
    362   Float_t estip2 = val2;
    363  
    364   if(effbinnum<=3)
    365     {
    366       fEstImpact = -1.;
    367     }
    368   if(chi2 > chi1)
    369     {
    370       fEstImpact = estip1;
    371       fChiArcPhi = chi1;
    372     }
    373   else
    374     {
    375       fEstImpact = estip2;
    376       fChiArcPhi = chi2;
    377     }
    378 }
    379 
    380 // --------------------------------------------------------------------------
    381 //
    382 //  Photon distribution of distance from the center of estimated ring is
    383 // fitted in order to get some more information such as ARC WIDTH which
    384 // can represent to the PSF of our reflector.
    385 //
    386 Float_t MMuonCalibPar::CalcWidth
    387 (const MGeomCam &geom, const MCerPhotEvt &evt,
    388  const MMuonSearchPar &musearch)
    389 {
    390   Float_t convbin2val = (fArcWidthHistEndVal - fArcWidthHistStartVal)
    391     /(Float_t)fArcWidthBinNum;
    392 
    393     // determination of fitting region
    394     Int_t maxbin = fHistWidth->GetMaximumBin();
    395     Float_t startfitval = 0.;
    396     Float_t endfitval   = 0.;
    397     Bool_t   IsInMaxim = kFALSE;
    398     Int_t    effbinnum = 0;
    399     for(Int_t i=1; i<fArcWidthBinNum+1; i++)
    400       {
    401         Float_t content = fHistWidth->GetBinContent(i);
    402         Float_t content_pre = fHistWidth->GetBinContent(i-1);
    403 
    404         if(content > fArcWidthThres)
    405           effbinnum++;
    406 
    407         if(content > fArcWidthThres && content_pre < fArcWidthThres)
    408           {
    409             startfitval = (Float_t)(i-4)*convbin2val + fArcWidthHistStartVal;
    410             if(startfitval<0.) startfitval = 0.;
    411           }
    412         if(i==maxbin)
    413           IsInMaxim = kTRUE;
    414 
    415         if(content < fArcWidthThres && IsInMaxim == kTRUE)
    416           {
    417             endfitval = (Float_t)(i+2)*convbin2val + fArcWidthHistStartVal;
    418             if(endfitval>180.) endfitval = 180.;
    419             break;
    420           }
    421         endfitval = fArcWidthHistEndVal;
    422       }
    423     effbinnum = (Int_t)((endfitval-startfitval)/convbin2val);
    424 
    425     TF1 f1("f1","gaus",startfitval,endfitval);
    426 
    427     fHistWidth->Fit("f1","QR","",startfitval,endfitval);
    428    
    429     if(effbinnum>3)
    430       fChiArcWidth = f1.GetChisquare()/((Float_t)(effbinnum-3));
    431 
    432     Double_t val,err;
    433     gMinuit->GetParameter(2,val,err); // get the sigma value
    434 
    435     return val;
    436 }
    437 
    438 // --------------------------------------------------------------------------
    439 //
    440 //  Calculation of muon parameters
    441 //
    442 Int_t MMuonCalibPar::Calc
    443 (const MGeomCam &geom, const MCerPhotEvt &evt,
    444  MMuonSearchPar &musearch, const Float_t *cuts)
    445 {
    446   // sanity check
    447   if(evt.GetNumPixels() < 3)
    448     return kCONTINUE;
    449 
    450   // If an event does not seem to be like muon, the calculation will be skipped.
    451   if(musearch.IsNoMuon())
    452     return kCONTINUE;
    453 
    454   // Pre Cuts 1
    455   if(!fDisablePreCuts)
    456     {
    457       if(musearch.GetRadius() < cuts[0] || musearch.GetRadius() > cuts[1])
    458         {
    459           musearch.SetNoMuon();
    460           return kCONTINUE;
    461         }
    462       if(musearch.GetDeviation() > cuts[2])
    463         {
    464           musearch.SetNoMuon();
    465           return kCONTINUE;
    466         }
    467     }
    468 
    469   Reset();
    470 
    471   FillHist(geom,evt,musearch);
    472 
    473   CalcPhi(geom,evt,musearch);
    474 
    475   // Pre Cuts 2
    476   if(!fDisablePreCuts)
    477     {
    478       if(fMuonSize < cuts[3] || fArcPhi < cuts[4])
    479         {
    480           musearch.SetNoMuon();
    481           return kCONTINUE;
    482         }
    483     }
    484 
    485   fArcWidth = CalcWidth(geom,evt,musearch);
    486  
    487   SetReadyToSave();
    488 
    489   return kTRUE;
    490 }
    491123
    492124void MMuonCalibPar::Print(Option_t *) const

  • trunk/MagicSoft/Mars/mmuon/MMuonCalibParCalc.cc

    r5325 r5332  
    3030// Task to calculate the muon parameters
    3131//
    32 //  This class allows you to get more muon information especially usefule for
     32//  This class allows you to get more muon information especially useful for
    3333// the calibration of our telescope. This class store the information into the
    34 // container of MMuonCalibPar. The actual calculation is done in the class of
    35 // MMuonCalibPar. Please see the manual.
     34// container of MMuonCalibPar.
    3635//
    3736//  In order to make this class work, we need the information of the arc
    … …  
    5352// the followings;
    5453//
     54//   mucalibcalc.SetMargin(60.);
     55//
     56//  You can retrieve the histogram (TH1F) using the function of GetHistPhi()
     57// (also GetHistWid()). Therefore, you can draw the histogram such as
     58//
    5559//   MParList  plist;
    5660//   MMuonCalibPar muparcalib;
    5761//   plist.AddToList(&muparcalib);.
    58 //   muparcalib.SetMargin(60.);
    59 //
    60 //  You can retrieve the histogram (TH1F) using the function of GetHistPhi()
    61 // (also GetHistWid()). Therefore, you can draw the histogram such as
    62 //
    6362//   muparcalib.GetHistPhi().Draw();.
    6463//
    … …  
    7372// namely;
    7473//
    75 //   muparcalib.EnableImpactCalc();.
     74//   mucalibcalc.EnableImpactCalc();.
    7675//
    7776//  In addition, for the faster comutation, pre cuts to select the candidates
    … …  
    8685// function of DisablePreCuts(), namely;
    8786// 
    88 //   muparcalib.DisablePreCuts();.
     87//   mucalibcalc.DisablePreCuts();.
    8988//
    9089//
    … …  
    9796// done correctly.
    9897//
    99 //  Because of this rough treatment of the error, it may be better to use
    100 // the image after the cleaning for the evaluation of the Arc Width. There
    101 // is a possibility to use the image after the cleaning for the Arc Width.
    102 // You can use the function of UseCleanForWid(). However, this treatment
    103 // should be temporal and if the problem of the error is fixed, this function
    104 // will be removed.
    105 //
    10698//
    10799//  Input Containers:
    … …  
    119111#include <fstream>
    120112
     113#include <TH1.h>
     114#include <TF1.h>
     115#include <TMinuit.h>
     116
    121117#include "MParList.h"
    122118
    123119#include "MGeomCam.h"
     120#include "MGeomPix.h"
    124121#include "MSrcPosCam.h"
    125122#include "MCerPhotEvt.h"
    … …  
    128125#include "MLog.h"
    129126#include "MLogManip.h"
     127#include "MBinning.h"
    130128
    131129using namespace std;
    … …  
    141139//
    142140MMuonCalibParCalc::MMuonCalibParCalc(const char *name, const char *title)
    143   : fCerPhotName("MCerPhotEvt")
     141  : fNameCerPhot("MCerPhotEvt")
    144142{
    145143    fName  = name  ? name  : gsDefName.Data();
    … …  
    151149    fPreCuts[3]     = 2000.;
    152150    fPreCuts[4]     = 150.;
     151
     152    fMargin = 60.;
     153    fArcPhiThres = 100.;
     154    fArcWidthThres = 100.;
     155
     156    fEnableImpactCalc = kFALSE; // By default the calculation of impact parameter is skipped.
     157    fDisablePreCuts =   kFALSE; // By default the pre cuts will be applied.
    153158}
    154159
    … …  
    157162Int_t MMuonCalibParCalc::PreProcess(MParList *pList)
    158163{
    159     fCerPhotEvt = (MCerPhotEvt*)pList->FindObject(AddSerialNumber(fCerPhotName), "MCerPhotEvt");
    160     if (!fCerPhotEvt)
    161     {
    162         *fLog << err << fCerPhotName << " [MCerPhotEvt] not found... aborting." << endl;
    163         return kFALSE;
    164     }
    165 
    166     fGeomCam = (MGeomCam*)pList->FindObject("MGeomCam");
    167     if (!fGeomCam)
    168     {
    169         *fLog << err << "MGeomCam not found... aborting." << 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;
    182 }
     164  fCerPhotEvt = (MCerPhotEvt*)pList->FindObject(AddSerialNumber(fNameCerPhot), "MCerPhotEvt");
     165  if (!fCerPhotEvt)
     166    {
     167      *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl;
     168      return kFALSE;
     169    }
     170 
     171  fGeomCam = (MGeomCam*)pList->FindObject("MGeomCam");
     172  if (!fGeomCam)
     173    {
     174      *fLog << dbginf << "MGeomCam (Camera Geometry) missing in Parameter List... aborting." << endl;
     175      return kFALSE;
     176    }
     177 
     178  fMuonCalibPar = (MMuonCalibPar*)pList->FindCreateObj("MMuonCalibPar", "MMuonCalibPar");
     179  if (!fMuonCalibPar)
     180    {
     181      *fLog << dbginf << "MMuonCalibPar missing in Parameter List... aborting." << endl;
     182      return kFALSE;
     183    }
     184 
     185  fMuonSearchPar = (MMuonSearchPar*)pList->FindCreateObj("MMuonSearchPar", "MMuonSearchPar");
     186  if (!fMuonSearchPar)
     187    {
     188      *fLog << dbginf << "MMuonSearchPar missing in Parameter List... aborting." << endl;
     189      return kFALSE;
     190    }
     191 
     192  return kTRUE;
     193}
     194
     195// --------------------------------------------------------------------------
     196//
     197//  This function fill the histograms in order to get muon parameters.
     198// For the evaluation of the Arc Width, we use only the signals in the inner
     199// part.
     200//
     201void MMuonCalibParCalc::FillHist()
     202{
     203  Float_t MuonSize = 0.;
     204
     205  Int_t binnumphi = fMuonCalibPar->fArcPhiBinNum;
     206  Int_t binnumwid = fMuonCalibPar->fArcWidthBinNum;
     207
     208  // preparation for a histgram
     209  MBinning binsphi;
     210  binsphi.SetEdges(binnumphi,
     211                   fMuonCalibPar->fArcPhiHistStartVal,
     212                   fMuonCalibPar->fArcPhiHistEndVal);
     213  binsphi.Apply(*(fMuonCalibPar->fHistPhi));
     214
     215  MBinning binswid;
     216  binswid.SetEdges(binnumwid,
     217                   fMuonCalibPar->fArcWidthHistStartVal,
     218                   fMuonCalibPar->fArcWidthHistEndVal);
     219  binswid.Apply(*(fMuonCalibPar->fHistWidth));
     220
     221  const Int_t entries = (*fCerPhotEvt).GetNumPixels();
     222
     223  // the position of the center of a muon ring
     224  const Float_t cenx = (*fMuonSearchPar).GetCenterX();
     225  const Float_t ceny = (*fMuonSearchPar).GetCenterY();
     226 
     227  for (Int_t i=0; i<entries; i++ )
     228    {
     229      MCerPhotPix &pix = (*fCerPhotEvt)[i];
     230     
     231      const MGeomPix &gpix = (*fGeomCam)[pix.GetPixId()];
     232     
     233      const Float_t dx = gpix.GetX() - cenx;
     234      const Float_t dy = gpix.GetY() - ceny;
     235     
     236     const Float_t dist = TMath::Sqrt(dx*dx+dy*dy);
     237     
     238      Float_t ang = TMath::ACos(dx/dist);
     239      if(dy>0)
     240        ang *= -1.0;
     241     
     242      // if the signal is not near the estimated circle, it is ignored.
     243      if(dist < (*fMuonSearchPar).GetRadius() + fMuonCalibPar->GetMargin()
     244         && dist > (*fMuonSearchPar).GetRadius() - fMuonCalibPar->GetMargin())
     245        {
     246          // check whether ummapped pixel is used or not.
     247          // if it is so, ingnore the pixel information since the pixels totally deteriorate the muon information.
     248          if(pix.IsPixelUnmapped())
     249            {
     250              fMuonCalibPar->SetUseUnmap();
     251              continue;
     252            }
     253          fMuonCalibPar->fHistPhi->Fill(ang*kRad2Deg, pix.GetNumPhotons());
     254          MuonSize += pix.GetNumPhotons();
     255        }
     256
     257      if(pix.GetPixId()>397)
     258        continue;  // use only the inner pixles
     259     
     260      fMuonCalibPar->fHistWidth
     261        ->Fill(dist*(*fGeomCam).GetConvMm2Deg(), pix.GetNumPhotons());
     262    }
     263
     264  fMuonCalibPar->SetMuonSize(MuonSize);
     265
     266  // error estimation (temporaly)
     267  //  The error is estimated from the signal. In order to do so, we have to
     268  // once convert the signal from ADC to photo-electron. Then we can get
     269  // the fluctuation such as F-factor*sqrt(phe).
     270  //  Up to now, the error of pedestal is not taken into accout. This is not
     271  // of course correct. We will include this soon.
     272    Double_t ADC2PhEl = 0.14;
     273    Double_t Ffactor = 1.4;
     274    for(Int_t i=0; i<binnumphi+1; i++)
     275      {
     276        Float_t rougherr  = TMath::Sqrt(TMath::Abs(fMuonCalibPar->fHistPhi->GetBinContent(i))*ADC2PhEl)/ADC2PhEl*Ffactor;
     277        {
     278          fMuonCalibPar->fHistPhi->SetBinError(i, rougherr);
     279        }
     280      }
     281    for(Int_t i=0; i<binnumwid+1; i++)
     282      {
     283        Float_t rougherr = TMath::Sqrt(TMath::Abs(fMuonCalibPar->fHistWidth->GetBinContent(i))*ADC2PhEl)/ADC2PhEl*Ffactor;
     284        {
     285          fMuonCalibPar->fHistWidth->SetBinError(i, rougherr);
     286        }
     287      }
     288}
     289
     290// --------------------------------------------------------------------------
     291//
     292//  Photon distribution along the estimated circle is fitted with theoritical
     293// function in order to get some more information such as Arc Phi and Arc
     294// Length.
     295//
     296void MMuonCalibParCalc::CalcPhi()
     297{
     298  Float_t thres = fMuonCalibPar->GetArcPhiThres();
     299  Float_t startval = fMuonCalibPar->fArcPhiHistStartVal;
     300  Float_t endval = fMuonCalibPar->fArcPhiHistEndVal;
     301  Int_t   binnum = fMuonCalibPar->fArcPhiBinNum;
     302
     303  Float_t convbin2val = (endval-startval)/(Float_t)binnum;
     304
     305    // adjust the peak to 0
     306    Float_t maxval = 0.;
     307    Int_t   maxbin = 0;
     308    maxval = fMuonCalibPar->fHistPhi->GetMaximum();
     309    maxbin = fMuonCalibPar->fHistPhi->GetMaximumBin();
     310    fMuonCalibPar->SetPeakPhi(180.-(Float_t)(maxbin-1.)*convbin2val);
     311    TArrayD tmp;
     312    tmp.Set(binnum+1);
     313    for(Int_t i=1; i<binnum+1; i++)
     314      {
     315        tmp[i] = fMuonCalibPar->fHistPhi->GetBinContent(i);
     316      }
     317    for(Int_t i=1; i<binnum+1; i++)
     318      {
     319        Int_t id;
     320        id = i + (maxbin-(Int_t)((Float_t)binnum/2.)-1);
     321        if(id>binnum)
     322          {
     323            id-=(binnum);
     324          }
     325        if(id<=0)
     326          {
     327            id+=(binnum);
     328          }
     329        fMuonCalibPar->fHistPhi->SetBinContent(i,tmp[id]);
     330      }
     331    maxbin = (Int_t)((Float_t)binnum/2.)+1;
     332
     333  // Determination of fitting region
     334  // The threshold is fixed with 100 [photons or ADC] in a bin. Therefore,
     335  // if you change the bin number, YOU HAVE TO CHANGE THIS VALUE!!!
     336  Float_t startfitval = 0.;
     337  Float_t endfitval   = 0.;
     338  Bool_t  IsInMaxim = kFALSE;
     339  Int_t effbinnum = 0;
     340  for(Int_t i=1; i<binnum+1; i++)
     341    {
     342      Float_t content = fMuonCalibPar->fHistPhi->GetBinContent(i);
     343      Float_t content_pre = fMuonCalibPar->fHistPhi->GetBinContent(i-1);
     344     
     345      if(content > thres && content_pre < thres)
     346        {
     347          startfitval = (Float_t)(i-1)*convbin2val+startval;
     348        }
     349      if(i==maxbin)
     350        IsInMaxim = kTRUE;
     351     
     352      if(content < thres && IsInMaxim == kTRUE)
     353        {
     354          endfitval = (Float_t)(i-1)*convbin2val+startval;
     355          break;
     356        }
     357      endfitval = endval;
     358    }
     359 
     360  effbinnum = (Int_t)((endfitval-startfitval)/convbin2val);
     361
     362  fMuonCalibPar->SetArcPhi(endfitval-startfitval);
     363
     364  fMuonCalibPar->SetArcLength( fMuonCalibPar->GetArcPhi()*TMath::DegToRad()*(*fMuonSearchPar).GetRadius()*(*fGeomCam).GetConvMm2Deg());
     365 
     366  if(fEnableImpactCalc)
     367      CalcImpact(effbinnum, startfitval, endfitval);
     368}
     369
     370// --------------------------------------------------------------------------
     371//
     372//  An impact parameter is calculated by fitting the histogram of photon
     373// distribution along the circle with a theoritical model.
     374// (See G. Vacanti et. al., Astroparticle Physics 2, 1994, 1-11.
     375// The function (6) is used here.)
     376//
     377//  By default this calculation is suppressed because this calculation is
     378// very time consuming. If you want to calculate an impact parameter,
     379// you can call the function of EnableImpactCalc().
     380//
     381void MMuonCalibParCalc::CalcImpact
     382(Int_t effbinnum, Float_t startfitval, Float_t endfitval)
     383{
     384  // Fit the distribution with Vacanti function. The function is different
     385  // for the impact parameter of inside or outside of our reflector.
     386  // Then two different functions are applied to the photon distribution,
     387  // and the one which give us smaller chisquare value is taken as a
     388  // proper one.
     389  Double_t val1,err1,val2,err2;
     390  // impact parameter inside mirror radius (8.5m)
     391  TString func1;
     392  Float_t tmpval = (*fMuonSearchPar).GetRadius()*(*fGeomCam).GetConvMm2Deg()*TMath::DegToRad();
     393  tmpval = sin(2.*tmpval)*8.5;
     394  func1 += "[0]*";
     395  func1 += tmpval;
     396  func1 += "*(sqrt(1.-([1]/8.5)**2*sin((x-[2])*3.1415926/180.)**2)+([1]/8.5)*cos((x-[2])*3.1415926/180.))";
     397 
     398  TF1 f1("f1",func1,startfitval,endfitval);
     399  f1.SetParameters(2000,3,0);
     400  f1.SetParLimits(1,0,8.5);
     401  f1.SetParLimits(2,-180.,180.);
     402 
     403  fMuonCalibPar->fHistPhi->Fit("f1","RQEM");
     404 
     405  Float_t chi1 = -1;
     406  Float_t chi2 = -1;
     407  if(effbinnum>3)
     408    chi1 = f1.GetChisquare()/((Float_t)(effbinnum-3));
     409 
     410  gMinuit->GetParameter(1,val1,err1);  // get the estimated IP
     411  Float_t estip1 = val1;
     412 
     413  // impact parameter beyond mirror area (8.5m)
     414  TString func2;
     415  Float_t tmpval2 = (*fMuonSearchPar).GetRadius()*(*fGeomCam).GetConvMm2Deg()*TMath::DegToRad();
     416  tmpval2 = sin(2.*tmpval2)*8.5*2.;
     417  func2 += "[0]*";
     418  func2 += tmpval2;
     419  func2 += "*sqrt(1.-(([1]/8.5)*sin((x-[2])*3.1415926/180.))**2)";
     420 
     421  TF1 f2("f2",func2,startfitval,endfitval);
     422  f2.SetParameters(2000,20,0);
     423  f2.SetParLimits(1,8.5,300.);
     424  f2.SetParLimits(2,-180.,180.);
     425 
     426  fMuonCalibPar->fHistPhi->Fit("f2","RQEM+");
     427 
     428  if(effbinnum>3)
     429    chi2 = f2.GetChisquare()/((Float_t)(effbinnum-3));
     430 
     431  gMinuit->GetParameter(1,val2,err2);  // get the estimated IP
     432  Float_t estip2 = val2;
     433 
     434  if(effbinnum<=3)
     435    {
     436      fMuonCalibPar->SetEstImpact(-1.);
     437    }
     438  if(chi2 > chi1)
     439    {
     440      fMuonCalibPar->SetEstImpact(estip1);
     441      fMuonCalibPar->SetChiArcPhi(chi1);
     442    }
     443  else
     444    {
     445      fMuonCalibPar->SetEstImpact(estip2);
     446      fMuonCalibPar->SetChiArcPhi(chi2);
     447    }
     448}
     449
     450// --------------------------------------------------------------------------
     451//
     452//  Photon distribution of distance from the center of estimated ring is
     453// fitted in order to get some more information such as ARC WIDTH which
     454// can represent to the PSF of our reflector.
     455//
     456Float_t MMuonCalibParCalc::CalcWidth()
     457{
     458  Float_t startval = fMuonCalibPar->fArcWidthHistStartVal;
     459  Float_t endval = fMuonCalibPar->fArcWidthHistEndVal;
     460  Int_t   binnum = fMuonCalibPar->fArcWidthBinNum;
     461  Float_t thres = fMuonCalibPar->GetArcWidthThres();
     462
     463  Float_t convbin2val = (endval - startval)
     464    /(Float_t)binnum;
     465
     466    // determination of fitting region
     467    Int_t maxbin = fMuonCalibPar->fHistWidth->GetMaximumBin();
     468    Float_t startfitval = 0.;
     469    Float_t endfitval   = 0.;
     470    Bool_t   IsInMaxim = kFALSE;
     471    Int_t    effbinnum = 0;
     472    for(Int_t i=1; i<binnum+1; i++)
     473      {
     474        Float_t content = fMuonCalibPar->fHistWidth->GetBinContent(i);
     475        Float_t content_pre = fMuonCalibPar->fHistWidth->GetBinContent(i-1);
     476
     477        if(content > thres)
     478          effbinnum++;
     479
     480        if(content > thres && content_pre < thres)
     481          {
     482            startfitval = (Float_t)(i-4)*convbin2val + startval;
     483            if(startfitval<0.) startfitval = 0.;
     484          }
     485        if(i==maxbin)
     486          IsInMaxim = kTRUE;
     487
     488        if(content < thres && IsInMaxim == kTRUE)
     489          {
     490            endfitval = (Float_t)(i+2)*convbin2val + startval;
     491            if(endfitval>180.) endfitval = 180.;
     492            break;
     493          }
     494        endfitval = endval;
     495      }
     496    effbinnum = (Int_t)((endfitval-startfitval)/convbin2val);
     497
     498    TF1 f1("f1","gaus",startfitval,endfitval);
     499
     500    fMuonCalibPar->fHistWidth->Fit("f1","QR","",startfitval,endfitval);
     501   
     502    if(effbinnum>3)
     503      fMuonCalibPar->SetChiArcWidth(f1.GetChisquare()/((Float_t)(effbinnum-3)));
     504
     505    Double_t val,err;
     506    gMinuit->GetParameter(2,val,err); // get the sigma value
     507
     508    return val;
     509}
     510
     511// --------------------------------------------------------------------------
     512//
     513//  Calculation of muon parameters
     514//
     515Int_t MMuonCalibParCalc::Calc(const Float_t *cuts)
     516{
     517  // sanity check
     518  if((*fCerPhotEvt).GetNumPixels() < 3)
     519    return kCONTINUE;
     520
     521  // If an event does not seem to be like muon, the calculation will be skipped.
     522  if((*fMuonSearchPar).IsNoMuon())
     523    return kCONTINUE;
     524
     525  // Pre Cuts 1
     526  if(!fDisablePreCuts)
     527    {
     528      if((*fMuonSearchPar).GetRadius() < cuts[0] || (*fMuonSearchPar).GetRadius() > cuts[1])
     529        {
     530          (*fMuonSearchPar).SetNoMuon();
     531          return kCONTINUE;
     532        }
     533      if((*fMuonSearchPar).GetDeviation() > cuts[2])
     534        {
     535          (*fMuonSearchPar).SetNoMuon();
     536          return kCONTINUE;
     537        }
     538    }
     539
     540  // initialization
     541  (*fMuonCalibPar).Reset();
     542
     543  // Fill signals to histograms
     544  FillHist();
     545
     546  // Calculation of Arc Phi etc...
     547  CalcPhi();
     548
     549  // Pre Cuts 2
     550  if(!fDisablePreCuts)
     551    {
     552      if(fMuonCalibPar->GetMuonSize() < cuts[3]
     553         || fMuonCalibPar->GetArcPhi() < cuts[4])
     554        {
     555          (*fMuonSearchPar).SetNoMuon();
     556          return kCONTINUE;
     557        }
     558    }
     559
     560  // Calculation of Arc Width etc...
     561  fMuonCalibPar->SetArcWidth(CalcWidth());
     562 
     563  return kTRUE;
     564}
     565
    183566
    184567// -------------------------------------------------------------------------
    … …  
    186569Int_t MMuonCalibParCalc::Process()
    187570{
    188 
    189   if(!fMuonCalibPar->Calc(*fGeomCam, *fCerPhotEvt, *fMuonSearchPar, fPreCuts))
     571  fMuonCalibPar->SetMargin(fMargin);
     572  fMuonCalibPar->SetArcPhiThres(fArcPhiThres);
     573  fMuonCalibPar->SetArcWidthThres(fArcWidthThres);
     574
     575  if(!Calc(fPreCuts))
    190576    return kCONTINUE;
    191577
    192578  return kTRUE;
    193 }
    194 
    195 void MMuonCalibParCalc::SetMargin(Float_t margin)
    196 {
    197   fMuonCalibPar->SetMargin(margin);
    198 }
    199 
    200 void MMuonCalibParCalc::EnableImpactCalc()
    201 {
    202   fMuonCalibPar->EnableImpactCalc();
    203 }
    204 
    205 void MMuonCalibParCalc::DisablePreCuts()
    206 {
    207   fMuonCalibPar->DisablePreCuts();
    208579}
    209580

  • trunk/MagicSoft/Mars/mmuon/MMuonCalibParCalc.h

    r5210 r5332  
    2020    MMuonSearchPar *fMuonSearchPar;
    2121
     22    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
     23    Float_t fArcPhiThres;     // The threshold value to define arc phi
     24    Float_t fArcWidthThres;   // The threshold value to define arc width
     25    Bool_t fEnableImpactCalc; // If true, the impact calculation will be done, which consumes a lot of time.
     26    Bool_t fDisablePreCuts;   // If true, the pre cuts to select muons for the calibration will be disabled.
     27
    2228    Float_t fPreCuts[5];  // The values for pre cuts.
    2329
    … …  
    2531    Int_t Process();
    2632
    27     TString fCerPhotName;
     33    TString fNameCerPhot;
    2834
    2935public:
    3036    MMuonCalibParCalc(const char *name=NULL, const char *title=NULL);
    3137
    32     void SetMargin(Float_t margin);
    33     void EnableImpactCalc();
    34     void DisablePreCuts();
     38    void SetMargin(Float_t margin)       { fMargin = margin; }
     39    void SetArcPhiThres(Float_t thres)   { fArcPhiThres = thres; }
     40    void SetArcWidthThres(Float_t thres) { fArcWidthThres = thres; }
     41    void EnableImpactCalc()              { fEnableImpactCalc = kTRUE; }
     42    void DisablePreCuts()                { fDisablePreCuts = kTRUE; }
    3543    void SetPreCuts(Float_t radcutlow, Float_t radcuthigh, Float_t devcuthigh,
    3644                    Float_t musizecutlow, Float_t arcphicutlow);
    3745
    38     void SetNameCerPhotEvt(const char *name) { fCerPhotName = name; }
     46    void SetNameCerPhotEvt(const char *name) { fNameCerPhot = name; }
     47
     48    void FillHist();
     49    void CalcPhi();
     50    void CalcImpact(Int_t effbinnum, Float_t startfitval, Float_t endfitval);
     51    Float_t CalcWidth();
     52    Int_t   Calc(const Float_t *cuts);
    3953
    4054    ClassDef(MMuonCalibParCalc, 0) // task to calculate muon parameters

  • trunk/MagicSoft/Mars/mmuon/Makefile

    r5172 r5332  
    2626
    2727SRCFILES = MMuonSearchPar.cc \
    28            MMuonCalibPar.cc
     28           MMuonSearchParCalc.cc \
     29           MMuonCalibPar.cc \
     30           MMuonCalibParCalc.cc
    2931
    3032############################################################
    … …  
    3739
    3840mrproper:       clean rmbak
     41
     42
     43
     44
     45
     46
     47

  • trunk/MagicSoft/Mars/mmuon/MuonLinkDef.h

    r5173 r5332  
    66
    77#pragma link C++ class MMuonSearchPar+;
     8#pragma link C++ class MMuonSearchParCalc+;
    89#pragma link C++ class MMuonCalibPar+;
     10#pragma link C++ class MMuonCalibParCalc+;
    911
    1012#endif
Note: See TracChangeset for help on using the changeset viewer.