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