event printing removed
[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 "AliGenHerwigEventHeader.h"
32 #include "AliHerwigRndm.h"
33 #include "AliMC.h"
34 #include "AliRun.h"
35 #include "driver.h"
36
37 using std::cerr;
38 using std::endl;
39
40 ClassImp(AliGenHerwig)
41
42
43   AliGenHerwig::AliGenHerwig() :
44     AliGenMC(),
45     fAutPDF("LHAPDF"),
46     fModPDF(19070),
47     fStrucFunc(kCTEQ5L),
48     fKeep(0),
49     fDecaysOff(1),
50     fTrigger(0),
51     fSelectAll(0),
52     fFlavor(0),
53     fMomentum1(7000),
54     fMomentum2(7000),
55     fKineBias(1),
56     fTrials(0),
57     fXsection(0),
58     fHerwig(0x0),
59     fProcess(0),
60     fPtHardMin(0.),
61     fPtHardMax(9999.),
62     fPtRMS(0.),
63     fMaxPr(10),
64     fMaxErrors(1000),
65     fEnSoft(1),
66     fEv1Pr(0),
67     fEv2Pr(0),
68     fFileName(0),
69     fEtaMinParton(-20.),     
70     fEtaMaxParton(20.),     
71     fPhiMinParton(0.),     
72     fPhiMaxParton(2.* TMath::Pi()),     
73     fEtaMinGamma(-20.),      
74     fEtaMaxGamma(20.),      
75     fPhiMinGamma(0.),      
76     fPhiMaxGamma(2. * TMath::Pi()),
77     fHeader(0)
78 {
79 // Constructor
80   fEnergyCMS = 14000;
81 }
82
83 AliGenHerwig::AliGenHerwig(Int_t npart)
84     :AliGenMC(npart),
85     fAutPDF("LHAPDF"),
86     fModPDF(19070),
87     fStrucFunc(kCTEQ5L),
88     fKeep(0),
89     fDecaysOff(1),
90     fTrigger(0),
91     fSelectAll(0),
92     fFlavor(0),
93     fMomentum1(7000),
94     fMomentum2(7000),
95     fKineBias(1),
96     fTrials(0),
97     fXsection(0),
98     fHerwig(0x0),
99     fProcess(0),
100     fPtHardMin(0.),
101     fPtHardMax(9999.),
102     fPtRMS(0.),
103     fMaxPr(10),
104     fMaxErrors(1000),
105     fEnSoft(1),
106     fEv1Pr(0),
107     fEv2Pr(0),
108     fFileName(0),
109     fEtaMinParton(-20.),     
110     fEtaMaxParton(20.),     
111     fPhiMinParton(0.),     
112     fPhiMaxParton(2.* TMath::Pi()),     
113     fEtaMinGamma(-20.),      
114     fEtaMaxGamma(20.),      
115     fPhiMinGamma(0.),
116     fPhiMaxGamma(2. * TMath::Pi()),
117     fHeader(0)
118 {
119 // Constructor
120     fEnergyCMS = 14000;
121     SetTarget();
122     SetProjectile();
123     // Set random number generator
124     AliHerwigRndm::SetHerwigRandom(GetRandom());
125 }
126
127 AliGenHerwig::~AliGenHerwig()
128 {
129 // Destructor
130 }
131
132 void AliGenHerwig::SetEventListRange(Int_t eventFirst, Int_t eventLast)
133 {
134   fEv1Pr = eventFirst;
135   fEv2Pr = eventLast;
136   if ( fEv2Pr == -1 ) fEv2Pr = fEv1Pr;
137 }
138
139 void AliGenHerwig::Init()
140 {
141 // Initialisation
142   fTarget.Resize(8);
143   fProjectile.Resize(8);
144   SetMC(new THerwig6());
145   fHerwig=(THerwig6*) fMCEvGen;
146   // initialize common blocks
147   fHerwig->Initialize(fProjectile.Data(), fTarget.Data(), fMomentum1, fMomentum2, fProcess); // 
148   // reset parameters according to user needs
149   InitPDF();
150   fHerwig->SetPTMIN(fPtHardMin);
151   fHerwig->SetPTMAX(fPtHardMax);
152   fHerwig->SetPTRMS(fPtRMS);
153   fHerwig->SetMAXPR(fMaxPr);
154   fHerwig->SetMAXER(fMaxErrors);
155   fHerwig->SetENSOF(fEnSoft);
156 // C---D,U,S,C,B,T QUARK AND GLUON MASSES (IN THAT ORDER)
157 //       RMASS(1)=0.32
158 //       RMASS(2)=0.32
159 //       RMASS(3)=0.5
160 //       RMASS(4)=1.55
161 //       RMASS(5)=4.75
162 //       RMASS(6)=174.3
163 //       RMASS(13)=0.75
164
165   fHerwig->SetRMASS(4,1.2);
166   fHerwig->SetRMASS(5,4.75);
167
168   if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(), 49);
169
170   //fHerwig->Hwusta("PI0     ");
171
172   // compute parameter dependent constants
173   fHerwig->PrepareRun();
174 }
175
176 void AliGenHerwig::InitJimmy()
177 {
178 // Initialisation
179   fTarget.Resize(8);
180   fProjectile.Resize(8);
181   SetMC(new THerwig6());
182   fHerwig=(THerwig6*) fMCEvGen;
183   // initialize common blocks
184   fHerwig->InitializeJimmy(fProjectile, fTarget, fMomentum1, fMomentum2, fProcess);
185   // reset parameters according to user needs
186   InitPDF();
187   fHerwig->SetPTMIN(fPtHardMin);
188   fHerwig->SetPTRMS(fPtRMS);
189   fHerwig->SetMAXPR(fMaxPr);
190   fHerwig->SetMAXER(fMaxErrors);
191   fHerwig->SetENSOF(fEnSoft);
192
193 // C---D,U,S,C,B,T QUARK AND GLUON MASSES (IN THAT ORDER)
194 //       RMASS(1)=0.32
195 //       RMASS(2)=0.32
196 //       RMASS(3)=0.5
197 //       RMASS(4)=1.55
198 //       RMASS(5)=4.75
199 //       RMASS(6)=174.3
200 //       RMASS(13)=0.75
201
202   fHerwig->SetRMASS(4,1.2);
203   fHerwig->SetRMASS(5,4.75);
204
205   if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(), 49);
206
207   //  fHerwig->Hwusta("PI0     ");
208
209   // compute parameter dependent constants
210   fHerwig->PrepareRunJimmy();
211 }
212
213 void AliGenHerwig::InitPDF()
214 {
215 // Initialize PDF
216   switch(fStrucFunc)
217     {
218 // ONLY USES LHAPDF STRUCTURE FUNCTIONS
219     case kGRVLO98:
220       fModPDF=80060;
221       fAutPDF="HWLHAPDF";
222       break;
223     case kCTEQ6:
224       fModPDF=10040;
225       fAutPDF="HWLHAPDF";
226       break;
227     case kCTEQ61:
228       fModPDF=10100;
229       fAutPDF="HWLHAPDF";
230       break;
231     case kCTEQ6m:
232       fModPDF=10050;
233       fAutPDF="HWLHAPDF";
234       break;
235     case kCTEQ6l:
236       fModPDF=10041;
237       fAutPDF="HWLHAPDF";
238       break;
239     case kCTEQ6ll:
240       fModPDF=10042;
241       fAutPDF="HWLHAPDF";
242       break;
243     case kCTEQ5M:
244       fModPDF=19050;
245       fAutPDF="HWLHAPDF";
246       break;
247     case kCTEQ5L:
248       fModPDF=19070;
249       fAutPDF="HWLHAPDF";
250       break;
251     case kCTEQ4M:
252       fModPDF=19150;
253       fAutPDF="HWLHAPDF";
254       break;
255     case kCTEQ4L:
256       fModPDF=19170;
257       fAutPDF="HWLHAPDF";
258       break;
259 //    case kMRST2004nlo:
260 //      fModPDF=20400;
261 //      fAutPDF="HWLHAPDF";
262 //      break;
263     default:
264       cerr << "This structure function is not inplemented " << fStrucFunc << endl;
265       break;
266     }
267   fAutPDF.Resize(20);
268   fHerwig->SetMODPDF(1,fModPDF);
269   fHerwig->SetMODPDF(2,fModPDF);
270   fHerwig->SetAUTPDF(1,fAutPDF);
271   fHerwig->SetAUTPDF(2,fAutPDF);
272 }
273
274 void AliGenHerwig::Generate()
275 {
276   // Generate one event
277
278   Float_t polar[3]  =   {0,0,0};
279   Float_t origin[3] =   {0,0,0};
280   Float_t p[4];
281
282   static TClonesArray *particles;
283   //  converts from mm/c to s
284   const Float_t kconv=0.001/2.999792458e8;
285   //
286   Int_t nt=0;
287   Int_t jev=0;
288   Int_t kf, ks, imo;
289   kf=0;
290
291   if(!particles) particles=new TClonesArray("TParticle",10000);
292
293   fTrials=0;
294
295   //  Set collision vertex position 
296   if (fVertexSmear == kPerEvent) Vertex();
297
298   while(1)
299     {
300         fHerwig->GenerateEvent();
301         fTrials++;
302         fHerwig->ImportParticles(particles,"All");
303         Int_t np = particles->GetEntriesFast()-1;
304         if (np == 0 ) continue;
305
306         //Check hard partons or direct gamma in kine range
307
308         if (fProcess == kHeJets || fProcess == kHeDirectGamma) {
309             TParticle* parton1 = (TParticle *) particles->At(6);
310             TParticle* parton2 = (TParticle *) particles->At(7);
311             if (!CheckParton(parton1, parton2))  continue ;
312         } 
313
314         // 
315         if (gAlice) {
316             if (gAlice->GetEvNumber()>=fEv1Pr &&
317                 gAlice->GetEvNumber()<=fEv2Pr) fHerwig->PrintEvt();
318
319         }
320
321         
322         Int_t nc = 0;
323         fNprimaries = 0;
324         
325         Int_t * newPos = new Int_t[np];
326         for (Int_t i = 0; i<np; i++) *(newPos+i)=-1;
327
328         for (Int_t i = 0; i<np; i++) {
329             TParticle *  iparticle       = (TParticle *) particles->At(i);
330             imo = iparticle->GetFirstMother();
331             kf        = iparticle->GetPdgCode();
332             ks        = iparticle->GetStatusCode();
333             if (ks != 3 &&
334                 KinematicSelection(iparticle,0))
335             {
336                 nc++;
337                 p[0]=iparticle->Px();
338                 p[1]=iparticle->Py();
339                 p[2]=iparticle->Pz();
340                 p[3]=iparticle->Energy();
341
342                 origin[0] = fVertex[0] + iparticle->Vx()/10; // [cm]
343                 origin[1] = fVertex[1] + iparticle->Vy()/10; // [cm]
344                 origin[2] = fVertex[2] + iparticle->Vz()/10; // [cm]
345
346                 Float_t tof = fTime + kconv*iparticle->T();
347                 Int_t   iparent = (imo > -1) ? newPos[imo] : -1;
348                 Int_t   trackIt = (ks == 1) && fTrackIt;
349                 PushTrack(trackIt, iparent, kf,
350                           p[0], p[1], p[2], p[3],
351                           origin[0], origin[1], origin[2],
352                           tof,
353                           polar[0], polar[1], polar[2],
354                           kPPrimary, nt, fHerwig->GetEVWGT(), ks);
355                 KeepTrack(nt);
356                 newPos[i]=nt;
357                 fNprimaries++;
358             } // end of if: selection of particle
359         } // end of for: particle loop
360         if (newPos) delete[] newPos;
361         if (nc > 0) {
362             jev+=nc;
363             if (jev >= fNpart || fNpart == -1) {
364                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
365                 break;
366             }
367         }
368     }
369 //
370   MakeHeader();
371 //  
372   SetHighWaterMark(nt);
373 //  adjust weight due to kinematic selection
374   AdjustWeights();
375 //  get cross-section
376   fXsection=fHerwig->GetAVWGT();
377   //printf(">> trials << %d\n",fTrials);
378 }
379
380 Bool_t AliGenHerwig::CheckParton(const TParticle* parton1, const TParticle* parton2)
381 {
382 // Check the kinematic trigger condition
383 //
384 //Select events with parton max energy
385     if(fPtHardMax < parton1->Pt()) return kFALSE;
386
387 // Select events within angular window
388     Double_t eta[2];
389     eta[0] = parton1->Eta();
390     eta[1] = parton2->Eta();
391     Double_t phi[2];
392     phi[0] = parton1->Phi();
393     phi[1] = parton2->Phi();
394     Int_t    pdg[2]; 
395     pdg[0] = parton1->GetPdgCode();
396     pdg[1] = parton2->GetPdgCode();   
397     printf("min %f, max %f\n",fPtHardMin, fPtHardMax);
398     printf("Parton 1: %s, pT= %2.2f, eta = %1.2f, phi = %2.2f\n", parton1->GetName(),parton1->Pt(), eta[0],phi[0]*TMath::RadToDeg());
399     printf("Parton 2: %s, pT= %2.2f, eta = %1.2f, phi = %2.2f\n", parton2->GetName(),parton2->Pt(), eta[1],phi[1]*TMath::RadToDeg());
400     
401     if (fProcess == kHeJets) {
402       //Check if one of the 2 outgoing partons are in the eta-phi window
403       for(Int_t i = 0; i < 2; i++)
404         if ((eta[i] < fEtaMaxParton  && eta[i] > fEtaMinParton) &&
405             (phi[i] < fPhiMaxParton  && phi[i] > fPhiMinParton)) return  kTRUE ;
406     }
407     
408     else {
409       //Check if the gamma and the jet  are in the eta-phi window
410       Int_t igj = 0;
411       Int_t ijj = 0;
412       if(pdg[0] == 22) ijj=1;
413       else igj=1;
414       if ((eta[ijj] < fEtaMaxParton   && eta[ijj] > fEtaMinParton) &&
415           (phi[ijj] < fPhiMaxParton   && phi[ijj] > fPhiMinParton)) {
416         
417         if ((eta[igj] < fEtaMaxGamma   && eta[igj] > fEtaMinGamma) &&
418             (phi[igj] < fPhiMaxGamma   && phi[igj] > fPhiMinGamma)) return  kTRUE;
419         
420       }
421     }
422
423     return kFALSE ;
424 }
425
426 void AliGenHerwig::AdjustWeights()
427 {
428 // Adjust the weights after generation of all events
429     TParticle *part;
430     Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
431     for (Int_t i=0; i<ntrack; i++) {
432         part= gAlice->GetMCApp()->Particle(i);
433         part->SetWeight(part->GetWeight()*fKineBias);
434     }
435 }
436
437
438 void AliGenHerwig::KeepFullEvent()
439 {
440     fKeep=1;
441 }
442
443 Bool_t AliGenHerwig::DaughtersSelection(const TParticle* iparticle, const TClonesArray* particles)
444 {
445 //
446 // Looks recursively if one of the daughters has been selected
447 //
448 //    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
449     Int_t imin=-1;
450     Int_t imax=-1;
451     Int_t i;
452     Bool_t hasDaughters= (iparticle->GetFirstDaughter() >=0);
453     Bool_t selected=kFALSE;
454     if (hasDaughters) {
455         imin=iparticle->GetFirstDaughter();
456         imax=iparticle->GetLastDaughter();
457         for (i=imin; i<= imax; i++){
458             TParticle *  jparticle       = (TParticle *) particles->At(i);
459             Int_t ip=jparticle->GetPdgCode();
460             if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
461                 selected=kTRUE; break;
462             }
463             if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
464         }
465     } else {
466         return kFALSE;
467     }
468
469     return selected;
470 }
471
472
473 Bool_t AliGenHerwig::SelectFlavor(Int_t pid) const
474 {
475 // Select flavor of particle
476 // 0: all
477 // 4: charm and beauty
478 // 5: beauty
479     if (fFlavor == 0) return kTRUE;
480
481     Int_t ifl=TMath::Abs(pid/100);
482     if (ifl > 10) ifl/=10;
483     return (fFlavor == ifl);
484 }
485
486 Bool_t AliGenHerwig::Stable(const TParticle*  particle) const
487 {
488 // Return true for a stable particle
489 //
490     Int_t kf = TMath::Abs(particle->GetPdgCode());
491
492     if ( (particle->GetFirstDaughter() < 0 ) || (kf == 1000*fFlavor+122))
493
494     {
495         return kTRUE;
496     } else {
497         return kFALSE;
498     }
499 }
500
501 void AliGenHerwig::FinishRun()
502 {
503   fHerwig->Hwefin();
504 }
505
506 void AliGenHerwig::FinishRunJimmy()
507 {
508   fHerwig->Hwefin();
509   fHerwig->Jmefin();
510
511 }
512
513
514 void AliGenHerwig::MakeHeader()
515 {
516 //
517 // Make header for the simulated event
518 // 
519   if (fHeader) delete fHeader;
520   fHeader = new AliGenHerwigEventHeader("Herwig");
521 //
522 // Event type  
523     ((AliGenHerwigEventHeader*) fHeader)->SetProcessType(fHerwig->GetIHPRO());
524 //
525 // Number of trials
526     ((AliGenHerwigEventHeader*) fHeader)->SetTrials(fTrials);
527 //
528 // Event weight (cross section)
529     ((AliGenHerwigEventHeader*) fHeader)->SetWeight(fHerwig->GetEVWGT());
530 //
531 // Event Vertex 
532     fHeader->SetPrimaryVertex(fVertex);
533     fHeader->SetInteractionTime(fTime);
534 //
535 // Number of primaries
536     fHeader->SetNProduced(fNprimaries);
537 //  Pass header
538 //
539     AddHeader(fHeader);
540     fHeader = 0x0;  
541 }