Getting AliGenHerwig in sync with rest of the framework.
[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 "AliHerwigRndm.h"
25 #include "AliRun.h"
26
27 #include <TParticle.h>
28 #include "THerwig6.h"
29
30 #include "Riostream.h"
31 #include "AliMC.h"
32
33 ClassImp(AliGenHerwig)
34
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=0.;
54     fPtRMS=0.0;
55     fEnSoft=1.0;
56     fMaxPr=10;
57     fMaxErrors=1000;
58     // Set random number generator   
59     AliHerwigRndm::SetHerwigRandom(GetRandom());
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         printf("Particles from Herwig %d \n", np);
183         
184         if (np == 0 ) continue;
185         
186         Int_t nc=0;
187         
188         Int_t * newPos = new Int_t[np];
189         for (Int_t i = 0; i<np; i++) *(newPos+i)=-1;
190         
191         for (Int_t i = 0; i<np; i++) {
192             TParticle *  iparticle       = (TParticle *) particles->At(i);
193             imo = iparticle->GetFirstMother();
194             kf        = iparticle->GetPdgCode();
195             ks        = iparticle->GetStatusCode();
196             if (ks != 3 && 
197                 KinematicSelection(iparticle,0)) 
198             {
199                 nc++;
200                 p[0]=iparticle->Px();
201                 p[1]=iparticle->Py();
202                 p[2]=iparticle->Pz();
203                 p[3]=iparticle->Energy();
204                 origin[0]=origin0[0]+iparticle->Vx()/10;
205                 origin[1]=origin0[1]+iparticle->Vy()/10;
206                 origin[2]=origin0[2]+iparticle->Vz()/10;
207                 Float_t tof = kconv*iparticle->T();
208                 Int_t   iparent = (imo > -1) ? newPos[imo] : -1;
209                 Int_t   trackIt = (ks == 1) && fTrackIt;
210                 PushTrack(trackIt, iparent, kf,
211                           p[0], p[1], p[2], p[3],
212                           origin[0], origin[1], origin[2], 
213                           tof,
214                           polar[0], polar[1], polar[2],
215                           kPPrimary, nt, 1., ks);
216                 KeepTrack(nt);
217                 newPos[i]=nt;
218             } // end of if: selection of particle
219         } // end of for: particle loop
220         if (newPos) delete[] newPos;
221         printf("\n I've put %i particles on the stack \n",nc);
222         //      MakeHeader();
223         printf("nc: %d %d\n", nc, fNpart);
224         
225         if (nc > 0) {
226             jev+=nc;
227             if (jev >= fNpart || fNpart == -1) {
228                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
229                 printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
230                 break;
231             }
232         }
233     }
234   SetHighWaterMark(nt);
235 //  adjust weight due to kinematic selection
236   AdjustWeights();
237 //  get cross-section
238   fXsection=fHerwig->GetAVWGT();
239 }
240
241 void AliGenHerwig::AdjustWeights()
242 {
243 // Adjust the weights after generation of all events
244     TParticle *part;
245     Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
246     for (Int_t i=0; i<ntrack; i++) {
247         part= gAlice->GetMCApp()->Particle(i);
248         part->SetWeight(part->GetWeight()*fKineBias);
249     }
250 }
251  
252
253 void AliGenHerwig::KeepFullEvent()
254 {
255     fKeep=1;
256 }
257
258 Bool_t AliGenHerwig::DaughtersSelection(TParticle* iparticle, TClonesArray* particles)
259 {
260 //
261 // Looks recursively if one of the daughters has been selected
262 //
263 //    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
264     Int_t imin=-1;
265     Int_t imax=-1;
266     Int_t i;
267     Bool_t hasDaughters= (iparticle->GetFirstDaughter() >=0);
268     Bool_t selected=kFALSE;
269     if (hasDaughters) {
270         imin=iparticle->GetFirstDaughter();
271         imax=iparticle->GetLastDaughter();       
272         for (i=imin; i<= imax; i++){
273             TParticle *  jparticle       = (TParticle *) particles->At(i);      
274             Int_t ip=jparticle->GetPdgCode();
275             if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
276                 selected=kTRUE; break;
277             }
278             if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
279         }
280     } else {
281         return kFALSE;
282     }
283
284     return selected;
285 }
286
287
288 Bool_t AliGenHerwig::SelectFlavor(Int_t pid)
289 {
290 // Select flavor of particle
291 // 0: all
292 // 4: charm and beauty
293 // 5: beauty
294     if (fFlavor == 0) return kTRUE;
295     
296     Int_t ifl=TMath::Abs(pid/100);
297     if (ifl > 10) ifl/=10;
298     return (fFlavor == ifl);
299 }
300
301 Bool_t AliGenHerwig::Stable(TParticle*  particle)
302 {
303 // Return true for a stable particle
304 //
305     Int_t kf = TMath::Abs(particle->GetPdgCode());
306     
307     if ( (particle->GetFirstDaughter() < 0 ) || (kf == 1000*fFlavor+122))
308          
309     {
310         return kTRUE;
311     } else {
312         return kFALSE;
313     }
314 }
315
316 void AliGenHerwig::FinishRun()
317 {
318   fHerwig->Hwefin();
319 }
320
321
322 AliGenHerwig& AliGenHerwig::operator=(const  AliGenHerwig& rhs)
323 {
324 // Assignment operator
325     rhs.Copy(*this);
326     return (*this);
327 }
328
329