obsolete/cc/FitHisto_Iterative.C

0001
0002 void FitHisto_Iterative(const char *fname, const char *hname, double peakEx, int nsigma){
0003   auto f = new TFile(fname);
0004   auto h = (TH1D *) f->Get(hname);
0005   if(!f || !h) printf("No File or Histo Found\n");
0006   float tolerance = 2;
0007   float startingvalue = h->GetXaxis()->GetBinCenter(1);
0008   printf("Starting value of bin: %f\n",startingvalue);
0009   double binwidth= h->GetBinWidth(1);
0010   double nbins = h->GetNbinsX();
0011   //Now search for the peak in the histo and test its goodness
0012   double peakLocation = h->GetMaximumBin();
0013   if(TMath::Abs((peakLocation*binwidth)+startingvalue - peakEx)>tolerance) {
0014     printf("Peak found in Unnatural location \n");
0015     printf("## %f compared to %f\n",(peakLocation*binwidth)+startingvalue, peakEx);
0016   }
0017   else{printf("peak found in good region; continuing ...\n");}
0018
0019   //Now the Fitting part; well at this level we don't care much of percision; Max bin --> +- half max
0020   // Get the first bin with diff = (BinContent - DesiredBinContent)<=maxdiff;
0021   int l=0, u=0;
0022   double difftol = 100;   // once should tune
0023   double halfcontent =(h->GetBinContent(h->GetMaximumBin())/2);
0024   printf("maxconetnt %f %f %f\n", h->GetBinContent(h->GetMaximumBin()),halfcontent, difftol );
0025   h->GetBinWithContent(halfcontent,l,1,peakLocation,difftol); //100 tolerance
0026   h->GetBinWithContent(halfcontent,u,peakLocation,nbins,difftol); //100 tolerance
0027   printf("First Estimate of Half Width %f %f\n ",h->GetBinContent(l), h->GetBinContent(u));
0028
0029   //Geting Bins to run first fit
0030   float lowedge = (l*binwidth)+startingvalue;
0031   float upperedge = (u*binwidth)+startingvalue;
0032
0033   auto f1 = new TF1("f1","gaus",200,300);
0034   h->Fit(f1,"R","",lowedge,upperedge);
0035   //auto f1 = h->GetFunction("gaus");
0036
0037   double mean =0;
0038   mean = f1->GetParameter(1); double sigma = 0;
0039   sigma = f1->GetParameter(2);
0040   f1->SetParameter(1,mean);
0041   f1->SetParameter(2,sigma);
0042
0043   h->Fit(f1,"R","",mean-2*sigma,mean+2*sigma);
0044   mean = f1->GetParameter(1);sigma = f1->GetParameter(2);
0045   f1->SetParameter(1,mean);
0046   f1->SetParameter(2,sigma);
0047   vector<double> SNR;
0048   vector<double> Sigma;
0049   auto cc1 = new TCanvas(); cc1->Divide(3,2);
0050   for(int i=1; i<=nsigma; i++){
0051     //cc1->cd(i);
0052     //h->Draw();
0053     float r1 =  mean - i*sigma;
0054     float r2 =  mean + i*sigma;
0055     printf("Fitting From %f %f \n",r1,r2);
0056     h->Fit(f1,"R","",r1,r2);
0057     mean = f1->GetParameter(1);sigma = f1->GetParameter(2);
0058     f1->SetParameter(1,mean);
0059     f1->SetParameter(2,sigma);
0060
0061     SNR.push_back(mean/sigma);
0062     Sigma.push_back(i);
0063   }
0064   auto cc = new TCanvas();
0065   auto G = new TGraph(5,&Sigma[0],&SNR[0]);
0066   G->Draw("AP*");
0067 }