using namespace RooFit; void Blind_fit_2011() { // gROOT->ProcessLine(".L RooCruijff.cxx+"); //Are we still blind? Only change to false after the OK from reviewers /*****************/ bool blind = true; /*****************/ // plot range: const double plotmin = 6.020; const double plotmax = 6.480; // blinded upper/lower sidebands: const double usb = 6.350; const double lsb = 6.150; //Defining Fit Model // --- Observable --- //RooRealVar mes("mB","m_{B} (GeV/c^{2})",6.00,6.500) ; RooRealVar s12("s12","s12 [GeV/c^{2}]^{2}",0.,50.) ; RooRealVar s31("s31","s31 [GeV/c^{2}]^{2}",0.,50.) ; RooRealVar BDT("BDT", "BDT", -20,20); //Define the mass observable RooRealVar mes("mB","m_{B} (GeV/c^{2})",plotmin,plotmax); // --- Parameters --- RooRealVar sigmean("#mu","B^{#pm} mass",6.280,6.270,6.290) ; RooRealVar sigwidth("#sigma","B^{#pm} width",14.3,1,50.) ; // ------- Double Crystal Ball PDF -------- // USING PARAMETERS FROM SIGNAL FIT 2011 RooRealVar meanCB("#mu_{CB}","muCB",6.27818); RooRealVar sigCB("#sigma_{CB}","sigCB",0.01578); RooRealVar nCB("n","n",1.274); RooRealVar alpha("#alpha","#alpha",1.494); RooRealVar meanCB2("#mu2_{CB}","muCB2",6.280,6.270,6.290); // Não é utilizado RooRealVar sigCB2("#sigma2_{CB}","sigCB",0.02242); RooRealVar nCB2("n2","n2",2.27); RooRealVar alpha2("#alpha2","#alpha2",-2.1097); RooRealVar frac("frac","frac",0.452); RooCBShape mycb("CB signal","CB function",mes,meanCB,sigCB,alpha,nCB); RooCBShape mycb2("CB signal 2","CB function 2",mes,meanCB,sigCB2,alpha2,nCB2); RooAddPdf signal("signal","signal model",RooArgList(mycb,mycb2),RooArgList(frac)); // --- Build Background PDF ---- // Polynomial shape \\ RooRealVar shape("P1","P1 Parameter",0,-1000,1000); RooRealVar shape2("P2","P2 Parameter",0,-100,100); // RooPolynomial Comb("Combinatorial","Combinatorial", mes,RooArgList(shape,shape2)); // Exponential shape \\ RooRealVar bkgslope( "bkgslope","background slope",0,-50,50); // aqui era -5,5 RooRealVar shift ("shift","shift",6.079); // MUDANÇA AQUI -> ERA 5.079 RooFormulaVar m_shifted("m_shifted","@0-@1",RooArgSet(mes,shift)); RooGenericPdf Comb("Combinatorial","Comb exp background","exp(@0*@1)", RooArgSet(bkgslope, m_shifted)); // Alternate PartReco, using CONVOLUTION // RooRealVar m0("m0","m_0 (GeV/c^{2})",6.1478) ; RooRealVar slope("slope","slope",-0.52) ; RooGaussian myGaus("resolution","gauss model",mes,RooFit::RooConst(0.0),RooFit::RooConst(0.0185)); RooArgusBG bkgArg("bkgArg","Argus Bkg",mes,m0,slope); RooFFTConvPdf PartReco( "PartReco1", "", mes, bkgArg, myGaus ); PartReco.setBufferFraction(0.5); //RooFFTConvPdf PartReco2( "PartReco2", "", mes, bkgArg, myGaus ); // PartReco2.setBufferFraction(0.3); //RooFFTConvPdf PartReco3( "PartReco3", "", mes, bkgArg, myGaus ); //PartReco3.setBufferFraction(0.3); // ********* READING DATA from NTUPLE and filling RooDataSet object ****** \\ TFile *f1 = TFile::Open("file1.root"); TTree *treedata1 = (TTree*)f1->Get("DecayTree"); RooDataSet data("data","data",RooArgSet(mes,s12,s31)); Double_t Bmass1,s12_DTF1,s31_DTF1; Bool_t B_L0HadronDecision_TOS1; Bool_t B_L0Global_TIS1; Bool_t B_Hlt1TrackAllL0Decision_TOS1; Bool_t B_Hlt2Topo2BodyBBDTDecision_TOS1; Bool_t B_Hlt2Topo3BodyBBDTDecision_TOS1; Bool_t B_Hlt2Topo4BodyBBDTDecision_TOS1; Double_t BDT1, K11, K12,P1; treedata1->SetBranchAddress("BDT",&BDT1); treedata1->SetBranchAddress("d1_MC12TuneV2_ProbNNk",&K11); treedata1->SetBranchAddress("d2_MC12TuneV2_ProbNNk",&K12); treedata1->SetBranchAddress("d3_MC12TuneV2_ProbNNpi",&P1); treedata1->SetBranchAddress("B_L0HadronDecision_TOS",&B_L0HadronDecision_TOS1); treedata1->SetBranchAddress("B_L0Global_TIS",&B_L0Global_TIS1); treedata1->SetBranchAddress("B_Hlt1TrackAllL0Decision_TOS",&B_Hlt1TrackAllL0Decision_TOS1); treedata1->SetBranchAddress("B_Hlt2Topo2BodyBBDTDecision_TOS",&B_Hlt2Topo2BodyBBDTDecision_TOS1); treedata1->SetBranchAddress("B_Hlt2Topo3BodyBBDTDecision_TOS",&B_Hlt2Topo3BodyBBDTDecision_TOS1); treedata1->SetBranchAddress("B_M",&Bmass1); treedata1->SetBranchAddress("s12_DTF",&s12_DTF1); treedata1->SetBranchAddress("s31_DTF",&s31_DTF1); Long64_t nentries1 = treedata1->GetEntries(); cout<<"Entries data: "<GetEntry(i); if (i%100000 == 0) { cout << "Finished processing " << i << " events." << endl; } if ((Bmass1 > 6020) && (Bmass1 < 6480) (B_L0HadronDecision_TOS1 == 1 || B_L0Global_TIS1 == 1) && (B_Hlt1TrackAllL0Decision_TOS1>0) && (B_Hlt2Topo2BodyBBDTDecision_TOS1>0 || B_Hlt2Topo3BodyBBDTDecision_TOS1>0 || B_Hlt2Topo4BodyBBDTDecision_TOS1>0) && (BDT1 > 0.04) && (K11 > 0.4) && (K12 > 0.4) && (P1 > 0.5)) { // With PID // if ((Bmass1 >= 6020) && (Bmass1 <= 6480) (B_L0HadronDecision_TOS1 == 1 || B_L0Global_TIS1 == 1) && (B_Hlt1TrackAllL0Decision_TOS1>0) && (B_Hlt2Topo2BodyBBDTDecision_TOS1>0 || B_Hlt2Topo3BodyBBDTDecision_TOS1>0 || B_Hlt2Topo4BodyBBDTDecision_TOS1>0) && (BDT1 >= 0.04)) { mes = Bmass1*0.001; s12 = s12_DTF1; s31 = s31_DTF1; data.add(RooArgSet(mes,s12,s31)); } } cout << "Finished first ntuple." << endl; TFile *f2 = TFile::Open("file2.root"); TTree *treedata2 = (TTree*)f2->Get("DecayTree"); Bool_t B_L0HadronDecision_TOS2; Bool_t B_L0Global_TIS2; Bool_t B_Hlt1TrackAllL0Decision_TOS2; Bool_t B_Hlt2Topo2BodyBBDTDecision_TOS2; Bool_t B_Hlt2Topo3BodyBBDTDecision_TOS2; Bool_t B_Hlt2Topo4BodyBBDTDecision_TOS2; Double_t Bmass2,s12_DTF2,s31_DTF2; Double_t BDT2, K21, K22,P2; treedata2->SetBranchAddress("BDT",&BDT2); treedata2->SetBranchAddress("d1_MC12TuneV2_ProbNNk",&K21); treedata2->SetBranchAddress("d2_MC12TuneV2_ProbNNk",&K22); treedata2->SetBranchAddress("d3_MC12TuneV2_ProbNNpi",&P2); treedata2->SetBranchAddress("B_M",&Bmass2); treedata2->SetBranchAddress("s12_DTF",&s12_DTF2); treedata2->SetBranchAddress("s31_DTF",&s31_DTF2); treedata2->SetBranchAddress("B_L0HadronDecision_TOS",&B_L0HadronDecision_TOS2); treedata2->SetBranchAddress("B_L0Global_TIS",&B_L0Global_TIS2); treedata2->SetBranchAddress("B_Hlt1TrackAllL0Decision_TOS",&B_Hlt1TrackAllL0Decision_TOS2); treedata2->SetBranchAddress("B_Hlt2Topo2BodyBBDTDecision_TOS",&B_Hlt2Topo2BodyBBDTDecision_TOS2); treedata2->SetBranchAddress("B_Hlt2Topo3BodyBBDTDecision_TOS",&B_Hlt2Topo3BodyBBDTDecision_TOS2); Long64_t nentries2 = treedata2->GetEntries(); cout<<"Entries data: "<GetEntry(i); if (i%100000 == 0) { cout << "Finished processing " << i << " events." << endl; } if ((Bmass2 >= 6020) && (Bmass2 <= 6480) && (B_L0HadronDecision_TOS2 == 1 || B_L0Global_TIS2 == 1) && (B_Hlt1TrackAllL0Decision_TOS2>0) && (B_Hlt2Topo2BodyBBDTDecision_TOS2>0 || B_Hlt2Topo3BodyBBDTDecision_TOS2>0 || B_Hlt2Topo4BodyBBDTDecision_TOS2>0)&& (BDT2 >= 0.04) && (K21 > 0.4) && (K22 > 0.4) && (P2 > 0.5)) { // with PID // if ((Bmass2 >= 6020) && (Bmass2 <= 6480) && (B_L0HadronDecision_TOS2 == 1 || B_L0Global_TIS2 == 1) && (B_Hlt1TrackAllL0Decision_TOS2>0) && (B_Hlt2Topo2BodyBBDTDecision_TOS2>0 || B_Hlt2Topo3BodyBBDTDecision_TOS2>0 || B_Hlt2Topo4BodyBBDTDecision_TOS2>0)&& (BDT2 >= 0.04)) { mes = Bmass2*0.001; s12 = s12_DTF2; s31 = s31_DTF2; data.add(RooArgSet(mes,s12,s31)); } } cout << "Finished second ntuple." << endl; // Yields (blinded) RooRealVar Nsig("Nsig","signal evts",0.5*data.numEntries(),-data.numEntries(),data.numEntries()); RooRealVar Ncomb("Ncomb","comb bkg evts",0.5*data.numEntries(),-data.numEntries(),data.numEntries()); RooRealVar Npartreco("Npartreco","partreco evts", 0.5*data.numEntries(),-data.numEntries(),data.numEntries()); //RooRealVar Npartreco("Npartreco","partreco evts", 0.5*data.numEntries(),0,data.numEntries()); // Make Blinding transforms RooCategory blindingstate("blindingstate","blinding state"); blindingstate.defineType("unblind", 0); blindingstate.defineType("blind", 1); //The blinding transforms. Always use different blinding strings for each parameter. RooUnblindPrecision UBNsig("UBNsig","Unblinded Nsig", "Mickey Mouse", 0.5*data.numEntries(), 2.0, Nsig, blindingstate); RooUnblindPrecision UBNbkg("UBNbkg","Unblinded Nbkg", "Donald Duck", 0.5*data.numEntries(), 2.0, Ncomb, blindingstate); RooUnblindPrecision UBNpartreco("UBNpartreco","Unblinded Npartreco", "Goofy", 0.5*data.numEntries(), 2.0, Npartreco, blindingstate); if(blind){ //Silence RooFit to prevent it printing the fit result: RooMsgService::instance().setSilentMode(kTRUE); //And suppress INFO messages for plot range yields: RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING) ; } // extended likelihood PDFs: RooExtendPdf sige("sige","sige",signal,UBNsig) ; RooExtendPdf Combe2("PartiallyReco","PartiallyReco",PartReco,UBNpartreco) ; RooExtendPdf Combe1("ExpComb","ExpComb",Comb,UBNbkg) ; RooAddPdf bkg("composite background", " ", RooArgList(Combe2,Combe1)); //Combine: //RooAddPdf sum("sum","g+p1",RooArgList(sige,Combe)); // Signal only: to be used for the matched Monte Carlo RooAddPdf sum("sum","dbc+Comb",RooArgList(sige,bkg)); // when using partreco+comb background //RooAddPdf sum("sum","dbc+Comb",RooArgList(sige,Combe1)); //when using only comb background // CALLING THE FITTER: RooFitResult* fitRes = sum.fitTo(data, RooFit::Save(true),RooFit::Strategy(2)); if(blind){ //If we are blind, we should print the fit status information //and the values of the non-blinded params: cout << "FIT IS BLINDED. ONLY PRINTING NUISANCE PARS AND FIT STATUS" << endl; cout << "Fit complete" << endl; cout << "covQual: " << fitRes->covQual() << endl; cout << "EDM: " << fitRes->edm() << endl; cout << "FCN at min: " << fitRes->minNll() << endl; RooArgSet notblinded(bkgslope); notblinded.Print("s"); // RooArgSet partreco(m0,slope); // partreco.Print("s"); } TCanvas *d = new TCanvas("goodblinding","goodblinding",1024,768); RooPlot *gbmassplot = mes.frame(); gbmassplot->SetTitle("2011 Blinded - 0.04 < BDT"); if(blind){ //define lower and upper sidebands, and normalisation: mes.setRange("lowersideband", 6.02, lsb); mes.setRange("uppersideband", usb, 6.48); double sidebands = data.sumEntries("1","uppersideband,lowersideband"); //Now plot only the sidebands: data.plotOn(gbmassplot,CutRange("uppersideband,lowersideband")); sum.plotOn(gbmassplot,Range("uppersideband,lowersideband"), Normalization(sidebands,RooAbsReal::NumEvent), LineColor(kRed)); gbmassplot->SetMinimum(1e-5); //Prevents drawing of 0 entry bins /***************************************************/ }else{ data.plotOn(gbmassplot); sum.plotOn(gbmassplot,Components("signal,Comb"),LineColor(kBlue)); sum.plotOn(gbmassplot,Components("Comb"),LineColor(kRed)); } gbmassplot->Draw(); d->Update(); if(blind){ d->Print("2011_blinded_massplot_goodblinding_004BDT_negative.pdf"); }else{ d->Print("2011_unblinded_massplot_goodblinding.pdf"); } }