Ticket #59: aliasdecays.patch

HADRONS++/Main/Hadrons.C

@@ -35 +35 @@
   msg_Tracking()<<"In Hadrons: ("<<path<<") "<<constfile<<std::endl;
   p_decaymap = new Hadron_Decay_Map(path, file, constfile);
   p_decaymap->ReadInConstants();
+  p_decaymap->ReadHadronAliases(path, "HadronAliases.dat");
   p_decaymap->SetHadronProperties();
-  p_decaymap->Read();
+  p_decaymap->Read(path, file);
+  p_decaymap->Read(path, "AliasDecays.dat");
   p_decaymap->Initialise();
   p_decaymap->ReadFixedTables();
   p_mixinghandler = new Mixing_Handler(p_decaymap);

HADRONS++/Main/Hadron_Decay_Map.H

@@ -39 +39 @@
               std::string file,
               std::string constfile );
     ~Hadron_Decay_Map();
-    void Read();
+    void Read(const std::string& path, const std::string& file);
+    void ReadHadronAliases(const std::string& path, const std::string& file);
     void Initialise();
     void ReadFixedTables();
     void ReadInConstants();

HADRONS++/Main/Tools.H

@@ -24 +24 @@
   typedef std::map<std::string,double> String_Map;
 
   struct GeneralModel: public String_Map {
+    std::map<kf_code, kf_code> m_aliases;
+    
   public:
     inline double operator()(const std::string &tag,const double &def) const 
     {

HADRONS++/Main/Hadron_Decay_Channel.C

@@ -27 +27 @@
   if(p_amplitudes) delete p_amplitudes; p_amplitudes=NULL;
   if(p_momenta)    delete [] p_momenta; p_momenta=NULL;
   if(p_flavours)   delete [] p_flavours; p_flavours=NULL;
+  if(p_physicalflavours) delete [] p_physicalflavours; p_physicalflavours=NULL;
   for(size_t i=0;i<m_mes.size();i++)      delete m_mes[i].second;
   for(size_t i=0;i<m_currents.size();i++) {
     delete m_currents[i].second.first;
@@ -65 +66 @@
     msg_Error()<<"  Total outgoing mass heavier than incoming particle. Will abort."<<endl;
     abort();
   }
-  p_amplitudes = new Spin_Amplitudes(Flavours(),GetN(),Complex(0.0,0.0));
+  p_physicalflavours=new Flavour[GetN()];
+  for (int i=0; i<GetN(); ++i) {
+    map<kf_code,kf_code>::const_iterator it =
+        startmd.m_aliases.find(p_flavours[i].Kfcode());
+    if (it!=startmd.m_aliases.end())
+      p_physicalflavours[i] = Flavour(it->second, p_flavours[i].IsAnti());
+    else
+      p_physicalflavours[i] = p_flavours[i];
+  }
+  p_amplitudes = new Spin_Amplitudes(p_physicalflavours,GetN(),Complex(0.0,0.0));
   p_ps = new HD_PS_Base(this);
   // check for identical particles
   Flavour refflav;
   double symfactor (1);         
   int l(0), lfac (1);              
   for( int i=0; i<NOut(); ++i ) {
-    refflav = p_flavours[i+1];
+    refflav = p_physicalflavours[i+1];
     l = 0;
     lfac = 1;
     for( int j=0; j<NOut(); ++j ) {
-      if( p_flavours[j+1]==refflav ) {
+      if( p_physicalflavours[j+1]==refflav ) {
         l++;
         lfac *= l;
       }
@@ -147 +157 @@
     for(int i=0;i<NOut()+1;i++) {
       decayindices[i]=i;
     }
-    m_mes.push_back(MEPair(1.0,new Generic(p_flavours,NOut()+1,
+    m_mes.push_back(MEPair(1.0,new Generic(p_physicalflavours,NOut()+1,
                                            decayindices,"Generic")));
     delete [] decayindices;
     p_ps->AddChannel( string("Isotropic"), 1., m_startmd );
@@ -266 +276 @@
     for(int i=0;i<NOut()+1;i++) {
       decayindices[i]=i;
     }
-    m_mes.push_back(MEPair(1.0,new Generic(p_flavours,NOut()+1,
+    m_mes.push_back(MEPair(1.0,new Generic(p_physicalflavours,NOut()+1,
                                            decayindices,"Generic")));
     delete [] decayindices;
   }
@@ -498 +508 @@
 // differential with random PS points; just for weight
 double Hadron_Decay_Channel::Differential()
 {
-  p_momenta[0] = Vec4D(p_flavours[0].HadMass(),0.,0.,0.);        // decay from rest
+  p_momenta[0] = Vec4D(p_physicalflavours[0].HadMass(),0.,0.,0.);        // decay from rest
   p_ps->GeneratePoint(p_momenta,(PHASIC::Cut_Data*)(NULL));     // generate a PS point
   p_ps->GenerateWeight(p_momenta,(PHASIC::Cut_Data*)(NULL));    // calculate its weight factor
   double weight = p_ps->Weight();                               // get weight factor
@@ -539 +549 @@
   fi.n=resultstrings.size()-1;
   int* indices = new int[fi.n]; // will get deleted by Current_Base
   for(int i=0; i<fi.n; i++) indices[i] = ToType<int>(resultstrings[i+1]);
-  fi.flavs=p_flavours; fi.indices=indices;
+  fi.flavs=p_physicalflavours; fi.indices=indices;
 
   Current_Base* current = Current_Getter_Function::GetObject(resultstrings[0],fi);
   if(current==NULL) {
@@ -573 +583 @@
   fi.n=NOut()+1;
   int* indices = new int[fi.n]; // will get deleted by HD_ME_Base
   for(int i=0; i<fi.n; i++) indices[i] = ToType<int>(resultstrings[i+1]);
-  fi.flavs=p_flavours; fi.indices=indices;
+  fi.flavs=p_physicalflavours; fi.indices=indices;
 
   HD_ME_Base* me = HD_ME_Getter_Function::GetObject(resultstrings[0],fi);
   if(me==NULL) {

HADRONS++/Main/Hadron_Decay_Channel.H

@@ -27 +27 @@
   class Hadron_Decay_Channel : public ATOOLS::Decay_Channel {
   private:
     ATOOLS::Flavour       * p_flavours;
+    ATOOLS::Flavour       * p_physicalflavours;
     ATOOLS::Vec4D         * p_momenta;
     std::vector<MEPair>     m_mes;    // MEs with their relative factors
     std::vector<CurrentsPair> m_currents;

HADRONS++/Main/Hadron_Decay_Map.C

@@ -72 +72 @@
   }
 }
 
-void Hadron_Decay_Map::Read()
+void Hadron_Decay_Map::ReadHadronAliases(const string& path, const string& file)
 {
+  Data_Reader reader = Data_Reader("->", ";", "#", "");
+  reader.AddWordSeparator("\t");
+  reader.SetAddCommandLine(false);
+  reader.AddComment("#");
+  reader.AddComment("//");
+  reader.SetInputPath(path);
+  reader.SetInputFile(file);
+  
+  vector<vector<string> > aliases;
+  reader.MatrixFromFile(aliases);
+  
+  for (size_t i=0;i<aliases.size();++i) {
+    if (aliases[i].size()!=2) {
+      msg_Error()<<METHOD<<": Wrong syntax in hadron alias file."<<endl
+          <<"  "<<aliases[i]<<endl;
+    }
+    kf_code alias = ToType<kf_code>(aliases[i][0]);
+    kf_code real = ToType<kf_code>(aliases[i][1]);
+    m_startmd.m_aliases[alias]=real;
+    Particle_Info* aliasinfo = new Particle_Info(*s_kftable[real]);
+    aliasinfo->m_kfc=alias;
+    s_kftable[alias]=aliasinfo;
+    msg_Info()<<METHOD<<" created alias "<<alias<<" for "<<Flavour(alias)<<endl;
+  }
+}
+
+void Hadron_Decay_Map::Read(const string& path, const string& file)
+{
   Data_Reader reader = Data_Reader(" ",";","!","->");
   reader.AddWordSeparator("\t");
   reader.SetAddCommandLine(false);
   reader.AddComment("#");
   reader.AddComment("//");
-  reader.SetInputPath(m_decaypath);
-  reader.SetInputFile(m_decayfile);
+  reader.SetInputPath(path);
+  reader.SetInputFile(file);
   
   vector<vector<string> > Decayers;
-  if(!reader.MatrixFromFile(Decayers)) {
-    msg_Error()<<METHOD<<" Warning:\n"
-      <<"   No decayers found in "<<m_decaypath<<m_decayfile<<"."<<endl;
-  }
-  else {
-    msg_Info()<<METHOD<<":"<<endl
-              <<"   Initializing hadron decay tables. This may take some time."
+  if(reader.MatrixFromFile(Decayers)) {
+    msg_Info()<<METHOD<<":"
+              <<"   Initializing from "<<file<<", this may take some time."
               <<endl;
   }