ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/AnalysisFramework/Plotting/Modules/SUSYScan.C
Revision: 1.6
Committed: Mon Jul 25 15:57:41 2011 UTC (13 years, 9 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.5: +38 -7 lines
Log Message:
Updated scripts to permit processing in parallel. This affects function brokers but not the functions themselves. When being called the standard way (locally), no changes are in effect

File Contents

# User Rev Content
1 buchmann 1.1 #include <iostream>
2     #include <vector>
3     #include <sys/stat.h>
4     #include <fstream>
5    
6     #include <TCut.h>
7     #include <TROOT.h>
8     #include <TCanvas.h>
9     #include <TMath.h>
10     #include <TColor.h>
11     #include <TPaveText.h>
12     #include <TRandom.h>
13     #include <TH1.h>
14     #include <TH2.h>
15     #include <TF1.h>
16     #include <TSQLResult.h>
17     #include <TProfile.h>
18    
19     //#include "TTbar_stuff.C"
20     using namespace std;
21    
22     using namespace PlottingSetup;
23    
24     void susy_scan_axis_labeling(TH2F *histo) {
25     histo->GetXaxis()->SetTitle("#Chi_{2}^{0}-LSP");
26     histo->GetXaxis()->CenterTitle();
27     histo->GetYaxis()->SetTitle("m_{#tilde{q}}");
28     histo->GetYaxis()->CenterTitle();
29     }
30    
31     void scan_susy_space(string mcjzb, string datajzb) {
32     TCanvas *c3 = new TCanvas("c3","c3");
33     vector<float> binning;
34     binning=allsamples.get_optimal_binsize(mcjzb,cutmass&&cutOSSF&&cutnJets,20,50,800);
35     float arrbinning[binning.size()];
36     for(int i=0;i<binning.size();i++) arrbinning[i]=binning[i];
37     TH1F *puredata = allsamples.Draw("puredata", datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
38     puredata->SetMarkerSize(DataMarkerSize);
39     TH1F *allbgs = allsamples.Draw("allbgs", "-"+datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
40     TH1F *allbgsb = allsamples.Draw("allbgsb", "-"+datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data,luminosity);
41     TH1F *allbgsc = allsamples.Draw("allbgsc", datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data,luminosity);
42     allbgs->Add(allbgsb,-1);
43     allbgs->Add(allbgsc);
44     int ndata=puredata->Integral();
45     ofstream myfile;
46     myfile.open ("susyscan_log.txt");
47     TFile *susyscanfile = new TFile("/scratch/fronga/SMS/T5z_SqSqToQZQZ_38xFall10.root");
48     TTree *suevents = (TTree*)susyscanfile->Get("events");
49     TH2F *exclusionmap = new TH2F("exclusionmap","",20,0,500,20,0,1000);
50     TH2F *exclusionmap1s = new TH2F("exclusionmap1s","",20,0,500,20,0,1000);
51     TH2F *exclusionmap2s = new TH2F("exclusionmap2s","",20,0,500,20,0,1000);
52     TH2F *exclusionmap3s = new TH2F("exclusionmap3s","",20,0,500,20,0,1000);
53    
54     susy_scan_axis_labeling(exclusionmap);
55     susy_scan_axis_labeling(exclusionmap1s);
56     susy_scan_axis_labeling(exclusionmap2s);
57     susy_scan_axis_labeling(exclusionmap3s);
58    
59     Int_t MyPalette[100];
60     Double_t r[] = {0., 0.0, 1.0, 1.0, 1.0};
61     Double_t g[] = {0., 0.0, 0.0, 1.0, 1.0};
62     Double_t b[] = {0., 1.0, 0.0, 0.0, 1.0};
63     Double_t stop[] = {0., .25, .50, .75, 1.0};
64     Int_t FI = TColor::CreateGradientColorTable(5, stop, r, g, b, 100);
65     for (int i=0;i<100;i++) MyPalette[i] = FI+i;
66    
67     gStyle->SetPalette(100, MyPalette);
68    
69     for(int m23=50;m23<500;m23+=25) {
70     for (int m0=(2*(m23-50)+150);m0<=1000;m0+=50)
71     {
72     c3->cd();
73     stringstream drawcondition;
74     drawcondition << "pfJetGoodNum>=3&&(TMath::Abs(masses[0]-"<<m0<<")<10&&TMath::Abs(masses[2]-masses[3]-"<<m23<<")<10)&&mll>5&&id1==id2";
75     TH1F *puresignal = new TH1F("puresignal","puresignal",binning.size()-1,arrbinning);
76     TH1F *puresignall= new TH1F("puresignall","puresignal",binning.size()-1,arrbinning);
77     stringstream drawvar,drawvar2;
78     drawvar<<mcjzb<<">>puresignal";
79     drawvar2<<"-"<<mcjzb<<">>puresignall";
80     suevents->Draw(drawvar.str().c_str(),drawcondition.str().c_str());
81     suevents->Draw(drawvar2.str().c_str(),drawcondition.str().c_str());
82     if(puresignal->Integral()<60) {
83     delete puresignal;
84     continue;
85     }
86     puresignal->Add(puresignall,-1);//we need to correct for the signal contamination - we effectively only see (JZB>0)-(JZB<0) !!
87     puresignal->Scale(ndata/(20*puresignal->Integral()));//normalizing it to 5% of the data
88     stringstream saveas;
89     saveas<<"Model_Scan/m0_"<<m0<<"__m23_"<<m23;
90     dout << "PLEASE KEEP IN MIND THAT SIGNAL CONTAMINATION IS NOT REALLY TAKEN CARE OF YET DUE TO LOW STATISTICS! SHOULD BE SOMETHING LIKE THIS : "<< endl;
91     // TH1F *signalpredlo = allsamples.Draw("signalpredlo", "-"+mcjzb, binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc, luminosity,allsamples.FindSample("LM4"));
92     // TH1F *signalpredro = allsamples.Draw("signalpredro", mcjzb, binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc, luminosity,allsamples.FindSample("LM4"));
93     // TH1F *puredata = allsamples.Draw("puredata", datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
94     // signalpred->Add(signalpredlo,-1);
95     // signalpred->Add(signalpredro);
96     // puresignal->Add(signalpred,-1);//subtracting signal contamination
97     //---------------------
98     // dout << "(m0,m23)=("<<m0<<","<<m23<<") contains " << puresignal->Integral() << endl;
99     // TH1F *puresignal = allsamples.Draw("puresignal",mcjzb, binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc, luminosity,allsamples.FindSample("LM4"));
100     vector<float> results=establish_upper_limits(puredata,allbgs,puresignal,saveas.str(),myfile);
101     if(results.size()==0) {
102     delete puresignal;
103     continue;
104     }
105     exclusionmap->Fill(m23,m0,results[0]);
106     exclusionmap1s->Fill(m23,m0,results[1]);
107     exclusionmap2s->Fill(m23,m0,results[2]);
108     exclusionmap3s->Fill(m23,m0,results[3]);
109     delete puresignal;
110     dout << "(m0,m23)=("<<m0<<","<<m23<<") : 3 sigma at " << results[3] << endl;
111     }
112     }//end of model scan for loop
113    
114     dout << "Exclusion Map contains" << exclusionmap->Integral() << " (integral) and entries: " << exclusionmap->GetEntries() << endl;
115     c3->cd();
116     exclusionmap->Draw("CONTZ");
117     CompleteSave(c3,"Model_Scan/CONT/Model_Scan_Mean_values");
118     exclusionmap->Draw("COLZ");
119     CompleteSave(c3,"Model_Scan/COL/Model_Scan_Mean_values");
120    
121     exclusionmap1s->Draw("CONTZ");
122     CompleteSave(c3,"Model_Scan/CONT/Model_Scan_1sigma_values");
123     exclusionmap1s->Draw("COLZ");
124     CompleteSave(c3,"Model_Scan/COL/Model_Scan_1sigma_values");
125    
126     exclusionmap2s->Draw("CONTZ");
127     CompleteSave(c3,"Model_Scan/CONT/Model_Scan_2sigma_values");
128     exclusionmap2s->Draw("COLZ");
129     CompleteSave(c3,"Model_Scan/COL/Model_Scan_2sigma_values");
130    
131     exclusionmap3s->Draw("CONTZ");
132     CompleteSave(c3,"Model_Scan/CONT/Model_Scan_3sigma_values");
133     exclusionmap3s->Draw("COLZ");
134     CompleteSave(c3,"Model_Scan/COL/Model_Scan_3sigma_values");
135    
136     TFile *exclusion_limits = new TFile("exclusion_limits.root","RECREATE");
137     exclusionmap->Write();
138     exclusionmap1s->Write();
139     exclusionmap2s->Write();
140     exclusionmap3s->Write();
141     exclusion_limits->Close();
142     susyscanfile->Close();
143    
144     myfile.close();
145     }
146    
147 buchmann 1.4 void efficiency_scan_in_susy_space(string mcjzb, string datajzb, bool requireZ, float peakerror) {
148     dout << "Doing efficiency scan using " << (scansample.collection)[0].filename << endl;
149 buchmann 1.5 geqleq="geq";
150 buchmann 1.1 automatized=true;
151 buchmann 1.2 mcjzbexpression=mcjzb;
152 buchmann 1.1
153     string massgluname="MassGlu";
154     string massLSPname="MassLSP";
155    
156     Int_t MyPalette[100];
157     Double_t r[] = {0., 0.0, 1.0, 1.0, 1.0};
158     Double_t g[] = {0., 0.0, 0.0, 1.0, 1.0};
159     Double_t b[] = {0., 1.0, 0.0, 0.0, 1.0};
160     Double_t stop[] = {0., .25, .50, .75, 1.0};
161 buchmann 1.2 Int_t FI = TColor::CreateGradientColorTable(5, stop, r, g, b, 110);
162 buchmann 1.1 for (int i=0;i<100;i++) MyPalette[i] = FI+i;
163    
164     gStyle->SetPalette(100, MyPalette);
165    
166     TH2F *efficiencymap = new TH2F("efficiencymap","",(mgluend-mglustart)/mglustep+1,mglustart-0.5*mglustep,mgluend+0.5*mglustep,(mLSPend-mLSPstart)/mLSPstep+1,mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
167 buchmann 1.4 TH2F *Neventsmap = new TH2F("Neventsmap","", (mgluend-mglustart)/mglustep+1,mglustart-0.5*mglustep,mgluend+0.5*mglustep,(mLSPend-mLSPstart)/mLSPstep+1,mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
168 buchmann 1.2
169 buchmann 1.1 float rightmargin=gStyle->GetPadRightMargin();
170 buchmann 1.2 gStyle->SetPadRightMargin(0.14);
171    
172 buchmann 1.1 TCanvas *effcanvas = new TCanvas("effcanvas","Efficiency canvas");
173 buchmann 1.2
174 buchmann 1.1 for(int mglu=mglustart;mglu<=mgluend;mglu+=mglustep) {
175     for (int mlsp=mLSPstart;mlsp<=mLSPend&&mlsp<=mglu;mlsp+=mLSPstep)
176     {
177     float result,resulterr;
178     stringstream addcut;
179     addcut << "(TMath::Abs("<<massgluname<<"-"<<mglu<<")<5)&&(TMath::Abs("<<massLSPname<<"-"<<mlsp<<")<5)";
180     MCefficiency((scansample.collection)[0].events,result,resulterr,mcjzb,requireZ,addcut.str());
181     if(result!=result||isanyinf(result)) {
182 buchmann 1.4 dout << "Watch out the result for mlsp=" << mlsp << " and mglu=" << mglu << " (" << result << ") has been detected as being not a number (nan) or infinity (inf) !" << endl;
183 buchmann 1.1 result=0; // kicking nan and inf errors
184     }
185 buchmann 1.5 else {
186     efficiencymap->Fill(mglu,mlsp,result);
187     dout << " ok! Added efficiency " << result << " for mlsp=" << mlsp << " and mglu=" << mglu << endl;
188     }
189 buchmann 1.2
190     (scansample.collection)[0].events->Draw(mcjzb.c_str(),addcut.str().c_str(),"goff");
191     float nevents = (scansample.collection)[0].events->GetSelectedRows();
192     Neventsmap->Fill(mglu,mlsp,nevents);
193 buchmann 1.1 }
194     }
195    
196     efficiencymap->GetXaxis()->SetTitle("m_{glu}");
197     efficiencymap->GetXaxis()->CenterTitle();
198     efficiencymap->GetYaxis()->SetTitle("m_{LSP}");
199     efficiencymap->GetYaxis()->CenterTitle();
200    
201 buchmann 1.2 Neventsmap->GetXaxis()->SetTitle("m_{glu}");
202     Neventsmap->GetXaxis()->CenterTitle();
203     Neventsmap->GetYaxis()->SetTitle("m_{LSP}");
204     Neventsmap->GetYaxis()->CenterTitle();
205 buchmann 1.4
206 buchmann 1.1 efficiencymap->Draw("COLZ");
207     TText *title = write_title("Efficiency in LSP-Glu plane");
208     title->Draw();
209     CompleteSave(effcanvas,"SUSYScan/Efficency");
210 buchmann 1.4
211 buchmann 1.2 Neventsmap->Draw("COLZ");
212     TText *title2 = write_title("Number of events in LSP-Glu plane");
213     title2->Draw();
214 buchmann 1.5 CompleteSave(effcanvas,"SUSYScan/Nevents");
215 buchmann 1.4 }
216    
217 buchmann 1.6 void scan_SUSY_parameter_space(string mcjzb,string datajzb,vector<float> jzb_cut,bool requireZ, float peakerror, int ibin,float njobs=-1, float jobnumber=-1) {
218    
219     bool runninglocally=true;
220     if(njobs>-1&&jobnumber>-1) {
221     runninglocally=false;
222     cout << "Running on the GRID (this is job " << jobnumber << "/" << njobs << ")" << endl;
223     } else {
224     cout << "Running locally " << endl;
225     }
226 buchmann 1.4
227     jzbSel=jzb_cut[ibin];
228     geqleq="geq";
229     automatized=true;
230     mcjzbexpression=mcjzb;
231    
232     string massgluname="MassGlu";
233     string massLSPname="MassLSP";
234    
235     Int_t MyPalette[100];
236     Double_t r[] = {0., 0.0, 1.0, 1.0, 1.0};
237     Double_t g[] = {0., 0.0, 0.0, 1.0, 1.0};
238     Double_t b[] = {0., 1.0, 0.0, 0.0, 1.0};
239     Double_t stop[] = {0., .25, .50, .75, 1.0};
240     Int_t FI = TColor::CreateGradientColorTable(5, stop, r, g, b, 110);
241     for (int i=0;i<100;i++) MyPalette[i] = FI+i;
242    
243     gStyle->SetPalette(100, MyPalette);
244    
245 buchmann 1.6 TH2F *limitmap = new TH2F(("limitmap"+any2string(jzbSel)).c_str(),"",(mgluend-mglustart)/mglustep+1,mglustart-0.5*mglustep,mgluend+0.5*mglustep,(mLSPend-mLSPstart)/mLSPstep+1,mLSPstart-0.5*mLSPstep,mLSPend+0.5*mLSPstep);
246 buchmann 1.4
247     float rightmargin=gStyle->GetPadRightMargin();
248     gStyle->SetPadRightMargin(0.15);
249    
250     TCanvas *limcanvas = new TCanvas("limcanvas","Limit canvas");
251    
252 buchmann 1.6
253     int Npoints=0;
254     for(int mglu=mglustart;mglu<=mgluend;mglu+=mglustep) {
255     for (int mlsp=mLSPstart;mlsp<=mLSPend&&mlsp<=mglu;mlsp+=mLSPstep) Npoints++;
256     }
257    
258     int ipoint=-1;
259 buchmann 1.4 for(int mglu=mglustart;mglu<=mgluend;mglu+=mglustep) {
260     for (int mlsp=mLSPstart;mlsp<=mLSPend&&mlsp<=mglu;mlsp+=mLSPstep)
261     {
262 buchmann 1.6 ipoint++;
263     if(!runninglocally&&!do_this_point(ipoint,Npoints,jobnumber,njobs)) continue;
264 buchmann 1.4 float result,resulterr;
265     stringstream addcut;
266     addcut << "(TMath::Abs("<<massgluname<<"-"<<mglu<<")<5)&&(TMath::Abs("<<massLSPname<<"-"<<mlsp<<")<5)";
267     vector<vector<float> > systematics;
268 buchmann 1.6 MCefficiency((scansample.collection)[0].events,result,resulterr,mcjzb,requireZ,addcut.str());
269     if(result!=result||resulterr!=resulterr) {
270     dout << "Limits can't be calculated for this configuration (mglu="<<mglu<<" , mlsp="<<mlsp << ") as either the efficiency or its error are not numbers! (result="<<result<<" and resulterr="<<resulterr<<")"<< endl;
271     continue;
272     } else {
273     dout << "Configuration (mglu="<<mglu<<" , mlsp="<<mlsp << ") passed the mcefficiency test; will establish limits now." << endl;
274     }
275 buchmann 1.4 do_systematics_for_one_file((scansample.collection)[0].events,"SUSY SCAN", systematics,mcjzb,datajzb,peakerror,requireZ, addcut.str());
276     float JZBcutat=systematics[0][0];
277     float mceff=systematics[0][1];
278     float toterr =systematics[0][4];
279     if(mceff!=mceff||toterr!=toterr) {
280 buchmann 1.6 dout << "Limits can't be calculated for this configuration (mglu="<<mglu<<" , mlsp="<<mlsp << ") as either the efficiency or its error are not numbers! (mceff="<<mceff<<" and toterr="<<toterr<<")"<< endl;
281 buchmann 1.4 continue;
282 buchmann 1.5 } else {
283 buchmann 1.6 dout << "Calculating limit now for m_{glu}="<<mglu << " and m_{lsp}="<<mlsp<<endl;
284 buchmann 1.5 vector<float> sigmas = compute_one_upper_limit(mceff,toterr,ibin,mcjzb);
285     limitmap->Fill(mglu,mlsp,sigmas[0]);
286     dout << "A limit has been added at " << sigmas[0] << " for m_{glu}="<<mglu << " and m_{lsp}="<<mlsp<<endl;
287 buchmann 1.4 }
288     }
289     }
290    
291     limitmap->GetXaxis()->SetTitle("m_{glu}");
292     limitmap->GetXaxis()->CenterTitle();
293     limitmap->GetYaxis()->SetTitle("m_{LSP}");
294     limitmap->GetYaxis()->CenterTitle();
295    
296     limitmap->GetXaxis()->SetTitle("m_{glu}");
297     limitmap->GetXaxis()->CenterTitle();
298     limitmap->GetYaxis()->SetTitle("m_{LSP}");
299     limitmap->GetYaxis()->CenterTitle();
300    
301     limcanvas->cd();
302     limitmap->Draw("COLZ");
303     TText *title = write_title("Limits in LSP-Glu plane");
304     title->Draw();
305 buchmann 1.6 if(runninglocally) {
306     CompleteSave(limcanvas,"SUSYScan/Limits_JZB_geq"+any2string(jzb_cut[ibin]));
307     } else {
308     TFile *outputfile=new TFile(("output/DistributedLimits_job"+string(any2string(jobnumber))+"_of_"+string(any2string(njobs))+".root").c_str(),"UPDATE");
309     limitmap->Write();
310     outputfile->Close();
311     }
312 buchmann 1.1 }
313 buchmann 1.4
314 buchmann 1.6 void scan_SUSY_parameter_space(string mcjzb,string datajzb,vector<float> jzb_cut,bool requireZ, float peakerror, float njobs=-1, float jobnumber=-1) {
315 buchmann 1.4 for(int ibin=0;ibin<jzb_cut.size();ibin++) {
316 buchmann 1.6 scan_SUSY_parameter_space(mcjzb,datajzb,jzb_cut,requireZ, peakerror, ibin, njobs, jobnumber);
317 buchmann 1.4 }
318 buchmann 1.6 }