Changeset 73e0386 in git for modules/Calorimeter.cc

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

File:

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

modules/Calorimeter.cc

@@ -43 +43 @@
   fECalResolutionFormula(0), fHCalResolutionFormula(0),
   fItParticleInputArray(0), fItTrackInputArray(0),
-  fTowerTrackArray(0), fItTowerTrackArray(0),
-  fTowerPhotonArray(0), fItTowerPhotonArray(0)
+  fTowerECalArray(0), fItTowerECalArray(0),
+  fTowerHCalArray(0), fItTowerHCalArray(0),
+  fTowerECalTrackArray(0), fItTowerECalTrackArray(0),
+  fTowerHCalTrackArray(0), fItTowerHCalTrackArray(0)
 {
   fECalResolutionFormula = new DelphesFormula;
   fHCalResolutionFormula = new DelphesFormula;
-  fTowerTrackArray = new TObjArray;
-  fItTowerTrackArray = fTowerTrackArray->MakeIterator();
-  fTowerPhotonArray = new TObjArray;
-  fItTowerPhotonArray = fTowerPhotonArray->MakeIterator();
+
+  fTowerECalArray = new TObjArray;
+  fItTowerECalArray = fTowerECalArray->MakeIterator();
+  fTowerHCalArray = new TObjArray;
+  fItTowerHCalArray = fTowerHCalArray->MakeIterator();
+
+  fTowerECalTrackArray = new TObjArray;
+  fItTowerECalTrackArray = fTowerECalTrackArray->MakeIterator();
+  fTowerHCalTrackArray = new TObjArray;
+  fItTowerHCalTrackArray = fTowerHCalTrackArray->MakeIterator();
 }
 
@@ -60 +68 @@
   if(fECalResolutionFormula) delete fECalResolutionFormula;
   if(fHCalResolutionFormula) delete fHCalResolutionFormula;
-  if(fTowerTrackArray) delete fTowerTrackArray;
-  if(fItTowerTrackArray) delete fItTowerTrackArray;
-  if(fTowerPhotonArray) delete fTowerPhotonArray;
-  if(fItTowerPhotonArray) delete fItTowerPhotonArray;
-}
+
+  if(fTowerECalArray) delete fTowerECalArray;
+  if(fItTowerECalArray) delete fItTowerECalArray;
+  if(fTowerHCalArray) delete fTowerHCalArray;
+  if(fItTowerHCalArray) delete fItTowerHCalArray;
+
+  if(fTowerECalTrackArray) delete fTowerECalTrackArray;
+  if(fItTowerECalTrackArray) delete fItTowerECalTrackArray;
+  if(fTowerHCalTrackArray) delete fTowerHCalTrackArray;
+  if(fItTowerHCalTrackArray) delete fItTowerHCalTrackArray;}
 
 //------------------------------------------------------------------------------
@@ -233 +246 @@
     phiBin = distance(phiBins->begin(), itPhiBin);
 
-    flags = (particle->Charge == 0);
-    flags |= (pdgCode == 22) << 1;
-    flags |= (pdgCode == 11) << 2;
+    flags = 0;
+    flags |= (ecalFraction >= 1.0E-9) << 1;
+    flags |= (hcalFraction >= 1.0E-9) << 2;
+    flags |= (pdgCode == 11 || pdgCode == 22) << 3;
 
     // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for particle number}
@@ -251 +265 @@
     ++number;
 
+    pdgCode = TMath::Abs(track->PID);
+
+    itFractionMap = fFractionMap.find(pdgCode);
+    if(itFractionMap == fFractionMap.end())
+    {
+      itFractionMap = fFractionMap.find(0);
+    }
+
+    ecalFraction = itFractionMap->second.first;
+    hcalFraction = itFractionMap->second.second;
+
     // find eta bin [1, fEtaBins.size - 1]
     itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), trackPosition.Eta());
@@ -264 +289 @@
     phiBin = distance(phiBins->begin(), itPhiBin);
 
+    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}
-    towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(1) << 27) | Long64_t(number);
+    towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
 
     fTowerHits.push_back(towerHit);
@@ -313 +342 @@
       fTowerHCalEnergy = 0.0;
 
-      fTowerECalNeutralEnergy = 0.0;
-      fTowerHCalNeutralEnergy = 0.0;
-
-      fTowerNeutralHits = 0;
       fTowerPhotonHits = 0;
-      fTowerElectronHits = 0;
-      fTowerTrackHits = 0;
+
       fTowerAllHits = 0;
-
-      fTowerTrackArray->Clear();
-      fTowerPhotonArray->Clear();
+      fTowerECalHits = 0;
+      fTowerHCalHits = 0;
+
+      fTowerTrackAllHits = 0;
+      fTowerECalTrackHits = 0;
+      fTowerHCalTrackHits = 0;
+
+      fTowerECalTrackArray->Clear();
+      fTowerHCalTrackArray->Clear();
+      fTowerECalArray->Clear();
+      fTowerHCalArray->Clear();
     }
 
     // check for track hits
-    if(flags & 8)
-    {
-      ++fTowerTrackHits;
+    if(flags & 1)
+    {
       track = static_cast<Candidate*>(fTrackInputArray->At(number));
+
+      ++fTowerTrackAllHits;
       fTowerTrackArray->Add(track);
+      if(flags & 2)
+      {
+        ++fTowerECalTrackHits;
+        fTowerECalTrackArray->Add(track);
+      }
+      if(flags & 4)
+      {
+        ++fTowerHCalTrackHits;
+        fTowerHCalTrackArray->Add(track);
+      }
       continue;
     }
+
+    ++fTowerAllHits;
+
+    // check for ECAL hits in current tower
+    if(flags & 2)
+    {
+      ++fTowerECalHits;
+      fTowerECalArray->Add(particle);
+    }
+
+    // check for HCAL hits in current tower
+    if(flags & 4)
+    {
+      ++fTowerHCalHits;
+      fTowerHCalArray->Add(particle);
+    }
+
+    // check for photon and electron hits in current tower
+    if(flags & 8) ++fTowerPhotonHits;
 
     particle = static_cast<Candidate*>(fParticleInputArray->At(number));
@@ -345 +407 @@
     fTowerHCalEnergy += hcalEnergy;
 
-    ++fTowerAllHits;
     fTower->AddCandidate(particle);
-
-    // check for neutral hits in current tower
-    if(flags & 1) ++fTowerNeutralHits;
-
-    // check for photon hits in current tower
-    if(flags & 2)
-    {
-      ++fTowerPhotonHits;
-      fTowerPhotonArray->Add(particle);
-    }
-
-    // check for electron hits in current tower
-    if(flags & 4) ++fTowerElectronHits;
   }
 
@@ -373 +421 @@
   Double_t energy, pt, eta, phi;
   Double_t ecalEnergy, hcalEnergy;
+  TIterator *iterator;
 
   if(!fTower) return;
@@ -409 +458 @@
   if(energy > 0.0)
   {
-    if((fTowerPhotonHits > 0 || fTowerElectronHits > 0) &&
-        fTowerTrackHits == 0)
+    if(fTowerPhotonHits > 0 && fTowerTrackAllHits == 0)
     {
       fPhotonOutputArray->Add(fTower);
@@ -419 +467 @@
 
   // fill energy flow candidates
-  if(fTowerTrackHits == fTowerAllHits)
+  if(fTowerTrackAllHits == fTowerAllHits)
   {
     fItTowerTrackArray->Reset();
@@ -427 +475 @@
     }
   }
-  else if(fTowerTrackHits > 0 &&
-          fTowerElectronHits == 0 &&
-          fTowerPhotonHits + fTowerTrackHits == fTowerAllHits)
-  {
-    fItTowerTrackArray->Reset();
-    while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
-    {
-      fEFlowTrackOutputArray->Add(track);
-    }
-
-    if(ecalEnergy > 0.0)
+  else if(fTowerTrackAllHits > 0 &&
+          fTowerECalHits + fTowerHCalHits == fTowerAllHits &&
+          (fTowerECalHits == fTowerECalTrackHits || fTowerHCalHits == fTowerHCalTrackHits))
+  {
+    if(fTowerECalHits == fTowerECalTrackHits)
+    {
+      fItTowerECalTrackArray->Reset();
+      while((track = static_cast<Candidate*>(fItTowerECalTrackArray->Next())))
+      {
+        fEFlowTrackOutputArray->Add(track);
+      }
+      energy = hcalEnergy;
+      iterator = fItTowerHCalArray;
+    }
+
+    if(fTowerHCalHits == fTowerHCalTrackHits)
+    {
+      fItTowerHCalTrackArray->Reset();
+      while((track = static_cast<Candidate*>(fItTowerHCalTrackArray->Next())))
+      {
+        fEFlowTrackOutputArray->Add(track);
+      }
+      energy = ecalEnergy;
+      iterator = fItTowerECalArray;
+    }
+
+    if(energy > 0.0 && iterator)
     {
       DelphesFactory *factory = GetFactory();
@@ -444 +508 @@
       tower = factory->NewCandidate();
 
-      fItTowerPhotonArray->Reset();
-      while((particle = static_cast<Candidate*>(fItTowerPhotonArray->Next())))
+      iterator->Reset();
+      while((particle = static_cast<Candidate*>(iterator->Next())))
       {
         tower->AddCandidate(particle);
       }
 
-      pt = ecalEnergy / TMath::CosH(eta);
+      pt = energy / TMath::CosH(eta);
 
       tower->Position.SetPtEtaPhiE(1.0, eta, phi, 0.0);
-      tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalEnergy);
-      tower->Eem = ecalEnergy;
+      tower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
+      tower->Eem = energy;
       tower->Ehad = 0.0;