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