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