add optional extra PFO collections in delphes2lcio

examples/cpp/delphes2lcio/README.md

@@ -216,7 +216,7 @@ The following maps are required:
 
 - *MCParticleMap, PFOMap, JetMap, MuonMap, ElectronMap, PhotonMap*.
 
-These maps for extra jet collections are optional:
+These maps for extra jet collections are optional (comment out from the configuration file, if not present for your delphes card) :
 
 - *ExtraJetMap2, ExtraJetMap3, ExtraJetMap4, ExtraJetMap5, ExtraJetMap6*.
 
@@ -225,6 +225,21 @@ all other names. For jet collections the parameter `useDelphes4Vec` defines weth
 jet (!=0) or wether it is computed from the jet constituent PFOs (default).
 
 
+Additional PFO maps can also be added. Their names have to start with `"ExtraPFO"`. For example
+this will create an extra PFO collection called `BCalPFOs` from the Delphes branch `BCalPhoton`
+and assign the photon PDG code:
+
+```
+ExtraPFOMapBCal
+branchName  BCalPhoton
+isCharged  0
+lcioName  BCalPFOs
+mass  0.
+pdg  22
+
+```
+
+
 
 #### Event Meta data
 

examples/cpp/delphes2lcio/examples/delphes2lcio.cfg

@@ -66,6 +66,12 @@ branchName  Photon
 lcioName  IsolatedPhotons
 pdg  22
 # --------------------------------
+ExtraPFOMapBCal
+branchName  BCalPhoton
+isCharged  0
+lcioName  BCalPFOs
+mass  0.
+pdg  22
 
 #--------------------------------------------------------------------
 # here one can set additional parameters copied to the event header

examples/cpp/delphes2lcio/include/DelphesLCIOConfig.h

@@ -56,6 +56,9 @@ class DelphesLCIOConfig{
   /// return list of map names for extra jet collection (containing substring 'ExtraJet' )
   std::vector<std::string> getExtraJetMapNames() const ;
 
+  /// return list of map names for extra jet collection (containing substring 'ExtraPFO' )
+  std::vector<std::string> getExtraPFOMapNames() const ;
+
   /// convert string to int
   int toInt(const std::string& val) const ;
 
@@ -140,8 +143,17 @@ class DelphesLCIOConfig{
     { "ExtraJetMap3" , { { "lcioName", "Durham3Jets" }, { "branchName", "Jet_N3" }, { "useDelphes4Vec" , "0" }, { "storeYMerge" , "0" } } },
     { "ExtraJetMap4" , { { "lcioName", "Durham4Jets" }, { "branchName", "Jet_N4" }, { "useDelphes4Vec" , "0" }, { "storeYMerge" , "0" } } },
     { "ExtraJetMap5" , { { "lcioName", "Durham5Jets" }, { "branchName", "Jet_N5" }, { "useDelphes4Vec" , "0" }, { "storeYMerge" , "0" } } },
-    { "ExtraJetMap6" , { { "lcioName", "Durham6Jets" }, { "branchName", "Jet_N6" }, { "useDelphes4Vec" , "0" }, { "storeYMerge" , "0" } } }
+    { "ExtraJetMap6" , { { "lcioName", "Durham6Jets" }, { "branchName", "Jet_N6" }, { "useDelphes4Vec" , "0" }, { "storeYMerge" , "0" } } },
 
+    { "ExtraPFOMapBCal" ,
+      {
+	{ "lcioName"          , "BCalPFOs" },
+	{ "branchName"        , "BCalPhoton" },
+	{ "pdg"               , "22" },
+	{ "mass"              , "0." },
+	{ "isCharged"         , "0" },
+      }
+    },
 
   } ;
 

examples/cpp/delphes2lcio/include/DelphesLCIOConverter.h

@@ -70,6 +70,18 @@ class DelphesLCIOConverter {
    */
   bool convertJetCollection( TClonesArray* tca, EVENT::LCCollection* col, int useDelphes4Vec, int storeYMerge ) ;
 
+  /** Helper function to convert extra charged (Track) particles to ReconstructedParticles (PFOs).
+   *  The new elements are added to the LCIO collection. They will not be linked to the MCParticles as is done
+   *  for the standard PFOs.
+   */
+  int convertExtraPFOsCharged(  TClonesArray* tca, EVENT::LCCollection* col, int pdg, double mass );
+
+  /** Helper function to convert extra neutral (Tower) particles to ReconstructedParticles (PFOs).
+   *  The new elements are added to the LCIO collection. They will not be linked to the MCParticles as is done
+   *  for the standard PFOs.
+   */
+  int convertExtraPFOsNeutral(  TClonesArray* tca, EVENT::LCCollection* col, int pdg, double mass );
+
 private:
   IO::LCWriter *_writer=nullptr;
   IMPL::LCCollectionVec *_evtSumCol=nullptr;

examples/cpp/delphes2lcio/src/DelphesLCIOConfig.cc

@@ -11,13 +11,26 @@ std::vector<std::string> DelphesLCIOConfig::getExtraJetMapNames() const {
   std::vector<std::string> ejmapNames ;
 
   for( auto& it : _maps ){
-      
+
     if( ( it.first.find( "ExtraJet" ) != std::string::npos) )
       ejmapNames.push_back( it.first ) ;
   }
   return ejmapNames ;
 }
 
+/// return list of map names for extra PFO collection (containing substring 'ExtraPFO' )
+std::vector<std::string> DelphesLCIOConfig::getExtraPFOMapNames() const {
+
+  std::vector<std::string> ejmapNames ;
+
+  for( auto& it : _maps ){
+
+    if( ( it.first.find( "ExtraPFO" ) != std::string::npos) )
+      ejmapNames.push_back( it.first ) ;
+  }
+  return ejmapNames ;
+}
+
 // ---------------------------------------
 /// convert string to int
 int DelphesLCIOConfig::toInt(const std::string& val) const {

examples/cpp/delphes2lcio/src/DelphesLCIOConverter.cc

@@ -168,6 +168,7 @@ void DelphesLCIOConverter::convertTree2LCIO( TTree *tree , lcio::LCEventImpl* ev
   evt->addCollection( photons, _cfg->getPhotonParameter("lcioName") ) ;
 
   lcio::PIDHandler pfopidH( pfos );
+
   int showerParamID = pfopidH.addAlgorithm( "ShowerParameters" , {"emFraction","hadFraction"} ) ;
   int trackParamID = pfopidH.addAlgorithm( "TrackParameters" , {"L","PT","D0","DZ","Phi","CtgTheta","ErrorPT","ErrorD0",
 								"ErrorDZ","ErrorPhi","ErrorCtgTheta"} ) ;
@@ -532,6 +533,38 @@ void DelphesLCIOConverter::convertTree2LCIO( TTree *tree , lcio::LCEventImpl* ev
     }
   }
 //======================================================================
+  // ----  convert extra pfo collections if requested:
+
+  const auto& epnames = _cfg->getExtraPFOMapNames() ;
+
+  for( auto& pName : epnames ){
+
+    auto* pcol = new lcio::LCCollectionVec( lcio::LCIO::RECONSTRUCTEDPARTICLE )  ;
+
+    evt->addCollection( pcol,  _cfg->getMapParameter("lcioName", pName ) ) ;
+
+    br = tree->GetBranch( _cfg->getMapParameter("branchName", pName ).c_str() ) ;
+
+    int pdg = _cfg->toInt( _cfg->getMapParameter("pdg", pName ) ) ;
+
+    int isCharged = _cfg->toInt( _cfg->getMapParameter("isCharged", pName ) ) ;
+
+    double mass = _cfg->toFloat( _cfg->getMapParameter("mass", pName ) ) ;
+
+    if( br != nullptr ){
+
+      TClonesArray* tca = *(TClonesArray**) br->GetAddress()  ;
+
+//      std::cout << " convertExtraPFOsCharged(" << pName <<", "  << pdg << ", " << mass << ")" << std::endl ;
+
+      if( isCharged )
+	convertExtraPFOsCharged( tca, pcol , pdg, mass ) ;
+      else
+	convertExtraPFOsNeutral( tca, pcol , pdg, mass ) ;
+    }
+  }
+//======================================================================
+
   br = tree->GetBranch( _cfg->getPhotonParameter("branchName").c_str() ) ;
 
   pdg = _cfg->toInt( _cfg->getPhotonParameter("pdg") ) ;
@@ -838,3 +871,96 @@ bool DelphesLCIOConverter::convertJetCollection(TClonesArray* tca, EVENT::LCColl
 }
 
 //======================================================================
+
+int DelphesLCIOConverter::convertExtraPFOsCharged(  TClonesArray* tca, EVENT::LCCollection* pfos,
+						    int pdgCh, double massCh ) {
+
+  int chPFONum = 0;
+
+  lcio::PIDHandler pfopidH( pfos );
+  int trackParamID = pfopidH.addAlgorithm( "TrackParameters" , {"L","PT","D0","DZ","Phi","CtgTheta","ErrorPT","ErrorD0",
+								"ErrorDZ","ErrorPhi","ErrorCtgTheta"} ) ;
+
+  int nnh = tca->GetEntriesFast();
+
+  for(int j = 0; j < nnh ; ++j){
+    Track* trk = static_cast<Track*> (tca->At(j));
+
+    auto* pfo = new lcio::ReconstructedParticleImpl ;
+    pfos->addElement( pfo ) ;
+    ++chPFONum;
+
+    pfo->setType( pdgCh * trk->Charge );
+
+    pfo->setMass( massCh ) ;
+
+    double p = trk->P ;
+    double e = sqrt( p*p + massCh*massCh ) ;
+    double th = 2.*atan( exp( - trk->Eta ) );
+    double ph = trk->Phi ;
+    double m[3] = { p * cos( ph ) * sin( th ) , p * sin( ph ) * sin( th ) ,   p * cos( th ) } ;
+
+    pfo->setEnergy( e ) ;
+    pfo->setMomentum( m ) ;
+//      pfo->setCovMatrix (const float *cov) ; //fixme
+    pfo->setCharge( trk->Charge ) ;
+
+    pfopidH.setParticleID( pfo, 0, trk->PID , 1.,
+			   trackParamID, { trk->L,trk->PT,trk->D0,trk->DZ,trk->Phi,trk->CtgTheta,
+					   trk->ErrorPT, trk->ErrorD0,trk->ErrorDZ,trk->ErrorPhi,
+					   trk->ErrorCtgTheta  } ) ;
+
+    pfopidH.setParticleIDUsed( pfo, trackParamID );
+
+    pfo->setGoodnessOfPID( 1. ) ;
+
+    float ref[3] = { trk->X,  trk->Y,  trk->Z  } ;
+    pfo->setReferencePoint ( ref ) ;
+
+  }
+  return chPFONum ;
+}
+//======================================================================
+int DelphesLCIOConverter::convertExtraPFOsNeutral(  TClonesArray* tca, EVENT::LCCollection* pfos,
+						    int pdgNH, double massNH ) {
+
+  int neuPFONum = 0;
+
+  lcio::PIDHandler pfopidH( pfos );
+  int showerParamID = pfopidH.addAlgorithm( "ShowerParameters" , {"emFraction","hadFraction"} ) ;
+
+  int nnh = tca->GetEntriesFast();
+
+  for(int j = 0; j < nnh ; ++j){
+
+    Tower* p = static_cast<Tower*> (tca->At(j));
+
+    auto* pfo = new lcio::ReconstructedParticleImpl ;
+    pfos->addElement( pfo ) ;
+    ++neuPFONum ;
+
+    pfo->setType( pdgNH );
+
+    double e = p->E ;
+    double th = 2.*atan( exp( - p->Eta ) );
+    double ph = p->Phi ;
+    double pp = sqrt( e * e - massNH * massNH ) ;
+
+    double m[3] = { pp * cos( ph ) * sin( th ) , pp * sin( ph ) * sin( th ) , pp * cos( th ) } ;
+
+    pfo->setEnergy( e ) ;
+    pfo->setMomentum( m ) ;
+    //pfo->setCovMatrix (const float *cov) ; //fixme
+    pfo->setMass( massNH ) ;
+    pfo->setCharge(0.) ;
+
+    pfopidH.setParticleID( pfo, 0, 130, 1. , showerParamID, { p->Eem  / p->E , p->Ehad / p->E   } ) ;
+    pfopidH.setParticleIDUsed( pfo, showerParamID );
+
+    pfo->setGoodnessOfPID( 1. ) ;
+  }
+
+
+  return neuPFONum ;
+}
+//======================================================================