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 time0 = 0.;
 116   Float_t p[3];
 117   Float_t tof;
 118
 119 //  converts from mm/c to s
 120   const Float_t kconv = 0.001/2.999792458e8;
 121 //
 122   Int_t nt  = 0;
 123   Int_t j, kf, ks, imo;
 124   kf = 0;
 125
 126   for (j = 0;j < 3; j++) origin0[j] = fOrigin[j];
 127   time0 = fTimeOrigin;
 128   if(fVertexSmear == kPerEvent) {
 129       Vertex();
 130       for (j=0; j < 3; j++) origin0[j] = fVertex[j];
 131       time0 = fTime;
 132   }
 133
 134   while(1)
 135   {
 136 //    Generate one event
 137       fHydjet->GenerateEvent();
 138       fHydjet->ImportParticles(&fParticles,"All");
 139
 140       Int_t np = fParticles.GetEntriesFast();
 141       Int_t nc = 0;
 142       if (np == 0 ) continue;
 143       Int_t i;
 144       Int_t* pSelected  = new Int_t[np];
 145
 146       for (i = 0; i < np; i++) {
 147            pSelected[i] = 0;
 148       }
 149
 150 //      Get event vertex
 151       fVertex[0] = origin0[0];
 152       fVertex[1] = origin0[1];
 153       fVertex[2] = origin0[2];
 154       fTime = time0;
 155
 156 //
 157 // Now select the final state particles
 158 //
 159
 160       for(i = 0; i<np; i++){
 161          TParticle *  iparticle = (TParticle *) fParticles.At(i);
 162          // Is this a final state particle ?
 163          if (!Stable(iparticle)) continue;
 164          Bool_t selected = kTRUE;
 165          kf = iparticle->GetPdgCode();
 166          ks = iparticle->GetStatusCode();
 167          if(!fSelectAll) selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
 168          // Put particle on the stack if selected
 169          if(selected) {
 170            nc++;
 171            pSelected[i] = 1;
 172          } // selected
 173       } // particle loop final state
 174
 175
 176 //
 177 // Write particles to stack
 178 //
 179
 180       for(i = 0; i<np; i++) {
 181          TParticle *  iparticle = (TParticle *) fParticles.At(i);
 182          Bool_t  hasDaughter = (iparticle->GetFirstDaughter() > 0);
 183          if(pSelected[i]){
 184            kf = iparticle->GetPdgCode();
 185            ks = iparticle->GetStatusCode();
 186            p[0] = iparticle->Px();
 187            p[1] = iparticle->Py();
 188            p[2] = iparticle->Pz();
 189            origin[0] = fVertex[0]+iparticle->Vx()/10;
 190            origin[1] = fVertex[1]+iparticle->Vy()/10;
 191            origin[2] = fVertex[2]+iparticle->Vz()/10;
 192            tof   = fTime + kconv*iparticle->T();
 193
 194            imo = -1;
 195            //imo = iparticle->GetFirstMother();
 196            //if(imo == 0) imo =-1;
 197            Bool_t tFlag = (fTrackIt && !hasDaughter);
 198            PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
 199            KeepTrack(nt);
 200          } // if selected
 201       } // particle loop
 202       delete[] pSelected;
 203
 204       AliInfo(Form("\n I've put %i particles on the stack \n",nc));
 205       if (nc > 0) break;
 206   } // event loop
 207   MakeHeader();
 208 }
 209
 210
 211 Bool_t AliGenHydjet::SelectFlavor(Int_t pid)
 212 {
 213 // Select flavor of particle
 214 // 0: all
 215 // 4: charm and beauty
 216 // 5: beauty
 217     Bool_t res = 0;
 218     if (fFlavor == 0) {
 219         res = kTRUE;
 220     } else {
 221         Int_t ifl = TMath::Abs(pid/100);
 222         if (ifl > 10) ifl/=10;
 223         res = (fFlavor == ifl);
 224     }
 225     return res;
 226 }
 227
 228 Bool_t AliGenHydjet::Stable(TParticle*  particle) const
 229 {
 230 // Return true for a stable particle
 231     if (particle->GetFirstDaughter() == 0 )
 232     {
 233         return kTRUE;
 234     } else {
 235         return kFALSE;
 236     }
 237 }
 238
 239
 240 ///////////////////////////////////
 241 void AliGenHydjet::MakeHeader()
 242 {
 243 // Builds the event header, to be called after each event
 244     AliGenEventHeader* header = new AliGenHydjetEventHeader("Hydjet");
 245     ((AliGenHydjetEventHeader*) header)->SetNProduced(fHydjet->GetN());
 246     ((AliGenHydjetEventHeader*) header)->SetImpactParameter(fHydjet->GetBGEN());
 247     ((AliGenHydjetEventHeader*) header)->SetNbcol(fHydjet->GetNBCOL());
 248     ((AliGenHydjetEventHeader*) header)->SetNpart(fHydjet->GetNPART());
 249     ((AliGenHydjetEventHeader*) header)->SetNpyt(fHydjet->GetNPYT());
 250     ((AliGenHydjetEventHeader*) header)->SetNhyd(fHydjet->GetNHYD());
 251   // draw generation values
 252   cout << "*********Hyflow ********** " << endl;
 253   cout << " YTFL " << fHydjet->GetYTFL() << endl;
 254   cout << " YLFL " << fHydjet->GetYLFL() << endl;
 255   cout << " FPART " << fHydjet->GetFPART() << endl;
 256   cout << "*********Hyjpar ********* " << endl;
 257   cout << " GetNHSEL() " << fHydjet->GetNHSEL() << endl;
 258   cout << " GetPTMIN() " << fHydjet->GetPTMIN() << endl;
 259   cout << " GetNJET() " << fHydjet->GetNJET() << endl;
 260   cout << "*********Hyfpar ********* " << endl;
 261   cout << " GetBGEN() " << fHydjet->GetBGEN() << endl;
 262   cout << " GetNBCOL() " << fHydjet->GetNBCOL() << endl;
 263   cout << " GetNPART() " << fHydjet->GetNPART() << endl;
 264   cout << " GetNPYT() " << fHydjet->GetNPYT() << endl;
 265   cout << " GetNHYD() " << fHydjet->GetNHYD() << endl;
 266
 267
 268 // Event Vertex
 269     header->SetPrimaryVertex(fVertex);
 270     header->SetInteractionTime(fTime);
 271     AddHeader(header);
 272     fCollisionGeometry = (AliGenHydjetEventHeader*)  header;
 273 }
 274
 275 void AliGenHydjet::AddHeader(AliGenEventHeader* header)
 276 {
 277     // Passes header either to the container or to gAlice
 278     if (fContainer) {
 279         fContainer->AddHeader(header);
 280     } else {
 281         gAlice->SetGenEventHeader(header);
 282     }
 283 }
 284
 285
 286 void AliGenHydjet::Copy(TObject &) const
 287 {
 288   Fatal("Copy","Not implemented!\n");
 289 }
 290
 291 AliGenHydjet& AliGenHydjet::operator=(const  AliGenHydjet& rhs)
 292 {
 293     rhs.Copy(*this);
 294     return (*this);
 295 }
 296