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