Index: /trunk/Mars/fact/analysis/gain/fit_spectra.C =================================================================== --- /trunk/Mars/fact/analysis/gain/fit_spectra.C (revision 17198) +++ /trunk/Mars/fact/analysis/gain/fit_spectra.C (revision 17199) @@ -18,4 +18,5 @@ #include "MArrayI.h" #include "MRawRunHeader.h" +#include "PixelMap.h" using namespace std; @@ -31,5 +32,5 @@ const Double_t cross = par[3]; const Double_t shift = par[4]; - const Double_t noise = par[5]; + const Double_t noise = par[5]<0 ? sigma : par[5]; const Double_t expo = par[6]; @@ -67,5 +68,5 @@ const Double_t sigma = func.GetParameter(2)*func.GetParameter(1); const Double_t cross = func.GetParameter(3); - const Double_t noise = func.GetParameter(5); + const Double_t noise = func.GetParameter(5)<0 ? sigma : func.GetParameter(5); const Double_t expo = func.GetParameter(6); @@ -94,5 +95,6 @@ cout << "Crosstalk: " << setw(8) << f.GetParameter(3) << " +/- " << f.GetParError(3) << endl; cout << "Pcrosstalk: " << setw(8) << xtalk(f) << endl; - cout << "Noise: " << setw(8) << f.GetParameter(5)/sqrt(integration_window) << " +/- " << f.GetParError(5)/sqrt(integration_window) << endl; + if (f.GetParameter(5)>=0) + cout << "Noise: " << setw(8) << f.GetParameter(5)/sqrt(integration_window) << " +/- " << f.GetParError(5)/sqrt(integration_window) << endl; cout << "Expo: " << setw(8) << f.GetParameter(6) << " +/- " << f.GetParError(6) << endl; //cout << "--------------------" << endl; @@ -104,7 +106,16 @@ // The parameters for the function are the filenam from the gain extraction // and the output filename -int fit_spectra(const char *filename = "20130222_018_018.root", - const char *outfile = "20120802_247_247-fit.root" ) +int fit_spectra(const char *filename = "20130222_018_018.root", + const char *outfile = "20130222_018_018-fit.root", + bool fixednoise=true) { + // Read the mapping between bias channels and hardware channels + PixelMap pmap; + if (!pmap.Read("FACTmap111030.txt")) + { + cout << "FACTmap111030.txt not found." << endl; + return 1; + } + // It is more correct to fit the integral, but this is super // slow, especially for 1440 pixel. To get that, one would have @@ -223,28 +234,42 @@ // Set name and title for the histograms - cRate.SetNameTitle ("Rate", "Dark count rate"); - cGain.SetNameTitle ("Gain", "Gain distribution"); - cRelSigma.SetNameTitle ("RelSigma", "Rel. Sigma"); - cCrosstalk.SetNameTitle("Crosstalk", "Crosstalk probability"); - cBaseline.SetNameTitle ("Baseline", "Baseline per sample"); - cNoise.SetNameTitle ("Noise", "Noise per sample"); - cChi2.SetNameTitle ("Chi2", "\\Chi^2"); - cNormGain.SetNameTitle ("NormGain", "Normalized gain"); - cFitProb.SetNameTitle ("FitProb", "Root's fit probability"); - cCrosstalkP.SetNameTitle("Pxtalk", "Crosstalk coeff. P"); - cCoeffR.SetNameTitle ("CoeffR", "Coefficient R"); + cRate.SetNameTitle ("Rate", "Dark count rate"); + cGain.SetNameTitle ("Gain", "Gain distribution"); + cRelSigma.SetNameTitle ("RelSigma", "Rel. Sigma"); + cCrosstalk.SetNameTitle ("Crosstalk", "Crosstalk probability"); + cBaseline.SetNameTitle ("Baseline", "Baseline per sample"); + cNoise.SetNameTitle ("Noise", "Noise per sample"); + cChi2.SetNameTitle ("Chi2", "\\Chi^2"); + cNormGain.SetNameTitle ("NormGain", "Normalized gain"); + cFitProb.SetNameTitle ("FitProb", "Root's fit probability"); + cCrosstalkP.SetNameTitle("Pxtalk", "Crosstalk coeff. P"); + cCoeffR.SetNameTitle ("CoeffR", "Coefficient R"); // Instantiate 1D histograms for the distributions - TH1F hRate ("Rate", "Dark count rate", 150, 0, 15); - TH1F hGain ("Gain", "Gain distribution", 100, 0, 400); - TH1F hRelSigma ("RelSigma", "Rel. Sigma", 160, 0, 0.40); - TH1F hCrosstalk("Crosstalk", "Crosstalk probability", 90, 0, 0.30); - TH1F hBaseline ("Baseline", "Baseline per sample", 75, -7.5, 7.5); - TH1F hNoise ("Noise", "Noise per sample", 60, 0, 30); - TH1F hChi2 ("Chi2", "\\Chi^2", 200, 0, 4); - TH1F hNormGain ("NormGain", "Normalized gain", 51, 0.5, 1.5); - TH1F hFitProb ("FitProb", "FitProb distribution", 100, 0, 1); - TH1F hCrosstalkP("Pxtalk", "Crosstalk coeff.", 90, 0, 0.3); - TH1F hCoeffR ("CoeffR", "Coefficient R", 90, -1, 2); + // including TM channels + TH1F hRate1 ("Rate1", "Dark count rate", 150, 0, 15); + TH1F hGain1 ("Gain1", "Gain distribution", 100, 0, 400); + TH1F hRelSigma1 ("RelSigma1", "Rel. Sigma", 160, 0, 0.40); + TH1F hCrosstalk1 ("Crosstalk1", "Crosstalk probability", 90, 0, 0.30); + TH1F hBaseline1 ("Baseline1", "Baseline per sample", 75, -7.5, 7.5); + TH1F hNoise1 ("Noise1", "Noise per sample", 60, 0, 30); + TH1F hChiSq1 ("ChiSq1", "\\Chi^2", 200, 0, 4); + TH1F hNormGain1 ("NormGain1", "Normalized gain", 51, 0.5, 1.5); + TH1F hFitProb1 ("FitProb1", "FitProb distribution", 100, 0, 1); + TH1F hCrosstalkP1("Pxtalk1", "Crosstalk coeff.", 90, 0, 0.3); + TH1F hCoeffR1 ("CoeffR1", "Coefficient R", 90, -1, 2); + + // excluding TM channels + TH1F hRate2 ("Rate2", "Dark count rate", 150, 0, 15); + TH1F hGain2 ("Gain2", "Gain distribution", 100, 0, 400); + TH1F hRelSigma2 ("RelSigma2", "Rel. Sigma", 160, 0, 0.40); + TH1F hCrosstalk2 ("Crosstalk2", "Crosstalk probability", 90, 0, 0.30); + TH1F hBaseline2 ("Baseline2", "Baseline per sample", 75, -7.5, 7.5); + TH1F hNoise2 ("Noise2", "Noise per sample", 60, 0, 30); + TH1F hChiSq2 ("ChiSq2", "\\Chi^2", 200, 0, 4); + TH1F hNormGain2 ("NormGain2", "Normalized gain", 51, 0.5, 1.5); + TH1F hFitProb2 ("FitProb2", "FitProb distribution", 100, 0, 1); + TH1F hCrosstalkP2("Pxtalk2", "Crosstalk coeff.", 90, 0, 0.3); + TH1F hCoeffR2 ("CoeffR2", "Coefficient R", 90, -1, 2); // Histigram for the sum of all spectrums @@ -257,11 +282,19 @@ // Histogram for the sum of all pixels excluding the ones with faulty fits - TH1F hSumClear("SumC", "Signal spectrum of all pixels", - hsignal->GetNbinsY(), hsignal->GetYaxis()->GetXmin(), hsignal->GetYaxis()->GetXmax()); - hSumClear.SetXTitle("pulse integral [mV sample]"); - hSumClear.SetYTitle("Counts"); - hSumClear.SetStats(false); - hSumClear.SetLineColor(kBlue); - hSumClear.Sumw2(); + TH1F hSumClear1("SumC1", "Signal spectrum of all pixels (incl TM)", + hsignal->GetNbinsY(), hsignal->GetYaxis()->GetXmin(), hsignal->GetYaxis()->GetXmax()); + hSumClear1.SetXTitle("pulse integral [mV sample]"); + hSumClear1.SetYTitle("Counts"); + hSumClear1.SetStats(false); + hSumClear1.SetLineColor(kBlue); + hSumClear1.Sumw2(); + + TH1F hSumClear2("SumC2", "Signal spectrum of all pixels (excp TM)", + hsignal->GetNbinsY(), hsignal->GetYaxis()->GetXmin(), hsignal->GetYaxis()->GetXmax()); + hSumClear2.SetXTitle("pulse integral [mV sample]"); + hSumClear2.SetYTitle("Counts"); + hSumClear2.SetStats(false); + hSumClear2.SetLineColor(kBlue); + hSumClear2.Sumw2(); // Arrival time spectrum of the extracted pulses @@ -278,9 +311,15 @@ // Spektrum for the normalized individual spectra - TH1F hSumScale("Sum2", "Signal spectrum of all pixels", 120, 0.05, 12.05); - hSumScale.SetXTitle("pulse integral [pe]"); - hSumScale.SetYTitle("Rate"); - hSumScale.SetStats(false); - hSumScale.Sumw2(); + TH1F hSumScale1("SumScale1", "Signal spectrum of all pixels (incl TM)", 120, 0.05, 12.05); + hSumScale1.SetXTitle("pulse integral [pe]"); + hSumScale1.SetYTitle("Rate"); + hSumScale1.SetStats(false); + hSumScale1.Sumw2(); + + TH1F hSumScale2("SumScale2", "Signal spectrum of all pixels (excl TM)", 120, 0.05, 12.05); + hSumScale2.SetXTitle("pulse integral [pe]"); + hSumScale2.SetYTitle("Rate"); + hSumScale2.SetStats(false); + hSumScale2.Sumw2(); // define fit function for Amplitudespectrum @@ -342,12 +381,17 @@ func.SetParLimits(2, 0, 1); // Sigma func.SetParLimits(3, 0, 1); // Crosstalk - func.SetParLimits(5, 0, 150); // Noise + if (!fixednoise) + func.SetParLimits(5, 0, 150); // Noise func.SetParameter(0, ampl); // Amplitude func.SetParameter(1, maxpos); // Gain func.SetParameter(2, 0.1); // Sigma - func.SetParameter(3, 0.15); // Crosstalk + func.SetParameter(3, 0.16); // Crosstalk func.SetParameter(4, 0*integration_window); // Baseline - func.SetParameter(5, 1.*sqrt(integration_window)); // Noise + if (fixednoise) + func.FixParameter(5, -1); // Noise + else + func.SetParameter(5, 0.1*maxpos); // Noise + func.SetParameter(6, 0.4); // Expo @@ -474,5 +518,5 @@ const float fCrosstlk = xtalk(func); const float fOffset = func.GetParameter(4); - const float fNoise = func.GetParameter(5)/sqrt(integration_window); + const float fNoise = func.GetParameter(5)<0 ? fGainRMS*fGain/sqrt(integration_window) : func.GetParameter(5)/sqrt(integration_window); const float fCoeffR = func.GetParameter(6); @@ -549,24 +593,46 @@ // Fill Parameters into histograms - hRate.Fill( fRate); - hGain.Fill( fGain); - hRelSigma.Fill( fGainRMS); - hCrosstalk.Fill( fCrosstlk); - hBaseline.Fill( fOffset/integration_window); - hNoise.Fill( fNoise); - hChi2.Fill( chi2); - hNormGain.Fill( fGain/gain); - hFitProb.Fill( fit_prob); - hCrosstalkP.Fill(fCrosstalkP); - hCoeffR.Fill( fCoeffR); + hRate1.Fill( fRate); + hGain1.Fill( fGain); + hRelSigma1.Fill( fGainRMS); + hCrosstalk1.Fill( fCrosstlk); + hBaseline1.Fill( fOffset/integration_window); + hNoise1.Fill( fNoise); + hChiSq1.Fill( chi2); + hNormGain1.Fill( fGain/gain); + hFitProb1.Fill( fit_prob); + hCrosstalkP1.Fill(fCrosstalkP); + hCoeffR1.Fill( fCoeffR); + + if (!pmap.index(pixel).isTM()) + { + hRate2.Fill( fRate); + hGain2.Fill( fGain); + hRelSigma2.Fill( fGainRMS); + hCrosstalk2.Fill( fCrosstlk); + hBaseline2.Fill( fOffset/integration_window); + hNoise2.Fill( fNoise); + hChiSq2.Fill( chi2); + hNormGain2.Fill( fGain/gain); + hFitProb2.Fill( fit_prob); + hCrosstalkP2.Fill(fCrosstalkP); + hCoeffR2.Fill( fCoeffR); + } // Fill sum spectrum for (int b=1; b<=hist->GetNbinsX(); b++) - hSumScale.Fill((hist->GetBinCenter(b)-fOffset)/fGain, hist->GetBinContent(b)); + hSumScale1.Fill((hist->GetBinCenter(b)-fOffset)/fGain, hist->GetBinContent(b)); + + if (!pmap.index(pixel).isTM()) + for (int b=1; b<=hist->GetNbinsX(); b++) + hSumScale2.Fill((hist->GetBinCenter(b)-fOffset)/fGain, hist->GetBinContent(b)); + delete hist; // Because of the rebinning... hist = hsignal->ProjectionY("proj", pixel+1, pixel+1); - hSumClear.Add(hist); + hSumClear1.Add(hist); + if (!pmap.index(pixel).isTM()) + hSumClear2.Add(hist); delete hist; @@ -647,24 +713,30 @@ canv->cd(1); - hh = hRate.DrawCopy(); + hh = hRate1.DrawCopy(); + hh = hRate2.DrawCopy("same"); canv->cd(2); - hh = hGain.DrawCopy(); + hh = hGain1.DrawCopy(); + hh = hGain2.DrawCopy("same"); hh->Fit("gaus"); canv->cd(3); - hh = hBaseline.DrawCopy(); + hh = hBaseline1.DrawCopy(); + hh = hBaseline2.DrawCopy("same"); hh->Fit("gaus"); canv->cd(4); - hh = hRelSigma.DrawCopy(); + hh = hRelSigma1.DrawCopy(); + hh = hRelSigma2.DrawCopy("same"); hh->Fit("gaus"); canv->cd(5); - hh = hCrosstalk.DrawCopy(); + hh = hCrosstalk1.DrawCopy(); + hh = hCrosstalk2.DrawCopy("same"); hh->Fit("gaus"); canv->cd(6); - hh = hNoise.DrawCopy(); + hh = hNoise1.DrawCopy(); + hh = hNoise2.DrawCopy("same"); hh->Fit("gaus"); @@ -676,16 +748,20 @@ canv->cd(1); gPad->SetLogy(); - hh = hFitProb.DrawCopy(); + hh = hFitProb1.DrawCopy(); + hh = hFitProb2.DrawCopy("same"); hh->Fit("gaus"); canv->cd(2); - hChi2.DrawCopy(); + hChiSq1.DrawCopy(); + hChiSq2.DrawCopy("same"); canv->cd(4); - hh = hCoeffR.DrawCopy(); + hh = hCoeffR1.DrawCopy(); + hh = hCoeffR2.DrawCopy("same"); hh->Fit("gaus"); canv->cd(5); - hh = hCrosstalkP.DrawCopy(); + hh = hCrosstalkP1.DrawCopy(); + hh = hCrosstalkP2.DrawCopy("same"); hh->Fit("gaus"); @@ -702,5 +778,6 @@ canv->cd(2); gPad->SetLogy(); - hh = hNormGain.DrawCopy(); + hh = hNormGain1.DrawCopy(); + hh = hNormGain2.DrawCopy("same"); hh->Fit("gaus"); @@ -708,10 +785,10 @@ gROOT->SetSelectedPad(0); c11.cd(); - hSumClear.DrawCopy("hist same"); + hSumClear1.DrawCopy("hist same"); //-------------------- fit gain corrected sum spectrum ------------------- gROOT->SetSelectedPad(0); - TCanvas &c11b = d->AddTab("CleanHist"); + TCanvas &c11b = d->AddTab("CleanHist1"); c11b.cd(); gPad->SetLogy(); @@ -719,6 +796,6 @@ gPad->SetGridy(); - const Int_t maxbin1 = hSumClear.GetMaximumBin(); - const Double_t ampl1 = hSumClear.GetBinContent(maxbin1); + const Int_t maxbin1 = hSumClear1.GetMaximumBin(); + const Double_t ampl1 = hSumClear1.GetBinContent(maxbin1); func.SetParameters(res_par); @@ -727,7 +804,7 @@ func.ReleaseParameter(6); - hSumClear.Fit(&func, fast?"LN0QSR":"LIN0QSR"); - - hSumClear.DrawCopy(); + hSumClear1.Fit(&func, fast?"LN0QSR":"LIN0QSR"); + + hSumClear1.DrawCopy(); func.DrawCopy("same"); @@ -739,5 +816,31 @@ gROOT->SetSelectedPad(0); - TCanvas &c12 = d->AddTab("GainHist"); + TCanvas &c11c = d->AddTab("CleanHist2"); + c11c.cd(); + gPad->SetLogy(); + gPad->SetGridx(); + gPad->SetGridy(); + + const Int_t maxbin1b = hSumClear2.GetMaximumBin(); + const Double_t ampl1b = hSumClear2.GetBinContent(maxbin1b); + + func.SetParameters(res_par); + func.SetParLimits(0, 0, 2*ampl1b); + func.SetParameter(0, ampl1b); + func.ReleaseParameter(6); + + hSumClear2.Fit(&func, fast?"LN0QSR":"LIN0QSR"); + + hSumClear2.DrawCopy(); + func.DrawCopy("same"); + + cout << "--------------------" << endl; + PrintFunc(func, integration_window); + cout << "--------------------" << endl; + + //-------------------- fit gain corrected sum spectrum ------------------- + + gROOT->SetSelectedPad(0); + TCanvas &c12 = d->AddTab("GainHist1"); c12.cd(); gPad->SetLogy(); @@ -745,11 +848,11 @@ gPad->SetGridy(); - const Int_t maxbin2 = hSumScale.GetMaximumBin(); - const Double_t ampl2 = hSumScale.GetBinContent(maxbin2); + const Int_t maxbin2 = hSumScale1.GetMaximumBin(); + const Double_t ampl2 = hSumScale1.GetBinContent(maxbin2); //Set fit parameters Double_t par2[7] = { - ampl2, 1, 0.1, res_par[3], 0, res_par[5]/res_par[1], res_par[6] + ampl2, 1, 0.1, res_par[3], 0, res_par[5]<0 ? -1 : res_par[5]/res_par[1], res_par[6] }; @@ -760,7 +863,40 @@ func.SetRange(0.62, 9); - hSumScale.Fit(&func, fast?"LN0QSR":"LIN0QSR"); - - hSumScale.DrawCopy(); + hSumScale1.Fit(&func, fast?"LN0QSR":"LIN0QSR"); + + hSumScale1.DrawCopy(); + func.DrawCopy("same"); + + cout << "--------------------" << endl; + PrintFunc(func, integration_window); + cout << "--------------------" << endl; + + //-------------------- fit gain corrected sum spectrum ------------------- + + gROOT->SetSelectedPad(0); + TCanvas &c12b = d->AddTab("GainHist2"); + c12b.cd(); + gPad->SetLogy(); + gPad->SetGridx(); + gPad->SetGridy(); + + const Int_t maxbin2b = hSumScale2.GetMaximumBin(); + const Double_t ampl2b = hSumScale2.GetBinContent(maxbin2b); + + //Set fit parameters + Double_t par2b[7] = + { + ampl2b, 1, 0.1, res_par[3], 0, res_par[5]<0 ? -1 : res_par[5]/res_par[1], res_par[6] + }; + + func.SetParameters(par2b); + func.SetParLimits(0, 0, 2*ampl2b); + func.FixParameter(1, 1); + func.FixParameter(4, 0); + + func.SetRange(0.62, 9); + hSumScale2.Fit(&func, fast?"LN0QSR":"LIN0QSR"); + + hSumScale2.DrawCopy(); func.DrawCopy("same"); @@ -794,3 +930,2 @@ return 0; } -