Fork me on GitHub

Changeset 04290b1 in git


Ignore:
Timestamp:
Jun 21, 2016, 4:17:30 PM (9 years ago)
Author:
Michele Selvaggi <michele.selvaggi@…>
Branches:
ImprovedOutputFile, Timing, dual_readout, llp, master
Children:
efac6f9, f688c89
Parents:
d97b2af
Message:

new particle-flow

Location:
modules
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • modules/Calorimeter.cc

    rd97b2af r04290b1  
    5858  fItParticleInputArray(0), fItTrackInputArray(0)
    5959{
    60   Int_t i;
    61 
     60 
    6261  fECalResolutionFormula = new DelphesFormula;
    6362  fHCalResolutionFormula = new DelphesFormula;
    6463
    65   for(i = 0; i < 2; ++i)
    66   {
    67     fECalTowerTrackArray[i] = new TObjArray;
    68     fItECalTowerTrackArray[i] = fECalTowerTrackArray[i]->MakeIterator();
    69 
    70     fHCalTowerTrackArray[i] = new TObjArray;
    71     fItHCalTowerTrackArray[i] = fHCalTowerTrackArray[i]->MakeIterator();
    72   }
     64  fECalTowerTrackArray = new TObjArray;
     65  fItECalTowerTrackArray = fECalTowerTrackArray->MakeIterator();
     66
     67  fHCalTowerTrackArray = new TObjArray;
     68  fItHCalTowerTrackArray = fHCalTowerTrackArray->MakeIterator();
     69 
    7370}
    7471
     
    7774Calorimeter::~Calorimeter()
    7875{
    79   Int_t i;
    80 
     76 
    8177  if(fECalResolutionFormula) delete fECalResolutionFormula;
    8278  if(fHCalResolutionFormula) delete fHCalResolutionFormula;
    8379
    84   for(i = 0; i < 2; ++i)
    85   {
    86     if(fECalTowerTrackArray[i]) delete fECalTowerTrackArray[i];
    87     if(fItECalTowerTrackArray[i]) delete fItECalTowerTrackArray[i];
    88 
    89     if(fHCalTowerTrackArray[i]) delete fHCalTowerTrackArray[i];
    90     if(fItHCalTowerTrackArray[i]) delete fItHCalTowerTrackArray[i];
    91   }
     80  if(fECalTowerTrackArray) delete fECalTowerTrackArray;
     81  if(fItECalTowerTrackArray) delete fItECalTowerTrackArray;
     82
     83  if(fHCalTowerTrackArray) delete fHCalTowerTrackArray;
     84  if(fItHCalTowerTrackArray) delete fItHCalTowerTrackArray;
     85 
    9286}
    9387
     
    218212  Double_t ecalFraction, hcalFraction;
    219213  Double_t ecalEnergy, hcalEnergy;
    220   Double_t ecalSigma, hcalSigma;
    221214  Int_t pdgCode;
    222215
     
    368361      fHCalTowerEnergy = 0.0;
    369362
    370       fECalTrackEnergy[0] = 0.0;
    371       fECalTrackEnergy[1] = 0.0;
    372 
    373       fHCalTrackEnergy[0] = 0.0;
    374       fHCalTrackEnergy[1] = 0.0;
    375 
     363      fECalTrackEnergy = 0.0;
     364      fHCalTrackEnergy = 0.0;
     365     
     366      fECalTrackSigma = 0.0;
     367      fHCalTrackSigma = 0.0;
     368     
    376369      fTowerTrackHits = 0;
    377370      fTowerPhotonHits = 0;
    378371
    379       fECalTowerTrackArray[0]->Clear();
    380       fECalTowerTrackArray[1]->Clear();
    381 
    382       fHCalTowerTrackArray[0]->Clear();
    383       fHCalTowerTrackArray[1]->Clear();
     372      fECalTowerTrackArray->Clear();
     373      fHCalTowerTrackArray->Clear();
     374   
    384375    }
    385376
     
    406397      if(fECalTrackFractions[number] > 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
    407398      {
    408         ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
    409         if(ecalSigma/momentum.E() < track->TrackResolution)
    410         {
    411           fECalTrackEnergy[0] += ecalEnergy;
    412           fECalTowerTrackArray[0]->Add(track);
    413         }
    414         else
    415         {
    416           fECalTrackEnergy[1] += ecalEnergy;
    417           fECalTowerTrackArray[1]->Add(track);
    418         }
     399        fECalTrackEnergy += ecalEnergy;
     400        fECalTrackSigma += (track->TrackResolution)*momentum.E()*(track->TrackResolution)*momentum.E();
     401        fECalTowerTrackArray->Add(track);
    419402      }
     403     
    420404      else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] > 1.0E-9)
    421405      {
    422         hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
    423         if(hcalSigma/momentum.E() < track->TrackResolution)
    424         {
    425           fHCalTrackEnergy[0] += hcalEnergy;
    426           fHCalTowerTrackArray[0]->Add(track);
    427         }
    428         else
    429         {
    430           fHCalTrackEnergy[1] += hcalEnergy;
    431           fHCalTowerTrackArray[1]->Add(track);
    432         }
     406        fHCalTrackEnergy += hcalEnergy;
     407        fHCalTrackSigma += (track->TrackResolution)*momentum.E()*(track->TrackResolution)*momentum.E();
     408        fHCalTowerTrackArray->Add(track);
    433409      }
     410     
    434411      else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
    435412      {
     
    476453  Double_t energy, pt, eta, phi;
    477454  Double_t ecalEnergy, hcalEnergy;
     455  Double_t ecalNeutralEnergy, hcalNeutralEnergy;
     456 
    478457  Double_t ecalSigma, hcalSigma;
    479 
     458  Double_t ecalNeutralSigma, hcalNeutralSigma;
     459
     460  Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor;
     461 
    480462  TLorentzVector momentum;
    481463  TFractionMap::iterator itFractionMap;
     
    484466
    485467  if(!fTower) return;
    486 
    487468
    488469  ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fECalTowerEnergy);
     
    554535    fTowerOutputArray->Add(fTower);
    555536  }
    556 
     537 
    557538  // fill energy flow candidates
    558 
    559   ecalEnergy -= fECalTrackEnergy[1];
    560   hcalEnergy -= fHCalTrackEnergy[1];
    561 
    562   fItECalTowerTrackArray[0]->Reset();
    563   while((track = static_cast<Candidate*>(fItECalTowerTrackArray[0]->Next())))
    564   {
    565     mother = track;
    566     track = static_cast<Candidate*>(track->Clone());
    567     track->AddCandidate(mother);
    568 
    569     track->Momentum *= ecalEnergy/fECalTrackEnergy[0];
    570 
    571     fEFlowTrackOutputArray->Add(track);
    572   }
    573 
    574   fItECalTowerTrackArray[1]->Reset();
    575   while((track = static_cast<Candidate*>(fItECalTowerTrackArray[1]->Next())))
    576   {
    577     mother = track;
    578     track = static_cast<Candidate*>(track->Clone());
    579     track->AddCandidate(mother);
    580 
    581     fEFlowTrackOutputArray->Add(track);
    582   }
    583 
    584   fItHCalTowerTrackArray[0]->Reset();
    585   while((track = static_cast<Candidate*>(fItHCalTowerTrackArray[0]->Next())))
    586   {
    587     mother = track;
    588     track = static_cast<Candidate*>(track->Clone());
    589     track->AddCandidate(mother);
    590 
    591     track->Momentum *= hcalEnergy/fHCalTrackEnergy[0];
    592 
    593     fEFlowTrackOutputArray->Add(track);
    594   }
    595 
    596   fItHCalTowerTrackArray[1]->Reset();
    597   while((track = static_cast<Candidate*>(fItHCalTowerTrackArray[1]->Next())))
    598   {
    599     mother = track;
    600     track = static_cast<Candidate*>(track->Clone());
    601     track->AddCandidate(mother);
    602 
    603     fEFlowTrackOutputArray->Add(track);
    604   }
    605 
    606   if(fECalTowerTrackArray[0]->GetEntriesFast() > 0) ecalEnergy = 0.0;
    607   if(fHCalTowerTrackArray[0]->GetEntriesFast() > 0) hcalEnergy = 0.0;
    608 
    609   ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy);
    610   hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy);
    611 
    612   if(ecalEnergy < fECalEnergyMin || ecalEnergy < fECalEnergySignificanceMin*ecalSigma) ecalEnergy = 0.0;
    613   if(hcalEnergy < fHCalEnergyMin || hcalEnergy < fHCalEnergySignificanceMin*hcalSigma) hcalEnergy = 0.0;
    614 
    615   energy = ecalEnergy + hcalEnergy;
    616 
    617   if(ecalEnergy > 0.0)
     539  fECalTrackSigma = TMath::Sqrt(fECalTrackSigma);
     540  fHCalTrackSigma = TMath::Sqrt(fHCalTrackSigma);
     541
     542  //compute neutral excesses
     543  ecalNeutralEnergy = max( (ecalEnergy - fECalTrackEnergy) , 0.0);
     544  hcalNeutralEnergy = max( (hcalEnergy - fHCalTrackEnergy) , 0.0);
     545 
     546  ecalNeutralSigma = ecalNeutralEnergy / TMath::Sqrt(fECalTrackSigma*fECalTrackSigma + ecalSigma*ecalSigma);
     547  hcalNeutralSigma = hcalNeutralEnergy / TMath::Sqrt(fHCalTrackSigma*fHCalTrackSigma + hcalSigma*hcalSigma);
     548 
     549   // if ecal neutral excess is significant, simply create neutral EflowPhoton tower and clone each track into eflowtrack
     550  if(ecalNeutralEnergy > fECalEnergyMin && ecalNeutralSigma > fECalEnergySignificanceMin)
    618551  {
    619552    // create new photon tower
    620553    tower = static_cast<Candidate*>(fTower->Clone());
    621 
    622     pt = ecalEnergy / TMath::CosH(eta);
    623 
    624     tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalEnergy);
    625     tower->Eem = ecalEnergy;
     554    pt =  ecalNeutralEnergy / TMath::CosH(eta);
     555   
     556    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalNeutralEnergy);
     557    tower->Eem = ecalNeutralEnergy;
    626558    tower->Ehad = 0.0;
    627     tower->PID = 22;
    628 
     559   
    629560    fEFlowPhotonOutputArray->Add(tower);
    630   }
    631   if(hcalEnergy > 0.0)
    632   {
    633     // create new neutral hadron tower
     561   
     562    //clone tracks
     563    fItECalTowerTrackArray->Reset();
     564    while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next())))
     565    {
     566      mother = track;
     567      track = static_cast<Candidate*>(track->Clone());
     568      track->AddCandidate(mother);
     569
     570      fEFlowTrackOutputArray->Add(track);
     571    }
     572 
     573  }
     574 
     575  // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the calorimeter and tracking
     576  else if(fECalTrackEnergy > 0.0)
     577  {
     578    weightTrack = (fECalTrackSigma > 0.0) ? 1 / (fECalTrackSigma*fECalTrackSigma) : 0.0;
     579    weightCalo  = (ecalSigma > 0.0) ? 1 / (ecalSigma*ecalSigma) : 0.0;
     580 
     581    bestEnergyEstimate = (weightTrack*fECalTrackEnergy + weightCalo*ecalEnergy) / (weightTrack + weightCalo);
     582    rescaleFactor = bestEnergyEstimate/fECalTrackEnergy;
     583
     584    //rescale tracks
     585    fItECalTowerTrackArray->Reset();
     586    while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next())))
     587    { 
     588      mother = track;
     589      track = static_cast<Candidate*>(track->Clone());
     590      track->AddCandidate(mother);
     591
     592      track->Momentum *= rescaleFactor;
     593
     594      fEFlowTrackOutputArray->Add(track);
     595    }
     596  }
     597
     598
     599  // if hcal neutral excess is significant, simply create neutral EflowNeutralHadron tower and clone each track into eflowtrack
     600  if(hcalNeutralEnergy > fHCalEnergyMin && hcalNeutralSigma > fHCalEnergySignificanceMin)
     601  {
     602    // create new photon tower
    634603    tower = static_cast<Candidate*>(fTower->Clone());
    635 
    636     pt = hcalEnergy / TMath::CosH(eta);
    637 
    638     tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalEnergy);
     604    pt =  hcalNeutralEnergy / TMath::CosH(eta);
     605   
     606    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalNeutralEnergy);
     607    tower->Ehad = hcalNeutralEnergy;
    639608    tower->Eem = 0.0;
    640     tower->Ehad = hcalEnergy;
    641 
     609   
    642610    fEFlowNeutralHadronOutputArray->Add(tower);
    643   }
     611   
     612    //clone tracks
     613    fItHCalTowerTrackArray->Reset();
     614    while((track = static_cast<Candidate*>(fItHCalTowerTrackArray->Next())))
     615    {
     616      mother = track;
     617      track = static_cast<Candidate*>(track->Clone());
     618      track->AddCandidate(mother);
     619
     620      fEFlowTrackOutputArray->Add(track);
     621    }
     622 
     623  }
     624 
     625  // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the calorimeter and tracking
     626  else if(fHCalTrackEnergy > 0.0)
     627  {
     628    weightTrack = (fHCalTrackSigma > 0.0) ? 1 / (fHCalTrackSigma*fHCalTrackSigma) : 0.0;
     629    weightCalo  = (hcalSigma > 0.0) ? 1 / (hcalSigma*hcalSigma) : 0.0;
     630 
     631    bestEnergyEstimate = (weightTrack*fHCalTrackEnergy + weightCalo*hcalEnergy) / (weightTrack + weightCalo);
     632    rescaleFactor = bestEnergyEstimate / fHCalTrackEnergy;
     633
     634    //rescale tracks
     635    fItECalTowerTrackArray->Reset();
     636    while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next())))
     637    { 
     638      mother = track;
     639      track = static_cast<Candidate*>(track->Clone());
     640      track->AddCandidate(mother);
     641
     642      track->Momentum *= rescaleFactor;
     643
     644      fEFlowTrackOutputArray->Add(track);
     645    }
     646  }
     647 
     648 
    644649}
    645650
  • modules/Calorimeter.h

    rd97b2af r04290b1  
    5858  Double_t fTowerEta, fTowerPhi, fTowerEdges[4];
    5959  Double_t fECalTowerEnergy, fHCalTowerEnergy;
    60   Double_t fECalTrackEnergy[2], fHCalTrackEnergy[2];
     60  Double_t fECalTrackEnergy, fHCalTrackEnergy;
    6161
    6262  Double_t fTimingEnergyMin;
     
    7070  Double_t fECalEnergySignificanceMin;
    7171  Double_t fHCalEnergySignificanceMin;
     72
     73  Double_t fECalTrackSigma;
     74  Double_t fHCalTrackSigma;
    7275
    7376  Bool_t fSmearTowerCenter;
     
    103106  TObjArray *fEFlowNeutralHadronOutputArray; //!
    104107
    105   TObjArray *fECalTowerTrackArray[2]; //!
    106   TIterator *fItECalTowerTrackArray[2]; //!
     108  TObjArray *fECalTowerTrackArray; //!
     109  TIterator *fItECalTowerTrackArray; //!
    107110
    108   TObjArray *fHCalTowerTrackArray[2]; //!
    109   TIterator *fItHCalTowerTrackArray[2]; //!
     111  TObjArray *fHCalTowerTrackArray; //!
     112  TIterator *fItHCalTowerTrackArray; //!
    110113
    111114  void FinalizeTower();
Note: See TracChangeset for help on using the changeset viewer.