Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 2789) +++ trunk/MagicSoft/Mars/Changelog (revision 2790) @@ -9,4 +9,12 @@ - make Stat Window in option 'proj' bigger - add comment line at beginning of Draw() + + * mcalib/MHCalibrationBlindPixel.[h.cc] + - add Function DrawClone() + - fix bug in calculation of integral (GetIntegral("width")) + - fit now with fixed normalization + - restrict parameters tighter + - declare fit as not valid, if less than 100 events + in single photo-electron peak Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.cc (revision 2789) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.cc (revision 2790) @@ -27,6 +27,10 @@ // MHCalibrationBlindPixel // // // -// Performs all the necessary fits to extract the mean number of photons // -// out of the derived light flux // +// Performs all the Single Photo-Electron Fit to extract // +// the mean number of photons and to derive the light flux // +// // +// The fit result is accepted under condition that: // +// 1) the Probability is greater than gkProbLimit (default 0.001 == 99.7%) // +// 2) at least 100 events are in the single Photo-electron peak // // // ////////////////////////////////////////////////////////////////////////////// @@ -83,4 +87,5 @@ fHBlindPixelCharge->SetYTitle("Nr. of events"); fHBlindPixelCharge->Sumw2(); + fHBlindPixelCharge->SetDirectory(NULL); Axis_t tfirst = -0.5; @@ -92,4 +97,5 @@ fHBlindPixelTime->SetYTitle("Nr. of events"); fHBlindPixelTime->Sumw2(); + fHBlindPixelTime->SetDirectory(NULL); // We define a reasonable number and later enlarge it if necessary @@ -101,6 +107,7 @@ fHBlindPixelChargevsN->SetXTitle("Event Nr."); fHBlindPixelChargevsN->SetYTitle("Sum of FADC slices"); - - fgSinglePheFitFunc = &gfKto8; + fHBlindPixelChargevsN->SetDirectory(NULL); + + fgSinglePheFitFunc = &gfKto4; fgSinglePheFitNPar = 5; } @@ -187,4 +194,26 @@ } + + +TObject *MHCalibrationBlindPixel::DrawClone(Option_t *opt) const +{ + + //TVirtualPad *pad = gROOT->GetSelectedPad(); + TVirtualPad *padsav = gPad; + //if (pad) pad->cd(); + + TObject *newobj = Clone(); + + if (!newobj) + return 0; + + newobj->Draw(opt); + + if (padsav) + padsav->cd(); + + return newobj; +} + @@ -317,48 +346,35 @@ // otherwise the fit goes gaga because of high number of dimensions ... // - const Stat_t entries = fHBlindPixelCharge->Integral(); + const Stat_t entries = fHBlindPixelCharge->Integral("width"); const Double_t lambda_guess = 0.5; const Double_t mu_0_guess = fHBlindPixelCharge->GetBinCenter(fHBlindPixelCharge->GetMaximumBin()); const Double_t si_0_guess = 20.; - const Double_t mu_1_guess = mu_0_guess + 100.; + const Double_t mu_1_guess = mu_0_guess + 50.; const Double_t si_1_guess = si_0_guess + si_0_guess; fSinglePheFit = new TF1("SinglePheFit",fgSinglePheFitFunc,rmin,rmax,fgSinglePheFitNPar+1); + // fSinglePheFit = new TF1("SinglePheFit",fgSinglePheFitFunc,rmin,rmax,fgSinglePheFitNPar+1); + // fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess); fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,entries); + // fSinglePheFit->SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1"); fSinglePheFit->SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1","area"); - fSinglePheFit->SetParLimits(0,0.,3.); - fSinglePheFit->SetParLimits(1,rmin,rmax); - fSinglePheFit->SetParLimits(2,rmin,rmax); - fSinglePheFit->SetParLimits(3,1.0,rmax-rmin); - fSinglePheFit->SetParLimits(4,1.7,rmax-rmin); - fSinglePheFit->SetParLimits(5,0.,1.5*entries); + fSinglePheFit->SetParLimits(0,0.,1.); + fSinglePheFit->SetParLimits(1,rmin,(rmax-rmin)/1.5); + fSinglePheFit->SetParLimits(2,(rmax-rmin)/2.,(rmax-0.05*(rmax-rmin))); + fSinglePheFit->SetParLimits(3,1.0,(rmax-rmin)/2.0); + fSinglePheFit->SetParLimits(4,1.0,(rmax-rmin)/2.5); + // fSinglePheFit->SetParLimits(5,entries/gkSq2Pi,entries/gkSq2Pi); + // fSinglePheFit->SetParLimits(5,0.,1.5*entries); // // Normalize the histogram to facilitate faster fitting of the area - // For speed reasons, FKto8 is normalized to Sqrt(2 pi). + // For speed reasons, gfKto4 is normalized to Sqrt(2 pi). // Therefore also normalize the histogram to that value // - // ROOT gives us another nice example of user-unfriendly behavior: - // Although the normalization of the function fSinglePhe and the - // Histogram fHBlindPixelCharge agree (!!), the fit does not normalize correctly INTERNALLY - // in the fitting procedure !!! - // - // This has to do with the fact that the internal function histogramming - // uses 100 bins and does not adapt to the binning of the fitted histogram, unlike PAW does - // (very important if you use Sumw2, see e.g. ROOTTALK: Mon May 26 1997 - 09:56:03 MEST) - // - // So, WE have to adapt to that internal flaw of ROOT: - // - // const Int_t npx = fSinglePheFit->GetNpx(); - // const Int_t bins = fHBlindPixelCharge->GetXaxis()->GetLast()-fHBlindPixelCharge->GetXaxis()->GetFirst(); // fHBlindPixelCharge->Scale(gkSq2Pi*(float)bins/npx/entries); - - // - // we need this, otherwise, ROOT does not calculate the area correctly - // don't ask me why it does not behave correctly, it's one of the nasty - // mysteries of ROOT which takes you a whole day to find out :-) - // - // fSinglePheFit->SetNpx(fChargenbins); - - fHBlindPixelCharge->Fit(fSinglePheFit,opt); + // fHBlindPixelCharge->Scale(gkSq2Pi/entries); + Float_t norm = entries/gkSq2Pi; + // norm *= (Float_t)fSinglePheFit->GetNpx()/(Float_t)fBlindPixelChargenbins; + fSinglePheFit->SetParLimits(5,norm-0.1,norm+0.1); + fHBlindPixelCharge->Fit(fSinglePheFit,opt); @@ -384,9 +400,5 @@ Double_t pedarea = fSinglePhePedFit->GetParameter(0)*gkSq2Pi*fSinglePhePedFit->GetParameter(2); - cout << "Parameter0: " << fSinglePhePedFit->GetParameter(0) << endl; - cout << "Parameter2: " << fSinglePhePedFit->GetParameter(2) << endl; - cout << "Pedarea: " << pedarea << endl; - cout << "entries: " << entries << endl; - fLambdaCheck = TMath::Log((double)entries/pedarea); + fLambdaCheck = TMath::Log(entries/pedarea); fLambdaCheckErr = fSinglePhePedFit->GetParError(0)/fSinglePhePedFit->GetParameter(0) + fSinglePhePedFit->GetParError(2)/fSinglePhePedFit->GetParameter(2); @@ -396,18 +408,29 @@ *fLog << "DoF: " << fNdf << endl; *fLog << "Probability: " << fProb << endl; - *fLog << "Integral: " << fSinglePheFit->Integral(rmin,rmax); - - // - // The fit result is accepted under condition - // The fit result is accepted under condition - // The Probability is greater than gkProbLimit (default 0.001 == 99.7%) + + // + // The fit result is accepted under condition that: + // 1) the Probability is greater than gkProbLimit (default 0.001 == 99.7%) + // 2) at least 100 events are in the single Photo-electron peak // if (fProb < gkProbLimit) { - *fLog << warn << "Prob: " << fProb << " is smaller than the allowed value: " << gkProbLimit << endl; + *fLog << err << "Prob: " << fProb << " is smaller than the allowed value: " << gkProbLimit << endl; fFitOK = kFALSE; return kFALSE; } - + + Float_t contSinglePhe = TMath::Exp(-1.0*fLambda)*fLambda*entries; + + if (contSinglePhe < 100.) + { + *fLog << err << "Statistics is too low: Only " << contSinglePhe + << " in the first electron peak " << endl; + fFitOK = kFALSE; + return kFALSE; + } + else + *fLog << GetDescriptor() << ": " << contSinglePhe + << " in first electron peak " << endl; fFitOK = kTRUE; @@ -420,5 +443,5 @@ { - Int_t nbins = 50; + Int_t nbins = 30; CutEdges(fHBlindPixelCharge,nbins); Index: trunk/MagicSoft/Mars/mhist/MHCamera.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHCamera.cc (revision 2789) +++ trunk/MagicSoft/Mars/mhist/MHCamera.cc (revision 2790) @@ -375,4 +375,5 @@ // 'content' Display the relative content aligned to GetMaximum() and // GeMinimum() ((val-min)/(max-min)) +// 'proj' Display the y-projection of the histogram // void MHCamera::Draw(Option_t *option) @@ -620,5 +621,11 @@ CreateProjection(); opt.ReplaceAll("proj", ""); + Float_t he = gStyle->GetStatH(); + Float_t wi = gStyle->GetStatH(); + gStyle->SetStatH(0.4); + gStyle->SetStatW(0.25); fYProj->Paint(opt); + gStyle->SetStatH(he); + gStyle->SetStatW(wi); return; }