Changes in modules/DualReadoutCalorimeter.cc [de2e39d:5eda6767] in git
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
modules/DualReadoutCalorimeter.cc
rde2e39d r5eda6767 24 24 // This implementation of dual calorimetry relies on several approximations: 25 25 // - If hadronic energy is found in the tower the energy resolution then the full tower enrgy is smeared according to hadronic resolution (pessimistic for (e,n) or (pi+,gamma)) 26 // - While e+ vs pi+ (or gamma vs n) separation is in principle possible for single particles (using C/S, PMT timing, lateral shower profile) it is not obvious it can be done overlapping particles. 26 // - While e+ vs pi+ (or gamma vs n) separation is in principle possible for single particles (using C/S, PMT timing, lateral shower profile) it is not obvious it can be done overlapping particles. 27 27 // Now we assume that regarless of the number of particle hits per tower we can always distinguish e+ vs pi+, which is probably not true in the case (e+,n) vs (pi+,gamma) without longitudinal segmentation. 28 28 … … 63 63 fItParticleInputArray(0), fItTrackInputArray(0) 64 64 { 65 65 66 66 fECalResolutionFormula = new DelphesFormula; 67 67 fHCalResolutionFormula = new DelphesFormula; … … 75 75 fTowerTrackArray = new TObjArray; 76 76 fItTowerTrackArray = fTowerTrackArray->MakeIterator(); 77 77 78 78 } 79 79 … … 82 82 DualReadoutCalorimeter::~DualReadoutCalorimeter() 83 83 { 84 84 85 85 if(fECalResolutionFormula) delete fECalResolutionFormula; 86 86 if(fHCalResolutionFormula) delete fHCalResolutionFormula; … … 94 94 if(fTowerTrackArray) delete fTowerTrackArray; 95 95 if(fItTowerTrackArray) delete fItTowerTrackArray; 96 96 97 97 } 98 98 … … 247 247 fItParticleInputArray->Reset(); 248 248 number = -1; 249 fTowerRmax=0.; 249 250 while((particle = static_cast<Candidate*>(fItParticleInputArray->Next()))) 250 251 { 251 252 const TLorentzVector &particlePosition = particle->Position; 252 253 ++number; 254 255 // compute maximum radius (needed in FinalizeTower to assess whether barrel or endcap tower) 256 if (particlePosition.Perp() > fTowerRmax) 257 fTowerRmax=particlePosition.Perp(); 253 258 254 259 pdgCode = TMath::Abs(particle->PID); … … 379 384 fHCalTrackEnergy = 0.0; 380 385 fTrackEnergy = 0.0; 381 386 382 387 fECalTrackSigma = 0.0; 383 388 fHCalTrackSigma = 0.0; 384 389 fTrackSigma = 0.0; 385 390 386 391 fTowerTrackHits = 0; 387 392 fTowerPhotonHits = 0; … … 390 395 fHCalTowerTrackArray->Clear(); 391 396 fTowerTrackArray->Clear(); 392 397 393 398 } 394 399 … … 414 419 } 415 420 416 417 /* 418 if(fECalTrackFractions[number] > 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9) 419 { 420 fECalTrackEnergy += ecalEnergy; 421 ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E()); 422 if(ecalSigma/momentum.E() < track->TrackResolution) energyGuess = ecalEnergy; 423 else energyGuess = momentum.E(); 424 425 fECalTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess; 426 fECalTowerTrackArray->Add(track); 427 } 428 429 else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] > 1.0E-9) 430 { 431 fHCalTrackEnergy += hcalEnergy; 432 hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E()); 433 if(hcalSigma/momentum.E() < track->TrackResolution) energyGuess = hcalEnergy; 434 else energyGuess = momentum.E(); 435 436 fHCalTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess; 437 fHCalTowerTrackArray->Add(track); 438 } 439 440 // muons 441 else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9) 442 { 443 fEFlowTrackOutputArray->Add(track); 444 } 445 */ 446 421 447 422 // in Dual Readout we do not care if tracks are ECAL of HCAL 448 423 if(fECalTrackFractions[number] > 1.0E-9 || fHCalTrackFractions[number] > 1.0E-9) 449 { 424 { 450 425 fTrackEnergy += energy; 451 // this sigma will be used to determine whether neutral excess is significant. We choose the resolution according to bthe higest deposited fraction (in practice had for charged hadrons and em for electrons) 426 // this sigma will be used to determine whether neutral excess is significant. We choose the resolution according to bthe higest deposited fraction (in practice had for charged hadrons and em for electrons) 452 427 sigma = 0.0; 453 428 if(fHCalTrackFractions[number] > 0) … … 455 430 else 456 431 sigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E()); 457 458 if(sigma/momentum.E() < track->TrackResolution) 459 energyGuess = ecalEnergy + hcalEnergy; 432 433 if(sigma/momentum.E() < track->TrackResolution) 434 energyGuess = ecalEnergy + hcalEnergy; 460 435 else 461 436 energyGuess = momentum.E(); 462 437 463 438 fTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess; 464 439 fTowerTrackArray->Add(track); 465 440 466 441 } 467 442 else … … 496 471 497 472 fTower->AddCandidate(particle); 473 fTower->Position = position; 498 474 } 499 475 … … 508 484 509 485 Candidate *track, *tower, *mother; 510 Double_t energy, pt, eta, phi ;486 Double_t energy, pt, eta, phi, r; 511 487 Double_t ecalEnergy, hcalEnergy; 512 488 Double_t ecalNeutralEnergy, hcalNeutralEnergy, neutralEnergy; 513 489 514 490 Double_t ecalSigma, hcalSigma, sigma; 515 491 Double_t ecalNeutralSigma, hcalNeutralSigma, neutralSigma; 516 492 517 493 Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor; 518 494 519 495 TLorentzVector momentum; 520 496 TFractionMap::iterator itFractionMap; … … 525 501 526 502 527 //if (fHCalTowerEnergy < 30 && fECalTowerEnergy < 30) return; 528 //cout<<"----------- New tower ---------"<<endl; 529 530 531 // here we change behaviour w.r.t to standard calorimeter. Start with worse case scenario. If fHCalTowerEnergy > 0, assume total energy smeared by HCAL resolution. 532 // For example, if overlapping charged pions and photons take hadronic resolution as combined measurement 533 534 // if no hadronic fraction at all, then use ECAL resolution 535 536 //cout<<"fECalTowerEnergy: "<<fECalTowerEnergy<<", fHCalTowerEnergy: "<<fHCalTowerEnergy<<", Eta: "<<fTowerEta<<endl; 537 538 // if no hadronic energy, use ECAL resolution 503 // if no hadronic energy, use ECAL resolution 539 504 if (fHCalTowerEnergy <= fHCalEnergyMin) 540 505 { … … 544 509 } 545 510 546 // if hadronic fraction > 0, use HCAL resolution 511 // if hadronic fraction > 0, use HCAL resolution 547 512 else 548 513 { … … 554 519 energy = LogNormal(energy, sigma); 555 520 //cout<<energy<<","<<ecalEnergy<<","<<hcalEnergy<<endl; 556 521 557 522 if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0; 558 523 … … 592 557 for(size_t i = 0; i < fTower->ECalEnergyTimePairs.size(); ++i) 593 558 { 594 weight = TMath:: Sqrt(fTower->ECalEnergyTimePairs[i].first);559 weight = TMath::Power((fTower->ECalEnergyTimePairs[i].first),2); 595 560 sumWeightedTime += weight * fTower->ECalEnergyTimePairs[i].second; 596 561 sumWeight += weight; … … 598 563 } 599 564 565 // check whether barrel or endcap tower 566 if (fTower->Position.Perp() < fTowerRmax && TMath::Abs(eta) > 0.) 567 r = fTower->Position.Z()/TMath::SinH(eta); 568 else 569 r = fTower->Position.Pt(); 570 600 571 if(sumWeight > 0.0) 601 572 { 602 fTower->Position.SetPtEtaPhiE( 1.0, eta, phi, sumWeightedTime/sumWeight);573 fTower->Position.SetPtEtaPhiE(r, eta, phi, sumWeightedTime/sumWeight); 603 574 } 604 575 else 605 576 { 606 fTower->Position.SetPtEtaPhiE( 1.0, eta, phi, 999999.9);577 fTower->Position.SetPtEtaPhiE(r, eta, phi, 999999.9); 607 578 } 608 579 … … 626 597 627 598 // fill energy flow candidates 628 599 629 600 fTrackSigma = TMath::Sqrt(fTrackSigma); 630 601 neutralEnergy = max( (energy - fTrackEnergy) , 0.0); … … 638 609 if(neutralEnergy > fEnergyMin && neutralSigma > fEnergySignificanceMin) 639 610 { 640 611 641 612 //cout<<"significant neutral excess found:"<<endl; 642 613 // create new photon tower … … 645 616 //cout<<"Creating tower with Pt, Eta, Phi, Energy: "<<pt<<","<<eta<<","<<phi<<","<<neutralEnergy<<endl; 646 617 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, neutralEnergy); 647 tower->Eem = neutralEnergy; 648 tower->Ehad = 0.0; 649 tower->PID = 22; 650 651 618 619 // if no hadronic energy, use ECAL resolution 620 if (fHCalTowerEnergy <= fHCalEnergyMin) 621 { 622 tower->Eem = neutralEnergy; 623 tower->Ehad = 0.0; 624 tower->PID = 22; 625 } 626 627 // if hadronic fraction > 0, use HCAL resolution 628 else 629 { 630 tower->Eem = 0; 631 tower->Ehad = neutralEnergy; 632 tower->PID = 130; 633 } 652 634 653 635 fEFlowPhotonOutputArray->Add(tower); … … 658 640 while((track = static_cast<Candidate*>(fItTowerTrackArray->Next()))) 659 641 { 660 //cout<<"looping over tracks"<<endl;661 642 mother = track; 662 643 track = static_cast<Candidate*>(track->Clone()); … … 666 647 } 667 648 668 649 669 650 // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the DualReadoutCalorimeter and tracking 670 651 else if(fTrackEnergy > 0.0) … … 674 655 weightCalo = (sigma > 0.0) ? 1 / (sigma*sigma) : 0.0; 675 656 676 bestEnergyEstimate = (weightTrack*fTrackEnergy + weightCalo*energy) / (weightTrack + weightCalo); 657 bestEnergyEstimate = (weightTrack*fTrackEnergy + weightCalo*energy) / (weightTrack + weightCalo); 677 658 rescaleFactor = bestEnergyEstimate/fTrackEnergy; 678 659 … … 680 661 fItTowerTrackArray->Reset(); 681 662 while((track = static_cast<Candidate*>(fItTowerTrackArray->Next()))) 682 { 663 { 683 664 mother = track; 684 track = static_cast<Candidate *>(track->Clone());665 track = static_cast<Candidate *>(track->Clone()); 685 666 track->AddCandidate(mother); 686 687 track->Momentum *= rescaleFactor; 688 667 track->Momentum.SetPtEtaPhiM(track->Momentum.Pt()*rescaleFactor, track->Momentum.Eta(), track->Momentum.Phi(), track->Momentum.M()); 689 668 fEFlowTrackOutputArray->Add(track); 690 669 } 691 670 } 692 671 693 672 694 673 }
Note:
See TracChangeset
for help on using the changeset viewer.