Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 3080) +++ trunk/MagicSoft/Mars/Changelog (revision 3081) @@ -8,4 +8,12 @@ * mcalib/MCalibrationCam.cc - fixed documentation + + * mhbase/MH.[h,cc] + - new function ProjectArray + + * mcalib/MHCalibrationPixel.[h,cc] + * mcalib/MHCalibrationBlindPixel.[h,cc] + - use ProjectArray from MH to plot the projection of the fourier + spectrum Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.cc (revision 3080) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.cc (revision 3081) @@ -88,6 +88,4 @@ const Double_t MHCalibrationBlindPixel::fgBlindPixelElectronicAmpError = 0.002; -const Axis_t MHCalibrationBlindPixel::fNyquistFreq = 1.0; -const Axis_t MHCalibrationBlindPixel::fMinFreq = 0.; const Int_t MHCalibrationBlindPixel::fPSDNbins = 30; @@ -266,29 +264,22 @@ fPSDHiGain = fourier.PowerSpectrumDensity(fHiGains); - Int_t entries; - TArrayF *array; - - fHPSD = new TH1F("HBlindPixelPSD", - "Power Spectrum Density Projection Blind Pixel", - fPSDNbins,fMinFreq,fNyquistFreq); - - array = fPSDHiGain; - entries = array->GetSize(); - - for (Int_t i=0;iFill(array->At(i)); - + fHPSD = MH::ProjectArray(fPSDHiGain, fPSDNbins, + "PSDProjectionBlindPixel", + "Power Spectrum Density Projection HiGain Blind Pixel"); // // First guesses for the fit (should be as close to reality as possible, // - const Double_t area_guess = entries*10.; - - fPSDExpFit = new TF1("PSDExpFit","[0]*exp(-[1]*x)",0.,1.); - - fPSDExpFit->SetParameters(entries,10.); - fPSDExpFit->SetParNames("Area","slope"); - fPSDExpFit->SetParLimits(0,0.,3.*area_guess); - fPSDExpFit->SetParLimits(1,0.,20.); - fPSDExpFit->SetRange(fMinFreq,fNyquistFreq); + const Double_t xmax = fHPSD->GetXaxis()->GetXmax(); + + fPSDExpFit = new TF1("PSDExpFit","exp([0]-[1]*x)",0.,xmax); + + const Double_t slope_guess = (TMath::Log(fHPSD->GetEntries())+1.)/xmax; + const Double_t offset_guess = slope_guess*xmax; + + fPSDExpFit->SetParameters(offset_guess, slope_guess); + fPSDExpFit->SetParNames("Offset","Slope"); + fPSDExpFit->SetParLimits(0,offset_guess/2.,2.*offset_guess); + fPSDExpFit->SetParLimits(1,slope_guess/2.,2.*slope_guess); + fPSDExpFit->SetRange(0.,xmax); fHPSD->Fit(fPSDExpFit,"RQL0"); Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.h (revision 3080) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.h (revision 3081) @@ -25,6 +25,4 @@ static const Double_t fgBlindPixelElectronicAmpError; - static const Axis_t fNyquistFreq; - static const Axis_t fMinFreq; static const Int_t fPSDNbins; @@ -40,5 +38,5 @@ TArrayF* fPSDLoGain; //-> Power spectrum density of fLoGains - TH1F* fHPSD; //-> + TH1I* fHPSD; //-> TF1* fPSDExpFit; //-> Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationPixel.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationPixel.cc (revision 3080) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationPixel.cc (revision 3081) @@ -61,6 +61,5 @@ const Int_t MHCalibrationPixel::fNDFLimit = 5; -const Axis_t MHCalibrationPixel::fMinFreq = 0.; -const Int_t MHCalibrationPixel::fPSDNbins = 50; +const Int_t MHCalibrationPixel::fPSDNbins = 30; // -------------------------------------------------------------------------- @@ -268,45 +267,28 @@ fPSDHiGain = fourier.PowerSpectrumDensity(fHiGains); - Int_t entries; - TArrayF *array; - - const Double_t max_freq = (fChargeSigma > 0.) ? fChargeSigma*0.2 : 2.; - const Double_t min_freq = 0.; - if (IsUseLoGain()) - { - fHPSD = new TH1F(Form("%s%d","HPSD",fPixId), - Form("%s%s","Power Spectrum Density Projection ","LoGain"), - fPSDNbins,min_freq,max_freq); - - array = fPSDLoGain; - } + fHPSD = MH::ProjectArray(fPSDLoGain, fPSDNbins, + Form("%s%d","PSDProjection",fPixId), + Form("%s%d","Power Spectrum Density Projection LoGain ",fPixId)); else - { - - fHPSD = new TH1F(Form("%s%d","HPSD",fPixId), - Form("%s%s","Power Spectrum Density Projection ","HiGain"), - fPSDNbins,min_freq,max_freq); - - array = fPSDHiGain; - } - - entries = array->GetSize(); - - for (Int_t i=0;iFill(array->At(i)); + fHPSD = MH::ProjectArray(fPSDHiGain, fPSDNbins, + Form("%s%d","PSDProjection",fPixId), + Form("%s%d","Power Spectrum Density Projection HiGain ",fPixId)); // // First guesses for the fit (should be as close to reality as possible, // - const Double_t area_guess = entries*10.; - - fPSDExpFit = new TF1(Form("%s%d","PSDExpFit",fPixId),"[0]*exp(-[1]*x)",0.,1.); - - fPSDExpFit->SetParameters(area_guess,10.); - fPSDExpFit->SetParNames("Area","slope"); - fPSDExpFit->SetParLimits(0,0.,20.*area_guess); - fPSDExpFit->SetParLimits(1,0.,30.); - fPSDExpFit->SetRange(min_freq,max_freq); + const Double_t xmax = fHPSD->GetXaxis()->GetXmax(); + + fPSDExpFit = new TF1(Form("%s%d","PSDExpFit",fPixId),"exp([0]-[1]*x)",0.,xmax); + + const Double_t slope_guess = (TMath::Log(fHPSD->GetEntries())+1.)/xmax; + const Double_t offset_guess = slope_guess*xmax; + + fPSDExpFit->SetParameters(offset_guess, slope_guess); + fPSDExpFit->SetParNames("Offset","Slope"); + fPSDExpFit->SetParLimits(0,offset_guess/2.,2.*offset_guess); + fPSDExpFit->SetParLimits(1,slope_guess/2.,2.*slope_guess); + fPSDExpFit->SetRange(0.,xmax); fHPSD->Fit(fPSDExpFit,"RQL0"); Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationPixel.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationPixel.h (revision 3080) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationPixel.h (revision 3081) @@ -26,5 +26,4 @@ static const Int_t fNDFLimit; - static const Axis_t fMinFreq; static const Int_t fPSDNbins; @@ -49,5 +48,5 @@ TF1* fChargeGausFit; //-> - TH1F* fHPSD; //-> + TH1I* fHPSD; //-> TF1* fPSDExpFit; //-> Index: trunk/MagicSoft/Mars/mhbase/MH.cc =================================================================== --- trunk/MagicSoft/Mars/mhbase/MH.cc (revision 3080) +++ trunk/MagicSoft/Mars/mhbase/MH.cc (revision 3081) @@ -1081,2 +1081,86 @@ return NULL; } + +TH1I* MH::ProjectArray(const TArrayF *array, Int_t nbins, const char* name, const char* title) +{ + + const Int_t size = array->GetSize(); + TH1I *h1=0; + + //check if histogram with identical name exist + TObject *h1obj = gROOT->FindObject(name); + if (h1obj && h1obj->InheritsFrom("TH1I")) { + h1 = (TH1I*)h1obj; + h1->Reset(); + } + + Double_t max = 0.; + Double_t min = 0.; + Int_t newbins = 0; + + // first loop over array to find the maximum: + for (Int_t i=0; iAt(i) > max) + max = array->At(i); + + FindGoodLimits(nbins, newbins, min, max, kFALSE); + + if (!h1) + { + h1 = new TH1I(name, title, newbins, min, max); + h1->SetXTitle(""); + h1->SetYTitle("Counts"); + h1->SetDirectory(NULL); + } + + // Second loop to fill the histogram + for (Int_t i=0;iFill(array->At(i)); + + h1->SetEntries(size); + + return h1; +} + +TH1I* MH::ProjectArray(const TArrayD *array, Int_t nbins, const char* name, const char* title) +{ + + const Int_t size = array->GetSize(); + + TH1I *h1=0; + + //check if histogram with identical name exist + TObject *h1obj = gROOT->FindObject(name); + if (h1obj && h1obj->InheritsFrom("TH1I")) { + h1 = (TH1I*)h1obj; + h1->Reset(); + } + + Double_t max = 0.; + Double_t min = 0.; + Int_t newbins = 0; + + // first loop over array to find the maximum: + for (Int_t i=0; iAt(i) > max) + max = array->At(i); + + FindGoodLimits( nbins, newbins, min, max, kFALSE); + + if (!h1) + { + h1 = new TH1I(name, title, newbins, min, max); + h1->SetXTitle(""); + h1->SetYTitle("Counts"); + h1->SetDirectory(NULL); + } + + // Second loop to fill the histogram + for (Int_t i=0;iFill(array->At(i)); + + h1->SetEntries(size); + + return h1; +} + Index: trunk/MagicSoft/Mars/mhbase/MH.h =================================================================== --- trunk/MagicSoft/Mars/mhbase/MH.h (revision 3080) +++ trunk/MagicSoft/Mars/mhbase/MH.h (revision 3081) @@ -7,8 +7,10 @@ class TH1; +class TH1I; class TH1D; class TH2; class TH3; class TAxis; +class TArrayF; class TArrayD; class TCanvas; @@ -89,4 +91,9 @@ static Double_t GetMinimumGT(const TH1 &h, Double_t gt=0); static Int_t CutEdges(TH1 *h, Int_t nbins); + + static TH1I* ProjectArray(const TArrayF *array, Int_t nbins=30, const char* name="ProjectArray", + const char* title="Projected Array"); + static TH1I* ProjectArray(const TArrayD *array, Int_t nbins=30, const char* name="ProjectArray", + const char* title="Projected Array"); static void ProjectionX(TH1D &dest, const TH2 &src, Int_t firstybin=-1, Int_t lastybin=9999);