Updates (N. Bastid)
[u/mrichter/AliRoot.git] / 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 }