diff --git a/fles/mcbm2018/tasks/CbmCheckTiming.cxx b/fles/mcbm2018/tasks/CbmCheckTiming.cxx index 65db603f799eb252c85fd5b8ecc7f513e70abe46..3d857df12738d5837d2e82b78dbbdc0fecd804c1 100644 --- a/fles/mcbm2018/tasks/CbmCheckTiming.cxx +++ b/fles/mcbm2018/tasks/CbmCheckTiming.cxx @@ -16,13 +16,14 @@ #include "CbmTofDigi.h" #include "CbmTrdDigi.h" +#include "FairLogger.h" #include "FairRootManager.h" #include "FairRunOnline.h" -#include <Logger.h> #include "TClonesArray.h" #include "TH1.h" #include "TH2.h" +#include "TF1.h" #include "THttpServer.h" #include <TFile.h> @@ -529,9 +530,10 @@ void CbmCheckTiming::Exec(Option_t* /*option*/) { if (fCheckTimeOrdering) CheckTimeOrder(); if (fCheckInterSystemOffset) CheckInterSystemOffset(); - + if (0 < fNrTs && 0 == fNrTs % fNrTsForFit) LOG(info) << "Fitting peaks for number of TS = " << fNrTs; + if (0 < fNrTs && 0 == fNrTs % fNrTsForFit) FitPeaks(); //Added AR if (0 < fNrTs && 0 == fNrTs % 2000) WriteHistos(); - + fNrTs++; } @@ -1076,6 +1078,165 @@ void CbmCheckTiming::Finish() { LOG(info) << Form("Checked %6d Timeslices", fNrTs); } + +void CbmCheckTiming::FitPeaks() { + +// Peak positions for differences wrt T0 +trd_peak_pos = fT0TrdDiff->GetMaximumBin()*fT0TrdDiff->GetBinWidth(1)+fT0TrdDiff->GetXaxis()->GetXmin(); +sts_peak_pos = fT0StsDiff->GetMaximumBin()*fT0StsDiff->GetBinWidth(1)+fT0StsDiff->GetXaxis()->GetXmin(); +much_peak_pos = fT0MuchDiff->GetMaximumBin()*fT0MuchDiff->GetBinWidth(1)+fT0MuchDiff->GetXaxis()->GetXmin(); +tof_peak_pos = fT0TofDiff->GetMaximumBin()*fT0TofDiff->GetBinWidth(1)+fT0TofDiff->GetXaxis()->GetXmin(); +rich_peak_pos = fT0RichDiff->GetMaximumBin()*fT0RichDiff->GetBinWidth(1)+fT0RichDiff->GetXaxis()->GetXmin(); +psd_peak_pos = fT0PsdDiff->GetMaximumBin()*fT0PsdDiff->GetBinWidth(1)+fT0PsdDiff->GetXaxis()->GetXmin(); + +// Peak positions for differences wrt T0 if coincidence with TOF +trd_coin_peak_pos = fSelT0TrdDiff->GetMaximumBin()*fSelT0TrdDiff->GetBinWidth(1)+fSelT0TrdDiff->GetXaxis()->GetXmin(); +sts_coin_peak_pos = fSelT0StsDiff->GetMaximumBin()*fSelT0StsDiff->GetBinWidth(1)+fSelT0StsDiff->GetXaxis()->GetXmin(); +much_coin_peak_pos = fSelT0MuchDiff->GetMaximumBin()*fSelT0MuchDiff->GetBinWidth(1)+fSelT0MuchDiff->GetXaxis()->GetXmin(); +//tof_coin_peak_pos = fT0TofDiff->GetMaximumBin()*fT0TofDiff->GetBinWidth(1)+fT0TofDiff->GetXaxis()->GetXmin(); +rich_coin_peak_pos = fSelT0RichDiff->GetMaximumBin()*fSelT0RichDiff->GetBinWidth(1)+fSelT0RichDiff->GetXaxis()->GetXmin(); +psd_coin_peak_pos = fSelT0PsdDiff->GetMaximumBin()*fSelT0PsdDiff->GetBinWidth(1)+fSelT0PsdDiff->GetXaxis()->GetXmin(); + +LOG(info) << "STS entries = " << fT0StsDiff->GetEntries(); +LOG(info) << "STS-T0 entries if T0 in coincidence with TOF ... " << fSelT0StsDiff->GetEntries(); +//LOG(info) << " "; +LOG(info) << "MUCH entries = " << fT0MuchDiff->GetEntries(); +LOG(info) << "MUCH-T0 entries if T0 in coincidence with TOF ... " << fSelT0MuchDiff->GetEntries(); +//LOG(info) << " "; +LOG(info) << "TRD entries = " << fT0TrdDiff->GetEntries(); +LOG(info) << "TRD-T0 entries if T0 in coincidence with TOF ... " << fSelT0TrdDiff->GetEntries(); +//LOG(info) << " "; +LOG(info) << "TOF entries = " << fT0TofDiff->GetEntries(); +//LOG(info) << " "; +LOG(info) << "RICH entries = " << fT0RichDiff->GetEntries(); +LOG(info) << "RICH-T0 entries if T0 in coincidence with TOF ... " << fSelT0RichDiff->GetEntries(); +//LOG(info) << " "; +LOG(info) << "PSD entries = " << fT0PsdDiff->GetEntries(); +LOG(info) << "PSD-T0 entries if T0 in coincidence with TOF ... " << fSelT0PsdDiff->GetEntries(); +//LOG(info) << " "; +LOG(info) << " "; +LOG(info) << "STS peak position [ns] = "<< sts_peak_pos; +LOG(info) << "STS peak position [ns] if T0 in coincidence with TOF ... "<< sts_coin_peak_pos; +//LOG(info) << " "; +LOG(info) << "MUCH peak position [ns] = "<< much_peak_pos; +LOG(info) << "MUCH peak position [ns] if T0 in coincidence with TOF ... "<< much_coin_peak_pos; +//LOG(info) << " "; +LOG(info) << "TRD peak position [ns] = "<< trd_peak_pos; +LOG(info) << "TRD peak position [ns] if T0 in coincidence with TOF ... "<< trd_coin_peak_pos; +//LOG(info) << " "; +LOG(info) << "TOF peak position [ns] = "<< tof_peak_pos; +//LOG(info) << " "; +LOG(info) << "RICH peak position [ns] = "<< rich_peak_pos; +LOG(info) << "RICH peak position [ns] if T0 in coincidence with TOF ... "<< rich_coin_peak_pos; +//LOG(info) << " "; +LOG(info) << "PSD peak position [ns] = "<< psd_peak_pos; +LOG(info) << "PSD peak position [ns] if T0 in coincidence with TOF ... "<< psd_coin_peak_pos; +LOG(info) << " "; + +//Average height of bins... +trd_average = fT0TrdDiff->Integral()/(fT0TrdDiff->GetNbinsX()); +sts_average = fT0StsDiff->Integral()/(fT0StsDiff->GetNbinsX()); //works +much_average = fT0MuchDiff->Integral()/(fT0MuchDiff->GetNbinsX()); +tof_average = fT0TofDiff->Integral()/(fT0TofDiff->GetNbinsX()); +rich_average = fT0RichDiff->Integral()/(fT0RichDiff->GetNbinsX()); +psd_average = fT0PsdDiff->Integral()/(fT0PsdDiff->GetNbinsX()); + +//LOG(info) << "STS average peak height = "<< sts_average; + +//Typical width in ns for Gauss +trd_width0_ns = 120; +sts_width0_ns = 30; +much_width0_ns = 100; +tof_width0_ns = 20; +rich_width0_ns = 40; +psd_width0_ns = 20; + +//TRD +if(trd_average > 0){ +//TCanvas* canvas_trd = new TCanvas ("canvas_trd"); + TF1 *gs_trd = new TF1("gs_trd", "gaus(0)+pol0(3)", trd_peak_pos-2*trd_width0_ns, trd_peak_pos+2*trd_width0_ns); + gs_trd->SetParameters(0.7*trd_average, trd_peak_pos, trd_width0_ns, trd_average); + fT0TrdDiff->Fit("gs_trd","R"); + TF1 *fitresult_trd = fT0TrdDiff->GetFunction("gs_trd"); + LOG(info) << "TRD parameters from Gauss fit = " << + fitresult_trd->GetParameter(0) << ", " << + fitresult_trd->GetParameter(1) << ", " << + fitresult_trd->GetParameter(2); + LOG(info)<< "TRD signal/background (p0/p3) " << (fitresult_trd->GetParameter(0))/(fitresult_trd->GetParameter(3)); +//LOG(info)<< "TRD chi2 per DOF= " << (fitresult_trd->GetChisquare())/(fitresult_trd->GetNDF()); +// fT0TrdDiff->Draw(); +} + +//STS +if(sts_average > 0){ + TF1 *gs_sts = new TF1("gs_sts", "gaus(0)+pol0(3)", sts_peak_pos-2*sts_width0_ns, sts_peak_pos+2*sts_width0_ns); + gs_sts->SetParameters(sts_average, sts_peak_pos, sts_width0_ns, sts_average); + fT0StsDiff->Fit("gs_sts","R"); + TF1 *fitresult_sts = fT0StsDiff->GetFunction("gs_sts"); + LOG(info) << "STS parameters from Gauss fit = " << + fitresult_sts->GetParameter(0) << ", " << + fitresult_sts->GetParameter(1) << ", " << + fitresult_sts->GetParameter(2); + LOG(info)<< "STS signal/background (p0/p3) " << (fitresult_sts->GetParameter(0))/(fitresult_sts->GetParameter(3)); +} + +//MUCH +if(much_average > 0){ + TF1 *gs_much = new TF1("gs_much", "gaus(0)+pol0(3)", much_peak_pos-2*much_width0_ns, much_peak_pos+2*much_width0_ns); + gs_much->SetParameters(much_average, much_peak_pos, much_width0_ns, much_average); + fT0MuchDiff->Fit("gs_much","R"); + TF1 *fitresult_much = fT0MuchDiff->GetFunction("gs_much"); + LOG(info) << "MUCH parameters from Gauss fit = " << + fitresult_much->GetParameter(0) << ", " << + fitresult_much->GetParameter(1) << ", " << + fitresult_much->GetParameter(2); + LOG(info)<< "MUCH signal/background (p0/p3) " << (fitresult_much->GetParameter(0))/(fitresult_much->GetParameter(3)); +} + +//TOF +if(tof_average > 0){ + TF1 *gs_tof = new TF1("gs_tof", "gaus(0)+pol0(3)", tof_peak_pos-2*tof_width0_ns, tof_peak_pos+2*tof_width0_ns); + gs_tof->SetParameters(tof_average, tof_peak_pos, tof_width0_ns, tof_average); + fT0TofDiff->Fit("gs_tof","R"); + TF1 *fitresult_tof = fT0TofDiff->GetFunction("gs_tof"); + LOG(info) << "TOF parameters from Gauss fit = " << + fitresult_tof->GetParameter(0) << ", " << + fitresult_tof->GetParameter(1) << ", " << + fitresult_tof->GetParameter(2); + LOG(info)<< "TOF signal/background (p0/p3) " << (fitresult_tof->GetParameter(0))/(fitresult_tof->GetParameter(3)); +} + +//RICH +if(rich_average > 0){ + TF1 *gs_rich = new TF1("gs_rich", "gaus(0)+pol0(3)", rich_peak_pos-2*rich_width0_ns, rich_peak_pos+2*rich_width0_ns); + gs_rich->SetParameters(0.5*rich_average, rich_peak_pos, rich_width0_ns, rich_average); + fT0RichDiff->Fit("gs_rich","R"); + TF1 *fitresult_rich = fT0RichDiff->GetFunction("gs_rich"); + LOG(info) << "RICH parameters from Gauss fit = " << + fitresult_rich->GetParameter(0) << ", " << + fitresult_rich->GetParameter(1) << ", " << + fitresult_rich->GetParameter(2); + LOG(info)<< "RICH signal/background (p0/p3) " << (fitresult_rich->GetParameter(0))/(fitresult_rich->GetParameter(3)); +} + +//PSD +if(psd_average > 0){ + TF1 *gs_psd = new TF1("gs_psd", "gaus(0)+pol0(3)", psd_peak_pos-2*psd_width0_ns, psd_peak_pos+2*psd_width0_ns); + gs_psd->SetParameters(psd_average, psd_peak_pos, psd_width0_ns, psd_average); + fT0PsdDiff->Fit("gs_psd","R"); + TF1 *fitresult_psd = fT0PsdDiff->GetFunction("gs_psd"); + LOG(info) << "PSD parameters from Gauss fit = " << + fitresult_psd->GetParameter(0) << ", " << + fitresult_psd->GetParameter(1) << ", " << + fitresult_psd->GetParameter(2); + LOG(info)<< "PSD signal/background (p0/p3) " << (fitresult_psd->GetParameter(0))/(fitresult_psd->GetParameter(3)); +} + + +} //end FitPeaks() + + + void CbmCheckTiming::WriteHistos() { TFile* old = gFile; TFile* outfile = TFile::Open(fOutFileName, "RECREATE"); diff --git a/fles/mcbm2018/tasks/CbmCheckTiming.h b/fles/mcbm2018/tasks/CbmCheckTiming.h index aa3024d5cca821585fccebf1e2d1a220e2b3aca3..4ef25a36954ad656590052700fec041cf1204563 100644 --- a/fles/mcbm2018/tasks/CbmCheckTiming.h +++ b/fles/mcbm2018/tasks/CbmCheckTiming.h @@ -72,7 +72,7 @@ public: void SetRichOffsetSearchRange(Int_t val = 1000) { fRichOffsetRange = val; } void SetPsdOffsetSearchRange(Int_t val = 1000) { fPsdOffsetRange = val; } - + void SetNrTsForFit(Int_t val = 200) { fNrTsForFit = val; } inline void SetT0PulserTotLimits(UInt_t uMin, UInt_t uMax) { fuMinTotPulserT0 = uMin; fuMaxTotPulserT0 = uMax; @@ -126,7 +126,7 @@ private: Int_t CalcNrBins(Int_t); void CreateHistos(); void WriteHistos(); - + void FitPeaks(); /** Input array from previous already existing data level **/ CbmDigiManager* fDigiMan = nullptr; //! @@ -135,6 +135,39 @@ private: const std::vector<CbmTofDigi>* fT0DigiVec = nullptr; //! TClonesArray* fT0DigiArr = nullptr; //! + + /** Peak position of time difference histograms **/ + Double_t trd_peak_pos; + Double_t sts_peak_pos; + Double_t much_peak_pos; + Double_t tof_peak_pos; + Double_t rich_peak_pos; + Double_t psd_peak_pos; + + /** Peak position of time difference histograms with TOF coincidence**/ + Double_t trd_coin_peak_pos; + Double_t sts_coin_peak_pos; + Double_t much_coin_peak_pos; + Double_t rich_coin_peak_pos; + Double_t psd_coin_peak_pos; + + /** Average height of bins **/ + Double_t trd_average; + Double_t sts_average; + Double_t much_average; + Double_t tof_average; + Double_t rich_average; + Double_t psd_average; + + /** Typical width from Gauss **/ + Double_t trd_width0_ns; + Double_t sts_width0_ns; + Double_t much_width0_ns; + Double_t tof_width0_ns; + Double_t rich_width0_ns; + Double_t psd_width0_ns; + + /// Variables to store the previous digi time Double_t fPrevTimeT0 = 0.; Double_t fPrevTimeSts = 0.; @@ -187,7 +220,7 @@ private: Int_t fTofOffsetRange = 1000; Int_t fRichOffsetRange = 1000; Int_t fPsdOffsetRange = 1000; - + Int_t fNrTsForFit = 100; Int_t fBinWidth = 1; TH1* fT0StsDiff = nullptr;