]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EVGEN/AliGenParam.cxx
New class added.
[u/mrichter/AliRoot.git] / EVGEN / AliGenParam.cxx
index 0fad3ae8c77f26765801fa8f1d676d37e98ca937..71847f8a7cfb9636f7639ea1d31b3cd7a4c89e90 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,76 @@ 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),
+       fSelectAll(kFALSE),
+       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, const AliGenLib * Library,  Int_t param, const 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),
+     fSelectAll(kFALSE),
+     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),
+    fSelectAll(kFALSE),
+    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 +130,35 @@ 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),
+     fSelectAll(kFALSE),
+     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 +168,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);
 }
 
 //____________________________________________________________
@@ -166,22 +190,25 @@ void AliGenParam::Init()
   */
   //End_Html
     char name[256];
-    sprintf(name, "pt-parameterisation for %s", GetName());
+    snprintf(name, 256, "pt-parameterisation for %s", GetName());
     
     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);
     
     fPtPara->SetNpx(npx);
 
-    sprintf(name, "y-parameterisation  for %s", GetName());
+    snprintf(name, 256, "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());
+    snprintf(name, 256, "pt-for-%s", GetName());
     TF1 ptPara(name ,fPtParaFunc, 0, 15, 0);
-    sprintf(name, "y-for-%s", GetName());
+    snprintf(name, 256, "y-for-%s", GetName());
     TF1 yPara(name, fYParaFunc, -6, 6, 0);
 
 //
@@ -192,9 +219,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;
@@ -235,6 +262,7 @@ void AliGenParam::Generate()
 //
   Float_t polar[3]= {0,0,0};  // Polarisation of the parent particle (for GEANT tracking)
   Float_t origin0[3];         // Origin of the generated parent particle (for GEANT tracking)
+  Float_t time0;              // Time0 of the generated parent particle
   Float_t pt, pl, ptot;       // Transverse, logitudinal and total momenta of the parent particle
   Float_t phi, theta;         // Phi and theta spherical angles of the parent particle momentum
   Float_t p[3], pc[3], 
@@ -253,14 +281,18 @@ void AliGenParam::Generate()
  
 // Calculating vertex position per event
   for (j=0;j<3;j++) origin0[j]=fOrigin[j];
+  time0 = fTimeOrigin;
   if(fVertexSmear==kPerEvent) {
       Vertex();
       for (j=0;j<3;j++) origin0[j]=fVertex[j];
+      time0 = fTime;
   }
   
   Int_t ipa=0;
   
 // Generating fNpart particles
+  fNprimaries = 0;
+  
   while (ipa<fNpart) {
       while(1) {
 //
@@ -296,7 +328,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;
@@ -314,6 +346,10 @@ void AliGenParam::Generate()
                      fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
                      TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
              }
+             Rndm(random,2);
+             time0 = fTimeOrigin + fOsigma[2]/TMath::Ccgs()*
+               TMath::Cos(2*random[0]*TMath::Pi())*
+               TMath::Sqrt(-2*TMath::Log(random[1]));
          }
          
 // Looking at fForceDecay : 
@@ -323,6 +359,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,7 +387,8 @@ void AliGenParam::Generate()
              }
              
              if (np >1) {
-                 TParticle* iparticle =  (TParticle *) particles->At(0);
+                 decayed = kTRUE;
+                 TParticle* iparticle =  0;
                  Int_t ipF, ipL;
                  for (i = 1; i<np ; i++) {
                      trackIt[i] = 1;
@@ -385,7 +424,7 @@ void AliGenParam::Generate()
 //
 // children
                      
-                     if (ChildSelected(TMath::Abs(kf)) || fForceDecay == kAll && trackIt[i])
+                     if ((ChildSelected(TMath::Abs(kf)) || fForceDecay == kAll || fSelectAll) && trackIt[i])
                      {
                          if (fCutOnChild) {
                              pc[0]=iparticle->Px();
@@ -412,17 +451,22 @@ void AliGenParam::Generate()
                  ipa++;
 //
 // Parent
-                 PushTrack(0, -1, iPart, p, origin0, polar, 0, kPPrimary, nt, wgtp);
+                 
+                 
+                 PushTrack(0, -1, iPart, p, origin0, polar, time0, 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 +475,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);
+                                          time0 + iparticle->T(), kPDecay, nt, wgtch, ksc);
                          pParent[i] = nt;
                          KeepTrack(nt); 
+                         fNprimaries++;
                      } // Selected
                  } // Particle loop 
              }  // Decays by Lujet
@@ -454,13 +499,21 @@ 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,time0,kPPrimary,nt,wgtp, 1);
             ipa++; 
+           fNprimaries++;
          }
          break;
     } // while
   } // event loop
+  
   SetHighWaterMark(nt);
+
+  AliGenEventHeader* header = new AliGenEventHeader("PARAM");
+  header->SetPrimaryVertex(fVertex);
+  header->SetInteractionTime(fTime);
+  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 +522,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 +545,6 @@ void AliGenParam::Draw( const char * /*opt*/)
      fYPara->GetHistogram()->SetXTitle("y");     
 }
 
-AliGenParam& AliGenParam::operator=(const  AliGenParam& rhs)
-{
-// Assignment operator
-    rhs.Copy(*this);
-    return *this;
-}