]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EVGEN/AliGenParam.cxx
When calculating a*a-b*b the form (a-b)*(a+b) is usually more numerically stable.
[u/mrichter/AliRoot.git] / EVGEN / AliGenParam.cxx
index 0fad3ae8c77f26765801fa8f1d676d37e98ca937..506347f5e663cc4ff4ca7a89faa14dd57545dd0d 100644 (file)
 // Decays are performed using Pythia.
 // andreas.morsch@cern.ch
 
-#include "AliGenParam.h"
-#include "AliDecayer.h"
-#include "AliGenMUONlib.h"
-#include "AliRun.h"
-#include <TParticle.h>
-#include <TParticlePDG.h>
+#include <TCanvas.h>
+#include <TClonesArray.h>
 #include <TDatabasePDG.h>
+#include <TF1.h>
+#include <TH1F.h>
 #include <TLorentzVector.h>
+#include <TMath.h>
+#include <TParticle.h>
+#include <TParticlePDG.h>
+#include <TROOT.h>
 #include <TVirtualMC.h>
 
-#include <TF1.h>
-#include <TCanvas.h>
-#include <TH1.h>
-#include <TMath.h>
+#include "AliDecayer.h"
+#include "AliGenMUONlib.h"
+#include "AliGenParam.h"
 #include "AliMC.h"
+#include "AliRun.h"
+#include "AliGenEventHeader.h"
 
 ClassImp(AliGenParam)
 
@@ -48,53 +51,73 @@ ClassImp(AliGenParam)
   //End_Html
 
 //____________________________________________________________
-AliGenParam::AliGenParam()
+    AliGenParam::AliGenParam():
+       fPtParaFunc(0),
+       fYParaFunc(0),
+       fIpParaFunc(0),
+       fPtPara(0),
+       fYPara(0),
+       fParam(0),
+       fdNdy0(0.),
+       fYWgt(0.),
+       fPtWgt(0.),
+       fBias(0.),
+       fTrials(0),
+       fDeltaPt(0.01),
+       fDecayer(0)
 {
-// Deafault constructor
-    fPtPara = 0;
-    fYPara  = 0;
-    fParam  = 0;
-    fAnalog = kAnalog;
-    SetDeltaPt();
-    fDecayer = 0;
-
-
+// Default constructor
 }
 //____________________________________________________________
-AliGenParam::AliGenParam(Int_t npart, AliGenLib * Library,  Int_t param, char* tname):AliGenMC(npart)
+AliGenParam::AliGenParam(Int_t npart, AliGenLib * Library,  Int_t param, char* tname)
+    :AliGenMC(npart),
+     fPtParaFunc(Library->GetPt(param, tname)),
+     fYParaFunc (Library->GetY (param, tname)),
+     fIpParaFunc(Library->GetIp(param, tname)),
+     fPtPara(0),
+     fYPara(0),
+     fParam(param),
+     fdNdy0(0.),
+     fYWgt(0.),
+     fPtWgt(0.),
+     fBias(0.),
+     fTrials(0),
+     fDeltaPt(0.01),
+     fDecayer(0)
 {
 // Constructor using number of particles parameterisation id and library
     fName = "Param";
     fTitle= "Particle Generator using pT and y parameterisation";
-    
-    fPtParaFunc = Library->GetPt(param, tname);
-    fYParaFunc  = Library->GetY (param, tname);
-    fIpParaFunc = Library->GetIp(param, tname);
-    
-    fPtPara = 0;
-    fYPara  = 0;
-    fParam  = param;
     fAnalog = kAnalog;
     SetForceDecay();
-    SetDeltaPt(); 
 }
 //____________________________________________________________
-AliGenParam::AliGenParam(Int_t npart, Int_t param, const char* tname, const char* name):AliGenMC(npart)
+AliGenParam::AliGenParam(Int_t npart, Int_t param, const char* tname, const char* name):
+    AliGenMC(npart),
+    fPtParaFunc(0),
+    fYParaFunc (0),
+    fIpParaFunc(0),
+    fPtPara(0),
+    fYPara(0),
+    fParam(param),
+    fdNdy0(0.),
+    fYWgt(0.),
+    fPtWgt(0.),
+    fBias(0.),
+    fTrials(0),
+    fDeltaPt(0.01),
+    fDecayer(0)
 {
 // Constructor using parameterisation id and number of particles
 //
-  fName = name;
-  fTitle= "Particle Generator using pT and y parameterisation";
+    fName = name;
+    fTitle= "Particle Generator using pT and y parameterisation";
       
     AliGenLib* pLibrary = new AliGenMUONlib();
     fPtParaFunc = pLibrary->GetPt(param, tname);
     fYParaFunc  = pLibrary->GetY (param, tname);
     fIpParaFunc = pLibrary->GetIp(param, tname);
     
-    fPtPara = 0;
-    fYPara  = 0;
-    fParam  = param;
     fAnalog = kAnalog;
     fChildSelect.Set(5);
     for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
@@ -104,28 +127,34 @@ AliGenParam::AliGenParam(Int_t npart, Int_t param, const char* tname, const char
     SetChildPtRange();
     SetChildPhiRange();
     SetChildThetaRange(); 
-    SetDeltaPt(); 
 }
 //____________________________________________________________
 
 AliGenParam::AliGenParam(Int_t npart, Int_t param,
-                         Double_t (*PtPara) (Double_t*, Double_t*),
-                         Double_t (*YPara ) (Double_t* ,Double_t*),
+                         Double_t (*PtPara) (const Double_t*, const Double_t*),
+                         Double_t (*YPara ) (const Double_t* ,const Double_t*),
                         Int_t    (*IpPara) (TRandom *))                 
-    :AliGenMC(npart)
+    :AliGenMC(npart),
+     
+     fPtParaFunc(PtPara),
+     fYParaFunc(YPara),
+     fIpParaFunc(IpPara),
+     fPtPara(0),
+     fYPara(0),
+     fParam(param),
+     fdNdy0(0.),
+     fYWgt(0.),
+     fPtWgt(0.),
+     fBias(0.),
+     fTrials(0),
+     fDeltaPt(0.01),
+     fDecayer(0)
 {
 // Constructor
 // Gines Martinez 1/10/99 
     fName = "Param";
     fTitle= "Particle Generator using pT and y parameterisation";
 
-    fPtParaFunc = PtPara; 
-    fYParaFunc  = YPara;  
-    fIpParaFunc = IpPara;
-//  
-    fPtPara = 0;
-    fYPara  = 0;
-    fParam  = param;
     fAnalog = kAnalog;
     fChildSelect.Set(5);
     for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
@@ -135,15 +164,6 @@ AliGenParam::AliGenParam(Int_t npart, Int_t param,
     SetChildPtRange();
     SetChildPhiRange();
     SetChildThetaRange();  
-    SetDeltaPt();
-}
-
-
-AliGenParam::AliGenParam(const AliGenParam & Param)
-    :AliGenMC(Param)
-{
-// Copy constructor
-    Param.Copy(*this);
 }
 
 //____________________________________________________________
@@ -170,6 +190,7 @@ void AliGenParam::Init()
     
     if (fPtPara) fPtPara->Delete();
     fPtPara = new TF1(name, fPtParaFunc, fPtMin, fPtMax,0);
+    gROOT->GetListOfFunctions()->Remove(fPtPara);
 //  Set representation precision to 10 MeV
     Int_t npx= Int_t((fPtMax - fPtMin) / fDeltaPt);
     
@@ -178,6 +199,8 @@ void AliGenParam::Init()
     sprintf(name, "y-parameterisation  for %s", GetName());
     if (fYPara) fYPara->Delete();
     fYPara  = new TF1(name, fYParaFunc, fYMin, fYMax, 0);
+    gROOT->GetListOfFunctions()->Remove(fYPara);
+
     
     sprintf(name, "pt-for-%s", GetName());
     TF1 ptPara(name ,fPtParaFunc, 0, 15, 0);
@@ -192,9 +215,9 @@ void AliGenParam::Init()
     fdNdy0=fYParaFunc(&y1,&y2);
 //
 // Integral over generation region
-    Float_t intYS  = yPara.Integral(fYMin, fYMax);
-    Float_t intPt0 = ptPara.Integral(0,15);
-    Float_t intPtS = ptPara.Integral(fPtMin,fPtMax);
+    Float_t intYS  = yPara.Integral(fYMin, fYMax,(Double_t*) 0x0,1.e-6);
+    Float_t intPt0 = ptPara.Integral(0,15,(Double_t *) 0x0,1.e-6);
+    Float_t intPtS = ptPara.Integral(fPtMin,fPtMax,(Double_t*) 0x0,1.e-6);
     Float_t phiWgt=(fPhiMax-fPhiMin)/2./TMath::Pi();
     if (fAnalog == kAnalog) {
        fYWgt  = intYS/fdNdy0;
@@ -261,6 +284,8 @@ void AliGenParam::Generate()
   Int_t ipa=0;
   
 // Generating fNpart particles
+  fNprimaries = 0;
+  
   while (ipa<fNpart) {
       while(1) {
 //
@@ -296,7 +321,7 @@ void AliGenParam::Generate()
                    "Division by 0: Please check you rapidity range !");
          }
          
-         pl=xmt*ty/sqrt(1.-ty*ty);
+         pl=xmt*ty/sqrt((1.-ty)*(1.+ty));
          theta=TMath::ATan2(pt,pl);
 // Cut on theta
          if(theta<fThetaMin || theta>fThetaMax) continue;
@@ -323,6 +348,8 @@ void AliGenParam::Generate()
 // if fForceDecay == none Primary particle is tracked by GEANT 
 // (In the latest, make sure that GEANT actually does all the decays you want)   
 //
+         Bool_t decayed = kFALSE;
+         
 
          if (fForceDecay != kNoDecay) {
 // Using lujet to decay particle
@@ -349,6 +376,7 @@ void AliGenParam::Generate()
              }
              
              if (np >1) {
+                 decayed = kTRUE;
                  TParticle* iparticle =  (TParticle *) particles->At(0);
                  Int_t ipF, ipL;
                  for (i = 1; i<np ; i++) {
@@ -385,7 +413,7 @@ void AliGenParam::Generate()
 //
 // children
                      
-                     if (ChildSelected(TMath::Abs(kf)) || fForceDecay == kAll && trackIt[i])
+                     if ((ChildSelected(TMath::Abs(kf)) || fForceDecay == kAll) && trackIt[i])
                      {
                          if (fCutOnChild) {
                              pc[0]=iparticle->Px();
@@ -412,17 +440,22 @@ void AliGenParam::Generate()
                  ipa++;
 //
 // Parent
-                 PushTrack(0, -1, iPart, p, origin0, polar, 0, kPPrimary, nt, wgtp);
+                 
+                 
+                 PushTrack(0, -1, iPart, p, origin0, polar, 0, kPPrimary, nt, wgtp, ((decayed)? 11 : 1));
                  pParent[0] = nt;
                  KeepTrack(nt); 
+                 fNprimaries++;
+                 
 //
 // Decay Products
 //               
                  for (i = 1; i < np; i++) {
                      if (pSelected[i]) {
                          TParticle* iparticle = (TParticle *) particles->At(i);
-                         Int_t kf  = iparticle->GetPdgCode();
-                         Int_t ipa = iparticle->GetFirstMother()-1;
+                         Int_t kf   = iparticle->GetPdgCode();
+                         Int_t ksc  = iparticle->GetStatusCode();
+                         Int_t jpa  = iparticle->GetFirstMother()-1;
                          
                          och[0] = origin0[0]+iparticle->Vx()/10;
                          och[1] = origin0[1]+iparticle->Vy()/10;
@@ -431,17 +464,18 @@ void AliGenParam::Generate()
                          pc[1]  = iparticle->Py();
                          pc[2]  = iparticle->Pz();
                          
-                         if (ipa > -1) {
-                             iparent = pParent[ipa];
+                         if (jpa > -1) {
+                             iparent = pParent[jpa];
                          } else {
                              iparent = -1;
                          }
                         
-                         PushTrack(fTrackIt*trackIt[i], iparent, kf,
+                         PushTrack(fTrackIt * trackIt[i], iparent, kf,
                                           pc, och, polar,
-                                          0, kPDecay, nt, wgtch);
+                                          0, kPDecay, nt, wgtch, ksc);
                          pParent[i] = nt;
                          KeepTrack(nt); 
+                         fNprimaries++;
                      } // Selected
                  } // Particle loop 
              }  // Decays by Lujet
@@ -454,13 +488,20 @@ void AliGenParam::Generate()
          else  // nodecay option, so parent will be tracked by GEANT (pions, kaons, eta, omegas, baryons)
          {
            gAlice->GetMCApp()->
-               PushTrack(fTrackIt,-1,iPart,p,origin0,polar,0,kPPrimary,nt,wgtp);
+               PushTrack(fTrackIt,-1,iPart,p,origin0,polar,0,kPPrimary,nt,wgtp, 1);
             ipa++; 
+           fNprimaries++;
          }
          break;
     } // while
   } // event loop
+  
   SetHighWaterMark(nt);
+
+  AliGenEventHeader* header = new AliGenEventHeader("PARAM");
+  header->SetPrimaryVertex(fVertex);
+  header->SetNProduced(fNprimaries);
+  AddHeader(header);
 }
 //____________________________________________________________________________________
 Float_t AliGenParam::GetRelativeArea(Float_t ptMin, Float_t ptMax, Float_t yMin, Float_t yMax, Float_t phiMin, Float_t phiMax)
@@ -469,8 +510,8 @@ Float_t AliGenParam::GetRelativeArea(Float_t ptMin, Float_t ptMax, Float_t yMin,
 // Normalisation for selected kinematic region
 //
   Float_t ratio =  
-    fPtPara->Integral(ptMin,ptMax) / fPtPara->Integral( fPtPara->GetXmin(), fPtPara->GetXmax()) *
-    fYPara->Integral(yMin,yMax)/fYPara->Integral(fYPara->GetXmin(),fYPara->GetXmax())   *
+    fPtPara->Integral(ptMin,ptMax,(Double_t *)0,1.e-6) / fPtPara->Integral( fPtPara->GetXmin(), fPtPara->GetXmax(),(Double_t *)0,1.e-6) *
+    fYPara->Integral(yMin,yMax,(Double_t *)0,1.e-6)/fYPara->Integral(fYPara->GetXmin(),fYPara->GetXmax(),(Double_t *)0,1.e-6)   *
     (phiMax-phiMin)/360.;
   return TMath::Abs(ratio);
 }
@@ -492,12 +533,6 @@ void AliGenParam::Draw( const char * /*opt*/)
      fYPara->GetHistogram()->SetXTitle("y");     
 }
 
-AliGenParam& AliGenParam::operator=(const  AliGenParam& rhs)
-{
-// Assignment operator
-    rhs.Copy(*this);
-    return *this;
-}