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

File:

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

modules/SimpleCalorimeter.cc

@@ -17 +17 @@
  */
 
+
 /** \class SimpleCalorimeter
  *
@@ -32 +33 @@
 #include "classes/DelphesFormula.h"
 
+#include "ExRootAnalysis/ExRootResult.h"
+#include "ExRootAnalysis/ExRootFilter.h"
 #include "ExRootAnalysis/ExRootClassifier.h"
-#include "ExRootAnalysis/ExRootFilter.h"
-#include "ExRootAnalysis/ExRootResult.h"
-
+
+#include "TMath.h"
+#include "TString.h"
+#include "TFormula.h"
+#include "TRandom3.h"
+#include "TObjArray.h"
 #include "TDatabasePDG.h"
-#include "TFormula.h"
 #include "TLorentzVector.h"
-#include "TMath.h"
-#include "TObjArray.h"
-#include "TRandom3.h"
-#include "TString.h"
 
 #include <algorithm>
+#include <stdexcept>
 #include <iostream>
 #include <sstream>
-#include <stdexcept>
 
 using namespace std;
@@ -57 +58 @@
   fItParticleInputArray(0), fItTrackInputArray(0)
 {
-
+  
   fResolutionFormula = new DelphesFormula;
   fTowerTrackArray = new TObjArray;
   fItTowerTrackArray = fTowerTrackArray->MakeIterator();
+ 
 }
 
@@ -67 +69 @@
 SimpleCalorimeter::~SimpleCalorimeter()
 {
-
+ 
   if(fResolutionFormula) delete fResolutionFormula;
   if(fTowerTrackArray) delete fTowerTrackArray;
   if(fItTowerTrackArray) delete fItTowerTrackArray;
+ 
 }
 
@@ -81 +84 @@
   Double_t fraction;
   TBinMap::iterator itEtaBin;
-  set<Double_t>::iterator itPhiBin;
-  vector<Double_t> *phiBins;
+  set< Double_t >::iterator itPhiBin;
+  vector< Double_t > *phiBins;
 
   // read eta and phi bins
@@ -90 +93 @@
   fEtaBins.clear();
   fPhiBins.clear();
-  for(i = 0; i < size / 2; ++i)
-  {
-    paramEtaBins = param[i * 2];
+  for(i = 0; i < size/2; ++i)
+  {
+    paramEtaBins = param[i*2];
     sizeEtaBins = paramEtaBins.GetSize();
-    paramPhiBins = param[i * 2 + 1];
+    paramPhiBins = param[i*2 + 1];
     sizePhiBins = paramPhiBins.GetSize();
 
@@ -111 +114 @@
   {
     fEtaBins.push_back(itEtaBin->first);
-    phiBins = new vector<double>(itEtaBin->second.size());
+    phiBins = new vector< double >(itEtaBin->second.size());
     fPhiBins.push_back(phiBins);
     phiBins->clear();
@@ -128 +131 @@
   fFractionMap[0] = 1.0;
 
-  for(i = 0; i < size / 2; ++i)
-  {
-    paramFractions = param[i * 2 + 1];
+  for(i = 0; i < size/2; ++i)
+  {
+    paramFractions = param[i*2 + 1];
     fraction = paramFractions[0].GetDouble();
-    fFractionMap[param[i * 2].GetInt()] = fraction;
+    fFractionMap[param[i*2].GetInt()] = fraction;
   }
 
@@ -167 +170 @@
 void SimpleCalorimeter::Finish()
 {
-  vector<vector<Double_t> *>::iterator itPhiBin;
+  vector< vector< Double_t >* >::iterator itPhiBin;
   if(fItParticleInputArray) delete fItParticleInputArray;
   if(fItTrackInputArray) delete fItTrackInputArray;
@@ -194 +197 @@
   TFractionMap::iterator itFractionMap;
 
-  vector<Double_t>::iterator itEtaBin;
-  vector<Double_t>::iterator itPhiBin;
-  vector<Double_t> *phiBins;
-
-  vector<Long64_t>::iterator itTowerHits;
+  vector< Double_t >::iterator itEtaBin;
+  vector< Double_t >::iterator itPhiBin;
+  vector< Double_t > *phiBins;
+
+  vector< Long64_t >::iterator itTowerHits;
 
   DelphesFactory *factory = GetFactory();
@@ -208 +211 @@
   fItParticleInputArray->Reset();
   number = -1;
-  while((particle = static_cast<Candidate *>(fItParticleInputArray->Next())))
+  while((particle = static_cast<Candidate*>(fItParticleInputArray->Next())))
   {
     const TLorentzVector &particlePosition = particle->Position;
@@ -251 +254 @@
   fItTrackInputArray->Reset();
   number = -1;
-  while((track = static_cast<Candidate *>(fItTrackInputArray->Next())))
+  while((track = static_cast<Candidate*>(fItTrackInputArray->Next())))
   {
     const TLorentzVector &trackPosition = track->Position;
@@ -300 +303 @@
     towerHit = (*itTowerHits);
     flags = (towerHit >> 24) & 0x00000000000000FFLL;
-    number = (towerHit)&0x0000000000FFFFFFLL;
+    number = (towerHit) & 0x0000000000FFFFFFLL;
     hitEtaPhi = towerHit >> 32;
 
@@ -321 +324 @@
 
       // calculate eta and phi of the tower's center
-      fTowerEta = 0.5 * (fEtaBins[etaBin - 1] + fEtaBins[etaBin]);
-      fTowerPhi = 0.5 * ((*phiBins)[phiBin - 1] + (*phiBins)[phiBin]);
+      fTowerEta = 0.5*(fEtaBins[etaBin - 1] + fEtaBins[etaBin]);
+      fTowerPhi = 0.5*((*phiBins)[phiBin - 1] + (*phiBins)[phiBin]);
 
       fTowerEdges[0] = fEtaBins[etaBin - 1];
@@ -343 +346 @@
 
       fTowerTrackArray->Clear();
-    }
+     }
 
     // check for track hits
@@ -350 +353 @@
       ++fTowerTrackHits;
 
-      track = static_cast<Candidate *>(fTrackInputArray->At(number));
+      track = static_cast<Candidate*>(fTrackInputArray->At(number));
       momentum = track->Momentum;
       position = track->Position;
@@ -356 +359 @@
       energy = momentum.E() * fTrackFractions[number];
 
-      fTrackTime += TMath::Sqrt(energy) * position.T();
+      fTrackTime += TMath::Sqrt(energy)*position.T();
       fTrackTimeWeight += TMath::Sqrt(energy);
 
       if(fTrackFractions[number] > 1.0E-9)
       {
-
-        // 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);
+             
+       // 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
       {
@@ -385 +386 @@
     if(flags & 2) ++fTowerPhotonHits;
 
-    particle = static_cast<Candidate *>(fParticleInputArray->At(number));
+    particle = static_cast<Candidate*>(fParticleInputArray->At(number));
     momentum = particle->Momentum;
     position = particle->Position;
@@ -394 +395 @@
     fTowerEnergy += energy;
 
-    fTowerTime += energy * position.T();
+    fTowerTime += energy*position.T();
     fTowerTimeWeight += energy;
 
@@ -409 +410 @@
 {
   Candidate *tower, *track, *mother;
-  Double_t energy, neutralEnergy, pt, eta, phi;
+  Double_t energy,neutralEnergy, pt, eta, phi;
   Double_t sigma, neutralSigma;
   Double_t time;
-
+   
   Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor;
 
@@ -424 +425 @@
   energy = LogNormal(fTowerEnergy, sigma);
 
-  time = (fTowerTimeWeight < 1.0E-09) ? 0.0 : fTowerTime / fTowerTimeWeight;
+  time = (fTowerTimeWeight < 1.0E-09 ) ? 0.0 : fTowerTime/fTowerTimeWeight;
 
   sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, energy);
 
-  if(energy < fEnergyMin || energy < fEnergySignificanceMin * sigma) energy = 0.0;
+  if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0;
+
 
   if(fSmearTowerCenter)
@@ -457 +459 @@
   if(energy > 0.0) fTowerOutputArray->Add(fTower);
 
+ 
   // e-flow candidates
 
   //compute neutral excess
-
+  
   fTrackSigma = TMath::Sqrt(fTrackSigma);
-  neutralEnergy = max((energy - fTrackEnergy), 0.0);
-
+  neutralEnergy = max( (energy - fTrackEnergy) , 0.0);
+  
   //compute sigma_trk total
-  neutralSigma = neutralEnergy / TMath::Sqrt(fTrackSigma * fTrackSigma + sigma * sigma);
-
+  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());
+    tower = static_cast<Candidate*>(fTower->Clone());
     pt = neutralEnergy / TMath::CosH(eta);
 
@@ -477 +480 @@
     tower->Ehad = (fIsEcal) ? 0 : neutralEnergy;
     tower->PID = (fIsEcal) ? 22 : 0;
-
+ 
     tower->Momentum.SetPtEtaPhiE(pt, eta, phi, neutralEnergy);
     fEFlowTowerOutputArray->Add(tower);
-
+    
     fItTowerTrackArray->Reset();
-    while((track = static_cast<Candidate *>(fItTowerTrackArray->Next())))
+    while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
     {
       mother = track;
-      track = static_cast<Candidate *>(track->Clone());
+      track = static_cast<Candidate*>(track->Clone());
       track->AddCandidate(mother);
 
@@ -491 +494 @@
     }
   }
-
+    
   // 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;
-
+  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())))
+    while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
     {
       mother = track;
-      track = static_cast<Candidate *>(track->Clone());
+      track = static_cast<Candidate*>(track->Clone());
       track->AddCandidate(mother);
 
@@ -513 +516 @@
     }
   }
+   
+ } 
+
+//------------------------------------------------------------------------------
+
+Double_t SimpleCalorimeter::LogNormal(Double_t mean, Double_t sigma)
+{
+  Double_t a, b;
+
+  if(mean > 0.0)
+  {
+    b = TMath::Sqrt(TMath::Log((1.0 + (sigma*sigma)/(mean*mean))));
+    a = TMath::Log(mean) - 0.5*b*b;
+
+    return TMath::Exp(a + b*gRandom->Gaus(0.0, 1.0));
+  }
+  else
+  {
+    return 0.0;
+  }
 }
-
-//------------------------------------------------------------------------------
-
-Double_t SimpleCalorimeter::LogNormal(Double_t mean, Double_t sigma)
-{
-  Double_t a, b;
-
-  if(mean > 0.0)
-  {
-    b = TMath::Sqrt(TMath::Log((1.0 + (sigma * sigma) / (mean * mean))));
-    a = TMath::Log(mean) - 0.5 * b * b;
-
-    return TMath::Exp(a + b * gRandom->Gaus(0.0, 1.0));
-  }
-  else
-  {
-    return 0.0;
-  }
-}