Updates.
[u/mrichter/AliRoot.git] / EVGEN / AliGenReaderEMD.cxx
index f6c7a87..8c8ec97 100644 (file)
@@ -35,36 +35,116 @@ AliGenReaderEMD::AliGenReaderEMD():
     fNcurrent(0),  
     fNparticle(0), 
     fTreeNtuple(0),
-    fIPSide(0),
     fPcToTrack(0),
+    fOffset(0),
     fNnAside(0),
     fEnAside(0),
+    fnPDGCode(0),
     fNnCside(0),
     fEnCside(0),
     fNpAside(0),
     fEtapAside(0),
+    fpPDGCode(0),
     fNpCside(0),
-    fEtapCside(0)
+    fEtapCside(0),
+    fNppAside(0),
+    fEtappAside(0),
+    fppPDGCode(0),
+    fNppCside(0),
+    fEtappCside(0),
+    fNpmAside(0),
+    fEtapmAside(0),
+    fpmPDGCode(0),
+    fNpmCside(0),
+    fEtapmCside(0),
+    fNp0Aside(0),
+    fEtap0Aside(0),
+    fp0PDGCode(0),
+    fNp0Cside(0),
+    fEtap0Cside(0),
+    fNetaAside(0),
+    fEtaetaAside(0),
+    fetaPDGCode(0),
+    fNetaCside(0),
+    fEtaetaCside(0),
+    fNomegaAside(0),
+    fEtaomegaAside(0),
+    fomegaPDGCode(0),
+    fNomegaCside(0),
+    fEtaomegaCside(0)
 {
-// Default constructor
+// Std constructor
+    for(int i=0; i<70; i++){
+       fPxnAside[i] = fPynAside[i] = fPznAside[i] = 0.;
+       fPxnCside[i] = fPynCside[i] = fPznCside[i] = 0.;
+       if(i<50){
+         fPxpAside[i] = fPypAside[i] = fPzpAside[i] = 0.;
+        fPxpCside[i] = fPypCside[i] = fPzpCside[i] = 0.;
+         if(i<30){
+           fPxppAside[i] = fPyppAside[i] = fPzppAside[i] = 0.;
+          fPxppCside[i] = fPyppCside[i] = fPzppCside[i] = 0.;
+           fPxpmAside[i] = fPypmAside[i] = fPzpmAside[i] = 0.;
+          fPxpmCside[i] = fPypmCside[i] = fPzpmCside[i] = 0.;
+           fPxp0Aside[i] = fPyp0Aside[i] = fPzp0Aside[i] = 0.;
+          fPxp0Cside[i] = fPyp0Cside[i] = fPzp0Cside[i] = 0.;
+          if(i<15){
+             fPxetaAside[i] = fPyetaAside[i] = fPzetaAside[i] = 0.;
+            fPxetaCside[i] = fPyetaCside[i] = fPzetaCside[i] = 0.;
+             fPxomegaAside[i] = fPyomegaAside[i] = fPzomegaAside[i] = 0.;
+            fPxomegaCside[i] = fPyomegaCside[i] = fPzomegaCside[i] = 0.;
+          }
+        }
+       }       
+    }
+    if(fPcToTrack==kAll) printf("\n\t   *** AliGenReaderEMD will track all produced particles \n\n");
+    else if(fPcToTrack==kNotNucleons) printf("\n\t   *** AliGenReaderEMD will track all produced particles except nucleons\n\n");
+    else if(fPcToTrack==kOnlyNucleons) printf("\n\t   *** AliGenReaderEMD will track only nucleons\n\n");
 }
 
+
 AliGenReaderEMD::AliGenReaderEMD(const AliGenReaderEMD &reader):
     AliGenReader(reader),
     fStartEvent(0),
     fNcurrent(0),  
     fNparticle(0), 
     fTreeNtuple(0),
-    fIPSide(0),
     fPcToTrack(0),
+    fOffset(0),
     fNnAside(0),
     fEnAside(0),
+    fnPDGCode(0),
     fNnCside(0),
     fEnCside(0),
     fNpAside(0),
     fEtapAside(0),
+    fpPDGCode(0),
     fNpCside(0),
-    fEtapCside(0)
+    fEtapCside(0),
+    fNppAside(0),
+    fEtappAside(0),
+    fppPDGCode(0),
+    fNppCside(0),
+    fEtappCside(0),
+    fNpmAside(0),
+    fEtapmAside(0),
+    fpmPDGCode(0),
+    fNpmCside(0),
+    fEtapmCside(0),
+    fNp0Aside(0),
+    fEtap0Aside(0),
+    fp0PDGCode(0),
+    fNp0Cside(0),
+    fEtap0Cside(0),
+    fNetaAside(0),
+    fEtaetaAside(0),
+    fetaPDGCode(0),
+    fNetaCside(0),
+    fEtaetaCside(0),
+    fNomegaAside(0),
+    fEtaomegaAside(0),
+    fomegaPDGCode(0),
+    fNomegaCside(0),
+    fEtaomegaCside(0)
 {
     // Copy Constructor
     reader.Copy(*this);
@@ -102,7 +182,7 @@ void AliGenReaderEMD::Init()
     if (!pFile) {
        pFile = new TFile(fFileName);
        pFile->cd();
-       printf("\n %s file opened\n\n", fFileName);
+       printf("\n %s file opened to read RELDIS EMD events\n\n", fFileName);
     }
     fTreeNtuple = (TTree*)gDirectory->Get("h2032");
     fNcurrent = fStartEvent;
@@ -113,6 +193,7 @@ void AliGenReaderEMD::Init()
     // **** neutrons
     Ntu->SetBranchAddress("Nleft",&fNnAside);
     Ntu->SetBranchAddress("Eleft",&fEnAside);
+    Ntu->SetBranchAddress("Ipdg_l_n",&fnPDGCode);
     Ntu->SetBranchAddress("Pxl",  fPxnAside);
     Ntu->SetBranchAddress("Pyl",  fPynAside);
     Ntu->SetBranchAddress("Pzl",  fPznAside);
@@ -124,6 +205,7 @@ void AliGenReaderEMD::Init()
     // **** protons
     Ntu->SetBranchAddress("Nleft_p",&fNpAside);
     Ntu->SetBranchAddress("Etaleft_p",&fEtapAside);
+    Ntu->SetBranchAddress("Ipdg_l_p",&fpPDGCode);
     Ntu->SetBranchAddress("Pxl_p",  fPxpAside);
     Ntu->SetBranchAddress("Pyl_p",  fPypAside);
     Ntu->SetBranchAddress("Pzl_p",  fPzpAside);
@@ -132,6 +214,66 @@ void AliGenReaderEMD::Init()
     Ntu->SetBranchAddress("Pxr_p",   fPxpCside);
     Ntu->SetBranchAddress("Pyr_p",   fPypCside);
     Ntu->SetBranchAddress("Pzr_p",   fPzpCside);
+    // **** pi+
+    Ntu->SetBranchAddress("Nleft_pp",&fNppAside);
+    Ntu->SetBranchAddress("Etaleft_pp",&fEtappAside);
+    Ntu->SetBranchAddress("Ipdg_l_pp",&fppPDGCode);
+    Ntu->SetBranchAddress("Pxl_pp",  fPxppAside);
+    Ntu->SetBranchAddress("Pyl_pp",  fPyppAside);
+    Ntu->SetBranchAddress("Pzl_pp",  fPzppAside);
+    Ntu->SetBranchAddress("Nright_pp",&fNppCside);
+    Ntu->SetBranchAddress("Etaright_pp",&fEtappCside);
+    Ntu->SetBranchAddress("Pxr_pp",   fPxppCside);
+    Ntu->SetBranchAddress("Pyr_pp",   fPyppCside);
+    Ntu->SetBranchAddress("Pzr_pp",   fPzppCside);
+    // **** pi-
+    Ntu->SetBranchAddress("Nleft_pm",&fNpmAside);
+    Ntu->SetBranchAddress("Etaleft_pm",&fEtapmAside);
+    Ntu->SetBranchAddress("Ipdg_l_pm",&fpmPDGCode);
+    Ntu->SetBranchAddress("Pxl_pm",  fPxpmAside);
+    Ntu->SetBranchAddress("Pyl_pm",  fPypmAside);
+    Ntu->SetBranchAddress("Pzl_pm",  fPzpmAside);
+    Ntu->SetBranchAddress("Nright_pm",&fNpmCside);
+    Ntu->SetBranchAddress("Etaright_pm",&fEtapmCside);
+    Ntu->SetBranchAddress("Pxr_pm",   fPxpmCside);
+    Ntu->SetBranchAddress("Pyr_pm",   fPypmCside);
+    Ntu->SetBranchAddress("Pzr_pm",   fPzpmCside);
+    // **** pi0
+    Ntu->SetBranchAddress("Nleft_p0",&fNp0Aside);
+    Ntu->SetBranchAddress("Etaleft_p0",&fEtap0Aside);
+    Ntu->SetBranchAddress("Ipdg_l_p0",&fp0PDGCode);
+    Ntu->SetBranchAddress("Pxl_p0",  fPxp0Aside);
+    Ntu->SetBranchAddress("Pyl_p0",  fPyp0Aside);
+    Ntu->SetBranchAddress("Pzl_p0",  fPzp0Aside);
+    Ntu->SetBranchAddress("Nright_p0",&fNp0Cside);
+    Ntu->SetBranchAddress("Etaright_p0",&fEtap0Cside);
+    Ntu->SetBranchAddress("Pxr_p0",   fPxp0Cside);
+    Ntu->SetBranchAddress("Pyr_p0",   fPyp0Cside);
+    Ntu->SetBranchAddress("Pzr_p0",   fPzp0Cside);
+    // **** eta
+    Ntu->SetBranchAddress("Nleft_et",&fNetaAside);
+    Ntu->SetBranchAddress("Etaleft_et",&fEtaetaAside);
+    Ntu->SetBranchAddress("Ipdg_l_et",&fetaPDGCode);
+    Ntu->SetBranchAddress("Pxl_et",  fPxetaAside);
+    Ntu->SetBranchAddress("Pyl_et",  fPyetaAside);
+    Ntu->SetBranchAddress("Pzl_et",  fPzetaAside);
+    Ntu->SetBranchAddress("Nright_et",&fNetaCside);
+    Ntu->SetBranchAddress("Etaright_et",&fEtaetaCside);
+    Ntu->SetBranchAddress("Pxr_et",   fPxetaCside);
+    Ntu->SetBranchAddress("Pyr_et",   fPyetaCside);
+    Ntu->SetBranchAddress("Pzr_et",   fPzetaCside);
+    // **** omega
+    Ntu->SetBranchAddress("Nleft_om",&fNomegaAside);
+    Ntu->SetBranchAddress("Etaleft_om",&fEtaomegaAside);
+    Ntu->SetBranchAddress("Ipdg_l_om",&fomegaPDGCode);
+    Ntu->SetBranchAddress("Pxl_om",  fPxomegaAside);
+    Ntu->SetBranchAddress("Pyl_om",  fPyomegaAside);
+    Ntu->SetBranchAddress("Pzl_om",  fPzomegaAside);
+    Ntu->SetBranchAddress("Nright_om",&fNomegaCside);
+    Ntu->SetBranchAddress("Etaright_om",&fEtaomegaCside);
+    Ntu->SetBranchAddress("Pxr_om",   fPxomegaCside);
+    Ntu->SetBranchAddress("Pyr_om",   fPyomegaCside);
+    Ntu->SetBranchAddress("Pzr_om",   fPzomegaCside);
 }
 
 // -----------------------------------------------------------------------------------
@@ -139,7 +281,7 @@ Int_t AliGenReaderEMD::NextEvent()
 {
     // Read the next event  
     Int_t nTracks=0;
-    fNparticle = 0;
+    fNparticle = 0; fOffset=0;
     
     TFile* pFile = fTreeNtuple->GetCurrentFile();
     pFile->cd();
@@ -148,41 +290,22 @@ Int_t AliGenReaderEMD::NextEvent()
     Int_t nentries = (Int_t) fTreeNtuple->GetEntries();
     if(fNcurrent < nentries) {
        fTreeNtuple->GetEvent(fNcurrent);
-       printf("\n *** Event %d read ***\n",fNcurrent);
-       //printf("\t A side-> %d neutrons and %d protons emitted\n", fNnAside, fNpAside);
-       //printf("\t C side-> %d neutrons and %d protons emitted\n\n", fNnCside, fNpCside);
+       if(fNcurrent%100 == 0) printf("\n *** Reading event %d ***\n",fNcurrent);
        //
-       // **** IPSIde = 0->both sides, 1->only Cside, 2->only A side  
-       // #### fPcToTrack = 0->neutrons, 1->protons
-       if(fPcToTrack==0){ // neutrons
-         if(fIPSide==0){
-           printf("\t Tracking %d+%d neutrons emitted on C+A side\n", fNnCside,fNnAside);
-           nTracks = fNnCside+fNnAside;
-         }
-         else if(fIPSide==1){ 
-           printf("\t Tracking %d neutrons emitted on C side\n", fNnCside);
-           nTracks    = fNnCside;
-         }
-         else if(fIPSide==2){
-           printf("\t Tracking %d neutrons emitted on A side\n", fNnAside);
-           nTracks    = fNnAside;
-         }
+       if(fPcToTrack==kAll || fPcToTrack==kOnlyNucleons){ // nucleons
+          nTracks = fNnCside+fNnAside+fNpCside+fNpAside;
        }
-       else if(fPcToTrack==1){ //protons
-         if(fIPSide==0){
-           printf("\t Tracking %d+%d protons emitted on C+A side\n", fNpCside,fNpAside);
-           nTracks    = fNpCside+fNpAside;
-         }
-         else if(fIPSide==1){
-           printf("\t Tracking %d protons emitted on C side\n", fNpCside);
-           nTracks    = fNpCside;
-         }
-         else if(fIPSide==2){
-           printf("\t Tracking %d protons emitted on A side\n", fNpAside);
-           nTracks    = fNpAside;
-         }
+       if(fPcToTrack==kAll || fPcToTrack==kNotNucleons){ //pions,eta,omega
+           nTracks += fNppCside+fNpmCside+fNppAside+fNpmAside+fNp0Aside+fNp0Cside+
+               fNetaAside+fNetaCside+fNomegaAside+fNomegaCside;
        }
        fNcurrent++;
+       printf("\t #### Putting %d particles in the stack\n", nTracks);
+       if(fPcToTrack==kAll || fPcToTrack==kOnlyNucleons) printf("\t\t  %d+%d neutrons, %d+%d protons\n", 
+               fNnAside,fNnCside, fNpAside,fNpCside);
+       if(fPcToTrack==kAll || fPcToTrack==kNotNucleons) printf("\t %d+%d pi+, %d+%d pi-, %d+%d pi0, %d+%d eta, %d+%d omega\n",
+               fNppAside,fNppCside,fNpmAside,fNpmCside, 
+               fNp0Aside,fNp0Cside,fNetaAside,fNetaCside, fNomegaAside,fNomegaCside);
        return nTracks;
     }
 
@@ -194,68 +317,127 @@ TParticle* AliGenReaderEMD::NextParticle()
 {
     // Read the next particle
     Float_t p[4]={0.,0.,0.,0.};
+    int pdgCode=0;
     
-    Int_t ipart=0;
-    if(fPcToTrack==0) ipart = kNeutron;
-    else  if(fPcToTrack==1) ipart = kProton;
-    Double_t amass = TDatabasePDG::Instance()->GetParticle(ipart)->Mass();
-
-    if(fPcToTrack==0){
-      if(fIPSide==0){
-        if(fNparticle<fNnCside){
-          p[0] = fPxnCside[fNparticle];
-          p[1] = fPynCside[fNparticle];
-          p[2] = fPznCside[fNparticle];
-       }
-       else{
-          p[0] = fPxnAside[fNparticle-fNnCside];
-          p[1] = fPynAside[fNparticle-fNnCside];
-          p[2] = fPznAside[fNparticle-fNnCside];
-       }
+    if(fPcToTrack==kAll || fPcToTrack==kOnlyNucleons){//***********************************************
+      if(fNparticle<fNnAside){
+        p[0] = fPxnAside[fNparticle];
+        p[1] = fPynAside[fNparticle];
+        p[2] = fPznAside[fNparticle];  
+       pdgCode = fnPDGCode;
+//    printf(" pc%d n sideA: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
       }
-      else if(fIPSide==1){
+      else if(fNparticle>=fNnAside && fNparticle<(fNnAside+fNnCside)){
         p[0] = fPxnCside[fNparticle];
         p[1] = fPynCside[fNparticle];
         p[2] = fPznCside[fNparticle];
+       pdgCode = fnPDGCode;
+//    printf(" pc%d n sideC: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
       }
-      else if(fIPSide==2){
-        p[0] = fPxnAside[fNparticle-fNpCside];
-        p[1] = fPynAside[fNparticle-fNpCside];
-        p[2] = fPznAside[fNparticle-fNpCside];
-      }
-    }
-    else if(fPcToTrack==1){
-      if(fIPSide==0){
-        if(fNparticle<fNnCside){
-          p[0] = fPxpCside[fNparticle];
-          p[1] = fPypCside[fNparticle];
-          p[2] = fPzpCside[fNparticle];
-       }
-       else{
-          p[0] = fPxpAside[fNparticle];
-          p[1] = fPypAside[fNparticle];
-          p[2] = fPzpAside[fNparticle];
-       }
+      else if(fNparticle>=fNnAside+fNnCside && fNparticle<(fNnAside+fNnCside+fNpAside)){
+        p[0] = fPxpAside[fNparticle];
+        p[1] = fPypAside[fNparticle];
+        p[2] = fPzpAside[fNparticle];
+       pdgCode = fpPDGCode;
+//    printf(" pc%d p sideA: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
       }
-      else if(fIPSide==1){
+      else if(fNparticle>=fNnAside+fNnCside+fNpAside && fNparticle<(fNnAside+fNnCside+fNpCside+fNpAside)){
         p[0] = fPxpCside[fNparticle];
         p[1] = fPypCside[fNparticle];
         p[2] = fPzpCside[fNparticle];
+       pdgCode = fpPDGCode;
+//    printf(" pc%d p sideC: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
       }
-      else if(fIPSide==2){
-        p[0] = fPxpAside[fNparticle];
-        p[1] = fPypAside[fNparticle];
-        p[2] = fPzpAside[fNparticle];
+      fOffset = fNnAside+fNnCside+fNpCside+fNpAside;
+    } //**********************************************************************************************
+    if(fPcToTrack==kAll || fPcToTrack==kNotNucleons){
+      if(fNparticle>=fOffset && fNparticle<fOffset+fNppAside){ // *** pi +
+        p[0] = fPxppAside[fNparticle];
+        p[1] = fPyppAside[fNparticle];
+        p[2] = fPzppAside[fNparticle];  
+       pdgCode = fppPDGCode;
+//    printf(" pc%d pi+ sideA: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
+      }
+      if(fNparticle>=fOffset+fNppAside && fNparticle<fOffset+fNppAside+fNppCside){
+        p[0] = fPxppCside[fNparticle];
+        p[1] = fPyppCside[fNparticle];
+        p[2] = fPzppCside[fNparticle];
+       pdgCode = fppPDGCode;
+//    printf(" pc%d pi+ sideC: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
+      }
+      if(fNparticle>=fOffset+fNppAside+fNppCside && fNparticle<fOffset+fNppAside+fNppCside+fNpmAside){ // *** pi -
+        p[0] = fPxpmAside[fNparticle];
+        p[1] = fPypmAside[fNparticle];
+        p[2] = fPzpmAside[fNparticle];
+       pdgCode = fpmPDGCode;
+//    printf(" pc%d pi- sideA: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
+      }
+      if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside && fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside){
+        p[0] = fPxpmCside[fNparticle];
+        p[1] = fPypmCside[fNparticle];
+        p[2] = fPzpmCside[fNparticle];
+       pdgCode = fpmPDGCode;
+//    printf(" pc%d pi- sideC: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
+      }
+      if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside && 
+         fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside){ // *** pi 0
+        p[0] = fPxp0Aside[fNparticle];
+        p[1] = fPyp0Aside[fNparticle];
+        p[2] = fPzp0Aside[fNparticle];
+       pdgCode = fp0PDGCode;
+//    printf(" pc%d pi0 sideA: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
+      }
+      if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside && 
+        fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside){
+        p[0] = fPxp0Cside[fNparticle];
+        p[1] = fPyp0Cside[fNparticle];
+        p[2] = fPzp0Cside[fNparticle];
+       pdgCode = fp0PDGCode;
+//    printf(" pc%d pi0 sideC: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
       }
+      if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside && 
+         fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside){ // *** eta
+        p[0] = fPxetaAside[fNparticle];
+        p[1] = fPyetaAside[fNparticle];
+        p[2] = fPzetaAside[fNparticle];
+       pdgCode = fetaPDGCode;
+//    printf(" pc%d eta sideA: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
+      }
+      if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside && 
+         fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside+fNetaCside){
+        p[0] = fPxetaCside[fNparticle];
+        p[1] = fPyetaCside[fNparticle];
+        p[2] = fPzetaCside[fNparticle];
+       pdgCode = fetaPDGCode;
+//    printf(" pc%d eta sideC: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
+      }
+      if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside+fNetaCside && 
+         fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside+fNetaCside+fNomegaAside){ // *** omega
+        p[0] = fPxomegaAside[fNparticle];
+        p[1] = fPyomegaAside[fNparticle];
+        p[2] = fPzomegaAside[fNparticle];
+       pdgCode = fomegaPDGCode;
+//    printf(" pc%d omega sideA: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
+      }
+      if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside+fNetaCside+fNomegaAside
+      && fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside+fNetaCside+fNomegaAside+fNomegaCside){
+        p[0] = fPxomegaCside[fNparticle];
+        p[1] = fPyomegaCside[fNparticle];
+        p[2] = fPzomegaCside[fNparticle];
+       pdgCode = fomegaPDGCode;
+//    printf(" pc%d omega sideC: PDG code %d,  momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
+      }
+      
     } 
+   
     Float_t ptot = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
+    Double_t amass = TDatabasePDG::Instance()->GetParticle(pdgCode)->Mass();
     p[3] = TMath::Sqrt(ptot*ptot+amass*amass);
-    //printf(" * pc %d: momentum (%f, %f, %f) \n", fNparticle, p[0],p[1],p[2]);
     
     if(p[3]<=amass) 
-      Warning("Generate","Particle %d  E = %f mass = %f \n",ipart,p[3],amass);
-
-    TParticle* particle = new TParticle(ipart, 0, -1, -1, -1, -1, 
+      Warning("Generate","Particle %d  E = %f GeV mass = %f GeV/c^2",pdgCode,p[3],amass);
+    
+    TParticle* particle = new TParticle(pdgCode, 0, -1, -1, -1, -1, 
        p[0], p[1], p[2], p[3], 0., 0., 0., 0.);
     particle->SetBit(kTransportBit);
     fNparticle++;