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