Changes in modules/Calorimeter.cc [38bf1ae:01f457a] 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
@@ -82 +86 @@
 {
   ExRootConfParam param, paramEtaBins, paramPhiBins, paramFractions;
-  Long_t i, j, k, size, sizeEtaBins, sizePhiBins;
+  Long_t i, j, k, size, sizeEtaBins, sizePhiBins, sizeFractions;
   Double_t ecalFraction, hcalFraction;
   TBinMap::iterator itEtaBin;
@@ -135 +139 @@
   {
     paramFractions = param[i*2 + 1];
+    sizeFractions = paramFractions.GetSize();
 
     ecalFraction = paramFractions[0].GetDouble();
@@ -141 +146 @@
     fFractionMap[param[i*2].GetInt()] = make_pair(ecalFraction, hcalFraction);
   }
-
 /*
   TFractionMap::iterator itFractionMap;
@@ -151 +155 @@
 
   // read min E value for towers to be saved
-  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);
-
+  fEcalEnergyMin = GetDouble("EcalTowerMinEnergy", 0.0); 
+  fHcalEnergyMin = GetDouble("HcalTowerMinEnergy", 0.0); 
+  
+  fEcalSigmaMin  = GetDouble("EcalTowerMinSignificance", 0.0); 
+  fHcalSigmaMin  = GetDouble("HcalTowerMinSignificance", 0.0); 
+
+ 
   // read resolution formulas
   fECalResolutionFormula->Compile(GetString("ECalResolutionFormula", "0"));
@@ -174 +176 @@
   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"));
+
+
 }
 
@@ -362 +366 @@
       fTrackHCalTime = 0.0;
 
-      fTowerECalTimeWeight = 0.0;
-      fTowerHCalTimeWeight = 0.0;
-
+      fTowerECalWeightTime = 0.0; 
+      fTowerHCalWeightTime = 0.0;
+     
       fTowerTrackHits = 0;
       fTowerPhotonHits = 0;
-
+      
       fTowerTrackArray->Clear();
     }
@@ -380 +384 @@
       position = track->Position;
 
-
+      
       ecalEnergy = momentum.E() * fTrackECalFractions[number];
       hcalEnergy = momentum.E() * fTrackHCalFractions[number];
@@ -386 +390 @@
       fTrackECalEnergy += ecalEnergy;
       fTrackHCalEnergy += hcalEnergy;
-
+      
       fTrackECalTime += TMath::Sqrt(ecalEnergy)*position.T();
       fTrackHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
-
-      fTrackECalTimeWeight += TMath::Sqrt(ecalEnergy);
-      fTrackHCalTimeWeight += TMath::Sqrt(hcalEnergy);
+       
+      fTrackECalWeightTime += TMath::Sqrt(ecalEnergy);
+      fTrackHCalWeightTime += TMath::Sqrt(hcalEnergy);
 
       fTowerTrackArray->Add(track);
@@ -397 +401 @@
       continue;
     }
-
+   
     // check for photon and electron hits in current tower
     if(flags & 2) ++fTowerPhotonHits;
-
+    
     particle = static_cast<Candidate*>(fParticleInputArray->At(number));
     momentum = particle->Momentum;
@@ -415 +419 @@
     fTowerHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
 
-    fTowerECalTimeWeight += TMath::Sqrt(ecalEnergy);
-    fTowerHCalTimeWeight += TMath::Sqrt(hcalEnergy);
-
+    fTowerECalWeightTime += TMath::Sqrt(ecalEnergy);
+    fTowerHCalWeightTime += TMath::Sqrt(hcalEnergy);
+    
 
     fTower->AddCandidate(particle);
@@ -437 +441 @@
 
   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;
+
+  ecalEnergy = LogNormal(fTowerECalEnergy, ecalSigma);
+  ecalTime = (fTowerECalWeightTime < 1.0E-09 ) ? 0 : fTowerECalTime/fTowerECalWeightTime;
+
   hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerHCalEnergy);
 
-  ecalEnergy = LogNormal(fTowerECalEnergy, ecalSigma);
+//  hcalEnergy = gRandom->Gaus(fTowerHCalEnergy, hcalSigma);
+//  if(hcalEnergy < 0.0) hcalEnergy = 0.0;
+
   hcalEnergy = LogNormal(fTowerHCalEnergy, hcalSigma);
-
-  ecalTime = (fTowerECalTimeWeight < 1.0E-09 ) ? 0.0 : fTowerECalTime/fTowerECalTimeWeight;
-  hcalTime = (fTowerHCalTimeWeight < 1.0E-09 ) ? 0.0 : fTowerHCalTime/fTowerHCalTimeWeight;
-
+  hcalTime = (fTowerHCalWeightTime < 1.0E-09 ) ? 0 : fTowerHCalTime/fTowerHCalWeightTime;
+
+  
   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;
+  ecalEnergy = (ecalEnergy < fEcalEnergyMin || ecalEnergy < fEcalSigmaMin*ecalSigma) ? 0 : ecalEnergy;
+  hcalEnergy = (hcalEnergy < fHcalEnergyMin || hcalEnergy < fHcalSigmaMin*hcalSigma) ? 0 : hcalEnergy;
 
   energy = ecalEnergy + hcalEnergy;
-  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;
-  }
+  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]);
 
   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);
@@ -479 +491 @@
   fTower->Edges[3] = fTowerEdges[3];
 
-  if(energy > 0.0)
+  if( energy > 0.0 )
   {
     if(fTowerPhotonHits > 0 && fTowerTrackHits == 0)
@@ -485 +497 @@
       fPhotonOutputArray->Add(fTower);
     }
-
+   
     fTowerOutputArray->Add(fTower);
   }
@@ -499 +511 @@
 
   ecalEnergy -= fTrackECalEnergy;
+  if(ecalEnergy < fEcalEnergyMin || ecalEnergy < fEcalSigmaMin*fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy)) ecalEnergy = 0.0;
+
   hcalEnergy -= fTrackHCalEnergy;
-
-  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;
+  if(hcalEnergy < fHcalEnergyMin || hcalEnergy < fHcalSigmaMin*fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy)) hcalEnergy = 0.0;
 
   energy = ecalEnergy + hcalEnergy;
@@ -518 +527 @@
     tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalEnergy);
     tower->Eem = ecalEnergy;
-    tower->Ehad = 0.0;
+    tower->Ehad = 0;
 
     fEFlowPhotonOutputArray->Add(tower);
@@ -530 +539 @@
 
     tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalEnergy);
-    tower->Eem = 0.0;
+    tower->Eem = 0;
     tower->Ehad = hcalEnergy;
 
     fEFlowNeutralHadronOutputArray->Add(tower);
   }
+
+
+
+
 }
 
@@ -548 +561 @@
     a = TMath::Log(mean) - 0.5*b*b;
 
-    return TMath::Exp(a + b*gRandom->Gaus(0.0, 1.0));
+    return TMath::Exp(a + b*gRandom->Gaus(0, 1));
   }
   else