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