EVGEN/AliGenMC.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /*
  17 $Log$
  18 Revision 1.6  2002/03/26 14:19:36  morsch
  19 Saver calculation of rapdity.
  20
  21 Revision 1.5  2002/03/12 17:02:20  morsch
  22 Change in calculation of rapidity, include case in which numerically e == pz.
  23
  24 Revision 1.4  2001/11/27 13:13:07  morsch
  25 Maximum lifetime for long-lived particles to be put on the stack is parameter.
  26 It can be set via SetMaximumLifetime(..).
  27
  28 Revision 1.3  2001/10/16 08:48:56  morsch
  29 Common vertex related code moved to base class AliGenerator.
  30
  31 Revision 1.2  2001/10/15 08:15:51  morsch
  32 Event vertex and vertex truncation setting moved into AliMC.
  33
  34 Revision 1.1  2001/07/13 10:56:00  morsch
  35 AliGenMC base class for AliGenParam and AliGenPythia commonalities.
  36
  37 */
  38
  39 // Base class for generators using external MC generators.
  40 // For example AliGenPythia using Pythia.
  41 // Provides basic functionality: setting of kinematic cuts on
  42 // decay products and particle selection.
  43 // andreas.morsch@cern.ch
  44
  45 #include "AliGenMC.h"
  46 #include "AliPDG.h"
  47 #include <TParticle.h>
  48
  49  ClassImp(AliGenMC)
  50
  51 AliGenMC::AliGenMC()
  52                  :AliGenerator()
  53 {
  54 // Default Constructor
  55     SetCutOnChild();
  56     SetChildMomentumRange();
  57     SetChildPtRange();
  58     SetChildPhiRange();
  59     SetChildThetaRange();
  60     SetChildYRange();
  61     SetMaximumLifetime();
  62 }
  63
  64 AliGenMC::AliGenMC(Int_t npart)
  65                  :AliGenerator(npart)
  66 {
  67 //  Constructor
  68     SetCutOnChild();
  69     SetChildMomentumRange();
  70     SetChildPtRange();
  71     SetChildPhiRange();
  72     SetChildThetaRange();
  73     SetChildYRange();
  74 //
  75     fParentSelect.Set(8);
  76     fChildSelect.Set(8);
  77     for (Int_t i=0; i<8; i++) fParentSelect[i]=fChildSelect[i]=0;
  78     SetMaximumLifetime();
  79 }
  80
  81 AliGenMC::AliGenMC(const AliGenMC & mc)
  82 {
  83 // copy constructor
  84 }
  85
  86 AliGenMC::~AliGenMC()
  87 {
  88 // Destructor
  89 }
  90
  91 void AliGenMC::Init()
  92 {
  93 //
  94 //  Initialization
  95     switch (fForceDecay) {
  96     case kSemiElectronic:
  97     case kDiElectron:
  98     case kBJpsiDiElectron:
  99     case kBPsiPrimeDiElectron:
 100         fChildSelect[0] = kElectron;
 101         break;
 102     case kSemiMuonic:
 103     case kDiMuon:
 104     case kBJpsiDiMuon:
 105     case kBPsiPrimeDiMuon:
 106     case kPiToMu:
 107     case kKaToMu:
 108         fChildSelect[0]=kMuonMinus;
 109         break;
 110     case kHadronicD:
 111         fChildSelect[0]=kPiPlus;
 112         fChildSelect[1]=kKPlus;
 113         break;
 114     case kOmega:
 115     case kAll:
 116     case kNoDecay:
 117         break;
 118     }
 119 }
 120
 121
 122 Bool_t AliGenMC::ParentSelected(Int_t ip) const
 123 {
 124 // True if particle is in list of parent particles to be selected
 125     for (Int_t i=0; i<8; i++)
 126     {
 127         if (fParentSelect.At(i) == ip) return kTRUE;
 128     }
 129     return kFALSE;
 130 }
 131
 132 Bool_t AliGenMC::ChildSelected(Int_t ip) const
 133 {
 134 // True if particle is in list of decay products to be selected
 135     for (Int_t i=0; i<5; i++)
 136     {
 137         if (fChildSelect.At(i) == ip) return kTRUE;
 138     }
 139     return kFALSE;
 140 }
 141
 142 Bool_t AliGenMC::KinematicSelection(TParticle *particle, Int_t flag) const
 143 {
 144 // Perform kinematic selection
 145     Float_t px    = particle->Px();
 146     Float_t py    = particle->Py();
 147     Float_t pz    = particle->Pz();
 148     Float_t  e    = particle->Energy();
 149     Float_t pt    = particle->Pt();
 150     Float_t p     = particle->P();
 151     Float_t theta = particle->Theta();
 152     Float_t mass  = particle->GetMass();
 153     Float_t mt2   = pt * pt + mass * mass;
 154
 155     Float_t phi   = Float_t(TMath::ATan2(Double_t(py),Double_t(px)));
 156     Double_t y, y0;
 157
 158     if (TMath::Abs(pz) <  e) {
 159         y = 0.5*TMath::Log((e+pz)/(e-pz));
 160     } else {
 161         y = 1.e10;
 162     }
 163
 164     if (mt2) {
 165         y0 = 0.5*TMath::Log((e+TMath::Abs(pz))*(e+TMath::Abs(pz))/mt2);
 166     } else {
 167         if (TMath::Abs(y) < 1.e10) {
 168             y0 = y;
 169         } else {
 170             y0 = 1.e10;
 171         }
 172     }
 173
 174     y = (pz < 0) ? -y0 : y0;
 175
 176     if (flag == 0) {
 177 //
 178 // Primary particle cuts
 179 //
 180 //  transverse momentum cut
 181         if (pt > fPtMax || pt < fPtMin) {
 182 //          printf("\n failed pt cut %f %f %f \n",pt,fPtMin,fPtMax);
 183             return kFALSE;
 184         }
 185 //
 186 // momentum cut
 187         if (p > fPMax || p < fPMin) {
 188 //          printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
 189             return kFALSE;
 190         }
 191 //
 192 // theta cut
 193         if (theta > fThetaMax || theta < fThetaMin) {
 194 //          printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
 195             return kFALSE;
 196         }
 197 //
 198 // rapidity cut
 199         if (y > fYMax || y < fYMin) {
 200 //          printf("\n failed y cut %f %f %f \n",y,fYMin,fYMax);
 201             return kFALSE;
 202         }
 203 //
 204 // phi cut
 205         if (phi > fPhiMax || phi < fPhiMin) {
 206 //          printf("\n failed phi cut %f %f %f \n",phi,fPhiMin,fPhiMax);
 207             return kFALSE;
 208         }
 209     } else {
 210 //
 211 // Decay product cuts
 212 //
 213 //  transverse momentum cut
 214         if (pt > fChildPtMax || pt < fChildPtMin) {
 215 //          printf("\n failed pt cut %f %f %f \n",pt,fChildPtMin,fChildPtMax);
 216             return kFALSE;
 217         }
 218 //
 219 // momentum cut
 220         if (p > fChildPMax || p < fChildPMin) {
 221 //          printf("\n failed p cut %f %f %f \n",p,fChildPMin,fChildPMax);
 222             return kFALSE;
 223         }
 224 //
 225 // theta cut
 226         if (theta > fChildThetaMax || theta < fChildThetaMin) {
 227 //          printf("\n failed theta cut %f %f %f \n",theta,fChildThetaMin,fChildThetaMax);
 228             return kFALSE;
 229         }
 230 //
 231 // rapidity cut
 232         if (y > fChildYMax || y < fChildYMin) {
 233 //          printf("\n failed y cut %f %f %f \n",y,fChildYMin,fChildYMax);
 234             return kFALSE;
 235         }
 236 //
 237 // phi cut
 238         if (phi > fChildPhiMax || phi < fChildPhiMin) {
 239 //          printf("\n failed phi cut %f %f %f \n",phi,fChildPhiMin,fChildPhiMax);
 240             return kFALSE;
 241         }
 242     }
 243
 244
 245
 246     return kTRUE;
 247 }
 248
 249 Int_t AliGenMC::CheckPDGCode(Int_t pdgcode) const
 250 {
 251 //
 252 //  If the particle is in a diffractive state, then take action accordingly
 253   switch (pdgcode) {
 254   case 91:
 255     return 92;
 256   case 110:
 257     //rho_diff0 -- difficult to translate, return rho0
 258     return 113;
 259   case 210:
 260     //pi_diffr+ -- change to pi+
 261     return 211;
 262   case 220:
 263     //omega_di0 -- change to omega0
 264     return 223;
 265   case 330:
 266     //phi_diff0 -- return phi0
 267     return 333;
 268   case 440:
 269     //J/psi_di0 -- return J/psi
 270     return 443;
 271   case 2110:
 272     //n_diffr -- return neutron
 273     return 2112;
 274   case 2210:
 275     //p_diffr+ -- return proton
 276     return 2212;
 277   }
 278   //non diffractive state -- return code unchanged
 279   return pdgcode;
 280 }
 281
 282 AliGenMC& AliGenMC::operator=(const  AliGenMC& rhs)
 283 {
 284 // Assignment operator
 285     return *this;
 286 }
 287