ECalCalibration/script/writeCALOROCToTxt.C

0001 #include <iostream>
0002 #include "TMath.h"
0003 #include "TTreeReader.h"
0004 #include "TTreeReaderArray.h"
0005 #include "TChain.h"
0006 #include "TH2.h"
0007 #include "TVector3.h"
0008 #include <fstream>
0009 #include <string>
0010 #include <fstream>
0011 #include <unordered_map>
0012
0013 void writeCALOROCToTxt(const std::string& input,
0014                        const std::string& outputTxt) {
0015     TChain chain("events");
0016     chain.Add(input.c_str());
0017
0018     TTreeReader reader(&chain);
0019
0020     TTreeReaderArray<int> genStatus(reader, "MCParticles.generatorStatus");
0021     TTreeReaderArray<double> momX(reader, "MCParticles.momentum.x");
0022     TTreeReaderArray<double> momY(reader, "MCParticles.momentum.y");
0023     TTreeReaderArray<double> momZ(reader, "MCParticles.momentum.z");
0024     TTreeReaderArray<int> PDG(reader, "MCParticles.PDG");
0025     TTreeReaderArray<double> mass(reader, "MCParticles.mass");
0026
0027     TTreeReaderArray<unsigned long> CellID(reader, "EcalBarrelScFiRecHits.cellID");
0028     TTreeReaderArray<float> CALOROCE(reader, "EcalBarrelScFiRecHits.energy");
0029     TTreeReaderArray<float> CALOROCX(reader, "EcalBarrelScFiRecHits.position.x");
0030     TTreeReaderArray<float> CALOROCY(reader, "EcalBarrelScFiRecHits.position.y");
0031     TTreeReaderArray<float> CALOROCZ(reader, "EcalBarrelScFiRecHits.position.z");
0032
0033     TTreeReaderArray<unsigned long> ScfiCellID(reader, "EcalBarrelScFiPNpeHits.cellID");
0034     TTreeReaderArray<unsigned int> ScfiBegin(reader, "EcalBarrelScFiPNpeHits.contributions_begin");
0035     TTreeReaderArray<unsigned int> ScfiEnd(reader, "EcalBarrelScFiPNpeHits.contributions_end");
0036     TTreeReaderArray<int> ScfiMCcId(reader, "_EcalBarrelScFiPNpeHits_contributions.index");
0037     TTreeReaderArray<int> ScfiMCId(reader, "_EcalBarrelScFiPAttenuatedHitContributions_particle.index");
0038
0039
0040     std::ofstream output(outputTxt.c_str());
0041     std::unordered_map<int, double> dp;
0042
0043     int nevent = 0;
0044     while(reader.Next()) {
0045         if(nevent % 10 == 0) std::cout << "Processing event " << nevent << std::endl;
0046         bool emptyEvent = true;
0047
0048         std::vector<int> Att2MC(ScfiCellID.GetSize(), -1);
0049         std::unordered_map<unsigned long, int> cellID2AttId;
0050
0051         for(size_t i = 0; i < ScfiCellID.GetSize(); ++i){
0052             for(size_t j = ScfiBegin[i]; j < ScfiEnd[i]; ++j){
0053                 int cid = ScfiMCcId[j];
0054                 Att2MC[i] = ScfiMCId[cid]; // BIG ASSUMPATIONS: ScfiMCId[j] all returns the same id for ScfiBegin[i] <= j < ScfiEnd[i]
0055             }
0056             cellID2AttId[ScfiCellID[i]] = i;
0057         }
0058
0059
0060         for(size_t i = 0; i < CellID.GetSize(); ++i) {
0061             auto it = cellID2AttId.find(CellID[i]);
0062             if(it == cellID2AttId.end()) continue;
0063
0064             int attId = it -> second;
0065             int mcId = Att2MC[attId];
0066             output << CALOROCX[i] << " "
0067                    << CALOROCY[i] << " "
0068                    << CALOROCZ[i] << " "
0069                    << CALOROCE[i] << " "
0070                    << CALOROCZ[i] << " "
0071                    << genStatus[mcId] << " "
0072                    << momX[mcId] << " "
0073                    << momY[mcId] << " "
0074                    << momZ[mcId] << " "
0075                    << PDG[mcId] << " "
0076                    << mcId << " "
0077                    << CALOROCE[i]
0078                    << " " << nevent << std::endl;
0079         }
0080         for(int i = 0; i < 12; ++i) output << "-9999 ";
0081         output << nevent << std::endl; // this serves as line break for next event
0082         ++nevent;
0083     }
0084 }