Changeset 61dccd3 in git for modules/DualReadoutCalorimeter.cc

Timestamp:

Aug 24, 2021, 11:57:44 AM (3 years ago)

Author:

michele <michele.selvaggi@…>

Branches:

master

Children:

83e77ee

Parents:

1ad8eca

Message:

added charged deposit to calo, fix time calc in DR, add path length to tower

File:

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

modules/DualReadoutCalorimeter.cc

@@ -248 +248 @@
   number = -1;
   fTowerRmax=0.;
+
+  //cout<<"--------- new event ---------- "<<endl;
+
   while((particle = static_cast<Candidate*>(fItParticleInputArray->Next())))
   {
@@ -391 +394 @@
       fTowerTrackHits = 0;
       fTowerPhotonHits = 0;
+
+      fTowerTime = 0.0;
+      fTowerTimeWeight = 0.0;
 
       fECalTowerTrackArray->Clear();
@@ -455 +461 @@
     position = particle->Position;
 
+
     // fill current tower
     ecalEnergy = momentum.E() * fECalTowerFractions[number];
@@ -462 +469 @@
     fHCalTowerEnergy += 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()));
-      }
-    }
+    // assume combined timing measurements in ECAL/HCAL sections
+    fTowerTime += (ecalEnergy + hcalEnergy) * position.T(); //sigma_t ~ 1/sqrt(E)
+    fTowerTimeWeight += ecalEnergy + hcalEnergy;
+    //fTowerTime += (hcalEnergy) * position.T(); //sigma_t ~ 1/sqrt(E)
+    //fTowerTimeWeight += hcalEnergy;
+    //fTowerTime +=  position.T(); //sigma_t ~ 1/sqrt(E)
+    //fTowerTimeWeight += 1;
+
+    //cout<<" tower particle PID, pt, eta, phi, l, tof:  "<<particle->PID<<", "<<momentum.E()<<", "<<momentum.Eta()<<", "<<momentum.Phi()<<", "<<position.Vect().Mag()<<", "<<position.T()/2.99792458E2<<endl;
+    //cout<<" tower particle time, weight:  "<<fTowerTime/2.99792458E2<<", "<<fTowerTimeWeight<<endl;
 
     fTower->AddCandidate(particle);
@@ -484 +494 @@
 
   Candidate *track, *tower, *mother;
-  Double_t energy, pt, eta, phi, r;
+  Double_t energy, pt, eta, phi, r, time;
   Double_t ecalEnergy, hcalEnergy;
   Double_t ecalNeutralEnergy, hcalNeutralEnergy, neutralEnergy;
@@ -529 +539 @@
   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);
@@ -535 +547 @@
   if(hcalEnergy < fHCalEnergyMin || hcalEnergy < fHCalEnergySignificanceMin*hcalSigma) hcalEnergy = 0.0;
 
-  //cout<<"Measured energy: "<<energy<<endl;
-
   if(fSmearTowerCenter)
   {
@@ -550 +560 @@
   pt = energy / TMath::CosH(eta);
 
-  // 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::Power((fTower->ECalEnergyTimePairs[i].first),2);
-    sumWeightedTime += weight * fTower->ECalEnergyTimePairs[i].second;
-    sumWeight += weight;
-    fTower->NTimeHits++;
-  }
-
   // check whether barrel or endcap tower
-  if (fTower->Position.Perp() < fTowerRmax && TMath::Abs(eta) > 0.)
+  if ((fTowerRmax - fTower->Position.Perp()) < 1.e-06 && TMath::Abs(eta) > 0.)
     r = fTower->Position.Z()/TMath::SinH(eta);
   else
     r = fTower->Position.Pt();
 
-  if(sumWeight > 0.0)
-  {
-    fTower->Position.SetPtEtaPhiE(r, eta, phi, sumWeightedTime/sumWeight);
-  }
-  else
-  {
-    fTower->Position.SetPtEtaPhiE(r, eta, phi, 999999.9);
-  }
+  fTower->Position.SetPtEtaPhiE(r, eta, phi, time);
+
+  fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
+  fTower->L = fTower->Position.Vect().Mag();
+  //cout<<"   tower pt, eta, phi, l, tof:  "<<fTower->Momentum.E()<<", "<<fTower->Momentum.Eta()<<", "<<fTower->Momentum.Phi()<<", "<<fTower->L<<", "<<fTower->Position.T()/2.99792458E2<<endl;
 
   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];
@@ -602 +596 @@
   neutralSigma = neutralEnergy / TMath::Sqrt(fTrackSigma*fTrackSigma + sigma*sigma);
 
-  //cout<<"trackEnergy: "<<fTrackEnergy<<", trackSigma: "<<fTrackSigma<<", Ntracks: "<<fTowerTrackArray->GetEntries()<<endl;
-
-  //cout<<"neutralEnergy: "<<neutralEnergy<<", neutralSigma: "<<neutralSigma<<", :fEnergyMin "<<fEnergyMin<<", fEnergySignificanceMin: "<<fEnergySignificanceMin<<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);