Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 4751) +++ trunk/MagicSoft/Mars/Changelog (revision 4752) @@ -19,4 +19,44 @@ -*-*- END OF LINE -*-*- + 2004/08/27: Thomas Bretz + + * Makefile: + - added comments how to link statically + + * callisto.cc: + - fixed some output + + * mbadpixels/Makefile: + - added a comment + + * mbase/BaseLinkDef.h, mbase/Makefile: + - added MArrayI + + * mbase/MArrayI.[h,cc]: + - added + + * mbase/MArrayD.cc: + - fixed some comments + + * mcalib/MCalibrateData.[h,cc]: + - unified CalibratePedestal and CalibrateData. Calling GetConvFactor twice + took a lot of time. + + * mjobs/MJCalibrateSignal.cc, mjobs/MJPedestal.cc: + - added two empty lines to output if finished + + * mpedestal/MPedPhotCam.cc: + - use faster MArrays in ReCalc + - accelerated GetPixelContent + + * msignal/MExtractTimeFastSpline.cc: + - accelerated a bit by defining + Float_t higainklo = fHiGainSecondDeriv[klo]; + Float_t higainkhi = fHiGainSecondDeriv[khi]; + instead of accesing the arrays many times inside the loops. + Somebody should do the same for logain. + + + 2004/08/26: Wolfgang Wittek Index: trunk/MagicSoft/Mars/Makefile =================================================================== --- trunk/MagicSoft/Mars/Makefile (revision 4751) +++ trunk/MagicSoft/Mars/Makefile (revision 4752) @@ -106,4 +106,9 @@ @echo " Linking $@ ..." $(CXX) $(CXXFLAGS) $(ROOTGLIBS) $(SOLIB) $@.o $(MARS_LIB) -o $@ + +# Use this to link the programs statically - for gprof +#$(PROGRAMS): $(OBJS) $(HEADERS) $(PROGRAMS:=.o) +# @echo " Linking $@ ..." +# $(CXX) $(CXXFLAGS) $(ROOTGLIBS) $(OBJS) $(SUBDIRS:=/*.o) $@.o $(MARS_LIB) -o $@ else $(DYLIB): $(LIBRARIES) $(OBJS) $(HEADERS) Index: trunk/MagicSoft/Mars/callisto.cc =================================================================== --- trunk/MagicSoft/Mars/callisto.cc (revision 4751) +++ trunk/MagicSoft/Mars/callisto.cc (revision 4752) @@ -325,6 +325,4 @@ return 1; } - - gLog << flush << all << endl << endl; } @@ -380,6 +378,4 @@ return 1; } - - gLog << flush << all << endl << endl; } Index: trunk/MagicSoft/Mars/mbase/BaseLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mbase/BaseLinkDef.h (revision 4751) +++ trunk/MagicSoft/Mars/mbase/BaseLinkDef.h (revision 4752) @@ -60,4 +60,5 @@ #pragma link C++ class MArrayS; #pragma link C++ class MArrayD; +#pragma link C++ class MArrayI; #pragma link C++ class MTime+; Index: trunk/MagicSoft/Mars/mbase/MArrayD.cc =================================================================== --- trunk/MagicSoft/Mars/mbase/MArrayD.cc (revision 4751) +++ trunk/MagicSoft/Mars/mbase/MArrayD.cc (revision 4752) @@ -26,7 +26,7 @@ ////////////////////////////////////////////////////////////////////////////// // -// MArrayS +// MArrayD // -// Array of UShort_t. It is almost the same than TArrayD but derives from +// Array of Double_t. It is almost the same than TArrayD but derives from // TObject // Index: trunk/MagicSoft/Mars/mbase/MArrayI.cc =================================================================== --- trunk/MagicSoft/Mars/mbase/MArrayI.cc (revision 4752) +++ trunk/MagicSoft/Mars/mbase/MArrayI.cc (revision 4752) @@ -0,0 +1,39 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Thomas Bretz 12/2000 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + + +////////////////////////////////////////////////////////////////////////////// +// +// MArrayI +// +// Array of Int_t. It is almost the same than TArrayI but derives from +// TObject +// +// Another advantage is: operator[] has no range check! +// +////////////////////////////////////////////////////////////////////////////// +#include "MArrayI.h" + +ClassImp(MArrayI); + Index: trunk/MagicSoft/Mars/mbase/MArrayI.h =================================================================== --- trunk/MagicSoft/Mars/mbase/MArrayI.h (revision 4752) +++ trunk/MagicSoft/Mars/mbase/MArrayI.h (revision 4752) @@ -0,0 +1,157 @@ +#ifndef MARS_MArrayI +#define MARS_MArrayI + +#ifndef MARS_MArray +#include "MArray.h" +#endif + +#include + +class MArrayI : public MArray +{ +private: + Int_t *fArray; //[fN] Array of fN chars + +public: + + MArrayI() + { + fN = 0; + fArray = NULL; + } + + MArrayI(UInt_t n) + { + fN = 0; + fArray = NULL; + Set(n); + } + + MArrayI(UInt_t n, Int_t *array) + { + // Create TArrayC object and initialize it with values of array. + fN = 0; + fArray = NULL; + Set(n, array); + } + + MArrayI(const MArrayI &array) + { + // Copy constructor. + fArray = NULL; + Set(array.fN, array.fArray); + } + + UInt_t GetSize() const + { + return fN; + } + + MArrayI &operator=(const MArrayI &rhs) + { + // TArrayC assignment operator. + if (this != &rhs) + Set(rhs.fN, rhs.fArray); + return *this; + } + + virtual ~MArrayI() + { + // Delete TArrayC object. + delete [] fArray; + fArray = NULL; + } + + void Adopt(UInt_t n, Int_t *array) + { + // Adopt array arr into TArrayC, i.e. don't copy arr but use it directly + // in TArrayC. User may not delete arr, TArrayC dtor will do it. + if (fArray) + delete [] fArray; + + fN = n; + fArray = array; + } + + void AddAt(Int_t c, UInt_t i) + { + // Add char c at position i. Check for out of bounds. + fArray[i] = c; + } + + Int_t At(UInt_t i) + { + return fArray[i]; + } + + Int_t *GetArray() const + { + return fArray; + } + + void Reset() + { + memset(fArray, 0, fN*sizeof(Int_t)); + } + + void Set(UInt_t n) + { + // Set size of this array to n chars. + // A new array is created, the old contents copied to the new array, + // then the old array is deleted. + + if (n==fN) + return; + + Int_t *temp = fArray; + if (n == 0) + fArray = NULL; + else + { + fArray = new Int_t[n]; + if (n < fN) + memcpy(fArray, temp, n*sizeof(Int_t)); + else + { + memcpy(fArray, temp, fN*sizeof(Int_t)); + memset(&fArray[fN], 0, (n-fN)*sizeof(Int_t)); + } + } + + if (fN) + delete [] temp; + + fN = n; + } + + void Set(UInt_t n, Int_t *array) + { + // Set size of this array to n chars and set the contents. + if (!array) + return; + + if (fArray && fN != n) + { + delete [] fArray; + fArray = 0; + } + fN = n; + + if (fN == 0) + return; + + if (!fArray) + fArray = new Int_t[fN]; + + memcpy(fArray, array, n*sizeof(Int_t)); + } + + Int_t &operator[](UInt_t i) + { + return fArray[i]; + } + + ClassDef(MArrayI, 1) //Array of Int_t +}; + +#endif Index: trunk/MagicSoft/Mars/mbase/Makefile =================================================================== --- trunk/MagicSoft/Mars/mbase/Makefile (revision 4751) +++ trunk/MagicSoft/Mars/mbase/Makefile (revision 4752) @@ -52,4 +52,5 @@ MArrayS.cc \ MArrayD.cc \ + MArrayI.cc \ MTime.cc \ MClone.cc \ Index: trunk/MagicSoft/Mars/mcalib/MCalibrateData.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrateData.cc (revision 4751) +++ trunk/MagicSoft/Mars/mcalib/MCalibrateData.cc (revision 4752) @@ -265,4 +265,11 @@ } + // Sizes might have changed + if (fPedestalFlag && (Int_t)fPedestal->GetSize() != fSignals->GetSize()) + { + *fLog << err << "Size mismatch of MPedestalCam and MCalibrationCam... abort." << endl; + return kFALSE; + } + if(fCalibrationMode == kBlindPixel && !fQEs->IsBlindPixelMethodValid()) { @@ -315,6 +322,5 @@ if (TestPedestalFlag(kRun)) - if (!CalibratePedestal()) - return kFALSE; + Calibrate(kFALSE, kTRUE); return kTRUE; @@ -323,63 +329,9 @@ // -------------------------------------------------------------------------- // -// Calibrate the Pedestal values -// -Bool_t MCalibrateData::CalibratePedestal() -{ - // - // Fill MPedPhot container using the informations from - // MPedestalCam, MExtractedSignalCam and MCalibrationCam - // - const Float_t slices = fSignals->GetNumUsedHiGainFADCSlices(); - const Float_t sqrtslices = TMath::Sqrt(slices); - - // is pixid equal to pixidx ? - if ((Int_t)fPedestal->GetSize() != fSignals->GetSize()) - { - *fLog << err << "Sizes of MPedestalCam and MCalibrationCam are different" << endl; - return kFALSE; - } - - const Int_t n = fPedestal->GetSize(); - for (Int_t pixid=0; pixidSetReadyToSave(); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// // Get conversion factor and its error from MCalibrationCam // Bool_t MCalibrateData::GetConversionFactor(UInt_t pixidx, Float_t &hiloconv, Float_t &hiloconverr, - Float_t &calibConv, Float_t &calibConvVar, Float_t &calibFFactor) + Float_t &calibConv, Float_t &calibConvVar, + Float_t &calibFFactor) const { // @@ -486,12 +438,17 @@ } -void MCalibrateData::CalibrateData() -{ - const UInt_t npix = fSignals->GetSize(); +void MCalibrateData::Calibrate(Bool_t data, Bool_t pedestal) const +{ + if (!data && !pedestal) + return; + + const UInt_t npix = fSignals->GetSize(); + const Float_t slices = fSignals->GetNumUsedHiGainFADCSlices(); + const Float_t sqrtslices = TMath::Sqrt(slices); Float_t hiloconv; Float_t hiloconverr; - Float_t calibrationConversionFactor; - Float_t calibrationConversionFactorErr; + Float_t calibConv; + Float_t calibConvErr; Float_t calibFFactor; @@ -499,47 +456,72 @@ { if (!GetConversionFactor(pixidx, hiloconv, hiloconverr, - calibrationConversionFactor, calibrationConversionFactorErr, calibFFactor)) + calibConv, calibConvErr, calibFFactor)) continue; - const MExtractedSignalPix &sig = (*fSignals)[pixidx]; - - Float_t signal = 0; - Float_t signalErr = 0.; - - if (sig.IsLoGainUsed()) - { - signal = sig.GetExtractedSignalLoGain()*hiloconv; - signalErr = signal*hiloconverr; - } - else - { - if (sig.GetExtractedSignalHiGain() <= 9999.) - signal = sig.GetExtractedSignalHiGain(); - } - - const Float_t nphot = signal*calibrationConversionFactor; - const Float_t nphotErr = calibFFactor*TMath::Sqrt(TMath::Abs(nphot)); - - // - // The following part is the commented first version of the error calculation - // Contact Markus Gaug for questions (or wait for the next documentation update...) - // - /* - nphotErr = signal > 0 ? signalErr*signalErr / (signal * signal) : 0. - + calibConv > 0 ? calibConvVar / (calibConv * calibConv ) : 0.; - nphotErr = TMath::Sqrt(nphotErr) * nphot; - */ - - MCerPhotPix *cpix = fCerPhotEvt->AddPixel(pixidx, nphot, nphotErr); - - if (sig.GetNumHiGainSaturated() > 0) - cpix->SetPixelHGSaturated(); - - if (sig.GetNumLoGainSaturated() > 0) - cpix->SetPixelSaturated(); - } - - fCerPhotEvt->FixSize(); - fCerPhotEvt->SetReadyToSave(); + if (data) + { + const MExtractedSignalPix &sig = (*fSignals)[pixidx]; + + Float_t signal = 0; + Float_t signalErr = 0.; + + if (sig.IsLoGainUsed()) + { + signal = sig.GetExtractedSignalLoGain()*hiloconv; + signalErr = signal*hiloconverr; + } + else + { + if (sig.GetExtractedSignalHiGain() <= 9999.) + signal = sig.GetExtractedSignalHiGain(); + } + + const Float_t nphot = signal*calibConv; + const Float_t nphotErr = calibFFactor*TMath::Sqrt(TMath::Abs(nphot)); + + // + // The following part is the commented first version of the error calculation + // Contact Markus Gaug for questions (or wait for the next documentation update...) + // + /* + nphotErr = signal > 0 ? signalErr*signalErr / (signal * signal) : 0. + + calibConv > 0 ? calibConvVar / (calibConv * calibConv ) : 0.; + nphotErr = TMath::Sqrt(nphotErr) * nphot; + */ + + MCerPhotPix *cpix = fCerPhotEvt->AddPixel(pixidx, nphot, nphotErr); + + if (sig.GetNumHiGainSaturated() > 0) + cpix->SetPixelHGSaturated(); + + if (sig.GetNumLoGainSaturated() > 0) + cpix->SetPixelSaturated(); + } + if (pedestal) + { + const MPedestalPix &ped = (*fPedestal)[pixidx]; + + // pedestals/(used FADC slices) in [ADC] counts + const Float_t pedes = ped.GetPedestal() * slices; + const Float_t pedrms = ped.GetPedestalRms() * sqrtslices; + + // + // pedestals/(used FADC slices) in [number of photons] + // + const Float_t pedphot = pedes * calibConv; + const Float_t pedphotrms = pedrms * calibConv; + + (*fPedPhot)[pixidx].Set(pedphot, pedphotrms); + } + } + + if (pedestal) + fPedPhot->SetReadyToSave(); + + if (data) + { + fCerPhotEvt->FixSize(); + fCerPhotEvt->SetReadyToSave(); + } } @@ -563,12 +545,6 @@ */ - if (fCalibrationMode!=kSkip) - CalibrateData(); - - if (TestPedestalFlag(kEvent)) - if (!CalibratePedestal()) - return kFALSE; - - return kTRUE; + Calibrate(fCalibrationMode!=kSkip, TestPedestalFlag(kEvent)); + return kTRUE; } Index: trunk/MagicSoft/Mars/mcalib/MCalibrateData.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrateData.h (revision 4751) +++ trunk/MagicSoft/Mars/mcalib/MCalibrateData.h (revision 4752) @@ -53,8 +53,7 @@ TString fNamePedPhotContainer; // name of fPedPhot - Bool_t CalibratePedestal(); - void CalibrateData(); + void Calibrate(Bool_t data, Bool_t pedestal) const; - Bool_t GetConversionFactor(UInt_t, Float_t &, Float_t &, Float_t &, Float_t &, Float_t &); + Bool_t GetConversionFactor(UInt_t, Float_t &, Float_t &, Float_t &, Float_t &, Float_t &) const; Int_t PreProcess(MParList *pList); Index: trunk/MagicSoft/Mars/mjobs/MJCalibrateSignal.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJCalibrateSignal.cc (revision 4751) +++ trunk/MagicSoft/Mars/mjobs/MJCalibrateSignal.cc (revision 4752) @@ -355,4 +355,5 @@ *fLog << all << GetDescriptor() << ": Done." << endl; + *fLog << endl << endl; return kTRUE; Index: trunk/MagicSoft/Mars/mjobs/MJPedestal.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJPedestal.cc (revision 4751) +++ trunk/MagicSoft/Mars/mjobs/MJPedestal.cc (revision 4752) @@ -650,4 +650,5 @@ *fLog << all << GetDescriptor() << ": Done." << endl; + *fLog << endl << endl; return kTRUE; Index: trunk/MagicSoft/Mars/mpedestal/MPedPhotCam.cc =================================================================== --- trunk/MagicSoft/Mars/mpedestal/MPedPhotCam.cc (revision 4751) +++ trunk/MagicSoft/Mars/mpedestal/MPedPhotCam.cc (revision 4752) @@ -39,7 +39,8 @@ #include "MPedPhotCam.h" -#include -#include #include + +#include "MArrayI.h" +#include "MArrayD.h" #include "MLog.h" @@ -233,10 +234,11 @@ const Int_t na = GetNumAreas(); - TArrayI acnt(na); - TArrayI scnt(ns); - TArrayD asumx(na); - TArrayD ssumx(ns); - TArrayD asumr(na); - TArrayD ssumr(ns); + // Using MArray instead of TArray because they don't do range checks + MArrayI acnt(na); + MArrayI scnt(ns); + MArrayD asumx(na); + MArrayD ssumx(ns); + MArrayD asumr(na); + MArrayD ssumr(ns); for (int i=0; i0) - { - pederr = rms/TMath::Sqrt((Float_t)nevt); - rmserr = rms/TMath::Sqrt((Float_t)nevt)/2.; - } - else - { - pederr = -1; - rmserr = -1; - } - switch (type) { case 0: - val = ped;// (*this)[idx].GetMean(); + val = (*this)[idx].GetMean(); break; case 1: - val = rms; // (*this)[idx].GetRms(); + val = (*this)[idx].GetRms(); break; case 2: - val = pederr; // new + val = (*this)[idx].GetNumEvents()>0 ? (*this)[idx].GetRms()/TMath::Sqrt((Float_t)(*this)[idx].GetNumEvents()) : -1; break; case 3: - val = rmserr; // new + val = (*this)[idx].GetNumEvents()>0 ? (*this)[idx].GetRms()/TMath::Sqrt((Float_t)(*this)[idx].GetNumEvents())/2. : -1; break; default: Index: trunk/MagicSoft/Mars/msignal/MExtractTimeFastSpline.cc =================================================================== --- trunk/MagicSoft/Mars/msignal/MExtractTimeFastSpline.cc (revision 4751) +++ trunk/MagicSoft/Mars/msignal/MExtractTimeFastSpline.cc (revision 4752) @@ -200,6 +200,6 @@ pp = fHiGainSecondDeriv[i-1] + 4.; fHiGainSecondDeriv[i] = -1.0/pp; - fHiGainFirstDeriv [i] = *(p+1) - 2.* *(p) + *(p-1); - fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp; + const Double_t deriv = *(p+1) - 2.* *(p) + *(p-1); + fHiGainFirstDeriv [i] = (6.0*deriv-fHiGainFirstDeriv[i-1])/pp; } @@ -232,4 +232,7 @@ // interval maxpos+1. // + + Float_t higainklo = fHiGainSecondDeriv[klo]; + Float_t higainkhi = fHiGainSecondDeriv[khi]; while (x abmax) @@ -265,4 +268,6 @@ khicont = (Float_t)*(first+khi); + higainklo = fHiGainSecondDeriv[klo]; + higainkhi = fHiGainSecondDeriv[khi]; while (x abmax) @@ -295,4 +300,6 @@ step = 0.04; // step size of 83 ps + higainklo = fHiGainSecondDeriv[klo]; + higainkhi = fHiGainSecondDeriv[khi]; while (x abmax) @@ -329,4 +336,6 @@ b = x - lower; + higainklo = fHiGainSecondDeriv[klo]; + higainkhi = fHiGainSecondDeriv[khi]; while (x>lo) { @@ -338,6 +347,6 @@ y = a* klocont + b* khicont - + (a*a*a-a)*fHiGainSecondDeriv[klo] - + (b*b*b-b)*fHiGainSecondDeriv[khi]; + + (a*a*a-a)*higainklo + + (b*b*b-b)*higainkhi; if (y > abmax)