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