THerwig/AliGenHerwig.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
  19 // Generator using Herwig as an external generator
  20 // The main Herwig options are accessable for the user through this interface.
  21 // Uses the THerwig implementation of TGenerator.
  22
  23 #include "AliGenHerwig.h"
  24 #include "AliRun.h"
  25
  26 #include <TParticle.h>
  27 #include "THerwig6.h"
  28
  29 #include "Riostream.h"
  30
  31  ClassImp(AliGenHerwig)
  32
  33 static TRandom * sRandom;
  34
  35 AliGenHerwig::AliGenHerwig()
  36 {
  37 // Constructor
  38 }
  39
  40 AliGenHerwig::AliGenHerwig(Int_t npart)
  41     :AliGenMC(npart)
  42 {
  43     SetBeamMomenta();
  44     SetTarget();
  45     SetProjectile();
  46     SetStrucFunc(kGRVLO98);
  47     fKeep=0;
  48     fTrigger=0;
  49     fDecaysOff=1;
  50     fSelectAll=0;
  51     fFlavor=0;
  52     fPtHardMin=10.;
  53     fPtRMS=0.0;
  54     fEnSoft=1.0;
  55     fMaxPr=1;
  56     fMaxErrors=1000;
  57 //  Set random number
  58     if (!sRandom) sRandom=fRandom;
  59 }
  60
  61 AliGenHerwig::AliGenHerwig(const AliGenHerwig & Herwig)
  62 {
  63 // copy constructor
  64 }
  65
  66
  67 AliGenHerwig::~AliGenHerwig()
  68 {
  69 // Destructor
  70 }
  71
  72 void AliGenHerwig::Init()
  73 {
  74 // Initialisation
  75   fTarget.Resize(8);
  76   fProjectile.Resize(8);
  77   SetMC(new THerwig6());
  78   fHerwig=(THerwig6*) fgMCEvGen;
  79   // initialize common blocks
  80   fHerwig->Initialize(fProjectile, fTarget, fMomentum1, fMomentum2, fProcess);
  81   // reset parameters according to user needs
  82   InitPDF();
  83   fHerwig->SetPTMIN(fPtHardMin);
  84   fHerwig->SetPTRMS(fPtRMS);
  85   fHerwig->SetMAXPR(fMaxPr);
  86   fHerwig->SetMAXER(fMaxErrors);
  87   fHerwig->SetENSOF(fEnSoft);
  88   // compute parameter dependent constants
  89   fHerwig->PrepareRun();
  90 }
  91
  92 void AliGenHerwig::InitPDF()
  93 {
  94   switch(fStrucFunc)
  95     {
  96     case kGRVLO:
  97       fModPDF=5;
  98       fAutPDF="GRV";
  99       break;
 100     case kGRVHO:
 101       fModPDF=6;
 102       fAutPDF="GRV";
 103       break;
 104     case kGRVLO98:
 105       fModPDF=12;
 106       fAutPDF="GRV";
 107       break;
 108     case kMRSDminus:
 109       fModPDF=31;
 110       fAutPDF="MRS";
 111       break;
 112     case kMRSD0:
 113       fModPDF=30;
 114       fAutPDF="MRS";
 115       break;
 116     case kMRSG:
 117       fModPDF=41;
 118       fAutPDF="MRS";
 119       break;
 120     case kMRSTcgLO:
 121       fModPDF=72;
 122       fAutPDF="MRS";
 123       break;
 124     case kCTEQ4M:
 125       fModPDF=34;
 126       fAutPDF="CTEQ";
 127       break;
 128     case kCTEQ5L:
 129       fModPDF=46;
 130       fAutPDF="CTEQ";
 131       break;
 132     default:
 133       cerr << "This structure function is not inplemented " << fStrucFunc << endl;
 134       break;
 135     }
 136   fAutPDF.Resize(20);
 137   fHerwig->SetMODPDF(1,fModPDF);
 138   fHerwig->SetMODPDF(2,fModPDF);
 139   fHerwig->SetAUTPDF(1,fAutPDF);
 140   fHerwig->SetAUTPDF(2,fAutPDF);
 141 }
 142
 143 void AliGenHerwig::Generate()
 144 {
 145   // Generate one event
 146
 147   Float_t polar[3] =   {0,0,0};
 148   Float_t origin[3]=   {0,0,0};
 149   Float_t origin0[3]=  {0,0,0};
 150   Float_t p[4], random[6];
 151
 152   static TClonesArray *particles;
 153   //  converts from mm/c to s
 154   const Float_t kconv=0.001/2.999792458e8;
 155   //
 156   Int_t nt=0;
 157   Int_t jev=0;
 158   Int_t j, kf, ks, imo;
 159   kf=0;
 160
 161   if(!particles) particles=new TClonesArray("TParticle",10000);
 162
 163   fTrials=0;
 164   for (j=0;j<3;j++) origin0[j]=fOrigin[j];
 165   if(fVertexSmear==kPerEvent) {
 166     Rndm(random,6);
 167     for (j=0;j<3;j++) {
 168       origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
 169         TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 170     }
 171   }
 172
 173   while(1)
 174     {
 175         fHerwig->GenerateEvent();
 176         fTrials++;
 177         fHerwig->ImportParticles(particles,"All");
 178         Int_t np = particles->GetEntriesFast()-1;
 179         if (np == 0 ) continue;
 180
 181         Int_t nc=0;
 182
 183         Int_t * newPos = new Int_t[np];
 184         for (Int_t i = 0; i<np; i++) *(newPos+i)=-1;
 185
 186         for (Int_t i = 0; i<np; i++) {
 187             TParticle *  iparticle       = (TParticle *) particles->At(i);
 188             imo = iparticle->GetFirstMother();
 189             kf        = iparticle->GetPdgCode();
 190             ks        = iparticle->GetStatusCode();
 191             if (ks != 3 &&
 192                 KinematicSelection(iparticle,0))
 193             {
 194                 nc++;
 195                 p[0]=iparticle->Px();
 196                 p[1]=iparticle->Py();
 197                 p[2]=iparticle->Pz();
 198                 p[3]=iparticle->Energy();
 199                 origin[0]=origin0[0]+iparticle->Vx()/10;
 200                 origin[1]=origin0[1]+iparticle->Vy()/10;
 201                 origin[2]=origin0[2]+iparticle->Vz()/10;
 202                 Float_t tof = kconv*iparticle->T();
 203                 Int_t   iparent = (imo > -1) ? newPos[imo] : -1;
 204                 Int_t   trackIt = (ks == 1) && fTrackIt;
 205                 gAlice->PushTrack(trackIt, iparent, kf,
 206                                  p[0], p[1], p[2], p[3],
 207                                  origin[0], origin[1], origin[2],
 208                                  tof,
 209                                  polar[0], polar[1], polar[2],
 210                                  kPPrimary, nt, 1., ks);
 211                 KeepTrack(nt);
 212                 newPos[i]=nt;
 213             } // end of if: selection of particle
 214         } // end of for: particle loop
 215         if (newPos) delete[] newPos;
 216         printf("\n I've put %i particles on the stack \n",nc);
 217         //      MakeHeader();
 218         printf("nc: %d %d\n", nc, fNpart);
 219
 220         if (nc > 0) {
 221             jev+=nc;
 222             if (jev >= fNpart || fNpart == -1) {
 223                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
 224                 printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
 225                 break;
 226             }
 227         }
 228     }
 229   SetHighWaterMark(nt);
 230 //  adjust weight due to kinematic selection
 231   AdjustWeights();
 232 //  get cross-section
 233   fXsection=fHerwig->GetAVWGT();
 234 }
 235
 236 void AliGenHerwig::AdjustWeights()
 237 {
 238 // Adjust the weights after generation of all events
 239     TParticle *part;
 240     Int_t ntrack=gAlice->GetNtrack();
 241     for (Int_t i=0; i<ntrack; i++) {
 242         part= gAlice->Particle(i);
 243         part->SetWeight(part->GetWeight()*fKineBias);
 244     }
 245 }
 246
 247
 248 void AliGenHerwig::KeepFullEvent()
 249 {
 250     fKeep=1;
 251 }
 252
 253 Bool_t AliGenHerwig::DaughtersSelection(TParticle* iparticle, TClonesArray* particles)
 254 {
 255 //
 256 // Looks recursively if one of the daughters has been selected
 257 //
 258 //    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
 259     Int_t imin=-1;
 260     Int_t imax=-1;
 261     Int_t i;
 262     Bool_t hasDaughters= (iparticle->GetFirstDaughter() >=0);
 263     Bool_t selected=kFALSE;
 264     if (hasDaughters) {
 265         imin=iparticle->GetFirstDaughter();
 266         imax=iparticle->GetLastDaughter();
 267         for (i=imin; i<= imax; i++){
 268             TParticle *  jparticle       = (TParticle *) particles->At(i);
 269             Int_t ip=jparticle->GetPdgCode();
 270             if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
 271                 selected=kTRUE; break;
 272             }
 273             if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
 274         }
 275     } else {
 276         return kFALSE;
 277     }
 278
 279     return selected;
 280 }
 281
 282
 283 Bool_t AliGenHerwig::SelectFlavor(Int_t pid)
 284 {
 285 // Select flavor of particle
 286 // 0: all
 287 // 4: charm and beauty
 288 // 5: beauty
 289     if (fFlavor == 0) return kTRUE;
 290
 291     Int_t ifl=TMath::Abs(pid/100);
 292     if (ifl > 10) ifl/=10;
 293     return (fFlavor == ifl);
 294 }
 295
 296 Bool_t AliGenHerwig::Stable(TParticle*  particle)
 297 {
 298 // Return true for a stable particle
 299 //
 300     Int_t kf = TMath::Abs(particle->GetPdgCode());
 301
 302     if ( (particle->GetFirstDaughter() < 0 ) || (kf == 1000*fFlavor+122))
 303
 304     {
 305         return kTRUE;
 306     } else {
 307         return kFALSE;
 308     }
 309 }
 310
 311 void AliGenHerwig::FinishRun()
 312 {
 313   fHerwig->Hwefin();
 314 }
 315
 316
 317 AliGenHerwig& AliGenHerwig::operator=(const  AliGenHerwig& rhs)
 318 {
 319 // Assignment operator
 320     return *this;
 321 }
 322
 323
 324 extern "C" {
 325   Double_t hwr_() {return sRandom->Rndm();}
 326 }