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 #include "AliGenReaderEMD.h"
  28 #include "AliStack.h"
  29
  30
  31 ClassImp(AliGenReaderEMD)
  32
  33 AliGenReaderEMD::AliGenReaderEMD():
  34     fStartEvent(0),
  35     fNcurrent(0),
  36     fNparticle(0),
  37     fTreeNtuple(0),
  38     fIPSide(0),
  39     fPcToTrack(0),
  40     fNnAside(0),
  41     fEnAside(0),
  42     fNnCside(0),
  43     fEnCside(0),
  44     fNpAside(0),
  45     fEtapAside(0),
  46     fNpCside(0),
  47     fEtapCside(0)
  48 {
  49 // Default constructor
  50 }
  51
  52 AliGenReaderEMD::AliGenReaderEMD(const AliGenReaderEMD &reader):
  53     AliGenReader(reader),
  54     fStartEvent(0),
  55     fNcurrent(0),
  56     fNparticle(0),
  57     fTreeNtuple(0),
  58     fIPSide(0),
  59     fPcToTrack(0),
  60     fNnAside(0),
  61     fEnAside(0),
  62     fNnCside(0),
  63     fEnCside(0),
  64     fNpAside(0),
  65     fEtapAside(0),
  66     fNpCside(0),
  67     fEtapCside(0)
  68 {
  69     // Copy Constructor
  70     reader.Copy(*this);
  71 }
  72   // -----------------------------------------------------------------------------------
  73 AliGenReaderEMD::~AliGenReaderEMD()
  74 {
  75     delete fTreeNtuple;
  76 }
  77
  78 // -----------------------------------------------------------------------------------
  79 AliGenReaderEMD& AliGenReaderEMD::operator=(const  AliGenReaderEMD& rhs)
  80 {
  81 // Assignment operator
  82     rhs.Copy(*this);
  83     return *this;
  84 }
  85
  86 // -----------------------------------------------------------------------------------
  87 void AliGenReaderEMD::Copy(TObject&) const
  88 {
  89     //
  90     // Copy
  91     //
  92     Fatal("Copy","Not implemented!\n");
  93 }
  94
  95 // -----------------------------------------------------------------------------------
  96 void AliGenReaderEMD::Init()
  97 {
  98 //
  99 // Reset the existing file environment and open a new root file
 100
 101     TFile *pFile=0;
 102     if (!pFile) {
 103         pFile = new TFile(fFileName);
 104         pFile->cd();
 105         printf("\n %s file opened\n\n", fFileName);
 106     }
 107     fTreeNtuple = (TTree*)gDirectory->Get("h2032");
 108     fNcurrent = fStartEvent;
 109
 110     TTree *Ntu=fTreeNtuple;
 111     //
 112     // Set branch addresses
 113     // **** neutrons
 114     Ntu->SetBranchAddress("Nleft",&fNnAside);
 115     Ntu->SetBranchAddress("Eleft",&fEnAside);
 116     Ntu->SetBranchAddress("Pxl",  fPxnAside);
 117     Ntu->SetBranchAddress("Pyl",  fPynAside);
 118     Ntu->SetBranchAddress("Pzl",  fPznAside);
 119     Ntu->SetBranchAddress("Nright",&fNnCside);
 120     Ntu->SetBranchAddress("Eright",&fEnCside);
 121     Ntu->SetBranchAddress("Pxr",   fPxnCside);
 122     Ntu->SetBranchAddress("Pyr",   fPynCside);
 123     Ntu->SetBranchAddress("Pzr",   fPznCside);
 124     // **** protons
 125     Ntu->SetBranchAddress("Nleft_p",&fNpAside);
 126     Ntu->SetBranchAddress("Etaleft_p",&fEtapAside);
 127     Ntu->SetBranchAddress("Pxl_p",  fPxpAside);
 128     Ntu->SetBranchAddress("Pyl_p",  fPypAside);
 129     Ntu->SetBranchAddress("Pzl_p",  fPzpAside);
 130     Ntu->SetBranchAddress("Nright_p",&fNpCside);
 131     Ntu->SetBranchAddress("Etaright_p",&fEtapCside);
 132     Ntu->SetBranchAddress("Pxr_p",   fPxpCside);
 133     Ntu->SetBranchAddress("Pyr_p",   fPypCside);
 134     Ntu->SetBranchAddress("Pzr_p",   fPzpCside);
 135 }
 136
 137 // -----------------------------------------------------------------------------------
 138 Int_t AliGenReaderEMD::NextEvent()
 139 {
 140     // Read the next event
 141     Int_t nTracks=0;
 142     fNparticle = 0;
 143
 144     TFile* pFile = fTreeNtuple->GetCurrentFile();
 145     pFile->cd();
 146
 147
 148     Int_t nentries = (Int_t) fTreeNtuple->GetEntries();
 149     if(fNcurrent < nentries) {
 150         fTreeNtuple->GetEvent(fNcurrent);
 151         printf("\n *** Event %d read ***\n",fNcurrent);
 152         //printf("\t A side-> %d neutrons and %d protons emitted\n", fNnAside, fNpAside);
 153         //printf("\t C side-> %d neutrons and %d protons emitted\n\n", fNnCside, fNpCside);
 154         //
 155         // **** IPSIde = 0->both sides, 1->only Cside, 2->only A side
 156         // #### fPcToTrack = 0->neutrons, 1->protons
 157         if(fPcToTrack==0){ // neutrons
 158           if(fIPSide==0){
 159             printf("\t Tracking %d+%d neutrons emitted on C+A side\n", fNnCside,fNnAside);
 160             nTracks = fNnCside+fNnAside;
 161           }
 162           else if(fIPSide==1){
 163             printf("\t Tracking %d neutrons emitted on C side\n", fNnCside);
 164             nTracks    = fNnCside;
 165           }
 166           else if(fIPSide==2){
 167             printf("\t Tracking %d neutrons emitted on A side\n", fNnAside);
 168             nTracks    = fNnAside;
 169           }
 170         }
 171         else if(fPcToTrack==1){ //protons
 172           if(fIPSide==0){
 173             printf("\t Tracking %d+%d protons emitted on C+A side\n", fNpCside,fNpAside);
 174             nTracks    = fNpCside+fNpAside;
 175           }
 176           else if(fIPSide==1){
 177             printf("\t Tracking %d protons emitted on C side\n", fNpCside);
 178             nTracks    = fNpCside;
 179           }
 180           else if(fIPSide==2){
 181             printf("\t Tracking %d protons emitted on A side\n", fNpAside);
 182             nTracks    = fNpAside;
 183           }
 184         }
 185         fNcurrent++;
 186         return nTracks;
 187     }
 188
 189     return 0;
 190 }
 191
 192 // -----------------------------------------------------------------------------------
 193 TParticle* AliGenReaderEMD::NextParticle()
 194 {
 195     // Read the next particle
 196     Float_t p[4]={0.,0.,0.,0.};
 197
 198     Int_t ipart=0;
 199     if(fPcToTrack==0) ipart = kNeutron;
 200     else  if(fPcToTrack==1) ipart = kProton;
 201     Double_t amass = TDatabasePDG::Instance()->GetParticle(ipart)->Mass();
 202
 203     if(fPcToTrack==0){
 204       if(fIPSide==0){
 205         if(fNparticle<fNnCside){
 206           p[0] = fPxnCside[fNparticle];
 207           p[1] = fPynCside[fNparticle];
 208           p[2] = fPznCside[fNparticle];
 209         }
 210         else{
 211           p[0] = fPxnAside[fNparticle-fNnCside];
 212           p[1] = fPynAside[fNparticle-fNnCside];
 213           p[2] = fPznAside[fNparticle-fNnCside];
 214         }
 215       }
 216       else if(fIPSide==1){
 217         p[0] = fPxnCside[fNparticle];
 218         p[1] = fPynCside[fNparticle];
 219         p[2] = fPznCside[fNparticle];
 220       }
 221       else if(fIPSide==2){
 222         p[0] = fPxnAside[fNparticle-fNpCside];
 223         p[1] = fPynAside[fNparticle-fNpCside];
 224         p[2] = fPznAside[fNparticle-fNpCside];
 225       }
 226     }
 227     else if(fPcToTrack==1){
 228       if(fIPSide==0){
 229         if(fNparticle<fNnCside){
 230           p[0] = fPxpCside[fNparticle];
 231           p[1] = fPypCside[fNparticle];
 232           p[2] = fPzpCside[fNparticle];
 233         }
 234         else{
 235           p[0] = fPxpAside[fNparticle];
 236           p[1] = fPypAside[fNparticle];
 237           p[2] = fPzpAside[fNparticle];
 238         }
 239       }
 240       else if(fIPSide==1){
 241         p[0] = fPxpCside[fNparticle];
 242         p[1] = fPypCside[fNparticle];
 243         p[2] = fPzpCside[fNparticle];
 244       }
 245       else if(fIPSide==2){
 246         p[0] = fPxpAside[fNparticle];
 247         p[1] = fPypAside[fNparticle];
 248         p[2] = fPzpAside[fNparticle];
 249       }
 250     }
 251     Float_t ptot = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
 252     p[3] = TMath::Sqrt(ptot*ptot+amass*amass);
 253     //printf(" * pc %d: momentum (%f, %f, %f) \n", fNparticle, p[0],p[1],p[2]);
 254
 255     if(p[3]<=amass)
 256       Warning("Generate","Particle %d  E = %f mass = %f \n",ipart,p[3],amass);
 257
 258     TParticle* particle = new TParticle(ipart, 0, -1, -1, -1, -1,
 259         p[0], p[1], p[2], p[3], 0., 0., 0., 0.);
 260     particle->SetBit(kTransportBit);
 261     fNparticle++;
 262     return particle;
 263 }