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