Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 3912) +++ trunk/MagicSoft/Mars/Changelog (revision 3913) @@ -46,4 +46,13 @@ * msignal/MExtractPINDiode.[h,cc] - now deriving from MExtractor + + * msignal/Makefile + * msignal/SignalLinkDef.h + * msignal/MExtractTimeSpline.[h,cc] + * msignal/MExtractTimeHighestIntegral.[h,cc] + - replacements for MArrivalTimeCalc and MArrivalTimeCalc2 + + * mjobs/MJCalibration.[h,cc] + - Set time extractor by pointer Index: trunk/MagicSoft/Mars/macros/calibration.C =================================================================== --- trunk/MagicSoft/Mars/macros/calibration.C (revision 3912) +++ trunk/MagicSoft/Mars/macros/calibration.C (revision 3913) @@ -55,24 +55,22 @@ // ///////////////////////////////////////////////////////////////////////////// -//static const TString inpath = "/home/rootdata/Calib/2004_04_16/"; -static const TString inpath = "/mnt/Data/rootdata/CrabNebula/2003_12_01/"; -static const Int_t pedrun = 3340; -static const Int_t calrun1 = 3341; +static const TString inpath = "/mnt/Data/rootdata/Miscellaneous/2004_03_03/"; +static const Int_t pedrun = 20123; +static const Int_t calrun1 = 20126; static const Int_t calrun2 = 0; -/* -static const Int_t pedrun = 22265; -static const Int_t calrun1 = 22300; -static const Int_t calrun2 = 0; -*/ + void calibration(const Int_t prun=pedrun, const Int_t crun1=calrun1, const Int_t crun2=calrun2) { MExtractFixedWindowPeakSearch extractor; - extractor.SetRange(0,6,0,6); + // extractor.SetRange(0,14,0,14); - /* - MExtractTimeSpline timeext; - timeext.SetRange(0,14,6,14); - */ + // MExtractSlidingWindow extractor; + // MExtractFixedWindow extractor; + + MExtractTimeHighestIntegral timeext; + + // MExtractTimeSpline timeext; + // timeext.SetRange(0,14,6,14); MRunIter pruns; @@ -137,6 +135,5 @@ calloop.SetRelTimeCalibration(); calloop.SetExtractor(&extractor); - // calloop.SetArrivalTimeLevel(1); - // calloop.SetTimeExtractor(&timeext); + calloop.SetTimeExtractor(&timeext); calloop.SetInput(&cruns); calloop.SetDisplay(display); Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimeCam.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimeCam.cc (revision 3912) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimeCam.cc (revision 3913) @@ -159,4 +159,5 @@ { + fCam = (MCalibrationCam*)pList->FindObject("MCalibrationRelTimeCam"); if (!fCam) @@ -172,4 +173,5 @@ } + MArrivalTimeCam *signal = (MArrivalTimeCam*)pList->FindObject("MArrivalTimeCam"); if (!signal) @@ -193,5 +195,4 @@ } } - if (fLoGainArray->GetEntries()==0) Index: trunk/MagicSoft/Mars/mjobs/MJCalibration.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJCalibration.cc (revision 3912) +++ trunk/MagicSoft/Mars/mjobs/MJCalibration.cc (revision 3913) @@ -39,5 +39,8 @@ // // Different signal extractors can be set with the command SetExtractor() -// Only extractor deriving from MExtractor can be set, default is MExtractSlidingWindow +// Only extractors deriving from MExtractor can be set, default is MExtractSlidingWindow +// +// Different arrival time extractors can be set with the command SetTimeExtractor() +// Only extractors deriving from MTimeExtractor can be set, default is MExtractTimeSpline // // At the end of the eventloop, plots and results are displayed, depending on @@ -116,8 +119,10 @@ #include "MBadPixelsMerge.h" #include "MBadPixelsCam.h" +#include "MTimeExtractor.h" #include "MExtractor.h" #include "MExtractPINDiode.h" #include "MExtractBlindPixel.h" #include "MExtractSlidingWindow.h" +#include "MExtractTimeSpline.h" #include "MFCosmics.h" #include "MContinue.h" @@ -139,10 +144,11 @@ // Default constructor. // -// Sets fRuns to 0, fExtractor to NULL, fColor to kNONE, fDisplay to kNormalDisplay, -// fRelTime to kFALSE, fDataCheck to kFALSE +// Sets fRuns to 0, fExtractor to NULL, fTimeExtractor to NULL, fColor to kNONE, +// fDisplay to kNormalDisplay, fRelTime to kFALSE, fDataCheck to kFALSE // MJCalibration::MJCalibration(const char *name, const char *title) - : fRuns(0), fExtractor(NULL), fColor(MCalibrationCam::kNONE), fDisplayType(kNormalDisplay), - fRelTimes(kFALSE), fDataCheck(kFALSE), fArrivalTimeLevel(2) + : fRuns(0), fExtractor(NULL), fTimeExtractor(NULL), + fColor(MCalibrationCam::kNONE), fDisplayType(kNormalDisplay), + fRelTimes(kFALSE), fDataCheck(kFALSE) { fName = name ? name : "MJCalibration"; @@ -761,5 +767,5 @@ // 5) MExtractPINDiode // 6) MExtractBlindPixel -// 7) MArrivalTimeCalc, MArrivalTimeCalc2, depending on fArrivalTimeLevel (only if flag fRelTimes is chosen) +// 7) MTimeExtractor (only if flag fRelTimes is chosen) // 8) MContinue(MFCosmics) // 9) MFillH("MHCalibrationChargePINDiode", "MExtractedSignalPINDiode") @@ -833,6 +839,5 @@ MExtractBlindPixel blindext; MExtractSlidingWindow extract2; - MArrivalTimeCalc tmecalc1; - MArrivalTimeCalc2 tmecalc2; + MExtractTimeSpline timespline; MCalibrationChargeCalc calcalc; @@ -879,14 +884,12 @@ { plist.AddToList(&fRelTimeCam); - if (fArrivalTimeLevel <= 1) - tlist.AddToList(&tmecalc1); - else if (fArrivalTimeLevel == 2) - tlist.AddToList(&tmecalc2); + + if (fTimeExtractor) + tlist.AddToList(fTimeExtractor); else { - *fLog << err << GetDescriptor() - << ": No valid Signal ArrivalTime Level has been chosen, have only: " << fArrivalTimeLevel - << " aborting..." << endl; - return kFALSE; + *fLog << warn << GetDescriptor() + << ": No extractor has been chosen, take default MTimeExtractSpine " << endl; + tlist.AddToList(×pline); } } @@ -900,5 +903,5 @@ if (fRelTimes) tlist.AddToList(&filltme); - + // Create and setup the eventloop MEvtLoop evtloop(fName); Index: trunk/MagicSoft/Mars/mjobs/MJCalibration.h =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJCalibration.h (revision 3912) +++ trunk/MagicSoft/Mars/mjobs/MJCalibration.h (revision 3913) @@ -22,4 +22,5 @@ class MPedestalCam; class MExtractor; +class MTimeExtractor; class MJCalibration : public MParContainer, public MGCamDisplays { @@ -30,6 +31,7 @@ TString fOutputPath; // Path to the output files - MRunIter *fRuns; // Calibration files - MExtractor *fExtractor; // Signal extractor + MRunIter *fRuns; // Calibration files + MExtractor *fExtractor; // Signal extractor + MTimeExtractor *fTimeExtractor; // Arrival Time extractor MBadPixelsCam fBadPixels; // Bad Pixels cam, can be set from previous runs @@ -46,5 +48,4 @@ Bool_t fRelTimes; // Flag if relative times have to be calibrated Bool_t fDataCheck; // Flag if the data check is run on raw data - UInt_t fArrivalTimeLevel; // Level arr. time extractor (e.g. MArrivalTimeCalc2) void DisplayResult(MParList &plist); @@ -69,4 +70,5 @@ void SetBadPixels(const MBadPixelsCam &bad) { bad.Copy(fBadPixels); } void SetExtractor(MExtractor* ext) { fExtractor = ext; } + void SetTimeExtractor(MTimeExtractor* ext) { fTimeExtractor = ext; } void SetQECam (const MCalibrationQECam &qe) { qe.Copy(fQECam); } void SetColor (const MCalibrationCam::PulserColor_t color) { fColor = color; } @@ -83,7 +85,4 @@ void SetDataCheck (const Bool_t b=kTRUE) { fDataCheck = b; SetDataCheckDisplay(); } - // Extractors - void SetArrivalTimeLevel (const UInt_t i=1 ) { fArrivalTimeLevel = i; } - Bool_t ReadCalibrationCam(); Bool_t ProcessFile( MPedestalCam &pedcam ); Index: trunk/MagicSoft/Mars/msignal/MExtractTimeHighestIntegral.cc =================================================================== --- trunk/MagicSoft/Mars/msignal/MExtractTimeHighestIntegral.cc (revision 3913) +++ trunk/MagicSoft/Mars/msignal/MExtractTimeHighestIntegral.cc (revision 3913) @@ -0,0 +1,298 @@ +/* ======================================================================== *\ +! +! * +! * 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, 02/2004 +! Hendrik Bartko, 02/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MExtractTimeHighestIntegral +// +// Calculates the arrival time as the mean time of the fWindowSize time slices +// which have the highest integral content. +// +// +////////////////////////////////////////////////////////////////////////////// +#include "MExtractTimeHighestIntegral.h" + +#include "MLog.h" +#include "MLogManip.h" +#include "MPedestalPix.h" + +ClassImp(MExtractTimeHighestIntegral); + +using namespace std; + +const Byte_t MExtractTimeHighestIntegral::fgHiGainFirst = 0; +const Byte_t MExtractTimeHighestIntegral::fgHiGainLast = 14; +const Byte_t MExtractTimeHighestIntegral::fgLoGainFirst = 3; +const Byte_t MExtractTimeHighestIntegral::fgLoGainLast = 14; +const Byte_t MExtractTimeHighestIntegral::fgHiGainWindowSize = 6; +const Byte_t MExtractTimeHighestIntegral::fgLoGainWindowSize = 6; +// -------------------------------------------------------------------------- +// +// Default constructor. +// +// Calls: +// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast) +// +MExtractTimeHighestIntegral::MExtractTimeHighestIntegral(const char *name, const char *title) + : fHiGainWindowSize(fgHiGainWindowSize), + fLoGainWindowSize(fgLoGainWindowSize) +{ + + fName = name ? name : "MExtractTimeHighestIntegral"; + fTitle = title ? title : "Task to extract the signal from the FADC slices"; + + SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast); +} + +// -------------------------------------------------------------------------- +// +// SetRange: +// +// Calls: +// - MExtractor::SetRange(hifirst,hilast,lofirst,lolast); +// +void MExtractTimeHighestIntegral::SetRange(Byte_t hifirst, Byte_t hilast, Byte_t lofirst, Byte_t lolast) +{ + + MExtractor::SetRange(hifirst,hilast,lofirst,lolast); + + + Int_t range = fHiGainLast - fHiGainFirst + 1; + + if (range < 2) + { + *fLog << warn << GetDescriptor() + << Form("%s%2i%s%2i%s",": Hi-Gain Extraction range [",(int)fHiGainFirst,"," + ,fHiGainLast,"] too small, ") << endl; + *fLog << warn << GetDescriptor() + << " will move higher limit to obtain 4 slices " << endl; + SetRange(fHiGainFirst, fHiGainLast+4-range,fLoGainFirst,fLoGainLast); + } + + range = fLoGainLast - fLoGainFirst + 1; + + if (range < 2) + { + *fLog << warn << GetDescriptor() + << Form("%s%2i%s%2i%s",": Lo-Gain Extraction range [",(int)fLoGainFirst,"," + ,fLoGainLast,"] too small, ") << endl; + *fLog << warn << GetDescriptor() + << " will move logher limit to obtain 4 slices " << endl; + SetRange(fHiGainFirst, fHiGainLast,fLoGainFirst,fLoGainLast+4-range); + } + + + SetWindowSize(fHiGainWindowSize,fLoGainWindowSize); + + fNumHiGainSamples = fHiGainLast-fHiGainFirst+1; + fNumLoGainSamples = fLoGainLast-fLoGainFirst+1; + +} + +// -------------------------------------------------------------------------- +// +// 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 +// +void MExtractTimeHighestIntegral::SetWindowSize(Byte_t windowh, Byte_t windowl) +{ + + fHiGainWindowSize = windowh & ~1; + fLoGainWindowSize = windowl & ~1; + + if (fHiGainWindowSize != windowh) + *fLog << warn << GetDescriptor() << ": Hi Gain window size has to be even, set to: " + << int(fHiGainWindowSize) << " samples " << endl; + + if (fLoGainWindowSize != windowl) + *fLog << warn << GetDescriptor() << ": Lo Gain window size has to be even, set to: " + << int(fLoGainWindowSize) << " samples " << endl; + + const Byte_t availhirange = (fHiGainLast-fHiGainFirst+1) & ~1; + const Byte_t availlorange = (fLoGainLast-fLoGainFirst+1) & ~1; + + if (fHiGainWindowSize > availhirange) + { + *fLog << warn << GetDescriptor() + << Form("%s%2i%s%2i%s%2i%s",": Hi Gain window size: ",(int)fHiGainWindowSize, + " is bigger than available range: [",(int)fHiGainFirst,",",(int)fHiGainLast,"]") << endl; + *fLog << warn << GetDescriptor() + << ": Will set window size to: " << (int)availhirange << endl; + fHiGainWindowSize = availhirange; + } + + if (fLoGainWindowSize > availlorange) + { + *fLog << warn << GetDescriptor() + << Form("%s%2i%s%2i%s%2i%s",": Lo Gain window size: ",(int)fLoGainWindowSize, + " is bigger than available range: [",(int)fLoGainFirst,",",(int)fLoGainLast,"]") << endl; + *fLog << warn << GetDescriptor() + << ": Will set window size to: " << (int)availlorange << endl; + fLoGainWindowSize= availlorange; + } + + fHiGainWindowSizeSqrt = TMath::Sqrt((Float_t)fHiGainWindowSize); + fLoGainWindowSizeSqrt = TMath::Sqrt((Float_t)fLoGainWindowSize); +} + + +void MExtractTimeHighestIntegral::FindTimeHiGain(Byte_t *ptr, Float_t &time, Float_t &deltatime, Byte_t &sat, const MPedestalPix &ped) const +{ + + const Byte_t *end = ptr + fHiGainLast - fHiGainFirst +1 ; + + Int_t sum=0; // integral content of the actual window + Int_t max = 0; // highest integral content of all windows + + // + // Calculate the sum of the first fWindowSize slices + // + sat = 0; + Byte_t *p = ptr; + + while (p= fSaturationLimit) + sat++; + } + + // + // Check for saturation in all other slices + // + while (p= fSaturationLimit) + sat++; + + // + // Calculate the i-th sum as + // sum_i+1 = sum_i + slice[i+8] - slice[i] + // This is fast and accurate (because we are using int's) + // + max=sum; + Byte_t *ptrmax=ptr; // pointer to the first slice of the maximum window + + for (p=ptr; p+fHiGainWindowSizemax) + { + max = sum; + ptrmax = p; + } + } + + // now calculate the time for the maximum window + Int_t timesignalsum = 0; + Int_t timesquaredsum =0; + Int_t timesum =0; + + for (p=ptrmax; p < ptrmax + fHiGainWindowSize; p++) + { + timesignalsum += *p*(p-ptr); + timesum += p-ptr; + timesquaredsum += (p-ptr)*(p-ptr); + } + + const Float_t pedes = ped.GetPedestal(); + const Float_t pedrms = ped.GetPedestalRms(); + const Float_t pedsubsum = max - fHiGainWindowSize*pedes; + const Float_t pedsubtimesignalsum = timesignalsum - timesum*pedes; + + time = pedsubsum != 0 ? pedsubtimesignalsum / pedsubsum : 1; + deltatime = pedsubsum != 0 ? pedrms / pedsubsum * sqrt(timesquaredsum - fHiGainWindowSize*time) : 1; +} + +void MExtractTimeHighestIntegral::FindTimeLoGain(Byte_t *ptr, Float_t &time, Float_t &deltatime, Byte_t &sat, const MPedestalPix &ped) const +{ + + const Byte_t *end = ptr + fLoGainLast - fLoGainFirst +1 ; + + Int_t sum=0; // integral content of the actual window + Int_t max = 0; // highest integral content of all windows + + // + // Calculate the sum of the first fWindowSize slices + // + sat = 0; + Byte_t *p = ptr; + + while (p= fSaturationLimit) + sat++; + } + + // + // Check for saturation in all other slices + // + while (p= fSaturationLimit) + sat++; + + // + // Calculate the i-th sum as + // sum_i+1 = sum_i + slice[i+8] - slice[i] + // This is fast and accurate (because we are using int's) + // + max=sum; + Byte_t *ptrmax=ptr; // pointer to the first slice of the maximum window + + for (p=ptr; p+fLoGainWindowSizemax) + { + max = sum; + ptrmax = p; + } + } + + // now calculate the time for the maximum window + Int_t timesignalsum = 0; + Int_t timesquaredsum =0; + Int_t timesum =0; + + for (p=ptrmax; p < ptrmax + fLoGainWindowSize; p++) + { + timesignalsum += *p*(p-ptr); + timesum += p-ptr; + timesquaredsum += (p-ptr)*(p-ptr); + } + + const Float_t pedes = ped.GetPedestal(); + const Float_t pedrms = ped.GetPedestalRms(); + const Float_t pedsubsum = max - fLoGainWindowSize*pedes; + const Float_t pedsubtimesignalsum = timesignalsum - timesum*pedes; + + time = pedsubsum != 0 ? pedsubtimesignalsum / pedsubsum : 1; + deltatime = pedsubsum != 0 ? pedrms / pedsubsum * sqrt(timesquaredsum - fLoGainWindowSize*time) : 1; +} + Index: trunk/MagicSoft/Mars/msignal/MExtractTimeHighestIntegral.h =================================================================== --- trunk/MagicSoft/Mars/msignal/MExtractTimeHighestIntegral.h (revision 3913) +++ trunk/MagicSoft/Mars/msignal/MExtractTimeHighestIntegral.h (revision 3913) @@ -0,0 +1,39 @@ +#ifndef MARS_MExtractTimeHighestIntegral +#define MARS_MExtractTimeHighestIntegral + +#ifndef MARS_MTimeExtractor +#include "MTimeExtractor.h" +#endif + +class MExtractTimeHighestIntegral : public MTimeExtractor +{ +private: + + static const Byte_t fgHiGainFirst; + static const Byte_t fgHiGainLast; + static const Byte_t fgLoGainFirst; + static const Byte_t fgLoGainLast; + static const Byte_t fgHiGainWindowSize; + static const Byte_t fgLoGainWindowSize; + + Byte_t fHiGainWindowSize; // Number of gains in window + Float_t fHiGainWindowSizeSqrt; // Sqaure root of number of gains in window + Byte_t fLoGainWindowSize; // Number of gains in window + Float_t fLoGainWindowSizeSqrt; // Sqaure root of number of gains in window + + void FindTimeHiGain(Byte_t *first, Float_t &time, Float_t &dtime, Byte_t &sat, const MPedestalPix &ped) const; + void FindTimeLoGain(Byte_t *first, Float_t &time, Float_t &dtime, Byte_t &sat, const MPedestalPix &ped) const; + +public: + MExtractTimeHighestIntegral(const char *name=NULL, const char *title=NULL); + + void SetRange(Byte_t hifirst=0, Byte_t hilast=0, Byte_t lofirst=0, Byte_t lolast=0); + void SetWindowSize(Byte_t windowh=fgHiGainWindowSize, + Byte_t windowl=fgLoGainWindowSize); + + ClassDef(MExtractTimeHighestIntegral, 0) // Task to Extract the Arrival Times As the mean time of the fWindowSize time slices +}; + +#endif + + Index: trunk/MagicSoft/Mars/msignal/Makefile =================================================================== --- trunk/MagicSoft/Mars/msignal/Makefile (revision 3912) +++ trunk/MagicSoft/Mars/msignal/Makefile (revision 3913) @@ -50,4 +50,5 @@ MTimeExtractor.cc \ MExtractTimeSpline.cc \ + MExtractTimeHighestIntegral.cc \ MArrivalTime.cc \ MArrivalTimeCalc.cc \ Index: trunk/MagicSoft/Mars/msignal/SignalLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/msignal/SignalLinkDef.h (revision 3912) +++ trunk/MagicSoft/Mars/msignal/SignalLinkDef.h (revision 3913) @@ -24,4 +24,5 @@ #pragma link C++ class MTimeExtractor+; #pragma link C++ class MExtractTimeSpline+; +#pragma link C++ class MExtractTimeHighestIntegral+; #pragma link C++ class MArrivalTimeCam; #pragma link C++ class MArrivalTimePix;