Changeset f6b9fec in git for modules/Calorimeter.cc

Timestamp:

Apr 16, 2014, 3:29:31 PM (11 years ago)

Author:

mselvaggi <mselvaggi@…>

Branches:

ImprovedOutputFile, Timing, dual_readout, llp, master

Children:

3c40083

Parents:

c526439

Message:

new calo module

File:

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

modules/Calorimeter.cc

@@ -41 +41 @@
 
 Calorimeter::Calorimeter() :
-  fECalResolutionFormula(0), fHCalResolutionFormula(0),
+  fResolutionFormula(0),
   fItParticleInputArray(0), fItTrackInputArray(0),
   fTowerTrackArray(0), fItTowerTrackArray(0)
 {
-  fECalResolutionFormula = new DelphesFormula;
-  fHCalResolutionFormula = new DelphesFormula;
-
+  fResolutionFormula = new DelphesFormula;
+ 
   fTowerTrackArray = new TObjArray;
   fItTowerTrackArray = fTowerTrackArray->MakeIterator();
@@ -56 +55 @@
 Calorimeter::~Calorimeter()
 {
-  if(fECalResolutionFormula) delete fECalResolutionFormula;
-  if(fHCalResolutionFormula) delete fHCalResolutionFormula;
-
+  if(fResolutionFormula) delete fResolutionFormula;
+ 
   if(fTowerTrackArray) delete fTowerTrackArray;
   if(fItTowerTrackArray) delete fItTowerTrackArray;
@@ -68 +66 @@
 {
   ExRootConfParam param, paramEtaBins, paramPhiBins, paramFractions;
-  Long_t i, j, k, size, sizeEtaBins, sizePhiBins, sizeFractions;
-  Double_t ecalFraction, hcalFraction;
+  Long_t i, j, k, size, sizeEtaBins, sizePhiBins;
+  Double_t fraction;
   TBinMap::iterator itEtaBin;
   set< Double_t >::iterator itPhiBin;
@@ -116 +114 @@
   // set default energy fractions values
   fFractionMap.clear();
-  fFractionMap[0] = make_pair(0.0, 1.0);
+  fFractionMap[0] = 1.0;
 
   for(i = 0; i < size/2; ++i)
   {
     paramFractions = param[i*2 + 1];
-    sizeFractions = paramFractions.GetSize();
-
-    ecalFraction = paramFractions[0].GetDouble();
-    hcalFraction = paramFractions[1].GetDouble();
-
-    fFractionMap[param[i*2].GetInt()] = make_pair(ecalFraction, hcalFraction);
+    fraction = paramFractions[0].GetDouble();
+    fFractionMap[param[i*2].GetInt()] = fraction;
   }
 /*
@@ -136 +130 @@
 */
   // read resolution formulas
-  fECalResolutionFormula->Compile(GetString("ECalResolutionFormula", "0"));
-  fHCalResolutionFormula->Compile(GetString("HCalResolutionFormula", "0"));
-
+  fResolutionFormula->Compile(GetString("ResolutionFormula", "0"));
+ 
   // import array with output from other modules
   fParticleInputArray = ImportArray(GetString("ParticleInputArray", "ParticlePropagator/particles"));
@@ -148 +141 @@
   // create output arrays
   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"));
-
-
+  fEFlowTowerOutputArray = ExportArray(GetString("EFlowTowerOutputArray", "eflowTowers"));
+ 
 }
 
@@ -179 +167 @@
   Int_t number;
   Long64_t towerHit, towerEtaPhi, hitEtaPhi;
-  Double_t ecalFraction, hcalFraction;
-  Double_t ecalEnergy, hcalEnergy;
+  Double_t fraction;
+  Double_t energy;
   Int_t pdgCode;
 
@@ -193 +181 @@
   DelphesFactory *factory = GetFactory();
   fTowerHits.clear();
-  fTowerECalFractions.clear();
-  fTowerHCalFractions.clear();
-  fTrackECalFractions.clear();
-  fTrackHCalFractions.clear();
-
+  fTowerFractions.clear();
+  fTrackFractions.clear();
+ 
   // loop over all particles
   fItParticleInputArray->Reset();
@@ -214 +200 @@
     }
 
-    ecalFraction = itFractionMap->second.first;
-    hcalFraction = itFractionMap->second.second;
-
-    fTowerECalFractions.push_back(ecalFraction);
-    fTowerHCalFractions.push_back(hcalFraction);
-
-    if(ecalFraction < 1.0E-9 && hcalFraction < 1.0E-9) continue;
+    fraction = itFractionMap->second;
+    fTowerFractions.push_back(fraction);
+   
+    if(fraction < 1.0E-9) continue;
 
     // find eta bin [1, fEtaBins.size - 1]
@@ -260 +243 @@
     }
 
-    ecalFraction = itFractionMap->second.first;
-    hcalFraction = itFractionMap->second.second;
-
-    fTrackECalFractions.push_back(ecalFraction);
-    fTrackHCalFractions.push_back(hcalFraction);
-
+    fraction = itFractionMap->second;
+  
+    fTrackFractions.push_back(fraction);
+  
     // find eta bin [1, fEtaBins.size - 1]
     itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), trackPosition.Eta());
@@ -327 +308 @@
       fTowerEdges[3] = (*phiBins)[phiBin];
 
-      fTowerECalEnergy = 0.0;
-      fTowerHCalEnergy = 0.0;
-
-      fTrackECalEnergy = 0.0;
-      fTrackHCalEnergy = 0.0;
-
-      fTowerECalTime = 0.0;
-      fTowerHCalTime = 0.0;
-
-      fTrackECalTime = 0.0;
-      fTrackHCalTime = 0.0;
-
-      fTowerECalWeightTime = 0.0; 
-      fTowerHCalWeightTime = 0.0;
+      fTowerEnergy = 0.0;
+      fTrackEnergy = 0.0;
+      
+      fTowerTime = 0.0;
+      fTrackTime = 0.0;
      
+      fTowerWeightTime = 0.0; 
+      
       fTowerTrackHits = 0;
       fTowerPhotonHits = 0;
@@ -356 +330 @@
       momentum = track->Momentum;
       position = track->Position;
-
+ 
+      energy = momentum.E() * fTrackFractions[number];
       
-      ecalEnergy = momentum.E() * fTrackECalFractions[number];
-      hcalEnergy = momentum.E() * fTrackHCalFractions[number];
-
-      fTrackECalEnergy += ecalEnergy;
-      fTrackHCalEnergy += hcalEnergy;
+      fTrackEnergy += energy;
       
-      fTrackECalTime += TMath::Sqrt(ecalEnergy)*position.T();
-      fTrackHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
-       
-      fTrackECalWeightTime += TMath::Sqrt(ecalEnergy);
-      fTrackHCalWeightTime += TMath::Sqrt(hcalEnergy);
-
+      fTrackTime += TMath::Sqrt(energy)*position.T();
+      fTrackWeightTime += TMath::Sqrt(energy);
+    
       fTowerTrackArray->Add(track);
 
@@ -383 +351 @@
 
     // fill current tower
-    ecalEnergy = momentum.E() * fTowerECalFractions[number];
-    hcalEnergy = momentum.E() * fTowerHCalFractions[number];
-
-    fTowerECalEnergy += ecalEnergy;
-    fTowerHCalEnergy += hcalEnergy;
-
-    fTowerECalTime += TMath::Sqrt(ecalEnergy)*position.T();
-    fTowerHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
-
-    fTowerECalWeightTime += TMath::Sqrt(ecalEnergy);
-    fTowerHCalWeightTime += TMath::Sqrt(hcalEnergy);
+    energy = momentum.E() * fTowerFractions[number];
+   
+    fTowerEnergy += energy;
     
-
+    fTowerTime += TMath::Sqrt(energy)*position.T();
+    fTowerWeightTime += TMath::Sqrt(energy);
+    
     fTower->AddCandidate(particle);
   }
@@ -407 +369 @@
 void Calorimeter::FinalizeTower()
 {
-  Candidate *track, *tower;
+  Candidate *tower;
   Double_t energy, pt, eta, phi;
-  Double_t ecalEnergy, hcalEnergy;
-  Double_t ecalSigma, hcalSigma;
-  Double_t ecalTime, hcalTime, time;
+  Double_t sigma;
+  Double_t time;
 
   if(!fTower) return;
 
-  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);
-
-//  hcalEnergy = gRandom->Gaus(fTowerHCalEnergy, hcalSigma);
-//  if(hcalEnergy < 0.0) hcalEnergy = 0.0;
-
-  hcalEnergy = LogNormal(fTowerHCalEnergy, hcalSigma);
-  hcalTime = (fTowerHCalWeightTime < 1.0E-09 ) ? 0 : fTowerHCalTime/fTowerHCalWeightTime;
-
-  energy = ecalEnergy + hcalEnergy;
-  time = (TMath::Sqrt(ecalEnergy)*ecalTime + TMath::Sqrt(hcalEnergy)*hcalTime)/(TMath::Sqrt(ecalEnergy) + TMath::Sqrt(hcalEnergy)); 
-
-//  eta = fTowerEta;
-//  phi = fTowerPhi;
+  sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerEnergy);
+
+//  energy = gRandom->Gaus(fTowerEnergy, sigma);
+//  if(energy < 0.0) energy = 0.0;
+
+  energy = LogNormal(fTowerEnergy, sigma);
+  time = (fTowerWeightTime < 1.0E-09 ) ? 0 : fTowerTime/fTowerWeightTime;
 
   eta = gRandom->Uniform(fTowerEdges[0], fTowerEdges[1]);
@@ -445 +392 @@
   fTower->Position.SetPtEtaPhiE(1.0, eta, phi, time);
   fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
-  fTower->Eem = ecalEnergy;
-  fTower->Ehad = hcalEnergy;
-
+ 
   fTower->Edges[0] = fTowerEdges[0];
   fTower->Edges[1] = fTowerEdges[1];
@@ -455 +400 @@
 
   // fill calorimeter towers
+  if(energy > 0.0) fTowerOutputArray->Add(fTower);
+
+  
+  // fill energy flow candidates
+  energy -= fTrackEnergy;
+  if(energy < 0.0) energy = 0.0;
+   
+  // save energy excess as an energy flow tower
   if(energy > 0.0)
-  {
-    if(fTowerPhotonHits > 0 && fTowerTrackHits == 0)
-    {
-      fPhotonOutputArray->Add(fTower);
-    }
-   
-    fTowerOutputArray->Add(fTower);
-  }
-
-  // fill energy flow candidates
-
-  // save all the tracks as energy flow tracks
-  fItTowerTrackArray->Reset();
-  while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
-  {
-    fEFlowTrackOutputArray->Add(track);
-  }
-
-  ecalEnergy -= fTrackECalEnergy;
-  if(ecalEnergy < 0.0) ecalEnergy = 0.0;
-
-  hcalEnergy -= fTrackHCalEnergy;
-  if(hcalEnergy < 0.0) hcalEnergy = 0.0;
-
-  energy = ecalEnergy + hcalEnergy;
-
-  
-  // save ECAL and/or HCAL energy excess as an energy flow tower
-  if(ecalEnergy > 0.0)
   {
     // create new photon tower
     tower = static_cast<Candidate*>(fTower->Clone());
-
-    pt = ecalEnergy / TMath::CosH(eta);
-
-    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalEnergy);
-    tower->Eem = ecalEnergy;
-    tower->Ehad = 0;
-
-    fEFlowPhotonOutputArray->Add(tower);
-  }
-
-  if(hcalEnergy > 0.0)
-  {
-    // create new neutral hadron tower
-    tower = static_cast<Candidate*>(fTower->Clone());
-
-    pt = hcalEnergy / TMath::CosH(eta);
-
-    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalEnergy);
-    tower->Eem = 0;
-    tower->Ehad = hcalEnergy;
-
-    fEFlowNeutralHadronOutputArray->Add(tower);
-  }
-
-
-
+    pt = energy / TMath::CosH(eta);
+
+    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
+    fEFlowTowerOutputArray->Add(tower);
+  }
 
 }