Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 3058) +++ /trunk/MagicSoft/Mars/Changelog (revision 3059) @@ -48,4 +48,7 @@ - removed times of GetPixelContent - added OscillationChecks in GetPixelContent + + * manalysis/MMcCalibrationUpdate.cc + - remove calpix.SkipTimeFits which is now obsolet Index: /trunk/MagicSoft/Mars/macros/calibration.C =================================================================== --- /trunk/MagicSoft/Mars/macros/calibration.C (revision 3058) +++ /trunk/MagicSoft/Mars/macros/calibration.C (revision 3059) @@ -29,5 +29,4 @@ const TString calfile = "/mnt/Data/rootdata/CrabNebula/2004_01_27/20040126_12525_C_Cab-On_E.root"; - //const TString pedfile = "/mnt/Data/rootdata/Miscellaneous/2003_12_19/20031218_03522_P_Park_E.root"; //const TString calfile = "/mnt/Data/rootdata/Miscellaneous/2003_12_19/20031218_03527_C_Park_E.root"; @@ -36,4 +35,8 @@ TString calname=calfile) { + + MStatusDisplay *display = new MStatusDisplay; + display->SetUpdateTime(3000); + display->Resize(850,700); // @@ -63,4 +66,5 @@ MPedCalcPedRun pedcalc; + pedcalc.SetNumHiGainSamples(14); MFillH fill("MPedestalCam", "MExtractedSignalCam"); @@ -81,4 +85,5 @@ MEvtLoop evtloop; evtloop.SetParList(&plist); + evtloop.SetDisplay(display); // @@ -92,21 +97,18 @@ // pedcam(17).DrawClone(); - MStatusDisplay *d1 = new MStatusDisplay; - d1->SetUpdateTime(3000); - d1->Resize(850,700); - - MHCamera dispped0 (geomcam, "MPedestalPix;Pedestal", "Mean per Slice"); - MHCamera dispped1 (geomcam, "MPedestalPix;PedestalErr", "Mean Error per Slice"); - MHCamera dispped2 (geomcam, "MPedestalPix;PedestalRms", "RMS per Slice"); - MHCamera dispped3 (geomcam, "MPedestalPix;PedestalRmsErr", "RMS Error per Slice"); - MHCamera dispped4 (geomcam, "MPedestalPix;Mean", "Fitted Mean per Slice"); - MHCamera dispped5 (geomcam, "MPedestalPix;MeanErr", "Fitted Error of Mean per Slice"); - MHCamera dispped6 (geomcam, "MPedestalPix;Sigma", "Fitted Sigma per Slice"); - MHCamera dispped7 (geomcam, "MPedestalPix;SigmaErr", "Fitted Error of Sigma per Slice"); - MHCamera dispped8 (geomcam, "MPedestalPix;Prob", "Probability of Fit"); - MHCamera dispped9 (geomcam, "MPedestalPix;DeltaPedestalMean", "Rel. Diff. Mean per Slice (Calc.-Fitte)"); - MHCamera dispped11 (geomcam, "MPedestalPix;DeltaPedestalMeanError", "Rel. Diff. Mean Error per Slice (Calc.-Fitted)"); - MHCamera dispped12 (geomcam, "MPedestalPix;DeltaRmsSigma", "Rel. Diff. RMS per Slice (Calc.-Fitted)"); - MHCamera dispped14 (geomcam, "MPedestalPix;DeltaRmsSigmaError", "Rel. Diff. RMS Error per Slice (Calc.-Fitted)"); + MHCamera dispped0 (geomcam, "Ped;Pedestal", "Mean per Slice"); + MHCamera dispped1 (geomcam, "Ped;PedestalErr", "Mean Error per Slice"); + MHCamera dispped2 (geomcam, "Ped;PedestalRms", "RMS per Slice"); + MHCamera dispped3 (geomcam, "Ped;PedestalRmsErr", "RMS Error per Slice"); + MHCamera dispped4 (geomcam, "Ped;Mean", "Fitted Mean per Slice"); + MHCamera dispped5 (geomcam, "Ped;MeanErr", "Fitted Error of Mean per Slice"); + MHCamera dispped6 (geomcam, "Ped;Sigma", "Fitted Sigma per Slice"); + MHCamera dispped7 (geomcam, "Ped;SigmaErr", "Fitted Error of Sigma per Slice"); + MHCamera dispped8 (geomcam, "Ped;Prob", "Probability of Fit"); + MHCamera dispped9 (geomcam, "Ped;DeltaPedestalMean", "Rel. Diff. Mean per Slice (Calc.-Fitte)"); + MHCamera dispped11 (geomcam, "Ped;DeltaPedestalMeanError", "Rel. Diff. Mean Error per Slice (Calc.-Fitted)"); + MHCamera dispped12 (geomcam, "Ped;DeltaRmsSigma", "Rel. Diff. RMS per Slice (Calc.-Fitted)"); + MHCamera dispped14 (geomcam, "Ped;DeltaRmsSigmaError", "Rel. Diff. RMS Error per Slice (Calc.-Fitted)"); + dispped0.SetCamContent(pedcam, 0); dispped0.SetCamError(pedcam, 1); @@ -147,5 +149,5 @@ // Histogram values - TCanvas &b1 = d1->AddTab("Direct Calculation"); + TCanvas &b1 = display->AddTab("Ped.Calc."); b1.Divide(4,3); @@ -156,5 +158,5 @@ // Fitted values - TCanvas &b2 = d1->AddTab("Fits"); + TCanvas &b2 = display->AddTab("Ped.Fit"); b2.Divide(4,3); @@ -166,5 +168,5 @@ // Fits Probability - TCanvas &b3 = d1->AddTab("Fit Probabilities"); + TCanvas &b3 = display->AddTab("Ped.Fit Prob."); b3.Divide(1,3); @@ -172,5 +174,5 @@ // Differences - TCanvas &c4 = d1->AddTab("Relative Difference Calculation-Fits"); + TCanvas &c4 = display->AddTab("Rel.Diff.Calc.-Fit"); c4.Divide(4,3); @@ -180,4 +182,5 @@ CamDraw(c4,dispped14,pedcam,4,4,1); + // // Create a empty Parameter List and an empty Task List @@ -201,5 +204,4 @@ // tlist2.AddToList(&geomapl); - // // Now setup the new tasks and tasklist for the calibration @@ -273,5 +275,5 @@ // of the arrival times using a spline, uncomment the next line // - tlist2.AddToList(&timecalc); + // tlist2.AddToList(&timecalc); tlist2.AddToList(&calcalc); @@ -281,7 +283,5 @@ MEvtLoop evtloop2; evtloop2.SetParList(&plist2); - - MProgressBar bar; - evtloop2.SetProgressBar(&bar); + evtloop2.SetDisplay(display); // @@ -301,80 +301,142 @@ // just one example how to get the plots of individual pixels // - // calcam[18].DrawClone(); - - MHCamera disp1 (geomcam, "MCalibrationPix;Charge", "Fitted Mean Charges"); - MHCamera disp3 (geomcam, "MCalibrationPix;SigmaCharge", "Sigma of Fitted Charges"); - MHCamera disp5 (geomcam, "MCalibrationPix;ChargeProb", "Probability of Fit"); - MHCamera disp6 (geomcam, "MCalibrationPix;Time", "Rel. Arrival Times"); - MHCamera disp7 (geomcam, "MCalibrationPix;SigmaTime", "Sigma of Rel. Arrival Times"); - MHCamera disp8 (geomcam, "MCalibrationPix;TimeProb", "Probability of Time Fit"); - MHCamera disp9 (geomcam, "MCalibrationPix;Ped", "Pedestals"); - MHCamera disp10 (geomcam, "MCalibrationPix;PedRms", "Pedestal RMS"); - MHCamera disp11 (geomcam, "MCalibrationPix;RSigma", "Reduced Sigmas"); - MHCamera disp12 (geomcam, "MCalibrationPix;PheFFactorMethod", "Nr. of Phe's (F-Factor Method)"); - MHCamera disp13 (geomcam, "MCalibrationPix;MeanConversionFFactorMethod", - "Conversion Factor (F-Factor Method)"); - MHCamera disp14 (geomcam, "MCalibrationPix;MeanPhotInsidePlexiglass", - "Nr. of Photons (Blind Pixel Method)"); - MHCamera disp15 (geomcam, "MCalibrationPix;MeanConversionBlindPixelMethod", - "Conversion Factor (Blind Pixel Method)"); - MHCamera disp16 (geomcam, "MCalibrationPix;RSigma/Charge", "Reduced Sigma per Charge"); - MHCamera disp17 (geomcam, "MCalibrationPix;AbsTimeMean", "Abs. Arrival Times"); - MHCamera disp18 (geomcam, "MCalibrationPix;AbsTimeRms", "RMS of Arrival Times"); - - disp1.SetCamContent(calcam, 0); - disp1.SetCamError(calcam,1); - - disp3.SetCamContent(calcam, 2); - disp3.SetCamError(calcam,3); - - disp5.SetCamContent(calcam, 4); - - disp6.SetCamContent(calcam, 5); - disp6.SetCamError(calcam, 6); - disp7.SetCamContent(calcam, 6); - disp8.SetCamContent(calcam, 7); - - disp9.SetCamContent(calcam, 8); - disp9.SetCamError(calcam, 9); - - disp10.SetCamContent(calcam, 9); - disp11.SetCamContent(calcam, 10); - - disp12.SetCamContent(calcam, 11); - disp12.SetCamError(calcam, 12); - - disp13.SetCamContent(calcam, 13); - disp13.SetCamError(calcam, 14); - - disp14.SetCamContent(calcam, 15); - disp15.SetCamContent(calcam, 16); - disp16.SetCamContent(calcam, 17); - - disp17.SetCamContent(calcam, 18); - disp17.SetCamError(calcam, 19); - disp18.SetCamContent(calcam, 20); - disp18.SetCamError(calcam, 21); - - disp1.SetYTitle("Charge [FADC counts]"); - disp3.SetYTitle("\\sigma_{Charge} [FADC counts]"); - disp5.SetYTitle("P_{Charge} [1]"); - disp6.SetYTitle("Rel. Arr. Time [ns]"); - disp7.SetYTitle("\\sigma_{Time} [ns]"); - disp8.SetYTitle("P_{Time} [1]"); - disp9.SetYTitle("Ped [FADC Counts ]"); - disp10.SetYTitle("RMS_{Ped} [FADC Counts ]"); - disp11.SetYTitle("\\sqrt{\\sigma^{2}_{Charge} - RMS^{2}_{Ped}} [FADC Counts]"); - 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]"); - disp17.SetYTitle("Mean Abs. Time [FADC slice]"); - disp18.SetYTitle("RMS Abs. Time [FADC slices]"); - - MStatusDisplay *d3 = new MStatusDisplay; - d3->SetUpdateTime(3000); - d3->Resize(850,700); + calcam[356].DrawClone(); + + // Create histograms to display + 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"); + + + MHCamera disp20 (geomcam, "cal;Ped", "Pedestals"); + MHCamera disp21 (geomcam, "cal;PedRms", "Pedestal RMS"); + + MHCamera disp22 (geomcam, "cal;Time", "Rel. Arrival Times"); + MHCamera disp23 (geomcam, "cal;SigmaTime", "Sigma of Rel. Arrival Times"); + MHCamera disp24 (geomcam, "cal;TimeProb", "Probability of Time Fit"); + + MHCamera disp25 (geomcam, "cal;AbsTimeMean", "Abs. Arrival Times"); + MHCamera disp26 (geomcam, "cal;AbsTimeRms", "RMS of Arrival Times"); + + + // Fitted charge means and sigmas + disp1.SetCamContent(calcam, 0); + disp1.SetCamError( calcam, 1); + disp2.SetCamContent(calcam, 2); + disp2.SetCamError( calcam, 3); + // Fit probabilities + disp3.SetCamContent(calcam, 4); + + // Reduced Sigmas and reduced sigmas per charge + disp4.SetCamContent(calcam, 5); + disp4.SetCamError( calcam, 6); + disp5.SetCamContent(calcam, 7); + disp5.SetCamError( calcam, 8); + + // F-Factor Method + disp6.SetCamContent(calcam, 9); + disp6.SetCamError( calcam, 10); + disp7.SetCamContent(calcam, 11); + disp7.SetCamError( calcam, 12); + disp8.SetCamContent(calcam, 13); + disp8.SetCamError( calcam, 14); + + /// Blind Pixel Method + disp9.SetCamContent(calcam, 15); + disp9.SetCamError( calcam, 16); + disp10.SetCamContent(calcam,17); + disp10.SetCamError( calcam,18); + disp11.SetCamContent(calcam,19); + disp11.SetCamError( calcam,20); + + // PIN Diode Method + disp12.SetCamContent(calcam,21); + disp12.SetCamError( calcam,22); + disp13.SetCamContent(calcam,23); + disp13.SetCamError( calcam,24); + disp14.SetCamContent(calcam,25); + disp14.SetCamError( calcam,26); + + // Pixels with defects + disp15.SetCamContent(calcam,27); + disp16.SetCamContent(calcam,28); + disp17.SetCamContent(calcam,29); + disp18.SetCamContent(calcam,30); + + // Lo Gain calibration + disp19.SetCamContent(calcam,31); + + + // Pedestals + disp20.SetCamContent(calcam,35); + disp20.SetCamError( calcam,36); + disp21.SetCamContent(calcam,37); + disp21.SetCamError( calcam,38); + + + // Relative Times + disp20.SetCamContent(calcam,39); + disp20.SetCamError( calcam,40); + disp21.SetCamContent(calcam,41); + disp21.SetCamError( calcam,42); + disp22.SetCamContent(calcam,43); + + // Absolute Times + disp23.SetCamContent(calcam,44); + disp23.SetCamError( calcam,45); + disp24.SetCamContent(calcam,46); + disp24.SetCamError( calcam,47); + + + disp1.SetYTitle("Charge [FADC units]"); + 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]"); + + disp19.SetYTitle("Ped [FADC Counts ]"); + disp20.SetYTitle("RMS_{Ped} [FADC Counts ]"); + + disp21.SetYTitle("Rel. Arr. Time [ns]"); + disp22.SetYTitle("\\sigma_{Time} [ns]"); + disp23.SetYTitle("P_{Time} [1]"); + + disp24.SetYTitle("Mean Abs. Time [FADC slice]"); + disp25.SetYTitle("RMS Abs. Time [FADC slices]"); + gStyle->SetOptStat(1111); @@ -382,59 +444,86 @@ // Charges - TCanvas &c1 = d3->AddTab("Fitted Charges"); - c1.Divide(2,3); - - CamDraw(c1,disp1,calcam,1,2,1); - CamDraw(c1,disp3,calcam,2,2,1); + TCanvas &c1 = display->AddTab("Fit.Charge"); + c1.Divide(2, 3); + + CamDraw(c1, disp1,calcam,1, 2 , 2); + CamDraw(c1, disp2,calcam,2, 2 , 2); // Fit Probability - TCanvas &c2 = d3->AddTab("Fit Prob."); + TCanvas &c2 = display->AddTab("Fit.Prob"); c2.Divide(1,3); - CamDraw(c2,disp5,calcam,1,1,3); + CamDraw(c2, disp3,calcam,1, 1 , 4); + + // Reduced Sigmas + TCanvas &c3 = display->AddTab("Red.Sigma"); + c3.Divide(2,3); + + CamDraw(c3, disp4,calcam,1, 2 , 2); + CamDraw(c3, disp5,calcam,2, 2 , 2); + + // F-Factor Method + TCanvas &c4 = display->AddTab("F-Factor"); + c4.Divide(3,3); + + CamDraw(c4, disp6,calcam,1, 3 , 2); + CamDraw(c4, disp7,calcam,2, 3 , 2); + CamDraw(c4, disp8,calcam,3, 3 , 2); + + // Blind Pixel Method + TCanvas &c5 = display->AddTab("BlindPix"); + c5.Divide(3, 3); + + CamDraw(c5, disp9,calcam,1, 3 , 9); + CamDraw(c5, disp10,calcam,2, 3 , 2); + CamDraw(c5, disp11,calcam,3, 3 , 2); + + // PIN Diode Method + TCanvas &c6 = display->AddTab("PINDiode"); + c6.Divide(3,3); + + CamDraw(c6, disp12,calcam,1, 3 , 9); + CamDraw(c6, disp13,calcam,2, 3 , 2); + CamDraw(c6, disp14,calcam,3, 3 , 2); + + // Defects + TCanvas &c7 = display->AddTab("Defects"); + c7.Divide(4,2); + + CamDraw(c7, disp15,calcam,1,4, 0); + CamDraw(c7, disp16,calcam,2,4, 0); + CamDraw(c7, disp17,calcam,3,4, 0); + CamDraw(c7, disp18,calcam,4,4,0); + + // Lo Gain Calibration + TCanvas &c8 = display->AddTab("LowGain"); + c8.Divide(1,3); + + CamDraw(c8, disp19,calcam,1,4,0); + + + // Pedestals + TCanvas &c9 = display->AddTab("Pedestals"); + c9.Divide(2,3); + + CamDraw(c9,disp20,calcam,1,3,1); + CamDraw(c9,disp21,calcam,2,3,2); + // Rel. Times - TCanvas &c3 = d3->AddTab("Fitted Rel. Times"); - c3.Divide(3,3); - - CamDraw(c3,disp6,calcam,1,3,1); - CamDraw(c3,disp7,calcam,2,3,1); - CamDraw(c3,disp8,calcam,3,3,0); + TCanvas &c10 = display->AddTab("Fitted Rel. Times"); + c10.Divide(3,3); + + CamDraw(c10,disp22,calcam,1,3,2); + CamDraw(c10,disp23,calcam,2,3,2); + CamDraw(c10,disp24,calcam,3,3,4); + // Abs. Times - TCanvas &c35 = d3->AddTab("Abs. Times"); - c35.Divide(2,3); - - CamDraw(c35,disp17,calcam,1,2,1); - CamDraw(c35,disp18,calcam,2,2,1); - - // Pedestals - TCanvas &c4 = d3->AddTab("Pedestals"); - c4.Divide(2,3); - - CamDraw(c4,disp9,calcam,1,2,0); - CamDraw(c4,disp10,calcam,2,2,1); - - // Reduced Sigmas - TCanvas &c5 = d3->AddTab("Reduced Sigmas"); - c5.Divide(2,3); - - // CamDraw(c5,disp11,calcam,1,2,1); - CamDraw(c5,disp11,calcam,1,2,2); - CamDraw(c5,disp16,calcam,2,2,1); - - // F-Factor Method - TCanvas &c6 = d3->AddTab("F-Factor Method"); - c6.Divide(2,3); - - CamDraw(c6,disp12,calcam,1,2,1); - CamDraw(c6,disp13,calcam,2,2,1); - - // Blind Pixel Method - TCanvas &c7 = d3->AddTab("Blind Pixel Method"); - c7.Divide(2, 3); - - CamDraw(c7,disp14,calcam,1,2,9); - CamDraw(c7,disp15,calcam,2,2,1); + TCanvas &c11 = display->AddTab("Abs. Times"); + c11.Divide(2,3); + + CamDraw(c11,disp25,calcam,1,2,2); + CamDraw(c11,disp26,calcam,2,2,2); } @@ -452,110 +541,101 @@ obj1->Draw(); ((MHCamera*)obj1)->SetPrettyPalette(); - - c.cd(i+2*j); - gPad->SetBorderMode(0); - 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.5066283; - const Double_t mean = obj2->GetMean(); - const Double_t rms = obj2->GetRMS(); - const Double_t width = max-min; - - if (rms == 0. || width == 0. ) - return; - - switch (fit) + + if (fit != 0) { - case 0: - TF1 *sgaus = new TF1("sgaus","gaus(0)",min,max); - sgaus->SetBit(kCanDelete); - sgaus->SetParNames("Area","#mu","#sigma"); - sgaus->SetParameters(integ/rms,mean,rms); - sgaus->SetParLimits(0,0.,integ); - sgaus->SetParLimits(1,min,max); - sgaus->SetParLimits(2,0,width/1.5); - obj2->Fit("sgaus","QLR"); - obj2->GetFunction("sgaus")->SetLineColor(kYellow); - break; - - case 1: - TString dgausform = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"; - dgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])"; - TF1 *dgaus = new TF1("dgaus",dgausform.Data(),min,max); - dgaus->SetBit(kCanDelete); - dgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}","A_{2}","#mu_{2}","#sigma_{2}"); - dgaus->SetParameters(integ,(min+mean)/2.,width/4., - integ/width/2.,(max+mean)/2.,width/4.); - // The left-sided Gauss - dgaus->SetParLimits(0,integ-1.5,integ+1.5); - dgaus->SetParLimits(1,min+(width/10.),mean); - dgaus->SetParLimits(2,0,width/2.); - // The right-sided Gauss - dgaus->SetParLimits(3,0,integ); - dgaus->SetParLimits(4,mean,max-(width/10.)); - dgaus->SetParLimits(5,0,width/2.); - obj2->Fit("dgaus","QLRM"); - obj2->GetFunction("dgaus")->SetLineColor(kYellow); - break; - - case 2: - TString tgausform = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"; - tgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])"; - tgausform += "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])"; - TF1 *tgaus = new TF1("tgaus",tgausform.Data(),min,max); - tgaus->SetBit(kCanDelete); - tgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}", - "A_{2}","#mu_{2}","#sigma_{2}", - "A_{3}","#mu_{3}","#sigma_{3}"); - tgaus->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 - tgaus->SetParLimits(0,integ-1.5,integ+1.5); - tgaus->SetParLimits(1,min+(width/10.),mean); - tgaus->SetParLimits(2,width/15.,width/2.); - // The right-sided Gauss - tgaus->SetParLimits(3,0.,integ); - tgaus->SetParLimits(4,mean,max-(width/10.)); - tgaus->SetParLimits(5,width/15.,width/2.); - // The Gauss describing the outliers - tgaus->SetParLimits(6,0.,integ); - tgaus->SetParLimits(7,min,max); - tgaus->SetParLimits(8,width/4.,width/1.5); - obj2->Fit("tgaus","QLRM"); - obj2->GetFunction("tgaus")->SetLineColor(kYellow); - break; - case 3: - obj2->Fit("pol0","Q"); - obj2->GetFunction("pol0")->SetLineColor(kYellow); - break; - case 9: - break; - default: - obj2->Fit("gaus","Q"); - obj2->GetFunction("gaus")->SetLineColor(kYellow); - break; + c.cd(i+2*j); + gPad->SetBorderMode(0); + 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.5066283; + const Double_t mean = obj2->GetMean(); + const Double_t rms = obj2->GetRMS(); + const Double_t width = max-min; + + if (rms == 0. || width == 0. ) + return; + + switch (fit) + { + case 1: + TF1 *sgaus = new TF1("sgaus","gaus(0)",min,max); + sgaus->SetBit(kCanDelete); + sgaus->SetParNames("Area","#mu","#sigma"); + sgaus->SetParameters(integ/rms,mean,rms); + sgaus->SetParLimits(0,0.,integ); + sgaus->SetParLimits(1,min,max); + sgaus->SetParLimits(2,0,width/1.5); + obj2->Fit("sgaus","QLR"); + obj2->GetFunction("sgaus")->SetLineColor(kYellow); + break; + + case 2: + TString dgausform = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"; + dgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])"; + TF1 *dgaus = new TF1("dgaus",dgausform.Data(),min,max); + dgaus->SetBit(kCanDelete); + dgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}","A_{2}","#mu_{2}","#sigma_{2}"); + dgaus->SetParameters(integ,(min+mean)/2.,width/4., + integ/width/2.,(max+mean)/2.,width/4.); + // The left-sided Gauss + dgaus->SetParLimits(0,integ-1.5,integ+1.5); + dgaus->SetParLimits(1,min+(width/10.),mean); + dgaus->SetParLimits(2,0,width/2.); + // The right-sided Gauss + dgaus->SetParLimits(3,0,integ); + dgaus->SetParLimits(4,mean,max-(width/10.)); + dgaus->SetParLimits(5,0,width/2.); + obj2->Fit("dgaus","QLRM"); + obj2->GetFunction("dgaus")->SetLineColor(kYellow); + break; + + case 3: + TString tgausform = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"; + tgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])"; + tgausform += "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])"; + TF1 *tgaus = new TF1("tgaus",tgausform.Data(),min,max); + tgaus->SetBit(kCanDelete); + tgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}", + "A_{2}","#mu_{2}","#sigma_{2}", + "A_{3}","#mu_{3}","#sigma_{3}"); + tgaus->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 + tgaus->SetParLimits(0,integ-1.5,integ+1.5); + tgaus->SetParLimits(1,min+(width/10.),mean); + tgaus->SetParLimits(2,width/15.,width/2.); + // The right-sided Gauss + tgaus->SetParLimits(3,0.,integ); + tgaus->SetParLimits(4,mean,max-(width/10.)); + tgaus->SetParLimits(5,width/15.,width/2.); + // The Gauss describing the outliers + tgaus->SetParLimits(6,0.,integ); + tgaus->SetParLimits(7,min,max); + tgaus->SetParLimits(8,width/4.,width/1.5); + obj2->Fit("tgaus","QLRM"); + obj2->GetFunction("tgaus")->SetLineColor(kYellow); + break; + case 4: + obj2->Fit("pol0","Q"); + obj2->GetFunction("pol0")->SetLineColor(kYellow); + break; + case 9: + break; + default: + obj2->Fit("gaus","Q"); + obj2->GetFunction("gaus")->SetLineColor(kYellow); + break; + } + + gPad->Modified(); + gPad->Update(); + } - - gPad->Modified(); - gPad->Update(); - - /* - TPaveStats *st = (TPaveStats*)gPad->GetPrimitive("stats"); - st->SetName("Projection"); - st->SetX1NDC(0.7); - st->SetX2NDC(1.2); - st->SetY1NDC(0.4); - st->SetY2NDC(1.2); - - gPad->Modified(); - gPad->Update(); - - */ } Index: /trunk/MagicSoft/Mars/mcalib/MCalibrationCam.cc =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MCalibrationCam.cc (revision 3058) +++ /trunk/MagicSoft/Mars/mcalib/MCalibrationCam.cc (revision 3059) @@ -35,6 +35,4 @@ // 3) It initializes a pointer to an MCalibrationPINDiode container // -// 4) -// ///////////////////////////////////////////////////////////////////////////// #include "MCalibrationCam.h" @@ -48,4 +46,5 @@ #include "MGeomCam.h" +#include "MGeomPix.h" #include "MCalibrationPix.h" @@ -62,5 +61,4 @@ const Int_t MCalibrationCam::gkBlindPixelId = 559; const Int_t MCalibrationCam::gkPINDiodeId = 9999; -const Float_t MCalibrationCam::gkTimeSliceWidth = 3.3; // -------------------------------------------------------------------------- @@ -383,98 +381,289 @@ return kFALSE; - if ( (!(*this)[idx].IsChargeFitValid()) || (*this)[idx].IsExcluded()) - return kFALSE; - - if (idx == gkBlindPixelId) - return kFALSE; - - if (idx == gkPINDiodeId) - return kFALSE; - switch (type) { case 0: + if ((*this)[idx].IsExcluded()) + return kFALSE; val = (*this)[idx].GetCharge(); break; case 1: + if ((*this)[idx].IsExcluded()) + return kFALSE; val = (*this)[idx].GetErrCharge(); break; case 2: + if ((*this)[idx].IsExcluded()) + return kFALSE; val = (*this)[idx].GetSigmaCharge(); break; case 3: + if ((*this)[idx].IsExcluded()) + return kFALSE; val = (*this)[idx].GetErrSigmaCharge(); break; case 4: + if ((*this)[idx].IsExcluded()) + return kFALSE; val = (*this)[idx].GetChargeProb(); break; case 5: - if (!(*this)[idx].IsTimeFitValid()) - return kFALSE; - val = (*this)[idx].GetMeanTimeOffset() * gkTimeSliceWidth; + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetRSigmaCharge(); break; case 6: - if (!(*this)[idx].IsTimeFitValid()) - return kFALSE; - val = (*this)[idx].GetTimingPrecision() * gkTimeSliceWidth; + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetErrRSigmaCharge(); break; case 7: - if (!(*this)[idx].IsTimeFitValid()) - return kFALSE; - val = (*this)[idx].GetTimeProb(); + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetRSigmaCharge() / (*this)[idx].GetCharge(); break; case 8: + if ((*this)[idx].IsExcluded()) + return kFALSE; + // relative error RsigmaCharge square + val = (*this)[idx].GetErrRSigmaCharge()* (*this)[idx].GetErrRSigmaCharge() + / ((*this)[idx].GetRSigmaCharge() * (*this)[idx].GetRSigmaCharge() ); + // relative error Charge square + val += (*this)[idx].GetErrCharge() * (*this)[idx].GetErrCharge() + / ((*this)[idx].GetCharge() * (*this)[idx].GetCharge() ); + // calculate relative error out of squares + val = TMath::Sqrt(val) ; + // multiply with value to get absolute error + val *= (*this)[idx].GetRSigmaCharge() / (*this)[idx].GetCharge(); + break; + case 9: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetPheFFactorMethod(); + break; + case 10: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetPheFFactorMethodError(); + break; + case 11: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetMeanConversionFFactorMethod(); + break; + case 12: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetErrorConversionFFactorMethod(); + break; + case 13: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetTotalFFactorFFactorMethod(); + break; + case 14: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetTotalFFactorErrorFFactorMethod(); + break; + case 15: + if ((*this)[idx].IsExcluded()) + return kFALSE; + if (cam[idx].IsInOuterRing()) + val = GetMeanPhotInsidePlexiglass()*gkCalibrationOutervsInnerPixelArea; + else + val = GetMeanPhotInsidePlexiglass(); + break; + case 16: + if ((*this)[idx].IsExcluded()) + return kFALSE; + if (cam[idx].IsInOuterRing()) + val = (double)fMeanPhotInsidePlexiglass*gkCalibrationOutervsInnerPixelArea; + else + val = (double)fMeanPhotInsidePlexiglass; + break; + case 17: + if ((*this)[idx].IsExcluded()) + return kFALSE; + if (cam[idx].IsInOuterRing()) + val = (*this)[idx].GetMeanConversionBlindPixelMethod()*gkCalibrationOutervsInnerPixelArea; + else + val = (*this)[idx].GetMeanConversionBlindPixelMethod(); + break; + case 18: + if ((*this)[idx].IsExcluded()) + return kFALSE; + if (cam[idx].IsInOuterRing()) + { + val = (*this)[idx].GetErrorConversionBlindPixelMethod()*(*this)[idx].GetErrorConversionBlindPixelMethod() + * gkCalibrationOutervsInnerPixelArea * gkCalibrationOutervsInnerPixelArea; + val += gkCalibrationOutervsInnerPixelAreaError * gkCalibrationOutervsInnerPixelAreaError + * (*this)[idx].GetMeanConversionBlindPixelMethod() *(*this)[idx].GetMeanConversionBlindPixelMethod(); + val = TMath::Sqrt(val); + } + else + val = (*this)[idx].GetErrorConversionBlindPixelMethod(); + break; + case 19: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetTotalFFactorBlindPixelMethod(); + break; + case 20: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetTotalFFactorErrorBlindPixelMethod(); + break; + case 21: + if ((*this)[idx].IsExcluded()) + return kFALSE; + if (cam[idx].IsInOuterRing()) + val = GetMeanPhotOutsidePlexiglass()*gkCalibrationOutervsInnerPixelArea; + else + val = GetMeanPhotOutsidePlexiglass(); + break; + case 22: + if ((*this)[idx].IsExcluded()) + return kFALSE; + if (cam[idx].IsInOuterRing()) + val = (double)fMeanPhotOutsidePlexiglass*gkCalibrationOutervsInnerPixelArea; + else + val = (double)fMeanPhotOutsidePlexiglass; + break; + case 23: + if ((*this)[idx].IsExcluded()) + return kFALSE; + if (cam[idx].IsInOuterRing()) + val = (*this)[idx].GetMeanConversionPINDiodeMethod()*gkCalibrationOutervsInnerPixelArea; + else + val = (*this)[idx].GetMeanConversionPINDiodeMethod(); + break; + case 24: + if ((*this)[idx].IsExcluded()) + return kFALSE; + if (cam[idx].IsInOuterRing()) + { + val = (*this)[idx].GetErrorConversionPINDiodeMethod()*(*this)[idx].GetErrorConversionPINDiodeMethod() + * gkCalibrationOutervsInnerPixelArea * gkCalibrationOutervsInnerPixelArea; + val += gkCalibrationOutervsInnerPixelAreaError * gkCalibrationOutervsInnerPixelAreaError + * (*this)[idx].GetMeanConversionPINDiodeMethod() *(*this)[idx].GetMeanConversionPINDiodeMethod(); + val = TMath::Sqrt(val); + } + else + val = (*this)[idx].GetErrorConversionPINDiodeMethod(); + break; + case 25: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetTotalFFactorBlindPixelMethod(); + break; + case 26: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = (*this)[idx].GetTotalFFactorErrorBlindPixelMethod(); + break; + case 27: + if ((*this)[idx].IsExcluded()) + val = 1.; + else + return kFALSE; + break; + case 28: + if ((*this)[idx].IsExcluded()) + return kFALSE; + if (!(*this)[idx].IsFitted()) + val = 1; + else + return kFALSE; + break; + case 29: + if ((*this)[idx].IsExcluded()) + return kFALSE; + if (!(*this)[idx].IsFitted()) + return kFALSE; + if (!(*this)[idx].IsChargeFitValid()) + val = 1; + else + return kFALSE; + break; + case 30: + if ((*this)[idx].IsExcluded()) + return kFALSE; + if ((*this)[idx].IsOscillating()) + val = 1; + else + return kFALSE; + break; + case 31: + if ((*this)[idx].IsExcluded()) + return kFALSE; + if ((*this)[idx].IsHiGainSaturation()) + val = 1; + else + return kFALSE; + break; + case 32: + if ((*this)[idx].IsExcluded()) + return kFALSE; + if ((*this)[idx].IsFFactorMethodValid()) + val = 1; + else + return kFALSE; + break; + case 33: + if ((*this)[idx].IsExcluded()) + return kFALSE; + if ((*this)[idx].IsBlindPixelMethodValid()) + val = 1; + else + return kFALSE; + break; + case 34: + if ((*this)[idx].IsExcluded()) + return kFALSE; + if ((*this)[idx].IsPINDiodeMethodValid()) + val = 1; + else + return kFALSE; + break; + case 35: + if ((*this)[idx].IsExcluded()) + return kFALSE; val = (*this)[idx].GetPed(); break; - case 9: + case 36: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = 1.; + // val = (*this)[idx].GetPedError(); + break; + case 37: + if ((*this)[idx].IsExcluded()) + return kFALSE; val = (*this)[idx].GetPedRms(); break; - case 10: - val = (*this)[idx].GetRSigmaCharge(); - break; - case 11: - val = (*this)[idx].GetPheFFactorMethod(); - break; - case 12: - val = (*this)[idx].GetPheFFactorMethodError(); - break; - case 13: - val = (*this)[idx].GetMeanConversionFFactorMethod(); - break; - case 14: - val = (*this)[idx].GetErrorConversionFFactorMethod(); - break; - case 15: - if (idx < 397) - val = (double)fMeanPhotInsidePlexiglass; - else - val = (double)fMeanPhotInsidePlexiglass*gkCalibrationOutervsInnerPixelArea; - break; - case 16: - if (idx < 397) - val = (*this)[idx].GetMeanConversionBlindPixelMethod(); - else - val = (*this)[idx].GetMeanConversionBlindPixelMethod()*gkCalibrationOutervsInnerPixelArea; - break; - case 17: - val = (*this)[idx].GetRSigmaCharge() / (*this)[idx].GetCharge(); - break; - case 18: - if (!(*this)[idx].IsTimeFitValid()) + case 38: + if ((*this)[idx].IsExcluded()) + return kFALSE; + val = 1.; + // val = (*this)[idx].GetPedRmsError(); + break; + case 39: + if ((*this)[idx].IsExcluded()) return kFALSE; val = (*this)[idx].GetAbsTimeMean(); break; - case 19: - if (!(*this)[idx].IsTimeFitValid()) + case 40: + if ((*this)[idx].IsExcluded()) return kFALSE; val = (*this)[idx].GetAbsTimeMeanErr(); break; - case 20: - if (!(*this)[idx].IsTimeFitValid()) + case 41: + if ((*this)[idx].IsExcluded()) return kFALSE; val = (*this)[idx].GetAbsTimeRms(); break; - case 21: - if (!(*this)[idx].IsTimeFitValid()) + case 42: + if ((*this)[idx].IsExcluded()) return kFALSE; val = (*this)[idx].GetAbsTimeMeanErr()/TMath::Sqrt(2.);