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

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.
  • [get | view] (2012-01-30 13:28:08, 4.3 KB) [[attachment:CalibrateData.C]]
  • [get | view] (2012-01-30 13:27:40, 12.0 KB) [[attachment:CalibrationFunctions.C]]
  • [get | view] (2014-06-20 09:25:03, 5.2 KB) [[attachment:CreateDataBackground.C]]
  • [get | view] (2014-06-20 09:21:28, 22.3 KB) [[attachment:CreateDataForSignalAndBackground.C]]
  • [get | view] (2014-06-20 09:24:20, 0.2 KB) [[attachment:CreateDataForSignalAndBackground.py]]
  • [get | view] (2012-01-30 12:32:43, 12.0 KB) [[attachment:FindMergeAngle.C]]
  • [get | view] (2012-01-30 13:22:33, 9.8 KB) [[attachment:PreprocessData.C]]
  • [get | view] (2012-01-30 13:51:59, 16.8 KB) [[attachment:Selection.C]]
  • [get | view] (2014-06-20 09:26:52, 13.2 KB) [[attachment:TSysLimit_searches.C]]
  • [get | view] (2014-06-20 09:28:25, 14.9 KB) [[attachment:TSysLimit_searches_onebin.C]]
  • [get | view] (2014-06-20 09:25:07, 12.9 KB) [[attachment:WIMP_functions.C]]
  • [get | view] (2014-06-20 12:52:10, 1.3 KB) [[attachment:runLambdaSearch.sh]]
  • [get | view] (2014-06-20 12:52:16, 0.2 KB) [[attachment:submitter.sh]]
 All files | Selected Files: delete move to page copy to page

You are not allowed to attach a file to this page.