Changes in modules/DualReadoutCalorimeter.cc [9a7ea36:a1b19ea] in git

File:

• modules/DualReadoutCalorimeter.cc (modified) (22 diffs)

modules/DualReadoutCalorimeter.cc

@@ -24 +24 @@
   // 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.
+  // - 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.
 
@@ -63 +63 @@
   fItParticleInputArray(0), fItTrackInputArray(0)
 {
-
+  
   fECalResolutionFormula = new DelphesFormula;
   fHCalResolutionFormula = new DelphesFormula;
@@ -75 +75 @@
   fTowerTrackArray = new TObjArray;
   fItTowerTrackArray = fTowerTrackArray->MakeIterator();
-
+  
 }
 
@@ -82 +82 @@
 DualReadoutCalorimeter::~DualReadoutCalorimeter()
 {
-
+  
   if(fECalResolutionFormula) delete fECalResolutionFormula;
   if(fHCalResolutionFormula) delete fHCalResolutionFormula;
@@ -94 +94 @@
   if(fTowerTrackArray) delete fTowerTrackArray;
   if(fItTowerTrackArray) delete fItTowerTrackArray;
-
+ 
 }
 
@@ -172 +172 @@
   fECalEnergyMin = GetDouble("ECalEnergyMin", 0.0);
   fHCalEnergyMin = GetDouble("HCalEnergyMin", 0.0);
-  fEnergyMin = GetDouble("EnergyMin", 0.0);
 
   fECalEnergySignificanceMin = GetDouble("ECalEnergySignificanceMin", 0.0);
   fHCalEnergySignificanceMin = GetDouble("HCalEnergySignificanceMin", 0.0);
-  fEnergySignificanceMin = GetDouble("EnergySignificanceMin", 0.0);
 
   // switch on or off the dithering of the center of DualReadoutCalorimeter towers
@@ -247 +245 @@
   fItParticleInputArray->Reset();
   number = -1;
-  fTowerRmax=0.;
-
-  //cout<<"--------- new event ---------- "<<endl;
-
   while((particle = static_cast<Candidate*>(fItParticleInputArray->Next())))
   {
     const TLorentzVector &particlePosition = particle->Position;
     ++number;
-
-    // compute maximum radius (needed in FinalizeTower to assess whether barrel or endcap tower)
-    if (particlePosition.Perp() > fTowerRmax)
-      fTowerRmax=particlePosition.Perp();
 
     pdgCode = TMath::Abs(particle->PID);
@@ -387 +377 @@
       fHCalTrackEnergy = 0.0;
       fTrackEnergy = 0.0;
-
+      
       fECalTrackSigma = 0.0;
       fHCalTrackSigma = 0.0;
       fTrackSigma = 0.0;
-
+      
       fTowerTrackHits = 0;
       fTowerPhotonHits = 0;
-
-      fTowerTime = 0.0;
-      fTowerTimeWeight = 0.0;
 
       fECalTowerTrackArray->Clear();
       fHCalTowerTrackArray->Clear();
       fTowerTrackArray->Clear();
-
+    
     }
 
@@ -425 +412 @@
       }
 
+      
+      /*
+      if(fECalTrackFractions[number] > 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
+      {
+        fECalTrackEnergy += ecalEnergy;
+        ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());        
+        if(ecalSigma/momentum.E() < track->TrackResolution) energyGuess = ecalEnergy;        
+        else energyGuess = momentum.E();
+
+        fECalTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess;
+        fECalTowerTrackArray->Add(track);
+      }
+     
+      else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] > 1.0E-9)
+      {
+        fHCalTrackEnergy += hcalEnergy;
+        hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
+        if(hcalSigma/momentum.E() < track->TrackResolution) energyGuess = hcalEnergy;
+        else energyGuess = momentum.E();
+
+        fHCalTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess;
+        fHCalTowerTrackArray->Add(track);
+      }
+      
+      // 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)
+        // 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)
@@ -435 +453 @@
         else
           sigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
-
-        if(sigma/momentum.E() < track->TrackResolution)
-          energyGuess = ecalEnergy + hcalEnergy;
+          
+        if(sigma/momentum.E() < track->TrackResolution) 
+          energyGuess = ecalEnergy + hcalEnergy;     
         else
           energyGuess = momentum.E();
-
+             
         fTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess;
         fTowerTrackArray->Add(track);
+      
       }
       else
@@ -459 +478 @@
     position = particle->Position;
 
-
     // fill current tower
     ecalEnergy = momentum.E() * fECalTowerFractions[number];
@@ -467 +485 @@
     fHCalTowerEnergy += hcalEnergy;
 
-    // assume combined timing measurements in ECAL/HCAL sections
-    fTowerTime += (ecalEnergy + hcalEnergy) * position.T(); //sigma_t ~ 1/sqrt(E)
-    fTowerTimeWeight += ecalEnergy + hcalEnergy;
+    if(ecalEnergy > fTimingEnergyMin && fTower)
+    {
+      if (abs(particle->PID) != 11 || !fElectronsFromTrack)
+      {
+        fTower->ECalEnergyTimePairs.push_back(make_pair<Float_t, Float_t>(ecalEnergy, particle->Position.T()));
+      }
+    }
 
     fTower->AddCandidate(particle);
-    fTower->Position = position;
   }
 
@@ -485 +506 @@
 
   Candidate *track, *tower, *mother;
-  Double_t energy, pt, eta, phi, r, time;
+  Double_t energy, pt, eta, phi;
   Double_t ecalEnergy, hcalEnergy;
   Double_t ecalNeutralEnergy, hcalNeutralEnergy, neutralEnergy;
-
+  
   Double_t ecalSigma, hcalSigma, sigma;
   Double_t ecalNeutralSigma, hcalNeutralSigma, neutralSigma;
 
   Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor;
-
+  
   TLorentzVector momentum;
   TFractionMap::iterator itFractionMap;
@@ -501 +522 @@
   if(!fTower) return;
 
-  // if no hadronic energy, use ECAL resolution
+
+  //if (fHCalTowerEnergy < 30 && fECalTowerEnergy < 30) return;
+  //cout<<"----------- New tower ---------"<<endl;
+
+
+  // here we change behaviour w.r.t to standard calorimeter. Start with worse case scenario. If fHCalTowerEnergy > 0, assume total energy smeared by HCAL resolution.
+  // For example, if overlapping charged pions and photons take hadronic resolution as combined measurement
+
+  // if no hadronic fraction at all, then use ECAL resolution
+
+  //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);
-  }
-
-  // if hadronic fraction > 0, use HCAL resolution
+    //cout<<"em energy"<<energy<<", sigma: "<<sigma<<endl;
+  }
+
+  // if hadronic fraction > 0, use HCAL resolution 
   else
   {
     energy = fECalTowerEnergy + fHCalTowerEnergy;
     sigma  = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, energy);
+    //cout<<"had energy: "<<energy<<", sigma: "<<sigma<<endl;
   }
 
   energy = LogNormal(energy, sigma);
-
+  //cout<<energy<<","<<ecalEnergy<<","<<hcalEnergy<<endl;
+  
   if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0;
 
@@ -526 +562 @@
   hcalEnergy = LogNormal(fHCalTowerEnergy, hcalSigma);
 
-  time = (fTowerTimeWeight < 1.0E-09) ? 0.0 : fTowerTime / fTowerTimeWeight;
-
   ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy);
   hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy);
@@ -534 +568 @@
   if(hcalEnergy < fHCalEnergyMin || hcalEnergy < fHCalEnergySignificanceMin*hcalSigma) hcalEnergy = 0.0;
 
+  //cout<<"Measured energy: "<<energy<<endl;
+
   if(fSmearTowerCenter)
   {
@@ -547 +583 @@
   pt = energy / TMath::CosH(eta);
 
-  // check whether barrel or endcap tower
-
-  // endcap
-  if (TMath::Abs(fTower->Position.Pt() - fTowerRmax) > 1.e-06 && TMath::Abs(eta) > 0.){
-    r = fTower->Position.Z()/TMath::SinH(eta);
-  }
-  // barrel
-  else {
-    r = fTower->Position.Pt();
-  }
-
-  fTower->Position.SetPtEtaPhiE(r, eta, phi, time);
-  fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
-  fTower->L = fTower->Position.Vect().Mag();
+  // Time calculation for tower
+  fTower->NTimeHits = 0;
+  sumWeightedTime = 0.0;
+  sumWeight = 0.0;
+
+  for(size_t i = 0; i < fTower->ECalEnergyTimePairs.size(); ++i)
+  {
+    weight = TMath::Sqrt(fTower->ECalEnergyTimePairs[i].first);
+    sumWeightedTime += weight * fTower->ECalEnergyTimePairs[i].second;
+    sumWeight += weight;
+    fTower->NTimeHits++;
+  }
+
+  if(sumWeight > 0.0)
+  {
+    fTower->Position.SetPtEtaPhiE(1.0, eta, phi, sumWeightedTime/sumWeight);
+  }
+  else
+  {
+    fTower->Position.SetPtEtaPhiE(1.0, eta, phi, 999999.9);
+  }
 
   fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
   fTower->Eem = ecalEnergy;
   fTower->Ehad = hcalEnergy;
-  fTower->Etrk = fTrackEnergy;
+
   fTower->Edges[0] = fTowerEdges[0];
   fTower->Edges[1] = fTowerEdges[1];
@@ -581 +624 @@
 
   // 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);
+    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);
-
-    // if no hadronic energy, use ECAL resolution
-    if (fHCalTowerEnergy <= fHCalEnergyMin)
-    {
-      tower->Eem = neutralEnergy;
-      tower->Ehad = 0.0;
-      tower->PID = 22;
-      fEFlowPhotonOutputArray->Add(tower);
-    }
-    // if hadronic fraction > 0, use HCAL resolution
-    else
-    {
-      tower->Eem = 0;
-      tower->Ehad = neutralEnergy;
-      tower->PID = 130;
-      fEFlowNeutralHadronOutputArray->Add(tower);
-    }
+    tower->Eem = neutralEnergy;
+    tower->Ehad = 0.0;
+    tower->PID = 22;
+
+    fEFlowPhotonOutputArray->Add(tower);
+
 
     //clone tracks
@@ -614 +654 @@
     while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
     {
+      //cout<<"looping over tracks"<<endl;
       mother = track;
       track = static_cast<Candidate*>(track->Clone());
@@ -621 +662 @@
   }
 
-
+  
   // 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)
@@ -629 +670 @@
     weightCalo  = (sigma > 0.0) ? 1 / (sigma*sigma) : 0.0;
 
-    bestEnergyEstimate = (weightTrack*fTrackEnergy + weightCalo*energy) / (weightTrack + weightCalo);
+    bestEnergyEstimate = (weightTrack*fTrackEnergy + weightCalo*energy) / (weightTrack + weightCalo); 
     rescaleFactor = bestEnergyEstimate/fTrackEnergy;
 
@@ -635 +676 @@
     fItTowerTrackArray->Reset();
     while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
-    {
+    {  
       mother = track;
-      track = static_cast<Candidate *>(track->Clone());
+      track = static_cast<Candidate*>(track->Clone());
       track->AddCandidate(mother);
-      track->Momentum.SetPtEtaPhiM(track->Momentum.Pt()*rescaleFactor, track->Momentum.Eta(), track->Momentum.Phi(), track->Momentum.M());
+
+      track->Momentum *= rescaleFactor;
+
       fEFlowTrackOutputArray->Add(track);
     }
   }
-
-
+  
+
+  /*
+  // 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)
+  {
+    // 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;
+    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 DualReadoutCalorimeter 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
+    tower = static_cast<Candidate*>(fTower->Clone());
+    pt =  hcalNeutralEnergy / TMath::CosH(eta);
+    
+    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalNeutralEnergy);
+    tower->Ehad = hcalNeutralEnergy;
+    tower->Eem = 0.0;
+    
+    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 DualReadoutCalorimeter 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);
+    }
+  }
+
+  */
 }