SplitDcmsmm.C

Go to the documentation of this file.

// Number of tokens in the header part
const Int_t nHead = 20;
const Int_t ntokh[nHead] = {7, 10, 6, 6, 8, 10, 16, 16, 16, 13, 13, 13, 13, 16, 2, 8};
 
// Read file header
void ReadHeader(ifstream& fin, vector<TString>& header);
// Read DCMSMM event
Int_t ReadEvent(ifstream& fin, vector<TString>& evBuff;);
// Write DCM-SMM event with the assigned number iev
void WriteEvent(ofstream& fout, vector<TString>& evBuff;);
// Get a string representation of the ratio: nevent/divider
TString GetRatioString(Int_t nevent, Int_t divider);
 
//__________________________________________________________________________
// Read DCMSMM file f12 and split it into pieces
void SplitDcmsmm(TString inFile = "./DCMSMM_XeCsI_3.8_AGeV_mb_50k_2.r12", Int_t nbefore = 10000, Int_t nafter = 2500)
{
    // trying to get input file from bash environment
    TString fileIn = TString(gSystem->Getenv("FILEIN"));
    if (fileIn != "")
        inFile = fileIn;
 
    gSystem->ExpandPathName(inFile);
    // Check if path exists
    if (gSystem->AccessPathName(inFile.Data(), kFileExists) == true) {
        cerr << "Error: Path does not exist: " << inFile << endl;
        return;
    }
    // Get detailed path information
    FileStat_t stat;
    if (gSystem->GetPathInfo(inFile, stat) != 0) {
        cerr << "Error: Cannot get path info: " << inFile << endl;
        return;
    }
 
    // Collect files to be processed
    vector<TString> filesToProcess;
    // Check if it is a FILE or DIRECTORY
    if (stat.fMode & S_IFREG) {
        // === CASE 1: It's a FILE ===
        // Check extension
        if (inFile.EndsWith(".r12")) {
            filesToProcess.push_back(inFile);
        } else {
            cerr << "Warning: Input file does not have .r12 extension: " << inFile << endl;
            return;
        }
    } else if (stat.fMode & S_IFDIR) {
        // === CASE 2: It's a DIRECTORY ===
        void* dirp = gSystem->OpenDirectory(inFile);
        if (!dirp) {
            cerr << "Error: Cannot open directory: " << inFile << endl;
            return;
        }
 
        const char* p_file_name = nullptr;
        while ((p_file_name = gSystem->GetDirEntry(dirp)) != nullptr) {
            TString fileName(p_file_name);
            // Check extension
            if (fileName.EndsWith(".r12")) {
                // Build full path
                TString fullPath = TString(inFile) + "/" + fileName;
                filesToProcess.push_back(fullPath);
            }
        }
        gSystem->FreeDirectory(dirp);
    } else {
        cerr << "Error: Input path is neither a file nor a directory: " << inFile << endl;
        return;
    }
 
    const int nfiles = nbefore / nafter;   // number of output files
    TString event_before_k = TString::Format("%sk", GetRatioString(nbefore, 1000).Data());
    TString event_after_k = TString::Format("%sk", GetRatioString(nafter, 1000).Data());
 
    // Process all input files
    cout << "Found " << filesToProcess.size() << " file(s) to process:" << endl << endl;
    for (const auto& input_file : filesToProcess) {
        cout << "Current input file: " << input_file << endl;
        ifstream fin;
        fin.open(input_file.Data());
 
        // read header of the current input file
        vector<TString> header ReadHeader(fin, header);
 
        int success = 0, ntot = 0;
        // Loop over file pieces
        for (int ifile = 0; ifile < nfiles; ++ifile) {
            // Change file name
            TString outFile = input_file;
            // search for the first position of the string containing number of events
            // int ip = outFile.First(event_before_k.Data());
            // search for the last position of the string containing number of events
            Int_t lastPos = -1, pos = 0;
            while ((pos = outFile.Index(event_before_k.Data(), pos, TString::kIgnoreCase)) >= 0) {
                lastPos = pos;
                pos += 1;   // Move forward to find next occurrence
            }
            if (lastPos >= 0)
                outFile.Replace(lastPos, event_before_k.Length(), event_after_k);
 
            // change number of the file according to the nfiles (number of output files for one input file)
            int idot = input_file.Last('.');
            int iund = input_file.Last('_');
            int last_number = (TString(input_file(iund + 1, idot - iund - 1))).Atoi();
            int new_last_number = (last_number - 1) * nfiles + ifile + 1;
            TString str_new_last_number = "";
            str_new_last_number += new_last_number;
            outFile.Replace(iund + 2, TString::Itoa(last_number, 10).Length(), str_new_last_number.Data());
            cout << " Current output file: " << outFile << endl;
 
            // Open output file
            ofstream fout;
            fout.open(outFile.Data());
            for (int j = 0; j < nHead; ++j)
                fout << header[j] << endl;
            // if (ifile == 0) for (int j = 0; j < nHead; ++j) fout << header[j] << endl;
            Int_t nsel = 0;
 
            for (Int_t iev = 0; iev < nbefore * 10; ++iev) {
                vector<TString> evBuff;
                success = ReadEvent(fin, evBuff);
                if (success == 99)
                    break;   // EOF
                if (success == 0) {
                    TString evNew;
                    ++nsel;
                    evNew.Form("%9d", nsel);
                    evBuff[0].Replace(0, 9, evNew);
                    WriteEvent(fout, evBuff);
                }
                if (iev % 1000 == 0)
                    cout << iev << " " << nsel << " " << success << endl;
                if (nsel >= nafter)
                    break;   // number of selected events reached
            }
            cout << " Number of selected events: " << nsel << endl;
            ntot += nsel;
            fout.close();
            // break;
            if (success == 99)
                break;
        }   // for (int ifile = 0;
 
        fin.close();
        cout << " Total number of events: " << ntot << endl;
    }   // cycle: process all input files
}
 
//__________________________________________________________________________
// Read file header
void ReadHeader(ifstream& fin, vector<TString>& header)
{
    TString chline;
    for (Int_t ihead = 0; ihead < nHead; ++ihead) {
        chline.ReadLine(fin, kFALSE);
        header.push_back(chline);
    }
}
 
//__________________________________________________________________________
// Read DCMSMM event
Int_t ReadEvent(ifstream& fin, vector<TString>& evBuff;)
{
    TString chline;
    TObjArray* tokens = NULL;
    TObjString* objstr = NULL;
    Int_t success = 1, nPart = 0, nNucl = 0, iok = 0, iev = 0;
    Float_t b, bx, by, beamA, beamZ;
    // success = 0; // check
 
    // for (Int_t ihead = 0; ihead < nHead; ++ihead) {
    for (Int_t ihead = 0; ihead > -1; ++ihead) {
        chline.ReadLine(fin, kFALSE);
        if (fin.eof())
            return 99;   // EOF
        tokens = chline.Tokenize(" -");
        Int_t ntok = tokens->GetEntriesFast();
        if (ntok <= 5 && ntok >= 3) {
            // Event header
            // cout << chline << endl;
            // objstr = (TObjString*) tokens->UncheckedAt(2);
            // if (!objstr->GetString().Contains("P")) {
            objstr = (TObjString*)tokens->UncheckedAt(0);
            TString token = objstr->GetString();
            if (token.IsDigit()) {
                sscanf(chline.Data(), "%d %d %f %f %f", &iev, &nPart, &b, &bx, &by);
                // Impact-parameter
                Double_t impPar = objstr->GetString().Atof();
                // bout << iev << " " << nPart << " " << b << " " << bx << " " << by << endl;
                ++iok;
            }
            evBuff.push_back(chline);
        } else if (ntok == 8) {
            // Beam remnants
            objstr = (TObjString*)tokens->UncheckedAt(0);
            TString token = objstr->GetString();
            if (token.IsDigit()) {
                sscanf(chline.Data(), "%d %f %f", &nNucl, &beamA, &beamZ);
                // tmpout << nNucl << " " << beamA << " " << beamZ << endl;
            }
            evBuff.push_back(chline);
            // Read fragments
 
            for (int j = 0; j < nNucl; ++j) {
                chline.ReadLine(fin, kFALSE);
                evBuff.push_back(chline);
            }
        } else if (ntok <= 1)
            evBuff.push_back(chline);   // line with number of particles
        if (tokens) {
            tokens->Delete();
            delete tokens;
            tokens = NULL;
        }
        if (ntok > 1)
            continue;
 
        // Read particles
        //*
        for (Int_t ipart = 0; ipart < nPart; ++ipart) {
            chline.ReadLine(fin, kFALSE);
            evBuff.push_back(chline);
        }
        //*/
 
        return 0;   // event processed
    }   // for (Int_t ihead = 0;
 
    return success;
}
 
//__________________________________________________________________________
// Write DCM-SMM event with the assigned number iev
void WriteEvent(ofstream& fout, vector<TString>& evBuff;)
{
    Int_t nlines = evBuff.size();
    for (int j = 0; j < nlines; ++j)
        fout << evBuff[j] << endl;
}
 
/**
 * @brief Get a string representation of the ratio: nevent/divider
 * If divider is zero, the result is "inf"
 * The ratio is formatted as string with integer (e.g. "21") or float (e.g. "2.5", without trailing zeros)
 * @param nevent The numerator of the ratio
 * @param divider The denominator of the ratio
 * @return A string representation of the ratio
 */
TString GetRatioString(Int_t nevent, Int_t divider)
{
    if (divider == 0)
        return "inf";
 
    Double_t ratio = ((double)nevent) / divider;
    // Format the ratio with a precision of 3
    TString str = Form("%.3f", ratio);
 
    // Remove trailing zeros after the decimal point
    str = str.Strip(TString::kTrailing, '0');
    // Guarantee that there is a decimal point, not a comma
    str.ReplaceAll(",", ".");
    // If there is a decimal point ath the end, remove it
    if (str.EndsWith("."))
        str.Remove(str.Length() - 1, 1);
 
    return str;
}