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 <TRandom.h>
#include <TF1.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;
  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<<","<<fPHI<<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); 
    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(frandom->Uniform(xmin,xmax),frandom->Uniform(ymin,ymax),frandom->Uniform(zmin,zmax)); }

    void GenDirection(Double_t thetamin, Double_t thetamax, Double_t phimin, Double_t phimax){
      SetDirection(ftheta->GetRandom(thetamin, thetamax), frandom->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); 
    void GenSpe_Flat(Double_t PEmin, Double_t PEmax, Double_t gamma, Bool_t isEnergy=kFALSE); 
    void GenSpe_3par(Double_t PEmin, Double_t PEmax, Double_t a, Double_t b, Double_t c);
    void GenSphericalPhiThe(); // cecilia // flusso istropo // sets phi the
    void GenSphPhiThe(Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax, 
                      Double_t zmin, Double_t zmax); // flusso isotropo entro S2 S3
    //                                                  sets position and phi the
    

    // 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; };

    //initialize random
    void SetRandom(TRandom* random){
      frandom = random;
    }

    //work with collection of primaries

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

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

  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; };
    Double_t function3par(Double_t xx, Double_t a, Double_t b, Double_t c){return  5.*pow((xx + b * exp(-c * sqrt(xx))),-a);};

    // data members
    TVirtualMCStack*  fStack;
    Int_t             fevno;
    PamVMCPrimary     fprim;
    Double_t          fmass;
    Double_t          fcharge;
    TClonesArray*     fprimColl; 
    TRandom*          frandom; // Class is not a owner of this object
    TF1*              ftheta; // To generate sherical distributhin in theta-angle
    
  ClassDef(PamVMCPrimaryGenerator,1)  //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	// $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 <TRandom.h>
13	#include <TF1.h>
14
15	#include "PamRootManager.h"
16
17	class TVirtualMCStack;
18
19	using std::cout;
20	using std::endl;
21	using TMath::Sqrt;
22
23	class PamVMCPrimary : public TObject
24	{
25
26	public:
27
28	PamVMCPrimary(): fPDG(0), fX0(0.), fY0(0.), fZ0(0.), fTHETA(0.), fPHI(0.), fP0(0.), fGOOD(kFALSE) { };
29
30
31	Int_t fPDG;
32	Double_t fX0, fY0, fZ0;
33	Double_t fTHETA, fPHI;
34	Double_t fP0;
35	Bool_t fGOOD;
36
37	void Clean() { fPDG=0; fX0=fY0=fZ0=fTHETA=fPHI=fP0=0.; fGOOD=kFALSE;}
38
39	void Clear(const Option_t * =""){ Clean(); }
40
41	void Print(const Option_t * ="") const
42	{ cout<<"PRIMARY particle infromation:"<<endl;
43	cout<<"Pdg="<<fPDG<<endl; //add name later
44	cout<<"Position: "<<"("<<fX0<<","<<fY0<<","<<fZ0<<")"<<endl;
45	cout<<"P0, Theta, Phi: "<<fP0<<","<<fTHETA<<","<<fPHI<<endl;
46	cout<<"GOOD Single Track:"<<fGOOD<<endl;
47	}
48
49	ClassDef(PamVMCPrimary,1);
50	};
51
52	PamVMCPrimary & operator+=(PamVMCPrimary &a, const PamVMCPrimary &b);
53
54	class PamVMCPrimaryGenerator : public TObject
55	{
56	public:
57	PamVMCPrimaryGenerator(TVirtualMCStack* stack);
58	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(frandom->Uniform(xmin,xmax),frandom->Uniform(ymin,ymax),frandom->Uniform(zmin,zmax)); }
83
84	void GenDirection(Double_t thetamin, Double_t thetamax, Double_t phimin, Double_t phimax){
85	SetDirection(ftheta->GetRandom(thetamin, thetamax), frandom->Uniform(phimin,phimax));
86	}
87	//flat spectra generator
88	void GenSpe(Double_t PEmin, Double_t PEmax, Bool_t isEnergy=kFALSE);
89	//power law spectra, gamma - differential spectral index
90	void GenSpe(Double_t PEmin, Double_t PEmax, Double_t gamma, Bool_t isEnergy=kFALSE);
91	void GenSpe_Flat(Double_t PEmin, Double_t PEmax, Double_t gamma, Bool_t isEnergy=kFALSE);
92	void GenSpe_3par(Double_t PEmin, Double_t PEmax, Double_t a, Double_t b, Double_t c);
93	void GenSphericalPhiThe(); // cecilia // flusso istropo // sets phi the
94	void GenSphPhiThe(Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax,
95	Double_t zmin, Double_t zmax); // flusso isotropo entro S2 S3
96	// sets position and phi the
97
98
99	// get methods
100	Int_t GetParticle(){ return fprim.fPDG; };
101	void GetPositon(TVector3 & v){ v.SetXYZ(fprim.fX0,fprim.fY0,fprim.fZ0); };
102	void GetPosition(Double_t X0, Double_t Y0,Double_t Z0)
103	{X0=fprim.fX0; Y0=fprim.fY0; Z0=fprim.fZ0;};
104
105	Double_t GetKinEnergy() { return MomentumToKinE(fprim.fP0); };
106	Double_t GetRigidity() { return MomentumToRig(fprim.fP0); };
107	Double_t GetMomentum() { return fprim.fP0; };
108
109	void GetDirection(TVector2 & v) {v.Set(fprim.fTHETA,fprim.fPHI);}
110
111	Bool_t Getgood(){ return fprim.fGOOD; };
112
113	//initialize random
114	void SetRandom(TRandom* random){
115	frandom = random;
116	}
117
118	//work with collection of primaries
119
120	void Register(){
121	PamRootManager::Instance()->
122	Register("PRIM","TClonesArray", &fprimColl);
123	}
124
125	void ClearPrimCol(){
126	fevno=0;
127	(*fprimColl).Clear();
128	}
129
130	private:
131	// methods
132
133	Double_t MomentumToKinE(Double_t P0);
134	Double_t KinEToMomentum(Double_t E0);
135	Double_t MomentumToRig(Double_t P0) { return P0/fcharge; };
136	Double_t RigToMomentum(Double_t R0){ return R0*fcharge; };
137	Double_t function3par(Double_t xx, Double_t a, Double_t b, Double_t c){return 5.pow((xx + b exp(-c * sqrt(xx))),-a);};
138
139	// data members
140	TVirtualMCStack* fStack;
141	Int_t fevno;
142	PamVMCPrimary fprim;
143	Double_t fmass;
144	Double_t fcharge;
145	TClonesArray* fprimColl;
146	TRandom* frandom; // Class is not a owner of this object
147	TF1* ftheta; // To generate sherical distributhin in theta-angle
148
149	ClassDef(PamVMCPrimaryGenerator,1) //PamVMCPrimaryGenerator
150	};
151
152	// inline functions
153
154
155	inline Double_t PamVMCPrimaryGenerator::MomentumToKinE(Double_t P0)
156	{ return Sqrt(P0P0+fmassfmass)-fmass; }
157
158	inline Double_t PamVMCPrimaryGenerator::KinEToMomentum(Double_t E0)
159	{ return Sqrt(E0E0+2.E0*fmass); }
160
161	#endif //PAMVMC_PRIMARY_GENERATOR_H
162