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