Changeset 3db5282 in git

Timestamp:

Jan 27, 2015, 5:32:59 PM (10 years ago)

Author:

Michele Selvaggi <michele.selvaggi@…>

Branches:

ImprovedOutputFile, Timing, dual_readout, llp, master

Children:

2862770

Parents:

2d494a6

Message:

included timing information in Calorimeter

Files:

• classes/DelphesClasses.cc (modified) (4 diffs)
• classes/DelphesClasses.h (modified) (3 diffs)
• modules/Calorimeter.cc (modified) (10 diffs)
• modules/Calorimeter.h (modified) (1 diff)
• modules/Isolation.cc (modified) (1 diff)
• modules/TreeWriter.cc (modified) (1 diff)

classes/DelphesClasses.cc

@@ -133 +133 @@
   MeanSqDeltaR(0),
   PTD(0),
+  Ntimes(-1),
   IsolationVar(-999),
   IsolationVarRhoCorr(-999),
@@ -255 +256 @@
   object.MeanSqDeltaR = MeanSqDeltaR;
   object.PTD = PTD;
-  
+  object.Ntimes = Ntimes; 
   object.IsolationVar = IsolationVar;
   object.IsolationVarRhoCorr = IsolationVarRhoCorr;
@@ -276 +277 @@
   object.fFactory = fFactory;
   object.fArray = 0;
+
+
+  // Copy cluster timing info
+  copy(Ecal_E_t.begin(),Ecal_E_t.end(),back_inserter(object.Ecal_E_t));
 
   if(fArray && fArray->GetEntriesFast() > 0)
@@ -326 +331 @@
   PTD = 0.0;
   
+  Ntimes = 0;
+  Ecal_E_t.clear();
+  
   IsolationVar = -999;
   IsolationVarRhoCorr = -999;

classes/DelphesClasses.h

@@ -419 +419 @@
   Float_t E; // calorimeter tower energy
 
-  Float_t T; //particle arrival time of flight
-
+  Float_t T; // ecal deposit time, averaged by sqrt(EM energy) over all particles, not smeared
+  Int_t   Ntimes; // number of hits contributing to time measurement
+  
   Float_t Eem; // calorimeter tower electromagnetic energy
   Float_t Ehad; // calorimeter tower hadronic energy
@@ -510 +511 @@
   Float_t  PTD;
   Float_t  FracPt[5];
+  
+  //Timing information
+  
+  Int_t    Ntimes;
+  std::vector<std::pair<Float_t,Float_t> > Ecal_E_t; 
 
   // Isolation variables
@@ -523 +529 @@
 
   Float_t Tau[5];
+  
+  
+  
 
   static CompBase *fgCompare; //!

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;

modules/Calorimeter.h

@@ -60 +60 @@
   Double_t fTrackECalEnergy, fTrackHCalEnergy;
 
-  Double_t fTowerECalTime, fTowerHCalTime;
-  Double_t fTrackECalTime, fTrackHCalTime;
-
-  Double_t fTowerECalTimeWeight, fTowerHCalTimeWeight;
-  Double_t fTrackECalTimeWeight, fTrackHCalTimeWeight;
-
+  Double_t fTimingEMin;
+  Bool_t fElectronsFromTrack; // for timing
+  
   Int_t fTowerTrackHits, fTowerPhotonHits;
 

modules/Isolation.cc

@@ -184 +184 @@
     // loop over all input tracks
     
+    sumNeutral = 0.0;
     sumCharged = 0.0;
     sumChargedPU = 0.0;

modules/TreeWriter.cc

@@ -362 +362 @@
 
     entry->T = position.T()*1.0E-3/c_light;
+    entry->Ntimes = candidate->Ntimes;    
 
     FillParticles(candidate, &entry->Particles);