Ignore:
Timestamp:
06/27/12 17:36:10 (12 years ago)
Author:
Jens Buss
Message:
titles for histograms and axes editted, fit-method reset to ch2, bug fix
in fit function, normalized gain distribution
File:
1 edited

Legend:

Unmodified
Added
Removed

  • fact/tools/marsmacros/singlepe.C

    r14241 r14242  
    153153
    154154        fSignal.SetName("Signal");
    155         fSignal.SetTitle("Title");
    156         fSignal.SetXTitle("X title");
    157         fSignal.SetYTitle("Y title");
     155        fSignal.SetTitle("pulse integral spectrum");
     156        fSignal.SetXTitle("pixel [SoftID]");
     157        fSignal.SetYTitle("time [a.u.]");
    158158        fSignal.SetDirectory(NULL);
    159159
    160160        fTime.SetName("Time");
    161         fTime.SetTitle("Title");
    162         fTime.SetXTitle("X title");
    163         fTime.SetYTitle("Y title");
     161        fTime.SetTitle("arival time spectrum");
     162        fTime.SetXTitle("pixel [SoftID]");
     163        fTime.SetYTitle("time [a.u.]");
    164164        fTime.SetDirectory(NULL);
    165165
    166166        fPulse.SetName("Pulse");
    167         fPulse.SetTitle("Title");
    168         fPulse.SetXTitle("X title");
    169         fPulse.SetYTitle("Y title");
     167        fPulse.SetTitle("average pulse shape spectrum");
     168        fPulse.SetXTitle("pixel [SoftID]");
     169        fPulse.SetYTitle("time [a.u.]");
    170170        fPulse.SetDirectory(NULL);
    171171    }
    … …  
    417417                             Int_t maxbin,  double fwhm, Double_t gain );
    418418
    419 int singlepe(const char *runfile, const char *drsfile, const char *outfile)
     419Double_t MedianOfH1 (
     420        TH1*            inputHisto
     421        );
     422
     423int singlepe(
     424//        const char *runfile, const char *drsfile, const char *outfile
     425        )
    420426{
     427
    421428    // ======================================================
    422429
    423     //const char *drsfile = "/fact/raw/2012/05/18/20120518_012.drs.fits.gz";
     430    const char *drsfile = "/fact/raw/2012/05/18/20120518_003.drs.fits.gz";
     431//    const char *drsfile = "/fact/raw/2012/06/25/20120625_112.drs.fits.gz";
    424432
    425433    MDirIter iter;
    426     //iter.AddDirectory("/fact/raw/2012/05/18", "20120518_005.fits.gz");
    427     iter.AddDirectory(gSystem->DirName(runfile), gSystem->BaseName(runfile));
     434    iter.AddDirectory("/fact/raw/2012/05/18", "20120518_005.fits.gz");
     435//    iter.AddDirectory(gSystem->DirName(runfile), gSystem->BaseName(runfile));
     436
     437//    TString filename;
     438//    for (UInt_t fileNr = 116; fileNr < 147; fileNr++)
     439//    {
     440//        filename    = "20120625_";
     441//        filename   += fileNr;
     442//        filename   += ".fits.gz";
     443//        iter.AddDirectory("/fact/raw/2012/06/25", filename);
     444//    }
     445
     446//     iter.AddDirectory("/fact/raw/2012/06/25", "20120625_116.fits.gz");
     447//     iter.AddDirectory("/fact/raw/2012/06/25", "20120625_117.fits.gz");
     448//     iter.AddDirectory("/fact/raw/2012/06/25", "20120625_118.fits.gz");
     449//     iter.AddDirectory("/fact/raw/2012/06/25", "20120625_119.fits.gz");
    428450
    429451    // ======================================================
    … …  
    570592    d->AddTab("Time");
    571593    TH1 *ht = htime->ProjectionY();
     594    ht->SetYTitle("Counts [a.u.]");
    572595    ht->Scale(1./1440);
    573596    ht->Draw();
    … …  
    575598    d->AddTab("Pulse");
    576599    TH1 *hp = hpulse->ProjectionY();
     600    hp->SetYTitle("Counts [a.u.]");
    577601    hp->Scale(1./1440);
    578602    hp->Draw();
    … …  
    582606    // AFTER this the Code for the spektrum fit follows
    583607    // access the spektrumhistogram by the pointer *hist
    584     UInt_t maxOrder     = 6;
     608    UInt_t maxOrder     = 8;
    585609
    586610
    587611    MGeomCamFACT fact;
    588612    MHCamera hcam(fact);
     613    MHCamera hcam2(fact);
    589614
    590615    //built an array of Amplitude spectra
    591616    TH1F hAmplitude ("Rate",      "Average number of dark counts per event",
    592617                     200,  0,  20);
     618    hAmplitude.SetXTitle("Amplitude [a.u.]");
     619    hAmplitude.SetYTitle("Rate");
     620
    593621    TH1F hGain      ("Gain",      "Gain distribution",
    594622                     50,  100,   300);
     623    hGain.SetXTitle("gain [a.u.]");
     624    hGain.SetYTitle("Rate");
     625
    595626    TH1F hGainRMS   ("RelSigma",  "Distribution of rel. Sigma",
    596627                     100,   0,   30);
     628    hGainRMS.SetXTitle("gainRMS [a.u.]");
     629    hGainRMS.SetYTitle("Rate");
     630    hGainRMS.SetStats(false);
     631
    597632    TH1F hCrosstlk  ("Crosstalk", "Distribution of crosstalk probability",
    598633                     100,   0,    25);
     634    hCrosstlk.SetXTitle("crosstalk [a.u.]");
     635    hCrosstlk.SetYTitle("Rate");
     636
    599637    TH1F hOffset    ("Offset",    "Baseline offset distribution",
    600638                     50, -7.5,  2.5);
     639    hOffset.SetXTitle("baseline offset [a.u.]");
     640    hOffset.SetYTitle("Rate");
     641
    601642    TH1F hFitProb   ("FitProb",   "FitProb distribution",
    602643                     50,   0,  100);
    603     TH1F hSum       ("Sum1",      "Sum of all pixel amplitude Spectra",
     644    hFitProb.SetXTitle("FitProb p-value [a.u.]");
     645    hFitProb.SetYTitle("Rate");
     646    hFitProb.SetStats(false);
     647
     648
     649    TH1F hSum       ("Sum1",      "Sum of all pixel's' integral spectra",
    604650                     100,    0,  2000);
     651    hSum.SetXTitle("pulse integral [a.u.]");
     652    hSum.SetYTitle("Rate");
     653    hSum.SetStats(false);
    605654    hSum.Sumw2();
    606655
     656
    607657    TH1F hSumScale  ("Sum2",
    608                      "Sum of all pixel amplitude Spectra (scaled with gain)",
     658                     "Sum of all pixel's' integral spectra (scaled with gain)",
    609659                     100,    0.05,   10.05);
     660    hSumScale.SetXTitle("pulse integral [pe]");
     661    hSumScale.SetYTitle("Rate");
     662    hSumScale.SetStats(false);
    610663    hSumScale.Sumw2();
    611664
    … …  
    616669
    617670    // define fit function for Amplitudespectrum
    618     TF1 func2("sum_spektrum", fcn_cross, 0, 2000, 6);
    619     func2.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset", "PowOfXtalk");
    620     func2.SetLineColor(kBlue);
     671    TF1 funcSum("sum_spektrum", fcn, 0, 2000, 5);
     672    funcSum.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset", "PowOfXtalk");
     673    funcSum.SetLineColor(kBlue);
    621674
    622675    // define fit function for Amplitudespectrum
    623     TF1 func3("gain_sum_spektrum", fcn, 0, 10, 5);
    624     func3.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset");
    625     func3.SetLineColor(kBlue);
     676    TF1 funcScaled("gain_sum_spektrum", fcn, 0, 10, 5);
     677    funcScaled.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset");
     678    funcScaled.SetLineColor(kBlue);
    626679
    627680
    … …  
    673726        }
    674727    }
    675     for(UInt_t bin = 1; bin < hSum.GetNbinsX()+1; bin++)
     728    for(UInt_t bin = 1; bin < (UInt_t)hSum.GetNbinsX()+1; bin++)
    676729    {
    677730        hSum.SetBinError(bin, hSum.GetBinContent(bin)*0.1);
    … …  
    680733    gROOT->SetSelectedPad(0);
    681734    TCanvas &c11 = d->AddTab("SumHist");
     735    c11.cd();
    682736    gPad->SetLogy();
    683737    gPad->SetGridx();
    … …  
    689743    const Double_t ampl     = hSum.GetBinContent(maxbin);
    690744
    691     double fwhm2 = 0;
     745    double fwhmSum = 0;
    692746
    693747    //Calculate full width half Maximum
    … …  
    695749        if (hSum.GetBinContent(maxbin-i)+hSum.GetBinContent(maxbin+i)<ampl*2)
    696750        {
    697             fwhm2 = (2*i+1)*binwidth;
     751            fwhmSum = (2*i+1)*binwidth;
    698752            break;
    699753        }
    700     if (fwhm2==0)
     754    if (fwhmSum==0)
    701755    {
    702756        cout << "Could not determine start value for sigma." << endl;
    … …  
    706760    Double_t sum_par[6] =
    707761    {
    708         ampl, hSum.GetBinCenter(maxbin), fwhm2*1.4, 0.1, -10, 1
     762        ampl, hSum.GetBinCenter(maxbin), fwhmSum*1.4, 0.1, -10, 1
    709763    };
    710     func2.SetParameters(sum_par);
     764    funcSum.SetParameters(sum_par);
    711765
    712766    //calculate fit range
    713     const double min = sum_par[1]-fwhm2*2;
    714     const double max = sum_par[1]*(maxOrder+1);
    715     func2.SetRange(min, max);
     767    const double min = sum_par[1]-fwhmSum*2;
     768    const double max = sum_par[1]*(maxOrder);
     769    funcSum.SetRange(min, max);
     770    funcSum.FixParameter(5,0.1);
    716771
    717772    //Fit and draw spectrum
    718     hSum.Fit(&func2, "NOS", "", min, max);
    719     func2.DrawCopy("SAME");
    720     func2.GetParameters(sum_par);
     773    hSum.Fit(&funcSum, "NOQS", "", min, max);
     774    funcSum.DrawCopy("SAME");
     775    funcSum.GetParameters(sum_par);
    721776
    722777    //Set Parameters for following fits
    723     const float  Amplitude      = sum_par[0];
     778//    const float  Amplitude      = sum_par[0];
    724779    const float  gain           = sum_par[1];
    725     const float  GainRMS        = sum_par[2];
     780//    const float  GainRMS        = sum_par[2];
    726781    const float  Crosstlk       = sum_par[3];
    727782    const float  Offset         = sum_par[4];
    728     const float  PowCrosstlk    = sum_par[5];
     783//    const float  PowCrosstlk    = sum_par[5];
    729784
    730785
    … …  
    733788    // create histo for height of Maximum amplitudes vs. pe
    734789    double min_bin  =   hSum.GetBinCenter(maxbin);
    735     double max_bin  =   hSum.GetBinCenter(maxOrder*(maxbin));
     790    double max_bin  =   hSum.GetBinCenter((maxOrder-2)*(maxbin));
    736791    double bin_width=   (max_bin - min_bin)/10;
    737792
    … …  
    741796    FitGaussiansToSpectrum
    742797            (
    743                 maxOrder,
     798                maxOrder-2,
    744799                hSum,
    745800                hMaxHeightsSum,
    746801                maxbin,
    747                 fwhm2,
     802                fwhmSum,
    748803                gain
    749804                );
    … …  
    757812    TF1 fExpo1( "fExpo1","expo", min, max);
    758813    fExpo1.SetLineColor(kRed);
    759     hMaxHeightsSum.Fit(&fExpo1, "WLN0S" );
    760     hMaxHeightsSum.DrawCopy("SAME");
    761     fExpo1.DrawCopy("SAME");
     814    hMaxHeightsSum.Fit(&fExpo1, "N0QS" );
     815//    hMaxHeightsSum.DrawCopy("SAME");
     816//    fExpo1.DrawCopy("SAME");
    762817
    763818    //  EOF:    Fill and calculate sum spectrum
    … …  
    781836        hist->SetDirectory(0);
    782837
    783         for(UInt_t bin = 1; bin < hSum.GetNbinsX()+1; bin++)
     838        for(UInt_t bin = 1; bin < (UInt_t)hSum.GetNbinsX()+1; bin++)
    784839        {
    785840            hist->SetBinError(bin, hSum.GetBinContent(bin)*0.1);
    … …  
    841896
    842897        //Fit Pixels spectrum
    843         const TFitResultPtr rc = hist->Fit(&func, "WLN0QS", "", fit_min, fit_max);
     898        const TFitResultPtr rc = hist->Fit(&func, "N0QS", "", fit_min, fit_max);
    844899
    845900        const bool ok = int(rc)==0;
    … …  
    853908//        const float  fCrossOrd  = func.GetParameter(5);
    854909
    855         if (suspicous[pixel] == true)
    856         {
    857             cout << "Amplitude, 1.Max, Sigma, fwhm:\t"
    858                  << fAmplitude << "\t"
    859                  << fGain << "\t"
    860                  << f2GainRMS << "\t"
    861                  << fwhm << "\t" << endl;
    862         }
    863 
    864910        hAmplitude.Fill(fAmplitude);
    865911        hGain.Fill(fGain);
    … …  
    871917
    872918        hcam.SetBinContent(pixel+1, fGain);
     919        hcam2.SetBinContent(pixel+1, 100*f2GainRMS/fGain);
    873920
    874921        usePixel[pixel] = ok;
    … …  
    884931        {
    885932//            hSum.Fill(hist->GetBinCenter(b), hist->GetBinContent(b));
    886             hSumScale.Fill((hist->GetBinCenter(b)-fOffset)/fGain, hist->GetBinContent(b));
     933            hSumScale.Fill((hist->GetBinCenter(b)-fOffset)/fGain, (hist->GetBinContent(b)));
    887934        }
    888935
    … …  
    895942            TCanvas &c = d->AddTab(Form("Pix%d", pixel));
    896943            c.cd();
     944            gPad->SetLogy();
     945            gPad->SetGridx();
     946            gPad->SetGridy();
     947
     948            hist->SetStats(false);
     949            hist->SetYTitle("Counts [a.u.]");
    897950            hist->SetName(Form("Pix%d", pixel));
    898951            hist->DrawCopy("hist")->SetDirectory(0);
    … …  
    911964
    912965    hcam.SetUsed(usePixel);
     966    hcam2.SetUsed(usePixel);
    913967
    914968    TCanvas &canv = d->AddTab("Gain");
    915969    canv.cd();
     970    canv.Divide(2,2);
     971
     972    canv.cd(1);
    916973    hcam.SetNameTitle( "gain","Gain distribution over whole camera");
    917974    hcam.DrawCopy();
     975
     976    canv.cd(2);
     977    hcam2.SetNameTitle( "gainRMS","GainRMS distribution over whole camera");
     978    hcam2.DrawCopy();
     979
     980    TCanvas &cam_canv = d->AddTab("Camera_Gain");
    918981
    919982    //------------------ Draw gain corrected sum specetrum -------------------
    920983    gROOT->SetSelectedPad(0);
    921984    TCanvas &c12 = d->AddTab("GainHist");
     985    c12.cd();
    922986    gPad->SetLogy();
    923987    gPad->SetGridx();
    924988    gPad->SetGridy();
    925     for(UInt_t bin = 1; bin < hSum.GetNbinsX()+1; bin++)
     989
     990    for(UInt_t bin = 1; bin < (UInt_t)hSum.GetNbinsX()+1; bin++)
    926991    {
    927992        hSumScale.SetBinError(bin, hSum.GetBinContent(bin)*0.1);
    … …  
    9361001
    9371002    //Calculate full width half Maximum
    938     fwhm2 = 0;
     1003    double fwhmScaled = 0;
    9391004    for (int i=1; i<maxbin; i++)
    9401005        if (hSumScale.GetBinContent(maxbin2-i)+hSumScale.GetBinContent(maxbin2+i)<ampl2*2)
    9411006        {
    942             fwhm2 = (2*i+1)*binwidth2;
     1007            fwhmScaled = (2*i+1)*binwidth2;
    9431008            break;
    9441009        }
    945     if (fwhm2==0)
     1010    if (fwhmScaled==0)
    9461011    {
    9471012        cout << "Could not determine start value for sigma." << endl;
    … …  
    9491014
    9501015    //Set fit parameters
    951     Double_t par2[5] =
    952     {
    953         ampl2, 1, fwhm2*1.4, Crosstlk, 0
     1016    Double_t par2[6] =
     1017    {
     1018        ampl2, 0.9, fwhmScaled*1.4, Crosstlk, 0, 1.1
    9541019    };
    9551020
    956     func3.SetParameters(par2);
    957     func3.FixParameter(1,1);
    958     func3.FixParameter(4,0);
    959 
    960     const double min2 = par2[1]-fwhm2*1.4;
    961     const double max2 = par2[1]*maxOrder;
    962 //    func.SetRange(min2-fwhm2, max);
    963     hSumScale.Fit(&func3, "WLN0QS", "", min2, max2);
    964     func3.DrawCopy("SAME");
    965 
     1021    funcScaled.SetParameters(par2);
     1022//    funcScaled.FixParameter(1,0.9);
     1023    funcScaled.FixParameter(4,0);
     1024    funcScaled.FixParameter(5,Crosstlk);
     1025
     1026    const double minScaled = par2[1]-fwhmScaled;
     1027    const double maxScaled = par2[1]*maxOrder;
     1028    cout << "par2[1]" << par2[1] << endl;
     1029    cout << "maxScaled\t"           << maxScaled
     1030         << " \t\t minScaled\t"     << minScaled << endl;
     1031
     1032    funcScaled.SetRange(minScaled, maxScaled);
     1033    hSumScale.Fit(&funcScaled, "N0QS", "", minScaled, maxScaled);
     1034    funcScaled.DrawCopy("SAME");
    9661035    //--------fit gausses to peaks of gain corrected sum specetrum -----------
    9671036
    … …  
    9721041    FitGaussiansToSpectrum
    9731042            (
    974                 maxOrder,
     1043                maxOrder-2,
    9751044                hSumScale,
    9761045                hMaxHeights,
    9771046                maxbin2,
    978                 fwhm2,
     1047                fwhmScaled,
    9791048                1
    9801049                );
    … …  
    9861055    hMaxHeights.SetLineColor(kRed);
    9871056//    hMaxHeights.DrawCopy("SAME");
    988     TF1 fExpo( "fExpo","expo", 0,  maxOrder-1);
     1057    TF1 fExpo( "fExpo","expo", 0,  maxOrder-2);
    9891058    fExpo.SetLineColor(kRed);
    9901059    hMaxHeights.Fit(&fExpo, "N0QS" );
    991     fExpo.DrawCopy("SAME");
     1060//    fExpo.DrawCopy("SAME");
    9921061//    c12.cd();
    9931062//    fExpo.DrawCopy("SAME");
    … …  
    10011070    c2.cd(2);
    10021071    gPad->SetLogy();
    1003     TF1 GainGaus( "GainGaus", "gaus");
     1072//    hGain.DrawCopy();
     1073
     1074
     1075    TF1 GainGaus( "GainGaus", "gaus", 0, hGain.GetXaxis()->GetXmax());
     1076    hGain.Rebin(2);
    10041077    hGain.Fit(&GainGaus, "N0QS");
     1078
     1079//    float gain_mean     = GainGaus.GetParameter(1);
     1080    float gain_median   = MedianOfH1(&hGain);
     1081    TH1F hNormGain          ("NormGain",      "median normalized Gain distribution",
     1082                            hGain.GetNbinsX(),
     1083                            (hGain.GetXaxis()->GetXmin())/gain_median,
     1084                            (hGain.GetXaxis()->GetXmax())/gain_median
     1085                            );
     1086    hNormGain.SetXTitle("gain [fraction of median gain value]");
     1087    hNormGain.SetYTitle("Counts");
     1088
     1089
     1090    hNormGain.Add(&hGain);
     1091//    hNormGain.SetStats(false);
     1092    hNormGain.GetXaxis()->SetRangeUser(
     1093                            1-(hNormGain.GetXaxis()->GetXmax()-1),
     1094                            hNormGain.GetXaxis()->GetXmax()
     1095                            );
     1096    hNormGain.DrawCopy();
     1097    hNormGain.Fit(&GainGaus, "N0QS");
     1098    GainGaus.DrawCopy("SAME");
     1099
     1100    //plot normailzed camera view
     1101    cam_canv.cd();
     1102    MHCamera hNormCam(fact);
     1103    hNormCam.SetStats(false);
     1104    for (UInt_t pixel =1; pixel <= 1440; pixel++)
     1105    {
     1106        hNormCam.SetBinContent(pixel, hcam.GetBinContent(pixel)/gain_median);
     1107    }
     1108
     1109    hNormCam.SetNameTitle( "gain","normalized [Median] Gain distribution over whole camera");
     1110    hNormCam.SetZTitle("Horst");
     1111    hNormCam.SetUsed(usePixel);
     1112    hNormCam.DrawCopy();
    10051113//    hGain.Scale(1/GainGaus.GetParameter(1));
    1006     hGain.DrawCopy();
    1007     GainGaus.DrawCopy("SAME");
     1114//    GainGaus.DrawCopy("SAME");
    10081115
    10091116    c2.cd(3);
    … …  
    10221129    gPad->SetLogy();
    10231130    hFitProb.DrawCopy();
     1131
     1132
     1133    TCanvas &c25 = d->AddTab("Spectra");
     1134    c25.Divide(2,1);
     1135    c25.cd(1);
     1136    gPad->SetLogy();
     1137    gPad->SetGrid();
     1138    hSum.DrawCopy("hist");
     1139    funcSum.DrawCopy("SAME");
     1140
     1141    c25.cd(2);
     1142    gPad->SetLogy();
     1143    gPad->SetGrid();
     1144    hSumScale.DrawCopy("hist");
     1145    funcScaled.DrawCopy("SAME");
     1146
    10241147
    10251148    /*
    … …  
    10291152    */
    10301153
    1031     d->SaveAs(outfile);
     1154//    d->SaveAs(outfile);
    10321155    return 0;
    10331156}
    … …  
    10731196    Double_t cross  = par[3];
    10741197    Double_t shift  = par[4];
    1075     Double_t powOfX = par[5];
     1198//    Double_t powOfX = par[5];
    10761199
    10771200    Double_t xTalk = 1;
    … …  
    10861209        y += xTalk*gauss;
    10871210
     1211//        xTalk = pow(cross, order*powOfX);
    10881212//        for (int j = 1; j < powOfX; j++)
    10891213//        {
    … …  
    11351259}
    11361260
     1261Double_t MedianOfH1 (
     1262        TH1*            inputHisto
     1263        )
     1264{
     1265   //compute the median for 1-d histogram h1
     1266   Int_t nbins = inputHisto->GetXaxis()->GetNbins();
     1267   Double_t *x = new Double_t[nbins];
     1268   Double_t *y = new Double_t[nbins];
     1269   for (Int_t i=0;i<nbins;i++) {
     1270      x[i] = inputHisto->GetXaxis()->GetBinCenter(i+1);
     1271      y[i] = inputHisto->GetBinContent(i+1);
     1272   }
     1273   Double_t median = TMath::Median(nbins,x,y);
     1274   delete [] x;
     1275   delete [] y;
     1276   return median;
     1277}
     1278// end of PlotMedianEachSliceOfPulse
     1279//----------------------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.