Changes in modules/Calorimeter.cc [3db5282:e33e6db] in git

File:

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

modules/Calorimeter.cc

@@ -150 +150 @@
 */
 
-  // read min E value for timing measurement in ECAL
-  fTimingEMin = GetDouble("TimingEMin",4.);
-  // For timing
-  // So far this flag needs to be false
-  // Curved extrapolation not supported
-  fElectronsFromTrack = false;
-
-  
   // read min E value for towers to be saved
   fECalEnergyMin = GetDouble("ECalEnergyMin", 0.0);
@@ -364 +356 @@
       fTrackHCalEnergy = 0.0;
 
+      fTowerECalTime = 0.0;
+      fTowerHCalTime = 0.0;
+
+      fTrackECalTime = 0.0;
+      fTrackHCalTime = 0.0;
+
+      fTowerECalTimeWeight = 0.0;
+      fTowerHCalTimeWeight = 0.0;
+
       fTowerTrackHits = 0;
       fTowerPhotonHits = 0;
@@ -379 +380 @@
       position = track->Position;
 
+
       ecalEnergy = momentum.E() * fTrackECalFractions[number];
       hcalEnergy = momentum.E() * fTrackHCalFractions[number];
@@ -385 +387 @@
       fTrackHCalEnergy += hcalEnergy;
 
-      bool dbg_scz = false;
-      if (dbg_scz) {
-        cout << "   Calorimeter input track has x y z t " << track->Position.X() << " " << track->Position.Y() << " " << track->Position.Z() << " " << track->Position.T() 
-             << endl;
-        Candidate *prt = static_cast<Candidate*>(track->GetCandidates()->Last());
-        const TLorentzVector &ini = prt->Position;
-
-        cout << "                and parent has x y z t " << ini.X() << " " << ini.Y() << " " << ini.Z() << " " << ini.T();
-
-      }
-      
-      if (ecalEnergy > fTimingEMin && fTower) {
-        if (fElectronsFromTrack) {
-          //      cout << " SCZ Debug pushing back track hit E=" << ecalEnergy << " T=" << track->Position.T() << " isPU=" << track->IsPU << " isRecoPU=" << track->IsRecoPU 
-          //           << " PID=" << track->PID << endl;
-          fTower->Ecal_E_t.push_back(std::make_pair<float,float>(ecalEnergy,track->Position.T()));
-        } else {
-          //      cout << " Skipping track hit E=" << ecalEnergy << " T=" << track->Position.T() << " isPU=" << track->IsPU << " isRecoPU=" << track->IsRecoPU 
-          //           << " PID=" << track->PID << endl;
-        }
-      }
-
+      fTrackECalTime += TMath::Sqrt(ecalEnergy)*position.T();
+      fTrackHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
+
+      fTrackECalTimeWeight += TMath::Sqrt(ecalEnergy);
+      fTrackHCalTimeWeight += TMath::Sqrt(hcalEnergy);
 
       fTowerTrackArray->Add(track);
@@ -412 +397 @@
       continue;
     }
-  
+
     // check for photon and electron hits in current tower
     if(flags & 2) ++fTowerPhotonHits;
@@ -427 +412 @@
     fTowerHCalEnergy += hcalEnergy;
 
-    if (ecalEnergy > fTimingEMin && fTower) {
-      if (abs(particle->PID) != 11 || !fElectronsFromTrack) {
-        //      cout << " SCZ Debug About to push back particle hit E=" << ecalEnergy << " T=" << particle->Position.T() << " isPU=" << particle->IsPU 
-        //   << " PID=" << particle->PID << endl;
-        fTower->Ecal_E_t.push_back(std::make_pair<Float_t,Float_t>(ecalEnergy,particle->Position.T()));
-      } else {
-        
-        // N.B. Only charged particle set to leave ecal energy is the electrons
-        //      cout << " SCZ Debug To avoid double-counting, skipping particle hit E=" << ecalEnergy << " T=" << particle->Position.T() << " isPU=" << particle->IsPU 
-        //           << " PID=" << particle->PID << endl;
-        
-      }
-    }
+    fTowerECalTime += TMath::Sqrt(ecalEnergy)*position.T();
+    fTowerHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
+
+    fTowerECalTimeWeight += TMath::Sqrt(ecalEnergy);
+    fTowerHCalTimeWeight += TMath::Sqrt(hcalEnergy);
+
 
     fTower->AddCandidate(particle);
@@ -456 +434 @@
   Double_t ecalEnergy, hcalEnergy;
   Double_t ecalSigma, hcalSigma;
- 
+  Double_t ecalTime, hcalTime, time;
+
   if(!fTower) return;
 
@@ -465 +444 @@
   hcalEnergy = LogNormal(fTowerHCalEnergy, hcalSigma);
 
+  ecalTime = (fTowerECalTimeWeight < 1.0E-09 ) ? 0.0 : fTowerECalTime/fTowerECalTimeWeight;
+  hcalTime = (fTowerHCalTimeWeight < 1.0E-09 ) ? 0.0 : fTowerHCalTime/fTowerHCalTimeWeight;
+
   ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy);
   hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy);
@@ -472 +454 @@
 
   energy = ecalEnergy + hcalEnergy;
-    
+  time = (TMath::Sqrt(ecalEnergy)*ecalTime + TMath::Sqrt(hcalEnergy)*hcalTime)/(TMath::Sqrt(ecalEnergy) + TMath::Sqrt(hcalEnergy));
+
   if(fSmearTowerCenter)
   {
@@ -486 +469 @@
   pt = energy / TMath::CosH(eta);
 
-  // Time calculation for tower
-  fTower->Ntimes = 0;
-  Float_t tow_sumT = 0;
-  Float_t tow_sumW = 0;
-  
-  for (Int_t i = 0 ; i < fTower->Ecal_E_t.size() ; i++) 
-  {
-    Float_t w = TMath::Sqrt(fTower->Ecal_E_t[i].first);
-    tow_sumT += w*fTower->Ecal_E_t[i].second;
-    tow_sumW += w;
-    fTower->Ntimes++;
-  }
-  
-  if (tow_sumW > 0.) {
-    fTower->Position.SetPtEtaPhiE(1.0, eta, phi,tow_sumT/tow_sumW);
-  } else {
-    fTower->Position.SetPtEtaPhiE(1.0,eta,phi,999999.);
-  }
-
-
+  fTower->Position.SetPtEtaPhiE(1.0, eta, phi, time);
   fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
   fTower->Eem = ecalEnergy;
@@ -517 +481 @@
   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;
-  hcalEnergy -= fTrackHCalEnergy;
-
-  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 = ecalEnergy / TMath::CosH(eta);
-
-    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalEnergy);
-    tower->Eem = ecalEnergy;
-    tower->Ehad = 0.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.0;
-    tower->Ehad = hcalEnergy;
-
-    fEFlowNeutralHadronOutputArray->Add(tower);
+     bool isCalPhoton = false;
+
+     if(fTowerTrackHits == 0)
+     {
+        // We have a CalPhoton when there are NOT tracks and there ARE photon hits
+        isCalPhoton = (fTowerPhotonHits > 0);
+     }
+     else
+     {
+         // save all the tracks as energy flow tracks
+         fItTowerTrackArray->Reset();
+         while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
+         {
+            fEFlowTrackOutputArray->Add(track);
+         }
+
+         ecalEnergy -= fTrackECalEnergy;
+         hcalEnergy -= fTrackHCalEnergy;
+
+         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 it's NOT a CalPhoton; add the tower as a whole entity.
+     // Otherwise, the tower will be split into its ECAL and HCAL components,
+     // then added to fPhotonArray and fTowerArray (respectively).
+     // By construction (no track hits), no track subtraction will occur
+     // for CalPhotons
+     if(!isCalPhoton)
+        fTowerOutputArray->Add(fTower);
+
+     // fill energy flow candidates
+
+     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.0;
+
+         fEFlowPhotonOutputArray->Add(tower);
+         if(isCalPhoton)
+                      fPhotonOutputArray->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.0;
+        tower->Ehad = hcalEnergy;
+
+        fEFlowNeutralHadronOutputArray->Add(tower);
+        if(isCalPhoton)
+           fTowerOutputArray->Add(tower);
+     }
   }
 }