Changes in modules/Calorimeter.cc [e33e6db:3db5282] 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);
@@ -356 +364 @@
       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;
@@ -380 +379 @@
       position = track->Position;
 
-
       ecalEnergy = momentum.E() * fTrackECalFractions[number];
       hcalEnergy = momentum.E() * fTrackHCalFractions[number];
@@ -387 +385 @@
       fTrackHCalEnergy += hcalEnergy;
 
-      fTrackECalTime += TMath::Sqrt(ecalEnergy)*position.T();
-      fTrackHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
-
-      fTrackECalTimeWeight += TMath::Sqrt(ecalEnergy);
-      fTrackHCalTimeWeight += TMath::Sqrt(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;
+        }
+      }
+
 
       fTowerTrackArray->Add(track);
@@ -397 +412 @@
       continue;
     }
-
+  
     // check for photon and electron hits in current tower
     if(flags & 2) ++fTowerPhotonHits;
@@ -412 +427 @@
     fTowerHCalEnergy += hcalEnergy;
 
-    fTowerECalTime += TMath::Sqrt(ecalEnergy)*position.T();
-    fTowerHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
-
-    fTowerECalTimeWeight += TMath::Sqrt(ecalEnergy);
-    fTowerHCalTimeWeight += TMath::Sqrt(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;
+        
+      }
+    }
 
     fTower->AddCandidate(particle);
@@ -434 +456 @@
   Double_t ecalEnergy, hcalEnergy;
   Double_t ecalSigma, hcalSigma;
-  Double_t ecalTime, hcalTime, time;
-
+ 
   if(!fTower) return;
 
@@ -444 +465 @@
   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);
@@ -454 +472 @@
 
   energy = ecalEnergy + hcalEnergy;
-  time = (TMath::Sqrt(ecalEnergy)*ecalTime + TMath::Sqrt(hcalEnergy)*hcalTime)/(TMath::Sqrt(ecalEnergy) + TMath::Sqrt(hcalEnergy));
-
+    
   if(fSmearTowerCenter)
   {
@@ -469 +486 @@
   pt = energy / TMath::CosH(eta);
 
-  fTower->Position.SetPtEtaPhiE(1.0, eta, phi, time);
+  // 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->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
   fTower->Eem = ecalEnergy;
@@ -481 +517 @@
   if(energy > 0.0)
   {
-     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);
-     }
+    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);
   }
 }