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

/**************************************************************************
 * Copyright(c) 1998-2007, 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.                  *
 **************************************************************************/


// Generator for muons according to kinematic parametrizations at ALICE
// (not at the surface).
// Origin: andrea.dainese@lnl.infn.it

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

#include "AliRun.h"
#include "AliESDtrack.h"
#include "AliESDVertex.h"
#include "AliGenCosmicsParam.h"

ClassImp(AliGenCosmicsParam)

//-----------------------------------------------------------------------------
AliGenCosmicsParam::AliGenCosmicsParam():
AliGenerator(),
fParamMI(kFALSE),
fParamACORDE(kFALSE),
fYOrigin(600.),
fMaxAngleWRTVertical(-99.),
fBkG(0.),
fTPC(kFALSE),
fITS(kFALSE),
fSPDouter(kFALSE),
fSPDinner(kFALSE),
fACORDE(kFALSE)
{
  //
  // Default constructor
  //
}
//-----------------------------------------------------------------------------
void AliGenCosmicsParam::Generate()
{
  //
  // Generate one muon
  //
  
  //
  Float_t origin[3];
  Float_t p[3];
  Int_t nt;
  Double_t ptot=0,pt=0,angleWRTVertical=0;
  Bool_t okMom=kFALSE,okAngle=kFALSE;
  //
  Float_t rtrigger=250.0,ztrigger=250.0; // TPC is default
  if(fITS)      { rtrigger=50.0; ztrigger=50.0; }
  if(fSPDouter) { rtrigger=6.5; ztrigger=14.0; }
  if(fSPDinner) { rtrigger=3.5; ztrigger=14.0; }


  // mu+ or mu-
  Float_t muMinusFraction = 4./9.; // mu+/mu- = 1.25
  Int_t ipart;
  if(gRandom->Rndm()<muMinusFraction) {
    ipart = 13; // mu-
  } else {
    ipart = -13; // mu+
  }

  if(fParamACORDE) { // extracted from AliGenACORDE events
    // sample total momentum only once (to speed up)
    TF1 *dNdpACORDE = new TF1("dNdpACORDE","x/(1.+(x/12.8)*(x/12.8))^1.96",fPMin,fPMax);
    ptot = (Double_t)dNdpACORDE->GetRandom();
    delete dNdpACORDE;
    dNdpACORDE = 0;
  }

  Int_t trials=0;

  while(1) {
    trials++;
    // origin
    origin[0]  = (fYOrigin*TMath::Tan(fMaxAngleWRTVertical)+rtrigger)*(-1.+2.*gRandom->Rndm());
    origin[1]  = fYOrigin;
    origin[2]  = (fYOrigin*TMath::Tan(fMaxAngleWRTVertical)+ztrigger)*(-1.+2.*gRandom->Rndm());

    // momentum
    while(1) {
      okMom=kFALSE; okAngle=kFALSE;

      if(fParamMI) { // parametrization by M.Ivanov of LEP cosmics data
	Float_t	pref  = 1. + gRandom->Exp(30.);
	p[1] = -pref; 
	p[0] = gRandom->Gaus(0.0,0.2)*pref;
	p[2] = gRandom->Gaus(0.0,0.2)*pref;
	if(gRandom->Rndm()>0.9) {
	  p[0] = gRandom->Gaus(0.0,0.4)*pref;
	  p[2] = gRandom->Gaus(0.0,0.4)*pref;
	}
	ptot=TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
	pt=TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
      } else if(fParamACORDE) { // extracted from AliGenACORDE events
	Float_t theta,phi;
	while(1) {
	  theta = gRandom->Gaus(0.5*TMath::Pi(),0.42);
	  if(TMath::Abs(theta-0.5*TMath::Pi())<fMaxAngleWRTVertical) break;
	}
	while(1) {
	  phi = gRandom->Gaus(-0.5*TMath::Pi(),0.42);
	  if(TMath::Abs(phi+0.5*TMath::Pi())<fMaxAngleWRTVertical) break;
	}
	pt = ptot*TMath::Sin(theta);
	p[0] = pt*TMath::Cos(phi); 
	p[1] = pt*TMath::Sin(phi); 
	p[2] = ptot*TMath::Cos(theta);
      } else {
	AliFatal("Parametrization not set: use SetParamMI or SetParamACORDE");
      }

      
      // check kinematic cuts
      if(TestBit(kMomentumRange)) {
	if(ptot>fPMin && ptot<fPMax) okMom=kTRUE;
      }

      angleWRTVertical=TMath::ACos(TMath::Abs(p[1])/ptot); // acos(|py|/ptot)
      if(angleWRTVertical<fMaxAngleWRTVertical) okAngle=kTRUE;
      
      if(okAngle&&okMom) break;
    }

    // acceptance trigger
    if(IntersectCylinder(rtrigger,ztrigger,ipart,origin,p)) {
      if(!fACORDE) {
	break; 
      } else {
	if(IntersectACORDE(ipart,origin,p)) break;
      }
    }
    //
  }

  Float_t polarization[3]= {0,0,0};
  PushTrack(fTrackIt,-1,ipart,p,origin,polarization,0,kPPrimary,nt);

  //printf("TRIALS %d\n",trials);
  

  return;
}
//-----------------------------------------------------------------------------
void AliGenCosmicsParam::Init()
{
  // 
  // Initialisation, check consistency of selected ranges
  //
  if(TestBit(kPtRange)) 
    AliFatal("You cannot set the pt range for this generator! Only momentum range");
  if(fPMin<8.) { 
    fPMin=8.; 
    if(TestBit(kMomentumRange)) 
      AliWarning("Minimum momentum cannot be < 8 GeV/c"); 
  }
  if(fMaxAngleWRTVertical<0.) 
    AliFatal("You must use SetMaxAngleWRTVertical() instead of SetThetaRange(), SetPhiRange()");

  printf("************ AliGenCosmicsParam ****************\n");
  printf("***** Muons generated at Y = %f cm\n",fYOrigin);
  printf("************************************************\n");

  return;
}
//-----------------------------------------------------------------------------
Bool_t AliGenCosmicsParam::IntersectCylinder(Float_t r,Float_t z,Int_t pdg,
					     Float_t o[3],Float_t p[3]) const
{
  //
  // Intersection between muon and cylinder [-z,+z] with radius r
  //

  Float_t en = TMath::Sqrt(0.105*0.105+p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
  TParticle part(pdg,0,0,0,0,0,p[0],p[1],p[2],en,o[0],o[1],o[2],0);
  AliESDtrack track(&part);
  Double_t pos[3]={0.,0.,0.},sigma[3]={0.,0.,0.};
  AliESDVertex origin(pos,sigma);

  track.RelateToVertex(&origin,fBkG,10000.);

  Float_t d0z0[2],covd0z0[3];
  track.GetImpactParameters(d0z0,covd0z0);

  // check rphi 
  if(TMath::Abs(d0z0[0])>r) return kFALSE;
  // check z
  if(TMath::Abs(d0z0[1])>z) return kFALSE;

  /*
    if(TMath::Abs(fB)<0.01) {  // NO FIELD
    Float_t drphi = TMath::Abs(o[1]-p[1]/p[0]*o[0])/
    TMath::Sqrt(p[1]*p[1]/p[0]/p[0]+1.);
    if(drphi>r) return kFALSE;
    Float_t dz = o[2]-p[2]/p[0]*o[0]+p[2]/p[0]*
    (p[1]*p[1]/p[0]/p[0]*o[0]-p[1]/p[0]*o[1])/(1.+p[1]*p[1]/p[0]/p[0]);
    if(TMath::Abs(dz)>z) return kFALSE;
    }
  */

  return kTRUE;
}
//-----------------------------------------------------------------------------
Bool_t AliGenCosmicsParam::IntersectACORDE(Int_t pdg,
					   Float_t o[3],Float_t p[3]) const
{
  //
  // Intersection between muon and ACORDE (very rough)
  //

  Float_t en = TMath::Sqrt(0.105*0.105+p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
  TParticle part(pdg,0,0,0,0,0,p[0],p[1],p[2],en,o[0],o[1],o[2],0);
  AliESDtrack track(&part);

  Float_t rACORDE=800.0,xACORDE=750.0,zACORDE=600.0;

  Double_t xyz[3];
  track.GetXYZAt(rACORDE,fBkG,xyz);

  // check global x 
  if(TMath::Abs(xyz[0]) > xACORDE) return kFALSE;
  // check global z
  if(TMath::Abs(xyz[2]) > zACORDE) return kFALSE;

  return kTRUE;
}
//-----------------------------------------------------------------------------
Commit	Line	Data
7b4a37c0	1	/**************************************************************************
	2	* Copyright(c) 1998-2007, 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
	17	// Generator for muons according to kinematic parametrizations at ALICE
	18	// (not at the surface).
	19	// Origin: andrea.dainese@lnl.infn.it
	20
	21	#include <TParticle.h>
	22	#include <TF1.h>
	23
	24	#include "AliRun.h"
	25	#include "AliESDtrack.h"
	26	#include "AliESDVertex.h"
	27	#include "AliGenCosmicsParam.h"
	28
	29	ClassImp(AliGenCosmicsParam)
	30
	31	//-----------------------------------------------------------------------------
	32	AliGenCosmicsParam::AliGenCosmicsParam():
	33	AliGenerator(),
	34	fParamMI(kFALSE),
	35	fParamACORDE(kFALSE),
	36	fYOrigin(600.),
	37	fMaxAngleWRTVertical(-99.),
	38	fBkG(0.),
	39	fTPC(kFALSE),
	40	fITS(kFALSE),
	41	fSPDouter(kFALSE),
fb2a6a7f	42	fSPDinner(kFALSE),
fb2a6a7f	43	fACORDE(kFALSE)
7b4a37c0	44	{
	45	//
	46	// Default constructor
	47	//
	48	}
	49	//-----------------------------------------------------------------------------
	50	void AliGenCosmicsParam::Generate()
	51	{
	52	//
	53	// Generate one muon
	54	//
	55
	56	//
	57	Float_t origin[3];
	58	Float_t p[3];
	59	Int_t nt;
	60	Double_t ptot=0,pt=0,angleWRTVertical=0;
	61	Bool_t okMom=kFALSE,okAngle=kFALSE;
	62	//
fb2a6a7f	63	Float_t rtrigger=250.0,ztrigger=250.0; // TPC is default
	64	if(fITS) { rtrigger=50.0; ztrigger=50.0; }
	65	if(fSPDouter) { rtrigger=6.5; ztrigger=14.0; }
	66	if(fSPDinner) { rtrigger=3.5; ztrigger=14.0; }
	67
7b4a37c0	68
7b4a37c0	69	// mu+ or mu-
fb2a6a7f	70	Float_t muMinusFraction = 4./9.; // mu+/mu- = 1.25
	71	Int_t ipart;
	72	if(gRandom->Rndm()<muMinusFraction) {
	73	ipart = 13; // mu-
	74	} else {
	75	ipart = -13; // mu+
	76	}
7b4a37c0	77
	78	if(fParamACORDE) { // extracted from AliGenACORDE events
	79	// sample total momentum only once (to speed up)
	80	TF1 dNdpACORDE = new TF1("dNdpACORDE","x/(1.+(x/12.8)(x/12.8))^1.96",fPMin,fPMax);
	81	ptot = (Double_t)dNdpACORDE->GetRandom();
	82	delete dNdpACORDE;
	83	dNdpACORDE = 0;
	84	}
	85
	86	Int_t trials=0;
	87
	88	while(1) {
	89	trials++;
	90	// origin
fb2a6a7f	91	origin[0] = (fYOriginTMath::Tan(fMaxAngleWRTVertical)+rtrigger)(-1.+2.*gRandom->Rndm());
7b4a37c0	92	origin[1] = fYOrigin;
fb2a6a7f	93	origin[2] = (fYOriginTMath::Tan(fMaxAngleWRTVertical)+ztrigger)(-1.+2.*gRandom->Rndm());
7b4a37c0	94
	95	// momentum
	96	while(1) {
	97	okMom=kFALSE; okAngle=kFALSE;
	98
	99	if(fParamMI) { // parametrization by M.Ivanov of LEP cosmics data
	100	Float_t pref = 1. + gRandom->Exp(30.);
	101	p[1] = -pref;
	102	p[0] = gRandom->Gaus(0.0,0.2)*pref;
	103	p[2] = gRandom->Gaus(0.0,0.2)*pref;
	104	if(gRandom->Rndm()>0.9) {
	105	p[0] = gRandom->Gaus(0.0,0.4)*pref;
	106	p[2] = gRandom->Gaus(0.0,0.4)*pref;
	107	}
	108	ptot=TMath::Sqrt(p[0]p[0]+p[1]p[1]+p[2]*p[2]);
	109	pt=TMath::Sqrt(p[0]p[0]+p[1]p[1]);
	110	} else if(fParamACORDE) { // extracted from AliGenACORDE events
	111	Float_t theta,phi;
	112	while(1) {
	113	theta = gRandom->Gaus(0.5*TMath::Pi(),0.42);
	114	if(TMath::Abs(theta-0.5*TMath::Pi())<fMaxAngleWRTVertical) break;
	115	}
	116	while(1) {
	117	phi = gRandom->Gaus(-0.5*TMath::Pi(),0.42);
	118	if(TMath::Abs(phi+0.5*TMath::Pi())<fMaxAngleWRTVertical) break;
	119	}
	120	pt = ptot*TMath::Sin(theta);
	121	p[0] = pt*TMath::Cos(phi);
	122	p[1] = pt*TMath::Sin(phi);
	123	p[2] = ptot*TMath::Cos(theta);
	124	} else {
	125	AliFatal("Parametrization not set: use SetParamMI or SetParamACORDE");
	126	}
	127
	128
	129	// check kinematic cuts
	130	if(TestBit(kMomentumRange)) {
	131	if(ptot>fPMin && ptot<fPMax) okMom=kTRUE;
	132	}
	133
	134	angleWRTVertical=TMath::ACos(TMath::Abs(p[1])/ptot); // acos(\|py\|/ptot)
	135	if(angleWRTVertical<fMaxAngleWRTVertical) okAngle=kTRUE;
	136
	137	if(okAngle&&okMom) break;
	138	}
	139
fb2a6a7f	140	// acceptance trigger
	141	if(IntersectCylinder(rtrigger,ztrigger,ipart,origin,p)) {
	142	if(!fACORDE) {
	143	break;
	144	} else {
	145	if(IntersectACORDE(ipart,origin,p)) break;
	146	}
	147	}
	148	//
7b4a37c0	149	}
	150
	151	Float_t polarization[3]= {0,0,0};
	152	PushTrack(fTrackIt,-1,ipart,p,origin,polarization,0,kPPrimary,nt);
	153
	154	//printf("TRIALS %d\n",trials);
	155
	156
	157	return;
	158	}
	159	//-----------------------------------------------------------------------------
	160	void AliGenCosmicsParam::Init()
	161	{
	162	//
	163	// Initialisation, check consistency of selected ranges
	164	//
	165	if(TestBit(kPtRange))
	166	AliFatal("You cannot set the pt range for this generator! Only momentum range");
	167	if(fPMin<8.) {
	168	fPMin=8.;
	169	if(TestBit(kMomentumRange))
	170	AliWarning("Minimum momentum cannot be < 8 GeV/c");
	171	}
	172	if(fMaxAngleWRTVertical<0.)
	173	AliFatal("You must use SetMaxAngleWRTVertical() instead of SetThetaRange(), SetPhiRange()");
	174
	175	printf("********** AliGenCosmicsParam **************\n");
	176	printf("***** Muons generated at Y = %f cm\n",fYOrigin);
	177	printf("************************************************\n");
	178
	179	return;
	180	}
	181	//-----------------------------------------------------------------------------
	182	Bool_t AliGenCosmicsParam::IntersectCylinder(Float_t r,Float_t z,Int_t pdg,
	183	Float_t o[3],Float_t p[3]) const
	184	{
	185	//
	186	// Intersection between muon and cylinder [-z,+z] with radius r
	187	//
	188
	189	Float_t en = TMath::Sqrt(0.1050.105+p[0]p[0]+p[1]p[1]+p[2]p[2]);
	190	TParticle part(pdg,0,0,0,0,0,p[0],p[1],p[2],en,o[0],o[1],o[2],0);
	191	AliESDtrack track(&part);
	192	Double_t pos[3]={0.,0.,0.},sigma[3]={0.,0.,0.};
	193	AliESDVertex origin(pos,sigma);
	194
	195	track.RelateToVertex(&origin,fBkG,10000.);
	196
	197	Float_t d0z0[2],covd0z0[3];
	198	track.GetImpactParameters(d0z0,covd0z0);
	199
	200	// check rphi
	201	if(TMath::Abs(d0z0[0])>r) return kFALSE;
	202	// check z
	203	if(TMath::Abs(d0z0[1])>z) return kFALSE;
	204
	205	/*
	206	if(TMath::Abs(fB)<0.01) { // NO FIELD
	207	Float_t drphi = TMath::Abs(o[1]-p[1]/p[0]*o[0])/
	208	TMath::Sqrt(p[1]*p[1]/p[0]/p[0]+1.);
	209	if(drphi>r) return kFALSE;
	210	Float_t dz = o[2]-p[2]/p[0]o[0]+p[2]/p[0]
	211	(p[1]p[1]/p[0]/p[0]o[0]-p[1]/p[0]o[1])/(1.+p[1]p[1]/p[0]/p[0]);
	212	if(TMath::Abs(dz)>z) return kFALSE;
213	}
214	*/
215
216	return kTRUE;
217	}
218	//-----------------------------------------------------------------------------
fb2a6a7f	219	Bool_t AliGenCosmicsParam::IntersectACORDE(Int_t pdg,
	220	Float_t o[3],Float_t p[3]) const
	221	{
	222	//
	223	// Intersection between muon and ACORDE (very rough)
	224	//
	225
	226	Float_t en = TMath::Sqrt(0.1050.105+p[0]p[0]+p[1]p[1]+p[2]p[2]);
	227	TParticle part(pdg,0,0,0,0,0,p[0],p[1],p[2],en,o[0],o[1],o[2],0);
	228	AliESDtrack track(&part);
	229
	230	Float_t rACORDE=800.0,xACORDE=750.0,zACORDE=600.0;
	231
	232	Double_t xyz[3];
	233	track.GetXYZAt(rACORDE,fBkG,xyz);
	234
	235	// check global x
	236	if(TMath::Abs(xyz[0]) > xACORDE) return kFALSE;
	237	// check global z
	238	if(TMath::Abs(xyz[2]) > zACORDE) return kFALSE;
	239
	240	return kTRUE;
	241	}
	242	//-----------------------------------------------------------------------------