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