]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - EVGEN/AliGenParam.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
update from Asako's tasks
[u/mrichter/AliRoot.git] / EVGEN / AliGenParam.cxx
... / ...
CommitLineData
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
43ClassImp(AliGenParam)
44
45//------------------------------------------------------------
46
47//Begin_Html
48/*
49 <img src="picts/AliGenParam.gif">
50*/
51//End_Html
52
53//____________________________________________________________
54AliGenParam::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//____________________________________________________________
76AliGenParam::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//____________________________________________________________
103AliGenParam::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
146AliGenParam::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//____________________________________________________________
188AliGenParam::~AliGenParam()
189{
190 // Destructor
191 delete fPtPara;
192 delete fYPara;
193 delete fV2Para;
194 delete fdNdPhi;
195}
196
197//____________________________________________________________
198void 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//____________________________________________________________
280void AliGenParam::Generate()
281{
282 //
283 // Generate 1 event (see Generate(Int_t ntimes) for details
284 //
285 GenerateN(1);
286}
287//____________________________________________________________
288void 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//____________________________________________________________________________________
570Float_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
591void 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