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