1 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
2 * See cxx source for full Copyright notice */
3 // Implementation of the interface for THBTprocessor
4 // Author: Piotr Krzysztof Skowronski <Piotr.Skowronski@cern.ch>
5 //////////////////////////////////////////////////////////////////////////////////
6 //Wrapper class for "hbt processor" after burner
7 //The origibal code is written in fortran by Lanny Ray
8 //and is put in the directory $ALICE_ROOT/HBTprocessor
9 //Detailed description is on the top of the file hbt_event_processor.f
11 //This class can be used ONLY with AliGenCocktailAfterBurner wrapper generator
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);
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
31 // A) HBT PROCESSOR NEEDS MORE THAN ONE EVENT TO WORK
32 // AS MORE AS IT BETTER WORKS
33 // B) IT IS ABLE TO "ADD" CORRELATIONS ONLY UP TO TWO PARTICLE TYPES AT ONES
34 //////////////////////////////////////////////////////////////////////////////////
36 // 11.11.2001 Piotr Skowronski
37 // Setting seed (date) in RNG in the constructor
39 // 09.10.2001 Piotr Skowronski
41 // Theta - Eta cohernecy facilities added:
42 // AliGenerator::SetThetaRange method overriden
43 // Static methods added
49 // Class description comments put on proper place
51 // 27.09.2001 Piotr Skowronski
52 // removing of redefinition of defaults velues
53 // in method's implementation.
57 #include "AliGenHBTprocessor.h"
64 #include "TParticle.h"
65 #include "THBTprocessor.h"
66 #include "AliGenCocktailAfterBurner.h"
70 ClassImp(AliGenHBTprocessor)
74 AliGenCocktailAfterBurner* GetGenerator();
75 /*******************************************************************/
76 AliGenHBTprocessor::AliGenHBTprocessor(const AliGenHBTprocessor& in)
79 AliGenHBTprocessor::AliGenHBTprocessor();
82 AliGenHBTprocessor::AliGenHBTprocessor() : AliGenerator()
85 // Standard constructor
86 // Sets default veues of all parameters
90 SetName("AliGenHBTprocessor");
91 SetTitle("AliGenHBTprocessor");
94 fHBTprocessor = new THBTprocessor();
99 SetTrackRejectionFactor();
119 SetSwitchCoherence();
121 SetSwitchFermiBose();
122 //SetMomentumRange();
135 SetNBins1DFineMesh();
136 SetBinSize1DFineMesh();
137 SetNBins1DCoarseMesh();
138 SetBinSize1DCoarseMesh();
139 SetNBins3DFineMesh();
140 SetBinSize3DFineMesh();
141 SetNBins3DCoarseMesh();
142 SetBinSize3DCoarseMesh();
143 SetNBins3DFineProjectMesh();
146 /*******************************************************************/
149 /*******************************************************************/
151 AliGenHBTprocessor::~AliGenHBTprocessor()
155 if (fHBTprocessor) delete fHBTprocessor; //delete generator
159 /*******************************************************************/
161 void AliGenHBTprocessor::InitStatusCodes()
163 //creates and inits status codes array to zero
164 AliGenCocktailAfterBurner *cab = GetGenerator();
166 if(!cab) Fatal("InitStatusCodes()","Can not find AliGenCocktailAfterBurner generator");
168 Int_t nev = cab->GetNumberOfEvents();
170 fHbtPStatCodes = new Int_t* [nev];
171 for( Int_t i =0; i<nev;i++)
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;
181 /*******************************************************************/
183 void AliGenHBTprocessor::CleanStatusCodes()
185 //Cleans up status codes
188 for (Int_t i =0; i<GetGenerator()->GetNumberOfEvents(); i++)
189 delete [] fHbtPStatCodes[i];
190 delete fHbtPStatCodes;
195 /*******************************************************************/
197 void AliGenHBTprocessor::Init()
199 //sets up parameters in generator
201 THBTprocessor *thbtp = fHBTprocessor;
204 thbtp->SetTrackRejectionFactor(fTrackRejectionFactor);
205 thbtp->SetRefControl(fReferenceControl);
207 if ((fPid[0] == fPid[1]) || (fPid[0] == 0) || (fPid[1] == 0))
210 thbtp->SetPIDs(IdFromPDG(fPid[1]) ,0);
212 thbtp->SetPIDs(IdFromPDG(fPid[0]) ,0);
213 thbtp->SetNPIDtypes(1);
216 Warning("AliGenHBTprocessor::Init","\nThere is only one particle type set,\n\
217 and Switch_Type differnt then 1,\n which does not make sense.\n\
218 Setting it to 1.\n");
224 thbtp->SetPIDs(IdFromPDG(fPid[0]) ,IdFromPDG(fPid[1]));
226 thbtp->SetSwitchType(fSwitchType);
229 thbtp->SetMaxIterations(fMaxit);
230 thbtp->SetDelChi(fDelchi);
231 thbtp->SetIRand(fIrand);
232 thbtp->SetLambda(fLambda);
234 thbtp->SetRSide(fRside);
235 thbtp->SetROut(fRout);
236 thbtp->SetRLong(fRlong);
237 thbtp->SetRPerp(fRperp);
238 thbtp->SetRParallel(fRparallel);
241 thbtp->SetSwitch1D(fSwitch1d);
242 thbtp->SetSwitch3D(fSwitch3d);
243 thbtp->SetSwitchType(fSwitchType);
244 thbtp->SetSwitchCoherence(fSwitchCoherence);
245 thbtp->SetSwitchCoulomb(fSwitchCoulomb);
246 thbtp->SetSwitchFermiBose(fSwitchFermiBose);
247 thbtp->SetPtRange(fPtMin,fPtMax);
248 thbtp->SetPxRange(fPxMin,fPxMax);
249 thbtp->SetPyRange(fPyMin,fPyMax);
250 thbtp->SetPzRange(fPzMin,fPzMax);
251 thbtp->SetPhiRange(fPhiMin*180./TMath::Pi(),fPhiMax*180./TMath::Pi());
252 thbtp->SetEtaRange(fEtaMin,fEtaMax);
253 thbtp->SetNPtBins(fNPtBins);
254 thbtp->SetNPhiBins(fNPhiBins);
255 thbtp->SetNEtaBins(fNEtaBins);
256 thbtp->SetNPxBins(fNPxBins);
257 thbtp->SetNPyBins(fNPyBins);
258 thbtp->SetNPzBins(fNPzBins);
259 thbtp->SetNBins1DFineMesh(fN1dFine);
260 thbtp->SetBinSize1DFineMesh(fBinsize1dFine);
261 thbtp->SetNBins1DCoarseMesh(fN1dCoarse);
262 thbtp->SetBinSize1DCoarseMesh(fBinsize1dCoarse);
263 thbtp->SetNBins3DFineMesh(fN3dFine);
264 thbtp->SetBinSize3DFineMesh(fBinsize3dFine);
265 thbtp->SetNBins3DCoarseMesh(fN3dCoarse);
266 thbtp->SetBinSize3DCoarseMesh(fBinsize3dCoarse);
267 thbtp->SetNBins3DFineProjectMesh(fN3dFineProject);
270 /*******************************************************************/
272 void AliGenHBTprocessor::Generate()
275 AliGenCocktailAfterBurner* cab = GetGenerator();
278 Fatal("Generate()","AliGenHBTprocessor needs AliGenCocktailAfterBurner to be main generator");
280 if (cab->GetNumberOfEvents() <2)
283 "HBT Processor needs more than 1 event to work properly,\
284 but there is only %d set", cab->GetNumberOfEvents());
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
294 InitStatusCodes(); //Init status codes
296 fHBTprocessor->GenerateEvent(); //Generates event
298 CleanStatusCodes(); //Clean Status codes - thet are not needed anymore
301 /*******************************************************************/
304 /*******************************************************************/
305 void AliGenHBTprocessor::GetParticles(TClonesArray * particles)
309 cout<<"Particles has to be initialized"<<endl;
312 fHBTprocessor->ImportParticles(particles);
315 /*******************************************************************/
317 Int_t AliGenHBTprocessor::GetHbtPStatusCode(Int_t part) const
319 //returns the status code of the given particle in the active event
320 //see SetActiveEvent in the bottom of AliGenHBTprocessor.cxx
321 //and in AliCocktailAfterBurner
322 Int_t activeEvent = GetGenerator()->GetActiveEventNumber();
323 return fHbtPStatCodes[activeEvent][part];
326 /*******************************************************************/
327 void AliGenHBTprocessor::SetHbtPStatusCode(Int_t hbtstatcode, Int_t part)
329 //Sets the given status code to given particle number (part) in the active event
330 Int_t activeEvent = GetGenerator()->GetActiveEventNumber();
331 fHbtPStatCodes[activeEvent][part] = hbtstatcode;
334 /*******************************************************************/
336 void AliGenHBTprocessor::SetTrackRejectionFactor(Float_t trf) //def 1.0
338 //Sets the Track Rejection Factor
339 //variates in range 0.0 <-> 1.0
340 //Describes the factor of particles rejected from the output.
341 //Used only in case of low muliplicity particles e.g. lambdas.
342 //Processor generates addisional particles and builds the
343 //correletions on such a statistics.
344 //At the end these particels are left in the event according
345 //to this factor: 1==all particles are left
346 // 0==all are removed
348 fTrackRejectionFactor=trf;
349 fHBTprocessor->SetTrackRejectionFactor(trf);
351 /*******************************************************************/
353 void AliGenHBTprocessor::SetRefControl(Int_t rc) //default 2
355 //Sets the Refernce Control Switch
356 //switch wether read reference histograms from file =1
357 // compute from input events =2 - default
358 fReferenceControl=rc;
359 fHBTprocessor->SetRefControl(rc);
361 /*******************************************************************/
363 void AliGenHBTprocessor::SetPIDs(Int_t pid1,Int_t pid2)
366 //Sets PDG Codes of particles to be processed, default \\Pi^{+} and \\Pi^{-}
367 //This method accepts PDG codes which is
368 //in opposite to THBBProcessor which accepts GEANT PID
369 if ( (pid1 == 0) && (pid2 == 0) )
371 Error("AliGenHBTprocessor::SetPIDs","Sensless Particle Types setting: 0 0, Ignoring\n");
379 fHBTprocessor->SetPIDs(IdFromPDG(pid1) ,0);
385 fHBTprocessor->SetPIDs(IdFromPDG(pid1) ,IdFromPDG(pid2));
389 /*******************************************************************/
391 void AliGenHBTprocessor::SetNPIDtypes(Int_t npidt)
393 //Number ofparticle types to be processed - default 2
394 //see AliGenHBTprocessor::SetPIDs
397 fHBTprocessor->SetNPIDtypes(npidt);
399 /*******************************************************************/
401 void AliGenHBTprocessor::SetDeltap(Float_t deltp)
404 //maximum range for random momentum shifts in GeV/c;
405 //px,py,pz independent; Default = 0.1 GeV/c.
407 fHBTprocessor->SetDeltap(deltp);
409 /*******************************************************************/
411 void AliGenHBTprocessor::SetMaxIterations(Int_t maxiter)
413 //maximum # allowed iterations thru all the
414 //tracks for each event. Default = 50.
415 //If maxit=0, then calculate the correlations
416 //for the input set of events without doing the
417 //track adjustment procedure.
420 fHBTprocessor->SetMaxIterations(maxiter);
423 /*******************************************************************/
424 void AliGenHBTprocessor::SetDelChi(Float_t dc)
426 //min % change in total chi-square for which
427 //the track adjustment iterations may stop,
431 fHBTprocessor->SetDelChi(dc);
433 /*******************************************************************/
435 void AliGenHBTprocessor::SetIRand(Int_t irnd)
437 //if fact dummy - used only for compatibility
438 //we are using AliRoot built in RNG
439 //dummy in fact since we are using aliroot build-in RNG
440 //Sets renaodom number generator
442 fHBTprocessor->SetIRand(irnd);
444 /*******************************************************************/
446 void AliGenHBTprocessor::SetLambda(Float_t lam)
449 // Sets Chaoticity Parameter
451 fHBTprocessor->SetLambda(lam);
453 /*******************************************************************/
454 void AliGenHBTprocessor::SetR1d(Float_t r)
457 //Sets Spherical source model radius (fm)
459 fHBTprocessor->SetR1d(r);
461 /*******************************************************************/
462 void AliGenHBTprocessor::SetRSide(Float_t rs)
465 //Rside,Rout,Rlong = Non-spherical Bertsch-Pratt source model (fm)
468 fHBTprocessor->SetRSide(rs);
470 /*******************************************************************/
471 void AliGenHBTprocessor::SetROut(Float_t ro)
474 //Rside,Rout,Rlong = Non-spherical Bertsch-Pratt source model (fm)
476 fHBTprocessor->SetROut(ro);
478 /*******************************************************************/
479 void AliGenHBTprocessor::SetRLong(Float_t rl)
482 //Rside,Rout,Rlong = Non-spherical Bertsch-Pratt source model (fm)
484 fHBTprocessor->SetRLong(rl);
486 /*******************************************************************/
487 void AliGenHBTprocessor::SetRPerp(Float_t rp)
490 //Rperp,Rparallel,R0= Non-spherical Yano-Koonin-Podgoretski source model (fm).
492 fHBTprocessor->SetRPerp(rp);
494 /*******************************************************************/
495 void AliGenHBTprocessor::SetRParallel(Float_t rprl)
498 //Rperp,Rparallel,R0= Non-spherical Yano-Koonin-Podgoretski source model (fm).
500 fHBTprocessor->SetRParallel(rprl);
502 /*******************************************************************/
503 void AliGenHBTprocessor::SetR0(Float_t r0)
506 //Rperp,Rparallel,R0= Non-spherical Yano-Koonin-Podgoretski source model (fm).
508 fHBTprocessor->SetR0(r0);
510 /*******************************************************************/
511 void AliGenHBTprocessor::SetQ0(Float_t q0)
514 //Sets Q0 = NA35 Coulomb parameter for finite source size in (GeV/c)
515 // if fSwitchCoulomb = 2
516 // = Spherical Coulomb source radius in (fm)
517 // if switchCoulomb = 3, used to interpolate the
518 // input Pratt/Cramer discrete Coulomb source
521 fHBTprocessor->SetQ0(q0);
524 /*******************************************************************/
525 void AliGenHBTprocessor::SetSwitch1D(Int_t s1d)
529 // = 0 to not compute the 1D two-body //orrelations.
530 // = 1 to compute this using Q-invariant
531 // = 2 to compute this using Q-total
532 // = 3 to compute this using Q-3-ve//tor
535 fHBTprocessor->SetSwitch1D(s1d);
537 /*******************************************************************/
538 void AliGenHBTprocessor::SetSwitch3D(Int_t s3d)
542 // = 0 to not compute the 3D two-body correlations.
543 // = 1 to compute this using the side-out-long form
544 // = 2 to compute this using the Yanno-Koonin-Pogoredskij form.
547 fHBTprocessor->SetSwitch3D(s3d);
549 /*******************************************************************/
550 void AliGenHBTprocessor::SetSwitchType(Int_t st)
553 // Sets switch_type = 1 to fit only the like pair correlations
554 // = 2 to fit only the unlike pair correlations
555 // = 3 to fit both the like and unlike pair correl.
556 //See SetPIDs and Init
557 //If only one particle type is set, unly==1 makes sens
560 fHBTprocessor->SetSwitchType(st);
562 /*******************************************************************/
563 void AliGenHBTprocessor::SetSwitchCoherence(Int_t sc)
566 // switchCoherence = 0 for purely incoherent source (but can have
568 // = 1 for mixed incoherent and coherent source
570 fSwitchCoherence = sc;
571 fHBTprocessor->SetSwitchCoherence(sc);
573 /*******************************************************************/
574 void AliGenHBTprocessor::SetSwitchCoulomb(Int_t scol)
577 // switchCoulomb = 0 no Coulomb correction
578 // = 1 Point source Gamow correction only
579 // = 2 NA35 finite source size correction
580 // = 3 Pratt/Cramer finite source size correction;
581 // interpolated from input tables.
582 fSwitchCoulomb =scol;
583 fHBTprocessor->SetSwitchCoulomb(scol);
585 /*******************************************************************/
586 void AliGenHBTprocessor::SetSwitchFermiBose(Int_t sfb)
589 // switchFermiBose = 1 Boson pairs
590 // = -1 Fermion pairs
592 fSwitchFermiBose = sfb;
593 fHBTprocessor->SetSwitchFermiBose(sfb);
595 /*******************************************************************/
596 void AliGenHBTprocessor::SetPtRange(Float_t ptmin, Float_t ptmax)
598 // default ptmin = 0.1, ptmax = 0.98
599 //Sets Pt range (GeV)
600 AliGenerator::SetPtRange(ptmin,ptmax);
601 fHBTprocessor->SetPtRange(ptmin,ptmax);
604 /*******************************************************************/
605 void AliGenHBTprocessor::SetPxRange(Float_t pxmin, Float_t pxmax)
607 //default pxmin = -1.0, pxmax = 1.0
611 fHBTprocessor->SetPxRange(pxmin,pxmax);
613 /*******************************************************************/
614 void AliGenHBTprocessor::SetPyRange(Float_t pymin, Float_t pymax)
616 //default pymin = -1.0, pymax = 1.0
620 fHBTprocessor->SetPyRange(pymin,pymax);
622 /*******************************************************************/
623 void AliGenHBTprocessor::SetPzRange(Float_t pzmin, Float_t pzmax)
625 //default pzmin = -3.6, pzmax = 3.6
629 fHBTprocessor->SetPzRange(pzmin,pzmax);
631 void AliGenHBTprocessor::SetMomentumRange(Float_t pmin, Float_t pmax)
633 //default pmin=0, pmax=0
634 //Do not use this method!
635 MayNotUse("AliGenHBTprocessor::SetMomentumRange Method is Dummy");
638 /*******************************************************************/
639 void AliGenHBTprocessor::SetPhiRange(Float_t phimin, Float_t phimax)
643 AliGenerator::SetPhiRange(phimin,phimax);
645 fHBTprocessor->SetPhiRange(phimin,phimax);
647 /*******************************************************************/
648 void AliGenHBTprocessor::SetEtaRange(Float_t etamin, Float_t etamax)
650 //default etamin = -1.5, etamax = 1.5
654 fHBTprocessor->SetEtaRange(etamin,etamax);
656 //set the azimothal angle range in the AliGeneraor -
657 //to keep coherency between azimuthal angle and pseudorapidity
658 //DO NOT CALL this->SetThetaRange, because it calls this method (where we are)
659 //which must cause INFINITE LOOP
660 AliGenerator::SetThetaRange(RadiansToDegrees(EtaToTheta(fEtaMin)),
661 RadiansToDegrees(EtaToTheta(fEtaMax)));
664 /*******************************************************************/
665 void AliGenHBTprocessor::SetThetaRange(Float_t thetamin, Float_t thetamax)
667 //default thetamin = 0, thetamax = 180
668 //Azimuthal angle, override AliGenerator method which uses widely (i.e. wrapper generators)
669 //core fortran HBTProcessor uses Eta (pseudorapidity)
670 //so these methods has to be synchronized
672 AliGenerator::SetThetaRange(thetamin,thetamax);
674 SetEtaRange( ThetaToEta(fThetaMin) , ThetaToEta(fThetaMax) );
678 /*******************************************************************/
679 void AliGenHBTprocessor::SetNPtBins(Int_t nptbin)
681 //default nptbin = 50
682 //set number of Pt bins
684 fHBTprocessor->SetNPtBins(nptbin);
686 /*******************************************************************/
687 void AliGenHBTprocessor::SetNPhiBins(Int_t nphibin)
689 //default nphibin = 50
690 //set number of Phi bins
692 fHBTprocessor->SetNPhiBins(nphibin);
694 /*******************************************************************/
695 void AliGenHBTprocessor::SetNEtaBins(Int_t netabin)
697 //default netabin = 50
698 //set number of Eta bins
700 fHBTprocessor->SetNEtaBins(netabin);
702 /*******************************************************************/
703 void AliGenHBTprocessor::SetNPxBins(Int_t npxbin)
705 //default npxbin = 20
706 //set number of Px bins
708 fHBTprocessor->SetNPxBins(npxbin);
710 /*******************************************************************/
711 void AliGenHBTprocessor::SetNPyBins(Int_t npybin)
713 //default npybin = 20
714 //set number of Py bins
716 fHBTprocessor->SetNPyBins(npybin);
718 /*******************************************************************/
719 void AliGenHBTprocessor::SetNPzBins(Int_t npzbin)
721 //default npzbin = 70
722 //set number of Pz bins
724 fHBTprocessor->SetNPzBins(npzbin);
726 /*******************************************************************/
727 void AliGenHBTprocessor::SetNBins1DFineMesh(Int_t n)
730 //Sets the number of bins in the 1D mesh
732 fHBTprocessor->SetNBins1DFineMesh(n);
735 /*******************************************************************/
736 void AliGenHBTprocessor::SetBinSize1DFineMesh(Float_t x)
739 //Sets the bin size in the 1D mesh
741 fHBTprocessor->SetBinSize1DFineMesh(x);
743 /*******************************************************************/
745 void AliGenHBTprocessor::SetNBins1DCoarseMesh(Int_t n)
748 //Sets the number of bins in the coarse 1D mesh
750 fHBTprocessor->SetNBins1DCoarseMesh(n);
752 /*******************************************************************/
753 void AliGenHBTprocessor::SetBinSize1DCoarseMesh(Float_t x)
756 //Sets the bin size in the coarse 1D mesh
758 fHBTprocessor->SetBinSize1DCoarseMesh(x);
760 /*******************************************************************/
762 void AliGenHBTprocessor::SetNBins3DFineMesh(Int_t n)
765 //Sets the number of bins in the 3D mesh
767 fHBTprocessor->SetNBins3DFineMesh(n);
769 /*******************************************************************/
770 void AliGenHBTprocessor::SetBinSize3DFineMesh(Float_t x)
773 //Sets the bin size in the 3D mesh
775 fHBTprocessor->SetBinSize3DFineMesh(x);
777 /*******************************************************************/
779 void AliGenHBTprocessor::SetNBins3DCoarseMesh(Int_t n)
782 //Sets the number of bins in the coarse 3D mesh
785 fHBTprocessor->SetNBins3DCoarseMesh(n);
787 /*******************************************************************/
788 void AliGenHBTprocessor::SetBinSize3DCoarseMesh(Float_t x)
791 //Sets the bin size in the coarse 3D mesh
792 fBinsize3dCoarse = x;
793 fHBTprocessor->SetBinSize3DCoarseMesh(x);
795 /*******************************************************************/
797 void AliGenHBTprocessor::SetNBins3DFineProjectMesh(Int_t n )
800 //Sets the number of bins in the fine project mesh
802 fHBTprocessor->SetNBins3DFineProjectMesh(n);
804 /*******************************************************************/
807 /*******************************************************************/
814 void AliGenHBTprocessor::DefineParticles()
817 // Load standard numbers for GEANT particles and PDG conversion
818 fNPDGCodes = 0; //this is done in the constructor - but in any case ...
820 fPDGCode[fNPDGCodes++]=-99; // 0 = unused location
821 fPDGCode[fNPDGCodes++]=22; // 1 = photon
822 fPDGCode[fNPDGCodes++]=-11; // 2 = positron
823 fPDGCode[fNPDGCodes++]=11; // 3 = electron
824 fPDGCode[fNPDGCodes++]=12; // 4 = neutrino e
825 fPDGCode[fNPDGCodes++]=-13; // 5 = muon +
826 fPDGCode[fNPDGCodes++]=13; // 6 = muon -
827 fPDGCode[fNPDGCodes++]=111; // 7 = pi0
828 fPDGCode[fNPDGCodes++]=211; // 8 = pi+
829 fPDGCode[fNPDGCodes++]=-211; // 9 = pi-
830 fPDGCode[fNPDGCodes++]=130; // 10 = Kaon Long
831 fPDGCode[fNPDGCodes++]=321; // 11 = Kaon +
832 fPDGCode[fNPDGCodes++]=-321; // 12 = Kaon -
833 fPDGCode[fNPDGCodes++]=2112; // 13 = Neutron
834 fPDGCode[fNPDGCodes++]=2212; // 14 = Proton
835 fPDGCode[fNPDGCodes++]=-2212; // 15 = Anti Proton
836 fPDGCode[fNPDGCodes++]=310; // 16 = Kaon Short
837 fPDGCode[fNPDGCodes++]=221; // 17 = Eta
838 fPDGCode[fNPDGCodes++]=3122; // 18 = Lambda
839 fPDGCode[fNPDGCodes++]=3222; // 19 = Sigma +
840 fPDGCode[fNPDGCodes++]=3212; // 20 = Sigma 0
841 fPDGCode[fNPDGCodes++]=3112; // 21 = Sigma -
842 fPDGCode[fNPDGCodes++]=3322; // 22 = Xi0
843 fPDGCode[fNPDGCodes++]=3312; // 23 = Xi-
844 fPDGCode[fNPDGCodes++]=3334; // 24 = Omega-
845 fPDGCode[fNPDGCodes++]=-2112; // 25 = Anti Proton
846 fPDGCode[fNPDGCodes++]=-3122; // 26 = Anti Proton
847 fPDGCode[fNPDGCodes++]=-3222; // 27 = Anti Sigma -
848 fPDGCode[fNPDGCodes++]=-3212; // 28 = Anti Sigma 0
849 fPDGCode[fNPDGCodes++]=-3112; // 29 = Anti Sigma 0
850 fPDGCode[fNPDGCodes++]=-3322; // 30 = Anti Xi 0
851 fPDGCode[fNPDGCodes++]=-3312; // 31 = Anti Xi +
852 fPDGCode[fNPDGCodes++]=-3334; // 32 = Anti Omega +
855 /*******************************************************************/
856 Int_t AliGenHBTprocessor::IdFromPDG(Int_t pdg) const
859 // Return Geant3 code from PDG and pseudo ENDF code
861 for(Int_t i=0;i<fNPDGCodes;++i)
862 if(pdg==fPDGCode[i]) return i;
865 Int_t AliGenHBTprocessor::PDGFromId(Int_t id) const
868 // Return PDG code and pseudo ENDF code from Geant3 code
870 if(id>0 && id<fNPDGCodes) return fPDGCode[id];
873 Double_t AliGenHBTprocessor::ThetaToEta(Double_t arg)
875 //converts etha(azimuthal angle) to Eta (pseudorapidity). Argument in radians
877 if(arg>= TMath::Pi()) return 709.0;//This number is the biggest wich not crashes exp(x)
878 if(arg<= 0.0) return -709.0;//
880 arg -= TMath::Pi()/2.;
883 return -TMath::Log( TMath::Tan(arg/2.)) ;
887 return TMath::Log( TMath::Tan(-arg/2.)) ;
891 /*******************************************************************/
892 /****** ROUTINES USED FOR COMMUNUCATION ********/
893 /******************** WITH FORTRAN ********************/
894 /*******************************************************************/
897 # define hbtpran hbtpran_
898 # define alihbtp_puttrack alihbtp_puttrack_
899 # define alihbtp_gettrack alihbtp_gettrack_
900 # define alihbtp_getnumberevents alihbtp_getnumberevents_
901 # define alihbtp_getnumbertracks alihbtp_getnumbertracks_
902 # define alihbtp_initialize alihbtp_initialize_
903 # define alihbtp_setactiveeventnumber alihbtp_setactiveeventnumber_
904 # define alihbtp_setparameters alihbtp_setparameters_
908 # define hbtpran HBTPRAN
909 # define alihbtp_puttrack ALIHBTP_PUTTRACK
910 # define alihbtp_gettrack ALIHBTP_GETTRACK
911 # define alihbtp_getnumberevents ALIHBTP_GETNUMBEREVENTS
912 # define alihbtp_getnumbertracks ALIHBTP_GETNUMBERTRACKS
913 # define alihbtp_initialize ALIHBTP_INITIALIZE
914 # define alihbtp_setactiveeventnumber ALIHBTP_SETACTIVEEVENTNUMBER
915 # define alihbtp_setparameters ALIHBTP_SETPARAMETERS
916 # define type_ofCall _stdcall
919 #include "AliGenCocktailAfterBurner.h"
921 /*******************************************************************/
923 AliGenCocktailAfterBurner* GetGenerator()
925 // This function has two tasks:
926 // Check if environment is OK (exist gAlice and generator)
927 // Returns pointer to genrator
928 //to be changed with TFolders
932 cout<<endl<<"ERROR: There is NO gALICE! Check what you are doing!"<<endl;
933 gROOT->Fatal("AliGenHBTprocessor.cxx: GetGenerator()",
934 "\nRunning HBT Processor without gAlice... Exiting \n");
937 AliGenerator * gen = gAlice->Generator();
941 gAlice->Fatal("AliGenHBTprocessor.cxx: GetGenerator()",
942 "\nThere is no generator in gAlice, exiting\n");
946 //we do not sure actual type of the genetator
947 //and simple casting is risky - we use ROOT machinery to do safe cast
949 TClass* cabclass = AliGenCocktailAfterBurner::Class(); //get AliGenCocktailAfterBurner TClass
950 TClass* genclass = gen->IsA();//get TClass of the generator we got from galice
951 //use casting implemented in TClass
952 //cast gen to cabclass
953 AliGenCocktailAfterBurner* cab=(AliGenCocktailAfterBurner*)genclass->DynamicCast(cabclass,gen);
955 if (cab == 0x0)//if generator that we got is not AliGenCocktailAfterBurner or its descendant we get null
956 { //then quit with error
957 gAlice->Fatal("AliGenHBTprocessor.cxx: GetGenerator()",
958 "\nThe main Generator is not a AliGenCocktailAfterBurner, exiting\n");
961 // cout<<endl<<"Got generator"<<endl;
965 /*******************************************************************/
967 AliGenHBTprocessor* GetAliGenHBTprocessor()
969 //returns pointer to the current instance of AliGenHBTprocessor in
970 //AliGenCocktailAfterBurner (can be more than one)
972 AliGenCocktailAfterBurner* gen = GetGenerator();
973 AliGenerator* g = gen->GetCurrentGenerator();
976 gAlice->Fatal("AliGenHBTprocessor.cxx: GetAliGenHBTprocessor()",
977 "Can not get the current generator. Exiting");
981 TClass* hbtpclass = AliGenHBTprocessor::Class(); //get AliGenCocktailAfterBurner TClass
982 TClass* gclass = g->IsA();//get TClass of the current generator we got from CAB
983 AliGenHBTprocessor* hbtp = (AliGenHBTprocessor*)gclass->DynamicCast(hbtpclass,g);//try to cast
986 gAlice->Fatal("AliGenHBTprocessor.cxx: GetAliGenHBTprocessor()",
987 "\nCurrernt generator in AliGenCocktailAfterBurner is not a AliGenHBTprocessor, exiting\n");
993 /*******************************************************************/
994 extern "C" void type_ofCall alihbtp_setparameters()
999 extern "C" void type_ofCall alihbtp_initialize()
1004 /*******************************************************************/
1006 extern "C" void type_ofCall alihbtp_getnumberevents(Int_t &nev)
1008 //passes number of events to the fortran
1009 if(gDebug) cout<<"alihbtp_getnumberevents("<<nev<<") ....";
1010 AliGenCocktailAfterBurner* gen = GetGenerator();
1017 nev = gen->GetNumberOfEvents();
1019 if(gDebug>5) cout<<"EXITED N Ev = "<<nev<<endl;
1023 /*******************************************************************/
1025 extern "C" void type_ofCall alihbtp_setactiveeventnumber(Int_t & nev)
1027 //sets active event in generator (AliGenCocktailAfterBurner)
1029 if(gDebug>5) cout<<"alihbtp_setactiveeventnumber("<<nev<<") ....";
1030 if(gDebug>0) cout<<"Asked for event "<<nev-1<<endl;
1031 AliGenCocktailAfterBurner* gen = GetGenerator();
1033 gen->SetActiveEventNumber(nev - 1); //fortran numerates events from 1 to N
1035 if(gDebug>5) cout<<"EXITED returned "<<nev<<endl;
1038 /*******************************************************************/
1040 extern "C" void type_ofCall alihbtp_getnumbertracks(Int_t &ntracks)
1042 //passes number of particles in active event to the fortran
1043 if(gDebug>5) cout<<"alihbtp_getnumbertracks("<<ntracks<<") ....";
1045 AliGenCocktailAfterBurner* gen = GetGenerator();
1052 ntracks = gen->GetActiveStack()->GetNprimary();
1053 if(gDebug>5) cout<<"EXITED Ntracks = "<<ntracks<<endl;
1056 /*******************************************************************/
1058 extern "C" void type_ofCall
1059 alihbtp_puttrack(Int_t & n,Int_t& flag, Float_t& px,
1060 Float_t& py, Float_t& pz, Int_t& geantpid)
1062 //sets new parameters (momenta) in track number n
1063 // in the active event
1064 // n - number of the track in active event
1065 // flag - flag of the track
1066 // px,py,pz - momenta
1067 // geantpid - type of the particle - Geant Particle ID
1069 if(gDebug>5) cout<<"alihbtp_puttrack("<<n<<") ....";
1071 AliGenCocktailAfterBurner* gen = GetGenerator();
1074 TParticle * track = gen->GetActiveStack()->Particle(n-1);
1076 AliGenHBTprocessor* g = GetAliGenHBTprocessor();
1078 //check to be deleted
1079 if (geantpid != (g->IdFromPDG( track->GetPdgCode() )))
1081 cerr<<endl<<" AliGenHBTprocessor.cxx: alihbtp_puttrack: SOMETHING IS GOING BAD:\n GEANTPIDS ARE NOT THE SAME"<<endl;
1085 if (px != track->Px())
1086 cout<<"Px diff. = "<<px - track->Px()<<endl;
1088 if(gDebug>3) cout<<" track->GetPdgCode() --> "<<track->GetPdgCode()<<endl;
1092 Float_t m = track->GetMass();
1093 Float_t e = TMath::Sqrt(m*m+px*px+py*py+pz*pz);
1094 track->SetMomentum(px,py,pz,e);
1096 g->SetHbtPStatusCode(flag,n-1);
1098 if(gDebug>5) cout<<"EXITED "<<endl;
1101 /*******************************************************************/
1103 extern "C" void type_ofCall
1104 alihbtp_gettrack(Int_t & n,Int_t & flag, Float_t & px,
1105 Float_t & py, Float_t & pz, Int_t & geantpid)
1108 //passes track parameters to the fortran
1109 // n - number of the track in active event
1110 // flag - flag of the track
1111 // px,py,pz - momenta
1112 // geantpid - type of the particle - Geant Particle ID
1114 if(gDebug>5) cout<<"alihbtp_gettrack("<<n<<") ....";
1115 AliGenCocktailAfterBurner* gen = GetGenerator();
1126 TParticle * track = gen->GetActiveStack()->Particle(n-1);
1127 AliGenHBTprocessor* g = GetAliGenHBTprocessor();
1129 flag = g->GetHbtPStatusCode(n-1);
1131 px = (Float_t)track->Px();
1132 py = (Float_t)track->Py();
1133 pz = (Float_t)track->Pz();
1135 geantpid = g->IdFromPDG( track->GetPdgCode() );
1137 if(gDebug>5) cout<<"EXITED "<<endl;
1140 /*******************************************************************/
1141 extern "C" Float_t type_ofCall hbtpran(Int_t &)
1143 //interface to the random number generator
1144 return sRandom->Rndm();
1147 /*******************************************************************/
1150 /*******************************************************************/