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