Fork me on GitHub

Changeset 0afc0fc in git for modules/RunPUPPI.cc


Ignore:
Timestamp:
Jun 1, 2016, 7:41:02 PM (8 years ago)
Author:
Philip <violatingcp@…>
Branches:
ImprovedOutputFile, Timing, dual_readout, llp, master
Children:
7b45ff5, ee0763a
Parents:
0a836f2
Message:

added the new Puppi

File:
1 edited

Legend:

Unmodified
Added
Removed

  • modules/RunPUPPI.cc

    r0a836f2 r0afc0fc  
    11#include "modules/RunPUPPI.h"
    22
    3 #include "PUPPI/puppiCleanContainer.hh"
    4 #include "PUPPI/RecoObj.hh"
    5 #include "PUPPI/puppiParticle.hh"
    6 #include "PUPPI/puppiAlgoBin.hh"
     3#include "PUPPI/RecoObj2.hh"
     4#include "PUPPI/AlgoObj.hh"
     5//#include "PUPPI/puppiParticle.hh"
     6//#include "PUPPI/puppiAlgoBin.hh"
    77
    88#include "classes/DelphesClasses.h"
    … …  
    3030
    3131void RunPUPPI::Init(){
    32 
    3332  // input collection
    3433  fTrackInputArray     = ImportArray(GetString("TrackInputArray", "Calorimeter/towers"));
    … …  
    3837  fPVInputArray        = ImportArray(GetString("PVInputArray", "PV"));
    3938  fPVItInputArray      = fPVInputArray->MakeIterator();
    40 
    41 
    42   // puppi parameters                                     
     39  // puppi parameters                                 
     40  fApplyNoLep     = GetBool("UseNoLep", true);   
    4341  fMinPuppiWeight = GetDouble("MinPuppiWeight", 0.01);
    4442  fUseExp         = GetBool("UseExp", false);
    45 
    46   // read eta min ranges                                                                                                                                                           
     43  // read eta min ranges
    4744  ExRootConfParam param = GetParam("EtaMinBin");
    4845  fEtaMinBin.clear();
    4946  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fEtaMinBin.push_back(param[iMap].GetDouble());
    50 
    51   // read eta max ranges                                                                                                                                                           
     47  // read eta max ranges
    5248  param = GetParam("EtaMaxBin");
    5349  fEtaMaxBin.clear();
    5450  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fEtaMaxBin.push_back(param[iMap].GetDouble());
    55 
    56   // read pt min value                                                                                                                                                           
     51  // read pt min value
    5752  param = GetParam("PtMinBin");
    5853  fPtMinBin.clear();
    5954  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fPtMinBin.push_back(param[iMap].GetDouble());
    60 
    6155  // read cone size                                                                                                                                                           
    6256  param = GetParam("ConeSizeBin");
    6357  fConeSizeBin.clear();
    6458  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fConeSizeBin.push_back(param[iMap].GetDouble());
    65 
    6659  // read RMS min pt                                                                                                                                             
    6760  param = GetParam("RMSPtMinBin");
    6861  fRMSPtMinBin.clear();
    6962  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fRMSPtMinBin.push_back(param[iMap].GetDouble());
    70 
    71   // read RMS scale factor                                                                                                                                                           
     63  // read RMS scale factor
    7264  param = GetParam("RMSScaleFactorBin");
    7365  fRMSScaleFactorBin.clear();
    7466  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fRMSScaleFactorBin.push_back(param[iMap].GetDouble());
    75 
    76   // read neutral pt min cut                                                                                                                                                           
     67  // read neutral pt min cut
    7768  param = GetParam("NeutralMinEBin");
    7869  fNeutralMinEBin.clear();
    7970  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fNeutralMinEBin.push_back(param[iMap].GetDouble());
    80 
    81   // read neutral pt min slope                                                                                                                                                           
     71  // read neutral pt min slope
    8272  param = GetParam("NeutralPtSlope");
    8373  fNeutralPtSlope.clear();
    8474  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fNeutralPtSlope.push_back(param[iMap].GetDouble());
    85 
    8675  // read apply chs                                                                                                                                                           
    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                                                                                                                                                           
     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
    9280  param = GetParam("UseCharged");
    9381  fUseCharged.clear();
    9482  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fUseCharged.push_back(param[iMap].GetBool());
    95 
    96   // read apply chs correction                                                                                                                                                           
     83  // read apply chs correction
    9784  param = GetParam("ApplyLowPUCorr");
    9885  fApplyLowPUCorr.clear();
    9986  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fApplyLowPUCorr.push_back(param[iMap].GetBool());
    100  
    10187  // read metric id                                                                                                                                                         
    10288  param = GetParam("MetricId");
    10389  fMetricId.clear();
    10490  for(int iMap = 0; iMap < param.GetSize(); ++iMap) fMetricId.push_back(param[iMap].GetInt());
    105 
     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());
    10695  // create output array
    10796  fOutputArray        = ExportArray(GetString("OutputArray", "puppiParticles"));
    10897  fOutputTrackArray   = ExportArray(GetString("OutputArrayTracks", "puppiTracks"));
    10998  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);
    110136}
    111137
    … …  
    124150  TLorentzVector momentum;
    125151
    126   DelphesFactory *factory = GetFactory();
     152  //DelphesFactory *factory = GetFactory();
    127153
    128154  // loop over input objects
    129   fItTrackInputArray->Reset();
    130   fItNeutralInputArray->Reset();
    131   fPVItInputArray->Reset();
     155  fItTrackInputArray   ->Reset();
     156  fItNeutralInputArray ->Reset();
     157  fPVItInputArray      ->Reset();
    132158
    133159  std::vector<Candidate *> InputParticles;
    … …  
    142168  std::vector<RecoObj> puppiInputVector;
    143169  puppiInputVector.clear();
    144 
     170  int lNBad  = 0;
    145171  // Loop on charge track candidate
    146172  while((candidate = static_cast<Candidate*>(fItTrackInputArray->Next()))){   
    147 
    148173      momentum = candidate->Momentum;
    149      
    150174      RecoObj curRecoObj;
    151175      curRecoObj.pt  = momentum.Pt();
    … …  
    154178      curRecoObj.m   = momentum.M(); 
    155179      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;
    156182      if (candidate->IsRecoPU and candidate->Charge !=0) { // if it comes fromPU vertexes after the resolution smearing and the dZ matching within resolution
     183        lNBad++;
    157184        curRecoObj.id    = 2;
    158         curRecoObj.vtxId = candidate->IsPU;
     185        curRecoObj.vtxId = 0.7*(fPVInputArray->GetEntries()); //Hack apply reco vtx efficiency of 70% for calibration
    159186        if(TMath::Abs(candidate->PID) == 11)      curRecoObj.pfType = 2;
    160187        else if(TMath::Abs(candidate->PID) == 13) curRecoObj.pfType = 3;
    … …  
    183210  // Loop on neutral calo cells
    184211  while((candidate = static_cast<Candidate*>(fItNeutralInputArray->Next()))){
    185 
    186212      momentum = candidate->Momentum;
    187 
    188213      RecoObj curRecoObj;
    189214      curRecoObj.pt  = momentum.Pt();
    … …  
    191216      curRecoObj.phi = momentum.Phi();
    192217      curRecoObj.m   = momentum.M();
     218      curRecoObj.charge = 0;
    193219      particle = static_cast<Candidate*>(candidate->GetCandidates()->Last());
    194 
    195220
    196221      if(candidate->Charge == 0){
    … …  
    210235      InputParticles.push_back(candidate);
    211236  }
    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 
    243237  // Create PUPPI container
    244   puppiCleanContainer curEvent(puppiInputVector,puppiAlgo,fMinPuppiWeight,fUseExp);
    245   std::vector<fastjet::PseudoJet> puppiParticles = curEvent.puppiEvent();
     238  fPuppi->initialize(puppiInputVector);
     239  fPuppi->puppiWeights();
     240  std::vector<fastjet::PseudoJet> puppiParticles = fPuppi->puppiParticles();
    246241
    247242  // Loop on final particles
Note: See TracChangeset for help on using the changeset viewer.