introducing SDD, SSD layer misal (Andrea Dainese)
[u/mrichter/AliRoot.git] / EVGEN / AliGenCosmicsParam.cxx
index a8f506d..3622700 100644 (file)
@@ -38,12 +38,20 @@ 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()
@@ -59,87 +67,116 @@ 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;
+    }
 
-    // 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) { // 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");
        }
-       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;
 }
 //-----------------------------------------------------------------------------
@@ -202,3 +239,62 @@ Bool_t AliGenCosmicsParam::IntersectCylinder(Float_t r,Float_t z,Int_t pdg,
   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;
+}
+//-----------------------------------------------------------------------------