Changes related to the initialization of random numbers generators. Now one can use...
[u/mrichter/AliRoot.git] / EVGEN / AliGenParam.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 $Log$
18 Revision 1.42  2003/03/15 19:48:01  morsch
19 AliDecayerPythia replaced by AliDecayer
20
21 Revision 1.41  2003/01/09 17:38:47  morsch
22 Draw() method added.
23
24 Revision 1.40  2002/10/14 14:55:35  hristov
25 Merging the VirtualMC branch to the main development branch (HEAD)
26
27 Revision 1.36.6.3  2002/10/10 16:40:08  hristov
28 Updating VirtualMC to v3-09-02
29
30 Revision 1.39  2002/09/16 08:21:16  morsch
31 Use TDatabasePDG::Instance();
32
33 Revision 1.38  2002/05/30 14:59:12  morsch
34 Check geometrical acceptance. (G. Martinez)
35
36 Revision 1.37  2002/04/17 10:20:44  morsch
37 Coding Rule violations corrected.
38
39 Revision 1.36  2002/02/08 16:50:50  morsch
40 Add name and title in constructor.
41
42 Revision 1.35  2002/01/21 10:02:40  morsch
43 ty is Double_t
44 Abort if too high rapidity causes numerical paroblem. User has to specify
45 meaningful y-range.
46
47 Revision 1.34  2001/11/27 13:13:07  morsch
48 Maximum lifetime for long-lived particles to be put on the stack is parameter.
49 It can be set via SetMaximumLifetime(..).
50
51 Revision 1.33  2001/10/21 18:35:56  hristov
52 Several pointers were set to zero in the default constructors to avoid memory management problems
53
54 Revision 1.32  2001/07/27 17:09:36  morsch
55 Use local SetTrack, KeepTrack and SetHighWaterMark methods
56 to delegate either to local stack or to stack owned by AliRun.
57 (Piotr Skowronski, A.M.)
58
59 Revision 1.31  2001/07/13 10:58:54  morsch
60 - Some coded moved to AliGenMC
61 - Improved handling of secondary vertices.
62
63 Revision 1.30  2001/06/15 07:55:04  morsch
64 Put only first generation decay products on the stack.
65
66 Revision 1.29  2001/03/27 10:58:41  morsch
67 Initialize decayer before generation. Important if run inside cocktail.
68
69 Revision 1.28  2001/03/09 13:01:41  morsch
70 - enum constants for paramterisation type (particle family) moved to AliGen*lib.h
71 - use AliGenGSIlib::kUpsilon, AliGenPHOSlib::kEtaPrime to access the constants
72
73 Revision 1.27  2001/02/02 15:21:10  morsch
74 Set high water mark after last particle.
75 Use Vertex() method for Vertex.
76
77 Revision 1.26  2000/12/21 16:24:06  morsch
78 Coding convention clean-up
79
80 Revision 1.25  2000/11/30 07:12:50  alibrary
81 Introducing new Rndm and QA classes
82
83 Revision 1.24  2000/10/18 19:11:27  hristov
84 Division by zero fixed
85
86 Revision 1.23  2000/10/02 21:28:06  fca
87 Removal of useless dependecies via forward declarations
88
89 Revision 1.22  2000/09/12 14:14:55  morsch
90 Call fDecayer->ForceDecay() at the beginning of Generate().
91
92 Revision 1.21  2000/09/08 15:39:01  morsch
93 Handle the case fForceDecay=all during the generation, i.e. select all secondaries.
94
95 Revision 1.20  2000/09/06 14:35:44  morsch
96 Use AliDecayerPythia for particle decays.
97
98 Revision 1.19  2000/07/11 18:24:56  fca
99 Coding convention corrections + few minor bug fixes
100
101 Revision 1.18  2000/06/29 21:08:27  morsch
102 All paramatrisation libraries derive from the pure virtual base class AliGenLib.
103 This allows to pass a pointer to a library directly to AliGenParam and avoids the
104 use of function pointers in Config.C.
105
106 Revision 1.17  2000/06/09 20:33:30  morsch
107 All coding rule violations except RS3 corrected
108
109 Revision 1.16  2000/05/02 07:51:31  morsch
110 - Control precision of pT sampling TF1::SetNpx(..)
111 - Correct initialisation of child-cuts in all constructors.
112 - Most coding rule violations corrected.
113
114 Revision 1.15  2000/04/03 15:42:12  morsch
115 Cuts on primary particles are separated from those on the decay products. Methods
116 SetChildMomentumRange, SetChildPtRange, SetChildPhiRange, SetChildThetaRange added.
117
118 Revision 1.14  1999/11/09 07:38:48  fca
119 Changes for compatibility with version 2.23 of ROOT
120
121 Revision 1.13  1999/11/04 11:30:31  fca
122 Correct the logics for SetForceDecay
123
124 Revision 1.12  1999/11/03 17:43:20  fca
125 New version from G.Martinez & A.Morsch
126
127 Revision 1.11  1999/09/29 09:24:14  fca
128 Introduction of the Copyright and cvs Log
129
130 */
131
132
133
134 // Class to generate particles from using paramtrized pT and y distributions.
135 // Distributions are obtained from pointer to object of type AliGenLib.
136 // (For example AliGenMUONlib)
137 // Decays are performed using Pythia.
138 // andreas.morsch@cern.ch
139
140 #include "AliGenParam.h"
141 #include "AliDecayer.h"
142 #include "AliGenMUONlib.h"
143 #include "AliRun.h"
144 #include <TParticle.h>
145 #include <TParticlePDG.h>
146 #include <TDatabasePDG.h>
147 #include <TLorentzVector.h>
148
149 #include <TF1.h>
150 #include <TCanvas.h>
151 #include <TH1.h>
152
153 ClassImp(AliGenParam)
154
155 //------------------------------------------------------------
156
157   //Begin_Html
158   /*
159     <img src="picts/AliGenParam.gif">
160   */
161   //End_Html
162
163 //____________________________________________________________
164 //____________________________________________________________
165 AliGenParam::AliGenParam()
166 {
167 // Deafault constructor
168     fPtPara = 0;
169     fYPara  = 0;
170     fParam  = 0;
171     fAnalog = kAnalog;
172     SetDeltaPt();
173     fDecayer = 0;
174
175
176 }
177
178 AliGenParam::AliGenParam(Int_t npart, AliGenLib * Library,  Int_t param, char* tname):AliGenMC(npart)
179 {
180 // Constructor using number of particles parameterisation id and library
181     fName = "Param";
182     fTitle= "Particle Generator using pT and y parameterisation";
183     
184     fPtParaFunc = Library->GetPt(param, tname);
185     fYParaFunc  = Library->GetY (param, tname);
186     fIpParaFunc = Library->GetIp(param, tname);
187     
188     fPtPara = 0;
189     fYPara  = 0;
190     fParam  = param;
191     fAnalog = kAnalog;
192     SetForceDecay();
193     SetDeltaPt(); 
194 }
195
196 //____________________________________________________________
197
198 AliGenParam::AliGenParam(Int_t npart, Int_t param, char* tname):AliGenMC(npart)
199 {
200 // Constructor using parameterisation id and number of particles
201 //
202     fName = "Param";
203     fTitle= "Particle Generator using pT and y parameterisation";
204       
205     AliGenLib* pLibrary = new AliGenMUONlib();
206  
207     fPtParaFunc = pLibrary->GetPt(param, tname);
208     fYParaFunc  = pLibrary->GetY (param, tname);
209     fIpParaFunc = pLibrary->GetIp(param, tname);
210     
211     fPtPara = 0;
212     fYPara  = 0;
213     fParam  = param;
214     fAnalog = kAnalog;
215     fChildSelect.Set(5);
216     for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
217     SetForceDecay();
218     SetCutOnChild();
219     SetChildMomentumRange();
220     SetChildPtRange();
221     SetChildPhiRange();
222     SetChildThetaRange(); 
223     SetDeltaPt(); 
224 }
225
226 AliGenParam::AliGenParam(Int_t npart, Int_t param,
227                          Double_t (*PtPara) (Double_t*, Double_t*),
228                          Double_t (*YPara ) (Double_t* ,Double_t*),
229                          Int_t    (*IpPara) (TRandom *))                 
230     :AliGenMC(npart)
231 {
232 // Constructor
233 // Gines Martinez 1/10/99 
234     fName = "Param";
235     fTitle= "Particle Generator using pT and y parameterisation";
236
237     fPtParaFunc = PtPara; 
238     fYParaFunc  = YPara;  
239     fIpParaFunc = IpPara;
240 //  
241     fPtPara = 0;
242     fYPara  = 0;
243     fParam  = param;
244     fAnalog = kAnalog;
245     fChildSelect.Set(5);
246     for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
247     SetForceDecay();
248     SetCutOnChild();
249     SetChildMomentumRange();
250     SetChildPtRange();
251     SetChildPhiRange();
252     SetChildThetaRange();  
253     SetDeltaPt();
254 }
255
256
257 AliGenParam::AliGenParam(const AliGenParam & Param)
258 {
259 // copy constructor
260     Param.Copy(*this);
261 }
262
263 //____________________________________________________________
264 AliGenParam::~AliGenParam()
265 {
266 // Destructor
267     delete  fPtPara;
268     delete  fYPara;
269 }
270
271 //____________________________________________________________
272 void AliGenParam::Init()
273 {
274 // Initialisation
275
276     if (gMC) fDecayer = gMC->GetDecayer();
277   //Begin_Html
278   /*
279     <img src="picts/AliGenParam.gif">
280   */
281   //End_Html
282  
283     fPtPara = new TF1("Pt-Parametrization",fPtParaFunc,fPtMin,fPtMax,0);
284 //  Set representation precision to 10 MeV
285     Int_t npx= Int_t((fPtMax-fPtMin)/fDeltaPt);
286     
287     fPtPara->SetNpx(npx);
288     
289     fYPara  = new TF1("Y -Parametrization",fYParaFunc,fYMin,fYMax,0);
290     TF1* ptPara = new TF1("Pt-Parametrization",fPtParaFunc,0,15,0);
291     TF1* yPara  = new TF1("Y -Parametrization",fYParaFunc,-6,6,0);
292
293 //
294 // dN/dy| y=0
295     Double_t y1=0;
296     Double_t y2=0;
297     
298     fdNdy0=fYParaFunc(&y1,&y2);
299 //
300 // Integral over generation region
301     Float_t intYS  = yPara ->Integral(fYMin, fYMax);
302     Float_t intPt0 = ptPara->Integral(0,15);
303     Float_t intPtS = ptPara->Integral(fPtMin,fPtMax);
304     Float_t phiWgt=(fPhiMax-fPhiMin)/2./TMath::Pi();
305     if (fAnalog == kAnalog) {
306         fYWgt  = intYS/fdNdy0;
307         fPtWgt = intPtS/intPt0;
308         fParentWeight = fYWgt*fPtWgt*phiWgt/fNpart;
309     } else {
310         fYWgt = intYS/fdNdy0;
311         fPtWgt = (fPtMax-fPtMin)/intPt0;
312         fParentWeight = fYWgt*fPtWgt*phiWgt/fNpart;
313     }
314 //
315 // particle decay related initialization
316     fDecayer->SetForceDecay(fForceDecay);
317     fDecayer->Init();
318
319 //
320     AliGenMC::Init();
321 }
322
323 //____________________________________________________________
324 void AliGenParam::Generate()
325 {
326 //
327 // Generate 'npart' of light and heavy mesons (J/Psi, upsilon or phi, Pion,
328 // Kaons, Etas, Omegas) and Baryons (proton, antiprotons, neutrons and 
329 // antineutrons in the the desired theta, phi and momentum windows; 
330 // Gaussian smearing on the vertex is done if selected. 
331 // The decay of heavy mesons is done using lujet, 
332 //    and the childern particle are tracked by GEANT
333 // However, light mesons are directly tracked by GEANT 
334 // setting fForceDecay = nodecay (SetForceDecay(nodecay)) 
335 //
336 //
337 //  Reinitialize decayer
338     fDecayer->Init();
339 //
340   Float_t polar[3]= {0,0,0};  // Polarisation of the parent particle (for GEANT tracking)
341   Float_t origin0[3];         // Origin of the generated parent particle (for GEANT tracking)
342   Float_t pt, pl, ptot;       // Transverse, logitudinal and total momenta of the parent particle
343   Float_t phi, theta;         // Phi and theta spherical angles of the parent particle momentum
344   Float_t p[3], pc[3], 
345           och[3];             // Momentum, polarisation and origin of the children particles from lujet
346   Double_t ty, xmt;
347   Int_t nt, i, j;
348   Float_t  wgtp, wgtch;
349   Double_t dummy;
350   static TClonesArray *particles;
351   //
352   if(!particles) particles = new TClonesArray("TParticle",1000);
353   
354   TDatabasePDG *pDataBase = TDatabasePDG::Instance();
355   //
356   Float_t random[6];
357  
358 // Calculating vertex position per event
359   for (j=0;j<3;j++) origin0[j]=fOrigin[j];
360   if(fVertexSmear==kPerEvent) {
361       Vertex();
362       for (j=0;j<3;j++) origin0[j]=fVertex[j];
363   }
364   
365   Int_t ipa=0;
366   
367 // Generating fNpart particles
368   while (ipa<fNpart) {
369       while(1) {
370 //
371 // particle type
372           Int_t iPart = fIpParaFunc(fRandom);
373           fChildWeight=(fDecayer->GetPartialBranchingRatio(iPart))*fParentWeight;          
374           TParticlePDG *particle = pDataBase->GetParticle(iPart);
375           Float_t am = particle->Mass();
376           
377           Rndm(random,2);
378 //
379 // phi
380           phi=fPhiMin+random[0]*(fPhiMax-fPhiMin);
381 //
382 // y
383           ty = TMath::TanH(fYPara->GetRandom());
384 //
385 // pT
386           if (fAnalog == kAnalog) {
387               pt=fPtPara->GetRandom();
388               wgtp=fParentWeight;
389               wgtch=fChildWeight;
390           } else {
391               pt=fPtMin+random[1]*(fPtMax-fPtMin);
392               Double_t ptd=pt;
393               wgtp=fParentWeight*fPtParaFunc(& ptd, &dummy);
394               wgtch=fChildWeight*fPtParaFunc(& ptd, &dummy);
395           }
396           xmt=sqrt(pt*pt+am*am);
397           if (TMath::Abs(ty)==1.) {
398               ty=0.;
399               Fatal("AliGenParam", 
400                     "Division by 0: Please check you rapidity range !");
401           }
402           
403           pl=xmt*ty/sqrt(1.-ty*ty);
404           theta=TMath::ATan2(pt,pl);
405 // Cut on theta
406           if(theta<fThetaMin || theta>fThetaMax) continue;
407           ptot=TMath::Sqrt(pt*pt+pl*pl);
408 // Cut on momentum
409           if(ptot<fPMin || ptot>fPMax) continue;
410 //
411           p[0]=pt*TMath::Cos(phi);
412           p[1]=pt*TMath::Sin(phi);
413           p[2]=pl;
414           if(fVertexSmear==kPerTrack) {
415               Rndm(random,6);
416               for (j=0;j<3;j++) {
417                   origin0[j]=
418                       fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
419                       TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
420               }
421           }
422           
423 // Looking at fForceDecay : 
424 // if fForceDecay != none Primary particle decays using 
425 // AliPythia and children are tracked by GEANT
426 //
427 // if fForceDecay == none Primary particle is tracked by GEANT 
428 // (In the latest, make sure that GEANT actually does all the decays you want)    
429 //
430
431           if (fForceDecay != kNoDecay) {
432 // Using lujet to decay particle
433               Float_t energy=TMath::Sqrt(ptot*ptot+am*am);
434               TLorentzVector pmom(p[0], p[1], p[2], energy);
435               fDecayer->Decay(iPart,&pmom);
436 //
437 // select decay particles
438               Int_t np=fDecayer->ImportParticles(particles);
439
440               //  Selecting  GeometryAcceptance for particles fPdgCodeParticleforAcceptanceCut;
441               if (fGeometryAcceptance) 
442                 if (!CheckAcceptanceGeometry(np,particles)) continue;
443               Int_t ncsel=0;
444               Int_t* pFlag      = new Int_t[np];
445               Int_t* pParent    = new Int_t[np];
446               Int_t* pSelected  = new Int_t[np];
447               Int_t* trackIt    = new Int_t[np];
448
449               for (i=0; i<np; i++) {
450                   pFlag[i]     =  0;
451                   pSelected[i] =  0;
452                   pParent[i]   = -1;
453               }
454               
455               if (np >1) {
456                   TParticle* iparticle =  (TParticle *) particles->At(0);
457                   Int_t ipF, ipL;
458                   for (i = 1; i<np ; i++) {
459                       trackIt[i] = 1;
460                       iparticle = (TParticle *) particles->At(i);
461                       Int_t kf = iparticle->GetPdgCode();
462                       Int_t ks = iparticle->GetStatusCode();
463 // flagged particle
464
465                       if (pFlag[i] == 1) {
466                           ipF = iparticle->GetFirstDaughter();
467                           ipL = iparticle->GetLastDaughter();   
468                           if (ipF > 0) for (j=ipF-1; j<ipL; j++) pFlag[j]=1;
469                           continue;
470                       }
471
472 // flag decay products of particles with long life-time (c tau > .3 mum)                      
473                       
474                       if (ks != 1) { 
475 //                        TParticlePDG *particle = pDataBase->GetParticle(kf);
476                           
477                           Double_t lifeTime = fDecayer->GetLifetime(kf);
478 //                        Double_t mass     = particle->Mass();
479 //                        Double_t width    = particle->Width();
480                           if (lifeTime > (Double_t) fMaxLifeTime) {
481                               ipF = iparticle->GetFirstDaughter();
482                               ipL = iparticle->GetLastDaughter();       
483                               if (ipF > 0) for (j=ipF-1; j<ipL; j++) pFlag[j]=1;
484                           } else{
485                               trackIt[i]     = 0;
486                               pSelected[i]   = 1;
487                           }
488                       } // ks==1 ?
489 //
490 // children
491                       
492                       if (ChildSelected(TMath::Abs(kf)) || fForceDecay == kAll && trackIt[i])
493                       {
494                           if (fCutOnChild) {
495                               pc[0]=iparticle->Px();
496                               pc[1]=iparticle->Py();
497                               pc[2]=iparticle->Pz();
498                               Bool_t  childok = KinematicSelection(iparticle, 1);
499                               if(childok) {
500                                   pSelected[i]  = 1;
501                                   ncsel++;
502                               } else {
503                                   ncsel=-1;
504                                   break;
505                               } // child kine cuts
506                           } else {
507                               pSelected[i]  = 1;
508                               ncsel++;
509                           } // if child selection
510                       } // select muon
511                   } // decay particle loop
512               } // if decay products
513               
514               Int_t iparent;
515               if ((fCutOnChild && ncsel >0) || !fCutOnChild){
516                   ipa++;
517 //
518 // Parent
519                   SetTrack(0, -1, iPart, p, origin0, polar, 0, kPPrimary, nt, wgtp);
520                   pParent[0] = nt;
521                   KeepTrack(nt); 
522 //
523 // Decay Products
524 //                
525                   for (i = 1; i < np; i++) {
526                       if (pSelected[i]) {
527                           TParticle* iparticle = (TParticle *) particles->At(i);
528                           Int_t kf  = iparticle->GetPdgCode();
529                           Int_t ipa = iparticle->GetFirstMother()-1;
530                           
531                           och[0] = origin0[0]+iparticle->Vx()/10;
532                           och[1] = origin0[1]+iparticle->Vy()/10;
533                           och[2] = origin0[2]+iparticle->Vz()/10;
534                           pc[0]  = iparticle->Px();
535                           pc[1]  = iparticle->Py();
536                           pc[2]  = iparticle->Pz();
537                           
538                           if (ipa > -1) {
539                               iparent = pParent[ipa];
540                           } else {
541                               iparent = -1;
542                           }
543                          
544                           SetTrack(fTrackIt*trackIt[i], iparent, kf,
545                                            pc, och, polar,
546                                            0, kPDecay, nt, wgtch);
547                           pParent[i] = nt;
548                           KeepTrack(nt); 
549                       } // Selected
550                   } // Particle loop 
551               }  // Decays by Lujet
552               particles->Clear();
553               if (pFlag)      delete[] pFlag;
554               if (pParent)    delete[] pParent;
555               if (pSelected)  delete[] pSelected;          
556               if (trackIt)    delete[] trackIt;
557           } // kinematic selection
558           else  // nodecay option, so parent will be tracked by GEANT (pions, kaons, eta, omegas, baryons)
559           {
560             gAlice->
561                 SetTrack(fTrackIt,-1,iPart,p,origin0,polar,0,kPPrimary,nt,wgtp);
562             ipa++; 
563           }
564           break;
565     } // while
566   } // event loop
567   SetHighWaterMark(nt);
568 }
569
570 void AliGenParam::Draw()
571 {
572     //
573     // Draw the pT and y Distributions
574     //
575      TCanvas *c0 = new TCanvas("c0","Canvas 0",400,10,600,700);
576      c0->Divide(2,1);
577      c0->cd(1);
578      fPtPara->Draw();
579      fPtPara->GetHistogram()->SetXTitle("p_{T} (GeV)");     
580      c0->cd(2);
581      fYPara->Draw();
582      fYPara->GetHistogram()->SetXTitle("y");     
583 }
584
585 AliGenParam& AliGenParam::operator=(const  AliGenParam& rhs)
586 {
587 // Assignment operator
588     return *this;
589 }
590
591
592