]> 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
Removing some meaningeless const (coverity)
[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 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//____________________________________________________________
274void 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();
508 och[1] = origin0[1]+iparticle->Vy();
509 och[2] = origin0[2]+iparticle->Vz();
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//____________________________________________________________________________________
555Float_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
569void 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