using namespace std;

#include "f1f2in07_.h"
//#include "Fit_P_RTPC.h"
#include "RTPCECorr.h"
#include "bonana.h"
#include "readPBtable.h"
#include "readSKtable.h"
#include "SetBranchAddress.h"
#include "createHistos.h"
#include "createCanvas.h"
#include "good_e_cut.h"
#include "callOsi.h"
#include "clas_pid.h"
#include "setRunCond.h"
#include "cos_th_pq.h"
#include "check_all_z.h"
#include "delta0_spectator.h"
//#include "momcorr_rev1.h"
//#include "inipar_rv.h"
#include "momcorr_rv.cc"
#include "GetRunQCFlag.cc"
#include "slava_costh_pq.h"
#include "dedxbb.h"
#include "dqdx_hedme.h"
#include "initCounters.h"
#include "makeRunPlots.h"
#include "countsBGsub.h"
#include "doRadCorr.h"
#include "writeAccCorr.h"
#include "readAccCorr.h"
#include "ACCEPTANCE.h"
#include "F2dF2pratio.h"
#include "sigd_sign_elastic.h"
#include "RTPC_dqdxCut.h"
#include "F2nF2dratio_vs_w.h"
#include "F2nF2dratio_vs_x.h"
#include "F2dF2pratio_vs_x.h"
#include "F2dF2pratio_vs_w.h"
#include "F2pF2dratio_vs_x.h"
#include "F2pF2dratio_vs_w.h"
#include "F2nF2pratio_vs_w.h"
#include "F2nF2pratio_vs_x.h"
#include "F2pF2nratio_vs_w.h"
#include "F2pF2nratio_vs_x.h"
#include "du_vs_w.h"
#include "du_vs_x.h"
#include "accCut.h"
#include "plot_rats.h"
#include "plotAccCorr.h"
#include "plotJXAccCorr.h"
#include "Conta.h"
#include "bbc.h"
#include "doPimCorr.h"
#include "doEpemCorr.h"
#include "doAccCorr.h"

#include <fstream>
#include <iostream>
#include <iomanip>
#include <cmath>
#include <string>
#include <string.h>
#include <vector>
#include <stdlib.h>
#include <vector>
#include <time.h>


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

  /* if( argc != 4)
    {
      fprintf(stderr,"Usage: deeps.exe [flags] [beam_e 1,2,4,5] [limit 0,1] [limit number]\n");
      exit(0);
      }*/
TApplication theApp("App", &argc, argv);
gROOT->Reset();
TStopwatch timer; 						//Start CPU timer
  timer.Start();
time_t start, end;
double elapsed;
time_t rawtime;
struct tm * timeinfo;
  //style preferences
  gStyle->SetTitle(0);
  gStyle->SetOptStat(0);
  gStyle->SetPalette(1);
  gROOT->SetStyle("Plain");
  
  //##########################config block start#############################
  
  int r_w_accCorr = 0; //0 means read the acc corr file and apply the correction
                       //1 means create the file
  
  int pbsk = 0;        //0 means use peter bosted's xsections
                       //1 means use sebastian kuhn xsections
  
 int skimmed = 1;    //0 means no skim, 1 means a skim was applied

 Int_t choose = 5;   //12 = 5+4+2+1 pass
                     //11 = 5+4+2 pass
                     //5 = 5pass
                     //4 = 4 pass 
                      //2 = 2 pass
                      //1 = 1 pass                    
                      //0 = test run
                      //-1 = simulation 
  Int_t tester = 0;
  Int_t truncate = 1000000;
  Int_t iEndIndex = 8;
  int iStartRun;
  int iEndRun;

  //choose   = atof(*++argv);
  //tester   = atof(*++argv);
  //truncate = atof(*++argv);

  int ibgsub = 1; //dz background subtraction
  int irc    = 1; //radiative correction
  int ipimc  = 0; //pi minus contamination correction
  int iepemc = 0; //positron-electron background correction
  int iac    = 0; //jixie's acceptance correction

  if(choose == 12) //5+4+2+1 pass!
	{
	  iStartRun=49490;
	  //iEndRun=50236;

	  //iStartRun=50063;
	  iEndRun=50349;
	}

      if(choose == 11) //5+4+2 pass!
	{
	  iStartRun=49670;
	  //iEndRun=50236;

	  //iStartRun=50063;
	  iEndRun=50331;

	}
  if(choose == 5) //5pass
    {
      iStartRun=49808;
      //iEndRun=50236;

      //iStartRun=50063;
      iEndRun=50268;

    }
  if(choose == 4) //4pass
    {
      iStartRun=49490;//start dirty
      //iEndRun=49629;  //end dirty

      //iStartRun=49620;//slava test
      //iEndRun=49629;  //slava test

      //iStartRun=49670;//start clean
      //iEndRun=49796;//end clean
      
      //iStartRun=49799;//4pass He
      iEndRun=49807;

    }

  if(choose == 2) //2pass
    {
      iStartRun=50269;
      //iStartRun=50291;//start D      
      //iEndRun=50310;
      iEndRun=50331; //end 2 pass
    }

  if(choose == 1) //1pass
    {
      iStartRun=50332;//start H
      iEndRun=50349;
    }

  if(choose == 0) //test
    {
      iStartRun=49793;  
      iEndRun=49800;
      tester = 0;
    }

  if(choose == -1) //sim
    {
      iStartRun=1;
      iEndRun=1;
      iEndIndex = 52;
    }


  char InPath[]="./ntuples";
  char ProgressPath[]="./progress_report"; 
#ifdef KEEP_DQDX_HISTO
 char OutPath[]="./ntuples/RTPCdqdx";
#endif

#ifdef WRITE_VIP_FILE
 char VIPPath[]="./ntuples/VIPfiles";
#endif

 //RTPC cut values

 Double_t deltazcut    = 15.0;
 Double_t chi2cut      = 4.0;
 Double_t sdist_lowcut  = -3.0;// -3.0,-9.0;
 Double_t sdist_highcut = 8.0;//8.0,1.0;
 Double_t edist_lowcut  = -9.0;//-3.0,-9.0;
 Double_t edist_highcut = 8.0;//8.0,5.0;
 Int_t    npadcut      = 5;
 Double_t radcut       = 0;
 Double_t zhighcut     = 100.0;
 Double_t zlowcut      = -60.0; 
 
 int iCounter, iPro, iPi, iPassDz = 0, iFailDz = 0;
 //these are to make squaring faster
 Double_t aa, bb, ccc, dd, ee;
 Int_t bin, kk, pad, qq, ww, ie;

////cuts!///////////////////////////////////////////////////////////////////
 bool vert_cut, clas_mom_cut, timing_cut, co4, cc_cut, elas_w_cut;
 bool phi_cut, idfail, dqdx_cut;

 //misc stuff
 Double_t dedx_hi, dedx_lo, predict_dqdx;
 Double_t xx1, yy1, xx2, yy2;
 Double_t mm, b_lo, b_hi;
 mm = 2.67; b_lo = -183.33; b_hi = -83.33;

 ////For Slava's Momentum Correction routine.......... 
 const Int_t NNN = 10;
 Int_t    charge[NNN];
 Double_t z_vertold[NNN];
 Double_t phiDC[NNN];
 Double_t px_old[NNN];
 Double_t py_old[NNN];
 Double_t pz_old[NNN];
 Int_t    sec_num[NNN];
 Double_t mass[NNN];
 Int_t    iClaspid[NNN];
 Int_t    my_pid[NNN];
 Int_t    iPartNum[NNN];
 Double_t px_new[NNN];
 Double_t py_new[NNN];
 Double_t pz_new[NNN];
 Double_t z_vertnew[NNN]; 
 Double_t zaver;
 Int_t    n_particles;

 //counters for the number of CLAS particles of different ids

 Int_t n_els, n_pros, n_pim, n_pip;

 Double_t p_trk;//represents the true p in the tracking region
 Double_t yy;   //typical nu/E variable

 cout<<"Analyzing runs "<<iStartRun<<" to "<<iEndRun<<endl;
 createHistos();
 createCanvas();
 initCounters();
 readPBtable();
 readSKtable();
 inipar();

 //// Class declarations

 eg1bCorr CorrMom; //Slava's Momentum Correction
 extern int GetRunQCFlag(int *RunQCFlag);
 int RunQCFlag[1000];
 GetRunQCFlag(RunQCFlag);
 //show the array
 //for(int i=0;i<100;i++) cout<<"RunQCFlag["<<49490+i<<"]="<<RunQCFlag[i]<<endl;


 if(r_w_accCorr == 0) readAccCorr();

 //now I'll make an array of gaussian functions to fit to any q.e. or e peaks

 TF1* egau[100];

  for(ii=0; ii<100; ii++)
    {
     sprintf(title, "egau%d",ii);
     egau[ii] = new TF1(title,"gaus",0.85,1.0);//0.924,0.954);
    }

  //to check on mom corrections....
for(ii=0; ii<6; ii++)
{
  sprintf(title,"uncorrected w sec%d",ii+1);
  hwb4corr[ii]= new TH1D(title, "", Wbplot, lowWplot, highWplot);
  sprintf(title,"corrected w sec%d",ii+1);
  hwafcorr[ii]= new TH1D(title, "", Wbplot, lowWplot, highWplot);
  sprintf(title,"uncorrected dp/p vs phi sec%d",ii+1);
  h2_dpVphisec[ii] = new TH2F(title,title,40,-30,30 ,40, -0.1,0.1);
  sprintf(title,"corrected dp/p vs phi sec%d",ii+1);
  h2_dpVphisec[ii+6] = new TH2F(title,title,40,-30,30 ,40, -0.1,0.1);
  //h2_dpVphisec[ii+6]->Sumw2();

}

//to check on rtpc dqdx cut
 for(ii=0; ii<20; ii++)
   {
     sprintf(title,"rtpc momentum %d - %d",ii*10+50, ii*10+60);
     dqdx_rat_pbinned[ii]= new TH1D(title, title, 150, 0, 4);
     dqdx_rat_pbinned[ii]->SetLineWidth(2); 
     dqdx_rat_pbinned[ii]->SetLineColor(kBlue);
     dqdx_rat_pbinned[ii]->GetXaxis()->SetTitle("measured/expected dqdx");
     dqdx_rat_pbinned[ii]->GetXaxis()->SetTitleSize(0.075);
   }

 //cut counters
 Int_t itotRTPC = 0;
 Int_t ipassvz = 0;
 Int_t ipassnpt = 0;
 Int_t ipassx2 = 0;
 Int_t ipasssdist = 0;
 Int_t ipassedist = 0;
 Int_t ipassr = 0;

 Double_t p_calc; //determined electron p from theta
 Double_t sin2;
 Double_t fill_phi, dp_over_p;
 Double_t cthpq, alpha_s, E_s; //cos(theta_pq), lightcone p fraction, energy of the spec
 Double_t xgtsum = 0.0;
Double_t q2gtsum = 0.0;
Double_t xltsum = 0.0;
Double_t q2ltsum = 0.0;
 int ltcount= 0;int gtcount = 0;

// - - - - - - - - - - Getting tree from the disk file - - - - - - - - - - - - 
//Get old file, old tree and set top branch address
 Int_t iRun=0,iIndex=0,runExists=0,totalNumberEntries=0, process;
char InFileName[200];
char OutFileName[200];
char NewRootFile[200];
char SaveFile[200];
 FILE* outfile, *elecfile, *vipfile, *statusfile, *ratfile;
iCounter=0;

//make a dir to store the dq/dx information (intermediate stage - delete later)
char cmd[255];
#ifdef KEEP_DQDX_HISTO
 sprintf(cmd,"mkdir %s",OutPath);
 system(cmd);
#endif

#ifdef WRITE_VIP_FILE
 sprintf(cmd,"mkdir %s",VIPPath);
 system(cmd);
#endif


#ifdef RUN_BY_RUN
 sprintf(OutFileName,"./tag_untag_rat.txt");
 //outfile.Open(OutFileName,"w");
 ratfile = fopen(OutFileName,"w");
#endif

 sprintf(cmd,"rm %s/*.done",ProgressPath);
 system(cmd);
 //fprintf(vipfile, "p(mev) theta(deg) phi(deg) w(gev) w*(gev) delta(w-w*)(gev)\n");



//---------first count the total number in the target files

for(iRun=iStartRun;iRun<=iEndRun;iRun++)
{
  if(skimmed) iEndIndex = 0;
  for(iIndex=0;iIndex <= iEndIndex;iIndex++)
    {
      if(choose >= 0 && skimmed) sprintf(InFileName,"%s/clasnt_0%d_skimmed.root",InPath,iRun);
      else if(choose >= 0 && !skimmed) sprintf(InFileName,"%s/clasnt_0%d_0%d.root",InPath,iRun, iIndex);
      else
	{
	  if(iIndex < 10)
	    sprintf(InFileName,"%s/Bonus_G4NG3Sim_%d_0%d_el_NoBrem.root",InPath,iRun, iIndex);
	  else sprintf(InFileName,"%s/Bonus_G4NG3Sim_%d_%d_el_NoBrem.root",InPath,iRun, iIndex);
	}	
      //check this file exist or not
      Long_t id,size,flags,mt;
      Int_t iFound = gSystem->GetPathInfo(InFileName,&id,&size,&flags,&mt);
      if(iFound==1 || RunQCFlag[iRun - 49490] != 2 )
	{
	  //printf("file %s DNE\n",InFileName);
	  continue; //File not exist
	}
      setRunCond(iRun);
#ifdef MY_MASS
      if(target_mass != MY_MASS) continue;
#endif
      
      printf("Target Mass = %.4f GeV, torus = %dA, E_BEAM = %.3f GeV\n",
	     target_mass, torus, e_beam);
      
      TFile *myfile = new TFile(InFileName);
      TTree *mytree = (TTree*)myfile->Get("h10"); //get the tree address	
      if (myfile->IsZombie()) {continue;} //File is a zombie and we don't want our brains eaten
      //SetBranchAddress(mytree);
      //mytree->SetBranchStatus("*",1);
      
      UInt_t nentries = (UInt_t)mytree->GetEntries();
      totalNumberEntries+=nentries;
      
      double eps = 0.01;
      
      ie = -1;
      if(fabs(e_beam - twopass) < eps) ie = 0;
      else if(fabs(e_beam - fourpass) < eps) ie = 1;
      else if(fabs(e_beam - fivepass) < eps) ie = 2;
      else {printf("WARNING: ATTEMPTING TO RUN A NON-EXISTENT BEAM ENERGY, e_beam = %f!!!!!!\n\n\n", e_beam);}
      
      delete mytree;
      delete myfile;
      printf("...Run(%d)=%d entries (%d accumulated)\n",iRun,nentries,totalNumberEntries);
    }
  if(tester && totalNumberEntries > truncate) break; //truncate to this number of events for testing
 }
 
 printf("\n \n --->You just asked me to analyze %d entries dude!!\n\n",totalNumberEntries);
 
 //----------Start processing
 time (&start);
 for(iRun=iStartRun;iRun<=iEndRun;iRun++)
   {
     runExists = 0;
     goodecount = 0;
     vipcount = 0;
     crapcount = 0;
     if(skimmed) iEndIndex = 0;     
     for(iIndex=0;iIndex<=iEndIndex;iIndex++)
       {
	 if(tester && iCounter > truncate) break; //truncate to this number of events for testing
	 
	 if(choose >= 0 && skimmed) sprintf(InFileName,"%s/clasnt_0%d_skimmed.root",InPath,iRun);
	 else if(choose >= 0 && ! skimmed) sprintf(InFileName,"%s/clasnt_0%d_0%d.root",InPath,iRun, iIndex);
	 else
	   {
	     if(iIndex < 10)
	       sprintf(InFileName,"%s/Bonus_G4NG3Sim_%d_0%d_el_NoBrem.root",InPath,iRun, iIndex);
	     else sprintf(InFileName,"%s/Bonus_G4NG3Sim_%d_%d_el_NoBrem.root",InPath,iRun, iIndex);
	   }
	 
	 //check this file exist or not
	 Long_t id,size,flags,mt;
	 Int_t iFound = gSystem->GetPathInfo(InFileName,&id,&size,&flags,&mt);
	 if(iFound==1 || setRunCond(iRun)==0 || RunQCFlag[iRun - 49490] != 2 )
	   {
	     //printf("file %s DNE\n",InFileName);
	     continue; //File not exist
	   }
#ifdef MY_MASS
	 if(target_mass != MY_MASS) continue;
#endif	
	 //e_beam -= 0.020;
	 
	 runExists = 1;
	 cout << "Analyzing "<<InFileName<<" ..."<<endl;
	 
	 TFile *myfile = new TFile(InFileName);
	 TTree *mytree = (TTree*)myfile->Get("h10"); //get the tree address
	 
#ifdef WRITE_ELEC 
	 
	 sprintf(OutFileName,"%s/elec_%d_A0%d.clas",InPath,iRun, iIndex);	
	 //outfile.Open(OutFileName,"w");
	 elecfile = fopen(OutFileName,"w");
	 fprintf(elecfile, "eid|total p(MeV)|theta(rad)|phi(rad)[0,2PI]|vertex x, y, z (cm, clas CS!)\n");	
#endif

#ifdef WRITE_MISS_PRO 
	 
	 sprintf(OutFileName,"%s/misspro_%d_A0%d.clas",InPath,iRun, iIndex);
	 sprintf(cmd,"rm %s",OutFileName);
	 system(cmd);	
	 //outfile.Open(OutFileName,"w");
	 mpfile = fopen(OutFileName,"w");
	 //fprintf(mpfile, "eid|total p(MeV)|theta(rad)|phi(rad)[0,2PI]|vertex x, y, z (cm, clas CS!)\n");
	 fprintf(mpfile, "predicted p_perp(MeV)|p_para(MeV)|measured p_perp(MeV)|p_para(MeV)\n");
#endif
	 
#ifdef WRITE_VIP_FILE
	 
	 sprintf(OutFileName,"%s/vips_%d.txt",VIPPath,iRun);
	 //outfile.Open(OutFileName,"w");
	 vipfile = fopen(OutFileName,"w");
#endif
	 
	 SetBranchAddress(mytree);
	 mytree->SetBranchStatus("*",1);
	 
	 UInt_t nentries = (UInt_t)mytree->GetEntries();
	 cout<<"Number of entries ="<<nentries<<endl;
	 
	 //oldtree->Print();
	 for (UInt_t ll=0;ll<nentries; ll++) 
	   {
	     //printf("______Entry %d______\n",ll);
	     if(tester && iCounter > truncate) break; //truncate to this number of events for testing
	     mytree->GetEntry(ll);
	     iCounter++;
	     if(!(ll%100000) && ll)
	       { 
		 time (&end);
		 elapsed = difftime (end,start);
		 //end = clock();
		 //elapsed = ((double) (end - start)) / CLOCKS_PER_SEC;
		 time ( &rawtime );
		 timeinfo = localtime ( &rawtime );
		 printf("%7d/ %d events have been processed. (%.2f min to process %.2f percent of total)\n",
			ll,nentries, elapsed/60.0, (double)iCounter*100/totalNumberEntries);
		 double minperevent = elapsed/(float)iCounter/60.0;
		 printf("Avg = %.2f min / mil events, estimate %.2f min left\n\n",
			1000000*minperevent, minperevent*(totalNumberEntries - iCounter));
	       }
	     hi_gcpart->Fill(gcpart);
	     hi_gpart->Fill(gpart);   
	     //if(gpart < 2) continue;
	     theta_e = acos(cz[0]);
	     /*printf("theta_vert = %.5f rad, vz[0]/100.0 = %.5f meters\n",
	       theta_e,vz[0]*CM2M);
	       if(gpart > 10) printf("FAIL: gpart > 10 \n\n");*/
	     

	     process = 0;
	     if(skimmed) process = 1;
	     else if(good_e_cut()
		     && gpart <= 10
		     //&& theta_z_cut(double(theta_e),double(vz[0]*CM2M))
		     ) process = 1;
	     else process = 0;
				    
	     if(process
		&& callOsi() > 0
		//&& vz[0] > -68
		//&& vz[0] < -48
		)
	       {

    //printf("event %d passed fid and good e cut\n",ll);
    //n_particles  = gpart;
     iClaspid[0]       = 11;
     idfail = false;
     iPro = -1; 
     iPi = -1;
     goodecount++;

     //printf("PASS (%d):Good Job electron!\n\n",goodecount);

     kk=0;

     for (jj=0; jj<gpart; jj++)
       {
	 iClaspid[jj] = clas_pid(jj);
	 if(iClaspid[jj]!=0)
	   {
	     iPartNum[kk] = jj;
	     kk++;
	   }
	 //if the particle can't be identified, break out of this event
	 //if(iClaspid[jj] == 0 || cl_stat[jj] < 1 || dc[jj] == 0) 
	 //  {
	 //printf("ID FAILED For part %d!!\n", jj);
	 //    idfail = true;
	 //break;
	 //  }
       }
     //if(idfail) continue;
     
     /////do Slava's Momentum Correction on the electron.////
     
     n_els=0; n_pros=0; n_pim=0; n_pip=0;
     for(jj=0; jj<kk; jj++)
       {
	 z_vertold[jj] = vz[jj];
	 phiDC[jj]     = (double)atan2(tl1_cy[dc[iPartNum[jj]]-1],tl1_cx[dc[iPartNum[jj]]-1]);
	 if(phiDC[jj]<0) phiDC[jj] += 2*PI;
	 px_old[jj]    = (double)p[iPartNum[jj]]*cx[iPartNum[jj]];
	 py_old[jj]    = (double)p[iPartNum[jj]]*cy[iPartNum[jj]];
	 pz_old[jj]    = (double)p[iPartNum[jj]]*cz[iPartNum[jj]];
	 sec_num[jj]   = dc_sect[dc[iPartNum[jj]]-1];
	 if(iClaspid[iPartNum[jj]] == 11)
	   {
	     my_pid[jj]    = iClaspid[iPartNum[jj]];
	     charge[jj]    = -1;
	     mass[jj]      = (double)M_EL;
	     n_els++;
	   }
	 if(iClaspid[iPartNum[jj]] == 2212)
	   {
	     my_pid[jj]    = iClaspid[iPartNum[jj]];
	     charge[jj]    = +1;
	     mass[jj]      = (double)M_PRO;
	     iPro          = iPartNum[jj];
	     n_pros++;
	   }
	 if(iClaspid[iPartNum[jj]] == 221)
	   {
	     my_pid[jj]    = iClaspid[iPartNum[jj]];
	     charge[jj]    = +1;
	     mass[jj]      = (double)M_PION;
	     n_pip++;
	   }
	 if(iClaspid[iPartNum[jj]] == -221)
	   {
	     my_pid[jj]    = iClaspid[iPartNum[jj]];
	     charge[jj]    = -1;
	     mass[jj]      = (double)M_PION;
	     iPi           = iPartNum[jj];
	     n_pim++;
	   }
	 if(iClaspid[iPartNum[jj]] == 0)
	   {
	     my_pid[jj]    = -1000;
	     charge[jj]    = -1000;
	     mass[jj]      = -1000;
	   }
	 
	 /*if(iPartNum[jj] != jj)
	   {
	     my_print_flag = 1;
	     printf("CLAS detected %d particles\n",gpart);
	     printf("event %d claspart %d my_part %d has iClaspid %d, phiDC %.2f, sec %d\n",
		    ll, iPartNum[jj], jj, iClaspid[iPartNum[jj]], phiDC[jj]*RAD2DEG,sec_num[jj]);
	     printf("dc[%d]-1 = %d\n",iPartNum[jj],dc[iPartNum[jj]]-1);
	   }
	 */
       }
     //getchar());
     //if(gpart > 1 && !check_all_z()) continue;
     n_particles  = kk;
     e_prime = sqrtf(M_EL*M_EL + p[0]*p[0]);
     Q2      = 4 * e_beam * e_prime * sin(theta_e/2.0)*sin(theta_e/2.0);
     w       = M_N*M_N + 2*M_N*(e_beam-e_prime) - Q2;
     w       = sqrtf(w);
     sin2    = sin(theta_e/2)*sin(theta_e/2);
     p_calc  = e_beam/(1 + 2*(e_beam/target_mass)*sin2);
     dp_over_p = (p[0]-p_calc)/e_prime;

     phi_clas = atan2(cy[0],cx[0]); 
     if(phi_clas<0) phi_clas += 2*PI;//puts phi_clas into the bonus coordinate system
	    
     h2_zVphi_nocorr->Fill(phi_clas*RAD2DEG,10.0*vz[0] + 580.0,1.);
     //if(idfail) n_particles = 1;

     for(jj=0; jj < n_particles; jj++)
       {
         px_new[jj] = px_old[jj];
         py_new[jj] = py_old[jj];
         pz_new[jj] = pz_old[jj];
	 z_vertnew[jj] = z_vertold[jj];
       }
     
     CorrMom.CallCorr(px_new, py_new, pz_new, mass, n_particles, iRun,
		      charge, sec_num, phiDC, z_vertnew, my_pid, e_beam);//my_ebeam);

     /*     old way of calling momentum correction routine

	    Loop(iRun, e_beam, torus, target_mass, n_particles, charge, z_vertold, 
	  phiDC, px_old, py_old, pz_old, sec_num, mass, my_pid, 
	  px_new, py_new, pz_new, z_vertnew, zaver);*/

     //     printf("__________________Event %d_________________________________\n", ll);
     for(jj=0; jj<kk; jj++)
       {
	 p[iPartNum[jj]] = sqrtf(px_new[jj]*px_new[jj] + py_new[jj]*py_new[jj] + pz_new[jj]*pz_new[jj]);
	 
	 cx[iPartNum[jj]] = px_new[jj]/p[iPartNum[jj]];
	 cy[iPartNum[jj]] = py_new[jj]/p[iPartNum[jj]];
	 cz[iPartNum[jj]] = pz_new[jj]/p[iPartNum[jj]];
	 vz[iPartNum[jj]] = z_vertnew[jj];//zaver;
     
#ifdef CLAS_MOM_DEBUG
	     printf("old(tracknum=%d) px py pz vz = %.3f %.3f %.3f %.3f\n",iPartNum[jj],px_old[jj],py_old[jj],pz_old[jj],z_vertold[jj]);
	     printf("new(clasnum=%d) px py pz vz= %.3f %.3f %.3f %.3f\n",iPartNum[jj],px_new[jj],py_new[jj],pz_new[jj],z_vertnew[jj]);
#endif
	 
       }

#ifdef CLAS_MOM_DEBUG
	 for(jj=0; jj<gpart; jj++)
	   { 
	     //printf("(clasnum=%d, pid=%d) = %.3f %.3f %.3f %.3f %.3f\n", 
	     //    jj, iClaspid[jj], p[jj], p[jj]*cx[jj], p[jj]*cy[jj], p[jj]*cz[jj], vz[jj]);
	   }
	 printf("\n\n\n");
	 getchar();
#endif	 
  
	      hi_beta->Fill(b[0]);                     
	      e_prime = sqrtf(M_EL*M_EL + p[0]*p[0]);
	      phi_clas = atan2(cy[0],cx[0]);// + 11/RAD2DEG;    
	      if(phi_clas<0) phi_clas += 2*PI;//puts phi_clas into the bonus coordinate system
	      
              theta_e = acos(cz[0]);
	      Q2    = 4 * e_beam * e_prime * sin(theta_e/2.0)*sin(theta_e/2.0);
	      w     = M_N*M_N + 2*M_N*(e_beam-e_prime) - Q2;
	      w     = sqrtf(w);
	      nu    = (e_beam - e_prime);
	      x_bj  = Q2/(2*M_N*nu);
	      sin2    = sin(theta_e/2)*sin(theta_e/2);
	      p_calc  = e_beam/(1 + 2*(e_beam/target_mass)*sin2);
	      dp_over_p = (p[0]-p_calc)/e_prime;
	      yy     = nu/e_beam;
  
	      //make the y cut to get rid of pion and e+e- backgrounds
	      //printf("y = %.4f\n", yy);
	      if(yy < 0 || yy > 0.85/*e_beam*/) 
		{
		  //printf("y = %.2f, event %d dn pass y cut\n",yy, ll);
		  continue;
		}
	      //printf("y=%.2f, passed y cut\n",yy);
	       hi_phi_clas->Fill(phi_clas*RAD2DEG);
	       h2_zVphi_corr->Fill(phi_clas*RAD2DEG,10.0*vz[0] + 580.0,1.);
	       
#ifdef WRITE_ELEC
	      
	       fprintf(elecfile, "%d %.5f %.5f %.5f %.5f %.5f %.5f\n",eid, p[0]*1000.0, 
		       theta_e, phi_clas, vx[0], vy[0], vz[0]);
#endif
	       
	       


	      if(w>=lowW && w<=highW)
		{
	      h2_zVw->Fill(w, 10.0*vz[0] + 580.0, 1.);
	      hi_Q2->Fill(Q2);
	      //Wbin=int(100*w);	       //W-bin defined (.01 GeV bins)
	      //Determine W bin:
	      for (bin=0;bin<numWbins;bin++)
		{
		  if ((w > Wmin[bin]) && (w <= Wmax[bin]))
		    {
		      Wbin=bin;
		    }
		}
	      //Determine Q^2 bin:
	      for (bin=0;bin<numQ2bins;bin++)
		{
		  if ((Q2 > Q2_min[bin]) && (Q2 <= Q2_max[bin]))
		    {
		      Q2_bin=bin;
		    }
		}
	      //printf("Wbin,Q2_bin = %d %d\n",Wbin,Q2_bin);
	      if (Wbin<numWbins && Wbin>0 && Q2_bin>0 && Q2_bin<numQ2bins)
		{
	      
	      if(Wbin > hiWbin) hiWbin = Wbin;
	      if(Wbin < loWbin) loWbin = Wbin;
	      if(Q2_bin > hiQ2bin) hiQ2bin = Q2_bin;
	      if(Q2_bin < loQ2bin) loQ2bin = Q2_bin;
	      

	      inclusiveCount[Wbin][Q2_bin] += 1.0; 
	      if(r_w_accCorr == 0) 
		{
		  acIncCount[Wbin][Q2_bin] += 1.0/epsilon[ie][Wbin][Q2_bin]; 
		  sumEps2Inc[Wbin][Q2_bin] += 1.0/(epsilon[ie][Wbin][Q2_bin]*epsilon[ie][Wbin][Q2_bin]);
		}
		  
	      else 
		{
		  acIncCount[Wbin][Q2_bin] += 1.0; 
		  sumEps2Inc[Wbin][Q2_bin] += 1.0;
		}
	      
	      runningW[Wbin][Q2_bin] += w;
	      runningQ2[Wbin][Q2_bin] += Q2;
	      runningx[Wbin][Q2_bin] += x_bj;
	      runningelp[Wbin][Q2_bin] += p[0];
	      runningthe[Wbin][Q2_bin] += theta_e;
	      vz[0] = 10.0*vz[0] + 580.0; // Turn into mm, and to BONUS coordinate system
		  
	      h2_wVQ2->Fill(Q2, w, 1.);
	      hi_clas_z->Fill(vz[0]);
	      //if(w > lowWplot && w<highWplot) hi_w->Fill(w);

#ifdef CHECK_DELTA0
	      elas_w_cut = false;
	      if(w > 0.938 - 0.07 && w < 0.938 + 0.07) elas_w_cut = true;
              if(n_els == 1
                 && n_pros == 1
		 && elas_w_cut)
		{ hi_dphi->Fill((atan2(cy[iPro],cx[iPro]) - atan2(cy[0],cx[0]))*RAD2DEG);
		}
	      if(n_els == 1
		 && n_pros == 1
		 && n_pim == 1
		 )
		{
		  //printf("calling delta0 mom x check routine\n");
		  delta0_spectator(-1,0,iPro, iPi);
		}
#endif
	      ////////////////////////////////////BEGIN LOOKING FOR SPECTATOR PROTON////////////////
	      for(int j1=0; j1<gcpart; j1++)
		{
		  if(int(q_tot[j1]) == 0) break;
		  //first fill come raw histograms
		  hi_bo_z->Fill(z[j1]);
		  hi_npt->Fill(npd_track[j1]);
		  hi_x2->Fill(x2[j1]);
		  hi_edist->Fill(edist[j1]); 
		  hi_sdist->Fill(sdist[j1]);		  
		  hi_r->Fill(r_0[j1]);
		  
		  //increment some cut counters
		  
		   itotRTPC++;
		   if(z[j1] > zlowcut && z[j1] < zhighcut)      ipassvz++;
		   if(npd_track[j1] >= npadcut)                 ipassnpt++;
		   if(x2[j1] < chi2cut)                         ipassx2++;
		   if(sdist[j1] < sdist_highcut 
		      && sdist[j1] > sdist_lowcut)               ipasssdist++;
		   if(edist[j1] < edist_highcut 
		      && edist[j1] > edist_lowcut)               ipassedist++;
		   if(r_0[j1] > radcut)                         ipassr++;

		  
		  if(z[j1] > zlowcut
		     && z[j1] < zhighcut
		     && vz[0] > zlowcut
		     && vz[0]< zhighcut
		     //&& r_0[j1] > 0
		     && npd_track[j1] >= npadcut
		     && x2[j1]  < chi2cut
		     )
		    {
		      //dqdx_cut = true;
		            //find sdist on left and right
		      int left  = 0;
		      int right = 0;
		      if(phi[j1] > PI/2 && phi[j1] < 3*PI/2)
			{
			  right = 1;
			  left  = 0;
			}  
		      else 
			{
			  right = 0;
			  left  = 1;
			}
		      if(theta[j1]*RAD2DEG > 90) bo_z_back->Fill(z[j1]);
		      //printf("%d\n",ll);
		      if(edist[j1] < edist_highcut 
			 && edist[j1] > edist_lowcut
			 && sdist[j1] < sdist_highcut
			 && sdist[j1] > sdist_lowcut
			 ) 
			{
			  timing_cut = 1;
			}
		      else timing_cut = 0;
		      vert_cut = fabs(z[j1]-vz[0]) < deltazcut; // In mm
		      
		      h2_ed_sd_dz->Fill(sdist[j1],edist[j1],fabs(z[j1]-vz[0]));
		      //for elastic protons elec_phi - prot_phi = PI 
		      dphi = phi_clas - phi[j1];
		      if(dphi < 0 ) dphi += 2*PI;
		      if(dphi > 2*PI) dphi -= 2*PI;
		      phi_cut = (fabs(phi_clas - phi[j1])*RAD2DEG < 15) ;
		      //&& (phi[j1]*RAD2DEG > 100)  
		      
		      //p_corr = (p_bonus + 132.41*exp(-0.0376*p_bonus)/sin(theta_bonus));
		      p_bonus_mev = sqrtf(px[j1]*px[j1] + py[j1]*py[j1] + pz[j1]*pz[j1]); 
		      p_bonus_gev = M2G*p_bonus_mev;
		      //p_corr = Fit_P_by_R((double)r_0[j1],(double)z[j1],(double)theta[j1]*RAD2DEG, 450)/p_bonus_mev;
		      int ret=RTPCECorr::DoPCorr( r_0[j1],z[j1],theta[j1],phi[j1], 450, p_corr);
		      //ret = 1; //I just turned off the momentum correction.
		      if(ret < 0)
			{
			  p_corr = 1;
			}
		      else{p_corr /= p_bonus_mev;} //1/0.783845 (to scale the uncorrected momentum)
		      
		      
		      //p_corr = (1 + 37.2209 / p_bonus_mev - 0.179484) ;
		      
		      E_s = sqrtf(M_PRO*M_PRO + p_corr*p_bonus_gev*p_corr*p_bonus_gev);
		      aa = e_beam + M_DEU - p[0] - E_s;
		      bb = p[0]*cx[0] + M2G*p_corr*px[j1];
		      ccc = p[0]*cy[0] + M2G*p_corr*py[j1];
		      dd = p[0]*cz[0] + M2G*p_corr*pz[j1] - e_beam;
		      wstar = aa*aa - bb*bb - ccc*ccc - dd*dd;
		      wstar = sqrtf(wstar);
#ifdef CALIBRATION_RUN
		      //temp check for electrons in min-I bonus run
		      if(right)
			{
			  hi_sdistR->Fill(sdist[j1]);
			  hi_edistR->Fill(edist[j1]);
			}  
		      if(left)
			{
			  hi_sdistL->Fill(sdist[j1]);
			  hi_edistL->Fill(edist[j1]);
			}

		      if(fabs(theta_e - theta[j1])*RAD2DEG < 20
			 //&& w_cut
			 && fabs(phi[j1] - phi_clas)*RAD2DEG < 20
			 && x2[j1] != 0
			 //&& theta_bonus*RAD2DEG > 10
			 && r_0[j1] != 0
			 //&& vert_cut
			 //&& timing_cut
			 )
			{
			  if(right)
			    {
			      hi_dzR->Fill(z[j1]-vz[0]);
			      hi_dthetaR->Fill((theta_e - theta[j1])*RAD2DEG);	
			      hi_dphiR->Fill((phi_clas - phi[j1])*RAD2DEG);
			    }
			  if(left)
			    {
			      hi_dzL->Fill(z[j1]-vz[0]);
			      hi_dthetaL->Fill((theta_e - theta[j1])*RAD2DEG);	
			      hi_dphiL->Fill((phi_clas - phi[j1])*RAD2DEG);
			    }
			  h2_zbVzc->Fill(vz[0], z[j1], 1.);
			}
#endif


	      if(wstar>=lowW 
		 && wstar<=highW
		 && r_0[j1] > radcut
		 )
		{	
		  predict_dqdx = dqdx_hedme(p_bonus_mev,M_PRO*G2M,1);	
#ifdef DO_RTPC_PID_CUT
		  dqdx_cut = false;
		  if(RTPC_dqdxCut(iRun, dedx[j1]/predict_dqdx)) dqdx_cut = true;  
#else
		  dqdx_cut = true;
		  
#endif
		  if(vert_cut 
		     && timing_cut 
		     )
		    {
		      hi_rat_dqdx->Fill(dedx[j1]/predict_dqdx);
		      hi_diff_dqdx->Fill(dedx[j1] - predict_dqdx);
		      h2_bethe->Fill(p_bonus_mev,dedx[j1],1.);
		      hi_p_raw->Fill(p_bonus_mev);
		      h2_avgadc_st->Fill(p_bonus_mev,vtc[j1]*sin(theta[j1])/npts[j1],1.);
                      h2_avgadc->Fill(p_bonus_mev,vtc[j1]/npts[j1],1.);
		      h2_nhits->Fill(p_bonus_mev,npts[j1],1.);
		      
		      if(p_bonus_mev >= 49.5 && p_bonus_mev <=250.5)
			{
			  int mybin = int((p_bonus_mev - 50)/10);
			  //printf("p = %f, mybin = %d\n", p_bonus_mev, mybin);
			  dqdx_rat_pbinned[mybin]->Fill(dedx[j1]/predict_dqdx);
			}
		      
		      if(dqdx_cut)
			{
			  hi_w_pro->Fill(wstar);
			  hi_x_bj_pro->Fill(x_bj);
			  hi_p_pro->Fill(p_bonus_mev);
			  hi_p_corr->Fill(p_corr*p_bonus_mev);
			  h2_bethepro->Fill(p_bonus_mev,dedx[j1],1.);
			  vipcount++;
			}
		      else
			{
			  hi_w_other->Fill(wstar);
			  hi_x_bj_other->Fill(x_bj);
			  hi_p_others->Fill(p_bonus_mev);
			  crapcount++;
			}
		    }
		  //printf("W, Wstar, Q2bin,Wbin = %f %f %d %d\n",w,wstar,Q2_bin,Wbin);
		  //if(Q2 > 2.43) getchar();
		      //Wstarbin=int(100*wstar); //W-bin defined (.01 GeV bins)
		      //Determine Wstar bin:
		      for (bin=0;bin<numWbins;bin++)
			{
			  if ((wstar > Wmin[bin]) && (wstar <= Wmax[bin]))
			    {
			      Wstarbin=bin;
			    }
			}
		   
		      cthpq = cos_th_pq(theta_e, phi_clas, theta[j1],
					phi[j1], e_prime, Q2, p_bonus_gev);
		      /*printf("Nate's cos_th_pq = %f\n",acos(cthpq)*RAD2DEG);
		      cthpq = slava_costh_pq(px[j1]/1000,py[j1]/1000,pz[j1]/1000,e_beam);
		      printf("Slava's cos_th_pq = %f\n\n\n",acos(cthpq)*RAD2DEG);
		      */
		      alpha_s = (E_s - M2G*p_corr*pz[j1])/M_PRO;
		      
#ifdef CHECK_DELTA0
                      if(//gpart == 3
                         n_els == 1
                         && n_pros == 1
                         && n_pim == 1
			 && vert_cut
			 && timing_cut
			 && dqdx_cut
			 //&& acos(cthpq)*RAD2DEG > 100
			 //&& p_corr*p_bonus/M2G < 180
			 )
                        {
                          //printf("calling delta0 mom x check routine\n");
                          delta0_spectator(j1,p_corr*p_bonus_mev,iPro, iPi);
                        }
#endif
		      
		      //printf("wbin = %d wstarbin = %d Wmin,max = %f %f\n",
		      //	     Wbin, Wstarbin, Wmin[Wbin],Wmax[Wbin]);
		      if(timing_cut 
			 && dqdx_cut
			 )
			{
			  hi_cutdiff_z->Fill(z[j1]-vz[0]);
			  hi_diff_z->Fill(z[j1] - 10.0*z_vertold[0]- 580.0);
			  if(!vert_cut) iFailDz++;			  
			  if(vert_cut)  iPassDz++;
			  if(!vert_cut
			     && cthpq < ctpqmax[nctpqbins-1] 
			     && p_corr*p_bonus_mev < psmax[npsbins-1]
			     )
			    {
			      tagIncBGCount[Wstarbin][Q2_bin] += 1.0;   
			      if(r_w_accCorr == 0)
				{
				  actagIncBGCount[Wstarbin][Q2_bin] += 1.0/epsilon[ie][Wbin][Q2_bin];   
				  sumEps2tagIncBG[Wstarbin][Q2_bin] += 1.0/(epsilon[ie][Wbin][Q2_bin]*epsilon[ie][Wbin][Q2_bin]);
				}
			      else
				{ 
				  actagIncBGCount[Wstarbin][Q2_bin] += 1.0;   
				  sumEps2tagIncBG[Wstarbin][Q2_bin] += 1.0;
				}
			    }
			}
		      
		      if(vert_cut 
			 && dqdx_cut
			 )
			{
			  //hi_edist->Fill(edist[j1]); 
			  //hi_sdist->Fill(sdist[j1]);
			  h2_ed_sd->Fill(sdist[j1],edist[j1],1.);
			  h2_sd->Fill(z[j1],phi[j1],1.0);
			  h2_sdVz->Fill(z[j1], sdist[j1],1.);
			  //////for eric c.//////////////////////////////////////////////
			  if(timing_cut)
			    {
			      if(w*w > 3.5)
				{
				  h2_psVcostpqgt->Fill(cthpq, p_corr*p_bonus_mev,1.);
				  h2_psValphagt->Fill(alpha_s, p_corr*p_bonus_mev,1.);
				  gtcount++;
				  xgtsum += x_bj;
				  q2gtsum += Q2;
				}
			      else
				{
				  h2_psVcostpq->Fill(cthpq, p_corr*p_bonus_mev,1.);
				  h2_psValpha->Fill(alpha_s, p_corr*p_bonus_mev,1.);
				  ltcount++;
				  xltsum += x_bj;
				  q2ltsum += Q2;
				}
			    }
			  ////////////////////////////////////////////////////////////////
			  if(timing_cut
			     && cthpq < ctpqmax[nctpqbins-1] 
			     && p_corr*p_bonus_mev < psmax[npsbins-1]
			     )
			    {
#ifdef WRITE_VIP_FILE
			      fprintf(vipfile, "%d %.5f %.5f %.5f %.5f %.5f %.5f\n",eid, M2G*p_corr*px[j1], M2G*p_corr*py[j1],M2G*p_corr*pz[j1],z[j1], w, wstar);
#endif
			      
			      
			      
			      double jxAcc = ACCEPTANCE::GetAcceptanceInc(iRun, double(wstar), double(Q2), double(acos(cthpq)*RAD2DEG), double(p_corr*p_bonus_mev*M2G));
			      double jxEff = ACCEPTANCE::GetRTPCEff(iRun,double(wstar),double(Q2));
			      /*
				
				printf("acc(wstar,q2,thpq,ps) = %f(%.2f,%.2f,%.2f,%.2f)\n",jxAcc,wstar,Q2,acos(cthpq)*RAD2DEG,p_corr*p_bonus_mev*M2G);
				printf("_______________________\n");
			      
			      if(jxAcc > 0.0)
				{
				  actagIncCount[Wstarbin][Q2_bin] += 1.0/jxAcc;
				  }*/

			      

			      tagIncCount[Wstarbin][Q2_bin] += 1.0;
			      if(r_w_accCorr == 0)
				{ 
				  actagIncCount[Wstarbin][Q2_bin] += 1.0/epsilon[ie][Wbin][Q2_bin];
				  sumEps2tagInc[Wstarbin][Q2_bin] += 1.0/(epsilon[ie][Wbin][Q2_bin]*epsilon[ie][Wbin][Q2_bin]);
				}			      
			      else 
				{
				  actagIncCount[Wstarbin][Q2_bin] += 1.0;
				  sumEps2tagInc[Wstarbin][Q2_bin] += 1.0;
				}
			      
			      runningACC[Wstarbin][Q2_bin] += jxAcc;
			      runningRTPCeff[Wstarbin][Q2_bin] += jxEff;
			      runningps[Wstarbin][Q2_bin] += p_corr*p_bonus_mev;
			      runningctpq[Wstarbin][Q2_bin] += cthpq;
			      runningdelta[Wstarbin][Q2_bin] += (w - wstar);
			      runningWstar[Wstarbin][Q2_bin] += wstar;
			      h2_wVwstar->Fill(wstar, w, 1.);
			      //if(wstar > lowWplot && wstar<highWplot)hi_wstar->Fill(wstar);
			      hi_bo_phi->Fill(phi[j1]*RAD2DEG); 
			      
			      hi_bo_theta->Fill(theta[j1]*RAD2DEG); 
			      
			      h2_phi->Fill(phi_clas,phi[j1],1.);
			      hi_phi_clas->Fill(phi_clas*RAD2DEG);
			      hi_theta_clas->Fill(theta_e*RAD2DEG);
			      
			      h2_dthVdphi->Fill((phi_clas - phi[j1])*RAD2DEG,(theta_e - theta[j1])*RAD2DEG,1.);
			      //printf("clas %.1f bo %.1f\n", phi_clas*RAD2DEG, phi[j1]*RAD2DEG);
			      
			      h2_phiVdth->Fill((theta_e - theta[j1])*RAD2DEG);
			      
			    }
			}
		      
		}//end Wstar range check
		      
		    }//end if(good bonus track)
		  
		}//end loop over bonus tracks
		}//end good bin check
		}//end of W range check
   }//end good first electron check

 			}//end loop over entries
			delete mytree;
			delete myfile;
#ifdef WRITE_ELEC
			fclose(elecfile);
#endif

#ifdef WRITE_VIP_FILE
                        fclose(vipfile);
#endif

#ifdef WRITE_MISS_PRO
			fclose(mpfile);
#endif
	    }//end subrun loop

#ifdef RUN_BY_RUN
	  if(runExists)
	    {
	      runnum[truns] = (double)iRun;
	      viprat[truns] = (double)vipcount/(double)goodecount;
	      eviprat[truns] = viprat[truns]*sqrtf(1/(double)vipcount + 1/(double)goodecount);
	      craprat[truns] = (double)crapcount/(double)goodecount;
	      ecraprat[truns] = craprat[truns]*sqrtf(1/(double)crapcount + 1/(double)goodecount);	      
	      targetMass[truns]=target_mass;
	      fprintf(ratfile, "%d %d %d %f %f\n",iRun,vipcount,goodecount,(double)vipcount/(double)goodecount,eviprat[truns]);
	      truns++;
	      if(truns >= NN)
		{
		  printf("you need to increase array size NN in deeps.C\n");
		  truns = NN-1;
		}
	    }
#endif
	
	      
#ifdef KEEP_DQDX_HISTO
	      sprintf(SaveFile,"%s/RTPC_%d_dqdx.root",OutPath,iRun);
	      TFile *f = new TFile(SaveFile,"new");
	      hi_rat_dqdx->Write();
	      hi_rat_dqdx->Reset();
	      delete f;
#endif
	      sprintf(OutFileName,"%s/run_%d.done",ProgressPath,iRun);
	      statusfile = fopen(OutFileName,"w");
	      fprintf(statusfile,"%d good electrons in this file (%.2f Percent finished)\n", 
		      goodecount, (double)iCounter*100/totalNumberEntries);
	      fclose(statusfile); 
	      
	}//end run loop



/////////////////start manipulation and drawing of histograms (should be in a new routine some day!)
/******CANVAS 0****************
 gStyle->SetPalette(1);
can[0]->cd(1); hi_pre_px->Draw(); hi_bo_px->Draw("same");
can[0]->cd(2); hi_pre_py->Draw(); hi_bo_py->Draw("same");
can[0]->cd(3); hi_pre_pz->Draw(); hi_bo_pz->Draw("same");

can[0]->cd(4); hi_dpx->Draw();
can[0]->cd(5); hi_dpy->Draw();
can[0]->cd(6); hi_dpz->Draw();
*/

/******CANVAS 1******************
//can[1]->cd(1); h2_zVphi_nocorr->Draw("colz");
//can[1]->cd(2); h2_zVphi_corr->Draw("colz");
can[1]->cd(1); hi_clas_z->Draw();
can[1]->cd(2); hi_bo_z->Draw();
can[1]->cd(3); bo_z_back->Draw();
can[1]->cd(4); hi_cutdiff_z->Draw();
hi_diff_z->Draw("same");
//can[1]->cd(3); hi_edist->Draw();
//can[1]->cd(4); hi_p_corr->Draw();

*/
/*******CANVAS 2******************
 */


/*******CANVAS 3******************
//can[3]->cd(3); hi_effVQ2->Draw();
*/

#ifdef RUN_BY_RUN
if(truns != 0)
{
  makeRunPlots(0,1);
}
fclose(ratfile);
#endif

//calculate avg values per bin...
for (int qq=0;qq<numQ2bins;qq++)
  {
    for (int ww=0;ww<numWbins;ww++)
      {
	if (inclusiveCount[ww][qq] > 0)
	  {
	    //first we set the relative statistical errors on each measurement
	    //stateUntag[ww][qq] = sqrtf(inclusiveCount[ww][qq])/inclusiveCount[ww][qq];
	    stateUntag[ww][qq] = sqrtf(sumEps2Inc[ww][qq])/acIncCount[ww][qq];

	    //now get the averages of the electron variables
	    Q2avg[ww][qq]=runningQ2[ww][qq]/double(inclusiveCount[ww][qq]);
	    Wavg[ww][qq]=runningW[ww][qq]/double(inclusiveCount[ww][qq]);
	    xavg[ww][qq]=runningx[ww][qq]/double(inclusiveCount[ww][qq]);
	    elpavg[ww][qq]=runningelp[ww][qq]/double(inclusiveCount[ww][qq]);
	    theavg[ww][qq]=runningthe[ww][qq]/double(inclusiveCount[ww][qq]);
	    //printf("qq, ww, q2avg, wavg = %d %d %f %f\n ", qq, ww,  Q2avg[ww][qq], Wavg[ww][qq]);
	  }
	if (tagIncCount[ww][qq] > 0)
	  {
	    //first we set the relative statistical errors on each measurement
	    //stateTag[ww][qq] = sqrtf(tagIncCount[ww][qq])/tagIncCount[ww][qq];
	    stateTag[ww][qq] = sqrtf(sumEps2tagInc[ww][qq])/actagIncCount[ww][qq];

	    //now get the averages of the tagging variables
	    psavg[ww][qq]=runningps[ww][qq]/double(tagIncCount[ww][qq]);
	    ctpqavg[ww][qq] = runningctpq[ww][qq]/double(tagIncCount[ww][qq]);
	    deltaavg[ww][qq] = runningdelta[ww][qq]/double(tagIncCount[ww][qq]);
	    Wstaravg[ww][qq] = runningWstar[ww][qq]/double(tagIncCount[ww][qq]);
	    accavg[ww][qq] = runningACC[ww][qq]/tagIncCount[ww][qq];
	    effavg[ww][qq] = runningRTPCeff[ww][qq]/tagIncCount[ww][qq];
	    
	    //if(psavg[ww][qq] < 100)
	    //printf("psavg[%d][%d] = %f MeV/c \n",ww,qq,psavg[ww][qq]);
	  }
      }
  }


//can[0]->cd(5);


can[2]->cd(1);h2_psVcostpq->Draw("colz");
can[2]->cd(2);h2_psValpha->Draw("colz");
can[2]->cd(3);h2_psVcostpqgt->Draw("colz");
can[2]->cd(4);h2_psValphagt->Draw("colz");
if(gtcount != 0 && ltcount != 0)
{
  printf("avg x, Q2 for w2 gt 3.5 = %.2f %.2f\n", xgtsum/double(gtcount), q2gtsum/double(gtcount));
  printf("avg x, Q2 for w2 lt 3.5 = %.2f %.2f\n", xltsum/double(ltcount), q2ltsum/double(ltcount));
}

/////determine what value to use to scale the background distributions

Double_t a1pars[2],a2pars[2],peakpars[3];
Double_t m1, m2, b1, b2, zm, zint;
Double_t a1, a2, a3, asum, Rbg;

can[3]->cd(2); 
gPad->SetLogz();
h2_ed_sd_dz->Divide(h2_ed_sd_dz,h2_ed_sd,1.0,1.0);
h2_ed_sd_dz->Draw("colz");
can[3]->cd(4); 
h2_ed_sd->Draw("colz");

can[3]->cd(1); hi_cutdiff_z->Draw();
TF1* p1 = new TF1("p1","pol1",deltazcut,100);
hi_cutdiff_z->Fit(p1, "R");
p1->GetParameters(&a2pars[0]);
p1->SetLineColor(kRed);
TF1* p2 = new TF1("p2","pol1",-100,-deltazcut);
hi_cutdiff_z->Fit(p2, "R");
p2->SetLineColor(kBlue);
p2->GetParameters(&a1pars[0]);
gStyle->SetOptFit(1111);
TF1* g3    = new TF1("g3","gaus",-deltazcut,deltazcut);
hi_cutdiff_z->Fit(g3,"R");
g3->GetParameters(&peakpars[0]);
p1->Draw("same");p2->Draw("same");

zm = peakpars[1];
m1 = a1pars[1];
m2 = a2pars[1];
b1 = a1pars[0];
b2 = a2pars[0];
zint = (b2 - b1) / (m1 - m2);
printf("-----------Start background dz subtraction calculation-----------\n");
printf("m1, m2, b1, b2, zm, zint = %f %f %f %f %f %f\n",m1,m2,b1,b2,zm,zint);

asum = 0.5*(m2*zint + b2)*(fabs(-b2/m2 - zint) + fabs(-b1/m1 - zint));
a2   = 0.5*(m2*(zm + deltazcut) + b2)*fabs(-b2/m2 - (zm + deltazcut));
a1   = 0.5*(m1*(zm - deltazcut) + b1)*fabs(-b1/m1 - (zm - deltazcut));
a3   = asum - a1 - a2;
Rbg  = a3/(a1+a2);
printf("a1, a2, a3, asum, Rbg = %f %f %f %f %f\n", a1,a2,a3,asum,Rbg);
printf("Number outside (%d) and inside (%d) dz cut.  outside/total = %f\n",
       iFailDz, iPassDz, (float)iFailDz/((float)(iPassDz + iFailDz)));
if(zm+deltazcut < zint || zm-deltazcut > zint) printf("WARNING bg calculation may be wrong!\n");
printf("-----------------------------------------------------------------\n");

if(ibgsub) countsBGsub(actagIncCount,stateTag,Rbg);
///////////////////////////////////////////////////////////////////////
double avgQ2, sumQ2;
double totCounts, totCorrCounts;
int plotbin;
double filler, oldcontent, olderror;
double newcontent, addstate, addsyse, newerror;

avgQ2 = 0.0;
sumQ2 = 0.0;
totCounts = 0.0;
totCorrCounts = 0.0;
 for(qq = 0; qq < numQ2bins; qq++)
   {
     for(ww=0;ww<numWbins;ww++)
       {
	 if(Wmin[ww] > lowWplot && Wmax[ww] < highWplot)
	   {
	     totCounts     += inclusiveCount[ww][qq];
	     totCorrCounts += tagIncCount[ww][qq];
	     sumQ2         += runningQ2[ww][qq];
	     
	     filler = inclusiveCount[ww][qq];
	     plotbin = hi_w->FindBin(Wavg[ww][qq]);
	     oldcontent = hi_w->GetBinContent(plotbin);
	     olderror   = hi_w->GetBinError(plotbin);
	     hi_w->AddBinContent(plotbin, filler);
	     newcontent = hi_w->GetBinContent(plotbin);
	     addstate = filler*stateUntag[ww][qq];
	     addsyse  = filler*syseUntag[ww][qq];
	     olderror = olderror*olderror;
	     addstate = addstate*addstate;
	     addsyse  = addsyse*addsyse;
	     newerror = sqrtf(olderror + addstate + addsyse);
	     hi_w->SetBinError(plotbin,newerror);

	     filler = tagIncCount[ww][qq];
	     plotbin = hi_wstar->FindBin(Wstaravg[ww][qq]);
	     oldcontent = hi_wstar->GetBinContent(plotbin);
	     olderror   = hi_wstar->GetBinError(plotbin);
	     hi_wstar->AddBinContent(plotbin, filler);
	     newcontent = hi_wstar->GetBinContent(plotbin);
	     addstate = filler*stateTag[ww][qq];
	     addsyse  = filler*syseTag[ww][qq];
	     olderror = olderror*olderror;
	     addstate = addstate*addstate;
	     addsyse  = addsyse*addsyse;
	     newerror = sqrtf(olderror + addstate + addsyse);
	     hi_wstar->SetBinError(plotbin,newerror);

	     {
	       //hi_wstarbkg->Fill((Wmax[ww]+Wmin[ww])/2);
	       hi_wstarbkg->Fill((Wmax[ww]+Wmin[ww])/2,tagIncBGCount[ww][qq]);
	     }
 
	   }
       }
   }

printf("______total inc = %.0f, total tagged = %.0f, overall ratio = %.4f_____\n",
       totCounts, totCorrCounts, totCorrCounts/totCounts);
printf("avg Q2 for tag over untagged plot:%.2f\n", sumQ2/totCounts);
 printf("///////////////CUT REPORT/////////////////////\n");
printf("Pass bonus z cut (%.0f,%.0f): %d (%.2f\%)\n",zlowcut, zhighcut, ipassvz, 100.0*double(ipassvz)/double(itotRTPC));
printf("Pass number pads/track >= %d: %d (%.2f\%)\n",npadcut, ipassnpt, 100.0*double(ipassnpt)/double(itotRTPC));
printf("Pass helix chi 2 < %.0f:      %d (%.2f\%)\n",chi2cut, ipassx2, 100.0*double(ipassx2)/double(itotRTPC));
printf("Pass sdist cut (%.0f,%.0f):   %d (%.2f\%)\n",sdist_lowcut, sdist_highcut, ipasssdist, 100.0*double(ipasssdist)/double(itotRTPC));
printf("Pass edist cut (%.0f,%.0f):   %d (%.2f\%)\n",edist_lowcut, edist_highcut, ipassedist, 100.0*double(ipassedist)/double(itotRTPC));
 printf("Pass positive R cut:          %d (%.2f\%)\n",ipassr, 100.0*double(ipassr)/double(itotRTPC));
 printf("//////////////////////////////////////////////\n");
 can[3]->cd(3); hi_w->DrawCopy();
//can[3]->cd(4); hi_wstar->DrawCopy();
//can[3]->cd(5); hi_wstarbkg->DrawCopy();
//hi_wstarcorr->Add(hi_wstar,hi_wstarbkg,1,-Rbg);
hi_wstar->Scale(totCounts/totCorrCounts);
//can[3]->cd(6); hi_wstar->DrawCopy();
  hi_wrat->Divide(hi_wstar,hi_w,1.,1.);

 hi_wrat->SetMinimum(0.00);
//hi_wrat->SetMaximum(0.04);
//can[3]->cd(6); hi_wrat->DrawCopy();

#ifdef PLOT_RATS
can[8]->cd(1); hi_wstar->DrawCopy(); hi_w->DrawCopy("same");
#endif

can[10]->cd(6); 
hi_wstar->DrawCopy(); 
hi_wstar->Fit(egau[0],"R");
hi_w->DrawCopy("same");
//egau[0]->DrawCopy("same");

hi_w->Reset();
hi_wstar->Reset();
hi_wstarbkg->Reset();
hi_wrat->Reset();

if(r_w_accCorr == 1) writeAccCorr(iStartRun, iEndRun, pbsk); //0 is pb, 1 is sk

/*
     for(ww=0;ww<numWbins;ww++)
       {
	 printf("%f %f ",tagIncCount[ww][28],inclusiveCount[ww][28]);
       }
     printf("\n ");
*/
if(irc) doRadCorr(actagIncCount,acIncCount, syseTag, syseUntag);
/*
     for(ww=0;ww<numWbins;ww++)
       {
	 printf("%f %f ",tagIncCount[ww][28],inclusiveCount[ww][28]);
       }
     printf("\n ");
*/
//pion cont
if(ipimc){
  doPimCorr(tagIncCount,inclusiveCount, syseTag,syseUntag);
can[3]->cd(5); 
prof_pimCorrVw = h2_pimCorrVw->ProfileX();
prof_pimCorrVw->Draw();
}
//ps bkgnd cont
if(iepemc){
  doEpemCorr(tagIncCount,inclusiveCount, syseTag,syseUntag);
can[3]->cd(6); 
prof_epemCorrVw = h2_epemCorrVw->ProfileX();
prof_epemCorrVw->Draw();
 }
can[3]->cd(5); h2_zVphi_nocorr->Draw("colz");
can[3]->cd(6); h2_zVphi_corr->Draw("colz");
//CANVAS 5******************
#ifdef CALIBRATION_RUN

can[5]->cd(1); hi_dzL->Draw();
TF1* g11    = new TF1("g11","gaus",-15,15);
hi_dzL->Fit(g11,"R","same"); 
hi_dzR->Draw("same");
TF1* g12    = new TF1("g12","gaus",-15,15);
hi_dzR->Fit(g12,"R","same");

can[5]->cd(2); h2_zbVzc->Draw("colz");
TLine * line;
line = new TLine(-60,-60,100,100);
line->SetLineColor(kBlue); line->Draw("SAME");

can[5]->cd(3); hi_dthetaL->Draw();
TF1* g13    = new TF1("g13","gaus",-10,10);
hi_dthetaL->Fit(g13,"R","same");
hi_dthetaR->Draw("same");
TF1* g14    = new TF1("g14","gaus",-10,10);
hi_dthetaR->Fit(g14,"R","same");

can[5]->cd(4); hi_dphiL->Draw();
TF1* g15    = new TF1("g15","gaus",-10,10);
   hi_dphiL->Fit(g15,"R","same");
hi_dphiR->Draw("same");
TF1* g16    = new TF1("g16","gaus",-10,10);
   hi_dphiR->Fit(g16,"R","same");

can[5]->cd(5); hi_edistL->Draw();hi_edistR->Draw("same");
can[5]->cd(6); hi_sdistL->Draw();hi_sdistR->Draw("same");

#endif

#ifdef PLOT_RATS

//this function takes in canvas#, pad#, start Q2 bin, end Q2 bin, acceptance correction flag.
//it makes plots on the given canvas# and the one after it.
/*plot_rats(4,1,0,23);
plot_rats(4,2,24,24);
plot_rats(4,3,25,25);
plot_rats(4,4,26,26);
plot_rats(4,5,27,27);
plot_rats(4,6,28,28);
*/
if(choose == 4)
  {
    plot_rats(tagIncCount,inclusiveCount,stateUntag,syseUntag,stateTag,syseTag,4,1,25,25,iac);
    plot_rats(tagIncCount,inclusiveCount,stateUntag,syseUntag,stateTag,syseTag,5,1,26,26,iac);
    plot_rats(tagIncCount,inclusiveCount,stateUntag,syseUntag,stateTag,syseTag,6,1,27,27,iac);
    plot_rats(tagIncCount,inclusiveCount,stateUntag,syseUntag,stateTag,syseTag,7,1,28,28,iac);
    plot_rats(tagIncCount,inclusiveCount,stateUntag,syseUntag,stateTag,syseTag,15,1,29,29,iac);
    plot_rats(tagIncCount,inclusiveCount,stateUntag,syseUntag,stateTag,syseTag,16,1,30,30,iac);
    plot_rats(tagIncCount,inclusiveCount,stateUntag,syseUntag,stateTag,syseTag,17,1,31,31,iac);
  
    plot_rats(tagIncCount,inclusiveCount,stateUntag,syseUntag,stateTag,syseTag,9,1,25,31,iac);
}
if(choose == 5)
  {

    plot_rats(actagIncCount,acIncCount,stateUntag,syseUntag,stateTag,syseTag,4,1,27,27,iac);
    plot_rats(actagIncCount,acIncCount,stateUntag,syseUntag,stateTag,syseTag,5,1,28,28,iac);
    plot_rats(actagIncCount,acIncCount,stateUntag,syseUntag,stateTag,syseTag,6,1,29,29,iac);
    plot_rats(actagIncCount,acIncCount,stateUntag,syseUntag,stateTag,syseTag,7,1,30,30,iac);
    
    plot_rats(actagIncCount,acIncCount,stateUntag,syseUntag,stateTag,syseTag,15,1,31,31,iac);
    plot_rats(actagIncCount,acIncCount,stateUntag,syseUntag,stateTag,syseTag,16,1,32,32,iac);
    plot_rats(actagIncCount,acIncCount,stateUntag,syseUntag,stateTag,syseTag,17,1,33,33,iac);
    plot_rats(actagIncCount,acIncCount,stateUntag,syseUntag,stateTag,syseTag,18,1,27,33,iac);
    
    plot_rats(actagIncCount,acIncCount,stateUntag,syseUntag,stateTag,syseTag,9,1,27,34,iac);
    
  }


    //plot_rats(12,1,27,40);
/*acccorr
if(iac) 
  {
    //doAccCorr(tagIncCount,inclusiveCount, syseTag, syseUntag);
    
    plot_rats(actagIncCount,inclusiveCount,stateUntag,syseUntag,stateTag,syseTag,6,1,27,27,1);
    plot_rats(actagIncCount,inclusiveCount,stateUntag,syseUntag,stateTag,syseTag,6,2,28,28,1);
    plot_rats(actagIncCount,inclusiveCount,stateUntag,syseUntag,stateTag,syseTag,6,3,29,29,1);
    plot_rats(actagIncCount,inclusiveCount,stateUntag,syseUntag,stateTag,syseTag,6,4,30,30,1);
    
    plot_rats(actagIncCount,inclusiveCount,stateUntag,syseUntag,stateTag,syseTag,7,1,31,31,1);
    plot_rats(actagIncCount,inclusiveCount,stateUntag,syseUntag,stateTag,syseTag,7,2,32,32,1);
    plot_rats(actagIncCount,inclusiveCount,stateUntag,syseUntag,stateTag,syseTag,7,3,33,33,1);
    plot_rats(actagIncCount,inclusiveCount,stateUntag,syseUntag,stateTag,syseTag,7,4,34,34,1);
  
    //plot_rats(12,1,27,40);
    plot_rats(actagIncCount,inclusiveCount,stateUntag,syseUntag,stateTag,syseTag,9,1,27,40,1);
  }
 else
   {
    plot_rats(tagIncCount,inclusiveCount,stateUntag,syseUntag,stateTag,syseTag,6,1,0,23,0);
    plot_rats(tagIncCount,inclusiveCount,stateUntag,syseUntag,stateTag,syseTag,6,2,24,24,0);
    plot_rats(tagIncCount,inclusiveCount,stateUntag,syseUntag,stateTag,syseTag,6,3,25,25,0);
    plot_rats(tagIncCount,inclusiveCount,stateUntag,syseUntag,stateTag,syseTag,6,4,26,26,0);
    
    plot_rats(tagIncCount,inclusiveCount,stateUntag,syseUntag,stateTag,syseTag,7,1,35,35,0);
    plot_rats(tagIncCount,inclusiveCount,stateUntag,syseUntag,stateTag,syseTag,7,2,36,36,0);
    plot_rats(tagIncCount,inclusiveCount,stateUntag,syseUntag,stateTag,syseTag,7,3,37,37,0);
    plot_rats(tagIncCount,inclusiveCount,stateUntag,syseUntag,stateTag,syseTag,7,4,38,40,0);
  
    //plot_rats(12,1,27,40);
    plot_rats(tagIncCount,inclusiveCount,stateUntag,syseUntag,stateTag,syseTag,9,1,27,40,0);

    }*/
#endif


 gStyle->SetPalette(1);

#ifdef CHECK_DELTA0
can[8]->cd(1); hi_p_rtpc->Draw(); hi_p_pre->Draw("same");
can[8]->cd(2);  hi_pt_rtpc->Draw(); hi_pt_m->Draw("same");
//can[8]->cd(3);h2_dpVp->Draw("colz");
can[8]->cd(3);hi_pre_pz->Draw(); hi_bo_pz->Draw("same");
can[8]->cd(4); hi_mmiss->Draw();
can[8]->cd(5); hi_wdelta0->Draw();
//can[8]->cd(4); hi_pmzVpmx->Draw("colz");
//can[8]->cd(5); hi_pmzVpmy->Draw("colz");

can[8]->cd(6);h2_mom->Draw(); 
 TLine * line;
 line = new TLine(60,60,200,200);
 line->SetLineColor(kBlue); line->Draw("SAME");


 gStyle->SetPalette(1);
//gStyle->SetOptStat(1111);
//gStyle->SetOptFit(1111);
//can[9]->cd(1); hi_mmiss->Draw();
//can[9]->cd(2); hi_wdelta0->Draw();
 can[9]->cd(1); hi_theta_meas->Draw();
 can[9]->cd(2); hi_phi_diff->Draw();
can[9]->cd(3);h2_dpVp->Draw();
prof_dpVp = h2_dpVp->ProfileX();
can[9]->cd(4);prof_dpVp->Draw();
can[9]->cd(5); hi_pt_diff->Draw();
TF1* g10    = new TF1("g10","gaus",-1,1);
hi_pt_diff->Fit(g10,"R");

can[9]->cd(6); hi_mom_diff->Draw();
TF1* g1    = new TF1("g1","gaus",-1,1);
hi_mom_diff->Fit(g1,"R");

can[3]->cd(4); hi_pe->Draw();
can[3]->cd(5); hi_ppro->Draw();
can[3]->cd(6); hi_ppim->Draw();

#else


#endif
//can[10]->cd(1);h2_wVwstar->Draw("colz");
//can[10]->cd(2);h2_deltaVps->Draw("colz");

can[10]->cd(1);
//gStyle->SetOptStat(0000);
//gStyle->SetOptFit(0000);
//h2_bethe->SetStats(0);
gPad->SetLogz();
h2_bethe->Draw("colz");
can[10]->cd(2); //gPad->SetLogz();
//h2_bethepro->Draw("colz");

const double   LO_FIT=20.0, HI_FIT=200;   // Fitting range

double par[5];
TF1* bbcc[5];

for(ii = 0; ii < 5; ii++)
  {
    sprintf(astring, "bbc%d", ii);
    bbcc[ii] = new TF1(astring,bbc,30,250,3);
  }
//TF1* myfunc = new TF1("myfunc", "[0]/x + [1]" , LO_FIT, HI_FIT, 5);
bbcc[0]->SetLineWidth(2); bbcc[0]->SetLineColor(kBlack);
gStyle->SetOptFit(1111);

// Initial values of parameters (proton):
bbcc[0]->FixParameter(0, M_PRO*G2M); //mass of particle
bbcc[0]->FixParameter(1, 1);     //charge of particle
//bbcc[0]->FixParameter(2, 2);     //Z of medium
//bbcc[0]->FixParameter(3, 4);     //A of medium

prof_bethe = h2_bethepro->ProfileX();
prof_bethe->Draw();
prof_bethe->Fit("bbc0","RV");
bbcc[0]->GetParameters(par);  
bbcc[0]->SetParameters(par); 
bbcc[0]->Draw("same");

can[10]->cd(1);
bbcc[0]->Draw("same");

//Deuteron
bbcc[1]->SetParameter(0, M_DEU*G2M); //mass of particle
bbcc[1]->SetParameter(1, 1);         //charge of particle
bbcc[1]->SetParameter(2, par[2]);    //offset
bbcc[1]->SetLineColor(kRed);
bbcc[1]->SetLineWidth(2);
bbcc[1]->Draw("same");

//Triton
bbcc[2]->SetParameter(0, M_TRI*G2M); //mass of particle
bbcc[2]->SetParameter(1, 1);         //charge of particle
bbcc[2]->SetParameter(2, par[2]);    //offset
bbcc[2]->SetLineColor(kBlue);
bbcc[2]->SetLineWidth(2);
bbcc[2]->Draw("same");

//He3
bbcc[3]->SetParameter(0, M_HE3*G2M); //mass of particle
bbcc[3]->SetParameter(1, 2);         //charge of particle
bbcc[3]->SetParameter(2, par[2]);    //offset
bbcc[3]->SetLineColor(6);
bbcc[3]->SetLineWidth(2);
bbcc[3]->Draw("same");

//He4
bbcc[4]->SetParameter(0, M_HE4*G2M); //mass of particle
bbcc[4]->SetParameter(1, 2);         //charge of particle
bbcc[4]->SetParameter(2, par[2]);    //offset
bbcc[4]->SetLineColor(kOrange);
bbcc[4]->SetLineWidth(2);
bbcc[4]->Draw("same");

/* yy1 = 0; yy2 = 600;
 xx1 = (yy1 - b_lo)/mm;
 xx2 = (yy2 - b_lo)/mm;
TLine* l1 = new TLine(xx1,yy1,xx2,yy2);
 l1->Draw("same");
 xx1 = (yy1 - b_hi)/mm;
 xx2 = (yy2 - b_hi)/mm;
TLine* l2 = new TLine(xx1,yy1,xx2,yy2);
 l2->Draw("same");
mm=-(3.0/2.0); b_hi = 300; 
 xx1 = (0 - b_hi)/mm;
 xx2 = (300 - b_hi)/mm;
TLine* l3 = new TLine(xx1,0,xx2,300);
l3->Draw("same");*/
//can[5]->cd(2); gPad->SetLogy(); hi_betheprof->Draw();
/*
gStyle->SetOptFit(1111);
TF1* g30    = new TF1("g30","gaus",0.85,1.05);
TF1* f1    = new TF1("f1","1.0",0.0,2.0);
hi_w_pro->Fit(g30,"R");
Double_t par[3];
Double_t aaa;
g30->GetParameters(&par[0]);
aaa = par[0];
//printf("par = %f\n",aaa);
TF1* g4    = new TF1("g4","gaus",0.85,1.05);
hi_w_other->Fit(g4,"R");
Double_t bbb;
g4->GetParameters(&par[0]);
bbb = par[0];
aaa /= bbb;
printf("par = %f\n",aaa);
hi_w_other->Multiply(f1,aaa);
hi_w_other->Draw();
can[10]->cd(4); hi_w_pro->Draw();

TF1* g101    = new TF1("g101","gaus",0.888,.988);
hi_w_pro->Fit(g101,"R");*/

can[10]->cd(3); hi_p_corr->Draw(); hi_p_pro->Draw("same"); hi_p_others->Draw("same"); hi_p_raw->Draw("same");
 can[10]->cd(4); gPad->SetLogz();
 h2_avgadc->Draw("colz");
can[10]->cd(5); gPad->SetLogz();
h2_avgadc_st->Draw("colz");
//h2_nhits->Draw("colz");

//can[10]->cd(5); hi_rat_dqdx->Draw(); //gPad->SetLogy(); hi_x_bj_pro->Draw();hi_x_bj_other->Draw("SAME");
//can[10]->cd(6); hi_diff_dqdx->Draw();

plotAccCorr(27, 11, 1, pbsk);//Q2 bin, can#, pad#, peter or sebastian's xs table
plotAccCorr(28, 11, 2, pbsk);
plotAccCorr(29, 11, 3, pbsk);
plotAccCorr(30, 11, 4, pbsk);

plotAccCorr(31, 12, 1, pbsk);
plotAccCorr(32, 12, 2, pbsk);
plotAccCorr(33, 12, 3, pbsk);
plotAccCorr(34, 12, 4, pbsk);

/*
plotJXAccCorr(24, 11, 1);//Q2 bin, can#, pad#, peter or sebastian's xs table
plotJXAccCorr(25, 11, 2);
plotJXAccCorr(26, 11, 3);
plotJXAccCorr(27, 11, 4);

plotJXAccCorr(28, 11, 5);
plotJXAccCorr(29, 11, 6);
plotJXAccCorr(30, 11, 7);
plotJXAccCorr(31, 11, 8);

plotJXAccCorr(32, 11, 9);
plotJXAccCorr(33, 11, 10);
plotJXAccCorr(34, 11, 11);
plotJXAccCorr(35, 11, 12);
*/
for(ii=0; ii<20; ii++)
  {
    can[13]->cd(ii+1);
    dqdx_rat_pbinned[ii]->Draw();
  }


can[14]->cd(1); hi_bo_z->Draw();
can[14]->cd(2); hi_npt->Draw();
can[14]->cd(3); hi_x2->Draw();
can[14]->cd(4); hi_edist->Draw();
can[14]->cd(5); hi_sdist->Draw();
can[14]->cd(6); hi_r->Draw();

/*double tag;

for(ww = 0; ww < numWbins; ww++)
  {
    printf("%d ",ww);
    for(int qq=0;qq<numQ2bins;qq++)
      {
	tag = ac_tagInc[ww][qq];
	printf("%.2e ", tag);
      }
    printf("\n");
  }
*/

sprintf(NewRootFile,"%d_%d.root",iStartRun,iEndRun);
//TFile *hfile = gROOT->FindObject(NewRootFile); if (hfile) hfile->Close();
//TFile *hfile = new TFile(NewRootFile,"RECREATE","ROOT file with histograms");
// TFile *hfile = 0;
TFile *f = new TFile(NewRootFile,"recreate");
gROOT->GetList()->Write();

for(ii=0;ii<numCans;ii++)
{
  if(ii == 13) gStyle->SetTitleFontSize(0.16);
  if(ii == 14) gStyle->SetTitleFontSize(0.1);
  gStyle->SetLabelSize(0.07,"x");
  gStyle->SetLabelSize(0.07,"y");
  gStyle->SetTitleSize(0.09,"x");
  gStyle->SetTitleSize(0.09,"y");
  gStyle->SetTitleOffset(1.3,"x");
  gStyle->SetTitleOffset(1.0,"y");
  if(printCan[ii])
    {
      can[ii]->Modified();  can[ii]->Update();
      sprintf(astring, "deeps_can%d.gif",ii);
      can[ii]->Print(astring);
      sprintf(astring, "deeps_can%d.eps",ii);
      can[ii]->Print(astring);
      //sprintf(astring, "deeps_can%d.C",ii);
      //can[ii]->SaveSource(astring);
      can[ii]->Write();
    }
}

delete f;
gDirectory->Clear();

// (TFile*)gROOT->FindObject(NewRootFile); if (hfile) hfile->Close();
// hfile = (TFile*)gROOT->FindObject(NewRootFile); if (hfile) hfile->Close();
// hfile = new TFile(NewRootFile,"RECREATE","Demo ROOT file with histograms");
/* TFile *f = new TFile("test.root", "recreate"); 
 f->cd();
 cout << "Writing file:" << endl;
 f->Write();
 gDirectory->Clear(); //<====== 
 f->Close();
*/


  timer.Stop();							//Stop CPU timer
printf("Real Time: %.2f min to process %d events\n",timer.RealTime()/60.0, iCounter);
printf("--->%.2f min per million events\n",timer.RealTime()*1000000/60.0/(double)iCounter );
//cout<<"CPU Time: "<<timer.CpuTime()<<endl;

//theApp.Run();

return 0;
} //  End of main(...)