Fork me on GitHub

Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed
  • modules/PileUpJetID.cc

    r6cdc544 r01f457a  
    22 *  Delphes: a framework for fast simulation of a generic collider experiment
    33 *  Copyright (C) 2012-2014  Universite catholique de Louvain (UCL), Belgium
    4  *
     4 * 
    55 *  This program is free software: you can redistribute it and/or modify
    66 *  it under the terms of the GNU General Public License as published by
    77 *  the Free Software Foundation, either version 3 of the License, or
    88 *  (at your option) any later version.
    9  *
     9 * 
    1010 *  This program is distributed in the hope that it will be useful,
    1111 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
    1212 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    1313 *  GNU General Public License for more details.
    14  *
     14 * 
    1515 *  You should have received a copy of the GNU General Public License
    1616 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
    1717 */
    1818
    19 
    2019/** \class PileUpJetID
    2120 *
     
    2322 *
    2423 *  \author S. Zenz, December 2013
     24 *
    2525 *
    2626 */
     
    5454
    5555PileUpJetID::PileUpJetID() :
    56   fItJetInputArray(0),fTrackInputArray(0),fNeutralInputArray(0),fItVertexInputArray(0)
     56  fItJetInputArray(0),fTrackInputArray(0),fNeutralInputArray(0),fItVertexInputArray(0) 
    5757{
    5858
     
    8686  fNeutralInputArray = ImportArray(GetString("NeutralInputArray", "Calorimeter/eflowTowers"));
    8787  fItNeutralInputArray = fNeutralInputArray->MakeIterator();
    88 
     88 
    8989  fVertexInputArray = ImportArray(GetString("VertexInputArray", "PileUpMerger/vertices"));
    9090  fItVertexInputArray = fVertexInputArray->MakeIterator();
    91 
     91 
    9292  fZVertexResolution  = GetDouble("ZVertexResolution", 0.005)*1.0E3;
    93 
    94   // create output array(s)
     93// create output array(s)
    9594
    9695  fOutputArray = ExportArray(GetString("OutputArray", "jets"));
     96
    9797}
    9898
     
    101101void PileUpJetID::Finish()
    102102{
     103
    103104  if(fItJetInputArray) delete fItJetInputArray;
    104105  if(fItTrackInputArray) delete fItTrackInputArray;
    105106  if(fItNeutralInputArray) delete fItNeutralInputArray;
    106107  if(fItVertexInputArray) delete fItVertexInputArray;
     108
    107109}
    108110
     
    113115  Candidate *candidate, *constituent;
    114116  TLorentzVector momentum, area;
    115   Int_t i, nc, nn;
    116   Double_t sumpt, sumptch, sumptchpv, sumptchpu, sumdrsqptsq, sumptsq;
    117   Double_t dr, pt, pt_ann[5];
    118   Double_t zvtx = 0.0;
    119 
    120   Candidate *track;
    121 
    122   // find z position of primary vertex
    123 
     117  Double_t zvtx=0;
     118
     119  Candidate *trk;
     120
     121 // find z position of primary vertex
     122   
    124123  fItVertexInputArray->Reset();
    125124  while((candidate = static_cast<Candidate*>(fItVertexInputArray->Next())))
     
    127126    if(!candidate->IsPU)
    128127    {
    129       zvtx = candidate->Position.Z();
    130       break;
     128    zvtx = candidate->Position.Z();
     129    break;
    131130    }
    132131  }
     
    139138    area = candidate->Area;
    140139
    141     sumpt = 0.0;
    142     sumptch = 0.0;
    143     sumptchpv = 0.0;
    144     sumptchpu = 0.0;
    145     sumdrsqptsq = 0.0;
    146     sumptsq = 0.0;
    147     nc = 0;
    148     nn = 0;
    149 
    150     for(i = 0; i < 5; ++i)
    151     {
    152       pt_ann[i] = 0.0;
    153     }
    154 
    155     if(fUseConstituents)
    156     {
     140    float sumpt = 0.;
     141    float sumptch = 0.;
     142    float sumptchpv = 0.;
     143    float sumptchpu = 0.;
     144    float sumdrsqptsq = 0.;
     145    float sumptsq = 0.;
     146    int nc = 0;
     147    int nn = 0;
     148    float pt_ann[5];
     149
     150    for (int i = 0 ; i < 5 ; i++) {
     151      pt_ann[i] = 0.;
     152    }
     153
     154    if (fUseConstituents) {
    157155      TIter itConstituents(candidate->GetCandidates());
    158       while((constituent = static_cast<Candidate*>(itConstituents.Next())))
    159       {
    160         pt = constituent->Momentum.Pt();
    161         dr = candidate->Momentum.DeltaR(constituent->Momentum);
     156      while((constituent = static_cast<Candidate*>(itConstituents.Next()))) {
     157        float pt = constituent->Momentum.Pt();
     158        float dr = candidate->Momentum.DeltaR(constituent->Momentum);
    162159        sumpt += pt;
    163160        sumdrsqptsq += dr*dr*pt*pt;
    164161        sumptsq += pt*pt;
    165         if(constituent->Charge == 0)
    166         {
    167           // neutrals
    168           ++nn;
    169         }
    170         else
    171         {
    172           // charged
    173           if(constituent->IsPU && TMath::Abs(constituent->Position.Z()-zvtx) > fZVertexResolution)
    174           {
    175             sumptchpu += pt;
    176           }
    177           else
    178           {
    179             sumptchpv += pt;
    180           }
    181           sumptch += pt;
    182           ++nc;
    183         }
    184         for(i = 0; i < 5; ++i)
    185         {
    186           if(dr > 0.1*i && dr < 0.1*(i + 1))
    187           {
    188             pt_ann[i] += pt;
    189           }
    190         }
    191       }
    192     }
    193     else
    194     {
     162        if (constituent->Charge == 0) {
     163          // neutrals
     164          nn++;
     165        } else {
     166          // charged
     167          if (constituent->IsPU && TMath::Abs(constituent->Position.Z()-zvtx) > fZVertexResolution) {
     168            sumptchpu += pt;
     169          } else {
     170            sumptchpv += pt;
     171          }
     172          sumptch += pt;
     173          nc++;
     174        }
     175        for (int i = 0 ; i < 5 ; i++) {
     176          if (dr > 0.1*i && dr < 0.1*(i+1)) {
     177            pt_ann[i] += pt;
     178          }
     179        }
     180      }
     181    } else {
    195182      // Not using constituents, using dr
    196183      fItTrackInputArray->Reset();
    197       while((track = static_cast<Candidate*>(fItTrackInputArray->Next())))
    198       {
    199         if(track->Momentum.DeltaR(candidate->Momentum) < fParameterR)
    200         {
    201           pt = track->Momentum.Pt();
    202           sumpt += pt;
    203           sumptch += pt;
    204           if(track->IsPU && TMath::Abs(track->Position.Z()-zvtx) > fZVertexResolution)
    205           {
    206             sumptchpu += pt;
    207           }
    208           else
    209           {
    210             sumptchpv += pt;
    211           }
    212           dr = candidate->Momentum.DeltaR(track->Momentum);
    213           sumdrsqptsq += dr*dr*pt*pt;
    214           sumptsq += pt*pt;
    215           nc++;
    216           for(i = 0; i < 5; ++i)
    217           {
    218             if(dr > 0.1*i && dr < 0.1*(i + 1))
    219             {
     184       while ((trk = static_cast<Candidate*>(fItTrackInputArray->Next()))) {
     185        if (trk->Momentum.DeltaR(candidate->Momentum) < fParameterR) {
     186          float pt = trk->Momentum.Pt();
     187          sumpt += pt;
     188          sumptch += pt;
     189          if (trk->IsPU && TMath::Abs(trk->Position.Z()-zvtx) > fZVertexResolution) {
     190            sumptchpu += pt;
     191          } else {
     192            sumptchpv += pt;
     193          }
     194          float dr = candidate->Momentum.DeltaR(trk->Momentum);
     195          sumdrsqptsq += dr*dr*pt*pt;
     196          sumptsq += pt*pt;
     197          nc++;
     198          for (int i = 0 ; i < 5 ; i++) {
     199            if (dr > 0.1*i && dr < 0.1*(i+1)) {
    220200              pt_ann[i] += pt;
    221             }
    222           }
    223         }
    224       }
    225 
     201            }
     202          }
     203        }
     204      }
    226205      fItNeutralInputArray->Reset();
    227       while ((constituent = static_cast<Candidate*>(fItNeutralInputArray->Next())))
    228       {
    229         if(constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR)
    230         {
    231           pt = constituent->Momentum.Pt();
    232           sumpt += pt;
    233           dr = candidate->Momentum.DeltaR(constituent->Momentum);
    234           sumdrsqptsq += dr*dr*pt*pt;
    235           sumptsq += pt*pt;
    236           nn++;
    237           for(i = 0; i < 5; ++i)
    238           {
    239             if(dr > 0.1*i && dr < 0.1*(i + 1))
    240             {
    241               pt_ann[i] += pt;
    242             }
    243           }
    244         }
    245       }
    246     }
    247 
    248     if(sumptch > 0.0)
    249     {
     206      while ((constituent = static_cast<Candidate*>(fItNeutralInputArray->Next()))) {
     207        if (constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR) {
     208          float pt = constituent->Momentum.Pt();
     209          sumpt += pt;
     210          float dr = candidate->Momentum.DeltaR(constituent->Momentum);
     211          sumdrsqptsq += dr*dr*pt*pt;
     212          sumptsq += pt*pt;
     213          nn++;
     214          for (int i = 0 ; i < 5 ; i++) {
     215            if (dr > 0.1*i && dr < 0.1*(i+1)) {
     216              pt_ann[i] += pt;
     217            }
     218          }
     219        }
     220      }
     221    }
     222         
     223    if (sumptch > 0.) {
    250224      candidate->Beta = sumptchpu/sumptch;
    251225      candidate->BetaStar = sumptchpv/sumptch;
    252     }
    253     else
    254     {
    255       candidate->Beta = -999.0;
    256       candidate->BetaStar = -999.0;
    257     }
    258     if(sumptsq > 0.0)
    259     {
     226    } else {
     227      candidate->Beta = -999.;
     228      candidate->BetaStar = -999.;
     229    }
     230    if (sumptsq > 0.) {
    260231      candidate->MeanSqDeltaR = sumdrsqptsq/sumptsq;
    261     }
    262     else
    263     {
    264       candidate->MeanSqDeltaR = -999.0;
     232    } else {
     233      candidate->MeanSqDeltaR = -999.;
    265234    }
    266235    candidate->NCharged = nc;
    267236    candidate->NNeutrals = nn;
    268     if(sumpt > 0.0)
    269     {
     237    if (sumpt > 0.) {
    270238      candidate->PTD = TMath::Sqrt(sumptsq) / sumpt;
    271       for(i = 0; i < 5; ++i)
    272       {
     239      for (int i = 0 ; i < 5 ; i++) {
    273240        candidate->FracPt[i] = pt_ann[i]/sumpt;
    274241      }
    275     }
    276     else
    277     {
    278       candidate->PTD = -999.0;
    279       for(i = 0; i < 5; ++i)
    280       {
    281         candidate->FracPt[i] = -999.0;
     242    } else {
     243      candidate->PTD = -999.;
     244      for (int i = 0 ; i < 5 ; i++) {
     245        candidate->FracPt[i] = -999.;
    282246      }
    283247    }
Note: See TracChangeset for help on using the changeset viewer.