Alifatal removed in GetWord. Put AliWarning (+AddMajorErrorLog)

[u/mrichter/AliRoot.git] / PYTHIA6 / AliGenPythia.cxx

diff --git a/PYTHIA6/AliGenPythia.cxx b/PYTHIA6/AliGenPythia.cxx
index d8238ce621d391878dd4dd68a1dead3cc2983b2f..175c7d89ac8718e76c7e1db5c1a37d46acd3a1ec 100644 (file)
--- a/PYTHIA6/AliGenPythia.cxx
+++ b/PYTHIA6/AliGenPythia.cxx
@@ -51,7 +51,6 @@ AliGenPythia::AliGenPythia():
     AliGenMC(),
     fProcess(kPyCharm),          
     fStrucFunc(kCTEQ5L), 
-    fEnergyCMS(5500.),
     fKineBias(0.),
     fTrials(0),
     fTrialsRun(0),
@@ -108,6 +107,8 @@ AliGenPythia::AliGenPythia():
     fHFoff(kFALSE),    
     fTriggerParticle(0),
     fTriggerEta(0.9),     
+    fTriggerMultiplicity(0),
+    fTriggerMultiplicityEta(0),
     fCountMode(kCountAll),      
     fHeader(0),  
     fRL(0),      
@@ -129,6 +130,7 @@ AliGenPythia::AliGenPythia():
 
 {
 // Default Constructor
+  fEnergyCMS = 5500.;
   SetNuclei(0,0);
   if (!AliPythiaRndm::GetPythiaRandom()) 
       AliPythiaRndm::SetPythiaRandom(GetRandom());
@@ -138,7 +140,6 @@ AliGenPythia::AliGenPythia(Int_t npart)
     :AliGenMC(npart),
      fProcess(kPyCharm),          
      fStrucFunc(kCTEQ5L), 
-     fEnergyCMS(5500.),
      fKineBias(0.),
      fTrials(0),
      fTrialsRun(0),
@@ -195,6 +196,8 @@ AliGenPythia::AliGenPythia(Int_t npart)
      fHFoff(kFALSE),    
      fTriggerParticle(0),
      fTriggerEta(0.9),     
+     fTriggerMultiplicity(0),
+     fTriggerMultiplicityEta(0),
      fCountMode(kCountAll),      
      fHeader(0),  
      fRL(0),      
@@ -218,13 +221,13 @@ AliGenPythia::AliGenPythia(Int_t npart)
 // semimuonic decay
 // structure function GRVHO
 //
+    fEnergyCMS = 5500.;
     fName = "Pythia";
     fTitle= "Particle Generator using PYTHIA";
     SetForceDecay();
     // Set random number generator 
     if (!AliPythiaRndm::GetPythiaRandom()) 
       AliPythiaRndm::SetPythiaRandom(GetRandom());
-    fParticles = new TClonesArray("TParticle",1000);
     SetNuclei(0,0);
  }
 
@@ -573,7 +576,7 @@ void AliGenPythia::Generate()
 //     printf("Calling hadronisation %d\n", fPythia->GetN());
        fPythia->Pyexec();      
        fTrials++;
-       fPythia->ImportParticles(fParticles,"All");
+       fPythia->ImportParticles(&fParticles,"All");
        Boost();
 //
 //
@@ -581,7 +584,7 @@ void AliGenPythia::Generate()
        Int_t i;
        
        fNprimaries = 0;
-       Int_t np = fParticles->GetEntriesFast();
+       Int_t np = fParticles.GetEntriesFast();
        
        if (np == 0) continue;
 //
@@ -611,7 +614,7 @@ void AliGenPythia::Generate()
            fProcess != kPyBeautyppMNRwmi) {
            
            for (i = 0; i < np; i++) {
-               TParticle* iparticle = (TParticle *) fParticles->At(i);
+               TParticle* iparticle = (TParticle *) fParticles.At(i);
                Int_t ks = iparticle->GetStatusCode();
                kf = CheckPDGCode(iparticle->GetPdgCode());
 // No initial state partons
@@ -638,20 +641,20 @@ void AliGenPythia::Generate()
                if (kf >= fFlavorSelect && kf <= 6) {
                    Int_t idau = iparticle->GetFirstDaughter() - 1;
                    if (idau > -1) {
-                       TParticle* daughter = (TParticle *) fParticles->At(idau);
+                       TParticle* daughter = (TParticle *) fParticles.At(idau);
                        Int_t pdgD = daughter->GetPdgCode();
                        if (pdgD == 91 || pdgD == 92) {
                            Int_t jmin = daughter->GetFirstDaughter() - 1;
                            Int_t jmax = daughter->GetLastDaughter()  - 1;                          
-                           for (Int_t j = jmin; j <= jmax; j++)
-                               ((TParticle *) fParticles->At(j))->SetFirstMother(i+1);
+                           for (Int_t jp = jmin; jp <= jmax; jp++)
+                               ((TParticle *) fParticles.At(jp))->SetFirstMother(i+1);
                        } // is string or cluster
                    } // has daughter
                } // heavy quark
                
 
                if (ipa > -1) {
-                   TParticle *  mother = (TParticle *) fParticles->At(ipa);
+                   TParticle *  mother = (TParticle *) fParticles.At(ipa);
                    kfMo = TMath::Abs(mother->GetPdgCode());
                }
                
@@ -745,7 +748,7 @@ void AliGenPythia::Generate()
            if (nc > 0) {
                for (i = 0; i<np; i++) {
                    if (!pSelected[i]) continue;
-                   TParticle *  iparticle = (TParticle *) fParticles->At(i);
+                   TParticle *  iparticle = (TParticle *) fParticles.At(i);
                    kf = CheckPDGCode(iparticle->GetPdgCode());
                    Int_t ks = iparticle->GetStatusCode();  
                    p[0] = iparticle->Px();
@@ -831,15 +834,15 @@ Int_t  AliGenPythia::GenerateMB()
     
 
     
-    Int_t np = (fHadronisation) ? fParticles->GetEntriesFast() : fNpartons;
+    Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
 
 
 
     Int_t* pParent = new Int_t[np];
     for (i=0; i< np; i++) pParent[i] = -1;
     if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
-       TParticle* jet1 = (TParticle *) fParticles->At(6);
-       TParticle* jet2 = (TParticle *) fParticles->At(7);
+       TParticle* jet1 = (TParticle *) fParticles.At(6);
+       TParticle* jet2 = (TParticle *) fParticles.At(7);
        if (!CheckTrigger(jet1, jet2)) {
          delete [] pParent;
          return 0;
@@ -856,11 +859,11 @@ Int_t  AliGenPythia::GenerateMB()
       else if (fPi0InCalo) pdg = 111 ; // Pi0
 
       for (i=0; i< np; i++) {
-       TParticle* iparticle = (TParticle *) fParticles->At(i);
+       TParticle* iparticle = (TParticle *) fParticles.At(i);
        if(iparticle->GetStatusCode()==1 && iparticle->GetPdgCode()==pdg && 
           iparticle->Pt() > fFragPhotonOrPi0MinPt){
          Int_t imother = iparticle->GetFirstMother() - 1;
-         TParticle* pmother = (TParticle *) fParticles->At(imother);
+         TParticle* pmother = (TParticle *) fParticles.At(imother);
          if(pdg == 111 || 
             (pdg == 22 && pmother->GetStatusCode() != 11))//No photon from hadron decay
            {
@@ -876,6 +879,40 @@ Int_t  AliGenPythia::GenerateMB()
        return 0;
     }
     
+    // Check for minimum multiplicity
+    if (fTriggerMultiplicity > 0) {
+      Int_t multiplicity = 0;
+      for (i = 0; i < np; i++) {
+       TParticle *  iparticle = (TParticle *) fParticles.At(i);
+       
+       Int_t statusCode = iparticle->GetStatusCode();
+       
+       // Initial state particle
+       if (statusCode > 20)
+         continue;
+       
+       // skip quarks and gluons
+       Int_t pdgCode = TMath::Abs(iparticle->GetPdgCode());
+       if (pdgCode <= 10 || pdgCode == 21)
+         continue;
+       
+       if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
+         continue;
+       
+       TParticlePDG* pdgPart = iparticle->GetPDG();
+       if (pdgPart && pdgPart->Charge() == 0)
+         continue;
+       
+       ++multiplicity;
+      }
+
+      if (multiplicity < fTriggerMultiplicity) {
+       delete [] pParent;
+       return 0;
+      }
+      
+      Printf("Triggered on event with multiplicity of %d > %d", multiplicity, fTriggerMultiplicity);
+    }    
     
      // Select events with a photon  pt > min pt going to PHOS eta acceptance or exactly PHOS eta phi
     if ((fProcess == kPyJets || fProcess == kPyDirectGamma) && fPhotonInCalo && (fCheckPHOSeta || fCheckPHOS)){
@@ -885,14 +922,14 @@ Int_t  AliGenPythia::GenerateMB()
       Int_t pdg  = 22; 
       Int_t iphcand = -1;
       for (i=0; i< np; i++) {
-        TParticle* iparticle = (TParticle *) fParticles->At(i);
+        TParticle* iparticle = (TParticle *) fParticles.At(i);
         Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
         Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax 
         
         if(iparticle->GetStatusCode() == 1 
            && iparticle->GetPdgCode() == pdg   
            && iparticle->Pt() > fPhotonMinPt    
-           && eta > fPHOSEta){                 
+           && eta < fPHOSEta){                 
            
            // first check if the photon is in PHOS phi
            if(IsInPHOS(phi,eta)){ 
@@ -914,7 +951,7 @@ Int_t  AliGenPythia::GenerateMB()
     if (fTriggerParticle) {
        Bool_t triggered = kFALSE;
        for (i = 0; i < np; i++) {
-           TParticle *  iparticle = (TParticle *) fParticles->At(i);
+           TParticle *  iparticle = (TParticle *) fParticles.At(i);
            kf = CheckPDGCode(iparticle->GetPdgCode());
            if (kf != fTriggerParticle) continue;
            if (iparticle->Pt() == 0.) continue;
@@ -937,7 +974,7 @@ Int_t  AliGenPythia::GenerateMB()
       Float_t yQ;  
       Int_t   pdgQ;
       for(i=0; i<np; i++) {
-       hvq = (TParticle*)fParticles->At(i);
+       hvq = (TParticle*)fParticles.At(i);
        pdgQ = hvq->GetPdgCode();  
        if(TMath::Abs(pdgQ) != fFlavorSelect) continue; 
        if(pdgQ>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
@@ -967,7 +1004,7 @@ Int_t  AliGenPythia::GenerateMB()
 
     for (i = 0; i < np; i++) {
        Int_t trackIt = 0;
-       TParticle *  iparticle = (TParticle *) fParticles->At(i);
+       TParticle *  iparticle = (TParticle *) fParticles.At(i);
        kf = CheckPDGCode(iparticle->GetPdgCode());
        Int_t ks = iparticle->GetStatusCode();
        Int_t km = iparticle->GetFirstMother();
@@ -1153,6 +1190,7 @@ void AliGenPythia::MakeHeader()
 //  Pass header
 //
     AddHeader(fHeader);
+    fHeader = 0x0;
 }
 
 Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
@@ -1214,10 +1252,10 @@ Bool_t AliGenPythia::CheckKinematicsOnChild(){
     Int_t nPartAcc = 0; //number of particles in the acceptance range
     Int_t numberOfAcceptedParticles = 1;
     if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
-    Int_t npart = fParticles->GetEntriesFast();
+    Int_t npart = fParticles.GetEntriesFast();
     
     for (j = 0; j<npart; j++) {
-       TParticle *  jparticle = (TParticle *) fParticles->At(j);
+       TParticle *  jparticle = (TParticle *) fParticles.At(j);
        kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
        ks = jparticle->GetStatusCode();
        km = jparticle->GetFirstMother(); 
@@ -1418,19 +1456,19 @@ void AliGenPythia::RotatePhi(Int_t iphcand, Bool_t& okdd)
   Double_t phiPHOS = gRandom->Uniform(phiPHOSmin,phiPHOSmax);
   
   //calculate deltaphi
-  TParticle* ph = (TParticle *) fParticles->At(iphcand);
+  TParticle* ph = (TParticle *) fParticles.At(iphcand);
   Double_t phphi = ph->Phi();
   Double_t deltaphi = phiPHOS - phphi;
 
   
   
   //loop for all particles and produce the phi rotation
-  Int_t np = (fHadronisation) ? fParticles->GetEntriesFast() : fNpartons;
+  Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
   Double_t oldphi, newphi;
   Double_t newVx, newVy, R, Vz, time; 
   Double_t newPx, newPy, pt, Pz, e;
   for(Int_t i=0; i< np; i++) {
-      TParticle* iparticle = (TParticle *) fParticles->At(i);
+      TParticle* iparticle = (TParticle *) fParticles.At(i);
       oldphi = iparticle->Phi();
       newphi = oldphi + deltaphi;
       if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle