JFCalibration/JFMC_CuFluoPeakFit.cpp

// file to fit fluo spectrum per pixel
// and save peak position and uncertainty

#include "/afs/psi.ch/project/mythen/sophie/sls_detector_calibration/jungfrauCommonHeader.h"
#include "/afs/psi.ch/project/mythen/sophie/sls_detector_calibration/jungfrauCommonFunctions.h"

#include "/afs/psi.ch/project/mythen/sophie/sls_detector_calibration/jungfrauFile.C"
#include "/afs/psi.ch/project/mythen/sophie/sls_detector_calibration/jungfrauPixelMask.C"
#include "/afs/psi.ch/project/mythen/sophie/sls_detector_calibration/jungfrauPedestal.C"

#include "/afs/psi.ch/project/mythen/sophie/sls_detector_calibration/energyCalibration.h"
#include "/afs/psi.ch/project/mythen/sophie/sls_detector_calibration/energyCalibration.cpp"

#include "TGraph.h"
#include "TGraphErrors.h"
#include "TF1.h"
#include "TFile.h"
#include "TPaveStats.h"

int main(int argc, char* argv[]) {
  jungfrauStyle();
  gStyle->SetOptFit(11);

  if (argc != 3) {
    cout << "Correct usage:" << endl;
    cout << "arg 1: specify module number" << endl;
    cout << "arg 2: specify data location" << endl;
    cout << "eg: ./JFMC_FluoCalibFit 006 /data_pool/Module_006_161116" << endl;
    cout << " " << endl;
    exit(1);
  }

  string module_str = argv[1];
  string data_loc = argv[2];

  char histoname[128];
  char savename[128];
  TCanvas* c1 = new TCanvas("c1","");

  sprintf(savename,"%s_ana_M%s/Fluo_comb.root", data_loc.c_str(), module_str.c_str());
  TFile* comb_file = new TFile((const char *)(savename),"READ");

  TH1F* peak_fit_pos = new TH1F("peak_fit_pos","",100,200,400);
  TH1F* peak_fit_poserr = new TH1F("peak_fit_poserr","",100,0,2);
  TH2F* peak_fit_pos_2d = new TH2F("peak_fit_pos_2d","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
  TH2F* peak_fit_poserr_2d = new TH2F("peak_fit_poserr_2d","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);

  for (int j=1; j<9; j++) {
    cout << "slice " << j << endl;
    sprintf(histoname,"adc2d_%i",j);
    TH2I* adc2d_j = (TH2I*)(comb_file->Get((const char *)(histoname)));
    cout << adc2d_j->GetEntries() << endl;
    adc2d_j->Draw("colz");
    c1->Update();

    for (int i=(65536*(j-1)); i<(65536*(j)); i++) {
      if (i%10000==0){cout << "another 10k" << endl;}

	TH1D* proj = adc2d_j->ProjectionX("bin1",i-(65536*(j-1))+1,i-(65536*(j-1))+1);
	proj->Draw();
	if (proj->Integral(300,600)!=0) {

	  proj->Rebin(4);
	  proj->GetXaxis()->SetRangeUser(100,400);
	  proj->SetStats(kTRUE);

	  //Double_t mypar[8];
	  Double_t mypar[6];
	  mypar[0] = 0.0;
	  mypar[1] = 0.0;
	  mypar[2] = proj->GetBinCenter(proj->GetMaximumBin());
	  //mypar[6] = 1.127;
	  //mypar[7] = 0.2;
	  mypar[3] = 19.;
	  mypar[4] = 60.;
	  mypar[5] = 0.18;
	  
	  //Double_t emypar[8];
	  Double_t emypar[6];
	  energyCalibration *thiscalibration = new energyCalibration();
	  thiscalibration->setScanSign(1); // new
	  thiscalibration->setStartParametersKb(mypar);
	  thiscalibration->fixParameter(0,0.); // no background
	  thiscalibration->fixParameter(1,0.); // new
	  //TF1* fittedfun = thiscalibration->fitSpectrumKb(proj,mypar,emypar);
	  TF1* fittedfun = thiscalibration->fitSpectrum(proj,mypar,emypar);

	  peak_fit_pos->Fill(mypar[2]);
	  peak_fit_poserr->Fill(emypar[2]);
	  peak_fit_pos_2d->Fill(i%NC,i/NC,mypar[2]);
	  peak_fit_poserr_2d->Fill(i%NC,i/NC,emypar[2]);

	  if (i > 10000 && i < 10010) {
	    proj->Draw();
	    fittedfun->Draw("same");
	    c1->Update();
	    proj->GetXaxis()->SetTitle("Pedestal corrected ADC [ADU]");
	    TPaveStats *st = (TPaveStats*)proj->FindObject("stats");
	    st->SetX1NDC(0.22);
	    st->SetX2NDC(0.55);
	    st->SetY1NDC(0.55);
	    st->SetY2NDC(0.94);
	    st->SetBorderSize(0);
	    sprintf(savename,"plots/M%s/CuFluo/proj_test%d.png",module_str.c_str(),i);
	    c1->SaveAs((const char *)(savename));
	    c1->SetLogy();
	    sprintf(savename,"plots/M%s/CuFluo/proj_test%d_log.png",module_str.c_str(),i);
	    c1->SaveAs((const char *)(savename));
	    c1->SetLogy(0);
	  }

	  delete thiscalibration;

	}
    }
  }

  TCanvas *mapcanvas = new TCanvas("mapcanvas","",150,10,800,400);
  mapcanvas->SetLeftMargin(0.1);
  mapcanvas->SetRightMargin(0.13);
  mapcanvas->SetTopMargin(0.08);
  mapcanvas->SetBottomMargin(0.15);

  peak_fit_pos_2d->GetXaxis()->SetTitle("Column");
  peak_fit_pos_2d->GetYaxis()->SetTitle("Row");
  peak_fit_pos_2d->GetYaxis()->SetTitleOffset(0.7);
  peak_fit_pos_2d->Draw("colz");
  peak_fit_pos_2d->GetZaxis()->SetRangeUser(200,400);
  sprintf(savename,"plots/M%s/CuFluo/peak_fit_pos_2d.png",module_str.c_str());
  mapcanvas->SaveAs((const char *)(savename));

  peak_fit_poserr_2d->GetXaxis()->SetTitle("Column");
  peak_fit_poserr_2d->GetYaxis()->SetTitle("Row");
  peak_fit_poserr_2d->GetYaxis()->SetTitleOffset(0.7);
  peak_fit_poserr_2d->Draw("colz");
  peak_fit_poserr_2d->GetZaxis()->SetRangeUser(0,3);
  sprintf(savename,"plots/M%s/CuFluo/peak_fit_poserr_2d.png",module_str.c_str());
  mapcanvas->SaveAs((const char *)(savename));

  c1->cd();
  peak_fit_pos->GetXaxis()->SetTitle("Peak position [ADU]");
  peak_fit_pos->Draw();
  sprintf(savename,"plots/M%s/CuFluo/peak_fit_pos.png",module_str.c_str());
  c1->SaveAs((const char *)(savename));

  peak_fit_poserr->GetXaxis()->SetTitle("Peak position uncert [ADU]");
  peak_fit_poserr->Draw();
  sprintf(savename,"plots/M%s/CuFluo/peak_fit_poserr.png",module_str.c_str());
  c1->SaveAs((const char *)(savename));

  sprintf(savename,"data/M%s/CuF_peakpos_M%s.root", module_str.c_str(), module_str.c_str());
  TFile* saved_file = new TFile((const char *)(savename),"RECREATE");
  peak_fit_pos_2d->Write();
  peak_fit_poserr_2d->Write();
  saved_file->Close();
  
}