Changes in modules/Calorimeter.cc [01f457a:38bf1ae] in git

File:

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

modules/Calorimeter.cc

@@ -2 +2 @@
  *  Delphes: a framework for fast simulation of a generic collider experiment
  *  Copyright (C) 2012-2014  Universite catholique de Louvain (UCL), Belgium
- *  
+ *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
  *  the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
- *  
+ *
  *  This program is distributed in the hope that it will be useful,
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *  GNU General Public License for more details.
- *  
+ *
  *  You should have received a copy of the GNU General Public License
  *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
@@ -22 +22 @@
  *  Fills calorimeter towers, performs calorimeter resolution smearing,
  *  and creates energy flow objects (tracks, photons, and neutral hadrons).
- *
- *  $Date$
- *  $Revision$
- *
  *
  *  \author P. Demin - UCL, Louvain-la-Neuve
@@ -86 +82 @@
 {
   ExRootConfParam param, paramEtaBins, paramPhiBins, paramFractions;
-  Long_t i, j, k, size, sizeEtaBins, sizePhiBins, sizeFractions;
+  Long_t i, j, k, size, sizeEtaBins, sizePhiBins;
   Double_t ecalFraction, hcalFraction;
   TBinMap::iterator itEtaBin;
@@ -139 +135 @@
   {
     paramFractions = param[i*2 + 1];
-    sizeFractions = paramFractions.GetSize();
 
     ecalFraction = paramFractions[0].GetDouble();
@@ -146 +141 @@
     fFractionMap[param[i*2].GetInt()] = make_pair(ecalFraction, hcalFraction);
   }
+
 /*
   TFractionMap::iterator itFractionMap;
@@ -155 +151 @@
 
   // read min E value for towers to be saved
-  fEcalEnergyMin = GetDouble("EcalTowerMinEnergy", 0.0); 
-  fHcalEnergyMin = GetDouble("HcalTowerMinEnergy", 0.0); 
-  
-  fEcalSigmaMin  = GetDouble("EcalTowerMinSignificance", 0.0); 
-  fHcalSigmaMin  = GetDouble("HcalTowerMinSignificance", 0.0); 
-
- 
+  fECalEnergyMin = GetDouble("ECalEnergyMin", 0.0);
+  fHCalEnergyMin = GetDouble("HCalEnergyMin", 0.0);
+
+  fECalEnergySignificanceMin = GetDouble("ECalEnergySignificanceMin", 0.0);
+  fHCalEnergySignificanceMin = GetDouble("HCalEnergySignificanceMin", 0.0);
+
+  // switch on or off the dithering of the center of calorimeter towers
+  fDitherTowerCenter = GetBool("DitherTowerCenter", true);
+
   // read resolution formulas
   fECalResolutionFormula->Compile(GetString("ECalResolutionFormula", "0"));
@@ -176 +174 @@
   fTowerOutputArray = ExportArray(GetString("TowerOutputArray", "towers"));
   fPhotonOutputArray = ExportArray(GetString("PhotonOutputArray", "photons"));
-  
+
   fEFlowTrackOutputArray = ExportArray(GetString("EFlowTrackOutputArray", "eflowTracks"));
   fEFlowPhotonOutputArray = ExportArray(GetString("EFlowPhotonOutputArray", "eflowPhotons"));
   fEFlowNeutralHadronOutputArray = ExportArray(GetString("EFlowNeutralHadronOutputArray", "eflowNeutralHadrons"));
-
-
 }
 
@@ -366 +362 @@
       fTrackHCalTime = 0.0;
 
-      fTowerECalWeightTime = 0.0; 
-      fTowerHCalWeightTime = 0.0;
-     
+      fTowerECalTimeWeight = 0.0;
+      fTowerHCalTimeWeight = 0.0;
+
       fTowerTrackHits = 0;
       fTowerPhotonHits = 0;
-      
+
       fTowerTrackArray->Clear();
     }
@@ -384 +380 @@
       position = track->Position;
 
-      
+
       ecalEnergy = momentum.E() * fTrackECalFractions[number];
       hcalEnergy = momentum.E() * fTrackHCalFractions[number];
@@ -390 +386 @@
       fTrackECalEnergy += ecalEnergy;
       fTrackHCalEnergy += hcalEnergy;
-      
+
       fTrackECalTime += TMath::Sqrt(ecalEnergy)*position.T();
       fTrackHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
-       
-      fTrackECalWeightTime += TMath::Sqrt(ecalEnergy);
-      fTrackHCalWeightTime += TMath::Sqrt(hcalEnergy);
+
+      fTrackECalTimeWeight += TMath::Sqrt(ecalEnergy);
+      fTrackHCalTimeWeight += TMath::Sqrt(hcalEnergy);
 
       fTowerTrackArray->Add(track);
@@ -401 +397 @@
       continue;
     }
-   
+
     // check for photon and electron hits in current tower
     if(flags & 2) ++fTowerPhotonHits;
-    
+
     particle = static_cast<Candidate*>(fParticleInputArray->At(number));
     momentum = particle->Momentum;
@@ -419 +415 @@
     fTowerHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
 
-    fTowerECalWeightTime += TMath::Sqrt(ecalEnergy);
-    fTowerHCalWeightTime += TMath::Sqrt(hcalEnergy);
-    
+    fTowerECalTimeWeight += TMath::Sqrt(ecalEnergy);
+    fTowerHCalTimeWeight += TMath::Sqrt(hcalEnergy);
+
 
     fTower->AddCandidate(particle);
@@ -441 +437 @@
 
   if(!fTower) return;
-//  cout<<"----------------------"<<endl;
-//  cout<<"Finalize Tower"<<endl;
-//  cout<<""<<endl;
-
 
   ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerECalEnergy);
-
-//  ecalEnergy = gRandom->Gaus(fTowerECalEnergy, ecalSigma);
-//  if(ecalEnergy < 0.0) ecalEnergy = 0.0;
+  hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerHCalEnergy);
 
   ecalEnergy = LogNormal(fTowerECalEnergy, ecalSigma);
-  ecalTime = (fTowerECalWeightTime < 1.0E-09 ) ? 0 : fTowerECalTime/fTowerECalWeightTime;
-
-  hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerHCalEnergy);
-
-//  hcalEnergy = gRandom->Gaus(fTowerHCalEnergy, hcalSigma);
-//  if(hcalEnergy < 0.0) hcalEnergy = 0.0;
-
   hcalEnergy = LogNormal(fTowerHCalEnergy, hcalSigma);
-  hcalTime = (fTowerHCalWeightTime < 1.0E-09 ) ? 0 : fTowerHCalTime/fTowerHCalWeightTime;
-
-  
+
+  ecalTime = (fTowerECalTimeWeight < 1.0E-09 ) ? 0.0 : fTowerECalTime/fTowerECalTimeWeight;
+  hcalTime = (fTowerHCalTimeWeight < 1.0E-09 ) ? 0.0 : fTowerHCalTime/fTowerHCalTimeWeight;
+
   ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy);
   hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy);
 
-  ecalEnergy = (ecalEnergy < fEcalEnergyMin || ecalEnergy < fEcalSigmaMin*ecalSigma) ? 0 : ecalEnergy;
-  hcalEnergy = (hcalEnergy < fHcalEnergyMin || hcalEnergy < fHcalSigmaMin*hcalSigma) ? 0 : hcalEnergy;
+  if(ecalEnergy < fECalEnergyMin || ecalEnergy < fECalEnergySignificanceMin*ecalSigma) ecalEnergy = 0.0;
+  if(hcalEnergy < fHCalEnergyMin || hcalEnergy < fHCalEnergySignificanceMin*hcalSigma) hcalEnergy = 0.0;
 
   energy = ecalEnergy + hcalEnergy;
-  time = (TMath::Sqrt(ecalEnergy)*ecalTime + TMath::Sqrt(hcalEnergy)*hcalTime)/(TMath::Sqrt(ecalEnergy) + TMath::Sqrt(hcalEnergy)); 
-
-//  eta = fTowerEta;
-//  phi = fTowerPhi;
-
-  eta = gRandom->Uniform(fTowerEdges[0], fTowerEdges[1]);
-  phi = gRandom->Uniform(fTowerEdges[2], fTowerEdges[3]);
+  time = (TMath::Sqrt(ecalEnergy)*ecalTime + TMath::Sqrt(hcalEnergy)*hcalTime)/(TMath::Sqrt(ecalEnergy) + TMath::Sqrt(hcalEnergy));
+
+  if(fDitherTowerCenter)
+  {
+    eta = gRandom->Uniform(fTowerEdges[0], fTowerEdges[1]);
+    phi = gRandom->Uniform(fTowerEdges[2], fTowerEdges[3]);
+  }
+  else
+  {
+    eta = fTowerEta;
+    phi = fTowerPhi;
+  }
 
   pt = energy / TMath::CosH(eta);
 
- // fTower->Position.SetXYZT(-time, 0.0, 0.0, time);
   fTower->Position.SetPtEtaPhiE(1.0, eta, phi, time);
   fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
@@ -491 +479 @@
   fTower->Edges[3] = fTowerEdges[3];
 
-  if( energy > 0.0 )
+  if(energy > 0.0)
   {
     if(fTowerPhotonHits > 0 && fTowerTrackHits == 0)
@@ -497 +485 @@
       fPhotonOutputArray->Add(fTower);
     }
-   
+
     fTowerOutputArray->Add(fTower);
   }
@@ -511 +499 @@
 
   ecalEnergy -= fTrackECalEnergy;
-  if(ecalEnergy < fEcalEnergyMin || ecalEnergy < fEcalSigmaMin*fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy)) ecalEnergy = 0.0;
-
   hcalEnergy -= fTrackHCalEnergy;
-  if(hcalEnergy < fHcalEnergyMin || hcalEnergy < fHcalSigmaMin*fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy)) hcalEnergy = 0.0;
+
+  ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy);
+  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;
 
   energy = ecalEnergy + hcalEnergy;
@@ -527 +518 @@
     tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalEnergy);
     tower->Eem = ecalEnergy;
-    tower->Ehad = 0;
+    tower->Ehad = 0.0;
 
     fEFlowPhotonOutputArray->Add(tower);
@@ -539 +530 @@
 
     tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalEnergy);
-    tower->Eem = 0;
+    tower->Eem = 0.0;
     tower->Ehad = hcalEnergy;
 
     fEFlowNeutralHadronOutputArray->Add(tower);
   }
-
-
-
-
 }
 
@@ -561 +548 @@
     a = TMath::Log(mean) - 0.5*b*b;
 
-    return TMath::Exp(a + b*gRandom->Gaus(0, 1));
+    return TMath::Exp(a + b*gRandom->Gaus(0.0, 1.0));
   }
   else