Updated version of the non-recursive Makefiles. See doc/README.FlatMakefile and build...
[u/mrichter/AliRoot.git] / EVGEN / AliGenHBTprocessor.cxx
CommitLineData
2b9786f4 1/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
2 * See cxx source for full Copyright notice */
3
4/* $Id$ */
5
6// Implementation of the interface for THBTprocessor
7// Author: Piotr Krzysztof Skowronski <Piotr.Skowronski@cern.ch>
8
9// 27.09.2001 Piotr Skowronski
10// removing of redefinition of defaults velues
11// in method's implementation.
12//
13
14#include "AliGenHBTprocessor.h"
15#include "TROOT.h"
16#include <iostream.h>
17
18#include "AliRun.h"
19#include "AliStack.h"
20#include "TParticle.h"
21#include "AliGenCocktailAfterBurner.h"
22
23
24
25const Int_t AliGenHBTprocessor::fgkHBTPMAXPART = 100000;
26
27ClassImp(AliGenHBTprocessor)
28
29//Wrapper class for "hbt processor" after burner
30//The origibal code is written in fortran by Lanny Ray
31//and is put in the directory $ALICE_ROOT/HBTprocessor
32//Detailed description is on the top of the file hbt_event_processor.f
33//
34//This class can be used ONLY with AliGenCocktailAfterBurner wrapper generator
35//i.e. (in Config.C)
36// ....
37// AliGenCocktailAfterBurner = gener = new AliGenCocktailAfterBurner();
38// gener->SetPhiRange(0, 360); //Set global parameters
39// gener->SetThetaRange(35., 145.);
40// AliGenHIJINGpara *hijing = new AliGenHIJINGpara(10000); //Add some generator
41// hijing->SetMomentumRange(0, 999);
42// gener->AddGenerator(hijing,"HIJING PARAMETRIZATION",1);
43//
44// AliGenHBTprocessor *hbtp = new AliGenHBTprocessor(); //create object
45// hbtp->SetRefControl(2); //set parameters
46// hbtp->SetSwitch1D(1);
47// hbtp->SetR0(6);//fm - radius
48// hbtp->SetLambda(0.7);//chaocity parameter
49// gener->AddAfterBurner(hbtp,"HBT PROCESSOR",1); //add to the main generator
50//
51// gener->Init();
52//
53//CAUTIONS:
54// A) HBT PROCESSOR NEEDS MORE THAN ONE EVENT TO WORK
55// AS MORE AS IT BETTER WORKS
56// B) IT IS ABLE TO "ADD" CORRELATIONS ONLY UP TO TWO PARTICLE TYPES AT ONES
57
58
59AliGenCocktailAfterBurner* GetGenerator();
60/*******************************************************************/
61
62AliGenHBTprocessor::AliGenHBTprocessor() : AliGenerator(-1)
63{
64 //
65 // Standard constructor
66 // Sets default veues of all parameters
67 SetName("AliGenHBTprocessor");
68 SetTitle("AliGenHBTprocessor");
69
70 sRandom = fRandom;
71 fHBTprocessor = new THBTprocessor();
72
73 fNPDGCodes = 0;
74 DefineParticles();
75
76 SetTrackRejectionFactor();
77 SetRefControl();
78 SetPIDs();
79 SetNPIDtypes();
80 SetDeltap();
81 SetMaxIterations();
82 SetDelChi();
83 SetIRand();
84 SetLambda();
85 SetR1d() ;
86 SetRSide();
87 SetROut() ;
88 SetRLong() ;
89 SetRPerp();
90 SetRParallel();
91 SetR0();
92 SetQ0();
93 SetSwitch1D();
94 SetSwitch3D();
95 SetSwitchType();
96 SetSwitchCoherence();
97 SetSwitchCoulomb();
98 SetSwitchFermiBose();
99 //SetMomentumRange();
100 SetPtRange();
101 SetPxRange();
102 SetPyRange();
103 SetPzRange();
104 SetPhiRange();
105 SetEtaRange();
106 SetNPtBins();
107 SetNPhiBins();
108 SetNEtaBins();
109 SetNPxBins();
110 SetNPyBins();
111 SetNPzBins();
112 SetNBins1DFineMesh();
113 SetBinSize1DFineMesh();
114 SetNBins1DCoarseMesh();
115 SetBinSize1DCoarseMesh();
116 SetNBins3DFineMesh();
117 SetBinSize3DFineMesh();
118 SetNBins3DCoarseMesh();
119 SetBinSize3DCoarseMesh();
120 SetNBins3DFineProjectMesh();
121}
122
123/*******************************************************************/
124
125
126/*******************************************************************/
127
128AliGenHBTprocessor::~AliGenHBTprocessor()
129{
130//destructor
131 CleanStatusCodes();
132 if (fHBTprocessor) delete fHBTprocessor; //delete generator
133
134}
135
136/*******************************************************************/
137
138void AliGenHBTprocessor::InitStatusCodes()
139{
140 //creates and inits status codes array to zero
141 AliGenCocktailAfterBurner *cab = GetGenerator();
142
143 if(!cab) Fatal("AliGenHBTprocessor::InitStatusCodes()","Can not find AliGenCocktailAfterBurner generator");
144
145 Int_t nev = cab->GetNumberOfEvents();
146
147 fHbtPStatCodes = new Int_t* [nev];
148 for( Int_t i =0; i<nev;i++)
149 {
150 Int_t nprim = cab->GetStack(i)->GetNprimary();
151 fHbtPStatCodes[i] = new Int_t[nprim];
152 for (Int_t k =0 ;k<nprim;k++)
153 fHbtPStatCodes[i][k] =0;
154
155 }
156
157}
158/*******************************************************************/
159
160void AliGenHBTprocessor::CleanStatusCodes()
161{//Cleans up status codes
162 if (fHbtPStatCodes)
163 {
164 for (Int_t i =0; i<GetGenerator()->GetNumberOfEvents(); i++)
165 delete [] fHbtPStatCodes[i];
166 delete fHbtPStatCodes;
167 fHbtPStatCodes = 0;
168 }
169
170}
171/*******************************************************************/
172
173void AliGenHBTprocessor::Init()
174 {
175 //sets up parameters in generator
176 THBTprocessor *thbtp = fHBTprocessor;
177
178 thbtp->SetTrackRejectionFactor(fTrackRejectionFactor);
179 thbtp->SetRefControl(fReferenceControl);
180
181 if ((fPid[0] == fPid[1]) || (fPid[0] == 0) || (fPid[1] == 0))
182 {
183 if (fPid[0] == 0)
184 thbtp->SetPIDs(IdFromPDG(fPid[1]) ,0);
185 else
186 thbtp->SetPIDs(IdFromPDG(fPid[0]) ,0);
187 thbtp->SetNPIDtypes(1);
188
189 if (fSwitch_type !=1)
190 Warning("AliGenHBTprocessor::Init","\nThere is only one particle type set,\n\
191 and Switch_Type differnt then 1,\n which does not make sense.\n\
192 Setting it to 1.\n");
193
194 SetSwitchType(1);
195 }
196 else
197 {
198 thbtp->SetPIDs(IdFromPDG(fPid[0]) ,IdFromPDG(fPid[1]));
199 SetNPIDtypes(2);
200 thbtp->SetSwitchType(fSwitch_type);
201 }
202
203 thbtp->SetMaxIterations(fMaxit);
204 thbtp->SetDelChi(fDelchi);
205 thbtp->SetIRand(fIrand);
206 thbtp->SetLambda(fLambda);
207 thbtp->SetR1d(fR_1d);
208 thbtp->SetRSide(fRside);
209 thbtp->SetROut(fRout);
210 thbtp->SetRLong(fRlong);
211 thbtp->SetRPerp(fRperp);
212 thbtp->SetRParallel(fRparallel);
213 thbtp->SetR0(fR0);
214 thbtp->SetQ0(fQ0);
215 thbtp->SetSwitch1D(fSwitch_1d);
216 thbtp->SetSwitch3D(fSwitch_3d);
217 thbtp->SetSwitchType(fSwitch_type);
218 thbtp->SetSwitchCoherence(fSwitch_coherence);
219 thbtp->SetSwitchCoulomb(fSwitch_coulomb);
220 thbtp->SetSwitchFermiBose(fSwitch_fermi_bose);
221 thbtp->SetPtRange(fPtMin,fPtMax);
222 thbtp->SetPxRange(fPx_min,fPx_max);
223 thbtp->SetPyRange(fPy_min,fPy_max);
224 thbtp->SetPzRange(fPz_min,fPz_max);
225 thbtp->SetPhiRange(fPhiMin*180/TMath::Pi(),fPhiMax*180/TMath::Pi());
226// cout<<"@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n";
227// cout<<fPhiMin<<fPhiMax<<endl;
228// cout<<"@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n";
229 thbtp->SetEtaRange(fEta_min,fEta_max);
230 thbtp->SetNPtBins(fN_pt_bins);
231 thbtp->SetNPhiBins(fN_phi_bins);
232 thbtp->SetNEtaBins(fN_eta_bins);
233 thbtp->SetNPxBins(fN_px_bins);
234 thbtp->SetNPyBins(fN_py_bins);
235 thbtp->SetNPzBins(fN_pz_bins);
236 thbtp->SetNBins1DFineMesh(fN_1d_fine);
237 thbtp->SetBinSize1DFineMesh(fBinsize_1d_fine);
238 thbtp->SetNBins1DCoarseMesh(fN_1d_coarse);
239 thbtp->SetBinSize1DCoarseMesh(fBinsize_1d_coarse);
240 thbtp->SetNBins3DFineMesh(fN_3d_fine);
241 thbtp->SetBinSize3DFineMesh(fBinsize_3d_fine);
242 thbtp->SetNBins3DCoarseMesh(fN_3d_coarse);
243 thbtp->SetBinSize3DCoarseMesh(fBinsize_3d_coarse);
244 thbtp->SetNBins3DFineProjectMesh(fN_3d_fine_project);
245
246 }
247/*******************************************************************/
248
249void AliGenHBTprocessor::Generate()
250 {
251 //starts processig
252
253 if (gAlice->GetEventsPerRun() <2)
254 {
255 Fatal("AliGenHBTprocessor::Generate()",
256 "HBT Processor needs more than 1 event to work properly,\
257 but there is only %d set", gAlice->GetEventsPerRun());
258 }
259
260 InitStatusCodes(); //Init status codes
261
262 fHBTprocessor->GenerateEvent(); //Generates event
263
264 CleanStatusCodes(); //Clean Status codes - thet are not needed anymore
265 }
266
267/*******************************************************************/
268
269
270/*******************************************************************/
271void AliGenHBTprocessor::GetParticles(TClonesArray * particles)
272 {//practically dumm
273 if (!particles)
274 {
275 cout<<"Particles has to be initialized"<<endl;
276 return;
277 }
278 fHBTprocessor->ImportParticles(particles);
279 }
280
281/*******************************************************************/
282
283Int_t AliGenHBTprocessor::GetHbtPStatusCode(Int_t part) const
284{
285//returns the status code of the given particle in the active event
286//see SetActiveEvent in the bottom of AliGenHBTprocessor.cxx
287//and in AliCocktailAfterBurner
288 Int_t ActiveEvent = GetGenerator()->GetActiveEventNumber();
289 return fHbtPStatCodes[ActiveEvent][part];
290}
291
292/*******************************************************************/
293void AliGenHBTprocessor::SetHbtPStatusCode(Int_t hbtstatcode, Int_t part)
294{
295 //Sets the given status code to given particle number (part) in the active event
296 Int_t ActiveEvent = GetGenerator()->GetActiveEventNumber();
297 fHbtPStatCodes[ActiveEvent][part] = hbtstatcode;
298}
299
300/*******************************************************************/
301
302void AliGenHBTprocessor::SetTrackRejectionFactor(Float_t trf) //def 1.0
303 {
304 //Sets the Track Rejection Factor
305 //variates in range 0.0 <-> 1.0
306 //Describes the factor of particles rejected from the output.
307 //Used only in case of low muliplicity particles e.g. lambdas.
308 //Processor generates addisional particles and builds the
309 //correletions on such a statistics.
310 //At the end these particels are left in the event according
311 //to this factor: 1==all particles are left
312 // 0==all are removed
313
314 fTrackRejectionFactor=trf;
315 fHBTprocessor->SetTrackRejectionFactor(trf);
316 }
317/*******************************************************************/
318
319void AliGenHBTprocessor::SetRefControl(Int_t rc) //default 2
320 {
321 //Sets the Refernce Control Switch
322 //switch wether read reference histograms from file =1
323 // compute from input events =2 - default
324 fReferenceControl=rc;
325 fHBTprocessor->SetRefControl(rc);
326 }
327/*******************************************************************/
328
329void AliGenHBTprocessor::SetPIDs(Int_t pid1,Int_t pid2)
330 {
331 //default pi+ pi-
332 //Sets PDG Codes of particles to be processed, default \\Pi^{+} and \\Pi^{-}
333 //This method accepts PDG codes which is
334 //in opposite to THBBProcessor which accepts GEANT PID
335 if ( (pid1 == 0) && (pid2 == 0) )
336 {
337 Error("AliGenHBTprocessor::SetPIDs","Sensless Particle Types setting: 0 0, Ignoring\n");
338 }
339
340 fPid[0]=pid1;
341 fPid[1]=pid2;
342
343 if(pid1 == pid2)
344 {
345 fHBTprocessor->SetPIDs(IdFromPDG(pid1) ,0);
346 SetNPIDtypes(1);
347 SetSwitchType(1);
348 }
349 else
350 {
351 fHBTprocessor->SetPIDs(IdFromPDG(pid1) ,IdFromPDG(pid2));
352 SetNPIDtypes(2);
353 }
354 }
355/*******************************************************************/
356
357void AliGenHBTprocessor::SetNPIDtypes(Int_t npidt)
358 {
359 //Number ofparticle types to be processed - default 2
360 //see AliGenHBTprocessor::SetPIDs
361
362 fNPidTypes = npidt;
363 fHBTprocessor->SetNPIDtypes(npidt);
364 }
365/*******************************************************************/
366
367void AliGenHBTprocessor::SetDeltap(Float_t deltp)
368 {
369 //default = 0.1 GeV
370 //maximum range for random momentum shifts in GeV/c;
371 //px,py,pz independent; Default = 0.1 GeV/c.
372 fDeltap=deltp;
373 fHBTprocessor->SetDeltap(deltp);
374 }
375/*******************************************************************/
376
377void AliGenHBTprocessor::SetMaxIterations(Int_t maxiter)
378 {
379 //maximum # allowed iterations thru all the
380 //tracks for each event. Default = 50.
381 //If maxit=0, then calculate the correlations
382 //for the input set of events without doing the
383 //track adjustment procedure.
384
385 fMaxit=maxiter;
386 fHBTprocessor->SetMaxIterations(maxiter);
387 }
388
389/*******************************************************************/
390void AliGenHBTprocessor::SetDelChi(Float_t dc)
391 {
392 //min % change in total chi-square for which
393 //the track adjustment iterations may stop,
394 //Default = 0.1%.
395
396 fDelchi=dc;
397 fHBTprocessor->SetDelChi(dc);
398 }
399/*******************************************************************/
400
401void AliGenHBTprocessor::SetIRand(Int_t irnd)
402 {
403 //if fact dummy - used only for compatibility
404 //we are using AliRoot built in RNG
405 //dummy in fact since we are using aliroot build-in RNG
406 //Sets renaodom number generator
407 fIrand=irnd;
408 fHBTprocessor->SetIRand(irnd);
409 }
410/*******************************************************************/
411
412void AliGenHBTprocessor::SetLambda(Float_t lam)
413 {
414 //default = 0.6
415 // Sets Chaoticity Parameter
416 fLambda = lam;
417 fHBTprocessor->SetLambda(lam);
418 }
419/*******************************************************************/
420void AliGenHBTprocessor::SetR1d(Float_t r)
421 {
422 //default 7.0
423 //Sets Spherical source model radius (fm)
424 fR_1d = r;
425 fHBTprocessor->SetR1d(r);
426 }
427/*******************************************************************/
428void AliGenHBTprocessor::SetRSide(Float_t rs)
429 {
430 //default rs = 6.0
431 //Rside,Rout,Rlong = Non-spherical Bertsch-Pratt source model (fm)
432
433 fRside = rs;
434 fHBTprocessor->SetRSide(rs);
435 }
436/*******************************************************************/
437void AliGenHBTprocessor::SetROut(Float_t ro)
438 {
439 //default ro = 7.0
440 //Rside,Rout,Rlong = Non-spherical Bertsch-Pratt source model (fm)
441 fRout = ro;
442 fHBTprocessor->SetROut(ro);
443 }
444/*******************************************************************/
445void AliGenHBTprocessor::SetRLong(Float_t rl)
446 {
447 //default rl = 4.0
448 //Rside,Rout,Rlong = Non-spherical Bertsch-Pratt source model (fm)
449 fRlong = rl;
450 fHBTprocessor->SetRLong(rl);
451 }
452/*******************************************************************/
453void AliGenHBTprocessor::SetRPerp(Float_t rp)
454 {
455 //default rp = 6.0
456 //Rperp,Rparallel,R0= Non-spherical Yano-Koonin-Podgoretski source model (fm).
457 fRperp = rp;
458 fHBTprocessor->SetRPerp(rp);
459 }
460/*******************************************************************/
461void AliGenHBTprocessor::SetRParallel(Float_t rprl)
462 {
463 //default rprl = 4.0
464 //Rperp,Rparallel,R0= Non-spherical Yano-Koonin-Podgoretski source model (fm).
465 fRparallel = rprl;
466 fHBTprocessor->SetRParallel(rprl);
467 }
468/*******************************************************************/
469void AliGenHBTprocessor::SetR0(Float_t r0)
470 {
471 //default r0 = 4.0
472 //Rperp,Rparallel,R0= Non-spherical Yano-Koonin-Podgoretski source model (fm).
473 fR0 = r0;
474 fHBTprocessor->SetR0(r0);
475 }
476/*******************************************************************/
477void AliGenHBTprocessor::SetQ0(Float_t q0)
478 {
479 //default q0 = 9.0
480 //Sets Q0 = NA35 Coulomb parameter for finite source size in (GeV/c)
481 // if fSwitchCoulomb = 2
482 // = Spherical Coulomb source radius in (fm)
483 // if switch_coulomb = 3, used to interpolate the
484 // input Pratt/Cramer discrete Coulomb source
485 // radii tables.
486 fQ0 = q0;
487 fHBTprocessor->SetQ0(q0);
488 }
489
490/*******************************************************************/
491void AliGenHBTprocessor::SetSwitch1D(Int_t s1d)
492 {
493//default s1d = 3
494// Sets fSwitch_3d
495// = 0 to not compute the 3D two-body correlations.
496// = 1 to compute this using the side-out-long form
497// = 2 to compute this using the Yanno-Koonin-Pogoredskij form.
498
499 fSwitch_1d = s1d;
500 fHBTprocessor->SetSwitch1D(s1d);
501 }
502/*******************************************************************/
503void AliGenHBTprocessor::SetSwitch3D(Int_t s3d)
504 {
505//default s3d = 0
506// Sets fSwitch_3d
507// = 0 to not compute the 1D two-body //orrelations.
508// = 1 to compute this using Q-invariant
509// = 2 to compute this using Q-total
510// = 3 to compute this using Q-3-ve//tor
511
512 fSwitch_3d = s3d;
513 fHBTprocessor->SetSwitch3D(s3d);
514 }
515/*******************************************************************/
516void AliGenHBTprocessor::SetSwitchType(Int_t st)
517 {
518//default st = 3
519// Sets switch_type = 1 to fit only the like pair correlations
520// = 2 to fit only the unlike pair correlations
521// = 3 to fit both the like and unlike pair correl.
522//See SetPIDs and Init
523//If only one particle type is set, unly==1 makes sens
524
525 fSwitch_type = st;
526 fHBTprocessor->SetSwitchType(st);
527 }
528/*******************************************************************/
529void AliGenHBTprocessor::SetSwitchCoherence(Int_t sc)
530 {
531// default sc = 0
532// switch_coherence = 0 for purely incoherent source (but can have
533// lambda < 1.0)
534// = 1 for mixed incoherent and coherent source
535
536 fSwitch_coherence = sc;
537 fHBTprocessor->SetSwitchCoherence(sc);
538 }
539/*******************************************************************/
540void AliGenHBTprocessor::SetSwitchCoulomb(Int_t scol)
541 {
542//default scol = 2
543// switch_coulomb = 0 no Coulomb correction
544// = 1 Point source Gamow correction only
545// = 2 NA35 finite source size correction
546// = 3 Pratt/Cramer finite source size correction;
547// interpolated from input tables.
548 fSwitch_coulomb =scol;
549 fHBTprocessor->SetSwitchCoulomb(scol);
550 }
551/*******************************************************************/
552void AliGenHBTprocessor::SetSwitchFermiBose(Int_t sfb)
553 {
554//default sfb = 1
555// switch_fermi_bose = 1 Boson pairs
556// = -1 Fermion pairs
557
558 fSwitch_fermi_bose = sfb;
559 fHBTprocessor->SetSwitchFermiBose(sfb);
560 }
561/*******************************************************************/
562void AliGenHBTprocessor::SetPtRange(Float_t ptmin, Float_t ptmax)
563 {
564// default ptmin = 0.1, ptmax = 0.98
565//Sets Pt range (GeV)
566 AliGenerator::SetPtRange(ptmin,ptmax);
567 fHBTprocessor->SetPtRange(ptmin,ptmax);
568 }
569
570/*******************************************************************/
571void AliGenHBTprocessor::SetPxRange(Float_t pxmin, Float_t pxmax)
572 {
573//default pxmin = -1.0, pxmax = 1.0
574//Sets Px range
575 fPx_min =pxmin;
576 fPx_max =pxmax;
577 fHBTprocessor->SetPxRange(pxmin,pxmax);
578 }
579/*******************************************************************/
580void AliGenHBTprocessor::SetPyRange(Float_t pymin, Float_t pymax)
581 {
582//default pymin = -1.0, pymax = 1.0
583//Sets Py range
584 fPy_min =pymin;
585 fPy_max =pymax;
586 fHBTprocessor->SetPyRange(pymin,pymax);
587 }
588/*******************************************************************/
589void AliGenHBTprocessor::SetPzRange(Float_t pzmin, Float_t pzmax)
590 {
591//default pzmin = -3.6, pzmax = 3.6
592//Sets Py range
593 fPz_min =pzmin;
594 fPz_max =pzmax;
595 fHBTprocessor->SetPzRange(pzmin,pzmax);
596 }
597void AliGenHBTprocessor::SetMomentumRange(Float_t pmin, Float_t pmax)
598 {
599 //default pmin=0, pmax=0
600 //Do not use this method!
601 MayNotUse("AliGenHBTprocessor::SetMomentumRange Method is Dummy");
602 }
603
604 /*******************************************************************/
605void AliGenHBTprocessor::SetPhiRange(Float_t phimin, Float_t phimax)
606 {
607//
608//Sets \\Phi range
609 AliGenerator::SetPhiRange(phimin,phimax);
610
611 fHBTprocessor->SetPhiRange(phimin,phimax);
612 }
613/*******************************************************************/
614void AliGenHBTprocessor::SetEtaRange(Float_t etamin, Float_t etamax)
615 {
616//default etamin = -1.5, etamax = 1.5
617//Sets \\Eta range
618 fEta_min= etamin;
619 fEta_max =etamax;
620 fHBTprocessor->SetEtaRange(etamin,etamax);
621 }
622/*******************************************************************/
623void AliGenHBTprocessor::SetNPtBins(Int_t nptbin)
624 {
625 //default nptbin = 50
626 //set number of Pt bins
627 fN_pt_bins= nptbin;
628 fHBTprocessor->SetNPtBins(nptbin);
629 }
630/*******************************************************************/
631void AliGenHBTprocessor::SetNPhiBins(Int_t nphibin)
632{
633 //default nphibin = 50
634 //set number of Phi bins
635 fN_phi_bins=nphibin;
636 fHBTprocessor->SetNPhiBins(nphibin);
637}
638/*******************************************************************/
639void AliGenHBTprocessor::SetNEtaBins(Int_t netabin)
640{
641 //default netabin = 50
642 //set number of Eta bins
643 fN_eta_bins = netabin;
644 fHBTprocessor->SetNEtaBins(netabin);
645}
646/*******************************************************************/
647void AliGenHBTprocessor::SetNPxBins(Int_t npxbin)
648{
649 //default npxbin = 20
650 //set number of Px bins
651 fN_px_bins = npxbin;
652 fHBTprocessor->SetNPxBins(npxbin);
653}
654/*******************************************************************/
655void AliGenHBTprocessor::SetNPyBins(Int_t npybin)
656{
657 //default npybin = 20
658 //set number of Py bins
659 fN_py_bins = npybin;
660 fHBTprocessor->SetNPyBins(npybin);
661}
662/*******************************************************************/
663void AliGenHBTprocessor::SetNPzBins(Int_t npzbin)
664{
665 //default npzbin = 70
666 //set number of Pz bins
667 fN_pz_bins = npzbin;
668 fHBTprocessor->SetNPzBins(npzbin);
669}
670/*******************************************************************/
671void AliGenHBTprocessor::SetNBins1DFineMesh(Int_t n)
672{
673//default n = 10
674//Sets the number of bins in the 1D mesh
675 fN_1d_fine =n;
676 fHBTprocessor->SetNBins1DFineMesh(n);
677
678}
679/*******************************************************************/
680void AliGenHBTprocessor::SetBinSize1DFineMesh(Float_t x)
681{
682//default x=0.01
683//Sets the bin size in the 1D mesh
684 fBinsize_1d_fine = x;
685 fHBTprocessor->SetBinSize1DFineMesh(x);
686}
687/*******************************************************************/
688
689void AliGenHBTprocessor::SetNBins1DCoarseMesh(Int_t n)
690{
691//default n =2
692//Sets the number of bins in the coarse 1D mesh
693 fN_1d_coarse =n;
694 fHBTprocessor->SetNBins1DCoarseMesh(n);
695}
696/*******************************************************************/
697void AliGenHBTprocessor::SetBinSize1DCoarseMesh(Float_t x)
698{
699//default x=0.05
700//Sets the bin size in the coarse 1D mesh
701 fBinsize_1d_coarse =x;
702 fHBTprocessor->SetBinSize1DCoarseMesh(x);
703}
704/*******************************************************************/
705
706void AliGenHBTprocessor::SetNBins3DFineMesh(Int_t n)
707{
708//default n = 8
709//Sets the number of bins in the 3D mesh
710 fN_3d_fine =n;
711 fHBTprocessor->SetNBins3DFineMesh(n);
712}
713/*******************************************************************/
714void AliGenHBTprocessor::SetBinSize3DFineMesh(Float_t x)
715{
716//default x=0.01
717//Sets the bin size in the 3D mesh
718 fBinsize_3d_fine =x;
719 fHBTprocessor->SetBinSize3DFineMesh(x);
720}
721/*******************************************************************/
722
723void AliGenHBTprocessor::SetNBins3DCoarseMesh(Int_t n)
724{
725//default n = 2
726//Sets the number of bins in the coarse 3D mesh
727
728 fN_3d_coarse = n;
729 fHBTprocessor->SetNBins3DCoarseMesh(n);
730}
731/*******************************************************************/
732void AliGenHBTprocessor::SetBinSize3DCoarseMesh(Float_t x)
733{
734//default x=0.08
735//Sets the bin size in the coarse 3D mesh
736 fBinsize_3d_coarse = x;
737 fHBTprocessor->SetBinSize3DCoarseMesh(x);
738}
739/*******************************************************************/
740
741void AliGenHBTprocessor::SetNBins3DFineProjectMesh(Int_t n )
742{
743//default n =3
744//Sets the number of bins in the fine project mesh
745 fN_3d_fine_project = n;
746 fHBTprocessor->SetNBins3DFineProjectMesh(n);
747}
748/*******************************************************************/
749
750
751/*******************************************************************/
752
753
754
755
756
757
758void AliGenHBTprocessor::DefineParticles()
759{
760 //
761 // Load standard numbers for GEANT particles and PDG conversion
762 fNPDGCodes = 0; //this is done in the constructor - but in any case ...
763
764 fPDGCode[fNPDGCodes++]=-99; // 0 = unused location
765 fPDGCode[fNPDGCodes++]=22; // 1 = photon
766 fPDGCode[fNPDGCodes++]=-11; // 2 = positron
767 fPDGCode[fNPDGCodes++]=11; // 3 = electron
768 fPDGCode[fNPDGCodes++]=12; // 4 = neutrino e
769 fPDGCode[fNPDGCodes++]=-13; // 5 = muon +
770 fPDGCode[fNPDGCodes++]=13; // 6 = muon -
771 fPDGCode[fNPDGCodes++]=111; // 7 = pi0
772 fPDGCode[fNPDGCodes++]=211; // 8 = pi+
773 fPDGCode[fNPDGCodes++]=-211; // 9 = pi-
774 fPDGCode[fNPDGCodes++]=130; // 10 = Kaon Long
775 fPDGCode[fNPDGCodes++]=321; // 11 = Kaon +
776 fPDGCode[fNPDGCodes++]=-321; // 12 = Kaon -
777 fPDGCode[fNPDGCodes++]=2112; // 13 = Neutron
778 fPDGCode[fNPDGCodes++]=2212; // 14 = Proton
779 fPDGCode[fNPDGCodes++]=-2212; // 15 = Anti Proton
780 fPDGCode[fNPDGCodes++]=310; // 16 = Kaon Short
781 fPDGCode[fNPDGCodes++]=221; // 17 = Eta
782 fPDGCode[fNPDGCodes++]=3122; // 18 = Lambda
783 fPDGCode[fNPDGCodes++]=3222; // 19 = Sigma +
784 fPDGCode[fNPDGCodes++]=3212; // 20 = Sigma 0
785 fPDGCode[fNPDGCodes++]=3112; // 21 = Sigma -
786 fPDGCode[fNPDGCodes++]=3322; // 22 = Xi0
787 fPDGCode[fNPDGCodes++]=3312; // 23 = Xi-
788 fPDGCode[fNPDGCodes++]=3334; // 24 = Omega-
789 fPDGCode[fNPDGCodes++]=-2112; // 25 = Anti Proton
790 fPDGCode[fNPDGCodes++]=-3122; // 26 = Anti Proton
791 fPDGCode[fNPDGCodes++]=-3222; // 27 = Anti Sigma -
792 fPDGCode[fNPDGCodes++]=-3212; // 28 = Anti Sigma 0
793 fPDGCode[fNPDGCodes++]=-3112; // 29 = Anti Sigma 0
794 fPDGCode[fNPDGCodes++]=-3322; // 30 = Anti Xi 0
795 fPDGCode[fNPDGCodes++]=-3312; // 31 = Anti Xi +
796 fPDGCode[fNPDGCodes++]=-3334; // 32 = Anti Omega +
797}
798
799/*******************************************************************/
800Int_t AliGenHBTprocessor::IdFromPDG(Int_t pdg) const
801{
802 //
803 // Return Geant3 code from PDG and pseudo ENDF code
804 //
805 for(Int_t i=0;i<fNPDGCodes;++i)
806 if(pdg==fPDGCode[i]) return i;
807 return -1;
808}
809Int_t AliGenHBTprocessor::PDGFromId(Int_t id) const
810{
811 //
812 // Return PDG code and pseudo ENDF code from Geant3 code
813 //
814 if(id>0 && id<fNPDGCodes) return fPDGCode[id];
815 else return -1;
816}
817/*******************************************************************/
818/****** ROUTINES USED FOR COMMUNUCATION ********/
819/******************** WITH FORTRAN ********************/
820/*******************************************************************/
821
822#ifndef WIN32
823 # define hbtpran hbtpran_
824 # define alihbtp_puttrack alihbtp_puttrack_
825 # define alihbtp_gettrack alihbtp_gettrack_
826 # define alihbtp_getnumberevents alihbtp_getnumberevents_
827 # define alihbtp_getnumbertracks alihbtp_getnumbertracks_
828 # define alihbtp_initialize alihbtp_initialize_
829 # define alihbtp_setactiveeventnumber alihbtp_setactiveeventnumber_
830 # define alihbtp_setparameters alihbtp_setparameters_
831 # define type_of_call
832
833#else
834 # define hbtpran HBTPRAN
835 # define alihbtp_puttrack ALIHBTP_PUTTRACK
836 # define alihbtp_gettrack ALIHBTP_GETTRACK
837 # define alihbtp_getnumberevents ALIHBTP_GETNUMBEREVENTS
838 # define alihbtp_getnumbertracks ALIHBTP_GETNUMBERTRACKS
839 # define alihbtp_initialize ALIHBTP_INITIALIZE
840 # define alihbtp_setactiveeventnumber ALIHBTP_SETACTIVEEVENTNUMBER
841 # define alihbtp_setparameters ALIHBTP_SETPARAMETERS
842 # define type_of_call _stdcall
843#endif
844
845#include "AliGenCocktailAfterBurner.h"
846#include <string.h>
847/*******************************************************************/
848
849AliGenCocktailAfterBurner* GetGenerator()
850 {
851 // This function has two tasks:
852 // Check if environment is OK (exist gAlice and generator)
853 // Returns pointer to genrator
854 //to be changed with TFolders
855
856 if(!gAlice)
857 {
858 cout<<endl<<"ERROR: There is NO gALICE! Check what you are doing!"<<endl;
859 gAlice->Fatal("AliGenHBTprocessor.cxx: alihbtp_getnumofev()",
860 "\nRunning HBT Processor without gAlice... Exiting \n");
861 return 0;
862 }
863 AliGenerator * gen = gAlice->Generator();
864
865 if (!gen)
866 {
867 gAlice->Fatal("AliGenHBTprocessor.cxx: GetGenerator()",
868 "\nThere is no generator in gAlice, exiting\n");
869 return 0;
870 }
871
872 const Char_t *genname = gen->GetName();
873 Char_t name[30];
874 strcpy(name,"AliGenCocktailAfterBurner");
875
876 if (strcmp(name,genname))
877 {
878 gAlice->Fatal("AliGenHBTprocessor.cxx: GetGenerator()",
879 "\nThe main Generator is not a AliGenCocktailAfterBurner, exiting\n");
880 return 0;
881 }
882
883 AliGenCocktailAfterBurner* CAB= (AliGenCocktailAfterBurner*) gen;
884
885 // cout<<endl<<"Got generator"<<endl;
886 return CAB;
887
888 }
889/*******************************************************************/
890
891AliGenHBTprocessor* GetAliGenHBTprocessor()
892{
893//returns pointer to the current instance of AliGenHBTprocessor in
894//AliGenCocktailAfterBurner (can be more than one)
895//
896 AliGenCocktailAfterBurner* gen = GetGenerator();
897 AliGenerator * g = gen->GetCurrentGenerator();
898
899 const Char_t *genname = g->GetName();
900 Char_t name[30];
901 strcpy(name,"AliGenHBTprocessor");
902
903 if (strcmp(name,genname))
904 {
905 gAlice->Fatal("AliGenHBTprocessor.cxx: GetAliGenHBTprocessor()",
906 "\nCurrernt generator in AliGenCocktailAfterBurner is not a AliGenHBTprocessor, exiting\n");
907 return 0;
908 }
909
910 AliGenHBTprocessor* hbtp = (AliGenHBTprocessor*)g;
911 return hbtp;
912}
913
914/*******************************************************************/
915extern "C" void type_of_call alihbtp_setparameters()
916 {
917 //dummy
918 }
919
920extern "C" void type_of_call alihbtp_initialize()
921 {
922 //dummy
923 }
924
925/*******************************************************************/
926
927extern "C" void type_of_call alihbtp_getnumberevents(Int_t &nev)
928 {
929 //passes number of events to the fortran
930 if(gDebug) cout<<"alihbtp_getnumberevents("<<nev<<") ....";
931 AliGenCocktailAfterBurner* gen = GetGenerator();
932 if(!gen)
933 {
934 nev = -1;
935 return;
936 }
937
938 nev = gen->GetNumberOfEvents();
939
940 if(gDebug>5) cout<<"EXITED N Ev = "<<nev<<endl;
941
942 }
943
944/*******************************************************************/
945
946extern "C" void type_of_call alihbtp_setactiveeventnumber(Int_t & nev)
947 {
948//sets active event in generator (AliGenCocktailAfterBurner)
949
950 if(gDebug>5) cout<<"alihbtp_setactiveeventnumber("<<nev<<") ....";
951 if(gDebug>0) cout<<"Asked for event "<<nev-1<<endl;
952 AliGenCocktailAfterBurner* gen = GetGenerator();
953 if(!gen) return;
954 gen->SetActiveEventNumber(nev - 1); //fortran numerates events from 1 to N
955
956 if(gDebug>5) cout<<"EXITED returned "<<nev<<endl;
957
958 }
959/*******************************************************************/
960
961extern "C" void type_of_call alihbtp_getnumbertracks(Int_t &ntracks)
962 {
963//passes number of particles in active event to the fortran
964 if(gDebug>5) cout<<"alihbtp_getnumbertracks("<<ntracks<<") ....";
965
966 AliGenCocktailAfterBurner* gen = GetGenerator();
967 if (!gen)
968 {
969 ntracks = -1;
970 return;
971 }
972
973 ntracks = gen->GetActiveStack()->GetNprimary();
974 if(gDebug>5) cout<<"EXITED Ntracks = "<<ntracks<<endl;
975 }
976
977/*******************************************************************/
978
979extern "C" void type_of_call
980 alihbtp_puttrack(Int_t & n,Int_t& flag, Float_t& px,
981 Float_t& py, Float_t& pz, Int_t& geantpid)
982 {
983//sets new parameters (momenta) in track number n
984// in the active event
985// n - number of the track in active event
986// flag - flag of the track
987// px,py,pz - momenta
988// geantpid - type of the particle - Geant Particle ID
989
990 if(gDebug>5) cout<<"alihbtp_puttrack("<<n<<") ....";
991
992 AliGenCocktailAfterBurner* gen = GetGenerator();
993 if(!gen) return;
994
995 TParticle * track = gen->GetActiveStack()->Particle(n-1);
996
997 AliGenHBTprocessor* g = GetAliGenHBTprocessor();
998
999 //check to be deleted
1000 if (geantpid != (g->IdFromPDG( track->GetPdgCode() )))
1001 {
1002 cout<<endl<<" AliGenHBTprocessor.cxx: alihbtp_puttrack: SOMETHING IS GOING BAD:\n GEANTPIDS ARE NOT THE SAME"<<endl;
1003 }
1004
1005 if(gDebug>0)
1006 if (px != track->Px())
1007 cout<<"Px diff. = "<<px - track->Px()<<endl;
1008
1009 if(gDebug>3) cout<<" track->GetPdgCode() --> "<<track->GetPdgCode()<<endl;
1010
1011
1012
5925692f 1013 Float_t m = track->GetMass();
1014 Float_t E = TMath::Sqrt(m*m+px*px+py*py+pz*pz);
1015 track->SetMomentum(px,py,pz,E);
2b9786f4 1016
1017 g->SetHbtPStatusCode(flag,n-1);
1018
1019 if(gDebug>5) cout<<"EXITED "<<endl;
1020 }
1021
1022/*******************************************************************/
1023
1024extern "C" void type_of_call
1025 alihbtp_gettrack(Int_t & n,Int_t & flag, Float_t & px,
1026 Float_t & py, Float_t & pz, Int_t & geantpid)
1027
1028 {
1029//passes track parameters to the fortran
1030// n - number of the track in active event
1031// flag - flag of the track
1032// px,py,pz - momenta
1033// geantpid - type of the particle - Geant Particle ID
1034
1035 if(gDebug>5) cout<<"alihbtp_gettrack("<<n<<") ....";
1036 AliGenCocktailAfterBurner* gen = GetGenerator();
1037
1038 if (!gen)
1039 {
1040 n = -1;
1041 flag =-1;
1042 px = py = pz = -1;
1043 geantpid = -1;
1044 return;
1045 }
1046
1047 TParticle * track = gen->GetActiveStack()->Particle(n-1);
1048 AliGenHBTprocessor* g = GetAliGenHBTprocessor();
1049
1050 flag = g->GetHbtPStatusCode(n-1);
1051
1052 px = (Float_t)track->Px();
1053 py = (Float_t)track->Py();
1054 pz = (Float_t)track->Pz();
1055
1056 geantpid = g->IdFromPDG( track->GetPdgCode() );
1057
1058 if(gDebug>5) cout<<"EXITED "<<endl;
1059 }
1060
1061/*******************************************************************/
1062extern "C" Float_t type_of_call hbtpran(Int_t &)
1063{
1064//interface to the random number generator
1065 return sRandom->Rndm();
1066}
1067
1068/*******************************************************************/
1069
1070
1071/*******************************************************************/