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