LCDG4: LCDG4readStdFile class Reference

00022 {
00023   debflg = 0;
00024 }

00027 {
00028   debflg=0;
00029   m_mcpartmgr = mcp;
00030 
00031   // Do trivial validity checks on input file
00032   G4String check = "test -f "+G4String(evfile);
00033   if( system(check.data()) ) {
00034     G4cout <<" ***** Error: input file does not exist--> "<< evfile << G4endl;
00035     exit(1);
00036   }
00037 
00038   reader = new lStdHep(evfile,false);
00039   if(reader->getError()) {
00040     G4cerr << "Could not open file <" << evfile << ">" << G4endl;
00041     G4Exception("LCDG4readStdFile: cannot open file.");
00042   }
00043 }

00046                                     {
00047   if(reader) delete reader;
00048 }

00385                                                                                                              {
00386 
00387   // insert the particle...
00388   aColl.push_back( parent );
00389 
00390   // ...and all its daughters recursively
00391   G4PrimaryParticle* ppp = parent->GetDaughter();
00392   while( ppp ) {
00393     AddFamilyToCollection(aColl, ppp);
00394     ppp = ppp->GetNext();
00395   }
00396 }

00340                                                                   {
00341   bool ok = true;
00342 
00343   //*** Build list of all particles and their descendants to be simulated
00344   list<G4PrimaryParticle*> toSim;
00345 
00346   // loop over primary vertices
00347   for(list<G4PrimaryVertex*>::iterator ivtx = vColl.begin();
00348       ivtx!=vColl.end(); ivtx++) {
00349 
00350     // loop over particles attached to this PV    
00351     G4int Nparts = (*ivtx)->GetNumberOfParticle();
00352     for(G4int j=0; j<Nparts; j++) {
00353       G4PrimaryParticle* ppp = (*ivtx)->GetPrimary(j);
00354       // add all its daughters
00355       AddFamilyToCollection( toSim, ppp );
00356     }
00357   }
00358 
00359   // Make sure that all FINAL STATEs will be simulated
00360   list<G4PrimaryParticle*>::iterator ip;
00361   for(unsigned int i=0; i<HPlist.size(); i++) {
00362     if( abs(HPlist[i]->GetISTHEP()) == 1 ) {
00363       G4PrimaryParticle* iptest = HPlist[i]->GetTheParticle();
00364       // check particle is in list to be simulated
00365       bool found=false;
00366       for(ip=toSim.begin(); ip!=toSim.end(); ip++) {
00367         if((*ip)==iptest) {
00368           found=true;
00369           break;
00370         }
00371       }
00372       if(!found) {
00373         // a FINAL STATE particle was not found on the list!!
00374         ok = false;
00375         G4cout<<" ***** ReadStdFile::CheckVertices problem:\n"
00376               <<" A FINAL STATE particle is not getting simulated.\n"
00377               <<" Please let Guilherme (lima@fnal.gov) know about it."
00378               << G4endl;
00379       }
00380     }
00381   }
00382   return ok;
00383 }

00188                                                          {
00189 
00190   // For debugging a specifig event ==================
00191   // create file debevt.dat with just an event number inside it (start from 1).
00192   // This number will be read into ievtdgb.
00193   // debugging flag in McPartMgr will be set just for event chosen
00194   static int ievt=0;
00195   static int ievtdbg, iskip;
00196   if(ievt==0) {
00197     ifstream debevtfile("debevt.dat");
00198     int iseed;
00199     debevtfile>> ievtdbg >> iseed >> iskip;
00200     debevtfile.close();
00201   }
00202 
00203   //===================================================
00204   // Skip events if requested
00205   if(iskip>0) {
00206     G4cout<< "\n***** Skipping "<< iskip << " events from input.\n"<< G4endl;
00207     for(;;) {
00208       if (ReadEvent() == 0) {
00209         G4Exception("End-Of-File : HEPEvt input file");
00210         return;
00211       }
00212       ievt++;
00213       if(ievt>=iskip) break;
00214     }
00215     ievt=0;
00216     iskip=0;
00217   }
00218   //===================================================
00219 
00220   ievt++;
00221   debflg=0;
00222   if(ievt==ievtdbg) debflg=1;
00223 
00224   // read MC particle information from the stdhep file
00225   if (ReadEvent() == 0) {
00226     G4Exception("End-Of-File : HEPEvt input file");
00227     return;
00228   }
00229 
00230   // Builds one LCDMcPart for each particle in the stdhep common block
00231   MakeMcParts();
00232 
00233 //   int nhep = m_mcpartmgr->MCparticles.size();
00234 // G4cout<<" LCDG4readStdFile::GeneratePrimVertex nhep=" << nhep << G4endl;
00235 
00236   int tag = 1;  // particle tags start at one, not zero
00237 
00238   vector<LCDMcPart*>::const_iterator part;
00239   vector<LCDMcPart*>& mclist = m_mcpartmgr->MCparticles;
00240   for(part = mclist.begin(); part != mclist.end(); part++) {
00241     LCDMcPart* mco   = (*part);
00242     G4int      isthep=mco->GetGenStatus() ;         // gen status code
00243     G4int      idhep =mco->GetParticleID();         // PDG code
00244     G4double   phep1 =mco->Get4MomentumPtr()->px(); // px
00245     G4double   phep2 =mco->Get4MomentumPtr()->py(); // py
00246     G4double   phep3 =mco->Get4MomentumPtr()->pz(); // pz
00247     G4double   phep5 =mco->Get4MomentumPtr()->e (); // energy
00248 
00249     // parentage info
00250     G4int firstMom, lastMom;
00251     mco->GetMotherTags(firstMom,lastMom);
00252     G4int firstDau, lastDau;
00253     mco->GetDaughterTags(firstDau, lastDau);
00254     if(debflg) {
00255       if(tag==1) G4cout<< " tag      id    istat    daus      moms" << G4endl;
00256       G4cout << myFormat(4,tag)
00257              << myFormat(8,idhep)
00258              << myFormat(8,isthep)
00259              << myFormat(6,firstDau) <<'-'<< myFormat(3,lastDau)
00260              << myFormat(6,firstMom) <<'-'<< myFormat(3,lastMom)
00261              <<"  "<< myFormat(10,phep5/GeV)
00262              << G4endl;
00263 //        G4cout << ", pvec(GeV)=("<< phep1 << ", "<< phep2 <<", "<< phep3 << ")\n";
00264 //        G4cout<< ", pos=" << *mco->GetPositionPtr()/mm
00265 //          << ", end=" << *mco->GetTermPositionPtr()/mm;
00266 //        G4cout << G4endl << G4endl;
00267     }
00268 
00269     // create G4PrimaryParticle object
00270     G4PrimaryParticle* particle
00271       = new G4PrimaryParticle(idhep,phep1,phep2,phep3);
00272     double mass = sqrt(phep5*phep5 - phep1*phep1 - phep2*phep2 - phep3*phep3);
00273     particle->SetMass(mass);
00274 
00275     // create G4HEPEvtParticle object
00276     G4HEPEvtParticle* hepParticle;
00277     if(firstDau>lastDau) {
00278       // no daughters
00279       hepParticle = new G4HEPEvtParticle( particle, isthep, 0, 0);
00280     }
00281     else {
00282       // one or more daughters
00283       hepParticle = new G4HEPEvtParticle( particle, isthep, firstDau, lastDau);
00284     }
00285 
00286     // Store
00287     HPlist.push_back(hepParticle);
00288     tag++;
00289   }
00290 
00291   int npart=HPlist.size();
00292 
00293   // check if there is at least one particle
00294   if(npart == 0) return;
00295 
00296   // Build links from mother to daughters
00297   for( int i=0; i<npart; i++ ) {
00298     G4PrimaryParticle* mother = HPlist[i]->GetTheParticle();
00299 //     G4cout<<"mother "<< i+1 <<", id="<< mother->GetPDGcode()
00300 //        <<", jdahep="<< HPlist[i]->GetJDAHEP1()
00301 //        <<' '<< HPlist[i]->GetJDAHEP2() << G4endl;
00302     if( HPlist[i]->GetJDAHEP1() > 0 ) { //  it has daughters
00303       G4int jda1 = HPlist[i]->GetJDAHEP1()-1; // FORTRAN index starts from 1
00304       G4int jda2 = HPlist[i]->GetJDAHEP2()-1; // but C++ starts from 0.
00305 
00306       for( G4int j=jda1; j<=jda2; j++ ) {
00307         G4PrimaryParticle* daughter = HPlist[j]->GetTheParticle();
00308 //      G4cout<<"daughter "<< j+1 <<", id="<< daughter->GetPDGcode()
00309 //            <<", isthep="<< HPlist[j]->GetISTHEP()
00310 //            << G4endl;
00311 
00312         if(HPlist[j]->GetISTHEP()>0) {
00313           mother->SetDaughter( daughter );
00314           HPlist[j]->Done();  // this sets ISTHEP *= -1 (negative)
00315         }
00316       }
00317     }
00318   }
00319 
00320   // Get vertices containing particles for simulation
00321   list<G4PrimaryVertex*> vColl = GetVertices();
00322 
00323   // Make sure that all final state particles will be handled
00324   if( CheckVertices(vColl) ) {
00325     // Insert the vertices into G4Event object
00326     list<G4PrimaryVertex*>::const_iterator ivtx = vColl.begin();
00327     for( ivtx=vColl.begin(); ivtx!=vColl.end(); ivtx++ ) {
00328       evt->AddPrimaryVertex(*ivtx);
00329     }
00330   }
00331 
00332   // clear G4HEPEvtParticles
00333   for(unsigned int i=0; i<HPlist.size(); i++) {
00334     delete HPlist[i];
00335   }
00336   HPlist.clear();
00337 }

00403                                                      {
00404   // create G4PrimaryVertex object
00405   G4double x0 = 0.;  // vertex point X
00406   G4double y0 = 0.;  // vertex point Y
00407   G4double z0 = 0.;  // vertex point Z
00408   G4double t0 = 0.;  // vertex time
00409   G4PrimaryVertex* vertex = new G4PrimaryVertex(x0,y0,z0,t0);
00410 
00411   // put initial particles to the vertex
00412   for (unsigned int i=0; i< HPlist.size(); i++ ) {
00413     int hepevtStatus = HPlist[i]->GetISTHEP();
00414     // ISTHEP of daughters had been made negative
00415     bool attach=false;
00416     if( hepevtStatus>0) {
00417       if( hepevtStatus!= 3 ) attach=true;
00418       else
00419         // hepevtStatus==3 here
00420         if( HPlist[i]->GetJDAHEP1()>0 && HPlist[i]->GetJDAHEP2()>0 ) {
00421           // fixes some Pandora-Pythia ee->WW->leptons and ee->leptons procs
00422           attach=true;
00423         }
00424     }
00425     if(attach) {
00426       if(debflg) {
00427         G4cout <<"Attaching to primary vertex: tag="<< i+1 <<", id="
00428                << HPlist[i]->GetTheParticle()->GetPDGcode() << G4endl;
00429       }
00430       G4PrimaryParticle* initialParticle = HPlist[i]->GetTheParticle();
00431       vertex->SetPrimary(initialParticle);
00432     }
00433   }
00434 
00435   list<G4PrimaryVertex*> vtxColl;
00436   vtxColl.push_back(vertex);
00437   return vtxColl;
00438 }

00070 {
00071   // Make McPart objects from the stdHEP common
00072   LCDMcPart*            part;
00073   Hep3Vector            position;
00074   HepLorentzVector      momentum4;
00075 
00076   int newTag = 0;
00077   double totFSE = 0.0;
00078   int npart  = reader->nTracks();
00079   for(int ihep=0; ihep<npart; ++ihep) {
00080     if(ihep > 3999) continue;
00081 
00082     // Possible status are: 1=FINALSTATE, 2=INTERMEDIATE, 3=DOCUMENTATION
00083     G4int status = reader->status(ihep);
00084     G4int id = reader->pid(ihep);
00085     if(id==99999999) continue;  // do not create an McPart for this one
00086 
00087     position.setX( reader->X(ihep)*mm );
00088     position.setY( reader->Y(ihep)*mm );
00089     position.setZ( reader->Z(ihep)*mm );
00090     //     G4cout <<"Primary position = "<< position << G4endl;
00091 
00092     momentum4.setPx(reader->Px(ihep)*GeV);
00093     momentum4.setPy(reader->Py(ihep)*GeV);
00094     momentum4.setPz(reader->Pz(ihep)*GeV);
00095     momentum4.setE(reader->E(ihep)*GeV);
00096     if(status==1) totFSE += reader->E(ihep)*GeV;
00097 
00098     // Parent and daughter information
00099 
00100     // mother and daughter pointers: convert from Fortran to C++
00101     int jmohep1 = reader->mother1(ihep) - 1;
00102     int jmohep2 = reader->mother2(ihep) - 1;
00103     int jdahep1 = reader->daughter1(ihep)%10000 - 1;
00104     int jdahep2 = reader->daughter2(ihep)%10000 - 1;
00105     int parentID = reader->pid(jmohep1);
00106 
00107     // look for pathological scenarios here
00108     if( jmohep1>=npart || jmohep2>=npart ) {
00114       G4cout<<"\n***** Inconsistent hierarchy in input file!\n";
00115       G4cout<<"  Daughter indices pointing outside HEPEVT block\n";
00116       G4cout<<"  Processing a single neutrino instead (no hits)\n";
00117       G4cout<< G4endl;
00118 
00119       // clear MCParticle arrays
00120       m_mcpartmgr->Clear();
00121 
00122       // Make a new LCDMcPart
00123       part = new LCDMcPart();
00124       part->SetGenStatus(4);
00125       part->SetParticleID(12);
00126       part->SetChargeFromPid(12);
00127       part->SetPosition( Hep3Vector(0,0,0) );
00128       part->Set4Momentum( HepLorentzVector(0,0,10,10) );
00129       part->SetParticleTag(1);  // start at one, not zero
00130       part->SetMotherTags(-1,-1);
00131       part->SetDaughterTags(0,0);
00132       m_mcpartmgr->MCparticles.push_back(part);
00133       return 1;
00134     }
00135 
00136     // Make a new LCDMcPart
00137     newTag++;  // tags start at 1, not zero
00138     part = new LCDMcPart();
00139     part->SetGenStatus(status);
00140     part->SetParticleID(id);
00141     part->SetChargeFromPid(id);
00142     part->SetPosition(position);
00143     part->Set4Momentum(momentum4);
00144     part->SetParticleTag(newTag);  // start at one, not zero
00145 
00146     // Now set the tags.  Tags start at 1.  Note that there might be
00147     // dropped particles in the input (id=99999999), so we need an offset
00148     // which is typically different for single particles and physics samples
00149     if(newTag==1) Hep2TagOffset = newTag - ihep;
00150     G4int firstMotherTag = jmohep1 + Hep2TagOffset;
00151     G4int lastMotherTag  = jmohep2 + Hep2TagOffset;
00152     part->SetMotherTags(firstMotherTag, lastMotherTag);
00153     G4int firstDaughterTag = jdahep1 + Hep2TagOffset;
00154     G4int lastDaughterTag = jdahep2 + Hep2TagOffset;
00155     part->SetDaughterTags(firstDaughterTag, lastDaughterTag);
00156 
00157     // some particles have more than one parent
00158     if( debflg && lastMotherTag>0 ) {
00159       printf("2-moms particle: tag=%d, par1=%d, par2=%d\n",
00160              newTag,firstMotherTag,lastMotherTag);
00161     }
00162 
00163     // set TermPosition for all parents, as Geant4 does not update some of them
00164     if(firstMotherTag > 0 && parentID!=99999999) {
00165       part->SetParentTag(firstMotherTag);
00166       LCDMcPart* parent = m_mcpartmgr->MCparticles[firstMotherTag-1];
00167       // parents die at the starting point of their children
00168       parent->SetTermPosition( part->GetPosition() );
00169     }
00170 
00171     m_mcpartmgr->MCparticles.push_back(part);
00172   }
00173 
00174   //   i = 0;
00175   //   vector<LCDMcPart*>::const_iterator piter;
00176   //   for(piter = m_mcpartmgr->MCparticles.begin();
00177   //       piter != m_mcpartmgr->MCparticles.end(); piter++) {
00178   //     part = *piter;
00179 
00180   //     i++;
00181   //   }
00182 
00183 //   G4cout<< " total energy from Final State particles: "
00184 //      << totFSE/GeV <<"GeV"<< G4endl;
00185   return npart;
00186 }

00052 {
00053   if(ifile) ; // make compiler happy
00054   // reading next event from input binary StdHep file
00055   if(reader->more()) {
00056     if(int errorcode = reader->readEvent()) {
00057       G4cout << "*** Error "<< errorcode <<" reading STDHEP file!"<< G4endl;
00058       return 0;
00059     }
00060   }
00061   else {
00062     G4Exception("*** End-Of-File: HEPEvt input file");
00063     return 0;
00064   }
00065 
00066   return reader->nTracks();
00067 }

00037 { debflg = dbg; }


Public Member Functions
	LCDG4readStdFile ()
	LCDG4readStdFile (char evfile, LCDG4McPartManager mcp)
	~LCDG4readStdFile ()
void	GeneratePrimaryVertex (G4Event *evt)
char *	GetFileName (int i)
int	MakeMcParts ()
int	ReadEvent (int ifile=0)
void	SetDebugFlag (int dbg)
bool	CheckVertices (std::list< G4PrimaryVertex * > &vColl)
std::list< G4PrimaryVertex * >	GetVertices ()
void	AddFamilyToCollection (std::list< G4PrimaryParticle * > &aColl, G4PrimaryParticle *parent) const
Private Member Functions
	LCDG4readStdFile (const LCDG4readStdFile &)
LCDG4readStdFile &	operator= (const LCDG4readStdFile &)
Private Attributes
lStdHep *	reader
std::vector< G4HEPEvtParticle * >	HPlist
LCDG4McPartManager *	m_mcpartmgr
G4int	debflg
G4int	Hep2TagOffset

LCDG4readStdFile Class Reference

Public Member Functions

Private Member Functions

Private Attributes

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation