Changes in modules/Calorimeter.cc [341014c:298734e] in git

File:

• modules/Calorimeter.cc (modified) (32 diffs)

modules/Calorimeter.cc

@@ -17 +17 @@
  */
 
+
 /** \class Calorimeter
  *
@@ -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)
 {
-
+  
   fECalResolutionFormula = new DelphesFormula;
   fHCalResolutionFormula = new DelphesFormula;
@@ -66 +67 @@
   fHCalTowerTrackArray = new TObjArray;
   fItHCalTowerTrackArray = fHCalTowerTrackArray->MakeIterator();
+  
 }
 
@@ -72 +74 @@
 Calorimeter::~Calorimeter()
 {
-
+  
   if(fECalResolutionFormula) delete fECalResolutionFormula;
   if(fHCalResolutionFormula) delete fHCalResolutionFormula;
@@ -81 +83 @@
   if(fHCalTowerTrackArray) delete fHCalTowerTrackArray;
   if(fItHCalTowerTrackArray) delete fItHCalTowerTrackArray;
+ 
 }
 
@@ -91 +94 @@
   Double_t ecalFraction, hcalFraction;
   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
@@ -100 +103 @@
   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();
 
@@ -121 +124 @@
   {
     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();
@@ -138 +141 @@
   fFractionMap[0] = make_pair(0.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];
 
     ecalFraction = paramFractions[0].GetDouble();
     hcalFraction = paramFractions[1].GetDouble();
 
-    fFractionMap[param[i * 2].GetInt()] = make_pair(ecalFraction, hcalFraction);
+    fFractionMap[param[i*2].GetInt()] = make_pair(ecalFraction, hcalFraction);
   }
 
   // read min E value for timing measurement in ECAL
-  fTimingEnergyMin = GetDouble("TimingEnergyMin", 4.);
+  fTimingEnergyMin = GetDouble("TimingEnergyMin",4.);
   // For timing
   // So far this flag needs to be false
@@ -189 +192 @@
 void Calorimeter::Finish()
 {
-  vector<vector<Double_t> *>::iterator itPhiBin;
+  vector< vector< Double_t >* >::iterator itPhiBin;
   if(fItParticleInputArray) delete fItParticleInputArray;
   if(fItTrackInputArray) delete fItTrackInputArray;
@@ -215 +218 @@
   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();
@@ -231 +234 @@
   fItParticleInputArray->Reset();
   number = -1;
-  while((particle = static_cast<Candidate *>(fItParticleInputArray->Next())))
+  while((particle = static_cast<Candidate*>(fItParticleInputArray->Next())))
   {
     const TLorentzVector &particlePosition = particle->Position;
@@ -277 +280 @@
   fItTrackInputArray->Reset();
   number = -1;
-  while((track = static_cast<Candidate *>(fItTrackInputArray->Next())))
+  while((track = static_cast<Candidate*>(fItTrackInputArray->Next())))
   {
     const TLorentzVector &trackPosition = track->Position;
@@ -328 +331 @@
     towerHit = (*itTowerHits);
     flags = (towerHit >> 24) & 0x00000000000000FFLL;
-    number = (towerHit)&0x0000000000FFFFFFLL;
+    number = (towerHit) & 0x0000000000FFFFFFLL;
     hitEtaPhi = towerHit >> 32;
 
@@ -349 +352 @@
 
       // 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];
@@ -362 +365 @@
       fECalTrackEnergy = 0.0;
       fHCalTrackEnergy = 0.0;
-
+      
       fECalTrackSigma = 0.0;
       fHCalTrackSigma = 0.0;
-
+      
       fTowerTrackHits = 0;
       fTowerPhotonHits = 0;
@@ -371 +374 @@
       fECalTowerTrackArray->Clear();
       fHCalTowerTrackArray->Clear();
+    
     }
 
@@ -378 +382 @@
       ++fTowerTrackHits;
 
-      track = static_cast<Candidate *>(fTrackInputArray->At(number));
+      track = static_cast<Candidate*>(fTrackInputArray->At(number));
       momentum = track->Momentum;
       position = track->Position;
@@ -396 +400 @@
       {
         fECalTrackEnergy += ecalEnergy;
-        ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
-        if(ecalSigma / momentum.E() < track->TrackResolution)
-          energyGuess = ecalEnergy;
-        else
-          energyGuess = momentum.E();
-
-        fECalTrackSigma += (track->TrackResolution) * energyGuess * (track->TrackResolution) * energyGuess;
+        ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());        
+        if(ecalSigma/momentum.E() < track->TrackResolution) energyGuess = ecalEnergy;        
+        else energyGuess = momentum.E();
+
+        fECalTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess;
         fECalTowerTrackArray->Add(track);
       }
-
+     
       else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] > 1.0E-9)
       {
         fHCalTrackEnergy += hcalEnergy;
         hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
-        if(hcalSigma / momentum.E() < track->TrackResolution)
-          energyGuess = hcalEnergy;
-        else
-          energyGuess = momentum.E();
-
-        fHCalTrackSigma += (track->TrackResolution) * energyGuess * (track->TrackResolution) * energyGuess;
+        if(hcalSigma/momentum.E() < track->TrackResolution) energyGuess = hcalEnergy;
+        else energyGuess = momentum.E();
+
+        fHCalTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess;
         fHCalTowerTrackArray->Add(track);
       }
-
+      
       else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
       {
@@ -430 +430 @@
     if(flags & 2) ++fTowerPhotonHits;
 
-    particle = static_cast<Candidate *>(fParticleInputArray->At(number));
+    particle = static_cast<Candidate*>(fParticleInputArray->At(number));
     momentum = particle->Momentum;
     position = particle->Position;
@@ -443 +443 @@
     if(ecalEnergy > fTimingEnergyMin && fTower)
     {
-      if(abs(particle->PID) != 11 || !fElectronsFromTrack)
+      if (abs(particle->PID) != 11 || !fElectronsFromTrack)
       {
         fTower->ECalEnergyTimePairs.push_back(make_pair<Float_t, Float_t>(ecalEnergy, particle->Position.T()));
@@ -464 +464 @@
   Double_t ecalEnergy, hcalEnergy;
   Double_t ecalNeutralEnergy, hcalNeutralEnergy;
-
+  
   Double_t ecalSigma, hcalSigma;
   Double_t ecalNeutralSigma, hcalNeutralSigma;
 
   Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor;
-
+  
   TLorentzVector momentum;
   TFractionMap::iterator itFractionMap;
@@ -486 +486 @@
   hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy);
 
-  if(ecalEnergy < fECalEnergyMin || ecalEnergy < fECalEnergySignificanceMin * ecalSigma) ecalEnergy = 0.0;
-  if(hcalEnergy < fHCalEnergyMin || hcalEnergy < fHCalEnergySignificanceMin * hcalSigma) hcalEnergy = 0.0;
+  if(ecalEnergy < fECalEnergyMin || ecalEnergy < fECalEnergySignificanceMin*ecalSigma) ecalEnergy = 0.0;
+  if(hcalEnergy < fHCalEnergyMin || hcalEnergy < fHCalEnergySignificanceMin*hcalSigma) hcalEnergy = 0.0;
 
   energy = ecalEnergy + hcalEnergy;
@@ -519 +519 @@
   if(sumWeight > 0.0)
   {
-    fTower->Position.SetPtEtaPhiE(1.0, eta, phi, sumWeightedTime / sumWeight);
+    fTower->Position.SetPtEtaPhiE(1.0, eta, phi, sumWeightedTime/sumWeight);
   }
   else
@@ -525 +525 @@
     fTower->Position.SetPtEtaPhiE(1.0, eta, phi, 999999.9);
   }
+
 
   fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
@@ -544 +545 @@
     fTowerOutputArray->Add(fTower);
   }
-
+  
   // fill energy flow candidates
   fECalTrackSigma = TMath::Sqrt(fECalTrackSigma);
@@ -550 +551 @@
 
   //compute neutral excesses
-  ecalNeutralEnergy = max((ecalEnergy - fECalTrackEnergy), 0.0);
-  hcalNeutralEnergy = max((hcalEnergy - fHCalTrackEnergy), 0.0);
-
-  ecalNeutralSigma = ecalNeutralEnergy / TMath::Sqrt(fECalTrackSigma * fECalTrackSigma + ecalSigma * ecalSigma);
-  hcalNeutralSigma = hcalNeutralEnergy / TMath::Sqrt(fHCalTrackSigma * fHCalTrackSigma + hcalSigma * hcalSigma);
-
-  // if ecal neutral excess is significant, simply create neutral EflowPhoton tower and clone each track into eflowtrack
+  ecalNeutralEnergy = max( (ecalEnergy - fECalTrackEnergy) , 0.0);
+  hcalNeutralEnergy = max( (hcalEnergy - fHCalTrackEnergy) , 0.0);
+  
+  ecalNeutralSigma = ecalNeutralEnergy / TMath::Sqrt(fECalTrackSigma*fECalTrackSigma + ecalSigma*ecalSigma);
+  hcalNeutralSigma = hcalNeutralEnergy / TMath::Sqrt(fHCalTrackSigma*fHCalTrackSigma + hcalSigma*hcalSigma);
+  
+   // if ecal neutral excess is significant, simply create neutral EflowPhoton tower and clone each track into eflowtrack
   if(ecalNeutralEnergy > fECalEnergyMin && ecalNeutralSigma > fECalEnergySignificanceMin)
   {
     // create new photon tower
-    tower = static_cast<Candidate *>(fTower->Clone());
-    pt = ecalNeutralEnergy / TMath::CosH(eta);
-
+    tower = static_cast<Candidate*>(fTower->Clone());
+    pt =  ecalNeutralEnergy / TMath::CosH(eta);
+    
     tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalNeutralEnergy);
     tower->Eem = ecalNeutralEnergy;
     tower->Ehad = 0.0;
     tower->PID = 22;
-
+    
     fEFlowPhotonOutputArray->Add(tower);
-
+   
     //clone tracks
     fItECalTowerTrackArray->Reset();
-    while((track = static_cast<Candidate *>(fItECalTowerTrackArray->Next())))
+    while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next())))
     {
       mother = track;
-      track = static_cast<Candidate *>(track->Clone());
+      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(fECalTrackEnergy > 0.0)
   {
-    weightTrack = (fECalTrackSigma > 0.0) ? 1 / (fECalTrackSigma * fECalTrackSigma) : 0.0;
-    weightCalo = (ecalSigma > 0.0) ? 1 / (ecalSigma * ecalSigma) : 0.0;
-
-    bestEnergyEstimate = (weightTrack * fECalTrackEnergy + weightCalo * ecalEnergy) / (weightTrack + weightCalo);
-    rescaleFactor = bestEnergyEstimate / fECalTrackEnergy;
+    weightTrack = (fECalTrackSigma > 0.0) ? 1 / (fECalTrackSigma*fECalTrackSigma) : 0.0;
+    weightCalo  = (ecalSigma > 0.0) ? 1 / (ecalSigma*ecalSigma) : 0.0;
+  
+    bestEnergyEstimate = (weightTrack*fECalTrackEnergy + weightCalo*ecalEnergy) / (weightTrack + weightCalo); 
+    rescaleFactor = bestEnergyEstimate/fECalTrackEnergy;
 
     //rescale tracks
     fItECalTowerTrackArray->Reset();
-    while((track = static_cast<Candidate *>(fItECalTowerTrackArray->Next())))
-    {
+    while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next())))
+    {  
       mother = track;
-      track = static_cast<Candidate *>(track->Clone());
+      track = static_cast<Candidate*>(track->Clone());
       track->AddCandidate(mother);
 
@@ -604 +606 @@
     }
   }
+
 
   // if hcal neutral excess is significant, simply create neutral EflowNeutralHadron tower and clone each track into eflowtrack
@@ -609 +612 @@
   {
     // create new photon tower
-    tower = static_cast<Candidate *>(fTower->Clone());
-    pt = hcalNeutralEnergy / TMath::CosH(eta);
-
+    tower = static_cast<Candidate*>(fTower->Clone());
+    pt =  hcalNeutralEnergy / TMath::CosH(eta);
+    
     tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalNeutralEnergy);
     tower->Ehad = hcalNeutralEnergy;
     tower->Eem = 0.0;
-
+    
     fEFlowNeutralHadronOutputArray->Add(tower);
-
+   
     //clone tracks
     fItHCalTowerTrackArray->Reset();
-    while((track = static_cast<Candidate *>(fItHCalTowerTrackArray->Next())))
+    while((track = static_cast<Candidate*>(fItHCalTowerTrackArray->Next())))
     {
       mother = track;
-      track = static_cast<Candidate *>(track->Clone());
+      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(fHCalTrackEnergy > 0.0)
   {
-    weightTrack = (fHCalTrackSigma > 0.0) ? 1 / (fHCalTrackSigma * fHCalTrackSigma) : 0.0;
-    weightCalo = (hcalSigma > 0.0) ? 1 / (hcalSigma * hcalSigma) : 0.0;
-
-    bestEnergyEstimate = (weightTrack * fHCalTrackEnergy + weightCalo * hcalEnergy) / (weightTrack + weightCalo);
+    weightTrack = (fHCalTrackSigma > 0.0) ? 1 / (fHCalTrackSigma*fHCalTrackSigma) : 0.0;
+    weightCalo  = (hcalSigma > 0.0) ? 1 / (hcalSigma*hcalSigma) : 0.0;
+  
+    bestEnergyEstimate = (weightTrack*fHCalTrackEnergy + weightCalo*hcalEnergy) / (weightTrack + weightCalo); 
     rescaleFactor = bestEnergyEstimate / fHCalTrackEnergy;
 
     //rescale tracks
     fItHCalTowerTrackArray->Reset();
-    while((track = static_cast<Candidate *>(fItHCalTowerTrackArray->Next())))
-    {
+    while((track = static_cast<Candidate*>(fItHCalTowerTrackArray->Next())))
+    {  
       mother = track;
-      track = static_cast<Candidate *>(track->Clone());
+      track = static_cast<Candidate*>(track->Clone());
       track->AddCandidate(mother);
 
@@ -652 +656 @@
     }
   }
+  
+  
 }
 
@@ -662 +668 @@
   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