Fork me on GitHub

Changeset 73e0386 in git


Ignore:
Timestamp:
Jul 22, 2013, 5:41:43 PM (11 years ago)
Author:
pavel <pavel@…>
Branches:
ImprovedOutputFile, Timing, dual_readout, llp, master
Children:
9493a0f
Parents:
d60d554
Message:

new energy flow

Location:
modules
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • modules/Calorimeter.cc

    rd60d554 r73e0386  
    4343  fECalResolutionFormula(0), fHCalResolutionFormula(0),
    4444  fItParticleInputArray(0), fItTrackInputArray(0),
    45   fTowerTrackArray(0), fItTowerTrackArray(0),
    46   fTowerPhotonArray(0), fItTowerPhotonArray(0)
     45  fTowerECalArray(0), fItTowerECalArray(0),
     46  fTowerHCalArray(0), fItTowerHCalArray(0),
     47  fTowerECalTrackArray(0), fItTowerECalTrackArray(0),
     48  fTowerHCalTrackArray(0), fItTowerHCalTrackArray(0)
    4749{
    4850  fECalResolutionFormula = new DelphesFormula;
    4951  fHCalResolutionFormula = new DelphesFormula;
    50   fTowerTrackArray = new TObjArray;
    51   fItTowerTrackArray = fTowerTrackArray->MakeIterator();
    52   fTowerPhotonArray = new TObjArray;
    53   fItTowerPhotonArray = fTowerPhotonArray->MakeIterator();
     52
     53  fTowerECalArray = new TObjArray;
     54  fItTowerECalArray = fTowerECalArray->MakeIterator();
     55  fTowerHCalArray = new TObjArray;
     56  fItTowerHCalArray = fTowerHCalArray->MakeIterator();
     57
     58  fTowerECalTrackArray = new TObjArray;
     59  fItTowerECalTrackArray = fTowerECalTrackArray->MakeIterator();
     60  fTowerHCalTrackArray = new TObjArray;
     61  fItTowerHCalTrackArray = fTowerHCalTrackArray->MakeIterator();
    5462}
    5563
     
    6068  if(fECalResolutionFormula) delete fECalResolutionFormula;
    6169  if(fHCalResolutionFormula) delete fHCalResolutionFormula;
    62   if(fTowerTrackArray) delete fTowerTrackArray;
    63   if(fItTowerTrackArray) delete fItTowerTrackArray;
    64   if(fTowerPhotonArray) delete fTowerPhotonArray;
    65   if(fItTowerPhotonArray) delete fItTowerPhotonArray;
    66 }
     70
     71  if(fTowerECalArray) delete fTowerECalArray;
     72  if(fItTowerECalArray) delete fItTowerECalArray;
     73  if(fTowerHCalArray) delete fTowerHCalArray;
     74  if(fItTowerHCalArray) delete fItTowerHCalArray;
     75
     76  if(fTowerECalTrackArray) delete fTowerECalTrackArray;
     77  if(fItTowerECalTrackArray) delete fItTowerECalTrackArray;
     78  if(fTowerHCalTrackArray) delete fTowerHCalTrackArray;
     79  if(fItTowerHCalTrackArray) delete fItTowerHCalTrackArray;}
    6780
    6881//------------------------------------------------------------------------------
     
    233246    phiBin = distance(phiBins->begin(), itPhiBin);
    234247
    235     flags = (particle->Charge == 0);
    236     flags |= (pdgCode == 22) << 1;
    237     flags |= (pdgCode == 11) << 2;
     248    flags = 0;
     249    flags |= (ecalFraction >= 1.0E-9) << 1;
     250    flags |= (hcalFraction >= 1.0E-9) << 2;
     251    flags |= (pdgCode == 11 || pdgCode == 22) << 3;
    238252
    239253    // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for particle number}
     
    251265    ++number;
    252266
     267    pdgCode = TMath::Abs(track->PID);
     268
     269    itFractionMap = fFractionMap.find(pdgCode);
     270    if(itFractionMap == fFractionMap.end())
     271    {
     272      itFractionMap = fFractionMap.find(0);
     273    }
     274
     275    ecalFraction = itFractionMap->second.first;
     276    hcalFraction = itFractionMap->second.second;
     277
    253278    // find eta bin [1, fEtaBins.size - 1]
    254279    itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), trackPosition.Eta());
     
    264289    phiBin = distance(phiBins->begin(), itPhiBin);
    265290
     291    flags = 1;
     292    flags |= (ecalFraction >= 1.0E-9) << 1;
     293    flags |= (hcalFraction >= 1.0E-9) << 2;
     294
    266295    // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for track number}
    267     towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(1) << 27) | Long64_t(number);
     296    towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
    268297
    269298    fTowerHits.push_back(towerHit);
     
    313342      fTowerHCalEnergy = 0.0;
    314343
    315       fTowerECalNeutralEnergy = 0.0;
    316       fTowerHCalNeutralEnergy = 0.0;
    317 
    318       fTowerNeutralHits = 0;
    319344      fTowerPhotonHits = 0;
    320       fTowerElectronHits = 0;
    321       fTowerTrackHits = 0;
     345
    322346      fTowerAllHits = 0;
    323 
    324       fTowerTrackArray->Clear();
    325       fTowerPhotonArray->Clear();
     347      fTowerECalHits = 0;
     348      fTowerHCalHits = 0;
     349
     350      fTowerTrackAllHits = 0;
     351      fTowerECalTrackHits = 0;
     352      fTowerHCalTrackHits = 0;
     353
     354      fTowerECalTrackArray->Clear();
     355      fTowerHCalTrackArray->Clear();
     356      fTowerECalArray->Clear();
     357      fTowerHCalArray->Clear();
    326358    }
    327359
    328360    // check for track hits
    329     if(flags & 8)
    330     {
    331       ++fTowerTrackHits;
     361    if(flags & 1)
     362    {
    332363      track = static_cast<Candidate*>(fTrackInputArray->At(number));
     364
     365      ++fTowerTrackAllHits;
    333366      fTowerTrackArray->Add(track);
     367      if(flags & 2)
     368      {
     369        ++fTowerECalTrackHits;
     370        fTowerECalTrackArray->Add(track);
     371      }
     372      if(flags & 4)
     373      {
     374        ++fTowerHCalTrackHits;
     375        fTowerHCalTrackArray->Add(track);
     376      }
    334377      continue;
    335378    }
     379
     380    ++fTowerAllHits;
     381
     382    // check for ECAL hits in current tower
     383    if(flags & 2)
     384    {
     385      ++fTowerECalHits;
     386      fTowerECalArray->Add(particle);
     387    }
     388
     389    // check for HCAL hits in current tower
     390    if(flags & 4)
     391    {
     392      ++fTowerHCalHits;
     393      fTowerHCalArray->Add(particle);
     394    }
     395
     396    // check for photon and electron hits in current tower
     397    if(flags & 8) ++fTowerPhotonHits;
    336398
    337399    particle = static_cast<Candidate*>(fParticleInputArray->At(number));
     
    345407    fTowerHCalEnergy += hcalEnergy;
    346408
    347     ++fTowerAllHits;
    348409    fTower->AddCandidate(particle);
    349 
    350     // check for neutral hits in current tower
    351     if(flags & 1) ++fTowerNeutralHits;
    352 
    353     // check for photon hits in current tower
    354     if(flags & 2)
    355     {
    356       ++fTowerPhotonHits;
    357       fTowerPhotonArray->Add(particle);
    358     }
    359 
    360     // check for electron hits in current tower
    361     if(flags & 4) ++fTowerElectronHits;
    362410  }
    363411
     
    373421  Double_t energy, pt, eta, phi;
    374422  Double_t ecalEnergy, hcalEnergy;
     423  TIterator *iterator;
    375424
    376425  if(!fTower) return;
     
    409458  if(energy > 0.0)
    410459  {
    411     if((fTowerPhotonHits > 0 || fTowerElectronHits > 0) &&
    412         fTowerTrackHits == 0)
     460    if(fTowerPhotonHits > 0 && fTowerTrackAllHits == 0)
    413461    {
    414462      fPhotonOutputArray->Add(fTower);
     
    419467
    420468  // fill energy flow candidates
    421   if(fTowerTrackHits == fTowerAllHits)
     469  if(fTowerTrackAllHits == fTowerAllHits)
    422470  {
    423471    fItTowerTrackArray->Reset();
     
    427475    }
    428476  }
    429   else if(fTowerTrackHits > 0 &&
    430           fTowerElectronHits == 0 &&
    431           fTowerPhotonHits + fTowerTrackHits == fTowerAllHits)
    432   {
    433     fItTowerTrackArray->Reset();
    434     while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
    435     {
    436       fEFlowTrackOutputArray->Add(track);
    437     }
    438 
    439     if(ecalEnergy > 0.0)
     477  else if(fTowerTrackAllHits > 0 &&
     478          fTowerECalHits + fTowerHCalHits == fTowerAllHits &&
     479          (fTowerECalHits == fTowerECalTrackHits || fTowerHCalHits == fTowerHCalTrackHits))
     480  {
     481    if(fTowerECalHits == fTowerECalTrackHits)
     482    {
     483      fItTowerECalTrackArray->Reset();
     484      while((track = static_cast<Candidate*>(fItTowerECalTrackArray->Next())))
     485      {
     486        fEFlowTrackOutputArray->Add(track);
     487      }
     488      energy = hcalEnergy;
     489      iterator = fItTowerHCalArray;
     490    }
     491
     492    if(fTowerHCalHits == fTowerHCalTrackHits)
     493    {
     494      fItTowerHCalTrackArray->Reset();
     495      while((track = static_cast<Candidate*>(fItTowerHCalTrackArray->Next())))
     496      {
     497        fEFlowTrackOutputArray->Add(track);
     498      }
     499      energy = ecalEnergy;
     500      iterator = fItTowerECalArray;
     501    }
     502
     503    if(energy > 0.0 && iterator)
    440504    {
    441505      DelphesFactory *factory = GetFactory();
     
    444508      tower = factory->NewCandidate();
    445509
    446       fItTowerPhotonArray->Reset();
    447       while((particle = static_cast<Candidate*>(fItTowerPhotonArray->Next())))
     510      iterator->Reset();
     511      while((particle = static_cast<Candidate*>(iterator->Next())))
    448512      {
    449513        tower->AddCandidate(particle);
    450514      }
    451515
    452       pt = ecalEnergy / TMath::CosH(eta);
     516      pt = energy / TMath::CosH(eta);
    453517
    454518      tower->Position.SetPtEtaPhiE(1.0, eta, phi, 0.0);
    455       tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalEnergy);
    456       tower->Eem = ecalEnergy;
     519      tower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
     520      tower->Eem = energy;
    457521      tower->Ehad = 0.0;
    458522
  • modules/Calorimeter.h

    rd60d554 r73e0386  
    4545  Double_t fTowerECalEnergy, fTowerHCalEnergy;
    4646  Double_t fTowerECalNeutralEnergy, fTowerHCalNeutralEnergy;
    47   Int_t fTowerNeutralHits, fTowerPhotonHits, fTowerElectronHits, fTowerTrackHits, fTowerAllHits;
     47  Int_t fTowerPhotonHits, fTowerECalHits, fTowerHCalHits, fTowerAllHits;
     48  Int_t fTowerECalTrackHits, fTowerHCalTrackHits, fTowerTrackAllHits;
    4849
    4950  TFractionMap fFractionMap; //!
     
    7374  TObjArray *fEFlowTowerOutputArray; //!
    7475
    75   TObjArray *fTowerTrackArray; //!
    76   TIterator *fItTowerTrackArray; //!
     76  TObjArray *fTowerECalTrackArray; //!
     77  TIterator *fItTowerECalTrackArray; //!
    7778
    78   TObjArray *fTowerPhotonArray; //!
    79   TIterator *fItTowerPhotonArray; //!
     79  TObjArray *fTowerHCalTrackArray; //!
     80  TIterator *fItTowerHCalTrackArray; //!
     81
     82  TObjArray *fTowerECalArray; //!
     83  TIterator *fItTowerECalArray; //!
     84
     85  TObjArray *fTowerHCalArray; //!
     86  TIterator *fItTowerHCalArray; //!
    8087
    8188  void FinalizeTower();
Note: See TracChangeset for help on using the changeset viewer.