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