SND@LHC Software
Loading...
Searching...
No Matches
Public Member Functions | Protected Member Functions | Protected Attributes | Private Member Functions | List of all members
NuageGenerator Class Reference

#include <NuageGenerator.h>

Inheritance diagram for NuageGenerator:
Inheritance graph
[legend]
Collaboration diagram for NuageGenerator:
Collaboration graph
[legend]

Public Member Functions

 NuageGenerator ()
 
virtual ~NuageGenerator ()
 
Bool_t ReadEvent (FairPrimaryGenerator *)
 
virtual Bool_t Init (const char *, int)
 
virtual Bool_t Init (const char *)
 
Int_t GetNevents ()
 
void NuOnly ()
 
void SetPositions (Double_t zTa, Double_t zS=-3400., Double_t zE=2650., Double_t xS=0, Double_t xE=0, Double_t yS=0, Double_t yE=0)
 
void AddBox (TVector3 dVec, TVector3 box)
 
void EnableExternalDecayer (Bool_t value)
 

Protected Member Functions

 ClassDef (NuageGenerator, 2)
 

Protected Attributes

Bool_t fExtDecayer
 
Double_t Yvessel
 
Double_t Xvessel
 
Double_t Lvessel
 
Double_t ztarget
 
Double_t startZ
 
Double_t endZ
 
Double_t startX
 
Double_t endX
 
Double_t startY
 
Double_t endY
 
Float_t Ev
 
Float_t pxv
 
Float_t pyv
 
Float_t pzv
 
Float_t El
 
Float_t pxl
 
Float_t pyl
 
Float_t pzl
 
Float_t vtxx
 
Float_t vtxy
 
Float_t vtxz
 
Float_t vtxt
 
Float_t vtxx2
 
Float_t vtxy2
 
Float_t vtxz2
 
Float_t vtxt2
 
Float_t vtxx4
 
Float_t vtxy3
 
Float_t vtxz3
 
Float_t vtxt3
 
Float_t vtxx3
 
Float_t vtxy4
 
Float_t vtxz4
 
Float_t vtxt4
 
Float_t pxf [500]
 
Float_t pyf [500]
 
Float_t pzf [500]
 
Float_t pdgf [500]
 
Float_t pxf2 [500]
 
Float_t pyf2 [500]
 
Float_t pzf2 [500]
 
Float_t pdgf2 [500]
 
Float_t pxf3 [500]
 
Float_t pyf3 [500]
 
Float_t pzf3 [500]
 
Float_t pdgf3 [500]
 
Float_t pxf4 [500]
 
Float_t pyf4 [500]
 
Float_t pzf4 [500]
 
Float_t pdgf4 [500]
 
std::vector< TVector3 > dVecs
 
std::vector< TVector3 > boxs
 
Int_t cc
 
Int_t nf
 
Int_t nf2
 
Int_t nf3
 
Int_t nf4
 
Int_t neu
 
Int_t parent2
 
Int_t parent3
 
Int_t parent4
 
Int_t fNvtx
 
FairLogger * fLogger
 
TFile * fInputFile
 don't make it persistent, magic ROOT command
 
TTree * fTree
 
int fNevents
 
Int_t fn
 
bool fFirst
 
bool fNuOnly
 
Double_t fznu
 
Double_t fXnu
 
Double_t fYnu
 
Double_t fEntrDz_inner
 
Double_t fEntrDz_outer
 
Double_t fEntrZ_inner
 
Double_t fEntrZ_outer
 
Double_t fEntrA
 
Double_t fEntrB
 
Double_t fL1z
 
Double_t fScintDz
 

Private Member Functions

std::vector< double > Rotate (Double_t x, Double_t y, Double_t z, Double_t px, Double_t py, Double_t pz)
 
Double_t MeanMaterialBudget (const Double_t *start, const Double_t *end, Double_t *mparam)
 

Detailed Description

Definition at line 14 of file NuageGenerator.h.

Constructor & Destructor Documentation

◆ NuageGenerator()

NuageGenerator::NuageGenerator ( )

default constructor

Definition at line 24 of file NuageGenerator.cxx.

24 {
25 cout<<endl;
26 cout<<"**************************************"<<endl;
27 cout<<"WARNING!!!!!!"<<endl;
28 cout<<"Bug on decays found and solved for GenieGenerator still to be fixed (25.08.2017)" <<endl;
29 cout<<"**************************************"<<endl;
30 cout<<endl;
31}

◆ ~NuageGenerator()

NuageGenerator::~NuageGenerator ( )
virtual

destructor

Definition at line 295 of file NuageGenerator.cxx.

296{
297 fInputFile->Close();
298 fInputFile->Delete();
299 delete fInputFile;
300}
NuageGenerator::fInputFile
TFile * fInputFile
don't make it persistent, magic ROOT command
Definition NuageGenerator.h:69

Member Function Documentation

◆ AddBox()

void NuageGenerator::AddBox ( TVector3  dVec,
TVector3  box 
)

◆ ClassDef()

NuageGenerator::ClassDef ( NuageGenerator  ,
2   
)
protected

◆ EnableExternalDecayer()

void NuageGenerator::EnableExternalDecayer ( Bool_t  value)
inline

Definition at line 40 of file NuageGenerator.h.

40{fExtDecayer=value;}

◆ GetNevents()

Int_t NuageGenerator::GetNevents ( )

Definition at line 596 of file NuageGenerator.cxx.

597{
598 return fNevents;
599}

◆ Init() [1/2]

Bool_t NuageGenerator::Init ( const char *  fileName)
virtual

Definition at line 34 of file NuageGenerator.cxx.

34 {
35 cout<<endl;
36 cout<<"**************************************"<<endl;
37 cout<<"WARNING!!!!!!"<<endl;
38 cout<<"Bug on decays found and solved for GenieGenerator still to be fixed (2\
395.08.2017)" <<endl;
40 cout<<"**************************************"<<endl;
41 cout<<endl;
42 return Init(fileName, 0);
43}
NuageGenerator::Init
virtual Bool_t Init(const char *, int)
Definition NuageGenerator.cxx:45

◆ Init() [2/2]

Bool_t NuageGenerator::Init ( const char *  fileName,
int  firstEvent 
)
virtual

Definition at line 45 of file NuageGenerator.cxx.

45 {
46 cout<<endl;
47 cout<<"**************************************"<<endl;
48 cout<<"WARNING!!!!!!"<<endl;
49 cout<<"Bug on decays found and solved for GenieGenerator still to be fixed (2\
505.08.2017)" <<endl;
51 cout<<"**************************************"<<endl;
52 cout<<endl;
53 fNuOnly = false;
54 LOGF(info, "Opening input file %s", fileName);
55 fInputFile = new TFile(fileName);
56 if (fInputFile->IsZombie()) {
57 LOG(FATAL) << "Error opening the Signal file";
58 }
59 fTree = (TTree *)fInputFile->Get("tree");
60 fNevents = fTree->GetEntries();
61 fn = firstEvent;
62 fTree->SetBranchAddress("Ev",&pxl); // incoming neutrino energy
63 fTree->SetBranchAddress("Nvtx",&fNvtx); //number of vertices
64 fTree->SetBranchAddress("pxv",&pxv);
65 fTree->SetBranchAddress("pyv",&pyv);
66 fTree->SetBranchAddress("pzv",&pzv);
67 fTree->SetBranchAddress("neu",&neu); // incoming neutrino PDG code
68 fTree->SetBranchAddress("cc",&cc); // Is it a CC event?
69 fTree->SetBranchAddress("vtxx",&vtxx); // vertex in SI units
70 fTree->SetBranchAddress("vtxy",&vtxy);
71 fTree->SetBranchAddress("vtxz",&vtxz);
72 fTree->SetBranchAddress("vtxt",&vtxt);
73 fTree->SetBranchAddress("El",&El); // outgoing lepton momentum
74 fTree->SetBranchAddress("pxl",&pxl);
75 fTree->SetBranchAddress("pyl",&pyl);
76 fTree->SetBranchAddress("pzl",&pzl);
77 fTree->SetBranchAddress("pxf",&pxf); // outgoing hadronic momenta
78 fTree->SetBranchAddress("pyf",&pyf);
79 fTree->SetBranchAddress("pzf",&pzf);
80 fTree->SetBranchAddress("nf",&nf); // nr of outgoing hadrons
81 fTree->SetBranchAddress("pdgf",&pdgf); // pdg code of hadron
82 fTree->SetBranchAddress("parent2",&parent2); //parent 2ry vtx tracks
83 fTree->SetBranchAddress("vtxx2",&vtxx2); // 2ry vertex
84 fTree->SetBranchAddress("vtxy2",&vtxy2);
85 fTree->SetBranchAddress("vtxz2",&vtxz2);
86 fTree->SetBranchAddress("vtxt2",&vtxt2);
87 fTree->SetBranchAddress("nf2",&nf2); // nr outgoing particles 2ry vtx
88 fTree->SetBranchAddress("pdgf2",&pdgf2); // pdg 2ry vtx particles
89 fTree->SetBranchAddress("pxf2",&pxf2); // outgoing hadronic momenta
90 fTree->SetBranchAddress("pyf2",&pyf2);
91 fTree->SetBranchAddress("pzf2",&pzf2);
92 fTree->SetBranchAddress("parent3",&parent3); //parent 3rd vtx tracks
93 fTree->SetBranchAddress("vtxx",&vtxx3); // 3rd vertex
94 fTree->SetBranchAddress("vtxy3",&vtxy3);
95 fTree->SetBranchAddress("vtxz3",&vtxz3);
96 fTree->SetBranchAddress("vtxt3",&vtxt3);
97 fTree->SetBranchAddress("nf3",&nf3); // nr outgoing particles 2ry vtx
98 fTree->SetBranchAddress("pdgf3",&pdgf3); // pdg 2ry vtx particles
99 fTree->SetBranchAddress("pxf3",&pxf3); // outgoing hadronic momenta
100 fTree->SetBranchAddress("pyf3",&pyf3);
101 fTree->SetBranchAddress("pzf3",&pzf3);
102 fTree->SetBranchAddress("parent4",&parent4); //parent 4th vtx tracks
103 fTree->SetBranchAddress("vtxx4",&vtxx4); // 4th vertex
104 fTree->SetBranchAddress("vtxy4",&vtxy4);
105 fTree->SetBranchAddress("vtxz4",&vtxz4);
106 fTree->SetBranchAddress("vtxt4",&vtxt4);
107 fTree->SetBranchAddress("nf4",&nf4); // nr outgoing particles 4t vtx
108 fTree->SetBranchAddress("pdgf4",&pdgf4); // pdg 4th vtx particles
109 fTree->SetBranchAddress("pxf4",&pxf4); // outgoing hadronic momenta
110 fTree->SetBranchAddress("pyf4",&pyf4);
111 fTree->SetBranchAddress("pzf4",&pzf4);
112
113 fFirst=kTRUE;
114 return kTRUE;
115}
NuageGenerator::pdgf4
Float_t pdgf4[500]
Definition NuageGenerator.h:61
NuageGenerator::pyf3
Float_t pyf3[500]
Definition NuageGenerator.h:58
NuageGenerator::pzf2
Float_t pzf2[500]
Definition NuageGenerator.h:56
NuageGenerator::pyf
Float_t pyf[500]
Definition NuageGenerator.h:54
NuageGenerator::pdgf3
Float_t pdgf3[500]
Definition NuageGenerator.h:59
NuageGenerator::pzf4
Float_t pzf4[500]
Definition NuageGenerator.h:60
NuageGenerator::pxf3
Float_t pxf3[500]
Definition NuageGenerator.h:58
NuageGenerator::pxf
Float_t pxf[500]
Definition NuageGenerator.h:54
NuageGenerator::pdgf
Float_t pdgf[500]
Definition NuageGenerator.h:55
NuageGenerator::pzf3
Float_t pzf3[500]
Definition NuageGenerator.h:58
NuageGenerator::pyf2
Float_t pyf2[500]
Definition NuageGenerator.h:56
NuageGenerator::pxf2
Float_t pxf2[500]
Definition NuageGenerator.h:56
NuageGenerator::pxf4
Float_t pxf4[500]
Definition NuageGenerator.h:60
NuageGenerator::pzf
Float_t pzf[500]
Definition NuageGenerator.h:54
NuageGenerator::pyf4
Float_t pyf4[500]
Definition NuageGenerator.h:60
NuageGenerator::pdgf2
Float_t pdgf2[500]
Definition NuageGenerator.h:57
MufluxReco.firstEvent
int firstEvent
Definition MufluxReco.py:13

◆ MeanMaterialBudget()

Double_t NuageGenerator::MeanMaterialBudget ( const Double_t *  start,
const Double_t *  end,
Double_t *  mparam 
)
private

Definition at line 118 of file NuageGenerator.cxx.

119{
120 //
121 // Calculate mean material budget and material properties between
122 // the points "start" and "end".
123 //
124 // "mparam" - parameters used for the energy and multiple scattering
125 // corrections:
126 //
127 // mparam[0] - mean density: sum(x_i*rho_i)/sum(x_i) [g/cm3]
128 // mparam[1] - equivalent rad length fraction: sum(x_i/X0_i) [adimensional]
129 // mparam[2] - mean A: sum(x_i*A_i)/sum(x_i) [adimensional]
130 // mparam[3] - mean Z: sum(x_i*Z_i)/sum(x_i) [adimensional]
131 // mparam[4] - length: sum(x_i) [cm]
132 // mparam[5] - Z/A mean: sum(x_i*Z_i/A_i)/sum(x_i) [adimensional]
133 // mparam[6] - number of boundary crosses
134 // mparam[7] - maximum density encountered (g/cm^3)
135 //
136 // Origin: Marian Ivanov, Marian.Ivanov@cern.ch
137 //
138 // Corrections and improvements by
139 // Andrea Dainese, Andrea.Dainese@lnl.infn.it,
140 // Andrei Gheata, Andrei.Gheata@cern.ch
141 //
142
143 mparam[0]=0; mparam[1]=1; mparam[2] =0; mparam[3] =0;
144 mparam[4]=0; mparam[5]=0; mparam[6]=0; mparam[7]=0;
145 //
146 Double_t bparam[6]; // total parameters
147 Double_t lparam[6]; // local parameters
148
149 for (Int_t i=0;i<6;i++) bparam[i]=0;
150
151 if (!gGeoManager) {
152 //AliFatalClass("No TGeo\n");
153 return 0.;
154 }
155 //
156 Double_t length;
157 Double_t dir[3];
158 length = TMath::Sqrt((end[0]-start[0])*(end[0]-start[0])+
159 (end[1]-start[1])*(end[1]-start[1])+
160 (end[2]-start[2])*(end[2]-start[2]));
161 mparam[4]=length;
162 if (length<TGeoShape::Tolerance()) return 0.0;
163 Double_t invlen = 1./length;
164 dir[0] = (end[0]-start[0])*invlen;
165 dir[1] = (end[1]-start[1])*invlen;
166 dir[2] = (end[2]-start[2])*invlen;
167
168 // Initialize start point and direction
169 TGeoNode *currentnode = 0;
170 TGeoNode *startnode = gGeoManager->InitTrack(start, dir);
171 if (!startnode) {
172 //AliErrorClass(Form("start point out of geometry: x %f, y %f, z %f",
173 // start[0],start[1],start[2]));
174 return 0.0;
175 }
176 TGeoMaterial *material = startnode->GetVolume()->GetMedium()->GetMaterial();
177 lparam[0] = material->GetDensity();
178 if (lparam[0] > mparam[7]) mparam[7]=lparam[0];
179 lparam[1] = material->GetRadLen();
180 lparam[2] = material->GetA();
181 lparam[3] = material->GetZ();
182 lparam[4] = length;
183 lparam[5] = lparam[3]/lparam[2];
184 if (material->IsMixture()) {
185 TGeoMixture * mixture = (TGeoMixture*)material;
186 lparam[5] =0;
187 Double_t sum =0;
188 for (Int_t iel=0;iel<mixture->GetNelements();iel++){
189 sum += mixture->GetWmixt()[iel];
190 lparam[5]+= mixture->GetZmixt()[iel]*mixture->GetWmixt()[iel]/mixture->GetAmixt()[iel];
191 }
192 lparam[5]/=sum;
193 }
194
195 // Locate next boundary within length without computing safety.
196 // Propagate either with length (if no boundary found) or just cross boundary
197 gGeoManager->FindNextBoundaryAndStep(length, kFALSE);
198 Double_t step = 0.0; // Step made
199 Double_t snext = gGeoManager->GetStep();
200 // If no boundary within proposed length, return current density
201 if (!gGeoManager->IsOnBoundary()) {
202 mparam[0] = lparam[0];
203 mparam[1] = lparam[4]/lparam[1];
204 mparam[2] = lparam[2];
205 mparam[3] = lparam[3];
206 mparam[4] = lparam[4];
207 return lparam[0];
208 }
209 // Try to cross the boundary and see what is next
210 Int_t nzero = 0;
211 while (length>TGeoShape::Tolerance()) {
212 currentnode = gGeoManager->GetCurrentNode();
213 if (snext<2.*TGeoShape::Tolerance()) nzero++;
214 else nzero = 0;
215 if (nzero>3) {
216 // This means navigation has problems on one boundary
217 // Try to cross by making a small step
218 //AliErrorClass("Cannot cross boundary\n");
219 mparam[0] = bparam[0]/step;
220 mparam[1] = bparam[1];
221 mparam[2] = bparam[2]/step;
222 mparam[3] = bparam[3]/step;
223 mparam[5] = bparam[5]/step;
224 mparam[4] = step;
225 mparam[0] = 0.; // if crash of navigation take mean density 0
226 mparam[1] = 1000000; // and infinite rad length
227 return bparam[0]/step;
228 }
229 mparam[6]+=1.;
230 step += snext;
231 bparam[1] += snext/lparam[1];
232 bparam[2] += snext*lparam[2];
233 bparam[3] += snext*lparam[3];
234 bparam[5] += snext*lparam[5];
235 bparam[0] += snext*lparam[0];
236
237 if (snext>=length) break;
238 if (!currentnode) break;
239 length -= snext;
240 material = currentnode->GetVolume()->GetMedium()->GetMaterial();
241 lparam[0] = material->GetDensity();
242 if (lparam[0] > mparam[7]) mparam[7]=lparam[0];
243 lparam[1] = material->GetRadLen();
244 lparam[2] = material->GetA();
245 lparam[3] = material->GetZ();
246 lparam[5] = lparam[3]/lparam[2];
247 if (material->IsMixture()) {
248 TGeoMixture * mixture = (TGeoMixture*)material;
249 lparam[5]=0;
250 Double_t sum =0;
251 for (Int_t iel=0;iel<mixture->GetNelements();iel++){
252 sum+= mixture->GetWmixt()[iel];
253 lparam[5]+= mixture->GetZmixt()[iel]*mixture->GetWmixt()[iel]/mixture->GetAmixt()[iel];
254 }
255 lparam[5]/=sum;
256 }
257 gGeoManager->FindNextBoundaryAndStep(length, kFALSE);
258 snext = gGeoManager->GetStep();
259 }
260 mparam[0] = bparam[0]/step;
261 mparam[1] = bparam[1];
262 mparam[2] = bparam[2]/step;
263 mparam[3] = bparam[3]/step;
264 mparam[5] = bparam[5]/step;
265 return bparam[0]/step;
266}
ShipAna.i
int i
Definition ShipAna.py:86
add_noise_to_field.step
step
Definition add_noise_to_field.py:69
charm-geometry_config.length
length
Definition charm-geometry_config.py:112
study_GammaConv.material
dict material
Definition study_GammaConv.py:24

◆ NuOnly()

void NuageGenerator::NuOnly ( )
inline

Definition at line 29 of file NuageGenerator.h.

29{fNuOnly = true;}

◆ ReadEvent()

Bool_t NuageGenerator::ReadEvent ( FairPrimaryGenerator *  cpg)

public method ReadEvent

Definition at line 305 of file NuageGenerator.cxx.

306{
307 fFirst = kFALSE;
308 //some start/end positions in z (emulsion to Tracker 1)
309 Double_t start[3]={startX,startY,startZ};
310 //cout << "startX = " << startX << " Y = " << startY << " Z = " << startZ << endl;
311 //cout << "endX = " << endX << " Y = " << endY << " Z = " << endZ << endl;
312
313 Double_t end[3]={endX,endY,endZ};
314 if (fFirst)
315 {
316 Double_t bparam=0.;
317 Double_t mparam[7];
318 bparam=MeanMaterialBudget(start, end, mparam);
319 cout << "Info NuageGenerator: MaterialBudget " << start[2] << " - "<< end[2] << endl;
320 cout << "Info NuageGenerator: MaterialBudget " << bparam << endl;
321 cout << "Info NuageGenerator: MaterialBudget 0 " << mparam[0] << endl;
322 cout << "Info NuageGenerator: MaterialBudget 1 " << mparam[1] << endl;
323 cout << "Info NuageGenerator: MaterialBudget 2 " << mparam[2] << endl;
324 cout << "Info NuageGenerator: MaterialBudget 3 " << mparam[3] << endl;
325 cout << "Info NuageGenerator: MaterialBudget 4 " << mparam[4] << endl;
326 cout << "Info NuageGenerator: MaterialBudget 5 " << mparam[5] << endl;
327 cout << "Info NuageGenerator: MaterialBudget 6 " << mparam[6] << endl;
328 cout << "Info NuageGenerator: MaterialBudget " << mparam[0]*mparam[4] << endl;
329
330 fFirst = kFALSE;
331 }
332 //cout << endl;
333 //cout << "*****************************************************************" << endl;
334
335 if (fn==fNevents)
336 {
337 LOG(WARNING) << "End of input file. Rewind.";
338 }
339 fTree->GetEntry(fn%fNevents);
340 fn++;
341 if (fn%100==0)
342 {
343 cout << "Info NuageGenerator: neutrino event-nr "<< fn << endl;
344 }
345
346 //Mean distance between start[2] and end[2]
347 Double_t zMean = (start[2]+end[2])/2;
348 //Distance from the target (Beam Dump (0,0));
349 Double_t zBD = zMean - ztarget;
350
351 //only accept neutrinos with a ThetaX in the ThetaXMin/Max range and same for ThetaY
352 Double_t ThetaXMax = startX/zBD;
353 Double_t ThetaXMin = endX/zBD;
354 Double_t ThetaYMax = startY/zBD;
355 Double_t ThetaYMin = endY/zBD;
356 // cout << "ThetaXMax = " << ThetaXMax << " ThetaXMin = " << ThetaXMin << endl;
357 //cout << "ThetaYMax = " << ThetaYMax << " ThetaYMin = " << ThetaYMin << endl;
358
359 // Incoming neutrino, get a random px,py
360 //cout << "Info NuageGenerator: neutrino " << neu << "p-in "<< pxv << pyv << pzv << " nf "<< nf << endl;
361 //cout << "Info NuageGenerator: ztarget " << ztarget << endl;
362 Double_t bparam=0.;
363 Double_t mparam[8];
364 Double_t pout[3] ;
365 Double_t txnu=0;
366 Double_t tynu=0;
367 //Does this neutrino fly through material? Otherwise draw another pt..
368 //cout << "Info NuageGenerator Start bparam while loop" << endl;
369 while (bparam<1.e-10)
370 {
371 txnu = gRandom->Uniform(ThetaXMax,ThetaXMin);
372 tynu = gRandom->Uniform(ThetaYMax,ThetaYMin);
373
374 pout[2] = pzv*pzv/(1+txnu*txnu+tynu*tynu);
375
376 if (pout[2]>0.)
377 {
378 pout[2]=TMath::Sqrt(pout[2]);
379 pout[0] = pout[2]*txnu;
380 pout[1] = pout[2]*tynu;
381 //cout << "txnu = " << txnu << " tynu = " << tynu << endl;
382
383 //cout << "Info NuageGenerator: neutrino pxyz " << pout[0] << ", " << pout[1] << ", " << pout[2] << endl;
384
385 start[0]=txnu*(start[2]-ztarget);
386 start[1]=tynu*(start[2]-ztarget);
387 //cout << "Info NuageGenerator: neutrino xyz-start " << start[0] << " - " << start[1] << " - " << start[2] << endl;
388
389 end[0]=txnu*(end[2]-ztarget);
390 end[1]=tynu*(end[2]-ztarget);
391 //cout << "Info NuageGenerator: neutrino xyz-end " << end[0] << " - " << end[1] << " - " << end[2] << endl;
392
393 //get material density between these two points
394 bparam=MeanMaterialBudget(start, end, mparam);
395
396 }
397 else bparam = 1.e-10;
398 }
399
400 //loop over trajectory between start and end to pick an interaction point
401 Double_t prob2int = 0.;
402 Double_t x=0.;
403 Double_t y=0.;
404 Double_t z=0.;
405 //Int_t count=0;
406 //cout << "Info NuageGenerator Start prob2int while loop, bparam= " << bparam << ", " << bparam*1.e8 <<endl;
407 //cout << "Info NuageGenerator What was maximum density, mparam[7]= " << mparam[7] << ", " << mparam[7]*1.e8 <<endl;
408 while (prob2int<gRandom->Uniform(0.,1.))
409 {
410 z=gRandom->Uniform(start[2],end[2]);
411 //x & y are computed so to be along the neutrino trajectory
412 x= txnu*(z-ztarget);
413 y= tynu*(z-ztarget);
414 //cout << "x = " << x << " y = " << y << " z " << z << endl;
415
416 //get local material at this point
417 TGeoNode *node = gGeoManager->FindNode(x,y,z);
418 TGeoMaterial *mat = 0;
419 if (node && !gGeoManager->IsOutside()) mat = node->GetVolume()->GetMaterial();
420 //cout << "Info NuageGenerator: mat " << count << ", " << mat->GetName() << ", " << mat->GetDensity() << endl;
421
422 //density relative to Prob largest density along this trajectory, i.e. use rho(Pt)
423 prob2int= mat->GetDensity()/mparam[7];
424
425 if (prob2int>1.)
426 cout << "***WARNING*** NuageGenerator: prob2int > Maximum density????" << prob2int << " maxrho:" << mparam[7] << " material: " << mat->GetName() << endl;
427 //count+=1;
428 }
429 //cout << "Info NuageGenerator: prob2int " << prob2int << ", " << count << endl;
430
431 //cout <<" Neutrino pdg = " << neu << endl;
432
433 Double_t zrelative=z-ztarget;
434 Double_t tof=TMath::Sqrt(x*x+y*y+zrelative*zrelative)/2.99792458e+6;
435 cpg->AddTrack(neu,pout[0],pout[1],pout[2],x,y,z,-1,false,TMath::Sqrt(pout[0]*pout[0]+pout[1]*pout[1]+pout[2]*pout[2]),tof,mparam[0]*mparam[4]);
436 if (not fNuOnly)
437 {
438 // second, outgoing lepton
439 std::vector<double> pp = Rotate(x,y,zrelative,pxl,pyl,pzl);
440 Int_t oLPdgCode = neu;
441 Int_t nAddTrk=0;
442 if (cc)
443 {
444 oLPdgCode = copysign(TMath::Abs(neu)-1,neu);
445 }
446 if(TMath::Abs(oLPdgCode)!=15)
447 {
448 cpg->AddTrack(oLPdgCode,pp[0],pp[1],pp[2],x,y,z,0,true,TMath::Sqrt(pp[0]*pp[0]+pp[1]*pp[1]+pp[2]*pp[2]),tof,mparam[0]*mparam[4]);
449 nAddTrk++;
450 }
451 //cout << "oLpdgCode " << oLPdgCode << " pz "<< pzl << endl;
452 if(TMath::Abs(oLPdgCode)==15)
453 {
454 if(fExtDecayer==0)
455 {
456 cpg->AddTrack(oLPdgCode,pp[0],pp[1],pp[2],x,y,z,0,false,TMath::Sqrt(pp[0]*pp[0]+pp[1]*pp[1]+pp[2]*pp[2]),tof,mparam[0]*mparam[4]);
457 nAddTrk++;
458 if(TMath::Abs(parent2)==15)
459 {
460 //Coordinate secondary vertex
461 Double_t x2=0., y2=0., z2=0.;
462 x2=x+vtxx2;
463 y2=y+vtxy2;
464 z2=z+vtxz2;
465 Double_t zrelative2 = z2-z;
466 Double_t tof2=TMath::Sqrt(x2*x2+y2*y2+zrelative2*zrelative2)/2.99792458e+6;
467 for(int j=0; j<nf2; j++)
468 {
469 pp = Rotate(x2,y2,zrelative2,pxf2[j],pyf2[j],pzf2[j]);
470 cpg->AddTrack(pdgf2[j],pp[0],pp[1],pp[2],x2,y2,z2,nAddTrk,true,TMath::Sqrt(pp[0]*pp[0]+pp[1]*pp[1]+pp[2]*pp[2]),tof2,mparam[0]*mparam[4]);
471 // cout << "pdgf2 " << pdgf2[j] << " p2 "<< TMath::Sqrt(pp[0]*pp[0]+pp[1]*pp[1]+pp[2]*pp[2]) << endl;
472 }
473 nAddTrk+=nf2;
474 }
475
476 if(TMath::Abs(parent3)==15)
477 {
478 //Coordinate third vertex
479 Double_t x3=0., y3=0., z3=0.;
480 x3=x+vtxx3;
481 y3=y+vtxy3;
482 z3=z+vtxz3;
483 Double_t zrelative3 = z3-z;
484 Double_t tof3=TMath::Sqrt(x3*x3+y3*y3+zrelative3*zrelative3)/2.99792458e+6;
485 for(int j=0; j<nf3; j++)
486 {
487 pp = Rotate(x3,y3,zrelative3,pxf3[j],pyf3[j],pzf3[j]);
488 cpg->AddTrack(pdgf3[j],pp[0],pp[1],pp[2],x3,y3,z3,nAddTrk,true,TMath::Sqrt(pp[0]*pp[0]+pp[1]*pp[1]+pp[2]*pp[2]),tof3,mparam[0]*mparam[4]);
489 }
490 nAddTrk+=nf3;
491 }
492 if(TMath::Abs(parent4)==15)
493 {
494 //Coordinate third vertex
495 Double_t x4=0., y4=0., z4=0.;
496 x4=x+vtxx4;
497 y4=y+vtxy4;
498 z4=z+vtxz4;
499 Double_t zrelative4 = z4-z;
500 Double_t tof4=TMath::Sqrt(x4*x4+y4*y4+zrelative4*zrelative4)/2.99792458e+6;
501 for(int j=0; j<nf4; j++)
502 {
503 pp = Rotate(x4,y4,zrelative4,pxf4[j],pyf4[j],pzf4[j]);
504 cpg->AddTrack(pdgf4[j],pp[0],pp[1],pp[2],x4,y4,z4,nAddTrk,true,TMath::Sqrt(pp[0]*pp[0]+pp[1]*pp[1]+pp[2]*pp[2]),tof4,mparam[0]*mparam[4]);
505 }
506 nAddTrk+=nf4;
507 }
508 }
509 else
510 cpg->AddTrack(oLPdgCode,pp[0],pp[1],pp[2],x,y,z,0,true,TMath::Sqrt(pp[0]*pp[0]+pp[1]*pp[1]+pp[2]*pp[2]),tof,mparam[0]*mparam[4]);
511 }
512
513 // last, all others
514 // cout<< "nf: " << nf<<endl;
515 for(int i=0; i<nf; i++)
516 {
517 pp = Rotate(x,y,zrelative,pxf[i],pyf[i],pzf[i]);
518 if(Abs(pdgf[i])!= 411&&Abs(pdgf[i])!=421&&Abs(pdgf[i])!=431&&Abs(pdgf[i])!=4122)
519 {
520 cpg->AddTrack(pdgf[i],pp[0],pp[1],pp[2],x,y,z,0,true,TMath::Sqrt(pp[0]*pp[0]+pp[1]*pp[1]+pp[2]*pp[2]),tof,mparam[0]*mparam[4]);
521 nAddTrk++;
522 }
523 //cout << "f " << pdgf[i] << " p "<< TMath::Sqrt(pp[0]*pp[0]+pp[1]*pp[1]+pp[2]*pp[2]) << endl;
524 if(Abs(pdgf[i]) == 411 || Abs(pdgf[i])==421 || Abs(pdgf[i])==431||Abs(pdgf[i])==4122)
525 cout << "charm particle " << endl;
526 {
527 if(fExtDecayer==0)
528 {
529 //disabilitated the tracking
530 cpg->AddTrack(pdgf[i],pp[0],pp[1],pp[2],x,y,z,0,false,TMath::Sqrt(pp[0]*pp[0]+pp[1]*pp[1]+pp[2]*pp[2]),tof,mparam[0]*mparam[4]);
531 nAddTrk++;
532 if(TMath::Abs(parent2)==pdgf[i])
533 {
534 //Coordinate secondary vertex
535 Double_t x2=0., y2=0., z2=0.;
536 x2=x+vtxx2;
537 y2=y+vtxy2;
538 z2=z+vtxz2;
539 Double_t zrelative2 = z2-z;
540 Double_t tof2=TMath::Sqrt(x2*x2+y2*y2+zrelative2*zrelative2)/2.99792458e+6;
541 for(int j=0; j<nf2; j++)
542 {
543 pp = Rotate(x2,y2,zrelative2,pxf2[j],pyf2[j],pzf2[j]);
544 cpg->AddTrack(pdgf2[j],pp[0],pp[1],pp[2],x2,y2,z2,nAddTrk,true,TMath::Sqrt(pp[0]*pp[0]+pp[1]*pp[1]+pp[2]*pp[2]),tof2,mparam[0]*mparam[4]);
545 // cout << "pdgf2 " << pdgf2[j] << " p2 "<< TMath::Sqrt(pp[0]*pp[0]+pp[1]*pp[1]+pp[2]*pp[2]) << endl;
546 }
547 nAddTrk+=nf2;
548 }
549
550 if(TMath::Abs(parent3)==pdgf[i])
551 {
552 //Coordinate third vertex
553 Double_t x3=0., y3=0., z3=0.;
554 x3=x+vtxx3;
555 y3=y+vtxy3;
556 z3=z+vtxz3;
557 Double_t zrelative3 = z3-z;
558 Double_t tof3=TMath::Sqrt(x3*x3+y3*y3+zrelative3*zrelative3)/2.99792458e+6;
559 for(int j=0; j<nf3; j++)
560 {
561 pp = Rotate(x3,y3,zrelative3,pxf3[j],pyf3[j],pzf3[j]);
562 cpg->AddTrack(pdgf3[j],pp[0],pp[1],pp[2],x3,y3,z3,nAddTrk,true,TMath::Sqrt(pp[0]*pp[0]+pp[1]*pp[1]+pp[2]*pp[2]),tof3,mparam[0]*mparam[4]);
563 }
564 nAddTrk+=nf3;
565 }
566 if(TMath::Abs(parent4)==pdgf[i])
567 {
568 //Coordinate third vertex
569 Double_t x4=0., y4=0., z4=0.;
570 x4=x+vtxx4;
571 y4=y+vtxy4;
572 z4=z+vtxz4;
573 Double_t zrelative4 = z4-z;
574 Double_t tof4=TMath::Sqrt(x4*x4+y4*y4+zrelative4*zrelative4)/2.99792458e+6;
575 for(int j=0; j<nf4; j++)
576 {
577 pp = Rotate(x4,y4,zrelative4,pxf4[j],pyf4[j],pzf4[j]);
578 cpg->AddTrack(pdgf4[j],pp[0],pp[1],pp[2],x4,y4,z4,nAddTrk,true,TMath::Sqrt(pp[0]*pp[0]+pp[1]*pp[1]+pp[2]*pp[2]),tof4,mparam[0]*mparam[4]);
579 }
580 nAddTrk+=nf4;
581 }
582 }
583 else
584 cpg->AddTrack(pdgf[i],pp[0],pp[1],pp[2],x,y,z,0,true,TMath::Sqrt(pp[0]*pp[0]+pp[1]*pp[1]+pp[2]*pp[2]),tof,mparam[0]*mparam[4]);
585 }
586
587
588 }
589 }
590
591 return kTRUE;
592}
NuageGenerator::MeanMaterialBudget
Double_t MeanMaterialBudget(const Double_t *start, const Double_t *end, Double_t *mparam)
Definition NuageGenerator.cxx:118
NuageGenerator::Rotate
std::vector< double > Rotate(Double_t x, Double_t y, Double_t z, Double_t px, Double_t py, Double_t pz)
Definition NuageGenerator.cxx:270
charm-geometry_config.z
z
Definition charm-geometry_config.py:7
conditionsDatabase.snd_conDB.j
int j
Definition snd_conDB.py:197
geometry_config.x
x
Definition geometry_config.py:180
geometry_config.y
y
Definition geometry_config.py:179
geometry_config.z4
int z4
Definition geometry_config.py:81

◆ Rotate()

std::vector< double > NuageGenerator::Rotate ( Double_t  x,
Double_t  y,
Double_t  z,
Double_t  px,
Double_t  py,
Double_t  pz 
)
private

Definition at line 270 of file NuageGenerator.cxx.

271{
272 //rotate vector px,py,pz to point at x,y,z at origin.
273 Double_t theta=atan(sqrt(x*x+y*y)/z);
274 Double_t c=cos(theta);
275 Double_t s=sin(theta);
276 //rotate around y-axis
277 Double_t px1=c*px+s*pz;
278 Double_t pzr=-s*px+c*pz;
279 Double_t phi=atan2(y,x);
280 c=cos(phi);
281 s=sin(phi);
282 //rotate around z-axis
283 Double_t pxr=c*px1-s*py;
284 Double_t pyr=s*px1+c*py;
285 std::vector<double> pout;
286 pout.push_back(pxr);
287 pout.push_back(pyr);
288 pout.push_back(pzr);
289 //cout << "Info NuageGenerator: rotated" << pout[0] << " " << pout[1] << " " << pout[2] << " " << x << " " << y << " " << z <<endl;
290 return pout;
291}
hnl.c
c
Definition hnl.py:100
study_GammaConv.s
s
Definition study_GammaConv.py:22

◆ SetPositions()

void NuageGenerator::SetPositions ( Double_t  zTa,
Double_t  zS = -3400.,
Double_t  zE = 2650.,
Double_t  xS = 0,
Double_t  xE = 0,
Double_t  yS = 0,
Double_t  yE = 0 
)
inline

Definition at line 30 of file NuageGenerator.h.

30 {
31 ztarget = zTa;
32 startZ = zS;
33 endZ = zE;
34 startX = xS;
35 endX = xE;
36 startY = yS;
37 endY = yE;
38 }

Member Data Documentation

◆ boxs

std::vector<TVector3> NuageGenerator::boxs
protected

Definition at line 63 of file NuageGenerator.h.

◆ cc

Int_t NuageGenerator::cc
protected

Definition at line 64 of file NuageGenerator.h.

◆ dVecs

std::vector<TVector3> NuageGenerator::dVecs
protected

Definition at line 62 of file NuageGenerator.h.

◆ El

Float_t NuageGenerator::El
protected

Definition at line 50 of file NuageGenerator.h.

◆ endX

Double_t NuageGenerator::endX
protected

Definition at line 49 of file NuageGenerator.h.

◆ endY

Double_t NuageGenerator::endY
protected

Definition at line 49 of file NuageGenerator.h.

◆ endZ

Double_t NuageGenerator::endZ
protected

Definition at line 49 of file NuageGenerator.h.

◆ Ev

Float_t NuageGenerator::Ev
protected

Definition at line 50 of file NuageGenerator.h.

◆ fEntrA

Double_t NuageGenerator::fEntrA
protected

Definition at line 76 of file NuageGenerator.h.

◆ fEntrB

Double_t NuageGenerator::fEntrB
protected

Definition at line 76 of file NuageGenerator.h.

◆ fEntrDz_inner

Double_t NuageGenerator::fEntrDz_inner
protected

Definition at line 76 of file NuageGenerator.h.

◆ fEntrDz_outer

Double_t NuageGenerator::fEntrDz_outer
protected

Definition at line 76 of file NuageGenerator.h.

◆ fEntrZ_inner

Double_t NuageGenerator::fEntrZ_inner
protected

Definition at line 76 of file NuageGenerator.h.

◆ fEntrZ_outer

Double_t NuageGenerator::fEntrZ_outer
protected

Definition at line 76 of file NuageGenerator.h.

◆ fExtDecayer

Bool_t NuageGenerator::fExtDecayer
protected

Definition at line 48 of file NuageGenerator.h.

◆ fFirst

bool NuageGenerator::fFirst
protected

Definition at line 73 of file NuageGenerator.h.

◆ fInputFile

TFile* NuageGenerator::fInputFile
protected

don't make it persistent, magic ROOT command

Definition at line 69 of file NuageGenerator.h.

◆ fL1z

Double_t NuageGenerator::fL1z
protected

Definition at line 76 of file NuageGenerator.h.

◆ fLogger

FairLogger* NuageGenerator::fLogger
protected

Definition at line 68 of file NuageGenerator.h.

◆ fn

Int_t NuageGenerator::fn
protected

Definition at line 72 of file NuageGenerator.h.

◆ fNevents

int NuageGenerator::fNevents
protected

Definition at line 71 of file NuageGenerator.h.

◆ fNuOnly

bool NuageGenerator::fNuOnly
protected

Definition at line 73 of file NuageGenerator.h.

◆ fNvtx

Int_t NuageGenerator::fNvtx
protected

Definition at line 67 of file NuageGenerator.h.

◆ fScintDz

Double_t NuageGenerator::fScintDz
protected

Definition at line 76 of file NuageGenerator.h.

◆ fTree

TTree* NuageGenerator::fTree
protected

Definition at line 70 of file NuageGenerator.h.

◆ fXnu

Double_t NuageGenerator::fXnu
protected

Definition at line 75 of file NuageGenerator.h.

◆ fYnu

Double_t NuageGenerator::fYnu
protected

Definition at line 75 of file NuageGenerator.h.

◆ fznu

Double_t NuageGenerator::fznu
protected

Definition at line 75 of file NuageGenerator.h.

◆ Lvessel

Double_t NuageGenerator::Lvessel
protected

Definition at line 49 of file NuageGenerator.h.

◆ neu

Int_t NuageGenerator::neu
protected

Definition at line 65 of file NuageGenerator.h.

◆ nf

Int_t NuageGenerator::nf
protected

Definition at line 65 of file NuageGenerator.h.

◆ nf2

Int_t NuageGenerator::nf2
protected

Definition at line 65 of file NuageGenerator.h.

◆ nf3

Int_t NuageGenerator::nf3
protected

Definition at line 65 of file NuageGenerator.h.

◆ nf4

Int_t NuageGenerator::nf4
protected

Definition at line 65 of file NuageGenerator.h.

◆ parent2

Int_t NuageGenerator::parent2
protected

Definition at line 66 of file NuageGenerator.h.

◆ parent3

Int_t NuageGenerator::parent3
protected

Definition at line 66 of file NuageGenerator.h.

◆ parent4

Int_t NuageGenerator::parent4
protected

Definition at line 66 of file NuageGenerator.h.

◆ pdgf

Float_t NuageGenerator::pdgf[500]
protected

Definition at line 55 of file NuageGenerator.h.

◆ pdgf2

Float_t NuageGenerator::pdgf2[500]
protected

Definition at line 57 of file NuageGenerator.h.

◆ pdgf3

Float_t NuageGenerator::pdgf3[500]
protected

Definition at line 59 of file NuageGenerator.h.

◆ pdgf4

Float_t NuageGenerator::pdgf4[500]
protected

Definition at line 61 of file NuageGenerator.h.

◆ pxf

Float_t NuageGenerator::pxf[500]
protected

Definition at line 54 of file NuageGenerator.h.

◆ pxf2

Float_t NuageGenerator::pxf2[500]
protected

Definition at line 56 of file NuageGenerator.h.

◆ pxf3

Float_t NuageGenerator::pxf3[500]
protected

Definition at line 58 of file NuageGenerator.h.

◆ pxf4

Float_t NuageGenerator::pxf4[500]
protected

Definition at line 60 of file NuageGenerator.h.

◆ pxl

Float_t NuageGenerator::pxl
protected

Definition at line 50 of file NuageGenerator.h.

◆ pxv

Float_t NuageGenerator::pxv
protected

Definition at line 50 of file NuageGenerator.h.

◆ pyf

Float_t NuageGenerator::pyf[500]
protected

Definition at line 54 of file NuageGenerator.h.

◆ pyf2

Float_t NuageGenerator::pyf2[500]
protected

Definition at line 56 of file NuageGenerator.h.

◆ pyf3

Float_t NuageGenerator::pyf3[500]
protected

Definition at line 58 of file NuageGenerator.h.

◆ pyf4

Float_t NuageGenerator::pyf4[500]
protected

Definition at line 60 of file NuageGenerator.h.

◆ pyl

Float_t NuageGenerator::pyl
protected

Definition at line 50 of file NuageGenerator.h.

◆ pyv

Float_t NuageGenerator::pyv
protected

Definition at line 50 of file NuageGenerator.h.

◆ pzf

Float_t NuageGenerator::pzf[500]
protected

Definition at line 54 of file NuageGenerator.h.

◆ pzf2

Float_t NuageGenerator::pzf2[500]
protected

Definition at line 56 of file NuageGenerator.h.

◆ pzf3

Float_t NuageGenerator::pzf3[500]
protected

Definition at line 58 of file NuageGenerator.h.

◆ pzf4

Float_t NuageGenerator::pzf4[500]
protected

Definition at line 60 of file NuageGenerator.h.

◆ pzl

Float_t NuageGenerator::pzl
protected

Definition at line 50 of file NuageGenerator.h.

◆ pzv

Float_t NuageGenerator::pzv
protected

Definition at line 50 of file NuageGenerator.h.

◆ startX

Double_t NuageGenerator::startX
protected

Definition at line 49 of file NuageGenerator.h.

◆ startY

Double_t NuageGenerator::startY
protected

Definition at line 49 of file NuageGenerator.h.

◆ startZ

Double_t NuageGenerator::startZ
protected

Definition at line 49 of file NuageGenerator.h.

◆ vtxt

Float_t NuageGenerator::vtxt
protected

Definition at line 50 of file NuageGenerator.h.

◆ vtxt2

Float_t NuageGenerator::vtxt2
protected

Definition at line 51 of file NuageGenerator.h.

◆ vtxt3

Float_t NuageGenerator::vtxt3
protected

Definition at line 52 of file NuageGenerator.h.

◆ vtxt4

Float_t NuageGenerator::vtxt4
protected

Definition at line 53 of file NuageGenerator.h.

◆ vtxx

Float_t NuageGenerator::vtxx
protected

Definition at line 50 of file NuageGenerator.h.

◆ vtxx2

Float_t NuageGenerator::vtxx2
protected

Definition at line 51 of file NuageGenerator.h.

◆ vtxx3

Float_t NuageGenerator::vtxx3
protected

Definition at line 53 of file NuageGenerator.h.

◆ vtxx4

Float_t NuageGenerator::vtxx4
protected

Definition at line 52 of file NuageGenerator.h.

◆ vtxy

Float_t NuageGenerator::vtxy
protected

Definition at line 50 of file NuageGenerator.h.

◆ vtxy2

Float_t NuageGenerator::vtxy2
protected

Definition at line 51 of file NuageGenerator.h.

◆ vtxy3

Float_t NuageGenerator::vtxy3
protected

Definition at line 52 of file NuageGenerator.h.

◆ vtxy4

Float_t NuageGenerator::vtxy4
protected

Definition at line 53 of file NuageGenerator.h.

◆ vtxz

Float_t NuageGenerator::vtxz
protected

Definition at line 50 of file NuageGenerator.h.

◆ vtxz2

Float_t NuageGenerator::vtxz2
protected

Definition at line 51 of file NuageGenerator.h.

◆ vtxz3

Float_t NuageGenerator::vtxz3
protected

Definition at line 52 of file NuageGenerator.h.

◆ vtxz4

Float_t NuageGenerator::vtxz4
protected

Definition at line 53 of file NuageGenerator.h.

◆ Xvessel

Double_t NuageGenerator::Xvessel
protected

Definition at line 49 of file NuageGenerator.h.

◆ Yvessel

Double_t NuageGenerator::Yvessel
protected

Definition at line 49 of file NuageGenerator.h.

◆ ztarget

Double_t NuageGenerator::ztarget
protected

Definition at line 49 of file NuageGenerator.h.

The documentation for this class was generated from the following files: