EVGEN/AliGenReaderEMD.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 // Class to read events from external (TNtupla) file
  18 // Events -> neutron removal by EM dissociation of Pb nuclei
  19 // Data from RELDIS code (by I. Pshenichov)
  20
  21 #include <TFile.h>
  22 #include <TParticle.h>
  23 #include <TTree.h>
  24 #include <TVirtualMC.h>
  25 #include <TDatabasePDG.h>
  26 #include <TPDGCode.h>
  27 #include "AliGenReaderEMD.h"
  28 #include "AliStack.h"
  29
  30
  31 ClassImp(AliGenReaderEMD)
  32
  33 AliGenReaderEMD::AliGenReaderEMD():
  34     fStartEvent(0),
  35     fNcurrent(0),
  36     fNparticle(0),
  37     fTreeNtuple(0),
  38     fPcToTrack(0),
  39     fOffset(0),
  40     fNnAside(0),
  41     fEnAside(0),
  42     fnPDGCode(0),
  43     fNnCside(0),
  44     fEnCside(0),
  45     fNpAside(0),
  46     fEtapAside(0),
  47     fpPDGCode(0),
  48     fNpCside(0),
  49     fEtapCside(0),
  50     fNppAside(0),
  51     fEtappAside(0),
  52     fppPDGCode(0),
  53     fNppCside(0),
  54     fEtappCside(0),
  55     fNpmAside(0),
  56     fEtapmAside(0),
  57     fpmPDGCode(0),
  58     fNpmCside(0),
  59     fEtapmCside(0),
  60     fNp0Aside(0),
  61     fEtap0Aside(0),
  62     fp0PDGCode(0),
  63     fNp0Cside(0),
  64     fEtap0Cside(0),
  65     fNetaAside(0),
  66     fEtaetaAside(0),
  67     fetaPDGCode(0),
  68     fNetaCside(0),
  69     fEtaetaCside(0),
  70     fNomegaAside(0),
  71     fEtaomegaAside(0),
  72     fomegaPDGCode(0),
  73     fNomegaCside(0),
  74     fEtaomegaCside(0)
  75 {
  76 // Std constructor
  77     for(int i=0; i<70; i++){
  78        fPxnAside[i] = fPynAside[i] = fPznAside[i] = 0.;
  79        fPxnCside[i] = fPynCside[i] = fPznCside[i] = 0.;
  80        if(i<50){
  81          fPxpAside[i] = fPypAside[i] = fPzpAside[i] = 0.;
  82          fPxpCside[i] = fPypCside[i] = fPzpCside[i] = 0.;
  83          if(i<30){
  84            fPxppAside[i] = fPyppAside[i] = fPzppAside[i] = 0.;
  85            fPxppCside[i] = fPyppCside[i] = fPzppCside[i] = 0.;
  86            fPxpmAside[i] = fPypmAside[i] = fPzpmAside[i] = 0.;
  87            fPxpmCside[i] = fPypmCside[i] = fPzpmCside[i] = 0.;
  88            fPxp0Aside[i] = fPyp0Aside[i] = fPzp0Aside[i] = 0.;
  89            fPxp0Cside[i] = fPyp0Cside[i] = fPzp0Cside[i] = 0.;
  90            if(i<15){
  91              fPxetaAside[i] = fPyetaAside[i] = fPzetaAside[i] = 0.;
  92              fPxetaCside[i] = fPyetaCside[i] = fPzetaCside[i] = 0.;
  93              fPxomegaAside[i] = fPyomegaAside[i] = fPzomegaAside[i] = 0.;
  94              fPxomegaCside[i] = fPyomegaCside[i] = fPzomegaCside[i] = 0.;
  95            }
  96          }
  97        }
  98     }
  99     if(fPcToTrack==kAll) printf("\n\t   *** AliGenReaderEMD will track all produced particles \n\n");
 100     else if(fPcToTrack==kNotNucleons) printf("\n\t   *** AliGenReaderEMD will track all produced particles except nucleons\n\n");
 101     else if(fPcToTrack==kOnlyNucleons) printf("\n\t   *** AliGenReaderEMD will track only nucleons\n\n");
 102 }
 103
 104
 105 AliGenReaderEMD::AliGenReaderEMD(const AliGenReaderEMD &reader):
 106     AliGenReader(reader),
 107     fStartEvent(0),
 108     fNcurrent(0),
 109     fNparticle(0),
 110     fTreeNtuple(0),
 111     fPcToTrack(0),
 112     fOffset(0),
 113     fNnAside(0),
 114     fEnAside(0),
 115     fnPDGCode(0),
 116     fNnCside(0),
 117     fEnCside(0),
 118     fNpAside(0),
 119     fEtapAside(0),
 120     fpPDGCode(0),
 121     fNpCside(0),
 122     fEtapCside(0),
 123     fNppAside(0),
 124     fEtappAside(0),
 125     fppPDGCode(0),
 126     fNppCside(0),
 127     fEtappCside(0),
 128     fNpmAside(0),
 129     fEtapmAside(0),
 130     fpmPDGCode(0),
 131     fNpmCside(0),
 132     fEtapmCside(0),
 133     fNp0Aside(0),
 134     fEtap0Aside(0),
 135     fp0PDGCode(0),
 136     fNp0Cside(0),
 137     fEtap0Cside(0),
 138     fNetaAside(0),
 139     fEtaetaAside(0),
 140     fetaPDGCode(0),
 141     fNetaCside(0),
 142     fEtaetaCside(0),
 143     fNomegaAside(0),
 144     fEtaomegaAside(0),
 145     fomegaPDGCode(0),
 146     fNomegaCside(0),
 147     fEtaomegaCside(0)
 148 {
 149     // Copy Constructor
 150     reader.Copy(*this);
 151 }
 152   // -----------------------------------------------------------------------------------
 153 AliGenReaderEMD::~AliGenReaderEMD()
 154 {
 155     delete fTreeNtuple;
 156 }
 157
 158 // -----------------------------------------------------------------------------------
 159 AliGenReaderEMD& AliGenReaderEMD::operator=(const  AliGenReaderEMD& rhs)
 160 {
 161 // Assignment operator
 162     rhs.Copy(*this);
 163     return *this;
 164 }
 165
 166 // -----------------------------------------------------------------------------------
 167 void AliGenReaderEMD::Copy(TObject&) const
 168 {
 169     //
 170     // Copy
 171     //
 172     Fatal("Copy","Not implemented!\n");
 173 }
 174
 175 // -----------------------------------------------------------------------------------
 176 void AliGenReaderEMD::Init()
 177 {
 178 //
 179 // Reset the existing file environment and open a new root file
 180
 181     TFile *pFile=0;
 182     if (!pFile) {
 183         pFile = new TFile(fFileName);
 184         pFile->cd();
 185         printf("\n %s file opened to read RELDIS EMD events\n\n", fFileName);
 186     }
 187     fTreeNtuple = (TTree*)gDirectory->Get("h2032");
 188     fNcurrent = fStartEvent;
 189
 190     TTree *Ntu=fTreeNtuple;
 191     //
 192     // Set branch addresses
 193     // **** neutrons
 194     Ntu->SetBranchAddress("Nleft",&fNnAside);
 195     Ntu->SetBranchAddress("Eleft",&fEnAside);
 196     Ntu->SetBranchAddress("Ipdg_l_n",&fnPDGCode);
 197     Ntu->SetBranchAddress("Pxl",  fPxnAside);
 198     Ntu->SetBranchAddress("Pyl",  fPynAside);
 199     Ntu->SetBranchAddress("Pzl",  fPznAside);
 200     Ntu->SetBranchAddress("Nright",&fNnCside);
 201     Ntu->SetBranchAddress("Eright",&fEnCside);
 202     Ntu->SetBranchAddress("Pxr",   fPxnCside);
 203     Ntu->SetBranchAddress("Pyr",   fPynCside);
 204     Ntu->SetBranchAddress("Pzr",   fPznCside);
 205     // **** protons
 206     Ntu->SetBranchAddress("Nleft_p",&fNpAside);
 207     Ntu->SetBranchAddress("Etaleft_p",&fEtapAside);
 208     Ntu->SetBranchAddress("Ipdg_l_p",&fpPDGCode);
 209     Ntu->SetBranchAddress("Pxl_p",  fPxpAside);
 210     Ntu->SetBranchAddress("Pyl_p",  fPypAside);
 211     Ntu->SetBranchAddress("Pzl_p",  fPzpAside);
 212     Ntu->SetBranchAddress("Nright_p",&fNpCside);
 213     Ntu->SetBranchAddress("Etaright_p",&fEtapCside);
 214     Ntu->SetBranchAddress("Pxr_p",   fPxpCside);
 215     Ntu->SetBranchAddress("Pyr_p",   fPypCside);
 216     Ntu->SetBranchAddress("Pzr_p",   fPzpCside);
 217     // **** pi+
 218     Ntu->SetBranchAddress("Nleft_pp",&fNppAside);
 219     Ntu->SetBranchAddress("Etaleft_pp",&fEtappAside);
 220     Ntu->SetBranchAddress("Ipdg_l_pp",&fppPDGCode);
 221     Ntu->SetBranchAddress("Pxl_pp",  fPxppAside);
 222     Ntu->SetBranchAddress("Pyl_pp",  fPyppAside);
 223     Ntu->SetBranchAddress("Pzl_pp",  fPzppAside);
 224     Ntu->SetBranchAddress("Nright_pp",&fNppCside);
 225     Ntu->SetBranchAddress("Etaright_pp",&fEtappCside);
 226     Ntu->SetBranchAddress("Pxr_pp",   fPxppCside);
 227     Ntu->SetBranchAddress("Pyr_pp",   fPyppCside);
 228     Ntu->SetBranchAddress("Pzr_pp",   fPzppCside);
 229     // **** pi-
 230     Ntu->SetBranchAddress("Nleft_pm",&fNpmAside);
 231     Ntu->SetBranchAddress("Etaleft_pm",&fEtapmAside);
 232     Ntu->SetBranchAddress("Ipdg_l_pm",&fpmPDGCode);
 233     Ntu->SetBranchAddress("Pxl_pm",  fPxpmAside);
 234     Ntu->SetBranchAddress("Pyl_pm",  fPypmAside);
 235     Ntu->SetBranchAddress("Pzl_pm",  fPzpmAside);
 236     Ntu->SetBranchAddress("Nright_pm",&fNpmCside);
 237     Ntu->SetBranchAddress("Etaright_pm",&fEtapmCside);
 238     Ntu->SetBranchAddress("Pxr_pm",   fPxpmCside);
 239     Ntu->SetBranchAddress("Pyr_pm",   fPypmCside);
 240     Ntu->SetBranchAddress("Pzr_pm",   fPzpmCside);
 241     // **** pi0
 242     Ntu->SetBranchAddress("Nleft_p0",&fNp0Aside);
 243     Ntu->SetBranchAddress("Etaleft_p0",&fEtap0Aside);
 244     Ntu->SetBranchAddress("Ipdg_l_p0",&fp0PDGCode);
 245     Ntu->SetBranchAddress("Pxl_p0",  fPxp0Aside);
 246     Ntu->SetBranchAddress("Pyl_p0",  fPyp0Aside);
 247     Ntu->SetBranchAddress("Pzl_p0",  fPzp0Aside);
 248     Ntu->SetBranchAddress("Nright_p0",&fNp0Cside);
 249     Ntu->SetBranchAddress("Etaright_p0",&fEtap0Cside);
 250     Ntu->SetBranchAddress("Pxr_p0",   fPxp0Cside);
 251     Ntu->SetBranchAddress("Pyr_p0",   fPyp0Cside);
 252     Ntu->SetBranchAddress("Pzr_p0",   fPzp0Cside);
 253     // **** eta
 254     Ntu->SetBranchAddress("Nleft_et",&fNetaAside);
 255     Ntu->SetBranchAddress("Etaleft_et",&fEtaetaAside);
 256     Ntu->SetBranchAddress("Ipdg_l_et",&fetaPDGCode);
 257     Ntu->SetBranchAddress("Pxl_et",  fPxetaAside);
 258     Ntu->SetBranchAddress("Pyl_et",  fPyetaAside);
 259     Ntu->SetBranchAddress("Pzl_et",  fPzetaAside);
 260     Ntu->SetBranchAddress("Nright_et",&fNetaCside);
 261     Ntu->SetBranchAddress("Etaright_et",&fEtaetaCside);
 262     Ntu->SetBranchAddress("Pxr_et",   fPxetaCside);
 263     Ntu->SetBranchAddress("Pyr_et",   fPyetaCside);
 264     Ntu->SetBranchAddress("Pzr_et",   fPzetaCside);
 265     // **** omega
 266     Ntu->SetBranchAddress("Nleft_om",&fNomegaAside);
 267     Ntu->SetBranchAddress("Etaleft_om",&fEtaomegaAside);
 268     Ntu->SetBranchAddress("Ipdg_l_om",&fomegaPDGCode);
 269     Ntu->SetBranchAddress("Pxl_om",  fPxomegaAside);
 270     Ntu->SetBranchAddress("Pyl_om",  fPyomegaAside);
 271     Ntu->SetBranchAddress("Pzl_om",  fPzomegaAside);
 272     Ntu->SetBranchAddress("Nright_om",&fNomegaCside);
 273     Ntu->SetBranchAddress("Etaright_om",&fEtaomegaCside);
 274     Ntu->SetBranchAddress("Pxr_om",   fPxomegaCside);
 275     Ntu->SetBranchAddress("Pyr_om",   fPyomegaCside);
 276     Ntu->SetBranchAddress("Pzr_om",   fPzomegaCside);
 277 }
 278
 279 // -----------------------------------------------------------------------------------
 280 Int_t AliGenReaderEMD::NextEvent()
 281 {
 282     // Read the next event
 283     Int_t nTracks=0;
 284     fNparticle = 0; fOffset=0;
 285
 286     TFile* pFile = fTreeNtuple->GetCurrentFile();
 287     pFile->cd();
 288
 289
 290     Int_t nentries = (Int_t) fTreeNtuple->GetEntries();
 291     if(fNcurrent < nentries) {
 292         fTreeNtuple->GetEvent(fNcurrent);
 293         if(fNcurrent%100 == 0) printf("\n *** Reading event %d ***\n",fNcurrent);
 294         //
 295         if(fPcToTrack==kAll || fPcToTrack==kOnlyNucleons){ // nucleons
 296            nTracks = fNnCside+fNnAside+fNpCside+fNpAside;
 297         }
 298         if(fPcToTrack==kAll || fPcToTrack==kNotNucleons){ //pions,eta,omega
 299             nTracks += fNppCside+fNpmCside+fNppAside+fNpmAside+fNp0Aside+fNp0Cside+
 300                 fNetaAside+fNetaCside+fNomegaAside+fNomegaCside;
 301         }
 302         fNcurrent++;
 303         printf("\t #### Putting %d particles in the stack\n", nTracks);
 304         /*if(fPcToTrack==kAll || fPcToTrack==kOnlyNucleons) printf("\t\t  %d+%d neutrons, %d+%d protons\n",
 305                 fNnAside,fNnCside, fNpAside,fNpCside);
 306         if(fPcToTrack==kAll || fPcToTrack==kNotNucleons) printf("\t %d+%d pi+, %d+%d pi-, %d+%d pi0, %d+%d eta, %d+%d omega\n",
 307                 fNppAside,fNppCside,fNpmAside,fNpmCside,
 308                 fNp0Aside,fNp0Cside,fNetaAside,fNetaCside, fNomegaAside,fNomegaCside);*/
 309         return nTracks;
 310     }
 311
 312     return 0;
 313 }
 314
 315 // -----------------------------------------------------------------------------------
 316 TParticle* AliGenReaderEMD::NextParticle()
 317 {
 318     // Read the next particle
 319     Float_t p[4]={0.,0.,0.,0.};
 320     int pdgCode=0;
 321
 322     if(fPcToTrack==kAll || fPcToTrack==kOnlyNucleons){//***********************************************
 323       if(fNparticle<fNnAside){
 324         p[0] = fPxnAside[fNparticle];
 325         p[1] = fPynAside[fNparticle];
 326         p[2] = fPznAside[fNparticle];
 327         pdgCode = fnPDGCode;
 328 //    printf(" pc%d n sideA: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
 329       }
 330       else if(fNparticle>=fNnAside && fNparticle<(fNnAside+fNnCside)){
 331         p[0] = fPxnCside[fNparticle];
 332         p[1] = fPynCside[fNparticle];
 333         p[2] = fPznCside[fNparticle];
 334         pdgCode = fnPDGCode;
 335 //    printf(" pc%d n sideC: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
 336       }
 337       else if(fNparticle>=fNnAside+fNnCside && fNparticle<(fNnAside+fNnCside+fNpAside)){
 338         p[0] = fPxpAside[fNparticle];
 339         p[1] = fPypAside[fNparticle];
 340         p[2] = fPzpAside[fNparticle];
 341         pdgCode = fpPDGCode;
 342 //    printf(" pc%d p sideA: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
 343       }
 344       else if(fNparticle>=fNnAside+fNnCside+fNpAside && fNparticle<(fNnAside+fNnCside+fNpCside+fNpAside)){
 345         p[0] = fPxpCside[fNparticle];
 346         p[1] = fPypCside[fNparticle];
 347         p[2] = fPzpCside[fNparticle];
 348         pdgCode = fpPDGCode;
 349 //    printf(" pc%d p sideC: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
 350       }
 351       fOffset = fNnAside+fNnCside+fNpCside+fNpAside;
 352     } //**********************************************************************************************
 353     if(fPcToTrack==kAll || fPcToTrack==kNotNucleons){
 354       if(fNparticle>=fOffset && fNparticle<fOffset+fNppAside){ // *** pi +
 355         p[0] = fPxppAside[fNparticle];
 356         p[1] = fPyppAside[fNparticle];
 357         p[2] = fPzppAside[fNparticle];
 358         pdgCode = fppPDGCode;
 359 //    printf(" pc%d pi+ sideA: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
 360       }
 361       if(fNparticle>=fOffset+fNppAside && fNparticle<fOffset+fNppAside+fNppCside){
 362         p[0] = fPxppCside[fNparticle];
 363         p[1] = fPyppCside[fNparticle];
 364         p[2] = fPzppCside[fNparticle];
 365         pdgCode = fppPDGCode;
 366 //    printf(" pc%d pi+ sideC: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
 367       }
 368       if(fNparticle>=fOffset+fNppAside+fNppCside && fNparticle<fOffset+fNppAside+fNppCside+fNpmAside){ // *** pi -
 369         p[0] = fPxpmAside[fNparticle];
 370         p[1] = fPypmAside[fNparticle];
 371         p[2] = fPzpmAside[fNparticle];
 372         pdgCode = fpmPDGCode;
 373 //    printf(" pc%d pi- sideA: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
 374       }
 375       if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside && fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside){
 376         p[0] = fPxpmCside[fNparticle];
 377         p[1] = fPypmCside[fNparticle];
 378         p[2] = fPzpmCside[fNparticle];
 379         pdgCode = fpmPDGCode;
 380 //    printf(" pc%d pi- sideC: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
 381       }
 382       if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside &&
 383          fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside){ // *** pi 0
 384         p[0] = fPxp0Aside[fNparticle];
 385         p[1] = fPyp0Aside[fNparticle];
 386         p[2] = fPzp0Aside[fNparticle];
 387         pdgCode = fp0PDGCode;
 388 //    printf(" pc%d pi0 sideA: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
 389       }
 390       if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside &&
 391         fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside){
 392         p[0] = fPxp0Cside[fNparticle];
 393         p[1] = fPyp0Cside[fNparticle];
 394         p[2] = fPzp0Cside[fNparticle];
 395         pdgCode = fp0PDGCode;
 396 //    printf(" pc%d pi0 sideC: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
 397       }
 398       if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside &&
 399          fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside){ // *** eta
 400         p[0] = fPxetaAside[fNparticle];
 401         p[1] = fPyetaAside[fNparticle];
 402         p[2] = fPzetaAside[fNparticle];
 403         pdgCode = fetaPDGCode;
 404 //    printf(" pc%d eta sideA: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
 405       }
 406       if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside &&
 407          fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside+fNetaCside){
 408         p[0] = fPxetaCside[fNparticle];
 409         p[1] = fPyetaCside[fNparticle];
 410         p[2] = fPzetaCside[fNparticle];
 411         pdgCode = fetaPDGCode;
 412 //    printf(" pc%d eta sideC: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
 413       }
 414       if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside+fNetaCside &&
 415          fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside+fNetaCside+fNomegaAside){ // *** omega
 416         p[0] = fPxomegaAside[fNparticle];
 417         p[1] = fPyomegaAside[fNparticle];
 418         p[2] = fPzomegaAside[fNparticle];
 419         pdgCode = fomegaPDGCode;
 420 //    printf(" pc%d omega sideA: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
 421       }
 422       if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside+fNetaCside+fNomegaAside
 423       && fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside+fNetaCside+fNomegaAside+fNomegaCside){
 424         p[0] = fPxomegaCside[fNparticle];
 425         p[1] = fPyomegaCside[fNparticle];
 426         p[2] = fPzomegaCside[fNparticle];
 427         pdgCode = fomegaPDGCode;
 428 //    printf(" pc%d omega sideC: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
 429       }
 430
 431     }
 432
 433     Float_t ptot = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
 434     Double_t amass = TDatabasePDG::Instance()->GetParticle(pdgCode)->Mass();
 435     p[3] = TMath::Sqrt(ptot*ptot+amass*amass);
 436
 437     if(p[3]<=amass){
 438        Warning("Generate","Particle %d  E = %f GeV mass = %f GeV ",pdgCode,p[3],amass);
 439     }
 440
 441     //printf("  Pc %d:  PDGcode %d  p(%1.2f, %1.2f, %1.2f, %1.3f)\n",
 442     //  fNparticle,pdgCode,p[0], p[1], p[2], p[3]);
 443
 444     TParticle* particle = new TParticle(pdgCode, 0, -1, -1, -1, -1,
 445         p[0], p[1], p[2], p[3], 0., 0., 0., 0.);
 446     if((p[0]*p[0]+p[1]*p[1]+p[2]*p[2])>1e-5) particle->SetBit(kTransportBit);
 447     fNparticle++;
 448     return particle;
 449 }