mbd/WAVEFORMS/find_th2ridge.C

0001 #ifndef MACRO_FIND_TH2RIDGE_C
0002 #define MACRO_FIND_TH2RIDGE_C
0003
0004 #include <TF1.h>
0005 #include <TFitResult.h>
0006 #include <TGraphErrors.h>
0007 #include <TH2.h>
0008 #include <TPad.h>
0009
0010 // find the ridge of the TH2
0011 TGraphErrors *find_th2ridge(const TH2* h2)
0012 {
0013   int nbinsx = h2->GetNbinsX();
0014   int nbinsy = h2->GetNbinsY();
0015 //  double min_xrange = h2->GetXaxis()->GetBinLowEdge(1);
0016 //  double max_xrange = h2->GetXaxis()->GetBinLowEdge(nbinsx+1);
0017   double min_yrange = h2->GetYaxis()->GetBinLowEdge(1);
0018   double max_yrange = h2->GetYaxis()->GetBinLowEdge(nbinsy+1);
0019
0020   TString name;
0021   TString title;
0022   name = "aaa";
0023   //title = "aaa";
0024   TGraphErrors *prof = new TGraphErrors();
0025   prof->SetName(name);
0026   //prof->SetTitle(title);
0027
0028   TH1 *h_projx = h2->ProjectionX("projx");
0029
0030   TF1 gaussian("gaussian","gaus",min_yrange,max_yrange);
0031   gaussian.SetLineColor(4);
0032
0033   TH1 *h_projy{nullptr};
0034   double adcmean = 0.;  // x-value
0035 //  double adcnum = 0.;
0036
0037   for (int ibin=1; ibin<=nbinsx; ibin++)
0038   {
0039     name = "hproj_"; name += ibin;
0040     if ( h_projy==nullptr )
0041     {
0042       h_projy = h2->ProjectionY(name,ibin,ibin);
0043       adcmean = h_projx->GetBinCenter(ibin);
0044 //      adcnum = 1.0;
0045     }
0046     /*
0047     else  // keep summing until we get enough statistics in bin
0048     {
0049       TH1 *h_projyadd =  h2->ProjectionY(name,ibin,ibin);
0050       h_projy->Add( h_projyadd );
0051       delete h_projyadd;
0052
0053       adcmean += h_projx->GetBinCenter(ibin);
0054       adcnum += 1.0;
0055     }
0056     */
0057
0058     std::cout << "adcmean " << adcmean << std::endl;
0059
0060     if ( h_projy->Integral()>100 || ibin==nbinsx )
0061     {
0062       //adcmean = adcmean/adcnum;
0063
0064       h_projy->Draw();
0065
0066       int maxbin = h_projy->GetMaximumBin();
0067       double xmax = h_projy->GetBinCenter( maxbin );
0068       double ymax = h_projy->GetBinContent( maxbin );
0069       gaussian.SetParameter(1,xmax);
0070       gaussian.SetParameter(0,ymax);
0071       gaussian.SetRange(xmax-0.6,xmax+0.6);
0072
0073       std::cout << "xmax " << xmax << "\t" << ymax << "\t" << h_projy->GetRMS() << std::endl;
0074       auto fitresult = h_projy->Fit("gaussian","RWWS");
0075
0076       int n = prof->GetN();
0077
0078       if ( fitresult->IsValid() )
0079       {
0080         double mean = gaussian.GetParameter(1);
0081         double meanerr = gaussian.GetParError(1);
0082
0083         prof->SetPoint(n,adcmean,mean);
0084         prof->SetPointError(n,0,meanerr);
0085         //std::cout << "mean and meanerr " << mean << "\t" << meanerr << std::endl;
0086       }
0087       else
0088       {
0089         prof->SetPoint(n,adcmean,h_projy->GetMean());
0090         prof->SetPointError(n,0,h_projy->GetRMS());
0091       }
0092
0093       gPad->Modified();
0094       gPad->Update();
0095
0096       /*
0097       string junk;
0098       cin >> junk;
0099       */
0100
0101       delete h_projy;
0102       h_projy = nullptr;
0103     }
0104
0105   }
0106
0107   delete h_projx;
0108
0109   prof->SetBit(TGraph::kIsSortedX);
0110   return prof;
0111 }
0112
0113 #endif