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

File:

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

modules/SimpleCalorimeter.cc

@@ -17 +17 @@
  */
 
-
 /** \class SimpleCalorimeter
  *
@@ -33 +32 @@
 #include "classes/DelphesFormula.h"
 
+#include "ExRootAnalysis/ExRootClassifier.h"
+#include "ExRootAnalysis/ExRootFilter.h"
 #include "ExRootAnalysis/ExRootResult.h"
-#include "ExRootAnalysis/ExRootFilter.h"
-#include "ExRootAnalysis/ExRootClassifier.h"
-
+
+#include "TDatabasePDG.h"
+#include "TFormula.h"
+#include "TLorentzVector.h"
 #include "TMath.h"
+#include "TObjArray.h"
+#include "TRandom3.h"
 #include "TString.h"
-#include "TFormula.h"
-#include "TRandom3.h"
-#include "TObjArray.h"
-#include "TDatabasePDG.h"
-#include "TLorentzVector.h"
 
 #include <algorithm>
-#include <stdexcept>
 #include <iostream>
 #include <sstream>
+#include <stdexcept>
 
 using namespace std;
@@ -58 +57 @@
   fItParticleInputArray(0), fItTrackInputArray(0)
 {
-  
+
   fResolutionFormula = new DelphesFormula;
   fTowerTrackArray = new TObjArray;
   fItTowerTrackArray = fTowerTrackArray->MakeIterator();
- 
 }
 
@@ -69 +67 @@
 SimpleCalorimeter::~SimpleCalorimeter()
 {
- 
+
   if(fResolutionFormula) delete fResolutionFormula;
   if(fTowerTrackArray) delete fTowerTrackArray;
   if(fItTowerTrackArray) delete fItTowerTrackArray;
- 
 }
 
@@ -84 +81 @@
   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
@@ -93 +90 @@
   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();
 
@@ -114 +111 @@
   {
     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();
@@ -131 +128 @@
   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;
   }
 
@@ -170 +167 @@
 void SimpleCalorimeter::Finish()
 {
-  vector< vector< Double_t >* >::iterator itPhiBin;
+  vector<vector<Double_t> *>::iterator itPhiBin;
   if(fItParticleInputArray) delete fItParticleInputArray;
   if(fItTrackInputArray) delete fItTrackInputArray;
@@ -197 +194 @@
   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();
@@ -211 +208 @@
   fItParticleInputArray->Reset();
   number = -1;
-  while((particle = static_cast<Candidate*>(fItParticleInputArray->Next())))
+  while((particle = static_cast<Candidate *>(fItParticleInputArray->Next())))
   {
     const TLorentzVector &particlePosition = particle->Position;
@@ -254 +251 @@
   fItTrackInputArray->Reset();
   number = -1;
-  while((track = static_cast<Candidate*>(fItTrackInputArray->Next())))
+  while((track = static_cast<Candidate *>(fItTrackInputArray->Next())))
   {
     const TLorentzVector &trackPosition = track->Position;
@@ -303 +300 @@
     towerHit = (*itTowerHits);
     flags = (towerHit >> 24) & 0x00000000000000FFLL;
-    number = (towerHit) & 0x0000000000FFFFFFLL;
+    number = (towerHit)&0x0000000000FFFFFFLL;
     hitEtaPhi = towerHit >> 32;
 
@@ -324 +321 @@
 
       // 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];
@@ -346 +343 @@
 
       fTowerTrackArray->Clear();
-     }
+    }
 
     // check for track hits
@@ -353 +350 @@
       ++fTowerTrackHits;
 
-      track = static_cast<Candidate*>(fTrackInputArray->At(number));
+      track = static_cast<Candidate *>(fTrackInputArray->At(number));
       momentum = track->Momentum;
       position = track->Position;
@@ -359 +356 @@
       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
       {
@@ -386 +385 @@
     if(flags & 2) ++fTowerPhotonHits;
 
-    particle = static_cast<Candidate*>(fParticleInputArray->At(number));
+    particle = static_cast<Candidate *>(fParticleInputArray->At(number));
     momentum = particle->Momentum;
     position = particle->Position;
@@ -395 +394 @@
     fTowerEnergy += energy;
 
-    fTowerTime += energy*position.T();
+    fTowerTime += energy * position.T();
     fTowerTimeWeight += energy;
 
@@ -410 +409 @@
 {
   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;
 
@@ -425 +424 @@
   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)
@@ -459 +457 @@
   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);
 
@@ -480 +477 @@
     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);
 
@@ -494 +491 @@
     }
   }
-    
+
   // 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);
 
@@ -516 +513 @@
     }
   }
-   
- } 
+}
 
 //------------------------------------------------------------------------------
@@ -527 +523 @@
   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));
+    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