// macro for creating "signal" and "background" samples // with separation by rec.variables, without use of MC info void CreateBkgSigSample(int uxcut=1, int uycut=0, int uzcut=0, int uthetacut=0, int upcut=0, int uallcut=0){ TPad foo; // never remove this line :-))) if(1){ gROOT->SetStyle("Plain"); const Int_t NRGBs = 5; const Int_t NCont = 255; Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 }; Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 }; Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 }; Double_t blue[NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 }; TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont); gStyle->SetNumberContours(NCont); gStyle->SetTitleFont(10*13+2,"xyz"); gStyle->SetTitleSize(0.06, "xyz"); gStyle->SetTitleOffset(1.3,"y"); gStyle->SetTitleOffset(1.3,"z"); gStyle->SetLabelFont(10*13+2,"xyz"); gStyle->SetLabelSize(0.06,"xyz"); gStyle->SetLabelOffset(0.009,"xyz"); gStyle->SetPadBottomMargin(0.16); gStyle->SetPadTopMargin(0.16); gStyle->SetPadLeftMargin(0.16); gStyle->SetPadRightMargin(0.16); gStyle->SetOptTitle(1); gStyle->SetOptStat(1); gROOT->ForceStyle(); gStyle->SetFrameFillColor(0); gStyle->SetFrameFillStyle(0); TGaxis::SetMaxDigits(3); } double pbeam = 15; gROOT->Macro("/panda/pandaroot/macro/lmd/Style_Imported_Style.C"); TString storePath = "/panda/myResults/DPM_pixel/HIMster/all_15/results13122012"; TString in=storePath+"/sumAll.root"; TFile *f = new TFile(in,"READ"); TTree *tAll_sum = (TTree *)f->Get("tAll"); TTree *tSig_sum = (TTree *)f->Get("tSig"); //Try to find some features of signal in general ------------------------------ //(assumes signal contrubution is dominated) //1) copy all important variables in hists TH1F *hpx = new TH1F("hpx","momentum; p_{x}, GeV/c",1e3,-0.1*pbeam,0.1*pbeam); tAll_sum->Project("hpx","px"); TH1F *hpy = new TH1F("hpy","momentum; p_{y}, GeV/c",1e3,-0.1*pbeam,0.1*pbeam); tAll_sum->Project("hpy","py"); TH1F *hpz = new TH1F("hpz","momentum; p_{z}, GeV/c",1e3,0.5*pbeam,1.5*pbeam); tAll_sum->Project("hpz","pz"); TH1F *hp = new TH1F("hp","momentum; p, GeV/c",1e3,0.5*pbeam,1.5*pbeam); tAll_sum->Project("hp","p"); TH1F *hxpca = new TH1F("hxpca","X_{PCA}; x, cm",1e3,-10.,10.); tAll_sum->Project("hxpca","x"); TH1F *hypca = new TH1F("hypca","Y_{PCA}; y, cm",1e3,-10.,10.); tAll_sum->Project("hypca","y"); TH1F *hzpca = new TH1F("hzpca","Z_{PCA}; z, cm",1e3,-0.01,0.01); tAll_sum->Project("hzpca","z"); TH1F *htheta = new TH1F("htheta","#theta; #theta, rad",1e3,0,0.05); tAll_sum->Project("htheta","theta"); TH1F *hphi = new TH1F("hphi","#phi; #phi, rad",1e3,-3.15,3.15); tAll_sum->Project("hphi","phi"); TH1F *hchi2 = new TH1F("hchi2","#chi^{2}; #chi^{2}",3e3,0,1e3); tAll_sum->Project("hchi2","chi2"); //2) fit hists: x, px TF1 *fgaus = new TF1("fgaus", "gaus"); hpx->Fit("fgaus"); double pxmean = fgaus->GetParameter(1); double pxsigma = fgaus->GetParameter(2); hxpca->Fit("fgaus"); double xmean = fgaus->GetParameter(1); double xsigma = fgaus->GetParameter(2); hypca->Fit("fgaus"); double ymean = fgaus->GetParameter(1); double ysigma = fgaus->GetParameter(2); hzpca->Fit("fgaus"); double zmean = fgaus->GetParameter(1); double zsigma = fgaus->GetParameter(2); cout<<"x: "<Divide(2,2); ccontrol->cd(1); hpx->Draw(); ccontrol->cd(2); hxpca->Draw(); ccontrol->cd(3); hpz->Draw(); // ccontrol->cd(4); // hxpca->Draw(); // ccontrol->cd(5); // hypca->Draw(); // ccontrol->cd(6); // hzpca->Draw(); // ccontrol->cd(7); // htheta->Draw(); // ccontrol->cd(8); // hp->Draw(); // ccontrol->cd(9); // hchi2->Draw(); //4) set up limit for approximate sig\bkg separation double xLowCut = xmean-3*xsigma; double xUpCut = xmean+3*xsigma; double yLowCut = ymean-3*ysigma; double yUpCut = ymean+3*ysigma; double zLowCut = zmean-3*zsigma; double zUpCut = zmean+3*zsigma; double thetaLowCut = 0.003; double thetaUpCut = 0.0085; double pLowCut = 0.99*pbeam; double pUpCut = 1.01*pbeam; // double pxLowCut = pxmean-3*pxsigma; // double pxUpCut = pxmean+3*pxsigma; // double pzLowCut = 0.98*pbeam; // double pzUpCut = 1.02*pbeam; cout<<"======= CUTS ======"<SetBranchAddress("theta",&glTheta); tAll_sum->SetBranchAddress("phi",&glPhi); tAll_sum->SetBranchAddress("x",&glXpca); tAll_sum->SetBranchAddress("y",&glYpca); tAll_sum->SetBranchAddress("z",&glZpca); tAll_sum->SetBranchAddress("p",&glP); tAll_sum->SetBranchAddress("px",&glPx); tAll_sum->SetBranchAddress("py",&glPy); tAll_sum->SetBranchAddress("pz",&glPz); tAll_sum->SetBranchAddress("id",&glID); tAll_sum->SetBranchAddress("sumid",&glSumID); tAll_sum->SetBranchAddress("motherid",&glMotherID); tAll_sum->SetBranchAddress("chi2",&glChi2); // ---- Output file ---------------------------------------------------------------- TString out=storePath+"/SigBkgTTree.root"; TFile *fout = new TFile(out,"RECREATE"); // --------------------------------------------------------------------------------- // TTree *tSig_sum = (TTree *)f->Get("tSig"); // TTree *tBkg_sum = (TTree *)f->Get("tBkg"); TTree *tSig_new = new TTree("tSig","Reconstructed variables (seems sig)"); TTree *tBkg_new= new TTree("tBkg","Reconstructed variables (seems bkg)"); TTree *tAll_new= new TTree("tAll","Reconstructed variables"); Double_t glThetaSig,glPhiSig,glXpcaSig,glYpcaSig,glZpcaSig,glPxSig,glPySig,glPzSig,glPSig; Int_t glIDSig,glSumIDSig,glMotherIDSig; Double_t glChi2Sig; int glSigFlag; tSig_new->Branch("theta",&glThetaSig); tSig_new->Branch("phi",&glPhiSig); tSig_new->Branch("x",&glXpcaSig); tSig_new->Branch("y",&glYpcaSig); tSig_new->Branch("z",&glZpcaSig); tSig_new->Branch("p",&glPSig); tSig_new->Branch("px",&glPxSig); tSig_new->Branch("py",&glPySig); tSig_new->Branch("pz",&glPzSig); tSig_new->Branch("id",&glIDSig); tSig_new->Branch("sumid",&glSumIDSig); tSig_new->Branch("motherid",&glMotherIDSig); tSig_new->Branch("chi2",&glChi2Sig); tSig_new->Branch("SigFlag",&glSigFlag); Double_t glThetaBkg,glPhiBkg,glXpcaBkg,glYpcaBkg,glZpcaBkg,glPxBkg,glPyBkg,glPzBkg,glPBkg; Int_t glIDBkg,glSumIDBkg,glMotherIDBkg; Double_t glChi2Bkg; int glBkgFlag; tBkg_new->Branch("theta",&glThetaBkg); tBkg_new->Branch("phi",&glPhiBkg); tBkg_new->Branch("x",&glXpcaBkg); tBkg_new->Branch("y",&glYpcaBkg); tBkg_new->Branch("z",&glZpcaBkg); tBkg_new->Branch("p",&glPBkg); tBkg_new->Branch("px",&glPxBkg); tBkg_new->Branch("py",&glPyBkg); tBkg_new->Branch("pz",&glPzBkg); tBkg_new->Branch("id",&glIDBkg); tBkg_new->Branch("sumid",&glSumIDBkg); tBkg_new->Branch("motherid",&glMotherIDBkg); tBkg_new->Branch("chi2",&glChi2Bkg); tBkg_new->Branch("BkgFlag",&glBkgFlag); Double_t glThetaAll,glXAll,glXBkg,glXSig; tAll_new->Branch("theta",&glThetaAll); tAll_new->Branch("x",&glXAll); tAll_new->Branch("xbkg",&glXBkg); tAll_new->Branch("xsig",&glXSig); int nEv = tAll_sum->GetEntries(); for(int i=0;iGetEvent(i); if(glSumID==4424 && glID==-2212){ glSigFlag = 1; glBkgFlag = -1; } else{ glSigFlag = -1; glBkgFlag = 1; } bool xcut = false ,ycut = false ,zcut = false ,pcut = false ,thetacut = false; if(glXpcaxUpCut) xcut = true; if(glYpcayUpCut) ycut = true; if(glZpcazUpCut) zcut = true; if(glThetathetaUpCut) thetacut = true; if(sqrt(glPx*glPx+glPy*glPy+glPz*glPz)pUpCut) pcut = true; // if(glPxpxUpCut) pxcut = true; // if(glPzpzUpCut) pzcut = true; // bool pxpycut = false; // double pxpyLow = 0.03*0.03; // double pxpyUp = 0.15*0.15; // if((glPx*glPx+glPy*glPy)pxpyUp) pxpycut = true; // if(xcut && pxcut && pzcut){// looks like bkg bool fincut = false; //int uxcut, int uycut, int uzcut, int uthetacut, int upcut, int uallcut if(uxcut>0) fincut = xcut || thetacut; if(uycut>0) fincut = ycut || thetacut; if(uzcut>0) fincut = zcut || thetacut; if(uthetacut>0) fincut = thetacut; if(upcut>0) fincut = pcut || thetacut; if(uallcut>0) fincut = (xcut || ycut || zcut || pcut || thetacut); if(fincut){// looks like bkg // if(thetacut){// looks like bkg glThetaBkg = glTheta; glPhiBkg = glPhi; glXpcaBkg = glXpca; glYpcaBkg = glYpca; glZpcaBkg = glZpca; glPxBkg = glPx; glPyBkg = glPy; glPzBkg = glPz; glPBkg = glP; glIDBkg = glID; glSumIDBkg = glSumID; glMotherIDBkg = glMotherID; glChi2Bkg = glChi2; tBkg_new->Fill(); glThetaAll = glTheta; glXAll = glXpca; glXBkg = glXpca; TBranch* brtheta = tAll_new->GetBranch("theta"); brtheta->Fill(); TBranch* brxall = tAll_new->GetBranch("x"); brxall->Fill(); TBranch* brxbkg = tAll_new->GetBranch("xbkg"); brxbkg->Fill(); } else{ //looks like sig glThetaSig = glTheta; glPhiSig = glPhi; glXpcaSig = glXpca; glYpcaSig = glYpca; glZpcaSig = glZpca; glPxSig = glPx; glPySig = glPy; glPzSig = glPz; glPSig = glP; glIDSig = glID; glSumIDSig = glSumID; glMotherIDSig = glMotherID; glChi2Sig = glChi2; tSig_new->Fill(); glThetaAll = glTheta; glXAll = glXpca; glXSig = glXpca; TBranch* brxsig = tAll_new->GetBranch("xsig"); brxsig->Fill(); TBranch* brtheta = tAll_new->GetBranch("theta"); brtheta->Fill(); TBranch* brxall = tAll_new->GetBranch("x"); brxall->Fill(); } } TH1F *htheta_sig_true = new TH1F("htheta_sig_true","#theta; #theta, mrad",1e2,0,10); tSig_sum->Project("htheta_sig_true","1000*theta"); htheta_sig_true->Write(); TH1F *htheta_sig_new = new TH1F("htheta_sig_new","#theta; #theta, mrad",1e2,0,10); tSig_new->Project("htheta_sig_new","1000*theta"); htheta_sig_new->Write(); TH1F *htheta_sig_only_new = new TH1F("htheta_sig_only_new","#theta; #theta, mrad",1e2,0,10); tSig_new->Project("htheta_sig_only_new","1000*theta","SigFlag>0"); htheta_sig_only_new->Write(); TH1F *htheta_bkg_only_new = new TH1F("htheta_bkg_only_new","#theta; #theta, mrad",1e2,0,10); tSig_new->Project("htheta_bkg_only_new","1000*theta","SigFlag<0"); htheta_bkg_only_new->Write(); TH1F *htheta_sig_new_ratio = new TH1F("htheta_sig_new_ratio","#theta; #theta, mrad",1e2,0,10); tSig_new->Project("htheta_sig_new_ratio","1000*theta"); htheta_sig_new_ratio->Divide(htheta_sig_true); htheta_sig_new_ratio->Write(); TH1F *htheta_sig_only_new_ratio = new TH1F("htheta_sig_only_new_ratio","#theta; #theta, mrad",1e2,0,10); tSig_new->Project("htheta_sig_only_new_ratio","1000*theta","SigFlag>0"); htheta_sig_only_new_ratio->Divide(htheta_sig_true); htheta_sig_only_new_ratio->Write(); TH1F *htheta_bkg_only_new_ratio = new TH1F("htheta_bkg_only_new_ratio","#theta; #theta, mrad",1e2,0,10); tSig_new->Project("htheta_bkg_only_new_ratio","1000*theta","SigFlag<0"); htheta_bkg_only_new_ratio->Divide(htheta_sig_true); htheta_bkg_only_new_ratio->Write(); // TH1F *htheta_diff = new TH1F("htheta_diff","#theta; #theta, rad; #theta_{true} - #theta_{sig}",1e2,0,0.01); // tSig_sum->Project("htheta_diff","theta"); // htheta_diff->Add(htheta_sig_new,-1); // htheta_diff->Write(); TCanvas* c0 = new TCanvas("thetaCan","",200,500,700,800); c0->SetFillColor(0); c0->SetBorderMode(0); c0->Divide(3,2); c0->cd(1); htheta_sig_true->SetLineColor(2); htheta_sig_true->Draw(); htheta_sig_new->SetLineColor(3); // htheta_sig_new->SetFillColor(3); htheta_sig_new->Draw("same"); c0->cd(4); htheta_sig_new_ratio->Draw(); c0->cd(2); htheta_sig_only_new->Draw(); c0->cd(5); htheta_sig_only_new_ratio->Draw(); c0->cd(3); htheta_bkg_only_new->Draw(); c0->cd(6); htheta_bkg_only_new_ratio->Draw(); c0->Write(); tSig_new->Write(); tBkg_new->Write(); tAll_new->Write(); fout->Close(); }