In This Package:

pmt_full_model_fit.h File Reference

#include <stdlib.h>
#include <iostream>
#include <TROOT.h>
#include <TApplication.h>
#include <TVirtualX.h>
#include <TGFileDialog.h>
#include <TPostScript.h>
#include <TStyle.h>
#include <TSystem.h>
#include <TRootHelpDialog.h>
#include <TMinuit.h>
#include <TVirtualFitter.h>
#include "TPaveStats.h"
#include "TText.h"
#include <fstream>
#include <iomanip>
#include <string.h>
#include <math.h>
#include <TMath.h>
#include "TF1.h"
#include "TCanvas.h"
#include <TH1F.h>
#include "Event/ReadoutPmtCrate.h"
#include "StatisticsSvc/IStatisticsSvc.h"
#include "DataSvc/ICableSvc.h"
#include "DataSvc/ICalibDataSvc.h"
#include "Context/ServiceMode.h"
#include "pmt_full_model.h"

Include dependency graph for pmt_full_model_fit.h:

Go to the source code of this file.

Defines
#define	FULL_PMT_RESPONSE_MODEL
#define	MAX_PE_NUMBER   20
Functions
Double_t	fitf (Double_t *x, Double_t *par)
int	pmt_full_model_fit (TH1F *histo1, double *par, double *par_err)
	--**---* main function *---**----
Variables
Double_t	x [4096]
Double_t	y [4096]
Float_t	num [4096]
Float_t	amplitude [4096]
Double_t	mu
Double_t	Q0
Double_t	Delta0
Double_t	Q1
Double_t	Delta1
Double_t	w
Double_t	A
Int_t	NON_ZERO_MAX_NUM
Double_t	EVENT_TOTLE_NUM
Double_t	N
Double_t	xlow
Double_t	xhigh
long	NNN

---

Define Documentation

#define FULL_PMT_RESPONSE_MODEL

Definition at line 54 of file pmt_full_model_fit.h.

#define MAX_PE_NUMBER   20

Definition at line 57 of file pmt_full_model_fit.h.

---

Function Documentation

Double_t fitf	(	Double_t *	x,
		Double_t *	par
	)

Definition at line 98 of file pmt_full_model_fit.h.

{
        Double_t RESULT=0;
        int Nmax=int(N);

        NNN++;

        if(x[0]<xlow||x[0]>xhigh)
                return 0;
        else
        {

                for(Int_t i=0;i<=Nmax;i++)
                {
                        #ifdef FULL_PMT_RESPONSE_MODEL
                        //SER_n
                        RESULT=RESULT+SER_n(i,x,par); // for full model
                        #endif
                }

                RESULT=RESULT*EVENT_TOTLE_NUM;
        }

        if(NNN%50000==0)
                cout<<"in fitf_IV..............."<<RESULT<<"  "<<NNN<<endl;
        return RESULT;
 }

int pmt_full_model_fit	(	TH1F *	histo1,
		double *	par,
		double *	par_err
	)

--**---* main function *---**----

Definition at line 135 of file pmt_full_model_fit.h.

{

        NNN=0;

        EVENT_TOTLE_NUM=par_err[0];
        if(EVENT_TOTLE_NUM<=0)
        {
                cout<<"total event number less than 0"<<endl;
                return -1;
        }

        for(int i=0;i<4096;i++)
        {
                x[i]=0;
                y[i]=0;
                num[i]=0;
                amplitude[i]=0;
        }
        Int_t max_nonzero_bin_X=0;
        Int_t min_nonzero_bin_X=0;

//************************************************************************************
        
        // get the embody fitting value for the channel
        //end to get the embody fitting value for the channel

        for(Int_t i=0;i<4096;i++)
        {     
                amplitude[i]=histo1->GetBinContent(i);
                num[i]=int(histo1->GetBinCenter(i)+1);          
                if(amplitude[i]!=0)
                        max_nonzero_bin_X=int(histo1->GetBinCenter(i)+1);
                x[i]=int(histo1->GetBinCenter(i)+1);
                y[i]=amplitude[i];
        }

//++++++++++++++++++++++++++

        for(int i=0;i<histo1->GetMaximumBin()-1;i++)
        {
                if(amplitude[i]==0)
                {
                        min_nonzero_bin_X=int(histo1->GetBinCenter(i)+1);
                }
                else
                {
                        break;
                }
        }



//---(------------------------------------------------------------------------------------------------------


        if(par[0]<=0.01||par[0]>=5)
        {
                cout<<" average_mu="<<par[0]<<" not suitable"<<endl;
                return -1;
        }
        else if(par[0]<=0.5)
        {
                N=6;
        }
        else if(par[0]<=1)
        {
                N=12;
        }
        else
        {
                N=MAX_PE_NUMBER;
        }


        double high_range=500;

        xlow=min_nonzero_bin_X-1;       // add 50 is for the FEE channel threshold effect
        xhigh=double(max_nonzero_bin_X-1);

        mu=par[0];
        double max_content_bin_X=histo1->GetBinCenter(histo1->GetMaximumBin());
        cout<<"the max_content_bin_X="<<max_content_bin_X<<endl;

        double max_content_X=histo1->GetMaximum();
        double bin_step=50;
        for(int t=0;t<10&&(max_content_bin_X+high_range<=max_nonzero_bin_X);t++)
        {
                if(histo1->GetBinContent(int(max_content_bin_X+high_range))<=0.01*max_content_X)
                {
                        break;
                }
                else
                {
                        high_range=high_range+bin_step;
                }
        }

        if(TMath::Abs(max_content_bin_X-histo1->GetMean())>100)
        {
                cout<<"max_content_bin_X="<<max_content_bin_X<<"        larger than the mean="<<histo1->GetMean()<<" more than 100"<<endl;
                return -1;
        }


        if(xhigh>max_content_bin_X-1+high_range)
                xhigh=max_content_bin_X-1+high_range;

        if(xlow>=max_content_bin_X-1)
        {
                cout<<"--------------+++++++++------------------------"<<endl;
                cout<<"xlow="<<xlow<<" larger than maxBin="<<max_content_bin_X-1<<endl;
                cout<<"--------------+++++++++------------------------"<<endl;
                return -1;
        }

        {
                Q0=par[1];
                if(Q0<0||Q0>xhigh||TMath::IsNaN(Q0))
                        Q0=0;                        /* the pedestal position */
                Delta0=par[2];
                if(Delta0<0||TMath::IsNaN(Delta0))
                        Delta0=Q0*0.02;          /* pedestal deviation */
                Q1=par[3];
                if(Q1<=0||Q1>xhigh||TMath::IsNaN(Q1)||Q1>(max_content_bin_X-1)*2)
                {
                        if(max_content_bin_X-1-Q0>0) //for the real spe position without pedestal
                        {
                                Q1=max_content_bin_X-1-Q0;
                        }
                        else
                        {
                                cout<<"--------------+++++++++------------------------"<<endl;
                                cout<<"max_content_bin_X="<<max_content_bin_X<<" lower than Q0="<<Q0<<endl;
                                cout<<"--------------+++++++++------------------------"<<endl;
                                return -1;
                        }
                }
                Delta1=par[4];
                if(Delta1<=0||TMath::IsNaN(Delta1))
                        Delta1=Q1*0.3;         /* the deviation of Gaussion of one p.e. */
                w=par[5];
                if(w<=0||TMath::IsNaN(w))
                        w=0.02;                       /* exponential part's probability */
                A=par[6];
                if(A<=0||TMath::IsNaN(A))
                        A=2.5*1/(Q1);                        /* exponential part's slope */
        }

        if(xlow<Q0+3*Delta0)
        {
                cout<<"xlow="<<xlow<<" larger than 3sigma("<<Delta0<<") of Q0("<<Q0<<")"<<endl;
                mu=mu-mu*(xlow-Q0)/Delta0*0.05;
                if(mu<0.5*par[0])
                {
                        mu=0.5*par[0];
                }
                cout<<"the corrected mu="<<mu<<"(the raw mu="<<par[0]<<")"<<endl;
        }

        if(xlow>=max_content_bin_X-1)
        {
                cout<<"--------------+++++++++------------------------"<<endl;
                cout<<"xlow="<<xlow<<" larger than maxBin="<<max_content_bin_X-1<<endl;
                cout<<"--------------+++++++++------------------------"<<endl;
                return -1;
        }
        if(xhigh>max_content_bin_X-1+high_range) // for the fitting maximum bin X should be less than the maximum X bin +300 (>5sigma)
        {
                xhigh=max_content_bin_X-1+high_range;
        }

cout<<"the data mu ="<<par[0]<<endl;
cout<<"the p.e. maximum for fitting is :"<<N<<endl;
        cout<<"the fitting range is:("<<xlow<<","<<xhigh<<")"<<endl;
        cout<<"mu="<<mu<<endl;
        cout<<"Q0="<<Q0<<endl;
        cout<<"Delta0="<<Delta0<<endl;
        cout<<"Q1="<<Q1<<endl;
        cout<<"Delta1="<<Delta1<<endl;
        cout<<"w="<<w<<endl;
        cout<<"A="<<A<<endl;

//++++++++++++++++++++++++++

        
                cout<<"=================================fitting=========================="<<endl;

                cout<<"=================================TF1=========================="<<endl;
                        TF1 *myfunc=new TF1("myfunc",fitf,0,4096,7);
                        histo1->SetAxisRange(xlow,xhigh);
                        myfunc->SetRange(xlow,xhigh);
                        myfunc->SetParameters(mu,Q0,Delta0,Q1,Delta1,w,A);



                        myfunc->SetParLimits(1,Q0+0.01,Q0+0.01);
                        myfunc->SetParLimits(2,Delta0,Delta0);
                        
                        histo1->Fit("myfunc","BNR"); //N: no draw; R: range; L: likilyhood
                                                        //B: Disable the automatic computation of the initial parameter values for the standard functions like poln, expo, and gaus.
                                                        //http://minos.phy.bnl.gov/~bviren/cern/root/rug/rug8.html

                        myfunc->GetParameters(par);
                        for(int i=0;i<7;i++)
                        {
                                par_err[i]=myfunc->GetParError(i);
                        }


        cout<<"               -----------------------Completed!-----------------------"<<endl;

        return 0;

}

---

Variable Documentation

Double_t x[4096]

Definition at line 61 of file pmt_full_model_fit.h.

Double_t y[4096]

Definition at line 62 of file pmt_full_model_fit.h.

Float_t num[4096]

Definition at line 63 of file pmt_full_model_fit.h.

Float_t amplitude[4096]

Definition at line 64 of file pmt_full_model_fit.h.

Double_t mu

Definition at line 73 of file pmt_full_model_fit.h.

Double_t Q0

Definition at line 74 of file pmt_full_model_fit.h.

Double_t Delta0

Definition at line 75 of file pmt_full_model_fit.h.

Double_t Q1

Definition at line 76 of file pmt_full_model_fit.h.

Double_t Delta1

Definition at line 77 of file pmt_full_model_fit.h.

Double_t w

Definition at line 78 of file pmt_full_model_fit.h.

Double_t A

Definition at line 79 of file pmt_full_model_fit.h.

Int_t NON_ZERO_MAX_NUM

Definition at line 31 of file pmt_full_model.h.

Double_t EVENT_TOTLE_NUM

Definition at line 32 of file pmt_full_model.h.

Double_t N

Definition at line 34 of file pmt_full_model.h.

Double_t xlow

Definition at line 36 of file pmt_full_model.h.

Double_t xhigh

Definition at line 37 of file pmt_full_model.h.

long NNN

Definition at line 39 of file pmt_full_model.h.

---