.so cleanup: removed from gSystem->Load()

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

diff --git a/EVGEN/AliGenParam.cxx b/EVGEN/AliGenParam.cxx
index 3e92b14c9c3e1f9db3d7060ddc42fcbf21288ec5..de448cb8368da514a0060095a2cd18fc5801694d 100644 (file)
--- a/EVGEN/AliGenParam.cxx
+++ b/EVGEN/AliGenParam.cxx
@@ -69,7 +69,9 @@ AliGenParam::AliGenParam()
        fDeltaPt(0.01),
        fSelectAll(kFALSE),
   fDecayer(0),
-  fForceConv(kFALSE)
+  fForceConv(kFALSE),
+  fKeepParent(kFALSE),
+  fKeepIfOneChildSelected(kFALSE)
 {
   // Default constructor
 }
@@ -93,7 +95,9 @@ AliGenParam::AliGenParam(Int_t npart, const AliGenLib * Library,  Int_t param, c
      fDeltaPt(0.01),
      fSelectAll(kFALSE),
    fDecayer(0),
-   fForceConv(kFALSE)
+   fForceConv(kFALSE),
+   fKeepParent(kFALSE),
+   fKeepIfOneChildSelected(kFALSE)
 {
   // Constructor using number of particles parameterisation id and library
     fName = "Param";
@@ -121,7 +125,9 @@ AliGenParam::AliGenParam(Int_t npart, Int_t param, const char* tname, const char
     fDeltaPt(0.01),
     fSelectAll(kFALSE),
   fDecayer(0),
-  fForceConv(kFALSE)
+  fForceConv(kFALSE),
+  fKeepParent(kFALSE),
+  fKeepIfOneChildSelected(kFALSE)
 {
   // Constructor using parameterisation id and number of particles
   //
@@ -170,7 +176,9 @@ AliGenParam::AliGenParam(Int_t npart, Int_t param,
      fDeltaPt(0.01),
      fSelectAll(kFALSE),
   fDecayer(0),
-  fForceConv(kFALSE)
+  fForceConv(kFALSE),
+  fKeepParent(kFALSE),
+  fKeepIfOneChildSelected(kFALSE)
 {
   // Constructor
   // Gines Martinez 1/10/99 
@@ -205,9 +213,9 @@ TVector3 AliGenParam::OrthogonalVector(TVector3 &inVec){
   int solvDim=0;
   double tmp=abc[0];
   for(int i=0; i<3; i++)
-    if(abs(abc[i])>tmp){
+    if(fabs(abc[i])>tmp){
       solvDim=i;
-      tmp=abs(abc[i]);
+      tmp=fabs(abc[i]);
     }
   xyz[solvDim]=(-abc[(1+solvDim)%3]-abc[(2+solvDim)%3])/abc[(0+solvDim)%3];
   
@@ -225,11 +233,6 @@ void AliGenParam::RotateVector(Double_t *pin, Double_t *pout, Double_t costheta,
   return;
 }
 
-Double_t AliGenParam::IntegratedKrollWada(Double_t mh){
-  if(mh<0.003) return 0;
-  return 2*log(mh/0.000511/exp(7.0/4.0))/411.0/TMath::Pi();
-}
-
 double AliGenParam::ScreenFunction1(double screenVariable){
   if(screenVariable>1)
     return 42.24 - 8.368 * log(screenVariable + 0.952);
@@ -321,8 +324,8 @@ Double_t AliGenParam::RandomMass(Double_t mh){
   while(true){
     double y=fRandom->Rndm();
     double mee=2*0.000511*TMath::Power(2*0.000511/mh,-y); //inverse of the enveloping cumulative distribution
-    double apxkw=2.0/3.0/137.036/TMath::Pi()/mee;
-    double val=fRandom->Uniform(0,apxkw);  //enveloping probability density0
+    double apxkw=2.0/3.0/137.036/TMath::Pi()/mee; //enveloping probability density
+    double val=fRandom->Uniform(0,apxkw);
     double kw=apxkw*sqrt(1-4*0.000511*0.000511/mee/mee)*(1+2*0.000511*0.000511/mee/mee)*1*1*TMath::Power(1-mee*mee/mh/mh,3);
     if(val<kw)
       return mee;
@@ -334,8 +337,8 @@ Int_t AliGenParam::VirtualGammaPairProduction(TClonesArray *particles, Int_t nPa
   Int_t nPartNew=nPart;
   for(int iPart=0; iPart<nPart; iPart++){      
     TParticle *gamma = (TParticle *) particles->At(iPart);
-    if(gamma->GetPdgCode()!=223000 || gamma->GetPdgCode()!=224000) continue;
-    //if(gamma->Energy()<0.001022) continue; //can never be
+    if(gamma->GetPdgCode()!=220000) continue;
+    if(gamma->Pt()<0.002941) continue;  //approximation of kw in AliGenEMlib is 0 below 0.002941
     double mass=RandomMass(gamma->Pt());
 
     // lepton pair kinematics in virtual photon rest frame 
@@ -378,6 +381,7 @@ Int_t AliGenParam::VirtualGammaPairProduction(TClonesArray *particles, Int_t nPa
     new((*particles)[nPartNew]) TParticle(-11, gamma->GetStatusCode(), iPart+1, -1, 0, 0, e2V4, vtx);
     nPartNew++;
   }
+  return nPartNew;
 }
 
 Int_t AliGenParam::ForceGammaConversion(TClonesArray *particles, Int_t nPart)
@@ -390,7 +394,6 @@ Int_t AliGenParam::ForceGammaConversion(TClonesArray *particles, Int_t nPart)
     TParticle *gamma = (TParticle *) particles->At(iPart);
     if(gamma->GetPdgCode()!=22) continue;
     if(gamma->Energy()<0.001022) continue;
-
     TVector3 gammaV3(gamma->Px(),gamma->Py(),gamma->Pz());
     double frac=RandomEnergyFraction(1,gamma->Energy());
     double Ee1=frac*gamma->Energy();
@@ -424,8 +427,7 @@ Int_t AliGenParam::ForceGammaConversion(TClonesArray *particles, Int_t nPart)
     new((*particles)[nPartNew]) TParticle(-11, gamma->GetStatusCode(), iPart+1, -1, 0, 0, TLorentzVector(e2V3,Ee2), vtx);
     nPartNew++;
   }
-  // particles->Compress();
-  return particles->GetEntriesFast();
+  return nPartNew;
 }
 
 //____________________________________________________________
@@ -433,7 +435,7 @@ void AliGenParam::Init()
 {
   // Initialisation
 
-    if (gMC) fDecayer = gMC->GetDecayer();
+    if (TVirtualMC::GetMC()) fDecayer = TVirtualMC::GetMC()->GetDecayer();
   //Begin_Html
   /*
     <img src="picts/AliGenParam.gif">
@@ -579,7 +581,7 @@ void AliGenParam::GenerateN(Int_t ntimes)
       Int_t iTemp = iPart;
 
       // custom pdg codes for to destinguish direct photons
-      if((iPart>=221000) && (iPart<=229000)) iPart=22;
+      if(iPart==220000) iPart=22;
 
          fChildWeight=(fDecayer->GetPartialBranchingRatio(iPart))*fParentWeight;          
          TParticlePDG *particle = pDataBase->GetParticle(iPart);
@@ -664,13 +666,12 @@ void AliGenParam::GenerateN(Int_t ntimes)
              Int_t np=fDecayer->ImportParticles(particles);
 
        iPart=iTemp;
-       if(fForceConv) np=ForceGammaConversion(particles,np);
-       if(iPart==223000 || iPart==224000){
-         // wgtp*=IntegratedKrollWada(pt);
-         // wgtch*=IntegratedKrollWada(pt);
-         // np=VirtualGammaPairProduction(particles,np)
+       if(iPart==220000){
+         TParticle *gamma = (TParticle *)particles->At(0);
+         gamma->SetPdgCode(iPart);
+         np=VirtualGammaPairProduction(particles,np);
        }
-
+       if(fForceConv) np=ForceGammaConversion(particles,np);
 
              //  Selecting  GeometryAcceptance for particles fPdgCodeParticleforAcceptanceCut;
              if (fGeometryAcceptance) 
@@ -736,8 +737,10 @@ void AliGenParam::GenerateN(Int_t ntimes)
                                  pSelected[i]  = 1;
                                  ncsel++;
                              } else {
-                                 ncsel=-1;
-                                 break;
+                                 if(!fKeepIfOneChildSelected){
+                                   ncsel=-1;
+                                   break;
+                                 }
                              } // child kine cuts
                          } else {
                              pSelected[i]  = 1;
@@ -748,8 +751,8 @@ void AliGenParam::GenerateN(Int_t ntimes)
              } // if decay products
              
              Int_t iparent;
-             if ((fCutOnChild && ncsel >0) || !fCutOnChild){
-                 ipa++;
+             
+             if (fKeepParent || (fCutOnChild && ncsel >0) || !fCutOnChild){
          //
          // Parent
                  
@@ -759,6 +762,11 @@ void AliGenParam::GenerateN(Int_t ntimes)
                  KeepTrack(nt); 
                  fNprimaries++;
                  
+                 //but count is as "generated" particle" only if it produced child(s) within cut
+                 if ((fCutOnChild && ncsel >0) || !fCutOnChild){
+                   ipa++;
+                 }
+                 
          //
          // Decay Products
          //