Changeset 17199 for trunk/Mars/fact


Ignore:
Timestamp:
10/04/13 18:55:41 (11 years ago)
Author:
tbretz
Message:
Implemented the option to fix the noise to sigma; fix noise to sigma as default; added new histograms which exclude the time marker channels; slightly modified start value for noise
File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/Mars/fact/analysis/gain/fit_spectra.C

    r17194 r17199  
    1818#include "MArrayI.h"
    1919#include "MRawRunHeader.h"
     20#include "PixelMap.h"
    2021
    2122using namespace std;
    … …  
    3132    const Double_t cross = par[3];
    3233    const Double_t shift = par[4];
    33     const Double_t noise = par[5];
     34    const Double_t noise = par[5]<0 ? sigma : par[5];
    3435    const Double_t expo  = par[6];
    3536
    … …  
    6768    const Double_t sigma = func.GetParameter(2)*func.GetParameter(1);
    6869    const Double_t cross = func.GetParameter(3);
    69     const Double_t noise = func.GetParameter(5);
     70    const Double_t noise = func.GetParameter(5)<0 ? sigma : func.GetParameter(5);
    7071    const Double_t expo  = func.GetParameter(6);
    7172
    … …  
    9495    cout << "Crosstalk:  " << setw(8) << f.GetParameter(3) << " +/- " << f.GetParError(3) << endl;
    9596    cout << "Pcrosstalk: " << setw(8) << xtalk(f) << endl;
    96     cout << "Noise:      " << setw(8) << f.GetParameter(5)/sqrt(integration_window) << " +/- " << f.GetParError(5)/sqrt(integration_window) << endl;
     97    if (f.GetParameter(5)>=0)
     98        cout << "Noise:      " << setw(8) << f.GetParameter(5)/sqrt(integration_window) << " +/- " << f.GetParError(5)/sqrt(integration_window) << endl;
    9799    cout << "Expo:       " << setw(8) << f.GetParameter(6) << " +/- " << f.GetParError(6) << endl;
    98100    //cout << "--------------------" << endl;
    … …  
    104106// The parameters for the function are the filenam from the gain extraction
    105107// and the output filename
    106 int fit_spectra(const char  *filename = "20130222_018_018.root",
    107                 const char  *outfile  = "20120802_247_247-fit.root" )
     108int fit_spectra(const char *filename = "20130222_018_018.root",
     109                const char *outfile  = "20130222_018_018-fit.root",
     110                bool fixednoise=true)
    108111{
     112    // Read the mapping between bias channels and hardware channels
     113    PixelMap pmap;
     114    if (!pmap.Read("FACTmap111030.txt"))
     115    {
     116        cout << "FACTmap111030.txt not found." << endl;
     117        return 1;
     118    }
     119
    109120    // It is more correct to fit the integral, but this is super
    110121    // slow, especially for 1440 pixel. To get that, one would have
    … …  
    223234
    224235    // Set name and title for the histograms
    225     cRate.SetNameTitle     ("Rate",      "Dark count rate");
    226     cGain.SetNameTitle     ("Gain",      "Gain distribution");
    227     cRelSigma.SetNameTitle ("RelSigma",  "Rel. Sigma");
    228     cCrosstalk.SetNameTitle("Crosstalk", "Crosstalk probability");
    229     cBaseline.SetNameTitle ("Baseline",  "Baseline per sample");
    230     cNoise.SetNameTitle    ("Noise",     "Noise per sample");
    231     cChi2.SetNameTitle     ("Chi2",      "\\Chi^2");
    232     cNormGain.SetNameTitle ("NormGain",  "Normalized gain");
    233     cFitProb.SetNameTitle  ("FitProb",   "Root's fit probability");
    234     cCrosstalkP.SetNameTitle("Pxtalk",   "Crosstalk coeff. P");
    235     cCoeffR.SetNameTitle   ("CoeffR",    "Coefficient R");
     236    cRate.SetNameTitle      ("Rate",      "Dark count rate");
     237    cGain.SetNameTitle      ("Gain",      "Gain distribution");
     238    cRelSigma.SetNameTitle  ("RelSigma",  "Rel. Sigma");
     239    cCrosstalk.SetNameTitle ("Crosstalk", "Crosstalk probability");
     240    cBaseline.SetNameTitle  ("Baseline",  "Baseline per sample");
     241    cNoise.SetNameTitle     ("Noise",     "Noise per sample");
     242    cChi2.SetNameTitle      ("Chi2",      "\\Chi^2");
     243    cNormGain.SetNameTitle  ("NormGain",  "Normalized gain");
     244    cFitProb.SetNameTitle   ("FitProb",   "Root's fit probability");
     245    cCrosstalkP.SetNameTitle("Pxtalk",    "Crosstalk coeff. P");
     246    cCoeffR.SetNameTitle    ("CoeffR",    "Coefficient R");
    236247
    237248    // Instantiate 1D histograms for the distributions
    238     TH1F hRate     ("Rate",      "Dark count rate",       150,  0,    15);
    239     TH1F hGain     ("Gain",      "Gain distribution",     100,  0,   400);
    240     TH1F hRelSigma ("RelSigma",  "Rel. Sigma",            160,  0,  0.40);
    241     TH1F hCrosstalk("Crosstalk", "Crosstalk probability",  90,  0,  0.30);
    242     TH1F hBaseline ("Baseline",  "Baseline per sample",    75, -7.5, 7.5);
    243     TH1F hNoise    ("Noise",     "Noise per sample",       60,  0,    30);
    244     TH1F hChi2     ("Chi2",      "\\Chi^2",               200,  0,     4);
    245     TH1F hNormGain ("NormGain",  "Normalized gain",        51,  0.5, 1.5);
    246     TH1F hFitProb  ("FitProb",   "FitProb distribution",  100,  0,     1);
    247     TH1F hCrosstalkP("Pxtalk",   "Crosstalk coeff.",       90,  0,   0.3);
    248     TH1F hCoeffR   ("CoeffR",    "Coefficient R",          90, -1,     2);
     249    // including TM channels
     250    TH1F hRate1      ("Rate1",      "Dark count rate",       150,  0,    15);
     251    TH1F hGain1      ("Gain1",      "Gain distribution",     100,  0,   400);
     252    TH1F hRelSigma1  ("RelSigma1",  "Rel. Sigma",            160,  0,  0.40);
     253    TH1F hCrosstalk1 ("Crosstalk1", "Crosstalk probability",  90,  0,  0.30);
     254    TH1F hBaseline1  ("Baseline1",  "Baseline per sample",    75, -7.5, 7.5);
     255    TH1F hNoise1     ("Noise1",     "Noise per sample",       60,  0,    30);
     256    TH1F hChiSq1     ("ChiSq1",     "\\Chi^2",               200,  0,     4);
     257    TH1F hNormGain1  ("NormGain1",  "Normalized gain",        51,  0.5, 1.5);
     258    TH1F hFitProb1   ("FitProb1",   "FitProb distribution",  100,  0,     1);
     259    TH1F hCrosstalkP1("Pxtalk1",    "Crosstalk coeff.",       90,  0,   0.3);
     260    TH1F hCoeffR1    ("CoeffR1",    "Coefficient R",          90, -1,     2);
     261
     262    // excluding TM channels
     263    TH1F hRate2      ("Rate2",      "Dark count rate",       150,  0,    15);
     264    TH1F hGain2      ("Gain2",      "Gain distribution",     100,  0,   400);
     265    TH1F hRelSigma2  ("RelSigma2",  "Rel. Sigma",            160,  0,  0.40);
     266    TH1F hCrosstalk2 ("Crosstalk2", "Crosstalk probability",  90,  0,  0.30);
     267    TH1F hBaseline2  ("Baseline2",  "Baseline per sample",    75, -7.5, 7.5);
     268    TH1F hNoise2     ("Noise2",     "Noise per sample",       60,  0,    30);
     269    TH1F hChiSq2     ("ChiSq2",     "\\Chi^2",               200,  0,     4);
     270    TH1F hNormGain2  ("NormGain2",  "Normalized gain",        51,  0.5, 1.5);
     271    TH1F hFitProb2   ("FitProb2",   "FitProb distribution",  100,  0,     1);
     272    TH1F hCrosstalkP2("Pxtalk2",    "Crosstalk coeff.",       90,  0,   0.3);
     273    TH1F hCoeffR2    ("CoeffR2",    "Coefficient R",          90, -1,     2);
    249274
    250275    // Histigram for the sum of all spectrums
    … …  
    257282
    258283    // Histogram for the sum of all pixels excluding the ones with faulty fits
    259     TH1F hSumClear("SumC", "Signal spectrum of all pixels",
    260                    hsignal->GetNbinsY(), hsignal->GetYaxis()->GetXmin(),  hsignal->GetYaxis()->GetXmax());
    261     hSumClear.SetXTitle("pulse integral [mV sample]");
    262     hSumClear.SetYTitle("Counts");
    263     hSumClear.SetStats(false);
    264     hSumClear.SetLineColor(kBlue);
    265     hSumClear.Sumw2();
     284    TH1F hSumClear1("SumC1", "Signal spectrum of all pixels (incl TM)",
     285                    hsignal->GetNbinsY(), hsignal->GetYaxis()->GetXmin(),  hsignal->GetYaxis()->GetXmax());
     286    hSumClear1.SetXTitle("pulse integral [mV sample]");
     287    hSumClear1.SetYTitle("Counts");
     288    hSumClear1.SetStats(false);
     289    hSumClear1.SetLineColor(kBlue);
     290    hSumClear1.Sumw2();
     291
     292    TH1F hSumClear2("SumC2", "Signal spectrum of all pixels (excp TM)",
     293                    hsignal->GetNbinsY(), hsignal->GetYaxis()->GetXmin(),  hsignal->GetYaxis()->GetXmax());
     294    hSumClear2.SetXTitle("pulse integral [mV sample]");
     295    hSumClear2.SetYTitle("Counts");
     296    hSumClear2.SetStats(false);
     297    hSumClear2.SetLineColor(kBlue);
     298    hSumClear2.Sumw2();
    266299
    267300    // Arrival time spectrum of the extracted pulses
    … …  
    278311
    279312    // Spektrum for the normalized individual spectra
    280     TH1F hSumScale("Sum2", "Signal spectrum of all pixels",  120, 0.05, 12.05);
    281     hSumScale.SetXTitle("pulse integral [pe]");
    282     hSumScale.SetYTitle("Rate");
    283     hSumScale.SetStats(false);
    284     hSumScale.Sumw2();
     313    TH1F hSumScale1("SumScale1", "Signal spectrum of all pixels (incl TM)",  120, 0.05, 12.05);
     314    hSumScale1.SetXTitle("pulse integral [pe]");
     315    hSumScale1.SetYTitle("Rate");
     316    hSumScale1.SetStats(false);
     317    hSumScale1.Sumw2();
     318
     319    TH1F hSumScale2("SumScale2", "Signal spectrum of all pixels (excl TM)",  120, 0.05, 12.05);
     320    hSumScale2.SetXTitle("pulse integral [pe]");
     321    hSumScale2.SetYTitle("Rate");
     322    hSumScale2.SetStats(false);
     323    hSumScale2.Sumw2();
    285324
    286325    // define fit function for Amplitudespectrum
    … …  
    342381    func.SetParLimits(2, 0, 1);            // Sigma
    343382    func.SetParLimits(3, 0, 1);            // Crosstalk
    344     func.SetParLimits(5, 0, 150);          // Noise
     383    if (!fixednoise)
     384        func.SetParLimits(5, 0, 150);          // Noise
    345385
    346386    func.SetParameter(0, ampl);                         // Amplitude
    347387    func.SetParameter(1, maxpos);                       // Gain
    348388    func.SetParameter(2, 0.1);                          // Sigma
    349     func.SetParameter(3, 0.15);                         // Crosstalk
     389    func.SetParameter(3, 0.16);                         // Crosstalk
    350390    func.SetParameter(4, 0*integration_window);         // Baseline
    351     func.SetParameter(5, 1.*sqrt(integration_window));  // Noise
     391    if (fixednoise)
     392        func.FixParameter(5, -1);                       // Noise
     393    else
     394        func.SetParameter(5, 0.1*maxpos);               // Noise
     395
    352396    func.SetParameter(6, 0.4);                          // Expo
    353397
    … …  
    474518        const float fCrosstlk  = xtalk(func);
    475519        const float fOffset    = func.GetParameter(4);
    476         const float fNoise     = func.GetParameter(5)/sqrt(integration_window);
     520        const float fNoise     = func.GetParameter(5)<0 ? fGainRMS*fGain/sqrt(integration_window) : func.GetParameter(5)/sqrt(integration_window);
    477521        const float fCoeffR    = func.GetParameter(6);
    478522
    … …  
    549593
    550594        // Fill Parameters into histograms
    551         hRate.Fill(      fRate);
    552         hGain.Fill(      fGain);
    553         hRelSigma.Fill(  fGainRMS);
    554         hCrosstalk.Fill( fCrosstlk);
    555         hBaseline.Fill(  fOffset/integration_window);
    556         hNoise.Fill(     fNoise);
    557         hChi2.Fill(      chi2);
    558         hNormGain.Fill(  fGain/gain);
    559         hFitProb.Fill(   fit_prob);
    560         hCrosstalkP.Fill(fCrosstalkP);
    561         hCoeffR.Fill(    fCoeffR);
     595        hRate1.Fill(      fRate);
     596        hGain1.Fill(      fGain);
     597        hRelSigma1.Fill(  fGainRMS);
     598        hCrosstalk1.Fill( fCrosstlk);
     599        hBaseline1.Fill(  fOffset/integration_window);
     600        hNoise1.Fill(     fNoise);
     601        hChiSq1.Fill(     chi2);
     602        hNormGain1.Fill(  fGain/gain);
     603        hFitProb1.Fill(   fit_prob);
     604        hCrosstalkP1.Fill(fCrosstalkP);
     605        hCoeffR1.Fill(    fCoeffR);
     606
     607        if (!pmap.index(pixel).isTM())
     608        {
     609            hRate2.Fill(      fRate);
     610            hGain2.Fill(      fGain);
     611            hRelSigma2.Fill(  fGainRMS);
     612            hCrosstalk2.Fill( fCrosstlk);
     613            hBaseline2.Fill(  fOffset/integration_window);
     614            hNoise2.Fill(     fNoise);
     615            hChiSq2.Fill(     chi2);
     616            hNormGain2.Fill(  fGain/gain);
     617            hFitProb2.Fill(   fit_prob);
     618            hCrosstalkP2.Fill(fCrosstalkP);
     619            hCoeffR2.Fill(    fCoeffR);
     620        }
    562621
    563622        // Fill sum spectrum
    564623        for (int b=1; b<=hist->GetNbinsX(); b++)
    565             hSumScale.Fill((hist->GetBinCenter(b)-fOffset)/fGain, hist->GetBinContent(b));
     624            hSumScale1.Fill((hist->GetBinCenter(b)-fOffset)/fGain, hist->GetBinContent(b));
     625
     626        if (!pmap.index(pixel).isTM())
     627            for (int b=1; b<=hist->GetNbinsX(); b++)
     628                hSumScale2.Fill((hist->GetBinCenter(b)-fOffset)/fGain, hist->GetBinContent(b));
     629
    566630        delete hist;
    567631
    568632        // Because of the rebinning...
    569633        hist = hsignal->ProjectionY("proj", pixel+1, pixel+1);
    570         hSumClear.Add(hist);
     634        hSumClear1.Add(hist);
     635        if (!pmap.index(pixel).isTM())
     636            hSumClear2.Add(hist);
    571637        delete hist;
    572638
    … …  
    647713
    648714    canv->cd(1);
    649     hh = hRate.DrawCopy();
     715    hh = hRate1.DrawCopy();
     716    hh = hRate2.DrawCopy("same");
    650717
    651718    canv->cd(2);
    652     hh = hGain.DrawCopy();
     719    hh = hGain1.DrawCopy();
     720    hh = hGain2.DrawCopy("same");
    653721    hh->Fit("gaus");
    654722
    655723    canv->cd(3);
    656     hh = hBaseline.DrawCopy();
     724    hh = hBaseline1.DrawCopy();
     725    hh = hBaseline2.DrawCopy("same");
    657726    hh->Fit("gaus");
    658727
    659728    canv->cd(4);
    660     hh = hRelSigma.DrawCopy();
     729    hh = hRelSigma1.DrawCopy();
     730    hh = hRelSigma2.DrawCopy("same");
    661731    hh->Fit("gaus");
    662732
    663733    canv->cd(5);
    664     hh = hCrosstalk.DrawCopy();
     734    hh = hCrosstalk1.DrawCopy();
     735    hh = hCrosstalk2.DrawCopy("same");
    665736    hh->Fit("gaus");
    666737
    667738    canv->cd(6);
    668     hh = hNoise.DrawCopy();
     739    hh = hNoise1.DrawCopy();
     740    hh = hNoise2.DrawCopy("same");
    669741    hh->Fit("gaus");
    670742
    … …  
    676748    canv->cd(1);
    677749    gPad->SetLogy();
    678     hh = hFitProb.DrawCopy();
     750    hh = hFitProb1.DrawCopy();
     751    hh = hFitProb2.DrawCopy("same");
    679752    hh->Fit("gaus");
    680753
    681754    canv->cd(2);
    682     hChi2.DrawCopy();
     755    hChiSq1.DrawCopy();
     756    hChiSq2.DrawCopy("same");
    683757
    684758    canv->cd(4);
    685     hh = hCoeffR.DrawCopy();
     759    hh = hCoeffR1.DrawCopy();
     760    hh = hCoeffR2.DrawCopy("same");
    686761    hh->Fit("gaus");
    687762
    688763    canv->cd(5);
    689     hh = hCrosstalkP.DrawCopy();
     764    hh = hCrosstalkP1.DrawCopy();
     765    hh = hCrosstalkP2.DrawCopy("same");
    690766    hh->Fit("gaus");
    691767
    … …  
    702778    canv->cd(2);
    703779    gPad->SetLogy();
    704     hh = hNormGain.DrawCopy();
     780    hh = hNormGain1.DrawCopy();
     781    hh = hNormGain2.DrawCopy("same");
    705782    hh->Fit("gaus");
    706783
    … …  
    708785    gROOT->SetSelectedPad(0);
    709786    c11.cd();
    710     hSumClear.DrawCopy("hist same");
     787    hSumClear1.DrawCopy("hist same");
    711788
    712789    //-------------------- fit gain corrected sum spectrum -------------------
    713790
    714791    gROOT->SetSelectedPad(0);
    715     TCanvas &c11b = d->AddTab("CleanHist");
     792    TCanvas &c11b = d->AddTab("CleanHist1");
    716793    c11b.cd();
    717794    gPad->SetLogy();
    … …  
    719796    gPad->SetGridy();
    720797
    721     const Int_t    maxbin1 = hSumClear.GetMaximumBin();
    722     const Double_t ampl1   = hSumClear.GetBinContent(maxbin1);
     798    const Int_t    maxbin1 = hSumClear1.GetMaximumBin();
     799    const Double_t ampl1   = hSumClear1.GetBinContent(maxbin1);
    723800
    724801    func.SetParameters(res_par);
    … …  
    727804    func.ReleaseParameter(6);
    728805
    729     hSumClear.Fit(&func, fast?"LN0QSR":"LIN0QSR");
    730 
    731     hSumClear.DrawCopy();
     806    hSumClear1.Fit(&func, fast?"LN0QSR":"LIN0QSR");
     807
     808    hSumClear1.DrawCopy();
    732809    func.DrawCopy("same");
    733810
    … …  
    739816
    740817    gROOT->SetSelectedPad(0);
    741     TCanvas &c12 = d->AddTab("GainHist");
     818    TCanvas &c11c = d->AddTab("CleanHist2");
     819    c11c.cd();
     820    gPad->SetLogy();
     821    gPad->SetGridx();
     822    gPad->SetGridy();
     823
     824    const Int_t    maxbin1b = hSumClear2.GetMaximumBin();
     825    const Double_t ampl1b   = hSumClear2.GetBinContent(maxbin1b);
     826
     827    func.SetParameters(res_par);
     828    func.SetParLimits(0, 0, 2*ampl1b);
     829    func.SetParameter(0, ampl1b);
     830    func.ReleaseParameter(6);
     831
     832    hSumClear2.Fit(&func, fast?"LN0QSR":"LIN0QSR");
     833
     834    hSumClear2.DrawCopy();
     835    func.DrawCopy("same");
     836
     837    cout << "--------------------" << endl;
     838    PrintFunc(func, integration_window);
     839    cout << "--------------------" << endl;
     840
     841    //-------------------- fit gain corrected sum spectrum -------------------
     842
     843    gROOT->SetSelectedPad(0);
     844    TCanvas &c12 = d->AddTab("GainHist1");
    742845    c12.cd();
    743846    gPad->SetLogy();
    … …  
    745848    gPad->SetGridy();
    746849
    747     const Int_t    maxbin2 = hSumScale.GetMaximumBin();
    748     const Double_t ampl2   = hSumScale.GetBinContent(maxbin2);
     850    const Int_t    maxbin2 = hSumScale1.GetMaximumBin();
     851    const Double_t ampl2   = hSumScale1.GetBinContent(maxbin2);
    749852
    750853    //Set fit parameters
    751854    Double_t par2[7] =
    752855    {
    753         ampl2, 1, 0.1, res_par[3], 0, res_par[5]/res_par[1], res_par[6]
     856        ampl2, 1, 0.1, res_par[3], 0, res_par[5]<0 ? -1 : res_par[5]/res_par[1], res_par[6]
    754857    };
    755858
    … …  
    760863
    761864    func.SetRange(0.62, 9);
    762     hSumScale.Fit(&func, fast?"LN0QSR":"LIN0QSR");
    763 
    764     hSumScale.DrawCopy();
     865    hSumScale1.Fit(&func, fast?"LN0QSR":"LIN0QSR");
     866
     867    hSumScale1.DrawCopy();
     868    func.DrawCopy("same");
     869
     870    cout << "--------------------" << endl;
     871    PrintFunc(func, integration_window);
     872    cout << "--------------------" << endl;
     873
     874    //-------------------- fit gain corrected sum spectrum -------------------
     875
     876    gROOT->SetSelectedPad(0);
     877    TCanvas &c12b = d->AddTab("GainHist2");
     878    c12b.cd();
     879    gPad->SetLogy();
     880    gPad->SetGridx();
     881    gPad->SetGridy();
     882
     883    const Int_t    maxbin2b = hSumScale2.GetMaximumBin();
     884    const Double_t ampl2b   = hSumScale2.GetBinContent(maxbin2b);
     885
     886    //Set fit parameters
     887    Double_t par2b[7] =
     888    {
     889        ampl2b, 1, 0.1, res_par[3], 0, res_par[5]<0 ? -1 : res_par[5]/res_par[1], res_par[6]
     890    };
     891
     892    func.SetParameters(par2b);
     893    func.SetParLimits(0, 0, 2*ampl2b);
     894    func.FixParameter(1, 1);
     895    func.FixParameter(4, 0);
     896
     897    func.SetRange(0.62, 9);
     898    hSumScale2.Fit(&func, fast?"LN0QSR":"LIN0QSR");
     899
     900    hSumScale2.DrawCopy();
    765901    func.DrawCopy("same");
    766902
    … …  
    794930    return 0;
    795931}
    796 
Note: See TracChangeset for help on using the changeset viewer.