EVGEN/AliGenCorrHF.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 // Class to generate correlated Heavy Flavor hadron pairs (one or several pairs
  19 // per event) using paramtrized kinematics of quark pairs from some generator
  20 // and quark fragmentation functions.
  21 // Is a generalisation of AliGenParam class for correlated pairs of hadrons.
  22 // In this version quark pairs and fragmentation functions are obtained from
  23 // Pythia6.124 using 100K events generated with kCharmppMNRwmi&kBeautyppMNRwmi
  24 // in pp collisions at 14 TeV.
  25 // Decays are performed by Pythia. Used AliRoot version: v4-04-Release
  26 // Author: S. Grigoryan, LPC Clermont-Fd & YerPhI, Smbat.Grigoryan@cern.ch
  27 // July 07: added quarks in the stack (B. Vulpescu)
  28 //-------------------------------------------------------------------------
  29 // How it works (for the given flavor):
  30 //
  31 // 1) Reads QQbar kinematical grid from the Input file and generates
  32 // quark pairs according to the weights of the cells.
  33 // It is a 5D grid in y1,y2,pt1,pt2 and deltaphi, with occupancy weights
  34 // of the cells obtained from Pythia (see details in GetQuarkPair).
  35 // 2) Reads "soft" and "hard" fragmentation functions (12 2D-histograms each,
  36 // for 12 pt bins) from the Input file, applies to quarks and produces hadrons
  37 // (only lower states, with proportions of species obtained from Pythia).
  38 // Fragmentation functions are the same for all hadron species and depend
  39 // on 2 variables - light cone energy-momentum fractions:
  40 //     z1=(E_H + Pz_H)/(E_Q + Pz_Q),  z2=(E_H - Pz_H)/(E_Q - Pz_Q).
  41 // "soft" & "hard" FFs correspond to "slower" & "faster" quark of a pair
  42 // (see details in GetHadronPair).
  43 // 3) Decays the hadrons and saves all the particles in the event stack in the
  44 // following order: HF hadron from Q, then its decay products, then HF hadron
  45 // from Qbar, then its decay productes, then next HF hadon pair (if any)
  46 // in the same way, etc...
  47 // 4) It is fast, e.g., generates the same number of events with a beauty pair
  48 //  ~15 times faster than AliGenPythia with kBeautyppMNRwmi (w/o tracking)
  49 //
  50 // An Input file for each quark flavor is included in EVGEN/dataCorrHF/
  51 // One can use also user-defined Input files.
  52 //
  53 // More details could be found in my presentation at DiMuonNet Workshop, Dec 2006:
  54 // http://www-dapnia.cea.fr/Sphn/Alice/DiMuonNet
  55 // and will be published in an Internal Note.
  56 //
  57 //-------------------------------------------------------------------------
  58 // How to use it:
  59 //
  60 // add the following typical lines in Config.C
  61 /*
  62   if (!strcmp(option,"corr")) {
  63     // Example for correlated charm or beauty hadron pair production
  64
  65     // AliGenCorrHF *gener = new AliGenCorrHF(1, 4);  // for charm, 1 pair per event
  66     AliGenCorrHF *gener = new AliGenCorrHF(1, 5);  // for beauty, 1 pair per event
  67
  68     gener->SetMomentumRange(0,9999);
  69     gener->SetCutOnChild(0);        // 1/0 means cuts on children enable/disable
  70     gener->SetChildThetaRange(171.0,178.0);
  71     gener->SetOrigin(0,0,0);          //vertex position
  72     gener->SetSigma(0,0,0);           //Sigma in (X,Y,Z) (cm) on IP position
  73     gener->SetForceDecay(kSemiMuonic);
  74     gener->SetTrackingFlag(0);
  75     gener->Init();
  76 }
  77 */
  78 // and in aliroot do e.g. gAlice->Run(10,"Config.C") to produce 10 events.
  79 // One can include AliGenCorrHF in an AliGenCocktail generator.
  80 //--------------------------------------------------------------------------
  81
  82 #include <Riostream.h>
  83 #include <TCanvas.h>
  84 #include <TClonesArray.h>
  85 #include <TDatabasePDG.h>
  86 #include <TFile.h>
  87 #include <TH2F.h>
  88 #include <TLorentzVector.h>
  89 #include <TMath.h>
  90 #include <TParticle.h>
  91 #include <TParticlePDG.h>
  92 #include <TROOT.h>
  93 #include <TRandom.h>
  94 #include <TTree.h>
  95 #include <TVirtualMC.h>
  96 #include <TVector3.h>
  97
  98 #include "AliGenCorrHF.h"
  99 #include "AliLog.h"
 100 #include "AliConst.h"
 101 #include "AliDecayer.h"
 102 #include "AliMC.h"
 103 #include "AliRun.h"
 104
 105 ClassImp(AliGenCorrHF)
 106
 107   //Begin_Html
 108   /*
 109     <img src="picts/AliGenCorrHF.gif">
 110   */
 111   //End_Html
 112
 113 Double_t AliGenCorrHF::fgdph[19] = {0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, 180};
 114 Double_t AliGenCorrHF::fgy[31] = {-10,-7, -6.5, -6, -5.5, -5, -4.5, -4, -3.5, -3, -2.5, -2,- 1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 5.5, 6, 6.5, 7, 10};
 115 Double_t AliGenCorrHF::fgpt[33] = {0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 5.5, 6, 6.6, 7.2, 7.8, 8.4, 9, 9.6, 10.3, 11.1, 12, 13.1, 14.3, 15.6, 17.1, 19, 21, 24, 28, 35, 50, 100};
 116 Int_t AliGenCorrHF::fgnptbins = 12;
 117 Double_t AliGenCorrHF::fgptbmin[12] = {0, 0.5, 1, 1.5, 2, 2.5, 3, 4, 5, 6, 7, 9};
 118 Double_t AliGenCorrHF::fgptbmax[12] = {0.5, 1, 1.5, 2, 2.5, 3, 4, 5, 6, 7, 9, 100};
 119
 120 Double_t* AliGenCorrHF::fgIntegral = 0;
 121
 122 //____________________________________________________________
 123     AliGenCorrHF::AliGenCorrHF():
 124         fFileName(0),
 125         fFile(0),
 126         fQuark(0),
 127         fBias(0.),
 128         fTrials(0),
 129         fDecayer(0)
 130 {
 131 // Default constructor
 132 }
 133
 134 //____________________________________________________________
 135 AliGenCorrHF::AliGenCorrHF(Int_t npart, Int_t param):
 136     AliGenMC(npart),
 137     fFileName(0),
 138     fFile(0),
 139     fQuark(param),
 140     fBias(0.),
 141     fTrials(0),
 142     //    fDecayer(new AliDecayerPythia())
 143     fDecayer(0)
 144 {
 145 // Constructor using number of particles, quark type & default InputFile
 146 //
 147     if (fQuark != 5) fQuark = 4;
 148     fFileName = "$ALICE_ROOT/EVGEN/dataCorrHF/CharmppMNRwmiCorr100K.root";
 149     if (fQuark == 5) fFileName = "$ALICE_ROOT/EVGEN/dataCorrHF/BeautyppMNRwmiCorr100K.root";
 150
 151     fName = "Default";
 152     fTitle= "Generator for correlated pairs of HF hadrons";
 153
 154     fChildSelect.Set(5);
 155     for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
 156     SetForceDecay();
 157     SetCutOnChild();
 158     SetChildMomentumRange();
 159     SetChildPtRange();
 160     SetChildPhiRange();
 161     SetChildThetaRange();
 162 }
 163
 164 //___________________________________________________________________
 165 AliGenCorrHF::AliGenCorrHF(char* tname, Int_t npart, Int_t param):
 166     AliGenMC(npart),
 167     fFileName(tname),
 168     fFile(0),
 169     fQuark(param),
 170     fBias(0.),
 171     fTrials(0),
 172     //    fDecayer(new AliDecayerPythia())
 173     fDecayer(0)
 174 {
 175 // Constructor using number of particles, quark type & user-defined InputFile
 176 //
 177     if (fQuark != 5) fQuark = 4;
 178     fName = "UserDefined";
 179     fTitle= "Generator for correlated pairs of HF hadrons";
 180
 181     fChildSelect.Set(5);
 182     for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
 183     SetForceDecay();
 184     SetCutOnChild();
 185     SetChildMomentumRange();
 186     SetChildPtRange();
 187     SetChildPhiRange();
 188     SetChildThetaRange();
 189 }
 190
 191 //____________________________________________________________
 192 AliGenCorrHF::~AliGenCorrHF()
 193 {
 194 // Destructor
 195   delete fFile;
 196 }
 197
 198 //____________________________________________________________
 199 void AliGenCorrHF::Init()
 200 {
 201 // Initialisation
 202
 203   AliInfo(Form(" QQbar kinematics and fragm. functions from:  %s",fFileName.Data() ));
 204     fFile = TFile::Open(fFileName.Data());
 205     if(!fFile->IsOpen()){
 206       AliError(Form("Could not open file %s",fFileName.Data() ));
 207     }
 208
 209     ComputeIntegral(fFile);
 210
 211     fParentWeight = 1./fNpart;   // fNpart is number of HF-hadron pairs
 212
 213 // particle decay related initialization
 214
 215     if (gMC) fDecayer = gMC->GetDecayer();
 216     fDecayer->SetForceDecay(fForceDecay);
 217     fDecayer->Init();
 218
 219 //
 220     AliGenMC::Init();
 221 }
 222
 223 //____________________________________________________________
 224 void AliGenCorrHF::Generate()
 225 {
 226 //
 227 // Generate fNpart of correlated HF hadron pairs per event
 228 // in the the desired theta and momentum windows (phi = 0 - 2pi).
 229 // Gaussian smearing on the vertex is done if selected.
 230 // The decay of heavy hadrons is done using lujet,
 231 //    and the childern particle are tracked by GEANT
 232 // However, light mesons are directly tracked by GEANT
 233 // setting fForceDecay = nodecay (SetForceDecay(nodecay))
 234 //
 235
 236
 237   Float_t polar[3]= {0,0,0};  // Polarisation of the parent particle (for GEANT tracking)
 238   Float_t origin0[3];         // Origin of the generated parent particle (for GEANT tracking)
 239   Float_t pt, pl, ptot;       // Transverse, logitudinal and total momenta of the parent particle
 240   Float_t phi, theta;         // Phi and theta spherical angles of the parent particle momentum
 241   Float_t p[3], pc[3], och[3];// Momentum, polarisation and origin of the children particles from lujet
 242
 243   Float_t pq1[3], pq2[3];     // Momenta of the two quarks
 244
 245   Double_t dphi=0, ptq[2], yq[2], pth[2], plh[2], ph[2], phih[2], phiq[2];
 246   Int_t i, j, ipair, ihadron[2];
 247   for (i=0; i<2; i++) {
 248     ptq[i]     =0;
 249     yq[i]      =0;
 250     pth[i]     =0;
 251     plh[i]     =0;
 252     phih[i]    =0;
 253     phiq[i]    =0;
 254     ihadron[i] =0;
 255   }
 256
 257   static TClonesArray *particles;
 258   //
 259   if(!particles) particles = new TClonesArray("TParticle",1000);
 260
 261   TDatabasePDG* pDataBase = TDatabasePDG::Instance();
 262
 263 // Calculating vertex position per event
 264   for (j=0;j<3;j++) origin0[j]=fOrigin[j];
 265   if(fVertexSmear==kPerEvent) {
 266       Vertex();
 267       for (j=0;j<3;j++) origin0[j]=fVertex[j];
 268   }
 269
 270   Float_t wgtp, wgtch, random[6];
 271   Int_t ipap = 0;
 272   Int_t nt   = 0;
 273   Int_t ntq1 = 0;
 274   Int_t ntq2 = 0;
 275
 276 // Generating fNpart HF-hadron pairs per event
 277   while(ipap<fNpart) {
 278
 279     while(1) {
 280
 281       GetQuarkPair(fFile, fgIntegral, yq[0], yq[1], ptq[0], ptq[1], dphi);
 282
 283 // Add the quarks in the stack
 284
 285       phiq[0] = Rndm()*k2PI;
 286       phiq[1] = phiq[0] + dphi*kDegrad;
 287       TVector3 qvect1 = TVector3();
 288       TVector3 qvect2 = TVector3();
 289       qvect1.SetPtEtaPhi(ptq[0],yq[0],phiq[0]);
 290       qvect2.SetPtEtaPhi(ptq[1],yq[1],phiq[1]);
 291       pq1[0] = qvect1.Px();
 292       pq1[1] = qvect1.Py();
 293       pq1[2] = qvect1.Pz();
 294       pq2[0] = qvect2.Px();
 295       pq2[1] = qvect2.Py();
 296       pq2[2] = qvect2.Pz();
 297
 298       wgtp  = fParentWeight;
 299
 300       PushTrack(0, -1, +fQuark, pq1, origin0, polar, 0,
 301                 kPPrimary, nt, wgtp);
 302       KeepTrack(nt);
 303       ntq1 = nt;
 304
 305       PushTrack(0, -1, -fQuark, pq2, origin0, polar, 0,
 306                 kPPrimary, nt, wgtp);
 307       KeepTrack(nt);
 308       ntq2 = nt;
 309
 310       nt = 2;
 311
 312       GetHadronPair(fFile, fQuark, yq[0], yq[1], ptq[0], ptq[1], ihadron[0], ihadron[1], plh[0], plh[1], pth[0], pth[1]);
 313
 314 // Here we assume that  |phi_H1 - phi_H2| = |phi_Q1 - phi_Q2| = dphi
 315 // which is a good approximation for heavy flavors in Pythia
 316
 317       /* // doesn't work if PhiMax < k2PI or PhiMin > 0, since dphi = 0 - 180
 318       phih[0] = fPhiMin + Rndm()*(fPhiMax-fPhiMin);
 319       phih[1] = phih[0] + dphi*kDegrad;
 320       if (phih[0] > fPhiMax/2.) phih[1] = phih[0] - dphi*kDegrad;
 321       */
 322       phih[0] = Rndm()*k2PI;
 323       phih[1] = phih[0] + dphi*kDegrad;
 324       if (phih[0] > TMath::Pi()) phih[1] = phih[0] - dphi*kDegrad;
 325
 326 // Cut on theta
 327       theta=TMath::ATan2(pth[0],plh[0]);
 328       if(theta<fThetaMin || theta>fThetaMax) continue;
 329       theta=TMath::ATan2(pth[1],plh[1]);
 330       if(theta<fThetaMin || theta>fThetaMax) continue;
 331
 332 // Cut on momentum
 333       ph[0]=TMath::Sqrt(pth[0]*pth[0]+plh[0]*plh[0]);
 334       if (ph[0]<fPMin || ph[0]>fPMax) continue;
 335       ph[1]=TMath::Sqrt(pth[1]*pth[1]+plh[1]*plh[1]);
 336       if (ph[1]<fPMin || ph[1]>fPMax) continue;
 337
 338 // Common origin for particles of the HF-hadron pair
 339       if(fVertexSmear==kPerTrack) {
 340           Rndm(random,6);
 341           for (j=0;j<3;j++) {
 342               origin0[j]=
 343                 fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
 344                 TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 345           }
 346       }
 347
 348       Int_t np1=0, kf1[100], select1[100], iparent1[100], trackIt1[100];
 349       Float_t  wgtch1=0, p1[3], pc1[100][3], och1[100][3];
 350
 351       for (j=0; j<3; j++)   p1[j] = 0;
 352       for (i=0; i<100; i++) {
 353         kf1[i]      = 0;
 354         select1[i]  = 0;
 355         iparent1[i] = 0;
 356         trackIt1[i] = 0;
 357         for (j=0; j<3; j++) {
 358           pc1[i][j]  = 0;
 359           och1[i][j] = 0;
 360         }
 361       }
 362
 363 //
 364 // Loop over particles of the HF-hadron pair
 365       Int_t nhadron = 0;
 366       for (ipair=0;ipair<2;ipair++) {
 367         phi  = phih[ipair];
 368         pl   = plh[ipair];
 369         pt   = pth[ipair];
 370         ptot = ph[ipair];
 371 //
 372 // particle type
 373           Int_t iPart = ihadron[ipair];
 374           Float_t am  = pDataBase->GetParticle(iPart)->Mass();
 375           fChildWeight=(fDecayer->GetPartialBranchingRatio(iPart))*fParentWeight;
 376
 377           wgtp  = fParentWeight;
 378           wgtch = fChildWeight;
 379
 380 //
 381           p[0]=pt*TMath::Cos(phi);
 382           p[1]=pt*TMath::Sin(phi);
 383           p[2]=pl;
 384
 385 // Looking at fForceDecay :
 386 // if fForceDecay != none Primary particle decays using
 387 // AliPythia and children are tracked by GEANT
 388 //
 389 // if fForceDecay == none Primary particle is tracked by GEANT
 390 // (In the latest, make sure that GEANT actually does all the decays you want)
 391 //
 392
 393           if (fForceDecay != kNoDecay) {
 394 // Using lujet to decay particle
 395               Float_t energy=TMath::Sqrt(ptot*ptot+am*am);
 396               TLorentzVector pmom(p[0], p[1], p[2], energy);
 397               fDecayer->Decay(iPart,&pmom);
 398 //
 399 // select decay particles
 400               Int_t np=fDecayer->ImportParticles(particles);
 401
 402 //  Selecting  GeometryAcceptance for particles fPdgCodeParticleforAcceptanceCut;
 403               if (fGeometryAcceptance)
 404                 if (!CheckAcceptanceGeometry(np,particles)) break;
 405               Int_t ncsel=0;
 406               Int_t* pParent    = new Int_t[np];
 407               Int_t* pSelected  = new Int_t[np];
 408               Int_t* trackIt    = new Int_t[np];
 409
 410               for (i=0; i<np; i++) {
 411                   pSelected[i] =  0;
 412                   pParent[i]   = -1;
 413               }
 414
 415               if (np >1) {
 416                   TParticle* iparticle =  (TParticle *) particles->At(0);
 417                   for (i=1; i<np; i++) {
 418                       trackIt[i] = 1;
 419                       iparticle = (TParticle *) particles->At(i);
 420                       Int_t kf = iparticle->GetPdgCode();
 421                       Int_t ks = iparticle->GetStatusCode();
 422
 423 // particles with long life-time (c tau > .3 mum)
 424                       if (ks != 1) {
 425                           Double_t lifeTime = fDecayer->GetLifetime(kf);
 426                           if (lifeTime <= (Double_t) fMaxLifeTime) {
 427                               trackIt[i]     = 0;
 428                               pSelected[i]   = 1;
 429                           }
 430                       } // ks==1 ?
 431 //
 432 // children, discard neutrinos
 433                       if (TMath::Abs(kf) == 12 || TMath::Abs(kf) == 14) continue;
 434                       if (trackIt[i])
 435                       {
 436                           if (fCutOnChild) {
 437                               pc[0]=iparticle->Px();
 438                               pc[1]=iparticle->Py();
 439                               pc[2]=iparticle->Pz();
 440                               Bool_t  childok = KinematicSelection(iparticle, 1);
 441                               if(childok) {
 442                                   pSelected[i]  = 1;
 443                                   ncsel++;
 444                               } else {
 445                                   ncsel=-1;
 446                                   break;
 447                               } // child kine cuts
 448                           } else {
 449                               pSelected[i]  = 1;
 450                               ncsel++;
 451                           } // if child selection
 452                       } // select muon
 453                   } // decay particle loop
 454               } // if decay products
 455
 456               Int_t iparent;
 457               if ((fCutOnChild && ncsel >0) || !fCutOnChild){
 458
 459                   nhadron++;
 460 //
 461 // Parents and Decay Products
 462                   if (ipair == 0) {
 463                       np1    = np;
 464                       wgtch1 = wgtch;
 465                       p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2];
 466                   } else {
 467                       ipap++;
 468                       PushTrack(0, ntq1, ihadron[0], p1, origin0, polar, 0,
 469                                 kPPrimary, nt, wgtp);
 470                       KeepTrack(nt);
 471                       for (i = 1; i < np1; i++) {
 472                         if (select1[i]) {
 473                           for (j=0; j<3; j++) {
 474                               och[j] = och1[i][j];
 475                               pc[j]  = pc1[i][j];
 476                           }
 477                           PushTrack(fTrackIt*trackIt1[i], iparent1[i], kf1[i], pc, och,
 478                                     polar, 0, kPDecay, nt, wgtch1);
 479                           KeepTrack(nt);
 480                           }
 481                       }
 482                       PushTrack(0, ntq2, iPart, p, origin0, polar, 0, kPPrimary, nt, wgtp);
 483                       KeepTrack(nt);
 484                   }
 485                   pParent[0] = nt;
 486 //
 487 // Decay Products
 488                   Int_t ntcount = 0;
 489                   for (i = 1; i < np; i++) {
 490                     if (pSelected[i]) {
 491                           TParticle* iparticle = (TParticle *) particles->At(i);
 492                           Int_t kf  = iparticle->GetPdgCode();
 493                           Int_t ipa = iparticle->GetFirstMother()-1;
 494
 495                           och[0] = origin0[0]+iparticle->Vx()/10;
 496                           och[1] = origin0[1]+iparticle->Vy()/10;
 497                           och[2] = origin0[2]+iparticle->Vz()/10;
 498                           pc[0]  = iparticle->Px();
 499                           pc[1]  = iparticle->Py();
 500                           pc[2]  = iparticle->Pz();
 501
 502                           if (ipa > -1) {
 503                               iparent = pParent[ipa];
 504                           } else {
 505                               iparent = -1;
 506                           }
 507
 508                           if (ipair == 0) {
 509                               kf1[i]      = kf;
 510                               select1[i]  = pSelected[i];
 511                               iparent1[i] = iparent;
 512                               trackIt1[i] = trackIt[i];
 513                               for (j=0; j<3; j++) {
 514                                   och1[i][j] = och[j];
 515                                   pc1[i][j]  = pc[j];
 516                               }
 517                               ntcount++;
 518                           } else {
 519                             PushTrack(fTrackIt*trackIt[i], iparent, kf, pc, och,
 520                                       polar, 0, kPDecay, nt, wgtch);
 521                             KeepTrack(nt);
 522                           }
 523                           pParent[i] = nt + ntcount;
 524                     } // Selected
 525                   } // Particle loop
 526               }  // Decays by Lujet
 527               particles->Clear();
 528               if (pParent)    delete[] pParent;
 529               if (pSelected)  delete[] pSelected;
 530               if (trackIt)    delete[] trackIt;
 531           } // kinematic selection
 532           else  // nodecay option, so parent will be tracked by GEANT
 533           {
 534             nhadron++;
 535             if (ipair == 0) {
 536                 p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2];
 537             } else {
 538                 ipap++;
 539                 gAlice->GetMCApp()->
 540                 PushTrack(fTrackIt,ntq1,ihadron[0],p1,origin0,polar,0,kPPrimary,nt,wgtp);
 541                 gAlice->GetMCApp()->
 542                 PushTrack(fTrackIt,ntq2,iPart,p,origin0,polar,0,kPPrimary,nt,wgtp);
 543             }
 544           }
 545           if (nhadron == 0) break;
 546       } // ipair loop
 547       if (nhadron != 2) continue;
 548       break;
 549     } // while(1)
 550     nt++;
 551   } // while(ipa<fNpart) --> event loop
 552
 553   SetHighWaterMark(nt);
 554 }
 555
 556 //____________________________________________________________________________________
 557 Int_t AliGenCorrHF::IpCharm(TRandom* ran)
 558 {
 559 // Composition of lower state charm hadrons, containing a c-quark
 560     Float_t random;
 561     Int_t ip;            // +- 411,421,431,4122,4132,4232,4332
 562     random = ran->Rndm();
 563 //  Rates from Pythia6.214 using 100Kevents with kPyCharmppMNRwmi at 14 TeV.
 564
 565     if (random < 0.6027) {
 566         ip=421;
 567     } else if (random < 0.7962) {
 568         ip=411;
 569     } else if (random < 0.9127) {
 570         ip=431;
 571     } else if (random < 0.9899) {
 572         ip=4122;
 573     } else if (random < 0.9948) {
 574         ip=4132;
 575     } else if (random < 0.9999) {
 576         ip=4232;
 577     } else {
 578         ip=4332;
 579     }
 580
 581     return ip;
 582 }
 583
 584 Int_t AliGenCorrHF::IpBeauty(TRandom* ran)
 585 {
 586 // Composition of lower state beauty hadrons, containing a b-quark
 587     Float_t random;
 588     Int_t ip;            // +- 511,521,531,5122,5132,5232,5332
 589     random = ran->Rndm();
 590 //  Rates from Pythia6.214 using 100Kevents with kPyBeautyppMNRwmi at 14 TeV.
 591                         // B-Bbar mixing will be done by Pythia at the decay point
 592  if (random < 0.3965) {
 593         ip=-511;
 594     } else if (random < 0.7930) {
 595         ip=-521;
 596     } else if (random < 0.9112) {
 597         ip=-531;
 598     } else if (random < 0.9887) {
 599         ip=5122;
 600     } else if (random < 0.9943) {
 601         ip=5132;
 602     } else if (random < 0.9999) {
 603         ip=5232;
 604     } else {
 605         ip=5332;
 606     }
 607
 608   return ip;
 609 }
 610
 611 //____________________________________________________________________________________
 612 Double_t AliGenCorrHF::ComputeIntegral(TFile* fG)       // needed by GetQuarkPair
 613 {
 614    // Read QQbar kinematical 5D grid's cell occupancy weights
 615    Int_t* cell = new Int_t[6];           // cell[6]={wght,iy1,iy2,ipt1,ipt2,idph}
 616    TTree* tG = (TTree*) fG->Get("tGqq");
 617    tG->GetBranch("cell")->SetAddress(cell);
 618    Int_t nbins = tG->GetEntries();
 619
 620    //   delete previously computed integral (if any)
 621    if(fgIntegral) delete [] fgIntegral;
 622
 623    fgIntegral = new Double_t[nbins+1];
 624    fgIntegral[0] = 0;
 625    Int_t bin;
 626    for(bin=0;bin<nbins;bin++) {
 627      tG->GetEvent(bin);
 628      fgIntegral[bin+1] = fgIntegral[bin] + cell[0];
 629    }
 630    //   Normalize integral to 1
 631    if (fgIntegral[nbins] == 0 ) {
 632       return 0;
 633    }
 634    for (bin=1;bin<=nbins;bin++)  fgIntegral[bin] /= fgIntegral[nbins];
 635
 636    return fgIntegral[nbins];
 637 }
 638
 639 //____________________________________________________________________________________
 640 void AliGenCorrHF::GetQuarkPair(TFile* fG, Double_t* fInt, Double_t &y1, Double_t &y2, Double_t &pt1, Double_t &pt2, Double_t &dphi)
 641                                  // modification of ROOT's TH3::GetRandom3 for 5D
 642 {
 643    // Read QQbar kinematical 5D grid's cell coordinates
 644    Int_t* cell = new Int_t[6];           // cell[6]={wght,iy1,iy2,ipt1,ipt2,idph}
 645    TTree* tG = (TTree*) fG->Get("tGqq");
 646    tG->GetBranch("cell")->SetAddress(cell);
 647    Int_t nbins = tG->GetEntries();
 648    Double_t rand[6];
 649    gRandom->RndmArray(6,rand);
 650    Int_t ibin = TMath::BinarySearch(nbins,fInt,rand[0]);
 651    tG->GetEvent(ibin);
 652    y1   = fgy[cell[1]]  + (fgy[cell[1]+1]-fgy[cell[1]])*rand[1];
 653    y2   = fgy[cell[2]]  + (fgy[cell[2]+1]-fgy[cell[2]])*rand[2];
 654    pt1  = fgpt[cell[3]] + (fgpt[cell[3]+1]-fgpt[cell[3]])*rand[3];
 655    pt2  = fgpt[cell[4]] + (fgpt[cell[4]+1]-fgpt[cell[4]])*rand[4];
 656    dphi = fgdph[cell[5]]+ (fgdph[cell[5]+1]-fgdph[cell[5]])*rand[5];
 657 }
 658
 659 //____________________________________________________________________________________
 660 void AliGenCorrHF::GetHadronPair(TFile* fG, Int_t idq, Double_t y1, Double_t y2, Double_t pt1, Double_t pt2, Int_t &id3, Int_t &id4, Double_t &pz3, Double_t &pz4, Double_t &pt3, Double_t &pt4)
 661 {
 662     // Generate a hadron pair
 663     Int_t (*fIpParaFunc )(TRandom*);//Pointer to particle type parametrisation function
 664     fIpParaFunc = IpCharm;
 665     Double_t mq = 1.2;              // c & b quark masses (used in AliPythia)
 666     if (idq == 5) {
 667       fIpParaFunc = IpBeauty;
 668       mq = 4.75;
 669     }
 670     Double_t z11, z12, z21, z22, pz1, pz2, e1, e2, mh, ptemp, rand[2];
 671     char tag[100];
 672     TH2F *h2h[12], *h2s[12];      // hard & soft Fragmentation Functions
 673     for (Int_t ipt = 0; ipt<fgnptbins; ipt++) {
 674       sprintf(tag,"h2h_pt%d",ipt);
 675       h2h[ipt] = (TH2F*) fG->Get(tag);
 676       sprintf(tag,"h2s_pt%d",ipt);
 677       h2s[ipt] = (TH2F*) fG->Get(tag);
 678     }
 679
 680        if (y1*y2 < 0) {
 681          for (Int_t ipt = 0; ipt<fgnptbins; ipt++) {
 682            if(pt1 >= fgptbmin[ipt] && pt1 < fgptbmax[ipt])
 683              h2h[ipt]->GetRandom2(z11, z21);
 684            if(pt2 >= fgptbmin[ipt] && pt2 < fgptbmax[ipt])
 685              h2h[ipt]->GetRandom2(z12, z22);
 686          }
 687        }
 688        else {
 689          if (TMath::Abs(y1) > TMath::Abs(y2)) {
 690            for (Int_t ipt = 0; ipt<fgnptbins; ipt++) {
 691              if(pt1 >= fgptbmin[ipt] && pt1 < fgptbmax[ipt])
 692                h2h[ipt]->GetRandom2(z11, z21);
 693              if(pt2 >= fgptbmin[ipt] && pt2 < fgptbmax[ipt])
 694                h2s[ipt]->GetRandom2(z12, z22);
 695            }
 696          }
 697          else {
 698            for (Int_t ipt = 0; ipt<fgnptbins; ipt++) {
 699              if(pt1 >= fgptbmin[ipt] && pt1 < fgptbmax[ipt])
 700                h2s[ipt]->GetRandom2(z11, z21);
 701              if(pt2 >= fgptbmin[ipt] && pt2 < fgptbmax[ipt])
 702                h2h[ipt]->GetRandom2(z12, z22);
 703            }
 704          }
 705        }
 706       gRandom->RndmArray(2,rand);
 707       ptemp = TMath::Sqrt(pt1*pt1 + mq*mq);
 708       pz1   = ptemp*TMath::SinH(y1);
 709       e1    = ptemp*TMath::CosH(y1);
 710       ptemp = TMath::Sqrt(pt2*pt2 + mq*mq);
 711       pz2   = ptemp*TMath::SinH(y2);
 712       e2    = ptemp*TMath::CosH(y2);
 713
 714       id3   = fIpParaFunc(gRandom);
 715       mh    = TDatabasePDG::Instance()->GetParticle(id3)->Mass();
 716       ptemp = z11*z21*(e1*e1-pz1*pz1) - mh*mh;
 717       pt3   = (idq-3)*rand[0];                // some smearing at low pt, try better
 718       if (ptemp > 0) pt3 = TMath::Sqrt(ptemp);
 719       if (pz1 > 0)   pz3 = (z11*(e1 + pz1) - z21*(e1 - pz1)) / 2;
 720       else           pz3 = (z21*(e1 + pz1) - z11*(e1 - pz1)) / 2;
 721       e1 = TMath::Sqrt(pz3*pz3 + pt3*pt3 + mh*mh);
 722
 723       id4   = - fIpParaFunc(gRandom);
 724       mh    = TDatabasePDG::Instance()->GetParticle(id4)->Mass();
 725       ptemp = z12*z22*(e2*e2-pz2*pz2) - mh*mh;
 726       pt4   = (idq-3)*rand[1];                // some smearing at low pt, try better
 727       if (ptemp > 0) pt4 = TMath::Sqrt(ptemp);
 728       if (pz2 > 0)   pz4 = (z12*(e2 + pz2) - z22*(e2 - pz2)) / 2;
 729       else           pz4 = (z22*(e2 + pz2) - z12*(e2 - pz2)) / 2;
 730       e2 = TMath::Sqrt(pz4*pz4 + pt4*pt4 + mh*mh);
 731
 732       // small corr. instead of using Frag. Func. depending on yQ (in addition to ptQ)
 733       Float_t ycorr = 0.2, y3, y4;
 734       gRandom->RndmArray(2,rand);
 735       y3 = 0.5 * TMath::Log((e1 + pz3 + 1.e-13)/(e1 - pz3 + 1.e-13));
 736       y4 = 0.5 * TMath::Log((e2 + pz4 + 1.e-13)/(e2 - pz4 + 1.e-13));
 737       if(TMath::Abs(y3)<ycorr && TMath::Abs(y4)<ycorr && rand[0]>0.5) {
 738         ptemp = TMath::Sqrt(e1*e1 - pz3*pz3);
 739         y3  = 4*(1 - 2*rand[1]);
 740         pz3 = ptemp*TMath::SinH(y3);
 741         pz4 = pz3;
 742       }
 743 }