Changeset df636c74 in git for modules/SimpleCalorimeter.cc

Timestamp:

Oct 8, 2015, 4:39:04 PM (9 years ago)

Author:

Pavel Demin <pavel.demin@…>

Branches:

ImprovedOutputFile, Timing, dual_readout, llp, master

Children:

d67b86b

Parents:

9da65a5

Message:

update eflow in SimpleCalorimeter

File:

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

modules/SimpleCalorimeter.cc

@@ -56 +56 @@
 SimpleCalorimeter::SimpleCalorimeter() :
   fResolutionFormula(0),
-  fItParticleInputArray(0), fItTrackInputArray(0),
-  fTowerTrackArray(0), fItTowerTrackArray(0)
-{
+  fItParticleInputArray(0), fItTrackInputArray(0)
+{
+  Int_t i;
+
   fResolutionFormula = new DelphesFormula;
 
-  fTowerTrackArray = new TObjArray;
-  fItTowerTrackArray = fTowerTrackArray->MakeIterator();
+  for(i = 0; i < 2; ++i)
+  {
+    fTowerTrackArray[i] = new TObjArray;
+    fItTowerTrackArray[i] = fTowerTrackArray[i]->MakeIterator();
+  }
 }
 
@@ -69 +73 @@
 SimpleCalorimeter::~SimpleCalorimeter()
 {
+  Int_t i;
+
   if(fResolutionFormula) delete fResolutionFormula;
 
-  if(fTowerTrackArray) delete fTowerTrackArray;
-  if(fItTowerTrackArray) delete fItTowerTrackArray;
+  for(i = 0; i < 2; ++i)
+  {
+    if(fTowerTrackArray[i]) delete fTowerTrackArray[i];
+    if(fItTowerTrackArray[i]) delete fItTowerTrackArray[i];
+  }
 }
 
@@ -137 +146 @@
   }
 
-/*
-  TFractionMap::iterator itFractionMap;
-  for(itFractionMap = fFractionMap.begin(); itFractionMap != fFractionMap.end(); ++itFractionMap)
-  {
-    cout << itFractionMap->first << "   " << itFractionMap->second.first  << "   " << itFractionMap->second.second << endl;
-  }
-*/
-
   // read min E value for towers to be saved
   fEnergyMin = GetDouble("EnergyMin", 0.0);
@@ -168 +169 @@
   // create output arrays
   fTowerOutputArray = ExportArray(GetString("TowerOutputArray", "towers"));
- 
+
   fEFlowTrackOutputArray = ExportArray(GetString("EFlowTrackOutputArray", "eflowTracks"));
   fEFlowTowerOutputArray = ExportArray(GetString("EFlowTowerOutputArray", "eflowTowers"));
@@ -197 +198 @@
   Double_t fraction;
   Double_t energy;
+  Double_t sigma;
   Int_t pdgCode;
 
@@ -337 +339 @@
 
       fTowerEnergy = 0.0;
-      fTrackEnergy = 0.0;
+
+      fTrackEnergy[0] = 0.0;
+      fTrackEnergy[1] = 0.0;
 
       fTowerTime = 0.0;
@@ -347 +351 @@
       fTowerPhotonHits = 0;
 
-      fTowerTrackArray->Clear();
+      fTowerTrackArray[0]->Clear();
+      fTowerTrackArray[1]->Clear();
     }
 
@@ -361 +366 @@
       energy = momentum.E() * fTrackFractions[number];
 
-      fTrackEnergy += energy;
-
       fTrackTime += TMath::Sqrt(energy)*position.T();
       fTrackTimeWeight += TMath::Sqrt(energy);
 
-      fTowerTrackArray->Add(track);
+      if(fTrackFractions[number] > 1.0E-9)
+      {
+        sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
+        if(sigma/momentum.E() < track->TrackResolution)
+        {
+          fTrackEnergy[0] += energy;
+          fTowerTrackArray[0]->Add(track);
+        }
+        else
+        {
+          fTrackEnergy[1] += energy;
+          fTowerTrackArray[1]->Add(track);
+        }
+      }
+      else if(fTrackFractions[number] < 1.0E-9)
+      {
+        fEFlowTrackOutputArray->Add(track);
+      }
 
       continue;
@@ -397 +417 @@
 void SimpleCalorimeter::FinalizeTower()
 {
-  Candidate *tower, *track;
+  Candidate *tower, *track, *mother;
   Double_t energy, pt, eta, phi;
   Double_t sigma;
   Double_t time;
 
-  Double_t trkSigma, fraction;
-  
-  Int_t pdgCode;
   TLorentzVector momentum;
   TFractionMap::iterator itFractionMap;
- 
+
   if(!fTower) return;
 
@@ -437 +454 @@
 
   fTower->Eem = (!fIsEcal) ? 0 : energy;
-  fTower->Ehad = (fIsEcal) ? 0 : energy; 
+  fTower->Ehad = (fIsEcal) ? 0 : energy;
 
   fTower->Edges[0] = fTowerEdges[0];
@@ -447 +464 @@
   if(energy > 0.0) fTowerOutputArray->Add(fTower);
 
-
-
-  // fill e-flow candidates 
-  fItTowerTrackArray->Reset();
-  while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
-  {
-     momentum = track->Momentum;
-   
-     pdgCode = TMath::Abs(track->PID);
-
-     itFractionMap = fFractionMap.find(pdgCode);
-     if(itFractionMap == fFractionMap.end())
-     {
-       itFractionMap = fFractionMap.find(0);
-     }
-
-     fraction = itFractionMap->second;
-   
-     // charged particle has to deposit either in ECAL or HCAL
-     if(fraction > 1.0E-9)
-     {
-       trkSigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());       
-       if(track->TrackResolution < trkSigma/momentum.E()) 
-       {
-         energy -= momentum.E();
-         fEFlowTrackOutputArray->Add(track);
-       }
-     }
-     //forward all tracks from ECAL to HCAL
-     else if(fIsEcal)
-     {
-       fEFlowTrackOutputArray->Add(track);
-     }
-     //store muons from HCAL 
-     else if(pdgCode == 13)
-     {
-       fEFlowTrackOutputArray->Add(track);
-     }
-  }
-
+  // fill e-flow candidates
+
+  energy -= fTrackEnergy[1];
+
+  fItTowerTrackArray[0]->Reset();
+  while((track = static_cast<Candidate*>(fItTowerTrackArray[0]->Next())))
+  {
+    mother = track;
+    track = static_cast<Candidate*>(track->Clone());
+    track->AddCandidate(mother);
+
+    track->Momentum *= energy/fTrackEnergy[0];
+
+    fEFlowTrackOutputArray->Add(track);
+  }
+
+  fItTowerTrackArray[1]->Reset();
+  while((track = static_cast<Candidate*>(fItTowerTrackArray[1]->Next())))
+  {
+    mother = track;
+    track = static_cast<Candidate*>(track->Clone());
+    track->AddCandidate(mother);
+
+    fEFlowTrackOutputArray->Add(track);
+  }
+
+  if(fTowerTrackArray[0]->GetEntriesFast() > 0) energy = 0.0;
 
   sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, energy);
@@ -499 +503 @@
 
     tower->Eem = (!fIsEcal) ? 0 : energy;
-    tower->Ehad = (fIsEcal) ? 0 : energy; 
+    tower->Ehad = (fIsEcal) ? 0 : energy;
 
     tower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);