Ignore:
Timestamp:
06/25/12 11:19:25 (13 years ago)
Author:
Jens Buss
Message:
add function to fit gaussians to spectrum to determine cross
talk, added list for plotting suspicious pixel, fit
praramter (gain) limitied by pos of 1.Max and fhwm,
sumspectrum fitted  in advance to get single pixel fit start
values from
File:
1 edited

Legend:

Unmodified
Added
Removed

  • fact/tools/marsmacros/singlepe.C

    r14182 r14230  
    88#include <TFitResultPtr.h>
    99#include <TFitResult.h>
     10#include <Getline.h>
     11
     12#include <cstdio>
     13#include <stdio.h>
     14#include <stdint.h>
    1015
    1116#include "MH.h"
    … …  
    407412
    408413Double_t fcn(Double_t *x, Double_t *par);
     414Double_t fcn_cross(Double_t *x, Double_t *par);
     415
     416void FitGaussiansToSpectrum( UInt_t nfunc,  TH1 &hSource,   TH1 &hDest,
     417                             Int_t maxbin,  double fwhm, Double_t gain );
    409418
    410419int singlepe()
    … …  
    412421    // ======================================================
    413422
    414     const char *drsfile = "/fact/raw/2012/06/01/20120601_003.drs.fits.gz";
     423    const char *drsfile = "/fact/raw/2012/03/09/20120309_012.drs.fits.gz";
     424//    const char *drsfile = "/fact/raw/2012/01/13/20120113_003.drs.fits.gz";
    415425
    416426    MDirIter iter;
    417     iter.AddDirectory("/fact/raw/2012/06/01", "20120601_005.fits.gz");
     427//    iter.AddDirectory("/fact/raw/2012/03/09", "20120309_014.fits.gz");
     428//    iter.AddDirectory("/fact/raw/2012/03/09", "20120309_032.fits.gz");
     429//    iter.AddDirectory("/fact/raw/2012/03/09", "20120309_038.fits.gz");
     430    iter.AddDirectory("/fact/raw/2012/03/09", "20120309_017.fits.gz");
     431//    iter.AddDirectory("/fact/raw/2012/03/09", "20120309_018.fits.gz");
     432//    iter.AddDirectory("/fact/raw/2012/01/13", "20120113_007.fits.gz");
    418433
    419434    // ======================================================
    … …  
    566581    // AFTER this the Code for the spektrum fit follows
    567582    // access the spektrumhistogram by the pointer *hist
     583    UInt_t maxOrder     = 6;
     584
    568585
    569586    MGeomCamFACT fact;
    … …  
    571588
    572589    //built an array of Amplitude spectra
    573     TH1F hAmplitude("Rate",      "Average number of dark counts per event",
    574                     100,  0,  10);
    575     TH1F hGain     ("Gain",      "Gain distribution",
    576                      50,  200,   350);
    577     TH1F hGainRMS  ("RelSigma",  "Distribution of rel. Sigma",
    578                     100,   0,   30);
    579     TH1F hCrosstlk ("Crosstalk", "Distribution of crosstalk probability",
    580                     100,   0,    25);
    581     TH1F hOffset   ("Offset",    "Baseline offset distribution",
    582                     50, -7.5,  2.5);
    583     TH1F hFitProb  ("FitProb",   "FitProb distribution",
     590    TH1F hAmplitude ("Rate",      "Average number of dark counts per event",
     591                     200,  0,  20);
     592    TH1F hGain      ("Gain",      "Gain distribution",
     593                     50,  100,   300);
     594    TH1F hGainRMS   ("RelSigma",  "Distribution of rel. Sigma",
     595                     100,   0,   30);
     596    TH1F hCrosstlk  ("Crosstalk", "Distribution of crosstalk probability",
     597                     100,   0,    25);
     598    TH1F hOffset    ("Offset",    "Baseline offset distribution",
     599                     50, -7.5,  2.5);
     600    TH1F hFitProb   ("FitProb",   "FitProb distribution",
    584601                     50,   0,  100);
    585     TH1F hSum1     ("Sum1",      "Sum",
    586                     100,    0,  2000);
    587     TH1F hSum2     ("Sum2",      "Sum",
    588                     100,    0,   10);
     602    TH1F hSum       ("Sum1",      "Sum of all pixel amplitude Spectra",
     603                     100,    0,  2000);
     604    TH1F hSumScale  ("Sum2",
     605                     "Sum of all pixel amplitude Spectra (scaled with gain)",
     606                     100,    0,   10);
    589607
    590608    // define fit function for Amplitudespectrum
    … …  
    594612
    595613    // define fit function for Amplitudespectrum
    596     TF1 func2("sum_spektrum", fcn, 0, 2000, 5);
    597     func2.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset");
     614    TF1 func2("sum_spektrum", fcn_cross, 0, 2000, 6);
     615    func2.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset", "PowOfXtalk");
    598616    func2.SetLineColor(kBlue);
    599617
    … …  
    623641    usePixel[1393] = 0;
    624642
    625     int count_ok = 0;
     643//    usePixel[990] = 0;
     644
     645    // Map for which pixel which were suspicous
     646    bool suspicous[1440];
     647    for (int i=0; i<1440; i++) suspicous[i]=false;
     648
     649    // List of Pixel that should be ignored in camera view
     650    suspicous[802]       = true;
     651    suspicous[543]     = true;
     652    suspicous[465]     = true;
     653    suspicous[1154]     = true;
     654    suspicous[342]     = true;
     655    suspicous[1319]      = true;
     656    suspicous[1401]      = true;
     657    suspicous[1318]      = true;
     658    suspicous[1400]     = true;
     659    suspicous[243]      = true;
     660    suspicous[1299]     = true;
     661    suspicous[764]     = true;
     662    suspicous[762]     = true;
     663    suspicous[1427]      = true;
     664    suspicous[1319]     = true;
     665    suspicous[118]     = true;
     666
     667    //------------------------------------------------------------------------
     668    //  Fill and calculate sum spectrum
    626669
    627670    // Begin of Loop over Pixels
    628671    for (UInt_t pixel = 0; pixel < singles.GetNumPixels(); pixel++)
    629672    {
     673        //jump to next pixel if the current is marked as faulty
     674        if (usePixel[pixel]==0)
     675            continue;
     676
     677        TH1D *hist = hsignal->ProjectionY("proj", pixel+1, pixel+1);
     678        hist->SetDirectory(0);
     679        //loop over slices
     680        for (int b=1; b<=hist->GetNbinsX(); b++)
     681        {
     682            //Fill sum spectrum
     683            hSum.Fill(hist->GetBinCenter(b), hist->GetBinContent(b));
     684        }
     685    }
     686
     687    gROOT->SetSelectedPad(0);
     688    TCanvas &c11 = d->AddTab("SumHist");
     689    gPad->SetLogy();
     690    hSum.DrawCopy();
     691    //------------------------fit sum spectrum-------------------------------
     692    const Int_t    maxbin   = hSum.GetMaximumBin();
     693    const Double_t binwidth = hSum.GetBinWidth(maxbin);
     694    const Double_t ampl     = hSum.GetBinContent(maxbin);
     695
     696    double fwhm2 = 0;
     697
     698    //Calculate full width half Maximum
     699    for (int i=1; i<maxbin; i++)
     700        if (hSum.GetBinContent(maxbin-i)+hSum.GetBinContent(maxbin+i)<ampl*2)
     701        {
     702            fwhm2 = (2*i+1)*binwidth;
     703            break;
     704        }
     705    if (fwhm2==0)
     706    {
     707        cout << "Could not determine start value for sigma." << endl;
     708    }
     709
     710    //Set fit parameters
     711    Double_t sum_par[6] =
     712    {
     713        ampl, hSum.GetBinCenter(maxbin), fwhm2*1.4, 0.1, -10, 1
     714    };
     715    func2.SetParameters(sum_par);
     716
     717    //calculate fit range
     718    const double min = sum_par[1]-fwhm2*1.4;
     719    const double max = sum_par[1]*(maxOrder+1);
     720    func2.SetRange(min-2*fwhm2, max);
     721
     722    //Fit and draw spectrum
     723    hSum.Fit(&func2, "NOS", "", min, max);
     724    func2.DrawCopy("SAME");
     725    func2.GetParameters(sum_par);
     726
     727    //Set Parameters for following fits
     728    const float  Amplitude      = sum_par[0];
     729    const float  gain           = sum_par[1];
     730    const float  GainRMS        = sum_par[2];
     731    const float  Crosstlk       = sum_par[3];
     732    const float  Offset         = sum_par[4];
     733    const float  PowCrosstlk    = sum_par[5];
     734
     735
     736    //--------fit gausses to peaks of sum specetrum -----------
     737
     738    // create histo for height of Maximum amplitudes vs. pe
     739    double min_bin  =   hSum.GetBinCenter(maxbin);
     740    double max_bin  =   hSum.GetBinCenter(maxOrder*(maxbin));
     741    double bin_width=   (max_bin - min_bin)/10;
     742
     743    TH1D hMaxHeightsSum("MaxHeightSum",      "peak heights",
     744                        10,  min_bin-bin_width, max_bin*2-bin_width/2 );
     745
     746    FitGaussiansToSpectrum
     747            (
     748                maxOrder,
     749                hSum,
     750                hMaxHeightsSum,
     751                maxbin,
     752                fwhm2,
     753                gain
     754                );
     755
     756    //EOF: fit sum spectrum-----------------------------
     757
     758    hMaxHeightsSum.SetLineColor(kRed);
     759//    hMaxHeightsSum.DrawCopy("SAME");
     760
     761    //Fit Peak heights
     762    TF1 fExpo1( "fExpo1","expo", min, max);
     763    fExpo1.SetLineColor(kRed);
     764    hMaxHeightsSum.Fit(&fExpo1, "N0S" );
     765    fExpo1.DrawCopy("SAME");
     766
     767    //  EOF:    Fill and calculate sum spectrum
     768    //------------------------------------------------------------------------
     769
     770
     771
     772    //------------------------------------------------------------------------
     773    //  Gain Calculation for each Pixel
     774
     775    int count_ok = 0;
     776    // Begin of Loop over Pixels
     777    for (UInt_t pixel = 0; pixel < singles.GetNumPixels(); pixel++)
     778    {
     779
    630780        if (usePixel[pixel]==0)
    631781            continue;
    … …  
    638788        const Double_t binwidth = hist->GetBinWidth(maxbin);
    639789        const Double_t ampl     = hist->GetBinContent(maxbin);
     790        const Double_t maxPos   = hist->GetBinCenter(maxbin);
    640791
    641792        double fwhm = 0;
    … …  
    657808        Double_t par[5] =
    658809        {
    659             ampl, hist->GetBinCenter(maxbin), fwhm*1.4, 0.1, -10
     810            ampl,
     811            hist->GetBinCenter(maxbin),
     812//            gain,
     813            fwhm*2,
     814//            GainRMS,
     815            Crosstlk,
     816            Offset
    660817        };
    661818
     819
    662820        func.SetParameters(par);
    663 
    664         const double min = par[1]-fwhm*1.4;
    665         const double max = par[1]*5;
    666 
     821        const double fit_min = par[1]-fwhm*1.4;
     822        const double fit_max = par[1]*maxOrder;
     823        func.SetRange(fit_min-fwhm*2, fit_max);
     824        func.SetParLimits(1, maxPos - fwhm, maxPos + fwhm);
     825
     826        if (suspicous[pixel] == true)
     827        {
     828            cout << "Maxbin\t" << maxbin << endl;
     829            cout << "min\t" << fit_min << endl;
     830            cout << "max\t" << fit_max << endl;
     831            cout << "Amplitude, 1.Max, Sigma, fwhm:\t"
     832                 << par[0] << "\t"
     833                 << par[1] << "\t"
     834                 << par[2] << "\t"
     835                 << fwhm << "\t" << endl;
     836        }
    667837        //Fit Pixels spectrum
    668         const TFitResultPtr rc = hist->Fit(&func, "N0QS", "", min, max);
     838        const TFitResultPtr rc = hist->Fit(&func, "N0QS", "", fit_min, fit_max);
    669839
    670840        const bool ok = int(rc)==0;
    … …  
    676846        const float  fCrosstlk  = func.GetParameter(3);
    677847        const float  fOffset    = func.GetParameter(4)/singles.GetIntegrationWindow();
     848//        const float  fCrossOrd  = func.GetParameter(5);
     849
     850        if (suspicous[pixel] == true)
     851        {
     852            cout << "Amplitude, 1.Max, Sigma, fwhm:\t"
     853                 << fAmplitude << "\t"
     854                 << fGain << "\t"
     855                 << f2GainRMS << "\t"
     856                 << fwhm << "\t" << endl;
     857        }
    678858
    679859        hAmplitude.Fill(fAmplitude);
    … …  
    695875        }
    696876
     877        //Fill sum spectrum
    697878        for (int b=1; b<=hist->GetNbinsX(); b++)
    698879        {
    699             hSum1.Fill(hist->GetBinCenter(b), hist->GetBinContent(b));
    700             hSum2.Fill((hist->GetBinCenter(b)-fOffset)/fGain, hist->GetBinContent(b));
    701         }
    702 
    703         if (count_ok==0)
     880            hSum.Fill(hist->GetBinCenter(b), hist->GetBinContent(b));
     881            hSumScale.Fill((hist->GetBinCenter(b)-fOffset)/fGain, hist->GetBinContent(b));
     882        }
     883
     884        //plott special pixel
     885        if (
     886                count_ok==0 ||
     887                suspicous[pixel]
     888                )
    704889        {
    705890            TCanvas &c = d->AddTab(Form("Pix%d", pixel));
    … …  
    707892            hist->SetName(Form("Pix%d", pixel));
    708893            hist->DrawCopy()->SetDirectory(0);
     894//            hist->Draw();
    709895            func.DrawCopy("SAME");
    710896        }
    711897
    712 
    713 
    714898        count_ok++;
    715899
    716900        delete hist;
    717     }
     901//        if (suspicous[pixel] == true) usePixel[pixel]=0;
     902    }
     903    // End of Loop over Pixel
    718904
    719905    //------------------fill histograms-----------------------
    … …  
    726912    hcam.DrawCopy();
    727913
    728     gROOT->SetSelectedPad(0);
    729     TCanvas &c11 = d->AddTab("SumHist");
    730     gPad->SetLogy();
    731     hSum1.DrawCopy();
    732 
    733     //------------------------fit sum spectrum-------------------------------
    734     const Int_t    maxbin   = hSum1.GetMaximumBin();
    735     const Double_t binwidth = hSum1.GetBinWidth(maxbin);
    736     const Double_t ampl     = hSum1.GetBinContent(maxbin);
    737 
    738     double fwhm2 = 0;
    739 
    740     //Calculate full width half Maximum
    741     for (int i=1; i<maxbin; i++)
    742         if (hSum1.GetBinContent(maxbin-i)+hSum1.GetBinContent(maxbin+i)<ampl*2)
    743         {
    744             fwhm2 = (2*i+1)*binwidth;
    745             break;
    746         }
    747     if (fwhm2==0)
    748     {
    749         cout << "Could not determine start value for sigma." << endl;
    750     }
    751 
    752     //Set fit parameters
    753     Double_t par[5] =
    754     {
    755         ampl, hSum1.GetBinCenter(maxbin), fwhm2*1.4, 0.1, -10
    756     };
    757     func2.SetParameters(par);
    758 
    759     //calculate fit range
    760     const double min = par[1]-fwhm2*1.4;
    761     const double max = par[1]*5;
    762 
    763     //Fit and draw spectrum
    764     hSum1.Fit(&func2, "NOQS", "", min, max);
    765     func2.DrawCopy("SAME");
    766 
    767     //-------------------END OF: fit sum spectrum-----------------------------
    768 
    769     // Draw gain corrected sum specetrum
     914    //------------------ Draw gain corrected sum specetrum -------------------
    770915    gROOT->SetSelectedPad(0);
    771916    TCanvas &c12 = d->AddTab("GainHist");
    772917    gPad->SetLogy();
    773     hSum2.DrawCopy();
    774     const Double_t fMaxPos  = hSum2.GetBinCenter(maxbin);
    775 
    776     //--------fit gausses to peaks of gain corrected sum specetrum -----------
    777     UInt_t nfunc    = 5;
    778 
    779     // create array for gaus funtions to fit to peaks of spectrum
    780     TF1 **gaus      = new TF1*[nfunc];
    781 
    782     // create histo for height of Maximum amplitudes vs. pe
    783     TH1D hMaxHeights("MaxHeights",      "peak heights",
    784                       20,  0,   20);
    785 
    786     // fit gauss functions to spectrum
    787     for (UInt_t nr = 0; nr<nfunc; nr++)
    788     {
    789         char fname[20];
    790         sprintf(fname,"gaus%d",nr);
    791 
    792         //create gauss functions
    793         gaus[nr] = new TF1( fname,"gaus",
    794                             fMaxPos*(nr+1)-0.6,     fMaxPos*(nr+1)+0.6);
    795         //fit and draw gauss functions
    796         hSum2.Fit(  gaus[nr],   "N0QS",     "",
    797                     fMaxPos*(nr+1)-0.6,     fMaxPos*(nr+1)+0.6);
    798         gaus[nr]->DrawCopy("SAME");
    799 
    800         //fill histo for height of Maximum amplitudes vs. pe
    801         hMaxHeights.SetBinContent(nr, gaus[nr]->GetParameter(0));
    802         delete gaus[nr];
    803     }
    804     delete[] gaus;
    805 
    806     //------------------fit gain corrected sum spectrum-----------------------
    807 
    808     const Int_t    maxbin2   = hSum2.GetMaximumBin();
    809     const Double_t binwidth2 = hSum2.GetBinWidth(maxbin);
    810     const Double_t ampl2     = hSum2.GetBinContent(maxbin);
     918    hSumScale.DrawCopy();
     919
     920    //-------------------- fit gain corrected sum spectrum -------------------
     921    const Int_t    maxbin2   = hSumScale.GetMaximumBin();
     922    const Double_t binwidth2 = hSumScale.GetBinWidth(maxbin2);
     923    const Double_t ampl2     = hSumScale.GetBinContent(maxbin2);
    811924
    812925    //Calculate full width half Maximum
    813926    fwhm2 = 0;
    814927    for (int i=1; i<maxbin; i++)
    815         if (hSum2.GetBinContent(maxbin2-i)+hSum2.GetBinContent(maxbin2+i)<ampl2*2)
     928        if (hSumScale.GetBinContent(maxbin2-i)+hSumScale.GetBinContent(maxbin2+i)<ampl2*2)
    816929        {
    817930            fwhm2 = (2*i+1)*binwidth2;
    … …  
    826939    Double_t par2[5] =
    827940    {
    828         ampl2, hSum2.GetBinCenter(maxbin2), fwhm2*1.4, 0.1, 0
     941        ampl2, hSumScale.GetBinCenter(maxbin2), fwhm2*1.4, 0.1, -0.2
    829942    };
    830943
    … …  
    832945
    833946    const double min2 = par2[1]-fwhm2*1.4;
    834     const double max2 = par2[1]*5;
    835 
    836     hSum2.Fit(&func3, "N0QS", "", min2, max2);
     947    const double max2 = par2[1]*maxOrder;
     948//    func.SetRange(min2-fwhm2, max);
     949    hSumScale.Fit(&func3, "N0QS", "", min2, max2);
    837950    func3.DrawCopy("SAME");
    838     //-------------------END OF: fit sum spectrum-----------------------------
    839 
    840     TCanvas &c15 = d->AddTab("PeakHeights");
    841     gPad->SetLogy();
    842     hMaxHeights.DrawCopy();
     951
     952    //--------fit gausses to peaks of gain corrected sum specetrum -----------
     953
     954    // create histo for height of Maximum amplitudes vs. pe
     955    TH1D hMaxHeights("MaxHeights",      "peak heights",
     956                      10,  hSumScale.GetBinCenter(maxbin-1)-0.5,   10+hSumScale.GetBinCenter(maxbin-1)-0.5);
     957
     958    FitGaussiansToSpectrum
     959            (
     960                maxOrder,
     961                hSumScale,
     962                hMaxHeights,
     963                maxbin2,
     964                fwhm2,
     965                1
     966                );
     967
     968
     969
     970//    TCanvas &c16 = d->AddTab("PeakHeights");
     971//    gPad->SetLogy();
     972    hMaxHeights.SetLineColor(kRed);
     973//    hMaxHeights.DrawCopy("SAME");
     974    TF1 fExpo( "fExpo","expo", 0,  maxOrder-1);
     975    fExpo.SetLineColor(kRed);
     976    hMaxHeights.Fit(&fExpo, "N0QS" );
     977    fExpo.DrawCopy("SAME");
     978//    c12.cd();
     979//    fExpo.DrawCopy("SAME");
    843980
    844981    TCanvas &c2 = d->AddTab("Result");
    … …  
    8801017    Double_t x     = xx[0];
    8811018
    882     Double_t ampl  = par[0];
    883     Double_t gain  = par[1];
    884     Double_t sigma = par[2];
    885     Double_t cross = par[3];
    886     Double_t shift = par[4];
     1019    Double_t ampl   = par[0];
     1020    Double_t gain   = par[1];
     1021    Double_t sigma  = par[2];
     1022    Double_t cross  = par[3];
     1023    Double_t shift  = par[4];
     1024//    Double_t k      = par[5];
    8871025
    8881026    Double_t xTalk = 1;
    … …  
    8971035        y += xTalk*gauss;
    8981036
    899         xTalk *= cross;
     1037//        for (int j = 1; j < k; j++)
     1038//        {
     1039            xTalk *= cross;
     1040//        }
    9001041    }
    9011042
    9021043    return ampl*y;
    9031044}
     1045
     1046Double_t fcn_cross(Double_t *xx, Double_t *par)
     1047{
     1048    Double_t x     = xx[0];
     1049
     1050    Double_t ampl   = par[0];
     1051    Double_t gain   = par[1];
     1052    Double_t sigma  = par[2];
     1053    Double_t cross  = par[3];
     1054    Double_t shift  = par[4];
     1055    Double_t powOfX = par[5];
     1056
     1057    Double_t xTalk = 1;
     1058
     1059    Double_t y = 0;
     1060    for (int order = 1; order < 14; order++)
     1061    {
     1062        Double_t arg   = (x - order*gain - shift)/sigma;
     1063
     1064        Double_t gauss = TMath::Exp(-0.5 * arg * arg/order);
     1065
     1066        y += xTalk*gauss;
     1067
     1068//        for (int j = 1; j < powOfX; j++)
     1069//        {
     1070//            xTalk *= cross;
     1071//        }
     1072        cross *= TMath::Exp(-powOfX*cross);
     1073        xTalk *= cross;
     1074    }
     1075
     1076    return ampl*y;
     1077}
     1078
     1079
     1080void
     1081FitGaussiansToSpectrum(
     1082        UInt_t      maxOrder,
     1083        TH1        &hSource,
     1084        TH1        &hDest,
     1085        Int_t       maxbin,
     1086        double      fwhm,
     1087        Double_t    gain
     1088        )
     1089{
     1090
     1091    Double_t    peakposition = hSource.GetBinCenter(maxbin);
     1092    char        fname[20];
     1093
     1094    // fit gauss functions to spectrum
     1095    for (UInt_t nr = 1; nr<maxOrder; nr++)
     1096    {
     1097        sprintf(fname,"gaus%d",nr);
     1098
     1099        //create gauss functions
     1100        TF1 gaus( fname,"gaus", peakposition-fwhm,     peakposition+fwhm);
     1101        gaus.SetParLimits(1, peakposition-fwhm, peakposition+fwhm);
     1102        gaus.SetParLimits(2, -fwhm, fwhm );
     1103        //fit and draw gauss functions
     1104        hSource.Fit(  &gaus,   "N0QS", "", peakposition-fwhm, peakposition+fwhm);
     1105//        gaus.DrawCopy("SAME");
     1106
     1107        peakposition = gaus.GetParameter(1);
     1108
     1109        //fill histo for height of Maximum amplitudes vs. pe
     1110        hDest.SetBinContent(nr, gaus.GetParameter(0));
     1111
     1112        peakposition += gain;
     1113    }
     1114    return;
     1115}
     1116
Note: See TracChangeset for help on using the changeset viewer.