Tspectrum1.C

#include "TCanvas.h"
#include "TMath.h"
#include "TGraph.h"
#include "TH1.h"
#include "TF1.h"
#include "TRandom.h"
#include "TSpectrum.h"
#include "TVirtualFitter.h"
#include "unistd.h"
#include "TNamed.h"
#include "TAttLine.h"
#include "TAttFill.h"
#include "TAttMarker.h"
#include "TVectorFfwd.h"
#include "TVectorDfwd.h"


Int_t npeaks = 3;
Double_t fpeaks(Double_t *x, Double_t *par) {
  Double_t result = par[0] + par[1]*x[0];
  for (Int_t p=0;p<npeaks;p++) {
    Double_t norm  = par[3*p+2];
    Double_t mean  = par[3*p+3];
    Double_t sigma = par[3*p+4];
    result += norm*TMath::Gaus(x[0],mean,sigma);
  }
  return result;
}
void Tspectrum1()
{


  Double_t StdDev;
  Double_t Mean;
  TString index_s;
  TString index_r;
  TString filename1 ="~/lmartin/agdaq_laser/ana/output02305.root"; 
  TString object = "h_pad_amp_per_row_cut_col";
  TFile *f = new TFile(filename1);
 
 
  
  TCanvas *c2 = new TCanvas("c2","A Graph ",200,10,700,500);
  TH1D* hdummy= new TH1D ("hdummy", "peaks graph", 1, 0., 576.);
  hdummy ->SetStats(kFALSE); //gets rid of the box 
  hdummy ->GetYaxis() ->SetRangeUser(0.,20000.);
  hdummy ->Draw();
  double maxpeak=0.;
  c2->SetFillColor(kGray);
  c2->SetGrid();
  

  for(Int_t i=0; i<33; i++){
    if ( i >5 && i<23){
      cout << "unwanted region" << endl;
    }
    else if (i==1){
      cout << "unwanted region" << endl;
    }
    else if (i==2){
      cout << "unwanted region" << endl;
    }
    else {
      index_r.Form("%i",i);
      // TString filename = object+index_r;
      //filename = TString::Format("%s%02i",object.Data(),i);
      // filename = TString::Format("h_pad_amp_per_row_cut_col%02i",i);
      TString filename = object+ TString::Format("%02i",i);

      TCanvas *c1 = new TCanvas("c1"+index_r,"Profile histogram example",1200,1000); 
      c1 ->Draw();
      TH2D *hfix = (TH2D*) f->FindObjectAny(filename); 
      cout<<filename<<endl;
      Int_t NX = hfix-> GetNbinsX();
      cout << "number of bins " << NX <<endl;
   

      vector<double> x, y;

      //hfix->GetXaxis()->SetRangeUser(0,400);
      //if ( i >6 && i<20){
      //cout << "unwanted region" << endl; 
      /*
	for (Int_t r=0; r < NX ;r++){
	index_s.Form("%i",r);
	TH1D *py = hfix->ProjectionY("py"+index_s+"col"+index_r,r,r,"d"); 
	py ->SetFillColor(kBlue-2);
	StdDev = py ->GetStdDev();
	Mean = py ->GetMean();

	if(StdDev < 500 && StdDev !=0 && py ->GetEntries()> 50 ){

	//Use TSpectrum to find the peak candidates
	x.push_back(r);
	y.push_back(Mean); 
	//printf( " %f %f \n" , x.back(), y.back()); 
	
	TSpectrum *s = new TSpectrum(2);
	Int_t nfound = s->Search(py,20,"",0.1);
	cout << "Found "<< nfound << " candidate peaks to fit." << endl;
	Double_t *xpeaks = s->GetPositionX();

	Double_t *ypeaks = s->GetPositionY();
	//cout << "xpeaks " << xpeaks <<endl;
	cout << "Now fitting: Be patient.\n";
	TF1 *f1 = new TF1("f1",fpeaks,0,1000,2+3*nfound);
	Double_t par[2000];
	par[0] = 0.1;
	par[1] = 1.;
	for (int p=0;p<nfound;p++) {
	Double_t xp = xpeaks[p];
	Double_t yp = ypeaks[p];
	par[3*p+2] = yp;
	par[3*p+3] = xp;
	par[3*p+4] = 50;
	}
      

	// TVirtualFitter::Fitter(hfix,10+3*nfound);
	f1->SetParameters(par);
	f1->FixParameter(1,0);
	//cout << "parameters " << f1->GetParameter(0) << ", " << f1->GetParameter(1) << ", " << f1->GetParameter(2) << ", " << f1->GetParameter(3) << ", " << f1->GetParameter(4) << endl;
	f1->SetNpx(1000);
	py ->Draw(); 
	py->Fit("f1");
	cout << "Now fitting: Be patient.\n";

	// cout<<index_s<<endl; 
	//cout<<StdDev<<endl;
	c1 ->Update();
	// sleep(1);

	} else if(py ->GetEntries()>1000 && StdDev != 0){
	cout << r << " : StdDev = " << StdDev << " py ->GetEntries() = " <<  py ->GetEntries() << endl;
	py ->Draw();
	c1 ->Update();
	//sleep(10);	
	}
	
	
      */
    
      //}
      //else{
      // for (Int_t r=300; r < NX ;r++){
      for (Int_t r=0; r < NX ;r++){
	index_s.Form("%i",r);
	TH1D *py = hfix->ProjectionY("py"+index_s+"col"+index_r,r,r,"d"); 
	py ->SetFillColor(kBlue-2);
	StdDev = py ->GetStdDev();
	Mean = py ->GetMean();

	if(StdDev < 500 && StdDev !=0 && py ->GetEntries()> 50 ){
	  
	  //Use TSpectrum to find the peak candidates
	  x.push_back(r);
	  y.push_back(Mean); 
	  printf( " %f %f \n" , x.back(), y.back()); 
	 
	  TSpectrum *s = new TSpectrum(2);
	  Int_t nfound = s->Search(py,20,"",0.1);
	  cout << "Found "<< nfound << " candidate peaks to fit." << endl;
	  Double_t *xpeaks = s->GetPositionX();
	  Double_t *ypeaks = s->GetPositionY();
	  //cout << "xpeaks " << xpeaks <<endl;
	  cout << "Now fitting: Be patient.\n";
	  TF1 *f1 = new TF1("f1",fpeaks,0,1000,2+3*nfound);
	  Double_t par[2000];
	  par[0] = 0.1;
	  par[1] = 1.;
	  for (int p=0;p<nfound;p++) {
	    Double_t xp = xpeaks[p];
	    Double_t yp = ypeaks[p];
	    par[3*p+2] = yp;
	    par[3*p+3] = xp;
	    par[3*p+4] = 50;
	  }
      

	  // TVirtualFitter::Fitter(hfix,10+3*nfound);
	  f1->SetParameters(par);
	  f1->FixParameter(1,0);
	  //cout << "parameters " << f1->GetParameter(0) << ", " << f1->GetParameter(1) << ", " << f1->GetParameter(2) << ", " << f1->GetParameter(3) << ", " << f1->GetParameter(4) << endl;
	  f1->SetNpx(1000);
	  py ->Draw(); 
	  py->Fit("f1");
	  cout << "Now fitting: Be patient.\n";

	  // cout<<index_s<<endl; 
	  //cout<<StdDev<<endl;
	  c1 ->Update();
	  //sleep(1);

	} else if(py ->GetEntries()>1000 && StdDev != 0){
	  cout << r << " : StdDev = " << StdDev << " py ->GetEntries() = " <<  py ->GetEntries() << endl;
	  py ->Draw();
	  c1 ->Update();
	  //sleep(10);	
	}
      }
      
      
    
      // TH1D* pad_peak_per_row;
      // pad_peak_per_row = new TH1D();
      // pad_peak_per_row -> Fill (x[n],y[n]);

      //new TCanvas;

      // TH1* h1 = new TH1D("h1", "peaks", 100, 0.0, 4.0);
      //h1 ->Fill(x.data(),y.data());
      // TCanvas *c2 = new TCanvas("c2","A Graph ",200,10,700,500);
      //c2->SetFillColor(42);
      //c2->SetGrid();



      TGraph *gr = new TGraph(x.size(),x.data(),y.data());
      gr ->SetName("col" +index_r);
      gr->SetLineColor(i-20);
      gr->SetLineWidth(1);
      gr->SetMarkerColor(i-20);
      gr->SetMarkerStyle(21);
      gr->SetTitle("col" +index_r);
      gr->GetXaxis()->SetTitle("pad rows");
      gr->GetYaxis()->SetTitle("amp");
      //maxpeak=maxpeak>*std::max_element(y.begin(),y.end())?maxpeak:*std::max_element(y.begin(),y.end());
  
      //c2 ->cd();
      //gr->Draw("CP SAME");
      //c2 ->Update();
  
   
      TH1* h1 = new TH1D("h1 col " +index_r, "peaks", NX, -0.5, NX-0.5);
      // h1 ->Fill(y.data());
      auto nPoints = gr->GetN();
      //cout << "number of Points is: " << nPoints << endl;
      for(int i=0; i < nPoints; ++i) {
	double x,y;
	gr->GetPoint(i, x, y);
	cout << x << y <<endl;
	h1->Fill(x,y); 
      }
      // h1 ->Draw();


      TSpectrum *s = new TSpectrum(9);
      Int_t nfound = s->Search(h1,1,"",0.05);
      cout << "Found "<< nfound << " candidate peaks to fit." << endl;
      Double_t *x1peaks = s->GetPositionX();
      Double_t *y1peaks = s->GetPositionY();
      cout << "x1peaks " << x1peaks <<endl;
      cout << "Now fitting: Be patient.\n";
      TF1 *f2 = new TF1("f2",fpeaks,0,1000,2+3*nfound);
      Double_t par[2000];
      par[0] = 0.1;
      par[1] = 1.;
      for (int p=0;p<nfound;p++) {
	Double_t x1p = x1peaks[p];
	Double_t y1p = y1peaks[p];
	par[3*p+2] = y1p;
	par[3*p+3] = x1p;
	par[3*p+4] = 2;
      }

      c2 ->cd();
      // new TCanvas; 
      f2->SetParameters(par);
      f2->FixParameter(1,0);
      f2->FixParameter(0,0);
      //cout << "parameters " << f1->GetParameter(0) << ", " << f1->GetParameter(1) << ", " << f1->GetParameter(2) << ", " << f1->GetParameter(3) << ", " << f1->GetParameter(4) << endl;
      f2->SetNpx(1000);
      //h1 ->Draw ();
      h1 ->SetStats(kFALSE); //gets rid of the box 
      h1 ->GetYaxis() ->SetRangeUser(0.,20000.);
      h1 ->Draw((i==0)?"ALP":"LP same");
      h1->Fit("f2");
      cout << "Now fitting: Be patient.\n";

      // cout<<index_s<<endl; 
      //cout<<StdDev<<endl;
      //c3 ->Update();
   
   
      //h1 ->Fit ("f1","R");
      //gr->Draw("CP SAME");
      //gr->Draw();
      //c2 ->Update();



    
      gr->Draw((i==0)?"ALP":"LP same");
      //gr->Draw("LP same");
  
  
      //c2->Update();
  
    }
  }
  //hdummy ->GetYaxis() ->SetRangeUser(0.,maxpeak*1.1);
  //c2->BuildLegend();
  //c2->Update();
}