Index: fact/tools/rootmacros/PulseTemplates/pulse.C =================================================================== --- fact/tools/rootmacros/PulseTemplates/pulse.C (revision 14475) +++ fact/tools/rootmacros/PulseTemplates/pulse.C (revision 14476) @@ -10,50 +10,76 @@ Pulse::Pulse() { + InitMembers(); +} + +Pulse::Pulse(TString name) +{ + InitMembers(); + mName = name; +} + +Pulse::Pulse(TString name, TH1* histo) +{ + InitMembers(); + mName = name; + mHisto = histo; + Fit(mName, mOptions, mType); +} + +Pulse::Pulse(TString name, TH1* histo, int type) +{ + InitMembers(); + mName = name; + mHisto = histo; + mType = type; + Fit(mName, mOptions, mType); +} + +Pulse::Pulse(TString name, TH1* histo, TString options) +{ + InitMembers(); + mName = name; + mHisto = histo; + mOptions += options; + Fit(mName, mOptions, mType); +} + +Pulse::Pulse(TString name, TH1* histo, TString options, int type) +{ + InitMembers(); + mName = name; + mHisto = histo; + mOptions += options; + mType = type; + Fit(mName, mOptions, mType); +} + + +Pulse::~Pulse() +{ + delete mModel; +} + +void +Pulse::InitMembers() +{ + mName = "pulse"; mHisto = NULL; - mOptions = ""; + mModel = NULL; + mOptions = "S"; mBsl = 0; mHeight = 0; - mStart = 0; - mRising = 0; - mTauRising = 0; - mTauFalling = 0; + mT0 = 0; + mT1 = 0; + mTau1 = 0; + mTau2 = 0; mIntegral = 0; mAmplitude = 0; mPhE = 0; -} - -Pulse::Pulse(TH1* histo) -{ - mHisto = histo; - mOptions = ""; - mBsl = 0; - mHeight = 0; - mStart = 0; - mRising = 0; - mTauRising = 0; - mTauFalling = 0; - mIntegral = 0; - mAmplitude = 0; - mPhE = 0; -} - -Pulse::Pulse(TH1* histo, TString options) -{ - mHisto = histo; - mOptions = options; - mBsl = 0; - mHeight = 0; - mStart = 0; - mRising = 0; - mTauRising = 0; - mTauFalling = 0; - mIntegral = 0; - mAmplitude = 0; - mPhE = 0; -} - -Pulse::~Pulse() -{ - + mType = 0; + mFitProb = 0; + mFitNCalls = 0; + mFitNdf = 0; + mChi2 = 0; } @@ -61,5 +87,6 @@ Pulse::Fit( TString fitName, - TString fitOptions + TString fitOptions, + int type ) { @@ -70,11 +97,45 @@ fitName, fitOptions, + type, fitMin, fitMax - ); + ); } void Pulse::Fit( + TString fitName, + TString fitOptions, + int type, + int fitMin, + int fitMax + ) +{ + fitOptions += "R"; + if (type == 0) + FitSectionWise( + "Model0", + fitOptions, + fitMin, + fitMax + ); + else if (type == 1) + FitContious( + "Model1", + fitOptions, + fitMin, + fitMax + ); + else{ + cout << "not a correct type number --> fitting skipped!" << endl; + return; + } + + CalculateParameters(); + return; +} + +void +Pulse::FitSectionWise( TString fitName, TString fitOptions, @@ -83,5 +144,7 @@ ) { - TF1 fit(fitName, shapeFunc, fitMin, fitMax, 6 ); + mFitMin = fitMin; + mFitMax = fitMax; + mModel = new TF1(fitName, shapeFunc, fitMin, fitMax, 6 ); // ====================================================================== @@ -92,13 +155,12 @@ int first_bin = mHisto->GetXaxis()->GetFirst(); - for (int i = 0; i < 10; i++){ + for (int i = 0; i < 30; i++){ bsl += mHisto->GetBinContent(first_bin+0); } - bsl = bsl/10; + bsl = bsl/30; double stop = mHisto->GetMaximumBin(); //pos of max double height = mHisto->GetBinContent(stop); double start = stop-10; //pos 10 slices before maximum - double rising = (stop - start)/tau; // ====================================================================== @@ -107,36 +169,20 @@ - fit.SetParameters(para); - fit.SetParNames("BSL", "Height", "Start", "Rising", "Tau1", "TauFalling"); - fit.SetLineColor(kRed); - mHisto->Fit(&fit, fitOptions); - - mBsl = fit.GetParameter(0); - mHeight = fit.GetParameter(1); - mStart = fit.GetParameter(2); - mRising = fit.GetParameter(3); - mTauRising = fit.GetParameter(4); - mTauFalling = fit.GetParameter(5); -} - -void -Pulse::Fit2( - TString fitName, - TString fitOptions - ) -{ - int fitMin = mHisto->GetXaxis()->GetFirst(); - int fitMax = mHisto->GetXaxis()->GetLast(); - - Fit2( - fitName, - fitOptions, - fitMin, - fitMax - ); -} - -void -Pulse::Fit2( + mModel->SetParameters(para); + mModel->SetParNames("BSL", "A0", "t0", "t1", "Tau1", "Tau2"); + mModel->SetLineColor(kRed); + + mFitResultPtr = mHisto->Fit(mModel, fitOptions); + + mBsl = mModel->GetParameter(0); + mHeight = mModel->GetParameter(1); + mT0 = mModel->GetParameter(2); + mT1 = mModel->GetParameter(3); + mTau1 = mModel->GetParameter(4); + mTau2 = mModel->GetParameter(5); +} + +void +Pulse::FitContious( TString fitName, TString fitOptions, @@ -145,5 +191,7 @@ ) { - TF1 fit(fitName, shapeFunc2, fitMin, fitMax, 5 ); + mFitMin = fitMin; + mFitMax = fitMax; + mModel = new TF1(fitName, shapeFunc2, fitMin, fitMax, 5 ); // ====================================================================== @@ -154,8 +202,8 @@ int first_bin = mHisto->GetXaxis()->GetFirst(); - for (int i = 0; i < 10; i++){ + for (int i = 0; i < 30; i++){ bsl += mHisto->GetBinContent(first_bin+0); } - bsl = bsl/10; + bsl = bsl/30; double stop = mHisto->GetMaximumBin(); //pos of max @@ -166,31 +214,44 @@ double para[] = {bsl, amplitude, start, tau, tau}; - fit.SetParameters(para); - fit.SetParNames("BSL", "Amplitude", "Start", "Tau1", "TauFalling"); - fit.SetLineColor(kBlue); - mHisto->Fit(&fit, fitOptions); - - mBsl = fit.GetParameter(0); - mAmplitude = fit.GetParameter(1); - mStart = fit.GetParameter(2); - mTauRising = fit.GetParameter(3); - mTauFalling = fit.GetParameter(4); + mModel->SetParameters(para); + mModel->SetParNames("BSL", "A0", "t0", "Tau1", "Tau2"); + mModel->SetLineColor(kBlue); + + mFitResultPtr = mHisto->Fit(mModel, fitOptions); + + mBsl = mModel->GetParameter(0); + mHeight = mModel->GetParameter(1); + mT0 = mModel->GetParameter(2); + mTau1 = mModel->GetParameter(3); + mTau2 = mModel->GetParameter(4); } double Pulse::GetBsl(){ return mBsl;} double Pulse::GetHeight(){ return mHeight;} -double Pulse::GetAvalancheStart(){ return mStart;} -double Pulse::GetAvalancheEnd(){ return mRising;} -double Pulse::GetTimeConstRising(){ return mTauRising;} -double Pulse::GetTimeConstFalling(){ return mTauFalling;} +double Pulse::GetAvalancheStart(){ return mT0;} +double Pulse::GetAvalancheEnd(){ return mT1;} +double Pulse::GetTimeConstRising(){ return mTau1;} +double Pulse::GetTimeConstFalling(){ return mTau2;} double Pulse::GetIntegral(){ return mIntegral;} double Pulse::GetAmplitude(){ return mAmplitude;} int Pulse::GetPE(){ return mPhE;} +void +Pulse::CalculateParameters() +{ + mIntegral = mModel->Integral(mFitMin, mFitMax); + mAmplitude = mModel->GetMaximum() - mBsl; + mFitProb = mFitResultPtr->Prob(); + mFitNCalls = mFitResultPtr->NCalls(); + mFitNdf = mFitResultPtr->Ndf(); + mChi2 = mFitResultPtr->Chi2(); +} + + //double int Heaviside(double val) { - if( val < 0) return 0; if( val >= 0) return 1; + else return 0; } @@ -214,17 +275,23 @@ // helper variables // double max = start + rising * tau1; - double e1 = height * (1 - Exp(-(t[0]-start)/tau1 ) ); + double e1 = (1 - Exp(-(t[0]-start)/tau1 ) ); // double e2 = (-1) * height * (1 - Exp(-(x[0]-max)/tau2 ) ); - double e2 = (-1) * height * (1 - Exp(-(t[0]-stop)/tau2 ) ); + double e2 = (-1) * (1 - Exp(-(t[0]-stop)/tau2 ) ); // calculate return value + returnval += Heaviside(t[0]-start)*e1; + returnval += Heaviside(t[0]-stop)*e2; + returnval *= height; returnval += bsl; - - if (t[0] > start) - returnval += e1; - if (t[0] > stop) - returnval += e2; +// if (t[0] > start) +// returnval += e1; +// if (t[0] > stop) +// returnval += e2; return returnval; -} + + +} + + /* @@ -282,3 +349,2 @@ return returnval; } -