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