STARLIGHT/starlight/utils/.svn/text-base/convertStarlightAsciiToTree.C.svn-base

   1 #include <iostream>
   2 #include <fstream>
   3 #include <sstream>
   4
   5 #include "TLorentzVector.h"
   6 #include "TClonesArray.h"
   7 #include "TTree.h"
   8 #include "TFile.h"
   9
  10
  11 using namespace std;
  12
  13
  14 void readStarlightAscii(const char* inFileName  = "slight.out",
  15                         const char* outFileName = "starlight.root")
  16 {
  17         const double daughterMass = 0.13957018; // charged pion [GeV]
  18
  19         // create output tree
  20         TFile* outFile = new TFile(outFileName, "RECREATE");
  21         if (!outFile) {
  22                 cerr << "    error: could not create output file '" << outFileName << "'" << endl;
  23                 return;
  24         }
  25         TTree*          outTree           = new TTree("starlightTree", "starlightTree");
  26         TLorentzVector* parentParticle    = new TLorentzVector();
  27         TClonesArray*   daughterParticles = new TClonesArray("TLorentzVector");
  28         outTree->Branch("parent",    "TLorentzVector", &parentParticle,    32000, -1);
  29         outTree->Branch("daughters", "TClonesArray",   &daughterParticles, 32000, -1);
  30
  31         ifstream inFile;
  32         inFile.open(inFileName);
  33         unsigned int countLines = 0;
  34         while (inFile.good()) {
  35                 string       line;
  36                 stringstream lineStream;
  37
  38                 // read EVENT
  39                 string label;
  40                 int    eventNmb, nmbTracks;
  41                 if (!getline(inFile, line))
  42                         break;
  43                 ++countLines;
  44                 lineStream.str(line);
  45                 assert(lineStream >> label >> eventNmb >> nmbTracks);
  46                 if (!(label == "EVENT:"))
  47                         continue;
  48
  49                 // read vertex
  50                 if (!getline(inFile, line))
  51                         break;
  52                 ++countLines;
  53                 lineStream.str(line);
  54                 assert(lineStream >> label);
  55                 assert(label == "VERTEX:");
  56
  57                 *parentParticle = TLorentzVector(0, 0, 0, 0);
  58                 for (int i = 0; i < nmbTracks; ++i) {
  59                         // read tracks
  60                         int    particleCode;
  61                         double momentum[3];
  62                         if (!getline(inFile, line))
  63                                 break;
  64                         ++countLines;
  65                         lineStream.str(line);
  66                         assert(lineStream >> label >> particleCode >> momentum[0] >> momentum[1] >> momentum[2]);
  67                         assert(label == "TRACK:");
  68                         const double E = sqrt(  momentum[0] * momentum[0] + momentum[1] * momentum[1]
  69                                               + momentum[2] * momentum[2] + daughterMass * daughterMass);
  70                         new ( (*daughterParticles)[i] ) TLorentzVector(momentum[0], momentum[1], momentum[2], E);
  71                         *parentParticle += *(static_cast<TLorentzVector*>(daughterParticles->At(i)));
  72                 }
  73
  74                 daughterParticles->Compress();
  75                 outTree->Fill();
  76         }
  77
  78         outTree->Write();
  79         if (outFile) {
  80                 outFile->Close();
  81                 delete outFile;
  82         }
  83 }