void plot_MM (void)
{
// 
// Read in ascii files for MM 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];
   Int_t j,jj;
   
   #define npts 120;

   // open valeri's file
   
   FILE *in;
   sprintf (filename,"vpk_511011.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 = 0.55;

   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);

   /*c1->cd(1);  
   c1_1->SetGridx();
   c1_1->SetGridy();
   c1_1->SetBorderMode(0);
   c1_1->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 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;    // 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;    // 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=20*sqrt(2*3.14159)*sig;       // gauss normalization   
   pulse_gaus->SetParameters(mu,sigma,gain,xshift,gnorm,mean,sig);


   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);

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

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

   valmKp->DrawClone("Ap");
   // valmKp->Fit("pulse_gaus","","",xmin,xmax);

   // open gagik's file
   
   FILE *in;
   sprintf (filename,"../gagik/g11finalGraphs/GraphSignal_04.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);
   // gagikmKp->Draw("Ap");
   // gagikmKp->Draw("psame");


   // 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.6;

   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 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);
   // toymcmKp->Draw("Ap");
   // toymcmKp->Draw("psame");

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

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

   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);

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

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

   toymcmKp->DrawClone("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.010;       // gauss sig  
   Double_t gnorm=40*sqrt(2*3.14159)*sig;       // gauss normalization  
   pulse_gaus->SetParameters(mu,sigma,gain,xshift,gnorm,mean,sig);
   // toymcmKp->Fit("pulse_gaus","","",xmin,xmax);

   // try some fits

   // plot pulse function

   Float_t xmin=1.435;
   Float_t xmax=1.7;


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

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

   gagikmKp->Draw("Ap");
   gagikmKp->Fit("pulse","","",xmin,xmax);

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

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

   gagikmKp->DrawClone("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.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#pm%.1f\n",Signif,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);

   c2->cd(1);  
   c2_1->SetGridx();
   c2_1->SetGridy();
   c2_1->SetBorderMode(0);
   c2_1->SetFillColor(0);
   valmKp->Draw("Ap");

   c2->cd(2);  
   c2_2->SetGridx();
   c2_2->SetGridy();
   c2_2->SetBorderMode(0);
   c2_2->SetFillColor(0);
   gagikmKp->DrawClone("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.005;       // gauss sig  
   Double_t gnorm=40*sqrt(2*3.14159)*sig;       // gauss normalization  
   pulse_gaus->SetParameters(mu,sigma,gain,xshift,gnorm,mean,sig);

   // fix to toy mc parameters

   // pulse_gaus->FixParameter(0,mu_toy);
   // pulse_gaus->FixParameter(1,sigma_toy);
   // pulse_gaus->FixParameter(2,gain_toy);
   // pulse_gaus->FixParameter(3,xshift_toy);

   // fix to valeri's parameters

   pulse_gaus->FixParameter(0,mu_val);
   pulse_gaus->FixParameter(1,sigma_val);
   // pulse_gaus->FixParameter(2,gain_val);
   pulse_gaus->FixParameter(3,xshift_val);

   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);

   sprintf(string,"Signal=%.1f#pm%.1f\n",signal,sigerr);
   t1 = new TLatex(0.15,0.96,string);
   t1->SetTextColor(1);
   t1->SetTextSize(0.04);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"Signif=%.1f\n",Signif);
   t1 = new TLatex(0.15,0.91,string);
   t1->SetTextColor(1);
   t1->SetTextSize(0.04);
   t1->SetNDC();
   t1->Draw();

   c2->cd(3);  
   c2_3->SetGridx();
   c2_3->SetGridy();
   c2_3->SetBorderMode(0);
   c2_3->SetFillColor(0);
   gagikmKp->Fit("pulse","","",xmin,xmax);
   gagikmKp->DrawClone("Ap");

   c2->cd(4);  
   c2_4->SetGridx();
   c2_4->SetGridy();
   c2_4->SetBorderMode(0);
   c2_4->SetFillColor(0);

   mu_gagik = pulse_gaus->GetParameter(0);
   sigma_gagik = pulse_gaus->GetParameter(1);
   gain_gagik = pulse_gaus->GetParameter(2);
   xshift_gagik = pulse_gaus->GetParameter(3);
   printf (" Gagik parameters mu=%f, sigma=%f, gain=%f, xshift=%f\n",mu_gagik,sigma_gagik,gain_gagik,xshift_gagik);

   // fix to gagik's parameters

   pulse_gaus->FixParameter(0,mu_gagik);
   pulse_gaus->FixParameter(1,sigma_gagik);
   // pulse_gaus->FixParameter(2,gain_gagik);
   pulse_gaus->FixParameter(3,xshift_gagik);

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

   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);

   sprintf(string,"Signal=%.1f#pm%.1f\n",signal,sigerr);
   t1 = new TLatex(0.15,0.96,string);
   t1->SetTextColor(1);
   t1->SetTextSize(0.04);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"Signif=%.1f\n",Signif);
   t1 = new TLatex(0.15,0.91,string);
   t1->SetTextColor(1);
   t1->SetTextSize(0.04);
   t1->SetNDC();
   t1->Draw();


   sprintf(filename,"plot_MM_c1.eps");
   c1->SaveAs(filename);
   sprintf(filename,"plot_MM_c2.eps");
   c2->SaveAs(filename);
   sprintf(filename,"plot_MM_c3.eps");
   c3->SaveAs(filename);
   sprintf(filename,"plot_MM_c4.eps");
   c4->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;
}