Changes in modules/SimpleCalorimeter.cc [e0f8f99:f8299bc] in git

File:

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

modules/SimpleCalorimeter.cc

@@ -58 +58 @@
   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];
+  }
 }
 
@@ -190 +198 @@
   Double_t fraction;
   Double_t energy;
+  Double_t sigma;
   Int_t pdgCode;
 
@@ -331 +340 @@
       fTowerEnergy = 0.0;
 
-      fTrackEnergy = 0.0;
-      fTrackSigma = 0.0;
+      fTrackEnergy[0] = 0.0;
+      fTrackEnergy[1] = 0.0;
 
       fTowerTime = 0.0;
@@ -342 +351 @@
       fTowerPhotonHits = 0;
 
-      fTowerTrackArray->Clear();
-     }
+      fTowerTrackArray[0]->Clear();
+      fTowerTrackArray[1]->Clear();
+    }
 
     // check for track hits
@@ -361 +371 @@
       if(fTrackFractions[number] > 1.0E-9)
       {
-             
-       // compute total charged energy   
-       fTrackEnergy += energy;
-       fTrackSigma += ((track->TrackResolution)*momentum.E())*((track->TrackResolution)*momentum.E());
-       
-       fTowerTrackArray->Add(track);
-      
+        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
       {
@@ -390 +403 @@
     fTowerEnergy += energy;
 
-    fTowerTime += energy*position.T();
-    fTowerTimeWeight += energy;
+    fTowerTime += TMath::Sqrt(energy)*position.T();
+    fTowerTimeWeight += TMath::Sqrt(energy);
 
     fTower->AddCandidate(particle);
@@ -405 +418 @@
 {
   Candidate *tower, *track, *mother;
-  Double_t energy,neutralEnergy, pt, eta, phi;
-  Double_t sigma, neutralSigma;
+  Double_t energy, pt, eta, phi;
+  Double_t sigma;
   Double_t time;
-   
-  Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor;
 
   TLorentzVector momentum;
@@ -425 +436 @@
 
   if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0;
-
 
   if(fSmearTowerCenter)
@@ -454 +464 @@
   if(energy > 0.0) fTowerOutputArray->Add(fTower);
 
- 
-  // e-flow candidates
-
-  //compute neutral excess
-  
-  fTrackSigma = TMath::Sqrt(fTrackSigma);
-  neutralEnergy = max( (energy - fTrackEnergy) , 0.0);
-  
-  //compute sigma_trk total
-  neutralSigma = neutralEnergy / TMath::Sqrt(fTrackSigma*fTrackSigma+ sigma*sigma);
-   
-  // if neutral excess is significant, simply create neutral Eflow tower and clone each track into eflowtrack
-  if(neutralEnergy > fEnergyMin && neutralSigma > fEnergySignificanceMin)
+  // 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);
+  if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0;
+
+  // save energy excess as an energy flow tower
+  if(energy > 0.0)
   {
     // create new photon tower
     tower = static_cast<Candidate*>(fTower->Clone());
-    pt = neutralEnergy / TMath::CosH(eta);
-
-    tower->Eem = (!fIsEcal) ? 0 : neutralEnergy;
-    tower->Ehad = (fIsEcal) ? 0 : neutralEnergy;
+    pt = energy / TMath::CosH(eta);
+
+    tower->Eem = (!fIsEcal) ? 0 : energy;
+    tower->Ehad = (fIsEcal) ? 0 : energy;
+    
+    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
+    
     tower->PID = (fIsEcal) ? 22 : 0;
- 
-    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, neutralEnergy);
+    
     fEFlowTowerOutputArray->Add(tower);
-    
-    fItTowerTrackArray->Reset();
-    while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
-    {
-      mother = track;
-      track = static_cast<Candidate*>(track->Clone());
-      track->AddCandidate(mother);
-
-      fEFlowTrackOutputArray->Add(track);
-    }
-  }
-    
-  // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the calorimeter and tracking
-  else if (fTrackEnergy > 0.0)
-  {
-    weightTrack = (fTrackSigma > 0.0) ? 1 / (fTrackSigma*fTrackSigma) : 0.0;
-    weightCalo  = (sigma > 0.0) ? 1 / (sigma*sigma) : 0.0;
-  
-    bestEnergyEstimate = (weightTrack*fTrackEnergy + weightCalo*energy) / (weightTrack + weightCalo); 
-    rescaleFactor = bestEnergyEstimate/fTrackEnergy;
-   
-    fItTowerTrackArray->Reset();
-    while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
-    {
-      mother = track;
-      track = static_cast<Candidate*>(track->Clone());
-      track->AddCandidate(mother);
-
-      track->Momentum *= rescaleFactor;
-
-      fEFlowTrackOutputArray->Add(track);
-    }
-  }
-   
- } 
+  }
+}
 
 //------------------------------------------------------------------------------