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