| 1 |
// C/C++ headers |
| 2 |
// |
| 3 |
#include <fstream> |
| 4 |
#include <string.h> |
| 5 |
#include <iostream> |
| 6 |
#include <stdio.h> |
| 7 |
#include <ctype.h> |
| 8 |
// |
| 9 |
// ROOT headers |
| 10 |
// |
| 11 |
|
| 12 |
#include <TTree.h> |
| 13 |
|
| 14 |
#include <TSystem.h> |
| 15 |
#include <TStopwatch.h> |
| 16 |
#include <TSystemDirectory.h> |
| 17 |
#include <TGraphAsymmErrors.h> |
| 18 |
#include <TString.h> |
| 19 |
#include <TFile.h> |
| 20 |
#include <TClass.h> |
| 21 |
#include <TGraph.h> |
| 22 |
#include <TLine.h> |
| 23 |
#include <TH1F.h> |
| 24 |
#include <TMath.h> |
| 25 |
#include <TRandom3.h> |
| 26 |
#include <TF1.h> |
| 27 |
// |
| 28 |
// RunInfo header |
| 29 |
// |
| 30 |
#include <RunInfo.h> |
| 31 |
|
| 32 |
#include <PamLevel2.h> |
| 33 |
|
| 34 |
#include <TClonesArray.h> |
| 35 |
#include <PamVMCPrimaryGenerator.h> |
| 36 |
|
| 37 |
using namespace std; |
| 38 |
using TMath::Pi; |
| 39 |
//This program generate trajectories, which are going in acceptance |
| 40 |
|
| 41 |
|
| 42 |
void Usage (); |
| 43 |
|
| 44 |
|
| 45 |
|
| 46 |
Bool_t IsInside(Trajectory *t, Double_t fiducial, |
| 47 |
float *xmin, float *xmax, float *ymin, float *ymax){ |
| 48 |
|
| 49 |
for(Int_t i = 0; i<t->npoint; i++){ |
| 50 |
if( !(((xmin[i]+fiducial)<t->x[i])*(t->x[i]<(xmax[i]-fiducial))) ) return kFALSE; |
| 51 |
if( !(((ymin[i]+fiducial)<t->y[i])*(t->y[i]<(ymax[i]-fiducial))) ) return kFALSE; |
| 52 |
//cout<<"XMIN:"<<xmin[i] |
| 53 |
//<<" XMAX:"<<xmax[i] |
| 54 |
//<<" YMIN:"<<ymin[i] |
| 55 |
//<<" YMAX:"<<ymax[i] |
| 56 |
//<<" Z:"<<TrkParams::zGF[i]<<endl; |
| 57 |
} |
| 58 |
//t->Dump(); |
| 59 |
|
| 60 |
return kTRUE; |
| 61 |
} |
| 62 |
|
| 63 |
|
| 64 |
|
| 65 |
void DoHistogramFromGraph(TGraph* g, TH1F* htemplate, Double_t factor){ |
| 66 |
//TH1F* hout = (TH1F*) htemplate->Clone(name); |
| 67 |
htemplate->Reset(); |
| 68 |
int n = htemplate->GetNbinsX(); |
| 69 |
for(int bin=1; bin<=n; bin++) { |
| 70 |
htemplate->SetBinContent(bin, g->Eval(htemplate->GetBinCenter(bin))*factor); |
| 71 |
Double_t *gx = g->GetX(); |
| 72 |
int p = 0; |
| 73 |
Double_t x = htemplate->GetBinCenter(bin); |
| 74 |
for(p=1; p<g->GetN(); p++) |
| 75 |
if((gx[p-1] <= x) && (gx[p] > x)) break; |
| 76 |
Double_t gey = g->GetErrorY(p-1)*factor; |
| 77 |
//cout<<"gx: "<< hout->GetBinCenter(bin)<<"gy: "<<hout->GetBinContent(bin)<<endl; |
| 78 |
htemplate->SetBinError(bin, gey); |
| 79 |
//cout<<"gx: "<< htemplate->GetBinCenter(bin)<<"gy: "<<htemplate->GetBinContent(bin)<<endl; |
| 80 |
} |
| 81 |
} |
| 82 |
|
| 83 |
void BinLogX(TH1*h){ |
| 84 |
// |
| 85 |
// Modify histogram h to create an histogram logarithmically binned |
| 86 |
// along the X axis. It works on normally defined, not |
| 87 |
// filled,histograms. |
| 88 |
// |
| 89 |
TAxis *axis = h->GetXaxis(); |
| 90 |
int bins = axis->GetNbins(); |
| 91 |
|
| 92 |
Axis_t xmin = axis->GetXmin(); |
| 93 |
Axis_t xmax = axis->GetXmax(); |
| 94 |
|
| 95 |
if (xmin==0) |
| 96 |
{ |
| 97 |
xmin = 1e-6; |
| 98 |
cout<<"Warning: low edge of histogram set from 0 to 1e-6"<<endl; |
| 99 |
} |
| 100 |
Axis_t from = TMath::Log10(xmin); |
| 101 |
Axis_t to = TMath::Log10(xmax); |
| 102 |
|
| 103 |
// Axis_t width = (xmax - xmin) / bins; |
| 104 |
Axis_t width = (to - from) / bins; |
| 105 |
Axis_t *new_bins = new Axis_t[bins + 1]; |
| 106 |
|
| 107 |
for (int i = 0; i <= bins; i++) { |
| 108 |
new_bins[i] = TMath::Power(10, from + i * width); |
| 109 |
} |
| 110 |
axis->Set(bins, new_bins); |
| 111 |
delete new_bins; |
| 112 |
} |
| 113 |
|
| 114 |
|
| 115 |
int main(int argc, char * argv[] ){ |
| 116 |
TStopwatch timer; |
| 117 |
timer.Start(); |
| 118 |
|
| 119 |
TString DataPATH = gSystem->Getenv("PWD"); |
| 120 |
|
| 121 |
|
| 122 |
TString inputfilename, outputfilename; |
| 123 |
TFile* rfin = 0; |
| 124 |
TFile* rfout = 0; |
| 125 |
TTree* primaries = 0; |
| 126 |
TClonesArray* PRIM = NULL; |
| 127 |
PamVMCPrimary *p = 0; |
| 128 |
Double_t fiducial, Z, Rmin, Rmax, R; |
| 129 |
Int_t PDG; |
| 130 |
fiducial=Z=Rmin=Rmax=R=0.; |
| 131 |
Int_t nev; |
| 132 |
TF1 *frigflat = 0; |
| 133 |
TH1F * flux = NULL; |
| 134 |
TH1F * evDist = NULL; |
| 135 |
|
| 136 |
enum mode {SPECTRA, FLATBIN, SINGLEPOINT}; |
| 137 |
mode md; |
| 138 |
|
| 139 |
TString innf(argv[5]); |
| 140 |
switch(argc){ |
| 141 |
case 7: |
| 142 |
md = SINGLEPOINT; |
| 143 |
break; |
| 144 |
case 8: |
| 145 |
md = FLATBIN; |
| 146 |
break; |
| 147 |
case 9: |
| 148 |
md = SPECTRA; |
| 149 |
break; |
| 150 |
default: |
| 151 |
Usage(); |
| 152 |
return 0; |
| 153 |
break; |
| 154 |
} |
| 155 |
|
| 156 |
nev = atoi(argv[1]); |
| 157 |
Z = atof(argv[2]); |
| 158 |
PDG = atoi(argv[3]); |
| 159 |
fiducial = atof(argv[4]); |
| 160 |
|
| 161 |
|
| 162 |
TSystemDirectory *workdir = 0; |
| 163 |
|
| 164 |
switch(md){ |
| 165 |
case SPECTRA: |
| 166 |
Rmin = atof(argv[5]); |
| 167 |
Rmax = atof(argv[6]); |
| 168 |
inputfilename = TString(argv[7]); |
| 169 |
outputfilename = TString(argv[8]); |
| 170 |
flux = new TH1F("flux", "flux", 2000, Rmin, Rmax); |
| 171 |
workdir=new TSystemDirectory("workdir",DataPATH); |
| 172 |
rfin = new TFile(inputfilename); |
| 173 |
if (rfin){ |
| 174 |
cout<<"Generating spectra from file: "<<inputfilename<<endl; |
| 175 |
cout << "Rmin="<<Rmin<<" GV,"<<" Rmax="<<Rmax<<" GV"<<endl; |
| 176 |
|
| 177 |
// TList * outlist = (TList*)rfin->Get("outlist"); |
| 178 |
//flux = (TH1F*)outlist->At(5); //CHANGE THIS IF NEED |
| 179 |
TGraphAsymmErrors *gr1 = (TGraphAsymmErrors*)rfin->Get("gflux"); |
| 180 |
//TH1F * temp = new TH1F("template", "template", 40, 0.4, 250.); |
| 181 |
BinLogX(flux); |
| 182 |
DoHistogramFromGraph(gr1, flux, 1000.); |
| 183 |
rfin->Close(); |
| 184 |
//cout<<"test"<<flux->GetBinCenter(10)<<endl; |
| 185 |
} |
| 186 |
evDist = (TH1F*) flux->Clone("evDist"); |
| 187 |
for(int ib=0; ib<evDist->GetNbinsX(); ib++) evDist->SetBinContent(ib+1, evDist->GetBinContent(ib+1)*evDist->GetBinWidth(ib+1)); |
| 188 |
break; |
| 189 |
case FLATBIN: |
| 190 |
Rmin = atof(argv[5]); |
| 191 |
Rmax = atof(argv[6]); |
| 192 |
if (Rmin == Rmax) cout<<"ERROR: Rmin = Rmax"<<endl; |
| 193 |
frigflat = new TF1("fmom","1",Rmin,Rmax); |
| 194 |
outputfilename = TString(argv[7]); |
| 195 |
cout<<"Generating flat spectra, Rmin="<<Rmin<<" GV,"<<" Rmax="<<Rmax<<" GV"<<endl; |
| 196 |
break; |
| 197 |
case SINGLEPOINT: |
| 198 |
R = atof(argv[5]); |
| 199 |
outputfilename = TString(argv[6]); |
| 200 |
cout<<"Generating fixed rigidity, R="<<R<<" GV"<<endl; |
| 201 |
break; |
| 202 |
default: |
| 203 |
break; |
| 204 |
} |
| 205 |
|
| 206 |
//Creating output root-file |
| 207 |
rfout = new TFile(outputfilename.Data(),"recreate"); |
| 208 |
primaries = new TTree("hits","hits"); |
| 209 |
TTree::SetMaxTreeSize(1000*Long64_t(2000000000)); |
| 210 |
PRIM = new TClonesArray("PamVMCPrimary"); |
| 211 |
if(!primaries->GetBranch("PRIM")) |
| 212 |
primaries ->Branch("PRIM","TClonesArray", &PRIM); |
| 213 |
else |
| 214 |
primaries->GetBranch("PRIM")->SetAddress(&PRIM); |
| 215 |
|
| 216 |
|
| 217 |
|
| 218 |
|
| 219 |
//Loading field, Initialization |
| 220 |
|
| 221 |
TRandom3 *ran = new TRandom3(0); |
| 222 |
TF1 *ftheta = new TF1("ftheta","sin(x)*cos(x)",0.,acos(-1.)/4.); |
| 223 |
ftheta->SetNpx(99999); |
| 224 |
|
| 225 |
cout<<"Load Magnetic Field... "<<endl; |
| 226 |
TrkTrack *track = new TrkTrack(); |
| 227 |
track->LoadField(getenv("PAM_FIELD")); |
| 228 |
track->SetTrackingMode(0); |
| 229 |
cout<<"Magnetic Field Loaded. "<<endl; |
| 230 |
|
| 231 |
//Start |
| 232 |
|
| 233 |
cout<<"Nevents: "<<nev<<", FI DUCIAL: "<<fiducial<<" (cm), Z:"<<Z<<", PDG:"<<PDG<<endl; |
| 234 |
|
| 235 |
#define ZCAL 13.05 //cm |
| 236 |
#define CALS3 10.639 //cm |
| 237 |
#define ZSPEBP 2.97 //cm |
| 238 |
#define HORMG 22.57 //cm |
| 239 |
|
| 240 |
|
| 241 |
#define XMINGEN -30. |
| 242 |
#define XMAXGEN 30. |
| 243 |
#define YMINGEN -30. |
| 244 |
#define YMAXGEN 30. |
| 245 |
#define ZINP 110. |
| 246 |
|
| 247 |
setprecision(8); |
| 248 |
|
| 249 |
float *xmin = TrkParams::xGF_min; |
| 250 |
float *xmax = TrkParams::xGF_max; |
| 251 |
float *ymin = TrkParams::yGF_min; |
| 252 |
float *ymax = TrkParams::yGF_max; |
| 253 |
|
| 254 |
Float_t zin[TrkParams::nGF]; |
| 255 |
for (Int_t igf=0; igf<TrkParams::nGF; igf++) zin[igf]=TrkParams::zGF[igf]; |
| 256 |
Trajectory* t = new Trajectory(TrkParams::nGF,zin); |
| 257 |
Float_t al[5]; |
| 258 |
|
| 259 |
//Main LOOP |
| 260 |
|
| 261 |
Int_t goodev = 0; |
| 262 |
Int_t totev = 0; |
| 263 |
Double_t mom = 0.; |
| 264 |
Double_t x,y,z,theta,phi; |
| 265 |
while (goodev<nev){ |
| 266 |
|
| 267 |
x = ran->Uniform(XMINGEN, XMAXGEN); |
| 268 |
y = ran->Uniform(YMINGEN, YMAXGEN); |
| 269 |
z = ZINP-ZCAL-CALS3-ZSPEBP-HORMG; |
| 270 |
theta = ftheta->GetRandom(); |
| 271 |
phi = ran->Uniform(0.,2.*acos(-1.)); |
| 272 |
|
| 273 |
switch(md){ |
| 274 |
case SPECTRA: |
| 275 |
mom = Z*evDist->GetRandom(); |
| 276 |
break; |
| 277 |
case FLATBIN: |
| 278 |
mom = Z*frigflat->GetRandom(); |
| 279 |
break; |
| 280 |
case SINGLEPOINT: |
| 281 |
mom = Z*R; |
| 282 |
break; |
| 283 |
default: |
| 284 |
break; |
| 285 |
} |
| 286 |
|
| 287 |
al[0] = x; |
| 288 |
al[1] = y; |
| 289 |
al[2] = sin(theta); |
| 290 |
al[3] = phi; |
| 291 |
al[4] = Z/mom; |
| 292 |
|
| 293 |
t->DoTrack(al, z); |
| 294 |
|
| 295 |
if(IsInside(t,fiducial, xmin, xmax, ymin, ymax)){ |
| 296 |
p = (PamVMCPrimary*) PRIM->New(0); |
| 297 |
|
| 298 |
p->fX0 = x; |
| 299 |
p->fY0 = y; |
| 300 |
p->fZ0 = ZINP; |
| 301 |
p->fTHETA = theta; |
| 302 |
p->fPHI = phi; |
| 303 |
p->fP0 = mom; |
| 304 |
p->fPDG=PDG; //He-4 |
| 305 |
|
| 306 |
primaries->Fill(); |
| 307 |
goodev++; |
| 308 |
|
| 309 |
PRIM->Delete(); |
| 310 |
} |
| 311 |
totev++; |
| 312 |
if( ((totev%10000)==0) || (goodev==nev) ) cout<<"-> "<<totev<<" "<<goodev<<endl; |
| 313 |
} |
| 314 |
primaries->Write(); |
| 315 |
if(md==SPECTRA){ |
| 316 |
rfout->WriteTObject(flux); |
| 317 |
rfout->WriteTObject(evDist); |
| 318 |
} |
| 319 |
rfout->Close(); |
| 320 |
|
| 321 |
Double_t gf = Double_t(goodev)/Double_t(totev)*(XMAXGEN-XMINGEN)*(YMAXGEN-YMINGEN)*TMath::Pi()/2.; |
| 322 |
Double_t err_gf = (Double_t)1./TMath::Sqrt(goodev)*gf; |
| 323 |
cout<<"GF="<<gf<<" +- "<<err_gf<<endl; |
| 324 |
|
| 325 |
//delete primaries; |
| 326 |
delete p; |
| 327 |
delete ran; |
| 328 |
delete ftheta; |
| 329 |
delete track; |
| 330 |
delete t; |
| 331 |
return 0; |
| 332 |
} |
| 333 |
|
| 334 |
void Usage(){ |
| 335 |
cout<<"NAME"<<endl; |
| 336 |
cout<<" sp_gen - spectra generator"<<endl; |
| 337 |
cout<<"SYNOPSIS"<<endl; |
| 338 |
cout<<" sp_gen NUMPAR Z PDG FIDUCIAL(cm) RIGIDITY(GV) OUTFILE"<<endl; |
| 339 |
cout<<" sp_gen NUMPAR Z PDG FIDUCIAL(cm) RIGMIN(GV) RIGMAX(GV) OUTFILE"<<endl; |
| 340 |
cout<<" sp_gen NUMPAR Z PDG FIDUCIAL(cm) RIGMIN(GV) RIGMAX(GV) INPUTFILE OUTFILE"<<endl; |
| 341 |
cout<<"DISCRIPTION"<<endl; |
| 342 |
cout<<" Generate tracks above the apparatus "<<endl; |
| 343 |
cout<<" isotropically and save initial kinematic parameters (X0,Y0,Z0=110., P, Theta, Phi)"<<endl; |
| 344 |
cout<<" good tracks in output root-file. After this, file could be used to do simulations."<<endl; |
| 345 |
cout<<" The condition to accept tracks is satisfaction to 14 acceptance planes, defined by"<<endl; |
| 346 |
cout<<" TrkParams in Level2. This tool is based on DoTrack method of Level2. "<<endl; |
| 347 |
cout<<" To save data in root-file PrimaryInfo class was developed and compiled in library. "<<endl; |
| 348 |
cout<<" If later you wish to develop your own script to work with data stored threre, "<<endl; |
| 349 |
cout<<" you need just load libPrimaryInfo.so shared library"<<endl; |
| 350 |
cout<<"NOTE"<<endl; |
| 351 |
cout<<" Use http://pdg.lbl.gov/2007/reviews/montecarlorpp.pdf document to define your"<<endl; |
| 352 |
cout<<" particle according PDG numeration scheme."<<endl; |
| 353 |
} |
Parent Directory
| 
Revision Log