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