Removing the fake copy constructors and assignment operator, moving their declaration...
[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 #include "driver.h"
34
35 ClassImp(AliGenHerwig)
36
37
38   AliGenHerwig::AliGenHerwig() :
39     AliGenMC(),
40     fAutPDF("LHAPDF"),
41     fModPDF(19070),
42     fStrucFunc(kCTEQ5L),
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),
60     fEnSoft(1),
61     fEv1Pr(0),
62     fEv2Pr(0),
63     fFileName(0)
64 {
65 // Constructor
66 }
67
68 AliGenHerwig::AliGenHerwig(Int_t npart)
69     :AliGenMC(npart),
70     fAutPDF("LHAPDF"),
71     fModPDF(19070),
72     fStrucFunc(kCTEQ5L),
73     fKeep(0),
74     fDecaysOff(1),
75     fTrigger(0),
76     fSelectAll(0),
77     fFlavor(0),
78     fEnergyCMS(14000),
79     fMomentum1(7000),
80     fMomentum2(7000),
81     fKineBias(1),
82     fTrials(0),
83     fXsection(0),
84     fHerwig(0x0),
85     fProcess(0),
86     fPtHardMin(0.),
87     fPtRMS(0.),
88     fMaxPr(10),
89     fMaxErrors(1000),
90     fEnSoft(1),
91     fEv1Pr(0),
92     fEv2Pr(0),
93     fFileName(0)
94 {
95     SetTarget();
96     SetProjectile();
97     // Set random number generator
98     AliHerwigRndm::SetHerwigRandom(GetRandom());
99 }
100
101 AliGenHerwig::~AliGenHerwig()
102 {
103 // Destructor
104 }
105
106 void AliGenHerwig::SetEventListRange(Int_t eventFirst, Int_t eventLast)
107 {
108   fEv1Pr = ++eventFirst;
109   fEv2Pr = ++eventLast;
110   if ( fEv2Pr == -1 ) fEv2Pr = fEv2Pr;
111 }
112
113 void AliGenHerwig::Init()
114 {
115 // Initialisation
116   fTarget.Resize(8);
117   fProjectile.Resize(8);
118   SetMC(new THerwig6());
119   fHerwig=(THerwig6*) fMCEvGen;
120   // initialize common blocks
121   fHerwig->Initialize(fProjectile, fTarget, fMomentum1, fMomentum2, fProcess);
122   // reset parameters according to user needs
123   InitPDF();
124   fHerwig->SetPTMIN(fPtHardMin);
125   fHerwig->SetPTRMS(fPtRMS);
126   fHerwig->SetMAXPR(fMaxPr);
127   fHerwig->SetMAXER(fMaxErrors);
128   fHerwig->SetENSOF(fEnSoft);
129
130   fHerwig->SetEV1PR(fEv1Pr);
131   fHerwig->SetEV2PR(fEv2Pr);
132
133 // C---D,U,S,C,B,T QUARK AND GLUON MASSES (IN THAT ORDER)
134 //       RMASS(1)=0.32
135 //       RMASS(2)=0.32
136 //       RMASS(3)=0.5
137 //       RMASS(4)=1.55
138 //       RMASS(5)=4.75
139 //       RMASS(6)=174.3
140 //       RMASS(13)=0.75
141
142   fHerwig->SetRMASS(4,1.2);
143   fHerwig->SetRMASS(5,4.75);
144
145   if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(),50);
146
147   fHerwig->Hwusta("PI0     ");
148
149   // compute parameter dependent constants
150   fHerwig->PrepareRun();
151 }
152
153 void AliGenHerwig::InitJimmy()
154 {
155 // Initialisation
156   fTarget.Resize(8);
157   fProjectile.Resize(8);
158   SetMC(new THerwig6());
159   fHerwig=(THerwig6*) fMCEvGen;
160   // initialize common blocks
161   fHerwig->InitializeJimmy(fProjectile, fTarget, fMomentum1, fMomentum2, fProcess);
162   // reset parameters according to user needs
163   InitPDF();
164   fHerwig->SetPTMIN(fPtHardMin);
165   fHerwig->SetPTRMS(fPtRMS);
166   fHerwig->SetMAXPR(fMaxPr);
167   fHerwig->SetMAXER(fMaxErrors);
168   fHerwig->SetENSOF(fEnSoft);
169
170   fHerwig->SetEV1PR(fEv1Pr);
171   fHerwig->SetEV2PR(fEv2Pr);
172
173 // C---D,U,S,C,B,T QUARK AND GLUON MASSES (IN THAT ORDER)
174 //       RMASS(1)=0.32
175 //       RMASS(2)=0.32
176 //       RMASS(3)=0.5
177 //       RMASS(4)=1.55
178 //       RMASS(5)=4.75
179 //       RMASS(6)=174.3
180 //       RMASS(13)=0.75
181
182   fHerwig->SetRMASS(4,1.2);
183   fHerwig->SetRMASS(5,4.75);
184
185   if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(),50);
186
187   fHerwig->Hwusta("PI0     ");
188
189   // compute parameter dependent constants
190   fHerwig->PrepareRunJimmy();
191 }
192
193 void AliGenHerwig::InitPDF()
194 {
195   switch(fStrucFunc)
196     {
197 // ONLY USES LHAPDF STRUCTURE FUNCTIONS
198     case kGRVLO98:
199       fModPDF=80060;
200       fAutPDF="HWLHAPDF";
201       break;
202     case kCTEQ6:
203       fModPDF=10040;
204       fAutPDF="HWLHAPDF";
205       break;
206     case kCTEQ61:
207       fModPDF=10100;
208       fAutPDF="HWLHAPDF";
209       break;
210     case kCTEQ6m:
211       fModPDF=10050;
212       fAutPDF="HWLHAPDF";
213       break;
214     case kCTEQ6l:
215       fModPDF=10041;
216       fAutPDF="HWLHAPDF";
217       break;
218     case kCTEQ6ll:
219       fModPDF=10042;
220       fAutPDF="HWLHAPDF";
221       break;
222     case kCTEQ5M:
223       fModPDF=19050;
224       fAutPDF="HWLHAPDF";
225       break;
226     case kCTEQ5L:
227       fModPDF=19070;
228       fAutPDF="HWLHAPDF";
229       break;
230     case kCTEQ4M:
231       fModPDF=19150;
232       fAutPDF="HWLHAPDF";
233       break;
234     case kCTEQ4L:
235       fModPDF=19170;
236       fAutPDF="HWLHAPDF";
237       break;
238 //    case kMRST2004nlo:
239 //      fModPDF=20400;
240 //      fAutPDF="HWLHAPDF";
241 //      break;
242     default:
243       cerr << "This structure function is not inplemented " << fStrucFunc << endl;
244       break;
245     }
246   fAutPDF.Resize(20);
247   fHerwig->SetMODPDF(1,fModPDF);
248   fHerwig->SetMODPDF(2,fModPDF);
249   fHerwig->SetAUTPDF(1,fAutPDF);
250   fHerwig->SetAUTPDF(2,fAutPDF);
251 }
252
253 void AliGenHerwig::Generate()
254 {
255   // Generate one event
256
257   Float_t polar[3] =   {0,0,0};
258   Float_t origin[3]=   {0,0,0};
259   Float_t origin0[3]=  {0,0,0};
260   Float_t p[4], random[6];
261
262   static TClonesArray *particles;
263   //  converts from mm/c to s
264   const Float_t kconv=0.001/2.999792458e8;
265   //
266   Int_t nt=0;
267   Int_t jev=0;
268   Int_t j, kf, ks, imo;
269   kf=0;
270
271   if(!particles) particles=new TClonesArray("TParticle",10000);
272
273   fTrials=0;
274   for (j=0;j<3;j++) origin0[j]=fOrigin[j];
275   if(fVertexSmear==kPerEvent) {
276     Rndm(random,6);
277     for (j=0;j<3;j++) {
278       origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
279         TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
280     }
281   }
282
283   while(1)
284     {
285         fHerwig->GenerateEvent();
286         fTrials++;
287         fHerwig->ImportParticles(particles,"All");
288         Int_t np = particles->GetEntriesFast()-1;
289         if (np == 0 ) continue;
290
291         Int_t nc=0;
292
293         Int_t * newPos = new Int_t[np];
294         for (Int_t i = 0; i<np; i++) *(newPos+i)=-1;
295
296         for (Int_t i = 0; i<np; i++) {
297             TParticle *  iparticle       = (TParticle *) particles->At(i);
298             imo = iparticle->GetFirstMother();
299             kf        = iparticle->GetPdgCode();
300             ks        = iparticle->GetStatusCode();
301             if (ks != 3 &&
302                 KinematicSelection(iparticle,0))
303             {
304                 nc++;
305                 p[0]=iparticle->Px();
306                 p[1]=iparticle->Py();
307                 p[2]=iparticle->Pz();
308                 p[3]=iparticle->Energy();
309                 origin[0]=origin0[0]+iparticle->Vx()/10;
310                 origin[1]=origin0[1]+iparticle->Vy()/10;
311                 origin[2]=origin0[2]+iparticle->Vz()/10;
312                 Float_t tof = kconv*iparticle->T();
313                 Int_t   iparent = (imo > -1) ? newPos[imo] : -1;
314                 Int_t   trackIt = (ks == 1) && fTrackIt;
315                 PushTrack(trackIt, iparent, kf,
316                           p[0], p[1], p[2], p[3],
317                           origin[0], origin[1], origin[2],
318                           tof,
319                           polar[0], polar[1], polar[2],
320                           kPPrimary, nt, fHerwig->GetEVWGT(), ks);
321                 KeepTrack(nt);
322                 newPos[i]=nt;
323             } // end of if: selection of particle
324         } // end of for: particle loop
325         if (newPos) delete[] newPos;
326         //      MakeHeader();
327         if (nc > 0) {
328             jev+=nc;
329             if (jev >= fNpart || fNpart == -1) {
330                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
331                 break;
332             }
333         }
334     }
335   SetHighWaterMark(nt);
336 //  adjust weight due to kinematic selection
337   AdjustWeights();
338 //  get cross-section
339   fXsection=fHerwig->GetAVWGT();
340 }
341
342 void AliGenHerwig::AdjustWeights()
343 {
344 // Adjust the weights after generation of all events
345     TParticle *part;
346     Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
347     for (Int_t i=0; i<ntrack; i++) {
348         part= gAlice->GetMCApp()->Particle(i);
349         part->SetWeight(part->GetWeight()*fKineBias);
350     }
351 }
352
353
354 void AliGenHerwig::KeepFullEvent()
355 {
356     fKeep=1;
357 }
358
359 Bool_t AliGenHerwig::DaughtersSelection(TParticle* iparticle, TClonesArray* particles)
360 {
361 //
362 // Looks recursively if one of the daughters has been selected
363 //
364 //    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
365     Int_t imin=-1;
366     Int_t imax=-1;
367     Int_t i;
368     Bool_t hasDaughters= (iparticle->GetFirstDaughter() >=0);
369     Bool_t selected=kFALSE;
370     if (hasDaughters) {
371         imin=iparticle->GetFirstDaughter();
372         imax=iparticle->GetLastDaughter();
373         for (i=imin; i<= imax; i++){
374             TParticle *  jparticle       = (TParticle *) particles->At(i);
375             Int_t ip=jparticle->GetPdgCode();
376             if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
377                 selected=kTRUE; break;
378             }
379             if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
380         }
381     } else {
382         return kFALSE;
383     }
384
385     return selected;
386 }
387
388
389 Bool_t AliGenHerwig::SelectFlavor(Int_t pid)
390 {
391 // Select flavor of particle
392 // 0: all
393 // 4: charm and beauty
394 // 5: beauty
395     if (fFlavor == 0) return kTRUE;
396
397     Int_t ifl=TMath::Abs(pid/100);
398     if (ifl > 10) ifl/=10;
399     return (fFlavor == ifl);
400 }
401
402 Bool_t AliGenHerwig::Stable(TParticle*  particle)
403 {
404 // Return true for a stable particle
405 //
406     Int_t kf = TMath::Abs(particle->GetPdgCode());
407
408     if ( (particle->GetFirstDaughter() < 0 ) || (kf == 1000*fFlavor+122))
409
410     {
411         return kTRUE;
412     } else {
413         return kFALSE;
414     }
415 }
416
417 void AliGenHerwig::FinishRun()
418 {
419   fHerwig->Hwefin();
420 }
421
422 void AliGenHerwig::FinishRunJimmy()
423 {
424   fHerwig->Hwefin();
425   fHerwig->Jmefin();
426
427 }
428