THijing/THijing.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 //                                                                            //
  20 // THijing                                                                    //
  21 //                                                                            //
  22 // THijing is an interface class to F77 version of Hijing 1.36                //
  23 // event generator, written by X.N. Wang and  M. Gyulassy.                    //
  24 // For details see http://nta2.lbl.gov/~xnwang                                //
  25 //                                                                            //
  26 //          **************************************************                //
  27 //          *     |             _______      /  ------/      *                //
  28 //          *   ----- ------     |_____|     /_/     /       *                //
  29 //          *    ||    /        |_____|      /    /          *                //
  30 //          *    /|   /_/       /_______    /_  /    _       *                //
  31 //          *   / |     / /     /  /  / |        -------     *                //
  32 //          *     |    / /       /  /  |     /     |         *                //
  33 //          *     |   / /       /  / _|    /   -------       *                //
  34 //          *                                                *                //
  35 //          **************************************************                //
  36 //                                HIJING                                      //
  37 //                 Heavy Ion Jet INteraction Generator                        //
  38 //                                  by                                        //
  39 //                     X. N. Wang  and  M. Gyulassy                           //
  40 //                     Lawrence Berkeley Laboratory                           //
  41 //****************************************************************************//
  42
  43
  44 #include <TClonesArray.h>
  45 #include <TObjArray.h>
  46 #include <TParticle.h>
  47 #include <TROOT.h>
  48
  49 #include "Hcommon.h"
  50 #include "THijing.h"
  51
  52 #ifndef WIN32
  53 # define hijset hijset_
  54 # define hijing hijing_
  55 # define profile profile_
  56 # define rluget_hijing rluget_hijing_
  57 # define rluset_hijing rluset_hijing_
  58 # define type_of_call
  59 #else
  60 # define hijset HIJSET
  61 # define hijing HIJING
  62 # define profile PROFILE
  63 # define rluget_hijing RLUGET_HIJING
  64 # define rluset_hijing RLUSET_HIJING
  65 # define type_of_call _stdcall
  66 #endif
  67
  68 #ifndef WIN32
  69 extern "C" void type_of_call hijset(Float_t & , const char *,
  70                                    const char *, const char *,
  71                                    Int_t & , Int_t &, Int_t &,
  72                                    Int_t &, const int,
  73                                    const int, const int);
  74 extern "C" float type_of_call profile(Float_t &);
  75 extern "C" void type_of_call hijing(const char *, Float_t  &,
  76                                    Float_t &, const int);
  77 extern "C" void type_of_call rluget_hijing(Int_t & lfn, Int_t & move);
  78
  79 extern "C" void type_of_call rluset_hijing(Int_t & lfn, Int_t & move);
  80
  81 #else
  82 #endif
  83
  84
  85
  86 ClassImp(THijing)
  87
  88
  89 THijing::THijing():
  90     TGenerator("Hijing","Hijing"),
  91     fEfrm(5500.),
  92     fFrame("CMS"),
  93     fProj("A"),
  94     fTarg("A"),
  95     fIap(208),
  96     fIzp(82),
  97     fIat(208),
  98     fIzt(82),
  99     fBmin(0.),
 100     fBmax(5.)
 101 {
 102 // Default constructor
 103 }
 104
 105 //______________________________________________________________________________
 106 THijing::THijing(Float_t efrm, const char *frame="CMS",
 107                  const char *proj="A", const char *targ="A",
 108                  Int_t iap=207, Int_t izp=82, Int_t iat=207, Int_t izt=82,
 109                  Float_t bmin=0, Float_t bmax=20):
 110     TGenerator("Hijing","Hijing"),
 111     fEfrm(efrm),
 112     fFrame(frame),
 113     fProj(proj),
 114     fTarg(targ),
 115     fIap(iap),
 116     fIzp(izp),
 117     fIat(iat),
 118     fIzt(izt),
 119     fBmin(bmin),
 120     fBmax(bmax)
 121 {
 122 // THijing constructor:
 123 // Note that there may be only one functional THijing object
 124 // at a time, so it's not use to create more than one instance of it.
 125 }
 126
 127 //______________________________________________________________________________
 128 THijing::~THijing()
 129 {
 130 // Destructor
 131 }
 132
 133
 134 TObjArray* THijing::ImportParticles(Option_t *option)
 135 {
 136 //
 137 //  Default primary creation method. It reads the /HEPEVT/ common block which
 138 //  has been filled by the GenerateEvent method. If the event generator does
 139 //  not use the HEPEVT common block, This routine has to be overloaded by
 140 //  the subclasses.
 141 //  The function loops on the generated particles and store them in
 142 //  the TClonesArray pointed by the argument particles.
 143 //  The default action is to store only the stable particles (ISTHEP = 1)
 144 //  This can be demanded explicitly by setting the option = "Final"
 145 //  If the option = "All", all the particles are stored.
 146 //
 147     fParticles->Clear();
 148     Int_t numpart = HIMAIN1.natt;
 149     printf("\n THijing: HIJING stack contains %d particles.", numpart);
 150     printf("\n THijing: Total energy:         %f           ", HIMAIN1.eatt);
 151     printf("\n THijing: Number of hard scatterings: %d     ", HIMAIN1.jatt);
 152     Int_t nump = 0;
 153     if (!strcmp(option,"") || !strcmp(option,"Final")) {
 154         for (Int_t i = 0; i < numpart; i++) {
 155
 156             if (HIMAIN2.katt[3][i] == 1) {
 157 //
 158 //  Use the common block values for the TParticle constructor
 159 //
 160                 nump++;
 161                 TParticle* p = new TParticle(
 162                     HIMAIN2.katt[0][i], HIMAIN2.katt[3][i] ,
 163                     -1, -1, -1, -1,
 164                     HIMAIN2.patt[0][i], HIMAIN2.patt[1][i], HIMAIN2.patt[2][i], HIMAIN2.patt[3][i] ,
 165                     HIMAIN2.vatt[0][i], HIMAIN2.vatt[1][i], HIMAIN2.vatt[2][i], HIMAIN2.vatt[3][i]
 166                   );
 167                 p->SetUniqueID(HIMAIN2.katt[1][i]);
 168                 fParticles->Add(p);
 169             }
 170         }
 171     }
 172     else if (!strcmp(option,"All")) {
 173         nump = numpart;
 174         for (Int_t i = 0; i < numpart; i++) {
 175
 176             Int_t iParent = HIMAIN2.katt[2][i]-1;
 177
 178             if (iParent >= 0) {
 179                 TParticle *mother = (TParticle*) (fParticles->UncheckedAt(iParent));
 180                 mother->SetLastDaughter(i);
 181                 if (mother->GetFirstDaughter()==-1)
 182                     mother->SetFirstDaughter(i);
 183             }
 184
 185             TParticle* p = new TParticle(
 186                 HIMAIN2.katt[0][i], HIMAIN2.katt[3][i], iParent,
 187                 -1, -1, -1,
 188                 HIMAIN2.patt[0][i], HIMAIN2.patt[1][i], HIMAIN2.patt[2][i], HIMAIN2.patt[3][i] ,
 189                 HIMAIN2.vatt[0][i], HIMAIN2.vatt[1][i], HIMAIN2.vatt[2][i], HIMAIN2.vatt[3][i]
 190                 );
 191             p->SetUniqueID(HIMAIN2.katt[1][i]);
 192             fParticles->Add(p);
 193         }
 194     }
 195     return fParticles;
 196 }
 197
 198 Int_t THijing::ImportParticles(TClonesArray *particles, Option_t *option)
 199 {
 200 //
 201 //  Default primary creation method. It reads the /HEPEVT/ common block which
 202 //  has been filled by the GenerateEvent method. If the event generator does
 203 //  not use the HEPEVT common block, This routine has to be overloaded by
 204 //  the subclasses.
 205 //  The function loops on the generated particles and store them in
 206 //  the TClonesArray pointed by the argument particles.
 207 //  The default action is to store only the stable particles (ISTHEP = 1)
 208 //  This can be demanded explicitly by setting the option = "Final"
 209 //  If the option = "All", all the particles are stored.
 210 //
 211   if (particles == 0) return 0;
 212   TClonesArray &particlesR = *particles;
 213   particlesR.Clear();
 214   Int_t numpart = HIMAIN1.natt;
 215   printf("\n THijing: HIJING stack contains %d particles.", numpart);
 216   printf("\n THijing: Total energy:         %f           ", HIMAIN1.eatt);
 217   printf("\n THijing: Number of hard scatterings: %d     ", HIMAIN1.jatt);
 218   Int_t nump = 0;
 219   if (!strcmp(option,"") || !strcmp(option,"Final")) {
 220       for (Int_t i = 0; i < numpart; i++) {
 221
 222           if (HIMAIN2.katt[3][i] == 1) {
 223 //
 224 //  Use the common block values for the TParticle constructor
 225 //
 226             nump++;
 227             new(particlesR[i]) TParticle(
 228                   HIMAIN2.katt[0][i] ,
 229                   HIMAIN2.katt[3][i] ,
 230                   -1 ,
 231                   -1,
 232                   -1,
 233                   -1,
 234
 235                   HIMAIN2.patt[0][i] ,
 236                   HIMAIN2.patt[1][i] ,
 237                   HIMAIN2.patt[2][i] ,
 238                   HIMAIN2.patt[3][i] ,
 239
 240                   HIMAIN2.vatt[0][i] ,
 241                   HIMAIN2.vatt[1][i] ,
 242                   HIMAIN2.vatt[2][i] ,
 243                   HIMAIN2.vatt[3][i]
 244                   );
 245             particlesR[i]->SetUniqueID(HIMAIN2.katt[1][i]);
 246           }
 247       }
 248   }
 249   else if (!strcmp(option,"All")) {
 250       nump = numpart;
 251       for (Int_t i = 0; i < numpart; i++) {
 252
 253           Int_t iParent = HIMAIN2.katt[2][i]-1;
 254
 255           if (iParent >= 0) {
 256               TParticle *mother = (TParticle*) (particlesR.UncheckedAt(iParent));
 257               mother->SetLastDaughter(i);
 258               if (mother->GetFirstDaughter()==-1)
 259                   mother->SetFirstDaughter(i);
 260           }
 261
 262           new(particlesR[i]) TParticle(
 263               HIMAIN2.katt[0][i] ,
 264               HIMAIN2.katt[3][i] ,
 265               iParent,
 266               -1,
 267               -1,
 268               -1,
 269
 270               HIMAIN2.patt[0][i] ,
 271               HIMAIN2.patt[1][i] ,
 272               HIMAIN2.patt[2][i] ,
 273               HIMAIN2.patt[3][i] ,
 274
 275               HIMAIN2.vatt[0][i] ,
 276               HIMAIN2.vatt[1][i] ,
 277               HIMAIN2.vatt[2][i] ,
 278               HIMAIN2.vatt[3][i]
 279               );
 280           particlesR[i]->SetUniqueID(HIMAIN2.katt[1][i]);
 281       }
 282   }
 283   return nump;
 284 }
 285
 286 //______________________________________________________________________________
 287 void THijing::SetEFRM(Float_t efrm)
 288 {
 289 // Set the centre of mass (CMS) or lab-energy (LAB)
 290    fEfrm=efrm;
 291 }
 292 //______________________________________________________________________________
 293 void THijing::SetFRAME(const char* frame)
 294 {
 295 // Set the frame type ("CMS" or "LAB")
 296    fFrame=frame;
 297 }
 298 //______________________________________________________________________________
 299 void THijing::SetPROJ(const char* proj)
 300 {
 301 // Set the projectile type
 302    fProj=proj;
 303 }
 304 //______________________________________________________________________________
 305 void THijing::SetTARG(const char* targ)
 306 {
 307 // Set the target type
 308    fTarg=targ;
 309 }
 310 //______________________________________________________________________________
 311 void THijing::SetIAP(Int_t iap)
 312 {
 313 // Set the projectile atomic number
 314    fIap=iap;
 315 }
 316 //______________________________________________________________________________
 317 void THijing::SetIZP(Int_t izp)
 318 {
 319 // Set the projectile charge number
 320    fIzp=izp;
 321 }
 322 //______________________________________________________________________________
 323 void THijing::SetIAT(Int_t iat)
 324 {
 325 // Set the target atomic number
 326    fIat=iat;
 327 }
 328 //______________________________________________________________________________
 329 void THijing::SetIZT(Int_t izt)
 330 {
 331 // Set the target charge number
 332    fIzt=izt;
 333 }
 334 //______________________________________________________________________________
 335 void THijing::SetBMIN(Float_t bmin)
 336 {
 337 // Set the minimum impact parameter
 338    fBmin=bmin;
 339 }
 340 //______________________________________________________________________________
 341 void THijing::SetBMAX(Float_t bmax)
 342 {
 343 // Set the maximum impact parameter
 344    fBmax=bmax;
 345 }
 346 //______________________________________________________________________________
 347 Float_t THijing::GetEFRM() const
 348 {
 349 // Get the centre of mass (CMS) or lab-energy (LAB)
 350    return fEfrm;
 351 }
 352 //______________________________________________________________________________
 353 const char* THijing::GetFRAME() const
 354 {
 355 // Get the frame type ("CMS" or "LAB")
 356    return fFrame.Data();
 357 }
 358 //______________________________________________________________________________
 359 const char* THijing::GetPROJ() const
 360 {
 361 // Get the projectile type
 362    return fProj;
 363 }
 364 //______________________________________________________________________________
 365 const char* THijing::GetTARG() const
 366 {
 367 // Set the target type
 368    return fTarg;
 369 }
 370 //______________________________________________________________________________
 371 Int_t THijing::GetIAP() const
 372 {
 373 // Get the projectile atomic number
 374    return fIap;
 375 }
 376 //______________________________________________________________________________
 377 Int_t THijing::GetIZP() const
 378 {
 379 // Get the projectile charge number
 380    return fIzp;
 381 }
 382 //______________________________________________________________________________
 383 Int_t THijing::GetIAT() const
 384 {
 385 // Get the target atomic number
 386    return fIat;
 387 }
 388 //______________________________________________________________________________
 389 Int_t THijing::GetIZT() const
 390 {
 391 // Get the target charge number
 392    return fIzt;
 393 }
 394 //______________________________________________________________________________
 395 Float_t THijing::GetBMIN() const
 396 {
 397 // Get the minimum impact parameter
 398    return fBmin;
 399 }
 400 //______________________________________________________________________________
 401 Float_t THijing::GetBMAX() const
 402 {
 403 // Get the maximum impact parameter
 404    return fBmax;
 405 }
 406
 407 //====================== access to common HIPARNT ===============================
 408
 409 //______________________________________________________________________________
 410 void THijing::SetHIPR1(Int_t key,Float_t value)
 411 {
 412 // Set the values of array HIPR1 in common HIPARNT
 413    if ( key<1 || key>100 ) {
 414       printf ("ERROR in THijing:SetHIPR1(key,value): \n ");
 415       printf ("      key=%i is out of range [1..100]!\n",key);
 416       return;
 417    }
 418
 419    HIPARNT.hipr1[key-1]=value;
 420
 421 }
 422
 423 //______________________________________________________________________________
 424 Float_t THijing::GetHIPR1(Int_t key) const
 425 {
 426 // Get the values of array HIPR1 in common HIPARNT
 427    if ( key<1 || key>100 ) {
 428       printf ("ERROR in THijing:GetHIPR1(key): \n ");
 429       printf ("      key=%i is out of range [1..100]!\n",key);
 430       return 0;
 431    }
 432
 433    return HIPARNT.hipr1[key-1];
 434
 435 }
 436
 437 //______________________________________________________________________________
 438 void THijing::SetIHPR2(Int_t key,Int_t value)
 439 {
 440 // Set the values of array HIPR2 in common HIPARNT
 441    if ( key<1 || key>50 ) {
 442       printf ("ERROR in THijing:SetIHPR2(key,value): \n ");
 443       printf ("      key=%i is out of range [1..50]!\n",key);
 444       return;
 445    }
 446
 447    HIPARNT.ihpr2[key-1]=value;
 448
 449 }
 450
 451 //______________________________________________________________________________
 452 Int_t THijing::GetIHPR2(Int_t key) const
 453 {
 454 // Get the values of array HIPR2 in common HIPARNT
 455    if ( key<1 || key>50 ) {
 456       printf ("ERROR in THijing:GetIHPR2(key): \n ");
 457       printf ("      key=%i is out of range [1..50]!\n",key);
 458       return 0;
 459    }
 460
 461    return HIPARNT.ihpr2[key-1];
 462
 463 }
 464
 465
 466 //______________________________________________________________________________
 467 Float_t THijing::GetHINT1(Int_t key) const
 468 {
 469 // Get the values of array HINT1 in common HIPARNT
 470    if ( key<1 || key>100 ) {
 471       printf ("ERROR in THijing:GetHINT1(key): \n ");
 472       printf ("      key=%i is out of range [1..100]!\n",key);
 473       return 0;
 474    }
 475
 476    return HIPARNT.hint1[key-1];
 477
 478 }
 479
 480
 481 //______________________________________________________________________________
 482 Int_t THijing::GetIHNT2(Int_t key) const
 483 {
 484 // Get the values of array HINT2 in common HIPARNT
 485    if ( key<1 || key>50 ) {
 486       printf ("ERROR in THijing:GetIHNT2(key): \n ");
 487       printf ("      key=%i is out of range [1..50]!\n",key);
 488       return 0;
 489    }
 490
 491    return HIPARNT.ihnt2[key-1];
 492
 493 }
 494
 495
 496 //====================== access to common HIMAIN1 ===============================
 497
 498 //______________________________________________________________________________
 499 Int_t  THijing::GetNATT() const
 500 {
 501 // Get the number of particles produces
 502    return HIMAIN1.natt;
 503
 504 }
 505
 506 //______________________________________________________________________________
 507 Float_t  THijing::GetEATT() const
 508 {
 509 // Get total energy of particles
 510
 511    return HIMAIN1.eatt;
 512
 513 }
 514
 515 //______________________________________________________________________________
 516 Int_t  THijing::GetJATT() const
 517 {
 518 // Get number of hard scatterings
 519
 520    return HIMAIN1.jatt;
 521
 522 }
 523
 524 //______________________________________________________________________________
 525 Int_t  THijing::GetNT() const
 526 {
 527 // Get number of target participants
 528
 529    return HIMAIN1.nt;
 530
 531 }
 532
 533 //______________________________________________________________________________
 534 Int_t  THijing::GetNP() const
 535 {
 536 // Get number of projectile participants
 537    return HIMAIN1.np;
 538
 539 }
 540
 541
 542 //______________________________________________________________________________
 543 Int_t  THijing::GetN0() const
 544 {
 545 // Get number of N-N collisions
 546    return HIMAIN1.n0;
 547
 548 }
 549 //______________________________________________________________________________
 550 Int_t  THijing::GetN01() const
 551 {
 552 // Get number of N-wounded collisions
 553
 554    return HIMAIN1.n01;
 555
 556 }
 557
 558 //______________________________________________________________________________
 559 Int_t  THijing::GetN10() const
 560 {
 561 // Get number of wounded-N collisions
 562
 563    return HIMAIN1.n10;
 564
 565 }
 566
 567 //______________________________________________________________________________
 568 Int_t  THijing::GetN11() const
 569 {
 570 // Get number of wounded-wounded collisions
 571
 572    return HIMAIN1.n11;
 573
 574 }
 575
 576 //______________________________________________________________________________
 577 Float_t  THijing::GetBB() const
 578 {
 579 // Get impact parameter
 580
 581    return HIMAIN1.bb;
 582
 583 }
 584
 585 //====================== access to common HIMAIN2 ===============================
 586
 587 //______________________________________________________________________________
 588 Int_t THijing::GetKATT(Int_t key1, Int_t key2) const
 589 {
 590 // Get values of array KATT in common HIMAIN2
 591    if ( key1<1 || key1>200000 ) {
 592       printf("ERROR in THijing::GetKATT(key1,key2):\n");
 593       printf("      key1=%i is out of range [1..200000]\n",key1);
 594       return 0;
 595    }
 596
 597    if ( key2<1 || key2>4 ) {
 598       printf("ERROR in THijing::GetKATT(key1,key2):\n");
 599       printf("      key2=%i is out of range [1..4]\n",key2);
 600       return 0;
 601    }
 602
 603    return   HIMAIN2.katt[key2-1][key1-1];
 604 }
 605
 606 //______________________________________________________________________________
 607 Float_t THijing::GetPATT(Int_t key1, Int_t key2) const
 608 {
 609 // Get values of array PATT in common HIMAIN2
 610    if ( key1<1 || key1>200000 ) {
 611       printf("ERROR in THijing::GetPATT(key1,key2):\n");
 612       printf("      key1=%i is out of range [1..130000]\n",key1);
 613       return 0;
 614    }
 615
 616    if ( key2<1 || key2>4 ) {
 617       printf("ERROR in THijing::GetPATT(key1,key2):\n");
 618       printf("      key2=%i is out of range [1..4]\n",key2);
 619       return 0;
 620    }
 621
 622    return   HIMAIN2.patt[key2-1][key1-1];
 623 }
 624
 625 Float_t THijing::GetVATT(Int_t key1, Int_t key2) const
 626 {
 627 // Get values of array VATT in common HIMAIN2
 628    if ( key1<1 || key1>200000 ) {
 629       printf("ERROR in THijing::GetVATT(key1,key2):\n");
 630       printf("      key1=%i is out of range [1..130000]\n",key1);
 631       return 0;
 632    }
 633
 634    if ( key2<1 || key2>4 ) {
 635       printf("ERROR in THijing::GetVATT(key1,key2):\n");
 636       printf("      key2=%i is out of range [1..4]\n",key2);
 637       return 0;
 638    }
 639
 640    return   HIMAIN2.vatt[key2-1][key1-1];
 641 }
 642
 643 //====================== access to common HIJJET1 ===============================
 644
 645 //______________________________________________________________________________
 646 Int_t THijing::GetNPJ(Int_t key) const
 647 {
 648 // Get values of array NPJ of common HIJJET1
 649    if ( key<1 || key>300 ) {
 650       printf("ERROR in THijing::GetNPJ(key):\n");
 651       printf("      key=%i is out of range [1..300]\n",key);
 652       return 0;
 653    }
 654    return HIJJET1.npj[key-1];
 655 }
 656
 657 //______________________________________________________________________________
 658 Int_t THijing::GetKFPJ(Int_t key1, Int_t key2) const
 659 {
 660 // Get values of array KFPJ in common HIJJET1
 661    if ( key1<1 || key1>300 ) {
 662       printf("ERROR in THijing::GetKFPJ(key1):\n");
 663       printf("      key1=%i is out of range [1..300]\n",key1);
 664       return 0;
 665    }
 666    if ( key2<1 || key2>500 ) {
 667       printf("ERROR in THijing::GetKFPJ(key1,key2):\n");
 668       printf("      key2=%i is out of range [1..500]\n",key2);
 669       return 0;
 670    }
 671
 672    return HIJJET1.kfpj[key2-1][key1-1];
 673 }
 674
 675 //______________________________________________________________________________
 676 Float_t THijing::GetPJPX(Int_t key1, Int_t key2) const
 677 {
 678 // Get values of array PJPX in common HIJJET1
 679    if ( key1<1 || key1>300 ) {
 680       printf("ERROR in THijing::GetPJPX(key1):\n");
 681       printf("      key1=%i is out of range [1..300]\n",key1);
 682       return 0;
 683    }
 684    if ( key2<1 || key2>500 ) {
 685       printf("ERROR in THijing::GetPJPX(key1,key2):\n");
 686       printf("      key2=%i is out of range [1..500]\n",key2);
 687       return 0;
 688    }
 689
 690    return HIJJET1.pjpx[key2-1][key1-1];
 691 }
 692
 693 //______________________________________________________________________________
 694 Float_t THijing::GetPJPY(Int_t key1, Int_t key2) const
 695 {
 696 // Get values of array PJPY in common HIJJET1
 697    if ( key1<1 || key1>300 ) {
 698       printf("ERROR in THijing::GetPJPY(key1):\n");
 699       printf("      key1=%i is out of range [1..300]\n",key1);
 700       return 0;
 701    }
 702    if ( key2<1 || key2>500 ) {
 703       printf("ERROR in THijing::GetPJPY(key1,key2):\n");
 704       printf("      key2=%i is out of range [1..500]\n",key2);
 705       return 0;
 706    }
 707
 708    return HIJJET1.pjpy[key2-1][key1-1];
 709 }
 710
 711 //______________________________________________________________________________
 712 Float_t THijing::GetPJPZ(Int_t key1, Int_t key2) const
 713 {
 714 // Get values of array PJPZ in common HIJJET1
 715    if ( key1<1 || key1>300 ) {
 716       printf("ERROR in THijing::GetPJPZ(key1):\n");
 717       printf("      key1=%i is out of range [1..300]\n",key1);
 718       return 0;
 719    }
 720    if ( key2<1 || key2>500 ) {
 721       printf("ERROR in THijing::GetPJPZ(key1,key2):\n");
 722       printf("      key2=%i is out of range [1..500]\n",key2);
 723       return 0;
 724    }
 725
 726    return HIJJET1.pjpz[key2-1][key1-1];
 727 }
 728
 729 //______________________________________________________________________________
 730 Float_t THijing::GetPJPE(Int_t key1, Int_t key2) const
 731 {
 732 // Get values of array PJPE in common HIJJET1
 733    if ( key1<1 || key1>300 ) {
 734       printf("ERROR in THijing::GetPJPE(key1):\n");
 735       printf("      key1=%i is out of range [1..300]\n",key1);
 736       return 0;
 737    }
 738    if ( key2<1 || key2>500 ) {
 739       printf("ERROR in THijing::GetPJPE(key1,key2):\n");
 740       printf("      key2=%i is out of range [1..500]\n",key2);
 741       return 0;
 742    }
 743
 744    return HIJJET1.pjpe[key2-1][key1-1];
 745 }
 746
 747 //______________________________________________________________________________
 748 Float_t THijing::GetPJPM(Int_t key1, Int_t key2) const
 749 {
 750 // Get values of array PJPM in common HIJJET1
 751    if ( key1<1 || key1>300 ) {
 752       printf("ERROR in THijing::GetPJPM(key1):\n");
 753       printf("      key1=%i is out of range [1..300]\n",key1);
 754       return 0;
 755    }
 756    if ( key2<1 || key2>500 ) {
 757       printf("ERROR in THijing::GetPJPM(key1,key2):\n");
 758       printf("      key2=%i is out of range [1..500]\n",key2);
 759       return 0;
 760    }
 761
 762    return HIJJET1.pjpm[key2-1][key1-1];
 763 }
 764
 765 //______________________________________________________________________________
 766 Int_t THijing::GetNTJ(Int_t key) const
 767 {
 768 // Get values of array NTJ in common HIJJET1
 769    if ( key<1 || key>300 ) {
 770       printf("ERROR in THijing::GetNTJ(key):\n");
 771       printf("      key=%i is out of range [1..300]\n",key);
 772       return 0;
 773    }
 774    return HIJJET1.ntj[key-1];
 775 }
 776
 777 //______________________________________________________________________________
 778 Int_t THijing::GetKFTJ(Int_t key1, Int_t key2) const
 779 {
 780 // Get values of array KFTJ in common HIJJET1
 781    if ( key1<1 || key1>300 ) {
 782       printf("ERROR in THijing::GetKFTJ(key1):\n");
 783       printf("      key1=%i is out of range [1..300]\n",key1);
 784       return 0;
 785    }
 786    if ( key2<1 || key2>500 ) {
 787       printf("ERROR in THijing::GetKFTJ(key1,key2):\n");
 788       printf("      key2=%i is out of range [1..500]\n",key2);
 789       return 0;
 790    }
 791
 792    return HIJJET1.kftj[key2-1][key1-1];
 793 }
 794
 795 //______________________________________________________________________________
 796 Float_t THijing::GetPJTX(Int_t key1, Int_t key2) const
 797 {
 798 // Get values of array PJTX in common HIJJET1
 799    if ( key1<1 || key1>300 ) {
 800       printf("ERROR in THijing::GetPJTX(key1):\n");
 801       printf("      key1=%i is out of range [1..300]\n",key1);
 802       return 0;
 803    }
 804    if ( key2<1 || key2>500 ) {
 805       printf("ERROR in THijing::GetPJTX(key1,key2):\n");
 806       printf("      key2=%i is out of range [1..500]\n",key2);
 807       return 0;
 808    }
 809
 810    return HIJJET1.pjtx[key2-1][key1-1];
 811 }
 812
 813 //______________________________________________________________________________
 814 Float_t THijing::GetPJTY(Int_t key1, Int_t key2) const
 815 {
 816 // Get values of array PJTY in common HIJJET1
 817    if ( key1<1 || key1>300 ) {
 818       printf("ERROR in THijing::GetPJTY(key1):\n");
 819       printf("      key1=%i is out of range [1..300]\n",key1);
 820       return 0;
 821    }
 822    if ( key2<1 || key2>500 ) {
 823       printf("ERROR in THijing::GetPJTY(key1,key2):\n");
 824       printf("      key2=%i is out of range [1..500]\n",key2);
 825       return 0;
 826    }
 827
 828    return HIJJET1.pjty[key2-1][key1-1];
 829 }
 830
 831 //______________________________________________________________________________
 832 Float_t THijing::GetPJTZ(Int_t key1, Int_t key2) const
 833 {
 834 // Get values of array PJTZ in common HIJJET1
 835    if ( key1<1 || key1>300 ) {
 836       printf("ERROR in THijing::GetPJTZ(key1):\n");
 837       printf("      key1=%i is out of range [1..300]\n",key1);
 838       return 0;
 839    }
 840    if ( key2<1 || key2>500 ) {
 841       printf("ERROR in THijing::GetPJTZ(key1,key2):\n");
 842       printf("      key2=%i is out of range [1..500]\n",key2);
 843       return 0;
 844    }
 845
 846    return HIJJET1.pjtz[key2-1][key1-1];
 847 }
 848
 849 //______________________________________________________________________________
 850 Float_t THijing::GetPJTE(Int_t key1, Int_t key2) const
 851 {
 852 // Get values of array PJTE in common HIJJET1
 853    if ( key1<1 || key1>300 ) {
 854       printf("ERROR in THijing::GetPJTE(key1):\n");
 855       printf("      key1=%i is out of range [1..300]\n",key1);
 856       return 0;
 857    }
 858    if ( key2<1 || key2>500 ) {
 859       printf("ERROR in THijing::GetPJTE(key1,key2):\n");
 860       printf("      key2=%i is out of range [1..500]\n",key2);
 861       return 0;
 862    }
 863
 864    return HIJJET1.pjte[key2-1][key1-1];
 865 }
 866
 867 //______________________________________________________________________________
 868 Float_t THijing::GetPJTM(Int_t key1, Int_t key2) const
 869 {
 870 // Get values of array PJTM in common HIJJET1
 871    if ( key1<1 || key1>300 ) {
 872       printf("ERROR in THijing::GetPJTM(key1):\n");
 873       printf("      key1=%i is out of range [1..300]\n",key1);
 874       return 0;
 875    }
 876    if ( key2<1 || key2>500 ) {
 877       printf("ERROR in THijing::GetPJTM(key1,key2):\n");
 878       printf("      key2=%i is out of range [1..500]\n",key2);
 879       return 0;
 880    }
 881
 882    return HIJJET1.pjtm[key2-1][key1-1];
 883 }
 884
 885 //====================== access to common HIJJET1 ===============================
 886
 887 //______________________________________________________________________________
 888 Int_t THijing::GetNSG() const
 889 {
 890 // Get value of NSG in common HIJJET2
 891    return HIJJET2.nsg;
 892 }
 893
 894 //______________________________________________________________________________
 895 Int_t THijing::GetNJSG(Int_t key) const
 896 {
 897 // Get values of array NJSG in common HIJJET2
 898    if ( key<1 || key>900 ) {
 899       printf ("ERROR in THijing:GetNJSG(key): \n ");
 900       printf ("      key=%i is out of range [1..900]!\n",key);
 901       return 0;
 902    }
 903
 904    return HIJJET2.njsg[key-1];
 905
 906 }
 907
 908 //______________________________________________________________________________
 909 Int_t THijing::GetIASG(Int_t key1, Int_t key2) const
 910 {
 911 // Get values of IASG in common HIJJET2
 912    if ( key1<1 || key1>900 ) {
 913       printf("ERROR in THijing::GetIASG(key1):\n");
 914       printf("      key1=%i is out of range [1..900]\n",key1);
 915       return 0;
 916    }
 917    if ( key2<1 || key2>3 ) {
 918       printf("ERROR in THijing::GetIASG(key1,key2):\n");
 919       printf("      key2=%i is out of range [1..3]\n",key2);
 920       return 0;
 921    }
 922
 923    return HIJJET2.iasg[key2-1][key1-1];
 924 }
 925
 926 //______________________________________________________________________________
 927 Int_t THijing::GetK1SG(Int_t key1, Int_t key2) const
 928 {
 929 // Get values of K1SG in common HIJJET2
 930    if ( key1<1 || key1>900 ) {
 931       printf("ERROR in THijing::GetK1SG(key1):\n");
 932       printf("      key1=%i is out of range [1..900]\n",key1);
 933       return 0;
 934    }
 935    if ( key2<1 || key2>100 ) {
 936       printf("ERROR in THijing::GetK1SG(key1,key2):\n");
 937       printf("      key2=%i is out of range [1..100]\n",key2);
 938       return 0;
 939    }
 940
 941    return HIJJET2.k1sg[key2-1][key1-1];
 942 }
 943
 944 //______________________________________________________________________________
 945 Int_t THijing::GetK2SG(Int_t key1, Int_t key2) const
 946 {
 947 // Get values of K2SG in common HIJJET2
 948    if ( key1<1 || key1>900 ) {
 949       printf("ERROR in THijing::GetK2SG(key1):\n");
 950       printf("      key1=%i is out of range [1..900]\n",key1);
 951       return 0;
 952    }
 953    if ( key2<1 || key2>100 ) {
 954       printf("ERROR in THijing::GetK2SG(key1,key2):\n");
 955       printf("      key2=%i is out of range [1..100]\n",key2);
 956       return 0;
 957    }
 958
 959    return HIJJET2.k2sg[key2-1][key1-1];
 960 }
 961
 962 //______________________________________________________________________________
 963 Float_t THijing::GetPXSG(Int_t key1, Int_t key2) const
 964 {
 965 // Get values of PXSG in common HIJJET2
 966    if ( key1<1 || key1>900 ) {
 967       printf("ERROR in THijing::GetPXSG(key1):\n");
 968       printf("      key1=%i is out of range [1..900]\n",key1);
 969       return 0;
 970    }
 971    if ( key2<1 || key2>100 ) {
 972       printf("ERROR in THijing::GetPXSG(key1,key2):\n");
 973       printf("      key2=%i is out of range [1..100]\n",key2);
 974       return 0;
 975    }
 976
 977    return HIJJET2.pxsg[key2-1][key1-1];
 978 }
 979
 980 //______________________________________________________________________________
 981 Float_t THijing::GetPYSG(Int_t key1, Int_t key2) const
 982 {
 983 // Get values of PYSG in common HIJJET2
 984    if ( key1<1 || key1>900 ) {
 985       printf("ERROR in THijing::GetPYSG(key1):\n");
 986       printf("      key1=%i is out of range [1..900]\n",key1);
 987       return 0;
 988    }
 989    if ( key2<1 || key2>100 ) {
 990       printf("ERROR in THijing::GetPYSG(key1,key2):\n");
 991       printf("      key2=%i is out of range [1..100]\n",key2);
 992       return 0;
 993    }
 994
 995    return HIJJET2.pysg[key2-1][key1-1];
 996 }
 997
 998 //______________________________________________________________________________
 999 Float_t THijing::GetPZSG(Int_t key1, Int_t key2) const
1000 {
1001 // Get values of PZSG in common HIJJET2
1002    if ( key1<1 || key1>900 ) {
1003       printf("ERROR in THijing::GetPZSG(key1):\n");
1004       printf("      key1=%i is out of range [1..900]\n",key1);
1005       return 0;
1006    }
1007    if ( key2<1 || key2>100 ) {
1008       printf("ERROR in THijing::GetPZSG(key1,key2):\n");
1009       printf("      key2=%i is out of range [1..100]\n",key2);
1010       return 0;
1011    }
1012
1013    return HIJJET2.pzsg[key2-1][key1-1];
1014 }
1015
1016 //______________________________________________________________________________
1017 Float_t THijing::GetPESG(Int_t key1, Int_t key2) const
1018 {
1019 // Get values of PESG in common HIJJET2
1020    if ( key1<1 || key1>900 ) {
1021       printf("ERROR in THijing::GetPESG(key1):\n");
1022       printf("      key1=%i is out of range [1..900]\n",key1);
1023       return 0;
1024    }
1025    if ( key2<1 || key2>100 ) {
1026       printf("ERROR in THijing::GetPESG(key1,key2):\n");
1027       printf("      key2=%i is out of range [1..100]\n",key2);
1028       return 0;
1029    }
1030
1031    return HIJJET2.pesg[key2-1][key1-1];
1032 }
1033
1034 //______________________________________________________________________________
1035 Float_t THijing::GetPMSG(Int_t key1, Int_t key2) const
1036 {
1037 // Get values of PMSG in common HIJJET2
1038    if ( key1<1 || key1>900 ) {
1039       printf("ERROR in THijing::GetPMSG(key1):\n");
1040       printf("      key1=%i is out of range [1..900]\n",key1);
1041       return 0;
1042    }
1043    if ( key2<1 || key2>100 ) {
1044       printf("ERROR in THijing::GetPMSG(key1,key2):\n");
1045       printf("      key2=%i is out of range [1..100]\n",key2);
1046       return 0;
1047    }
1048
1049    return HIJJET2.pmsg[key2-1][key1-1];
1050 }
1051
1052 //====================== access to common HISTRNG ===============================
1053
1054 //______________________________________________________________________________
1055 Int_t THijing::GetNFP(Int_t key1, Int_t key2) const
1056 {
1057 // Get values of array NFP in common HISTRNG
1058    if ( key1<1 || key1>300 ) {
1059       printf("ERROR in THijing::GetNFP(key1):\n");
1060       printf("      key1=%i is out of range [1..300]\n",key1);
1061       return 0;
1062    }
1063    if ( key2<1 || key2>15 ) {
1064       printf("ERROR in THijing::GetNFP(key1,key2):\n");
1065       printf("      key2=%i is out of range [1..15]\n",key2);
1066       return 0;
1067    }
1068
1069    return HISTRNG.nfp[key2-1][key1-1];
1070 }
1071
1072 //______________________________________________________________________________
1073 Float_t THijing::GetPP(Int_t key1, Int_t key2) const
1074 {
1075 // Get values of array PP in common HISTRNG
1076    if ( key1<1 || key1>300 ) {
1077       printf("ERROR in THijing::GetPP(key1):\n");
1078       printf("      key1=%i is out of range [1..300]\n",key1);
1079       return 0;
1080    }
1081    if ( key2<1 || key2>15 ) {
1082       printf("ERROR in THijing::GetPP(key1,key2):\n");
1083       printf("      key2=%i is out of range [1..15]\n",key2);
1084       return 0;
1085    }
1086
1087    return HISTRNG.pp[key2-1][key1-1];
1088 }
1089
1090 //______________________________________________________________________________
1091 Int_t THijing::GetNFT(Int_t key1, Int_t key2) const
1092 {
1093 // Get values of array NFT in common HISTRNG
1094    if ( key1<1 || key1>300 ) {
1095       printf("ERROR in THijing::GetNFT(key1):\n");
1096       printf("      key1=%i is out of range [1..300]\n",key1);
1097       return 0;
1098    }
1099    if ( key2<1 || key2>15 ) {
1100       printf("ERROR in THijing::GetNFT(key1,key2):\n");
1101       printf("      key2=%i is out of range [1..15]\n",key2);
1102       return 0;
1103    }
1104
1105    return HISTRNG.nft[key2-1][key1-1];
1106 }
1107
1108 //______________________________________________________________________________
1109 Float_t THijing::GetPT(Int_t key1, Int_t key2) const
1110 {
1111 // Get values of array PT in common HISTRNG
1112    if ( key1<1 || key1>300 ) {
1113       printf("ERROR in THijing::GetPT(key1):\n");
1114       printf("      key1=%i is out of range [1..300]\n",key1);
1115       return 0;
1116    }
1117    if ( key2<1 || key2>15 ) {
1118       printf("ERROR in THijing::GetPT(key1,key2):\n");
1119       printf("      key2=%i is out of range [1..15]\n",key2);
1120       return 0;
1121    }
1122
1123    return HISTRNG.pt[key2-1][key1-1];
1124 }
1125
1126 void THijing::SetPARJ(Int_t key, Float_t parm)
1127 {
1128 // Set values of array PARJ in common HISTRNG
1129     if ( key < 1 || key > 200) {
1130         printf("ERROR in THijing::SetPARJ(key,parm):\n");
1131         printf("      key=%i is out of range [1..200]\n",key);
1132     } else {
1133       LUDAT1_HIJING.parj[key-1] = parm;
1134     }
1135 }
1136
1137
1138 void THijing::SetMSTJ(Int_t key, Int_t parm)
1139 {
1140 // Set values of array MSTJ in common HISTRNG
1141     if ( key < 1 || key > 200) {
1142         printf("ERROR in THijing::SetMSTJ(key,parm):\n");
1143         printf("      key=%i is out of range [1..200]\n",key);
1144     } else {
1145       LUDAT1_HIJING.mstj[key-1] = parm;
1146     }
1147 }
1148
1149
1150 //====================== access to Hijing subroutines =========================
1151
1152
1153 //______________________________________________________________________________
1154 void THijing::Initialize()
1155 {
1156 //////////////////////////////////////////////////////////////////////////////////
1157 // Calls Hijset with the either default parameters or the ones set by the user  //
1158 // via SetEFRM, SetFRAME, SetPROJ, SetTARG, SetIAP, SetIZP, SetIAT, SetIZT      //
1159 //////////////////////////////////////////////////////////////////////////////////
1160
1161    if ( (!strcmp(fFrame.Data(), "CMS     "  )) &&
1162         (!strcmp(fFrame.Data(), "LAB     "  ))){
1163       printf("WARNING! In THijing:Initialize():\n");
1164       printf(" specified frame=%s is neither CMS or LAB\n",fFrame.Data());
1165       printf(" resetting to default \"CMS\" .");
1166       fFrame="CMS";
1167    }
1168
1169    if ( (!strcmp(fProj.Data(), "A       "     )) &&
1170         (!strcmp(fProj.Data(), "P       "     )) &&
1171         (!strcmp(fProj.Data(), "PBAR    "  ))){
1172       printf("WARNING! In THijing:Initialize():\n");
1173       printf(" specified projectile=%s is neither A, P or PBAR\n",fProj.Data());
1174       printf(" resetting to default \"A\" .");
1175       fProj="A";
1176    }
1177
1178    if ( (!strcmp(fTarg.Data(), "A       "     )) &&
1179         (!strcmp(fTarg.Data(), "P       "     )) &&
1180         (!strcmp(fTarg.Data(), "PBAR    "  ))){
1181       printf("WARNING! In THijing:Initialize():\n");
1182       printf(" specified target=%s is neither A, P or PBAR\n",fTarg.Data());
1183       printf(" resetting to default \"A\" .");
1184       fTarg="A";
1185    }
1186
1187    printf(" %s-%s at %f GeV \n",fProj.Data(),fTarg.Data(),fEfrm);
1188
1189    Hijset(fEfrm,fFrame.Data(),fProj.Data(),fTarg.Data(),fIap,fIzp,fIat,fIzt);
1190
1191    printf(" %s-%s at %f GeV \n",fProj.Data(),fTarg.Data(),fEfrm);
1192 }
1193
1194
1195 //______________________________________________________________________________
1196 void THijing::GenerateEvent()
1197 {
1198 // Generates one event;
1199
1200    Hijing(fFrame.Data(),fBmin,fBmax);
1201
1202 }
1203 //______________________________________________________________________________
1204 void THijing::Hijset(float efrm, const char *frame, const char *proj,
1205                      const char *targ, int iap, int izp, int iat, int izt)
1206 {
1207 // Call HIJING routine HIJSET passing the parameters in a way accepted by
1208 // Fortran routines
1209
1210   int s1 = strlen(frame);
1211   int s2 = strlen(proj);
1212   int s3 = strlen(targ);
1213   printf("s1 = %d s2 = %d s3 = %d\n",s1,s2,s3);
1214 #ifndef WIN32
1215   hijset(efrm, frame, proj, targ, iap, izp, iat, izt, s1, s2, s3);
1216 #else
1217   hijset(efrm, frame, s1, proj, s2, targ, s3, iap, izp, iat, izt);
1218 #endif
1219 }
1220 //______________________________________________________________________________
1221 void THijing::Hijing(const char *frame, float bmin, float bmax)
1222 {
1223 // Call HIJING routine HIJSET passing the parameters in a way accepted by
1224 // Fortran routines
1225
1226   int s1 = strlen(frame);
1227
1228 #ifndef WIN32
1229   hijing(frame, bmin, bmax, s1);
1230 #else
1231   hijing(frame, s1, bmin, bmax);
1232 #endif
1233 }
1234
1235
1236 Float_t  THijing::Profile(float b)
1237 {
1238 // Call HIJING routine PROFILE
1239   return profile(b);
1240 }
1241
1242
1243 void  THijing::Rluget(Int_t lfn, Int_t move)
1244 {
1245 // write seed to file
1246   rluget_hijing(lfn, move);
1247 }
1248
1249
1250 void  THijing::Rluset(Int_t lfn, Int_t move)
1251 {
1252 // read seed from file
1253   rluset_hijing(lfn, move);
1254 }
1255