Minor corrections

[u/mrichter/AliRoot.git] / ZDC / AliGenZDC.cxx

diff --git a/ZDC/AliGenZDC.cxx b/ZDC/AliGenZDC.cxx
index c17bebe9cc7f24ec0e7a93464177573b24bcc7a5..13ff20efdb42ec00e0a19f763f69f42c301ee854 100644 (file)
--- a/ZDC/AliGenZDC.cxx
+++ b/ZDC/AliGenZDC.cxx
@@ -15,20 +15,51 @@
 
 /*
 $Log$
 
 /*
 $Log$

+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
+

 Revision 1.7  2000/01/19 17:17:40  fca
 
 Revision 1.6  1999/09/29 09:24:35  fca
 Introduction of the Copyright and cvs Log
 
 */
 Revision 1.7  2000/01/19 17:17:40  fca
 
 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 <TRandom.h>
 #include <TLorentzVector.h>
 #include <TVector3.h>

+#include "TDatabasePDG.h"

 
 #include "AliGenZDC.h"
 #include "AliConst.h"
 #include "AliPDG.h"
 #include "AliRun.h"
 
 #include "AliGenZDC.h"
 #include "AliConst.h"
 #include "AliPDG.h"
 #include "AliRun.h"

+#include "AliMCProcess.h"

  
 ClassImp(AliGenZDC)
  
  
 ClassImp(AliGenZDC)
  


@@ -56,20 +87,32 @@ AliGenZDC::AliGenZDC(Int_t npart)
   fCosy  = 0.;
   fCosz  = 1.;
   fPseudoRapidity = 0.;
   fCosy  = 0.;
   fCosz  = 1.;
   fPseudoRapidity = 0.;

+  

   fFermiflag = 1;
   // LHC values for beam divergence and crossing angle
   fBeamDiv = 0.000032;
   fBeamCrossAngle = 0.0001;
   fBeamCrossPlane = 2;
   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()
 {
 }
 
 //_____________________________________________________________________________
 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);
          "= %f, Crossing Plane = %d\n\n", fFermiflag, fBeamDiv, fBeamCrossAngle,
         fBeamCrossPlane);

+

   //Initialize Fermi momentum distributions for Pb-Pb
   FermiTwoGaussian(207.,82.,fPp,fProbintp,fProbintn);
 }  
   //Initialize Fermi momentum distributions for Pb-Pb
   FermiTwoGaussian(207.,82.,fPp,fProbintp,fProbintn);
 }  


@@ -80,11 +123,13 @@ void AliGenZDC::Generate()
   //
   // Generate one trigger (n or p)
   //
   //
   // Generate one trigger (n or p)
   //

-  Double_t mass, pLab[3], balp0, balp[3], ddp[3], dddp0, dddp[3];
-  Float_t ptot = fPMin;
+  Int_t i;
+
+  Double_t Mass, pLab[3], fP0, fP[3], fBoostP[3], ddp[3], dddp0, dddp[3]; 
+  Float_t  fPTrack[3], ptot = fPMin;

   Int_t nt;
   
   Int_t nt;
   

-  if(fPseudoRapidity==0.){
+  if(fPseudoRapidity==0.){ 

     pLab[0] = ptot*fCosx;
     pLab[1] = ptot*fCosy;
     pLab[2] = ptot*fCosz;
     pLab[0] = ptot*fCosx;
     pLab[1] = ptot*fCosy;
     pLab[2] = ptot*fCosz;


@@ -95,63 +140,70 @@ void AliGenZDC::Generate()
     pLab[1] = 0.;
     pLab[2] = ptot*TMath::Cos(scang);
   }
     pLab[1] = 0.;
     pLab[2] = ptot*TMath::Cos(scang);
   }

-  for(Int_t i=0; i<=2; i++){
+  for(i=0; i<=2; i++){

      fP[i] = pLab[i];
   }
   
      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 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(Int_t 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);
-    for(Int_t i=0; i<=2; i++){
+    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];
     }
        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);
 
     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);
     pFermi.Boost(b);

-
-//    printf(" Boosted momentum -> px = %f, py = %f, pz = %f\n",
-//          pFermi[0], pFermi[1], pFermi[2]);
-    for(Int_t i=0; i<=2; i++){
+    for(i=0; i<=2; i++){

        fBoostP[i] = pFermi[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};
   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
 //
 {
 //
 // 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;
    fProbintp[0] = 0;
    fProbintn[0] = 0;
    Double_t sig1 = 0.113;


@@ -162,34 +214,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;
    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;
       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];
       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
 //
 {
 //
 // Compute Fermi momentum for spectator nucleons
 //

+  

   Int_t i;
   Float_t xx = gRandom->Rndm();
   Int_t i;
   Float_t xx = gRandom->Rndm();

+  assert ( id==kProton || id==kNeutron );

   if(id==kProton){
   if(id==kProton){

-    for(i=0; i<=200; i++){
+    for(i=1; i<=200; i++){

        if((xx>=fProbintp[i-1]) && (xx<fProbintp[i])) break;
        }
   }
        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;
        }
     for(i=0; i<=200; i++){
        if((xx>=fProbintn[i-1]) && (xx<fProbintn[i])) break;
        }


@@ -201,24 +258,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;
         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, 
 }
 
 //_____________________________________________________________________________
 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;
   Double_t rvec;

-  
+
+  Int_t i;

   Double_t pmq = 0.;
   Double_t pmq = 0.;

-  for(int i=0; i<=2; i++){
-     dplab[i] = pLab[i];
+  for(i=0; i<=2; i++){

      pmq = pmq+pLab[i]*pLab[i];
   }
   Double_t pmod = TMath::Sqrt(pmq);
      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);
   if(icross==0){
     rvec = gRandom->Gaus(0.0,1.0);
     tetdiv = fBeamDiv * TMath::Abs(rvec);


@@ -238,16 +298,15 @@ void AliGenZDC::BeamDivCross(Int_t icross, Float_t fBeamDiv, Float_t fBeamCrossA
       fidiv = k2PI/4.;
     }
   }
       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;}
   }
   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;
   tetdiv = tetdiv*kRaddeg;
   fidiv = fidiv*kRaddeg;
   tetpart = tetpart*kRaddeg;


@@ -255,22 +314,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];
   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);
   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(Int_t 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,
   }
 }
   
 //_____________________________________________________________________________
 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;
 {
   Double_t temp, conv, cx, cy, cz, ct1, st1, ct2, st2, cp1, sp1, cp2, sp2;
   Double_t rtetsum, tetsum, fisum;


@@ -301,7 +360,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;}
   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;
 }  
   angleSum[0] = tetsum;
   angleSum[1] = fisum;
 }