Context Navigation

resolutions.C@ 733

Last change on this file since 733 was 503, checked in by Xavier Rouby, 15 years ago
added CMS expectation band for etmis
File size: 14.7 KB

Rev	Line
[443]	1	/***********************************************************************
	2
	3	/----------------------------------------------\
	4	\| Delphes, a framework for the fast simulation \|
	5	\| of a generic collider experiment \|
	6	\------------- arXiv:0903.2225v1 ------------/
	7
	8
	9	This package uses:
	10	------------------
	11	ROOT: Nucl. Inst. & Meth. in Phys. Res. A389 (1997) 81-86
	12	FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210]
	13	Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2]
	14	FROG: [hep-ex/0901.2718v1]
	15	HepMC: Comput. Phys. Commun.134 (2001) 41
	16
	17	------------------------------------------------------------------
	18
	19	Main authors:
	20	-------------
	21
	22	Severine Ovyn Xavier Rouby
	23	severine.ovyn@uclouvain.be xavier.rouby@cern
	24
	25	Center for Particle Physics and Phenomenology (CP3)
	26	Universite catholique de Louvain (UCL)
	27	Louvain-la-Neuve, Belgium
	28
	29	Copyright (C) 2008-2009,
	30	All rights reserved.
	31
	32	***********************************************************************/
	33
[9]	34	#include "TROOT.h"
	35	#include "TFile.h"
	36	#include "TTree.h"
	37	#include "TCanvas.h"
	38	#include "TProfile.h"
	39	#include "TF1.h"
	40	#include "TGraph.h"
[503]	41	#include "TMultiGraph.h"
	42	#include "TGraphErrors.h"
[153]	43	#include "TLegend.h"
[9]	44
	45	#include "interface/FuncDef.h"
	46
[27]	47	void JetResol()
[9]	48	{
[39]	49	setTDRStyle();
	50	gROOT->Reset();
	51
[258]	52	TFile *f1 = new TFile("JET2.root","read");
[473]	53	if(!f1->IsOpen()) { cout << "could not open "<< f1->GetName() << ". Exiting..." << endl; return;}
	54	if(!f1->FindKey("Analysis")) { cout << "Bad input file, could not find the \"Analysis\" tree. Exiting..." << endl; return;}
[39]	55	TTree Analyze = (TTree)f1->Get("Analysis");
	56
[153]	57	const Int_t numBin=16;
	58	double bins[numBin]={0,10,20,30,40,50,60,70,80,100,120,140,180,220,300,1200};
[39]	59	TProfile *ETSoverET = new TProfile("ETSoverET","Jet resolution: E_{T}^{rec}/E_{T}^{mc}",(numBin-1),bins,-10,10);
	60
	61	double rms[numBin];
	62	double mean[numBin];
	63
	64	TCanvas *c1 = new TCanvas("c1","JET resol",0,0,1000,650);
	65	c1->cd(); int frame=0;
	66	c1->Divide(6,2);
	67
	68	float x[numBin-1];
	69	float y[numBin-1];
	70	float ex[numBin-1];
	71	float ey[numBin-1];
[503]	72	float erec_over_etruth[numBin-1];
[39]	73
	74	float finval=0;//valeur a remplir puis a fitter
	75
	76	for ( int i=0; i<(numBin-1); i++) // premiÃšre bin : i ==1 et pas i == 0
	77	{
	78	TAxis *xaxis = ETSoverET->GetXaxis();
	79	float binCenter = xaxis->GetBinCenter(i+1);
	80	int binMin = (int)xaxis->GetBinLowEdge(i+1);
	81	int binMax = (int)xaxis->GetBinUpEdge(i+1);
	82	char tempMin[500];
	83	if(i==0)binMin=5;
[91]	84	sprintf(tempMin,"JetPTResol.PT > %d",binMin);
[39]	85	string mystringMin(tempMin);
	86	char tempMax[500];
[91]	87	sprintf(tempMax,"JetPTResol.PT < %d",binMax);
[39]	88	string mystringMax(tempMax);
	89	char tempName[500];
[91]	90	sprintf(tempName,"(JetPTResol.SmearedPT)>>hETSoverET%d",i);
[39]	91	string mystringName(tempName);
	92
	93	c1->cd(++frame);
	94	GaussValues(Analyze,tempName,rms[i],mean[i],mystringMin,mystringMax);
	95	x[i]=binCenter;
	96	finval=rms[i]/mean[i];
	97	y[i]=(finval*100);
	98	ex[i]=0;
	99	ey[i]=0;
[503]	100	erec_over_etruth[i]=mean[i]+1;
[39]	101	}
[503]	102
	103
	104	TCanvas *c3 = new TCanvas("c3","JET linearity",100,100,600,450);
	105	c3->cd();
	106	//for(int i=0; i<numBin-1; i++) cout << x[i] << "\t\t" << erec_over_etruth[i] << endl;
	107	TGraph *linearity = new TGraph(numBin-1,x,erec_over_etruth);
	108	linearity->Draw("AP");
	109	linearity->SetTitle("");
	110	linearity->GetXaxis()->SetTitle("E^{Truth} [GeV]");
	111	linearity->GetYaxis()->SetTitle("<E^{Reco}/E^{Truth}>");
	112	linearity->GetXaxis()->SetRangeUser(30,1000);
	113	linearity->SetMinimum(0.85);
	114	linearity->SetMaximum(1.20);
	115	gPad->SetLogx();
	116
	117
[39]	118	TCanvas *c2 = new TCanvas("c2","JET resol",100,100,600,450);
	119	c2->cd();
	120
[153]	121	TF1 *fitfun = new TF1("user","sqrt(pow([0]/x,2)+pow([1]/sqrt(x),2)+pow([2],2))",10,800);
	122	TF1 fitfunCMS = new TF1("userCMS","sqrt(pow(1.05100/sqrt(x),2)+pow(0.06100,2)+pow(0100/x,2))",10,800);
	123	//TF1 fitfunCMS = new TF1("userCMS","sqrt(pow(1.17100/sqrt(x),2)+pow(0.039*100,2))",10,800);
[39]	124
[153]	125	TGraph *gr11 = new TGraph((numBin-1),x,y);
[39]	126	gr11->Draw("AP");
	127	gr11->SetTitle("");
	128	gr11->GetXaxis()->SetTitle("E_{T}^{MC} [GeV]");
	129	gr11->GetYaxis()->SetRangeUser(0,50);
	130	gr11->GetYaxis()->SetTitle("#sigma(E_{T}^{rec}/E_{T}^{MC})_{fit}/<E_{T}^{rec}/E_{T}^{MC}>_{fit}");
	131
	132	Double_t* params = fitfun->GetParameters();
[153]	133
	134	fitfun->SetLineColor(6);
[39]	135	gr11->Fit("user","QR");
[153]	136	gr11->Fit("user","QR");
	137	gr11->Fit("user","QR");
	138	gr11->Fit("user","QR");
	139	char tempResol[100];
	140	sprintf(tempResol,"Delphes resolution: #frac{#sigma(E_{T}^{rec}/E_{T}^{MC})}{<E_{T}^{rec}/E_{T}^{MC}>} = #frac{%f}{#sqrt{E_{T}^{MC}}} #oplus %f",params[1],params[2]);
	141	char tempResol2[100];
	142	sprintf(tempResol2,"sqrt(pow(%f/sqrt(x),2)+pow(%f,2))",params[1],params[2]);
[27]	143
[153]	144
	145	TF1 *fitfunDelphes = new TF1("userDelphes",tempResol2,7,1000);
	146	fitfunDelphes->SetLineColor(6);
	147	fitfunDelphes->Draw("same");
	148
	149	//CMS values
	150	fitfunCMS->SetLineColor(9);
	151	fitfunCMS->Draw("same");
	152
[60]	153	TPaveText *events = MakeTPave(0.2,0.75,0.35,0.8,"Events: pp #rightarrow gg ");
	154	events->Draw();
	155
[39]	156	TPaveText *Delphes = MakeTPave(0.2,0.15,0.35,0.2,"MG/ME + Delphes");
	157	Delphes->Draw();
[153]	158
	159	TLegend *legend = new TLegend(0.2,0.6,0.9,0.85,NULL,"NDC");
	160	legend->AddEntry(fitfunCMS,"CMS resolution","l");
	161	legend->AddEntry(fitfunDelphes,tempResol,"l");
	162	legend->SetFillColor(10);
	163	legend->SetBorderSize(0);
	164	legend->Draw();
	165
[39]	166	delete fitfun;
	167	}
[27]	168
[39]	169	void ElecResol()
	170	{
[9]	171
[39]	172	setTDRStyle();
	173	gROOT->Reset();
	174
[153]	175	TFile *f1 = new TFile("ETMIS2.root","read");
[473]	176	if(!f1->IsOpen()) { cout << "could not open "<< f1->GetName() << ". Exiting..." << endl; return;}
	177	if(!f1->FindKey("Analysis")) { cout << "Bad input file, could not find the \"Analysis\" tree. Exiting..." << endl; return;}
[39]	178	TTree Analyze = (TTree)f1->Get("Analysis");
	179
[153]	180	const Int_t numBin=11;
	181	double bins[numBin]={0,10,20,30,40,50,60,70,80,100,120};
[39]	182	TProfile *ETSminusET = new TProfile("ETSminusET","Electron resolution: E_{T}^{rec}/E_{T}^{mc}",(numBin-1),bins,-10,10);
	183
	184	double rms[numBin];
	185	double mean[numBin];
	186
	187	TCanvas *c3 = new TCanvas("c3","ELEC resol",0,0,1000,650);
	188	c3->cd(); int frame=0;
	189	c3->Divide(6,2);
	190
	191	float x[numBin-1];
	192	float y[numBin-1];
	193	float ex[numBin-1];
	194	float ey[numBin-1];
	195
	196	float finval=0;//valeur a remplir puis a fitter
	197
	198	for ( int i=0; i<(numBin-1); i++) // premiÃšre bin : i ==1 et pas i == 0
	199	{
	200	TAxis *xaxis = ETSminusET->GetXaxis();
	201	float binCenter = xaxis->GetBinCenter(i+1);
	202	int binMin = (int)xaxis->GetBinLowEdge(i+1);
	203	int binMax = (int)xaxis->GetBinUpEdge(i+1);
	204	char tempMin[500];
	205	if(i==0)binMin=5;
[89]	206	sprintf(tempMin,"ElecEResol.E > %d",binMin);
[39]	207	string mystringMin(tempMin);
	208	char tempMax[500];
[89]	209	sprintf(tempMax,"ElecEResol.E < %d",binMax);
[39]	210	string mystringMax(tempMax);
	211	char tempName[500];
[89]	212	sprintf(tempName,"(ElecEResol.E-ElecEResol.SmearedE)>>hETSoverET%d",i);
[39]	213	string mystringName(tempName);
	214
	215	c3->cd(++frame);
[153]	216	GaussValuesElec(Analyze,tempName,rms[i],mean[i],mystringMin,mystringMax);
	217	//GaussValues(Analyze,tempName,rms[i],mean[i],mystringMin,mystringMax);
[39]	218	x[i]=binCenter;
	219	finval=rms[i]/binCenter;
	220	y[i]=(finval*100);
	221	ex[i]=0;
	222	ey[i]=0;
	223	}
	224
	225	TCanvas *c4 = new TCanvas("c4","ELEC resol",100,100,600,450);
	226	c4->cd();
	227
	228	TF1 *fitfun = new TF1("user","sqrt(pow([0],2)+pow([1]/sqrt(x),2)+pow([2]/x,2))",1,400);
	229
	230	TGraphErrors *gr11 = new TGraphErrors((numBin-1),x,y,ex,ey);
	231	gr11->Draw("AP");
	232	gr11->SetTitle("");
	233	gr11->GetXaxis()->SetTitle("E [GeV]");
	234	gr11->GetYaxis()->SetRangeUser(0,5);
	235	gr11->GetYaxis()->SetTitle("#sigma/E");
	236
	237	Double_t* params = fitfun->GetParameters();
	238
	239	gr11->Fit("user","QR");
	240	gr11->Fit("user","QRI");
	241	gr11->Fit("user","QRI");
	242	gr11->Fit("user","QRI");
	243	char tempResol[500];
	244	sprintf(tempResol,"#frac{#sigma}{E} = #frac{%f}{#sqrt{E}} #oplus #frac{%f}{E} #oplus %f",params[1]/100,params[2]/100,params[0]/100);
	245
	246	TPaveText *leg1 = MakeTPave(0.4,0.6,0.8,0.65,tempResol);
	247	leg1->Draw();
	248
	249	TPaveText *Delphes = MakeTPave(0.2,0.15,0.35,0.2,"MG/ME + Delphes");
	250	Delphes->Draw();
	251
[60]	252	TPaveText *events = MakeTPave(0.2,0.75,0.35,0.8,"Events: WHq'#rightarrow W#tau#tauq'#rightarrowjjl#tauq', m_{H}=150 GeV ");
[39]	253	events->Draw();
	254
	255	delete fitfun;
	256	}
[9]	257
[39]	258	void TauJetInfo()
	259	{
	260
	261	setTDRStyle();
	262	gROOT->Reset();
	263
[258]	264	TFile *f1 = new TFile("TAUJET2.root","read");
[473]	265	if(!f1->IsOpen()) { cout << "could not open "<< f1->GetName() << ". Exiting..." << endl; return;}
	266	if(!f1->FindKey("Analysis")) { cout << "Bad input file, could not find the \"Analysis\" tree. Exiting..." << endl; return;}
[39]	267	TTree Analyze = (TTree)f1->Get("Analysis");
	268
[60]	269	TCanvas *ct1 = new TCanvas("ct1","Tau information",100,100,600,450);
	270	ct1->cd();
[39]	271
	272	TH1F *tauEnergy =MakeNormTH1F(20,0.8,1,Analyze,"TauJetPTResol.EnergieCen>>tauEnergy",1, 0, 1,2,false);
	273	tauEnergy->Draw();
	274	tauEnergy->SetTitle("");
	275	tauEnergy->GetYaxis()->SetTitle("Fraction of events");
	276	tauEnergy->GetXaxis()->SetTitle("C_{#tau}");
	277
	278
	279	TPaveText *Delphes1 = MakeTPave(0.3,0.85,0.45,0.9,"MG/ME + Delphes");
	280	Delphes1->Draw();
	281	TPaveText *ctau = MakeTPave(0.3,0.75,0.45,0.8,"C_{#tau} = #frac{#sum E^{towers} (#DeltaR = 0.15)}{E^{jet}}");
	282	ctau->Draw();
[153]	283	TPaveText *events1 = MakeTPave(0.3,0.65,0.45,0.7,"Events: WHq'#rightarrow WWWq'#rightarrow lllq', m_{H}=150 GeV ");
[39]	284	events1->Draw();
	285
	286
[60]	287	TCanvas *ct2 = new TCanvas("ct2","Tau information",100,100,600,450);
	288	ct2->cd();
[39]	289	TH1F *NumTrack =MakeNormTH1F(6,0,6,Analyze,"TauJetPTResol.NumTrack>>NumTrack",1, 0, 1,2,false);
	290	NumTrack->Draw();
	291	NumTrack->SetTitle("");
	292	NumTrack->GetYaxis()->SetTitle("Fraction of events");
	293	NumTrack->GetXaxis()->SetTitle("N^{tracks}");
	294
	295	TPaveText *Delphes = MakeTPave(0.6,0.85,0.85,0.9,"MG/ME + Delphes");
	296	Delphes->Draw();
	297	TPaveText *numtracks = MakeTPave(0.6,0.75,0.85,0.8,"#DeltaR < 0.4, p_{T}^{track} > 2 GeV");
	298	numtracks->Draw();
[153]	299	TPaveText *events = MakeTPave(0.6,0.65,0.85,0.7,"Events: WHq'#rightarrow WWWq'#rightarrow lllq', m_{H}=150 GeV ");
[39]	300	events->Draw();
	301
[9]	302	}
[23]	303
[27]	304	void ETmisResol()
[23]	305	{
[39]	306
	307	setTDRStyle();
	308	gROOT->Reset();
	309
[473]	310	//TFile *f1 = new TFile("JET2.root","read");
	311	TFile *f1 = new TFile("JET2_cms_midpoint.root","read");
	312
	313	if(!f1->IsOpen()) { cout << "could not open "<< f1->GetName() << ". Exiting..." << endl; return;}
	314	if(!f1->FindKey("Analysis")) { cout << "Bad input file, could not find the \"Analysis\" tree. Exiting..." << endl; return;}
[39]	315	TTree Analyze = (TTree)f1->Get("Analysis");
	316
	317	TF1 fitfun = new TF1("user","[0]sqrt(x)",0,600);
[503]	318	fitfun->SetLineColor(kBlack);
	319	fitfun->SetLineStyle(2);
[153]	320	const Int_t numBin=6;
	321	double bins[numBin]={180,260,340,420,500,580};
[39]	322	TProfile *ETSoverET = new TProfile("ETSoverET","ETmis resol",(numBin-1),bins,-1000,1000);
	323
	324
	325	double rms[numBin-1];
	326
	327	TCanvas *c5 = new TCanvas("c5","PTmis resol",0,0,1000,650);
	328	c5->cd(); int frame=0;
	329	c5->Divide(6,2);
	330
	331	double x[numBin];
	332	double y[numBin];
	333
	334	for ( int i=0; i<(numBin-1); i++) // premiÃšre bin : i ==1 et pas i == 0
	335	{
	336	TAxis *xaxis = ETSoverET->GetXaxis();
	337	float binCenter = xaxis->GetBinCenter(i+1);
	338	int binMin = (int)xaxis->GetBinLowEdge(i+1);
	339	int binMax = (int)xaxis->GetBinUpEdge(i+1);
	340	char tempMin[500];
	341	sprintf(tempMin,"ETmisResol.SEt>%d",binMin);
	342	string mystringMin(tempMin);
	343	char tempMax[500];
	344	sprintf(tempMax,"ETmisResol.SEt<%d",binMax);
	345	string mystringMax(tempMax);
	346
	347	char tempName[500];
	348	sprintf(tempName,"ETmisResol.ExSmeare>>hETSoverET%d",i);
	349	string mystringName(tempName);
	350
	351	c5->cd(++frame);
	352	GaussValuesETmis(Analyze,tempName,rms[i],mystringMin,mystringMax);
	353	x[i]=binCenter;
	354	y[i]=rms[i];
	355
	356	}
	357
	358	TCanvas *c6 = new TCanvas("c6","ETmis resolution",100,100,600,450);
	359	c6->cd();
[503]	360
	361	TF1 fdo = new TF1("fdo","0.6sqrt(x)",1,600);
	362	TF1 fup = new TF1("fup","0.7sqrt(x)",1,600);
	363	int N = 100;
	364	float xx[N], yy[N], eex[N], eey[N];
	365	for (int i=0; i<N; i++) {
	366	xx[i] = 6*i;
	367	yy[i] = (fup->Eval(xx[i]) + fdo->Eval(xx[i]))/2.;
	368	eex[i] = 0;
	369	eey[i] = (fup->Eval(xx[i]) - fdo->Eval(xx[i]))/2.;
	370	}
	371	TGraphErrors *gr = new TGraphErrors(N,xx,yy,eex,eey);
	372	gr->SetFillColor(18);
	373
[39]	374
	375	x[numBin]=0;
	376	y[numBin]=0;
[473]	377	TGraph *gr11 = new TGraph(numBin,x,rms);
	378	gr11->SetTitle("");
	379	gr11->SetMarkerStyle(3);
[503]	380	gr11->SetLineStyle(2);
[39]	381	gr11->Fit("user","RQ");
	382	gr11->Fit("user","RQ");
	383	gr11->Fit("user","RQ");
	384	gr11->Fit("user","RQ");
[23]	385
[503]	386
	387	TMultiGraph *mgr = new TMultiGraph("mgr","");
	388	mgr->Add(gr,"3");
	389	mgr->Add(gr11,"P");
	390	mgr->Draw("A");
	391	mgr->GetXaxis()->SetTitle("#Sigma E_{T} [GeV]");
	392	mgr->GetYaxis()->SetTitle("Resolution of x-component of MET [GeV]");
	393	mgr->GetYaxis()->SetRangeUser(0,60);
	394	mgr->GetXaxis()->SetRangeUser(1,600);
	395	mgr->SetMaximum(16);
	396
[39]	397	Double_t* params = fitfun->GetParameters();
[503]	398
[39]	399	char tempResol[500];
[473]	400	sprintf(tempResol,"E_{x}^{mis} resolution: %.2f #times #sqrt{E_{T}}",params[0]);
[23]	401
[473]	402	TPaveText *leg1 = MakeTPave(0.26,0.84,0.43,0.89,tempResol);
[39]	403	leg1->Draw();
[23]	404
[473]	405	TPaveText *leg2 = MakeTPave(0.15,0.75,0.45,0.8,"Events: pp #rightarrow ggX");
[39]	406	leg2->Draw();
[27]	407
[473]	408	TPaveText *Delphes = MakeTPave(0.65,0.19,0.8,0.24,"MG/ME + Pythia + Delphes");
[39]	409	Delphes->Draw();
[473]	410
[503]	411	TPaveText *cms = MakeTPave(0.65,0.25,0.8,0.30,"CMS-like detector");
	412	cms->Draw();
	413
	414	TLegend *legend = new TLegend(0.69,0.6,0.94,0.74,NULL,"NDC");
	415	legend->AddEntry(gr,"CMS","f");
	416	legend->AddEntry(gr11,"Delphes","l");
	417	legend->SetFillColor(10);
	418	legend->SetBorderSize(0);
	419	legend->Draw();
	420
[473]	421
[39]	422	delete fitfun;
[23]	423	}
	424
	425	void General()
	426	{
[264]	427	JetResol();
	428	ElecResol();
	429	ETmisResol();
	430	//TauJetInfo();
[39]	431
[27]	432	}
[23]	433
	434

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: svn/trunk/routines/resolutions.C@ 733

Download in other formats: