Fork me on GitHub

Changeset 341014c in git for modules/PileUpJetID.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/PileUpJetID.cc

    r45e58be r341014c  
    1414#include "classes/DelphesFormula.h"
    1515
     16#include "ExRootAnalysis/ExRootClassifier.h"
     17#include "ExRootAnalysis/ExRootFilter.h"
    1618#include "ExRootAnalysis/ExRootResult.h"
    17 #include "ExRootAnalysis/ExRootFilter.h"
    18 #include "ExRootAnalysis/ExRootClassifier.h"
    19 
     19
     20#include "TFormula.h"
    2021#include "TMath.h"
     22#include "TObjArray.h"
     23#include "TRandom3.h"
    2124#include "TString.h"
    22 #include "TFormula.h"
    23 #include "TRandom3.h"
    24 #include "TObjArray.h"
    2525//#include "TDatabasePDG.h"
    2626#include "TLorentzVector.h"
    2727
    2828#include <algorithm>
    29 #include <stdexcept>
    3029#include <iostream>
    3130#include <sstream>
     31#include <stdexcept>
    3232
    3333using namespace std;
     
    3636
    3737PileUpJetID::PileUpJetID() :
    38   fItJetInputArray(0),fTrackInputArray(0),fNeutralInputArray(0)
    39 {
    40 
     38  fItJetInputArray(0), fTrackInputArray(0), fNeutralInputArray(0)
     39{
    4140}
    4241
     
    4544PileUpJetID::~PileUpJetID()
    4645{
    47 
    4846}
    4947
     
    5755  fUseConstituents = GetInt("UseConstituents", 0);
    5856
    59   fMeanSqDeltaRMaxBarrel = GetDouble("MeanSqDeltaRMaxBarrel",0.1);
    60   fBetaMinBarrel = GetDouble("BetaMinBarrel",0.1);
    61   fMeanSqDeltaRMaxEndcap = GetDouble("MeanSqDeltaRMaxEndcap",0.1);
    62   fBetaMinEndcap = GetDouble("BetaMinEndcap",0.1);
    63   fMeanSqDeltaRMaxForward = GetDouble("MeanSqDeltaRMaxForward",0.1);
     57  fMeanSqDeltaRMaxBarrel = GetDouble("MeanSqDeltaRMaxBarrel", 0.1);
     58  fBetaMinBarrel = GetDouble("BetaMinBarrel", 0.1);
     59  fMeanSqDeltaRMaxEndcap = GetDouble("MeanSqDeltaRMaxEndcap", 0.1);
     60  fBetaMinEndcap = GetDouble("BetaMinEndcap", 0.1);
     61  fMeanSqDeltaRMaxForward = GetDouble("MeanSqDeltaRMaxForward", 0.1);
    6462  fJetPTMinForNeutrals = GetDouble("JetPTMinForNeutrals", 20.0);
    6563  fNeutralPTMin = GetDouble("NeutralPTMin", 2.0);
     
    8280  fOutputArray = ExportArray(GetString("OutputArray", "jets"));
    8381
    84   fNeutralsInPassingJets = ExportArray(GetString("NeutralsInPassingJets","eflowtowers"));
    85 
     82  fNeutralsInPassingJets = ExportArray(GetString("NeutralsInPassingJets", "eflowtowers"));
    8683}
    8784
     
    9592  if(fItTrackInputArray) delete fItTrackInputArray;
    9693  if(fItNeutralInputArray) delete fItNeutralInputArray;
    97 
    9894}
    9995
     
    109105  // loop over all input candidates
    110106  fItJetInputArray->Reset();
    111   while((candidate = static_cast<Candidate*>(fItJetInputArray->Next())))
     107  while((candidate = static_cast<Candidate *>(fItJetInputArray->Next())))
    112108  {
    113109    momentum = candidate->Momentum;
     
    133129    float pt_ann[5];
    134130
    135     for (int i = 0 ; i < 5 ; i++) {
     131    for(int i = 0; i < 5; i++)
     132    {
    136133      pt_ann[i] = 0.;
    137134    }
    138135
    139     if (fUseConstituents) {
     136    if(fUseConstituents)
     137    {
    140138      TIter itConstituents(candidate->GetCandidates());
    141       while((constituent = static_cast<Candidate*>(itConstituents.Next()))) {
     139      while((constituent = static_cast<Candidate *>(itConstituents.Next())))
     140      {
    142141        float pt = constituent->Momentum.Pt();
    143142        float dr = candidate->Momentum.DeltaR(constituent->Momentum);
    144         //      cout << " There exists a constituent with dr=" << dr << endl;
    145         sumpt += pt;
    146         sumdrsqptsq += dr*dr*pt*pt;
    147         sumptsq += pt*pt;
    148         if (constituent->Charge == 0) {
    149           nn++;
    150         } else {
    151           if (constituent->IsRecoPU) {
    152             sumptchpu += pt;
    153           } else {
    154             sumptchpv += pt;
    155           }
    156           sumptch += pt;
    157           nc++;
    158         }
    159         for (int i = 0 ; i < 5 ; i++) {
    160           if (dr > 0.1*i && dr < 0.1*(i+1)) {
    161             pt_ann[i] += pt;
    162           }
    163         }
    164         float tow_sumT = 0;
    165         float tow_sumW = 0;
    166         for (int i = 0 ; i < constituent->ECalEnergyTimePairs.size() ; i++) {
    167           float w = TMath::Sqrt(constituent->ECalEnergyTimePairs[i].first);
    168           if (fAverageEachTower) {
    169             tow_sumT += w*constituent->ECalEnergyTimePairs[i].second;
     143        //      cout << " There exists a constituent with dr=" << dr << endl;
     144        sumpt += pt;
     145        sumdrsqptsq += dr * dr * pt * pt;
     146        sumptsq += pt * pt;
     147        if(constituent->Charge == 0)
     148        {
     149          nn++;
     150        }
     151        else
     152        {
     153          if(constituent->IsRecoPU)
     154          {
     155            sumptchpu += pt;
     156          }
     157          else
     158          {
     159            sumptchpv += pt;
     160          }
     161          sumptch += pt;
     162          nc++;
     163        }
     164        for(int i = 0; i < 5; i++)
     165        {
     166          if(dr > 0.1 * i && dr < 0.1 * (i + 1))
     167          {
     168            pt_ann[i] += pt;
     169          }
     170        }
     171        float tow_sumT = 0;
     172        float tow_sumW = 0;
     173        for(int i = 0; i < constituent->ECalEnergyTimePairs.size(); i++)
     174        {
     175          float w = TMath::Sqrt(constituent->ECalEnergyTimePairs[i].first);
     176          if(fAverageEachTower)
     177          {
     178            tow_sumT += w * constituent->ECalEnergyTimePairs[i].second;
    170179            tow_sumW += w;
    171           } else {
    172             sumT0 += w*constituent->ECalEnergyTimePairs[i].second;
    173             sumT1 += w*gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second,0.001);
    174             sumT10 += w*gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second,0.010);
    175             sumT20 += w*gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second,0.020);
    176             sumT30 += w*gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second,0.030);
    177             sumT40 += w*gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second,0.040);
    178             sumWeightsForT += w;
    179             candidate->NTimeHits++;
    180           }
    181         }
    182         if (fAverageEachTower && tow_sumW > 0.) {
    183           sumT0 += tow_sumT;
    184           sumT1 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.001);
    185           sumT10 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0010);
    186           sumT20 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0020);
    187           sumT30 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0030);
    188           sumT40 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0040);
    189           sumWeightsForT += tow_sumW;
    190           candidate->NTimeHits++;
    191         }
    192       }
    193     } else {
     180          }
     181          else
     182          {
     183            sumT0 += w * constituent->ECalEnergyTimePairs[i].second;
     184            sumT1 += w * gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second, 0.001);
     185            sumT10 += w * gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second, 0.010);
     186            sumT20 += w * gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second, 0.020);
     187            sumT30 += w * gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second, 0.030);
     188            sumT40 += w * gRandom->Gaus(constituent->ECalEnergyTimePairs[i].second, 0.040);
     189            sumWeightsForT += w;
     190            candidate->NTimeHits++;
     191          }
     192        }
     193        if(fAverageEachTower && tow_sumW > 0.)
     194        {
     195          sumT0 += tow_sumT;
     196          sumT1 += tow_sumW * gRandom->Gaus(tow_sumT / tow_sumW, 0.001);
     197          sumT10 += tow_sumW * gRandom->Gaus(tow_sumT / tow_sumW, 0.0010);
     198          sumT20 += tow_sumW * gRandom->Gaus(tow_sumT / tow_sumW, 0.0020);
     199          sumT30 += tow_sumW * gRandom->Gaus(tow_sumT / tow_sumW, 0.0030);
     200          sumT40 += tow_sumW * gRandom->Gaus(tow_sumT / tow_sumW, 0.0040);
     201          sumWeightsForT += tow_sumW;
     202          candidate->NTimeHits++;
     203        }
     204      }
     205    }
     206    else
     207    {
    194208      // Not using constituents, using dr
    195209      fItTrackInputArray->Reset();
    196       while ((trk = static_cast<Candidate*>(fItTrackInputArray->Next()))) {
    197         if (trk->Momentum.DeltaR(candidate->Momentum) < fParameterR) {
    198           float pt = trk->Momentum.Pt();
    199           sumpt += pt;
    200           sumptch += pt;
    201           if (trk->IsRecoPU) {
    202             sumptchpu += pt;
    203           } else {
    204             sumptchpv += pt;
    205           }
    206           float dr = candidate->Momentum.DeltaR(trk->Momentum);
    207           sumdrsqptsq += dr*dr*pt*pt;
    208           sumptsq += pt*pt;
    209           nc++;
    210           for (int i = 0 ; i < 5 ; i++) {
    211             if (dr > 0.1*i && dr < 0.1*(i+1)) {
    212               pt_ann[i] += pt;
    213             }
    214           }
    215         }
     210      while((trk = static_cast<Candidate *>(fItTrackInputArray->Next())))
     211      {
     212        if(trk->Momentum.DeltaR(candidate->Momentum) < fParameterR)
     213        {
     214          float pt = trk->Momentum.Pt();
     215          sumpt += pt;
     216          sumptch += pt;
     217          if(trk->IsRecoPU)
     218          {
     219            sumptchpu += pt;
     220          }
     221          else
     222          {
     223            sumptchpv += pt;
     224          }
     225          float dr = candidate->Momentum.DeltaR(trk->Momentum);
     226          sumdrsqptsq += dr * dr * pt * pt;
     227          sumptsq += pt * pt;
     228          nc++;
     229          for(int i = 0; i < 5; i++)
     230          {
     231            if(dr > 0.1 * i && dr < 0.1 * (i + 1))
     232            {
     233              pt_ann[i] += pt;
     234            }
     235          }
     236        }
    216237      }
    217238      fItNeutralInputArray->Reset();
    218       while ((constituent = static_cast<Candidate*>(fItNeutralInputArray->Next()))) {
    219         if (constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR) {
    220           float pt = constituent->Momentum.Pt();
    221           sumpt += pt;
    222           float dr = candidate->Momentum.DeltaR(constituent->Momentum);
    223           sumdrsqptsq += dr*dr*pt*pt;
    224           sumptsq += pt*pt;
    225           nn++;
    226           for (int i = 0 ; i < 5 ; i++) {
    227             if (dr > 0.1*i && dr < 0.1*(i+1)) {
    228             pt_ann[i] += pt;
    229             }
    230           }
    231         }
    232       }
    233     }
    234 
    235     if (sumptch > 0.) {
    236       candidate->Beta = sumptchpv/sumptch;
    237       candidate->BetaStar = sumptchpu/sumptch;
    238     } else {
     239      while((constituent = static_cast<Candidate *>(fItNeutralInputArray->Next())))
     240      {
     241        if(constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR)
     242        {
     243          float pt = constituent->Momentum.Pt();
     244          sumpt += pt;
     245          float dr = candidate->Momentum.DeltaR(constituent->Momentum);
     246          sumdrsqptsq += dr * dr * pt * pt;
     247          sumptsq += pt * pt;
     248          nn++;
     249          for(int i = 0; i < 5; i++)
     250          {
     251            if(dr > 0.1 * i && dr < 0.1 * (i + 1))
     252            {
     253              pt_ann[i] += pt;
     254            }
     255          }
     256        }
     257      }
     258    }
     259
     260    if(sumptch > 0.)
     261    {
     262      candidate->Beta = sumptchpv / sumptch;
     263      candidate->BetaStar = sumptchpu / sumptch;
     264    }
     265    else
     266    {
    239267      candidate->Beta = -999.;
    240268      candidate->BetaStar = -999.;
    241269    }
    242     if (sumptsq > 0.) {
    243       candidate->MeanSqDeltaR = sumdrsqptsq/sumptsq;
    244     } else {
     270    if(sumptsq > 0.)
     271    {
     272      candidate->MeanSqDeltaR = sumdrsqptsq / sumptsq;
     273    }
     274    else
     275    {
    245276      candidate->MeanSqDeltaR = -999.;
    246277    }
    247278    candidate->NCharged = nc;
    248279    candidate->NNeutrals = nn;
    249     if (sumpt > 0.) {
     280    if(sumpt > 0.)
     281    {
    250282      candidate->PTD = TMath::Sqrt(sumptsq) / sumpt;
    251       for (int i = 0 ; i < 5 ; i++) {
    252         candidate->FracPt[i] = pt_ann[i]/sumpt;
    253       }
    254     } else {
     283      for(int i = 0; i < 5; i++)
     284      {
     285        candidate->FracPt[i] = pt_ann[i] / sumpt;
     286      }
     287    }
     288    else
     289    {
    255290      candidate->PTD = -999.;
    256       for (int i = 0 ; i < 5 ; i++) {
     291      for(int i = 0; i < 5; i++)
     292      {
    257293        candidate->FracPt[i] = -999.;
    258294      }
     
    271307
    272308    bool passId = false;
    273     if (candidate->Momentum.Pt() > fJetPTMinForNeutrals && candidate->MeanSqDeltaR > -0.1) {
    274       if (fabs(candidate->Momentum.Eta())<1.5) {
    275         passId = ((candidate->Beta > fBetaMinBarrel) && (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxBarrel));
    276       } else if (fabs(candidate->Momentum.Eta())<4.0) {
    277         passId = ((candidate->Beta > fBetaMinEndcap) && (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxEndcap));
    278       } else {
    279         passId = (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxForward);
    280       }
    281     }
    282 
    283     //    cout << " Pt Eta MeanSqDeltaR Beta PassId " << candidate->Momentum.Pt()
     309    if(candidate->Momentum.Pt() > fJetPTMinForNeutrals && candidate->MeanSqDeltaR > -0.1)
     310    {
     311      if(fabs(candidate->Momentum.Eta()) < 1.5)
     312      {
     313        passId = ((candidate->Beta > fBetaMinBarrel) && (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxBarrel));
     314      }
     315      else if(fabs(candidate->Momentum.Eta()) < 4.0)
     316      {
     317        passId = ((candidate->Beta > fBetaMinEndcap) && (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxEndcap));
     318      }
     319      else
     320      {
     321        passId = (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxForward);
     322      }
     323    }
     324
     325    //    cout << " Pt Eta MeanSqDeltaR Beta PassId " << candidate->Momentum.Pt()
    284326    //   << " " << candidate->Momentum.Eta() << " " << candidate->MeanSqDeltaR << " " << candidate->Beta << " " << passId << endl;
    285327
    286     if (passId) {
    287       if (fUseConstituents) {
    288         TIter itConstituents(candidate->GetCandidates());
    289         while((constituent = static_cast<Candidate*>(itConstituents.Next()))) {
    290           if (constituent->Charge == 0 && constituent->Momentum.Pt() > fNeutralPTMin) {
    291             fNeutralsInPassingJets->Add(constituent);
    292             //      cout << "    Constitutent added Pt Eta Charge " << constituent->Momentum.Pt() << " " << constituent->Momentum.Eta() << " " << constituent->Charge << endl;
    293           }
    294         }
    295       } else { // use DeltaR
    296         fItNeutralInputArray->Reset();
    297         while ((constituent = static_cast<Candidate*>(fItNeutralInputArray->Next()))) {
    298           if (constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR && constituent->Momentum.Pt() > fNeutralPTMin) {
    299             fNeutralsInPassingJets->Add(constituent);
    300             //            cout << "    Constitutent added Pt Eta Charge " << constituent->Momentum.Pt() << " " << constituent->Momentum.Eta() << " " << constituent->Charge << endl;
    301           }
    302         }
    303       }
    304     }
    305 
    306 
     328    if(passId)
     329    {
     330      if(fUseConstituents)
     331      {
     332        TIter itConstituents(candidate->GetCandidates());
     333        while((constituent = static_cast<Candidate *>(itConstituents.Next())))
     334        {
     335          if(constituent->Charge == 0 && constituent->Momentum.Pt() > fNeutralPTMin)
     336          {
     337            fNeutralsInPassingJets->Add(constituent);
     338            //      cout << "    Constitutent added Pt Eta Charge " << constituent->Momentum.Pt() << " " << constituent->Momentum.Eta() << " " << constituent->Charge << endl;
     339          }
     340        }
     341      }
     342      else
     343      { // use DeltaR
     344        fItNeutralInputArray->Reset();
     345        while((constituent = static_cast<Candidate *>(fItNeutralInputArray->Next())))
     346        {
     347          if(constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR && constituent->Momentum.Pt() > fNeutralPTMin)
     348          {
     349            fNeutralsInPassingJets->Add(constituent);
     350            //            cout << "    Constitutent added Pt Eta Charge " << constituent->Momentum.Pt() << " " << constituent->Momentum.Eta() << " " << constituent->Charge << endl;
     351          }
     352        }
     353      }
     354    }
    307355  }
    308356}
Note: See TracChangeset for help on using the changeset viewer.