AlpgenInterface/src/AlpgenEventRecordFixes.cc

0001 #include "GeneratorInterface/AlpgenInterface/interface/AlpgenEventRecordFixes.h"
0002 #include <cstdlib>
0003
0004 // An important peculiarity of this code is that, inside the
0005 // HEPEUP, particles are numbered from 1 (a la FORTRAN), but
0006 // the array that holds HEPEUP is acessed starting from 0
0007 // (a la C/C++). This has to be given due consideration.
0008 // For instance: to set a particle w/ index iup status code:
0009 //
0010 //     hepeup.ISTUP[iup] = 2
0011 //
0012 // but to set that particle as mother of a particle idau:
0013 //
0014 //     hepeup.MOTHUP[idau].first = iup+1
0015 //                                ^^^^^^^ See the difference?
0016 // Be careful with this.
0017
0018 math::XYZTLorentzVector alpgen::vectorFromHepeup(const lhef::HEPEUP &hepeup, int index) {
0019   return math::XYZTLorentzVector(
0020       hepeup.PUP[index][0], hepeup.PUP[index][1], hepeup.PUP[index][2], hepeup.PUP[index][3]);
0021 }
0022
0023 void alpgen::fixEventWZ(lhef::HEPEUP &hepeup) {
0024   // Nomenclature...
0025   int nup = hepeup.NUP;
0026
0027   // Open up space for the vector boson.
0028   hepeup.resize(nup + 1);
0029
0030   // Assignments specific to individual hard processes.
0031   // This one fixes the Event Record for W and Z.
0032   int iwch = 0;
0033
0034   // The last two particles in the record make the boson.
0035   double bosonPx = 0.;
0036   double bosonPy = 0.;
0037   double bosonPz = 0.;
0038   double bosonE = 0.;
0039   for (int i = nup - 2; i != nup; ++i) {
0040     bosonPx += hepeup.PUP[i][0];
0041     bosonPy += hepeup.PUP[i][1];
0042     bosonPz += hepeup.PUP[i][2];
0043     bosonE += hepeup.PUP[i][3];
0044     hepeup.MOTHUP[i].first = nup + 1;
0045     hepeup.MOTHUP[i].second = 0;
0046     iwch = (iwch - hepeup.IDUP[i] % 2);
0047   }
0048   // electron+nubar -> 11 + (-12) => -(1)+0 = -1  => W-
0049   // positron+nu    -> -11+ 12    => -(-1)+0 = +1 => W+
0050   // u dbar -> 2 -1  => 0 -(-1) = 1 => W+
0051   // c dbar -> 4 -1  => W+
0052   // etc.
0053
0054   // Boson ID.
0055   int bosonIndex = nup;
0056   int bosonId = 23;
0057   double bosonMass = std::sqrt(bosonE * bosonE - (bosonPx * bosonPx + bosonPy * bosonPy + bosonPz * bosonPz));
0058   if (iwch > 0)
0059     bosonId = 24;
0060   if (iwch < 0)
0061     bosonId = -24;
0062
0063   // Boson in the Event Record.
0064   hepeup.IDUP[bosonIndex] = bosonId;
0065   hepeup.ISTUP[bosonIndex] = 2;
0066   hepeup.MOTHUP[bosonIndex].first = 1;
0067   hepeup.MOTHUP[bosonIndex].second = 2;
0068   hepeup.PUP[bosonIndex][0] = bosonPx;
0069   hepeup.PUP[bosonIndex][1] = bosonPy;
0070   hepeup.PUP[bosonIndex][2] = bosonPz;
0071   hepeup.PUP[bosonIndex][3] = bosonE;
0072   hepeup.PUP[bosonIndex][4] = bosonMass;
0073   hepeup.ICOLUP[bosonIndex].first = 0;
0074   hepeup.ICOLUP[bosonIndex].second = 0;
0075
0076   hepeup.AQEDUP = hepeup.AQCDUP = -1.0;  // alphas are not saved by Alpgen
0077   for (int i = 0; i < hepeup.NUP; i++)
0078     hepeup.SPINUP[i] = -9;  // Alpgen does not store spin information
0079 }
0080
0081 void alpgen::fixEventMultiBoson(lhef::HEPEUP &hepeup) {
0082   int nup = hepeup.NUP;
0083
0084   // find first gauge bosons
0085   int ivstart = 0;
0086   int ivend = 0;
0087   for (int i = 0; i != nup; ++i) {
0088     if (std::abs(hepeup.IDUP[i]) == 24 || hepeup.IDUP[i] == 23) {
0089       hepeup.ISTUP[i] = 2;
0090       if (ivstart == 0)
0091         ivstart = i;
0092       ivend = i + 1;
0093     }
0094   }
0095   int nvb = ivend - ivstart;
0096
0097   // decay products pointers, starting from the end
0098   for (int i = 0; i != nvb; ++i) {
0099     hepeup.MOTHUP[nup - 2 * i - 1].first = ivend - i;
0100     hepeup.MOTHUP[nup - 2 * i - 2].first = ivend - i;
0101     hepeup.MOTHUP[nup - 2 * i - 1].second = 0;
0102     hepeup.MOTHUP[nup - 2 * i - 2].second = 0;
0103   }
0104
0105   hepeup.AQEDUP = hepeup.AQCDUP = -1.0;  // alphas are not saved by Alpgen
0106   for (int i = 0; i < hepeup.NUP; i++)
0107     hepeup.SPINUP[i] = -9;  // Alpgen does not store spin information
0108 }
0109
0110 void alpgen::fixEventTTbar(lhef::HEPEUP &hepeup) {
0111   using namespace math;
0112
0113   int nup = hepeup.NUP;
0114
0115   // Assert top is in the third position.
0116   int thirdID = hepeup.IDUP[2];
0117   if (std::abs(thirdID) != 6) {
0118     std::cout << "Top is NOT in the third position - no need to fix." << std::endl;
0119     return;
0120   }
0121
0122   // Open up space for 2 W bosons and two b quarks.
0123   hepeup.resize(nup + 4);
0124
0125   // reset top status codes
0126   hepeup.ISTUP[2] = 2;
0127   hepeup.ISTUP[3] = 2;
0128   int it;
0129   int itbar;
0130   if (thirdID == 6) {
0131     it = 2;
0132     itbar = 3;
0133   } else {
0134     it = 3;
0135     itbar = 2;
0136   }
0137
0138   // Reconstruct W's from decay product, fix mother-daughter relations.
0139   // Glossary: iwdec is the location of the first W decay product;
0140   //           iwup is the location where the W will be put.
0141   //           ibup is the location where the b will be put.
0142   for (int iw = 0; iw != 2; ++iw) {
0143     int iwdec = nup - 4 + 2 * iw;
0144     int iwup = nup + iw;
0145     int ibup = iwup + 2;
0146     int iwch = 0;
0147     for (int iup = iwdec; iup != iwdec + 2; ++iup) {
0148       hepeup.MOTHUP[iup].first = iwup + 1;
0149       hepeup.MOTHUP[iup].second = 0;
0150       iwch = (iwch - hepeup.IDUP[iup] % 2);
0151     }
0152     if (iwch > 0) {
0153       hepeup.IDUP[iwup] = 24;
0154       hepeup.IDUP[ibup] = 5;
0155       hepeup.MOTHUP[iwup].first = it + 1;
0156       hepeup.MOTHUP[iwup].second = 0;
0157       hepeup.MOTHUP[ibup].first = it + 1;
0158       hepeup.MOTHUP[ibup].second = 0;
0159     }
0160     if (iwch < 0) {
0161       hepeup.IDUP[iwup] = -24;
0162       hepeup.IDUP[ibup] = -5;
0163       hepeup.MOTHUP[iwup].first = itbar + 1;
0164       hepeup.MOTHUP[iwup].second = 0;
0165       hepeup.MOTHUP[ibup].first = itbar + 1;
0166       hepeup.MOTHUP[ibup].second = 0;
0167     }
0168     hepeup.ISTUP[iwup] = 2;
0169     hepeup.ISTUP[ibup] = 1;
0170
0171     // Reconstruct W momentum from its children.
0172     XYZTLorentzVector child1;
0173     child1 = vectorFromHepeup(hepeup, iwdec);
0174     XYZTLorentzVector child2;
0175     child2 = vectorFromHepeup(hepeup, iwdec + 1);
0176
0177     XYZTLorentzVector bosonW;
0178     bosonW = (child1 + child2);
0179     hepeup.PUP[iwup][0] = bosonW.Px();
0180     hepeup.PUP[iwup][1] = bosonW.Py();
0181     hepeup.PUP[iwup][2] = bosonW.Pz();
0182     hepeup.PUP[iwup][3] = bosonW.E();
0183     hepeup.PUP[iwup][4] = bosonW.M();
0184
0185     // Reconstruct b momentum from W and t.
0186     int topIndex = (hepeup.MOTHUP[iwup].first) - 1;
0187     XYZTLorentzVector topQuark;
0188     topQuark = vectorFromHepeup(hepeup, topIndex);
0189
0190     XYZTLorentzVector bottomQuark;
0191     bottomQuark = (topQuark - bosonW);
0192     hepeup.PUP[ibup][0] = bottomQuark.Px();
0193     hepeup.PUP[ibup][1] = bottomQuark.Py();
0194     hepeup.PUP[ibup][2] = bottomQuark.Pz();
0195     hepeup.PUP[ibup][3] = bottomQuark.E();
0196     hepeup.PUP[ibup][4] = bottomQuark.M();
0197
0198     // Set color labels.
0199     hepeup.ICOLUP[iwup].first = 0;
0200     hepeup.ICOLUP[iwup].second = 0;
0201     hepeup.ICOLUP[ibup].first = hepeup.ICOLUP[(hepeup.MOTHUP[iwup].first) - 1].first;
0202     hepeup.ICOLUP[ibup].second = hepeup.ICOLUP[(hepeup.MOTHUP[iwup].first) - 1].second;
0203   }
0204
0205   hepeup.AQEDUP = hepeup.AQCDUP = -1.0;  // alphas are not saved by Alpgen
0206   for (int i = 0; i < hepeup.NUP; i++)
0207     hepeup.SPINUP[i] = -9;  // Alpgen does not store spin information
0208 }
0209
0210 void alpgen::fixEventHiggsTTbar(lhef::HEPEUP &hepeup) {
0211   using namespace math;
0212
0213   int nup = hepeup.NUP;
0214
0215   // Assert top is in the fourth position.
0216   int fourthID = hepeup.IDUP[3];
0217   if (std::abs(fourthID) != 6) {
0218     std::cout << "Top is NOT in the fourth position - no need to fix." << std::endl;
0219     return;
0220   }
0221
0222   // Open up space for 2 W bosons and two b quarks.
0223   hepeup.resize(nup + 4);
0224
0225   // reset top status codes
0226   hepeup.ISTUP[3] = 2;
0227   hepeup.ISTUP[4] = 2;
0228   int it;
0229   int itbar;
0230   if (fourthID == 6) {
0231     it = 3;
0232     itbar = 4;
0233   } else {
0234     it = 4;
0235     itbar = 3;
0236   }
0237
0238   // Reconstruct W's from decay product, fix mother-daughter relations.
0239   // Glossary: iwdec is the location of the first W decay product;
0240   //           iwup is the location where the W will be put.
0241   //           ibup is the location where the b will be put.
0242   for (int iw = 0; iw != 2; ++iw) {
0243     int iwdec = nup - 4 + 2 * iw;
0244     int iwup = nup + iw;
0245     int ibup = iwup + 2;
0246     int iwch = 0;
0247     for (int iup = iwdec; iup != iwdec + 2; ++iup) {
0248       hepeup.MOTHUP[iup].first = iwup + 1;
0249       hepeup.MOTHUP[iup].second = 0;
0250       iwch = (iwch - hepeup.IDUP[iup] % 2);
0251     }
0252     if (iwch > 0) {
0253       hepeup.IDUP[iwup] = 24;
0254       hepeup.IDUP[ibup] = 5;
0255       hepeup.MOTHUP[iwup].first = it + 1;
0256       hepeup.MOTHUP[iwup].second = 0;
0257       hepeup.MOTHUP[ibup].first = it + 1;
0258       hepeup.MOTHUP[ibup].second = 0;
0259     }
0260     if (iwch < 0) {
0261       hepeup.IDUP[iwup] = -24;
0262       hepeup.IDUP[ibup] = -5;
0263       hepeup.MOTHUP[iwup].first = itbar + 1;
0264       hepeup.MOTHUP[iwup].second = 0;
0265       hepeup.MOTHUP[ibup].first = itbar + 1;
0266       hepeup.MOTHUP[ibup].second = 0;
0267     }
0268     hepeup.ISTUP[iwup] = 2;
0269     hepeup.ISTUP[ibup] = 1;
0270
0271     // Reconstruct W momentum from its children.
0272     XYZTLorentzVector child1;
0273     child1 = vectorFromHepeup(hepeup, iwdec);
0274     XYZTLorentzVector child2;
0275     child2 = vectorFromHepeup(hepeup, iwdec + 1);
0276
0277     XYZTLorentzVector bosonW;
0278     bosonW = (child1 + child2);
0279     hepeup.PUP[iwup][0] = bosonW.Px();
0280     hepeup.PUP[iwup][1] = bosonW.Py();
0281     hepeup.PUP[iwup][2] = bosonW.Pz();
0282     hepeup.PUP[iwup][3] = bosonW.E();
0283     hepeup.PUP[iwup][4] = bosonW.M();
0284
0285     // Reconstruct b momentum from W and t.
0286     int topIndex = (hepeup.MOTHUP[iwup].first) - 1;
0287     XYZTLorentzVector topQuark;
0288     topQuark = vectorFromHepeup(hepeup, topIndex);
0289
0290     XYZTLorentzVector bottomQuark;
0291     bottomQuark = (topQuark - bosonW);
0292     hepeup.PUP[ibup][0] = bottomQuark.Px();
0293     hepeup.PUP[ibup][1] = bottomQuark.Py();
0294     hepeup.PUP[ibup][2] = bottomQuark.Pz();
0295     hepeup.PUP[ibup][3] = bottomQuark.E();
0296     hepeup.PUP[ibup][4] = bottomQuark.M();
0297
0298     // Set color labels.
0299     hepeup.ICOLUP[iwup].first = 0;
0300     hepeup.ICOLUP[iwup].second = 0;
0301     hepeup.ICOLUP[ibup].first = hepeup.ICOLUP[(hepeup.MOTHUP[iwup].first) - 1].first;
0302     hepeup.ICOLUP[ibup].second = hepeup.ICOLUP[(hepeup.MOTHUP[iwup].first) - 1].second;
0303   }
0304
0305   hepeup.AQEDUP = hepeup.AQCDUP = -1.0;  // alphas are not saved by Alpgen
0306   for (int i = 0; i < hepeup.NUP; i++)
0307     hepeup.SPINUP[i] = -9;  // Alpgen does not store spin information
0308 }
0309
0310 void alpgen::fixEventSingleTop(lhef::HEPEUP &hepeup, double mb, int itopprc) {
0311   using namespace math;
0312
0313   int nup = hepeup.NUP;
0314
0315   // Number of W bosons.
0316   int nw = 1;
0317   // Fix number of Ws in final state.
0318   if (itopprc >= 3)
0319     nw = 2;
0320
0321   // Open up space for W bosons and b quarks.
0322   hepeup.resize(nup + 2);
0323
0324   // Assign mass to the incoming bottom quark, if required.
0325   for (int i = 0; i != 2; ++i) {
0326     if (std::abs(hepeup.IDUP[i]) == 5) {
0327       hepeup.PUP[i][4] = mb;
0328       double energyb = hepeup.PUP[i][3];
0329       hepeup.PUP[i][3] = std::sqrt(energyb * energyb + mb * mb);
0330     }
0331   }
0332
0333   // Identify the top.
0334   hepeup.ISTUP[2] = 2;
0335   int it = 0;
0336   int itbar = 0;
0337   if (hepeup.IDUP[2] == 6)
0338     it = 2;
0339   else if (hepeup.IDUP[2] == -6)
0340     itbar = 2;
0341   else {
0342     std::cout << "Wrong assumption about top position, stop." << std::endl;
0343     // FIXME: Should throw an exception here.
0344     return;
0345   }
0346
0347   // TOP DECAY
0348   // Reconstruct W's from decay products.
0349   // iwdec: 1st W decay product.
0350   int iwdec = 0;
0351   if (nw == 1)
0352     iwdec = nup - 2;
0353   else if (nw == 2)
0354     iwdec = nup - 4;
0355
0356   // Put W and b at the end.
0357   int iwup = nup;
0358   int ibup = iwup + 1;
0359   int iwch = 0;
0360   for (int iup = iwdec; iup != iwdec + 2; ++iup) {
0361     hepeup.MOTHUP[iup].first = iwup + 1;
0362     hepeup.MOTHUP[iup].second = 0;
0363     iwch = (iwch - hepeup.IDUP[iup] % 2);
0364   }
0365
0366   if (iwch > 0) {
0367     hepeup.IDUP[iwup] = 24;
0368     hepeup.IDUP[ibup] = 5;
0369     hepeup.MOTHUP[iwup].first = it + 1;
0370     hepeup.MOTHUP[iwup].second = 0;
0371     hepeup.MOTHUP[ibup].first = it + 1;
0372     hepeup.MOTHUP[ibup].second = 0;
0373   }
0374   if (iwch < 0) {
0375     hepeup.IDUP[iwup] = -24;
0376     hepeup.IDUP[ibup] = -5;
0377     hepeup.MOTHUP[iwup].first = itbar + 1;
0378     hepeup.MOTHUP[iwup].second = 0;
0379     hepeup.MOTHUP[ibup].first = itbar + 1;
0380     hepeup.MOTHUP[ibup].second = 0;
0381   }
0382   hepeup.ISTUP[iwup] = 2;
0383   hepeup.ISTUP[ibup] = 1;
0384
0385   // Reconstruct W momentum from its children.
0386   XYZTLorentzVector child1;
0387   child1 = vectorFromHepeup(hepeup, iwdec);
0388   XYZTLorentzVector child2;
0389   child2 = vectorFromHepeup(hepeup, iwdec + 1);
0390
0391   XYZTLorentzVector bosonW;
0392   bosonW = (child1 + child2);
0393   hepeup.PUP[iwup][0] = bosonW.Px();
0394   hepeup.PUP[iwup][1] = bosonW.Py();
0395   hepeup.PUP[iwup][2] = bosonW.Pz();
0396   hepeup.PUP[iwup][3] = bosonW.E();
0397   hepeup.PUP[iwup][4] = bosonW.M();
0398
0399   // Reconstruct b momentum from W and t.
0400   int topIndex = (hepeup.MOTHUP[iwup].first) - 1;
0401   XYZTLorentzVector topQuark;
0402   topQuark = vectorFromHepeup(hepeup, topIndex);
0403
0404   XYZTLorentzVector bottomQuark;
0405   bottomQuark = (topQuark - bosonW);
0406   hepeup.PUP[ibup][0] = bottomQuark.Px();
0407   hepeup.PUP[ibup][1] = bottomQuark.Py();
0408   hepeup.PUP[ibup][2] = bottomQuark.Pz();
0409   hepeup.PUP[ibup][3] = bottomQuark.E();
0410   hepeup.PUP[ibup][4] = bottomQuark.M();
0411
0412   // Set color labels.
0413   hepeup.ICOLUP[iwup].first = 0;
0414   hepeup.ICOLUP[iwup].second = 0;
0415   hepeup.ICOLUP[ibup].first = hepeup.ICOLUP[(hepeup.MOTHUP[iwup].first) - 1].first;
0416   hepeup.ICOLUP[ibup].second = hepeup.ICOLUP[(hepeup.MOTHUP[iwup].first) - 1].second;
0417
0418   nup = nup + 2;
0419
0420   // If there is a second W, we deal with it now.
0421   // Unlike the first W, this one appears EXPLICITLY in the Event Record.
0422   if (nw == 2) {
0423     // W DECAY
0424     // iwdec: 1st W decay product.
0425     iwdec = nup - 4;
0426     // iwup: location of the W in the event record.
0427     iwup = nup - 7;
0428     iwch = 0;
0429     for (int iup = iwdec; iup != iwdec + 2; ++iup) {
0430       hepeup.MOTHUP[iup].first = iwup + 1;
0431       hepeup.MOTHUP[iup].second = 0;
0432       iwch = (iwch - hepeup.IDUP[iup] % 2);
0433     }
0434     hepeup.ISTUP[iwup] = 2;
0435     hepeup.ICOLUP[iwup].first = 0;
0436     hepeup.ICOLUP[iwup].second = 0;
0437   }
0438 }