#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)
{
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)
{
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.;
}
}
}
}
reader.Copy(*this);
}
AliGenReaderEMD::~AliGenReaderEMD()
{
delete fTreeNtuple;
}
AliGenReaderEMD& AliGenReaderEMD::operator=(const AliGenReaderEMD& rhs)
{
rhs.Copy(*this);
return *this;
}
void AliGenReaderEMD::Copy(TObject&) const
{
Fatal("Copy","Not implemented!\n");
}
void AliGenReaderEMD::Init()
{
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;
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);
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);
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);
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);
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);
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);
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()
{
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){
nTracks = fNnCside+fNnAside+fNpCside+fNpAside;
}
if(fPcToTrack==kAll || fPcToTrack==kNotNucleons){
nTracks += fNppCside+fNpmCside+fNppAside+fNpmAside+fNp0Aside+fNp0Cside+
fNetaAside+fNetaCside+fNomegaAside+fNomegaCside;
}
fNcurrent++;
printf("\t #### Putting %d particles in the stack\n", nTracks);
return nTracks;
}
return 0;
}
TParticle* AliGenReaderEMD::NextParticle()
{
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;
}
else if(fNparticle>=fNnAside && fNparticle<(fNnAside+fNnCside)){
p[0] = fPxnCside[fNparticle];
p[1] = fPynCside[fNparticle];
p[2] = fPznCside[fNparticle];
pdgCode = fnPDGCode;
}
else if(fNparticle>=fNnAside+fNnCside && fNparticle<(fNnAside+fNnCside+fNpAside)){
p[0] = fPxpAside[fNparticle];
p[1] = fPypAside[fNparticle];
p[2] = fPzpAside[fNparticle];
pdgCode = fpPDGCode;
}
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;
}
fOffset = fNnAside+fNnCside+fNpCside+fNpAside;
}
if(fPcToTrack==kAll || fPcToTrack==kNotNucleons){
if(fNparticle>=fOffset && fNparticle<fOffset+fNppAside){
p[0] = fPxppAside[fNparticle];
p[1] = fPyppAside[fNparticle];
p[2] = fPzppAside[fNparticle];
pdgCode = fppPDGCode;
}
if(fNparticle>=fOffset+fNppAside && fNparticle<fOffset+fNppAside+fNppCside){
p[0] = fPxppCside[fNparticle];
p[1] = fPyppCside[fNparticle];
p[2] = fPzppCside[fNparticle];
pdgCode = fppPDGCode;
}
if(fNparticle>=fOffset+fNppAside+fNppCside && fNparticle<fOffset+fNppAside+fNppCside+fNpmAside){
p[0] = fPxpmAside[fNparticle];
p[1] = fPypmAside[fNparticle];
p[2] = fPzpmAside[fNparticle];
pdgCode = fpmPDGCode;
}
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;
}
if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside &&
fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside){
p[0] = fPxp0Aside[fNparticle];
p[1] = fPyp0Aside[fNparticle];
p[2] = fPzp0Aside[fNparticle];
pdgCode = fp0PDGCode;
}
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;
}
if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside &&
fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside){
p[0] = fPxetaAside[fNparticle];
p[1] = fPyetaAside[fNparticle];
p[2] = fPzetaAside[fNparticle];
pdgCode = fetaPDGCode;
}
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;
}
if(fNparticle>=fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside+fNetaCside &&
fNparticle<fOffset+fNppAside+fNppCside+fNpmAside+fNpmCside+fNp0Aside+fNp0Cside+fNetaAside+fNetaCside+fNomegaAside){
p[0] = fPxomegaAside[fNparticle];
p[1] = fPyomegaAside[fNparticle];
p[2] = fPzomegaAside[fNparticle];
pdgCode = fomegaPDGCode;
}
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;
}
}
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 ",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.);
if((p[0]*p[0]+p[1]*p[1]+p[2]*p[2])>1e-5) particle->SetBit(kTransportBit);
fNparticle++;
return particle;
}