EVGEN/AliGenPythiaJets.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /*
  17 $Log$
  18 */
  19
  20 //
  21 // Generator using the TPythia interface (via AliPythia)
  22 // to generate pp collisions.
  23 // Using SetNuclei() also nuclear modifications to the structure functions
  24 // can be taken into account. This makes, of course, only sense for the
  25 // generation of the products of hard processes (heavy flavor, jets ...)
  26 //
  27 // andreas.morsch@cern.ch
  28 //
  29
  30
  31 #include "AliGenPythiaJets.h"
  32 #include "AliRun.h"
  33 ClassImp(AliGenPythiaJets)
  34
  35 AliGenPythiaJets::AliGenPythiaJets()
  36                  :AliGenPythia()
  37 {
  38 // Default Constructor
  39 }
  40
  41 AliGenPythiaJets::AliGenPythiaJets(Int_t npart)
  42                  :AliGenPythia(npart)
  43 {
  44     fName = "PythiaJets";
  45     fTitle= "Jet Generator using PYTHIA";
  46 }
  47
  48 AliGenPythiaJets::AliGenPythiaJets(const AliGenPythiaJets & Pythia)
  49 {
  50 // copy constructor
  51     Pythia.Copy(*this);
  52 }
  53
  54 AliGenPythiaJets::~AliGenPythiaJets()
  55 {
  56 // Destructor
  57 }
  58
  59 void AliGenPythiaJets::Init()
  60 {
  61 // Initialization
  62 //
  63     printf("AliGenPythiaJets::Init() \n");
  64
  65     AliGenPythia::Init();
  66
  67     if (fQuench > 0.) {
  68          fEtMinJetQ[0]  = fEtMinJet  * fQuench;
  69          fEtMaxJetQ[0]  = fEtMaxJet  * fQuench;
  70          fEtMinJetQ[1]  = fEtMinJet  * (1. - fQuench);
  71          fEtMaxJetQ[1]  = fEtMaxJet  * (1. - fQuench);
  72          fPtHardMinQ[0] = fPtHardMin * fQuench;
  73          fPtHardMaxQ[0] = fPtHardMax * fQuench;
  74          fPtHardMinQ[1] = fPtHardMin * (1. - fQuench);
  75          fPtHardMaxQ[1] = fPtHardMax * (1. - fQuench);
  76     }
  77 }
  78
  79 void AliGenPythiaJets::Generate()
  80 {
  81 // Generate one event
  82     fDecayer->ForceDecay();
  83
  84     Float_t polar[3]   =   {0,0,0};
  85     Float_t origin[3]  =   {0,0,0};
  86     Float_t p[3];
  87 //  converts from mm/c to s
  88     const Float_t kconv=0.001/2.999792458e8;
  89 //
  90     Int_t nt  = 0;
  91     Int_t nc  = 0;
  92     Int_t jev = 0;
  93     Int_t j, kf, iparent;
  94     fTrials=0;
  95 //
  96 //  Set collision vertex position
  97     if(fVertexSmear==kPerEvent) {
  98         fPythia->SetMSTP(151,1);
  99         for (j=0;j<3;j++) {
 100             fPythia->SetPARP(151+j, fOsigma[j]*10.);
 101         }
 102     } else if (fVertexSmear==kPerTrack) {
 103         fPythia->SetMSTP(151,0);
 104     }
 105 //  Event loop
 106     while(1)
 107     {
 108         if (fQuench > 0.) {
 109             fPythia->SetCKIN(3,fPtHardMinQ[jev]);
 110             fPythia->SetCKIN(4,fPtHardMaxQ[jev]);
 111             fEtMinJet = fEtMinJetQ[jev];
 112             fEtMaxJet = fEtMaxJetQ[jev];
 113         }
 114
 115         fPythia->Pyevnt();
 116         if (gAlice->GetEvNumber()>=fDebugEventFirst &&
 117             gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1);
 118         fTrials++;
 119         //
 120         // Has this jet triggered
 121         //
 122         if ((fEtMinJet != -1) && ! CheckTrigger()) continue;
 123 //
 124         fPythia->ImportParticles(fParticles,"All");
 125         Int_t i;
 126         Int_t np = fParticles->GetEntriesFast();
 127         if (np == 0 ) continue;
 128 // Get event vertex and discard the event if the z coord. is too big
 129         TParticle *iparticle = (TParticle *) fParticles->At(0);
 130         Float_t distz = iparticle->Vz()/10.;
 131         if(TMath::Abs(distz)>fCutVertexZ*fOsigma[2]) continue;
 132 //
 133 //
 134         fEventVertex[0] = iparticle->Vx()/10.+fOrigin.At(0);
 135         fEventVertex[1] = iparticle->Vy()/10.+fOrigin.At(1);
 136         fEventVertex[2] = iparticle->Vz()/10.+fOrigin.At(2);
 137
 138         Int_t* pParent = new Int_t[np];
 139         for (i=0; i< np; i++) pParent[i] = -1;
 140
 141
 142         //
 143         for (i = 0; i<np; i++) {
 144             Int_t trackIt = 0;
 145             TParticle *  iparticle = (TParticle *) fParticles->At(i);
 146             kf = CheckPDGCode(iparticle->GetPdgCode());
 147             Int_t ks = iparticle->GetStatusCode();
 148             Int_t km = iparticle->GetFirstMother();
 149             if ((ks == 1  && kf !=0 && KinematicSelection(iparticle, 0)) ||
 150                 (ks != 1) ||
 151                 (fProcess == kPyJets && ks == 21 && km == 0 && i > 1)) {
 152                 nc++;
 153                 if (ks == 1) trackIt = 1;
 154                 Int_t ipa = iparticle->GetFirstMother() - 1;
 155                 iparent = (ipa > -1) ? pParent[ipa] : -1;
 156 //
 157 // Store track information
 158 //
 159                 p[0] = iparticle->Px();
 160                 p[1] = iparticle->Py();
 161                 p[2] = iparticle->Pz();
 162                 origin[0] = fOrigin[0]+iparticle->Vx()/10.;
 163                 origin[1] = fOrigin[1]+iparticle->Vy()/10.;
 164                 origin[2] = fOrigin[2]+iparticle->Vz()/10.;
 165                 Float_t tof=kconv*iparticle->T();
 166                 SetTrack(fTrackIt*trackIt, iparent, kf, p, origin, polar,
 167                          tof, kPPrimary, nt, 1., ks);
 168                 KeepTrack(nt);
 169                 pParent[i] = nt;
 170             } // select particle
 171         } // particle loop
 172
 173         if (pParent) delete[] pParent;
 174         printf("\n AliGenPythiaJets: I've put %i particles on the stack \n",nc);
 175         if (nc > 0) {
 176             jev += 1;
 177             if ((fQuench <= 0.) || (fQuench > 0. && jev == 2)) {
 178                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
 179                 printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
 180                 fNev++;
 181                 MakeHeader();
 182                 break;
 183             }
 184         }
 185     }
 186     SetHighWaterMark(nt);
 187 //  Get cross-section
 188     fXsection=fPythia->GetPARI(1);
 189 }
 190
 191 Bool_t AliGenPythiaJets::CheckTrigger()
 192 {
 193 // Check the kinematic trigger condition
 194 //
 195     Bool_t   triggered = kFALSE;
 196     Int_t njets = 0;
 197     Int_t ntrig = 0;
 198     Float_t jets[4][10];
 199 //
 200 // Use Pythia clustering on parton level to determine jet axis
 201 //
 202     GetJets(njets, ntrig, jets);
 203
 204     if (ntrig) {
 205         triggered = kTRUE;
 206         Float_t px   = jets[0][0];
 207         Float_t py   = jets[1][0];
 208         Float_t pz   = jets[2][0];
 209         Float_t e    = jets[3][0];
 210         Float_t beta = pz/e;
 211         Float_t phi  =  TMath::Pi()+TMath::ATan2(-py,-px);
 212         TransformEvent(beta, -TMath::Pi()/2. + phi);
 213     }
 214     return triggered;
 215 }
 216
 217 AliGenPythiaJets& AliGenPythiaJets::operator=(const  AliGenPythiaJets& rhs)
 218 {
 219 // Assignment operator
 220     return *this;
 221 }
 222
 223 void  AliGenPythiaJets::TransformEvent(Float_t beta, Float_t phi)
 224 {
 225 //
 226 // Perform Lorentz Transformation and Rotation
 227 //
 228     printf("Transforming \n");
 229
 230
 231     Float_t gamma = 1./TMath::Sqrt(1. - beta * beta);
 232     Int_t npart = (fPythia->GetPyjets())->N;
 233
 234     for (Int_t part = 0; part < npart; part++) {
 235         Float_t px   =  (fPythia->GetPyjets())->P[0][part];
 236         Float_t py   =  (fPythia->GetPyjets())->P[1][part];
 237         Float_t pz   =  (fPythia->GetPyjets())->P[2][part];
 238         Float_t e    =  (fPythia->GetPyjets())->P[3][part];
 239         //
 240         // Lorentz Transform
 241         //
 242         Float_t pzt =  gamma * pz        - gamma * beta * e;
 243         Float_t et  = -gamma * beta * pz + gamma        * e;
 244         //
 245         // Rotation
 246         //
 247         Float_t pxt =   TMath::Cos(phi) * px +  TMath::Sin(phi) * py;
 248         Float_t pyt = - TMath::Sin(phi) * px +  TMath::Cos(phi) * py;
 249         //
 250         //
 251         (fPythia->GetPyjets())->P[0][part] = pxt;
 252         (fPythia->GetPyjets())->P[1][part] = pyt;
 253         (fPythia->GetPyjets())->P[2][part] = pzt;
 254         (fPythia->GetPyjets())->P[3][part] = et;
 255     }
 256 }
 257