// wc_res.C cbc 2003/04/02
//   originally resolution.C
// plots momentum, theta, phi, z0 resolution
//
// root wc_res.C analdir=$DATADIR/root/ 3060-3080
// root wc_res.C +p ntuple=flr 3800-3810
//
{
  bool p=false, t=false, f=false, z=false;
  gROOT->ProcessLine(".x lr.C");
  if (!(p|t|f|z)) { p=t=f=z=true; }

  gStyle->SetOptStat(10);
  gStyle->SetOptFit(6);
  gStyle->SetStatW(.25);
  gStyle->SetStatH(.06);
  float mean, rms;

  // elastic kinematics (momentum, theta)
  TString en(".85/(1+2*.85/.938*(sin(twl/57.3/2)*sin(twl/57.3/2)))");
  TString tp("57.3*asin(("+en+")*sin(twl/57.3)/sqrt((.85-("+en
	     +"))*(.85-("+en+"))+2*.938*(.85-("+en+"))))");

  // momemtum resolution
  if (p) {
    TPad* pp=InitCanvas("recon_p","Reconstructed Momentum, runs "
			+runlist[0],3,2);
    TString cp("0<twl&twl<90 & 0<pwl&pwl<1.2"
	       " & xwl<100000&xwr<100000 & qwr==1&qwl==-1");
    pp->cd(1);
    chain[0]->Draw("pwl:twl>>p_th", cp);
    ((TH1*)gPad->GetPrimitive("p_th"))
      ->SetTitle("p_{left} : #theta_{left}");
    chain[0]->Draw(en+":twl","","same");
    pp->cd(2);
    chain[0]->Draw("pwl-"+en+">>p_res(100,-.15,.15)", cp);
    mean = p_res->GetMean();
    rms = p_res->GetRMS();
    p_res->Fit("gaus","","",mean-1.5*rms, mean+1.5*rms);
    p_res->GetFunction("gaus")->SetLineColor(kBlue);
    p_res->SetTitle("p_{recon} - p_{calc}");
 
    pp->cd(4);
    chain[0]->Draw("pwl-"+en+":pwl", cp);
    TProfile *pf = new TProfile("pf","pf",14,.2,.9, -.6,.4);
 
    pp->cd(5);    
    chain[0]->Draw("pwl-"+en+":pwl >> pf", cp);
    pf->Fit("pol1","same");
    
    pp->cd(3);
    chain[0]->Draw("pwl*0.34+0.5-"+en+">>p_res2(100,-.15,.15)", cp);
    mean = p_res2->GetMean();
    rms = p_res2->GetRMS();
    p_res2->Fit("gaus","","",mean-1.5*rms, mean+1.5*rms);
    p_res2->GetFunction("gaus")->SetLineColor(kBlue);
    p_res2->SetTitle("p_{recon} - p_{calc}");

    pp->cd(6);
    chain[0]->Draw("(pwl*0.34+0.5-"+en+")*100.0/"+en+">>p_res4(100,-15,15)", cp);
    mean = p_res4->GetMean();
    rms = p_res4->GetRMS();
    p_res4->Fit("gaus","","",mean-1.5*rms, mean+1.5*rms);
    p_res4->GetFunction("gaus")->SetLineColor(kBlue);
    p_res4->SetTitle("p_{recon} - p_{calc}");


  }

}