Changeset 14166 for fact


Ignore:
Timestamp:
06/13/12 12:46:17 (12 years ago)
Author:
Jens Buss
Message:
with thomas changes
File:
1 edited

Legend:

Unmodified
Added
Removed

  • fact/tools/marsmacros/singlepe.C

    r14165 r14166  
    33#include <TSystem.h>
    44#include <TF1.h>
     5#include <TProfile.h>
     6#include <TProfile2D.h>
    57#include <TMath.h>
    68#include <TFitResultPtr.h>
    … …  
    2830#include "MHCamera.h"
    2931#include "MGeomCamFACT.h"
     32#include "MRawRunHeader.h"
    3033
    3134using namespace std;
    … …  
    3639struct Single
    3740{
    38     float fSignal;
    39     float fTime;
     41    float    fSignal;
     42    float    fTime;
    4043};
    4144
    4245class MSingles : public MParContainer, public MCamEvent
    4346{
     47    Int_t fIntegrationWindow;
     48
    4449    vector<vector<Single>> fData;
    4550
    4651public:
    47     MSingles(const char *name=NULL, const char *title=NULL)
     52    MSingles(const char *name=NULL, const char *title=NULL) : fIntegrationWindow(0)
    4853    {
    4954        fName = name ? name : "MSingles";
    … …  
    6065    UInt_t GetNumPixels() const { return fData.size(); }
    6166
     67    void SetIntegrationWindow(Int_t w) { fIntegrationWindow = w; }
     68    Int_t GetIntegrationWindow() const { return fIntegrationWindow; }
     69
    6270    Bool_t GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const
    6371    {
     72        return kTRUE;
    6473    }
    6574    void   DrawPixelContent(Int_t num) const { }
    … …  
    6877};
    6978
     79class MH2F : public TH2F
     80{
     81public:
     82    Int_t Fill(Double_t x,Double_t y)
     83    {
     84        Int_t binx, biny, bin;
     85        fEntries++;
     86        binx = TMath::Nint(x+1);
     87        biny = fYaxis.FindBin(y);
     88        bin  = biny*(fXaxis.GetNbins()+2) + binx;
     89        AddBinContent(bin);
     90        if (fSumw2.fN) ++fSumw2.fArray[bin];
     91        if (biny == 0 || biny > fYaxis.GetNbins()) {
     92            if (!fgStatOverflows) return -1;
     93        }
     94        ++fTsumw;
     95        ++fTsumw2;
     96        fTsumwy  += y;
     97        fTsumwy2 += y*y;
     98        return bin;
     99    }
     100    ClassDef(MH2F, 1);
     101};
     102
     103class MProfile2D : public TProfile2D
     104{
     105public:
     106    Int_t Fill(Double_t x, Double_t y, Double_t z)
     107    {
     108        Int_t bin,binx,biny;
     109        fEntries++;
     110        binx =TMath::Nint(x+1);
     111        biny =fYaxis.FindBin(y);
     112        bin = GetBin(binx, biny);
     113        fArray[bin] += z;
     114        fSumw2.fArray[bin] += z*z;
     115        fBinEntries.fArray[bin] += 1;
     116        if (fBinSumw2.fN)  fBinSumw2.fArray[bin] += 1;
     117        if (biny == 0 || biny > fYaxis.GetNbins()) {
     118            if (!fgStatOverflows) return -1;
     119        }
     120        ++fTsumw;
     121        ++fTsumw2;
     122        fTsumwy  += y;
     123        fTsumwy2 += y*y;
     124        fTsumwz  += z;
     125        fTsumwz2 += z*z;
     126        return bin;
     127    }
     128    ClassDef(MProfile2D, 1);
     129};
     130
     131
    70132
    71133class MHSingles : public MH
    72134{
    73     MSingles *fSingles;
    74 
    75     TH2F fHist;
     135    MH2F       fSignal;
     136    MH2F       fTime;
     137    MProfile2D fPulse;
     138
     139    UInt_t fNumEntries;
     140
     141    MSingles               *fSingles;   //!
     142    MPedestalSubtractedEvt *fEvent;     //!
    76143
    77144public:
    78     MHSingles() : fSingles(0)
     145    MHSingles() : fNumEntries(0), fSingles(0), fEvent(0)
    79146    {
    80147        fName = "MHSingles";
    81148
    82         fHist.SetName("Name");
    83         fHist.SetTitle("Title");
    84         fHist.SetXTitle("X title");
    85         fHist.SetYTitle("Y title");
    86         fHist.SetDirectory(NULL);
     149        fSignal.SetName("Signal");
     150        fSignal.SetTitle("Title");
     151        fSignal.SetXTitle("X title");
     152        fSignal.SetYTitle("Y title");
     153        fSignal.SetDirectory(NULL);
     154
     155        fTime.SetName("Time");
     156        fTime.SetTitle("Title");
     157        fTime.SetXTitle("X title");
     158        fTime.SetYTitle("Y title");
     159        fTime.SetDirectory(NULL);
     160
     161        fPulse.SetName("Pulse");
     162        fPulse.SetTitle("Title");
     163        fPulse.SetXTitle("X title");
     164        fPulse.SetYTitle("Y title");
     165        fPulse.SetDirectory(NULL);
    87166    }
    88167
    … …  
    96175        }
    97176
     177        fEvent = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt");
     178        if (!fEvent)
     179        {
     180            *fLog << /* err << */ "MPedestalSubtractedEvt not found... abort." << endl;
     181            return kFALSE;
     182        }
     183
     184        fNumEntries = 0;
     185
    98186        return kTRUE;
    99187    }
    100     Bool_t ReInit(MParList *)
    101     {
     188    Bool_t ReInit(MParList *plist)
     189    {
     190        const MRawRunHeader *header = (MRawRunHeader*)plist->FindObject("MRawRunHeader");
     191        if (!header)
     192        {
     193            *fLog << /* err << */ "MRawRunHeader not found... abort." << endl;
     194            return kFALSE;
     195        }
     196
     197        const Int_t w = fSingles->GetIntegrationWindow();
     198
    102199        MBinning binsx, binsy;
    103200        binsx.SetEdges(fSingles->GetNumPixels(), -0.5, fSingles->GetNumPixels()-0.5);
    104         binsy.SetEdges(300, -25, 600-25);
    105 
    106         MH::SetBinning(fHist, binsx, binsy);
     201        binsy.SetEdges(300, -7.5*w, (60-7.5)*w);
     202
     203        MH::SetBinning(fSignal, binsx, binsy);
     204
     205        const UShort_t roi = header->GetNumSamples();
     206
     207        // Correct for DRS timing!!
     208        MBinning binst(roi, -0.5, roi-0.5);
     209        MH::SetBinning(fTime, binsx, binst);
     210
     211        MBinning binspy(2*roi, -roi-0.5, roi-0.5);
     212        MH::SetBinning(fPulse, binsx, binspy);
    107213
    108214        return kTRUE;
    … …  
    111217    Int_t Fill(const MParContainer *par, const Stat_t weight=1)
    112218    {
     219        const Float_t *ptr = fEvent->GetSamples(0);
     220        const Short_t  roi = fEvent->GetNumSamples();
     221
    113222        for (unsigned int i=0; i<fSingles->GetNumPixels(); i++)
    114223        {
    … …  
    116225
    117226            for (unsigned int j=0; j<result.size(); j++)
    118                 fHist.Fill(i, result[j].fSignal);
    119         }
     227            {
     228                fSignal.Fill(i, result[j].fSignal);
     229                fTime.Fill  (i, result[j].fTime);
     230
     231                if (!ptr)
     232                    continue;
     233
     234                const Short_t tm = floor(result[j].fTime);
     235
     236                for (int s=0; s<roi; s++)
     237                    fPulse.Fill(i, s-tm, ptr[s]);
     238            }
     239
     240            ptr += roi;
     241        }
     242
     243        fNumEntries++;
    120244
    121245        return kTRUE;
    122246    }
    123247
    124     TH2 *GetHist() { return &fHist; }
     248    TH2 *GetSignal() { return &fSignal; }
     249    TH2 *GetTime()   { return &fTime; }
     250    TH2 *GetPulse()  { return &fPulse; }
     251
     252    UInt_t GetNumEntries() const { return fNumEntries; }
    125253
    126254    void Draw(Option_t *opt)
    … …  
    130258        AppendPad("");
    131259
    132         pad->Divide(1,1);
     260        pad->Divide(2,2);
    133261
    134262        pad->cd(1);
    135 
    136263        gPad->SetBorderMode(0);
    137264        gPad->SetFrameBorderMode(0);
    138 
    139         fHist.Draw("colz");
     265        fSignal.Draw("colz");
     266
     267        pad->cd(2);
     268        gPad->SetBorderMode(0);
     269        gPad->SetFrameBorderMode(0);
     270        fTime.Draw("colz");
     271
     272        pad->cd(3);
     273        gPad->SetBorderMode(0);
     274        gPad->SetFrameBorderMode(0);
     275        fPulse.Draw("colz");
    140276    }
    141277
    … …  
    143279};
    144280
     281MStatusDisplay *d = new MStatusDisplay;
     282
    145283MSingles *fSingles;
     284
     285TH1F *fluct1 = new TH1F("", "", 100, -50, 50);
     286TH1F *fluct2 = new TH1F("", "", 100, -50, 50);
     287
     288TGraph weights("weights.txt");
    146289
    147290Int_t PreProcess(MParList *plist)
    … …  
    158301    }
    159302
     303    d->AddTab("Fluct");
     304    fluct2->Draw();
     305    fluct1->Draw("same");
     306
     307    fSingles->SetIntegrationWindow(20);
     308
     309    //if (weights.GetN()==0)
     310    //    return kFALSE;
     311
    160312    return kTRUE;
    161313}
    … …  
    163315Int_t Process()
    164316{
    165 
    166     for (int i=0; i<fEvent->GetNumPixels(); i++)
    167     {
    168         vector<Single> &result = fSingles->GetVector(i);
     317    const UInt_t roi = fEvent->GetNumSamples();
     318
     319    const size_t navg             = 10;
     320    const float  threshold        = 5;
     321    const UInt_t start_skip       = 20;
     322    const UInt_t end_skip         = 10;
     323    const UInt_t integration_size = 10*3;
     324
     325    vector<float> val(roi-navg);//result of first sliding average
     326    for (int pix=0; pix<fEvent->GetNumPixels(); pix++)
     327    {
     328        vector<Single> &result = fSingles->GetVector(pix);
    169329        result.clear();
    170330
    171         const UInt_t roi = fEvent->GetNumSamples();
    172 
    173         const Float_t *ptr = fEvent->GetSamples(i);
    174 
    175     //---------------------------
    176         // ...Extract pulses here....
    177     //---------------------------
    178 
    179         //sliding average parameter (center weight = 3)
    180     int avg1 = 8;
    181         vector<float> Vslide(roi-2*avg1);//result of first sliding average
     331        const Float_t *ptr = fEvent->GetSamples(pix);
    182332
    183333        //Sliding window
    184         for (UInt_t k=avg1; k<roi-avg1; k++)
    185         {
    186             Vslide[k-avg1]=ptr[k];
    187             for(unsigned int offset=0; offset<avg1; offset++)
    188             {
    189                 Vslide[k-avg1] += ptr[k+offset];
    190                 Vslide[k-avg1] += ptr[k-offset];
    191             }
    192             Vslide[k-avg1] /= float(2*avg1+1);
    193         }
    194 
    195     //peak finding via threshold
    196         float threshold = 5;
    197         Single single;
    198         UInt_t start_sample = 5;
    199         UInt_t end_sample = roi-2*avg1-21;
    200         UInt_t min_dist = 20;
    201         UInt_t integration_size = 5;
    202         UInt_t integration_delay = 8;
    203 
    204         for (UInt_t k=start_sample; k<end_sample; k++)
    205         {
    206       //search for threshold crossings
    207             if(     (Vslide[k-start_sample]<threshold)
    208                     &&(Vslide[k-1]<threshold)
    209                     &&(Vslide[k]>threshold)
    210                     &&(Vslide[k+5]>threshold)
    211                     )
    212             {
    213                 //average baseline before crossing
    214                 //double baseline = (Vslide[k-start_sample]+Vslide[k-start_sample+1]+Vslide[k-start_sample-1])/3.;
    215                 single.fSignal = 0;
    216                 single.fTime = k+integration_delay;
    217                 for(UInt_t l=integration_delay; l<integration_delay+integration_size; l++)
    218                 {
    219                     single.fSignal+=Vslide[k+l];
    220                 }
    221                 //single.fSignal-=baseline;//baseline subtraction
    222                 result.push_back(single);
    223                 k+=min_dist;
    224             }
     334        for (size_t i=0; i<roi-navg; i++)
     335        {
     336            val[i] = 0;
     337            for (size_t j=i; j<i+navg; j++)
     338                val[i] += ptr[j];
     339            val[i] /= navg;
     340        }
     341
     342        //peak finding via threshold
     343        for (UInt_t i=start_skip; i<roi-navg-end_skip-10; i++)
     344        {
     345            //search for threshold crossings
     346            if (val[i+0]>threshold ||
     347                val[i+4]>threshold ||
     348
     349                val[i+5]<threshold ||
     350                val[i+9]<threshold)
     351                continue;
     352
     353            // Evaluate arrival time more precisely!!!
     354            // Get a better integration window
     355
     356            Single single;
     357            single.fSignal = 0;
     358            single.fTime   = i;
     359
     360            // Crossing of the threshold is at 5
     361            for (UInt_t j=i+5; j<i+5+integration_size; j++)
     362                single.fSignal += ptr[j];
     363
     364            result.push_back(single);
     365
     366            i += 5+integration_size;
    225367        }
    226368    }
    … …  
    233375}
    234376
    235 Double_t
    236 spek_fit(
    237         Double_t*   x,
    238         Double_t*   par
    239         );
    240 
    241 int singlepe_final()
     377Double_t fcn(Double_t *x, Double_t *par);
     378
     379int singlepe_final2()
    242380{
    243381    // ======================================================
    244382
    245     const char *drsfile = "/fact/raw/2012/02/25/20120225_004.drs.fits.gz";
     383    const char *drsfile = "/fact/raw/2012/06/01/20120601_003.drs.fits.gz";
    246384
    247385    MDirIter iter;
    … …  
    258396    // map file to use (get that from La Palma!)
    259397    const char *map = usemap ? "/scratch_nfs/FACTmap111030.txt" : NULL;
    260     gStyle->SetPalette(1,0);
     398
    261399    // ------------------------------------------------------
    262400
    … …  
    295433    // ======================================================
    296434
    297     MStatusDisplay *d = new MStatusDisplay;
     435    //MStatusDisplay *d = new MStatusDisplay;
    298436
    299437    MBadPixelsCam badpixels;
    … …  
    316454    //MDrsCalibrationTime timecam;
    317455
    318     gStyle->SetOptFit(kTRUE);
     456    MH::SetPalette();
    319457
    320458    // ======================================================
    … …  
    326464
    327465    MSingles singles;
     466    MRawRunHeader header;
    328467
    329468    MParList plist1;
    … …  
    332471    plist1.AddToList(&badpixels);
    333472    plist1.AddToList(&singles);
     473    plist1.AddToList(&header);
    334474
    335475    MEvtLoop loop1("DetermineCalConst");
    … …  
    371511        return 11;
    372512
     513    gStyle->SetOptFit(kTRUE);
     514
    373515    MHSingles* h = (MHSingles*) plist1.FindObject("MHSingles");
    374516    if (!h)
    375517        return 0;
    376518
    377     TH2 *hist = h->GetHist();
     519    TH2 *hsignal = h->GetSignal();
     520    TH2 *htime   = h->GetTime();
     521    TH2 *hpulse  = h->GetPulse();
     522
     523    d->AddTab("Time");
     524    TH1 *ht = htime->ProjectionY();
     525    ht->Scale(1./1440);
     526    ht->Draw();
     527
     528    d->AddTab("Pulse");
     529    TH1 *hp = hpulse->ProjectionY();
     530    hp->Scale(1./1440);
     531    hp->Draw();
     532
    378533
    379534    // ------------------ Spectrum Fit ---------------
    … …  
    381536    // access the spektrumhistogram by the pointer *hist
    382537
    383     hist->SetTitle("Amplitude Spectrum vs. Pixel ID");
    384     hist->SetXTitle("Pixel SoftID");
    385     hist->SetYTitle("Amplitude [mV]");
    386     hist->SetZTitle("counts");
    387 
    388538    MGeomCamFACT fact;
    389539    MHCamera hcam(fact);
    390540
    391     struct fit
    392     {
    393         float fAmplitude;
    394         float fGain;
    395         float f2GainRMS;
    396         float fCrosstlk;
    397         float fOffset;
    398     };
    399     fit fitParameters[singles.GetNumPixels()];
    400 
    401     TCanvas &canv1 = d->AddTab("Fits");
    402     TCanvas &canv2 = d->AddTab("Distributions");
    403     canv2.Divide(3,2);
    404     TCanvas &gain_canv = d->AddTab("Gain Distribution");
    405 
    406541    //built an array of Amplitude spectra
    407     TH1D*       hAmplSpek[singles.GetNumPixels()];
    408     TString     histotitle;
    409 
    410     TF1 gaus4("g4","gaus");
    411     gaus4.SetLineColor(2);
    412 
    413     TH1F        hAmplitude(
    414                     "hAmplitude",
    415                     "Amplitude of Single Preak",
    416                     200,
    417                     50,
    418                     250
    419                     );
    420     hAmplitude.GetXaxis()->SetTitle("Amplitude (from Integral 5 Slices) [mV]");
    421     hAmplitude.GetYaxis()->SetTitle("counts");
    422 
    423     TH1F        hGain(
    424                     "hGain",
    425                     "Gain",
    426                     40*5,
    427                     30,
    428                     70
    429                     );
    430     hGain.GetXaxis()->SetTitle("Gain (from Integral 5 Slices) [mV]");
    431     hGain.GetYaxis()->SetTitle("counts");
    432 
    433     TH1F        hGainRMS(
    434                     "hSigma",
    435                     "Sigma of gauss peak in gainfit",
    436                     20*5,
    437                     0,
    438                     40
    439                     );
    440     hGainRMS.GetXaxis()->SetTitle("Sigma (from Integral 5 Slices) [mV]");
    441     hGainRMS.GetYaxis()->SetTitle("counts");
    442 
    443     TH1F        hCrosstlk(
    444                     "hCrosstlk",
    445                     "Parameter proportional to Crosstalk",
    446                     30*5,
    447                     0,
    448                     4
    449                     );
    450     hCrosstlk.GetXaxis()->SetTitle("~Crosstalk (from Integral 5 Slices) [mV]");
    451     hCrosstlk.GetYaxis()->SetTitle("counts");
    452 
    453     TH1F        hOffset(
    454                     "hOffset",
    455                     "X-Offset of Spectrum",
    456                     100*2,
    457                     0,
    458                     4
    459                     );
    460     hOffset.GetXaxis()->SetTitle(" X-Offset of Spectrum (from Integral 5 Slices) [mV]");
    461     hOffset.GetYaxis()->SetTitle("counts");
    462 
    463     TH1F        hFitProb(
    464                     "hFitProb",
    465                     "Fit Probability",
    466                     100*2,
    467                     0,
    468                     100
    469                     );
    470     hFitProb.GetXaxis()->SetTitle("fit-probability [%]");
    471     hFitProb.GetYaxis()->SetTitle("counts");
     542    TH1F hAmplitude("Rate",      "Average number of dark counts per event",
     543                    100,  0,  10);
     544    TH1F hGain     ("Gain",      "Gain distribution",
     545                     50,  100,   250);
     546    TH1F hGainRMS  ("RelSigma",  "Distribution of rel. Sigma",
     547                    100,   0,   30);
     548    TH1F hCrosstlk ("Crosstalk", "Distribution of crosstalk probability",
     549                    100,   0,    25);
     550    TH1F hOffset   ("Offset",    "Baseline offset distribution",
     551                    50, -7.5,  2.5);
     552    TH1F hFitProb  ("FitProb",   "FitProb distribution",
     553                     50,   0,  100);
     554    TH1F hSum1     ("Sum1",      "Sum",
     555                    100,    0,  2000);
     556    TH1F hSum2     ("Sum2",      "Sum",
     557                    100,    0,   10);
    472558
    473559    // define fit function for Amplitudespectrum
    474     TF1         fspek_fit("fspek_fit", spek_fit, 10, 600, 5);
    475     fspek_fit.SetParNames("Ampl 1.Max", "Gain", "Sigma", "Ampl 2.Max", "Offset");
    476 
    477     // fit parameters
    478     Double_t    par[5];
     560    TF1 func("spektrum", fcn, 0, 1000, 5);
     561    func.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset");
     562    func.SetLineColor(kBlue);
    479563
    480564    // Map for which pixel shall be plotted and which not
    481565    TArrayC usePixel(1440);
     566    memset(usePixel.GetArray(), 1, 1440);
    482567
    483568    // List of Pixel that should be ignored in camera view
    484     Int_t doNotUse[12]= {
    485         424,
    486         923,
    487         1208,
    488         583,
    489         830,
    490         1399,
    491         113,
    492         115,
    493         354,
    494         423,
    495         1195,
    496         1393
    497     };
     569    usePixel[424]  = 0;
     570    usePixel[923]  = 0;
     571    usePixel[1208] = 0;
     572    usePixel[583]  = 0;
     573    usePixel[830]  = 0;
     574    usePixel[1399] = 0;
     575    usePixel[113]  = 0;
     576    usePixel[115]  = 0;
     577    usePixel[354]  = 0;
     578    usePixel[423]  = 0;
     579    usePixel[1195] = 0;
     580    usePixel[1393] = 0;
     581
     582    int count_ok = 0;
    498583
    499584    // Begin of Loop over Pixels
    500585    for (UInt_t pixel = 0; pixel < singles.GetNumPixels(); pixel++)
    501586    {
    502         usePixel[pixel] = 0;
    503 
    504         //define projection histogram title
    505         histotitle = "hPixelPro_";
    506         histotitle += pixel;
     587        if (usePixel[pixel]==0)
     588            continue;
    507589
    508590        //Projectipon of a certain Pixel out of Ampl.Spectrum
    509         hAmplSpek[pixel] = hist->ProjectionY( histotitle , pixel+1, pixel+1);
    510 
    511         //set histogram and Axis title
    512         histotitle = "Amplitude Spectrum of Pixel #";
    513         histotitle += pixel;
    514         hAmplSpek[pixel]->SetTitle(histotitle);
    515 
    516         // define appearance of fit
    517         fspek_fit.SetLineColor(kBlue);
    518         fspek_fit.SetLineStyle(2);
    519 
    520         par[0] = hAmplSpek[pixel]->GetBinContent(hAmplSpek[pixel]->GetMaximumBin());
     591        TH1D *hist = hsignal->ProjectionY("proj", pixel+1, pixel+1);
     592        hist->SetDirectory(0);
     593
     594        const Int_t    maxbin   = hist->GetMaximumBin();
     595        const Double_t binwidth = hist->GetBinWidth(maxbin);
     596        const Double_t ampl     = hist->GetBinContent(maxbin);
     597
     598        double fwhm = 0;
     599        for (int i=1; i<maxbin; i++)
     600            if (hist->GetBinContent(maxbin-i)+hist->GetBinContent(maxbin+i)<ampl*2)
     601            {
     602                fwhm = (2*i+1)*binwidth;
     603                break;
     604            }
     605
     606        if (fwhm==0)
     607        {
     608            cout << "Could not determine start value for sigma." << endl;
     609            continue;
     610        }
     611
     612        // Amplitude, Gain, Sigma, Crosstalk Prob, Shift
    521613        // par[1] + par[4] is the position of the first maximum
    522         par[1] = 0.8*hAmplSpek[pixel]->GetBinCenter(hAmplSpek[pixel]->GetMaximumBin());
    523         par[2] = 0.2*par[1];
    524         par[3] = 0.1*par[0];
    525         par[4] = 0.2*hAmplSpek[pixel]->GetBinCenter(hAmplSpek[pixel]->GetMaximumBin());
    526         fspek_fit.SetParameters(par);
     614        Double_t par[5] =
     615        {
     616            ampl, hist->GetBinCenter(maxbin), fwhm*1.4, 0.1, -10
     617        };
     618
     619        func.SetParameters(par);
     620
     621        const double min = par[1]-fwhm*1.4;
     622        const double max = par[1]*5;
    527623
    528624        //Fit Pixels spectrum
    529         canv1.cd();
    530         TFitResultPtr fitStatus(hAmplSpek[ pixel ]->Fit(&fspek_fit,"+S","",0.8*par[1],5*par[1]));
    531         double fit_prob = fitStatus->Prob();
    532 
    533         //cout << "pixel: " << pixel << "\t status: " << fitStatus << endl;
    534 
    535         fitParameters[pixel].fAmplitude = fspek_fit.GetParameter(0);
    536         fitParameters[pixel].fGain      = fspek_fit.GetParameter(1);
    537         fitParameters[pixel].f2GainRMS  = fspek_fit.GetParameter(2);
    538         fitParameters[pixel].fCrosstlk  = fspek_fit.GetParameter(3)/fspek_fit.GetParameter(0);
    539         fitParameters[pixel].fOffset    = fspek_fit.GetParameter(4);
    540 
    541         gain_canv.cd();
    542         gPad->SetLogy(1);
    543         hGain.Fill(fitParameters[pixel].fGain);
    544         hGain.Fit(&gaus4);
    545         hGain.DrawCopy();
    546 
    547         canv2.cd(2);
    548         gPad->SetLogy(1);
    549         hGain.DrawCopy();
    550 
    551         canv2.cd(3);
    552         gPad->SetLogy(1);
    553         hGainRMS.Fill(fitParameters[pixel].f2GainRMS);
    554         hGainRMS.DrawCopy();
    555 
    556         canv2.cd(1);
    557         gPad->SetLogy(1);
    558         hAmplitude.Fill(fitParameters[pixel].fAmplitude);
    559         hAmplitude.DrawCopy();
    560 
    561         canv2.cd(4);
    562         gPad->SetLogy(1);
    563         hCrosstlk.Fill(fitParameters[pixel].fCrosstlk);
    564         hCrosstlk.DrawCopy();
    565 
    566         canv2.cd(5);
    567         gPad->SetLogy(1);
    568         hOffset.Fill(fitParameters[pixel].fOffset);
    569         hOffset.DrawCopy();
    570 
    571         canv2.cd(6);
    572         gPad->SetLogy(1);
     625        const TFitResultPtr rc = hist->Fit(&func, "N0QS", "", min, max);
     626
     627        const bool ok = int(rc)==0;
     628
     629        const double fit_prob   = rc->Prob();
     630        const float  fGain      = func.GetParameter(1);
     631        const float  fAmplitude = func.Integral(0, 7*fGain)/h->GetNumEntries(); //func.GetParameter(0);
     632        const float  f2GainRMS  = func.GetParameter(2);
     633        const float  fCrosstlk  = func.GetParameter(3);
     634        const float  fOffset    = func.GetParameter(4)/singles.GetIntegrationWindow();
     635
     636        hAmplitude.Fill(fAmplitude);
     637        hGain.Fill(fGain);
     638        //hGainRMS.Fill(f2GainRMS);
     639        hCrosstlk.Fill(fCrosstlk*100);
     640        hOffset.Fill(fOffset);
    573641        hFitProb.Fill(100*fit_prob);
    574         hFitProb.DrawCopy();
    575 
    576 
    577         hcam.SetBinContent(pixel+1, fitParameters[pixel].fGain);
    578         if (!fitStatus)
    579         {
    580             usePixel[pixel] = 1;
    581             cout << "Pixel: " << pixel << " success!" << endl << endl;
    582         }
    583         else
    584         {
    585             cout << "Pixel: " << pixel << " failed!" << endl << endl;
    586         }
    587 
    588     }
    589 //    TCanvas &canvBroken[12];
    590 //    TString name;
    591     for  (int i = 0; i < 12; i++)
    592     {
    593 //        name    = "Pixel";
    594 //        name   += doNotUse[i];
    595 //        canvBroken[i] = d->AddTab(name);
    596 //        hAmplSpek[doNotUse[ i ]]->Draw();
    597 
    598         usePixel[doNotUse[ i ]] = 0;
     642        hGainRMS.Fill(100*f2GainRMS/fGain);
     643
     644        hcam.SetBinContent(pixel+1, fGain);
     645
     646        usePixel[pixel] = ok;
     647
     648        if (!ok)
     649        {
     650            cout << "Pixel " << setw(4) << pixel << ": failed!" << endl;
     651            continue;
     652        }
     653
     654        for (int b=1; b<=hist->GetNbinsX(); b++)
     655        {
     656            hSum1.Fill(hist->GetBinCenter(b), hist->GetBinContent(b));
     657            hSum2.Fill((hist->GetBinCenter(b)-fOffset)/fGain, hist->GetBinContent(b));
     658        }
     659
     660        if (count_ok==0)
     661        {
     662            TCanvas &c = d->AddTab(Form("Pix%d", pixel));
     663            c.cd();
     664            hist->SetName(Form("Pix%d", pixel));
     665            hist->DrawCopy()->SetDirectory(0);
     666            //func.DrawCopy("same");
     667        }
     668
     669        count_ok++;
     670
     671        delete hist;
    599672    }
    600673
    … …  
    606679    hcam.DrawCopy();
    607680
    608 
    609     canv1.cd();
     681    TCanvas &c11 = d->AddTab("SumHist");
     682    gPad->SetLogy();
     683    hSum1.DrawCopy();
     684
     685    TCanvas &c12 = d->AddTab("GainHist");
     686    gPad->SetLogy();
     687    hSum2.DrawCopy();
     688
     689    TCanvas &c2 = d->AddTab("Result");
     690    c2.Divide(3,2);
     691    c2.cd(1);
     692    gPad->SetLogy();
     693    hAmplitude.DrawCopy();
     694
     695    c2.cd(2);
     696    gPad->SetLogy();
     697    hGain.DrawCopy();
     698
     699    c2.cd(3);
     700    gPad->SetLogy();
     701    hGainRMS.DrawCopy();
     702
     703    c2.cd(4);
     704    gPad->SetLogy();
     705    hCrosstlk.DrawCopy();
     706
     707    c2.cd(5);
     708    gPad->SetLogy();
     709    hOffset.DrawCopy();
     710
     711    c2.cd(6);
     712    gPad->SetLogy();
     713    hFitProb.DrawCopy();
    610714
    611715    /*
    612     TString title = "--  Calibrated signal #";
    613     title += seq.GetSequence();
    614     title += " (";
    615     title += drsfile;
    616     title += ")  --";
    617     d->SetTitle(title, kFALSE);
    618 
    619     TString path;
    620     path += Form("%s/20%6d_%03d-calibration.root", outpath,
    621                  seq.GetSequence()/1000, seq.GetSequence()%1000);
    622 
    623     d->SaveAs(path);
    624     */
    625 
     716    TCanvas &c3 = d->AddTab("Separation");
     717    gPad->SetLogy();
     718    hSep.DrawCopy();
     719     */
    626720    return 0;
    627721}
    628722
    629 Double_t
    630 spek_fit(
    631         Double_t*   x,
    632         Double_t*   par
    633         )
    634 {
    635     Double_t arg        = 0;
    636     Double_t cross      = 1;
    637     Double_t peak       = 0;
    638 //    Double_t arg2     = 0;
    639 //    Double_t arg3     = 0;
    640 //    Double_t arg4     = 0;
    641 
    642 //  if (par[0] != 0) cross   = par[3]/par[0];
    643 //  if (par[2] != 0) arg     = (x[0] - par[1] + par[4])/par[2];
    644 //  if (par[2] != 0) arg2    = (x[0] - 2*par[1] + par[4])/par[2];
    645 //  if (par[2] != 0) arg3    = (x[0] - 3*par[1] + par[4])/par[2];
    646 //  if (par[2] != 0) arg4    = (x[0] - 4*par[1] + par[4])/par[2];
    647 //    if (par[0] = 0) return 0;
    648 //    if (par[2] = 0) return 0;
    649 
    650     if (par[0] != 0) cross   = par[3]/par[0];
     723Double_t fcn(Double_t *xx, Double_t *par)
     724{
     725    Double_t x     = xx[0];
     726
     727    Double_t ampl  = par[0];
     728    Double_t gain  = par[1];
     729    Double_t sigma = par[2];
     730    Double_t cross = par[3];
     731    Double_t shift = par[4];
     732
     733    Double_t xTalk = 1;
     734
     735    Double_t y = 0;
    651736    for (int order = 1; order < 5; order++)
    652737    {
    653         arg = (x[0] - order*par[1] - par[4] )/par[2];
    654         Double_t xTalk    = TMath::Power(cross, order-1);
    655 //        Double_t gauss    = TMath::Exp(-0.5 * arg * arg/order);
    656         Double_t gauss    = TMath::Exp(-0.5 * arg * arg);
    657         peak += xTalk*par[0]*gauss;
    658     }
    659 
    660 //  peak += par[0]*TMath::Exp(-0.5 * arg * arg);
    661 //  peak += cross*par[0]*TMath::Exp(-0.5 * arg2 * arg2 /2);
    662 //  peak += cross*cross*par[0]*TMath::Exp(-0.5 * arg3 * arg3 /3);
    663 //  peak += cross*cross*cross*par[0]*TMath::Exp(-0.5 * arg4 * arg4 /4);
    664 
    665     return peak;
     738        Double_t arg   = (x - order*gain - shift)/sigma;
     739
     740        Double_t gauss = TMath::Exp(-0.5 * arg * arg/order);
     741
     742        y += xTalk*gauss;
     743
     744        xTalk *= cross;
     745    }
     746
     747    return ampl*y;
    666748}
Note: See TracChangeset for help on using the changeset viewer.