AliGenerator(),
fParamMI(kFALSE),
fParamACORDE(kFALSE),
+fParamDataTPC(kTRUE),
fYOrigin(600.),
fMaxAngleWRTVertical(-99.),
fBkG(0.),
fTPC(kFALSE),
fITS(kFALSE),
+fSPDinner(kFALSE),
fSPDouter(kFALSE),
-fSPDinner(kFALSE)
+fSDDinner(kFALSE),
+fSDDouter(kFALSE),
+fSSDinner(kFALSE),
+fSSDouter(kFALSE),
+fACORDE(kFALSE),
+fACORDE4ITS(kFALSE),
+fBottomScintillator(kFALSE)
{
//
// Default constructor
//
+ SetNumberParticles(1);
}
//-----------------------------------------------------------------------------
void AliGenCosmicsParam::Generate()
Double_t ptot=0,pt=0,angleWRTVertical=0;
Bool_t okMom=kFALSE,okAngle=kFALSE;
//
+ Float_t rtrigger=1000.0,ztrigger=600.0;
+ if(fTPC) { rtrigger=250.0; ztrigger=250.0; }
+ if(fITS) { rtrigger=50.0; ztrigger=50.0; }
+ if(fSPDinner) { rtrigger=3.5; ztrigger=14.0; }
+ if(fSPDouter) { rtrigger=6.5; ztrigger=14.0; }
+ if(fSDDinner) { rtrigger=14.0; ztrigger=21.0; }
+ if(fSDDouter) { rtrigger=23.0; ztrigger=29.0; }
+ if(fSSDinner) { rtrigger=37.0; ztrigger=42.0; }
+ if(fSSDouter) { rtrigger=42.0; ztrigger=48.0; }
+
// mu+ or mu-
- Int_t ipart=13;
- if(gRandom->Rndm()<0.5) ipart *= -1;
-
- 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;
- }
+ Float_t muMinusFraction = 4./9.; // mu+/mu- = 1.25
+ Int_t ipart;
Int_t trials=0;
+ Int_t npart=0;
+
+ while (npart<fNpart) {
- while(1) {
- trials++;
- // origin
- origin[0] = fYOrigin*TMath::Tan(fMaxAngleWRTVertical)*(-1.+2.*gRandom->Rndm());
- origin[1] = fYOrigin;
- origin[2] = fYOrigin*TMath::Tan(fMaxAngleWRTVertical)*(-1.+2.*gRandom->Rndm());
+ 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;
+ } else if(fParamDataTPC) { // extracted from cosmics in TPC (Summer 08)
+ // sample total momentum only once (to speed up)
+ TF1 *dNdpTPC = new TF1("dNdpTPC","x/(1.+(x/3.)*(x/3.))^1.",fPMin,fPMax);
+ ptot = (Double_t)dNdpTPC->GetRandom();
+ delete dNdpTPC;
+ dNdpTPC = 0;
+ }
- // 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;
+ 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 || fParamDataTPC) {
+ 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 SetParamDataTPC, SetParamMI, or SetParamACORDE");
}
- while(1) {
- phi = gRandom->Gaus(-0.5*TMath::Pi(),0.42);
- if(TMath::Abs(phi+0.5*TMath::Pi())<fMaxAngleWRTVertical) break;
+
+
+ // check kinematic cuts
+ if(TestBit(kMomentumRange)) {
+ if(ptot>fPMin && ptot<fPMax) okMom=kTRUE;
+ } else {
+ okMom=kTRUE;
}
- 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;
}
- 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 && !fBottomScintillator) {
+ if(IntersectACORDE(ipart,origin,p)) break;
+ } else if(!fACORDE && fBottomScintillator) {
+ if(IntersectBottomScintillator(ipart,origin,p)) break;
+ } else if(fACORDE && fBottomScintillator) {
+ if(IntersectACORDE(ipart,origin,p) &&
+ IntersectBottomScintillator(ipart,origin,p)) break;
+ } else { // !fACORDE && !fBottomScintillator
+ break;
+ }
+ }
+ //
}
- // acceptance
- if(!fTPC && !fITS && !fSPDinner && !fSPDouter) break;
- if(fTPC) if(IntersectCylinder(250.,250.,ipart,origin,p)) break;
- if(fITS) if(IntersectCylinder(50.,50.,ipart,origin,p)) break;
- if(fSPDouter) if(IntersectCylinder(6.5,14.,ipart,origin,p)) break;
- if(fSPDinner) if(IntersectCylinder(3.5,14.0,ipart,origin,p)) break;
+ Float_t polarization[3]= {0,0,0};
+ PushTrack(fTrackIt,-1,ipart,p,origin,polarization,0,kPPrimary,nt);
+ npart++;
+ //printf("TRIALS %d\n",trials);
}
- Float_t polarization[3]= {0,0,0};
- PushTrack(fTrackIt,-1,ipart,p,origin,polarization,0,kPPrimary,nt);
-
- //printf("TRIALS %d\n",trials);
-
-
return;
}
//-----------------------------------------------------------------------------
//
if(TestBit(kPtRange))
AliFatal("You cannot set the pt range for this generator! Only momentum range");
- if(fPMin<8.) {
- fPMin=8.;
+ Double_t pmin=8.; // fParamACORDE
+ if(fParamDataTPC) pmin=0.5;
+ if(fPMin<pmin) {
+ fPMin=pmin;
if(TestBit(kMomentumRange))
- AliWarning("Minimum momentum cannot be < 8 GeV/c");
+ AliWarning(Form("Minimum momentum cannot be < %f GeV/c",pmin));
}
if(fMaxAngleWRTVertical<0.)
AliFatal("You must use SetMaxAngleWRTVertical() instead of SetThetaRange(), SetPhiRange()");
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/*250.0*/,zACORDE=500.0;
+ if(fACORDE4ITS) { xACORDE=100.0; zACORDE=100.0; }
+
+ Double_t planepoint[3]={0.,rACORDE,0.};
+ Double_t planenorm[3]={0.,1.,0.};
+
+ if(!track.Intersect(planepoint,planenorm,fBkG)) return kFALSE;
+
+ Double_t xyz[3]={planepoint[0],planepoint[1],planepoint[2]};
+ //printf("XYZ = %f %f %f\n",xyz[0],xyz[1],xyz[2]);
+
+ // check global x
+ if(TMath::Abs(xyz[0]) > xACORDE) return kFALSE;
+ // check global z
+ if(TMath::Abs(xyz[2]) > zACORDE) return kFALSE;
+
+ return kTRUE;
+}
+//-----------------------------------------------------------------------------
+Bool_t AliGenCosmicsParam::IntersectBottomScintillator(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);
+
+ Double_t xSc=40.,ySc=-350.,zSc=40.;
+
+ Double_t planepoint[3]={0.,ySc,0.};
+ Double_t planenorm[3]={0.,1.,0.};
+
+ if(!track.Intersect(planepoint,planenorm,fBkG)) return kFALSE;
+
+ Double_t xyz[3]={planepoint[0],planepoint[1],planepoint[2]};
+ //printf("XYZ = %f %f %f\n",xyz[0],xyz[1],xyz[2]);
+
+ // check global x
+ if(TMath::Abs(xyz[0]) > xSc) return kFALSE;
+ // check global z
+ if(TMath::Abs(xyz[2]) > zSc) return kFALSE;
+
+ return kTRUE;
+}
+//-----------------------------------------------------------------------------
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff