Fork me on GitHub

Changeset 341014c in git for modules/SimpleCalorimeter.cc


Ignore:
Timestamp:
Feb 12, 2019, 9:29:17 PM (6 years ago)
Author:
Pavel Demin <pavel-demin@…>
Branches:
ImprovedOutputFile, Timing, llp, master
Children:
6455202
Parents:
45e58be
Message:

apply .clang-format to all .h, .cc and .cpp files

File:
1 edited

Legend:

Unmodified
Added
Removed
  • modules/SimpleCalorimeter.cc

    r45e58be r341014c  
    1717 */
    1818
    19 
    2019/** \class SimpleCalorimeter
    2120 *
     
    3332#include "classes/DelphesFormula.h"
    3433
     34#include "ExRootAnalysis/ExRootClassifier.h"
     35#include "ExRootAnalysis/ExRootFilter.h"
    3536#include "ExRootAnalysis/ExRootResult.h"
    36 #include "ExRootAnalysis/ExRootFilter.h"
    37 #include "ExRootAnalysis/ExRootClassifier.h"
    38 
     37
     38#include "TDatabasePDG.h"
     39#include "TFormula.h"
     40#include "TLorentzVector.h"
    3941#include "TMath.h"
     42#include "TObjArray.h"
     43#include "TRandom3.h"
    4044#include "TString.h"
    41 #include "TFormula.h"
    42 #include "TRandom3.h"
    43 #include "TObjArray.h"
    44 #include "TDatabasePDG.h"
    45 #include "TLorentzVector.h"
    4645
    4746#include <algorithm>
    48 #include <stdexcept>
    4947#include <iostream>
    5048#include <sstream>
     49#include <stdexcept>
    5150
    5251using namespace std;
     
    5857  fItParticleInputArray(0), fItTrackInputArray(0)
    5958{
    60  
     59
    6160  fResolutionFormula = new DelphesFormula;
    6261  fTowerTrackArray = new TObjArray;
    6362  fItTowerTrackArray = fTowerTrackArray->MakeIterator();
    64  
    6563}
    6664
     
    6967SimpleCalorimeter::~SimpleCalorimeter()
    7068{
    71  
     69
    7270  if(fResolutionFormula) delete fResolutionFormula;
    7371  if(fTowerTrackArray) delete fTowerTrackArray;
    7472  if(fItTowerTrackArray) delete fItTowerTrackArray;
    75  
    7673}
    7774
     
    8481  Double_t fraction;
    8582  TBinMap::iterator itEtaBin;
    86   set< Double_t >::iterator itPhiBin;
    87   vector< Double_t > *phiBins;
     83  set<Double_t>::iterator itPhiBin;
     84  vector<Double_t> *phiBins;
    8885
    8986  // read eta and phi bins
     
    9390  fEtaBins.clear();
    9491  fPhiBins.clear();
    95   for(i = 0; i < size/2; ++i)
    96   {
    97     paramEtaBins = param[i*2];
     92  for(i = 0; i < size / 2; ++i)
     93  {
     94    paramEtaBins = param[i * 2];
    9895    sizeEtaBins = paramEtaBins.GetSize();
    99     paramPhiBins = param[i*2 + 1];
     96    paramPhiBins = param[i * 2 + 1];
    10097    sizePhiBins = paramPhiBins.GetSize();
    10198
     
    114111  {
    115112    fEtaBins.push_back(itEtaBin->first);
    116     phiBins = new vector< double >(itEtaBin->second.size());
     113    phiBins = new vector<double>(itEtaBin->second.size());
    117114    fPhiBins.push_back(phiBins);
    118115    phiBins->clear();
     
    131128  fFractionMap[0] = 1.0;
    132129
    133   for(i = 0; i < size/2; ++i)
    134   {
    135     paramFractions = param[i*2 + 1];
     130  for(i = 0; i < size / 2; ++i)
     131  {
     132    paramFractions = param[i * 2 + 1];
    136133    fraction = paramFractions[0].GetDouble();
    137     fFractionMap[param[i*2].GetInt()] = fraction;
     134    fFractionMap[param[i * 2].GetInt()] = fraction;
    138135  }
    139136
     
    170167void SimpleCalorimeter::Finish()
    171168{
    172   vector< vector< Double_t >* >::iterator itPhiBin;
     169  vector<vector<Double_t> *>::iterator itPhiBin;
    173170  if(fItParticleInputArray) delete fItParticleInputArray;
    174171  if(fItTrackInputArray) delete fItTrackInputArray;
     
    197194  TFractionMap::iterator itFractionMap;
    198195
    199   vector< Double_t >::iterator itEtaBin;
    200   vector< Double_t >::iterator itPhiBin;
    201   vector< Double_t > *phiBins;
    202 
    203   vector< Long64_t >::iterator itTowerHits;
     196  vector<Double_t>::iterator itEtaBin;
     197  vector<Double_t>::iterator itPhiBin;
     198  vector<Double_t> *phiBins;
     199
     200  vector<Long64_t>::iterator itTowerHits;
    204201
    205202  DelphesFactory *factory = GetFactory();
     
    211208  fItParticleInputArray->Reset();
    212209  number = -1;
    213   while((particle = static_cast<Candidate*>(fItParticleInputArray->Next())))
     210  while((particle = static_cast<Candidate *>(fItParticleInputArray->Next())))
    214211  {
    215212    const TLorentzVector &particlePosition = particle->Position;
     
    254251  fItTrackInputArray->Reset();
    255252  number = -1;
    256   while((track = static_cast<Candidate*>(fItTrackInputArray->Next())))
     253  while((track = static_cast<Candidate *>(fItTrackInputArray->Next())))
    257254  {
    258255    const TLorentzVector &trackPosition = track->Position;
     
    303300    towerHit = (*itTowerHits);
    304301    flags = (towerHit >> 24) & 0x00000000000000FFLL;
    305     number = (towerHit) & 0x0000000000FFFFFFLL;
     302    number = (towerHit)&0x0000000000FFFFFFLL;
    306303    hitEtaPhi = towerHit >> 32;
    307304
     
    324321
    325322      // calculate eta and phi of the tower's center
    326       fTowerEta = 0.5*(fEtaBins[etaBin - 1] + fEtaBins[etaBin]);
    327       fTowerPhi = 0.5*((*phiBins)[phiBin - 1] + (*phiBins)[phiBin]);
     323      fTowerEta = 0.5 * (fEtaBins[etaBin - 1] + fEtaBins[etaBin]);
     324      fTowerPhi = 0.5 * ((*phiBins)[phiBin - 1] + (*phiBins)[phiBin]);
    328325
    329326      fTowerEdges[0] = fEtaBins[etaBin - 1];
     
    346343
    347344      fTowerTrackArray->Clear();
    348      }
     345    }
    349346
    350347    // check for track hits
     
    353350      ++fTowerTrackHits;
    354351
    355       track = static_cast<Candidate*>(fTrackInputArray->At(number));
     352      track = static_cast<Candidate *>(fTrackInputArray->At(number));
    356353      momentum = track->Momentum;
    357354      position = track->Position;
     
    359356      energy = momentum.E() * fTrackFractions[number];
    360357
    361       fTrackTime += TMath::Sqrt(energy)*position.T();
     358      fTrackTime += TMath::Sqrt(energy) * position.T();
    362359      fTrackTimeWeight += TMath::Sqrt(energy);
    363360
    364361      if(fTrackFractions[number] > 1.0E-9)
    365362      {
    366              
    367        // compute total charged energy   
    368        fTrackEnergy += energy;
    369        sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
    370        if(sigma/momentum.E() < track->TrackResolution) energyGuess = energy;
    371        else energyGuess = momentum.E();
    372 
    373        fTrackSigma += ((track->TrackResolution)*energyGuess)*((track->TrackResolution)*energyGuess);
    374        fTowerTrackArray->Add(track);
     363
     364        // compute total charged energy
     365        fTrackEnergy += energy;
     366        sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
     367        if(sigma / momentum.E() < track->TrackResolution)
     368          energyGuess = energy;
     369        else
     370          energyGuess = momentum.E();
     371
     372        fTrackSigma += ((track->TrackResolution) * energyGuess) * ((track->TrackResolution) * energyGuess);
     373        fTowerTrackArray->Add(track);
    375374      }
    376        
     375
    377376      else
    378377      {
     
    386385    if(flags & 2) ++fTowerPhotonHits;
    387386
    388     particle = static_cast<Candidate*>(fParticleInputArray->At(number));
     387    particle = static_cast<Candidate *>(fParticleInputArray->At(number));
    389388    momentum = particle->Momentum;
    390389    position = particle->Position;
     
    395394    fTowerEnergy += energy;
    396395
    397     fTowerTime += energy*position.T();
     396    fTowerTime += energy * position.T();
    398397    fTowerTimeWeight += energy;
    399398
     
    410409{
    411410  Candidate *tower, *track, *mother;
    412   Double_t energy,neutralEnergy, pt, eta, phi;
     411  Double_t energy, neutralEnergy, pt, eta, phi;
    413412  Double_t sigma, neutralSigma;
    414413  Double_t time;
    415    
     414
    416415  Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor;
    417416
     
    425424  energy = LogNormal(fTowerEnergy, sigma);
    426425
    427   time = (fTowerTimeWeight < 1.0E-09 ) ? 0.0 : fTowerTime/fTowerTimeWeight;
     426  time = (fTowerTimeWeight < 1.0E-09) ? 0.0 : fTowerTime / fTowerTimeWeight;
    428427
    429428  sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, energy);
    430429
    431   if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0;
    432 
     430  if(energy < fEnergyMin || energy < fEnergySignificanceMin * sigma) energy = 0.0;
    433431
    434432  if(fSmearTowerCenter)
     
    459457  if(energy > 0.0) fTowerOutputArray->Add(fTower);
    460458
    461  
    462459  // e-flow candidates
    463460
    464461  //compute neutral excess
    465  
     462
    466463  fTrackSigma = TMath::Sqrt(fTrackSigma);
    467   neutralEnergy = max( (energy - fTrackEnergy) , 0.0);
    468  
     464  neutralEnergy = max((energy - fTrackEnergy), 0.0);
     465
    469466  //compute sigma_trk total
    470   neutralSigma = neutralEnergy / TMath::Sqrt(fTrackSigma*fTrackSigma+ sigma*sigma);
    471    
     467  neutralSigma = neutralEnergy / TMath::Sqrt(fTrackSigma * fTrackSigma + sigma * sigma);
     468
    472469  // if neutral excess is significant, simply create neutral Eflow tower and clone each track into eflowtrack
    473470  if(neutralEnergy > fEnergyMin && neutralSigma > fEnergySignificanceMin)
    474471  {
    475472    // create new photon tower
    476     tower = static_cast<Candidate*>(fTower->Clone());
     473    tower = static_cast<Candidate *>(fTower->Clone());
    477474    pt = neutralEnergy / TMath::CosH(eta);
    478475
     
    480477    tower->Ehad = (fIsEcal) ? 0 : neutralEnergy;
    481478    tower->PID = (fIsEcal) ? 22 : 0;
    482  
     479
    483480    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, neutralEnergy);
    484481    fEFlowTowerOutputArray->Add(tower);
    485    
     482
    486483    fItTowerTrackArray->Reset();
    487     while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
     484    while((track = static_cast<Candidate *>(fItTowerTrackArray->Next())))
    488485    {
    489486      mother = track;
    490       track = static_cast<Candidate*>(track->Clone());
     487      track = static_cast<Candidate *>(track->Clone());
    491488      track->AddCandidate(mother);
    492489
     
    494491    }
    495492  }
    496    
     493
    497494  // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the calorimeter and tracking
    498   else if (fTrackEnergy > 0.0)
    499   {
    500     weightTrack = (fTrackSigma > 0.0) ? 1 / (fTrackSigma*fTrackSigma) : 0.0;
    501     weightCalo  = (sigma > 0.0) ? 1 / (sigma*sigma) : 0.0;
    502  
    503     bestEnergyEstimate = (weightTrack*fTrackEnergy + weightCalo*energy) / (weightTrack + weightCalo);
    504     rescaleFactor = bestEnergyEstimate/fTrackEnergy;
    505    
     495  else if(fTrackEnergy > 0.0)
     496  {
     497    weightTrack = (fTrackSigma > 0.0) ? 1 / (fTrackSigma * fTrackSigma) : 0.0;
     498    weightCalo = (sigma > 0.0) ? 1 / (sigma * sigma) : 0.0;
     499
     500    bestEnergyEstimate = (weightTrack * fTrackEnergy + weightCalo * energy) / (weightTrack + weightCalo);
     501    rescaleFactor = bestEnergyEstimate / fTrackEnergy;
     502
    506503    fItTowerTrackArray->Reset();
    507     while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
     504    while((track = static_cast<Candidate *>(fItTowerTrackArray->Next())))
    508505    {
    509506      mother = track;
    510       track = static_cast<Candidate*>(track->Clone());
     507      track = static_cast<Candidate *>(track->Clone());
    511508      track->AddCandidate(mother);
    512509
     
    516513    }
    517514  }
    518    
    519  }
     515}
    520516
    521517//------------------------------------------------------------------------------
     
    527523  if(mean > 0.0)
    528524  {
    529     b = TMath::Sqrt(TMath::Log((1.0 + (sigma*sigma)/(mean*mean))));
    530     a = TMath::Log(mean) - 0.5*b*b;
    531 
    532     return TMath::Exp(a + b*gRandom->Gaus(0.0, 1.0));
     525    b = TMath::Sqrt(TMath::Log((1.0 + (sigma * sigma) / (mean * mean))));
     526    a = TMath::Log(mean) - 0.5 * b * b;
     527
     528    return TMath::Exp(a + b * gRandom->Gaus(0.0, 1.0));
    533529  }
    534530  else
Note: See TracChangeset for help on using the changeset viewer.