| 1 | /***********************************************************************
|
|---|
| 2 | ** **
|
|---|
| 3 | ** /----------------------------------------------\ **
|
|---|
| 4 | ** | Delphes, a framework for the fast simulation | **
|
|---|
| 5 | ** | of a generic collider experiment | **
|
|---|
| 6 | ** \------------- arXiv:0903.2225v1 ------------/ **
|
|---|
| 7 | ** **
|
|---|
| 8 | ** **
|
|---|
| 9 | ** This package uses: **
|
|---|
| 10 | ** ------------------ **
|
|---|
| 11 | ** ROOT: Nucl. Inst. & Meth. in Phys. Res. A389 (1997) 81-86 **
|
|---|
| 12 | ** FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210] **
|
|---|
| 13 | ** Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2] **
|
|---|
| 14 | ** FROG: [hep-ex/0901.2718v1] **
|
|---|
| 15 | ** HepMC: Comput. Phys. Commun.134 (2001) 41 **
|
|---|
| 16 | ** **
|
|---|
| 17 | ** ------------------------------------------------------------------ **
|
|---|
| 18 | ** **
|
|---|
| 19 | ** Main authors: **
|
|---|
| 20 | ** ------------- **
|
|---|
| 21 | ** **
|
|---|
| 22 | ** Severine Ovyn Xavier Rouby **
|
|---|
| 23 | ** severine.ovyn@uclouvain.be xavier.rouby@cern **
|
|---|
| 24 | ** **
|
|---|
| 25 | ** Center for Particle Physics and Phenomenology (CP3) **
|
|---|
| 26 | ** Universite catholique de Louvain (UCL) **
|
|---|
| 27 | ** Louvain-la-Neuve, Belgium **
|
|---|
| 28 | ** **
|
|---|
| 29 | ** Copyright (C) 2008-2009, **
|
|---|
| 30 | ** All rights reserved. **
|
|---|
| 31 | ** **
|
|---|
| 32 | ***********************************************************************/
|
|---|
| 33 |
|
|---|
| 34 | /// \file Resolution_ATLAS.cpp
|
|---|
| 35 | /// \brief Resolution for Atlas
|
|---|
| 36 |
|
|---|
| 37 | #include "TChain.h"
|
|---|
| 38 | #include "TApplication.h"
|
|---|
| 39 | #include "TFile.h"
|
|---|
| 40 |
|
|---|
| 41 | #include "ExRootTreeReader.h"
|
|---|
| 42 | #include "ExRootTreeWriter.h"
|
|---|
| 43 | #include "ExRootTreeBranch.h"
|
|---|
| 44 | #include "TreeClasses.h"
|
|---|
| 45 |
|
|---|
| 46 | #include "SmearUtil.h"
|
|---|
| 47 | #include "JetsUtil.h"
|
|---|
| 48 | #include "BFieldProp.h"
|
|---|
| 49 |
|
|---|
| 50 | #include<vector>
|
|---|
| 51 | #include<iostream>
|
|---|
| 52 |
|
|---|
| 53 | using namespace std;
|
|---|
| 54 |
|
|---|
| 55 | //------------------------------------------------------------------------------
|
|---|
| 56 |
|
|---|
| 57 | // //********************************** PYTHIA INFORMATION*********************************
|
|---|
| 58 |
|
|---|
| 59 | TSimpleArray<TRootGenParticle> TauHadr(const TClonesArray *GEN)
|
|---|
| 60 | {
|
|---|
| 61 | TIter it((TCollection*)GEN);
|
|---|
| 62 | it.Reset();
|
|---|
| 63 | TRootGenParticle *gen1;
|
|---|
| 64 | TSimpleArray<TRootGenParticle> array,array2;
|
|---|
| 65 |
|
|---|
| 66 | while((gen1 = (TRootGenParticle*) it.Next()))
|
|---|
| 67 | {
|
|---|
| 68 | array.Add(gen1);
|
|---|
| 69 | }
|
|---|
| 70 | it.Reset();
|
|---|
| 71 | bool tauhad;
|
|---|
| 72 | while((gen1 = (TRootGenParticle*) it.Next()))
|
|---|
| 73 | {
|
|---|
| 74 | tauhad=false;
|
|---|
| 75 | if(abs(gen1->PID)==15)
|
|---|
| 76 | {
|
|---|
| 77 | int d1=gen1->D1;
|
|---|
| 78 | int d2=gen1->D2;
|
|---|
| 79 | if((d1 < array.GetEntries()) && (d1 > 0) && (d2 < array.GetEntries()) && (d2 > 0))
|
|---|
| 80 | {
|
|---|
| 81 | tauhad=true;
|
|---|
| 82 | for(int d=d1; d < d2+1; d++)
|
|---|
| 83 | {
|
|---|
| 84 | if(abs(array[d]->PID)== pE || abs(array[d]->PID)== pMU)tauhad=false;
|
|---|
| 85 | }
|
|---|
| 86 | }
|
|---|
| 87 | }
|
|---|
| 88 | if(tauhad)array2.Add(gen1);
|
|---|
| 89 | }
|
|---|
| 90 | return array2;
|
|---|
| 91 | }
|
|---|
| 92 |
|
|---|
| 93 | double EnergySmallCone(const vector<TLorentzVector> &towers, const float eta, const float phi,float energy_scone,float JET_seed) {
|
|---|
| 94 | double Energie=0;
|
|---|
| 95 | for(unsigned int i=0; i < towers.size(); i++) {
|
|---|
| 96 | if(towers[i].Pt() < JET_seed) continue;
|
|---|
| 97 | if((DeltaR(phi,eta,towers[i].Phi(),towers[i].Eta()) < energy_scone)) {
|
|---|
| 98 | Energie += towers[i].E();
|
|---|
| 99 | }
|
|---|
| 100 | }
|
|---|
| 101 | return Energie;
|
|---|
| 102 | }
|
|---|
| 103 |
|
|---|
| 104 |
|
|---|
| 105 | void PairingJet(TLorentzVector &JETSm, const TLorentzVector &JET, const TClonesArray *branchJet, const float dR=0.25)
|
|---|
| 106 | {
|
|---|
| 107 | JETSm.SetPxPyPzE(0,0,0,0);
|
|---|
| 108 | float deltaRtest=5000;
|
|---|
| 109 | TIter itJet((TCollection*)branchJet);
|
|---|
| 110 | TRootJet *jet;
|
|---|
| 111 | itJet.Reset();
|
|---|
| 112 | while( (jet = (TRootJet*) itJet.Next()) )
|
|---|
| 113 | {
|
|---|
| 114 | TLorentzVector Att;
|
|---|
| 115 | Att.SetPtEtaPhiE(jet->PT,jet->Eta,jet->Phi,jet->E);
|
|---|
| 116 | if(DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta()) < deltaRtest)
|
|---|
| 117 | {
|
|---|
| 118 | deltaRtest = DeltaR(JET.Phi(),JET.Eta(),Att.Phi(),Att.Eta());
|
|---|
| 119 | if(deltaRtest < dR)
|
|---|
| 120 | {
|
|---|
| 121 | JETSm = Att;
|
|---|
| 122 | }
|
|---|
| 123 | }
|
|---|
| 124 | }
|
|---|
| 125 | }
|
|---|
| 126 |
|
|---|
| 127 | void PairingElec(TLorentzVector &ELECSm, const TLorentzVector &ELEC, const TClonesArray *branchElec)
|
|---|
| 128 | {
|
|---|
| 129 | ELECSm.SetPxPyPzE(0,0,0,0);
|
|---|
| 130 | float deltaRtest=5000;
|
|---|
| 131 | TIter itElec((TCollection*)branchElec);
|
|---|
| 132 | TRootElectron *elec;
|
|---|
| 133 | itElec.Reset();
|
|---|
| 134 | while( (elec = (TRootElectron*) itElec.Next()) )
|
|---|
| 135 | {
|
|---|
| 136 | TLorentzVector Att;
|
|---|
| 137 | Att.SetPtEtaPhiE(elec->PT,elec->Eta,elec->Phi,elec->E);
|
|---|
| 138 | if(DeltaR(ELEC.Phi(),ELEC.Eta(),Att.Phi(),Att.Eta()) < deltaRtest)
|
|---|
| 139 | {
|
|---|
| 140 | deltaRtest = DeltaR(ELEC.Phi(),ELEC.Eta(),Att.Phi(),Att.Eta());
|
|---|
| 141 | if(deltaRtest < 0.025)
|
|---|
| 142 | {
|
|---|
| 143 | ELECSm = Att;
|
|---|
| 144 | }
|
|---|
| 145 | }
|
|---|
| 146 | }
|
|---|
| 147 | }
|
|---|
| 148 |
|
|---|
| 149 | void PairingMuon(TLorentzVector &MUONSm, const TLorentzVector &MUON, const TClonesArray *branchMuon)
|
|---|
| 150 | {
|
|---|
| 151 | MUONSm.SetPxPyPzE(0,0,0,0);
|
|---|
| 152 | float deltaRtest=5000;
|
|---|
| 153 | TIter itMuon((TCollection*)branchMuon);
|
|---|
| 154 | TRootMuon *muon;
|
|---|
| 155 | itMuon.Reset();
|
|---|
| 156 | while( (muon = (TRootMuon*) itMuon.Next()) )
|
|---|
| 157 | {
|
|---|
| 158 | TLorentzVector Att;
|
|---|
| 159 | Att.SetPxPyPzE(muon->Px,muon->Py,muon->Pz,muon->E);
|
|---|
| 160 | if(DeltaR(MUON.Phi(),MUON.Eta(),Att.Phi(),Att.Eta()) < deltaRtest)
|
|---|
| 161 | {
|
|---|
| 162 | deltaRtest = DeltaR(MUON.Phi(),MUON.Eta(),Att.Phi(),Att.Eta());
|
|---|
| 163 | if(deltaRtest < 0.025)
|
|---|
| 164 | {
|
|---|
| 165 | MUONSm = Att;
|
|---|
| 166 | }
|
|---|
| 167 | }
|
|---|
| 168 | }
|
|---|
| 169 | }
|
|---|
| 170 |
|
|---|
| 171 | unsigned int NumTracks(const TClonesArray *branchTracks, const float pt_track, const float eta, const float phi,float track_scone) {
|
|---|
| 172 | unsigned int numtrack=0;
|
|---|
| 173 | TIter itTrack((TCollection*)branchTracks);
|
|---|
| 174 | TRootTracks *track;
|
|---|
| 175 | itTrack.Reset();
|
|---|
| 176 | while( (track = (TRootTracks*) itTrack.Next()) )
|
|---|
| 177 | {
|
|---|
| 178 | if((track->PT < pt_track )||
|
|---|
| 179 | (DeltaR(phi,eta,track->Phi,track->Eta) > track_scone)
|
|---|
| 180 | )continue;
|
|---|
| 181 | numtrack++;
|
|---|
| 182 | }
|
|---|
| 183 | return numtrack;
|
|---|
| 184 | }
|
|---|
| 185 |
|
|---|
| 186 |
|
|---|
| 187 |
|
|---|
| 188 | int main(int argc, char *argv[])
|
|---|
| 189 | {
|
|---|
| 190 | int appargc = 2;
|
|---|
| 191 | char *appName = "Resolution";
|
|---|
| 192 | char *appargv[] = {appName, "-b"};
|
|---|
| 193 | TApplication app(appName, &appargc, appargv);
|
|---|
| 194 |
|
|---|
| 195 | if(argc != 3) {
|
|---|
| 196 | cout << " Usage: " << argv[0] << " input_file" << " output_file" << endl;
|
|---|
| 197 | cout << " input_file - input file in Delphes-root format," << endl;
|
|---|
| 198 | cout << " output_file - output file." << endl;
|
|---|
| 199 | exit(1);
|
|---|
| 200 | }
|
|---|
| 201 |
|
|---|
| 202 | srand (time (NULL)); /* Initialisation du générateur */
|
|---|
| 203 |
|
|---|
| 204 | //read the input TROOT file
|
|---|
| 205 | string inputfilename(argv[1]), outputfilename(argv[2]);
|
|---|
| 206 |
|
|---|
| 207 | if(outputfilename.find(".root") > outputfilename.length() ) {
|
|---|
| 208 | cout << "output_file should be a .root file!\n";
|
|---|
| 209 | return -1;
|
|---|
| 210 | }
|
|---|
| 211 |
|
|---|
| 212 |
|
|---|
| 213 |
|
|---|
| 214 | TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE");// Creates the file, but should be closed just after
|
|---|
| 215 | outputFile->Close();
|
|---|
| 216 |
|
|---|
| 217 | TChain chainGEN("GEN");
|
|---|
| 218 | chainGEN.Add(inputfilename.c_str());
|
|---|
| 219 | ExRootTreeReader *treeReaderGEN = new ExRootTreeReader(&chainGEN);
|
|---|
| 220 | TChain chain("Analysis");
|
|---|
| 221 | chain.Add(inputfilename.c_str());
|
|---|
| 222 | ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
|
|---|
| 223 | const TClonesArray *branchJet = treeReader->UseBranch("Jet");
|
|---|
| 224 | const TClonesArray *branchElec = treeReader->UseBranch("Electron");
|
|---|
| 225 | const TClonesArray *branchMuon = treeReader->UseBranch("Muon");
|
|---|
| 226 | const TClonesArray *branchTracks = treeReader->UseBranch("Tracks");
|
|---|
| 227 | const TClonesArray *branchTowers = treeReader->UseBranch("CaloTower");
|
|---|
| 228 | const TClonesArray *branchGen = treeReaderGEN->UseBranch("Particle");
|
|---|
| 229 | TIter itGen((TCollection*)branchGen);
|
|---|
| 230 |
|
|---|
| 231 | //write the output root file
|
|---|
| 232 | ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputfilename, "Analysis");
|
|---|
| 233 | ExRootTreeBranch *branchjet = treeWriter->NewBranch("JetPTResol", RESOLJET::Class());
|
|---|
| 234 | ExRootTreeBranch *branchelec = treeWriter->NewBranch("ElecEResol", RESOLELEC::Class());
|
|---|
| 235 | ExRootTreeBranch *branchmuon = treeWriter->NewBranch("MuonPTResol", RESOLMUON::Class());
|
|---|
| 236 | ExRootTreeBranch *branchtaujet = treeWriter->NewBranch("TauJetPTResol", TAUHAD::Class());
|
|---|
| 237 | ExRootTreeBranch *branchetmis = treeWriter->NewBranch("ETmisResol",ETMIS::Class());
|
|---|
| 238 |
|
|---|
| 239 | TRootGenParticle *particle;
|
|---|
| 240 |
|
|---|
| 241 | RESOLELEC * elementElec;
|
|---|
| 242 | RESOLMUON *elementMuon;
|
|---|
| 243 | RESOLJET *elementJet;
|
|---|
| 244 | TAUHAD *elementTaujet;
|
|---|
| 245 | ETMIS *elementEtmis;
|
|---|
| 246 |
|
|---|
| 247 | int numTau=0;
|
|---|
| 248 | int numTauRec=0;
|
|---|
| 249 |
|
|---|
| 250 | RESOLution *DET = new RESOLution();
|
|---|
| 251 | /*
|
|---|
| 252 | string detectorcard = "data/DetectorCard_CMS.dat";
|
|---|
| 253 | const float dR_jetpairing = 0.25;
|
|---|
| 254 | const float jet_pt_cut = 1;
|
|---|
| 255 | */
|
|---|
| 256 | string detectorcard = "data/DetectorCard_ATLAS.dat";
|
|---|
| 257 | const float dR_jetpairing = 0.2;
|
|---|
| 258 | const float jet_pt_cut = 7;
|
|---|
| 259 | DET->ReadDataCard(detectorcard);
|
|---|
| 260 |
|
|---|
| 261 |
|
|---|
| 262 | //Jet information
|
|---|
| 263 | JetsUtil *JETRUN = new JetsUtil(detectorcard);
|
|---|
| 264 |
|
|---|
| 265 | TLorentzVector genMomentum(0,0,0,0);//TLorentzVector containing generator level information
|
|---|
| 266 | TLorentzVector recoMomentum(0,0,0,0);//TLorentzVector containing generator level information
|
|---|
| 267 | LorentzVector jetMomentum;
|
|---|
| 268 |
|
|---|
| 269 | vector<fastjet::PseudoJet> input_particlesGEN;//for FastJet algorithm
|
|---|
| 270 | vector<fastjet::PseudoJet> sorted_jetsGEN;
|
|---|
| 271 | vector<int> NTrackJet;
|
|---|
| 272 | vector<TLorentzVector> towers;
|
|---|
| 273 |
|
|---|
| 274 | // Loop over all events
|
|---|
| 275 | Long64_t entry, allEntries = treeReader->GetEntries();
|
|---|
| 276 | cout << "** Chain contains " << allEntries << " events" << endl;
|
|---|
| 277 | for(entry = 0; entry < allEntries; ++entry)
|
|---|
| 278 | {
|
|---|
| 279 | TLorentzVector PTmisReco(0,0,0,0);
|
|---|
| 280 | TLorentzVector PTmisGEN(0,0,0,0);
|
|---|
| 281 | treeReader->ReadEntry(entry);
|
|---|
| 282 | treeReaderGEN->ReadEntry(entry);
|
|---|
| 283 | treeWriter->Clear();
|
|---|
| 284 | if((entry % 100) == 0 && entry > 0 ) cout << "** Processing element # " << entry << endl;
|
|---|
| 285 |
|
|---|
| 286 | TSimpleArray<TRootGenParticle> bGen;
|
|---|
| 287 | itGen.Reset();
|
|---|
| 288 | TSimpleArray<TRootGenParticle> NFCentralQ;
|
|---|
| 289 |
|
|---|
| 290 | input_particlesGEN.clear();
|
|---|
| 291 | towers.clear();
|
|---|
| 292 |
|
|---|
| 293 | // Loop over all particles in event
|
|---|
| 294 | while( (particle = (TRootGenParticle*) itGen.Next()) )
|
|---|
| 295 | {
|
|---|
| 296 | genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
|
|---|
| 297 |
|
|---|
| 298 | int pid = abs(particle->PID);
|
|---|
| 299 | float eta = fabs(particle->Eta);
|
|---|
| 300 |
|
|---|
| 301 | //input generator level particle for jet algorithm
|
|---|
| 302 | if(particle->Status == 1 && eta < DET->CEN_max_calo_fwd)
|
|---|
| 303 | {
|
|---|
| 304 | input_particlesGEN.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
|
|---|
| 305 | }
|
|---|
| 306 |
|
|---|
| 307 | //Calculate ETMIS from generated particles
|
|---|
| 308 | if((pid == pNU1) || (pid == pNU2) || (pid == pNU3))PTmisGEN = PTmisGEN + genMomentum;
|
|---|
| 309 |
|
|---|
| 310 | //Electrons and muons
|
|---|
| 311 | if( (particle->Status == 1) &&
|
|---|
| 312 | ((pid != pNU1) && (pid != pNU2) && (pid != pNU3)) &&
|
|---|
| 313 | (fabs(particle->Eta) < DET->CEN_max_calo_fwd)
|
|---|
| 314 | )
|
|---|
| 315 | {
|
|---|
| 316 | eta=fabs(genMomentum.Eta());
|
|---|
| 317 |
|
|---|
| 318 | switch(pid) {
|
|---|
| 319 |
|
|---|
| 320 | case pE: // all electrons with eta < DET->MAX_CALO_FWD
|
|---|
| 321 | PairingElec(recoMomentum,genMomentum,branchElec);
|
|---|
| 322 | if(recoMomentum.E()!=0){
|
|---|
| 323 | elementElec=(RESOLELEC*) branchelec->NewEntry();
|
|---|
| 324 | elementElec->E = genMomentum.E();
|
|---|
| 325 | elementElec->SmearedE = recoMomentum.E();}
|
|---|
| 326 | break; // case pE
|
|---|
| 327 | case pMU: // all muons with eta < DET->MAX_MU
|
|---|
| 328 | PairingMuon(recoMomentum,genMomentum,branchMuon);
|
|---|
| 329 | if(recoMomentum.E()!=0){
|
|---|
| 330 | elementMuon = (RESOLMUON*) branchmuon->NewEntry();
|
|---|
| 331 | elementMuon->OverPT = 1./genMomentum.Pt();
|
|---|
| 332 | elementMuon->OverSmearedPT = 1./recoMomentum.Pt();}
|
|---|
| 333 | break; // case pMU
|
|---|
| 334 | default:
|
|---|
| 335 | break;
|
|---|
| 336 | } // switch (pid)
|
|---|
| 337 | }
|
|---|
| 338 |
|
|---|
| 339 | } // while
|
|---|
| 340 |
|
|---|
| 341 | //compute missing transverse energy from calo towers
|
|---|
| 342 | TIter itCalo((TCollection*)branchTowers);
|
|---|
| 343 | TRootCalo *calo;
|
|---|
| 344 | itCalo.Reset();
|
|---|
| 345 | TLorentzVector Att(0.,0.,0.,0.);
|
|---|
| 346 | float ScalarEt=0;
|
|---|
| 347 | while( (calo = (TRootCalo*) itCalo.Next()) )
|
|---|
| 348 | {
|
|---|
| 349 | if(calo->E !=0){
|
|---|
| 350 | Att.SetPtEtaPhiE(calo->getET(),calo->Eta,calo->Phi,calo->E);
|
|---|
| 351 | towers.push_back(Att);
|
|---|
| 352 | if(fabs(Att.Eta()) < DET->CEN_max_calo_fwd)
|
|---|
| 353 | {
|
|---|
| 354 | ScalarEt = ScalarEt + calo->getET();
|
|---|
| 355 | PTmisReco = PTmisReco + Att;
|
|---|
| 356 | }
|
|---|
| 357 | }
|
|---|
| 358 | }
|
|---|
| 359 | elementEtmis= (ETMIS*) branchetmis->NewEntry();
|
|---|
| 360 | elementEtmis->Et = (PTmisGEN).Pt();
|
|---|
| 361 | elementEtmis->Ex = (-PTmisGEN).Px();
|
|---|
| 362 | elementEtmis->SEt = ScalarEt;
|
|---|
| 363 | elementEtmis->EtSmeare = (PTmisReco).Pt()-(PTmisGEN).Pt();
|
|---|
| 364 | elementEtmis->ExSmeare = (-PTmisReco).Px()-(PTmisGEN).Px();
|
|---|
| 365 |
|
|---|
| 366 | //*****************************
|
|---|
| 367 | sorted_jetsGEN=JETRUN->RunJetsResol(input_particlesGEN);
|
|---|
| 368 |
|
|---|
| 369 | TSimpleArray<TRootGenParticle> TausHadr = TauHadr(branchGen);
|
|---|
| 370 | TLorentzVector JETreco(0,0,0,0);
|
|---|
| 371 | for (unsigned int i = 0; i < sorted_jetsGEN.size(); i++) {
|
|---|
| 372 | TLorentzVector JETgen(0,0,0,0);
|
|---|
| 373 | JETgen.SetPxPyPzE(sorted_jetsGEN[i].px(),sorted_jetsGEN[i].py(),sorted_jetsGEN[i].pz(),sorted_jetsGEN[i].E());
|
|---|
| 374 | PairingJet(JETreco,JETgen,branchJet);
|
|---|
| 375 | if(JETreco.Pt()>jet_pt_cut)
|
|---|
| 376 | {
|
|---|
| 377 | elementJet= (RESOLJET*) branchjet->NewEntry();
|
|---|
| 378 | elementJet->PT = JETgen.Et();
|
|---|
| 379 | elementJet->SmearedPT = JETreco.Et()/JETgen.Et();
|
|---|
| 380 | elementJet->E = JETgen.E();
|
|---|
| 381 | elementJet->dE = (JETreco.E()-JETgen.E())/JETgen.E() ;
|
|---|
| 382 | elementJet->dE2 = pow( (JETreco.E()-JETgen.E())/JETgen.E() , 2.) ;
|
|---|
| 383 | }
|
|---|
| 384 | }
|
|---|
| 385 | numTau = numTau+TausHadr.GetEntries();
|
|---|
| 386 |
|
|---|
| 387 | TIter itJet((TCollection*)branchJet);
|
|---|
| 388 | TRootJet *jet;
|
|---|
| 389 | itJet.Reset();
|
|---|
| 390 | while( (jet = (TRootJet*) itJet.Next()) )
|
|---|
| 391 | {
|
|---|
| 392 | TLorentzVector JETT(0,0,0,0);
|
|---|
| 393 | JETT.SetPxPyPzE(jet->Px,jet->Py,jet->Pz,jet->E);
|
|---|
| 394 | if(fabs(JETT.Eta()) < (DET->CEN_max_tracker - DET->TAU_track_scone))
|
|---|
| 395 | {
|
|---|
| 396 | for(Int_t i=0; i<TausHadr.GetEntries();i++)
|
|---|
| 397 | {
|
|---|
| 398 | if(DeltaR(TausHadr[i]->Phi,TausHadr[i]->Eta,JETT.Phi(),JETT.Eta())<0.1)
|
|---|
| 399 | {
|
|---|
| 400 | elementTaujet= (TAUHAD*) branchtaujet->NewEntry();
|
|---|
| 401 | elementTaujet->EnergieCen = EnergySmallCone(towers,JETT.Eta(),JETT.Phi(),DET->TAU_energy_scone,DET->JET_seed)/JETT.E();
|
|---|
| 402 | elementTaujet->NumTrack = NumTracks(branchTracks,DET->TAU_track_pt,JETT.Eta(),JETT.Phi(),DET->TAU_track_scone);
|
|---|
| 403 | if( (EnergySmallCone(towers,JETT.Eta(),JETT.Phi(),DET->TAU_energy_scone,DET->JET_seed)/JETT.E()) > 0.95
|
|---|
| 404 | && (NumTracks(branchTracks,DET->TAU_track_pt,JETT.Eta(),JETT.Phi(),DET->TAU_track_scone))==1)numTauRec++;
|
|---|
| 405 | }
|
|---|
| 406 | }
|
|---|
| 407 | }
|
|---|
| 408 |
|
|---|
| 409 |
|
|---|
| 410 | } // for itJet : loop on all jets
|
|---|
| 411 |
|
|---|
| 412 | treeWriter->Fill();
|
|---|
| 413 | } // Loop over all events
|
|---|
| 414 | treeWriter->Write();
|
|---|
| 415 |
|
|---|
| 416 | cout << detectorcard << " has been used.\n";
|
|---|
| 417 | cout << "** Exiting..." << endl;
|
|---|
| 418 |
|
|---|
| 419 | delete treeWriter;
|
|---|
| 420 | delete treeReader;
|
|---|
| 421 | delete DET;
|
|---|
| 422 | }
|
|---|
| 423 |
|
|---|
