Results saved in output root files
[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)
265 Float_t pt, pl, ptot; // Transverse, logitudinal and total momenta of the parent particle
266 Float_t phi, theta; // Phi and theta spherical angles of the parent particle momentum
267 Float_t p[3], pc[3],
fff02fee 268 och[3]; // Momentum, polarisation and origin of the children particles from lujet
2d7a47be 269 Double_t ty, xmt;
fff02fee 270 Int_t nt, i, j;
fe4da5cc 271 Float_t wgtp, wgtch;
272 Double_t dummy;
09fd3ea2 273 static TClonesArray *particles;
fe4da5cc 274 //
fff02fee 275 if(!particles) particles = new TClonesArray("TParticle",1000);
8b31bfac 276
349be858 277 TDatabasePDG *pDataBase = TDatabasePDG::Instance();
fe4da5cc 278 //
279 Float_t random[6];
28337bc1 280
281// Calculating vertex position per event
fe4da5cc 282 for (j=0;j<3;j++) origin0[j]=fOrigin[j];
aee8290b 283 if(fVertexSmear==kPerEvent) {
d9ea0e3b 284 Vertex();
285 for (j=0;j<3;j++) origin0[j]=fVertex[j];
fe4da5cc 286 }
d9ea0e3b 287
21aaa175 288 Int_t ipa=0;
fff02fee 289
28337bc1 290// Generating fNpart particles
7cf36cae 291 fNprimaries = 0;
292
21aaa175 293 while (ipa<fNpart) {
cc5d764c 294 while(1) {
fe4da5cc 295//
2ad38d56 296// particle type
65fb704d 297 Int_t iPart = fIpParaFunc(fRandom);
fa480368 298 fChildWeight=(fDecayer->GetPartialBranchingRatio(iPart))*fParentWeight;
675e9664 299 TParticlePDG *particle = pDataBase->GetParticle(iPart);
fa480368 300 Float_t am = particle->Mass();
8b31bfac 301
65fb704d 302 Rndm(random,2);
fe4da5cc 303//
304// phi
305 phi=fPhiMin+random[0]*(fPhiMax-fPhiMin);
306//
307// y
2d7a47be 308 ty = TMath::TanH(fYPara->GetRandom());
fe4da5cc 309//
310// pT
34f60c01 311 if (fAnalog == kAnalog) {
fe4da5cc 312 pt=fPtPara->GetRandom();
313 wgtp=fParentWeight;
314 wgtch=fChildWeight;
315 } else {
316 pt=fPtMin+random[1]*(fPtMax-fPtMin);
317 Double_t ptd=pt;
318 wgtp=fParentWeight*fPtParaFunc(& ptd, &dummy);
319 wgtch=fChildWeight*fPtParaFunc(& ptd, &dummy);
320 }
321 xmt=sqrt(pt*pt+am*am);
2d7a47be 322 if (TMath::Abs(ty)==1.) {
323 ty=0.;
324 Fatal("AliGenParam",
325 "Division by 0: Please check you rapidity range !");
326 }
327
60e55aee 328 pl=xmt*ty/sqrt((1.-ty)*(1.+ty));
fe4da5cc 329 theta=TMath::ATan2(pt,pl);
cc5d764c 330// Cut on theta
fe4da5cc 331 if(theta<fThetaMin || theta>fThetaMax) continue;
332 ptot=TMath::Sqrt(pt*pt+pl*pl);
cc5d764c 333// Cut on momentum
fe4da5cc 334 if(ptot<fPMin || ptot>fPMax) continue;
fff02fee 335//
fe4da5cc 336 p[0]=pt*TMath::Cos(phi);
337 p[1]=pt*TMath::Sin(phi);
338 p[2]=pl;
aee8290b 339 if(fVertexSmear==kPerTrack) {
65fb704d 340 Rndm(random,6);
fe4da5cc 341 for (j=0;j<3;j++) {
342 origin0[j]=
343 fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
344 TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
345 }
346 }
28337bc1 347
348// Looking at fForceDecay :
cc5d764c 349// if fForceDecay != none Primary particle decays using
350// AliPythia and children are tracked by GEANT
28337bc1 351//
cc5d764c 352// if fForceDecay == none Primary particle is tracked by GEANT
353// (In the latest, make sure that GEANT actually does all the decays you want)
fe4da5cc 354//
cd716030 355 Bool_t decayed = kFALSE;
356
fff02fee 357
34f60c01 358 if (fForceDecay != kNoDecay) {
28337bc1 359// Using lujet to decay particle
886b6f73 360 Float_t energy=TMath::Sqrt(ptot*ptot+am*am);
fa480368 361 TLorentzVector pmom(p[0], p[1], p[2], energy);
362 fDecayer->Decay(iPart,&pmom);
fe4da5cc 363//
cc5d764c 364// select decay particles
fa480368 365 Int_t np=fDecayer->ImportParticles(particles);
1242532d 366
367 // Selecting GeometryAcceptance for particles fPdgCodeParticleforAcceptanceCut;
368 if (fGeometryAcceptance)
369 if (!CheckAcceptanceGeometry(np,particles)) continue;
886b6f73 370 Int_t ncsel=0;
fff02fee 371 Int_t* pFlag = new Int_t[np];
372 Int_t* pParent = new Int_t[np];
373 Int_t* pSelected = new Int_t[np];
374 Int_t* trackIt = new Int_t[np];
375
376 for (i=0; i<np; i++) {
377 pFlag[i] = 0;
378 pSelected[i] = 0;
379 pParent[i] = -1;
380 }
381
6ba00c52 382 if (np >1) {
cd716030 383 decayed = kTRUE;
08bffa4d 384 TParticle* iparticle = 0;
fff02fee 385 Int_t ipF, ipL;
386 for (i = 1; i<np ; i++) {
387 trackIt[i] = 1;
6ba00c52 388 iparticle = (TParticle *) particles->At(i);
389 Int_t kf = iparticle->GetPdgCode();
fff02fee 390 Int_t ks = iparticle->GetStatusCode();
391// flagged particle
392
393 if (pFlag[i] == 1) {
fff02fee 394 ipF = iparticle->GetFirstDaughter();
395 ipL = iparticle->GetLastDaughter();
396 if (ipF > 0) for (j=ipF-1; j<ipL; j++) pFlag[j]=1;
397 continue;
398 }
399
400// flag decay products of particles with long life-time (c tau > .3 mum)
401
402 if (ks != 1) {
2d7a47be 403// TParticlePDG *particle = pDataBase->GetParticle(kf);
fff02fee 404
405 Double_t lifeTime = fDecayer->GetLifetime(kf);
2d7a47be 406// Double_t mass = particle->Mass();
407// Double_t width = particle->Width();
47fc6bd5 408 if (lifeTime > (Double_t) fMaxLifeTime) {
fff02fee 409 ipF = iparticle->GetFirstDaughter();
410 ipL = iparticle->GetLastDaughter();
411 if (ipF > 0) for (j=ipF-1; j<ipL; j++) pFlag[j]=1;
412 } else{
413 trackIt[i] = 0;
414 pSelected[i] = 1;
415 }
416 } // ks==1 ?
fe4da5cc 417//
418// children
47fc6bd5 419
18e09c20 420 if ((ChildSelected(TMath::Abs(kf)) || fForceDecay == kAll || fSelectAll) && trackIt[i])
6ba00c52 421 {
6ba00c52 422 if (fCutOnChild) {
fff02fee 423 pc[0]=iparticle->Px();
424 pc[1]=iparticle->Py();
425 pc[2]=iparticle->Pz();
426 Bool_t childok = KinematicSelection(iparticle, 1);
427 if(childok) {
428 pSelected[i] = 1;
6ba00c52 429 ncsel++;
430 } else {
431 ncsel=-1;
432 break;
433 } // child kine cuts
434 } else {
fff02fee 435 pSelected[i] = 1;
886b6f73 436 ncsel++;
6ba00c52 437 } // if child selection
438 } // select muon
439 } // decay particle loop
440 } // if decay products
441
886b6f73 442 Int_t iparent;
443 if ((fCutOnChild && ncsel >0) || !fCutOnChild){
444 ipa++;
21aaa175 445//
fff02fee 446// Parent
cd716030 447
448
449 PushTrack(0, -1, iPart, p, origin0, polar, 0, kPPrimary, nt, wgtp, ((decayed)? 11 : 1));
fff02fee 450 pParent[0] = nt;
a99cf51f 451 KeepTrack(nt);
7cf36cae 452 fNprimaries++;
453
fff02fee 454//
455// Decay Products
456//
457 for (i = 1; i < np; i++) {
458 if (pSelected[i]) {
459 TParticle* iparticle = (TParticle *) particles->At(i);
ea79897e 460 Int_t kf = iparticle->GetPdgCode();
cd716030 461 Int_t ksc = iparticle->GetStatusCode();
ea79897e 462 Int_t jpa = iparticle->GetFirstMother()-1;
fff02fee 463
464 och[0] = origin0[0]+iparticle->Vx()/10;
465 och[1] = origin0[1]+iparticle->Vy()/10;
466 och[2] = origin0[2]+iparticle->Vz()/10;
467 pc[0] = iparticle->Px();
468 pc[1] = iparticle->Py();
469 pc[2] = iparticle->Pz();
470
ea79897e 471 if (jpa > -1) {
472 iparent = pParent[jpa];
fff02fee 473 } else {
474 iparent = -1;
475 }
1242532d 476
ea79897e 477 PushTrack(fTrackIt * trackIt[i], iparent, kf,
fff02fee 478 pc, och, polar,
cd716030 479 0, kPDecay, nt, wgtch, ksc);
fff02fee 480 pParent[i] = nt;
a99cf51f 481 KeepTrack(nt);
7cf36cae 482 fNprimaries++;
8b31bfac 483 } // Selected
484 } // Particle loop
28337bc1 485 } // Decays by Lujet
8b31bfac 486 particles->Clear();
fff02fee 487 if (pFlag) delete[] pFlag;
488 if (pParent) delete[] pParent;
489 if (pSelected) delete[] pSelected;
8b31bfac 490 if (trackIt) delete[] trackIt;
21aaa175 491 } // kinematic selection
28337bc1 492 else // nodecay option, so parent will be tracked by GEANT (pions, kaons, eta, omegas, baryons)
493 {
5d12ce38 494 gAlice->GetMCApp()->
cd716030 495 PushTrack(fTrackIt,-1,iPart,p,origin0,polar,0,kPPrimary,nt,wgtp, 1);
28337bc1 496 ipa++;
7cf36cae 497 fNprimaries++;
28337bc1 498 }
fe4da5cc 499 break;
21aaa175 500 } // while
fe4da5cc 501 } // event loop
7cf36cae 502
a99cf51f 503 SetHighWaterMark(nt);
7cf36cae 504
505 AliGenEventHeader* header = new AliGenEventHeader("PARAM");
506 header->SetPrimaryVertex(fVertex);
507 header->SetNProduced(fNprimaries);
508 AddHeader(header);
fe4da5cc 509}
2ad38d56 510//____________________________________________________________________________________
511Float_t AliGenParam::GetRelativeArea(Float_t ptMin, Float_t ptMax, Float_t yMin, Float_t yMax, Float_t phiMin, Float_t phiMax)
512{
513//
514// Normalisation for selected kinematic region
515//
516 Float_t ratio =
389a9124 517 fPtPara->Integral(ptMin,ptMax,(Double_t *)0,1.e-6) / fPtPara->Integral( fPtPara->GetXmin(), fPtPara->GetXmax(),(Double_t *)0,1.e-6) *
518 fYPara->Integral(yMin,yMax,(Double_t *)0,1.e-6)/fYPara->Integral(fYPara->GetXmin(),fYPara->GetXmax(),(Double_t *)0,1.e-6) *
2ad38d56 519 (phiMax-phiMin)/360.;
520 return TMath::Abs(ratio);
521}
522
523//____________________________________________________________________________________
fe4da5cc 524
dc1d768c 525void AliGenParam::Draw( const char * /*opt*/)
5bd39445 526{
527 //
528 // Draw the pT and y Distributions
529 //
530 TCanvas *c0 = new TCanvas("c0","Canvas 0",400,10,600,700);
531 c0->Divide(2,1);
532 c0->cd(1);
533 fPtPara->Draw();
534 fPtPara->GetHistogram()->SetXTitle("p_{T} (GeV)");
535 c0->cd(2);
536 fYPara->Draw();
537 fYPara->GetHistogram()->SetXTitle("y");
538}
539
f87cfe57 540
fe4da5cc 541
542