// file to calculate pedestal correction of fluo data for all storage cells, only G0
// make correction and save spectrum per pixel
// then fit fluo spectrum per pixel
// and save peak position and uncertainty
// changes by VH 210906: to eliminate hardcoded absolute paths, uses location of the analysis root files as additional input argument (accordingly changed in filename_creator.sh)

#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 "TObjArray.h"

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

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

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

  /*
  if (argc != 6) {
    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" << endl;
    cout << "arg 4: specify pede file" << endl;
    cout << "arg 5: specify data file" << endl;
    cout << " " << endl;
    exit(1);
  }
  */ //uncomment for SR

  if (argc != 7) {
    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 << " " << endl;
    exit(1);
  } //uncomment for VH 210906

  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]; //uncomment for VH 210906

  bool save_to_afs = false;

  int createHistoFile = 1;

  //char histoname[128];
  //char savename[128];
  char histoname[256]; // VH 210902
  char savename[256];  // VH 210902
  char plotfolder[256];
  char rootdatafolder[256];
  int filen = 352;
  int pedefilen = 5;

  // 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);
  // /mnt/pcmoench_jungfrau_data/jungfrau_ana_sophie/Mxxx_CalibAna
  //sprintf(savename,"/mnt/sls_det_storage/jungfrau_data1/jungfrau_ana_sophie/M%s_CalibAna", module_str.c_str()); //uncomment for SR
  sprintf(savename,"%s", anadata_loc.c_str()); //uncomment for VH 210906
  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();

    vector<jungfrauPedestal*> pedestalObject_SC; //( NSC, new jungfrauPedestal() ); //should call default constructor 16 times
    pedestalObject_SC.reserve(NSC);

    vector< vector<uint16_t> > pedestals16_G0_start( NSC, vector<uint16_t>(NCH) ); //I suppose, this is to track pedestal shifting over the course of data taking
    vector< vector<double> > pedeRMS16_G0( NSC, vector<double>(NCH) );

    { //scope to limit the following 2 vectors
      vector<TH2F*> pedestalsG0;
      pedestalsG0.reserve(NSC);
      vector<TH2F*> pedeRMSG0;
      pedeRMSG0.reserve(NSC);
      for (int sci = 0; sci < NSC; ++sci) {
	pedestalObject_SC.push_back( new jungfrauPedestal() );
	pedestalObject_SC.at(sci)->pedestalSetNFrames(100);
	
	Char_t *pedehistoname = new Char_t[50];
	snprintf( pedehistoname, 50, "pedestalsG0_sc%d", sci );
	pedestalsG0.push_back( new TH2F( pedehistoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 ) );
	snprintf( pedehistoname, 50, "pedeRMSG0_sc%d", sci );
	pedeRMSG0.push_back( new TH2F( pedehistoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 ) );
      }

      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() ); //VH: note, this adds a double slash in the filepath
	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
	  pedestalObject_SC.at( thisfile->currentSCnumber() )->addFrameToPedestalCalculation( thisfile->getFrameDataHandle() );
	}
	
	thisfile->close();
      
      } // end of loops over files
      
      for (int sci = 0; sci < NSC; ++sci) {
	for (int i = 0; i < NCH; ++i) {
	  if (pixel_mask[i] == true) {
	    pedestalsG0.at(sci)->Fill( i%NC, i/NC, pedestalObject_SC[sci]->pedestalOfChannel(i) );
	    pedeRMSG0.at(sci)->Fill( i%NC, i/NC, pedestalObject_SC[sci]->rmsOfChannel(i) );
	  }
	}
	pedestalObject_SC.at(sci)->pedestalData((uint16_t*)(pedestals16_G0_start.at(sci).data())); //this loads the pedestal data into the array pedestals16_G0_start
	pedestalObject_SC.at(sci)->pedestalRMSData(pedeRMS16_G0[sci].data()); //same here
	pedestalObject_SC.at(sci)->pedestalResetUpdates();
	//pedestalObject[sci]->pedestalClear(); //I don't need to clear if I only calculate pedestal once.
      }

      sprintf( savename, "%s/pixelmask_%s_M%s.png", plotfolder, gain_str.c_str(), module_str.c_str() );
      pixelMaskObject->plotPixelMask( pixel_mask, savename );
      cout << "after chip mask, n masked pixels is " << pixelMaskObject->getNMasked(pixel_mask) << endl;
      
      mapcanvas->cd();
      
      for (int sci = 0; sci < 16; ++sci) {
	
	pedestalsG0[sci]->GetXaxis()->SetTitle("Column");
	pedestalsG0[sci]->GetYaxis()->SetTitle("Row");
	pedestalsG0[sci]->GetYaxis()->SetTitleOffset(0.7);
	pedestalsG0[sci]->Draw("colz");
	sprintf(savename,"%s/pedeG0_sc%d_%s_M%s.png", plotfolder, sci, gain_str.c_str(), module_str.c_str());
	mapcanvas->SaveAs((const char *)(savename));

	pedestalsG0[sci]->Delete(); //not needed hereafter (try to save memory)
	
	pedeRMSG0[sci]->GetXaxis()->SetTitle("Column");
	pedeRMSG0[sci]->GetYaxis()->SetTitle("Row");
	pedeRMSG0[sci]->GetYaxis()->SetTitleOffset(0.7);
	pedeRMSG0[sci]->GetZaxis()->SetRangeUser(0,30);
	pedeRMSG0[sci]->Draw("colz");
	sprintf(savename,"%s/pedeRMSG0_sc%d_%s_M%s.png", plotfolder, sci, gain_str.c_str(), module_str.c_str());
	mapcanvas->SaveAs((const char *)(savename));

	pedeRMSG0[sci]->Delete(); //not needed hereafter (try to save memory)
	
      }

      //Make sure histograms are deleted

    } // end of scope to limit pedestalsG0 and pedeRMSG0
    
    int adc2d_nbin= 1200;
    if (isJF11) adc2d_nbin= 1600;
    
    vector<TH2I*> adc2d_1; //(8, vector <TH2I*>(16) ); //declare
    vector<TH2I*> adc2d_2;
    vector<TH2I*> adc2d_3;
    vector<TH2I*> adc2d_4;
    vector<TH2I*> adc2d_5;
    vector<TH2I*> adc2d_6;
    vector<TH2I*> adc2d_7;
    vector<TH2I*> adc2d_8;
    //TObjArray *adc2d_8 = new TObjArray(NSC);
    //adc2d.reserve(8);
    /*
      TH2I *adc2d_2[16]; = new TH2I("adc2d_2","",adc2d_nbin,-200-0.5,adc2d_nbin-200-0.5,65536,(65536*1-0.5),(65536*2-0.5));
      TH2I *adc2d_3[16]; = new TH2I("adc2d_3","",adc2d_nbin,-200-0.5,adc2d_nbin-200-0.5,65536,(65536*2-0.5),(65536*3-0.5));
      TH2I *adc2d_4[16]; = new TH2I("adc2d_4","",adc2d_nbin,-200-0.5,adc2d_nbin-200-0.5,65536,(65536*3-0.5),(65536*4-0.5));
      TH2I *adc2d_5[14]; = new TH2I("adc2d_5","",adc2d_nbin,-200-0.5,adc2d_nbin-200-0.5,65536,(65536*4-0.5),(65536*5-0.5));
      TH2I *adc2d_6[16]; = new TH2I("adc2d_6","",adc2d_nbin,-200-0.5,adc2d_nbin-200-0.5,65536,(65536*5-0.5),(65536*6-0.5));
      TH2I *adc2d_7[16]; = new TH2I("adc2d_7","",adc2d_nbin,-200-0.5,adc2d_nbin-200-0.5,65536,(65536*6-0.5),(65536*7-0.5));
      TH2I *adc2d_8[16]; = new TH2I("adc2d_8","",adc2d_nbin,-200-0.5,adc2d_nbin-200-0.5,65536,(65536*7-0.5),(65536*8-0.5));
    */
    
    cout << "Initialize adc2d vectors..." << endl;
    
    //initialize
    //for ( int adci = 0; adci < 8; ++adci ) {
    //vector<TH2I*> v;
    //v.reserve(NSC);
    for ( int sci = 0; sci < NSC; ++sci ) {
      Char_t *_histoname = new Char_t[50];
      snprintf( _histoname, 50, "adc2d_1_sc%d", sci );
      adc2d_1.push_back( new TH2I( _histoname, "", adc2d_nbin, -200-0.5, adc2d_nbin-200-0.5, 65536, ( 65536*0-0.5 ), ( 65536*( 0+1 )-0.5 ) ) );
      snprintf( _histoname, 50, "adc2d_2_sc%d", sci );
      adc2d_2.push_back( new TH2I( _histoname, "", adc2d_nbin, -200-0.5, adc2d_nbin-200-0.5, 65536, ( 65536*1-0.5 ), ( 65536*( 1+1 )-0.5 ) ) );
      snprintf( _histoname, 50, "adc2d_3_sc%d", sci );
      adc2d_3.push_back( new TH2I( _histoname, "", adc2d_nbin, -200-0.5, adc2d_nbin-200-0.5, 65536, ( 65536*2-0.5 ), ( 65536*( 2+1 )-0.5 ) ) );
      snprintf( _histoname, 50, "adc2d_4_sc%d", sci );
      adc2d_4.push_back( new TH2I( _histoname, "", adc2d_nbin, -200-0.5, adc2d_nbin-200-0.5, 65536, ( 65536*3-0.5 ), ( 65536*( 3+1 )-0.5 ) ) );
      snprintf( _histoname, 50, "adc2d_5_sc%d", sci );
      adc2d_5.push_back( new TH2I( _histoname, "", adc2d_nbin, -200-0.5, adc2d_nbin-200-0.5, 65536, ( 65536*4-0.5 ), ( 65536*( 4+1 )-0.5 ) ) );
      snprintf( _histoname, 50, "adc2d_6_sc%d", sci );
      adc2d_6.push_back( new TH2I( _histoname, "", adc2d_nbin, -200-0.5, adc2d_nbin-200-0.5, 65536, ( 65536*5-0.5 ), ( 65536*( 5+1 )-0.5 ) ) );
      snprintf( _histoname, 50, "adc2d_7_sc%d", sci );
      adc2d_7.push_back( new TH2I( _histoname, "", adc2d_nbin, -200-0.5, adc2d_nbin-200-0.5, 65536, ( 65536*6-0.5 ), ( 65536*( 6+1 )-0.5 ) ) );	
      snprintf( _histoname, 50, "adc2d_8_sc%d", sci ); 
      adc2d_8.push_back( new TH2I( _histoname, "", adc2d_nbin, -200-0.5, adc2d_nbin-200-0.5, 65536, ( 65536*7-0.5 ), ( 65536*( 7+1 )-0.5 ) ) );
    }
      //adc2d.push_back(v);
      //}
    
    cout << "Done." << endl;
    
    //declare
    vector<TH1D*> adcpc_spec;
    adcpc_spec.reserve(NSC);
    vector<TH2F*> pede_updates;
    pede_updates.reserve(NSC);
    vector<TH2F*> pede_diff;
    pede_diff.reserve(NSC);

    cout << "Initialize vectors of histograms..." << endl;
    
    //initialize
    for ( int sci = 0; sci < NSC; ++sci ) {
      Char_t *_histoname = new Char_t[50];
      snprintf( _histoname, 50, "adcpc_spec_sc%d", sci );
      adcpc_spec.push_back( new TH1D( _histoname,"",300,0,3000 ) ); //spectrum for every sc, every 10000 frames
      snprintf( _histoname, 50, "pede_updates_sc%d", sci );
      pede_updates.push_back( new TH2F( _histoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 ) );
      snprintf( _histoname, 50, "pede_diff_sc%d", sci );
      pede_diff.push_back( new TH2F( _histoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 ) );
    }

    cout << "Done." << endl;

    for (int filei = 0; filei < filen; filei++) {

      // open data file
      sprintf(savename,"%s/%s_%%6.6d.dat", data_loc.c_str(), data_file.c_str()); //VH: note, this adds a double slash in the filepath
      thisfile->open((char*)savename, filei);
      
      while (thisfile->readNextFrame()) {
	
	uint16_t* imagedptr = thisfile->getFrameDataHandle();
	uint64_t scnumber = thisfile->currentSCnumber();
	
	pedestalObject_SC.at(scnumber)->addG0FrameToPedestalCalculationWThreshold( imagedptr, pedestalObject_SC.at(scnumber), pedeRMS16_G0.at(scnumber).data() );

	for (int i = 0; i < NCH; i++) {
	    
	  uint16_t gain = (imagedptr[i]&0xc000) >> 14;
	    
	  if (gain == 0) {
	      
            int adcpc = (imagedptr[i]&0x3fff) - pedestalObject_SC.at(scnumber)->pedestalOfChannel(i);
	    adcpc_spec.at(scnumber)->Fill(adcpc);
	      
	    if (i < (65536*1)) {
	      adc2d_1[scnumber]->Fill(adcpc,i);
	    } else if (i < (65536*2)) {
	      adc2d_2[scnumber]->Fill(adcpc,i);
	    } else if (i < (65536*3)) {
	      adc2d_3[scnumber]->Fill(adcpc,i);
	    } else if (i < (65536*4)) {
	      adc2d_4[scnumber]->Fill(adcpc,i);
	    } else if (i < (65536*5)) {
	      adc2d_5[scnumber]->Fill(adcpc,i);
	    } else if (i < (65536*6)) {
	      adc2d_6[scnumber]->Fill(adcpc,i);
	    } else if (i < (65536*7)) {
	      adc2d_7[scnumber]->Fill(adcpc,i);
	    } else if (i < (65536*8)) {
    //TH2I *adc2d_now = dynamic_cast<TH2I*>(adc2d_8->At(scnumber));
              adc2d_8[scnumber]->Fill(adcpc,i);
	    }
	  }
	}
      }
      thisfile->close();

      if ( filei%16 == 0 ) { // every 16th file equals every 10000 frames for each storage cell
      
	for ( int sci = 0; sci < 16; ++sci ) {

	  adcpc_spec[sci]->GetXaxis()->SetTitle("Pedestal corrected ADC [ADU]");
	  adcpc_spec[sci]->Draw();
	  mapcanvas->SetLogy();
	  sprintf(savename,"%s/adcpc_spec_sc%d_slice%d_%s_M%s.png", plotfolder, sci, ( filei/16 ), gain_str.c_str(), module_str.c_str());
	  mapcanvas->SaveAs((const char *)(savename));
	  mapcanvas->SetLogy(0);
	  adcpc_spec[sci]->Reset();

	  pede_updates[sci]->Reset();
	  pede_diff[sci]->Reset();
	  for (int i = 0; i < NCH; i++) {
	    pede_updates[sci]->Fill( i%NC, i/NC, pedestalObject_SC[sci]->pedestalUpdates(i) );
	    pede_diff[sci]->Fill( i%NC, i/NC, pedestalObject_SC[sci]->pedestalOfChannel(i) - pedestals16_G0_start[sci][i] );
	    pedestals16_G0_start[sci][i] = pedestalObject_SC[sci]->pedestalOfChannel(i);
	  }

	  pede_updates[sci]->GetXaxis()->SetTitle("Column");
	  pede_updates[sci]->GetYaxis()->SetTitle("Row");
	  pede_updates[sci]->GetYaxis()->SetTitleOffset(0.7);
	  pede_updates[sci]->GetZaxis()->SetRangeUser(0,10000);
	  pede_updates[sci]->Draw("colz");
	  sprintf(savename,"%s/pede_updates_sc%d_slice%d_%s_M%s.png", plotfolder, sci, ( filei/16 ), gain_str.c_str(), module_str.c_str());
	  mapcanvas->SaveAs((const char *)(savename));

	  pede_diff[sci]->GetXaxis()->SetTitle("Column");
	  pede_diff[sci]->GetYaxis()->SetTitle("Row");
	  pede_diff[sci]->GetYaxis()->SetTitleOffset(0.7);
	  pede_diff[sci]->GetZaxis()->SetRangeUser(-40,40);
	  pede_diff[sci]->Draw("colz");
	  sprintf(savename,"%s/pede_diff_sc%d_slice%d_%s_M%s.png", plotfolder, sci, ( filei/16 ), gain_str.c_str(), module_str.c_str());
	  mapcanvas->SaveAs((const char *)(savename));
	  
	  pedestalObject_SC[sci]->pedestalResetUpdates();
	
	} //end of storage cell loop
      
      } //end of if filei%16

    } // end of file loop

    //sprintf(savename,"/mnt/sls_det_storage/jungfrau_data1/jungfrau_ana_sophie/M%s_CalibAna/CuFluo_%s_file0to%d.root", module_str.c_str(), gain_str.c_str(), filen-1);
    sprintf(savename,"%s/CuFluo_%s_file0to%d.root", anadata_loc.c_str(), gain_str.c_str(), filen-1); //uncomment for VH 210906
    TFile* saved_file = new TFile((const char *)(savename),"RECREATE");
    //for ( int adci = 0; adci < 8; ++adci ) {
      for ( int sci = 0; sci < 16; ++sci ) {
	adc2d_1[sci]->Write();
	adc2d_2[sci]->Write();
	adc2d_3[sci]->Write();
	adc2d_4[sci]->Write();
	adc2d_5[sci]->Write();
	adc2d_6[sci]->Write();
	adc2d_7[sci]->Write();
	adc2d_8[sci]->Write();
	//TH2I *adc2d_now = dynamic_cast<TH2I*>(adc2d_8->At(sci));

	//All pede updates have been written, so they can be destroyed.
	pede_updates[sci]->Delete();
	pede_diff[sci]->Delete();

	//In theory, I believe, adc2d can now also be destroyed since it has already been written to file
	adc2d_1[sci]->Delete();
	adc2d_2[sci]->Delete();
	adc2d_3[sci]->Delete();
	adc2d_4[sci]->Delete();
	adc2d_5[sci]->Delete();
	adc2d_6[sci]->Delete();
	adc2d_7[sci]->Delete();
	adc2d_8[sci]->Delete();
      }
      //}
    saved_file->Close();
    
  } // end if i creatHistoFile

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

  c1->cd();

  //sprintf(savename,"/mnt/sls_det_storage/jungfrau_data1/jungfrau_ana_sophie/M%s_CalibAna/CuFluo_%s_file0to%d.root", module_str.c_str(), gain_str.c_str(), filen-1);
  sprintf(savename,"%s/CuFluo_%s_file0to%d.root", anadata_loc.c_str(), gain_str.c_str(), 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
  vector<TH1F*> fit_par3;
  fit_par3.reserve(NSC);
  vector<TH1F*> fit_par4;
  fit_par4.reserve(NSC);
  vector<TH1F*> fit_par5;
  fit_par5.reserve(NSC);
  vector<TH1F*> fit_par6;
  fit_par6.reserve(NSC);
  vector<TH1F*> fit_par7;
  fit_par7.reserve(NSC);

  vector<TH2F*> fit_par3_2d;
  fit_par3_2d.reserve(NSC);
  vector<TH2F*> fit_par4_2d;
  fit_par4_2d.reserve(NSC);
  vector<TH2F*> fit_par5_2d;
  fit_par5_2d.reserve(NSC);
  vector<TH2F*> fit_par6_2d;
  fit_par6_2d.reserve(NSC);
  vector<TH2F*> fit_par7_2d;
  fit_par7_2d.reserve(NSC);

  vector<TH1F*> peak_fit_pos;
  peak_fit_pos.reserve(NSC);
  vector<TH1F*> peak_fit_poserr;
  peak_fit_poserr.reserve(NSC);
  vector<TH2F*> peak_fit_pos_2d;
  peak_fit_pos_2d.reserve(NSC);
  vector<TH2F*> peak_fit_poserr_2d;
  peak_fit_poserr_2d.reserve(NSC);

  vector<TH1F*> noise_fit_pos;
  noise_fit_pos.reserve(NSC);
  vector<TH1F*> noise_fit_poserr;
  noise_fit_poserr.reserve(NSC);
  vector<TH2F*> noise_fit_pos_2d;
  noise_fit_pos_2d.reserve(NSC);
  vector<TH2F*> noise_fit_poserr_2d;
  noise_fit_poserr_2d.reserve(NSC);

  vector<TH1F*> gain_fit;
  gain_fit.reserve(NSC);
  vector<TH1F*> gain_fiterr;
  gain_fiterr.reserve(NSC);
  vector<TH2F*> gain_fit_2d;
  gain_fit_2d.reserve(NSC);
  vector<TH2F*> gain_fiterr_2d;
  gain_fiterr_2d.reserve(NSC);

  vector<TH1F*> gain_fit_isEdge;
  gain_fit_isEdge.reserve(NSC);
  vector<TH1F*> gain_fit_isInnerEdge;
  gain_fit_isInnerEdge.reserve(NSC);
  vector<TH1F*> gain_fit_isDouble;
  gain_fit_isDouble.reserve(NSC);
  vector<TH1F*> gain_fit_isNextToDouble;
  gain_fit_isNextToDouble.reserve(NSC); 
  vector<TH1F*> gain_fit_isQuad;
  gain_fit_isQuad.reserve(NSC);
  vector<TH1F*> gain_fit_isBulk;
  gain_fit_isBulk.reserve(NSC);

  vector<TH2F*> gain_ADUper1keV_2d;
  gain_ADUper1keV_2d.reserve(NSC);
  vector<TH2F*> gainerr_ADUper1keV_2d;
  gainerr_ADUper1keV_2d.reserve(NSC);

  //initialize
  for ( int sci = 0; sci < NSC; ++sci ) {

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

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

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

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

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

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

  } //end histo creating for loop

  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 sci = 0; sci < NSC; ++sci ) {

    for ( int j = 1; j < 9; ++j ) {
      cout << "SC = " << sci << " slice " << j << endl;
      sprintf( histoname, "adc2d_%d_sc%d", j, sci );
      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[sci]->Fill(fit->GetParameter(1));
	  noise_fit_pos_2d[sci]->Fill(i%NC,i/NC,fit->GetParameter(1));
	  noise_fit_poserr[sci]->Fill(fit->GetParError(1));
	  noise_fit_poserr_2d[sci]->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[sci]->Fill( mypar[3] );
	  fit_par4[sci]->Fill( mypar[4] );
	  fit_par5[sci]->Fill( mypar[5] );
	  fit_par6[sci]->Fill( mypar[6] );
	  fit_par7[sci]->Fill( mypar[7] );

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

	  peak_fit_pos[sci]->Fill( mypar[2] );
	  peak_fit_poserr[sci]->Fill( emypar[2] );
	  peak_fit_pos_2d[sci]->Fill( i%NC, i/NC, mypar[2] );
	  peak_fit_poserr_2d[sci]->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(), sci, 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(), sci, i, gain_str.c_str(), module_str.c_str() );
	    c1->SaveAs( (const char *)(savename) );
	  }
	  
	  // gain
	  gain_fit[sci]->Fill( mypar[2] - fit->GetParameter(1) );
	  gain_fiterr[sci]->Fill( sqrt( pow(emypar[2],2) + pow(fit->GetParError(1),2) ) );
	  gain_fit_2d[sci]->Fill( i%NC, i/NC, mypar[2] - fit->GetParameter(1) );
	  gain_fiterr_2d[sci]->Fill( i%NC, i/NC, sqrt( pow(emypar[2],2) + pow(fit->GetParError(1),2) ) );
	  gain_ADUper1keV_2d[sci]->Fill( i%NC, i/NC, ( mypar[2] - fit->GetParameter(1) ) / 8.0 );
	  gainerr_ADUper1keV_2d[sci]->Fill( i%NC, i/NC, sqrt( pow( emypar[2],2 ) + pow( fit->GetParError(1),2) ) / 8.0 );
	  
	  if (isEdge(i)) {
	    gain_fit_isEdge[sci]->Fill( mypar[2] - fit->GetParameter(1) );
	  }
	  if (isInnerEdge(i)) {
	    gain_fit_isInnerEdge[sci]->Fill( mypar[2] - fit->GetParameter(1) );
	  }
	  if (isDouble(i)) {
	    gain_fit_isDouble[sci]->Fill( mypar[2] - fit->GetParameter(1) );
	  }
	  if (isNextToDouble(i)) {
	    gain_fit_isNextToDouble[sci]->Fill( mypar[2] - fit->GetParameter(1) );
	  }
	  if (isQuad(i)) {
	    gain_fit_isQuad[sci]->Fill( mypar[2] - fit->GetParameter(1) );
	  }
	  if (isBulk(i)) {
	    gain_fit_isBulk[sci]->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, sci, 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[sci]->GetXaxis()->SetTitle("Column");
    fit_par3_2d[sci]->GetYaxis()->SetTitle("Row");
    fit_par3_2d[sci]->GetYaxis()->SetTitleOffset(0.7);
    fit_par3_2d[sci]->Draw("colz");
    fit_par3_2d[sci]->GetZaxis()->SetRangeUser(0,50);
    sprintf( savename, "%s/fit_par3_2d_sc%d_%s_M%s.png", plotfolder, sci, gain_str.c_str(), module_str.c_str() );
    mapcanvas->SaveAs( (const char *)(savename) );

    fit_par3_2d[sci]->Delete();

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

    fit_par4_2d[sci]->Delete();

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

    fit_par5_2d[sci]->Delete();

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

    fit_par6_2d[sci]->Delete();

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

    peak_fit_pos_2d[sci]->GetXaxis()->SetTitle("Column");
    peak_fit_pos_2d[sci]->GetYaxis()->SetTitle("Row");
    peak_fit_pos_2d[sci]->GetYaxis()->SetTitleOffset(0.7);
    peak_fit_pos_2d[sci]->Draw("colz");
    peak_fit_pos_2d[sci]->GetZaxis()->SetRangeUser( low_ADU_peak, high_ADU_peak );
    sprintf( savename, "%s/peak_fit_pos_2d_sc%d_%s_M%s.png", plotfolder, sci, gain_str.c_str(), module_str.c_str() );
    mapcanvas->SaveAs( (const char *)(savename) );

    peak_fit_pos_2d[sci]->Delete();

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

    peak_fit_poserr_2d[sci]->Delete();

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

    noise_fit_pos_2d[sci]->Delete();

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

    noise_fit_poserr_2d[sci]->Delete();

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

    gain_fit_2d[sci]->Delete();

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

    gain_fiterr_2d[sci]->Delete();

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

    //gain_ADUper1keV_2d[sci]->Delete();

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

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

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

    fit_par4[sci]->Delete();

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

    fit_par5[sci]->Delete();

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

    fit_par6[sci]->Delete();

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

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

    peak_fit_pos[sci]->Delete();

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

    peak_fit_poserr[sci]->Delete();

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

    noise_fit_pos[sci]->Delete();

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

    noise_fit_poserr[sci]->Delete();

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

    //gain_fit[sci]->Delete();

    gain_fit[sci]->GetXaxis()->SetRangeUser( low_ADU_peak+30, high_ADU_peak );
    gain_fit[sci]->Fit("gaus");
    gain_fit[sci]->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[sci]->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[sci]->SetStats(kFALSE);
    sprintf( savename, "%s/gain_fit_fit_sc%d_%s_M%s.png", plotfolder, sci, gain_str.c_str(), module_str.c_str() );
    c1->SaveAs( (const char *)(savename) );

    gain_fit[sci]->Delete();

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

    gain_fiterr[sci]->Delete();

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

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

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

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

    gain_fit_isDouble[sci]->Delete();
    gain_fit_isEdge[sci]->Delete();
    gain_fit_isInnerEdge[sci]->Delete();
    gain_fit_isNextToDouble[sci]->Delete();
    gain_fit_isBulk[sci]->Delete();

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

    gain_ADUper1keV_2d[sci]->Delete();
    gainerr_ADUper1keV_2d[sci]->Delete();
  
  } // end of sc loop

}