Changes for
[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 #if ROOT_VERSION_CODE < ROOT_VERSION(5,99,0)
252     Float_t intYS  = yPara.Integral(fYMin, fYMax,(Double_t*) 0x0,1.e-6);
253     Float_t intPt0 = ptPara.Integral(0,15,(Double_t *) 0x0,1.e-6);
254     Float_t intPtS = ptPara.Integral(fPtMin,fPtMax,(Double_t*) 0x0,1.e-6);
255 #else
256     Float_t intYS  = yPara.Integral(fYMin, fYMax,1.e-6);
257     Float_t intPt0 = ptPara.Integral(0,15,1.e-6);
258     Float_t intPtS = ptPara.Integral(fPtMin,fPtMax,1.e-6);
259 #endif
260   Float_t phiWgt=(fPhiMax-fPhiMin)/2./TMath::Pi();    //TR: should probably be done differently in case of anisotropic phi...
261     if (fAnalog == kAnalog) {
262         fYWgt  = intYS/fdNdy0;
263         fPtWgt = intPtS/intPt0;
264         fParentWeight = fYWgt*fPtWgt*phiWgt/fNpart;
265     } else {
266         fYWgt = intYS/fdNdy0;
267         fPtWgt = (fPtMax-fPtMin)/intPt0;
268         fParentWeight = fYWgt*fPtWgt*phiWgt/fNpart;
269     }
270   //
271   // particle decay related initialization
272     fDecayer->SetForceDecay(fForceDecay);
273     fDecayer->Init();
274
275   //
276     AliGenMC::Init();
277 }
278
279 //____________________________________________________________
280 void AliGenParam::Generate()
281 {
282   // 
283   // Generate 1 event (see Generate(Int_t ntimes) for details
284   //
285   GenerateN(1);
286 }
287 //____________________________________________________________
288 void AliGenParam::GenerateN(Int_t ntimes)
289 {
290   //
291   // Generate ntimes*'npart' light and heavy mesons (J/Psi, upsilon or phi, Pion,
292   // Kaons, Etas, Omegas) and Baryons (proton, antiprotons, neutrons and 
293   // antineutrons in the the desired theta, phi and momentum windows; 
294   // Gaussian smearing on the vertex is done if selected. 
295   // The decay of heavy mesons is done using lujet, 
296   //    and the childern particle are tracked by GEANT
297   // However, light mesons are directly tracked by GEANT 
298   // setting fForceDecay = nodecay (SetForceDecay(nodecay)) 
299   //
300   //
301   //  Reinitialize decayer
302   fDecayer->SetForceDecay(fForceDecay);
303   fDecayer->Init();
304
305   //
306   Float_t polar[3]= {0,0,0};  // Polarisation of the parent particle (for GEANT tracking)
307   Float_t origin0[3];         // Origin of the generated parent particle (for GEANT tracking)
308   Float_t time0;              // Time0 of the generated parent particle
309   Float_t pt, pl, ptot;       // Transverse, logitudinal and total momenta of the parent particle
310   Float_t phi, theta;         // Phi and theta spherical angles of the parent particle momentum
311   Float_t p[3], pc[3], 
312           och[3];             // Momentum, polarisation and origin of the children particles from lujet
313   Double_t ty, xmt;
314   Int_t nt, i, j;
315   Float_t  wgtp, wgtch;
316   Double_t dummy;
317   static TClonesArray *particles;
318   //
319   if(!particles) particles = new TClonesArray("TParticle",1000);
320   
321   TDatabasePDG *pDataBase = TDatabasePDG::Instance();
322   //
323   Float_t random[6];
324  
325   // Calculating vertex position per event
326   for (j=0;j<3;j++) origin0[j]=fOrigin[j];
327   time0 = fTimeOrigin;
328   if(fVertexSmear==kPerEvent) {
329       Vertex();
330       for (j=0;j<3;j++) origin0[j]=fVertex[j];
331       time0 = fTime;
332   }
333   
334   Int_t ipa=0;
335   
336   // Generating fNpart particles
337   fNprimaries = 0;
338   
339   Int_t nGen = fNpart*ntimes;
340   while (ipa<nGen) {
341       while(1) {
342       //
343       // particle type 
344           Int_t iPart = fIpParaFunc(fRandom);
345           fChildWeight=(fDecayer->GetPartialBranchingRatio(iPart))*fParentWeight;          
346           TParticlePDG *particle = pDataBase->GetParticle(iPart);
347           Float_t am = particle->Mass();
348           
349           Rndm(random,2);
350       //
351       // y
352           ty = TMath::TanH(fYPara->GetRandom());
353       //
354       // pT
355           if (fAnalog == kAnalog) {
356               pt=fPtPara->GetRandom();
357               wgtp=fParentWeight;
358               wgtch=fChildWeight;
359           } else {
360               pt=fPtMin+random[1]*(fPtMax-fPtMin);
361               Double_t ptd=pt;
362               wgtp=fParentWeight*fPtParaFunc(& ptd, &dummy);
363               wgtch=fChildWeight*fPtParaFunc(& ptd, &dummy);
364           }
365           xmt=sqrt(pt*pt+am*am);
366           if (TMath::Abs(ty)==1.) {
367               ty=0.;
368               Fatal("AliGenParam", 
369                     "Division by 0: Please check you rapidity range !");
370           }
371       //
372       // phi
373       // if(!ipa)
374       //        phi=fEvPlane; //align first particle of each event with event plane
375       // else{
376         double v2 = fV2Para->Eval(pt);
377         fdNdPhi->SetParameter(0,v2);
378         fdNdPhi->SetParameter(1,fEvPlane);
379         phi=fdNdPhi->GetRandom();
380       // }
381           
382           pl=xmt*ty/sqrt((1.-ty)*(1.+ty));
383           theta=TMath::ATan2(pt,pl);
384       // Cut on theta
385           if(theta<fThetaMin || theta>fThetaMax) continue;
386           ptot=TMath::Sqrt(pt*pt+pl*pl);
387       // Cut on momentum
388           if(ptot<fPMin || ptot>fPMax) continue;
389       //
390           p[0]=pt*TMath::Cos(phi);
391           p[1]=pt*TMath::Sin(phi);
392           p[2]=pl;
393           if(fVertexSmear==kPerTrack) {
394               Rndm(random,6);
395               for (j=0;j<3;j++) {
396                   origin0[j]=
397                       fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
398                       TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
399               }
400               Rndm(random,2);
401               time0 = fTimeOrigin + fOsigma[2]/TMath::Ccgs()*
402                 TMath::Cos(2*random[0]*TMath::Pi())*
403                 TMath::Sqrt(-2*TMath::Log(random[1]));
404           }
405           
406       // Looking at fForceDecay : 
407       // if fForceDecay != none Primary particle decays using 
408       // AliPythia and children are tracked by GEANT
409       //
410       // if fForceDecay == none Primary particle is tracked by GEANT 
411       // (In the latest, make sure that GEANT actually does all the decays you want)      
412       //
413           Bool_t decayed = kFALSE;
414           
415
416           if (fForceDecay != kNoDecay) {
417         // Using lujet to decay particle
418               Float_t energy=TMath::Sqrt(ptot*ptot+am*am);
419               TLorentzVector pmom(p[0], p[1], p[2], energy);
420               fDecayer->Decay(iPart,&pmom);
421         //
422         // select decay particles
423               Int_t np=fDecayer->ImportParticles(particles);
424
425               //  Selecting  GeometryAcceptance for particles fPdgCodeParticleforAcceptanceCut;
426               if (fGeometryAcceptance) 
427                 if (!CheckAcceptanceGeometry(np,particles)) continue;
428               Int_t ncsel=0;
429               Int_t* pFlag      = new Int_t[np];
430               Int_t* pParent    = new Int_t[np];
431               Int_t* pSelected  = new Int_t[np];
432               Int_t* trackIt    = new Int_t[np];
433
434               for (i=0; i<np; i++) {
435                   pFlag[i]     =  0;
436                   pSelected[i] =  0;
437                   pParent[i]   = -1;
438               }
439               
440               if (np >1) {
441                   decayed = kTRUE;
442                   TParticle* iparticle =  0;
443                   Int_t ipF, ipL;
444                   for (i = 1; i<np ; i++) {
445                       trackIt[i] = 1;
446                       iparticle = (TParticle *) particles->At(i);
447                       Int_t kf = iparticle->GetPdgCode();
448                       Int_t ks = iparticle->GetStatusCode();
449             // flagged particle
450
451                       if (pFlag[i] == 1) {
452                           ipF = iparticle->GetFirstDaughter();
453                           ipL = iparticle->GetLastDaughter();   
454                           if (ipF > 0) for (j=ipF-1; j<ipL; j++) pFlag[j]=1;
455                           continue;
456                       }
457
458             // flag decay products of particles with long life-time (c tau > .3 mum)                  
459                       
460                       if (ks != 1) { 
461               //                          TParticlePDG *particle = pDataBase->GetParticle(kf);
462                           
463                           Double_t lifeTime = fDecayer->GetLifetime(kf);
464               //                          Double_t mass     = particle->Mass();
465               //                          Double_t width    = particle->Width();
466                           if (lifeTime > (Double_t) fMaxLifeTime) {
467                               ipF = iparticle->GetFirstDaughter();
468                               ipL = iparticle->GetLastDaughter();       
469                               if (ipF > 0) for (j=ipF-1; j<ipL; j++) pFlag[j]=1;
470                           } else{
471                               trackIt[i]     = 0;
472                               pSelected[i]   = 1;
473                           }
474                       } // ks==1 ?
475             //
476             // children
477                       
478                       if ((ChildSelected(TMath::Abs(kf)) || fForceDecay == kAll || fSelectAll) && trackIt[i])
479                       {
480                           if (fCutOnChild) {
481                               pc[0]=iparticle->Px();
482                               pc[1]=iparticle->Py();
483                               pc[2]=iparticle->Pz();
484                               Bool_t  childok = KinematicSelection(iparticle, 1);
485                               if(childok) {
486                                   pSelected[i]  = 1;
487                                   ncsel++;
488                               } else {
489                                   ncsel=-1;
490                                   break;
491                               } // child kine cuts
492                           } else {
493                               pSelected[i]  = 1;
494                               ncsel++;
495                           } // if child selection
496                       } // select muon
497                   } // decay particle loop
498               } // if decay products
499               
500               Int_t iparent;
501               if ((fCutOnChild && ncsel >0) || !fCutOnChild){
502                   ipa++;
503           //
504           // Parent
505                   
506                   
507                   PushTrack(0, -1, iPart, p, origin0, polar, time0, kPPrimary, nt, wgtp, ((decayed)? 11 : 1));
508                   pParent[0] = nt;
509                   KeepTrack(nt); 
510                   fNprimaries++;
511                   
512           //
513           // Decay Products
514           //              
515                   for (i = 1; i < np; i++) {
516                       if (pSelected[i]) {
517                           TParticle* iparticle = (TParticle *) particles->At(i);
518                           Int_t kf   = iparticle->GetPdgCode();
519                           Int_t ksc  = iparticle->GetStatusCode();
520                           Int_t jpa  = iparticle->GetFirstMother()-1;
521                           
522                           och[0] = origin0[0]+iparticle->Vx();
523                           och[1] = origin0[1]+iparticle->Vy();
524                           och[2] = origin0[2]+iparticle->Vz();
525                           pc[0]  = iparticle->Px();
526                           pc[1]  = iparticle->Py();
527                           pc[2]  = iparticle->Pz();
528                           
529                           if (jpa > -1) {
530                               iparent = pParent[jpa];
531                           } else {
532                               iparent = -1;
533                           }
534                          
535                           PushTrack(fTrackIt * trackIt[i], iparent, kf,
536                                            pc, och, polar,
537                                            time0 + iparticle->T(), kPDecay, nt, wgtch, ksc);
538                           pParent[i] = nt;
539                           KeepTrack(nt); 
540                           fNprimaries++;
541                       } // Selected
542                   } // Particle loop 
543               }  // Decays by Lujet
544               particles->Clear();
545               if (pFlag)      delete[] pFlag;
546               if (pParent)    delete[] pParent;
547               if (pSelected)  delete[] pSelected;          
548               if (trackIt)    delete[] trackIt;
549           } // kinematic selection
550           else  // nodecay option, so parent will be tracked by GEANT (pions, kaons, eta, omegas, baryons)
551           {
552             gAlice->GetMCApp()->
553                 PushTrack(fTrackIt,-1,iPart,p,origin0,polar,time0,kPPrimary,nt,wgtp, 1);
554             ipa++; 
555             fNprimaries++;
556           }
557           break;
558     } // while
559   } // event loop
560   
561   SetHighWaterMark(nt);
562
563   AliGenEventHeader* header = new AliGenEventHeader("PARAM");
564   header->SetPrimaryVertex(fVertex);
565   header->SetInteractionTime(fTime);
566   header->SetNProduced(fNprimaries);
567   AddHeader(header);
568 }
569 //____________________________________________________________________________________
570 Float_t AliGenParam::GetRelativeArea(Float_t ptMin, Float_t ptMax, Float_t yMin, Float_t yMax, Float_t phiMin, Float_t phiMax)
571 {
572   //
573   // Normalisation for selected kinematic region
574   //
575 #if ROOT_VERSION_CODE < ROOT_VERSION(5,99,0)
576   Float_t ratio =  
577     fPtPara->Integral(ptMin,ptMax,(Double_t *)0,1.e-6) / fPtPara->Integral( fPtPara->GetXmin(), fPtPara->GetXmax(),(Double_t *)0,1.e-6) *
578     fYPara->Integral(yMin,yMax,(Double_t *)0,1.e-6)/fYPara->Integral(fYPara->GetXmin(),fYPara->GetXmax(),(Double_t *)0,1.e-6)   *
579     (phiMax-phiMin)/360.;
580 #else
581   Float_t ratio =  
582     fPtPara->Integral(ptMin,ptMax,1.e-6) / fPtPara->Integral( fPtPara->GetXmin(), fPtPara->GetXmax(),1.e-6) *
583     fYPara->Integral(yMin,yMax,1.e-6)/fYPara->Integral(fYPara->GetXmin(),fYPara->GetXmax(),1.e-6)   *
584     (phiMax-phiMin)/360.;
585 #endif
586   return TMath::Abs(ratio);
587 }
588
589 //____________________________________________________________________________________
590
591 void AliGenParam::Draw( const char * /*opt*/)
592 {
593     //
594     // Draw the pT and y Distributions
595     //
596      TCanvas *c0 = new TCanvas("c0","Canvas 0",400,10,600,700);
597      c0->Divide(2,1);
598      c0->cd(1);
599      fPtPara->Draw();
600      fPtPara->GetHistogram()->SetXTitle("p_{T} (GeV)");     
601      c0->cd(2);
602      fYPara->Draw();
603      fYPara->GetHistogram()->SetXTitle("y");     
604 }
605
606
607
608