]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TDPMjet/AliGenDPMjet.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Added more functions and fixed blast wave
[u/mrichter/AliRoot.git] / TDPMjet / AliGenDPMjet.cxx
CommitLineData
1a52e0ed 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
17// Generator using DPMJET as an external generator
18// The main DPMJET options are accessable for the user through this interface.
19// Uses the TDPMjet implementation of TGenerator.
20
21#include <TDPMjet.h>
22#include <TRandom.h>
23#include <TArrayI.h>
24#include <TParticle.h>
25#include <TGraph.h>
26#include <TDatabasePDG.h>
27#include <TParticlePDG.h>
28#include <TParticleClassPDG.h>
29#include <TPDGCode.h>
30#include <TLorentzVector.h>
07c71c2e 31#include <TClonesArray.h>
9dcc0beb 32#include "AliRunLoader.h"
1a52e0ed 33#include "AliGenDPMjet.h"
34#include "AliGenDPMjetEventHeader.h"
1a52e0ed 35#include "AliRun.h"
f2ba02a4 36#include "AliDpmJetRndm.h"
9dcc0beb 37#include "AliHeader.h"
38#include "AliStack.h"
39#include "AliMC.h"
1a52e0ed 40
f97d4a8e 41ClassImp(AliGenDPMjet)
1a52e0ed 42
43//______________________________________________________________________________
44AliGenDPMjet::AliGenDPMjet()
f97d4a8e 45 :AliGenMC(),
46 fBeamEn(2750.),
f97d4a8e 47 fMinImpactParam(0.),
48 fMaxImpactParam(5.),
49 fICentr(0),
50 fSelectAll(0),
51 fFlavor(0),
52 fTrials(0),
53 fSpectators(1),
54 fSpecn(0),
55 fSpecp(0),
56 fDPMjet(0),
f97d4a8e 57 fNoGammas(0),
58 fLHC(0),
9c97cedd 59 fPi0Decay(1),
717765ce 60 fDecayAll(0),
f97d4a8e 61 fGenImpPar(0.),
3b50c39e 62 fProcess(kDpmMb),
63 fTriggerMultiplicity(0),
64 fTriggerMultiplicityEta(0),
19896767 65 fTriggerMultiplicityPtMin(0),
66 fkTuneForDiff(0),
67 fProcDiff(0)
1a52e0ed 68{
69// Constructor
e7c989e4 70 fEnergyCMS = 5500.;
f2ba02a4 71 AliDpmJetRndm::SetDpmJetRandom(GetRandom());
1a52e0ed 72}
73
74
75//______________________________________________________________________________
76AliGenDPMjet::AliGenDPMjet(Int_t npart)
f97d4a8e 77 :AliGenMC(npart),
78 fBeamEn(2750.),
f97d4a8e 79 fMinImpactParam(0.),
80 fMaxImpactParam(5.),
81 fICentr(0),
82 fSelectAll(0),
83 fFlavor(0),
84 fTrials(0),
85 fSpectators(1),
86 fSpecn(0),
87 fSpecp(0),
88 fDPMjet(0),
f97d4a8e 89 fNoGammas(0),
90 fLHC(0),
9c97cedd 91 fPi0Decay(1),
717765ce 92 fDecayAll(0),
f97d4a8e 93 fGenImpPar(0.),
3b50c39e 94 fProcess(kDpmMb),
95 fTriggerMultiplicity(0),
96 fTriggerMultiplicityEta(0),
19896767 97 fTriggerMultiplicityPtMin(0),
98 fkTuneForDiff(0),
99 fProcDiff(0)
1a52e0ed 100{
101// Default PbPb collisions at 5. 5 TeV
102//
e7c989e4 103 fEnergyCMS = 5500.;
1a52e0ed 104 fName = "DPMJET";
105 fTitle= "Particle Generator using DPMJET";
1a52e0ed 106 SetTarget();
107 SetProjectile();
f2ba02a4 108 fVertex.Set(3);
f2ba02a4 109 AliDpmJetRndm::SetDpmJetRandom(GetRandom());
1a52e0ed 110}
111
ba758f5a 112AliGenDPMjet::AliGenDPMjet(const AliGenDPMjet &/*Dpmjet*/)
f97d4a8e 113 :AliGenMC(),
114 fBeamEn(2750.),
f97d4a8e 115 fMinImpactParam(0.),
116 fMaxImpactParam(5.),
117 fICentr(0),
118 fSelectAll(0),
119 fFlavor(0),
120 fTrials(0),
121 fSpectators(1),
122 fSpecn(0),
123 fSpecp(0),
124 fDPMjet(0),
f97d4a8e 125 fNoGammas(0),
126 fLHC(0),
9c97cedd 127 fPi0Decay(1),
717765ce 128 fDecayAll(0),
f97d4a8e 129 fGenImpPar(0.),
3b50c39e 130 fProcess(kDpmMb),
131 fTriggerMultiplicity(0),
132 fTriggerMultiplicityEta(0),
19896767 133 fTriggerMultiplicityPtMin(0),
134 fkTuneForDiff(0),
135 fProcDiff(0)
ba758f5a 136{
f97d4a8e 137 // Dummy copy constructor
e7c989e4 138 fEnergyCMS = 5500.;
ba758f5a 139}
1a52e0ed 140
141//______________________________________________________________________________
142AliGenDPMjet::~AliGenDPMjet()
143{
144// Destructor
1a52e0ed 145}
1a52e0ed 146//______________________________________________________________________________
147void AliGenDPMjet::Init()
148{
149// Initialization
150
e67e0fb7 151 SetMC(new TDPMjet(fProcess, fAProjectile, fZProjectile, fATarget, fZTarget,
1a52e0ed 152 fBeamEn,fEnergyCMS));
153
35c87559 154 fDPMjet=(TDPMjet*) fMCEvGen;
1a52e0ed 155 //
156 // **** Flag to force central production
157 // fICentr=1. central production forced
158 // fICentr<0 && fICentr>-100 -> bmin = fMinImpactParam, bmax = fMaxImpactParam
159 // fICentr<-99 -> fraction of x-sec. = XSFRAC
160 // fICentr=-1. -> evaporation/fzc suppressed
161 // fICentr<-1. -> evaporation/fzc suppressed
cd19c617 162 if (fAProjectile == 1 && TMath::Abs(fZProjectile == 1)) fDPMjet->SetfIdp(1);
ba758f5a 163
1a52e0ed 164 fDPMjet->SetfFCentr(fICentr);
e67e0fb7 165 fDPMjet->SetbRange(fMinImpactParam, fMaxImpactParam);
f97d4a8e 166 fDPMjet->SetPi0Decay(fPi0Decay);
717765ce 167 fDPMjet->SetDecayAll(fDecayAll);
1a52e0ed 168//
169// Initialize DPMjet
170//
171 fDPMjet->Initialize();
1a52e0ed 172}
173
174
175//______________________________________________________________________________
176void AliGenDPMjet::Generate()
177{
178// Generate one event
179
c0d429ef 180 Double_t polar[3] = {0,0,0};
181 Double_t origin[3] = {0,0,0};
182 Double_t p[4] = {0};
1a52e0ed 183 Float_t tof;
184
185// converts from mm/c to s
186 const Float_t kconv = 0.001/2.999792458e8;
187 Int_t nt = 0;
188 Int_t jev = 0;
e67e0fb7 189 Int_t kf, ks, imo;
1a52e0ed 190 kf = 0;
1a52e0ed 191 fTrials = 0;
e67e0fb7 192 // Set collision vertex position
193 if (fVertexSmear == kPerEvent) Vertex();
1a52e0ed 194
195 while(1)
196 {
197// Generate one event
198// --------------------------------------------------------------------------
199 fSpecn = 0;
200 fSpecp = 0;
201// --------------------------------------------------------------------------
202 fDPMjet->GenerateEvent();
9c97cedd 203
1a52e0ed 204 fTrials++;
205
8507138f 206 fDPMjet->ImportParticles(&fParticles,"All");
1a52e0ed 207 if (fLHC) Boost();
208
209 // Temporaneo
210 fGenImpPar = fDPMjet->GetBImpac();
41cc7c18 211
212 if(TMath::Abs(fXingAngleY) > 1.e-10) BeamCrossAngle();
213
8507138f 214 Int_t np = fParticles.GetEntriesFast();
3b50c39e 215 //
216 // Multiplicity Trigger
217 if (fTriggerMultiplicity > 0) {
218 Int_t multiplicity = 0;
219 for (Int_t i = 0; i < np; i++) {
220 TParticle * iparticle = (TParticle *) fParticles.At(i);
221
222 Int_t statusCode = iparticle->GetStatusCode();
223
224 // Initial state particle
225 if (statusCode != 1)
226 continue;
227 // eta cut
228 if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
229 continue;
230 // pt cut
231 if (iparticle->Pt() < fTriggerMultiplicityPtMin)
232 continue;
233
234 TParticlePDG* pdgPart = iparticle->GetPDG();
235 if (pdgPart && pdgPart->Charge() == 0)
236 continue;
237 ++multiplicity;
238 }
239 //
240 //
241 if (multiplicity < fTriggerMultiplicity) continue;
242 Printf("Triggered on event with multiplicity of %d >= %d", multiplicity, fTriggerMultiplicity);
243 }
244
19896767 245
246 if(fkTuneForDiff && (TMath::Abs(fEnergyCMS - 900) < 1)) {
247 if(!CheckDiffraction() ) continue;
248 }
249
250
3b50c39e 251 Int_t nc = 0;
252 if (np == 0) continue;
253
1a52e0ed 254 Int_t i;
255 Int_t* newPos = new Int_t[np];
256 Int_t* pSelected = new Int_t[np];
257
258 for (i = 0; i<np; i++) {
259 newPos[i] = i;
260 pSelected[i] = 0;
261 }
262
1a52e0ed 263// First select parent particles
264
265 for (i = 0; i<np; i++) {
8507138f 266 TParticle *iparticle = (TParticle *) fParticles.At(i);
1a52e0ed 267
268// Is this a parent particle ?
269
270 if (Stable(iparticle)) continue;
271
272 Bool_t selected = kTRUE;
273 Bool_t hasSelectedDaughters = kFALSE;
274
275 kf = iparticle->GetPdgCode();
cd19c617 276 if (kf == 92 || kf == 99999) continue;
1a52e0ed 277 ks = iparticle->GetStatusCode();
278// No initial state partons
279 if (ks==21) continue;
280
281 if (!fSelectAll) selected = KinematicSelection(iparticle, 0) &&
282 SelectFlavor(kf);
9c97cedd 283
284
1a52e0ed 285 hasSelectedDaughters = DaughtersSelection(iparticle);
286
9c97cedd 287
1a52e0ed 288// Put particle on the stack if it is either selected or
289// it is the mother of at least one seleted particle
290
291 if (selected || hasSelectedDaughters) {
292 nc++;
293 pSelected[i] = 1;
294 } // selected
295 } // particle loop parents
296
297// Now select the final state particles
298
299
300 for (i=0; i<np; i++) {
8507138f 301 TParticle *iparticle = (TParticle *) fParticles.At(i);
1a52e0ed 302
303// Is this a final state particle ?
304
305 if (!Stable(iparticle)) continue;
306
307 Bool_t selected = kTRUE;
308 kf = iparticle->GetPdgCode();
309 ks = iparticle->GetStatusCode();
310
311// --------------------------------------------------------------------------
312// Count spectator neutrons and protons (ks == 13, 14)
313 if(ks == 13 || ks == 14){
314 if(kf == kNeutron) fSpecn += 1;
315 if(kf == kProton) fSpecp += 1;
316 }
317// --------------------------------------------------------------------------
318
319 if (!fSelectAll) {
320 selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
321 if (!fSpectators && selected) selected = (ks == 13 || ks == 14);
322 }
323
324// Put particle on the stack if selected
325
326 if (selected) {
327 nc++;
328 pSelected[i] = 1;
329 } // selected
330 } // particle loop final state
331
332// Write particles to stack
333
334 for (i = 0; i<np; i++) {
8507138f 335 TParticle * iparticle = (TParticle *) fParticles.At(i);
1a52e0ed 336 Bool_t hasMother = (iparticle->GetFirstMother()>=0);
1a52e0ed 337 if (pSelected[i]) {
338
339 kf = iparticle->GetPdgCode();
1a52e0ed 340 ks = iparticle->GetStatusCode();
341
342 p[0] = iparticle->Px();
343 p[1] = iparticle->Py();
344 p[2] = iparticle->Pz();
c0d429ef 345 p[3] = iparticle->Energy();
e67e0fb7 346 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
347 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
348 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
349
21391258 350 tof = fTime + kconv*iparticle->T();
1a52e0ed 351
352 imo = -1;
353 TParticle* mother = 0;
354 if (hasMother) {
355 imo = iparticle->GetFirstMother();
8507138f 356 mother = (TParticle *) fParticles.At(imo);
cd19c617 357 imo = (mother->GetPdgCode() != 92 && mother->GetPdgCode() != 99999) ? newPos[imo] : -1;
1a52e0ed 358 } // if has mother
359
9c97cedd 360
361
1a52e0ed 362 Bool_t tFlag = (fTrackIt && (ks == 1));
c0d429ef 363 PushTrack(tFlag, imo, kf,
364 p[0], p[1], p[2], p[3],
365 origin[0], origin[1], origin[2], tof,
366 polar[0], polar[1], polar[2],
367 kPNoProcess, nt, 1., ks);
1a52e0ed 368 KeepTrack(nt);
369 newPos[i] = nt;
370 } // if selected
371 } // particle loop
372 delete[] newPos;
373 delete[] pSelected;
1a52e0ed 374 if (nc>0) {
375 jev += nc;
376 if (jev >= fNpart || fNpart == -1) {
1a52e0ed 377 break;
378 }
1a52e0ed 379 }
380 } // event loop
381 MakeHeader();
382 SetHighWaterMark(nt);
383}
384
1a52e0ed 385//______________________________________________________________________________
386Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle)
387{
388//
389// Looks recursively if one of the daughters has been selected
390//
391// printf("\n Consider daughters %d:",iparticle->GetPdgCode());
392 Int_t imin = -1;
393 Int_t imax = -1;
394 Int_t i;
395 Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
396 Bool_t selected = kFALSE;
397 if (hasDaughters) {
398 imin = iparticle->GetFirstDaughter();
399 imax = iparticle->GetLastDaughter();
400 for (i = imin; i <= imax; i++){
8507138f 401 TParticle * jparticle = (TParticle *) fParticles.At(i);
1a52e0ed 402 Int_t ip = jparticle->GetPdgCode();
403 if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
404 selected=kTRUE; break;
405 }
406 if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
407 }
408 } else {
409 return kFALSE;
410 }
411 return selected;
412}
413
414
415
416//______________________________________________________________________________
417Bool_t AliGenDPMjet::SelectFlavor(Int_t pid)
418{
419// Select flavor of particle
420// 0: all
421// 4: charm and beauty
422// 5: beauty
423 Bool_t res = 0;
424
425 if (fFlavor == 0) {
426 res = kTRUE;
427 } else {
428 Int_t ifl = TMath::Abs(pid/100);
429 if (ifl > 10) ifl/=10;
430 res = (fFlavor == ifl);
431 }
432//
433// This part if gamma writing is inhibited
434 if (fNoGammas)
435 res = res && (pid != kGamma && pid != kPi0);
436//
437 return res;
438}
439
440//______________________________________________________________________________
441Bool_t AliGenDPMjet::Stable(TParticle* particle)
442{
443// Return true for a stable particle
444//
445
446// if (particle->GetFirstDaughter() < 0 ) return kTRUE;
447 if (particle->GetStatusCode() == 1) return kTRUE;
448 else return kFALSE;
449
450}
451
1a52e0ed 452//______________________________________________________________________________
453void AliGenDPMjet::MakeHeader()
454{
19896767 455// printf("MakeHeader %13.3f \n", fDPMjet->GetBImpac());
1a52e0ed 456// Builds the event header, to be called after each event
457 AliGenEventHeader* header = new AliGenDPMjetEventHeader("DPMJET");
458 ((AliGenDPMjetEventHeader*) header)->SetNProduced(fDPMjet->GetNumStablePc());
459 ((AliGenDPMjetEventHeader*) header)->SetImpactParameter(fDPMjet->GetBImpac());
460 ((AliGenDPMjetEventHeader*) header)->SetTotalEnergy(fDPMjet->GetTotEnergy());
3c349853 461 ((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetProjParticipants(),
462 fDPMjet->GetTargParticipants());
19896767 463
464 if(fProcDiff>0){
465 ((AliGenDPMjetEventHeader*) header)->SetProcessType(fProcDiff);
466 }
467 else
468 ((AliGenDPMjetEventHeader*) header)->SetProcessType(fDPMjet->GetProcessCode());
469
3c349853 470 // Bookkeeping for kinematic bias
1a52e0ed 471 ((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
3c349853 472 // Event Vertex
f2ba02a4 473 header->SetPrimaryVertex(fVertex);
21391258 474 header->SetInteractionTime(fTime);
1a52e0ed 475 gAlice->SetGenEventHeader(header);
d2593e7a 476 AddHeader(header);
1a52e0ed 477}
478
9dcc0beb 479void AliGenDPMjet::AddHeader(AliGenEventHeader* header)
480{
481 // Add header to container or runloader
482 if (fContainer) {
483 fContainer->AddHeader(header);
484 } else {
33c3c91a 485 AliRunLoader::Instance()->GetHeader()->SetGenEventHeader(header);
9dcc0beb 486 }
487}
1a52e0ed 488
489
490//______________________________________________________________________________
ba758f5a 491AliGenDPMjet& AliGenDPMjet::operator=(const AliGenDPMjet& /*rhs*/)
1a52e0ed 492{
493// Assignment operator
494 return *this;
495}
496
497
d2593e7a 498void AliGenDPMjet::FinishRun()
499{
500 // Print run statistics
501 fDPMjet->Dt_Dtuout();
502}
503
19896767 504
505
506Bool_t AliGenDPMjet::CheckDiffraction()
507{
508
509 // printf("AAA\n");
510
511 Int_t np = fParticles.GetEntriesFast();
512
513 Int_t iPart1=-1;
514 Int_t iPart2=-1;
515
516 Double_t y1 = 1e10;
517 Double_t y2 = -1e10;
518
519 const Int_t kNstable=20;
520 const Int_t pdgStable[20] = {
521 22, // Photon
522 11, // Electron
523 12, // Electron Neutrino
524 13, // Muon
525 14, // Muon Neutrino
526 15, // Tau
527 16, // Tau Neutrino
528 211, // Pion
529 321, // Kaon
530 311, // K0
531 130, // K0s
532 310, // K0l
533 2212, // Proton
534 2112, // Neutron
535 3122, // Lambda_0
536 3112, // Sigma Minus
537 3222, // Sigma Plus
538 3312, // Xsi Minus
539 3322, // Xsi0
540 3334 // Omega
541 };
d2593e7a 542
19896767 543 for (Int_t i = 0; i < np; i++) {
544 TParticle * part = (TParticle *) fParticles.At(i);
545
546 Int_t statusCode = part->GetStatusCode();
547
548 // Initial state particle
549 if (statusCode != 1)
550 continue;
551
552 Int_t pdg = TMath::Abs(part->GetPdgCode());
553 Bool_t isStable = kFALSE;
554 for (Int_t i1 = 0; i1 < kNstable; i1++) {
555 if (pdg == pdgStable[i1]) {
556 isStable = kTRUE;
557 break;
558 }
559 }
560 if(!isStable)
561 continue;
562
563 Double_t y = part->Y();
564
565 if (y < y1)
566 {
567 y1 = y;
568 iPart1 = i;
569 }
570 if (y > y2)
571 {
572 y2 = y;
573 iPart2 = i;
574 }
575 }
576
577 if(iPart1<0 || iPart2<0) return kFALSE;
578
579 y1=TMath::Abs(y1);
580 y2=TMath::Abs(y2);
581
582 TParticle * part1 = (TParticle *) fParticles.At(iPart1);
583 TParticle * part2 = (TParticle *) fParticles.At(iPart2);
584
585 Int_t pdg1 = part1->GetPdgCode();
586 Int_t pdg2 = part2->GetPdgCode();
587
588
589 Int_t iPart = -1;
590 if (pdg1 == 2212 && pdg2 == 2212)
591 {
592 if(y1 > y2)
593 iPart = iPart1;
594 else if(y1 < y2)
595 iPart = iPart2;
596 else {
597 iPart = iPart1;
598 if((AliDpmJetRndm::GetDpmJetRandom())->Uniform(0.,1.)>0.5) iPart = iPart2;
599 }
600 }
601 else if (pdg1 == 2212)
602 iPart = iPart1;
603 else if (pdg2 == 2212)
604 iPart = iPart2;
605
606
607
608
609
610 Double_t M=-1.;
611 if(iPart>0) {
612 TParticle * part = (TParticle *) fParticles.At(iPart);
613 Double_t E= part->Energy();
614 Double_t P= part->P();
615 M= TMath::Sqrt((fEnergyCMS-E-P)*(fEnergyCMS-E+P));
616 }
617
618const Int_t nbin=120;
619Double_t bin[]={
6201.080000, 1.274258, 1.468516, 1.662773, 1.857031, 2.051289,
6212.245547, 2.439805, 2.634062, 2.828320, 3.022578, 3.216836,
6223.411094, 3.605352, 3.799609, 3.993867, 4.188125, 4.382383,
6234.576641, 4.770898, 4.965156, 5.547930, 6.130703, 6.713477,
6247.296250, 7.879023, 8.461797, 9.044570, 9.627344, 10.210117,
62510.792891, 11.375664, 11.958437, 12.541211, 13.123984, 13.706758,
62614.289531, 14.872305, 15.455078, 16.037852, 16.620625, 17.203398,
62717.786172, 18.368945, 18.951719, 19.534492, 20.117266, 20.700039,
62821.282812, 21.865586, 22.448359, 23.031133, 23.613906, 24.196680,
62924.779453, 25.362227, 25.945000, 26.527773, 27.110547, 27.693320,
63028.276094, 28.858867, 29.441641, 30.024414, 30.607187, 31.189961,
63131.772734, 32.355508, 32.938281, 33.521055, 34.103828, 34.686602,
63235.269375, 35.852148, 36.434922, 37.017695, 37.600469, 38.183242,
63338.766016, 39.348789, 39.931562, 40.514336, 41.097109, 41.679883,
63442.262656, 42.845430, 43.428203, 44.010977, 44.593750, 45.176523,
63545.759297, 46.342070, 46.924844, 47.507617, 48.090391, 48.673164,
63649.255937, 49.838711, 50.421484, 57.220508, 64.019531, 70.818555,
63777.617578, 84.416602, 91.215625, 98.014648, 104.813672, 111.612695,
638118.411719, 125.210742, 132.009766, 138.808789, 145.607812, 152.406836,
639159.205859, 166.004883, 172.803906, 179.602930, 186.401953, 193.200977,
640200.000000};
641Double_t w[]={
6421.000000, 0.367136, 0.239268, 0.181139, 0.167470, 0.160072,
6430.147832, 0.162765, 0.176103, 0.156382, 0.146040, 0.143375,
6440.134038, 0.126747, 0.123152, 0.119424, 0.113839, 0.109433,
6450.107180, 0.104690, 0.096427, 0.090603, 0.083706, 0.077206,
6460.074603, 0.069698, 0.067315, 0.064980, 0.063560, 0.059573,
6470.058712, 0.057581, 0.055944, 0.055442, 0.053272, 0.051769,
6480.051672, 0.049284, 0.048980, 0.048797, 0.047434, 0.047039,
6490.046395, 0.046227, 0.044288, 0.044743, 0.043772, 0.043902,
6500.042771, 0.043232, 0.042222, 0.041668, 0.041988, 0.040858,
6510.039672, 0.040069, 0.040274, 0.039438, 0.039903, 0.039083,
6520.038741, 0.038182, 0.037664, 0.038610, 0.038759, 0.038688,
6530.038039, 0.038220, 0.038145, 0.037445, 0.036765, 0.037333,
6540.036753, 0.036405, 0.036339, 0.037659, 0.036139, 0.036706,
6550.035393, 0.037136, 0.036570, 0.035234, 0.036832, 0.035560,
6560.035509, 0.035579, 0.035100, 0.035471, 0.035421, 0.034494,
6570.035596, 0.034935, 0.035810, 0.034324, 0.035355, 0.034323,
6580.033486, 0.034622, 0.034805, 0.034419, 0.033946, 0.033927,
6590.034224, 0.033942, 0.034088, 0.034190, 0.034620, 0.035294,
6600.035650, 0.035378, 0.036028, 0.035933, 0.036753, 0.037171,
6610.037528, 0.037985, 0.039589, 0.039359, 0.040269, 0.040755};
662
663 Double_t wSD=1.;
664 Double_t wDD=0.100418;
665 Double_t wND=0.050277;
666
667 if(M>-1 && M<bin[0]) return kFALSE;
668 if(M>bin[nbin]) M=-1;
669
670 Int_t procType=fDPMjet->GetProcessCode();//fPythia->GetMSTI(1);
671 Int_t proc0=2;
672 if(procType== 7) proc0=1;
673 if(procType== 5 || procType== 6) proc0=0;
674
675
676 // printf("M = %f bin[nbin] = %f\n",M, bin[nbin]);
677
678 Int_t proc=2;
679 if(M>0) proc=0;
680 else if(proc0==1) proc=1;
681
682 if(proc==0 && (AliDpmJetRndm::GetDpmJetRandom())->Uniform(0.,1.) > wSD) return kFALSE;
683 if(proc==1 && (AliDpmJetRndm::GetDpmJetRandom())->Uniform(0.,1.) > wDD) return kFALSE;
684 if(proc==2 && (AliDpmJetRndm::GetDpmJetRandom())->Uniform(0.,1.) > wND) return kFALSE;
685
686
687 // if(proc==1 || proc==2) return kFALSE;
688
689 if(proc!=0) {
690 if(proc0!=0) fProcDiff = procType;
691 else fProcDiff = 1;
692 return kTRUE;
693 }
694
695 Int_t ibin=nbin-1;
696 for(Int_t i=1; i<=nbin; i++)
697 if(M<=bin[i]) {
698 ibin=i-1;
699 // printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
700 break;
701 }
702
703 // printf("w[ibin] = %f\n", w[ibin]);
704
705 if((AliDpmJetRndm::GetDpmJetRandom())->Uniform(0.,1.)> w[ibin]) return kFALSE;
706
707 // printf("iPart = %d\n", iPart);
708
709 if(iPart==iPart1) fProcDiff=5;
710 else if(iPart==iPart2) fProcDiff=6;
711 else {
712 printf("EROOR: iPart!=iPart1 && iPart!=iPart2\n");
713
714 }
715
716 return kTRUE;
717}
718
1a52e0ed 719
720//______________________________________________________________________________