Changeset 2dab783 in git for modules/Calorimeter.cc

Timestamp:

Sep 3, 2013, 5:54:56 PM (11 years ago)

Author:

pavel <pavel@…>

Branches:

ImprovedOutputFile, Timing, dual_readout, llp, master

Children:

2d28bda

Parents:

2eb3782

Message:

new energy flow

File:

• modules/Calorimeter.cc (modified) (14 diffs)

modules/Calorimeter.cc

@@ -43 +43 @@
   fECalResolutionFormula(0), fHCalResolutionFormula(0),
   fItParticleInputArray(0), fItTrackInputArray(0),
-  fTowerECalArray(0), fItTowerECalArray(0),
-  fTowerHCalArray(0), fItTowerHCalArray(0),
-  fTowerTrackArray(0), fItTowerTrackArray(0),
-  fTowerECalTrackArray(0), fItTowerECalTrackArray(0),
-  fTowerHCalTrackArray(0), fItTowerHCalTrackArray(0)
+  fTowerTrackArray(0), fItTowerTrackArray(0)
 {
   fECalResolutionFormula = new DelphesFormula;
   fHCalResolutionFormula = new DelphesFormula;
 
-  fTowerECalArray = new TObjArray;
-  fItTowerECalArray = fTowerECalArray->MakeIterator();
-  fTowerHCalArray = new TObjArray;
-  fItTowerHCalArray = fTowerHCalArray->MakeIterator();
-
   fTowerTrackArray = new TObjArray;
   fItTowerTrackArray = fTowerTrackArray->MakeIterator();
-  fTowerECalTrackArray = new TObjArray;
-  fItTowerECalTrackArray = fTowerECalTrackArray->MakeIterator();
-  fTowerHCalTrackArray = new TObjArray;
-  fItTowerHCalTrackArray = fTowerHCalTrackArray->MakeIterator();
 }
 
@@ -72 +59 @@
   if(fHCalResolutionFormula) delete fHCalResolutionFormula;
 
-  if(fTowerECalArray) delete fTowerECalArray;
-  if(fItTowerECalArray) delete fItTowerECalArray;
-  if(fTowerHCalArray) delete fTowerHCalArray;
-  if(fItTowerHCalArray) delete fItTowerHCalArray;
-
   if(fTowerTrackArray) delete fTowerTrackArray;
   if(fItTowerTrackArray) delete fItTowerTrackArray;
-  if(fTowerECalTrackArray) delete fTowerECalTrackArray;
-  if(fItTowerECalTrackArray) delete fItTowerECalTrackArray;
-  if(fTowerHCalTrackArray) delete fTowerHCalTrackArray;
-  if(fItTowerHCalTrackArray) delete fItTowerHCalTrackArray;}
+}
 
 //------------------------------------------------------------------------------
@@ -179 +158 @@
 void Calorimeter::Finish()
 {
-  vector< vector< Double_t>* >::iterator itPhiBin;
+  vector< vector< Double_t >* >::iterator itPhiBin;
   if(fItParticleInputArray) delete fItParticleInputArray;
   if(fItTrackInputArray) delete fItTrackInputArray;
@@ -211 +190 @@
   DelphesFactory *factory = GetFactory();
   fTowerHits.clear();
-  fECalFractions.clear();
-  fHCalFractions.clear();
+  fTowerECalFractions.clear();
+  fTowerHCalFractions.clear();
+  fTrackECalFractions.clear();
+  fTrackHCalFractions.clear();
 
   // loop over all particles
@@ -233 +214 @@
     hcalFraction = itFractionMap->second.second;
 
-    fECalFractions.push_back(ecalFraction);
-    fHCalFractions.push_back(hcalFraction);
+    fTowerECalFractions.push_back(ecalFraction);
+    fTowerHCalFractions.push_back(hcalFraction);
 
     if(ecalFraction < 1.0E-9 && hcalFraction < 1.0E-9) continue;
@@ -252 +233 @@
 
     flags = 0;
-    flags |= (ecalFraction >= 1.0E-9) << 1;
-    flags |= (hcalFraction >= 1.0E-9) << 2;
-    flags |= (pdgCode == 11 || pdgCode == 22) << 3;
+    flags |= (pdgCode == 11 || pdgCode == 22) << 1;
 
     // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for particle number}
@@ -281 +260 @@
     hcalFraction = itFractionMap->second.second;
 
+    fTrackECalFractions.push_back(ecalFraction);
+    fTrackHCalFractions.push_back(hcalFraction);
+
+    if(ecalFraction < 1.0E-9 && hcalFraction < 1.0E-9) continue;
+
     // find eta bin [1, fEtaBins.size - 1]
     itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), trackPosition.Eta());
@@ -295 +279 @@
 
     flags = 1;
-    flags |= (ecalFraction >= 1.0E-9) << 1;
-    flags |= (hcalFraction >= 1.0E-9) << 2;
 
     // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for track number}
@@ -347 +329 @@
       fTowerHCalEnergy = 0.0;
 
+      fTrackECalEnergy = 0.0;
+      fTrackHCalEnergy = 0.0;
+
+      fTowerTrackHits = 0;
       fTowerPhotonHits = 0;
 
-      fTowerAllHits = 0;
-      fTowerECalHits = 0;
-      fTowerHCalHits = 0;
-
-      fTowerTrackAllHits = 0;
-      fTowerECalTrackHits = 0;
-      fTowerHCalTrackHits = 0;
-
-      fTowerECalArray->Clear();
-      fTowerHCalArray->Clear();
-
       fTowerTrackArray->Clear();
-      fTowerECalTrackArray->Clear();
-      fTowerHCalTrackArray->Clear();
     }
 
@@ -368 +341 @@
     if(flags & 1)
     {
+      ++fTowerTrackHits;
+
       track = static_cast<Candidate*>(fTrackInputArray->At(number));
-
-      ++fTowerTrackAllHits;
+      momentum = track->Momentum;
+
+      ecalEnergy = momentum.E() * fTrackECalFractions[number];
+      hcalEnergy = momentum.E() * fTrackHCalFractions[number];
+
+      fTrackECalEnergy += ecalEnergy;
+      fTrackHCalEnergy += hcalEnergy;
+
       fTowerTrackArray->Add(track);
 
-      // check for track ECAL hits
-      if(flags & 2)
-      {
-        ++fTowerECalTrackHits;
-        fTowerECalTrackArray->Add(track);
-      }
-
-      // check for track HCAL hits
-      if(flags & 4)
-      {
-        ++fTowerHCalTrackHits;
-        fTowerHCalTrackArray->Add(track);
-      }
       continue;
     }
 
-    ++fTowerAllHits;
-
-    // check for ECAL hits
-    if(flags & 2)
-    {
-      ++fTowerECalHits;
-      fTowerECalArray->Add(particle);
-    }
-
-    // check for HCAL hits
-    if(flags & 4)
-    {
-      ++fTowerHCalHits;
-      fTowerHCalArray->Add(particle);
-    }
-
     // check for photon and electron hits in current tower
-    if(flags & 8) ++fTowerPhotonHits;
+    if(flags & 2) ++fTowerPhotonHits;
 
     particle = static_cast<Candidate*>(fParticleInputArray->At(number));
@@ -412 +364 @@
 
     // fill current tower
-    ecalEnergy = momentum.E() * fECalFractions[number];
-    hcalEnergy = momentum.E() * fHCalFractions[number];
+    ecalEnergy = momentum.E() * fTowerECalFractions[number];
+    hcalEnergy = momentum.E() * fTowerHCalFractions[number];
 
     fTowerECalEnergy += ecalEnergy;
@@ -429 +381 @@
 void Calorimeter::FinalizeTower()
 {
-  Candidate *particle, *track, *tower;
+  Candidate *track, *tower;
   Double_t energy, pt, eta, phi;
   Double_t ecalEnergy, hcalEnergy;
-  TIterator *itTowerTrackArray;
 
   if(!fTower) return;
@@ -469 +420 @@
   if(energy > 0.0)
   {
-    if(fTowerPhotonHits > 0 && fTowerTrackAllHits == 0)
+    if(fTowerPhotonHits > 0 && fTowerTrackHits == 0)
     {
       fPhotonOutputArray->Add(fTower);
@@ -478 +429 @@
 
   // fill energy flow candidates
-  if(fTowerTrackAllHits == fTowerAllHits)
-  {
-    fItTowerTrackArray->Reset();
-    while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
-    {
-      fEFlowTrackOutputArray->Add(track);
-    }
-  }
-  else if(fTowerTrackAllHits > 0 &&
-          fTowerECalHits + fTowerHCalHits == fTowerAllHits)
-  {
-    if(fTowerECalHits == fTowerECalTrackHits &&
-       fTowerHCalHits == fTowerHCalTrackHits)
-    {
-      itTowerTrackArray = fItTowerTrackArray;
-    }
-    else if(fTowerECalHits == fTowerECalTrackHits)
-    {
-      itTowerTrackArray = fItTowerECalTrackArray;
-
-      if(hcalEnergy > 0.0)
-      {
-        DelphesFactory *factory = GetFactory();
-
-        // create new tower
-        tower = factory->NewCandidate();
-
-        fItTowerHCalArray->Reset();
-        while((particle = static_cast<Candidate*>(fItTowerHCalArray->Next())))
-        {
-          tower->AddCandidate(particle);
-        }
-
-        pt = hcalEnergy / TMath::CosH(eta);
-
-        tower->Position.SetPtEtaPhiE(1.0, eta, phi, 0.0);
-        tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalEnergy);
-        tower->Eem = 0.0;
-        tower->Ehad = hcalEnergy;
-
-        tower->Edges[0] = fTowerEdges[0];
-        tower->Edges[1] = fTowerEdges[1];
-        tower->Edges[2] = fTowerEdges[2];
-        tower->Edges[3] = fTowerEdges[3];
-
-        fEFlowTowerOutputArray->Add(tower);
-      }
-    }
-    else if(fTowerHCalHits == fTowerHCalTrackHits)
-    {
-      itTowerTrackArray = fItTowerHCalTrackArray;
-
-      if(ecalEnergy > 0.0)
-      {
-        DelphesFactory *factory = GetFactory();
-
-        // create new tower
-        tower = factory->NewCandidate();
-
-        fItTowerECalArray->Reset();
-        while((particle = static_cast<Candidate*>(fItTowerECalArray->Next())))
-        {
-          tower->AddCandidate(particle);
-        }
-
-        pt = ecalEnergy / TMath::CosH(eta);
-
-        tower->Position.SetPtEtaPhiE(1.0, eta, phi, 0.0);
-        tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalEnergy);
-        tower->Eem = ecalEnergy;
-        tower->Ehad = 0.0;
-
-        tower->Edges[0] = fTowerEdges[0];
-        tower->Edges[1] = fTowerEdges[1];
-        tower->Edges[2] = fTowerEdges[2];
-        tower->Edges[3] = fTowerEdges[3];
-
-        fEFlowTowerOutputArray->Add(tower);
-      }
-    }
-    else
-    {
-      itTowerTrackArray = 0;
-      fEFlowTowerOutputArray->Add(fTower);
-    }
-
-    if(itTowerTrackArray)
-    {
-      itTowerTrackArray->Reset();
-      while((track = static_cast<Candidate*>(itTowerTrackArray->Next())))
-      {
-        fEFlowTrackOutputArray->Add(track);
-      }
-    }
-  }
-  else if(energy > 0.0)
-  {
-    fEFlowTowerOutputArray->Add(fTower);
+
+  // save all the tracks as energy flow tracks
+  fItTowerTrackArray->Reset();
+  while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
+  {
+    fEFlowTrackOutputArray->Add(track);
+  }
+
+  ecalEnergy -= fTrackECalEnergy;
+  if(ecalEnergy < 0.0) ecalEnergy = 0.0;
+
+  hcalEnergy -= fTrackHCalEnergy;
+  if(hcalEnergy < 0.0) hcalEnergy = 0.0;
+
+  energy = ecalEnergy + hcalEnergy;
+
+  // save ECAL and/or HCAL energy excess as an energy flow tower
+  if(energy > 0.0)
+  {
+    // create new tower
+    tower = static_cast<Candidate*>(fTower->Clone());
+
+    pt = energy / TMath::CosH(eta);
+
+    tower->Position.SetPtEtaPhiE(1.0, eta, phi, 0.0);
+    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
+    tower->Eem = ecalEnergy;
+    tower->Ehad = hcalEnergy;
+
+    tower->Edges[0] = fTowerEdges[0];
+    tower->Edges[1] = fTowerEdges[1];
+    tower->Edges[2] = fTowerEdges[2];
+    tower->Edges[3] = fTowerEdges[3];
+
+    fEFlowTowerOutputArray->Add(tower);
   }
 }