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