Fork me on GitHub

Context Navigation

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

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

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

Download in other formats: