Changeset 1951


Ignore:
Timestamp:
04/12/03 16:40:23 (21 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
19 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/manalysis/MApplyPadding.cc

    r1713 r1951  
    7878// Default constructor.
    7979//
    80 MApplyPadding::MApplyPadding(const char *name, const char *title) : fRunType(0), fGroup(0), fUseHistogram(kTRUE), fFixedSigmabar(0.0)
     80MApplyPadding::MApplyPadding(const char *name, const char *title) : fRunType(0), fGroup(0), fUseHistogram(kTRUE), fFixedSigmabar(0.0), fRnd(0)
    8181{
    8282  fName  = name  ? name  : "MApplyPadding";
    … …  
    8787// --------------------------------------------------------------------------
    8888//
    89 // Destructor.
    90 //
    91 MApplyPadding::~MApplyPadding()
    92 {
    93   //nothing yet
    94 }
    95 
    96 // --------------------------------------------------------------------------
    97 //
    9889// You can provide a TH1D* histogram containing the target Sigmabar in
    9990// bins of theta. Be sure to use the same binning as for the analysis
    10091//
    101 Bool_t MApplyPadding::SetDefiningHistogram(TH1D *histo)
     92void MApplyPadding::SetDefiningHistogram(TH1D *histo)
    10293{
    10394  fHSigmabarMax = histo;
    104   return kTRUE;
    10595}
    10696
    … …  
    113103Bool_t MApplyPadding::PreProcess(MParList *pList)
    114104{
    115   fRnd = new TRandom3(0);
    116 
    117105  fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
    118106  if (!fMcEvt)
    … …  
    163151   if ((!fUseHistogram) && (fHSigmabarMax==NULL)) {
    164152     
    165      fHSigmabarMax = new TH1D();
     153     fHSigmabarMax = new TH1D;
    166154     fHSigmabarMax->SetNameTitle("fHSigmabarMax","Sigmabarmax for this analysis");
    167      TAxis &x = *fHSigmabarMax->GetXaxis();
    168 #if ROOT_VERSION_CODE < ROOT_VERSION(3,03,03)
    169      TString xtitle = x.GetTitle();
    170 #endif
    171      fHSigmabarMax->SetBins(binstheta->GetNumBins(), 0, 1);
    172      // Set the binning of the current histogram to the binning
    173      // in one of the two given histograms
    174      x.Set(binstheta->GetNumBins(), binstheta->GetEdges());
    175 #if ROOT_VERSION_CODE < ROOT_VERSION(3,03,03)
    176      x.SetTitle(xtitle);
    177 #endif
    178      
     155
     156     MH::SetBinning(fHSigmabarMax, binstheta);
     157
    179158     // -------------------------------------------------
    180159     // read in SigmabarParams
    … …  
    229208
    230209  // Get sigmabar which we have to pad to
    231   Double_t otherSig;
     210  Double_t otherSig = fFixedSigmabar;
    232211  if (fUseHistogram) {
    233     Int_t binNumber = fHSigmabarMax->GetXaxis()->FindBin(fMcEvt->GetTheta()*kRad2Deg);
     212    Int_t binNumber = fHSigmabarMax->GetXaxis()->FindBin(fMcEvt->GetTelescopeTheta()*kRad2Deg);
    234213    otherSig = fHSigmabarMax->GetBinContent(binNumber);
    235   } else {
    236     otherSig = fFixedSigmabar;
    237214  }
    238215
    239216  // Determine quadratic difference other-mine
    240   Double_t quadraticDiff = otherSig*otherSig - mySig*mySig;
    241 
     217  const Double_t quadraticDiff = otherSig*otherSig - mySig*mySig;
     218
     219  // crosscheck, should never happen
    242220  if (quadraticDiff < 0) {
    243221    *fLog << err << dbginf << "Event has higher Sigmabar="<<mySig<<" than Sigmabarmax="<<otherSig << " ...Skipping this event" <<endl;
    … …  
    248226 
    249227  // Pad if quadratic difference > 0
    250   if (quadraticDiff > 0) {   
    251 
    252228   MPedestalCam newPed;
    253229   newPed.InitSize(fPed->GetSize());
    … …  
    257233
    258234   for (UInt_t i=0; i<npix; i++) {
    259      MCerPhotPix pix = fEvt->operator[](i);     
     235     MCerPhotPix &pix = (*fEvt)[i];
    260236     if (!pix.IsPixelUsed())
    261237       continue;
    262238     pix.SetNumPhotons(pix.GetNumPhotons() +
    263239                       sqrt(quadraticDiff)*
    264                        fRnd->Gaus(0.0, 1.0)/
     240                       fRnd.Gaus(0.0, 1.0)/
    265241                       fCam->GetPixRatio(pix.GetPixId())
    266242                       );
    … …  
    268244     Double_t error = pix.GetErrorPhot();
    269245     pix.SetErrorPhot(sqrt(error*error + quadraticDiff));
    270      MPedestalPix ppix = fPed->operator[](i);
    271      MPedestalPix npix;
    272      npix.SetSigma(sqrt(ppix.GetSigma()*ppix.GetSigma() + quadraticDiff));
    273      newPed[i]=npix;
     246     MPedestalPix &ppix = (*fPed)[i];
     247     newPed[i].SetSigma(sqrt(ppix.GetSigma()*ppix.GetSigma() + quadraticDiff));
    274248    } //for
    275249   // Calculate Sigmabar again and crosscheck
    … …  
    278252   // fTest->Fill(otherSig,mySig);
    279253   return kTRUE;
    280   } //if
    281   return kFALSE;
    282 }
    283 
    284 Bool_t MApplyPadding::PostProcess()
    285 {
    286   //  fTest->DrawClone();
    287   return kTRUE;
    288 }
     254}
     255

  • trunk/MagicSoft/Mars/manalysis/MApplyPadding.h

    r1682 r1951  
    2424    MCerPhotEvt *fEvt;
    2525    MSigmabar *fSigmabar;
    26     TRandom3 *fRnd;
     26    TRandom3 fRnd;
    2727    Int_t fRunType;
    2828    Int_t fGroup;
    … …  
    3737public:
    3838    MApplyPadding(const char *name=NULL, const char *title=NULL);
    39     ~MApplyPadding();
    4039
    4140    Bool_t PreProcess(MParList *pList);
    4241    Bool_t Process();
    43     Bool_t PostProcess();
    4442   
    4543    void SetRunType(Int_t runtype) { fRunType =  runtype; }
    … …  
    4745    void SetDatabaseFile(char *filename) { fDatabaseFilename = filename; }
    4846    void SetTargetLevel(Double_t sigmabar) { fFixedSigmabar = sigmabar; fUseHistogram=kFALSE; }
    49     Bool_t SetDefiningHistogram(TH1D *histo);
     47    void SetDefiningHistogram(TH1D *histo);
    5048
    5149    ClassDef(MApplyPadding, 1)    // task for applying padding

  • trunk/MagicSoft/Mars/manalysis/MCerPhotPix.h

    r1715 r1951  
    1515    Bool_t  fIsCore;      // the pixel is a Core pixel          --> kTRUE
    1616
     17    UShort_t fRing;       // NT: number of analyzed rings around the core pixels
    1718    Float_t fPhot;        // The number of Cerenkov photons
    1819    Float_t fErrPhot;     // the error of fPhot
    … …  
    3132
    3233    Bool_t  IsPixelUsed() const      { return fIsUsed;   }
    33     void    SetPixelUnused()         { fIsUsed = kFALSE; }
    34     void    SetPixelUsed()           { fIsUsed = kTRUE;  }
     34    void    SetPixelUnused()         { fIsUsed = kFALSE; fRing=0; }
     35    void    SetPixelUsed()           { fIsUsed = kTRUE;  fRing=1;  }
     36
     37    void    SetRing(Short_t r)       { fRing = r; fIsUsed = (r>0); }
     38
     39    Short_t GetRing() const          { return fRing;}
    3540
    3641    void    SetPixelCore()           { fIsCore = kTRUE; }

  • trunk/MagicSoft/Mars/manalysis/MPadSchweizer.cc

    r1885 r1951  
    145145   if (!fPed)
    146146     {
    147        *fLog << dbginf << "MPedestalCam not found... aborting." << endl;
     147       *fLog << err << "MPedestalCam not found... aborting." << endl;
    148148       return kFALSE;
    149149     }
    … …  
    152152   if (!fCam)
    153153     {
    154        *fLog << dbginf << "MGeomCam not found (no geometry information available)... aborting." << endl;
     154       *fLog << err << "MGeomCam not found (no geometry information available)... aborting." << endl;
    155155       return kFALSE;
    156156     }
    … …  
    159159   if (!fEvt)
    160160     {
    161        *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl;
     161       *fLog << err << "MCerPhotEvt not found... aborting." << endl;
    162162       return kFALSE;
    163163     }
    … …  
    166166   if (!fSigmabar)
    167167     {
    168        *fLog << dbginf << "MSigmabar not found... aborting." << endl;
     168       *fLog << err << "MSigmabar not found... aborting." << endl;
    169169       return kFALSE;
    170170     }
    171171   
    172172   // Get Theta Binning 
    173    MBinning* binstheta  = (MBinning*)pList->FindObject("BinningTheta");
     173   MBinning* binstheta = (MBinning*)pList->FindObject("BinningTheta");
    174174   if (!binstheta)
    175175     {
    176        *fLog << err << dbginf << "BinningTheta not found... aborting." << endl;
     176       *fLog << err << "BinningTheta not found... aborting." << endl;
    177177       return kFALSE;     
    178178     }
    179179
    180180   // Get Sigma Binning 
    181    MBinning* binssigma  = (MBinning*)pList->FindObject("BinningSigmabar");
     181   MBinning* binssigma = (MBinning*)pList->FindObject("BinningSigmabar");
    182182   if (!binssigma)
    183183     {
    184        *fLog << err << dbginf << "BinningSigmabar not found... aborting." << endl;
     184       *fLog << err << "BinningSigmabar not found... aborting." << endl;
    185185       return kFALSE;     
    186186     }
    187187
    188188   // Get binning for (sigma^2-sigmabar^2) 
    189    MBinning* binsdiff  = (MBinning*)pList->FindObject("BinningDiffsigma2");
     189   MBinning* binsdiff = (MBinning*)pList->FindObject("BinningDiffsigma2");
    190190   if (!binsdiff)
    191191     {
    192        *fLog << err << dbginf << "BinningDiffsigma2 not found... aborting." << endl;
     192       *fLog << err << "BinningDiffsigma2 not found... aborting." << endl;
    193193       return kFALSE;     
    194194     }

  • trunk/MagicSoft/Mars/manalysis/MPadding.cc

    r1885 r1951  
    1616!
    1717!
    18 !   Author(s): Robert Wagner   <mailto:magicsoft@rwagner.de> 10/2002
    19 !              Wolfgang Wittek <mailto:wittek@mppmu.mpg.de>  01/2003
     18!   Author(s): Robert Wagner, 10/2002   <mailto:magicsoft@rwagner.de>
     19!   Author(s): Wolfgang Wittek, 01/2003 <mailto:wittek@mppmu.mpg.de>
     20!   Author(s): Thomas Bretz, 04/2003    <mailto:tbretz@astro.uni-wuerzburg.de>
    2021!
    2122!   Copyright: MAGIC Software Development, 2000-2003
    … …  
    6162//  missing. It is not the FINAL MAGIC VERSION.                            //
    6263//                                                                         //
     64//  For random numbers gRandom is used.                                    //
     65//                                                                         //
    6366/////////////////////////////////////////////////////////////////////////////
    6467#include "MPadding.h"
    … …  
    6669#include <stdio.h>
    6770
    68 #include "TH1.h"
    69 #include "TH2.h"
    70 #include "TH3.h"
    71 #include "TProfile.h"
    72 #include "TRandom.h"
    73 #include "TCanvas.h"
    74 
     71#include <TH1.h>
     72#include <TH2.h>
     73#include <TH3.h>
     74#include <TRandom.h>
     75#include <TCanvas.h>
     76#include <TProfile.h>
     77
     78#include "MH.h"
    7579#include "MBinning.h"
     80
    7681#include "MSigmabar.h"
     82
    7783#include "MMcEvt.hxx"
     84
     85#include "MParList.h"
     86
    7887#include "MLog.h"
    7988#include "MLogManip.h"
    80 #include "MParList.h"
     89
    8190#include "MGeomCam.h"
     91
     92#include "MCerPhotEvt.h"
    8293#include "MCerPhotPix.h"
    83 #include "MCerPhotEvt.h"
     94
     95#include "MPedestalCam.h"
    8496#include "MPedestalPix.h"
    8597
    … …  
    90102// Default constructor.
    91103//
    92 MPadding::MPadding(const char *name, const char *title) : fRunType(0), fGroup(0), fFixedSigmabar(0.0)
    93 {
    94   fName  = name  ? name  : "MPadding";
    95   fTitle = title ? title : "Task for the padding";
    96 
    97   fFixedSigmabar = 0.0;
    98   fHSigmabarMax  = NULL;
    99   fHSigmaTheta   = NULL;
    100 
    101   Print();
     104MPadding::MPadding(const char *name, const char *title)
     105    : fRunType(0), fGroup(0), fFixedSigmabar(0), fHSigMaxAllocated(kFALSE), fHSigmabarMax(NULL), fHSigmaTheta(NULL)
     106{
     107    fName  = name  ? name  : "MPadding";
     108    fTitle = title ? title : "Task for the padding";
     109
     110    //--------------------------------------------------------------------
     111    // plot pedestal sigmas for testing purposes
     112    fHSigmaPedestal = new TH2D("SigPed", "Padded vs orig. sigma",
     113                               100, 0.0, 5.0, 100, 0.0, 5.0);
     114   fHSigmaPedestal->SetXTitle("Orig. Pedestal sigma");
     115   fHSigmaPedestal->SetYTitle("Padded Pedestal sigma");
     116
     117   // plot no.of photons (before vs. after padding) for testing purposes
     118   fHPhotons = new TH2D("Photons", "Photons after vs.before padding",
     119                        100, -10.0, 90.0, 100, -10, 90);
     120   fHPhotons->SetXTitle("No.of photons before padding");
     121   fHPhotons->SetYTitle("No.of photons after padding");
     122
     123   // plot of added NSB
     124   fHNSB = new TH1D("NSB", "Distribution of added NSB", 100, -10.0, 10.0);
     125   fHNSB->SetXTitle("No.of added NSB photons");
     126   fHNSB->SetYTitle("No.of pixels");
     127
     128   fHSigmaOld = new TH2D;
     129   fHSigmaOld->SetNameTitle("fHSigmaOld", "Sigma before padding");
     130   fHSigmaOld->SetXTitle("Theta");
     131   fHSigmaOld->SetYTitle("Sigma");
     132
     133   fHSigmaNew = new TH2D;
     134   fHSigmaNew->SetNameTitle("fHSigmaNew", "Sigma after padding");
     135   fHSigmaNew->SetXTitle("Theta");
     136   fHSigmaNew->SetYTitle("Sigma");
    102137}
    103138
    … …  
    108143MPadding::~MPadding()
    109144{
    110   //nothing yet
     145   delete fHSigmaPedestal;
     146   delete fHPhotons;
     147   delete fHNSB;
     148   delete fHSigmaOld;
     149   delete fHSigmaNew;
     150   if (fHSigMaxAllocated && fHSigmabarMax)
     151       delete fHSigmabarMax;
    111152}
    112153
    … …  
    119160Bool_t MPadding::SetDefiningHistogram(TH1D *histo)
    120161{
    121   fHSigmabarMax  = histo;
    122 
    123   fFixedSigmabar = 0.0;
    124   fHSigmaTheta   = NULL;
    125 
    126   *fLog << "MPadding::SetDefiningHistogram" << endl;
    127 
    128   return kTRUE;
     162    if (fHSigmabarMax)
     163    {
     164        *fLog << warn << "MPadding - SigmabarMax already set.";
     165        return kFALSE;
     166    }
     167
     168    fHSigmabarMax  = histo;
     169
     170    fFixedSigmabar = 0;
     171    fHSigmaTheta   = NULL;
     172
     173    *fLog << inf << "MPadding - Use Defining Histogram.";
     174
     175    return kTRUE;
    129176}
    130177
    … …  
    136183Bool_t MPadding::SetSigmaThetaHist(TH2D *histo)
    137184{
    138   fHSigmaTheta   = histo;
    139 
    140   fFixedSigmabar = 0.0;
    141   fHSigmabarMax  = NULL;
    142 
    143   *fLog << "MPadding::SetSigmaThetaHist" << endl;
    144 
    145   return kTRUE;
     185    if (fHSigmaTheta)
     186    {
     187        *fLog << warn << "MPadding - SigmaTheta already set.";
     188        return kFALSE;
     189    }
     190
     191    fHSigmaTheta   = histo;
     192
     193    fFixedSigmabar = 0;
     194    if (fHSigMaxAllocated)
     195    {
     196        fHSigMaxAllocated = kFALSE;
     197        delete fHSigmabarMax;
     198    }
     199    fHSigmabarMax  = NULL;
     200
     201    *fLog << inf << "MPadding - Use Sigma Theta Histogram.";
     202
     203    return kTRUE;
    146204}
    147205
    … …  
    151209void MPadding::SetTargetLevel(Double_t sigmabar)
    152210{
    153   fFixedSigmabar = sigmabar;
    154 
    155   fHSigmabarMax  = NULL;
    156   fHSigmaTheta   = NULL;
    157 
    158   *fLog << "MPadding::SetTargetLevel; use fixed sigmabar : fFixedSigmabar = "
    159         << fFixedSigmabar << endl;
     211    fFixedSigmabar = sigmabar;
     212
     213    fHSigmaTheta   = NULL;
     214    if (fHSigMaxAllocated)
     215    {
     216        fHSigMaxAllocated = kFALSE;
     217        delete fHSigmabarMax;
     218    }
     219    fHSigmabarMax = NULL;
     220
     221    *fLog << inf << "MPadding - Use fixed sigmabar: fFixedSigmabar = ";
     222    *fLog << fFixedSigmabar << endl;
    160223}
    161224
    … …  
    167230Bool_t MPadding::PreProcess(MParList *pList)
    168231{
    169   fRnd = new TRandom3(0);
    170 
    171232  fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
    172233  if (!fMcEvt)
    173234    {
    174        *fLog << err << dbginf << "MPadding::PreProcess : MMcEvt not found... aborting." << endl;
     235       *fLog << err << dbginf << "MMcEvt not found... aborting." << endl;
    175236       return kFALSE;
    176237     }
    … …  
    179240   if (!fPed)
    180241     {
    181        *fLog << dbginf << "MPadding::PreProcess : MPedestalCam not found... aborting." << endl;
     242       *fLog << err << dbginf << "MPedestalCam not found... aborting." << endl;
    182243       return kFALSE;
    183244     }
    … …  
    186247   if (!fCam)
    187248     {
    188        *fLog << dbginf << "MPadding::PreProcess : MGeomCam not found (no geometry information available)... aborting." << endl;
     249       *fLog << err << dbginf << "MGeomCam not found (no geometry information available)... aborting." << endl;
    189250       return kFALSE;
    190251     }
    … …  
    193254   if (!fEvt)
    194255     {
    195        *fLog << dbginf << "MPadding::PreProcess : MCerPhotEvt not found... aborting." << endl;
     256       *fLog << err << dbginf << "MCerPhotEvt not found... aborting." << endl;
    196257       return kFALSE;
    197258     }
    … …  
    200261   if (!fSigmabar)
    201262     {
    202        *fLog << dbginf << "MPadding::PreProcess : MSigmabar not found... aborting." << endl;
     263       *fLog << err << dbginf << "MSigmabar not found... aborting." << endl;
    203264       return kFALSE;
    204265     }
    … …  
    208269   if (!binstheta)
    209270     {
    210        *fLog << err << dbginf << "MPadding::PreProcess : BinningTheta not found... aborting." << endl;
     271       *fLog << err << dbginf << "BinningTheta not found... aborting." << endl;
    211272       return kFALSE;     
    212273     }
    … …  
    216277   if (!binssigma)
    217278     {
    218        *fLog << err << dbginf << "MPadding::PreProcess : BinningSigmabar not found... aborting." << endl;
     279       *fLog << err << dbginf << "BinningSigmabar not found... aborting." << endl;
    219280       return kFALSE;     
    220281     }
    221282
    222    //--------------------------------------------------------------------
    223    // plot pedestal sigmas for testing purposes
    224    fHSigmaPedestal = new TH2D("SigPed","padded vs orig. sigma",
    225                      100, 0.0, 5.0, 100, 0.0, 5.0);
    226    fHSigmaPedestal->SetXTitle("orig. Pedestal sigma");
    227    fHSigmaPedestal->SetYTitle("padded Pedestal sigma");
    228 
    229    // plot no.of photons (before vs. after padding) for testing purposes
    230    fHPhotons = new TH2D("Photons","Photons after vs.before padding",
    231                         100, -10.0, 90.0, 100, -10, 90);
    232    fHPhotons->SetXTitle("no.of photons before padding");
    233    fHPhotons->SetYTitle("no.of photons after padding");
    234 
    235    // plot of added NSB
    236    fHNSB = new TH1D("NSB","Distribution of added NSB", 100, -10.0, 10.0);
    237    fHNSB->SetXTitle("no.of added NSB photons");
    238    fHNSB->SetYTitle("no.of pixels");
    239 
    240    fHSigmaOld = new TH2D();
    241    fHSigmaOld->SetNameTitle("fHSigmaOld","Sigma before padding");
    242283   MH::SetBinning(fHSigmaOld, binstheta, binssigma);
    243    fHSigmaOld->SetXTitle("Theta");
    244    fHSigmaOld->SetYTitle("Sigma");
    245 
    246    fHSigmaNew = new TH2D();
    247    fHSigmaNew->SetNameTitle("fHSigmaNew","Sigma after padding");
    248284   MH::SetBinning(fHSigmaNew, binstheta, binssigma);
    249    fHSigmaNew->SetXTitle("Theta");
    250    fHSigmaNew->SetYTitle("Sigma");
    251 
    252285
    253286   //************************************************************************
    … …  
    256289   // provided)
    257290   //
    258    if ( (fFixedSigmabar==0.0)   && (fHSigmabarMax==NULL)
    259         && (fHSigmaTheta==NULL) )
     291   if (fFixedSigmabar==0 && !fHSigmabarMax && !fHSigmaTheta)
    260292   {
    261         *fLog << "MPadding::PreProcess() : fFixedSigmabar, fHSigmabarMax = "
    262           <<  fFixedSigmabar << ",  " << fHSigmabarMax << endl;
    263      *fLog << "          create fSigmabarMax histogram" << endl;
    264 
    265      fHSigmabarMax = new TH1D();
    266      fHSigmabarMax->SetNameTitle("fHSigmabarMax","Sigmabarmax for this analysis");
    267      TAxis &x = *fHSigmabarMax->GetXaxis();
    268      fHSigmabarMax->SetBins(binstheta->GetNumBins(), 0, 1);
    269      // Set the binning of the current histogram to the binning
    270      // in one of the two given histograms
    271      x.Set(binstheta->GetNumBins(), binstheta->GetEdges());
    272      x.SetTitle("Theta");
    273      TAxis &y = *fHSigmabarMax->GetYaxis();     
    274      y.SetTitle("SigmabarMax");
     293       *fLog << inf << "MPadding - Creating fSigmabarMax histogram: ";
     294       *fLog << "fFixedSigmabar=" << fFixedSigmabar << ", ";
     295       *fLog << "fHSigmabarMax = " << fHSigmabarMax << endl;
     296
     297       // FIXME: Not deleted
     298     fHSigmabarMax = new TH1D;
     299     fHSigmabarMax->SetNameTitle("fHSigmabarMax", "Sigmabarmax for this analysis");
     300
     301     fHSigMaxAllocated = kTRUE;
     302
     303     MH::SetBinning(fHSigmabarMax, binstheta);
    275304
    276305     // -------------------------------------------------
    … …  
    279308     // -------------------------------------------------
    280309     
    281      FILE *f;
    282      if( !(f =fopen(fDatabaseFilename, "r")) ) {
    283        *fLog << err << dbginf << "Database file " << fDatabaseFilename
    284              << "was not found... (specify with MPadding::SetDatabaseFile) aborting." << endl;
    285        return kFALSE; 
    286      }
    287      char line[80];   
    288      Float_t sigmabarMin, sigmabarMax, thetaMin, thetaMax, ra, dec2;
    289      Int_t type, group, mjd, nr;
     310     FILE *f=fopen(fDatabaseFilename, "r");
     311     if(!f) {
     312         *fLog << err << dbginf << "Database file '" << fDatabaseFilename;
     313         *fLog << "' was not found (specified by MPadding::SetDatabaseFile) ...aborting." << endl;
     314         return kFALSE;
     315     }
     316
     317     TAxis &axe = *fHSigmabarMax->GetXaxis();
     318
     319     char line[80];
    290320     while ( fgets(line, sizeof(line), f) != NULL) {
    291        if ((line[0]!='#')) {   
    292          sscanf(line,"%d %d %f %f %d %d %f %f %f %f",
    293                 &type, &group, &ra, &dec2, &mjd, &nr,
    294                 &sigmabarMin,&sigmabarMax,&thetaMin,&thetaMax);
    295          if ((group==fGroup)||(type==1)) //selected ON group or OFF         
    296            {   
    297              // find out which bin(s) we have to look at
    298              for (Int_t i=fHSigmabarMax->GetXaxis()->FindBin(thetaMin);
    299                i<fHSigmabarMax->GetXaxis()->FindBin(thetaMax)+1; i++)
    300              if (sigmabarMax > fHSigmabarMax->GetBinContent(i))
    301                fHSigmabarMax->SetBinContent(i, sigmabarMax);     
    302            }
    303        }
     321         if (line[0]=='#')
     322             continue;
     323
     324         Float_t sigmabarMin, sigmabarMax, thetaMin, thetaMax, ra, dec2;
     325         Int_t type, group, mjd, nr;
     326
     327         sscanf(line,"%d %d %f %f %d %d %f %f %f %f",
     328                &type, &group, &ra, &dec2, &mjd, &nr,
     329                &sigmabarMin, &sigmabarMax, &thetaMin, &thetaMax);
     330
     331         if (group!=fGroup && type!=1) //selected ON group or OFF
     332             continue;
     333
     334         const Int_t from = axe.FindFixBin(thetaMin);
     335         const Int_t to   = axe.FindFixBin(thetaMax);
     336
     337         // find out which bin(s) we have to look at
     338         for (Int_t i=from; i<to+1; i++)
     339             if (sigmabarMax > fHSigmabarMax->GetBinContent(i))
     340                 fHSigmabarMax->SetBinContent(i, sigmabarMax);
    304341     }//while
    305342 
    … …  
    307344   //************************************************************************
    308345
     346   if (!fHSigmabarMax && !fHSigmaTheta && fFixedSigmabar==0)
     347   {
     348       *fLog << err << "ERROR: Sigmabar for padding not defined... aborting." << endl;
     349       return kFALSE;
     350   }
     351
    309352   return kTRUE;
    310353}
     354
     355Double_t MPadding::CalcOtherSig(const Double_t mySig, const Double_t theta) const
     356{
     357  //
     358  // Get sigmabar which we have to pad to
     359  //
     360  const TAxis   &axe     = *fHSigmabarMax->GetXaxis();
     361  const Int_t    binnum  =  axe.FindFixBin(theta);
     362  const Bool_t   inrange =  theta>=axe.GetXmin() && theta<=axe.GetXmax();
     363
     364  if ((fHSigmabarMax || fHSigmaTheta) && !inrange)
     365  {
     366      *fLog << err << dbginf;
     367      *fLog << "Theta of current event is beyond the limits, Theta = ";
     368      *fLog << theta << " ...skipping." <<endl;
     369      return -1;
     370  }
     371
     372
     373  //
     374  // get target sigma for the current Theta from the histogram fHSigmabarMax
     375  //
     376  if (fHSigmabarMax != NULL)
     377      return fHSigmabarMax->GetBinContent(binnum);
     378
     379  //
     380  // for the current Theta,
     381  // generate a sigma according to the histogram fHSigmaTheta
     382  //
     383  if (fHSigmaTheta != NULL)
     384  {
     385      Double_t otherSig = -1;
     386
     387      TH1D* fHSigma = fHSigmaTheta->ProjectionY("", binnum, binnum, "");
     388
     389      if (fHSigma->GetEntries()>0)
     390          otherSig = fHSigma->GetRandom();
     391
     392      delete fHSigma;
     393
     394      return otherSig;
     395  }
     396
     397  //
     398  // use a fixed target sigma
     399  //
     400  return fFixedSigmabar;
     401}
     402
     403// --------------------------------------------------------------------------
     404//
     405// Do the padding  (mySig ==> otherSig)
     406//
     407Bool_t MPadding::Padding(const Double_t quadraticDiff, const Double_t theta)
     408{
     409   const UInt_t npix = fEvt->GetNumPixels();
     410
     411   // pad only pixels   - which are used (before image cleaning)
     412   //                   - and for which the no.of photons is != 0.0
     413   //
     414   for (UInt_t i=0; i<npix; i++)
     415   {   
     416     MCerPhotPix &pix = (*fEvt)[i];
     417     if ( !pix.IsPixelUsed() )
     418       continue;
     419/*
     420     if ( pix.GetNumPhotons() == 0)
     421     {
     422       *fLog << "MPadding::Process(); no.of photons is 0 for used pixel"
     423             << endl;
     424       continue;
     425     }
     426*/
     427     const Double_t area = fCam->GetPixRatio(pix.GetPixId());
     428
     429     // add additional NSB to the number of photons
     430     const Double_t NSB = sqrt(quadraticDiff*area) * gRandom->Gaus(0, 1);
     431     const Double_t oldphotons = pix.GetNumPhotons();
     432     const Double_t newphotons = oldphotons + NSB;
     433     pix.SetNumPhotons( newphotons );
     434
     435     fHNSB->Fill( NSB/sqrt(area) );
     436     fHPhotons->Fill( newphotons/sqrt(area), oldphotons/sqrt(area) );
     437
     438     // error: add sigma of padded noise quadratically
     439     const Double_t olderror = pix.GetErrorPhot();
     440     const Double_t newerror = sqrt( olderror*olderror + quadraticDiff*area );
     441     pix.SetErrorPhot( newerror );
     442
     443     MPedestalPix &ppix = (*fPed)[i];
     444
     445     ppix.SetMeanRms(0);
     446
     447     const Double_t oldsigma = ppix.GetMeanRms();
     448     const Double_t newsigma = sqrt( oldsigma*oldsigma + quadraticDiff*area );
     449     ppix.SetMeanRms( newsigma );
     450
     451     fHSigmaPedestal->Fill( oldsigma, newsigma );
     452     fHSigmaOld->Fill( theta, oldsigma );
     453     fHSigmaNew->Fill( theta, newsigma );
     454   } //for
     455
     456   return kTRUE;
     457}
     458
    311459
    312460// --------------------------------------------------------------------------
    … …  
    321469Bool_t MPadding::Process()
    322470{
    323   //-------------------------------------------
    324   // Calculate sigmabar of event
    325   //
    326   Double_t mySig = fSigmabar->Calc(*fCam, *fPed, *fEvt);
    327   //fSigmabar->Print("");
    328 
    329   //$$$$$$$$$$$$$$$$$$$$$$$$$$
    330   mySig = 0.0;
    331   //$$$$$$$$$$$$$$$$$$$$$$$$$$
    332 
    333 
    334   // Get sigmabar which we have to pad to
    335   //
    336   Double_t otherSig;
    337   Double_t Theta  = kRad2Deg*fMcEvt->GetTelescopeTheta();
    338 
    339   // *fLog << "Theta = " << Theta << endl;
    340 
    341   //-------------------------------------------
    342   // get target sigma for the current Theta from the histogram fHSigmabarMax
    343   //
    344 
    345   if (fHSigmabarMax != NULL)
    346   {
    347     Int_t binNumber = fHSigmabarMax->GetXaxis()->FindBin(Theta);
    348     if (binNumber < 1  || binNumber > fHSigmabarMax->GetNbinsX())
    349     {
    350       *fLog << err << dbginf
    351             << "Theta of current event is beyond the limits, Theta = " 
    352             << kRad2Deg*fMcEvt->GetTelescopeTheta()
    353             << "   ...Skipping this event" <<endl;
    354       return kCONTINUE;
    355     }
    356     else
    357     {
    358       otherSig = fHSigmabarMax->GetBinContent(binNumber);
    359 
    360       //*fLog << "Theta, binNumber, otherSig = "
    361       //      << kRad2Deg*fMcEvt->GetTelescopeTheta() << ",  "
    362       //      << binNumber << ",  " << otherSig << endl;
    363     }
    364   }
    365 
    366   //-------------------------------------------
    367   // for the current Theta,
    368   // generate a sigma according to the histogram fHSigmaTheta
    369   //
    370   else if (fHSigmaTheta != NULL)
    371   {
    372     Int_t binNumber = fHSigmaTheta->GetXaxis()->FindBin(Theta);
    373 
    374     if ( binNumber < 1  ||  binNumber > fHSigmaTheta->GetNbinsX() )
    375     {
    376       //      *fLog << "MPadding::Process(); binNumber out of range, binNumber = "
    377       //      << binNumber << ",  Skip event " << endl;
    378       return kCONTINUE;
    379     }
    380     else
    381     {
    382       TH1D* fHSigma =
    383             fHSigmaTheta->ProjectionY("", binNumber, binNumber, "");
    384       if ( fHSigma->GetEntries() == 0.0 )
    385       {
    386         //*fLog << "MPadding::Process(); projection for Theta bin " << binNumber
    387         //      << " has no entries,  Skip event" << endl;
     471    const Double_t theta = kRad2Deg*fMcEvt->GetTelescopeTheta();
     472
     473    //
     474    // Calculate sigmabar of event
     475    //
     476    Double_t mySig = fSigmabar->Calc(*fCam, *fPed, *fEvt);
     477
     478    //$$$$$$$$$$$$$$$$$$$$$$$$$$
     479    mySig = 0.0;  // FIXME?
     480    //$$$$$$$$$$$$$$$$$$$$$$$$$$
     481
     482    const Double_t otherSig = CalcOtherSig(mySig, theta);
     483
     484    // Skip event if target sigma is zero
     485    if (otherSig<=0)
    388486        return kCONTINUE;
    389       }
    390       else
    391       {
    392         otherSig = fHSigma->GetRandom();
    393 
    394         //*fLog << "Theta, bin number = " << Theta << ",  " << binNumber
    395         //      << ",  otherSig = " << otherSig << endl;
    396       }
    397       delete fHSigma;
    398     }
    399   }
    400 
    401   //-------------------------------------------
    402   // use a fixed target sigma
    403   //
    404   else if (fFixedSigmabar != 0.0)
    405   {
    406     otherSig = fFixedSigmabar;
    407   }
    408 
    409   //-------------------------------------------
    410   else
    411   {
    412     *fLog << "MPadding::Process(); sigmabar for padding not defined" << endl;
    413     return kFALSE;
    414   }
    415 
    416   //-------------------------------------------
    417   //
    418 
    419   //*fLog << "MPadding::Process(); mySig, otherSig = " << mySig << ",  "
    420   //      << otherSig << endl;
    421 
    422   // Skip event if target sigma is zero
    423   if (otherSig == 0.0)
    424   {
    425     return kCONTINUE;     
    426   }
    427 
    428   // Determine quadratic difference other-mine
    429   Double_t quadraticDiff = otherSig*otherSig - mySig*mySig;
    430 
    431   if (quadraticDiff < 0) {
    432     *fLog << err << dbginf << "Event has higher Sigmabar=" << mySig
    433           <<" than Sigmabarmax=" << otherSig << "; Theta ="
    434           << kRad2Deg*fMcEvt->GetTelescopeTheta()
    435           << "   ...Skipping this event" <<endl;
    436     return kCONTINUE; //skip
    437   }
    438 
    439   if (quadraticDiff == 0) return kTRUE; //no padding necessary.
    440  
    441 
    442   //-------------------------------------------
    443   // quadratic difference is > 0, do the padding;
    444   //
    445   Padding(quadraticDiff);
    446 
    447   // Calculate Sigmabar again and crosscheck
    448   mySig = fSigmabar->Calc(*fCam, *fPed, *fEvt);
    449 
    450   //fSigmabar->Print("");
    451 
    452   return kTRUE;
    453 }
    454 
    455 // --------------------------------------------------------------------------
    456 //
    457 // Do the padding  (mySig ==> otherSig)
    458 //
    459 Bool_t MPadding::Padding(Double_t quadraticDiff)
    460 {
    461    const UInt_t npix = fEvt->GetNumPixels();
    462 
    463    // pad only pixels   - which are used (before image cleaning)
    464    //                   - and for which the no.of photons is != 0.0
    465    //
    466    for (UInt_t i=0; i<npix; i++)
    467    {   
    468      MCerPhotPix &pix = fEvt->operator[](i);     
    469      if ( !pix.IsPixelUsed() )
    470        continue;
    471 
    472      if ( pix.GetNumPhotons() == 0.0)
    473      {
    474        *fLog << "MPadding::Process(); no.of photons is 0 for used pixel"
    475              << endl;
    476        continue;
    477      }
    478 
    479      Int_t j = pix.GetPixId();
    480      Double_t Area = fCam->GetPixRatio(j);
    481 
    482      // add additional NSB to the number of photons
    483      Double_t oldphotons = pix.GetNumPhotons();
    484      Double_t NSB = sqrt(quadraticDiff*Area) * fRnd->Gaus(0.0, 1.0);
    485      Double_t newphotons = oldphotons + NSB;
    486      pix.SetNumPhotons( newphotons );
    487 
    488      fHNSB->Fill( NSB/sqrt(Area) );
    489      fHPhotons->Fill( newphotons/sqrt(Area), oldphotons/sqrt(Area) );
    490 
    491 
    492      // error: add sigma of padded noise quadratically
    493      Double_t olderror = pix.GetErrorPhot();
    494      Double_t newerror = sqrt(  olderror*olderror + quadraticDiff*Area );
    495      pix.SetErrorPhot( newerror );
    496 
    497 
    498      MPedestalPix &ppix = fPed->operator[](j);
    499 
    500      //$$$$$$$$$$$$$$$$$$$$$$$$$$
    501      ppix.SetMeanRms(0.0);
    502      //$$$$$$$$$$$$$$$$$$$$$$$$$$
    503 
    504      Double_t oldsigma = ppix.GetMeanRms();
    505      Double_t newsigma = sqrt(  oldsigma*oldsigma + quadraticDiff*Area );
    506      ppix.SetMeanRms( newsigma );
    507 
    508      fHSigmaPedestal->Fill( oldsigma, newsigma );
    509      fHSigmaOld->Fill( kRad2Deg*fMcEvt->GetTelescopeTheta(), oldsigma );
    510      fHSigmaNew->Fill( kRad2Deg*fMcEvt->GetTelescopeTheta(), newsigma );
    511    } //for
    512 
    513    return kTRUE;
    514 }
    515 
    516 // --------------------------------------------------------------------------
    517 //
     487
     488    // Determine quadratic difference other-mine
     489    const Double_t quadraticDiff = otherSig*otherSig - mySig*mySig;
     490
     491    if (quadraticDiff < 0) {
     492        *fLog << err << "ERROR - MPadding: Event has higher Sigmabar=" << mySig;
     493        *fLog << " than Sigmabarmax=" << otherSig << " @ Theta =" << theta;
     494        *fLog << " ...skipping." << endl;
     495        return kCONTINUE; //skip
     496    }
     497
     498    if (quadraticDiff == 0)
     499        return kTRUE; //no padding necessary.
     500
     501    //
     502    // quadratic difference is > 0, do the padding;
     503    //
     504    Padding(quadraticDiff, theta);
     505
     506    // Calculate Sigmabar again and crosscheck
     507    //mySig = fSigmabar->Calc(*fCam, *fPed, *fEvt);
     508
     509    return kTRUE;
     510}
     511
     512// --------------------------------------------------------------------------
     513//
     514// Draws some histograms if IsGraphicalOutputEnabled
    518515//
    519516Bool_t MPadding::PostProcess()
    520517{
    521     TCanvas &c = *(MH::MakeDefCanvas("Padding", "", 600, 900));
     518    if (!IsGraphicalOutputEnabled())
     519        return kTRUE;
     520
     521    TCanvas &c = *MH::MakeDefCanvas("Padding", "", 600, 900);
    522522    c.Divide(2,3);
    523523    gROOT->SetSelectedPad(NULL);
    524524
    525 
    526525    if (fHSigmabarMax != NULL)
    527526    {
    528       c.cd(1);
    529       fHSigmabarMax->DrawClone();     
    530       fHSigmabarMax->SetBit(kCanDelete);     
     527        c.cd(1);
     528        fHSigmabarMax->DrawClone();
    531529    }
    532530    else if (fHSigmaTheta != NULL)
    533531    {
    534       c.cd(1);
    535       fHSigmaTheta->DrawClone();     
    536       fHSigmaTheta->SetBit(kCanDelete);     
     532        c.cd(1);
     533        fHSigmaTheta->DrawClone();
    537534    }
    538535
    539536    c.cd(3);
    540537    fHSigmaPedestal->DrawClone();
    541     fHSigmaPedestal->SetBit(kCanDelete);   
    542538
    543539    c.cd(5);
    544540    fHPhotons->DrawClone();
    545     fHPhotons->SetBit(kCanDelete);   
    546541
    547542    c.cd(2);
    548543    fHNSB->DrawClone();
    549     fHNSB->SetBit(kCanDelete);   
    550544
    551545    c.cd(4);
    552546    fHSigmaOld->DrawClone();     
    553     fHSigmaOld->SetBit(kCanDelete);     
    554547
    555548    c.cd(6);
    556549    fHSigmaNew->DrawClone();     
    557     fHSigmaNew->SetBit(kCanDelete);     
    558 
    559 
    560   return kTRUE;
    561 }
    562 
    563 // --------------------------------------------------------------------------
    564 
    565 
    566 
    567 
    568 
    569 
     550
     551    return kTRUE;
     552}

  • trunk/MagicSoft/Mars/manalysis/MPadding.h

    r1768 r1951  
    66#endif
    77
    8 #ifndef MARS_MH
    9 #include "MH.h"
    10 #endif
    11 
    12 #include "TRandom3.h"
    13 #include "TH1.h"
    14 #include "TH2.h"
    15 #include "TH3.h"
    16 #include "TProfile.h"
    17 
    18 
     8class TH1D;
     9class TH2D;
    1910class MGeomCam;
    2011class MCerPhotEvt;
    … …  
    2819{
    2920private:
    30     MGeomCam *fCam;
    31     MCerPhotEvt    *fEvt;
    32     MSigmabar      *fSigmabar;
    33     MMcEvt         *fMcEvt;
    34     MPedestalCam   *fPed;
     21    MGeomCam     *fCam;
     22    MCerPhotEvt  *fEvt;
     23    MSigmabar    *fSigmabar;
     24    MMcEvt       *fMcEvt;
     25    MPedestalCam *fPed;
    3526
    36     TRandom3       *fRnd;
     27    Int_t     fRunType;
     28    Int_t     fGroup;
    3729
    38     Int_t          fRunType;
    39     Int_t          fGroup;
     30    TString   fDatabaseFilename; // data file used for generating fHSigmabarMax histogram
     31    Double_t  fFixedSigmabar;    // fixed sigmabar value
    4032
    41     TH1D           *fHSigmabarMax;   // histogram (sigmabarmax vs. Theta)
    42     char           *fDatabaseFilename; // data file used for generating
    43                                      //   fHSigmabarMax histogram
     33    Bool_t    fHSigMaxAllocated; // flag whether MPadding allocated it
     34    TH1D     *fHSigmabarMax;     // histogram (sigmabarmax vs. Theta)
     35    TH2D     *fHSigmaTheta;      // 2D-histogram (sigmabar vs. Theta)
     36    TH2D     *fHSigmaPedestal;   //-> for testing: plot of padded vs orig. pedestal sigmas
     37    TH2D     *fHPhotons;         //-> for testing: no.of photons after versus before padding
     38    TH2D     *fHSigmaOld;        //-> histogram (sigma vs. Theta) before padding
     39    TH2D     *fHSigmaNew;        //-> histogram (sigma vs. Theta) after padding
     40    TH1D     *fHNSB;             //-> histogram of added NSB
    4441
    45     TH2D           *fHSigmaTheta;    // 2D-histogram (sigmabar vs. Theta)
    46 
    47     Double_t       fFixedSigmabar;   // fixed sigmabar value
    48 
    49     TH2D           *fHSigmaPedestal; // for testing : plot of padded vs
    50                                      //               orig. pedestal sigmas
    51     TH2D           *fHPhotons;       // for testing : no.of photons after
    52                                      //               versus before padding
    53     TH2D           *fHSigmaOld;      // histogram (sigma vs. Theta)
    54                                      // before padding
    55 
    56     TH2D           *fHSigmaNew   ;   // histogram (sigma vs. Theta)
    57                                      // after padding
    58     TH1D           *fHNSB;           // histogram of added NSB
    59 
     42    Double_t CalcOtherSig(const Double_t mySig, const Double_t theta) const;
     43    Bool_t   Padding(Double_t quadDiff, const Double_t theta);
    6044
    6145public:
    … …  
    6751    Bool_t PostProcess();
    6852   
    69     Bool_t Padding(Double_t quadDiff);
    70 
    7153    void SetRunType(Int_t runtype) { fRunType =  runtype; }
    7254    void SetGroup(Int_t group)     { fGroup   =  group; }
    … …  
    7961    void SetTargetLevel(Double_t sigmabar);
    8062
    81     ClassDef(MPadding, 1)    // task for the padding
     63    ClassDef(MPadding, 0)   // task for the padding
    8264};
    8365

  • trunk/MagicSoft/Mars/manalysis/MParameters.h

    r1884 r1951  
    3131    Int_t GetVal() const { return fVal; }
    3232
    33     ClassDef(MParameterI, 1) // Container to hold a generalized parameters (double)
     33    ClassDef(MParameterI, 1) // Container to hold a generalized parameters (integer)
    3434};
     35/*
     36class MParameters : public MParContainer
     37{
     38private:
     39    TObjArray fList;
     40    TObjArray fNames;
    3541
     42public:
     43    MParameters(const char *name=NULL, const char *title=NULL)
     44    {
     45        fName  = name  ? name  : "MParameters";
     46        fTitle = title ? title : "Additional temporary parameters";
     47
     48        SetReadyToSave();
     49    }
     50
     51    MParamaterI &AddInteger(const TString name, const TString title, Int_t val=0)
     52    {
     53        MParameterI &p = *new MParameterI(name, title);
     54        p.SetValue(val);
     55
     56        fList.Add(&p);
     57
     58        TNamed &n = *new TNamed(name, title);
     59        fNames.Add(&n);
     60
     61        return p;
     62    }
     63
     64    MParameterD &AddDouble(const TString name, const TString title, Double_t val=0)
     65    {
     66        MParameterD &p = *new MParameterD(name, title);
     67        p.SetValue(val);
     68
     69        fList.Add(&p);
     70
     71        TNamed &n = *new TNamed(name, title);
     72        fNames.Add(&n);
     73
     74        return p;
     75    }
     76
     77    const TObjArray &GetList()
     78    {
     79        fList.SetNames(&fNames);
     80        return fList;
     81    }
     82
     83    MParameterD *GetParameterD(const TString &name)
     84    {
     85        fList.SetNames(&fNames);
     86        return (MParamaterD*)fList.FindObject(name);
     87    }
     88
     89    MParameterI *GetParameterI(const TString &name)
     90    {
     91        fList.SetNames(&fNames);
     92        return (MParameterI*)fList.FindObject(name);
     93    }
     94
     95    ClassDef(MParameters, 1) // List to hold generalized parameters (MParameterD/I)
     96    }
     97    */
    3698#endif

  • trunk/MagicSoft/Mars/manalysis/MPointingCorr.h

    r1888 r1951  
    1818class MParameterD;
    1919
    20 
    2120class MPointingCorr : public MTask
    2221{
    … …  
    2726    MParameterD  *fHourAngle;
    2827
    29     Float_t      fMm2Deg;
     28    Float_t       fMm2Deg;
    3029
    3130public:

  • trunk/MagicSoft/Mars/manalysis/MSigmabar.cc

    r1885 r1951  
    1616!
    1717!
    18 !   Author(s): Robert Wagner   10/2002 <mailto:magicsoft@rwagner.de>
    19 !              Wolfgang Wittek 01/2003 <mailto:wittek@mppmu.mpg.de>
     18!   Author(s): Robert Wagner, 10/2002 <mailto:magicsoft@rwagner.de>
     19!   Author(s): Wolfgang Wittek, 01/2003 <mailto:wittek@mppmu.mpg.de>
     20!   Author(s): Thomas Bretz, 04/2003 <mailto:tbretz@astro.uni-wuerzburg.de>
    2021!
    2122!   Copyright: MAGIC Software Development, 2003
    … …  
    4041#include "MLog.h"
    4142#include "MLogManip.h"
     43
    4244#include "MParList.h"
     45
    4346#include "MGeomCam.h"
     47#include "MGeomPix.h"
     48
     49#include "MCerPhotEvt.h"
     50#include "MCerPhotPix.h"
     51
    4452#include "MPedestalCam.h"
    4553#include "MPedestalPix.h"
    46 #include "MGeomPix.h"
    47 #include "MCerPhotEvt.h"
    48 #include "MCerPhotPix.h"
    4954
    5055ClassImp(MSigmabar);
    … …  
    5257// --------------------------------------------------------------------------
    5358//
    54 MSigmabar::MSigmabar(const char *name, const char *title)
     59MSigmabar::MSigmabar(const char *name, const char *title) : fCalcPixNum(kTRUE)
    5560{
    5661    fName  = name  ? name  : "MSigmabar";
    5762    fTitle = title ? title : "Storage container for Sigmabar";
    58    
    59 
    60     fCalcPixNum=kTRUE;
    6163}
    6264
    … …  
    6668{
    6769  // do nothing special.
     70}
     71
     72void MSigmabar::Reset()
     73{
     74    fSigmabar      = -1;
     75    fInnerPixels   = -1;
     76    fOuterPixels   = -1;
     77    fSigmabarInner = -1;
     78    fSigmabarOuter = -1;
     79
     80    memset(fSigmabarSector, 0, sizeof(fSigmabarSector));
     81    memset(fSigmabarInnerSector, 0, sizeof(fSigmabarInnerSector));
     82    memset(fSigmabarOuterSector, 0, sizeof(fSigmabarOuterSector));
    6883}
    6984
    … …  
    8196                        const MCerPhotEvt &evt)
    8297{
    83   fSigmabar      = 0.0;
    84   fSigmabarInner = 0.0;
    85   fSigmabarOuter = 0.0;
    86 
    87 
    88   Float_t innerSquaredSum[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
    89   Float_t outerSquaredSum[6] = {0.0, 0.0, 0.0, 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 
    93   // sum up sigma**2 for each sector, separately for inner and outer region;
    94   // all pixels are renormalized to the area of pixel 0
    95   //
    96   // consider all pixels with Cherenkov photon information
    97   // and require "Used"
    98   //
    99 
    100   const UInt_t npix = evt.GetNumPixels();
    101 
    102   for (UInt_t i=0; i<npix; i++)
    103   {
    104       MCerPhotPix &cerpix = evt.operator[](i);
    105       if (!cerpix.IsPixelUsed())
    106         continue;
    107  
    108       if ( cerpix.GetNumPhotons() == 0 )
    109       {
    110         *fLog << "MSigmabar::Calc(); no.of photons is 0 for used pixel"
    111               << endl;
    112         continue;
    113       }
    114 
    115       Int_t j = cerpix.GetPixId();
    116       Double_t area = geom.GetPixRatio(j);
    117 
    118       const MGeomPix    &gpix = geom[j];
    119       const MPedestalPix &pix =  ped[j];
    120 
    121       //angle = 6.0*atan2(gpix.GetX(),gpix.GetY()) / (2.0*TMath::Pi()) - 0.5;
    122       //if (angle<0.0) angle+=6.0;
    123 
    124       Float_t angle = atan2(gpix.GetY(),gpix.GetX())*6 / (TMath::Pi()*2);
    125       if (angle<0) angle+=6;
    126 
    127       Int_t currentSector=(Int_t)angle;
    128        
    129       // count only those pixels which have a sigma != 0.0
    130       Float_t sigma = pix.GetMeanRms();
    131 
    132       if ( sigma != 0.0 )
    133       { 
     98    Int_t innerPixels[6];
     99    Int_t outerPixels[6];
     100    Float_t innerSquaredSum[6];
     101    Float_t outerSquaredSum[6];
     102
     103    memset(innerPixels, 0, sizeof(innerPixels));
     104    memset(outerPixels, 0, sizeof(outerPixels));
     105    memset(outerSquaredSum, 0, sizeof(outerSquaredSum));
     106    memset(outerSquaredSum, 0, sizeof(outerSquaredSum));
     107
     108    //
     109    // sum up sigma^2 for each sector, separately for inner and outer region;
     110    // all pixels are renormalized to the area of pixel 0
     111    //
     112    // consider all pixels with Cherenkov photon information
     113    // and require "Used"
     114    //
     115
     116    const UInt_t npix = evt.GetNumPixels();
     117
     118    for (UInt_t i=0; i<npix; i++)
     119    {
     120        MCerPhotPix &cerpix = evt.operator[](i);
     121        if (!cerpix.IsPixelUsed())
     122            continue;
     123        /*
     124         if ( cerpix.GetNumPhotons() == 0 )
     125         {
     126         *fLog << "MSigmabar::Calc(); no.of photons is 0 for used pixel"
     127         << endl;
     128         continue;
     129         }
     130         */
     131        const Int_t id = cerpix.GetPixId();
     132        const Double_t area = geom.GetPixRatio(id);
     133
     134        const MGeomPix    &gpix = geom[id];
     135        const MPedestalPix &pix =  ped[id];
     136
     137        Int_t sector = (Int_t)(atan2(gpix.GetY(),gpix.GetX())*6 / (TMath::Pi()*2));
     138        if (sector<0)
     139            sector+=6;
     140
     141        // count only those pixels which have a sigma != 0.0
     142        const Float_t sigma = pix.GetMeanRms();
     143
     144        if ( sigma <= 0 )
     145            continue;
     146
    134147        if (area < 1.5)
    135148        {
    136           innerPixels[currentSector]++;
    137           innerSquaredSum[currentSector]+= sigma*sigma / area;
     149            innerPixels[sector]++;
     150            innerSquaredSum[sector]+= sigma*sigma / area;
    138151        }
    139152        else
    140153        {
    141           outerPixels[currentSector]++;
    142           outerSquaredSum[currentSector]+= sigma*sigma / area;
     154            outerPixels[sector]++;
     155            outerSquaredSum[sector]+= sigma*sigma / area;
    143156        }
    144       }
    145   }
    146  
    147   fSigmabarInner=0; fInnerPixels=0;
    148   fSigmabarOuter=0; fOuterPixels=0;
    149   for (UInt_t i=0; i<6; i++) {
    150     fSigmabarInner += innerSquaredSum[i];
    151     fInnerPixels   += innerPixels[i];
    152     fSigmabarOuter += outerSquaredSum[i];
    153     fOuterPixels   += outerPixels[i];
    154   }
    155 
    156 
    157   // this is the sqrt of the average sigma**2;
     157    }
     158
     159    fInnerPixels   = 0;
     160    fOuterPixels   = 0;
     161    fSigmabarInner = 0;
     162    fSigmabarOuter = 0;
     163    for (UInt_t i=0; i<6; i++) {
     164        fSigmabarInner += innerSquaredSum[i];
     165        fInnerPixels   += innerPixels[i];
     166        fSigmabarOuter += outerSquaredSum[i];
     167        fOuterPixels   += outerPixels[i];
     168    }
     169
     170    //
     171    // this is the sqrt of the average sigma^2;
     172    //
     173    fSigmabar=fInnerPixels+fOuterPixels<=0?0:sqrt((fSigmabarInner+fSigmabarOuter)/(fInnerPixels+fOuterPixels));
     174
     175    if (fInnerPixels > 0) fSigmabarInner /= fInnerPixels;
     176    if (fOuterPixels > 0) fSigmabarOuter /= fOuterPixels;
     177
    158178  //
    159   if ( (fInnerPixels+fOuterPixels) > 0)
    160     fSigmabar=sqrt(   (fSigmabarInner + fSigmabarOuter)
    161                     / (fInnerPixels   + fOuterPixels) );
    162 
    163   if (fInnerPixels > 0) fSigmabarInner /= fInnerPixels;
    164   if (fOuterPixels > 0) fSigmabarOuter /= fOuterPixels;
    165  
    166   // this is the sqrt of the average sigma**2
     179  // this is the sqrt of the average sigma^2
    167180  // for the inner and outer pixels respectively
    168181  //
    169182  fSigmabarInner = sqrt( fSigmabarInner );
    170183  fSigmabarOuter = sqrt( fSigmabarOuter );
    171  
    172184
    173185  for (UInt_t i=0; i<6; i++) {
    174     fSigmabarSector[i]         = 0.0;
    175     fSigmabarInnerSector[i]    = 0.0;
    176     fSigmabarOuterSector[i]    = 0.0;
    177 
    178     if ( (innerPixels[i]+outerPixels[i]) > 0)
    179       fSigmabarSector[i] = sqrt(  (innerSquaredSum[i] + outerSquaredSum[i])
    180                                 / (innerPixels[i]     + outerPixels[i]    ) );
    181     if ( innerPixels[i] > 0)
    182       fSigmabarInnerSector[i] = innerSquaredSum[i] / innerPixels[i];
    183     if ( outerPixels[i] > 0)
    184       fSigmabarOuterSector[i] = outerSquaredSum[i] / outerPixels[i];
    185 
    186 
    187     fSigmabarInnerSector[i] = sqrt( fSigmabarInnerSector[i] );
    188     fSigmabarOuterSector[i] = sqrt( fSigmabarOuterSector[i] );
     186      fSigmabarSector[i] = innerPixels[i]+outerPixels[i]<=0?0:sqrt((innerSquaredSum[i]+outerSquaredSum[i])/(innerPixels[i]+outerPixels[i]));
     187
     188      const Double_t is = innerPixels[i]<=0?0:innerSquaredSum[i]/innerPixels[i];
     189      const Double_t os = outerPixels[i]<=0?0:outerSquaredSum[i]/outerPixels[i];
     190
     191      fSigmabarInnerSector[i] = sqrt( is );
     192      fSigmabarOuterSector[i] = sqrt( os );
    189193  }
    190194   
    … …  
    196200void MSigmabar::Print(Option_t *) const
    197201{
    198   *fLog << endl;
     202  *fLog << all << endl;
    199203  *fLog << "Total number of inner pixels is " << fInnerPixels << endl;
    200204  *fLog << "Total number of outer pixels is " << fOuterPixels << endl;
    … …  
    203207  *fLog << "Sigmabar     Overall : " << fSigmabar << "   ";
    204208  *fLog << " Sectors: ";
    205   for (Int_t i=0;i<6;i++) *fLog << fSigmabarSector[i] << ", ";
     209  for (Int_t i=0;i<6;i++)
     210      *fLog << fSigmabarSector[i] << ", ";
    206211  *fLog << endl;
    207212
    … …  
    209214  *fLog << "Sigmabar     Inner   : " << fSigmabarInner << "   ";
    210215  *fLog << " Sectors: ";
    211   for (Int_t i=0;i<6;i++) *fLog << fSigmabarInnerSector[i] << ", ";
     216  for (Int_t i=0;i<6;i++)
     217      *fLog << fSigmabarInnerSector[i] << ", ";
    212218  *fLog << endl;
    213219
    … …  
    215221  *fLog << "Sigmabar     Outer   : " << fSigmabarOuter << "   ";
    216222  *fLog << " Sectors: ";
    217   for (Int_t i=0;i<6;i++) *fLog << fSigmabarOuterSector[i] << ", ";
    218   *fLog << endl;
    219 
    220 }
    221 // --------------------------------------------------------------------------
    222 
    223 
    224 
    225 
     223  for (Int_t i=0;i<6;i++)
     224      *fLog << fSigmabarOuterSector[i] << ", ";
     225  *fLog << endl;
     226
     227}

  • trunk/MagicSoft/Mars/manalysis/MSigmabar.h

    r1748 r1951  
    66#endif
    77
    8 #ifndef MARS_MParList
    9 #include "MParList.h"
    10 #endif
    11 
    12 #ifndef MARS_MGeomCam
    13 #include "MGeomCam.h"
    14 #endif
    15 
    16 #ifndef MARS_MPedestalCam
    17 #include "MPedestalCam.h"
    18 #endif
    19 
     8class MGeomCam;
     9class MParList;
    2010class MCerPhotEvt;
     11class MPedestalCam;
    2112
    2213class MSigmabar : public MParContainer
    … …  
    3122    Float_t fSigmabarOuterSector[6];
    3223
    33     UInt_t  fInnerPixels;       // Overall number of inner pixels
    34     UInt_t  fOuterPixels;       // Overall number of outer pixels
     24    Int_t  fInnerPixels;       // Overall number of inner pixels
     25    Int_t  fOuterPixels;       // Overall number of outer pixels
    3526
    3627    Bool_t  fCalcPixNum;
    … …  
    4132    ~MSigmabar();
    4233   
     34    void Reset();
     35
    4336    void Print(Option_t *) const;
    4437 

  • trunk/MagicSoft/Mars/manalysis/MSigmabarCalc.cc

    r1748 r1951  
    1616!
    1717!
    18 !   Author(s): Robert Wagner <rwagner@mppmu.mpg.de> 10/2002
     18!   Author(s): Robert Wagner, 10/2002 <rwagner@mppmu.mpg.de>
     19!   Author(s): Thomas Bretz, 4/2003 <tbretz@astro.uni-wuerzburg.de>
    1920!
    2021!   Copyright: MAGIC Software Development, 2000-2003
    … …  
    3940
    4041#include "MParList.h"
     42
    4143#include "MGeomCam.h"
    4244#include "MPedestalCam.h"
     45
    4346#include "MSigmabar.h"
    4447#include "MSigmabarParam.h"
     48
    4549#include "MMcEvt.hxx"
    4650
    … …  
    5357MSigmabarCalc::MSigmabarCalc(const char *name, const char *title)
    5458{
    55   fName  = name  ? name  : "MSigmabarCalc";
    56   fTitle = title ? title : "Task to calculate Sigmabar";
     59    fName  = name  ? name  : "MSigmabarCalc";
     60    fTitle = title ? title : "Task to calculate Sigmabar";
    5761
    58   Reset();
     62    Reset();
    5963}
    6064
    … …  
    7478    if (!fCam)
    7579    {
    76         *fLog << dbginf << "MGeomCam not found (no geometry information available)... aborting." << endl;
     80        *fLog << err << "MGeomCam not found (no geometry information available)... aborting." << endl;
    7781        return kFALSE;
    7882    }
    … …  
    8185    if (!fPed)
    8286    {
    83         *fLog << dbginf << "MPedestalCam not found... aborting." << endl;
     87        *fLog << err << "MPedestalCam not found... aborting." << endl;
    8488        return kFALSE;
    8589    }
    … …  
    8892    if (!fSig)
    8993    {
    90         *fLog << dbginf << "MSigmabar not found... aborting." << endl;
     94        *fLog << err << "MSigmabar not found... aborting." << endl;
    9195        return kFALSE;
    9296    }
    … …  
    9599    if (!fSigParam)
    96100    {
    97         *fLog << dbginf << "MSigmabarParam not found... aborting." << endl;
     101        *fLog << err << "MSigmabarParam not found... aborting." << endl;
    98102        return kFALSE;
    99103    }
    … …  
    102106    if (!fRun)
    103107    {
    104         *fLog << dbginf << "MRawRunHeader not found... aborting." << endl;
     108        *fLog << err << "MRawRunHeader not found... aborting." << endl;
    105109        return kFALSE;
    106110    }
    107111   
     112    // This is needed for determining min/max Theta
    108113    fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
    109     // This is needed for determining min/max Theta
    110114    if (!fMcEvt)
    111115      {
    112         *fLog << err << dbginf << "MHSigmabaCalc : MMcEvt not found... aborting." << endl;
     116        *fLog << err << "MMcEvt not found... aborting." << endl;
    113117        return kFALSE;
    114118      }
    … …  
    117121    if (!fEvt)
    118122      {
    119         *fLog << err << dbginf << "MHSigmabarCalc : MCerPhotEvt not found... aborting." << endl;
     123        *fLog << err << "MCerPhotEvt not found... aborting." << endl;
    120124        return kFALSE;
    121125      }
    … …  
    133137Bool_t MSigmabarCalc::Process()
    134138{
    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);
     139    const Double_t rc = fSig->Calc(*fCam, *fPed, *fEvt);
    138140
    139   if (fMcEvt->GetTelescopeTheta()*kRad2Deg < 120)
    140     fThetaMax    = TMath::Max(fMcEvt->GetTelescopeTheta()*kRad2Deg, fThetaMax);
    141   fThetaMin    = TMath::Min(fMcEvt->GetTelescopeTheta()*kRad2Deg, fThetaMin);
     141    fSigmabarMax = TMath::Max(rc, fSigmabarMax);
     142    fSigmabarMin = TMath::Min(rc, fSigmabarMin);
    142143
    143   return kTRUE;
     144    const Double_t theta = fMcEvt->GetTelescopeTheta()*kRad2Deg;
     145
     146    fThetaMax = TMath::Max(theta, fThetaMax);
     147    fThetaMin = TMath::Min(theta, fThetaMin);
     148
     149    return kTRUE;
    144150}
    145151
    … …  
    151157Bool_t MSigmabarCalc::ReInit(MParList *pList)
    152158{
    153    
    154   fSigParam->SetParams(1, fSigmabarMin, fSigmabarMax, fThetaMin, fThetaMax);
    155   fSigParam->SetRunNumber(fRun->GetRunNumber());
    156  
    157   Reset();
    158  
    159   return kTRUE;
     159    fSigParam->SetParams(1, fSigmabarMin, fSigmabarMax, fThetaMin, fThetaMax);
     160    fSigParam->SetRunNumber(fRun->GetRunNumber());
     161
     162    Reset();
     163
     164    return kTRUE;
    160165}
    161166
    162167void MSigmabarCalc::Reset()
    163168{
    164   fThetaMin = 200;    //there must be a function which gives me the hightest
    165   fThetaMax = 0;     // value allowed for a certain type!
    166   fSigmabarMin = 200;
    167   fSigmabarMax = 0;
     169    fThetaMin = 200;    // there must be a function which gives me the hightest
     170    fThetaMax = 0;      // value allowed for a certain type!
     171    fSigmabarMin = 200;
     172    fSigmabarMax = 0;
    168173}
    169174

  • trunk/MagicSoft/Mars/manalysis/MSigmabarCalc.h

    r1768 r1951  
    3333{
    3434private:
    35     MGeomCam      *fCam; 
    36     MPedestalCam  *fPed; 
    37     MRawRunHeader *fRun;
    38     MSigmabar    *fSig;
     35    MMcEvt         *fMcEvt;
     36    MCerPhotEvt    *fEvt;
     37    MGeomCam       *fCam;
     38    MPedestalCam   *fPed;
     39    MRawRunHeader  *fRun;
     40    MSigmabar      *fSig;
     41    MSigmabarParam *fSigParam;
     42
    3943    Double_t fSigmabarMin; // Parametrization
    4044    Double_t fSigmabarMax;
    4145    Double_t fThetaMin;
    4246    Double_t fThetaMax;
    43     MSigmabarParam *fSigParam;
    44     MMcEvt *fMcEvt;
    45     MCerPhotEvt *fEvt;
     47
    4648    void Reset();
    4749
    … …  
    5456    Bool_t Process();
    5557
    56     ClassDef(MSigmabarCalc, 2)    // task for calculating sigmabar
     58    ClassDef(MSigmabarCalc, 0) // task for calculating sigmabar
    5759};
    5860

  • trunk/MagicSoft/Mars/mhist/MHSigmaPixel.cc

    r1682 r1951  
    1717!
    1818!   Author(s): Robert Wagner 10/2002 <mailto:magicsoft@rwagner.de>
    19 !   Copyright: MAGIC Software Development, 2000-2002
     19!
     20!   Copyright: MAGIC Software Development, 2000-2003
    2021!
    2122!
    … …  
    7980  if (!fPedestalCam)
    8081    {
    81       *fLog << err << dbginf << "MHSigmaPixel: MPedestalCam not found... aborting." << endl;
     82      *fLog << err << "MPedestalCam not found... aborting." << endl;
    8283      return kFALSE;
    8384    }
    8485 
    8586  MBinning* binssigma = (MBinning*)plist->FindObject("BinningSigma");
    86   MBinning* binspixel = new MBinning();
    87   binspixel->SetEdges(fPedestalCam->GetSize(), -0.5, -0.5+fPedestalCam->GetSize());
    88  
    8987  if (!binssigma)
    9088    {
    91       *fLog << err << dbginf << "MHSigmaPixel: BinningSigma not found... aborting." << endl;
     89      *fLog << err << "BinningSigma [MBinning] not found... aborting." << endl;
    9290      return kFALSE;     
    93    }
     91    }
    9492
    95    SetBinning(&fHist, binspixel, binssigma);
     93  const Int_t n = fPedestalCam->GetSize();
    9694
    97    fHist.Sumw2();
     95  MBinning binspixel;
     96  binspixel.SetEdges(n, -0.5, -0.5+n);
    9897
    99    return kTRUE;
     98  SetBinning(&fHist, &binspixel, binssigma);
     99
     100  fHist.Sumw2();
     101
     102  return kTRUE;
    100103}
    101104
    … …  
    106109Bool_t MHSigmaPixel::Fill(const MParContainer *par)
    107110{
    108   MPedestalCam &ped = *(MPedestalCam*)par;
     111    const MPedestalCam &ped = *(MPedestalCam*)par;
    109112    for (Int_t i=0;i<(ped.GetSize());i++)
    110113    {
    111       const MPedestalPix pix = ped[i];       
    112       fHist.Fill(i, pix.GetSigma());
     114        const MPedestalPix pix = ped[i];
     115        fHist.Fill(i, pix.GetSigma());
    113116    }
    114117    return kTRUE;

  • trunk/MagicSoft/Mars/mhist/MHSigmaPixel.h

    r1682 r1951  
    1919{
    2020private:
    21     MPedestalCam *fPedestalCam;
     21    MPedestalCam *fPedestalCam; //!
     22
    2223    TH2D fHist;
    2324

  • trunk/MagicSoft/Mars/mhist/MHSigmaTheta.cc

    r1809 r1951  
    9696  if (!fMcEvt)
    9797    {
    98        *fLog << err << dbginf << "MHSigmaTheta::SetupFill : MMcEvt not found... aborting." << endl;
     98       *fLog << err << "MMcEvt not found... aborting." << endl;
    9999       return kFALSE;
    100100     }
    … …  
    103103   if (!fPed)
    104104     {
    105        *fLog << dbginf << "MHSigmaTheta::SetupFill : MPedestalCam not found... aborting." << endl;
     105       *fLog << err << "MPedestalCam not found... aborting." << endl;
    106106       return kFALSE;
    107107     }
    … …  
    110110   if (!fCam)
    111111     {
    112        *fLog << dbginf << "MHSigmaTheta::SetupFill : MGeomCam not found (no geometry information available)... aborting." << endl;
     112       *fLog << err << "MGeomCam not found (no geometry information available)... aborting." << endl;
    113113       return kFALSE;
    114114     }
    … …  
    117117   if (!fEvt)
    118118     {
    119        *fLog << dbginf << "MHSigmaTheta::SetupFill : MCerPhotEvt not found... aborting." << endl;
     119       *fLog << err << "MCerPhotEvt not found... aborting." << endl;
    120120       return kFALSE;
    121121     }
    … …  
    124124   if (!fSigmabar)
    125125     {
    126        *fLog << dbginf << "MHSigmaTheta::SetupFill : MSigmabar not found... aborting." << endl;
     126       *fLog << err << "MSigmabar not found... aborting." << endl;
    127127       return kFALSE;
    128128     }
    … …  
    132132   if (!binstheta)
    133133     {
    134        *fLog << err << dbginf << "MHSigmaTheta::SetupFill : BinningTheta not found... aborting." << endl;
     134       *fLog << err << "BinningTheta [MBinning] not found... aborting." << endl;
    135135       return kFALSE;     
    136136     }
    137137
    138138   // Get Sigmabar binning 
    139    MBinning* binssigma  = (MBinning*)plist->FindObject("BinningSigmabar");
     139   MBinning* binssigma = (MBinning*)plist->FindObject("BinningSigmabar");
    140140   if (!binssigma)
    141141     {
    142        *fLog << err << dbginf << "MHSigmaTheta::SetupFill : BinningSigmabar not found... aborting." << endl;
     142       *fLog << err << "BinningSigmabar [MBinning] not found... aborting." << endl;
    143143       return kFALSE;     
    144144     }
    … …  
    148148   if (!binsdiff)
    149149     {
    150        *fLog << err << dbginf << "MHSigmaTheta::SetupFill : BinningDiffsigma2 not found... aborting." << endl;
     150       *fLog << err << "BinningDiffsigma2 [MBinning] not found... aborting." << endl;
    151151       return kFALSE;     
    152152     }
    … …  
    154154
    155155   // Set binnings in histograms
    156    SetBinning(&fSigmaTheta,    binstheta, binssigma);
     156   SetBinning(&fSigmaTheta, binstheta, binssigma);
    157157
    158158   // Get binning for pixel number
    159    UInt_t npix = fPed->GetSize();
     159   const UInt_t npix = fPed->GetSize();
     160
    160161   MBinning binspix("BinningPixel");
    161    MBinning* binspixel = &binspix;
    162    binspixel->SetEdges(npix, -0.5, ((float)npix)-0.5 );
    163 
    164    SetBinning(&fSigmaPixTheta, binstheta, binspixel, binssigma);
    165    SetBinning(&fDiffPixTheta,  binstheta, binspixel, binsdiff);
     162   binspix.SetEdges(npix, -0.5, -0.5+npix );
     163
     164   SetBinning(&fSigmaPixTheta, binstheta, &binspix, binssigma);
     165   SetBinning(&fDiffPixTheta,  binstheta, &binspix, binsdiff);
    166166
    167167   return kTRUE;
    … …  
    176176  //*fLog << "entry Fill" << endl;
    177177
    178     Double_t Theta = fMcEvt->GetTelescopeTheta()*kRad2Deg;
     178    Double_t theta = fMcEvt->GetTelescopeTheta()*kRad2Deg;
    179179    Double_t mySig = fSigmabar->Calc(*fCam, *fPed, *fEvt);
    180     fSigmaTheta.Fill(Theta, mySig);
    181 
    182     //*fLog << "Theta, mySig = " << Theta << ",  " << mySig << endl;
     180
     181    fSigmaTheta.Fill(theta, mySig);
    183182
    184183    const UInt_t npix = fEvt->GetNumPixels();
    … …  
    189188        continue;
    190189
     190      /*
    191191      if (cerpix.GetNumPhotons() == 0)
    192192        continue;
    193 
    194       Int_t j = cerpix.GetPixId();
    195       const MPedestalPix pix = fPed->operator[](j);
    196 
    197       Double_t Sigma = pix.GetMeanRms();
    198       Double_t Area  = fCam->GetPixRatio(j);
    199 
    200       fSigmaPixTheta.Fill(Theta, (Double_t)j, Sigma);
    201 
    202       Double_t Diff = Sigma*Sigma/Area - mySig*mySig;
    203       fDiffPixTheta.Fill (Theta, (Double_t)j, Diff);
     193        */
     194
     195      const Int_t id = cerpix.GetPixId();
     196      const MPedestalPix &pix = (*fPed)[id];
     197
     198      const Double_t sigma = pix.GetMeanRms();
     199      const Double_t area  = fCam->GetPixRatio(id);
     200
     201      fSigmaPixTheta.Fill(theta, (Double_t)id, sigma);
     202
     203      const Double_t diff = sigma*sigma/area - mySig*mySig;
     204      fDiffPixTheta.Fill(theta, (Double_t)id, diff);
    204205    }
    205 
    206     return kTRUE;
    207 }
    208 
    209 // --------------------------------------------------------------------------
    210 //
    211 //  Plot the results
    212 //
    213 Bool_t MHSigmaTheta::Finalize()
    214 {
    215     DrawClone();
    216206
    217207    return kTRUE;
    … …  
    329319{
    330320    if (!gPad)
    331         MakeDefCanvas("SigmaTheta", "Sigmabar vs. Theta",
    332                        600, 600);
     321        MakeDefCanvas("SigmaTheta", "Sigmabar vs. Theta", 600, 600);
    333322
    334323    TH1D *h;
    … …  
    358347    fSigmaTheta.DrawCopy(opt);
    359348
    360 
    361 
    362349    gPad->Modified();
    363350    gPad->Update();

  • trunk/MagicSoft/Mars/mhist/MHSigmaTheta.h

    r1748 r1951  
    1414#endif
    1515
    16 #ifndef ROOT_TProfile2D
    17 #include "TProfile2D.h"
    18 #endif
    19 
    2016class MGeomCam;
    2117class MCerPhotEvt;
    … …  
    3026{
    3127private:
    32     const MGeomCam *fCam;
    33     MPedestalCam   *fPed;
    34     MCerPhotEvt    *fEvt;
    35     MSigmabar      *fSigmabar;
    36     MMcEvt         *fMcEvt;
     28    const MGeomCam *fCam;        //!
     29    MPedestalCam   *fPed;        //!
     30    MCerPhotEvt    *fEvt;        //!
     31    MSigmabar      *fSigmabar;   //!
     32    MMcEvt         *fMcEvt;      //!
    3733 
    38     TH2D           fSigmaTheta;   // 2D-distribution sigmabar versus Theta;
    39                                   // sigmabar is the average pedestasl sigma
    40                                   // in an event
    41     TH3D           fSigmaPixTheta;// 3D-distr.:Theta, pixel, pedestal sigma
    42     TH3D           fDiffPixTheta; // 3D-distr.:Theta, pixel, sigma^2-sigmabar^2
     34    TH2D fSigmaTheta;    // 2D-distribution sigmabar versus Theta; sigmabar is the average pedestasl sigma in an event
     35    TH3D fSigmaPixTheta; // 3D-distr.:Theta, pixel, pedestal sigma
     36    TH3D fDiffPixTheta;  // 3D-distr.:Theta, pixel, sigma^2-sigmabar^2
    4337
    4438
    … …  
    4640    MHSigmaTheta(const char *name=NULL, const char *title=NULL);
    4741
    48     virtual Bool_t SetupFill(const MParList *plist);
    49     virtual Bool_t Fill(const MParContainer *par);
    50     virtual Bool_t Finalize();
     42    Bool_t SetupFill(const MParList *plist);
     43    Bool_t Fill(const MParContainer *par);
    5144
    5245    const TH2D *GetSigmaTheta() { return &fSigmaTheta; }

  • trunk/MagicSoft/Mars/mhist/MHSigmabarTheta.cc

    r1682 r1951  
    7676   if (!fMcEvt)
    7777   {
    78        *fLog << err << dbginf << "MHSigmabarTheta : MMcEvt not found... aborting." << endl;
     78       *fLog << err << "MMcEvt not found... aborting." << endl;
    7979       return kFALSE;
    8080   }
    … …  
    8383   if (!fSigmabar)
    8484   {
    85        *fLog << err << dbginf << "MHSigmabarTheta : MSigmabar not found... aborting." << endl;
     85       *fLog << err << "MSigmabar not found... aborting." << endl;
    8686       return kFALSE;
    8787   }
    8888
     89   MBinning* binstheta  = (MBinning*)plist->FindObject("BinningTheta");
     90   if (!binstheta)
     91   {
     92       *fLog << err << "BinningTheta [MBinning] not found... aborting." << endl;
     93       return kFALSE;     
     94   }
    8995   MBinning* binssigmabar = (MBinning*)plist->FindObject("BinningSigmabar");
    90    MBinning* binstheta  = (MBinning*)plist->FindObject("BinningTheta");
    91    if (!binssigmabar || !binstheta)
     96   if (!binssigmabar)
    9297   {
    93        *fLog << err << dbginf << "MHSigmabarTheta : At least one MBinning not found... aborting." << endl;
     98       *fLog << err << "BinningSigmabar [MBinning] not found... aborting." << endl;
    9499       return kFALSE;     
    95100   }

  • trunk/MagicSoft/Mars/mhist/MHSigmabarTheta.h

    r1682 r1951  
    1919{
    2020private:
    21     MMcEvt    *fMcEvt;
    22     MSigmabar *fSigmabar;
     21    MMcEvt    *fMcEvt;        //!
     22    MSigmabar *fSigmabar;     //!
    2323
    2424    TH2D fHist;

  • trunk/MagicSoft/Mars/mhist/MHStarMap.cc

    r1867 r1951  
    156156    const float sind = sqrt(1.0-cosd*cosd);
    157157
    158 
    159     float t = h.GetMeanY() - m*h.GetMeanX();
     158    const float t = h.GetMeanY() - m*h.GetMeanX();
    160159
    161160    if (!fUseMmScale)
    … …  
    163162
    164163    // get step size ds along the main axis of the ellipse
    165     TAxis &axe = *fStarMap->GetXaxis();
     164    const TAxis &axe = *fStarMap->GetXaxis();
    166165    const int   N    = axe.GetNbins();
    167     const float xmin = axe.GetBinLowEdge(1);
    168     const float xmax = axe.GetBinLowEdge(N+1);
     166    const float xmin = axe.GetXmin();
     167    const float xmax = axe.GetXmax();
     168    // FIXME: Fixed number?
    169169    const float ds   = (xmax-xmin) / 200.0;
    170170
    171171    if (m>-1 && m<1)
    172172    {
    173         float dx = ds * cosd;
    174         float  x = xmin + dx/2.0;
    175         int   N1 = (int) ((xmax-xmin)/dx+1.0);
     173        const float dx = ds * cosd;
     174        const float  x = xmin + dx/2.0;
     175        const int   N1 = (int) ((xmax-xmin)/dx+1.0);
    176176
    177177        for (int i=0; i<N1; i++)
    … …  
    184184    else
    185185    {
    186         TAxis &axe = *fStarMap->GetYaxis();
     186        const TAxis &axe = *fStarMap->GetYaxis();
    187187        const int   M    = axe.GetNbins();
    188         const float ymin = axe.GetBinLowEdge(1);
    189         const float ymax = axe.GetBinLowEdge(M+1);
    190 
    191         float dy = ds * sind;
    192         float  y = ymin + dy/2.0;
    193         int   M1 = (int) ((ymax-ymin)/dy+1.0);
     188        const float ymin = axe.GetXmin();
     189        const float ymax = axe.GetXmax();
     190
     191        const float dy = ds * sind;
     192        const float  y = ymin + dy/2.0;
     193        const int   M1 = (int) ((ymax-ymin)/dy+1.0);
    194194
    195195        for (int i=0; i<M1; i++)
Note: See TracChangeset for help on using the changeset viewer.