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