EVGEN/AliGenExtFile.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 $Log$
  18 Revision 1.15  2001/01/23 13:29:37  morsch
  19 Add method SetParticleCode and enum type Code_t to handle both PDG (new ntuples)
  20 and GEANT3 codes (old ntuples) in input file.
  21
  22 Revision 1.14  2000/12/21 16:24:06  morsch
  23 Coding convention clean-up
  24
  25 Revision 1.13  2000/11/30 07:12:50  alibrary
  26 Introducing new Rndm and QA classes
  27
  28 Revision 1.12  2000/10/27 13:54:45  morsch
  29 Remove explicite reference to input file from constuctor.
  30
  31 Revision 1.11  2000/10/02 21:28:06  fca
  32 Removal of useless dependecies via forward declarations
  33
  34 Revision 1.10  2000/07/11 18:24:55  fca
  35 Coding convention corrections + few minor bug fixes
  36
  37 Revision 1.9  2000/06/14 15:20:09  morsch
  38 Include clean-up (IH)
  39
  40 Revision 1.8  2000/06/09 20:36:44  morsch
  41 All coding rule violations except RS3 corrected
  42
  43 Revision 1.7  2000/02/16 14:56:27  morsch
  44 Convert geant particle code into pdg code before putting particle on the stack.
  45
  46 Revision 1.6  1999/11/09 07:38:48  fca
  47 Changes for compatibility with version 2.23 of ROOT
  48
  49 Revision 1.5  1999/09/29 09:24:12  fca
  50 Introduction of the Copyright and cvs Log
  51
  52 */
  53
  54
  55 // Event generator that can read the old ALICE event format based on CW-ntuples
  56 // http://consult.cern.ch/alice/Internal_Notes/1995/32/abstract
  57 // .cwn file have to be converted to .root using h2root
  58 // Use SetFileName(file) to read from "file"
  59 // Author: andreas.morsch@cern.ch
  60
  61 #include <iostream.h>
  62
  63 #include "AliGenExtFile.h"
  64 #include "AliMC.h"
  65 #include "AliRun.h"
  66
  67 #include <TDirectory.h>
  68 #include <TDatabasePDG.h>
  69 #include <TFile.h>
  70 #include "TTree.h"
  71 #include <stdlib.h>
  72
  73  ClassImp(AliGenExtFile)
  74      AliGenExtFile::AliGenExtFile()
  75          :AliGenerator(-1)
  76 {
  77 //  Constructor
  78     fName="ExtFile";
  79     fTitle="Primaries from ext. File";
  80     fFileName="";
  81     fTreeNtuple=0;
  82     fNcurrent=0;
  83     fCode = kGEANT3;
  84 //
  85 //  Read all particles
  86     fNpart=-1;
  87 }
  88
  89 AliGenExtFile::AliGenExtFile(Int_t npart)
  90     :AliGenerator(npart)
  91 {
  92 //  Constructor
  93     fName="ExtFile";
  94     fTitle="Primaries from ext. File";
  95     fFileName="";
  96     fCode = kGEANT3;
  97     fTreeNtuple=0;
  98     fNcurrent=0;
  99 }
 100
 101 AliGenExtFile::AliGenExtFile(const AliGenExtFile & ExtFile)
 102 {
 103 // copy constructor
 104 }
 105 //____________________________________________________________
 106 AliGenExtFile::~AliGenExtFile()
 107 {
 108 // Destructor
 109     delete fTreeNtuple;
 110 }
 111
 112 //____________________________________________________________
 113 void AliGenExtFile::NtupleInit()
 114 {
 115 //
 116 // reset the existing file environment and open a new root file if
 117 // the pointer to the Fluka tree is null
 118
 119     TFile *pFile=0;
 120     if (fTreeNtuple==0) {
 121         if (!pFile) {
 122             pFile = new TFile(fFileName);
 123             pFile->cd();
 124             cout<<"I have opened "<<fFileName<<" file "<<endl;
 125         }
 126 // get the tree address in the Fluka boundary source file
 127         fTreeNtuple = (TTree*)gDirectory->Get("h888");
 128     } else {
 129         pFile = fTreeNtuple->GetCurrentFile();
 130         pFile->cd();
 131     }
 132
 133     TTree *h2=fTreeNtuple;
 134 //Set branch addresses
 135     h2->SetBranchAddress("Nihead",&fNihead);
 136     h2->SetBranchAddress("Ihead",fIhead);
 137     h2->SetBranchAddress("Nrhead",&fNrhead);
 138     h2->SetBranchAddress("Rhead",fRhead);
 139     h2->SetBranchAddress("Idpart",&fIdpart);
 140     h2->SetBranchAddress("Theta",&fTheta);
 141     h2->SetBranchAddress("Phi",&fPhi);
 142     h2->SetBranchAddress("P",&fP);
 143     h2->SetBranchAddress("E",&fE);
 144 }
 145
 146
 147 //____________________________________________________________
 148 void AliGenExtFile::Generate()
 149 {
 150 // Generate particles
 151
 152   Float_t polar[3]= {0,0,0};
 153   //
 154   Float_t origin[3]={0,0,0};
 155   Float_t p[3];
 156   Float_t random[6];
 157   Float_t prwn;
 158   Int_t i, j, nt, nTracks=0;
 159   //
 160   NtupleInit();
 161   TTree *h2=fTreeNtuple;
 162   Int_t nentries = (Int_t) h2->GetEntries();
 163   // loop over number of particles
 164   Int_t nb = (Int_t)h2->GetEvent(fNcurrent);
 165   Int_t i5=fIhead[4];
 166   Int_t i6=fIhead[5];
 167
 168   for (j=0;j<3;j++) origin[j]=fOrigin[j];
 169   if(fVertexSmear==kPerTrack) {
 170     Rndm(random,6);
 171     for (j=0;j<3;j++) {
 172         origin[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
 173             TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 174     }
 175   }
 176
 177   if (fNcurrent >= nentries) {
 178       printf("\n No more entries !!! !\n");
 179       return;
 180   }
 181
 182
 183   if (i5==0) {
 184       printf("\n This should never happen !\n");
 185   } else {
 186       printf("\n Next event contains %d tracks! \n", i6);
 187       nTracks=i6;
 188   }
 189   for (i=0; i<nTracks; i++) {
 190       if (fCode == kGEANT3) fIdpart=gMC->PDGFromId(fIdpart);
 191       Double_t amass = TDatabasePDG::Instance()->GetParticle(fIdpart)->Mass();
 192       if(fE<=amass) {
 193         Warning("Generate","Particle %d no %d E = %f mass = %f %f %f \n",
 194                 fIdpart,i,fE,amass, fPhi, fTheta);
 195         prwn=0;
 196       } else {
 197         prwn=sqrt((fE+amass)*(fE-amass));
 198       }
 199
 200       fTheta *= TMath::Pi()/180.;
 201       fPhi    = (fPhi-180)*TMath::Pi()/180.;
 202       if(fTheta<fThetaMin || fTheta>fThetaMax ||
 203          fPhi<fPhiMin || fPhi>fPhiMax         ||
 204          prwn<fPMin || prwn>fPMax)
 205       {
 206           ;
 207       } else {
 208           p[0]=prwn*TMath::Sin(fTheta)*TMath::Cos(fPhi);
 209           p[1]=prwn*TMath::Sin(fTheta)*TMath::Sin(fPhi);
 210           p[2]=prwn*TMath::Cos(fTheta);
 211
 212           if(fVertexSmear==kPerTrack) {
 213               Rndm(random,6);
 214               for (j=0;j<3;j++) {
 215                   origin[j]=fOrigin[j]
 216                       +fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
 217                       TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 218               }
 219           }
 220           SetTrack(fTrackIt,-1,fIdpart,p,origin,polar,0,kPPrimary,nt);
 221       }
 222       fNcurrent++;
 223       nb = (Int_t)h2->GetEvent(fNcurrent);
 224   }
 225
 226     TFile *pFile=0;
 227 // Get AliRun object or create it
 228     if (!gAlice) {
 229         gAlice = (AliRun*)pFile->Get("gAlice");
 230         if (gAlice) printf("AliRun object found on file\n");
 231         if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
 232     }
 233     TTree *fAli=gAlice->TreeK();
 234     if (fAli) pFile =fAli->GetCurrentFile();
 235     pFile->cd();
 236 }
 237
 238
 239 AliGenExtFile& AliGenExtFile::operator=(const  AliGenExtFile& rhs)
 240 {
 241 // Assignment operator
 242     return *this;
 243 }
 244
 245
 246
 247
 248
 249
 250
 251