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