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