EVGEN/AliGenReaderEMD.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 // Class to read events from external (TNtupla) file
  18 // Events -> neutron removal by EM dissociation of Pb nuclei
  19 // Data from RELDIS code (by I. Pshenichov)
  20
  21 #include <TFile.h>
  22 #include <TParticle.h>
  23 #include <TTree.h>
  24 #include <TVirtualMC.h>
  25 #include <TDatabasePDG.h>
  26 #include <TPDGCode.h>
  27
  28 #include "AliGenReaderEMD.h"
  29
  30 ClassImp(AliGenReaderEMD)
  31
  32 AliGenReaderEMD::AliGenReaderEMD():
  33     fStartEvent(0),
  34     fNcurrent(0),
  35     fNparticle(0),
  36     fTreeNtuple(0),
  37     fIPSide(0),
  38     fPcToTrack(0),
  39     fNnLeft(0),
  40     fEnLeft(0),
  41     fNnRight(0),
  42     fEnRight(0),
  43     fNpLeft(0),
  44     fEtapLeft(0),
  45     fNpRight(0),
  46     fEtapRight(0)
  47 {
  48 // Default constructor
  49 }
  50
  51 AliGenReaderEMD::AliGenReaderEMD(const AliGenReaderEMD &reader):
  52     AliGenReader(reader),
  53     fStartEvent(0),
  54     fNcurrent(0),
  55     fNparticle(0),
  56     fTreeNtuple(0),
  57     fIPSide(0),
  58     fPcToTrack(0),
  59     fNnLeft(0),
  60     fEnLeft(0),
  61     fNnRight(0),
  62     fEnRight(0),
  63     fNpLeft(0),
  64     fEtapLeft(0),
  65     fNpRight(0),
  66     fEtapRight(0)
  67 {
  68     // Copy Constructor
  69     reader.Copy(*this);
  70 }
  71   // -----------------------------------------------------------------------------------
  72 AliGenReaderEMD::~AliGenReaderEMD()
  73 {
  74     delete fTreeNtuple;
  75 }
  76
  77 // -----------------------------------------------------------------------------------
  78 AliGenReaderEMD& AliGenReaderEMD::operator=(const  AliGenReaderEMD& rhs)
  79 {
  80 // Assignment operator
  81     rhs.Copy(*this);
  82     return *this;
  83 }
  84
  85 // -----------------------------------------------------------------------------------
  86 void AliGenReaderEMD::Copy(TObject&) const
  87 {
  88     //
  89     // Copy
  90     //
  91     Fatal("Copy","Not implemented!\n");
  92 }
  93
  94 // -----------------------------------------------------------------------------------
  95 void AliGenReaderEMD::Init()
  96 {
  97 //
  98 // Reset the existing file environment and open a new root file
  99
 100     TFile *pFile=0;
 101     if (!pFile) {
 102         pFile = new TFile(fFileName);
 103         pFile->cd();
 104         printf("\n %s file opened\n\n", fFileName);
 105     }
 106     fTreeNtuple = (TTree*)gDirectory->Get("h2032");
 107     fNcurrent = fStartEvent;
 108
 109     TTree *Ntu=fTreeNtuple;
 110     //
 111     // Set branch addresses
 112     // **** neutrons
 113     Ntu->SetBranchAddress("Nleft",&fNnLeft);
 114     Ntu->SetBranchAddress("Eleft",&fEnLeft);
 115     Ntu->SetBranchAddress("Pxl",  fPxnLeft);
 116     Ntu->SetBranchAddress("Pyl",  fPynLeft);
 117     Ntu->SetBranchAddress("Pzl",  fPznLeft);
 118     Ntu->SetBranchAddress("Nright",&fNnRight);
 119     Ntu->SetBranchAddress("Eright",&fEnRight);
 120     Ntu->SetBranchAddress("Pxr",   fPxnRight);
 121     Ntu->SetBranchAddress("Pyr",   fPynRight);
 122     Ntu->SetBranchAddress("Pzr",   fPznRight);
 123     // **** protons
 124     Ntu->SetBranchAddress("Nleft_p",&fNpLeft);
 125     Ntu->SetBranchAddress("Etaleft_p",&fEtapLeft);
 126     Ntu->SetBranchAddress("Pxl_p",  fPxpLeft);
 127     Ntu->SetBranchAddress("Pyl_p",  fPypLeft);
 128     Ntu->SetBranchAddress("Pzl_p",  fPzpLeft);
 129     Ntu->SetBranchAddress("Nright_p",&fNpRight);
 130     Ntu->SetBranchAddress("Etaright_p",&fEtapRight);
 131     Ntu->SetBranchAddress("Pxr_p",   fPxpRight);
 132     Ntu->SetBranchAddress("Pyr_p",   fPypRight);
 133     Ntu->SetBranchAddress("Pzr_p",   fPzpRight);
 134 }
 135
 136 // -----------------------------------------------------------------------------------
 137 Int_t AliGenReaderEMD::NextEvent()
 138 {
 139     // Read the next event
 140     Int_t nTracks=0;
 141
 142     TFile* pFile = fTreeNtuple->GetCurrentFile();
 143     pFile->cd();
 144
 145
 146     Int_t nentries = (Int_t) fTreeNtuple->GetEntries();
 147     if(fNcurrent < nentries) {
 148         fTreeNtuple->GetEvent(fNcurrent);
 149         printf("\n *** Event %d read ***\n\n",fNcurrent);
 150         fNcurrent++;
 151         //printf("\n \t \t LEFT-> %d neutrons and %d protons emitted\n", fNnLeft, fNpLeft);
 152         //printf("\t \t RIGHT-> %d neutrons and %d protons emitted\n\n", fNnRight, fNpRight);
 153         //
 154         // **** IPSIde          =0->RIGHT side, =1->LEFT side
 155         // #### fPcToTrack      =0->neutrons, =1->protons
 156         if(fIPSide==0){
 157           if(fPcToTrack==0){
 158             printf("\n \t \t Tracking %d neutrons emitted on RIGHT side\n\n", fNnRight);
 159             nTracks    = fNnRight;
 160           }
 161           else if(fPcToTrack==1){
 162             printf("\n \t \t Tracking %d protons emitted on RIGHT side\n\n", fNpRight);
 163             nTracks    = fNpRight;
 164           }
 165         }
 166         else if(fIPSide==1){
 167           if(fPcToTrack==0){
 168             printf("\n \t \t Tracking %d neutrons emitted on LEFT side\n", fNnLeft);
 169             nTracks    = fNnLeft;
 170           }
 171           else if(fPcToTrack==1){
 172             printf("\n \t \t Tracking %d protons emitted on LEFT side\n", fNpLeft);
 173             nTracks    = fNpLeft;
 174           }
 175         }
 176         fNparticle = 0;
 177         return nTracks;
 178     }
 179
 180     return 0;
 181 }
 182
 183 // -----------------------------------------------------------------------------------
 184 TParticle* AliGenReaderEMD::NextParticle()
 185 {
 186     // Read the next particle
 187     Float_t p[4];
 188     Int_t ipart = kNeutron;
 189     Double_t amass = TDatabasePDG::Instance()->GetParticle(kNeutron)->Mass();
 190     p[0] = fPxnRight[fNparticle];
 191     p[1] = fPynRight[fNparticle];
 192     p[2] = fPznRight[fNparticle];
 193     Float_t ptot = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
 194     p[3] = TMath::Sqrt(ptot*ptot+amass*amass);
 195     //printf("\t Nucleon momentum: px = %f, py = %f, pz = %f\n", p[0],p[1],p[2]);
 196
 197     if(p[3]<=amass)
 198       Warning("Generate","Particle %d  E = %f mass = %f \n",ipart,p[3],amass);
 199
 200     TParticle* particle = new TParticle(ipart, 0, -1, -1, -1, -1,
 201         p[0], p[1], p[2], p[3], 0., 0., 0., 0.);
 202     //fNcurrent++;
 203     fNparticle++;
 204     return particle;
 205 }