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