#include "TFile.h"
#include "TTree.h"
#include "TH1.h"
#include "TH2F.h"
#include "TF1.h"
#include "TROOT.h"
#include "TString.h"
#include <stdio.h>
#include <stdlib.h>
#include <bitset>

int plot3(){
   const int n = 17051;
    Double_t x[n]=0,y[n]=0,ex[n]=0,ey[n]=0;
   TCanvas *c1 = new TCanvas("c1","multipads",900,700);
   c1->SetGrid();
   
   //gStyle->SetOptStat(0);c1->Divide(2,2);
   Double_t DoseDeposited;
   Double_t Polarization;
   Int_t ExperimentNum;
   
    TFile *f = new TFile("PDPout.root");

   f->Get("PDPvar");
   // stuff
   PDPvar->SetBranchAddress("DoseDeposited",&DoseDeposited);
   PDPvar->SetBranchAddress("Polarization",&Polarization);  
   PDPvar->SetBranchAddress("ExperimentNum",&ExperimentNum);
   
   /*
   ExperimentNum coding
   # - Material - Temperature - Beam Voltage - Magnetic Field - Experiment Name
   1 - LiD - 1K - 48 GeV - N/A - Nucl. instr. and meth res. manual p35
   2 - NH3 - N/A - N/A - N/A - Annealing, Proceeding of the 4th... p28
   3 - NH3 - 1K - 2.5 GeV - 2.5 T - Proceeding of the 4th... p29 graph 3
   4 - NH3 - 1K - 2.5 GeV - 2.5 T - Proceeding of the 4th... p30 graph 1
   5 - pre-irradiated NH3 - 1K - 2.5 GeV - 2.5 T - Proceeding of the 4th... p30 graph 1   
   6 - pre-irradiated ND3 - 1K - 2.5 GeV - 2.5 T - Proceeding of the 4th... p30 graph 2
   7 - ND3 - N/A - N/A - 2.5 T - Annealing - Proceeding of the 4th... p31 graph 2
   http://p25ext.lanl.gov/~ming/SantaFe-DY/talks/Monday/Maxwell_DY.pdf
   8 - NH3 - 1K - 6.0 GeV - 5 T - SANE 72417 -  negative pol.
   9 - NH3 - 1K - 6.0 GeV - 5 T - SANE 72417 - positive pol.
   10- NH3 - 1K - 6.0 GeV - 5 T - SANE 72428 - negative pol.
   11- NH3 - 1K - 6.0 GeV - 5 T - SANE 72428 - positive pol.
   12- NH3 - 1K - 6.0 GeV - 5 T - SANE 72986 - negative pol.
   13- NH3 - 1K - 6.0 GeV - 5 T - SANE 72986 - positive pol.
   14- run 72162 positive pol.
   15- run 72162 negative pol.
   16- run 72164 negative pol.
   17- run 72164 positive pol.
   18- run 72657 negative pol.
   19- run 72657 positive pol.
   20- run 72669 negative pol.
   21- run 72669 posiitve pol.
   22- run 72824 negative pol.
   23- run 72824 positive pol.
   24- run 72837 negative pol.
   25- run 72837 positive pol.
   26- run 72984 negative pol.
   27- run 72984 positive pol.
   28- run 72986 negative pol.
   29- run 72986 positive pol.
   30- jlab2001 proton polarization
   31- eg1 dvcs dose top negative pol.
   32- eg1 dvcs dose top positive pol.
   33- eg1 dvcs dose bottom negative pol.
   34- eg1 dvcs dose bottom positive pol.
    */
   Int_t nentries = (Int_t)PDPvar->GetEntries();   
   //cout<<nentries<<endl;

  
   for (Int_t j=0;j<nentries;j++) { PDPvar->GetEntry(j);
       if (ExperimentNum == 31)
       {
           //start of te_10
           //TString = gSystem->GetFromPipe("date -d@'"+DateTime+"' +%s");
           //if(time>0 && time<30){
           x[j]  = DoseDeposited;
           
           y[j]  = Polarization;
      
	   if (x[j] > .5)
           {
	   if (y[j] < y[j-1]) 
	   {	
		y[j] = y[j-1];
	   }
	   else if (y[j] - y[j-1] > 1.3)
	   {
  
		y[j] = y[j-1];
	    }
	   
           }
	   
	   cout<<y[j]<<endl;
           //y[j]  = He4Press;
           //ex[j] = 0;
           //ey[j] = -1*0.017*Polarization; //uncertainty in TE fit
           //ey[j] = 0.01*He4Press;
       }
       //else
       //{
       //}
   }

   TGraphErrors *gra = new TGraphErrors(j,x,y,ex,ey);
   gra->SetTitle("TGraphErrors Example");
   gra->SetMarkerColor(4);
   gra->SetMarkerStyle(20);
   gra->SetMarkerSize(0.8);
   gra->Draw("AP");     
   c1->Update();
   
   //TF1 *a = new TF1("f", "log(x*[0])");
   //gra->Fit(a);
   //gra->Draw("AL");

   //Double_t par[9];
   //TF1 *g1 = new TF1("g1","guas",0,x[j]/3);
   //TF1 *g2 = new TF1("g1","guas",x[j]/3, x[j]*2/3);
   //TF1 *g1 = new TF1("g1","guas",2*x[j]/3,x[j]);
   //TF1 *total = new TF1("total","guas(0)+guas(3)+guas(6)",0,x[j]);
   //total->SetLineColor(2);
   //h->Fit(g1,"R");
   //h->Fit(g2,"R+");
   //h->Fit(g3,"R+");
   //g1->GetParameters(&par[0]);
   //g2->GetParameters(&par[3]);
   //g3->GetParameters(&par[6]);
   //total->SetParameters(par);
   //h->Fit(total,"R+"); 

    //RooPlot* frame = x.frame() ;
    //data->plotOn(frame,Cut("fabs(y)<2")) ; 
	
}