Ticket #39: partonic-decay-vertex-position.patch

SHERPA/Single_Events/Hadron_Decays.C

@@ -314 +314 @@
 
           return true;
         }
-        DEBUG_INFO("did not find exclusive decay channel, continue.");
+        else {
+          DEBUG_INFO("did not find exclusive decay channel, continue.");
+          Vec4D vertex_position=decayblob->Position();
+          showerblob->SetPosition(vertex_position);
+          fragmentationblob->SetPosition(vertex_position);
+        }
       }
     }
   }
@@ -344 +349 @@
 void Hadron_Decays::SetPosition(ATOOLS::Blob* blob)
 {
   Particle* inpart = blob->InParticle(0);
-  if(inpart->Flav().Kfcode()==kf_K) return;
+  if(inpart->Flav().Kfcode()==kf_K) {
+    blob->SetPosition(inpart->XProd());
+    return;
+  }
   
   // boost lifetime into lab
   double gamma = 1./rpa.gen.Accu();

SHERPA/SoftPhysics/Fragmentation_Handler.C

@@ -116 +116 @@
 {
   Particle  * part(NULL);
   Blob      * blob(NULL);
-  Part_List * plist = new Part_List;
+  std::vector<Part_List*> plists;
+  plists.push_back(new Part_List);
   for (Blob_List::iterator blit=bloblist->begin();blit!=bloblist->end();++blit) {
     if ((*blit)->Has(blob_status::needs_hadronization)) {
+      // If not coming from hadron decays, fill default plists[0].
+      // If from hadron decays, create separate plist for each blob
+      // such that there is a one-to-one correspondence between
+      // fragmentation outcome and hadron decay. This is needed for setting
+      // the correct vertex position, and to reject exclusive final states
+      Part_List* plist(plists[0]);
+      Blob* upstream_blob=(*blit)->UpstreamBlob();
+      if (upstream_blob && upstream_blob->Type()==btp::Hadron_Decay) {
+        plist=new Part_List;
+        plists.push_back(plist);
+      }
+
       std::vector<Particle*> taus;
       for (int i=0;i<(*blit)->NOutP();i++) {
         part = (*blit)->OutParticle(i); 
@@ -149 +162 @@
       (*blit)->UnsetStatus(blob_status::needs_hadronization);
     }
   }
-  if (plist->empty()) {
+  if (plists[0]->empty() && plists.size()<2) {
     msg_Debugging()<<"WARNING in Lund_Interface::PrepareFragmentationBlob:"<<endl
                    <<"   No coloured particle found leaving shower blobs."<<endl;
-    delete plist;
+    delete plists[0];
     return Return_Value::Nothing;
   }
 
 
+  Return_Value::code ret(Return_Value::Success);
+  for (size_t i=0; i<plists.size(); ++i) {
+  if (plists[i]->empty()) continue;
+  Part_List* plist=plists[i];
   int  col1, col2;
   bool hit1, hit2;
   Part_List * pli(NULL);
@@ -233 +250 @@
       delete partlists[i];
     }
     delete plist;
-    return Return_Value::Success;
+    ret=Return_Value::Success;
+  }
+  else {
+    for(size_t i=0;i<partlists.size();i++) {
+      delete partlists[i];
+    }
+    delete plist;
+    ret=Return_Value::Error;
+    break;
   }
-  for(size_t i=0;i<partlists.size();i++) {
-    delete partlists[i];
   }
-  delete plist;
-  return Return_Value::Error;
+  return ret;
 }