diff --git a/BarbaraTP/eff_plotTP.py b/BarbaraTP/eff_plotTP.py new file mode 100644 index 0000000..f5229a6 --- /dev/null +++ b/BarbaraTP/eff_plotTP.py @@ -0,0 +1,177 @@ +from ROOT import * + +gROOT.SetBatch(kTRUE) + +nbins = 20 +gROOT.ProcessLine(".x mystyle.C") +#from ShipStyle import * +#lhcbstyleSetup() + + +# legTuple (title, h, hcolor) +def makeLegend(legTupl, pos=(0.66,0.78,0.90,0.90)): + leg = TLegend(*pos); + leg.SetLineColor(kWhite) + for lt in legTuple: + leg.AddEntry(lt[0],lt[1],"l"); + + return leg + + +#f = TFile("/afs/cern.ch/work/a/austyuzh/public/4barbara/newGen_15340000_nuBg.root")#newGen/ShipAna_Barbara_nu_th45MeV_LSSegments45MeV_FINALTEST.root") +f = TFile("../../DATA/nuFromYandex/newGen_15340000_nuBg.root") +t = f.Get("t") + + + +#h_veto = TH1F("veto","veto",nbins,0,nbins) +#h_effOff = TH1F("effOff","effOff",nbins,0,nbins) +#h_tot = TH1F("tot","tot",nbins,0,nbins) + + +#flux_nu = 1.09e+12 * 2.e+20/ 4.e+13 +#flux_anti = 6.98e+11 * 2.e+20 / 4.e+13 +#entries_nu = 1.345e+7#200000-2 +#entries_anti = 200000-2 + +#1.61x10^11, +flux_nu = 1.05e+11 * 2.e+20/ 5.e+13#1.61e+11 * 2.e+20/ 4.e+13 +flux_anti = 7.38e+10 * 2.e+20 / 5.e+13#1.02e+11 * 2.e+20 / 4.e+13 +entries_nu = 1.534e+7#200000-2 +entries_anti =175*10e+3 +weight_nu = "(weight * %s / %s)"%(flux_nu,entries_nu) +''' +cutEllipse = "( Sum$(NoB_vtxxSqr/(245.*245.)) + Sum$(NoB_vtxySqr/(490.*490.)) ) < 1." + +docaVtxIp = "Sum$(NoB_DOCA < 30.) && Sum$(NoB_IP0 < 250.) && Sum$(NoB_vtxz > -1958.) && Sum$(NoB_vtxz < 2568.)" + +goodTracks = "Sum$(DaughtersFitConverged==1) && Sum$(DaughtersAlwaysIn==1)" + +trackQ = "Sum$(NoB_MaxDaughtersRedChi2<5.)" + +ndf = "Sum$(NoB_MinDaughtersNdf>10.)" + +hasMuon = "Sum$(Has1Muon1==1) && Sum$(Has1Muon2==1)" + +hasCalo = "Sum$(HasEcal==1)" + +vtxz = "recoed==1"#"Sum$(NoB_vtxz<2500)" +offline_cut = "%s && %s && %s && %s && %s && %s && %s && %s"%(docaVtxIp,goodTracks,trackQ,hasMuon,cutEllipse,hasCalo,ndf,vtxz) +''' + +''' +trackQ = "Sum$(NoB_MaxDaughtersRedChi2<5.)" +cutEllipse = "Sum$(( (NoB_vtxxSqr)/(245.*245.) + (NoB_vtxySqr)/(490.*490.) ) < 1.)" +docaVtxIp = "Sum$(NoB_DOCA < 30.) && Sum$(NoB_IP0 < 250.) && Sum$(NoB_vtxz > -1958. && NoB_vtxz < 2568.)" +vtxz = "fabs(Max$(NoB_vtxz)-Min$(NoB_vtxz))<30" +offline_cut = "%s && %s && %s && %s"%(docaVtxIp,vtxz,trackQ,cutEllipse) +''' + +c1 = "Sum$(NoB_MaxDaughtersRedChi2<5.)" +c2 = "Sum$(( (NoB_vtxxSqr)/(245.*245.) + (NoB_vtxySqr)/(490.*490.) ) < 1.)" +c3 = "Sum$(NoB_DOCA < 30.)" +c4 = "Sum$(NoB_IP0 < 250.)" +c5 = "Sum$(NoB_vtxz > -1958. && NoB_vtxz < 2568.)" +c6 = "TMath::Abs(Max$(NoB_vtxz)-Min$(NoB_vtxz))<30" +offline_cut = "%s && %s && %s && %s && %s && %s"%(c1,c2,c3,c4,c5,c6) + +f_dummy = TFile("dummy.root","RECREATE") + +cuts = {"not-vetoed":"(nRecoed>0 && numberOfHitLSSegments == 0 && strawVetoAny_eff==0 && upstreamVetoAny_eff==0 && RPCany_eff==0 && TrackSyst==1)", + "selected" :"%s"%offline_cut, + "vetoed":"!(nRecoed>0 && numberOfHitLSSegments == 0 && strawVetoAny_eff==0 && upstreamVetoAny_eff==0 && RPCany_eff==0 && TrackSyst==1 && Sum$(NoB_IP0 < 1000.))", + "not-selected" :"!(%s)"%offline_cut, + "recoed":"nRecoed>0"} + +c = {"not-vetoed":kRed, + "vetoed":kBlue, + "selected":kRed, + 'killed':kBlue, + "reco":kBlack} + +n = {'not-vetoed':0, + "vetoed":0, + "selected":0, + 'not-selected':0, + "recoed":0} + +h = {} + +s = 0 +types = ["nChargedPart_fromNu","nNeutrPart_fromNu","nPart_fromNu"] +titles = ["charged", "neutral", ""] +canv = {} +for (ty,tt) in zip(types,titles): + + h[ty] = {"not-vetoed":TH1F("%s_not-vetoed"%ty,"not-vetoed",nbins,0,nbins), + "selected" :TH1F("%s_selected"%ty,"selected",nbins,0,nbins), + "vetoed":TH1F("%s_vetoed"%ty,"vetoed",nbins,0,nbins), + "not-selected" :TH1F("%s_not-selected"%ty,"not-selected",nbins,0,nbins), + "recoed":TH1F("%s_recoed"%ty,"recoed",nbins,0,nbins)} + + + for k in n: + t.Draw("%s>>%s"%(ty,h[ty][k].GetName()),"(%s)*%s"%(cuts[k],weight_nu)) + + h[ty]['res'] = h[ty]['vetoed'].Clone() + h[ty]['res'].Divide(h[ty]['recoed']) + h[ty]['res'].SetLineColor(kBlack) + h[ty]['res'].SetMarkerColor(kBlack) + + yaxis = h[ty]['res'].GetYaxis() + xaxis = h[ty]['res'].GetXaxis() + yaxis.SetRangeUser(0.96,1.003) + yaxis.SetDecimals() + yaxis.SetTitle("efficiency") + yaxis.SetTitleOffset(1.1) + xaxis.SetTitle("%s particles / event"%tt) + h[ty]['res2'] = h[ty]['not-selected'].Clone() + h[ty]['res2'].Divide(h[ty]['recoed']) + h[ty]['res2'].SetLineColor(kBlue) + h[ty]['res2'].SetMarkerColor(kBlue) + h[ty]['res2'].SetLineStyle(2) + + + canv[ty] = TCanvas("c_efficiency_%s"%ty,"c_efficiency_%s"%ty) + h[ty]['res'].Draw('') + h[ty]['res2'].Draw('same') + + legTuple = [(h[ty]['res'],'veto',kBlack),(h[ty]['res2'],'rejection','kBlue')] + leg = makeLegend(legTuple, pos=(0.30,0.34,0.56,0.47))#(0.34,31,58,43)) + leg.Draw("same") + + canv[ty].Print("eff/"+canv[ty].GetName()+".pdf") +assert(False) + + + +#h['res'].Multiply(h['veto']) +#h['res'].Draw() +#s = h['res'].GetSum() +#print s + + + + +''' +s = 0 +for i in xrange(nbins): + for k in n: + print k, i + print t.Draw("event>>h_%s_%s"%(k,i),"(%s && nChargedPart_fromNu==%s)*%s"%(cuts[k],i,weight_nu)) + #t.Draw("event","(%s && nChargedPart_fromNu==%s)*%s"%(cuts[k],i,weight_nu)) + #n[k] = float(t.GetHistogram().GetSum()) + print i, k, t.GetHistogram().GetSum() + print t.GetHistogram().GetSum() + print t.GetHistogram().GetSum() + assert(False) + h_veto.SetBinContent(i+1,n['veto']) + if n['reco']>0: + print n["off"]/n["reco"], n['veto']," -> ", n["off"]/n["reco"]*n["veto"] + # h_effOff.SetBinContent(i+1,n["off"]/n["reco"]) + #h_tot.SetBinContent(i+1,n["off"]/n["reco"]*n["veto"]) + s += n["off"]/n["reco"]*n["veto"] +#print h_tot.GetSum() +#h_tot.Draw() +print s +''' diff --git a/BarbaraTP/int_plotTP.py b/BarbaraTP/int_plotTP.py new file mode 100644 index 0000000..de63b4b --- /dev/null +++ b/BarbaraTP/int_plotTP.py @@ -0,0 +1,230 @@ +from ROOT import * + +gROOT.SetBatch(kTRUE) + +'''f = TFile("/afs/cern.ch/work/a/austyuzh/public/4barbara/newGen_15340000_nuBg.root")#newGen/ShipAna_Barbara_nu_th45MeV_LSSegments45MeV_FINALTEST.root") +t = f.Get("t") + +flux_nu = 1.05e+11 * 2.e+20/ 5.e+13#1.61e+11 * 2.e+20/ 4.e+13 +flux_anti = 7.38e+10 * 2.e+20 / 5.e+13#1.02e+11 * 2.e+20 / 4.e+13 +entries_nu = 1.345e+7#200000-2 +entries_anti = 200000-2 +weight_nu = "(weight * %s / %s)"%(flux_nu,entries_nu) +''' + +what = "nu" + + +if what == "anti": + f = TFile("/afs/cern.ch/work/a/austyuzh/public/4barbara/newGen_15340000_anti-nuBg.root") + t = f.Get("t") +else: + #f = TFile("/afs/cern.ch/work/a/austyuzh/public/4barbara/newGen_15340000_nuBg.root")#newGen/ShipAna_Barbara_nu_th45MeV_LSSegments45MeV_FINALTEST.root") + f = TFile("../../DATA/nuFromYandex/newGen_15340000_nuBg.root") + t = f.Get("t") + + + +#b = TBrowser() + +flux_nu = 1.05e+11 * 2.e+20/ 5.e+13#1.61e+11 * 2.e+20/ 4.e+13 +flux_anti = 7.38e+10 * 2.e+20 / 5.e+13#1.02e+11 * 2.e+20 / 4.e+13 +entries_nu = 1.534e+7#1.345e+7#200000-2 +entries_anti =175*10e+4 +if what=="anti": + + weight_nu = "(weight * %s / %s)"%(flux_anti,entries_anti) + +else: + weight_nu = "(weight * %s / %s)"%(flux_nu,entries_nu) + + +res = {'recoed':{}, 'not-vetoed':{}, 'selected':{}} +entries = t.GetEntries() +bbb = {} +for entry in xrange(entries): + t.GetEntry(entry) + intEl = t.nuIntNumSimpl + if intEl in [12,8,11,0,9,16]: + inter = 'nuDet' + #res['recoed'][inter]=0 + #res['not-vetoed'][inter]=0 + #res['selected'][inter]=0 + elif intEl in [4,10,7,2]: + inter = "vesselWalls" + #res['recoed'][inter]=0 + #res['not-vetoed'][inter]=0 + #res['selected'][inter]=0 + elif intEl in [14,18,19,17]: + inter = 'TrackSyst' + #res['recoed'][inter]=0 + #res['not-vetoed'][inter]=0 + #res['selected'][inter]=0 + elif intEl in [15,]: + inter = 'strawVeto' + #res['recoed'][inter]=0 + #res['not-vetoed'][inter]=0 + #res['selected'][inter]=0 + elif intEl in [5,1]: + inter = 'cave' + #res['recoed'][inter]=0 + #res['not-vetoed'][inter]=0 + #res['selected'][inter]=0 + elif intEl in [6,3,13]: + inter = 'vesselLids' + #res['recoed'][inter]=0 + #res['not-vetoed'][inter]=0 + #res['selected'][inter]=0 + + val = (1*t.weight*flux_nu/entries_nu) + + if inter in res['recoed']: + res['recoed'][inter]+=val + else: + res['recoed'][inter]=val + + if (t.numberOfHitLSSegments == 0 and t.strawVetoAny_eff==0 and t.upstreamVetoAny_eff==0 and t.RPCany_eff==0 ): + flag = False + for p in t.NoB_IP0: + if p < 1000.: + flag = True + if flag: + if inter in res['not-vetoed']: + res['not-vetoed'][inter]+=val + else: + res['not-vetoed'][inter]=val + + track_f = False + cutElipse_f = False + docaVtxIp_f = False + doca_f = False + ip_f = False + vessel_f = False + vtx_f = False + for a in t.NoB_MaxDaughtersRedChi2: + if a<5.: + track_f=True + if '1' in bbb: + bbb['1']+=val + else: + bbb['1']=val + assert(len(t.NoB_vtxxSqr)==len(t.NoB_vtxySqr)) + for (a,b) in zip(t.NoB_vtxxSqr,t.NoB_vtxySqr): + if (a/(245.*245.) + b/(490.*490.) ) < 1.: + cutElipse_f = True + #if t.event == 453 and t.jobID==2750601: + #print t.event, t.jobID, t.event+t.jobID, a,b + #print "all" + #for (a,b) in zip(t.NoB_vtxxSqr,t.NoB_vtxySqr): + # print a, b + #assert(False) + if '2' in bbb: + bbb['2']+=val + else: + bbb['2']=val + + for a in t.NoB_DOCA: + if a < 30.: + doca_f=True + if '3' in bbb: + bbb['3']+=val + else: + bbb['3']=val + for a in t.NoB_IP0: + if a< 250.: + ip_f = True + if '4' in bbb: + bbb['4']+=val + else: + bbb['4']=val + for a in t.NoB_vtxz: + if a > -1958. and a<2568: + vessel_f=True + if '5' in bbb: + bbb['5']+=val + else: + bbb['5']=val + if fabs(max(t.NoB_vtxz)-min(t.NoB_vtxz))<30: + vtx_f = True + if '6' in bbb: + bbb['6']+=val + else: + bbb['6']=val + + if track_f and cutElipse_f and doca_f and ip_f and vessel_f and vtx_f: + if inter in res['selected']: + res['selected'][inter]+=val + else: + res['selected'][inter]=val + +''' +trackQ = "Sum$(NoB_MaxDaughtersRedChi2<5.)" +cutEllipse = "Sum$(( (NoB_vtxxSqr)/(245.*245.) + (NoB_vtxySqr)/(490.*490.) ) < 1.))" +docaVtxIp = "Sum$(NoB_DOCA < 30.) && Sum$(NoB_IP0 < 250.) && Sum$(NoB_vtxz > -1958. && NoB_vtxz < 2568.)" +vtxz = "fabs(Max$(NoB_vtxz)-Min$(NoB_vtxz))<30" +offline_cut = "%s && %s && %s && %s"%(docaVtxIp,vtxz,trackQ,cutEllipse) +''' + +c1 = "Sum$(NoB_MaxDaughtersRedChi2<5.)" +c2 = "Sum$(( (NoB_vtxxSqr)/(245.*245.) + (NoB_vtxySqr)/(490.*490.) ) < 1.)" +c3 = "Sum$(NoB_DOCA < 30.)" +c4 = "Sum$(NoB_IP0 < 250.)" +c5 = "Sum$(NoB_vtxz > -1958. && NoB_vtxz < 2568.)" +c6 = "TMath::Abs(Max$(NoB_vtxz)-Min$(NoB_vtxz))<30" +offline_cut = "%s && %s && %s && %s && %s && %s"%(c1,c2,c3,c4,c5,c6) + +t.Draw("event","(%s)*%s"%(c1,weight_nu)) +print bbb['1'],t.GetHistogram().GetSum() + +t.Draw("event","(%s)*%s"%(c2,weight_nu)) +print bbb['2'],t.GetHistogram().GetSum() + +t.Draw("event","(%s)*%s"%(c3,weight_nu)) +print bbb['3'],t.GetHistogram().GetSum() + +t.Draw("event","(%s)*%s"%(c4,weight_nu)) +print bbb['4'],t.GetHistogram().GetSum() + +t.Draw("event","(%s)*%s"%(c5,weight_nu)) +print bbb['5'],t.GetHistogram().GetSum() + +t.Draw("event","(%s)*%s"%(c6,weight_nu)) +print bbb['6'],t.GetHistogram().GetSum() + + +print res['selected'] + +a = 0 +tot = {} +for tt in res: + a = 0 + for k in res[tt]: + a+=res[tt][k] + tot[tt]=a + +t.Draw("event","(%s)*%s"%(offline_cut,weight_nu)) +print tot['selected'],t.GetHistogram().GetSum() + +for k in res['recoed']: + if not k in res['not-vetoed']: + res['not-vetoed'][k] = 0 + if not k in res['selected']: + res['selected'][k] = 0 + + +print "\\begin{table}" +print " \\centering" +print " \\begin{tabular}{lrrr}" +print " \\hline" +print " detector & reconstructed (\\%) & not vetoed (\\%) & selected (\\%) \\\\" +print " \\hline" +print " $\\nu$ detector & %.0f (%.1f) & %.0f (%.1f)& %.0f (%.1f)\\\\"%(res['recoed']['nuDet'],res['recoed']['nuDet']/tot['recoed']*100.,res['not-vetoed']['nuDet'],res['not-vetoed']['nuDet']/tot['not-vetoed']*100.,res['selected']['nuDet'],res['selected']['nuDet']/tot['selected']*100.) +print " vessel lids & %.0f (%.1f)& %.0f (%.1f)& %.0f (%.1f)\\\\"%(res['recoed']['vesselLids'],res['recoed']['vesselLids']/tot['recoed']*100.,res['not-vetoed']['vesselLids'],res['not-vetoed']['vesselLids']/tot['not-vetoed']*100.,res['selected']['vesselLids'],res['selected']['vesselLids']/tot['selected']*100.) +print " vessel walls & %.0f (%.1f) & %.0f (%.1f)& %.0f (%.1f)\\\\"%(res['recoed']["vesselWalls"],res['recoed']["vesselWalls"]/tot['recoed']*100.,res['not-vetoed']["vesselWalls"],res['not-vetoed']["vesselWalls"]/tot['not-vetoed']*100.,res['selected']["vesselWalls"],res['selected']["vesselWalls"]/tot['selected']*100.) +#print " straw veto & %.0f (%.0f) & %.0f (%.0f)& %.0f (%.0f)\\\\"%(res['recoed']['strawVeto'],res['recoed']['strawVeto']/tot['recoed'],res['not-vetoed']['strawVeto'],res['not-vetoed']['strawVeto']/tot['not-vetoed'],res['selected']['strawVeto'],res['selected']['strawVeto']/tot['selected']) +print " tracking system & %.0f (%.1f) & %.1f (%.1f)& %.0f (%.1f)\\\\"%(res['recoed']['TrackSyst'],res['recoed']['TrackSyst']/tot['recoed']*100.,res['not-vetoed']['TrackSyst'],res['not-vetoed']['TrackSyst']/tot['not-vetoed']*100.,res['selected']['TrackSyst'],res['selected']['vesselLids']/tot['selected']*100.) +print " cave & %.0f (%.1f) & %.0f (%.1f)& %.0f (%.1f)\\\\"%(res['recoed']['cave'],res['recoed']['cave']/tot['recoed']*100.,res['not-vetoed']['cave'],res['not-vetoed']['cave']/tot['not-vetoed']*100.,res['selected']['cave'],res['selected']['cave']/tot['selected']*100.) +print " \\hline" +print " \\textbf{Total} & \\textbf{%.0f} & \\textbf{%.0f} & \\textbf{%.0f} \\\\"%(tot['recoed'],tot['not-vetoed'],tot['selected']) +print " \\end{tabular}" +print "\\end{table}" diff --git a/BarbaraTP/mystyle.C b/BarbaraTP/mystyle.C new file mode 100644 index 0000000..a903c2f --- /dev/null +++ b/BarbaraTP/mystyle.C @@ -0,0 +1,179 @@ +{ +//////////////////////////////////////////////////////////////////// +// PURPOSE: +// +// This macro defines a reasonable style for (black-and-white) +// "publication quality" ROOT plots. The default settings contain +// many features that are either not desirable for printing on white +// paper or impair the general readibility of plots. +// +// USAGE: +// +// Simply include the line +// gROOT->ProcessLine(".x $LHCBSTYLE/root/lhcbstyle.C"); +// at the beginning of your root macro. +// +// SOME COMMENTS: +// +// Statistics and fit boxes: +// +// "Decorative" items around the histogram are kept to a minimum. +// In particular there is no box with statistics or fit information. +// You can easily change this either by editing your private copy +// of this style file or by calls to "gStyle" in your macro. +// For example, +// gStyle->SetOptFit(1011); +// will add some fit information. +// +// Font: +// +// The font is chosen to be 62, i.e.helvetica-bold-r-normal with +// precision 2. Font is of course a matter of taste, but most people +// will probably agree that Helvetica bold gives close to optimal +// readibility in presentations. It appears to be the ROOT default, +// and since there are still some features in ROOT that simply won't +// respond to any font requests, it is the wise choice to avoid +// ugly font mixtures on the same plot... The precision of the font (2) +// is chosen in order to have a rotatable and scalable font. Be sure +// to use true-type fonts! I.e. +// Unix.*.Root.UseTTFonts: true in your .rootrc file. +// +// "Landscape histograms": +// +// The style here is designed for more or less quadratic plots. +// For very long histograms, adjustements are needed. For instance, +// for a canvas with 1x5 histograms: +// TCanvas* c1 = new TCanvas("c1", "L0 muons", 600, 800); +// c1->Divide(1,5); +// adaptions like the following will be needed: +// gStyle->SetTickLength(0.05,"x"); +// gStyle->SetTickLength(0.01,"y"); +// gStyle->SetLabelSize(0.15,"x"); +// gStyle->SetLabelSize(0.1,"y"); +// gStyle->SetStatW(0.15); +// gStyle->SetStatH(0.5); +// +//////////////////////////////////////////////////////////////////// + +gROOT->Reset(); + +cout << "executing lhcbStyle.C:" << endl; +cout << " " << endl; +cout << " " << endl; +cout << " $ $ $ $$$ $ " << endl; +cout << " $ $ $ $ $ " << endl; +cout << " $ $$$$$ $ $$$ " << endl; +cout << " $ $ $ $ $ $ " << endl; +cout << " $$$$$ $ $ $$$ $$$ " << endl; +cout << " " << endl; +cout << " LHCb ROOT style file " << endl; +cout << " " << endl; +cout << +" Problems, suggestions, contributions to Thomas.Schietinger@cern.ch" + << endl; +cout << " " << endl; + +TStyle *lhcbStyle= new TStyle("lhcbStyle","LHCb official plots style"); + +// use helvetica-bold-r-normal, precision 2 (rotatable) + Int_t lhcbFont = 42; //62; +// line thickness +Double_t lhcbWidth = 2; + +// use plain black on white colors +lhcbStyle->SetFrameBorderMode(0); +lhcbStyle->SetCanvasBorderMode(0); +lhcbStyle->SetPadBorderMode(0); +lhcbStyle->SetPadColor(0); +lhcbStyle->SetCanvasColor(0); +lhcbStyle->SetStatColor(0); +lhcbStyle->SetPalette(1); +lhcbStyle->SetTitleColor(1); + lhcbStyle->SetFillColor(1);//0 + lhcbStyle->SetFillStyle(0); + +// set the paper & margin sizes +lhcbStyle->SetPaperSize(20,26); +lhcbStyle->SetPadTopMargin(0.05); +lhcbStyle->SetPadRightMargin(0.05); // increase for colz plots!! +lhcbStyle->SetPadBottomMargin(0.16); +lhcbStyle->SetPadLeftMargin(0.14); + +// use large fonts +lhcbStyle->SetTextFont(lhcbFont); +lhcbStyle->SetTextSize(0.08); +lhcbStyle->SetLabelFont(lhcbFont,"x"); +lhcbStyle->SetLabelFont(lhcbFont,"y"); +lhcbStyle->SetLabelFont(lhcbFont,"z"); +lhcbStyle->SetLabelSize(0.05,"x"); +lhcbStyle->SetLabelSize(0.05,"y"); +lhcbStyle->SetLabelSize(0.05,"z"); +lhcbStyle->SetTitleFont(lhcbFont); +lhcbStyle->SetTitleFont(lhcbFont,"y"); +lhcbStyle->SetTitleFont(lhcbFont,"x"); +lhcbStyle->SetTitleFont(lhcbFont,"z"); +lhcbStyle->SetTitleSize(0.06,"x"); +lhcbStyle->SetTitleSize(0.06,"y"); +lhcbStyle->SetTitleSize(0.06,"z"); + +// use bold lines and markers +lhcbStyle->SetLineWidth(lhcbWidth); +lhcbStyle->SetFrameLineWidth(lhcbWidth); +lhcbStyle->SetHistLineWidth(lhcbWidth); +lhcbStyle->SetFuncWidth(lhcbWidth); +lhcbStyle->SetGridWidth(lhcbWidth); +lhcbStyle->SetLineStyleString(2,"[12 12]"); // postscript dashes +lhcbStyle->SetMarkerStyle(8); +lhcbStyle->SetMarkerSize(1.5); + +// label offsets +lhcbStyle->SetLabelOffset(0.015); + +// by default, do not display histogram decorations: +lhcbStyle->SetOptStat(0); +//lhcbStyle->SetOptStat(1110); // show only nent, mean, rms +lhcbStyle->SetOptTitle(0); +lhcbStyle->SetOptFit(0); +//lhcbStyle->SetOptFit(1011); // show probability, parameters and errors + +// look of the statistics box: +lhcbStyle->SetStatBorderSize(1); +lhcbStyle->SetStatFont(lhcbFont); +lhcbStyle->SetStatFontSize(0.05); +lhcbStyle->SetStatX(0.9); +lhcbStyle->SetStatY(0.9); +lhcbStyle->SetStatW(0.25); +lhcbStyle->SetStatH(0.15); + +// put tick marks on top and RHS of plots +lhcbStyle->SetPadTickX(1); +lhcbStyle->SetPadTickY(1); + +// histogram divisions: only 5 in x to avoid label overlaps +lhcbStyle->SetNdivisions(505,"x"); +lhcbStyle->SetNdivisions(510,"y"); + +gROOT->SetStyle("lhcbStyle"); +gROOT->ForceStyle(); + +TPaveText *lhcbName = new TPaveText(0.65,0.8,0.9,0.9,"BRNDC"); +lhcbName->SetFillColor(0); +lhcbName->SetTextAlign(12); +lhcbName->SetBorderSize(0); +lhcbName->AddText("LHCb"); + +TText *lhcbLabel = new TText(); +lhcbLabel->SetTextFont(lhcbFont); +lhcbLabel->SetTextColor(1); +lhcbLabel->SetTextSize(0.04); +lhcbLabel->SetTextAlign(12); + +TLatex *lhcbLatex = new TLatex(); +lhcbLatex->SetTextFont(lhcbFont); +lhcbLatex->SetTextColor(1); +lhcbLatex->SetTextSize(0.04); +lhcbLatex->SetTextAlign(12); +lhcbStyle->SetLegendFillColor(kWhite); + +} + diff --git a/BarbaraTP/weight_plotTP.py b/BarbaraTP/weight_plotTP.py new file mode 100644 index 0000000..65035c9 --- /dev/null +++ b/BarbaraTP/weight_plotTP.py @@ -0,0 +1,103 @@ +from ROOT import * + +gROOT.SetBatch(kTRUE) + +pdg = TDatabasePDG.Instance() +gROOT.ProcessLine(".x mystyle.C") +#from ShipStyle import * +#lhcbstyleSetup() +#gROOT.SetBatch(kTRUE) + + +what = "nu" + + +if what == "anti": + f = TFile("/afs/cern.ch/work/a/austyuzh/public/4barbara/newGen_15340000_anti-nuBg.root") + t = f.Get("t") +else: + #f = TFile("/afs/cern.ch/work/a/austyuzh/public/4barbara/newGen_15340000_nuBg.root")#newGen/ShipAna_Barbara_nu_th45MeV_LSSegments45MeV_FINALTEST.root") + f = TFile("../../DATA/nuFromYandex/newGen_15340000_nuBg.root") + t = f.Get("t") + + + +#b = TBrowser() + +flux_nu = 1.05e+11 * 2.e+20/ 5.e+13#1.61e+11 * 2.e+20/ 4.e+13 +flux_anti = 7.38e+10 * 2.e+20 / 5.e+13#1.02e+11 * 2.e+20 / 4.e+13 +entries_nu = 1.534e+7#200000-2 +entries_anti =175*10e+3 +if what=="anti": + + weight_nu = "(weight * %s / %s)"%(flux_anti,entries_anti) + +else: + weight_nu = "(weight * %s / %s)"%(flux_nu,entries_nu) + +'''flux_nu = 1.61e+11 * 2.e+20/ 4.e+13 +flux_anti = 1.02e+11 * 2.e+20 / 4.e+13 +entries_nu = 1.345e+7#200000-2 +entries_anti = 200000-2 +weight_nu = "(weight * %s / %s)"%(flux_nu,entries_nu) +''' +n = t.Draw("event", "(nuE>10)*%s"%weight_nu) +print "nuE>10: entries=%.0f weighted=%.0f"%(n, t.GetHistogram().GetSum()) +n = t.Draw("event", "(nuE>4 && nuE<10)*%s"%weight_nu) +print "4>h3", "(nuE<4)*%s"%weight_nu) +print "nuE<4: entries=%.0f weighted=%.0f"%(n, t.GetHistogram().GetSum()) + + +assert(False) + +def makeText(textes): + pave = TPaveText(0.63, 0.72, 0.92, 0.84,"TNDC") + for tx in textes: + pave.AddText("%s"%tx) + pave.SetTextAlign(11) + return pave + +def makeLegend(legTuple, pos=(0.66,0.78,0.90,0.90)): + leg = TLegend(*pos); + leg.SetLineColor(kWhite) + for lt in legTuple: + leg.AddEntry(lt[0],lt[1],"l"); + + return leg + + +def makePlot(name,var,cut,logY=False,legend=True): + c = TCanvas("c_%s"%name,"c_%s"%name) + c.SetLogy(logY) + t.Draw("%s"%var, "%s"%cut) + h = t.GetHistogram() + if legend: + pave = makeLegend([[h,cut],],pos=(0.56,0.69,0.90,0.90)) + pave.Draw("same") + c.Print(c.GetName()+".pdf") + return c,h,pave + +c0,h0,p0 = makePlot("recoed",weight_nu,"",True) +c1,h1,p1 = makePlot("ndf_25",weight_nu,"NoB_MinDaughtersNdf>25",True) +c2,h2,p2 = makePlot("acceptance",weight_nu,"NoB_MinDaughtersNdf>25 && Sum$( (NoB_vtxxSqr)/(245.*245.) + (NoB_vtxySqr)/(490.*490.) ) < 1. && Sum$(NoB_vtxz > -1958. && NoB_vtxz < 2568.)",True) +c3,h3,p3 = makePlot("IP0<10m",weight_nu,"NoB_IP0<1000",True) +c8, h8, p8 = makePlot("vtxSec",weight_nu,"TMath::Abs(Max$(NoB_vtxz)-Min$(NoB_vtxz))<30",True) + +c4,h4,p4 = makePlot("ndf_weight","%s:NoB_MinDaughtersNdf"%weight_nu,"",False,False) +c5,h5,p5 = makePlot("vtxz_weight","%s:vtxz"%weight_nu,"vtxz>-100",False,False) +c6,h6,p6 = makePlot("IP_weight","%s:NoB_IP0"%weight_nu,"",False,False) +c7,h7,p7 = makePlot("nrecoed_weight","%s:nRecoed"%weight_nu,"nRecoed<30",False,False) +c9, h9, p9 = makePlot("vtxSec_weight","%s:TMath::Abs(Max$(NoB_vtxz)-Min$(NoB_vtxz))"%weight_nu,"TMath::Abs(Max$(NoB_vtxz)-Min$(NoB_vtxz))",False,False) + +c10, h10, p10 = makePlot("notVetoed",weight_nu,"(nRecoed>0 && numberOfHitLSSegments == 0 && strawVetoAny_eff==0 && upstreamVetoAny_eff==0 && RPCany_eff==0 && TrackSyst==1)",True) + +c11, h11, p11 = makePlot("notVetoed_IP10",weight_nu,"(nRecoed>0 && numberOfHitLSSegments == 0 && strawVetoAny_eff==0 && upstreamVetoAny_eff==0 && RPCany_eff==0 && TrackSyst==1 && NoB_IP0<1000)",True) + +c12, h12, p12 = makePlot("notVetoed_acceptance_IP",weight_nu,"(nRecoed>0 && numberOfHitLSSegments == 0 && strawVetoAny_eff==0 && upstreamVetoAny_eff==0 && RPCany_eff==0 && TrackSyst==1 && NoB_IP0<1000) && NoB_MinDaughtersNdf>25 && Sum$( (NoB_vtxxSqr)/(245.*245.) + (NoB_vtxySqr)/(490.*490.) ) < 1. && Sum$(NoB_vtxz > -1958. && NoB_vtxz < 2568.)",True) + +c13, h13, p13 = makePlot("notVetoed_acceptance_IP",weight_nu,"(nRecoed>0 && numberOfHitLSSegments == 0 && strawVetoAny_eff==0 && upstreamVetoAny_eff==0 && RPCany_eff==0 && TrackSyst==1 ) && NoB_MinDaughtersNdf>25 && Sum$( (NoB_vtxxSqr)/(245.*245.) + (NoB_vtxySqr)/(490.*490.) ) < 1. && Sum$(NoB_vtxz > -1958. && NoB_vtxz < 2568.)",True) + +c13, h13, p13 = makePlot("notVetoed_acceptance",weight_nu,"(nRecoed>0 && numberOfHitLSSegments == 0 && strawVetoAny_eff==0 && upstreamVetoAny_eff==0 && RPCany_eff==0 && TrackSyst==1 && NoB_IP0<1000) && NoB_MinDaughtersNdf>25 && Sum$( (NoB_vtxxSqr)/(245.*245.) + (NoB_vtxySqr)/(490.*490.) ) < 1. && Sum$(NoB_vtxz > -1958. && NoB_vtxz < 2568.)",True) + +assert(False) diff --git a/TP.py b/TP.py new file mode 100644 index 0000000..0389708 --- /dev/null +++ b/TP.py @@ -0,0 +1,70 @@ +from distribsForHNLandBG_byEvent import * + + + +## non funziona +#from ShipStyle import * +#lhcbstyleSetup() + +r.gROOT.ProcessLine(".x newGen/mystyle.C") +r.gROOT.SetBatch(r.kTRUE) + +#accPlotsTP('pimu') + +if not os.path.exists('accPlotsTP'): + os.system('mkdir accPlotsTP') + +cSaver = {} + +# BG stuff +bg_t = studies['nu']['data'] +bg_geo = studies['nu']['geo'] +bg_ntot = studies['nu']['ntot'] +bg_tc = cutsWithDraw() +bg_tc.setNoB('NoB_') +bg_zmin = bg_geo['Veto_5']['z']['pos']+bg_geo['Veto_5']['z']['dim'] +bg_zmax = bg_geo['Tr1_1']['z']['pos']-bg_geo['Tr1_1']['z']['dim'] + +for study in ['pimu', 'mumunu', 'tiny']: + #r.gROOT.ProcessLine(".x newGen/mystyle.C") + t = studies[study]['data'] + geo = studies[study]['geo'] + ntot = studies[study]['ntot'] + tc = cutsWithDraw() + zmin = geo['Veto_5']['z']['pos']+geo['Veto_5']['z']['dim'] + zmax = geo['Tr1_1']['z']['pos']-geo['Tr1_1']['z']['dim'] + # Latest requests... making Hans and Thomas happy :) + additionalCuts = "&&"+tc.ndf(25)[0]+"&&"+tc.redChi2(5.)[0]#+"&&"+tc.fiducial(zmin, zmax, 245., 495.)[0]+"&&"+tc.doca(2.)[0]#+"&&"+tc.nRecoed(1)[0] + prefix = 'NoB_' + # Not vetoed + ip = plotSigBg(cSaver, study+"-IP0", t, ntot, bg_t, prefix+"IP0", tc.notVetoed(1000.)[0]+additionalCuts, "Impact parameter to target [cm]", (0.27, 0.75, 0.50, 0.92, "Not vetoed"), ("y")) + iplin = plotSigBg(cSaver, study+"-IP0-lin", t, ntot, bg_t, prefix+"IP0", tc.notVetoed(1000.)[0]+additionalCuts, "Impact parameter to target [cm]", (0.64, 0.75, 0.87, 0.92, "Not vetoed"), ()) + chi2 = plotSigBg(cSaver, study+"-chi2", t, ntot, bg_t, prefix+"MaxDaughtersRedChi2", tc.notVetoed(1000.)[0], "Maximum daughters #chi^{2}/ndf", (0.64, 0.72, 0.87, 0.89, "Not vetoed"), ("y")) + doca = plotSigBg(cSaver, study+"-doca", t, ntot, bg_t, prefix+"DOCA", tc.notVetoed(1000.)[0]+additionalCuts, "DOCA of daughter tracks [cm]", (0.64, 0.72, 0.87, 0.89, "Not vetoed"), ("y")) + z = plotSigBg(cSaver, study+"-z", t, ntot, bg_t, prefix+"vtxz", tc.notVetoed(1000.)[0]+additionalCuts, "Reco vertex z position [cm]", (0.165, 0.72, 0.395, 0.89, "Not vetoed"), ("y")) + ndf = plotSigBg(cSaver, study+"-ndf", t, ntot, bg_t, prefix+"MinDaughtersNdf", tc.notVetoed(1000.)[0], "Minimum daughters ndf", (0.50, 0.74, 0.73, 0.92, "Not vetoed"), ()) + mass = plotSigBg(cSaver, study+"-mass", t, ntot, bg_t, prefix+"Mass", tc.notVetoed(1000.)[0]+additionalCuts, "Candidate mass [GeV]", (0.64, 0.72, 0.87, 0.89, "Not vetoed"), ("y")) + mass = plotSigBg(cSaver, study+"-mass-no-bg-mass", t, ntot, bg_t, prefix+"Mass", tc.notVetoed(1000.)[0]+additionalCuts, "Candidate mass [GeV]", (0.64, 0.72, 0.87, 0.89, "Not vetoed"), ("y")) + #mass2 = plotSigBg(cSaver, study+"-mass2", t, ntot, bg_t, prefix+"Mass", tc.notVetoed(1000.)[0]+'&&'+tc.ndf(25)[0]+'&&'+tc.redChi2(5.)[0]+additionalCuts, "Candidate mass [GeV]", (0.57, 0.72, 0.95, 0.89, "Not vetoed, #chi^{2}, ndf"), ("y")) + # Reconstructed + ip_all = plotSigBg(cSaver, study+"-IP0-all", t, ntot, bg_t, prefix+"IP0", tc.recoed()[0]+additionalCuts, "Impact parameter to target [cm]", (0.27, 0.75, 0.58, 0.92, "Reconstructed"), ("y")) + iplin_all = plotSigBg(cSaver, study+"-IP0-lin-all", t, ntot, bg_t, prefix+"IP0", tc.recoed()[0]+additionalCuts, "Impact parameter to target [cm]", (0.64, 0.75, 0.96, 0.92, "Reconstructed"), ()) + chi2_all = plotSigBg(cSaver, study+"-chi2-all", t, ntot, bg_t, prefix+"MaxDaughtersRedChi2", tc.recoed()[0], "Maximum daughters #chi^{2}/ndf", (0.64, 0.72, 0.95, 0.89, "Reconstructed"), ("y")) + doca_all = plotSigBg(cSaver, study+"-doca-all", t, ntot, bg_t, prefix+"DOCA", tc.recoed()[0]+additionalCuts, "DOCA of daughter tracks [cm]", (0.64, 0.72, 0.95, 0.89, "Reconstructed"), ("y")) + z_all = plotSigBg(cSaver, study+"-z-all", t, ntot, bg_t, prefix+"vtxz", tc.recoed()[0]+additionalCuts, "Reco vertex z position [cm]", (0.165, 0.72, 0.475, 0.89, "Reconstructed"), ("y")) + ndf_all = plotSigBg(cSaver, study+"-ndf-all", t, ntot, bg_t, prefix+"MinDaughtersNdf", tc.recoed()[0], "Minimum daughters ndf", (0.43, 0.74, 0.74, 0.92, "Reconstructed"), ()) + mass_all = plotSigBg(cSaver, study+"-mass-all", t, ntot, bg_t, prefix+"Mass", tc.recoed()[0]+additionalCuts, "Candidate mass [GeV]", (0.64, 0.72, 0.95, 0.89, "Reconstructed"), ("y")) + mass_all = plotSigBg(cSaver, study+"-mass-no-bg-mass-all", t, ntot, bg_t, prefix+"Mass", tc.recoed()[0]+additionalCuts, "Candidate mass [GeV]", (0.64, 0.72, 0.95, 0.89, "Reconstructed"), ("y")) + #mass2_all = plotSigBg(cSaver, study+"-mass2-all", t, ntot, bg_t, prefix+"Mass", tc.recoed()[0]+'&&'+tc.ndf(25)[0]+'&&'+tc.redChi2(5.)[0]+additionalCuts, "Candidate mass [GeV]", (0.57, 0.72, 0.98, 0.89, "Reconstructed, #chi^{2}, ndf"), ("y")) + +for name in cSaver.keys(): + cSaver[name][0].Print("accPlotsTP/"+name+".pdf") + cSaver[name][0].Print("accPlotsTP/"+name+".png") + + + +#tex = raw_input("\n\nWould you like me to print some latex acceptance tables?\n") +#if (str(tex) == 'yes'): +res = tableWithDraw() +makeTex(res) + diff --git a/distribsForHNLandBG_byEvent.py b/distribsForHNLandBG_byEvent.py index 8f0ea20..cb8a57a 100755 --- a/distribsForHNLandBG_byEvent.py +++ b/distribsForHNLandBG_byEvent.py @@ -1,9 +1,9 @@ -import sys +import sys,os import ROOT as r import numpy as np import funcsByBarbara as tools -r.gROOT.SetBatch(r.kTRUE) +#r.gROOT.SetBatch(r.kTRUE) r.gROOT.ProcessLine(".X lhcbStyle.C") plots = {} pdg = r.TDatabasePDG.Instance() @@ -80,6 +80,12 @@ pimu = fpimu.Get('t') pimu_geo = tools.searchForNodes3_xyz_dict('../DATA/NewPIMU/geofile_full.10.0.Pythia8-TGeant4-1.0GeV-008.root')#'../DATA/MUPI/geofile_full.10.0.Pythia8-TGeant4-1.0GeV-008-2985556.root') +""" +fpimu = r.TFile('../DATA/signalFromThomas/ShipAna_newGen.root','read')#r.TFile('../DATA/MUPI/ShipAna.root','read') +pimu = fpimu.Get('t') +pimu_geo = tools.searchForNodes3_xyz_dict('../DATA/signalFromThomas/geofile_full.10.0.Pythia8-TGeant4_50k.root')#'../DATA/MUPI/geofile_full.10.0.Pythia8-TGeant4-1.0GeV-008-2985556.root') +""" + fmumunu = r.TFile('../DATA/NewMUMUNU/ShipAna_newGen.root','read')#r.TFile('../DATA/MUMUNU/ShipAna.root','read') mumunu = fmumunu.Get('t') mumunu_geo = tools.searchForNodes3_xyz_dict('../DATA/NewMUMUNU/geofile_full.10.0.Pythia8-TGeant4-1.0GeV-008-mumunu.root')#'../DATA/MUMUNU/geofile_full.10.0.Pythia8-TGeant4-1.0GeV-008-TOTAL.root') @@ -108,9 +114,10 @@ # For a correct normalisation, we subtract the events that have an HNL vertex produced outside of the 60m fiducial volume studies = { - 'pimu': {'data':pimu, 'geo':pimu_geo, 'seen':[], 'ntot':20000-136},#1000}, - 'mumunu': {'data':mumunu, 'geo':mumunu_geo, 'seen':[], 'ntot':20000-289},#7000}, - 'tiny': {'data': tiny, 'geo':tiny_geo, 'seen':[], 'ntot':20000-66}, + 'pimu': {'data':pimu, 'geo':pimu_geo, 'seen':[], 'ntot':20000-125},#136}, + #'pimu': {'data':pimu, 'geo':pimu_geo, 'seen':[], 'ntot':50000-125},#136}, + 'mumunu': {'data':mumunu, 'geo':mumunu_geo, 'seen':[], 'ntot':20000-280},#289}, + 'tiny': {'data': tiny, 'geo':tiny_geo, 'seen':[], 'ntot':20000-43},#66}, 'nu': {'data': nu, 'geo':nu_geo, 'seen':[], 'ntot':1.534e7}, #'cosmics': {'data':cosmics, 'geo':cosmics_geo, 'seen':[], 'ntot':cosmics.GetEntriesFast()}, #'nubar': {'data': nubar, 'geo':nubar_geo, 'seen':[], 'ntot':99999*2},#99999*2}, @@ -126,10 +133,12 @@ # - DOCA < 10 cm # - has no veto? -flux_nu = (1.09e+12 * 2.e+20)/ 5.e+13 -flux_anti = 6.98e+11 * 2.e+20 / 5.e+13 +flux_nu = 1.05e+11 * 2.e+20/ 5.e+13#1.61e+11 * 2.e+20/ 4.e+13 +flux_anti = 7.38e+10 * 2.e+20 / 5.e+13#1.02e+11 * 2.e+20 / 4.e+13 +#flux_nu = (1.09e+12 * 2.e+20)/ 5.e+13 +#flux_anti = 6.98e+11 * 2.e+20 / 5.e+13 entries_nu = 1.534e7#200000-2 -entries_anti = 200000-2 +entries_anti =175*10e+3#entries_anti = 200000-2 weight_nu = "(weight * %s / %s)"%(flux_nu,entries_nu) weight_anti ="(weight * %s / %s)"%(flux_anti,entries_anti) @@ -464,51 +473,115 @@ class cutsWithDraw(object): def __init__(self): - self.noB = '' + self.noB = 'NoB_'#'' self.opt = '' # '_eff' # 'Ethr' + self.nu = False def setNoB(self, noB): # to ignore or take into account Katerina's modifications self.noB = noB + def recoed(self): + cut = " recoed==1 && Sum$(BadTruthVtx==0) " + if self.nu: + cut = " recoed==1 " + string = "Event reconstructed" + return cut, string + def nRecoed(self, num): + cut = " nRecoed==%s "%num + if num == 1: + strn = "1 HNL candidate" + elif num>1: + strn = "%s HNL candidates"%num + string = strn + return cut, string def doca(self, th, noB=None): if not noB: noB = self.noB - return " Sum$(%sDOCA < %s) "%(noB, th) + cut = " Sum$(%sDOCA < %s) "%(noB, th) + if th == int(round(th)): + th = int(round(th)) + string = "DOCA $<$ %s~cm"%th + return cut, string def ip(self, th, noB=None): if not noB: noB = self.noB - return " Sum$(%sIP0 < %s) "%(noB, th) + cut = " Sum$(%sIP0 < %s) "%(noB, th) + string = "IP $<$ %.1f m"%(th/100.) + return cut, string def z(self, zmin, zmax, noB=None): if not noB: noB = self.noB - return " Sum$(%svtxz > %s) && Sum$(%svtxz < %s) "%(noB, zmin, noB, zmax) + cut = " Sum$(%svtxz > %s && %svtxz < %s) "%(noB, zmin, noB, zmax) + string = "$z_{straw\, veto} < z_{vtx} < z_{trk1}$" + return cut, string def ellipse(self, a, b, noB=None): if not noB: noB = self.noB - return " Sum$(( (%svtxx*%svtxx)/(%s*%s) + (%svtxy*%svtxy)/(%s*%s) ) < 1.) "%(noB,noB, a,a, noB,noB, b,b) + cut = " Sum$(( (%svtxx*%svtxx)/(%s*%s) + (%svtxy*%svtxy)/(%s*%s) ) < 1.) "%(noB,noB, a,a, noB,noB, b,b) + string = "Vertex $\\in $ ellipse" + return cut, string + def fiducial(self, zmin, zmax, a, b, noB=None): + cutz = self.z(zmin, zmax, noB)[0] + cute = self.ellipseSq(a, b, noB)[0] + string = "Vtx in fiducial vol." + cut = cutz+"&&"+cute + return cut, string def ellipseSq(self, a, b, noB=None): if not noB: noB = self.noB - return " Sum$(( (%svtxxSqr)/(%s*%s) + (%svtxySqr)/(%s*%s) ) < 1.) "%(noB, a,a, noB, b,b) + cut = " Sum$(( (%svtxxSqr)/(%s*%s) + (%svtxySqr)/(%s*%s) ) < 1.) "%(noB, a,a, noB, b,b) + string = "Vertex $\\in $ ellipse" + return cut, string def converged(self): - return " Sum$(DaughtersFitConverged==1) " + cut = " Sum$(DaughtersFitConverged==1) " + string = "Track fit converged" + return cut, string def goodtracks(self): - return " Sum$(DaughtersAlwaysIn==1) " + cut = " Sum$(DaughtersAlwaysIn==1) " + string = "Tracks in fiducial vol."# $\\in$ ellipse" + return cut, string def redChi2(self, th, noB=None): if not noB: noB = self.noB - return " Sum$(%sMaxDaughtersRedChi2<%s) "%(noB, th) + cut = " Sum$(%sMaxDaughtersRedChi2<%s) "%(noB, th) + if th == int(round(th)): + th = int(round(th)) + string = "$\\chi^2/\\text{N.d.f.} < %s $"%th + return cut, string def ndf(self, th, noB=None): if not noB: noB = self.noB - return " Sum$(%sMinDaughtersNdf>%s) "%(noB, th) + cut = " Sum$(%sMinDaughtersNdf>%s) "%(noB, th) + string = "$\\text{N.d.f.} > %s$"%th + return cut, string def muon1(self, num): - return " Sum$(Has%sMuon1==1) "%num + cut = " Sum$(Has%sMuon1==1) "%num + if num==1: + strn = "1 muon" + elif num==2: + strn = "2 muons" + string = strn+" in 1$^{\\text{st}}$ muon station" + return cut, string def muon2(self, num): - return " Sum$(Has%sMuon2==1) "%num - def notVetoed(self, opt=None): + cut = " Sum$(Has%sMuon2==1) "%num + if num==1: + strn = "1 muon" + elif num==2: + strn = "2 muons" + string = strn+" in 2$^{\\text{nd}}$ muon station" + return cut, string + def ecal(self): + cut = " Sum$(HasEcal==1) " + string = "150~MeV in Ecal" + return cut, string + def notVetoed(self, th, noB=None, opt=None): + if not noB: + noB = self.noB if not opt: opt = self.opt - return " Sum$(numberOfHitLSSegments == 0) && Sum$(strawVetoAny%s==0) && Sum$(upstreamVetoAny%s==0) && Sum$(RPCany%s==0) "%(opt, opt, opt) + cut = " Sum$(%sIP0 < %s) &&"%(noB, th)+" Sum$(numberOfHitLSSegments == 0) && Sum$(strawVetoAny%s==0) && Sum$(upstreamVetoAny%s==0) && Sum$(RPCany%s==0) "%(opt, opt, opt) + string = "Event not vetoed" + return cut, string + def countWithDraw(study): @@ -516,26 +589,177 @@ geo = studies[study]['geo'] ntot = studies[study]['ntot'] tc = cutsWithDraw() - tc.setNoB('NoB_') + if study == 'nu': + tc.setNoB('NoB_') + tc.nu = True zmin = geo['Veto_5']['z']['pos']+geo['Veto_5']['z']['dim'] zmax = geo['Tr1_1']['z']['pos']-geo['Tr1_1']['z']['dim'] - cuts = [ tc.notVetoed(), tc.converged(), tc.redChi2(5.), tc.z(zmin, zmax), tc.ellipseSq(250., 500.), tc.goodtracks(), - tc.muon1(1), tc.muon2(1), tc.doca(30.) ] - cuts.extend([tc.ip(float(i)) for i in range(100, 1000, 50)][::-1]) - cutnames = [ "not vetoed", "converged", "chi2/n<5", "z", "ellipse", "good tracks", "muon1", "muon2", "doca<30" ] - cutnames.extend(["IP<%s"%i for i in range(100, 1000, 50)][::-1]) + if (study == 'pimu'): + cuts = [ tc.recoed()[0], tc.notVetoed(1000.)[0], tc.redChi2(5.)[0], tc.ndf(25)[0], tc.fiducial(zmin, zmax, 250., 500.)[0], tc.goodtracks()[0], + tc.ecal()[0], tc.muon1(1)[0], tc.muon2(1)[0], tc.doca(30.)[0], tc.ip(250.)[0] ] + cutnames = [ tc.recoed()[1], tc.notVetoed(1000.)[1], tc.redChi2(5.)[1], tc.ndf(25)[1], tc.fiducial(zmin, zmax, 250., 500.)[1], tc.goodtracks()[1], + tc.ecal()[1], tc.muon1(1)[1], tc.muon2(1)[1], tc.doca(30.)[1], tc.ip(250.)[1] ] + elif (study == 'nu'): + cuts = [ tc.recoed()[0], tc.notVetoed(1000.)[0], tc.redChi2(5.)[0], tc.ndf(25)[0], tc.fiducial(zmin, zmax, 250., 500.)[0], tc.goodtracks()[0], + tc.ecal()[0], tc.muon1(1)[0], tc.muon2(1)[0], tc.doca(30.)[0], tc.ip(250.)[0] ] + cutnames = [ tc.recoed()[1], tc.notVetoed(1000.)[1], tc.redChi2(5.)[1], tc.ndf(25)[1], tc.fiducial(zmin, zmax, 250., 500.)[1], tc.goodtracks()[1], + tc.ecal()[1], tc.muon1(1)[1], tc.muon2(1)[1], tc.doca(30.)[1], tc.ip(250.)[1] ] + elif (study == 'mumunu'): + cuts = [ tc.recoed()[0], tc.notVetoed(1000.)[0], tc.redChi2(5.)[0], tc.ndf(25)[0], tc.fiducial(zmin, zmax, 250., 500.)[0], tc.goodtracks()[0], + tc.muon1(2)[0], tc.muon2(2)[0], tc.doca(30.)[0], tc.ip(250.)[0] ] + cutnames = [ tc.recoed()[1], tc.notVetoed(1000.)[1], tc.redChi2(5.)[1], tc.ndf(25)[1], tc.fiducial(zmin, zmax, 250., 500.)[1], tc.goodtracks()[1], + tc.muon1(2)[1], tc.muon2(2)[1], tc.doca(30.)[1], tc.ip(250.)[1] ] + elif (study == 'tiny'): + cuts = [ tc.recoed()[0], tc.notVetoed(1000.)[0], tc.redChi2(5.)[0], tc.ndf(25)[0], tc.fiducial(zmin, zmax, 250., 500.)[0], tc.goodtracks()[0], + tc.ecal()[0], tc.doca(30.)[0], tc.ip(250.)[0] ] + cutnames = [ tc.recoed()[1], tc.notVetoed(1000.)[1], tc.redChi2(5.)[1], tc.ndf(25)[1], tc.fiducial(zmin, zmax, 250., 500.)[1], tc.goodtracks()[1], + tc.ecal()[1], tc.doca(30.)[1], tc.ip(250.)[1] ] + # if not cutnames: + # cutnames = [ "converged", "not vetoed", "chi2/n<5", "ndf>15", "z", "ellipse", "good tracks", "ecal", "muon1", "muon2", "doca<30", "IP<250" ] + # cutStrings = [ "Track fit converged", + # "Event not vetoed", + # "$\\xi^2/ndf < 5 $", + # "$ndf > 15$", + # "$z_{straw veto} < z_{vtx} < z_{Trk1}$", + # "Vertex $\\in $ ellipse", + # "Tracks $\\in$ ellipse", + # "150~MeV in Ecal", + # "1$^{\\text{st}}$ muon station", + # "2$^{\\text{nd}}$ muon station", + # "DOCA $<$ 30~cm", + # "IP $<$ 2.5~m", + # ] entries = [] weights = [] print r.gROOT.SetBatch(r.kTRUE) for i,cut in enumerate(cuts): n = t.Draw("event", "&&".join([c for c in cuts[:i+1]])) - t.Draw("event","weight*%s/%s("%(flux_nu, entries_nu) + str("&&".join([c for c in cuts[:i+1]])) + ")") - w = t.GetHistogram().GetSum() / ntot + #if study == 'nu': + # t.Draw("1.>>hw","(weight*%s/%s)*("%(flux_nu, entries_nu) + str("&&".join([c for c in cuts[:i+1]])) + ")") + #elif study == 'nubar': + # t.Draw("1.>>hw","(weight*%s/%s)*("%(flux_anti, entries_anti) + str("&&".join([c for c in cuts[:i+1]])) + ")") + #else: + # t.Draw("1.>>hw","(weight/%s)*("%studies[study]['ntot'] + str("&&".join([c for c in cuts[:i+1]])) + ")") + if study == 'nu': + t.Draw("1.>>hw","(weight*%s/%s)*("%(flux_nu, entries_nu) + str("&&".join([c for c in cuts[:i+1]])) + ")") + elif study == 'nubar': + t.Draw("1.>>hw","(weight*%s/%s)*("%(flux_anti, entries_anti) + str("&&".join([c for c in cuts[:i+1]])) + ")") + else: + t.Draw("1.>>hw","(weight/%s)*("%ntot + str("&&".join([c for c in cuts[:i+1]])) + ")") + w = t.GetHistogram().GetSumOfWeights()#Integral()#Integral('width')#GetSum() + #w = r.gDirectory.Get('hw').Integral() entries.append(n) weights.append(w) - print cutnames[i], '\t\t', n, '\t', w + print cutnames[i], '\t\t\t', n, '\t', w, cuts[i] print + r.gROOT.SetBatch(r.kFALSE) + return entries, weights, cutnames + +def tableWithDraw(): + results = {} + for study in ['nu', 'pimu', 'mumunu', 'tiny']: + results[study] = {} + numbers, numbersWithWeights, cuts = countWithDraw(study) + results[study]['numbers'] = numbers + results[study]['numbersWithWeights'] = numbersWithWeights + results[study]['cuts'] = cuts + return results + + +def plotSigBg(cSaver, name, sig_tree, sig_ntot, bg_tree, what, cuts, xaxis, legSize, log, sig_cuts=None, bg_cuts=None): + if not sig_cuts: + sig_cuts = cuts + if not bg_cuts: + bg_cuts = cuts.replace("&& Sum$(BadTruthVtx==0) ", "") + # colors + bgcolor = r.kRed + sigcolor = r.kBlack + sigfill = r.kGray + bgstyle = 1#9 + r.gStyle.SetOptTitle(r.kFALSE) + # draw normalized + c1 = r.TCanvas(name, name) + whatbg = what + #if what in ["DOCA", "vtxz"]: + # whatbg = "NoB_"+what + if "z" in name: + sig_tree.Draw(what+'>>h_sig_%s(100,-15500.,4500.)'%name, "(weight/%s)*("%sig_ntot + sig_cuts + ")", 'norm') + elif "pimu-IP0-lin" in name: + sig_tree.Draw(what+'>>h_sig_%s(100,0.,20.)'%name, "(weight/%s)*("%sig_ntot + sig_cuts + ")", 'norm') + elif "IP0-lin" in name: + sig_tree.Draw(what+'>>h_sig_%s(100,0.,250.)'%name, "(weight/%s)*("%sig_ntot + sig_cuts + ")", 'norm') + else: + sig_tree.Draw(what+'>>h_sig_%s'%name, "(weight/%s)*("%sig_ntot + sig_cuts + ")", 'norm') + hsig = r.gDirectory.Get('h_sig_%s'%name) + hsig.SetMarkerColor(sigcolor) + hsig.SetLineColor(sigcolor) + hsig.SetFillColor(sigfill) + hsig.GetXaxis().SetTitle(xaxis) + hsig.GetYaxis().SetTitle("a. u.") + hsig.GetYaxis().SetDecimals() + if (("mumunu-mass" in name) or ("mumunu-IP0" in name) or ("tiny-mass" in name) or ("tiny-IP0" in name)) and not ("lin" in name): + hsig.GetYaxis().SetRangeUser(0.0001,5.) + #if "z" in name: + # print name + # hsig.GetXaxis().SetLimits(-17000.,4500.) + # hsig.Draw() + hbg = None + if not "no-bg-mass" in name: + bg_tree.Draw(whatbg+'>>h_bg_%s'%name, "(weight*%s/%s)*("%(flux_nu, entries_nu) + bg_cuts + ")", 'norm same') + hbg = r.gDirectory.Get('h_bg_%s'%name) + hbg.SetMarkerColor(bgcolor) + hbg.SetLineColor(bgcolor) + hbg.SetLineStyle(bgstyle) + if "y" in log: + c1.SetLogy() + if "x" in log: + c1.SetLogx() + newLegSize = list(legSize) + if ("tiny" in name) and not ("2" in name): + newLegSize[2] += 0.02 + newLegSize[3] += 0.02 + leg = r.TLegend(*newLegSize) + leg.SetTextFont(42) + leg.SetLineColor(r.kWhite) + leg.SetBorderSize(0) + sigstring = "HNL #rightarrow #pi#mu" + if ('eenu' in name) or ('tiny' in name): + sigstring = sigstring.replace("#pi#mu", "ee#nu") + if ('mumu' in name): + sigstring = sigstring.replace("#pi#mu", "#mu#mu#nu") + leg.AddEntry(hsig, sigstring, 'lf' ) + if not "no-bg-mass" in name: + leg.AddEntry(hbg, "Neutrino BG", 'l' ) + leg.Draw("same") + c1.Modified(); c1.Update() + print "\n", name, hsig.GetMean(), hsig.GetRMS() + cSaver[name] = (c1, leg, hsig, hbg) + return c1, leg, hsig, hbg + +cSaver = [] +def accPlotsTP(study): + if not os.path.exists('accPlotsTP'): + os.system('mkdir accPlotsTP') + r.gROOT.ProcessLine(".x newGen/mystyle.C") + t = studies[study]['data'] + geo = studies[study]['geo'] + ntot = studies[study]['ntot'] + tc = cutsWithDraw() + zmin = geo['Veto_5']['z']['pos']+geo['Veto_5']['z']['dim'] + zmax = geo['Tr1_1']['z']['pos']-geo['Tr1_1']['z']['dim'] + # BG stuff + bg_t = studies['nu']['data'] + bg_geo = studies['nu']['geo'] + bg_ntot = studies['nu']['ntot'] + bg_tc = cutsWithDraw() + bg_tc.setNoB('NoB_') + bg_zmin = geo['Veto_5']['z']['pos']+geo['Veto_5']['z']['dim'] + bg_zmax = geo['Tr1_1']['z']['pos']-geo['Tr1_1']['z']['dim'] + plot_ip = plotSigBg(study+"-IP0", t, ntot, bg_t, "IP0", tc.notVetoed(1000.)[0], "Impact parameter to target [cm]")[0] + cSaver.append(plot_ip) + raw_input('') + @@ -629,7 +853,40 @@ print return results +header = { 'pimu': ("$HNL \\to \\pi \\mu$", "tab:recpimu"), + 'mumunu': ("$HNL \\to \\mu \\mu \\nu$", "tab:recmmnu"), + 'tiny': ("$HNL \\to e e \\nu$", "tab:receenu"), + 'nu': ("neutrino-induced background", "tab:recnubkg") + } + def makeTex(results): + for study in results.keys(): + weighted_string = "Acceptance" + if study == 'nu': + weighted_string = "Events / 5 years" + results[study]['efficiencies'] = [ results[study]['numbersWithWeights'][i]/results[study]['numbersWithWeights'][i-1] for i in xrange(1, len(results[study]['numbersWithWeights']))] + results[study]['effStrings'] = ['-'] + [ "{:.1f} \\%".format(float(e) * 100) for e in results[study]['efficiencies'] ] + #for s in results[study]['effStrings']: + # print s + print '\n\n\n' + print '\t\\begin{table}[!!!!h]' + print '\t\\begin{center}' + print '\t\t\\caption{ Effect of the offline selection on %s \\label{%s}}'%(header[study][0],header[study][1]) + print '\t\t\\begin{tabular}{|l|c|c|c|}' + print '\t\t\t\\hline' + print '\t\t\t\\textbf{Selection} & \\textbf{Entries} & \\textbf{%s} & \\textbf{Selection efficiency} \\\\'%weighted_string + print '\t\t\t\\hline' + assert(len(results[study]['cuts']) == len(results[study]['numbers']) == len(results[study]['numbersWithWeights'])) + for i in xrange(len(results[study]['cuts'])): + print '\t\t\t', results[study]['cuts'][i], ' & ', '%s'%int(results[study]['numbers'][i]), ' & ', '%s'%latex_float(results[study]['numbersWithWeights'][i]), ' & ', results[study]['effStrings'][i], '\\\\' + print '\t\t\t\\hline' + print '\t\t\\end{tabular}' + print '\t\\end{center}' + print '\t\\end{table}' + print '\n\n' + + +def makeTexOld(results, cuts=None): print print 'Background:' print '\t\\begin{table}' @@ -638,8 +895,9 @@ print '\t\t\t\\hline' print '\t\t\tCut & Weighted $\\nu$ & $\\nu$ entries & Weighted $\\bar{\\nu}$ & $\\bar{\\nu}$ entries & Weighted cosmics & cosmics entries \\\\' print '\t\t\t\\hline' - cuts = ['Track fit converged', '$\chi^2/n < 5$', '$z_{in}+5~\\text{m} < z_{vtx} < z_{out}$', '$z_{vtx} \in $ ellipse', 'Tracks $\in$ ellipse', 'DOCA $<$ 50 cm', - 'IP $<$ 5 m', 'Muon1 \& Muon2', 'DOCA $<$ 30 cm', 'IP $<$ 2.5 m', '$N_{candidates} = 1$'] + if not cuts: + cuts = ['Track fit converged', '$\chi^2/n < 5$', '$z_{in}+5~\\text{m} < z_{vtx} < z_{out}$', '$z_{vtx} \in $ ellipse', 'Tracks $\in$ ellipse', 'DOCA $<$ 50 cm', + 'IP $<$ 5 m', 'Muon1 \& Muon2', 'DOCA $<$ 30 cm', 'IP $<$ 2.5 m', '$N_{candidates} = 1$'] for i in xrange(len(results['nu']['numbers'])): print '\t\t\t', cuts[i], ' & ', '%s'%latex_float(results['nu']['numbersWithWeights'][i]), ' & ', '%s'%results['nu']['numbers'][i], ' & ',\ '%s'%latex_float(results['nubar']['numbersWithWeights'][i]), ' & ', '%s'%results['nubar']['numbers'][i], ' & ',\ @@ -655,8 +913,8 @@ print '\t\t\t\\hline' print '\t\t\tCut & $\\mathcal{A}(\\pi\\mu)$ & $\\pi\\mu$ entries & $\\mathcal{A}(\\mu\\mu\\nu)$ & $\\mu\\mu\\nu$ entries & $\\mathcal{A}(ee\\nu)$ & $ee\\nu$ entries \\\\' print '\t\t\t\\hline' - cuts = ['Track fit converged', '$\chi^2/n < 5$', '$z_{in}+5~\\text{m} < z_{vtx} < z_{out}$', '$z_{vtx} \in $ ellipse', 'Tracks $\in$ ellipse', 'DOCA $<$ 50 cm', - 'IP $<$ 5 m', 'Muon1 \& Muon2', 'DOCA $<$ 30 cm', 'IP $<$ 2.5 m', '$N_{candidates} = 1$'] + #cuts = ['Track fit converged', '$\chi^2/n < 5$', '$z_{in}+5~\\text{m} < z_{vtx} < z_{out}$', '$z_{vtx} \in $ ellipse', 'Tracks $\in$ ellipse', 'DOCA $<$ 50 cm', + # 'IP $<$ 5 m', 'Muon1 \& Muon2', 'DOCA $<$ 30 cm', 'IP $<$ 2.5 m', '$N_{candidates} = 1$'] for i in xrange(len(results['nu']['numbers'])): print '\t\t\t', cuts[i], ' & ', '%s'%latex_float(results['pimu']['numbersWithWeights'][i]), ' & ', '%s'%results['pimu']['numbers'][i], ' & ',\ '%s'%latex_float(results['mumunu']['numbersWithWeights'][i]), ' & ', '%s'%results['mumunu']['numbers'][i], ' & ',\ diff --git a/newGen/mystyle.C b/newGen/mystyle.C new file mode 100644 index 0000000..a903c2f --- /dev/null +++ b/newGen/mystyle.C @@ -0,0 +1,179 @@ +{ +//////////////////////////////////////////////////////////////////// +// PURPOSE: +// +// This macro defines a reasonable style for (black-and-white) +// "publication quality" ROOT plots. The default settings contain +// many features that are either not desirable for printing on white +// paper or impair the general readibility of plots. +// +// USAGE: +// +// Simply include the line +// gROOT->ProcessLine(".x $LHCBSTYLE/root/lhcbstyle.C"); +// at the beginning of your root macro. +// +// SOME COMMENTS: +// +// Statistics and fit boxes: +// +// "Decorative" items around the histogram are kept to a minimum. +// In particular there is no box with statistics or fit information. +// You can easily change this either by editing your private copy +// of this style file or by calls to "gStyle" in your macro. +// For example, +// gStyle->SetOptFit(1011); +// will add some fit information. +// +// Font: +// +// The font is chosen to be 62, i.e.helvetica-bold-r-normal with +// precision 2. Font is of course a matter of taste, but most people +// will probably agree that Helvetica bold gives close to optimal +// readibility in presentations. It appears to be the ROOT default, +// and since there are still some features in ROOT that simply won't +// respond to any font requests, it is the wise choice to avoid +// ugly font mixtures on the same plot... The precision of the font (2) +// is chosen in order to have a rotatable and scalable font. Be sure +// to use true-type fonts! I.e. +// Unix.*.Root.UseTTFonts: true in your .rootrc file. +// +// "Landscape histograms": +// +// The style here is designed for more or less quadratic plots. +// For very long histograms, adjustements are needed. For instance, +// for a canvas with 1x5 histograms: +// TCanvas* c1 = new TCanvas("c1", "L0 muons", 600, 800); +// c1->Divide(1,5); +// adaptions like the following will be needed: +// gStyle->SetTickLength(0.05,"x"); +// gStyle->SetTickLength(0.01,"y"); +// gStyle->SetLabelSize(0.15,"x"); +// gStyle->SetLabelSize(0.1,"y"); +// gStyle->SetStatW(0.15); +// gStyle->SetStatH(0.5); +// +//////////////////////////////////////////////////////////////////// + +gROOT->Reset(); + +cout << "executing lhcbStyle.C:" << endl; +cout << " " << endl; +cout << " " << endl; +cout << " $ $ $ $$$ $ " << endl; +cout << " $ $ $ $ $ " << endl; +cout << " $ $$$$$ $ $$$ " << endl; +cout << " $ $ $ $ $ $ " << endl; +cout << " $$$$$ $ $ $$$ $$$ " << endl; +cout << " " << endl; +cout << " LHCb ROOT style file " << endl; +cout << " " << endl; +cout << +" Problems, suggestions, contributions to Thomas.Schietinger@cern.ch" + << endl; +cout << " " << endl; + +TStyle *lhcbStyle= new TStyle("lhcbStyle","LHCb official plots style"); + +// use helvetica-bold-r-normal, precision 2 (rotatable) + Int_t lhcbFont = 42; //62; +// line thickness +Double_t lhcbWidth = 2; + +// use plain black on white colors +lhcbStyle->SetFrameBorderMode(0); +lhcbStyle->SetCanvasBorderMode(0); +lhcbStyle->SetPadBorderMode(0); +lhcbStyle->SetPadColor(0); +lhcbStyle->SetCanvasColor(0); +lhcbStyle->SetStatColor(0); +lhcbStyle->SetPalette(1); +lhcbStyle->SetTitleColor(1); + lhcbStyle->SetFillColor(1);//0 + lhcbStyle->SetFillStyle(0); + +// set the paper & margin sizes +lhcbStyle->SetPaperSize(20,26); +lhcbStyle->SetPadTopMargin(0.05); +lhcbStyle->SetPadRightMargin(0.05); // increase for colz plots!! +lhcbStyle->SetPadBottomMargin(0.16); +lhcbStyle->SetPadLeftMargin(0.14); + +// use large fonts +lhcbStyle->SetTextFont(lhcbFont); +lhcbStyle->SetTextSize(0.08); +lhcbStyle->SetLabelFont(lhcbFont,"x"); +lhcbStyle->SetLabelFont(lhcbFont,"y"); +lhcbStyle->SetLabelFont(lhcbFont,"z"); +lhcbStyle->SetLabelSize(0.05,"x"); +lhcbStyle->SetLabelSize(0.05,"y"); +lhcbStyle->SetLabelSize(0.05,"z"); +lhcbStyle->SetTitleFont(lhcbFont); +lhcbStyle->SetTitleFont(lhcbFont,"y"); +lhcbStyle->SetTitleFont(lhcbFont,"x"); +lhcbStyle->SetTitleFont(lhcbFont,"z"); +lhcbStyle->SetTitleSize(0.06,"x"); +lhcbStyle->SetTitleSize(0.06,"y"); +lhcbStyle->SetTitleSize(0.06,"z"); + +// use bold lines and markers +lhcbStyle->SetLineWidth(lhcbWidth); +lhcbStyle->SetFrameLineWidth(lhcbWidth); +lhcbStyle->SetHistLineWidth(lhcbWidth); +lhcbStyle->SetFuncWidth(lhcbWidth); +lhcbStyle->SetGridWidth(lhcbWidth); +lhcbStyle->SetLineStyleString(2,"[12 12]"); // postscript dashes +lhcbStyle->SetMarkerStyle(8); +lhcbStyle->SetMarkerSize(1.5); + +// label offsets +lhcbStyle->SetLabelOffset(0.015); + +// by default, do not display histogram decorations: +lhcbStyle->SetOptStat(0); +//lhcbStyle->SetOptStat(1110); // show only nent, mean, rms +lhcbStyle->SetOptTitle(0); +lhcbStyle->SetOptFit(0); +//lhcbStyle->SetOptFit(1011); // show probability, parameters and errors + +// look of the statistics box: +lhcbStyle->SetStatBorderSize(1); +lhcbStyle->SetStatFont(lhcbFont); +lhcbStyle->SetStatFontSize(0.05); +lhcbStyle->SetStatX(0.9); +lhcbStyle->SetStatY(0.9); +lhcbStyle->SetStatW(0.25); +lhcbStyle->SetStatH(0.15); + +// put tick marks on top and RHS of plots +lhcbStyle->SetPadTickX(1); +lhcbStyle->SetPadTickY(1); + +// histogram divisions: only 5 in x to avoid label overlaps +lhcbStyle->SetNdivisions(505,"x"); +lhcbStyle->SetNdivisions(510,"y"); + +gROOT->SetStyle("lhcbStyle"); +gROOT->ForceStyle(); + +TPaveText *lhcbName = new TPaveText(0.65,0.8,0.9,0.9,"BRNDC"); +lhcbName->SetFillColor(0); +lhcbName->SetTextAlign(12); +lhcbName->SetBorderSize(0); +lhcbName->AddText("LHCb"); + +TText *lhcbLabel = new TText(); +lhcbLabel->SetTextFont(lhcbFont); +lhcbLabel->SetTextColor(1); +lhcbLabel->SetTextSize(0.04); +lhcbLabel->SetTextAlign(12); + +TLatex *lhcbLatex = new TLatex(); +lhcbLatex->SetTextFont(lhcbFont); +lhcbLatex->SetTextColor(1); +lhcbLatex->SetTextSize(0.04); +lhcbLatex->SetTextAlign(12); +lhcbStyle->SetLegendFillColor(kWhite); + +} + diff --git a/newGen/offlineForBarbara.py b/newGen/offlineForBarbara.py index 9fe706e..954d478 100644 --- a/newGen/offlineForBarbara.py +++ b/newGen/offlineForBarbara.py @@ -472,6 +472,22 @@ # By particle global NoB_vtxx, NoB_vtxy, NoB_vtxz global NoB_IP0, NoB_DOCA + HNLPos = ROOT.TLorentzVector() + particle.ProductionVertex(HNLPos) + xv, yv, zv, doca = HNLPos.X(),HNLPos.Y(),HNLPos.Z(),HNLPos.T() + tools.Push(NoB_DOCA, doca) + tools.Push(NoB_vtxx, xv*xv); tools.Push(NoB_vtxy, yv*yv); tools.Push(NoB_vtxz, zv) + # impact parameter to target + HNLMom = ROOT.TLorentzVector() + particle.Momentum(HNLMom) + tr = ROOT.TVector3(0,0,ShipGeo.target.z0) + t = 0 + for i in range(3): t += HNLMom(i)/HNLMom.P()*(tr(i)-HNLPos(i)) + ip = 0 + for i in range(3): ip += (tr(i)-HNLPos(i)-t*HNLMom(i)/HNLMom.P())**2 + ip = ROOT.TMath.Sqrt(ip) + tools.Push(NoB_IP0, ip) + """ t1,t2 = particle.GetDaughter(0),particle.GetDaughter(1) PosDir = {} for tr in [t1,t2]: @@ -493,6 +509,7 @@ for i in range(3): ip += (tr(i)-HNLPos(i)-t*HNLMom(i)/HNLMom.P())**2 ip = ROOT.TMath.Sqrt(ip) tools.Push(NoB_IP0, ip) + """ def kinematics(tree, particle, ntr, nref): diff --git a/pleaseRunMe/elenaOutput.txt b/pleaseRunMe/elenaOutput.txt new file mode 100644 index 0000000..4557603 --- /dev/null +++ b/pleaseRunMe/elenaOutput.txt @@ -0,0 +1,45 @@ +[INFO ] Media file used : /home/ubuntu/ShipSoft/FairShip/geometry/media.geo +Refitting 2 tracks in event 1 +genfit failed to fit track +genfit failed to fit track +Refitting 2 tracks in event 4 +genfit failed to fit track +genfit failed to fit track +Refitting 2 tracks in event 8 +genfit failed to fit track +genfit failed to fit track +Refitting 3 tracks in event 15 +genfit failed to fit track +genfit failed to fit track +genfit failed to fit track +Refitting 3 tracks in event 17 +genfit failed to fit track +genfit failed to fit track +genfit failed to fit track +Refitting 2 tracks in event 22 +genfit failed to fit track +genfit failed to fit track +Refitting 2 tracks in event 28 +genfit failed to fit track +genfit failed to fit track +Refitting 2 tracks in event 31 +genfit failed to fit track +genfit failed to fit track +Refitting 2 tracks in event 33 +genfit failed to fit track +genfit failed to fit track +Refitting 2 tracks in event 38 +genfit failed to fit track +genfit failed to fit track +Refitting 2 tracks in event 42 +genfit failed to fit track +genfit failed to fit track +1 +genfit failed to fit track +0 +2 +genfit failed to fit track +genfit failed to fit track +0 +0 +0 diff --git a/plotsForKaterina.py b/plotsForKaterina.py new file mode 100644 index 0000000..912c1ba --- /dev/null +++ b/plotsForKaterina.py @@ -0,0 +1,44 @@ +import os +import ROOT as r + +if not os.path.exists('PlotsForKaterina'): + os.system('mkdir PlotsForKaterina') + +cSaver = {} +f, t = {}, {} + +r.gStyle.SetOptStat('') +r.gStyle.SetOptTitle(r.kFALSE) + +for data in ['PIMU', 'MUMUNU']: + cSaver[data] = [] + f[data] = r.TFile("../DATA/New%s/ShipAna_newGen.root"%data, "read") + t[data] = f[data].Get('t') + c1 = r.TCanvas(data+'-'+"vtxz",data+'-'+"vtxz") + t[data].Draw('vtxz>>histo') + r.gDirectory.Get('histo').SetMarkerColor(r.kRed) + r.gDirectory.Get('histo').SetLineColor(r.kRed) + t[data].Draw('NoB_vtxz', '', 'same') + c1.BuildLegend(0.3, 0.67, 0.58, 0.88) + cSaver[data].append(c1) + c2 = r.TCanvas(data+'-'+"doca",data+'-'+"doca") + t[data].Draw('DOCA>>histo2') + r.gDirectory.Get('histo2').SetMarkerColor(r.kRed) + r.gDirectory.Get('histo2').SetLineColor(r.kRed) + t[data].Draw('NoB_DOCA','','same') + c2.BuildLegend() + c2.SetLogy() + cSaver[data].append(c2) + c3 = r.TCanvas(data+'-'+"vtxxy",data+'-'+"vtxxy") + t[data].SetMarkerSize(0.3) + t[data].SetMarkerColor(r.kRed) + t[data].SetLineColor(r.kRed) + t[data].Draw('(TMath::Sqrt(vtxySqr)) : (TMath::Sqrt(vtxxSqr)) >> histo3') + t[data].SetMarkerColor(r.kBlack) + t[data].SetLineColor(r.kBlack) + t[data].Draw('(TMath::Sqrt(NoB_vtxySqr)) : (TMath::Sqrt(NoB_vtxxSqr))','','same') + cSaver[data].append(c3) + + for canvas in cSaver[data]: + canvas.Modified(); canvas.Update() + canvas.Print('PlotsForKaterina/'+canvas.GetName()+'.pdf') \ No newline at end of file