pamVMC/include/PamVMCPrimaryGenerator.h

// $Id: PamVMCPrimaryGenerator.h,v 1.0 2007/06/03 

#ifndef PAMVMC_PRIMARY_GENERATOR_H
#define PAMVMC_PRIMARY_GENERATOR_H

#include <iostream>
#include <TVirtualMCApplication.h>
#include <TClonesArray.h>
#include <TVector3.h>
#include <TVector2.h>
#include <TMath.h>
#include <TRandom3.h>

#include "PamRootManager.h"

class TVirtualMCStack;

using std::cout;
using std::endl;
using TMath::Sqrt;

class PamVMCPrimary : public TObject
{

 public:
  
  PamVMCPrimary(): fPDG(0), fX0(0.), fY0(0.), fZ0(0.), fTHETA(0.), fPHI(0.), fP0(0.), fGOOD(kFALSE) { };


  Int_t fPDG;
  Double_t fX0, fY0, fZ0;
  Double_t fTHETA, fPHI;  // RADIANTS
  Double_t fP0;
  Bool_t fGOOD;

  void Clean() { fPDG=0; fX0=fY0=fZ0=fTHETA=fPHI=fP0=0.; fGOOD=kFALSE;}

  void Clear(const Option_t * =""){ Clean(); }

  void Print(const Option_t * ="") const
    { cout<<"PRIMARY particle infromation:"<<endl;
      cout<<"Pdg="<<fPDG<<endl; //add name later
      cout<<"Position: "<<"("<<fX0<<","<<fY0<<","<<fZ0<<")"<<endl;
      cout<<"P0, Theta, Phi: "<<fP0<<","<<fTHETA*180./3.1415<<","<<fPHI*180./3.1415<<endl;
      cout<<"GOOD Single Track:"<<fGOOD<<endl;
    }

  ClassDef(PamVMCPrimary,1);
};

PamVMCPrimary & operator+=(PamVMCPrimary &a, const PamVMCPrimary &b);

class PamVMCPrimaryGenerator : public TObject
{
  public:
    PamVMCPrimaryGenerator(TVirtualMCStack* stack, UInt_t seed); 
    //PamVMCPrimaryGenerator(UInt_t seed);
    PamVMCPrimaryGenerator();

    virtual ~PamVMCPrimaryGenerator();

    // methods
    
    void GeneratePrimary();
    // set methods
    void SetParticle(Int_t pdg);
    void SetPosition(Double_t vx, Double_t vy, Double_t vz) 
         {fprim.fX0=vx; fprim.fY0=vy; fprim.fZ0=vz;};
   
    void SetKinEnergy(Double_t kinEnergy) {fprim.fP0=KinEToMomentum(kinEnergy);};
    void SetRigidity(Double_t rigidity) {fprim.fP0=RigToMomentum(rigidity);};
    void SetMomentum(Double_t momentum) {fprim.fP0=momentum; };

    void SetDirection(Double_t theta, Double_t phi) {fprim.fTHETA=theta; fprim.fPHI=phi; };
    void SetMomentum(Double_t px, Double_t py, Double_t pz);

    void SetGood(Bool_t isgood) {fprim.fGOOD=isgood; };

    // gen methods

    void GenPosition(Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax, Double_t zmin, Double_t zmax){
      SetPosition(frnd->Uniform(xmin,xmax),frnd->Uniform(ymin,ymax),frnd->Uniform(zmin,zmax)); }
        //Spherical distribution in radius, thetamin, thetamax, phimin, phimax.
    void GenSphDist(Double_t, Double_t, Double_t, Double_t, Double_t);
    
    void GenDirection(Double_t thetamin, Double_t thetamax, Double_t phimin, Double_t phimax){
        SetDirection(frnd->Uniform(thetamin,thetamax),frnd->Uniform(phimin,phimax)); }

    //flat spectra generator
    void GenSpe(Double_t PEmin, Double_t PEmax, Bool_t isEnergy=kFALSE); 
    //power law spectra, gamma - differential spectral index
    void GenSpe(Double_t PEmin, Double_t PEmax, Double_t gamma, Bool_t isEnergy=kFALSE); 

    //spherical dist. from radius r, allow only those particles going to cross 3 TOFs
    void GenSphDist(Double_t);


    void GenSphPhiThe(Double_t , Double_t , Double_t , Double_t , 
                      Double_t , Double_t ); // cecilia

    // get methods
    Int_t GetParticle(){ return fprim.fPDG; };
    void  GetPositon(TVector3 & v){ v.SetXYZ(fprim.fX0,fprim.fY0,fprim.fZ0); };
    void  GetPosition(Double_t X0, Double_t Y0,Double_t Z0)
          {X0=fprim.fX0; Y0=fprim.fY0; Z0=fprim.fZ0;};

    Double_t GetKinEnergy() { return MomentumToKinE(fprim.fP0); };
    Double_t GetRigidity() { return MomentumToRig(fprim.fP0); };
    Double_t GetMomentum() { return fprim.fP0; };
    
    void GetDirection(TVector2 & v) {v.Set(fprim.fTHETA,fprim.fPHI);}
    
    Bool_t Getgood(){ return fprim.fGOOD; };

    //work with collection of primaries

    void Register(){
      PamRootManager::Instance()->
        Register("PRIM","TClonesArray", &fprimColl);
    }

    void ClearPrimCol(){
      fevno=0;
      (*fprimColl).Clear();
    }

        void Print(const Option_t * ="") const {fprim.Print();};

  private:
    // methods
    
    Double_t MomentumToKinE(Double_t P0);
    Double_t KinEToMomentum(Double_t E0);
    Double_t MomentumToRig(Double_t P0) { return P0/fcharge; };
    Double_t RigToMomentum(Double_t R0){ return R0*fcharge; };

    // data members
    TVirtualMCStack*  fStack;
    Int_t             fevno;
    PamVMCPrimary     fprim;
    Double_t          fmass;
    Double_t          fcharge;
    TClonesArray*     fprimColl; 

    TRandom3 *        frnd;
    
  ClassDef(PamVMCPrimaryGenerator,2)  //PamVMCPrimaryGenerator
};

// inline functions


inline Double_t PamVMCPrimaryGenerator::MomentumToKinE(Double_t P0)
{ return  Sqrt(P0*P0+fmass*fmass)-fmass; }

inline Double_t PamVMCPrimaryGenerator::KinEToMomentum(Double_t E0)
{ return  Sqrt(E0*E0+2.*E0*fmass); }

#endif //PAMVMC_PRIMARY_GENERATOR_H

1	pamelats	1.1	// $Id: PamVMCPrimaryGenerator.h,v 1.0 2007/06/03
2
3			#ifndef PAMVMC_PRIMARY_GENERATOR_H
4			#define PAMVMC_PRIMARY_GENERATOR_H
5
6			#include <iostream>
7			#include <TVirtualMCApplication.h>
8			#include <TClonesArray.h>
9			#include <TVector3.h>
10			#include <TVector2.h>
11			#include <TMath.h>
12			#include <TRandom3.h>
13
14			#include "PamRootManager.h"
15
16			class TVirtualMCStack;
17
18			using std::cout;
19			using std::endl;
20			using TMath::Sqrt;
21
22			class PamVMCPrimary : public TObject
23			{
24
25			public:
26
27			PamVMCPrimary(): fPDG(0), fX0(0.), fY0(0.), fZ0(0.), fTHETA(0.), fPHI(0.), fP0(0.), fGOOD(kFALSE) { };
28
29
30			Int_t fPDG;
31			Double_t fX0, fY0, fZ0;
32	pizzolot	1.2	Double_t fTHETA, fPHI; // RADIANTS
33	pamelats	1.1	Double_t fP0;
34			Bool_t fGOOD;
35
36			void Clean() { fPDG=0; fX0=fY0=fZ0=fTHETA=fPHI=fP0=0.; fGOOD=kFALSE;}
37
38			void Clear(const Option_t * =""){ Clean(); }
39
40			void Print(const Option_t * ="") const
41			{ cout<<"PRIMARY particle infromation:"<<endl;
42			cout<<"Pdg="<<fPDG<<endl; //add name later
43			cout<<"Position: "<<"("<<fX0<<","<<fY0<<","<<fZ0<<")"<<endl;
44	pizzolot	1.2	cout<<"P0, Theta, Phi: "<<fP0<<","<<fTHETA180./3.1415<<","<<fPHI180./3.1415<<endl;
45	pamelats	1.1	cout<<"GOOD Single Track:"<<fGOOD<<endl;
46			}
47
48			ClassDef(PamVMCPrimary,1);
49			};
50
51			PamVMCPrimary & operator+=(PamVMCPrimary &a, const PamVMCPrimary &b);
52
53			class PamVMCPrimaryGenerator : public TObject
54			{
55			public:
56	pizzolot	1.2	PamVMCPrimaryGenerator(TVirtualMCStack* stack, UInt_t seed);
57			//PamVMCPrimaryGenerator(UInt_t seed);
58	pamelats	1.1	PamVMCPrimaryGenerator();
59
60			virtual ~PamVMCPrimaryGenerator();
61
62			// methods
63
64			void GeneratePrimary();
65			// set methods
66			void SetParticle(Int_t pdg);
67			void SetPosition(Double_t vx, Double_t vy, Double_t vz)
68			{fprim.fX0=vx; fprim.fY0=vy; fprim.fZ0=vz;};
69
70			void SetKinEnergy(Double_t kinEnergy) {fprim.fP0=KinEToMomentum(kinEnergy);};
71			void SetRigidity(Double_t rigidity) {fprim.fP0=RigToMomentum(rigidity);};
72			void SetMomentum(Double_t momentum) {fprim.fP0=momentum; };
73
74			void SetDirection(Double_t theta, Double_t phi) {fprim.fTHETA=theta; fprim.fPHI=phi; };
75			void SetMomentum(Double_t px, Double_t py, Double_t pz);
76
77			void SetGood(Bool_t isgood) {fprim.fGOOD=isgood; };
78
79			// gen methods
80
81			void GenPosition(Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax, Double_t zmin, Double_t zmax){
82			SetPosition(frnd->Uniform(xmin,xmax),frnd->Uniform(ymin,ymax),frnd->Uniform(zmin,zmax)); }
83			//Spherical distribution in radius, thetamin, thetamax, phimin, phimax.
84			void GenSphDist(Double_t, Double_t, Double_t, Double_t, Double_t);
85
86			void GenDirection(Double_t thetamin, Double_t thetamax, Double_t phimin, Double_t phimax){
87	pizzolot	1.2	SetDirection(frnd->Uniform(thetamin,thetamax),frnd->Uniform(phimin,phimax)); }
88	pamelats	1.1
89			//flat spectra generator
90			void GenSpe(Double_t PEmin, Double_t PEmax, Bool_t isEnergy=kFALSE);
91			//power law spectra, gamma - differential spectral index
92			void GenSpe(Double_t PEmin, Double_t PEmax, Double_t gamma, Bool_t isEnergy=kFALSE);
93
94	pizzolot	1.2	//spherical dist. from radius r, allow only those particles going to cross 3 TOFs
95			void GenSphDist(Double_t);
96
97
98			void GenSphPhiThe(Double_t , Double_t , Double_t , Double_t ,
99			Double_t , Double_t ); // cecilia
100	pamelats	1.1
101			// get methods
102			Int_t GetParticle(){ return fprim.fPDG; };
103			void GetPositon(TVector3 & v){ v.SetXYZ(fprim.fX0,fprim.fY0,fprim.fZ0); };
104			void GetPosition(Double_t X0, Double_t Y0,Double_t Z0)
105			{X0=fprim.fX0; Y0=fprim.fY0; Z0=fprim.fZ0;};
106
107			Double_t GetKinEnergy() { return MomentumToKinE(fprim.fP0); };
108			Double_t GetRigidity() { return MomentumToRig(fprim.fP0); };
109			Double_t GetMomentum() { return fprim.fP0; };
110
111			void GetDirection(TVector2 & v) {v.Set(fprim.fTHETA,fprim.fPHI);}
112
113			Bool_t Getgood(){ return fprim.fGOOD; };
114
115			//work with collection of primaries
116
117			void Register(){
118			PamRootManager::Instance()->
119			Register("PRIM","TClonesArray", &fprimColl);
120			}
121
122			void ClearPrimCol(){
123			fevno=0;
124			(*fprimColl).Clear();
125			}
126
127			void Print(const Option_t * ="") const {fprim.Print();};
128
129			private:
130			// methods
131
132			Double_t MomentumToKinE(Double_t P0);
133			Double_t KinEToMomentum(Double_t E0);
134			Double_t MomentumToRig(Double_t P0) { return P0/fcharge; };
135			Double_t RigToMomentum(Double_t R0){ return R0*fcharge; };
136
137			// data members
138			TVirtualMCStack* fStack;
139			Int_t fevno;
140			PamVMCPrimary fprim;
141			Double_t fmass;
142			Double_t fcharge;
143			TClonesArray* fprimColl;
144
145			TRandom3 * frnd;
146
147	pizzolot	1.2	ClassDef(PamVMCPrimaryGenerator,2) //PamVMCPrimaryGenerator
148	pamelats	1.1	};
149
150			// inline functions
151
152
153			inline Double_t PamVMCPrimaryGenerator::MomentumToKinE(Double_t P0)
154			{ return Sqrt(P0P0+fmassfmass)-fmass; }
155
156			inline Double_t PamVMCPrimaryGenerator::KinEToMomentum(Double_t E0)
157			{ return Sqrt(E0E0+2.E0*fmass); }
158
159			#endif //PAMVMC_PRIMARY_GENERATOR_H
160