TDPMjet/AliGenDPMjet.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
  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>
  31 #include <TClonesArray.h>
  32 #include "AliRunLoader.h"
  33 #include "AliGenDPMjet.h"
  34 #include "AliGenDPMjetEventHeader.h"
  35 #include "AliRun.h"
  36 #include "AliDpmJetRndm.h"
  37 #include "AliHeader.h"
  38 #include "AliStack.h"
  39 #include "AliMC.h"
  40
  41 ClassImp(AliGenDPMjet)
  42
  43 //______________________________________________________________________________
  44 AliGenDPMjet::AliGenDPMjet()
  45     :AliGenMC(),
  46      fBeamEn(2750.),
  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),
  57      fNoGammas(0),
  58      fLHC(0),
  59      fPi0Decay(1),
  60      fDecayAll(0),
  61      fGenImpPar(0.),
  62      fProcess(kDpmMb),
  63      fTriggerMultiplicity(0),
  64      fTriggerMultiplicityEta(0),
  65      fTriggerMultiplicityPtMin(0),
  66      fkTuneForDiff(0),
  67      fProcDiff(0)
  68 {
  69 // Constructor
  70     fEnergyCMS = 5500.;
  71     AliDpmJetRndm::SetDpmJetRandom(GetRandom());
  72 }
  73
  74
  75 //______________________________________________________________________________
  76 AliGenDPMjet::AliGenDPMjet(Int_t npart)
  77     :AliGenMC(npart),
  78      fBeamEn(2750.),
  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),
  89      fNoGammas(0),
  90      fLHC(0),
  91      fPi0Decay(1),
  92      fDecayAll(0),
  93      fGenImpPar(0.),
  94      fProcess(kDpmMb),
  95      fTriggerMultiplicity(0),
  96      fTriggerMultiplicityEta(0),
  97      fTriggerMultiplicityPtMin(0),
  98      fkTuneForDiff(0),
  99      fProcDiff(0)
 100 {
 101 // Default PbPb collisions at 5. 5 TeV
 102 //
 103     fEnergyCMS = 5500.;
 104     fName = "DPMJET";
 105     fTitle= "Particle Generator using DPMJET";
 106     SetTarget();
 107     SetProjectile();
 108     fVertex.Set(3);
 109     AliDpmJetRndm::SetDpmJetRandom(GetRandom());
 110 }
 111
 112 AliGenDPMjet::AliGenDPMjet(const AliGenDPMjet &/*Dpmjet*/)
 113     :AliGenMC(),
 114      fBeamEn(2750.),
 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),
 125      fNoGammas(0),
 126      fLHC(0),
 127      fPi0Decay(1),
 128      fDecayAll(0),
 129      fGenImpPar(0.),
 130      fProcess(kDpmMb),
 131      fTriggerMultiplicity(0),
 132      fTriggerMultiplicityEta(0),
 133      fTriggerMultiplicityPtMin(0),
 134      fkTuneForDiff(0),
 135      fProcDiff(0)
 136 {
 137     // Dummy copy constructor
 138     fEnergyCMS = 5500.;
 139 }
 140
 141 //______________________________________________________________________________
 142 AliGenDPMjet::~AliGenDPMjet()
 143 {
 144 // Destructor
 145 }
 146 //______________________________________________________________________________
 147 void AliGenDPMjet::Init()
 148 {
 149 // Initialization
 150
 151     SetMC(new TDPMjet(fProcess, fAProjectile, fZProjectile, fATarget, fZTarget,
 152                       fBeamEn,fEnergyCMS));
 153
 154     fDPMjet=(TDPMjet*) fMCEvGen;
 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
 162     if (fAProjectile == 1 && TMath::Abs(fZProjectile == 1)) fDPMjet->SetfIdp(1);
 163
 164     fDPMjet->SetfFCentr(fICentr);
 165     fDPMjet->SetbRange(fMinImpactParam, fMaxImpactParam);
 166     fDPMjet->SetPi0Decay(fPi0Decay);
 167     fDPMjet->SetDecayAll(fDecayAll);
 168 //
 169 //  Initialize DPMjet
 170 //
 171     fDPMjet->Initialize();
 172 }
 173
 174
 175 //______________________________________________________________________________
 176 void AliGenDPMjet::Generate()
 177 {
 178 // Generate one event
 179
 180   Double_t polar[3]    =   {0,0,0};
 181   Double_t origin[3]   =   {0,0,0};
 182   Double_t p[4]        =   {0};
 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;
 189   Int_t kf, ks, imo;
 190   kf = 0;
 191   fTrials = 0;
 192   //  Set collision vertex position
 193   if (fVertexSmear == kPerEvent) Vertex();
 194
 195   while(1)
 196   {
 197 //    Generate one event
 198 // --------------------------------------------------------------------------
 199       fSpecn = 0;
 200       fSpecp = 0;
 201 // --------------------------------------------------------------------------
 202       fDPMjet->GenerateEvent();
 203
 204       fTrials++;
 205
 206       fDPMjet->ImportParticles(&fParticles,"All");
 207       if (fLHC) Boost();
 208
 209       // Temporaneo
 210       fGenImpPar = fDPMjet->GetBImpac();
 211
 212       if(TMath::Abs(fXingAngleY) > 1.e-10) BeamCrossAngle();
 213
 214       Int_t np = fParticles.GetEntriesFast();
 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
 245
 246       if(fkTuneForDiff && (TMath::Abs(fEnergyCMS - 900) < 1)) {
 247         if(!CheckDiffraction() ) continue;
 248       }
 249
 250
 251       Int_t nc = 0;
 252       if (np == 0) continue;
 253
 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
 263 //      First select parent particles
 264
 265       for (i = 0; i<np; i++) {
 266           TParticle *iparticle = (TParticle *) fParticles.At(i);
 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();
 276           if (kf == 92 || kf == 99999) continue;
 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);
 283
 284
 285           hasSelectedDaughters = DaughtersSelection(iparticle);
 286
 287
 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++) {
 301           TParticle *iparticle = (TParticle *) fParticles.At(i);
 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++) {
 335           TParticle *  iparticle = (TParticle *) fParticles.At(i);
 336           Bool_t  hasMother   = (iparticle->GetFirstMother()>=0);
 337           if (pSelected[i]) {
 338
 339               kf   = iparticle->GetPdgCode();
 340               ks   = iparticle->GetStatusCode();
 341
 342               p[0] = iparticle->Px();
 343               p[1] = iparticle->Py();
 344               p[2] = iparticle->Pz();
 345               p[3] = iparticle->Energy();
 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
 350               tof = fTime + kconv*iparticle->T();
 351
 352               imo = -1;
 353               TParticle* mother = 0;
 354               if (hasMother) {
 355                   imo = iparticle->GetFirstMother();
 356                   mother = (TParticle *) fParticles.At(imo);
 357                   imo = (mother->GetPdgCode() != 92 && mother->GetPdgCode() != 99999) ? newPos[imo] : -1;
 358               } // if has mother
 359
 360
 361
 362               Bool_t tFlag = (fTrackIt && (ks == 1));
 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);
 368               KeepTrack(nt);
 369               newPos[i] = nt;
 370           } // if selected
 371       } // particle loop
 372       delete[] newPos;
 373       delete[] pSelected;
 374       if (nc>0) {
 375           jev += nc;
 376           if (jev >= fNpart || fNpart == -1) {
 377               break;
 378           }
 379       }
 380   } // event loop
 381   MakeHeader();
 382   SetHighWaterMark(nt);
 383 }
 384
 385 //______________________________________________________________________________
 386 Bool_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++){
 401             TParticle *  jparticle = (TParticle *) fParticles.At(i);
 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 //______________________________________________________________________________
 417 Bool_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 //______________________________________________________________________________
 441 Bool_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
 452 //______________________________________________________________________________
 453 void AliGenDPMjet::MakeHeader()
 454 {
 455 //  printf("MakeHeader %13.3f \n", fDPMjet->GetBImpac());
 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());
 461     ((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetProjParticipants(),
 462                                                          fDPMjet->GetTargParticipants());
 463
 464     if(fProcDiff>0){
 465     ((AliGenDPMjetEventHeader*) header)->SetProcessType(fProcDiff);
 466     }
 467     else
 468       ((AliGenDPMjetEventHeader*) header)->SetProcessType(fDPMjet->GetProcessCode());
 469
 470     // Bookkeeping for kinematic bias
 471     ((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
 472     // Event Vertex
 473     header->SetPrimaryVertex(fVertex);
 474     header->SetInteractionTime(fTime);
 475     gAlice->SetGenEventHeader(header);
 476     AddHeader(header);
 477 }
 478
 479 void AliGenDPMjet::AddHeader(AliGenEventHeader* header)
 480 {
 481     // Add header to container or runloader
 482     if (fContainer) {
 483         fContainer->AddHeader(header);
 484     } else {
 485         AliRunLoader::Instance()->GetHeader()->SetGenEventHeader(header);
 486     }
 487 }
 488
 489
 490 //______________________________________________________________________________
 491 AliGenDPMjet& AliGenDPMjet::operator=(const  AliGenDPMjet& /*rhs*/)
 492 {
 493 // Assignment operator
 494     return *this;
 495 }
 496
 497
 498 void AliGenDPMjet::FinishRun()
 499 {
 500     // Print run statistics
 501     fDPMjet->Dt_Dtuout();
 502 }
 503
 504
 505
 506 Bool_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   };
 542
 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
 618 const Int_t nbin=120;
 619 Double_t bin[]={
 620 1.080000, 1.274258, 1.468516, 1.662773, 1.857031, 2.051289,
 621 2.245547, 2.439805, 2.634062, 2.828320, 3.022578, 3.216836,
 622 3.411094, 3.605352, 3.799609, 3.993867, 4.188125, 4.382383,
 623 4.576641, 4.770898, 4.965156, 5.547930, 6.130703, 6.713477,
 624 7.296250, 7.879023, 8.461797, 9.044570, 9.627344, 10.210117,
 625 10.792891, 11.375664, 11.958437, 12.541211, 13.123984, 13.706758,
 626 14.289531, 14.872305, 15.455078, 16.037852, 16.620625, 17.203398,
 627 17.786172, 18.368945, 18.951719, 19.534492, 20.117266, 20.700039,
 628 21.282812, 21.865586, 22.448359, 23.031133, 23.613906, 24.196680,
 629 24.779453, 25.362227, 25.945000, 26.527773, 27.110547, 27.693320,
 630 28.276094, 28.858867, 29.441641, 30.024414, 30.607187, 31.189961,
 631 31.772734, 32.355508, 32.938281, 33.521055, 34.103828, 34.686602,
 632 35.269375, 35.852148, 36.434922, 37.017695, 37.600469, 38.183242,
 633 38.766016, 39.348789, 39.931562, 40.514336, 41.097109, 41.679883,
 634 42.262656, 42.845430, 43.428203, 44.010977, 44.593750, 45.176523,
 635 45.759297, 46.342070, 46.924844, 47.507617, 48.090391, 48.673164,
 636 49.255937, 49.838711, 50.421484, 57.220508, 64.019531, 70.818555,
 637 77.617578, 84.416602, 91.215625, 98.014648, 104.813672, 111.612695,
 638 118.411719, 125.210742, 132.009766, 138.808789, 145.607812, 152.406836,
 639 159.205859, 166.004883, 172.803906, 179.602930, 186.401953, 193.200977,
 640 200.000000};
 641 Double_t w[]={
 642 1.000000, 0.367136, 0.239268, 0.181139, 0.167470, 0.160072,
 643 0.147832, 0.162765, 0.176103, 0.156382, 0.146040, 0.143375,
 644 0.134038, 0.126747, 0.123152, 0.119424, 0.113839, 0.109433,
 645 0.107180, 0.104690, 0.096427, 0.090603, 0.083706, 0.077206,
 646 0.074603, 0.069698, 0.067315, 0.064980, 0.063560, 0.059573,
 647 0.058712, 0.057581, 0.055944, 0.055442, 0.053272, 0.051769,
 648 0.051672, 0.049284, 0.048980, 0.048797, 0.047434, 0.047039,
 649 0.046395, 0.046227, 0.044288, 0.044743, 0.043772, 0.043902,
 650 0.042771, 0.043232, 0.042222, 0.041668, 0.041988, 0.040858,
 651 0.039672, 0.040069, 0.040274, 0.039438, 0.039903, 0.039083,
 652 0.038741, 0.038182, 0.037664, 0.038610, 0.038759, 0.038688,
 653 0.038039, 0.038220, 0.038145, 0.037445, 0.036765, 0.037333,
 654 0.036753, 0.036405, 0.036339, 0.037659, 0.036139, 0.036706,
 655 0.035393, 0.037136, 0.036570, 0.035234, 0.036832, 0.035560,
 656 0.035509, 0.035579, 0.035100, 0.035471, 0.035421, 0.034494,
 657 0.035596, 0.034935, 0.035810, 0.034324, 0.035355, 0.034323,
 658 0.033486, 0.034622, 0.034805, 0.034419, 0.033946, 0.033927,
 659 0.034224, 0.033942, 0.034088, 0.034190, 0.034620, 0.035294,
 660 0.035650, 0.035378, 0.036028, 0.035933, 0.036753, 0.037171,
 661 0.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
 719
 720 //______________________________________________________________________________