// fit fluo spectrum per pixel for a given storage cell, only G0
// (can be called multiple times to handle multi-storage-cell data, each time with different Sc number as argument, i.e. supports manual multi-threading on the same data set)
// and save peak position and uncertainty
// Takes pedestal-corrected data (in root file) as input
// Call CuFluo_analysis_sc_singlethread_data before to analyse raw data files and create root file

#include "../sls_detector_calibration/jungfrauCommonHeader.h"
#include "../sls_detector_calibration/jungfrauCommonFunctions.h"

#include "../sls_detector_calibration/jungfrauFile.C"
#include "../sls_detector_calibration/jungfrauPixelMask.C"
#include "../sls_detector_calibration/jungfrauPedestal.C"

#include "../sls_detector_calibration/energyCalibration.h"
#include "../sls_detector_calibration/energyCalibration.cpp"

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

#include <sys/stat.h>
#include <vector>
#include <iostream>
#include <string>

int main(int argc, char* argv[]) {

  jungfrauStyle();
  gStyle->SetOptFit(11);

  if (argc != 8) {
    cout << "Correct usage:" << endl;
    cout << "arg 1: specify module number" << endl;
    cout << "arg 2: specify HG0 or G0" << endl;
    cout << "arg 3: specify data location folder" << endl;
    cout << "arg 4: specify pede file prefix" << endl;
    cout << "arg 5: specify data file prefix" << endl;
    cout << "arg 6: specify location of analysis root files" << endl;
    cout << "arg 7: specify storage cell number" << endl;
    cout << " " << endl;
    exit(1);
  }

  string module_str = argv[1];
  string gain_str = argv[2];
  bool isJF11=false;
  if (gain_str == "HG0JF11") {
    gain_str = "HG0";
    isJF11=true;
    }

  string data_loc = argv[3];
  string pede_file = argv[4];
  string data_file = argv[5];
  string anadata_loc = argv[6];

  bool save_to_afs = false;

  int createHistoFile = 1;

  char histoname[256];
  char savename[256];
  char plotfolder[256];
  char rootdatafolder[256];
  int filen = 352;
  int pedefilen = 5;
  int this_storagecell = stoi(argv[7]); // needs compiler flag -std=c++11

  // create necessary directories with permissions drwxrwxr-x
  // data/Mxxx
  sprintf(savename,"data/M%s", module_str.c_str());
  mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
  // plots/Mxxx
  sprintf(savename,"plots/M%s", module_str.c_str());
  mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
  // plots/Mxxx/CuFluo
  sprintf(savename,"plots/M%s/CuFluo", module_str.c_str());
  mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
  // plots/Mxxx/CuFluo/G0
  sprintf(savename,"plots/M%s/CuFluo/G0", module_str.c_str());
  mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
  // plots/Mxxx/CuFluo/HG0
  sprintf(savename,"plots/M%s/CuFluo/HG0", module_str.c_str());
  mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
  sprintf(savename,"%s", anadata_loc.c_str());
  mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);

  if ( save_to_afs ) {

    cout << "Saving to afs." << endl;

    sprintf( plotfolder, "plots/M%s/CuFluo/%s", module_str.c_str(), gain_str.c_str() );
    sprintf( rootdatafolder, "data/M%s", module_str.c_str() );

  } else {

    cout << "Not saving to afs." << endl;

    // data/Mxxx
    sprintf( savename,"%s/data", anadata_loc.c_str() );
    mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
    // plots/Mxxx
    sprintf( savename,"%s/plots", anadata_loc.c_str() );
    mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
    // plots/Mxxx/CuFluo
    sprintf( savename,"%s/plots/CuFluo", anadata_loc.c_str() );
    mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
    // plots/Mxxx/CuFluo/G0
    sprintf( savename,"%s/plots/CuFluo/G0", anadata_loc.c_str() );
    mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
    // plots/Mxxx/CuFluo/HG0
    sprintf( savename,"%s/plots/CuFluo/HG0", anadata_loc.c_str() );
    mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);

    sprintf( plotfolder, "%s/plots/CuFluo/%s", anadata_loc.c_str(), gain_str.c_str() );
    sprintf( rootdatafolder, "%s/data", anadata_loc.c_str() );

  }

  jungfrauPixelMask *pixelMaskObject = new jungfrauPixelMask();
  bool pixel_mask [NCH];
  pixelMaskObject->initialisePixelMask(pixel_mask);

  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);

  TCanvas* c1 = new TCanvas("c1","");

  /*****************************************************************************************************************/
  /* NOTE: Since arrays, as they were used before, tend to cause memory leaks if they grow too large on the stack, */
  /*       it might be smarter to realize all array objects with std::vector instead.                              */
  /*****************************************************************************************************************/

  if (createHistoFile) {

    jungfrauFile *thisfile = new jungfrauFile();

    jungfrauPedestal* pedestalObject_SC = new jungfrauPedestal();
    pedestalObject_SC->pedestalSetNFrames(100);

    for (int pedefilei = 0; pedefilei < pedefilen; ++pedefilei) {
    
      // open pede file
      sprintf( savename, "%s/%s_%%6.6d.dat", data_loc.c_str(), pede_file.c_str() );
      thisfile->open( (char*)savename, pedefilei );
	
      //count events in file
      int nevents = 0;
      while ( thisfile->readNextFrame() ) {
	nevents++;
      }
      thisfile->rewind();
      cout << "read " << nevents << " events" << endl;
	
      while ( thisfile->readNextFrame() ) {
	// calculate pixel mask
	pixelMaskObject->maskIfGainNot( 0, thisfile->getFrameDataHandle(), pixel_mask );
	// caluclate pedestals
	if( thisfile->currentSCnumber() == (uint64_t)this_storagecell ) {
	  pedestalObject_SC->addFrameToPedestalCalculation( thisfile->getFrameDataHandle() );
	}
      }
	
      thisfile->close();
      
    } // end of loops over files
      
    cout << "after chip mask, n masked pixels is " << pixelMaskObject->getNMasked(pixel_mask) << endl;

  } // end if i creatHistoFile

  // We have saved the data into root-file by first calling CuFluo_analysis_sc_singlethread_data.cpp. Now, we read from that same root-file and perform the fitting.

  c1->cd();

  sprintf(savename,"%s/CuFluo_%s_sc%d_file0to%d.root", anadata_loc.c_str(), gain_str.c_str(), this_storagecell, filen-1); //uncomment for VH 210906
  TFile* comb_file = new TFile((const char *)(savename),"READ");

  pixelMaskObject->initialisePixelMask(pixel_mask);

  int low_ADU_peak = 0;
  int high_ADU_peak = 0;
  if (gain_str == "HG0") {
    
    
    low_ADU_peak = 700;
    high_ADU_peak = 900;
    if  (isJF11) {
        low_ADU_peak = 850;
    high_ADU_peak = 1350;

    }

  } else if (gain_str == "G0") {
    low_ADU_peak = 200; // was 250
    high_ADU_peak = 400;
  }

  //declare
  TH1F* fit_par3;
  TH1F* fit_par4;
  TH1F* fit_par5;
  TH1F* fit_par6;
  TH1F* fit_par7;

  TH2F* fit_par3_2d;
  TH2F* fit_par4_2d;
  TH2F* fit_par5_2d;
  TH2F* fit_par6_2d;
  TH2F* fit_par7_2d;

  TH1F* peak_fit_pos;
  TH1F* peak_fit_poserr;
  TH2F* peak_fit_pos_2d;
  TH2F* peak_fit_poserr_2d;

  TH1F* noise_fit_pos;
  TH1F* noise_fit_poserr;
  TH2F* noise_fit_pos_2d;
  TH2F* noise_fit_poserr_2d;

  TH1F* gain_fit;
  TH1F* gain_fiterr;
  TH2F* gain_fit_2d;
  TH2F* gain_fiterr_2d;

  TH1F* gain_fit_isEdge;
  TH1F* gain_fit_isInnerEdge;
  TH1F* gain_fit_isDouble;
  TH1F* gain_fit_isNextToDouble; 
  TH1F* gain_fit_isQuad;
  TH1F* gain_fit_isBulk;

  TH2F* gain_ADUper1keV_2d;
  TH2F* gainerr_ADUper1keV_2d;

  //initialize
  Char_t *fithistoname = new Char_t[50];
  snprintf( fithistoname, 50, "fit_par3_sc%d", this_storagecell );
  fit_par3 = new TH1F( fithistoname, "", 100, 0, 50 );
  snprintf( fithistoname, 50, "fit_par4_sc%d", this_storagecell );
  fit_par4 = new TH1F( fithistoname, "", 100, 0, 500 );
  snprintf( fithistoname, 50, "fit_par5_sc%d", this_storagecell );
  fit_par5 = new TH1F( fithistoname, "", 100, 0, 0.5 );
  snprintf( fithistoname, 50, "fit_par6_sc%d", this_storagecell );
  fit_par6 = new TH1F( fithistoname, "", 100, 1.05, 1.3 ); // was 1.05, 1.25
  snprintf( fithistoname, 50, "fit_par7_sc%d", this_storagecell );
  fit_par7 = new TH1F( fithistoname, "", 100, 0, 0.4 );

  snprintf( fithistoname, 50, "fit_par3_2d_sc%d", this_storagecell );
  fit_par3_2d = new TH2F( fithistoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 );
  snprintf( fithistoname, 50, "fit_par4_2d_sc%d", this_storagecell );
  fit_par4_2d = new TH2F( fithistoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 );
  snprintf( fithistoname, 50, "fit_par5_2d_sc%d", this_storagecell );
  fit_par5_2d = new TH2F( fithistoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 );
  snprintf( fithistoname, 50, "fit_par6_2d_sc%d", this_storagecell );
  fit_par6_2d = new TH2F( fithistoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 );
  snprintf( fithistoname, 50, "fit_par7_2d_sc%d", this_storagecell );
  fit_par7_2d = new TH2F( fithistoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 );
   
  snprintf( fithistoname, 50, "peak_fit_pos_sc%d", this_storagecell );
  peak_fit_pos       = new TH1F( fithistoname, "", 100, low_ADU_peak, high_ADU_peak );
  snprintf( fithistoname, 50, "peak_fit_poserr_sc%d", this_storagecell );
  peak_fit_poserr    = new TH1F( fithistoname, "", 100, 0, 2 );
  snprintf( fithistoname, 50, "peak_fit_pos_2d_sc%d", this_storagecell );
  peak_fit_pos_2d    = new TH2F( fithistoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 );
  snprintf( fithistoname, 50, "peak_fit_poserr_2d_sc%d", this_storagecell );
  peak_fit_poserr_2d = new TH2F( fithistoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 );

  snprintf( fithistoname, 50, "noise_fit_pos_sc%d", this_storagecell );
  noise_fit_pos       = new TH1F( fithistoname, "", 100, -10, 10 );
  snprintf( fithistoname, 50, "noise_fit_poserr_sc%d", this_storagecell );
  noise_fit_poserr    = new TH1F( fithistoname, "", 100, 0, 0.15 ); // was 0, 0.1
  snprintf( fithistoname, 50, "noise_fit_pos_2d_sc%d", this_storagecell );
  noise_fit_pos_2d    = new TH2F( fithistoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 );
  snprintf( fithistoname, 50, "noise_fit_poserr_2d_sc%d", this_storagecell );
  noise_fit_poserr_2d = new TH2F( fithistoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 );

  snprintf( fithistoname, 50, "gain_fit_sc%d", this_storagecell );
  gain_fit       = new TH1F( fithistoname, "", 100, low_ADU_peak, high_ADU_peak );
  snprintf( fithistoname, 50, "gain_fiterr_sc%d", this_storagecell );
  gain_fiterr    = new TH1F( fithistoname, "", 100, 0, 2 );
  snprintf( fithistoname, 50, "gain_fit_2d_sc%d", this_storagecell );
  gain_fit_2d    = new TH2F( fithistoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 );
  snprintf( fithistoname, 50, "gain_fiterr_2d_sc%d", this_storagecell ); 
  gain_fiterr_2d = new TH2F( fithistoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 );

  snprintf( fithistoname, 50, "gain_fit_isEdge_sc%d", this_storagecell );
  gain_fit_isEdge         = new TH1F( fithistoname, "", 100, low_ADU_peak, high_ADU_peak );
  snprintf( fithistoname, 50, "gain_fit_isInnerEdge_sc%d", this_storagecell );
  gain_fit_isInnerEdge    = new TH1F( fithistoname, "", 100, low_ADU_peak, high_ADU_peak );
  snprintf( fithistoname, 50, "gain_fit_isDouble_sc%d", this_storagecell );
  gain_fit_isDouble       = new TH1F( fithistoname, "", 100, low_ADU_peak, high_ADU_peak );
  snprintf( fithistoname, 50, "gain_fit_isNextToDouble_sc%d", this_storagecell );
  gain_fit_isNextToDouble = new TH1F( fithistoname, "", 100, low_ADU_peak, high_ADU_peak );
  snprintf( fithistoname, 50, "gain_fit_isQuad_sc%d", this_storagecell );
  gain_fit_isQuad         = new TH1F( fithistoname, "", 100, low_ADU_peak, high_ADU_peak );
  snprintf( fithistoname, 50, "gain_fit_isBulk_sc%d", this_storagecell );
  gain_fit_isBulk         = new TH1F( fithistoname, "", 100, low_ADU_peak, high_ADU_peak );

  snprintf( fithistoname, 50, "gain_ADUper1keV_2d_sc%d", this_storagecell );
  gain_ADUper1keV_2d    = new TH2F( fithistoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 );
  snprintf( fithistoname, 50, "gainerr_ADUper1keV_2d_sc%d", this_storagecell );
  gainerr_ADUper1keV_2d = new TH2F( fithistoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 );

  int low_bin_noise = 101;
  int high_bin_noise = 301;
  int low_bin_peak = 0;
  int high_bin_peak = 0;
  if (gain_str == "HG0") {
    low_bin_peak = 701;
    high_bin_peak = 1200;
    if  (isJF11) {
      low_bin_peak = 801;
      high_bin_peak = 1451;
    }

  } else if (gain_str == "G0") {
    low_bin_peak = 301;
    high_bin_peak = 651;
  }

  for ( int j = 1; j < 9; ++j ) {
    cout << "SC = " << this_storagecell << " slice " << j << endl;
    sprintf( histoname, "adc2d_%d_sc%d", j, this_storagecell );
    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);

      if ( proj->Integral( low_bin_noise, high_bin_noise ) != 0 && proj->Integral( low_bin_peak, high_bin_peak ) != 0 ) {

	// noise
	TH1D *proj_noise = dynamic_cast<TH1D*>(proj->Rebin(4,"proj_noise"));
	proj_noise->SetStats(kTRUE);
	proj_noise->GetXaxis()->SetRangeUser(proj->GetBinLowEdge(low_bin_noise),proj->GetBinLowEdge(high_bin_noise+1));
	proj_noise->Fit("gaus","Q");
	TF1 *fit = proj_noise->GetFunction("gaus");

	noise_fit_pos->Fill(fit->GetParameter(1));
	noise_fit_pos_2d->Fill(i%NC,i/NC,fit->GetParameter(1));
	noise_fit_poserr->Fill(fit->GetParError(1));
	noise_fit_poserr_2d->Fill(i%NC,i/NC,fit->GetParError(1));

	// peak
	TH1D *proj_peak = dynamic_cast<TH1D*>(proj->Rebin(4,"proj_peak"));
	proj_peak->SetStats(kTRUE);
	proj_peak->GetXaxis()->SetRangeUser(proj->GetBinLowEdge(low_bin_peak),proj->GetBinLowEdge(high_bin_peak+1));

	Double_t mypar[8];
	mypar[0] = 0.0;
	mypar[1] = 0.0;
	mypar[2] = proj_peak->GetBinCenter(proj_peak->GetMaximumBin());
	if (gain_str == "G0") {
	  mypar[3] = 16.;
	} else if (gain_str == "HG0") {
	  mypar[3] = 29.;
	}
	mypar[4] = proj_peak->GetBinContent(proj_peak->GetMaximumBin());
	if (gain_str == "G0") {
	  mypar[5] = 0.17;
	} else if (gain_str == "HG0") {
	  mypar[5] = 0.14;
	}
	mypar[6] = 1.12;
	if (gain_str == "G0") {
	  mypar[7] = 0.12;
	} else if (gain_str == "HG0") {
	  mypar[7] = 0.14;
	}

	Double_t emypar[8];
	energyCalibration *thiscalibration = new energyCalibration();
	thiscalibration->setScanSign(1);
	thiscalibration->setStartParametersKb(mypar);
	thiscalibration->fixParameter(0,0.); // no background
	thiscalibration->fixParameter(1,0.);
	TF1* fittedfun = thiscalibration->fitSpectrumKb( proj_peak, mypar, emypar );

	fit_par3->Fill( mypar[3] );
	fit_par4->Fill( mypar[4] );
	fit_par5->Fill( mypar[5] );
	fit_par6->Fill( mypar[6] );
	fit_par7->Fill( mypar[7] );

	fit_par3_2d->Fill( i%NC, i/NC, mypar[3] );
	fit_par4_2d->Fill( i%NC, i/NC, mypar[4] );
	fit_par5_2d->Fill( i%NC, i/NC, mypar[5] );
	fit_par6_2d->Fill( i%NC, i/NC, mypar[6] );
	fit_par7_2d->Fill( i%NC, i/NC, mypar[7] );

	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 >= 58000 && i < 58000+10 ) || // bulk
	     ( i >= 10 && i < 10+10 ) || // edge
	     ( i >= 1024+10 && i < 1024+10+10 ) || // inner edge
	     ( i >= ( 256*1024 )+10 && i < ( 256*1024 )+10+10 ) || // double
	     ( i >= ( 257*1024 )+10 && i < ( 257*1024 )+10+10 ) || // next to double
	     ( i == ( 255*1024 )+255 ) // quad
	     ) {

	  string pixel_type = "x";
	  if ( i >= 58000 && i < 58000+10 ) {
	    pixel_type = "b";
	  } else if ( i >= 10 && i < 10+10 ) {
	    pixel_type = "e";
	  } else if ( i >= 1024+10 && i < 1024+10+10 ) {
	    pixel_type = "ie";
	  } else if ( i >= ( 256*1024 )+10 && i < ( 256*1024 )+10+10 ) {
	    pixel_type = "d";
	  } else if ( i >= ( 257*1024 )+10 && i < ( 257*1024 )+10+10 ) {
	    pixel_type = "ntd";
	  } else if ( i == ( 255*1024 )+255 ) {
	    pixel_type = "q";
	  }
	    
	  proj_noise->Draw();
	  c1->Update();
	  proj_noise->GetXaxis()->SetTitle( "Pedestal corrected ADC [ADU]" );
	  proj_noise->GetXaxis()->SetRangeUser( -100, 150 );
	  fit->SetParNames( "N_{#gamma}", "Peak pos", "Noise RMS" );
	  TPaveStats *st0 = (TPaveStats*)proj_noise->FindObject( "stats" );
	  st0->SetX1NDC(0.53);
	  st0->SetX2NDC(0.94);
	  st0->SetY1NDC(0.75);
	  st0->SetY2NDC(0.94);
	  st0->SetBorderSize(0);
	  st0->SetTextSize(0.04);
	  sprintf( savename, "%s/noise_%s_sc%d_%d_%s_M%s.png", plotfolder, pixel_type.c_str(), this_storagecell, i, gain_str.c_str(), module_str.c_str() );
	  c1->SaveAs( (const char *)(savename) );
	    
	  TF1 *gaus_Ka = new TF1( "gaus_Ka", "gaus", proj->GetBinLowEdge(low_bin_peak), proj->GetBinLowEdge(high_bin_peak+1) );
	  gaus_Ka->SetParameters( mypar[4], mypar[2], mypar[3] );
	  gaus_Ka->SetLineColor( kBlue );
	    
	  TF1 *erfc_Ka = new TF1( "erfc_Ka", "[0]/2.*(TMath::Erfc(([1]*(x-[2])/[3])/(TMath::Sqrt(2.))))", proj->GetBinLowEdge(low_bin_peak), proj->GetBinLowEdge(high_bin_peak+1) );
	  erfc_Ka->SetParameters( mypar[4]*mypar[5], 1, mypar[2], mypar[3] );
	  erfc_Ka->SetLineColor( kOrange );
	  
	  TF1 *gaus_Kb = new TF1( "gaus_Kb", "gaus", proj->GetBinLowEdge(low_bin_peak), proj->GetBinLowEdge(high_bin_peak+1) );
	  gaus_Kb->SetParameters( mypar[4]*mypar[7], mypar[6]*mypar[2], mypar[3] );
	  gaus_Kb->SetLineColor( kGreen+2 );
	    
	  TF1 *erfc_Kb = new TF1( "erfc_Kb", "[0]/2.*(TMath::Erfc(([1]*(x-[2])/[3])/(TMath::Sqrt(2.))))", proj->GetBinLowEdge(low_bin_peak), proj->GetBinLowEdge(high_bin_peak+1) );
	  erfc_Kb->SetParameters( mypar[4]*mypar[7]*mypar[5], 1, mypar[6]*mypar[2], mypar[3] );
	  erfc_Kb->SetLineColor( kOrange+7 );
	    
	  proj_peak->Draw();
	  erfc_Kb->Draw("same");
	  erfc_Ka->Draw("same");
	  gaus_Kb->Draw("same");
	  gaus_Ka->Draw("same");
	  fittedfun->Draw("same");
	  c1->Update();
	  proj_peak->GetXaxis()->SetTitle( "Pedestal corrected ADC [ADU]" );
	  fittedfun->SetParNames( "Bkg height", "Bkg grad", "K_{#alpha} pos", "Noise RMS", "K_{#alpha} height", "CS", "K_{#beta}/K_{#alpha} pos", "K_{#beta} frac" );
	  TPaveStats *st = (TPaveStats*)proj_peak->FindObject( "stats" );
	  st->SetX1NDC(0.22);
	  st->SetX2NDC(0.62);
	  st->SetY1NDC(0.7);
	  st->SetY2NDC(0.94);
	  st->SetBorderSize(0);
	  st->SetTextSize(0.04);
	  sprintf( savename, "%s/peak_%s_sc%d_%d_%s_M%s.png", plotfolder, pixel_type.c_str(), this_storagecell, i, gain_str.c_str(), module_str.c_str() );
	  c1->SaveAs( (const char *)(savename) );
	}
	  
	// gain
	gain_fit->Fill( mypar[2] - fit->GetParameter(1) );
	gain_fiterr->Fill( sqrt( pow(emypar[2],2) + pow(fit->GetParError(1),2) ) );
	gain_fit_2d->Fill( i%NC, i/NC, mypar[2] - fit->GetParameter(1) );
	gain_fiterr_2d->Fill( i%NC, i/NC, sqrt( pow(emypar[2],2) + pow(fit->GetParError(1),2) ) );
	gain_ADUper1keV_2d->Fill( i%NC, i/NC, ( mypar[2] - fit->GetParameter(1) ) / 8.0 );
	gainerr_ADUper1keV_2d->Fill( i%NC, i/NC, sqrt( pow( emypar[2],2 ) + pow( fit->GetParError(1),2) ) / 8.0 );
	  
	if (isEdge(i)) {
	  gain_fit_isEdge->Fill( mypar[2] - fit->GetParameter(1) );
	}
	if (isInnerEdge(i)) {
	  gain_fit_isInnerEdge->Fill( mypar[2] - fit->GetParameter(1) );
	}
	if (isDouble(i)) {
	  gain_fit_isDouble->Fill( mypar[2] - fit->GetParameter(1) );
	}
	if (isNextToDouble(i)) {
	    gain_fit_isNextToDouble->Fill( mypar[2] - fit->GetParameter(1) );
	}
	if (isQuad(i)) {
	  gain_fit_isQuad->Fill( mypar[2] - fit->GetParameter(1) );
	}
	if (isBulk(i)) {
	  gain_fit_isBulk->Fill( mypar[2] - fit->GetParameter(1) );
	}
	  
	delete thiscalibration;
	delete proj_noise;
	delete proj_peak;
	delete proj;
	  
      } else {
	pixel_mask[i] = false;
      } //end if can be fitted

    } // end for loop pixel numbers (i)
    
  } // end for loop sensor slices (j)

  sprintf( savename, "%s/pixelmask_afterfit_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() );
  pixelMaskObject->plotPixelMask( pixel_mask, savename );

  TPaveText *pave = new TPaveText( 0.86, 0.95, 0.91, 0.98, "blNDC" );
  pave->SetBorderSize(0);
  pave->SetFillStyle(0);
  pave->SetTextSize(0.06);
  pave->SetTextAlign(32);
  
  mapcanvas->cd();

  fit_par3_2d->GetXaxis()->SetTitle("Column");
  fit_par3_2d->GetYaxis()->SetTitle("Row");
  fit_par3_2d->GetYaxis()->SetTitleOffset(0.7);
  fit_par3_2d->Draw("colz");
  fit_par3_2d->GetZaxis()->SetRangeUser(0,50);
  sprintf( savename, "%s/fit_par3_2d_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() );
  mapcanvas->SaveAs( (const char *)(savename) );

  fit_par3_2d->Delete();

  fit_par4_2d->GetXaxis()->SetTitle("Column");
  fit_par4_2d->GetYaxis()->SetTitle("Row");
  fit_par4_2d->GetYaxis()->SetTitleOffset(0.7);
  fit_par4_2d->Draw("colz");
  fit_par4_2d->GetZaxis()->SetRangeUser(0,500);
  sprintf( savename, "%s/fit_par4_2d_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() );
  mapcanvas->SaveAs( (const char *)(savename) );

  fit_par4_2d->Delete();

  fit_par5_2d->GetXaxis()->SetTitle("Column");
  fit_par5_2d->GetYaxis()->SetTitle("Row");
  fit_par5_2d->GetYaxis()->SetTitleOffset(0.7);
  fit_par5_2d->Draw("colz");
  fit_par5_2d->GetZaxis()->SetRangeUser(0,0.5);
  sprintf( savename, "%s/fit_par5_2d_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() );
  mapcanvas->SaveAs( (const char *)(savename) );

  fit_par5_2d->Delete();

  fit_par6_2d->GetXaxis()->SetTitle("Column");
  fit_par6_2d->GetYaxis()->SetTitle("Row");
  fit_par6_2d->GetYaxis()->SetTitleOffset(0.7);
  fit_par6_2d->Draw("colz");
  fit_par6_2d->GetZaxis()->SetRangeUser(1.0,1.3); //was (1.0,1.25)
  sprintf( savename, "%s/fit_par6_2d_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() );
  mapcanvas->SaveAs( (const char *)(savename) );

  fit_par6_2d->Delete();

  fit_par7_2d->GetXaxis()->SetTitle("Column");
  fit_par7_2d->GetYaxis()->SetTitle("Row");
  fit_par7_2d->GetYaxis()->SetTitleOffset(0.7);
  fit_par7_2d->Draw("colz");
  fit_par7_2d->GetZaxis()->SetRangeUser(0.,0.4);
  sprintf( savename, "%s/fit_par7_2d_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() );
  mapcanvas->SaveAs( (const char *)(savename) );

  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( low_ADU_peak, high_ADU_peak );
  sprintf( savename, "%s/peak_fit_pos_2d_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() );
  mapcanvas->SaveAs( (const char *)(savename) );

  peak_fit_pos_2d->Delete();

  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,2);
  sprintf( savename, "%s/peak_fit_poserr_2d_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() );
  mapcanvas->SaveAs( (const char *)(savename) );

  peak_fit_poserr_2d->Delete();

  noise_fit_pos_2d->GetXaxis()->SetTitle("Column");
  noise_fit_pos_2d->GetYaxis()->SetTitle("Row");
  noise_fit_pos_2d->GetYaxis()->SetTitleOffset(0.7);
  noise_fit_pos_2d->Draw("colz");
  noise_fit_pos_2d->GetZaxis()->SetRangeUser(-5,5);
  sprintf( savename, "%s/noise_fit_pos_2d_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() );
  mapcanvas->SaveAs( (const char *)(savename) );

  noise_fit_pos_2d->Delete();

  noise_fit_poserr_2d->GetXaxis()->SetTitle("Column");
  noise_fit_poserr_2d->GetYaxis()->SetTitle("Row");
  noise_fit_poserr_2d->GetYaxis()->SetTitleOffset(0.7);
  noise_fit_poserr_2d->Draw("colz");
  if (gain_str == "HG0") {
    noise_fit_poserr_2d->GetZaxis()->SetRangeUser(0,0.1);
  } else if (gain_str == "G0") {
    noise_fit_poserr_2d->GetZaxis()->SetRangeUser(0,0.1); // was 0,0.05
  }
  sprintf( savename, "%s/noise_fit_poserr_2d_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() );
  mapcanvas->SaveAs( (const char *)(savename) );

  noise_fit_poserr_2d->Delete();

  gain_fit_2d->GetXaxis()->SetTitle("Column");
  gain_fit_2d->GetYaxis()->SetTitle("Row");
  gain_fit_2d->GetYaxis()->SetTitleOffset(0.7);
  gain_fit_2d->Draw("colz");
  sprintf( savename, "%s [ADU/8 keV]", gain_str.c_str() );
  pave->AddText( (const char *)(savename) );
  pave->Draw();
  gain_fit_2d->GetZaxis()->SetRangeUser( low_ADU_peak, high_ADU_peak );
  sprintf( savename, "%s/gain_fit_2d_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() );
  mapcanvas->SaveAs( (const char *)(savename) );
    
  gain_fit_2d->Delete();

  gain_fiterr_2d->GetXaxis()->SetTitle("Column");
  gain_fiterr_2d->GetYaxis()->SetTitle("Row");
  gain_fiterr_2d->GetYaxis()->SetTitleOffset(0.7);
  gain_fiterr_2d->Draw("colz");
  gain_fiterr_2d->GetZaxis()->SetRangeUser(0,2);
  sprintf( savename, "%s/gain_fiterr_2d_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() );
  mapcanvas->SaveAs( (const char *)(savename) );

  gain_fiterr_2d->Delete();

  gain_ADUper1keV_2d->GetXaxis()->SetTitle("Column");
  gain_ADUper1keV_2d->GetYaxis()->SetTitle("Row");
  gain_ADUper1keV_2d->GetYaxis()->SetTitleOffset(0.7);
  gain_ADUper1keV_2d->Draw("colz");
  if (gain_str == "HG0") {
    gain_ADUper1keV_2d->GetZaxis()->SetRangeUser(80,120);
  } else if (gain_str == "G0") {
    gain_ADUper1keV_2d->GetZaxis()->SetRangeUser(25,45); //was 35,50
  }
  sprintf( savename, "%s/gain_ADUper1keV_2d_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() );
  mapcanvas->SaveAs( (const char *)(savename) );
  
  //gain_ADUper1keV_2d->Delete();

  gainerr_ADUper1keV_2d->GetXaxis()->SetTitle("Column");
  gainerr_ADUper1keV_2d->GetYaxis()->SetTitle("Row");
  gainerr_ADUper1keV_2d->GetYaxis()->SetTitleOffset(0.7);
  gainerr_ADUper1keV_2d->Draw("colz");
  if (gain_str == "HG0") {
    gainerr_ADUper1keV_2d->GetZaxis()->SetRangeUser(0,0.5);
  } else if (gain_str == "G0") {
    gainerr_ADUper1keV_2d->GetZaxis()->SetRangeUser(0,0.25);
  }
  sprintf( savename, "%s/gainerr_ADUper1keV_2d_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() );
  mapcanvas->SaveAs( (const char *)(savename) );

  //gainerr_ADUper1keV_2d[this_storagecell]->Delete();
    
  c1->cd();
  
  fit_par3->GetXaxis()->SetTitle("Fit par 3");
  fit_par3->Draw();
  sprintf( savename, "%s/fit_par3_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() );
  c1->SaveAs( (const char *)(savename) );
    
  fit_par3->Delete();

  fit_par4->GetXaxis()->SetTitle("Fit par 4");
  fit_par4->Draw();
  sprintf( savename, "%s/fit_par4_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() );
  c1->SaveAs( (const char *)(savename) );

  fit_par4->Delete();

  fit_par5->GetXaxis()->SetTitle("Fit par 5");
  fit_par5->Draw();
  sprintf( savename, "%s/fit_par5_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() );
  c1->SaveAs( (const char *)(savename) );

  fit_par5->Delete();

  fit_par6->GetXaxis()->SetTitle("Fit par 6");
  fit_par6->Draw();
  sprintf( savename, "%s/fit_par6_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() );
  c1->SaveAs( (const char *)(savename) );

  fit_par6->Delete();

  fit_par7->GetXaxis()->SetTitle("Fit par 7");
  fit_par7->Draw();
  sprintf( savename, "%s/fit_par7_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() );
  c1->SaveAs( (const char *)(savename) );

  peak_fit_pos->GetXaxis()->SetTitle("Peak position [ADU]");
  peak_fit_pos->Draw();
  sprintf( savename, "%s/peak_fit_pos_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() );
  c1->SaveAs( (const char *)(savename) );

  peak_fit_pos->Delete();

  peak_fit_poserr->GetXaxis()->SetTitle("Peak position uncert [ADU]");
  peak_fit_poserr->Draw();
  sprintf( savename, "%s/peak_fit_poserr_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() );
  c1->SaveAs( (const char *)(savename) );

  peak_fit_poserr->Delete();

  noise_fit_pos->GetXaxis()->SetTitle("Noise position [ADU]");
  noise_fit_pos->Draw();
  sprintf( savename, "%s/noise_fit_pos_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() );
  c1->SaveAs( (const char *)(savename) );

  noise_fit_pos->Delete();

  noise_fit_poserr->GetXaxis()->SetTitle("Noise position uncert [ADU]");
  noise_fit_poserr->Draw();
  sprintf( savename, "%s/noise_fit_poserr_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() );
  c1->SaveAs( (const char *)(savename) );

  noise_fit_poserr->Delete();

  sprintf( savename, "Gain %s [ADU / 8 keV]", gain_str.c_str() );
  gain_fit->GetXaxis()->SetTitle( (const char *)(savename) );
  gain_fit->Draw();
  sprintf( savename, "%s/gain_fit_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() );
  c1->SaveAs( (const char *)(savename) );

  //gain_fit->Delete();

  gain_fit->GetXaxis()->SetRangeUser( low_ADU_peak+30, high_ADU_peak );
  gain_fit->Fit("gaus");
  gain_fit->Draw();
  c1->Update();
  TPaveText *pave2 = new TPaveText( 0.6, 0.8, 0.94, 0.94, "blNDC" );
  pave2->SetBorderSize(0);
  pave2->SetFillStyle(0);
  pave2->SetTextSize(0.04);
  pave2->SetTextAlign(32);
  TF1* gain_fit_gaus = gain_fit->GetFunction("gaus");
  sprintf( savename, "Mean %0.2f #pm %0.2f", gain_fit_gaus->GetParameter(1), gain_fit_gaus->GetParError(1) );
  pave2->AddText( (const char *)(savename) );
  sprintf( savename, "Sigma %0.2f #pm %0.2f", gain_fit_gaus->GetParameter(2), gain_fit_gaus->GetParError(2) );
  pave2->AddText( (const char *)(savename) );
  pave2->Draw();
  gain_fit->SetStats(kFALSE);
  sprintf( savename, "%s/gain_fit_fit_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() );
  c1->SaveAs( (const char *)(savename) );
  
  gain_fit->Delete();

  gain_fiterr->GetXaxis()->SetTitle("Gain uncert [ADU / 8 keV]");
  gain_fiterr->Draw();
  sprintf( savename, "%s/gain_fiterr_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() );
  c1->SaveAs( (const char *)(savename) );

  gain_fiterr->Delete();

  gain_fit_isEdge->SetLineColor(kBlue);
  gain_fit_isInnerEdge->SetLineColor(kCyan);
  gain_fit_isDouble->SetLineColor(kGreen+2);
  gain_fit_isNextToDouble->SetLineColor(kRed);
  gain_fit_isQuad->SetLineColor(kOrange);

  gain_fit_isEdge->Scale( 1./gain_fit_isEdge->GetEntries() );
  gain_fit_isInnerEdge->Scale( 1./gain_fit_isInnerEdge->GetEntries() );
  gain_fit_isDouble->Scale( 1./gain_fit_isDouble->GetEntries() );
  gain_fit_isNextToDouble->Scale( 1./gain_fit_isNextToDouble->GetEntries() );
  gain_fit_isQuad->Scale( 1./gain_fit_isQuad->GetEntries() );
  gain_fit_isBulk->Scale( 1./gain_fit_isBulk->GetEntries() );

  TLegend *leg = new TLegend( 0.62, 0.6, 0.93, 0.93 );
  leg->AddEntry( gain_fit_isBulk, "Normal", "l" );
  leg->AddEntry( gain_fit_isDouble, "Double", "l" );
  leg->AddEntry( gain_fit_isNextToDouble, "Next to D", "l" );
  leg->AddEntry( gain_fit_isEdge, "Edge", "l" );
  leg->AddEntry( gain_fit_isInnerEdge, "Inner E", "l" );

  sprintf( savename, "Gain %s [ADU / 8 keV]", gain_str.c_str() );
  gain_fit_isDouble->GetXaxis()->SetTitle( (const char *)(savename) );
  gain_fit_isDouble->GetYaxis()->SetTitle("Normalised");
  gain_fit_isDouble->GetYaxis()->SetTitleOffset(1.3);
  gain_fit_isDouble->SetMinimum(0.0);
  gain_fit_isDouble->SetMaximum(0.16);
  gain_fit_isDouble->Draw();
  gain_fit_isEdge->Draw("same");
  gain_fit_isInnerEdge->Draw("same");
  gain_fit_isNextToDouble->Draw("same");
  gain_fit_isBulk->Draw("same");
  leg->Draw("same");
  sprintf( savename, "%s/gain_fit_perType_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() );
  c1->SaveAs( (const char *)(savename) );

  gain_fit_isDouble->Delete();
  gain_fit_isEdge->Delete();
  gain_fit_isInnerEdge->Delete();
  gain_fit_isNextToDouble->Delete();
  gain_fit_isBulk->Delete();

  sprintf( savename, "%s/CuFluo_gain_sc%d_%s_M%s.root", rootdatafolder, this_storagecell, gain_str.c_str(), module_str.c_str() );
  TFile* saved_file = new TFile( (const char *)(savename), "RECREATE" );
  gain_ADUper1keV_2d->Write();
  gainerr_ADUper1keV_2d->Write();
  saved_file->Close();

  gain_ADUper1keV_2d->Delete();
  gainerr_ADUper1keV_2d->Delete();

}