ECalCalibration/script/AttCor.C

0001 #include <cmath>
0002 #include <fstream>
0003 #include <iostream>
0004 #include "TMath.h"
0005 #include "TDatabasePDG.h"
0006 #include "TVector3.h"
0007 #include <unordered_map>
0008
0009 void AttCor(const std::string& side,
0010             const std::string& inputTxt,
0011             const std::string& outputTxt,
0012             double thetaMin = 0,
0013             double thetaMax = 180) {
0014     std::ifstream input(inputTxt.c_str());
0015     std::ofstream output(outputTxt.c_str());
0016
0017     // read input
0018     double x, y, temp, npe, z, momx, momy, momz;
0019     int status, mcId, pdg, nevent;
0020
0021     // convert pdg to mass
0022     auto pdgDB = TDatabasePDG::Instance();
0023
0024     // sum energies
0025     std::unordered_map<int, double> mcId2TotE, mcId2Pdg, mcId2Mom, mcIdX, mcIdY, mcIdZ;
0026     std::unordered_map<int, int> mcStatus;
0027     double prevMomx = 0, prevMomy = 0, prevMomz = 0;
0028     while((input >> x >> y >> temp >> npe >> z >> status >> momx >> momy >> momz >> pdg >> mcId >> temp >> nevent)) {
0029          if(temp < -999) {
0030               if(nevent % 100 == 0) std::cout << "Processing event " << nevent << std::endl;
0031               // end of event line
0032               // summarized the event
0033               for(const auto& key : mcId2TotE) {
0034                   TVector3 vec(mcIdX[key.first], mcIdY[key.first], mcIdZ[key.first]);
0035                   if(thetaMin <= vec.Theta()*TMath::RadToDeg() && vec.Theta()*TMath::RadToDeg() <= thetaMax) {
0036                       auto particle = pdgDB -> GetParticle(mcId2Pdg[key.first]);
0037                       double mass = 0;
0038                       if(particle) mass = particle -> Mass(); // in GeV/c^2
0039                       double totE = std::sqrt(mcId2Mom[key.first]*mcId2Mom[key.first] + mass*mass);
0040                       output << key.first << " "
0041                              << mcId2Pdg[key.first] << " "
0042                              << vec.Theta()*TMath::RadToDeg() << " "
0043                              << key.second << " "
0044                              << totE << " "
0045                              << vec.Phi() * TMath::RadToDeg() << " "
0046                              << vec.Eta() << " "
0047                              << vec.Phi() * TMath::RadToDeg() << " "
0048                              << vec.Eta() << " "
0049                              << mcStatus[key.first] << " "
0050                              << nevent << std::endl;
0051                   }
0052               }
0053               output << "-9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999 -9999 " << nevent << std::endl;
0054               mcId2TotE.clear();
0055               mcId2Pdg.clear();
0056               mcIdX.clear();
0057               mcIdY.clear();
0058               mcIdZ.clear();
0059           }
0060           else {
0061               mcIdX[mcId] += x;
0062               mcIdY[mcId] += y;
0063               mcIdZ[mcId] += z;
0064               mcId2Mom[mcId] = TVector3(momx, momy, momz).Mag();
0065               mcId2TotE[mcId] += npe;
0066               mcStatus[mcId] = status;
0067               mcId2Pdg[mcId] = pdg;
0068           }
0069      }
0070
0071 }