Changeset 1369 in svn for trunk/modules/Calorimeter.cc

Timestamp:

Apr 16, 2014, 5:17:35 PM (10 years ago)

Author:

Michele Selvaggi

Message:

added simplecalo

File:

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

trunk/modules/Calorimeter.cc

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