SND@LHC Software
Loading...
Searching...
No Matches
NuageGenerator.cxx
Go to the documentation of this file.
1#include <math.h>
2#include "TROOT.h"
3#include "TMath.h"
4#include "TFile.h"
5#include "TRandom.h"
6#include "FairPrimaryGenerator.h"
7#include "NuageGenerator.h"
8#include "TGeoVolume.h"
9#include "TGeoNode.h"
10#include "TGeoManager.h"
11#include "TGeoEltu.h"
12#include "TGeoCompositeShape.h"
13#include "TVirtualMC.h"
14
15using std::cout;
16using std::endl;
17using namespace TMath;
18// read events from ntuples produced with Nuage
19// Nuage momentum GeV
20// Vertex in SI units, assume this means m
21// important to read back number of events to give to FairRoot
22
23// ----- Default constructor -------------------------------------------
24NuageGenerator::NuageGenerator() {
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}
32// -------------------------------------------------------------------------
33// ----- Default constructor -------------------------------------------
34Bool_t NuageGenerator::Init(const char* fileName) {
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}
44// ----- Default constructor -------------------------------------------
45Bool_t NuageGenerator::Init(const char* fileName, const int firstEvent) {
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}
116// -------------------------------------------------------------------------
117
118Double_t NuageGenerator::MeanMaterialBudget(const Double_t *start, const Double_t *end, Double_t *mparam)
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}
267
268//----------------------------------------------------------------------------------------
269
270std::vector<double> NuageGenerator::Rotate(Double_t x, Double_t y, Double_t z, Double_t px, Double_t py, Double_t pz)
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}
292
293
294// ----- Destructor ----------------------------------------------------
295NuageGenerator::~NuageGenerator()
296{
297 fInputFile->Close();
298 fInputFile->Delete();
299 delete fInputFile;
300}
301// -------------------------------------------------------------------------
302
303
304// ----- Passing the event ---------------------------------------------
305Bool_t NuageGenerator::ReadEvent(FairPrimaryGenerator* cpg)
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}
593
594
595// -------------------------------------------------------------------------
596Int_t NuageGenerator::GetNevents()
597{
598 return fNevents;
599}
600
601
602ClassImp(NuageGenerator)
NuageGenerator::pdgf4
Float_t pdgf4[500]
Definition NuageGenerator.h:61
NuageGenerator::pyf3
Float_t pyf3[500]
Definition NuageGenerator.h:58
NuageGenerator::Init
virtual Bool_t Init(const char *, int)
Definition NuageGenerator.cxx:45
NuageGenerator::pzf2
Float_t pzf2[500]
Definition NuageGenerator.h:56
NuageGenerator::pyf
Float_t pyf[500]
Definition NuageGenerator.h:54
NuageGenerator::ReadEvent
Bool_t ReadEvent(FairPrimaryGenerator *)
Definition NuageGenerator.cxx:305
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::fInputFile
TFile * fInputFile
don't make it persistent, magic ROOT command
Definition NuageGenerator.h:69
NuageGenerator::pxf
Float_t pxf[500]
Definition NuageGenerator.h:54
NuageGenerator::pdgf
Float_t pdgf[500]
Definition NuageGenerator.h:55
NuageGenerator::~NuageGenerator
virtual ~NuageGenerator()
Definition NuageGenerator.cxx:295
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::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
NuageGenerator::pyf4
Float_t pyf4[500]
Definition NuageGenerator.h:60
NuageGenerator::pdgf2
Float_t pdgf2[500]
Definition NuageGenerator.h:57
ClassImp
ClassImp(ecalContFact) ecalContFact
Definition ecalContFact.cxx:8