EvtDtoKSpietaPlot.cc

Go to the documentation of this file.

//--------------------------------------------------------------------------
//
// Environment:
//      This software is part of models developed at BES collaboration
//      based on the EvtGen framework.  If you use all or part
//      of it, please give an appropriate acknowledgement.
//
// Copyright Information: See EvtGen/BesCopyright
//      Copyright (A) 2006      Ping Rong-Gang @IHEP
//
// Module: EvtDtoKSpietaPlot.cc
//
// Modification history:
//
//    Liaoyuan Dong    August, 2022       Module created
//
//------------------------------------------------------------------------
#include "EvtDtoKSpietaPlot.hh"
#include "../EvtGenBase/EvtGenKine.hh"
#include "../EvtGenBase/EvtPDL.hh"
#include "../EvtGenBase/EvtParticle.hh"
#include "../EvtGenBase/EvtPatches.hh"
#include "../EvtGenBase/EvtRandom.hh"
#include "../EvtGenBase/EvtReport.hh"
#include <stdlib.h>
#include <string>
 
#include "TApplication.h"
#include "TAxis.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TROOT.h"
using std::endl;
 
EvtDtoKSpietaPlot::~EvtDtoKSpietaPlot() {}
 
void EvtDtoKSpietaPlot::getName( std::string& model_name ) { model_name = "DtoKSpietaPlot"; }
 
EvtDecayBase* EvtDtoKSpietaPlot::clone() { return new EvtDtoKSpietaPlot; }
 
void EvtDtoKSpietaPlot::init() {
 
  checkNArg( 0 );
 
  bool idN = getDaugs()[0] == EvtPDL::getId( std::string( "K_S0" ) ) ||
             getDaugs()[0] == EvtPDL::getId( std::string( "K_L0" ) );
  bool idKs = getDaugs()[1] == EvtPDL::getId( std::string( "pi+" ) ) ||
              getDaugs()[1] == EvtPDL::getId( std::string( "pi-" ) );
  bool idPi = getDaugs()[2] == EvtPDL::getId( std::string( "eta" ) );
  if ( !( idN && idKs && idPi ) )
  {
    std::cout << "EvtDtoKSpietaPlot: the daughter sequence should be Ks/Kl pi+ eta"
              << std::endl;
    abort();
  }
  EvtSpinType::spintype parenttype = EvtPDL::getSpinType( getParentId() );
 
  Xmin                     = 0.406;
  Xmax                     = 1.747;
  Xwid                     = 0.1341;
  Ymin                     = 0.473;
  Ymax                     = 1.882;
  Ywid                     = 0.1409;
  avm1                     = 0.240741;
  double HisPDFtmp[12][12] = { { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
                               { 0, 1e-07, 0.0119599, 0.0296515, 0.0444187, 0.0196774,
                                 0.00679957, 0.00950145, 0.00704225, 0.0833333, 1e-07, 0 },
                               { 0, 0.00265756, 0.00594999, 0.04354, 0.0266106, 0.00888889,
                                 0.0104167, 0.00903452, 0.0110731, 0.00418803, 1e-07, 0 },
                               { 0, 0.00629828, 0.0145145, 0.0280041, 0.0236976, 0.0160962,
                                 0.00604595, 0.0225225, 0.00838338, 0.00828402, 0.0102041, 0 },
                               { 0, 0.0058566, 0.0214311, 0.0238095, 0.0421477, 0.0116178,
                                 0.0188645, 0.0262201, 0.00917431, 0.00817757, 0.00310856, 0 },
                               { 0, 0.00679852, 0.0101833, 0.0106561, 0.0188282, 0.0231757,
                                 0.0246205, 0.0146294, 0.0124968, 0.0152171, 0.00802826, 0 },
                               { 0, 0.01116, 0.00961875, 0.0153356, 0.0245832, 0.0207115,
                                 0.0180312, 0.01606, 0.0126437, 0.00146566, 0.00326797, 0 },
                               { 0, 0.005919, 0.00744745, 0.0114328, 0.0330096, 0.0306159,
                                 0.0280269, 0.0182468, 0.0111238, 0.00387888, 1e-07, 0 },
                               { 0, 0.00185185, 0.00144012, 0.0189673, 0.0478907, 0.0233768,
                                 0.0250272, 0.0130943, 0.0166364, 1e-07, 1e-07, 0 },
                               { 0, 1e-07, 0.00462388, 0.0170154, 0.0512206, 0.0407905,
                                 0.0289623, 0.0307394, 1e-07, 1e-07, 1e-07, 0 },
                               { 0, 1e-07, 0.00657895, 0.0152134, 0.0479368, 0.0486111,
                                 0.0318021, 0.166667, 1e-07, 1e-07, 1e-07, 0 },
                               { 0, 0, 0, 0, 0.0689655, 0.240741, 0, 0, 0, 0, 0, 0 } };
 
  for ( int i = 0; i < 12; i++ )
  {
    for ( int j = 0; j < 12; j++ ) { HisPDF[i][j] = HisPDFtmp[i][j]; }
  }
}
void EvtDtoKSpietaPlot::initProbMax() { noProbMax(); }
 
void EvtDtoKSpietaPlot::decay( EvtParticle* p ) {
 
loop:
  p->initializePhaseSpace( getNDaug(), getDaugs() );
 
  EvtParticle *id1, *id2, *id3;
  EvtVector4R  pd1, pd2, pd3;
  double       xmass13, xmass12, xmass23;
 
  id1 = p->getDaug( 0 );
  id2 = p->getDaug( 1 );
  id3 = p->getDaug( 2 );
 
  pd1 = id1->getP4Lab();
  pd2 = id2->getP4Lab();
  pd3 = id3->getP4Lab();
 
  xmass12 = ( pd1 + pd2 ).mass() * ( pd1 + pd2 ).mass(); // M_ksopi
  //        xmass13=(pd1+pd3).mass()*(pd1+pd3).mass();  // M_ksow
  xmass23 = ( pd2 + pd3 ).mass() * ( pd2 + pd3 ).mass(); // M_piw
 
  int    xbin     = FindXBin( xmass12 );
  int    ybin     = FindYBin( xmass23 );
  double xratio12 = HisPDF[xbin][ybin] / avm1;
 
  if ( xratio12 <= 0 ) goto loop;
 
  double rd12 = EvtRandom::Flat( 0.0, 1.0 );
  if ( rd12 > xratio12 ) goto loop;
 
  return;
}
 
int EvtDtoKSpietaPlot::FindXBin( double mass2 ) {
  if ( mass2 < Xmin ) { return 0; }
  else if ( mass2 >= Xmax ) { return 11; }
  else { return int( ( mass2 - Xmin ) / Xwid ) + 1; }
}
 
int EvtDtoKSpietaPlot::FindYBin( double mass2 ) {
  if ( mass2 < Ymin ) { return 0; }
  else if ( mass2 >= Ymax ) { return 11; }
  else { return int( ( mass2 - Ymin ) / Ywid ) + 1; }
}