THydjet/AliGenHydjet.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 // Generator using Hydjet as an external generator
  19 // The main Hydjet options are accessable for the user through this interface.
  20 // Uses the THydjet implementation of TGenerator.
  21 // Author:
  22 // Rafael Diaz Valdes    (rafael.diaz.valdes@cern.ch)
  23 //
  24 #include <Riostream.h>
  25 #include <THydjet.h>
  26 #include <TParticle.h>
  27 #include <TClonesArray.h>
  28
  29 #include "AliGenHydjet.h"
  30 #include "AliLog.h"
  31 #include "AliGenHydjetEventHeader.h"
  32 #include "AliRun.h"
  33 #include "AliPythiaRndm.h"
  34
  35 ClassImp(AliGenHydjet)
  36
  37 AliGenHydjet::AliGenHydjet(Int_t npart) :
  38   AliGenMC(npart),
  39     //initial parameters
  40   fFrame("CMS"),      // Reference frame
  41   fAtomicWeigth(207),        // Projectile-Target atomic weight
  42   fIfbtype(0),          // centrality type
  43   fFixImpactParam(0),          // fixed impact parameter
  44   fMinImpactParam(0.),         // minimum impact parameter
  45   fMaxImpactParam(1.),         // maximum impact parameter
  46   fSoftMult(20000),      // mean soft multiplicity
  47   //hydro parameters
  48   fJetProd(2),          // flag Jet production (nhsel)
  49   fYflow(5.),         // max longitudinal flow
  50   fTflow(1.),         // max transverse flow
  51   fSoftFract(1.),         // Soft multiplicity fraction
  52   // PYTHIA parameters
  53   fDijetProd(1),          // flag dijet production
  54   fMinPtHard(10.),        // min pt hard
  55   fStructFunction(7),          // Structure Function (default CTEQ5M)
  56   fMultipleInt(0),          // flag multiple interaction
  57
  58   fHydjet(0),
  59
  60   fSelectAll(0),
  61   fFlavor(0)
  62
  63 {
  64 // Default PbPb collisions at 5. 5 TeV
  65 //
  66   fEnergyCMS = 5500.;      //Energy cms
  67   fName = "Hydjet";
  68   fTitle = "Particle Generator using Hydjet";
  69   // Set random number generator
  70   if (!AliPythiaRndm::GetPythiaRandom())
  71     AliPythiaRndm::SetPythiaRandom(GetRandom());
  72 }
  73
  74 AliGenHydjet::~AliGenHydjet()
  75 {
  76 // Destructor
  77 }
  78
  79 void AliGenHydjet::Init()
  80 {
  81 // Initialisation
  82
  83     SetMC(new THydjet(fEnergyCMS, fFrame,fAtomicWeigth,
  84                       fIfbtype,fMinImpactParam, fMaxImpactParam,
  85                       fFixImpactParam,fSoftMult));
  86
  87     fHydjet=(THydjet*) fMCEvGen;
  88
  89     // init hydro common blocks
  90     fHydjet->SetNHSEL(fJetProd);
  91     fHydjet->SetYLFL(fYflow);
  92     fHydjet->SetYTFL(fTflow);
  93     fHydjet->SetFPART(fSoftFract);
  94     // init PYTHIA common block
  95     fHydjet->SetMSEL(fDijetProd);
  96     fHydjet->SetPTMIN(fMinPtHard);
  97     fHydjet->SetCKIN(3,fMinPtHard);
  98     fHydjet->SetMSTP(51,fStructFunction);
  99     fHydjet->SetMSTP(81,0);
 100     fHydjet->SetMSTU(21,1);
 101     fHydjet->SetPARU(14,1.);
 102
 103 //  Initialize Hydjet
 104     fHydjet->Initialize();
 105
 106 }
 107
 108 void AliGenHydjet::Generate()
 109 {
 110
 111 ///////////////////////////////////////////////////////////////////////////////////////////
 112   Float_t polar[3]    =   {0,0,0};
 113   Float_t origin[3]   =   {0,0,0};
 114   Float_t origin0[3]  =   {0,0,0};
 115   Float_t p[3];
 116   Float_t tof;
 117
 118 //  converts from mm/c to s
 119   const Float_t kconv = 0.001/2.999792458e8;
 120 //
 121   Int_t nt  = 0;
 122   Int_t j, kf, ks, imo;
 123   kf = 0;
 124
 125   for (j = 0;j < 3; j++) origin0[j] = fOrigin[j];
 126   if(fVertexSmear == kPerEvent) {
 127       Vertex();
 128       for (j=0; j < 3; j++) origin0[j] = fVertex[j];
 129   }
 130
 131   while(1)
 132   {
 133 //    Generate one event
 134       fHydjet->GenerateEvent();
 135       fHydjet->ImportParticles(&fParticles,"All");
 136
 137       Int_t np = fParticles.GetEntriesFast();
 138       Int_t nc = 0;
 139       if (np == 0 ) continue;
 140       Int_t i;
 141       Int_t* pSelected  = new Int_t[np];
 142
 143       for (i = 0; i < np; i++) {
 144            pSelected[i] = 0;
 145       }
 146
 147 //      Get event vertex
 148       fVertex[0] = origin0[0];
 149       fVertex[1] = origin0[1];
 150       fVertex[2] = origin0[2];
 151
 152 //
 153 // Now select the final state particles
 154 //
 155
 156       for(i = 0; i<np; i++){
 157          TParticle *  iparticle = (TParticle *) fParticles.At(i);
 158          // Is this a final state particle ?
 159          if (!Stable(iparticle)) continue;
 160          Bool_t selected = kTRUE;
 161          kf = iparticle->GetPdgCode();
 162          ks = iparticle->GetStatusCode();
 163          if(!fSelectAll) selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
 164          // Put particle on the stack if selected
 165          if(selected) {
 166            nc++;
 167            pSelected[i] = 1;
 168          } // selected
 169       } // particle loop final state
 170
 171
 172 //
 173 // Write particles to stack
 174 //
 175
 176       for(i = 0; i<np; i++) {
 177          TParticle *  iparticle = (TParticle *) fParticles.At(i);
 178          Bool_t  hasDaughter = (iparticle->GetFirstDaughter() > 0);
 179          if(pSelected[i]){
 180            kf = iparticle->GetPdgCode();
 181            ks = iparticle->GetStatusCode();
 182            p[0] = iparticle->Px();
 183            p[1] = iparticle->Py();
 184            p[2] = iparticle->Pz();
 185            origin[0] = fVertex[0]+iparticle->Vx()/10;
 186            origin[1] = fVertex[1]+iparticle->Vy()/10;
 187            origin[2] = fVertex[2]+iparticle->Vz()/10;
 188            tof   = kconv*iparticle->T();
 189
 190            imo = -1;
 191            //imo = iparticle->GetFirstMother();
 192            //if(imo == 0) imo =-1;
 193            Bool_t tFlag = (fTrackIt && !hasDaughter);
 194            PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
 195            KeepTrack(nt);
 196          } // if selected
 197       } // particle loop
 198       delete[] pSelected;
 199
 200       AliInfo(Form("\n I've put %i particles on the stack \n",nc));
 201       if (nc > 0) break;
 202   } // event loop
 203   MakeHeader();
 204 }
 205
 206
 207 Bool_t AliGenHydjet::SelectFlavor(Int_t pid)
 208 {
 209 // Select flavor of particle
 210 // 0: all
 211 // 4: charm and beauty
 212 // 5: beauty
 213     Bool_t res = 0;
 214     if (fFlavor == 0) {
 215         res = kTRUE;
 216     } else {
 217         Int_t ifl = TMath::Abs(pid/100);
 218         if (ifl > 10) ifl/=10;
 219         res = (fFlavor == ifl);
 220     }
 221     return res;
 222 }
 223
 224 Bool_t AliGenHydjet::Stable(TParticle*  particle) const
 225 {
 226 // Return true for a stable particle
 227     if (particle->GetFirstDaughter() == 0 )
 228     {
 229         return kTRUE;
 230     } else {
 231         return kFALSE;
 232     }
 233 }
 234
 235
 236 ///////////////////////////////////
 237 void AliGenHydjet::MakeHeader()
 238 {
 239 // Builds the event header, to be called after each event
 240     AliGenEventHeader* header = new AliGenHydjetEventHeader("Hydjet");
 241     ((AliGenHydjetEventHeader*) header)->SetNProduced(fHydjet->GetN());
 242     ((AliGenHydjetEventHeader*) header)->SetImpactParameter(fHydjet->GetBGEN());
 243     ((AliGenHydjetEventHeader*) header)->SetNbcol(fHydjet->GetNBCOL());
 244     ((AliGenHydjetEventHeader*) header)->SetNpart(fHydjet->GetNPART());
 245     ((AliGenHydjetEventHeader*) header)->SetNpyt(fHydjet->GetNPYT());
 246     ((AliGenHydjetEventHeader*) header)->SetNhyd(fHydjet->GetNHYD());
 247   // draw generation values
 248   cout << "*********Hyflow ********** " << endl;
 249   cout << " YTFL " << fHydjet->GetYTFL() << endl;
 250   cout << " YLFL " << fHydjet->GetYLFL() << endl;
 251   cout << " FPART " << fHydjet->GetFPART() << endl;
 252   cout << "*********Hyjpar ********* " << endl;
 253   cout << " GetNHSEL() " << fHydjet->GetNHSEL() << endl;
 254   cout << " GetPTMIN() " << fHydjet->GetPTMIN() << endl;
 255   cout << " GetNJET() " << fHydjet->GetNJET() << endl;
 256   cout << "*********Hyfpar ********* " << endl;
 257   cout << " GetBGEN() " << fHydjet->GetBGEN() << endl;
 258   cout << " GetNBCOL() " << fHydjet->GetNBCOL() << endl;
 259   cout << " GetNPART() " << fHydjet->GetNPART() << endl;
 260   cout << " GetNPYT() " << fHydjet->GetNPYT() << endl;
 261   cout << " GetNHYD() " << fHydjet->GetNHYD() << endl;
 262
 263
 264 // Event Vertex
 265     header->SetPrimaryVertex(fVertex);
 266     AddHeader(header);
 267     fCollisionGeometry = (AliGenHydjetEventHeader*)  header;
 268 }
 269
 270 void AliGenHydjet::AddHeader(AliGenEventHeader* header)
 271 {
 272     // Passes header either to the container or to gAlice
 273     if (fContainer) {
 274         fContainer->AddHeader(header);
 275     } else {
 276         gAlice->SetGenEventHeader(header);
 277     }
 278 }
 279
 280
 281 void AliGenHydjet::Copy(TObject &) const
 282 {
 283   Fatal("Copy","Not implemented!\n");
 284 }
 285
 286 AliGenHydjet& AliGenHydjet::operator=(const  AliGenHydjet& rhs)
 287 {
 288     rhs.Copy(*this);
 289     return (*this);
 290 }
 291