Add const to SolenoidField() method for correct overwriting from the parent class
[u/mrichter/AliRoot.git] / EVGEN / AliGenHBTprocessor.cxx
CommitLineData
2b9786f4 1/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
2 * See cxx source for full Copyright notice */
2b9786f4 3// Implementation of the interface for THBTprocessor
4// Author: Piotr Krzysztof Skowronski <Piotr.Skowronski@cern.ch>
4c90f2a0 5//////////////////////////////////////////////////////////////////////////////////
2b9786f4 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
4c90f2a0 34//////////////////////////////////////////////////////////////////////////////////
35
ea6982d2 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
4c90f2a0 57#include "AliGenHBTprocessor.h"
58#include "TROOT.h"
59#include <iostream.h>
20dddfab 60#include <TFile.h>
61#include <TTree.h>
4c90f2a0 62#include "AliRun.h"
63#include "AliStack.h"
64#include "TParticle.h"
20dddfab 65#include "THBTprocessor.h"
4c90f2a0 66#include "AliGenCocktailAfterBurner.h"
67
68
69
4c90f2a0 70ClassImp(AliGenHBTprocessor)
71
2b9786f4 72
73
74AliGenCocktailAfterBurner* GetGenerator();
75/*******************************************************************/
20dddfab 76AliGenHBTprocessor::AliGenHBTprocessor(const AliGenHBTprocessor& in)
77{
78//copy contructor
635be086 79 // AliGenHBTprocessor::AliGenHBTprocessor();
20dddfab 80}
2b9786f4 81
84b7985c 82AliGenHBTprocessor::AliGenHBTprocessor() : AliGenerator()
2b9786f4 83{
84 //
85 // Standard constructor
86 // Sets default veues of all parameters
20dddfab 87 fHbtPStatCodes = 0x0;
88 fHBTprocessor = 0x0;
89
2b9786f4 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
151AliGenHBTprocessor::~AliGenHBTprocessor()
152{
153//destructor
154 CleanStatusCodes();
155 if (fHBTprocessor) delete fHBTprocessor; //delete generator
156
157}
158
159/*******************************************************************/
160
161void AliGenHBTprocessor::InitStatusCodes()
162{
163 //creates and inits status codes array to zero
164 AliGenCocktailAfterBurner *cab = GetGenerator();
165
84b7985c 166 if(!cab) Fatal("InitStatusCodes()","Can not find AliGenCocktailAfterBurner generator");
2b9786f4 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
183void AliGenHBTprocessor::CleanStatusCodes()
20dddfab 184{
185 //Cleans up status codes
2b9786f4 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
197void AliGenHBTprocessor::Init()
198 {
199 //sets up parameters in generator
84b7985c 200
2b9786f4 201 THBTprocessor *thbtp = fHBTprocessor;
202
4c90f2a0 203
2b9786f4 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
20dddfab 215 if (fSwitchType !=1)
2b9786f4 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);
20dddfab 226 thbtp->SetSwitchType(fSwitchType);
2b9786f4 227 }
228
229 thbtp->SetMaxIterations(fMaxit);
230 thbtp->SetDelChi(fDelchi);
231 thbtp->SetIRand(fIrand);
232 thbtp->SetLambda(fLambda);
20dddfab 233 thbtp->SetR1d(fR1d);
2b9786f4 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);
20dddfab 241 thbtp->SetSwitch1D(fSwitch1d);
242 thbtp->SetSwitch3D(fSwitch3d);
243 thbtp->SetSwitchType(fSwitchType);
244 thbtp->SetSwitchCoherence(fSwitchCoherence);
245 thbtp->SetSwitchCoulomb(fSwitchCoulomb);
246 thbtp->SetSwitchFermiBose(fSwitchFermiBose);
2b9786f4 247 thbtp->SetPtRange(fPtMin,fPtMax);
20dddfab 248 thbtp->SetPxRange(fPxMin,fPxMax);
249 thbtp->SetPyRange(fPyMin,fPyMax);
250 thbtp->SetPzRange(fPzMin,fPzMax);
4c90f2a0 251 thbtp->SetPhiRange(fPhiMin*180./TMath::Pi(),fPhiMax*180./TMath::Pi());
20dddfab 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);
2b9786f4 268
269 }
270/*******************************************************************/
271
272void AliGenHBTprocessor::Generate()
273 {
274 //starts processig
84b7985c 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)
2b9786f4 281 {
84b7985c 282 Fatal("Generate()",
2b9786f4 283 "HBT Processor needs more than 1 event to work properly,\
84b7985c 284 but there is only %d set", cab->GetNumberOfEvents());
2b9786f4 285 }
4c90f2a0 286
84b7985c 287
4c90f2a0 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
2b9786f4 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/*******************************************************************/
ff6aba1e 305void AliGenHBTprocessor::GetParticles(TClonesArray * particles) const
306 {
307 //practically dumm
2b9786f4 308 if (!particles)
309 {
310 cout<<"Particles has to be initialized"<<endl;
311 return;
312 }
313 fHBTprocessor->ImportParticles(particles);
314 }
315
316/*******************************************************************/
317
318Int_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
20dddfab 323 Int_t activeEvent = GetGenerator()->GetActiveEventNumber();
324 return fHbtPStatCodes[activeEvent][part];
2b9786f4 325}
326
327/*******************************************************************/
328void AliGenHBTprocessor::SetHbtPStatusCode(Int_t hbtstatcode, Int_t part)
329{
330 //Sets the given status code to given particle number (part) in the active event
20dddfab 331 Int_t activeEvent = GetGenerator()->GetActiveEventNumber();
332 fHbtPStatCodes[activeEvent][part] = hbtstatcode;
2b9786f4 333}
334
335/*******************************************************************/
336
337void 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
354void 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
364void 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);
4c90f2a0 381 SetNPIDtypes(1);
382 SetSwitchType(1);
2b9786f4 383 }
384 else
385 {
386 fHBTprocessor->SetPIDs(IdFromPDG(pid1) ,IdFromPDG(pid2));
4c90f2a0 387 SetNPIDtypes(2);
2b9786f4 388 }
389 }
390/*******************************************************************/
391
392void 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
402void 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
412void 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/*******************************************************************/
425void 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
436void 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
447void AliGenHBTprocessor::SetLambda(Float_t lam)
448 {
449 //default = 0.6
450 // Sets Chaoticity Parameter
451 fLambda = lam;
452 fHBTprocessor->SetLambda(lam);
453 }
454/*******************************************************************/
455void AliGenHBTprocessor::SetR1d(Float_t r)
456 {
457 //default 7.0
458 //Sets Spherical source model radius (fm)
20dddfab 459 fR1d = r;
2b9786f4 460 fHBTprocessor->SetR1d(r);
461 }
462/*******************************************************************/
463void 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/*******************************************************************/
472void 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/*******************************************************************/
480void 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/*******************************************************************/
488void 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/*******************************************************************/
496void 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/*******************************************************************/
504void 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/*******************************************************************/
512void 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)
20dddfab 518 // if switchCoulomb = 3, used to interpolate the
2b9786f4 519 // input Pratt/Cramer discrete Coulomb source
520 // radii tables.
521 fQ0 = q0;
522 fHBTprocessor->SetQ0(q0);
523 }
524
525/*******************************************************************/
526void AliGenHBTprocessor::SetSwitch1D(Int_t s1d)
527 {
528//default s1d = 3
20dddfab 529// Sets fSwitch1d
4c90f2a0 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
2b9786f4 534
20dddfab 535 fSwitch1d = s1d;
2b9786f4 536 fHBTprocessor->SetSwitch1D(s1d);
537 }
538/*******************************************************************/
539void AliGenHBTprocessor::SetSwitch3D(Int_t s3d)
540 {
541//default s3d = 0
20dddfab 542// Sets fSwitch3d
4c90f2a0 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.
2b9786f4 546
20dddfab 547 fSwitch3d = s3d;
2b9786f4 548 fHBTprocessor->SetSwitch3D(s3d);
549 }
550/*******************************************************************/
551void 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
20dddfab 560 fSwitchType = st;
2b9786f4 561 fHBTprocessor->SetSwitchType(st);
562 }
563/*******************************************************************/
564void AliGenHBTprocessor::SetSwitchCoherence(Int_t sc)
565 {
566// default sc = 0
20dddfab 567// switchCoherence = 0 for purely incoherent source (but can have
2b9786f4 568// lambda < 1.0)
569// = 1 for mixed incoherent and coherent source
570
20dddfab 571 fSwitchCoherence = sc;
2b9786f4 572 fHBTprocessor->SetSwitchCoherence(sc);
573 }
574/*******************************************************************/
575void AliGenHBTprocessor::SetSwitchCoulomb(Int_t scol)
576 {
577//default scol = 2
20dddfab 578// switchCoulomb = 0 no Coulomb correction
2b9786f4 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.
20dddfab 583 fSwitchCoulomb =scol;
2b9786f4 584 fHBTprocessor->SetSwitchCoulomb(scol);
585 }
586/*******************************************************************/
587void AliGenHBTprocessor::SetSwitchFermiBose(Int_t sfb)
588 {
589//default sfb = 1
20dddfab 590// switchFermiBose = 1 Boson pairs
2b9786f4 591// = -1 Fermion pairs
592
20dddfab 593 fSwitchFermiBose = sfb;
2b9786f4 594 fHBTprocessor->SetSwitchFermiBose(sfb);
595 }
596/*******************************************************************/
597void 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/*******************************************************************/
606void AliGenHBTprocessor::SetPxRange(Float_t pxmin, Float_t pxmax)
607 {
608//default pxmin = -1.0, pxmax = 1.0
609//Sets Px range
20dddfab 610 fPxMin =pxmin;
611 fPxMax =pxmax;
2b9786f4 612 fHBTprocessor->SetPxRange(pxmin,pxmax);
613 }
614/*******************************************************************/
615void AliGenHBTprocessor::SetPyRange(Float_t pymin, Float_t pymax)
616 {
617//default pymin = -1.0, pymax = 1.0
618//Sets Py range
20dddfab 619 fPyMin =pymin;
620 fPyMax =pymax;
2b9786f4 621 fHBTprocessor->SetPyRange(pymin,pymax);
622 }
623/*******************************************************************/
624void AliGenHBTprocessor::SetPzRange(Float_t pzmin, Float_t pzmax)
625 {
626//default pzmin = -3.6, pzmax = 3.6
627//Sets Py range
20dddfab 628 fPzMin =pzmin;
629 fPzMax =pzmax;
2b9786f4 630 fHBTprocessor->SetPzRange(pzmin,pzmax);
631 }
632void 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 /*******************************************************************/
640void 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/*******************************************************************/
649void AliGenHBTprocessor::SetEtaRange(Float_t etamin, Float_t etamax)
650 {
651//default etamin = -1.5, etamax = 1.5
652//Sets \\Eta range
20dddfab 653 fEtaMin= etamin;
654 fEtaMax =etamax;
2b9786f4 655 fHBTprocessor->SetEtaRange(etamin,etamax);
4c90f2a0 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
20dddfab 661 AliGenerator::SetThetaRange(RadiansToDegrees(EtaToTheta(fEtaMin)),
662 RadiansToDegrees(EtaToTheta(fEtaMax)));
4c90f2a0 663
2b9786f4 664 }
665/*******************************************************************/
4c90f2a0 666void 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/*******************************************************************/
2b9786f4 680void AliGenHBTprocessor::SetNPtBins(Int_t nptbin)
681 {
682 //default nptbin = 50
683 //set number of Pt bins
20dddfab 684 fNPtBins= nptbin;
2b9786f4 685 fHBTprocessor->SetNPtBins(nptbin);
686 }
687/*******************************************************************/
688void AliGenHBTprocessor::SetNPhiBins(Int_t nphibin)
689{
690 //default nphibin = 50
691 //set number of Phi bins
20dddfab 692 fNPhiBins=nphibin;
2b9786f4 693 fHBTprocessor->SetNPhiBins(nphibin);
694}
695/*******************************************************************/
696void AliGenHBTprocessor::SetNEtaBins(Int_t netabin)
697{
698 //default netabin = 50
699 //set number of Eta bins
20dddfab 700 fNEtaBins = netabin;
2b9786f4 701 fHBTprocessor->SetNEtaBins(netabin);
702}
703/*******************************************************************/
704void AliGenHBTprocessor::SetNPxBins(Int_t npxbin)
705{
706 //default npxbin = 20
707 //set number of Px bins
20dddfab 708 fNPxBins = npxbin;
2b9786f4 709 fHBTprocessor->SetNPxBins(npxbin);
710}
711/*******************************************************************/
712void AliGenHBTprocessor::SetNPyBins(Int_t npybin)
713{
714 //default npybin = 20
715 //set number of Py bins
20dddfab 716 fNPyBins = npybin;
2b9786f4 717 fHBTprocessor->SetNPyBins(npybin);
718}
719/*******************************************************************/
720void AliGenHBTprocessor::SetNPzBins(Int_t npzbin)
721{
722 //default npzbin = 70
723 //set number of Pz bins
20dddfab 724 fNPzBins = npzbin;
2b9786f4 725 fHBTprocessor->SetNPzBins(npzbin);
726}
727/*******************************************************************/
728void AliGenHBTprocessor::SetNBins1DFineMesh(Int_t n)
729{
730//default n = 10
731//Sets the number of bins in the 1D mesh
20dddfab 732 fN1dFine =n;
2b9786f4 733 fHBTprocessor->SetNBins1DFineMesh(n);
734
735}
736/*******************************************************************/
737void AliGenHBTprocessor::SetBinSize1DFineMesh(Float_t x)
738{
739//default x=0.01
740//Sets the bin size in the 1D mesh
20dddfab 741 fBinsize1dFine = x;
2b9786f4 742 fHBTprocessor->SetBinSize1DFineMesh(x);
743}
744/*******************************************************************/
745
746void AliGenHBTprocessor::SetNBins1DCoarseMesh(Int_t n)
747{
748//default n =2
749//Sets the number of bins in the coarse 1D mesh
20dddfab 750 fN1dCoarse =n;
2b9786f4 751 fHBTprocessor->SetNBins1DCoarseMesh(n);
752}
753/*******************************************************************/
754void AliGenHBTprocessor::SetBinSize1DCoarseMesh(Float_t x)
755{
756//default x=0.05
757//Sets the bin size in the coarse 1D mesh
20dddfab 758 fBinsize1dCoarse =x;
2b9786f4 759 fHBTprocessor->SetBinSize1DCoarseMesh(x);
760}
761/*******************************************************************/
762
763void AliGenHBTprocessor::SetNBins3DFineMesh(Int_t n)
764{
765//default n = 8
766//Sets the number of bins in the 3D mesh
20dddfab 767 fN3dFine =n;
2b9786f4 768 fHBTprocessor->SetNBins3DFineMesh(n);
769}
770/*******************************************************************/
771void AliGenHBTprocessor::SetBinSize3DFineMesh(Float_t x)
772{
773//default x=0.01
774//Sets the bin size in the 3D mesh
20dddfab 775 fBinsize3dFine =x;
2b9786f4 776 fHBTprocessor->SetBinSize3DFineMesh(x);
777}
778/*******************************************************************/
779
780void AliGenHBTprocessor::SetNBins3DCoarseMesh(Int_t n)
781{
782//default n = 2
783//Sets the number of bins in the coarse 3D mesh
784
20dddfab 785 fN3dCoarse = n;
2b9786f4 786 fHBTprocessor->SetNBins3DCoarseMesh(n);
787}
788/*******************************************************************/
789void AliGenHBTprocessor::SetBinSize3DCoarseMesh(Float_t x)
790{
791//default x=0.08
792//Sets the bin size in the coarse 3D mesh
20dddfab 793 fBinsize3dCoarse = x;
2b9786f4 794 fHBTprocessor->SetBinSize3DCoarseMesh(x);
795}
796/*******************************************************************/
797
798void AliGenHBTprocessor::SetNBins3DFineProjectMesh(Int_t n )
799{
800//default n =3
801//Sets the number of bins in the fine project mesh
20dddfab 802 fN3dFineProject = n;
2b9786f4 803 fHBTprocessor->SetNBins3DFineProjectMesh(n);
804}
805/*******************************************************************/
806
807
808/*******************************************************************/
809
810
811
812
813
814
815void 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/*******************************************************************/
857Int_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}
866Int_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;
4c90f2a0 873}
874Double_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
2b9786f4 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_
20dddfab 906 # define type_ofCall
2b9786f4 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
20dddfab 917 # define type_ofCall _stdcall
2b9786f4 918#endif
919
920#include "AliGenCocktailAfterBurner.h"
921#include <string.h>
922/*******************************************************************/
923
924AliGenCocktailAfterBurner* 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;
84b7985c 934 gROOT->Fatal("AliGenHBTprocessor.cxx: GetGenerator()",
2b9786f4 935 "\nRunning HBT Processor without gAlice... Exiting \n");
84b7985c 936 return 0x0;
2b9786f4 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");
84b7985c 944 return 0x0;
2b9786f4 945 }
84b7985c 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
2b9786f4 949
84b7985c 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
20dddfab 954 AliGenCocktailAfterBurner* cab=(AliGenCocktailAfterBurner*)genclass->DynamicCast(cabclass,gen);
84b7985c 955
20dddfab 956 if (cab == 0x0)//if generator that we got is not AliGenCocktailAfterBurner or its descendant we get null
84b7985c 957 { //then quit with error
2b9786f4 958 gAlice->Fatal("AliGenHBTprocessor.cxx: GetGenerator()",
84b7985c 959 "\nThe main Generator is not a AliGenCocktailAfterBurner, exiting\n");
960 return 0x0;
2b9786f4 961 }
2b9786f4 962 // cout<<endl<<"Got generator"<<endl;
20dddfab 963 return cab;
2b9786f4 964
965 }
966/*******************************************************************/
967
968AliGenHBTprocessor* GetAliGenHBTprocessor()
969{
970//returns pointer to the current instance of AliGenHBTprocessor in
971//AliGenCocktailAfterBurner (can be more than one)
972//
973 AliGenCocktailAfterBurner* gen = GetGenerator();
84b7985c 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)
2b9786f4 986 {
987 gAlice->Fatal("AliGenHBTprocessor.cxx: GetAliGenHBTprocessor()",
988 "\nCurrernt generator in AliGenCocktailAfterBurner is not a AliGenHBTprocessor, exiting\n");
84b7985c 989 return 0x0;
2b9786f4 990 }
2b9786f4 991 return hbtp;
992}
993
994/*******************************************************************/
20dddfab 995extern "C" void type_ofCall alihbtp_setparameters()
2b9786f4 996 {
997 //dummy
998 }
999
20dddfab 1000extern "C" void type_ofCall alihbtp_initialize()
2b9786f4 1001 {
1002 //dummy
1003 }
1004
1005/*******************************************************************/
1006
20dddfab 1007extern "C" void type_ofCall alihbtp_getnumberevents(Int_t &nev)
2b9786f4 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
20dddfab 1026extern "C" void type_ofCall alihbtp_setactiveeventnumber(Int_t & nev)
2b9786f4 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
20dddfab 1041extern "C" void type_ofCall alihbtp_getnumbertracks(Int_t &ntracks)
2b9786f4 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
20dddfab 1059extern "C" void type_ofCall
2b9786f4 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 {
ea6982d2 1082 cerr<<endl<<" AliGenHBTprocessor.cxx: alihbtp_puttrack: SOMETHING IS GOING BAD:\n GEANTPIDS ARE NOT THE SAME"<<endl;
2b9786f4 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
5925692f 1093 Float_t m = track->GetMass();
20dddfab 1094 Float_t e = TMath::Sqrt(m*m+px*px+py*py+pz*pz);
1095 track->SetMomentum(px,py,pz,e);
2b9786f4 1096
1097 g->SetHbtPStatusCode(flag,n-1);
1098
1099 if(gDebug>5) cout<<"EXITED "<<endl;
1100 }
1101
1102/*******************************************************************/
1103
20dddfab 1104extern "C" void type_ofCall
2b9786f4 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/*******************************************************************/
20dddfab 1142extern "C" Float_t type_ofCall hbtpran(Int_t &)
2b9786f4 1143{
1144//interface to the random number generator
1145 return sRandom->Rndm();
1146}
1147
1148/*******************************************************************/
1149
1150
1151/*******************************************************************/