Changes in modules/Calorimeter.cc [b6e6d36:f8299bc] in git

File:

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

modules/Calorimeter.cc

@@ -58 +58 @@
   fItParticleInputArray(0), fItTrackInputArray(0)
 {
-  
+  Int_t i;
+
   fECalResolutionFormula = new DelphesFormula;
   fHCalResolutionFormula = new DelphesFormula;
 
-  fECalTowerTrackArray = new TObjArray;
-  fItECalTowerTrackArray = fECalTowerTrackArray->MakeIterator();
-
-  fHCalTowerTrackArray = new TObjArray;
-  fItHCalTowerTrackArray = fHCalTowerTrackArray->MakeIterator();
-  
+  for(i = 0; i < 2; ++i)
+  {
+    fECalTowerTrackArray[i] = new TObjArray;
+    fItECalTowerTrackArray[i] = fECalTowerTrackArray[i]->MakeIterator();
+
+    fHCalTowerTrackArray[i] = new TObjArray;
+    fItHCalTowerTrackArray[i] = fHCalTowerTrackArray[i]->MakeIterator();
+  }
 }
 
@@ -74 +77 @@
 Calorimeter::~Calorimeter()
 {
-  
+  Int_t i;
+
   if(fECalResolutionFormula) delete fECalResolutionFormula;
   if(fHCalResolutionFormula) delete fHCalResolutionFormula;
 
-  if(fECalTowerTrackArray) delete fECalTowerTrackArray;
-  if(fItECalTowerTrackArray) delete fItECalTowerTrackArray;
-
-  if(fHCalTowerTrackArray) delete fHCalTowerTrackArray;
-  if(fItHCalTowerTrackArray) delete fItHCalTowerTrackArray;
- 
+  for(i = 0; i < 2; ++i)
+  {
+    if(fECalTowerTrackArray[i]) delete fECalTowerTrackArray[i];
+    if(fItECalTowerTrackArray[i]) delete fItECalTowerTrackArray[i];
+
+    if(fHCalTowerTrackArray[i]) delete fHCalTowerTrackArray[i];
+    if(fItHCalTowerTrackArray[i]) delete fItHCalTowerTrackArray[i];
+  }
 }
 
@@ -212 +218 @@
   Double_t ecalFraction, hcalFraction;
   Double_t ecalEnergy, hcalEnergy;
+  Double_t ecalSigma, hcalSigma;
   Int_t pdgCode;
 
@@ -361 +368 @@
       fHCalTowerEnergy = 0.0;
 
-      fECalTrackEnergy = 0.0;
-      fHCalTrackEnergy = 0.0;
-      
-      fECalTrackSigma = 0.0;
-      fHCalTrackSigma = 0.0;
-      
+      fECalTrackEnergy[0] = 0.0;
+      fECalTrackEnergy[1] = 0.0;
+
+      fHCalTrackEnergy[0] = 0.0;
+      fHCalTrackEnergy[1] = 0.0;
+
       fTowerTrackHits = 0;
       fTowerPhotonHits = 0;
 
-      fECalTowerTrackArray->Clear();
-      fHCalTowerTrackArray->Clear();
-    
+      fECalTowerTrackArray[0]->Clear();
+      fECalTowerTrackArray[1]->Clear();
+
+      fHCalTowerTrackArray[0]->Clear();
+      fHCalTowerTrackArray[1]->Clear();
     }
 
@@ -397 +406 @@
       if(fECalTrackFractions[number] > 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
       {
-        fECalTrackEnergy += ecalEnergy;
-        fECalTrackSigma += (track->TrackResolution)*momentum.E()*(track->TrackResolution)*momentum.E();
-        fECalTowerTrackArray->Add(track);
+        ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
+        if(ecalSigma/momentum.E() < track->TrackResolution)
+        {
+          fECalTrackEnergy[0] += ecalEnergy;
+          fECalTowerTrackArray[0]->Add(track);
+        }
+        else
+        {
+          fECalTrackEnergy[1] += ecalEnergy;
+          fECalTowerTrackArray[1]->Add(track);
+        }
       }
-     
       else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] > 1.0E-9)
       {
-        fHCalTrackEnergy += hcalEnergy;
-        fHCalTrackSigma += (track->TrackResolution)*momentum.E()*(track->TrackResolution)*momentum.E();
-        fHCalTowerTrackArray->Add(track);
+        hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
+        if(hcalSigma/momentum.E() < track->TrackResolution)
+        {
+          fHCalTrackEnergy[0] += hcalEnergy;
+          fHCalTowerTrackArray[0]->Add(track);
+        }
+        else
+        {
+          fHCalTrackEnergy[1] += hcalEnergy;
+          fHCalTowerTrackArray[1]->Add(track);
+        }
       }
-      
       else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
       {
@@ -453 +476 @@
   Double_t energy, pt, eta, phi;
   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;
@@ -466 +484 @@
 
   if(!fTower) return;
+
 
   ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fECalTowerEnergy);
@@ -535 +554 @@
     fTowerOutputArray->Add(fTower);
   }
-  
+
   // fill energy flow candidates
-  fECalTrackSigma = TMath::Sqrt(fECalTrackSigma);
-  fHCalTrackSigma = TMath::Sqrt(fHCalTrackSigma);
-
-  //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
-  if(ecalNeutralEnergy > fECalEnergyMin && ecalNeutralSigma > fECalEnergySignificanceMin)
+
+  ecalEnergy -= fECalTrackEnergy[1];
+  hcalEnergy -= fHCalTrackEnergy[1];
+
+  fItECalTowerTrackArray[0]->Reset();
+  while((track = static_cast<Candidate*>(fItECalTowerTrackArray[0]->Next())))
+  {
+    mother = track;
+    track = static_cast<Candidate*>(track->Clone());
+    track->AddCandidate(mother);
+
+    track->Momentum *= ecalEnergy/fECalTrackEnergy[0];
+
+    fEFlowTrackOutputArray->Add(track);
+  }
+
+  fItECalTowerTrackArray[1]->Reset();
+  while((track = static_cast<Candidate*>(fItECalTowerTrackArray[1]->Next())))
+  {
+    mother = track;
+    track = static_cast<Candidate*>(track->Clone());
+    track->AddCandidate(mother);
+
+    fEFlowTrackOutputArray->Add(track);
+  }
+
+  fItHCalTowerTrackArray[0]->Reset();
+  while((track = static_cast<Candidate*>(fItHCalTowerTrackArray[0]->Next())))
+  {
+    mother = track;
+    track = static_cast<Candidate*>(track->Clone());
+    track->AddCandidate(mother);
+
+    track->Momentum *= hcalEnergy/fHCalTrackEnergy[0];
+
+    fEFlowTrackOutputArray->Add(track);
+  }
+
+  fItHCalTowerTrackArray[1]->Reset();
+  while((track = static_cast<Candidate*>(fItHCalTowerTrackArray[1]->Next())))
+  {
+    mother = track;
+    track = static_cast<Candidate*>(track->Clone());
+    track->AddCandidate(mother);
+
+    fEFlowTrackOutputArray->Add(track);
+  }
+
+  if(fECalTowerTrackArray[0]->GetEntriesFast() > 0) ecalEnergy = 0.0;
+  if(fHCalTowerTrackArray[0]->GetEntriesFast() > 0) 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;
+
+  if(ecalEnergy > 0.0)
   {
     // create new photon tower
     tower = static_cast<Candidate*>(fTower->Clone());
-    pt =  ecalNeutralEnergy / TMath::CosH(eta);
-    
-    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalNeutralEnergy);
-    tower->Eem = ecalNeutralEnergy;
+
+    pt = ecalEnergy / TMath::CosH(eta);
+
+    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalEnergy);
+    tower->Eem = ecalEnergy;
     tower->Ehad = 0.0;
     tower->PID = 22;
-    
+
     fEFlowPhotonOutputArray->Add(tower);
-   
-    //clone tracks
-    fItECalTowerTrackArray->Reset();
-    while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next())))
-    {
-      mother = track;
-      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;
-
-    //rescale tracks
-    fItECalTowerTrackArray->Reset();
-    while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next())))
-    {  
-      mother = track;
-      track = static_cast<Candidate*>(track->Clone());
-      track->AddCandidate(mother);
-
-      track->Momentum *= rescaleFactor;
-
-      fEFlowTrackOutputArray->Add(track);
-    }
-  }
-
-
-  // if hcal neutral excess is significant, simply create neutral EflowNeutralHadron tower and clone each track into eflowtrack
-  if(hcalNeutralEnergy > fHCalEnergyMin && hcalNeutralSigma > fHCalEnergySignificanceMin)
-  {
-    // create new photon tower
+  }
+  if(hcalEnergy > 0.0)
+  {
+    // create new neutral hadron tower
     tower = static_cast<Candidate*>(fTower->Clone());
-    pt =  hcalNeutralEnergy / TMath::CosH(eta);
-    
-    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalNeutralEnergy);
-    tower->Ehad = hcalNeutralEnergy;
+
+    pt = hcalEnergy / TMath::CosH(eta);
+
+    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalEnergy);
     tower->Eem = 0.0;
-    
+    tower->Ehad = hcalEnergy;
+
     fEFlowNeutralHadronOutputArray->Add(tower);
-   
-    //clone tracks
-    fItHCalTowerTrackArray->Reset();
-    while((track = static_cast<Candidate*>(fItHCalTowerTrackArray->Next())))
-    {
-      mother = track;
-      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); 
-    rescaleFactor = bestEnergyEstimate / fHCalTrackEnergy;
-
-    //rescale tracks
-    fItHCalTowerTrackArray->Reset();
-    while((track = static_cast<Candidate*>(fItHCalTowerTrackArray->Next())))
-    {  
-      mother = track;
-      track = static_cast<Candidate*>(track->Clone());
-      track->AddCandidate(mother);
-
-      track->Momentum *= rescaleFactor;
-
-      fEFlowTrackOutputArray->Add(track);
-    }
-  }
-  
-  
+  }
 }