theFit.h File Reference

#include "TH1.h"
#include "TF1.h"
#include "TMath.h"

Go to the source code of this file.

Functions

void fit_histo (class TH1F *histo, Double_t *params, Double_t *parerrs, Bool_t drawIt)
void fit_histo (class TH1F *histo)
Double_t myfped_1gauss (Double_t *x, Double_t *par)
Double_t mypefit (Double_t *x, Double_t *par)
Double_t my_gamma (Double_t *x, Double_t *par)

Function Documentation

void fit_histo ( class TH1F *  histo  ) 

Definition at line 11 of file theFit.h.

References fit_histo().

Referenced by fit_histo().

00011                                   {
00012    Double_t params[9]={0};
00013    Double_t parerrs[9]={0};
00014    fit_histo(histo,params,parerrs,kTRUE);
00015 }

void fit_histo ( class TH1F *  histo,
Double_t *  params,
Double_t *  parerrs,
Bool_t  drawIt = kFALSE 
)

Definition at line 21 of file theFit.h.

References myfped_1gauss(), mypefit(), and nentries.

00022 {
00023    
00024    TF1 *peddy = new TF1("peddy",myfped_1gauss,-500,16500,3);
00025    TF1 *fitty = new TF1("fitty",mypefit,-500,16500,8);
00026    fitty->SetNpx(10000);
00027             
00028    Int_t peakbin=histo->GetMaximumBin();
00029    Int_t nentries=(int)histo->GetEntries();
00030    Int_t lastbin=0;
00031    Int_t nbinsx = histo->GetNbinsX();
00032    for(Int_t j = 0 ; j < nbinsx-1 ; j++) {
00033       Stat_t tempbc =histo->GetBinContent(j+1);
00034       if(tempbc > 3 && lastbin<j) {
00035          lastbin=j;
00036       }
00037       
00038    } 
00039 
00040    Int_t startbin= (int)histo->GetBinCenter(0);
00041    Int_t finishbin= (int) histo->GetBinCenter(nbinsx);
00042 
00043    Int_t zerobin=-startbin+1;
00044    
00045    Int_t minbin=peakbin;
00046    Int_t maxbin=16000; 
00047    Bool_t foundped=kTRUE;
00048    Int_t tempval;
00049    for(Int_t i = 0 ; i < peakbin ; i++) {
00050       if((peakbin-i)==1&&i<4) {
00051          minbin=0;
00052          break;
00053       }     
00054       if((peakbin-i)==1&&i>3) {
00055          minbin=i;
00056          break;
00057       }
00058       tempval=(int)histo->GetBinContent(i+1);
00059       if(tempval>0&&(peakbin-i)<10) {
00060          minbin=i;
00061          break;
00062       }
00063       
00064       if(tempval>10) {
00065          minbin=i;
00066          foundped=kFALSE;
00067          break;
00068       }
00069       
00070    }
00071 
00072    if(!foundped) {
00073       //       cout << "Pedestal not found!" << endl;
00074       Int_t tempvali,tempvalh,tempvalg,tempvalj,tempvalk;
00075       
00076 
00077       for(Int_t i = 2 ; i < peakbin ; i++) {
00078          tempvalg=(int)histo->GetBinContent(i-1);
00079          tempvalh=(int)histo->GetBinContent(i);
00080          tempvali=(int)histo->GetBinContent(i+1);
00081          tempvalj=(int)histo->GetBinContent(i+2);
00082          tempvalk=(int)histo->GetBinContent(i+3);
00083          if((tempvali>tempvalh)&&(tempvali>tempvalg)&(tempvali>tempvalj)&&(tempvali>tempvalk)&&tempvali>5) {
00084             peakbin=i;
00085             break;
00086          }
00087       }
00088       for(Int_t i = 0 ; i < peakbin ; i++) {
00089          if((peakbin-i)==1&&i<4) {
00090             minbin=0;
00091             break;
00092          }
00093          if((peakbin-i)==1&&i>3) {
00094             minbin=i;
00095             break;
00096          }
00097          
00098          if(histo->GetBinContent(i)>0&&(peakbin-i)<10) {
00099             minbin=i;
00100             break;
00101          }
00102       }
00103    }
00104    
00105    if(peakbin+(peakbin-minbin)<16000) {
00106       maxbin=peakbin+(peakbin-minbin);
00107       Int_t maxtemp=maxbin;
00108       for(Int_t i=maxbin+1; i < maxbin+10; i++) {
00109          Int_t temp1=(int)histo->GetBinContent(maxtemp);
00110          Int_t temp2=(int)histo->GetBinContent(i);
00111          if(temp2==0) break;
00112          //cout << temp1 << "       "<< temp2 << endl;
00113          if((Float_t(temp1)/Float_t(temp2))>1.3) {
00114             maxtemp=i;
00115          }
00116          else break;
00117       }    
00118       maxbin=maxtemp;
00119    }
00120    
00121    
00122   
00123             
00124    //cout << peakbin << "\t" << minbin << "\t" << maxbin << endl; 
00125    
00126    Int_t plotmin=(int)histo->GetBinCenter(minbin)-20;
00127    Int_t plotmax=(int)histo->GetBinCenter(lastbin)+20;
00128   
00129    if(drawIt)
00130       histo->SetAxisRange(plotmin,plotmax);
00131    
00132    int temppedmean=(int)histo->GetBinCenter(peakbin);
00133 
00134    //   cout << "Pedestal " << histo->GetBinCenter(peakbin) << endl;
00135    peddy->SetParameters(temppedmean,2,nentries);
00136    peddy->SetParNames("Position","width","No. of entries");
00137    peddy->SetParLimits(0,temppedmean-5,temppedmean+5);
00138    peddy->SetParLimits(1,0,10);
00139    peddy->SetParLimits(2,0,nentries);
00140    
00141    peddy->SetRange(temppedmean-10,temppedmean+10);
00142    if(drawIt)
00143        histo->Fit("peddy","NR");
00144    else
00145        histo->Fit("peddy","NRQ");
00146    Int_t ntotped=0;
00147    Double_t pedsig=peddy->GetParameter(1);
00148    //Int_t minbin=int(500-pedsig*3);
00149    //Int_t maxbin=int(500+pedsig*3);
00150    
00151    ntotped=(int)histo->Integral(minbin,maxbin);
00152    Double_t meanpe=0;
00153    Double_t pedpos=peddy->GetParameter(0);
00154    params[0]=peddy->GetParameter(0);
00155    parerrs[0]=peddy->GetParError(0);
00156    params[5]=peddy->GetParameter(1);
00157    parerrs[5]=peddy->GetParError(1);
00158    params[8]=peddy->GetParameter(2);
00159    parerrs[8]=peddy->GetParError(2);
00160    if(ntotped!=0) {
00161       meanpe=TMath::Log(Float_t(nentries)/Float_t(ntotped));
00162       //cout << "meanpe " << meanpe << endl;
00163    }
00164    else meanpe=0;
00165    Double_t spe=60;
00166    if(meanpe>0) {
00167      spe=histo->GetMean()/meanpe;
00168    }
00169    if(spe<40 || spe>90) spe=60;
00170    Double_t spewidth=0.51;
00171    Double_t dynfrac=0.05;
00172    Double_t dynscale=4;
00173    
00174                   
00175    fitty->SetParameters(pedpos,meanpe,spe,spewidth,dynfrac,pedsig,dynscale,nentries);
00176    fitty->SetParNames("pedpos","meanpe","spepos","spewidth","dynodefrac","pedsigma","dynscale","numentries");
00177    fitty->SetParLimits(0,1,1);
00178    fitty->SetParLimits(1,0,10);
00179    fitty->SetParLimits(2,20,120);
00180    fitty->SetParLimits(3,0,1);
00181    fitty->SetParLimits(4,0,0.2);
00182    fitty->SetParLimits(6,2,9);
00183    fitty->SetParLimits(7,1,1);
00184    fitty->SetParLimits(5,1,1);
00185    
00186    fitty->SetRange(histo->GetBinCenter(minbin),histo->GetBinCenter(lastbin));
00187    if(meanpe>0.04) {
00188        if(drawIt)
00189            histo->Fit("fitty","NR");
00190        else
00191            histo->Fit("fitty","QNR");
00192       for(Int_t parno=1;parno<5;parno++) {
00193          params[parno]=fitty->GetParameter(parno);
00194          parerrs[parno]=fitty->GetParError(parno);
00195       }
00196       for(Int_t parno=6;parno<8;parno++) {
00197          params[parno]=fitty->GetParameter(parno);
00198          parerrs[parno]=fitty->GetParError(parno);
00199       }
00200    }
00201    else {
00202       for(Int_t fred=0; fred<9;fred++) {
00203          params[fred]=-100;
00204          parerrs[fred]=-100;
00205       }
00206    }
00207    if(drawIt) {
00208        peddy->SetLineColor(kBlue);
00209        histo->Draw();
00210        peddy->Draw("SAME");
00211        
00212        if(meanpe>0.1) {
00213            fitty->SetLineColor(kRed);
00214            fitty->SetRange(plotmin,plotmax);
00215            fitty->Draw("SAME");
00216        }
00217    }
00218    else {
00219        delete peddy;
00220        delete fitty;  
00221    }
00222 /*    cout << histo->GetName() << "\t" */
00223 /*      << histo->GetMean() << "\t" */
00224 /*      << histo->GetRMS() << "\t" */
00225 /*      << histo->GetEntries() << "\t"; */
00226 /*    for(int i=0;i<8;i++) { */
00227 /*        if(i==0) { */
00228 /*         cout << peddy->GetParameter(0) << "\t" */
00229 /*              << peddy->GetParError(0) << "\t"; */
00230 /*        } */
00231 /*        else if(i==5) { */
00232 /*         cout << peddy->GetParameter(1) << "\t" */
00233 /*              << peddy->GetParError(1) << "\t"; */
00234 /*        } */
00235 /*        else if(i==7) { */
00236 /*         cout << peddy->GetParameter(2) << "\t" */
00237 /*              << peddy->GetParError(2) << "\t"; */
00238 /*        } */
00239 /*        else { */
00240 /*         cout << fitty->GetParameter(i) << "\t" */
00241 /*              << fitty->GetParError(i) << "\t"; */
00242 /*        } */
00243 /*    }     */
00244 /*    cout << endl; */
00245 
00246    
00247 }

Double_t my_gamma ( Double_t *  x,
Double_t *  par 
)

Definition at line 301 of file theFit.h.

00302 {
00303    Double_t top= pow(x[0],par[1]-1)*exp(-x[0]/par[2]);
00304    Double_t bottom= pow(par[2],par[1])*TMath::Gamma(par[1]);
00305 
00306 
00307 
00308    return par[0]*top /bottom;
00309 }

Double_t myfped_1gauss ( Double_t *  x,
Double_t *  par 
)

Definition at line 288 of file theFit.h.

Referenced by fit_histo().

00289 {
00290    static double r2pi = 2.506628275;
00291    static double r2 = 1.414213562;
00292    
00293    Double_t s = r2*par[1];
00294    Double_t y = 0.5*par[2]*(TMath::Erf(Double_t(x[0]-par[0]+.5)/s)
00295                             -TMath::Erf(Double_t(x[0]-par[0]-.5)/s));
00296    return y;       
00297          
00298 }

Double_t mypefit ( Double_t *  x,
Double_t *  par 
)

Definition at line 249 of file theFit.h.

References MuELoss::a.

Referenced by fit_histo().

00250 {
00251    static double r2pi = 2.506628275;
00252    static double r2 = 1.414213562;
00253 
00254    double a,s,y;
00255   
00256    /* Pedestal p.d.f.: double Gaussian. */
00257    a = par[7];
00258    a = a*exp(-par[1]);
00259    s = r2*par[5];
00260    y = 0.5*a*(TMath::Erf(Double_t(x[0]-par[0]+.5)/s)
00261                                   -TMath::Erf(Double_t(x[0]-par[0]-.5)/s));
00262 
00263    
00264    /*Photoelectron p.d.f.: combination of Poisson and Gaussian.*/
00265    a=par[7]*exp(-par[1])*(1-par[4]);
00266    a = a/r2pi;
00267    s=par[5];
00268    for (int j=1; j<12; j++)
00269       {
00270          a = a*par[1]/double(j);
00271          s = sqrt(pow(par[3]*par[2],2)+pow(s,2));
00272          y += a/s*exp(-.5*pow(((double(x[0])-double(j)*par[2]-par[0])/s),2));
00273       }
00274    
00275    
00276    /* Dynode PE: assumed to be a gaussian with the mean and the width
00277       1/10 of that of single PE. */
00278    a = par[7]*par[4]*(1-exp(-par[1]));
00279    s = par[5];
00280    s = sqrt( pow(par[3]*par[2]/par[6],2) + pow(s,2) );
00281 
00282    y += a/(s*r2pi)*exp(-.5*pow(((double(x[0])-par[2]/par[6]-par[0])/s),2));
00283    return y;
00284    
00285 }


Generated on 2 Nov 2017 for loon by  doxygen 1.6.1