effc++ changes
[u/mrichter/AliRoot.git] / PYTHIA6 / AliGenPythia.cxx
CommitLineData
8d2cd130 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
7cdba479 16/* $Id$ */
8d2cd130 17
18//
19// Generator using the TPythia interface (via AliPythia)
20// to generate pp collisions.
21// Using SetNuclei() also nuclear modifications to the structure functions
22// can be taken into account. This makes, of course, only sense for the
23// generation of the products of hard processes (heavy flavor, jets ...)
24//
25// andreas.morsch@cern.ch
26//
27
28#include <TDatabasePDG.h>
29#include <TParticle.h>
30#include <TPDGCode.h>
31#include <TSystem.h>
32#include <TTree.h>
8d2cd130 33#include "AliConst.h"
34#include "AliDecayerPythia.h"
35#include "AliGenPythia.h"
5fa4b20b 36#include "AliHeader.h"
8d2cd130 37#include "AliGenPythiaEventHeader.h"
38#include "AliPythia.h"
7cdba479 39#include "AliPythiaRndm.h"
8d2cd130 40#include "AliRun.h"
7ea3ea5b 41#include "AliStack.h"
42#include "AliRunLoader.h"
5d12ce38 43#include "AliMC.h"
7c21f297 44#include "pyquenCommon.h"
8d2cd130 45
014a9521 46ClassImp(AliGenPythia)
8d2cd130 47
48AliGenPythia::AliGenPythia()
49 :AliGenMC()
50{
51// Default Constructor
52 fParticles = 0;
53 fPythia = 0;
e5c87a3d 54 fHeader = 0;
5fffd0fe 55 fReadFromFile = 0;
f913ec4f 56 fEventTime = 0.;
8ec89232 57 fInteractionRate = 0.;
9d7108a7 58 fTimeWindow = 0.;
59 fEventsTime = 0;
60 fCurSubEvent = 0;
8d2cd130 61 fDecayer = new AliDecayerPythia();
62 SetEventListRange();
63 SetJetPhiRange();
64 SetJetEtaRange();
65 SetJetEtRange();
66 SetGammaPhiRange();
67 SetGammaEtaRange();
68 SetPtKick();
7ea3ea5b 69 SetQuench();
5fa4b20b 70 SetHadronisation();
e0e89f40 71 SetTriggerParticle();
2d9f9817 72 SetNuclei(0,0);
1d568bc2 73 fSetNuclei = kFALSE;
beac474c 74 fNewMIS = kFALSE;
75 fHFoff = kFALSE;
76d6ba9a 76 fGinit = 1;
77 fGfinal = 1;
78
7cdba479 79 if (!AliPythiaRndm::GetPythiaRandom())
80 AliPythiaRndm::SetPythiaRandom(GetRandom());
8d2cd130 81}
82
83AliGenPythia::AliGenPythia(Int_t npart)
84 :AliGenMC(npart)
85{
86// default charm production at 5. 5 TeV
87// semimuonic decay
88// structure function GRVHO
89//
90 fName = "Pythia";
91 fTitle= "Particle Generator using PYTHIA";
92 fXsection = 0.;
5fffd0fe 93 fReadFromFile = 0;
f913ec4f 94 fEventTime = 0.;
8ec89232 95 fInteractionRate = 0.;
9d7108a7 96 fTimeWindow = 0.;
97 fEventsTime = 0;
98 fCurSubEvent = 0;
8d2cd130 99 SetProcess();
100 SetStrucFunc();
101 SetForceDecay();
102 SetPtHard();
103 SetYHard();
104 SetEnergyCMS();
105 fDecayer = new AliDecayerPythia();
106 // Set random number generator
7cdba479 107 if (!AliPythiaRndm::GetPythiaRandom())
108 AliPythiaRndm::SetPythiaRandom(GetRandom());
8d2cd130 109 fFlavorSelect = 0;
110 // Produced particles
111 fParticles = new TClonesArray("TParticle",1000);
e5c87a3d 112 fHeader = 0;
8d2cd130 113 SetEventListRange();
114 SetJetPhiRange();
115 SetJetEtaRange();
116 SetJetEtRange();
117 SetGammaPhiRange();
118 SetGammaEtaRange();
119 SetJetReconstructionMode();
7ea3ea5b 120 SetQuench();
5fa4b20b 121 SetHadronisation();
8d2cd130 122 SetPtKick();
e0e89f40 123 SetTriggerParticle();
2d9f9817 124 SetNuclei(0,0);
8d2cd130 125 // Options determining what to keep in the stack (Heavy flavour generation)
126 fStackFillOpt = kFlavorSelection; // Keep particle with selected flavor
127 fFeedDownOpt = kTRUE; // allow feed down from higher family
128 // Fragmentation on/off
129 fFragmentation = kTRUE;
130 // Default counting mode
131 fCountMode = kCountAll;
592f8307 132 // Pycel
133 SetPycellParameters();
1d568bc2 134 fSetNuclei = kFALSE;
beac474c 135 fNewMIS = kFALSE;
136 fHFoff = kFALSE;
76d6ba9a 137 fGinit = 1;
138 fGfinal = 1;
139 }
8d2cd130 140
141AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
f936b8ea 142 :AliGenMC(Pythia)
8d2cd130 143{
144// copy constructor
145 Pythia.Copy(*this);
146}
147
148AliGenPythia::~AliGenPythia()
149{
150// Destructor
9d7108a7 151 if(fEventsTime) delete fEventsTime;
152}
153
154void AliGenPythia::SetInteractionRate(Float_t rate,Float_t timewindow)
155{
156// Generate pileup using user specified rate
157 fInteractionRate = rate;
158 fTimeWindow = timewindow;
159 GeneratePileup();
160}
161
162void AliGenPythia::GeneratePileup()
163{
164// Generate sub events time for pileup
165 fEventsTime = 0;
166 if(fInteractionRate == 0.) {
167 Warning("GeneratePileup","Zero interaction specified. Skipping pileup generation.\n");
168 return;
169 }
170
171 Int_t npart = NumberParticles();
172 if(npart < 0) {
173 Warning("GeneratePileup","Negative number of particles. Skipping pileup generation.\n");
174 return;
175 }
176
177 if(fEventsTime) delete fEventsTime;
178 fEventsTime = new TArrayF(npart);
179 TArrayF &array = *fEventsTime;
180 for(Int_t ipart = 0; ipart < npart; ipart++)
181 array[ipart] = 0.;
182
183 Float_t eventtime = 0.;
184 while(1)
185 {
186 eventtime += (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
187 if(eventtime > fTimeWindow) break;
188 array.Set(array.GetSize()+1);
189 array[array.GetSize()-1] = eventtime;
190 }
191
192 eventtime = 0.;
193 while(1)
194 {
195 eventtime -= (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
196 if(TMath::Abs(eventtime) > fTimeWindow) break;
197 array.Set(array.GetSize()+1);
198 array[array.GetSize()-1] = eventtime;
199 }
200
201 SetNumberParticles(fEventsTime->GetSize());
8d2cd130 202}
203
592f8307 204void AliGenPythia::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
205 Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
206{
207// Set pycell parameters
208 fPycellEtaMax = etamax;
209 fPycellNEta = neta;
210 fPycellNPhi = nphi;
211 fPycellThreshold = thresh;
212 fPycellEtSeed = etseed;
213 fPycellMinEtJet = minet;
214 fPycellMaxRadius = r;
215}
216
217
218
8d2cd130 219void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
220{
221 // Set a range of event numbers, for which a table
222 // of generated particle will be printed
223 fDebugEventFirst = eventFirst;
224 fDebugEventLast = eventLast;
225 if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
226}
227
228void AliGenPythia::Init()
229{
230// Initialisation
231
232 SetMC(AliPythia::Instance());
b88f5cea 233 fPythia=(AliPythia*) fMCEvGen;
e2bddf81 234
8d2cd130 235//
236 fParentWeight=1./Float_t(fNpart);
237//
238// Forward Paramters to the AliPythia object
239 fDecayer->SetForceDecay(fForceDecay);
240 fDecayer->Init();
241
242
243 fPythia->SetCKIN(3,fPtHardMin);
244 fPythia->SetCKIN(4,fPtHardMax);
245 fPythia->SetCKIN(7,fYHardMin);
246 fPythia->SetCKIN(8,fYHardMax);
247
1a626d4e 248 if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget);
8d2cd130 249 // Fragmentation?
250 if (fFragmentation) {
251 fPythia->SetMSTP(111,1);
252 } else {
253 fPythia->SetMSTP(111,0);
254 }
255
256
257// initial state radiation
258 fPythia->SetMSTP(61,fGinit);
259// final state radiation
260 fPythia->SetMSTP(71,fGfinal);
261// pt - kick
262 if (fPtKick > 0.) {
263 fPythia->SetMSTP(91,1);
264 fPythia->SetPARP(91,fPtKick);
265 } else {
266 fPythia->SetMSTP(91,0);
267 }
268
5fa4b20b 269
270 if (fReadFromFile) {
271 fRL = AliRunLoader::Open(fFileName, "Partons");
272 fRL->LoadKinematics();
273 fRL->LoadHeader();
274 } else {
275 fRL = 0x0;
276 }
beac474c 277// Switch off Heavy Flavors on request
278 if (fHFoff) {
279 fPythia->SetMSTP(58, 3);
280 fPythia->SetMSTJ(45, 3);
281 for (Int_t i = 156; i <= 160; i++) fPythia->SetMDME(i, 1, 0);
282 }
8d2cd130 283 //
284 fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
285
286// Parent and Children Selection
287 switch (fProcess)
288 {
65f2626c 289 case kPyOldUEQ2ordered:
290 case kPyOldUEQ2ordered2:
291 case kPyOldPopcorn:
292 break;
8d2cd130 293 case kPyCharm:
294 case kPyCharmUnforced:
adf4d898 295 case kPyCharmPbPbMNR:
aabc7187 296 case kPyCharmpPbMNR:
e0e89f40 297 case kPyCharmppMNR:
298 case kPyCharmppMNRwmi:
8d2cd130 299 fParentSelect[0] = 411;
300 fParentSelect[1] = 421;
301 fParentSelect[2] = 431;
302 fParentSelect[3] = 4122;
303 fFlavorSelect = 4;
304 break;
adf4d898 305 case kPyD0PbPbMNR:
306 case kPyD0pPbMNR:
307 case kPyD0ppMNR:
8d2cd130 308 fParentSelect[0] = 421;
309 fFlavorSelect = 4;
310 break;
90d7b703 311 case kPyDPlusPbPbMNR:
312 case kPyDPluspPbMNR:
313 case kPyDPlusppMNR:
314 fParentSelect[0] = 411;
315 fFlavorSelect = 4;
316 break;
e0e89f40 317 case kPyDPlusStrangePbPbMNR:
318 case kPyDPlusStrangepPbMNR:
319 case kPyDPlusStrangeppMNR:
320 fParentSelect[0] = 431;
321 fFlavorSelect = 4;
322 break;
8d2cd130 323 case kPyBeauty:
adf4d898 324 case kPyBeautyPbPbMNR:
325 case kPyBeautypPbMNR:
326 case kPyBeautyppMNR:
e0e89f40 327 case kPyBeautyppMNRwmi:
8d2cd130 328 fParentSelect[0]= 511;
329 fParentSelect[1]= 521;
330 fParentSelect[2]= 531;
331 fParentSelect[3]= 5122;
332 fParentSelect[4]= 5132;
333 fParentSelect[5]= 5232;
334 fParentSelect[6]= 5332;
335 fFlavorSelect = 5;
336 break;
337 case kPyBeautyUnforced:
338 fParentSelect[0] = 511;
339 fParentSelect[1] = 521;
340 fParentSelect[2] = 531;
341 fParentSelect[3] = 5122;
342 fParentSelect[4] = 5132;
343 fParentSelect[5] = 5232;
344 fParentSelect[6] = 5332;
345 fFlavorSelect = 5;
346 break;
347 case kPyJpsiChi:
348 case kPyJpsi:
349 fParentSelect[0] = 443;
350 break;
351 case kPyMb:
352 case kPyMbNonDiffr:
d7de4a67 353 case kPyMbMSEL1:
8d2cd130 354 case kPyJets:
355 case kPyDirectGamma:
356 break;
589380c6 357 case kPyW:
0f6ee828 358 case kPyZ:
589380c6 359 break;
8d2cd130 360 }
361//
592f8307 362//
363// JetFinder for Trigger
364//
365// Configure detector (EMCAL like)
366//
d7de4a67 367 fPythia->SetPARU(51, fPycellEtaMax);
368 fPythia->SetMSTU(51, fPycellNEta);
369 fPythia->SetMSTU(52, fPycellNPhi);
592f8307 370//
371// Configure Jet Finder
372//
d7de4a67 373 fPythia->SetPARU(58, fPycellThreshold);
374 fPythia->SetPARU(52, fPycellEtSeed);
375 fPythia->SetPARU(53, fPycellMinEtJet);
376 fPythia->SetPARU(54, fPycellMaxRadius);
377 fPythia->SetMSTU(54, 2);
592f8307 378//
8d2cd130 379// This counts the total number of calls to Pyevnt() per run.
380 fTrialsRun = 0;
381 fQ = 0.;
382 fX1 = 0.;
383 fX2 = 0.;
384 fNev = 0 ;
385//
1d568bc2 386//
387//
8d2cd130 388 AliGenMC::Init();
1d568bc2 389//
390//
391//
392 if (fSetNuclei) {
393 fDyBoost = 0;
394 Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
395 }
d7de4a67 396
32d6ef7d 397 if (fQuench) {
5fa4b20b 398 fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
32d6ef7d 399 }
8d2cd130 400}
401
402void AliGenPythia::Generate()
403{
404// Generate one event
e5c87a3d 405
8d2cd130 406 fDecayer->ForceDecay();
407
408 Float_t polar[3] = {0,0,0};
409 Float_t origin[3] = {0,0,0};
a920faf9 410 Float_t p[4];
8d2cd130 411// converts from mm/c to s
412 const Float_t kconv=0.001/2.999792458e8;
413//
414 Int_t nt=0;
415 Int_t jev=0;
416 Int_t j, kf;
417 fTrials=0;
f913ec4f 418 fEventTime = 0.;
419
2590ccf9 420
8d2cd130 421
422 // Set collision vertex position
2590ccf9 423 if (fVertexSmear == kPerEvent) Vertex();
424
8d2cd130 425// event loop
426 while(1)
427 {
32d6ef7d 428//
5fa4b20b 429// Produce event
32d6ef7d 430//
32d6ef7d 431//
5fa4b20b 432// Switch hadronisation off
433//
434 fPythia->SetMSTJ(1, 0);
32d6ef7d 435//
5fa4b20b 436// Either produce new event or read partons from file
437//
438 if (!fReadFromFile) {
beac474c 439 if (!fNewMIS) {
440 fPythia->Pyevnt();
441 } else {
442 fPythia->Pyevnw();
443 }
5fa4b20b 444 fNpartons = fPythia->GetN();
445 } else {
446 printf("Loading Event %d\n",AliRunLoader::GetRunLoader()->GetEventNumber());
447 fRL->GetEvent(AliRunLoader::GetRunLoader()->GetEventNumber());
448 fPythia->SetN(0);
449 LoadEvent(fRL->Stack(), 0 , 1);
450 fPythia->Pyedit(21);
451 }
452
32d6ef7d 453//
454// Run quenching routine
455//
5fa4b20b 456 if (fQuench == 1) {
457 fPythia->Quench();
458 } else if (fQuench == 2){
459 fPythia->Pyquen(208., 0, 0.);
460 }
32d6ef7d 461//
5fa4b20b 462// Switch hadronisation on
32d6ef7d 463//
aea21c57 464 fPythia->SetMSTJ(1, 1);
5fa4b20b 465//
466// .. and perform hadronisation
aea21c57 467// printf("Calling hadronisation %d\n", fPythia->GetN());
5fa4b20b 468 fPythia->Pyexec();
8d2cd130 469 fTrials++;
8d2cd130 470 fPythia->ImportParticles(fParticles,"All");
1d568bc2 471 Boost();
8d2cd130 472//
473//
474//
475 Int_t i;
476
5fa4b20b 477
8d2cd130 478 Int_t np = fParticles->GetEntriesFast();
5fa4b20b 479
7c21f297 480 if (np == 0) continue;
8d2cd130 481//
2590ccf9 482
8d2cd130 483//
484 Int_t* pParent = new Int_t[np];
485 Int_t* pSelected = new Int_t[np];
486 Int_t* trackIt = new Int_t[np];
5fa4b20b 487 for (i = 0; i < np; i++) {
8d2cd130 488 pParent[i] = -1;
489 pSelected[i] = 0;
490 trackIt[i] = 0;
491 }
492
493 Int_t nc = 0; // Total n. of selected particles
494 Int_t nParents = 0; // Selected parents
495 Int_t nTkbles = 0; // Trackable particles
496 if (fProcess != kPyMb && fProcess != kPyJets &&
497 fProcess != kPyDirectGamma &&
589380c6 498 fProcess != kPyMbNonDiffr &&
d7de4a67 499 fProcess != kPyMbMSEL1 &&
e0e89f40 500 fProcess != kPyW && fProcess != kPyZ &&
501 fProcess != kPyCharmppMNRwmi && fProcess != kPyBeautyppMNRwmi) {
8d2cd130 502
5fa4b20b 503 for (i = 0; i < np; i++) {
2590ccf9 504 TParticle* iparticle = (TParticle *) fParticles->At(i);
8d2cd130 505 Int_t ks = iparticle->GetStatusCode();
506 kf = CheckPDGCode(iparticle->GetPdgCode());
507// No initial state partons
508 if (ks==21) continue;
509//
510// Heavy Flavor Selection
511//
512 // quark ?
513 kf = TMath::Abs(kf);
514 Int_t kfl = kf;
9ff6c04c 515 // Resonance
f913ec4f 516
9ff6c04c 517 if (kfl > 100000) kfl %= 100000;
183a5ca9 518 if (kfl > 10000) kfl %= 10000;
8d2cd130 519 // meson ?
520 if (kfl > 10) kfl/=100;
521 // baryon
522 if (kfl > 10) kfl/=10;
8d2cd130 523 Int_t ipa = iparticle->GetFirstMother()-1;
524 Int_t kfMo = 0;
f913ec4f 525//
526// Establish mother daughter relation between heavy quarks and mesons
527//
528 if (kf >= fFlavorSelect && kf <= 6) {
529 Int_t idau = iparticle->GetFirstDaughter() - 1;
530 if (idau > -1) {
531 TParticle* daughter = (TParticle *) fParticles->At(idau);
532 Int_t pdgD = daughter->GetPdgCode();
533 if (pdgD == 91 || pdgD == 92) {
534 Int_t jmin = daughter->GetFirstDaughter() - 1;
535 Int_t jmax = daughter->GetLastDaughter() - 1;
536 for (Int_t j = jmin; j <= jmax; j++)
537 ((TParticle *) fParticles->At(j))->SetFirstMother(i+1);
538 } // is string or cluster
539 } // has daughter
540 } // heavy quark
8d2cd130 541
f913ec4f 542
8d2cd130 543 if (ipa > -1) {
544 TParticle * mother = (TParticle *) fParticles->At(ipa);
545 kfMo = TMath::Abs(mother->GetPdgCode());
546 }
f913ec4f 547
8d2cd130 548 // What to keep in Stack?
549 Bool_t flavorOK = kFALSE;
550 Bool_t selectOK = kFALSE;
551 if (fFeedDownOpt) {
32d6ef7d 552 if (kfl >= fFlavorSelect) flavorOK = kTRUE;
8d2cd130 553 } else {
32d6ef7d 554 if (kfl > fFlavorSelect) {
555 nc = -1;
556 break;
557 }
558 if (kfl == fFlavorSelect) flavorOK = kTRUE;
8d2cd130 559 }
560 switch (fStackFillOpt) {
561 case kFlavorSelection:
32d6ef7d 562 selectOK = kTRUE;
563 break;
8d2cd130 564 case kParentSelection:
32d6ef7d 565 if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
566 break;
8d2cd130 567 }
568 if (flavorOK && selectOK) {
569//
570// Heavy flavor hadron or quark
571//
572// Kinematic seletion on final state heavy flavor mesons
573 if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
574 {
9ff6c04c 575 continue;
8d2cd130 576 }
577 pSelected[i] = 1;
578 if (ParentSelected(kf)) ++nParents; // Update parent count
579// printf("\n particle (HF) %d %d %d", i, pSelected[i], kf);
580 } else {
581// Kinematic seletion on decay products
582 if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
9ff6c04c 583 && !KinematicSelection(iparticle, 1))
8d2cd130 584 {
9ff6c04c 585 continue;
8d2cd130 586 }
587//
588// Decay products
589// Select if mother was selected and is not tracked
590
591 if (pSelected[ipa] &&
592 !trackIt[ipa] && // mother will be tracked ?
593 kfMo != 5 && // mother is b-quark, don't store fragments
594 kfMo != 4 && // mother is c-quark, don't store fragments
595 kf != 92) // don't store string
596 {
597//
598// Semi-stable or de-selected: diselect decay products:
599//
600//
601 if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > fMaxLifeTime)
602 {
603 Int_t ipF = iparticle->GetFirstDaughter();
604 Int_t ipL = iparticle->GetLastDaughter();
605 if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
606 }
607// printf("\n particle (decay) %d %d %d", i, pSelected[i], kf);
608 pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
609 }
610 }
611 if (pSelected[i] == -1) pSelected[i] = 0;
612 if (!pSelected[i]) continue;
613 // Count quarks only if you did not include fragmentation
614 if (fFragmentation && kf <= 10) continue;
9ff6c04c 615
8d2cd130 616 nc++;
617// Decision on tracking
618 trackIt[i] = 0;
619//
620// Track final state particle
621 if (ks == 1) trackIt[i] = 1;
622// Track semi-stable particles
d25cfd65 623 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime)) trackIt[i] = 1;
8d2cd130 624// Track particles selected by process if undecayed.
625 if (fForceDecay == kNoDecay) {
626 if (ParentSelected(kf)) trackIt[i] = 1;
627 } else {
628 if (ParentSelected(kf)) trackIt[i] = 0;
629 }
630 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
631//
632//
633
634 } // particle selection loop
635 if (nc > 0) {
636 for (i = 0; i<np; i++) {
637 if (!pSelected[i]) continue;
638 TParticle * iparticle = (TParticle *) fParticles->At(i);
639 kf = CheckPDGCode(iparticle->GetPdgCode());
640 Int_t ks = iparticle->GetStatusCode();
641 p[0] = iparticle->Px();
642 p[1] = iparticle->Py();
643 p[2] = iparticle->Pz();
a920faf9 644 p[3] = iparticle->Energy();
645
2590ccf9 646 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
647 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
648 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
649
8d2cd130 650 Float_t tof = kconv*iparticle->T();
651 Int_t ipa = iparticle->GetFirstMother()-1;
652 Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
a920faf9 653
654 PushTrack(fTrackIt*trackIt[i], iparent, kf,
655 p[0], p[1], p[2], p[3],
656 origin[0], origin[1], origin[2], tof,
657 polar[0], polar[1], polar[2],
658 kPPrimary, nt, 1., ks);
8d2cd130 659 pParent[i] = nt;
660 KeepTrack(nt);
642f15cf 661 } // PushTrack loop
8d2cd130 662 }
663 } else {
664 nc = GenerateMB();
665 } // mb ?
f913ec4f 666
667 GetSubEventTime();
8d2cd130 668
234f6d32 669 delete[] pParent;
670 delete[] pSelected;
671 delete[] trackIt;
8d2cd130 672
673 if (nc > 0) {
674 switch (fCountMode) {
675 case kCountAll:
676 // printf(" Count all \n");
677 jev += nc;
678 break;
679 case kCountParents:
680 // printf(" Count parents \n");
681 jev += nParents;
682 break;
683 case kCountTrackables:
684 // printf(" Count trackable \n");
685 jev += nTkbles;
686 break;
687 }
688 if (jev >= fNpart || fNpart == -1) {
689 fKineBias=Float_t(fNpart)/Float_t(fTrials);
e0e89f40 690
8d2cd130 691 fQ += fPythia->GetVINT(51);
692 fX1 += fPythia->GetVINT(41);
693 fX2 += fPythia->GetVINT(42);
694 fTrialsRun += fTrials;
695 fNev++;
696 MakeHeader();
697 break;
698 }
699 }
700 } // event loop
701 SetHighWaterMark(nt);
702// adjust weight due to kinematic selection
b88f5cea 703// AdjustWeights();
8d2cd130 704// get cross-section
705 fXsection=fPythia->GetPARI(1);
706}
707
708Int_t AliGenPythia::GenerateMB()
709{
710//
711// Min Bias selection and other global selections
712//
713 Int_t i, kf, nt, iparent;
714 Int_t nc = 0;
bf950da8 715 Float_t p[4];
8d2cd130 716 Float_t polar[3] = {0,0,0};
717 Float_t origin[3] = {0,0,0};
718// converts from mm/c to s
719 const Float_t kconv=0.001/2.999792458e8;
720
e0e89f40 721
722
5fa4b20b 723 Int_t np = (fHadronisation) ? fParticles->GetEntriesFast() : fNpartons;
aea21c57 724
5fa4b20b 725
e0e89f40 726
8d2cd130 727 Int_t* pParent = new Int_t[np];
728 for (i=0; i< np; i++) pParent[i] = -1;
729 if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
730 TParticle* jet1 = (TParticle *) fParticles->At(6);
731 TParticle* jet2 = (TParticle *) fParticles->At(7);
234f6d32 732 if (!CheckTrigger(jet1, jet2)) {
733 delete [] pParent;
734 return 0;
735 }
8d2cd130 736 }
e0e89f40 737
7ce1321b 738 if (fTriggerParticle) {
739 Bool_t triggered = kFALSE;
740 for (i = 0; i < np; i++) {
741 TParticle * iparticle = (TParticle *) fParticles->At(i);
742 kf = CheckPDGCode(iparticle->GetPdgCode());
38db0ee6 743 if (kf != fTriggerParticle) continue;
7ce1321b 744 if (iparticle->Pt() == 0.) continue;
745 if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
746 triggered = kTRUE;
747 break;
748 }
234f6d32 749 if (!triggered) {
750 delete [] pParent;
751 return 0;
752 }
7ce1321b 753 }
e0e89f40 754
755
756 // Check if there is a ccbar or bbbar pair with at least one of the two
757 // in fYMin < y < fYMax
758 if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi) {
759 TParticle *hvq;
760 Bool_t theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
761 Float_t yQ;
762 Int_t pdgQ;
763 for(i=0; i<np; i++) {
764 hvq = (TParticle*)fParticles->At(i);
765 pdgQ = hvq->GetPdgCode();
766 if(TMath::Abs(pdgQ) != fFlavorSelect) continue;
767 if(pdgQ>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
768 yQ = 0.5*TMath::Log((hvq->Energy()+hvq->Pz()+1.e-13)/
769 (hvq->Energy()-hvq->Pz()+1.e-13));
770 if(yQ>fYMin && yQ<fYMax) inYcut=kTRUE;
771 }
772 if (!theQ || !theQbar || !inYcut) {
234f6d32 773 delete[] pParent;
e0e89f40 774 return 0;
775 }
776 }
aea21c57 777
0f6ee828 778 //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
779 if ( (fProcess == kPyW || fProcess == kPyZ || fProcess == kPyMb || fProcess == kPyMbNonDiffr)
780 && (fCutOnChild == 1) ) {
781 if ( !CheckKinematicsOnChild() ) {
234f6d32 782 delete[] pParent;
0f6ee828 783 return 0;
784 }
aea21c57 785 }
786
787
f913ec4f 788 for (i = 0; i < np; i++) {
8d2cd130 789 Int_t trackIt = 0;
790 TParticle * iparticle = (TParticle *) fParticles->At(i);
791 kf = CheckPDGCode(iparticle->GetPdgCode());
792 Int_t ks = iparticle->GetStatusCode();
793 Int_t km = iparticle->GetFirstMother();
794 if ((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) ||
795 (ks != 1) ||
796 (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
797 nc++;
798 if (ks == 1) trackIt = 1;
799 Int_t ipa = iparticle->GetFirstMother()-1;
800
801 iparent = (ipa > -1) ? pParent[ipa] : -1;
802
803//
804// store track information
805 p[0] = iparticle->Px();
806 p[1] = iparticle->Py();
807 p[2] = iparticle->Pz();
a920faf9 808 p[3] = iparticle->Energy();
1406f599 809
a920faf9 810
2590ccf9 811 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
812 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
813 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
814
f913ec4f 815 Float_t tof = fEventTime + kconv * iparticle->T();
a920faf9 816
817 PushTrack(fTrackIt*trackIt, iparent, kf,
818 p[0], p[1], p[2], p[3],
819 origin[0], origin[1], origin[2], tof,
820 polar[0], polar[1], polar[2],
821 kPPrimary, nt, 1., ks);
5fa4b20b 822 //
823 // Special Treatment to store color-flow
824 //
825 if (ks == 3 || ks == 13 || ks == 14) {
826 TParticle* particle = 0;
827 if (fStack) {
828 particle = fStack->Particle(nt);
829 } else {
830 particle = gAlice->Stack()->Particle(nt);
831 }
832 particle->SetFirstDaughter(fPythia->GetK(2, i));
833 particle->SetLastDaughter(fPythia->GetK(3, i));
834 }
835
8d2cd130 836 KeepTrack(nt);
837 pParent[i] = nt;
f913ec4f 838 SetHighWaterMark(nt);
839
8d2cd130 840 } // select particle
841 } // particle loop
842
234f6d32 843 delete[] pParent;
e0e89f40 844
f913ec4f 845 return 1;
8d2cd130 846}
847
848
849void AliGenPythia::FinishRun()
850{
851// Print x-section summary
852 fPythia->Pystat(1);
2779fc64 853
854 if (fNev > 0.) {
855 fQ /= fNev;
856 fX1 /= fNev;
857 fX2 /= fNev;
858 }
859
8d2cd130 860 printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
861 printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
8d2cd130 862}
863
7184e472 864void AliGenPythia::AdjustWeights() const
8d2cd130 865{
866// Adjust the weights after generation of all events
867//
e2bddf81 868 if (gAlice) {
869 TParticle *part;
870 Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
871 for (Int_t i=0; i<ntrack; i++) {
872 part= gAlice->GetMCApp()->Particle(i);
873 part->SetWeight(part->GetWeight()*fKineBias);
874 }
8d2cd130 875 }
876}
877
878void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
879{
880// Treat protons as inside nuclei with mass numbers a1 and a2
1d568bc2 881
1a626d4e 882 fAProjectile = a1;
883 fATarget = a2;
1d568bc2 884 fSetNuclei = kTRUE;
8d2cd130 885}
886
887
888void AliGenPythia::MakeHeader()
889{
7184e472 890//
891// Make header for the simulated event
892//
183a5ca9 893 if (gAlice) {
894 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
f913ec4f 895 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
183a5ca9 896 }
897
8d2cd130 898// Builds the event header, to be called after each event
e5c87a3d 899 if (fHeader) delete fHeader;
900 fHeader = new AliGenPythiaEventHeader("Pythia");
8d2cd130 901//
902// Event type
e5c87a3d 903 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
8d2cd130 904//
905// Number of trials
e5c87a3d 906 ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
8d2cd130 907//
908// Event Vertex
d25cfd65 909 fHeader->SetPrimaryVertex(fVertex);
8d2cd130 910//
911// Jets that have triggered
f913ec4f 912
8d2cd130 913 if (fProcess == kPyJets)
914 {
915 Int_t ntrig, njet;
916 Float_t jets[4][10];
917 GetJets(njet, ntrig, jets);
9ff6c04c 918
8d2cd130 919
920 for (Int_t i = 0; i < ntrig; i++) {
e5c87a3d 921 ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
8d2cd130 922 jets[3][i]);
923 }
924 }
5fa4b20b 925//
926// Copy relevant information from external header, if present.
927//
928 Float_t uqJet[4];
929
930 if (fRL) {
931 AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
932 for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
933 {
934 printf("Adding Jet %d %d \n", i, exHeader->NTriggerJets());
935
936
937 exHeader->TriggerJet(i, uqJet);
938 ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
939 }
940 }
941//
942// Store quenching parameters
943//
944 if (fQuench){
945 Double_t z[4];
946 Double_t xp, yp;
7c21f297 947 if (fQuench == 1) {
948 // Pythia::Quench()
949 fPythia->GetQuenchingParameters(xp, yp, z);
950 } else {
951 // Pyquen
952 Double_t r1 = PARIMP.rb1;
953 Double_t r2 = PARIMP.rb2;
954 Double_t b = PARIMP.b1;
955 Double_t r = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
956 Double_t phi = PARIMP.psib1;
957 xp = r * TMath::Cos(phi);
958 yp = r * TMath::Sin(phi);
959
960 }
961 ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
962 ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
963 }
beac474c 964//
965// Store pt^hard
966 ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
5fa4b20b 967//
cf57b268 968// Pass header
5fa4b20b 969//
cf57b268 970 AddHeader(fHeader);
8d2cd130 971}
cf57b268 972
973void AliGenPythia::AddHeader(AliGenEventHeader* header)
974{
975 // Add header to container or runloader
976 if (fContainer) {
977 fContainer->AddHeader(header);
978 } else {
979 AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(header);
980 }
981}
982
8d2cd130 983
984Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
985{
986// Check the kinematic trigger condition
987//
988 Double_t eta[2];
989 eta[0] = jet1->Eta();
990 eta[1] = jet2->Eta();
991 Double_t phi[2];
992 phi[0] = jet1->Phi();
993 phi[1] = jet2->Phi();
994 Int_t pdg[2];
995 pdg[0] = jet1->GetPdgCode();
996 pdg[1] = jet2->GetPdgCode();
997 Bool_t triggered = kFALSE;
998
999 if (fProcess == kPyJets) {
1000 Int_t njets = 0;
1001 Int_t ntrig = 0;
1002 Float_t jets[4][10];
1003//
1004// Use Pythia clustering on parton level to determine jet axis
1005//
1006 GetJets(njets, ntrig, jets);
1007
76d6ba9a 1008 if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
8d2cd130 1009//
1010 } else {
1011 Int_t ij = 0;
1012 Int_t ig = 1;
1013 if (pdg[0] == kGamma) {
1014 ij = 1;
1015 ig = 0;
1016 }
1017 //Check eta range first...
1018 if ((eta[ij] < fEtaMaxJet && eta[ij] > fEtaMinJet) &&
1019 (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1020 {
1021 //Eta is okay, now check phi range
1022 if ((phi[ij] < fPhiMaxJet && phi[ij] > fPhiMinJet) &&
1023 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1024 {
1025 triggered = kTRUE;
1026 }
1027 }
1028 }
1029 return triggered;
1030}
aea21c57 1031
1032
aea21c57 1033
7184e472 1034Bool_t AliGenPythia::CheckKinematicsOnChild(){
1035//
1036//Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1037//
aea21c57 1038 Bool_t checking = kFALSE;
1039 Int_t j, kcode, ks, km;
1040 Int_t nPartAcc = 0; //number of particles in the acceptance range
1041 Int_t numberOfAcceptedParticles = 1;
1042 if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
1043 Int_t npart = fParticles->GetEntriesFast();
1044
0f6ee828 1045 for (j = 0; j<npart; j++) {
aea21c57 1046 TParticle * jparticle = (TParticle *) fParticles->At(j);
1047 kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1048 ks = jparticle->GetStatusCode();
1049 km = jparticle->GetFirstMother();
1050
1051 if( (ks == 1) && (kcode == fPdgCodeParticleforAcceptanceCut) && (KinematicSelection(jparticle,1)) ){
1052 nPartAcc++;
1053 }
0f6ee828 1054 if( numberOfAcceptedParticles <= nPartAcc){
1055 checking = kTRUE;
1056 break;
1057 }
aea21c57 1058 }
0f6ee828 1059
aea21c57 1060 return checking;
aea21c57 1061}
1062
8d2cd130 1063
1064AliGenPythia& AliGenPythia::operator=(const AliGenPythia& rhs)
1065{
1066// Assignment operator
014a9521 1067 rhs.Copy(*this);
8d2cd130 1068 return *this;
1069}
1070
5fa4b20b 1071void AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
8d2cd130 1072{
1073//
1074// Load event into Pythia Common Block
1075//
5fa4b20b 1076
32d6ef7d 1077 Int_t npart = stack -> GetNprimary();
1078 Int_t n0 = 0;
1079
1080 if (!flag) {
1081 (fPythia->GetPyjets())->N = npart;
1082 } else {
1083 n0 = (fPythia->GetPyjets())->N;
1084 (fPythia->GetPyjets())->N = n0 + npart;
1085 }
1086
1087
8d2cd130 1088 for (Int_t part = 0; part < npart; part++) {
7184e472 1089 TParticle *mPart = stack->Particle(part);
32d6ef7d 1090
7184e472 1091 Int_t kf = mPart->GetPdgCode();
1092 Int_t ks = mPart->GetStatusCode();
1093 Int_t idf = mPart->GetFirstDaughter();
1094 Int_t idl = mPart->GetLastDaughter();
5fa4b20b 1095
1096 if (reHadr) {
1097 if (ks == 11 || ks == 12) {
1098 ks -= 10;
1099 idf = -1;
1100 idl = -1;
1101 }
1102 }
32d6ef7d 1103
7184e472 1104 Float_t px = mPart->Px();
1105 Float_t py = mPart->Py();
1106 Float_t pz = mPart->Pz();
1107 Float_t e = mPart->Energy();
1108 Float_t m = mPart->GetCalcMass();
8d2cd130 1109
1110
32d6ef7d 1111 (fPythia->GetPyjets())->P[0][part+n0] = px;
1112 (fPythia->GetPyjets())->P[1][part+n0] = py;
1113 (fPythia->GetPyjets())->P[2][part+n0] = pz;
1114 (fPythia->GetPyjets())->P[3][part+n0] = e;
1115 (fPythia->GetPyjets())->P[4][part+n0] = m;
8d2cd130 1116
32d6ef7d 1117 (fPythia->GetPyjets())->K[1][part+n0] = kf;
1118 (fPythia->GetPyjets())->K[0][part+n0] = ks;
5fa4b20b 1119 (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1120 (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
7184e472 1121 (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
8d2cd130 1122 }
1123}
1124
5fa4b20b 1125
014a9521 1126void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
8d2cd130 1127{
1128//
1129// Calls the Pythia jet finding algorithm to find jets in the current event
1130//
1131//
8d2cd130 1132//
1133// Save jets
1134 Int_t n = fPythia->GetN();
1135
1136//
1137// Run Jet Finder
1138 fPythia->Pycell(njets);
1139 Int_t i;
1140 for (i = 0; i < njets; i++) {
1141 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1142 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1143 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1144 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1145
1146 jets[0][i] = px;
1147 jets[1][i] = py;
1148 jets[2][i] = pz;
1149 jets[3][i] = e;
1150 }
1151}
1152
1153
1154
1155void AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1156{
1157//
1158// Calls the Pythia clustering algorithm to find jets in the current event
1159//
1160 Int_t n = fPythia->GetN();
1161 nJets = 0;
1162 nJetsTrig = 0;
1163 if (fJetReconstruction == kCluster) {
1164//
1165// Configure cluster algorithm
1166//
1167 fPythia->SetPARU(43, 2.);
1168 fPythia->SetMSTU(41, 1);
1169//
1170// Call cluster algorithm
1171//
1172 fPythia->Pyclus(nJets);
1173//
1174// Loading jets from common block
1175//
1176 } else {
592f8307 1177
8d2cd130 1178//
1179// Run Jet Finder
1180 fPythia->Pycell(nJets);
1181 }
1182
1183 Int_t i;
1184 for (i = 0; i < nJets; i++) {
1185 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
1186 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
1187 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
1188 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
1189 Float_t pt = TMath::Sqrt(px * px + py * py);
a920faf9 1190 Float_t phi = TMath::Pi() + TMath::ATan2(-py, -px);
8d2cd130 1191 Float_t theta = TMath::ATan2(pt,pz);
1192 Float_t et = e * TMath::Sin(theta);
1193 Float_t eta = -TMath::Log(TMath::Tan(theta / 2.));
8d2cd130 1194 if (
1195 eta > fEtaMinJet && eta < fEtaMaxJet &&
675eb105 1196 phi > fPhiMinJet && phi < fPhiMaxJet &&
8d2cd130 1197 et > fEtMinJet && et < fEtMaxJet
1198 )
1199 {
1200 jets[0][nJetsTrig] = px;
1201 jets[1][nJetsTrig] = py;
1202 jets[2][nJetsTrig] = pz;
1203 jets[3][nJetsTrig] = e;
1204 nJetsTrig++;
5fa4b20b 1205// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
8d2cd130 1206 } else {
1207// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1208 }
1209 }
1210}
1211
f913ec4f 1212void AliGenPythia::GetSubEventTime()
1213{
1214 // Calculates time of the next subevent
9d7108a7 1215 fEventTime = 0.;
1216 if (fEventsTime) {
1217 TArrayF &array = *fEventsTime;
1218 fEventTime = array[fCurSubEvent++];
1219 }
1220 // printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1221 return;
f913ec4f 1222}
8d2cd130 1223
1224#ifdef never
1225void AliGenPythia::Streamer(TBuffer &R__b)
1226{
1227 // Stream an object of class AliGenPythia.
1228
1229 if (R__b.IsReading()) {
1230 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1231 AliGenerator::Streamer(R__b);
1232 R__b >> (Int_t&)fProcess;
1233 R__b >> (Int_t&)fStrucFunc;
1234 R__b >> (Int_t&)fForceDecay;
1235 R__b >> fEnergyCMS;
1236 R__b >> fKineBias;
1237 R__b >> fTrials;
1238 fParentSelect.Streamer(R__b);
1239 fChildSelect.Streamer(R__b);
1240 R__b >> fXsection;
1241// (AliPythia::Instance())->Streamer(R__b);
1242 R__b >> fPtHardMin;
1243 R__b >> fPtHardMax;
1244// if (fDecayer) fDecayer->Streamer(R__b);
1245 } else {
1246 R__b.WriteVersion(AliGenPythia::IsA());
1247 AliGenerator::Streamer(R__b);
1248 R__b << (Int_t)fProcess;
1249 R__b << (Int_t)fStrucFunc;
1250 R__b << (Int_t)fForceDecay;
1251 R__b << fEnergyCMS;
1252 R__b << fKineBias;
1253 R__b << fTrials;
1254 fParentSelect.Streamer(R__b);
1255 fChildSelect.Streamer(R__b);
1256 R__b << fXsection;
1257// R__b << fPythia;
1258 R__b << fPtHardMin;
1259 R__b << fPtHardMax;
1260 // fDecayer->Streamer(R__b);
1261 }
1262}
1263#endif
1264
90d7b703 1265
589380c6 1266