[u/mrichter/AliRoot.git] / EVGEN / AliGenMUONCocktail.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 MUON coktail for physics in the Alice muon spectrometer
// Gines Martinez, jan 2004, Nantes martinez@in2p3.fr


//

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

#include "AliGenParam.h"
#include "AliGenMUONlib.h"
#include "AliGenMUONCocktail.h"
#include "AliGenCocktailEntry.h"
#include "AliCollisionGeometry.h"
#include "AliRun.h"
#include "AliMC.h"
#include "AliStack.h"

ClassImp(AliGenMUONCocktail)
  
  
//________________________________________________________________________
AliGenMUONCocktail::AliGenMUONCocktail()
  :AliGenCocktail()
{
// Constructor
  fTotalRate =0;  
  fNSucceded=0; 
  fNGenerated=0; 
  fMuonMultiplicity=1;
  fMuonPtCut= 1.;
  fMuonThetaMinCut=171.; 
  fMuonThetaMaxCut=178.;
  fNumberOfCollisions = 400; // Minimum bias Pb+Pb collisions
    //
}
//_________________________________________________________________________
AliGenMUONCocktail::AliGenMUONCocktail(const AliGenMUONCocktail & cocktail):
    AliGenCocktail(cocktail)
{
// Copy constructor
  fTotalRate =0; 
  fNSucceded=0; 
  fNGenerated=0; 
  fMuonMultiplicity=1;
  fMuonPtCut= 1.;
  fMuonThetaMinCut=171.; 
  fMuonThetaMaxCut=178.;
  fNumberOfCollisions = 400; // Minimum bias Pb+Pb collisions
    //
}
//_________________________________________________________________________
AliGenMUONCocktail::~AliGenMUONCocktail()
{
// Destructor
    delete fEntries;
}

//_________________________________________________________________________
void AliGenMUONCocktail::Init()
{
  // Defining MUON physics cocktail

  // Kinematical limits for particle generation
  Float_t ptMin  = fPtMin;
  Float_t ptMax  = fPtMax;
  Float_t yMin   = fYMin;;
  Float_t yMax   = fYMax;;
  Float_t phiMin = fPhiMin*180./TMath::Pi();
  Float_t phiMax = fPhiMax*180./TMath::Pi();
  printf(">>> Kinematical range pT:%f:%f  y:%f:%f Phi:%f:%f\n",ptMin,ptMax,yMin,yMax,phiMin,phiMax);

  Float_t sigma_reaction =   0.072;   // MINB pp at LHC energies 72 mb
  
  // Generating J/Psi Physics
  AliGenParam * gen_jpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "Vogt", "Jpsi");
  // 4pi Generation 
  gen_jpsi->SetPtRange(0,100.);
  gen_jpsi->SetYRange(-8.,8);
  gen_jpsi->SetPhiRange(0.,360.);
  gen_jpsi->SetForceDecay(kDiMuon);
  gen_jpsi->SetTrackingFlag(1);
  // Calculation of the paritcle multiplicity per event in the muonic channel
  Float_t ratio_jpsi; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
  Float_t sigma_jpsi     = 31.0e-6 * 0.437;   //   section "6.7 Quarkonia Production" table 6.5 for pp  times shadowing
  Float_t br_jpsi        = 0.0588;           // Branching Ratio for JPsi
  gen_jpsi->Init();  // Generating pT and Y parametrsation for the 4pi
  ratio_jpsi = sigma_jpsi * br_jpsi * fNumberOfCollisions / sigma_reaction * gen_jpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
  printf(">>> ratio jpsi %g et %g Ncol %g sigma %g\n",ratio_jpsi,gen_jpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax),fNumberOfCollisions, sigma_jpsi );
  // Generation in the kinematical limits of AliGenMUONCocktail
  gen_jpsi->SetPtRange(ptMin, ptMax);  
  gen_jpsi->SetYRange(yMin, yMax);
  gen_jpsi->SetPhiRange(phiMin, phiMax);
  gen_jpsi->Init(); // Generating pT and Y parametrsation in the desired kinematic range
  // Adding Generator 
  AddGenerator(gen_jpsi, "Jpsi", ratio_jpsi);
  fTotalRate+=ratio_jpsi;

// Generating Upsilon Physics
  AliGenParam * gen_upsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "Vogt", "Upsilon");  
  gen_upsilon->SetPtRange(0,100.);  
  gen_upsilon->SetYRange(-8.,8);
  gen_upsilon->SetPhiRange(0.,360.);
  gen_upsilon->SetForceDecay(kDiMuon);
  gen_upsilon->SetTrackingFlag(1);
  Float_t ratio_upsilon; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
  Float_t sigma_upsilon     = 0.501e-6 * 0.674;   //  section "6.7 Quarkonia Production" table 6.5 for pp  times shadowing
  Float_t br_upsilon        = 0.0248;  // Branching Ratio for Upsilon
  gen_upsilon->Init();  // Generating pT and Y parametrsation for the 4pi
  ratio_upsilon = sigma_upsilon * br_upsilon * fNumberOfCollisions / sigma_reaction * gen_upsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
  printf(">>> ratio upsilon %g et %g\n",ratio_upsilon, gen_upsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax));
  gen_upsilon->SetPtRange(ptMin, ptMax);  
  gen_upsilon->SetYRange(yMin, yMax);
  gen_upsilon->SetPhiRange(phiMin, phiMax);
  gen_upsilon->Init(); // Generating pT and Y parametrsation in the desired kinematic range
  AddGenerator(gen_upsilon,"Upsilon", ratio_upsilon);
  fTotalRate+=ratio_upsilon;

// Generating Charm Physics 
  AliGenParam * gen_charm = new AliGenParam(1, AliGenMUONlib::kCharm, "Vogt", "Charm");  
  gen_charm->SetPtRange(0,100.);  
  gen_charm->SetYRange(-8.,8);
  gen_charm->SetPhiRange(0.,360.);
  gen_charm->SetForceDecay(kSemiMuonic);
  gen_charm->SetTrackingFlag(1);
  Float_t ratio_charm; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
  Float_t sigma_charm     = 2. * 6.64e-3 * 0.65 ;   // 
  Float_t br_charm        = 0.12;  // Branching Ratio for Charm
  gen_charm->Init();  // Generating pT and Y parametrsation for the 4pi
  ratio_charm = sigma_charm * br_charm * fNumberOfCollisions / sigma_reaction * 
    gen_charm->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
  gen_charm->SetPtRange(ptMin, ptMax);  
  gen_charm->SetYRange(yMin, yMax);
  gen_charm->SetPhiRange(phiMin, phiMax);
  gen_charm->Init(); // Generating pT and Y parametrsation in the desired kinematic range

  printf(">>> ratio charm %f\n",ratio_charm);
  AddGenerator(gen_charm,"Charm", ratio_charm);
  fTotalRate+=ratio_charm;

// Generating Beauty Physics "Correlated Pairs"
  AliGenParam * gen_beauty = new AliGenParam(2, AliGenMUONlib::kBeauty, "Vogt", "Beauty");  
  gen_beauty->SetPtRange(0,100.);  
  gen_beauty->SetYRange(-8.,8);
  gen_beauty->SetPhiRange(0.,360.);
  gen_beauty->SetForceDecay(kSemiMuonic);
  gen_beauty->SetTrackingFlag(1);
  Float_t ratio_beauty; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
  Float_t sigma_beauty     = 2. * 0.210e-3 * 0.84;   // 
  Float_t br_beauty        = 0.15;  // Branching Ratio for Beauty
  gen_beauty->Init();  // Generating pT and Y parametrsation for the 4pi
  ratio_beauty = sigma_beauty * br_beauty * fNumberOfCollisions / sigma_reaction * 
    gen_beauty->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
  gen_beauty->SetPtRange(ptMin, ptMax);  
  gen_beauty->SetYRange(yMin, yMax);
  gen_beauty->SetPhiRange(phiMin, phiMax);
  gen_beauty->Init(); // Generating pT and Y parametrisation in the desired kinematic range

  printf(">>> ratio beauty %f\n",ratio_beauty);
  AddGenerator(gen_beauty,"Beauty", ratio_beauty);
  fTotalRate+=ratio_beauty;

// Generating Pion Physics
  AliGenParam * gen_pion = new AliGenParam(1, AliGenMUONlib::kPion, "Vogt", "Pion");  
  gen_pion->SetPtRange(0,100.);  
  gen_pion->SetYRange(-8.,8);
  gen_pion->SetPhiRange(0.,360.);
  gen_pion->SetForceDecay(kPiToMu);
  gen_pion->SetTrackingFlag(1);
  Float_t ratio_pion; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
  Float_t sigma_pion     = 1.6e-2;   // A ojo TO be studied in detail.  
  Float_t br_pion        = 0.9999;  // Branching Ratio for Pion 
  gen_pion->Init();  // Generating pT and Y parametrsation for the 4pi
  ratio_pion = sigma_pion * br_pion * fNumberOfParticipants / sigma_reaction * gen_pion->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
  gen_pion->SetPtRange(ptMin, ptMax);  
  gen_pion->SetYRange(yMin, yMax);
  gen_pion->SetPhiRange(phiMin, phiMax);
  gen_pion->Init(); // Generating pT and Y parametrsation in the desired kinematic range
  printf(">>> ratio pion %f\n",ratio_pion);
  AddGenerator(gen_pion,"Pion", ratio_pion);
  fTotalRate+=ratio_pion;

// Generating Kaon Physics
  AliGenParam * gen_kaon = new AliGenParam(1, AliGenMUONlib::kKaon, "Vogt", "Kaon");  
  gen_kaon->SetPtRange(0,100.);  
  gen_kaon->SetYRange(-8.,8);
  gen_kaon->SetPhiRange(0.,360.);
  gen_kaon->SetForceDecay(kKaToMu);
  gen_kaon->SetTrackingFlag(1);
  Float_t ratio_kaon; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
  Float_t sigma_kaon     = 1.8e-4;   // OJO 
  Float_t br_kaon        = 0.6351 ;  // Branching Ratio for Kaon 
  gen_kaon->Init();  // Generating pT and Y parametrsation for the 4pi
  ratio_kaon = sigma_kaon * br_kaon * fNumberOfParticipants/ sigma_reaction * gen_kaon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
  gen_kaon->SetPtRange(ptMin, ptMax);  
  gen_kaon->SetYRange(yMin, yMax);
  gen_kaon->SetPhiRange(phiMin, phiMax);
  gen_kaon->Init(); // Generating pT and Y parametrsation in the desired kinematic range
  printf(">>> ratio kaon %f\n",ratio_kaon);
  AddGenerator(gen_kaon,"Kaon", ratio_kaon);
  fTotalRate+=ratio_kaon;
}

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

    TObjArray *partArray = gAlice->GetMCApp()->Particles();

//
//  Generate the vertex position used by all generators
//    
    if(fVertexSmear == kPerEvent) Vertex();
    Bool_t PrimordialTrigger = kFALSE;

    while(!PrimordialTrigger) {

      //Reseting stack
      AliRunLoader * runloader = gAlice->GetRunLoader();
      if (runloader)
	if (runloader->Stack())
	  runloader->Stack()->Reset();
      //
      // Loop over generators and generate events
      Int_t igen=0;
      Int_t npart =0;
      
      while((entry = (AliGenCocktailEntry*)next())) {
	gen = entry->Generator();
	gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
	if ( (npart = gRandom->Poisson(entry->Rate())) >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();

      // Tesitng primordial trigger : Muon  pair in the MUON spectrometer acceptance 171,178 and pTCut
      Int_t iPart;
      fNGenerated++;
      Int_t numberOfMuons=0;
      //      printf(">>>fNGenerated is %d\n",fNGenerated);
      for(iPart=0; iPart<partArray->GetEntriesFast(); iPart++){      
	//	gAlice->GetMCApp()->Particle(iPart)->Print();
	if ( (TMath::Abs(gAlice->GetMCApp()->Particle(iPart)->GetPdgCode())==13)  &&
	     (gAlice->GetMCApp()->Particle(iPart)->Theta()*180./TMath::Pi()>fMuonThetaMinCut) &&
	     (gAlice->GetMCApp()->Particle(iPart)->Theta()*180./TMath::Pi()<fMuonThetaMaxCut) &&
	     (gAlice->GetMCApp()->Particle(iPart)->Pt()>fMuonPtCut)                             ) { 
	  gAlice->GetMCApp()->Particle(iPart)->SetProductionVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2), 0.);   
	  numberOfMuons++;
	}
      }
      //  printf(">>> Number of Muons is %d \n", numberOfMuons);
      if (numberOfMuons >= fMuonMultiplicity ) PrimordialTrigger = kTRUE;
    }
    //printf(">>> Trigger Accepted!!! \n");
    fNSucceded++;
    //   Float_t Ratio = (Float_t) fNSucceded/fNGenerated;
    //    printf("Generated %d, Succeded %d and Ratio %f\n",fNGenerated, fNSucceded,Ratio);
}
Commit	Line	Data
84954c47	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	// Classe to create the MUON coktail for physics in the Alice muon spectrometer
	19	// Gines Martinez, jan 2004, Nantes martinez@in2p3.fr
	20
	21
	22	//
	23
	24	#include <TList.h>
	25	#include <TObjArray.h>
	26	#include <TF1.h>
	27	#include <TParticle.h>
	28
	29	#include "AliGenParam.h"
	30	#include "AliGenMUONlib.h"
	31	#include "AliGenMUONCocktail.h"
	32	#include "AliGenCocktailEntry.h"
	33	#include "AliCollisionGeometry.h"
	34	#include "AliRun.h"
	35	#include "AliMC.h"
	36	#include "AliStack.h"
	37
	38	ClassImp(AliGenMUONCocktail)
	39
	40
	41	//________________________________________________________________________
	42	AliGenMUONCocktail::AliGenMUONCocktail()
	43	:AliGenCocktail()
	44	{
	45	// Constructor
	46	fTotalRate =0;
	47	fNSucceded=0;
	48	fNGenerated=0;
	49	fMuonMultiplicity=1;
	50	fMuonPtCut= 1.;
	51	fMuonThetaMinCut=171.;
	52	fMuonThetaMaxCut=178.;
	53	fNumberOfCollisions = 400; // Minimum bias Pb+Pb collisions
	54	//
	55	}
	56	//_________________________________________________________________________
	57	AliGenMUONCocktail::AliGenMUONCocktail(const AliGenMUONCocktail & cocktail):
	58	AliGenCocktail(cocktail)
	59	{
	60	// Copy constructor
	61	fTotalRate =0;
	62	fNSucceded=0;
	63	fNGenerated=0;
	64	fMuonMultiplicity=1;
65	fMuonPtCut= 1.;
66	fMuonThetaMinCut=171.;
67	fMuonThetaMaxCut=178.;
68	fNumberOfCollisions = 400; // Minimum bias Pb+Pb collisions
69	//
70	}
71	//_________________________________________________________________________
72	AliGenMUONCocktail::~AliGenMUONCocktail()
73	{
74	// Destructor
75	delete fEntries;
76	}
77
78	//_________________________________________________________________________
79	void AliGenMUONCocktail::Init()
80	{
81	// Defining MUON physics cocktail
82
83	// Kinematical limits for particle generation
84	Float_t ptMin = fPtMin;
85	Float_t ptMax = fPtMax;
86	Float_t yMin = fYMin;;
87	Float_t yMax = fYMax;;
88	Float_t phiMin = fPhiMin*180./TMath::Pi();
89	Float_t phiMax = fPhiMax*180./TMath::Pi();
90	printf(">>> Kinematical range pT:%f:%f y:%f:%f Phi:%f:%f\n",ptMin,ptMax,yMin,yMax,phiMin,phiMax);
91
92	Float_t sigma_reaction = 0.072; // MINB pp at LHC energies 72 mb
93
94	// Generating J/Psi Physics
95	AliGenParam * gen_jpsi = new AliGenParam(1, AliGenMUONlib::kJpsi, "Vogt", "Jpsi");
96	// 4pi Generation
97	gen_jpsi->SetPtRange(0,100.);
98	gen_jpsi->SetYRange(-8.,8);
99	gen_jpsi->SetPhiRange(0.,360.);
100	gen_jpsi->SetForceDecay(kDiMuon);
101	gen_jpsi->SetTrackingFlag(1);
102	// Calculation of the paritcle multiplicity per event in the muonic channel
103	Float_t ratio_jpsi; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
104	Float_t sigma_jpsi = 31.0e-6 * 0.437; // section "6.7 Quarkonia Production" table 6.5 for pp times shadowing
105	Float_t br_jpsi = 0.0588; // Branching Ratio for JPsi
106	gen_jpsi->Init(); // Generating pT and Y parametrsation for the 4pi
107	ratio_jpsi = sigma_jpsi * br_jpsi * fNumberOfCollisions / sigma_reaction * gen_jpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
108	printf(">>> ratio jpsi %g et %g Ncol %g sigma %g\n",ratio_jpsi,gen_jpsi->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax),fNumberOfCollisions, sigma_jpsi );
109	// Generation in the kinematical limits of AliGenMUONCocktail
110	gen_jpsi->SetPtRange(ptMin, ptMax);
111	gen_jpsi->SetYRange(yMin, yMax);
112	gen_jpsi->SetPhiRange(phiMin, phiMax);
113	gen_jpsi->Init(); // Generating pT and Y parametrsation in the desired kinematic range
114	// Adding Generator
115	AddGenerator(gen_jpsi, "Jpsi", ratio_jpsi);
116	fTotalRate+=ratio_jpsi;
117
118	// Generating Upsilon Physics
119	AliGenParam * gen_upsilon = new AliGenParam(1, AliGenMUONlib::kUpsilon, "Vogt", "Upsilon");
120	gen_upsilon->SetPtRange(0,100.);
121	gen_upsilon->SetYRange(-8.,8);
122	gen_upsilon->SetPhiRange(0.,360.);
123	gen_upsilon->SetForceDecay(kDiMuon);
124	gen_upsilon->SetTrackingFlag(1);
125	Float_t ratio_upsilon; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
126	Float_t sigma_upsilon = 0.501e-6 * 0.674; // section "6.7 Quarkonia Production" table 6.5 for pp times shadowing
127	Float_t br_upsilon = 0.0248; // Branching Ratio for Upsilon
128	gen_upsilon->Init(); // Generating pT and Y parametrsation for the 4pi
129	ratio_upsilon = sigma_upsilon * br_upsilon * fNumberOfCollisions / sigma_reaction * gen_upsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
130	printf(">>> ratio upsilon %g et %g\n",ratio_upsilon, gen_upsilon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax));
131	gen_upsilon->SetPtRange(ptMin, ptMax);
132	gen_upsilon->SetYRange(yMin, yMax);
133	gen_upsilon->SetPhiRange(phiMin, phiMax);
134	gen_upsilon->Init(); // Generating pT and Y parametrsation in the desired kinematic range
135	AddGenerator(gen_upsilon,"Upsilon", ratio_upsilon);
136	fTotalRate+=ratio_upsilon;
137
138	// Generating Charm Physics
139	AliGenParam * gen_charm = new AliGenParam(1, AliGenMUONlib::kCharm, "Vogt", "Charm");
140	gen_charm->SetPtRange(0,100.);
141	gen_charm->SetYRange(-8.,8);
142	gen_charm->SetPhiRange(0.,360.);
143	gen_charm->SetForceDecay(kSemiMuonic);
144	gen_charm->SetTrackingFlag(1);
145	Float_t ratio_charm; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
146	Float_t sigma_charm = 2. * 6.64e-3 * 0.65 ; //
147	Float_t br_charm = 0.12; // Branching Ratio for Charm
148	gen_charm->Init(); // Generating pT and Y parametrsation for the 4pi
149	ratio_charm = sigma_charm * br_charm * fNumberOfCollisions / sigma_reaction *
150	gen_charm->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
151	gen_charm->SetPtRange(ptMin, ptMax);
152	gen_charm->SetYRange(yMin, yMax);
153	gen_charm->SetPhiRange(phiMin, phiMax);
154	gen_charm->Init(); // Generating pT and Y parametrsation in the desired kinematic range
155
156	printf(">>> ratio charm %f\n",ratio_charm);
157	AddGenerator(gen_charm,"Charm", ratio_charm);
158	fTotalRate+=ratio_charm;
159
160	// Generating Beauty Physics "Correlated Pairs"
161	AliGenParam * gen_beauty = new AliGenParam(2, AliGenMUONlib::kBeauty, "Vogt", "Beauty");
162	gen_beauty->SetPtRange(0,100.);
163	gen_beauty->SetYRange(-8.,8);
164	gen_beauty->SetPhiRange(0.,360.);
165	gen_beauty->SetForceDecay(kSemiMuonic);
166	gen_beauty->SetTrackingFlag(1);
167	Float_t ratio_beauty; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
168	Float_t sigma_beauty = 2. * 0.210e-3 * 0.84; //
169	Float_t br_beauty = 0.15; // Branching Ratio for Beauty
170	gen_beauty->Init(); // Generating pT and Y parametrsation for the 4pi
171	ratio_beauty = sigma_beauty * br_beauty * fNumberOfCollisions / sigma_reaction *
172	gen_beauty->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
173	gen_beauty->SetPtRange(ptMin, ptMax);
174	gen_beauty->SetYRange(yMin, yMax);
175	gen_beauty->SetPhiRange(phiMin, phiMax);
176	gen_beauty->Init(); // Generating pT and Y parametrisation in the desired kinematic range
177
178	printf(">>> ratio beauty %f\n",ratio_beauty);
179	AddGenerator(gen_beauty,"Beauty", ratio_beauty);
180	fTotalRate+=ratio_beauty;
181
182	// Generating Pion Physics
183	AliGenParam * gen_pion = new AliGenParam(1, AliGenMUONlib::kPion, "Vogt", "Pion");
184	gen_pion->SetPtRange(0,100.);
185	gen_pion->SetYRange(-8.,8);
186	gen_pion->SetPhiRange(0.,360.);
187	gen_pion->SetForceDecay(kPiToMu);
188	gen_pion->SetTrackingFlag(1);
189	Float_t ratio_pion; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
190	Float_t sigma_pion = 1.6e-2; // A ojo TO be studied in detail.
191	Float_t br_pion = 0.9999; // Branching Ratio for Pion
192	gen_pion->Init(); // Generating pT and Y parametrsation for the 4pi
193	ratio_pion = sigma_pion * br_pion * fNumberOfParticipants / sigma_reaction * gen_pion->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
194	gen_pion->SetPtRange(ptMin, ptMax);
195	gen_pion->SetYRange(yMin, yMax);
196	gen_pion->SetPhiRange(phiMin, phiMax);
197	gen_pion->Init(); // Generating pT and Y parametrsation in the desired kinematic range
198	printf(">>> ratio pion %f\n",ratio_pion);
199	AddGenerator(gen_pion,"Pion", ratio_pion);
200	fTotalRate+=ratio_pion;
201
202	// Generating Kaon Physics
203	AliGenParam * gen_kaon = new AliGenParam(1, AliGenMUONlib::kKaon, "Vogt", "Kaon");
204	gen_kaon->SetPtRange(0,100.);
205	gen_kaon->SetYRange(-8.,8);
206	gen_kaon->SetPhiRange(0.,360.);
207	gen_kaon->SetForceDecay(kKaToMu);
208	gen_kaon->SetTrackingFlag(1);
209	Float_t ratio_kaon; // Ratio with respect to the reaction cross-section for the muonic channel in the kinematics limit of the MUONCocktail
210	Float_t sigma_kaon = 1.8e-4; // OJO
211	Float_t br_kaon = 0.6351 ; // Branching Ratio for Kaon
212	gen_kaon->Init(); // Generating pT and Y parametrsation for the 4pi
213	ratio_kaon = sigma_kaon * br_kaon * fNumberOfParticipants/ sigma_reaction * gen_kaon->GetRelativeArea(ptMin,ptMax,yMin,yMax,phiMin,phiMax);
214	gen_kaon->SetPtRange(ptMin, ptMax);
215	gen_kaon->SetYRange(yMin, yMax);
216	gen_kaon->SetPhiRange(phiMin, phiMax);
217	gen_kaon->Init(); // Generating pT and Y parametrsation in the desired kinematic range
218	printf(">>> ratio kaon %f\n",ratio_kaon);
219	AddGenerator(gen_kaon,"Kaon", ratio_kaon);
220	fTotalRate+=ratio_kaon;
221	}
222
223	//_________________________________________________________________________
224	void AliGenMUONCocktail::Generate()
225	{
226	//
227	// Generate event
228	TIter next(fEntries);
229	AliGenCocktailEntry *entry = 0;
230	AliGenCocktailEntry *preventry = 0;
231	AliGenerator* gen = 0;
232
233	TObjArray *partArray = gAlice->GetMCApp()->Particles();
234
235	//
236	// Generate the vertex position used by all generators
237	//
238	if(fVertexSmear == kPerEvent) Vertex();
239	Bool_t PrimordialTrigger = kFALSE;
240
241	while(!PrimordialTrigger) {
242
243	//Reseting stack
244	AliRunLoader * runloader = gAlice->GetRunLoader();
245	if (runloader)
246	if (runloader->Stack())
247	runloader->Stack()->Reset();
248	//
249	// Loop over generators and generate events
250	Int_t igen=0;
251	Int_t npart =0;
252
253	while((entry = (AliGenCocktailEntry*)next())) {
254	gen = entry->Generator();
255	gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
256	if ( (npart = gRandom->Poisson(entry->Rate())) >0 ) {
257	igen++;
258	if (igen ==1) entry->SetFirst(0);
259	else entry->SetFirst((partArray->GetEntriesFast())+1);
260	gen->SetNumberParticles(npart);
261	gen->Generate();
262	entry->SetLast(partArray->GetEntriesFast());
263	preventry = entry;
264	}
265	}
266	next.Reset();
267
268	// Tesitng primordial trigger : Muon pair in the MUON spectrometer acceptance 171,178 and pTCut
269	Int_t iPart;
270	fNGenerated++;
271	Int_t numberOfMuons=0;
272	// printf(">>>fNGenerated is %d\n",fNGenerated);
273	for(iPart=0; iPart<partArray->GetEntriesFast(); iPart++){
274	// gAlice->GetMCApp()->Particle(iPart)->Print();
275	if ( (TMath::Abs(gAlice->GetMCApp()->Particle(iPart)->GetPdgCode())==13) &&
276	(gAlice->GetMCApp()->Particle(iPart)->Theta()*180./TMath::Pi()>fMuonThetaMinCut) &&
277	(gAlice->GetMCApp()->Particle(iPart)->Theta()*180./TMath::Pi()<fMuonThetaMaxCut) &&
278	(gAlice->GetMCApp()->Particle(iPart)->Pt()>fMuonPtCut) ) {
279	gAlice->GetMCApp()->Particle(iPart)->SetProductionVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2), 0.);
280	numberOfMuons++;
281	}
282	}
283	// printf(">>> Number of Muons is %d \n", numberOfMuons);
284	if (numberOfMuons >= fMuonMultiplicity ) PrimordialTrigger = kTRUE;
285	}
286	//printf(">>> Trigger Accepted!!! \n");
287	fNSucceded++;
288	// Float_t Ratio = (Float_t) fNSucceded/fNGenerated;
289	// printf("Generated %d, Succeded %d and Ratio %f\n",fNGenerated, fNSucceded,Ratio);
290	}
291
292
293
294
295
296