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