Changeset 3889


Ignore:
Timestamp:
04/29/04 15:48:38 (21 years ago)
Author:
gaug
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/Changelog

    r3888 r3889  
    3535     - add the possibility to set ranges and extraction windows. Default
    3636       is what has always been
     37
     38
     39   * mjobs/MJPedestals.[h,cc]
     40     - add the possibility to set a signal extractor which gives the
     41       extraction ranges to MPedCalcPedRun
     42     - derive from MHCamDisplays
    3743
    3844

  • trunk/MagicSoft/Mars/mjobs/MGCamDisplays.cc

    r3875 r3889  
    121121// 3: Triple Gauss (for distributions with inner, outer pixels and outliers)
    122122// 4: flat         (for the probability distributions)
     123// (1-4:) Moreover, sectors 6,1 and 2 of the camera and sectors 3,4 and 5 are
     124//        drawn separately, for inner and outer pixels.
    123125// 5: Fit Inner and Outer pixels separately by a single Gaussian
    124126//                 (only for MAGIC cameras)
    125 //
    126 // Moreover, sectors 6,1 and 2 of the camera and sectors 3,4 and 5 are
    127 // drawn separately, for inner and outer pixels.
     127// 6: Fit Inner and Outer pixels separately by a single Gaussian and display
     128//                 additionally the two camera halfs separately (for MAGIC camera)
     129//
    128130//
    129131void MGCamDisplays::DrawProjection(MHCamera *obj, Int_t fit) const
    … …  
    136138  outer[0] = 1;
    137139         
    138   if (fit==5)
     140  if (fit==5 || fit==6)
    139141    {
    140142     
    … …  
    149151          s0[5] = 5;
    150152
     153          TArrayI s1(3);
     154          s1[0] = 6;
     155          s1[1] = 1;
     156          s1[2] = 2;
     157         
     158          TArrayI s2(3);
     159          s2[0] = 3;
     160          s2[1] = 4;
     161          s2[2] = 5;
     162
    151163          gPad->Clear();
    152164          TVirtualPad *pad = gPad;
    153165          pad->Divide(2,1);
    154166         
    155           TH1D *half[2];
    156           half[0] = obj->ProjectionS(s0, inner, "Inner");
    157           half[1] = obj->ProjectionS(s0, outer, "Outer");
    158          
    159           half[0]->SetDirectory(NULL);
    160           half[1]->SetDirectory(NULL);
     167          TH1D *inout[2];
     168          inout[0] = obj->ProjectionS(s0, inner, "Inner");
     169          inout[1] = obj->ProjectionS(s0, outer, "Outer");
     170         
     171          inout[0]->SetDirectory(NULL);
     172          inout[1]->SetDirectory(NULL);
    161173         
    162174          for (int i=0; i<2; i++)
    163175            {
    164176              pad->cd(i+1);
    165               half[i]->SetLineColor(kRed+i);
    166               half[i]->SetBit(kCanDelete);
    167               half[i]->Draw();
    168               half[i]->Fit("gaus","Q");
     177              inout[i]->SetLineColor(kRed+i);
     178              inout[i]->SetBit(kCanDelete);
     179              inout[i]->Draw();
     180              inout[i]->Fit("gaus","Q");
     181
     182              if (fit == 6)
     183                {
     184                  TH1D *half[2];
     185                  half[0] = obj->ProjectionS(s1, i==0 ? inner : outer , "Sector 6-1-2");
     186                  half[1] = obj->ProjectionS(s2, i==0 ? inner : outer , "Sector 3-4-5");
     187                 
     188                  for (int j=0; j<2; j++)
     189                    {
     190                      half[j]->SetLineColor(kRed+i+j);
     191                      half[j]->SetDirectory(0);
     192                      half[j]->SetBit(kCanDelete);
     193                      half[j]->Draw("same");
     194                    }
     195                }
     196             
    169197            }
    170198         
    171199          gLog << all << obj->GetName()
    172200               << Form("%s%5.3f%s%3.2f"," Mean: Inner Pixels: ",
    173                        half[0]->GetFunction("gaus")->GetParameter(1),"+-",
    174                        half[0]->GetFunction("gaus")->GetParError(1));
     201                       inout[0]->GetFunction("gaus")->GetParameter(1),"+-",
     202                       inout[0]->GetFunction("gaus")->GetParError(1));
    175203          gLog << Form("%s%5.3f%s%3.2f","  Outer Pixels: ",
    176                        half[1]->GetFunction("gaus")->GetParameter(1),"+-",
    177                        half[1]->GetFunction("gaus")->GetParError(1)) << endl;
     204                       inout[1]->GetFunction("gaus")->GetParameter(1),"+-",
     205                       inout[1]->GetFunction("gaus")->GetParError(1)) << endl;
    178206          gLog << all << obj->GetName()
    179207               << Form("%s%5.3f%s%3.2f"," Sigma: Inner Pixels: ",
    180                        half[0]->GetFunction("gaus")->GetParameter(2),"+-",
    181                        half[0]->GetFunction("gaus")->GetParError(2));
     208                       inout[0]->GetFunction("gaus")->GetParameter(2),"+-",
     209                       inout[0]->GetFunction("gaus")->GetParError(2));
    182210          gLog << Form("%s%5.3f%s%3.2f","  Outer Pixels: ",
    183                        half[1]->GetFunction("gaus")->GetParameter(2),"+-",
    184                        half[1]->GetFunction("gaus")->GetParError(2)) << endl;
    185          
     211                       inout[1]->GetFunction("gaus")->GetParameter(2),"+-",
     212                       inout[1]->GetFunction("gaus")->GetParError(2)) << endl;
     213
    186214        }
    187215      return;

  • trunk/MagicSoft/Mars/mjobs/MJPedestal.cc

    r3810 r3889  
    4545#include "MTaskList.h"
    4646#include "MEvtLoop.h"
     47#include "MExtractor.h"
    4748
    4849#include "MStatusDisplay.h"
    … …  
    6263
    6364using namespace std;
    64 
     65// --------------------------------------------------------------------------
     66//
     67// Default constructor.
     68//
     69// Sets fRuns to 0, fExtractor to NULL, fDataCheck to kFALSE
     70//
    6571MJPedestal::MJPedestal(const char *name, const char *title)
    66     : fRuns(0), fDataCheck(kFALSE)
    67 
     72    : fRuns(0), fExtractor(NULL), fDataCheck(kFALSE)
    6873{
    6974    fName  = name  ? name  : "MJPedestal";
    … …  
    115120}
    116121
    117 void MJPedestal::DrawProjection(MHCamera *obj1, Int_t fit) const
    118 {
    119     TH1D *obj2 = (TH1D*)obj1->Projection(obj1->GetName());
    120     obj2->SetDirectory(0);
    121     obj2->Draw();
    122     obj2->SetBit(kCanDelete);
    123     gPad->SetLogy();
    124 
    125     if (obj1->GetGeomCam().InheritsFrom("MGeomCamMagic"))
    126     {
    127         TArrayI s0(3);
    128         s0[0] = 6;
    129         s0[1] = 1;
    130         s0[2] = 2;
    131 
    132         TArrayI s1(3);
    133         s1[0] = 3;
    134         s1[1] = 4;
    135         s1[2] = 5;
    136 
    137         TArrayI inner(1);
    138         inner[0] = 0;
    139 
    140         TArrayI outer(1);
    141         outer[0] = 1;
    142 
    143         // Just to get the right (maximum) binning
    144         TH1D *half[4];
    145         half[0] = obj1->ProjectionS(s0, inner, "Sector 6-1-2 Inner");
    146         half[1] = obj1->ProjectionS(s1, inner, "Sector 3-4-5 Inner");
    147         half[2] = obj1->ProjectionS(s0, outer, "Sector 6-1-2 Outer");
    148         half[3] = obj1->ProjectionS(s1, outer, "Sector 3-4-5 Outer");
    149 
    150         for (int i=0; i<4; i++)
    151         {
    152             half[i]->SetLineColor(kRed+i);
    153             half[i]->SetDirectory(0);
    154             half[i]->SetBit(kCanDelete);
    155             half[i]->Draw("same");
    156         }
    157     }
    158 
    159     const Double_t min   = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
    160     const Double_t max   = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
    161     const Double_t integ = obj2->Integral("width")/2.5066283;
    162     const Double_t mean  = obj2->GetMean();
    163     const Double_t rms   = obj2->GetRMS();
    164     const Double_t width = max-min;
    165 
    166     if (rms==0 || width==0)
    167         return;
    168 
    169     const TString dgausformula("([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"
    170                                "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])");
    171 
    172     TF1 *f=0;
    173     switch (fit)
    174     {
    175         // FIXME: MAYBE add function to TH1?
    176     case 0:
    177         f = new TF1("sgaus", "gaus(0)", min, max);
    178         f->SetLineColor(kYellow);
    179         f->SetBit(kCanDelete);
    180         f->SetParNames("Area", "#mu", "#sigma");
    181         f->SetParameters(integ/rms, mean, rms);
    182         f->SetParLimits(0, 0,   integ);
    183         f->SetParLimits(1, min, max);
    184         f->SetParLimits(2, 0,   width/1.5);
    185         obj2->Fit(f, "QLRM");
    186         break;
    187 
    188     case 1:
    189         f = new TF1("dgaus", dgausformula, min, max);
    190         f->SetLineColor(kYellow);
    191         f->SetBit(kCanDelete);
    192         f->SetParNames("A_{tot}", "#mu_{1}", "#sigma_{1}", "A_{2}", "#mu_{2}", "#sigma_{2}");
    193         f->SetParameters(integ,         (min+mean)/2, width/4,
    194                          integ/width/2, (max+mean)/2, width/4);
    195         // The left-sided Gauss
    196         f->SetParLimits(0, integ-1.5,      integ+1.5);
    197         f->SetParLimits(1, min+(width/10), mean);
    198         f->SetParLimits(2, 0,              width/2);
    199         // The right-sided Gauss
    200         f->SetParLimits(3, 0,    integ);
    201         f->SetParLimits(4, mean, max-(width/10));
    202         f->SetParLimits(5, 0,    width/2);
    203         obj2->Fit(f, "QLRM");
    204         break;
    205 
    206     default:
    207         obj2->Fit("gaus", "Q");
    208         obj2->GetFunction("gaus")->SetLineColor(kYellow);
    209         break;
    210     }
    211 }
    212 
    213 void MJPedestal::CamDraw(TCanvas &c, Int_t x, Int_t y, MHCamera &cam1, Int_t fit)
    214 {
    215     c.cd(x);
    216     gPad->SetBorderMode(0);
    217     MHCamera *obj1=(MHCamera*)cam1.DrawCopy("hist");
    218 
    219     obj1->AddNotify(&fPedestalCam);
    220 
    221     c.cd(x+y);
    222     gPad->SetBorderMode(0);
    223     obj1->Draw();
    224 
    225     c.cd(x+2*y);
    226     gPad->SetBorderMode(0);
    227     DrawProjection(obj1, fit);
    228 }
    229122
    230123void MJPedestal::DisplayResult(MParList &plist)
    … …  
    251144    // Create container to display
    252145    //
    253     MHCamera disp0(geomcam, "MPedestalCam;ped", "Pedestals");
    254     MHCamera disp1(geomcam, "MPedestalCam;var", "Sigma Pedestal");
     146    MHCamera disp0(geomcam, "MPedestalCam;ped", "Pedestal");
     147    MHCamera disp1(geomcam, "MPedestalCam;RMS", "Ped. RMS");
    255148
    256149    disp0.SetCamContent(fPedestalCam, 0);
    … …  
    269162    c3.Divide(2,3);
    270163
    271     CamDraw(c3, 1, 2, disp0, 0);
    272     CamDraw(c3, 2, 2, disp1, 1);
     164    CamDraw(c3, 1, 2, disp0, 1);
     165    CamDraw(c3, 2, 2, disp1, 6);
    273166}
    274167
    … …  
    368261    MPedCalcPedRun  pedcalc;
    369262
     263    if (fExtractor)
     264      {
     265        pedcalc.SetWindowSize((Int_t)fExtractor->GetNumHiGainSamples());
     266        pedcalc.SetRange(fExtractor->GetHiGainFirst(), fExtractor->GetHiGainLast());
     267      }
     268    else
     269      {
     270        *fLog << warn << GetDescriptor()
     271            << ": No extractor has been chosen, take default number of FADC samples " << endl;
     272      }
     273
    370274    tlist.AddToList(&geomapl);
    371275    //    tlist.AddToList(&merge);

  • trunk/MagicSoft/Mars/mjobs/MJPedestal.h

    r3810 r3889  
    88#include "MBadPixelsCam.h"
    99#endif
    10 
     10#ifndef MARS_MGCamDisplays
     11#include "MGCamDisplays.h"
     12#endif
    1113
    1214class TCanvas;
    … …  
    1416class MRunIter;
    1517class MHCamera;
    16 
    17 class MJPedestal : public MParContainer
     18class MExtractor;
     19class MJPedestal : public MParContainer, public MGCamDisplays 
    1820{
    1921private:
    20     TString fOutputPath;
    2122
    22     MRunIter *fRuns;
     23  TString fOutputPath;
     24 
     25  MRunIter   *fRuns;
     26  MExtractor *fExtractor;                         // Signal extractor, used to find the nr. of used FADC slices
    2327
    24     MPedestalCam  fPedestalCam;
    25     MBadPixelsCam fBadPixels;
     28  MPedestalCam  fPedestalCam;
     29  MBadPixelsCam fBadPixels;
     30 
     31  Bool_t fDataCheck;                              // Flag if the data check is run on raw data
     32 
     33  Bool_t ReadPedestalCam();
     34  Bool_t WriteResult();
     35 
     36  void   DisplayResult(MParList &plist);
     37 
     38public:
    2639
    27     Bool_t fDataCheck;                              // Flag if the data check is run on raw data
     40  MJPedestal(const char *name=NULL, const char *title=NULL);
    2841
    29     Bool_t ReadPedestalCam();
    30     Bool_t WriteResult();
    31    
    32     void   DrawProjection(MHCamera *obj1, Int_t fit) const;
    33     void   CamDraw(TCanvas &c, Int_t x, Int_t y, MHCamera &cam1, Int_t fit);
    34     void   DisplayResult(MParList &plist);
    35 
    36 public:
    37     MJPedestal(const char *name=NULL, const char *title=NULL);
    38 
    39     void SetOutputPath(const char *path=".");
    40     void SetInput(MRunIter *iter) { fRuns=iter; }
    41 
    42     // Data Check
    43     void SetDataCheck(const Bool_t b=kTRUE) { fDataCheck = b; }
    44 
    45     TString GetOutputFile() const;
    46 
    47     MPedestalCam &GetPedestalCam() { return fPedestalCam; }
    48     const MBadPixelsCam &GetBadPixels() const { return fBadPixels; }
    49 
    50     void SetBadPixels(const MBadPixelsCam &bad) { bad.Copy(fBadPixels); }
    51 
    52     Bool_t ProcessFile();
    53     Bool_t Process();
    54 
    55     ClassDef(MJPedestal, 0) // Tool to create a pedestal file (MPedestalCam)
     42  MPedestalCam &GetPedestalCam()            { return fPedestalCam; }
     43  const MBadPixelsCam &GetBadPixels() const { return fBadPixels;   }
     44 
     45  TString GetOutputFile() const;
     46 
     47  Bool_t Process();
     48  Bool_t ProcessFile();
     49 
     50  void SetBadPixels ( const MBadPixelsCam &bad) { bad.Copy(fBadPixels); }
     51  void SetDataCheck ( const Bool_t b=kTRUE    ) { fDataCheck = b;       }
     52  void SetExtractor (       MExtractor* ext   ) { fExtractor = ext;     }
     53  void SetInput     (       MRunIter *iter    ) { fRuns      = iter;    }
     54  void SetOutputPath( const char *path="."    );
     55 
     56  ClassDef(MJPedestal, 0) // Tool to create a pedestal file (MPedestalCam)
    5657};
    5758

  • trunk/MagicSoft/Mars/mpedestal/MPedCalcPedRun.cc

    r3804 r3889  
    2626!
    2727\* ======================================================================== */
    28 
    2928/////////////////////////////////////////////////////////////////////////////
    3029//
    … …  
    7069// (and moreover asymmetric) and that they are also correlated.)
    7170//
     71//  Call: SetRange(higainfirst, higainlast, logainfirst, logainlast)
     72//  to modify the ranges in which the window is allowed to move.
     73//  Defaults are:
     74//
     75//   fHiGainFirst =  fgHiGainFirst =  0
     76//   fHiGainLast  =  fgHiGainLast  =  14
     77//   fLoGainFirst =  fgLoGainFirst =  0
     78//   fLoGainLast  =  fgLoGainLast  =  14
     79//
     80//  Call: SetWindowSize(windowhigain, windowlogain)
     81//  to modify the sliding window widths. Windows have to be an even number.
     82//  In case of odd numbers, the window will be modified.
     83//
     84//  Defaults are:
     85//
     86//   fHiGainWindowSize = fgHiGainWindowSize = 14
     87//   fLoGainWindowSize = fgLoGainWindowSize = 0
     88//
    7289//
    7390//  Input Containers:
    7491//   MRawEvtData
     92//   MRawRunHeader
     93//   MGeomCam
    7594//
    7695//  Output Containers:
    … …  
    8099/////////////////////////////////////////////////////////////////////////////
    81100#include "MPedCalcPedRun.h"
     101#include "MExtractor.h"
    82102
    83103#include "MParList.h"
    … …  
    93113#include "MPedestalCam.h"
    94114
    95 #include "MExtractedSignalPix.h"
    96 #include "MExtractedSignalCam.h"
    97 
    98115#include "MGeomPix.h"
    99116#include "MGeomCam.h"
    … …  
    105122using namespace std;
    106123
    107 // --------------------------------------------------------------------------
    108 //
    109 // default constructor
     124const Byte_t MPedCalcPedRun::fgHiGainFirst      = 0;
     125const Byte_t MPedCalcPedRun::fgHiGainLast       = 14;
     126const Byte_t MPedCalcPedRun::fgLoGainFirst      = 0;
     127const Byte_t MPedCalcPedRun::fgLoGainLast       = 14;
     128const Byte_t MPedCalcPedRun::fgHiGainWindowSize = 14;
     129const Byte_t MPedCalcPedRun::fgLoGainWindowSize = 0;
     130// --------------------------------------------------------------------------
     131//
     132// Default constructor:
     133//
     134// Sets:
     135// - fWindowSizeHiGain to fgHiGainWindowSize
     136// - fWindowSizeLoGain to fgLoGainWindowSize
     137//
     138// Calls:
     139// - AddToBranchList("fHiGainPixId");
     140// - AddToBranchList("fHiGainFadcSamples");
     141// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast)
     142// - Clear()
    110143//
    111144MPedCalcPedRun::MPedCalcPedRun(const char *name, const char *title)
    112     : fRawEvt(NULL), fPedestals(NULL)
    113 {
    114     fName  = name  ? name  : "MPedCalcPedRun";
    115     fTitle = title ? title : "Task to calculate pedestals from pedestal runs raw data";
    116 
    117     AddToBranchList("fHiGainPixId");
    118     AddToBranchList("fHiGainFadcSamples");
    119 
    120     fNumHiGainSamples = 0;
    121     Clear();
    122 }
    123 
     145    : fWindowSizeHiGain(fgHiGainWindowSize),
     146      fWindowSizeLoGain(fgLoGainWindowSize)
     147{
     148  fName  = name  ? name  : "MPedCalcPedRun";
     149  fTitle = title ? title : "Task to calculate pedestals from pedestal runs raw data";
     150 
     151  AddToBranchList("fHiGainPixId");
     152  AddToBranchList("fHiGainFadcSamples");
     153 
     154  SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
     155  Clear();
     156}
     157
     158// --------------------------------------------------------------------------
     159//
     160// Sets:
     161// - fNumSamplesTot to 0
     162// - fRawEvt to NULL
     163// - fRunHeader to NULL
     164// - fPedestals to NULL
     165//
    124166void MPedCalcPedRun::Clear(const Option_t *o)
    125167{
    … …  
    128170
    129171  fRawEvt    = NULL;
     172  fRunHeader = NULL;
    130173  fPedestals = NULL;
    131174}
    132175
    133 
     176// --------------------------------------------------------------------------
     177//
     178// SetRange:
     179//
     180// Calls:
     181// - MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
     182// - SetWindowSize(fWindowSizeHiGain,fWindowSizeLoGain);
     183//
     184void MPedCalcPedRun::SetRange(Byte_t hifirst, Byte_t hilast, Byte_t lofirst, Byte_t lolast)
     185{
     186
     187  MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
     188
     189  //
     190  // Redo the checks if the window is still inside the ranges
     191  //
     192  SetWindowSize(fWindowSizeHiGain,fWindowSizeLoGain);
     193 
     194}
     195
     196
     197// --------------------------------------------------------------------------
     198//
     199// Checks:
     200// - if a window is odd, subtract one
     201// - if a window is bigger than the one defined by the ranges, set it to the available range
     202// - if a window is smaller than 2, set it to 2
     203//
     204// Sets:
     205// - fNumHiGainSamples to: (Float_t)fWindowSizeHiGain
     206// - fNumLoGainSamples to: (Float_t)fWindowSizeLoGain
     207// - fSqrtHiGainSamples to: TMath::Sqrt(fNumHiGainSamples)
     208// - fSqrtLoGainSamples to: TMath::Sqrt(fNumLoGainSamples) 
     209// 
     210void MPedCalcPedRun::SetWindowSize(Byte_t windowh, Byte_t windowl)
     211{
     212 
     213  fWindowSizeHiGain = windowh & ~1;
     214  fWindowSizeLoGain = windowl & ~1;
     215
     216  if (fWindowSizeHiGain != windowh)
     217    *fLog << warn << GetDescriptor() << ": Hi Gain window size has to be even, set to: "
     218          << int(fWindowSizeHiGain) << " samples " << endl;
     219 
     220  if (fWindowSizeLoGain != windowl)
     221    *fLog << warn << GetDescriptor() << ": Lo Gain window size has to be even, set to: "
     222          << int(fWindowSizeLoGain) << " samples " << endl;
     223   
     224  const Byte_t availhirange = (fHiGainLast-fHiGainFirst+1) & ~1;
     225  const Byte_t availlorange = (fLoGainLast-fLoGainFirst+1) & ~1;
     226
     227  if (fWindowSizeHiGain > availhirange)
     228    {
     229      *fLog << warn << GetDescriptor()
     230            << Form("%s%2i%s%2i%s%2i%s",": Hi Gain window size: ",(int)fWindowSizeHiGain,
     231                    " is bigger than available range: [",(int)fHiGainFirst,",",(int)fHiGainLast,"]") << endl;
     232      *fLog << warn << GetDescriptor()
     233            << ": Will set window size to: " << (int)availhirange << endl;
     234      fWindowSizeHiGain = availhirange;
     235    }
     236 
     237  if (fWindowSizeLoGain > availlorange)
     238    {
     239      *fLog << warn << GetDescriptor()
     240            << Form("%s%2i%s%2i%s%2i%s",": Lo Gain window size: ",(int)fWindowSizeLoGain,
     241                    " is bigger than available range: [",(int)fLoGainFirst,",",(int)fLoGainLast,"]") << endl;
     242      *fLog << warn << GetDescriptor()
     243            << ": Will set window size to: " << (int)availlorange << endl;
     244      fWindowSizeLoGain = availlorange;
     245    }
     246 
     247   
     248  fNumHiGainSamples = (Float_t)fWindowSizeHiGain;
     249  fNumLoGainSamples = (Float_t)fWindowSizeLoGain;
     250 
     251  fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
     252  fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
     253
     254}
     255
     256 
     257   
    134258// --------------------------------------------------------------------------
    135259//
    … …  
    137261//
    138262//  - MRawEvtData
     263//  - MRawRunHeader
     264//  - MGeomCam
    139265//
    140266// The following output containers are also searched and created if
    … …  
    155281    }
    156282 
     283  fRunHeader = (MRawRunHeader*)pList->FindObject(AddSerialNumber("MRawRunHeader"));
     284  if (!fRunHeader)
     285    {
     286      *fLog << err << AddSerialNumber("MRawRunHeader") << " not found... aborting." << endl;
     287      return kFALSE;
     288    }
     289   
    157290  fGeom   =  (MGeomCam*)pList->FindObject("MGeomCam");
    158291  if (!fGeom)
    … …  
    165298  if (!fPedestals)
    166299    return kFALSE;
    167  
     300
    168301  return kTRUE;
    169302}
    … …  
    171304// --------------------------------------------------------------------------
    172305//
    173 // The ReInit searches for the following input containers:
    174 //  - MRawRunHeader
    175 //
    176 // It also initializes the data arrays fSumx and fSumx2
    177 // (only for the first read file)
    178 //
     306// The ReInit searches for:
     307// -  MRawRunHeader::GetNumSamplesHiGain()
     308// -  MRawRunHeader::GetNumSamplesLoGain()
     309//
     310// In case that the variables fHiGainLast and fLoGainLast are smaller than
     311// the even part of the number of samples obtained from the run header, a
     312// warning is given an the range is set back accordingly. A call to: 
     313// - SetRange(fHiGainFirst, fHiGainLast-diff, fLoGainFirst, fLoGainLast) or
     314// - SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast-diff)
     315// is performed in that case. The variable diff means here the difference
     316// between the requested range (fHiGainLast) and the available one. Note that
     317// the functions SetRange() are mostly overloaded and perform more checks,
     318// modifying the ranges again, if necessary.
     319//
    179320Bool_t MPedCalcPedRun::ReInit(MParList *pList)
    180321{
    181     const MRawRunHeader *runheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
    182     if (!runheader)
    183     {
    184         *fLog << warn << dbginf;
    185         *fLog << "Warning - cannot check file type, MRawRunHeader not found." << endl;
    186     }
    187     else
    188       if (runheader->IsMonteCarloRun())
    189         return kTRUE;
    190 
    191     Int_t npixels  = fPedestals->GetSize();
    192     Int_t areas    = fPedestals->GetAverageAreas();
    193     Int_t sectors  = fPedestals->GetAverageSectors();
    194 
    195     if (fSumx.GetSize()==0)
    196     {
    197         fSumx. Set(npixels);
    198         fSumx2.Set(npixels);
    199 
    200         fAreaSumx. Set(areas);
    201         fAreaSumx2.Set(areas);
    202         fAreaValid.Set(areas);
    203 
    204         fSectorSumx. Set(sectors);
    205         fSectorSumx2.Set(sectors);
    206         fSectorValid.Set(sectors);
    207 
    208         fSumx.Reset();
    209         fSumx2.Reset();
    210     }
    211 
    212     // Calculate an even number for the hi gain samples to avoid
    213     // biases due to the fluctuation in pedestal from one slice to
    214     // the other one
    215     fNumHiGainSamples = runheader->GetNumSamplesHiGain() & ~1;
    216 
    217     return kTRUE;
     322 
     323  Int_t lastdesired   = (Int_t)fHiGainLast;
     324  Int_t lastavailable = (Int_t)fRunHeader->GetNumSamplesHiGain()-1;
     325 
     326  if (lastdesired > lastavailable)
     327    {
     328      const Int_t diff = lastdesired - lastavailable;
     329      *fLog << endl;
     330      *fLog << warn << GetDescriptor()
     331            << Form("%s%2i%s%2i%s%2i%s",": Selected Hi Gain FADC Window [",
     332                    (int)fHiGainFirst,",",lastdesired,
     333                    "] ranges out of the available limits: [0,",lastavailable,"].") << endl;
     334      *fLog << GetDescriptor() << ": Will reduce the upper edge to " << (int)(fHiGainLast - diff) << endl;
     335      SetRange(fHiGainFirst, fHiGainLast-diff, fLoGainFirst, fLoGainLast);
     336    }
     337 
     338  lastdesired   = (Int_t)(fLoGainLast);
     339  lastavailable = (Int_t)fRunHeader->GetNumSamplesLoGain()-1;
     340 
     341  if (lastdesired > lastavailable)
     342    {
     343      const Int_t diff = lastdesired - lastavailable;
     344      *fLog << endl;
     345      *fLog << warn << GetDescriptor()
     346            << Form("%s%2i%s%2i%s%2i%s",": Selected Lo Gain FADC Window [",
     347                    (int)fLoGainFirst,",",lastdesired,
     348                    "] ranges out of the available limits: [0,",lastavailable,"].") << endl;
     349      *fLog << GetDescriptor() << ": Will reduce the upper edge to " << (int)(fLoGainLast - diff) << endl;
     350      SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast-diff);
     351    }
     352
     353  Int_t npixels  = fPedestals->GetSize();
     354  Int_t areas    = fPedestals->GetAverageAreas();
     355  Int_t sectors  = fPedestals->GetAverageSectors();
     356 
     357  if (fSumx.GetSize()==0)
     358    {
     359      fSumx. Set(npixels);
     360      fSumx2.Set(npixels);
     361     
     362      fAreaSumx. Set(areas);
     363      fAreaSumx2.Set(areas);
     364      fAreaValid.Set(areas);
     365     
     366      fSectorSumx. Set(sectors);
     367      fSectorSumx2.Set(sectors);
     368      fSectorValid.Set(sectors);
     369     
     370      fSumx.Reset();
     371      fSumx2.Reset();
     372    }
     373 
     374  if (fWindowSizeHiGain == 0 && fWindowSizeLoGain == 0)
     375    {
     376      *fLog << err << GetDescriptor()
     377            << ": Number of extracted Slices is 0, abort ... " << endl;
     378      return kFALSE;
     379    }
     380 
     381
     382  *fLog << endl;
     383  *fLog << inf << GetDescriptor() << ": Taking " << (int)fWindowSizeHiGain
     384        << " HiGain FADC samples starting with slice: " << (int)fHiGainFirst << endl;
     385  *fLog << inf << GetDescriptor() << ": Taking " << (int)fWindowSizeLoGain
     386        << " LoGain FADC samples starting with slice: " << (int)fLoGainFirst << endl;
     387 
     388  return kTRUE;
     389 
    218390}
    219391// --------------------------------------------------------------------------
    … …  
    235407      const UInt_t sector = (*fGeom)[idx].GetSector();     
    236408
    237       Byte_t *ptr = pixel.GetHiGainSamples();
    238       const Byte_t *end = ptr + fNumHiGainSamples;
     409      Byte_t *ptr = pixel.GetHiGainSamples() + fHiGainFirst;
     410      Byte_t *end = ptr + fWindowSizeHiGain;
    239411     
    240412      UInt_t sum = 0;
    … …  
    247419        }
    248420      while (++ptr != end);
     421
     422      if (fWindowSizeLoGain != 0)
     423        {
     424         
     425          ptr = pixel.GetLoGainSamples() + fLoGainFirst;
     426          end = ptr + fWindowSizeLoGain;
     427     
     428          do
     429            {
     430              sum += *ptr;
     431              sqr += *ptr * *ptr;
     432            }
     433          while (++ptr != end);
     434         
     435        }
    249436     
    250437      const Float_t msum = (Float_t)sum;
    … …  
    254441      // If anybody needs them, please contact me!!
    255442      //
    256       //        const Float_t higainped = msum/fNumHiGainSamples;
    257       //        const Float_t higainrms = TMath::Sqrt((msqr-msum*msum/fNumHiGainSamples)/(fNumHiGainSamples-1.));
     443      //        const Float_t higainped = msum/fNumHiGainSlices;
     444      //        const Float_t higainrms = TMath::Sqrt((msqr-msum*msum/fNumHiGainSlices)/(fNumHiGainSlices-1.));
    258445      //        (*fPedestals)[idx].Set(higainped, higainrms);
    259446     
    … …  
    276463 
    277464  fPedestals->SetReadyToSave();
    278   fNumSamplesTot += fNumHiGainSamples;
    279  
     465  fNumSamplesTot += fWindowSizeHiGain + fWindowSizeLoGain;
    280466 
    281467  return kTRUE;
    … …  
    319505      Float_t higainVar = (sum2-sum*sum/nevts)/(nevts-1.);
    320506      // 2. Scale the variance to the number of slices:
    321       higainVar /= (Float_t)fNumHiGainSamples;
     507      higainVar /= (Float_t)(fWindowSizeHiGain+fWindowSizeLoGain);
    322508      // 3. Calculate the RMS from the Variance:
    323509      const Float_t higainrms = TMath::Sqrt(higainVar);
    … …  
    344530      Float_t higainVar = (sum2-sum*sum/aevts)/(aevts-1.);
    345531      // 2. Scale the variance to the number of slices:
    346       higainVar /= (Float_t)fNumHiGainSamples;
     532      higainVar /= fWindowSizeHiGain+fWindowSizeLoGain;
    347533      // 3. Calculate the RMS from the Variance:
    348534      Float_t higainrms = TMath::Sqrt(higainVar);
    … …  
    371557      Float_t higainVar = (sum2-sum*sum/sevts)/(sevts-1.);
    372558      // 2. Scale the variance to the number of slices:
    373       higainVar /= (Float_t)fNumHiGainSamples;
     559      higainVar /= fWindowSizeHiGain+fWindowSizeLoGain;
    374560      // 3. Calculate the RMS from the Variance:
    375561      Float_t higainrms = TMath::Sqrt(higainVar);
Note: See TracChangeset for help on using the changeset viewer.