TEvtGen/EvtGenModels/EvtHelAmp.cpp

   1 //--------------------------------------------------------------------------
   2 //
   3 // Environment:
   4 //      This software is part of the EvtGen package developed jointly
   5 //      for the BaBar and CLEO collaborations.  If you use all or part
   6 //      of it, please give an appropriate acknowledgement.
   7 //
   8 // Copyright Information: See EvtGen/COPYRIGHT
   9 //      Copyright (C) 1998      Caltech, UCSB
  10 //
  11 // Module: EvtHelAmp.cc
  12 //
  13 // Description: Decay model for implementation of generic 2 body
  14 //              decay specified by the helicity amplitudes
  15 //
  16 //
  17 // Modification history:
  18 //
  19 //    RYD       March 14, 1999       Module created
  20 //
  21 //------------------------------------------------------------------------
  22 //
  23 #include "EvtGenBase/EvtPatches.hh"
  24 #include <stdlib.h>
  25 #include "EvtGenBase/EvtParticle.hh"
  26 #include "EvtGenBase/EvtGenKine.hh"
  27 #include "EvtGenBase/EvtPDL.hh"
  28 #include "EvtGenBase/EvtReport.hh"
  29 #include "EvtGenModels/EvtHelAmp.hh"
  30 #include "EvtGenBase/EvtId.hh"
  31 #include <string>
  32 #include "EvtGenBase/EvtConst.hh"
  33 #include "EvtGenBase/EvtEvalHelAmp.hh"
  34 using std::endl;
  35
  36
  37 EvtHelAmp::~EvtHelAmp() {
  38
  39   delete _evalHelAmp;
  40
  41 }
  42
  43 std::string EvtHelAmp::getName(){
  44
  45   return "HELAMP";
  46
  47 }
  48
  49
  50 EvtDecayBase* EvtHelAmp::clone(){
  51
  52   return new EvtHelAmp;
  53
  54 }
  55
  56 void EvtHelAmp::init(){
  57
  58   checkNDaug(2);
  59
  60
  61   //find out how many states each particle have
  62   int _nA=EvtSpinType::getSpinStates(EvtPDL::getSpinType(getParentId()));
  63   int _nB=EvtSpinType::getSpinStates(EvtPDL::getSpinType(getDaug(0)));
  64   int _nC=EvtSpinType::getSpinStates(EvtPDL::getSpinType(getDaug(1)));
  65
  66   if (verbose()){
  67     report(INFO,"EvtGen")<<"_nA,_nB,_nC:"
  68                          <<_nA<<","<<_nB<<","<<_nC<<endl;
  69   }
  70
  71   //find out what 2 times the spin is
  72   int _JA2=EvtSpinType::getSpin2(EvtPDL::getSpinType(getParentId()));
  73   int _JB2=EvtSpinType::getSpin2(EvtPDL::getSpinType(getDaug(0)));
  74   int _JC2=EvtSpinType::getSpin2(EvtPDL::getSpinType(getDaug(1)));
  75
  76   if (verbose()){
  77     report(INFO,"EvtGen")<<"_JA2,_JB2,_JC2:"
  78                          <<_JA2<<","<<_JB2<<","<<_JC2<<endl;
  79   }
  80
  81   //allocate memory
  82   int* _lambdaA2=new int[_nA];
  83   int* _lambdaB2=new int[_nB];
  84   int* _lambdaC2=new int[_nC];
  85
  86   EvtComplexPtr* _HBC=new EvtComplexPtr[_nB];
  87   for(int ib=0;ib<_nB;ib++){
  88     _HBC[ib]=new EvtComplex[_nC];
  89   }
  90
  91   int i;
  92   //find the allowed helicities (actually 2*times the helicity!)
  93
  94   fillHelicity(_lambdaA2,_nA,_JA2,getParentId());
  95   fillHelicity(_lambdaB2,_nB,_JB2,getDaug(0));
  96   fillHelicity(_lambdaC2,_nC,_JC2,getDaug(1));
  97
  98   if (verbose()){
  99     report(INFO,"EvtGen")<<"Helicity states of particle A:"<<endl;
 100     for(i=0;i<_nA;i++){
 101       report(INFO,"EvtGen")<<_lambdaA2[i]<<endl;
 102     }
 103
 104     report(INFO,"EvtGen")<<"Helicity states of particle B:"<<endl;
 105     for(i=0;i<_nB;i++){
 106       report(INFO,"EvtGen")<<_lambdaB2[i]<<endl;
 107     }
 108
 109     report(INFO,"EvtGen")<<"Helicity states of particle C:"<<endl;
 110     for(i=0;i<_nC;i++){
 111       report(INFO,"EvtGen")<<_lambdaC2[i]<<endl;
 112     }
 113   }
 114
 115   //now read in the helicity amplitudes
 116
 117   int argcounter=0;
 118
 119   for(int ib=0;ib<_nB;ib++){
 120     for(int ic=0;ic<_nC;ic++){
 121       _HBC[ib][ic]=0.0;
 122       if (abs(_lambdaB2[ib]-_lambdaC2[ic])<=_JA2) argcounter+=2;
 123     }
 124   }
 125
 126   checkNArg(argcounter);
 127
 128   argcounter=0;
 129
 130   for(int ib=0;ib<_nB;ib++){
 131     for(int ic=0;ic<_nC;ic++){
 132       if (abs(_lambdaB2[ib]-_lambdaC2[ic])<=_JA2) {
 133         _HBC[ib][ic]=getArg(argcounter)*exp(EvtComplex(0.0,getArg(argcounter+1)));;
 134         argcounter+=2;
 135         if (verbose()){
 136           report(INFO,"EvtGen")<<"_HBC["<<ib<<"]["<<ic<<"]="
 137                                <<_HBC[ib][ic]<<endl;
 138         }
 139       }
 140     }
 141   }
 142
 143   _evalHelAmp=new EvtEvalHelAmp(getParentId(),
 144                                 getDaug(0),
 145                                 getDaug(1),
 146                                 _HBC);
 147
 148   // Note: these are not class data members but local variables.
 149   delete [] _lambdaA2;
 150   delete [] _lambdaB2;
 151   delete [] _lambdaC2;
 152   for(int ib=0;ib<_nB;ib++){
 153     delete [] _HBC[ib];
 154   }
 155   delete [] _HBC;  // _HBC is copied in ctor of EvtEvalHelAmp above.
 156
 157 }
 158
 159
 160 void EvtHelAmp::initProbMax(){
 161
 162   double maxprob=_evalHelAmp->probMax();
 163
 164   if (verbose()){
 165     report(INFO,"EvtGen")<<"Calculated probmax"<<maxprob<<endl;
 166   }
 167
 168   setProbMax(maxprob);
 169
 170 }
 171
 172
 173 void EvtHelAmp::decay( EvtParticle *p){
 174
 175   //first generate simple phase space
 176   p->initializePhaseSpace(getNDaug(),getDaugs());
 177
 178   _evalHelAmp->evalAmp(p,_amp2);
 179
 180   return ;
 181
 182 }
 183
 184
 185 void EvtHelAmp::fillHelicity(int* lambda2,int n,int J2, EvtId id){
 186
 187   int i;
 188
 189   //photon is special case!
 190   if (n==2&&J2==2) {
 191     lambda2[0]=2;
 192     lambda2[1]=-2;
 193     return;
 194   }
 195
 196   //and so is the neutrino!
 197   if (n==1&&J2==1) {
 198     if (EvtPDL::getStdHep(id)>0){
 199         //particle i.e. lefthanded
 200         lambda2[0]=-1;
 201     }else{
 202         //anti particle i.e. righthanded
 203         lambda2[0]=1;
 204     }
 205     return;
 206   }
 207
 208   assert(n==J2+1);
 209
 210   for(i=0;i<n;i++){
 211     lambda2[i]=n-i*2-1;
 212   }
 213
 214   return;
 215
 216 }
 217
 218
 219
 220
 221
 222
 223