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("GRV"),
  41     fModPDF(5),
  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("GRV"),
  71     fModPDF(5),
  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("GRV"),
 104     fModPDF(5),
 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.95
 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::InitPDF()
 186 {
 187   switch(fStrucFunc)
 188     {
 189 //    case kGRVLO:
 190 //      fModPDF=5;
 191 //      fAutPDF="GRV";
 192 //      break;
 193 //    case kGRVHO:
 194 //      fModPDF=6;
 195 //      fAutPDF="GRV";
 196 //      break;
 197     case kGRVLO98:
 198       fModPDF=12;
 199       fAutPDF="GRV";
 200       break;
 201 //    case kMRSDminus:
 202 //      fModPDF=31;
 203 //      fAutPDF="MRS";
 204 //      break;
 205 //    case kMRSD0:
 206 //      fModPDF=30;
 207 //      fAutPDF="MRS";
 208 //      break;
 209 //    case kMRSG:
 210 //      fModPDF=41;
 211 //      fAutPDF="MRS";
 212 //      break;
 213 //    case kMRSTcgLO:
 214 //      fModPDF=72;
 215 //      fAutPDF="MRS";
 216 //      break;
 217     case kCTEQ4M:
 218       fModPDF=34;
 219       fAutPDF="CTEQ";
 220       break;
 221     case kCTEQ5L:
 222       fModPDF=46;
 223       fAutPDF="CTEQ";
 224       break;
 225 //    case kCTEQ5M:
 226 //      fModPDF=48;
 227 //      fAutPDF="CTEQ";
 228 //      break;
 229     default:
 230       cerr << "This structure function is not inplemented " << fStrucFunc << endl;
 231       break;
 232     }
 233   fAutPDF.Resize(20);
 234   fHerwig->SetMODPDF(1,fModPDF);
 235   fHerwig->SetMODPDF(2,fModPDF);
 236   fHerwig->SetAUTPDF(1,fAutPDF);
 237   fHerwig->SetAUTPDF(2,fAutPDF);
 238 }
 239
 240 void AliGenHerwig::Generate()
 241 {
 242   // Generate one event
 243
 244   Float_t polar[3] =   {0,0,0};
 245   Float_t origin[3]=   {0,0,0};
 246   Float_t origin0[3]=  {0,0,0};
 247   Float_t p[4], random[6];
 248
 249   static TClonesArray *particles;
 250   //  converts from mm/c to s
 251   const Float_t kconv=0.001/2.999792458e8;
 252   //
 253   Int_t nt=0;
 254   Int_t jev=0;
 255   Int_t j, kf, ks, imo;
 256   kf=0;
 257
 258   if(!particles) particles=new TClonesArray("TParticle",10000);
 259
 260   fTrials=0;
 261   for (j=0;j<3;j++) origin0[j]=fOrigin[j];
 262   if(fVertexSmear==kPerEvent) {
 263     Rndm(random,6);
 264     for (j=0;j<3;j++) {
 265       origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
 266         TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 267     }
 268   }
 269
 270   while(1)
 271     {
 272         fHerwig->GenerateEvent();
 273         fTrials++;
 274         fHerwig->ImportParticles(particles,"All");
 275         Int_t np = particles->GetEntriesFast()-1;
 276         if (np == 0 ) continue;
 277
 278         Int_t nc=0;
 279
 280         Int_t * newPos = new Int_t[np];
 281         for (Int_t i = 0; i<np; i++) *(newPos+i)=-1;
 282
 283         for (Int_t i = 0; i<np; i++) {
 284             TParticle *  iparticle       = (TParticle *) particles->At(i);
 285             imo = iparticle->GetFirstMother();
 286             kf        = iparticle->GetPdgCode();
 287             ks        = iparticle->GetStatusCode();
 288             if (ks != 3 &&
 289                 KinematicSelection(iparticle,0))
 290             {
 291                 nc++;
 292                 p[0]=iparticle->Px();
 293                 p[1]=iparticle->Py();
 294                 p[2]=iparticle->Pz();
 295                 p[3]=iparticle->Energy();
 296                 origin[0]=origin0[0]+iparticle->Vx()/10;
 297                 origin[1]=origin0[1]+iparticle->Vy()/10;
 298                 origin[2]=origin0[2]+iparticle->Vz()/10;
 299                 Float_t tof = kconv*iparticle->T();
 300                 Int_t   iparent = (imo > -1) ? newPos[imo] : -1;
 301                 Int_t   trackIt = (ks == 1) && fTrackIt;
 302                 PushTrack(trackIt, iparent, kf,
 303                           p[0], p[1], p[2], p[3],
 304                           origin[0], origin[1], origin[2],
 305                           tof,
 306                           polar[0], polar[1], polar[2],
 307                           kPPrimary, nt, fHerwig->GetEVWGT(), ks);
 308                 KeepTrack(nt);
 309                 newPos[i]=nt;
 310             } // end of if: selection of particle
 311         } // end of for: particle loop
 312         if (newPos) delete[] newPos;
 313         //      MakeHeader();
 314         if (nc > 0) {
 315             jev+=nc;
 316             if (jev >= fNpart || fNpart == -1) {
 317                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
 318                 break;
 319             }
 320         }
 321     }
 322   SetHighWaterMark(nt);
 323 //  adjust weight due to kinematic selection
 324   AdjustWeights();
 325 //  get cross-section
 326   fXsection=fHerwig->GetAVWGT();
 327 }
 328
 329 void AliGenHerwig::AdjustWeights()
 330 {
 331 // Adjust the weights after generation of all events
 332     TParticle *part;
 333     Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
 334     for (Int_t i=0; i<ntrack; i++) {
 335         part= gAlice->GetMCApp()->Particle(i);
 336         part->SetWeight(part->GetWeight()*fKineBias);
 337     }
 338 }
 339
 340
 341 void AliGenHerwig::KeepFullEvent()
 342 {
 343     fKeep=1;
 344 }
 345
 346 Bool_t AliGenHerwig::DaughtersSelection(TParticle* iparticle, TClonesArray* particles)
 347 {
 348 //
 349 // Looks recursively if one of the daughters has been selected
 350 //
 351 //    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
 352     Int_t imin=-1;
 353     Int_t imax=-1;
 354     Int_t i;
 355     Bool_t hasDaughters= (iparticle->GetFirstDaughter() >=0);
 356     Bool_t selected=kFALSE;
 357     if (hasDaughters) {
 358         imin=iparticle->GetFirstDaughter();
 359         imax=iparticle->GetLastDaughter();
 360         for (i=imin; i<= imax; i++){
 361             TParticle *  jparticle       = (TParticle *) particles->At(i);
 362             Int_t ip=jparticle->GetPdgCode();
 363             if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
 364                 selected=kTRUE; break;
 365             }
 366             if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
 367         }
 368     } else {
 369         return kFALSE;
 370     }
 371
 372     return selected;
 373 }
 374
 375
 376 Bool_t AliGenHerwig::SelectFlavor(Int_t pid)
 377 {
 378 // Select flavor of particle
 379 // 0: all
 380 // 4: charm and beauty
 381 // 5: beauty
 382     if (fFlavor == 0) return kTRUE;
 383
 384     Int_t ifl=TMath::Abs(pid/100);
 385     if (ifl > 10) ifl/=10;
 386     return (fFlavor == ifl);
 387 }
 388
 389 Bool_t AliGenHerwig::Stable(TParticle*  particle)
 390 {
 391 // Return true for a stable particle
 392 //
 393     Int_t kf = TMath::Abs(particle->GetPdgCode());
 394
 395     if ( (particle->GetFirstDaughter() < 0 ) || (kf == 1000*fFlavor+122))
 396
 397     {
 398         return kTRUE;
 399     } else {
 400         return kFALSE;
 401     }
 402 }
 403
 404 void AliGenHerwig::FinishRun()
 405 {
 406   fHerwig->Hwefin();
 407 }
 408
 409
 410 AliGenHerwig& AliGenHerwig::operator=(const  AliGenHerwig& rhs)
 411 {
 412 // Assignment operator
 413     rhs.Copy(*this);
 414     return (*this);
 415 }
 416
 417