Fork me on GitHub

Changeset 380 in svn


Ignore:
Timestamp:
May 12, 2009, 9:47:12 AM (15 years ago)
Author:
Xavier Rouby
Message:

new PDG table

Location:
trunk
Files:
1 added
18 edited

Legend:

Unmodified
Added
Removed

  • trunk/CHANGELOG

    r373 r380  
    66        - E_cut for n/gamma reconstruction
    77   - time resolution also for RP220,RP420
    8    - Distinction between forward detectors with (TRootForwardTagger) /without (TRootRomanPot) time measurement
     8   - new : PDG tables for particles
    99     
    1010

  • trunk/Delphes.cpp

    r372 r380  
    213213  //Smearing information
    214214  RESOLution *DET = new RESOLution();
     215
    215216  cout <<"**                                                                 **"<< endl;
    216217  cout <<"**        ####### Start reading DETECTOR parameters #######        **"<< endl;
    … …  
    224225       << right << setw(2) <<"**"<<""<<endl;
    225226  cout <<"**                                                                 **"<< endl;
     227  DET->ReadParticleDataGroupTable();
     228  //DET->PDGtable.print();
    226229 
    227230  //Trigger information
    … …  
    260263      cout <<"**                 StdHEP file format detected                     **"<<endl;
    261264      cout <<"**                This can take several minutes                    **"<< endl;
    262       STDHEPConverter converter(inputFileList,outputfilename);//case ntpl file in input list
     265      STDHEPConverter converter(inputFileList,outputfilename,DET->PDGtable);//case ntpl file in input list
    263266    }
    264267  else if(line.length() == 1+line.find_last_of(".lhe"))
    … …  
    266269      cout <<"**                   LHEF file format detected                     **"<<endl;
    267270      cout <<"**                 This can take several minutes                   **"<< endl;
    268       LHEFConverter converter(inputFileList,outputfilename);//case ntpl file in input list
     271      LHEFConverter converter(inputFileList,outputfilename,DET->PDGtable);//case ntpl file in input list
    269272    }
    270273  else if(line.length() == 1+line.find_last_of(".root"))
    … …  
    272275      cout <<"**                   h2root file format detected                   **"<<endl;
    273276      cout <<"**                  This can take several minutes                  **"<< endl;
    274       HEPEVTConverter converter(inputFileList,outputfilename);//case ntpl file in input list
     277      HEPEVTConverter converter(inputFileList,outputfilename,DET->PDGtable);//case ntpl file in input list
    275278    }
    276279  else if(line.length() == 1+line.find_last_of(".hepmc"))
    … …  
    278281      cout <<"**                HepMC ASCII file format detected                 **"<<endl;
    279282      cout <<"**                  This can take several minutes                  **"<< endl;
    280       HepMCConverter converter(inputFileList,outputfilename);
     283      HepMCConverter converter(inputFileList,outputfilename,DET->PDGtable);
    281284    }
    282285  else {
    … …  
    309312  ExRootTreeBranch *branchZDC = treeWriter->NewBranch("ZDChits", TRootZdcHits::Class());
    310313  ExRootTreeBranch *branchRP220 = treeWriter->NewBranch("RP220hits", TRootRomanPotHits::Class());
    311   ExRootTreeBranch *branchFP420 = treeWriter->NewBranch("FP420hits", TRootForwardTaggerHits::Class());
     314  //ExRootTreeBranch *branchFP420 = treeWriter->NewBranch("FP420hits", TRootForwardTaggerHits::Class());
     315  ExRootTreeBranch *branchFP420 = treeWriter->NewBranch("FP420hits", TRootRomanPotHits::Class());
    312316 
    313317  TRootETmis *elementEtmis;
    … …  
    397401        {
    398402          TRootGenParticle *particle = new TRootGenParticle(particleG);
     403          PdgParticle pdg_part = DET->PDGtable[particle->PID];
     404          particle->Charge  = pdg_part.charge();
     405          particle->M  = pdg_part.mass();
     406          //particle->Charge=ChargeVal(particle->PID);
     407          particle->setFractions(); // init
    399408          int pid = abs(particle->PID);
    400           particle->Charge=ChargeVal(particle->PID);
    401           particle->setFractions(); // init
    402409         
    403410         

  • trunk/Examples/Convertors_Only.cpp

    r266 r380  
    4545#include "DataConverter.h"
    4646#include "HEPEVTConverter.h"
     47#include "HepMCConverter.h"
    4748#include "LHEFConverter.h"
    4849#include "STDHEPConverter.h"
    … …  
    6061
    6162using namespace std;
    62 
    63 //------------------------------------------------------------------------------
    64 void todo(string filename) {
    65   ifstream infile(filename.c_str());
    66   cout << "** TODO list ..." << endl;
    67   while(infile.good()) {
    68     string temp;
    69     getline(infile,temp);
    70     cout << "*" << temp << endl;
    71   }
    72   cout << "** done...\n";
    73 }
    7463
    7564//------------------------------------------------------------------------------
    … …  
    8877  delete [] appOpt;
    8978
    90   if(argc != 3) {
    91     cout << " Usage: " << argv[0] << " input_file output_file " << endl;
     79  if(argc != 4) {
     80    cout << " Usage: " << argv[0] << " input_file output_file PDG_Table" << endl;
    9281    cout << " input_list - list of files in Ntpl, StdHep of LHEF format," << endl;
    93     cout << " output_file - output file." << endl;
     82    cout << " output_file - output file," << endl;
     83    cout << " PDG_Table - file with PDG particle table" << endl;
    9484    exit(1);
    9585  }
    … …  
    160150 
    161151  //read the output TROOT file
    162   string inputFileList(argv[1]), outputfilename(argv[2]);
     152  string inputFileList(argv[1]), outputfilename(argv[2]), pdgtablefilename(argv[3]);
    163153  if(outputfilename.find(".root") > outputfilename.length()) {
    164154    cout <<"**    ERROR: 'output_file' should be a .root file. Exiting...      **"<< endl;
    165155    exit(1);
    166156  }
    167  
     157 
     158  RESOLution DET;
     159  DET.PdgTableFilename = pdgtablefilename;
     160  DET.ReadParticleDataGroupTable(); //  DET.PDGtable.print();
     161
    168162  TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE"); // Creates the file, but should be closed just after
    169163  outputFile->Close();
    … …  
    174168 
    175169  // data converters
    176   DataConverter *converter=NULL;
    177170  cout <<"**                                                                 **"<<endl;
    178171  cout <<"**        ####### Start convertion to TRoot format ########        **"<< endl;
    179172
    180   if(strstr(line.c_str(),".hep"))
    181     {                           
     173  if(line.length() == 1+line.find_last_of(".hep"))
     174    {
    182175      cout <<"**                 StdHEP file format detected                     **"<<endl;
    183176      cout <<"**                This can take several minutes                    **"<< endl;
    184       converter = new STDHEPConverter(inputFileList,outputfilename);//case ntpl file in input list
    185     }
    186   else if(strstr(line.c_str(),".lhe"))
     177      STDHEPConverter converter(inputFileList,outputfilename,DET.PDGtable);//case ntpl file in input list
     178    }
     179  else if(line.length() == 1+line.find_last_of(".lhe"))
    187180    {
    188181      cout <<"**                   LHEF file format detected                     **"<<endl;
    189182      cout <<"**                 This can take several minutes                   **"<< endl;
    190       converter = new LHEFConverter(inputFileList,outputfilename);//case ntpl file in input list
    191     }
    192   else if(strstr(line.c_str(),".root"))
     183      LHEFConverter converter(inputFileList,outputfilename,DET.PDGtable);//case ntpl file in input list
     184    }
     185  else if(line.length() == 1+line.find_last_of(".root"))
    193186    {
    194187      cout <<"**                   h2root file format detected                   **"<<endl;
    195188      cout <<"**                  This can take several minutes                  **"<< endl;
    196       converter = new HEPEVTConverter(inputFileList,outputfilename);//case ntpl file in input list
     189      HEPEVTConverter converter(inputFileList,outputfilename,DET.PDGtable);//case ntpl file in input list
     190    }
     191  else if(line.length() == 1+line.find_last_of(".hepmc"))
     192    {
     193      cout <<"**                HepMC ASCII file format detected                 **"<<endl;
     194      cout <<"**                  This can take several minutes                  **"<< endl;
     195      HepMCConverter converter(inputFileList,outputfilename,DET.PDGtable);
    197196    }
    198197  else {
    199       cerr << left  << setw(4) <<"**  "<<""
    200            << left  << setw(63) << line.c_str() <<""
    201            << right << setw(2) <<"**"<<endl;
    202       cerr <<"**         ERROR: File format not identified --  Exiting...        **"<< endl;
    203       cout <<"**                                                                 **"<< endl;
    204       cout <<"*********************************************************************"<< endl;
    205       return -1;};
     198    cerr << left  << setw(4) <<"**  "<<""
     199         << left  << setw(63) << line.c_str() <<""
     200         << right << setw(2) <<"**"<<endl;
     201    cerr <<"**         ERROR: File format not identified --  Exiting...        **"<< endl;
     202    cout <<"**                                                                 **"<< endl;
     203    cout <<"*********************************************************************"<< endl;
     204    return -1;};
     205
    206206  cout <<"**                       Exiting conversion...                     **"<< endl;
    207207  cout <<"*********************************************************************"<< endl;
    … …  
    209209
    210210
    211   delete converter;
    212211 
    213212}

  • trunk/data/DetectorCard.dat

    r374 r380  
    120120NEvents_Frog      100
    121121
     122# input PDG tables
     123PdgTableFilename  data/particle.tbl     // table with particle pid,mass,charge,...
     124

  • trunk/data/DetectorCard_ATLAS.dat

    r374 r380  
    120120NEvents_Frog      100
    121121
     122# input PDG tables
     123PdgTableFilename  data/particle.tbl     // table with particle pid,mass,charge,...
     124

  • trunk/data/DetectorCard_CMS.dat

    r374 r380  
    121121NEvents_Frog      100
    122122
     123# input PDG tables
     124PdgTableFilename  data/particle.tbl     // table with particle pid,mass,charge,...
     125

  • trunk/interface/DataConverter.h

    r357 r380  
    3838#include <string>
    3939
     40#include "PdgParticle.h"
    4041#include "Utilities/ExRootAnalysis/interface/ExRootTreeBranch.h"
    4142#include "Utilities/ExRootAnalysis/interface/LHEF.h"
    … …  
    4849
    4950  DataConverter() : Nevt(-1) {};
    50   DataConverter(const int &N) : Nevt(N) {};
     51  DataConverter(const PdgTable& table, const int &N) : Nevt(N), pdg_table(table) {};
    5152  virtual ~DataConverter() {};
    5253
    5354private:
    5455  int Nevt; // number of events to read; -1 means "all of them"
     56
     57public:
     58  PdgTable pdg_table;   // table of particles from PDG
    5559};
    5660

  • trunk/interface/HEPEVTConverter.h

    r357 r380  
    8888public:
    8989
    90   HEPEVTConverter(const string& inputFileList, const string& outputFileName, const int& Nevents=-1);
     90  HEPEVTConverter(const string& inputFileList, const string& outputFileName, const PdgTable& pdg, const int& Nevents=-1);
    9191
    9292};

  • trunk/interface/HepMCConverter.h

    r368 r380  
    4646class HepMCConverter : public DataConverter {
    4747  public:
    48     HepMCConverter(const string& inputFileList, const string& outputFileName, const int& Nevents=-1);
     48    HepMCConverter(const string& inputFileList, const string& outputFileName, const PdgTable& pdg, const int& Nevents=-1);
    4949    virtual ~HepMCConverter();
    5050

  • trunk/interface/LHEFConverter.h

    r357 r380  
    4949public:
    5050
    51   LHEFConverter(const string& inputFileList, const string& outputFileName, const int& Nevents=-1);
     51  LHEFConverter(const string& inputFileList, const string& outputFileName, const PdgTable& pdg, const int& Nevents=-1);
    5252  virtual ~LHEFConverter();
    5353

  • trunk/interface/STDHEPConverter.h

    r357 r380  
    4141class STDHEPConverter : public DataConverter {
    4242  public:
    43     STDHEPConverter(const string& inputFileList, const string& outputFileName, const int& Nevents=-1);
     43    STDHEPConverter(const string& inputFileList, const string& outputFileName, const PdgTable& pdg, const int& Nevents=-1);
    4444    virtual ~STDHEPConverter();
    4545

  • trunk/interface/SmearUtil.h

    r374 r380  
    4545#include "TSimpleArray.h"
    4646#include "PhysicsTower.hh"
     47#include "PdgParticle.h"
    4748
    4849using namespace std;
    … …  
    213214 
    214215 int NEvents_Frog;
    215   float PT_QUARKS_MIN;    // minimal pt needed for quarks to reach the tracker, in GeV
     216 float PT_QUARKS_MIN;    // minimal pt needed for quarks to reach the tracker, in GeV
     217
     218 string PdgTableFilename;
     219 //map<int,PdgParticle> PdgTable;
     220 PdgTable PDGtable;
    216221
    217222  // to sort a vector
    … …  
    221226  /// Reads the data card for the initialisation of the parameters
    222227  void ReadDataCard(const string datacard);
     228
     229  /// Reads the PDG table
     230  void ReadParticleDataGroupTable();
    223231 
    224232  /// Create the output log file

  • trunk/src/BFieldProp.cc

    r294 r380  
    3131
    3232#include "BFieldProp.h"
     33#include "PdgParticle.h"
    3334#include "SystemOfUnits.h"
    3435#include "PhysicalConstants.h"
    … …  
    131132  if (!DET->FLAG_bfield ) return;
    132133
    133   q  = ChargeVal(Part->PID) *eplus; // in units of 'e'
     134  double M;
     135  //int q1  = ChargeVal(Part->PID) *eplus; // in units of 'e'
     136  if(Part->M < -999) { // unitialised!
     137     PdgParticle pdg_part = DET->PDGtable[Part->PID];
     138     q  = pdg_part.charge() *eplus; // in units of 'e'
     139     M  = pdg_part.mass(); // GeV
     140  } else  { q = Part->Charge; M = Part->M; }
     141
    134142  if(q==0) return;
    135143  if(R_max==0) { cout << "ERROR: magnetic field has no lateral extention\n"; return;}
    … …  
    165173    double PT = Part->PT;
    166174    double E  = Part->E;              // [GeV]
    167     double M  = Part->M;              // [GeV]/c²
     175    //double M  = Part->M;            // [GeV]/c²
    168176      double Phi = UNDEFINED;
    169177
    … …  
    477485  double p  = sqrt(pz*pz + pt*pt); //sqrt(px*px+py*py+pz*pz);
    478486 
    479   double M  = Part->M;
     487  //double M  = Part->M; // see above
    480488  double vx = px/M;
    481489  double vy = py/M;

  • trunk/src/HEPEVTConverter.cc

    r352 r380  
    169169
    170170
    171 HEPEVTConverter::HEPEVTConverter(const string& inputFileList, const string& outputFileName, const int& Nevents) : DataConverter(Nevents)
     171HEPEVTConverter::HEPEVTConverter(const string& inputFileList, const string& outputFileName, const PdgTable& pdg, const int& Nevents) : DataConverter(pdg,Nevents)
    172172{
    173173  string buffer;

  • trunk/src/HepMCConverter.cc

    r376 r380  
    3636#include "BlockClasses.h"
    3737
    38 #include "SmearUtil.h"
     38#include "PdgParticle.h"
    3939#include "ExRootTreeWriter.h"
    4040#include "ExRootTreeBranch.h"
    … …  
    197197      element->D1 = da1;
    198198      element->D2 = da2;
    199       element->Charge = ChargeVal(pid);
    200199
    201200      element->E = index_to_particle[n]->momentum().e();
    … …  
    203202      element->Py = index_to_particle[n]->momentum().py();
    204203      element->Pz = index_to_particle[n]->momentum().pz();
    205       element->M =  index_to_particle[n]->momentum().m();
    206 
    207       float PT = sqrt(pow(element->Px,2)+pow(element->Py,2));
    208       element->PT = PT;
     204
     205
     206      //cout << "element->PID = " << pid << "\t";
     207      //PdgParticle pdg_part(PdgID[pid]);
     208      //element->M =  index_to_particle[n]->momentum().m(); // this is the particle virtuality, not its rest mass
     209      //element->M      = pdg_part.mass();
     210      //element->Charge = pdg_part.charge();
     211      //cout << "element->M = " << element->M << " \t element->Charge = " << element->Charge << endl;
     212
     213      element->PT = sqrt(pow(element->Px,2)+pow(element->Py,2));
    209214
    210215      momentum.SetPxPyPzE(element->Px, element->Py, element->Pz, element->E);
    … …  
    256261//------------------------------------------------------------------------------
    257262
    258 HepMCConverter::HepMCConverter(const string& inputFileList, const string& outputFileName, const int& Nevents) : DataConverter(Nevents)
     263HepMCConverter::HepMCConverter(const string& inputFileList, const string& outputFileName, const PdgTable& pdg, const int& Nevents) : DataConverter(pdg,Nevents)
    259264{
    260265 

  • trunk/src/LHEFConverter.cc

    r369 r380  
    3939#include "BlockClasses.h"
    4040#include "LHEFConverter.h"
     41#include "SmearUtil.h"
     42#include "PdgParticle.h"
    4143
    4244using namespace std;
    … …  
    8587    element->Pz = hepeup.PUP[particle][2];
    8688    element->E = hepeup.PUP[particle][3];
    87     element->M = hepeup.PUP[particle][4];
    88 
     89    //element->M = hepeup.PUP[particle][4];
     90    //element->Charge = ChargeVal(element->PID);
     91   
     92    PdgParticle pdg_part = pdg_table[element->PID];
     93    element->Charge = pdg_part.charge();
     94    element->M = pdg_part.mass();
     95   
    8996    momentum.SetPxPyPzE(element->Px, element->Py, element->Pz, element->E);
    9097    element->PT = momentum.Perp();
    … …  
    105112//------------------------------------------------------------------------------
    106113//Nevents not yet implemented! 08.03.2009
    107 LHEFConverter::LHEFConverter(const string& inputFileList, const string& outputFileName, const int& Nevents) :
    108  DataConverter(Nevents) {
     114LHEFConverter::LHEFConverter(const string& inputFileList, const string& outputFileName, const PdgTable& pdg, const int& Nevents) :
     115 DataConverter(pdg,Nevents) {
    109116
    110117  ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputFileName, "GEN");

  • trunk/src/STDHEPConverter.cc

    r376 r380  
    4141#include "stdhep_declarations.h"
    4242#include "STDHEPConverter.h"
     43#include "PdgParticle.h"
    4344
    4445using namespace std;
    … …  
    7980    element->D1 = myhepevt.jdahep[number][0] - 1;
    8081    element->D2 = myhepevt.jdahep[number][1] - 1;
    81 //    element->Charge = myhepevt.hepchg(element->PID);
    82     element->Charge = ChargeVal(element->PID);
     82    PdgParticle pdg_part = pdg_table[element->PID];
     83    element->Charge = pdg_part.charge();
     84    element->M = pdg_part.mass();
    8385
    8486    element->E = myhepevt.phep[number][3];
    … …  
    8688    element->Py = myhepevt.phep[number][1];
    8789    element->Pz = myhepevt.phep[number][2];
    88     element->M = myhepevt.phep[number][4];
     90    //element->M = myhepevt.phep[number][4];
    8991
    9092    momentum.SetPxPyPzE(element->Px, element->Py, element->Pz, element->E);
    … …  
    111113//------------------------------------------------------------------------------
    112114
    113 STDHEPConverter::STDHEPConverter(const string& inputFileList, const string& outputFileName, const int& Nevents) : DataConverter(Nevents)
     115STDHEPConverter::STDHEPConverter(const string& inputFileList, const string& outputFileName, const PdgTable& pdg, const int& Nevents) : DataConverter(pdg, Nevents)
    114116{
    115117  int ierr, entryType;

  • trunk/src/SmearUtil.cc

    r376 r380  
    4242#include <sstream>
    4343#include <iomanip>
     44#include <map>
    4445using namespace std;
    4546
    … …  
    195196  RP_cross_y        = 0.0;
    196197  RP_cross_ang      = 142.5;
     198
     199  PdgTableFilename  = "data/particle.tbl";
    197200 
    198201}
    … …  
    332335
    333336  PT_QUARKS_MIN    = DET.PT_QUARKS_MIN;
     337  PdgTableFilename = DET.PdgTableFilename;
     338  PDGtable         = DET.PDGtable;
    334339}
    335340
    … …  
    468473
    469474  PT_QUARKS_MIN    = DET.PT_QUARKS_MIN;
     475
     476  PdgTableFilename = DET.PdgTableFilename;
     477  PDGtable         = DET.PDGtable;
    470478  return *this;
    471479}
    … …  
    577585    else if(strstr(temp_string.c_str(),"FLAG_lhco"))        {curstring >> varname >> ivalue; FLAG_lhco        = ivalue;}
    578586    else if(strstr(temp_string.c_str(),"NEvents_Frog"))     {curstring >> varname >> ivalue; NEvents_Frog     = ivalue;}
     587
     588    else if(strstr(temp_string.c_str(),"PdgTableFilename"))     {curstring >> varname >> svalue; PdgTableFilename = svalue;}
    579589  }
    580590 
    … …  
    657667  f_out <<"**********************************************************************"<< endl;
    658668  f_out <<"**                                                                  **"<< endl;
    659   f_out<<"#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"<<"\n";
     669  f_out<<"#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"<<"\n";
     670  f_out<<"#  Input PDG table : " << PdgTableFilename << "                               *"<<"\n";
     671  f_out<<"*                                                                    *"<<"\n";
    660672  f_out<<"*                                                                    *"<<"\n";
    661673  f_out<<"#********************************                                    *"<<"\n";
    … …  
    12801292int ChargeVal(const int pid)
    12811293{
     1294  cout << "ChargeVal :: deprecated function, do not use it anymore" << endl;
    12821295  int charge;
    12831296  if(
    … …  
    12951308           
    12961309}
     1310
     1311//------------------------------------------------------------------------------
     1312void RESOLution::ReadParticleDataGroupTable() {
     1313
     1314  string temp_string;
     1315  istringstream curstring;
     1316
     1317  ifstream fichier_a_lire(PdgTableFilename.c_str());
     1318  if(!fichier_a_lire.good()) {
     1319    cout <<"**        ERROR: PDG Table ("<< PdgTableFilename
     1320         <<  ") not found! exit.                        **" << endl;
     1321    exit(1);
     1322    return;
     1323  }
     1324  // first three lines of the file are useless
     1325  getline(fichier_a_lire,temp_string);
     1326  getline(fichier_a_lire,temp_string);
     1327  getline(fichier_a_lire,temp_string);
     1328
     1329
     1330  while (getline(fichier_a_lire,temp_string)) {
     1331    curstring.clear(); // needed when using several times istringstream::str(string)
     1332    curstring.str(temp_string);
     1333    int ID; string name; int charge; float mass; float width; float lifetime;
     1334    // ID name   chg       mass    total width   lifetime
     1335    //  1 d      -1      0.33000     0.00000   0.00000E+00
     1336    curstring >> ID >> name >> charge >> mass >> width >> lifetime;
     1337    PdgParticle particle(ID,name,mass,charge/3.,width,lifetime);
     1338    PDGtable.insert(ID,particle);
     1339    //PdgTable.insert(pair<int,PdgParticle>(ID,particle));
     1340    //cout << PDGtable[ID].name() <<  "\t" << PDGtable[ID].mass() << "\t" << PDGtable[ID].charge() << endl;
     1341  }
     1342 
     1343} // ReadParticleDataGroupTable
Note: See TracChangeset for help on using the changeset viewer.