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