- Timestamp:
- Jun 21, 2016, 4:17:30 PM (9 years ago)
- Branches:
- ImprovedOutputFile, Timing, dual_readout, llp, master
- Children:
- efac6f9, f688c89
- Parents:
- d97b2af
- Location:
- modules
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
modules/Calorimeter.cc
rd97b2af r04290b1 58 58 fItParticleInputArray(0), fItTrackInputArray(0) 59 59 { 60 Int_t i; 61 60 62 61 fECalResolutionFormula = new DelphesFormula; 63 62 fHCalResolutionFormula = new DelphesFormula; 64 63 65 for(i = 0; i < 2; ++i) 66 { 67 fECalTowerTrackArray[i] = new TObjArray; 68 fItECalTowerTrackArray[i] = fECalTowerTrackArray[i]->MakeIterator(); 69 70 fHCalTowerTrackArray[i] = new TObjArray; 71 fItHCalTowerTrackArray[i] = fHCalTowerTrackArray[i]->MakeIterator(); 72 } 64 fECalTowerTrackArray = new TObjArray; 65 fItECalTowerTrackArray = fECalTowerTrackArray->MakeIterator(); 66 67 fHCalTowerTrackArray = new TObjArray; 68 fItHCalTowerTrackArray = fHCalTowerTrackArray->MakeIterator(); 69 73 70 } 74 71 … … 77 74 Calorimeter::~Calorimeter() 78 75 { 79 Int_t i; 80 76 81 77 if(fECalResolutionFormula) delete fECalResolutionFormula; 82 78 if(fHCalResolutionFormula) delete fHCalResolutionFormula; 83 79 84 for(i = 0; i < 2; ++i) 85 { 86 if(fECalTowerTrackArray[i]) delete fECalTowerTrackArray[i]; 87 if(fItECalTowerTrackArray[i]) delete fItECalTowerTrackArray[i]; 88 89 if(fHCalTowerTrackArray[i]) delete fHCalTowerTrackArray[i]; 90 if(fItHCalTowerTrackArray[i]) delete fItHCalTowerTrackArray[i]; 91 } 80 if(fECalTowerTrackArray) delete fECalTowerTrackArray; 81 if(fItECalTowerTrackArray) delete fItECalTowerTrackArray; 82 83 if(fHCalTowerTrackArray) delete fHCalTowerTrackArray; 84 if(fItHCalTowerTrackArray) delete fItHCalTowerTrackArray; 85 92 86 } 93 87 … … 218 212 Double_t ecalFraction, hcalFraction; 219 213 Double_t ecalEnergy, hcalEnergy; 220 Double_t ecalSigma, hcalSigma;221 214 Int_t pdgCode; 222 215 … … 368 361 fHCalTowerEnergy = 0.0; 369 362 370 fECalTrackEnergy [0]= 0.0;371 f ECalTrackEnergy[1]= 0.0;372 373 f HCalTrackEnergy[0]= 0.0;374 fHCalTrack Energy[1]= 0.0;375 363 fECalTrackEnergy = 0.0; 364 fHCalTrackEnergy = 0.0; 365 366 fECalTrackSigma = 0.0; 367 fHCalTrackSigma = 0.0; 368 376 369 fTowerTrackHits = 0; 377 370 fTowerPhotonHits = 0; 378 371 379 fECalTowerTrackArray[0]->Clear(); 380 fECalTowerTrackArray[1]->Clear(); 381 382 fHCalTowerTrackArray[0]->Clear(); 383 fHCalTowerTrackArray[1]->Clear(); 372 fECalTowerTrackArray->Clear(); 373 fHCalTowerTrackArray->Clear(); 374 384 375 } 385 376 … … 406 397 if(fECalTrackFractions[number] > 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9) 407 398 { 408 ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E()); 409 if(ecalSigma/momentum.E() < track->TrackResolution) 410 { 411 fECalTrackEnergy[0] += ecalEnergy; 412 fECalTowerTrackArray[0]->Add(track); 413 } 414 else 415 { 416 fECalTrackEnergy[1] += ecalEnergy; 417 fECalTowerTrackArray[1]->Add(track); 418 } 399 fECalTrackEnergy += ecalEnergy; 400 fECalTrackSigma += (track->TrackResolution)*momentum.E()*(track->TrackResolution)*momentum.E(); 401 fECalTowerTrackArray->Add(track); 419 402 } 403 420 404 else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] > 1.0E-9) 421 405 { 422 hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E()); 423 if(hcalSigma/momentum.E() < track->TrackResolution) 424 { 425 fHCalTrackEnergy[0] += hcalEnergy; 426 fHCalTowerTrackArray[0]->Add(track); 427 } 428 else 429 { 430 fHCalTrackEnergy[1] += hcalEnergy; 431 fHCalTowerTrackArray[1]->Add(track); 432 } 406 fHCalTrackEnergy += hcalEnergy; 407 fHCalTrackSigma += (track->TrackResolution)*momentum.E()*(track->TrackResolution)*momentum.E(); 408 fHCalTowerTrackArray->Add(track); 433 409 } 410 434 411 else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9) 435 412 { … … 476 453 Double_t energy, pt, eta, phi; 477 454 Double_t ecalEnergy, hcalEnergy; 455 Double_t ecalNeutralEnergy, hcalNeutralEnergy; 456 478 457 Double_t ecalSigma, hcalSigma; 479 458 Double_t ecalNeutralSigma, hcalNeutralSigma; 459 460 Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor; 461 480 462 TLorentzVector momentum; 481 463 TFractionMap::iterator itFractionMap; … … 484 466 485 467 if(!fTower) return; 486 487 468 488 469 ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fECalTowerEnergy); … … 554 535 fTowerOutputArray->Add(fTower); 555 536 } 556 537 557 538 // fill energy flow candidates 558 559 ecalEnergy -= fECalTrackEnergy[1]; 560 hcalEnergy -= fHCalTrackEnergy[1]; 561 562 fItECalTowerTrackArray[0]->Reset(); 563 while((track = static_cast<Candidate*>(fItECalTowerTrackArray[0]->Next()))) 564 { 565 mother = track; 566 track = static_cast<Candidate*>(track->Clone()); 567 track->AddCandidate(mother); 568 569 track->Momentum *= ecalEnergy/fECalTrackEnergy[0]; 570 571 fEFlowTrackOutputArray->Add(track); 572 } 573 574 fItECalTowerTrackArray[1]->Reset(); 575 while((track = static_cast<Candidate*>(fItECalTowerTrackArray[1]->Next()))) 576 { 577 mother = track; 578 track = static_cast<Candidate*>(track->Clone()); 579 track->AddCandidate(mother); 580 581 fEFlowTrackOutputArray->Add(track); 582 } 583 584 fItHCalTowerTrackArray[0]->Reset(); 585 while((track = static_cast<Candidate*>(fItHCalTowerTrackArray[0]->Next()))) 586 { 587 mother = track; 588 track = static_cast<Candidate*>(track->Clone()); 589 track->AddCandidate(mother); 590 591 track->Momentum *= hcalEnergy/fHCalTrackEnergy[0]; 592 593 fEFlowTrackOutputArray->Add(track); 594 } 595 596 fItHCalTowerTrackArray[1]->Reset(); 597 while((track = static_cast<Candidate*>(fItHCalTowerTrackArray[1]->Next()))) 598 { 599 mother = track; 600 track = static_cast<Candidate*>(track->Clone()); 601 track->AddCandidate(mother); 602 603 fEFlowTrackOutputArray->Add(track); 604 } 605 606 if(fECalTowerTrackArray[0]->GetEntriesFast() > 0) ecalEnergy = 0.0; 607 if(fHCalTowerTrackArray[0]->GetEntriesFast() > 0) hcalEnergy = 0.0; 608 609 ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy); 610 hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy); 611 612 if(ecalEnergy < fECalEnergyMin || ecalEnergy < fECalEnergySignificanceMin*ecalSigma) ecalEnergy = 0.0; 613 if(hcalEnergy < fHCalEnergyMin || hcalEnergy < fHCalEnergySignificanceMin*hcalSigma) hcalEnergy = 0.0; 614 615 energy = ecalEnergy + hcalEnergy; 616 617 if(ecalEnergy > 0.0) 539 fECalTrackSigma = TMath::Sqrt(fECalTrackSigma); 540 fHCalTrackSigma = TMath::Sqrt(fHCalTrackSigma); 541 542 //compute neutral excesses 543 ecalNeutralEnergy = max( (ecalEnergy - fECalTrackEnergy) , 0.0); 544 hcalNeutralEnergy = max( (hcalEnergy - fHCalTrackEnergy) , 0.0); 545 546 ecalNeutralSigma = ecalNeutralEnergy / TMath::Sqrt(fECalTrackSigma*fECalTrackSigma + ecalSigma*ecalSigma); 547 hcalNeutralSigma = hcalNeutralEnergy / TMath::Sqrt(fHCalTrackSigma*fHCalTrackSigma + hcalSigma*hcalSigma); 548 549 // if ecal neutral excess is significant, simply create neutral EflowPhoton tower and clone each track into eflowtrack 550 if(ecalNeutralEnergy > fECalEnergyMin && ecalNeutralSigma > fECalEnergySignificanceMin) 618 551 { 619 552 // create new photon tower 620 553 tower = static_cast<Candidate*>(fTower->Clone()); 621 622 pt = ecalEnergy / TMath::CosH(eta); 623 624 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalEnergy); 625 tower->Eem = ecalEnergy; 554 pt = ecalNeutralEnergy / TMath::CosH(eta); 555 556 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalNeutralEnergy); 557 tower->Eem = ecalNeutralEnergy; 626 558 tower->Ehad = 0.0; 627 tower->PID = 22; 628 559 629 560 fEFlowPhotonOutputArray->Add(tower); 630 } 631 if(hcalEnergy > 0.0) 632 { 633 // create new neutral hadron tower 561 562 //clone tracks 563 fItECalTowerTrackArray->Reset(); 564 while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next()))) 565 { 566 mother = track; 567 track = static_cast<Candidate*>(track->Clone()); 568 track->AddCandidate(mother); 569 570 fEFlowTrackOutputArray->Add(track); 571 } 572 573 } 574 575 // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the calorimeter and tracking 576 else if(fECalTrackEnergy > 0.0) 577 { 578 weightTrack = (fECalTrackSigma > 0.0) ? 1 / (fECalTrackSigma*fECalTrackSigma) : 0.0; 579 weightCalo = (ecalSigma > 0.0) ? 1 / (ecalSigma*ecalSigma) : 0.0; 580 581 bestEnergyEstimate = (weightTrack*fECalTrackEnergy + weightCalo*ecalEnergy) / (weightTrack + weightCalo); 582 rescaleFactor = bestEnergyEstimate/fECalTrackEnergy; 583 584 //rescale tracks 585 fItECalTowerTrackArray->Reset(); 586 while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next()))) 587 { 588 mother = track; 589 track = static_cast<Candidate*>(track->Clone()); 590 track->AddCandidate(mother); 591 592 track->Momentum *= rescaleFactor; 593 594 fEFlowTrackOutputArray->Add(track); 595 } 596 } 597 598 599 // if hcal neutral excess is significant, simply create neutral EflowNeutralHadron tower and clone each track into eflowtrack 600 if(hcalNeutralEnergy > fHCalEnergyMin && hcalNeutralSigma > fHCalEnergySignificanceMin) 601 { 602 // create new photon tower 634 603 tower = static_cast<Candidate*>(fTower->Clone()); 635 636 pt = hcalEnergy / TMath::CosH(eta);637 638 tower-> Momentum.SetPtEtaPhiE(pt, eta, phi, hcalEnergy);604 pt = hcalNeutralEnergy / TMath::CosH(eta); 605 606 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalNeutralEnergy); 607 tower->Ehad = hcalNeutralEnergy; 639 608 tower->Eem = 0.0; 640 tower->Ehad = hcalEnergy; 641 609 642 610 fEFlowNeutralHadronOutputArray->Add(tower); 643 } 611 612 //clone tracks 613 fItHCalTowerTrackArray->Reset(); 614 while((track = static_cast<Candidate*>(fItHCalTowerTrackArray->Next()))) 615 { 616 mother = track; 617 track = static_cast<Candidate*>(track->Clone()); 618 track->AddCandidate(mother); 619 620 fEFlowTrackOutputArray->Add(track); 621 } 622 623 } 624 625 // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the calorimeter and tracking 626 else if(fHCalTrackEnergy > 0.0) 627 { 628 weightTrack = (fHCalTrackSigma > 0.0) ? 1 / (fHCalTrackSigma*fHCalTrackSigma) : 0.0; 629 weightCalo = (hcalSigma > 0.0) ? 1 / (hcalSigma*hcalSigma) : 0.0; 630 631 bestEnergyEstimate = (weightTrack*fHCalTrackEnergy + weightCalo*hcalEnergy) / (weightTrack + weightCalo); 632 rescaleFactor = bestEnergyEstimate / fHCalTrackEnergy; 633 634 //rescale tracks 635 fItECalTowerTrackArray->Reset(); 636 while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next()))) 637 { 638 mother = track; 639 track = static_cast<Candidate*>(track->Clone()); 640 track->AddCandidate(mother); 641 642 track->Momentum *= rescaleFactor; 643 644 fEFlowTrackOutputArray->Add(track); 645 } 646 } 647 648 644 649 } 645 650 -
modules/Calorimeter.h
rd97b2af r04290b1 58 58 Double_t fTowerEta, fTowerPhi, fTowerEdges[4]; 59 59 Double_t fECalTowerEnergy, fHCalTowerEnergy; 60 Double_t fECalTrackEnergy [2], fHCalTrackEnergy[2];60 Double_t fECalTrackEnergy, fHCalTrackEnergy; 61 61 62 62 Double_t fTimingEnergyMin; … … 70 70 Double_t fECalEnergySignificanceMin; 71 71 Double_t fHCalEnergySignificanceMin; 72 73 Double_t fECalTrackSigma; 74 Double_t fHCalTrackSigma; 72 75 73 76 Bool_t fSmearTowerCenter; … … 103 106 TObjArray *fEFlowNeutralHadronOutputArray; //! 104 107 105 TObjArray *fECalTowerTrackArray [2]; //!106 TIterator *fItECalTowerTrackArray [2]; //!108 TObjArray *fECalTowerTrackArray; //! 109 TIterator *fItECalTowerTrackArray; //! 107 110 108 TObjArray *fHCalTowerTrackArray [2]; //!109 TIterator *fItHCalTowerTrackArray [2]; //!111 TObjArray *fHCalTowerTrackArray; //! 112 TIterator *fItHCalTowerTrackArray; //! 110 113 111 114 void FinalizeTower();
Note:
See TracChangeset
for help on using the changeset viewer.