Adding Domenico Colella as responsible for SPD part in TRI pp
[u/mrichter/AliRoot.git] / 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 "DPMcommon.h"
23 #include <TRandom.h>
24 #include <TArrayI.h>
25 #include <TParticle.h>
26 #include <TGraph.h>
27 #include <TDatabasePDG.h>
28 #include <TParticlePDG.h>
29 #include <TParticleClassPDG.h>
30 #include <TPDGCode.h>
31 #include <TLorentzVector.h>
32 #include <TClonesArray.h>
33 #include "AliRunLoader.h"
34 #include "AliGenDPMjet.h"
35 #include "AliGenDPMjetEventHeader.h"
36 #include "AliRun.h"
37 #include "AliDpmJetRndm.h"
38 #include "AliIonPDGCodes.h"
39 #include "AliHeader.h"
40 #include "AliStack.h"
41 #include "AliMC.h"
42 #include "AliLog.h"
43
44 ClassImp(AliGenDPMjet)
45
46 //______________________________________________________________________________
47 AliGenDPMjet::AliGenDPMjet()
48     :AliGenMC(), 
49      fBeamEn(0.),
50      fMinImpactParam(0.),
51      fMaxImpactParam(100.),
52      fICentr(0),
53      fSelectAll(0),
54      fFlavor(0),
55      fTrials(0),
56      fNprimaries(0),
57      fSpectators(1),
58      fSpecn(0),
59      fSpecp(0),
60      fDPMjet(0),
61      fNoGammas(0),
62      fLHC(0),
63      fPi0Decay(1),
64      fDecayAll(0),
65      fGenImpPar(0.),
66      fProcess(kDpmMb),
67      fTriggerParticle(0),
68      fTriggerEta(0.9),     
69      fTriggerMinPt(-1),  
70      fTriggerMaxPt(1000),  
71      fTriggerMultiplicity(0),
72      fTriggerMultiplicityEta(0),
73      fTriggerMultiplicityPtMin(0),
74      fkTuneForDiff(0),
75      fProcDiff(0),
76      fFragmentation(kFALSE),
77      fHeader(AliGenDPMjetEventHeader("DPMJET"))
78
79 {
80 // Constructor
81     fEnergyCMS = 5500.;
82     AliDpmJetRndm::SetDpmJetRandom(GetRandom());
83 }
84
85
86 //______________________________________________________________________________
87 AliGenDPMjet::AliGenDPMjet(Int_t npart)
88     :AliGenMC(npart),
89      fBeamEn(0.),
90      fMinImpactParam(0.),
91      fMaxImpactParam(100.),
92      fICentr(0),
93      fSelectAll(0),
94      fFlavor(0),
95      fTrials(0),
96      fNprimaries(0),
97      fSpectators(1),
98      fSpecn(0),
99      fSpecp(0),
100      fDPMjet(0),
101      fNoGammas(0),
102      fLHC(0),
103      fPi0Decay(1),
104      fDecayAll(0),
105      fGenImpPar(0.),
106      fProcess(kDpmMb),
107      fTriggerParticle(0),
108      fTriggerEta(0.9),     
109      fTriggerMinPt(-1),  
110      fTriggerMaxPt(1000),  
111      fTriggerMultiplicity(0),
112      fTriggerMultiplicityEta(0),
113      fTriggerMultiplicityPtMin(0),
114      fkTuneForDiff(0),
115      fProcDiff(0),
116      fFragmentation(kFALSE),
117      fHeader(AliGenDPMjetEventHeader("DPMJET"))
118
119
120 {
121 // Default PbPb collisions at 5. 5 TeV
122 //
123     fEnergyCMS = 5500.;
124     fName = "DPMJET";
125     fTitle= "Particle Generator using DPMJET";
126     SetTarget();
127     SetProjectile();
128     fVertex.Set(3);
129     AliDpmJetRndm::SetDpmJetRandom(GetRandom());
130 }
131
132 AliGenDPMjet::AliGenDPMjet(const AliGenDPMjet &/*Dpmjet*/)
133     :AliGenMC(),
134      fBeamEn(0.),
135      fMinImpactParam(0.),
136      fMaxImpactParam(100.),
137      fICentr(0),
138      fSelectAll(0),
139      fFlavor(0),
140      fTrials(0),
141      fNprimaries(0),
142      fSpectators(1),
143      fSpecn(0),
144      fSpecp(0),
145      fDPMjet(0),
146      fNoGammas(0),
147      fLHC(0),
148      fPi0Decay(1),
149      fDecayAll(0),
150      fGenImpPar(0.),
151      fProcess(kDpmMb),
152      fTriggerParticle(0),
153      fTriggerEta(0.9),     
154      fTriggerMinPt(-1),  
155      fTriggerMaxPt(1000),  
156      fTriggerMultiplicity(0),
157      fTriggerMultiplicityEta(0),
158      fTriggerMultiplicityPtMin(0),
159      fkTuneForDiff(0),
160      fProcDiff(0),
161      fFragmentation(kFALSE),
162      fHeader(0x0)
163
164
165 {
166     // Dummy copy constructor
167     fEnergyCMS = 5500.;
168 }
169
170 //______________________________________________________________________________
171 AliGenDPMjet::~AliGenDPMjet()
172 {
173 // Destructor
174 }
175 //______________________________________________________________________________
176 void AliGenDPMjet::Init()
177 {
178 // Initialization
179     
180     if(fEnergyCMS>0. && fBeamEn<0.1) fBeamEn = fEnergyCMS/2;
181     SetMC(new TDPMjet(fProcess, fAProjectile, fZProjectile, fATarget, fZTarget, 
182                       fBeamEn,fEnergyCMS));
183
184     fDPMjet=(TDPMjet*) fMCEvGen;
185     if (fAProjectile == 1 && TMath::Abs(fZProjectile == 1)) fDPMjet->SetfIdp(1);
186     
187     // **** Flag to force central production
188     // fICentr=1. central production forced 
189     // fICentr<0 && fICentr>-100 -> bmin = fMinImpactParam, bmax = fMaxImpactParam        
190     // fICentr<-99 -> fraction of x-sec. = XSFRAC                 
191     // fICentr=-1. -> evaporation/fzc suppressed                  
192     // fICentr<-1. -> evaporation/fzc allowed                 
193     fDPMjet->SetfFCentr(fICentr);  
194     
195     fDPMjet->SetbRange(fMinImpactParam, fMaxImpactParam); 
196     fDPMjet->SetPi0Decay(fPi0Decay);
197     fDPMjet->SetDecayAll(fDecayAll);
198     fDPMjet->SetFragmentProd(fFragmentation);
199     
200 //
201 //  Initialize DPMjet  
202 //    
203     fDPMjet->Initialize();
204 }
205
206
207 //______________________________________________________________________________
208 void AliGenDPMjet::Generate()
209 {
210 // Generate one event
211
212   Double_t polar[3]    =   {0,0,0};
213   Double_t origin[3]   =   {0,0,0};
214   Double_t p[4]        =   {0};
215   Float_t tof;
216
217 //  converts from mm/c to s
218   const Float_t kconv = 0.001/2.999792458e8;
219   Int_t nt  = 0;
220   Int_t jev = 0;
221   Int_t kf, ks, imo;
222   kf = 0;
223   fTrials = 0;
224   fNprimaries = 0;
225   //  Set collision vertex position 
226   if (fVertexSmear == kPerEvent) Vertex();
227   
228   while(1)
229   {
230 //    Generate one event
231 // --------------------------------------------------------------------------
232       fSpecn = 0;  
233       fSpecp = 0;
234 // --------------------------------------------------------------------------
235       fDPMjet->GenerateEvent();
236       
237       fTrials++;
238
239       fDPMjet->ImportParticles(&fParticles,"All");      
240       if (fLHC) Boost();
241
242       // Temporaneo
243       fGenImpPar = fDPMjet->GetBImpac();
244
245       if(TMath::Abs(fXingAngleY) > 1.e-10) BeamCrossAngle();
246
247       Int_t np = fParticles.GetEntriesFast();
248       //
249       // Multiplicity Trigger
250       if (fTriggerMultiplicity > 0) {
251         Int_t multiplicity = 0;
252         for (Int_t i = 0; i < np; i++) {
253           TParticle *  iparticle = (TParticle *) fParticles.At(i);
254         
255           Int_t statusCode = iparticle->GetStatusCode();
256         
257           // Initial state particle
258           if (statusCode != 1)
259             continue;
260           // eta cut
261           if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
262             continue;
263           // pt cut
264           if (iparticle->Pt() < fTriggerMultiplicityPtMin) 
265             continue;
266           
267           TParticlePDG* pdgPart = iparticle->GetPDG();
268           if (pdgPart && pdgPart->Charge() == 0)
269             continue;
270           ++multiplicity;
271         }
272         //
273         //
274         if (multiplicity < fTriggerMultiplicity) continue;
275         Printf("Triggered on event with multiplicity of %d >= %d", multiplicity, fTriggerMultiplicity);
276       }    
277
278   //Trigger on the presence of a given particle in some phase space
279     if (fTriggerParticle) {
280         Bool_t triggered = kFALSE;
281             for (Long_t i = 0; i < np; i++) {
282                 TParticle *  iparticle = (TParticle *) fParticles.At(i);
283                 kf = CheckPDGCode(iparticle->GetPdgCode());
284                 if (kf != fTriggerParticle) continue;
285                 if (iparticle->Pt() == 0.) continue;
286                 if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
287                 if ( iparticle->Pt() > fTriggerMaxPt || iparticle->Pt() < fTriggerMinPt ) continue;
288                 triggered = kTRUE;
289                 break;
290             }
291       if (!triggered) continue; 
292     }
293
294       if(fkTuneForDiff && ( (TMath::Abs(fEnergyCMS - 900) < 1) || (TMath::Abs(fEnergyCMS - 2760) < 1) || (TMath::Abs(fEnergyCMS - 7000) < 1)) ) {
295         if(!CheckDiffraction() ) continue;
296       }
297
298
299       Int_t nc = 0;
300       if (np == 0) continue;
301
302       Int_t i;
303       Int_t* newPos     = new Int_t[np];
304       Int_t* pSelected  = new Int_t[np];
305
306       for (i = 0; i<np; i++) {
307           newPos[i]    = i;
308           pSelected[i] = 0;
309       }
310       
311 //      First select parent particles
312
313       for (i = 0; i<np; i++) {
314           TParticle *iparticle = (TParticle *) fParticles.At(i);
315
316 // Is this a parent particle ?
317
318           if (Stable(iparticle)) continue;
319
320           Bool_t  selected             =  kTRUE;
321           Bool_t  hasSelectedDaughters =  kFALSE;
322           
323           kf = iparticle->GetPdgCode();
324           if (kf == 92 || kf == 99999) continue;
325           ks = iparticle->GetStatusCode();
326 // No initial state partons
327           if (ks==21) continue;
328           if (!fSelectAll) selected = KinematicSelection(iparticle, 0) && 
329                                SelectFlavor(kf);
330
331           
332           hasSelectedDaughters = DaughtersSelection(iparticle);
333
334
335 // Put particle on the stack if it is either selected or 
336 // it is the mother of at least one seleted particle
337
338           if (selected || hasSelectedDaughters) {
339               nc++;
340               pSelected[i] = 1;
341           } // selected
342       } // particle loop parents
343
344 // Now select the final state particles
345
346
347       for (i=0; i<np; i++) {
348           TParticle *iparticle = (TParticle *) fParticles.At(i);
349
350 // Is this a final state particle ?
351
352           if (!Stable(iparticle)) continue;
353       
354           Bool_t  selected =  kTRUE;
355           kf = iparticle->GetPdgCode();
356           ks = iparticle->GetStatusCode();
357
358 // --------------------------------------------------------------------------
359 // Count spectator neutrons and protons (ks == 13, 14)
360           if(ks == 13 || ks == 14){
361               if(kf == kNeutron) fSpecn += 1;
362               if(kf == kProton)  fSpecp += 1;
363           }
364 // --------------------------------------------------------------------------
365
366           if (!fSelectAll) {
367               selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
368               if (!fSpectators && selected) selected = (ks == 13 || ks == 14);
369           }
370
371 // Put particle on the stack if selected
372
373           if (selected) {
374               nc++;
375               pSelected[i] = 1;
376           } // selected
377       } // particle loop final state
378
379 // Write particles to stack
380
381       for (i = 0; i<np; i++) {
382           TParticle *  iparticle = (TParticle *) fParticles.At(i);
383           Bool_t  hasMother   = (iparticle->GetFirstMother()>=0);
384           if (pSelected[i]) {
385               
386               kf   = iparticle->GetPdgCode();         
387               ks   = iparticle->GetStatusCode();              
388               
389               p[0] = iparticle->Px();
390               p[1] = iparticle->Py();
391               p[2] = iparticle->Pz();
392               p[3] = iparticle->Energy();
393               origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
394               origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
395               origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
396                     
397               tof = fTime + kconv*iparticle->T();
398               
399               imo = -1;
400               TParticle* mother = 0;
401               if (hasMother) {
402                   imo = iparticle->GetFirstMother();
403                   mother = (TParticle *) fParticles.At(imo);
404                   imo = (mother->GetPdgCode() != 92 && mother->GetPdgCode() != 99999) ? newPos[imo] : -1;
405               } // if has mother   
406
407
408               
409               Bool_t tFlag = (fTrackIt && (ks==1 || ks==-1));
410               //printf(" AliGemDPMJet->PushTrack: kf %d  ks %d  flag %d\n",kf,ks,tFlag);
411               if(kf>10000 && (ks==-1 || ks==1000 || ks==1001)) kf += 1000000000;
412               PushTrack(tFlag, imo, kf, 
413                         p[0], p[1], p[2], p[3], 
414                         origin[0], origin[1], origin[2], tof,
415                         polar[0], polar[1], polar[2],
416                         kPNoProcess, nt, 1., ks);
417               KeepTrack(nt);
418               fNprimaries++;
419               newPos[i] = nt;
420           } // if selected
421       } // particle loop
422       delete[] newPos;
423       delete[] pSelected;
424       if (nc>0) {
425           jev += nc;
426           if (jev >= fNpart || fNpart == -1) {
427               break;
428           }
429       }
430   } // event loop
431   MakeHeader();
432   SetHighWaterMark(nt);
433 }
434
435 //______________________________________________________________________________
436 Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle)
437 {
438 //
439 // Looks recursively if one of the daughters has been selected
440 //
441 //    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
442     Int_t imin = -1;
443     Int_t imax = -1;
444     Int_t i;
445     Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
446     Bool_t selected = kFALSE;
447     if (hasDaughters) {
448         imin = iparticle->GetFirstDaughter();
449         imax = iparticle->GetLastDaughter();       
450         for (i = imin; i <= imax; i++){
451             TParticle *  jparticle = (TParticle *) fParticles.At(i);    
452             Int_t ip = jparticle->GetPdgCode();
453             if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
454                 selected=kTRUE; break;
455             }
456             if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
457         }
458     } else {
459         return kFALSE;
460     }
461     return selected;
462 }
463
464
465
466 //______________________________________________________________________________
467 Bool_t AliGenDPMjet::SelectFlavor(Int_t pid)
468 {
469 // Select flavor of particle
470 // 0: all
471 // 4: charm and beauty
472 // 5: beauty
473     Bool_t res = 0;
474     
475     if (fFlavor == 0) {
476         res = kTRUE;
477     } else {
478         Int_t ifl = TMath::Abs(pid/100);
479         if (ifl > 10) ifl/=10;
480         res = (fFlavor == ifl);
481     }
482 //
483 //  This part if gamma writing is inhibited
484     if (fNoGammas) 
485         res = res && (pid != kGamma && pid != kPi0);
486 //
487     return res;
488 }
489
490 //______________________________________________________________________________
491 Bool_t AliGenDPMjet::Stable(TParticle*  particle)
492 {
493 // Return true for a stable particle
494 //
495     int st = particle->GetStatusCode();
496     if(st == 1 || st == -1) return kTRUE;
497     else return kFALSE;
498
499 }
500
501 //______________________________________________________________________________
502 void AliGenDPMjet::MakeHeader()
503 {
504 // Builds the event header, to be called after each event
505     fHeader.SetNProduced(fNprimaries);
506     fHeader.SetImpactParameter(fDPMjet->GetBImpac());
507     fHeader.SetTotalEnergy(fDPMjet->GetTotEnergy());
508     fHeader.SetParticipants(fDPMjet->GetProjParticipants(), 
509                             fDPMjet->GetTargParticipants());
510
511     fHeader.SetCollisions(DTGLCP.ncp, DTGLCP.nct, 
512         fDPMjet->GetProjWounded(),fDPMjet->GetTargWounded());
513                 
514     if(fProcDiff>0)  fHeader.SetProcessType(fProcDiff);
515     else fHeader.SetProcessType(fDPMjet->GetProcessCode());
516
517     // Bookkeeping for kinematic bias
518     fHeader.SetTrials(fTrials);
519     // Event Vertex
520     fHeader.SetPrimaryVertex(fVertex);
521     fHeader.SetInteractionTime(fTime);
522     fHeader.SetNDiffractive(POEVT1.nsd1, POEVT1.nsd2, POEVT1.ndd);
523 //    gAlice->SetGenEventHeader(fHeader);    
524     AddHeader(&fHeader);
525     fCollisionGeometry = &fHeader;
526 }
527
528 //______________________________________________________________________________
529 /*void AliGenDPMjet::AddHeader(AliGenEventHeader* fHeader)
530 {
531     // Add fHeader to container or runloader
532     if (fContainer) {
533         fContainer->AddHeader(fHeader);
534     } else {
535         AliRunLoader::Instance()->GetHeader()->SetGenEventHeader(fHeader);
536     }
537 }*/
538
539
540 //______________________________________________________________________________
541 AliGenDPMjet& AliGenDPMjet::operator=(const  AliGenDPMjet& /*rhs*/)
542 {
543 // Assignment operator
544     return *this;
545 }
546
547
548 //______________________________________________________________________________
549 void AliGenDPMjet::FinishRun()
550 {
551     // Print run statistics
552     fDPMjet->Dt_Dtuout();
553 }
554
555
556 //______________________________________________________________________________
557 Bool_t AliGenDPMjet::CheckDiffraction()
558 {
559
560   //  printf("AAA\n");
561
562    Int_t np = fParticles.GetEntriesFast();
563
564    Int_t iPart1=-1;
565    Int_t iPart2=-1;
566
567    Double_t y1 = 1e10;
568    Double_t y2 = -1e10;
569
570   const Int_t kNstable=20;
571   const Int_t pdgStable[20] = {
572     22,             // Photon
573     11,             // Electron
574     12,             // Electron Neutrino 
575     13,             // Muon 
576     14,             // Muon Neutrino
577     15,             // Tau 
578     16,             // Tau Neutrino
579     211,            // Pion
580     321,            // Kaon
581     311,            // K0
582     130,            // K0s
583     310,            // K0l
584     2212,           // Proton 
585     2112,           // Neutron
586     3122,           // Lambda_0
587     3112,           // Sigma Minus
588     3222,           // Sigma Plus
589     3312,           // Xsi Minus 
590     3322,           // Xsi0
591     3334            // Omega
592   };
593     
594      for (Int_t i = 0; i < np; i++) {
595         TParticle *  part = (TParticle *) fParticles.At(i);
596         
597         Int_t statusCode = part->GetStatusCode();
598         
599         // Initial state particle
600         if (statusCode != 1)
601           continue;
602
603         Int_t pdg = TMath::Abs(part->GetPdgCode());
604         Bool_t isStable = kFALSE;
605         for (Int_t i1 = 0; i1 < kNstable; i1++) {
606           if (pdg == pdgStable[i1]) {
607             isStable = kTRUE;
608             break;
609           }
610         }
611         if(!isStable) 
612           continue;
613
614         Double_t y = part->Y();
615
616         if (y < y1)
617           {
618             y1 = y;
619             iPart1 = i;
620           }
621         if (y > y2)
622         {
623           y2 = y;
624           iPart2 = i;
625         }
626      }
627
628      if(iPart1<0 || iPart2<0) return kFALSE;
629
630      y1=TMath::Abs(y1);
631      y2=TMath::Abs(y2);
632
633      TParticle *  part1 = (TParticle *) fParticles.At(iPart1);
634      TParticle *  part2 = (TParticle *) fParticles.At(iPart2);
635
636      Int_t pdg1 = part1->GetPdgCode();
637      Int_t pdg2 = part2->GetPdgCode();
638
639
640      Int_t iPart = -1;
641      if (pdg1 == 2212 && pdg2 == 2212)
642        {
643          if(y1 > y2) 
644            iPart = iPart1;
645          else if(y1 < y2) 
646            iPart = iPart2;
647          else {
648            iPart = iPart1;
649            if((AliDpmJetRndm::GetDpmJetRandom())->Uniform(0.,1.)>0.5) iPart = iPart2;
650          }
651        }
652      else if (pdg1 == 2212)
653        iPart = iPart1;
654      else if (pdg2 == 2212)
655        iPart = iPart2;
656
657
658
659
660
661      Double_t M=-1.;
662      if(iPart>0) {
663        TParticle *  part = (TParticle *) fParticles.At(iPart);
664        Double_t E= part->Energy();
665        Double_t P= part->P();
666        Double_t M2 = (fEnergyCMS-E-P)*(fEnergyCMS-E+P);
667        if(M2<0)  return kFALSE;
668        M= TMath::Sqrt(M2);
669      }
670
671      Double_t Mmin, Mmax, wSD, wDD, wND;
672
673      if(!GetWeightsDiffraction(M, Mmin, Mmax, wSD, wDD, wND)) return kFALSE;
674      if(M>-1 && M<Mmin) return kFALSE;
675      if(M>Mmax) M=-1;
676
677
678  Int_t procType=fDPMjet->GetProcessCode();//fPythia->GetMSTI(1);
679  if(procType== 4) return kFALSE;
680  Int_t proc0=2;
681  if(procType== 7) proc0=1;
682  if(procType== 5 || procType== 6) proc0=0;
683
684
685  Int_t proc=2;
686  if(M>0) proc=0;
687  else if(proc0==1) proc=1;
688
689  if(proc==1 && (AliDpmJetRndm::GetDpmJetRandom())->Uniform(0.,1.) > wDD) return kFALSE;
690  if(proc==2 && (AliDpmJetRndm::GetDpmJetRandom())->Uniform(0.,1.) > wND) return kFALSE;
691
692
693
694     if(proc!=0) {
695       if(proc0!=0) fProcDiff = procType;
696       else       fProcDiff = 1;
697       return kTRUE;
698     }
699
700     if(wSD<0)  AliError("wSD<0 ! \n");
701
702     if((AliDpmJetRndm::GetDpmJetRandom())->Uniform(0.,1.)> wSD) return kFALSE;
703
704     //    printf("iPart = %d\n", iPart);
705
706     if(iPart==iPart1) fProcDiff=5;
707     else if(iPart==iPart2) fProcDiff=6;
708     else {
709       printf("EROOR:  iPart!=iPart1 && iPart!=iPart2\n");
710
711     }
712
713     return kTRUE;
714 }
715
716 // -------------------------------------------------------
717 void AliGenDPMjet::SetIonPDGCodes()
718 {
719    // Defining PDG codes for the ions
720    AliIonPDGCodes *pdgcodes = new AliIonPDGCodes();
721    pdgcodes->AddParticlesToPdgDataBase();
722 }
723
724 // -------------------------------------------------------
725 Bool_t AliGenDPMjet::GetWeightsDiffraction(Double_t M, Double_t &Mmin, Double_t &Mmax, 
726                                                        Double_t &wSD, Double_t &wDD, Double_t &wND)
727 {
728
729   if(TMath::Abs(fEnergyCMS-900)<1 ){
730 const Int_t nbin=400;
731 Double_t bin[]={
732 1.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500, 
733 4.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300, 
734 7.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100, 
735 10.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900, 
736 13.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700, 
737 15.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500, 
738 18.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300, 
739 21.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100, 
740 24.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900, 
741 27.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700, 
742 30.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500, 
743 33.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300, 
744 36.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100, 
745 39.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900, 
746 42.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700, 
747 45.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500, 
748 48.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300, 
749 51.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100, 
750 54.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900, 
751 57.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700, 
752 60.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500, 
753 63.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300, 
754 66.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100, 
755 69.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900, 
756 72.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700, 
757 75.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500, 
758 78.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300, 
759 81.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100, 
760 84.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900, 
761 87.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700, 
762 90.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500, 
763 93.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300, 
764 96.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100, 
765 99.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900, 
766 102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700, 
767 105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500, 
768 108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300, 
769 111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100, 
770 114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900, 
771 117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700, 
772 120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500, 
773 123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300, 
774 126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100, 
775 129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900, 
776 132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700, 
777 135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500, 
778 138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300, 
779 141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100, 
780 144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900, 
781 147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700, 
782 150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500, 
783 153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300, 
784 156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100, 
785 159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900, 
786 162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700, 
787 165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500, 
788 168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300, 
789 171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100, 
790 174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900, 
791 177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700, 
792 180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500, 
793 183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300, 
794 186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100, 
795 189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900, 
796 192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700, 
797 195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500, 
798 198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
799 Double_t w[]={
800 1.000000, 0.389079, 0.326612, 0.356363, 0.310726, 0.264037, 
801 0.253810, 0.224655, 0.207990, 0.198149, 0.186803, 0.180209, 
802 0.178928, 0.161772, 0.162314, 0.158409, 0.148379, 0.143557, 
803 0.140357, 0.136272, 0.136957, 0.136606, 0.129376, 0.127671, 
804 0.128406, 0.132239, 0.119311, 0.130227, 0.130233, 0.123312, 
805 0.115257, 0.120105, 0.119733, 0.117116, 0.110096, 0.117509, 
806 0.109149, 0.114372, 0.100322, 0.106227, 0.108696, 0.110352, 
807 0.106748, 0.101475, 0.101837, 0.098585, 0.094433, 0.100148, 
808 0.096505, 0.100155, 0.103526, 0.098161, 0.093647, 0.100121, 
809 0.097426, 0.093414, 0.090241, 0.097228, 0.098098, 0.098340, 
810 0.096779, 0.099265, 0.095101, 0.090342, 0.097469, 0.090663, 
811 0.093422, 0.093103, 0.089648, 0.096923, 0.088079, 0.089360, 
812 0.097772, 0.092932, 0.093554, 0.085784, 0.096395, 0.096304, 
813 0.090183, 0.089255, 0.090265, 0.086262, 0.087044, 0.088965, 
814 0.085825, 0.095073, 0.085933, 0.091005, 0.091707, 0.092428, 
815 0.091689, 0.091224, 0.095256, 0.083407, 0.087983, 0.090320, 
816 0.081580, 0.088077, 0.083478, 0.091309, 0.083734, 0.089906, 
817 0.100155, 0.092728, 0.086542, 0.078091, 0.085261, 0.094302, 
818 0.078063, 0.090070, 0.086566, 0.095020, 0.082124, 0.084791, 
819 0.090624, 0.090236, 0.086344, 0.085706, 0.085913, 0.083107, 
820 0.092557, 0.081144, 0.091254, 0.096139, 0.079458, 0.088124, 
821 0.088777, 0.076652, 0.085168, 0.082326, 0.080435, 0.085022, 
822 0.085693, 0.093957, 0.088134, 0.082347, 0.082139, 0.090980, 
823 0.078397, 0.093595, 0.091562, 0.081789, 0.085026, 0.080868, 
824 0.083455, 0.084590, 0.084124, 0.081486, 0.087242, 0.080928, 
825 0.096392, 0.078324, 0.092093, 0.074267, 0.082925, 0.090717, 
826 0.081467, 0.080435, 0.078825, 0.075843, 0.088940, 0.081451, 
827 0.096047, 0.088102, 0.091710, 0.095208, 0.086160, 0.091305, 
828 0.085406, 0.094432, 0.078227, 0.099870, 0.094140, 0.080130, 
829 0.090707, 0.083268, 0.082222, 0.088767, 0.084477, 0.088069, 
830 0.089382, 0.086164, 0.089123, 0.092799, 0.076710, 0.090727, 
831 0.077097, 0.099905, 0.089733, 0.088101, 0.093705, 0.089215, 
832 0.085110, 0.086032, 0.083719, 0.079693, 0.088116, 0.090519, 
833 0.091150, 0.090855, 0.090547, 0.077773, 0.081914, 0.080864, 
834 0.082935, 0.092952, 0.079390, 0.080255, 0.091123, 0.091331, 
835 0.095160, 0.089343, 0.083353, 0.087445, 0.094036, 0.084719, 
836 0.099665, 0.085104, 0.100912, 0.085958, 0.083972, 0.098284, 
837 0.078318, 0.082042, 0.088007, 0.085469, 0.089984, 0.081181, 
838 0.098850, 0.086409, 0.089070, 0.081579, 0.089622, 0.089396, 
839 0.093922, 0.089472, 0.090806, 0.075034, 0.090346, 0.083871, 
840 0.086931, 0.089207, 0.094425, 0.087830, 0.079537, 0.097316, 
841 0.096513, 0.092264, 0.082211, 0.083841, 0.081861, 0.092591, 
842 0.077785, 0.079646, 0.093721, 0.090735, 0.086910, 0.075837, 
843 0.090729, 0.092800, 0.098704, 0.092441, 0.086404, 0.089344, 
844 0.083650, 0.082569, 0.085753, 0.091430, 0.085460, 0.095210, 
845 0.083226, 0.094882, 0.084856, 0.093350, 0.084579, 0.096744, 
846 0.099211, 0.085142, 0.085499, 0.083395, 0.088352, 0.110418, 
847 0.093788, 0.096140, 0.082758, 0.097145, 0.089170, 0.090720, 
848 0.084708, 0.095927, 0.090330, 0.074239, 0.089532, 0.090514, 
849 0.089823, 0.089709, 0.101840, 0.082676, 0.087157, 0.079221, 
850 0.096460, 0.108192, 0.088904, 0.084510, 0.096624, 0.099242, 
851 0.094470, 0.092473, 0.091745, 0.090439, 0.095316, 0.087963, 
852 0.096400, 0.089479, 0.094880, 0.085964, 0.092775, 0.083200, 
853 0.095133, 0.090079, 0.088828, 0.087600, 0.094123, 0.089135, 
854 0.082617, 0.085109, 0.114091, 0.096331, 0.108465, 0.080318, 
855 0.106576, 0.089671, 0.092023, 0.090722, 0.086603, 0.091788, 
856 0.098375, 0.082712, 0.102681, 0.098869, 0.089051, 0.109972, 
857 0.080440, 0.097860, 0.093074, 0.097028, 0.107826, 0.117152, 
858 0.090968, 0.096790, 0.096725, 0.094641, 0.089535, 0.092561, 
859 0.095828, 0.084320, 0.089942, 0.105476, 0.087495, 0.089805, 
860 0.092238, 0.094141, 0.104390, 0.082958, 0.097449, 0.099594, 
861 0.092640, 0.097332, 0.093223, 0.100183, 0.092511, 0.087035, 
862 0.088741, 0.097856, 0.116682, 0.091732, 0.097753, 0.094283, 
863 0.094235, 0.082016, 0.098370, 0.085676, 0.104529, 0.087319, 
864 0.090289, 0.105314, 0.103634, 0.101261, 0.092764, 0.098217, 
865 0.098939, 0.096071, 0.096071, 0.094027, 0.092713, 0.089542, 
866 0.112293, 0.112293, 0.089531, 0.087752};
867 wDD = 0.188264;
868 wND = 0.102252;
869 wSD = -1;
870
871     Mmin = bin[0];
872     Mmax = bin[nbin];
873     if(M<Mmin || M>Mmax) return kTRUE;
874
875     Int_t ibin=nbin-1;
876     for(Int_t i=1; i<=nbin; i++) 
877       if(M<=bin[i]) {
878         ibin=i-1;
879         //      printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
880         break;
881       }
882     wSD=w[ibin];
883     return kTRUE;
884
885   }
886   else if(TMath::Abs(fEnergyCMS-2760)<1 ){
887    
888 const Int_t nbin=400;
889 Double_t bin[]={
890 1.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500, 
891 4.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300, 
892 7.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100, 
893 10.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900, 
894 13.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700, 
895 15.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500, 
896 18.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300, 
897 21.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100, 
898 24.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900, 
899 27.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700, 
900 30.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500, 
901 33.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300, 
902 36.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100, 
903 39.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900, 
904 42.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700, 
905 45.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500, 
906 48.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300, 
907 51.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100, 
908 54.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900, 
909 57.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700, 
910 60.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500, 
911 63.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300, 
912 66.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100, 
913 69.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900, 
914 72.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700, 
915 75.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500, 
916 78.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300, 
917 81.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100, 
918 84.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900, 
919 87.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700, 
920 90.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500, 
921 93.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300, 
922 96.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100, 
923 99.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900, 
924 102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700, 
925 105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500, 
926 108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300, 
927 111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100, 
928 114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900, 
929 117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700, 
930 120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500, 
931 123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300, 
932 126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100, 
933 129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900, 
934 132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700, 
935 135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500, 
936 138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300, 
937 141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100, 
938 144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900, 
939 147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700, 
940 150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500, 
941 153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300, 
942 156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100, 
943 159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900, 
944 162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700, 
945 165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500, 
946 168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300, 
947 171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100, 
948 174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900, 
949 177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700, 
950 180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500, 
951 183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300, 
952 186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100, 
953 189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900, 
954 192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700, 
955 195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500, 
956 198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
957 Double_t w[]={
958 1.000000, 0.402402, 0.347976, 0.386866, 0.304413, 0.275746, 
959 0.256941, 0.250439, 0.228486, 0.219225, 0.203368, 0.196124, 
960 0.180551, 0.169230, 0.159021, 0.157136, 0.154743, 0.160522, 
961 0.153993, 0.143448, 0.148062, 0.136337, 0.130853, 0.127798, 
962 0.137576, 0.137530, 0.124175, 0.120225, 0.120774, 0.123584, 
963 0.121883, 0.114522, 0.118338, 0.119587, 0.106017, 0.125802, 
964 0.106256, 0.101585, 0.097737, 0.100817, 0.102279, 0.110935, 
965 0.101429, 0.107629, 0.105802, 0.095688, 0.096617, 0.100985, 
966 0.101521, 0.100450, 0.111836, 0.088828, 0.102264, 0.101012, 
967 0.097673, 0.097634, 0.091720, 0.096167, 0.084537, 0.104443, 
968 0.091966, 0.097204, 0.089456, 0.088273, 0.098220, 0.095151, 
969 0.088201, 0.085570, 0.088431, 0.077625, 0.088314, 0.084582, 
970 0.091545, 0.083774, 0.086183, 0.089195, 0.086255, 0.077167, 
971 0.096118, 0.089881, 0.086206, 0.086968, 0.079827, 0.082002, 
972 0.083081, 0.076587, 0.081820, 0.092382, 0.082964, 0.080901, 
973 0.093512, 0.082907, 0.086691, 0.093639, 0.089010, 0.082857, 
974 0.090795, 0.082403, 0.078602, 0.087284, 0.079638, 0.082905, 
975 0.077189, 0.084329, 0.075900, 0.082559, 0.084210, 0.070053, 
976 0.068453, 0.083369, 0.077659, 0.073286, 0.075396, 0.073946, 
977 0.083293, 0.076235, 0.082184, 0.080104, 0.073755, 0.077369, 
978 0.068281, 0.084593, 0.074923, 0.069467, 0.072094, 0.077973, 
979 0.067034, 0.080314, 0.087441, 0.071728, 0.072302, 0.076783, 
980 0.079931, 0.075557, 0.076318, 0.079029, 0.075126, 0.076859, 
981 0.080253, 0.074344, 0.083387, 0.087553, 0.084437, 0.076322, 
982 0.077748, 0.087559, 0.075649, 0.076615, 0.077098, 0.074559, 
983 0.075540, 0.085883, 0.079269, 0.079866, 0.068922, 0.069770, 
984 0.068024, 0.077783, 0.086620, 0.086769, 0.068433, 0.078059, 
985 0.084463, 0.072790, 0.073889, 0.079090, 0.081759, 0.069576, 
986 0.069160, 0.071695, 0.078569, 0.075727, 0.088055, 0.082395, 
987 0.079915, 0.075150, 0.073580, 0.069968, 0.071141, 0.065823, 
988 0.077270, 0.076229, 0.071735, 0.079271, 0.078209, 0.069503, 
989 0.064732, 0.076312, 0.088579, 0.087271, 0.080566, 0.073527, 
990 0.071376, 0.068785, 0.071812, 0.072680, 0.067252, 0.059543, 
991 0.072818, 0.071856, 0.083427, 0.080664, 0.077315, 0.066060, 
992 0.078772, 0.065442, 0.067255, 0.077676, 0.068306, 0.071578, 
993 0.059728, 0.073526, 0.078932, 0.070316, 0.067620, 0.067323, 
994 0.078316, 0.080366, 0.091194, 0.066739, 0.076238, 0.070382, 
995 0.077118, 0.073822, 0.072830, 0.062947, 0.077375, 0.076284, 
996 0.070951, 0.086127, 0.071074, 0.076621, 0.070502, 0.070897, 
997 0.058764, 0.065322, 0.073620, 0.062051, 0.074548, 0.083770, 
998 0.072502, 0.066823, 0.074192, 0.081570, 0.067340, 0.081896, 
999 0.080680, 0.072835, 0.069675, 0.078881, 0.072826, 0.069616, 
1000 0.082744, 0.069138, 0.073333, 0.080881, 0.069091, 0.064266, 
1001 0.062849, 0.072846, 0.067462, 0.067229, 0.071319, 0.071099, 
1002 0.069381, 0.073799, 0.076106, 0.075006, 0.073074, 0.092589, 
1003 0.074277, 0.061324, 0.069788, 0.070310, 0.063160, 0.073077, 
1004 0.075447, 0.081013, 0.074102, 0.076565, 0.060192, 0.058326, 
1005 0.071508, 0.072974, 0.065098, 0.060432, 0.062077, 0.075644, 
1006 0.075417, 0.066947, 0.066744, 0.065111, 0.075468, 0.078185, 
1007 0.073216, 0.066502, 0.079582, 0.065405, 0.069828, 0.072165, 
1008 0.071112, 0.064656, 0.080396, 0.066505, 0.063419, 0.066128, 
1009 0.080616, 0.071463, 0.064867, 0.074782, 0.063103, 0.068179, 
1010 0.064911, 0.073029, 0.075746, 0.062264, 0.099126, 0.067890, 
1011 0.074940, 0.070090, 0.073583, 0.056904, 0.063745, 0.063780, 
1012 0.076876, 0.066519, 0.067164, 0.068668, 0.075065, 0.070983, 
1013 0.075751, 0.066317, 0.077520, 0.073139, 0.073961, 0.085492, 
1014 0.074629, 0.059555, 0.081026, 0.059274, 0.055448, 0.084269, 
1015 0.066222, 0.078874, 0.062021, 0.075156, 0.090478, 0.066265, 
1016 0.067845, 0.076653, 0.086671, 0.082362, 0.067246, 0.075975, 
1017 0.072594, 0.068570, 0.095034, 0.065644, 0.070945, 0.066190, 
1018 0.075906, 0.055154, 0.075575, 0.058078, 0.071004, 0.070850, 
1019 0.064374, 0.058251, 0.069430, 0.072293, 0.065484, 0.084197, 
1020 0.090119, 0.091619, 0.067620, 0.079462, 0.063025, 0.068128, 
1021 0.056927, 0.076351, 0.073869, 0.061597, 0.083741, 0.063762, 
1022 0.064489, 0.074269, 0.068832, 0.058648, 0.069536, 0.074824, 
1023 0.081009, 0.073389, 0.076792, 0.084855, 0.075280, 0.061473, 
1024 0.063840, 0.062891, 0.071328, 0.065250};
1025 wDD = 0.077845;
1026 wND = 0.066355;
1027 wSD = -1;
1028
1029     Mmin = bin[0];
1030     Mmax = bin[nbin];
1031     if(M<Mmin || M>Mmax) return kTRUE;
1032
1033     Int_t ibin=nbin-1;
1034     for(Int_t i=1; i<=nbin; i++) 
1035       if(M<=bin[i]) {
1036         ibin=i-1;
1037         //      printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
1038         break;
1039       }
1040     wSD=w[ibin];
1041     return kTRUE;
1042    
1043   }
1044   else if(TMath::Abs(fEnergyCMS-7000)<1 ){
1045     
1046 const Int_t nbin=400;
1047 Double_t bin[]={
1048 1.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500, 
1049 4.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300, 
1050 7.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100, 
1051 10.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900, 
1052 13.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700, 
1053 15.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500, 
1054 18.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300, 
1055 21.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100, 
1056 24.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900, 
1057 27.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700, 
1058 30.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500, 
1059 33.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300, 
1060 36.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100, 
1061 39.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900, 
1062 42.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700, 
1063 45.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500, 
1064 48.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300, 
1065 51.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100, 
1066 54.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900, 
1067 57.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700, 
1068 60.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500, 
1069 63.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300, 
1070 66.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100, 
1071 69.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900, 
1072 72.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700, 
1073 75.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500, 
1074 78.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300, 
1075 81.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100, 
1076 84.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900, 
1077 87.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700, 
1078 90.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500, 
1079 93.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300, 
1080 96.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100, 
1081 99.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900, 
1082 102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700, 
1083 105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500, 
1084 108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300, 
1085 111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100, 
1086 114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900, 
1087 117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700, 
1088 120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500, 
1089 123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300, 
1090 126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100, 
1091 129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900, 
1092 132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700, 
1093 135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500, 
1094 138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300, 
1095 141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100, 
1096 144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900, 
1097 147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700, 
1098 150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500, 
1099 153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300, 
1100 156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100, 
1101 159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900, 
1102 162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700, 
1103 165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500, 
1104 168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300, 
1105 171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100, 
1106 174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900, 
1107 177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700, 
1108 180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500, 
1109 183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300, 
1110 186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100, 
1111 189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900, 
1112 192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700, 
1113 195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500, 
1114 198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
1115 Double_t w[]={
1116 1.000000, 0.526293, 0.446686, 0.437789, 0.366854, 0.333320, 
1117 0.291931, 0.266464, 0.253870, 0.248706, 0.232788, 0.220736, 
1118 0.209886, 0.202741, 0.188617, 0.182767, 0.178748, 0.169039, 
1119 0.175911, 0.169098, 0.171256, 0.146728, 0.144543, 0.159470, 
1120 0.153171, 0.151883, 0.144693, 0.136307, 0.140226, 0.135388, 
1121 0.141317, 0.151121, 0.131209, 0.144039, 0.124688, 0.128020, 
1122 0.119122, 0.121868, 0.122332, 0.119561, 0.115744, 0.102504, 
1123 0.114726, 0.109518, 0.115418, 0.114860, 0.110026, 0.107693, 
1124 0.103005, 0.115985, 0.108629, 0.105937, 0.101056, 0.101228, 
1125 0.113305, 0.110302, 0.104696, 0.107447, 0.099095, 0.107378, 
1126 0.103090, 0.111384, 0.090821, 0.109806, 0.093639, 0.096655, 
1127 0.110416, 0.104446, 0.098530, 0.095105, 0.106641, 0.101804, 
1128 0.091798, 0.094775, 0.104130, 0.088436, 0.119692, 0.099160, 
1129 0.086163, 0.089250, 0.101898, 0.091786, 0.087346, 0.095312, 
1130 0.107186, 0.085671, 0.093283, 0.091992, 0.114654, 0.086172, 
1131 0.084559, 0.097862, 0.079417, 0.094737, 0.089678, 0.105265, 
1132 0.099092, 0.087901, 0.086828, 0.084756, 0.077769, 0.084254, 
1133 0.083262, 0.092935, 0.088858, 0.093377, 0.083569, 0.084771, 
1134 0.084957, 0.084033, 0.095007, 0.076690, 0.087726, 0.097913, 
1135 0.084119, 0.076361, 0.084172, 0.091068, 0.094564, 0.091406, 
1136 0.088405, 0.090355, 0.085388, 0.084586, 0.085752, 0.084298, 
1137 0.091294, 0.081491, 0.075024, 0.082543, 0.088048, 0.074621, 
1138 0.078514, 0.077249, 0.078378, 0.092945, 0.075859, 0.077602, 
1139 0.074592, 0.078181, 0.081527, 0.080865, 0.078854, 0.078237, 
1140 0.085506, 0.089823, 0.072738, 0.078430, 0.077162, 0.081495, 
1141 0.090878, 0.099417, 0.100966, 0.087960, 0.085556, 0.081661, 
1142 0.078066, 0.089059, 0.073564, 0.092704, 0.073148, 0.098057, 
1143 0.067906, 0.079452, 0.090567, 0.082454, 0.077851, 0.079694, 
1144 0.091272, 0.078628, 0.096906, 0.098779, 0.088906, 0.070174, 
1145 0.083822, 0.084241, 0.093237, 0.071062, 0.075771, 0.096405, 
1146 0.098441, 0.086007, 0.069599, 0.078400, 0.083481, 0.081054, 
1147 0.085552, 0.069582, 0.071336, 0.073207, 0.087913, 0.080125, 
1148 0.075189, 0.067217, 0.073509, 0.099694, 0.080781, 0.073943, 
1149 0.062822, 0.076611, 0.087672, 0.077632, 0.063959, 0.077754, 
1150 0.084651, 0.097348, 0.063909, 0.055053, 0.087616, 0.074428, 
1151 0.101165, 0.078446, 0.070683, 0.071162, 0.091516, 0.069641, 
1152 0.070969, 0.069778, 0.083318, 0.066703, 0.074027, 0.067589, 
1153 0.070620, 0.081307, 0.094406, 0.076188, 0.061663, 0.077561, 
1154 0.076159, 0.071851, 0.074417, 0.076136, 0.069073, 0.075450, 
1155 0.070218, 0.089709, 0.079974, 0.082077, 0.076979, 0.070497, 
1156 0.076296, 0.087347, 0.073508, 0.088073, 0.067186, 0.059898, 
1157 0.085495, 0.100152, 0.079468, 0.093783, 0.082680, 0.077272, 
1158 0.081995, 0.068240, 0.071881, 0.073737, 0.072346, 0.079034, 
1159 0.078721, 0.067518, 0.068196, 0.081738, 0.082814, 0.082480, 
1160 0.069072, 0.066853, 0.081523, 0.073588, 0.082308, 0.091641, 
1161 0.072747, 0.084317, 0.072190, 0.079372, 0.067424, 0.077450, 
1162 0.062343, 0.090708, 0.065470, 0.086588, 0.071172, 0.066705, 
1163 0.070670, 0.070407, 0.096111, 0.066738, 0.081249, 0.072882, 
1164 0.075144, 0.060331, 0.074589, 0.076968, 0.085913, 0.072561, 
1165 0.064645, 0.078742, 0.075670, 0.065984, 0.080932, 0.069898, 
1166 0.065303, 0.096856, 0.057690, 0.065720, 0.066545, 0.068521, 
1167 0.068278, 0.069245, 0.086643, 0.063401, 0.070933, 0.070752, 
1168 0.066978, 0.058891, 0.070073, 0.078031, 0.082691, 0.101414, 
1169 0.075814, 0.072790, 0.057622, 0.093002, 0.084660, 0.079216, 
1170 0.070371, 0.070141, 0.076944, 0.067285, 0.078016, 0.077807, 
1171 0.066668, 0.066459, 0.059962, 0.062774, 0.083450, 0.064554, 
1172 0.067887, 0.064165, 0.072782, 0.067285, 0.052710, 0.096824, 
1173 0.071931, 0.064190, 0.074442, 0.082647, 0.055797, 0.078632, 
1174 0.061116, 0.063092, 0.049131, 0.074517, 0.069915, 0.079021, 
1175 0.088656, 0.101045, 0.090432, 0.076460, 0.067451, 0.071381, 
1176 0.069790, 0.077330, 0.077115, 0.084403, 0.065138, 0.056436, 
1177 0.088024, 0.069893, 0.055985, 0.089655, 0.062911, 0.075311, 
1178 0.086679, 0.093645, 0.068632, 0.064498, 0.057677, 0.081391, 
1179 0.062781, 0.075467, 0.061314, 0.073394, 0.084462, 0.068470, 
1180 0.071267, 0.060556, 0.072487, 0.063785, 0.079164, 0.070406, 
1181 0.073394, 0.063168, 0.066968, 0.064125, 0.063954, 0.072504, 
1182 0.058948, 0.057740, 0.083383, 0.067262};
1183 wDD = 0.098024;
1184 wND = 0.053452;
1185 wSD = -1;
1186
1187     Mmin = bin[0];
1188     Mmax = bin[nbin];
1189     if(M<Mmin || M>Mmax) return kTRUE;
1190
1191     Int_t ibin=nbin-1;
1192     for(Int_t i=1; i<=nbin; i++) 
1193       if(M<=bin[i]) {
1194         ibin=i-1;
1195         //      printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
1196         break;
1197       }
1198     wSD=w[ibin];
1199     return kTRUE;
1200    
1201   }
1202
1203   return kFALSE;
1204 }
1205
1206
1207 //______________________________________________________________________________