[u/mrichter/AliRoot.git] / EVGEN / AliGenMUONCocktailpp.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

//
// Classe to create the coktail for physics with muons for pp at 14 TeV
// The followoing sources: 
// jpsi,psiP, upsilon, upsilonP, upsilonPP are added to Pythia
// 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 
// Aug. 07: added trigger cut on total momentum

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

#include "AliGenCocktailEventHeader.h"

#include "AliGenCocktailEntry.h"
#include "AliGenMUONCocktailpp.h"
#include "AliGenMUONlib.h"
#include "AliGenParam.h"
#include "AliMC.h"
#include "AliRun.h"
#include "AliStack.h"
#include "AliDecayer.h"
#include "AliLog.h"
#include "AliGenCorrHF.h"

ClassImp(AliGenMUONCocktailpp)  
  
//________________________________________________________________________
AliGenMUONCocktailpp::AliGenMUONCocktailpp()
    :AliGenCocktail(),
     fDecayer(0),
     fDecayModeResonance(kAll),
     fDecayModePythia(kAll),
     fMuonMultiplicity(0),
     fMuonPtCut(0.),
     fMuonPCut(0.),
     fMuonThetaMinCut(0.),
     fMuonThetaMaxCut(180.),
     fMuonOriginCut(-999.),
     fNSucceded(0),
     fNGenerated(0)
{
// Constructor
}

//_________________________________________________________________________
AliGenMUONCocktailpp::~AliGenMUONCocktailpp()
// Destructor
{
    
}

//_________________________________________________________________________
void AliGenMUONCocktailpp::CreateCocktail()
{
// MinBias NN cross section @ pp 14 TeV -PR  Vol II p:473
    Double_t sigmaReaction =   0.070; 
    
// These limits are only used for renormalization of quarkonia cross section
// Pythia events are generated in 4pi  
    Double_t ptMin  = fPtMin;
    Double_t ptMax  = fPtMax;
    Double_t yMin   = fYMin;;
    Double_t yMax   = fYMax;;
    Double_t phiMin = fPhiMin*180./TMath::Pi();
    Double_t phiMax = fPhiMax*180./TMath::Pi();
    AliInfo(Form("Ranges pT:%4.1f : %4.1f GeV/c, y:%4.2f : %4.2f, Phi:%5.1f : %5.1f degres",ptMin,ptMax,yMin,yMax,phiMin,phiMax));
    
// Ratios with respect to the reaction cross-section in the 
// kinematics limit of the MUONCocktail
    Double_t ratiojpsi;
    Double_t ratiopsiP;
    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 

    Double_t sigmajpsi = 49.44e-6;  
    Double_t sigmapsiP = 7.67e-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;

    Bool_t showSIGJPSI = kTRUE;
    if (gSystem->Getenv("COCKTAIL_SIGJPSI")) {
      sigmajpsi = atof(gSystem->Getenv("COCKTAIL_SIGJPSI"));
      AliInfo("Cross-section for JPsi set from external value");
      if (!gSystem->Getenv("COCKTAIL_SIGSHOW")) showSIGJPSI = kFALSE;
    }
    Bool_t showSIGPSIP = kTRUE;
    if (gSystem->Getenv("COCKTAIL_SIGPSIP")) {
      sigmapsiP = atof(gSystem->Getenv("COCKTAIL_SIGPSIP"));
      AliInfo("Cross-section for Psi-prime set from external value");
      if (!gSystem->Getenv("COCKTAIL_SIGSHOW")) showSIGPSIP = kFALSE;
    }
    Bool_t showSIGUPSI = kTRUE;
    if (gSystem->Getenv("COCKTAIL_SIGUPSI")) {
      sigmaupsilon = atof(gSystem->Getenv("COCKTAIL_SIGUPSI"));
      AliInfo("Cross-section for Upsilon 1S set from external value");
      if (!gSystem->Getenv("COCKTAIL_SIGSHOW")) showSIGUPSI = kFALSE;
    }
    Bool_t showSIGUPSIP = kTRUE;
    if (gSystem->Getenv("COCKTAIL_SIGUPSIP")) {
      sigmaupsilonP = atof(gSystem->Getenv("COCKTAIL_SIGUPSIP"));
      AliInfo("Cross-section for Upsilon 2S set from external value");
      if (!gSystem->Getenv("COCKTAIL_SIGSHOW")) showSIGUPSIP = kFALSE;
    }
    Bool_t showSIGUPSIPP = kTRUE;
    if (gSystem->Getenv("COCKTAIL_SIGUPSIPP")) {
      sigmaupsilonPP = atof(gSystem->Getenv("COCKTAIL_SIGUPSIPP"));
      AliInfo("Cross-section for Upsilon 3S set from external value");
      if (!gSystem->Getenv("COCKTAIL_SIGSHOW")) showSIGUPSIPP = kFALSE;
    }
    Bool_t showSIGCCBAR = kTRUE;
    if (gSystem->Getenv("COCKTAIL_SIGCCBAR")) {
      sigmaccbar = atof(gSystem->Getenv("COCKTAIL_SIGCCBAR"));
      AliInfo("Cross-section for c-cbar set from external value");
      if (!gSystem->Getenv("COCKTAIL_SIGSHOW")) showSIGCCBAR = kFALSE;
    }
    Bool_t showSIGBBBAR = kTRUE;
    if (gSystem->Getenv("COCKTAIL_SIGBBBAR")) {
      sigmabbbar = atof(gSystem->Getenv("COCKTAIL_SIGBBBAR"));
      AliInfo("Cross-section for b-bbar set from external value");
      if (!gSystem->Getenv("COCKTAIL_SIGSHOW")) showSIGBBBAR = kFALSE;
    }
   
    AliInfo(Form("the parametrised resonances uses the decay mode %d",fDecayModeResonance));

// Generation using CDF scaled distribution for pp@14TeV (from D. Stocco)
//----------------------------------------------------------------------
// Jpsi generator
    AliGenParam * genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp", "Jpsi");
// first step: generation in 4pi
    genjpsi->SetPtRange(0.,100.);
    genjpsi->SetYRange(-8.,8.);
    genjpsi->SetPhiRange(0.,360.);
    genjpsi->SetForceDecay(fDecayModeResonance);
    if (!gMC) genjpsi->SetDecayer(fDecayer);

    genjpsi->Init();  // generation in 4pi
    ratiojpsi = sigmajpsi / sigmaReaction * genjpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); // get weight
    if (showSIGJPSI) {
      AliInfo(Form("jpsi prod. cross-section in pp(14 TeV) %5.3g b",sigmajpsi));
      AliInfo(Form("jpsi prod. probability per collision in acceptance %5.3g",ratiojpsi));
    }
// second step: generation in selected kinematical range
    genjpsi->SetPtRange(ptMin, ptMax);  
    genjpsi->SetYRange(yMin, yMax);
    genjpsi->SetPhiRange(phiMin, phiMax);
    genjpsi->Init(); // generation in selected kinematic range
    AddGenerator(genjpsi, "Jpsi", ratiojpsi); // Adding Generator
//------------------------------------------------------------------
// Psi prime generator
    AliGenParam * genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp", "PsiP");
// first step: generation in 4pi
    genpsiP->SetPtRange(0.,100.);
    genpsiP->SetYRange(-8.,8.);
    genpsiP->SetPhiRange(0.,360.);
    genpsiP->SetForceDecay(fDecayModeResonance);
    if (!gMC) genpsiP->SetDecayer(fDecayer);

    genpsiP->Init();  // generation in 4pi
    ratiopsiP = sigmapsiP / sigmaReaction * genpsiP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
    if (showSIGPSIP) {
      AliInfo(Form("psiP prod. cross-section in pp(14 TeV) %5.3g b",sigmapsiP));
      AliInfo(Form("psiP prod. probability per collision in acceptance %5.3g",ratiopsiP));
    }
// second step: generation in selected kinematical range
    genpsiP->SetPtRange(ptMin, ptMax);  
    genpsiP->SetYRange(yMin, yMax);
    genpsiP->SetPhiRange(phiMin, phiMax);
    genpsiP->Init(); // generation in selected kinematic range
    AddGenerator(genpsiP, "PsiP", ratiopsiP); // Adding Generator
//------------------------------------------------------------------
// Upsilon 1S generator
    AliGenParam * genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp", "Upsilon");
// first step: generation in 4pi
    genupsilon->SetPtRange(0.,100.); 
    genupsilon->SetYRange(-8.,8.);
    genupsilon->SetPhiRange(0.,360.);
    genupsilon->SetForceDecay(fDecayModeResonance);
    if (!gMC) genupsilon->SetDecayer(fDecayer);
    genupsilon->Init();  // generation in 4pi
    ratioupsilon = sigmaupsilon / sigmaReaction * genupsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
    if (showSIGUPSI) {
      AliInfo(Form("Upsilon 1S prod. cross-section in pp(14 TeV) %5.3g b",sigmaupsilon));
      AliInfo(Form("Upsilon 1S prod. probability per collision in acceptance %5.3g",ratioupsilon));
    }
// second step: generation in selected kinematical range
    genupsilon->SetPtRange(ptMin, ptMax);  
    genupsilon->SetYRange(yMin, yMax);
    genupsilon->SetPhiRange(phiMin, phiMax);
    genupsilon->Init(); // generation in selected kinematical range
    AddGenerator(genupsilon,"Upsilon", ratioupsilon); // Adding Generator
//------------------------------------------------------------------
// Upsilon 2S generator
    AliGenParam * genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF pp", "UpsilonP");
// first step: generation in 4pi
    genupsilonP->SetPtRange(0.,100.);
    genupsilonP->SetYRange(-8.,8.);
    genupsilonP->SetPhiRange(0.,360.);
    genupsilonP->SetForceDecay(fDecayModeResonance);
    if (!gMC) genupsilonP->SetDecayer(fDecayer);  
    genupsilonP->Init();  // generation in 4pi
    ratioupsilonP = sigmaupsilonP / sigmaReaction * genupsilonP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
    if (showSIGUPSIP) {
      AliInfo(Form("Upsilon 2S prod. cross-section in pp(14 TeV) %5.3g b",sigmaupsilonP));
      AliInfo(Form("Upsilon 2S prod. probability per collision in acceptance %5.3g",ratioupsilonP));
    }
// second step: generation in the kinematical range
    genupsilonP->SetPtRange(ptMin, ptMax);  
    genupsilonP->SetYRange(yMin, yMax);
    genupsilonP->SetPhiRange(phiMin, phiMax);
    genupsilonP->Init(); // generation in selected kinematical range
    AddGenerator(genupsilonP,"UpsilonP", ratioupsilonP); // Adding Generator
    
//------------------------------------------------------------------
// Upsilon 3S generator
    AliGenParam * genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp", "UpsilonPP");
// first step: generation in 4pi
    genupsilonPP->SetPtRange(0.,100.); 
    genupsilonPP->SetYRange(-8.,8.);
    genupsilonPP->SetPhiRange(0.,360.);
    genupsilonPP->SetForceDecay(fDecayModeResonance);
    if (!gMC) genupsilonPP->SetDecayer(fDecayer);  
    genupsilonPP->Init();  // generation in 4pi
    ratioupsilonPP = sigmaupsilonPP / sigmaReaction * genupsilonPP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
    if (showSIGUPSIPP) {
      AliInfo(Form("Upsilon 3S prod. cross-section in pp(14 TeV) %5.3g b",sigmaupsilonPP));
      AliInfo(Form("Upsilon 3S prod. probability per collision in acceptance %5.3g",ratioupsilonPP));
    }
// second step: generation in selected kinematical range
    genupsilonPP->SetPtRange(ptMin, ptMax);  
    genupsilonPP->SetYRange(yMin, yMax);
    genupsilonPP->SetPhiRange(phiMin, phiMax);
    genupsilonPP->Init(); // generation in selected kinematical range
    AddGenerator(genupsilonPP,"UpsilonPP", ratioupsilonPP); // Adding Generator

//------------------------------------------------------------------
// Generator of charm
    
    AliGenCorrHF *gencharm = new AliGenCorrHF(1, 4);  
    gencharm->SetMomentumRange(0,9999);
    gencharm->SetForceDecay(kAll);
    ratioccbar = sigmaccbar/sigmaReaction;
    if (!gMC) gencharm->SetDecayer(fDecayer);  
    gencharm->Init();
    if (showSIGCCBAR) {
      AliInfo(Form("c-cbar prod. cross-section in pp(14 TeV) %5.3g b",sigmaccbar));
      AliInfo(Form("c-cbar prod. probability per collision in acceptance %5.3g",ratioccbar));
    }
    AddGenerator(gencharm,"CorrHFCharm",ratioccbar);

//------------------------------------------------------------------
// Generator of beauty

    AliGenCorrHF *genbeauty = new AliGenCorrHF(1, 5);  
    genbeauty->SetMomentumRange(0,9999);
    genbeauty->SetForceDecay(kAll);
    ratiobbbar = sigmabbbar/sigmaReaction;
    if (!gMC) genbeauty->SetDecayer(fDecayer);  
    genbeauty->Init();
    if (showSIGBBBAR) {
      AliInfo(Form("b-bbar prod. cross-section in pp(14 TeV) %5.3g b",sigmabbbar));
      AliInfo(Form("b-bbar prod. probability per collision in acceptance %5.3g",ratiobbbar));
    }
    AddGenerator(genbeauty,"CorrHFBeauty",ratiobbbar); 

//-------------------------------------------------------------------
// Pythia generator
//
// This has to go into the Config.C
//
//    AliGenPythia *pythia = new AliGenPythia(1);
//    pythia->SetProcess(kPyMbMSEL1);
//    pythia->SetStrucFunc(kCTEQ5L);
//    pythia->SetEnergyCMS(14000.);
//    AliInfo(Form("\n\npythia uses the decay mode %d", GetDecayModePythia()));
//    Decay_t dt = gener->GetDecayModePythia();
//    pythia->SetForceDecay(dt);
//    pythia->SetPtRange(0.,100.);
//    pythia->SetYRange(-8.,8.);
//    pythia->SetPhiRange(0.,360.);
//    pythia->SetPtHard(2.76,-1.0);
//    pythia->SwitchHFOff();
//    pythia->Init(); 
//    AddGenerator(pythia,"Pythia",1);
}

void AliGenMUONCocktailpp::Init()
{
    //
    // Initialisation
    TIter next(fEntries);
    AliGenCocktailEntry *entry;
    if (fStack) {
	while((entry = (AliGenCocktailEntry*)next())) {
	    entry->Generator()->SetStack(fStack);
	}
    }
}

//_________________________________________________________________________
void AliGenMUONCocktailpp::Generate()
{
// Generate event 
    TIter next(fEntries);
    AliGenCocktailEntry *entry = 0;
    AliGenCocktailEntry *preventry = 0;
    AliGenerator* gen = 0;

    const TObjArray *partArray = gAlice->GetMCApp()->Particles();
    
// Generate the vertex position used by all generators    
    if(fVertexSmear == kPerEvent) Vertex();

// Loop on primordialTrigger: 
// minimum muon multiplicity above a pt cut in a theta acceptance region

    Bool_t primordialTrigger = kFALSE;
    while(!primordialTrigger) {
	//Reseting stack
	AliRunLoader * runloader = AliRunLoader::Instance();
	if (runloader)
	    if (runloader->Stack())
		runloader->Stack()->Clean();
	// Loop over generators and generate events
	Int_t igen = 0;
	Int_t npart = 0;	
	const char* genName = 0;
	while((entry = (AliGenCocktailEntry*)next())) {
	    gen = entry->Generator();
	    genName = entry->GetName();
	    gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));

	    npart = (strcmp(genName,"Pythia") == 0) ? 1 :
		gRandom->Poisson(entry->Rate());

	    if (npart > 0) {
		igen++;	
		if (igen == 1) entry->SetFirst(0);		
		else  entry->SetFirst((partArray->GetEntriesFast())+1);
		
		gen->SetNumberParticles(npart);		
		gen->Generate();
		entry->SetLast(partArray->GetEntriesFast());
		preventry = entry;
	    }
	}  
	next.Reset();

// Testing primordial trigger: Single muons or dimuons with Pt above a Pt cut 
// in the muon spectrometer acceptance
	Int_t iPart;
	fNGenerated++;
	Int_t numberOfMuons=0;Int_t maxPart = partArray->GetEntriesFast();
	for(iPart=0; iPart<maxPart; iPart++){      
	    
	  TParticle *part = gAlice->GetMCApp()->Particle(iPart);
	  if ( TMath::Abs(part->GetPdgCode()) == 13 ){  
	    if((part->Vz() > fMuonOriginCut) && //take only the muons that decayed before the abs + 1 int. length in C abs
	       (part->Theta()*180./TMath::Pi()>fMuonThetaMinCut) &&
	       (part->Theta()*180./TMath::Pi()<fMuonThetaMaxCut) &&
	       (part->Pt()>fMuonPtCut) &&
	       (part->P()>fMuonPCut)) {
	      numberOfMuons++;
	    }
	  }
	}	
	if (numberOfMuons >= fMuonMultiplicity) primordialTrigger = kTRUE;
    }
    fNSucceded++;

//     AliInfo(Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
    AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
}
Commit	Line	Data
103ac317	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	//
	19	// Classe to create the coktail for physics with muons for pp at 14 TeV
	20	// The followoing sources:
	21	// jpsi,psiP, upsilon, upsilonP, upsilonPP are added to Pythia
	22	// The free parameeters are :
	23	// pp reaction cross-section
	24	// production cross-sections in pp collisions
b44c3901	25	// July 07:added heavy quark production from AliGenCorrHF and heavy quark
b44c3901	26	// production switched off in Pythia
f81a4aca	27	// Aug. 07: added trigger cut on total momentum
103ac317	28
	29	#include <TObjArray.h>
	30	#include <TParticle.h>
	31	#include <TF1.h>
aaa95f22	32	#include <TVirtualMC.h>
103ac317	33
	34	#include "AliGenCocktailEventHeader.h"
	35
	36	#include "AliGenCocktailEntry.h"
	37	#include "AliGenMUONCocktailpp.h"
	38	#include "AliGenMUONlib.h"
	39	#include "AliGenParam.h"
	40	#include "AliMC.h"
	41	#include "AliRun.h"
	42	#include "AliStack.h"
aaa95f22	43	#include "AliDecayer.h"
103ac317	44	#include "AliLog.h"
b44c3901	45	#include "AliGenCorrHF.h"
103ac317	46
	47	ClassImp(AliGenMUONCocktailpp)
	48
	49	//________________________________________________________________________
	50	AliGenMUONCocktailpp::AliGenMUONCocktailpp()
1c56e311	51	:AliGenCocktail(),
aaa95f22	52	fDecayer(0),
	53	fDecayModeResonance(kAll),
	54	fDecayModePythia(kAll),
1c56e311	55	fMuonMultiplicity(0),
1c56e311	56	fMuonPtCut(0.),
f81a4aca	57	fMuonPCut(0.),
18bc0c8e	58	fMuonThetaMinCut(0.),
f81a4aca	59	fMuonThetaMaxCut(180.),
aaa95f22	60	fMuonOriginCut(-999.),
18bc0c8e	61	fNSucceded(0),
1c56e311	62	fNGenerated(0)
103ac317	63	{
103ac317	64	// Constructor
103ac317	65	}
93a2041b	66
103ac317	67	//_________________________________________________________________________
	68	AliGenMUONCocktailpp::~AliGenMUONCocktailpp()
	69	// Destructor
	70	{
	71
	72	}
	73
	74	//_________________________________________________________________________
8058ead1	75	void AliGenMUONCocktailpp::CreateCocktail()
103ac317	76	{
	77	// MinBias NN cross section @ pp 14 TeV -PR Vol II p:473
	78	Double_t sigmaReaction = 0.070;
	79
	80	// These limits are only used for renormalization of quarkonia cross section
	81	// Pythia events are generated in 4pi
	82	Double_t ptMin = fPtMin;
	83	Double_t ptMax = fPtMax;
	84	Double_t yMin = fYMin;;
	85	Double_t yMax = fYMax;;
	86	Double_t phiMin = fPhiMin*180./TMath::Pi();
	87	Double_t phiMax = fPhiMax*180./TMath::Pi();
	88	AliInfo(Form("Ranges pT:%4.1f : %4.1f GeV/c, y:%4.2f : %4.2f, Phi:%5.1f : %5.1f degres",ptMin,ptMax,yMin,yMax,phiMin,phiMax));
	89
	90	// Ratios with respect to the reaction cross-section in the
	91	// kinematics limit of the MUONCocktail
	92	Double_t ratiojpsi;
	93	Double_t ratiopsiP;
	94	Double_t ratioupsilon;
	95	Double_t ratioupsilonP;
	96	Double_t ratioupsilonPP;
b44c3901	97	Double_t ratioccbar;
	98	Double_t ratiobbbar;
	99
18bc0c8e	100	// Cross sections in barns (from PPR Vol. II p: 552) pp - 14 TeV and
18bc0c8e	101	// corrected from feed down of higher resonances
103ac317	102
18bc0c8e	103	Double_t sigmajpsi = 49.44e-6;
103ac317	104	Double_t sigmapsiP = 7.67e-6;
18bc0c8e	105	Double_t sigmaupsilon = 0.989e-6;
18bc0c8e	106	Double_t sigmaupsilonP = 0.502e-6;
103ac317	107	Double_t sigmaupsilonPP = 0.228e-6;
b44c3901	108	Double_t sigmaccbar = 11.2e-3;
b44c3901	109	Double_t sigmabbbar = 0.51e-3;
3285b419	110
	111	Bool_t showSIGJPSI = kTRUE;
	112	if (gSystem->Getenv("COCKTAIL_SIGJPSI")) {
	113	sigmajpsi = atof(gSystem->Getenv("COCKTAIL_SIGJPSI"));
	114	AliInfo("Cross-section for JPsi set from external value");
	115	if (!gSystem->Getenv("COCKTAIL_SIGSHOW")) showSIGJPSI = kFALSE;
	116	}
	117	Bool_t showSIGPSIP = kTRUE;
	118	if (gSystem->Getenv("COCKTAIL_SIGPSIP")) {
	119	sigmapsiP = atof(gSystem->Getenv("COCKTAIL_SIGPSIP"));
	120	AliInfo("Cross-section for Psi-prime set from external value");
	121	if (!gSystem->Getenv("COCKTAIL_SIGSHOW")) showSIGPSIP = kFALSE;
	122	}
	123	Bool_t showSIGUPSI = kTRUE;
	124	if (gSystem->Getenv("COCKTAIL_SIGUPSI")) {
	125	sigmaupsilon = atof(gSystem->Getenv("COCKTAIL_SIGUPSI"));
	126	AliInfo("Cross-section for Upsilon 1S set from external value");
	127	if (!gSystem->Getenv("COCKTAIL_SIGSHOW")) showSIGUPSI = kFALSE;
	128	}
	129	Bool_t showSIGUPSIP = kTRUE;
	130	if (gSystem->Getenv("COCKTAIL_SIGUPSIP")) {
	131	sigmaupsilonP = atof(gSystem->Getenv("COCKTAIL_SIGUPSIP"));
	132	AliInfo("Cross-section for Upsilon 2S set from external value");
	133	if (!gSystem->Getenv("COCKTAIL_SIGSHOW")) showSIGUPSIP = kFALSE;
	134	}
	135	Bool_t showSIGUPSIPP = kTRUE;
	136	if (gSystem->Getenv("COCKTAIL_SIGUPSIPP")) {
	137	sigmaupsilonPP = atof(gSystem->Getenv("COCKTAIL_SIGUPSIPP"));
	138	AliInfo("Cross-section for Upsilon 3S set from external value");
	139	if (!gSystem->Getenv("COCKTAIL_SIGSHOW")) showSIGUPSIPP = kFALSE;
	140	}
	141	Bool_t showSIGCCBAR = kTRUE;
	142	if (gSystem->Getenv("COCKTAIL_SIGCCBAR")) {
	143	sigmaccbar = atof(gSystem->Getenv("COCKTAIL_SIGCCBAR"));
	144	AliInfo("Cross-section for c-cbar set from external value");
	145	if (!gSystem->Getenv("COCKTAIL_SIGSHOW")) showSIGCCBAR = kFALSE;
	146	}
	147	Bool_t showSIGBBBAR = kTRUE;
	148	if (gSystem->Getenv("COCKTAIL_SIGBBBAR")) {
	149	sigmabbbar = atof(gSystem->Getenv("COCKTAIL_SIGBBBAR"));
	150	AliInfo("Cross-section for b-bbar set from external value");
	151	if (!gSystem->Getenv("COCKTAIL_SIGSHOW")) showSIGBBBAR = kFALSE;
	152	}
	153
aaa95f22	154	AliInfo(Form("the parametrised resonances uses the decay mode %d",fDecayModeResonance));
aaa95f22	155
18bc0c8e	156	// Generation using CDF scaled distribution for pp@14TeV (from D. Stocco)
18bc0c8e	157	//----------------------------------------------------------------------
103ac317	158	// Jpsi generator
18bc0c8e	159	AliGenParam * genjpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "CDF pp", "Jpsi");
103ac317	160	// first step: generation in 4pi
	161	genjpsi->SetPtRange(0.,100.);
	162	genjpsi->SetYRange(-8.,8.);
	163	genjpsi->SetPhiRange(0.,360.);
aaa95f22	164	genjpsi->SetForceDecay(fDecayModeResonance);
	165	if (!gMC) genjpsi->SetDecayer(fDecayer);
	166
103ac317	167	genjpsi->Init(); // generation in 4pi
103ac317	168	ratiojpsi = sigmajpsi / sigmaReaction * genjpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax); // get weight
3285b419	169	if (showSIGJPSI) {
	170	AliInfo(Form("jpsi prod. cross-section in pp(14 TeV) %5.3g b",sigmajpsi));
	171	AliInfo(Form("jpsi prod. probability per collision in acceptance %5.3g",ratiojpsi));
	172	}
103ac317	173	// second step: generation in selected kinematical range
	174	genjpsi->SetPtRange(ptMin, ptMax);
	175	genjpsi->SetYRange(yMin, yMax);
	176	genjpsi->SetPhiRange(phiMin, phiMax);
	177	genjpsi->Init(); // generation in selected kinematic range
	178	AddGenerator(genjpsi, "Jpsi", ratiojpsi); // Adding Generator
103ac317	179	//------------------------------------------------------------------
103ac317	180	// Psi prime generator
18bc0c8e	181	AliGenParam * genpsiP = new AliGenParam(1, AliGenMUONlib::kPsiP, "CDF pp", "PsiP");
103ac317	182	// first step: generation in 4pi
	183	genpsiP->SetPtRange(0.,100.);
	184	genpsiP->SetYRange(-8.,8.);
	185	genpsiP->SetPhiRange(0.,360.);
aaa95f22	186	genpsiP->SetForceDecay(fDecayModeResonance);
	187	if (!gMC) genpsiP->SetDecayer(fDecayer);
	188
103ac317	189	genpsiP->Init(); // generation in 4pi
103ac317	190	ratiopsiP = sigmapsiP / sigmaReaction * genpsiP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
3285b419	191	if (showSIGPSIP) {
	192	AliInfo(Form("psiP prod. cross-section in pp(14 TeV) %5.3g b",sigmapsiP));
	193	AliInfo(Form("psiP prod. probability per collision in acceptance %5.3g",ratiopsiP));
	194	}
103ac317	195	// second step: generation in selected kinematical range
	196	genpsiP->SetPtRange(ptMin, ptMax);
	197	genpsiP->SetYRange(yMin, yMax);
	198	genpsiP->SetPhiRange(phiMin, phiMax);
	199	genpsiP->Init(); // generation in selected kinematic range
	200	AddGenerator(genpsiP, "PsiP", ratiopsiP); // Adding Generator
103ac317	201	//------------------------------------------------------------------
103ac317	202	// Upsilon 1S generator
18bc0c8e	203	AliGenParam * genupsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "CDF pp", "Upsilon");
103ac317	204	// first step: generation in 4pi
	205	genupsilon->SetPtRange(0.,100.);
	206	genupsilon->SetYRange(-8.,8.);
	207	genupsilon->SetPhiRange(0.,360.);
aaa95f22	208	genupsilon->SetForceDecay(fDecayModeResonance);
aaa95f22	209	if (!gMC) genupsilon->SetDecayer(fDecayer);
103ac317	210	genupsilon->Init(); // generation in 4pi
103ac317	211	ratioupsilon = sigmaupsilon / sigmaReaction * genupsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
3285b419	212	if (showSIGUPSI) {
	213	AliInfo(Form("Upsilon 1S prod. cross-section in pp(14 TeV) %5.3g b",sigmaupsilon));
	214	AliInfo(Form("Upsilon 1S prod. probability per collision in acceptance %5.3g",ratioupsilon));
	215	}
103ac317	216	// second step: generation in selected kinematical range
	217	genupsilon->SetPtRange(ptMin, ptMax);
	218	genupsilon->SetYRange(yMin, yMax);
	219	genupsilon->SetPhiRange(phiMin, phiMax);
	220	genupsilon->Init(); // generation in selected kinematical range
	221	AddGenerator(genupsilon,"Upsilon", ratioupsilon); // Adding Generator
103ac317	222	//------------------------------------------------------------------
103ac317	223	// Upsilon 2S generator
18bc0c8e	224	AliGenParam * genupsilonP = new AliGenParam(1, AliGenMUONlib::kUpsilonP, "CDF pp", "UpsilonP");
103ac317	225	// first step: generation in 4pi
	226	genupsilonP->SetPtRange(0.,100.);
	227	genupsilonP->SetYRange(-8.,8.);
	228	genupsilonP->SetPhiRange(0.,360.);
aaa95f22	229	genupsilonP->SetForceDecay(fDecayModeResonance);
aaa95f22	230	if (!gMC) genupsilonP->SetDecayer(fDecayer);
103ac317	231	genupsilonP->Init(); // generation in 4pi
103ac317	232	ratioupsilonP = sigmaupsilonP / sigmaReaction * genupsilonP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
3285b419	233	if (showSIGUPSIP) {
	234	AliInfo(Form("Upsilon 2S prod. cross-section in pp(14 TeV) %5.3g b",sigmaupsilonP));
	235	AliInfo(Form("Upsilon 2S prod. probability per collision in acceptance %5.3g",ratioupsilonP));
	236	}
103ac317	237	// second step: generation in the kinematical range
	238	genupsilonP->SetPtRange(ptMin, ptMax);
	239	genupsilonP->SetYRange(yMin, yMax);
	240	genupsilonP->SetPhiRange(phiMin, phiMax);
	241	genupsilonP->Init(); // generation in selected kinematical range
	242	AddGenerator(genupsilonP,"UpsilonP", ratioupsilonP); // Adding Generator
103ac317	243
	244	//------------------------------------------------------------------
	245	// Upsilon 3S generator
18bc0c8e	246	AliGenParam * genupsilonPP = new AliGenParam(1, AliGenMUONlib::kUpsilonPP, "CDF pp", "UpsilonPP");
103ac317	247	// first step: generation in 4pi
	248	genupsilonPP->SetPtRange(0.,100.);
	249	genupsilonPP->SetYRange(-8.,8.);
	250	genupsilonPP->SetPhiRange(0.,360.);
aaa95f22	251	genupsilonPP->SetForceDecay(fDecayModeResonance);
aaa95f22	252	if (!gMC) genupsilonPP->SetDecayer(fDecayer);
103ac317	253	genupsilonPP->Init(); // generation in 4pi
103ac317	254	ratioupsilonPP = sigmaupsilonPP / sigmaReaction * genupsilonPP->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
3285b419	255	if (showSIGUPSIPP) {
	256	AliInfo(Form("Upsilon 3S prod. cross-section in pp(14 TeV) %5.3g b",sigmaupsilonPP));
	257	AliInfo(Form("Upsilon 3S prod. probability per collision in acceptance %5.3g",ratioupsilonPP));
	258	}
103ac317	259	// second step: generation in selected kinematical range
	260	genupsilonPP->SetPtRange(ptMin, ptMax);
	261	genupsilonPP->SetYRange(yMin, yMax);
	262	genupsilonPP->SetPhiRange(phiMin, phiMax);
	263	genupsilonPP->Init(); // generation in selected kinematical range
	264	AddGenerator(genupsilonPP,"UpsilonPP", ratioupsilonPP); // Adding Generator
b44c3901	265
	266	//------------------------------------------------------------------
	267	// Generator of charm
	268
	269	AliGenCorrHF *gencharm = new AliGenCorrHF(1, 4);
8058ead1	270	gencharm->SetMomentumRange(0,9999);
	271	gencharm->SetForceDecay(kAll);
	272	ratioccbar = sigmaccbar/sigmaReaction;
	273	if (!gMC) gencharm->SetDecayer(fDecayer);
	274	gencharm->Init();
3285b419	275	if (showSIGCCBAR) {
	276	AliInfo(Form("c-cbar prod. cross-section in pp(14 TeV) %5.3g b",sigmaccbar));
	277	AliInfo(Form("c-cbar prod. probability per collision in acceptance %5.3g",ratioccbar));
	278	}
8058ead1	279	AddGenerator(gencharm,"CorrHFCharm",ratioccbar);
b44c3901	280
103ac317	281	//------------------------------------------------------------------
b44c3901	282	// Generator of beauty
b44c3901	283
8058ead1	284	AliGenCorrHF *genbeauty = new AliGenCorrHF(1, 5);
	285	genbeauty->SetMomentumRange(0,9999);
	286	genbeauty->SetForceDecay(kAll);
	287	ratiobbbar = sigmabbbar/sigmaReaction;
	288	if (!gMC) genbeauty->SetDecayer(fDecayer);
	289	genbeauty->Init();
3285b419	290	if (showSIGBBBAR) {
	291	AliInfo(Form("b-bbar prod. cross-section in pp(14 TeV) %5.3g b",sigmabbbar));
	292	AliInfo(Form("b-bbar prod. probability per collision in acceptance %5.3g",ratiobbbar));
	293	}
8058ead1	294	AddGenerator(genbeauty,"CorrHFBeauty",ratiobbbar);
b44c3901	295
8058ead1	296	//-------------------------------------------------------------------
103ac317	297	// Pythia generator
8058ead1	298	//
	299	// This has to go into the Config.C
	300	//
	301	// AliGenPythia *pythia = new AliGenPythia(1);
	302	// pythia->SetProcess(kPyMbMSEL1);
	303	// pythia->SetStrucFunc(kCTEQ5L);
	304	// pythia->SetEnergyCMS(14000.);
	305	// AliInfo(Form("\n\npythia uses the decay mode %d", GetDecayModePythia()));
	306	// Decay_t dt = gener->GetDecayModePythia();
	307	// pythia->SetForceDecay(dt);
	308	// pythia->SetPtRange(0.,100.);
	309	// pythia->SetYRange(-8.,8.);
	310	// pythia->SetPhiRange(0.,360.);
	311	// pythia->SetPtHard(2.76,-1.0);
	312	// pythia->SwitchHFOff();
	313	// pythia->Init();
	314	// AddGenerator(pythia,"Pythia",1);
	315	}
103ac317	316
8058ead1	317	void AliGenMUONCocktailpp::Init()
	318	{
	319	//
	320	// Initialisation
aaa95f22	321	TIter next(fEntries);
	322	AliGenCocktailEntry *entry;
	323	if (fStack) {
	324	while((entry = (AliGenCocktailEntry*)next())) {
	325	entry->Generator()->SetStack(fStack);
	326	}
	327	}
103ac317	328	}
	329
	330	//_________________________________________________________________________
	331	void AliGenMUONCocktailpp::Generate()
	332	{
	333	// Generate event
	334	TIter next(fEntries);
	335	AliGenCocktailEntry *entry = 0;
	336	AliGenCocktailEntry *preventry = 0;
	337	AliGenerator* gen = 0;
	338
d94af0c1	339	const TObjArray *partArray = gAlice->GetMCApp()->Particles();
103ac317	340
	341	// Generate the vertex position used by all generators
	342	if(fVertexSmear == kPerEvent) Vertex();
	343
	344	// Loop on primordialTrigger:
	345	// minimum muon multiplicity above a pt cut in a theta acceptance region
	346
	347	Bool_t primordialTrigger = kFALSE;
	348	while(!primordialTrigger) {
	349	//Reseting stack
33c3c91a	350	AliRunLoader * runloader = AliRunLoader::Instance();
103ac317	351	if (runloader)
103ac317	352	if (runloader->Stack())
34da8494	353	runloader->Stack()->Clean();
103ac317	354	// Loop over generators and generate events
	355	Int_t igen = 0;
	356	Int_t npart = 0;
	357	const char* genName = 0;
	358	while((entry = (AliGenCocktailEntry*)next())) {
	359	gen = entry->Generator();
	360	genName = entry->GetName();
	361	gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
	362
	363	npart = (strcmp(genName,"Pythia") == 0) ? 1 :
	364	gRandom->Poisson(entry->Rate());
	365
	366	if (npart > 0) {
	367	igen++;
	368	if (igen == 1) entry->SetFirst(0);
	369	else entry->SetFirst((partArray->GetEntriesFast())+1);
	370
	371	gen->SetNumberParticles(npart);
	372	gen->Generate();
	373	entry->SetLast(partArray->GetEntriesFast());
	374	preventry = entry;
	375	}
	376	}
	377	next.Reset();
	378
	379	// Testing primordial trigger: Single muons or dimuons with Pt above a Pt cut
	380	// in the muon spectrometer acceptance
	381	Int_t iPart;
	382	fNGenerated++;
b44c3901	383	Int_t numberOfMuons=0;Int_t maxPart = partArray->GetEntriesFast();
aaa95f22	384	for(iPart=0; iPart<maxPart; iPart++){
103ac317	385
aaa95f22	386	TParticle *part = gAlice->GetMCApp()->Particle(iPart);
	387	if ( TMath::Abs(part->GetPdgCode()) == 13 ){
	388	if((part->Vz() > fMuonOriginCut) && //take only the muons that decayed before the abs + 1 int. length in C abs
	389	(part->Theta()*180./TMath::Pi()>fMuonThetaMinCut) &&
	390	(part->Theta()*180./TMath::Pi()<fMuonThetaMaxCut) &&
f81a4aca	391	(part->Pt()>fMuonPtCut) &&
f81a4aca	392	(part->P()>fMuonPCut)) {
aaa95f22	393	numberOfMuons++;
103ac317	394	}
aaa95f22	395	}
aaa95f22	396	}
103ac317	397	if (numberOfMuons >= fMuonMultiplicity) primordialTrigger = kTRUE;
	398	}
	399	fNSucceded++;
	400
aaa95f22	401	// AliInfo(Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
103ac317	402	AliDebug(5,Form("Generated Events are %d and Succeeded Events are %d",fNGenerated,fNSucceded));
	403	}
	404
	405