Fork me on GitHub

Changeset 1273 in svn for trunk


Ignore:
Timestamp:
Sep 3, 2013, 5:54:56 PM (11 years ago)
Author:
Pavel Demin
Message:

new energy flow

Location:
trunk/modules
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/modules/Calorimeter.cc

    r1246 r1273  
    4343  fECalResolutionFormula(0), fHCalResolutionFormula(0),
    4444  fItParticleInputArray(0), fItTrackInputArray(0),
    45   fTowerECalArray(0), fItTowerECalArray(0),
    46   fTowerHCalArray(0), fItTowerHCalArray(0),
    47   fTowerTrackArray(0), fItTowerTrackArray(0),
    48   fTowerECalTrackArray(0), fItTowerECalTrackArray(0),
    49   fTowerHCalTrackArray(0), fItTowerHCalTrackArray(0)
     45  fTowerTrackArray(0), fItTowerTrackArray(0)
    5046{
    5147  fECalResolutionFormula = new DelphesFormula;
    5248  fHCalResolutionFormula = new DelphesFormula;
    5349
    54   fTowerECalArray = new TObjArray;
    55   fItTowerECalArray = fTowerECalArray->MakeIterator();
    56   fTowerHCalArray = new TObjArray;
    57   fItTowerHCalArray = fTowerHCalArray->MakeIterator();
    58 
    5950  fTowerTrackArray = new TObjArray;
    6051  fItTowerTrackArray = fTowerTrackArray->MakeIterator();
    61   fTowerECalTrackArray = new TObjArray;
    62   fItTowerECalTrackArray = fTowerECalTrackArray->MakeIterator();
    63   fTowerHCalTrackArray = new TObjArray;
    64   fItTowerHCalTrackArray = fTowerHCalTrackArray->MakeIterator();
    6552}
    6653
     
    7259  if(fHCalResolutionFormula) delete fHCalResolutionFormula;
    7360
    74   if(fTowerECalArray) delete fTowerECalArray;
    75   if(fItTowerECalArray) delete fItTowerECalArray;
    76   if(fTowerHCalArray) delete fTowerHCalArray;
    77   if(fItTowerHCalArray) delete fItTowerHCalArray;
    78 
    7961  if(fTowerTrackArray) delete fTowerTrackArray;
    8062  if(fItTowerTrackArray) delete fItTowerTrackArray;
    81   if(fTowerECalTrackArray) delete fTowerECalTrackArray;
    82   if(fItTowerECalTrackArray) delete fItTowerECalTrackArray;
    83   if(fTowerHCalTrackArray) delete fTowerHCalTrackArray;
    84   if(fItTowerHCalTrackArray) delete fItTowerHCalTrackArray;}
     63}
    8564
    8665//------------------------------------------------------------------------------
     
    179158void Calorimeter::Finish()
    180159{
    181   vector< vector< Double_t>* >::iterator itPhiBin;
     160  vector< vector< Double_t >* >::iterator itPhiBin;
    182161  if(fItParticleInputArray) delete fItParticleInputArray;
    183162  if(fItTrackInputArray) delete fItTrackInputArray;
     
    211190  DelphesFactory *factory = GetFactory();
    212191  fTowerHits.clear();
    213   fECalFractions.clear();
    214   fHCalFractions.clear();
     192  fTowerECalFractions.clear();
     193  fTowerHCalFractions.clear();
     194  fTrackECalFractions.clear();
     195  fTrackHCalFractions.clear();
    215196
    216197  // loop over all particles
     
    233214    hcalFraction = itFractionMap->second.second;
    234215
    235     fECalFractions.push_back(ecalFraction);
    236     fHCalFractions.push_back(hcalFraction);
     216    fTowerECalFractions.push_back(ecalFraction);
     217    fTowerHCalFractions.push_back(hcalFraction);
    237218
    238219    if(ecalFraction < 1.0E-9 && hcalFraction < 1.0E-9) continue;
     
    252233
    253234    flags = 0;
    254     flags |= (ecalFraction >= 1.0E-9) << 1;
    255     flags |= (hcalFraction >= 1.0E-9) << 2;
    256     flags |= (pdgCode == 11 || pdgCode == 22) << 3;
     235    flags |= (pdgCode == 11 || pdgCode == 22) << 1;
    257236
    258237    // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for particle number}
     
    281260    hcalFraction = itFractionMap->second.second;
    282261
     262    fTrackECalFractions.push_back(ecalFraction);
     263    fTrackHCalFractions.push_back(hcalFraction);
     264
     265    if(ecalFraction < 1.0E-9 && hcalFraction < 1.0E-9) continue;
     266
    283267    // find eta bin [1, fEtaBins.size - 1]
    284268    itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), trackPosition.Eta());
     
    295279
    296280    flags = 1;
    297     flags |= (ecalFraction >= 1.0E-9) << 1;
    298     flags |= (hcalFraction >= 1.0E-9) << 2;
    299281
    300282    // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for track number}
     
    347329      fTowerHCalEnergy = 0.0;
    348330
     331      fTrackECalEnergy = 0.0;
     332      fTrackHCalEnergy = 0.0;
     333
     334      fTowerTrackHits = 0;
    349335      fTowerPhotonHits = 0;
    350336
    351       fTowerAllHits = 0;
    352       fTowerECalHits = 0;
    353       fTowerHCalHits = 0;
    354 
    355       fTowerTrackAllHits = 0;
    356       fTowerECalTrackHits = 0;
    357       fTowerHCalTrackHits = 0;
    358 
    359       fTowerECalArray->Clear();
    360       fTowerHCalArray->Clear();
    361 
    362337      fTowerTrackArray->Clear();
    363       fTowerECalTrackArray->Clear();
    364       fTowerHCalTrackArray->Clear();
    365338    }
    366339
     
    368341    if(flags & 1)
    369342    {
     343      ++fTowerTrackHits;
     344
    370345      track = static_cast<Candidate*>(fTrackInputArray->At(number));
    371 
    372       ++fTowerTrackAllHits;
     346      momentum = track->Momentum;
     347
     348      ecalEnergy = momentum.E() * fTrackECalFractions[number];
     349      hcalEnergy = momentum.E() * fTrackHCalFractions[number];
     350
     351      fTrackECalEnergy += ecalEnergy;
     352      fTrackHCalEnergy += hcalEnergy;
     353
    373354      fTowerTrackArray->Add(track);
    374355
    375       // check for track ECAL hits
    376       if(flags & 2)
    377       {
    378         ++fTowerECalTrackHits;
    379         fTowerECalTrackArray->Add(track);
    380       }
    381 
    382       // check for track HCAL hits
    383       if(flags & 4)
    384       {
    385         ++fTowerHCalTrackHits;
    386         fTowerHCalTrackArray->Add(track);
    387       }
    388356      continue;
    389357    }
    390358
    391     ++fTowerAllHits;
    392 
    393     // check for ECAL hits
    394     if(flags & 2)
    395     {
    396       ++fTowerECalHits;
    397       fTowerECalArray->Add(particle);
    398     }
    399 
    400     // check for HCAL hits
    401     if(flags & 4)
    402     {
    403       ++fTowerHCalHits;
    404       fTowerHCalArray->Add(particle);
    405     }
    406 
    407359    // check for photon and electron hits in current tower
    408     if(flags & 8) ++fTowerPhotonHits;
     360    if(flags & 2) ++fTowerPhotonHits;
    409361
    410362    particle = static_cast<Candidate*>(fParticleInputArray->At(number));
     
    412364
    413365    // fill current tower
    414     ecalEnergy = momentum.E() * fECalFractions[number];
    415     hcalEnergy = momentum.E() * fHCalFractions[number];
     366    ecalEnergy = momentum.E() * fTowerECalFractions[number];
     367    hcalEnergy = momentum.E() * fTowerHCalFractions[number];
    416368
    417369    fTowerECalEnergy += ecalEnergy;
     
    429381void Calorimeter::FinalizeTower()
    430382{
    431   Candidate *particle, *track, *tower;
     383  Candidate *track, *tower;
    432384  Double_t energy, pt, eta, phi;
    433385  Double_t ecalEnergy, hcalEnergy;
    434   TIterator *itTowerTrackArray;
    435386
    436387  if(!fTower) return;
     
    469420  if(energy > 0.0)
    470421  {
    471     if(fTowerPhotonHits > 0 && fTowerTrackAllHits == 0)
     422    if(fTowerPhotonHits > 0 && fTowerTrackHits == 0)
    472423    {
    473424      fPhotonOutputArray->Add(fTower);
     
    478429
    479430  // fill energy flow candidates
    480   if(fTowerTrackAllHits == fTowerAllHits)
    481   {
    482     fItTowerTrackArray->Reset();
    483     while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
    484     {
    485       fEFlowTrackOutputArray->Add(track);
    486     }
    487   }
    488   else if(fTowerTrackAllHits > 0 &&
    489           fTowerECalHits + fTowerHCalHits == fTowerAllHits)
    490   {
    491     if(fTowerECalHits == fTowerECalTrackHits &&
    492        fTowerHCalHits == fTowerHCalTrackHits)
    493     {
    494       itTowerTrackArray = fItTowerTrackArray;
    495     }
    496     else if(fTowerECalHits == fTowerECalTrackHits)
    497     {
    498       itTowerTrackArray = fItTowerECalTrackArray;
    499 
    500       if(hcalEnergy > 0.0)
    501       {
    502         DelphesFactory *factory = GetFactory();
    503 
    504         // create new tower
    505         tower = factory->NewCandidate();
    506 
    507         fItTowerHCalArray->Reset();
    508         while((particle = static_cast<Candidate*>(fItTowerHCalArray->Next())))
    509         {
    510           tower->AddCandidate(particle);
    511         }
    512 
    513         pt = hcalEnergy / TMath::CosH(eta);
    514 
    515         tower->Position.SetPtEtaPhiE(1.0, eta, phi, 0.0);
    516         tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalEnergy);
    517         tower->Eem = 0.0;
    518         tower->Ehad = hcalEnergy;
    519 
    520         tower->Edges[0] = fTowerEdges[0];
    521         tower->Edges[1] = fTowerEdges[1];
    522         tower->Edges[2] = fTowerEdges[2];
    523         tower->Edges[3] = fTowerEdges[3];
    524 
    525         fEFlowTowerOutputArray->Add(tower);
    526       }
    527     }
    528     else if(fTowerHCalHits == fTowerHCalTrackHits)
    529     {
    530       itTowerTrackArray = fItTowerHCalTrackArray;
    531 
    532       if(ecalEnergy > 0.0)
    533       {
    534         DelphesFactory *factory = GetFactory();
    535 
    536         // create new tower
    537         tower = factory->NewCandidate();
    538 
    539         fItTowerECalArray->Reset();
    540         while((particle = static_cast<Candidate*>(fItTowerECalArray->Next())))
    541         {
    542           tower->AddCandidate(particle);
    543         }
    544 
    545         pt = ecalEnergy / TMath::CosH(eta);
    546 
    547         tower->Position.SetPtEtaPhiE(1.0, eta, phi, 0.0);
    548         tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalEnergy);
    549         tower->Eem = ecalEnergy;
    550         tower->Ehad = 0.0;
    551 
    552         tower->Edges[0] = fTowerEdges[0];
    553         tower->Edges[1] = fTowerEdges[1];
    554         tower->Edges[2] = fTowerEdges[2];
    555         tower->Edges[3] = fTowerEdges[3];
    556 
    557         fEFlowTowerOutputArray->Add(tower);
    558       }
    559     }
    560     else
    561     {
    562       itTowerTrackArray = 0;
    563       fEFlowTowerOutputArray->Add(fTower);
    564     }
    565 
    566     if(itTowerTrackArray)
    567     {
    568       itTowerTrackArray->Reset();
    569       while((track = static_cast<Candidate*>(itTowerTrackArray->Next())))
    570       {
    571         fEFlowTrackOutputArray->Add(track);
    572       }
    573     }
    574   }
    575   else if(energy > 0.0)
    576   {
    577     fEFlowTowerOutputArray->Add(fTower);
     431
     432  // save all the tracks as energy flow tracks
     433  fItTowerTrackArray->Reset();
     434  while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
     435  {
     436    fEFlowTrackOutputArray->Add(track);
     437  }
     438
     439  ecalEnergy -= fTrackECalEnergy;
     440  if(ecalEnergy < 0.0) ecalEnergy = 0.0;
     441
     442  hcalEnergy -= fTrackHCalEnergy;
     443  if(hcalEnergy < 0.0) hcalEnergy = 0.0;
     444
     445  energy = ecalEnergy + hcalEnergy;
     446
     447  // save ECAL and/or HCAL energy excess as an energy flow tower
     448  if(energy > 0.0)
     449  {
     450    // create new tower
     451    tower = static_cast<Candidate*>(fTower->Clone());
     452
     453    pt = energy / TMath::CosH(eta);
     454
     455    tower->Position.SetPtEtaPhiE(1.0, eta, phi, 0.0);
     456    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
     457    tower->Eem = ecalEnergy;
     458    tower->Ehad = hcalEnergy;
     459
     460    tower->Edges[0] = fTowerEdges[0];
     461    tower->Edges[1] = fTowerEdges[1];
     462    tower->Edges[2] = fTowerEdges[2];
     463    tower->Edges[3] = fTowerEdges[3];
     464
     465    fEFlowTowerOutputArray->Add(tower);
    578466  }
    579467}
  • trunk/modules/Calorimeter.h

    r1235 r1273  
    4444  Double_t fTowerEta, fTowerPhi, fTowerEdges[4];
    4545  Double_t fTowerECalEnergy, fTowerHCalEnergy;
    46   Double_t fTowerECalNeutralEnergy, fTowerHCalNeutralEnergy;
    47   Int_t fTowerPhotonHits, fTowerECalHits, fTowerHCalHits, fTowerAllHits;
    48   Int_t fTowerECalTrackHits, fTowerHCalTrackHits, fTowerTrackAllHits;
     46  Double_t fTrackECalEnergy, fTrackHCalEnergy;
     47  Int_t fTowerTrackHits, fTowerPhotonHits;
    4948
    5049  TFractionMap fFractionMap; //!
     
    5655  std::vector < Long64_t > fTowerHits;
    5756
    58   std::vector < Double_t > fECalFractions;
    59   std::vector < Double_t > fHCalFractions;
     57  std::vector < Double_t > fTowerECalFractions;
     58  std::vector < Double_t > fTowerHCalFractions;
     59
     60  std::vector < Double_t > fTrackECalFractions;
     61  std::vector < Double_t > fTrackHCalFractions;
    6062
    6163  DelphesFormula *fECalResolutionFormula; //!
     
    7476  TObjArray *fEFlowTowerOutputArray; //!
    7577
    76   TObjArray *fTowerECalArray; //!
    77   TIterator *fItTowerECalArray; //!
    78 
    79   TObjArray *fTowerHCalArray; //!
    80   TIterator *fItTowerHCalArray; //!
    81 
    8278  TObjArray *fTowerTrackArray; //!
    8379  TIterator *fItTowerTrackArray; //!
    84 
    85   TObjArray *fTowerECalTrackArray; //!
    86   TIterator *fItTowerECalTrackArray; //!
    87 
    88   TObjArray *fTowerHCalTrackArray; //!
    89   TIterator *fItTowerHCalTrackArray; //!
    9080
    9181  void FinalizeTower();
Note: See TracChangeset for help on using the changeset viewer.