| 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 |
|
|---|
| 25 | Double_t acut = 10; // alpha region of excess
|
|---|
| 26 | Bool_t aplot = 1; // make only alpha plot (1) , do sig. optimization and alpha plot (0)
|
|---|
| 27 |
|
|---|
| 28 | //--------------------------------------------------------------------------------------------------------------
|
|---|
| 29 | //--------------------------------------------------------------------------------------------------------------
|
|---|
| 30 |
|
|---|
| 31 | const TMatrix *matOn;
|
|---|
| 32 | const TMatrix *matOff;
|
|---|
| 33 |
|
|---|
| 34 | Int_t isize = 0;
|
|---|
| 35 | Int_t idist = 1;
|
|---|
| 36 | Int_t iwidth = 2;
|
|---|
| 37 | Int_t ilength = 3;
|
|---|
| 38 | Int_t ialpha = 4;
|
|---|
| 39 |
|
|---|
| 40 | //------------------------------------------------------------------------------
|
|---|
| 41 | //------------------------------------------------------------------------------
|
|---|
| 42 |
|
|---|
| 43 | Double_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 |
|
|---|
| 47 | Double_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 |
|
|---|
| 52 | Double_t fitFunction(Double_t *x, Double_t *par){
|
|---|
| 53 | return signal(x,par) + background(x,&par[3]);
|
|---|
| 54 | }
|
|---|
| 55 |
|
|---|
| 56 | //------------------------------------------------------------------------------
|
|---|
| 57 | //------------------------------------------------------------------------------
|
|---|
| 58 | void 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 | //------------------------------------------------------------------------------
|
|---|
| 120 | void 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 |
|
|---|
| 265 | APlot:
|
|---|
| 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 |
|
|---|
