EffC++ warning corrected.
[u/mrichter/AliRoot.git] / THerwig / AliGenHerwig.cxx
CommitLineData
618e1dc0 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
7cdba479 16/* $Id$ */
618e1dc0 17
18
19// Generator using Herwig as an external generator
20// The main Herwig options are accessable for the user through this interface.
21// Uses the THerwig implementation of TGenerator.
22
23#include "AliGenHerwig.h"
93831bcf 24#include "AliHerwigRndm.h"
618e1dc0 25#include "AliRun.h"
26
27#include <TParticle.h>
28#include "THerwig6.h"
29
30#include "Riostream.h"
5d12ce38 31#include "AliMC.h"
618e1dc0 32
e2054d85 33#include "driver.h"
34
014a9521 35ClassImp(AliGenHerwig)
618e1dc0 36
7cdba479 37
a1d47d99 38 AliGenHerwig::AliGenHerwig() :
39 AliGenMC(),
40 fAutPDF("GRV"),
41 fModPDF(5),
2ad1e1c2 42 fStrucFunc(kCTEQ5L),
a1d47d99 43 fKeep(0),
44 fDecaysOff(1),
45 fTrigger(0),
46 fSelectAll(0),
47 fFlavor(0),
48 fEnergyCMS(14000),
49 fMomentum1(7000),
50 fMomentum2(7000),
51 fKineBias(1),
52 fTrials(0),
53 fXsection(0),
54 fHerwig(0x0),
55 fProcess(0),
56 fPtHardMin(0),
57 fPtRMS(0),
58 fMaxPr(10),
59 fMaxErrors(1000),
e2054d85 60 fEnSoft(1),
61 fEv1Pr(0),
62 fEv2Pr(0)
618e1dc0 63{
64// Constructor
65}
66
67AliGenHerwig::AliGenHerwig(Int_t npart)
68 :AliGenMC(npart)
69{
70 SetBeamMomenta();
71 SetTarget();
72 SetProjectile();
2ad1e1c2 73 SetStrucFunc(kCTEQ5L);
618e1dc0 74 fKeep=0;
75 fTrigger=0;
76 fDecaysOff=1;
77 fSelectAll=0;
78 fFlavor=0;
93831bcf 79 fPtHardMin=0.;
618e1dc0 80 fPtRMS=0.0;
81 fEnSoft=1.0;
93831bcf 82 fMaxPr=10;
618e1dc0 83 fMaxErrors=1000;
e2054d85 84 fEv1Pr=0;
85 fEv2Pr=0;
93831bcf 86 // Set random number generator
87 AliHerwigRndm::SetHerwigRandom(GetRandom());
618e1dc0 88}
89
90AliGenHerwig::AliGenHerwig(const AliGenHerwig & Herwig)
014a9521 91 :AliGenMC(Herwig)
618e1dc0 92{
014a9521 93// Copy constructor
94 Herwig.Copy(*this);
618e1dc0 95}
96
97
98AliGenHerwig::~AliGenHerwig()
99{
100// Destructor
101}
102
e2054d85 103void AliGenHerwig::SetEventListRange(Int_t eventFirst, Int_t eventLast)
104{
105 fEv1Pr = ++eventFirst;
106 fEv2Pr = ++eventLast;
107 if ( fEv2Pr == -1 ) fEv2Pr = fEv2Pr;
108}
109
618e1dc0 110void AliGenHerwig::Init()
111{
112// Initialisation
113 fTarget.Resize(8);
114 fProjectile.Resize(8);
115 SetMC(new THerwig6());
d274cd60 116 fHerwig=(THerwig6*) fMCEvGen;
618e1dc0 117 // initialize common blocks
118 fHerwig->Initialize(fProjectile, fTarget, fMomentum1, fMomentum2, fProcess);
119 // reset parameters according to user needs
120 InitPDF();
121 fHerwig->SetPTMIN(fPtHardMin);
122 fHerwig->SetPTRMS(fPtRMS);
123 fHerwig->SetMAXPR(fMaxPr);
124 fHerwig->SetMAXER(fMaxErrors);
125 fHerwig->SetENSOF(fEnSoft);
e2054d85 126
127 fHerwig->SetEV1PR(fEv1Pr);
128 fHerwig->SetEV2PR(fEv2Pr);
129
130// C---D,U,S,C,B,T QUARK AND GLUON MASSES (IN THAT ORDER)
131// RMASS(1)=0.32
132// RMASS(2)=0.32
133// RMASS(3)=0.5
134// RMASS(4)=1.55
135// RMASS(5)=4.95
136// RMASS(6)=174.3
137// RMASS(13)=0.75
138
139 fHerwig->SetRMASS(4,1.2);
140 fHerwig->SetRMASS(5,4.75);
141
142 if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(),50);
143
144 fHerwig->Hwusta("PI0 ");
145
618e1dc0 146 // compute parameter dependent constants
147 fHerwig->PrepareRun();
148}
149
150void AliGenHerwig::InitPDF()
151{
152 switch(fStrucFunc)
153 {
618e1dc0 154 case kGRVLO98:
155 fModPDF=12;
156 fAutPDF="GRV";
157 break;
618e1dc0 158 case kCTEQ4M:
159 fModPDF=34;
160 fAutPDF="CTEQ";
161 break;
162 case kCTEQ5L:
163 fModPDF=46;
164 fAutPDF="CTEQ";
165 break;
e2054d85 166 case kCTEQ5M:
167 fModPDF=48;
168 fAutPDF="CTEQ";
169 break;
618e1dc0 170 default:
171 cerr << "This structure function is not inplemented " << fStrucFunc << endl;
172 break;
173 }
174 fAutPDF.Resize(20);
175 fHerwig->SetMODPDF(1,fModPDF);
176 fHerwig->SetMODPDF(2,fModPDF);
177 fHerwig->SetAUTPDF(1,fAutPDF);
178 fHerwig->SetAUTPDF(2,fAutPDF);
179}
180
181void AliGenHerwig::Generate()
182{
183 // Generate one event
184
185 Float_t polar[3] = {0,0,0};
186 Float_t origin[3]= {0,0,0};
187 Float_t origin0[3]= {0,0,0};
188 Float_t p[4], random[6];
189
190 static TClonesArray *particles;
191 // converts from mm/c to s
192 const Float_t kconv=0.001/2.999792458e8;
193 //
194 Int_t nt=0;
195 Int_t jev=0;
196 Int_t j, kf, ks, imo;
197 kf=0;
198
199 if(!particles) particles=new TClonesArray("TParticle",10000);
200
201 fTrials=0;
202 for (j=0;j<3;j++) origin0[j]=fOrigin[j];
203 if(fVertexSmear==kPerEvent) {
204 Rndm(random,6);
205 for (j=0;j<3;j++) {
206 origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
207 TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
208 }
209 }
210
211 while(1)
212 {
213 fHerwig->GenerateEvent();
214 fTrials++;
215 fHerwig->ImportParticles(particles,"All");
216 Int_t np = particles->GetEntriesFast()-1;
217 if (np == 0 ) continue;
218
219 Int_t nc=0;
220
221 Int_t * newPos = new Int_t[np];
222 for (Int_t i = 0; i<np; i++) *(newPos+i)=-1;
223
224 for (Int_t i = 0; i<np; i++) {
225 TParticle * iparticle = (TParticle *) particles->At(i);
226 imo = iparticle->GetFirstMother();
227 kf = iparticle->GetPdgCode();
228 ks = iparticle->GetStatusCode();
229 if (ks != 3 &&
230 KinematicSelection(iparticle,0))
231 {
232 nc++;
233 p[0]=iparticle->Px();
234 p[1]=iparticle->Py();
235 p[2]=iparticle->Pz();
236 p[3]=iparticle->Energy();
237 origin[0]=origin0[0]+iparticle->Vx()/10;
238 origin[1]=origin0[1]+iparticle->Vy()/10;
239 origin[2]=origin0[2]+iparticle->Vz()/10;
240 Float_t tof = kconv*iparticle->T();
241 Int_t iparent = (imo > -1) ? newPos[imo] : -1;
242 Int_t trackIt = (ks == 1) && fTrackIt;
93831bcf 243 PushTrack(trackIt, iparent, kf,
244 p[0], p[1], p[2], p[3],
245 origin[0], origin[1], origin[2],
246 tof,
247 polar[0], polar[1], polar[2],
e2054d85 248 kPPrimary, nt, fHerwig->GetEVWGT(), ks);
618e1dc0 249 KeepTrack(nt);
250 newPos[i]=nt;
251 } // end of if: selection of particle
252 } // end of for: particle loop
253 if (newPos) delete[] newPos;
618e1dc0 254 // MakeHeader();
618e1dc0 255 if (nc > 0) {
256 jev+=nc;
257 if (jev >= fNpart || fNpart == -1) {
258 fKineBias=Float_t(fNpart)/Float_t(fTrials);
618e1dc0 259 break;
260 }
261 }
262 }
263 SetHighWaterMark(nt);
264// adjust weight due to kinematic selection
265 AdjustWeights();
266// get cross-section
267 fXsection=fHerwig->GetAVWGT();
268}
269
270void AliGenHerwig::AdjustWeights()
271{
272// Adjust the weights after generation of all events
273 TParticle *part;
5d12ce38 274 Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
618e1dc0 275 for (Int_t i=0; i<ntrack; i++) {
5d12ce38 276 part= gAlice->GetMCApp()->Particle(i);
618e1dc0 277 part->SetWeight(part->GetWeight()*fKineBias);
278 }
279}
280
281
282void AliGenHerwig::KeepFullEvent()
283{
284 fKeep=1;
285}
286
287Bool_t AliGenHerwig::DaughtersSelection(TParticle* iparticle, TClonesArray* particles)
288{
289//
290// Looks recursively if one of the daughters has been selected
291//
292// printf("\n Consider daughters %d:",iparticle->GetPdgCode());
293 Int_t imin=-1;
294 Int_t imax=-1;
295 Int_t i;
296 Bool_t hasDaughters= (iparticle->GetFirstDaughter() >=0);
297 Bool_t selected=kFALSE;
298 if (hasDaughters) {
299 imin=iparticle->GetFirstDaughter();
300 imax=iparticle->GetLastDaughter();
301 for (i=imin; i<= imax; i++){
302 TParticle * jparticle = (TParticle *) particles->At(i);
303 Int_t ip=jparticle->GetPdgCode();
304 if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
305 selected=kTRUE; break;
306 }
307 if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
308 }
309 } else {
310 return kFALSE;
311 }
312
313 return selected;
314}
315
316
317Bool_t AliGenHerwig::SelectFlavor(Int_t pid)
318{
319// Select flavor of particle
320// 0: all
321// 4: charm and beauty
322// 5: beauty
323 if (fFlavor == 0) return kTRUE;
324
325 Int_t ifl=TMath::Abs(pid/100);
326 if (ifl > 10) ifl/=10;
327 return (fFlavor == ifl);
328}
329
330Bool_t AliGenHerwig::Stable(TParticle* particle)
331{
332// Return true for a stable particle
333//
334 Int_t kf = TMath::Abs(particle->GetPdgCode());
335
336 if ( (particle->GetFirstDaughter() < 0 ) || (kf == 1000*fFlavor+122))
337
338 {
339 return kTRUE;
340 } else {
341 return kFALSE;
342 }
343}
344
345void AliGenHerwig::FinishRun()
346{
347 fHerwig->Hwefin();
348}
349
350
351AliGenHerwig& AliGenHerwig::operator=(const AliGenHerwig& rhs)
352{
353// Assignment operator
014a9521 354 rhs.Copy(*this);
355 return (*this);
618e1dc0 356}
357
358