/////////////this block was in plot_rats().  the purpose was to perform an acceptance correction using 
//peter bosted's model.  i needed to also control the error bars, this is probably the only thing i will
//get out of this code block
	     
	     //determine the delta at this ps and cthpq bin 
	     /*psbin = ctpqbin = -1;
	       for (int bin=0;bin<npsbins;bin++)
	       {
	       if ((psavg[ww][qq] > psmin[bin]) && (psavg[ww][qq] <= psmax[bin]))
	       {
	       psbin=bin;
	       }
	       }
	       for (int bin=0;bin<nctpqbins;bin++)
	       {
	       if ((ctpqavg[ww][qq] > ctpqmin[bin]) && (ctpqavg[ww][qq] <= ctpqmax[bin]))
	       {
	       ctpqbin=bin;
	       }
	       }
	     
	     //if(psbin > -1 && ctpqbin > -1)
	     //{
	     Wbin = ww;
	     wstar_pl_delta = Wstaravg[ww][qq] + deltaavg[ww][qq];
	     for (int bin=0;bin<numWbins;bin++)
	       {
		 if ((wstar_pl_delta > Wmin[bin]) && (wstar_pl_delta <= Wmax[bin]))
		   {
		     Wbin=bin;
		   }
	       }
	     if(rc_untagInc[ww][qq] != 0)
	       {
		 plotbin   = hi_x->FindBin(xfill);
		 oldcontent = hi_x->GetBinContent(plotbin);
		 olderror   = hi_x->GetBinError(plotbin);
		 //hi_x->AddBinContent(plotbin,d2mrat[ww][qq]*rc_untagInc[ww][qq] );
		 hi_x->AddBinContent(plotbin,rc_untagInc[ww][qq] );
		 newcontent = hi_x->GetBinContent(plotbin);
		 if(oldcontent != 0)
		   {
		     newerror = newcontent*sqrtf(olderror*olderror/(oldcontent*oldcontent) 
						 + 1/rc_untagInc[ww][qq]); 
		     //+ d2mratPerErr[ww][qq]*d2mratPerErr[ww][qq]);
		   }
		 else newerror = newcontent*sqrtf(1/rc_untagInc[ww][qq]);// + d2mratPerErr[ww][qq]*d2mratPerErr[ww][qq]);
		 hi_x->SetBinError(plotbin,newerror);
		 
		 //printf("x %d, old %f+-%f, new %f+-%f\n", plotbin, oldcontent, olderror, newcontent, newerror);
		 
		 plotbin   = hi_w->FindBin(Wavg[ww][qq]);
		 oldcontent = hi_w->GetBinContent(plotbin);
		 olderror   = hi_w->GetBinError(plotbin);
		 //hi_w->AddBinContent(plotbin,d2mrat[ww][qq]*rc_untagInc[ww][qq] );
		 hi_w->AddBinContent(plotbin,rc_untagInc[ww][qq] );
		 newcontent = hi_w->GetBinContent(plotbin);
		 if(oldcontent != 0)
		   {
		     newerror = newcontent*sqrtf(olderror*olderror/(oldcontent*oldcontent) 
						 + 1/rc_untagInc[ww][qq]); 
		     //+ d2mratPerErr[ww][qq]*d2mratPerErr[ww][qq]);
		   }
		 else newerror = newcontent*sqrtf(1/rc_untagInc[ww][qq]);//+ d2mratPerErr[ww][qq]*d2mratPerErr[ww][qq]);
		 hi_w->SetBinError(plotbin,newerror);
		 
		 //printf("w %d, old %f+-%f, new %f+-%f\n", plotbin, oldcontent, olderror, newcontent, newerror);
	       }
	     
	     if(rc_bgc_tagInc[ww][qq] != 0)
	       {
		 plotbin   = hi_xstar->FindBin(xfill);
		 oldcontent = hi_xstar->GetBinContent(plotbin);
		 olderror   = hi_xstar->GetBinError(plotbin);
		 //hi_xstar->AddBinContent(plotbin,d2mrat[Wbin][qq]*rc_bgc_tagInc[ww][qq] );
		 hi_xstar->AddBinContent(plotbin,rc_bgc_tagInc[ww][qq] );
		 newcontent = hi_xstar->GetBinContent(plotbin);
		 if(oldcontent != 0)
		   {
		     newerror = newcontent*sqrtf(olderror*olderror/(oldcontent*oldcontent) 
						 + 1/rc_bgc_tagInc[ww][qq]); 
		     //+ d2mratPerErr[Wbin][qq]*d2mratPerErr[Wbin][qq]);
		   }
		 else newerror = sqrtf(1/rc_bgc_tagInc[ww][qq]);// + d2mratPerErr[Wbin][qq]*d2mratPerErr[Wbin][qq]);
		 hi_xstar->SetBinError(plotbin,newerror);
		 
		 //printf("xstar %d, old %f+-%f, new %f+-%f\n", plotbin, oldcontent, olderror, newcontent, newerror);
		 
		 plotbin   = hi_wstar->FindBin(Wavg[ww][qq]);
		 oldcontent = hi_wstar->GetBinContent(plotbin);
		 olderror   = hi_wstar->GetBinError(plotbin);
		 //hi_wstar->AddBinContent(plotbin,d2mrat[Wbin][qq]*rc_bgc_tagInc[ww][qq] );
		 hi_wstar->AddBinContent(plotbin,rc_bgc_tagInc[ww][qq] );
		 newcontent = hi_wstar->GetBinContent(plotbin);
		 if(oldcontent != 0)
		   {
		     newerror = newcontent*sqrtf(olderror*olderror/(oldcontent*oldcontent) 
						 + 1/rc_bgc_tagInc[ww][qq]); 
		     //+ d2mratPerErr[Wbin][qq]*d2mratPerErr[Wbin][qq]);
		   }
		 else newerror = sqrtf(1/rc_bgc_tagInc[ww][qq]);//+ d2mratPerErr[Wbin][qq]*d2mratPerErr[Wbin][qq]);
		 hi_wstar->SetBinError(plotbin,newerror);
		 
		 //printf("wstar %d, old %f+-%f, new %f+-%f\n", plotbin, oldcontent, olderror, newcontent, newerror);
	       }
	     */


 /* else
   {
     for (int qq=minbin;qq<maxbin;qq++)
       {
	 for (int ww=0;ww<numWbins;ww++)
	   {
	     if (rc_untagInc[ww][qq]!=0 
		 && rc_bgc_tagInc[ww][qq]!=0
		 && Wavg[ww][qq] > 2.0
		 )
	       { 
		 F2n_F2p = double(rc_bgc_tagInc[ww][qq])*dScaler*F2dF2pratio(Wavg[ww][qq],Q2avg[ww][qq])/double(rc_untagInc[ww][qq]);
		 //printf("F2d/F2p = %f, F2n/F2p = %f\n", F2d/F2p, F2n_F2p);
		 h2_ratVx->Fill(xavg[ww][qq],F2n_F2p, 1.0);
	       }
	   }
       }
     
     prof_ratVx = h2_ratVx->ProfileX();
     prof_ratVx->Draw();
     can[cc]->cd(pad+1);
     h2_ratVx->Draw("colz");

     }*/

////////////////////////////////////////////////////////////////
//these things are for making a graph of sigma_n/sigma_d vs p_spec
//this is quite old and reallllly inefficient.  if i do it again it 
//should be redone from scratch, using newer arrays that replace these ones

 Double_t ps[npsbins],noverd[npsbins],enod[npsbins], eps[npsbins];
 Int_t totalD[npsbins], totalN[npsbins], psbin;
 Double_t pen, ped, ndrat;

 Double_t noverd2d[numWbins][npsbins];
 Int_t tD2d[numWbins][npsbins], tN2d[numWbins][npsbins];
 
 for(jj = 0; jj < npsbins; jj++)
   {
     for (int ww=0;ww<numWbins;ww++)
       {
	 tD2d[ww][jj] = 0;
	 tN2d[ww][jj] = 0;
       }
     totalD[jj] = 0;
     totalN[jj] = 0;
   }

for (int qq=0;qq<numQ2bins;qq++)
  {
    for (int ww=0;ww<numWbins;ww++)
      {
	//Determine ps bin:
	for (bin=0;bin<npsbins;bin++)
	  {
	    if ((psavg[ww][qq] > psmin[bin]) && (psavg[ww][qq] <= psmax[bin]))
	      {
		psbin=bin;
	      }
	  }
	if(psbin>=0 && psbin<npsbins)
	  {
	    totalN[psbin] += tagIncCount[ww][qq];
	    totalD[psbin] += inclusiveCount[ww][qq];
	    tN2d[ww][psbin] += tagIncCount[ww][qq];
	    tD2d[ww][psbin] += inclusiveCount[ww][qq];
	  }
      }
  }

 for (int t=0;t<npsbins;t++)
   {
     for (int ww=0;ww<numWbins;ww++)
       {
	 ndrat = (double)tN2d[ww][t]/(double)tD2d[ww][t];
	 if(ndrat > 0 && ndrat < 0.03)
	   {
	     h2_ps_w_sig->Fill((Wmax[ww]+Wmin[ww])/2,
			       (psmin[t] + psmax[t])/2, ndrat
			       );
	   }
       }
   }
 gStyle->SetPalette(1);
//can[3]->cd(5); h2_ps_w_sig->SetContour(20);h2_ps_w_sig->Draw("colz");

  for (int t=0;t<npsbins;t++)
    {
      ps[t] = (psmin[t] + psmax[t])/2;
      if(totalD[t]>0) noverd[t] = double(totalN[t])/double(totalD[t]);
      else noverd[t] = 0;
      eps[t]    = double(dps)/double(npsbins)/2;
      pen       = sqrtf((double)totalN[t])/double(totalN[t]);
      ped       = sqrtf((double)totalD[t])/double(totalD[t]);
      enod[t]   = noverd[t]*sqrtf(pen*pen + ped*ped);
      //printf("%d %f %f %f %f\n",t,ps[t],noverd[t],eps[t],enod[t]);
					 
    }

TGraphErrors* gr1 = new TGraphErrors(npsbins,ps,noverd,eps,enod);
 gr1->SetTitle("all W");
gr1->SetMarkerStyle(22);
gr1->SetMarkerColor(kBlue);
gr1->GetXaxis()->SetTitle("p_{spec} (MeV/c)");
gr1->GetYaxis()->SetTitle("#sigma_{n}/#sigma_{d}");
gr1->GetXaxis()->CenterTitle();
gr1->GetYaxis()->CenterTitle();
//gr1->SetMinimum(0.0);
//gr1->SetMaximum(0.04);
//can[0]->cd(1);gr1->Draw("AP");

for(jj = 0; jj < npsbins; jj++)
   {
     totalD[jj] = 0;
     totalN[jj] = 0;
   }

for (int qq=0;qq<numQ2bins;qq++)
  {
    for (int ww=0;ww<numWbins;ww++)
      {
	//Determine ps bin:
	for (bin=0;bin<npsbins;bin++)
	  {
	    if ((psavg[ww][qq] > psmin[bin]) && (psavg[ww][qq] <= psmax[bin]))
	      {
		psbin=bin;
	      }
	  }
	if(psbin>=0 && psbin<npsbins && Wmin[ww] > 2.0)
	  {
	    totalN[psbin] += tagIncCount[ww][qq];
	    totalD[psbin] += inclusiveCount[ww][qq];
	  }
      }
  }

  for (int t=0;t<npsbins;t++)
    {
      ps[t] = (psmin[t] + psmax[t])/2;
      if(totalD[t]>0) noverd[t] = double(totalN[t])/double(totalD[t]);
      else noverd[t] = 0;
      eps[t]    = double(dps)/double(npsbins)/2;
      pen       = sqrtf((double)totalN[t])/double(totalN[t]);
      ped       = sqrtf((double)totalD[t])/double(totalD[t]);
      enod[t]   = noverd[t]*sqrtf(pen*pen + ped*ped);
      //printf("%d %f %f %f %f\n",t,ps[t],noverd[t],eps[t],enod[t]);
					 
    }
TGraphErrors* gr2 = new TGraphErrors(npsbins,ps,noverd,eps,enod);
 gr2->SetTitle("W > 2.0 GeV");
gr2->SetMarkerStyle(22);
gr2->SetMarkerColor(kBlue);
gr2->GetXaxis()->SetTitle("p_{spec} (MeV/c)");
gr2->GetYaxis()->SetTitle("#sigma_{n}/#sigma_{d}");
gr2->GetXaxis()->CenterTitle();
gr2->GetYaxis()->CenterTitle();
 //can[0]->cd(2);gr2->Draw("AP");

for(jj = 0; jj < npsbins; jj++)
   {
     totalD[jj] = 0;
     totalN[jj] = 0;
   }

for (int qq=0;qq<numQ2bins;qq++)
  {
    for (int ww=0;ww<numWbins;ww++)
      {
	//Determine ps bin:
	for (bin=0;bin<npsbins;bin++)
	  {
	    if ((psavg[ww][qq] > psmin[bin]) && (psavg[ww][qq] <= psmax[bin]))
	      {
		psbin=bin;
	      }
	  }
	if(psbin>=0 && psbin<npsbins && Wmin[ww] < 2.0 && Wmin[ww] > 1.07)
	  {
	    totalN[psbin] += tagIncCount[ww][qq];
	    totalD[psbin] += inclusiveCount[ww][qq];
	  }
      }
  }

  for (int t=0;t<npsbins;t++)
    {
      ps[t] = (psmin[t] + psmax[t])/2;
      if(totalD[t]>0) noverd[t] = double(totalN[t])/double(totalD[t]);
      else noverd[t] = 0;
      eps[t]    = double(dps)/double(npsbins)/2;
      pen       = sqrtf((double)totalN[t])/double(totalN[t]);
      ped       = sqrtf((double)totalD[t])/double(totalD[t]);
      enod[t]   = noverd[t]*sqrtf(pen*pen + ped*ped);
      //printf("%d %f %f %f %f\n",t,ps[t],noverd[t],eps[t],enod[t]);
					 
    }
TGraphErrors* gr3 = new TGraphErrors(npsbins,ps,noverd,eps,enod);
 gr3->SetTitle("1.07 GeV < W < 2.0 GeV");
gr3->SetMarkerStyle(22);
gr3->SetMarkerColor(kBlue);
gr3->GetXaxis()->SetTitle("p_{spec} (MeV/c)");
gr3->GetYaxis()->SetTitle("#sigma_{n}/#sigma_{d}");
gr3->GetXaxis()->CenterTitle();
gr3->GetYaxis()->CenterTitle();
 //can[0]->cd(3);gr3->Draw("AP");

for(jj = 0; jj < npsbins; jj++)
   {
     totalD[jj] = 0;
     totalN[jj] = 0;
   }

for (int qq=0;qq<numQ2bins;qq++)
  {
    for (int ww=0;ww<numWbins;ww++)
      {
	//Determine ps bin:
	for (bin=0;bin<npsbins;bin++)
	  {
	    if ((psavg[ww][qq] > psmin[bin]) && (psavg[ww][qq] <= psmax[bin]))
	      {
		psbin=bin;
	      }
	  }
	if(psbin>=0 && psbin<npsbins && Wmin[ww] < 1.07)
	  {
	    totalN[psbin] += tagIncCount[ww][qq];
	    totalD[psbin] += inclusiveCount[ww][qq];
	  }
      }
  }

  for (int t=0;t<npsbins;t++)
    {
      ps[t] = (psmin[t] + psmax[t])/2;
      if(totalD[t]>0) noverd[t] = double(totalN[t])/double(totalD[t]);
      else noverd[t] = 0;
      eps[t]    = double(dps)/double(npsbins)/2;
      pen       = sqrtf((double)totalN[t])/double(totalN[t]);
      ped       = sqrtf((double)totalD[t])/double(totalD[t]);
      enod[t]   = noverd[t]*sqrtf(pen*pen + ped*ped);
      //printf("%d %f %f %f %f\n",t,ps[t],noverd[t],eps[t],enod[t]);
					 
    }
TGraphErrors* gr4 = new TGraphErrors(npsbins,ps,noverd,eps,enod);
 gr4->SetTitle("W < 1.07 GeV");
gr4->SetMarkerStyle(22);
gr4->SetMarkerColor(kBlue);
gr4->GetXaxis()->SetTitle("p_{spec} (MeV/c)");
gr4->GetYaxis()->SetTitle("#sigma_{n}/#sigma_{d}");
gr4->GetXaxis()->CenterTitle();
gr4->GetYaxis()->CenterTitle();
 //can[0]->cd(4);gr4->Draw("AP");


//////////this stuff was to check slava's mom correction when he was gone
/*
can[2]->cd(1); hi_wstar->Draw();
hi_wstar->Fit(egau[1],"R");

 for(ii=0;ii<6;ii++)
   {
     can[3]->cd(ii+1);
     hwb4corr[ii]->Draw();
     hwb4corr[ii]->Fit(egau[ii+2],"R");
   }
 for(ii=0;ii<6;ii++)
   {
     can[4]->cd(ii+1);
     hwafcorr[ii]->Draw();
     hwafcorr[ii]->Fit(egau[ii+7],"R");
   }
 for(ii=0;ii<3;ii++)
   {
     can[5]->cd(ii+1);
     //h2_dpVphisec[ii]->Draw("colz");
     prof[ii] = h2_dpVphisec[ii]->ProfileX();
     prof[ii]->SetMinimum(-0.02);
     prof[ii]->SetMaximum(0.02);
     prof[ii]->Fit("pol1", "", "", -20, 20);
     can[5]->cd(ii+4);
     //h2_dpVphisec[ii+6]->Draw("colz");
     prof[ii+4] = h2_dpVphisec[ii+6]->ProfileX();
     prof[ii+4]->SetMinimum(-0.02);
     prof[ii+4]->SetMaximum(0.02);
     prof[ii+4]->Fit("pol1", "", "", -20, 20);
     
   }
jj=3;
 for(ii=0;ii<3;ii++)
   {
     can[6]->cd(ii+1);
     //h2_dpVphisec[jj]->Draw("colz");
     prof[jj] = h2_dpVphisec[jj]->ProfileX();
     prof[jj]->SetMinimum(-0.02);
     prof[jj]->SetMaximum(0.02);
     prof[jj]->Fit("pol1", "", "", -20, 20);
     can[6]->cd(ii+4);
     //h2_dpVphisec[jj+6]->Draw("colz");
     prof[jj+6] = h2_dpVphisec[jj+6]->ProfileX();
     prof[jj+6]->SetMinimum(-0.02);
     prof[jj+6]->SetMaximum(0.02);
     prof[jj+6]->Fit("pol1", "", "", -20, 20);
     jj++;
   }

*/