loon: theFit.h File Reference

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 }

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)
theFit.h File Reference

Functions

Function Documentation
