Use ParticlesImport with TParticle instead than TMCParticle
[u/mrichter/AliRoot.git] / EVGEN / AliGenPythia.cxx
index 34f3ee7de7da74af6e5997fd2a7f387df2e9f046..2329c4dcc7eaf84338ce4866fddc27bec8f8f722 100644 (file)
@@ -1,6 +1,5 @@
 #include "AliGenerator.h"
 #include "AliGenPythia.h"
-#include "TGeant3.h"
 #include "AliRun.h"
 #include "AliPythia.h"
 #include <TDirectory.h>
@@ -8,8 +7,8 @@
 #include <TTree.h>
 #include <stdlib.h>
 #include <AliPythia.h>
-#include <TMCParticle.h>
-#include <GParticle.h>
+#include <TParticle.h>
+//#include <GParticle.h>
  ClassImp(AliGenPythia)
 
 AliGenPythia::AliGenPythia()
@@ -48,10 +47,10 @@ void AliGenPythia::Init()
 //
 //  Forward Paramters to the AliPythia object
     fPythia->DefineParticles();
-    fPythia->ForceDecay(fForceDecay);
     fPythia->SetCKIN(3,fPtHardMin);
     fPythia->SetCKIN(4,fPtHardMax);    
     fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
+    fPythia->ForceDecay(fForceDecay);
     fPythia->LuList(0);
     fPythia->PyStat(2);
 //  Parent and Children Selection
@@ -87,6 +86,8 @@ void AliGenPythia::Init()
     case jpsi:
        fParentSelect[0]=443;
        break;
+    case mb:
+       break;
     }
 
     switch (fForceDecay) 
@@ -101,6 +102,8 @@ void AliGenPythia::Init()
     case dimuon:
     case b_jpsi_dimuon:
     case b_psip_dimuon:
+    case pitomu:
+    case katomu:
        fChildSelect[0]=13;
        break;
     }
@@ -108,24 +111,29 @@ void AliGenPythia::Init()
 
 void AliGenPythia::Generate()
 {
-    AliMC* pMC = AliMC::GetMC();
 
-    Float_t polar[3] = {0,0,0};
-    Float_t origin[3]= {0,0,0};
-    Float_t origin0[3]= {0,0,0};
-    Float_t p[3], random[6];
+    Float_t polar[3] =   {0,0,0};
+    Float_t origin[3]=   {0,0,0};
+    Float_t origin_p[3]= {0,0,0};
+    Float_t origin0[3]=  {0,0,0};
+    Float_t p[3], p_p[4], random[6];
+    static TClonesArray *particles;
+//  converts from mm/c to s
+    const Float_t kconv=0.001/2.999792458e8;
     
     
 //
-    printf("Generate event");
     Int_t nt=0;
+    Int_t nt_p=0;
     Int_t jev=0;
     Int_t j;
+
+    if(!particles) particles=new TClonesArray("TParticle",1000);
     
     fTrials=0;
     for (j=0;j<3;j++) origin0[j]=fOrigin[j];
     if(fVertexSmear==perEvent) {
-       pMC->Rndm(random,6);
+       gMC->Rndm(random,6);
        for (j=0;j<3;j++) {
            origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
                TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
@@ -137,74 +145,113 @@ void AliGenPythia::Generate()
            fPythia->SetPARP(151+j, fOsigma[j]*10.);
        }
     }
-    
-
     while(1)
     {
        fPythia->PyEvnt();
        fTrials++;
-       TObjArray* particles = fPythia->GetPrimaries() ;
+       fPythia->ImportParticles(particles);
        Int_t np = particles->GetEntriesFast();
        printf("\n **************************************************%d\n",np);
        Int_t nc=0;
        if (np == 0 ) continue;
-       for (Int_t i = 0; i<np; i++) {
-           TMCParticle *  iparticle = (TMCParticle *) particles->At(i);
-           Int_t kf = iparticle->GetKF();
-           fChildWeight=(fPythia->GetBraPart(kf))*fParentWeight;         
+       if (fProcess != mb) {
+           for (Int_t i = 0; i<np; i++) {
+               TParticle *  iparticle = (TParticle *) particles->At(i);
+               Int_t kf = iparticle->GetPdgCode();
+               fChildWeight=(fPythia->GetBraPart(kf))*fParentWeight;     
 //
 // Parent
-           if (ParentSelected(TMath::Abs(kf))) {
-               
-               if (KinematicSelection(iparticle)) {
-                   nc++;
+               if (ParentSelected(TMath::Abs(kf))) {
+                   if (KinematicSelection(iparticle)) {
+                       if (nc==0) {
+//
+// Store information concerning the hard scattering process
+//
+                           Float_t mass_p = fPythia->GetPARI(13);
+                           Float_t   pt_p = fPythia->GetPARI(17);
+                           Float_t    y_p = fPythia->GetPARI(37);
+                           Float_t  xmt_p = sqrt(pt_p*pt_p+mass_p*mass_p);
+                           Float_t     ty = Float_t(TMath::TanH(y_p));
+                           p_p[0] = pt_p;
+                           p_p[1] = 0;
+                           p_p[2] = xmt_p*ty/sqrt(1.-ty*ty);
+                           p_p[3] = mass_p;
+                           gAlice->SetTrack(0,-1,-1,
+                                            p_p,origin_p,polar,
+                                            0,"Hard Scat.",nt_p,fParentWeight);
+                           gAlice->KeepTrack(nt_p);
+                       }
+                       nc++;
 //
 // store parent track information
-                   p[0]=iparticle->GetPx();
-                   p[1]=iparticle->GetPy();
-                   p[2]=iparticle->GetPz();
-                   origin[0]=origin0[0]+iparticle->GetVx()/10;
-                   origin[1]=origin0[1]+iparticle->GetVy()/10;
-                   origin[2]=origin0[2]+iparticle->GetVz()/10;
+                       p[0]=iparticle->Px();
+                       p[1]=iparticle->Py();
+                       p[2]=iparticle->Pz();
+                       origin[0]=origin0[0]+iparticle->Vx()/10;
+                       origin[1]=origin0[1]+iparticle->Vy()/10;
+                       origin[2]=origin0[2]+iparticle->Vz()/10;
 
-                   Int_t ifch=iparticle->GetFirstChild();
-                   Int_t ilch=iparticle->GetLastChild();       
-                   if (ifch !=0 && ilch !=0) {
-                       gAlice->SetTrack(0,-1,fPythia->GetGeantCode(kf),
-                                        p,origin,polar,
-                                        0,"Primary",nt,fParentWeight);
-                       Int_t iparent = nt;
+                       Int_t ifch=iparticle->GetFirstDaughter();
+                       Int_t ilch=iparticle->GetLastDaughter();        
+                       if (ifch !=0 && ilch !=0) {
+                           gAlice->SetTrack(0,nt_p,kf,
+                                            p,origin,polar,
+                                            0,"Primary",nt,fParentWeight);
+                           gAlice->KeepTrack(nt);
+                           Int_t iparent = nt;
 //
 // Children        
 
-                       for (Int_t j=ifch; j<=ilch; j++)
-                       {
-                           TMCParticle *  ichild = 
-                               (TMCParticle *) particles->At(j-1);
-                           Int_t kf = ichild->GetKF();
+                           for (Int_t j=ifch; j<=ilch; j++)
+                           {
+                               TParticle *  ichild = 
+                                   (TParticle *) particles->At(j-1);
+                               Int_t kf = ichild->GetPdgCode();
 //
 // 
-                           if (ChildSelected(TMath::Abs(kf))) {
-                               Int_t kg=fPythia->GetGeantCode(kf);
-                               origin[0]=ichild->GetVx();
-                               origin[1]=ichild->GetVy();
-                               origin[2]=ichild->GetVz();              
-                               p[0]=ichild->GetPx();
-                               p[1]=ichild->GetPy();
-                               p[2]=ichild->GetPz();
-                               Float_t tof=ichild->GetTime();
-                               gAlice->SetTrack(1, iparent, kg,
-                                                p,origin,polar,
-                                                tof,"Decay",nt,fChildWeight);
-                               gAlice->KeepTrack(nt);
-                           } // select child
-                       } // child loop
-                   }
-               } // kinematic selection
-           } // select particle 
-       } // particles
+                               if (ChildSelected(TMath::Abs(kf))) {
+                                   origin[0]=ichild->Vx();
+                                   origin[1]=ichild->Vy();
+                                   origin[2]=ichild->Vz();             
+                                   p[0]=ichild->Px();
+                                   p[1]=ichild->Py();
+                                   p[2]=ichild->Pz();
+                                   Float_t tof=kconv*ichild->T();
+                                   gAlice->SetTrack(fTrackIt, iparent, kf,
+                                                    p,origin,polar,
+                                                    tof,"Decay",nt,fChildWeight);
+                                   gAlice->KeepTrack(nt);
+                               } // select child
+                           } // child loop
+                       }
+                   } // kinematic selection
+               } // select particle
+           } // particle loop
+       } else {
+           for (Int_t i = 0; i<np; i++) {
+               TParticle *  iparticle = (TParticle *) particles->At(i);
+               Int_t kf = iparticle->GetPdgCode();
+               Int_t ks = iparticle->GetStatusCode();
+               if (ks==1 && kf!=0 && KinematicSelection(iparticle)) {
+                       nc++;
+//
+// store track information
+                       p[0]=iparticle->Px();
+                       p[1]=iparticle->Py();
+                       p[2]=iparticle->Pz();
+                       origin[0]=origin0[0]+iparticle->Vx()/10;
+                       origin[1]=origin0[1]+iparticle->Vy()/10;
+                       origin[2]=origin0[2]+iparticle->Vz()/10;
+                       Float_t tof=kconv*iparticle->T();
+                       gAlice->SetTrack(fTrackIt,-1,kf,p,origin,polar,
+                                        tof,"Primary",nt);
+                       gAlice->KeepTrack(nt);
+               } // select particle
+           } // particle loop 
+           printf("\n I've put %i particles on the stack \n",nc);
+       } // mb ?
        if (nc > 0) {
-           jev++;
+           jev+=nc;
            if (jev >= fNpart) {
                fKineBias=Float_t(fNpart)/Float_t(fTrials);
                printf("\n Trials: %i\n",fTrials);
@@ -236,12 +283,12 @@ Bool_t AliGenPythia::ChildSelected(Int_t ip)
     return kFALSE;
 }
 
-Bool_t AliGenPythia::KinematicSelection(TMCParticle *particle)
+Bool_t AliGenPythia::KinematicSelection(TParticle *particle)
 {
-    Float_t px=particle->GetPx();
-    Float_t py=particle->GetPy();
-    Float_t pz=particle->GetPz();
-    Float_t  e=particle->GetEnergy();
+    Float_t px=particle->Px();
+    Float_t py=particle->Py();
+    Float_t pz=particle->Pz();
+    Float_t  e=particle->Energy();
 
 //
 //  transverse momentum cut    
@@ -262,7 +309,7 @@ Bool_t AliGenPythia::KinematicSelection(TMCParticle *particle)
     
 //
 // theta cut
-    Float_t  theta = Float_t(TMath::ATan2(Double_t(pt),Double_t(p)));
+    Float_t  theta = Float_t(TMath::ATan2(Double_t(pt),Double_t(pz)));
     if (theta > fThetaMax || theta < fThetaMin) 
     {
 //     printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
@@ -292,11 +339,11 @@ Bool_t AliGenPythia::KinematicSelection(TMCParticle *particle)
 void AliGenPythia::AdjustWeights()
 {
     TClonesArray *PartArray = gAlice->Particles();
-    GParticle *Part;
+    TParticle *Part;
     Int_t ntrack=gAlice->GetNtrack();
     for (Int_t i=0; i<ntrack; i++) {
-       Part= (GParticle*) PartArray->UncheckedAt(i);
-       Part->SetWgt(Part->GetWgt()*fKineBias);
+       Part= (TParticle*) PartArray->UncheckedAt(i);
+       Part->SetWeight(Part->GetWeight()*fKineBias);
     }
 }
 
@@ -304,3 +351,7 @@ void AliGenPythia::AdjustWeights()
 
 
 
+
+
+
+