]> 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
add dphi histograms for deta > 0.8
[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//____________________________________________________________
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),
67 fSelectAll(kFALSE),
68 fDecayer(0)
69{
70// Default constructor
71}
72//____________________________________________________________
73AliGenParam::AliGenParam(Int_t npart, const AliGenLib * Library, Int_t param, const char* tname)
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),
87 fSelectAll(kFALSE),
88 fDecayer(0)
89{
90// Constructor using number of particles parameterisation id and library
91 fName = "Param";
92 fTitle= "Particle Generator using pT and y parameterisation";
93 fAnalog = kAnalog;
94 SetForceDecay();
95}
96//____________________________________________________________
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),
111 fSelectAll(kFALSE),
112 fDecayer(0)
113{
114// Constructor using parameterisation id and number of particles
115//
116 fName = name;
117 fTitle= "Particle Generator using pT and y parameterisation";
118
119 AliGenLib* pLibrary = new AliGenMUONlib();
120 fPtParaFunc = pLibrary->GetPt(param, tname);
121 fYParaFunc = pLibrary->GetY (param, tname);
122 fIpParaFunc = pLibrary->GetIp(param, tname);
123
124 fAnalog = kAnalog;
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();
133}
134//____________________________________________________________
135
136AliGenParam::AliGenParam(Int_t npart, Int_t param,
137 Double_t (*PtPara) (const Double_t*, const Double_t*),
138 Double_t (*YPara ) (const Double_t* ,const Double_t*),
139 Int_t (*IpPara) (TRandom *))
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),
154 fSelectAll(kFALSE),
155 fDecayer(0)
156{
157// Constructor
158// Gines Martinez 1/10/99
159 fName = "Param";
160 fTitle= "Particle Generator using pT and y parameterisation";
161
162 fAnalog = kAnalog;
163 fChildSelect.Set(5);
164 for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
165 SetForceDecay();
166 SetCutOnChild();
167 SetChildMomentumRange();
168 SetChildPtRange();
169 SetChildPhiRange();
170 SetChildThetaRange();
171}
172
173//____________________________________________________________
174AliGenParam::~AliGenParam()
175{
176// Destructor
177 delete fPtPara;
178 delete fYPara;
179}
180
181//____________________________________________________________
182void AliGenParam::Init()
183{
184// Initialisation
185
186 if (gMC) fDecayer = gMC->GetDecayer();
187 //Begin_Html
188 /*
189 <img src="picts/AliGenParam.gif">
190 */
191 //End_Html
192 char name[256];
193 snprintf(name, 256, "pt-parameterisation for %s", GetName());
194
195 if (fPtPara) fPtPara->Delete();
196 fPtPara = new TF1(name, fPtParaFunc, fPtMin, fPtMax,0);
197 gROOT->GetListOfFunctions()->Remove(fPtPara);
198// Set representation precision to 10 MeV
199 Int_t npx= Int_t((fPtMax - fPtMin) / fDeltaPt);
200
201 fPtPara->SetNpx(npx);
202
203 snprintf(name, 256, "y-parameterisation for %s", GetName());
204 if (fYPara) fYPara->Delete();
205 fYPara = new TF1(name, fYParaFunc, fYMin, fYMax, 0);
206 gROOT->GetListOfFunctions()->Remove(fYPara);
207
208
209 snprintf(name, 256, "pt-for-%s", GetName());
210 TF1 ptPara(name ,fPtParaFunc, 0, 15, 0);
211 snprintf(name, 256, "y-for-%s", GetName());
212 TF1 yPara(name, fYParaFunc, -6, 6, 0);
213
214//
215// dN/dy| y=0
216 Double_t y1=0;
217 Double_t y2=0;
218
219 fdNdy0=fYParaFunc(&y1,&y2);
220//
221// Integral over generation region
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);
225 Float_t phiWgt=(fPhiMax-fPhiMin)/2./TMath::Pi();
226 if (fAnalog == kAnalog) {
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 }
235//
236// particle decay related initialization
237 fDecayer->SetForceDecay(fForceDecay);
238 fDecayer->Init();
239
240//
241 AliGenMC::Init();
242}
243
244//____________________________________________________________
245void AliGenParam::Generate()
246{
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.
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))
256//
257//
258// Reinitialize decayer
259 fDecayer->SetForceDecay(fForceDecay);
260 fDecayer->Init();
261
262//
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 time0; // Time0 of the generated parent particle
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],
269 och[3]; // Momentum, polarisation and origin of the children particles from lujet
270 Double_t ty, xmt;
271 Int_t nt, i, j;
272 Float_t wgtp, wgtch;
273 Double_t dummy;
274 static TClonesArray *particles;
275 //
276 if(!particles) particles = new TClonesArray("TParticle",1000);
277
278 TDatabasePDG *pDataBase = TDatabasePDG::Instance();
279 //
280 Float_t random[6];
281
282// Calculating vertex position per event
283 for (j=0;j<3;j++) origin0[j]=fOrigin[j];
284 time0 = fTimeOrigin;
285 if(fVertexSmear==kPerEvent) {
286 Vertex();
287 for (j=0;j<3;j++) origin0[j]=fVertex[j];
288 time0 = fTime;
289 }
290
291 Int_t ipa=0;
292
293// Generating fNpart particles
294 fNprimaries = 0;
295
296 while (ipa<fNpart) {
297 while(1) {
298//
299// particle type
300 Int_t iPart = fIpParaFunc(fRandom);
301 fChildWeight=(fDecayer->GetPartialBranchingRatio(iPart))*fParentWeight;
302 TParticlePDG *particle = pDataBase->GetParticle(iPart);
303 Float_t am = particle->Mass();
304
305 Rndm(random,2);
306//
307// phi
308 phi=fPhiMin+random[0]*(fPhiMax-fPhiMin);
309//
310// y
311 ty = TMath::TanH(fYPara->GetRandom());
312//
313// pT
314 if (fAnalog == kAnalog) {
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);
325 if (TMath::Abs(ty)==1.) {
326 ty=0.;
327 Fatal("AliGenParam",
328 "Division by 0: Please check you rapidity range !");
329 }
330
331 pl=xmt*ty/sqrt((1.-ty)*(1.+ty));
332 theta=TMath::ATan2(pt,pl);
333// Cut on theta
334 if(theta<fThetaMin || theta>fThetaMax) continue;
335 ptot=TMath::Sqrt(pt*pt+pl*pl);
336// Cut on momentum
337 if(ptot<fPMin || ptot>fPMax) continue;
338//
339 p[0]=pt*TMath::Cos(phi);
340 p[1]=pt*TMath::Sin(phi);
341 p[2]=pl;
342 if(fVertexSmear==kPerTrack) {
343 Rndm(random,6);
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 }
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]));
353 }
354
355// Looking at fForceDecay :
356// if fForceDecay != none Primary particle decays using
357// AliPythia and children are tracked by GEANT
358//
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)
361//
362 Bool_t decayed = kFALSE;
363
364
365 if (fForceDecay != kNoDecay) {
366// Using lujet to decay particle
367 Float_t energy=TMath::Sqrt(ptot*ptot+am*am);
368 TLorentzVector pmom(p[0], p[1], p[2], energy);
369 fDecayer->Decay(iPart,&pmom);
370//
371// select decay particles
372 Int_t np=fDecayer->ImportParticles(particles);
373
374 // Selecting GeometryAcceptance for particles fPdgCodeParticleforAcceptanceCut;
375 if (fGeometryAcceptance)
376 if (!CheckAcceptanceGeometry(np,particles)) continue;
377 Int_t ncsel=0;
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
389 if (np >1) {
390 decayed = kTRUE;
391 TParticle* iparticle = 0;
392 Int_t ipF, ipL;
393 for (i = 1; i<np ; i++) {
394 trackIt[i] = 1;
395 iparticle = (TParticle *) particles->At(i);
396 Int_t kf = iparticle->GetPdgCode();
397 Int_t ks = iparticle->GetStatusCode();
398// flagged particle
399
400 if (pFlag[i] == 1) {
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) {
410// TParticlePDG *particle = pDataBase->GetParticle(kf);
411
412 Double_t lifeTime = fDecayer->GetLifetime(kf);
413// Double_t mass = particle->Mass();
414// Double_t width = particle->Width();
415 if (lifeTime > (Double_t) fMaxLifeTime) {
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 ?
424//
425// children
426
427 if ((ChildSelected(TMath::Abs(kf)) || fForceDecay == kAll || fSelectAll) && trackIt[i])
428 {
429 if (fCutOnChild) {
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;
436 ncsel++;
437 } else {
438 ncsel=-1;
439 break;
440 } // child kine cuts
441 } else {
442 pSelected[i] = 1;
443 ncsel++;
444 } // if child selection
445 } // select muon
446 } // decay particle loop
447 } // if decay products
448
449 Int_t iparent;
450 if ((fCutOnChild && ncsel >0) || !fCutOnChild){
451 ipa++;
452//
453// Parent
454
455
456 PushTrack(0, -1, iPart, p, origin0, polar, time0, kPPrimary, nt, wgtp, ((decayed)? 11 : 1));
457 pParent[0] = nt;
458 KeepTrack(nt);
459 fNprimaries++;
460
461//
462// Decay Products
463//
464 for (i = 1; i < np; i++) {
465 if (pSelected[i]) {
466 TParticle* iparticle = (TParticle *) particles->At(i);
467 Int_t kf = iparticle->GetPdgCode();
468 Int_t ksc = iparticle->GetStatusCode();
469 Int_t jpa = iparticle->GetFirstMother()-1;
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
478 if (jpa > -1) {
479 iparent = pParent[jpa];
480 } else {
481 iparent = -1;
482 }
483
484 PushTrack(fTrackIt * trackIt[i], iparent, kf,
485 pc, och, polar,
486 time0 + iparticle->T(), kPDecay, nt, wgtch, ksc);
487 pParent[i] = nt;
488 KeepTrack(nt);
489 fNprimaries++;
490 } // Selected
491 } // Particle loop
492 } // Decays by Lujet
493 particles->Clear();
494 if (pFlag) delete[] pFlag;
495 if (pParent) delete[] pParent;
496 if (pSelected) delete[] pSelected;
497 if (trackIt) delete[] trackIt;
498 } // kinematic selection
499 else // nodecay option, so parent will be tracked by GEANT (pions, kaons, eta, omegas, baryons)
500 {
501 gAlice->GetMCApp()->
502 PushTrack(fTrackIt,-1,iPart,p,origin0,polar,time0,kPPrimary,nt,wgtp, 1);
503 ipa++;
504 fNprimaries++;
505 }
506 break;
507 } // while
508 } // event loop
509
510 SetHighWaterMark(nt);
511
512 AliGenEventHeader* header = new AliGenEventHeader("PARAM");
513 header->SetPrimaryVertex(fVertex);
514 header->SetInteractionTime(fTime);
515 header->SetNProduced(fNprimaries);
516 AddHeader(header);
517}
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 =
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) *
527 (phiMax-phiMin)/360.;
528 return TMath::Abs(ratio);
529}
530
531//____________________________________________________________________________________
532
533void AliGenParam::Draw( const char * /*opt*/)
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
548
549
550