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 /* $Id$ */
  17
  18 // Event generator that using an instance of type AliGenReader
  19 // reads particles from a file and applies cuts.
  20 // Example: In your Config.C you can include the following lines
  21 //  AliGenExtFile *gener = new AliGenExtFile(-1);
  22 //  gener->SetMomentumRange(0,999);
  23 //  gener->SetPhiRange(-180.,180.);
  24 //  gener->SetThetaRange(0,180);
  25 //  gener->SetYRange(-999,999);
  26 //  AliGenReaderTreeK * reader = new AliGenReaderTreeK();
  27 //  reader->SetFileName("myFileWithTreeK.root");
  28 //  gener->SetReader(reader);
  29 //  gener->Init();
  30
  31
  32 #include <Riostream.h>
  33
  34 #include "AliGenExtFile.h"
  35 #include "AliRunLoader.h"
  36 #include "AliHeader.h"
  37 #include "AliGenEventHeader.h"
  38
  39 #include <TParticle.h>
  40 #include <TFile.h>
  41 #include <TTree.h>
  42
  43
  44 ClassImp(AliGenExtFile)
  45
  46 AliGenExtFile::AliGenExtFile()
  47     :AliGenMC(),
  48      fFileName(0),
  49      fReader(0)
  50 {
  51 //  Constructor
  52 //
  53 //  Read all particles
  54     fNpart  =- 1;
  55 }
  56
  57 AliGenExtFile::AliGenExtFile(Int_t npart)
  58     :AliGenMC(npart),
  59      fFileName(0),
  60      fReader(0)
  61 {
  62 //  Constructor
  63     fName   = "ExtFile";
  64     fTitle  = "Primaries from ext. File";
  65 }
  66
  67 //____________________________________________________________
  68 AliGenExtFile::~AliGenExtFile()
  69 {
  70 // Destructor
  71     delete fReader;
  72 }
  73
  74 //___________________________________________________________
  75 void AliGenExtFile::Init()
  76 {
  77 // Initialize
  78     if (fReader) fReader->Init();
  79 }
  80
  81
  82 void AliGenExtFile::Generate()
  83 {
  84 // Generate particles
  85
  86   Float_t polar[3]  = {0,0,0};
  87   //
  88   Float_t origin[3] = {0,0,0};
  89   Float_t p[3];
  90   Float_t random[6];
  91   Int_t i = 0, j, nt;
  92   //
  93   //
  94   if (fVertexSmear == kPerEvent) Vertex();
  95
  96   while(1) {
  97     Int_t nTracks = fReader->NextEvent();
  98     if (nTracks == 0) {
  99       // printf("\n No more events !!! !\n");
 100       Warning("AliGenExtFile::Generate","\nNo more events in external file!!!\nLast event may be empty or incomplete.\n");
 101       return;
 102     }
 103
 104     //
 105     // Particle selection loop
 106     //
 107     // The selection criterium for the external file generator is as follows:
 108     //
 109     // 1) All tracks are subject to the cuts defined by AliGenerator, i.e.
 110     //    fThetaMin, fThetaMax, fPhiMin, fPhiMax, fPMin, fPMax, fPtMin, fPtMax,
 111     //    fYMin, fYMax.
 112     //    If the particle does not satisfy these cuts, it is not put on the
 113     //    stack.
 114     // 2) If fCutOnChild and some specific child is selected (e.g. if
 115     //    fForceDecay==kSemiElectronic) the event is rejected if NOT EVEN ONE
 116     //    child falls into the child-cuts.
 117     TParticle* iparticle = 0x0;
 118
 119     if(fCutOnChild) {
 120       // Count the selected children
 121       Int_t nSelected = 0;
 122       while ((iparticle=fReader->NextParticle()) ) {
 123           Int_t kf = CheckPDGCode(iparticle->GetPdgCode());
 124           kf = TMath::Abs(kf);
 125           if (ChildSelected(kf) && KinematicSelection(iparticle, 1)) {
 126               nSelected++;
 127           }
 128       }
 129       if (!nSelected) continue;    // No particle selected:  Go to next event
 130       fReader->RewindEvent();
 131     }
 132
 133     //
 134     // Stack filling loop
 135     //
 136     fNprimaries = 0;
 137     for (i = 0; i < nTracks; i++) {
 138         TParticle* jparticle = fReader->NextParticle();
 139         Bool_t selected = KinematicSelection(jparticle,0);
 140         if (!selected) continue;
 141         p[0] = jparticle->Px();
 142         p[1] = jparticle->Py();
 143         p[2] = jparticle->Pz();
 144         Int_t idpart = iparticle->GetPdgCode();
 145         if(fVertexSmear==kPerTrack)
 146         {
 147             Rndm(random,6);
 148             for (j = 0; j < 3; j++) {
 149                 origin[j]=fOrigin[j]+
 150                     fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
 151                     TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 152             }
 153         } else {
 154             origin[0] = fVertex[0] + jparticle->Vx();
 155             origin[1] = fVertex[1] + jparticle->Vy();
 156             origin[2] = fVertex[2] + jparticle->Vz();
 157         }
 158
 159         Int_t decayed    = jparticle->GetFirstDaughter();
 160         Int_t doTracking = fTrackIt && (decayed < 0) && (TMath::Abs(idpart) > 10) && selected;
 161         Int_t parent     = jparticle->GetFirstMother();
 162
 163         PushTrack(doTracking, parent, idpart, p, origin, polar, 0, kPPrimary, nt);
 164         KeepTrack(nt);
 165         fNprimaries++;
 166     } // track loop
 167
 168     // Generated event header
 169
 170     AliGenEventHeader * header = new AliGenEventHeader();
 171     header->SetNProduced(fNprimaries);
 172     header->SetPrimaryVertex(fVertex);
 173     AddHeader(header);
 174     break;
 175
 176   } // event loop
 177
 178   SetHighWaterMark(nt);
 179   CdEventFile();
 180 }
 181
 182 void AliGenExtFile::CdEventFile()
 183 {
 184 // CD back to the event file
 185     (AliRunLoader::GetRunLoader())->CdGAFile();
 186 }
 187
 188
 189
 190
 191