STARLIGHT/AliStarLight/AliGenStarLight.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 #include <Riostream.h>
  18 #include "AliLog.h"
  19 #include "TStarLight.h"
  20 #include "starlight.h"
  21 #include "upcevent.h"
  22 #include "AliRunLoader.h"
  23 #include "AliGenEventHeader.h"
  24 #include "AliGenStarLight.h"
  25
  26 ClassImp(AliGenStarLight);
  27
  28 AliGenStarLight::AliGenStarLight()
  29   : AliGenMC()
  30   , fRapidityMotherMin( 1) // Max < Min: no cut
  31   , fRapidityMotherMax(-1)
  32   , fEtaChildMin( 1)       // Max < Min: no cut
  33   , fEtaChildMax(-1)
  34   , fSLgenerator(NULL)
  35 {
  36 }
  37 //----------------------------------------------------------------------
  38 AliGenStarLight::AliGenStarLight(Int_t npart)
  39   : AliGenMC(npart)
  40   , fRapidityMotherMin( 1) // Max < Min: no cut
  41   , fRapidityMotherMax(-1)
  42   , fEtaChildMin( 1)       // Max < Min: no cut
  43   , fEtaChildMax(-1)
  44   , fSLgenerator(new TStarLight("TStarLight",
  45                                 "StarLight UPC Event Generator",
  46                                 ""))  // no config file name
  47 {
  48   //
  49 }
  50 //----------------------------------------------------------------------
  51 AliGenStarLight::~AliGenStarLight() {
  52   if (NULL != fSLgenerator) delete fSLgenerator;
  53   fSLgenerator = NULL;
  54 }
  55 void AliGenStarLight::ImportConfigurationFromFile(const char* filename) {
  56   if (NULL == fSLgenerator) {
  57     AliFatal("AliGenStarLight class not constructed properly. ");
  58     return;
  59   }
  60   fSLgenerator->ImportConfigurationFromFile(filename);
  61 }
  62 void AliGenStarLight::SetParameter(const char* line) {
  63   if (NULL == fSLgenerator) {
  64     AliFatal("AliGenStarLight class not constructed properly. ");
  65     return;
  66   }
  67   fSLgenerator->SetParameter(line);
  68 }
  69 //----------------------------------------------------------------------
  70 void AliGenStarLight::Init() {
  71   if (NULL == fSLgenerator) {
  72     AliFatal("AliGenStarLight class not constructed properly. ");
  73     return;
  74   }
  75   fSLgenerator->InitStarLight();
  76 }
  77 //----------------------------------------------------------------------
  78 void AliGenStarLight::Generate() {
  79   Float_t vpos[4] = { 0, 0, 0, 0 };
  80   if (fVertexSmear == kPerEvent) {
  81     Vertex(); // get vertex
  82     for (Int_t i=0; i<3; ++i)
  83       vpos[i] = fVertex[i];
  84     vpos[3] = fTime;
  85   }
  86
  87   Int_t   nt(0);     // number of tracks
  88   Bool_t  genOK(kFALSE);
  89   // generate events until all constraints are fulfilled
  90   for (Int_t trials=0; !genOK && trials < 100*1000; ++trials) {
  91     fSLgenerator->GenerateEvent();
  92     fSLgenerator->BoostEvent();
  93     fSLgenerator->ImportParticles(&fParticles, "ALL");
  94
  95     TLorentzVector vSum;
  96     genOK = kTRUE;
  97     const Long64_t n(fParticles.GetEntries());
  98     if (n == 0) {
  99       AliFatal("no particles generated");
 100       return;
 101     }
 102     for (Long64_t i(0); i<n; ++i) {
 103       TParticle *part(dynamic_cast<TParticle*>(fParticles.At(i)));
 104       if (NULL == part) {
 105         AliFatal("NULL == part");
 106         return;
 107       }
 108       genOK =
 109         (part->Theta() >= fThetaMin) && (part->Theta() <  fThetaMax) &&
 110         (part->Phi()   >= fPhiMin)   && (part->Phi()   <  fPhiMax)   &&
 111         (part->Y()     >= fYMin)     && (part->Y()     <  fYMax)     &&
 112         (part->P()     >= fPMin)     && (part->P()     <  fPMax)     &&
 113         (part->Pt()    >= fPtMin)    && (part->Pt()    <  fPtMax);
 114       if (fEtaChildMin <= fEtaChildMax) // no cut if Max < Min
 115         genOK = genOK && (part->Eta() >= fEtaChildMin &&
 116                           part->Eta() <  fEtaChildMax);
 117       if (kFALSE == genOK)
 118         break;
 119
 120       TLorentzVector v;
 121       part->Momentum(v);
 122       vSum += v;
 123     }
 124
 125     if (fRapidityMotherMin <= fRapidityMotherMax) // no cut if Max < Min
 126       genOK = (genOK
 127                ? (vSum.Rapidity() > fRapidityMotherMin &&
 128                   vSum.Rapidity() < fRapidityMotherMax)
 129                : kFALSE);
 130
 131     if (kFALSE == genOK) continue;
 132
 133     fNprimaries = 0;
 134     for (Long64_t i(0), n(fParticles.GetEntries()); i<n; ++i) {
 135       const TParticle *part(dynamic_cast<TParticle*>(fParticles.At(i)));
 136       if (NULL == part) {
 137         AliFatal("NULL == part");
 138         return;
 139       }
 140       const Int_t   iparent(-1);
 141       const Float_t polar[3] = { 0, 0, 0 };
 142       const Float_t weight(trials+1);
 143       PushTrack(fTrackIt, iparent, part->GetPdgCode(),
 144                 part->Px(), part->Py(), part->Pz(), part->Energy(),
 145                 vpos[0],    vpos[1],    vpos[2],    vpos[3],
 146                 polar[0],   polar[1],   polar[2],
 147                 kPPrimary, nt, weight, part->GetStatusCode());
 148       AliInfo(Form("weight=%.0f nt=%d fTrackIt=%d statusCode=%d",
 149                    weight, nt, fTrackIt, part->GetStatusCode()));
 150       part->Print();
 151       KeepTrack(nt);
 152       ++fNprimaries;
 153     }
 154     fParticles.Clear();
 155   }
 156   if (kFALSE == genOK)
 157     AliFatal("Maximum number of trials reached");
 158
 159   AliGenEventHeader* header = new AliGenEventHeader("SL");
 160   const TArrayF vertexPosition(3, vpos);
 161   header->SetPrimaryVertex(vertexPosition);
 162   header->SetInteractionTime(vpos[3]);
 163   header->SetNProduced(nt);
 164   AddHeader(header);
 165   SetHighWaterMark(nt);
 166   AliRunLoader::Instance()->CdGAFile();
 167 }