Changeset 1345 in svn for trunk/modules
- Timestamp:
- Dec 21, 2013, 3:00:11 PM (11 years ago)
- Location:
- trunk/modules
- Files:
-
- 8 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/modules/Calorimeter.cc
r1280 r1345 330 330 fTrackHCalEnergy = 0.0; 331 331 332 fTowerECalTime = 0.0; 333 fTowerHCalTime = 0.0; 334 335 fTrackECalTime = 0.0; 336 fTrackHCalTime = 0.0; 337 338 fTowerECalWeightTime = 0.0; 339 fTowerHCalWeightTime = 0.0; 340 332 341 fTowerTrackHits = 0; 333 342 fTowerPhotonHits = 0; … … 343 352 track = static_cast<Candidate*>(fTrackInputArray->At(number)); 344 353 momentum = track->Momentum; 345 354 position = track->Position; 355 356 346 357 ecalEnergy = momentum.E() * fTrackECalFractions[number]; 347 358 hcalEnergy = momentum.E() * fTrackHCalFractions[number]; … … 349 360 fTrackECalEnergy += ecalEnergy; 350 361 fTrackHCalEnergy += hcalEnergy; 362 363 fTrackECalTime += TMath::Sqrt(ecalEnergy)*position.T(); 364 fTrackHCalTime += TMath::Sqrt(hcalEnergy)*position.T(); 365 366 fTrackECalWeightTime += TMath::Sqrt(ecalEnergy); 367 fTrackHCalWeightTime += TMath::Sqrt(hcalEnergy); 351 368 352 369 fTowerTrackArray->Add(track); … … 360 377 particle = static_cast<Candidate*>(fParticleInputArray->At(number)); 361 378 momentum = particle->Momentum; 379 position = particle->Position; 362 380 363 381 // fill current tower … … 367 385 fTowerECalEnergy += ecalEnergy; 368 386 fTowerHCalEnergy += hcalEnergy; 387 388 fTowerECalTime += TMath::Sqrt(ecalEnergy)*position.T(); 389 fTowerHCalTime += TMath::Sqrt(hcalEnergy)*position.T(); 390 391 fTowerECalWeightTime += TMath::Sqrt(ecalEnergy); 392 fTowerHCalWeightTime += TMath::Sqrt(hcalEnergy); 393 369 394 370 395 fTower->AddCandidate(particle); … … 383 408 Double_t ecalEnergy, hcalEnergy; 384 409 Double_t ecalSigma, hcalSigma; 410 Double_t ecalTime, hcalTime, time; 385 411 386 412 if(!fTower) return; … … 392 418 393 419 ecalEnergy = LogNormal(fTowerECalEnergy, ecalSigma); 420 ecalTime = (fTowerECalWeightTime < 1.0E-09 ) ? 0 : fTowerECalTime/fTowerECalWeightTime; 394 421 395 422 hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerHCalEnergy); … … 399 426 400 427 hcalEnergy = LogNormal(fTowerHCalEnergy, hcalSigma); 428 hcalTime = (fTowerHCalWeightTime < 1.0E-09 ) ? 0 : fTowerHCalTime/fTowerHCalWeightTime; 401 429 402 430 energy = ecalEnergy + hcalEnergy; 431 time = (TMath::Sqrt(ecalEnergy)*ecalTime + TMath::Sqrt(hcalEnergy)*hcalTime)/(TMath::Sqrt(ecalEnergy) + TMath::Sqrt(hcalEnergy)); 403 432 404 433 // eta = fTowerEta; … … 410 439 pt = energy / TMath::CosH(eta); 411 440 412 fTower->Position.SetPtEtaPhiE(1.0, eta, phi, 0.0); 441 // fTower->Position.SetXYZT(-time, 0.0, 0.0, time); 442 fTower->Position.SetPtEtaPhiE(1.0, eta, phi, time); 413 443 fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy); 414 444 fTower->Eem = ecalEnergy; … … 420 450 fTower->Edges[3] = fTowerEdges[3]; 421 451 452 422 453 // fill calorimeter towers and photon candidates 423 454 if(energy > 0.0) -
trunk/modules/Calorimeter.h
r1273 r1345 45 45 Double_t fTowerECalEnergy, fTowerHCalEnergy; 46 46 Double_t fTrackECalEnergy, fTrackHCalEnergy; 47 48 Double_t fTowerECalTime, fTowerHCalTime; 49 Double_t fTrackECalTime, fTrackHCalTime; 50 51 Double_t fTowerECalWeightTime, fTowerHCalWeightTime; 52 Double_t fTrackECalWeightTime, fTrackHCalWeightTime; 53 47 54 Int_t fTowerTrackHits, fTowerPhotonHits; 48 55 -
trunk/modules/FastJetFinder.cc
r1335 r1345 193 193 TLorentzVector momentum; 194 194 Double_t deta, dphi, detaMax, dphiMax; 195 Double_t time, weightTime, avTime; 195 196 Int_t number; 196 197 Double_t rho = 0; … … 259 260 candidate = factory->NewCandidate(); 260 261 262 time=0; 263 weightTime=0; 264 261 265 inputList.clear(); 262 266 inputList = sequence->constituents(*itOutputList); … … 269 273 if(deta > detaMax) detaMax = deta; 270 274 if(dphi > dphiMax) dphiMax = dphi; 271 275 276 time += TMath::Sqrt(constituent->Momentum.E())*(constituent->Position.T()); 277 weightTime += TMath::Sqrt(constituent->Momentum.E()); 278 272 279 candidate->AddCandidate(constituent); 273 280 } 281 282 avTime = time/weightTime; 274 283 275 284 candidate->Momentum = momentum; 285 candidate->Position.SetT(avTime); 276 286 candidate->Area.SetPxPyPzE(area.px(), area.py(), area.pz(), area.E()); 277 287 -
trunk/modules/ModulesLinkDef.h
r1177 r1345 19 19 #include "modules/EnergySmearing.h" 20 20 #include "modules/MomentumSmearing.h" 21 #include "modules/TimeSmearing.h" 21 22 #include "modules/Calorimeter.h" 22 23 #include "modules/Isolation.h" … … 50 51 #pragma link C++ class EnergySmearing+; 51 52 #pragma link C++ class MomentumSmearing+; 53 #pragma link C++ class TimeSmearing+; 52 54 #pragma link C++ class Calorimeter+; 53 55 #pragma link C++ class Isolation+; -
trunk/modules/PileUpMerger.cc
r1323 r1345 16 16 #include "classes/DelphesClasses.h" 17 17 #include "classes/DelphesFactory.h" 18 #include "classes/Delphes Formula.h"18 #include "classes/DelphesTF2.h" 19 19 #include "classes/DelphesPileUpReader.h" 20 20 … … 41 41 42 42 PileUpMerger::PileUpMerger() : 43 fReader(0), fItInputArray(0) 44 { 45 } 43 fFunction(0), fReader(0), fItInputArray(0) 44 { 45 fFunction = new DelphesTF2; 46 } 47 46 48 47 49 //------------------------------------------------------------------------------ … … 49 51 PileUpMerger::~PileUpMerger() 50 52 { 53 delete fFunction; 51 54 } 52 55 … … 60 63 61 64 fMeanPileUp = GetDouble("MeanPileUp", 10); 62 fZVertexSpread = GetDouble("ZVertexSpread", 0.05)*1.0E3; 63 65 66 fZVertexSpread = GetDouble("ZVertexSpread", 0.15); 67 fTVertexSpread = GetDouble("TVertexSpread", 1.5E-09); 68 69 // read vertex smearing formula 70 71 fFunction->Compile(GetString("VertexDistributionFormula", "0.0")); 72 fFunction->SetRange(-fZVertexSpread,-fTVertexSpread,fZVertexSpread,fTVertexSpread); 73 64 74 fileName = GetString("PileUpFile", "MinBias.pileup"); 65 75 fReader = new DelphesPileUpReader(fileName); … … 90 100 Float_t x, y, z, t; 91 101 Float_t px, py, pz, e; 92 Double_t dz, dphi ;102 Double_t dz, dphi, dt; 93 103 Int_t numberOfEvents, event; 94 104 Long64_t allEntries, entry; 95 Candidate *candidate ;105 Candidate *candidate, *vertexcandidate; 96 106 DelphesFactory *factory; 107 108 const Double_t c_light = 2.99792458E8; 97 109 98 110 fItInputArray->Reset(); 111 112 // --- Deal with Primary vertex first ------ 113 114 fFunction->GetRandom2(dz,dt); 115 116 dt *= c_light*1.0E3; //necessary in order to make t in mm/c 117 dz *= 1.0E3; //necessary in order to make z in mm 118 99 119 while((candidate = static_cast<Candidate*>(fItInputArray->Next()))) 100 120 { 121 candidate->Position.SetXYZT(x, y, z+dz, t+dt); 101 122 fParticleOutputArray->Add(candidate); 102 123 } 103 124 104 125 factory = GetFactory(); 105 126 127 vertexcandidate = factory->NewCandidate(); 128 vertexcandidate->Position.SetXYZT(0.0, 0.0, dz, dt); 129 fVertexOutputArray->Add(vertexcandidate); 130 131 132 // --- Then with pile-up vertices ------ 133 106 134 switch(fPileUpDistribution) 107 135 { … … 118 146 119 147 allEntries = fReader->GetEntries(); 120 148 121 149 for(event = 0; event < numberOfEvents; ++event) 122 150 { … … 128 156 129 157 fReader->ReadEntry(entry); 130 131 dz = gRandom->Gaus(0.0, fZVertexSpread); 158 159 // --- Pile-up vertex smearing 160 161 fFunction->GetRandom2(dz,dt); 162 163 dt *= c_light*1.0E3; //necessary in order to make t in mm/c 164 dz *= 1.0E3; //necessary in order to make z in mm 165 132 166 dphi = gRandom->Uniform(-TMath::Pi(), TMath::Pi()); 133 167 134 candidate = factory->NewCandidate(); 135 candidate->Position.SetXYZT(0.0, 0.0, dz, 0.0); 136 fVertexOutputArray->Add(candidate); 168 vertexcandidate = factory->NewCandidate(); 169 vertexcandidate->Position.SetXYZT(0.0, 0.0, dz, dt); 170 vertexcandidate->IsPU = 1; 171 172 fVertexOutputArray->Add(vertexcandidate); 137 173 138 174 while(fReader->ReadParticle(pid, x, y, z, t, px, py, pz, e)) … … 149 185 150 186 candidate->IsPU = 1; 151 187 152 188 candidate->Momentum.SetPxPyPzE(px, py, pz, e); 153 189 candidate->Momentum.RotateZ(dphi); 154 190 155 candidate->Position.SetXYZT(x, y, z + dz,t);191 candidate->Position.SetXYZT(x, y, z+dz, t+dt); 156 192 candidate->Position.RotateZ(dphi); 157 193 158 194 fParticleOutputArray->Add(candidate); 159 195 } -
trunk/modules/PileUpMerger.h
r1323 r1345 19 19 class TObjArray; 20 20 class DelphesPileUpReader; 21 class DelphesTF2; 21 22 22 23 class PileUpMerger: public DelphesModule … … 35 36 Int_t fPileUpDistribution; 36 37 Double_t fMeanPileUp; 38 37 39 Double_t fZVertexSpread; 40 Double_t fTVertexSpread; 38 41 39 42 DelphesPileUpReader *fReader; 40 43 44 DelphesTF2 *fFunction; //! 45 41 46 TIterator *fItInputArray; //! 42 47 -
trunk/modules/TrackPileUpSubtractor.cc
r1314 r1345 53 53 void TrackPileUpSubtractor::Init() 54 54 { 55 // import input array 56 57 fVertexInputArray = ImportArray(GetString("VertexInputArray", "PileUpMerger/vertices")); 58 fItVertexInputArray = fVertexInputArray->MakeIterator(); 59 55 60 fZVertexResolution = GetDouble("ZVertexResolution", 0.005)*1.0E3; 56 61 … … 85 90 if(iterator) delete iterator; 86 91 } 92 93 if(fItVertexInputArray) delete fItVertexInputArray; 87 94 } 88 95 … … 95 102 TIterator *iterator; 96 103 TObjArray *array; 97 Double_t z; 104 Double_t z, zvtx; 105 106 107 // find z position of primary vertex 108 109 fItVertexInputArray->Reset(); 110 while((candidate = static_cast<Candidate*>(fItVertexInputArray->Next()))) 111 { 112 if(candidate->IsPU == 0) 113 { 114 zvtx = candidate->Position.Z(); 115 break; 116 } 117 } 118 98 119 99 120 // loop over all input arrays … … 112 133 // apply pile-up subtraction 113 134 // assume perfect pile-up subtraction for tracks outside fZVertexResolution 114 if(candidate->IsPU && TMath::Abs(z ) > fZVertexResolution) continue;135 if(candidate->IsPU && TMath::Abs(z-zvtx) > fZVertexResolution) continue; 115 136 116 137 array->Add(candidate); -
trunk/modules/TreeWriter.cc
r1325 r1345 95 95 continue; 96 96 } 97 97 98 98 itClassMap = fClassMap.find(branchClass); 99 99 if(itClassMap == fClassMap.end()) … … 115 115 void TreeWriter::Finish() 116 116 { 117 118 117 } 119 118 … … 161 160 GenParticle *entry = 0; 162 161 Double_t pt, signPz, cosTheta, eta, rapidity; 163 162 163 const Double_t c_light = 2.99792458E8; 164 164 165 // loop over all particles 165 166 iterator.Reset(); … … 208 209 entry->Y = position.Y(); 209 210 entry->Z = position.Z(); 210 entry->T = position.T() ;211 entry->T = position.T()*1.0E-3/c_light; 211 212 } 212 213 } … … 219 220 Candidate *candidate = 0; 220 221 Vertex *entry = 0; 222 223 const Double_t c_light = 2.99792458E8; 221 224 222 225 // loop over all vertices … … 231 234 entry->Y = position.Y(); 232 235 entry->Z = position.Z(); 233 entry->T = position.T() ;236 entry->T = position.T()*1.0E-3/c_light; 234 237 } 235 238 } … … 244 247 Track *entry = 0; 245 248 Double_t pt, signz, cosTheta, eta, rapidity; 246 249 const Double_t c_light = 2.99792458E8; 250 247 251 // loop over all tracks 248 252 iterator.Reset(); … … 271 275 entry->YOuter = position.Y(); 272 276 entry->ZOuter = position.Z(); 277 entry->TOuter = position.T()*1.0E-3/c_light; 273 278 274 279 const TLorentzVector &momentum = candidate->Momentum; … … 290 295 entry->Y = initialPosition.Y(); 291 296 entry->Z = initialPosition.Z(); 297 entry->T = initialPosition.T()*1.0E-3/c_light; 292 298 293 299 entry->Particle = particle; … … 303 309 Tower *entry = 0; 304 310 Double_t pt, signPz, cosTheta, eta, rapidity; 305 311 const Double_t c_light = 2.99792458E8; 312 306 313 // loop over all towers 307 314 iterator.Reset(); … … 309 316 { 310 317 const TLorentzVector &momentum = candidate->Momentum; 311 318 const TLorentzVector &position = candidate->Position; 319 312 320 pt = momentum.Pt(); 313 321 cosTheta = TMath::Abs(momentum.CosTheta()); … … 331 339 entry->Edges[2] = candidate->Edges[2]; 332 340 entry->Edges[3] = candidate->Edges[3]; 333 341 342 entry->T = position.T()*1.0E-3/c_light; 343 334 344 FillParticles(candidate, &entry->Particles); 335 345 } … … 344 354 Photon *entry = 0; 345 355 Double_t pt, signPz, cosTheta, eta, rapidity; 346 356 const Double_t c_light = 2.99792458E8; 357 347 358 array->Sort(); 348 359 … … 353 364 TIter it1(candidate->GetCandidates()); 354 365 const TLorentzVector &momentum = candidate->Momentum; 366 const TLorentzVector &position = candidate->Position; 367 355 368 356 369 pt = momentum.Pt(); … … 366 379 entry->PT = pt; 367 380 entry->E = momentum.E(); 368 381 382 entry->T = position.T()*1.0E-3/c_light; 383 369 384 entry->EhadOverEem = candidate->Eem > 0.0 ? candidate->Ehad/candidate->Eem : 999.9; 370 385 … … 381 396 Electron *entry = 0; 382 397 Double_t pt, signPz, cosTheta, eta, rapidity; 383 398 const Double_t c_light = 2.99792458E8; 399 384 400 array->Sort(); 385 401 … … 389 405 { 390 406 const TLorentzVector &momentum = candidate->Momentum; 391 407 const TLorentzVector &position = candidate->Position; 408 392 409 pt = momentum.Pt(); 393 410 cosTheta = TMath::Abs(momentum.CosTheta()); … … 401 418 entry->Phi = momentum.Phi(); 402 419 entry->PT = pt; 403 420 421 entry->T = position.T()*1.0E-3/c_light; 422 404 423 entry->Charge = candidate->Charge; 405 424 … … 418 437 Muon *entry = 0; 419 438 Double_t pt, signPz, cosTheta, eta, rapidity; 420 439 440 const Double_t c_light = 2.99792458E8; 441 421 442 array->Sort(); 422 443 … … 426 447 { 427 448 const TLorentzVector &momentum = candidate->Momentum; 449 const TLorentzVector &position = candidate->Position; 450 428 451 429 452 pt = momentum.Pt(); … … 442 465 entry->PT = pt; 443 466 467 entry->T = position.T()*1.0E-3/c_light; 468 469 // cout<<entry->PT<<","<<entry->T<<endl; 470 444 471 entry->Charge = candidate->Charge; 445 472 … … 457 484 Double_t pt, signPz, cosTheta, eta, rapidity; 458 485 Double_t ecalEnergy, hcalEnergy; 459 486 const Double_t c_light = 2.99792458E8; 487 460 488 array->Sort(); 461 489 … … 465 493 { 466 494 TIter itConstituents(candidate->GetCandidates()); 467 const TLorentzVector &momentum = candidate->Momentum; 468 495 496 const TLorentzVector &momentum = candidate->Momentum; 497 const TLorentzVector &position = candidate->Position; 498 469 499 pt = momentum.Pt(); 470 500 cosTheta = TMath::Abs(momentum.CosTheta()); … … 479 509 entry->PT = pt; 480 510 511 entry->T = position.T()*1.0E-3/c_light; 512 481 513 entry->Mass = momentum.M(); 482 514
Note:
See TracChangeset
for help on using the changeset viewer.