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