monitor dQdl (PID) at TRD entrance
[u/mrichter/AliRoot.git] / PYTHIA6 / AliGenPythiaPlus.cxx
index 35ba99a..4a19b5f 100644 (file)
@@ -52,7 +52,6 @@ AliGenPythiaPlus::AliGenPythiaPlus():
     fPythia(0),
     fProcess(kPyCharm),          
     fStrucFunc(kCTEQ5L), 
-    fEnergyCMS(5500.),
     fKineBias(0.),
     fTrials(0),
     fTrialsRun(0),
@@ -129,6 +128,7 @@ AliGenPythiaPlus::AliGenPythiaPlus():
 
 {
 // Default Constructor
+  fEnergyCMS = 5500.;
   SetNuclei(0,0);
   if (!AliPythiaRndm::GetPythiaRandom()) 
       AliPythiaRndm::SetPythiaRandom(GetRandom());
@@ -139,7 +139,6 @@ AliGenPythiaPlus::AliGenPythiaPlus(AliPythiaBase* pythia)
      fPythia(pythia),
      fProcess(kPyCharm),          
      fStrucFunc(kCTEQ5L), 
-     fEnergyCMS(5500.),
      fKineBias(0.),
      fTrials(0),
      fTrialsRun(0),
@@ -218,13 +217,13 @@ AliGenPythiaPlus::AliGenPythiaPlus(AliPythiaBase* pythia)
 // 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);
  }
 
@@ -376,6 +375,9 @@ void AliGenPythiaPlus::Init()
        fParentSelect[1] =   421;
        fParentSelect[2] =   431;
        fParentSelect[3] =  4122;
+       fParentSelect[4] =  4232;
+       fParentSelect[5] =  4132;
+       fParentSelect[6] =  4332;
        fFlavorSelect    =  4;  
        break;
     case kPyD0PbPbMNR:
@@ -396,7 +398,12 @@ void AliGenPythiaPlus::Init()
        fParentSelect[0] =   431;
        fFlavorSelect    =   4; 
        break;
+    case kPyLambdacppMNR:
+       fParentSelect[0] =  4122;
+       fFlavorSelect    =   4; 
+       break;      
     case kPyBeauty:
+    case kPyBeautyJets:
     case kPyBeautyPbPbMNR:
     case kPyBeautypPbMNR:
     case kPyBeautyppMNR:
@@ -424,8 +431,10 @@ void AliGenPythiaPlus::Init()
     case kPyJpsi:
        fParentSelect[0] = 443;
        break;
+    case kPyMbAtlasTuneMC09:
     case kPyMbDefault:
     case kPyMb:
+    case kPyMbWithDirectPhoton:
     case kPyMbNonDiffr:
     case kPyMbMSEL1:
     case kPyJets:
@@ -488,9 +497,9 @@ void AliGenPythiaPlus::Generate()
     
     fDecayer->ForceDecay();
 
-    Float_t polar[3]   =   {0,0,0};
-    Float_t origin[3]  =   {0,0,0};
-    Float_t p[4];
+    Double_t polar[3]   =   {0,0,0};
+    Double_t origin[3]  =   {0,0,0};
+    Double_t p[4];
 //  converts from mm/c to s
     const Float_t kconv=0.001/2.999792458e8;
 //
@@ -526,8 +535,8 @@ void AliGenPythiaPlus::Generate()
            }
            fNpartons = fPythia->GetNumberOfParticles();
        } else {
-           printf("Loading Event %d\n",AliRunLoader::GetRunLoader()->GetEventNumber());
-           fRL->GetEvent(AliRunLoader::GetRunLoader()->GetEventNumber());
+           printf("Loading Event %d\n",AliRunLoader::Instance()->GetEventNumber());
+           fRL->GetEvent(AliRunLoader::Instance()->GetEventNumber());
            fPythia->SetNumberOfParticles(0);
            fPythia->LoadEvent(fRL->Stack(), 0 , 1);
            fPythia->EditEventList(21);
@@ -553,7 +562,7 @@ void AliGenPythiaPlus::Generate()
 //     printf("Calling hadronisation %d\n", fPythia->GetN());
 //     fPythia->HadronizeEvent();      
        fTrials++;
-       fPythia->GetParticles(fParticles);
+       fPythia->GetParticles(&fParticles);
        Boost();
 //
 //
@@ -561,7 +570,7 @@ void AliGenPythiaPlus::Generate()
        Int_t i;
        
        fNprimaries = 0;
-       Int_t np = fParticles->GetEntriesFast();
+       Int_t np = fParticles.GetEntriesFast();
        
        if (np == 0) continue;
 //
@@ -581,6 +590,7 @@ void AliGenPythiaPlus::Generate()
        Int_t nTkbles = 0;   // Trackable particles
        if (fProcess != kPyMbDefault && 
            fProcess != kPyMb && 
+           fProcess != kPyMbWithDirectPhoton && 
            fProcess != kPyJets && 
            fProcess != kPyDirectGamma &&
            fProcess != kPyMbNonDiffr  &&
@@ -591,7 +601,7 @@ void AliGenPythiaPlus::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
@@ -618,21 +628,21 @@ void AliGenPythiaPlus::Generate()
                if (kf >= fFlavorSelect && kf <= 6) {
                    Int_t idau = (fPythia->Version() == 6) ? (iparticle->GetFirstDaughter() - 1) :(iparticle->GetFirstDaughter());
                    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 = (fPythia->Version() == 6) ? (daughter->GetFirstDaughter() - 1) : (daughter->GetFirstDaughter());
                            Int_t jmax = (fPythia->Version() == 6) ? (daughter->GetLastDaughter() - 1)  : (daughter->GetLastDaughter());
 
-                           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());
                }
                
@@ -726,7 +736,7 @@ void AliGenPythiaPlus::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();
@@ -738,7 +748,7 @@ void AliGenPythiaPlus::Generate()
                    origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
                    origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
                    
-                   Float_t tof   = kconv*iparticle->T();
+                   Float_t tof   = fTime + kconv*iparticle->T();
                    Int_t ipa = (fPythia->Version() == 6) ? (iparticle->GetFirstMother() - 1) :(iparticle->GetFirstMother()) ;
                    Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
  
@@ -807,13 +817,13 @@ Int_t  AliGenPythiaPlus::GenerateMB()
 //  converts from mm/c to s
     const Float_t kconv = 0.001 / 2.999792458e8;
     
-    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;
@@ -830,11 +840,11 @@ Int_t  AliGenPythiaPlus::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 = (fPythia->Version() == 6) ? (iparticle->GetFirstMother() - 1) :(iparticle->GetFirstMother()) ;
-         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
            {
@@ -846,8 +856,10 @@ Int_t  AliGenPythiaPlus::GenerateMB()
            }
        }
       }
-      if(!ok)
-       return 0;
+      if(!ok){
+         delete [] pParent;
+         return 0;
+      }
     }
     
     
@@ -859,7 +871,7 @@ Int_t  AliGenPythiaPlus::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 
         
@@ -881,14 +893,16 @@ Int_t  AliGenPythiaPlus::GenerateMB()
       if(!okd && iphcand != -1) // execute rotation in phi 
           RotatePhi(iphcand,okd);
       
-      if(!okd)
-       return 0;
+      if(!okd) {
+         delete[] pParent;
+         return 0;
+      }
     }
     
     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;
@@ -906,23 +920,45 @@ Int_t  AliGenPythiaPlus::GenerateMB()
     // Check if there is a ccbar or bbbar pair with at least one of the two
     // in fYMin < y < fYMax
     if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi) {
-      TParticle *hvq;
+      TParticle *partCheck;
+      TParticle *mother;
       Bool_t  theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
-      Float_t yQ;  
-      Int_t   pdgQ;
+      Bool_t  theChild=kFALSE;
+      Float_t y;  
+      Int_t   pdg,mpdg,mpdgUpperFamily;
       for(i=0; i<np; i++) {
-       hvq = (TParticle*)fParticles->At(i);
-       pdgQ = hvq->GetPdgCode();  
-       if(TMath::Abs(pdgQ) != fFlavorSelect) continue; 
-       if(pdgQ>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
-       yQ = 0.5*TMath::Log((hvq->Energy()+hvq->Pz()+1.e-13)/
-                           (hvq->Energy()-hvq->Pz()+1.e-13));
-       if(yQ>fYMin && yQ<fYMax) inYcut=kTRUE;
+       partCheck = (TParticle*)fParticles.At(i);
+       pdg = partCheck->GetPdgCode();  
+       if(TMath::Abs(pdg) == fFlavorSelect) { // quark  
+         if(pdg>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
+         y = 0.5*TMath::Log((partCheck->Energy()+partCheck->Pz()+1.e-13)/
+                            (partCheck->Energy()-partCheck->Pz()+1.e-13));
+         if(y>fYMin && y<fYMax) inYcut=kTRUE;
+       }
+
+       if(fCutOnChild && TMath::Abs(pdg) == fPdgCodeParticleforAcceptanceCut) {
+         Int_t mi = partCheck->GetFirstMother() - 1;
+         if(mi<0) continue;
+         mother = (TParticle*)fParticles.At(mi);
+         mpdg=TMath::Abs(mother->GetPdgCode());
+         mpdgUpperFamily=(mpdg>1000 ? mpdg+1000 : mpdg+100); // keep e from c from b
+         if ( ParentSelected(mpdg) || 
+             (fFlavorSelect==5 && ParentSelected(mpdgUpperFamily))) {
+           if (KinematicSelection(partCheck,1)) {
+             theChild=kTRUE;
+           }
+         }
+       }
       }
-      if (!theQ || !theQbar || !inYcut) {
+      if (!theQ || !theQbar || !inYcut) { // one of the c/b conditions not satisfied
        delete[] pParent;
        return 0;
       }
+      if (fCutOnChild && !theChild) { // one of the child conditions not satisfied
+       delete[] pParent;
+       return 0;       
+      }
+
     }
 
     //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
@@ -930,6 +966,7 @@ Int_t  AliGenPythiaPlus::GenerateMB()
          fProcess == kPyZ ||
          fProcess == kPyMbDefault ||
          fProcess == kPyMb ||
+         fProcess == kPyMbWithDirectPhoton ||
          fProcess == kPyMbNonDiffr)  
         && (fCutOnChild == 1) ) {
       if ( !CheckKinematicsOnChild() ) {
@@ -941,7 +978,7 @@ Int_t  AliGenPythiaPlus::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();
@@ -967,7 +1004,7 @@ Int_t  AliGenPythiaPlus::GenerateMB()
            origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
            origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
            
-           Float_t tof = fEventTime + kconv * iparticle->T();
+           Float_t tof = fTime + fEventTime + kconv * iparticle->T();
 
            PushTrack(fTrackIt*trackIt, iparent, kf, 
                      p[0], p[1], p[2], p[3], 
@@ -980,17 +1017,18 @@ Int_t  AliGenPythiaPlus::GenerateMB()
            //
            // Special Treatment to store color-flow
            //
+           /*
            if (ks == 3 || ks == 13 || ks == 14) {
                TParticle* particle = 0;
                if (fStack) {
                    particle = fStack->Particle(nt);
                } else {
-                   particle = gAlice->Stack()->Particle(nt);
+                   particle = AliRunLoader::Instance()->Stack()->Particle(nt);
                }
-//             particle->SetFirstDaughter(fPythia->GetK(2, i));
-//             particle->SetLastDaughter(fPythia->GetK(3, i));         
+               particle->SetFirstDaughter(fPythia->GetK(2, i));
+               particle->SetLastDaughter(fPythia->GetK(3, i));         
            }
-           
+           */  
            KeepTrack(nt);
            pParent[i] = nt;
            SetHighWaterMark(nt);
@@ -1065,7 +1103,7 @@ void AliGenPythiaPlus::MakeHeader()
 //
 // Event Vertex 
     fHeader->SetPrimaryVertex(fVertex);
-    
+    fHeader->SetInteractionTime(fTime+fEventTime);    
 //
 // Number of primaries
     fHeader->SetNProduced(fNprimaries);
@@ -1104,8 +1142,9 @@ void AliGenPythiaPlus::MakeHeader()
 // Store quenching parameters
 //
     if (fQuench){
-       Double_t z[4];
-       Double_t xp, yp;
+        Double_t z[4] = {0.};
+       Double_t xp = 0.;
+       Double_t yp = 0.;
        if (fQuench == 1) {
            // Pythia::Quench()
            fPythia->GetQuenchingParameters(xp, yp, z);
@@ -1133,7 +1172,7 @@ void AliGenPythiaPlus::MakeHeader()
     fHeader = 0x0;
 }
 
-Bool_t AliGenPythiaPlus::CheckTrigger(TParticle* jet1, TParticle* jet2)
+Bool_t AliGenPythiaPlus::CheckTrigger(const TParticle* jet1, const TParticle* jet2)
 {
 // Check the kinematic trigger condition
 //
@@ -1192,10 +1231,10 @@ Bool_t AliGenPythiaPlus::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(); 
@@ -1302,7 +1341,7 @@ void AliGenPythiaPlus::GetSubEventTime()
 
 
 
-Bool_t AliGenPythiaPlus::IsInEMCAL(Float_t phi, Float_t eta)
+Bool_t AliGenPythiaPlus::IsInEMCAL(Float_t phi, Float_t eta) const
 {
   // Is particle in EMCAL acceptance? 
   // phi in degrees, etamin=-etamax
@@ -1313,7 +1352,7 @@ Bool_t AliGenPythiaPlus::IsInEMCAL(Float_t phi, Float_t eta)
     return kFALSE;
 }
 
-Bool_t AliGenPythiaPlus::IsInPHOS(Float_t phi, Float_t eta)
+Bool_t AliGenPythiaPlus::IsInPHOS(Float_t phi, Float_t eta) const
 {
   // Is particle in PHOS acceptance? 
   // Acceptance slightly larger considered.
@@ -1333,19 +1372,19 @@ void AliGenPythiaPlus::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