]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TDPMjet/AliGenDPMjet.cxx
Coverity warning 23946
[u/mrichter/AliRoot.git] / TDPMjet / AliGenDPMjet.cxx
index 7e3900e31fe9aa168b81e04e010ff7d89fa90db2..25a4a0fcf82c5970cfaf5bb2a432c170ad939131 100644 (file)
@@ -19,6 +19,7 @@
 // Uses the TDPMjet implementation of TGenerator.
 
 #include <TDPMjet.h>
+#include "DPMcommon.h"
 #include <TRandom.h>
 #include <TArrayI.h>
 #include <TParticle.h>
@@ -34,6 +35,7 @@
 #include "AliGenDPMjetEventHeader.h"
 #include "AliRun.h"
 #include "AliDpmJetRndm.h"
+#include "AliIonPDGCodes.h"
 #include "AliHeader.h"
 #include "AliStack.h"
 #include "AliMC.h"
@@ -44,13 +46,14 @@ ClassImp(AliGenDPMjet)
 //______________________________________________________________________________
 AliGenDPMjet::AliGenDPMjet()
     :AliGenMC(), 
-     fBeamEn(2750.),
+     fBeamEn(0.),
      fMinImpactParam(0.),
-     fMaxImpactParam(5.),
+     fMaxImpactParam(100.),
      fICentr(0),
      fSelectAll(0),
      fFlavor(0),
      fTrials(0),
+     fNprimaries(0),
      fSpectators(1),
      fSpecn(0),
      fSpecp(0),
@@ -61,11 +64,18 @@ AliGenDPMjet::AliGenDPMjet()
      fDecayAll(0),
      fGenImpPar(0.),
      fProcess(kDpmMb),
+     fTriggerParticle(0),
+     fTriggerEta(0.9),     
+     fTriggerMinPt(-1),  
+     fTriggerMaxPt(1000),  
      fTriggerMultiplicity(0),
      fTriggerMultiplicityEta(0),
      fTriggerMultiplicityPtMin(0),
      fkTuneForDiff(0),
-     fProcDiff(0)
+     fProcDiff(0),
+     fFragmentation(kFALSE),
+     fHeader(AliGenDPMjetEventHeader("DPMJET"))
+
 {
 // Constructor
     fEnergyCMS = 5500.;
@@ -76,13 +86,14 @@ AliGenDPMjet::AliGenDPMjet()
 //______________________________________________________________________________
 AliGenDPMjet::AliGenDPMjet(Int_t npart)
     :AliGenMC(npart),
-     fBeamEn(2750.),
+     fBeamEn(0.),
      fMinImpactParam(0.),
-     fMaxImpactParam(5.),
+     fMaxImpactParam(100.),
      fICentr(0),
      fSelectAll(0),
      fFlavor(0),
      fTrials(0),
+     fNprimaries(0),
      fSpectators(1),
      fSpecn(0),
      fSpecp(0),
@@ -93,11 +104,19 @@ AliGenDPMjet::AliGenDPMjet(Int_t npart)
      fDecayAll(0),
      fGenImpPar(0.),
      fProcess(kDpmMb),
+     fTriggerParticle(0),
+     fTriggerEta(0.9),     
+     fTriggerMinPt(-1),  
+     fTriggerMaxPt(1000),  
      fTriggerMultiplicity(0),
      fTriggerMultiplicityEta(0),
      fTriggerMultiplicityPtMin(0),
      fkTuneForDiff(0),
-     fProcDiff(0)
+     fProcDiff(0),
+     fFragmentation(kFALSE),
+     fHeader(AliGenDPMjetEventHeader("DPMJET"))
+
+
 {
 // Default PbPb collisions at 5. 5 TeV
 //
@@ -112,13 +131,14 @@ AliGenDPMjet::AliGenDPMjet(Int_t npart)
 
 AliGenDPMjet::AliGenDPMjet(const AliGenDPMjet &/*Dpmjet*/)
     :AliGenMC(),
-     fBeamEn(2750.),
+     fBeamEn(0.),
      fMinImpactParam(0.),
-     fMaxImpactParam(5.),
+     fMaxImpactParam(100.),
      fICentr(0),
      fSelectAll(0),
      fFlavor(0),
      fTrials(0),
+     fNprimaries(0),
      fSpectators(1),
      fSpecn(0),
      fSpecp(0),
@@ -129,11 +149,19 @@ AliGenDPMjet::AliGenDPMjet(const AliGenDPMjet &/*Dpmjet*/)
      fDecayAll(0),
      fGenImpPar(0.),
      fProcess(kDpmMb),
+     fTriggerParticle(0),
+     fTriggerEta(0.9),     
+     fTriggerMinPt(-1),  
+     fTriggerMaxPt(1000),  
      fTriggerMultiplicity(0),
      fTriggerMultiplicityEta(0),
      fTriggerMultiplicityPtMin(0),
      fkTuneForDiff(0),
-     fProcDiff(0)
+     fProcDiff(0),
+     fFragmentation(kFALSE),
+     fHeader(0x0)
+
+
 {
     // Dummy copy constructor
     fEnergyCMS = 5500.;
@@ -149,25 +177,25 @@ void AliGenDPMjet::Init()
 {
 // Initialization
     
+    if(fEnergyCMS>0. && fBeamEn<0.1) fBeamEn = fEnergyCMS/2;
     SetMC(new TDPMjet(fProcess, fAProjectile, fZProjectile, fATarget, fZTarget, 
                      fBeamEn,fEnergyCMS));
 
     fDPMjet=(TDPMjet*) fMCEvGen;
-    //
+    if (fAProjectile == 1 && TMath::Abs(fZProjectile == 1)) fDPMjet->SetfIdp(1);
+    
     // **** Flag to force central production
     // fICentr=1. central production forced 
     // fICentr<0 && fICentr>-100 -> bmin = fMinImpactParam, bmax = fMaxImpactParam       
     // fICentr<-99 -> fraction of x-sec. = XSFRAC                
     // fICentr=-1. -> evaporation/fzc suppressed                 
-    // fICentr<-1. -> evaporation/fzc allowed            
-    if (fAProjectile == 1 && TMath::Abs(fZProjectile == 1)) fDPMjet->SetfIdp(1);
-    
+    // fICentr<-1. -> evaporation/fzc allowed                
     fDPMjet->SetfFCentr(fICentr);  
+    
     fDPMjet->SetbRange(fMinImpactParam, fMaxImpactParam); 
     fDPMjet->SetPi0Decay(fPi0Decay);
     fDPMjet->SetDecayAll(fDecayAll);
-
-    AliGenMC::Init();
+    fDPMjet->SetFragmentProd(fFragmentation);
     
 //
 //  Initialize DPMjet  
@@ -193,6 +221,7 @@ void AliGenDPMjet::Generate()
   Int_t kf, ks, imo;
   kf = 0;
   fTrials = 0;
+  fNprimaries = 0;
   //  Set collision vertex position 
   if (fVertexSmear == kPerEvent) Vertex();
   
@@ -246,6 +275,21 @@ void AliGenDPMjet::Generate()
        Printf("Triggered on event with multiplicity of %d >= %d", multiplicity, fTriggerMultiplicity);
       }    
 
+  //Trigger on the presence of a given particle in some phase space
+    if (fTriggerParticle) {
+       Bool_t triggered = kFALSE;
+           for (Long_t i = 0; i < np; i++) {
+               TParticle *  iparticle = (TParticle *) fParticles.At(i);
+               kf = CheckPDGCode(iparticle->GetPdgCode());
+               if (kf != fTriggerParticle) continue;
+               if (iparticle->Pt() == 0.) continue;
+               if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
+               if ( iparticle->Pt() > fTriggerMaxPt || iparticle->Pt() < fTriggerMinPt ) continue;
+               triggered = kTRUE;
+               break;
+           }
+      if (!triggered) continue; 
+    }
 
       if(fkTuneForDiff && ( (TMath::Abs(fEnergyCMS - 900) < 1) || (TMath::Abs(fEnergyCMS - 2760) < 1) || (TMath::Abs(fEnergyCMS - 7000) < 1)) ) {
        if(!CheckDiffraction() ) continue;
@@ -362,13 +406,16 @@ void AliGenDPMjet::Generate()
 
 
              
-             Bool_t tFlag = (fTrackIt && (ks == 1));
+             Bool_t tFlag = (fTrackIt && (ks==1 || ks==-1));
+             //printf(" AliGemDPMJet->PushTrack: kf %d  ks %d  flag %d\n",kf,ks,tFlag);
+             if(kf>10000 && (ks==-1 || ks==1000 || ks==1001)) kf += 1000000000;
              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);
+             fNprimaries++;
              newPos[i] = nt;
          } // if selected
       } // particle loop
@@ -445,9 +492,8 @@ Bool_t AliGenDPMjet::Stable(TParticle*  particle)
 {
 // Return true for a stable particle
 //
-    
-//    if (particle->GetFirstDaughter() < 0 ) return kTRUE;
-    if (particle->GetStatusCode() == 1) return kTRUE;
+    int st = particle->GetStatusCode();
+    if(st == 1 || st == -1) return kTRUE;
     else return kFALSE;
 
 }
@@ -455,39 +501,40 @@ 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->GetProjParticipants(), 
-                                                        fDPMjet->GetTargParticipants());
-
-    if(fProcDiff>0){
-    ((AliGenDPMjetEventHeader*) header)->SetProcessType(fProcDiff);
-    }
-    else 
-      ((AliGenDPMjetEventHeader*) header)->SetProcessType(fDPMjet->GetProcessCode());
+    fHeader.SetNProduced(fNprimaries);
+    fHeader.SetImpactParameter(fDPMjet->GetBImpac());
+    fHeader.SetTotalEnergy(fDPMjet->GetTotEnergy());
+    fHeader.SetParticipants(fDPMjet->GetProjParticipants(), 
+                           fDPMjet->GetTargParticipants());
+
+    fHeader.SetCollisions(DTGLCP.ncp, DTGLCP.nct, 
+       fDPMjet->GetProjWounded(),fDPMjet->GetTargWounded());
+               
+    if(fProcDiff>0)  fHeader.SetProcessType(fProcDiff);
+    else fHeader.SetProcessType(fDPMjet->GetProcessCode());
 
     // Bookkeeping for kinematic bias
-    ((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
+    fHeader.SetTrials(fTrials);
     // Event Vertex
-    header->SetPrimaryVertex(fVertex);
-    header->SetInteractionTime(fTime);
-    gAlice->SetGenEventHeader(header);    
-    AddHeader(header);
+    fHeader.SetPrimaryVertex(fVertex);
+    fHeader.SetInteractionTime(fTime);
+    fHeader.SetNDiffractive(POEVT1.nsd1, POEVT1.nsd2, POEVT1.ndd);
+//    gAlice->SetGenEventHeader(fHeader);    
+    AddHeader(&fHeader);
+    fCollisionGeometry = &fHeader;
 }
 
-void AliGenDPMjet::AddHeader(AliGenEventHeader* header)
+//______________________________________________________________________________
+/*void AliGenDPMjet::AddHeader(AliGenEventHeader* fHeader)
 {
-    // Add header to container or runloader
+    // Add fHeader to container or runloader
     if (fContainer) {
-        fContainer->AddHeader(header);
+        fContainer->AddHeader(fHeader);
     } else {
-        AliRunLoader::Instance()->GetHeader()->SetGenEventHeader(header);
+        AliRunLoader::Instance()->GetHeader()->SetGenEventHeader(fHeader);
     }
-}
+}*/
 
 
 //______________________________________________________________________________
@@ -616,7 +663,9 @@ Bool_t AliGenDPMjet::CheckDiffraction()
        TParticle *  part = (TParticle *) fParticles.At(iPart);
        Double_t E= part->Energy();
        Double_t P= part->P();
-       M= TMath::Sqrt((fEnergyCMS-E-P)*(fEnergyCMS-E+P));
+       Double_t M2 = (fEnergyCMS-E-P)*(fEnergyCMS-E+P);
+       if(M2<0)  return kFALSE;
+       M= TMath::Sqrt(M2);
      }
 
      Double_t Mmin, Mmax, wSD, wDD, wND;
@@ -664,8 +713,15 @@ Bool_t AliGenDPMjet::CheckDiffraction()
     return kTRUE;
 }
 
+// -------------------------------------------------------
+void AliGenDPMjet::SetIonPDGCodes()
+{
+   // Defining PDG codes for the ions
+   AliIonPDGCodes *pdgcodes = new AliIonPDGCodes();
+   pdgcodes->AddParticlesToPdgDataBase();
+}
 
-
+// -------------------------------------------------------
 Bool_t AliGenDPMjet::GetWeightsDiffraction(Double_t M, Double_t &Mmin, Double_t &Mmax, 
                                                        Double_t &wSD, Double_t &wDD, Double_t &wND)
 {