void plot_MM2 (void)
{
// 
// Read in ascii files for MM2 phi-theta mass distributions Valeri and Moskov
//      7/31/09
//

#include <TMath.h>
#include <TRandom.h>
#include <TTree.h>

gROOT->Reset();

gStyle->SetPalette(1,0);
gStyle->SetOptStat(kFALSE);;
// gStyle->SetOptStat(kTRUE);
// gStyle->SetOptFit(kFALSE);
gStyle->SetOptFit(1111);
gStyle->SetPadRightMargin(0.15);
gStyle->SetPadLeftMargin(0.15);
gStyle->SetPadBottomMargin(0.15);
// gStyle->SetFillColor(0);
//

   char string[256];
   char filename[132];
   char dataname[80];
   Int_t j,jj;
   
   #define npts 120;

   // decide on data set

   // sprintf (dataname,"gagik");
   sprintf (dataname,"val");
   // sprintf (dataname,"toy");
   printf ("Data selection=%s\n",dataname);

   // open valeri's file

   FILE *in;
   sprintf (filename,"vpk_510011.dat");
   in = fopen (filename,"r");
   
   // plot measured mass distribution

   Float_t val_counts[npts];
   Float_t val_xerrors[npts];
   Float_t val_yerrors[npts];
   Float_t val_mKp[npts];
   Float_t val_delta=0.005;
   Float_t val_offset=1.4;
   j=0;

   // normalize valeri's data by factor norm
   Float_t norm = 1.;

   while(fgets(string, 256, in)!=NULL) {
      // printf ("string=%s",string);
      Float_t x, y, yerr;
      sscanf (string,"%f %f %f",&x,&y,&yerr);
      // printf ("j=%d y=%f, y=%f, yerr=%f\n",j,x,y,yerr);
      
      val_mKp[j]=x;
      val_xerrors[j] = 0.0001;
      val_counts[j] = y*norm;
      val_yerrors[j] = yerr*norm;
      printf ("j=%d, val_mKp=%f, val_counts=%f, val_yerrors=%f\n",j,val_mKp[j],val_counts[j],val_yerrors[j]);
      j++;
      }
   val_mKp[0] = 1.3;      // to plot from 1.4
   val_mKp[npts-1] = 2.2; // to plot up to 2.0
   fclose (in);

   TGraph *valmKp = new TGraphErrors(npts,val_mKp,val_counts,val_xerrors,val_yerrors);

   TCanvas *c1 = new TCanvas("c1","c1 Plot_MM Canvas",200,10,700,700);
   c1->SetBorderMode(0);
   c1->SetFillColor(0);
   // c1->Divide(2,2);

   c1->SetGridx();
   c1->SetGridy();
   c1->SetBorderMode(0);
   c1->SetFillColor(0);

   Float_t x1=1.4;
   Float_t x2=2.0;
   Float_t y1=0;
   Float_t y2=140;
   valmKp->SetTitle("");
   valmKp->GetXaxis()->SetRangeUser(x1,x2);
   valmKp->GetYaxis()->SetRangeUser(y1,y2);
   valmKp->GetXaxis()->SetTitleSize(0.05);
   valmKp->GetYaxis()->SetTitleSize(0.05);
   valmKp->GetYaxis()->SetTitleOffset(1.5);
   valmKp->GetXaxis()->SetTitle("m_{pK} (GeV)");
   valmKp->GetYaxis()->SetTitle("Counts");
   valmKp->GetXaxis()->SetNdivisions(5);
   // valmKp->GetYaxis()->SetNdivisions(5);
   valmKp->SetMarkerColor(2);
   valmKp->SetMarkerStyle(21);
   valmKp->SetMarkerSize(0.7);

   // Float_t xmin=1.435;
   Float_t xmin=1.45;
   Float_t xmax=1.7;

   TF1 *pulse = new TF1("pulse",pulse,xmin,xmax,4);
   Double_t mu=6;    // unis of 1/x
   Double_t sigma=0.04;     // units of x  
   Double_t gain=30;    // original 
   Double_t gain=300;    // arbitrary
   Double_t xshift=1.47;  // x translation  
   pulse->SetParameter(0,mu);    
   pulse->SetParameter(1,sigma);    
   pulse->SetParameter(2,gain);    
   pulse->SetParameter(3,xshift);
   // pulse->DrawCopy("same");

   TF1 *pulse_gaus = new TF1("pulse_gaus",pulse_gaus,xmin,xmax,7);
   Double_t mu=6;    // unis of 1/x
   Double_t sigma=0.04;     // units of x  
   Double_t gain=30;    // original
   Double_t xshift=1.47;  // x translation 
   Double_t mean=1.54;       // gauss mean  
   Double_t sig=0.010;       // gauss sig   
   Double_t gnorm=20*sqrt(2*3.14159)*sig;       // gauss normalization   
   pulse_gaus->SetParameters(mu,sigma,gain,xshift,gnorm,mean,sig);


  if (!strcmp(dataname,"val")) {
      valmKp->Draw("Ap");
      valmKp->Fit("pulse","","",xmin,xmax); 
   }
   mu_val = pulse->GetParameter(0);
   sigma_val = pulse->GetParameter(1);
   gain_val = pulse->GetParameter(2);
   xshift_val = pulse->GetParameter(3);

   // open gagik's file
   
   FILE *in;
   sprintf (filename,"../gagik2/g11finalGraphs/GraphSignal_11.data");
   in = fopen (filename,"r");
   
   // plot measured mass distribution

   Float_t gagik_counts[npts];
   Float_t gagik_xerrors[npts];
   Float_t gagik_yerrors[npts];
   Float_t gagik_mKp[npts];
   Float_t gagik_delta=0.005;
   Float_t gagik_offset=1.4;
   j=0;
   Float_t binsize;

   while(fgets(string, 256, in)!=NULL) {
      // printf ("string=%s",string);
      Float_t x, y, bin, yerr;
      binsize = bin;
      sscanf (string,"%f %f %f %f",&x,&y,&bin,&yerr);
      // printf ("j=%d x=%f, y=%f, bin=%f, yerr=%f\n",j,x,y,bin,yerr);
      
      gagik_mKp[j]=x;
      gagik_xerrors[j] = 0.0001;
      gagik_counts[j] = y;
      gagik_yerrors[j] = yerr;
      // printf ("j=%d, gagik_mKp=%f, gagik_counts=%f, gagik_yerrors=%f\n",j,gagik_mKp[j],gagik_counts[j],gagik_yerrors[j]);
      j++;
      }
   gagik_mKp[0] = 1.3;      // to plot from 1.4
   gagik_mKp[99] = 2.2; // to plot up to 2.0
   fclose (in);

   TGraph *gagikmKp = new TGraphErrors(npts,gagik_mKp,gagik_counts,gagik_xerrors,gagik_yerrors);

   // now get data from Moskov

   gagikmKp->SetTitle("");
   gagikmKp->GetXaxis()->SetRangeUser(x1,x2);
   gagikmKp->GetYaxis()->SetRangeUser(y1,y2);
   gagikmKp->GetXaxis()->SetTitleSize(0.05);
   gagikmKp->GetYaxis()->SetTitleSize(0.05);
   gagikmKp->GetYaxis()->SetTitleOffset(1.5);
   gagikmKp->GetXaxis()->SetTitle("m_{pK} (GeV)");
   gagikmKp->GetYaxis()->SetTitle("Counts");
   gagikmKp->GetXaxis()->SetNdivisions(5);
   gagikmKp->SetMarkerColor(4);
   gagikmKp->SetMarkerStyle(20);

   // open toy MC's file
   
   FILE *in;
   sprintf (filename,"../calcs/tslope4/ToyMC_.dat_1500-2550_tcut45.dat");
   in = fopen (filename,"r");
   
   // plot measured mass distribution

   Float_t toymc_counts[npts];
   Float_t toymc_xerrors[npts];
   Float_t toymc_yerrors[npts];
   Float_t toymc_mKp[npts];
   Float_t toymc_delta=0.005;
   Float_t toymc_offset=1.4;
   j=0;
   Float_t binsize;
   Double_t mc_norm=1.1;

   while(fgets(string, 256, in)!=NULL) {
      // printf ("string=%s",string);
      Float_t x, y, yerr;
      binsize = bin;
      sscanf (string,"%f %f %f %f",&x,&y,&yerr);
      // printf ("j=%d x=%f, y=%f, yerr=%f\n",j,x,y,bin,yerr);
      
      toymc_mKp[j]=x;
      toymc_xerrors[j] = 0.0001;
      toymc_counts[j] = mc_norm*y;
      toymc_yerrors[j] = mc_norm*yerr;
      // printf ("j=%d, toymc_mKp=%f, toymc_counts=%f, toymc_yerrors=%f\n",j,toymc_mKp[j],toymc_counts[j],toymc_yerrors[j]);
      j++;
      }
   toymc_mKp[0] = 1.3;      // to plot from 1.4
   toymc_mKp[99] = 2.2; // to plot up to 2.0
   fclose (in);

   TGraph *toymcmKp = new TGraphErrors(npts,toymc_mKp,toymc_counts,toymc_xerrors,toymc_yerrors);

   // Float_t xmin=1.435;
   Float_t xmin=1.45;
   Float_t xmax=1.7;

   Float_t x1=1.4;
   Float_t x2=2.0;
   Float_t y1=0;
   Float_t y2=140;
   toymcmKp->SetTitle("");
   toymcmKp->GetXaxis()->SetRangeUser(x1,x2);
   toymcmKp->GetYaxis()->SetRangeUser(y1,y2);
   toymcmKp->GetXaxis()->SetTitleSize(0.05);
   toymcmKp->GetYaxis()->SetTitleSize(0.05);
   toymcmKp->GetYaxis()->SetTitleOffset(1.5);
   toymcmKp->GetXaxis()->SetTitle("m_{pK} (GeV)");
   toymcmKp->GetYaxis()->SetTitle("Counts");
   toymcmKp->GetXaxis()->SetNdivisions(5);
   toymcmKp->SetMarkerColor(3);
   toymcmKp->SetMarkerStyle(21);


   if (!strcmp(dataname,"toy")) {
      toymcmKp->Draw("Ap");
      toymcmKp->Fit("pulse","","",xmin,xmax);
   } 
   mu_toy = pulse->GetParameter(0);
   sigma_toy = pulse->GetParameter(1);
   gain_toy = pulse->GetParameter(2);
   xshift_toy = pulse->GetParameter(3);
   printf (" toy par mu=%f, sigma=%f, gain=%f, xshift=%f\n",mu_toy,sigma_toy,gain_toy,xshift_toy);

   // Float_t xmin=1.435;
   Float_t xmin=1.45;
   Float_t xmax=1.7;

   if (!strcmp(dataname,"gagik")) {
      gagikmKp->Draw("Ap");
      gagikmKp->Fit("pulse","","",xmin,xmax);
   } 
   mu_gagik = pulse->GetParameter(0);
   sigma_gagik = pulse->GetParameter(1);
   gain_gagik = pulse->GetParameter(2);
   xshift_gagik = pulse->GetParameter(3);
   printf (" gagik par mu=%f, sigma=%f, gain=%f, xshift=%f\n",mu_gagik,sigma_gagik,gain_gagik,xshift_gagik);

   Double_t mu=6;    // unis of 1/x
   Double_t sigma=0.04;     // units of x  
   Double_t gain=30;    // arbitrary
   Double_t xshift=1.47;  // x translation
   Double_t mean=1.54;       // gauss mean  
   Double_t sig=0.010;       // gauss sig  
   Double_t gnorm=40*sqrt(2*3.14159)*sig;       // gauss normalization  
   pulse_gaus->SetParameters(mu,sigma,gain,xshift,gnorm,mean,sig);
   // gagikmKp->Fit("pulse_gaus","","",xmin,xmax); 

   mu = pulse_gaus->GetParameter(0);
   sigma = pulse_gaus->GetParameter(1);
   gain = pulse_gaus->GetParameter(2);
   xshift = pulse_gaus->GetParameter(3);
   gnorm = fabs(pulse_gaus->GetParameter(4));
   mean = pulse_gaus->GetParameter(5);
   sig = fabs(pulse_gaus->GetParameter(6));
   printf ("mu=%f, sigma=%f, gain=%f, xshift=%f, gnorm=%f, mean=%f, sig=%f, binsize=%f\n",mu,sigma,gain,xshift,gnorm,mean,sig,binsize);
   Double_t signal = (1-0.0455)*gnorm/binsize;  // area outside 2sigma
   gnormerr = pulse_gaus->GetParError(4);
   gainerr = pulse_gaus->GetParError(2);
   Double_t sigerr = (1-0.0455)*gnormerr/binsize;
   Double_t back = pulse_gaus->Integral(mean-2*sig,mean+2*sig)/binsize - signal;
   Double_t backerr = back*gainerr/gain;
   if (backerr < sqrt(back)) backerr = sqrt(back);
   Double_t Signif = signal/backerr;

   
   /*sprintf(string,"signal=%.1f#pm%.1f\n",signal,sigerr);
   t1 = new TLatex(0.5,0.6,string);
   t1->SetTextColor(1);
   t1->SetTextSize(0.04);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"back=%.1f#pm%.1f\n",back,backerr);
   t1 = new TLatex(0.5,0.5,string);
   t1->SetTextColor(1);
   t1->SetTextSize(0.04);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"Signif=%.1f\n",Signif);
   t1 = new TLatex(0.5,0.4,string);
   t1->SetTextColor(1);
   t1->SetTextSize(0.04);
   t1->SetNDC();
   t1->Draw();*/


   TCanvas *c2 = new TCanvas("c2","c2 Plot_MM Canvas",200,10,700,700);
   c2->SetBorderMode(0);
   c2->SetFillColor(0);
   // c2->Divide(2,2);

   c2->SetGridx();
   c2->SetGridy();
   c2->SetBorderMode(0);
   c2->SetFillColor(0);

   // toymcmKp->Draw("Ap");

   gagikmKp->Draw("Ap");
   Double_t mu=6;    // unis of 1/x
   Double_t sigma=0.04;     // units of x  
   Double_t gain=30;    // arbitrary
   Double_t xshift=1.47;  // x translation
   Double_t mean=1.54;       // gauss mean  
   Double_t sig=0.007;       // gauss sig  
   Double_t gnorm=40*sqrt(2*3.14159)*sig;       // gauss normalization  
   pulse_gaus->SetParameters(mu,sigma,gain,xshift,gnorm,mean,sig);


   Double_t backnorm;
   if (!strcmp(dataname,"toy")) {
      pulse_gaus->FixParameter(0,mu_toy);
      pulse_gaus->FixParameter(1,sigma_toy);
      pulse_gaus->FixParameter(3,xshift_toy);
      pulse_gaus->FixParameter(2,gain_toy);
      backnorm = gain_toy;
      }
   elseif (!strcmp(dataname,"val")) {
      pulse_gaus->FixParameter(0,mu_val);
      pulse_gaus->FixParameter(1,sigma_val);
      pulse_gaus->FixParameter(3,xshift_val);
      backnorm = gain_val;
      }
   elseif (!strcmp(dataname,"gagik")) {
      pulse_gaus->FixParameter(0,mu_gagik);
      pulse_gaus->FixParameter(1,sigma_gagik);
      pulse_gaus->FixParameter(3,xshift_gagik);
      backnorm = gain_gagik;
      }
   else {
      printf ("*** plot_MM2 - illegal selection=%s\n",dataname);
      }

   Float_t xmin=1.45;
   Float_t xmax=1.7;

   gagikmKp->Fit("pulse_gaus","","",xmin,xmax); 

   mu = pulse_gaus->GetParameter(0);
   sigma = pulse_gaus->GetParameter(1);
   gain = pulse_gaus->GetParameter(2);
   xshift = pulse_gaus->GetParameter(3);
   gnorm = fabs(pulse_gaus->GetParameter(4));
   mean = pulse_gaus->GetParameter(5);
   sig = fabs(pulse_gaus->GetParameter(6));
   printf ("mu=%f, sigma=%f, gain=%f, xshift=%f, gnorm=%f, mean=%f, sig=%f, binsize=%f\n",mu,sigma,gain,xshift,gnorm,mean,sig,binsize);
   Double_t signal = (1-0.0455)*gnorm/binsize;  // area outside 2sigma
   gnormerr = pulse_gaus->GetParError(4);
   gainerr = pulse_gaus->GetParError(2);
   Double_t sigerr = (1-0.0455)*gnormerr/binsize;
   Double_t back = pulse_gaus->Integral(mean-2*sig,mean+2*sig)/binsize - signal;
   Double_t backerr = back*gainerr/gain;
   if (backerr < sqrt(back)) backerr = sqrt(back);
   Double_t Signif = signal/backerr;
   printf ("signal=%f, gnormerr=%f, gainerr=%f, sigerr=%f, back=%f, backerr=%f, Signif=%f\n",signal,gnormerr,gainerr,sigerr,back,backerr,Signif);

   // compute log likelihood

   Double_t logLB = 0;
   Double_t logLBS = 0;

   // find normalization factor

   Double_t datanorm=0;
   for (j=0;j<npts;j++) {
      Double_t datax = gagik_mKp[j];
      if (datax>xmin && datax<xmax) {
         Double_t datay = gagik_counts[j];
         datanorm = datanorm + datay;
         }
      }
   Double_t funcnorm = pulse->Integral(xmin,xmax)/binsize;
   printf (" datanorm=%f, funcnorm=%f ratio=%f\n",datanorm,funcnorm,datanorm/funcnorm);
      
   for (j=0;j<npts;j++) {
      Double_t datax = gagik_mKp[j];
      if (datax>xmin && datax<xmax) {
         Double_t datay = gagik_counts[j];
         Double_t toty = pulse_gaus->Eval(datax);
         // Double_t backy = pulse->Eval(datax)*gain/backnorm;
         // Double_t backy = pulse->Eval(datax);
         Double_t backy = pulse->Eval(datax)*datanorm/funcnorm;
         printf ("datax=%f, datay=%f, backy=%f, toty=%f\n",datax,datay,backy,toty);
         logLB = logLB -backy + datay*log(backy);
         logLBS = logLBS -toty + datay*log(toty);
         }
      }
   printf ("backnorm=%f, gnorm=%f\n",backnorm,gain);
   Double_t Signif_log = sqrt(2*(logLBS-logLB));
   printf ("Signif_log=%f, logLB =%f, logLBS=%f\n",Signif_log,logLB,logLBS);

   sprintf(string,"Signal=%.1f#pm%.1f\n",signal,sigerr);
   t1 = new TLatex(0.55,0.55,string);
   t1->SetTextColor(1);
   t1->SetTextSize(0.04);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"Back=%.1f#pm%.1f\n",back,backerr);
   t1 = new TLatex(0.55,0.50,string);
   t1->SetTextColor(1);
   t1->SetTextSize(0.04);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"Naive Sig=%.1f\n",Signif);
   t1 = new TLatex(0.55,0.45,string);
   t1->SetTextColor(1);
   t1->SetTextSize(0.04);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"ln(LB)=%.1f\n",logLB);
   t1 = new TLatex(0.55,0.40,string);
   t1->SetTextColor(1);
   t1->SetTextSize(0.04);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"ln(LB+S)=%.1f\n",logLBS);
   t1 = new TLatex(0.55,0.35,string);
   t1->SetTextColor(1);
   t1->SetTextSize(0.04);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"Log Sig=%.1f\n",Signif_log);
   t1 = new TLatex(0.55,0.30,string);
   t1->SetTextColor(1);
   t1->SetTextSize(0.04);
   t1->SetNDC();
   t1->Draw();



   sprintf(filename,"plot_MM2_%s_c1.eps",dataname);
   printf ("file selection =%s\n",filename);
   c1->SaveAs(filename);
   sprintf(filename,"plot_MM2_%s_c2.eps",dataname);
   c2->SaveAs(filename);


}
Double_t pulse (Double_t *x, Double_t *par) {

   Double_t mu=par[0];
   Double_t sigma=par[1];
   Double_t gain=par[2];
   Double_t xshift=par[3];
   Double_t x1=x[0]-xshift;

   char string[256];

   Double_t arg = (mu*sigma*sigma-x1)/(sqrt(2)*sigma);
   Double_t amplitude = (mu/2)* exp(-mu*x1 +(mu*sigma)*(mu*sigma)/2) *  TMath::Erfc(arg);

   amplitude = gain*amplitude;
   
   /*sprintf (string,"amplitude=%f mu=%f sigma=%f arg=%f\n",amplitude,mu,sigma,arg);
   printf ("string=%s",string);*/

   return amplitude;
}
Double_t pulse_gaus (Double_t *x, Double_t *par) {

   Double_t mu=par[0];
   Double_t sigma=par[1];
   Double_t gain=par[2];
   Double_t xshift=par[3];
   Double_t norm=par[4];
   Double_t mean=par[5];
   Double_t sig=par[6];

   Double_t x1=x[0]-xshift;
   Double_t x2=x[0];         // use absolute x for gaussian

   char string[256];

   Double_t arg = (mu*sigma*sigma-x1)/(sqrt(2)*sigma);
   Double_t amplitude = (mu/2)* exp(-mu*x1 +(mu*sigma)*(mu*sigma)/2) *  TMath::Erfc(arg);
   Double_t gaus = norm*exp(-(x2-mean)*(x2-mean)/(2*sig*sig))/(sqrt(2*3.14159)*sig);

   amplitude = gain*amplitude + gaus;
   
   /*sprintf (string,"amplitude=%f mu=%f sigma=%f arg=%f\n",amplitude,mu,sigma,arg);
   printf ("string=%s",string);*/

   return amplitude;
}