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
36 // Artificial particle dennsity enhancment feature
37 // HBT Processor is unable to process correctly low multiplicity particles
38 // even if the high statistics (more events) is supplied.
39 // For that reason was implemented artificial multiplicity enhancement
40 // - see also comments in HBT Processor source file (HBTP/hbt_event_processor.f)
41 // or web page (http://alisoft.cern.ch/people/skowron/hbtprocessor/hbteventprocessor.html)
42 // concerning trk_accep or SetTrackRejectionFactor method in this class
43 // Fortran is cheated by masking several events a single one. Number is defined by fEventMerge
44 // variable. In case it is equal to 1, there is no masking and the feature is off.
45 // But, for example if it is 5, and we have 100 events, number of events passed to fortran is 100/5=20
46 // When fortran asks about multiplcity of event, let sey, 1 - multiplicity of events 5 to 9 is summed up
50 //////////////////////////////////////////////////////////////////////////////////
52 // 11.11.2001 Piotr Skowronski
53 // Setting seed (date) in RNG in the constructor
55 // 09.10.2001 Piotr Skowronski
57 // Theta - Eta cohernecy facilities added:
58 // AliGenerator::SetThetaRange method overriden
59 // Static methods added
65 // Class description comments put on proper place
67 // 27.09.2001 Piotr Skowronski
68 // removing of redefinition of defaults velues
69 // in method's implementation.
73 #include "AliGenHBTprocessor.h"
75 #include <TParticle.h>
76 #include "THBTprocessor.h"
81 #include "AliGenCocktailAfterBurner.h"
86 ClassImp(AliGenHBTprocessor)
88 Int_t AliGenHBTprocessor::fgDebug = 0;
89 static TRandom* gAliRandom;//RNG to be used by the fortran code
91 AliGenCocktailAfterBurner* GetGenerator();
92 /*******************************************************************/
94 AliGenHBTprocessor::AliGenHBTprocessor():
99 fTrackRejectionFactor(0),
100 fReferenceControl(0),
111 fEventLineCounter(0),
165 // Standard constructor
166 // Sets default veues of all parameters
168 SetName("AliGenHBTprocessor");
169 SetTitle("AliGenHBTprocessor");
171 gAliRandom = fRandom;
172 fHBTprocessor = new THBTprocessor();
177 SetTrackRejectionFactor();
197 SetSwitchCoherence();
199 SetSwitchFermiBose();
200 //SetMomentumRange();
213 SetNBins1DFineMesh();
214 SetBinSize1DFineMesh();
215 SetNBins1DCoarseMesh();
216 SetBinSize1DCoarseMesh();
217 SetNBins3DFineMesh();
218 SetBinSize3DFineMesh();
219 SetNBins3DCoarseMesh();
220 SetBinSize3DCoarseMesh();
221 SetNBins3DFineProjectMesh();
224 /*******************************************************************/
227 /*******************************************************************/
229 AliGenHBTprocessor::~AliGenHBTprocessor()
233 if (fHBTprocessor) delete fHBTprocessor; //delete generator
237 /*******************************************************************/
239 void AliGenHBTprocessor::InitStatusCodes()
241 //creates and inits status codes array to zero
242 AliGenCocktailAfterBurner *cab = GetGenerator();
245 Fatal("InitStatusCodes()","Can not find AliGenCocktailAfterBurner generator");
249 Int_t nev = cab->GetNumberOfEvents();
251 fHbtPStatCodes = new Int_t* [nev];
252 for( Int_t i =0; i<nev;i++)
254 Int_t nprim = cab->GetStack(i)->GetNprimary();
255 fHbtPStatCodes[i] = new Int_t[nprim];
256 for (Int_t k =0 ;k<nprim;k++)
257 fHbtPStatCodes[i][k] =0;
262 /*******************************************************************/
264 void AliGenHBTprocessor::CleanStatusCodes()
266 //Cleans up status codes
269 for (Int_t i =0; i<GetGenerator()->GetNumberOfEvents(); i++)
270 delete [] fHbtPStatCodes[i];
271 delete [] fHbtPStatCodes;
276 /*******************************************************************/
278 void AliGenHBTprocessor::Init()
280 //sets up parameters in generator
282 THBTprocessor *thbtp = fHBTprocessor;
285 thbtp->SetTrackRejectionFactor(fTrackRejectionFactor);
286 thbtp->SetRefControl(fReferenceControl);
288 if ((fPid[0] == fPid[1]) || (fPid[0] == 0) || (fPid[1] == 0))
291 thbtp->SetPIDs(IdFromPDG(fPid[1]) ,0);
293 thbtp->SetPIDs(IdFromPDG(fPid[0]) ,0);
294 thbtp->SetNPIDtypes(1);
297 Warning("AliGenHBTprocessor::Init","\nThere is only one particle type set,\n\
298 and Switch_Type differnt then 1,\n which does not make sense.\n\
299 Setting it to 1.\n");
305 thbtp->SetPIDs(IdFromPDG(fPid[0]) ,IdFromPDG(fPid[1]));
307 thbtp->SetSwitchType(fSwitchType);
310 thbtp->SetMaxIterations(fMaxit);
311 thbtp->SetDelChi(fDelchi);
312 thbtp->SetIRand(fIrand);
313 thbtp->SetLambda(fLambda);
315 thbtp->SetRSide(fRside);
316 thbtp->SetROut(fRout);
317 thbtp->SetRLong(fRlong);
318 thbtp->SetRPerp(fRperp);
319 thbtp->SetRParallel(fRparallel);
322 thbtp->SetSwitch1D(fSwitch1d);
323 thbtp->SetSwitch3D(fSwitch3d);
324 thbtp->SetSwitchType(fSwitchType);
325 thbtp->SetSwitchCoherence(fSwitchCoherence);
326 thbtp->SetSwitchCoulomb(fSwitchCoulomb);
327 thbtp->SetSwitchFermiBose(fSwitchFermiBose);
328 thbtp->SetPtRange(fPtMin,fPtMax);
329 thbtp->SetPxRange(fPxMin,fPxMax);
330 thbtp->SetPyRange(fPyMin,fPyMax);
331 thbtp->SetPzRange(fPzMin,fPzMax);
332 thbtp->SetPhiRange(fPhiMin*180./((Float_t)TMath::Pi())+180.0, //casting is because if fPhiMin = 180.0 then
333 fPhiMax*180./((Float_t)TMath::Pi())+180.0);//TMath::Pi() != TMath::Pi()*fPhiMin/180.0,
334 thbtp->SetEtaRange(fEtaMin,fEtaMax);
335 thbtp->SetNPtBins(fNPtBins);
336 thbtp->SetNPhiBins(fNPhiBins);
337 thbtp->SetNEtaBins(fNEtaBins);
338 thbtp->SetNPxBins(fNPxBins);
339 thbtp->SetNPyBins(fNPyBins);
340 thbtp->SetNPzBins(fNPzBins);
341 thbtp->SetNBins1DFineMesh(fN1dFine);
342 thbtp->SetBinSize1DFineMesh(fBinsize1dFine);
343 thbtp->SetNBins1DCoarseMesh(fN1dCoarse);
344 thbtp->SetBinSize1DCoarseMesh(fBinsize1dCoarse);
345 thbtp->SetNBins3DFineMesh(fN3dFine);
346 thbtp->SetBinSize3DFineMesh(fBinsize3dFine);
347 thbtp->SetNBins3DCoarseMesh(fN3dCoarse);
348 thbtp->SetBinSize3DCoarseMesh(fBinsize3dCoarse);
349 thbtp->SetNBins3DFineProjectMesh(fN3dFineProject);
351 thbtp->SetPrintFull(fPrintFull);
354 /*******************************************************************/
356 void AliGenHBTprocessor::Generate()
359 AliGenCocktailAfterBurner* cab = GetGenerator();
362 Fatal("Generate()","AliGenHBTprocessor needs AliGenCocktailAfterBurner to be main generator");
365 if (cab->GetNumberOfEvents() <2)
368 "HBT Processor needs more than 1 event to work properly,\
369 but there is only %d set", cab->GetNumberOfEvents());
373 fHBTprocessor->Initialize(); //reset all fields of common blocks
374 //in case there are many HBT processors
375 //run one after one (i.e. in AliCocktailAfterBurner)
376 //between Init() called and Generate there might
377 Init(); //be different instance running - be sure that we have our settings
379 InitStatusCodes(); //Init status codes
381 fHBTprocessor->GenerateEvent(); //Generates event
383 CleanStatusCodes(); //Clean Status codes - thet are not needed anymore
386 /*******************************************************************/
389 /*******************************************************************/
390 void AliGenHBTprocessor::GetParticles(TClonesArray * particles) const
395 Error("GetParticles","Particles has to be initialized");
398 fHBTprocessor->ImportParticles(particles);
401 /*******************************************************************/
403 Int_t AliGenHBTprocessor::GetHbtPStatusCode(Int_t part) const
405 //returns the status code of the given particle in the active event
406 //see SetActiveEvent in the bottom of AliGenHBTprocessor.cxx
407 //and in AliCocktailAfterBurner
409 GetTrackEventIndex(part,ev,idx);
410 if ( (ev<0) || (idx<0) )
412 Error("GetHbtPStatusCode","GetTrackEventIndex returned error");
415 return fHbtPStatCodes[ev][idx];
419 /*******************************************************************/
420 void AliGenHBTprocessor::SetHbtPStatusCode(Int_t hbtstatcode, Int_t part)
422 //Sets the given status code to given particle number (part) in the active event
424 GetTrackEventIndex(part,ev,idx);
425 if ( (ev<0) || (idx<0) )
427 Error("GetHbtPStatusCode","GetTrackEventIndex returned error");
430 else fHbtPStatCodes[ev][idx] = hbtstatcode;
433 /*******************************************************************/
435 void AliGenHBTprocessor::SetTrackRejectionFactor(Float_t trf) //def 1.0
437 //Sets the Track Rejection Factor
438 //variates in range 0.0 <-> 1.0
439 //Describes the factor of particles rejected from the output.
440 //Used only in case of low muliplicity particles e.g. lambdas.
441 //Processor generates addisional particles and builds the
442 //correletions on such a statistics.
443 //At the end these particels are left in the event according
444 //to this factor: 1==all particles are left
445 // 0==all are removed
447 fTrackRejectionFactor=trf;
448 fHBTprocessor->SetTrackRejectionFactor(trf);
450 /*******************************************************************/
452 void AliGenHBTprocessor::SetRefControl(Int_t rc) //default 2
454 //Sets the Refernce Control Switch
455 //switch wether read reference histograms from file =1
456 // compute from input events =2 - default
457 fReferenceControl=rc;
458 fHBTprocessor->SetRefControl(rc);
460 /*******************************************************************/
462 void AliGenHBTprocessor::SetPIDs(Int_t pid1,Int_t pid2)
465 //Sets PDG Codes of particles to be processed, default \\Pi^{+} and \\Pi^{-}
466 //This method accepts PDG codes which is
467 //in opposite to THBBProcessor which accepts GEANT PID
468 if ( (pid1 == 0) && (pid2 == 0) )
470 Error("AliGenHBTprocessor::SetPIDs","Sensless Particle Types setting: 0 0, Ignoring\n");
478 fHBTprocessor->SetPIDs(IdFromPDG(pid1) ,0);
484 fHBTprocessor->SetPIDs(IdFromPDG(pid1) ,IdFromPDG(pid2));
488 /*******************************************************************/
490 void AliGenHBTprocessor::SetNPIDtypes(Int_t npidt)
492 //Number ofparticle types to be processed - default 2
493 //see AliGenHBTprocessor::SetPIDs
496 fHBTprocessor->SetNPIDtypes(npidt);
498 /*******************************************************************/
500 void AliGenHBTprocessor::SetDeltap(Float_t deltp)
503 //maximum range for random momentum shifts in GeV/c;
504 //px,py,pz independent; Default = 0.1 GeV/c.
506 fHBTprocessor->SetDeltap(deltp);
508 /*******************************************************************/
510 void AliGenHBTprocessor::SetMaxIterations(Int_t maxiter)
512 //maximum # allowed iterations thru all the
513 //tracks for each event. Default = 50.
514 //If maxit=0, then calculate the correlations
515 //for the input set of events without doing the
516 //track adjustment procedure.
519 fHBTprocessor->SetMaxIterations(maxiter);
522 /*******************************************************************/
523 void AliGenHBTprocessor::SetDelChi(Float_t dc)
525 //min % change in total chi-square for which
526 //the track adjustment iterations may stop,
530 fHBTprocessor->SetDelChi(dc);
532 /*******************************************************************/
534 void AliGenHBTprocessor::SetIRand(Int_t irnd)
536 //if fact dummy - used only for compatibility
537 //we are using AliRoot built in RNG
538 //dummy in fact since we are using aliroot build-in RNG
539 //Sets renaodom number generator
541 fHBTprocessor->SetIRand(irnd);
543 /*******************************************************************/
545 void AliGenHBTprocessor::SetLambda(Float_t lam)
548 // Sets Chaoticity Parameter
550 fHBTprocessor->SetLambda(lam);
552 /*******************************************************************/
553 void AliGenHBTprocessor::SetR1d(Float_t r)
556 //Sets Spherical source model radius (fm)
558 fHBTprocessor->SetR1d(r);
560 /*******************************************************************/
561 void AliGenHBTprocessor::SetRSide(Float_t rs)
564 //Rside,Rout,Rlong = Non-spherical Bertsch-Pratt source model (fm)
567 fHBTprocessor->SetRSide(rs);
569 /*******************************************************************/
570 void AliGenHBTprocessor::SetROut(Float_t ro)
573 //Rside,Rout,Rlong = Non-spherical Bertsch-Pratt source model (fm)
575 fHBTprocessor->SetROut(ro);
577 /*******************************************************************/
578 void AliGenHBTprocessor::SetRLong(Float_t rl)
581 //Rside,Rout,Rlong = Non-spherical Bertsch-Pratt source model (fm)
583 fHBTprocessor->SetRLong(rl);
585 /*******************************************************************/
586 void AliGenHBTprocessor::SetRPerp(Float_t rp)
589 //Rperp,Rparallel,R0= Non-spherical Yano-Koonin-Podgoretski source model (fm).
591 fHBTprocessor->SetRPerp(rp);
593 /*******************************************************************/
594 void AliGenHBTprocessor::SetRParallel(Float_t rprl)
597 //Rperp,Rparallel,R0= Non-spherical Yano-Koonin-Podgoretski source model (fm).
599 fHBTprocessor->SetRParallel(rprl);
601 /*******************************************************************/
602 void AliGenHBTprocessor::SetR0(Float_t r0)
605 //Rperp,Rparallel,R0= Non-spherical Yano-Koonin-Podgoretski source model (fm).
607 fHBTprocessor->SetR0(r0);
609 /*******************************************************************/
610 void AliGenHBTprocessor::SetQ0(Float_t q0)
613 //Sets Q0 = NA35 Coulomb parameter for finite source size in (GeV/c)
614 // if fSwitchCoulomb = 2
615 // = Spherical Coulomb source radius in (fm)
616 // if switchCoulomb = 3, used to interpolate the
617 // input Pratt/Cramer discrete Coulomb source
620 fHBTprocessor->SetQ0(q0);
623 /*******************************************************************/
624 void AliGenHBTprocessor::SetSwitch1D(Int_t s1d)
628 // = 0 to not compute the 1D two-body //orrelations.
629 // = 1 to compute this using Q-invariant
630 // = 2 to compute this using Q-total
631 // = 3 to compute this using Q-3-ve//tor
634 fHBTprocessor->SetSwitch1D(s1d);
636 /*******************************************************************/
637 void AliGenHBTprocessor::SetSwitch3D(Int_t s3d)
641 // = 0 to not compute the 3D two-body correlations.
642 // = 1 to compute this using the side-out-long form
643 // = 2 to compute this using the Yanno-Koonin-Pogoredskij form.
646 fHBTprocessor->SetSwitch3D(s3d);
648 /*******************************************************************/
649 void AliGenHBTprocessor::SetSwitchType(Int_t st)
652 // Sets switch_type = 1 to fit only the like pair correlations
653 // = 2 to fit only the unlike pair correlations
654 // = 3 to fit both the like and unlike pair correl.
655 //See SetPIDs and Init
656 //If only one particle type is set, unly==1 makes sens
659 fHBTprocessor->SetSwitchType(st);
661 /*******************************************************************/
662 void AliGenHBTprocessor::SetSwitchCoherence(Int_t sc)
665 // switchCoherence = 0 for purely incoherent source (but can have
667 // = 1 for mixed incoherent and coherent source
669 fSwitchCoherence = sc;
670 fHBTprocessor->SetSwitchCoherence(sc);
672 /*******************************************************************/
673 void AliGenHBTprocessor::SetSwitchCoulomb(Int_t scol)
676 // switchCoulomb = 0 no Coulomb correction
677 // = 1 Point source Gamow correction only
678 // = 2 NA35 finite source size correction
679 // = 3 Pratt/Cramer finite source size correction;
680 // interpolated from input tables.
681 fSwitchCoulomb =scol;
682 fHBTprocessor->SetSwitchCoulomb(scol);
684 /*******************************************************************/
685 void AliGenHBTprocessor::SetSwitchFermiBose(Int_t sfb)
688 // switchFermiBose = 1 Boson pairs
689 // = -1 Fermion pairs
691 fSwitchFermiBose = sfb;
692 fHBTprocessor->SetSwitchFermiBose(sfb);
694 /*******************************************************************/
695 void AliGenHBTprocessor::SetPtRange(Float_t ptmin, Float_t ptmax)
697 // default ptmin = 0.1, ptmax = 0.98
698 //Sets Pt range (GeV)
699 AliGenerator::SetPtRange(ptmin,ptmax);
700 fHBTprocessor->SetPtRange(ptmin,ptmax);
703 /*******************************************************************/
704 void AliGenHBTprocessor::SetPxRange(Float_t pxmin, Float_t pxmax)
706 //default pxmin = -1.0, pxmax = 1.0
710 fHBTprocessor->SetPxRange(pxmin,pxmax);
712 /*******************************************************************/
713 void AliGenHBTprocessor::SetPyRange(Float_t pymin, Float_t pymax)
715 //default pymin = -1.0, pymax = 1.0
719 fHBTprocessor->SetPyRange(pymin,pymax);
721 /*******************************************************************/
722 void AliGenHBTprocessor::SetPzRange(Float_t pzmin, Float_t pzmax)
724 //default pzmin = -3.6, pzmax = 3.6
728 fHBTprocessor->SetPzRange(pzmin,pzmax);
730 void AliGenHBTprocessor::SetMomentumRange(Float_t /*pmin*/, Float_t /*pmax*/)
732 //default pmin=0, pmax=0
733 //Do not use this method!
734 MayNotUse("AliGenHBTprocessor::SetMomentumRange Method is Dummy");
737 /*******************************************************************/
738 void AliGenHBTprocessor::SetPhiRange(Float_t phimin, Float_t phimax)
742 AliGenerator::SetPhiRange(phimin,phimax);
744 fHBTprocessor->SetPhiRange(phimin+180.0,phimax+180.0);
746 /*******************************************************************/
747 void AliGenHBTprocessor::SetEtaRange(Float_t etamin, Float_t etamax)
749 //default etamin = -1.5, etamax = 1.5
753 fHBTprocessor->SetEtaRange(etamin,etamax);
755 //set the azimothal angle range in the AliGeneraor -
756 //to keep coherency between azimuthal angle and pseudorapidity
757 //DO NOT CALL this->SetThetaRange, because it calls this method (where we are)
758 //which must cause INFINITE LOOP
759 AliGenerator::SetThetaRange(RadiansToDegrees(EtaToTheta(fEtaMin)),
760 RadiansToDegrees(EtaToTheta(fEtaMax)));
763 /*******************************************************************/
764 void AliGenHBTprocessor::SetThetaRange(Float_t thetamin, Float_t thetamax)
766 //default thetamin = 0, thetamax = 180
767 //Azimuthal angle, override AliGenerator method which uses widely (i.e. wrapper generators)
768 //core fortran HBTProcessor uses Eta (pseudorapidity)
769 //so these methods has to be synchronized
771 AliGenerator::SetThetaRange(thetamin,thetamax);
773 SetEtaRange( ThetaToEta(fThetaMin) , ThetaToEta(fThetaMax) );
777 /*******************************************************************/
778 void AliGenHBTprocessor::SetNPtBins(Int_t nptbin)
780 //default nptbin = 50
781 //set number of Pt bins
783 fHBTprocessor->SetNPtBins(nptbin);
785 /*******************************************************************/
786 void AliGenHBTprocessor::SetNPhiBins(Int_t nphibin)
788 //default nphibin = 50
789 //set number of Phi bins
791 fHBTprocessor->SetNPhiBins(nphibin);
793 /*******************************************************************/
794 void AliGenHBTprocessor::SetNEtaBins(Int_t netabin)
796 //default netabin = 50
797 //set number of Eta bins
799 fHBTprocessor->SetNEtaBins(netabin);
801 /*******************************************************************/
802 void AliGenHBTprocessor::SetNPxBins(Int_t npxbin)
804 //default npxbin = 20
805 //set number of Px bins
807 fHBTprocessor->SetNPxBins(npxbin);
809 /*******************************************************************/
810 void AliGenHBTprocessor::SetNPyBins(Int_t npybin)
812 //default npybin = 20
813 //set number of Py bins
815 fHBTprocessor->SetNPyBins(npybin);
817 /*******************************************************************/
818 void AliGenHBTprocessor::SetNPzBins(Int_t npzbin)
820 //default npzbin = 70
821 //set number of Pz bins
823 fHBTprocessor->SetNPzBins(npzbin);
825 /*******************************************************************/
826 void AliGenHBTprocessor::SetNBins1DFineMesh(Int_t n)
829 //Sets the number of bins in the 1D mesh
831 fHBTprocessor->SetNBins1DFineMesh(n);
834 /*******************************************************************/
835 void AliGenHBTprocessor::SetBinSize1DFineMesh(Float_t x)
838 //Sets the bin size in the 1D mesh
840 fHBTprocessor->SetBinSize1DFineMesh(x);
842 /*******************************************************************/
844 void AliGenHBTprocessor::SetNBins1DCoarseMesh(Int_t n)
847 //Sets the number of bins in the coarse 1D mesh
849 fHBTprocessor->SetNBins1DCoarseMesh(n);
851 /*******************************************************************/
852 void AliGenHBTprocessor::SetBinSize1DCoarseMesh(Float_t x)
855 //Sets the bin size in the coarse 1D mesh
857 fHBTprocessor->SetBinSize1DCoarseMesh(x);
859 /*******************************************************************/
861 void AliGenHBTprocessor::SetNBins3DFineMesh(Int_t n)
864 //Sets the number of bins in the 3D mesh
866 fHBTprocessor->SetNBins3DFineMesh(n);
868 /*******************************************************************/
869 void AliGenHBTprocessor::SetBinSize3DFineMesh(Float_t x)
872 //Sets the bin size in the 3D mesh
874 fHBTprocessor->SetBinSize3DFineMesh(x);
876 /*******************************************************************/
878 void AliGenHBTprocessor::SetNBins3DCoarseMesh(Int_t n)
881 //Sets the number of bins in the coarse 3D mesh
884 fHBTprocessor->SetNBins3DCoarseMesh(n);
886 /*******************************************************************/
887 void AliGenHBTprocessor::SetBinSize3DCoarseMesh(Float_t x)
890 //Sets the bin size in the coarse 3D mesh
891 fBinsize3dCoarse = x;
892 fHBTprocessor->SetBinSize3DCoarseMesh(x);
894 /*******************************************************************/
896 void AliGenHBTprocessor::SetNBins3DFineProjectMesh(Int_t n )
899 //Sets the number of bins in the fine project mesh
901 fHBTprocessor->SetNBins3DFineProjectMesh(n);
903 /*******************************************************************/
904 void AliGenHBTprocessor::SetPrintFull(Int_t flag)
906 //sets the print full flag
908 fHBTprocessor->SetPrintFull(flag);
912 /*******************************************************************/
914 Int_t AliGenHBTprocessor::GetNumberOfEvents()
916 //returns number of available events
917 AliGenerator* g = gAlice->GetMCApp()->Generator();
918 AliGenCocktailAfterBurner* cab = (g)?dynamic_cast<AliGenCocktailAfterBurner*>(g):0x0;
921 Fatal("GetNumberOfEvents","Master Generator is not an AliGenCocktailAfterBurner");
924 return (Int_t)TMath::Ceil(cab->GetNumberOfEvents()/((Float_t)fEventMerge));
927 /*******************************************************************/
929 void AliGenHBTprocessor::SetActiveEventNumber(Int_t n)
931 //sets the active event
932 fActiveStack = n*fEventMerge;
933 AliDebug(1,Form("Settimg active event %d passed %d",fActiveStack,n));
935 /*******************************************************************/
937 Int_t AliGenHBTprocessor::GetNumberOfTracks()
939 //returns number of tracks in active event
940 AliGenerator* g = gAlice->GetMCApp()->Generator();
941 AliGenCocktailAfterBurner* cab = (g)?dynamic_cast<AliGenCocktailAfterBurner*>(g):0x0;
944 Fatal("GetNumberOfEvents","Master Generator is not an AliGenCocktailAfterBurner");
948 for (Int_t i = fActiveStack;i < fActiveStack+fEventMerge; i++)
950 if (i >= GetNumberOfEvents()) break; //protection not to overshoot nb of events
951 AliStack* stack = cab->GetStack(i);
953 Error("GetNumberOfTracks","There is no stack %d",i);
956 n+=stack->GetNprimary();
960 /*******************************************************************/
962 void AliGenHBTprocessor::SetNEventsToMerge(Int_t nev)
964 //sets number of events to merge
965 if (nev > 0 ) fEventMerge = nev;
967 /*******************************************************************/
969 TParticle* AliGenHBTprocessor::GetTrack(Int_t n)
971 //returns track that hbtp thinks is n in active event
972 AliDebug(5,Form("n = %d",n));
973 AliGenerator* g = gAlice->GetMCApp()->Generator();
974 AliGenCocktailAfterBurner* cab = (g)?dynamic_cast<AliGenCocktailAfterBurner*>(g):0x0;
977 Fatal("GetTrackEventIndex","Master Generator is not an AliGenCocktailAfterBurner");
982 GetTrackEventIndex(n,ev,idx);
983 AliDebug(5,Form("Event = %d Particle = %d",ev,idx));
984 if ( (ev<0) || (idx<0) )
986 Error("GetTrack","GetTrackEventIndex returned error");
989 AliDebug(5,Form("Number of Tracks in Event(%d) = %d",ev,cab->GetStack(ev)->GetNprimary()));
990 return cab->GetStack(ev)->Particle(idx);//safe - in case stack does not exist
992 /*******************************************************************/
994 void AliGenHBTprocessor::GetTrackEventIndex(Int_t n, Int_t &evno, Int_t &index) const
996 //returns event(stack) number and particle index
997 AliGenerator* g = gAlice->GetMCApp()->Generator();
998 AliGenCocktailAfterBurner* cab = (g)?dynamic_cast<AliGenCocktailAfterBurner*>(g):0x0;
1001 Fatal("GetTrackEventIndex","Master Generator is not an AliGenCocktailAfterBurner");
1007 for (Int_t i = fActiveStack;i < fActiveStack+fEventMerge; i++)
1009 AliStack* stack = cab->GetStack(i);
1012 Error("GetTrackEventIndex","There is no stack %d",i);
1016 Int_t ntracks = stack->GetNprimary();
1017 AliDebug(10,Form("Event %d has %d tracks. n = %d",i,ntracks,n));
1031 Error("GetTrackEventIndex","Could not find given track");
1034 /*******************************************************************/
1035 /*******************************************************************/
1036 /*******************************************************************/
1042 void AliGenHBTprocessor::DefineParticles()
1045 // Load standard numbers for GEANT particles and PDG conversion
1046 fNPDGCodes = 0; //this is done in the constructor - but in any case ...
1048 fPDGCode[fNPDGCodes++]=-99; // 0 = unused location
1049 fPDGCode[fNPDGCodes++]=22; // 1 = photon
1050 fPDGCode[fNPDGCodes++]=-11; // 2 = positron
1051 fPDGCode[fNPDGCodes++]=11; // 3 = electron
1052 fPDGCode[fNPDGCodes++]=12; // 4 = neutrino e
1053 fPDGCode[fNPDGCodes++]=-13; // 5 = muon +
1054 fPDGCode[fNPDGCodes++]=13; // 6 = muon -
1055 fPDGCode[fNPDGCodes++]=111; // 7 = pi0
1056 fPDGCode[fNPDGCodes++]=211; // 8 = pi+
1057 fPDGCode[fNPDGCodes++]=-211; // 9 = pi-
1058 fPDGCode[fNPDGCodes++]=130; // 10 = Kaon Long
1059 fPDGCode[fNPDGCodes++]=321; // 11 = Kaon +
1060 fPDGCode[fNPDGCodes++]=-321; // 12 = Kaon -
1061 fPDGCode[fNPDGCodes++]=2112; // 13 = Neutron
1062 fPDGCode[fNPDGCodes++]=2212; // 14 = Proton
1063 fPDGCode[fNPDGCodes++]=-2212; // 15 = Anti Proton
1064 fPDGCode[fNPDGCodes++]=310; // 16 = Kaon Short
1065 fPDGCode[fNPDGCodes++]=221; // 17 = Eta
1066 fPDGCode[fNPDGCodes++]=3122; // 18 = Lambda
1067 fPDGCode[fNPDGCodes++]=3222; // 19 = Sigma +
1068 fPDGCode[fNPDGCodes++]=3212; // 20 = Sigma 0
1069 fPDGCode[fNPDGCodes++]=3112; // 21 = Sigma -
1070 fPDGCode[fNPDGCodes++]=3322; // 22 = Xi0
1071 fPDGCode[fNPDGCodes++]=3312; // 23 = Xi-
1072 fPDGCode[fNPDGCodes++]=3334; // 24 = Omega-
1073 fPDGCode[fNPDGCodes++]=-2112; // 25 = Anti Proton
1074 fPDGCode[fNPDGCodes++]=-3122; // 26 = Anti Proton
1075 fPDGCode[fNPDGCodes++]=-3222; // 27 = Anti Sigma -
1076 fPDGCode[fNPDGCodes++]=-3212; // 28 = Anti Sigma 0
1077 fPDGCode[fNPDGCodes++]=-3112; // 29 = Anti Sigma 0
1078 fPDGCode[fNPDGCodes++]=-3322; // 30 = Anti Xi 0
1079 fPDGCode[fNPDGCodes++]=-3312; // 31 = Anti Xi +
1080 fPDGCode[fNPDGCodes++]=-3334; // 32 = Anti Omega +
1083 /*******************************************************************/
1084 Int_t AliGenHBTprocessor::IdFromPDG(Int_t pdg) const
1087 // Return Geant3 code from PDG and pseudo ENDF code
1089 for(Int_t i=0;i<fNPDGCodes;++i)
1090 if(pdg==fPDGCode[i]) return i;
1093 Int_t AliGenHBTprocessor::PDGFromId(Int_t id) const
1096 // Return PDG code and pseudo ENDF code from Geant3 code
1098 if(id>0 && id<fNPDGCodes) return fPDGCode[id];
1101 Double_t AliGenHBTprocessor::ThetaToEta(Double_t arg)
1103 //converts etha(azimuthal angle) to Eta (pseudorapidity). Argument in radians
1105 if(arg>= TMath::Pi()) return 708.0;//This number is the biggest wich not crashes exp(x)
1106 if(arg<= 0.0) return -708.0;//
1107 if(arg == TMath::Pi()/2.) return 0.0;//
1111 return TMath::Log( TMath::Tan(arg/2.)) ;
1115 return -TMath::Log( TMath::Tan(-arg/2.)) ;
1119 /*******************************************************************/
1120 /****** ROUTINES USED FOR COMMUNUCATION ********/
1121 /******************** WITH FORTRAN ********************/
1122 /*******************************************************************/
1125 # define hbtpran hbtpran_
1126 # define alihbtp_puttrack alihbtp_puttrack_
1127 # define alihbtp_gettrack alihbtp_gettrack_
1128 # define alihbtp_getnumberevents alihbtp_getnumberevents_
1129 # define alihbtp_getnumbertracks alihbtp_getnumbertracks_
1130 # define alihbtp_initialize alihbtp_initialize_
1131 # define alihbtp_setactiveeventnumber alihbtp_setactiveeventnumber_
1132 # define alihbtp_setparameters alihbtp_setparameters_
1133 # define type_ofCall
1136 # define hbtpran HBTPRAN
1137 # define alihbtp_puttrack ALIHBTP_PUTTRACK
1138 # define alihbtp_gettrack ALIHBTP_GETTRACK
1139 # define alihbtp_getnumberevents ALIHBTP_GETNUMBEREVENTS
1140 # define alihbtp_getnumbertracks ALIHBTP_GETNUMBERTRACKS
1141 # define alihbtp_initialize ALIHBTP_INITIALIZE
1142 # define alihbtp_setactiveeventnumber ALIHBTP_SETACTIVEEVENTNUMBER
1143 # define alihbtp_setparameters ALIHBTP_SETPARAMETERS
1144 # define type_ofCall _stdcall
1147 /*******************************************************************/
1149 AliGenCocktailAfterBurner* GetGenerator()
1151 // This function has two tasks:
1152 // Check if environment is OK (exist gAlice and generator)
1153 // Returns pointer to genrator
1154 //to be changed with TFolders
1158 ::Error("AliGenHBTprocessor.cxx: GetGenerator()",
1159 "There is NO gALICE! Check what you are doing!");
1160 ::Fatal("AliGenHBTprocessor.cxx: GetGenerator()",
1161 "Running HBT Processor without gAlice... Exiting \n");
1162 return 0x0;//pro forma
1164 AliGenerator * gen = gAlice->GetMCApp()->Generator();
1168 ::Fatal("AliGenHBTprocessor.cxx: GetGenerator()",
1169 "There is no generator in gAlice, exiting\n");
1173 //we do not sure actual type of the genetator
1174 //and simple casting is risky - we use ROOT machinery to do safe cast
1176 TClass* cabclass = AliGenCocktailAfterBurner::Class(); //get AliGenCocktailAfterBurner TClass
1177 TClass* genclass = gen->IsA();//get TClass of the generator we got from galice
1178 //use casting implemented in TClass
1179 //cast gen to cabclass
1180 AliGenCocktailAfterBurner* cab=(AliGenCocktailAfterBurner*)genclass->DynamicCast(cabclass,gen);
1182 if (cab == 0x0)//if generator that we got is not AliGenCocktailAfterBurner or its descendant we get null
1183 { //then quit with error
1184 ::Fatal("AliGenHBTprocessor.cxx: GetGenerator()",
1185 "The main Generator is not a AliGenCocktailAfterBurner, exiting\n");
1190 /*******************************************************************/
1192 AliGenHBTprocessor* GetAliGenHBTprocessor()
1194 //returns pointer to the current instance of AliGenHBTprocessor in
1195 //AliGenCocktailAfterBurner (can be more than one)
1197 AliGenCocktailAfterBurner* gen = GetGenerator();
1198 AliGenerator* g = gen->GetCurrentGenerator();
1201 ::Fatal("AliGenHBTprocessor.cxx: GetAliGenHBTprocessor()",
1202 "Can not get the current generator. Exiting");
1203 return 0x0;//pro forma
1206 TClass* hbtpclass = AliGenHBTprocessor::Class(); //get AliGenCocktailAfterBurner TClass
1207 TClass* gclass = g->IsA();//get TClass of the current generator we got from CAB
1208 AliGenHBTprocessor* hbtp = (AliGenHBTprocessor*)gclass->DynamicCast(hbtpclass,g);//try to cast
1211 ::Fatal("AliGenHBTprocessor.cxx: GetAliGenHBTprocessor()",
1212 "\nCurrernt generator in AliGenCocktailAfterBurner is not a AliGenHBTprocessor, exiting\n");
1218 /*******************************************************************/
1219 extern "C" void type_ofCall alihbtp_setparameters()
1224 extern "C" void type_ofCall alihbtp_initialize()
1229 /*******************************************************************/
1231 extern "C" void type_ofCall alihbtp_getnumberevents(Int_t &nev)
1233 //passes number of events to the fortran
1234 if(AliGenHBTprocessor::GetDebug())
1235 AliInfoGeneral("alihbtp_getnumberevents",Form("(%d) ....",nev));
1237 AliGenHBTprocessor* gen = GetAliGenHBTprocessor();//we dont check because it is done in function
1238 nev = gen->GetNumberOfEvents();
1240 if(AliGenHBTprocessor::GetDebug()>5)
1241 AliInfoGeneral("alihbtp_getnumberevents",Form("EXITED N Ev = %d",nev));
1244 /*******************************************************************/
1246 extern "C" void type_ofCall alihbtp_setactiveeventnumber(Int_t & nev)
1248 //sets active event in generator (AliGenCocktailAfterBurner)
1250 if(AliGenHBTprocessor::GetDebug()>5)
1251 ::Info("AliGenHBTpocessor.cxx: alihbtp_setactiveeventnumber","(%d)",nev);
1252 if(AliGenHBTprocessor::GetDebug()>0)
1253 ::Info("AliGenHBTpocessor.cxx: alihbtp_setactiveeventnumber","Asked for event %d",nev-1);
1255 AliGenHBTprocessor* gen = GetAliGenHBTprocessor();//we dont check because it is done in function
1257 gen->SetActiveEventNumber(nev - 1); //fortran numerates events from 1 to N
1259 if(AliGenHBTprocessor::GetDebug()>5)
1260 ::Info("AliGenHBTpocessor.cxx: alihbtp_setactiveeventnumber","EXITED returned %d",nev);
1262 /*******************************************************************/
1264 extern "C" void type_ofCall alihbtp_getnumbertracks(Int_t &ntracks)
1266 //passes number of particles in active event to the fortran
1267 if(AliGenHBTprocessor::GetDebug()>5)
1268 ::Info("AliGenHBTpocessor.cxx: alihbtp_getnumbertracks","(%d)",ntracks);
1270 AliGenHBTprocessor* gen = GetAliGenHBTprocessor();//we dont check because it is done in function
1272 ntracks = gen->GetNumberOfTracks();
1273 if(AliGenHBTprocessor::GetDebug()>5)
1274 ::Info("AliGenHBTpocessor.cxx: alihbtp_getnumbertracks","EXITED Ntracks = %d",ntracks);
1277 /*******************************************************************/
1279 extern "C" void type_ofCall
1280 alihbtp_puttrack(Int_t & n,Int_t& flag, Float_t& px,
1281 Float_t& py, Float_t& pz, Int_t& geantpid)
1283 //sets new parameters (momenta) in track number n
1284 // in the active event
1285 // n - number of the track in active event
1286 // flag - flag of the track
1287 // px,py,pz - momenta
1288 // geantpid - type of the particle - Geant Particle ID
1290 if(AliGenHBTprocessor::GetDebug()>5)
1291 ::Info("AliGenHBTpocessor.cxx: alihbtp_puttrack","(%d)",n);
1293 AliGenHBTprocessor* gen = GetAliGenHBTprocessor();//we dont check because it is done in function
1295 TParticle * track = gen->GetTrack(n-1);
1298 ::Error("AliGenHBTprocessor.cxx","Can not get track from AliGenHBTprocessor");
1302 //check to be deleted
1303 if (geantpid != (gen->IdFromPDG( track->GetPdgCode() )))
1305 ::Error("AliGenHBTprocessor.cxx: alihbtp_puttrack",
1306 "SOMETHING IS GOING BAD:\n GEANTPIDS ARE NOT THE SAME");
1309 if(AliGenHBTprocessor::GetDebug()>0)
1310 if (px != track->Px())
1311 ::Info("AliGenHBTprocessor.cxx: alihbtp_puttrack",
1312 "Px diff. = %f", px - track->Px());
1314 if(AliGenHBTprocessor::GetDebug()>3)
1315 ::Info("AliGenHBTprocessor.cxx: alihbtp_puttrack",
1316 "track->GetPdgCode() --> %d",track->GetPdgCode());
1320 Float_t m = track->GetMass();
1321 Float_t e = TMath::Sqrt(m*m+px*px+py*py+pz*pz);
1322 track->SetMomentum(px,py,pz,e);
1324 gen->SetHbtPStatusCode(flag,n-1);
1326 if(AliGenHBTprocessor::GetDebug()>5) ::Info("AliGenHBTprocessor.cxx: alihbtp_puttrack","EXITED ");
1329 /*******************************************************************/
1331 extern "C" void type_ofCall
1332 alihbtp_gettrack(Int_t & n,Int_t & flag, Float_t & px,
1333 Float_t & py, Float_t & pz, Int_t & geantpid)
1336 //passes track parameters to the fortran
1337 // n - number of the track in active event
1338 // flag - flag of the track
1339 // px,py,pz - momenta
1340 // geantpid - type of the particle - Geant Particle ID
1342 if(AliGenHBTprocessor::GetDebug()>5) ::Info("AliGenHBTprocessor.cxx: alihbtp_gettrack","(%d)",n);
1343 AliGenHBTprocessor* g = GetAliGenHBTprocessor();
1345 TParticle * track = g->GetTrack(n-1);
1347 flag = g->GetHbtPStatusCode(n-1);
1349 px = (Float_t)track->Px();
1350 py = (Float_t)track->Py();
1351 pz = (Float_t)track->Pz();
1353 geantpid = g->IdFromPDG( track->GetPdgCode() );
1355 if(AliGenHBTprocessor::GetDebug()>5) ::Info("AliGenHBTprocessor.cxx: alihbtp_gettrack","EXITED");
1358 /*******************************************************************/
1359 extern "C" Float_t type_ofCall hbtpran(Int_t &)
1361 //interface to the random number generator
1362 return gAliRandom->Rndm();
1365 /*******************************************************************/
1368 /*******************************************************************/