%%%%%%%%%%%%%%%%%%%%%%%%%%
DSSSD_Energy_Calibration.C
%%%%%%%%%%%%%%%%%%%%%%%%%%

#include "TStrip.h"
#include "TStrip.cxx"

#include "TFile.h"
#include "TROOT.h"
#include "TH1.h"
#include "TKey.h"
#include "TTree.h"
#include "TTree.h"
#include "TCanvas.h"
#include "TH2.h"
#include "TF1.h"
#include "TStyle.h"
#include "TSystem.h"
#include "TSpectrum.h"
#include "TKey.h"
#include "TCanvas.h"
#include "TCutG.h"
#include "TColor.h"
#include "TObject.h"
#include "TMath.h"
#include <iostream>
#include <sstream>

using namespace std;

TCanvas* Results_print;
TH1F* Strip_calibrated;


int DSSSD_Energy_Calibration(TString root_file="DSSSD_energy_calibration.root",string Strip_cuts_from_file="1 15 3 5 2 11 10 16")
    {
     TString file=Read_root_file(root_file);
    
     if (!gSystem->AccessPathName(file.Data(),kFileExists))
         {
          printf("**************************************\n");
          printf("Reading %s\n",root_file.Data());
          printf("**************************************\n");
          TFile *File = new TFile(file.Data(),"READ");
         
          Int_t Strip_cut_file_number[16]={0};
          extractIntegerFromString(Strip_cuts_from_file,Strip_cut_file_number);


          printf("*********************************************\n");
          printf("Creating %s\n","Energy_Calibrated_DSSSD.root");
          printf("*********************************************\n");
          TFile* Calibrated_DSSSD = new TFile("Energy_Calibrated_DSSSD.root","RECREATE");

          TStrip* strip = new TStrip();
          TTree* DSSSD_tree = Create_DSSSD_tree(strip);
     
          TH2F* Strip_vs_PAD_histo;
          TString Strip_histo_name;
            
          for (int strip_number=10; strip_number<11; strip_number++)
              {
               Strip_histo_name=Form("BidiP%d",strip_number+1); 
               Results_print = new TCanvas(Strip_histo_name,Strip_histo_name);
               Results_print->Divide(2,3);
              
               if(Strip_cut_file_number[strip_number]==100+strip_number+1)
                 {
                  Results_print->cd(2);
                  printf("Calibrating strip %d \n",strip_number+1);
                  Strip_vs_PAD_histo=(TH2F*)File->Get(Strip_histo_name);
                  Strip_vs_PAD_histo->Draw("COLZ");
                 
                  Results_print->cd(4);
                  Float_t x[3];

                  for (int i=0;i<3;i++)
                      {
                       printf("*********************************************************\n");
                       printf("Reading strip cut BidiP%d_%d_user.C defined by user.\n",strip_number+1,i+1);
                       printf("*********************************************************\n");
                       gROOT->ProcessLine(Form(".x BidiP%d_%d_user.C",strip_number+1,i+1));
                       TH1D* histo_proj_Y = Strip_vs_PAD_histo->ProjectionY("Projection_Y",0,500000,Form("[BidiP%d_%d_user]",strip_number+1,i+1)); 

                       TF1* gauss_fit = new TF1("Gauss_fit","gaus",150000,350000);
                       if (i==0){Results_print->cd(1);}
                       if (i==1){Results_print->cd(3);}
                       if (i==2){Results_print->cd(5);}          
           
                       histo_proj_Y->Fit(gauss_fit,"R+");
                       strip->Set_channel_alpha_E(gauss_fit->GetParameter(1),i+1);
                       x[i]=gauss_fit->GetParameter(1);
                       }
                  
                  Results_print->cd(6);       
                  TGraph* Linear_fit = new TGraph(3,x,strip->dE);;
                  Linear_fit->Draw("a*");
                  TF1* Linear_function = new TF1("Linear_fit","pol1");
                  Linear_fit->Fit("Linear_fit");
       
                  Fill_strip_parmeters(strip,strip_number,Linear_function);                  
                  }
               else
                  {
                   Results_print = new TCanvas(Strip_histo_name,Strip_histo_name);
                   Results_print->Divide(2,3);
                   printf("Calibrating strip %d \n",strip_number+1);
                   Strip_vs_PAD_histo=Transf_Strip_vs_PAD_histo(File,Strip_histo_name,strip);
                   Calibrate_strip(strip,Strip_vs_PAD_histo,Strip_histo_name,strip_number,File,Strip_cut_file_number);
                   }
               DSSSD_tree->Fill();
               //DSSSD_tree->Print();
               }
     
           DSSSD_tree->Write();
           Results_print->Write();
           Strip_calibrated->Write();  
           //PAD_results_print->Write();
           //PAD_calibrated->Write();
    
           File->Close();
  
           return ;
           }
     else
         {
          printf("********************************************************************\n");
          printf("File %s do not exist. Aborting analysis.\n",root_file.Data());
          printf("********************************************************************\n");          
          return ;
          } 
     }


TString Read_root_file(TString root_file)
     {
      TString dir = gSystem->UnixPathName(gInterpreter->GetCurrentMacroName());
      dir.ReplaceAll("DSSSD_Energy_Calibration.C","");
      dir.ReplaceAll("/./","/");
      TString file=dir+root_file;
      return file;
      }


TTree* Create_DSSSD_tree(TStrip* strip)
       {
        TTree* DSSSD_tree = new TTree("DSSSD_Energy_Calibrated_Tree","DSSSD_Energy_Calibrated");
        DSSSD_tree->Branch("Strip_number",&strip->strip_number,"Strip_number/I"); 
        DSSSD_tree->Branch("Slope",&strip->slope,"Slope/F");
        DSSSD_tree->Branch("Slope_error",&strip->slope_error,"Slope_error/F");
        DSSSD_tree->Branch("Intercept",&strip->intercept,"Intercept/F");
        DSSSD_tree->Branch("Intercept_error",&strip->intercept_error,"Intercept_error/F");
        DSSSD_tree->Branch("Chi2_red",&strip->chi2_red,"Chi2_red/F");
        
        //DSSSD_tree->Print();
        return DSSSD_tree;
        }

TH2F* Transf_Strip_vs_PAD_histo(TFile* File, TString Strip_histo_name,TStrip* strip)
      {
       TH2F* Strip_vs_PAD_histo=(TH2F*)File->Get(Strip_histo_name);
       TH2F* Strip_vs_PAD_histo_transformed = new TH2F(Strip_histo_name,Strip_histo_name,Strip_vs_PAD_histo->GetNbinsX(),Strip_vs_PAD_histo->GetXaxis()->GetXmin(),Strip_vs_PAD_histo->GetXaxis()->GetXmax(),Strip_vs_PAD_histo->GetNbinsY(),Strip_vs_PAD_histo->GetYaxis()->GetXmin(),Strip_vs_PAD_histo->GetYaxis()->GetXmax());
       for (int i=1;i<=Strip_vs_PAD_histo->GetNbinsY();i++)
           {
            for (int m=1;m<=Strip_vs_PAD_histo->GetNbinsX();m++)
                {
                 Strip_vs_PAD_histo_transformed->SetBinContent(m+i,m,Strip_vs_PAD_histo->GetBinContent(i,m));
                 }
            }
       Results_print->cd(2);
       Strip_vs_PAD_histo_transformed->Draw("COLZ");
       return Strip_vs_PAD_histo_transformed;
       }




void Calibrate_strip(TStrip* strip,TH2F* Strip_vs_PAD_histo,TString Strip_histo_name,Int_t strip_number,TFile* File,Int_t Strip_cut_file_number[])
     {
      TSpectrum *Strip_vs_PAD_proy_PAD = new TSpectrum();
      TH1D* histo_proj_X = Strip_vs_PAD_histo->ProjectionX();
      Results_print->cd(4);
      Results_print->Update();
      Int_t No_found_peaks = Strip_vs_PAD_proy_PAD->Search(histo_proj_X,2);
      Float_t* xpeaks = Strip_vs_PAD_proy_PAD->GetPositionX();
      Results_print->Update();

      Float_t* xpeaks=Sort(xpeaks,No_found_peaks);
      Float_t x[3]={0};
      Float_t cut_with=0.8;

      for (int i=0;i<3;i++)
          {
           TF1 *alpha_line = new TF1("Alpha_line_fit","gaus",0.95*xpeaks[i],1.05*xpeaks[i]);
           Results_print->cd(4);
           histo_proj_X->Fit(alpha_line,"R+");
           Results_print->Update();
           Float_t mean = alpha_line->GetParameter(1);
           Float_t x_P1 = mean-cut_with*2.355*(alpha_line->GetParameter(2));
           Float_t y_P1 = Strip_vs_PAD_histo->GetYaxis()->GetXmin();
           Float_t y_P2 = Strip_vs_PAD_histo->GetXaxis()->GetXmax();
           Float_t x_P3 = mean+cut_with*2.355*(alpha_line->GetParameter(2));
        
           TCutG* CUT=Create_CUTS(Strip_histo_name,x_P1,y_P1,y_P2,x_P3,Strip_vs_PAD_histo,i,strip,strip_number);
           TH1D* histo_proj_Y = Strip_vs_PAD_histo->ProjectionY("Projection_Y",0,500000,Form("[BidiP%d_%d]",strip_number+1,i+1)); 
           
           TF1* gauss_fit = new TF1("Gauss_fit","gaus",150000,350000);
           if (i==0){Results_print->cd(1); histo_proj_Y->Fit(gauss_fit,"R");}
           if (i==1){Results_print->cd(3); histo_proj_Y->Fit(gauss_fit,"R");}
           if (i==2){Results_print->cd(5); histo_proj_Y->Fit(gauss_fit,"R");}          
                      
           strip->Set_channel_alpha_E(gauss_fit->GetParameter(1),i+1);
           x[i]=gauss_fit->GetParameter(1);
           }    
                  
      TGraph* Linear_fit = new TGraph(No_found_peaks,x,strip->dE);
      Results_print->cd(6); 
      Linear_fit->Draw("a*");
      TF1* Linear_function = new TF1("Linear_fit","pol1");
      Linear_fit->Fit("Linear_fit");
       
      Fill_strip_parmeters(strip,strip_number,Linear_function); 
      Calibrate_strip(File,strip_number,strip);
      } 
     

Float_t* Sort(Float_t *xpeaks,Int_t No_found_peaks)
         {
          Float_t temp;
          for (int i=0;i < (No_found_peaks-1);i++)
              {
               for (int j=i+1;j < No_found_peaks;j++)
                   {
                    if (xpeaks[j] < xpeaks[i])
                       {
                        temp=xpeaks[j];
                        xpeaks[j]=xpeaks[i];
                        xpeaks[i]=temp; 
                        }
                    }
               }
          return xpeaks;
          }

   

TCutG* Create_CUTS(TString Strip_histo_name,Float_t x_P1,Float_t y_P1,Float_t y_P2,Float_t x_P3,TH2F* Strip_vs_PAD_histo,int i,TStrip* strip,Int_t strip_number)
     {
      ofstream myfile;
      myfile.open(Form("BidiP%d_%d.C",strip_number+1,i+1));
      
      TCutG* cutg = new TCutG(Form("BidiP%d_%d",strip_number+1,i+1),5);
      cutg->SetLineColor(i+1);
      cutg->SetPoint(0,x_P1,y_P1);
      cutg->SetPoint(1,x_P1,y_P2);
      cutg->SetPoint(2,x_P3,y_P2);
      cutg->SetPoint(3,x_P3,y_P1);
      cutg->SetPoint(4,x_P1,y_P1);
      cutg->Draw("SAME");
      Results_print->cd(2);
      cutg->Draw("SAME");
               
      myfile << Form(" TCutG* cutg = new TCutG(\"%s\",5);     \
           \n cutg->SetPoint(0,%f,%f);                        \
           \n cutg->SetPoint(1,%f,%f);                        \
           \n cutg->SetPoint(2,%f,%f);                        \
           \n cutg->SetPoint(3,%f,%f);                        \
           \n cutg->SetPoint(4,%f,%f);                        \
           \n cutg->Draw(\"same\");",Form("BidiP%d_%d",strip_number+1,i+1),x_P1,y_P1,x_P1,y_P2,x_P3,y_P2,x_P3,y_P1,x_P1,y_P1);
      return cutg; 
      }


void extractIntegerFromString(string Strip_cuts_from_file,Int_t Strip_cut_file_number[])
     {
      stringstream ss(Strip_cuts_from_file);
      string tempStr;
      int found;
      while (ss >> tempStr) 
            { 
             if (stringstream(tempStr) >> found)
                { 
                 Strip_cut_file_number[found-1]=found;
                 }
             }
      }

void Fill_strip_parmeters(TStrip* strip,Int_t strip_number,TF1* Linear_function)
     {
      strip->strip_number=strip_number+1;    
      strip->slope = Linear_function->GetParameter(0);
      strip->slope_error = Linear_function->GetParError(0);
      strip->intercept = Linear_function->GetParameter(1);
      strip->intercept_error = Linear_function->GetParError(1);
      strip->chi2_red = Linear_function->GetChisquare();
      }

void Calibrate_strip(TFile* File,Int_t strip_number,TStrip* strip)
     {
      TH1F* Strip_uncalibrated = (TH1F*)File->Get(Form("HistoP%d",strip_number+1)); 
      Strip_calibrated = new TH1F(Form("HistoP%d",strip_number+1),Form("HistoP%d",strip_number+1),Strip_uncalibrated->GetNbinsX(),Strip_uncalibrated->GetXaxis()->GetXmin(),strip->slope*Strip_uncalibrated->GetXaxis()->GetXmax()+strip->intercept);
      for (int i=0;i<Strip_uncalibrated->GetNbinsX();i++)
          {
           Strip_calibrated->SetBinContent(i,Strip_uncalibrated->GetBinContent(i));
           }
      new TCanvas(Form("HistoP%d",strip_number+1),Form("HistoP%d",strip_number+1));
      Strip_calibrated->Draw();
      }


%%%%%%%%%
TStrip.h
%%%%%%%%%

#ifndef ROOT_TStrip
#define ROOT_TStrip

#include "TObject.h"


class TStrip : public TObject 
     {
      public:
              Float_t channel_alpha_E1;
              Float_t channel_alpha_E2;
              Float_t channel_alpha_E3;
              Float_t* dE;
              Float_t slope;
              Float_t slope_error;
              Float_t intercept;
              Float_t intercept_error;
              Float_t chi2_red;
              Int_t strip_number;
                          
      public:
             TStrip(Float_t line1_dE=3.812,Float_t line2_dE=3.522,Float_t line3_dE=3.294)
                   {
                    dE = new Float_t[3];
                    dE[0]=line1_dE;
                    dE[1]=line2_dE;
                    dE[2]=line3_dE;
                    }
             Int_t Get_channel_alpha_E(Int_t alpha_line);
             const Float_t Get_dE(Int_t alpha_line) const  {return dE[alpha_line-1];} 
             const Float_t Get_slope(Int_t strip_number) const {return slope;}
             const Float_t Get_slope_error(Int_t strip_number) const {return slope_error;}
             const Float_t Get_intercept(Int_t strip_number) const {return intercept;}
             const Float_t Get_intercept_error(Int_t strip_number) const {return intercept_error;}
             const Float_t Get_chi2_red(Int_t strip_number) const {return chi2_red;}
             const Int_t Get_strip_number() const {return strip_number;}
             void Set_channel_alpha_E(Float_t alpha_channel, Int_t alpha_line);
             void Set_strip_number(Int_t number){strip_number = number;}
           
           
      ClassDef(TStrip,1)
      };
 
#endif 

%%%%%%%%%%%
TStrip.cxx
%%%%%%%%%%%

#include "TStrip.h"


ClassImp(TStrip);

Float_t TStrip::Get_channel_alpha_E(Int_t alpha_line)
              {
               switch (alpha_line)
                      {
                       case 1:
                              return channel_alpha_E1;
                              break;
                       case 2:
                              return channel_alpha_E2;
                              break;
                       case 3:
                              return channel_alpha_E3;
                              break;
                       }
               }
 
void TStrip::Set_channel_alpha_E(Float_t alpha_channel, Int_t alpha_line)
           {
            switch (alpha_line)
                   {
                    case 1:
                           channel_alpha_E1 = alpha_channel;
                           break;
                    case 2:
                           channel_alpha_E2 = alpha_channel;
                           break;
                    case 3:
                           channel_alpha_E3 = alpha_channel;
                           break;
                    }
            }