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
  24 #include <Riostream.h>
  25 #include <TClonesArray.h>
  26 #include <TParticle.h>
  27
  28 #include <THerwig6.h>
  29
  30 #include "AliGenHerwig.h"
  31 #include "AliGenHerwigEventHeader.h"
  32 #include "AliHerwigRndm.h"
  33 #include "AliMC.h"
  34 #include "AliRun.h"
  35 #include "driver.h"
  36
  37 using std::cerr;
  38 using std::endl;
  39
  40 ClassImp(AliGenHerwig)
  41
  42
  43   AliGenHerwig::AliGenHerwig() :
  44     AliGenMC(),
  45     fAutPDF("LHAPDF"),
  46     fModPDF(19070),
  47     fStrucFunc(kCTEQ5L),
  48     fKeep(0),
  49     fDecaysOff(1),
  50     fTrigger(0),
  51     fSelectAll(0),
  52     fFlavor(0),
  53     fMomentum1(7000),
  54     fMomentum2(7000),
  55     fKineBias(1),
  56     fTrials(0),
  57     fXsection(0),
  58     fHerwig(0x0),
  59     fProcess(0),
  60     fPtHardMin(0.),
  61     fPtHardMax(9999.),
  62     fPtRMS(0.),
  63     fMaxPr(10),
  64     fMaxErrors(1000),
  65     fEnSoft(1),
  66     fEv1Pr(0),
  67     fEv2Pr(0),
  68     fFileName(0),
  69     fEtaMinParton(-20.),
  70     fEtaMaxParton(20.),
  71     fPhiMinParton(0.),
  72     fPhiMaxParton(2.* TMath::Pi()),
  73     fEtaMinGamma(-20.),
  74     fEtaMaxGamma(20.),
  75     fPhiMinGamma(0.),
  76     fPhiMaxGamma(2. * TMath::Pi()),
  77     fHeader(0)
  78 {
  79 // Constructor
  80   fEnergyCMS = 14000;
  81 }
  82
  83 AliGenHerwig::AliGenHerwig(Int_t npart)
  84     :AliGenMC(npart),
  85     fAutPDF("LHAPDF"),
  86     fModPDF(19070),
  87     fStrucFunc(kCTEQ5L),
  88     fKeep(0),
  89     fDecaysOff(1),
  90     fTrigger(0),
  91     fSelectAll(0),
  92     fFlavor(0),
  93     fMomentum1(7000),
  94     fMomentum2(7000),
  95     fKineBias(1),
  96     fTrials(0),
  97     fXsection(0),
  98     fHerwig(0x0),
  99     fProcess(0),
 100     fPtHardMin(0.),
 101     fPtHardMax(9999.),
 102     fPtRMS(0.),
 103     fMaxPr(10),
 104     fMaxErrors(1000),
 105     fEnSoft(1),
 106     fEv1Pr(0),
 107     fEv2Pr(0),
 108     fFileName(0),
 109     fEtaMinParton(-20.),
 110     fEtaMaxParton(20.),
 111     fPhiMinParton(0.),
 112     fPhiMaxParton(2.* TMath::Pi()),
 113     fEtaMinGamma(-20.),
 114     fEtaMaxGamma(20.),
 115     fPhiMinGamma(0.),
 116     fPhiMaxGamma(2. * TMath::Pi()),
 117     fHeader(0)
 118 {
 119 // Constructor
 120     fEnergyCMS = 14000;
 121     SetTarget();
 122     SetProjectile();
 123     // Set random number generator
 124     AliHerwigRndm::SetHerwigRandom(GetRandom());
 125 }
 126
 127 AliGenHerwig::~AliGenHerwig()
 128 {
 129 // Destructor
 130 }
 131
 132 void AliGenHerwig::SetEventListRange(Int_t eventFirst, Int_t eventLast)
 133 {
 134   fEv1Pr = eventFirst;
 135   fEv2Pr = eventLast;
 136   if ( fEv2Pr == -1 ) fEv2Pr = fEv1Pr;
 137 }
 138
 139 void AliGenHerwig::Init()
 140 {
 141 // Initialisation
 142   fTarget.Resize(8);
 143   fProjectile.Resize(8);
 144   SetMC(new THerwig6());
 145   fHerwig=(THerwig6*) fMCEvGen;
 146   // initialize common blocks
 147   fHerwig->Initialize(fProjectile.Data(), fTarget.Data(), fMomentum1, fMomentum2, fProcess); //
 148   // reset parameters according to user needs
 149   InitPDF();
 150   fHerwig->SetPTMIN(fPtHardMin);
 151   fHerwig->SetPTMAX(fPtHardMax);
 152   fHerwig->SetPTRMS(fPtRMS);
 153   fHerwig->SetMAXPR(fMaxPr);
 154   fHerwig->SetMAXER(fMaxErrors);
 155   fHerwig->SetENSOF(fEnSoft);
 156 // C---D,U,S,C,B,T QUARK AND GLUON MASSES (IN THAT ORDER)
 157 //       RMASS(1)=0.32
 158 //       RMASS(2)=0.32
 159 //       RMASS(3)=0.5
 160 //       RMASS(4)=1.55
 161 //       RMASS(5)=4.75
 162 //       RMASS(6)=174.3
 163 //       RMASS(13)=0.75
 164
 165   fHerwig->SetRMASS(4,1.2);
 166   fHerwig->SetRMASS(5,4.75);
 167
 168   if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(), 49);
 169
 170   //fHerwig->Hwusta("PI0     ");
 171
 172   // compute parameter dependent constants
 173   fHerwig->PrepareRun();
 174 }
 175
 176 void AliGenHerwig::InitJimmy()
 177 {
 178 // Initialisation
 179   fTarget.Resize(8);
 180   fProjectile.Resize(8);
 181   SetMC(new THerwig6());
 182   fHerwig=(THerwig6*) fMCEvGen;
 183   // initialize common blocks
 184   fHerwig->InitializeJimmy(fProjectile, fTarget, fMomentum1, fMomentum2, fProcess);
 185   // reset parameters according to user needs
 186   InitPDF();
 187   fHerwig->SetPTMIN(fPtHardMin);
 188   fHerwig->SetPTRMS(fPtRMS);
 189   fHerwig->SetMAXPR(fMaxPr);
 190   fHerwig->SetMAXER(fMaxErrors);
 191   fHerwig->SetENSOF(fEnSoft);
 192
 193 // C---D,U,S,C,B,T QUARK AND GLUON MASSES (IN THAT ORDER)
 194 //       RMASS(1)=0.32
 195 //       RMASS(2)=0.32
 196 //       RMASS(3)=0.5
 197 //       RMASS(4)=1.55
 198 //       RMASS(5)=4.75
 199 //       RMASS(6)=174.3
 200 //       RMASS(13)=0.75
 201
 202   fHerwig->SetRMASS(4,1.2);
 203   fHerwig->SetRMASS(5,4.75);
 204
 205   if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(), 49);
 206
 207   //  fHerwig->Hwusta("PI0     ");
 208
 209   // compute parameter dependent constants
 210   fHerwig->PrepareRunJimmy();
 211 }
 212
 213 void AliGenHerwig::InitPDF()
 214 {
 215 // Initialize PDF
 216   switch(fStrucFunc)
 217     {
 218 // ONLY USES LHAPDF STRUCTURE FUNCTIONS
 219     case kGRVLO98:
 220       fModPDF=80060;
 221       fAutPDF="HWLHAPDF";
 222       break;
 223     case kCTEQ6:
 224       fModPDF=10040;
 225       fAutPDF="HWLHAPDF";
 226       break;
 227     case kCTEQ61:
 228       fModPDF=10100;
 229       fAutPDF="HWLHAPDF";
 230       break;
 231     case kCTEQ6m:
 232       fModPDF=10050;
 233       fAutPDF="HWLHAPDF";
 234       break;
 235     case kCTEQ6l:
 236       fModPDF=10041;
 237       fAutPDF="HWLHAPDF";
 238       break;
 239     case kCTEQ6ll:
 240       fModPDF=10042;
 241       fAutPDF="HWLHAPDF";
 242       break;
 243     case kCTEQ5M:
 244       fModPDF=19050;
 245       fAutPDF="HWLHAPDF";
 246       break;
 247     case kCTEQ5L:
 248       fModPDF=19070;
 249       fAutPDF="HWLHAPDF";
 250       break;
 251     case kCTEQ4M:
 252       fModPDF=19150;
 253       fAutPDF="HWLHAPDF";
 254       break;
 255     case kCTEQ4L:
 256       fModPDF=19170;
 257       fAutPDF="HWLHAPDF";
 258       break;
 259 //    case kMRST2004nlo:
 260 //      fModPDF=20400;
 261 //      fAutPDF="HWLHAPDF";
 262 //      break;
 263     default:
 264       cerr << "This structure function is not inplemented " << fStrucFunc << endl;
 265       break;
 266     }
 267   fAutPDF.Resize(20);
 268   fHerwig->SetMODPDF(1,fModPDF);
 269   fHerwig->SetMODPDF(2,fModPDF);
 270   fHerwig->SetAUTPDF(1,fAutPDF);
 271   fHerwig->SetAUTPDF(2,fAutPDF);
 272 }
 273
 274 void AliGenHerwig::Generate()
 275 {
 276   // Generate one event
 277
 278   Float_t polar[3]  =   {0,0,0};
 279   Float_t origin[3] =   {0,0,0};
 280   Float_t p[4];
 281
 282   static TClonesArray *particles;
 283   //  converts from mm/c to s
 284   const Float_t kconv=0.001/2.999792458e8;
 285   //
 286   Int_t nt=0;
 287   Int_t jev=0;
 288   Int_t kf, ks, imo;
 289   kf=0;
 290
 291   if(!particles) particles=new TClonesArray("TParticle",10000);
 292
 293   fTrials=0;
 294
 295   //  Set collision vertex position
 296   if (fVertexSmear == kPerEvent) Vertex();
 297
 298   while(1)
 299     {
 300         fHerwig->GenerateEvent();
 301         fTrials++;
 302         fHerwig->ImportParticles(particles,"All");
 303         Int_t np = particles->GetEntriesFast()-1;
 304         if (np == 0 ) continue;
 305
 306         //Check hard partons or direct gamma in kine range
 307
 308         if (fProcess == kHeJets || fProcess == kHeDirectGamma) {
 309             TParticle* parton1 = (TParticle *) particles->At(6);
 310             TParticle* parton2 = (TParticle *) particles->At(7);
 311             if (!CheckParton(parton1, parton2))  continue ;
 312         }
 313
 314         //
 315         if (gAlice) {
 316             if (gAlice->GetEvNumber()>=fEv1Pr &&
 317                 gAlice->GetEvNumber()<=fEv2Pr) fHerwig->PrintEvt();
 318
 319         }
 320
 321
 322         Int_t nc = 0;
 323         fNprimaries = 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
 342                 origin[0] = fVertex[0] + iparticle->Vx()/10; // [cm]
 343                 origin[1] = fVertex[1] + iparticle->Vy()/10; // [cm]
 344                 origin[2] = fVertex[2] + iparticle->Vz()/10; // [cm]
 345
 346                 Float_t tof = fTime + kconv*iparticle->T();
 347                 Int_t   iparent = (imo > -1) ? newPos[imo] : -1;
 348                 Int_t   trackIt = (ks == 1) && fTrackIt;
 349                 PushTrack(trackIt, iparent, kf,
 350                           p[0], p[1], p[2], p[3],
 351                           origin[0], origin[1], origin[2],
 352                           tof,
 353                           polar[0], polar[1], polar[2],
 354                           kPPrimary, nt, fHerwig->GetEVWGT(), ks);
 355                 KeepTrack(nt);
 356                 newPos[i]=nt;
 357                 fNprimaries++;
 358             } // end of if: selection of particle
 359         } // end of for: particle loop
 360         if (newPos) delete[] newPos;
 361         if (nc > 0) {
 362             jev+=nc;
 363             if (jev >= fNpart || fNpart == -1) {
 364                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
 365                 break;
 366             }
 367         }
 368     }
 369 //
 370   MakeHeader();
 371 //
 372   SetHighWaterMark(nt);
 373 //  adjust weight due to kinematic selection
 374   AdjustWeights();
 375 //  get cross-section
 376   fXsection=fHerwig->GetAVWGT();
 377   //printf(">> trials << %d\n",fTrials);
 378 }
 379
 380 Bool_t AliGenHerwig::CheckParton(const TParticle* parton1, const TParticle* parton2)
 381 {
 382 // Check the kinematic trigger condition
 383 //
 384 //Select events with parton max energy
 385     if(fPtHardMax < parton1->Pt()) return kFALSE;
 386
 387 // Select events within angular window
 388     Double_t eta[2];
 389     eta[0] = parton1->Eta();
 390     eta[1] = parton2->Eta();
 391     Double_t phi[2];
 392     phi[0] = parton1->Phi();
 393     phi[1] = parton2->Phi();
 394     Int_t    pdg[2];
 395     pdg[0] = parton1->GetPdgCode();
 396     pdg[1] = parton2->GetPdgCode();
 397     printf("min %f, max %f\n",fPtHardMin, fPtHardMax);
 398     printf("Parton 1: %s, pT= %2.2f, eta = %1.2f, phi = %2.2f\n", parton1->GetName(),parton1->Pt(), eta[0],phi[0]*TMath::RadToDeg());
 399     printf("Parton 2: %s, pT= %2.2f, eta = %1.2f, phi = %2.2f\n", parton2->GetName(),parton2->Pt(), eta[1],phi[1]*TMath::RadToDeg());
 400
 401     if (fProcess == kHeJets) {
 402       //Check if one of the 2 outgoing partons are in the eta-phi window
 403       for(Int_t i = 0; i < 2; i++)
 404         if ((eta[i] < fEtaMaxParton  && eta[i] > fEtaMinParton) &&
 405             (phi[i] < fPhiMaxParton  && phi[i] > fPhiMinParton)) return  kTRUE ;
 406     }
 407
 408     else {
 409       //Check if the gamma and the jet  are in the eta-phi window
 410       Int_t igj = 0;
 411       Int_t ijj = 0;
 412       if(pdg[0] == 22) ijj=1;
 413       else igj=1;
 414       if ((eta[ijj] < fEtaMaxParton   && eta[ijj] > fEtaMinParton) &&
 415           (phi[ijj] < fPhiMaxParton   && phi[ijj] > fPhiMinParton)) {
 416
 417         if ((eta[igj] < fEtaMaxGamma   && eta[igj] > fEtaMinGamma) &&
 418             (phi[igj] < fPhiMaxGamma   && phi[igj] > fPhiMinGamma)) return  kTRUE;
 419
 420       }
 421     }
 422
 423     return kFALSE ;
 424 }
 425
 426 void AliGenHerwig::AdjustWeights()
 427 {
 428 // Adjust the weights after generation of all events
 429     TParticle *part;
 430     Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
 431     for (Int_t i=0; i<ntrack; i++) {
 432         part= gAlice->GetMCApp()->Particle(i);
 433         part->SetWeight(part->GetWeight()*fKineBias);
 434     }
 435 }
 436
 437
 438 void AliGenHerwig::KeepFullEvent()
 439 {
 440     fKeep=1;
 441 }
 442
 443 Bool_t AliGenHerwig::DaughtersSelection(const TParticle* iparticle, const TClonesArray* particles)
 444 {
 445 //
 446 // Looks recursively if one of the daughters has been selected
 447 //
 448 //    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
 449     Int_t imin=-1;
 450     Int_t imax=-1;
 451     Int_t i;
 452     Bool_t hasDaughters= (iparticle->GetFirstDaughter() >=0);
 453     Bool_t selected=kFALSE;
 454     if (hasDaughters) {
 455         imin=iparticle->GetFirstDaughter();
 456         imax=iparticle->GetLastDaughter();
 457         for (i=imin; i<= imax; i++){
 458             TParticle *  jparticle       = (TParticle *) particles->At(i);
 459             Int_t ip=jparticle->GetPdgCode();
 460             if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
 461                 selected=kTRUE; break;
 462             }
 463             if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
 464         }
 465     } else {
 466         return kFALSE;
 467     }
 468
 469     return selected;
 470 }
 471
 472
 473 Bool_t AliGenHerwig::SelectFlavor(Int_t pid) const
 474 {
 475 // Select flavor of particle
 476 // 0: all
 477 // 4: charm and beauty
 478 // 5: beauty
 479     if (fFlavor == 0) return kTRUE;
 480
 481     Int_t ifl=TMath::Abs(pid/100);
 482     if (ifl > 10) ifl/=10;
 483     return (fFlavor == ifl);
 484 }
 485
 486 Bool_t AliGenHerwig::Stable(const TParticle*  particle) const
 487 {
 488 // Return true for a stable particle
 489 //
 490     Int_t kf = TMath::Abs(particle->GetPdgCode());
 491
 492     if ( (particle->GetFirstDaughter() < 0 ) || (kf == 1000*fFlavor+122))
 493
 494     {
 495         return kTRUE;
 496     } else {
 497         return kFALSE;
 498     }
 499 }
 500
 501 void AliGenHerwig::FinishRun()
 502 {
 503   fHerwig->Hwefin();
 504 }
 505
 506 void AliGenHerwig::FinishRunJimmy()
 507 {
 508   fHerwig->Hwefin();
 509   fHerwig->Jmefin();
 510
 511 }
 512
 513
 514 void AliGenHerwig::MakeHeader()
 515 {
 516 //
 517 // Make header for the simulated event
 518 //
 519   if (fHeader) delete fHeader;
 520   fHeader = new AliGenHerwigEventHeader("Herwig");
 521 //
 522 // Event type
 523     ((AliGenHerwigEventHeader*) fHeader)->SetProcessType(fHerwig->GetIHPRO());
 524 //
 525 // Number of trials
 526     ((AliGenHerwigEventHeader*) fHeader)->SetTrials(fTrials);
 527 //
 528 // Event weight (cross section)
 529     ((AliGenHerwigEventHeader*) fHeader)->SetWeight(fHerwig->GetEVWGT());
 530 //
 531 // Event Vertex
 532     fHeader->SetPrimaryVertex(fVertex);
 533     fHeader->SetInteractionTime(fTime);
 534 //
 535 // Number of primaries
 536     fHeader->SetNProduced(fNprimaries);
 537 //  Pass header
 538 //
 539     AddHeader(fHeader);
 540     fHeader = 0x0;
 541 }