]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TDPMjet/AliGenDPMjet.cxx
Fix Coverity defect
[u/mrichter/AliRoot.git] / TDPMjet / AliGenDPMjet.cxx
index d528ab28ce0beddb9e2e1f3d0796473041f93e74..68ddec32abc0abbfd012df9871f1afba7d47e947 100644 (file)
 #include <TParticleClassPDG.h>
 #include <TPDGCode.h>
 #include <TLorentzVector.h>
-
+#include <TClonesArray.h>
+#include "AliRunLoader.h"
 #include "AliGenDPMjet.h"
 #include "AliGenDPMjetEventHeader.h"
-#include "AliPythia.h"
 #include "AliRun.h"
 #include "AliDpmJetRndm.h"
+#include "AliHeader.h"
+#include "AliStack.h"
+#include "AliMC.h"
 
 ClassImp(AliGenDPMjet)
 
@@ -41,7 +44,6 @@ ClassImp(AliGenDPMjet)
 AliGenDPMjet::AliGenDPMjet()
     :AliGenMC(), 
      fBeamEn(2750.),
-     fEnergyCMS(5500.),
      fMinImpactParam(0.),
      fMaxImpactParam(5.),
      fICentr(0),
@@ -52,14 +54,18 @@ AliGenDPMjet::AliGenDPMjet()
      fSpecn(0),
      fSpecp(0),
      fDPMjet(0),
-     fParticles(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());
 }
 
@@ -68,7 +74,6 @@ AliGenDPMjet::AliGenDPMjet()
 AliGenDPMjet::AliGenDPMjet(Int_t npart)
     :AliGenMC(npart),
      fBeamEn(2750.),
-     fEnergyCMS(5500.),
      fMinImpactParam(0.),
      fMaxImpactParam(5.),
      fICentr(0),
@@ -79,29 +84,30 @@ AliGenDPMjet::AliGenDPMjet(Int_t npart)
      fSpecn(0),
      fSpecp(0),
      fDPMjet(0),
-     fParticles(new TClonesArray("TParticle",10000)),
      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();
     SetProjectile();
     fVertex.Set(3);
-    // Instance AliPythia
-    AliPythia::Instance(); 
     AliDpmJetRndm::SetDpmJetRandom(GetRandom());
 }
 
 AliGenDPMjet::AliGenDPMjet(const AliGenDPMjet &/*Dpmjet*/)
     :AliGenMC(),
      fBeamEn(2750.),
-     fEnergyCMS(5500.),
      fMinImpactParam(0.),
      fMaxImpactParam(5.),
      fICentr(0),
@@ -112,21 +118,24 @@ AliGenDPMjet::AliGenDPMjet(const AliGenDPMjet &/*Dpmjet*/)
      fSpecn(0),
      fSpecp(0),
      fDPMjet(0),
-     fParticles(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.;
 }
 
 //______________________________________________________________________________
 AliGenDPMjet::~AliGenDPMjet()
 {
 // Destructor
-    delete fParticles;
 }
 //______________________________________________________________________________
 void AliGenDPMjet::Init()
@@ -144,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  
 //    
@@ -163,7 +173,6 @@ void AliGenDPMjet::Generate()
 
   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];
   Float_t tof;
 
@@ -185,18 +194,51 @@ void AliGenDPMjet::Generate()
       fSpecp = 0;
 // --------------------------------------------------------------------------
       fDPMjet->GenerateEvent();
+      
       fTrials++;
 
-      fDPMjet->ImportParticles(fParticles,"All");      
+      fDPMjet->ImportParticles(&fParticles,"All");      
       if (fLHC) Boost();
 
       // Temporaneo
       fGenImpPar = fDPMjet->GetBImpac();
-      
-      Int_t np = fParticles->GetEntriesFast();
-      printf("\n **************************************************%d\n",np);
-      Int_t nc=0;
-      if (np==0) continue;
+
+      if(TMath::Abs(fXingAngleY) > 1.e-10) BeamCrossAngle();
+
+      Int_t np = fParticles.GetEntriesFast();
+      //
+      // 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];
@@ -206,16 +248,10 @@ void AliGenDPMjet::Generate()
          pSelected[i] = 0;
       }
       
-//      Get event vertex
-      
-      fVertex[0] = origin0[0];
-      fVertex[1] = origin0[1]; 
-      fVertex[2] = origin0[2];
-      
 //      First select parent particles
 
       for (i = 0; i<np; i++) {
-         TParticle *iparticle = (TParticle *) fParticles->At(i);
+         TParticle *iparticle = (TParticle *) fParticles.At(i);
 
 // Is this a parent particle ?
 
@@ -225,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
 
@@ -247,7 +286,7 @@ void AliGenDPMjet::Generate()
 
 
       for (i=0; i<np; i++) {
-         TParticle *iparticle = (TParticle *) fParticles->At(i);
+         TParticle *iparticle = (TParticle *) fParticles.At(i);
 
 // Is this a final state particle ?
 
@@ -281,7 +320,7 @@ void AliGenDPMjet::Generate()
 // Write particles to stack
 
       for (i = 0; i<np; i++) {
-         TParticle *  iparticle = (TParticle *) fParticles->At(i);
+         TParticle *  iparticle = (TParticle *) fParticles.At(i);
          Bool_t  hasMother   = (iparticle->GetFirstMother()>=0);
          if (pSelected[i]) {
              
@@ -301,12 +340,13 @@ void AliGenDPMjet::Generate()
              TParticle* mother = 0;
              if (hasMother) {
                  imo = iparticle->GetFirstMother();
-                 mother = (TParticle *) fParticles->At(imo);
-                 imo = (mother->GetPdgCode() != 92) ? imo = newPos[imo] : -1;
+                 mother = (TParticle *) fParticles.At(imo);
+                 imo = (mother->GetPdgCode() != 92 && mother->GetPdgCode() != 99999) ? newPos[imo] : -1;
              } // if has mother   
 
-             Bool_t tFlag = (fTrackIt && (ks == 1));
+
              
+             Bool_t tFlag = (fTrackIt && (ks == 1));
              PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
              KeepTrack(nt);
              newPos[i] = nt;
@@ -341,7 +381,7 @@ Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle)
        imin = iparticle->GetFirstDaughter();
        imax = iparticle->GetLastDaughter();       
        for (i = imin; i <= imax; i++){
-           TParticle *  jparticle = (TParticle *) fParticles->At(i);   
+           TParticle *  jparticle = (TParticle *) fParticles.At(i);    
            Int_t ip = jparticle->GetPdgCode();
            if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
                selected=kTRUE; break;
@@ -402,14 +442,24 @@ void AliGenDPMjet::MakeHeader()
     ((AliGenDPMjetEventHeader*) header)->SetTotalEnergy(fDPMjet->GetTotEnergy());
     ((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetfIp(), 
                                                         fDPMjet->GetfIt());
-
+ ((AliGenDPMjetEventHeader*) header)->SetProcessType(fDPMjet->GetProcessCode());
 // Bookkeeping for kinematic bias
     ((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
 // Event Vertex
     header->SetPrimaryVertex(fVertex);
     gAlice->SetGenEventHeader(header);    
+    AddHeader(header);
 }
 
+void AliGenDPMjet::AddHeader(AliGenEventHeader* header)
+{
+    // Add header to container or runloader
+    if (fContainer) {
+        fContainer->AddHeader(header);
+    } else {
+        AliRunLoader::Instance()->GetHeader()->SetGenEventHeader(header);
+    }
+}
 
 
 //______________________________________________________________________________
@@ -420,5 +470,12 @@ AliGenDPMjet& AliGenDPMjet::operator=(const  AliGenDPMjet& /*rhs*/)
 }
 
 
+void AliGenDPMjet::FinishRun()
+{
+    // Print run statistics
+    fDPMjet->Dt_Dtuout();
+}
+
+    
 
 //______________________________________________________________________________