Changes in modules/DualReadoutCalorimeter.cc [5eda6767:de2e39d] in git
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
modules/DualReadoutCalorimeter.cc
r5eda6767 rde2e39d 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.;250 249 while((particle = static_cast<Candidate*>(fItParticleInputArray->Next()))) 251 250 { 252 251 const TLorentzVector &particlePosition = particle->Position; 253 252 ++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();258 253 259 254 pdgCode = TMath::Abs(particle->PID); … … 384 379 fHCalTrackEnergy = 0.0; 385 380 fTrackEnergy = 0.0; 386 381 387 382 fECalTrackSigma = 0.0; 388 383 fHCalTrackSigma = 0.0; 389 384 fTrackSigma = 0.0; 390 385 391 386 fTowerTrackHits = 0; 392 387 fTowerPhotonHits = 0; … … 395 390 fHCalTowerTrackArray->Clear(); 396 391 fTowerTrackArray->Clear(); 397 392 398 393 } 399 394 … … 419 414 } 420 415 421 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 422 447 // in Dual Readout we do not care if tracks are ECAL of HCAL 423 448 if(fECalTrackFractions[number] > 1.0E-9 || fHCalTrackFractions[number] > 1.0E-9) 424 { 449 { 425 450 fTrackEnergy += energy; 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) 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) 427 452 sigma = 0.0; 428 453 if(fHCalTrackFractions[number] > 0) … … 430 455 else 431 456 sigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E()); 432 433 if(sigma/momentum.E() < track->TrackResolution) 434 energyGuess = ecalEnergy + hcalEnergy; 457 458 if(sigma/momentum.E() < track->TrackResolution) 459 energyGuess = ecalEnergy + hcalEnergy; 435 460 else 436 461 energyGuess = momentum.E(); 437 462 438 463 fTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess; 439 464 fTowerTrackArray->Add(track); 440 465 441 466 } 442 467 else … … 471 496 472 497 fTower->AddCandidate(particle); 473 fTower->Position = position;474 498 } 475 499 … … 484 508 485 509 Candidate *track, *tower, *mother; 486 Double_t energy, pt, eta, phi , r;510 Double_t energy, pt, eta, phi; 487 511 Double_t ecalEnergy, hcalEnergy; 488 512 Double_t ecalNeutralEnergy, hcalNeutralEnergy, neutralEnergy; 489 513 490 514 Double_t ecalSigma, hcalSigma, sigma; 491 515 Double_t ecalNeutralSigma, hcalNeutralSigma, neutralSigma; 492 516 493 517 Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor; 494 518 495 519 TLorentzVector momentum; 496 520 TFractionMap::iterator itFractionMap; … … 501 525 502 526 503 // if no hadronic energy, use ECAL resolution 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 504 539 if (fHCalTowerEnergy <= fHCalEnergyMin) 505 540 { … … 509 544 } 510 545 511 // if hadronic fraction > 0, use HCAL resolution 546 // if hadronic fraction > 0, use HCAL resolution 512 547 else 513 548 { … … 519 554 energy = LogNormal(energy, sigma); 520 555 //cout<<energy<<","<<ecalEnergy<<","<<hcalEnergy<<endl; 521 556 522 557 if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0; 523 558 … … 557 592 for(size_t i = 0; i < fTower->ECalEnergyTimePairs.size(); ++i) 558 593 { 559 weight = TMath:: Power((fTower->ECalEnergyTimePairs[i].first),2);594 weight = TMath::Sqrt(fTower->ECalEnergyTimePairs[i].first); 560 595 sumWeightedTime += weight * fTower->ECalEnergyTimePairs[i].second; 561 596 sumWeight += weight; … … 563 598 } 564 599 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); 600 if(sumWeight > 0.0) 601 { 602 fTower->Position.SetPtEtaPhiE(1.0, eta, phi, sumWeightedTime/sumWeight); 603 } 568 604 else 569 r = fTower->Position.Pt(); 570 571 if(sumWeight > 0.0) 572 { 573 fTower->Position.SetPtEtaPhiE(r, eta, phi, sumWeightedTime/sumWeight); 574 } 575 else 576 { 577 fTower->Position.SetPtEtaPhiE(r, eta, phi, 999999.9); 605 { 606 fTower->Position.SetPtEtaPhiE(1.0, eta, phi, 999999.9); 578 607 } 579 608 … … 597 626 598 627 // fill energy flow candidates 599 628 600 629 fTrackSigma = TMath::Sqrt(fTrackSigma); 601 630 neutralEnergy = max( (energy - fTrackEnergy) , 0.0); … … 609 638 if(neutralEnergy > fEnergyMin && neutralSigma > fEnergySignificanceMin) 610 639 { 611 640 612 641 //cout<<"significant neutral excess found:"<<endl; 613 642 // create new photon tower … … 616 645 //cout<<"Creating tower with Pt, Eta, Phi, Energy: "<<pt<<","<<eta<<","<<phi<<","<<neutralEnergy<<endl; 617 646 tower->Momentum.SetPtEtaPhiE(pt, eta, phi, neutralEnergy); 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 } 647 tower->Eem = neutralEnergy; 648 tower->Ehad = 0.0; 649 tower->PID = 22; 650 651 634 652 635 653 fEFlowPhotonOutputArray->Add(tower); … … 640 658 while((track = static_cast<Candidate*>(fItTowerTrackArray->Next()))) 641 659 { 660 //cout<<"looping over tracks"<<endl; 642 661 mother = track; 643 662 track = static_cast<Candidate*>(track->Clone()); … … 647 666 } 648 667 649 668 650 669 // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the DualReadoutCalorimeter and tracking 651 670 else if(fTrackEnergy > 0.0) … … 655 674 weightCalo = (sigma > 0.0) ? 1 / (sigma*sigma) : 0.0; 656 675 657 bestEnergyEstimate = (weightTrack*fTrackEnergy + weightCalo*energy) / (weightTrack + weightCalo); 676 bestEnergyEstimate = (weightTrack*fTrackEnergy + weightCalo*energy) / (weightTrack + weightCalo); 658 677 rescaleFactor = bestEnergyEstimate/fTrackEnergy; 659 678 … … 661 680 fItTowerTrackArray->Reset(); 662 681 while((track = static_cast<Candidate*>(fItTowerTrackArray->Next()))) 663 { 682 { 664 683 mother = track; 665 track = static_cast<Candidate 684 track = static_cast<Candidate*>(track->Clone()); 666 685 track->AddCandidate(mother); 667 track->Momentum.SetPtEtaPhiM(track->Momentum.Pt()*rescaleFactor, track->Momentum.Eta(), track->Momentum.Phi(), track->Momentum.M()); 686 687 track->Momentum *= rescaleFactor; 688 668 689 fEFlowTrackOutputArray->Add(track); 669 690 } 670 691 } 671 692 672 693 673 694 }
Note:
See TracChangeset
for help on using the changeset viewer.