Index: /fact/tools/marsmacros/singlepe.C =================================================================== --- /fact/tools/marsmacros/singlepe.C (revision 14244) +++ /fact/tools/marsmacros/singlepe.C (revision 14245) @@ -44,6 +44,6 @@ struct Single { - float fSignal; - float fTime; + float fSignal; + float fTime; }; @@ -54,26 +54,41 @@ vector> fData; -public: + public: MSingles(const char *name=NULL, const char *title=NULL) : fIntegrationWindow(0) { - fName = name ? name : "MSingles"; + 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; } + 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; + return kTRUE; } void DrawPixelContent(Int_t num) const { } @@ -84,22 +99,26 @@ class MH2F : public TH2F { -public: + 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; + 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); @@ -108,26 +127,30 @@ class MProfile2D : public TProfile2D { -public: + 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; + 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); @@ -147,136 +170,155 @@ MPedestalSubtractedEvt *fEvent; //! -public: + 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); + 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; + 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; + 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; jGetSamples(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"); + TVirtualPad *pad = gPad ? gPad : MakeDefCanvas(this); + + AppendPad(""); + + pad->Divide(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"); } @@ -295,118 +337,130 @@ 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; + 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() { - const UInt_t roi = fEvent->GetNumSamples(); - - const size_t navg = 10; - const float threshold = 5; - const UInt_t start_skip = 20; - const UInt_t end_skip = 10; - const UInt_t integration_size = 10*3; - const UInt_t max_search_window = 30; - - vector val(roi-navg);//result of first sliding average - for (int pix=0; pixGetNumPixels(); pix++) - { - vector &result = fSingles->GetVector(pix); - result.clear(); - - const Float_t *ptr = fEvent->GetSamples(pix); - - //Sliding window - for (size_t i=0; ithreshold || - val[i+4]>threshold || - - val[i+5]GetNumSamples(); + + const size_t navg = 10; + const float threshold = 5; + const UInt_t start_skip = 20; + const UInt_t end_skip = 10; + const UInt_t integration_size = 10*3; + const UInt_t max_search_window = 30; + + vector val(roi-navg);//result of first sliding average + + for (int pix=0; pixGetNumPixels(); pix++) + { + vector &result = fSingles->GetVector(pix); + result.clear(); + + const Float_t *ptr = fEvent->GetSamples(pix); + + //Sliding window + for (size_t i=0; ithreshold || + val[i+4]>threshold || + + val[i+5] val[k_max]) - k_max = k; + if (val[k] > val[k_max]) + k_max = k; } - if (k_max == i+5 || k_max == i + max_search_window-1) continue; - - //search for half maximum before maximum - UInt_t k_half_max = 0; - for (UInt_t k=k_max; k>k_max-25; k--) + if (k_max == i+5 || k_max == i + max_search_window-1) continue; + + //search for half maximum before maximum + UInt_t k_half_max = 0; + + for (UInt_t k=k_max; k>k_max-25; k--) { - if (val[k-1] < val[k_max]/2 && - val[k] >= val[k_max]/2 ) + if (val[k-1] < val[k_max]/2 && + val[k] >= val[k_max]/2 ) { - k_half_max = k; - break; + k_half_max = k; + break; } } - if (k_half_max > roi-navg-end_skip-max_search_window - integration_size) - continue; - if (k_half_max == 0) - continue; - if (k_max - k_half_max > 14) - continue; - - fluct2->Fill(k_max - k_half_max); - - // Evaluate arrival time more precisely!!! - // Get a better integration window - - Single single; - single.fSignal = 0; - single.fTime = k_half_max; - - - // Crossing of the threshold is at 5 - for (UInt_t j=k_half_max; j roi-navg-end_skip-max_search_window - integration_size) + continue; + + if (k_half_max == 0) + continue; + + if (k_max - k_half_max > 14) + continue; + + fluct2->Fill(k_max - k_half_max); + + // Evaluate arrival time more precisely!!! + // Get a better integration window + + Single single; + single.fSignal = 0; + single.fTime = k_half_max; + + + // Crossing of the threshold is at 5 + for (UInt_t j=k_half_max; jDirName(runfile), gSystem->BaseName(runfile)); @@ -449,603 +503,629 @@ // iter.AddDirectory("/fact/raw/2012/06/25", "20120625_119.fits.gz"); - // ====================================================== - - // 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; - - // ------------------------------------------------------ - - Long_t max1 = 0; - - // ====================================================== - - if (map && gSystem->AccessPathName(map, kFileExists)) - { - gLog << "ERROR - Cannot access mapping file '" << map << "'" << endl; - return 1; - } - if (gSystem->AccessPathName(drsfile, kFileExists)) - { - gLog << "ERROR - Cannot access drsfile file '" << drsfile << "'" << endl; - return 2; - } - - // -------------------------------------------------------------------------------- - - gLog.Separator("Singles"); - gLog << "Extract singles with '" << drsfile << "'" << endl; - gLog << endl; - - // ------------------------------------------------------ - - MDrsCalibration drscalib300; - if (!drscalib300.ReadFits(drsfile)) - return 4; - - // ------------------------------------------------------ - - iter.Sort(); - iter.Print(); - - // ====================================================== - - //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 - - //MDrsCalibrationTime timecam; - - MH::SetPalette(); - - // ====================================================== - - gLog << endl; - gLog.Separator("Processing external light pulser run"); - - MTaskList tlist1; - - MSingles singles; - MRawRunHeader header; - - MParList plist1; - plist1.AddToList(&tlist1); - plist1.AddToList(&drscalib300); - plist1.AddToList(&badpixels); - plist1.AddToList(&singles); - plist1.AddToList(&header); - - TString Title; - Title = iter.Next(); - iter.Reset(); - Title += "; "; - Title += max1; - - MEvtLoop loop1(Title); - loop1.SetDisplay(d); - loop1.SetParList(&plist1); - - // ------------------ Setup the tasks --------------- - - MRawFitsRead read1; - read1.LoadMap(map); - read1.AddFiles(iter); - - MContinue cont1("(MRawEvtHeader.GetTriggerID&0xff00)!=0x100", "SelectCal"); - - MGeomApply apply1; - - MDrsCalibApply drsapply1; - - MTaskInteractive mytask; - mytask.SetPreProcess(PreProcess); - mytask.SetProcess(Process); - mytask.SetPostProcess(PostProcess); - - MFillH fill("MHSingles"); - - // ------------------ Setup eventloop and run analysis --------------- - - tlist1.AddToList(&read1); + // ====================================================== + + // 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; + + // ------------------------------------------------------ + + Long_t max1 = 0; + + // ====================================================== + + if (map && gSystem->AccessPathName(map, kFileExists)) + { + gLog << "ERROR - Cannot access mapping file '" << map << "'" << endl; + return 1; + } + + if (gSystem->AccessPathName(drsfile, kFileExists)) + { + gLog << "ERROR - Cannot access drsfile file '" << drsfile << "'" << endl; + return 2; + } + + // -------------------------------------------------------------------------------- + + gLog.Separator("Singles"); + gLog << "Extract singles with '" << drsfile << "'" << endl; + gLog << endl; + + // ------------------------------------------------------ + + MDrsCalibration drscalib300; + + if (!drscalib300.ReadFits(drsfile)) + return 4; + + // ------------------------------------------------------ + + iter.Sort(); + iter.Print(); + + // ====================================================== + + //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 + + //MDrsCalibrationTime timecam; + + MH::SetPalette(); + + // ====================================================== + + gLog << endl; + gLog.Separator("Processing external light pulser run"); + + MTaskList tlist1; + + MSingles singles; + MRawRunHeader header; + + MParList plist1; + plist1.AddToList(&tlist1); + plist1.AddToList(&drscalib300); + plist1.AddToList(&badpixels); + plist1.AddToList(&singles); + plist1.AddToList(&header); + + TString Title; + Title = iter.Next(); + iter.Reset(); + Title += "; "; + Title += max1; + + MEvtLoop loop1(Title); + loop1.SetDisplay(d); + loop1.SetParList(&plist1); + + // ------------------ Setup the tasks --------------- + + MRawFitsRead read1; + read1.LoadMap(map); + read1.AddFiles(iter); + + MContinue cont1("(MRawEvtHeader.GetTriggerID&0xff00)!=0x100", "SelectCal"); + + MGeomApply apply1; + + MDrsCalibApply drsapply1; + + MTaskInteractive mytask; + mytask.SetPreProcess(PreProcess); + mytask.SetProcess(Process); + mytask.SetPostProcess(PostProcess); + + MFillH fill("MHSingles"); + + // ------------------ Setup eventloop and run analysis --------------- + + tlist1.AddToList(&read1); // tlist1.AddToList(&cont1); - tlist1.AddToList(&apply1); - tlist1.AddToList(&drsapply1); - tlist1.AddToList(&mytask); - tlist1.AddToList(&fill); - - if (!loop1.Eventloop(max1)) - return 10; - - if (!loop1.GetDisplay()) - return 11; - - gStyle->SetOptFit(kTRUE); - - MHSingles* h = (MHSingles*) plist1.FindObject("MHSingles"); - if (!h) - return 0; - - TH2 *hsignal = h->GetSignal(); - TH2 *htime = h->GetTime(); - TH2 *hpulse = h->GetPulse(); - - d->AddTab("Time"); - TH1 *ht = htime->ProjectionY(); - ht->SetYTitle("Counts [a.u.]"); - ht->Scale(1./1440); - ht->Draw(); - - d->AddTab("Pulse"); - TH1 *hp = hpulse->ProjectionY(); - hp->SetYTitle("Counts [a.u.]"); - hp->Scale(1./1440); - hp->Draw(); - - - // ------------------ Spectrum Fit --------------- - // AFTER this the Code for the spektrum fit follows - // access the spektrumhistogram by the pointer *hist - UInt_t maxOrder = 8; - - - MGeomCamFACT fact; - MHCamera hcam(fact); - MHCamera hcam2(fact); - - //built an array of Amplitude spectra - TH1F hAmplitude ("Rate", "Average number of dark counts per event", - 200, 0, 20); - hAmplitude.SetXTitle("Amplitude [a.u.]"); - hAmplitude.SetYTitle("Rate"); - - TH1F hGain ("Gain", "Gain distribution", - 50, 100, 300); - hGain.SetXTitle("gain [a.u.]"); - hGain.SetYTitle("Rate"); - - TH1F hGainRMS ("RelSigma", "Distribution of rel. Sigma", - 100, 0, 30); - hGainRMS.SetXTitle("gainRMS [a.u.]"); - hGainRMS.SetYTitle("Rate"); - hGainRMS.SetStats(false); - - TH1F hCrosstlk ("Crosstalk", "Distribution of crosstalk probability", - 100, 0, 25); - hCrosstlk.SetXTitle("crosstalk [a.u.]"); - hCrosstlk.SetYTitle("Rate"); - - TH1F hOffset ("Offset", "Baseline offset distribution", - 50, -7.5, 2.5); - hOffset.SetXTitle("baseline offset [a.u.]"); - hOffset.SetYTitle("Rate"); - - TH1F hFitProb ("FitProb", "FitProb distribution", - 50, 0, 100); - hFitProb.SetXTitle("FitProb p-value [a.u.]"); - hFitProb.SetYTitle("Rate"); - hFitProb.SetStats(false); - - - TH1F hSum ("Sum1", "Sum of all pixel's' integral spectra", - 100, 0, 2000); - hSum.SetXTitle("pulse integral [a.u.]"); - hSum.SetYTitle("Rate"); - hSum.SetStats(false); - hSum.Sumw2(); - - - TH1F hSumScale ("Sum2", - "Sum of all pixel's' integral spectra (scaled with gain)", - 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, 1000, 5); - func.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset"); - func.SetLineColor(kBlue); - - // define fit function for Amplitudespectrum - TF1 funcSum("sum_spektrum", fcn, 0, 2000, 5); - funcSum.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset", "PowOfXtalk"); - funcSum.SetLineColor(kBlue); - - // define fit function for Amplitudespectrum - TF1 funcScaled("gain_sum_spektrum", fcn, 0, 10, 5); - funcScaled.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset"); - funcScaled.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; + tlist1.AddToList(&apply1); + tlist1.AddToList(&drsapply1); + tlist1.AddToList(&mytask); + tlist1.AddToList(&fill); + + if (!loop1.Eventloop(max1)) + return 10; + + if (!loop1.GetDisplay()) + return 11; + + gStyle->SetOptFit(kTRUE); + + MHSingles* h = (MHSingles*) plist1.FindObject("MHSingles"); + + if (!h) + return 0; + + TH2 *hsignal = h->GetSignal(); + TH2 *htime = h->GetTime(); + TH2 *hpulse = h->GetPulse(); + + d->AddTab("Time"); + TH1 *ht = htime->ProjectionY(); + ht->SetYTitle("Counts [a.u.]"); + ht->Scale(1./1440); + ht->Draw(); + + d->AddTab("Pulse"); + TH1 *hp = hpulse->ProjectionY(); + hp->SetYTitle("Counts [a.u.]"); + hp->Scale(1./1440); + hp->Draw(); + + + // ------------------ Spectrum Fit --------------- + // AFTER this the Code for the spektrum fit follows + // access the spektrumhistogram by the pointer *hist + UInt_t maxOrder = 8; + + + MGeomCamFACT fact; + MHCamera hcam(fact); + MHCamera hcam2(fact); + + //built an array of Amplitude spectra + TH1F hAmplitude ("Rate", "Average number of dark counts per event", + 200, 0, 20); + hAmplitude.SetXTitle("Amplitude [a.u.]"); + hAmplitude.SetYTitle("Rate"); + + TH1F hGain ("Gain", "Gain distribution", + 400, 100, 300); + hGain.SetXTitle("gain [a.u.]"); + hGain.SetYTitle("Rate"); + + TH1F hGainRMS ("RelSigma", "Distribution of rel. Sigma", + 100, 0, 30); + hGainRMS.SetXTitle("gainRMS [a.u.]"); + hGainRMS.SetYTitle("Rate"); + hGainRMS.SetStats(false); + + TH1F hCrosstlk ("Crosstalk", "Distribution of crosstalk probability", + 100, 0, 25); + hCrosstlk.SetXTitle("crosstalk [a.u.]"); + hCrosstlk.SetYTitle("Rate"); + + TH1F hOffset ("Offset", "Baseline offset distribution", + 50, -7.5, 2.5); + hOffset.SetXTitle("baseline offset [a.u.]"); + hOffset.SetYTitle("Rate"); + + TH1F hFitProb ("FitProb", "FitProb distribution", + 50, 0, 100); + hFitProb.SetXTitle("FitProb p-value [a.u.]"); + hFitProb.SetYTitle("Rate"); + hFitProb.SetStats(false); + + + TH1F hSum ("Sum1", "Sum of all pixel's' integral spectra", + 100, 0, 2000); + hSum.SetXTitle("pulse integral [a.u.]"); + hSum.SetYTitle("Rate"); + hSum.SetStats(false); + hSum.Sumw2(); + + + TH1F hSumScale ("Sum2", + "Sum of all pixel's' integral spectra (scaled with gain)", + 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, 1000, 5); + func.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset"); + func.SetLineColor(kBlue); + + // define fit function for Amplitudespectrum + TF1 funcSum("sum_spektrum", fcn, 0, 2000, 5); + funcSum.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset", "PowOfXtalk"); + funcSum.SetLineColor(kBlue); + + // define fit function for Amplitudespectrum + TF1 funcScaled("gain_sum_spektrum", fcn, 0, 10, 5); + funcScaled.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset"); + funcScaled.SetLineColor(kBlue); + + TF1 funcScaledBL("gain_sum_spektrum", fcn, 0, 10, 5); + funcScaledBL.SetParNames("Maximum", "Gain", "Sigma", "XtalkProb", "Offset"); + funcScaledBL.SetLineColor(kRed); + + + + // 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; // usePixel[990] = 0; - // Map for which pixel which were suspicous - bool suspicous[1440]; - for (int i=0; i<1440; i++) suspicous[i]=false; - - // List of Pixel that should be ignored in camera view - // suspicous[802] = true; - - //------------------------------------------------------------------------ - // Fill and calculate sum spectrum - - // Begin of Loop over Pixels - for (UInt_t pixel = 0; pixel < singles.GetNumPixels(); pixel++) - { - //jump to next pixel if the current is marked as faulty - if (usePixel[pixel]==0) - continue; - - TH1D *hist = hsignal->ProjectionY("proj", pixel+1, pixel+1); - hist->SetDirectory(0); - //loop over slices - for (int b=1; b<=hist->GetNbinsX(); b++) - { - //Fill sum spectrum - hSum.Fill(hist->GetBinCenter(b), hist->GetBinContent(b)); - } - } - for(UInt_t bin = 1; bin < (UInt_t)hSum.GetNbinsX()+1; bin++) - { - hSum.SetBinError(bin, hSum.GetBinContent(bin)*0.1); - } - - gROOT->SetSelectedPad(0); - TCanvas &c11 = d->AddTab("SumHist"); - c11.cd(); - gPad->SetLogy(); - gPad->SetGridx(); - gPad->SetGridy(); - hSum.DrawCopy("hist"); - //------------------------fit sum spectrum------------------------------- - const Int_t maxbin = hSum.GetMaximumBin(); - const Double_t binwidth = hSum.GetBinWidth(maxbin); - const Double_t ampl = hSum.GetBinContent(maxbin); - - double fwhmSum = 0; - - //Calculate full width half Maximum - for (int i=1; iProjectionY("proj", pixel+1, pixel+1); + hist->SetDirectory(0); + + //loop over slices + for (int b=1; b<=hist->GetNbinsX(); b++) + { + //Fill sum spectrum + hSum.Fill(hist->GetBinCenter(b), hist->GetBinContent(b)); + } + } + + for (UInt_t bin = 1; bin < (UInt_t)hSum.GetNbinsX()+1; bin++) + { + hSum.SetBinError(bin, hSum.GetBinContent(bin)*0.1); + } + + gROOT->SetSelectedPad(0); + TCanvas &c11 = d->AddTab("SumHist"); + c11.cd(); + gPad->SetLogy(); + gPad->SetGridx(); + gPad->SetGridy(); + hSum.DrawCopy("hist"); + //------------------------fit sum spectrum------------------------------- + const Int_t maxbin = hSum.GetMaximumBin(); + const Double_t binwidth = hSum.GetBinWidth(maxbin); + const Double_t ampl = hSum.GetBinContent(maxbin); + + double fwhmSum = 0; + + //Calculate full width half Maximum + for (int i=1; iProjectionY("proj", pixel+1, pixel+1); - hist->SetDirectory(0); - - for(UInt_t bin = 1; bin < (UInt_t)hSum.GetNbinsX()+1; bin++) - { - hist->SetBinError(bin, hSum.GetBinContent(bin)*0.1); - } - - const Int_t maxbin = hist->GetMaximumBin(); - const Double_t binwidth = hist->GetBinWidth(maxbin); - const Double_t ampl = hist->GetBinContent(maxbin); - const Double_t maxPos = hist->GetBinCenter(maxbin); - - double fwhm = 0; - for (int i=1; iGetBinContent(maxbin-i)+hist->GetBinContent(maxbin+i)GetBinCenter(maxbin), + // EOF: Fill and calculate sum spectrum + //------------------------------------------------------------------------ + + + + //------------------------------------------------------------------------ + // Gain Calculation for each Pixel + + int count_ok = 0; + + // Begin of Loop over Pixels + for (UInt_t pixel = 0; pixel < singles.GetNumPixels(); pixel++) + { + + 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); + + for (UInt_t bin = 1; bin < (UInt_t)hSum.GetNbinsX()+1; bin++) + { + hist->SetBinError(bin, hSum.GetBinContent(bin)*0.1); + } + + const Int_t maxbin = hist->GetMaximumBin(); + + const Double_t binwidth = hist->GetBinWidth(maxbin); + + const Double_t ampl = hist->GetBinContent(maxbin); + + const Double_t maxPos = hist->GetBinCenter(maxbin); + + double fwhm = 0; + + for (int i=1; iGetBinContent(maxbin-i)+hist->GetBinContent(maxbin+i)GetBinCenter(maxbin), // gain, - fwhm*2, + fwhm*2, // GainRMS, - Crosstlk, - Offset - }; - - - func.SetParameters(par); - const double fit_min = par[1]-fwhm*1.4; - const double fit_max = par[1]*maxOrder; - func.SetRange(fit_min-fwhm*2, fit_max); - func.SetParLimits(1, maxPos - fwhm, maxPos + fwhm); - - if (suspicous[pixel] == true) - { - cout << "Maxbin\t" << maxbin << endl; - cout << "min\t" << fit_min << endl; - cout << "max\t" << fit_max << endl; - cout << "Amplitude, 1.Max, Sigma, fwhm:\t" - << par[0] << "\t" - << par[1] << "\t" - << par[2] << "\t" - << fwhm << "\t" << endl; - } - - //Rebin Projection - hist->Rebin(2); - - //Fit Pixels spectrum - const TFitResultPtr rc = hist->Fit(&func, "N0QS", "", fit_min, fit_max); - - const bool ok = int(rc)==0; - - const double fit_prob = rc->Prob(); - const float fGain = func.GetParameter(1); - const float fAmplitude = func.Integral(0, 7*fGain)/h->GetNumEntries(); //func.GetParameter(0); - const float f2GainRMS = func.GetParameter(2); - const float fCrosstlk = func.GetParameter(3); - const float fOffset = func.GetParameter(4)/singles.GetIntegrationWindow(); + Crosstlk, + Offset + }; + + + func.SetParameters(par); + const double fit_min = par[1]-fwhm*1.4; + const double fit_max = par[1]*maxOrder; + func.SetRange(fit_min-fwhm*2, fit_max); + func.SetParLimits(1, maxPos - fwhm, maxPos + fwhm); + + if (suspicous[pixel] == true) + { + cout << "Maxbin\t" << maxbin << endl; + cout << "min\t" << fit_min << endl; + cout << "max\t" << fit_max << endl; + cout << "Amplitude, 1.Max, Sigma, fwhm:\t" + << par[0] << "\t" + << par[1] << "\t" + << par[2] << "\t" + << fwhm << "\t" << endl; + } + + //Rebin Projection + hist->Rebin(2); + + //Fit Pixels spectrum + const TFitResultPtr rc = hist->Fit(&func, "N0QS", "", fit_min, fit_max); + + const bool ok = int(rc)==0; + + const double fit_prob = rc->Prob(); + const float fGain = func.GetParameter(1); + const float fAmplitude = func.Integral(0, 7*fGain)/h->GetNumEntries(); //func.GetParameter(0); + const float f2GainRMS = func.GetParameter(2); + const float fCrosstlk = func.GetParameter(3); + const float fOffset = func.GetParameter(4)/singles.GetIntegrationWindow(); // const float fCrossOrd = func.GetParameter(5); - hAmplitude.Fill(fAmplitude); - hGain.Fill(fGain); - //hGainRMS.Fill(f2GainRMS); - hCrosstlk.Fill(fCrosstlk*100); - hOffset.Fill(fOffset); - hFitProb.Fill(100*fit_prob); - hGainRMS.Fill(100*f2GainRMS/fGain); - - hcam.SetBinContent(pixel+1, fGain); - hcam2.SetBinContent(pixel+1, 100*f2GainRMS/fGain); - - usePixel[pixel] = ok; - - if (!ok) - { - cout << "Pixel " << setw(4) << pixel << ": failed!" << endl; - continue; - } - - //Fill sum spectrum - for (int b=1; b<=hist->GetNbinsX(); b++) + hAmplitude.Fill(fAmplitude); + hGain.Fill(fGain); + //hGainRMS.Fill(f2GainRMS); + hCrosstlk.Fill(fCrosstlk*100); + hOffset.Fill(fOffset); + hFitProb.Fill(100*fit_prob); + hGainRMS.Fill(100*f2GainRMS/fGain); + + hcam.SetBinContent(pixel+1, fGain); + hcam2.SetBinContent(pixel+1, 100*f2GainRMS/fGain); + + usePixel[pixel] = ok; + + if (!ok) + { + cout << "Pixel " << setw(4) << pixel << ": failed!" << endl; + continue; + } + + //Fill sum spectrum + for (int b=1; b<=hist->GetNbinsX(); b++) { // hSum.Fill(hist->GetBinCenter(b), hist->GetBinContent(b)); - hSumScale.Fill((hist->GetBinCenter(b)-fOffset)/fGain, (hist->GetBinContent(b))); - } - - //plott special pixel - if ( - count_ok==0 || - suspicous[pixel] - ) - { - TCanvas &c = d->AddTab(Form("Pix%d", pixel)); - c.cd(); - gPad->SetLogy(); - gPad->SetGridx(); - gPad->SetGridy(); - - hist->SetStats(false); - hist->SetYTitle("Counts [a.u.]"); - hist->SetName(Form("Pix%d", pixel)); - hist->DrawCopy("hist")->SetDirectory(0); + hSumScale.Fill((hist->GetBinCenter(b)-fOffset)/fGain, (hist->GetBinContent(b))); + } + + //plott special pixel + if ( + count_ok==0 || + suspicous[pixel] + ) + { + TCanvas &c = d->AddTab(Form("Pix%d", pixel)); + c.cd(); + gPad->SetLogy(); + gPad->SetGridx(); + gPad->SetGridy(); + + hist->SetStats(false); + hist->SetYTitle("Counts [a.u.]"); + hist->SetName(Form("Pix%d", pixel)); + hist->DrawCopy("hist")->SetDirectory(0); // hist->Draw(); - func.DrawCopy("SAME"); - } - - count_ok++; - - delete hist; + func.DrawCopy("SAME"); + } + + count_ok++; + + delete hist; // if (suspicous[pixel] == true) usePixel[pixel]=0; } - // End of Loop over Pixel - - //------------------fill histograms----------------------- - - 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(); - - for(UInt_t bin = 1; bin < (UInt_t)hSum.GetNbinsX()+1; bin++) - { - hSumScale.SetBinError(bin, hSum.GetBinContent(bin)*0.1); - } - hSumScale.DrawCopy("hist"); - - - //-------------------- fit gain corrected sum spectrum ------------------- - const Int_t maxbin2 = hSumScale.GetMaximumBin(); - const Double_t binwidth2 = hSumScale.GetBinWidth(maxbin2); - const Double_t ampl2 = hSumScale.GetBinContent(maxbin2); - - //Calculate full width half Maximum - double fwhmScaled = 0; - for (int i=1; iAddTab("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(); + + for (UInt_t bin = 1; bin < (UInt_t)hSum.GetNbinsX()+1; bin++) + { + hSumScale.SetBinError(bin, hSum.GetBinContent(bin)*0.1); + } + + hSumScale.DrawCopy("hist"); + + + //-------------------- fit gain corrected sum spectrum ------------------- + const Int_t maxbin2 = hSumScale.GetMaximumBin(); + const Double_t binwidth2 = hSumScale.GetBinWidth(maxbin2); + const Double_t ampl2 = hSumScale.GetBinContent(maxbin2); + + //Calculate full width half Maximum + double fwhmScaled = 0; + + for (int i=1; iAddTab("PeakHeights"); // gPad->SetLogy(); - hMaxHeights.SetLineColor(kRed); + hMaxHeights.SetLineColor(kRed); // hMaxHeights.DrawCopy("SAME"); - TF1 fExpo( "fExpo","expo", 0, maxOrder-2); - fExpo.SetLineColor(kRed); - hMaxHeights.Fit(&fExpo, "N0QS" ); + TF1 fExpo( "fExpo","expo", 0, maxOrder-2); + fExpo.SetLineColor(kRed); + hMaxHeights.Fit(&fExpo, "N0QS" ); // fExpo.DrawCopy("SAME"); // c12.cd(); // fExpo.DrawCopy("SAME"); - TCanvas &c2 = d->AddTab("Result"); - c2.Divide(3,2); - c2.cd(1); - gPad->SetLogy(); - hAmplitude.DrawCopy(); - - c2.cd(2); - gPad->SetLogy(); + TCanvas &c2 = d->AddTab("Result"); + c2.Divide(3,2); + c2.cd(1); + gPad->SetLogy(); + hAmplitude.DrawCopy(); + + c2.cd(2); + gPad->SetLogy(); // hGain.DrawCopy(); - TF1 GainGaus( "GainGaus", "gaus", 0, hGain.GetXaxis()->GetXmax()); - hGain.Rebin(2); - hGain.Fit(&GainGaus, "N0QS"); + TF1 GainGaus( "GainGaus", "gaus", 0, hGain.GetXaxis()->GetXmax()); + hGain.Rebin(2); + hGain.Fit(&GainGaus, "N0QS"); // float gain_mean = GainGaus.GetParameter(1); - float gain_median = MedianOfH1(&hGain); - TH1F hNormGain ("NormGain", "median normalized Gain distribution", - hGain.GetNbinsX(), - (hGain.GetXaxis()->GetXmin())/gain_median, - (hGain.GetXaxis()->GetXmax())/gain_median - ); - hNormGain.SetXTitle("gain [fraction of median gain value]"); - hNormGain.SetYTitle("Counts"); - - - hNormGain.Add(&hGain); + float gain_median = MedianOfH1(&hGain); + TH1F hNormGain ("NormGain", "median normalized Gain distribution", + hGain.GetNbinsX(), + (hGain.GetXaxis()->GetXmin())/gain_median, + (hGain.GetXaxis()->GetXmax())/gain_median + ); + hNormGain.SetXTitle("gain [fraction of median gain value]"); + hNormGain.SetYTitle("Counts"); + + + hNormGain.Add(&hGain); // hNormGain.SetStats(false); - hNormGain.GetXaxis()->SetRangeUser( - 1-(hNormGain.GetXaxis()->GetXmax()-1), - hNormGain.GetXaxis()->GetXmax() - ); - hNormGain.DrawCopy(); - hNormGain.Fit(&GainGaus, "N0QS"); - GainGaus.DrawCopy("SAME"); - - //plot normailzed camera view - cam_canv.cd(); - MHCamera hNormCam(fact); - hNormCam.SetStats(false); - for (UInt_t pixel =1; pixel <= 1440; pixel++) - { - hNormCam.SetBinContent(pixel, hcam.GetBinContent(pixel)/gain_median); - } - - hNormCam.SetNameTitle( "gain","normalized [Median] Gain distribution over whole camera"); - hNormCam.SetZTitle("Horst"); - hNormCam.SetUsed(usePixel); - hNormCam.DrawCopy(); + hNormGain.GetXaxis()->SetRangeUser( + 1-(hNormGain.GetXaxis()->GetXmax()-1), + hNormGain.GetXaxis()->GetXmax() + ); + hNormGain.DrawCopy(); + hNormGain.Fit(&GainGaus, "N0QS"); + GainGaus.DrawCopy("SAME"); + + //plot normailzed camera view + cam_canv.cd(); + MHCamera hNormCam(fact); + hNormCam.SetStats(false); + + for (UInt_t pixel =1; pixel <= 1440; pixel++) + { + hNormCam.SetBinContent(pixel, hcam.GetBinContent(pixel)/gain_median); + } + + hNormCam.SetNameTitle( "gain","normalized [Median] Gain distribution over whole camera"); + hNormCam.SetZTitle("Horst"); + hNormCam.SetUsed(usePixel); + hNormCam.DrawCopy(); // hGain.Scale(1/GainGaus.GetParameter(1)); // GainGaus.DrawCopy("SAME"); - 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(); - hFitProb.DrawCopy(); - - - TCanvas &c25 = d->AddTab("Spectra"); - c25.Divide(2,1); - c25.cd(1); - gPad->SetLogy(); - gPad->SetGrid(); - hSum.DrawCopy("hist"); - funcSum.DrawCopy("SAME"); - - c25.cd(2); - gPad->SetLogy(); - gPad->SetGrid(); - hSumScale.DrawCopy("hist"); - funcScaled.DrawCopy("SAME"); - - - /* - TCanvas &c3 = d->AddTab("Separation"); - gPad->SetLogy(); - hSep.DrawCopy(); - */ - -// d->SaveAs(outfile); - return 0; + 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(); +// hFitProb.DrawCopy(); + hGain.DrawCopy(); + + + TCanvas &c25 = d->AddTab("Spectra"); + c25.Divide(2,1); + c25.cd(1); + gPad->SetLogy(); + gPad->SetGrid(); + hSum.DrawCopy("hist"); + funcSum.DrawCopy("SAME"); + + c25.cd(2); + gPad->SetLogy(); + gPad->SetGrid(); + hSumScale.DrawCopy("hist"); + funcScaled.DrawCopy("SAME"); + funcScaledBL.DrawCopy("SAME"); + + + /* + TCanvas &c3 = d->AddTab("Separation"); + gPad->SetLogy(); + hSep.DrawCopy(); + */ + + d->SaveAs("/home_nfs/isdc/jbbuss/singlepe.root"); + 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 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 k = par[5]; - Double_t xTalk = 1; - - Double_t y = 0; - for (int order = 1; order < 14; order++) - { - Double_t arg = (x - order*gain - shift)/sigma; - - Double_t gauss = TMath::Exp(-0.5 * arg * arg/order); - - y += xTalk*gauss; + Double_t xTalk = 1; + + Double_t y = 0; + + for (int order = 1; order < 14; order++) + { + Double_t arg = (x - order*gain - shift)/sigma; + + Double_t gauss = TMath::Exp(-0.5 * arg * arg/order); + + y += xTalk*gauss; // for (int j = 1; j < k; j++) // { - xTalk *= cross; + xTalk *= cross; // } } - return ampl*y; + return ampl*y; } Double_t fcn_cross(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 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 powOfX = par[5]; - Double_t xTalk = 1; - - Double_t y = 0; - for (int order = 1; order < 14; order++) - { - Double_t arg = (x - order*gain - shift)/sigma; - - Double_t gauss = TMath::Exp(-0.5 * arg * arg/order); - - y += xTalk*gauss; + Double_t xTalk = 1; + + Double_t y = 0; + + for (int order = 1; order < 14; order++) + { + Double_t arg = (x - order*gain - shift)/sigma; + + Double_t gauss = TMath::Exp(-0.5 * arg * arg/order); + + y += xTalk*gauss; // xTalk = pow(cross, order*powOfX); @@ -1215,8 +1300,8 @@ // } // cross *= TMath::Exp(-powOfX*cross); - xTalk *= cross; - } - - return ampl*y; + xTalk *= cross; + } + + return ampl*y; } @@ -1224,55 +1309,59 @@ void FitGaussiansToSpectrum( - UInt_t maxOrder, - TH1 &hSource, - TH1 &hDest, - Int_t maxbin, - double fwhm, - Double_t gain - ) + UInt_t maxOrder, + TH1 &hSource, + TH1 &hDest, + Int_t maxbin, + double fwhm, + Double_t gain +) { - Double_t peakposition = hSource.GetBinCenter(maxbin); - char fname[20]; - - // fit gauss functions to spectrum - for (UInt_t nr = 1; nrGetXaxis()->GetNbins(); - Double_t *x = new Double_t[nbins]; - Double_t *y = new Double_t[nbins]; - for (Int_t i=0;iGetXaxis()->GetNbins(); + Double_t *x = new Double_t[nbins]; + Double_t *y = new Double_t[nbins]; + + for (Int_t i=0; iGetXaxis()->GetBinCenter(i+1); y[i] = inputHisto->GetBinContent(i+1); - } - Double_t median = TMath::Median(nbins,x,y); - delete [] x; - delete [] y; - return median; + } + + Double_t median = TMath::Median(nbins,x,y); + delete [] x; + delete [] y; + return median; } // end of PlotMedianEachSliceOfPulse