C++ Analysis::applyCuts方法代码示例

RooDataSet * getDataAndFrac(TString name, TString q2name, TreeReader * mydata, TCut cut, RooRealVar * MM, double * frac, Str2VarMap jpsiPars, double *outnsig)
{
	RooRealVar * cosThetaL = new RooRealVar("cosThetaL","cosThetaL",0.,-1.,1.);
	RooRealVar * cosThetaB = new RooRealVar("cosThetaB","cosThetaB",0.,-1.,1.);
	TCut massCut = "Lb_MassConsLambda > 5590 && Lb_MassConsLambda < 5650";

	Analysis * ana = new Analysis(name+"_mass"+q2name,"Lb",mydata,&cut,MM);
	ana->AddVariable("J_psi_1S_MM");
	ana->AddVariable(cosThetaL);
	ana->AddVariable(cosThetaB);
	RooAbsPdf * mysig = stringToPdf("Gauss","sig",MM,jpsiPars);
	RooAbsPdf * mybkg = stringToPdf("Exp","bkgM",MM);
	RooRealVar * mynsig = new RooRealVar("mynsig","mynsig",50,0,100000);
	RooRealVar * mynbkg = new RooRealVar("mynbkg","mynbkg",10,0,100000);
	RooAbsPdf * Mmodel = new RooAddPdf("MassModel","MassModel",RooArgSet(*mysig,*mybkg),RooArgSet(*mynsig,*mynbkg));
	ana->applyCuts(&cut);
	RooDataSet * data = ana->GetDataSet("-recalc");
	Mmodel->fitTo(*data,Extended(kTRUE));
	
	double sigBkg = mybkg->createIntegral(*MM,NormSet(*MM),Range("Signal"))->getVal();
	double sig = mysig->createIntegral(*MM,NormSet(*MM),Range("Signal"))->getVal();
	double nsig = mynsig->getVal();
	double nbkg = mynbkg->getVal();
	if(frac)
	{
		frac[0] = nsig*sig/(nsig*sig+nbkg*sigBkg);
		frac[1] = frac[0]*TMath::Sqrt( TMath::Power(mynsig->getError()/nsig,2) + TMath::Power(mynbkg->getError()/nbkg,2) );
	}
	TCut mycut = cut + massCut;
	ana->applyCuts(&mycut);

	TCanvas * cc = new TCanvas();
	GetFrame(MM,Mmodel,data,"-nochi2-plotAllComp",30,NULL,0,"cos#theta_{#Lambda}")->Draw();
	cc->Print("M_"+name+"_"+q2name+".pdf");
	if(*outnsig) *outnsig = nsig;
	return ana->GetDataSet("-recalc");
}

void buildBkgPdfs(double q2min, double q2max, TString name, TCut cut, RooAbsPdf ** bkg, RooAbsPdf ** bkgB, TString model = "Poly2")
{
	RooRealVar * cosThetaL = new RooRealVar("cosThetaL","cosThetaL",0.,-1.,1.);
	RooRealVar * cosThetaB = new RooRealVar("cosThetaB","cosThetaB",0.,-1.,1.);

	TString q2name = ((TString)Form("_q2_%4.2f_%4.2f",q2min,q2max)).ReplaceAll(".","");
	string datafile = "/afs/cern.ch/work/p/pluca/weighted/Lmumu/trainingSamples.root";

	string treename = "bkgTestSample";
	if(q2min == 8 && q2max == 11)
	{
		q2name = "_jpsi";
		datafile = "/afs/cern.ch/work/p/pluca/weighted/Lmumu/Lb2Lmumu_CL_NBweighted.root";
		treename = "tree";
	}
	

	Analysis * ana = new Analysis("data"+name+q2name,"Lb",treename,datafile);
	ana->AddVariable(cosThetaL);
	ana->AddVariable(cosThetaB);
	ana->AddVariable("J_psi_1S_MM");
	ana->AddVariable("Lb_MM");

	TCut sideBand = "(Lb_MM > 5690 || Lb_MM < 5550) && weight > 0.8";
	if(q2min == 8 && q2max == 11) sideBand = "Lb_MM > 6000 && weight > 0.8 && TMath::Abs(J_psi_1S_MM - 3096.916) < 92.9";
	TCut q2cut = (TCut)Form("TMath::Power(J_psi_1S_MM/1000,2) > %e && TMath::Power(J_psi_1S_MM/1000,2) < %e",q2min,q2max);
	TCut mycut = sideBand + q2cut + cut;

	ana->applyCuts(&mycut);
	RooDataSet * data = ana->GetDataSet("-recalc");

	if(model=="RooKeyPdf")
	{
		(*bkg) = new RooKeysPdf("bkg"+name,"bkg"+name,*cosThetaL,*data,RooKeysPdf::MirrorBoth,1);
		(*bkgB) = new RooKeysPdf("bkg"+name+"B","bkg"+name+"B",*cosThetaB,*data,RooKeysPdf::MirrorBoth,1);
	}
	else
	{
		(*bkg) = stringToPdf(model,"bkg"+name, cosThetaL);
		(*bkg)->fitTo(*data,PrintLevel(-1));
		fixParams((*bkg),cosThetaL);
		(*bkgB) = stringToPdf(model,"bkg"+name+"B", cosThetaB); 
		(*bkgB)->fitTo(*data,PrintLevel(-1));
		fixParams((*bkgB),cosThetaB);
	}
}

int main(int argc, char **argv)
{
	TString analysis = "Lb2Lmumu";
	bool MC = false;
	TString base = "/afs/cern.ch/work/p/pluca/weighted/Lmumu/";

	if(argc > 1)
	{
		string arg = argv[1];
		if(arg == "MC") MC = true;
	}

	vector< string > novar;
	novar.push_back("Lb_MassCons");
	novar.push_back("Lb_MassConsLambda");
	novar.push_back("Lb_MassConsJpsiLambda");
	novar.push_back("cosTheta");
	novar.push_back("cosThetaL");
	novar.push_back("cosThetaB");
	novar.push_back("phiL");
	novar.push_back("phiB");
	novar.push_back("dphi");
	novar.push_back("cosTheta_TRUE");
	novar.push_back("cosThetaL_TRUE");
	novar.push_back("cosThetaB_TRUE");
	novar.push_back("phiL_TRUE");
	novar.push_back("phiB_TRUE");
	novar.push_back("dphi_TRUE");

	TCut cutJpsi = CutsDef::cutJpsi;
	TCut cutMuMu = CutsDef::cutMuMu_veto;

	TreeReader* treeReader = new TreeReader("tree");

	TString namefile = base + "candLb";
	if(MC) namefile += "_MC";
	namefile += ".root";

	TFile * candFile = new TFile(namefile,"recreate");

	if(!MC) treeReader->AddFile(base+analysis+"_CL_NBweighted.root");
	else treeReader->AddFile(base+analysis+"_MC_Pythia8_NBweighted.root");
	treeReader->Initialize(novar,"except");

	Analysis * anaLbMuMu = new Analysis("Lb2Lmumu","Lb",treeReader,&cutMuMu);

	candFile->cd();
	TTree * candLbMuMu = anaLbMuMu->applyCuts(&addVariables);
	candLbMuMu->Write();
	string tnameMuMu = candLbMuMu->GetName();

	candFile->Close();
	candFile = TFile::Open(namefile,"update");
	TTree * singleCand_LbMuMu = anaLbMuMu->checkMultiple("weight",namefile,tnameMuMu,&randomKill);
	singleCand_LbMuMu->Write();



	if(MC)
	{
		treeReader = new TreeReader("tree");
		treeReader->AddFile(base+"Lb2JpsiL_MC_Pythia8_NBweighted.root");
		treeReader->Initialize(novar,"except");
	}
	Analysis * anaLbJpsi = new Analysis("Lb2JpsiL","Lb",treeReader,&cutJpsi);

	candFile->cd();
	TTree * candLbJpsi = anaLbJpsi->applyCuts(&addVariables);
	candLbJpsi->Write();
	string tnameJpsi = candLbJpsi->GetName();

	candFile->Close();
	candFile = TFile::Open(namefile,"update");	
	TTree * singleCand_LbJpsi = anaLbJpsi->checkMultiple("weight",namefile,tnameJpsi,&randomKill);
	singleCand_LbJpsi->Write();

	candFile->cd();
	TTree * candLbJpsi_reduced = anaLbJpsi->applyCuts(&addVariables,300);
	candLbJpsi_reduced->SetName("candLb2JpsiL_reduced");
	candLbJpsi_reduced->Write();



	if(MC)
	{
		candFile->cd();
		TCut jpsiSwap = cutJpsi + CutsDef::jpsiSwapID;
		TCut mumuSwap = cutMuMu + CutsDef::mumuSwapID;
		TTree * mumuSwapTree = anaLbMuMu->applyCuts(&mumuSwap, false,&addVariables);
		mumuSwapTree->SetName("candLmumuSwap");
		mumuSwapTree->Write();
		TTree * jpsiSwapTree = anaLbJpsi->applyCuts(&jpsiSwap, false, &addVariables);
		jpsiSwapTree->SetName("candJpsiLSwap");
		jpsiSwapTree->Write();

		TreeReader * KSReader = new TreeReader("tree");
		KSReader->AddFile(base+"Bd2JpsiKS_MC12_NBweighted.root");
		KSReader->Initialize(novar,"except");
		TCut cutBdLL = cutJpsi + CutsDef::LLcut;
		TCut cutBdDD = cutJpsi + CutsDef::DDcut;
//.........这里部分代码省略.........