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