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

#include <MuDISGenerator.h>

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

Public Member Functions

 MuDISGenerator ()
 
virtual ~MuDISGenerator ()
 
Bool_t ReadEvent (FairPrimaryGenerator *)
 
virtual Bool_t Init (const char *, int)
 
virtual Bool_t Init (const char *)
 
Int_t GetNevents ()
 
void SetPositions (Double_t zTa, Double_t zS=-3400., Double_t zE=2650.)
 

Protected Member Functions

 ClassDef (MuDISGenerator, 1)
 

Protected Attributes

Double_t startZ
 
Double_t endZ
 
TClonesArray * iMuon
 
TClonesArray * dPart
 
FairLogger * fLogger
 
TFile * fInputFile
 don't make it persistent, magic ROOT command
 
TTree * fTree
 
int fNevents
 
int fn
 
bool fFirst
 

Private Member Functions

Double_t MeanMaterialBudget (const Double_t *start, const Double_t *end, Double_t *mparam)
 

Detailed Description

Definition at line 15 of file MuDISGenerator.h.

Constructor & Destructor Documentation

◆ MuDISGenerator()

MuDISGenerator::MuDISGenerator ( )

default constructor

Definition at line 26 of file MuDISGenerator.cxx.

26{}

◆ ~MuDISGenerator()

MuDISGenerator::~MuDISGenerator ( )
virtual

destructor

Definition at line 207 of file MuDISGenerator.cxx.

208{
209 fInputFile->Close();
210 fInputFile->Delete();
211 delete fInputFile;
212}
MuDISGenerator::fInputFile
TFile * fInputFile
don't make it persistent, magic ROOT command
Definition MuDISGenerator.h:45

Member Function Documentation

◆ ClassDef()

MuDISGenerator::ClassDef ( MuDISGenerator  ,
1   
)
protected

◆ GetNevents()

Int_t MuDISGenerator::GetNevents ( )

Definition at line 292 of file MuDISGenerator.cxx.

293{
294 return fNevents;
295}

◆ Init() [1/2]

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

Definition at line 29 of file MuDISGenerator.cxx.

29 {
30 return Init(fileName, 0);
31}
MuDISGenerator::Init
virtual Bool_t Init(const char *, int)
Definition MuDISGenerator.cxx:33

◆ Init() [2/2]

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

Definition at line 33 of file MuDISGenerator.cxx.

33 {
34 LOG(INFO) << "Opening input file"<< fileName;
35
36 iMuon = 0;
37 dPart = 0;
38 if (0 == strncmp("/eos",fileName,4) ) {
39 TString tmp = gSystem->Getenv("EOSSHIP");
40 tmp+=fileName;
41 fInputFile = TFile::Open(tmp);
42 LOG(INFO) << "Open external file on eos:"<<tmp.Data();
43 }else{
44 fInputFile = new TFile(fileName);
45 }
46 if (fInputFile->IsZombie() or !fInputFile) {
47 LOG(FATAL) << "Error opening input file";
48 return kFALSE; }
49 fTree = (TTree *)fInputFile->Get("DIS");
50 fNevents = fTree->GetEntries();
51 fn = firstEvent;
52 fTree->SetBranchAddress("InMuon",&iMuon); // incoming muon
53 fTree->SetBranchAddress("Particles",&dPart);
54 return kTRUE;
55}
MuDISGenerator::iMuon
TClonesArray * iMuon
Definition MuDISGenerator.h:42
MuDISGenerator::dPart
TClonesArray * dPart
Definition MuDISGenerator.h:43
MufluxReco.firstEvent
int firstEvent
Definition MufluxReco.py:13

◆ MeanMaterialBudget()

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

Definition at line 56 of file MuDISGenerator.cxx.

57{
58 //
59 // Calculate mean material budget and material properties between
60 // the points "start" and "end".
61 //
62 // "mparam" - parameters used for the energy and multiple scattering
63 // corrections:
64 //
65 // mparam[0] - mean density: sum(x_i*rho_i)/sum(x_i) [g/cm3]
66 // mparam[1] - equivalent rad length fraction: sum(x_i/X0_i) [adimensional]
67 // mparam[2] - mean A: sum(x_i*A_i)/sum(x_i) [adimensional]
68 // mparam[3] - mean Z: sum(x_i*Z_i)/sum(x_i) [adimensional]
69 // mparam[4] - length: sum(x_i) [cm]
70 // mparam[5] - Z/A mean: sum(x_i*Z_i/A_i)/sum(x_i) [adimensional]
71 // mparam[6] - number of boundary crosses
72 // mparam[7] - maximum density encountered (g/cm^3)
73 //
74 // Origin: Marian Ivanov, Marian.Ivanov@cern.ch
75 //
76 // Corrections and improvements by
77 // Andrea Dainese, Andrea.Dainese@lnl.infn.it,
78 // Andrei Gheata, Andrei.Gheata@cern.ch
79 //
80
81 mparam[0]=0; mparam[1]=1; mparam[2] =0; mparam[3] =0;
82 mparam[4]=0; mparam[5]=0; mparam[6]=0; mparam[7]=0;
83 //
84 Double_t bparam[6]; // total parameters
85 Double_t lparam[6]; // local parameters
86
87 for (Int_t i=0;i<6;i++) bparam[i]=0;
88
89 if (!gGeoManager) {
90 //AliFatalClass("No TGeo\n");
91 return 0.;
92 }
93 //
94 Double_t length;
95 Double_t dir[3];
96 length = TMath::Sqrt((end[0]-start[0])*(end[0]-start[0])+
97 (end[1]-start[1])*(end[1]-start[1])+
98 (end[2]-start[2])*(end[2]-start[2]));
99 mparam[4]=length;
100 if (length<TGeoShape::Tolerance()) return 0.0;
101 Double_t invlen = 1./length;
102 dir[0] = (end[0]-start[0])*invlen;
103 dir[1] = (end[1]-start[1])*invlen;
104 dir[2] = (end[2]-start[2])*invlen;
105
106 // Initialize start point and direction
107 TGeoNode *currentnode = 0;
108 TGeoNode *startnode = gGeoManager->InitTrack(start, dir);
109 if (!startnode) {
110 //AliErrorClass(Form("start point out of geometry: x %f, y %f, z %f",
111 // start[0],start[1],start[2]));
112 return 0.0;
113 }
114 TGeoMaterial *material = startnode->GetVolume()->GetMedium()->GetMaterial();
115 lparam[0] = material->GetDensity();
116 if (lparam[0] > mparam[7]) mparam[7]=lparam[0];
117 lparam[1] = material->GetRadLen();
118 lparam[2] = material->GetA();
119 lparam[3] = material->GetZ();
120 lparam[4] = length;
121 lparam[5] = lparam[3]/lparam[2];
122 if (material->IsMixture()) {
123 TGeoMixture * mixture = (TGeoMixture*)material;
124 lparam[5] =0;
125 Double_t sum =0;
126 for (Int_t iel=0;iel<mixture->GetNelements();iel++){
127 sum += mixture->GetWmixt()[iel];
128 lparam[5]+= mixture->GetZmixt()[iel]*mixture->GetWmixt()[iel]/mixture->GetAmixt()[iel];
129 }
130 lparam[5]/=sum;
131 }
132
133 // Locate next boundary within length without computing safety.
134 // Propagate either with length (if no boundary found) or just cross boundary
135 gGeoManager->FindNextBoundaryAndStep(length, kFALSE);
136 Double_t step = 0.0; // Step made
137 Double_t snext = gGeoManager->GetStep();
138 // If no boundary within proposed length, return current density
139 if (!gGeoManager->IsOnBoundary()) {
140 mparam[0] = lparam[0];
141 mparam[1] = lparam[4]/lparam[1];
142 mparam[2] = lparam[2];
143 mparam[3] = lparam[3];
144 mparam[4] = lparam[4];
145 return lparam[0];
146 }
147 // Try to cross the boundary and see what is next
148 Int_t nzero = 0;
149 while (length>TGeoShape::Tolerance()) {
150 currentnode = gGeoManager->GetCurrentNode();
151 if (snext<2.*TGeoShape::Tolerance()) nzero++;
152 else nzero = 0;
153 if (nzero>3) {
154 // This means navigation has problems on one boundary
155 // Try to cross by making a small step
156 //AliErrorClass("Cannot cross boundary\n");
157 mparam[0] = bparam[0]/step;
158 mparam[1] = bparam[1];
159 mparam[2] = bparam[2]/step;
160 mparam[3] = bparam[3]/step;
161 mparam[5] = bparam[5]/step;
162 mparam[4] = step;
163 mparam[0] = 0.; // if crash of navigation take mean density 0
164 mparam[1] = 1000000; // and infinite rad length
165 return bparam[0]/step;
166 }
167 mparam[6]+=1.;
168 step += snext;
169 bparam[1] += snext/lparam[1];
170 bparam[2] += snext*lparam[2];
171 bparam[3] += snext*lparam[3];
172 bparam[5] += snext*lparam[5];
173 bparam[0] += snext*lparam[0];
174
175 if (snext>=length) break;
176 if (!currentnode) break;
177 length -= snext;
178 material = currentnode->GetVolume()->GetMedium()->GetMaterial();
179 lparam[0] = material->GetDensity();
180 if (lparam[0] > mparam[7]) mparam[7]=lparam[0];
181 lparam[1] = material->GetRadLen();
182 lparam[2] = material->GetA();
183 lparam[3] = material->GetZ();
184 lparam[5] = lparam[3]/lparam[2];
185 if (material->IsMixture()) {
186 TGeoMixture * mixture = (TGeoMixture*)material;
187 lparam[5]=0;
188 Double_t sum =0;
189 for (Int_t iel=0;iel<mixture->GetNelements();iel++){
190 sum+= mixture->GetWmixt()[iel];
191 lparam[5]+= mixture->GetZmixt()[iel]*mixture->GetWmixt()[iel]/mixture->GetAmixt()[iel];
192 }
193 lparam[5]/=sum;
194 }
195 gGeoManager->FindNextBoundaryAndStep(length, kFALSE);
196 snext = gGeoManager->GetStep();
197 }
198 mparam[0] = bparam[0]/step;
199 mparam[1] = bparam[1];
200 mparam[2] = bparam[2]/step;
201 mparam[3] = bparam[3]/step;
202 mparam[5] = bparam[5]/step;
203 return bparam[0]/step;
204}
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

◆ ReadEvent()

Bool_t MuDISGenerator::ReadEvent ( FairPrimaryGenerator *  cpg)

public method ReadEvent

Definition at line 214 of file MuDISGenerator.cxx.

215{
216 if (fn==fNevents) {LOG(WARNING) << "End of input file. Rewind.";}
217 fTree->GetEntry(fn%fNevents);
218 fn++;
219 if (fn%10000==0) {
220 LOG(INFO) << " heartbeat event-nr "<< fn;
221 }
222 int nf = dPart->GetEntries();
223 LOG(DEBUG) << "*********************************************************";
224 LOG(DEBUG) << " " << iMuon->GetEntries()<< " "<< iMuon->AddrAt(0) << " nf "<< nf << " fn=" << fn ;
225
226 //some start/end positions in z (f.i. emulsion to Tracker 1)
227 Double_t start[3]={0.,0.,startZ};
228 Double_t end[3]={0.,0.,endZ};
229
230 TParticle* mu = dynamic_cast<TParticle*>(iMuon->AddrAt(0));
231 LOG(DEBUG) << " in muon " << mu->GetPdgCode()<< endl;
232 Double_t x = mu->Vx();
233 Double_t y = mu->Vy();
234 Double_t z = mu->Vz();
235 Double_t w = mu->GetWeight();
236// calculate start/end positions along this muon, and amout of material in between
237 Double_t txmu=mu->Px()/mu->Pz();
238 Double_t tymu=mu->Py()/mu->Pz();
239 start[0]=x-(z-start[2])*txmu;
240 start[1]=y-(z-start[2])*tymu;
241 end[0] =x-(z-end[2])*txmu;
242 end[1] =y-(z-end[2])*tymu;
243 LOG(DEBUG) << " mu xyz position " << x << ", " << y << ", " << z << endl;
244 LOG(DEBUG) << " mu pxyz " << mu->Px() << ", " << mu->Py() << ", " << mu->Pz() ;
245 LOG(DEBUG) << " mu weight " << w ;
246 LOG(DEBUG) << " start position " << start[0] << ", " << start[1] << ", " << start[2] ;
247 LOG(DEBUG) << " end position " << end[0] << ", " << end[1] << ", " << end[2] ;
248 Double_t bparam;
249 Double_t mparam[8];
250 bparam=MeanMaterialBudget(start, end, mparam);
251 //loop over trajectory between start and end to pick an interaction point
252 Double_t prob2int = 0.;
253 Double_t xmu;
254 Double_t ymu;
255 Double_t zmu;
256 Int_t count=0;
257 LOG(DEBUG) << " Start prob2int while loop, bparam= " << bparam << ", " << bparam*1.e8 ;
258 LOG(DEBUG) << " What was maximum density, mparam[7]= " << mparam[7] << ", " << mparam[7]*1.e8 ;
259 while (prob2int<gRandom->Uniform(0.,1.)) {
260 zmu=gRandom->Uniform(start[2],end[2]);
261 xmu=x-(z-zmu)*txmu;
262 ymu=y-(z-zmu)*tymu;
263 //get local material at this point
264 TGeoNode *node = gGeoManager->FindNode(xmu,ymu,zmu);
265 TGeoMaterial *mat = 0;
266 if (node && !gGeoManager->IsOutside()) mat = node->GetVolume()->GetMaterial();
267 LOG(DEBUG) << "mat " << count << ", " << mat->GetName() << ", " << mat->GetDensity();
268 //density relative to Prob largest density along this trajectory, i.e. use rho(Pt)
269 prob2int= mat->GetDensity()/mparam[7];
270 if (prob2int>1.) LOG(WARNING) << "***WARNING*** MuDISGenerator: prob2int > Maximum density????" << prob2int << " maxrho:" << mparam[7] << " material: " << mat->GetName();
271 count+=1;
272 }
273 LOG(DEBUG) << "prob2int " << prob2int << ", " << count ;
274
275 LOG(DEBUG) << "put position " << xmu << ", " << ymu << ", " << zmu;
276 //modify weight, by multiplying with average densiy along track??
277 cpg->AddTrack(mu->GetPdgCode(),mu->Px(),mu->Py(),mu->Pz(),xmu,ymu,zmu,-1,false,mu->Energy(),mu->T(),w);
278 // in order to have a simulation of possible veto detector hits, let Geant4 follow the muon backward
279 // due to dE/dx, as soon as the muon travers thick material, this approximation will become bad.
280 // a proper implementation however would be to have this better integrated in Geant4, follow the muon, call DIS event, continue
281 cpg->AddTrack(mu->GetPdgCode(),-(mu->Px()),-(mu->Py()),-(mu->Pz()),xmu,ymu,zmu,0,true,mu->Energy(),mu->T(),w);
282// outgoing particles, [did,dpx,dpy,dpz,E], put density along trajectory as weight, g/cm^2
283 w=mparam[0]*mparam[4];
284 for(int i=0; i<nf; i++)
285 {
286 TParticle* Part = dynamic_cast<TParticle*>(dPart->AddrAt(i));
287 cpg->AddTrack(Part->GetPdgCode(),Part->Px(),Part->Py(),Part->Pz(),xmu,ymu,zmu,0,true,Part->Energy(),mu->T(),w);
288 }
289 return kTRUE;
290}
MuDISGenerator::MeanMaterialBudget
Double_t MeanMaterialBudget(const Double_t *start, const Double_t *end, Double_t *mparam)
Definition MuDISGenerator.cxx:56
charm-geometry_config.z
z
Definition charm-geometry_config.py:7
conditionsDatabase.getconditions_condDB.count
int count
Definition getconditions_condDB.py:39
geometry_config.x
x
Definition geometry_config.py:180
geometry_config.y
y
Definition geometry_config.py:179

◆ SetPositions()

void MuDISGenerator::SetPositions ( Double_t  zTa,
Double_t  zS = -3400.,
Double_t  zE = 2650. 
)
inline

Definition at line 31 of file MuDISGenerator.h.

31 {
32 startZ = zS;
33 endZ = zE;
34 }

Member Data Documentation

◆ dPart

TClonesArray* MuDISGenerator::dPart
protected

Definition at line 43 of file MuDISGenerator.h.

◆ endZ

Double_t MuDISGenerator::endZ
protected

Definition at line 41 of file MuDISGenerator.h.

◆ fFirst

bool MuDISGenerator::fFirst
protected

Definition at line 49 of file MuDISGenerator.h.

◆ fInputFile

TFile* MuDISGenerator::fInputFile
protected

don't make it persistent, magic ROOT command

Definition at line 45 of file MuDISGenerator.h.

◆ fLogger

FairLogger* MuDISGenerator::fLogger
protected

Definition at line 44 of file MuDISGenerator.h.

◆ fn

int MuDISGenerator::fn
protected

Definition at line 48 of file MuDISGenerator.h.

◆ fNevents

int MuDISGenerator::fNevents
protected

Definition at line 47 of file MuDISGenerator.h.

◆ fTree

TTree* MuDISGenerator::fTree
protected

Definition at line 46 of file MuDISGenerator.h.

◆ iMuon

TClonesArray* MuDISGenerator::iMuon
protected

Definition at line 42 of file MuDISGenerator.h.

◆ startZ

Double_t MuDISGenerator::startZ
protected

Definition at line 41 of file MuDISGenerator.h.

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