Introducing Id
[u/mrichter/AliRoot.git] / ZDC / AliGenZDC.cxx
index ac87c05..96b9f05 100644 (file)
 
 /*
 $Log$
+Revision 1.11  2001/11/11 17:27:50  hristov
+Minor corrections
+
+Revision 1.10  2001/09/26 16:00:47  coppedis
+Minor change
+
+Revision 1.8  2001/04/20 10:10:25  coppedis
+Minor changes
+
+Revision 1.7  2001/03/15 16:13:28  coppedis
+Code review
+
+Revision 1.6  2000/11/30 17:16:14  coppedis
+Changes suggested by fca
+
+Revision 1.5  2000/11/22 11:30:12  coppedis
+Major code revision
+
+Revision 1.4  2000/10/05 08:02:47  fca
+Correction of the generator direction
+
+Revision 1.3  2000/10/02 21:28:20  fca
+Removal of useless dependecies via forward declarations
+
+Revision 1.2  2000/07/11 11:12:34  fca
+Some syntax corrections for non standard HP aCC
+
 Revision 1.1  2000/07/10 13:58:01  fca
 New version of ZDC from E.Scomparin & C.Oppedisano
 
@@ -24,9 +51,13 @@ Revision 1.6  1999/09/29 09:24:35  fca
 Introduction of the Copyright and cvs Log
 
 */
+#include <assert.h>
+
 #include <TRandom.h>
 #include <TLorentzVector.h>
 #include <TVector3.h>
+#include "TDatabasePDG.h"
+#include "TMCProcess.h"
 
 #include "AliGenZDC.h"
 #include "AliConst.h"
@@ -59,20 +90,32 @@ AliGenZDC::AliGenZDC(Int_t npart)
   fCosy  = 0.;
   fCosz  = 1.;
   fPseudoRapidity = 0.;
+  
   fFermiflag = 1;
   // LHC values for beam divergence and crossing angle
   fBeamDiv = 0.000032;
   fBeamCrossAngle = 0.0001;
   fBeamCrossPlane = 2;
+  
+  Int_t i, j;
+  for(i=0; i<201; i++){
+     fProbintp[i] = 0;
+     fProbintn[i] = 0;
+  }
+  for(j=0; j<3; j++){
+     fPp[i] = 0;
+  }
+  fDebugOpt = 0;
 }
 
 //_____________________________________________________________________________
 void AliGenZDC::Init()
 {
-  printf("             Initializing AliGenZDC\n");
-  printf("     Fermi flag = %d, Beam Divergence = %f, Crossing Angle "
+  printf("\n\n         AliGenZDC initialized with:\n");
+  printf("   Fermi flag = %d, Beam Divergence = %f, Crossing Angle "
          "= %f, Crossing Plane = %d\n\n", fFermiflag, fBeamDiv, fBeamCrossAngle,
         fBeamCrossPlane);
+
   //Initialize Fermi momentum distributions for Pb-Pb
   FermiTwoGaussian(207.,82.,fPp,fProbintp,fProbintn);
 }  
@@ -85,11 +128,11 @@ void AliGenZDC::Generate()
   //
   Int_t i;
 
-  Double_t mass, pLab[3], balp0, balp[3], ddp[3], dddp0, dddp[3];
-  Float_t ptot = fPMin;
+  Double_t Mass, pLab[3], fP0, fP[3], fBoostP[3], ddp[3], dddp0, dddp[3]; 
+  Float_t  fPTrack[3], ptot = fPMin;
   Int_t nt;
   
-  if(fPseudoRapidity==0.){
+  if(fPseudoRapidity==0.){ 
     pLab[0] = ptot*fCosx;
     pLab[1] = ptot*fCosy;
     pLab[2] = ptot*fCosz;
@@ -104,59 +147,66 @@ void AliGenZDC::Generate()
      fP[i] = pLab[i];
   }
   
-  // Beam divergence and crossing angle
-  if(fBeamDiv!=0.) {BeamDivCross(0,fBeamDiv,fBeamCrossAngle,fBeamCrossPlane,pLab);}
-  if(fBeamCrossAngle!=0.) {BeamDivCross(1,fBeamDiv,fBeamCrossAngle,fBeamCrossPlane,pLab);}
   
+  // Beam divergence and crossing angle
+  if(fBeamCrossAngle!=0.) {
+    BeamDivCross(1,fBeamDiv,fBeamCrossAngle,fBeamCrossPlane,pLab);
+    for(i=0; i<=2; i++){
+       fP[i] = pLab[i];
+    }
+  }
+  if(fBeamDiv!=0.) {
+    BeamDivCross(0,fBeamDiv,fBeamCrossAngle,fBeamCrossPlane,pLab);
+    for(i=0; i<=2; i++){
+       fP[i] = pLab[i];
+    }
+  }
+
   // If required apply the Fermi momentum
   if(fFermiflag==1){
     if((fIpart==kProton) || (fIpart==kNeutron)){
       ExtractFermi(fIpart,fPp,fProbintp,fProbintn,ddp);
     }
-    if(fIpart==kProton) {mass = 0.93956563;}
-    if(fIpart==kNeutron) {mass = 0.93827231;}
-//  printf(" pLABx = %f  pLABy = %f  pLABz = %f \n",pLab[0],pLab[1],pLab[2]); 
-    for(i=0; i<=2; i++){
-       balp[i] = -pLab[i];
-    }
-    balp0 = TMath::Sqrt(pLab[0]*pLab[0]+pLab[1]*pLab[1]+pLab[2]*pLab[2]+mass*mass);
+    Mass=gAlice->PDGDB()->GetParticle(fIpart)->Mass();
+    fP0 = TMath::Sqrt(fP[0]*fP[0]+fP[1]*fP[1]+fP[2]*fP[2]+Mass*Mass);
     for(i=0; i<=2; i++){
        dddp[i] = ddp[i];
     }
-    dddp0 = TMath::Sqrt(dddp[0]*dddp[0]+dddp[1]*dddp[1]+dddp[2]*dddp[2]+mass*mass);
+    dddp0 = TMath::Sqrt(dddp[0]*dddp[0]+dddp[1]*dddp[1]+dddp[2]*dddp[2]+Mass*Mass);
     
-    TVector3 b(balp[0]/balp0, balp[1]/balp0, balp[2]/balp0);
+    TVector3 b(fP[0]/fP0, fP[1]/fP0, fP[2]/fP0);
     TLorentzVector pFermi(dddp[0], dddp[1], dddp[2], dddp0);
 
-//    printf(" pmu -> pLABx = %f  pLABy = %f  pLABz = %f  E = %f\n",
-//           balp[0],balp[1],balp[2],balp0); 
-//    printf(" Beta -> bx = %f  by = %f  bz = %f\n", b[0], b[1], b[2]);  
-//    printf(" pFermi -> px = %f, py = %f, pz = %f\n", pFermi[0], pFermi[1], pFermi[2]);
-    
     pFermi.Boost(b);
-
-//    printf(" Boosted momentum -> px = %f, py = %f, pz = %f\n",
-//          pFermi[0], pFermi[1], pFermi[2]);
     for(i=0; i<=2; i++){
        fBoostP[i] = pFermi[i];
+       fP[i] = pFermi[i];
     }
 
   }
-    
+  
+  for(i=0; i<=2; i++){
+     fPTrack[i] = fP[i];
+  }
+      
   Float_t polar[3] = {0,0,0};
-  gAlice->SetTrack(fTrackIt,-1,fIpart,fBoostP,fOrigin.GetArray(),polar,0,
-                  "Primary",nt);
+  gAlice->SetTrack(fTrackIt,-1,fIpart,fPTrack,fOrigin.GetArray(),polar,0,
+                  kPPrimary,nt);
+  if(fDebugOpt == 1){
+    printf("\n\n               Track momentum:\n");
+    printf("\n  fPTrack = %f, %f, %f \n",fPTrack[0],fPTrack[1],fPTrack[2]);
+  }
 }
 
 //_____________________________________________________________________________
-void AliGenZDC::FermiTwoGaussian(Double_t A, Float_t Z, Double_t* fPp, Double_t*
-               fProbintp, Double_t* fProbintn)
+void AliGenZDC::FermiTwoGaussian(Float_t A, Float_t Z, Double_t *fPp, 
+                Double_t *fProbintp, Double_t *fProbintn)
 {
 //
 // Momenta distributions according to the "double-gaussian"
 // distribution (Ilinov) - equal for protons and neutrons
 //
-//   printf("          Initialization of Fermi momenta distribution\n");
+
    fProbintp[0] = 0;
    fProbintn[0] = 0;
    Double_t sig1 = 0.113;
@@ -167,34 +217,39 @@ void AliGenZDC::FermiTwoGaussian(Double_t A, Float_t Z, Double_t* fPp, Double_t*
    for(Int_t i=1; i<=200; i++){
       Double_t p = i*0.005;
       fPp[i] = p;
-//      printf(" fPp[%d] = %f\n",i,fPp[i]);
       Double_t e1 = (p*p)/(2.*sig1*sig1);
       Double_t e2 = (p*p)/(2.*sig2*sig2);
       Double_t f1 = TMath::Exp(-(e1));
       Double_t f2 = TMath::Exp(-(e2));
       Double_t probp = xk*p*p*(f1/(TMath::Power(sig1,3.))+
                       alfa*f2/(TMath::Power(sig2,3.)))*0.005;
-//      printf("       probp = %f\n",probp);
       fProbintp[i] = fProbintp[i-1] + probp;
       fProbintn[i] = fProbintp[i];
-//      printf(" fProbintp[%d] = %f, fProbintp[%d] = %f\n",i,fProbintp[i],i,fProbintn[i]);
+   }
+   if(fDebugOpt == 1){
+     printf("\n\n              Initialization of Fermi momenta distribution \n");
+     for(Int_t i=0; i<=200; i++){
+        printf(" fProbintp[%d] = %f, fProbintn[%d] = %f\n",i,fProbintp[i],i,fProbintn[i]);
+     }
    }
 } 
 //_____________________________________________________________________________
-void AliGenZDC::ExtractFermi(Int_t id, Double_t* fPp, Double_t* fProbintp,
-                Double_t* fProbintn, Double_t* ddp)
+void AliGenZDC::ExtractFermi(Int_t id, Double_t *fPp, Double_t *fProbintp,
+                Double_t *fProbintn, Double_t *ddp)
 {
 //
 // Compute Fermi momentum for spectator nucleons
 //
+  
   Int_t i;
   Float_t xx = gRandom->Rndm();
+  assert ( id==kProton || id==kNeutron );
   if(id==kProton){
-    for(i=0; i<=200; i++){
+    for(i=1; i<=200; i++){
        if((xx>=fProbintp[i-1]) && (xx<fProbintp[i])) break;
        }
   }
-  else if(id==kNeutron){
+  else {
     for(i=0; i<=200; i++){
        if((xx>=fProbintn[i-1]) && (xx<fProbintn[i])) break;
        }
@@ -206,26 +261,27 @@ void AliGenZDC::ExtractFermi(Int_t id, Double_t* fPp, Double_t* fProbintp,
         ddp[0] = pext*TMath::Sin(tet)*TMath::Cos(phi);
         ddp[1] = pext*TMath::Sin(tet)*TMath::Sin(phi);
         ddp[2] = pext*cost;
+
+  if(fDebugOpt == 1){
+    printf("\n\n               Extraction of Fermi momentum\n");
+    printf("\n         pxFermi = %f  pyFermi = %f  pzFermi = %f \n",ddp[0],ddp[1],ddp[2]); 
+  }
 }
 
 //_____________________________________________________________________________
 void AliGenZDC::BeamDivCross(Int_t icross, Float_t fBeamDiv, Float_t fBeamCrossAngle, 
-                Int_t fBeamCrossPlane, Double_t* pLab)
+                Int_t fBeamCrossPlane, Double_t *pLab)
 {
-  Double_t tetpart, fipart, tetdiv, fidiv, angleSum[2], tetsum, fisum, dplab[3];
+  Double_t tetpart, fipart, tetdiv=0, fidiv=0, angleSum[2], tetsum, fisum;
   Double_t rvec;
 
   Int_t i;
-  
   Double_t pmq = 0.;
   for(i=0; i<=2; i++){
-     dplab[i] = pLab[i];
      pmq = pmq+pLab[i]*pLab[i];
   }
   Double_t pmod = TMath::Sqrt(pmq);
-//  printf("   pmod = %f\n",pmod);
 
-//  printf("   icross = %d, fBeamDiv = %f\n",icross,fBeamDiv);
   if(icross==0){
     rvec = gRandom->Gaus(0.0,1.0);
     tetdiv = fBeamDiv * TMath::Abs(rvec);
@@ -245,16 +301,15 @@ void AliGenZDC::BeamDivCross(Int_t icross, Float_t fBeamDiv, Float_t fBeamCrossA
       fidiv = k2PI/4.;
     }
   }
-//  printf("   tetdiv = %f, fidiv = %f\n",tetdiv,fidiv);
-  tetpart = TMath::ATan(TMath::Sqrt(dplab[0]*dplab[0]+dplab[1]*dplab[1])/dplab[2]);
-  if(dplab[1]!=0. || dplab[0]!=0.){
-    fipart = TMath::ATan2(dplab[1],dplab[0]);
+
+  tetpart = TMath::ATan(TMath::Sqrt(pLab[0]*pLab[0]+pLab[1]*pLab[1])/pLab[2]);
+  if(pLab[1]!=0. || pLab[0]!=0.){
+    fipart = TMath::ATan2(pLab[1],pLab[0]);
   }
   else{
     fipart = 0.;
   }
   if(fipart<0.) {fipart = fipart+k2PI;}
-//  printf("   tetpart = %f, fipart = %f\n",tetpart,fipart);
   tetdiv = tetdiv*kRaddeg;
   fidiv = fidiv*kRaddeg;
   tetpart = tetpart*kRaddeg;
@@ -262,22 +317,22 @@ void AliGenZDC::BeamDivCross(Int_t icross, Float_t fBeamDiv, Float_t fBeamCrossA
   AddAngle(tetpart,fipart,tetdiv,fidiv,angleSum);
   tetsum = angleSum[0];
   fisum  = angleSum[1];
-//  printf("   tetsum = %f, fisum = %f\n",tetsum,fisum);
   tetsum = tetsum*kDegrad;
   fisum = fisum*kDegrad;
   pLab[0] = pmod*TMath::Sin(tetsum)*TMath::Cos(fisum);
   pLab[1] = pmod*TMath::Sin(tetsum)*TMath::Sin(fisum);
   pLab[2] = pmod*TMath::Cos(tetsum);
-//  printf("   pLab[0] = %f pLab[1] = %f pLab[2] = %f \n\n",
-//         pLab[0],pLab[1],pLab[2]);
-  for(i=0; i<=2; i++){
-     fDivP[i] = pLab[i];
+  if(fDebugOpt == 1){
+    printf("\n\n               Beam divergence and crossing angle\n");
+    for(i=0; i<=2; i++){
+       printf("        pLab[%d] = %f\n",i,pLab[i]);
+    }
   }
 }
   
 //_____________________________________________________________________________
 void  AliGenZDC::AddAngle(Double_t theta1, Double_t phi1, Double_t theta2,
-              Double_t phi2, Double_t* angleSum)
+              Double_t phi2, Double_t *angleSum)
 {
   Double_t temp, conv, cx, cy, cz, ct1, st1, ct2, st2, cp1, sp1, cp2, sp2;
   Double_t rtetsum, tetsum, fisum;
@@ -308,7 +363,6 @@ void  AliGenZDC::AddAngle(Double_t theta1, Double_t phi1, Double_t theta2,
   if(temp<-1.) temp=-1.;
   fisum = conv*TMath::ACos(temp);
   if(cy<0) {fisum = 360.-fisum;}
-//  printf("   AddAngle -> tetsum = %f, fisum = %f\n",tetsum, fisum); 
   angleSum[0] = tetsum;
   angleSum[1] = fisum;
 }