Attachment 'TSysLimit_searches_onebin.C'
Download 1 #include <iostream>
2 #include <TFile.h>
3 #include <TH1D.h>
4 #include <TRandom3.h>
5 #include <TMath.h>
6 #include <TCanvas.h>
7 #include <TColor.h>
8 #include <TLimit.h>
9 #include <TConfidenceLevel.h>
10 #include <TSysLimitChannel.h>
11 #include <TSysLimit.h>
12 #include <TSysLimitResult.h>
13 #include <TSysLimitScan.h>
14 #include "TTree.h"
15
16 using namespace std;
17 #define Systematics
18
19 int main(int argc, char const *argv[]) {
20
21
22 // TH1::SetDefaultSumw2(kTRUE);
23
24 //#########################################################
25 // Input data
26
27 Double_t lumi = 500.0;
28 Double_t Pe = 0.8;
29 Double_t Pp = -0.3;
30 char filename[256];
31 char _operator[20];
32 Double_t _mass=1;
33 Double_t _spin=0.5;
34 Double_t _lambda =1580;
35 Double_t ExpCL = 0.997; //0.9
36
37
38 if(argc>1) {
39 cout << argv[1] << endl;
40 _mass=atof(argv[1]);
41 } else { exit(1); }
42
43
44 cout<<" Pe = "<<Pe<<"| Pp = "<<Pp<<endl;
45
46 // sprintf(filename,"Results/Sensitivity_Pe%.1f_Pp%.1f_Lumi%.0f_Mass%.0f_3Sigma_CLTSys_onebin_3permile.root",Pe,Pp,lumi,_mass);
47 sprintf(filename,"Results/Sensitivity_Pe%.1f_Pp%.1f_Lumi%.0f_Mass%.0f_3Sigma_CLTSys_onebin_with_errors.root",Pe,Pp,lumi,_mass);
48 // sprintf(filename,"Results/Sensitivity_%s_Pe%.1f_Pp%.1f_Lumi%.0f_Mass%.0f_90CLTSys_lumi_pole_polp_withsys_corel.root",_operator,Pe,Pp,lumi,_mass);
49
50 cout<<filename<<endl;
51 TFile * f= new TFile(filename,"recreate");
52 TList* list = new TList();
53 TH1D* Sensitivity;
54 TTree * t1=new TTree("Parameters","All parameters used in calculation");
55 t1->Branch("Pe", &Pe, "new_v/D");
56 t1->Branch("Pp", &Pp, "new_v/D");
57 t1->Branch("Lumi", &lumi, "new_v/D");
58 t1->Branch("ExpCL", &ExpCL, "new_v/D");
59
60
61 for (int l=0;l<1;l++){
62
63 if (l == 0){
64 sprintf(_operator,"Vector");
65 } else
66 if (l == 1){
67 sprintf(_operator,"Scalar,s-channel");
68 } else
69 if (l == 2){
70 sprintf(_operator,"Axial-vector");
71 } else {cout<<"###########################FALD ############################### "<<endl;}
72
73 sprintf(filename,"Sensitivity %s",_operator);
74 Sensitivity = new TH1D(filename,filename,250,1,250);
75
76
77
78
79
80 cout<<" Mass = "<<_mass<<endl;
81
82
83 sprintf(filename,"/afs/desy.de/group/flc/pool/achaus/WIMPs/data/Templates_Andrii/Lumi%.0f_%s_full/1D_Data_Signal_Pe%.1f_Pp%.1f_Lumi%.0f.root",lumi,_operator,Pe,Pp,lumi);
84 TFile * infile_sig = new TFile(filename,"READ");
85 infile_sig->cd();
86
87
88 sprintf(filename,"M%.0f_S%.1f_%s_Lambda%.0f",_mass,_spin,_operator,_lambda);
89 TH1D* sh=(TH1D*)infile_sig->Get(filename);
90
91
92 sprintf(filename,"/afs/desy.de/group/flc/pool/achaus/WIMPs/data/Templates_Andrii/Data_Background_Pe%.1f_Pp%.1f_Lumi%.0f.root",Pe,Pp,lumi);
93 // sprintf(filename,"/afs/desy.de/group/flc/pool/achaus/WIMPs/data/Templates_Andrii/Background_only_nunugamma/Data_Background_Pe%.1f_Pp%.1f_Lumi%.0f_nunugamma.root",Pe,Pp,lumi);
94 TFile * infile_bg = new TFile(filename,"READ");
95 infile_bg->cd();
96
97 sprintf(filename,"Background_Pe%.1f_Pp%.1f_Lumi%.0f",Pe,Pp,lumi);
98 TH1D* bh=(TH1D*)infile_bg->Get(filename);
99
100
101
102 // cout<<" Integral after 1/lambda^4 :"<<sh->Integral()<<endl;
103 // cout<<" Integral Background :"<<bh->Integral()<<endl;
104
105 sh->Rebin(300);
106 bh->Rebin(300);
107 cout<<" Nbins :"<<bh->GetNbinsX()<<endl;
108
109
110
111 // // consider NSYS systematic sources
112 // // create histograms with relative errors
113 // Int_t NSYS=1;
114 // TString sys_names[NSYS];
115 // sys_names[0]="luminosity";
116 // sys_names[1]="polp";
117 // sys_names[2]="pole";
118 // sys_names[3]="beam";
119 // TH1D *histBgrSys[NSYS][2];
120 // TH1D *histSignalSys[NSYS][2];
121 // for(Int_t s=0;s<NSYS;s++) {
122 // for(Int_t ds=0;ds<2;ds++) {
123 // TString name("Sys_error_");
124 // name += sys_names[s];
125 //
126 // TString name1,name2;
127 //
128 // name2 =name+"_background";
129 // name2 +=ds;
130 // cout<<name2<<endl;
131 // histBgrSys[s][ds]=new TH1D(name2,"Background errors",1,0.0,300.0);
132 //
133 // name1 =name+"_signal";
134 // name1 +=ds;
135 // histSignalSys[s][ds]=new TH1D(name1,"Signal errors",1,0.0,300.0);
136 // cout<<name1<<endl;
137 // }
138 // }
139 //
140 // sh->Scale(pow(_lambda,4));
141 //
142 // // Fill errors histograms
143 // cout<<"BINS = "<<histBgrSys[0][0]->GetNbinsX()<<endl;
144 //
145 // //FILL LUMINOSITY Systematic errors
146 // Int_t s=0;
147 //
148 // for(Int_t i=1;i<=histBgrSys[s][0]->GetNbinsX();i++){
149 // Double_t bg_sys=0.003,sig_sys=0.000;
150 // histBgrSys[s][0]->Fill(i,bg_sys);
151 // histSignalSys[s][0]->Fill(i,sig_sys);
152 // }
153 //
154 //
155 // // cout<<" Integral before :"<<sh->Integral()<<endl;
156 // // sh->Scale(pow(_lambda,4));
157 // // cout<<" Integral after lambda^4 :"<<sh->Integral()<<endl;
158 //
159 //
160 // cout<<"Before add channel"<<endl;
161 // TSysLimitChannel *channel=new TSysLimitChannel();
162 // channel->SetData(bh);
163 // channel->SetPrediction(sh,bh,0); // without errors (sh,bh,0)
164 //
165 //
166 // //################# NUMBER OF SYSTEMATIC ERRORS ########################
167 // NSYS=0;
168 // // ADD SYSTEMATIC ERRORS
169 // for(Int_t m=0;m<NSYS;m++) {
170 // TString name("Sys_error_");
171 // if(m==0){name += "luminosity";
172 // cout<< "DONE "<<name<<endl;
173 // channel->AddSysRelative(histSignalSys[m][0],histBgrSys[m][0],name,1); // correlative add 1
174 // }else if (m==3) {
175 // name += "beam spectra";
176 // cout<< "DONE "<<name<<endl;
177 // channel->AddSysUpDown(histSignalSys[m][0],histSignalSys[m][1],histBgrSys[m][0],histBgrSys[m][1],name,1); // correlative add 1
178 // }
179 // else {name += "polarization";
180 // name +=m;
181 // cout<< "DONE "<<name<<endl;
182 // // channel->AddSysRelative(histSignalSys[m],histBgrSys[m],name); // correlative add 1
183 // channel->AddSysUpDown(histSignalSys[m][0],histSignalSys[m][1],histBgrSys[m][0],histBgrSys[m][1],name,1);
184 // }
185 //
186 //
187 //
188 // }
189
190 #ifdef Systematics
191 cout<<" Calculation with systematic errors :"<<endl;
192 // consider NSYS systematic sources
193 // create histograms with relative errors
194 Int_t NSYS=5;
195 TString sys_names[NSYS];
196 sys_names[0]="luminosity";
197 sys_names[1]="polp";
198 sys_names[2]="pole";
199 sys_names[3]="beam";
200 sys_names[4]="select";
201 TH1D *histBgrSys[NSYS][2];
202 TH1D *histSignalSys[NSYS][2];
203 for(Int_t s=0;s<NSYS;s++) {
204 for(Int_t ds=0;ds<2;ds++) {
205 TString name("Sys_error_");
206 name += sys_names[s];
207
208 TString name1,name2;
209
210 name2 =name+"_background";
211 name2 +=ds;
212 cout<<name2<<endl;
213 histBgrSys[s][ds]=new TH1D(name2,"Background errors",1,0.0,300.0); //300 instead of 1
214
215 name1 =name+"_signal";
216 name1 +=ds;
217 histSignalSys[s][ds]=new TH1D(name1,"Signal errors",1,0.0,300.0); //300 instead of 1
218 cout<<name1<<endl;
219 }
220 }
221 #endif
222 sh->Scale(pow(_lambda,4));
223
224 #ifdef Systematics
225 // Fill errors histograms
226 cout<<"BINS = "<<histBgrSys[0][0]->GetNbinsX()<<endl;
227
228 //FILL LUMINOSITY Systematic errors
229 Int_t s=0;
230
231 for(Int_t i=1;i<=histBgrSys[s][0]->GetNbinsX();i++){
232 Double_t bg_sys=0.0011,sig_sys=0.0011;
233 if(i==1){
234 t1->Branch("bg_sys_lumi", &bg_sys);
235 t1->Branch("sig_sys_lumi", &sig_sys);}
236 histBgrSys[s][0]->Fill(i,bg_sys);
237 histSignalSys[s][0]->Fill(i,sig_sys);
238 }
239
240 //FILL Pp Systematic errors
241 s=1;
242
243 for(Int_t i=1;i<=histBgrSys[s][0]->GetNbinsX();i++){
244 Double_t error_size= 0.0025;
245 // Double_t error_size= 0.001;
246 if(i==1){cout<<"POL P eroor ="<<error_size<<endl;
247 t1->Branch("Pp_Systematic_error", &error_size);}
248 Double_t signal_bin = sh->GetBinContent(i);
249 Double_t bg_bin = bh->GetBinContent(i);
250 // cout<<"signal_bin = "<<signal_bin<<"| bg_bin = "<<bg_bin<<" ;"<<endl;
251
252 Double_t bg_sysup=bg_bin*(1.0+error_size);
253 Double_t bg_sysdown=bg_bin*(1.0-error_size);
254 Double_t sig_sysup=signal_bin*(1.0+error_size);
255 Double_t sig_sysdown=signal_bin*(1.0-error_size);
256
257 histBgrSys[s][0]->Fill(i,bg_sysup);
258 histBgrSys[s][1]->Fill(i,bg_sysdown);
259 histSignalSys[s][0]->Fill(i,sig_sysup);
260 histSignalSys[s][1]->Fill(i,sig_sysdown);
261
262 }
263
264
265 //FILL Pe Systematic errors
266 s=2;
267 for(Int_t i=1;i<=histBgrSys[s][0]->GetNbinsX();i++){
268 Double_t error_size= 0.0025;
269 // Double_t error_size= 0.001;
270 if(i==1){cout<<"POL E eroor ="<<error_size<<endl;
271 t1->Branch("Pe_Systematic_error", &error_size);}
272 Double_t signal_bin = sh->GetBinContent(i);
273 Double_t bg_bin = bh->GetBinContent(i);
274 // cout<<"signal_bin = "<<signal_bin<<"| bg_bin = "<<bg_bin<<" ;"<<endl;
275
276 Double_t bg_sysup=bg_bin*(1.0+error_size);
277 Double_t bg_sysdown=bg_bin*(1.0-error_size);
278 Double_t sig_sysup=signal_bin*(1.0+error_size);
279 Double_t sig_sysdown=signal_bin*(1.0-error_size);
280
281 // cout<<"sig_sys = "<<sig_sys<<"| bg_sys = "<<bg_sys<<" ;"<<endl;
282
283 histBgrSys[s][0]->Fill(i,bg_sysup);
284 histBgrSys[s][1]->Fill(i,bg_sysdown);
285 histSignalSys[s][0]->Fill(i,sig_sysup);
286 histSignalSys[s][1]->Fill(i,sig_sysdown);
287
288 }
289
290 //##########################################################################
291 //FILL Beam spectra parameters
292
293 //FILL Beam spectra parameters errors
294 s=3;
295 Int_t q=1;
296 Double_t error_size= -0.03566; //-0.3566
297 Double_t signal_bin = sh->GetBinContent(q);
298 Double_t bg_bin = bh->GetBinContent(q);
299
300
301 Double_t bg_sysup=bg_bin*(1.0+error_size);
302 Double_t bg_sysdown=bg_bin*(1.0-error_size);
303 Double_t sig_sysup=signal_bin*(1.0+error_size);
304 Double_t sig_sysdown=signal_bin*(1.0-error_size);
305
306
307
308
309
310
311 histBgrSys[s][0]->Fill(q,bg_sysup);
312 histBgrSys[s][1]->Fill(q,bg_sysdown);
313 histSignalSys[s][0]->Fill(q,sig_sysup);
314 histSignalSys[s][1]->Fill(q,sig_sysdown);
315
316
317
318
319
320 //FILL Selection Systematic errors
321 s=4;
322
323 for(Int_t i=1;i<=histBgrSys[s][0]->GetNbinsX();i++){
324 Double_t bg_sys=0.0043,sig_sys=0.0043;
325 if(i==1){
326 t1->Branch("Selection_Systematic_errors_bg", &bg_sys);
327 t1->Branch("Selection_Systematic_errors_sig`", &sig_sys);}
328 histBgrSys[s][0]->Fill(i,bg_sys);
329 histSignalSys[s][0]->Fill(i,sig_sys);
330 }
331
332 #endif
333
334
335
336
337
338
339 cout<<" Integral before :"<<sh->Integral()<<endl;
340 // sh->Scale(pow(_lambda,4));
341 cout<<" Integral after lambda^4 :"<<sh->Integral()<<endl;
342
343
344 cout<<"Before add channel"<<endl;
345 TSysLimitChannel *channel=new TSysLimitChannel();
346 channel->SetData(bh);
347 channel->SetPrediction(sh,bh,0); // without errors (sh,bh,0)
348
349 #ifdef Systematics
350 //################# NUMBER OF SYSTEMATIC ERRORS ########################
351 NSYS=5;
352 // ADD SYSTEMATIC ERRORS
353 for(Int_t m=0;m<NSYS;m++) {
354 TString name("Sys_error_");
355 if(m==0){name += "luminosity";
356 cout<< "DONE "<<name<<endl;
357 channel->AddSysRelative(histSignalSys[m][0],histBgrSys[m][0],name,1); // correlative add 1
358 }else if (m==3) {
359 name += "beam spectra";
360 cout<< "DONE "<<name<<endl;
361 channel->AddSysUpDown(histSignalSys[m][0],histSignalSys[m][1],histBgrSys[m][0],histBgrSys[m][1],name,1); // correlative add 1
362 } else if (m==4) {name += "selection";
363 name +=m;
364 cout<< "DONE "<<name<<endl;
365 channel->AddSysRelative(histSignalSys[m][0],histBgrSys[m][0],name,1); // correlative add 1
366 // channel->AddSysUpDown(histSignalSys[m][0],histSignalSys[m][1],histBgrSys[m][0],histBgrSys[m][1],name,1);
367 }else {name += "polarization";
368 name +=m;
369 cout<< "DONE "<<name<<endl;
370 channel->AddSysUpDown(histSignalSys[m][0],histSignalSys[m][1],histBgrSys[m][0],histBgrSys[m][1],name,1);
371 }
372
373
374
375 }
376
377 #endif
378
379
380 TSysLimit limit(0);
381 limit.AddChannel(channel);
382
383
384 Double_t CL=0.0;
385 Double_t k = 100.0;
386 Double_t k_last = 100.0;
387 Double_t stepwidth = 512;
388 Double_t factor =0;
389 int counter =0;
390
391 while(stepwidth>0.0044)
392
393 {
394
395 CL=0.0;
396 while( CL < (1-ExpCL))
397 {
398
399
400 factor =1/pow(k,4);
401
402
403
404 limit.CalculateCL(factor,50000);
405 TSysLimitResult const *result=limit.GetResult(3); // Log(S+B/B)
406 result->Print();
407 CL=result->GetDataCLs();
408 cout<<"CLs expected = "<<result->GetExpectedCLs()<<", CLs observed "<<CL<<"| k ="<<k<<endl;
409
410
411
412 /* for(Int_t method=0;method<limit.GetNMethod();method++) {
413 TSysLimitResult const *result=limit.GetResult(method);
414 //TSysLimitResult const *result=limit.GetResult();
415 Double_t_t cl;
416 cout<<"=========================================\n";
417 cout<<"Method: "<<method;
418 if(method==limit.GetIndexPBock()) cout<<" (the default method)";
419 cout<<"\n";
420 result->Print();
421
422 cl=result->GetDataCLs();
423 cout<<"CLs = "<<cl;
424 cl=result->GetExpectedCLs();
425 cout<<" expected: "<<cl<<"\n";
426
427 cl=result->GetDataCLb();
428 cout<<"CLb = "<<cl;
429 cl=result->GetExpectedCLb();
430 cout<<" expected: "<<cl<<"\n";
431
432 cl=result->GetDataCLsb();
433 cout<<"CLsb = "<<cl;
434 cl=result->GetExpectedCLsb();
435 cout<<" expected: "<<cl<<"\n";
436
437 } */
438
439
440 if (CL < (1-ExpCL))
441 k_last = k;
442
443 k += stepwidth;
444
445 counter++;
446 cout<<" COUNTER = "<<counter<<endl;
447 cout<<" stepwidth = "<<stepwidth<<endl;
448 if (counter>40)break;
449
450 }
451
452 k = k_last;
453 stepwidth *= 0.5;
454 }
455 cout<<" COUNTER AFTER = "<<counter<<endl;
456
457 Sensitivity->Fill(_mass,k);
458 list->Add(Sensitivity);
459 cout <<"Lambda out "<<k<<endl;
460 cout <<"CLs = "<<CL<<endl;
461
462
463 channel->Delete();
464 // histSignalSys[0]->Delete();
465 // histBgrSys[0]->Delete();
466
467
468
469
470
471 sh->Delete();
472 bh->Delete();
473
474
475 // for(Int_t method=0;method<limit.GetNMethod();method++) {
476 // TSysLimitResult const *result=limit.GetResult(method);
477 // //TSysLimitResult const *result=limit.GetResult();
478 // Double_t_t cl;
479 // cout<<"=========================================\n";
480 // cout<<"Method: "<<method;
481 // if(method==limit.GetIndexPBock()) cout<<" (the default method)";
482 // cout<<"\n";
483 // result->Print();
484 //
485 // cl=result->GetDataCLs();
486 // cout<<"CLs = "<<cl;
487 // cl=result->GetExpectedCLs();
488 // cout<<" expected: "<<cl<<"\n";
489 //
490 // cl=result->GetDataCLb();
491 // cout<<"CLb = "<<cl;
492 // cl=result->GetExpectedCLb();
493 // cout<<" expected: "<<cl<<"\n";
494 //
495 // cl=result->GetDataCLsb();
496 // cout<<"CLsb = "<<cl;
497 // cl=result->GetExpectedCLsb();
498 // cout<<" expected: "<<cl<<"\n";
499 //
500 // /* Double_t_t x,dx;
501 // x=result->GetAvgXData(&dx);
502 // cout<<"X(data) = "<<x<<" +/- "<<dx<<"\n";
503 // x=result->GetAvgXB(&dx);
504 // cout<<"X(B) = "<<x<<" +/- "<<dx<<"\n";
505 // x=result->GetAvgXSB(&dx);
506 // cout<<"X(SB) = "<<x<<" +/- "<<dx<<"\n"; */
507 //
508 // }
509
510 /* cout<<"=========================================\n";
511 cout<<"run TLIMIT without systematic error\n";
512 TConfidenceLevel *tcl=
513 TLimit::ComputeLimit(sh,bh,bh,50000);
514 cout<<" CLb="<<tcl->CLb()<<"\n";
515 cout<<" CLs="<<tcl->CLs()<<"\n";
516 cout<<" CLsb="<<tcl->CLsb()<<"\n";
517 cout<<" Expected CLs="<<tcl->GetExpectedCLs_b()<<" "<<"\n";
518 cout<<" Expected CLsb="<<tcl->GetExpectedCLsb_b()<<" "<<"\n";
519 cout<<" Expected CLb="<<tcl->GetExpectedCLb_b()<<" "<<"\n";
520
521 */
522
523 infile_sig->Close();
524 infile_bg->Close();
525
526 }
527
528 f->cd();
529 list->Print();
530 list->Write();
531 // beamspec->Write();
532 cout<<"END LOOP"<<endl;
533 f->Close();
534 cout<<"SAVE FILE"<<endl;
535
536
537
538
539
540 return 0;
541 }
Attached Files
To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.You are not allowed to attach a file to this page.