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