Import All particles from Pythia and add CutOnChild
[u/mrichter/AliRoot.git] / EVGEN / AliGenParam.cxx
index 7e38529ff020fe65eb0396948a33d31b9cb1f096..12566847dd94520f3d95d63dff1f1b82dbebfa3e 100644 (file)
@@ -1,13 +1,12 @@
 #include "AliGenParam.h"
 #include "AliGenMUONlib.h"
 #include "AliRun.h"
-#include "TGeant3.h"
 #include "AliPythia.h"
 #include <TDirectory.h>
 #include <TFile.h>
 #include <TTree.h>
 #include <stdlib.h>
-#include <TMCParticle.h>
+#include <TParticle.h>
 
 ClassImp(AliGenParam)
 
@@ -15,7 +14,7 @@ ClassImp(AliGenParam)
 
   //Begin_Html
   /*
-    <img src="gif/AliGenParam.gif">
+    <img src="picts/AliGenParam.gif">
   */
   //End_Html
 
@@ -28,6 +27,7 @@ AliGenParam::AliGenParam()
   fYPara  = 0;
   fParam  = jpsi_p;
   fAnalog = analog;
+  SetCutOnChild();
 }
 
 //____________________________________________________________
@@ -50,6 +50,7 @@ AliGenParam::AliGenParam(Int_t npart, Param_t param)
   fChildSelect.Set(5);
   for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
   ForceDecay();
+  SetCutOnChild();
 }
 
 //____________________________________________________________
@@ -68,7 +69,7 @@ void AliGenParam::Init()
 //  End of the test !!!
   //Begin_Html
   /*
-    <img src="gif/AliGenParam.gif">
+    <img src="picts/AliGenParam.gif">
   */
   //End_Html
  
@@ -135,37 +136,41 @@ void AliGenParam::Generate()
 // the desired theta, phi and momentum windows; Gaussian smearing 
 // on the vertex is done if selected
 
-  AliMC* pMC = AliMC::GetMC();
+
+  //printf("Generate !!!!!!!!!!!!!\n");
 
   Float_t polar[3]= {0,0,0};
   //
   Float_t origin[3], origin0[3];
   Float_t pt, pl, ptot;
   Float_t phi, theta;
-  Float_t p[3];
+  Float_t p[3], pc[3], och[3], pch[10][3];
   Float_t ty, xmt;
-  Int_t i, nt, j;
+  Int_t nt, i, j, kfch[10];
   Float_t  wgtp, wgtch;
   Double_t dummy;
+  static TClonesArray *particles;
   //
+  if(!particles) particles=new TClonesArray("TParticle",1000);
   //
   Float_t random[6];
   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]));
       }
   }
-  for(i=0;i<fNpart;i++) {
-      while(1) {
+  Int_t ipa=0;
+  while (ipa<fNpart) {
+    while(1) {
 //
 // particle type
          Int_t Ipart = fIpParaFunc();
          fChildWeight=(fPythia->GetBraPart(Ipart))*fParentWeight;        
          Float_t am=fPythia->GetPMAS(fPythia->LuComp(Ipart),1);
-         pMC->Rndm(random,2);
+         gMC->Rndm(random,2);
 //
 // phi
          phi=fPhiMin+random[0]*(fPhiMax-fPhiMin);
@@ -194,53 +199,89 @@ void AliGenParam::Generate()
          p[1]=pt*TMath::Sin(phi);
          p[2]=pl;
          if(fVertexSmear==perTrack) {
-             pMC->Rndm(random,6);
+             gMC->Rndm(random,6);
              for (j=0;j<3;j++) {
                  origin0[j]=
                      fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
                      TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
              }
          }
-         
-         Int_t kg=fPythia->GetGeantCode(Ipart);
-//
-// parent
-
-         gAlice->
-             SetTrack(0,-1,kg,p,origin,polar,0,"Primary",nt,wgtp);
-         Int_t iparent=nt;
 //
 // use lujet to decay particle
 
          Float_t energy=TMath::Sqrt(ptot*ptot+am*am);
          fPythia->DecayParticle(Ipart,energy,theta,phi);
-//       fPythia->LuList(1);
-         
+         //      fPythia->LuList(1);
+
+
+         //printf("origin0 %f %f %f\n",origin0[0],origin0[1],origin0[2]);
+         //printf("fCutOnChild %d \n",fCutOnChild);
 //
 // select muons
-         TObjArray* particles = fPythia->GetPrimaries() ;
-         Int_t np = particles->GetEntriesFast();
-         for (Int_t i = 0; i<np; i++) {
-             TMCParticle *  iparticle = (TMCParticle *) particles->At(i);
-             Int_t kf = iparticle->GetKF();
+         Int_t np=fPythia->ImportParticles(particles,"All");
+          //printf("np     %d \n",np);
+         Int_t ncsel=0;
+         for (i = 1; i<np; i++) {
+             TParticle *  iparticle = (TParticle *) particles->At(i);
+             Int_t kf = iparticle->GetPdgCode();
+              //printf("kf %d\n",kf);
 //
 // children
              if (ChildSelected(TMath::Abs(kf)))
              {
-                 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;
-                 gAlice->SetTrack(fTrackIt,iparent,fPythia->GetGeantCode(kf),
-                                  p,origin,polar,
-                                  0,"Decay",nt,wgtch);
-                 gAlice->KeepTrack(nt);
+                 pc[0]=iparticle->Px();
+                 pc[1]=iparticle->Py();
+                 pc[2]=iparticle->Pz();
+                 och[0]=origin0[0]+iparticle->Vx()/10;
+                 och[1]=origin0[1]+iparticle->Vy()/10;
+                 och[2]=origin0[2]+iparticle->Vz()/10;
+                 if (fCutOnChild) {
+                   Float_t PtChild=TMath::Sqrt(pc[0]*pc[0]+pc[1]*pc[1]);
+                   Float_t PChild=TMath::Sqrt(PtChild*PtChild+pc[2]*pc[2]);
+                   Float_t ThetaChild=TMath::ATan2(PtChild,pc[2]);
+                   Float_t PhiChild=TMath::ATan2(pc[1],pc[0])+TMath::Pi();
+                   if ((PtChild   > fPtMin   && PtChild   <fPtMax)      &&
+                       (PChild    > fPMin    && PChild    <fPMax)       &&
+                       (ThetaChild>fThetaMin && ThetaChild<fThetaMax)   &&
+                       (PhiChild  >  fPhiMin && PhiChild  <fPhiMax))
+                     {
+                       pch[ncsel][0]=pc[0];
+                       pch[ncsel][1]=pc[1];
+                       pch[ncsel][2]=pc[2];
+                       kfch[ncsel]=kf;
+                       ncsel++;
+                     } else {
+                       ncsel=-1;
+                       break;
+                     } // child kine cuts
+                 } else {
+                   pch[ncsel][0]=pc[0];
+                   pch[ncsel][1]=pc[1];
+                   pch[ncsel][2]=pc[2];
+                   kfch[ncsel]=kf;
+                   ncsel++;
+                 } // if child selection
              } // select muon
          } // decay particle loop
+         Int_t iparent;
+         if ((fCutOnChild && ncsel >0) || !fCutOnChild){
+             ipa++;
+//
+// parent
+             gAlice->
+                 SetTrack(0,-1,Ipart,p,origin,polar,0,"Primary",nt,wgtp);
+             iparent=nt;
+
+             for (i=0; i< ncsel; i++) {
+                 gAlice->SetTrack(fTrackIt,iparent,kfch[i],
+                                  &pch[i][0],och,polar,
+                                  0,"Decay",nt,wgtch);
+                 gAlice->KeepTrack(nt); 
+             }
+             
+         } // kinematic selection
          break;
-      } // kinematic selection
+    } // while
   } // event loop
 }
 
@@ -253,11 +294,11 @@ Bool_t AliGenParam::ChildSelected(Int_t ip)
     return kFALSE;
 }
 
-Bool_t AliGenParam::KinematicSelection(TMCParticle *particle)
+Bool_t AliGenParam::KinematicSelection(TParticle *particle)
 {
-    Float_t px=particle->GetPx();
-    Float_t py=particle->GetPy();
-    Float_t pz=particle->GetPz();
+    Float_t px=particle->Px();
+    Float_t py=particle->Py();
+    Float_t pz=particle->Pz();
 //
 // momentum cut
     Float_t p=TMath::Sqrt(px*px+py*py+pz*pz);