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     for(int i=0; i<70; i++){
 151        fPxnAside[i] = fPynAside[i] = fPznAside[i] = 0.;
 152        fPxnCside[i] = fPynCside[i] = fPznCside[i] = 0.;
 153        if(i<50){
 154          fPxpAside[i] = fPypAside[i] = fPzpAside[i] = 0.;
 155          fPxpCside[i] = fPypCside[i] = fPzpCside[i] = 0.;
 156          if(i<30){
 157            fPxppAside[i] = fPyppAside[i] = fPzppAside[i] = 0.;
 158            fPxppCside[i] = fPyppCside[i] = fPzppCside[i] = 0.;
 159            fPxpmAside[i] = fPypmAside[i] = fPzpmAside[i] = 0.;
 160            fPxpmCside[i] = fPypmCside[i] = fPzpmCside[i] = 0.;
 161            fPxp0Aside[i] = fPyp0Aside[i] = fPzp0Aside[i] = 0.;
 162            fPxp0Cside[i] = fPyp0Cside[i] = fPzp0Cside[i] = 0.;
 163            if(i<15){
 164              fPxetaAside[i] = fPyetaAside[i] = fPzetaAside[i] = 0.;
 165              fPxetaCside[i] = fPyetaCside[i] = fPzetaCside[i] = 0.;
 166              fPxomegaAside[i] = fPyomegaAside[i] = fPzomegaAside[i] = 0.;
 167              fPxomegaCside[i] = fPyomegaCside[i] = fPzomegaCside[i] = 0.;
 168            }
 169          }
 170        }
 171     }
 172     reader.Copy(*this);
 173 }
 174   // -----------------------------------------------------------------------------------
 175 AliGenReaderEMD::~AliGenReaderEMD()
 176 {
 177     delete fTreeNtuple;
 178 }
 179
 180 // -----------------------------------------------------------------------------------
 181 AliGenReaderEMD& AliGenReaderEMD::operator=(const  AliGenReaderEMD& rhs)
 182 {
 183 // Assignment operator
 184     rhs.Copy(*this);
 185     return *this;
 186 }
 187
 188 // -----------------------------------------------------------------------------------
 189 void AliGenReaderEMD::Copy(TObject&) const
 190 {
 191     //
 192     // Copy
 193     //
 194     Fatal("Copy","Not implemented!\n");
 195 }
 196
 197 // -----------------------------------------------------------------------------------
 198 void AliGenReaderEMD::Init()
 199 {
 200 //
 201 // Reset the existing file environment and open a new root file
 202
 203     TFile *pFile=0;
 204     if (!pFile) {
 205         pFile = new TFile(fFileName);
 206         pFile->cd();
 207         printf("\n %s file opened to read RELDIS EMD events\n\n", fFileName);
 208     }
 209     fTreeNtuple = (TTree*)gDirectory->Get("h2032");
 210     fNcurrent = fStartEvent;
 211
 212     TTree *Ntu=fTreeNtuple;
 213     //
 214     // Set branch addresses
 215     // **** neutrons
 216     Ntu->SetBranchAddress("Nleft",&fNnAside);
 217     Ntu->SetBranchAddress("Eleft",&fEnAside);
 218     Ntu->SetBranchAddress("Ipdg_l_n",&fnPDGCode);
 219     Ntu->SetBranchAddress("Pxl",  fPxnAside);
 220     Ntu->SetBranchAddress("Pyl",  fPynAside);
 221     Ntu->SetBranchAddress("Pzl",  fPznAside);
 222     Ntu->SetBranchAddress("Nright",&fNnCside);
 223     Ntu->SetBranchAddress("Eright",&fEnCside);
 224     Ntu->SetBranchAddress("Pxr",   fPxnCside);
 225     Ntu->SetBranchAddress("Pyr",   fPynCside);
 226     Ntu->SetBranchAddress("Pzr",   fPznCside);
 227     // **** protons
 228     Ntu->SetBranchAddress("Nleft_p",&fNpAside);
 229     Ntu->SetBranchAddress("Etaleft_p",&fEtapAside);
 230     Ntu->SetBranchAddress("Ipdg_l_p",&fpPDGCode);
 231     Ntu->SetBranchAddress("Pxl_p",  fPxpAside);
 232     Ntu->SetBranchAddress("Pyl_p",  fPypAside);
 233     Ntu->SetBranchAddress("Pzl_p",  fPzpAside);
 234     Ntu->SetBranchAddress("Nright_p",&fNpCside);
 235     Ntu->SetBranchAddress("Etaright_p",&fEtapCside);
 236     Ntu->SetBranchAddress("Pxr_p",   fPxpCside);
 237     Ntu->SetBranchAddress("Pyr_p",   fPypCside);
 238     Ntu->SetBranchAddress("Pzr_p",   fPzpCside);
 239     // **** pi+
 240     Ntu->SetBranchAddress("Nleft_pp",&fNppAside);
 241     Ntu->SetBranchAddress("Etaleft_pp",&fEtappAside);
 242     Ntu->SetBranchAddress("Ipdg_l_pp",&fppPDGCode);
 243     Ntu->SetBranchAddress("Pxl_pp",  fPxppAside);
 244     Ntu->SetBranchAddress("Pyl_pp",  fPyppAside);
 245     Ntu->SetBranchAddress("Pzl_pp",  fPzppAside);
 246     Ntu->SetBranchAddress("Nright_pp",&fNppCside);
 247     Ntu->SetBranchAddress("Etaright_pp",&fEtappCside);
 248     Ntu->SetBranchAddress("Pxr_pp",   fPxppCside);
 249     Ntu->SetBranchAddress("Pyr_pp",   fPyppCside);
 250     Ntu->SetBranchAddress("Pzr_pp",   fPzppCside);
 251     // **** pi-
 252     Ntu->SetBranchAddress("Nleft_pm",&fNpmAside);
 253     Ntu->SetBranchAddress("Etaleft_pm",&fEtapmAside);
 254     Ntu->SetBranchAddress("Ipdg_l_pm",&fpmPDGCode);
 255     Ntu->SetBranchAddress("Pxl_pm",  fPxpmAside);
 256     Ntu->SetBranchAddress("Pyl_pm",  fPypmAside);
 257     Ntu->SetBranchAddress("Pzl_pm",  fPzpmAside);
 258     Ntu->SetBranchAddress("Nright_pm",&fNpmCside);
 259     Ntu->SetBranchAddress("Etaright_pm",&fEtapmCside);
 260     Ntu->SetBranchAddress("Pxr_pm",   fPxpmCside);
 261     Ntu->SetBranchAddress("Pyr_pm",   fPypmCside);
 262     Ntu->SetBranchAddress("Pzr_pm",   fPzpmCside);
 263     // **** pi0
 264     Ntu->SetBranchAddress("Nleft_p0",&fNp0Aside);
 265     Ntu->SetBranchAddress("Etaleft_p0",&fEtap0Aside);
 266     Ntu->SetBranchAddress("Ipdg_l_p0",&fp0PDGCode);
 267     Ntu->SetBranchAddress("Pxl_p0",  fPxp0Aside);
 268     Ntu->SetBranchAddress("Pyl_p0",  fPyp0Aside);
 269     Ntu->SetBranchAddress("Pzl_p0",  fPzp0Aside);
 270     Ntu->SetBranchAddress("Nright_p0",&fNp0Cside);
 271     Ntu->SetBranchAddress("Etaright_p0",&fEtap0Cside);
 272     Ntu->SetBranchAddress("Pxr_p0",   fPxp0Cside);
 273     Ntu->SetBranchAddress("Pyr_p0",   fPyp0Cside);
 274     Ntu->SetBranchAddress("Pzr_p0",   fPzp0Cside);
 275     // **** eta
 276     Ntu->SetBranchAddress("Nleft_et",&fNetaAside);
 277     Ntu->SetBranchAddress("Etaleft_et",&fEtaetaAside);
 278     Ntu->SetBranchAddress("Ipdg_l_et",&fetaPDGCode);
 279     Ntu->SetBranchAddress("Pxl_et",  fPxetaAside);
 280     Ntu->SetBranchAddress("Pyl_et",  fPyetaAside);
 281     Ntu->SetBranchAddress("Pzl_et",  fPzetaAside);
 282     Ntu->SetBranchAddress("Nright_et",&fNetaCside);
 283     Ntu->SetBranchAddress("Etaright_et",&fEtaetaCside);
 284     Ntu->SetBranchAddress("Pxr_et",   fPxetaCside);
 285     Ntu->SetBranchAddress("Pyr_et",   fPyetaCside);
 286     Ntu->SetBranchAddress("Pzr_et",   fPzetaCside);
 287     // **** omega
 288     Ntu->SetBranchAddress("Nleft_om",&fNomegaAside);
 289     Ntu->SetBranchAddress("Etaleft_om",&fEtaomegaAside);
 290     Ntu->SetBranchAddress("Ipdg_l_om",&fomegaPDGCode);
 291     Ntu->SetBranchAddress("Pxl_om",  fPxomegaAside);
 292     Ntu->SetBranchAddress("Pyl_om",  fPyomegaAside);
 293     Ntu->SetBranchAddress("Pzl_om",  fPzomegaAside);
 294     Ntu->SetBranchAddress("Nright_om",&fNomegaCside);
 295     Ntu->SetBranchAddress("Etaright_om",&fEtaomegaCside);
 296     Ntu->SetBranchAddress("Pxr_om",   fPxomegaCside);
 297     Ntu->SetBranchAddress("Pyr_om",   fPyomegaCside);
 298     Ntu->SetBranchAddress("Pzr_om",   fPzomegaCside);
 299 }
 300
 301 // -----------------------------------------------------------------------------------
 302 Int_t AliGenReaderEMD::NextEvent()
 303 {
 304     // Read the next event
 305     Int_t nTracks=0;
 306     fNparticle = 0; fOffset=0;
 307
 308     TFile* pFile = fTreeNtuple->GetCurrentFile();
 309     pFile->cd();
 310
 311
 312     Int_t nentries = (Int_t) fTreeNtuple->GetEntries();
 313     if(fNcurrent < nentries) {
 314         fTreeNtuple->GetEvent(fNcurrent);
 315         if(fNcurrent%100 == 0) printf("\n *** Reading event %d ***\n",fNcurrent);
 316         //
 317         if(fPcToTrack==kAll || fPcToTrack==kOnlyNucleons){ // nucleons
 318            nTracks = fNnCside+fNnAside+fNpCside+fNpAside;
 319         }
 320         if(fPcToTrack==kAll || fPcToTrack==kNotNucleons){ //pions,eta,omega
 321             nTracks += fNppCside+fNpmCside+fNppAside+fNpmAside+fNp0Aside+fNp0Cside+
 322                 fNetaAside+fNetaCside+fNomegaAside+fNomegaCside;
 323         }
 324         fNcurrent++;
 325         printf("\t #### Putting %d particles in the stack\n", nTracks);
 326         /*if(fPcToTrack==kAll || fPcToTrack==kOnlyNucleons) printf("\t\t  %d+%d neutrons, %d+%d protons\n",
 327                 fNnAside,fNnCside, fNpAside,fNpCside);
 328         if(fPcToTrack==kAll || fPcToTrack==kNotNucleons) printf("\t %d+%d pi+, %d+%d pi-, %d+%d pi0, %d+%d eta, %d+%d omega\n",
 329                 fNppAside,fNppCside,fNpmAside,fNpmCside,
 330                 fNp0Aside,fNp0Cside,fNetaAside,fNetaCside, fNomegaAside,fNomegaCside);*/
 331         return nTracks;
 332     }
 333
 334     return 0;
 335 }
 336
 337 // -----------------------------------------------------------------------------------
 338 TParticle* AliGenReaderEMD::NextParticle()
 339 {
 340     // Read the next particle
 341     Float_t p[4]={0.,0.,0.,0.};
 342     int pdgCode=0;
 343
 344     if(fPcToTrack==kAll || fPcToTrack==kOnlyNucleons){//***********************************************
 345       if(fNparticle<fNnAside){
 346         p[0] = fPxnAside[fNparticle];
 347         p[1] = fPynAside[fNparticle];
 348         p[2] = fPznAside[fNparticle];
 349         pdgCode = fnPDGCode;
 350 //    printf(" pc%d n sideA: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
 351       }
 352       else if(fNparticle>=fNnAside && fNparticle<(fNnAside+fNnCside)){
 353         p[0] = fPxnCside[fNparticle];
 354         p[1] = fPynCside[fNparticle];
 355         p[2] = fPznCside[fNparticle];
 356         pdgCode = fnPDGCode;
 357 //    printf(" pc%d n sideC: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
 358       }
 359       else if(fNparticle>=fNnAside+fNnCside && fNparticle<(fNnAside+fNnCside+fNpAside)){
 360         p[0] = fPxpAside[fNparticle];
 361         p[1] = fPypAside[fNparticle];
 362         p[2] = fPzpAside[fNparticle];
 363         pdgCode = fpPDGCode;
 364 //    printf(" pc%d p sideA: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
 365       }
 366       else if(fNparticle>=fNnAside+fNnCside+fNpAside && fNparticle<(fNnAside+fNnCside+fNpCside+fNpAside)){
 367         p[0] = fPxpCside[fNparticle];
 368         p[1] = fPypCside[fNparticle];
 369         p[2] = fPzpCside[fNparticle];
 370         pdgCode = fpPDGCode;
 371 //    printf(" pc%d p sideC: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
 372       }
 373       fOffset = fNnAside+fNnCside+fNpCside+fNpAside;
 374     } //**********************************************************************************************
 375     if(fPcToTrack==kAll || fPcToTrack==kNotNucleons){
 376       if(fNparticle>=fOffset && fNparticle<fOffset+fNppAside){ // *** pi +
 377         p[0] = fPxppAside[fNparticle];
 378         p[1] = fPyppAside[fNparticle];
 379         p[2] = fPzppAside[fNparticle];
 380         pdgCode = fppPDGCode;
 381 //    printf(" pc%d pi+ sideA: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
 382       }
 383       if(fNparticle>=fOffset+fNppAside && fNparticle<fOffset+fNppAside+fNppCside){
 384         p[0] = fPxppCside[fNparticle];
 385         p[1] = fPyppCside[fNparticle];
 386         p[2] = fPzppCside[fNparticle];
 387         pdgCode = fppPDGCode;
 388 //    printf(" pc%d pi+ sideC: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
 389       }
 390       if(fNparticle>=fOffset+fNppAside+fNppCside && fNparticle<fOffset+fNppAside+fNppCside+fNpmAside){ // *** pi -
 391         p[0] = fPxpmAside[fNparticle];
 392         p[1] = fPypmAside[fNparticle];
 393         p[2] = fPzpmAside[fNparticle];
 394         pdgCode = fpmPDGCode;
 395 //    printf(" pc%d pi- sideA: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
 396       }
 397       if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside && fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside){
 398         p[0] = fPxpmCside[fNparticle];
 399         p[1] = fPypmCside[fNparticle];
 400         p[2] = fPzpmCside[fNparticle];
 401         pdgCode = fpmPDGCode;
 402 //    printf(" pc%d pi- sideC: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
 403       }
 404       if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside &&
 405          fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside){ // *** pi 0
 406         p[0] = fPxp0Aside[fNparticle];
 407         p[1] = fPyp0Aside[fNparticle];
 408         p[2] = fPzp0Aside[fNparticle];
 409         pdgCode = fp0PDGCode;
 410 //    printf(" pc%d pi0 sideA: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
 411       }
 412       if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside &&
 413         fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside){
 414         p[0] = fPxp0Cside[fNparticle];
 415         p[1] = fPyp0Cside[fNparticle];
 416         p[2] = fPzp0Cside[fNparticle];
 417         pdgCode = fp0PDGCode;
 418 //    printf(" pc%d pi0 sideC: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
 419       }
 420       if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside &&
 421          fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside){ // *** eta
 422         p[0] = fPxetaAside[fNparticle];
 423         p[1] = fPyetaAside[fNparticle];
 424         p[2] = fPzetaAside[fNparticle];
 425         pdgCode = fetaPDGCode;
 426 //    printf(" pc%d eta sideA: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
 427       }
 428       if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside &&
 429          fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside+fNetaCside){
 430         p[0] = fPxetaCside[fNparticle];
 431         p[1] = fPyetaCside[fNparticle];
 432         p[2] = fPzetaCside[fNparticle];
 433         pdgCode = fetaPDGCode;
 434 //    printf(" pc%d eta sideC: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
 435       }
 436       if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside+fNetaCside &&
 437          fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside+fNetaCside+fNomegaAside){ // *** omega
 438         p[0] = fPxomegaAside[fNparticle];
 439         p[1] = fPyomegaAside[fNparticle];
 440         p[2] = fPzomegaAside[fNparticle];
 441         pdgCode = fomegaPDGCode;
 442 //    printf(" pc%d omega sideA: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
 443       }
 444       if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside+fNetaCside+fNomegaAside
 445       && fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside+fNetaCside+fNomegaAside+fNomegaCside){
 446         p[0] = fPxomegaCside[fNparticle];
 447         p[1] = fPyomegaCside[fNparticle];
 448         p[2] = fPzomegaCside[fNparticle];
 449         pdgCode = fomegaPDGCode;
 450 //    printf(" pc%d omega sideC: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
 451       }
 452
 453     }
 454
 455     Float_t ptot = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
 456     Double_t amass = TDatabasePDG::Instance()->GetParticle(pdgCode)->Mass();
 457     p[3] = TMath::Sqrt(ptot*ptot+amass*amass);
 458
 459     if(p[3]<=amass){
 460        Warning("Generate","Particle %d  E = %f GeV mass = %f GeV ",pdgCode,p[3],amass);
 461     }
 462
 463     //printf("  Pc %d:  PDGcode %d  p(%1.2f, %1.2f, %1.2f, %1.3f)\n",
 464     //  fNparticle,pdgCode,p[0], p[1], p[2], p[3]);
 465
 466     TParticle* particle = new TParticle(pdgCode, 0, -1, -1, -1, -1,
 467         p[0], p[1], p[2], p[3], 0., 0., 0., 0.);
 468     if((p[0]*p[0]+p[1]*p[1]+p[2]*p[2])>1e-5) particle->SetBit(kTransportBit);
 469     fNparticle++;
 470     return particle;
 471 }