[2] | 1 | /*
|
---|
| 2 | ---- Delphes ----
|
---|
| 3 | A Fast Simulator for general purpose LHC detector
|
---|
| 4 | S. Ovyn ~~~~ severine.ovyn@uclouvain.be
|
---|
| 5 |
|
---|
| 6 | Center for Particle Physics and Phenomenology (CP3)
|
---|
| 7 | Universite Catholique de Louvain (UCL)
|
---|
| 8 | Louvain-la-Neuve, Belgium
|
---|
| 9 | */
|
---|
| 10 |
|
---|
| 11 | /// \file Delphes.cpp
|
---|
| 12 | /// \brief executable for the Delphes
|
---|
| 13 |
|
---|
| 14 | #include "TChain.h"
|
---|
| 15 | #include "TApplication.h"
|
---|
| 16 |
|
---|
| 17 | #include "Utilities/ExRootAnalysis/interface/ExRootTreeReader.h"
|
---|
| 18 | #include "Utilities/ExRootAnalysis/interface/ExRootTreeWriter.h"
|
---|
| 19 | #include "Utilities/ExRootAnalysis/interface/ExRootTreeBranch.h"
|
---|
| 20 |
|
---|
| 21 | #include "H_BeamParticle.h"
|
---|
| 22 | #include "H_BeamLine.h"
|
---|
| 23 | #include "H_RomanPot.h"
|
---|
| 24 |
|
---|
| 25 | #include "interface/DataConverter.h"
|
---|
| 26 | #include "interface/HEPEVTConverter.h"
|
---|
| 27 | #include "interface/LHEFConverter.h"
|
---|
| 28 | #include "interface/STDHEPConverter.h"
|
---|
| 29 |
|
---|
| 30 | #include "interface/SmearUtil.h"
|
---|
[11] | 31 | #include "Utilities/Fastjet/include/fastjet/PseudoJet.hh"
|
---|
| 32 | #include "Utilities/Fastjet/include/fastjet/ClusterSequence.hh"
|
---|
[2] | 33 |
|
---|
[11] | 34 | // get info on how fastjet was configured
|
---|
| 35 | #include "Utilities/Fastjet/include/fastjet/config.h"
|
---|
| 36 |
|
---|
| 37 | // make sure we have what is needed
|
---|
| 38 | #ifdef ENABLE_PLUGIN_SISCONE
|
---|
| 39 | # include "Utilities/Fastjet/plugins/SISCone/SISConePlugin.hh"
|
---|
| 40 | #endif
|
---|
| 41 | #ifdef ENABLE_PLUGIN_CDFCONES
|
---|
| 42 | # include "Utilities/Fastjet/plugins/CDFCones/CDFMidPointPlugin.hh"
|
---|
| 43 | # include "Utilities/Fastjet/plugins/CDFCones/CDFJetCluPlugin.hh"
|
---|
| 44 | #endif
|
---|
| 45 |
|
---|
| 46 | #include<vector>
|
---|
| 47 | #include<iostream>
|
---|
| 48 |
|
---|
| 49 |
|
---|
| 50 |
|
---|
[2] | 51 | using namespace std;
|
---|
| 52 |
|
---|
| 53 | //------------------------------------------------------------------------------
|
---|
| 54 | void todo(string filename) {
|
---|
| 55 | ifstream infile(filename.c_str());
|
---|
| 56 | cout << "** TODO list ..." << endl;
|
---|
| 57 | while(infile.good()) {
|
---|
| 58 | string temp;
|
---|
| 59 | getline(infile,temp);
|
---|
| 60 | cout << "*" << temp << endl;
|
---|
| 61 | }
|
---|
| 62 | cout << "** done...\n";
|
---|
| 63 | }
|
---|
| 64 |
|
---|
| 65 | //------------------------------------------------------------------------------
|
---|
| 66 |
|
---|
| 67 | int main(int argc, char *argv[])
|
---|
| 68 | {
|
---|
| 69 | int appargc = 2;
|
---|
[29] | 70 | char *appName = "JetsSim";
|
---|
[2] | 71 | char *appargv[] = {appName, "-b"};
|
---|
| 72 | TApplication app(appName, &appargc, appargv);
|
---|
| 73 |
|
---|
| 74 | if(argc != 4 && argc != 3) {
|
---|
| 75 | cout << " Usage: " << argv[0] << " input_file" << " output_file" << " data_card " << endl;
|
---|
| 76 | cout << " input_list - list of files in Ntpl, StdHep of LHEF format," << endl;
|
---|
| 77 | cout << " output_file - output file." << endl;
|
---|
| 78 | cout << " data_card - Datacard containing resolution variables for the detector simulation (optional) "<<endl;
|
---|
| 79 | exit(1);
|
---|
| 80 | }
|
---|
| 81 |
|
---|
| 82 | srand (time (NULL)); /* Initialisation du générateur */
|
---|
| 83 |
|
---|
| 84 | //read the input TROOT file
|
---|
| 85 | string inputFileList(argv[1]), outputfilename(argv[2]);
|
---|
| 86 | if(outputfilename.find(".root") > outputfilename.length() ) {
|
---|
| 87 | cout << "output_file should be a .root file!\n";
|
---|
| 88 | exit(1);
|
---|
| 89 | }
|
---|
| 90 |
|
---|
| 91 | TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE"); // Creates the file, but should be closed just after
|
---|
| 92 | outputFile->Close();
|
---|
| 93 |
|
---|
| 94 | string line;
|
---|
| 95 | ifstream infile(inputFileList.c_str());
|
---|
| 96 | infile >> line; // the first line determines the type of input files
|
---|
| 97 |
|
---|
| 98 | DataConverter *converter=0;
|
---|
| 99 |
|
---|
| 100 | if(strstr(line.c_str(),".hep"))
|
---|
| 101 | {
|
---|
| 102 | cout<<"*************************************************************************"<<endl;
|
---|
| 103 | cout<<"************ StdHEP file format detected **************"<<endl;
|
---|
| 104 | cout<<"************ Starting convertion to TRoot format **************"<<endl;
|
---|
| 105 | cout<<"*************************************************************************"<<endl;
|
---|
| 106 | converter = new STDHEPConverter(inputFileList,outputfilename);//case ntpl file in input list
|
---|
| 107 | }
|
---|
| 108 | else if(strstr(line.c_str(),".lhe"))
|
---|
| 109 | {
|
---|
| 110 | cout<<"*************************************************************************"<<endl;
|
---|
| 111 | cout<<"************ LHEF file format detected **************"<<endl;
|
---|
| 112 | cout<<"************ Starting convertion to TRoot format **************"<<endl;
|
---|
| 113 | cout<<"*************************************************************************"<<endl;
|
---|
| 114 | converter = new LHEFConverter(inputFileList,outputfilename);//case ntpl file in input list
|
---|
| 115 | }
|
---|
| 116 | else if(strstr(line.c_str(),".root"))
|
---|
| 117 | {
|
---|
| 118 | cout<<"*************************************************************************"<<endl;
|
---|
| 119 | cout<<"************ h2root file format detected **************"<<endl;
|
---|
| 120 | cout<<"************ Starting convertion to TRoot format **************"<<endl;
|
---|
| 121 | cout<<"*************************************************************************"<<endl;
|
---|
| 122 | converter = new HEPEVTConverter(inputFileList,outputfilename);//case ntpl file in input list
|
---|
| 123 | }
|
---|
| 124 | else { cout << "*** " << line.c_str() << "\n*** file format not identified\n*** Exiting\n"; return -1;};
|
---|
| 125 |
|
---|
| 126 | TChain chain("GEN");
|
---|
| 127 | chain.Add(outputfilename.c_str());
|
---|
| 128 | ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
|
---|
| 129 | const TClonesArray *branchGen = treeReader->UseBranch("Particle");
|
---|
| 130 | TIter itGen((TCollection*)branchGen);
|
---|
| 131 |
|
---|
| 132 | //write the output root file
|
---|
| 133 | ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputfilename, "Analysis");
|
---|
| 134 | ExRootTreeBranch *branchJet = treeWriter->NewBranch("Jet", TRootJet::Class());
|
---|
| 135 | ExRootTreeBranch *branchTauJet = treeWriter->NewBranch("TauJet", TRootTauJet::Class());
|
---|
| 136 | ExRootTreeBranch *branchElectron = treeWriter->NewBranch("Electron", TRootElectron::Class());
|
---|
| 137 | ExRootTreeBranch *branchMuon = treeWriter->NewBranch("Muon", TRootMuon::Class());
|
---|
| 138 | ExRootTreeBranch *branchPhoton = treeWriter->NewBranch("Photon", TRootPhoton::Class());
|
---|
| 139 | ExRootTreeBranch *branchTracks = treeWriter->NewBranch("Tracks", TRootTracks::Class());
|
---|
| 140 | ExRootTreeBranch *branchETmis = treeWriter->NewBranch("ETmis", TRootETmis::Class());
|
---|
| 141 | ExRootTreeBranch *branchCalo = treeWriter->NewBranch("CaloTower", TRootCalo::Class());
|
---|
| 142 | ExRootTreeBranch *branchZDC = treeWriter->NewBranch("ZDChits", TRootZdcHits::Class());
|
---|
| 143 | ExRootTreeBranch *branchRP220 = treeWriter->NewBranch("RP220hits", TRootRomanPotHits::Class());
|
---|
| 144 | ExRootTreeBranch *branchFP420 = treeWriter->NewBranch("FP420hits", TRootRomanPotHits::Class());
|
---|
| 145 |
|
---|
| 146 |
|
---|
| 147 | TRootGenParticle *particle;
|
---|
| 148 | TRootETmis *elementEtmis;
|
---|
| 149 | TRootElectron *elementElec;
|
---|
| 150 | TRootMuon *elementMu;
|
---|
| 151 | TRootPhoton *elementPhoton;
|
---|
| 152 | TRootJet *elementJet;
|
---|
| 153 | TRootTauJet *elementTauJet;
|
---|
| 154 | TRootTracks *elementTracks;
|
---|
| 155 | TRootCalo *elementCalo;
|
---|
| 156 | TRootZdcHits *elementZdc;
|
---|
| 157 | TRootRomanPotHits *elementRP220, *elementFP420;
|
---|
| 158 |
|
---|
| 159 | //read the datacard input file
|
---|
| 160 | string DetDatacard("");
|
---|
| 161 | if(argc==4) DetDatacard =argv[3];
|
---|
| 162 | RESOLution *DET = new RESOLution();
|
---|
| 163 | DET->ReadDataCard(DetDatacard);
|
---|
| 164 |
|
---|
| 165 |
|
---|
| 166 | TLorentzVector genMomentum(0,0,0,0);
|
---|
| 167 | LorentzVector jetMomentum;
|
---|
| 168 | vector<TLorentzVector> TrackCentral;
|
---|
| 169 | vector<PhysicsTower> towers;
|
---|
| 170 |
|
---|
[11] | 171 | vector<fastjet::PseudoJet> input_particles;//for FastJet algorithm
|
---|
| 172 | vector<fastjet::PseudoJet> inclusive_jets;
|
---|
| 173 | vector<fastjet::PseudoJet> sorted_jets;
|
---|
| 174 |
|
---|
[2] | 175 | //Initialisation of Hector
|
---|
| 176 | extern bool relative_energy;
|
---|
| 177 | relative_energy = true; // should always be true
|
---|
| 178 | extern int kickers_on;
|
---|
| 179 | kickers_on = 1; // should always be 1
|
---|
[11] | 180 |
|
---|
| 181 | // user should provide : (1) optics file for each beamline, and IPname,
|
---|
| 182 | // and offset data (s,x) for optical elements
|
---|
[2] | 183 | H_BeamLine* beamline1 = new H_BeamLine(1,500.);
|
---|
| 184 | beamline1->fill("data/LHCB1IR5_v6.500.tfs",1,"IP5");
|
---|
| 185 | beamline1->offsetElements(120,-0.097);
|
---|
| 186 | H_RomanPot * rp220_1 = new H_RomanPot("rp220_1",220,2000); // RP 220m, 2mm, beam 1
|
---|
| 187 | H_RomanPot * rp420_1 = new H_RomanPot("rp420_1",420,4000); // RP 420m, 4mm, beam 1
|
---|
| 188 | beamline1->add(rp220_1);
|
---|
| 189 | beamline1->add(rp420_1);
|
---|
[11] | 190 |
|
---|
[2] | 191 | H_BeamLine* beamline2 = new H_BeamLine(1,500.);
|
---|
| 192 | beamline2->fill("data/LHCB1IR5_v6.500.tfs",-1,"IP5");
|
---|
| 193 | beamline2->offsetElements(120,+0.097);
|
---|
| 194 | H_RomanPot * rp220_2 = new H_RomanPot("rp220_2",220,2000);// RP 220m, 2mm, beam 2
|
---|
| 195 | H_RomanPot * rp420_2 = new H_RomanPot("rp420_2",420,4000);// RP 420m, 4mm, beam 2
|
---|
| 196 | beamline2->add(rp220_2);
|
---|
| 197 | beamline2->add(rp420_2);
|
---|
| 198 |
|
---|
[19] | 199 | // we will have four jet definitions, and the first three will be
|
---|
| 200 | // plugins
|
---|
| 201 | fastjet::JetDefinition jet_def;
|
---|
| 202 | fastjet::JetDefinition::Plugin * plugins;
|
---|
| 203 |
|
---|
| 204 | switch(DET->JETALGO) {
|
---|
| 205 | default:
|
---|
| 206 | case 1: {
|
---|
| 207 |
|
---|
| 208 | // set up a CDF midpoint jet definition
|
---|
| 209 | #ifdef ENABLE_PLUGIN_CDFCONES
|
---|
| 210 | plugins = new fastjet::CDFJetCluPlugin(DET->C_SEEDTHRESHOLD,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
|
---|
| 211 | jet_def = fastjet::JetDefinition(plugins);
|
---|
| 212 | #else
|
---|
| 213 | plugins = NULL;
|
---|
| 214 | #endif
|
---|
| 215 | }
|
---|
| 216 | break;
|
---|
| 217 |
|
---|
| 218 | case 2: {
|
---|
| 219 |
|
---|
| 220 | // set up a CDF midpoint jet definition
|
---|
| 221 | #ifdef ENABLE_PLUGIN_CDFCONES
|
---|
| 222 | plugins = new fastjet::CDFMidPointPlugin (DET->M_SEEDTHRESHOLD,DET->CONERADIUS,DET->M_CONEAREAFRACTION,DET->M_MAXPAIRSIZE,DET->M_MAXPAIRSIZE,DET->C_OVERLAPTHRESHOLD);
|
---|
| 223 | jet_def = fastjet::JetDefinition(plugins);
|
---|
| 224 | #else
|
---|
| 225 | plugins = NULL;
|
---|
| 226 | #endif
|
---|
| 227 | }
|
---|
| 228 | break;
|
---|
| 229 | case 3: {
|
---|
| 230 | // set up a siscone jet definition
|
---|
| 231 | #ifdef ENABLE_PLUGIN_SISCONE
|
---|
| 232 | int npass = 0; // do infinite number of passes
|
---|
| 233 | double protojet_ptmin = 0.0; // use all protojets
|
---|
| 234 | plugins = new fastjet::SISConePlugin (DET->CONERADIUS,DET->C_OVERLAPTHRESHOLD,npass, protojet_ptmin);
|
---|
| 235 | jet_def = fastjet::JetDefinition(plugins);
|
---|
| 236 | #else
|
---|
| 237 | plugins = NULL;
|
---|
| 238 | #endif
|
---|
| 239 | }
|
---|
| 240 | break;
|
---|
| 241 |
|
---|
| 242 | case 4: {
|
---|
| 243 | jet_def = fastjet::JetDefinition(fastjet::kt_algorithm, DET->CONERADIUS);
|
---|
| 244 | //jet_defs[4] = fastjet::JetDefinition(fastjet::cambridge_algorithm,jet_radius);
|
---|
| 245 | //jet_defs[5] = fastjet::JetDefinition(fastjet::antikt_algorithm,jet_radius);
|
---|
| 246 | }
|
---|
| 247 | break;
|
---|
| 248 | }
|
---|
[2] | 249 |
|
---|
[19] | 250 | // Loop over all events
|
---|
[2] | 251 | Long64_t entry, allEntries = treeReader->GetEntries();
|
---|
| 252 | cout << "** Chain contains " << allEntries << " events" << endl;
|
---|
| 253 | for(entry = 0; entry < allEntries; ++entry)
|
---|
| 254 | {
|
---|
| 255 | TLorentzVector PTmis(0,0,0,0);
|
---|
| 256 | treeReader->ReadEntry(entry);
|
---|
| 257 | treeWriter->Clear();
|
---|
| 258 |
|
---|
| 259 | if((entry % 100) == 0 && entry > 0 ) cout << "** Processing element # " << entry << endl;
|
---|
| 260 |
|
---|
[29] | 261 | TSimpleArray<TRootGenParticle> electron;
|
---|
| 262 | TSimpleArray<TRootGenParticle> muon;
|
---|
[2] | 263 | itGen.Reset();
|
---|
| 264 | TrackCentral.clear();
|
---|
| 265 | towers.clear();
|
---|
[11] | 266 | input_particles.clear();
|
---|
[2] | 267 | TSimpleArray<TRootGenParticle> NFCentralQ;
|
---|
[15] | 268 | inclusive_jets.clear();
|
---|
| 269 | sorted_jets.clear();
|
---|
[2] | 270 |
|
---|
| 271 | // Loop over all particles in event
|
---|
| 272 | while( (particle = (TRootGenParticle*) itGen.Next()) )
|
---|
| 273 | {
|
---|
| 274 | genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
|
---|
| 275 |
|
---|
| 276 | int pid = abs(particle->PID);
|
---|
| 277 | float eta = fabs(particle->Eta);
|
---|
| 278 |
|
---|
| 279 | //subarray of the quarks (i.e. not final particles, i.e status not equal to 1)
|
---|
| 280 | // in the tracker (i.e. for those we know the tracks)
|
---|
| 281 | // with enough p_T
|
---|
| 282 | //// This subarray is needed for the B-jet algorithm
|
---|
| 283 | // optimization for speed : put first PID condition, then ETA condition, then either pt or status
|
---|
| 284 | if( (pid <= pB || pid == pGLUON) &&// is it a light quark or a gluon, i.e. is it one of these : u,d,c,s,b,g ?
|
---|
| 285 | eta < DET->MAX_TRACKER &&
|
---|
| 286 | particle->Status != 1 &&
|
---|
| 287 | particle->PT > DET->PT_QUARKS_MIN ) {
|
---|
| 288 | NFCentralQ.Add(particle);
|
---|
| 289 | }
|
---|
| 290 |
|
---|
| 291 |
|
---|
| 292 | // keeps only final particles, visible by the central detector, including the fiducial volume
|
---|
| 293 | // the ordering of conditions have been optimised for speed : put first the STATUS condition
|
---|
| 294 | if( (particle->Status == 1) &&
|
---|
| 295 | (
|
---|
| 296 | (pid == pMU && eta < DET->MAX_MU) ||
|
---|
| 297 | (pid != pMU && (pid != pNU1) && (pid != pNU2) && (pid != pNU3) && eta < DET->MAX_CALO_FWD)
|
---|
| 298 | )
|
---|
| 299 | ) {
|
---|
| 300 | switch(pid) {
|
---|
| 301 |
|
---|
| 302 | case pE: // all electrons with eta < DET->MAX_CALO_FWD
|
---|
| 303 | DET->SmearElectron(genMomentum);
|
---|
[29] | 304 | electron.Add(particle);
|
---|
[2] | 305 | break; // case pE
|
---|
| 306 | case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD
|
---|
| 307 | DET->SmearElectron(genMomentum);
|
---|
| 308 | if(genMomentum.E()!=0 && eta < DET->MAX_TRACKER) {
|
---|
| 309 | elementPhoton = (TRootPhoton*) branchPhoton->NewEntry();
|
---|
| 310 | elementPhoton->Set(genMomentum);
|
---|
| 311 | }
|
---|
| 312 | break; // case pGAMMA
|
---|
| 313 | case pMU: // all muons with eta < DET->MAX_MU
|
---|
| 314 | DET->SmearMu(genMomentum);
|
---|
[29] | 315 | muon.Add(particle);
|
---|
[2] | 316 | break; // case pMU
|
---|
| 317 | case pLAMBDA: // all lambdas with eta < DET->MAX_CALO_FWD
|
---|
| 318 | case pK0S: // all K0s with eta < DET->MAX_CALO_FWD
|
---|
| 319 | DET->SmearHadron(genMomentum, 0.7);
|
---|
| 320 | break; // case hadron
|
---|
| 321 | default: // all other final particles with eta < DET->MAX_CALO_FWD
|
---|
| 322 | DET->SmearHadron(genMomentum, 1.0);
|
---|
| 323 | break;
|
---|
| 324 | } // switch (pid)
|
---|
[11] | 325 |
|
---|
[2] | 326 | // all final particles but muons and neutrinos
|
---|
| 327 | // for calorimetric towers and mission PT
|
---|
| 328 | if(genMomentum.E()!=0) {
|
---|
[11] | 329 | PTmis = PTmis + genMomentum;//ptmis
|
---|
| 330 | if(pid !=pMU) {
|
---|
| 331 | towers.push_back(PhysicsTower(LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())));
|
---|
| 332 | // create a fastjet::PseudoJet with these components and put it onto
|
---|
| 333 | // back of the input_particles vector
|
---|
| 334 | input_particles.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
|
---|
| 335 | elementCalo = (TRootCalo*) branchCalo->NewEntry();
|
---|
| 336 | elementCalo->Set(genMomentum);
|
---|
| 337 | }
|
---|
[2] | 338 | }
|
---|
[11] | 339 |
|
---|
| 340 |
|
---|
[2] | 341 | // all final charged particles
|
---|
| 342 | if(
|
---|
[11] | 343 | ((rand()%100) < DET->TRACKING_EFF) &&
|
---|
| 344 | (genMomentum.E()!=0) &&
|
---|
| 345 | (fabs(particle->Eta) < DET->MAX_TRACKER) &&
|
---|
| 346 | (genMomentum.Pt() > DET->PT_TRACKS_MIN ) && // pt too small to be taken into account
|
---|
| 347 | (pid != pGAMMA) &&
|
---|
| 348 | (pid != pPI0) &&
|
---|
| 349 | (pid != pK0L) &&
|
---|
| 350 | (pid != pN) &&
|
---|
| 351 | (pid != pSIGMA0) &&
|
---|
| 352 | (pid != pDELTA0) &&
|
---|
| 353 | (pid != pK0S) // not charged particles : invisible by tracker
|
---|
| 354 | )
|
---|
| 355 | {
|
---|
| 356 | elementTracks = (TRootTracks*) branchTracks->NewEntry();
|
---|
| 357 | elementTracks->Set(genMomentum);
|
---|
| 358 | TrackCentral.push_back(genMomentum);
|
---|
| 359 | }
|
---|
| 360 | } // switch
|
---|
[29] | 361 |
|
---|
| 362 | //branchElectron->Reset();
|
---|
| 363 |
|
---|
| 364 | if(electron.GetEntries()>0) {
|
---|
| 365 | for(int i=0; i < electron.GetEntries();i++) {
|
---|
| 366 | genMomentum.SetPxPyPzE(electron[i]->Px,electron[i]->Py,electron[i]->Pz,electron[i]->E);
|
---|
| 367 | if(genMomentum.E()!=0 && fabs(genMomentum.Eta()) < DET->MAX_TRACKER) {
|
---|
| 368 | elementElec = (TRootElectron*) branchElectron->NewEntry();
|
---|
| 369 | elementElec->Set(genMomentum);
|
---|
| 370 | elementElec->Charge = sign(electron[i]->PID);
|
---|
| 371 | if(fabs(genMomentum.Eta()))elementElec->IsolFlag = DET->Isolation(genMomentum.Phi(),genMomentum.Eta(),TrackCentral,2.0);
|
---|
| 372 | }
|
---|
| 373 | }}
|
---|
| 374 | if(muon.GetEntries()>0) {
|
---|
| 375 | for(int i=0; i < muon.GetEntries();i++) {
|
---|
| 376 | genMomentum.SetPxPyPzE(muon[i]->Px,muon[i]->Py,muon[i]->Pz,muon[i]->E);
|
---|
| 377 | if(genMomentum.E()!=0 && fabs(genMomentum.Eta()) < DET->MAX_MU) {
|
---|
| 378 | elementMu = (TRootMuon*) branchMuon->NewEntry();
|
---|
| 379 | elementMu->Charge = sign(particle->PID);
|
---|
| 380 | elementMu->Set(genMomentum);
|
---|
| 381 | elementMu->IsolFlag = DET->Isolation(genMomentum.Phi(),genMomentum.Eta(),TrackCentral,2.0);
|
---|
| 382 | }
|
---|
| 383 | }}
|
---|
| 384 |
|
---|
[2] | 385 | // Forward particles in CASTOR ?
|
---|
[11] | 386 | /* if (particle->Status == 1 && (fabs(particle->Eta) > DET->MIN_CALO_VFWD)
|
---|
| 387 | && (fabs(particle->Eta) < DET->MAX_CALO_VFWD)) {
|
---|
| 388 |
|
---|
| 389 |
|
---|
| 390 | } // CASTOR
|
---|
| 391 | */
|
---|
[2] | 392 | // Zero degree calorimeter, for forward neutrons and photons
|
---|
| 393 | if (particle->Status ==1 && (pid == pN || pid == pGAMMA ) && eta > DET->MIN_ZDC ) {
|
---|
[11] | 394 | // !!!!!!!!! vérifier que particle->Z est bien en micromÚtres!!!
|
---|
| 395 | // !!!!!!!!! vérifier que particle->T est bien en secondes!!!
|
---|
| 396 | // !!!!!!!!! pas de smearing ! on garde trop d'info !
|
---|
| 397 | elementZdc = (TRootZdcHits*) branchZDC->NewEntry();
|
---|
| 398 | elementZdc->Set(genMomentum);
|
---|
| 399 |
|
---|
| 400 | // time of flight t is t = T + d/[ cos(theta) v ]
|
---|
| 401 | //double tx = acos(particle->Px/particle->Pz);
|
---|
| 402 | //double ty = acos(particle->Py/particle->Pz);
|
---|
| 403 | //double theta = (1E-6)*sqrt( pow(tx,2) + pow(ty,2) );
|
---|
| 404 | //double flight_distance = (DET->ZDC_S - particle->Z*(1E-6))/cos(theta) ; // assumes that Z is in micrometers
|
---|
| 405 | double flight_distance = (DET->ZDC_S - particle->Z*(1E-6));
|
---|
| 406 | // assumes also that the emission angle is so small that 1/(cos theta) = 1
|
---|
| 407 | elementZdc->T = 0*particle->T + flight_distance/speed_of_light; // assumes highly relativistic particles
|
---|
| 408 | elementZdc->side = sign(particle->Eta);
|
---|
| 409 |
|
---|
[2] | 410 | }
|
---|
[11] | 411 |
|
---|
[2] | 412 | // if forward proton
|
---|
| 413 | if( (pid == pP) && (particle->Status == 1) && (fabs(particle->Eta) > DET->MAX_CALO_FWD) )
|
---|
| 414 | {
|
---|
[11] | 415 | // !!!!!!!! put here particle->CHARGE and particle->MASS
|
---|
| 416 | H_BeamParticle p1; /// put here particle->CHARGE and particle->MASS
|
---|
| 417 | p1.smearAng();
|
---|
| 418 | p1.smearPos();
|
---|
| 419 | p1.setPosition(p1.getX()-500.,p1.getY(),p1.getTX()-1*kickers_on*CRANG,p1.getTY(),0);
|
---|
| 420 | p1.set4Momentum(particle->Px,particle->Py,particle->Pz,particle->E);
|
---|
| 421 |
|
---|
| 422 | H_BeamLine *beamline;
|
---|
| 423 | if(particle->Eta >0) beamline = beamline1;
|
---|
| 424 | else beamline = beamline2;
|
---|
| 425 |
|
---|
| 426 | p1.computePath(beamline,1);
|
---|
| 427 |
|
---|
| 428 | if(p1.stopped(beamline)) {
|
---|
| 429 | if (p1.getStoppingElement()->getName()=="rp220_1" || p1.getStoppingElement()->getName()=="rp220_2") {
|
---|
| 430 | p1.propagate(DET->RP220_S);
|
---|
| 431 | elementRP220 = (TRootRomanPotHits*) branchRP220->NewEntry();
|
---|
| 432 | elementRP220->X = (1E-6)*p1.getX(); // [m]
|
---|
| 433 | elementRP220->Y = (1E-6)*p1.getY(); // [m]
|
---|
| 434 | elementRP220->Tx = (1E-6)*p1.getTX(); // [rad]
|
---|
| 435 | elementRP220->Ty = (1E-6)*p1.getTY(); // [rad]
|
---|
| 436 | elementRP220->S = p1.getS(); // [m]
|
---|
| 437 | elementRP220->T = -1; // not yet implemented
|
---|
| 438 | elementRP220->E = p1.getE(); // not yet implemented
|
---|
| 439 | elementRP220->q2 = -1; // not yet implemented
|
---|
| 440 | elementRP220->side = sign(particle->Eta);
|
---|
| 441 |
|
---|
| 442 | } else if (p1.getStoppingElement()->getName()=="rp420_1" || p1.getStoppingElement()->getName()=="rp420_2") {
|
---|
| 443 | p1.propagate(DET->FP420_S);
|
---|
| 444 | elementFP420 = (TRootRomanPotHits*) branchFP420->NewEntry();
|
---|
| 445 | elementFP420->X = (1E-6)*p1.getX(); // [m]
|
---|
| 446 | elementFP420->Y = (1E-6)*p1.getY(); // [m]
|
---|
| 447 | elementFP420->Tx = (1E-6)*p1.getTX(); // [rad]
|
---|
| 448 | elementFP420->Ty = (1E-6)*p1.getTY(); // [rad]
|
---|
| 449 | elementFP420->S = p1.getS(); // [m]
|
---|
| 450 | elementFP420->T = -1; // not yet implemented
|
---|
| 451 | elementFP420->E = p1.getE(); // not yet implemented
|
---|
| 452 | elementFP420->q2 = -1; // not yet implemented
|
---|
| 453 | elementFP420->side = sign(particle->Eta);
|
---|
[2] | 454 | }
|
---|
[11] | 455 | }
|
---|
| 456 |
|
---|
| 457 | // if(p1.stopped(beamline) && (p1.getStoppingElement()->getS() > 100))
|
---|
| 458 | // cout << "Eloss =" << 7000.-p1.getE() << " ; " << p1.getStoppingElement()->getName() << endl;
|
---|
[2] | 459 | } // if forward proton
|
---|
[11] | 460 |
|
---|
[2] | 461 | } // while
|
---|
[11] | 462 |
|
---|
[2] | 463 | // computes the Missing Transverse Momentum
|
---|
| 464 | elementEtmis = (TRootETmis*) branchETmis->NewEntry();
|
---|
| 465 | elementEtmis->ET = (-PTmis).Pt();
|
---|
| 466 | elementEtmis->Phi = (-PTmis).Phi();
|
---|
| 467 | elementEtmis->Px = (-PTmis).Px();
|
---|
| 468 | elementEtmis->Py = (-PTmis).Py();
|
---|
[11] | 469 |
|
---|
| 470 | //*****************************
|
---|
[2] | 471 |
|
---|
[11] | 472 | // run the jet clustering with the above jet definition
|
---|
| 473 | if(input_particles.size()!=0)
|
---|
| 474 | {
|
---|
| 475 | fastjet::ClusterSequence clust_seq(input_particles, jet_def);
|
---|
| 476 |
|
---|
| 477 |
|
---|
| 478 | // extract the inclusive jets with pt > 5 GeV
|
---|
| 479 | double ptmin = 5.0;
|
---|
| 480 | inclusive_jets = clust_seq.inclusive_jets(ptmin);
|
---|
| 481 |
|
---|
| 482 | // sort jets into increasing pt
|
---|
[15] | 483 | sorted_jets = sorted_by_pt(inclusive_jets);
|
---|
[11] | 484 | }
|
---|
[15] | 485 | for (unsigned int i = 0; i < sorted_jets.size(); i++) {
|
---|
[11] | 486 | elementJet = (TRootJet*) branchJet->NewEntry();
|
---|
| 487 | TLorentzVector JET;
|
---|
[15] | 488 | JET.SetPxPyPzE(sorted_jets[i].px(),sorted_jets[i].py(),sorted_jets[i].pz(),sorted_jets[i].E());
|
---|
[21] | 489 | //cout<<"Jet.Pt() "<<JET.Pt()<<endl;
|
---|
[11] | 490 | elementJet->Set(JET);
|
---|
| 491 | // b-jets
|
---|
| 492 | bool btag=false;
|
---|
| 493 | if((fabs(JET.Eta()) < DET->MAX_TRACKER && DET->Btaggedjet(JET, NFCentralQ)))btag=true;
|
---|
| 494 | elementJet->Btag = btag;
|
---|
| 495 |
|
---|
| 496 | // Tau jet identification : 1! track and electromagnetic collimation
|
---|
| 497 | if(fabs(JET.Eta()) < (DET->MAX_TRACKER - DET->TAU_CONE_TRACKS)) {
|
---|
| 498 | double Energie_tau_central = DET->EnergySmallCone(towers,JET.Eta(),JET.Phi());
|
---|
| 499 | if(
|
---|
| 500 | ( Energie_tau_central/JET.E() > DET->TAU_EM_COLLIMATION ) &&
|
---|
| 501 | ( DET->NumTracks(TrackCentral,DET->PT_TRACK_TAU,JET.Eta(),JET.Phi()) == 1 )
|
---|
| 502 | ) {
|
---|
| 503 | elementTauJet = (TRootTauJet*) branchTauJet->NewEntry();
|
---|
| 504 | elementTauJet->Set(JET);
|
---|
| 505 | } // if tau jet
|
---|
| 506 | } // if JET.eta < tracker - tau_cone : Tau jet identification
|
---|
| 507 | } // for itJet : loop on all jets
|
---|
| 508 |
|
---|
[2] | 509 | treeWriter->Fill();
|
---|
| 510 | // Add here the trigger
|
---|
| 511 | // Should test all the trigger table on the event, based on reconstructed objects
|
---|
| 512 | } // Loop over all events
|
---|
| 513 | treeWriter->Write();
|
---|
| 514 |
|
---|
| 515 | cout << "** Exiting..." << endl;
|
---|
| 516 |
|
---|
| 517 | delete treeWriter;
|
---|
| 518 | delete treeReader;
|
---|
| 519 | delete DET;
|
---|
| 520 | if(converter) delete converter;
|
---|
| 521 |
|
---|
| 522 | todo("TODO");
|
---|
| 523 | }
|
---|
| 524 |
|
---|