Changeset a1b19ea in git

Timestamp:

Jan 11, 2019, 4:39:56 PM (6 years ago)

Author:

Michele Selvaggi <michele.selvaggi@…>

Branches:

ImprovedOutputFile, Timing, dual_readout, master

Children:

e39abb4

Parents:

ef8a06d

Message:

first attempt at energy flow with Dual Readout

Files:

• cards/delphes_card_IDEA.tcl (modified) (2 diffs)
• modules/DualReadoutCalorimeter.cc (modified) (16 diffs)
• modules/DualReadoutCalorimeter.h (modified) (4 diffs)

cards/delphes_card_IDEA.tcl

@@ -211 +211 @@
   set ECalEnergyMin 0.5
   set HCalEnergyMin 1.0
-  set EnergyMin 0.0
+  set EnergyMin 0.5
 
   set ECalEnergySignificanceMin 1.0
@@ -279 +279 @@
 
   # set HCalResolutionFormula {resolution formula as a function of eta and energy}
-  set HCalResolutionFormula {                  (abs(eta) <= 3.0) * sqrt(energy^2*0.050^2 + energy*1.50^2) +
+  set HCalResolutionFormula {                  (abs(eta) <= 3.0) * sqrt(energy^2*0.050^2 + energy*0.30^2) +
                              (abs(eta) > 3.0 && abs(eta) <= 5.0) * sqrt(energy^2*0.130^2 + energy*2.70^2)}
 }

modules/DualReadoutCalorimeter.cc

@@ -20 +20 @@
 /** \class DualReadoutCalorimeter
  *
- *  Fills DualReadoutCalorimeter towers, performs DualReadoutCalorimeter resolution smearing,
- *  and creates energy flow objects (tracks, photons, and neutral hadrons).
+
+  // ============  TODO  =========================================================
+  // This implementation of dual calorimetry relies on several approximations:
+  // - If hadronic energy is found in the tower the energy resolution then the full tower enrgy is smeared according to hadronic resolution (pessimistic for (e,n) or (pi+,gamma))
+  // - While e+ vs pi+ (or gamma vs n) separation is in principle possible for single particles (using C/S, PMT timing, lateral shower profile) it is not obvious it can be done overlapping particles. 
+  //   Now we assume that regarless of the number of particle hits per tower we can always distinguish e+ vs pi+, which is probably not true in the case (e+,n) vs (pi+,gamma) without longitudinal segmentation.
+
  *
- *  \author P. Demin - UCL, Louvain-la-Neuve
+ *  \author M. Selvaggi - CERN
  *
  */
@@ -67 +72 @@
   fHCalTowerTrackArray = new TObjArray;
   fItHCalTowerTrackArray = fHCalTowerTrackArray->MakeIterator();
+
+  fTowerTrackArray = new TObjArray;
+  fItTowerTrackArray = fTowerTrackArray->MakeIterator();
   
 }
@@ -83 +91 @@
   if(fHCalTowerTrackArray) delete fHCalTowerTrackArray;
   if(fItHCalTowerTrackArray) delete fItHCalTowerTrackArray;
+
+  if(fTowerTrackArray) delete fTowerTrackArray;
+  if(fItTowerTrackArray) delete fItTowerTrackArray;
  
 }
@@ -212 +223 @@
   Double_t ecalFraction, hcalFraction;
   Double_t ecalEnergy, hcalEnergy;
-  Double_t ecalSigma, hcalSigma;
-  Double_t energyGuess;
+  Double_t ecalSigma, hcalSigma, sigma;
+  Double_t energyGuess, energy;
   Int_t pdgCode;
 
@@ -365 +376 @@
       fECalTrackEnergy = 0.0;
       fHCalTrackEnergy = 0.0;
+      fTrackEnergy = 0.0;
       
       fECalTrackSigma = 0.0;
       fHCalTrackSigma = 0.0;
+      fTrackSigma = 0.0;
       
       fTowerTrackHits = 0;
@@ -374 +387 @@
       fECalTowerTrackArray->Clear();
       fHCalTowerTrackArray->Clear();
+      fTowerTrackArray->Clear();
     
     }
@@ -388 +402 @@
       ecalEnergy = momentum.E() * fECalTrackFractions[number];
       hcalEnergy = momentum.E() * fHCalTrackFractions[number];
+      energy = ecalEnergy + hcalEnergy;
 
       if(ecalEnergy > fTimingEnergyMin && fTower)
@@ -397 +412 @@
       }
 
+      
+      /*
       if(fECalTrackFractions[number] > 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
       {
@@ -419 +436 @@
       }
       
+      // muons
       else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
       {
         fEFlowTrackOutputArray->Add(track);
       }
+      */
+      
+      // in Dual Readout we do not care if tracks are ECAL of HCAL
+      if(fECalTrackFractions[number] > 1.0E-9 || fHCalTrackFractions[number] > 1.0E-9)
+      { 
+        fTrackEnergy += energy;
+        // this sigma will be used to determine whether neutral excess is significant. We choose the resolution according to bthe higest deposited fraction (in practice had for charged hadrons and em for electrons)      
+        sigma = 0.0;
+        if(fHCalTrackFractions[number] > 0)
+          sigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
+        else
+          sigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
+          
+        if(sigma/momentum.E() < track->TrackResolution) 
+          energyGuess = ecalEnergy + hcalEnergy;     
+        else
+          energyGuess = momentum.E();
+             
+        fTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess;
+        fTowerTrackArray->Add(track);
+      
+      }
+      else
+      {
+        fEFlowTrackOutputArray->Add(track);
+      }
 
       continue;
@@ -461 +505 @@
 {
 
-
-
-
   Candidate *track, *tower, *mother;
   Double_t energy, pt, eta, phi;
   Double_t ecalEnergy, hcalEnergy;
-  Double_t ecalNeutralEnergy, hcalNeutralEnergy;
+  Double_t ecalNeutralEnergy, hcalNeutralEnergy, neutralEnergy;
   
   Double_t ecalSigma, hcalSigma, sigma;
-  Double_t ecalNeutralSigma, hcalNeutralSigma;
+  Double_t ecalNeutralSigma, hcalNeutralSigma, neutralSigma;
 
   Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor;
@@ -491 +532 @@
   // if no hadronic fraction at all, then use ECAL resolution
 
-  //cout<<"fECalTowerEnergy: "<<fECalTowerEnergy<<", fHCalTowerEnergy: "<<fHCalTowerEnergy<<endl;
-
+  //cout<<"fECalTowerEnergy: "<<fECalTowerEnergy<<", fHCalTowerEnergy: "<<fHCalTowerEnergy<<", Eta: "<<fTowerEta<<endl;
+
+  // if no hadronic energy, use ECAL resolution 
   if (fHCalTowerEnergy <= fHCalEnergyMin)
   {
     energy = fECalTowerEnergy;
     sigma  = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, energy);
-    //cout<<"em: "<<energy<<","<<sigma<<endl;
+    //cout<<"em energy"<<energy<<", sigma: "<<sigma<<endl;
   }
 
@@ -505 +547 @@
     energy = fECalTowerEnergy + fHCalTowerEnergy;
     sigma  = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, energy);
-    //cout<<"had: "<<energy<<","<<sigma<<endl;
+    //cout<<"had energy: "<<energy<<", sigma: "<<sigma<<endl;
   }
 
@@ -525 +567 @@
   if(ecalEnergy < fECalEnergyMin || ecalEnergy < fECalEnergySignificanceMin*ecalSigma) ecalEnergy = 0.0;
   if(hcalEnergy < fHCalEnergyMin || hcalEnergy < fHCalEnergySignificanceMin*hcalSigma) hcalEnergy = 0.0;
-  //cout<<energy<<","<<ecalEnergy<<","<<hcalEnergy<<endl;
-
+
+  //cout<<"Measured energy: "<<energy<<endl;
 
   if(fSmearTowerCenter)
@@ -563 +605 @@
   }
 
-
   fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
   fTower->Eem = ecalEnergy;
@@ -579 +620 @@
       fPhotonOutputArray->Add(fTower);
     }
-
     fTowerOutputArray->Add(fTower);
   }
-  
+
+  // fill energy flow candidates
+  
+  fTrackSigma = TMath::Sqrt(fTrackSigma);
+  neutralEnergy = max( (energy - fTrackEnergy) , 0.0);
+  neutralSigma = neutralEnergy / TMath::Sqrt(fTrackSigma*fTrackSigma + sigma*sigma);
+
+  //cout<<"trackEnergy: "<<fTrackEnergy<<", trackSigma: "<<fTrackSigma<<", Ntracks: "<<fTowerTrackArray->GetEntries()<<endl;
+
+  //cout<<"neutralEnergy: "<<neutralEnergy<<", neutralSigma: "<<neutralSigma<<endl;
+
+   // For now, if neutral excess is significant, simply create neutral EflowPhoton tower and clone each track into eflowtrack !!! -> Creating only photons !! EFlowNeutralHadron collection will be empy!!! TO BE FIXED
+  if(neutralEnergy > fEnergyMin && neutralSigma > fEnergySignificanceMin)
+  {
+    
+    //cout<<"significant neutral excess found:"<<endl;
+    // create new photon tower
+    tower = static_cast<Candidate*>(fTower->Clone());
+    pt =  neutralEnergy / TMath::CosH(eta);
+    //cout<<"Creating tower with Pt, Eta, Phi, Energy: "<<pt<<","<<eta<<","<<phi<<","<<neutralEnergy<<endl;
+    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, neutralEnergy);
+    tower->Eem = neutralEnergy;
+    tower->Ehad = 0.0;
+    tower->PID = 22;
+
+    fEFlowPhotonOutputArray->Add(tower);
+
+
+    //clone tracks
+    fItTowerTrackArray->Reset();
+    while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
+    {
+      //cout<<"looping over tracks"<<endl;
+      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 DualReadoutCalorimeter and tracking
+  else if(fTrackEnergy > 0.0)
+  {
+    //cout<<"no significant neutral excess found:"<<endl;
+    weightTrack = (fTrackSigma > 0.0) ? 1 / (fTrackSigma*fTrackSigma) : 0.0;
+    weightCalo  = (sigma > 0.0) ? 1 / (sigma*sigma) : 0.0;
+
+    bestEnergyEstimate = (weightTrack*fTrackEnergy + weightCalo*energy) / (weightTrack + weightCalo); 
+    rescaleFactor = bestEnergyEstimate/fTrackEnergy;
+
+    //rescale tracks
+    fItTowerTrackArray->Reset();
+    while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
+    {  
+      mother = track;
+      track = static_cast<Candidate*>(track->Clone());
+      track->AddCandidate(mother);
+
+      track->Momentum *= rescaleFactor;
+
+      fEFlowTrackOutputArray->Add(track);
+    }
+  }
+  
+
+  /*
   // fill energy flow candidates
   fECalTrackSigma = TMath::Sqrt(fECalTrackSigma);
@@ -693 +799 @@
     }
   }
-  
-  
+
+  */
 }
 

modules/DualReadoutCalorimeter.h

@@ -25 +25 @@
  *  and creates energy flow objects (tracks, photons, and neutral hadrons).
  *
- *  \author P. Demin - UCL, Louvain-la-Neuve
+ *  \author M. Selvaggi - CERN
  *
  */
@@ -59 +59 @@
   Double_t fECalTowerEnergy, fHCalTowerEnergy;
   Double_t fECalTrackEnergy, fHCalTrackEnergy;
+  Double_t fTrackEnergy;
 
   Double_t fTimingEnergyMin;
@@ -75 +76 @@
   Double_t fECalTrackSigma;
   Double_t fHCalTrackSigma;
+  Double_t fTrackSigma;
 
   Bool_t fSmearTowerCenter;
@@ -114 +116 @@
   TIterator *fItHCalTowerTrackArray; //!
 
+  TObjArray *fTowerTrackArray; //!
+  TIterator *fItTowerTrackArray; //!
+
   void FinalizeTower();
   Double_t LogNormal(Double_t mean, Double_t sigma);