Fork me on GitHub

Changeset 4dd5358 in git for modules/RunPUPPI.cc


Ignore:
Timestamp:
Jun 5, 2016, 11:13:04 AM (8 years ago)
Author:
Michele Selvaggi <michele.selvaggi@…>
Branches:
ImprovedOutputFile, Timing, dual_readout, llp, master
Children:
355a7d7
Parents:
a9c67b3a (diff), 7d55e69a (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'master' into dev_01

File:
1 edited

Legend:

Unmodified
Added
Removed

  • modules/RunPUPPI.cc

    ra9c67b3a r4dd5358  
    11#include "modules/RunPUPPI.h"
    22
    3 #include "PUPPI/RecoObj2.hh"
    4 #include "PUPPI/AlgoObj.hh"
    5 //#include "PUPPI/puppiParticle.hh"
    6 //#include "PUPPI/puppiAlgoBin.hh"
     3#include "PUPPI/puppiCleanContainer.hh"
     4#include "PUPPI/RecoObj.hh"
     5#include "PUPPI/puppiParticle.hh"
     6#include "PUPPI/puppiAlgoBin.hh"
    77
    88#include "classes/DelphesClasses.h"
    … …  
    3030
    3131void RunPUPPI::Init(){
     32
    3233  // input collection
    3334  fTrackInputArray     = ImportArray(GetString("TrackInputArray", "Calorimeter/towers"));
    … …  
    3738  fPVInputArray        = ImportArray(GetString("PVInputArray", "PV"));
    3839  fPVItInputArray      = fPVInputArray->MakeIterator();
    39   // puppi parameters                                 
    40   fApplyNoLep     = GetBool("UseNoLep", true);   
     40
     41
     42  // puppi parameters                                     
    4143  fMinPuppiWeight = GetDouble("MinPuppiWeight", 0.01);
    4244  fUseExp         = GetBool("UseExp", false);
    43   // read eta min ranges
     45
     46  // read eta min ranges                                                                                                                                                           
    4447  ExRootConfParam param = GetParam("EtaMinBin");
    4548  fEtaMinBin.clear();
    4649  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fEtaMinBin.push_back(param[iMap].GetDouble());
    47   // read eta max ranges
     50
     51  // read eta max ranges                                                                                                                                                           
    4852  param = GetParam("EtaMaxBin");
    4953  fEtaMaxBin.clear();
    5054  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fEtaMaxBin.push_back(param[iMap].GetDouble());
    51   // read pt min value
     55
     56  // read pt min value                                                                                                                                                           
    5257  param = GetParam("PtMinBin");
    5358  fPtMinBin.clear();
    5459  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fPtMinBin.push_back(param[iMap].GetDouble());
     60
    5561  // read cone size                                                                                                                                                           
    5662  param = GetParam("ConeSizeBin");
    5763  fConeSizeBin.clear();
    5864  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fConeSizeBin.push_back(param[iMap].GetDouble());
     65
    5966  // read RMS min pt                                                                                                                                             
    6067  param = GetParam("RMSPtMinBin");
    6168  fRMSPtMinBin.clear();
    6269  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fRMSPtMinBin.push_back(param[iMap].GetDouble());
    63   // read RMS scale factor
     70
     71  // read RMS scale factor                                                                                                                                                           
    6472  param = GetParam("RMSScaleFactorBin");
    6573  fRMSScaleFactorBin.clear();
    6674  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fRMSScaleFactorBin.push_back(param[iMap].GetDouble());
    67   // read neutral pt min cut
     75
     76  // read neutral pt min cut                                                                                                                                                           
    6877  param = GetParam("NeutralMinEBin");
    6978  fNeutralMinEBin.clear();
    7079  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fNeutralMinEBin.push_back(param[iMap].GetDouble());
    71   // read neutral pt min slope
     80
     81  // read neutral pt min slope                                                                                                                                                           
    7282  param = GetParam("NeutralPtSlope");
    7383  fNeutralPtSlope.clear();
    7484  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fNeutralPtSlope.push_back(param[iMap].GetDouble());
     85
    7586  // read apply chs                                                                                                                                                           
    76   //param = GetParam("ApplyCHS");
    77   //fApplyCHS.clear();
    78   //for(int iMap = 0; iMap < param.GetSize(); ++iMap) fApplyCHS.push_back(param[iMap].GetBool());
    79   // read use charged
     87  param = GetParam("ApplyCHS");
     88  fApplyCHS.clear();
     89  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fApplyCHS.push_back(param[iMap].GetBool());
     90
     91  // read use charged                                                                                                                                                           
    8092  param = GetParam("UseCharged");
    8193  fUseCharged.clear();
    8294  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fUseCharged.push_back(param[iMap].GetBool());
    83   // read apply chs correction
     95
     96  // read apply chs correction                                                                                                                                                           
    8497  param = GetParam("ApplyLowPUCorr");
    8598  fApplyLowPUCorr.clear();
    8699  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fApplyLowPUCorr.push_back(param[iMap].GetBool());
     100 
    87101  // read metric id                                                                                                                                                         
    88102  param = GetParam("MetricId");
    89103  fMetricId.clear();
    90104  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fMetricId.push_back(param[iMap].GetInt());
    91   // scheme for combining
    92   param = GetParam("CombId");
    93   fCombId.clear();
    94   for(int iMap = 0; iMap < param.GetSize(); ++iMap) fCombId.push_back(param[iMap].GetInt());
     105
    95106  // create output array
    96107  fOutputArray        = ExportArray(GetString("OutputArray", "puppiParticles"));
    97108  fOutputTrackArray   = ExportArray(GetString("OutputArrayTracks", "puppiTracks"));
    98109  fOutputNeutralArray = ExportArray(GetString("OutputArrayNeutrals", "puppiNeutrals"));
    99   // Create algorithm list for puppi
    100   std::vector<AlgoObj> puppiAlgo;
    101   if(puppiAlgo.empty()){
    102     if(!(fEtaMinBin.size() == fEtaMaxBin.size() and fEtaMinBin.size() == fPtMinBin.size() and fEtaMinBin.size() == fConeSizeBin.size() and fEtaMinBin.size() == fRMSPtMinBin.size()
    103          and fEtaMinBin.size() == fRMSScaleFactorBin.size() and fEtaMinBin.size() == fNeutralMinEBin.size() and  fEtaMinBin.size() == fNeutralPtSlope.size()
    104          and fEtaMinBin.size() == fUseCharged.size()
    105          and fEtaMinBin.size() == fApplyLowPUCorr.size() and fEtaMinBin.size() == fMetricId.size())) {
    106       std::cerr<<" Error in PUPPI configuration, algo info should have the same size --> exit from the code"<<std::endl;
    107       std::exit(EXIT_FAILURE);
    108     }
    109   }
    110   for( size_t iAlgo =  0 ; iAlgo < fEtaMinBin.size() ; iAlgo++){
    111     AlgoObj algoTmp ;
    112     algoTmp.etaMin            = fEtaMinBin.at(iAlgo);
    113     algoTmp.etaMax            = fEtaMaxBin.at(iAlgo);
    114     algoTmp.ptMin             = fPtMinBin.at(iAlgo);
    115     algoTmp.minNeutralPt      = fNeutralMinEBin.at(iAlgo);
    116     algoTmp.minNeutralPtSlope = fNeutralPtSlope.at(iAlgo);
    117     //Eta Extrapolation stuff is missing
    118     //Loop through file requiring algos for same bins to be adjacent
    119     while(iAlgo < fEtaMinBin.size() and algoTmp.etaMin == fEtaMinBin.at(iAlgo) and algoTmp.etaMax == fEtaMaxBin.at(iAlgo)) {
    120       AlgoSubObj algoSubTmp;
    121       algoSubTmp.metricId          = fMetricId.at(iAlgo);
    122       algoSubTmp.useCharged        = fUseCharged.at(iAlgo);
    123       algoSubTmp.applyLowPUCorr    = fApplyLowPUCorr.at(iAlgo);
    124       algoSubTmp.combId            = fCombId.at(iAlgo);
    125       algoSubTmp.coneSize          = fConeSizeBin.at(iAlgo);
    126       algoSubTmp.rmsPtMin          = fRMSPtMinBin.at(iAlgo);
    127       algoSubTmp.rmsScaleFactor    = fRMSScaleFactorBin.at(iAlgo);
    128       algoTmp.subAlgos.push_back(algoSubTmp);
    129       iAlgo++;
    130     }
    131     iAlgo--;
    132     //if(std::find(puppiAlgo.begin(),puppiAlgo.end(),algoTmp) != puppiAlgo.end()) continue;   
    133     puppiAlgo.push_back(algoTmp);     
    134   }
    135   fPuppi  = new PuppiContainer(true,fUseExp,fMinPuppiWeight,puppiAlgo);
    136110}
    137111
    … …  
    150124  TLorentzVector momentum;
    151125
    152   //DelphesFactory *factory = GetFactory();
     126  DelphesFactory *factory = GetFactory();
    153127
    154128  // loop over input objects
    155   fItTrackInputArray   ->Reset();
    156   fItNeutralInputArray ->Reset();
    157   fPVItInputArray      ->Reset();
     129  fItTrackInputArray->Reset();
     130  fItNeutralInputArray->Reset();
     131  fPVItInputArray->Reset();
    158132
    159133  std::vector<Candidate *> InputParticles;
    … …  
    168142  std::vector<RecoObj> puppiInputVector;
    169143  puppiInputVector.clear();
    170   int lNBad  = 0;
     144
    171145  // Loop on charge track candidate
    172146  while((candidate = static_cast<Candidate*>(fItTrackInputArray->Next()))){   
     147
    173148      momentum = candidate->Momentum;
     149     
    174150      RecoObj curRecoObj;
    175151      curRecoObj.pt  = momentum.Pt();
    … …  
    178154      curRecoObj.m   = momentum.M(); 
    179155      particle = static_cast<Candidate*>(candidate->GetCandidates()->Last());
    180       //if(fApplyNoLep && TMath::Abs(candidate->PID) == 11) continue; //Dumb cut to minimize the nolepton on electron
    181       //if(fApplyNoLep && TMath::Abs(candidate->PID) == 13) continue;
    182156      if (candidate->IsRecoPU and candidate->Charge !=0) { // if it comes fromPU vertexes after the resolution smearing and the dZ matching within resolution
    183         lNBad++;
    184157        curRecoObj.id    = 2;
    185         curRecoObj.vtxId = 0.7*(fPVInputArray->GetEntries()); //Hack apply reco vtx efficiency of 70% for calibration
     158        curRecoObj.vtxId = candidate->IsPU;
    186159        if(TMath::Abs(candidate->PID) == 11)      curRecoObj.pfType = 2;
    187160        else if(TMath::Abs(candidate->PID) == 13) curRecoObj.pfType = 3;
    … …  
    210183  // Loop on neutral calo cells
    211184  while((candidate = static_cast<Candidate*>(fItNeutralInputArray->Next()))){
     185
    212186      momentum = candidate->Momentum;
     187
    213188      RecoObj curRecoObj;
    214189      curRecoObj.pt  = momentum.Pt();
    … …  
    216191      curRecoObj.phi = momentum.Phi();
    217192      curRecoObj.m   = momentum.M();
    218       curRecoObj.charge = 0;
    219193      particle = static_cast<Candidate*>(candidate->GetCandidates()->Last());
     194
    220195
    221196      if(candidate->Charge == 0){
    … …  
    235210      InputParticles.push_back(candidate);
    236211  }
     212
     213  // Create algorithm list for puppi
     214  std::vector<puppiAlgoBin> puppiAlgo;
     215  if(puppiAlgo.empty()){
     216   if(!(fEtaMinBin.size() == fEtaMaxBin.size() and fEtaMinBin.size() == fPtMinBin.size() and fEtaMinBin.size() == fConeSizeBin.size() and fEtaMinBin.size() == fRMSPtMinBin.size()
     217       and fEtaMinBin.size() == fRMSScaleFactorBin.size() and fEtaMinBin.size() == fNeutralMinEBin.size() and  fEtaMinBin.size() == fNeutralPtSlope.size()
     218       and fEtaMinBin.size() == fApplyCHS.size()  and fEtaMinBin.size() == fUseCharged.size()
     219       and fEtaMinBin.size() == fApplyLowPUCorr.size() and fEtaMinBin.size() == fMetricId.size())) {
     220    std::cerr<<" Error in PUPPI configuration, algo info should have the same size --> exit from the code"<<std::endl;
     221    std::exit(EXIT_FAILURE);
     222   }
     223
     224   for( size_t iAlgo =  0 ; iAlgo < fEtaMinBin.size() ; iAlgo++){
     225    puppiAlgoBin algoTmp ;
     226    algoTmp.fEtaMin_ = fEtaMinBin.at(iAlgo);
     227    algoTmp.fEtaMax_ = fEtaMaxBin.at(iAlgo);
     228    algoTmp.fPtMin_  = fPtMinBin.at(iAlgo);
     229    algoTmp.fConeSize_        = fConeSizeBin.at(iAlgo);
     230    algoTmp.fRMSPtMin_        = fRMSPtMinBin.at(iAlgo);
     231    algoTmp.fRMSScaleFactor_  = fRMSScaleFactorBin.at(iAlgo);
     232    algoTmp.fNeutralMinE_     = fNeutralMinEBin.at(iAlgo);
     233    algoTmp.fNeutralPtSlope_  = fNeutralPtSlope.at(iAlgo);
     234    algoTmp.fApplyCHS_        = fApplyCHS.at(iAlgo);
     235    algoTmp.fUseCharged_      = fUseCharged.at(iAlgo);
     236    algoTmp.fApplyLowPUCorr_  = fApplyLowPUCorr.at(iAlgo);
     237    algoTmp.fMetricId_        = fMetricId.at(iAlgo);
     238    if(std::find(puppiAlgo.begin(),puppiAlgo.end(),algoTmp) != puppiAlgo.end()) continue;   
     239    puppiAlgo.push_back(algoTmp);     
     240   }
     241  } 
     242
    237243  // Create PUPPI container
    238   fPuppi->initialize(puppiInputVector);
    239   fPuppi->puppiWeights();
    240   std::vector<fastjet::PseudoJet> puppiParticles = fPuppi->puppiParticles();
     244  puppiCleanContainer curEvent(puppiInputVector,puppiAlgo,fMinPuppiWeight,fUseExp);
     245  std::vector<fastjet::PseudoJet> puppiParticles = curEvent.puppiEvent();
    241246
    242247  // Loop on final particles
Note: See TracChangeset for help on using the changeset viewer.