// MACRO FOR DRAWING NORMALIZED DISTRIBUTIONS.
// FIRST DEFINES THE CANVAS WITH  FOUR PADS 
// OPEN 14 FILES 
// READS NORMALISATION FACTORS FOR THE FILES AND THE CUT VARIABLE FROM FILE NORM.XX
// CREATES THE CUT
// FILLS THE HISTOGRAMS WITH DATA FROM THE FILES (WITH CUT AND NORMALISATION)
// BY SIGVE HAUG SEPTEMBER 2001.


// YOU HAVE TO SET DMAX !!!!!


#include <fstream>


void DrawTrees(char varname[],char Title[],char cutstring[],float norm[],TTree *t0,TTree *t1,TTree *t2,TTree *t3,TTree *t4,TTree *t5,TTree *t6,TTree *t7,TTree *t8,TTree *t9,TTree *t10,TTree *t11,TTree *t12,TTree *t13,int nbin,float fbin,float lbin)
{
   TH1F *bkg=new TH1F("bkg","Background",nbin,fbin,lbin);
   TH1F *sig=new TH1F("sig","Signal",nbin,fbin,lbin);
   TH1F *zgp=new TH1F("zgp","Signal",nbin,fbin,lbin);
   TH1F *wwe=new TH1F("wwe","Signal",nbin,fbin,lbin);
   TH1F *h=new TH1F("h","h",nbin,fbin,lbin);
// Add the normalized backgrounds to background histo 
   t0 ->Draw(varname,cutstring);
   if (h->Integral()>0) sig->Add(h,norm[0]/h->Integral());
   t1 ->Draw(varname,cutstring);   
   if (h->Integral()>0) bkg->Add(h,norm[1]/h->Integral());
   if (h->Integral()>0) zgp->Add(h,norm[1]/h->Integral());   
   if (h->Integral()>0) wwe->Add(h,norm[1]/h->Integral());  
   t2 ->Draw(varname,cutstring);   
   if (h->Integral()>0) bkg->Add(h,norm[2]/h->Integral());
   t3 ->Draw(varname,cutstring);   
   if (h->Integral()>0) bkg->Add(h,norm[3]/h->Integral());
   t4 ->Draw(varname,cutstring);   
   if (h->Integral()>0) bkg->Add(h,norm[4]/h->Integral());
   t5 ->Draw(varname,cutstring);   
   if (h->Integral()>0) bkg->Add(h,norm[5]/h->Integral());
   if (h->Integral()>0) wwe->Add(h,norm[5]/h->Integral());
   t6 ->Draw(varname,cutstring);   
   if (h->Integral()>0) bkg->Add(h,norm[6]/h->Integral());
   t7 ->Draw(varname,cutstring);   
   if (h->Integral()>0) bkg->Add(h,norm[7]/h->Integral());
   t8 ->Draw(varname,cutstring);   
   if (h->Integral()>0) bkg->Add(h,norm[8]/h->Integral());
   t9 ->Draw(varname,cutstring);   
   if (h->Integral()>0) bkg->Add(h,norm[9]/h->Integral());
   t10->Draw(varname,cutstring);   
   if (h->Integral()>0) bkg->Add(h,norm[10]/h->Integral());
   t11->Draw(varname,cutstring);   
   if (h->Integral()>0) bkg->Add(h,norm[11]/h->Integral());
   t12->Draw(varname,cutstring);   
   if (h->Integral()>0) bkg->Add(h,norm[12]/h->Integral());
// Draw bkg with data and scaled signal
   h->SetTitle();  
   h->GetXaxis()->SetTitleOffset(1.5);
   h->GetXaxis()->SetTitle(Title);
   h->GetYaxis()->SetTitle("Events");
   t13->Draw(varname,cutstring,"E");
   bkg->SetFillColor(41);
   bkg->DrawCopy("SAME");
   wwe->SetFillColor(42);
   wwe->DrawCopy("SAME");
   zgp->SetFillColor(43);
   zgp->DrawCopy("SAME");
   t13->Draw(varname,cutstring,"ESAME");
   sig->Scale(10.);
   sig->SetFillColor(2);
   sig->DrawCopy("SAME");
   cout<<"Exp bgd events "<<bkg->Integral()<<" Data "<<h->Integral()<<endl;
}



main(){
   float dmax=12.; //   float dmax=12.5; // =28.;
// Set style which is defined in style.cc
   gROOT->ProcessLine(".x style.cc");
// Open files 
   TFile *f0 =new TFile("/mn/susy/data/delphi/sigh/btag/s115.root","READ");
   TFile *f1 =new TFile("/mn/susy/data/delphi/sigh/btag/zgpy.root","READ");
   TFile *f2 =new TFile("/mn/susy/data/delphi/sigh/btag/bhfo.root","READ");
   TFile *f3 =new TFile("/mn/susy/data/delphi/sigh/btag/kora.root","READ");
   TFile *f4 =new TFile("/mn/susy/data/delphi/sigh/btag/mumu.root","READ");
   TFile *f5 =new TFile("/mn/susy/data/delphi/sigh/btag/wwex.root","READ");
   TFile *f6 =new TFile("/mn/susy/data/delphi/sigh/btag/qqnn.root","READ");
   TFile *f7 =new TFile("/mn/susy/data/delphi/sigh/btag/evqq.root","READ");
   TFile *f8 =new TFile("/mn/susy/data/delphi/sigh/btag/qqee.root","READ");
   TFile *f9 =new TFile("/mn/susy/data/delphi/sigh/btag/qqmm.root","READ");
   TFile *f10=new TFile("/mn/susy/data/delphi/sigh/btag/qqtt.root","READ");
   TFile *f11=new TFile("/mn/susy/data/delphi/sigh/btag/ggm1.root","READ");
   TFile *f12=new TFile("/mn/susy/data/delphi/sigh/btag/ggtt.root","READ");
   TFile *f13=new TFile("/mn/susy/data/delphi/sigh/btag/xeds.root","READ");
// Get the normalisation factors and the discriminant cut value
   ifstream normfile("/mn/susy/data/delphi/sigh/prf206/btag/ef86/norm.2"); 
   float norm[13];
   float nevtorg[13];
   float ntxsecorg[13];
   float dcut,lumo;
   normfile>>dcut; cout <<"DCut value : "<<dcut<<endl; 
   for (int i=0;i<13;i++){
     normfile>>norm[i]; cout<<"Normfactor "<<i<<" "<<norm[i]<<endl;
     normfile>>nevtorg[i];
     normfile>>ntxsecorg[i];
   }
   normfile>>lumo; cout << "Lumo : "<<lumo<<endl;
   normfile.close();
// Create the cut
   char cutstring[20];  
   sprintf(cutstring,"Sigtp > %f",dcut);
   cout << cutstring<<endl;
// Get tree h31 which was a hbook ntuple   
   TTree *t0 = (TTree*)f0 ->Get("h31");
   TTree *t1 = (TTree*)f1 ->Get("h31");
   TTree *t2 = (TTree*)f2 ->Get("h31");
   TTree *t3 = (TTree*)f3 ->Get("h31");
   TTree *t4 = (TTree*)f4 ->Get("h31");
   TTree *t5 = (TTree*)f5 ->Get("h31");
   TTree *t6 = (TTree*)f6 ->Get("h31");
   TTree *t7 = (TTree*)f7 ->Get("h31");
   TTree *t8 = (TTree*)f8 ->Get("h31");
   TTree *t9 = (TTree*)f9 ->Get("h31");
   TTree *t10= (TTree*)f10->Get("h31");
   TTree *t11= (TTree*)f11->Get("h31");
   TTree *t12= (TTree*)f12->Get("h31");
   TTree *t13= (TTree*)f13->Get("h31");




// Create histos
   c1_1->cd();
   char varname[15];
   strcpy(varname,"Mcon");
   strcat(varname," >> h");
   char title[]="Mass / GeV c^{-2}";
   DrawTrees(varname,title,cutstring,norm,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,60,-5.,119.5);

   /* Draw some text 
   TLatex *tex = new TLatex(9.6,9.7,"Data from 2000");
   tex->SetTextSize(0.03);
   tex->SetLineWidth(2);
   tex->Draw();
   TLatex *tex = new TLatex(9.6,8.7,"Analysis at 206 GeV");
   tex->SetTextSize(0.03);
   tex->SetLineWidth(2);
   tex->Draw();
   TLatex *tex = new TLatex(9.6,7.7,"30 % Higgs efficiency");
   tex->SetTextSize(0.03);
   tex->SetLineWidth(2);
   tex->Draw();
   */
// Create histos
   c1_2->cd();
   char varname[15];
   strcpy(varname,"Sigtp");
   strcat(varname," >> h");
   char tit2[]="Discriminating variable";
   DrawTrees(varname,tit2,cutstring,norm,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,40,8.,dmax);



   
// Draw efficiency curves (temporary pasta code)
// Fill tmp histo with fixed bins
   gROOT->ProcessLine(".x style2.cc");
   c2_1->cd();
   int nbin=75;
   TH1F *tmp=new TH1F("tmp","Temporary histo",nbin,-0.005,0.745);
   TH1F *dattmp=new TH1F("dattmp","Temporary histo",nbin,-0.005,0.745);
   float dmin=dcut;
   float x,y;
   for (int i=0;i<100;i++){
     sprintf(cutstring,"Sigtp > %f",dcut);
     t0->Draw("Mcon",cutstring);
     x=t0->GetSelectedRows()/nevtorg[0];
     if (tmp->GetBinContent((int)x)==0){
       t1->Draw("Mcon",cutstring);
       if (nevtorg[1]>0){ 
	 y=t1->GetSelectedRows()*lumo*ntxsecorg[1]/(float)nevtorg[1]);
       }
       tmp->Fill(x,y);
       t2->Draw("Mcon",cutstring);
       if (nevtorg[2]>0){ 
	 y=t2->GetSelectedRows()*lumo*ntxsecorg[2]/(float)nevtorg[2]);
       }
       tmp->Fill(x,y);
       t3->Draw("Mcon",cutstring);
       if (nevtorg[3]>0){ 
	 y=t3->GetSelectedRows()*lumo*ntxsecorg[3]/(float)nevtorg[3]);
       }
       tmp->Fill(x,y);
       t4->Draw("Mcon",cutstring);
       if (nevtorg[4]>0){ 
	 y=t4->GetSelectedRows()*lumo*ntxsecorg[4]/(float)nevtorg[4]);
       }
       tmp->Fill(x,y);
       t5->Draw("Mcon",cutstring);
       if (nevtorg[5]>0){ 
	 y=t5->GetSelectedRows()*lumo*ntxsecorg[5]/(float)nevtorg[5]);
       }
       tmp->Fill(x,y);
       t6->Draw("Mcon",cutstring);
       if (nevtorg[6]>0){ 
	 y=t6->GetSelectedRows()*lumo*ntxsecorg[6]/(float)nevtorg[6]);
       }
       tmp->Fill(x,y);
       t7->Draw("Mcon",cutstring);
       if (nevtorg[7]>0){ 
	 y=t7->GetSelectedRows()*lumo*ntxsecorg[7]/(float)nevtorg[7]);
       }
       tmp->Fill(x,y);
       t8->Draw("Mcon",cutstring);
       if (nevtorg[8]>0){ 
	 y=t8->GetSelectedRows()*lumo*ntxsecorg[8]/(float)nevtorg[8]);
       }
       tmp->Fill(x,y);
       t9->Draw("Mcon",cutstring);
       if (nevtorg[9]>0){ 
	 y=t9->GetSelectedRows()*lumo*ntxsecorg[9]/(float)nevtorg[9]);
       }
       tmp->Fill(x,y);
       t10->Draw("Mcon",cutstring);
       if (nevtorg[10]>0){ 
	 y=t10->GetSelectedRows()*lumo*ntxsecorg[10]/(float)nevtorg[10]);
       }
       tmp->Fill(x,y);
       t11->Draw("Mcon",cutstring);
       if (nevtorg[11]>0){ 
	 y=t11->GetSelectedRows()*lumo*ntxsecorg[11]/(float)nevtorg[11]);
       }
       tmp->Fill(x,y);
       t12->Draw("Mcon",cutstring);
       if (nevtorg[12]>0){ 
	 y=t12->GetSelectedRows()*lumo*ntxsecorg[12]/(float)nevtorg[12]);
       }
       tmp->Fill(x,y);
       t13->Draw("Mcon",cutstring);
       y=t13->GetSelectedRows();
       dattmp->Fill(x,y);
     }
     dcut=dcut+(dmax-dmin)/20.;
   }
// Make array with low bin edges for variable bin size histo  
   int i=0;
   float xbins[75];
   cout<<"NBINS = "<<nbin<<endl;
   for (int ibin=1;ibin<nbin;ibin++){
     if (tmp->GetBinContent(ibin)!=0){
       i++;
       xbins[i]=tmp->GetBinLowEdge(ibin);
     cout<<"BIN "<<i<<" "<<xbins[i]<<" "<<tmp->GetBinContent(ibin)<<endl;   
     }
   } 
// Create and fill variable bin size histo eff
   TH1F *eff=new TH1F("eff","Events vs signal efficiency ",i,xbins);
   TH1F *dat=new TH1F("dat","Events vs signal efficiency ",i,xbins);
   for (int ibin=1;ibin<nbin;ibin++){
     if (tmp->GetBinContent(ibin)!=0){
       eff->Fill(tmp->GetBinCenter(ibin),tmp->GetBinContent(ibin));
       dat->Fill(dattmp->GetBinCenter(ibin),dattmp->GetBinContent(ibin)); 
     }
   }   
   c2_1->SetLogy();
   c2_1->SetGridy();
   c2_1->SetGridx();   
   dat->SetTitle();  
   dat->GetXaxis()->SetTitleOffset(1.5);
   dat->GetXaxis()->SetTitle("Signal efficiency");
   dat->GetYaxis()->SetTitle("Events");
   dat->SetLineColor(2);
   dat->SetLineWidth(2);
   dat->Draw("E");
   eff->SetLineColor(1);
   eff->Draw("LSAME");

// without log scale
   c2_2->cd();
   c2_2->SetLogy(0);
   c2_2->SetGridy();
   c2_2->SetGridx();   
   dat->SetTitle();  
   dat->GetXaxis()->SetTitleOffset(1.5);
   dat->GetXaxis()->SetTitle("Signal efficiency");
   dat->GetYaxis()->SetTitle("Events");
   dat->SetLineColor(2);
   dat->SetLineWidth(2);
   dat->Draw("E");
   eff->SetLineColor(1);
   eff->Draw("LSAME");

}