Setter for calling Generate n-times
[u/mrichter/AliRoot.git] / EVGEN / AliGenHIJINGpara.cxx
index 9d5cb51..737b0cd 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-/*
-$Log$
-Revision 1.4  2000/11/30 07:12:50  alibrary
-Introducing new Rndm and QA classes
-
-Revision 1.3  2000/10/02 21:28:06  fca
-Removal of useless dependecies via forward declarations
-
-Revision 1.2  2000/07/11 18:24:55  fca
-Coding convention corrections + few minor bug fixes
-
-Revision 1.1  2000/06/09 20:20:30  morsch
-Same class as previously in AliSimpleGen.cxx
-All coding rule violations except RS3 corrected (AM)
-
-*/
+/* $Id$ */
 
 // Parameterisation of pi and K, eta and pt distributions
 // used for the ALICE TDRs.
@@ -52,44 +37,54 @@ All coding rule violations except RS3 corrected (AM)
 //                                                               //
 ///////////////////////////////////////////////////////////////////
 
+#include <TArrayF.h>
+#include <TCanvas.h>
+#include <TClonesArray.h>
+#include <TDatabasePDG.h>
+#include <TF1.h>
+#include <TH1.h>
+#include <TPDGCode.h>
+#include <TParticle.h>
+#include <TROOT.h>
+#include <TVirtualMC.h>
+
+#include "AliConst.h"
+#include "AliDecayer.h"
+#include "AliGenEventHeader.h"
 #include "AliGenHIJINGpara.h"
-#include "TF1.h"
+#include "AliLog.h"
 #include "AliRun.h"
-#include "AliConst.h"
-#include "AliPDG.h"
 
 ClassImp(AliGenHIJINGpara)
 
-AliGenHIJINGpara::AliGenHIJINGpara(const AliGenHIJINGpara & para)
-{
-// copy constructor
-}
+
 
 //_____________________________________________________________________________
-static Double_t ptpi(Double_t *px, Double_t *)
+static Double_t ptpi(const Double_t *px, const Double_t *)
 {
   //
   //     PT-PARAMETERIZATION CDF, PRL 61(88) 1819
   //     POWER LAW FOR PT > 500 MEV
   //     MT SCALING BELOW (T=160 MEV)
   //
-  const Double_t kp0 = 1.3;
-  const Double_t kxn = 8.28;
-  const Double_t kxlim=0.5;
-  const Double_t kt=0.160;
-  const Double_t kxmpi=0.139;
-  const Double_t kb=1.;
+  const Double_t kp0    = 1.3;
+  const Double_t kxn    = 8.28;
+  const Double_t kxlim  = 0.5;
+  const Double_t kt     = 0.160;
+  const Double_t kxmpi  = 0.139;
+  const Double_t kb     = 1.;
   Double_t y, y1, xmpi2, ynorm, a;
-  Double_t x=*px;
+  Double_t x = *px;
   //
-  y1=TMath::Power(kp0/(kp0+kxlim),kxn);
-  xmpi2=kxmpi*kxmpi;
-  ynorm=kb*(TMath::Exp(-sqrt(kxlim*kxlim+xmpi2)/kt));
-  a=ynorm/y1;
+  y1 = TMath::Power(kp0 / (kp0 + kxlim), kxn);
+  xmpi2 = kxmpi * kxmpi;
+  ynorm = kb * (TMath::Exp(-sqrt(kxlim * kxlim + xmpi2) / kt ));
+  a = ynorm / y1;
   if (x > kxlim)
-    y=a*TMath::Power(kp0/(kp0+x),kxn);
+    y = a * TMath::Power(kp0 / (kp0 + x), kxn);
   else
-    y=kb*TMath::Exp(-sqrt(x*x+xmpi2)/kt);
+    y = kb* TMath::Exp(-sqrt(x * x + xmpi2) / kt);
+  
   return y*x;
 }
 
@@ -163,30 +158,46 @@ static Double_t etakac( Double_t *py, Double_t *)
 
 //_____________________________________________________________________________
 AliGenHIJINGpara::AliGenHIJINGpara()
-  :AliGenerator()
+    :AliGenerator(),
+       fNt(-1),
+       fNpartProd(0),
+       fPi0Decays(kFALSE),
+       fPtWgtPi(0.),
+       fPtWgtKa(0.),
+       fPtpi(0),
+       fPtka(0),
+       fETApic(0),
+       fETAkac(0),
+       fDecayer(0)
 {
     //
     // Default constructor
     //
-    fPtpi = 0;
-    fPtka = 0;
-    fETApic = 0;
-    fETAkac = 0;
+    SetCutVertexZ();
+    SetPtRange();
 }
 
 //_____________________________________________________________________________
 AliGenHIJINGpara::AliGenHIJINGpara(Int_t npart)
-  :AliGenerator(npart)
+    :AliGenerator(npart),
+       fNt(-1),
+       fNpartProd(npart),
+       fPi0Decays(kFALSE),
+       fPtWgtPi(0.),
+       fPtWgtKa(0.),
+       fPtpi(0),
+       fPtka(0),
+       fETApic(0),
+       fETAkac(0),
+       fDecayer(0)
 {
   // 
   // Standard constructor
   //
-    fName="HIGINGpara";
+    fName="HIJINGpara";
     fTitle="HIJING Parametrisation Particle Generator";
-    fPtpi = 0;
-    fPtka = 0;
-    fETApic = 0;
-    fETAkac = 0;
+    SetCutVertexZ();
+    SetPtRange();
 }
 
 //_____________________________________________________________________________
@@ -211,33 +222,75 @@ void AliGenHIJINGpara::Init()
        TMath::Min((Double_t)fThetaMax/2,TMath::Pi()/2-1.e-10)));
     Float_t etaMax = -TMath::Log(TMath::Tan(
        TMath::Max((Double_t)fThetaMin/2,1.e-10)));
-    fPtpi = new TF1("ptpi",&ptpi,0,20,0);
-    fPtka = new TF1("ptka",&ptka,0,20,0);
+    fPtpi   = new TF1("ptpi",&ptpi,0,20,0);
+    gROOT->GetListOfFunctions()->Remove(fPtpi);
+    fPtka   = new TF1("ptka",&ptka,0,20,0);
+    gROOT->GetListOfFunctions()->Remove(fPtka);
+    fPtpi->SetNpx(1000);
+    fPtka->SetNpx(1000);
     fETApic = new TF1("etapic",&etapic,etaMin,etaMax,0);
+    gROOT->GetListOfFunctions()->Remove(fETApic);
     fETAkac = new TF1("etakac",&etakac,etaMin,etaMax,0);
-    TF1 *etaPic0 = new TF1("etapic",&etapic,-7,7,0);
-    TF1 *etaKac0 = new TF1("etakac",&etakac,-7,7,0);
-    Float_t intETApi  = etaPic0->Integral(-0.5, 0.5);
-    Float_t intETAka  = etaKac0->Integral(-0.5, 0.5);
-    Float_t scalePi=7316/(intETApi/1.5);
-    Float_t scaleKa= 684/(intETAka/2.0);
+    gROOT->GetListOfFunctions()->Remove(fETAkac);
+
+    TF1 etaPic0("etaPic0",&etapic,-7,7,0);
+    TF1 etaKac0("etaKac0",&etakac,-7,7,0);
+
+    TF1 ptPic0("ptPic0",&ptpi,0.,15.,0);
+    TF1 ptKac0("ptKac0",&ptka,0.,15.,0);
+
+    Float_t intETApi  = etaPic0.Integral(-0.5, 0.5);
+    Float_t intETAka  = etaKac0.Integral(-0.5, 0.5);
+    Float_t scalePi   = 7316/(intETApi/1.5);
+    Float_t scaleKa   =  684/(intETAka/2.0);
+
+//  Fraction of events corresponding to the selected pt-range    
+    Float_t intPt    = (0.877*ptPic0.Integral(0, 15)+
+                       0.123*ptKac0.Integral(0, 15));
+    Float_t intPtSel = (0.877*ptPic0.Integral(fPtMin, fPtMax)+
+                       0.123*ptKac0.Integral(fPtMin, fPtMax));
+    Float_t ptFrac   = intPtSel/intPt;
+
+//  Fraction of events corresponding to the selected eta-range    
+    Float_t intETASel  = (scalePi*etaPic0.Integral(etaMin, etaMax)+
+                         scaleKa*etaKac0.Integral(etaMin, etaMax));
+//  Fraction of events corresponding to the selected phi-range    
+    Float_t phiFrac    = (fPhiMax-fPhiMin)/2/TMath::Pi();
+
     
-    Float_t intPt  = (0.877*etaPic0->Integral(0, 15)+
-                     0.123*etaKac0->Integral(0, 15));
-    Float_t intPtSel = (0.877*etaPic0->Integral(fPtMin, fPtMax)+
-                       0.123*etaKac0->Integral(fPtMin, fPtMax));
-    Float_t ptFrac = intPtSel/intPt;
+    fParentWeight = (intETASel*ptFrac*phiFrac) / Float_t(fNpart);
     
+    if (fAnalog != 0) {
+       fPtWgtPi = (fPtMax - fPtMin) / fPtpi->Integral(0., 20.);
+       fPtWgtKa = (fPtMax - fPtMin) / fPtka->Integral(0., 20.);
+       fParentWeight = (intETASel*phiFrac) / Float_t(fNpart);
+    }
     
-    Float_t intETASel  = (scalePi*etaPic0->Integral(etaMin, etaMax)+
-                         scaleKa*etaKac0->Integral(etaMin, etaMax));
-    Float_t phiFrac = (fPhiMax-fPhiMin)/2/TMath::Pi();
-    fParentWeight = Float_t(fNpart)/intETASel*ptFrac*phiFrac;
     
-    printf("\n The number of particles in the selected kinematic region corresponds to %f percent of a full event\n ", 100.*fParentWeight);
+    AliInfo(Form("The number of particles in the selected kinematic region corresponds to %f percent of a full event", 
+                100./ fParentWeight));
+
+// Issue warning message if etaMin or etaMax are outside the alowed range 
+// of the parametrization
+    if (etaMin < -8.001 || etaMax > 8.001) {
+       AliWarning("\nYOU ARE USING THE PARAMETERISATION OUTSIDE ");    
+       AliWarning("THE ALLOWED PSEUDORAPIDITY RANGE (-8. - 8.)");          
+       AliWarning(Form("YOUR LIMITS: %f %f \n ", etaMin, etaMax));
+    }
+//
+//
+    if (fPi0Decays && gMC)
+       fDecayer = gMC->GetDecayer();
+
+    if (fPi0Decays)
+    {
+       fDecayer->SetForceDecay(kNeutralPion);
+       fDecayer->Init();
+    }
     
 }
 
+
 //_____________________________________________________________________________
 void AliGenHIJINGpara::Generate()
 {
@@ -254,10 +307,11 @@ void AliGenHIJINGpara::Generate()
     const Int_t kKaons[4] = {kK0Long, kK0Short, kKPlus, kKMinus};
     //
     Float_t origin[3];
-    Float_t pt, pl, ptot;
+    Float_t time;
+    Float_t pt, pl, ptot, wgt;
     Float_t phi, theta;
     Float_t p[3];
-    Int_t i, part, nt, j;
+    Int_t i, part, j;
     //
     TF1 *ptf;
     TF1 *etaf;
@@ -265,29 +319,45 @@ void AliGenHIJINGpara::Generate()
     Float_t random[6];
     //
     for (j=0;j<3;j++) origin[j]=fOrigin[j];
-    if(fVertexSmear==kPerEvent) {
-       Rndm(random,6);
-       for (j=0;j<3;j++) {
-           origin[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
-               TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
-       }
-    }
+    time = fTimeOrigin;
+
+    if(fVertexSmear == kPerEvent) {
+       Vertex();
+       for (j=0; j < 3; j++) origin[j] = fVertex[j];
+       time = fTime;
+    } // if kPerEvent
+    TArrayF eventVertex;
+    eventVertex.Set(3);
+    eventVertex[0] = origin[0];
+    eventVertex[1] = origin[1];
+    eventVertex[2] = origin[2];
+    Float_t eventTime = time;
+     
     for(i=0;i<fNpart;i++) {
        while(1) {
-           Rndm(random,3);
+           Rndm(random,4);
            if(random[0]<kBorne) {
                part=kPions[Int_t (random[1]*3)];
                ptf=fPtpi;
-             etaf=fETApic;
+               etaf=fETApic;
+               wgt = fPtWgtPi;
            } else {
                part=kKaons[Int_t (random[1]*4)];
                ptf=fPtka;
                etaf=fETAkac;
+               wgt = fPtWgtKa;
            }
            phi=fPhiMin+random[2]*(fPhiMax-fPhiMin);
            theta=2*TMath::ATan(TMath::Exp(-etaf->GetRandom()));
            if(theta<fThetaMin || theta>fThetaMax) continue;
-           pt=ptf->GetRandom();
+        
+           if (fAnalog == 0) { 
+               pt = ptf->GetRandom();
+           } else {
+               pt = fPtMin + random[3] * (fPtMax - fPtMin);
+           }
+           
+           
            pl=pt/TMath::Tan(theta);
            ptot=TMath::Sqrt(pt*pt+pl*pl);
            if(ptot<fPMin || ptot>fPMax) continue;
@@ -300,17 +370,100 @@ void AliGenHIJINGpara::Generate()
                    origin[j]=fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
                        TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
                }
+               Rndm(random,2);
+               time = fTimeOrigin + fOsigma[2]/TMath::Ccgs()*
+                 TMath::Cos(2*random[0]*TMath::Pi())*
+                 TMath::Sqrt(-2*TMath::Log(random[1]));
+           }
+
+           if (fAnalog == 0) { 
+               wgt = fParentWeight;
+           } else {
+               wgt *= (fParentWeight * ptf->Eval(pt));
+           }
+           
+           
+           if (part == kPi0 && fPi0Decays){
+//
+//          Decay pi0 if requested
+               PushTrack(0,-1,part,p,origin,polar,time,kPPrimary,fNt,wgt);
+               KeepTrack(fNt);
+               DecayPi0(origin, p, time);
+           } else {
+      // printf("fNt %d", fNt);
+               PushTrack(fTrackIt,-1,part,p,origin,polar,time,kPPrimary,fNt,wgt);
+
+               KeepTrack(fNt);
            }
-           gAlice->SetTrack(fTrackIt,-1,part,p,origin,polar,0,kPPrimary,nt,fParentWeight);
+
            break;
        }
+       SetHighWaterMark(fNt);
     }
+//
+
+// Header
+    AliGenEventHeader* header = new AliGenEventHeader("HIJINGparam");
+// Event Vertex
+    header->SetPrimaryVertex(eventVertex);
+    header->SetInteractionTime(eventTime);
+    header->SetNProduced(fNpartProd);
+    gAlice->SetGenEventHeader(header); 
+}
+
+void AliGenHIJINGpara::SetPtRange(Float_t ptmin, Float_t ptmax) {
+    AliGenerator::SetPtRange(ptmin, ptmax);
 }
 
-AliGenHIJINGpara& AliGenHIJINGpara::operator=(const  AliGenHIJINGpara& rhs)
+void AliGenHIJINGpara::DecayPi0(Float_t* orig, Float_t * p, Float_t time) 
 {
-// Assignment operator
-    return *this;
+//
+//    Decay the pi0
+//    and put decay products on the stack
+//
+    static TClonesArray *particles;
+    if(!particles) particles = new TClonesArray("TParticle",1000);
+//    
+    const Float_t kMass = TDatabasePDG::Instance()->GetParticle(kPi0)->Mass();
+    Float_t       e     = TMath::Sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2]+ kMass * kMass);
+//
+//  Decay the pi0    
+    TLorentzVector pmom(p[0], p[1], p[2], e);
+    fDecayer->Decay(kPi0, &pmom);
+    
+//
+// Put decay particles on the stack
+//
+    Float_t polar[3] = {0., 0., 0.};
+    Int_t np = fDecayer->ImportParticles(particles);
+    fNpartProd += (np-1);
+    Int_t nt = 0;    
+    for (Int_t i = 1; i < np; i++)
+    {
+       TParticle* iParticle =  (TParticle *) particles->At(i);
+       p[0] = iParticle->Px();
+       p[1] = iParticle->Py();
+       p[2] = iParticle->Pz();
+       Int_t part = iParticle->GetPdgCode();
+
+       PushTrack(fTrackIt, fNt, part, p, orig, polar, time, kPDecay, nt, fParentWeight);
+       KeepTrack(nt);
+    }
+    fNt = nt;
 }
 
+void AliGenHIJINGpara::Draw( const char * /*opt*/)
+{
+    //
+    // Draw the pT and y Distributions
+    //
+     TCanvas *c0 = new TCanvas("c0","Canvas 0",400,10,600,700);
+     c0->Divide(2,1);
+     c0->cd(1);
+     fPtpi->Draw();
+     fPtpi->GetHistogram()->SetXTitle("p_{T} (GeV)");     
+     c0->cd(2);
+     fPtka->Draw();
+     fPtka->GetHistogram()->SetXTitle("p_{T} (GeV)");     
 
+}