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