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