Fork me on GitHub

Context Navigation

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

Last change on this file since 455 was 443, checked in by Xavier Rouby, 15 years ago
new header in all files
File size: 12.2 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"
[153]	41	#include "TLegend.h"
[9]	42
	43	#include "interface/FuncDef.h"
	44
[27]	45	void JetResol()
[9]	46	{
[39]	47	setTDRStyle();
	48	gROOT->Reset();
	49
[258]	50	TFile *f1 = new TFile("JET2.root","read");
[39]	51	TTree Analyze = (TTree)f1->Get("Analysis");
	52
[153]	53	const Int_t numBin=16;
	54	double bins[numBin]={0,10,20,30,40,50,60,70,80,100,120,140,180,220,300,1200};
[39]	55	TProfile *ETSoverET = new TProfile("ETSoverET","Jet resolution: E_{T}^{rec}/E_{T}^{mc}",(numBin-1),bins,-10,10);
	56
	57	double rms[numBin];
	58	double mean[numBin];
	59
	60	TCanvas *c1 = new TCanvas("c1","JET resol",0,0,1000,650);
	61	c1->cd(); int frame=0;
	62	c1->Divide(6,2);
	63
	64	float x[numBin-1];
	65	float y[numBin-1];
	66	float ex[numBin-1];
	67	float ey[numBin-1];
	68
	69	float finval=0;//valeur a remplir puis a fitter
	70
	71	for ( int i=0; i<(numBin-1); i++) // premiÃšre bin : i ==1 et pas i == 0
	72	{
	73	TAxis *xaxis = ETSoverET->GetXaxis();
	74	float binCenter = xaxis->GetBinCenter(i+1);
	75	int binMin = (int)xaxis->GetBinLowEdge(i+1);
	76	int binMax = (int)xaxis->GetBinUpEdge(i+1);
	77	char tempMin[500];
	78	if(i==0)binMin=5;
[91]	79	sprintf(tempMin,"JetPTResol.PT > %d",binMin);
[39]	80	string mystringMin(tempMin);
	81	char tempMax[500];
[91]	82	sprintf(tempMax,"JetPTResol.PT < %d",binMax);
[39]	83	string mystringMax(tempMax);
	84	char tempName[500];
[91]	85	sprintf(tempName,"(JetPTResol.SmearedPT)>>hETSoverET%d",i);
[39]	86	string mystringName(tempName);
	87
	88	c1->cd(++frame);
	89	GaussValues(Analyze,tempName,rms[i],mean[i],mystringMin,mystringMax);
	90	x[i]=binCenter;
	91	finval=rms[i]/mean[i];
	92	y[i]=(finval*100);
	93	ex[i]=0;
	94	ey[i]=0;
	95	}
	96
	97	TCanvas *c2 = new TCanvas("c2","JET resol",100,100,600,450);
	98	c2->cd();
	99
[153]	100	TF1 *fitfun = new TF1("user","sqrt(pow([0]/x,2)+pow([1]/sqrt(x),2)+pow([2],2))",10,800);
	101	TF1 fitfunCMS = new TF1("userCMS","sqrt(pow(1.05100/sqrt(x),2)+pow(0.06100,2)+pow(0100/x,2))",10,800);
	102	//TF1 fitfunCMS = new TF1("userCMS","sqrt(pow(1.17100/sqrt(x),2)+pow(0.039*100,2))",10,800);
[39]	103
[153]	104	TGraph *gr11 = new TGraph((numBin-1),x,y);
[39]	105	gr11->Draw("AP");
	106	gr11->SetTitle("");
	107	gr11->GetXaxis()->SetTitle("E_{T}^{MC} [GeV]");
	108	gr11->GetYaxis()->SetRangeUser(0,50);
	109	gr11->GetYaxis()->SetTitle("#sigma(E_{T}^{rec}/E_{T}^{MC})_{fit}/<E_{T}^{rec}/E_{T}^{MC}>_{fit}");
	110
	111	Double_t* params = fitfun->GetParameters();
[153]	112
	113	fitfun->SetLineColor(6);
[39]	114	gr11->Fit("user","QR");
[153]	115	gr11->Fit("user","QR");
	116	gr11->Fit("user","QR");
	117	gr11->Fit("user","QR");
	118	char tempResol[100];
	119	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]);
	120	char tempResol2[100];
	121	sprintf(tempResol2,"sqrt(pow(%f/sqrt(x),2)+pow(%f,2))",params[1],params[2]);
[27]	122
[153]	123
	124	TF1 *fitfunDelphes = new TF1("userDelphes",tempResol2,7,1000);
	125	fitfunDelphes->SetLineColor(6);
	126	fitfunDelphes->Draw("same");
	127
	128	//CMS values
	129	fitfunCMS->SetLineColor(9);
	130	fitfunCMS->Draw("same");
	131
[60]	132	TPaveText *events = MakeTPave(0.2,0.75,0.35,0.8,"Events: pp #rightarrow gg ");
	133	events->Draw();
	134
[39]	135	TPaveText *Delphes = MakeTPave(0.2,0.15,0.35,0.2,"MG/ME + Delphes");
	136	Delphes->Draw();
[153]	137
	138	TLegend *legend = new TLegend(0.2,0.6,0.9,0.85,NULL,"NDC");
	139	legend->AddEntry(fitfunCMS,"CMS resolution","l");
	140	legend->AddEntry(fitfunDelphes,tempResol,"l");
	141	legend->SetFillColor(10);
	142	legend->SetBorderSize(0);
	143	legend->Draw();
	144
[39]	145	delete fitfun;
	146	}
[27]	147
[39]	148	void ElecResol()
	149	{
[9]	150
[39]	151	setTDRStyle();
	152	gROOT->Reset();
	153
[153]	154	TFile *f1 = new TFile("ETMIS2.root","read");
[39]	155	TTree Analyze = (TTree)f1->Get("Analysis");
	156
[153]	157	const Int_t numBin=11;
	158	double bins[numBin]={0,10,20,30,40,50,60,70,80,100,120};
[39]	159	TProfile *ETSminusET = new TProfile("ETSminusET","Electron resolution: E_{T}^{rec}/E_{T}^{mc}",(numBin-1),bins,-10,10);
	160
	161	double rms[numBin];
	162	double mean[numBin];
	163
	164	TCanvas *c3 = new TCanvas("c3","ELEC resol",0,0,1000,650);
	165	c3->cd(); int frame=0;
	166	c3->Divide(6,2);
	167
	168	float x[numBin-1];
	169	float y[numBin-1];
	170	float ex[numBin-1];
	171	float ey[numBin-1];
	172
	173	float finval=0;//valeur a remplir puis a fitter
	174
	175	for ( int i=0; i<(numBin-1); i++) // premiÃšre bin : i ==1 et pas i == 0
	176	{
	177	TAxis *xaxis = ETSminusET->GetXaxis();
	178	float binCenter = xaxis->GetBinCenter(i+1);
	179	int binMin = (int)xaxis->GetBinLowEdge(i+1);
	180	int binMax = (int)xaxis->GetBinUpEdge(i+1);
	181	char tempMin[500];
	182	if(i==0)binMin=5;
[89]	183	sprintf(tempMin,"ElecEResol.E > %d",binMin);
[39]	184	string mystringMin(tempMin);
	185	char tempMax[500];
[89]	186	sprintf(tempMax,"ElecEResol.E < %d",binMax);
[39]	187	string mystringMax(tempMax);
	188	char tempName[500];
[89]	189	sprintf(tempName,"(ElecEResol.E-ElecEResol.SmearedE)>>hETSoverET%d",i);
[39]	190	string mystringName(tempName);
	191
	192	c3->cd(++frame);
[153]	193	GaussValuesElec(Analyze,tempName,rms[i],mean[i],mystringMin,mystringMax);
	194	//GaussValues(Analyze,tempName,rms[i],mean[i],mystringMin,mystringMax);
[39]	195	x[i]=binCenter;
	196	finval=rms[i]/binCenter;
	197	y[i]=(finval*100);
	198	ex[i]=0;
	199	ey[i]=0;
	200	}
	201
	202	TCanvas *c4 = new TCanvas("c4","ELEC resol",100,100,600,450);
	203	c4->cd();
	204
	205	TF1 *fitfun = new TF1("user","sqrt(pow([0],2)+pow([1]/sqrt(x),2)+pow([2]/x,2))",1,400);
	206
	207	TGraphErrors *gr11 = new TGraphErrors((numBin-1),x,y,ex,ey);
	208	gr11->Draw("AP");
	209	gr11->SetTitle("");
	210	gr11->GetXaxis()->SetTitle("E [GeV]");
	211	gr11->GetYaxis()->SetRangeUser(0,5);
	212	gr11->GetYaxis()->SetTitle("#sigma/E");
	213
	214	Double_t* params = fitfun->GetParameters();
	215
	216	gr11->Fit("user","QR");
	217	gr11->Fit("user","QRI");
	218	gr11->Fit("user","QRI");
	219	gr11->Fit("user","QRI");
	220	char tempResol[500];
	221	sprintf(tempResol,"#frac{#sigma}{E} = #frac{%f}{#sqrt{E}} #oplus #frac{%f}{E} #oplus %f",params[1]/100,params[2]/100,params[0]/100);
	222
	223	TPaveText *leg1 = MakeTPave(0.4,0.6,0.8,0.65,tempResol);
	224	leg1->Draw();
	225
	226	TPaveText *Delphes = MakeTPave(0.2,0.15,0.35,0.2,"MG/ME + Delphes");
	227	Delphes->Draw();
	228
[60]	229	TPaveText *events = MakeTPave(0.2,0.75,0.35,0.8,"Events: WHq'#rightarrow W#tau#tauq'#rightarrowjjl#tauq', m_{H}=150 GeV ");
[39]	230	events->Draw();
	231
	232	delete fitfun;
	233	}
[9]	234
[39]	235	void TauJetInfo()
	236	{
	237
	238	setTDRStyle();
	239	gROOT->Reset();
	240
[258]	241	TFile *f1 = new TFile("TAUJET2.root","read");
[39]	242	TTree Analyze = (TTree)f1->Get("Analysis");
	243
[60]	244	TCanvas *ct1 = new TCanvas("ct1","Tau information",100,100,600,450);
	245	ct1->cd();
[39]	246
	247	TH1F *tauEnergy =MakeNormTH1F(20,0.8,1,Analyze,"TauJetPTResol.EnergieCen>>tauEnergy",1, 0, 1,2,false);
	248	tauEnergy->Draw();
	249	tauEnergy->SetTitle("");
	250	tauEnergy->GetYaxis()->SetTitle("Fraction of events");
	251	tauEnergy->GetXaxis()->SetTitle("C_{#tau}");
	252
	253
	254	TPaveText *Delphes1 = MakeTPave(0.3,0.85,0.45,0.9,"MG/ME + Delphes");
	255	Delphes1->Draw();
	256	TPaveText *ctau = MakeTPave(0.3,0.75,0.45,0.8,"C_{#tau} = #frac{#sum E^{towers} (#DeltaR = 0.15)}{E^{jet}}");
	257	ctau->Draw();
[153]	258	TPaveText *events1 = MakeTPave(0.3,0.65,0.45,0.7,"Events: WHq'#rightarrow WWWq'#rightarrow lllq', m_{H}=150 GeV ");
[39]	259	events1->Draw();
	260
	261
[60]	262	TCanvas *ct2 = new TCanvas("ct2","Tau information",100,100,600,450);
	263	ct2->cd();
[39]	264	TH1F *NumTrack =MakeNormTH1F(6,0,6,Analyze,"TauJetPTResol.NumTrack>>NumTrack",1, 0, 1,2,false);
	265	NumTrack->Draw();
	266	NumTrack->SetTitle("");
	267	NumTrack->GetYaxis()->SetTitle("Fraction of events");
	268	NumTrack->GetXaxis()->SetTitle("N^{tracks}");
	269
	270	TPaveText *Delphes = MakeTPave(0.6,0.85,0.85,0.9,"MG/ME + Delphes");
	271	Delphes->Draw();
	272	TPaveText *numtracks = MakeTPave(0.6,0.75,0.85,0.8,"#DeltaR < 0.4, p_{T}^{track} > 2 GeV");
	273	numtracks->Draw();
[153]	274	TPaveText *events = MakeTPave(0.6,0.65,0.85,0.7,"Events: WHq'#rightarrow WWWq'#rightarrow lllq', m_{H}=150 GeV ");
[39]	275	events->Draw();
	276
[9]	277	}
[23]	278
[27]	279	void ETmisResol()
[23]	280	{
[39]	281
	282	setTDRStyle();
	283	gROOT->Reset();
	284
[258]	285	TFile *f1 = new TFile("JET2.root","read");
[39]	286	TTree Analyze = (TTree)f1->Get("Analysis");
	287
	288	TF1 fitfun = new TF1("user","[0]sqrt(x)",0,600);
[153]	289	const Int_t numBin=6;
	290	double bins[numBin]={180,260,340,420,500,580};
[39]	291	TProfile *ETSoverET = new TProfile("ETSoverET","ETmis resol",(numBin-1),bins,-1000,1000);
	292
	293
	294	double rms[numBin-1];
	295
	296	TCanvas *c5 = new TCanvas("c5","PTmis resol",0,0,1000,650);
	297	c5->cd(); int frame=0;
	298	c5->Divide(6,2);
	299
	300	double x[numBin];
	301	double y[numBin];
	302
	303	for ( int i=0; i<(numBin-1); i++) // premiÃšre bin : i ==1 et pas i == 0
	304	{
	305	TAxis *xaxis = ETSoverET->GetXaxis();
	306	float binCenter = xaxis->GetBinCenter(i+1);
	307	int binMin = (int)xaxis->GetBinLowEdge(i+1);
	308	int binMax = (int)xaxis->GetBinUpEdge(i+1);
	309	char tempMin[500];
	310	sprintf(tempMin,"ETmisResol.SEt>%d",binMin);
	311	string mystringMin(tempMin);
	312	char tempMax[500];
	313	sprintf(tempMax,"ETmisResol.SEt<%d",binMax);
	314	string mystringMax(tempMax);
	315
	316	char tempName[500];
	317	sprintf(tempName,"ETmisResol.ExSmeare>>hETSoverET%d",i);
	318	string mystringName(tempName);
	319
	320	c5->cd(++frame);
	321	GaussValuesETmis(Analyze,tempName,rms[i],mystringMin,mystringMax);
	322	x[i]=binCenter;
	323	y[i]=rms[i];
	324
	325	}
	326
	327	TCanvas *c6 = new TCanvas("c6","ETmis resolution",100,100,600,450);
	328	c6->cd();
	329
	330	x[numBin]=0;
	331	y[numBin]=0;
	332	TGraph *gr11 = new TGraph((numBin),x,rms);
	333	gr11->Draw("AP");
	334	gr11->GetXaxis()->SetTitle("Offline sum of E_{T} [GeV]");
	335	gr11->GetYaxis()->SetTitle("Resolution of x-component of MET [GeV]");
	336	gr11->Fit("user","RQ");
	337	gr11->Fit("user","RQ");
	338	gr11->Fit("user","RQ");
	339	gr11->Fit("user","RQ");
[153]	340	gr11->GetYaxis()->SetRangeUser(0,60);
[39]	341	gr11->GetXaxis()->SetRangeUser(1,600);
[23]	342
[39]	343	Double_t* params = fitfun->GetParameters();
	344
	345	char tempResol[500];
	346	sprintf(tempResol,"%f * #sqrt{E_{T}}",params[0]);
[23]	347
[39]	348	TPaveText *leg1 = MakeTPave(0.4,0.6,0.8,0.65,tempResol);
	349	leg1->Draw();
[23]	350
[153]	351	TPaveText *leg2 = MakeTPave(0.2,0.8,0.8,0.85,"WHq'#rightarrow WWWq'#rightarrow lllq', m_{H}=150 GeV ");
[39]	352	leg2->Draw();
[27]	353
[39]	354	TPaveText *Delphes = MakeTPave(0.2,0.15,0.35,0.2,"MG/ME + Delphes");
	355	Delphes->Draw();
	356
	357
	358	delete fitfun;
[23]	359	}
	360
	361	void General()
	362	{
[264]	363	JetResol();
	364	ElecResol();
	365	ETmisResol();
	366	//TauJetInfo();
[39]	367
[27]	368	}
[23]	369
	370

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

Download in other formats: