newly implemented model for slow nucleon production

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

diff --git a/EVGEN/AliGenCosmicsParam.cxx b/EVGEN/AliGenCosmicsParam.cxx
index 07a29789b074b42d1b832a078c245428561016f9..a40f84ff2e20cf13ebf804d06f911d1ca50f2ddb 100644 (file)
--- a/EVGEN/AliGenCosmicsParam.cxx
+++ b/EVGEN/AliGenCosmicsParam.cxx
@@ -33,6 +33,7 @@ AliGenCosmicsParam::AliGenCosmicsParam():
 AliGenerator(),
 fParamMI(kFALSE),
 fParamACORDE(kFALSE),
+fParamDataTPC(kTRUE),
 fYOrigin(600.),
 fMaxAngleWRTVertical(-99.),
 fBkG(0.),
@@ -45,7 +46,8 @@ fSDDouter(kFALSE),
 fSSDinner(kFALSE),
 fSSDouter(kFALSE),
 fACORDE(kFALSE),
-fACORDE4ITS(kFALSE)
+fACORDE4ITS(kFALSE),
+fBottomScintillator(kFALSE)
 {
   //
   // Default constructor
@@ -98,6 +100,12 @@ void AliGenCosmicsParam::Generate()
       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;
     }
 
     while(1) {
@@ -122,7 +130,7 @@ void AliGenCosmicsParam::Generate()
          }
          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
+       } else if(fParamACORDE || fParamDataTPC) {
          Float_t theta,phi;
          while(1) {
            theta = gRandom->Gaus(0.5*TMath::Pi(),0.42);
@@ -137,7 +145,7 @@ void AliGenCosmicsParam::Generate()
          p[1] = pt*TMath::Sin(phi); 
          p[2] = ptot*TMath::Cos(theta);
        } else {
-         AliFatal("Parametrization not set: use SetParamMI or SetParamACORDE");
+         AliFatal("Parametrization not set: use SetParamDataTPC, SetParamMI, or SetParamACORDE");
        }
        
        
@@ -156,10 +164,15 @@ void AliGenCosmicsParam::Generate()
 
       // acceptance trigger
       if(IntersectCylinder(rtrigger,ztrigger,ipart,origin,p)) {
-       if(!fACORDE) {
-         break; 
-       } else {
+       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;
        }
       }
       //
@@ -181,10 +194,12 @@ void AliGenCosmicsParam::Init()
   //
   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()");
@@ -241,14 +256,19 @@ Bool_t AliGenCosmicsParam::IntersectACORDE(Int_t pdg,
   //
 
   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);
+  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=500.0;
+  Float_t rACORDE=800.0,xACORDE=750.0/*250.0*/,zACORDE=500.0;
   if(fACORDE4ITS) { xACORDE=100.0; zACORDE=100.0; }
 
-  Double_t xyz[3];
-  track.GetXYZAt(rACORDE,fBkG,xyz);
+  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;
@@ -258,3 +278,32 @@ Bool_t AliGenCosmicsParam::IntersectACORDE(Int_t pdg,
   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;
+}
+//-----------------------------------------------------------------------------