[u/mrichter/AliRoot.git] / EVGEN / AliGenReaderEMD.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/


// Class to read events from external (TNtupla) file
// Events -> neutron removal by EM dissociation of Pb nuclei
// Data from RELDIS code (by I. Pshenichov)

#include <TFile.h>
#include <TParticle.h>
#include <TTree.h>
#include <TVirtualMC.h>
#include <TDatabasePDG.h>
#include <TPDGCode.h>
#include "AliGenReaderEMD.h"
#include "AliStack.h"


ClassImp(AliGenReaderEMD)

AliGenReaderEMD::AliGenReaderEMD():
    fStartEvent(0),
    fNcurrent(0),  
    fNparticle(0), 
    fTreeNtuple(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),
    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)
{
// 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),
    fPcToTrack(0),
    fOffset(0),
    fNnAside(0),
    fEnAside(0),
    fnPDGCode(0),
    fNnCside(0),
    fEnCside(0),
    fNpAside(0),
    fEtapAside(0),
    fpPDGCode(0),
    fNpCside(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);
}
  // -----------------------------------------------------------------------------------
AliGenReaderEMD::~AliGenReaderEMD()
{
    delete fTreeNtuple;
}

// -----------------------------------------------------------------------------------
AliGenReaderEMD& AliGenReaderEMD::operator=(const  AliGenReaderEMD& rhs)
{
// Assignment operator
    rhs.Copy(*this);
    return *this;
}

// -----------------------------------------------------------------------------------
void AliGenReaderEMD::Copy(TObject&) const
{
    //
    // Copy 
    //
    Fatal("Copy","Not implemented!\n");
}

// -----------------------------------------------------------------------------------
void AliGenReaderEMD::Init() 
{
//
// Reset the existing file environment and open a new root file
    
    TFile *pFile=0;
    if (!pFile) {
	pFile = new TFile(fFileName);
	pFile->cd();
	printf("\n %s file opened to read RELDIS EMD events\n\n", fFileName);
    }
    fTreeNtuple = (TTree*)gDirectory->Get("h2032");
    fNcurrent = fStartEvent;

    TTree *Ntu=fTreeNtuple;
    //
    // Set branch addresses
    // **** 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);
    Ntu->SetBranchAddress("Nright",&fNnCside);
    Ntu->SetBranchAddress("Eright",&fEnCside);
    Ntu->SetBranchAddress("Pxr",   fPxnCside);
    Ntu->SetBranchAddress("Pyr",   fPynCside);
    Ntu->SetBranchAddress("Pzr",   fPznCside);
    // **** 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);
    Ntu->SetBranchAddress("Nright_p",&fNpCside);
    Ntu->SetBranchAddress("Etaright_p",&fEtapCside);
    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);
}

// -----------------------------------------------------------------------------------
Int_t AliGenReaderEMD::NextEvent() 
{
    // Read the next event  
    Int_t nTracks=0;
    fNparticle = 0; fOffset=0;
    
    TFile* pFile = fTreeNtuple->GetCurrentFile();
    pFile->cd();
    

    Int_t nentries = (Int_t) fTreeNtuple->GetEntries();
    if(fNcurrent < nentries) {
	fTreeNtuple->GetEvent(fNcurrent);
	if(fNcurrent%100 == 0) printf("\n *** Reading event %d ***\n",fNcurrent);
	//
	if(fPcToTrack==kAll || fPcToTrack==kOnlyNucleons){ // nucleons
	   nTracks = fNnCside+fNnAside+fNpCside+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;
    }

    return 0;
}

// -----------------------------------------------------------------------------------
TParticle* AliGenReaderEMD::NextParticle() 
{
    // Read the next particle
    Float_t p[4]={0.,0.,0.,0.};
    int pdgCode=0;
    
    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(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(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(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]);
      }
      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);
    
    if(p[3]<=amass) 
      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++;
    return particle;
}
Commit	Line	Data
f1d1affa	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>
f1d1affa	27	#include "AliGenReaderEMD.h"
2078b0f7	28	#include "AliStack.h"
2078b0f7	29
f1d1affa	30
706938e6	31	ClassImp(AliGenReaderEMD)
f1d1affa	32
1c56e311	33	AliGenReaderEMD::AliGenReaderEMD():
	34	fStartEvent(0),
	35	fNcurrent(0),
	36	fNparticle(0),
	37	fTreeNtuple(0),
1c56e311	38	fPcToTrack(0),
a16b785f	39	fOffset(0),
2078b0f7	40	fNnAside(0),
2078b0f7	41	fEnAside(0),
a16b785f	42	fnPDGCode(0),
2078b0f7	43	fNnCside(0),
	44	fEnCside(0),
	45	fNpAside(0),
	46	fEtapAside(0),
a16b785f	47	fpPDGCode(0),
2078b0f7	48	fNpCside(0),
a16b785f	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)
f1d1affa	75	{
a16b785f	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");
f1d1affa	102	}
f1d1affa	103
a16b785f	104
1c56e311	105	AliGenReaderEMD::AliGenReaderEMD(const AliGenReaderEMD &reader):
	106	AliGenReader(reader),
	107	fStartEvent(0),
	108	fNcurrent(0),
	109	fNparticle(0),
	110	fTreeNtuple(0),
1c56e311	111	fPcToTrack(0),
a16b785f	112	fOffset(0),
2078b0f7	113	fNnAside(0),
2078b0f7	114	fEnAside(0),
a16b785f	115	fnPDGCode(0),
2078b0f7	116	fNnCside(0),
	117	fEnCside(0),
	118	fNpAside(0),
	119	fEtapAside(0),
a16b785f	120	fpPDGCode(0),
2078b0f7	121	fNpCside(0),
a16b785f	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)
1c56e311	148	{
	149	// Copy Constructor
	150	reader.Copy(*this);
	151	}
f1d1affa	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();
a16b785f	185	printf("\n %s file opened to read RELDIS EMD events\n\n", fFileName);
f1d1affa	186	}
	187	fTreeNtuple = (TTree*)gDirectory->Get("h2032");
	188	fNcurrent = fStartEvent;
	189
	190	TTree *Ntu=fTreeNtuple;
	191	//
	192	// Set branch addresses
	193	// **** neutrons
2078b0f7	194	Ntu->SetBranchAddress("Nleft",&fNnAside);
2078b0f7	195	Ntu->SetBranchAddress("Eleft",&fEnAside);
a16b785f	196	Ntu->SetBranchAddress("Ipdg_l_n",&fnPDGCode);
2078b0f7	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);
f1d1affa	205	// **** protons
2078b0f7	206	Ntu->SetBranchAddress("Nleft_p",&fNpAside);
2078b0f7	207	Ntu->SetBranchAddress("Etaleft_p",&fEtapAside);
a16b785f	208	Ntu->SetBranchAddress("Ipdg_l_p",&fpPDGCode);
2078b0f7	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);
a16b785f	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);
f1d1affa	277	}
	278
	279	// -----------------------------------------------------------------------------------
	280	Int_t AliGenReaderEMD::NextEvent()
	281	{
	282	// Read the next event
	283	Int_t nTracks=0;
a16b785f	284	fNparticle = 0; fOffset=0;
f1d1affa	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);
a16b785f	293	if(fNcurrent%100 == 0) printf("\n * Reading event %d *\n",fNcurrent);
f1d1affa	294	//
a16b785f	295	if(fPcToTrack==kAll \|\| fPcToTrack==kOnlyNucleons){ // nucleons
a16b785f	296	nTracks = fNnCside+fNnAside+fNpCside+fNpAside;
f1d1affa	297	}
a16b785f	298	if(fPcToTrack==kAll \|\| fPcToTrack==kNotNucleons){ //pions,eta,omega
	299	nTracks += fNppCside+fNpmCside+fNppAside+fNpmAside+fNp0Aside+fNp0Cside+
	300	fNetaAside+fNetaCside+fNomegaAside+fNomegaCside;
f1d1affa	301	}
2078b0f7	302	fNcurrent++;
a16b785f	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);
f1d1affa	309	return nTracks;
	310	}
	311
	312	return 0;
	313	}
	314
	315	// -----------------------------------------------------------------------------------
	316	TParticle* AliGenReaderEMD::NextParticle()
	317	{
	318	// Read the next particle
2078b0f7	319	Float_t p[4]={0.,0.,0.,0.};
a16b785f	320	int pdgCode=0;
2078b0f7	321
a16b785f	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]);
2078b0f7	329	}
a16b785f	330	else if(fNparticle>=fNnAside && fNparticle<(fNnAside+fNnCside)){
2078b0f7	331	p[0] = fPxnCside[fNparticle];
	332	p[1] = fPynCside[fNparticle];
	333	p[2] = fPznCside[fNparticle];
a16b785f	334	pdgCode = fnPDGCode;
a16b785f	335	// printf(" pc%d n sideC: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
2078b0f7	336	}
a16b785f	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]);
2078b0f7	343	}
a16b785f	344	else if(fNparticle>=fNnAside+fNnCside+fNpAside && fNparticle<(fNnAside+fNnCside+fNpCside+fNpAside)){
2078b0f7	345	p[0] = fPxpCside[fNparticle];
	346	p[1] = fPypCside[fNparticle];
	347	p[2] = fPzpCside[fNparticle];
a16b785f	348	pdgCode = fpPDGCode;
a16b785f	349	// printf(" pc%d p sideC: PDG code %d, momentum (%f, %f, %f) \n", fNparticle, pdgCode, p[0],p[1],p[2]);
2078b0f7	350	}
a16b785f	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]);
2078b0f7	397	}
a16b785f	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
2078b0f7	431	}
a16b785f	432
f1d1affa	433	Float_t ptot = TMath::Sqrt(p[0]p[0]+p[1]p[1]+p[2]*p[2]);
a16b785f	434	Double_t amass = TDatabasePDG::Instance()->GetParticle(pdgCode)->Mass();
f1d1affa	435	p[3] = TMath::Sqrt(ptotptot+amassamass);
f1d1affa	436
f1d1affa	437	if(p[3]<=amass)
a16b785f	438	Warning("Generate","Particle %d E = %f GeV mass = %f GeV/c^2",pdgCode,p[3],amass);
	439
	440	TParticle* particle = new TParticle(pdgCode, 0, -1, -1, -1, -1,
f1d1affa	441	p[0], p[1], p[2], p[3], 0., 0., 0., 0.);
2078b0f7	442	particle->SetBit(kTransportBit);
f1d1affa	443	fNparticle++;
	444	return particle;
	445	}