ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/AnalysisFramework/various_assignments/ExclusionLineExperimentation/newSMS_utils.C
Revision: 1.1
Committed: Thu Dec 15 15:27:31 2011 UTC (13 years, 4 months ago) by buchmann
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Log Message: 
added utilities

File Contents

# User Rev Content
1 buchmann 1.1 #include "TFile.h"
2     #include "TLine.h"
3     #include "TGraph.h"
4     #include "TH1D.h"
5     #include "TH2D.h"
6     #include "TH3D.h"
7     #include "TText.h"
8     #include <TF1.h>
9     //#include "TLatex.h"
10     //#include "THStack.h"
11     //#include "TDirectory.h"
12     #include "TLegend.h"
13     #include "TCanvas.h"
14     //#include "TMarker.h"
15     #include "TMath.h"
16     #include <vector>
17     #include <iostream>
18     #include <fstream>
19     #include <sstream>
20    
21     using namespace std;
22    
23     bool fail(double xs, double xsLimit) {return xsLimit>1 or !xsLimit;}
24    
25     TH2F *getrelativeXS(TH2F* limit, string myType, double refMult) {
26     //written by Mariarosaria D'Alfonso (dalfonso@mail.cern.ch)
27     //
28     //returns relative cross section (i.e. upper limit / reference cross section)
29     // output can then be used to obtain the exclusion line (get_exclusion_line)
30     string mass="";
31    
32     if((myType=="T1") || (myType=="T1bbbb") || (myType=="T1lnu") || (myType=="T1lh") || (myType=="T5zz") || (myType=="GMSB")) mass="gluino";
33     if(myType=="T2") mass="squark";
34     if(myType=="T2bb") mass="sbottom";
35     if(myType=="T2tt") mass="stop";
36     // if(myType=="T2tt") mass="gluino";
37     if(myType=="T1tttt") mass="gluino";
38    
39     string fileName="reference_xSec.root";
40     if(myType=="T2tt") fileName="/afs/cern.ch/user/d/dalfonso/public/reference_xSec_stop.root";
41    
42     TFile *_file = TFile::Open((char *) fileName.c_str());
43     TH1 *hRef= (TH2F*) _file->Get((char *) mass.c_str());
44    
45     TH2F* limit_ref = (TH2F*) limit->Clone();
46    
47     Int_t nBinX= limit->GetXaxis()->GetNbins();
48     double xLow=limit->GetXaxis()->GetBinLowEdge(1);
49     double xHigh=limit->GetXaxis()->GetBinLowEdge(nBinX+1);
50    
51     Int_t nBinY= limit->GetYaxis()->GetNbins();
52     double yLow=limit->GetYaxis()->GetBinLowEdge(1);
53     double yHigh=limit->GetYaxis()->GetBinLowEdge(nBinY+1);
54    
55     TH2F * xSec_ref= new TH2F("xSec_ref","xSec_ref",nBinX, xLow ,xHigh , nBinY , yLow , yHigh);
56    
57     for(int i=1; i<(nBinX+1); i++) {
58     for(int j=1; j<(nBinY+1); j++) {
59     if((!(limit->GetBinContent(i,j)==0)) && ((myType=="T1") || (myType=="T1bbbb") || (myType=="T1lnu") || (myType=="T1lh") || (myType=="T5zz") || (myType=="GMSB"))) xSec_ref->Fill(limit->GetBinCenter(i),limit->GetBinCenter(j),refMult*hRef->GetBinContent(hRef->FindBin(limit->GetXaxis()->GetBinCenter(i))));
60     if((!(limit->GetBinContent(i,j)==0)) && ((myType=="T2"))) xSec_ref->Fill(limit->GetBinCenter(i),limit->GetBinCenter(j),refMult*(4./5)*hRef->GetBinContent(hRef->FindBin(limit->GetXaxis()->GetBinCenter(i))));
61     if((!(limit->GetBinContent(i,j)==0)) && ((myType=="T2bb"))) xSec_ref->Fill(limit->GetBinCenter(i),limit->GetBinCenter(j),refMult*(1./5)*hRef->GetBinContent(hRef->FindBin(limit->GetXaxis()->GetBinCenter(i))));
62     if((!(limit->GetBinContent(i,j)==0)) && (myType=="T2tt")) xSec_ref->Fill(limit->GetBinCenter(i),limit->GetBinCenter(j),refMult*hRef->GetBinContent(hRef->FindBin(limit->GetXaxis()->GetBinCenter(i))));
63     }
64     }
65    
66     limit_ref->Divide(xSec_ref);
67     xSec_ref->SetDirectory(0);
68     xSec_ref->Delete();
69     hRef->Delete();
70     return limit_ref;
71     }
72    
73     //----------------------------------------------------------------------------------------
74    
75     bool dosmoothing=false; //smooth (neighbors)
76     bool dosupersmoothing=false; //smooth (NN)
77    
78     void project_to_last_point_on_line(float x2, float y2, float x1, float y1, float &x, float &y, string type) {
79     //written by Marco - Andrea Buchmann, ETH Zürich, Switzerland (marco.andrea.buchmann@cern.ch)
80     //currently implemented for GMSB and T5zz
81     float b = 0;
82     if(type=="T5zz") b = -75;
83     x = floor((x2+y2-b)/2.0);
84     y = floor(x + b);
85     if(type!="GMSB"&&type!="T5zz") {
86     //rest not implemented - just return the last value.
87     x=x1;
88     y=y1;
89     }
90     }
91    
92     pair <float,float> find_point(TH2F* histo, int a) {
93     //written by Marco - Andrea Buchmann, ETH Zürich, Switzerland (marco.andrea.buchmann@cern.ch)
94     int nbins=histo->GetNbinsX();
95     bool foundpoint=false;
96     pair <float,float> newpoint;
97     for(int i=nbins;i>0&&!foundpoint;i--) {
98     int b = i + a;
99     if(b>nbins) continue;
100     float value=(histo->GetBinContent(b,i));
101     if(value>0 && value<1) {
102     newpoint.first=histo->GetXaxis()->GetBinLowEdge(b)+histo->GetXaxis()->GetBinWidth(b);
103     newpoint.second=histo->GetYaxis()->GetBinLowEdge(i)+histo->GetYaxis()->GetBinWidth(i);
104     foundpoint=true;
105     }
106     }
107     return newpoint;
108     }
109    
110     bool saveeachexclusionline=true;
111     int plotcounter=0;
112    
113     pair<float,float> pointzero(TH2F *histo) {
114     //written by Marco - Andrea Buchmann, ETH Zürich, Switzerland (marco.andrea.buchmann@cern.ch)
115     bool done=false;
116     pair<float,float> point;
117     for(int j=1;j<=histo->GetNbinsY()&&!done;j++) {
118     for(int i=histo->GetNbinsX();i>0&&!done;i--) {
119     if(histo->GetBinContent(i,j)>0&&histo->GetBinContent(i,j)<1) {
120     point.first=histo->GetXaxis()->GetBinCenter(i);
121     point.second=histo->GetYaxis()->GetBinCenter(j);
122     done=true;
123     }
124     }
125     }
126     return point;
127     }
128    
129    
130     TH2F* get_excluded_region(TH2F *ehisto) {
131     //written by Marco - Andrea Buchmann, ETH Zürich, Switzerland (marco.andrea.buchmann@cern.ch)
132     TH2F *excludedregion = (TH2F*)ehisto->Clone("excludedregion");
133     for(int i=1;i<excludedregion->GetNbinsX();i++) {
134     for(int j=1;j<excludedregion->GetNbinsY();j++) {
135     float val = excludedregion->GetBinContent(i,j);
136     if(val<1&&val>0) excludedregion->SetBinContent(i,j,1);
137     else excludedregion->SetBinContent(i,j,0.0);
138     }
139     }
140     return excludedregion;
141     }
142    
143     TGraph* get_exclusion_line(TH2F *exclusionhisto, string type) {
144     //written by Marco - Andrea Buchmann, ETH Zürich, Switzerland (marco.andrea.buchmann@cern.ch)
145     vector<pair <float,float> > points;
146     points.push_back(pointzero(exclusionhisto));
147     for(int i=exclusionhisto->GetNbinsX();i>0;i--) {
148     pair<float,float> anything = find_point(exclusionhisto,i);
149     if(anything.second>0&&anything.second>0) points.push_back(anything);
150     }
151     Double_t xpoints[points.size()];
152     Double_t ypoints[points.size()];
153    
154     TGraph *graph = new TGraph(points.size()+1);
155     graph->SetLineColor(kBlack);
156     graph->SetLineWidth(2);
157    
158     int pointcounter=0;
159     float lastx=0;
160     float lasty=0;
161     float lastx2=0;
162     float lasty2=0;
163    
164     for(int i=0;i<points.size();i++) {
165     xpoints[i]=points[i].first;
166     ypoints[i]=points[i].second;
167     if(dosmoothing) {
168     if(i>1&&i<points.size()-1) ypoints[i]=(1.0/3.0)*(points[i-1].second+points[i].second+points[i+1].second);
169     if(i>2&&i<points.size()-2) ypoints[i]=(1.0/5.0)*(points[i-2].second+points[i-1].second+points[i].second+points[i+1].second+points[i+2].second);
170     if(i>3&&i<points.size()-3&&dosupersmoothing) ypoints[i]=(1.0/7.0)*(points[i-3].second+points[i-2].second+points[i-1].second+points[i].second+points[i+1].second+points[i+2].second+points[i+3].second);
171    
172     if(i>1&&i<points.size()-1) xpoints[i]=(1.0/3.0)*(points[i-1].first+points[i].first+points[i+1].first);
173     if(i>2&&i<points.size()-2) xpoints[i]=(1.0/5.0)*(points[i-2].first+points[i-1].first+points[i].first+points[i+1].first+points[i+2].first);
174     if(i>3&&i<points.size()-3&&dosupersmoothing) xpoints[i]=(1.0/7.0)*(points[i-3].first+points[i-2].first+points[i-1].first+points[i].first+points[i+1].first+points[i+2].first+points[i+3].first);
175     }
176    
177     graph->SetPoint(pointcounter,xpoints[i],ypoints[i]);
178     lastx2=lastx;
179     lasty2=lasty;
180     lastx=xpoints[i];
181     lasty=ypoints[i];
182     pointcounter++;
183     }
184    
185     for(int i=pointcounter;i<points.size();i++) graph->SetPoint(i,lastx,lasty);
186     if(type=="GMSB"||type=="T5zz") {
187     //The final point will be a continuation of the last one, towards the diagonal
188     //at the moment this is only implemented for GMSB and T5zz. Feel free to add more! (see project_to_last_point_on_line)
189     float x,y;
190     project_to_last_point_on_line(lastx,lasty,lastx2,lasty2,x,y,type);
191     if(x>0&&y>0) graph->SetPoint(points.size(),x,y);
192     else graph->SetPoint(points.size(),lastx,lasty);
193     }
194    
195     if(saveeachexclusionline) {
196     //if you'd like to see how the exclusion line compares to the shape you're excluding, set saveeachexclusionline to true
197     TCanvas *can = new TCanvas("cannie","cannie");
198     exclusionhisto->GetZaxis()->SetRangeUser(0,20);
199     exclusionhisto->SetStats(0);
200     TH2F *exlregion = get_excluded_region(exclusionhisto);
201     exclusionhisto->Draw("COLZ");
202     exlregion->Draw("COLZ");
203     graph->Draw("L*");
204     stringstream graphsaveas;
205     graphsaveas << "new/ExclusionLine___" << plotcounter++ << ".png";
206     can->SaveAs(graphsaveas.str().c_str());
207     delete can;
208     }
209     return graph;
210     }
211    
212     TGraph *get_Marias_exclusion_line(TH2F *histo, string name, float factor) {
213     //written by Marco - Andrea Buchmann, ETH Zürich, Switzerland (marco.andrea.buchmann@cern.ch)
214     return get_exclusion_line(getrelativeXS(histo,name,factor),name);
215     }
216    
217    
218     void draw_Marias_line(string name, string type) {
219     //written by Marco - Andrea Buchmann, ETH Zürich, Switzerland (marco.andrea.buchmann@cern.ch)
220     TFile *f = new TFile("all_bestlimits_PAPER.root");
221     TH2F *hist = (TH2F*)f->Get(name.c_str());
222     if(!hist) {
223     cout << "Whoops no histo found called " << name << endl;
224     f->ls();
225     return;
226     }
227     TGraph *line = get_Marias_exclusion_line(hist,type,1.0);
228     TGraph *lined3 = get_Marias_exclusion_line(hist,type,1.0/3);
229     TGraph *linet3 = get_Marias_exclusion_line(hist,type,3.0);
230     linet3->SetLineColor(kBlue);
231     lined3->SetLineColor(kRed);
232     TCanvas *can = new TCanvas("can","can");
233     can->SetLogz(1);
234     hist->Draw("COLZ");
235     line->Draw();
236     linet3->Draw();
237     lined3->Draw();
238     stringstream saveas;
239     saveas << "new/Exclusionline_NEW_Maria__" << name << "S.png";
240     can->SaveAs(saveas.str().c_str());
241     delete can;
242     }
243    
244     int main() {
245     //written by Marco - Andrea Buchmann, ETH Zürich, Switzerland (marco.andrea.buchmann@cern.ch)
246     cout << "hello" << endl;
247     draw_Marias_line("T5zz","T5zz");
248     draw_Marias_line("T5zzl","T5zz");
249     draw_Marias_line("GMSB","GMSB");
250     return 0;
251     }