yangyong/combineICEE/effSigma.C

//******************************************//
// author:
//
//  This function evaluates the effective 
//  sigma of a histogram. 
//
//
//  TODO list:
//
//
//
//  Modified by:
//
//******************************************//


//Double_t effSigma(TH1 * hist)

void effSigma(TH1 * hist,double res[])

{

  TAxis *xaxis = hist->GetXaxis();
  Int_t nb = xaxis->GetNbins();
  if(nb < 10) {
    cout << "effsigma: Not a valid histo. nbins = " << nb << endl;
    return ; 
  }
  
  Double_t bwid = xaxis->GetBinWidth(1);
  if(bwid == 0) {
    cout << "effsigma: Not a valid histo. bwid = " << bwid << endl;
    return ; 
  }
  Double_t xmax = xaxis->GetXmax();
  Double_t xmin = xaxis->GetXmin();
  Double_t ave = hist->GetMean();
  Double_t rms = hist->GetRMS();

  Double_t total=0.;
  for(Int_t i=0; i<nb+2; i++) {
    total+=hist->GetBinContent(i);
  }
  //  if(total < 100.) {
  if(total < 50.) {
    cout << "effsigma: Too few entries " << total << endl;
    return ; 
  }
  Int_t ierr=0;
  Int_t ismin=999;
  
  Double_t rlim=0.683*total;
  Int_t nrms=int(rms/(bwid));    // Set scan size to +/- rms
  if(nrms > nb/10) nrms=nb/10; // Could be tuned...

  Double_t widmin=9999999.;
  
  // cout<<"iscan: "<< xmin <<" "<< xmax <<" "<< ave <<" "<< rms <<" "<< rlim <<" "<< nrms <<endl; 
  
  int binlow = 1; 
  int binhigh = nb; 
  
  double mean_min = ave; 
  
  int binlow2 = 1; 
  int binhigh2 = nb; 
  
  
  for(Int_t iscan=-nrms;iscan<nrms+1;iscan++) { // Scan window centre
    Int_t ibm=(ave-xmin)/bwid+1+iscan;
    Double_t x=(ibm-0.5)*bwid+xmin;
    Double_t xj=x;
    Double_t xk=x;
    Int_t jbm=ibm;
    Int_t kbm=ibm;
    Double_t bin=hist->GetBinContent(ibm);
    total=bin;

    //  cout<<"ibm: "<< iscan <<" "<< ibm <<" "<< xj<<" "<< xk <<" "<<jbm <<" "<<kbm <<endl; 
    
    double mean = bin * x; 
    
    
    for(Int_t j=1;j<nb;j++){
      if(jbm < nb) {
        jbm++;
        xj+=bwid;
        bin=hist->GetBinContent(jbm);
        total+=bin;
        
        if(total > rlim) break;
      }
      else ierr=1;
      if(kbm > 0) {
        kbm--;
        xk-=bwid;
        bin=hist->GetBinContent(kbm);
        total+=bin;
        
        if(total > rlim) break;
      }
      else ierr=1;
    }
    Double_t dxf=(total-rlim)*bwid/bin;
    Double_t wid=(xj-xk+bwid-dxf)*0.5;

    
    //cout<<"dxf: "<< total <<" "<< rlim <<" "<< bwid <<" "<< bin <<" "<< wid <<" "<< xj<<" "<<xk<<" "<< dxf <<endl; 

    if(wid < widmin) {
      widmin=wid;
      ismin=iscan;
      
      binlow = kbm; 
      binhigh = jbm; 


      binlow2 = ibm - 2*(ibm - kbm);  
      binhigh2 = ibm + 2*(-ibm + jbm);  
      
    }
    
    // cout<<"widmin: "<< wid<<" "<< widmin <<" "<< ismin <<" "<<ibm <<" "<< jbm<<" "<< kbm <<" "<<x <<" "<< xj <<" "<< xk <<" "<<(ibm-0.5)*bwid+xmin<<" "<<
    //" "<< (jbm-0.5)*bwid+xmin <<" "<<(kbm-0.5)*bwid+xmin <<endl; 
    
  }
  if(ismin == nrms || ismin == -nrms) ierr=3;
  if(ierr != 0) cout << "effsigma: Error of type " << ierr << endl;
  
  //return widmin;
  
  double mean = 0; 
  double n =0; 
  for(int b = binlow ; b<= binhigh; b++){
    
    if( b<1 || b> nb) continue; 

    double x= (b-0.5)*bwid+xmin;
    mean += hist->GetBinContent(b) * x;
    n += hist->GetBinContent(b); 
  }
  
  mean /= n; 
  
  double mean2 = 0; 
  double n2 =0; 
  for(int b = binlow2 ; b<= binhigh2; b++){

    if( b<1 || b> nb) continue; 
    
    double x= (b-0.5)*bwid+xmin;
    mean2 += hist->GetBinContent(b) * x;
    n2 += hist->GetBinContent(b); 
  }
  
  mean2 /= n2; 
    

  //cout<<widmin <<" "<<mean <<endl; 
  
  res[0] = mean; 
  res[1] = mean2; 
  
  res[2] = widmin; 
  
  
}
Revision:	1.1
Committed:	Tue Apr 10 20:03:11 2012 UTC (13 years, 1 month ago) by yangyong
Content type:	text/plain
Branch:	MAIN
CVS Tags:	V2012comb, pi0etaee_laser20111122_eefloatalpha, HEAD
Log Message:	LaserTag: EcalLaserAPDPNRatios_data_20111122_158851_180363 and alpha Tag ( for Endcap Only) EcalLaserAlphas_lto420-620_progr_data_20111122
#	User	Rev	Content
1	yangyong	1.1	//******************************************//
2			// author:
3			//
4			// This function evaluates the effective
5			// sigma of a histogram.
6			//
7			//
8			// TODO list:
9			//
10			//
11			//
12			// Modified by:
13			//
14			//******************************************//
15
16
17
18
19			//Double_t effSigma(TH1 * hist)
20
21			void effSigma(TH1 * hist,double res[])
22
23			{
24
25			TAxis *xaxis = hist->GetXaxis();
26			Int_t nb = xaxis->GetNbins();
27			if(nb < 10) {
28			cout << "effsigma: Not a valid histo. nbins = " << nb << endl;
29			return ;
30			}
31
32			Double_t bwid = xaxis->GetBinWidth(1);
33			if(bwid == 0) {
34			cout << "effsigma: Not a valid histo. bwid = " << bwid << endl;
35			return ;
36			}
37			Double_t xmax = xaxis->GetXmax();
38			Double_t xmin = xaxis->GetXmin();
39			Double_t ave = hist->GetMean();
40			Double_t rms = hist->GetRMS();
41
42			Double_t total=0.;
43			for(Int_t i=0; i<nb+2; i++) {
44			total+=hist->GetBinContent(i);
45			}
46			// if(total < 100.) {
47			if(total < 50.) {
48			cout << "effsigma: Too few entries " << total << endl;
49			return ;
50			}
51			Int_t ierr=0;
52			Int_t ismin=999;
53
54			Double_t rlim=0.683*total;
55			Int_t nrms=int(rms/(bwid)); // Set scan size to +/- rms
56			if(nrms > nb/10) nrms=nb/10; // Could be tuned...
57
58			Double_t widmin=9999999.;
59
60			// cout<<"iscan: "<< xmin <<" "<< xmax <<" "<< ave <<" "<< rms <<" "<< rlim <<" "<< nrms <<endl;
61
62			int binlow = 1;
63			int binhigh = nb;
64
65			double mean_min = ave;
66
67			int binlow2 = 1;
68			int binhigh2 = nb;
69
70
71
72			for(Int_t iscan=-nrms;iscan<nrms+1;iscan++) { // Scan window centre
73			Int_t ibm=(ave-xmin)/bwid+1+iscan;
74			Double_t x=(ibm-0.5)*bwid+xmin;
75			Double_t xj=x;
76			Double_t xk=x;
77			Int_t jbm=ibm;
78			Int_t kbm=ibm;
79			Double_t bin=hist->GetBinContent(ibm);
80			total=bin;
81
82			// cout<<"ibm: "<< iscan <<" "<< ibm <<" "<< xj<<" "<< xk <<" "<<jbm <<" "<<kbm <<endl;
83
84			double mean = bin * x;
85
86
87			for(Int_t j=1;j<nb;j++){
88			if(jbm < nb) {
89			jbm++;
90			xj+=bwid;
91			bin=hist->GetBinContent(jbm);
92			total+=bin;
93
94			if(total > rlim) break;
95			}
96			else ierr=1;
97			if(kbm > 0) {
98			kbm--;
99			xk-=bwid;
100			bin=hist->GetBinContent(kbm);
101			total+=bin;
102
103			if(total > rlim) break;
104			}
105			else ierr=1;
106			}
107			Double_t dxf=(total-rlim)*bwid/bin;
108			Double_t wid=(xj-xk+bwid-dxf)*0.5;
109
110
111			//cout<<"dxf: "<< total <<" "<< rlim <<" "<< bwid <<" "<< bin <<" "<< wid <<" "<< xj<<" "<<xk<<" "<< dxf <<endl;
112
113			if(wid < widmin) {
114			widmin=wid;
115			ismin=iscan;
116
117			binlow = kbm;
118			binhigh = jbm;
119
120
121			binlow2 = ibm - 2*(ibm - kbm);
122			binhigh2 = ibm + 2*(-ibm + jbm);
123
124			}
125
126			// cout<<"widmin: "<< wid<<" "<< widmin <<" "<< ismin <<" "<<ibm <<" "<< jbm<<" "<< kbm <<" "<<x <<" "<< xj <<" "<< xk <<" "<<(ibm-0.5)*bwid+xmin<<" "<<
127			//" "<< (jbm-0.5)bwid+xmin <<" "<<(kbm-0.5)bwid+xmin <<endl;
128
129			}
130			if(ismin == nrms \|\| ismin == -nrms) ierr=3;
131			if(ierr != 0) cout << "effsigma: Error of type " << ierr << endl;
132
133			//return widmin;
134
135			double mean = 0;
136			double n =0;
137			for(int b = binlow ; b<= binhigh; b++){
138
139			if( b<1 \|\| b> nb) continue;
140
141			double x= (b-0.5)*bwid+xmin;
142			mean += hist->GetBinContent(b) * x;
143			n += hist->GetBinContent(b);
144			}
145
146			mean /= n;
147
148			double mean2 = 0;
149			double n2 =0;
150			for(int b = binlow2 ; b<= binhigh2; b++){
151
152			if( b<1 \|\| b> nb) continue;
153
154			double x= (b-0.5)*bwid+xmin;
155			mean2 += hist->GetBinContent(b) * x;
156			n2 += hist->GetBinContent(b);
157			}
158
159			mean2 /= n2;
160
161
162			//cout<<widmin <<" "<<mean <<endl;
163
164			res[0] = mean;
165			res[1] = mean2;
166
167			res[2] = widmin;
168
169
170			}