| 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 |
//#include <TObjectTable.h> |
| 15 |
|
| 16 |
#include <PamLevel2.h> |
| 17 |
|
| 18 |
#include "PamRootManager.h" |
| 19 |
#include "PamVMC2PiAux.h" |
| 20 |
#include "PamVMCPrimaryGF.h" |
| 21 |
|
| 22 |
class TVirtualMCStack; |
| 23 |
|
| 24 |
using std::cout; |
| 25 |
using std::endl; |
| 26 |
using TMath::Sqrt; |
| 27 |
|
| 28 |
class PamVMCPrimary : public TObject |
| 29 |
{ |
| 30 |
|
| 31 |
public: |
| 32 |
|
| 33 |
PamVMCPrimary(): fID(0), fPDG(0), fX0(0.), fY0(0.), fZ0(0.), fTHETA(0.), fPHI(0.), fP0(0.), fGOOD(kFALSE), fXTOL(-100.), fYTOL(-100.){ |
| 34 |
|
| 35 |
|
| 36 |
for (Int_t j = 0; j<14; j++){ |
| 37 |
fCR_X[j] = -100.; |
| 38 |
fCR_Y[j] = -100.; |
| 39 |
} |
| 40 |
}; |
| 41 |
|
| 42 |
virtual ~PamVMCPrimary(){ /*gObjectTable->Remove(this);*/}; |
| 43 |
|
| 44 |
UInt_t fID; |
| 45 |
Int_t fPDG; |
| 46 |
Double_t fX0, fY0, fZ0; |
| 47 |
Double_t fTHETA, fPHI; |
| 48 |
Double_t fP0; |
| 49 |
Bool_t fGOOD; |
| 50 |
Double_t fXTOL, fYTOL; |
| 51 |
Double_t fCR_X [14]; |
| 52 |
Double_t fCR_Y [14]; |
| 53 |
|
| 54 |
|
| 55 |
void Clean() { fID=0; fPDG=0; fX0=fY0=fZ0=fTHETA=fPHI=fP0=0.; fGOOD=kFALSE; fXTOL=fYTOL=-100.; |
| 56 |
for (Int_t j = 0; j<14; j++){ |
| 57 |
fCR_X[j] = -100.; |
| 58 |
fCR_Y[j] = -100.; |
| 59 |
} |
| 60 |
|
| 61 |
} |
| 62 |
|
| 63 |
void Clear(const Option_t * =""){ this->~PamVMCPrimary(); } |
| 64 |
|
| 65 |
void Print(const Option_t * ="") const |
| 66 |
{ cout<<"PRIMARY particle "<<fID<<" infromation:"<<endl; |
| 67 |
cout<<"Pdg="<<fPDG<<endl; |
| 68 |
cout<<"Position: "<<"("<<fX0<<","<<fY0<<","<<fZ0<<")"<<endl; |
| 69 |
cout<<"P0, Theta, Phi: "<<fP0<<","<<fTHETA<<","<<fPHI<<endl; |
| 70 |
cout<<"GOOD Single Track:"<<fGOOD<<endl; |
| 71 |
cout<<"Xtol= "<<fXTOL<<" Ytol= "<<fYTOL<<endl; |
| 72 |
cout<<" Coordinates of acceptance planes crossing:"<<endl; |
| 73 |
for(Int_t i = 0; i< 14; i++){ |
| 74 |
cout<<"i="<<i<<" X="<<fCR_X[i]<<" Y="<<fCR_Y[i]<<endl; |
| 75 |
} |
| 76 |
} |
| 77 |
|
| 78 |
ClassDef(PamVMCPrimary,1); |
| 79 |
}; |
| 80 |
|
| 81 |
PamVMCPrimary & operator+=(PamVMCPrimary &a, const PamVMCPrimary &b); |
| 82 |
|
| 83 |
class PamVMCPrimaryGenerator : public TObject |
| 84 |
{ |
| 85 |
public: |
| 86 |
PamVMCPrimaryGenerator(TVirtualMCStack* stack); |
| 87 |
PamVMCPrimaryGenerator(); |
| 88 |
|
| 89 |
virtual ~PamVMCPrimaryGenerator(); |
| 90 |
|
| 91 |
// methods |
| 92 |
|
| 93 |
void GeneratePrimary(); |
| 94 |
// set methods |
| 95 |
void SetParticle(Int_t pdg); |
| 96 |
void SetParticleID(UInt_t id){ fprim.fID = id; }; |
| 97 |
void SetPosition(Double_t vx, Double_t vy, Double_t vz) |
| 98 |
{fprim.fX0=vx; fprim.fY0=vy; fprim.fZ0=vz;}; |
| 99 |
|
| 100 |
void SetKinEnergy(Double_t kinEnergy) {fprim.fP0=KinEToMomentum(kinEnergy);}; |
| 101 |
void SetRigidity(Double_t rigidity) {fprim.fP0=RigToMomentum(rigidity);}; |
| 102 |
void SetMomentum(Double_t momentum) {fprim.fP0=momentum; }; |
| 103 |
|
| 104 |
void SetDirection(Double_t theta, Double_t phi) {fprim.fTHETA=theta; fprim.fPHI=phi; }; |
| 105 |
void SetMomentum(Double_t px, Double_t py, Double_t pz); |
| 106 |
|
| 107 |
// gen methods |
| 108 |
|
| 109 |
void GenPosition(Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax, Double_t zmin, Double_t zmax){ |
| 110 |
SetPosition(frandom->Uniform(xmin,xmax),frandom->Uniform(ymin,ymax),frandom->Uniform(zmin,zmax)); } |
| 111 |
|
| 112 |
void GenDirection(Double_t thetamin, Double_t thetamax, Double_t phimin, Double_t phimax){ |
| 113 |
Double_t theta = TMath::ASin(Sqrt( frandom->Uniform( TMath::ASin(thetamin), TMath::ASin(thetamax)))); |
| 114 |
Double_t phi = TMath::TwoPi() * frandom->Uniform( phimin, phimax ); |
| 115 |
SetDirection(theta, phi); |
| 116 |
} |
| 117 |
|
| 118 |
void GenPosAng(TVector3 &pos, TVector2 &ang) { f2PiAux->GenPosAng(pos,ang); } |
| 119 |
|
| 120 |
void Gen2PiPosAng(){ |
| 121 |
TVector3 pos; |
| 122 |
TVector2 ang; |
| 123 |
GenPosAng(pos,ang); |
| 124 |
SetPosition( pos.X(), pos.Y(), pos.Z() ); |
| 125 |
SetDirection( ang.X(), ang.Y() ); |
| 126 |
} |
| 127 |
|
| 128 |
void ExtrapolateTrajectory(){ //Added by V. Formato on 2012/06/19 |
| 129 |
Float_t zin[TrkParams::nGF]; |
| 130 |
for (Int_t igf=0; igf<TrkParams::nGF; igf++) zin[igf]=TrkParams::zGF[igf]; |
| 131 |
Trajectory* t = new Trajectory(TrkParams::nGF,zin); |
| 132 |
Float_t al[5]; |
| 133 |
al[0] = fprim.fX0; |
| 134 |
al[1] = fprim.fY0; |
| 135 |
al[2] = sin(fprim.fTHETA); |
| 136 |
al[3] = fprim.fPHI; |
| 137 |
al[4] = fcharge/fprim.fP0; |
| 138 |
|
| 139 |
t->DoTrack(al, fprim.fZ0 - 13.05 - 10.639 - 2.97 - 22.57); |
| 140 |
// t->DoTrack(al, fprim.fZ0); |
| 141 |
|
| 142 |
for(Int_t igf=0; igf<NGF; igf++){ |
| 143 |
fprim.fCR_X[igf] = t->x[igf]; |
| 144 |
fprim.fCR_Y[igf] = t->y[igf]; |
| 145 |
} |
| 146 |
|
| 147 |
delete t; |
| 148 |
|
| 149 |
} |
| 150 |
|
| 151 |
Bool_t ToBeSimulated(){ |
| 152 |
|
| 153 |
float *xmin = TrkParams::xGF_min; |
| 154 |
float *xmax = TrkParams::xGF_max; |
| 155 |
float *ymin = TrkParams::yGF_min; |
| 156 |
float *ymax = TrkParams::yGF_max; |
| 157 |
|
| 158 |
for(Int_t i=0; i<NGF; i++){ |
| 159 |
if( !((xmin[i]<fprim.fCR_X[i])&&(fprim.fCR_X[i]<xmax[i])) ) return kFALSE; |
| 160 |
if( !((ymin[i]<fprim.fCR_Y[i])&&(fprim.fCR_Y[i]<ymax[i])) ) return kFALSE; |
| 161 |
} |
| 162 |
|
| 163 |
return kTRUE; |
| 164 |
} |
| 165 |
|
| 166 |
|
| 167 |
//flat spectra generator |
| 168 |
void GenSpe(Double_t PEmin, Double_t PEmax, Bool_t isEnergy=kFALSE); |
| 169 |
//power law spectra, gamma - differential spectral index |
| 170 |
void GenSpe(Double_t PEmin, Double_t PEmax, Double_t gamma, Bool_t isEnergy=kFALSE); |
| 171 |
void GenSpe_Flat(Double_t PEmin, Double_t PEmax, Double_t gamma, Bool_t isEnergy=kFALSE); |
| 172 |
void GenSpe_3par(Double_t PEmin, Double_t PEmax, Double_t a, Double_t b, Double_t c); |
| 173 |
void GenSphericalPhiThe(); // cecilia // flusso istropo // sets phi the |
| 174 |
void GenSphPhiThe(Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax, |
| 175 |
Double_t zmin, Double_t zmax); // flusso isotropo entro S2 S3 |
| 176 |
// sets position and phi the |
| 177 |
|
| 178 |
|
| 179 |
// get methods |
| 180 |
Int_t GetParticleID(){ return fprim.fID; }; |
| 181 |
Int_t GetParticle(){ return fprim.fPDG; }; |
| 182 |
void GetPositon(TVector3 & v){ v.SetXYZ(fprim.fX0,fprim.fY0,fprim.fZ0); }; |
| 183 |
void GetPosition(Double_t & X0, Double_t & Y0,Double_t & Z0) |
| 184 |
{X0=fprim.fX0; Y0=fprim.fY0; Z0=fprim.fZ0;}; |
| 185 |
|
| 186 |
Double_t GetKinEnergy() { return MomentumToKinE(fprim.fP0); }; |
| 187 |
Double_t GetRigidity() { return MomentumToRig(fprim.fP0); }; |
| 188 |
Double_t GetMomentum() { return fprim.fP0; }; |
| 189 |
|
| 190 |
Double_t Get2PiArea() { return f2PiAux->GetArea(); }; |
| 191 |
|
| 192 |
void GetDirection(TVector2 & v) {v.Set(fprim.fTHETA,fprim.fPHI);} |
| 193 |
void GetDirection(Double_t &theta, Double_t & phi) { theta = fprim.fTHETA; phi = fprim.fPHI;} |
| 194 |
|
| 195 |
Bool_t Getgood(){ return fprim.fGOOD; }; |
| 196 |
|
| 197 |
void GetTol(Double_t& xtol, Double_t& ytol){ xtol=fprim.fXTOL; ytol=fprim.fYTOL; }; |
| 198 |
|
| 199 |
void GetCR_GF(Double_t& xcr, Double_t& ycr, Int_t id_GF){ |
| 200 |
if(id_GF<14){ |
| 201 |
xcr=fprim.fCR_X[id_GF]; |
| 202 |
ycr=fprim.fCR_Y[id_GF]; |
| 203 |
} |
| 204 |
} |
| 205 |
|
| 206 |
Int_t GetNgood(){ return fngood; } |
| 207 |
|
| 208 |
//initialize random |
| 209 |
void SetRandom(TRandom* random){ |
| 210 |
frandom = random; |
| 211 |
f2PiAux->SetRandom(random); |
| 212 |
} |
| 213 |
|
| 214 |
//work with collection of primaries |
| 215 |
|
| 216 |
void Register(){ |
| 217 |
PamRootManager::Instance()-> |
| 218 |
Register("PRIM","TClonesArray", &fprimColl); |
| 219 |
} |
| 220 |
|
| 221 |
void ClearPrimCol(){ |
| 222 |
fevno=0; |
| 223 |
(*fprimColl).Clear("C"); |
| 224 |
} |
| 225 |
//all methods to operate with PamVMCPrimaryGF object |
| 226 |
//after each step of primary |
| 227 |
void SetZeroStep(Double_t z_curr, Double_t x_curr, Double_t y_curr){ |
| 228 |
fzprev = z_curr; |
| 229 |
fxprev = x_curr; |
| 230 |
fyprev = y_curr; |
| 231 |
} |
| 232 |
void CheckCrossing(Double_t z_curr, Double_t x_curr, Double_t y_curr){ |
| 233 |
fpgf->CheckCrossing(fzprev,fxprev,fyprev, z_curr, x_curr, y_curr); |
| 234 |
SetZeroStep(z_curr,x_curr, y_curr); |
| 235 |
} |
| 236 |
//after last step of primary |
| 237 |
void CheckInsideAcceptance(){ |
| 238 |
// for(Int_t i = 0; i<NGF; i++) fprim.fCR_X[i] = fprim.fCR_Y[i] = fprim.fXTOL = fprim.fYTOL = -100.; //zeroing |
| 239 |
fprim.fGOOD = fpgf->IsInsideAcceptance(fprim.fCR_X, fprim.fCR_Y, fprim.fXTOL, fprim.fYTOL); |
| 240 |
fpgf->Clear(); |
| 241 |
if(fprim.fGOOD) fngood++; |
| 242 |
//Here I'm getting what I have in collection and asume that I have only 1 primary here |
| 243 |
//If I wish to launch more in 1 event I need change this stuff |
| 244 |
PamVMCPrimary* tmp = (PamVMCPrimary*)fprimColl->At(0); |
| 245 |
if(tmp) *tmp = fprim; |
| 246 |
} |
| 247 |
|
| 248 |
Double_t MomentumToKinE(Double_t P0); |
| 249 |
Double_t KinEToMomentum(Double_t E0); |
| 250 |
Double_t MomentumToRig(Double_t P0) { return TMath::Abs(P0/fcharge); }; |
| 251 |
Double_t RigToMomentum(Double_t R0){ return TMath::Abs(R0*fcharge); }; |
| 252 |
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);}; |
| 253 |
|
| 254 |
private: |
| 255 |
// methods |
| 256 |
|
| 257 |
|
| 258 |
|
| 259 |
// data members |
| 260 |
TVirtualMCStack* fStack; |
| 261 |
Int_t fevno; |
| 262 |
Int_t fngood; //particles inside nominal acceptance |
| 263 |
PamVMCPrimary fprim; |
| 264 |
Double_t fmass; |
| 265 |
Double_t fcharge; |
| 266 |
Double_t fxprev, fyprev, fzprev; // coordinates of primary on prev. step |
| 267 |
TClonesArray* fprimColl; |
| 268 |
TRandom* frandom; // Class is not a owner of this object |
| 269 |
TF1* ftheta; // To generate sherical distributhin in theta-angle |
| 270 |
PamVMC2PiAux* f2PiAux; //2 Pi generator aux stuff |
| 271 |
PamVMCPrimaryGF* fpgf; // Object which |
| 272 |
|
| 273 |
|
| 274 |
ClassDef(PamVMCPrimaryGenerator,1) //PamVMCPrimaryGenerator |
| 275 |
}; |
| 276 |
|
| 277 |
// inline functions |
| 278 |
|
| 279 |
|
| 280 |
inline Double_t PamVMCPrimaryGenerator::MomentumToKinE(Double_t P0) |
| 281 |
{ return Sqrt(P0*P0+fmass*fmass)-fmass; } |
| 282 |
|
| 283 |
inline Double_t PamVMCPrimaryGenerator::KinEToMomentum(Double_t E0) |
| 284 |
{ return Sqrt(E0*E0+2.*E0*fmass); } |
| 285 |
|
| 286 |
#endif //PAMVMC_PRIMARY_GENERATOR_H |
| 287 |
|
Parent Directory
| 
Revision Log