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