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 "AliLog.h"
  35 #include "AliGenExtFile.h"
  36 #include "AliRunLoader.h"
  37 #include "AliHeader.h"
  38 #include "AliStack.h"
  39 #include "AliGenEventHeader.h"
  40 #include "AliGenReader.h"
  41
  42 #include <TParticle.h>
  43 #include <TFile.h>
  44 #include <TTree.h>
  45
  46
  47 using std::cout;
  48 using std::endl;
  49 using std::map;
  50
  51 ClassImp(AliGenExtFile)
  52
  53 AliGenExtFile::AliGenExtFile()
  54     :AliGenMC(),
  55      fFileName(0),
  56      fReader(0),
  57      fStartEvent(0)
  58 {
  59 //  Constructor
  60 //
  61 //  Read all particles
  62     fNpart  = -1;
  63 }
  64
  65 AliGenExtFile::AliGenExtFile(Int_t npart)
  66     :AliGenMC(npart),
  67      fFileName(0),
  68      fReader(0),
  69      fStartEvent(0)
  70 {
  71 //  Constructor
  72     fName   = "ExtFile";
  73     fTitle  = "Primaries from ext. File";
  74 }
  75
  76 //____________________________________________________________
  77 AliGenExtFile::~AliGenExtFile()
  78 {
  79 // Destructor
  80     delete fReader;
  81 }
  82
  83 //___________________________________________________________
  84 void AliGenExtFile::Init()
  85 {
  86 // Initialize
  87     if (fReader) fReader->Init();
  88 }
  89
  90 //___________________________________________________________
  91 void AliGenExtFile::Generate()
  92 {
  93 // Generate particles
  94
  95   Double_t polar[3]  = {0,0,0};
  96   //
  97   Double_t origin[3] = {0,0,0};
  98   Double_t time = 0.;
  99   Double_t p[4];
 100   Float_t random[6];
 101   Int_t i=0, j, nt;
 102   //
 103   //
 104   if (fVertexSmear == kPerEvent) Vertex();
 105
 106   // Fast forward up to start Event
 107   for (Int_t ie=0; ie<fStartEvent; ++ie ) {
 108     Int_t nTracks = fReader->NextEvent();
 109     if (nTracks == 0) {
 110       // printf("\n No more events !!! !\n");
 111       Warning("AliGenExtFile::Generate","\nNo more events in external file!!!\nLast event may be empty or incomplete.\n");
 112       return;
 113     }
 114     for (i = 0; i < nTracks; ++i) {
 115       if (NULL == fReader->NextParticle())
 116         AliFatal("Error while skipping tracks");
 117     }
 118     cout << "Skipping event " << ie << endl;
 119   }
 120   fStartEvent = 0; // do not skip events the second time
 121
 122   while(1) {
 123     Int_t nTracks = fReader->NextEvent();
 124     if (nTracks == 0) {
 125       // printf("\n No more events !!! !\n");
 126       Warning("AliGenExtFile::Generate","\nNo more events in external file!!!\nLast event may be empty or incomplete.\n");
 127       return;
 128     }
 129
 130     //
 131     // Particle selection loop
 132     //
 133     // The selection criterium for the external file generator is as follows:
 134     //
 135     // 1) All tracks are subject to the cuts defined by AliGenerator, i.e.
 136     //    fThetaMin, fThetaMax, fPhiMin, fPhiMax, fPMin, fPMax, fPtMin, fPtMax,
 137     //    fYMin, fYMax.
 138     //    If the particle does not satisfy these cuts, it is not put on the
 139     //    stack.
 140     // 2) If fCutOnChild and some specific child is selected (e.g. if
 141     //    fForceDecay==kSemiElectronic) the event is rejected if NOT EVEN ONE
 142     //    child falls into the child-cuts.
 143     TParticle* iparticle = 0x0;
 144
 145     if(fCutOnChild) {
 146       // Count the selected children
 147       Int_t nSelected = 0;
 148       while ((iparticle=fReader->NextParticle()) ) {
 149           Int_t kf = CheckPDGCode(iparticle->GetPdgCode());
 150           kf = TMath::Abs(kf);
 151           if (ChildSelected(kf) && KinematicSelection(iparticle, 1)) {
 152               nSelected++;
 153           }
 154       }
 155       if (!nSelected) continue;    // No particle selected:  Go to next event
 156       fReader->RewindEvent();
 157     }
 158
 159     //
 160     // Stack selection loop
 161     //
 162     class SelectorLogic { // need to do recursive back tracking, requires a "nested" function
 163     private:
 164        Int_t idCount;
 165        map<Int_t,Int_t> firstMotherMap;
 166        map<Int_t,Int_t> secondMotherMap;
 167        map<Int_t,Bool_t> selectedIdMap;
 168        map<Int_t,Int_t> newIdMap;
 169        void selectMothersToo(Int_t particleId) {
 170          Int_t mum1 = firstMotherMap[particleId];
 171           if (mum1 > -1 && !selectedIdMap[mum1]) {
 172              selectedIdMap[mum1] = true;
 173              selectMothersToo(mum1);
 174           }
 175           Int_t mum2 = secondMotherMap[particleId];
 176           if (mum2 > -1 && !selectedIdMap[mum2]) {
 177              selectedIdMap[mum2] = true;
 178              selectMothersToo(mum2);
 179           }
 180        }
 181     public:
 182       SelectorLogic():idCount(0), firstMotherMap(), secondMotherMap(), selectedIdMap(), newIdMap() {}
 183       void init() {
 184           idCount = 0;
 185        }
 186        void setData(Int_t id, Int_t mum1, Int_t mum2, Bool_t selected) {
 187           idCount++; // we know that this function is called in succession of ids, so counting is fine to determine max id
 188           firstMotherMap[id] = mum1;
 189           secondMotherMap[id] = mum2;
 190           selectedIdMap[id] = selected;
 191        }
 192        void reselectCuttedMothersAndRemapIDs() {
 193           for (Int_t id = 0; id < idCount; ++id) {
 194              if (selectedIdMap[id]) {
 195                 selectMothersToo(id);
 196              }
 197           }
 198           Int_t newId0 = 0;
 199           for (Int_t id = 0; id < idCount; id++) {
 200              if (selectedIdMap[id]) {
 201                 newIdMap[id] = newId0; ++newId0;
 202              } else {
 203                 newIdMap[id] = -1;
 204              }
 205           }
 206        }
 207        Bool_t isSelected(Int_t id) {
 208           return selectedIdMap[id];
 209        }
 210        Int_t newId(Int_t id) {
 211           if (id == -1) return -1;
 212           return newIdMap[id];
 213        }
 214     };
 215     SelectorLogic selector;
 216     selector.init();
 217     for (i = 0; i < nTracks; i++) {
 218        TParticle* jparticle = fReader->NextParticle();
 219        selector.setData(i,
 220              jparticle->GetFirstMother(),
 221              jparticle->GetSecondMother(),
 222              KinematicSelection(jparticle,0));
 223     }
 224     selector.reselectCuttedMothersAndRemapIDs();
 225     fReader->RewindEvent();
 226
 227     //
 228     // Stack filling loop
 229     //
 230     fNprimaries = 0;
 231     for (i = 0; i < nTracks; i++) {
 232        TParticle* jparticle = fReader->NextParticle();
 233        Bool_t selected = selector.isSelected(i);
 234        if (!selected) {
 235           continue;
 236        }
 237        Int_t parent = selector.newId(jparticle->GetFirstMother());
 238 //       printf("particle %d -> %d, with mother %d -> %d\n", i, selector.newId(i), jparticle->GetFirstMother(), parent);
 239
 240         p[0] = jparticle->Px();
 241         p[1] = jparticle->Py();
 242         p[2] = jparticle->Pz();
 243         p[3] = jparticle->Energy();
 244
 245         Int_t idpart = jparticle->GetPdgCode();
 246         if(fVertexSmear==kPerTrack)
 247         {
 248             Rndm(random,6);
 249             for (j = 0; j < 3; j++) {
 250                 origin[j]=fOrigin[j]+
 251                     fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
 252                     TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 253             }
 254             Rndm(random,2);
 255             time = fTimeOrigin + fOsigma[2]/TMath::Ccgs()*
 256               TMath::Cos(2*random[0]*TMath::Pi())*
 257               TMath::Sqrt(-2*TMath::Log(random[1]));
 258         } else {
 259             origin[0] = fVertex[0] + jparticle->Vx();
 260             origin[1] = fVertex[1] + jparticle->Vy();
 261             origin[2] = fVertex[2] + jparticle->Vz();
 262             time = fTime + jparticle->T();
 263         }
 264         Int_t doTracking = fTrackIt && selected && (jparticle->TestBit(kTransportBit));
 265
 266         PushTrack(doTracking, parent, idpart,
 267                   p[0], p[1], p[2], p[3], origin[0], origin[1], origin[2], time,
 268                   polar[0], polar[1], polar[2],
 269                   kPPrimary, nt, 1., jparticle->GetStatusCode());
 270
 271         KeepTrack(nt);
 272         fNprimaries++;
 273     } // track loop
 274
 275     // Generated event header
 276     AliGenEventHeader * header = fReader->GetGenEventHeader();
 277     if (!header) header = new AliGenEventHeader();
 278     header->SetNProduced(fNprimaries);
 279     header->SetPrimaryVertex(fVertex);
 280     header->SetInteractionTime(fTime);
 281     AddHeader(header);
 282     break;
 283
 284   } // event loop
 285
 286   SetHighWaterMark(nt);
 287   CdEventFile();
 288 }
 289
 290 //___________________________________________________________
 291 void AliGenExtFile::CdEventFile()
 292 {
 293 // CD back to the event file
 294   AliRunLoader::Instance()->CdGAFile();
 295 }
 296
 297
 298
 299
 300