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