Context Navigation

SmearUtil.cc@ 752

Last change on this file since 752 was 671, checked in by cp3-support, 12 years ago
corrected b-tagging by a.mertens
File size: 96.2 KB

Rev	Line
[260]	1	/***********************************************************************
	2
	3	/----------------------------------------------\
	4	\| Delphes, a framework for the fast simulation \|
	5	\| of a generic collider experiment \|
[443]	6	\------------- arXiv:0903.2225v1 ------------/
[260]	7
	8
	9	This package uses:
	10	------------------
[443]	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]
[260]	14	FROG: [hep-ex/0901.2718v1]
[443]	15	HepMC: Comput. Phys. Commun.134 (2001) 41
[260]	16
	17	------------------------------------------------------------------
	18
	19	Main authors:
	20	-------------
	21
[443]	22	Severine Ovyn Xavier Rouby
	23	severine.ovyn@uclouvain.be xavier.rouby@cern
[260]	24
[443]	25	Center for Particle Physics and Phenomenology (CP3)
	26	Universite catholique de Louvain (UCL)
	27	Louvain-la-Neuve, Belgium
	28
[567]	29	Copyright (C) 2008-2010,
[443]	30	All rights reserved.
[260]	31
	32	***********************************************************************/
[2]	33
	34	/// \file SmearUtil.cc
	35	/// \brief RESOLution class, and some generic definitions
	36
	37
[219]	38	#include "SmearUtil.h"
[399]	39	#include "TStopwatch.h"
[558]	40	#include "TF2.h"
[2]	41
	42	#include <iostream>
[219]	43	#include <fstream>
[2]	44	#include <sstream>
[44]	45	#include <iomanip>
[380]	46	#include <map>
[454]	47	#include <vector>
[526]	48	#include <cmath>
[537]	49	#include <cstdlib> // for exit()
[219]	50	using namespace std;
[44]	51
[2]	52	//------------------------------------------------------------------------------
	53
	54	RESOLution::RESOLution() {
	55
[94]	56	// Detector characteristics
	57	CEN_max_tracker = 2.5; // Maximum tracker coverage
[494]	58	CEN_max_calo_cen = 1.7; // central calorimeter coverage
	59	CEN_max_calo_ec = 3.0; // calorimeter endcap coverage
[94]	60	CEN_max_calo_fwd = 5.0; // forward calorimeter pseudorapidity coverage
	61	CEN_max_mu = 2.4; // muon chambers pseudorapidity coverage
	62
	63	// Energy resolution for electron/photon
	64	// \sigma/E = C + N/E + S/\sqrt{E}
	65	ELG_Scen = 0.05; // S term for central ECAL
[494]	66	ELG_Ncen = 0.25; // N term
	67	ELG_Ccen = 0.005; // C term
	68	ELG_Sec = 0.05; // S term for central ECAL endcap
	69	ELG_Nec = 0.25; // S term
	70	ELG_Cec = 0.005; // S term
[257]	71	ELG_Sfwd = 2.084; // S term for FCAL
[494]	72	ELG_Nfwd = 0.0; // N term
	73	ELG_Cfwd = 0.107; // C term
[374]	74	ELG_Szdc = 0.70; // S term for ZDC
[494]	75	ELG_Nzdc = 0.0; // N term
	76	ELG_Czdc = 0.08; // C term
[2]	77
[494]	78	// Energy resolution for hadrons in ecal/hcal/fwd
[94]	79	// \sigma/E = C + N/E + S/\sqrt{E}
[494]	80	HAD_Scen = 1.5; // S term for central HCAL
	81	HAD_Ncen = 0.; // N term
	82	HAD_Ccen = 0.05; // C term
	83	HAD_Sec = 1.5; // S term for HCAL endcap
	84	HAD_Nec = 0.; // N term
	85	HAD_Cec = 0.05; // C term
	86	HAD_Sfwd = 2.7; // S term for FCAL
	87	HAD_Nfwd = 0.; // N term
	88	HAD_Cfwd = 0.13; // C term
[374]	89	HAD_Szdc = 1.38; // S term for ZDC
[494]	90	HAD_Nzdc = 0.; // N term
	91	HAD_Czdc = 0.13; // C term
[2]	92
[94]	93	// Muon smearing
	94	MU_SmearPt = 0.01;
[2]	95
[374]	96	// time resolution
	97	ZDC_T_resolution = 0; // resolution for time measurement [s]
	98	RP220_T_resolution = 0;
	99	RP420_T_resolution = 0;
	100
[94]	101	// Tracking efficiencies
	102	TRACK_ptmin = 0.9; // minimal pt needed to reach the calorimeter in GeV
	103	TRACK_eff = 100; // efficiency associated to the tracking
[2]	104
[94]	105	// Calorimetric towers
	106	TOWER_number = 40;
	107	const float tower_eta_edges[41] = {
	108	0., 0.087, 0.174, 0.261, 0.348, 0.435, 0.522, 0.609, 0.696, 0.783, 0.870, 0.957, 1.044, 1.131, 1.218, 1.305, 1.392, 1.479, 1.566,
	109	1.653, 1.740, 1.830, 1.930, 2.043, 2.172, 2.322, 2.500, 2.650, 2.868, 2.950, 3.125, 3.300, 3.475, 3.650, 3.825, 4.000, 4.175,
	110	4.350, 4.525, 4.700, 5.000}; // temporary object
	111	TOWER_eta_edges = new float[TOWER_number+1];
	112	for(unsigned int i=0; i<TOWER_number +1; i++) TOWER_eta_edges[i] = tower_eta_edges[i];
	113
	114	const float tower_dphi[40] = {
	115	5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 10,
	116	10,10,10,10,10, 10,10,10,10,10, 10,10,10,10,10, 10,10,10,20, 20 }; // temporary object
	117	TOWER_dphi = new float[TOWER_number];
	118	for(unsigned int i=0; i<TOWER_number; i++) TOWER_dphi[i] = tower_dphi[i];
[2]	119
	120
[374]	121	// Thresholds for reconstructed objetcs (GeV)
[94]	122	PTCUT_elec = 10.0;
	123	PTCUT_muon = 10.0;
	124	PTCUT_jet = 20.0;
	125	PTCUT_gamma = 10.0;
	126	PTCUT_taujet = 10.0;
[33]	127
[374]	128	ZDC_gamma_E = 20; // GeV
	129	ZDC_n_E = 50; // GeV
	130
[321]	131	// Isolation
[557]	132	ISOL_trk_PT = 2.0; //minimal pt of tracks for isolation criteria
	133	ISOL_trk_Cone = 0.5; //Cone for isolation criteria
	134	ISOL_calo_Cone = 0.5; //Cone for isolation criteria
	135	ISOL_calo_ET = 1E99; //minimal tower energy for isolation criteria. Default off = 1E99
	136	ISOL_calo_Grid = 3; //Grid size (N x N) for calorimetric isolation -- should be odd
[305]	137
[94]	138	// General jet variable
	139	JET_coneradius = 0.7; // generic jet radius ; not for tau's !!!
[567]	140	JET_jetalgo = 2; // 1 for Cone algorithm, 2 for MidPoint algorithm, 3 for SIScone algorithm, 4 for kt algorithm
[94]	141	JET_seed = 1.0; // minimum seed to start jet reconstruction
[567]	142	JET_Eflow = 0; // 1 for Energy flow in jets reco ; 0 if not
[33]	143
[94]	144	// Tagging definition
[558]	145	BTAG_func_b = "0.4 + 0x + 0y"; // efficiency functions (x=Pt, y=eta) in ROOT::TF2 format
	146	BTAG_func_c = "0.1 + 0x + 0y";
	147	BTAG_func_l = "0.01 + 0x + 0y";
	148	BTAG_func_g = "0.01 + 0x + 0y";
[2]	149
[558]	150
[94]	151	// FLAGS
	152	FLAG_bfield = 1; //1 to run the bfield propagation else 0
	153	FLAG_vfd = 1; //1 to run the very forward detectors else 0
[307]	154	FLAG_RP = 1; //1 to run the zero degree calorimeter else 0
[94]	155	FLAG_trigger = 1; //1 to run the trigger selection else 0
	156	FLAG_frog = 1; //1 to run the FROG event display
[307]	157	FLAG_lhco = 1;
[2]	158
[94]	159	// In case BField propagation allowed
	160	TRACK_radius = 129; //radius of the BField coverage
	161	TRACK_length = 300; //length of the BField coverage
	162	TRACK_bfield_x = 0; //X composant of the BField
	163	TRACK_bfield_y = 0; //Y composant of the BField
	164	TRACK_bfield_z = 3.8; //Z composant of the BField
[2]	165
[94]	166	// In case Very forward detectors allowed
	167	VFD_min_calo_vfd = 5.2; // very forward calorimeter (if any) like CASTOR
	168	VFD_max_calo_vfd = 6.6;
	169	VFD_min_zdc = 8.3;
	170	VFD_s_zdc = 140; // distance of the Zero Degree Calorimeter, from the Interaction poin, in [m]
[2]	171
[94]	172	RP_220_s = 220; // distance of the RP to the IP, in meters
	173	RP_220_x = 0.002; // distance of the RP to the beam, in meters
	174	RP_420_s = 420; // distance of the RP to the IP, in meters
	175	RP_420_x = 0.004; // distance of the RP to the beam, in meters
[257]	176	RP_IP_name = "IP5";
[252]	177	RP_beam1Card = "data/LHCB1IR5_v6.500.tfs";
	178	RP_beam2Card = "data/LHCB1IR5_v6.500.tfs";
[2]	179
[94]	180	// In case FROG event display allowed
	181	NEvents_Frog = 10;
[2]	182
[422]	183	// Number of events to be processed
	184	NEvents = -1;
	185
[94]	186	//********************************************
	187	//jet stuffs not defined in the input datacard
	188	//********************************************
	189
	190	JET_overlap = 0.75;
	191	// MidPoint algorithm definition
	192	JET_M_coneareafraction = 0.25;
	193	JET_M_maxpairsize = 2;
	194	JET_M_maxiterations = 100;
	195	// Define Cone algorithm.
	196	JET_C_adjacencycut = 2;
	197	JET_C_maxiterations = 100;
	198	JET_C_iratch = 1;
	199	//Define SISCone algorithm.
	200	JET_S_npass = 0;
	201	JET_S_protojet_ptmin= 0.0;
	202
	203	//For Tau-jet definition
	204	TAU_energy_scone = 0.15; // radius R of the cone for tau definition, based on energy threshold
	205	TAU_track_scone = 0.4; // radius R of the cone for tau definition, based on track number
	206	TAU_track_pt = 2; // minimal pt [GeV] for tracks to be considered in tau definition
	207	TAU_energy_frac = 0.95; // fraction of energy required in the central part of the cone, for tau jets
	208
	209	PT_QUARKS_MIN = 2.0 ; // minimal pt needed by quarks to do b-tag
[252]	210
	211	//for very forward detectors
[399]	212	RP_offsetEl_s = 120; // distance of beam separation point, from IP
[404]	213	RP_offsetEl_x = -0.097; // half distance of separation of beams in horizontal plan, in m
	214	RP_offsetEl_y = 0; // half distance of separation of beams in vertical plan, in m
[399]	215	RP_cross_x = -500; // IP offset in horizontal plane, in micrometers
	216	RP_cross_y = 0.0; // IP offset in vertical plane, in micrometers
	217	RP_cross_ang_x = 142.5; // half-crossing angle in horizontal plane, in microrad
	218	RP_cross_ang_y = 0.0; // half-crossing angle in vertical plane, in microrad
[380]	219
[399]	220
[380]	221	PdgTableFilename = "data/particle.tbl";
[454]	222	inputfilelist = "";
	223	detectorcard = "";
	224	triggercard = "";
[526]	225
	226	grandom = new TRandom3(0); // a new seed is set everytime Delphes is run
[2]	227	}
	228
[219]	229
	230	RESOLution::RESOLution(const RESOLution & DET) {
	231	// Detector characteristics
	232	CEN_max_tracker = DET.CEN_max_tracker;
	233	CEN_max_calo_cen = DET.CEN_max_calo_cen;
[494]	234	CEN_max_calo_ec = DET.CEN_max_calo_ec;
[219]	235	CEN_max_calo_fwd = DET.CEN_max_calo_fwd;
	236	CEN_max_mu = DET.CEN_max_mu;
	237
	238	// Energy resolution for electron/photon
	239	ELG_Scen = DET.ELG_Scen;
	240	ELG_Ncen = DET.ELG_Ncen;
	241	ELG_Ccen = DET.ELG_Ccen;
[494]	242	ELG_Sec = DET.ELG_Sec;
	243	ELG_Nec = DET.ELG_Nec;
	244	ELG_Cec = DET.ELG_Cec;
[219]	245	ELG_Cfwd = DET.ELG_Cfwd;
	246	ELG_Sfwd = DET.ELG_Sfwd;
	247	ELG_Nfwd = DET.ELG_Nfwd;
[374]	248	ELG_Czdc = DET.ELG_Czdc;
	249	ELG_Szdc = DET.ELG_Szdc;
	250	ELG_Nzdc = DET.ELG_Nzdc;
[219]	251
[494]	252	// Energy resolution for hadrons in ecal/hcal/fwd/zdc
	253	HAD_Scen = DET.HAD_Scen;
	254	HAD_Ncen = DET.HAD_Ncen;
	255	HAD_Ccen = DET.HAD_Ccen;
	256	HAD_Sec = DET.HAD_Sec;
	257	HAD_Nec = DET.HAD_Nec;
	258	HAD_Cec = DET.HAD_Cec;
	259	HAD_Sfwd = DET.HAD_Sfwd;
	260	HAD_Nfwd = DET.HAD_Nfwd;
	261	HAD_Cfwd = DET.HAD_Cfwd;
[374]	262	HAD_Szdc = DET.HAD_Szdc;
	263	HAD_Nzdc = DET.HAD_Nzdc;
	264	HAD_Czdc = DET.HAD_Czdc;
[219]	265
[374]	266	// time resolution
	267	ZDC_T_resolution = DET.ZDC_T_resolution; // resolution for time measurement [s]
	268	RP220_T_resolution = DET.RP220_T_resolution;
	269	RP420_T_resolution = DET.RP420_T_resolution;
	270
[219]	271	// Muon smearing
	272	MU_SmearPt = DET.MU_SmearPt;
	273
	274	// Tracking efficiencies
	275	TRACK_ptmin = DET.TRACK_ptmin;
	276	TRACK_eff = DET.TRACK_eff;
	277
	278	// Calorimetric towers
	279	TOWER_number = DET.TOWER_number;
	280	TOWER_eta_edges = new float[TOWER_number+1];
	281	for(unsigned int i=0; i<TOWER_number +1; i++) TOWER_eta_edges[i] = DET.TOWER_eta_edges[i];
	282
	283	TOWER_dphi = new float[TOWER_number];
	284	for(unsigned int i=0; i<TOWER_number; i++) TOWER_dphi[i] = DET.TOWER_dphi[i];
	285
	286	// Thresholds for reconstructed objetcs
	287	PTCUT_elec = DET.PTCUT_elec;
	288	PTCUT_muon = DET.PTCUT_muon;
	289	PTCUT_jet = DET.PTCUT_jet;
	290	PTCUT_gamma = DET.PTCUT_gamma;
	291	PTCUT_taujet = DET.PTCUT_taujet;
	292
[374]	293	ZDC_gamma_E = DET.ZDC_gamma_E;
	294	ZDC_n_E = DET.ZDC_n_E;
	295
[321]	296	// Isolation
[557]	297	ISOL_trk_PT = DET.ISOL_trk_PT; // tracking isolation
	298	ISOL_trk_Cone = DET.ISOL_trk_Cone;
	299	ISOL_calo_Cone = DET.ISOL_calo_Cone;
	300	ISOL_calo_ET = DET.ISOL_calo_ET; // calorimeter isolation, defaut off
	301	ISOL_calo_Grid = DET.ISOL_calo_Grid;
[305]	302
	303
[219]	304	// General jet variable
	305	JET_coneradius = DET.JET_coneradius;
	306	JET_jetalgo = DET.JET_jetalgo;
	307	JET_seed = DET.JET_seed;
[383]	308	JET_Eflow = DET.JET_Eflow;
[219]	309
	310	// Tagging definition
[558]	311	BTAG_func_b = DET.BTAG_func_b;
	312	BTAG_func_c = DET.BTAG_func_c;
	313	BTAG_func_l = DET.BTAG_func_l;
	314	BTAG_func_g = DET.BTAG_func_g;
[219]	315
[558]	316
[219]	317	// FLAGS
	318	FLAG_bfield = DET.FLAG_bfield;
	319	FLAG_vfd = DET.FLAG_vfd;
[306]	320	FLAG_RP = DET.FLAG_RP;
[219]	321	FLAG_trigger = DET.FLAG_trigger;
	322	FLAG_frog = DET.FLAG_frog;
[307]	323	FLAG_lhco = DET.FLAG_lhco;
[219]	324
	325	// In case BField propagation allowed
	326	TRACK_radius = DET.TRACK_radius;
	327	TRACK_length = DET.TRACK_length;
	328	TRACK_bfield_x = DET.TRACK_bfield_x;
	329	TRACK_bfield_y = DET.TRACK_bfield_y;
	330	TRACK_bfield_z = DET.TRACK_bfield_z;
	331
	332	// In case Very forward detectors allowed
	333	VFD_min_calo_vfd = DET.VFD_min_calo_vfd;
	334	VFD_max_calo_vfd = DET.VFD_max_calo_vfd;
	335	VFD_min_zdc = DET.VFD_min_zdc;
	336	VFD_s_zdc = DET.VFD_s_zdc;
	337
	338	RP_220_s = DET.RP_220_s;
	339	RP_220_x = DET.RP_220_x;
	340	RP_420_s = DET.RP_420_s;
	341	RP_420_x = DET.RP_420_x;
[252]	342	RP_beam1Card = DET.RP_beam1Card;
	343	RP_beam2Card = DET.RP_beam2Card;
	344	RP_offsetEl_s = DET.RP_offsetEl_s;
	345	RP_offsetEl_x = DET.RP_offsetEl_x;
[404]	346	RP_offsetEl_y = DET.RP_offsetEl_y;
[254]	347	RP_cross_x = DET.RP_cross_x;
	348	RP_cross_y = DET.RP_cross_y;
[399]	349	RP_cross_ang_x = DET.RP_cross_ang_x;
	350	RP_cross_ang_y = DET.RP_cross_ang_y;
[257]	351	RP_IP_name = DET.RP_IP_name;
[219]	352
	353	// In case FROG event display allowed
	354	NEvents_Frog = DET.NEvents_Frog;
	355
[422]	356	// Number of events to be processed
	357	NEvents = DET.NEvents;
	358
[219]	359	JET_overlap = DET.JET_overlap;
	360	// MidPoint algorithm definition
	361	JET_M_coneareafraction = DET.JET_M_coneareafraction;
	362	JET_M_maxpairsize = DET.JET_M_maxpairsize;
	363	JET_M_maxiterations = DET.JET_M_maxiterations;
	364	// Define Cone algorithm.
	365	JET_C_adjacencycut = DET.JET_C_adjacencycut;
	366	JET_C_maxiterations = DET.JET_C_maxiterations;
	367	JET_C_iratch = DET.JET_C_iratch;
	368	//Define SISCone algorithm.
	369	JET_S_npass = DET.JET_S_npass;
	370	JET_S_protojet_ptmin = DET.JET_S_protojet_ptmin;
	371
	372	//For Tau-jet definition
	373	TAU_energy_scone = DET.TAU_energy_scone;
	374	TAU_track_scone = DET.TAU_track_scone;
	375	TAU_track_pt = DET.TAU_track_pt;
	376	TAU_energy_frac = DET.TAU_energy_frac;
	377
	378	PT_QUARKS_MIN = DET.PT_QUARKS_MIN;
[380]	379	PdgTableFilename = DET.PdgTableFilename;
	380	PDGtable = DET.PDGtable;
[454]	381	inputfilelist = DET.inputfilelist;
	382	detectorcard = DET.detectorcard;
	383	triggercard = DET.triggercard;
[526]	384
	385	grandom = new TRandom3(*(DET.grandom));
[219]	386	}
	387
	388	RESOLution& RESOLution::operator=(const RESOLution& DET) {
	389	if(this==&DET) return *this;
	390	// Detector characteristics
	391	CEN_max_tracker = DET.CEN_max_tracker;
	392	CEN_max_calo_cen = DET.CEN_max_calo_cen;
[494]	393	CEN_max_calo_ec = DET.CEN_max_calo_ec;
[219]	394	CEN_max_calo_fwd = DET.CEN_max_calo_fwd;
	395	CEN_max_mu = DET.CEN_max_mu;
	396
	397	// Energy resolution for electron/photon
	398	ELG_Scen = DET.ELG_Scen;
	399	ELG_Ncen = DET.ELG_Ncen;
	400	ELG_Ccen = DET.ELG_Ccen;
[494]	401	ELG_Sec = DET.ELG_Sec;
	402	ELG_Nec = DET.ELG_Nec;
	403	ELG_Cec = DET.ELG_Cec;
[219]	404	ELG_Cfwd = DET.ELG_Cfwd;
	405	ELG_Sfwd = DET.ELG_Sfwd;
	406	ELG_Nfwd = DET.ELG_Nfwd;
[374]	407	ELG_Czdc = DET.ELG_Czdc;
	408	ELG_Szdc = DET.ELG_Szdc;
	409	ELG_Nzdc = DET.ELG_Nzdc;
[219]	410
[494]	411	// Energy resolution for hadrons in ecal/hcal/fwd/zdc
	412	HAD_Scen = DET.HAD_Scen ;
	413	HAD_Ncen = DET.HAD_Ncen;
	414	HAD_Ccen = DET.HAD_Ccen;
	415	HAD_Sec = DET.HAD_Sec;
	416	HAD_Nec = DET.HAD_Nec;
	417	HAD_Cec = DET.HAD_Cec;
	418	HAD_Sfwd = DET.HAD_Sfwd;
	419	HAD_Nfwd = DET.HAD_Nfwd;
	420	HAD_Cfwd = DET.HAD_Cfwd;
[374]	421	HAD_Szdc = DET.HAD_Szdc;
	422	HAD_Nzdc = DET.HAD_Nzdc;
	423	HAD_Czdc = DET.HAD_Czdc;
[219]	424
[374]	425	// time resolution
	426	ZDC_T_resolution = DET.ZDC_T_resolution; // resolution for time measurement [s]
	427	RP220_T_resolution = DET.RP220_T_resolution;
	428	RP420_T_resolution = DET.RP420_T_resolution;
	429
[219]	430	// Muon smearing
	431	MU_SmearPt = DET.MU_SmearPt;
	432
	433	// Tracking efficiencies
	434	TRACK_ptmin = DET.TRACK_ptmin;
	435	TRACK_eff = DET.TRACK_eff;
	436
	437	// Calorimetric towers
	438	TOWER_number = DET.TOWER_number;
	439	TOWER_eta_edges = new float[TOWER_number+1];
	440	for(unsigned int i=0; i<TOWER_number +1; i++) TOWER_eta_edges[i] = DET.TOWER_eta_edges[i];
	441
	442	TOWER_dphi = new float[TOWER_number];
	443	for(unsigned int i=0; i<TOWER_number; i++) TOWER_dphi[i] = DET.TOWER_dphi[i];
	444
	445	// Thresholds for reconstructed objetcs
	446	PTCUT_elec = DET.PTCUT_elec;
	447	PTCUT_muon = DET.PTCUT_muon;
	448	PTCUT_jet = DET.PTCUT_jet;
	449	PTCUT_gamma = DET.PTCUT_gamma;
	450	PTCUT_taujet = DET.PTCUT_taujet;
	451
[374]	452	ZDC_gamma_E = DET.ZDC_gamma_E;
	453	ZDC_n_E = DET.ZDC_n_E;
	454
[321]	455	// Isolation
[557]	456	ISOL_trk_PT = DET.ISOL_trk_PT; // tracking isolation
	457	ISOL_trk_Cone = DET.ISOL_trk_Cone;
	458	ISOL_calo_Cone = DET.ISOL_calo_Cone;
	459	ISOL_calo_ET = DET.ISOL_calo_ET; // calorimeter isolation, defaut off
	460	ISOL_calo_Grid = DET.ISOL_calo_Grid;
[305]	461
[219]	462	// General jet variable
	463	JET_coneradius = DET.JET_coneradius;
	464	JET_jetalgo = DET.JET_jetalgo;
	465	JET_seed = DET.JET_seed;
[383]	466	JET_Eflow = DET.JET_Eflow;
[219]	467
	468	// Tagging definition
[558]	469	BTAG_func_b = DET.BTAG_func_b;
	470	BTAG_func_c = DET.BTAG_func_c;
	471	BTAG_func_l = DET.BTAG_func_l;
	472	BTAG_func_g = DET.BTAG_func_g;
[219]	473
	474	// FLAGS
	475	FLAG_bfield = DET.FLAG_bfield;
	476	FLAG_vfd = DET.FLAG_vfd;
[306]	477	FLAG_RP = DET.FLAG_RP;
[219]	478	FLAG_trigger = DET.FLAG_trigger;
	479	FLAG_frog = DET.FLAG_frog;
[307]	480	FLAG_lhco = DET.FLAG_lhco;
[219]	481
	482	// In case BField propagation allowed
	483	TRACK_radius = DET.TRACK_radius;
	484	TRACK_length = DET.TRACK_length;
	485	TRACK_bfield_x = DET.TRACK_bfield_x;
	486	TRACK_bfield_y = DET.TRACK_bfield_y;
	487	TRACK_bfield_z = DET.TRACK_bfield_z;
	488
	489	// In case Very forward detectors allowed
	490	VFD_min_calo_vfd = DET.VFD_min_calo_vfd;
	491	VFD_max_calo_vfd = DET.VFD_max_calo_vfd;
	492	VFD_min_zdc = DET.VFD_min_zdc;
	493	VFD_s_zdc = DET.VFD_s_zdc;
	494
	495	RP_220_s = DET.RP_220_s;
	496	RP_220_x = DET.RP_220_x;
	497	RP_420_s = DET.RP_420_s;
	498	RP_420_x = DET.RP_420_x;
[252]	499	RP_offsetEl_s = DET.RP_offsetEl_s;
	500	RP_offsetEl_x = DET.RP_offsetEl_x;
[404]	501	RP_offsetEl_y = DET.RP_offsetEl_y;
[252]	502	RP_beam1Card = DET.RP_beam1Card;
	503	RP_beam2Card = DET.RP_beam2Card;
[254]	504	RP_cross_x = DET.RP_cross_x;
	505	RP_cross_y = DET.RP_cross_y;
[399]	506	RP_cross_ang_x = DET.RP_cross_ang_x;
	507	RP_cross_ang_y = DET.RP_cross_ang_y;
[257]	508	RP_IP_name = DET.RP_IP_name;
[219]	509
[252]	510
[219]	511	// In case FROG event display allowed
	512	NEvents_Frog = DET.NEvents_Frog;
	513
[422]	514	// Number of events to be processed
	515	NEvents = DET.NEvents;
	516
[219]	517	JET_overlap = DET.JET_overlap;
	518	// MidPoint algorithm definition
	519	JET_M_coneareafraction = DET.JET_M_coneareafraction;
	520	JET_M_maxpairsize = DET.JET_M_maxpairsize;
	521	JET_M_maxiterations = DET.JET_M_maxiterations;
	522	// Define Cone algorithm.
	523	JET_C_adjacencycut = DET.JET_C_adjacencycut;
	524	JET_C_maxiterations = DET.JET_C_maxiterations;
	525	JET_C_iratch = DET.JET_C_iratch;
	526	//Define SISCone algorithm.
	527	JET_S_npass = DET.JET_S_npass;
	528	JET_S_protojet_ptmin = DET.JET_S_protojet_ptmin;
	529
	530	//For Tau-jet definition
	531	TAU_energy_scone = DET.TAU_energy_scone;
	532	TAU_track_scone = DET.TAU_track_scone;
	533	TAU_track_pt = DET.TAU_track_pt;
	534	TAU_energy_frac = DET.TAU_energy_frac;
	535
	536	PT_QUARKS_MIN = DET.PT_QUARKS_MIN;
[380]	537
	538	PdgTableFilename = DET.PdgTableFilename;
	539	PDGtable = DET.PDGtable;
[454]	540
	541	inputfilelist = DET.inputfilelist;
	542	detectorcard = DET.detectorcard;
	543	triggercard = DET.triggercard;
	544
[526]	545	grandom = new TRandom3(*(DET.grandom));
	546
[219]	547	return *this;
	548	}
	549
[454]	550	void RESOLution::setNames(const string& list, const string& det, const string& trig) {
	551	inputfilelist = list;
	552	detectorcard = det;
	553	triggercard = trig;
	554	}
[219]	555
[2]	556	//------------------------------------------------------------------------------
	557	void RESOLution::ReadDataCard(const string datacard) {
	558
	559	string temp_string;
	560	istringstream curstring;
	561
	562	ifstream fichier_a_lire(datacard.c_str());
	563	if(!fichier_a_lire.good()) {
[249]	564	cout <<" WARNING: Datadard not found, use default values " << endl;
[94]	565	return;
[2]	566	}
[494]	567	bool CEN_max_calo_ec_flag = false;
[94]	568
[2]	569	while (getline(fichier_a_lire,temp_string)) {
	570	curstring.clear(); // needed when using several times istringstream::str(string)
	571	curstring.str(temp_string);
	572	string varname;
[252]	573	float value; int ivalue; string svalue;
[494]	574
[2]	575	if(strstr(temp_string.c_str(),"#")) { }
[94]	576	else if(strstr(temp_string.c_str(),"CEN_max_tracker")) {curstring >> varname >> value; CEN_max_tracker = value;}
	577	else if(strstr(temp_string.c_str(),"CEN_max_calo_cen")) {curstring >> varname >> value; CEN_max_calo_cen = value;}
[494]	578	else if(strstr(temp_string.c_str(),"CEN_max_calo_ec")) {CEN_max_calo_ec_flag=true; curstring >> varname >> value; CEN_max_calo_ec = value;}
[94]	579	else if(strstr(temp_string.c_str(),"CEN_max_calo_fwd")) {curstring >> varname >> value; CEN_max_calo_fwd = value;}
	580	else if(strstr(temp_string.c_str(),"CEN_max_mu")) {curstring >> varname >> value; CEN_max_mu = value;}
[494]	581
[94]	582	else if(strstr(temp_string.c_str(),"VFD_min_calo_vfd")) {curstring >> varname >> value; VFD_min_calo_vfd = value;}
	583	else if(strstr(temp_string.c_str(),"VFD_max_calo_vfd")) {curstring >> varname >> value; VFD_max_calo_vfd = value;}
	584	else if(strstr(temp_string.c_str(),"VFD_min_zdc")) {curstring >> varname >> value; VFD_min_zdc = value;}
	585	else if(strstr(temp_string.c_str(),"VFD_s_zdc")) {curstring >> varname >> value; VFD_s_zdc = value;}
	586
	587	else if(strstr(temp_string.c_str(),"RP_220_s")) {curstring >> varname >> value; RP_220_s = value;}
	588	else if(strstr(temp_string.c_str(),"RP_220_x")) {curstring >> varname >> value; RP_220_x = value;}
	589	else if(strstr(temp_string.c_str(),"RP_420_s")) {curstring >> varname >> value; RP_420_s = value;}
	590	else if(strstr(temp_string.c_str(),"RP_420_x")) {curstring >> varname >> value; RP_420_x = value;}
[257]	591	else if(strstr(temp_string.c_str(),"RP_beam1Card")) {curstring >> varname >> svalue;RP_beam1Card = svalue;}
	592	else if(strstr(temp_string.c_str(),"RP_beam2Card")) {curstring >> varname >> svalue;RP_beam2Card = svalue;}
	593	else if(strstr(temp_string.c_str(),"RP_IP_name")) {curstring >> varname >> svalue;RP_IP_name = svalue;}
[399]	594
	595	else if(strstr(temp_string.c_str(),"RP_offsetEl_s")) {curstring >> varname >> value; RP_offsetEl_s = value;}
	596	else if(strstr(temp_string.c_str(),"RP_offsetEl_x")) {curstring >> varname >> value; RP_offsetEl_x = value;}
[404]	597	else if(strstr(temp_string.c_str(),"RP_offsetEl_y")) {curstring >> varname >> value; RP_offsetEl_y = value;}
[399]	598	else if(strstr(temp_string.c_str(),"RP_cross_x")) {curstring >> varname >> value; RP_cross_x = value;}
	599	else if(strstr(temp_string.c_str(),"RP_cross_y")) {curstring >> varname >> value; RP_cross_y = value;}
	600	else if(strstr(temp_string.c_str(),"RP_cross_ang_x")) {curstring >> varname >> value; RP_cross_ang_x = value;}
	601	else if(strstr(temp_string.c_str(),"RP_cross_ang_y")) {curstring >> varname >> value; RP_cross_ang_y = value;}
[94]	602
	603	else if(strstr(temp_string.c_str(),"ELG_Scen")) {curstring >> varname >> value; ELG_Scen = value;}
	604	else if(strstr(temp_string.c_str(),"ELG_Ncen")) {curstring >> varname >> value; ELG_Ncen = value;}
	605	else if(strstr(temp_string.c_str(),"ELG_Ccen")) {curstring >> varname >> value; ELG_Ccen = value;}
[494]	606	else if(strstr(temp_string.c_str(),"ELG_Sec")) {curstring >> varname >> value; ELG_Sec = value;}
	607	else if(strstr(temp_string.c_str(),"ELG_Nec")) {curstring >> varname >> value; ELG_Nec = value;}
	608	else if(strstr(temp_string.c_str(),"ELG_Cec")) {curstring >> varname >> value; ELG_Cec = value;}
[94]	609	else if(strstr(temp_string.c_str(),"ELG_Sfwd")) {curstring >> varname >> value; ELG_Sfwd = value;}
	610	else if(strstr(temp_string.c_str(),"ELG_Cfwd")) {curstring >> varname >> value; ELG_Cfwd = value;}
	611	else if(strstr(temp_string.c_str(),"ELG_Nfwd")) {curstring >> varname >> value; ELG_Nfwd = value;}
[374]	612	else if(strstr(temp_string.c_str(),"ELG_Szdc")) {curstring >> varname >> value; ELG_Szdc = value;}
	613	else if(strstr(temp_string.c_str(),"ELG_Czdc")) {curstring >> varname >> value; ELG_Czdc = value;}
	614	else if(strstr(temp_string.c_str(),"ELG_Nzdc")) {curstring >> varname >> value; ELG_Nzdc = value;}
	615
[494]	616	else if(strstr(temp_string.c_str(),"HAD_Shcal")) {warning("HAD_Shcal","HAD_Scen"); curstring >> varname >> value; HAD_Scen = value;}
	617	else if(strstr(temp_string.c_str(),"HAD_Nhcal")) {warning("HAD_Nhcal","HAD_Ncen"); curstring >> varname >> value; HAD_Ncen = value;}
	618	else if(strstr(temp_string.c_str(),"HAD_Chcal")) {warning("HAD_Chcal","HAD_Ccen"); curstring >> varname >> value; HAD_Ccen = value;}
	619	else if(strstr(temp_string.c_str(),"HAD_Shf")) {warning("HAD_Shf","HAD_Sfwd"); curstring >> varname >> value; HAD_Sfwd = value;}
	620	else if(strstr(temp_string.c_str(),"HAD_Nhf")) {warning("HAD_Nhf","HAD_Nfwd"); curstring >> varname >> value; HAD_Nfwd = value;}
	621	else if(strstr(temp_string.c_str(),"HAD_Chf")) {warning("HAD_Chf","HAD_Cfwd"); curstring >> varname >> value; HAD_Cfwd = value;}
	622
	623
	624
	625	else if(strstr(temp_string.c_str(),"HAD_Scen")) {curstring >> varname >> value; HAD_Scen = value;}
	626	else if(strstr(temp_string.c_str(),"HAD_Ncen")) {curstring >> varname >> value; HAD_Ncen = value;}
	627	else if(strstr(temp_string.c_str(),"HAD_Ccen")) {curstring >> varname >> value; HAD_Ccen = value;}
	628	else if(strstr(temp_string.c_str(),"HAD_Sec")) {curstring >> varname >> value; HAD_Sec = value;}
	629	else if(strstr(temp_string.c_str(),"HAD_Nec")) {curstring >> varname >> value; HAD_Nec = value;}
	630	else if(strstr(temp_string.c_str(),"HAD_Cec")) {curstring >> varname >> value; HAD_Cec = value;}
	631	else if(strstr(temp_string.c_str(),"HAD_Sfwd")) {curstring >> varname >> value; HAD_Sfwd = value;}
	632	else if(strstr(temp_string.c_str(),"HAD_Nfwd")) {curstring >> varname >> value; HAD_Nfwd = value;}
	633	else if(strstr(temp_string.c_str(),"HAD_Cfwd")) {curstring >> varname >> value; HAD_Cfwd = value;}
[374]	634	else if(strstr(temp_string.c_str(),"HAD_Szdc")) {curstring >> varname >> value; HAD_Szdc = value;}
	635	else if(strstr(temp_string.c_str(),"HAD_Nzdc")) {curstring >> varname >> value; HAD_Nzdc = value;}
	636	else if(strstr(temp_string.c_str(),"HAD_Czdc")) {curstring >> varname >> value; HAD_Czdc = value;}
[494]	637
[374]	638	else if(strstr(temp_string.c_str(),"ZDC_T_resolution")) {curstring >> varname >> value; ZDC_T_resolution = value;}
	639	else if(strstr(temp_string.c_str(),"RP220_T_resolution")) {curstring >> varname >> value; RP220_T_resolution = value;}
	640	else if(strstr(temp_string.c_str(),"RP420_T_resolution")) {curstring >> varname >> value; RP420_T_resolution = value;}
[94]	641	else if(strstr(temp_string.c_str(),"MU_SmearPt")) {curstring >> varname >> value; MU_SmearPt = value;}
	642
	643	else if(strstr(temp_string.c_str(),"TRACK_radius")) {curstring >> varname >> ivalue;TRACK_radius = ivalue;}
	644	else if(strstr(temp_string.c_str(),"TRACK_length")) {curstring >> varname >> ivalue;TRACK_length = ivalue;}
	645	else if(strstr(temp_string.c_str(),"TRACK_bfield_x")) {curstring >> varname >> value; TRACK_bfield_x = value;}
	646	else if(strstr(temp_string.c_str(),"TRACK_bfield_y")) {curstring >> varname >> value; TRACK_bfield_y = value;}
	647	else if(strstr(temp_string.c_str(),"TRACK_bfield_z")) {curstring >> varname >> value; TRACK_bfield_z = value;}
	648	else if(strstr(temp_string.c_str(),"FLAG_bfield")) {curstring >> varname >> ivalue; FLAG_bfield = ivalue;}
	649	else if(strstr(temp_string.c_str(),"TRACK_ptmin")) {curstring >> varname >> value; TRACK_ptmin = value;}
[531]	650	else if(strstr(temp_string.c_str(),"TRACK_eff")) {curstring >> varname >> value; TRACK_eff = value;}
[33]	651
[94]	652	else if(strstr(temp_string.c_str(),"TOWER_number")) {curstring >> varname >> ivalue;TOWER_number = ivalue;}
	653	else if(strstr(temp_string.c_str(),"TOWER_eta_edges")){
	654	curstring >> varname; for(unsigned int i=0; i<TOWER_number+1; i++) {curstring >> value; TOWER_eta_edges[i] = value;} }
	655	else if(strstr(temp_string.c_str(),"TOWER_dphi")){
	656	curstring >> varname; for(unsigned int i=0; i<TOWER_number; i++) {curstring >> value; TOWER_dphi[i] = value;} }
[2]	657
[94]	658	else if(strstr(temp_string.c_str(),"PTCUT_elec")) {curstring >> varname >> value; PTCUT_elec = value;}
	659	else if(strstr(temp_string.c_str(),"PTCUT_muon")) {curstring >> varname >> value; PTCUT_muon = value;}
	660	else if(strstr(temp_string.c_str(),"PTCUT_jet")) {curstring >> varname >> value; PTCUT_jet = value;}
	661	else if(strstr(temp_string.c_str(),"PTCUT_gamma")) {curstring >> varname >> value; PTCUT_gamma = value;}
	662	else if(strstr(temp_string.c_str(),"PTCUT_taujet")) {curstring >> varname >> value; PTCUT_taujet = value;}
[374]	663	else if(strstr(temp_string.c_str(),"ZDC_gamma_E")) {curstring >> varname >> value; ZDC_gamma_E = value;}
	664	else if(strstr(temp_string.c_str(),"ZDC_n_E")) {curstring >> varname >> value; ZDC_n_E = value;}
[43]	665
[557]	666	else if(strstr(temp_string.c_str(),"ISOL_trk_PT")) {curstring >> varname >> value; ISOL_trk_PT = value;}
	667	else if(strstr(temp_string.c_str(),"ISOL_trk_Cone")) {curstring >> varname >> value; ISOL_trk_Cone = value;}
	668	else if(strstr(temp_string.c_str(),"ISOL_calo_Cone")) {curstring >> varname >> value; ISOL_calo_Cone = value;}
	669	else if(strstr(temp_string.c_str(),"ISOL_calo_ET")) {curstring >> varname >> value; ISOL_calo_ET = value;}
	670	else if(strstr(temp_string.c_str(),"ISOL_calo_Grid")) {curstring >> varname >> ivalue; ISOL_calo_Grid = ivalue;}
[305]	671
[94]	672	else if(strstr(temp_string.c_str(),"JET_coneradius")) {curstring >> varname >> value; JET_coneradius = value;}
	673	else if(strstr(temp_string.c_str(),"JET_jetalgo")) {curstring >> varname >> ivalue;JET_jetalgo = ivalue;}
	674	else if(strstr(temp_string.c_str(),"JET_seed")) {curstring >> varname >> value; JET_seed = value;}
[384]	675	else if(strstr(temp_string.c_str(),"JET_Eflow")) {curstring >> varname >> ivalue; JET_Eflow = ivalue;}
[94]	676
[558]	677	else if(strstr(temp_string.c_str(),"BTAG_func_b")) {curstring >> varname >> svalue; BTAG_func_b = svalue;}
	678	else if(strstr(temp_string.c_str(),"BTAG_func_c")) {curstring >> varname >> svalue; BTAG_func_c = svalue;}
	679	else if(strstr(temp_string.c_str(),"BTAG_func_l")) {curstring >> varname >> svalue; BTAG_func_l = svalue;}
	680	else if(strstr(temp_string.c_str(),"BTAG_func_g")) {curstring >> varname >> svalue; BTAG_func_g = svalue;}
[2]	681
[94]	682	else if(strstr(temp_string.c_str(),"FLAG_vfd")) {curstring >> varname >> ivalue; FLAG_vfd = ivalue;}
[526]	683	else if(strstr(temp_string.c_str(),"FLAG_RP")) {curstring >> varname >> ivalue; FLAG_RP = ivalue;}
[94]	684	else if(strstr(temp_string.c_str(),"FLAG_trigger")) {curstring >> varname >> ivalue; FLAG_trigger = ivalue;}
	685	else if(strstr(temp_string.c_str(),"FLAG_frog")) {curstring >> varname >> ivalue; FLAG_frog = ivalue;}
[307]	686	else if(strstr(temp_string.c_str(),"FLAG_lhco")) {curstring >> varname >> ivalue; FLAG_lhco = ivalue;}
[94]	687	else if(strstr(temp_string.c_str(),"NEvents_Frog")) {curstring >> varname >> ivalue; NEvents_Frog = ivalue;}
[422]	688	else if(strstr(temp_string.c_str(),"NEvents")) {curstring >> varname >> ivalue; NEvents = ivalue;}
[380]	689
	690	else if(strstr(temp_string.c_str(),"PdgTableFilename")) {curstring >> varname >> svalue; PdgTableFilename = svalue;}
[94]	691	}
[494]	692
	693	// for compatibility with old data cards
	694	if(!CEN_max_calo_ec_flag) {
	695	cout << " Warning \'CEN_max_calo_ec\' not found in datacard. "<< endl;
	696	cout << " Same values will be applied for calorimeter endcaps "<< endl;
	697	cout << " as for central calorimeters "<< endl;
	698	string text = "** Please update your card ("+ datacard +")";
	699	cout << left << setw(67) << text << right << setw(2) << "**" << endl;
	700	cout << " This change is 100\% backward compatibly with older DetectorCard. " << endl;
	701	cout << " However, please update your DetectorCard. " << endl;
	702	CEN_max_calo_ec = CEN_max_calo_cen;
	703	CEN_max_calo_cen = CEN_max_calo_cen/2;
	704	ELG_Sec = ELG_Scen;
	705	ELG_Nec = ELG_Ncen;
	706	ELG_Cec = ELG_Ccen;
	707	HAD_Sec = HAD_Scen;
	708	HAD_Nec = HAD_Ncen;
	709	HAD_Cec = HAD_Ccen;
	710	}
[392]	711
[557]	712	if(ISOL_calo_Grid%2 ==0) {
	713	ISOL_calo_Grid++;
	714	cout <<" WARNING: ISOL_Calo_Grid is not odd. Set it to "<< ISOL_calo_Grid << " " << endl;
[392]	715	}
	716
[94]	717	//jet stuffs not defined in the input datacard
	718	JET_overlap = 0.75;
	719	// MidPoint algorithm definition
	720	JET_M_coneareafraction = 0.25;
	721	JET_M_maxpairsize = 2;
	722	JET_M_maxiterations = 100;
	723	// Define Cone algorithm.
	724	JET_C_adjacencycut = 2;
	725	JET_C_maxiterations = 100;
	726	JET_C_iratch = 1;
	727	//Define SISCone algorithm.
	728	JET_S_npass = 0;
	729	JET_S_protojet_ptmin= 0.0;
	730
	731	//For Tau-jet definition
	732	TAU_energy_scone = 0.15; // radius R of the cone for tau definition, based on energy threshold
	733	TAU_track_scone = 0.4; // radius R of the cone for tau definition, based on track number
	734	TAU_track_pt = 2; // minimal pt [GeV] for tracks to be considered in tau definition
	735	TAU_energy_frac = 0.95; // fraction of energy required in the central part of the cone, for tau jets
	736
[2]	737	}
	738
[219]	739	void RESOLution::Logfile(const string& LogName) {
[454]	740
	741	// creates the list of good input files
	742	// this list is vector<string> inputfiles.
	743	ifstream infile(inputfilelist.c_str());
	744	vector<string> inputfiles;
	745	string filename;
	746	while(1) {
	747	infile >> filename; // reads the first line of the list
	748	if(!infile.good()) break; // quits when at the end of the list
	749	ifstream checking_the_file(filename.c_str()); // try to open the file
	750	if(!checking_the_file.good()) continue; // skips bad/unknown files
	751	else checking_the_file.close(); // close file if found
	752	inputfiles.push_back(filename); // append the name to the vector
	753	}
	754	infile.close();
	755
[44]	756	ofstream f_out(LogName.c_str());
[260]	757
	758	f_out <<"**********************************************************************"<< endl;
	759	f_out <<"**********************************************************************"<< endl;
	760	f_out <<" "<< endl;
	761	f_out <<" Welcome to "<< endl;
	762	f_out <<" "<< endl;
	763	f_out <<" "<< endl;
	764	f_out <<" .ddddddd- lL hH "<< endl;
	765	f_out <<" -Dd` `dD: Ll hH` "<< endl;
	766	f_out <<" dDd dDd eeee. lL .pp+pp Hh+hhh` -eeee- `sssss "<< endl;
	767	f_out <<" -Dd `DD ee. ee Ll .Pp. PP Hh. HH. ee. ee sSs "<< endl;
	768	f_out <<" dD` dDd eEeee: lL. pP. pP hH hH` eEeee:` -sSSSs. "<< endl;
	769	f_out <<" .Dd :dd eE. LlL PpppPP Hh Hh eE sSS "<< endl;
	770	f_out <<" dddddd:. eee+: lL. pp. hh. hh eee+ sssssS "<< endl;
	771	f_out <<" Pp "<< endl;
	772	f_out <<" "<< endl;
	773	f_out <<" Delphes, a framework for the fast simulation "<< endl;
	774	f_out <<" of a generic collider experiment "<< endl;
	775	f_out <<" "<< endl;
[567]	776	f_out <<" --- Version 1.9 of Delphes --- "<< endl;
	777	f_out <<" Last date of change: 7 May 2010 "<< endl;
[260]	778	f_out <<" "<< endl;
	779	f_out <<" "<< endl;
	780	f_out <<" This package uses: "<< endl;
	781	f_out <<" ------------------ "<< endl;
	782	f_out <<" FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210] "<< endl;
	783	f_out <<" Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2] "<< endl;
	784	f_out <<" FROG: L. Quertenmont, V. Roberfroid [hep-ex/0901.2718v1] "<< endl;
	785	f_out <<" "<< endl;
	786	f_out <<" ---------------------------------------------------------------- "<< endl;
	787	f_out <<" "<< endl;
	788	f_out <<" Main authors: "<< endl;
	789	f_out <<" ------------- "<< endl;
	790	f_out <<" "<< endl;
	791	f_out <<" SÃ©verine Ovyn Xavier Rouby "<< endl;
	792	f_out <<" severine.ovyn@uclouvain.be xavier.rouby@cern "<< endl;
	793	f_out <<" Center for Particle Physics and Phenomenology (CP3) "<< endl;
	794	f_out <<" Universite Catholique de Louvain (UCL) "<< endl;
	795	f_out <<" Louvain-la-Neuve, Belgium "<< endl;
	796	f_out <<" "<< endl;
	797	f_out <<" ---------------------------------------------------------------- "<< endl;
	798	f_out <<" "<< endl;
	799	f_out <<" Former Delphes versions and documentation can be found on : "<< endl;
	800	f_out <<" http://www.fynu.ucl.ac.be/delphes.html "<< endl;
	801	f_out <<" "<< endl;
	802	f_out <<" "<< endl;
[528]	803	f_out <<" Disclaimer: this program comes without guarantees. "<< endl;
	804	f_out <<" Beware of errors and please give us "<< endl;
	805	f_out <<" your feedbacks about potential bugs. "<< endl;
[260]	806	f_out <<" "<< endl;
	807	f_out <<"**********************************************************************"<< endl;
	808	f_out <<" "<< endl;
[380]	809	f_out<<"#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"<<"\n";
[454]	810	f_out <<"* *"<< endl;
	811	f_out <<"#******************************** *"<<"\n";
	812	f_out <<"# Input files *"<<"\n";
	813	f_out <<"#******************************** *"<<"\n";
	814	f_out << left << setw(22) <<"* Input list "<<""
	815	<< left << setw(39) << inputfilelist << "" << right << setw(9) << "*"<<"\n";
	816	for (unsigned int i =0; i<inputfiles.size(); i++) {
	817	f_out << left << setw(22) <<"* - file "<<""
	818	<< left << setw(43) << inputfiles[i] << "" << right << setw(5) << "*"<<"\n";
	819	}
	820	if(detectorcard != "")
	821	f_out << left << setw(22) <<"* Detector card "<<""
	822	<< left << setw(39) << detectorcard << "" << right << setw(9) << "*"<<"\n";
	823	if(triggercard != "")
	824	f_out << left << setw(22) <<"* Trigger card "<<""
	825	<< left << setw(39) << triggercard << "" << right << setw(9) << "*"<<"\n";
	826	f_out<<"* Beam optics : *"<<"\n";
	827	f_out << left << setw(22) <<"* - beam 1 "<<""
	828	<< left << setw(33) << RP_beam1Card << "" << right << setw(15) << "*"<<"\n";
	829	f_out << left << setw(22) <<"* - beam 2 "<<""
	830	<< left << setw(33) << RP_beam2Card << "" << right << setw(15) << "*"<<"\n";
	831	f_out << left << setw(22) <<"* Input PDG table " << ""
	832	<< left << setw(39) << PdgTableFilename << "" << right << setw(9) << "*"<<"\n";
	833
[44]	834	f_out<<"* *"<<"\n";
[380]	835	f_out<<"* *"<<"\n";
[44]	836	f_out<<"#******************************** *"<<"\n";
	837	f_out<<"# Central detector caracteristics *"<<"\n";
	838	f_out<<"#******************************** *"<<"\n";
	839	f_out<<"* *"<<"\n";
	840	f_out << left << setw(30) <<"* Maximum tracking system: "<<""
[94]	841	<< left << setw(10) <<CEN_max_tracker <<""<< right << setw(15)<<"*"<<"\n";
[44]	842	f_out << left << setw(30) <<"* Maximum central calorimeter: "<<""
[94]	843	<< left << setw(10) <<CEN_max_calo_cen <<""<< right << setw(15)<<"*"<<"\n";
[494]	844	f_out << left << setw(30) <<"* Maximum endcap calorimeter: "<<""
	845	<< left << setw(10) <<CEN_max_calo_ec <<""<< right << setw(15)<<"*"<<"\n";
	846	f_out << left << setw(30) <<"* Maximum central calorimeter: "<<""
[94]	847	<< left << setw(10) <<CEN_max_calo_fwd <<""<< right << setw(15)<<"*"<<"\n";
[44]	848	f_out << left << setw(30) <<"* Muon chambers coverage: "<<""
[94]	849	<< left << setw(10) <<CEN_max_mu <<""<< right << setw(15)<<"*"<<"\n";
[44]	850	f_out<<"* *"<<"\n";
[306]	851	if(FLAG_RP==1){
	852	f_out<<"#************************************ *"<<"\n";
	853	f_out<<"# Very forward Roman Pots switched on *"<<"\n";
	854	f_out<<"#************************************ *"<<"\n";
[94]	855	f_out<<"* *"<<"\n";
[306]	856	f_out << left << setw(55) <<"* Distance of the 220 RP to the IP in meters:"<<""
[94]	857	<< left << setw(5) <<RP_220_s <<""<< right << setw(10)<<"*"<<"\n";
[306]	858	f_out << left << setw(55) <<"* Distance of the 220 RP to the beam in meters:"<<""
[94]	859	<< left << setw(5) <<RP_220_x <<""<< right << setw(10)<<"*"<<"\n";
[306]	860	f_out << left << setw(55) <<"* Distance of the 420 RP to the IP in meters:"<<""
[94]	861	<< left << setw(5) <<RP_420_s <<""<< right << setw(10)<<"*"<<"\n";
[306]	862	f_out << left << setw(55) <<"* Distance of the 420 RP to the beam in meters:"<<""
[94]	863	<< left << setw(5) <<RP_420_x <<""<< right << setw(10)<<"*"<<"\n";
[257]	864	f_out << left << setw(55) <<"* Interaction point at the LHC named: "<<""
	865	<< left << setw(5) <<RP_IP_name <<""<< right << setw(10)<<"*"<<"\n";
[252]	866	f_out << left << setw(35) <<"* Datacard for beam 1: "<<""
	867	<< left << setw(25) <<RP_beam1Card <<""<< right << setw(10)<<"*"<<"\n";
	868	f_out << left << setw(35) <<"* Datacard for beam 2: "<<""
	869	<< left << setw(25) <<RP_beam2Card <<""<< right << setw(10)<<"*"<<"\n";
[494]	870	f_out << left << setw(54) <<"* Beam separation, in meters(hor):"<<""
[399]	871	<< left << setw(6) << RP_offsetEl_x <<""<< right << setw(10)<<"*"<<"\n";
[494]	872	f_out << left << setw(54) <<"* Beam separation, in meters(ver):"<<""
[404]	873	<< left << setw(6) << RP_offsetEl_y <<""<< right << setw(10)<<"*"<<"\n";
[494]	874	f_out << left << setw(54) <<"* Distance from IP for Beam separation (m):"<<""
[399]	875	<< left << setw(6) <<RP_offsetEl_s <<""<< right << setw(10)<<"*"<<"\n";
[494]	876	f_out << left << setw(54) <<"* X offset of beam crossing in micrometers:"<<""
[399]	877	<< left << setw(6) <<RP_cross_x <<""<< right << setw(10)<<"*"<<"\n";
[494]	878	f_out << left << setw(54) <<"* Y offset of beam crossing in micrometers:"<<""
[399]	879	<< left << setw(6) <<RP_cross_y <<""<< right << setw(10)<<"*"<<"\n";
[494]	880	f_out << left << setw(54) <<"* X Angle of beam crossing:"<<""
[399]	881	<< left << setw(6) <<RP_cross_ang_x <<""<< right << setw(10)<<"*"<<"\n";
[494]	882	f_out << left << setw(54) <<"* Y Angle of beam crossing:"<<""
[399]	883	<< left << setw(6) <<RP_cross_ang_y <<""<< right << setw(10)<<"*"<<"\n";
[94]	884	f_out<<"* *"<<"\n";
	885	}
	886	else {
[306]	887	f_out<<"#************************************* *"<<"\n";
	888	f_out<<"# Very forward Roman Pots switched off *"<<"\n";
	889	f_out<<"#************************************* *"<<"\n";
[94]	890	f_out<<"* *"<<"\n";
	891	}
[306]	892	if(FLAG_vfd==1){
	893	f_out<<"#************************************** *"<<"\n";
	894	f_out<<"# Very forward calorimeters switched on *"<<"\n";
	895	f_out<<"#************************************** *"<<"\n";
	896	f_out<<"* *"<<"\n";
	897	f_out << left << setw(55) <<"* Minimum very forward calorimeter: "<<""
	898	<< left << setw(5) <<VFD_min_calo_vfd <<""<< right << setw(10)<<"*"<<"\n";
	899	f_out << left << setw(55) <<"* Maximum very forward calorimeter: "<<""
	900	<< left << setw(5) <<VFD_max_calo_vfd <<""<< right << setw(10)<<"*"<<"\n";
	901	f_out << left << setw(55) <<"* Minimum coverage zero_degree calorimeter "<<""
	902	<< left << setw(5) <<VFD_min_zdc <<""<< right << setw(10)<<"*"<<"\n";
	903	f_out << left << setw(55) <<"* Distance of the ZDC to the IP, in meters: "<<""
	904	<< left << setw(5) <<VFD_s_zdc <<""<< right << setw(10)<<"*"<<"\n";
	905	f_out<<"* *"<<"\n";
	906	}
	907	else {
	908	f_out<<"#*************************************** *"<<"\n";
	909	f_out<<"# Very forward calorimeters switched off *"<<"\n";
	910	f_out<<"#*************************************** *"<<"\n";
	911	f_out<<"* *"<<"\n";
	912	}
	913
[44]	914	f_out<<"#************************************ *"<<"\n";
	915	f_out<<"# Electromagnetic smearing parameters *"<<"\n";
	916	f_out<<"#************************************ *"<<"\n";
	917	f_out<<"* *"<<"\n";
	918	//# \sigma/E = C + N/E + S/\sqrt{E}
	919	f_out << left << setw(30) <<"* S term for central ECAL: "<<""
	920	<< left << setw(30) <<ELG_Scen <<""<< right << setw(10)<<"*"<<"\n";
	921	f_out << left << setw(30) <<"* N term for central ECAL: "<<""
	922	<< left << setw(30) <<ELG_Ncen <<""<< right << setw(10)<<"*"<<"\n";
	923	f_out << left << setw(30) <<"* C term for central ECAL: "<<""
	924	<< left << setw(30) <<ELG_Ccen <<""<< right << setw(10)<<"*"<<"\n";
[494]	925	f_out << left << setw(30) <<"* S term for ECAL end-cap: "<<""
	926	<< left << setw(30) <<ELG_Sec <<""<< right << setw(10)<<"*"<<"\n";
	927	f_out << left << setw(30) <<"* N term for ECAL end-cap: "<<""
	928	<< left << setw(30) <<ELG_Nec <<""<< right << setw(10)<<"*"<<"\n";
	929	f_out << left << setw(30) <<"* C term for ECAL end-cap: "<<""
	930	<< left << setw(30) <<ELG_Cec <<""<< right << setw(10)<<"*"<<"\n";
	931
	932
[257]	933	f_out << left << setw(30) <<"* S term for FCAL: "<<""
[44]	934	<< left << setw(30) <<ELG_Sfwd <<""<< right << setw(10)<<"*"<<"\n";
[257]	935	f_out << left << setw(30) <<"* N term for FCAL: "<<""
[44]	936	<< left << setw(30) <<ELG_Nfwd <<""<< right << setw(10)<<"*"<<"\n";
[257]	937	f_out << left << setw(30) <<"* C term for FCAL: "<<""
[44]	938	<< left << setw(30) <<ELG_Cfwd <<""<< right << setw(10)<<"*"<<"\n";
[374]	939	f_out << left << setw(30) <<"* S term for ZDC: "<<""
	940	<< left << setw(30) <<ELG_Szdc <<""<< right << setw(10)<<"*"<<"\n";
	941	f_out << left << setw(30) <<"* N term for ZDC: "<<""
	942	<< left << setw(30) <<ELG_Nzdc <<""<< right << setw(10)<<"*"<<"\n";
	943	f_out << left << setw(30) <<"* C term for ZDC: "<<""
	944	<< left << setw(30) <<ELG_Czdc <<""<< right << setw(10)<<"*"<<"\n";
	945
[44]	946	f_out<<"* *"<<"\n";
	947	f_out<<"#***************************** *"<<"\n";
	948	f_out<<"# Hadronic smearing parameters *"<<"\n";
	949	f_out<<"#***************************** *"<<"\n";
	950	f_out<<"* *"<<"\n";
	951	f_out << left << setw(30) <<"* S term for central HCAL: "<<""
[494]	952	<< left << setw(30) <<HAD_Scen <<""<< right << setw(10)<<"*"<<"\n";
[44]	953	f_out << left << setw(30) <<"* N term for central HCAL: "<<""
[494]	954	<< left << setw(30) <<HAD_Ncen <<""<< right << setw(10)<<"*"<<"\n";
[44]	955	f_out << left << setw(30) <<"* C term for central HCAL: "<<""
[494]	956	<< left << setw(30) <<HAD_Ccen <<""<< right << setw(10)<<"*"<<"\n";
	957	f_out << left << setw(30) <<"* S term for HCAL endcap: "<<""
	958	<< left << setw(30) <<HAD_Sec <<""<< right << setw(10)<<"*"<<"\n";
	959	f_out << left << setw(30) <<"* N term for HCAL endcap: "<<""
	960	<< left << setw(30) <<HAD_Nec <<""<< right << setw(10)<<"*"<<"\n";
	961	f_out << left << setw(30) <<"* C term for HCAL endcap: "<<""
	962	<< left << setw(30) <<HAD_Cec <<""<< right << setw(10)<<"*"<<"\n";
[257]	963	f_out << left << setw(30) <<"* S term for FCAL: "<<""
[494]	964	<< left << setw(30) <<HAD_Sfwd <<""<< right << setw(10)<<"*"<<"\n";
[257]	965	f_out << left << setw(30) <<"* N term for FCAL: "<<""
[494]	966	<< left << setw(30) <<HAD_Nfwd <<""<< right << setw(10)<<"*"<<"\n";
[257]	967	f_out << left << setw(30) <<"* C term for FCAL: "<<""
[494]	968	<< left << setw(30) <<HAD_Cfwd <<""<< right << setw(10)<<"*"<<"\n";
[374]	969	f_out << left << setw(30) <<"* S term for ZDC: "<<""
	970	<< left << setw(30) <<HAD_Szdc <<""<< right << setw(10)<<"*"<<"\n";
	971	f_out << left << setw(30) <<"* N term for ZDC: "<<""
	972	<< left << setw(30) <<HAD_Nzdc <<""<< right << setw(10)<<"*"<<"\n";
	973	f_out << left << setw(30) <<"* C term for ZDC: "<<""
	974	<< left << setw(30) <<HAD_Czdc <<""<< right << setw(10)<<"*"<<"\n";
	975
[44]	976	f_out<<"* *"<<"\n";
	977	f_out<<"#************************* *"<<"\n";
[374]	978	f_out<<"# Time smearing parameters *"<<"\n";
	979	f_out<<"#************************* *"<<"\n";
	980	f_out<<"* *"<<"\n";
	981	f_out << left << setw(55) <<"* Time resolution for ZDC : "<<""
	982	<< left << setw(5) <<ZDC_T_resolution <<""<< right << setw(10)<<"*"<<"\n";
	983	f_out << left << setw(55) <<"* Time resolution for RP220 : "<<""
	984	<< left << setw(5) <<RP220_T_resolution <<""<< right << setw(10)<<"*"<<"\n";
	985	f_out << left << setw(55) <<"* Time resolution for RP420 : "<<""
	986	<< left << setw(5) <<RP420_T_resolution <<""<< right << setw(10)<<"*"<<"\n";
	987	f_out<<"* *"<<"\n";
	988
	989	f_out<<"* *"<<"\n";
	990	f_out<<"#************************* *"<<"\n";
[44]	991	f_out<<"# Muon smearing parameters *"<<"\n";
	992	f_out<<"#************************* *"<<"\n";
	993	f_out<<"* *"<<"\n";
[94]	994	f_out << left << setw(55) <<"* PT resolution for muons : "<<""
	995	<< left << setw(5) <<MU_SmearPt <<""<< right << setw(10)<<"*"<<"\n";
[44]	996	f_out<<"* *"<<"\n";
[94]	997	if(FLAG_bfield==1){
	998	f_out<<"#*************************** *"<<"\n";
[264]	999	f_out<<"# Magnetic field switched on *"<<"\n";
[94]	1000	f_out<<"#*************************** *"<<"\n";
	1001	f_out<<"* *"<<"\n";
	1002	f_out << left << setw(55) <<"* Radius of the BField coverage: "<<""
	1003	<< left << setw(5) <<TRACK_radius <<""<< right << setw(10)<<"*"<<"\n";
	1004	f_out << left << setw(55) <<"* Length of the BField coverage: "<<""
	1005	<< left << setw(5) <<TRACK_length <<""<< right << setw(10)<<"*"<<"\n";
	1006	f_out << left << setw(55) <<"* BField X component: "<<""
	1007	<< left << setw(5) <<TRACK_bfield_x <<""<< right << setw(10)<<"*"<<"\n";
	1008	f_out << left << setw(55) <<"* BField Y component: "<<""
	1009	<< left << setw(5) <<TRACK_bfield_y <<""<< right << setw(10)<<"*"<<"\n";
	1010	f_out << left << setw(55) <<"* BField Z component: "<<""
	1011	<< left << setw(5) <<TRACK_bfield_z <<""<< right << setw(10)<<"*"<<"\n";
	1012	f_out << left << setw(55) <<"* Minimal pT needed to reach the calorimeter [GeV]: "<<""
	1013	<< left << setw(10) <<TRACK_ptmin <<""<< right << setw(5)<<"*"<<"\n";
	1014	f_out << left << setw(55) <<"* Efficiency associated to the tracking: "<<""
	1015	<< left << setw(10) <<TRACK_eff <<""<< right << setw(5)<<"*"<<"\n";
	1016	f_out<<"* *"<<"\n";
	1017	}
	1018	else {
	1019	f_out<<"#**************************** *"<<"\n";
[264]	1020	f_out<<"# Magnetic field switched off *"<<"\n";
[94]	1021	f_out<<"#**************************** *"<<"\n";
	1022	f_out << left << setw(55) <<"* Minimal pT needed to reach the calorimeter [GeV]: "<<""
	1023	<< left << setw(10) <<TRACK_ptmin <<""<< right << setw(5)<<"*"<<"\n";
	1024	f_out << left << setw(55) <<"* Efficiency associated to the tracking: "<<""
	1025	<< left << setw(10) <<TRACK_eff <<""<< right << setw(5)<<"*"<<"\n";
	1026	f_out<<"* *"<<"\n";
	1027	}
	1028	f_out<<"#******************** *"<<"\n";
	1029	f_out<<"# Calorimetric Towers *"<<"\n";
	1030	f_out<<"#******************** *"<<"\n";
	1031	f_out << left << setw(55) <<"* Number of calorimetric towers in eta, for eta>0: "<<""
	1032	<< left << setw(5) << TOWER_number <<""<< right << setw(10)<<"*"<<"\n";
	1033	f_out << left << setw(55) <<"* Tower edges in eta, for eta>0: "<<"" << right << setw(15)<<"*"<<"\n";
	1034	f_out << "* ";
	1035	for (unsigned int i=0; i<TOWER_number+1; i++) {
	1036	f_out << left << setw(7) << TOWER_eta_edges[i];
	1037	if(!( (i+1) %9 )) f_out << right << setw(3) << "" << "\n" << " ";
	1038	}
	1039	for (unsigned int i=(TOWER_number+1)%9; i<9; i++) f_out << left << setw(7) << "";
	1040	f_out << right << setw(3)<<"*"<<"\n";
	1041	f_out << left << setw(55) <<"* Tower sizes in phi, for eta>0 [degree]:"<<"" << right << setw(15)<<"*"<<"\n";
	1042	f_out << "* ";
	1043	for (unsigned int i=0; i<TOWER_number; i++) {
	1044	f_out << left << setw(7) << TOWER_dphi[i];
	1045	if(!( (i+1) %9 )) f_out << right << setw(3) << "" << "\n" << " ";
	1046	}
	1047	for (unsigned int i=(TOWER_number)%9; i<9; i++) f_out << left << setw(7) << "";
	1048	f_out << right << setw(3)<<"*"<<"\n";
[44]	1049	f_out<<"* *"<<"\n";
	1050	f_out<<"#******************* *"<<"\n";
	1051	f_out<<"# Minimum pT's [GeV] *"<<"\n";
	1052	f_out<<"#******************* *"<<"\n";
	1053	f_out<<"* *"<<"\n";
	1054	f_out << left << setw(40) <<"* Minimum pT for electrons: "<<""
[94]	1055	<< left << setw(20) <<PTCUT_elec <<""<< right << setw(10)<<"*"<<"\n";
[44]	1056	f_out << left << setw(40) <<"* Minimum pT for muons: "<<""
[94]	1057	<< left << setw(20) <<PTCUT_muon <<""<< right << setw(10)<<"*"<<"\n";
[44]	1058	f_out << left << setw(40) <<"* Minimum pT for jets: "<<""
[94]	1059	<< left << setw(20) <<PTCUT_jet <<""<< right << setw(10)<<"*"<<"\n";
[44]	1060	f_out << left << setw(40) <<"* Minimum pT for Tau-jets: "<<""
[94]	1061	<< left << setw(20) <<PTCUT_taujet <<""<< right << setw(10)<<"*"<<"\n";
[74]	1062	f_out << left << setw(40) <<"* Minimum pT for photons: "<<""
[94]	1063	<< left << setw(20) <<PTCUT_gamma <<""<< right << setw(10)<<"*"<<"\n";
[374]	1064	f_out << left << setw(40) <<"* Minimum E for photons in ZDC: "<<""
	1065	<< left << setw(20) <<ZDC_gamma_E <<""<< right << setw(10)<<"*"<<"\n";
	1066	f_out << left << setw(40) <<"* Minimum E for neutrons in ZDC: "<<""
	1067	<< left << setw(20) <<ZDC_n_E <<""<< right << setw(10)<<"*"<<"\n";
	1068
[44]	1069	f_out<<"* *"<<"\n";
[305]	1070	f_out<<"#******************* *"<<"\n";
	1071	f_out<<"# Isolation criteria *"<<"\n";
	1072	f_out<<"#******************* *"<<"\n";
	1073	f_out<<"* *"<<"\n";
[557]	1074	f_out << left << setw(40) <<"* Minimum pT for tracking isolation [GeV]: "<<""
	1075	<< left << setw(20) <<ISOL_trk_PT <<""<< right << setw(10)<<"*"<<"\n";
	1076	f_out << left << setw(40) <<"* Cone for tracking isolation criteria: "<<""
	1077	<< left << setw(20) <<ISOL_trk_Cone <<""<< right << setw(10)<<"*"<<"\n";
	1078	f_out << left << setw(40) <<"* Cone for calorimetric isolation criteria: "<<""
	1079	<< left << setw(20) <<ISOL_calo_Cone <<""<< right << setw(10)<<"*"<<"\n";
[321]	1080
[557]	1081
	1082	if(ISOL_calo_ET > 1E98) f_out<<"# No Calorimetric isolation applied *"<<"\n";
[321]	1083	else {
	1084	f_out << left << setw(40) <<"* Minimum ET for towers [GeV]: "<<""
[557]	1085	<< left << setw(20) <<ISOL_calo_ET <<""<< right << setw(10)<<"*"<<"\n";
[321]	1086	f_out << left << setw(40) <<"* Grid size (NxN) for calorimetric isolation: "<<""
[557]	1087	<< left << setw(20) <<ISOL_calo_Grid <<""<< right << setw(4)<<"*"<<"\n";
[321]	1088	}
	1089
	1090
[305]	1091	f_out<<"* *"<<"\n";
[44]	1092	f_out<<"#*************** *"<<"\n";
	1093	f_out<<"# Jet definition *"<<"\n";
	1094	f_out<<"#*************** *"<<"\n";
[383]	1095	if(JET_Eflow)
	1096	{
	1097	f_out<<"#*************** *"<<"\n";
	1098	f_out<<"#* Running considering perfect energy flow on the tracker coverage *"<<"\n";
	1099	}
	1100	else
	1101	{
	1102	f_out<<"#* Running considering no energy flow on the tracker coverage *"<<"\n";
	1103	f_out<<"#* --> jet algo applied on the calorimetric towers *"<<"\n";
	1104	}
[44]	1105	f_out<<"* *"<<"\n";
[49]	1106	f_out<<"* Six algorithms are currently available: *"<<"\n";
	1107	f_out<<"* - 1) CDF cone algorithm, *"<<"\n";
	1108	f_out<<"* - 2) CDF MidPoint algorithm, *"<<"\n";
	1109	f_out<<"* - 3) SIScone algorithm, *"<<"\n";
	1110	f_out<<"* - 4) kt algorithm, *"<<"\n";
	1111	f_out<<"* - 5) Cambrigde/Aachen algorithm, *"<<"\n";
	1112	f_out<<"* - 6) Anti-kt algorithm. *"<<"\n";
	1113	f_out<<"* *"<<"\n";
	1114	f_out<<"* You have chosen *"<<"\n";
[94]	1115	switch(JET_jetalgo) {
[44]	1116	default:
	1117	case 1: {
[94]	1118	f_out<<"* CDF JetClu jet algorithm with parameters: *"<<"\n";
	1119	f_out << left << setw(40) <<"* - Seed threshold: "<<""
	1120	<< left << setw(10) <<JET_seed <<""<< right << setw(20)<<"! not in datacard *"<<"\n";
	1121	f_out << left << setw(40) <<"* - Cone radius: "<<""
	1122	<< left << setw(10) <<JET_coneradius <<""<< right << setw(20)<<"*"<<"\n";
	1123	f_out << left << setw(40) <<"* - Adjacency cut: "<<""
	1124	<< left << setw(10) <<JET_C_adjacencycut <<""<< right << setw(20)<<"! not in datacard *"<<"\n";
	1125	f_out << left << setw(40) <<"* - Max iterations: "<<""
	1126	<< left << setw(10) <<JET_C_maxiterations <<""<< right << setw(20)<<"! not in datacard *"<<"\n";
	1127	f_out << left << setw(40) <<"* - Iratch: "<<""
	1128	<< left << setw(10) <<JET_C_iratch <<""<< right << setw(20)<<"! not in datacard *"<<"\n";
	1129	f_out << left << setw(40) <<"* - Overlap threshold: "<<""
	1130	<< left << setw(10) <<JET_overlap <<""<< right << setw(20)<<"! not in datacard *"<<"\n";
[44]	1131	}
	1132	break;
	1133	case 2: {
[94]	1134	f_out<<"* CDF midpoint jet algorithm with parameters: *"<<"\n";
	1135	f_out << left << setw(40) <<"* - Seed threshold: "<<""
	1136	<< left << setw(20) <<JET_seed <<""<< right << setw(10)<<"! not in datacard *"<<"\n";
	1137	f_out << left << setw(40) <<"* - Cone radius: "<<""
	1138	<< left << setw(20) <<JET_coneradius <<""<< right << setw(10)<<"*"<<"\n";
	1139	f_out << left << setw(40) <<"* - Cone area fraction:"<<""
	1140	<< left << setw(20) <<JET_M_coneareafraction <<""<< right << setw(10)<<"! not in datacard *"<<"\n";
	1141	f_out << left << setw(40) <<"* - Maximum pair size: "<<""
	1142	<< left << setw(20) <<JET_M_maxpairsize <<""<< right << setw(10)<<"! not in datacard *"<<"\n";
	1143	f_out << left << setw(40) <<"* - Max iterations: "<<""
	1144	<< left << setw(20) <<JET_M_maxiterations <<""<< right << setw(10)<<"! not in datacard *"<<"\n";
	1145	f_out << left << setw(40) <<"* - Overlap threshold: "<<""
	1146	<< left << setw(20) <<JET_overlap <<""<< right << setw(10)<<"! not in datacard *"<<"\n";
[44]	1147	}
	1148	break;
	1149	case 3: {
[94]	1150	f_out <<"* SISCone jet algorithm with parameters: *"<<"\n";
	1151	f_out << left << setw(40) <<"* - Cone radius: "<<""
	1152	<< left << setw(20) <<JET_coneradius <<""<< right << setw(10)<<"*"<<"\n";
	1153	f_out << left << setw(40) <<"* - Overlap threshold: "<<""
	1154	<< left << setw(20) <<JET_overlap <<""<< right << setw(10)<<"! not in datacard *"<<"\n";
	1155	f_out << left << setw(40) <<"* - Number pass max: "<<""
	1156	<< left << setw(20) <<JET_S_npass <<""<< right << setw(10)<<"! not in datacard *"<<"\n";
	1157	f_out << left << setw(40) <<"* - Minimum pT for protojet: "<<""
	1158	<< left << setw(20) <<JET_S_protojet_ptmin <<""<< right << setw(10)<<"! not in datacard *"<<"\n";
[44]	1159	}
	1160	break;
	1161	case 4: {
[94]	1162	f_out <<"* KT jet algorithm with parameters: *"<<"\n";
	1163	f_out << left << setw(40) <<"* - Cone radius: "<<""
	1164	<< left << setw(20) <<JET_coneradius <<""<< right << setw(10)<<"*"<<"\n";
[44]	1165	}
	1166	break;
[49]	1167	case 5: {
[94]	1168	f_out <<"* Cambridge/Aachen jet algorithm with parameters: *"<<"\n";
	1169	f_out << left << setw(40) <<"* - Cone radius: "<<""
	1170	<< left << setw(20) <<JET_coneradius <<""<< right << setw(10)<<"*"<<"\n";
[44]	1171	}
[49]	1172	break;
	1173	case 6: {
[94]	1174	f_out <<"* Anti-kt jet algorithm with parameters: *"<<"\n";
	1175	f_out << left << setw(40) <<"* - Cone radius: "<<""
	1176	<< left << setw(20) <<JET_coneradius <<""<< right << setw(10)<<"*"<<"\n";
[49]	1177	}
	1178	break;
	1179	}
[44]	1180	f_out<<"* *"<<"\n";
[94]	1181	f_out<<"#****************************** *"<<"\n";
	1182	f_out<<"# Tau-jet definition parameters *"<<"\n";
	1183	f_out<<"#****************************** *"<<"\n";
	1184	f_out<<"* *"<<"\n";
	1185	f_out << left << setw(45) <<"* Cone radius for calorimeter tagging: "<<""
	1186	<< left << setw(5) <<TAU_energy_scone <<""<< right << setw(20)<<"*"<<"\n";
	1187	f_out << left << setw(45) <<"* Fraction of energy in the small cone: "<<""
	1188	<< left << setw(5) <<TAU_energy_frac100 <<""<< right << setw(20)<<"! not in datacard "<<"\n";
	1189	f_out << left << setw(45) <<"* Cone radius for tracking tagging: "<<""
	1190	<< left << setw(5) <<TAU_track_scone <<""<< right << setw(20)<<"*"<<"\n";
	1191	f_out << left << setw(45) <<"* Minimum track pT [GeV]: "<<""
	1192	<< left << setw(5) <<TAU_track_pt <<""<< right << setw(20)<<"*"<<"\n";
	1193	f_out<<"* *"<<"\n";
[558]	1194	f_out<<"#******************************* *"<<"\n";
	1195	f_out<<"# B-tagging efficiency functions *"<<"\n";
	1196	f_out<<"#******************************* *"<<"\n";
[94]	1197	f_out<<"* *"<<"\n";
[558]	1198	f_out << left << setw(50) <<"* Tagging a \"b\" as a b-jet: "<<""
	1199	<< left << setw(18) <<BTAG_func_b <<""<< right << setw(2)<<"*"<<"\n";
	1200	f_out << left << setw(50) <<"* Mistagging a c-jet as a b-jet: "<<""
	1201	<< left << setw(18) <<BTAG_func_c <<""<< right << setw(2)<<"*"<<"\n";
	1202	f_out << left << setw(50) <<"* Mistagging a gluon-jet as a b-jet: "<<""
	1203	<< left << setw(18) <<BTAG_func_g <<""<< right << setw(2)<<"*"<<"\n";
	1204	f_out << left << setw(50) <<"* Mistagging a light jet as a b-jet: "<<""
	1205	<< left << setw(18) <<BTAG_func_l <<""<< right << setw(2)<<"*"<<"\n";
	1206
	1207
[94]	1208	f_out<<"* *"<<"\n";
	1209	f_out<<"* *"<<"\n";
[44]	1210	f_out<<"#....................................................................*"<<"\n";
	1211	f_out<<"#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"<<"\n";
[399]	1212
	1213	f_out.close();
[44]	1214	}
	1215
[2]	1216	// ********Provides the smeared TLorentzVector for the electrons******
	1217	// Smears the electron energy, and changes the 4-momentum accordingly
	1218	// different smearing if the electron is central (eta < 2.5) or forward
	1219	void RESOLution::SmearElectron(TLorentzVector &electron) {
	1220	// the 'electron' variable will be changed by the function
	1221	float energy = electron.E(); // before smearing
	1222	float energyS = 0.0; // after smearing // \sigma/E = C + N/E + S/\sqrt{E}
[71]	1223
[94]	1224	if(fabs(electron.Eta()) < CEN_max_tracker) { // if the electron is inside the tracker
[526]	1225	energyS = grandom->Gaus(energy, sqrt(
[2]	1226	pow(ELG_Ncen,2) +
	1227	pow(ELG_Ccen*energy,2) +
[22]	1228	pow(ELG_Scen*sqrt(energy),2) ));
[55]	1229	}
[494]	1230	else if(fabs(electron.Eta()) >= CEN_max_tracker && fabs(electron.Eta()) < CEN_max_calo_ec){
[526]	1231	energyS = grandom->Gaus(energy, sqrt(
[494]	1232	pow(ELG_Nec,2) +
	1233	pow(ELG_Cec*energy,2) +
	1234	pow(ELG_Sec*sqrt(energy),2) ) );
	1235	}
	1236	else if(fabs(electron.Eta()) >= CEN_max_calo_ec && fabs(electron.Eta()) < CEN_max_calo_fwd){
[526]	1237	energyS = grandom->Gaus(energy, sqrt(
[2]	1238	pow(ELG_Nfwd,2) +
	1239	pow(ELG_Cfwd*energy,2) +
	1240	pow(ELG_Sfwd*sqrt(energy),2) ) );
	1241	}
	1242	electron.SetPtEtaPhiE(energyS/cosh(electron.Eta()), electron.Eta(), electron.Phi(), energyS);
	1243	if(electron.E() < 0)electron.SetPxPyPzE(0,0,0,0); // no negative values after smearing !
	1244	}
	1245
	1246
	1247	// ********Provides the smeared TLorentzVector for the muons******
	1248	// Smears the muon pT and changes the 4-momentum accordingly
	1249	void RESOLution::SmearMu(TLorentzVector &muon) {
	1250	// the 'muon' variable will be changed by the function
	1251	float pt = muon.Pt(); // before smearing
[61]	1252	float ptS=pt;
	1253
[94]	1254	if(fabs(muon.Eta()) < CEN_max_mu )
[61]	1255	{
[526]	1256	ptS = grandom->Gaus(pt, MU_SmearPt*pt ); // after smearing // \sigma/E = C + N/E + S/\sqrt{E}
[61]	1257	}
	1258	muon.SetPtEtaPhiE(ptS, muon.Eta(), muon.Phi(), ptS*cosh(muon.Eta()));
[2]	1259
	1260	if(muon.E() < 0)muon.SetPxPyPzE(0,0,0,0); // no negative values after smearing !
	1261	}
	1262
	1263
	1264	// ********Provides the smeared TLorentzVector for the hadrons******
	1265	// Smears the hadron 4-momentum
	1266	void RESOLution::SmearHadron(TLorentzVector &hadron, const float frac)
	1267	// the 'hadron' variable will be changed by the function
	1268	// the 'frac' variable describes the long-living particles. Should be 0.7 for K0S and Lambda, 1. otherwise
	1269	{
	1270	float energy = hadron.E(); // before smearing
	1271	float energyS = 0.0; // after smearing // \sigma/E = C + N/E + S/\sqrt{E}
	1272	float energy_ecal = (1.0 - frac)*energy; // electromagnetic calorimeter
	1273	float energy_hcal = frac*energy; // hadronic calorimeter
	1274	// frac takes into account the decay of long-living particles, that decay in the calorimeters
	1275	// some of the particles decay mostly in the ecal, some mostly in the hcal
	1276
[31]	1277	float energyS1,energyS2;
[94]	1278	if(fabs(hadron.Eta()) < CEN_max_calo_cen) {
[526]	1279	energyS1 = grandom->Gaus(energy_hcal, sqrt(
[494]	1280	pow(HAD_Ncen,2) +
	1281	pow(HAD_Ccen*energy_hcal,2) +
	1282	pow(HAD_Scen*sqrt(energy_hcal),2) )) ;
[9]	1283
[526]	1284	energyS2 = grandom->Gaus(energy_ecal, sqrt(
[32]	1285	pow(ELG_Ncen,2) +
	1286	pow(ELG_Ccen*energy_ecal,2) +
	1287	pow(ELG_Scen*sqrt(energy_ecal),2) ) );
[9]	1288
[10]	1289	energyS = ((energyS1>0)?energyS1:0) + ((energyS2>0)?energyS2:0);
[55]	1290	}
[494]	1291	else if(fabs(hadron.Eta()) >= CEN_max_calo_cen && fabs(hadron.Eta()) < CEN_max_calo_ec) {
[526]	1292	energyS1 = grandom->Gaus(energy_hcal, sqrt(
[494]	1293	pow(HAD_Nec,2) +
	1294	pow(HAD_Cec*energy_hcal,2) +
	1295	pow(HAD_Sec*sqrt(energy_hcal),2) )) ;
	1296
[526]	1297	energyS2 = grandom->Gaus(energy_ecal, sqrt(
[494]	1298	pow(ELG_Nec,2) +
	1299	pow(ELG_Cec*energy_ecal,2) +
	1300	pow(ELG_Sec*sqrt(energy_ecal),2) ) );
	1301
	1302	energyS = ((energyS1>0)?energyS1:0) + ((energyS2>0)?energyS2:0);
	1303	}
	1304	else if(fabs(hadron.Eta()) >= CEN_max_calo_ec && fabs(hadron.Eta()) < CEN_max_calo_fwd){
[526]	1305	energyS = grandom->Gaus(energy, sqrt(
[494]	1306	pow(HAD_Nfwd,2) +
	1307	pow(HAD_Cfwd*energy,2) +
	1308	pow(HAD_Sfwd*sqrt(energy),2) ));
[55]	1309	}
	1310
[10]	1311
	1312
[2]	1313	hadron.SetPtEtaPhiE(energyS/cosh(hadron.Eta()),hadron.Eta(), hadron.Phi(), energyS);
	1314
	1315	if(hadron.E() < 0)hadron.SetPxPyPzE(0,0,0,0);
	1316	}
	1317
[74]	1318	//******************************************************************************************
	1319
[264]	1320	//void RESOLution::SortedVector(vector<ParticleUtil> &vect)
	1321	void RESOLution::SortedVector(vector<D_Particle> &vect)
[74]	1322	{
	1323	int i,j = 0;
	1324	TLorentzVector tmp;
	1325	bool en_desordre = true;
	1326	int entries=vect.size();
	1327	for(i = 0 ; (i < entries) && en_desordre; i++)
	1328	{
	1329	en_desordre = false;
	1330	for(j = 1 ; j < entries - i ; j++)
	1331	{
	1332	if ( vect[j].Pt() > vect[j-1].Pt() )
	1333	{
[264]	1334	//ParticleUtil tmp = vect[j-1];
	1335	D_Particle tmp = vect[j-1];
[74]	1336	vect[j-1] = vect[j];
	1337	vect[j] = tmp;
	1338	en_desordre = true;
	1339	}
	1340	}
	1341	}
	1342	}
	1343
[2]	1344	// ********Provides the energy in the cone of radius TAU_CONE_ENERGY for the tau identification******
	1345	// to be taken into account, a calo tower should
	1346	// 1) have a transverse energy \f$ E_T = \sqrt{E_X^2 + E_Y^2} \f$ above a given threshold
	1347	// 2) be inside a cone with a radius R and the axis defined by (eta,phi)
	1348	double RESOLution::EnergySmallCone(const vector<PhysicsTower> &towers, const float eta, const float phi) {
	1349	double Energie=0;
	1350	for(unsigned int i=0; i < towers.size(); i++) {
[94]	1351	if(towers[i].fourVector.pt() < JET_seed) continue;
	1352	if((DeltaR(phi,eta,towers[i].fourVector.phi(),towers[i].fourVector.eta()) < TAU_energy_scone)) {
[2]	1353	Energie += towers[i].fourVector.E;
	1354	}
	1355	}
	1356	return Energie;
	1357	}
	1358
	1359
	1360	// ********Provides the number of tracks in the cone of radius TAU_CONE_TRACKS for the tau identification******
	1361	// to be taken into account, a track should
	1362	// 1) avec a transverse momentum \$f p_T \$ above a given threshold
	1363	// 2) be inside a cone with a radius R and the axis defined by (eta,phi)
	1364	// IMPORTANT REMARK !!!!!
[287]	1365	// NEW : "charge" will contain the sum of all charged tracks in the cone TAU_track_scone
	1366	unsigned int RESOLution::NumTracks(float& charge, const vector<TRootTracks> &tracks, const float pt_track, const float eta, const float phi) {
	1367	unsigned int numbtrack=0; // number of track in the tau-jet cone, which is smaller than R;
	1368	charge=0;
[2]	1369	for(unsigned int i=0; i < tracks.size(); i++) {
[287]	1370	if(tracks[i].PT < pt_track ) continue;
[319]	1371	//float dr = DeltaR(phi,eta,tracks[i].PhiOuter,tracks[i].EtaOuter);
[287]	1372	float dr = DeltaR(phi,eta,tracks[i].Phi,tracks[i].Eta);
	1373	if (dr > TAU_track_scone) continue;
	1374	numbtrack++;
	1375	charge += tracks[i].Charge; // total charge in the cone for Tau-jet
[2]	1376	}
[287]	1377	return numbtrack;
[2]	1378	}
	1379
	1380	//* Returns the PID of the particle with the highest energy, in a cone with a radius CONERADIUS and an axis (eta,phi) *******
	1381	//used by Btaggedjet
	1382	///// Attention : bug removed => CONERADIUS/2 -> CONERADIUS !!
[350]	1383	int RESOLution::Bjets(const TSimpleArray<TRootC::GenParticle> &subarray, const float& eta, const float& phi) {
[671]	1384	int flav_max=0;
[2]	1385	int Ppid=0;
	1386	if(subarray.GetEntries()>0) {
	1387	for(int i=0; i < subarray.GetEntries();i++) { // should have pt>PT_JETMIN and a small cone radius (r<CONE_JET)
[671]	1388	float genDeltaR = DeltaR(subarray[i]->Phi,subarray[i]->Eta,phi,eta);
	1389	if(genDeltaR < JET_coneradius && (abs( subarray[i]->PID) < 6 && abs(subarray[i]->PID) > flav_max)) {
	1390	flav_max=abs(subarray[i]->PID);
[2]	1391	Ppid=abs(subarray[i]->PID);
	1392	}
	1393	}
	1394	}
[671]	1395	return Ppid;
[2]	1396	}
	1397
[558]	1398	//***************** Simulates the b-tagging efficiency for real bjet, or the misendentification for other jets**************
[526]	1399
[350]	1400	bool RESOLution::Btaggedjet(const TLorentzVector &JET, const TSimpleArray<TRootC::GenParticle> &subarray) {
[526]	1401
[558]	1402	const float ptmax = 9E6; // GeV/c
	1403	bool tag = false;
	1404	string efficiency_formula="";
	1405
	1406	// selects the correct formula
[671]	1407	int jet_id = Bjets(subarray,JET.Eta(),JET.Phi()); // gets the flavor of the jet
[558]	1408	switch (jet_id) {
[671]	1409	case 5: efficiency_formula = BTAG_func_b; break;
	1410	case 4: efficiency_formula = BTAG_func_c; break;
[558]	1411	default: efficiency_formula = BTAG_func_l; break;
[526]	1412	}
	1413
[558]	1414	TF2 efficiency("efficiency",efficiency_formula.c_str(),0,ptmax,-CEN_max_tracker,CEN_max_tracker); // efficiency function wrt Pt
	1415	float x = grandom->Uniform(); // randomisation
	1416	tag = (x < efficiency.Eval(JET.Pt(),JET.Eta())) ? true : false;
[526]	1417
[558]	1418	return tag;
[2]	1419	}
	1420
[526]	1421
[31]	1422	//*********************Isolation criteria*********************
	1423	//****************************************************************
[557]	1424	bool RESOLution::Isolation(const D_Particle& part, const vector<TRootTracks> &tracks, const float& pt_second_track, float& ptiso )
[31]	1425	{
	1426	bool isolated = false;
[321]	1427	ptiso = 0; // sum of all track pt in isolation cone
	1428	float deltar=1E99; // Initial value; should be high; no further repercussion
	1429
	1430	// loop on all tracks, with p_t above threshold, close enough from the charged lepton
	1431	for(unsigned int i=0; i < tracks.size(); i++) {
	1432	if(tracks[i].PT < pt_second_track) continue; // ptcut on tracks
	1433	float genDeltaR = DeltaR(part.Phi(),part.Eta(),tracks[i].Phi,tracks[i].Eta);
[31]	1434	if(
[321]	1435	(genDeltaR==0) // rejets the track of the particle itself
[31]	1436	) continue ;
[665]	1437	if (genDeltaR < deltar) deltar=genDeltaR; // finds the closest track
[321]	1438
	1439	// as long as (genDeltaR==0) is put above, the particle itself is not taken into account
[557]	1440	if( genDeltaR < ISOL_trk_Cone && tracks[i].PT>ISOL_trk_PT) ptiso += tracks[i].PT; // dR cut on tracks
[31]	1441	}
[557]	1442	if(deltar > ISOL_trk_Cone) isolated = true;
[31]	1443	return isolated;
	1444	}
	1445
[556]	1446	float RESOLution::IsolationSumEt(const float iPhi, const float iEta, D_CaloTowerList & towers)
	1447	{
	1448	float etiso = 0; // sum of all track pt in isolation cone
	1449	for (unsigned int i=0; i< towers.size(); i++) {
[558]	1450	if( (DeltaR(iPhi,iEta,towers[i].getPhi(),towers[i].getEta()) < ISOL_calo_Cone) && (towers[i].getET() > ISOL_calo_ET) ) {
[556]	1451	etiso += towers[i].getET(); // dR cut on tracks
	1452	}
	1453	}
	1454	return etiso;
	1455	}
	1456
	1457
	1458
[558]	1459	// ******* Calorimetric isolation -- for LHCO only
[550]	1460	float RESOLution::CaloIsolation(const D_Particle& part, const D_CaloTowerList & towers, const float iPhi, const float iEta, const int iTower) {
	1461	// iPhi and iEta are the coordinates of the center of the tower crossed by the particle
	1462	// iPhi is in [-pi ; pi] and iEta is in [-etamax ; etamax]
	1463	// iTower is the index of the tower, in [0, n_tower]. iTower points only towers in positive range
	1464
[392]	1465	float et_sum=0;
	1466
[557]	1467	unsigned int N = ISOL_calo_Grid;
[550]	1468	if(N%2 ==0 \|\| N==0) { cout << "Warning: ISOL_Calo_Grid must be odd. ISOL_Calo_Grid = 3 will be used\n"; N=3;}
	1469	int index= iTower; // index of the central tower of the grid in TOWER_eta_edges[.];
	1470	// !! TOWER_eta_edges is only with eta>0
	1471
[551]	1472	//cout << "iEta= " << iEta << "\tiPhi= " << iPhi << endl;
[392]	1473	if(index != iUNDEFINED) {
[550]	1474	int range = (int) (N-1)/2; // the (int) conversion is needed, as N is "unsigned int"
	1475	for(int i= -range; i<= range; i++) { // shift of index in eta: e.g.: i = -1, 0, 1 if N=3;
	1476	// goal : to identify the center of each cell in the NxN grid.
	1477	// the eta/phi coordinates of each center will be in (eta_ith_tower,phi_ith_tower)
[551]	1478	unsigned int i_eta = (index+i>=0) ? index + i : (int) -1*(index + i +1); // i_eta is the index in [0; index_eta_max]
	1479	if(i_eta>=TOWER_number) continue; // avoid going too far: TOWER_dphi[TOWER_number], TOWER_eta_edges[TOWER_number+1]
	1480	//cout << "i_eta" << i_eta << "\t TOWER_number = " << TOWER_number << endl;
[392]	1481
[550]	1482	float dphi = TOWER_dphi[i_eta]*pi/180.; // in rad // size in phi of the cells for this eta
	1483	float eta_ith_tower = (TOWER_eta_edges[i_eta]+TOWER_eta_edges[i_eta+1])/2.0;
	1484	if(iEta<0) eta_ith_tower *= -1; // needed if the central tower is in negative region
[551]	1485	if(index+i<0) eta_ith_tower *= -1; // needed if the "eta=0"-axis is crossed by the grid
[550]	1486	//cout << "pour eta_ith_tower=" << eta_ith_tower << ", on va dphi = " << dphi << endl;
	1487	// !!! at this point, eta_ith_tower is the value in eta of the center of the current calo cell
	1488
	1489	for(int j= -range; j<= range; j++) { // shift of index in phi. e.g.: j = -1, 0, 1 if N=3;
	1490	float phi_ith_tower = iPhi + j*dphi; // in radian
	1491	//cout << "eta_ith_tower= " << eta_ith_tower << "\tphi_ith_tower= " << phi_ith_tower << "\tj= " << j << "\tdphi=" << dphi << endl;
	1492	if(phi_ith_tower > pi) phi_ith_tower -= 2.*pi;
	1493	else if (phi_ith_tower < -pi) phi_ith_tower += 2.*pi;
	1494
	1495	D_CaloTower calMuon(towers.getElement(eta_ith_tower,phi_ith_tower));
[551]	1496	//cout << "eta_ith_tower= " << eta_ith_tower << "\tphi_ith_tower= " << phi_ith_tower <<"\t"
	1497	// << "calMuon.getEta= " << calMuon.getEta() << "\tcalMuon.getPhi()= " << calMuon.getPhi() <<"\t";
[550]	1498
[558]	1499	//if(calMuon.getEta() != UNDEFINED && calMuon.getET() > ISOL_calo_ET) {
	1500	if(calMuon.getEta() != UNDEFINED && calMuon.getET() > 0.0) {
[551]	1501	et_sum += calMuon.getET();
	1502	//cout << calMuon.getET() << " GeV";
[550]	1503	}
[551]	1504	//else cout << " not active";
	1505	} // j-loop (phi indices)
	1506	} // i-loop (eta indices)
[550]	1507	} // undefined index
	1508	else {
	1509	if (CEN_max_mu < CEN_max_calo_fwd)
[392]	1510	cout << " ERROR in RESOLution::CaloIsolation: 'muon'-tower not found! " << endl;
[551]	1511	} // should never happen ! this would be a bug
	1512	return (part.Pt()==0)? UNDEFINED: et_sum/part.Pt();
[321]	1513	}
[31]	1514
[321]	1515
[71]	1516	//******** returns a segmented value for eta and phi, for calo towers ***
[550]	1517	int RESOLution::BinEtaPhi(const float phi, const float eta, float& iPhi, float& iEta){
[423]	1518	iEta = UNDEFINED;
	1519	int index= iUNDEFINED;
	1520	for (unsigned int i=1; i< TOWER_number+1; i++) {
	1521	if(fabs(eta)>TOWER_eta_edges[i-1] && fabs(eta)<=TOWER_eta_edges[i]) {
	1522	iEta = (eta>0) ? ((TOWER_eta_edges[i-1]+TOWER_eta_edges[i])/2.0) : -((TOWER_eta_edges[i-1]+TOWER_eta_edges[i])/2.0);
	1523	index = i-1;
	1524	break;
	1525	}
	1526	}
	1527	iPhi = UNDEFINED;
[550]	1528	if(index==iUNDEFINED) return index;
	1529
[423]	1530	float dphi = TOWER_dphi[index]*pi/180.;
	1531	for (unsigned int i=1; i < 360/TOWER_dphi[index]; i++ ) {
	1532	float low = -pi+(i-1)*dphi;
	1533	float high= low+dphi;
	1534	if(phi > low && phi <= high ){
	1535	iPhi = (low+high)/2.0;
	1536	break;
	1537	}
	1538	}
	1539	if (phi > pi-dphi) iPhi = pi-dphi;
[550]	1540	return index;
[423]	1541	}
[71]	1542
[2]	1543	//************************** Returns the delta Phi **************************
	1544	float DeltaPhi(const float phi1, const float phi2) {
[244]	1545	float deltaphi=phi1-phi2; // in here, -pi < phi < pi
	1546	if(fabs(deltaphi) > pi) {
	1547	deltaphi=2.*pi -fabs(deltaphi);// put deltaphi between 0 and pi
[219]	1548	}
[2]	1549	else deltaphi=fabs(deltaphi);
	1550
	1551	return deltaphi;
	1552	}
	1553
	1554	//************************** Returns the delta R**************************
	1555	float DeltaR(const float phi1, const float eta1, const float phi2, const float eta2) {
	1556	return sqrt(pow(DeltaPhi(phi1,phi2),2) + pow(eta1-eta2,2));
	1557	}
	1558
	1559	int sign(const int myint) {
	1560	if (myint >0) return 1;
	1561	else if (myint <0) return -1;
	1562	else return 0;
	1563	}
	1564
	1565	int sign(const float myfloat) {
	1566	if (myfloat >0) return 1;
	1567	else if (myfloat <0) return -1;
	1568	else return 0;
	1569	}
	1570
[270]	1571	int ChargeVal(const int pid)
[55]	1572	{
[380]	1573	cout << "ChargeVal :: deprecated function, do not use it anymore" << endl;
[55]	1574	int charge;
	1575	if(
	1576	(pid == pGAMMA) \|\|
	1577	(pid == pPI0) \|\|
	1578	(pid == pK0L) \|\|
	1579	(pid == pN) \|\|
	1580	(pid == pSIGMA0) \|\|
	1581	(pid == pDELTA0) \|\|
	1582	(pid == pK0S) // not charged particles : invisible by tracker
	1583	)
	1584	charge = 0;
[376]	1585	else charge = sign(pid);
[55]	1586	return charge;
	1587
[2]	1588	}
[380]	1589
	1590	//------------------------------------------------------------------------------
	1591	void RESOLution::ReadParticleDataGroupTable() {
	1592
	1593	string temp_string;
	1594	istringstream curstring;
	1595
	1596	ifstream fichier_a_lire(PdgTableFilename.c_str());
	1597	if(!fichier_a_lire.good()) {
	1598	cout <<"** ERROR: PDG Table ("<< PdgTableFilename
	1599	<< ") not found! exit. **" << endl;
	1600	exit(1);
	1601	return;
	1602	}
	1603	// first three lines of the file are useless
	1604	getline(fichier_a_lire,temp_string);
	1605	getline(fichier_a_lire,temp_string);
	1606	getline(fichier_a_lire,temp_string);
	1607
	1608
	1609	while (getline(fichier_a_lire,temp_string)) {
	1610	curstring.clear(); // needed when using several times istringstream::str(string)
	1611	curstring.str(temp_string);
[469]	1612	long int ID; std::string name; int charge; float mass; float width; float lifetime;
[380]	1613	// ID name chg mass total width lifetime
	1614	// 1 d -1 0.33000 0.00000 0.00000E+00
[404]	1615	// in the table, the charge is in units of e+/3
	1616	// the total width is in GeV
	1617	// the lifetime is ctau in mm
[380]	1618	curstring >> ID >> name >> charge >> mass >> width >> lifetime;
[404]	1619	PdgParticle particle(ID,name,mass,charge/3.,width,lifetime/1000.);
[380]	1620	PDGtable.insert(ID,particle);
	1621	//PdgTable.insert(pair<int,PdgParticle>(ID,particle));
	1622	//cout << PDGtable[ID].name() << "\t" << PDGtable[ID].mass() << "\t" << PDGtable[ID].charge() << endl;
	1623	}
	1624
	1625	} // ReadParticleDataGroupTable
[399]	1626
	1627
	1628	// to be improved in order to avoid code repetition
	1629	// sorry, no time to do it right now (XR, 19/05/2009)
	1630	void time_report(const TStopwatch& global,const TStopwatch& loop,const TStopwatch& trigger,const TStopwatch& frog,const TStopwatch& lhco, const int flag_frog, const int flag_trigger, const int flag_lhco, const string& LogName, const Long64_t allEntries) {
	1631
	1632	TStopwatch globalwatch(global), loopwatch(loop), triggerwatch(trigger), frogwatch(frog), lhcowatch(lhco);
	1633
	1634	cout <<" "<< endl;
	1635	cout <<" ################## Time report ################# "<< endl;
	1636	cout << left << setw(32) <<"** Time report for "<<""
	1637	<< left << setw(15) << allEntries <<""
	1638	<< right << setw(22) <<"events **"<<endl;
	1639	cout <<" "<< endl;
	1640	cout << left << setw(10) <<"**"<<""
	1641	<< left << setw(15) <<"Time (s):"<<""
	1642	<< right << setw(15) <<"CPU"<<""
	1643	<< right << setw(15) <<"Real"<<""
	1644	<< right << setw(14) <<"**"<<endl;
	1645	cout << left << setw(10) <<"**"<<""
	1646	<< left << setw(15) <<" + Global:"<<""
	1647	<< right << setw(15) <<globalwatch.CpuTime()<<""
	1648	<< right << setw(15) <<globalwatch.RealTime()<<""
	1649	<< right << setw(14) <<"**"<<endl;
	1650	cout << left << setw(10) <<"**"<<""
	1651	<< left << setw(15) <<" + Events:"<<""
	1652	<< right << setw(15) <<loopwatch.CpuTime()<<""
	1653	<< right << setw(15) <<loopwatch.RealTime()<<""
	1654	<< right << setw(14) <<"**"<<endl;
	1655	if(flag_trigger == 1)
	1656	{
	1657	cout << left << setw(10) <<"**"<<""
	1658	<< left << setw(15) <<" + Trigger:"<<""
	1659	<< right << setw(15) <<triggerwatch.CpuTime()<<""
	1660	<< right << setw(15) <<triggerwatch.RealTime()<<""
	1661	<< right << setw(14) <<"**"<<endl;
	1662	}
	1663	if(flag_frog == 1)
	1664	{
	1665	cout << left << setw(10) <<"**"<<""
	1666	<< left << setw(15) <<" + Frog:"<<""
	1667	<< right << setw(15) <<frogwatch.CpuTime()<<""
	1668	<< right << setw(15) <<frogwatch.RealTime()<<""
	1669	<< right << setw(14) <<"**"<<endl;
	1670	}
	1671	if(flag_lhco == 1)
	1672	{
	1673	cout << left << setw(10) <<"**"<<""
	1674	<< left << setw(15) <<" + LHCO:"<<""
	1675	<< right << setw(15) <<lhcowatch.CpuTime()<<""
	1676	<< right << setw(15) <<lhcowatch.RealTime()<<""
	1677	<< right << setw(14) <<"**"<<endl;
	1678	}
	1679
	1680
	1681	ofstream f_out(LogName.c_str(),ios_base::app);
	1682
	1683	f_out <<"** *"<< endl;
	1684	f_out <<"** ################## Time report ################# *"<< endl;
	1685	f_out << left << setw(32) <<"** Time report for "<<""
	1686	<< left << setw(15) << allEntries <<""
[472]	1687	<< right << setw(23) <<"events *"<<endl;
[399]	1688	f_out <<"** *"<< endl;
	1689	f_out << left << setw(10) <<"**"<<""
	1690	<< left << setw(15) <<"Time (s):"<<""
	1691	<< right << setw(15) <<"CPU"<<""
	1692	<< right << setw(15) <<"Real"<<""
	1693	<< right << setw(15) <<" *"<<endl;
	1694	f_out << left << setw(10) <<"**"<<""
	1695	<< left << setw(15) <<" + Global:"<<""
	1696	<< right << setw(15) <<globalwatch.CpuTime()<<""
	1697	<< right << setw(15) <<globalwatch.RealTime()<<""
	1698	<< right << setw(15) <<" *"<<endl;
	1699	f_out << left << setw(10) <<"**"<<""
	1700	<< left << setw(15) <<" + Events:"<<""
	1701	<< right << setw(15) <<loopwatch.CpuTime()<<""
	1702	<< right << setw(15) <<loopwatch.RealTime()<<""
	1703	<< right << setw(15) <<" *"<<endl;
	1704	if(flag_trigger == 1)
	1705	{
	1706	f_out << left << setw(10) <<"**"<<""
	1707	<< left << setw(15) <<" + Trigger:"<<""
	1708	<< right << setw(15) <<triggerwatch.CpuTime()<<""
	1709	<< right << setw(15) <<triggerwatch.RealTime()<<""
	1710	<< right << setw(15) <<" *"<<endl;
	1711	}
	1712	if(flag_frog == 1)
	1713	{
	1714	f_out << left << setw(10) <<"**"<<""
	1715	<< left << setw(15) <<" + Frog:"<<""
	1716	<< right << setw(15) <<frogwatch.CpuTime()<<""
	1717	<< right << setw(15) <<frogwatch.RealTime()<<""
	1718	<< right << setw(15) <<" *"<<endl;
	1719	}
	1720	if(flag_lhco == 1)
	1721	{
	1722	f_out << left << setw(10) <<"**"<<""
	1723	<< left << setw(15) <<" + LHCO:"<<""
	1724	<< right << setw(15) <<lhcowatch.CpuTime()<<""
	1725	<< right << setw(15) <<lhcowatch.RealTime()<<""
	1726	<< right << setw(15) <<" *"<<endl;
	1727	}
[472]	1728	f_out << left << setw(16) << "** " << ""
	1729	<< left << setw(52) << get_time_date() << "**" << endl;
	1730
	1731
[399]	1732	f_out<<"* *"<<"\n";
	1733	f_out<<"* *"<<"\n";
	1734	f_out<<"#....................................................................*"<<"\n";
	1735	f_out<<"#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"<<"\n";
	1736
	1737	f_out.close();
	1738
	1739	}
[404]	1740
	1741	void print_header() {
	1742	cout << endl << endl;
	1743
	1744	cout <<"*********************************************************************"<< endl;
	1745	cout <<"*********************************************************************"<< endl;
	1746	cout <<" "<< endl;
	1747	cout <<" Welcome to "<< endl;
	1748	cout <<" "<< endl;
	1749	cout <<" "<< endl;
	1750	cout <<" .ddddddd- lL hH "<< endl;
	1751	cout <<" -Dd` `dD: Ll hH` "<< endl;
	1752	cout <<" dDd dDd eeee. lL .pp+pp Hh+hhh` -eeee- `sssss "<< endl;
	1753	cout <<" -Dd `DD ee. ee Ll .Pp. PP Hh. HH. ee. ee sSs "<< endl;
	1754	cout <<" dD` dDd eEeee: lL. pP. pP hH hH` eEeee:` -sSSSs. "<< endl;
	1755	cout <<" .Dd :dd eE. LlL PpppPP Hh Hh eE sSS "<< endl;
	1756	cout <<" dddddd:. eee+: lL. pp. hh. hh eee+ sssssS "<< endl;
	1757	cout <<" Pp "<< endl;
	1758	cout <<" "<< endl;
	1759	cout <<" Delphes, a framework for the fast simulation "<< endl;
	1760	cout <<" of a generic collider experiment "<< endl;
	1761	cout <<" arXiv:0903.2225v1 [hep-ph] "<< endl;
	1762	cout <<" "<< endl;
[567]	1763	cout <<" --- Version 1.9 of Delphes --- "<< endl;
	1764	cout <<" Last date of change: 7 May 2010 "<< endl;
[404]	1765	cout <<" "<< endl;
	1766	cout <<" "<< endl;
	1767	cout <<" This package uses: "<< endl;
	1768	cout <<" ------------------ "<< endl;
	1769	cout <<" FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210] "<< endl;
	1770	cout <<" Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2] "<< endl;
	1771	cout <<" FROG: [hep-ex/0901.2718v1] "<< endl;
	1772	cout <<" "<< endl;
	1773	cout <<"-----------------------------------------------------------------"<< endl;
	1774	cout <<" "<< endl;
	1775	cout <<" Main authors: "<< endl;
	1776	cout <<" ------------- "<< endl;
	1777	cout <<" "<< endl;
	1778	cout <<" SÃ©verine Ovyn Xavier Rouby "<< endl;
	1779	cout <<" severine.ovyn@uclouvain.be xavier.rouby@cern "<< endl;
	1780	cout <<" Center for Particle Physics and Phenomenology (CP3) "<< endl;
	1781	cout <<" Universite Catholique de Louvain (UCL) "<< endl;
	1782	cout <<" Louvain-la-Neuve, Belgium "<< endl;
	1783	cout <<" "<< endl;
	1784	cout <<"-----------------------------------------------------------------"<< endl;
	1785	cout <<" "<< endl;
	1786	cout <<" Former Delphes versions and documentation can be found on : "<< endl;
	1787	cout <<" http://www.fynu.ucl.ac.be/delphes.html "<< endl;
	1788	cout <<" "<< endl;
	1789	cout <<" "<< endl;
[528]	1790	cout <<" Disclaimer: this program comes without guarantees. "<< endl;
	1791	cout <<" Beware of errors and please give us "<< endl;
	1792	cout <<" your feedbacks about potential bugs. "<< endl;
[404]	1793	cout <<" "<< endl;
	1794	cout <<"*********************************************************************"<< endl;
	1795	cout <<"*********************************************************************"<< endl;
	1796
	1797	}
[465]	1798
	1799	string get_time_date() {
	1800	time_t rawtime;
	1801	struct tm * timeinfo;
	1802
	1803	time ( &rawtime );
	1804	timeinfo = localtime ( &rawtime );
	1805
	1806	char temp[100];
[553]	1807	//sprintf(temp,"%i/%i/%i %i:%i:%i",timeinfo->tm_mday,timeinfo->tm_mon+1,timeinfo->tm_year+1900,timeinfo->tm_hour,timeinfo->tm_min,timeinfo->tm_sec);
	1808	sprintf(temp,"%i",timeinfo->tm_mday);
	1809	if(timeinfo->tm_mon+1<10) sprintf(temp,"%s/0%d",temp,timeinfo->tm_mon+1);else sprintf(temp,"%s/%d",temp,timeinfo->tm_mon+1);
	1810	sprintf(temp,"%s/%d %d",temp,timeinfo->tm_year+1900,timeinfo->tm_hour);
	1811	if(timeinfo->tm_min<10) sprintf(temp,"%s:0%d",temp,timeinfo->tm_min); else sprintf(temp,"%s:%d",temp,timeinfo->tm_min);
	1812	if(timeinfo->tm_sec<10) sprintf(temp,"%s:0%d",temp,timeinfo->tm_sec); else sprintf(temp,"%s:%d",temp,timeinfo->tm_sec);
[465]	1813	string tempstring(temp);
	1814	return tempstring;
	1815	}
	1816
[494]	1817	void warning(const string oldname, const string newname) {
	1818	string text = "** Warning in datacard: " + oldname + " deprecated. Use " + newname + " instead.";
	1819	cout << left << setw(67) << text
	1820	<< right << setw(2) <<"**" << endl;
	1821	}

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

Context Navigation

source: svn/trunk/src/SmearUtil.cc@ 752

Download in other formats: