Fork me on GitHub

Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • modules/RunPUPPI.cc

    r341014c r509d1b4  
    11#include "modules/RunPUPPI.h"
    22
     3#include "PUPPI/RecoObj2.hh"
    34#include "PUPPI/AlgoObj.hh"
    45#include "PUPPI/PuppiContainer.hh"
    5 #include "PUPPI/RecoObj2.hh"
    66
    77#include "fastjet/PseudoJet.hh"
    … …  
    1111#include "classes/DelphesFormula.h"
    1212
    13 #include <algorithm>
     13#include <algorithm>
     14#include <stdexcept>
    1415#include <iostream>
    1516#include <sstream>
    16 #include <stdexcept>
    1717#include <vector>
    1818
    … …  
    2222//------------------------------------------------------------------------------
    2323RunPUPPI::RunPUPPI() :
    24   fItTrackInputArray(0),
     24  fItTrackInputArray(0), 
    2525  fItNeutralInputArray(0)
    26 {
    27 }
    28 
    29 //------------------------------------------------------------------------------
    30 RunPUPPI::~RunPUPPI() {}
    31 
    32 //------------------------------------------------------------------------------
    33 
    34 void RunPUPPI::Init()
    35 {
     26{}
     27
     28//------------------------------------------------------------------------------
     29RunPUPPI::~RunPUPPI(){}
     30
     31//------------------------------------------------------------------------------
     32
     33void RunPUPPI::Init(){
    3634  // input collection
    37   fTrackInputArray = ImportArray(GetString("TrackInputArray", "Calorimeter/towers"));
    38   fItTrackInputArray = fTrackInputArray->MakeIterator();
    39   fNeutralInputArray = ImportArray(GetString("NeutralInputArray", "Calorimeter/towers"));
     35  fTrackInputArray     = ImportArray(GetString("TrackInputArray", "Calorimeter/towers"));
     36  fItTrackInputArray   = fTrackInputArray->MakeIterator();
     37  fNeutralInputArray   = ImportArray(GetString("NeutralInputArray", "Calorimeter/towers"));
    4038  fItNeutralInputArray = fNeutralInputArray->MakeIterator();
    41   fPVInputArray = ImportArray(GetString("PVInputArray", "PV"));
    42   fPVItInputArray = fPVInputArray->MakeIterator();
    43   // puppi parameters
    44   fApplyNoLep = GetBool("UseNoLep", true);
     39  fPVInputArray        = ImportArray(GetString("PVInputArray", "PV"));
     40  fPVItInputArray      = fPVInputArray->MakeIterator();
     41  // puppi parameters                                 
     42  fApplyNoLep     = GetBool("UseNoLep", true);   
    4543  fMinPuppiWeight = GetDouble("MinPuppiWeight", 0.01);
    46   fUseExp = GetBool("UseExp", false);
     44  fUseExp         = GetBool("UseExp", false);
    4745  // read eta min ranges
    4846  ExRootConfParam param = GetParam("EtaMinBin");
    … …  
    5755  fPtMinBin.clear();
    5856  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fPtMinBin.push_back(param[iMap].GetDouble());
    59   // read cone size
     57  // read cone size                                                                                                                                                           
    6058  param = GetParam("ConeSizeBin");
    6159  fConeSizeBin.clear();
    6260  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fConeSizeBin.push_back(param[iMap].GetDouble());
    63   // read RMS min pt
     61  // read RMS min pt                                                                                                                                             
    6462  param = GetParam("RMSPtMinBin");
    6563  fRMSPtMinBin.clear();
    … …  
    7775  fNeutralPtSlope.clear();
    7876  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fNeutralPtSlope.push_back(param[iMap].GetDouble());
    79   // read apply chs
     77  // read apply chs                                                                                                                                                           
    8078  //param = GetParam("ApplyCHS");
    8179  //fApplyCHS.clear();
    … …  
    8987  fApplyLowPUCorr.clear();
    9088  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fApplyLowPUCorr.push_back(param[iMap].GetBool());
    91   // read metric id
     89  // read metric id                                                                                                                                                         
    9290  param = GetParam("MetricId");
    9391  fMetricId.clear();
    … …  
    9896  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fCombId.push_back(param[iMap].GetInt());
    9997  // create output array
    100   fOutputArray = ExportArray(GetString("OutputArray", "puppiParticles"));
    101   fOutputTrackArray = ExportArray(GetString("OutputArrayTracks", "puppiTracks"));
     98  fOutputArray        = ExportArray(GetString("OutputArray", "puppiParticles"));
     99  fOutputTrackArray   = ExportArray(GetString("OutputArrayTracks", "puppiTracks"));
    102100  fOutputNeutralArray = ExportArray(GetString("OutputArrayNeutrals", "puppiNeutrals"));
    103101  // Create algorithm list for puppi
    104102  std::vector<AlgoObj> puppiAlgo;
    105   if(puppiAlgo.empty())
    106   {
     103  if(puppiAlgo.empty()){
    107104    if(!(fEtaMinBin.size() == fEtaMaxBin.size() and fEtaMinBin.size() == fPtMinBin.size() and fEtaMinBin.size() == fConeSizeBin.size() and fEtaMinBin.size() == fRMSPtMinBin.size()
    108          and fEtaMinBin.size() == fRMSScaleFactorBin.size() and fEtaMinBin.size() == fNeutralMinEBin.size() and fEtaMinBin.size() == fNeutralPtSlope.size()
    109          and fEtaMinBin.size() == fUseCharged.size()
    110          and fEtaMinBin.size() == fApplyLowPUCorr.size() and fEtaMinBin.size() == fMetricId.size()))
    111     {
    112       std::cerr << " Error in PUPPI configuration, algo info should have the same size --> exit from the code" << std::endl;
     105         and fEtaMinBin.size() == fRMSScaleFactorBin.size() and fEtaMinBin.size() == fNeutralMinEBin.size() and  fEtaMinBin.size() == fNeutralPtSlope.size()
     106         and fEtaMinBin.size() == fUseCharged.size()
     107         and fEtaMinBin.size() == fApplyLowPUCorr.size() and fEtaMinBin.size() == fMetricId.size())) {
     108      std::cerr<<" Error in PUPPI configuration, algo info should have the same size --> exit from the code"<<std::endl;
    113109      std::exit(EXIT_FAILURE);
    114110    }
    115111  }
    116   for(size_t iAlgo = 0; iAlgo < fEtaMinBin.size(); iAlgo++)
    117   {
    118     AlgoObj algoTmp;
    119     algoTmp.etaMin = fEtaMinBin.at(iAlgo);
    120     algoTmp.etaMax = fEtaMaxBin.at(iAlgo);
    121     algoTmp.ptMin = fPtMinBin.at(iAlgo);
    122     algoTmp.minNeutralPt = fNeutralMinEBin.at(iAlgo);
     112  for( size_t iAlgo =  0 ; iAlgo < fEtaMinBin.size() ; iAlgo++){
     113    AlgoObj algoTmp ;
     114    algoTmp.etaMin            = fEtaMinBin.at(iAlgo);
     115    algoTmp.etaMax            = fEtaMaxBin.at(iAlgo);
     116    algoTmp.ptMin             = fPtMinBin.at(iAlgo);
     117    algoTmp.minNeutralPt      = fNeutralMinEBin.at(iAlgo);
    123118    algoTmp.minNeutralPtSlope = fNeutralPtSlope.at(iAlgo);
    124119    //Eta Extrapolation stuff is missing
    125120    //Loop through file requiring algos for same bins to be adjacent
    126     while(iAlgo < fEtaMinBin.size() and algoTmp.etaMin == fEtaMinBin.at(iAlgo) and algoTmp.etaMax == fEtaMaxBin.at(iAlgo))
    127     {
     121    while(iAlgo < fEtaMinBin.size() and algoTmp.etaMin == fEtaMinBin.at(iAlgo) and algoTmp.etaMax == fEtaMaxBin.at(iAlgo)) {
    128122      AlgoSubObj algoSubTmp;
    129       algoSubTmp.metricId = fMetricId.at(iAlgo);
    130       algoSubTmp.useCharged = fUseCharged.at(iAlgo);
    131       algoSubTmp.applyLowPUCorr = fApplyLowPUCorr.at(iAlgo);
    132       algoSubTmp.combId = fCombId.at(iAlgo);
    133       algoSubTmp.coneSize = fConeSizeBin.at(iAlgo);
    134       algoSubTmp.rmsPtMin = fRMSPtMinBin.at(iAlgo);
    135       algoSubTmp.rmsScaleFactor = fRMSScaleFactorBin.at(iAlgo);
     123      algoSubTmp.metricId          = fMetricId.at(iAlgo);
     124      algoSubTmp.useCharged        = fUseCharged.at(iAlgo);
     125      algoSubTmp.applyLowPUCorr    = fApplyLowPUCorr.at(iAlgo);
     126      algoSubTmp.combId            = fCombId.at(iAlgo);
     127      algoSubTmp.coneSize          = fConeSizeBin.at(iAlgo);
     128      algoSubTmp.rmsPtMin          = fRMSPtMinBin.at(iAlgo);
     129      algoSubTmp.rmsScaleFactor    = fRMSScaleFactorBin.at(iAlgo);
    136130      algoTmp.subAlgos.push_back(algoSubTmp);
    137131      iAlgo++;
    138132    }
    139133    iAlgo--;
    140     //if(std::find(puppiAlgo.begin(),puppiAlgo.end(),algoTmp) != puppiAlgo.end()) continue;
    141     puppiAlgo.push_back(algoTmp);
    142   }
    143   fPuppi = new PuppiContainer(true, fUseExp, fMinPuppiWeight, puppiAlgo);
     134    //if(std::find(puppiAlgo.begin(),puppiAlgo.end(),algoTmp) != puppiAlgo.end()) continue;   
     135    puppiAlgo.push_back(algoTmp);     
     136  }
     137  fPuppi  = new PuppiContainer(true,fUseExp,fMinPuppiWeight,puppiAlgo);
    144138}
    145139
    146140//------------------------------------------------------------------------------
    147141
    148 void RunPUPPI::Finish()
    149 {
    150   if(fItTrackInputArray) delete fItTrackInputArray;
     142void RunPUPPI::Finish(){
     143  if(fItTrackInputArray)   delete fItTrackInputArray;
    151144  if(fItNeutralInputArray) delete fItNeutralInputArray;
    152145  if(fPuppi) delete fPuppi;
    … …  
    155148//------------------------------------------------------------------------------
    156149
    157 void RunPUPPI::Process()
    158 {
     150void RunPUPPI::Process(){
    159151
    160152  Candidate *candidate, *particle;
    … …  
    164156
    165157  // loop over input objects
    166   fItTrackInputArray->Reset();
    167   fItNeutralInputArray->Reset();
    168   fPVItInputArray->Reset();
     158  fItTrackInputArray   ->Reset();
     159  fItNeutralInputArray ->Reset();
     160  fPVItInputArray      ->Reset();
    169161
    170162  std::vector<Candidate *> InputParticles;
    171163  InputParticles.clear();
    172164
    173   // take the leading vertex
     165  // take the leading vertex 
    174166  float PVZ = 0.;
    175   Candidate *pv = static_cast<Candidate *>(fPVItInputArray->Next());
    176   if(pv) PVZ = pv->Position.Z();
     167  Candidate *pv = static_cast<Candidate*>(fPVItInputArray->Next());
     168  if (pv) PVZ = pv->Position.Z();
    177169  // Fill input particles for puppi
    178170  std::vector<RecoObj> puppiInputVector;
    179171  puppiInputVector.clear();
    180   int lNBad = 0;
     172  int lNBad  = 0;
    181173  // Loop on charge track candidate
    182   while((candidate = static_cast<Candidate *>(fItTrackInputArray->Next())))
    183   {
    184     momentum = candidate->Momentum;
    185     RecoObj curRecoObj;
    186     curRecoObj.pt = momentum.Pt();
    187     curRecoObj.eta = momentum.Eta();
    188     curRecoObj.phi = momentum.Phi();
    189     curRecoObj.m = momentum.M();
    190     particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0)); //if(fApplyNoLep && TMath::Abs(candidate->PID) == 11) continue; //Dumb cut to minimize the nolepton on electron
    191     //if(fApplyNoLep && TMath::Abs(candidate->PID) == 13) continue;
    192     if(candidate->IsRecoPU and candidate->Charge != 0)
    193     { // if it comes fromPU vertexes after the resolution smearing and the dZ matching within resolution
    194       lNBad++;
    195       curRecoObj.id = 2;
    196       curRecoObj.vtxId = 0.7 * (fPVInputArray->GetEntries()); //Hack apply reco vtx efficiency of 70% for calibration
    197       if(TMath::Abs(candidate->PID) == 11)
    198         curRecoObj.pfType = 2;
    199       else if(TMath::Abs(candidate->PID) == 13)
    200         curRecoObj.pfType = 3;
    201       else if(TMath::Abs(candidate->PID) == 22)
    202         curRecoObj.pfType = 4;
    203       else
    204         curRecoObj.pfType = 1;
    205       curRecoObj.dZ = particle->Position.Z() - PVZ;
    206     }
    207     else if(!candidate->IsRecoPU && candidate->Charge != 0)
    208     {
    209       curRecoObj.id = 1; // charge from LV
    210       curRecoObj.vtxId = 1; // from PV
    211       if(TMath::Abs(candidate->PID) == 11)
    212         curRecoObj.pfType = 2;
    213       else if(TMath::Abs(candidate->PID) == 13)
    214         curRecoObj.pfType = 3;
    215       else if(TMath::Abs(candidate->PID) == 22)
    216         curRecoObj.pfType = 4;
    217       else
    218         curRecoObj.pfType = 1;
    219       curRecoObj.dZ = particle->Position.Z() - PVZ;
    220     }
    221     else
    222     {
    223       std::cerr << " RunPUPPI: problem with a charged track --> it has charge 0 " << std::endl;
    224       continue;
    225     }
    226 
    227     puppiInputVector.push_back(curRecoObj);
    228     InputParticles.push_back(candidate);
    229   }
    230 
    231   // Loop on neutral calo cells
    232   while((candidate = static_cast<Candidate *>(fItNeutralInputArray->Next())))
    233   {
    234     momentum = candidate->Momentum;
    235     RecoObj curRecoObj;
    236     curRecoObj.pt = momentum.Pt();
    237     curRecoObj.eta = momentum.Eta();
    238     curRecoObj.phi = momentum.Phi();
    239     curRecoObj.m = momentum.M();
    240     curRecoObj.charge = 0;
    241     particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
    242     if(candidate->Charge == 0)
    243     {
    244       curRecoObj.id = 0; // neutrals have id==0
    245       curRecoObj.vtxId = 0; // neutrals have vtxId==0
    246       if(TMath::Abs(candidate->PID) == 11)
    247         curRecoObj.pfType = 2;
    248       else if(TMath::Abs(candidate->PID) == 13)
    249         curRecoObj.pfType = 3;
    250       else if(TMath::Abs(candidate->PID) == 22)
    251         curRecoObj.pfType = 4;
    252       else
    253         curRecoObj.pfType = 5;
    254       curRecoObj.dZ = particle->Position.Z() - PVZ;
    255     }
    256     else
    257     {
    258       std::cerr << " RunPUPPI: problem with a neutrals cells --> it has charge !=0 " << std::endl;
    259       continue;
    260     }
    261     puppiInputVector.push_back(curRecoObj);
    262     InputParticles.push_back(candidate);
     174  while((candidate = static_cast<Candidate*>(fItTrackInputArray->Next()))){   
     175      momentum = candidate->Momentum;
     176      RecoObj curRecoObj;
     177      curRecoObj.pt  = momentum.Pt();
     178      curRecoObj.eta = momentum.Eta();
     179      curRecoObj.phi = momentum.Phi();
     180      curRecoObj.m   = momentum.M(); 
     181      particle = static_cast<Candidate*>(candidate->GetCandidates()->At(0));//if(fApplyNoLep && TMath::Abs(candidate->PID) == 11) continue; //Dumb cut to minimize the nolepton on electron
     182      //if(fApplyNoLep && TMath::Abs(candidate->PID) == 13) continue;
     183      if (candidate->IsRecoPU and candidate->Charge !=0) { // if it comes fromPU vertexes after the resolution smearing and the dZ matching within resolution
     184        lNBad++;
     185        curRecoObj.id    = 2;
     186        curRecoObj.vtxId = 0.7*(fPVInputArray->GetEntries()); //Hack apply reco vtx efficiency of 70% for calibration
     187        if(TMath::Abs(candidate->PID) == 11)      curRecoObj.pfType = 2;
     188        else if(TMath::Abs(candidate->PID) == 13) curRecoObj.pfType = 3;
     189        else if(TMath::Abs(candidate->PID) == 22) curRecoObj.pfType = 4;
     190        else curRecoObj.pfType = 1;
     191        curRecoObj.dZ = particle->Position.Z()-PVZ;
     192      }
     193      else if(!candidate->IsRecoPU && candidate->Charge !=0) {
     194        curRecoObj.id    = 1;  // charge from LV
     195        curRecoObj.vtxId = 1; // from PV
     196        if(TMath::Abs(candidate->PID) == 11)      curRecoObj.pfType = 2;
     197        else if(TMath::Abs(candidate->PID) == 13) curRecoObj.pfType = 3;
     198        else if(TMath::Abs(candidate->PID) == 22) curRecoObj.pfType = 4;
     199        else curRecoObj.pfType = 1;
     200        curRecoObj.dZ = particle->Position.Z()-PVZ;
     201      }
     202      else {
     203        std::cerr<<" RunPUPPI: problem with a charged track --> it has charge 0 "<<std::endl;
     204        continue;
     205      }
     206
     207      puppiInputVector.push_back(curRecoObj);
     208      InputParticles.push_back(candidate);
     209  }
     210
     211  // Loop on neutral calo cells
     212  while((candidate = static_cast<Candidate*>(fItNeutralInputArray->Next()))){
     213      momentum = candidate->Momentum;
     214      RecoObj curRecoObj;
     215      curRecoObj.pt  = momentum.Pt();
     216      curRecoObj.eta = momentum.Eta();
     217      curRecoObj.phi = momentum.Phi();
     218      curRecoObj.m   = momentum.M();
     219      curRecoObj.charge = 0;
     220      particle = static_cast<Candidate*>(candidate->GetCandidates()->At(0));
     221      if(candidate->Charge == 0){
     222        curRecoObj.id    = 0; // neutrals have id==0
     223        curRecoObj.vtxId = 0; // neutrals have vtxId==0
     224        if(TMath::Abs(candidate->PID) == 11)      curRecoObj.pfType = 2;
     225        else if(TMath::Abs(candidate->PID) == 13) curRecoObj.pfType = 3;
     226        else if(TMath::Abs(candidate->PID) == 22) curRecoObj.pfType = 4;
     227        else curRecoObj.pfType = 5;
     228        curRecoObj.dZ = particle->Position.Z()-PVZ;
     229      }
     230      else{
     231        std::cerr<<" RunPUPPI: problem with a neutrals cells --> it has charge !=0 "<<std::endl;
     232        continue;
     233      }
     234      puppiInputVector.push_back(curRecoObj);
     235      InputParticles.push_back(candidate);
    263236  }
    264237  // Create PUPPI container
    … …  
    268241
    269242  // Loop on final particles
    270   for(std::vector<PseudoJet>::iterator it = puppiParticles.begin(); it != puppiParticles.end(); it++)
    271   {
    272     if(it->user_index() <= int(InputParticles.size()))
    273     {
     243  for (std::vector<PseudoJet>::iterator it = puppiParticles.begin() ; it != puppiParticles.end() ; it++) {
     244    if(it->user_index() <= int(InputParticles.size())){     
    274245      candidate = static_cast<Candidate *>(InputParticles.at(it->user_index())->Clone());
    275       candidate->Momentum.SetPxPyPzE(it->px(), it->py(), it->pz(), it->e());
     246      candidate->Momentum.SetPxPyPzE(it->px(),it->py(),it->pz(),it->e());
    276247      fOutputArray->Add(candidate);
    277       if(puppiInputVector.at(it->user_index()).id == 1 or puppiInputVector.at(it->user_index()).id == 2)
    278         fOutputTrackArray->Add(candidate);
    279       else if(puppiInputVector.at(it->user_index()).id == 0)
    280         fOutputNeutralArray->Add(candidate);
     248      if( puppiInputVector.at(it->user_index()).id == 1 or puppiInputVector.at(it->user_index()).id == 2) fOutputTrackArray->Add(candidate);
     249      else if (puppiInputVector.at(it->user_index()).id == 0) fOutputNeutralArray->Add(candidate);
    281250    }
    282     else
    283     {
    284       std::cerr << " particle not found in the input Array --> skip " << std::endl;
     251    else{
     252      std::cerr<<" particle not found in the input Array --> skip "<<std::endl;
    285253      continue;
    286     }
    287   }
     254    }     
     255  }
     256
    288257}
Note: See TracChangeset for help on using the changeset viewer.