Fork me on GitHub

Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • modules/RunPUPPI.cc

    r509d1b4 r341014c  
    11#include "modules/RunPUPPI.h"
    22
    3 #include "PUPPI/RecoObj2.hh"
    43#include "PUPPI/AlgoObj.hh"
    54#include "PUPPI/PuppiContainer.hh"
     5#include "PUPPI/RecoObj2.hh"
    66
    77#include "fastjet/PseudoJet.hh"
    … …  
    1111#include "classes/DelphesFormula.h"
    1212
    13 #include <algorithm>
    14 #include <stdexcept>
     13#include <algorithm>
    1514#include <iostream>
    1615#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 RunPUPPI::~RunPUPPI(){}
    30 
    31 //------------------------------------------------------------------------------
    32 
    33 void RunPUPPI::Init(){
     26{
     27}
     28
     29//------------------------------------------------------------------------------
     30RunPUPPI::~RunPUPPI() {}
     31
     32//------------------------------------------------------------------------------
     33
     34void RunPUPPI::Init()
     35{
    3436  // input collection
    35   fTrackInputArray     = ImportArray(GetString("TrackInputArray", "Calorimeter/towers"));
    36   fItTrackInputArray   = fTrackInputArray->MakeIterator();
    37   fNeutralInputArray   = ImportArray(GetString("NeutralInputArray", "Calorimeter/towers"));
     37  fTrackInputArray = ImportArray(GetString("TrackInputArray", "Calorimeter/towers"));
     38  fItTrackInputArray = fTrackInputArray->MakeIterator();
     39  fNeutralInputArray = ImportArray(GetString("NeutralInputArray", "Calorimeter/towers"));
    3840  fItNeutralInputArray = fNeutralInputArray->MakeIterator();
    39   fPVInputArray        = ImportArray(GetString("PVInputArray", "PV"));
    40   fPVItInputArray      = fPVInputArray->MakeIterator();
    41   // puppi parameters                                 
    42   fApplyNoLep     = GetBool("UseNoLep", true);   
     41  fPVInputArray = ImportArray(GetString("PVInputArray", "PV"));
     42  fPVItInputArray = fPVInputArray->MakeIterator();
     43  // puppi parameters
     44  fApplyNoLep = GetBool("UseNoLep", true);
    4345  fMinPuppiWeight = GetDouble("MinPuppiWeight", 0.01);
    44   fUseExp         = GetBool("UseExp", false);
     46  fUseExp = GetBool("UseExp", false);
    4547  // read eta min ranges
    4648  ExRootConfParam param = GetParam("EtaMinBin");
    … …  
    5557  fPtMinBin.clear();
    5658  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fPtMinBin.push_back(param[iMap].GetDouble());
    57   // read cone size                                                                                                                                                           
     59  // read cone size
    5860  param = GetParam("ConeSizeBin");
    5961  fConeSizeBin.clear();
    6062  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fConeSizeBin.push_back(param[iMap].GetDouble());
    61   // read RMS min pt                                                                                                                                             
     63  // read RMS min pt
    6264  param = GetParam("RMSPtMinBin");
    6365  fRMSPtMinBin.clear();
    … …  
    7577  fNeutralPtSlope.clear();
    7678  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fNeutralPtSlope.push_back(param[iMap].GetDouble());
    77   // read apply chs                                                                                                                                                           
     79  // read apply chs
    7880  //param = GetParam("ApplyCHS");
    7981  //fApplyCHS.clear();
    … …  
    8789  fApplyLowPUCorr.clear();
    8890  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fApplyLowPUCorr.push_back(param[iMap].GetBool());
    89   // read metric id                                                                                                                                                         
     91  // read metric id
    9092  param = GetParam("MetricId");
    9193  fMetricId.clear();
    … …  
    9698  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fCombId.push_back(param[iMap].GetInt());
    9799  // create output array
    98   fOutputArray        = ExportArray(GetString("OutputArray", "puppiParticles"));
    99   fOutputTrackArray   = ExportArray(GetString("OutputArrayTracks", "puppiTracks"));
     100  fOutputArray = ExportArray(GetString("OutputArray", "puppiParticles"));
     101  fOutputTrackArray = ExportArray(GetString("OutputArrayTracks", "puppiTracks"));
    100102  fOutputNeutralArray = ExportArray(GetString("OutputArrayNeutrals", "puppiNeutrals"));
    101103  // Create algorithm list for puppi
    102104  std::vector<AlgoObj> puppiAlgo;
    103   if(puppiAlgo.empty()){
     105  if(puppiAlgo.empty())
     106  {
    104107    if(!(fEtaMinBin.size() == fEtaMaxBin.size() and fEtaMinBin.size() == fPtMinBin.size() and fEtaMinBin.size() == fConeSizeBin.size() and fEtaMinBin.size() == fRMSPtMinBin.size()
    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;
     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;
    109113      std::exit(EXIT_FAILURE);
    110114    }
    111115  }
    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);
     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);
    118123    algoTmp.minNeutralPtSlope = fNeutralPtSlope.at(iAlgo);
    119124    //Eta Extrapolation stuff is missing
    120125    //Loop through file requiring algos for same bins to be adjacent
    121     while(iAlgo < fEtaMinBin.size() and algoTmp.etaMin == fEtaMinBin.at(iAlgo) and algoTmp.etaMax == fEtaMaxBin.at(iAlgo)) {
     126    while(iAlgo < fEtaMinBin.size() and algoTmp.etaMin == fEtaMinBin.at(iAlgo) and algoTmp.etaMax == fEtaMaxBin.at(iAlgo))
     127    {
    122128      AlgoSubObj algoSubTmp;
    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);
     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);
    130136      algoTmp.subAlgos.push_back(algoSubTmp);
    131137      iAlgo++;
    132138    }
    133139    iAlgo--;
    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);
    138 }
    139 
    140 //------------------------------------------------------------------------------
    141 
    142 void RunPUPPI::Finish(){
    143   if(fItTrackInputArray)   delete fItTrackInputArray;
     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);
     144}
     145
     146//------------------------------------------------------------------------------
     147
     148void RunPUPPI::Finish()
     149{
     150  if(fItTrackInputArray) delete fItTrackInputArray;
    144151  if(fItNeutralInputArray) delete fItNeutralInputArray;
    145152  if(fPuppi) delete fPuppi;
    … …  
    148155//------------------------------------------------------------------------------
    149156
    150 void RunPUPPI::Process(){
     157void RunPUPPI::Process()
     158{
    151159
    152160  Candidate *candidate, *particle;
    … …  
    156164
    157165  // loop over input objects
    158   fItTrackInputArray   ->Reset();
    159   fItNeutralInputArray ->Reset();
    160   fPVItInputArray      ->Reset();
     166  fItTrackInputArray->Reset();
     167  fItNeutralInputArray->Reset();
     168  fPVItInputArray->Reset();
    161169
    162170  std::vector<Candidate *> InputParticles;
    163171  InputParticles.clear();
    164172
    165   // take the leading vertex 
     173  // take the leading vertex
    166174  float PVZ = 0.;
    167   Candidate *pv = static_cast<Candidate*>(fPVItInputArray->Next());
    168   if (pv) PVZ = pv->Position.Z();
     175  Candidate *pv = static_cast<Candidate *>(fPVItInputArray->Next());
     176  if(pv) PVZ = pv->Position.Z();
    169177  // Fill input particles for puppi
    170178  std::vector<RecoObj> puppiInputVector;
    171179  puppiInputVector.clear();
    172   int lNBad  = 0;
     180  int lNBad = 0;
    173181  // Loop on charge track 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);
     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);
    236263  }
    237264  // Create PUPPI container
    … …  
    241268
    242269  // Loop on final particles
    243   for (std::vector<PseudoJet>::iterator it = puppiParticles.begin() ; it != puppiParticles.end() ; it++) {
    244     if(it->user_index() <= int(InputParticles.size())){     
     270  for(std::vector<PseudoJet>::iterator it = puppiParticles.begin(); it != puppiParticles.end(); it++)
     271  {
     272    if(it->user_index() <= int(InputParticles.size()))
     273    {
    245274      candidate = static_cast<Candidate *>(InputParticles.at(it->user_index())->Clone());
    246       candidate->Momentum.SetPxPyPzE(it->px(),it->py(),it->pz(),it->e());
     275      candidate->Momentum.SetPxPyPzE(it->px(), it->py(), it->pz(), it->e());
    247276      fOutputArray->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);
    250     }
    251     else{
    252       std::cerr<<" particle not found in the input Array --> skip "<<std::endl;
     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);
     281    }
     282    else
     283    {
     284      std::cerr << " particle not found in the input Array --> skip " << std::endl;
    253285      continue;
    254     }     
    255   }
    256 
    257 }
     286    }
     287  }
     288}
Note: See TracChangeset for help on using the changeset viewer.