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