Index: /fact/tools/marsmacros/recalcSinglepe.C =================================================================== --- /fact/tools/marsmacros/recalcSinglepe.C (revision 14277) +++ /fact/tools/marsmacros/recalcSinglepe.C (revision 14277) @@ -0,0 +1,1055 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include // TGHProgressBar +#include +#include + +#include +#include +#include + +#include "MH.h" +#include "MLog.h" +#include "MLogManip.h" +#include "MDirIter.h" +#include "MFillH.h" +#include "MEvtLoop.h" +#include "MCamEvent.h" +#include "MGeomApply.h" +#include "MTaskList.h" +#include "MParList.h" +#include "MContinue.h" +#include "MBinning.h" +#include "MHDrsCalib.h" +#include "MStatusArray.h" +#include "MDrsCalibApply.h" +#include "MRawFitsRead.h" +#include "MBadPixelsCam.h" +#include "MStatusDisplay.h" +#include "MTaskInteractive.h" +#include "MPedestalSubtractedEvt.h" +#include "MHCamera.h" +#include "MGeomCamFACT.h" +#include "MRawRunHeader.h" + +using namespace std; +using namespace TMath; + +MPedestalSubtractedEvt *fEvent = 0; +MLog *fLog = &gLog; + +struct Single +{ + float fSignal; + float fTime; +}; + +class MSingles : public MParContainer, public MCamEvent +{ + Int_t fIntegrationWindow; + + vector> fData; + +public: + MSingles(const char *name=NULL, const char *title=NULL) : fIntegrationWindow(0) + { + fName = name ? name : "MSingles"; + } + + void InitSize(const UInt_t i) + { + fData.resize(i); + } + + vector &operator[](UInt_t i) { return fData[i]; } + vector &GetVector(UInt_t i) { return fData[i]; } + + UInt_t GetNumPixels() const { return fData.size(); } + + void SetIntegrationWindow(Int_t w) { fIntegrationWindow = w; } + Int_t GetIntegrationWindow() const { return fIntegrationWindow; } + + Bool_t GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const + { + return kTRUE; + } + void DrawPixelContent(Int_t num) const { } + + ClassDef(MSingles, 1) +}; +/* +class MH2F : public TH2F +{ +public: + Int_t Fill(Double_t x,Double_t y) + { + Int_t binx, biny, bin; + fEntries++; + binx = TMath::Nint(x+1); + biny = fYaxis.FindBin(y); + bin = biny*(fXaxis.GetNbins()+2) + binx; + AddBinContent(bin); + if (fSumw2.fN) ++fSumw2.fArray[bin]; + if (biny == 0 || biny > fYaxis.GetNbins()) { + if (!fgStatOverflows) return -1; + } + ++fTsumw; + ++fTsumw2; + fTsumwy += y; + fTsumwy2 += y*y; + return bin; + } + ClassDef(MH2F, 1); +}; + +class MProfile2D : public TProfile2D +{ +public: + Int_t Fill(Double_t x, Double_t y, Double_t z) + { + Int_t bin,binx,biny; + fEntries++; + binx =TMath::Nint(x+1); + biny =fYaxis.FindBin(y); + bin = GetBin(binx, biny); + fArray[bin] += z; + fSumw2.fArray[bin] += z*z; + fBinEntries.fArray[bin] += 1; + if (fBinSumw2.fN) fBinSumw2.fArray[bin] += 1; + if (biny == 0 || biny > fYaxis.GetNbins()) { + if (!fgStatOverflows) return -1; + } + ++fTsumw; + ++fTsumw2; + fTsumwy += y; + fTsumwy2 += y*y; + fTsumwz += z; + fTsumwz2 += z*z; + return bin; + } + ClassDef(MProfile2D, 1); +}; +*/ + + +class MHSingles : public MH +{ + TH2F fSignal; + TH2F fTime; + TProfile2D fPulse; + + UInt_t fNumEntries; + + MSingles *fSingles; //! + MPedestalSubtractedEvt *fEvent; //! + +public: + MHSingles() : fNumEntries(0), fSingles(0), fEvent(0) + { + fName = "MHSingles"; + + fSignal.SetName("Signal"); + fSignal.SetTitle("pulse integral spectrum"); + fSignal.SetXTitle("pixel [SoftID]"); + fSignal.SetYTitle("time [a.u.]"); + fSignal.SetDirectory(NULL); + + fTime.SetName("Time"); + fTime.SetTitle("arival time spectrum"); + fTime.SetXTitle("pixel [SoftID]"); + fTime.SetYTitle("time [a.u.]"); + fTime.SetDirectory(NULL); + + fPulse.SetName("Pulse"); + fPulse.SetTitle("average pulse shape spectrum"); + fPulse.SetXTitle("pixel [SoftID]"); + fPulse.SetYTitle("time [a.u.]"); + fPulse.SetDirectory(NULL); + } + + Bool_t SetupFill(const MParList *plist) + { + fSingles = (MSingles*)plist->FindObject("MSingles"); + if (!fSingles) + { + *fLog << /* err << */ "MSingles not found... abort." << endl; + return kFALSE; + } + + fEvent = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt"); + if (!fEvent) + { + *fLog << /* err << */ "MPedestalSubtractedEvt not found... abort." << endl; + return kFALSE; + } + + fNumEntries = 0; + + return kTRUE; + } + Bool_t ReInit(MParList *plist) + { + const MRawRunHeader *header = (MRawRunHeader*)plist->FindObject("MRawRunHeader"); + if (!header) + { + *fLog << /* err << */ "MRawRunHeader not found... abort." << endl; + return kFALSE; + } + + const Int_t w = fSingles->GetIntegrationWindow(); + + MBinning binsx, binsy; + binsx.SetEdges(fSingles->GetNumPixels(), -0.5, fSingles->GetNumPixels()-0.5); + binsy.SetEdges(2*150, -2*7.5*w, 2*(60-7.5)*w); + + MH::SetBinning(fSignal, binsx, binsy); + + const UShort_t roi = header->GetNumSamples(); + + // Correct for DRS timing!! + MBinning binst(roi, -0.5, roi-0.5); + MH::SetBinning(fTime, binsx, binst); + + MBinning binspy(2*roi, -roi-0.5, roi-0.5); + MH::SetBinning(fPulse, binsx, binspy); + + return kTRUE; + } + + Int_t Fill(const MParContainer *par, const Stat_t weight=1) + { + const Float_t *ptr = fEvent->GetSamples(0); + const Short_t roi = fEvent->GetNumSamples(); + + for (unsigned int i=0; iGetNumPixels(); i++) + { + const vector &result = fSingles->GetVector(i); + + for (unsigned int j=0; jDivide(2,2); + + pad->cd(1); + gPad->SetBorderMode(0); + gPad->SetFrameBorderMode(0); + fSignal.Draw("colz"); + + pad->cd(2); + gPad->SetBorderMode(0); + gPad->SetFrameBorderMode(0); + fTime.Draw("colz"); + + pad->cd(3); + gPad->SetBorderMode(0); + gPad->SetFrameBorderMode(0); + fPulse.Draw("colz"); + } + + void DrawCopy(Option_t *opt) + { + TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); + + AppendPad(""); + + pad->Divide(2,2); + + pad->cd(1); + gPad->SetBorderMode(0); + gPad->SetFrameBorderMode(0); + fSignal.DrawCopy("colz"); + + pad->cd(2); + gPad->SetBorderMode(0); + gPad->SetFrameBorderMode(0); + fTime.DrawCopy("colz"); + + pad->cd(3); + gPad->SetBorderMode(0); + gPad->SetFrameBorderMode(0); + fPulse.DrawCopy("colz"); + } + + ClassDef(MHSingles, 1) +}; + +MStatusDisplay *d = new MStatusDisplay; + +MSingles *fSingles; + +Int_t PreProcess(MParList *plist) +{ +// fSingles = (MSingles*)plist->FindCreateObj("MSingles"); +// if (!fSingles) +// return kFALSE; + +// fEvent = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt"); +// if (!fEvent) +// { +// *fLog << /* err << */ "MPedestalSubtractedEvt not found... abort." << endl; +// return kFALSE; +// } + +// d->AddTab("Fluct"); +// fluct2->Draw(); +// fluct1->Draw("same"); + +// fSingles->SetIntegrationWindow(20); + +// //if (weights.GetN()==0) +// // return kFALSE; + + return kTRUE; +} + +Int_t Process() +{ + return kTRUE; +} + +Int_t PostProcess() +{ + return kTRUE; +} + +Double_t fcn(Double_t *x, Double_t *par); + +int recalcSinglepe(const char *file, const char *outfile="", TString option="") +{ + bool save = 0; + + // ====================================================== + //save options + + if ( option.Contains("S") ){ + save = 1; + cout << "...will save" << endl; + } + + //Maximum pe order to which the sepctrum will be fitted + Int_t maxOrder = 5; + + // ====================================================== + + cout << file << endl; + TFile fi(file, "READ"); + if (!fi.IsOpen()) + { + cout << "ERROR - Could not find file " << file << endl; + return 2; + } + + MStatusArray arr; + if (arr.Read()<=0) + { + cout << "ERROR - could not read MStatusArray." << endl; + return 2; + } + + // ====================================================== + + // Load histogramm for fluctuations of all Pixel + TH1 *fluct = (TH1*)arr.FindObjectInCanvas("", "TH1F", "Fluct"); + if (!fluct) + { + cout << "WARNING - Reading of fluct failed." << endl; + return 2; + } + + d->AddTab("Fluct"); + fluct->DrawCopy(); + + + MHSingles* singles = (MHSingles*) arr.FindObjectInCanvas("MHSingles", "MHSingles", "MHSingles"); + if (!singles) + { + cout << "WARNING - Reading of singles failed." << endl; + return 2; + } + + d->AddTab("MHSingles"); + singles->DrawCopy("colz"); + + TH1 *time = (TH1*)arr.FindObjectInCanvas("Time_py", "TH1D", "Time"); + if (!fluct) + { + cout << "WARNING - Reading of time failed." << endl; + return 2; + } + + d->AddTab("Time"); + time->DrawCopy(); + + TH1 *pulse = (TH1*)arr.FindObjectInCanvas("Pulse_py", "TH1D", "Pulse"); + if (!fluct) + { + cout << "WARNING - Reading of pulse failed." << endl; + return 2; + } + + d->AddTab("Pulse"); + pulse->DrawCopy(); + + // ====================================================== + + // Load histogramm for Sum Spectrum of all Pixel + TH1 *hSum = (TH1*)arr.FindObjectInCanvas("Sum1", "TH1F", "SumHist"); + if (!hSum) + { + cout << "WARNING - Reading of hsum failed." << endl; + return 2; + } + + // Load fit function for Sum Spectrum of all Pixel + TF1 *func1 = (TF1*)arr.FindObjectInCanvas("spektrum", "TF1", "SumHist"); + if (!func1) + { + cout << "WARNING - Reading of func failed." << endl; + return 2; + } + + //Draw histogramm for Sum Spectrum of all Pixel +// d->AddTab("SumHist"); +// gPad->SetLogy(); +// gPad->SetGrid(); +// hSum->DrawCopy("hist"); +// func1->DrawCopy("SAME"); + + // Load TH2 histogramm for Pixel Spectra + TH2F *hsignal = (TH2F*)arr.FindObjectInCanvas("Signal", "TH2F", "MHSingles"); + if (!hsignal) + { + cout << "WARNING - Reading of hSpectra failed." << endl; + return 2; + } + + // Draw TH2 histogramm for Pixel Spectra + d->AddTab("Signal"); + hsignal->SetYTitle("Counts [a.u.]"); + hsignal->DrawCopy("colz"); + + + // ====================================================== + + // true: Display correctly mapped pixels in the camera displays + // but the value-vs-index plot is in software/spiral indices + // false: Display pixels in hardware/linear indices, + // but the order is the camera display is distorted + bool usemap = true; + + // map file to use (get that from La Palma!) + const char *map = usemap ? "/scratch_nfs/FACTmap111030.txt" : NULL; + + + + // ====================================================== + + if (map && gSystem->AccessPathName(map, kFileExists)) + { + gLog << "ERROR - Cannot access mapping file '" << map << "'" << endl; + return 1; + } + + // ====================================================== + + //MStatusDisplay *d = new MStatusDisplay; + + MBadPixelsCam badpixels; + badpixels.InitSize(1440); + badpixels[ 424].SetUnsuitable(MBadPixelsPix::kUnsuitable); + badpixels[ 583].SetUnsuitable(MBadPixelsPix::kUnsuitable); + badpixels[ 830].SetUnsuitable(MBadPixelsPix::kUnsuitable); + badpixels[ 923].SetUnsuitable(MBadPixelsPix::kUnsuitable); + badpixels[1208].SetUnsuitable(MBadPixelsPix::kUnsuitable); + badpixels[1399].SetUnsuitable(MBadPixelsPix::kUnsuitable); + + // Twin pixel + // 113 + // 115 + // 354 + // 423 + // 1195 + // 1393 + + // ====================================================== + + gLog << endl; + gLog.Separator("recalculation of gain from pixel spectra"); + + MTaskList tlist1; + + MParList plist1; + plist1.AddToList(&tlist1); + plist1.AddToList(&badpixels); + + MEvtLoop loop1(file); + loop1.SetDisplay(d); + loop1.SetParList(&plist1); + + // ------------------ Setup the tasks --------------- + + MRawFitsRead read1; + read1.LoadMap(map); + + MGeomApply apply1; + + MTaskInteractive mytask; + mytask.SetPreProcess(PreProcess); + mytask.SetProcess(Process); + mytask.SetPostProcess(PostProcess); + + // ------------------ Setup eventloop and run analysis --------------- + + tlist1.AddToList(&mytask); + +// if (!loop1.Eventloop(max1)) +// return 10; + + if (!loop1.GetDisplay()) + return 11; + + gStyle->SetOptFit(kTRUE); + + // ====================================================== + // ----------------------- Spectrum Fit ----------------- + // AFTER this the Code for the spektrum fit follows + // access the spektrumhistogram by the pointer *hist + + MGeomCamFACT fact; + MHCamera hcam(fact); + MHCamera hcam2(fact); + MHCamera hNormCam(fact); + + // ====================================================== + // create histograms to plot spectra and parameter distributions + + TH1F hAmplitude ("Rate", "Average number of dark counts per event", + 200, 0, 20); + TH1F hGain ("Gain", "Gain distribution", + 50, 150, 350); + TH1F hGain2 ("NormGain", "Normalized gain distribution", + 51, 0.5, 1.5); + TH1F hGainRMS ("RelSigma", "Distribution of rel. Sigma", + 100, 0, 30); + TH1F hCrosstlk ("Crosstalk", "Distribution of crosstalk probability", + 100, 0, 25); + TH1F hOffset ("Offset", "Baseline offset distribution", + 50, -7.5, 2.5); + TH1F hFitProb ("FitProb", "FitProb distribution", + 50, 0, 100); + TH1F hNoise ("Noise", "Noise distribution", + 80, -10, 30); + TH1F hXtalkPot ("XtalkPot", "Crosstalk Potential Distribution", + 200, 0, 2); + + hAmplitude.SetXTitle("Amplitude [a.u.]"); + hGain.SetXTitle("gain [a.u.]"); + hGain2.SetXTitle("gain [a.u.]"); + hGainRMS.SetXTitle("gainRMS [a.u.]"); + //hGainRMS.SetStats(false); + hCrosstlk.SetXTitle("crosstalk [a.u.]"); + hOffset.SetXTitle("baseline offset [a.u.]"); + hFitProb.SetXTitle("FitProb p-value [a.u.]"); + //hFitProb.SetStats(false); + + TH1F hSumScale("Sum2", "Signal spectrum of all pixels", 100, 0.05, 10.05); + hSumScale.SetXTitle("pulse integral [pe]"); + hSumScale.SetYTitle("Rate"); + hSumScale.SetStats(false); + hSumScale.Sumw2(); + + // define fit function for Amplitudespectrum + TF1 func("spektrum", fcn, 0, hSum->GetXaxis()->GetXmax(), 7); + func.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset", "Noise", "XTalkPot"); + func.SetLineColor(kBlue); + + // ====================================================== + + // Map for which pixel shall be plotted and which not + TArrayC usePixel(1440); + memset(usePixel.GetArray(), 1, 1440); + + // List of Pixel that should be ignored in camera view + usePixel[424] = 0; + usePixel[923] = 0; + usePixel[1208] = 0; + usePixel[583] = 0; + usePixel[830] = 0; + usePixel[1399] = 0; + usePixel[113] = 0; + usePixel[115] = 0; + usePixel[354] = 0; + usePixel[423] = 0; + usePixel[1195] = 0; + usePixel[1393] = 0; + + // ====================================================== + //-------------------_fit sum spectrum------------------- + + gLog << endl; + gLog.Separator("Fitting Sum Spectrum"); + + //Set Status Display Output and Prograssbar position + d->SetStatusLine( "...fitting Sum Spectrum", 0); + d->SetProgressBarPosition( 0.0 , 1); + + const Int_t maxbin = hSum->GetMaximumBin(); + const Double_t maxpos = hSum->GetBinCenter(maxbin); + const Double_t binwidth = hSum->GetBinWidth(maxbin); + const Double_t ampl = hSum->GetBinContent(maxbin); + + // ====================================================== + + // full width half maximum of sumspectra + double fwhmSum = 0; + + //Calculate full width half Maximum + for (int i=1; iGetBinContent(maxbin-i)+hSum->GetBinContent(maxbin+i)GetXaxis()->GetXmax(); + + //set status display prograssbar position + d->SetProgressBarPosition(0.05 , 1); + + //Fit and draw spectrum + func.SetParameters(sum_par); + func.SetRange(fitmin, fitmax); + hSum->Fit(&func, "INOQSR", "", fitmin, fitmax); + func.GetParameters(sum_par); + + //Save Parameter Errors for output + for (int i = 0; i < 7; i++){ + sum_par_err[i] = func.GetParError(i); + } + + //Set Parameters for following fits + const float gain = sum_par[1]; + const float GainRMS = sum_par[2]; + const float Crosstlk = sum_par[3]; + const float Offset = sum_par[4]; + const float Noise = sum_par[5]; + const float XtalkPot = sum_par[6]; + + cout << "--------------------" << endl; + cout << "MaxPos: " << maxpos << " \t +/- " << sum_par_err[0] << endl; + cout << "FWHM: " << fwhmSum << endl; + cout << "GAIN: " << gain << " \t +/- " << sum_par_err[1] << endl; + cout << "RMS: " << GainRMS << " \t +/- " << sum_par_err[2] << endl; + cout << "xTalk: " << Crosstlk << " \t +/- " << sum_par_err[3] << endl; + cout << "XtPot: " << XtalkPot << " \t +/- " << sum_par_err[6] << endl; + cout << "Noise: " << Noise << " \t +/- " << sum_par_err[5] << endl; + cout << "Offset: " << Offset << " \t +/- " << sum_par_err[4] << endl; + cout << "--------------------" << endl; + + gROOT->SetSelectedPad(0); + + //set status display prograssbar position + d->SetProgressBarPosition(0.1 , 1); + + // ====================================================== + + TCanvas &c11 = d->AddTab("SumHist"); + c11.cd(); + gPad->SetLogy(); + gPad->SetGridx(); + gPad->SetGridy(); + hSum->DrawCopy("hist"); + func.DrawCopy("SAME"); + + // ====================================================== + //----------Gain Calculation for each Pixel-------------- + + gLog << endl; + gLog.Separator("pixelwise gain calculation"); + + // counter for number of processed pixel + int count_ok = 0; + + // marker for pixel with abnormal behavior + bool suspicous = 0; + + // number of pixels to process + UInt_t nPixels = hsignal->GetNbinsX(); + + //set status display status + d->SetStatusLine( "...fitting pixels", 0); + + // ====================================================== + + // Loop over Pixels + for (UInt_t pixel = 0; pixel < nPixels; pixel++) + { + //set status display prograssbar position + d->SetProgressBarPosition(0.1 + 0.8*( (double)(pixel+1)/(double)nPixels ), 1); + suspicous = 0; + + // jump to next pixel if the current is marked as broken + if (usePixel[pixel]==0) + continue; + + //Projectipon of a certain Pixel out of Ampl.Spectrum + TH1D *hist = hsignal->ProjectionY("proj", pixel+1, pixel+1); + hist->SetDirectory(0); + + //Rebin Projection + hist->Rebin(2); + + // Callculate values of the bin with maximum entries + const Int_t maxbin = hist->GetMaximumBin(); + const Double_t maxPos = hist->GetBinCenter(maxbin); + + // Calculate range for the fit + const double fit_min = maxPos-GainRMS*2; + const double fit_max = maxPos+gain*(maxOrder-0.5); + + // ====================================================== + + // Set fit parameter start values + sum_par[0] = hist->GetBinContent(maxbin); +// sum_par[2] = first_gaus_RMS; +// sum_par[3] = second_gaus_ampl/first_gaus_ampl; +// sum_par[4] = -1*(second_gaus_ampl-2*first_gaus_ampl); + func.SetParameters(sum_par); +// func.FixParameter(0, first_gaus_ampl); +// func.SetParLimits(0, 0.8*sum_par[0], 1.1*sum_par[0]); +// func.SetParLimits(1, maxPos - 2*GainRMS, maxPos + 2*GainRMS); +// func.SetParLimits(2, -first_gaus_RMS, +first_gaus_RMS); +// func.FixParameter(5, sum_par[5]); +// func.SetParLimits(5, 0.7, 1.1); +// func.FixParameter(6, sum_par[6]); +// func.SetParLimits(6, 0.1, 1.9); + + // ====================================================== + // -------------- Fit Pixels spectrum ------------------- + + // Save fit result to a variable + const TFitResultPtr rc = hist->Fit(&func, "LLN0QS", "", fit_min, fit_max); + + // marker to mark if fit was successfull or not + const bool ok = int(rc)==0; + + // Save Parametervalues for statistics and bug fixing + const double fit_prob = rc->Prob(); +// const float fMaxAmpl = func.GetParameter(0); + const float fGain = func.GetParameter(1); + const float fAmplitude = func.Integral(Offset, maxOrder*fGain+Offset); //func.GetParameter(0); + const float f2GainRMS = func.GetParameter(2); + const float fCrosstlk = func.GetParameter(3); + const float fOffset = func.GetParameter(4)/30; + const float fNoise = func.GetParameter(5); + const float fXtlkPot = func.GetParameter(6); + + // Fill Parameter value hgistograms + hAmplitude.Fill(fAmplitude); + hGain.Fill(fGain); + hGain2.Fill(fGain/gain); + hCrosstlk.Fill(fCrosstlk*100); + hOffset.Fill(fOffset); + hFitProb.Fill(100*fit_prob); + hGainRMS.Fill(100*f2GainRMS/fGain); + hNoise.Fill(fNoise); + hXtalkPot.Fill(fXtlkPot); + + // Plot calculated Gain values to Camera Display + hcam.SetBinContent(pixel+1, fGain); + hcam2.SetBinContent(pixel+1, 100*f2GainRMS/fGain); + hNormCam.SetBinContent(pixel+1, fGain/gain); + + // ====================================================== + + // cancel out pixel where the fit was not succsessfull + usePixel[pixel] = ok; + + // mark pixels suspicious with negative noise + if ( fNoise < 0) + suspicous = 1; + + // mark pixels suspicious with negative GainRMS + if ( f2GainRMS < 0) + suspicous = 1; + + // mark pixels suspicious with with large deviation from mean gain + if ( fGain < gain - 1.4*GainRMS || fGain > gain + 1.4*GainRMS) + suspicous = 1; + + //plott especialy marked pixel + if (count_ok<3 || suspicous == 1 || !ok) + { + hist->SetName("Pix%d"); + TCanvas &c = d->AddTab(Form("Pix%d", pixel)); + c.cd(); + gPad->SetLogy(); + gPad->SetGridx(); + gPad->SetGridy(); + + rc->Print("V"); + +// hist->SetStats(false); + hist->SetYTitle("Counts [a.u.]"); + hist->SetName(Form("Pix%d", pixel)); + hist->DrawCopy()->SetDirectory(0); + func.DrawCopy("SAME")->SetRange(fit_min, fit_max); + } + + if (!ok) + { + cout << "Pixel " << setw(4) << pixel << ": failed!" << endl; + delete hist; + continue; + } + + // ====================================================== + + //Fill sum spectrum + for (int b=1; b<=hist->GetNbinsX(); b++) + hSumScale.Fill( + (hist->GetBinCenter(b)-fOffset)/fGain, + hist->GetBinContent(b) + ); + + count_ok++; + + delete hist; + } + + // define a general gaus function + TF1 fgaus("fgaus", "gaus"); + hGain.Fit(&fgaus, "N0QS"); + + // ====================================================== + //------------------Draw histograms---------------------- + + gLog << endl; + gLog.Separator("Drawing Histograms"); + + d->SetStatusLine( "...drawing Histograms", 0); + d->SetProgressBarPosition(0.9 , 1); + + // cancel out pixel in camera view + hcam.SetUsed(usePixel); + hcam2.SetUsed(usePixel); + + // ====================================================== + + TCanvas &canv = d->AddTab("Gain"); + canv.cd(); + canv.Divide(2,2); + + canv.cd(1); + hcam.SetNameTitle("gain", "Gain distribution over whole camera"); + hcam.DrawCopy(); + + canv.cd(2); + hcam2.SetNameTitle("gainRMS", "GainRMS distribution over whole camera"); + hcam2.DrawCopy(); + + TCanvas &cam_canv = d->AddTab("Camera_Gain"); + + // ====================================================== + //---------- Draw gain corrected sum specetrum ---------- + + gROOT->SetSelectedPad(0); + TCanvas &c12 = d->AddTab("GainHist"); + c12.cd(); + gPad->SetLogy(); + gPad->SetGridx(); + gPad->SetGridy(); + + hSumScale.DrawCopy("hist"); + + //---------fit gain corrected sum spectrum --------------- + + // Callculate values of the bin with maximum entries + const Int_t maxbin2 = hSumScale.GetMaximumBin(); + const Double_t ampl2 = hSumScale.GetBinContent(maxbin2); + + //Set fit start parameters + Double_t par2[7] = + { + ampl2, 1, GainRMS/fgaus.GetParameter(1), Crosstlk, 0, 0, 1 + }; + func.SetParameters(par2); + + // fix gain to 1 as it was normalized allready + func.FixParameter(1, 1); + + // fix offset parameter + // func.FixParameter(4, 0); + + // set progressbar possition + d->SetStatusLine( "...fitting corr. spectrum", 0); + + // Fit corrected Spectrum with Likelyhood Method + const TFitResultPtr rc = hSumScale.Fit(&func, "LLEN0QS", "", 0.75, 10); + + func.DrawCopy("SAME")->SetRange(0.5, 10); + + // Print fit results + rc->Print("V"); + + // set progressbar possition + d->SetProgressBarPosition(0.95 , 1); + + // ====================================================== + + TCanvas &c2 = d->AddTab("Result"); + c2.Divide(3,2); + c2.cd(1); + gPad->SetLogy(); + hXtalkPot.DrawCopy(); + + c2.cd(2); + gPad->SetLogy(); + hGain.DrawCopy(); + + //plot normailzed camera view + cam_canv.cd(); + hNormCam.SetStats(false); + + hNormCam.SetNameTitle("GainCam", "Normalized gain distribution"); + hNormCam.SetUsed(usePixel); + hNormCam.DrawCopy(); + + c2.cd(3); + gPad->SetLogy(); + hGainRMS.DrawCopy(); + + c2.cd(4); + gPad->SetLogy(); + hCrosstlk.DrawCopy(); + + c2.cd(5); + gPad->SetLogy(); + hOffset.DrawCopy(); + + c2.cd(6); + gPad->SetLogy(); + hNoise.DrawCopy(); + + // ====================================================== + + TCanvas &c3 =d->AddTab("Prop"); + c3.Divide(2,1); + c3.cd(1); + gPad->SetGrid(); + gPad->SetLogy(); + hFitProb.DrawCopy(); + + c3.cd(2); + gPad->SetGrid(); + gPad->SetLogy(); + hAmplitude.DrawCopy(); + + // ====================================================== + + d->AddTab("Norm"); + gPad->SetGrid(); + gPad->SetLogy(); + TH1 *hh = hGain2.DrawCopy(); + hh->Fit("gaus"); + + // ====================================================== + + //save results to file + if (save){ + d->SetStatusLine( "...saving results to rootfile", 0); + d->SaveAs(outfile); + cout << "..success!" << endl; + } + d->SetProgressBarPosition(1 , 1); + d->SetStatusLine( "...done", 0); + return 0; +} + +Double_t fcn(Double_t *xx, Double_t *par) +{ + Double_t x = xx[0]; + + Double_t ampl = par[0]; + Double_t gain = par[1]; + Double_t sigma = par[2]; + Double_t cross = par[3]; + Double_t shift = par[4]; + Double_t noise = par[5]; + Double_t crossPot = par[6]; + + Double_t xTalk = 1; + + Double_t y = 0; + for (int order = 1; order < 14; order++) + { + Double_t arg = (x - order*gain - shift)/(order==1?sigma+noise:sigma); + + Double_t gauss = TMath::Exp(-0.5 * arg * arg/order); + + y += xTalk*gauss; + xTalk = Power(cross, Power(order, crossPot) ); + } + + return ampl*y; +} + +// end of PlotMedianEachSliceOfPulse +//----------------------------------------------------------------------------