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