Updates.
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 22 Feb 2011 17:18:12 +0000 (17:18 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 22 Feb 2011 17:18:12 +0000 (17:18 +0000)
Chiara Oppedisano

EVGEN/AliGenReaderEMD.cxx
EVGEN/AliGenReaderEMD.h

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++;
index 68152dc..0b6b5c3 100644 (file)
@@ -4,7 +4,7 @@
  * See cxx source for full Copyright notice                               */
 
 // Class to read events from external (TNtupla) file
-// Events -> neutron removal by EM dissociation of Pb nuclei
+// Events -> EM dissociation of Pb nuclei
 // Data from RELDIS code (by I. Pshenichov)
 //
 #include "AliGenReader.h"
@@ -14,6 +14,8 @@
 class AliGenReaderEMD : public AliGenReader
 {
  public:
+    enum TrackedPc{kAll=0, kOnlyNucleons=1, kNotNucleons=2};
+    
     AliGenReaderEMD();
     
     AliGenReaderEMD(const AliGenReaderEMD &reader);
@@ -21,13 +23,14 @@ class AliGenReaderEMD : public AliGenReader
     AliGenReaderEMD & operator=(const AliGenReaderEMD & rhs);
     // Initialise 
     virtual void Init();
-    // Read
+    // Reader
     virtual Int_t NextEvent();
     virtual TParticle*  NextParticle();
     virtual void RewindEvent(){;}
-    // Setter
-    void SetIPSide(Int_t iside) {fIPSide = iside;}
-    void SetPcToTrack(Int_t ipart) {fPcToTrack = ipart;}
+    // Setters
+    void TrackNotNucleons()  {fPcToTrack = kNotNucleons;}
+    void TrackOnlyNucleons() {fPcToTrack = kOnlyNucleons;}
+    void TrackAllParticles() {fPcToTrack = kAll;}
     void SetStartEvent(Int_t nev) {fStartEvent = nev;}
     
  protected:
@@ -36,39 +39,111 @@ class AliGenReaderEMD : public AliGenReader
     Int_t           fNparticle;        // number of particles
     TTree           *fTreeNtuple;      // pointer to the TTree
     //
-    Int_t          fIPSide;            // ZDC arm relative to IP
     Int_t          fPcToTrack;         // flag for particles to be tracked
+    Int_t           fOffset;           // Needed to correctly read next particle
     //
     // --- Declaration of leaves types
     // **** neutrons
     Int_t           fNnAside;          // No. of neutrons emitted on left side            
-    Float_t         fEnAside;          // Forward energy Aside side
-    Float_t        fPxnAside[70];      // momentum x component - left side
-    Float_t        fPynAside[70];      // momentum y component - left side     
-    Float_t        fPznAside[70];      // momentum z component - left side   
+    Float_t         fEnAside;          // Forward energy A side
+    Int_t           fnPDGCode;         // PDG code
+    Float_t        fPxnAside[70];      // momentum x component A side
+    Float_t        fPynAside[70];      // momentum y component A side    
+    Float_t        fPznAside[70];      // momentum z component A side  
     //  
     Int_t           fNnCside;          // No. of neutrons emitted on right side            
-    Float_t         fEnCside;          // Forward energy Cside side
-    Float_t        fPxnCside[70];      // momentum x component - right side
-    Float_t        fPynCside[70];      // momentum y component - right side     
-    Float_t        fPznCside[70];      // momentum z component - right side     
+    Float_t         fEnCside;          // Forward energy C side
+    Float_t        fPxnCside[70];      // momentum x component C side
+    Float_t        fPynCside[70];      // momentum y component C side    
+    Float_t        fPznCside[70];      // momentum z component C side    
     //
     // **** protons
     Int_t           fNpAside;          // No. of protons emitted on left side            
-    Float_t         fEtapAside;                // Forward energy Aside side
-    Float_t        fPxpAside[70];      // momentum x component - left side
-    Float_t        fPypAside[70];      // momentum y component - left side     
-    Float_t        fPzpAside[70];      // momentum z component - left side   
+    Float_t         fEtapAside;        // Forward energy A side
+    Int_t           fpPDGCode;         // PDG code
+    Float_t        fPxpAside[50];      // momentum x component A side
+    Float_t        fPypAside[50];      // momentum y component A side    
+    Float_t        fPzpAside[50];      // momentum z component A side  
     //  
     Int_t           fNpCside;          // No. of protons emitted on right side            
-    Float_t         fEtapCside;        // Forward energy Cside side
-    Float_t        fPxpCside[70];      // momentum x component - right side
-    Float_t        fPypCside[70];      // momentum y component - right side     
-    Float_t        fPzpCside[70];      // momentum z component - right side
-    
+    Float_t         fEtapCside;        // Forward energy C side
+    Float_t        fPxpCside[50];      // momentum x component C side
+    Float_t        fPypCside[50];      // momentum y component C side    
+    Float_t        fPzpCside[50];      // momentum z component C side
+    //
+    // **** pi +
+    Int_t           fNppAside;         // No. of pi+ emitted pi+ on A side            
+    Float_t         fEtappAside;       // Forward energy pi+ A side
+    Int_t           fppPDGCode;        // PDG code
+    Float_t        fPxppAside[30];     // momentum x component pi+ A side
+    Float_t        fPyppAside[30];     // momentum y component pi+ A side        
+    Float_t        fPzppAside[30];     // momentum z component pi+ A side      
+    //  
+    Int_t           fNppCside;         // No. of pi+ emitted on C side            
+    Float_t         fEtappCside;       // Forward energy pi+ C side
+    Float_t        fPxppCside[30];     // momentum x component pi+ C side
+    Float_t        fPyppCside[30];     // momentum y component pi+ C side        
+    Float_t        fPzppCside[30];     // momentum z component pi+ C side
+    //
+    // **** pi -
+    Int_t           fNpmAside;         // No. of pi- emitted on A side            
+    Float_t         fEtapmAside;       // Forward energy pi- A side
+    Int_t           fpmPDGCode;        // PDG code
+    Float_t        fPxpmAside[30];     // momentum x component pi- A side
+    Float_t        fPypmAside[30];     // momentum y component pi- A side        
+    Float_t        fPzpmAside[30];     // momentum z component pi- A side      
+    //  
+    Int_t           fNpmCside;         // No. of pi- emitted on C side            
+    Float_t         fEtapmCside;       // Forward energy pi- C side
+    Float_t        fPxpmCside[30];     // momentum x component pi- C side
+    Float_t        fPypmCside[30];     // momentum y component pi- C side        
+    Float_t        fPzpmCside[30];     // momentum z component pi- C side
+    //
+    // **** pi0
+    Int_t           fNp0Aside;         // No. of pi0 emitted on A side            
+    Float_t         fEtap0Aside;       // Forward energy pi0 A side
+    Int_t           fp0PDGCode;        // PDG code
+    Float_t        fPxp0Aside[30];     // momentum x component pi0 A side
+    Float_t        fPyp0Aside[30];     // momentum y component pi0 A side        
+    Float_t        fPzp0Aside[30];     // momentum z component pi0 A side      
+    //  
+    Int_t           fNp0Cside;         // No. of pi0 emitted on C side            
+    Float_t         fEtap0Cside;       // Forward energy pi0 C side
+    Float_t        fPxp0Cside[30];     // momentum x component pi0 C side
+    Float_t        fPyp0Cside[30];     // momentum y component pi0 C side        
+    Float_t        fPzp0Cside[30];     // momentum z component pi0 C side
+    //
+    // **** eta
+    Int_t           fNetaAside;                // No. of eta emitted on A side            
+    Float_t         fEtaetaAside;      // Forward energy eta A side
+    Int_t           fetaPDGCode;       // PDG code
+    Float_t        fPxetaAside[15];    // momentum x component eta A side
+    Float_t        fPyetaAside[15];    // momentum y component eta A side       
+    Float_t        fPzetaAside[15];    // momentum z component eta A side     
+    //  
+    Int_t           fNetaCside;                // No. of eta emitted on C side            
+    Float_t         fEtaetaCside;      // Forward energy eta C side
+    Float_t        fPxetaCside[15];    // momentum x component eta C side
+    Float_t        fPyetaCside[15];    // momentum y component eta C side       
+    Float_t        fPzetaCside[15];    // momentum z component eta C side
+    //
+    // **** omega
+    Int_t           fNomegaAside;      // No. of omega emitted on A side            
+    Float_t         fEtaomegaAside;    // Forward energy omega A side
+    Int_t           fomegaPDGCode;     // PDG code
+    Float_t        fPxomegaAside[15];  // momentum x component omega A side
+    Float_t        fPyomegaAside[15];  // momentum y component omega A side   
+    Float_t        fPzomegaAside[15];  // momentum z component omega A side     
+    //  
+    Int_t           fNomegaCside;      // No. of omega emitted on C side            
+    Float_t         fEtaomegaCside;    // Forward energy omega C side
+    Float_t        fPxomegaCside[15];  // momentum x component omega C side
+    Float_t        fPyomegaCside[15];  // momentum y component omega C side     
+    Float_t        fPzomegaCside[15];  // momentum z component omega C side
+   
  private:
     void Copy(TObject&) const;
-    ClassDef(AliGenReaderEMD,1) // Class to read EMD data
+    ClassDef(AliGenReaderEMD, 2) // Class to read EMD data
 };
 #endif