Changeset 1748


Ignore:
Timestamp:
02/07/03 13:26:50 (22 years ago)
Author:
wittek
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
4 added
4 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Mars/manalysis/MSigmabar.cc

    r1682 r1748  
    1616!
    1717!
    18 !   Author(s): Robert Wagner  10/2002 <mailto:magicsoft@rwagner.de>
     18!   Author(s): Robert Wagner   10/2002 <mailto:magicsoft@rwagner.de>
     19!              Wolfgang Wittek 01/2003 <mailto:wittek@mppmu.mpg.de>
    1920!
    2021!   Copyright: MAGIC Software Development, 2002
     
    2728// MSigmabar                                                               //
    2829//                                                                         //
    29 // This is the storage container to hold information about the mean sigma  //
    30 // (aka Sigmabar) of all pedestals                                         //
     30// This is the storage container to hold information about                 //
     31// "average" pedestal sigmas                                               //
     32//                                                                         //
     33// In calculating averages all sigmas are normalized to the area of pixel 0//
    3134//                                                                         //
    3235/////////////////////////////////////////////////////////////////////////////
     
    3740#include "MLog.h"
    3841#include "MLogManip.h"
    39 
     42#include "MParList.h"
    4043#include "MGeomCam.h"
    4144#include "MPedestalCam.h"
     45#include "MPedestalPix.h"
    4246#include "MGeomPix.h"
    43 #include "MPedestalPix.h"
     47#include "MCerPhotEvt.h"
     48#include "MCerPhotPix.h"
    4449
    4550ClassImp(MSigmabar);
    4651
     52// --------------------------------------------------------------------------
     53//
    4754MSigmabar::MSigmabar(const char *name, const char *title)
    4855{
     
    5057    fTitle = title ? title : "Storage container for Sigmabar";
    5158   
    52     fSigmabar = 0.0;
    53     fSigmabarInner = 0.0;
    54     fSigmabarOuter = 0.0;
    55     fRatioA = 0.0;
    5659
    5760    fCalcPixNum=kTRUE;
    5861}
    5962
     63// --------------------------------------------------------------------------
     64//
    6065MSigmabar::~MSigmabar()
    6166{
    6267  // do nothing special.
    6368}
    64 
    6569
    6670// --------------------------------------------------------------------------
     
    6872// Actual calculation of sigmabar. This is done for each of the six sectors
    6973// separately due to their possibly different HV behavior. Also inner and
    70 // outer pixels are treated separately
     74// outer pixels are treated separately.
    7175//
    7276// Preliminary! Works for CT1 test, for real MAGIC crosschecks still have
     
    7478// determination of sector to which a respective pixel belongs
    7579//
    76 Bool_t MSigmabar::Calc(const MGeomCam &geom, const MPedestalCam &ped)
    77 {
    78   const UInt_t npix = ped.GetSize();
     80Float_t MSigmabar::Calc(const MGeomCam &geom, const MPedestalCam &ped,
     81                        const MCerPhotEvt &evt)
     82{
     83  fSigmabar      = 0.0;
     84  fSigmabarInner = 0.0;
     85  fSigmabarOuter = 0.0;
     86
     87
    7988  Float_t innerSquaredSum[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
    8089  Float_t outerSquaredSum[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
    81   Int_t innerPixels[6] = {0,0,0,0,0,0};
    82   Int_t outerPixels[6] = {0,0,0,0,0,0};
     90  Int_t innerPixels[6]       = {0,0,0,0,0,0};
     91  Int_t outerPixels[6]       = {0,0,0,0,0,0};
     92
    8393  Int_t currentSector;
    8494  Float_t angle;
    8595 
    86   for (UInt_t i=0; i<npix;i++)
    87     {
    88       const MGeomPix gpix = geom[i];
    89       angle=((6*atan2(gpix.GetX(),gpix.GetY())/(2*TMath::Pi()))-0.5);
    90       if (angle<0) angle+=6;
    91       currentSector=(Int_t)angle;         
    92       geom.GetPixRatio(i) == 1 ? innerPixels[currentSector]++ : outerPixels[currentSector]++;   
    93      
    94       const MPedestalPix pix = ped[i];
    95       geom.GetPixRatio(i) == 1 ? innerSquaredSum[currentSector]+=(pix.GetSigma()*pix.GetSigma()) : outerSquaredSum[currentSector]+=((pix.GetSigma()*pix.GetSigma()) / geom.GetPixRatio(i));
    96 
    97       // Get once and forever the ratio of areas outer/inner pixels
    98       if (fRatioA && (geom.GetPixRatio(i) != 1)) fRatioA=geom.GetPixRatio(i);
    99     }
    100 
    101   // Overall Sigma
     96  // sum up sigma**2 for each sector, separately for inner and outer region;
     97  // all pixels are renormalized to the area of pixel 0
     98  //
     99  // consider all pixels with Cherenkov photon information
     100  // and require "Used"
     101  //
     102
     103  const UInt_t npix = evt.GetNumPixels();
     104
     105  for (UInt_t i=0; i<npix; i++)
     106  {
     107      MCerPhotPix &cerpix = evt.operator[](i);
     108      if (!cerpix.IsPixelUsed())
     109        continue;
     110 
     111      if ( cerpix.GetNumPhotons() == 0 )
     112      {
     113        *fLog << "MSigmabar::Calc(); no.of photons is 0 for used pixel"
     114              << endl;
     115        continue;
     116      }
     117
     118      Int_t j = cerpix.GetPixId();
     119      Double_t Area = geom.GetPixRatio(j);
     120
     121      const MGeomPix    &gpix = geom[j];
     122      const MPedestalPix &pix =  ped[j];
     123
     124      //angle = 6.0*atan2(gpix.GetX(),gpix.GetY()) / (2.0*TMath::Pi()) - 0.5;
     125      //if (angle<0.0) angle+=6.0;
     126
     127      angle = 6.0*atan2(gpix.GetY(),gpix.GetX()) / (2.0*TMath::Pi());
     128      if (angle<0.0) angle+=6.0;
     129      currentSector=(Int_t)angle;
     130       
     131      // count only those pixels which have a sigma != 0.0
     132      Float_t sigma = pix.GetMeanRms();
     133
     134      if ( sigma != 0.0 )
     135      { 
     136        if (Area < 1.5)
     137        {
     138          innerPixels[currentSector]++;
     139          innerSquaredSum[currentSector]+= sigma*sigma / Area;
     140        }
     141        else
     142        {
     143          outerPixels[currentSector]++;
     144          outerSquaredSum[currentSector]+= sigma*sigma / Area;
     145        }
     146      }
     147  }
     148 
    102149  fSigmabarInner=0; fInnerPixels=0;
    103150  fSigmabarOuter=0; fOuterPixels=0;
    104151  for (UInt_t i=0; i<6; i++) {
    105     fSigmabarInner+=innerSquaredSum[i];
    106     fInnerPixels  +=innerPixels[i];
    107     fSigmabarOuter+=outerSquaredSum[i];
    108     fOuterPixels  +=outerPixels[i];
     152    fSigmabarInner += innerSquaredSum[i];
     153    fInnerPixels   += innerPixels[i];
     154    fSigmabarOuter += outerSquaredSum[i];
     155    fOuterPixels   += outerPixels[i];
    109156  }
    110157
    111   fSigmabarInner/=fInnerPixels;
    112   if (fSigmabarOuter != 0) fSigmabarOuter/=fOuterPixels;
    113   fSigmabar=sqrt(fSigmabarInner + fSigmabarOuter/( fOuterPixels==0 ? 1 : fRatioA));
    114   fSigmabarInner=sqrt(fSigmabarInner);
    115   fSigmabarOuter=sqrt(fSigmabarOuter);
     158
     159  // this is the sqrt of the average sigma**2;
     160  //
     161  if ( (fInnerPixels+fOuterPixels) > 0)
     162    fSigmabar=sqrt(   (fSigmabarInner + fSigmabarOuter)
     163                    / (fInnerPixels   + fOuterPixels) );
     164
     165  if (fInnerPixels > 0) fSigmabarInner /= fInnerPixels;
     166  if (fOuterPixels > 0) fSigmabarOuter /= fOuterPixels;
     167 
     168  // this is the sqrt of the average sigma**2
     169  // for the inner and outer pixels respectively
     170  //
     171  fSigmabarInner = sqrt( fSigmabarInner );
     172  fSigmabarOuter = sqrt( fSigmabarOuter );
    116173 
     174
    117175  for (UInt_t i=0; i<6; i++) {
    118     fSigmabarInnerSector[i]=innerSquaredSum[i]/innerPixels[i];
    119     fSigmabarOuterSector[i]=outerSquaredSum[i]/( outerSquaredSum[i]==0 ? 1: outerPixels[i] );
    120 
    121     fSigmabarSector[i]=sqrt(fSigmabarInnerSector[i] + fSigmabarOuterSector[i]*fRatioA);
    122     fSigmabarInnerSector[i]=sqrt(fSigmabarInnerSector[i]);
    123     fSigmabarOuterSector[i]=sqrt(fSigmabarOuterSector[i]);
     176    fSigmabarSector[i]         = 0.0;
     177    fSigmabarInnerSector[i]    = 0.0;
     178    fSigmabarOuterSector[i]    = 0.0;
     179
     180    if ( (innerPixels[i]+outerPixels[i]) > 0)
     181      fSigmabarSector[i] = sqrt(  (innerSquaredSum[i] + outerSquaredSum[i])
     182                                / (innerPixels[i]     + outerPixels[i]    ) );
     183    if ( innerPixels[i] > 0)
     184      fSigmabarInnerSector[i] = innerSquaredSum[i] / innerPixels[i];
     185    if ( outerPixels[i] > 0)
     186      fSigmabarOuterSector[i] = outerSquaredSum[i] / outerPixels[i];
     187
     188
     189    fSigmabarInnerSector[i] = sqrt( fSigmabarInnerSector[i] );
     190    fSigmabarOuterSector[i] = sqrt( fSigmabarOuterSector[i] );
    124191  }
    125192   
    126   // Did all our calculations work? fOuterPixels==0 could happen, however.
    127   return (fInnerPixels!=0);
    128 }
    129 
    130 
     193  return (fSigmabar);
     194}
     195
     196// --------------------------------------------------------------------------
     197//
    131198void MSigmabar::Print(Option_t *) const
    132199{
    133   *fLog << all;
    134   *fLog << "Sigmabar     Overall " << fSigmabar;
     200  *fLog << endl;
     201  *fLog << "Total number of inner pixels is " << fInnerPixels << endl;
     202  *fLog << "Total number of outer pixels is " << fOuterPixels << endl;
     203  *fLog << endl;
     204
     205  *fLog << "Sigmabar     Overall : " << fSigmabar << "   ";
    135206  *fLog << " Sectors: ";
    136   for (Int_t i=0;i<6;i++) *fLog << fSigmabarSector[i] << " ";
    137   *fLog << endl;
    138   *fLog << "Sigmabar     Inner   " << fSigmabarInner;
     207  for (Int_t i=0;i<6;i++) *fLog << fSigmabarSector[i] << ", ";
     208  *fLog << endl;
     209
     210
     211  *fLog << "Sigmabar     Inner   : " << fSigmabarInner << "   ";
    139212  *fLog << " Sectors: ";
    140   for (Int_t i=0;i<6;i++) *fLog << fSigmabarInnerSector[i] << " ";
    141   *fLog << endl;
    142   *fLog << "Sigmabar     Outer   " << fSigmabarOuter;
     213  for (Int_t i=0;i<6;i++) *fLog << fSigmabarInnerSector[i] << ", ";
     214  *fLog << endl;
     215
     216
     217  *fLog << "Sigmabar     Outer   : " << fSigmabarOuter << "   ";
    143218  *fLog << " Sectors: ";
    144   for (Int_t i=0;i<6;i++) *fLog << fSigmabarOuterSector[i] << " ";
    145   *fLog << endl;
    146   *fLog << "Total number of inner pixels found is " << fInnerPixels << endl;
    147   *fLog << "Total number of outer pixels found is " << fOuterPixels << endl;
    148   *fLog << "Ratio of areas outer/inner pixels found is " << fRatioA << endl;
    149   *fLog << endl;
    150 }
     219  for (Int_t i=0;i<6;i++) *fLog << fSigmabarOuterSector[i] << ", ";
     220  *fLog << endl;
     221
     222}
     223// --------------------------------------------------------------------------
     224
     225
     226
     227
  • trunk/MagicSoft/Mars/manalysis/MSigmabar.h

    r1693 r1748  
    44#ifndef MARS_MParContainer
    55#include "MParContainer.h"
     6#endif
     7
     8#ifndef MARS_MParList
     9#include "MParList.h"
    610#endif
    711
     
    1418#endif
    1519
     20class MCerPhotEvt;
     21
    1622class MSigmabar : public MParContainer
    1723{
    1824private:
    19     Float_t fSigmabar; // Sigmabar (mean standard deviation) of pedestal
     25    Float_t fSigmabar;          // Sigmabar ("average" RMS) of pedestal
     26    Float_t fSigmabarInner;     // --only for inner pixels
     27    Float_t fSigmabarOuter;     // --only for outer pixels 
     28
    2029    Float_t fSigmabarSector[6]; // --for the 6 sectors of the camera
    2130    Float_t fSigmabarInnerSector[6];
    2231    Float_t fSigmabarOuterSector[6];
    23     Float_t fSigmabarInner; // --only for inner pixels
    24     Float_t fSigmabarOuter; // --only for outer pixels 
    25     UInt_t  fInnerPixels; // Overall number of inner pixels
    26     UInt_t  fOuterPixels; // Overall number of outer pixels
    27     Float_t fRatioA; // Ratio of areas (outer/inner pixels)
    28     Bool_t fCalcPixNum;
    2932
     33    UInt_t  fInnerPixels;       // Overall number of inner pixels
     34    UInt_t  fOuterPixels;       // Overall number of outer pixels
     35
     36    Bool_t  fCalcPixNum;
     37   
    3038public:
    3139
     
    4654    //    void SetSigmabarOuter(Float_t f) { fSigmabarOuter = f; }   
    4755
    48     Bool_t Calc(const MGeomCam &geom, const MPedestalCam &ped);
     56    Float_t Calc(const MGeomCam &geom, const MPedestalCam &ped, const MCerPhotEvt &evt);
    4957     
    5058    ClassDef(MSigmabar, 1)  // Storage Container for Sigmabar
     
    5361#endif
    5462
     63
     64
     65
     66
  • trunk/MagicSoft/Mars/manalysis/MSigmabarCalc.cc

    r1737 r1748  
    110110    if (!fMcEvt)
    111111      {
    112         *fLog << err << dbginf << "MHSigmabarTheta : MMcEvt not found... aborting." << endl;
     112        *fLog << err << dbginf << "MHSigmabaCalc : MMcEvt not found... aborting." << endl;
     113        return kFALSE;
     114      }
     115
     116    fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt");
     117    if (!fEvt)
     118      {
     119        *fLog << err << dbginf << "MHSigmabarCalc : MCerPhotEvt not found... aborting." << endl;
    113120        return kFALSE;
    114121      }
     
    126133Bool_t MSigmabarCalc::Process()
    127134{
    128   Bool_t rc = fSig->Calc(*fCam, *fPed);   
    129   fSigmabarMax = TMath::Max((Double_t)fSig->GetSigmabar(), fSigmabarMax);
    130   fSigmabarMin = TMath::Min((Double_t)fSig->GetSigmabar(), fSigmabarMin);
     135  Float_t rc = fSig->Calc(*fCam, *fPed, *fEvt);   
     136  fSigmabarMax = TMath::Max((Double_t)rc, fSigmabarMax);
     137  fSigmabarMin = TMath::Min((Double_t)rc, fSigmabarMin);
    131138
    132   if (fMcEvt->GetTelescopeTheta()*kRad2Deg < 90)
     139  if (fMcEvt->GetTelescopeTheta()*kRad2Deg < 120)
    133140    fThetaMax    = TMath::Max(fMcEvt->GetTelescopeTheta()*kRad2Deg, fThetaMax);
    134141  fThetaMin    = TMath::Min(fMcEvt->GetTelescopeTheta()*kRad2Deg, fThetaMin);
    135142
    136   return rc;
     143  return kTRUE;
    137144}
    138145
  • trunk/MagicSoft/Mars/manalysis/MSigmabarCalc.h

    r1737 r1748  
    4343    MSigmabarParam *fSigParam;
    4444    MMcEvt *fMcEvt;
     45    MCerPhotEvt *fEvt;
    4546    void Reset();
    4647
     
    5859#endif
    5960
     61
     62
Note: See TracChangeset for help on using the changeset viewer.