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