#include<iostream>
#include<fstream>
using namespace std;   
void CanvasPartition(TCanvas *C,const Int_t Nx = 2,const Int_t Ny = 2,
                     Float_t lMargin = 0.15, Float_t rMargin = 0.05,
                     Float_t bMargin = 0.15, Float_t tMargin = 0.05);
void multipad()
{
   gStyle->SetOptStat(0);
   TCanvas *C = (TCanvas*) gROOT->FindObject("C");
   if (C) delete C;
   C = new TCanvas("C","canvas",1024,640);
   C->SetFillStyle(4000);

   ifstream f1;
   ifstream f2;
   ifstream f3;
   double Evt, q0, Pvis, Q2, q3;
   f1.open ("DUNECCQE.dat");
   f2.open ("DUNECCQERPA.dat");
   f3.open ("DUNECC2p2h.dat");

   // Number of PADS
   const Int_t Nx = 3;
   const Int_t Ny = 2;
   // Margins
   Float_t lMargin = 0.12;
   Float_t rMargin = 0.05;
   Float_t bMargin = 0.15;
   Float_t tMargin = 0.05;
   // Canvas setup
   CanvasPartition(C,Nx,Ny,lMargin,rMargin,bMargin,tMargin);
   // Dummy histogram.

  TH1F *h1 = new TH1F("h1", "0 < q_{3}/GeV <0.2", 25,0.0,.5);  // here p1 u hv already used as variable so can nt use it as histogram name 
  TH1F *h2 = new TH1F("h2", "0.2 < q_{3}/GeV <0.3", 25,0.0,.5);
  TH1F *h3 = new TH1F("h3", "0.3 < q_{3}/GeV <0.4", 25,0.0,.5);  // here p1 u hv already used as variable so can nt use it as histogram name 
  TH1F *h4 = new TH1F("h4", "0.4 < q_{3}/GeV <0.5", 25,0.0,.5);
  TH1F *h5 = new TH1F("h5", "0.5 < q_{3}/GeV <0.6", 25,0.0,.5);  // here p1 u hv already used as variable so can nt use it as histogram name 
  TH1F *h6 = new TH1F("h6", "0.6 < q_{3}/GeV <0.8", 25,0.0,.5);
  TH1F *h11 = new TH1F("h11", "0 < q_{3}/GeV <0.2", 25,0.0,.5);  // here p1 u hv already used as variable so can nt use it as histogram name 
  TH1F *h22 = new TH1F("h22", "0.2 < q_{3}/GeV <0.3", 25,0.0,.5);
  TH1F *h33 = new TH1F("h33", "0.3 < q_{3}/GeV <0.4", 25,0.0,.5);  // here p1 u hv already used as variable so can nt use it as histogram name 
  TH1F *h44 = new TH1F("h44", "0.4 < q_{3}/GeV <0.5", 25,0.0,.5);
  TH1F *h55 = new TH1F("h55", "0.5 < q_{3}/GeV <0.6", 25,0.0,.5);  // here p1 u hv already used as variable so can nt use it as histogram name 
  TH1F *h66 = new TH1F("h66", "0.6 < q_{3}/GeV <0.8", 25,0.0,.5);
  TH1F *h11a = new TH1F("h11a", "", 25,0.0,.5);  // here p1 u hv already used as variable so can nt use it as histogram name 
  TH1F *h22a = new TH1F("h22a", "", 25,0.0,.5);
  TH1F *h33a = new TH1F("h33a", "", 25,0.0,.5);  // here p1 u hv already used as variable so can nt use it as histogram name 
  TH1F *h44a = new TH1F("h44a", "", 25,0.0,.5);
  TH1F *h55a = new TH1F("h55a", "", 25,0.0,.5);  // here p1 u hv already used as variable so can nt use it as histogram name 
  TH1F *h66a = new TH1F("h66a", "", 25,0.0,.5);
 

  while(!f1.eof())
       {
        f1>>Evt>>q0>>Pvis>>Q2>>q3;
       if(q3>0.0 && q3<0.2) {h1->Fill(q0);}
       if(q3>0.2 && q3<0.3) {h2->Fill(q0);}
       if(q3>0.3 && q3<0.4) {h3->Fill(q0);}
       if(q3>0.4 && q3<0.5) {h4->Fill(q0);}
       if(q3>0.5 && q3<0.6) {h5->Fill(q0);}
       if(q3>0.6 && q3<0.8) {h6->Fill(q0);}
       }
       
      while(!f2.eof())
       {
        f2>>Evt>>q0>>Pvis>>Q2>>q3;
       if(q3>0.0 && q3<0.2) {h11->Fill(q0);}
       if(q3>0.2 && q3<0.3) {h22->Fill(q0);}
       if(q3>0.3 && q3<0.4) {h33->Fill(q0);}
       if(q3>0.4 && q3<0.5) {h44->Fill(q0);}
       if(q3>0.5 && q3<0.6) {h55->Fill(q0);}
       if(q3>0.6 && q3<0.8) {h66->Fill(q0);}
       } 
       
       while(!f3.eof())
       {
        f3>>Evt>>q0>>Pvis>>Q2>>q3;
       if(q3>0.0 && q3<0.2) {h11a->Fill(q0);}
       if(q3>0.2 && q3<0.3) {h22a->Fill(q0);}
       if(q3>0.3 && q3<0.4) {h33a->Fill(q0);}
       if(q3>0.4 && q3<0.5) {h44a->Fill(q0);}
       if(q3>0.5 && q3<0.6) {h55a->Fill(q0);}
       if(q3>0.6 && q3<0.8) {h66a->Fill(q0);}
       }  
  
     TPad *pad[Nx][Ny];
     for (Int_t i=0;i<Nx;i++) {
     for (Int_t j=0;j<Ny;j++) {
     C->cd(0);
         // Get the pads previously created.
         char pname[16];
         sprintf(pname,"pad_%i_%i",i,j);
         pad[i][j] = (TPad*) gROOT->FindObject(pname);
         pad[i][j]->Draw();
         pad[i][j]->SetFillStyle(4000);
         pad[i][j]->SetFrameFillStyle(4000);
         pad[i][j]->cd();
         // Size factors
         Float_t xFactor = pad[0][0]->GetAbsWNDC()/pad[i][j]->GetAbsWNDC();
         Float_t yFactor = pad[0][0]->GetAbsHNDC()/pad[i][j]->GetAbsHNDC();
    
      }
   }
         C->cd();

         gStyle->SetOptStat(0);
      h1->Scale(1/8.665e25); 
      h11->Scale(1/8.665e25); 
      h11a->Scale(1/8.665e25);
      h11->GetXaxis()->SetTitle("Reconstructed available energy(GeV)");
      h11->GetYaxis()->SetTitle("d^{2}#sigma/dE_{avail}dq_{3}(10^{-42}cm^{2}/GeV^{2})");
      h11->GetYaxis()->CenterTitle(); 
      h11->GetYaxis()->SetTitleOffset(1.0);
      h11->GetXaxis()->CenterTitle();
      h11->SetLineWidth(1); 
      h11->SetLineWidth(1); 
      h11->SetLineColor(3); 
      h11a->SetLineColor(2);
      h1->SetMarkerStyle(21);
      h1->GetXaxis()->SetTitleSize(0.04);
      h1->GetYaxis()->SetTitleSize(0.04);
      h1->SetMaximum(0.5e-21);
       
      TH1F *h1c = (TH1F*) h1->Clone();
      h11->Add(h1c);
      h11a->Add(h11);
      h11->GetYaxis()->SetRangeUser(0.0,0.5e-21);
    
      h2->Scale(1/8.665e25);   // Here we have devide the events by No of Target* Ave Flux
      h22->Scale(1/8.665e25); 
      h22a->Scale(1/8.665e25); 
    //  h2->GetXaxis()->SetTitle("Reconstructed available energy(GeV)");
    //  h2->GetYaxis()->SetTitle("Events");
   //   h2->GetYaxis()->CenterTitle(); 
   //   h2->GetYaxis()->SetTitleOffset(1.0);
   //   h2->GetXaxis()->CenterTitle();
      h2->SetLineWidth(1); 
      h22->SetLineWidth(1); 
      h22->SetLineColor(3); 
      h22a->SetLineColor(2);
    //  h2->GetXaxis()->SetTitleSize(0.04);
    //  h2->GetYaxis()->SetTitleSize(0.04);
      h2->SetMaximum(0.5e-21);         

      TH1F *h2c = (TH1F*) h2->Clone();
      h22->Add(h2c);
      h22a->Add(h22);
      h22->GetYaxis()->SetRangeUser(0.0,0.5e-21);
     
      h3->Scale(1/8.665e25); 
      h33->Scale(1/8.665e25);
      h33a->Scale(1/8.665e25); 
     // h3->GetYaxis()->CenterTitle(); 
     // h3->GetYaxis()->SetTitleOffset(1.0);
     // h3->GetXaxis()->CenterTitle();
      h3->SetLineWidth(1); 
      h33->SetLineWidth(1); 
      h33->SetLineColor(3); 
      h33a->SetLineColor(2);
    //  h3->GetXaxis()->SetTitleSize(0.04);
    //  h3->GetYaxis()->SetTitleSize(0.04);
      h3->SetMaximum(0.5e-21);
         
      TH1F *h3c = (TH1F*) h3->Clone();
      h33->Add(h3c);
      h33a->Add(h33);
      h33->GetYaxis()->SetRangeUser(0.0,0.5e-21);

     h4->Scale(1/8.665e25); 
     h44->Scale(1/8.665e25);
     h44a->Scale(1/8.665e25);
   //  h4->GetYaxis()->CenterTitle(); 
   //  h4->GetYaxis()->SetTitleOffset(1.0);
   //  h4->GetXaxis()->CenterTitle();
     h4->SetLineWidth(1); 
     h44->SetLineWidth(1); 
     h44->SetLineColor(3); 
     h44a->SetLineColor(2);
  //   h4->GetXaxis()->SetTitleSize(0.04);
  //   h4->GetYaxis()->SetTitleSize(0.04);
     h4->SetMaximum(0.5e-21);
     
     TH1F *h4c = (TH1F*) h4->Clone();
     h44->Add(h4c);
     h44a->Add(h44);
     h44->GetYaxis()->SetRangeUser(0.0,0.5e-21);
       
     h5->Scale(1/8.665e25); 
     h55->Scale(1/8.665e25); 
     h55a->Scale(1/8.665e25); 
  //   h5->GetYaxis()->CenterTitle(); 
  //   h5->GetYaxis()->SetTitleOffset(1.0);
  //   h5->GetXaxis()->CenterTitle();
     h5->SetLineWidth(1);
     h55->SetLineColor(3); 
     h55a->SetLineColor(2); 
   //  h5->GetXaxis()->SetTitleSize(0.04);
   //  h5->GetYaxis()->SetTitleSize(0.04);
     h5->SetMaximum(0.5e-21);
         
     TH1F *h5c = (TH1F*) h5->Clone();
     h55->Add(h5c);
     h55a->Add(h55);
     h55->GetYaxis()->SetRangeUser(0.0,0.5e-21);
      
    h6->Scale(1/8.665e25); 
    h66->Scale(1/8.665e25); 
    h66a->Scale(1/8.665e25);
 //   h6->GetYaxis()->CenterTitle(); 
 //   h6->GetYaxis()->SetTitleOffset(1.0);
 //   h6->GetXaxis()->CenterTitle();
    h6->SetLineWidth(1); 
    h66->SetLineColor(3); 
    h66a->SetLineColor(2);
  //  h6->GetXaxis()->SetTitleSize(0.04);
  //  h6->GetYaxis()->SetTitleSize(0.04);
    h6->SetMaximum(0.5e-21);

    TH1F *h6c = (TH1F*) h6->Clone();
    h66->Add(h6c);
    h66a->Add(h66);
    h66->GetYaxis()->SetRangeUser(0.0,0.5e-21);
   
      pad[0][0]->cd(); h44->Draw();  h4->Draw("same");  h44a->Draw("same");
      pad[1][0]->cd(); h55->Draw();  h5->Draw("same");  h55a->Draw("same");
      pad[2][0]->cd(); h66->Draw();  h6->Draw("same");  h66a->Draw("same");


      pad[0][1]->cd(); h11->Draw();  h1->Draw("same");  h11a->Draw("same");
      leg = new TLegend(0.9,0.7,0.48,0.9);
      leg->SetHeader("GENIE-2.12.6");
      leg->AddEntry(h1,"Default","l");
      leg->AddEntry(h11,"+RPA","l");
      leg->AddEntry(h11a,"+RPA+2p2h","l");
      leg->Draw();
      pad[1][1]->cd(); h22->Draw();  h2->Draw("same");  h22a->Draw("same");
      pad[2][1]->cd(); h33->Draw();  h3->Draw("same");  h33a->Draw("same");
     
   
  
}
void CanvasPartition(TCanvas *C,const Int_t Nx,const Int_t Ny,
                     Float_t lMargin, Float_t rMargin,
                     Float_t bMargin, Float_t tMargin)
{
   if (!C) return;
   // Setup Pad layout:
   Float_t vSpacing = 0.0;
   Float_t vStep  = (1.- bMargin - tMargin - (Ny-1) * vSpacing) / Ny;
   Float_t hSpacing = 0.0;
   Float_t hStep  = (1.- lMargin - rMargin - (Nx-1) * hSpacing) / Nx;
   Float_t vposd,vposu,vmard,vmaru,vfactor;
   Float_t hposl,hposr,hmarl,hmarr,hfactor;
   for (Int_t i=0;i<Nx;i++) {
      if (i==0) {
         hposl = 0.0;
         hposr = lMargin + hStep;
         hfactor = hposr-hposl;
         hmarl = lMargin / hfactor;
         hmarr = 0.0;
      } else if (i == Nx-1) {
         hposl = hposr + hSpacing;
         hposr = hposl + hStep + rMargin;
         hfactor = hposr-hposl;
         hmarl = 0.0;
         hmarr = rMargin / (hposr-hposl);
      } else {
         hposl = hposr + hSpacing;
         hposr = hposl + hStep;
         hfactor = hposr-hposl;
         hmarl = 0.0;
         hmarr = 0.0;
      }
      for (Int_t j=0;j<Ny;j++) {
         if (j==0) {
            vposd = 0.0;
            vposu = bMargin + vStep;
            vfactor = vposu-vposd;
            vmard = bMargin / vfactor;
            vmaru = 0.0;
         } else if (j == Ny-1) {
            vposd = vposu + vSpacing;
            vposu = vposd + vStep + tMargin;
            vfactor = vposu-vposd;
            vmard = 0.0;
            vmaru = tMargin / (vposu-vposd);
         } else {
            vposd = vposu + vSpacing;
            vposu = vposd + vStep;
            vfactor = vposu-vposd;
            vmard = 0.0;
            vmaru = 0.0;
         }
         C->cd(0);
         char name[16];
         sprintf(name,"pad_%i_%i",i,j);
         TPad *pad = (TPad*) gROOT->FindObject(name);
         if (pad) delete pad;
         pad = new TPad(name,"",hposl,vposd,hposr,vposu);
         pad->SetLeftMargin(hmarl);
         pad->SetRightMargin(hmarr);
         pad->SetBottomMargin(vmard);
         pad->SetTopMargin(vmaru);
         pad->SetFrameBorderMode(0);
         pad->SetBorderMode(0);
         pad->SetBorderSize(0);
         pad->Draw();
      }
   }
}