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