Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 3024) +++ trunk/MagicSoft/Mars/Changelog (revision 3025) @@ -4,4 +4,11 @@ -*-*- END OF LINE -*-*- + 2004/02/05: Markus Gaug + * mcalib/MHCalibrationPixel.[h,cc] + - put an NDF limit for the fit + - removed all TStrings and replace them by Form(..) + - Check for nan's of all returned fit values + + 2004/02/05: Abelardo Moralejo Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationPixel.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationPixel.cc (revision 3024) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationPixel.cc (revision 3025) @@ -52,12 +52,13 @@ const Int_t MHCalibrationPixel::fChargeNbinsHiGain = 2100; const Int_t MHCalibrationPixel::fChargeNbinsLoGain = 1010; +const Int_t MHCalibrationPixel::fChargevsNbins = 5000; const Int_t MHCalibrationPixel::fAbsTimeNbins = 32; -const Int_t MHCalibrationPixel::fRelTimeNbins = 480; -const Int_t MHCalibrationPixel::fChargevsNbins = 5000; const Axis_t MHCalibrationPixel::fAbsTimeFirst = - 0.25; const Axis_t MHCalibrationPixel::fAbsTimeLast = 15.75; -const Axis_t MHCalibrationPixel::fRelTimeFirst = -11.; -const Axis_t MHCalibrationPixel::fRelTimeLast = 12.; +const Int_t MHCalibrationPixel::fRelTimeNbins = 900; +const Axis_t MHCalibrationPixel::fRelTimeFirst = -13.; +const Axis_t MHCalibrationPixel::fRelTimeLast = 13.; const Float_t MHCalibrationPixel::fProbLimit = 0.001; +const Int_t MHCalibrationPixel::fNDFLimit = 5; // -------------------------------------------------------------------------- @@ -287,10 +288,10 @@ Bool_t MHCalibrationPixel::FillAbsTimeLoGain(Float_t t) { - return (fHAbsTimeLoGain->Fill(t) > -1); + return (fHAbsTimeLoGain->Fill(t)> -1); } Bool_t MHCalibrationPixel::FillRelTimeLoGain(Float_t t) { - return (fHRelTimeLoGain->Fill(t) > -1); + return (fHRelTimeLoGain->Fill(t) > -1); } @@ -323,37 +324,43 @@ { - fPixId = id; - - // - // Names Hi gain Histograms - // - fHChargeHiGain->SetName( Form("%s%d",fHChargeHiGain->GetName(), id)); - fHAbsTimeHiGain->SetName( Form("%s%d",fHAbsTimeHiGain->GetName(), id)); - fHRelTimeHiGain->SetName( Form("%s%d",fHRelTimeHiGain->GetName(), id)); - fHChargevsNHiGain->SetName(Form("%s%d",fHChargevsNHiGain->GetName(),id)); - - // - // Title Hi gain Histograms - // - fHChargeHiGain->SetTitle( Form("%s%d",fHChargeHiGain->GetTitle(), id)); - fHAbsTimeHiGain->SetTitle( Form("%s%d",fHAbsTimeHiGain->GetTitle(), id)); - fHRelTimeHiGain->SetTitle( Form("%s%d",fHRelTimeHiGain->GetTitle(), id)); - fHChargevsNHiGain->SetTitle( Form("%s%d",fHChargevsNHiGain->GetTitle(),id)); - - // - // Names Low Gain Histograms - // - fHChargeLoGain->SetName( Form("%s%d",fHChargeLoGain->GetName(),id)); - fHAbsTimeLoGain->SetName( Form("%s%d",fHAbsTimeLoGain->GetName(),id)); - fHRelTimeLoGain->SetName( Form("%s%d",fHRelTimeLoGain->GetName(),id)); - fHChargevsNLoGain->SetName( Form("%s%d",fHChargevsNLoGain->GetName(),id)); - - // - // Titles Low Gain Histograms - // - fHChargeLoGain->SetTitle( Form("%s%d",fHChargeLoGain->GetTitle(),id)); - fHAbsTimeLoGain->SetTitle( Form("%s%d",fHAbsTimeLoGain->GetTitle(),id)); - fHRelTimeLoGain->SetTitle( Form("%s%d",fHRelTimeLoGain->GetTitle(),id)); - fHChargevsNLoGain->SetTitle( Form("%s%d",fHChargevsNLoGain->GetTitle(),id)); + // Change only if the names have not yet been changed + if (fPixId == -1) + { + + // + // Names Hi gain Histograms + // + fHChargeHiGain->SetName( Form("%s%d",fHChargeHiGain->GetName(), id)); + fHAbsTimeHiGain->SetName( Form("%s%d",fHAbsTimeHiGain->GetName(), id)); + fHRelTimeHiGain->SetName( Form("%s%d",fHRelTimeHiGain->GetName(), id)); + fHChargevsNHiGain->SetName(Form("%s%d",fHChargevsNHiGain->GetName(),id)); + + // + // Title Hi gain Histograms + // + fHChargeHiGain->SetTitle( Form("%s%d",fHChargeHiGain->GetTitle(), id)); + fHAbsTimeHiGain->SetTitle( Form("%s%d",fHAbsTimeHiGain->GetTitle(), id)); + fHRelTimeHiGain->SetTitle( Form("%s%d",fHRelTimeHiGain->GetTitle(), id)); + fHChargevsNHiGain->SetTitle( Form("%s%d",fHChargevsNHiGain->GetTitle(),id)); + + // + // Names Low Gain Histograms + // + fHChargeLoGain->SetName( Form("%s%d",fHChargeLoGain->GetName(),id)); + fHAbsTimeLoGain->SetName( Form("%s%d",fHAbsTimeLoGain->GetName(),id)); + fHRelTimeLoGain->SetName( Form("%s%d",fHRelTimeLoGain->GetName(),id)); + fHChargevsNLoGain->SetName( Form("%s%d",fHChargevsNLoGain->GetName(),id)); + + // + // Titles Low Gain Histograms + // + fHChargeLoGain->SetTitle( Form("%s%d",fHChargeLoGain->GetTitle(),id)); + fHAbsTimeLoGain->SetTitle( Form("%s%d",fHAbsTimeLoGain->GetTitle(),id)); + fHRelTimeLoGain->SetTitle( Form("%s%d",fHRelTimeLoGain->GetTitle(),id)); + fHChargevsNLoGain->SetTitle( Form("%s%d",fHChargevsNLoGain->GetTitle(),id)); + + fPixId = id; + } + } @@ -673,5 +680,5 @@ fHivsLoGain->Fit("pol1","rq","",fHivsLoGain->GetMean()-2.5*fHivsLoGain->GetRMS(), - fHivsLoGain->GetMean()+2.0*fHivsLoGain->GetRMS()); + fHivsLoGain->GetMean()+2.5*fHivsLoGain->GetRMS()); fOffset = fHivsLoGain->GetFunction("pol1")->GetParameter(0); @@ -744,6 +751,6 @@ if (fRelTimeGausFit->GetProb() < 0.001) { - rmin = fRelTimeGausFit->GetParameter(1) - 2.*fRelTimeGausFit->GetParameter(2); - rmax = fRelTimeGausFit->GetParameter(1) + 2.*fRelTimeGausFit->GetParameter(2); + rmin = fRelTimeGausFit->GetParameter(1) - 1.5*fRelTimeGausFit->GetParameter(2); + rmax = fRelTimeGausFit->GetParameter(1) + 1.5*fRelTimeGausFit->GetParameter(2); fRelTimeGausFit->SetRange(rmin,rmax); hist->Fit(fRelTimeGausFit,option); @@ -765,5 +772,5 @@ } - if (TMath::IsNaN(fRelTimeChisquare) || (fRelTimeProb < gkProbLimit)) + if (TMath::IsNaN(fRelTimeChisquare) || (fRelTimeProb < fProbLimit)) { CLRBIT(fFlags,kTimeFitOK); @@ -788,28 +795,23 @@ // Axis_t rmin = fChargeFirstHiGain; + Axis_t rmax = fChargeLastHiGain; + TH1F *hist = fHChargeHiGain; + if (TESTBIT(fFlags,kUseLoGain)) - rmin = fChargeFirstLoGain; - - Axis_t rmax = fChargeLastHiGain; - if (TESTBIT(fFlags,kUseLoGain)) - rmax = fChargeLastLoGain; - - TH1F *hist = fHChargeHiGain; - if (TESTBIT(fFlags,kUseLoGain)) - hist = fHChargeLoGain; - + { + rmin = fChargeFirstLoGain; + rmax = fChargeLastLoGain; + hist = fHChargeLoGain; + } + // // First guesses for the fit (should be as close to reality as possible, - // otherwise the fit goes gaga because of high number of dimensions ... - // - const Stat_t entries = hist->Integral(); - const Double_t area_guess = entries/gkSq2Pi; + // + const Stat_t entries = hist->Integral("width"); const Double_t mu_guess = hist->GetBinCenter(hist->GetMaximumBin()); - const Double_t sigma_guess = mu_guess/15.; - - TString name = TString("ChargeGausFit"); - name += fPixId; - - fChargeGausFit = new TF1(name.Data(),"gaus",rmin,rmax); + const Double_t sigma_guess = (rmax-rmin)/2.; + const Double_t area_guess = entries/gkSq2Pi/sigma_guess; + + fChargeGausFit = new TF1(Form("%s%d","ChargeGausFit",fPixId),"gaus",rmin,rmax); if (!fChargeGausFit) @@ -831,13 +833,13 @@ // If we are not able to fit, try once again // - if (fChargeGausFit->GetProb() < gkProbLimit) - { - - Axis_t rtry = fChargeGausFit->GetParameter(1) - 3.0*fChargeGausFit->GetParameter(2); + if (fChargeGausFit->GetProb() < fProbLimit) + { + + Axis_t rtry = fChargeGausFit->GetParameter(1) - 2.5*fChargeGausFit->GetParameter(2); rmin = (rtry < rmin ? rmin : rtry); - rmax = fChargeGausFit->GetParameter(1) + 3.0*fChargeGausFit->GetParameter(2); + rmax = fChargeGausFit->GetParameter(1) + 2.5*fChargeGausFit->GetParameter(2); + fChargeGausFit->SetRange(rmin,rmax); - - fHChargeHiGain->Fit(fChargeGausFit,option); + hist->Fit(fChargeGausFit,option); } @@ -852,19 +854,28 @@ // // The fit result is accepted under condition: - // The Results are not nan's - // The Probability is greater than gkProbLimit (default 0.001 == 99.9%) - // - if (TMath::IsNaN(fChargeMean) || TMath::IsNaN(fChargeMeanErr)) - { + // 1) The results are not nan's + // 2) The NDF is not smaller than fNDFLimit (5) + // 3) The Probability is greater than fProbLimit (default 0.001 == 99.9%) + // + // Otherwise take means and RMS of the histograms + // + if ( TMath::IsNaN(fChargeMean) + || TMath::IsNaN(fChargeMeanErr) + || TMath::IsNaN(fChargeProb) + || TMath::IsNaN(fChargeSigma) + || TMath::IsNaN(fChargeSigmaErr) + || (fChargeNdf < fNDFLimit) + || (fChargeProb < fProbLimit) ) + { + + fChargeMean = hist->GetMean(); + fChargeMeanErr = hist->GetRMS()/hist->GetEntries(); + fChargeSigma = hist->GetRMS(); + fChargeSigmaErr = fChargeMeanErr/2.; + CLRBIT(fFlags,kChargeFitOK); return kFALSE; } - if ((fChargeProb < gkProbLimit) || (TMath::IsNaN(fChargeProb))) - { - CLRBIT(fFlags,kChargeFitOK); - return kFALSE; - } - SETBIT(fFlags,kChargeFitOK); return kTRUE; @@ -877,5 +888,5 @@ { - Int_t nbins = 30; + Int_t nbins = 50; CutEdges(fHChargeHiGain,nbins); @@ -891,5 +902,5 @@ +fHChargeLoGain->GetBinWidth(0); - CutEdges(fHRelTimeHiGain,nbins); + CutEdges(fHRelTimeHiGain,0); fRelTimeLowerFitRangeHiGain = fHRelTimeHiGain->GetBinLowEdge(fHRelTimeHiGain->GetXaxis()->GetFirst()); @@ -897,5 +908,5 @@ +fHRelTimeHiGain->GetBinWidth(0); - CutEdges(fHRelTimeLoGain,nbins); + CutEdges(fHRelTimeLoGain,0); fRelTimeLowerFitRangeLoGain = fHRelTimeLoGain->GetBinLowEdge(fHRelTimeLoGain->GetXaxis()->GetFirst()); @@ -903,5 +914,5 @@ +fHRelTimeLoGain->GetBinWidth(0); - CutEdges(fHAbsTimeHiGain,nbins); + CutEdges(fHAbsTimeHiGain,0); fAbsTimeFirstHiGain = fHAbsTimeHiGain->GetBinLowEdge(fHAbsTimeHiGain->GetXaxis()->GetFirst()); @@ -909,5 +920,5 @@ +fHAbsTimeHiGain->GetBinWidth(0); - CutEdges(fHAbsTimeLoGain,nbins); + CutEdges(fHAbsTimeLoGain,0); fAbsTimeFirstLoGain = fHAbsTimeLoGain->GetBinLowEdge(fHAbsTimeLoGain->GetXaxis()->GetFirst()); Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationPixel.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationPixel.h (revision 3024) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationPixel.h (revision 3025) @@ -20,12 +20,13 @@ static const Int_t fChargeNbinsHiGain; static const Int_t fChargeNbinsLoGain; + static const Int_t fChargevsNbins; static const Int_t fAbsTimeNbins; - static const Int_t fRelTimeNbins; - static const Int_t fChargevsNbins; static const Axis_t fAbsTimeFirst; static const Axis_t fAbsTimeLast; + static const Int_t fRelTimeNbins; static const Axis_t fRelTimeFirst; static const Axis_t fRelTimeLast; static const Float_t fProbLimit; + static const Int_t fNDFLimit; Int_t fPixId; // Pixel Nr