MultivariateAnalysis/root/plot_diff.C

{ gROOT->SetStyle("Plain");
  gStyle->SetOptStat(00000); 
  gStyle->SetPalette(1);
  const int nvar =14 ;
  int color;
  TString var[nvar]={"Mt","ht","M123inv","centrality","sphericity","htPlusLepton","getHt3","ht2p","ktMinPrime","DphiJMET","aplanarity", "dPhiLMet","W_MT","metHtPlusLepton"};
  TString lepton_ID = "TM2DCompatibilityLoose_caloIso";
  TString  tmp = "TM2DCompatibilityLoose_caloIso.root";
  TFile * f1 = new TFile(tmp);

  TH1F *hnomsig[nvar];
  TH1F *hrevsig[nvar];
  TH1F *hdiff[nvar];
  Double_t normnomsig, normrevsig;
  Double_t N1[1500][20],N2[1500][20]; 
  for(int j=0;j<nvar;j++){

    hnomsig[j] = (TH1F*)f1->Get("h"+var[j]+"_nominal");
    hrevsig[j] = (TH1F*)f1->Get("h"+var[j]+"_reversed");
    tmp = var[j]+"_diff";
    hdiff[j] = (TH1F*) hnomsig[j]->Clone(tmp);
    hdiff[j]->Reset();
 
    int nbinssig = hnomsig[j]->GetNbinsX()+1 ;
    std::cout<< "var " << var[j] << "nbins " << nbinssig << std::endl;  
 
 
    for(int i = 1; i < nbinssig; i++){
      // std::cout<< var[j] << std::endl;
      N1[i][j] = hnomsig[j]->GetBinContent(i);
      N2[i][j] = hrevsig[j]->GetBinContent(i);
    }
    normnomsig  = hnomsig[j]->Integral();
    if(normnomsig > 0) 
      normnomsig = normnomsig;
    else
      normnomsig = 1;
          

    normrevsig  = hrevsig[j]->Integral();
    if(normrevsig > 0) 
      normrevsig = normrevsig;
    else
      normrevsig = 1;  
    
    hnomsig[j]->Scale(1./normnomsig);
    hrevsig[j]->Scale(1./normrevsig);
  }
 ///////////////////////  DRAWING /////////////////////////////
 
 
    hnomsig[0]->GetYaxis()->SetTitle("Normalized");
    TCanvas *c1[nvar];
    TCanvas *c2[nvar];
    Double_t SigNom, SigRev, Diff, Err;
    Double_t maxnom,maxrev,range;
    Int_t maxsig_bin,maxrev_bin;
for(j=0;j<nvar;j++){
    leg = new TLegend(0.4,0.6,0.65,0.88);
        c1[j] = new TCanvas(var[j],var[j],900,600);
        tmp = var[j]+"_diff";
        c2[j] = new TCanvas(tmp,tmp,900,600);
        c1[j]->cd();
        hnomsig[j]->GetXaxis()->SetTitle(var[j]);
        hnomsig[j]->SetTitle(var[j]);
        hnomsig[j]->SetLineColor(1);
        hrevsig[j]->SetLineColor(2);
        maxsig_bin = hnomsig[j]->GetMaximumBin();
        maxnom =  hnomsig[j]->GetBinContent(maxsig_bin);
        maxrev_bin = hrevsig[j]->GetMaximumBin();
        maxrev =  hrevsig[j]->GetBinContent(maxrev_bin);
        if(maxnom>maxrev)
          range  = maxnom+.1*maxnom;
        else
          range = maxrev+.1*maxrev;
        hnomsig[j]->SetAxisRange(0,range,"Y");
        hnomsig[j]->Draw("HIST");    
        hrevsig[j]->Draw("HIST same");
        tmp = lepton_ID;
        leg->AddEntry(hnomsig[j],tmp,"L");
        tmp = lepton_ID+" Reversed";
        leg->AddEntry(hrevsig[j],tmp,"L");
        leg->SetFillColor(0);
        leg->SetLineColor(0);
        leg->Draw();
        tmp= var[j]+"_"+lepton_ID+".png";
        // c1[j]->Print(tmp);
        int nbins = hnomsig[j]->GetNbinsX();
        double chi2 = 0;
        double chi = 0;
        int ndof=0;
        for(int i = 1; i < nbins+1; i++){
                SigNom = hnomsig[j]->GetBinContent(i);
                SigRev = hrevsig[j]->GetBinContent(i);
                Diff = SigRev - SigNom;
                if( N1[i][j] >0 &&  N2[i][j] > 0)
                  {     
                    Err  =(1./ sqrt(N1[i][j]))*SigNom + (1./ sqrt(N2[i][j]))*SigRev;
                  }
                else
                  {
                    //  std::cout << "diff = " << Diff <<std::endl;
                    Err = 0;}

                if(SigNom > 0)
                  {
                    chi = (Diff*Diff)/SigNom;
                    chi2 = chi2 + chi;
                    ndof++;
                  }
        
                hdiff[j]->SetBinContent(i,Diff);
                hdiff[j]->SetBinError(i,Err);
        }
        std::cout << "Variable = " << var[j] << " Chi2 = " << chi2<< " ndof = " << ndof-1 <<std::endl;
                
        c2[j]->cd();
        hdiff[j]->SetTitle(var[j]+" difference");
        hdiff[j]->SetXTitle(var[j]);
        hdiff[j]->SetYTitle("difference");
        hdiff[j]->Draw();
        tmp= var[j]+"_"+lepton_ID+"_diff.png";
        
        //c2[j]->Print(tmp);


        }
}

 
Revision:	1.4
Committed:	Fri Feb 13 00:53:03 2009 UTC (16 years, 2 months ago) by kukartse
Content type:	text/plain
Branch:	MAIN
CVS Tags:	V00-03-01, ZMorph_BASE_20100408, gak040610_morphing, V00-02-02, gak011410, gak010310, ejterm2010_25nov2009, V00-02-01, V00-02-00, gak112409, CMSSW_22X_branch_base, segala101609, V00-01-15, V00-01-14, V00-01-13, V00-01-12, V00-01-11, V00-01-10, gak031009, gak030509, gak022309, gak021209, HEAD
Branch point for:	ZMorph-V00-03-01, CMSSW_22X_branch
Changes since 1.3:	+25 -8 lines
Log Message:	* empty log message *
#	User	Rev	Content
1	jindal	1.1	{ gROOT->SetStyle("Plain");
2	jindal	1.2	gStyle->SetOptStat(00000);
3	jindal	1.1	gStyle->SetPalette(1);
4	jindal	1.2	const int nvar =14 ;
5	jindal	1.1	int color;
6	jindal	1.2	TString var[nvar]={"Mt","ht","M123inv","centrality","sphericity","htPlusLepton","getHt3","ht2p","ktMinPrime","DphiJMET","aplanarity", "dPhiLMet","W_MT","metHtPlusLepton"};
7	kukartse	1.4	TString lepton_ID = "TM2DCompatibilityLoose_caloIso";
8			TString tmp = "TM2DCompatibilityLoose_caloIso.root";
9	jindal	1.1	TFile * f1 = new TFile(tmp);
10	jindal	1.2
11			TH1F *hnomsig[nvar];
12			TH1F *hrevsig[nvar];
13	jindal	1.1	TH1F *hdiff[nvar];
14	jindal	1.2	Double_t normnomsig, normrevsig;
15			Double_t N1[1500][20],N2[1500][20];
16	jindal	1.1	for(int j=0;j<nvar;j++){
17
18	jindal	1.2	hnomsig[j] = (TH1F*)f1->Get("h"+var[j]+"_nominal");
19			hrevsig[j] = (TH1F*)f1->Get("h"+var[j]+"_reversed");
20			tmp = var[j]+"_diff";
21			hdiff[j] = (TH1F*) hnomsig[j]->Clone(tmp);
22			hdiff[j]->Reset();
23
24			int nbinssig = hnomsig[j]->GetNbinsX()+1 ;
25			std::cout<< "var " << var[j] << "nbins " << nbinssig << std::endl;
26
27
28			for(int i = 1; i < nbinssig; i++){
29			// std::cout<< var[j] << std::endl;
30			N1[i][j] = hnomsig[j]->GetBinContent(i);
31			N2[i][j] = hrevsig[j]->GetBinContent(i);
32			}
33			normnomsig = hnomsig[j]->Integral();
34			if(normnomsig > 0)
35			normnomsig = normnomsig;
36			else
37			normnomsig = 1;
38	jindal	1.1
39
40	jindal	1.2	normrevsig = hrevsig[j]->Integral();
41			if(normrevsig > 0)
42			normrevsig = normrevsig;
43			else
44			normrevsig = 1;
45
46			hnomsig[j]->Scale(1./normnomsig);
47			hrevsig[j]->Scale(1./normrevsig);
48	jindal	1.1	}
49			/////////////////////// DRAWING /////////////////////////////
50
51
52
53	jindal	1.2	hnomsig[0]->GetYaxis()->SetTitle("Normalized");
54			TCanvas *c1[nvar];
55			TCanvas *c2[nvar];
56			Double_t SigNom, SigRev, Diff, Err;
57			Double_t maxnom,maxrev,range;
58			Int_t maxsig_bin,maxrev_bin;
59	jindal	1.1	for(j=0;j<nvar;j++){
60			leg = new TLegend(0.4,0.6,0.65,0.88);
61			c1[j] = new TCanvas(var[j],var[j],900,600);
62			tmp = var[j]+"_diff";
63			c2[j] = new TCanvas(tmp,tmp,900,600);
64			c1[j]->cd();
65	jindal	1.2	hnomsig[j]->GetXaxis()->SetTitle(var[j]);
66			hnomsig[j]->SetTitle(var[j]);
67			hnomsig[j]->SetLineColor(1);
68			hrevsig[j]->SetLineColor(2);
69			maxsig_bin = hnomsig[j]->GetMaximumBin();
70			maxnom = hnomsig[j]->GetBinContent(maxsig_bin);
71			maxrev_bin = hrevsig[j]->GetMaximumBin();
72			maxrev = hrevsig[j]->GetBinContent(maxrev_bin);
73			if(maxnom>maxrev)
74			range = maxnom+.1*maxnom;
75			else
76			range = maxrev+.1*maxrev;
77			hnomsig[j]->SetAxisRange(0,range,"Y");
78			hnomsig[j]->Draw("HIST");
79			hrevsig[j]->Draw("HIST same");
80			tmp = lepton_ID;
81			leg->AddEntry(hnomsig[j],tmp,"L");
82			tmp = lepton_ID+" Reversed";
83			leg->AddEntry(hrevsig[j],tmp,"L");
84	jindal	1.1	leg->SetFillColor(0);
85			leg->SetLineColor(0);
86			leg->Draw();
87	jindal	1.2	tmp= var[j]+"_"+lepton_ID+".png";
88	kukartse	1.4	// c1[j]->Print(tmp);
89	jindal	1.2	int nbins = hnomsig[j]->GetNbinsX();
90	kukartse	1.4	double chi2 = 0;
91			double chi = 0;
92			int ndof=0;
93	jindal	1.1	for(int i = 1; i < nbins+1; i++){
94	jindal	1.2	SigNom = hnomsig[j]->GetBinContent(i);
95			SigRev = hrevsig[j]->GetBinContent(i);
96	kukartse	1.4	Diff = SigRev - SigNom;
97	jindal	1.1	if( N1[i][j] >0 && N2[i][j] > 0)
98	kukartse	1.4	{
99			Err =(1./ sqrt(N1[i][j]))SigNom + (1./ sqrt(N2[i][j]))SigRev;
100			}
101	jindal	1.1	else
102	kukartse	1.4	{
103			// std::cout << "diff = " << Diff <<std::endl;
104			Err = 0;}
105
106			if(SigNom > 0)
107			{
108			chi = (Diff*Diff)/SigNom;
109			chi2 = chi2 + chi;
110			ndof++;
111			}
112
113	jindal	1.1	hdiff[j]->SetBinContent(i,Diff);
114	kukartse	1.4	hdiff[j]->SetBinError(i,Err);
115			}
116			std::cout << "Variable = " << var[j] << " Chi2 = " << chi2<< " ndof = " << ndof-1 <<std::endl;
117	jindal	1.1
118			c2[j]->cd();
119			hdiff[j]->SetTitle(var[j]+" difference");
120			hdiff[j]->SetXTitle(var[j]);
121			hdiff[j]->SetYTitle("difference");
122			hdiff[j]->Draw();
123	jindal	1.3	tmp= var[j]+"_"+lepton_ID+"_diff.png";
124	jindal	1.1
125	kukartse	1.4	//c2[j]->Print(tmp);
126	jindal	1.1
127
128
129
130
131
132
133
134			}
135			}
136
137