Ahadic.C 8.69 KB
Newer Older
1 2 3
#include <cassert>
#include "AHADIC++/Main/Ahadic.H"
#include "AHADIC++/Tools/Soft_Cluster_Handler.H"
4
#include "AHADIC++/Tools/Dipole.H"
5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28
#include "AHADIC++/Tools/Cluster.H"
#include "ATOOLS/Org/Message.H"
#include "ATOOLS/Org/Run_Parameter.H"
#include "AHADIC++/Tools/Hadron_Init.H"
#include "AHADIC++/Tools/Hadronisation_Parameters.H"
#include "AHADIC++/Formation/Cluster_Formation_Handler.H"
#include "AHADIC++/Decays/Cluster_Decay_Handler.H"
#include "ATOOLS/Org/Data_Reader.H"

using namespace AHADIC;
using namespace ATOOLS;
using namespace std;


Ahadic::Ahadic(string path,string file,bool ana)  :
  m_fullinfo(false), m_maxtrials(3), m_clulist()
{
  
  Data_Reader dr(" ",";","!","=");
  dr.AddComment("#");
  dr.AddWordSeparator("\t");
  dr.SetInputPath(path);
  dr.SetInputFile(file);

29 30
  hadpars =  new Hadronisation_Parameters();
  hadpars->Init(path,file);
31 32
  ana=false;

33 34
  p_cformhandler = new Cluster_Formation_Handler(&m_clulist,hadpars->AnaOn());
  p_cdechandler  = new Cluster_Decay_Handler(&m_clulist,hadpars->AnaOn());
35 36 37 38 39 40 41 42
  msg_Tracking()<<"Initialisation of Ahadic complete."<<endl;
}

Ahadic::~Ahadic() 
{
  CleanUp();
  if (p_cdechandler)  { delete p_cdechandler;  p_cdechandler=NULL;  }
  if (p_cformhandler) { delete p_cformhandler; p_cformhandler=NULL; }
43
  delete hadpars;
44 45 46 47 48
}

Return_Value::code Ahadic::Hadronize(Blob_List * blobs)
{
  static std::string mname(METHOD);
49
  Return_Value::IncCall(mname);
50 51
  Blob * blob(NULL);
  bool moveon(false);
52
  double norm2(sqr(rpa->gen.Ecms()));
53 54 55 56 57
  Return_Value::code result;
  for (Blob_List::iterator blit=blobs->begin();blit!=blobs->end();) {
    if ((*blit)->Has(blob_status::needs_hadronization) &&
	(*blit)->Type()==btp::Fragmentation) {
      blob   = (*blit);
58 59 60 61
      // msg_Out()<<"====================================================\n"
      // 	  <<"====================================================\n"
      // 	  <<"====================================================\n"
      // 	  <<(*blob)<<"\n";
62 63 64 65 66 67
      moveon = false;
      Reset();
      for (short int i=0;i<m_maxtrials;i++) {
	try {
	  result = Hadronize(blob,i);
	} catch (Return_Value::code ret) {
68 69
	  msg_Error()<<"ERROR in "<<METHOD<<" : \n"
		     <<"   Hadronization for blob "<<(*blob)<<"\n"
70 71
		     <<"("<<blob<<"; "<<blob->NInP()<<" -> "
		     <<blob->NOutP()<<") "
72 73
		     <<"did not work out,\n"
		     <<"   will trigger retrying the event.\n";
74 75 76
	  CleanUp(blob);
	  return Return_Value::Retry_Event;
	}
77
	
78 79 80 81 82 83
	switch (result) {
	case Return_Value::Success : 
	  ++blit;
	  moveon = true;
	  break;
	case Return_Value::Retry_Event : 
84
	  {
85
	    blobs->ColorConservation();
86
	    msg_Tracking()<<"ERROR in "<<METHOD<<" :\n"
87 88 89 90
			  <<"   Hadronization for blob "
			  <<"("<<blob<<"; "<<blob->NInP()<<" -> "
			  <<blob->NOutP()<<") "
			  <<"did not work out,"<<std::endl
91
			  <<"   will trigger retrying the event.\n";
92 93 94
	    CleanUp(blob);
	    if (rpa->gen.Beam1().IsLepton() ||
	        rpa->gen.Beam2().IsLepton()) {
95 96
	      msg_Tracking()<<METHOD<<": "
			    <<"Non-hh collision. Request new event instead.\n";
97 98 99
	      return Return_Value::New_Event;
	    }
	    return result;
100
	  }
101 102 103
	case Return_Value::Retry_Method :
	  msg_Tracking()<<"Warning in "<<METHOD<<" : "<<std::endl
			<<"   Hadronization for blob "
104 105 106 107
			<<"("<<blob<<"; "<<blob->NInP()<<" -> "
			<<blob->NOutP()<<") "
			<<"did not work out properly in the "
			<<(i+1)<<"th attempt,"<<std::endl
108
			<<"   retry it "<<m_maxtrials<<" times."<<std::endl;
109
	  Return_Value::IncRetryMethod(mname);
110 111 112 113 114
	  CleanUp(blob);
	  break;
	case Return_Value::Nothing :
	default:
	  msg_Tracking()<<"Warning in "<<METHOD<<":"<<std::endl
115 116
			<<"   Calling Hadronization for Blob("<<blob<<") "
			<<"yields "<<int(result)<<"."<<std::endl
117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137
			<<"   Continue and hope for the best."<<std::endl;
	  moveon = true;
	  ++blit;
	  break;
	}
	if (moveon) break;
      }

      CleanUp();
      moveon = moveon && SanityCheck(blob,norm2);	    
      if (moveon) {
	blob->SetStatus(blob_status::needs_hadrondecays);
	blob->SetType(btp::Fragmentation);
      }
      else {
	CleanUp(blob);
	return Return_Value::Retry_Event;
      }
    }
    else blit++;
  }
Frank Krauss's avatar
Frank Krauss committed
138 139 140 141 142 143 144 145 146 147 148 149
  // for (size_t i(0);i<blob->NOutP();i++) {
  //   if (blob->OutParticle(i)->Momentum().Nan()) {
  //     msg_Out()<<"=======================================================\n"
  // 	       <<"=======================================================\n"
  // 	       <<"=======================================================\n"
  // 	       <<"=======================================================\n"
  // 	       <<"=======================================================\n"
  // 	       <<(*blob)<<"\n";
  //     msg_Out()<<" nan momentum: "<<blob->OutParticle(i)->Momentum()<<"\n";
  //     exit(0);
  //   }
  // }
150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172
  return Return_Value::Success;
}  



Return_Value::code Ahadic::Hadronize(Blob * blob,int retry) {
  assert(m_clulist.empty());
  blob->SetType(btp::Cluster_Formation);
  blob->SetTypeSpec("AHADIC-1.0");


  switch (p_cformhandler->FormClusters(blob)) {
  case -1 : 
    msg_Tracking()<<"ERROR in "<<METHOD<<" :"<<std::endl
		  <<"   Will retry event."<<std::endl;
    p_cformhandler->Reset();
    return Return_Value::Retry_Event;
  case  0 :
    msg_Tracking()<<"ERROR in "<<METHOD<<" :"<<std::endl
		  <<"   Will retry method."<<std::endl;
    p_cformhandler->Reset();
    return Return_Value::Retry_Method;
  case 1 :
173 174
    if (retry>=0) 
      msg_Tracking()<<"   Passed cluster form now ("<<retry<<"th trial).\n";
175 176 177 178 179 180 181 182 183 184 185 186 187
    break;
  }
  
  switch (p_cdechandler->DecayClusters(blob)) {
  case -1 : 
    msg_Tracking()<<"ERROR in "<<METHOD<<" :"<<std::endl
		  <<"   Will retry event."<<std::endl;
    return Return_Value::Retry_Event;
  case  0 :
    msg_Tracking()<<"ERROR in "<<METHOD<<" :"<<std::endl
		  <<"   Will retry method."<<std::endl;
    return Return_Value::Retry_Method;
  case  1 :
188 189
    if (retry>=0) 
      msg_Tracking()<<"   Passed cluster decays now ("<<retry<<"th trial).\n";
190 191 192 193 194 195 196 197 198 199 200 201
    break;
  }


  
  if (!m_clulist.empty()) {
    msg_Error()<<"Error in "<<METHOD<<":"<<std::endl
	       <<"   "<<m_clulist.size()<<" clusters undeleted:"<<std::endl
	       <<m_clulist<<std::endl;
  }
  assert(m_clulist.empty());

202 203 204
  if (blob->CheckMomentumConservation().Abs2()>1.e-6) {
    msg_Error()<<"Error in "<<METHOD<<": blob seem to be fishy: "
	       <<blob->CheckMomentumConservation()<<"\n"
205
	       <<(*blob)<<"\n";
206 207
    exit(1);
  }
208 209 210 211 212 213 214
  return Return_Value::Success;
}


void Ahadic::Reset() {
  if (Cluster::RemainingClusters()) {
    msg_Error()<<"Error in "<<METHOD<<":"<<std::endl
215
	       <<"   "<<Cluster::RemainingClusters()<<" clusters undeleted.\n";
216 217 218 219 220 221 222
  }
  assert(!Cluster::RemainingClusters());
  Cluster::ResetClusterNumber();
  control::s_AHAparticles=0;
}

bool Ahadic::SanityCheck(Blob * blob,double norm2) {
223 224
  Vec4D checkmom(blob->CheckMomentumConservation());
  if (dabs(checkmom.Abs2())/norm2>1.e-12 ||
225 226
      (norm2<0. && norm2>0.) ||
      Cluster::RemainingClusters()!=0 ||
227 228
      control::s_AHAparticles!=blob->NOutP()
      /* || control::s_AHAprotoparticles!=0*/) {
229
    msg_Error()<<"ERROR in "<<METHOD<<" : "<<endl
230 231
	       <<"   Momentum/particle-blob number violation for "
	       <<(Cluster::RemainingClusters())
232 233 234
	       <<" remaining clusters (norm2 = "<<norm2<<")."<<endl
	       <<"   Protoparticles = "<<control::s_AHAprotoparticles
	       <<"/ parts = "<<control::s_AHAparticles<<" vs. "<<blob->NOutP()
235
	       <<"   : "<<checkmom<<" ("<<sqrt(Max(0.,checkmom.Abs2()))<<")\n"
236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279
	       <<(*blob)<<endl;
    return false;
  }
  msg_Tracking()<<"Passed "<<METHOD<<" with "
		<<"protoparticles = "<<control::s_AHAprotoparticles
		<<"/ parts = "<<control::s_AHAparticles<<" vs. "<<blob->NOutP()
		<<"   : "<<blob->CheckMomentumConservation()<<endl;
  msg_Debugging()<<(*blob)<<endl;
  return true;
}

void Ahadic::CleanUp(Blob * blob) {
  if (Cluster::RemainingActives()>0) {
    msg_Tracking()<<METHOD<<": "<<Cluster::RemainingActives()
		  <<" remaining Clusters found:"<<std::endl;
    //Cluster::PrintActives(msg_Tracking());
    Cluster::DeleteActives();
  }
  if (!m_clulist.empty()) m_clulist.clear();
  
  if (Dipole::RemainingActives()>0) {
    msg_Tracking()<<METHOD<<": "<<Dipole::RemainingActives()
		  <<" remaining Dipoles found:"<<std::endl;
    //Dipole::PrintActives(msg_Tracking());
    Dipole::DeleteActives();
  }

  if (Proto_Particle::RemainingActives()>0) {
    msg_Tracking()<<METHOD<<": "<<Proto_Particle::RemainingActives()
		  <<" remaining Proto_Particles found:"<<std::endl;
    //Proto_Particle::PrintActives(msg_Tracking());
    Proto_Particle::DeleteActives();
  }

  if (Proto_Particle_List::RemainingActives()>0) {
    msg_Tracking()<<METHOD<<": "<<Proto_Particle_List::RemainingActives()
		  <<" remaining Proto_Particle_Lists found:"<<std::endl;
    //Proto_Particle_List::PrintActives(msg_Tracking());
    Proto_Particle_List::DeleteActives();
  }


  if(blob) blob->DeleteOutParticles(0);
}