]>
Commit | Line | Data |
---|---|---|
cb3eff76 | 1 | // module identifier line... |
2 | // Author: Brian Thorsbro, 24/6-2014 | |
3 | ||
4 | #include <iostream> | |
5 | #include <sstream> | |
6 | #include <string> | |
7 | #include <list> | |
8 | #include <set> | |
9 | #include <time.h> | |
10 | ||
11 | #include "THepMCParser.h" | |
cb3eff76 | 12 | #include "TObject.h" |
13 | #include "TTree.h" | |
14 | #include "TClonesArray.h" | |
15 | #include "TFile.h" | |
16 | #include "TParticle.h" | |
17 | #include "TDatabasePDG.h" | |
71c3b2be | 18 | #include "HepMC/IO_GenEvent.h" |
cb3eff76 | 19 | |
cedc6926 | 20 | |
cb3eff76 | 21 | using namespace std; |
22 | ||
23 | ||
24 | THepMCParser::THepMCParser(const char * infile) : fTree(0) | |
25 | { | |
26 | HepMC::IO_BaseClass * events = new HepMC::IO_GenEvent(infile, std::ios::in); | |
27 | init(events); | |
28 | delete events; | |
29 | events = 0; // nullptr | |
30 | } | |
31 | THepMCParser::THepMCParser(HepMC::IO_BaseClass * events) : fTree(0) | |
32 | { | |
33 | init(events); | |
34 | } | |
35 | void THepMCParser::init(HepMC::IO_BaseClass * events) | |
36 | { | |
37 | int particlecount = 0; | |
38 | fTree = new TTree("treeEPOS","Tree EPOS"); | |
39 | TClonesArray * array = new TClonesArray("TParticle"); | |
40 | // array->BypassStreamer(); | |
41 | fTree->Branch("Particles",&array); // more flags? | |
42 | THepMCParser::HeavyIonHeader_t heavyIonHeader; | |
43 | fTree->Branch("HeavyIonInfo", &heavyIonHeader, THepMCParser::fgHeavyIonHeaderBranchString); | |
44 | THepMCParser::PdfHeader_t pdfHeader; | |
45 | fTree->Branch("PdfInfo", &pdfHeader, THepMCParser::fgPdfHeaderBranchString); | |
46 | HepMC::GenEvent* evt = 0; // nullptr | |
47 | while ((evt = events->read_next_event())) { | |
48 | string errMsg1 = ParseGenEvent2TCloneArray(evt,array); | |
49 | string errMsg2 = ParseGenEvent2HeaderStructs(evt,heavyIonHeader,pdfHeader); | |
50 | if (errMsg1.length() == 0 && errMsg2.length() == 0) { | |
51 | fTree->Fill(); | |
52 | } else { | |
53 | if (errMsg1.length() != 0) cerr << errMsg1 << endl; | |
54 | if (errMsg2.length() != 0) cerr << errMsg2 << endl; | |
55 | } | |
56 | particlecount += array->Capacity(); | |
57 | } | |
58 | // array->Clear(); | |
59 | // delete array; | |
60 | // array = 0; // nullptr | |
61 | cout << " parsed " << particlecount << " particles" << endl; | |
62 | } | |
63 | ||
64 | ||
65 | TTree * THepMCParser::GetTTree() | |
66 | { | |
67 | return fTree; | |
68 | } | |
69 | void THepMCParser::WriteTTreeToFile(const char *outfile) | |
70 | { | |
71 | TFile * f = new TFile(outfile, "recreate"); | |
72 | fTree->Write(); | |
73 | delete f; | |
74 | f = 0; // nullptr | |
75 | } | |
76 | ||
77 | ||
78 | ||
79 | // Based on a validator written by Peter Hristov, CERN | |
80 | bool THepMCParser::IsValidMotherDaughtersConsitency(bool useStdErr, bool requireSecondMotherBeforeDaughters) | |
81 | { | |
82 | bool valid = true; | |
83 | TClonesArray * array = new TClonesArray("TParticle"); | |
84 | TBranch* branch = fTree->GetBranch("Particles"); | |
85 | branch->SetAddress(&array); | |
86 | Int_t count = branch->GetEntries(); | |
87 | for (Int_t idx=0; idx<count; ++idx) { | |
88 | array->Clear(); | |
89 | branch->GetEntry(idx); // "fill" the array | |
90 | Int_t nkeep = array->GetEntriesFast(); | |
91 | for (Int_t i=0; i<nkeep; i++) { | |
92 | TParticle * part = (TParticle*)array->AddrAt(i); | |
93 | Int_t mum1 = part->GetFirstMother(); | |
94 | Int_t mum2 = part->GetSecondMother(); | |
95 | Int_t fd = part->GetFirstDaughter(); | |
96 | Int_t ld = part->GetLastDaughter(); | |
97 | if (mum1>-1 && i<mum1) { | |
98 | valid = false; | |
99 | if (useStdErr) cerr << "Problem: first_mother(" << mum1 << ") after daughter(" << i << ")" << endl; | |
100 | } | |
101 | if (mum2>-1 && i<mum2 && requireSecondMotherBeforeDaughters) { | |
102 | valid = false; | |
103 | if (useStdErr) cerr << "Problem: second_mother(" << mum2 << ") after daughter(" << i << ")" << endl; | |
104 | } | |
105 | if (fd > ld ) { | |
106 | valid = false; | |
107 | if (useStdErr) cerr << "Problem: first_daughter(" << fd << ") > last_daughter(" << ld << ")" << endl; | |
108 | } | |
109 | for (Int_t id=TMath::Max(fd,0); id<=ld; id++) { | |
110 | TParticle * daughter = (TParticle*)array->AddrAt(id); | |
111 | if (daughter->GetFirstMother() != i && daughter->GetSecondMother() != i) { | |
112 | valid = false; | |
113 | if (useStdErr) cerr << "Problem: mother("<< i << ") not registered as mother for either first_daughter(" | |
114 | << daughter->GetFirstMother() << ") or second_daughter(" | |
115 | << daughter->GetSecondMother() << ")" << endl; | |
116 | } | |
117 | } | |
118 | } | |
119 | } | |
120 | delete array; | |
121 | array = 0; | |
122 | return valid; | |
123 | } | |
124 | ||
125 | bool THepMCParser::IsValidParticleInvariantMass(bool useStdErr, bool includeStatusCode2Particles) | |
126 | { | |
127 | bool valid = true; | |
128 | TClonesArray *array = new TClonesArray("TParticle"); | |
129 | TBranch* branch = fTree->GetBranch("Particles"); | |
130 | branch->SetAddress(&array); | |
131 | Int_t count = branch->GetEntries(); | |
132 | for (Int_t idx=0; idx<count; ++idx) { | |
133 | array->Clear(); | |
134 | branch->GetEntry(idx); | |
135 | Int_t nkeep = array->GetEntries(); | |
136 | for (Int_t i=0; i<nkeep; i++) { | |
137 | TParticle * parton = (TParticle*)array->AddrAt(i); | |
138 | if (parton->GetStatusCode()==2 && !includeStatusCode2Particles) { | |
139 | continue; | |
140 | } | |
141 | TLorentzVector v; | |
142 | parton->Momentum(v); | |
143 | Double_t m1 = v.M(); // invariant mass from the particle in the HepMC event | |
144 | TParticlePDG *dbParton = parton->GetPDG(); | |
145 | if (!dbParton) { | |
146 | if (useStdErr) cerr << "Warning: could not look up particle with PDG Code: " << parton->GetPdgCode() << endl << endl; | |
147 | continue; | |
148 | } | |
149 | Double_t m2 = dbParton->Mass(); | |
150 | bool checkok; | |
151 | if (m2 == 0) { | |
152 | checkok = abs(m1) < 0.0001; // no such thing as negative mass... | |
153 | } else { | |
154 | checkok = abs(1 - m1/m2) < 0.01; | |
155 | } | |
156 | if (!checkok && useStdErr) { | |
157 | cerr << "Problem: " << GetParticleName(parton) << " HepMC:" << m1 << endl; | |
158 | cerr << ListReactionChain(array,i); | |
159 | cerr << endl; | |
160 | } | |
161 | if (!checkok) | |
162 | valid = false; | |
163 | } | |
164 | } | |
165 | delete array; | |
166 | array = 0; | |
167 | return valid; | |
168 | } | |
169 | ||
170 | bool THepMCParser::IsValidVertexInvariantMass(bool useStdErr, bool includeStatusCode2Particles) | |
171 | { | |
172 | bool valid = true; | |
173 | TClonesArray * array = new TClonesArray("TParticle"); | |
174 | TBranch* branch = fTree->GetBranch("Particles"); | |
175 | branch->SetAddress(&array); | |
176 | Int_t count = branch->GetEntries(); | |
177 | for (Int_t idx=0; idx<count; ++idx) { | |
178 | array->Clear(); | |
179 | branch->GetEntry(idx); // "fill" the array | |
180 | TLorentzVector v_st1; | |
181 | TLorentzVector v_st4; | |
182 | Int_t nkeep = array->GetEntriesFast(); | |
183 | for (Int_t i=0; i<nkeep; i++) { | |
184 | TParticle * parton = (TParticle*)array->AddrAt(i); | |
185 | TLorentzVector v_in; | |
186 | parton->Momentum(v_in); | |
187 | if (parton->GetStatusCode()==4) { | |
188 | v_st4 += v_in; | |
189 | } else if (parton->GetStatusCode()==1) { | |
190 | v_st1 += v_in; | |
191 | } | |
192 | if (!includeStatusCode2Particles) { // only check beam particles vs final particles | |
193 | continue; | |
194 | } | |
195 | Int_t fd = parton->GetFirstDaughter(); | |
196 | Int_t ld = parton->GetLastDaughter(); | |
197 | if (fd == -1) continue; // no daughters, continue loop | |
198 | Int_t mother2 = -1; | |
199 | TLorentzVector v_out; | |
200 | bool oneok = false; bool allok = true; bool agreemother2 = true; | |
201 | ostringstream daughterpdg; | |
202 | ostringstream motherpdg; | |
203 | for (Int_t id=TMath::Max(fd,0); id<=ld; id++) { | |
204 | TParticle * daughter = (TParticle*)array->AddrAt(id); | |
205 | if (fd==id) { | |
206 | daughter->Momentum(v_out); | |
207 | mother2 = daughter->GetSecondMother(); | |
208 | } else { | |
209 | TLorentzVector d; | |
210 | daughter->Momentum(d); | |
211 | v_out += d; | |
212 | if (daughter->GetSecondMother() != mother2) agreemother2 = false; | |
213 | } | |
214 | if (daughter->GetFirstMother() == i) { | |
215 | oneok = true; | |
216 | } else { | |
217 | allok = false; | |
218 | } | |
219 | daughterpdg << " " << daughter->GetPdgCode(); | |
220 | } | |
221 | motherpdg << " " << parton->GetPdgCode(); | |
222 | if (mother2 > -1 && agreemother2) { | |
223 | TParticle * m2 = (TParticle*)array->AddrAt(mother2); | |
224 | TLorentzVector m2v; | |
225 | m2->Momentum(m2v); | |
226 | v_in += m2v; | |
227 | motherpdg << " " << m2->GetPdgCode(); | |
228 | } | |
229 | if (oneok && allok && agreemother2) { | |
230 | bool checkok = abs(1 - v_in.M()/v_out.M()) < 0.1; | |
231 | if (!checkok) valid=false; | |
232 | if (!checkok && useStdErr) { | |
233 | // cerr << "Problem: " << i << "[" << fd << "," << ld << "] PDG:" << motherpdg.str() << " ->" << daughterpdg.str() << " Inv.mass: " << v_in.M() << " -> " << v_out.M() << endl; | |
234 | cerr << ListReactionChain(array,i); | |
235 | cerr << endl; | |
236 | // cerr << v_in.Px() << " " << v_in.Py() << " " << v_in.Pz() << " " << v_in.E() << endl; | |
237 | } else if (useStdErr) { | |
238 | // cerr << "OK: " << i << "[" << fd << "," << ld << "] PDG:" << motherpdg.str() << " ->" << daughterpdg.str() << " Inv.mass: " << v_in.M() << " -> " << v_out.M() << endl; | |
239 | } | |
240 | } | |
241 | } | |
242 | bool checkok = abs(1 - v_st4.M()/v_st1.M()) < 0.001; | |
243 | if (!checkok) valid=false; | |
244 | if (!checkok && useStdErr) { | |
245 | cerr << " BEAM PARTICLES -> FINAL PARTICLES " << endl; | |
246 | cerr << " status 4 (" << v_st4.M() << ") -> status 1 (" << v_st1.M() << ")" << endl << endl; | |
247 | } | |
248 | } | |
249 | delete array; | |
250 | array = 0; | |
251 | return valid; | |
252 | } | |
253 | ||
254 | string THepMCParser::GetParticleName(TParticle * thisPart) | |
255 | { | |
256 | TParticlePDG *dbPart = thisPart->GetPDG(); | |
257 | ostringstream name; | |
258 | if (dbPart) { | |
259 | name << dbPart->GetName() << "(" << dbPart->PdgCode() << "," << dbPart->Mass() << ")"; | |
260 | } else { | |
261 | name << thisPart->GetPdgCode() << "(NoDBinfo)"; | |
262 | } | |
263 | return name.str(); | |
264 | } | |
265 | ||
266 | string THepMCParser::ListReactionChain(TClonesArray * particles, Int_t particleId) | |
267 | { | |
268 | ostringstream output; | |
269 | ||
270 | TParticle * part = (TParticle*)particles->AddrAt(particleId); | |
271 | Int_t m1id = part->GetFirstMother(); | |
272 | Int_t m2id = part->GetSecondMother(); | |
273 | if (m1id > 1) { // ignore the initial collision with beam particles | |
274 | ostringstream inStr; | |
275 | ostringstream outStr; | |
276 | TParticle * m1 = (TParticle*)particles->AddrAt(m1id); | |
277 | TLorentzVector v_in; | |
278 | m1->Momentum(v_in); | |
279 | inStr << GetParticleName(m1) << "[" << v_in.M() << "]"; | |
280 | if (m2id > 1) { | |
281 | TParticle * m2 = (TParticle*)particles->AddrAt(m2id); | |
282 | TLorentzVector v_m2; | |
283 | m2->Momentum(v_m2); | |
284 | v_in += v_m2; | |
285 | inStr << " " << GetParticleName(m2) << "[" << v_m2.M() << "]"; | |
286 | } | |
287 | Int_t fd = m1->GetFirstDaughter(); | |
288 | Int_t ld = m1->GetLastDaughter(); | |
289 | TLorentzVector v_out; | |
290 | part->Momentum(v_out); | |
291 | outStr << GetParticleName(part) << "[" << v_out.M() << "]"; | |
292 | for (Int_t i=fd; i<=ld; ++i) { | |
293 | if (i!=particleId) { | |
294 | TParticle * d = (TParticle*)particles->AddrAt(i); | |
295 | TLorentzVector v_d; | |
296 | d->Momentum(v_d); | |
297 | v_out += v_d; | |
298 | outStr << " " << GetParticleName(d) << "[" << v_d.M() << "]"; | |
299 | } | |
300 | } | |
301 | output << "Parent reaction, inv mass: " << v_in.M() << " -> " << v_out.M() << endl; | |
302 | output << " - partons: " << inStr.str() << " -> " << outStr.str() << endl; | |
303 | } | |
304 | Int_t fd = part->GetFirstDaughter(); | |
305 | Int_t ld = part->GetLastDaughter(); | |
306 | if (fd > -1) { | |
307 | ostringstream inStr; | |
308 | ostringstream outStr; | |
309 | TLorentzVector v_in; | |
310 | part->Momentum(v_in); | |
311 | inStr << GetParticleName(part) << "[" << v_in.M() << "]"; | |
312 | ||
313 | TParticle * f = (TParticle*)particles->AddrAt(fd); | |
314 | m2id = f->GetSecondMother(); | |
315 | if (m2id == particleId) { | |
316 | m2id = f->GetFirstMother(); | |
317 | } | |
318 | if (m2id > -1) { | |
319 | TParticle * m2 = (TParticle*)particles->AddrAt(m2id); | |
320 | TLorentzVector v_m2; | |
321 | m2->Momentum(v_m2); | |
322 | v_in += v_m2; | |
323 | inStr << " " << GetParticleName(m2) << "[" << v_m2.M() << "]"; | |
324 | } | |
325 | TLorentzVector v_out; | |
326 | f->Momentum(v_out); | |
327 | outStr << GetParticleName(f) << "[" << v_out.M() << "]"; | |
328 | for (Int_t i=fd+1; i<=ld; ++i) { | |
329 | TParticle * d = (TParticle*)particles->AddrAt(i); | |
330 | TLorentzVector v_d; | |
331 | d->Momentum(v_d); | |
332 | v_out += v_d; | |
333 | outStr << " " << GetParticleName(d) << "[" << v_d.M() << "]"; | |
334 | } | |
335 | output << "Child reaction, inv mass: " << v_in.M() << " -> " << v_out.M() << endl; | |
336 | output << " - partons: " << inStr.str() << " -> " << outStr.str() << endl; | |
337 | } else { | |
338 | output << "Child reaction" << endl << " - none" << endl; | |
339 | } | |
340 | ||
341 | return output.str(); | |
342 | } | |
343 | ||
344 | ||
cedc6926 | 345 | string THepMCParser::ParseGenEvent2TCloneArray(HepMC::GenEvent * genEvent, TClonesArray * array, string momUnit, string lenUnit, bool requireSecondMotherBeforeDaughters) |
cb3eff76 | 346 | { |
347 | if (requireSecondMotherBeforeDaughters) { | |
348 | return "requireSecondMotherBeforeDaughters not implemented yet!"; | |
349 | } | |
cedc6926 | 350 | genEvent->use_units(momUnit, lenUnit); |
cb3eff76 | 351 | array->Clear(); |
352 | ostringstream errMsgStream; | |
353 | map<int,Int_t> partonMemory; // unordered_map in c++11 - but probably not much performance gain from that: log(n) vs log(1) where constant can be high | |
354 | ||
355 | // Check event with HepMC's internal validation algorithm | |
356 | if (!genEvent->is_valid()) { | |
357 | errMsgStream << "Error with event id: " << genEvent->event_number() << ", event is not valid!"; | |
358 | return errMsgStream.str(); | |
359 | } | |
360 | ||
361 | // Pull out the beam particles from the event | |
362 | const pair<HepMC::GenParticle *,HepMC::GenParticle *> beamparts = genEvent->beam_particles(); | |
363 | ||
364 | // Four sanity checks: | |
365 | // - Beam particles exists and are not the same | |
366 | // - Both beam particles should have no production vertices, they come from the beams | |
367 | // - Both beam particles should have defined end vertices, as they both should contribute | |
368 | // - Both beam particles should have the exact same end vertex | |
369 | if (!beamparts.first || !beamparts.second || beamparts.first->barcode()==beamparts.second->barcode()) { | |
370 | errMsgStream << "Error with event id: " << genEvent->event_number() << ", beam particles doesn't exists or are the same"; | |
371 | return errMsgStream.str(); | |
372 | } | |
373 | if (beamparts.first->production_vertex() || beamparts.second->production_vertex()) { | |
374 | errMsgStream << "Error with event id: " << genEvent->event_number() << ", beam particles have production vertex/vertices..."; | |
375 | return errMsgStream.str(); | |
376 | } | |
377 | if (!beamparts.first->end_vertex() || !beamparts.second->end_vertex()) { | |
378 | errMsgStream << "Error with event id: " << genEvent->event_number() << ", beam particles have undefined end vertex/vertices..."; | |
379 | return errMsgStream.str(); | |
380 | } | |
381 | if (beamparts.first->end_vertex() != beamparts.second->end_vertex()) { | |
382 | errMsgStream << "Error with event id: " << genEvent->event_number() << ", beam particles do not collide in the same end vertex."; | |
383 | return errMsgStream.str(); | |
384 | } | |
385 | ||
386 | // Set the array to hold the number of particles in the event | |
387 | array->Expand(genEvent->particles_size()); | |
388 | ||
389 | // Create a TParticle for each beam particle | |
390 | new((*array)[0]) TParticle( | |
391 | beamparts.first->pdg_id(), | |
392 | beamparts.first->status(), // check if status has the same meaning | |
393 | -1, // no mother1 | |
394 | -1, // no mother2 | |
395 | -1, // first daughter not known yet | |
396 | -1, // last daughter not known yet | |
397 | beamparts.first->momentum().px(), | |
398 | beamparts.first->momentum().py(), | |
399 | beamparts.first->momentum().pz(), | |
400 | beamparts.first->momentum().e(), | |
401 | 0, // no production vertex, so zero? | |
402 | 0, | |
403 | 0, | |
404 | 0 | |
405 | ); | |
406 | partonMemory[beamparts.first->barcode()] = 0; | |
407 | new((*array)[1]) TParticle( | |
408 | beamparts.second->pdg_id(), | |
409 | beamparts.second->status(), | |
410 | -1, // no mother1 | |
411 | -1, // no mother2 | |
412 | -1, // first daughter not known yet | |
413 | -1, // last daughter not known yet | |
414 | beamparts.second->momentum().px(), | |
415 | beamparts.second->momentum().py(), | |
416 | beamparts.second->momentum().pz(), | |
417 | beamparts.second->momentum().e(), | |
418 | 0, // no production vertex, so zero? | |
419 | 0, | |
420 | 0, | |
421 | 0 | |
422 | ); | |
423 | partonMemory[beamparts.second->barcode()] = 1; | |
424 | Int_t arrayID = 2; // start counting IDs after the beam particles | |
425 | ||
426 | // Do first vertex as a special case for performance, since its often has the most daughters and both mothers are known | |
427 | Int_t firstDaughter = arrayID; | |
428 | for (HepMC::GenVertex::particles_out_const_iterator iter = beamparts.first->end_vertex()->particles_out_const_begin(); | |
429 | iter != beamparts.first->end_vertex()->particles_out_const_end(); | |
430 | ++iter) { | |
431 | new((*array)[arrayID]) TParticle( | |
432 | (*iter)->pdg_id(), | |
433 | (*iter)->status(), | |
434 | 0, // beam particle 1 | |
435 | 1, // beam particle 2 | |
436 | -1, // first daughter not known yet | |
437 | -1, // last daughter not known yet | |
438 | (*iter)->momentum().px(), | |
439 | (*iter)->momentum().py(), | |
440 | (*iter)->momentum().pz(), | |
441 | (*iter)->momentum().e(), | |
442 | beamparts.first->end_vertex()->position().x(), | |
443 | beamparts.first->end_vertex()->position().y(), | |
444 | beamparts.first->end_vertex()->position().z(), | |
445 | beamparts.first->end_vertex()->position().t() | |
446 | ); | |
447 | partonMemory[(*iter)->barcode()] = arrayID; | |
448 | ++arrayID; | |
449 | } | |
450 | Int_t lastDaughter = arrayID-1; | |
451 | ((TParticle*)array->AddrAt(0))->SetFirstDaughter(firstDaughter); // beam particle 1 | |
452 | ((TParticle*)array->AddrAt(0))->SetLastDaughter(lastDaughter); | |
453 | ((TParticle*)array->AddrAt(1))->SetFirstDaughter(firstDaughter); // beam particle 2 | |
454 | ((TParticle*)array->AddrAt(1))->SetLastDaughter(lastDaughter); | |
455 | ||
456 | // Build vertex list by exploring tree and sorting such that daughters comes after mothers | |
457 | list<HepMC::GenVertex*> vertexList; | |
458 | set<int> vertexSearchSet; | |
459 | ExploreVertex(beamparts.first->end_vertex(),vertexList,vertexSearchSet,requireSecondMotherBeforeDaughters); | |
460 | ||
461 | // Analyze each vertex | |
462 | for (list<HepMC::GenVertex*>::iterator i = vertexList.begin(); i != vertexList.end(); ++i) { | |
463 | HepMC::GenVertex * vertex = (*i); | |
464 | ||
465 | // first establish mother relations | |
466 | HepMC::GenVertex::particles_in_const_iterator iter = vertex->particles_in_const_begin(); | |
467 | if (iter == vertex->particles_in_const_end()) { | |
468 | return "Particle without a mother, and its not a beam particle!"; | |
469 | } | |
470 | int motherA = partonMemory[(*iter)->barcode()]; | |
471 | if (((TParticle*)array->AddrAt(motherA))->GetFirstDaughter() > -1) { | |
472 | return "Trying to assign new daughters to a particle that already has daughters defined!"; | |
473 | } | |
474 | ++iter; | |
475 | int motherB = -1; | |
476 | if (iter != vertex->particles_in_const_end()) { | |
477 | motherB = partonMemory[(*iter)->barcode()]; | |
478 | if (((TParticle*)array->AddrAt(motherB))->GetFirstDaughter() > -1) { | |
479 | return "Trying to assign new daughters to a particle that already has daughters defined!"; | |
480 | } | |
481 | ++iter; | |
482 | if (iter != vertex->particles_in_const_end()) { | |
483 | return "Particle with more than two mothers!"; | |
484 | } | |
485 | } | |
486 | if (motherB > -1 && motherB < motherA) { | |
487 | int swap = motherA; motherA = motherB; motherB = swap; | |
488 | } | |
489 | ||
490 | // add the particles to the array, important that they are add in succession with respect to arrayID | |
491 | firstDaughter = arrayID; | |
492 | for (iter = vertex->particles_out_const_begin(); | |
493 | iter != vertex->particles_out_const_end(); | |
494 | ++iter) { | |
495 | new((*array)[arrayID]) TParticle( | |
496 | (*iter)->pdg_id(), | |
497 | (*iter)->status(), | |
498 | motherA, // mother 1 | |
499 | motherB, // mother 2 | |
500 | -1, // first daughter, if applicable, not known yet | |
501 | -1, // last daughter, if applicable, not known yet | |
502 | (*iter)->momentum().px(), | |
503 | (*iter)->momentum().py(), | |
504 | (*iter)->momentum().pz(), | |
505 | (*iter)->momentum().e(), | |
506 | vertex->position().x(), | |
507 | vertex->position().y(), | |
508 | vertex->position().z(), | |
509 | vertex->position().t() | |
510 | ); | |
511 | partonMemory[(*iter)->barcode()] = arrayID; | |
512 | ++arrayID; | |
513 | } | |
514 | lastDaughter = arrayID-1; | |
515 | if (lastDaughter < firstDaughter) { | |
516 | return "Vertex with no out particles, should not be possible!"; | |
517 | } | |
518 | // update mother with daughter interval | |
519 | ((TParticle*)array->AddrAt(motherA))->SetFirstDaughter(firstDaughter); | |
520 | ((TParticle*)array->AddrAt(motherA))->SetLastDaughter(lastDaughter); | |
521 | if (motherB > -1) { | |
522 | ((TParticle*)array->AddrAt(motherB))->SetFirstDaughter(firstDaughter); | |
523 | ((TParticle*)array->AddrAt(motherB))->SetLastDaughter(lastDaughter); | |
524 | } | |
525 | } | |
526 | ||
527 | return ""; | |
528 | } | |
529 | ||
530 | ||
531 | void THepMCParser::ExploreVertex(HepMC::GenVertex * vertex, list<HepMC::GenVertex*> & vertexList, set<int> & vertexSearchSet, bool requireSecondMotherBeforeDaughters) | |
532 | { | |
533 | for (HepMC::GenVertex::particles_out_const_iterator iter = vertex->particles_out_const_begin(); | |
534 | iter != vertex->particles_out_const_end(); | |
535 | ++iter) { | |
536 | HepMC::GenVertex * testVertex = (*iter)->end_vertex(); | |
537 | if (testVertex) { | |
538 | bool foundVertex = vertexSearchSet.find((*iter)->end_vertex()->barcode()) != vertexSearchSet.end(); | |
539 | if (requireSecondMotherBeforeDaughters) { | |
540 | // redo this algorithem to move subtree instead of node.... | |
541 | // its not completely error proof in its current implementation even though the error is extremely rare | |
542 | ||
543 | if (foundVertex) for (list<HepMC::GenVertex*>::iterator i = vertexList.begin(); i != vertexList.end(); ++i) { | |
544 | if ((*i)->barcode() == testVertex->barcode()) { | |
545 | vertexList.erase(i); | |
546 | cout << " it happened, the vertex parsing order had to be changed " << endl; | |
547 | break; | |
548 | } | |
549 | } else { | |
550 | vertexSearchSet.insert((*iter)->end_vertex()->barcode()); | |
551 | } | |
552 | vertexList.push_back(testVertex); | |
553 | if (!foundVertex) ExploreVertex(testVertex,vertexList,vertexSearchSet,requireSecondMotherBeforeDaughters); | |
554 | ||
555 | } else { | |
556 | if (!foundVertex) { | |
557 | vertexSearchSet.insert((*iter)->end_vertex()->barcode()); | |
558 | vertexList.push_back(testVertex); | |
559 | ExploreVertex(testVertex,vertexList,vertexSearchSet,requireSecondMotherBeforeDaughters); | |
560 | } | |
561 | } | |
562 | } | |
563 | } | |
564 | } | |
565 | ||
566 | ||
567 | ||
568 | const char * THepMCParser::fgHeavyIonHeaderBranchString = "Ncoll_hard/I:Npart_proj:Npart_targ:Ncoll:spectator_neutrons:spectator_protons:N_Nwounded_collisions:Nwounded_N_collisions:Nwounded_Nwounded_collisions:impact_parameter/F:event_plane_angle:eccentricity:sigma_inel_NN"; | |
569 | const char * THepMCParser::fgPdfHeaderBranchString = "id1/I:id2:pdf_id1:pdf_id2:x1/D:x2:scalePDF:pdf1:pdf2"; | |
570 | ||
571 | string THepMCParser::ParseGenEvent2HeaderStructs(HepMC::GenEvent * genEvent, HeavyIonHeader_t & heavyIonHeader, PdfHeader_t & pdfHeader, bool fillZeroOnMissingHeavyIon, bool fillZeroOnMissingPdf) | |
572 | { | |
573 | HepMC::HeavyIon * heavyIon = genEvent->heavy_ion(); | |
574 | HepMC::PdfInfo * pdfInfo = genEvent->pdf_info(); | |
575 | if ((!heavyIon && !fillZeroOnMissingHeavyIon) || (!pdfInfo && !fillZeroOnMissingPdf)) { | |
576 | return "HeavyIonInfo and/or PdfInfo not defined for this event, did you read it with IO_GenEvent?"; | |
577 | } | |
578 | if (heavyIon) { | |
579 | heavyIonHeader.Ncoll_hard = heavyIon->Ncoll_hard(); | |
580 | heavyIonHeader.Npart_proj = heavyIon->Npart_proj(); | |
581 | heavyIonHeader.Npart_targ = heavyIon->Npart_targ(); | |
582 | heavyIonHeader.Ncoll = heavyIon->Ncoll(); | |
583 | heavyIonHeader.spectator_neutrons = heavyIon->spectator_neutrons(); | |
584 | heavyIonHeader.spectator_protons = heavyIon->spectator_protons(); | |
585 | heavyIonHeader.N_Nwounded_collisions = heavyIon->N_Nwounded_collisions(); | |
586 | heavyIonHeader.Nwounded_N_collisions = heavyIon->Nwounded_N_collisions(); | |
587 | heavyIonHeader.Nwounded_Nwounded_collisions = heavyIon->Nwounded_Nwounded_collisions(); | |
588 | heavyIonHeader.impact_parameter = heavyIon->impact_parameter(); | |
589 | heavyIonHeader.event_plane_angle = heavyIon->event_plane_angle(); | |
590 | heavyIonHeader.eccentricity = heavyIon->eccentricity(); | |
591 | heavyIonHeader.sigma_inel_NN = heavyIon->sigma_inel_NN(); | |
592 | } else { | |
593 | heavyIonHeader.Ncoll_hard = 0; | |
594 | heavyIonHeader.Npart_proj = 0; | |
595 | heavyIonHeader.Npart_targ = 0; | |
596 | heavyIonHeader.Ncoll = 0; | |
597 | heavyIonHeader.spectator_neutrons = 0; | |
598 | heavyIonHeader.spectator_protons = 0; | |
599 | heavyIonHeader.N_Nwounded_collisions = 0; | |
600 | heavyIonHeader.Nwounded_N_collisions = 0; | |
601 | heavyIonHeader.Nwounded_Nwounded_collisions = 0; | |
602 | heavyIonHeader.impact_parameter = 0.0; | |
603 | heavyIonHeader.event_plane_angle = 0.0; | |
604 | heavyIonHeader.eccentricity = 0.0; | |
605 | heavyIonHeader.sigma_inel_NN = 0.0; | |
606 | } | |
607 | if (pdfInfo) { | |
608 | pdfHeader.id1 = pdfInfo->id1(); | |
609 | pdfHeader.id2 = pdfInfo->id2(); | |
610 | pdfHeader.pdf_id1 = pdfInfo->pdf_id1(); | |
611 | pdfHeader.pdf_id2 = pdfInfo->pdf_id2(); | |
612 | pdfHeader.x1 = pdfInfo->x1(); | |
613 | pdfHeader.x2 = pdfInfo->x2(); | |
614 | pdfHeader.scalePDF = pdfInfo->scalePDF(); | |
615 | pdfHeader.pdf1 = pdfInfo->pdf1(); | |
616 | pdfHeader.pdf2 = pdfInfo->pdf2(); | |
617 | } else { | |
618 | pdfHeader.id1 = 0; | |
619 | pdfHeader.id2 = 0; | |
620 | pdfHeader.pdf_id1 = 0; | |
621 | pdfHeader.pdf_id2 = 0; | |
622 | pdfHeader.x1 = 0.0; | |
623 | pdfHeader.x2 = 0.0; | |
624 | pdfHeader.scalePDF = 0.0; | |
625 | pdfHeader.pdf1 = 0.0; | |
626 | pdfHeader.pdf2 = 0.0; | |
627 | } | |
628 | return ""; | |
629 | } | |
630 | ||
631 | ||
632 | ||
633 | ||
634 | ||
635 | ||
636 |