Draw() method added.
[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.40  2002/10/14 14:55:35  hristov
19 Merging the VirtualMC branch to the main development branch (HEAD)
20
21 Revision 1.36.6.3  2002/10/10 16:40:08  hristov
22 Updating VirtualMC to v3-09-02
23
24 Revision 1.39  2002/09/16 08:21:16  morsch
25 Use TDatabasePDG::Instance();
26
27 Revision 1.38  2002/05/30 14:59:12  morsch
28 Check geometrical acceptance. (G. Martinez)
29
30 Revision 1.37  2002/04/17 10:20:44  morsch
31 Coding Rule violations corrected.
32
33 Revision 1.36  2002/02/08 16:50:50  morsch
34 Add name and title in constructor.
35
36 Revision 1.35  2002/01/21 10:02:40  morsch
37 ty is Double_t
38 Abort if too high rapidity causes numerical paroblem. User has to specify
39 meaningful y-range.
40
41 Revision 1.34  2001/11/27 13:13:07  morsch
42 Maximum lifetime for long-lived particles to be put on the stack is parameter.
43 It can be set via SetMaximumLifetime(..).
44
45 Revision 1.33  2001/10/21 18:35:56  hristov
46 Several pointers were set to zero in the default constructors to avoid memory management problems
47
48 Revision 1.32  2001/07/27 17:09:36  morsch
49 Use local SetTrack, KeepTrack and SetHighWaterMark methods
50 to delegate either to local stack or to stack owned by AliRun.
51 (Piotr Skowronski, A.M.)
52
53 Revision 1.31  2001/07/13 10:58:54  morsch
54 - Some coded moved to AliGenMC
55 - Improved handling of secondary vertices.
56
57 Revision 1.30  2001/06/15 07:55:04  morsch
58 Put only first generation decay products on the stack.
59
60 Revision 1.29  2001/03/27 10:58:41  morsch
61 Initialize decayer before generation. Important if run inside cocktail.
62
63 Revision 1.28  2001/03/09 13:01:41  morsch
64 - enum constants for paramterisation type (particle family) moved to AliGen*lib.h
65 - use AliGenGSIlib::kUpsilon, AliGenPHOSlib::kEtaPrime to access the constants
66
67 Revision 1.27  2001/02/02 15:21:10  morsch
68 Set high water mark after last particle.
69 Use Vertex() method for Vertex.
70
71 Revision 1.26  2000/12/21 16:24:06  morsch
72 Coding convention clean-up
73
74 Revision 1.25  2000/11/30 07:12:50  alibrary
75 Introducing new Rndm and QA classes
76
77 Revision 1.24  2000/10/18 19:11:27  hristov
78 Division by zero fixed
79
80 Revision 1.23  2000/10/02 21:28:06  fca
81 Removal of useless dependecies via forward declarations
82
83 Revision 1.22  2000/09/12 14:14:55  morsch
84 Call fDecayer->ForceDecay() at the beginning of Generate().
85
86 Revision 1.21  2000/09/08 15:39:01  morsch
87 Handle the case fForceDecay=all during the generation, i.e. select all secondaries.
88
89 Revision 1.20  2000/09/06 14:35:44  morsch
90 Use AliDecayerPythia for particle decays.
91
92 Revision 1.19  2000/07/11 18:24:56  fca
93 Coding convention corrections + few minor bug fixes
94
95 Revision 1.18  2000/06/29 21:08:27  morsch
96 All paramatrisation libraries derive from the pure virtual base class AliGenLib.
97 This allows to pass a pointer to a library directly to AliGenParam and avoids the
98 use of function pointers in Config.C.
99
100 Revision 1.17  2000/06/09 20:33:30  morsch
101 All coding rule violations except RS3 corrected
102
103 Revision 1.16  2000/05/02 07:51:31  morsch
104 - Control precision of pT sampling TF1::SetNpx(..)
105 - Correct initialisation of child-cuts in all constructors.
106 - Most coding rule violations corrected.
107
108 Revision 1.15  2000/04/03 15:42:12  morsch
109 Cuts on primary particles are separated from those on the decay products. Methods
110 SetChildMomentumRange, SetChildPtRange, SetChildPhiRange, SetChildThetaRange added.
111
112 Revision 1.14  1999/11/09 07:38:48  fca
113 Changes for compatibility with version 2.23 of ROOT
114
115 Revision 1.13  1999/11/04 11:30:31  fca
116 Correct the logics for SetForceDecay
117
118 Revision 1.12  1999/11/03 17:43:20  fca
119 New version from G.Martinez & A.Morsch
120
121 Revision 1.11  1999/09/29 09:24:14  fca
122 Introduction of the Copyright and cvs Log
123
124 */
125
126
127
128 // Class to generate particles from using paramtrized pT and y distributions.
129 // Distributions are obtained from pointer to object of type AliGenLib.
130 // (For example AliGenMUONlib)
131 // Decays are performed using Pythia.
132 // andreas.morsch@cern.ch
133
134 #include "AliGenParam.h"
135 #include "AliDecayerPythia.h"
136 #include "AliGenMUONlib.h"
137 #include "AliRun.h"
138 #include <TParticle.h>
139 #include <TParticlePDG.h>
140 #include <TDatabasePDG.h>
141 #include <TLorentzVector.h>
142
143 #include <TF1.h>
144 #include <TCanvas.h>
145 #include <TH1.h>
146
147 ClassImp(AliGenParam)
148
149 //------------------------------------------------------------
150
151   //Begin_Html
152   /*
153     <img src="picts/AliGenParam.gif">
154   */
155   //End_Html
156
157 //____________________________________________________________
158 //____________________________________________________________
159 AliGenParam::AliGenParam()
160 {
161 // Deafault constructor
162     fPtPara = 0;
163     fYPara  = 0;
164     fParam  = 0;
165     fAnalog = kAnalog;
166     SetDeltaPt();
167 //
168 // Set random number generator   
169     sRandom = fRandom;
170     fDecayer = 0;
171
172
173 }
174
175 AliGenParam::AliGenParam(Int_t npart, AliGenLib * Library,  Int_t param, char* tname):AliGenMC(npart)
176 {
177 // Constructor using number of particles parameterisation id and library
178     fName = "Param";
179     fTitle= "Particle Generator using pT and y parameterisation";
180     
181     fPtParaFunc = Library->GetPt(param, tname);
182     fYParaFunc  = Library->GetY (param, tname);
183     fIpParaFunc = Library->GetIp(param, tname);
184     
185     fPtPara = 0;
186     fYPara  = 0;
187     fParam  = param;
188     fAnalog = kAnalog;
189     SetForceDecay();
190     SetDeltaPt(); 
191 //
192 // Set random number generator   
193     sRandom = fRandom;
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     fDecayer = new AliDecayerPythia();
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