TAmpt/AliGenAmpt.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 // Generator using AMPT as an external generator
  19
  20 #include "AliGenAmpt.h"
  21
  22 #include <TClonesArray.h>
  23 #include <TGraph.h>
  24 #include <TAmpt.h>
  25 #include <TLorentzVector.h>
  26 #include <TPDGCode.h>
  27 #include <TParticle.h>
  28 #include <TVirtualMC.h>
  29 #include <TParticlePDG.h>
  30 #include "AliGenHijingEventHeader.h"
  31 #define AliGenAmptEventHeader AliGenHijingEventHeader
  32 #include "AliAmptRndm.h"
  33 #include "AliLog.h"
  34 #include "AliRun.h"
  35 #include "AliDecayer.h"
  36
  37 ClassImp(AliGenAmpt)
  38
  39 AliGenAmpt::AliGenAmpt()
  40   : AliGenMC(),
  41     fDecayer(NULL),
  42     fFrame("CMS"),
  43     fMinImpactParam(0.),
  44     fMaxImpactParam(5.),
  45     fKeep(0),
  46     fQuench(0),
  47     fShadowing(1),
  48     fDecaysOff(1),
  49     fTrigger(0),
  50     fEvaluate(0),
  51     fSelectAll(0),
  52     fFlavor(0),
  53     fKineBias(0.),
  54     fTrials(0),
  55     fXsection(0.),
  56     fAmpt(0),
  57     fPtHardMin(2.0),
  58     fPtHardMax(-1),
  59     fSpectators(1),
  60     fDsigmaDb(0),
  61     fDnDb(0),
  62     fPtMinJet(-2.5),
  63     fEtaMinJet(-20.),
  64     fEtaMaxJet(+20.),
  65     fPhiMinJet(0.),
  66     fPhiMaxJet(TMath::TwoPi()),
  67     fRadiation(3),
  68     fSimpleJet(kFALSE),
  69     fNoGammas(kFALSE),
  70     fProjectileSpecn(0),
  71     fProjectileSpecp(0),
  72     fTargetSpecn(0),
  73     fTargetSpecp(0),
  74     fLHC(kFALSE),
  75     fRandomPz(kFALSE),
  76     fNoHeavyQuarks(kFALSE),
  77     fIsoft(4),
  78     fNtMax(150),
  79     fIpop(1),
  80     fXmu(3.2264),
  81     fAlpha(1./3),
  82     fStringA(0.5),
  83     fStringB(0.9),
  84     fEventTime(0.),
  85     fHeader(new AliGenAmptEventHeader("Ampt")),
  86     fDecay(kTRUE),
  87     fRotating(kFALSE)
  88 {
  89   // Constructor
  90   fEnergyCMS = 2760.;
  91   AliAmptRndm::SetAmptRandom(GetRandom());
  92 }
  93
  94 AliGenAmpt::AliGenAmpt(Int_t npart)
  95   : AliGenMC(npart),
  96     fDecayer(NULL),
  97     fFrame("CMS"),
  98     fMinImpactParam(0.),
  99     fMaxImpactParam(5.),
 100     fKeep(0),
 101     fQuench(0),
 102     fShadowing(1),
 103     fDecaysOff(1),
 104     fTrigger(0),
 105     fEvaluate(0),
 106     fSelectAll(0),
 107     fFlavor(0),
 108     fKineBias(0.),
 109     fTrials(0),
 110     fXsection(0.),
 111     fAmpt(0),
 112     fPtHardMin(2.0),
 113     fPtHardMax(-1),
 114     fSpectators(1),
 115     fDsigmaDb(0),
 116     fDnDb(0),
 117     fPtMinJet(-2.5),
 118     fEtaMinJet(-20.),
 119     fEtaMaxJet(+20.),
 120     fPhiMinJet(0.),
 121     fPhiMaxJet(2. * TMath::Pi()),
 122     fRadiation(3),
 123     fSimpleJet(kFALSE),
 124     fNoGammas(kFALSE),
 125     fProjectileSpecn(0),
 126     fProjectileSpecp(0),
 127     fTargetSpecn(0),
 128     fTargetSpecp(0),
 129     fLHC(kFALSE),
 130     fRandomPz(kFALSE),
 131     fNoHeavyQuarks(kFALSE),
 132     fIsoft(1),
 133     fNtMax(150),
 134     fIpop(1),
 135     fXmu(3.2264),
 136     fAlpha(1./3),
 137     fStringA(0.5),
 138     fStringB(0.9),
 139     fEventTime(0.),
 140     fHeader(new AliGenAmptEventHeader("Ampt")),
 141     fDecay(kTRUE),
 142     fRotating(kFALSE)
 143 {
 144   // Default PbPb collisions at 2.76 TeV
 145
 146   fEnergyCMS = 2760.;
 147   fName = "Ampt";
 148   fTitle= "Particle Generator using AMPT";
 149   AliAmptRndm::SetAmptRandom(GetRandom());
 150 }
 151
 152 AliGenAmpt::~AliGenAmpt()
 153 {
 154   // Destructor
 155   if ( fDsigmaDb) delete fDsigmaDb;
 156   if ( fDnDb)     delete fDnDb;
 157   if ( fHeader)   delete fHeader;
 158 }
 159
 160 void AliGenAmpt::Init()
 161 {
 162   // Initialisation
 163
 164   fFrame.Resize(8);
 165   fTarget.Resize(8);
 166   fProjectile.Resize(8);
 167
 168   fAmpt = new TAmpt(fEnergyCMS, fFrame, fProjectile, fTarget,
 169                     fAProjectile, fZProjectile, fATarget, fZTarget,
 170                     fMinImpactParam, fMaxImpactParam);
 171   SetMC(fAmpt);
 172
 173   fAmpt->SetIHPR2(2,  fRadiation);
 174   fAmpt->SetIHPR2(3,  fTrigger);
 175   fAmpt->SetIHPR2(6,  fShadowing);
 176   fAmpt->SetIHPR2(12, fDecaysOff);
 177   fAmpt->SetIHPR2(21, fKeep);
 178   fAmpt->SetHIPR1(8,  fPtHardMin);
 179   fAmpt->SetHIPR1(9,  fPtHardMax);
 180   fAmpt->SetHIPR1(10, fPtMinJet);
 181   fAmpt->SetHIPR1(50, fSimpleJet);
 182
 183   //  Quenching
 184   //  fQuench = 0:  no quenching
 185   //  fQuench = 1:  Hijing default
 186   //  fQuench = 2:  new LHC  parameters for HIPR1(11) and HIPR1(14)
 187   //  fQuench = 3:  new RHIC parameters for HIPR1(11) and HIPR1(14)
 188   //  fQuench = 4:  new LHC  parameters with log(e) dependence
 189   //  fQuench = 5:  new RHIC parameters with log(e) dependence
 190   fAmpt->SetIHPR2(50, 0);
 191   if (fQuench > 0)
 192     fAmpt->SetIHPR2(4,  1);
 193   else
 194     fAmpt->SetIHPR2(4,  0);
 195
 196   if (fQuench == 2) {
 197     fAmpt->SetHIPR1(14, 1.1);
 198     fAmpt->SetHIPR1(11, 3.7);
 199   } else if (fQuench == 3) {
 200     fAmpt->SetHIPR1(14, 0.20);
 201     fAmpt->SetHIPR1(11, 2.5);
 202   } else if (fQuench == 4) {
 203     fAmpt->SetIHPR2(50, 1);
 204     fAmpt->SetHIPR1(14, 4.*0.34);
 205     fAmpt->SetHIPR1(11, 3.7);
 206   } else if (fQuench == 5) {
 207     fAmpt->SetIHPR2(50, 1);
 208     fAmpt->SetHIPR1(14, 0.34);
 209     fAmpt->SetHIPR1(11, 2.5);
 210   }
 211
 212   // Heavy quarks
 213   if (fNoHeavyQuarks) {
 214     fAmpt->SetIHPR2(49, 1);
 215   } else {
 216     fAmpt->SetIHPR2(49, 0);
 217   }
 218
 219   // Ampt specific
 220   fAmpt->SetIsoft(fIsoft);
 221   fAmpt->SetNtMax(fNtMax);
 222   fAmpt->SetIpop(fIpop);
 223   fAmpt->SetXmu(fXmu);
 224   fAmpt->SetAlpha(fAlpha);
 225   fAmpt->SetStringFrag(fStringA, fStringB);
 226
 227   AliGenMC::Init();
 228
 229   // Initialize Ampt
 230   fAmpt->Initialize();
 231   if (fEvaluate)
 232     EvaluateCrossSections();
 233
 234   fAmpt->SetReactionPlaneAngle(0.0);
 235   fRotating=kFALSE;
 236 }
 237
 238 void AliGenAmpt::Generate()
 239 {
 240   // Generate one event
 241
 242   Float_t polar[3]    =   {0,0,0};
 243   Float_t origin[3]   =   {0,0,0};
 244   Float_t origin0[3]  =   {0,0,0};
 245   Float_t time0 = 0.;
 246   Float_t p[3];
 247   Float_t tof;
 248
 249   Int_t nt  = 0;
 250   Int_t jev = 0;
 251   Int_t j, kf, ks, ksp, imo;
 252   kf = 0;
 253
 254   fTrials = 0;
 255   for (j = 0;j < 3; j++)
 256     origin0[j] = fOrigin[j];
 257   //time0 = fTimeOrigin;
 258
 259   if(fVertexSmear == kPerEvent) {
 260     Vertex();
 261     for (j=0; j < 3; j++)
 262       origin0[j] = fVertex[j];
 263     //time0 = fTime;
 264   }
 265
 266   Float_t sign = (fRandomPz && (Rndm() < 0.5))? -1. : 1.;
 267
 268   while(1) {
 269
 270     // Generate random reaction plane angle if requested
 271     if( fRotating ) {
 272       TRandom *r=AliAmptRndm::GetAmptRandom();
 273       fAmpt->SetReactionPlaneAngle(TMath::TwoPi()*r->Rndm());
 274     }
 275
 276     // Generate one event
 277     Int_t fpemask = gSystem->GetFPEMask();
 278     gSystem->SetFPEMask(0);
 279     fAmpt->GenerateEvent();
 280     gSystem->SetFPEMask(fpemask);
 281     fTrials++;
 282     fNprimaries = 0;
 283
 284
 285     fAmpt->ImportParticles(&fParticles,"All");
 286     Int_t np = fParticles.GetEntriesFast();
 287     if (np == 0 )
 288       continue;
 289     //
 290     //RS>>: Decayers now returns cm and sec. Since TAmpt returns mm and mm/c, convert its result to cm and sec here
 291     const Float_t kconvT=0.001/2.999792458e8; // mm/c to seconds conversion
 292     const Float_t kconvL=1./10; // mm to cm conversion
 293     for (int ip=np;ip--;) {
 294       TParticle* part = (TParticle*)fParticles[ip];
 295       if (!part) continue;
 296       part->SetProductionVertex(part->Vx()*kconvL,part->Vy()*kconvL,part->Vz()*kconvL,kconvT*part->T());
 297     }
 298     // RS<<
 299     //
 300     if (fTrigger != kNoTrigger) {
 301       if (!CheckTrigger())
 302         continue;
 303     }
 304
 305     AliDecayer *decayer = 0;
 306     //if (gMC)
 307     //  decayer = gMC->GetDecayer();
 308     decayer = fDecayer; //AMPT does not do the strong decays per dafault
 309
 310     if (decayer&&fDecay) {
 311       TClonesArray arr("TParticle",100);
 312       for( Int_t nLoop=0; nLoop!=2; ++nLoop) { // In order to produce more than one generation of decays: NumberOfNestedLoops set to 2
 313         Int_t np2 = np;
 314             for (Int_t i = 0; i < np; i++) {
 315                 TParticle *iparticle = (TParticle *)fParticles.At(i);
 316                 if (!Stable(iparticle)) // true if particle has daughters already
 317                   continue;
 318                 kf = TMath::Abs(iparticle->GetPdgCode());
 319                 if (kf==92)
 320                   continue;
 321           if( !IsThisAKnownParticle(iparticle) ) continue; // skip undesired particles
 322                 /*
 323                 if (0) { // this turned out to be too cumbersome!
 324                   if (kf!=331&&kf!=3114&&kf!=3114&&kf!=411&&kf!=-4122&&kf!=-3324&&kf!=-3312&&kf!=-3114&&
 325                       kf!=-311&&kf!=3214&&kf!=-3214&&kf!=-433&&kf!=413&&kf!=3122&&kf!=-3122&&kf!=-413&&
 326                       kf!=-421&&kf!=-423&&kf!=3324&&kf!=-313&&kf!=213&&kf!=-213&&kf!=3314&&kf!=3222&&
 327                       kf!=-3222&&kf!=3224&&kf!=-3224&&kf!=-4212&&kf!=4212&&kf!=433&&kf!=423&&kf!=-3322&&
 328                       kf!=3322&&kf!=-3314)
 329                     continue; //decay eta',Sigma*+,Sigma*-,D+,Lambda_c-,Xi*0_bar,Xi-_bar,Sigma*-,
 330                               //      K0_bar,Sigma*0,Sigma*0_bar,D*_s-,D*+,Lambda0,Lambda0_bar,D*-
 331                               //      D0_bar,D*0_bar,Xi*0,K*0_bar,rho+,rho-,Xi*-,Sigma-,
 332                             //      Sigma+,Sigma*+,Sigma*-,Sigma_c-,Sigma_c+,D*_s+,D*0,Xi0_bar
 333                             //      Xi0,Xi*+
 334               //} else { // really only decay particles if there are not known to Geant3
 335               //  if (gMC->IdFromPDG(kf)>0)
 336               //    continue;
 337               }
 338               if (0) { // defining the particle for Geant3 leads to a floating point exception.
 339                 TParticlePDG *pdg = iparticle->GetPDG(1);
 340                 //pdg->Print(); printf("%s\n",pdg->ParticleClass());
 341                 TString ptype(pdg->ParticleClass());
 342                 TMCParticleType mctype(kPTUndefined);
 343                 if (ptype=="Baryon" || ptype=="Meson")
 344                   mctype = kPTHadron;
 345                 gMC->DefineParticle(pdg->PdgCode(), pdg->GetName(), mctype, pdg->Mass(), pdg->Charge(), pdg->Lifetime(),
 346                                     ptype,pdg->Width(), (Int_t)pdg->Spin(), (Int_t)pdg->Parity(), 0,
 347                                     (Int_t)pdg->Isospin(), 0, 0, 0, 0, pdg->Stable());
 348                 gMC->SetUserDecay(pdg->PdgCode());
 349                 continue;
 350               }
 351               */
 352                 TLorentzVector pmom(iparticle->Px(),iparticle->Py(),iparticle->Pz(),iparticle->Energy());
 353                 decayer->Decay(kf,&pmom);
 354                 decayer->ImportParticles(&arr);
 355                 Int_t ndecayed = arr.GetEntries();
 356                 if (ndecayed>1) {
 357                   if (np2+ndecayed>fParticles.GetSize())
 358                     fParticles.Expand(2*fParticles.GetSize());
 359                   //arr.Print();
 360                   // iparticle->SetStatusCode(2);  to be compatible with Hijing
 361                   iparticle->SetFirstDaughter(np2);
 362                   for (Int_t jj = 1; jj < ndecayed; jj++) {
 363                   TParticle *jp = (TParticle *)arr.At(jj);
 364                     if (jp->GetFirstMother()!=1)
 365                       continue;
 366
 367                     TParticle *newp = new(fParticles[np2]) TParticle(jp->GetPdgCode(),
 368                                                                    0, //1,  //to be compatible with Hijing
 369                                                                    i,
 370                                                                    -1,
 371                                                                    -1,
 372                                                                    -1,
 373                                                                    jp->Px(),jp->Py(),jp->Pz(),jp->Energy(),
 374                                                                      jp->Vx(),jp->Vy(),jp->Vz(),jp->T());
 375                     //take care of the phi
 376                     //if((kf == 333)||(kf == 313)) {
 377                     if(IsThisAKnownParticle(iparticle)) {
 378                       //Printf("=============PANOS===================");
 379                       //Printf("Phi detected - daughet is: %d",jp->GetPdgCode());
 380                       newp->SetUniqueID(4);
 381                     }
 382                     else newp->SetUniqueID( jp->GetStatusCode() );
 383                   np2++;
 384                 } // end of jj->nDecayedParticles
 385                 iparticle->SetLastDaughter(np2-1);
 386               } // end of nDecayedPrticles>1
 387             } // end of i->np
 388         np = fParticles.GetEntries();
 389         if (np!=np2) {
 390           AliError(Form("Something is fishy: %d %d\n", np,np2));
 391         }
 392       } // end of nLoop->NumberOfNestedLoops
 393     } else {
 394       if (fDecay)
 395         AliError("No decayer found, but fDecay==kTRUE!");
 396     }
 397
 398     if (fLHC)
 399       Boost();
 400
 401     Int_t nc = 0;
 402     Int_t* newPos     = new Int_t[np];
 403     Int_t* pSelected  = new Int_t[np];
 404
 405     for (Int_t i = 0; i < np; i++) {
 406       newPos[i]    = i;
 407       pSelected[i] = 0;
 408     }
 409
 410     // Get event vertex
 411     //TParticle *  iparticle = (TParticle *) fParticles.At(0);
 412     fVertex[0] = origin0[0];
 413     fVertex[1] = origin0[1];
 414     fVertex[2] = origin0[2];
 415     //fTime = time0;
 416
 417     // First select parent particles
 418     for (Int_t i = 0; i < np; i++) {
 419       TParticle *iparticle = (TParticle *) fParticles.At(i);
 420
 421     // Is this a parent particle ?
 422       if (Stable(iparticle)) continue;  // quit if particle has no daughters
 423       Bool_t  selected             =  kTRUE;
 424       Bool_t  hasSelectedDaughters =  kFALSE;
 425       kf = iparticle->GetPdgCode();
 426       ks = iparticle->GetStatusCode();
 427       if (kf == 92)
 428         continue;
 429
 430       if (!fSelectAll)
 431         selected = KinematicSelection(iparticle, 0) && SelectFlavor(kf);
 432       hasSelectedDaughters = DaughtersSelection(iparticle);
 433
 434       // Put particle on the stack if it is either selected or
 435       // it is the mother of at least one seleted particle
 436       if (selected || hasSelectedDaughters) {
 437         nc++;
 438         pSelected[i] = 1;
 439       } // selected
 440     } // particle loop parents
 441
 442     // Now select the final state particles
 443     fProjectileSpecn    = 0;
 444     fProjectileSpecp    = 0;
 445     fTargetSpecn        = 0;
 446     fTargetSpecp        = 0;
 447     for (Int_t i = 0; i<np; i++) {
 448       TParticle *iparticle = (TParticle *) fParticles.At(i);
 449       // Is this a final state particle ?
 450       if (!Stable(iparticle)) continue;  // quit if particle has daughters
 451       Bool_t  selected =  kTRUE;
 452       kf = iparticle->GetPdgCode();
 453       if (kf == 92)
 454         continue;
 455       ks  = iparticle->GetStatusCode();
 456       ksp = iparticle->GetUniqueID();
 457
 458       // --------------------------------------------------------------------------
 459       // Count spectator neutrons and protons
 460       if(ksp == 0 || ksp == 1) {
 461         if(kf == kNeutron) fProjectileSpecn += 1;
 462         if(kf == kProton)  fProjectileSpecp += 1;
 463       } else if(ksp == 10 || ksp == 11) {
 464         if(kf == kNeutron) fTargetSpecn += 1;
 465         if(kf == kProton)  fTargetSpecp += 1;
 466       }
 467       // --------------------------------------------------------------------------
 468       if (!fSelectAll) {
 469         selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
 470         if (!fSpectators && selected)
 471           selected = (ksp != 0 && ksp != 1 && ksp != 10 && ksp != 11);
 472       }
 473
 474      // Put particle on the stack if selected
 475       if (selected) {
 476         nc++;
 477         pSelected[i] = 1;
 478         if (0) printf("---> %d %d %d %s\n",i,nc,kf,iparticle->GetName());
 479       } // selected
 480     } // particle loop final state
 481
 482     // Write particles to stack
 483     for (Int_t i = 0; i<np; i++) {
 484       if (pSelected[i]) {
 485               TParticle *iparticle = (TParticle *) fParticles.At(i);
 486               Bool_t  hasMother   = (iparticle->GetFirstMother()     >=0);
 487               Bool_t  hasDaughter = (iparticle->GetFirstDaughter()   >=0);
 488         kf   = iparticle->GetPdgCode();
 489         ks   = iparticle->GetStatusCode();
 490         p[0] = iparticle->Px();
 491         p[1] = iparticle->Py();
 492         p[2] = iparticle->Pz() * sign;
 493         origin[0] = origin0[0]+iparticle->Vx();
 494         origin[1] = origin0[1]+iparticle->Vy();
 495         origin[2] = origin0[2]+iparticle->Vz();
 496         tof = time0 + iparticle->T();
 497
 498         imo = -1;
 499         TParticle* mother = 0;
 500         TMCProcess procID = (TMCProcess) iparticle->GetUniqueID();
 501         if (hasMother) {
 502           imo = iparticle->GetFirstMother();
 503           mother = (TParticle *) fParticles.At(imo);
 504           imo = (mother->GetPdgCode() != 92) ? newPos[imo] : -1;
 505         } else { // if has no mothers then it was created by AMPT
 506           if(procID==999)
 507             procID = kPPrimary; // reseting to ALIROOT convention
 508           else
 509             procID = kPNoProcess; // for expectators
 510         } // if has mother
 511         Bool_t tFlag = (fTrackIt && !hasDaughter);
 512         PushTrack(tFlag,imo,kf,p,origin,polar,tof,procID,nt, 1., ks);
 513         fNprimaries++;
 514         KeepTrack(nt);
 515         newPos[i] = nt;
 516       } // if selected
 517     } // particle loop
 518     delete[] newPos;
 519     delete[] pSelected;
 520
 521     AliInfo(Form("\n I've put %i particles on the stack \n",nc));
 522     if (nc > 0) {
 523       jev += nc;
 524       if (jev >= fNpart || fNpart == -1) {
 525         fKineBias = Float_t(fNpart)/Float_t(fTrials);
 526         AliInfo(Form("\n Trials: %i %i %i\n",fTrials, fNpart, jev));
 527         break;
 528       }
 529     }
 530   } // event loop
 531   MakeHeader();
 532   SetHighWaterMark(nt);
 533 }
 534
 535 Bool_t AliGenAmpt::IsThisAKnownParticle(TParticle *thisGuy)
 536 {
 537   // In order to prevent AMPT to introduce weird particles into the decayer and transporter
 538   // blame cperez@cern.ch for this method
 539
 540   Int_t pdgcode = TMath::Abs( thisGuy->GetPdgCode() );
 541
 542   Int_t myFavoriteParticles[ 38] = { 3322, 3314, 3312, 3224, 3222,  // Xi0       Xi*+-   Xi+-    Sigma*-+ Sigma-+
 543                                      3214, 3212, 3122, 3114, 3112,  // Sigma*0   Sigma0  Lambda0 Sigma*+- Sigma+-
 544                                      2224, 2214, 2212, 2114, 2112,  // Delta--++ Delta-+ proton  Delta0   neutron
 545                                      1114,  323,  321,  313,  311,  // Delta+-   K*-+    K-+     K*0      K0
 546                                       213,  211,   11,   22,  111,  // rho-+     pi-+    e+-     gamma    pi0
 547                                       113,  130,  221,  223,  310,  // rho0      K_L0    eta     omega    K_S0
 548                                       331,  333, 3324,  431,  421,  // eta'      phi     Xi*0    Ds-+     D0
 549                                       411,  413,   13               // D-+       D*-+    mu+-
 550                                     };
 551
 552   Bool_t found = kFALSE;
 553   for(Int_t i=0; i!=38; ++i)
 554     if( myFavoriteParticles[i] == pdgcode ) {
 555       found = kTRUE;
 556       break;
 557     }
 558
 559   return found;
 560 }
 561
 562 void AliGenAmpt::EvaluateCrossSections()
 563 {
 564   // Glauber Calculation of geometrical x-section
 565
 566   Float_t xTot       = 0.;          // barn
 567   Float_t xTotHard   = 0.;          // barn
 568   Float_t xPart      = 0.;          // barn
 569   Float_t xPartHard  = 0.;          // barn
 570   Float_t sigmaHard  = 0.1;         // mbarn
 571   Float_t bMin       = 0.;
 572   Float_t bMax       = fAmpt->GetHIPR1(34)+fAmpt->GetHIPR1(35);
 573   const Float_t kdib = 0.2;
 574   Int_t   kMax       = Int_t((bMax-bMin)/kdib)+1;
 575
 576   printf("\n Projectile Radius (fm): %f \n",fAmpt->GetHIPR1(34));
 577   printf("\n Target     Radius (fm): %f \n",fAmpt->GetHIPR1(35));
 578
 579   Int_t i;
 580   Float_t oldvalue= 0.;
 581   Float_t* b   = new Float_t[kMax]; memset(b,0,kMax*sizeof(Float_t));
 582   Float_t* si1 = new Float_t[kMax]; memset(si1,0,kMax*sizeof(Float_t));
 583   Float_t* si2 = new Float_t[kMax]; memset(si2,0,kMax*sizeof(Float_t));
 584   for (i = 0; i < kMax; i++) {
 585     Float_t xb  = bMin+i*kdib;
 586     Float_t ov=fAmpt->Profile(xb);
 587     Float_t gb  =  2.*0.01*fAmpt->GetHIPR1(40)*kdib*xb*(1.-TMath::Exp(-fAmpt->GetHINT1(12)*ov));
 588     Float_t gbh =  2.*0.01*fAmpt->GetHIPR1(40)*kdib*xb*sigmaHard*ov;
 589     xTot+=gb;
 590     xTotHard += gbh;
 591     printf("profile %f %f %f\n", xb, ov, fAmpt->GetHINT1(12));
 592
 593     if (xb > fMinImpactParam && xb < fMaxImpactParam) {
 594       xPart += gb;
 595       xPartHard += gbh;
 596     }
 597
 598     if ((oldvalue) && ((xTot-oldvalue)/oldvalue<0.0001))
 599       break;
 600     oldvalue = xTot;
 601     printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot);
 602     printf("\n Hard  cross section (barn): %d %f %f \n\n",i, xb, xTotHard);
 603     if (i>0) {
 604       si1[i] = gb/kdib;
 605       si2[i] = gbh/gb;
 606       b[i]  = xb;
 607     }
 608   }
 609
 610   printf("\n Total cross section (barn): %f \n",xTot);
 611   printf("\n Hard  cross section (barn): %f \n \n",xTotHard);
 612   printf("\n Partial       cross section (barn): %f %f \n",xPart, xPart/xTot*100.);
 613   printf("\n Partial  hard cross section (barn): %f %f \n",xPartHard, xPartHard/xTotHard*100.);
 614
 615   //  Store result as a graph
 616   b[0] = 0;
 617   si1[0] = 0;
 618   si2[0]=si2[1];
 619   delete fDsigmaDb;
 620   fDsigmaDb  = new TGraph(i, b, si1);
 621   delete fDnDb;
 622   fDnDb      = new TGraph(i, b, si2);
 623 }
 624
 625 Bool_t AliGenAmpt::DaughtersSelection(TParticle* iparticle)
 626 {
 627   // Looks recursively if one of the daughters has been selected
 628   //printf("\n Consider daughters %d:",iparticle->GetPdgCode());
 629   Int_t imin = -1;
 630   Int_t imax = -1;
 631   Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
 632   Bool_t selected = kFALSE;
 633   if (hasDaughters) {
 634     imin = iparticle->GetFirstDaughter();
 635     imax = iparticle->GetLastDaughter();
 636     for (Int_t i = imin; i <= imax; i++){
 637       TParticle *  jparticle = (TParticle *) fParticles.At(i);
 638       Int_t ip = jparticle->GetPdgCode();
 639       if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
 640         selected=kTRUE; break;
 641       }
 642       if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
 643     }
 644   } else {
 645     return kFALSE;
 646   }
 647   return selected;
 648 }
 649
 650 Bool_t AliGenAmpt::SelectFlavor(Int_t pid)
 651 {
 652   // Select flavor of particle
 653   // 0: all
 654   // 4: charm and beauty
 655   // 5: beauty
 656   Bool_t res = 0;
 657
 658   if (fFlavor == 0) {
 659     res = kTRUE;
 660   } else {
 661     Int_t ifl = TMath::Abs(pid/100);
 662     if (ifl > 10) ifl/=10;
 663     res = (fFlavor == ifl);
 664   }
 665
 666   //  This part if gamma writing is inhibited
 667   if (fNoGammas)
 668     res = res && (pid != kGamma && pid != kPi0);
 669
 670   return res;
 671 }
 672
 673 Bool_t AliGenAmpt::Stable(TParticle* particle) const
 674 {
 675   // Return true for a stable particle
 676
 677   if (!particle)
 678     return kFALSE;
 679   if (particle->GetFirstDaughter() < 0 )
 680     return kTRUE;
 681   return kFALSE;
 682
 683   /// ADD LIST
 684
 685 }
 686
 687 void AliGenAmpt::MakeHeader()
 688 {
 689   // Fills the event header, to be called after each event
 690
 691   fHeader->SetNProduced(fNprimaries);
 692   fHeader->SetImpactParameter(fAmpt->GetHINT1(19));
 693   fHeader->SetTotalEnergy(fAmpt->GetEATT());
 694   fHeader->SetHardScatters(fAmpt->GetJATT());
 695   fHeader->SetParticipants(fAmpt->GetNP(), fAmpt->GetNT());
 696   fHeader->SetCollisions(fAmpt->GetN0(),
 697                         fAmpt->GetN01(),
 698                         fAmpt->GetN10(),
 699                         fAmpt->GetN11());
 700   fHeader->SetSpectators(fProjectileSpecn, fProjectileSpecp,
 701                         fTargetSpecn,fTargetSpecp);
 702   //fHeader->SetReactionPlaneAngle(fAmpt->GetHINT1(20));
 703   fHeader->SetReactionPlaneAngle(fAmpt->GetReactionPlaneAngle());
 704   //printf("Impact Parameter %13.3f \n", fAmpt->GetHINT1(19));
 705
 706   // 4-momentum vectors of the triggered jets.
 707   // Before final state gluon radiation.
 708   TLorentzVector* jet1 = new TLorentzVector(fAmpt->GetHINT1(21),
 709                                             fAmpt->GetHINT1(22),
 710                                             fAmpt->GetHINT1(23),
 711                                             fAmpt->GetHINT1(24));
 712
 713   TLorentzVector* jet2 = new TLorentzVector(fAmpt->GetHINT1(31),
 714                                             fAmpt->GetHINT1(32),
 715                                             fAmpt->GetHINT1(33),
 716                                             fAmpt->GetHINT1(34));
 717   // After final state gluon radiation.
 718   TLorentzVector* jet3 = new TLorentzVector(fAmpt->GetHINT1(26),
 719                                             fAmpt->GetHINT1(27),
 720                                             fAmpt->GetHINT1(28),
 721                                             fAmpt->GetHINT1(29));
 722
 723   TLorentzVector* jet4 = new TLorentzVector(fAmpt->GetHINT1(36),
 724                                             fAmpt->GetHINT1(37),
 725                                             fAmpt->GetHINT1(38),
 726                                             fAmpt->GetHINT1(39));
 727   fHeader->SetJets(jet1, jet2, jet3, jet4);
 728   // Bookkeeping for kinematic bias
 729   fHeader->SetTrials(fTrials);
 730   // Event Vertex
 731   fHeader->SetPrimaryVertex(fVertex);
 732   fHeader->SetInteractionTime(fEventTime);
 733
 734   fCollisionGeometry = fHeader;
 735   AddHeader(fHeader);
 736 }
 737
 738
 739 Bool_t AliGenAmpt::CheckTrigger()
 740 {
 741   // Check the kinematic trigger condition
 742
 743   Bool_t   triggered = kFALSE;
 744
 745   if (fTrigger == 1) {
 746     // jet-jet Trigger
 747     TLorentzVector* jet1 = new TLorentzVector(fAmpt->GetHINT1(26),
 748                                               fAmpt->GetHINT1(27),
 749                                               fAmpt->GetHINT1(28),
 750                                               fAmpt->GetHINT1(29));
 751
 752     TLorentzVector* jet2 = new TLorentzVector(fAmpt->GetHINT1(36),
 753                                               fAmpt->GetHINT1(37),
 754                                               fAmpt->GetHINT1(38),
 755                                               fAmpt->GetHINT1(39));
 756     Double_t eta1      = jet1->Eta();
 757     Double_t eta2      = jet2->Eta();
 758     Double_t phi1      = jet1->Phi();
 759     Double_t phi2      = jet2->Phi();
 760     //printf("\n Trigger: %f %f %f %f", fEtaMinJet, fEtaMaxJet, fPhiMinJet, fPhiMaxJet);
 761     if ( (eta1 < fEtaMaxJet && eta1 > fEtaMinJet &&
 762           phi1 < fPhiMaxJet && phi1 > fPhiMinJet)
 763          ||
 764          (eta2 < fEtaMaxJet && eta2 > fEtaMinJet &&
 765           phi2 < fPhiMaxJet && phi2 > fPhiMinJet)
 766       )
 767       triggered = kTRUE;
 768   } else if (fTrigger == 2) {
 769   // Gamma Jet
 770     Int_t np = fParticles.GetEntriesFast();
 771     for (Int_t i = 0; i < np; i++) {
 772       TParticle* part = (TParticle*) fParticles.At(i);
 773       Int_t kf = part->GetPdgCode();
 774       Int_t ksp = part->GetUniqueID();
 775       if (kf == 22 && ksp == 40) {
 776         Float_t phi = part->Phi();
 777         Float_t eta = part->Eta();
 778         if  (eta < fEtaMaxJet &&
 779              eta > fEtaMinJet &&
 780              phi < fPhiMaxJet &&
 781              phi > fPhiMinJet) {
 782           triggered = 1;
 783           break;
 784         } // check phi,eta within limits
 785       } // direct gamma ?
 786     } // particle loop
 787   } // fTrigger == 2
 788   return triggered;
 789 }