]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TDPMjet/AliGenDPMjet.cxx
double for 4-momentum vector in Generate()
[u/mrichter/AliRoot.git] / TDPMjet / AliGenDPMjet.cxx
index 4692ccf07907d226d265136157c061822d725050..c3d15078b162032fdeab3dbda1aa73d3c2e3cdf7 100644 (file)
@@ -44,7 +44,6 @@ ClassImp(AliGenDPMjet)
 AliGenDPMjet::AliGenDPMjet()
     :AliGenMC(), 
      fBeamEn(2750.),
-     fEnergyCMS(5500.),
      fMinImpactParam(0.),
      fMaxImpactParam(5.),
      fICentr(0),
@@ -57,11 +56,16 @@ AliGenDPMjet::AliGenDPMjet()
      fDPMjet(0),
      fNoGammas(0),
      fLHC(0),
-     fPi0Decay(0),
+     fPi0Decay(1),
+     fDecayAll(0),
      fGenImpPar(0.),
-     fProcess(kDpmMb)
+     fProcess(kDpmMb),
+     fTriggerMultiplicity(0),
+     fTriggerMultiplicityEta(0),
+     fTriggerMultiplicityPtMin(0)
 {
 // Constructor
+    fEnergyCMS = 5500.;
     AliDpmJetRndm::SetDpmJetRandom(GetRandom());
 }
 
@@ -70,7 +74,6 @@ AliGenDPMjet::AliGenDPMjet()
 AliGenDPMjet::AliGenDPMjet(Int_t npart)
     :AliGenMC(npart),
      fBeamEn(2750.),
-     fEnergyCMS(5500.),
      fMinImpactParam(0.),
      fMaxImpactParam(5.),
      fICentr(0),
@@ -83,12 +86,17 @@ AliGenDPMjet::AliGenDPMjet(Int_t npart)
      fDPMjet(0),
      fNoGammas(0),
      fLHC(0),
-     fPi0Decay(0),
+     fPi0Decay(1),
+     fDecayAll(0),
      fGenImpPar(0.),
-     fProcess(kDpmMb)
+     fProcess(kDpmMb),
+     fTriggerMultiplicity(0),
+     fTriggerMultiplicityEta(0),
+     fTriggerMultiplicityPtMin(0)
 {
 // Default PbPb collisions at 5. 5 TeV
 //
+    fEnergyCMS = 5500.;
     fName = "DPMJET";
     fTitle= "Particle Generator using DPMJET";
     SetTarget();
@@ -100,7 +108,6 @@ AliGenDPMjet::AliGenDPMjet(Int_t npart)
 AliGenDPMjet::AliGenDPMjet(const AliGenDPMjet &/*Dpmjet*/)
     :AliGenMC(),
      fBeamEn(2750.),
-     fEnergyCMS(5500.),
      fMinImpactParam(0.),
      fMaxImpactParam(5.),
      fICentr(0),
@@ -113,11 +120,16 @@ AliGenDPMjet::AliGenDPMjet(const AliGenDPMjet &/*Dpmjet*/)
      fDPMjet(0),
      fNoGammas(0),
      fLHC(0),
-     fPi0Decay(0),
+     fPi0Decay(1),
+     fDecayAll(0),
      fGenImpPar(0.),
-     fProcess(kDpmMb)
+     fProcess(kDpmMb),
+     fTriggerMultiplicity(0),
+     fTriggerMultiplicityEta(0),
+     fTriggerMultiplicityPtMin(0)
 {
     // Dummy copy constructor
+    fEnergyCMS = 5500.;
 }
 
 //______________________________________________________________________________
@@ -141,11 +153,12 @@ void AliGenDPMjet::Init()
     // fICentr<-99 -> fraction of x-sec. = XSFRAC                
     // fICentr=-1. -> evaporation/fzc suppressed                 
     // fICentr<-1. -> evaporation/fzc suppressed                 
-    if (fAProjectile == 1 && fZProjectile == 1) fDPMjet->SetfIdp(1);
+    if (fAProjectile == 1 && TMath::Abs(fZProjectile == 1)) fDPMjet->SetfIdp(1);
     
     fDPMjet->SetfFCentr(fICentr);  
     fDPMjet->SetbRange(fMinImpactParam, fMaxImpactParam); 
     fDPMjet->SetPi0Decay(fPi0Decay);
+    fDPMjet->SetDecayAll(fDecayAll);
 //
 //  Initialize DPMjet  
 //    
@@ -158,9 +171,9 @@ void AliGenDPMjet::Generate()
 {
 // Generate one event
 
-  Float_t polar[3]    =   {0,0,0};
-  Float_t origin[3]   =   {0,0,0};
-  Float_t p[3];
+  Double_t polar[3]    =   {0,0,0};
+  Double_t origin[3]   =   {0,0,0};
+  Double_t p[4]        =   {0};
   Float_t tof;
 
 //  converts from mm/c to s
@@ -181,6 +194,7 @@ void AliGenDPMjet::Generate()
       fSpecp = 0;
 // --------------------------------------------------------------------------
       fDPMjet->GenerateEvent();
+      
       fTrials++;
 
       fDPMjet->ImportParticles(&fParticles,"All");      
@@ -188,11 +202,43 @@ void AliGenDPMjet::Generate()
 
       // Temporaneo
       fGenImpPar = fDPMjet->GetBImpac();
-      
+
+      if(TMath::Abs(fXingAngleY) > 1.e-10) BeamCrossAngle();
+
       Int_t np = fParticles.GetEntriesFast();
-      printf("\n **************************************************%d\n",np);
-      Int_t nc=0;
-      if (np==0) continue;
+      //
+      // Multiplicity Trigger
+      if (fTriggerMultiplicity > 0) {
+       Int_t multiplicity = 0;
+       for (Int_t i = 0; i < np; i++) {
+         TParticle *  iparticle = (TParticle *) fParticles.At(i);
+       
+         Int_t statusCode = iparticle->GetStatusCode();
+       
+         // Initial state particle
+         if (statusCode != 1)
+           continue;
+         // eta cut
+         if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
+           continue;
+         // pt cut
+         if (iparticle->Pt() < fTriggerMultiplicityPtMin) 
+           continue;
+         
+         TParticlePDG* pdgPart = iparticle->GetPDG();
+         if (pdgPart && pdgPart->Charge() == 0)
+           continue;
+         ++multiplicity;
+       }
+       //
+       //
+       if (multiplicity < fTriggerMultiplicity) continue;
+       Printf("Triggered on event with multiplicity of %d >= %d", multiplicity, fTriggerMultiplicity);
+      }    
+
+      Int_t nc = 0;
+      if (np == 0) continue;
+
       Int_t i;
       Int_t* newPos     = new Int_t[np];
       Int_t* pSelected  = new Int_t[np];
@@ -215,15 +261,18 @@ void AliGenDPMjet::Generate()
          Bool_t  hasSelectedDaughters =  kFALSE;
          
          kf = iparticle->GetPdgCode();
-         if (kf == 92) continue;
+         if (kf == 92 || kf == 99999) continue;
          ks = iparticle->GetStatusCode();
 // No initial state partons
           if (ks==21) continue;
            
          if (!fSelectAll) selected = KinematicSelection(iparticle, 0) && 
                               SelectFlavor(kf);
+
+         
          hasSelectedDaughters = DaughtersSelection(iparticle);
 
+
 // Put particle on the stack if it is either selected or 
 // it is the mother of at least one seleted particle
 
@@ -281,6 +330,7 @@ void AliGenDPMjet::Generate()
              p[0] = iparticle->Px();
              p[1] = iparticle->Py();
              p[2] = iparticle->Pz();
+             p[3] = iparticle->Energy();
              origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
              origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
              origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
@@ -292,12 +342,17 @@ void AliGenDPMjet::Generate()
              if (hasMother) {
                  imo = iparticle->GetFirstMother();
                  mother = (TParticle *) fParticles.At(imo);
-                 imo = (mother->GetPdgCode() != 92) ? newPos[imo] : -1;
+                 imo = (mother->GetPdgCode() != 92 && mother->GetPdgCode() != 99999) ? newPos[imo] : -1;
              } // if has mother   
 
-             Bool_t tFlag = (fTrackIt && (ks == 1));
+
              
-             PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
+             Bool_t tFlag = (fTrackIt && (ks == 1));
+             PushTrack(tFlag, imo, kf, 
+                       p[0], p[1], p[2], p[3], 
+                       origin[0], origin[1], origin[2], tof,
+                       polar[0], polar[1], polar[2],
+                       kPNoProcess, nt, 1., ks);
              KeepTrack(nt);
              newPos[i] = nt;
          } // if selected
@@ -385,17 +440,18 @@ Bool_t AliGenDPMjet::Stable(TParticle*  particle)
 //______________________________________________________________________________
 void AliGenDPMjet::MakeHeader()
 {
+  printf("MakeHeader %13.3f \n", fDPMjet->GetBImpac());
 // Builds the event header, to be called after each event
     AliGenEventHeader* header = new AliGenDPMjetEventHeader("DPMJET");
     ((AliGenDPMjetEventHeader*) header)->SetNProduced(fDPMjet->GetNumStablePc());
     ((AliGenDPMjetEventHeader*) header)->SetImpactParameter(fDPMjet->GetBImpac());
     ((AliGenDPMjetEventHeader*) header)->SetTotalEnergy(fDPMjet->GetTotEnergy());
-    ((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetfIp(), 
-                                                        fDPMjet->GetfIt());
- ((AliGenDPMjetEventHeader*) header)->SetProcessType(fDPMjet->GetProcessCode());
-// Bookkeeping for kinematic bias
+    ((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetProjParticipants(), 
+                                                        fDPMjet->GetTargParticipants());
    ((AliGenDPMjetEventHeader*) header)->SetProcessType(fDPMjet->GetProcessCode());
+    // Bookkeeping for kinematic bias
     ((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
-// Event Vertex
+    // Event Vertex
     header->SetPrimaryVertex(fVertex);
     gAlice->SetGenEventHeader(header);    
     AddHeader(header);
@@ -407,7 +463,7 @@ void AliGenDPMjet::AddHeader(AliGenEventHeader* header)
     if (fContainer) {
         fContainer->AddHeader(header);
     } else {
-        AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(header);
+        AliRunLoader::Instance()->GetHeader()->SetGenEventHeader(header);
     }
 }