source: trunk/MagicSoft/Mars/mtemp/mberlin/macros/CutOptim.C@ 5138

Last change on this file since 5138 was 4137, checked in by hengsteb, 21 years ago
*** empty log message ***
File size: 13.4 KB
Line 
1#include "TMinuit.h"
2#include "TMatrix.h"
3#include "TString.h"
4#include "TCanvas.h"
5#include "TH1F.h"
6#include "TF1.h"
7#include "TLegend.h"
8#include "TFile.h"
9#include "TArrayF.h"
10#include "TGraphErrors.h"
11#include "MRawRunHeader.h"
12#include "TPad.h"
13#include "TPaveText.h"
14#include "iostream.h"
15#include "TObjArray.h"
16
17#include "MHMatrix.h"
18
19//--------------------------------------------------------------------------------------
20// MAIN INPUT/OUTPUT
21
22#include "IOMkn421.h"
23//--------------------------------------------------------------------------------------
24
25Double_t acut = 10; // alpha region of excess
26Bool_t aplot = 1; // make only alpha plot (1) , do sig. optimization and alpha plot (0)
27
28//--------------------------------------------------------------------------------------------------------------
29//--------------------------------------------------------------------------------------------------------------
30
31const TMatrix *matOn;
32const TMatrix *matOff;
33
34Int_t isize = 0;
35Int_t idist = 1;
36Int_t iwidth = 2;
37Int_t ilength = 3;
38Int_t ialpha = 4;
39
40//------------------------------------------------------------------------------
41//------------------------------------------------------------------------------
42
43Double_t signal(Double_t *x, Double_t *par){
44 return par[0]/(TMath::Sqrt(2.*TMath::Pi())*par[1])*TMath::Exp(-(x[0]-par[2])*(x[0]-par[2])/2./par[1]/par[1]);
45}
46
47Double_t background(Double_t *x, Double_t *par){
48 return par[0]+par[1]*x[0]+par[2]*pow(x[0],2.)+par[3]*x[0]*x[0]*x[0];
49 //return par[0]+par[1]*x[0]+par[2]*x[0]*x[0]+par[3]*x[0]*x[0]*x[0];
50}
51
52Double_t fitFunction(Double_t *x, Double_t *par){
53 return signal(x,par) + background(x,&par[3]);
54}
55
56//------------------------------------------------------------------------------
57//------------------------------------------------------------------------------
58void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
59{
60 Double_t w1=par[0];
61 Double_t w2=par[1];
62 Double_t l1=par[2];
63 Double_t l2=par[3];
64 Double_t d1=par[4];
65 Double_t d2=par[5];
66
67 Double_t non=0;
68 Double_t noff=0;
69
70 TH1F hon("hon","",18,0.,90.);//5 deg bins
71 for(Int_t i=0;i<matOn->GetNrows();i++)
72 {
73 if( (*matOn)(i,isize)<SizeLo) continue;
74 if(((*matOn)(i,idist)>d1) ||((*matOn)(i,idist)<d2 )) continue;
75
76 if(((*matOn)(i,iwidth)>w1) ||((*matOn)(i,iwidth)<w2)) continue;
77 if(((*matOn)(i,ilength)>l1)||((*matOn)(i,ilength)<l2)) continue;
78
79 Float_t alpha=(*matOn)(i,ialpha);
80 hon.Fill(TMath::Abs(alpha));
81 }
82
83 for(Int_t i=1;i<=4;i++)
84 non+=hon.GetBinContent(i);
85
86 //-----------------------------------------------------------------
87 /*TH1F hoff("hoff","",18,0.,90.);//5 deg bins
88 for(Int_t i=0;i<matOff->GetNrows();i++)
89 {
90 if( (*matOff)(i,isize)<SizeLo) continue;
91 if(((*matOff)(i,idist)>d1) ||((*matOff)(i,idist)<d2 )) continue;
92
93 if(((*matOff)(i,iwidth)>w1) ||((*matOff)(i,iwidth)<w2)) continue;
94 if(((*matOff)(i,ilength)>l1)||((*matOff)(i,ilength)<l2)) continue;
95
96 Float_t alpha=(*matOff)(i,ialpha);
97 hoff.Fill(TMath::Abs(alpha));
98 }
99
100 for(Int_t i=1;i<=4;i++)
101 noff+=hoff.GetBinContent(i);*/
102 //-----------------------------------------------------------------
103
104 // flat background region
105 for(Int_t i=7;i<=15;i++)
106 noff+=hon.GetBinContent(i);
107
108 Float_t offscale=4./9.;
109 noff*=offscale;
110
111 if(non+noff<0.001)
112 f=0.;
113 else
114 f = -(non-noff)/(TMath::Sqrt(non+noff));
115
116 cout<<"Non="<<non<<" Noff="<<noff*offscale<<" Sig="<<-f<<endl;
117}
118//------------------------------------------------------------------------------
119//------------------------------------------------------------------------------
120void CutOptim()
121{
122 //-----------------------------------------------------------------------
123 // read matrices
124
125 MHMatrix hmatOn;
126 TFile *fileOn=new TFile(fileMatOn.Data());
127 hmatOn.Read("Matrix");
128 delete fileOn;
129 const TMatrix &mOn=hmatOn.GetM();
130 matOn=(const TMatrix*)&mOn;
131
132 MHMatrix hmatOff;
133 TFile *fileOff=new TFile(fileMatOff.Data());
134 hmatOff.Read("Matrix");
135 delete fileOff;
136 const TMatrix &mOff=hmatOff.GetM();
137 matOff=(const TMatrix*)&mOff;
138
139 //-----------------------------------------------------------------------
140
141 Int_t nrowOn = matOn->GetNrows();
142 Int_t nrowOff = matOff->GetNrows();
143
144 Int_t ncolOn = matOn->GetNcols();
145 Int_t ncolOff = matOff->GetNcols();
146
147 cout<<"matrix On: ncol="<<ncolOn<<" nrow="<<nrowOn<<endl;
148 cout<<"matrix Off: ncol="<<ncolOff<<" nrow="<<nrowOff<<endl;
149
150 TCanvas *c = new TCanvas("c1","Alpha distribution after g/h separation",200,10,800,500);
151 c->SetFillColor(0);
152
153 //-----------------------------------------------------------------------
154 // optimize cuts
155
156 Int_t iflag;
157 Double_t val;
158 Double_t vstart[6];
159
160 Double_t *gin=NULL;
161
162 Double_t arglist[10];
163 Int_t ierflg = 0;
164
165 TMinuit *gMinuit = new TMinuit(10);
166 TGraph* gr=NULL;
167
168 if(!aplot)
169 c->Divide(2,1);
170 else
171 goto APlot;
172
173 // start values
174 vstart[0]=WidthUp;
175 vstart[1]=WidthLo;
176 vstart[2]=LengthUp;
177 vstart[3]=LengthLo;
178 vstart[4]=DistUp;
179 vstart[5]=DistLo;
180
181 gMinuit->SetFCN(fcn);
182
183 arglist[0] = 1;
184 gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);
185
186 // Set starting values and step sizes for parameters
187 static Double_t step[6];
188 const Double_t def_step=0.2;
189 step[0]=def_step;
190 step[1]=def_step;
191 step[2]=def_step;
192 step[3]=def_step;
193 step[4]=def_step;
194 step[5]=def_step;
195
196 gMinuit->mnparm(0,"w1", vstart[0], step[0], 0.,0.,ierflg);
197 gMinuit->mnparm(1,"w2", vstart[1], step[1], 0.,0.,ierflg);
198 gMinuit->mnparm(2,"l1", vstart[2], step[2], 0.,0.,ierflg);
199 gMinuit->mnparm(3,"l2", vstart[3], step[3], 0.,0.,ierflg);
200 gMinuit->mnparm(4,"d1", vstart[4], step[4], 0.2,1.,ierflg);
201 gMinuit->mnparm(5,"d2", vstart[5], step[5], 0.5,1.5,ierflg);
202
203 //gMinuit->FixParameter(0); // WidthUp
204 gMinuit->FixParameter(1); // WidthLo
205 //gMinuit->FixParameter(2); // LengthUp
206 gMinuit->FixParameter(3); // LengthLo
207 gMinuit->FixParameter(4); // DistUp
208 gMinuit->FixParameter(5); // DistLo
209
210 gMinuit->Eval(6,gin,val,vstart,iflag);
211
212 arglist[0] = 5000;
213 arglist[1] = 1.;
214
215 //gMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);
216 gMinuit->mnexcm("SIM", arglist ,2,ierflg);
217
218 // Print results
219 Double_t amin,edm,errdef;
220 Int_t nvpar,nparx,icstat;
221 gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
222 gMinuit->mnprin(3,amin);
223
224 Double_t fpar[6];
225 Double_t fparerr[6];
226 for(Int_t i=0;i<6;i++)
227 gMinuit->GetParameter(i,fpar[i],fparerr[i]);
228
229 WidthUp = fpar[0];
230 WidthLo = fpar[1];
231 LengthUp = fpar[2];
232 LengthLo = fpar[3];
233 DistUp = fpar[4];
234 DistLo = fpar[5];
235
236
237 // parameter's sigma plots
238 c->cd(2);
239 // 2 sigma contour
240 /*
241 gMinuit->SetErrorDef(4);
242 gr = (TGraph*) gMinuit->Contour(100,0,2);
243 gr->DrawClone("AL");*/
244
245 // 1 sigma contour
246 gMinuit->SetErrorDef(1);
247 gr = (TGraph*) gMinuit->Contour(100,0,2);
248 gr->DrawClone("AL");
249 /*
250 // contour of min. function
251 c->cd(2);
252 arglist[0] = 1; // WidthUp
253 gMinuit->mnexcm("SCA", arglist ,1,ierflg);
254 TGraph* gr = (TGraph*)gMinuit->GetPlot();
255 gr->DrawClone("al");
256
257 c->cd(4);
258 arglist[0] = 4; // LengthLo
259 gMinuit->mnexcm("SCA", arglist ,1,ierflg);
260 gr = (TGraph*)gMinuit->GetPlot();
261 gr->DrawClone("al");//*/
262
263 c->cd(1);
264
265APlot:
266 //-----------------------------------------------------------------------
267 // plot alpha histogram with current cuts and calc significance
268
269 // alpha histogram
270 TH1F *halpha=new TH1F("halpha","",18,0.,90.);
271 TH1F *halphaOff=new TH1F("halphaOff","",18,0.,90.);
272
273 halpha->SetStats(kFALSE);
274 halphaOff->SetStats(kFALSE);
275
276 for(Int_t i=0;i<nrowOn;i++)
277 {
278 if (mOn(i,isize)<SizeLo) continue;
279 if((mOn(i,idist)<DistLo)||(mOn(i,idist)>DistUp)) continue;
280
281 if((mOn(i,iwidth) < WidthLo) ||(mOn(i,iwidth) > WidthUp)) continue;
282 if((mOn(i,ilength)< LengthLo)||(mOn(i,ilength)> LengthUp)) continue;
283
284 halpha->Fill(TMath::Abs(mOn(i,ialpha)));
285 }
286 halpha->SetLineWidth(2);
287 halpha->DrawClone();
288
289 for(Int_t i=0;i<nrowOff;i++)
290 {
291 if (mOff(i,isize)<SizeLo) continue;
292 if((mOff(i,idist)<DistLo)||(mOff(i,idist)>DistUp)) continue;
293
294 if((mOff(i,iwidth) < WidthLo) ||(mOff(i,iwidth) > WidthUp)) continue;
295 if((mOff(i,ilength)< LengthLo)||(mOff(i,ilength)> LengthUp)) continue;
296
297 halphaOff->Fill(TMath::Abs(mOff(i,ialpha)));
298
299 }
300
301 // scale Off histogram to On data using bins 6-18
302 Double_t nsample_on=0;
303 Double_t nsample_off=0;
304 for(Int_t i=6;i<18;i++)
305 nsample_on+=halpha->GetBinContent(i);
306 for(Int_t i=6;i<18;i++)
307 nsample_off+=halphaOff->GetBinContent(i);
308
309 Double_t scal=nsample_on/nsample_off;
310
311 halphaOff->SetLineColor(2);
312 halphaOff->Scale(scal);
313 halphaOff->Draw("same");
314
315 gPad->SetGridx();
316 gPad->SetGridy();
317
318 // define functions
319 TF1 *backFcn = new TF1("backFcn",background,20.,89.5,4);
320 backFcn->SetLineWidth(2);
321 backFcn->SetLineColor(kRed);
322
323 TF1 *signalFcn = new TF1("signalFcn",signal,0.,89.5,3);
324 signalFcn->SetLineWidth(2);
325 signalFcn->SetLineColor(kBlue);
326 signalFcn->SetNpx(500);
327
328 TF1 *fitFcn = new TF1("fitFcn",fitFunction,0.,89.5,7);
329 fitFcn->SetNpx(500);
330 fitFcn->SetLineWidth(2);
331 fitFcn->SetLineColor(kMagenta);
332
333 // start-parameters
334 Double_t p3start=halpha->GetBinContent(halpha->GetNbinsX()-1);
335 Double_t p4start=0.;
336 Double_t p5start=0.;
337 Double_t p6start=0.;
338
339 Double_t p0start=(halpha->GetEntries()*halpha->GetBinWidth(1)-p3start)*0.05;
340 Double_t p1start=5.;
341 Double_t p2start=0.;
342
343 cout<<"start values:"<<endl;
344
345 cout<<" gaussian:"<<endl;
346 cout<<" p0 = "<<p0start<<endl;
347 cout<<" p1 = "<<p1start<<endl;
348 cout<<" p2 = "<<p2start<<endl;
349
350 cout<<" pol3:"<<endl;
351 cout<<" p0 = "<<p3start<<endl;
352 cout<<" p1 = "<<p4start<<endl;
353 cout<<" p2 = "<<p5start<<endl;
354 cout<<" p3 = "<<p6start<<endl<<endl;
355
356
357 // background fit
358 Double_t parb[4];
359 backFcn->SetParameters(p3start,p4start,p5start,p6start);
360 backFcn->FixParameter(1,0); // deriv. at zero equal zero
361 //backFcn->FixParameter(2,0); // no 2nd order term
362 backFcn->FixParameter(3,0); // no 3rd order term
363 //halphaOff->Fit ("backFcn","RN");
364 halpha->Fit ("backFcn","RN");
365 backFcn->GetParameters(parb);
366
367 // total fit
368 Double_t par[7];
369 fitFcn->SetParameters(p0start,p1start,p2start,
370 parb[0],parb[1],parb[2],parb[3]);
371
372 fitFcn->FixParameter(2,0.);
373 fitFcn->FixParameter(3,parb[0]);
374 fitFcn->FixParameter(4,parb[1]);
375 fitFcn->FixParameter(5,parb[2]);
376 fitFcn->FixParameter(6,parb[3]);
377 halpha->Fit ("fitFcn","RN");
378
379 // draw the legend
380 TLegend *legend=new TLegend(0.4,0.5,0.93,0.65);
381 legend->SetTextFont(40);
382 legend->SetTextSize(0.04);
383 legend->SetFillColor(19);
384 legend->AddEntry(halpha,"Data","lp");
385 TString strpol("");
386 strpol=Form("Polynom");
387 legend->AddEntry(backFcn,strpol.Data(),"l");
388 legend->AddEntry(fitFcn,"Global Fit","l");
389 //legend->Draw();
390
391 // drawing the functions
392 fitFcn->GetParameters(par);
393 fitFcn->SetRange(-90.,90.);
394 fitFcn->SetLineWidth(2);
395 fitFcn->Draw("same");
396
397 signalFcn->SetParameters(par);
398 signalFcn->SetLineWidth(2);
399 //signalFcn->Draw("same");
400
401 backFcn->SetParameters(&par[3]);
402 backFcn->SetRange(-90,90);
403 backFcn->SetLineWidth(2);
404 backFcn->Draw("same");
405
406
407 //-----------------------------------------------------------------------
408 // calculate significance
409
410 // analytical method
411 Double_t nSig=signalFcn->Integral(0.,acut)/halpha->GetBinWidth(1);
412 Double_t nOffScaled=backFcn->Integral(0.,acut)/halpha->GetBinWidth(1);
413
414 Double_t nOn=nSig+nOffScaled;
415
416 // significance after Li/Ma
417 Double_t theta=1.; // should be = Ton/Toff ???
418 nSig=nOn-nOffScaled;
419 Double_t nOff=nOffScaled/theta;
420
421 Double_t siglima = sqrt(2*(nOn*TMath::Log((1+theta)*nOn/(theta*(nOn+nOff)))
422 +nOff*TMath::Log((1+theta)*(nOff/(nOn+nOff)))));
423
424 // "standard" significance
425 Double_t sigstd = nSig/sqrt(nOn+nOffScaled);
426
427
428 TPaveText *pt = new TPaveText(0.25,0.75,0.56,0.98,"brNDC");
429
430 pt->SetTextFont(60);
431 pt->SetTextSize(0.03);
432 pt->SetTextAlign(12);
433
434 TString strcut1(Form("cuts:"));
435 TString strcut3(Form(" |alpha| < %4.2f ",acut));
436 pt->AddText(strcut3.Data());
437 TString str("");
438 pt->AddText(str);
439
440 TString strsig(Form("NExcess = %d",Int_t(nSig)));
441 pt->AddText(strsig.Data());
442 TString strbak(Form("NOff = %d",Int_t(nOffScaled)));
443 pt->AddText(strbak.Data());
444 pt->AddText(str);
445
446 TString strsg(Form("Significance (Li/Ma) = %4.2f",siglima));
447 pt->AddText(strsg.Data());
448
449 cout<<"results for events with alpha < "<<acut<<" degrees :"<<endl;
450 cout<<" nbackground = "<<nOffScaled<<endl;
451 cout<<" theta = "<<theta<<endl;
452 cout<<" nsignal = "<<nSig<<endl;
453 cout<<" sig. (Li/Ma) = "<<siglima<<endl;
454 cout<<" sig.-std = "<<sigstd<<endl;
455 pt->Draw();
456 gPad->SetGridx();
457 gPad->SetGridy();
458
459 c->Update();
460
461 cout<<endl<<"Optimized Cuts:"<<endl;
462 cout<<" "<<DistLo<<" < dist <"<<DistUp<<endl;
463 cout<<" "<<WidthLo<<" < width <"<<WidthUp<<endl;
464 cout<<" "<<LengthLo<<" < length <"<<LengthUp<<endl;
465
466 return;
467}
468
Note: See TracBrowser for help on using the repository browser.