717dae0f7d8c02ca97d44e8b1af0739d4b96d7d8
[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 /* $Id$ */
17
18 // Class to generate particles from using paramtrized pT and y distributions.
19 // Distributions are obtained from pointer to object of type AliGenLib.
20 // (For example AliGenMUONlib)
21 // Decays are performed using Pythia.
22 // andreas.morsch@cern.ch
23
24 #include <TCanvas.h>
25 #include <TClonesArray.h>
26 #include <TDatabasePDG.h>
27 #include <TF1.h>
28 #include <TH1F.h>
29 #include <TLorentzVector.h>
30 #include <TMath.h>
31 #include <TParticle.h>
32 #include <TParticlePDG.h>
33 #include <TROOT.h>
34 #include <TVirtualMC.h>
35
36 #include "AliDecayer.h"
37 #include "AliGenMUONlib.h"
38 #include "AliGenParam.h"
39 #include "AliMC.h"
40 #include "AliRun.h"
41 #include "AliGenEventHeader.h"
42
43 ClassImp(AliGenParam)
44
45 //------------------------------------------------------------
46
47 //Begin_Html
48 /*
49     <img src="picts/AliGenParam.gif">
50 */
51 //End_Html
52
53 //____________________________________________________________
54 AliGenParam::AliGenParam()
55 : fPtParaFunc(0),
56         fYParaFunc(0),
57         fIpParaFunc(0),
58   fV2ParaFunc(0),
59         fPtPara(0),
60         fYPara(0),
61   fV2Para(0),
62   fdNdPhi(0),
63         fParam(0),
64         fdNdy0(0.),
65         fYWgt(0.),
66         fPtWgt(0.),
67         fBias(0.),
68         fTrials(0),
69         fDeltaPt(0.01),
70         fSelectAll(kFALSE),
71         fDecayer(0)
72 {
73   // Default constructor
74 }
75 //____________________________________________________________
76 AliGenParam::AliGenParam(Int_t npart, const AliGenLib * Library,  Int_t param, const char* tname)
77     :AliGenMC(npart),
78      fPtParaFunc(Library->GetPt(param, tname)),
79      fYParaFunc (Library->GetY (param, tname)),
80      fIpParaFunc(Library->GetIp(param, tname)),
81    fV2ParaFunc(Library->GetV2(param, tname)),
82      fPtPara(0),
83      fYPara(0),
84    fV2Para(0),
85    fdNdPhi(0),
86      fParam(param),
87      fdNdy0(0.),
88      fYWgt(0.),
89      fPtWgt(0.),
90      fBias(0.),
91      fTrials(0),
92      fDeltaPt(0.01),
93      fSelectAll(kFALSE),
94      fDecayer(0)
95 {
96   // Constructor using number of particles parameterisation id and library
97     fName = "Param";
98     fTitle= "Particle Generator using pT and y parameterisation";
99     fAnalog = kAnalog;
100     SetForceDecay();
101 }
102 //____________________________________________________________
103 AliGenParam::AliGenParam(Int_t npart, Int_t param, const char* tname, const char* name):
104     AliGenMC(npart),
105     fPtParaFunc(0),
106     fYParaFunc (0),
107     fIpParaFunc(0),
108   fV2ParaFunc(0),
109     fPtPara(0),
110     fYPara(0),
111   fV2Para(0),
112   fdNdPhi(0),
113     fParam(param),
114     fdNdy0(0.),
115     fYWgt(0.),
116     fPtWgt(0.),
117     fBias(0.),
118     fTrials(0),
119     fDeltaPt(0.01),
120     fSelectAll(kFALSE),
121     fDecayer(0)
122 {
123   // Constructor using parameterisation id and number of particles
124   //
125     fName = name;
126     fTitle= "Particle Generator using pT and y parameterisation";
127       
128     AliGenLib* pLibrary = new AliGenMUONlib();
129     fPtParaFunc = pLibrary->GetPt(param, tname);
130     fYParaFunc  = pLibrary->GetY (param, tname);
131     fIpParaFunc = pLibrary->GetIp(param, tname);
132   fV2ParaFunc = pLibrary->GetV2(param, tname);
133     
134     fAnalog = kAnalog;
135     fChildSelect.Set(5);
136     for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
137     SetForceDecay();
138     SetCutOnChild();
139     SetChildMomentumRange();
140     SetChildPtRange();
141     SetChildPhiRange();
142     SetChildThetaRange(); 
143 }
144 //____________________________________________________________
145
146 AliGenParam::AliGenParam(Int_t npart, Int_t param,
147                          Double_t (*PtPara) (const Double_t*, const Double_t*),
148                          Double_t (*YPara ) (const Double_t* ,const Double_t*),
149                          Double_t (*V2Para) (const Double_t* ,const Double_t*),
150                          Int_t    (*IpPara) (TRandom *))                 
151     :AliGenMC(npart),
152      
153      fPtParaFunc(PtPara),
154      fYParaFunc(YPara),
155      fIpParaFunc(IpPara),
156    fV2ParaFunc(V2Para),
157      fPtPara(0),
158      fYPara(0),
159    fV2Para(0),
160    fdNdPhi(0),
161      fParam(param),
162      fdNdy0(0.),
163      fYWgt(0.),
164      fPtWgt(0.),
165      fBias(0.),
166      fTrials(0),
167      fDeltaPt(0.01),
168      fSelectAll(kFALSE),
169      fDecayer(0)
170 {
171   // Constructor
172   // Gines Martinez 1/10/99 
173     fName = "Param";
174     fTitle= "Particle Generator using pT and y parameterisation";
175
176     fAnalog = kAnalog;
177     fChildSelect.Set(5);
178     for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
179     SetForceDecay();
180     SetCutOnChild();
181     SetChildMomentumRange();
182     SetChildPtRange();
183     SetChildPhiRange();
184     SetChildThetaRange();  
185 }
186
187 //____________________________________________________________
188 AliGenParam::~AliGenParam()
189 {
190   // Destructor
191     delete  fPtPara;
192     delete  fYPara;
193   delete  fV2Para;
194   delete  fdNdPhi;
195 }
196
197 //____________________________________________________________
198 void AliGenParam::Init()
199 {
200   // Initialisation
201
202     if (gMC) fDecayer = gMC->GetDecayer();
203   //Begin_Html
204   /*
205     <img src="picts/AliGenParam.gif">
206   */
207   //End_Html
208     char name[256];
209     snprintf(name, 256, "pt-parameterisation for %s", GetName());
210     
211     if (fPtPara) fPtPara->Delete();
212     fPtPara = new TF1(name, fPtParaFunc, fPtMin, fPtMax,0);
213     gROOT->GetListOfFunctions()->Remove(fPtPara);
214   //  Set representation precision to 10 MeV
215     Int_t npx= Int_t((fPtMax - fPtMin) / fDeltaPt);
216     
217     fPtPara->SetNpx(npx);
218
219     snprintf(name, 256, "y-parameterisation  for %s", GetName());
220     if (fYPara) fYPara->Delete();
221     fYPara  = new TF1(name, fYParaFunc, fYMin, fYMax, 0);
222     gROOT->GetListOfFunctions()->Remove(fYPara);
223
224   snprintf(name, 256, "v2-parameterisation for %s", GetName());
225   if (fV2Para) fV2Para->Delete();
226   fV2Para  = new TF1(name, fV2ParaFunc, fPtMin, fPtMax, 0);
227   // fV2Para  = new TF1(name, "2*[0]/(1+TMath::Exp([1]*([2]-x)))-[0]", fPtMin, fPtMax);
228   // fV2Para->SetParameter(0, 0.236910);
229   // fV2Para->SetParameter(1, 1.71122);
230   // fV2Para->SetParameter(2, 0.0827617);
231   //gROOT->GetListOfFunctions()->Remove(fV2Para);  //TR: necessary?
232     
233   snprintf(name, 256, "dNdPhi for %s", GetName());
234   if (fdNdPhi) fdNdPhi->Delete();
235   fdNdPhi  = new TF1(name, "1+2*[0]*TMath::Cos(2*(x-[1]))", fPhiMin, fPhiMax);
236   //gROOT->GetListOfFunctions()->Remove(fdNdPhi);  //TR: necessary?
237     
238     snprintf(name, 256, "pt-for-%s", GetName());
239     TF1 ptPara(name ,fPtParaFunc, 0, 15, 0);
240     snprintf(name, 256, "y-for-%s", GetName());
241     TF1 yPara(name, fYParaFunc, -6, 6, 0);
242
243   //
244   // dN/dy| y=0
245     Double_t y1=0;
246     Double_t y2=0;
247     
248     fdNdy0=fYParaFunc(&y1,&y2);
249   //
250   // Integral over generation region
251     Float_t intYS  = yPara.Integral(fYMin, fYMax,(Double_t*) 0x0,1.e-6);
252     Float_t intPt0 = ptPara.Integral(0,15,(Double_t *) 0x0,1.e-6);
253     Float_t intPtS = ptPara.Integral(fPtMin,fPtMax,(Double_t*) 0x0,1.e-6);
254   Float_t phiWgt=(fPhiMax-fPhiMin)/2./TMath::Pi();    //TR: should probably be done differently in case of anisotropic phi...
255     if (fAnalog == kAnalog) {
256         fYWgt  = intYS/fdNdy0;
257         fPtWgt = intPtS/intPt0;
258         fParentWeight = fYWgt*fPtWgt*phiWgt/fNpart;
259     } else {
260         fYWgt = intYS/fdNdy0;
261         fPtWgt = (fPtMax-fPtMin)/intPt0;
262         fParentWeight = fYWgt*fPtWgt*phiWgt/fNpart;
263     }
264   //
265   // particle decay related initialization
266     fDecayer->SetForceDecay(fForceDecay);
267     fDecayer->Init();
268
269   //
270     AliGenMC::Init();
271 }
272
273 //____________________________________________________________
274 void AliGenParam::Generate()
275 {
276   // 
277   // Generate 1 event (see Generate(Int_t ntimes) for details
278   //
279   GenerateN(1);
280 }
281 //____________________________________________________________
282 void AliGenParam::GenerateN(Int_t ntimes)
283 {
284   //
285   // Generate ntimes*'npart' light and heavy mesons (J/Psi, upsilon or phi, Pion,
286   // Kaons, Etas, Omegas) and Baryons (proton, antiprotons, neutrons and 
287   // antineutrons in the the desired theta, phi and momentum windows; 
288   // Gaussian smearing on the vertex is done if selected. 
289   // The decay of heavy mesons is done using lujet, 
290   //    and the childern particle are tracked by GEANT
291   // However, light mesons are directly tracked by GEANT 
292   // setting fForceDecay = nodecay (SetForceDecay(nodecay)) 
293   //
294   //
295   //  Reinitialize decayer
296   fDecayer->SetForceDecay(fForceDecay);
297   fDecayer->Init();
298
299   //
300   Float_t polar[3]= {0,0,0};  // Polarisation of the parent particle (for GEANT tracking)
301   Float_t origin0[3];         // Origin of the generated parent particle (for GEANT tracking)
302   Float_t time0;              // Time0 of the generated parent particle
303   Float_t pt, pl, ptot;       // Transverse, logitudinal and total momenta of the parent particle
304   Float_t phi, theta;         // Phi and theta spherical angles of the parent particle momentum
305   Float_t p[3], pc[3], 
306           och[3];             // Momentum, polarisation and origin of the children particles from lujet
307   Double_t ty, xmt;
308   Int_t nt, i, j;
309   Float_t  wgtp, wgtch;
310   Double_t dummy;
311   static TClonesArray *particles;
312   //
313   if(!particles) particles = new TClonesArray("TParticle",1000);
314   
315   TDatabasePDG *pDataBase = TDatabasePDG::Instance();
316   //
317   Float_t random[6];
318  
319   // Calculating vertex position per event
320   for (j=0;j<3;j++) origin0[j]=fOrigin[j];
321   time0 = fTimeOrigin;
322   if(fVertexSmear==kPerEvent) {
323       Vertex();
324       for (j=0;j<3;j++) origin0[j]=fVertex[j];
325       time0 = fTime;
326   }
327   
328   Int_t ipa=0;
329   
330   // Generating fNpart particles
331   fNprimaries = 0;
332   
333   Int_t nGen = fNpart*ntimes;
334   while (ipa<nGen) {
335       while(1) {
336       //
337       // particle type 
338           Int_t iPart = fIpParaFunc(fRandom);
339           fChildWeight=(fDecayer->GetPartialBranchingRatio(iPart))*fParentWeight;          
340           TParticlePDG *particle = pDataBase->GetParticle(iPart);
341           Float_t am = particle->Mass();
342           
343           Rndm(random,2);
344       //
345       // y
346           ty = TMath::TanH(fYPara->GetRandom());
347       //
348       // pT
349           if (fAnalog == kAnalog) {
350               pt=fPtPara->GetRandom();
351               wgtp=fParentWeight;
352               wgtch=fChildWeight;
353           } else {
354               pt=fPtMin+random[1]*(fPtMax-fPtMin);
355               Double_t ptd=pt;
356               wgtp=fParentWeight*fPtParaFunc(& ptd, &dummy);
357               wgtch=fChildWeight*fPtParaFunc(& ptd, &dummy);
358           }
359           xmt=sqrt(pt*pt+am*am);
360           if (TMath::Abs(ty)==1.) {
361               ty=0.;
362               Fatal("AliGenParam", 
363                     "Division by 0: Please check you rapidity range !");
364           }
365       //
366       // phi
367       // if(!ipa)
368       //        phi=fEvPlane; //align first particle of each event with event plane
369       // else{
370         double v2 = fV2Para->Eval(pt);
371         fdNdPhi->SetParameter(0,v2);
372         fdNdPhi->SetParameter(1,fEvPlane);
373         phi=fdNdPhi->GetRandom();
374       // }
375           
376           pl=xmt*ty/sqrt((1.-ty)*(1.+ty));
377           theta=TMath::ATan2(pt,pl);
378       // Cut on theta
379           if(theta<fThetaMin || theta>fThetaMax) continue;
380           ptot=TMath::Sqrt(pt*pt+pl*pl);
381       // Cut on momentum
382           if(ptot<fPMin || ptot>fPMax) continue;
383       //
384           p[0]=pt*TMath::Cos(phi);
385           p[1]=pt*TMath::Sin(phi);
386           p[2]=pl;
387           if(fVertexSmear==kPerTrack) {
388               Rndm(random,6);
389               for (j=0;j<3;j++) {
390                   origin0[j]=
391                       fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
392                       TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
393               }
394               Rndm(random,2);
395               time0 = fTimeOrigin + fOsigma[2]/TMath::Ccgs()*
396                 TMath::Cos(2*random[0]*TMath::Pi())*
397                 TMath::Sqrt(-2*TMath::Log(random[1]));
398           }
399           
400       // Looking at fForceDecay : 
401       // if fForceDecay != none Primary particle decays using 
402       // AliPythia and children are tracked by GEANT
403       //
404       // if fForceDecay == none Primary particle is tracked by GEANT 
405       // (In the latest, make sure that GEANT actually does all the decays you want)      
406       //
407           Bool_t decayed = kFALSE;
408           
409
410           if (fForceDecay != kNoDecay) {
411         // Using lujet to decay particle
412               Float_t energy=TMath::Sqrt(ptot*ptot+am*am);
413               TLorentzVector pmom(p[0], p[1], p[2], energy);
414               fDecayer->Decay(iPart,&pmom);
415         //
416         // select decay particles
417               Int_t np=fDecayer->ImportParticles(particles);
418
419               //  Selecting  GeometryAcceptance for particles fPdgCodeParticleforAcceptanceCut;
420               if (fGeometryAcceptance) 
421                 if (!CheckAcceptanceGeometry(np,particles)) continue;
422               Int_t ncsel=0;
423               Int_t* pFlag      = new Int_t[np];
424               Int_t* pParent    = new Int_t[np];
425               Int_t* pSelected  = new Int_t[np];
426               Int_t* trackIt    = new Int_t[np];
427
428               for (i=0; i<np; i++) {
429                   pFlag[i]     =  0;
430                   pSelected[i] =  0;
431                   pParent[i]   = -1;
432               }
433               
434               if (np >1) {
435                   decayed = kTRUE;
436                   TParticle* iparticle =  0;
437                   Int_t ipF, ipL;
438                   for (i = 1; i<np ; i++) {
439                       trackIt[i] = 1;
440                       iparticle = (TParticle *) particles->At(i);
441                       Int_t kf = iparticle->GetPdgCode();
442                       Int_t ks = iparticle->GetStatusCode();
443             // flagged particle
444
445                       if (pFlag[i] == 1) {
446                           ipF = iparticle->GetFirstDaughter();
447                           ipL = iparticle->GetLastDaughter();   
448                           if (ipF > 0) for (j=ipF-1; j<ipL; j++) pFlag[j]=1;
449                           continue;
450                       }
451
452             // flag decay products of particles with long life-time (c tau > .3 mum)                  
453                       
454                       if (ks != 1) { 
455               //                          TParticlePDG *particle = pDataBase->GetParticle(kf);
456                           
457                           Double_t lifeTime = fDecayer->GetLifetime(kf);
458               //                          Double_t mass     = particle->Mass();
459               //                          Double_t width    = particle->Width();
460                           if (lifeTime > (Double_t) fMaxLifeTime) {
461                               ipF = iparticle->GetFirstDaughter();
462                               ipL = iparticle->GetLastDaughter();       
463                               if (ipF > 0) for (j=ipF-1; j<ipL; j++) pFlag[j]=1;
464                           } else{
465                               trackIt[i]     = 0;
466                               pSelected[i]   = 1;
467                           }
468                       } // ks==1 ?
469             //
470             // children
471                       
472                       if ((ChildSelected(TMath::Abs(kf)) || fForceDecay == kAll || fSelectAll) && trackIt[i])
473                       {
474                           if (fCutOnChild) {
475                               pc[0]=iparticle->Px();
476                               pc[1]=iparticle->Py();
477                               pc[2]=iparticle->Pz();
478                               Bool_t  childok = KinematicSelection(iparticle, 1);
479                               if(childok) {
480                                   pSelected[i]  = 1;
481                                   ncsel++;
482                               } else {
483                                   ncsel=-1;
484                                   break;
485                               } // child kine cuts
486                           } else {
487                               pSelected[i]  = 1;
488                               ncsel++;
489                           } // if child selection
490                       } // select muon
491                   } // decay particle loop
492               } // if decay products
493               
494               Int_t iparent;
495               if ((fCutOnChild && ncsel >0) || !fCutOnChild){
496                   ipa++;
497           //
498           // Parent
499                   
500                   
501                   PushTrack(0, -1, iPart, p, origin0, polar, time0, kPPrimary, nt, wgtp, ((decayed)? 11 : 1));
502                   pParent[0] = nt;
503                   KeepTrack(nt); 
504                   fNprimaries++;
505                   
506           //
507           // Decay Products
508           //              
509                   for (i = 1; i < np; i++) {
510                       if (pSelected[i]) {
511                           TParticle* iparticle = (TParticle *) particles->At(i);
512                           Int_t kf   = iparticle->GetPdgCode();
513                           Int_t ksc  = iparticle->GetStatusCode();
514                           Int_t jpa  = iparticle->GetFirstMother()-1;
515                           
516                           och[0] = origin0[0]+iparticle->Vx();
517                           och[1] = origin0[1]+iparticle->Vy();
518                           och[2] = origin0[2]+iparticle->Vz();
519                           pc[0]  = iparticle->Px();
520                           pc[1]  = iparticle->Py();
521                           pc[2]  = iparticle->Pz();
522                           
523                           if (jpa > -1) {
524                               iparent = pParent[jpa];
525                           } else {
526                               iparent = -1;
527                           }
528                          
529                           PushTrack(fTrackIt * trackIt[i], iparent, kf,
530                                            pc, och, polar,
531                                            time0 + iparticle->T(), kPDecay, nt, wgtch, ksc);
532                           pParent[i] = nt;
533                           KeepTrack(nt); 
534                           fNprimaries++;
535                       } // Selected
536                   } // Particle loop 
537               }  // Decays by Lujet
538               particles->Clear();
539               if (pFlag)      delete[] pFlag;
540               if (pParent)    delete[] pParent;
541               if (pSelected)  delete[] pSelected;          
542               if (trackIt)    delete[] trackIt;
543           } // kinematic selection
544           else  // nodecay option, so parent will be tracked by GEANT (pions, kaons, eta, omegas, baryons)
545           {
546             gAlice->GetMCApp()->
547                 PushTrack(fTrackIt,-1,iPart,p,origin0,polar,time0,kPPrimary,nt,wgtp, 1);
548             ipa++; 
549             fNprimaries++;
550           }
551           break;
552     } // while
553   } // event loop
554   
555   SetHighWaterMark(nt);
556
557   AliGenEventHeader* header = new AliGenEventHeader("PARAM");
558   header->SetPrimaryVertex(fVertex);
559   header->SetInteractionTime(fTime);
560   header->SetNProduced(fNprimaries);
561   AddHeader(header);
562 }
563 //____________________________________________________________________________________
564 Float_t AliGenParam::GetRelativeArea(Float_t ptMin, Float_t ptMax, Float_t yMin, Float_t yMax, Float_t phiMin, Float_t phiMax)
565 {
566   //
567   // Normalisation for selected kinematic region
568   //
569   Float_t ratio =  
570     fPtPara->Integral(ptMin,ptMax,(Double_t *)0,1.e-6) / fPtPara->Integral( fPtPara->GetXmin(), fPtPara->GetXmax(),(Double_t *)0,1.e-6) *
571     fYPara->Integral(yMin,yMax,(Double_t *)0,1.e-6)/fYPara->Integral(fYPara->GetXmin(),fYPara->GetXmax(),(Double_t *)0,1.e-6)   *
572     (phiMax-phiMin)/360.;
573   return TMath::Abs(ratio);
574 }
575
576 //____________________________________________________________________________________
577
578 void AliGenParam::Draw( const char * /*opt*/)
579 {
580     //
581     // Draw the pT and y Distributions
582     //
583      TCanvas *c0 = new TCanvas("c0","Canvas 0",400,10,600,700);
584      c0->Divide(2,1);
585      c0->cd(1);
586      fPtPara->Draw();
587      fPtPara->GetHistogram()->SetXTitle("p_{T} (GeV)");     
588      c0->cd(2);
589      fYPara->Draw();
590      fYPara->GetHistogram()->SetXTitle("y");     
591 }
592
593
594
595