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

File:

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

modules/SimpleCalorimeter.cc

@@ -58 +58 @@
   fItParticleInputArray(0), fItTrackInputArray(0)
 {
-  Int_t i;
-
+  
   fResolutionFormula = new DelphesFormula;
-
-  for(i = 0; i < 2; ++i)
-  {
-    fTowerTrackArray[i] = new TObjArray;
-    fItTowerTrackArray[i] = fTowerTrackArray[i]->MakeIterator();
-  }
+  fTowerTrackArray = new TObjArray;
+  fItTowerTrackArray = fTowerTrackArray->MakeIterator();
+ 
 }
 
@@ -73 +69 @@
 SimpleCalorimeter::~SimpleCalorimeter()
 {
-  Int_t i;
-
+ 
   if(fResolutionFormula) delete fResolutionFormula;
-
-  for(i = 0; i < 2; ++i)
-  {
-    if(fTowerTrackArray[i]) delete fTowerTrackArray[i];
-    if(fItTowerTrackArray[i]) delete fItTowerTrackArray[i];
-  }
+  if(fTowerTrackArray) delete fTowerTrackArray;
+  if(fItTowerTrackArray) delete fItTowerTrackArray;
+ 
 }
 
@@ -199 +191 @@
   Double_t energy;
   Double_t sigma;
+  Double_t energyGuess;
+
   Int_t pdgCode;
 
@@ -340 +334 @@
       fTowerEnergy = 0.0;
 
-      fTrackEnergy[0] = 0.0;
-      fTrackEnergy[1] = 0.0;
+      fTrackEnergy = 0.0;
+      fTrackSigma = 0.0;
 
       fTowerTime = 0.0;
@@ -351 +345 @@
       fTowerPhotonHits = 0;
 
-      fTowerTrackArray[0]->Clear();
-      fTowerTrackArray[1]->Clear();
-    }
+      fTowerTrackArray->Clear();
+     }
 
     // check for track hits
@@ -371 +364 @@
       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);
-        }
+             
+       // compute total charged energy   
+       fTrackEnergy += energy;
+       sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
+       if(sigma/momentum.E() < track->TrackResolution) energyGuess = energy;
+       else energyGuess = momentum.E();
+
+       fTrackSigma += ((track->TrackResolution)*energyGuess)*((track->TrackResolution)*energyGuess);
+       fTowerTrackArray->Add(track);
       }
+       
       else
       {
@@ -403 +395 @@
     fTowerEnergy += energy;
 
-    fTowerTime += TMath::Sqrt(energy)*position.T();
-    fTowerTimeWeight += TMath::Sqrt(energy);
+    fTowerTime += energy*position.T();
+    fTowerTimeWeight += energy;
 
     fTower->AddCandidate(particle);
@@ -418 +410 @@
 {
   Candidate *tower, *track, *mother;
-  Double_t energy, pt, eta, phi;
-  Double_t sigma;
+  Double_t energy,neutralEnergy, pt, eta, phi;
+  Double_t sigma, neutralSigma;
   Double_t time;
+   
+  Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor;
 
   TLorentzVector momentum;
@@ -436 +430 @@
 
   if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0;
+
 
   if(fSmearTowerCenter)
@@ -464 +459 @@
   if(energy > 0.0) fTowerOutputArray->Add(fTower);
 
-  // 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)
+ 
+  // 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)
   {
     // create new photon tower
     tower = static_cast<Candidate*>(fTower->Clone());
-    pt = energy / TMath::CosH(eta);
-
-    tower->Eem = (!fIsEcal) ? 0 : energy;
-    tower->Ehad = (fIsEcal) ? 0 : energy;
+    pt = neutralEnergy / TMath::CosH(eta);
+
+    tower->Eem = (!fIsEcal) ? 0 : neutralEnergy;
+    tower->Ehad = (fIsEcal) ? 0 : neutralEnergy;
+    tower->PID = (fIsEcal) ? 22 : 0;
+ 
+    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, neutralEnergy);
+    fEFlowTowerOutputArray->Add(tower);
     
-    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
+    fItTowerTrackArray->Reset();
+    while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
+    {
+      mother = track;
+      track = static_cast<Candidate*>(track->Clone());
+      track->AddCandidate(mother);
+
+      fEFlowTrackOutputArray->Add(track);
+    }
+  }
     
-    tower->PID = (fIsEcal) ? 22 : 0;
-    
-    fEFlowTowerOutputArray->Add(tower);
-  }
-}
+  // 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);
+    }
+  }
+   
+ } 
 
 //------------------------------------------------------------------------------