Index: trunk/MagicSoft/Mars/mjobs/MJCalibration.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJCalibration.cc (revision 3054) +++ trunk/MagicSoft/Mars/mjobs/MJCalibration.cc (revision 3055) @@ -52,4 +52,5 @@ #include "MGeomApply.h" #include "MExtractSignal.h" +#include "MExtractSignal2.h" #include "MCalibrationCalc.h" @@ -58,5 +59,4 @@ ClassImp(MJCalibration); - using namespace std; @@ -69,83 +69,100 @@ void MJCalibration::DrawProjection(MHCamera *obj1, Int_t fit) const { - TH1D *obj2 = (TH1D*)obj1->Projection(); - obj2->Draw(); - obj2->SetBit(kCanDelete); - - const Double_t min = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst()); - const Double_t max = obj2->GetBinCenter(obj2->GetXaxis()->GetLast()); - const Double_t integ = obj2->Integral("width")/2.5; - const Double_t mean = obj2->GetMean(); - const Double_t rms = obj2->GetRMS(); - const Double_t width = max-min; - - TF1 *f=0; - switch (fit) - { - case 0: - f = new TF1("sgaus", "gaus(0)", min, max); - f->SetLineColor(kYellow); - f->SetBit(kCanDelete); - f->SetParNames("Area", "#mu", "#sigma"); - f->SetParameters(integ/rms, mean, rms); - f->SetParLimits(0, 0, integ); - f->SetParLimits(1, min, max); - f->SetParLimits(2, 0, width/1.5); - - obj2->Fit(f, "QLR"); - break; - + + TH1D *obj2 = (TH1D*)obj1->Projection(); + obj2->Draw(); + obj2->SetBit(kCanDelete); + + const Double_t min = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst()); + const Double_t max = obj2->GetBinCenter(obj2->GetXaxis()->GetLast()); + const Double_t integ = obj2->Integral("width")/2.5; + const Double_t mean = obj2->GetMean(); + const Double_t rms = obj2->GetRMS(); + const Double_t width = max-min; + + const TString dgausformula = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])" + "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])"; + + const TString tgausformula = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])" + "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])" + "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])"; + TF1 *f=0; + switch (fit) + { case 1: - f = new TF1("dgaus", "gaus(0)+gaus(3)", min, max); - f->SetLineColor(kYellow); - f->SetBit(kCanDelete); - f->SetParNames("A1", "#mu1", "#sigma1", "A2", "#mu2", "#sigma2"); - f->SetParameters(integ/width, max-width/6, width/4, integ/width, min+width/6, width/4); - f->SetParLimits(0, 0, integ); - f->SetParLimits(1, min, max); - f->SetParLimits(2, 0, width/3); - f->SetParLimits(3, 0, integ); - f->SetParLimits(4, min, max); - f->SetParLimits(5, 0, width/3); - - obj2->Fit(f, "QLR"); - break; - + f = new TF1("sgaus", "gaus(0)", min, max); + f->SetLineColor(kYellow); + f->SetBit(kCanDelete); + f->SetParNames("Area", "#mu", "#sigma"); + f->SetParameters(integ/rms, mean, rms); + f->SetParLimits(0, 0, integ); + f->SetParLimits(1, min, max); + f->SetParLimits(2, 0, width/1.5); + + obj2->Fit(f, "QLR"); + break; + case 2: - f = new TF1("tgaus", "gaus(0)+gaus(3)+gaus(6)", min, max); - f->SetLineColor(kYellow); - f->SetBit(kCanDelete); - f->SetParNames("A1", "#mu1", "#sigma1", "A2", "#mu2", "#sigma2", "A3", "#mu3", "#sigma3"); - f->SetParameters(integ/width, max-width/6, width/4, integ/width, min+width/6, width/4, integ/width,min+width/6, width/2); - f->SetParLimits(0, 0, integ); - f->SetParLimits(1, min, max); - f->SetParLimits(2, 0, width/4); - f->SetParLimits(3, 0, integ); - f->SetParLimits(4, min, max); - f->SetParLimits(5, 0, width/4); - f->SetParLimits(6, 0, integ); - f->SetParLimits(7, min, max); - f->SetParLimits(8, 0, width/2); - - obj2->Fit(f, "QLR"); - break; - + f = new TF1("dgaus",dgausformula.Data(),min,max); + f->SetLineColor(kYellow); + f->SetBit(kCanDelete); + f->SetParNames("A_{tot}", "#mu1", "#sigma1", "A2", "#mu2", "#sigma2"); + f->SetParameters(integ,(min+mean)/2.,width/4., + integ/width/2.,(max+mean)/2.,width/4.); + // The left-sided Gauss + f->SetParLimits(0,integ-1.5 , integ+1.5); + f->SetParLimits(1,min+(width/10.), mean); + f->SetParLimits(2,0 , width/2.); + // The right-sided Gauss + f->SetParLimits(3,0 , integ); + f->SetParLimits(4,mean, max-(width/10.)); + f->SetParLimits(5,0 , width/2.); + obj2->Fit(f,"QLRM"); + break; + case 3: - obj2->Fit("pol0", "Q"); - obj2->GetFunction("pol0")->SetLineColor(kYellow); - break; - + f = new TF1("tgaus",tgausformula.Data(),min,max); + f->SetLineColor(kYellow); + f->SetBit(kCanDelete); + f->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}", + "A_{2}","#mu_{2}","#sigma_{2}", + "A_{3}","#mu_{3}","#sigma_{3}"); + f->SetParameters(integ,(min+mean)/2,width/4., + integ/width/3.,(max+mean)/2.,width/4., + integ/width/3.,mean,width/2.); + // The left-sided Gauss + f->SetParLimits(0,integ-1.5,integ+1.5); + f->SetParLimits(1,min+(width/10.),mean); + f->SetParLimits(2,width/15.,width/2.); + // The right-sided Gauss + f->SetParLimits(3,0.,integ); + f->SetParLimits(4,mean,max-(width/10.)); + f->SetParLimits(5,width/15.,width/2.); + // The Gauss describing the outliers + f->SetParLimits(6,0.,integ); + f->SetParLimits(7,min,max); + f->SetParLimits(8,width/4.,width/1.5); + obj2->Fit(f,"QLRM"); + break; + + case 4: + obj2->Fit("pol0", "Q"); + obj2->GetFunction("pol0")->SetLineColor(kYellow); + break; + case 9: - break; - + break; + default: - obj2->Fit("gaus", "Q"); - obj2->GetFunction("gaus")->SetLineColor(kYellow); - break; - } -} - -void MJCalibration::CamDraw(TCanvas &c, Int_t x, Int_t y, MHCamera &cam1, Int_t fit) -{ + obj2->Fit("gaus", "Q"); + obj2->GetFunction("gaus")->SetLineColor(kYellow); + break; + } +} + +void MJCalibration::CamDraw(TCanvas &c, const Int_t x, const Int_t y, const MHCamera &cam1, const Int_t fit) +{ + + c.cd(x); gPad->SetBorderMode(0); @@ -157,16 +174,12 @@ obj1->Draw(); - c.cd(x+2*y); - gPad->SetBorderMode(0); - DrawProjection(obj1, fit); - - // tbretz: Using gStyle in this context is totally useless! - //const Float_t he = gStyle->GetStatH(); - //const Float_t wi = gStyle->GetStatH(); - //gStyle->SetStatH(0.4); - //gStyle->SetStatW(0.25); - //gStyle->SetStatH(he); - //gStyle->SetStatW(wi); -} + if (fit) + { + c.cd(x+2*y); + gPad->SetBorderMode(0); + DrawProjection(obj1, fit); + } +} + void MJCalibration::DisplayResult(MParList &plist) @@ -190,118 +203,158 @@ // Create histograms to display - MHCamera disp1 (geomcam, "Cal;Charge", "Fitted Mean Charges"); - MHCamera disp3 (geomcam, "Cal;SigmaCharge", "Sigma of Fitted Charges"); - MHCamera disp5 (geomcam, "Cal;FitProb", "Probability of Fit"); - /* - MHCamera disp6 (geomcam, "Cal;Time", "Arrival Times"); - MHCamera disp7 (geomcam, "Cal;SigmaTime", "Sigma of Arrival Times"); - MHCamera disp8 (geomcam, "Cal;TimeChiSquare", "Chi Square of Time Fit"); - */ -// MHCamera disp9 (geomcam, "Cal;Ped", "Pedestals"); -// MHCamera disp10(geomcam, "Cal;PedRms", "Pedestal RMS"); -// MHCamera disp11(geomcam, "Cal;RSigma", "Reduced Sigmas"); - MHCamera disp12(geomcam, "Cal;FFactorPhe", "Nr. of Phe's (F-Factor Method)"); - MHCamera disp13(geomcam, "Cal;FFactorConv", "Conversion Factor (F-Factor Method)"); - MHCamera disp14(geomcam, "Cal;BlindPixPhe", "Nr. of Photons inside plexiglass (Blind Pixel Method)"); - MHCamera disp15(geomcam, "Cal;BlindPixConv", "Conversion Factor (Blind Pixel Method)"); -// MHCamera disp16(geomcam, "Cal;RSigma/Charge", "Reduced Sigma per Charge"); - - disp1.SetCamContent(fCalibrationCam, 0); - disp1.SetCamError(fCalibrationCam, 1); - - disp3.SetCamContent(fCalibrationCam, 2); - disp3.SetCamError(fCalibrationCam, 3); - - disp5.SetCamContent(fCalibrationCam, 4); - - /* - disp6.SetCamContent(fCalibrationCam, 5); - disp6.SetCamError(fCalibrationCam, 6); - disp7.SetCamContent(fCalibrationCam, 6); - disp8.SetCamContent(fCalibrationCam, 7); - */ -/* - disp9.SetCamContent(calcam, 8); - disp9.SetCamError(calcam, 9); - disp10.SetCamContent(calcam, 9); - disp11.SetCamContent(fCalibrationCam, 10); - */ - disp12.SetCamContent(fCalibrationCam, 11); - disp12.SetCamError(fCalibrationCam, 12); - disp13.SetCamContent(fCalibrationCam, 13); - disp13.SetCamError(fCalibrationCam, 14); - - disp14.SetCamContent(fCalibrationCam, 15); - disp15.SetCamContent(fCalibrationCam, 16); - // disp16.SetCamContent(fCalibrationCam, 17); + MHCamera disp1 (geomcam, "Cal;Charge", "Fitted Mean Charges"); + MHCamera disp2 (geomcam, "Cal;SigmaCharge", "Sigma of Fitted Charges"); + MHCamera disp3 (geomcam, "Cal;FitProb", "Probability of Fit"); + MHCamera disp4 (geomcam, "Cal;RSigma", "Reduced Sigmas"); + MHCamera disp5 (geomcam, "Cal;RSigma/Charge", "Reduced Sigma per Charge"); + MHCamera disp6 (geomcam, "Cal;FFactorPhe", "Nr. of Phe's (F-Factor Method)"); + MHCamera disp7 (geomcam, "Cal;FFactorConv", "Conversion Factor (F-Factor Method)"); + MHCamera disp8 (geomcam, "Cal;FFactorFFactor", "Total F-Factor (F-Factor Method)"); + MHCamera disp9 (geomcam, "Cal;BlindPixPh", "Nr. of Photons inside plexiglass (Blind Pixel Method)"); + MHCamera disp10 (geomcam, "Cal;BlindPixConv", "Conversion Factor (Blind Pixel Method)"); + MHCamera disp11 (geomcam, "Cal;BlindPixFFactor","Total F-Factor (Blind Pixel Method)"); + MHCamera disp12 (geomcam, "Cal;PINDiodePh", "Nr. of Photons outside plexiglass (PIN Diode Method)"); + MHCamera disp13 (geomcam, "Cal;PINDiodeConv", "Conversion Factor (PIN Diode Method)"); + MHCamera disp14 (geomcam, "Cal;PINDiodeFFactor","Total F-Factor (PIN Diode Method)"); + MHCamera disp15 (geomcam, "Cal;Excluded", "Pixels previously excluded"); + MHCamera disp16 (geomcam, "Cal;NotFitted", "Pixels that could not be fitted"); + MHCamera disp17 (geomcam, "Cal;NotFitValid", "Pixels with not valid fit results"); + MHCamera disp18 (geomcam, "Cal;Oscillating", "Oscillating Pixels"); + MHCamera disp19 (geomcam, "Cal;Saturation", "Pixels with saturated Hi Gain"); + + + // Fitted charge means and sigmas + disp1.SetCamContent(fCalibrationCam, 0); + disp1.SetCamError( fCalibrationCam, 1); + disp2.SetCamContent(fCalibrationCam, 2); + disp2.SetCamError( fCalibrationCam, 3); + // Fit probabilities + disp3.SetCamContent(fCalibrationCam, 4); + + // Reduced Sigmas and reduced sigmas per charge + disp4.SetCamContent(fCalibrationCam, 5); + disp4.SetCamError( fCalibrationCam, 6); + disp5.SetCamContent(fCalibrationCam, 7); + disp5.SetCamError( fCalibrationCam, 8); + + // F-Factor Method + disp6.SetCamContent(fCalibrationCam, 9); + disp6.SetCamError( fCalibrationCam, 10); + disp7.SetCamContent(fCalibrationCam, 11); + disp7.SetCamError( fCalibrationCam, 12); + disp8.SetCamContent(fCalibrationCam, 13); + disp8.SetCamError( fCalibrationCam, 14); + + /// Blind Pixel Method + disp9.SetCamContent(fCalibrationCam, 15); + disp9.SetCamError( fCalibrationCam, 16); + disp10.SetCamContent(fCalibrationCam,17); + disp10.SetCamError( fCalibrationCam,18); + disp11.SetCamContent(fCalibrationCam,19); + disp11.SetCamError( fCalibrationCam,20); + + // PIN Diode Method + disp12.SetCamContent(fCalibrationCam,21); + disp12.SetCamError( fCalibrationCam,22); + disp13.SetCamContent(fCalibrationCam,23); + disp13.SetCamError( fCalibrationCam,24); + disp14.SetCamContent(fCalibrationCam,25); + disp14.SetCamError( fCalibrationCam,26); + + // Pixels with defects + disp15.SetCamContent(fCalibrationCam,27); + disp16.SetCamContent(fCalibrationCam,28); + disp17.SetCamContent(fCalibrationCam,29); + disp18.SetCamContent(fCalibrationCam,30); + + // Lo Gain calibration + disp19.SetCamContent(fCalibrationCam,31); + disp1.SetYTitle("Charge [FADC units]"); - disp3.SetYTitle("\\sigma_{Charge} [FADC units]"); - disp5.SetYTitle("P_{Charge} [1]"); - /* - disp6.SetYTitle("Arr. Time [FADC slice Nr.]"); - disp7.SetYTitle("\\sigma_{Time} [FADC slices]"); - disp8.SetYTitle("\\chi^{2}_{Time} [1]"); - disp9.SetYTitle("Ped [FADC Counts ]"); - disp10.SetYTitle("RMS_{Ped} [FADC Counts ]"); - disp11.SetYTitle("\\sqrt{\\sigma^{2}_{Charge} - RMS^{2}_{Ped}} [FADC units]"); - */ - disp12.SetYTitle("Nr. Photo-Electrons [1]"); - disp13.SetYTitle("Conversion Factor [PhE/FADC Count]"); - disp14.SetYTitle("Nr. Photons [1]"); - disp15.SetYTitle("Conversion Factor [Phot/FADC Count]"); -// disp16.SetYTitle("Reduced Sigma / Charge [1]"); + disp2.SetYTitle("\\sigma_{Charge} [FADC units]"); + disp3.SetYTitle("P_{Charge} [1]"); + + disp4.SetYTitle("\\sqrt{\\sigma^{2}_{Charge} - RMS^{2}_{Ped}} [FADC Counts]"); + disp5.SetYTitle("Reduced Sigma / Mean Charge [1]"); + + disp6.SetYTitle("Nr. Photo-Electrons [1]"); + disp7.SetYTitle("Conversion Factor [PhE/FADC Count]"); + disp8.SetYTitle("\\sqrt{N_{PhE}}*\\sigma_{Charge}/\\mu_{Charge} [1]"); + + disp9.SetYTitle("Nr. Photons [1]"); + disp10.SetYTitle("Conversion Factor [Phot/FADC Count]"); + disp11.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]"); + + disp12.SetYTitle("Nr. Photons [1]"); + disp13.SetYTitle("Conversion Factor [Phot/FADC Count]"); + disp14.SetYTitle("\\sqrt{N_{Ph}}*\\sigma_{Charge}/\\mu_{Charge} [1]"); + + disp15.SetYTitle("[1]"); + disp16.SetYTitle("[1]"); + disp17.SetYTitle("[1]"); + disp18.SetYTitle("[1]"); gStyle->SetOptStat(1111); gStyle->SetOptFit(); - TCanvas &c1 = fDisplay->AddTab("FitCharge"); + // Charges + TCanvas &c1 = fDisplay->AddTab("Fit.Charge"); c1.Divide(2, 3); - CamDraw(c1, 1, 2, disp1, 1); - CamDraw(c1, 2, 2, disp3, 2); + CamDraw(c1, 1, 2, disp1, 2); + CamDraw(c1, 2, 2, disp2, 2); // Fit Probability - TCanvas &c2 = fDisplay->AddTab("FitProb"); + TCanvas &c2 = fDisplay->AddTab("Fit.Prob"); c2.Divide(1,3); - CamDraw(c2, 1, 1, disp5, 3); -/* - // Times - TCanvas &c3 = fDisplay->AddTab("FitTimes"); - c3.Divide(3,3); - - CamDraw(c3, 1, 3, disp6, 1); - CamDraw(c3, 2, 3, disp7, 0); - CamDraw(c3, 3, 3, disp8, 0); - - // Pedestals - TCanvas &c4 = d3->AddTab("Pedestals"); - c4.Divide(2,3); - - CamDraw(c4, 1, 2, disp9, 0); - CamDraw(c4, 2, 2, disp10, 1); + CamDraw(c2, 1, 1, disp3, 4); // Reduced Sigmas - TCanvas &c5 = fDisplay->AddTab("RedSigma"); - c5.Divide(2,3); - - CamDraw(c5, 1, 2, disp11, 2); - CamDraw(c5, 2, 2, disp16, 2); - */ + TCanvas &c3 = fDisplay->AddTab("Red.Sigma"); + c3.Divide(2,3); + + CamDraw(c3, 1, 2, disp4, 2); + CamDraw(c3, 2, 2, disp5, 2); // F-Factor Method - TCanvas &c6 = fDisplay->AddTab("F-Factor"); - c6.Divide(2,3); - - CamDraw(c6, 1, 2, disp12, 1); - CamDraw(c6, 2, 2, disp13, 1); + TCanvas &c4 = fDisplay->AddTab("F-Factor"); + c4.Divide(3,3); + + CamDraw(c4, 1, 3, disp6, 2); + CamDraw(c4, 2, 3, disp7, 2); + CamDraw(c4, 3, 3, disp8, 2); // Blind Pixel Method - TCanvas &c7 = fDisplay->AddTab("BlindPix"); - c7.Divide(2, 3); - - CamDraw(c7, 1, 2, disp14, 9); - CamDraw(c7, 2, 2, disp15, 1); + TCanvas &c5 = fDisplay->AddTab("BlindPix"); + c5.Divide(3, 3); + + CamDraw(c5, 1, 3, disp9, 9); + CamDraw(c5, 2, 3, disp10, 2); + CamDraw(c5, 3, 3, disp11, 2); + + // PIN Diode Method + TCanvas &c6 = fDisplay->AddTab("PINDiode"); + c6.Divide(3,3); + + CamDraw(c6, 1, 3, disp12, 9); + CamDraw(c6, 2, 3, disp13, 2); + CamDraw(c6, 3, 3, disp14, 2); + + // Defects + TCanvas &c7 = fDisplay->AddTab("Defects"); + c7.Divide(4,2); + + CamDraw(c7, 1, 4, disp15, 0); + CamDraw(c7, 2, 4, disp16, 0); + CamDraw(c7, 3, 4, disp17, 0); + CamDraw(c7, 4, 4, disp18, 0); + + // Lo Gain Calibration + TCanvas &c8 = fDisplay->AddTab("LowGain"); + c8.Divide(1,3); + + CamDraw(c8, 1, 4, disp19, 0); } @@ -413,5 +466,5 @@ MGeomApply apply; - MExtractSignal extract; + MExtractSignal2 extract; MCalibrationCalc calcalc;