EvtD0toKSpi0etaPlot.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: EvtD0toKspi0etaplot.cc
//
// Modification history:
//
//    Liaoyuan Dong    August, 2022       Module created
//
//------------------------------------------------------------------------
#include "EvtGenBase/EvtPatches.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtGenKine.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenModels/EvtD0toKSpi0etaPlot.hh"
#include "EvtGenBase/EvtRandom.hh"
#include <stdlib.h>
#include <string>
 
#include "TFile.h"  
#include "TApplication.h"
#include "TROOT.h"
#include "TH1.h"
#include "TAxis.h"
#include "TH2.h"
using std::endl;
 
EvtD0toKSpi0etaPlot::~EvtD0toKSpi0etaPlot() {}
 
void EvtD0toKSpi0etaPlot::getName(std::string& model_name){
        model_name="D0toKSpi0etaPlot";     
}
 
EvtDecayBase* EvtD0toKSpi0etaPlot::clone(){
        return new EvtD0toKSpi0etaPlot;
}
 
void EvtD0toKSpi0etaPlot::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("pi0"));
    bool idPi = getDaugs()[2]==EvtPDL::getId(std::string("eta"));
    if(!(idN && idKs && idPi ) ){std::cout<<"EvtD0toKSpi0etaPlot: the daughter sequence should be Ks/Kl pi0 eta"<<std::endl;abort();}
 
    EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getParentId());
 
Xmin = 0.4;
Xmax = 1.734;
Xwid = 0.1334;
Ymin = 0.466;
Ymax = 1.869;
Ywid = 0.1403;
avm1  = 0.125;
double HisPDFtmp[12][12] =
{
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0.00527169, 0.0113872, 0.0301807, 0.0804598, 0.0617421, 0.0427827, 0.0540283, 0.0165289, 0, 0, 0},
{0, 0.020155, 0.0159653, 0.0296171, 0.0808365, 0.0416933, 0.0265935, 0.0193752, 0.0441176, 0.0259259, 0.0287037, 0},
{0, 0.0417141, 0.0513983, 0.0479749, 0.0659032, 0.0323741, 0.0175809, 0.0304969, 0.0254219, 0.054984, 0.0801812, 0.0493827},
{0, 0.0388568, 0.00614035, 0.0154096, 0.028749, 0.0281987, 0.0182064, 0.00844951, 0.00682977, 0.0027885, 0.00239316, -0.0173611},
{0, 0.0100519, 0.0169912, 0.0200723, 0.0272782, 0.0163995, 0.0231309, 0.0170477, 0.013863, 0.000892459, 0.00420168, 0},
{0, 0.0361953, 0.0251366, 0.0173922, 0.0399432, 0.00877736, 0.0105188, 0.0110497, 0.0107359, 0.00134903, 0.00448934, 0},
{0, 0.0078219, 0.0234282, 0.0324155, 0.0451149, 0.011002, 0.0149195, 0.00172265, 0.00821524, 0, 0, 0},
{0, 0.021721, 0.0280722, 0.0341469, 0.0507545, 0.0136752, 0.011148, 0.0100806, 0.00383142, 0, 0, 0},
{0, 0, 0.0237679, 0.0316915, 0.0674883, 0.018785, 0.0192308, 0.0037594, 0.047619, 0, 0, 0},
{0, 0, 0.0307692, 0.0461138, 0.0755956, 0.0250696, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0.0938697, 0.125, 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 EvtD0toKSpi0etaPlot::initProbMax(){
        noProbMax();
}
 
void EvtD0toKSpi0etaPlot::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 EvtD0toKSpi0etaPlot::FindXBin(double mass2) {
  if (mass2 < Xmin) { return 0; }
  else if (mass2>=Xmax) { return 11; }
  else { return int((mass2-Xmin)/Xwid)+1; }
}
 
int EvtD0toKSpi0etaPlot::FindYBin(double mass2) {
  if (mass2 < Ymin) { return 0; }
  else if (mass2>=Ymax) { return 11; }
  else { return int((mass2-Ymin)/Ywid)+1; }
}