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