Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 3888) +++ trunk/MagicSoft/Mars/Changelog (revision 3889) @@ -35,4 +35,10 @@ - add the possibility to set ranges and extraction windows. Default is what has always been + + + * mjobs/MJPedestals.[h,cc] + - add the possibility to set a signal extractor which gives the + extraction ranges to MPedCalcPedRun + - derive from MHCamDisplays Index: trunk/MagicSoft/Mars/mjobs/MGCamDisplays.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MGCamDisplays.cc (revision 3888) +++ trunk/MagicSoft/Mars/mjobs/MGCamDisplays.cc (revision 3889) @@ -121,9 +121,11 @@ // 3: Triple Gauss (for distributions with inner, outer pixels and outliers) // 4: flat (for the probability distributions) +// (1-4:) Moreover, sectors 6,1 and 2 of the camera and sectors 3,4 and 5 are +// drawn separately, for inner and outer pixels. // 5: Fit Inner and Outer pixels separately by a single Gaussian // (only for MAGIC cameras) -// -// Moreover, sectors 6,1 and 2 of the camera and sectors 3,4 and 5 are -// drawn separately, for inner and outer pixels. +// 6: Fit Inner and Outer pixels separately by a single Gaussian and display +// additionally the two camera halfs separately (for MAGIC camera) +// // void MGCamDisplays::DrawProjection(MHCamera *obj, Int_t fit) const @@ -136,5 +138,5 @@ outer[0] = 1; - if (fit==5) + if (fit==5 || fit==6) { @@ -149,39 +151,65 @@ s0[5] = 5; + TArrayI s1(3); + s1[0] = 6; + s1[1] = 1; + s1[2] = 2; + + TArrayI s2(3); + s2[0] = 3; + s2[1] = 4; + s2[2] = 5; + gPad->Clear(); TVirtualPad *pad = gPad; pad->Divide(2,1); - TH1D *half[2]; - half[0] = obj->ProjectionS(s0, inner, "Inner"); - half[1] = obj->ProjectionS(s0, outer, "Outer"); - - half[0]->SetDirectory(NULL); - half[1]->SetDirectory(NULL); + TH1D *inout[2]; + inout[0] = obj->ProjectionS(s0, inner, "Inner"); + inout[1] = obj->ProjectionS(s0, outer, "Outer"); + + inout[0]->SetDirectory(NULL); + inout[1]->SetDirectory(NULL); for (int i=0; i<2; i++) { pad->cd(i+1); - half[i]->SetLineColor(kRed+i); - half[i]->SetBit(kCanDelete); - half[i]->Draw(); - half[i]->Fit("gaus","Q"); + inout[i]->SetLineColor(kRed+i); + inout[i]->SetBit(kCanDelete); + inout[i]->Draw(); + inout[i]->Fit("gaus","Q"); + + if (fit == 6) + { + TH1D *half[2]; + half[0] = obj->ProjectionS(s1, i==0 ? inner : outer , "Sector 6-1-2"); + half[1] = obj->ProjectionS(s2, i==0 ? inner : outer , "Sector 3-4-5"); + + for (int j=0; j<2; j++) + { + half[j]->SetLineColor(kRed+i+j); + half[j]->SetDirectory(0); + half[j]->SetBit(kCanDelete); + half[j]->Draw("same"); + } + } + } gLog << all << obj->GetName() << Form("%s%5.3f%s%3.2f"," Mean: Inner Pixels: ", - half[0]->GetFunction("gaus")->GetParameter(1),"+-", - half[0]->GetFunction("gaus")->GetParError(1)); + inout[0]->GetFunction("gaus")->GetParameter(1),"+-", + inout[0]->GetFunction("gaus")->GetParError(1)); gLog << Form("%s%5.3f%s%3.2f"," Outer Pixels: ", - half[1]->GetFunction("gaus")->GetParameter(1),"+-", - half[1]->GetFunction("gaus")->GetParError(1)) << endl; + inout[1]->GetFunction("gaus")->GetParameter(1),"+-", + inout[1]->GetFunction("gaus")->GetParError(1)) << endl; gLog << all << obj->GetName() << Form("%s%5.3f%s%3.2f"," Sigma: Inner Pixels: ", - half[0]->GetFunction("gaus")->GetParameter(2),"+-", - half[0]->GetFunction("gaus")->GetParError(2)); + inout[0]->GetFunction("gaus")->GetParameter(2),"+-", + inout[0]->GetFunction("gaus")->GetParError(2)); gLog << Form("%s%5.3f%s%3.2f"," Outer Pixels: ", - half[1]->GetFunction("gaus")->GetParameter(2),"+-", - half[1]->GetFunction("gaus")->GetParError(2)) << endl; - + inout[1]->GetFunction("gaus")->GetParameter(2),"+-", + inout[1]->GetFunction("gaus")->GetParError(2)) << endl; + } return; Index: trunk/MagicSoft/Mars/mjobs/MJPedestal.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJPedestal.cc (revision 3888) +++ trunk/MagicSoft/Mars/mjobs/MJPedestal.cc (revision 3889) @@ -45,4 +45,5 @@ #include "MTaskList.h" #include "MEvtLoop.h" +#include "MExtractor.h" #include "MStatusDisplay.h" @@ -62,8 +63,12 @@ using namespace std; - +// -------------------------------------------------------------------------- +// +// Default constructor. +// +// Sets fRuns to 0, fExtractor to NULL, fDataCheck to kFALSE +// MJPedestal::MJPedestal(const char *name, const char *title) - : fRuns(0), fDataCheck(kFALSE) - + : fRuns(0), fExtractor(NULL), fDataCheck(kFALSE) { fName = name ? name : "MJPedestal"; @@ -115,116 +120,4 @@ } -void MJPedestal::DrawProjection(MHCamera *obj1, Int_t fit) const -{ - TH1D *obj2 = (TH1D*)obj1->Projection(obj1->GetName()); - obj2->SetDirectory(0); - obj2->Draw(); - obj2->SetBit(kCanDelete); - gPad->SetLogy(); - - if (obj1->GetGeomCam().InheritsFrom("MGeomCamMagic")) - { - TArrayI s0(3); - s0[0] = 6; - s0[1] = 1; - s0[2] = 2; - - TArrayI s1(3); - s1[0] = 3; - s1[1] = 4; - s1[2] = 5; - - TArrayI inner(1); - inner[0] = 0; - - TArrayI outer(1); - outer[0] = 1; - - // Just to get the right (maximum) binning - TH1D *half[4]; - half[0] = obj1->ProjectionS(s0, inner, "Sector 6-1-2 Inner"); - half[1] = obj1->ProjectionS(s1, inner, "Sector 3-4-5 Inner"); - half[2] = obj1->ProjectionS(s0, outer, "Sector 6-1-2 Outer"); - half[3] = obj1->ProjectionS(s1, outer, "Sector 3-4-5 Outer"); - - for (int i=0; i<4; i++) - { - half[i]->SetLineColor(kRed+i); - half[i]->SetDirectory(0); - half[i]->SetBit(kCanDelete); - half[i]->Draw("same"); - } - } - - const Double_t min = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst()); - const Double_t max = obj2->GetBinCenter(obj2->GetXaxis()->GetLast()); - const Double_t integ = obj2->Integral("width")/2.5066283; - const Double_t mean = obj2->GetMean(); - const Double_t rms = obj2->GetRMS(); - const Double_t width = max-min; - - if (rms==0 || width==0) - return; - - const TString dgausformula("([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])" - "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])"); - - TF1 *f=0; - switch (fit) - { - // FIXME: MAYBE add function to TH1? - case 0: - f = new TF1("sgaus", "gaus(0)", min, max); - f->SetLineColor(kYellow); - f->SetBit(kCanDelete); - f->SetParNames("Area", "#mu", "#sigma"); - f->SetParameters(integ/rms, mean, rms); - f->SetParLimits(0, 0, integ); - f->SetParLimits(1, min, max); - f->SetParLimits(2, 0, width/1.5); - obj2->Fit(f, "QLRM"); - break; - - case 1: - f = new TF1("dgaus", dgausformula, min, max); - f->SetLineColor(kYellow); - f->SetBit(kCanDelete); - f->SetParNames("A_{tot}", "#mu_{1}", "#sigma_{1}", "A_{2}", "#mu_{2}", "#sigma_{2}"); - f->SetParameters(integ, (min+mean)/2, width/4, - integ/width/2, (max+mean)/2, width/4); - // The left-sided Gauss - f->SetParLimits(0, integ-1.5, integ+1.5); - f->SetParLimits(1, min+(width/10), mean); - f->SetParLimits(2, 0, width/2); - // The right-sided Gauss - f->SetParLimits(3, 0, integ); - f->SetParLimits(4, mean, max-(width/10)); - f->SetParLimits(5, 0, width/2); - obj2->Fit(f, "QLRM"); - break; - - default: - obj2->Fit("gaus", "Q"); - obj2->GetFunction("gaus")->SetLineColor(kYellow); - break; - } -} - -void MJPedestal::CamDraw(TCanvas &c, Int_t x, Int_t y, MHCamera &cam1, Int_t fit) -{ - c.cd(x); - gPad->SetBorderMode(0); - MHCamera *obj1=(MHCamera*)cam1.DrawCopy("hist"); - - obj1->AddNotify(&fPedestalCam); - - c.cd(x+y); - gPad->SetBorderMode(0); - obj1->Draw(); - - c.cd(x+2*y); - gPad->SetBorderMode(0); - DrawProjection(obj1, fit); -} void MJPedestal::DisplayResult(MParList &plist) @@ -251,6 +144,6 @@ // Create container to display // - MHCamera disp0(geomcam, "MPedestalCam;ped", "Pedestals"); - MHCamera disp1(geomcam, "MPedestalCam;var", "Sigma Pedestal"); + MHCamera disp0(geomcam, "MPedestalCam;ped", "Pedestal"); + MHCamera disp1(geomcam, "MPedestalCam;RMS", "Ped. RMS"); disp0.SetCamContent(fPedestalCam, 0); @@ -269,6 +162,6 @@ c3.Divide(2,3); - CamDraw(c3, 1, 2, disp0, 0); - CamDraw(c3, 2, 2, disp1, 1); + CamDraw(c3, 1, 2, disp0, 1); + CamDraw(c3, 2, 2, disp1, 6); } @@ -368,4 +261,15 @@ MPedCalcPedRun pedcalc; + if (fExtractor) + { + pedcalc.SetWindowSize((Int_t)fExtractor->GetNumHiGainSamples()); + pedcalc.SetRange(fExtractor->GetHiGainFirst(), fExtractor->GetHiGainLast()); + } + else + { + *fLog << warn << GetDescriptor() + << ": No extractor has been chosen, take default number of FADC samples " << endl; + } + tlist.AddToList(&geomapl); // tlist.AddToList(&merge); Index: trunk/MagicSoft/Mars/mjobs/MJPedestal.h =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJPedestal.h (revision 3888) +++ trunk/MagicSoft/Mars/mjobs/MJPedestal.h (revision 3889) @@ -8,5 +8,7 @@ #include "MBadPixelsCam.h" #endif - +#ifndef MARS_MGCamDisplays +#include "MGCamDisplays.h" +#endif class TCanvas; @@ -14,44 +16,43 @@ class MRunIter; class MHCamera; - -class MJPedestal : public MParContainer +class MExtractor; +class MJPedestal : public MParContainer, public MGCamDisplays { private: - TString fOutputPath; - MRunIter *fRuns; + TString fOutputPath; + + MRunIter *fRuns; + MExtractor *fExtractor; // Signal extractor, used to find the nr. of used FADC slices - MPedestalCam fPedestalCam; - MBadPixelsCam fBadPixels; + MPedestalCam fPedestalCam; + MBadPixelsCam fBadPixels; + + Bool_t fDataCheck; // Flag if the data check is run on raw data + + Bool_t ReadPedestalCam(); + Bool_t WriteResult(); + + void DisplayResult(MParList &plist); + +public: - Bool_t fDataCheck; // Flag if the data check is run on raw data + MJPedestal(const char *name=NULL, const char *title=NULL); - Bool_t ReadPedestalCam(); - Bool_t WriteResult(); - - void DrawProjection(MHCamera *obj1, Int_t fit) const; - void CamDraw(TCanvas &c, Int_t x, Int_t y, MHCamera &cam1, Int_t fit); - void DisplayResult(MParList &plist); - -public: - MJPedestal(const char *name=NULL, const char *title=NULL); - - void SetOutputPath(const char *path="."); - void SetInput(MRunIter *iter) { fRuns=iter; } - - // Data Check - void SetDataCheck(const Bool_t b=kTRUE) { fDataCheck = b; } - - TString GetOutputFile() const; - - MPedestalCam &GetPedestalCam() { return fPedestalCam; } - const MBadPixelsCam &GetBadPixels() const { return fBadPixels; } - - void SetBadPixels(const MBadPixelsCam &bad) { bad.Copy(fBadPixels); } - - Bool_t ProcessFile(); - Bool_t Process(); - - ClassDef(MJPedestal, 0) // Tool to create a pedestal file (MPedestalCam) + MPedestalCam &GetPedestalCam() { return fPedestalCam; } + const MBadPixelsCam &GetBadPixels() const { return fBadPixels; } + + TString GetOutputFile() const; + + Bool_t Process(); + Bool_t ProcessFile(); + + void SetBadPixels ( const MBadPixelsCam &bad) { bad.Copy(fBadPixels); } + void SetDataCheck ( const Bool_t b=kTRUE ) { fDataCheck = b; } + void SetExtractor ( MExtractor* ext ) { fExtractor = ext; } + void SetInput ( MRunIter *iter ) { fRuns = iter; } + void SetOutputPath( const char *path="." ); + + ClassDef(MJPedestal, 0) // Tool to create a pedestal file (MPedestalCam) }; Index: trunk/MagicSoft/Mars/mpedestal/MPedCalcPedRun.cc =================================================================== --- trunk/MagicSoft/Mars/mpedestal/MPedCalcPedRun.cc (revision 3888) +++ trunk/MagicSoft/Mars/mpedestal/MPedCalcPedRun.cc (revision 3889) @@ -26,5 +26,4 @@ ! \* ======================================================================== */ - ///////////////////////////////////////////////////////////////////////////// // @@ -70,7 +69,27 @@ // (and moreover asymmetric) and that they are also correlated.) // +// Call: SetRange(higainfirst, higainlast, logainfirst, logainlast) +// to modify the ranges in which the window is allowed to move. +// Defaults are: +// +// fHiGainFirst = fgHiGainFirst = 0 +// fHiGainLast = fgHiGainLast = 14 +// fLoGainFirst = fgLoGainFirst = 0 +// fLoGainLast = fgLoGainLast = 14 +// +// Call: SetWindowSize(windowhigain, windowlogain) +// to modify the sliding window widths. Windows have to be an even number. +// In case of odd numbers, the window will be modified. +// +// Defaults are: +// +// fHiGainWindowSize = fgHiGainWindowSize = 14 +// fLoGainWindowSize = fgLoGainWindowSize = 0 +// // // Input Containers: // MRawEvtData +// MRawRunHeader +// MGeomCam // // Output Containers: @@ -80,4 +99,5 @@ ///////////////////////////////////////////////////////////////////////////// #include "MPedCalcPedRun.h" +#include "MExtractor.h" #include "MParList.h" @@ -93,7 +113,4 @@ #include "MPedestalCam.h" -#include "MExtractedSignalPix.h" -#include "MExtractedSignalCam.h" - #include "MGeomPix.h" #include "MGeomCam.h" @@ -105,21 +122,46 @@ using namespace std; -// -------------------------------------------------------------------------- -// -// default constructor +const Byte_t MPedCalcPedRun::fgHiGainFirst = 0; +const Byte_t MPedCalcPedRun::fgHiGainLast = 14; +const Byte_t MPedCalcPedRun::fgLoGainFirst = 0; +const Byte_t MPedCalcPedRun::fgLoGainLast = 14; +const Byte_t MPedCalcPedRun::fgHiGainWindowSize = 14; +const Byte_t MPedCalcPedRun::fgLoGainWindowSize = 0; +// -------------------------------------------------------------------------- +// +// Default constructor: +// +// Sets: +// - fWindowSizeHiGain to fgHiGainWindowSize +// - fWindowSizeLoGain to fgLoGainWindowSize +// +// Calls: +// - AddToBranchList("fHiGainPixId"); +// - AddToBranchList("fHiGainFadcSamples"); +// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast) +// - Clear() // MPedCalcPedRun::MPedCalcPedRun(const char *name, const char *title) - : fRawEvt(NULL), fPedestals(NULL) -{ - fName = name ? name : "MPedCalcPedRun"; - fTitle = title ? title : "Task to calculate pedestals from pedestal runs raw data"; - - AddToBranchList("fHiGainPixId"); - AddToBranchList("fHiGainFadcSamples"); - - fNumHiGainSamples = 0; - Clear(); -} - + : fWindowSizeHiGain(fgHiGainWindowSize), + fWindowSizeLoGain(fgLoGainWindowSize) +{ + fName = name ? name : "MPedCalcPedRun"; + fTitle = title ? title : "Task to calculate pedestals from pedestal runs raw data"; + + AddToBranchList("fHiGainPixId"); + AddToBranchList("fHiGainFadcSamples"); + + SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast); + Clear(); +} + +// -------------------------------------------------------------------------- +// +// Sets: +// - fNumSamplesTot to 0 +// - fRawEvt to NULL +// - fRunHeader to NULL +// - fPedestals to NULL +// void MPedCalcPedRun::Clear(const Option_t *o) { @@ -128,8 +170,90 @@ fRawEvt = NULL; + fRunHeader = NULL; fPedestals = NULL; } - +// -------------------------------------------------------------------------- +// +// SetRange: +// +// Calls: +// - MExtractor::SetRange(hifirst,hilast,lofirst,lolast); +// - SetWindowSize(fWindowSizeHiGain,fWindowSizeLoGain); +// +void MPedCalcPedRun::SetRange(Byte_t hifirst, Byte_t hilast, Byte_t lofirst, Byte_t lolast) +{ + + MExtractor::SetRange(hifirst,hilast,lofirst,lolast); + + // + // Redo the checks if the window is still inside the ranges + // + SetWindowSize(fWindowSizeHiGain,fWindowSizeLoGain); + +} + + +// -------------------------------------------------------------------------- +// +// Checks: +// - if a window is odd, subtract one +// - if a window is bigger than the one defined by the ranges, set it to the available range +// - if a window is smaller than 2, set it to 2 +// +// Sets: +// - fNumHiGainSamples to: (Float_t)fWindowSizeHiGain +// - fNumLoGainSamples to: (Float_t)fWindowSizeLoGain +// - fSqrtHiGainSamples to: TMath::Sqrt(fNumHiGainSamples) +// - fSqrtLoGainSamples to: TMath::Sqrt(fNumLoGainSamples) +// +void MPedCalcPedRun::SetWindowSize(Byte_t windowh, Byte_t windowl) +{ + + fWindowSizeHiGain = windowh & ~1; + fWindowSizeLoGain = windowl & ~1; + + if (fWindowSizeHiGain != windowh) + *fLog << warn << GetDescriptor() << ": Hi Gain window size has to be even, set to: " + << int(fWindowSizeHiGain) << " samples " << endl; + + if (fWindowSizeLoGain != windowl) + *fLog << warn << GetDescriptor() << ": Lo Gain window size has to be even, set to: " + << int(fWindowSizeLoGain) << " samples " << endl; + + const Byte_t availhirange = (fHiGainLast-fHiGainFirst+1) & ~1; + const Byte_t availlorange = (fLoGainLast-fLoGainFirst+1) & ~1; + + if (fWindowSizeHiGain > availhirange) + { + *fLog << warn << GetDescriptor() + << Form("%s%2i%s%2i%s%2i%s",": Hi Gain window size: ",(int)fWindowSizeHiGain, + " is bigger than available range: [",(int)fHiGainFirst,",",(int)fHiGainLast,"]") << endl; + *fLog << warn << GetDescriptor() + << ": Will set window size to: " << (int)availhirange << endl; + fWindowSizeHiGain = availhirange; + } + + if (fWindowSizeLoGain > availlorange) + { + *fLog << warn << GetDescriptor() + << Form("%s%2i%s%2i%s%2i%s",": Lo Gain window size: ",(int)fWindowSizeLoGain, + " is bigger than available range: [",(int)fLoGainFirst,",",(int)fLoGainLast,"]") << endl; + *fLog << warn << GetDescriptor() + << ": Will set window size to: " << (int)availlorange << endl; + fWindowSizeLoGain = availlorange; + } + + + fNumHiGainSamples = (Float_t)fWindowSizeHiGain; + fNumLoGainSamples = (Float_t)fWindowSizeLoGain; + + fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples); + fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples); + +} + + + // -------------------------------------------------------------------------- // @@ -137,4 +261,6 @@ // // - MRawEvtData +// - MRawRunHeader +// - MGeomCam // // The following output containers are also searched and created if @@ -155,4 +281,11 @@ } + fRunHeader = (MRawRunHeader*)pList->FindObject(AddSerialNumber("MRawRunHeader")); + if (!fRunHeader) + { + *fLog << err << AddSerialNumber("MRawRunHeader") << " not found... aborting." << endl; + return kFALSE; + } + fGeom = (MGeomCam*)pList->FindObject("MGeomCam"); if (!fGeom) @@ -165,5 +298,5 @@ if (!fPedestals) return kFALSE; - + return kTRUE; } @@ -171,49 +304,88 @@ // -------------------------------------------------------------------------- // -// The ReInit searches for the following input containers: -// - MRawRunHeader -// -// It also initializes the data arrays fSumx and fSumx2 -// (only for the first read file) -// +// The ReInit searches for: +// - MRawRunHeader::GetNumSamplesHiGain() +// - MRawRunHeader::GetNumSamplesLoGain() +// +// In case that the variables fHiGainLast and fLoGainLast are smaller than +// the even part of the number of samples obtained from the run header, a +// warning is given an the range is set back accordingly. A call to: +// - SetRange(fHiGainFirst, fHiGainLast-diff, fLoGainFirst, fLoGainLast) or +// - SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast-diff) +// is performed in that case. The variable diff means here the difference +// between the requested range (fHiGainLast) and the available one. Note that +// the functions SetRange() are mostly overloaded and perform more checks, +// modifying the ranges again, if necessary. +// Bool_t MPedCalcPedRun::ReInit(MParList *pList) { - const MRawRunHeader *runheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); - if (!runheader) - { - *fLog << warn << dbginf; - *fLog << "Warning - cannot check file type, MRawRunHeader not found." << endl; - } - else - if (runheader->IsMonteCarloRun()) - return kTRUE; - - Int_t npixels = fPedestals->GetSize(); - Int_t areas = fPedestals->GetAverageAreas(); - Int_t sectors = fPedestals->GetAverageSectors(); - - if (fSumx.GetSize()==0) - { - fSumx. Set(npixels); - fSumx2.Set(npixels); - - fAreaSumx. Set(areas); - fAreaSumx2.Set(areas); - fAreaValid.Set(areas); - - fSectorSumx. Set(sectors); - fSectorSumx2.Set(sectors); - fSectorValid.Set(sectors); - - fSumx.Reset(); - fSumx2.Reset(); - } - - // Calculate an even number for the hi gain samples to avoid - // biases due to the fluctuation in pedestal from one slice to - // the other one - fNumHiGainSamples = runheader->GetNumSamplesHiGain() & ~1; - - return kTRUE; + + Int_t lastdesired = (Int_t)fHiGainLast; + Int_t lastavailable = (Int_t)fRunHeader->GetNumSamplesHiGain()-1; + + if (lastdesired > lastavailable) + { + const Int_t diff = lastdesired - lastavailable; + *fLog << endl; + *fLog << warn << GetDescriptor() + << Form("%s%2i%s%2i%s%2i%s",": Selected Hi Gain FADC Window [", + (int)fHiGainFirst,",",lastdesired, + "] ranges out of the available limits: [0,",lastavailable,"].") << endl; + *fLog << GetDescriptor() << ": Will reduce the upper edge to " << (int)(fHiGainLast - diff) << endl; + SetRange(fHiGainFirst, fHiGainLast-diff, fLoGainFirst, fLoGainLast); + } + + lastdesired = (Int_t)(fLoGainLast); + lastavailable = (Int_t)fRunHeader->GetNumSamplesLoGain()-1; + + if (lastdesired > lastavailable) + { + const Int_t diff = lastdesired - lastavailable; + *fLog << endl; + *fLog << warn << GetDescriptor() + << Form("%s%2i%s%2i%s%2i%s",": Selected Lo Gain FADC Window [", + (int)fLoGainFirst,",",lastdesired, + "] ranges out of the available limits: [0,",lastavailable,"].") << endl; + *fLog << GetDescriptor() << ": Will reduce the upper edge to " << (int)(fLoGainLast - diff) << endl; + SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast-diff); + } + + Int_t npixels = fPedestals->GetSize(); + Int_t areas = fPedestals->GetAverageAreas(); + Int_t sectors = fPedestals->GetAverageSectors(); + + if (fSumx.GetSize()==0) + { + fSumx. Set(npixels); + fSumx2.Set(npixels); + + fAreaSumx. Set(areas); + fAreaSumx2.Set(areas); + fAreaValid.Set(areas); + + fSectorSumx. Set(sectors); + fSectorSumx2.Set(sectors); + fSectorValid.Set(sectors); + + fSumx.Reset(); + fSumx2.Reset(); + } + + if (fWindowSizeHiGain == 0 && fWindowSizeLoGain == 0) + { + *fLog << err << GetDescriptor() + << ": Number of extracted Slices is 0, abort ... " << endl; + return kFALSE; + } + + + *fLog << endl; + *fLog << inf << GetDescriptor() << ": Taking " << (int)fWindowSizeHiGain + << " HiGain FADC samples starting with slice: " << (int)fHiGainFirst << endl; + *fLog << inf << GetDescriptor() << ": Taking " << (int)fWindowSizeLoGain + << " LoGain FADC samples starting with slice: " << (int)fLoGainFirst << endl; + + return kTRUE; + } // -------------------------------------------------------------------------- @@ -235,6 +407,6 @@ const UInt_t sector = (*fGeom)[idx].GetSector(); - Byte_t *ptr = pixel.GetHiGainSamples(); - const Byte_t *end = ptr + fNumHiGainSamples; + Byte_t *ptr = pixel.GetHiGainSamples() + fHiGainFirst; + Byte_t *end = ptr + fWindowSizeHiGain; UInt_t sum = 0; @@ -247,4 +419,19 @@ } while (++ptr != end); + + if (fWindowSizeLoGain != 0) + { + + ptr = pixel.GetLoGainSamples() + fLoGainFirst; + end = ptr + fWindowSizeLoGain; + + do + { + sum += *ptr; + sqr += *ptr * *ptr; + } + while (++ptr != end); + + } const Float_t msum = (Float_t)sum; @@ -254,6 +441,6 @@ // If anybody needs them, please contact me!! // - // const Float_t higainped = msum/fNumHiGainSamples; - // const Float_t higainrms = TMath::Sqrt((msqr-msum*msum/fNumHiGainSamples)/(fNumHiGainSamples-1.)); + // const Float_t higainped = msum/fNumHiGainSlices; + // const Float_t higainrms = TMath::Sqrt((msqr-msum*msum/fNumHiGainSlices)/(fNumHiGainSlices-1.)); // (*fPedestals)[idx].Set(higainped, higainrms); @@ -276,6 +463,5 @@ fPedestals->SetReadyToSave(); - fNumSamplesTot += fNumHiGainSamples; - + fNumSamplesTot += fWindowSizeHiGain + fWindowSizeLoGain; return kTRUE; @@ -319,5 +505,5 @@ Float_t higainVar = (sum2-sum*sum/nevts)/(nevts-1.); // 2. Scale the variance to the number of slices: - higainVar /= (Float_t)fNumHiGainSamples; + higainVar /= (Float_t)(fWindowSizeHiGain+fWindowSizeLoGain); // 3. Calculate the RMS from the Variance: const Float_t higainrms = TMath::Sqrt(higainVar); @@ -344,5 +530,5 @@ Float_t higainVar = (sum2-sum*sum/aevts)/(aevts-1.); // 2. Scale the variance to the number of slices: - higainVar /= (Float_t)fNumHiGainSamples; + higainVar /= fWindowSizeHiGain+fWindowSizeLoGain; // 3. Calculate the RMS from the Variance: Float_t higainrms = TMath::Sqrt(higainVar); @@ -371,5 +557,5 @@ Float_t higainVar = (sum2-sum*sum/sevts)/(sevts-1.); // 2. Scale the variance to the number of slices: - higainVar /= (Float_t)fNumHiGainSamples; + higainVar /= fWindowSizeHiGain+fWindowSizeLoGain; // 3. Calculate the RMS from the Variance: Float_t higainrms = TMath::Sqrt(higainVar);