Pythia8Interface/plugins/Py8toJetInput.cc

0001 #include "GeneratorInterface/Pythia8Interface/plugins/Py8toJetInput.h"
0002
0003 #include "SimDataFormats/GeneratorProducts/interface/LHECommonBlocks.h"
0004 #include "GeneratorInterface/LHEInterface/interface/LHERunInfo.h"
0005 #include "GeneratorInterface/LHEInterface/interface/LHEEvent.h"
0006
0007 const std::vector<fastjet::PseudoJet> Py8toJetInput::fillJetAlgoInput(const Event& event,
0008                                                                       const Event& workEvent,
0009                                                                       const lhef::LHEEvent* lhee,
0010                                                                       const std::vector<int>* typeIdx) {
0011   fJetInput.clear();
0012
0013   Event workEventJet = workEvent;
0014
0015   const lhef::HEPEUP& hepeup = *lhee->getHEPEUP();
0016
0017   std::set<int> typeSet[3];
0018
0019   // FIXME !!!
0020   // This is not safe to assume it's 3 because we're passing in a pointer
0021   // and we do NOT know what the actuial size is. I'll have to improve it.
0022   //
0023   for (size_t i = 0; i < 3; i++) {
0024     typeSet[i].clear();
0025     for (size_t j = 0; j < typeIdx[i].size(); j++) {
0026       // HEPEUP and Py8 event record are shifted by 3
0027       // (system particle + 2 beam particles in Py8 event)
0028       // ... EXCEPT FOR THE DECAY PRODUCTS IF DONE AT THE ME LEVEL !!!!!
0029       // ... in such case, the decay productes get gutted out of the sequence
0030       // and get placed in Py8 Event later on...
0031       // so we need to figure out the shift
0032       int pos = typeIdx[i][j];
0033       int shift = 3;
0034       // skip the first 2 entrirs in HEPEUP - they're incoming partons
0035       for (int ip = 2; ip < pos; ip++) {
0036         // alternative can be: moth1 != 1 && moth2 !=2...
0037         // but moth1 == moth2 means pointer to the same mother t
0038         // that can only be a resonance, unless moth1==moth2==0
0039         //
0040         if (hepeup.MOTHUP[ip].first == hepeup.MOTHUP[ip].second) {
0041           shift--;
0042         }
0043       }
0044       pos += shift;
0045       // typeSet[i].insert( event[pos].daughter1() );
0046       typeSet[i].insert(pos);
0047     }
0048   }
0049
0050   // int iTypeEnd = (typeIdx[1].empty()) ? 1 : 2;
0051
0052   // --> FIXME !!!
0053   int iType = 0;     // only LIGHT jets for now
0054   int jetAllow = 0;  // hardcoded for now for the same value as is set in Steve's example
0055                      // at present, not even in use...
0056                      // int jetMatch = 0; // hardcoded for now for the same value as is set in Steve's example
0057
0058   // Loop over particles and decide what to pass to the jet algorithm
0059   for (int i = 0; i < workEventJet.size(); ++i) {
0060     if (!workEventJet[i].isFinal())
0061       continue;
0062
0063     // jetAllow option to disallow certain particle types
0064     if (jetAllow == 1) {
0065       // Original AG+Py6 algorithm explicitly excludes tops,
0066       // leptons and photons.
0067       int id = workEventJet[i].idAbs();
0068       if ((id >= 11 && id <= 16) || id == ID_TOP || id == ID_PHOTON) {
0069         workEventJet[i].statusNeg();
0070         continue;
0071       }
0072     }
0073
0074     // Get the index of this particle in original event
0075     int idx = workEventJet[i].daughter1();
0076
0077     // Start with particle idx, and afterwards track mothers
0078     while (true) {
0079       // Light jets
0080       if (iType == 0) {
0081         // Do not include if originates from heavy jet or 'other'
0082         if (typeSet[1].find(idx) != typeSet[1].end() || typeSet[2].find(idx) != typeSet[2].end()) {
0083           workEventJet[i].statusNeg();
0084           break;
0085         }
0086
0087         // Made it to start of event record so done
0088         if (idx == 0) {
0089           break;
0090         }
0091         // Otherwise next mother and continue
0092         idx = event[idx].mother1();
0093
0094         // Heavy jets
0095       } else if (iType == 1) {
0096         // Only include if originates from heavy jet
0097         if (typeSet[1].find(idx) != typeSet[1].end())
0098           break;
0099
0100         // Made it to start of event record with no heavy jet mother,
0101         // so DO NOT include particle
0102         if (idx == 0) {
0103           workEventJet[i].statusNeg();
0104           break;
0105         }
0106
0107         // Otherwise next mother and continue
0108         idx = event[idx].mother1();
0109
0110       }  // if (iType)
0111     }    // while (true)
0112   }      // for (i)
0113
0114   // For jetMatch = 2, insert ghost particles corresponding to
0115   // each hard parton in the original process
0116   /*
0117   if (jetMatch > 0)
0118   {
0119
0120     for (int i = 0; i < int(typeIdx[iType].size()); i++)
0121     {
0122       // Get y/phi of the parton
0123       Vec4   pIn = eventProcess[typeIdx[iType][i]].p();
0124       double y   = Vec4y(pIn);
0125       double phi = pIn.phi();
0126
0127       // Create a ghost particle and add to the workEventJet
0128       double e   = MG5_GHOSTENERGY;
0129       double e2y = exp(2. * y);
0130       double pz  = e * (e2y - 1.) / (e2y + 1.);
0131       double pt  = sqrt(e*e - pz*pz);
0132       double px  = pt * cos(phi);
0133       double py  = pt * sin(phi);
0134       workEventJet.append(Particle(ID_GLUON, 99, 0, 0, 0, 0, 0, 0,
0135                                 px, py, pz, e, 0., 0, 9.));
0136
0137     } // for (i)
0138   } // if (jetMatch == 2)
0139 */
0140
0141   for (int i = 0; i < workEventJet.size(); i++) {
0142     // fisrt, weed out all entries marked with statusNeg();
0143     //
0144     if (workEventJet[i].status() < 0)
0145       continue;
0146
0147     // now weed out all entries above etaMax
0148     // in principle, we can use etaJetMaxAlgo because it's set equal to etaJetMax
0149     // as for etaJetMax, it gets set to memain_.etaclmax
0150     //
0151     if (fabs(workEventJet[i].eta()) > fJetEtaMax)
0152       continue;
0153
0154     // need to double check if native FastJet understands Py8 Event format
0155     // in general, PseudoJet gets formed from (px,py,pz,E)
0156     //
0157     fastjet::PseudoJet partTmp = workEventJet[i];
0158     fJetInput.push_back(partTmp);
0159   }
0160
0161   return fJetInput;
0162 }
0163
0164 int Py8toJetInput::getAncestor(int pos, const Event& fullEvent, const Event& workEvent) {
0165   int parentId = fullEvent[pos].mother1();
0166   int parentPrevId = 0;
0167   int counter = pos;
0168
0169   while (parentId > 0) {
0170     if (parentId == fullEvent[counter].mother2())  // carbon copy, keep walking up
0171     {
0172       parentPrevId = parentId;
0173       counter = parentId;
0174       parentId = fullEvent[parentPrevId].mother1();
0175       continue;
0176     }
0177
0178     // we get here if not a carbon copy
0179
0180     // let's check if it's a normal process, etc.
0181     //
0182     if ((parentId < parentPrevId) || parentId < fullEvent[counter].mother2())  // normal process
0183     {
0184       // first of all, check if hard block
0185       if (abs(fullEvent[counter].status()) == 22 || abs(fullEvent[counter].status()) == 23) {
0186         // yes, it's the hard block
0187         // we got what we want, and can exit now !
0188         parentId = counter;
0189         break;
0190       } else {
0191         parentPrevId = parentId;
0192         parentId = fullEvent[parentPrevId].mother1();
0193       }
0194     } else if (parentId > parentPrevId ||
0195                parentId > pos)  // "circular"/"forward-pointing parent" - intermediate process
0196     {
0197       parentId = -1;
0198       break;
0199     }
0200
0201     // additional checks... although we shouldn't be geeting here all that much...
0202     //
0203     if (abs(fullEvent[parentId].status()) == 22 || abs(fullEvent[parentId].status()) == 23)  // hard block
0204     {
0205       break;
0206     }
0207     if (abs(fullEvent[parentId].status()) < 22)  // incoming
0208     {
0209       parentId = -1;
0210       break;
0211     }
0212   }
0213
0214   return parentId;
0215 }
0216
0217 #include "HepMC/HEPEVT_Wrapper.h"
0218 #include <cassert>
0219
0220 const std::vector<fastjet::PseudoJet> Py8toJetInputHEPEVT::fillJetAlgoInput(const Event& event,
0221                                                                             const Event& workEvent,
0222                                                                             const lhef::LHEEvent* lhee,
0223                                                                             const std::vector<int>*) {
0224   fJetInput.clear();
0225
0226   HepMC::HEPEVT_Wrapper::zero_everything();
0227
0228   // service container for further mother-daughters links
0229   //
0230   std::vector<int> Py8PartonIdx;  // position of original (LHE) partons in Py8::Event
0231   Py8PartonIdx.clear();
0232   std::vector<int> HEPEVTPartonIdx;  // position of LHE partons in HEPEVT (incl. ME-generated decays)
0233   HEPEVTPartonIdx.clear();
0234
0235   // general counter
0236   //
0237   int index = 0;
0238
0239   int Py8PartonCounter = 0;
0240   int HEPEVTPartonCounter = 0;
0241
0242   // find the fisrt parton that comes from LHE (ME-generated)
0243   // skip the incoming particles/partons
0244   for (int iprt = 1; iprt < event.size(); iprt++) {
0245     const Particle& part = event[iprt];
0246     if (abs(part.status()) < 22)
0247       continue;  // below 10 is "service"
0248                  // 11-19 are beam particles; below 10 is "service"
0249                  // 21 is incoming partons
0250     Py8PartonCounter = iprt;
0251     break;
0252   }
0253
0254   const lhef::HEPEUP& hepeup = *lhee->getHEPEUP();
0255   // start the counter from 2, because we don't want the incoming particles/oartons !
0256   for (int iprt = 2; iprt < hepeup.NUP; iprt++) {
0257     index++;
0258     HepMC::HEPEVT_Wrapper::set_id(index, hepeup.IDUP[iprt]);
0259     HepMC::HEPEVT_Wrapper::set_status(index, 2);
0260     HepMC::HEPEVT_Wrapper::set_momentum(
0261         index, hepeup.PUP[iprt][0], hepeup.PUP[iprt][1], hepeup.PUP[iprt][2], hepeup.PUP[iprt][4]);
0262     HepMC::HEPEVT_Wrapper::set_mass(index, hepeup.PUP[iprt][4]);
0263     // --> FIXME HepMC::HEPEVT_Wrapper::set_position( index, part.xProd(), part.yProd(), part.zProd(), part.tProd() );
0264     HepMC::HEPEVT_Wrapper::set_parents(index, 0, 0);  // NO, not anymore to the "system particle"
0265     HepMC::HEPEVT_Wrapper::set_children(index, 0, 0);
0266     if (hepeup.MOTHUP[iprt].first > 2 &&
0267         hepeup.MOTHUP[iprt].second > 2)  // decay from LHE, will NOT show at the start of Py8 event !!!
0268     {
0269       HEPEVTPartonCounter++;
0270       continue;
0271     }
0272     Py8PartonIdx.push_back(Py8PartonCounter);
0273     Py8PartonCounter++;
0274     HEPEVTPartonIdx.push_back(HEPEVTPartonCounter);
0275     HEPEVTPartonCounter++;
0276   }
0277
0278   HepMC::HEPEVT_Wrapper::set_number_entries(index);
0279
0280   // now that the initial partons are in, attach parton-level from Pythia8
0281   // do NOT reset index as we need to *add* more particles sequentially
0282   //
0283   for (int iprt = 1; iprt < workEvent.size(); iprt++)  // from 0-entry (system) or from 1st entry ???
0284   {
0285     const Particle& part = workEvent[iprt];
0286
0287     //      if ( part.status() != 62 ) continue;
0288     if (part.status() < 51)
0289       continue;
0290     index++;
0291     HepMC::HEPEVT_Wrapper::set_id(index, part.id());
0292
0293     // HepMC::HEPEVT_Wrapper::set_status( index, event.statusHepMC(iprt) );
0294     HepMC::HEPEVT_Wrapper::set_status(index, 1);
0295     HepMC::HEPEVT_Wrapper::set_momentum(index, part.px(), part.py(), part.pz(), part.e());
0296     HepMC::HEPEVT_Wrapper::set_mass(index, part.m());
0297     HepMC::HEPEVT_Wrapper::set_position(index, part.xProd(), part.yProd(), part.zProd(), part.tProd());
0298     HepMC::HEPEVT_Wrapper::set_parents(index, 0, 0);  // just set to 0 like in Py6...
0299                                                       // although for some, mother will need to be re-set properly !
0300     HepMC::HEPEVT_Wrapper::set_children(index, 0, 0);
0301
0302     // now refine mother-daughters links, where applicable
0303
0304     int parentId = getAncestor(part.daughter1(), event, workEvent);
0305
0306     if (parentId <= 0)
0307       continue;
0308
0309     for (int idx = 0; idx < (int)Py8PartonIdx.size(); idx++) {
0310       if (parentId == Py8PartonIdx[idx]) {
0311         int idx1 = HEPEVTPartonIdx[idx];
0312         HepMC::HEPEVT_Wrapper::set_parents(index, idx1 + 1, idx1 + 1);
0313         break;
0314       }
0315     }
0316   }
0317
0318   HepMC::HEPEVT_Wrapper::set_number_entries(index);
0319
0320   // --> FIXME   HepMC::HEPEVT_Wrapper::set_event_number( fEventNumber ); // well, if you know it... well, it's one of the counters...
0321
0322   //   HepMC::HEPEVT_Wrapper::print_hepevt();
0323
0324   return fJetInput;
0325 }