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
  33   // -----------------------------------------------------------------------------------
  34 AliGenReaderEMD::AliGenReaderEMD()
  35 {
  36 // Default constructor
  37     fStartEvent  = 0;
  38     fTreeNtuple  = 0;
  39     fIPSide      = 0;
  40     fPcToTrack = 0;
  41 }
  42
  43   // -----------------------------------------------------------------------------------
  44 AliGenReaderEMD::~AliGenReaderEMD()
  45 {
  46     delete fTreeNtuple;
  47 }
  48
  49 // -----------------------------------------------------------------------------------
  50 AliGenReaderEMD& AliGenReaderEMD::operator=(const  AliGenReaderEMD& rhs)
  51 {
  52 // Assignment operator
  53     rhs.Copy(*this);
  54     return *this;
  55 }
  56
  57 // -----------------------------------------------------------------------------------
  58 void AliGenReaderEMD::Copy(TObject&) const
  59 {
  60     //
  61     // Copy
  62     //
  63     Fatal("Copy","Not implemented!\n");
  64 }
  65
  66 // -----------------------------------------------------------------------------------
  67 void AliGenReaderEMD::Init()
  68 {
  69 //
  70 // Reset the existing file environment and open a new root file
  71
  72     TFile *pFile=0;
  73     if (!pFile) {
  74         pFile = new TFile(fFileName);
  75         pFile->cd();
  76         printf("\n %s file opened\n\n", fFileName);
  77     }
  78     fTreeNtuple = (TTree*)gDirectory->Get("h2032");
  79     fNcurrent = fStartEvent;
  80
  81     TTree *Ntu=fTreeNtuple;
  82     //
  83     // Set branch addresses
  84     // **** neutrons
  85     Ntu->SetBranchAddress("Nleft",&fNnLeft);
  86     Ntu->SetBranchAddress("Eleft",&fEnLeft);
  87     Ntu->SetBranchAddress("Pxl",  fPxnLeft);
  88     Ntu->SetBranchAddress("Pyl",  fPynLeft);
  89     Ntu->SetBranchAddress("Pzl",  fPznLeft);
  90     Ntu->SetBranchAddress("Nright",&fNnRight);
  91     Ntu->SetBranchAddress("Eright",&fEnRight);
  92     Ntu->SetBranchAddress("Pxr",   fPxnRight);
  93     Ntu->SetBranchAddress("Pyr",   fPynRight);
  94     Ntu->SetBranchAddress("Pzr",   fPznRight);
  95     // **** protons
  96     Ntu->SetBranchAddress("Nleft_p",&fNpLeft);
  97     Ntu->SetBranchAddress("Etaleft_p",&fEtapLeft);
  98     Ntu->SetBranchAddress("Pxl_p",  fPxpLeft);
  99     Ntu->SetBranchAddress("Pyl_p",  fPypLeft);
 100     Ntu->SetBranchAddress("Pzl_p",  fPzpLeft);
 101     Ntu->SetBranchAddress("Nright_p",&fNpRight);
 102     Ntu->SetBranchAddress("Etaright_p",&fEtapRight);
 103     Ntu->SetBranchAddress("Pxr_p",   fPxpRight);
 104     Ntu->SetBranchAddress("Pyr_p",   fPypRight);
 105     Ntu->SetBranchAddress("Pzr_p",   fPzpRight);
 106 }
 107
 108 // -----------------------------------------------------------------------------------
 109 Int_t AliGenReaderEMD::NextEvent()
 110 {
 111     // Read the next event
 112     Int_t nTracks=0;
 113
 114     TFile* pFile = fTreeNtuple->GetCurrentFile();
 115     pFile->cd();
 116
 117
 118     Int_t nentries = (Int_t) fTreeNtuple->GetEntries();
 119     if(fNcurrent < nentries) {
 120         fTreeNtuple->GetEvent(fNcurrent);
 121         printf("\n *** Event %d read ***\n\n",fNcurrent);
 122         fNcurrent++;
 123         //printf("\n \t \t LEFT-> %d neutrons and %d protons emitted\n", fNnLeft, fNpLeft);
 124         //printf("\t \t RIGHT-> %d neutrons and %d protons emitted\n\n", fNnRight, fNpRight);
 125         //
 126         // **** IPSIde          =0->RIGHT side, =1->LEFT side
 127         // #### fPcToTrack      =0->neutrons, =1->protons
 128         if(fIPSide==0){
 129           if(fPcToTrack==0){
 130             printf("\n \t \t Tracking %d neutrons emitted on RIGHT side\n\n", fNnRight);
 131             nTracks    = fNnRight;
 132           }
 133           else if(fPcToTrack==1){
 134             printf("\n \t \t Tracking %d protons emitted on RIGHT side\n\n", fNpRight);
 135             nTracks    = fNpRight;
 136           }
 137         }
 138         else if(fIPSide==1){
 139           if(fPcToTrack==0){
 140             printf("\n \t \t Tracking %d neutrons emitted on LEFT side\n", fNnLeft);
 141             nTracks    = fNnLeft;
 142           }
 143           else if(fPcToTrack==1){
 144             printf("\n \t \t Tracking %d protons emitted on LEFT side\n", fNpLeft);
 145             nTracks    = fNpLeft;
 146           }
 147         }
 148         fNparticle = 0;
 149         return nTracks;
 150     }
 151
 152     return 0;
 153 }
 154
 155 // -----------------------------------------------------------------------------------
 156 TParticle* AliGenReaderEMD::NextParticle()
 157 {
 158     // Read the next particle
 159     Float_t p[4];
 160     Int_t ipart = kNeutron;
 161     Double_t amass = TDatabasePDG::Instance()->GetParticle(kNeutron)->Mass();
 162     p[0] = fPxnRight[fNparticle];
 163     p[1] = fPynRight[fNparticle];
 164     p[2] = fPznRight[fNparticle];
 165     Float_t ptot = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
 166     p[3] = TMath::Sqrt(ptot*ptot+amass*amass);
 167     //printf("\t Nucleon momentum: px = %f, py = %f, pz = %f\n", p[0],p[1],p[2]);
 168
 169     if(p[3]<=amass)
 170       Warning("Generate","Particle %d  E = %f mass = %f \n",ipart,p[3],amass);
 171
 172     TParticle* particle = new TParticle(ipart, 0, -1, -1, -1, -1,
 173         p[0], p[1], p[2], p[3], 0., 0., 0., 0.);
 174     //fNcurrent++;
 175     fNparticle++;
 176     return particle;
 177 }