Changes related to the initialization of random numbers generators. Now one can use...
[u/mrichter/AliRoot.git] / THerwig / AliGenHerwig.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 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"
24 #include "AliRun.h"
25
26 #include <TParticle.h>
27 #include "THerwig6.h"
28
29 #include "Riostream.h"
30
31  ClassImp(AliGenHerwig)
32
33 static TRandom * sRandom;
34
35 AliGenHerwig::AliGenHerwig()
36 {
37 // Constructor
38 }
39
40 AliGenHerwig::AliGenHerwig(Int_t npart)
41     :AliGenMC(npart)
42 {
43     SetBeamMomenta();
44     SetTarget();
45     SetProjectile();
46     SetStrucFunc(kGRVLO98);
47     fKeep=0;
48     fTrigger=0;
49     fDecaysOff=1;
50     fSelectAll=0;
51     fFlavor=0;
52     fPtHardMin=10.;
53     fPtRMS=0.0;
54     fEnSoft=1.0;
55     fMaxPr=1;
56     fMaxErrors=1000;
57 //  Set random number
58     if (!sRandom) sRandom=fRandom;
59 }
60
61 AliGenHerwig::AliGenHerwig(const AliGenHerwig & Herwig)
62 {
63 // copy constructor
64 }
65
66
67 AliGenHerwig::~AliGenHerwig()
68 {
69 // Destructor
70 }
71
72 void AliGenHerwig::Init()
73 {
74 // Initialisation
75   fTarget.Resize(8);
76   fProjectile.Resize(8);
77   SetMC(new THerwig6());
78   fHerwig=(THerwig6*) fgMCEvGen;
79   // initialize common blocks
80   fHerwig->Initialize(fProjectile, fTarget, fMomentum1, fMomentum2, fProcess);
81   // reset parameters according to user needs
82   InitPDF();
83   fHerwig->SetPTMIN(fPtHardMin);
84   fHerwig->SetPTRMS(fPtRMS);
85   fHerwig->SetMAXPR(fMaxPr);
86   fHerwig->SetMAXER(fMaxErrors);
87   fHerwig->SetENSOF(fEnSoft);
88   // compute parameter dependent constants
89   fHerwig->PrepareRun();
90 }
91
92 void AliGenHerwig::InitPDF()
93 {
94   switch(fStrucFunc)
95     {
96     case kGRVLO:
97       fModPDF=5;
98       fAutPDF="GRV";
99       break;
100     case kGRVHO:
101       fModPDF=6;
102       fAutPDF="GRV";
103       break;
104     case kGRVLO98:
105       fModPDF=12;
106       fAutPDF="GRV";
107       break;
108     case kMRSDminus:
109       fModPDF=31;
110       fAutPDF="MRS";
111       break;
112     case kMRSD0:
113       fModPDF=30;
114       fAutPDF="MRS";
115       break;
116     case kMRSG:
117       fModPDF=41;
118       fAutPDF="MRS";
119       break;
120     case kMRSTcgLO:
121       fModPDF=72;
122       fAutPDF="MRS";
123       break;
124     case kCTEQ4M:
125       fModPDF=34;
126       fAutPDF="CTEQ";
127       break;
128     case kCTEQ5L:
129       fModPDF=46;
130       fAutPDF="CTEQ";
131       break;
132     default:
133       cerr << "This structure function is not inplemented " << fStrucFunc << endl;
134       break;
135     }
136   fAutPDF.Resize(20);      
137   fHerwig->SetMODPDF(1,fModPDF);
138   fHerwig->SetMODPDF(2,fModPDF);
139   fHerwig->SetAUTPDF(1,fAutPDF);
140   fHerwig->SetAUTPDF(2,fAutPDF);
141 }
142
143 void AliGenHerwig::Generate()
144 {
145   // Generate one event
146
147   Float_t polar[3] =   {0,0,0};
148   Float_t origin[3]=   {0,0,0};
149   Float_t origin0[3]=  {0,0,0};
150   Float_t p[4], random[6];
151   
152   static TClonesArray *particles;
153   //  converts from mm/c to s
154   const Float_t kconv=0.001/2.999792458e8;
155   //
156   Int_t nt=0;
157   Int_t jev=0;
158   Int_t j, kf, ks, imo;
159   kf=0;
160   
161   if(!particles) particles=new TClonesArray("TParticle",10000);
162   
163   fTrials=0;
164   for (j=0;j<3;j++) origin0[j]=fOrigin[j];
165   if(fVertexSmear==kPerEvent) {
166     Rndm(random,6);
167     for (j=0;j<3;j++) {
168       origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
169         TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
170     }
171   }
172   
173   while(1)
174     {
175         fHerwig->GenerateEvent();
176         fTrials++;
177         fHerwig->ImportParticles(particles,"All");
178         Int_t np = particles->GetEntriesFast()-1;
179         if (np == 0 ) continue;
180         
181         Int_t nc=0;
182         
183         Int_t * newPos = new Int_t[np];
184         for (Int_t i = 0; i<np; i++) *(newPos+i)=-1;
185         
186         for (Int_t i = 0; i<np; i++) {
187             TParticle *  iparticle       = (TParticle *) particles->At(i);
188             imo = iparticle->GetFirstMother();
189             kf        = iparticle->GetPdgCode();
190             ks        = iparticle->GetStatusCode();
191             if (ks != 3 && 
192                 KinematicSelection(iparticle,0)) 
193             {
194                 nc++;
195                 p[0]=iparticle->Px();
196                 p[1]=iparticle->Py();
197                 p[2]=iparticle->Pz();
198                 p[3]=iparticle->Energy();
199                 origin[0]=origin0[0]+iparticle->Vx()/10;
200                 origin[1]=origin0[1]+iparticle->Vy()/10;
201                 origin[2]=origin0[2]+iparticle->Vz()/10;
202                 Float_t tof = kconv*iparticle->T();
203                 Int_t   iparent = (imo > -1) ? newPos[imo] : -1;
204                 Int_t   trackIt = (ks == 1) && fTrackIt;
205                 gAlice->SetTrack(trackIt, iparent, kf,
206                                  p[0], p[1], p[2], p[3],
207                                  origin[0], origin[1], origin[2], 
208                                  tof,
209                                  polar[0], polar[1], polar[2],
210                                  kPPrimary, nt, 1., ks);
211                 KeepTrack(nt);
212                 newPos[i]=nt;
213             } // end of if: selection of particle
214         } // end of for: particle loop
215         if (newPos) delete[] newPos;
216         printf("\n I've put %i particles on the stack \n",nc);
217         //      MakeHeader();
218         printf("nc: %d %d\n", nc, fNpart);
219         
220         if (nc > 0) {
221             jev+=nc;
222             if (jev >= fNpart || fNpart == -1) {
223                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
224                 printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
225                 break;
226             }
227         }
228     }
229   SetHighWaterMark(nt);
230 //  adjust weight due to kinematic selection
231   AdjustWeights();
232 //  get cross-section
233   fXsection=fHerwig->GetAVWGT();
234 }
235
236 void AliGenHerwig::AdjustWeights()
237 {
238 // Adjust the weights after generation of all events
239     TParticle *part;
240     Int_t ntrack=gAlice->GetNtrack();
241     for (Int_t i=0; i<ntrack; i++) {
242         part= gAlice->Particle(i);
243         part->SetWeight(part->GetWeight()*fKineBias);
244     }
245 }
246  
247
248 void AliGenHerwig::KeepFullEvent()
249 {
250     fKeep=1;
251 }
252
253 Bool_t AliGenHerwig::DaughtersSelection(TParticle* iparticle, TClonesArray* particles)
254 {
255 //
256 // Looks recursively if one of the daughters has been selected
257 //
258 //    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
259     Int_t imin=-1;
260     Int_t imax=-1;
261     Int_t i;
262     Bool_t hasDaughters= (iparticle->GetFirstDaughter() >=0);
263     Bool_t selected=kFALSE;
264     if (hasDaughters) {
265         imin=iparticle->GetFirstDaughter();
266         imax=iparticle->GetLastDaughter();       
267         for (i=imin; i<= imax; i++){
268             TParticle *  jparticle       = (TParticle *) particles->At(i);      
269             Int_t ip=jparticle->GetPdgCode();
270             if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
271                 selected=kTRUE; break;
272             }
273             if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
274         }
275     } else {
276         return kFALSE;
277     }
278
279     return selected;
280 }
281
282
283 Bool_t AliGenHerwig::SelectFlavor(Int_t pid)
284 {
285 // Select flavor of particle
286 // 0: all
287 // 4: charm and beauty
288 // 5: beauty
289     if (fFlavor == 0) return kTRUE;
290     
291     Int_t ifl=TMath::Abs(pid/100);
292     if (ifl > 10) ifl/=10;
293     return (fFlavor == ifl);
294 }
295
296 Bool_t AliGenHerwig::Stable(TParticle*  particle)
297 {
298 // Return true for a stable particle
299 //
300     Int_t kf = TMath::Abs(particle->GetPdgCode());
301     
302     if ( (particle->GetFirstDaughter() < 0 ) || (kf == 1000*fFlavor+122))
303          
304     {
305         return kTRUE;
306     } else {
307         return kFALSE;
308     }
309 }
310
311 void AliGenHerwig::FinishRun()
312 {
313   fHerwig->Hwefin();
314 }
315
316
317 AliGenHerwig& AliGenHerwig::operator=(const  AliGenHerwig& rhs)
318 {
319 // Assignment operator
320     return *this;
321 }
322
323
324 extern "C" {
325   Double_t hwr_() {return sRandom->Rndm();}
326 }