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 "AliHerwigRndm.h"
  25 #include "AliRun.h"
  26
  27 #include <TParticle.h>
  28 #include "THerwig6.h"
  29
  30 #include "Riostream.h"
  31 #include "AliMC.h"
  32
  33 #include "driver.h"
  34
  35 ClassImp(AliGenHerwig)
  36
  37
  38   AliGenHerwig::AliGenHerwig() :
  39     AliGenMC(),
  40     fAutPDF("LHAPDF"),
  41     fModPDF(19070),
  42     fStrucFunc(kCTEQ5L),
  43     fKeep(0),
  44     fDecaysOff(1),
  45     fTrigger(0),
  46     fSelectAll(0),
  47     fFlavor(0),
  48     fEnergyCMS(14000),
  49     fMomentum1(7000),
  50     fMomentum2(7000),
  51     fKineBias(1),
  52     fTrials(0),
  53     fXsection(0),
  54     fHerwig(0x0),
  55     fProcess(0),
  56     fPtHardMin(0.),
  57     fPtRMS(0.),
  58     fMaxPr(10),
  59     fMaxErrors(1000),
  60     fEnSoft(1),
  61     fEv1Pr(0),
  62     fEv2Pr(0),
  63     fFileName(0)
  64 {
  65 // Constructor
  66 }
  67
  68 AliGenHerwig::AliGenHerwig(Int_t npart)
  69     :AliGenMC(npart),
  70     fAutPDF("LHAPDF"),
  71     fModPDF(19070),
  72     fStrucFunc(kCTEQ5L),
  73     fKeep(0),
  74     fDecaysOff(1),
  75     fTrigger(0),
  76     fSelectAll(0),
  77     fFlavor(0),
  78     fEnergyCMS(14000),
  79     fMomentum1(7000),
  80     fMomentum2(7000),
  81     fKineBias(1),
  82     fTrials(0),
  83     fXsection(0),
  84     fHerwig(0x0),
  85     fProcess(0),
  86     fPtHardMin(0.),
  87     fPtRMS(0.),
  88     fMaxPr(10),
  89     fMaxErrors(1000),
  90     fEnSoft(1),
  91     fEv1Pr(0),
  92     fEv2Pr(0),
  93     fFileName(0)
  94 {
  95     SetTarget();
  96     SetProjectile();
  97     // Set random number generator
  98     AliHerwigRndm::SetHerwigRandom(GetRandom());
  99 }
 100
 101 AliGenHerwig::AliGenHerwig(const AliGenHerwig & Herwig)
 102     :AliGenMC(Herwig),
 103     fAutPDF("LHAPDF"),
 104     fModPDF(19070),
 105     fStrucFunc(kCTEQ5L),
 106     fKeep(0),
 107     fDecaysOff(1),
 108     fTrigger(0),
 109     fSelectAll(0),
 110     fFlavor(0),
 111     fEnergyCMS(14000),
 112     fMomentum1(7000),
 113     fMomentum2(7000),
 114     fKineBias(1),
 115     fTrials(0),
 116     fXsection(0),
 117     fHerwig(0x0),
 118     fProcess(0),
 119     fPtHardMin(0.),
 120     fPtRMS(0.),
 121     fMaxPr(10),
 122     fMaxErrors(1000),
 123     fEnSoft(1),
 124     fEv1Pr(0),
 125     fEv2Pr(0),
 126     fFileName(0)
 127 {
 128 // Copy constructor
 129     Herwig.Copy(*this);
 130 }
 131
 132
 133 AliGenHerwig::~AliGenHerwig()
 134 {
 135 // Destructor
 136 }
 137
 138 void AliGenHerwig::SetEventListRange(Int_t eventFirst, Int_t eventLast)
 139 {
 140   fEv1Pr = ++eventFirst;
 141   fEv2Pr = ++eventLast;
 142   if ( fEv2Pr == -1 ) fEv2Pr = fEv2Pr;
 143 }
 144
 145 void AliGenHerwig::Init()
 146 {
 147 // Initialisation
 148   fTarget.Resize(8);
 149   fProjectile.Resize(8);
 150   SetMC(new THerwig6());
 151   fHerwig=(THerwig6*) fMCEvGen;
 152   // initialize common blocks
 153   fHerwig->Initialize(fProjectile, fTarget, fMomentum1, fMomentum2, fProcess);
 154   // reset parameters according to user needs
 155   InitPDF();
 156   fHerwig->SetPTMIN(fPtHardMin);
 157   fHerwig->SetPTRMS(fPtRMS);
 158   fHerwig->SetMAXPR(fMaxPr);
 159   fHerwig->SetMAXER(fMaxErrors);
 160   fHerwig->SetENSOF(fEnSoft);
 161
 162   fHerwig->SetEV1PR(fEv1Pr);
 163   fHerwig->SetEV2PR(fEv2Pr);
 164
 165 // C---D,U,S,C,B,T QUARK AND GLUON MASSES (IN THAT ORDER)
 166 //       RMASS(1)=0.32
 167 //       RMASS(2)=0.32
 168 //       RMASS(3)=0.5
 169 //       RMASS(4)=1.55
 170 //       RMASS(5)=4.75
 171 //       RMASS(6)=174.3
 172 //       RMASS(13)=0.75
 173
 174   fHerwig->SetRMASS(4,1.2);
 175   fHerwig->SetRMASS(5,4.75);
 176
 177   if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(),50);
 178
 179   fHerwig->Hwusta("PI0     ");
 180
 181   // compute parameter dependent constants
 182   fHerwig->PrepareRun();
 183 }
 184
 185 void AliGenHerwig::InitJimmy()
 186 {
 187 // Initialisation
 188   fTarget.Resize(8);
 189   fProjectile.Resize(8);
 190   SetMC(new THerwig6());
 191   fHerwig=(THerwig6*) fMCEvGen;
 192   // initialize common blocks
 193   fHerwig->InitializeJimmy(fProjectile, fTarget, fMomentum1, fMomentum2, fProcess);
 194   // reset parameters according to user needs
 195   InitPDF();
 196   fHerwig->SetPTMIN(fPtHardMin);
 197   fHerwig->SetPTRMS(fPtRMS);
 198   fHerwig->SetMAXPR(fMaxPr);
 199   fHerwig->SetMAXER(fMaxErrors);
 200   fHerwig->SetENSOF(fEnSoft);
 201
 202   fHerwig->SetEV1PR(fEv1Pr);
 203   fHerwig->SetEV2PR(fEv2Pr);
 204
 205 // C---D,U,S,C,B,T QUARK AND GLUON MASSES (IN THAT ORDER)
 206 //       RMASS(1)=0.32
 207 //       RMASS(2)=0.32
 208 //       RMASS(3)=0.5
 209 //       RMASS(4)=1.55
 210 //       RMASS(5)=4.75
 211 //       RMASS(6)=174.3
 212 //       RMASS(13)=0.75
 213
 214   fHerwig->SetRMASS(4,1.2);
 215   fHerwig->SetRMASS(5,4.75);
 216
 217   if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(),50);
 218
 219   fHerwig->Hwusta("PI0     ");
 220
 221   // compute parameter dependent constants
 222   fHerwig->PrepareRunJimmy();
 223 }
 224
 225 void AliGenHerwig::InitPDF()
 226 {
 227   switch(fStrucFunc)
 228     {
 229 // ONLY USES LHAPDF STRUCTURE FUNCTIONS
 230     case kGRVLO98:
 231       fModPDF=80060;
 232       fAutPDF="HWLHAPDF";
 233       break;
 234     case kCTEQ6:
 235       fModPDF=10040;
 236       fAutPDF="HWLHAPDF";
 237       break;
 238     case kCTEQ61:
 239       fModPDF=10100;
 240       fAutPDF="HWLHAPDF";
 241       break;
 242     case kCTEQ6m:
 243       fModPDF=10050;
 244       fAutPDF="HWLHAPDF";
 245       break;
 246     case kCTEQ6l:
 247       fModPDF=10041;
 248       fAutPDF="HWLHAPDF";
 249       break;
 250     case kCTEQ6ll:
 251       fModPDF=10042;
 252       fAutPDF="HWLHAPDF";
 253       break;
 254     case kCTEQ5M:
 255       fModPDF=19050;
 256       fAutPDF="HWLHAPDF";
 257       break;
 258     case kCTEQ5L:
 259       fModPDF=19070;
 260       fAutPDF="HWLHAPDF";
 261       break;
 262     case kCTEQ4M:
 263       fModPDF=19150;
 264       fAutPDF="HWLHAPDF";
 265       break;
 266     case kCTEQ4L:
 267       fModPDF=19170;
 268       fAutPDF="HWLHAPDF";
 269       break;
 270 //    case kMRST2004nlo:
 271 //      fModPDF=20400;
 272 //      fAutPDF="HWLHAPDF";
 273 //      break;
 274     default:
 275       cerr << "This structure function is not inplemented " << fStrucFunc << endl;
 276       break;
 277     }
 278   fAutPDF.Resize(20);
 279   fHerwig->SetMODPDF(1,fModPDF);
 280   fHerwig->SetMODPDF(2,fModPDF);
 281   fHerwig->SetAUTPDF(1,fAutPDF);
 282   fHerwig->SetAUTPDF(2,fAutPDF);
 283 }
 284
 285 void AliGenHerwig::Generate()
 286 {
 287   // Generate one event
 288
 289   Float_t polar[3] =   {0,0,0};
 290   Float_t origin[3]=   {0,0,0};
 291   Float_t origin0[3]=  {0,0,0};
 292   Float_t p[4], random[6];
 293
 294   static TClonesArray *particles;
 295   //  converts from mm/c to s
 296   const Float_t kconv=0.001/2.999792458e8;
 297   //
 298   Int_t nt=0;
 299   Int_t jev=0;
 300   Int_t j, kf, ks, imo;
 301   kf=0;
 302
 303   if(!particles) particles=new TClonesArray("TParticle",10000);
 304
 305   fTrials=0;
 306   for (j=0;j<3;j++) origin0[j]=fOrigin[j];
 307   if(fVertexSmear==kPerEvent) {
 308     Rndm(random,6);
 309     for (j=0;j<3;j++) {
 310       origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
 311         TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 312     }
 313   }
 314
 315   while(1)
 316     {
 317         fHerwig->GenerateEvent();
 318         fTrials++;
 319         fHerwig->ImportParticles(particles,"All");
 320         Int_t np = particles->GetEntriesFast()-1;
 321         if (np == 0 ) continue;
 322
 323         Int_t nc=0;
 324
 325         Int_t * newPos = new Int_t[np];
 326         for (Int_t i = 0; i<np; i++) *(newPos+i)=-1;
 327
 328         for (Int_t i = 0; i<np; i++) {
 329             TParticle *  iparticle       = (TParticle *) particles->At(i);
 330             imo = iparticle->GetFirstMother();
 331             kf        = iparticle->GetPdgCode();
 332             ks        = iparticle->GetStatusCode();
 333             if (ks != 3 &&
 334                 KinematicSelection(iparticle,0))
 335             {
 336                 nc++;
 337                 p[0]=iparticle->Px();
 338                 p[1]=iparticle->Py();
 339                 p[2]=iparticle->Pz();
 340                 p[3]=iparticle->Energy();
 341                 origin[0]=origin0[0]+iparticle->Vx()/10;
 342                 origin[1]=origin0[1]+iparticle->Vy()/10;
 343                 origin[2]=origin0[2]+iparticle->Vz()/10;
 344                 Float_t tof = kconv*iparticle->T();
 345                 Int_t   iparent = (imo > -1) ? newPos[imo] : -1;
 346                 Int_t   trackIt = (ks == 1) && fTrackIt;
 347                 PushTrack(trackIt, iparent, kf,
 348                           p[0], p[1], p[2], p[3],
 349                           origin[0], origin[1], origin[2],
 350                           tof,
 351                           polar[0], polar[1], polar[2],
 352                           kPPrimary, nt, fHerwig->GetEVWGT(), ks);
 353                 KeepTrack(nt);
 354                 newPos[i]=nt;
 355             } // end of if: selection of particle
 356         } // end of for: particle loop
 357         if (newPos) delete[] newPos;
 358         //      MakeHeader();
 359         if (nc > 0) {
 360             jev+=nc;
 361             if (jev >= fNpart || fNpart == -1) {
 362                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
 363                 break;
 364             }
 365         }
 366     }
 367   SetHighWaterMark(nt);
 368 //  adjust weight due to kinematic selection
 369   AdjustWeights();
 370 //  get cross-section
 371   fXsection=fHerwig->GetAVWGT();
 372 }
 373
 374 void AliGenHerwig::AdjustWeights()
 375 {
 376 // Adjust the weights after generation of all events
 377     TParticle *part;
 378     Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
 379     for (Int_t i=0; i<ntrack; i++) {
 380         part= gAlice->GetMCApp()->Particle(i);
 381         part->SetWeight(part->GetWeight()*fKineBias);
 382     }
 383 }
 384
 385
 386 void AliGenHerwig::KeepFullEvent()
 387 {
 388     fKeep=1;
 389 }
 390
 391 Bool_t AliGenHerwig::DaughtersSelection(TParticle* iparticle, TClonesArray* particles)
 392 {
 393 //
 394 // Looks recursively if one of the daughters has been selected
 395 //
 396 //    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
 397     Int_t imin=-1;
 398     Int_t imax=-1;
 399     Int_t i;
 400     Bool_t hasDaughters= (iparticle->GetFirstDaughter() >=0);
 401     Bool_t selected=kFALSE;
 402     if (hasDaughters) {
 403         imin=iparticle->GetFirstDaughter();
 404         imax=iparticle->GetLastDaughter();
 405         for (i=imin; i<= imax; i++){
 406             TParticle *  jparticle       = (TParticle *) particles->At(i);
 407             Int_t ip=jparticle->GetPdgCode();
 408             if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
 409                 selected=kTRUE; break;
 410             }
 411             if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
 412         }
 413     } else {
 414         return kFALSE;
 415     }
 416
 417     return selected;
 418 }
 419
 420
 421 Bool_t AliGenHerwig::SelectFlavor(Int_t pid)
 422 {
 423 // Select flavor of particle
 424 // 0: all
 425 // 4: charm and beauty
 426 // 5: beauty
 427     if (fFlavor == 0) return kTRUE;
 428
 429     Int_t ifl=TMath::Abs(pid/100);
 430     if (ifl > 10) ifl/=10;
 431     return (fFlavor == ifl);
 432 }
 433
 434 Bool_t AliGenHerwig::Stable(TParticle*  particle)
 435 {
 436 // Return true for a stable particle
 437 //
 438     Int_t kf = TMath::Abs(particle->GetPdgCode());
 439
 440     if ( (particle->GetFirstDaughter() < 0 ) || (kf == 1000*fFlavor+122))
 441
 442     {
 443         return kTRUE;
 444     } else {
 445         return kFALSE;
 446     }
 447 }
 448
 449 void AliGenHerwig::FinishRun()
 450 {
 451   fHerwig->Hwefin();
 452 }
 453
 454
 455 AliGenHerwig& AliGenHerwig::operator=(const  AliGenHerwig& rhs)
 456 {
 457 // Assignment operator
 458     rhs.Copy(*this);
 459     return (*this);
 460 }
 461
 462 void AliGenHerwig::FinishRunJimmy()
 463 {
 464   fHerwig->Hwefin();
 465   fHerwig->Jmefin();
 466
 467 }
 468