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