Updates (N. Bastid)

diff --git a/EVGEN/AliGenCorrHF.cxx b/EVGEN/AliGenCorrHF.cxx
index 5dfdcee33dc151a3eeb5bc2354e12dbb82a300f2..41e7c3cfe8833bf8dad3deb1bdb309b4864bb0ca 100644 (file)
--- a/EVGEN/AliGenCorrHF.cxx
+++ b/EVGEN/AliGenCorrHF.cxx
@@ -20,11 +20,11 @@
 // and quark fragmentation functions.
 // Is a generalisation of AliGenParam class for correlated pairs of hadrons.
 // In this version quark pairs and fragmentation functions are obtained from
 // and quark fragmentation functions.
 // Is a generalisation of AliGenParam class for correlated pairs of hadrons.
 // In this version quark pairs and fragmentation functions are obtained from

-// Pythia6.124 using 100K events generated with kCharmppMNRwmi & kBeautyppMNRwmi 
+// Pythia6.124 using 100K events generated with kCharmppMNRwmi&kBeautyppMNRwmi 

 // in pp collisions at 14 TeV.
 // Decays are performed by Pythia. Used AliRoot version: v4-04-Release
 // Author: S. Grigoryan, LPC Clermont-Fd & YerPhI, Smbat.Grigoryan@cern.ch
 // in pp collisions at 14 TeV.
 // Decays are performed by Pythia. Used AliRoot version: v4-04-Release
 // Author: S. Grigoryan, LPC Clermont-Fd & YerPhI, Smbat.Grigoryan@cern.ch

-//
+// July 07: added quarks in the stack (B. Vulpescu)

 //-------------------------------------------------------------------------
 // How it works (for the given flavor):
 //
 //-------------------------------------------------------------------------
 // How it works (for the given flavor):
 //


@@ -93,6 +93,7 @@
 #include <TRandom.h>
 #include <TTree.h>
 #include <TVirtualMC.h>
 #include <TRandom.h>
 #include <TTree.h>
 #include <TVirtualMC.h>

+#include <TVector3.h>

 
 #include "AliGenCorrHF.h"
 #include "AliLog.h"
 
 #include "AliGenCorrHF.h"
 #include "AliLog.h"


@@ -206,7 +207,7 @@ void AliGenCorrHF::Init()
     }
 
     ComputeIntegral(fFile);
     }
 
     ComputeIntegral(fFile);

-
+    

     fParentWeight = 1./fNpart;   // fNpart is number of HF-hadron pairs
 
 // particle decay related initialization
     fParentWeight = 1./fNpart;   // fNpart is number of HF-hadron pairs
 
 // particle decay related initialization


@@ -239,14 +240,17 @@ void AliGenCorrHF::Generate()
   Float_t phi, theta;         // Phi and theta spherical angles of the parent particle momentum
   Float_t p[3], pc[3], och[3];// Momentum, polarisation and origin of the children particles from lujet
 
   Float_t phi, theta;         // Phi and theta spherical angles of the parent particle momentum
   Float_t p[3], pc[3], och[3];// Momentum, polarisation and origin of the children particles from lujet
 

+  Float_t pq1[3], pq2[3];     // Momenta of the two quarks

 
 

-  Double_t dphi=0, ptq[2], yq[2], pth[2], plh[2], ph[2], phih[2];
+  Double_t dphi=0, ptq[2], yq[2], pth[2], plh[2], ph[2], phih[2], phiq[2];

   Int_t i, j, ipair, ihadron[2];
   for (i=0; i<2; i++) { 
   Int_t i, j, ipair, ihadron[2];
   for (i=0; i<2; i++) { 

-    ptq[i] =0; 
-    yq[i]  =0; 
-    pth[i] =0; 
-    plh[i] =0; 
+    ptq[i]     =0; 
+    yq[i]      =0; 
+    pth[i]     =0; 
+    plh[i]     =0;
+    phih[i]    =0; 
+    phiq[i]    =0;

     ihadron[i] =0; 
   }
 
     ihadron[i] =0; 
   }
 


@@ -266,6 +270,8 @@ void AliGenCorrHF::Generate()
   Float_t wgtp, wgtch, random[6];
   Int_t ipap = 0;
   Int_t nt   = 0;
   Float_t wgtp, wgtch, random[6];
   Int_t ipap = 0;
   Int_t nt   = 0;

+  Int_t ntq1 = 0;
+  Int_t ntq2 = 0;

 
 // Generating fNpart HF-hadron pairs per event
   while(ipap<fNpart) {
 
 // Generating fNpart HF-hadron pairs per event
   while(ipap<fNpart) {


@@ -274,6 +280,35 @@ void AliGenCorrHF::Generate()
 
       GetQuarkPair(fFile, fgIntegral, yq[0], yq[1], ptq[0], ptq[1], dphi);
 
 
       GetQuarkPair(fFile, fgIntegral, yq[0], yq[1], ptq[0], ptq[1], dphi);
 

+// Add the quarks in the stack
+
+      phiq[0] = Rndm()*k2PI;
+      phiq[1] = phiq[0] + dphi*kDegrad; 
+      TVector3 qvect1 = TVector3();
+      TVector3 qvect2 = TVector3();
+      qvect1.SetPtEtaPhi(ptq[0],yq[0],phiq[0]);
+      qvect2.SetPtEtaPhi(ptq[1],yq[1],phiq[1]);
+      pq1[0] = qvect1.Px();
+      pq1[1] = qvect1.Py();
+      pq1[2] = qvect1.Pz();
+      pq2[0] = qvect2.Px();
+      pq2[1] = qvect2.Py();
+      pq2[2] = qvect2.Pz();
+      
+      wgtp  = fParentWeight;
+
+      PushTrack(0, -1, +fQuark, pq1, origin0, polar, 0, 
+               kPPrimary, nt, wgtp);
+      KeepTrack(nt);
+      ntq1 = nt;
+
+      PushTrack(0, -1, -fQuark, pq2, origin0, polar, 0, 
+               kPPrimary, nt, wgtp);
+      KeepTrack(nt);
+      ntq2 = nt;
+
+      nt = 2;
+      

       GetHadronPair(fFile, fQuark, yq[0], yq[1], ptq[0], ptq[1], ihadron[0], ihadron[1], plh[0], plh[1], pth[0], pth[1]);
 
 // Here we assume that  |phi_H1 - phi_H2| = |phi_Q1 - phi_Q2| = dphi
       GetHadronPair(fFile, fQuark, yq[0], yq[1], ptq[0], ptq[1], ihadron[0], ihadron[1], plh[0], plh[1], pth[0], pth[1]);
 
 // Here we assume that  |phi_H1 - phi_H2| = |phi_Q1 - phi_Q2| = dphi


@@ -430,7 +465,7 @@ void AliGenCorrHF::Generate()
                      p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2];
                  } else {
                      ipap++;
                      p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2];
                  } else {
                      ipap++;

-                     PushTrack(0, -1, ihadron[0], p1, origin0, polar, 0, 
+                     PushTrack(0, ntq1, ihadron[0], p1, origin0, polar, 0, 

                                kPPrimary, nt, wgtp);
                      KeepTrack(nt);
                      for (i = 1; i < np1; i++) {
                                kPPrimary, nt, wgtp);
                      KeepTrack(nt);
                      for (i = 1; i < np1; i++) {


@@ -444,7 +479,7 @@ void AliGenCorrHF::Generate()
                          KeepTrack(nt);
                          }
                      }
                          KeepTrack(nt);
                          }
                      }

-                     PushTrack(0, -1, iPart, p, origin0, polar, 0, kPPrimary, nt, wgtp);
+                     PushTrack(0, ntq2, iPart, p, origin0, polar, 0, kPPrimary, nt, wgtp);

                      KeepTrack(nt);
                  } 
                  pParent[0] = nt;
                      KeepTrack(nt);
                  } 
                  pParent[0] = nt;


@@ -502,9 +537,9 @@ void AliGenCorrHF::Generate()
            } else {
                ipap++;
                gAlice->GetMCApp()->
            } else {
                ipap++;
                gAlice->GetMCApp()->

-               PushTrack(fTrackIt,-1,ihadron[0],p1,origin0,polar,0,kPPrimary,nt,wgtp);
+               PushTrack(fTrackIt,ntq1,ihadron[0],p1,origin0,polar,0,kPPrimary,nt,wgtp);

                gAlice->GetMCApp()->
                gAlice->GetMCApp()->

-               PushTrack(fTrackIt,-1,iPart,p,origin0,polar,0,kPPrimary,nt,wgtp);
+               PushTrack(fTrackIt,ntq2,iPart,p,origin0,polar,0,kPPrimary,nt,wgtp);

            }
          }
          if (nhadron == 0) break;
            }
          }
          if (nhadron == 0) break;

diff --git a/EVGEN/AliGenMUONCocktailpp.cxx b/EVGEN/AliGenMUONCocktailpp.cxx
index e96d1e4a69c2544d4b6f7211218008a015b36d4e..b8f95538fc60eb668b30cafe5ef5640c4aa0ffac 100644 (file)
--- a/EVGEN/AliGenMUONCocktailpp.cxx
+++ b/EVGEN/AliGenMUONCocktailpp.cxx
@@ -22,7 +22,8 @@
 // The free parameeters are :
 //      pp reaction cross-section
 //      production cross-sections in pp collisions 
 // The free parameeters are :
 //      pp reaction cross-section
 //      production cross-sections in pp collisions 

-// 
+// July 07:added heavy quark production from AliGenCorrHF and heavy quark 
+//         production switched off in Pythia 

 
 #include <TObjArray.h>
 #include <TParticle.h>
 
 #include <TObjArray.h>
 #include <TParticle.h>


@@ -41,7 +42,7 @@
 #include "AliDecayer.h"
 #include "AliLog.h"
 #include "AliGenPythia.h"
 #include "AliDecayer.h"
 #include "AliLog.h"
 #include "AliGenPythia.h"

-
+#include "AliGenCorrHF.h"

 
 ClassImp(AliGenMUONCocktailpp)  
   
 
 ClassImp(AliGenMUONCocktailpp)  
   


@@ -93,7 +94,9 @@ void AliGenMUONCocktailpp::Init()
     Double_t ratioupsilon;
     Double_t ratioupsilonP;
     Double_t ratioupsilonPP;
     Double_t ratioupsilon;
     Double_t ratioupsilonP;
     Double_t ratioupsilonPP;

-    
+    Double_t ratioccbar;
+    Double_t ratiobbbar;
+

 // Cross sections in barns (from PPR Vol. II p: 552) pp - 14 TeV and 
 // corrected from feed down of higher resonances 
 
 // Cross sections in barns (from PPR Vol. II p: 552) pp - 14 TeV and 
 // corrected from feed down of higher resonances 
 


@@ -102,6 +105,8 @@ void AliGenMUONCocktailpp::Init()
     Double_t sigmaupsilon = 0.989e-6;  
     Double_t sigmaupsilonP = 0.502e-6;  
     Double_t sigmaupsilonPP = 0.228e-6;
     Double_t sigmaupsilon = 0.989e-6;  
     Double_t sigmaupsilonP = 0.502e-6;  
     Double_t sigmaupsilonPP = 0.228e-6;

+    Double_t sigmaccbar = 11.2e-3;
+    Double_t sigmabbbar = 0.51e-3;

     
     AliInfo(Form("the parametrised resonances uses the decay mode %d",fDecayModeResonance));
 
     
     AliInfo(Form("the parametrised resonances uses the decay mode %d",fDecayModeResonance));
 


@@ -212,7 +217,30 @@ void AliGenMUONCocktailpp::Init()
     genupsilonPP->Init(); // generation in selected kinematical range
     AddGenerator(genupsilonPP,"UpsilonPP", ratioupsilonPP); // Adding Generator
     fTotalRate+=ratioupsilonPP;
     genupsilonPP->Init(); // generation in selected kinematical range
     AddGenerator(genupsilonPP,"UpsilonPP", ratioupsilonPP); // Adding Generator
     fTotalRate+=ratioupsilonPP;

+
+//------------------------------------------------------------------
+// Generator of charm
+    
+    AliGenCorrHF *gencharm = new AliGenCorrHF(1, 4);  
+          gencharm->SetMomentumRange(0,9999);
+          gencharm->SetForceDecay(kAll);
+         ratioccbar = sigmaccbar/sigmaReaction;
+          gencharm->Init();
+          AddGenerator(gencharm,"CorrHFCharm",ratioccbar);
+         fTotalRate+=ratioccbar;
+

 //------------------------------------------------------------------
 //------------------------------------------------------------------

+// Generator of beauty
+
+         AliGenCorrHF *genbeauty = new AliGenCorrHF(1, 5);  
+          genbeauty->SetMomentumRange(0,9999);
+          genbeauty->SetForceDecay(kAll);
+         ratiobbbar = sigmabbbar/sigmaReaction;
+         genbeauty->Init();
+          AddGenerator(genbeauty,"CorrHFBeauty",ratiobbbar); 
+         fTotalRate+=ratiobbbar;
+
+//--------------------------t----------------------------------------

 // Pythia generator
     AliGenPythia *pythia = new AliGenPythia(1);
     pythia->SetProcess(kPyMbMSEL1);
 // Pythia generator
     AliGenPythia *pythia = new AliGenPythia(1);
     pythia->SetProcess(kPyMbMSEL1);


@@ -224,6 +252,7 @@ void AliGenMUONCocktailpp::Init()
     pythia->SetYRange(-8.,8.);
     pythia->SetPhiRange(0.,360.);
     pythia->SetPtHard(2.76,-1.0);
     pythia->SetYRange(-8.,8.);
     pythia->SetPhiRange(0.,360.);
     pythia->SetPtHard(2.76,-1.0);

+    pythia->SwitchHFOff();

     pythia->Init(); 
     AddGenerator(pythia,"Pythia",1);
     fTotalRate+=1.;
     pythia->Init(); 
     AddGenerator(pythia,"Pythia",1);
     fTotalRate+=1.;


@@ -290,9 +319,7 @@ void AliGenMUONCocktailpp::Generate()
 // in the muon spectrometer acceptance
        Int_t iPart;
        fNGenerated++;
 // in the muon spectrometer acceptance
        Int_t iPart;
        fNGenerated++;

-       Int_t numberOfMuons=0;  
-
-       Int_t maxPart = partArray->GetEntriesFast();
+       Int_t numberOfMuons=0;Int_t maxPart = partArray->GetEntriesFast();

        for(iPart=0; iPart<maxPart; iPart++){      
            
          TParticle *part = gAlice->GetMCApp()->Particle(iPart);
        for(iPart=0; iPart<maxPart; iPart++){      
            
          TParticle *part = gAlice->GetMCApp()->Particle(iPart);