From 150359420637fbeb32ad74342f55ccd189ff3b5c Mon Sep 17 00:00:00 2001 From: bspage912 Date: Fri, 3 May 2024 17:24:32 -0400 Subject: [PATCH 01/10] Initial commit to branch may24LutGeneration which will be used for pfRICH LUT generation for the May campaign. Scripts to run simulation and generate analysis trees are added along with a README with basic instructions. Changes to be committed: new file: CondorSubmitter.csh new file: README.BP new file: hepmc_writer.C new file: runSimu.sh new file: scripts/reco-epic-LUT.C --- CondorSubmitter.csh | 47 +++++++++++++ README.BP | 19 ++++++ hepmc_writer.C | 121 +++++++++++++++++++++++++++++++++ runSimu.sh | 61 +++++++++++++++++ scripts/reco-epic-LUT.C | 147 ++++++++++++++++++++++++++++++++++++++++ 5 files changed, 395 insertions(+) create mode 100755 CondorSubmitter.csh create mode 100644 README.BP create mode 100644 hepmc_writer.C create mode 100755 runSimu.sh create mode 100644 scripts/reco-epic-LUT.C diff --git a/CondorSubmitter.csh b/CondorSubmitter.csh new file mode 100755 index 0000000..e361d37 --- /dev/null +++ b/CondorSubmitter.csh @@ -0,0 +1,47 @@ +#! /bin/env csh + +## Variabe Exec should be set to the absolute path of the runSimu.sh script (should be in same dir as this script) +## set Exec = /yourDIR/runSimu.sh +set Exec = + +####### Initialize condor file +echo "" > CondorFile +echo "Universe = vanilla" >> CondorFile +echo "Executable = ${Exec}" >> CondorFile +echo "getenv = true" >> CondorFile + +## Output Directory +## Output should point to the dir you wan the error and stdOut files to go (again, likely the same dir as this script) +set Output = + +####### Set number of jobs and events per job +set NUMBER = 1 +set LIMIT = 10 +set NEVENTS = 100000 + +####### Loop over jobs +while ( "$NUMBER" <= "$LIMIT" ) + + set OutFile = ${Output}/submitSimu_$NUMBER.out + set ErrFile = ${Output}/submitSimu_$NUMBER.err + + set Args = ( $NUMBER $NEVENTS ) + echo "" >> CondorFile + echo "Output = ${OutFile}" >> CondorFile + echo "Error = ${ErrFile}" >> CondorFile + echo "Arguments = ${Args}" >> CondorFile + echo "Queue" >> CondorFile + + echo Submitting: + echo $Exec $Args + echo "Logging output to " $OutFile + echo "Logging errors to " $ErrFile + + @ NUMBER++ + +end + + +#condor_submit CondorFile + + diff --git a/README.BP b/README.BP new file mode 100644 index 0000000..3ef4517 --- /dev/null +++ b/README.BP @@ -0,0 +1,19 @@ +Four scripts are needed to run mass amounts of simulation to create LUTs at RCF: + +hepmc_writer.C: This script produces single particle hepmc events over the perscribed kinematic ranges. Input to the GEANT simulator. Should not need to be modified + +scripts/reco-epic-LUT.C: This script reads in the output of the GEANT simulation and generates a root tree from which we can generate LUTs and diagnostic plots + +runSimu.sh: Master script which will generate hepmc files, run the pfRICH GEANT simulation, and generate the simple analysis root tree. User needs to modify to specify paths to the various scripts + +CondorSubmitter.csh: This script generates the condor files which are submitted to the batch system at RCF. The user will need to specify the path to runSimu.sh and the location for log files and also the number of jobs to submit and the number of events to simulate per job (NEVENTS specifies the number of events per species so setting NEVENTS = 1000 will generate 1000 electron events, 1000 pion, events, 1000 kaon events, and 1000 proton events. + +The below instructions assume you have a working pfRICH simulation environment at RCF. The output will go directly to the directory you submit from, so this needs to be set up in a location with ample storage (100s to 1000s of GBs). + +1. Modify CondorSubmitter.csh and runSimu.sh as described above + +2. Run CondorSubmitter.csh to generate a CondorFile + +3. Submit the CondorFile to the batch system with the command: condor_submit CondorFile + +4. You can check progress with the command: condor_q userName (Seems a job with NEVENTS = 100000 should take between 5.5 adn 6.5 hours \ No newline at end of file diff --git a/hepmc_writer.C b/hepmc_writer.C new file mode 100644 index 0000000..8dcb128 --- /dev/null +++ b/hepmc_writer.C @@ -0,0 +1,121 @@ +// +// export LD_LIBRARY_PATH=/home/eic/lib64:${LD_LIBRARY_PATH} +// +// root -l +// +// root [0] gSystem->AddIncludePath("-I/home/eic/include"); +// root [1] .x hepmc_writer_two_particles.C+("out.hepmc", 1000) +// + +#include "HepMC3/GenEvent.h" +#include "HepMC3/ReaderAscii.h" +#include "HepMC3/WriterAscii.h" +#include "HepMC3/Print.h" + +#include +#include +#include +#include + +#include +#include +#include + +using namespace HepMC3; + +void hepmc_writer(const char *out_fname, int part, Double_t thMin, Double_t thMax, Double_t pMin, Double_t pMax, int n_events) +{ + auto *DatabasePDG = new TDatabasePDG(); + auto *electron = DatabasePDG->GetParticle(11); + auto *pion = DatabasePDG->GetParticle(211); // + auto *kaon = DatabasePDG->GetParticle(321); + auto *proton = DatabasePDG->GetParticle(2212); + //const char * out_fname= "out.hepmc"; + WriterAscii hepmc_output(out_fname); + int events_parsed = 0; + GenEvent evt(Units::GEV, Units::MM); + + // Random number generator + //TRandom *rdmn_gen = new TRandom(0x12345678); + TRandom3 *rdmn_gen = new TRandom3(0); + + for (events_parsed = 0; events_parsed < n_events; events_parsed++) { + // type 4 is beam + GenParticlePtr p1 = + std::make_shared(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4); + GenParticlePtr p2 = std::make_shared( + FourVector(0.0, 0.0, 100.0, 100.004), 2212, 4); + + GenVertexPtr v1 = std::make_shared(); + v1->add_particle_in(p1); + v1->add_particle_in(p2); + + // Set Kinematics + //Double_t eta = rdmn_gen->Uniform(etaMin, etaMax); + //Double_t th = 2*std::atan(exp(-eta)); + Double_t th = rdmn_gen->Uniform(thMin, thMax); + Double_t p = rdmn_gen->Uniform(pMin, pMax); + Double_t phi = rdmn_gen->Uniform(0.0*M_PI/180.0, 360.0*M_PI/180.); + //Double_t phi = rdmn_gen->Uniform(85.0*M_PI/180.0, 95.0*M_PI/180.); + + Double_t px = p * std::cos(phi) * std::sin(th); + Double_t py = p * std::sin(phi) * std::sin(th); + Double_t pz = p * std::cos(th); + + /* + if(events_parsed < 50) + { + cout << events_parsed << " " << p << " " << eta << " " << th << " " << phi << " " << px << " " << py << " " << pz << endl; + } + */ + + TParticlePDG *particle; + switch(part) + { + case 0: + particle = electron; + break; + case 1: + particle = pion; + break; + case 2: + particle = kaon; + break; + case 3: + particle = proton; + break; + default: + cout << "Invalid Particle Selection - Default to Pion" << endl; + particle = pion; + break; + } + + //auto particle = pion; + //if(part == 0) particle = electron; + //if(part == 1) particle = pion; + //if(part == 2) particle = kaon; + //if(part == 3) particle = proton; + + GenParticlePtr pq = std::make_shared(FourVector( + px, py, pz, + sqrt(p*p + pow(particle->Mass(), 2))), + particle->PdgCode(), 1); + v1->add_particle_out(pq); + + evt.add_vertex(v1); + + if (events_parsed == 0) { + std::cout << "First event: " << std::endl; + Print::listing(evt); + } + + hepmc_output.write_event(evt); + if (events_parsed % 10000 == 0) { + std::cout << "Event: " << events_parsed << std::endl; + } + evt.clear(); + } + hepmc_output.close(); + std::cout << "Events parsed and written: " << events_parsed << std::endl; + exit(0); +} diff --git a/runSimu.sh b/runSimu.sh new file mode 100755 index 0000000..1a18562 --- /dev/null +++ b/runSimu.sh @@ -0,0 +1,61 @@ +#!/bin/bash + +RUN=$1 +NUM=$2 + +date + +fileNameElectron=electronEvents_${RUN}.hepMC +fileNamePion=pionEvents_${RUN}.hepMC +fileNameKaon=kaonEvents_${RUN}.hepMC +fileNameProton=protonEvents_${RUN}.hepMC + +echo "" +echo "Generating HepMC Files" +echo ${fileNameElectron} +echo ${fileNamePion} +echo ${fileNameKaon} +echo ${fileNameProton} + +## Supply absolute path to the hepmc_writer.C script +root -b -l -q "/yourDir/hepmc_writer.C+(\"${fileNameElectron}\",0,2.70,3.10,0.1,15.0,${NUM})" +root -b -l -q "/yourDir/hepmc_writer.C+(\"${fileNamePion}\",1,2.70,3.10,0.1,15.0,${NUM})" +root -b -l -q "/yourDir/hepmc_writer.C+(\"${fileNameKaon}\",2,2.70,3.10,0.1,15.0,${NUM})" +root -b -l -q "/yourDir/hepmc_writer.C+(\"${fileNameProton}\",3,2.70,3.10,0.1,15.0,${NUM})" + +simuNameElectron=pfRICH-epic_electron_${RUN}.root +simuNamePion=pfRICH-epic_pion_${RUN}.root +simuNameKaon=pfRICH-epic_kaon_${RUN}.root +simuNameProton=pfRICH-epic_proton_${RUN}.root + +echo "" +echo "Generating Simulation Files" +echo ${simuNameElectron} +echo ${simuNamePion} +echo ${simuNameKaon} +echo ${simuNameProton} + +## Supply absolute path to the pfRICH build dir +/yourDir/build/pfrich-epic -i ${fileNameElectron} -o ${simuNameElectron} -s ${NUM} +/yourDir/build/pfrich-epic -i ${fileNamePion} -o ${simuNamePion} -s ${NUM} +/yourDir/build/pfrich-epic -i ${fileNameKaon} -o ${simuNameKaon} -s ${NUM} +/yourDir/build/pfrich-epic -i ${fileNameProton} -o ${simuNameProton} -s ${NUM} + +rootNameElectron=outElectron_${RUN}.root +rootNamePion=outPion_${RUN}.root +rootNameKaon=outKaon_${RUN}.root +rootNameProton=outProton_${RUN}.root + +echo "" +echo "Generating Root Trees" +echo ${rootNameElectron} +echo ${rootNamePion} +echo ${rootNameKaon} +echo ${rootNameProton} + +root -b -l -q "scripts/reco-epic-LUT.C(\"${simuNameElectron}\",\"${rootNameElectron}\")" +root -b -l -q "scripts/reco-epic-LUT.C(\"${simuNamePion}\",\"${rootNamePion}\")" +root -b -l -q "scripts/reco-epic-LUT.C(\"${simuNameKaon}\",\"${rootNameKaon}\")" +root -b -l -q "scripts/reco-epic-LUT.C(\"${simuNameProton}\",\"${rootNameProton}\")" + +date diff --git a/scripts/reco-epic-LUT.C b/scripts/reco-epic-LUT.C new file mode 100644 index 0000000..4fd4917 --- /dev/null +++ b/scripts/reco-epic-LUT.C @@ -0,0 +1,147 @@ + + +void reco_epic_LUT(const char *dfname, const char *output, const char *cfname = 0) +{ + auto *reco = new ReconstructionFactory(dfname, cfname, "pfRICH"); + + //const char* output="test.hist.root"; + TFile *ofile = TFile::Open(output,"recreate"); + + TTree *tree = new TTree("T","Tree of Quantities"); + + Int_t partPDG; + Double_t partMom; + Double_t partEta; + Double_t partPhi; + Double_t partVertX; + Double_t partVertY; + Double_t partVertZ; + + Int_t recoPDG; + Int_t recoPhotons; + Int_t recoCherenkovHitCount; + Double_t recoTrackCherenkovAngle; + + tree->Branch("partPDG",&partPDG,"partPDG/I"); + tree->Branch("partMom",&partMom,"partMom/D"); + tree->Branch("partEta",&partEta,"partEta/D"); + tree->Branch("partPhi",&partPhi,"partPhi/D"); + tree->Branch("partVertX",&partVertX,"partVertX/D"); + tree->Branch("partVertY",&partVertY,"partVertY/D"); + tree->Branch("partVertZ",&partVertZ,"partVertZ/D"); + + tree->Branch("recoPDG",&recoPDG,"recoPDG/I"); + tree->Branch("recoPhotons",&recoPhotons,"recoPhotons/I"); + tree->Branch("recoCherenkovHitCount",&recoCherenkovHitCount,"recoCherenkovHitCount/I"); + tree->Branch("recoTrackCherenkovAngle",&recoTrackCherenkovAngle,"recoTrackCherenkovAngle/D"); + + // + // Factory configuration part; + // + //reco->IgnoreTimingInChiSquare(); + //reco->IgnorePoissonTermInChiSquare(); + //reco->SetSingleHitCCDFcut(0.005); + //reco->RemoveAmbiguousHits(); + // This only affects the calibration stage; + //reco->SetDefaultSinglePhotonThetaResolution(0.0040); + // Sensor active area pixelated will be pixellated NxN in digitization; + //reco->SetSensorActiveAreaPixellation(24); + // [rad] (should match SPE sigma) & [ns]; + //auto *a1 = reco->UseRadiator("Aerogel225", 0.0040); + + auto *a1 = reco->UseRadiator("BelleIIAerogel3"); + //auto *a1 = reco->UseRadiator("BelleIIAerogel3"); + + + //reco->SetSinglePhotonTimingResolution(0.030); + //reco->SetQuietMode(); + //reco->AddHypothesis(-11); + reco->AddHypothesis(11); + reco->AddHypothesis("pi+"); + reco->AddHypothesis(321); + reco->AddHypothesis(2212); + //reco->IgnoreMcTruthPhotonDirectionSeed(); + + // Mark all pads hit by "calibration" (above the QE curve) photons as "useless"; + reco->AddBlackoutRadiator("QuartzWindow"); + reco->AddBlackoutRadiator("Acrylic"); + reco->AddBlackoutRadiator("GasVolume"); + // Carelessly remove (0x1 << n)x(0x1 << n) square area "around" these hits; + reco->SetBlackoutBlowupValue(3); + + //auto hEta = new TH1D("hEta","",100,-5.,5.); + //auto hMom = new TH1D("hMom","",100,0.,10.); + + //auto factoryPhotons = new TH1D("factoryPhotons","",50,0.,50.); + + //auto hmatch = new TH1D("hmatch", "PID evaluation correctness", 2, 0, 2); + //auto hthtr1 = new TH1D("thtr1", "Cherenkov angle (track)", 200, 220, 320); + //auto hthtr1 = new TH1D("thtr1", "Cherenkov angle (track)", 40, 270, 290); + //auto hthtr1 = new TH1D("thtr1", "Cherenkov angle (track)", 580, 0., 290.); + // For a dual aerogel configuration; + //auto hthtr2 = new TH1D("thtr2", "Cherenkov angle (track)", 200, 220, 320); + + reco->UseAutomaticCalibration(); + // This call is mandatory; second argument: statistics (default: all events); + reco->PerformCalibration(200); + + //int numElec = 0; + //int numPion = 0; + //int numKaon = 0; + //int numProt = 0; + { + CherenkovEvent *event; + + // Event loop; + while((event = reco->GetNextEvent())) { + for(auto mcparticle: event->ChargedParticles()) { + if (!mcparticle->IsPrimary()) continue; + + /* + if (mcparticle->GetPDG() == mcparticle->GetRecoPdgCode()) { + hmatch->Fill(0.5); + } else { + hmatch->Fill(1.5); + } //if + */ + + // GetRecoCherenkovPhotonCount + + //hEta->Fill(mcparticle->GetVertexMomentum().PseudoRapidity()); + //hMom->Fill(mcparticle->GetVertexMomentum().Mag()); + + partPDG = mcparticle->GetPDG(); + partMom = mcparticle->GetVertexMomentum().Mag(); + partEta = mcparticle->GetVertexMomentum().PseudoRapidity(); + partPhi = mcparticle->GetVertexMomentum().Phi(); + partVertX = mcparticle->GetVertexPosition().X(); + partVertY = mcparticle->GetVertexPosition().Y(); + partVertZ = mcparticle->GetVertexPosition().Z(); + + recoPDG = mcparticle->GetRecoPdgCode(); + recoPhotons = mcparticle->GetRecoCherenkovPhotonCount(); + recoCherenkovHitCount = mcparticle->GetRecoCherenkovHitCount(); + recoTrackCherenkovAngle = 1000.0*mcparticle->GetRecoCherenkovAverageTheta(a1); + + tree->Fill(); + + //if(TMath::Abs(mcparticle->GetRecoPdgCode()) == 11) numElec++; + //if(TMath::Abs(mcparticle->GetRecoPdgCode()) == 211) numPion++; + //if(TMath::Abs(mcparticle->GetRecoPdgCode()) == 321) numKaon++; + //if(TMath::Abs(mcparticle->GetRecoPdgCode()) == 2212) numProt++; + + //hthtr1->Fill(1000*mcparticle->GetRecoCherenkovAverageTheta(a1)); + //hthtr2->Fill(1000*mcparticle->GetRecoCherenkovAverageTheta(a2)); + } //for mcparticle + } //while + } + + //factoryPhotons = (TH1D *)reco->hnpetr()->Clone("factoryPhotons"); + + ///hEta->Write(); + //hMom->Write(); + //factoryPhotons->Write(); + tree->Write(); + ofile->Close(); + +} From a35033c6b753354635fbc096b2f5d0296c1ba72f Mon Sep 17 00:00:00 2001 From: Youqi Song Date: Fri, 3 May 2024 20:49:17 -0400 Subject: [PATCH 02/10] Added tracking smearing and new Tree branches --- scripts/reco-epic-LUT.C | 152 +++++++++++++++------------------------- 1 file changed, 57 insertions(+), 95 deletions(-) diff --git a/scripts/reco-epic-LUT.C b/scripts/reco-epic-LUT.C index 4fd4917..e1f2a0d 100644 --- a/scripts/reco-epic-LUT.C +++ b/scripts/reco-epic-LUT.C @@ -1,17 +1,16 @@ - - void reco_epic_LUT(const char *dfname, const char *output, const char *cfname = 0) { auto *reco = new ReconstructionFactory(dfname, cfname, "pfRICH"); - //const char* output="test.hist.root"; - TFile *ofile = TFile::Open(output,"recreate"); + // const char* output="test.hist.root"; + TFile *ofile = TFile::Open(output, "recreate"); - TTree *tree = new TTree("T","Tree of Quantities"); + TTree *tree = new TTree("T", "Tree of Quantities"); Int_t partPDG; Double_t partMom; Double_t partEta; + Double_t partTheta; Double_t partPhi; Double_t partVertX; Double_t partVertY; @@ -19,48 +18,45 @@ void reco_epic_LUT(const char *dfname, const char *output, const char *cfname = Int_t recoPDG; Int_t recoPhotons; - Int_t recoCherenkovHitCount; Double_t recoTrackCherenkovAngle; - tree->Branch("partPDG",&partPDG,"partPDG/I"); - tree->Branch("partMom",&partMom,"partMom/D"); - tree->Branch("partEta",&partEta,"partEta/D"); - tree->Branch("partPhi",&partPhi,"partPhi/D"); - tree->Branch("partVertX",&partVertX,"partVertX/D"); - tree->Branch("partVertY",&partVertY,"partVertY/D"); - tree->Branch("partVertZ",&partVertZ,"partVertZ/D"); + tree->Branch("partPDG", &partPDG, "partPDG/I"); + tree->Branch("partMom", &partMom, "partMom/D"); + tree->Branch("partEta", &partEta, "partEta/D"); + tree->Branch("partTheta", &partTheta, "partTheta/D"); + tree->Branch("partPhi", &partPhi, "partPhi/D"); + tree->Branch("partVertX", &partVertX, "partVertX/D"); + tree->Branch("partVertY", &partVertY, "partVertY/D"); + tree->Branch("partVertZ", &partVertZ, "partVertZ/D"); - tree->Branch("recoPDG",&recoPDG,"recoPDG/I"); - tree->Branch("recoPhotons",&recoPhotons,"recoPhotons/I"); - tree->Branch("recoCherenkovHitCount",&recoCherenkovHitCount,"recoCherenkovHitCount/I"); - tree->Branch("recoTrackCherenkovAngle",&recoTrackCherenkovAngle,"recoTrackCherenkovAngle/D"); + tree->Branch("recoPDG", &recoPDG, "recoPDG/I"); + tree->Branch("recoPhotons", &recoPhotons, "recoPhotons/I"); + tree->Branch("recoTrackCherenkovAngle", &recoTrackCherenkovAngle, "recoTrackCherenkovAngle/D"); - // // Factory configuration part; // - //reco->IgnoreTimingInChiSquare(); - //reco->IgnorePoissonTermInChiSquare(); - //reco->SetSingleHitCCDFcut(0.005); - //reco->RemoveAmbiguousHits(); + // reco->IgnoreTimingInChiSquare(); + // reco->IgnorePoissonTermInChiSquare(); + // reco->SetSingleHitCCDFcut(0.005); + // reco->RemoveAmbiguousHits(); // This only affects the calibration stage; - //reco->SetDefaultSinglePhotonThetaResolution(0.0040); + // reco->SetDefaultSinglePhotonThetaResolution(0.0040); // Sensor active area pixelated will be pixellated NxN in digitization; - //reco->SetSensorActiveAreaPixellation(24); + // reco->SetSensorActiveAreaPixellation(24); // [rad] (should match SPE sigma) & [ns]; - //auto *a1 = reco->UseRadiator("Aerogel225", 0.0040); + // auto *a1 = reco->UseRadiator("Aerogel225", 0.0040); auto *a1 = reco->UseRadiator("BelleIIAerogel3"); - //auto *a1 = reco->UseRadiator("BelleIIAerogel3"); + // auto *a1 = reco->UseRadiator("BelleIIAerogel3"); - - //reco->SetSinglePhotonTimingResolution(0.030); - //reco->SetQuietMode(); - //reco->AddHypothesis(-11); + // reco->SetSinglePhotonTimingResolution(0.030); + // reco->SetQuietMode(); + // reco->AddHypothesis(-11); reco->AddHypothesis(11); reco->AddHypothesis("pi+"); reco->AddHypothesis(321); reco->AddHypothesis(2212); - //reco->IgnoreMcTruthPhotonDirectionSeed(); + // reco->IgnoreMcTruthPhotonDirectionSeed(); // Mark all pads hit by "calibration" (above the QE curve) photons as "useless"; reco->AddBlackoutRadiator("QuartzWindow"); @@ -68,80 +64,46 @@ void reco_epic_LUT(const char *dfname, const char *output, const char *cfname = reco->AddBlackoutRadiator("GasVolume"); // Carelessly remove (0x1 << n)x(0x1 << n) square area "around" these hits; reco->SetBlackoutBlowupValue(3); - - //auto hEta = new TH1D("hEta","",100,-5.,5.); - //auto hMom = new TH1D("hMom","",100,0.,10.); - - //auto factoryPhotons = new TH1D("factoryPhotons","",50,0.,50.); - - //auto hmatch = new TH1D("hmatch", "PID evaluation correctness", 2, 0, 2); - //auto hthtr1 = new TH1D("thtr1", "Cherenkov angle (track)", 200, 220, 320); - //auto hthtr1 = new TH1D("thtr1", "Cherenkov angle (track)", 40, 270, 290); - //auto hthtr1 = new TH1D("thtr1", "Cherenkov angle (track)", 580, 0., 290.); - // For a dual aerogel configuration; - //auto hthtr2 = new TH1D("thtr2", "Cherenkov angle (track)", 200, 220, 320); - + reco->ImportTrackingSmearing("./database/dtheta_seed_param.reformatted.dat", "./database/dphi_seed_param.reformatted.dat"); + reco->UseAutomaticCalibration(); // This call is mandatory; second argument: statistics (default: all events); reco->PerformCalibration(200); - //int numElec = 0; - //int numPion = 0; - //int numKaon = 0; - //int numProt = 0; + // int numElec = 0; + // int numPion = 0; + // int numKaon = 0; + // int numProt = 0; { CherenkovEvent *event; // Event loop; - while((event = reco->GetNextEvent())) { - for(auto mcparticle: event->ChargedParticles()) { - if (!mcparticle->IsPrimary()) continue; - - /* - if (mcparticle->GetPDG() == mcparticle->GetRecoPdgCode()) { - hmatch->Fill(0.5); - } else { - hmatch->Fill(1.5); - } //if - */ - - // GetRecoCherenkovPhotonCount - - //hEta->Fill(mcparticle->GetVertexMomentum().PseudoRapidity()); - //hMom->Fill(mcparticle->GetVertexMomentum().Mag()); - - partPDG = mcparticle->GetPDG(); - partMom = mcparticle->GetVertexMomentum().Mag(); - partEta = mcparticle->GetVertexMomentum().PseudoRapidity(); - partPhi = mcparticle->GetVertexMomentum().Phi(); - partVertX = mcparticle->GetVertexPosition().X(); - partVertY = mcparticle->GetVertexPosition().Y(); - partVertZ = mcparticle->GetVertexPosition().Z(); - - recoPDG = mcparticle->GetRecoPdgCode(); - recoPhotons = mcparticle->GetRecoCherenkovPhotonCount(); - recoCherenkovHitCount = mcparticle->GetRecoCherenkovHitCount(); - recoTrackCherenkovAngle = 1000.0*mcparticle->GetRecoCherenkovAverageTheta(a1); - - tree->Fill(); - - //if(TMath::Abs(mcparticle->GetRecoPdgCode()) == 11) numElec++; - //if(TMath::Abs(mcparticle->GetRecoPdgCode()) == 211) numPion++; - //if(TMath::Abs(mcparticle->GetRecoPdgCode()) == 321) numKaon++; - //if(TMath::Abs(mcparticle->GetRecoPdgCode()) == 2212) numProt++; - - //hthtr1->Fill(1000*mcparticle->GetRecoCherenkovAverageTheta(a1)); - //hthtr2->Fill(1000*mcparticle->GetRecoCherenkovAverageTheta(a2)); - } //for mcparticle - } //while + while ((event = reco->GetNextEvent())) + { + for (auto mcparticle : event->ChargedParticles()) + { + if (!mcparticle->IsPrimary()) + continue; + + partPDG = mcparticle->GetPDG(); + partMom = mcparticle->GetVertexMomentum().Mag(); + partEta = mcparticle->GetVertexMomentum().PseudoRapidity(); + partPhi = mcparticle->GetVertexMomentum().Phi(); + partTheta = mcparticle->GetVertexMomentum().Theta(); + partVertX = mcparticle->GetVertexPosition().X(); + partVertY = mcparticle->GetVertexPosition().Y(); + partVertZ = mcparticle->GetVertexPosition().Z(); + + recoPDG = mcparticle->GetRecoPdgCode(); + recoPhotons = mcparticle->GetRecoCherenkovPhotonCount(); + recoTrackCherenkovAngle = 1000.0 * mcparticle->GetRecoCherenkovAverageTheta(a1); + + tree->Fill(); + + } // for mcparticle + } // while } - //factoryPhotons = (TH1D *)reco->hnpetr()->Clone("factoryPhotons"); - - ///hEta->Write(); - //hMom->Write(); - //factoryPhotons->Write(); tree->Write(); ofile->Close(); - } From b80fbaa723f4816217e6d9bed928381b64cc201c Mon Sep 17 00:00:00 2001 From: Youqi Song Date: Fri, 3 May 2024 21:53:03 -0400 Subject: [PATCH 03/10] Added scripts to read trees and output LUT --- readTree.C | 244 ++++++++++++++++++++++++++++++++++++ readTreeQA.C | 347 +++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 591 insertions(+) create mode 100644 readTree.C create mode 100644 readTreeQA.C diff --git a/readTree.C b/readTree.C new file mode 100644 index 0000000..9085192 --- /dev/null +++ b/readTree.C @@ -0,0 +1,244 @@ +void readTree(const char *input, const char *output, const char *output_txt) +{ + + cout << "running readTree.C" << endl; + + TFile *fa = new TFile(input); + + TTree *ta = (TTree *)fa->Get("T"); + + TFile *ofile = TFile::Open(output, "recreate"); + + Int_t partPDG, recoPDG, recoPhotons, recoCherenkovHitCount; + Double_t partMom, partEta, partTheta, partPhi, partVertX, partVertY, partVertZ, recoTrackCherenkovAngle; + Double_t zero = pow(10., -50.); + + ta->SetBranchAddress("partPDG", &partPDG); + ta->SetBranchAddress("partMom", &partMom); + ta->SetBranchAddress("partEta", &partEta); + ta->SetBranchAddress("partTheta", &partTheta); + ta->SetBranchAddress("partPhi", &partPhi); + ta->SetBranchAddress("partVertX", &partVertX); + ta->SetBranchAddress("partVertY", &partVertY); + ta->SetBranchAddress("partVertZ", &partVertZ); + ta->SetBranchAddress("recoPDG", &recoPDG); + ta->SetBranchAddress("recoPhotons", &recoPhotons); + ta->SetBranchAddress("recoTrackCherenkovAngle", &recoTrackCherenkovAngle); + + TH1D *h_p = new TH1D("p", "", 15, 0.1, 15.1); + TH1D *h_t = new TH1D("theta", "", 20, 2.7, 3.1); + TH1D *h_phi = new TH1D("phi", "", 360, -3.1415926, 3.1415926); + + static Int_t fee[15][20][360], fepi[15][20][360], fek[15][20][360], fep[15][20][360], feb[15][20][360]; + static Int_t fpie[15][20][360], fpipi[15][20][360], fpik[15][20][360], fpip[15][20][360], fpib[15][20][360]; + static Int_t fke[15][20][360], fkpi[15][20][360], fkk[15][20][360], fkp[15][20][360], fkb[15][20][360]; + static Int_t fpe[15][20][360], fppi[15][20][360], fpk[15][20][360], fpp[15][20][360], fpb[15][20][360]; + + for (int pbin = 0; pbin < 15; pbin++) + { + for (int tbin = 0; tbin < 20; tbin++) + { + for (int phibin = 0; phibin < 360; phibin++) + { + fee[pbin][tbin][phibin] = fepi[pbin][tbin][phibin] = fek[pbin][tbin][phibin] = fep[pbin][tbin][phibin] = feb[pbin][tbin][phibin] = 0; + fpie[pbin][tbin][phibin] = fpipi[pbin][tbin][phibin] = fpik[pbin][tbin][phibin] = fpip[pbin][tbin][phibin] = fpib[pbin][tbin][phibin] = 0; + fke[pbin][tbin][phibin] = fkpi[pbin][tbin][phibin] = fkk[pbin][tbin][phibin] = fkp[pbin][tbin][phibin] = fkb[pbin][tbin][phibin] = 0; + fpe[pbin][tbin][phibin] = fppi[pbin][tbin][phibin] = fpk[pbin][tbin][phibin] = fpp[pbin][tbin][phibin] = fpb[pbin][tbin][phibin] = 0; + } + } + } + + Int_t nEntries = (Int_t)ta->GetEntries(); + for (int i = 0; i < nEntries; i++) + { + if (i % 100000 == 0) + cout << "Processing Event " << i << endl; + ta->GetEntry(i); + + Int_t pbin = h_p->FindBin(partMom) - 1; + Int_t tbin = h_t->FindBin(partTheta) - 1; + Int_t phibin = h_phi->FindBin(partPhi) - 1; + + // if (partMom < 1) cout << pbin << endl;// -> bin number starts from 1 + + if (partPDG == 11) // Look at Thrown Electrons + { + if (recoPhotons < 4) + { + feb[pbin][tbin][phibin]++; + } + else + { + if (recoPDG == 11) + fee[pbin][tbin][phibin]++; + if (recoPDG == 211) + fepi[pbin][tbin][phibin]++; + if (recoPDG == 321) + fek[pbin][tbin][phibin]++; + if (recoPDG == 2212) + fep[pbin][tbin][phibin]++; + } + } + + if (partPDG == 211) // Look at Thrown Pions + { + if (recoPhotons < 4) + { + fpib[pbin][tbin][phibin]++; + } + else + { + if (recoPDG == 11) + fpie[pbin][tbin][phibin]++; + if (recoPDG == 211) + fpipi[pbin][tbin][phibin]++; + if (recoPDG == 321) + fpik[pbin][tbin][phibin]++; + if (recoPDG == 2212) + fpip[pbin][tbin][phibin]++; + } + } + + if (partPDG == 321) // Look at Thrown Kaons + { + if (recoPhotons < 4) + { + fkb[pbin][tbin][phibin]++; + } + else + { + if (recoPDG == 11) + fke[pbin][tbin][phibin]++; + if (recoPDG == 211) + fkpi[pbin][tbin][phibin]++; + if (recoPDG == 321) + fkk[pbin][tbin][phibin]++; + if (recoPDG == 2212) + fkp[pbin][tbin][phibin]++; + } + } + + if (partPDG == 2212) // Look at Thrown Protons + { + if (recoPhotons < 4) + { + fpb[pbin][tbin][phibin]++; + } + else + { + if (recoPDG == 11) + fpe[pbin][tbin][phibin]++; + if (recoPDG == 211) + fppi[pbin][tbin][phibin]++; + if (recoPDG == 321) + fpk[pbin][tbin][phibin]++; + if (recoPDG == 2212) + fpp[pbin][tbin][phibin]++; + } + } + } + + TH2D *h[15][20][360]; + for (int pbin = 0; pbin < 15; pbin++) + { + for (int tbin = 0; tbin < 20; tbin++) + { + for (int phibin = 0; phibin < 360; phibin++) + { + h[pbin][tbin][phibin] = new TH2D(Form("p%it%iphi%i", pbin, tbin, phibin), Form("p%it%iphi%i;reco;truth", pbin, tbin, phibin), 5, 0., 5., 4, 0., 4.); + } + } + } + + FILE *filePointer; + filePointer = fopen(output_txt, "w"); + + for (int pbin = 0; pbin < 15; pbin++) + { + for (int tbin = 0; tbin < 20; tbin++) + { + for (int phibin = 0; phibin < 360; phibin++) + { + + Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1); + Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1); + Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); + + Double_t mc_e_rec_e = (1.0 * fee[pbin][tbin][phibin]) / (1.0 * (fee[pbin][tbin][phibin] + fepi[pbin][tbin][phibin] + fek[pbin][tbin][phibin] + fep[pbin][tbin][phibin] + feb[pbin][tbin][phibin] + zero)); + Double_t mc_e_rec_pi = (1.0 * fepi[pbin][tbin][phibin]) / (1.0 * (fee[pbin][tbin][phibin] + fepi[pbin][tbin][phibin] + fek[pbin][tbin][phibin] + fep[pbin][tbin][phibin] + feb[pbin][tbin][phibin] + zero)); + Double_t mc_e_rec_k = (1.0 * fek[pbin][tbin][phibin]) / (1.0 * (fee[pbin][tbin][phibin] + fepi[pbin][tbin][phibin] + fek[pbin][tbin][phibin] + fep[pbin][tbin][phibin] + feb[pbin][tbin][phibin] + zero)); + Double_t mc_e_rec_p = (1.0 * fep[pbin][tbin][phibin]) / (1.0 * (fee[pbin][tbin][phibin] + fepi[pbin][tbin][phibin] + fek[pbin][tbin][phibin] + fep[pbin][tbin][phibin] + feb[pbin][tbin][phibin] + zero)); + Double_t mc_e_rec_b = (1.0 * feb[pbin][tbin][phibin]) / (1.0 * (fee[pbin][tbin][phibin] + fepi[pbin][tbin][phibin] + fek[pbin][tbin][phibin] + fep[pbin][tbin][phibin] + feb[pbin][tbin][phibin] + zero)); + + Double_t mc_pi_rec_e = (1.0 * fpie[pbin][tbin][phibin]) / (1.0 * (fpie[pbin][tbin][phibin] + fpipi[pbin][tbin][phibin] + fpik[pbin][tbin][phibin] + fpip[pbin][tbin][phibin] + fpib[pbin][tbin][phibin] + zero)); + Double_t mc_pi_rec_pi = (1.0 * fpipi[pbin][tbin][phibin]) / (1.0 * (fpie[pbin][tbin][phibin] + fpipi[pbin][tbin][phibin] + fpik[pbin][tbin][phibin] + fpip[pbin][tbin][phibin] + fpib[pbin][tbin][phibin] + zero)); + Double_t mc_pi_rec_k = (1.0 * fpik[pbin][tbin][phibin]) / (1.0 * (fpie[pbin][tbin][phibin] + fpipi[pbin][tbin][phibin] + fpik[pbin][tbin][phibin] + fpip[pbin][tbin][phibin] + fpib[pbin][tbin][phibin] + zero)); + Double_t mc_pi_rec_p = (1.0 * fpip[pbin][tbin][phibin]) / (1.0 * (fpie[pbin][tbin][phibin] + fpipi[pbin][tbin][phibin] + fpik[pbin][tbin][phibin] + fpip[pbin][tbin][phibin] + fpib[pbin][tbin][phibin] + zero)); + Double_t mc_pi_rec_b = (1.0 * fpib[pbin][tbin][phibin]) / (1.0 * (fpie[pbin][tbin][phibin] + fpipi[pbin][tbin][phibin] + fpik[pbin][tbin][phibin] + fpip[pbin][tbin][phibin] + fpib[pbin][tbin][phibin] + zero)); + + Double_t mc_k_rec_e = (1.0 * fke[pbin][tbin][phibin]) / (1.0 * (fke[pbin][tbin][phibin] + fkpi[pbin][tbin][phibin] + fkk[pbin][tbin][phibin] + fkp[pbin][tbin][phibin] + fkb[pbin][tbin][phibin] + zero)); + Double_t mc_k_rec_pi = (1.0 * fkpi[pbin][tbin][phibin]) / (1.0 * (fke[pbin][tbin][phibin] + fkpi[pbin][tbin][phibin] + fkk[pbin][tbin][phibin] + fkp[pbin][tbin][phibin] + fkb[pbin][tbin][phibin] + zero)); + Double_t mc_k_rec_k = (1.0 * fkk[pbin][tbin][phibin]) / (1.0 * (fke[pbin][tbin][phibin] + fkpi[pbin][tbin][phibin] + fkk[pbin][tbin][phibin] + fkp[pbin][tbin][phibin] + fkb[pbin][tbin][phibin] + zero)); + Double_t mc_k_rec_p = (1.0 * fkp[pbin][tbin][phibin]) / (1.0 * (fke[pbin][tbin][phibin] + fkpi[pbin][tbin][phibin] + fkk[pbin][tbin][phibin] + fkp[pbin][tbin][phibin] + fkb[pbin][tbin][phibin] + zero)); + Double_t mc_k_rec_b = (1.0 * fkb[pbin][tbin][phibin]) / (1.0 * (fke[pbin][tbin][phibin] + fkpi[pbin][tbin][phibin] + fkk[pbin][tbin][phibin] + fkp[pbin][tbin][phibin] + fkb[pbin][tbin][phibin] + zero)); + + Double_t mc_p_rec_e = (1.0 * fpe[pbin][tbin][phibin]) / (1.0 * (fpe[pbin][tbin][phibin] + fppi[pbin][tbin][phibin] + fpk[pbin][tbin][phibin] + fpp[pbin][tbin][phibin] + fpb[pbin][tbin][phibin] + zero)); + Double_t mc_p_rec_pi = (1.0 * fppi[pbin][tbin][phibin]) / (1.0 * (fpe[pbin][tbin][phibin] + fppi[pbin][tbin][phibin] + fpk[pbin][tbin][phibin] + fpp[pbin][tbin][phibin] + fpb[pbin][tbin][phibin] + zero)); + Double_t mc_p_rec_k = (1.0 * fpk[pbin][tbin][phibin]) / (1.0 * (fpe[pbin][tbin][phibin] + fppi[pbin][tbin][phibin] + fpk[pbin][tbin][phibin] + fpp[pbin][tbin][phibin] + fpb[pbin][tbin][phibin] + zero)); + Double_t mc_p_rec_p = (1.0 * fpp[pbin][tbin][phibin]) / (1.0 * (fpe[pbin][tbin][phibin] + fppi[pbin][tbin][phibin] + fpk[pbin][tbin][phibin] + fpp[pbin][tbin][phibin] + fpb[pbin][tbin][phibin] + zero)); + Double_t mc_p_rec_b = (1.0 * fpb[pbin][tbin][phibin]) / (1.0 * (fpe[pbin][tbin][phibin] + fppi[pbin][tbin][phibin] + fpk[pbin][tbin][phibin] + fpp[pbin][tbin][phibin] + fpb[pbin][tbin][phibin] + zero)); + + h[pbin][tbin][phibin]->SetBinContent(1, 4, mc_e_rec_e); + h[pbin][tbin][phibin]->SetBinContent(2, 4, mc_e_rec_pi); + h[pbin][tbin][phibin]->SetBinContent(3, 4, mc_e_rec_k); + h[pbin][tbin][phibin]->SetBinContent(4, 4, mc_e_rec_p); + h[pbin][tbin][phibin]->SetBinContent(5, 4, mc_e_rec_b); + + h[pbin][tbin][phibin]->SetBinContent(1, 3, mc_pi_rec_e); + h[pbin][tbin][phibin]->SetBinContent(2, 3, mc_pi_rec_pi); + h[pbin][tbin][phibin]->SetBinContent(3, 3, mc_pi_rec_k); + h[pbin][tbin][phibin]->SetBinContent(4, 3, mc_pi_rec_p); + h[pbin][tbin][phibin]->SetBinContent(5, 3, mc_pi_rec_b); + + h[pbin][tbin][phibin]->SetBinContent(1, 2, mc_k_rec_e); + h[pbin][tbin][phibin]->SetBinContent(2, 2, mc_k_rec_pi); + h[pbin][tbin][phibin]->SetBinContent(3, 2, mc_k_rec_k); + h[pbin][tbin][phibin]->SetBinContent(4, 2, mc_k_rec_p); + h[pbin][tbin][phibin]->SetBinContent(5, 2, mc_k_rec_b); + + h[pbin][tbin][phibin]->SetBinContent(1, 1, mc_p_rec_e); + h[pbin][tbin][phibin]->SetBinContent(2, 1, mc_p_rec_pi); + h[pbin][tbin][phibin]->SetBinContent(3, 1, mc_p_rec_k); + h[pbin][tbin][phibin]->SetBinContent(4, 1, mc_p_rec_p); + h[pbin][tbin][phibin]->SetBinContent(5, 1, mc_p_rec_b); + + // h[pbin][tbin][phibin]->Write(); + + fprintf(filePointer, "11 %f %f %f %f %f %f %f %f\n", p, t, phi, mc_e_rec_e, mc_e_rec_pi, mc_e_rec_k, mc_e_rec_p, mc_e_rec_b); + fprintf(filePointer, "211 %f %f %f %f %f %f %f %f\n", p, t, phi, mc_pi_rec_e, mc_pi_rec_pi, mc_pi_rec_k, mc_pi_rec_p, mc_pi_rec_b); + fprintf(filePointer, "321 %f %f %f %f %f %f %f %f\n", p, t, phi, mc_k_rec_e, mc_k_rec_pi, mc_k_rec_k, mc_k_rec_p, mc_k_rec_b); + fprintf(filePointer, "2212 %f %f %f %f %f %f %f %f\n", p, t, phi, mc_p_rec_e, mc_p_rec_pi, mc_p_rec_k, mc_p_rec_p, mc_p_rec_b); + + if (pbin == 0 && tbin == 0 && phibin == 0) + { + cout << "truth PID, p, theta, phi, prob(e), prob(pi), prob(k), prob(p), prob(fail)" << endl; + cout << "11 " << p << " " << t << " " << phi << " " << mc_e_rec_e << " " << mc_e_rec_pi << " " << mc_e_rec_k << " " << mc_e_rec_p << " " << mc_e_rec_b << endl; + cout << "211 " << p << " " << t << " " << phi << " " << mc_pi_rec_e << " " << mc_pi_rec_pi << " " << mc_pi_rec_k << " " << mc_pi_rec_p << " " << mc_pi_rec_b << endl; + cout << "321 " << p << " " << t << " " << phi << " " << mc_k_rec_e << " " << mc_k_rec_pi << " " << mc_k_rec_k << " " << mc_k_rec_p << " " << mc_k_rec_b << endl; + cout << "2212 " << p << " " << t << " " << phi << " " << mc_p_rec_e << " " << mc_p_rec_pi << " " << mc_p_rec_k << " " << mc_p_rec_p << " " << mc_p_rec_b << endl; + } + if (pbin == 14 && tbin == 9 && phibin == 11) + { + cout << "11 " << p << " " << t << " " << phi << " " << mc_e_rec_e << " " << mc_e_rec_pi << " " << mc_e_rec_k << " " << mc_e_rec_p << " " << mc_e_rec_b << endl; + cout << "211 " << p << " " << t << " " << phi << " " << mc_pi_rec_e << " " << mc_pi_rec_pi << " " << mc_pi_rec_k << " " << mc_pi_rec_p << " " << mc_pi_rec_b << endl; + cout << "321 " << p << " " << t << " " << phi << " " << mc_k_rec_e << " " << mc_k_rec_pi << " " << mc_k_rec_k << " " << mc_k_rec_p << " " << mc_k_rec_b << endl; + cout << "2212 " << p << " " << t << " " << phi << " " << mc_p_rec_e << " " << mc_p_rec_pi << " " << mc_p_rec_k << " " << mc_p_rec_p << " " << mc_p_rec_b << endl; + } + } + } + } + + fclose(filePointer); + ofile->Write(); + ofile->Close(); +} diff --git a/readTreeQA.C b/readTreeQA.C new file mode 100644 index 0000000..b561200 --- /dev/null +++ b/readTreeQA.C @@ -0,0 +1,347 @@ +void readTreeQA(const char *input, const char *output, const char *output_txt) +{ + + TFile *fa = new TFile(input); + + TTree *ta = (TTree *)fa->Get("T"); + + TFile *ofile = TFile::Open(output, "recreate"); + + Int_t partPDG, recoPDG, recoPhotons, recoCherenkovHitCount; + Double_t partMom, partEta, partTheta, partPhi, partVertX, partVertY, partVertZ, recoTrackCherenkovAngle; + Double_t zero = pow(10., -50.); + + ta->SetBranchAddress("partPDG", &partPDG); + ta->SetBranchAddress("partMom", &partMom); + ta->SetBranchAddress("partEta", &partEta); + ta->SetBranchAddress("partTheta", &partTheta); + ta->SetBranchAddress("partPhi", &partPhi); + ta->SetBranchAddress("partVertX", &partVertX); + ta->SetBranchAddress("partVertY", &partVertY); + ta->SetBranchAddress("partVertZ", &partVertZ); + ta->SetBranchAddress("recoPDG", &recoPDG); + ta->SetBranchAddress("recoPhotons", &recoPhotons); + ta->SetBranchAddress("recoTrackCherenkovAngle", &recoTrackCherenkovAngle); + + TH1D *h_p = new TH1D("p", "", 15, 0.1, 15.1); + TH1D *h_t = new TH1D("theta", "", 16, 2.7, 3.1); + TH1D *h_phi = new TH1D("phi", "", 360, -3.1415926, 3.1415926); + + TProfile2D *recoPhiVsEtaElecPhoton[4]; + TH2D *partPhiVsEtaElec[4]; + TH2D *partPhiVsEtaElecMatch[4]; + TH2D *partPhiVsEtaElecUnMatch[4]; + TH2D *partPhiVsEtaElecMatchPhoton[4]; + + TProfile2D *recoPhiVsEtaPionPhoton[4]; + TH2D *partPhiVsEtaPion[4]; + TH2D *partPhiVsEtaPionMatch[4]; + TH2D *partPhiVsEtaPionUnMatch[4]; + + TProfile2D *recoPhiVsEtaKaonPhoton[4]; + TH2D *partPhiVsEtaKaon[4]; + TH2D *partPhiVsEtaKaonMatch[4]; + TH2D *partPhiVsEtaKaonUnMatch[4]; + + TProfile2D *recoPhiVsEtaProtonPhoton[4]; + TH2D *partPhiVsEtaProton[4]; + TH2D *partPhiVsEtaProtonMatch[4]; + TH2D *partPhiVsEtaProtonUnMatch[4]; + + for (int i = 0; i < 4; i++) + { + recoPhiVsEtaElecPhoton[i] = new TProfile2D(Form("recoPhiVsEtaElecPhoton_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + partPhiVsEtaElec[i] = new TH2D(Form("partPhiVsEtaElecMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + partPhiVsEtaElecMatch[i] = new TH2D(Form("partPhiVsEtaElecMatchMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + partPhiVsEtaElecUnMatch[i] = new TH2D(Form("partPhiVsEtaElecUnMatchMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + partPhiVsEtaElecMatchPhoton[i] = new TH2D(Form("partPhiVsEtaElecMatchPhotonMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + + recoPhiVsEtaPionPhoton[i] = new TProfile2D(Form("recoPhiVsEtaPionPhoton_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + partPhiVsEtaPion[i] = new TH2D(Form("partPhiVsEtaPionMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + partPhiVsEtaPionMatch[i] = new TH2D(Form("partPhiVsEtaPionMatchMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + partPhiVsEtaPionUnMatch[i] = new TH2D(Form("partPhiVsEtaPionUnMatchMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + + recoPhiVsEtaKaonPhoton[i] = new TProfile2D(Form("recoPhiVsEtaKaonPhoton_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + partPhiVsEtaKaon[i] = new TH2D(Form("partPhiVsEtaKaonMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + partPhiVsEtaKaonMatch[i] = new TH2D(Form("partPhiVsEtaKaonMatchMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + partPhiVsEtaKaonUnMatch[i] = new TH2D(Form("partPhiVsEtaKaonUnMatchMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + + recoPhiVsEtaProtonPhoton[i] = new TProfile2D(Form("recoPhiVsEtaProtonPhoton_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + partPhiVsEtaProton[i] = new TH2D(Form("partPhiVsEtaProtonMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + partPhiVsEtaProtonMatch[i] = new TH2D(Form("partPhiVsEtaProtonMatchMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + partPhiVsEtaProtonUnMatch[i] = new TH2D(Form("partPhiVsEtaProtonUnMatchMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + } + + TH2D *numPhotonsVsPhiEta25Mom10Pion = new TH2D("numPhotonsVsPhiEta25Mom10Pion", "", 60, -TMath::Pi(), TMath::Pi(), 25, 0., 25.); + TH2D *numPhotonsVsPhiEta25Mom4Pion = new TH2D("numPhotonsVsPhiEta25Mom4Pion", "", 60, -TMath::Pi(), TMath::Pi(), 25, 0., 25.); + + Int_t fee[15][16][360], fepi[15][16][360], fek[15][16][360], fep[15][16][360], feb[15][16][360]; + Int_t fpie[15][16][360], fpipi[15][16][360], fpik[15][16][360], fpip[15][16][360], fpib[15][16][360]; + Int_t fke[15][16][360], fkpi[15][16][360], fkk[15][16][360], fkp[15][16][360], fkb[15][16][360]; + Int_t fpe[15][16][360], fppi[15][16][360], fpk[15][16][360], fpp[15][16][360], fpb[15][16][360]; + + for (int pbin = 0; pbin < 15; pbin++) + { + for (int tbin = 0; tbin < 16; tbin++) + { + for (int phibin = 0; phibin < 360; phibin++) + { + fee[pbin][tbin][phibin] = fepi[pbin][tbin][phibin] = fek[pbin][tbin][phibin] = fep[pbin][tbin][phibin] = feb[pbin][tbin][phibin] = 0; + fpie[pbin][tbin][phibin] = fpipi[pbin][tbin][phibin] = fpik[pbin][tbin][phibin] = fpip[pbin][tbin][phibin] = fpib[pbin][tbin][phibin] = 0; + fke[pbin][tbin][phibin] = fkpi[pbin][tbin][phibin] = fkk[pbin][tbin][phibin] = fkp[pbin][tbin][phibin] = fkb[pbin][tbin][phibin] = 0; + fpe[pbin][tbin][phibin] = fppi[pbin][tbin][phibin] = fpk[pbin][tbin][phibin] = fpp[pbin][tbin][phibin] = fpb[pbin][tbin][phibin] = 0; + } + } + } + + Int_t nEntries = (Int_t)ta->GetEntries(); + for (int i = 0; i < nEntries; i++) + { + if (i % 100000 == 0) + cout << "Processing Event " << i << endl; + ta->GetEntry(i); + + Int_t pbin = h_p->FindBin(partMom) - 1; + Int_t tbin = h_t->FindBin(partTheta) - 1; + Int_t phibin = h_phi->FindBin(partPhi) - 1; + + Int_t momSelect = -1; + if (partMom > 1.0 && partMom < 3.0) + momSelect = 0; + if (partMom > 3.0 && partMom < 5.0) + momSelect = 1; + if (partMom > 5.0 && partMom < 7.0) + momSelect = 2; + if (partMom > 7.0 && partMom < 9.0) + momSelect = 3; + + // if (partMom < 1) cout << pbin << endl;// -> bin number starts from 1 + + if (partPDG == 11) // Look at Thrown Electrons + { + if (recoPhotons < 4) + { + feb[pbin][tbin][phibin]++; + } + else + { + if (recoPDG == 11) + fee[pbin][tbin][phibin]++; + if (recoPDG == 211) + fepi[pbin][tbin][phibin]++; + if (recoPDG == 321) + fek[pbin][tbin][phibin]++; + if (recoPDG == 2212) + fep[pbin][tbin][phibin]++; + } + if (momSelect != -1) + { + recoPhiVsEtaElecPhoton[momSelect]->Fill(partEta, partPhi, recoPhotons); + partPhiVsEtaElec[momSelect]->Fill(partEta, partPhi); + if (recoPDG == 11) + partPhiVsEtaElecMatch[momSelect]->Fill(partEta, partPhi); + if (recoPDG != 11) + partPhiVsEtaElecUnMatch[momSelect]->Fill(partEta, partPhi); + if (recoPDG == 11 && recoPhotons > 0) + partPhiVsEtaElecMatchPhoton[momSelect]->Fill(partEta, partPhi); + } + } + + if (partPDG == 211) // Look at Thrown Pions + { + if (recoPhotons < 4) + { + fpib[pbin][tbin][phibin]++; + } + else + { + if (recoPDG == 11) + fpie[pbin][tbin][phibin]++; + if (recoPDG == 211) + fpipi[pbin][tbin][phibin]++; + if (recoPDG == 321) + fpik[pbin][tbin][phibin]++; + if (recoPDG == 2212) + fpip[pbin][tbin][phibin]++; + } + if (momSelect != -1) + { + recoPhiVsEtaPionPhoton[momSelect]->Fill(partEta, partPhi, recoPhotons); + partPhiVsEtaPion[momSelect]->Fill(partEta, partPhi); + if (recoPDG == 211) + partPhiVsEtaPionMatch[momSelect]->Fill(partEta, partPhi); + if (recoPDG != 211) + partPhiVsEtaPionUnMatch[momSelect]->Fill(partEta, partPhi); + } + if (partEta < -2.49 && partEta > -2.51) + { + if (partMom > 9.5 && partMom < 10.5) + { + numPhotonsVsPhiEta25Mom10Pion->Fill(partPhi, recoPhotons); + } + if (partMom > 3.5 && partMom < 4.5) + { + numPhotonsVsPhiEta25Mom4Pion->Fill(partPhi, recoPhotons); + } + } + } + + if (partPDG == 321) // Look at Thrown Kaons + { + if (recoPhotons < 4) + { + fkb[pbin][tbin][phibin]++; + } + else + { + if (recoPDG == 11) + fke[pbin][tbin][phibin]++; + if (recoPDG == 211) + fkpi[pbin][tbin][phibin]++; + if (recoPDG == 321) + fkk[pbin][tbin][phibin]++; + if (recoPDG == 2212) + fkp[pbin][tbin][phibin]++; + } + if (momSelect != -1) + { + recoPhiVsEtaKaonPhoton[momSelect]->Fill(partEta, partPhi, recoPhotons); + partPhiVsEtaKaon[momSelect]->Fill(partEta, partPhi); + if (recoPDG == 321) + partPhiVsEtaKaonMatch[momSelect]->Fill(partEta, partPhi); + if (recoPDG != 321) + partPhiVsEtaKaonUnMatch[momSelect]->Fill(partEta, partPhi); + } + } + + if (partPDG == 2212) // Look at Thrown Protons + { + if (recoPhotons < 4) + { + fpb[pbin][tbin][phibin]++; + } + else + { + if (recoPDG == 11) + fpe[pbin][tbin][phibin]++; + if (recoPDG == 211) + fppi[pbin][tbin][phibin]++; + if (recoPDG == 321) + fpk[pbin][tbin][phibin]++; + if (recoPDG == 2212) + fpp[pbin][tbin][phibin]++; + } + if (momSelect != -1) + { + recoPhiVsEtaProtonPhoton[momSelect]->Fill(partEta, partPhi, recoPhotons); + partPhiVsEtaProton[momSelect]->Fill(partEta, partPhi); + if (recoPDG == 2212) + partPhiVsEtaProtonMatch[momSelect]->Fill(partEta, partPhi); + if (recoPDG != 2212) + partPhiVsEtaProtonUnMatch[momSelect]->Fill(partEta, partPhi); + } + } + } + + TH2D *h[15][16][360]; + for (int pbin = 0; pbin < 15; pbin++) + { + for (int tbin = 0; tbin < 16; tbin++) + { + for (int phibin = 0; phibin < 360; phibin++) + { + h[pbin][tbin][phibin] = new TH2D(Form("p%it%iphi%i", pbin, tbin, phibin), Form("p%it%iphi%i;reco;truth", pbin, tbin, phibin), 5, 0., 5., 4, 0., 4.); + } + } + } + + FILE *filePointer; + filePointer = fopen(output_txt, "w"); + + for (int pbin = 0; pbin < 15; pbin++) + { + for (int tbin = 0; tbin < 16; tbin++) + { + for (int phibin = 0; phibin < 360; phibin++) + { + + Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1); + Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1); + Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); + + Double_t mc_e_rec_e = (1.0 * fee[pbin][tbin][phibin]) / (1.0 * (fee[pbin][tbin][phibin] + fepi[pbin][tbin][phibin] + fek[pbin][tbin][phibin] + fep[pbin][tbin][phibin] + feb[pbin][tbin][phibin] + zero)); + Double_t mc_e_rec_pi = (1.0 * fepi[pbin][tbin][phibin]) / (1.0 * (fee[pbin][tbin][phibin] + fepi[pbin][tbin][phibin] + fek[pbin][tbin][phibin] + fep[pbin][tbin][phibin] + feb[pbin][tbin][phibin] + zero)); + Double_t mc_e_rec_k = (1.0 * fek[pbin][tbin][phibin]) / (1.0 * (fee[pbin][tbin][phibin] + fepi[pbin][tbin][phibin] + fek[pbin][tbin][phibin] + fep[pbin][tbin][phibin] + feb[pbin][tbin][phibin] + zero)); + Double_t mc_e_rec_p = (1.0 * fep[pbin][tbin][phibin]) / (1.0 * (fee[pbin][tbin][phibin] + fepi[pbin][tbin][phibin] + fek[pbin][tbin][phibin] + fep[pbin][tbin][phibin] + feb[pbin][tbin][phibin] + zero)); + Double_t mc_e_rec_b = (1.0 * feb[pbin][tbin][phibin]) / (1.0 * (fee[pbin][tbin][phibin] + fepi[pbin][tbin][phibin] + fek[pbin][tbin][phibin] + fep[pbin][tbin][phibin] + feb[pbin][tbin][phibin] + zero)); + + Double_t mc_pi_rec_e = (1.0 * fpie[pbin][tbin][phibin]) / (1.0 * (fpie[pbin][tbin][phibin] + fpipi[pbin][tbin][phibin] + fpik[pbin][tbin][phibin] + fpip[pbin][tbin][phibin] + fpib[pbin][tbin][phibin] + zero)); + Double_t mc_pi_rec_pi = (1.0 * fpipi[pbin][tbin][phibin]) / (1.0 * (fpie[pbin][tbin][phibin] + fpipi[pbin][tbin][phibin] + fpik[pbin][tbin][phibin] + fpip[pbin][tbin][phibin] + fpib[pbin][tbin][phibin] + zero)); + Double_t mc_pi_rec_k = (1.0 * fpik[pbin][tbin][phibin]) / (1.0 * (fpie[pbin][tbin][phibin] + fpipi[pbin][tbin][phibin] + fpik[pbin][tbin][phibin] + fpip[pbin][tbin][phibin] + fpib[pbin][tbin][phibin] + zero)); + Double_t mc_pi_rec_p = (1.0 * fpip[pbin][tbin][phibin]) / (1.0 * (fpie[pbin][tbin][phibin] + fpipi[pbin][tbin][phibin] + fpik[pbin][tbin][phibin] + fpip[pbin][tbin][phibin] + fpib[pbin][tbin][phibin] + zero)); + Double_t mc_pi_rec_b = (1.0 * fpib[pbin][tbin][phibin]) / (1.0 * (fpie[pbin][tbin][phibin] + fpipi[pbin][tbin][phibin] + fpik[pbin][tbin][phibin] + fpip[pbin][tbin][phibin] + fpib[pbin][tbin][phibin] + zero)); + + Double_t mc_k_rec_e = (1.0 * fke[pbin][tbin][phibin]) / (1.0 * (fke[pbin][tbin][phibin] + fkpi[pbin][tbin][phibin] + fkk[pbin][tbin][phibin] + fkp[pbin][tbin][phibin] + fkb[pbin][tbin][phibin] + zero)); + Double_t mc_k_rec_pi = (1.0 * fkpi[pbin][tbin][phibin]) / (1.0 * (fke[pbin][tbin][phibin] + fkpi[pbin][tbin][phibin] + fkk[pbin][tbin][phibin] + fkp[pbin][tbin][phibin] + fkb[pbin][tbin][phibin] + zero)); + Double_t mc_k_rec_k = (1.0 * fkk[pbin][tbin][phibin]) / (1.0 * (fke[pbin][tbin][phibin] + fkpi[pbin][tbin][phibin] + fkk[pbin][tbin][phibin] + fkp[pbin][tbin][phibin] + fkb[pbin][tbin][phibin] + zero)); + Double_t mc_k_rec_p = (1.0 * fkp[pbin][tbin][phibin]) / (1.0 * (fke[pbin][tbin][phibin] + fkpi[pbin][tbin][phibin] + fkk[pbin][tbin][phibin] + fkp[pbin][tbin][phibin] + fkb[pbin][tbin][phibin] + zero)); + Double_t mc_k_rec_b = (1.0 * fkb[pbin][tbin][phibin]) / (1.0 * (fke[pbin][tbin][phibin] + fkpi[pbin][tbin][phibin] + fkk[pbin][tbin][phibin] + fkp[pbin][tbin][phibin] + fkb[pbin][tbin][phibin] + zero)); + + Double_t mc_p_rec_e = (1.0 * fpe[pbin][tbin][phibin]) / (1.0 * (fpe[pbin][tbin][phibin] + fppi[pbin][tbin][phibin] + fpk[pbin][tbin][phibin] + fpp[pbin][tbin][phibin] + fpb[pbin][tbin][phibin] + zero)); + Double_t mc_p_rec_pi = (1.0 * fppi[pbin][tbin][phibin]) / (1.0 * (fpe[pbin][tbin][phibin] + fppi[pbin][tbin][phibin] + fpk[pbin][tbin][phibin] + fpp[pbin][tbin][phibin] + fpb[pbin][tbin][phibin] + zero)); + Double_t mc_p_rec_k = (1.0 * fpk[pbin][tbin][phibin]) / (1.0 * (fpe[pbin][tbin][phibin] + fppi[pbin][tbin][phibin] + fpk[pbin][tbin][phibin] + fpp[pbin][tbin][phibin] + fpb[pbin][tbin][phibin] + zero)); + Double_t mc_p_rec_p = (1.0 * fpp[pbin][tbin][phibin]) / (1.0 * (fpe[pbin][tbin][phibin] + fppi[pbin][tbin][phibin] + fpk[pbin][tbin][phibin] + fpp[pbin][tbin][phibin] + fpb[pbin][tbin][phibin] + zero)); + Double_t mc_p_rec_b = (1.0 * fpb[pbin][tbin][phibin]) / (1.0 * (fpe[pbin][tbin][phibin] + fppi[pbin][tbin][phibin] + fpk[pbin][tbin][phibin] + fpp[pbin][tbin][phibin] + fpb[pbin][tbin][phibin] + zero)); + + h[pbin][tbin][phibin]->SetBinContent(1, 4, mc_e_rec_e); + h[pbin][tbin][phibin]->SetBinContent(2, 4, mc_e_rec_pi); + h[pbin][tbin][phibin]->SetBinContent(3, 4, mc_e_rec_k); + h[pbin][tbin][phibin]->SetBinContent(4, 4, mc_e_rec_p); + h[pbin][tbin][phibin]->SetBinContent(5, 4, mc_e_rec_b); + + h[pbin][tbin][phibin]->SetBinContent(1, 3, mc_pi_rec_e); + h[pbin][tbin][phibin]->SetBinContent(2, 3, mc_pi_rec_pi); + h[pbin][tbin][phibin]->SetBinContent(3, 3, mc_pi_rec_k); + h[pbin][tbin][phibin]->SetBinContent(4, 3, mc_pi_rec_p); + h[pbin][tbin][phibin]->SetBinContent(5, 3, mc_pi_rec_b); + + h[pbin][tbin][phibin]->SetBinContent(1, 2, mc_k_rec_e); + h[pbin][tbin][phibin]->SetBinContent(2, 2, mc_k_rec_pi); + h[pbin][tbin][phibin]->SetBinContent(3, 2, mc_k_rec_k); + h[pbin][tbin][phibin]->SetBinContent(4, 2, mc_k_rec_p); + h[pbin][tbin][phibin]->SetBinContent(5, 2, mc_k_rec_b); + + h[pbin][tbin][phibin]->SetBinContent(1, 1, mc_p_rec_e); + h[pbin][tbin][phibin]->SetBinContent(2, 1, mc_p_rec_pi); + h[pbin][tbin][phibin]->SetBinContent(3, 1, mc_p_rec_k); + h[pbin][tbin][phibin]->SetBinContent(4, 1, mc_p_rec_p); + h[pbin][tbin][phibin]->SetBinContent(5, 1, mc_p_rec_b); + + fprintf(filePointer, "11 %f %f %f %f %f %f %f %f\n", p, t, phi, mc_e_rec_e, mc_e_rec_pi, mc_e_rec_k, mc_e_rec_p, mc_e_rec_b); + fprintf(filePointer, "211 %f %f %f %f %f %f %f %f\n", p, t, phi, mc_pi_rec_e, mc_pi_rec_pi, mc_pi_rec_k, mc_pi_rec_p, mc_pi_rec_b); + fprintf(filePointer, "321 %f %f %f %f %f %f %f %f\n", p, t, phi, mc_k_rec_e, mc_k_rec_pi, mc_k_rec_k, mc_k_rec_p, mc_k_rec_b); + fprintf(filePointer, "2212 %f %f %f %f %f %f %f %f\n", p, t, phi, mc_p_rec_e, mc_p_rec_pi, mc_p_rec_k, mc_p_rec_p, mc_p_rec_b); + + if (pbin == 0 && tbin == 0 && phibin == 0) + { + cout << "truth PID, p, theta, phi, prob(e), prob(pi), prob(k), prob(p), prob(fail)" << endl; + cout << "11 " << p << " " << t << " " << phi << " " << mc_e_rec_e << " " << mc_e_rec_pi << " " << mc_e_rec_k << " " << mc_e_rec_p << " " << mc_e_rec_b << endl; + cout << "211 " << p << " " << t << " " << phi << " " << mc_pi_rec_e << " " << mc_pi_rec_pi << " " << mc_pi_rec_k << " " << mc_pi_rec_p << " " << mc_pi_rec_b << endl; + cout << "321 " << p << " " << t << " " << phi << " " << mc_k_rec_e << " " << mc_k_rec_pi << " " << mc_k_rec_k << " " << mc_k_rec_p << " " << mc_k_rec_b << endl; + cout << "2212 " << p << " " << t << " " << phi << " " << mc_p_rec_e << " " << mc_p_rec_pi << " " << mc_p_rec_k << " " << mc_p_rec_p << " " << mc_p_rec_b << endl; + } + if (pbin == 14 && tbin == 9 && phibin == 11) + { + cout << "11 " << p << " " << t << " " << phi << " " << mc_e_rec_e << " " << mc_e_rec_pi << " " << mc_e_rec_k << " " << mc_e_rec_p << " " << mc_e_rec_b << endl; + cout << "211 " << p << " " << t << " " << phi << " " << mc_pi_rec_e << " " << mc_pi_rec_pi << " " << mc_pi_rec_k << " " << mc_pi_rec_p << " " << mc_pi_rec_b << endl; + cout << "321 " << p << " " << t << " " << phi << " " << mc_k_rec_e << " " << mc_k_rec_pi << " " << mc_k_rec_k << " " << mc_k_rec_p << " " << mc_k_rec_b << endl; + cout << "2212 " << p << " " << t << " " << phi << " " << mc_p_rec_e << " " << mc_p_rec_pi << " " << mc_p_rec_k << " " << mc_p_rec_p << " " << mc_p_rec_b << endl; + } + } + } + } + + fclose(filePointer); + ofile->Write(); + ofile->Close(); +} From ee8305f6f154ded51cac7e065e0f62b1c238b73e Mon Sep 17 00:00:00 2001 From: Youqi Song Date: Sat, 4 May 2024 16:00:46 -0400 Subject: [PATCH 04/10] Reordred the output by PID --- readTree.C | 210 +++++++++++++++++++++++++++++++++-------------------- 1 file changed, 132 insertions(+), 78 deletions(-) diff --git a/readTree.C b/readTree.C index 9085192..e6a273d 100644 --- a/readTree.C +++ b/readTree.C @@ -27,12 +27,12 @@ void readTree(const char *input, const char *output, const char *output_txt) TH1D *h_p = new TH1D("p", "", 15, 0.1, 15.1); TH1D *h_t = new TH1D("theta", "", 20, 2.7, 3.1); - TH1D *h_phi = new TH1D("phi", "", 360, -3.1415926, 3.1415926); + TH1D *h_phi = new TH1D("phi", "", 360, -1.0*TMath::Pi(), TMath::Pi()); - static Int_t fee[15][20][360], fepi[15][20][360], fek[15][20][360], fep[15][20][360], feb[15][20][360]; - static Int_t fpie[15][20][360], fpipi[15][20][360], fpik[15][20][360], fpip[15][20][360], fpib[15][20][360]; - static Int_t fke[15][20][360], fkpi[15][20][360], fkk[15][20][360], fkp[15][20][360], fkb[15][20][360]; - static Int_t fpe[15][20][360], fppi[15][20][360], fpk[15][20][360], fpp[15][20][360], fpb[15][20][360]; + static Float_t fee[15][20][360], fepi[15][20][360], fek[15][20][360], fep[15][20][360], feb[15][20][360]; + static Float_t fpie[15][20][360], fpipi[15][20][360], fpik[15][20][360], fpip[15][20][360], fpib[15][20][360]; + static Float_t fke[15][20][360], fkpi[15][20][360], fkk[15][20][360], fkp[15][20][360], fkb[15][20][360]; + static Float_t fpe[15][20][360], fppi[15][20][360], fpk[15][20][360], fpp[15][20][360], fpb[15][20][360]; for (int pbin = 0; pbin < 15; pbin++) { @@ -40,10 +40,10 @@ void readTree(const char *input, const char *output, const char *output_txt) { for (int phibin = 0; phibin < 360; phibin++) { - fee[pbin][tbin][phibin] = fepi[pbin][tbin][phibin] = fek[pbin][tbin][phibin] = fep[pbin][tbin][phibin] = feb[pbin][tbin][phibin] = 0; - fpie[pbin][tbin][phibin] = fpipi[pbin][tbin][phibin] = fpik[pbin][tbin][phibin] = fpip[pbin][tbin][phibin] = fpib[pbin][tbin][phibin] = 0; - fke[pbin][tbin][phibin] = fkpi[pbin][tbin][phibin] = fkk[pbin][tbin][phibin] = fkp[pbin][tbin][phibin] = fkb[pbin][tbin][phibin] = 0; - fpe[pbin][tbin][phibin] = fppi[pbin][tbin][phibin] = fpk[pbin][tbin][phibin] = fpp[pbin][tbin][phibin] = fpb[pbin][tbin][phibin] = 0; + fee[pbin][tbin][phibin] = fepi[pbin][tbin][phibin] = fek[pbin][tbin][phibin] = fep[pbin][tbin][phibin] = feb[pbin][tbin][phibin] = 0.; + fpie[pbin][tbin][phibin] = fpipi[pbin][tbin][phibin] = fpik[pbin][tbin][phibin] = fpip[pbin][tbin][phibin] = fpib[pbin][tbin][phibin] = 0.; + fke[pbin][tbin][phibin] = fkpi[pbin][tbin][phibin] = fkk[pbin][tbin][phibin] = fkp[pbin][tbin][phibin] = fkb[pbin][tbin][phibin] = 0.; + fpe[pbin][tbin][phibin] = fppi[pbin][tbin][phibin] = fpk[pbin][tbin][phibin] = fpp[pbin][tbin][phibin] = fpb[pbin][tbin][phibin] = 0.; } } } @@ -150,6 +150,81 @@ void readTree(const char *input, const char *output, const char *output_txt) } } + for (int pbin = 0; pbin < 15; pbin++) + { + for (int tbin = 0; tbin < 20; tbin++) + { + for (int phibin = 0; phibin < 360; phibin++) + { + + Double_t fe_tot = fee[pbin][tbin][phibin] + fepi[pbin][tbin][phibin] + fek[pbin][tbin][phibin] + fep[pbin][tbin][phibin] + feb[pbin][tbin][phibin] + zero; + fee[pbin][tbin][phibin] /= fe_tot; + fepi[pbin][tbin][phibin] /= fe_tot; + fek[pbin][tbin][phibin] /= fe_tot; + fep[pbin][tbin][phibin] /= fe_tot; + feb[pbin][tbin][phibin] /= fe_tot; + + Double_t fpi_tot = fpie[pbin][tbin][phibin] + fpipi[pbin][tbin][phibin] + fpik[pbin][tbin][phibin] + fpip[pbin][tbin][phibin] + fpib[pbin][tbin][phibin] + zero; + fpie[pbin][tbin][phibin] /= fpi_tot; + fpipi[pbin][tbin][phibin] /= fpi_tot; + fpik[pbin][tbin][phibin] /= fpi_tot; + fpip[pbin][tbin][phibin] /= fpi_tot; + fpib[pbin][tbin][phibin] /= fpi_tot; + + Double_t fk_tot = fke[pbin][tbin][phibin] + fkpi[pbin][tbin][phibin] + fkk[pbin][tbin][phibin] + fkp[pbin][tbin][phibin] + fkb[pbin][tbin][phibin] + zero; + fke[pbin][tbin][phibin] /= fk_tot; + fkpi[pbin][tbin][phibin] /= fk_tot; + fkk[pbin][tbin][phibin] /= fk_tot; + fkp[pbin][tbin][phibin] /= fk_tot; + fkb[pbin][tbin][phibin] /= fk_tot; + + Double_t fp_tot = fpe[pbin][tbin][phibin] + fppi[pbin][tbin][phibin] + fpk[pbin][tbin][phibin] + fpp[pbin][tbin][phibin] + fpb[pbin][tbin][phibin] + zero; + fpe[pbin][tbin][phibin] /= fp_tot; + fppi[pbin][tbin][phibin] /= fp_tot; + fpk[pbin][tbin][phibin] /= fp_tot; + fpp[pbin][tbin][phibin] /= fp_tot; + fpb[pbin][tbin][phibin] /= fp_tot; + + h[pbin][tbin][phibin]->SetBinContent(1, 4, fee[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(2, 4, fepi[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(3, 4, fek[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(4, 4, fep[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(5, 4, feb[pbin][tbin][phibin]); + + h[pbin][tbin][phibin]->SetBinContent(1, 3, fpie[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(2, 3, fpipi[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(3, 3, fpik[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(4, 3, fpip[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(5, 3, fpib[pbin][tbin][phibin]); + + h[pbin][tbin][phibin]->SetBinContent(1, 2, fke[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(2, 2, fkpi[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(3, 2, fkk[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(4, 2, fkp[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(5, 2, fkb[pbin][tbin][phibin]); + + h[pbin][tbin][phibin]->SetBinContent(1, 1, fpe[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(2, 1, fppi[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(3, 1, fpk[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(4, 1, fpp[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(5, 1, fpb[pbin][tbin][phibin]); + + Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1); + Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1); + Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); + + if (pbin == 0 && tbin == 0 && phibin == 0) + { + cout << "truth PID, charge, p, theta, phi, prob(e), prob(pi), prob(k), prob(p), prob(fail)" << endl; + cout << "11 -1 " << p << " " << t << " " << phi << " " << fee[pbin][tbin][phibin] << " " << fepi[pbin][tbin][phibin] << " " << fek[pbin][tbin][phibin] << " " << fep[pbin][tbin][phibin] << " " << feb[pbin][tbin][phibin] << endl; + cout << "211 1 " << p << " " << t << " " << phi << " " << fpie[pbin][tbin][phibin] << " " << fpipi[pbin][tbin][phibin] << " " << fpik[pbin][tbin][phibin] << " " << fpip[pbin][tbin][phibin] << " " << fpib[pbin][tbin][phibin] << endl; + cout << "321 1 " << p << " " << t << " " << phi << " " << fke[pbin][tbin][phibin] << " " << fkpi[pbin][tbin][phibin] << " " << fkk[pbin][tbin][phibin] << " " << fkp[pbin][tbin][phibin] << " " << fkb[pbin][tbin][phibin] << endl; + cout << "2212 1 " << p << " " << t << " " << phi << " " << fpe[pbin][tbin][phibin] << " " << fppi[pbin][tbin][phibin] << " " << fpk[pbin][tbin][phibin] << " " << fpp[pbin][tbin][phibin] << " " << fpb[pbin][tbin][phibin] << endl; + } + } + } + } + FILE *filePointer; filePointer = fopen(output_txt, "w"); @@ -159,84 +234,63 @@ void readTree(const char *input, const char *output, const char *output_txt) { for (int phibin = 0; phibin < 360; phibin++) { + Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1); + Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1); + Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); + fprintf(filePointer, "11 -1 %f %f %f %f %f %f %f %f\n", p, t, phi, fee[pbin][tbin][phibin], fepi[pbin][tbin][phibin], fek[pbin][tbin][phibin], fep[pbin][tbin][phibin], feb[pbin][tbin][phibin]); + + } + } + } + + for (int pbin = 0; pbin < 15; pbin++) + { + for (int tbin = 0; tbin < 20; tbin++) + { + for (int phibin = 0; phibin < 360; phibin++) + { Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1); Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1); Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); - Double_t mc_e_rec_e = (1.0 * fee[pbin][tbin][phibin]) / (1.0 * (fee[pbin][tbin][phibin] + fepi[pbin][tbin][phibin] + fek[pbin][tbin][phibin] + fep[pbin][tbin][phibin] + feb[pbin][tbin][phibin] + zero)); - Double_t mc_e_rec_pi = (1.0 * fepi[pbin][tbin][phibin]) / (1.0 * (fee[pbin][tbin][phibin] + fepi[pbin][tbin][phibin] + fek[pbin][tbin][phibin] + fep[pbin][tbin][phibin] + feb[pbin][tbin][phibin] + zero)); - Double_t mc_e_rec_k = (1.0 * fek[pbin][tbin][phibin]) / (1.0 * (fee[pbin][tbin][phibin] + fepi[pbin][tbin][phibin] + fek[pbin][tbin][phibin] + fep[pbin][tbin][phibin] + feb[pbin][tbin][phibin] + zero)); - Double_t mc_e_rec_p = (1.0 * fep[pbin][tbin][phibin]) / (1.0 * (fee[pbin][tbin][phibin] + fepi[pbin][tbin][phibin] + fek[pbin][tbin][phibin] + fep[pbin][tbin][phibin] + feb[pbin][tbin][phibin] + zero)); - Double_t mc_e_rec_b = (1.0 * feb[pbin][tbin][phibin]) / (1.0 * (fee[pbin][tbin][phibin] + fepi[pbin][tbin][phibin] + fek[pbin][tbin][phibin] + fep[pbin][tbin][phibin] + feb[pbin][tbin][phibin] + zero)); - - Double_t mc_pi_rec_e = (1.0 * fpie[pbin][tbin][phibin]) / (1.0 * (fpie[pbin][tbin][phibin] + fpipi[pbin][tbin][phibin] + fpik[pbin][tbin][phibin] + fpip[pbin][tbin][phibin] + fpib[pbin][tbin][phibin] + zero)); - Double_t mc_pi_rec_pi = (1.0 * fpipi[pbin][tbin][phibin]) / (1.0 * (fpie[pbin][tbin][phibin] + fpipi[pbin][tbin][phibin] + fpik[pbin][tbin][phibin] + fpip[pbin][tbin][phibin] + fpib[pbin][tbin][phibin] + zero)); - Double_t mc_pi_rec_k = (1.0 * fpik[pbin][tbin][phibin]) / (1.0 * (fpie[pbin][tbin][phibin] + fpipi[pbin][tbin][phibin] + fpik[pbin][tbin][phibin] + fpip[pbin][tbin][phibin] + fpib[pbin][tbin][phibin] + zero)); - Double_t mc_pi_rec_p = (1.0 * fpip[pbin][tbin][phibin]) / (1.0 * (fpie[pbin][tbin][phibin] + fpipi[pbin][tbin][phibin] + fpik[pbin][tbin][phibin] + fpip[pbin][tbin][phibin] + fpib[pbin][tbin][phibin] + zero)); - Double_t mc_pi_rec_b = (1.0 * fpib[pbin][tbin][phibin]) / (1.0 * (fpie[pbin][tbin][phibin] + fpipi[pbin][tbin][phibin] + fpik[pbin][tbin][phibin] + fpip[pbin][tbin][phibin] + fpib[pbin][tbin][phibin] + zero)); - - Double_t mc_k_rec_e = (1.0 * fke[pbin][tbin][phibin]) / (1.0 * (fke[pbin][tbin][phibin] + fkpi[pbin][tbin][phibin] + fkk[pbin][tbin][phibin] + fkp[pbin][tbin][phibin] + fkb[pbin][tbin][phibin] + zero)); - Double_t mc_k_rec_pi = (1.0 * fkpi[pbin][tbin][phibin]) / (1.0 * (fke[pbin][tbin][phibin] + fkpi[pbin][tbin][phibin] + fkk[pbin][tbin][phibin] + fkp[pbin][tbin][phibin] + fkb[pbin][tbin][phibin] + zero)); - Double_t mc_k_rec_k = (1.0 * fkk[pbin][tbin][phibin]) / (1.0 * (fke[pbin][tbin][phibin] + fkpi[pbin][tbin][phibin] + fkk[pbin][tbin][phibin] + fkp[pbin][tbin][phibin] + fkb[pbin][tbin][phibin] + zero)); - Double_t mc_k_rec_p = (1.0 * fkp[pbin][tbin][phibin]) / (1.0 * (fke[pbin][tbin][phibin] + fkpi[pbin][tbin][phibin] + fkk[pbin][tbin][phibin] + fkp[pbin][tbin][phibin] + fkb[pbin][tbin][phibin] + zero)); - Double_t mc_k_rec_b = (1.0 * fkb[pbin][tbin][phibin]) / (1.0 * (fke[pbin][tbin][phibin] + fkpi[pbin][tbin][phibin] + fkk[pbin][tbin][phibin] + fkp[pbin][tbin][phibin] + fkb[pbin][tbin][phibin] + zero)); - - Double_t mc_p_rec_e = (1.0 * fpe[pbin][tbin][phibin]) / (1.0 * (fpe[pbin][tbin][phibin] + fppi[pbin][tbin][phibin] + fpk[pbin][tbin][phibin] + fpp[pbin][tbin][phibin] + fpb[pbin][tbin][phibin] + zero)); - Double_t mc_p_rec_pi = (1.0 * fppi[pbin][tbin][phibin]) / (1.0 * (fpe[pbin][tbin][phibin] + fppi[pbin][tbin][phibin] + fpk[pbin][tbin][phibin] + fpp[pbin][tbin][phibin] + fpb[pbin][tbin][phibin] + zero)); - Double_t mc_p_rec_k = (1.0 * fpk[pbin][tbin][phibin]) / (1.0 * (fpe[pbin][tbin][phibin] + fppi[pbin][tbin][phibin] + fpk[pbin][tbin][phibin] + fpp[pbin][tbin][phibin] + fpb[pbin][tbin][phibin] + zero)); - Double_t mc_p_rec_p = (1.0 * fpp[pbin][tbin][phibin]) / (1.0 * (fpe[pbin][tbin][phibin] + fppi[pbin][tbin][phibin] + fpk[pbin][tbin][phibin] + fpp[pbin][tbin][phibin] + fpb[pbin][tbin][phibin] + zero)); - Double_t mc_p_rec_b = (1.0 * fpb[pbin][tbin][phibin]) / (1.0 * (fpe[pbin][tbin][phibin] + fppi[pbin][tbin][phibin] + fpk[pbin][tbin][phibin] + fpp[pbin][tbin][phibin] + fpb[pbin][tbin][phibin] + zero)); - - h[pbin][tbin][phibin]->SetBinContent(1, 4, mc_e_rec_e); - h[pbin][tbin][phibin]->SetBinContent(2, 4, mc_e_rec_pi); - h[pbin][tbin][phibin]->SetBinContent(3, 4, mc_e_rec_k); - h[pbin][tbin][phibin]->SetBinContent(4, 4, mc_e_rec_p); - h[pbin][tbin][phibin]->SetBinContent(5, 4, mc_e_rec_b); - - h[pbin][tbin][phibin]->SetBinContent(1, 3, mc_pi_rec_e); - h[pbin][tbin][phibin]->SetBinContent(2, 3, mc_pi_rec_pi); - h[pbin][tbin][phibin]->SetBinContent(3, 3, mc_pi_rec_k); - h[pbin][tbin][phibin]->SetBinContent(4, 3, mc_pi_rec_p); - h[pbin][tbin][phibin]->SetBinContent(5, 3, mc_pi_rec_b); - - h[pbin][tbin][phibin]->SetBinContent(1, 2, mc_k_rec_e); - h[pbin][tbin][phibin]->SetBinContent(2, 2, mc_k_rec_pi); - h[pbin][tbin][phibin]->SetBinContent(3, 2, mc_k_rec_k); - h[pbin][tbin][phibin]->SetBinContent(4, 2, mc_k_rec_p); - h[pbin][tbin][phibin]->SetBinContent(5, 2, mc_k_rec_b); - - h[pbin][tbin][phibin]->SetBinContent(1, 1, mc_p_rec_e); - h[pbin][tbin][phibin]->SetBinContent(2, 1, mc_p_rec_pi); - h[pbin][tbin][phibin]->SetBinContent(3, 1, mc_p_rec_k); - h[pbin][tbin][phibin]->SetBinContent(4, 1, mc_p_rec_p); - h[pbin][tbin][phibin]->SetBinContent(5, 1, mc_p_rec_b); - - // h[pbin][tbin][phibin]->Write(); - - fprintf(filePointer, "11 %f %f %f %f %f %f %f %f\n", p, t, phi, mc_e_rec_e, mc_e_rec_pi, mc_e_rec_k, mc_e_rec_p, mc_e_rec_b); - fprintf(filePointer, "211 %f %f %f %f %f %f %f %f\n", p, t, phi, mc_pi_rec_e, mc_pi_rec_pi, mc_pi_rec_k, mc_pi_rec_p, mc_pi_rec_b); - fprintf(filePointer, "321 %f %f %f %f %f %f %f %f\n", p, t, phi, mc_k_rec_e, mc_k_rec_pi, mc_k_rec_k, mc_k_rec_p, mc_k_rec_b); - fprintf(filePointer, "2212 %f %f %f %f %f %f %f %f\n", p, t, phi, mc_p_rec_e, mc_p_rec_pi, mc_p_rec_k, mc_p_rec_p, mc_p_rec_b); + fprintf(filePointer, "211 1 %f %f %f %f %f %f %f %f\n", p, t, phi, fpie[pbin][tbin][phibin], fpipi[pbin][tbin][phibin], fpik[pbin][tbin][phibin], fpip[pbin][tbin][phibin], fpib[pbin][tbin][phibin]); + + } + } + } + + for (int pbin = 0; pbin < 15; pbin++) + { + for (int tbin = 0; tbin < 20; tbin++) + { + for (int phibin = 0; phibin < 360; phibin++) + { + Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1); + Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1); + Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); - if (pbin == 0 && tbin == 0 && phibin == 0) - { - cout << "truth PID, p, theta, phi, prob(e), prob(pi), prob(k), prob(p), prob(fail)" << endl; - cout << "11 " << p << " " << t << " " << phi << " " << mc_e_rec_e << " " << mc_e_rec_pi << " " << mc_e_rec_k << " " << mc_e_rec_p << " " << mc_e_rec_b << endl; - cout << "211 " << p << " " << t << " " << phi << " " << mc_pi_rec_e << " " << mc_pi_rec_pi << " " << mc_pi_rec_k << " " << mc_pi_rec_p << " " << mc_pi_rec_b << endl; - cout << "321 " << p << " " << t << " " << phi << " " << mc_k_rec_e << " " << mc_k_rec_pi << " " << mc_k_rec_k << " " << mc_k_rec_p << " " << mc_k_rec_b << endl; - cout << "2212 " << p << " " << t << " " << phi << " " << mc_p_rec_e << " " << mc_p_rec_pi << " " << mc_p_rec_k << " " << mc_p_rec_p << " " << mc_p_rec_b << endl; - } - if (pbin == 14 && tbin == 9 && phibin == 11) - { - cout << "11 " << p << " " << t << " " << phi << " " << mc_e_rec_e << " " << mc_e_rec_pi << " " << mc_e_rec_k << " " << mc_e_rec_p << " " << mc_e_rec_b << endl; - cout << "211 " << p << " " << t << " " << phi << " " << mc_pi_rec_e << " " << mc_pi_rec_pi << " " << mc_pi_rec_k << " " << mc_pi_rec_p << " " << mc_pi_rec_b << endl; - cout << "321 " << p << " " << t << " " << phi << " " << mc_k_rec_e << " " << mc_k_rec_pi << " " << mc_k_rec_k << " " << mc_k_rec_p << " " << mc_k_rec_b << endl; - cout << "2212 " << p << " " << t << " " << phi << " " << mc_p_rec_e << " " << mc_p_rec_pi << " " << mc_p_rec_k << " " << mc_p_rec_p << " " << mc_p_rec_b << endl; - } + fprintf(filePointer, "321 1 %f %f %f %f %f %f %f %f\n", p, t, phi, fke[pbin][tbin][phibin], fkpi[pbin][tbin][phibin], fkk[pbin][tbin][phibin], fkp[pbin][tbin][phibin], fkb[pbin][tbin][phibin]); + } } } + + for (int pbin = 0; pbin < 15; pbin++) + { + for (int tbin = 0; tbin < 20; tbin++) + { + for (int phibin = 0; phibin < 360; phibin++) + { + Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1); + Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1); + Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); + + fprintf(filePointer, "2212 1 %f %f %f %f %f %f %f %f\n", p, t, phi, fpe[pbin][tbin][phibin], fppi[pbin][tbin][phibin], fpk[pbin][tbin][phibin], fpp[pbin][tbin][phibin], fpb[pbin][tbin][phibin]); + + } + } + } fclose(filePointer); ofile->Write(); From 18f328046bc065ca5344ecc9cedc57af9a5eed32 Mon Sep 17 00:00:00 2001 From: Youqi Song Date: Sat, 4 May 2024 22:34:13 -0400 Subject: [PATCH 05/10] Added new plots for QA --- readTree.C | 14 +- readTreeQA.C | 708 ++++++++++++++++++++++++++++----------------------- 2 files changed, 391 insertions(+), 331 deletions(-) diff --git a/readTree.C b/readTree.C index e6a273d..8f9c221 100644 --- a/readTree.C +++ b/readTree.C @@ -27,7 +27,7 @@ void readTree(const char *input, const char *output, const char *output_txt) TH1D *h_p = new TH1D("p", "", 15, 0.1, 15.1); TH1D *h_t = new TH1D("theta", "", 20, 2.7, 3.1); - TH1D *h_phi = new TH1D("phi", "", 360, -1.0*TMath::Pi(), TMath::Pi()); + TH1D *h_phi = new TH1D("phi", "", 360, -1.0 * TMath::Pi(), TMath::Pi()); static Float_t fee[15][20][360], fepi[15][20][360], fek[15][20][360], fep[15][20][360], feb[15][20][360]; static Float_t fpie[15][20][360], fpipi[15][20][360], fpik[15][20][360], fpip[15][20][360], fpib[15][20][360]; @@ -212,7 +212,7 @@ void readTree(const char *input, const char *output, const char *output_txt) Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1); Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1); Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); - + if (pbin == 0 && tbin == 0 && phibin == 0) { cout << "truth PID, charge, p, theta, phi, prob(e), prob(pi), prob(k), prob(p), prob(fail)" << endl; @@ -239,7 +239,6 @@ void readTree(const char *input, const char *output, const char *output_txt) Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); fprintf(filePointer, "11 -1 %f %f %f %f %f %f %f %f\n", p, t, phi, fee[pbin][tbin][phibin], fepi[pbin][tbin][phibin], fek[pbin][tbin][phibin], fep[pbin][tbin][phibin], feb[pbin][tbin][phibin]); - } } } @@ -255,11 +254,10 @@ void readTree(const char *input, const char *output, const char *output_txt) Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); fprintf(filePointer, "211 1 %f %f %f %f %f %f %f %f\n", p, t, phi, fpie[pbin][tbin][phibin], fpipi[pbin][tbin][phibin], fpik[pbin][tbin][phibin], fpip[pbin][tbin][phibin], fpib[pbin][tbin][phibin]); - } } } - + for (int pbin = 0; pbin < 15; pbin++) { for (int tbin = 0; tbin < 20; tbin++) @@ -271,11 +269,10 @@ void readTree(const char *input, const char *output, const char *output_txt) Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); fprintf(filePointer, "321 1 %f %f %f %f %f %f %f %f\n", p, t, phi, fke[pbin][tbin][phibin], fkpi[pbin][tbin][phibin], fkk[pbin][tbin][phibin], fkp[pbin][tbin][phibin], fkb[pbin][tbin][phibin]); - } } } - + for (int pbin = 0; pbin < 15; pbin++) { for (int tbin = 0; tbin < 20; tbin++) @@ -287,10 +284,9 @@ void readTree(const char *input, const char *output, const char *output_txt) Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); fprintf(filePointer, "2212 1 %f %f %f %f %f %f %f %f\n", p, t, phi, fpe[pbin][tbin][phibin], fppi[pbin][tbin][phibin], fpk[pbin][tbin][phibin], fpp[pbin][tbin][phibin], fpb[pbin][tbin][phibin]); - } } - } + } fclose(filePointer); ofile->Write(); diff --git a/readTreeQA.C b/readTreeQA.C index b561200..b34ca6e 100644 --- a/readTreeQA.C +++ b/readTreeQA.C @@ -1,347 +1,411 @@ void readTreeQA(const char *input, const char *output, const char *output_txt) { - TFile *fa = new TFile(input); - - TTree *ta = (TTree *)fa->Get("T"); - - TFile *ofile = TFile::Open(output, "recreate"); - - Int_t partPDG, recoPDG, recoPhotons, recoCherenkovHitCount; - Double_t partMom, partEta, partTheta, partPhi, partVertX, partVertY, partVertZ, recoTrackCherenkovAngle; - Double_t zero = pow(10., -50.); - - ta->SetBranchAddress("partPDG", &partPDG); - ta->SetBranchAddress("partMom", &partMom); - ta->SetBranchAddress("partEta", &partEta); - ta->SetBranchAddress("partTheta", &partTheta); - ta->SetBranchAddress("partPhi", &partPhi); - ta->SetBranchAddress("partVertX", &partVertX); - ta->SetBranchAddress("partVertY", &partVertY); - ta->SetBranchAddress("partVertZ", &partVertZ); - ta->SetBranchAddress("recoPDG", &recoPDG); - ta->SetBranchAddress("recoPhotons", &recoPhotons); - ta->SetBranchAddress("recoTrackCherenkovAngle", &recoTrackCherenkovAngle); - - TH1D *h_p = new TH1D("p", "", 15, 0.1, 15.1); - TH1D *h_t = new TH1D("theta", "", 16, 2.7, 3.1); - TH1D *h_phi = new TH1D("phi", "", 360, -3.1415926, 3.1415926); - - TProfile2D *recoPhiVsEtaElecPhoton[4]; - TH2D *partPhiVsEtaElec[4]; - TH2D *partPhiVsEtaElecMatch[4]; - TH2D *partPhiVsEtaElecUnMatch[4]; - TH2D *partPhiVsEtaElecMatchPhoton[4]; - - TProfile2D *recoPhiVsEtaPionPhoton[4]; - TH2D *partPhiVsEtaPion[4]; - TH2D *partPhiVsEtaPionMatch[4]; - TH2D *partPhiVsEtaPionUnMatch[4]; - - TProfile2D *recoPhiVsEtaKaonPhoton[4]; - TH2D *partPhiVsEtaKaon[4]; - TH2D *partPhiVsEtaKaonMatch[4]; - TH2D *partPhiVsEtaKaonUnMatch[4]; - - TProfile2D *recoPhiVsEtaProtonPhoton[4]; - TH2D *partPhiVsEtaProton[4]; - TH2D *partPhiVsEtaProtonMatch[4]; - TH2D *partPhiVsEtaProtonUnMatch[4]; - - for (int i = 0; i < 4; i++) - { - recoPhiVsEtaElecPhoton[i] = new TProfile2D(Form("recoPhiVsEtaElecPhoton_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); - partPhiVsEtaElec[i] = new TH2D(Form("partPhiVsEtaElecMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); - partPhiVsEtaElecMatch[i] = new TH2D(Form("partPhiVsEtaElecMatchMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); - partPhiVsEtaElecUnMatch[i] = new TH2D(Form("partPhiVsEtaElecUnMatchMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); - partPhiVsEtaElecMatchPhoton[i] = new TH2D(Form("partPhiVsEtaElecMatchPhotonMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); - - recoPhiVsEtaPionPhoton[i] = new TProfile2D(Form("recoPhiVsEtaPionPhoton_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); - partPhiVsEtaPion[i] = new TH2D(Form("partPhiVsEtaPionMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); - partPhiVsEtaPionMatch[i] = new TH2D(Form("partPhiVsEtaPionMatchMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); - partPhiVsEtaPionUnMatch[i] = new TH2D(Form("partPhiVsEtaPionUnMatchMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); - - recoPhiVsEtaKaonPhoton[i] = new TProfile2D(Form("recoPhiVsEtaKaonPhoton_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); - partPhiVsEtaKaon[i] = new TH2D(Form("partPhiVsEtaKaonMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); - partPhiVsEtaKaonMatch[i] = new TH2D(Form("partPhiVsEtaKaonMatchMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); - partPhiVsEtaKaonUnMatch[i] = new TH2D(Form("partPhiVsEtaKaonUnMatchMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); - - recoPhiVsEtaProtonPhoton[i] = new TProfile2D(Form("recoPhiVsEtaProtonPhoton_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); - partPhiVsEtaProton[i] = new TH2D(Form("partPhiVsEtaProtonMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); - partPhiVsEtaProtonMatch[i] = new TH2D(Form("partPhiVsEtaProtonMatchMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); - partPhiVsEtaProtonUnMatch[i] = new TH2D(Form("partPhiVsEtaProtonUnMatchMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); - } - - TH2D *numPhotonsVsPhiEta25Mom10Pion = new TH2D("numPhotonsVsPhiEta25Mom10Pion", "", 60, -TMath::Pi(), TMath::Pi(), 25, 0., 25.); - TH2D *numPhotonsVsPhiEta25Mom4Pion = new TH2D("numPhotonsVsPhiEta25Mom4Pion", "", 60, -TMath::Pi(), TMath::Pi(), 25, 0., 25.); - - Int_t fee[15][16][360], fepi[15][16][360], fek[15][16][360], fep[15][16][360], feb[15][16][360]; - Int_t fpie[15][16][360], fpipi[15][16][360], fpik[15][16][360], fpip[15][16][360], fpib[15][16][360]; - Int_t fke[15][16][360], fkpi[15][16][360], fkk[15][16][360], fkp[15][16][360], fkb[15][16][360]; - Int_t fpe[15][16][360], fppi[15][16][360], fpk[15][16][360], fpp[15][16][360], fpb[15][16][360]; - - for (int pbin = 0; pbin < 15; pbin++) - { - for (int tbin = 0; tbin < 16; tbin++) + TFile *fa = new TFile(input); + + TTree *ta = (TTree *)fa->Get("T"); + + TFile *ofile = TFile::Open(output, "recreate"); + + Int_t partPDG, recoPDG, recoPhotons, recoCherenkovHitCount; + Double_t partMom, partEta, partTheta, partPhi, partVertX, partVertY, partVertZ, recoTrackCherenkovAngle; + Double_t zero = pow(10., -50.); + + ta->SetBranchAddress("partPDG", &partPDG); + ta->SetBranchAddress("partMom", &partMom); + ta->SetBranchAddress("partEta", &partEta); + ta->SetBranchAddress("partTheta", &partTheta); + ta->SetBranchAddress("partPhi", &partPhi); + ta->SetBranchAddress("partVertX", &partVertX); + ta->SetBranchAddress("partVertY", &partVertY); + ta->SetBranchAddress("partVertZ", &partVertZ); + ta->SetBranchAddress("recoPDG", &recoPDG); + ta->SetBranchAddress("recoPhotons", &recoPhotons); + ta->SetBranchAddress("recoTrackCherenkovAngle", &recoTrackCherenkovAngle); + + TH1D *h_p = new TH1D("p", "", 15, 0.1, 15.1); + TH1D *h_t = new TH1D("theta", "", 10, 2.7, 3.1); + TH1D *h_phi = new TH1D("phi", "", 12, -1.0 * TMath::Pi(), TMath::Pi()); + + TProfile2D *recoPhiVsEtaElecPhoton[4]; + TH2D *partPhiVsEtaElec[4]; + TH2D *partPhiVsEtaElecMatch[4]; + TH2D *partPhiVsEtaElecUnMatch[4]; + TH2D *partPhiVsEtaElecMatchPhoton[4]; + TH1D *recoElecPhoton[4]; + + TProfile2D *recoPhiVsEtaPionPhoton[4]; + TH2D *partPhiVsEtaPion[4]; + TH2D *partPhiVsEtaPionMatch[4]; + TH2D *partPhiVsEtaPionUnMatch[4]; + TH1D *recoPionPhoton[4]; + + TProfile2D *recoPhiVsEtaKaonPhoton[4]; + TH2D *partPhiVsEtaKaon[4]; + TH2D *partPhiVsEtaKaonMatch[4]; + TH2D *partPhiVsEtaKaonUnMatch[4]; + TH1D *recoKaonPhoton[4]; + + TProfile2D *recoPhiVsEtaProtonPhoton[4]; + TH2D *partPhiVsEtaProton[4]; + TH2D *partPhiVsEtaProtonMatch[4]; + TH2D *partPhiVsEtaProtonUnMatch[4]; + TH1D *recoProtonPhoton[4]; + + for (int i = 0; i < 4; i++) + { + recoPhiVsEtaElecPhoton[i] = new TProfile2D(Form("recoPhiVsEtaElecPhoton_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + partPhiVsEtaElec[i] = new TH2D(Form("partPhiVsEtaElecMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + partPhiVsEtaElecMatch[i] = new TH2D(Form("partPhiVsEtaElecMatchMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + partPhiVsEtaElecUnMatch[i] = new TH2D(Form("partPhiVsEtaElecUnMatchMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + partPhiVsEtaElecMatchPhoton[i] = new TH2D(Form("partPhiVsEtaElecMatchPhotonMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + recoElecPhoton[i] = new TH1D(Form("recoElecPhoton_%d", i), "", 25, -0.1, 24.9); + + recoPhiVsEtaPionPhoton[i] = new TProfile2D(Form("recoPhiVsEtaPionPhoton_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + partPhiVsEtaPion[i] = new TH2D(Form("partPhiVsEtaPionMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + partPhiVsEtaPionMatch[i] = new TH2D(Form("partPhiVsEtaPionMatchMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + partPhiVsEtaPionUnMatch[i] = new TH2D(Form("partPhiVsEtaPionUnMatchMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + recoPionPhoton[i] = new TH1D(Form("recoPionPhoton_%d", i), "", 25, -0.1, 24.9); + + recoPhiVsEtaKaonPhoton[i] = new TProfile2D(Form("recoPhiVsEtaKaonPhoton_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + partPhiVsEtaKaon[i] = new TH2D(Form("partPhiVsEtaKaonMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + partPhiVsEtaKaonMatch[i] = new TH2D(Form("partPhiVsEtaKaonMatchMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + partPhiVsEtaKaonUnMatch[i] = new TH2D(Form("partPhiVsEtaKaonUnMatchMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + recoKaonPhoton[i] = new TH1D(Form("recoKaonPhoton_%d", i), "", 25, -0.1, 24.9); + + recoPhiVsEtaProtonPhoton[i] = new TProfile2D(Form("recoPhiVsEtaProtonPhoton_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + partPhiVsEtaProton[i] = new TH2D(Form("partPhiVsEtaProtonMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + partPhiVsEtaProtonMatch[i] = new TH2D(Form("partPhiVsEtaProtonMatchMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + partPhiVsEtaProtonUnMatch[i] = new TH2D(Form("partPhiVsEtaProtonUnMatchMom_%d", i), "", 50, -4., -1., 100, -TMath::Pi(), TMath::Pi()); + recoProtonPhoton[i] = new TH1D(Form("recoProtonPhoton_%d", i), "", 25, -0.1, 24.9); + } + + TH2D *numPhotonsVsPhiEta25Mom10Pion = new TH2D("numPhotonsVsPhiEta25Mom10Pion", "", 60, -TMath::Pi(), TMath::Pi(), 25, 0., 25.); + TH2D *numPhotonsVsPhiEta25Mom4Pion = new TH2D("numPhotonsVsPhiEta25Mom4Pion", "", 60, -TMath::Pi(), TMath::Pi(), 25, 0., 25.); + + static Float_t fee[15][10][12], fepi[15][10][12], fek[15][10][12], fep[15][10][12], feb[15][10][12]; + static Float_t fpie[15][10][12], fpipi[15][10][12], fpik[15][10][12], fpip[15][10][12], fpib[15][10][12]; + static Float_t fke[15][10][12], fkpi[15][10][12], fkk[15][10][12], fkp[15][10][12], fkb[15][10][12]; + static Float_t fpe[15][10][12], fppi[15][10][12], fpk[15][10][12], fpp[15][10][12], fpb[15][10][12]; + + for (int pbin = 0; pbin < 15; pbin++) + { + for (int tbin = 0; tbin < 10; tbin++) + { + for (int phibin = 0; phibin < 12; phibin++) { - for (int phibin = 0; phibin < 360; phibin++) - { - fee[pbin][tbin][phibin] = fepi[pbin][tbin][phibin] = fek[pbin][tbin][phibin] = fep[pbin][tbin][phibin] = feb[pbin][tbin][phibin] = 0; - fpie[pbin][tbin][phibin] = fpipi[pbin][tbin][phibin] = fpik[pbin][tbin][phibin] = fpip[pbin][tbin][phibin] = fpib[pbin][tbin][phibin] = 0; - fke[pbin][tbin][phibin] = fkpi[pbin][tbin][phibin] = fkk[pbin][tbin][phibin] = fkp[pbin][tbin][phibin] = fkb[pbin][tbin][phibin] = 0; - fpe[pbin][tbin][phibin] = fppi[pbin][tbin][phibin] = fpk[pbin][tbin][phibin] = fpp[pbin][tbin][phibin] = fpb[pbin][tbin][phibin] = 0; - } + fee[pbin][tbin][phibin] = fepi[pbin][tbin][phibin] = fek[pbin][tbin][phibin] = fep[pbin][tbin][phibin] = feb[pbin][tbin][phibin] = 0; + fpie[pbin][tbin][phibin] = fpipi[pbin][tbin][phibin] = fpik[pbin][tbin][phibin] = fpip[pbin][tbin][phibin] = fpib[pbin][tbin][phibin] = 0; + fke[pbin][tbin][phibin] = fkpi[pbin][tbin][phibin] = fkk[pbin][tbin][phibin] = fkp[pbin][tbin][phibin] = fkb[pbin][tbin][phibin] = 0; + fpe[pbin][tbin][phibin] = fppi[pbin][tbin][phibin] = fpk[pbin][tbin][phibin] = fpp[pbin][tbin][phibin] = fpb[pbin][tbin][phibin] = 0; } - } - - Int_t nEntries = (Int_t)ta->GetEntries(); - for (int i = 0; i < nEntries; i++) - { - if (i % 100000 == 0) - cout << "Processing Event " << i << endl; - ta->GetEntry(i); - - Int_t pbin = h_p->FindBin(partMom) - 1; - Int_t tbin = h_t->FindBin(partTheta) - 1; - Int_t phibin = h_phi->FindBin(partPhi) - 1; - - Int_t momSelect = -1; - if (partMom > 1.0 && partMom < 3.0) - momSelect = 0; - if (partMom > 3.0 && partMom < 5.0) - momSelect = 1; - if (partMom > 5.0 && partMom < 7.0) - momSelect = 2; - if (partMom > 7.0 && partMom < 9.0) - momSelect = 3; - - // if (partMom < 1) cout << pbin << endl;// -> bin number starts from 1 - - if (partPDG == 11) // Look at Thrown Electrons + } + } + + Int_t nEntries = (Int_t)ta->GetEntries(); + for (int i = 0; i < nEntries; i++) + { + if (i % 100000 == 0) + cout << "Processing Event " << i << endl; + ta->GetEntry(i); + + Int_t pbin = h_p->FindBin(partMom) - 1; + Int_t tbin = h_t->FindBin(partTheta) - 1; + Int_t phibin = h_phi->FindBin(partPhi) - 1; + + Int_t momSelect = -1; + if (partMom > 1.0 && partMom < 3.0) + momSelect = 0; + if (partMom > 3.0 && partMom < 5.0) + momSelect = 1; + if (partMom > 5.0 && partMom < 7.0) + momSelect = 2; + if (partMom > 7.0 && partMom < 9.0) + momSelect = 3; + + // if (partMom < 1) cout << pbin << endl;// -> bin number starts from 1 + + if (partPDG == 11) // Look at Thrown Electrons + { + if (recoPhotons < 4) { - if (recoPhotons < 4) - { - feb[pbin][tbin][phibin]++; - } - else - { - if (recoPDG == 11) - fee[pbin][tbin][phibin]++; - if (recoPDG == 211) - fepi[pbin][tbin][phibin]++; - if (recoPDG == 321) - fek[pbin][tbin][phibin]++; - if (recoPDG == 2212) - fep[pbin][tbin][phibin]++; - } - if (momSelect != -1) - { - recoPhiVsEtaElecPhoton[momSelect]->Fill(partEta, partPhi, recoPhotons); - partPhiVsEtaElec[momSelect]->Fill(partEta, partPhi); - if (recoPDG == 11) - partPhiVsEtaElecMatch[momSelect]->Fill(partEta, partPhi); - if (recoPDG != 11) - partPhiVsEtaElecUnMatch[momSelect]->Fill(partEta, partPhi); - if (recoPDG == 11 && recoPhotons > 0) - partPhiVsEtaElecMatchPhoton[momSelect]->Fill(partEta, partPhi); - } + feb[pbin][tbin][phibin]++; } + else + { + if (recoPDG == 11) + fee[pbin][tbin][phibin]++; + if (recoPDG == 211) + fepi[pbin][tbin][phibin]++; + if (recoPDG == 321) + fek[pbin][tbin][phibin]++; + if (recoPDG == 2212) + fep[pbin][tbin][phibin]++; + } + if (momSelect != -1 && partEta > -2.6 && partEta < -2.4 && partPhi > 2 && partPhi < 2.5) + { + recoElecPhoton[momSelect]->Fill(recoPhotons); + recoPhiVsEtaElecPhoton[momSelect]->Fill(partEta, partPhi, recoPhotons); + partPhiVsEtaElec[momSelect]->Fill(partEta, partPhi); + if (recoPDG == 11) + partPhiVsEtaElecMatch[momSelect]->Fill(partEta, partPhi); + if (recoPDG != 11) + partPhiVsEtaElecUnMatch[momSelect]->Fill(partEta, partPhi); + if (recoPDG == 11 && recoPhotons > 0) + partPhiVsEtaElecMatchPhoton[momSelect]->Fill(partEta, partPhi); + } + } - if (partPDG == 211) // Look at Thrown Pions + if (partPDG == 211) // Look at Thrown Pions + { + if (recoPhotons < 4) + { + fpib[pbin][tbin][phibin]++; + } + else + { + if (recoPDG == 11) + fpie[pbin][tbin][phibin]++; + if (recoPDG == 211) + fpipi[pbin][tbin][phibin]++; + if (recoPDG == 321) + fpik[pbin][tbin][phibin]++; + if (recoPDG == 2212) + fpip[pbin][tbin][phibin]++; + } + if (momSelect != -1 && partEta > -2.6 && partEta < -2.4 && partPhi > 2 && partPhi < 2.5) { - if (recoPhotons < 4) - { - fpib[pbin][tbin][phibin]++; - } - else - { - if (recoPDG == 11) - fpie[pbin][tbin][phibin]++; - if (recoPDG == 211) - fpipi[pbin][tbin][phibin]++; - if (recoPDG == 321) - fpik[pbin][tbin][phibin]++; - if (recoPDG == 2212) - fpip[pbin][tbin][phibin]++; - } - if (momSelect != -1) - { - recoPhiVsEtaPionPhoton[momSelect]->Fill(partEta, partPhi, recoPhotons); - partPhiVsEtaPion[momSelect]->Fill(partEta, partPhi); - if (recoPDG == 211) - partPhiVsEtaPionMatch[momSelect]->Fill(partEta, partPhi); - if (recoPDG != 211) - partPhiVsEtaPionUnMatch[momSelect]->Fill(partEta, partPhi); - } - if (partEta < -2.49 && partEta > -2.51) - { - if (partMom > 9.5 && partMom < 10.5) - { - numPhotonsVsPhiEta25Mom10Pion->Fill(partPhi, recoPhotons); - } - if (partMom > 3.5 && partMom < 4.5) - { - numPhotonsVsPhiEta25Mom4Pion->Fill(partPhi, recoPhotons); - } - } + recoPionPhoton[momSelect]->Fill(recoPhotons); + recoPhiVsEtaPionPhoton[momSelect]->Fill(partEta, partPhi, recoPhotons); + partPhiVsEtaPion[momSelect]->Fill(partEta, partPhi); + if (recoPDG == 211) + partPhiVsEtaPionMatch[momSelect]->Fill(partEta, partPhi); + if (recoPDG != 211) + partPhiVsEtaPionUnMatch[momSelect]->Fill(partEta, partPhi); } + if (partEta < -2.49 && partEta > -2.51) + { + if (partMom > 9.5 && partMom < 10.5) + { + numPhotonsVsPhiEta25Mom10Pion->Fill(partPhi, recoPhotons); + } + if (partMom > 3.5 && partMom < 4.5) + { + numPhotonsVsPhiEta25Mom4Pion->Fill(partPhi, recoPhotons); + } + } + } - if (partPDG == 321) // Look at Thrown Kaons + if (partPDG == 321) // Look at Thrown Kaons + { + if (recoPhotons < 4) + { + fkb[pbin][tbin][phibin]++; + } + else { - if (recoPhotons < 4) - { - fkb[pbin][tbin][phibin]++; - } - else - { - if (recoPDG == 11) - fke[pbin][tbin][phibin]++; - if (recoPDG == 211) - fkpi[pbin][tbin][phibin]++; - if (recoPDG == 321) - fkk[pbin][tbin][phibin]++; - if (recoPDG == 2212) - fkp[pbin][tbin][phibin]++; - } - if (momSelect != -1) - { - recoPhiVsEtaKaonPhoton[momSelect]->Fill(partEta, partPhi, recoPhotons); - partPhiVsEtaKaon[momSelect]->Fill(partEta, partPhi); - if (recoPDG == 321) - partPhiVsEtaKaonMatch[momSelect]->Fill(partEta, partPhi); - if (recoPDG != 321) - partPhiVsEtaKaonUnMatch[momSelect]->Fill(partEta, partPhi); - } + if (recoPDG == 11) + fke[pbin][tbin][phibin]++; + if (recoPDG == 211) + fkpi[pbin][tbin][phibin]++; + if (recoPDG == 321) + fkk[pbin][tbin][phibin]++; + if (recoPDG == 2212) + fkp[pbin][tbin][phibin]++; } + if (momSelect != -1 && partEta > -2.6 && partEta < -2.4 && partPhi > 2 && partPhi < 2.5) + { + recoKaonPhoton[momSelect]->Fill(recoPhotons); + recoPhiVsEtaKaonPhoton[momSelect]->Fill(partEta, partPhi, recoPhotons); + partPhiVsEtaKaon[momSelect]->Fill(partEta, partPhi); + if (recoPDG == 321) + partPhiVsEtaKaonMatch[momSelect]->Fill(partEta, partPhi); + if (recoPDG != 321) + partPhiVsEtaKaonUnMatch[momSelect]->Fill(partEta, partPhi); + } + } - if (partPDG == 2212) // Look at Thrown Protons + if (partPDG == 2212) // Look at Thrown Protons + { + if (recoPhotons < 4) + { + fpb[pbin][tbin][phibin]++; + } + else + { + if (recoPDG == 11) + fpe[pbin][tbin][phibin]++; + if (recoPDG == 211) + fppi[pbin][tbin][phibin]++; + if (recoPDG == 321) + fpk[pbin][tbin][phibin]++; + if (recoPDG == 2212) + fpp[pbin][tbin][phibin]++; + } + if (momSelect != -1 && partEta > -2.6 && partEta < -2.4 && partPhi > 2. && partPhi < 2.5) + { + recoProtonPhoton[momSelect]->Fill(recoPhotons); + recoPhiVsEtaProtonPhoton[momSelect]->Fill(partEta, partPhi, recoPhotons); + partPhiVsEtaProton[momSelect]->Fill(partEta, partPhi); + if (recoPDG == 2212) + partPhiVsEtaProtonMatch[momSelect]->Fill(partEta, partPhi); + if (recoPDG != 2212) + partPhiVsEtaProtonUnMatch[momSelect]->Fill(partEta, partPhi); + } + } + } + + TH2D *h[15][10][12]; + for (int pbin = 0; pbin < 15; pbin++) + { + for (int tbin = 0; tbin < 10; tbin++) + { + for (int phibin = 0; phibin < 12; phibin++) + { + h[pbin][tbin][phibin] = new TH2D(Form("p%it%iphi%i", pbin, tbin, phibin), Form("p%it%iphi%i;reco;truth", pbin, tbin, phibin), 5, 0., 5., 4, 0., 4.); + } + } + } + + for (int pbin = 0; pbin < 15; pbin++) + { + for (int tbin = 0; tbin < 10; tbin++) + { + for (int phibin = 0; phibin < 12; phibin++) { - if (recoPhotons < 4) - { - fpb[pbin][tbin][phibin]++; - } - else - { - if (recoPDG == 11) - fpe[pbin][tbin][phibin]++; - if (recoPDG == 211) - fppi[pbin][tbin][phibin]++; - if (recoPDG == 321) - fpk[pbin][tbin][phibin]++; - if (recoPDG == 2212) - fpp[pbin][tbin][phibin]++; - } - if (momSelect != -1) - { - recoPhiVsEtaProtonPhoton[momSelect]->Fill(partEta, partPhi, recoPhotons); - partPhiVsEtaProton[momSelect]->Fill(partEta, partPhi); - if (recoPDG == 2212) - partPhiVsEtaProtonMatch[momSelect]->Fill(partEta, partPhi); - if (recoPDG != 2212) - partPhiVsEtaProtonUnMatch[momSelect]->Fill(partEta, partPhi); - } + + Double_t fe_tot = fee[pbin][tbin][phibin] + fepi[pbin][tbin][phibin] + fek[pbin][tbin][phibin] + fep[pbin][tbin][phibin] + feb[pbin][tbin][phibin] + zero; + fee[pbin][tbin][phibin] /= fe_tot; + fepi[pbin][tbin][phibin] /= fe_tot; + fek[pbin][tbin][phibin] /= fe_tot; + fep[pbin][tbin][phibin] /= fe_tot; + feb[pbin][tbin][phibin] /= fe_tot; + + Double_t fpi_tot = fpie[pbin][tbin][phibin] + fpipi[pbin][tbin][phibin] + fpik[pbin][tbin][phibin] + fpip[pbin][tbin][phibin] + fpib[pbin][tbin][phibin] + zero; + fpie[pbin][tbin][phibin] /= fpi_tot; + fpipi[pbin][tbin][phibin] /= fpi_tot; + fpik[pbin][tbin][phibin] /= fpi_tot; + fpip[pbin][tbin][phibin] /= fpi_tot; + fpib[pbin][tbin][phibin] /= fpi_tot; + + Double_t fk_tot = fke[pbin][tbin][phibin] + fkpi[pbin][tbin][phibin] + fkk[pbin][tbin][phibin] + fkp[pbin][tbin][phibin] + fkb[pbin][tbin][phibin] + zero; + fke[pbin][tbin][phibin] /= fk_tot; + fkpi[pbin][tbin][phibin] /= fk_tot; + fkk[pbin][tbin][phibin] /= fk_tot; + fkp[pbin][tbin][phibin] /= fk_tot; + fkb[pbin][tbin][phibin] /= fk_tot; + + Double_t fp_tot = fpe[pbin][tbin][phibin] + fppi[pbin][tbin][phibin] + fpk[pbin][tbin][phibin] + fpp[pbin][tbin][phibin] + fpb[pbin][tbin][phibin] + zero; + fpe[pbin][tbin][phibin] /= fp_tot; + fppi[pbin][tbin][phibin] /= fp_tot; + fpk[pbin][tbin][phibin] /= fp_tot; + fpp[pbin][tbin][phibin] /= fp_tot; + fpb[pbin][tbin][phibin] /= fp_tot; + + h[pbin][tbin][phibin]->SetBinContent(1, 4, fee[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(2, 4, fepi[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(3, 4, fek[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(4, 4, fep[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(5, 4, feb[pbin][tbin][phibin]); + + h[pbin][tbin][phibin]->SetBinContent(1, 3, fpie[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(2, 3, fpipi[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(3, 3, fpik[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(4, 3, fpip[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(5, 3, fpib[pbin][tbin][phibin]); + + h[pbin][tbin][phibin]->SetBinContent(1, 2, fke[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(2, 2, fkpi[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(3, 2, fkk[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(4, 2, fkp[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(5, 2, fkb[pbin][tbin][phibin]); + + h[pbin][tbin][phibin]->SetBinContent(1, 1, fpe[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(2, 1, fppi[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(3, 1, fpk[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(4, 1, fpp[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(5, 1, fpb[pbin][tbin][phibin]); + + Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1); + Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1); + Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); + + if (pbin == 0 && tbin == 0 && phibin == 0) + { + cout << "truth PID, charge, p, theta, phi, prob(e), prob(pi), prob(k), prob(p), prob(fail)" << endl; + cout << "11 -1 " << p << " " << t << " " << phi << " " << fee[pbin][tbin][phibin] << " " << fepi[pbin][tbin][phibin] << " " << fek[pbin][tbin][phibin] << " " << fep[pbin][tbin][phibin] << " " << feb[pbin][tbin][phibin] << endl; + cout << "211 1 " << p << " " << t << " " << phi << " " << fpie[pbin][tbin][phibin] << " " << fpipi[pbin][tbin][phibin] << " " << fpik[pbin][tbin][phibin] << " " << fpip[pbin][tbin][phibin] << " " << fpib[pbin][tbin][phibin] << endl; + cout << "321 1 " << p << " " << t << " " << phi << " " << fke[pbin][tbin][phibin] << " " << fkpi[pbin][tbin][phibin] << " " << fkk[pbin][tbin][phibin] << " " << fkp[pbin][tbin][phibin] << " " << fkb[pbin][tbin][phibin] << endl; + cout << "2212 1 " << p << " " << t << " " << phi << " " << fpe[pbin][tbin][phibin] << " " << fppi[pbin][tbin][phibin] << " " << fpk[pbin][tbin][phibin] << " " << fpp[pbin][tbin][phibin] << " " << fpb[pbin][tbin][phibin] << endl; + } } - } + } + } - TH2D *h[15][16][360]; - for (int pbin = 0; pbin < 15; pbin++) - { - for (int tbin = 0; tbin < 16; tbin++) + FILE *filePointer; + filePointer = fopen(output_txt, "w"); + + for (int pbin = 0; pbin < 15; pbin++) + { + for (int tbin = 0; tbin < 10; tbin++) + { + for (int phibin = 0; phibin < 12; phibin++) { - for (int phibin = 0; phibin < 360; phibin++) - { - h[pbin][tbin][phibin] = new TH2D(Form("p%it%iphi%i", pbin, tbin, phibin), Form("p%it%iphi%i;reco;truth", pbin, tbin, phibin), 5, 0., 5., 4, 0., 4.); - } + Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1); + Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1); + Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); + + fprintf(filePointer, "11 -1 %f %f %f %f %f %f %f %f\n", p, t, phi, fee[pbin][tbin][phibin], fepi[pbin][tbin][phibin], fek[pbin][tbin][phibin], fep[pbin][tbin][phibin], feb[pbin][tbin][phibin]); } - } + } + } + + for (int pbin = 0; pbin < 15; pbin++) + { + for (int tbin = 0; tbin < 10; tbin++) + { + for (int phibin = 0; phibin < 12; phibin++) + { + Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1); + Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1); + Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); - FILE *filePointer; - filePointer = fopen(output_txt, "w"); + fprintf(filePointer, "211 1 %f %f %f %f %f %f %f %f\n", p, t, phi, fpie[pbin][tbin][phibin], fpipi[pbin][tbin][phibin], fpik[pbin][tbin][phibin], fpip[pbin][tbin][phibin], fpib[pbin][tbin][phibin]); + } + } + } + + for (int pbin = 0; pbin < 15; pbin++) + { + for (int tbin = 0; tbin < 10; tbin++) + { + for (int phibin = 0; phibin < 12; phibin++) + { + Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1); + Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1); + Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); - for (int pbin = 0; pbin < 15; pbin++) - { - for (int tbin = 0; tbin < 16; tbin++) + fprintf(filePointer, "321 1 %f %f %f %f %f %f %f %f\n", p, t, phi, fke[pbin][tbin][phibin], fkpi[pbin][tbin][phibin], fkk[pbin][tbin][phibin], fkp[pbin][tbin][phibin], fkb[pbin][tbin][phibin]); + } + } + } + + for (int pbin = 0; pbin < 15; pbin++) + { + for (int tbin = 0; tbin < 10; tbin++) + { + for (int phibin = 0; phibin < 12; phibin++) { - for (int phibin = 0; phibin < 360; phibin++) - { - - Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1); - Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1); - Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); - - Double_t mc_e_rec_e = (1.0 * fee[pbin][tbin][phibin]) / (1.0 * (fee[pbin][tbin][phibin] + fepi[pbin][tbin][phibin] + fek[pbin][tbin][phibin] + fep[pbin][tbin][phibin] + feb[pbin][tbin][phibin] + zero)); - Double_t mc_e_rec_pi = (1.0 * fepi[pbin][tbin][phibin]) / (1.0 * (fee[pbin][tbin][phibin] + fepi[pbin][tbin][phibin] + fek[pbin][tbin][phibin] + fep[pbin][tbin][phibin] + feb[pbin][tbin][phibin] + zero)); - Double_t mc_e_rec_k = (1.0 * fek[pbin][tbin][phibin]) / (1.0 * (fee[pbin][tbin][phibin] + fepi[pbin][tbin][phibin] + fek[pbin][tbin][phibin] + fep[pbin][tbin][phibin] + feb[pbin][tbin][phibin] + zero)); - Double_t mc_e_rec_p = (1.0 * fep[pbin][tbin][phibin]) / (1.0 * (fee[pbin][tbin][phibin] + fepi[pbin][tbin][phibin] + fek[pbin][tbin][phibin] + fep[pbin][tbin][phibin] + feb[pbin][tbin][phibin] + zero)); - Double_t mc_e_rec_b = (1.0 * feb[pbin][tbin][phibin]) / (1.0 * (fee[pbin][tbin][phibin] + fepi[pbin][tbin][phibin] + fek[pbin][tbin][phibin] + fep[pbin][tbin][phibin] + feb[pbin][tbin][phibin] + zero)); - - Double_t mc_pi_rec_e = (1.0 * fpie[pbin][tbin][phibin]) / (1.0 * (fpie[pbin][tbin][phibin] + fpipi[pbin][tbin][phibin] + fpik[pbin][tbin][phibin] + fpip[pbin][tbin][phibin] + fpib[pbin][tbin][phibin] + zero)); - Double_t mc_pi_rec_pi = (1.0 * fpipi[pbin][tbin][phibin]) / (1.0 * (fpie[pbin][tbin][phibin] + fpipi[pbin][tbin][phibin] + fpik[pbin][tbin][phibin] + fpip[pbin][tbin][phibin] + fpib[pbin][tbin][phibin] + zero)); - Double_t mc_pi_rec_k = (1.0 * fpik[pbin][tbin][phibin]) / (1.0 * (fpie[pbin][tbin][phibin] + fpipi[pbin][tbin][phibin] + fpik[pbin][tbin][phibin] + fpip[pbin][tbin][phibin] + fpib[pbin][tbin][phibin] + zero)); - Double_t mc_pi_rec_p = (1.0 * fpip[pbin][tbin][phibin]) / (1.0 * (fpie[pbin][tbin][phibin] + fpipi[pbin][tbin][phibin] + fpik[pbin][tbin][phibin] + fpip[pbin][tbin][phibin] + fpib[pbin][tbin][phibin] + zero)); - Double_t mc_pi_rec_b = (1.0 * fpib[pbin][tbin][phibin]) / (1.0 * (fpie[pbin][tbin][phibin] + fpipi[pbin][tbin][phibin] + fpik[pbin][tbin][phibin] + fpip[pbin][tbin][phibin] + fpib[pbin][tbin][phibin] + zero)); - - Double_t mc_k_rec_e = (1.0 * fke[pbin][tbin][phibin]) / (1.0 * (fke[pbin][tbin][phibin] + fkpi[pbin][tbin][phibin] + fkk[pbin][tbin][phibin] + fkp[pbin][tbin][phibin] + fkb[pbin][tbin][phibin] + zero)); - Double_t mc_k_rec_pi = (1.0 * fkpi[pbin][tbin][phibin]) / (1.0 * (fke[pbin][tbin][phibin] + fkpi[pbin][tbin][phibin] + fkk[pbin][tbin][phibin] + fkp[pbin][tbin][phibin] + fkb[pbin][tbin][phibin] + zero)); - Double_t mc_k_rec_k = (1.0 * fkk[pbin][tbin][phibin]) / (1.0 * (fke[pbin][tbin][phibin] + fkpi[pbin][tbin][phibin] + fkk[pbin][tbin][phibin] + fkp[pbin][tbin][phibin] + fkb[pbin][tbin][phibin] + zero)); - Double_t mc_k_rec_p = (1.0 * fkp[pbin][tbin][phibin]) / (1.0 * (fke[pbin][tbin][phibin] + fkpi[pbin][tbin][phibin] + fkk[pbin][tbin][phibin] + fkp[pbin][tbin][phibin] + fkb[pbin][tbin][phibin] + zero)); - Double_t mc_k_rec_b = (1.0 * fkb[pbin][tbin][phibin]) / (1.0 * (fke[pbin][tbin][phibin] + fkpi[pbin][tbin][phibin] + fkk[pbin][tbin][phibin] + fkp[pbin][tbin][phibin] + fkb[pbin][tbin][phibin] + zero)); - - Double_t mc_p_rec_e = (1.0 * fpe[pbin][tbin][phibin]) / (1.0 * (fpe[pbin][tbin][phibin] + fppi[pbin][tbin][phibin] + fpk[pbin][tbin][phibin] + fpp[pbin][tbin][phibin] + fpb[pbin][tbin][phibin] + zero)); - Double_t mc_p_rec_pi = (1.0 * fppi[pbin][tbin][phibin]) / (1.0 * (fpe[pbin][tbin][phibin] + fppi[pbin][tbin][phibin] + fpk[pbin][tbin][phibin] + fpp[pbin][tbin][phibin] + fpb[pbin][tbin][phibin] + zero)); - Double_t mc_p_rec_k = (1.0 * fpk[pbin][tbin][phibin]) / (1.0 * (fpe[pbin][tbin][phibin] + fppi[pbin][tbin][phibin] + fpk[pbin][tbin][phibin] + fpp[pbin][tbin][phibin] + fpb[pbin][tbin][phibin] + zero)); - Double_t mc_p_rec_p = (1.0 * fpp[pbin][tbin][phibin]) / (1.0 * (fpe[pbin][tbin][phibin] + fppi[pbin][tbin][phibin] + fpk[pbin][tbin][phibin] + fpp[pbin][tbin][phibin] + fpb[pbin][tbin][phibin] + zero)); - Double_t mc_p_rec_b = (1.0 * fpb[pbin][tbin][phibin]) / (1.0 * (fpe[pbin][tbin][phibin] + fppi[pbin][tbin][phibin] + fpk[pbin][tbin][phibin] + fpp[pbin][tbin][phibin] + fpb[pbin][tbin][phibin] + zero)); - - h[pbin][tbin][phibin]->SetBinContent(1, 4, mc_e_rec_e); - h[pbin][tbin][phibin]->SetBinContent(2, 4, mc_e_rec_pi); - h[pbin][tbin][phibin]->SetBinContent(3, 4, mc_e_rec_k); - h[pbin][tbin][phibin]->SetBinContent(4, 4, mc_e_rec_p); - h[pbin][tbin][phibin]->SetBinContent(5, 4, mc_e_rec_b); - - h[pbin][tbin][phibin]->SetBinContent(1, 3, mc_pi_rec_e); - h[pbin][tbin][phibin]->SetBinContent(2, 3, mc_pi_rec_pi); - h[pbin][tbin][phibin]->SetBinContent(3, 3, mc_pi_rec_k); - h[pbin][tbin][phibin]->SetBinContent(4, 3, mc_pi_rec_p); - h[pbin][tbin][phibin]->SetBinContent(5, 3, mc_pi_rec_b); - - h[pbin][tbin][phibin]->SetBinContent(1, 2, mc_k_rec_e); - h[pbin][tbin][phibin]->SetBinContent(2, 2, mc_k_rec_pi); - h[pbin][tbin][phibin]->SetBinContent(3, 2, mc_k_rec_k); - h[pbin][tbin][phibin]->SetBinContent(4, 2, mc_k_rec_p); - h[pbin][tbin][phibin]->SetBinContent(5, 2, mc_k_rec_b); - - h[pbin][tbin][phibin]->SetBinContent(1, 1, mc_p_rec_e); - h[pbin][tbin][phibin]->SetBinContent(2, 1, mc_p_rec_pi); - h[pbin][tbin][phibin]->SetBinContent(3, 1, mc_p_rec_k); - h[pbin][tbin][phibin]->SetBinContent(4, 1, mc_p_rec_p); - h[pbin][tbin][phibin]->SetBinContent(5, 1, mc_p_rec_b); - - fprintf(filePointer, "11 %f %f %f %f %f %f %f %f\n", p, t, phi, mc_e_rec_e, mc_e_rec_pi, mc_e_rec_k, mc_e_rec_p, mc_e_rec_b); - fprintf(filePointer, "211 %f %f %f %f %f %f %f %f\n", p, t, phi, mc_pi_rec_e, mc_pi_rec_pi, mc_pi_rec_k, mc_pi_rec_p, mc_pi_rec_b); - fprintf(filePointer, "321 %f %f %f %f %f %f %f %f\n", p, t, phi, mc_k_rec_e, mc_k_rec_pi, mc_k_rec_k, mc_k_rec_p, mc_k_rec_b); - fprintf(filePointer, "2212 %f %f %f %f %f %f %f %f\n", p, t, phi, mc_p_rec_e, mc_p_rec_pi, mc_p_rec_k, mc_p_rec_p, mc_p_rec_b); - - if (pbin == 0 && tbin == 0 && phibin == 0) - { - cout << "truth PID, p, theta, phi, prob(e), prob(pi), prob(k), prob(p), prob(fail)" << endl; - cout << "11 " << p << " " << t << " " << phi << " " << mc_e_rec_e << " " << mc_e_rec_pi << " " << mc_e_rec_k << " " << mc_e_rec_p << " " << mc_e_rec_b << endl; - cout << "211 " << p << " " << t << " " << phi << " " << mc_pi_rec_e << " " << mc_pi_rec_pi << " " << mc_pi_rec_k << " " << mc_pi_rec_p << " " << mc_pi_rec_b << endl; - cout << "321 " << p << " " << t << " " << phi << " " << mc_k_rec_e << " " << mc_k_rec_pi << " " << mc_k_rec_k << " " << mc_k_rec_p << " " << mc_k_rec_b << endl; - cout << "2212 " << p << " " << t << " " << phi << " " << mc_p_rec_e << " " << mc_p_rec_pi << " " << mc_p_rec_k << " " << mc_p_rec_p << " " << mc_p_rec_b << endl; - } - if (pbin == 14 && tbin == 9 && phibin == 11) - { - cout << "11 " << p << " " << t << " " << phi << " " << mc_e_rec_e << " " << mc_e_rec_pi << " " << mc_e_rec_k << " " << mc_e_rec_p << " " << mc_e_rec_b << endl; - cout << "211 " << p << " " << t << " " << phi << " " << mc_pi_rec_e << " " << mc_pi_rec_pi << " " << mc_pi_rec_k << " " << mc_pi_rec_p << " " << mc_pi_rec_b << endl; - cout << "321 " << p << " " << t << " " << phi << " " << mc_k_rec_e << " " << mc_k_rec_pi << " " << mc_k_rec_k << " " << mc_k_rec_p << " " << mc_k_rec_b << endl; - cout << "2212 " << p << " " << t << " " << phi << " " << mc_p_rec_e << " " << mc_p_rec_pi << " " << mc_p_rec_k << " " << mc_p_rec_p << " " << mc_p_rec_b << endl; - } - } + Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1); + Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1); + Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); + + fprintf(filePointer, "2212 1 %f %f %f %f %f %f %f %f\n", p, t, phi, fpe[pbin][tbin][phibin], fppi[pbin][tbin][phibin], fpk[pbin][tbin][phibin], fpp[pbin][tbin][phibin], fpb[pbin][tbin][phibin]); } - } + } + } - fclose(filePointer); - ofile->Write(); - ofile->Close(); + fclose(filePointer); + ofile->Write(); + ofile->Close(); } From 4abeb1801d0f0ed61b8bd5ac47027d6616db7a8c Mon Sep 17 00:00:00 2001 From: bspage912 Date: Sun, 5 May 2024 00:26:04 -0400 Subject: [PATCH 06/10] Changed generated particle theta range from 2.70 < theta < 3.10 to 2.65 < theta < 3.10 (eta -1.5 to -1.4) in runSimu.sh --- runSimu.sh | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/runSimu.sh b/runSimu.sh index 1a18562..fbe17f7 100755 --- a/runSimu.sh +++ b/runSimu.sh @@ -18,10 +18,10 @@ echo ${fileNameKaon} echo ${fileNameProton} ## Supply absolute path to the hepmc_writer.C script -root -b -l -q "/yourDir/hepmc_writer.C+(\"${fileNameElectron}\",0,2.70,3.10,0.1,15.0,${NUM})" -root -b -l -q "/yourDir/hepmc_writer.C+(\"${fileNamePion}\",1,2.70,3.10,0.1,15.0,${NUM})" -root -b -l -q "/yourDir/hepmc_writer.C+(\"${fileNameKaon}\",2,2.70,3.10,0.1,15.0,${NUM})" -root -b -l -q "/yourDir/hepmc_writer.C+(\"${fileNameProton}\",3,2.70,3.10,0.1,15.0,${NUM})" +root -b -l -q "/yourDir/hepmc_writer.C+(\"${fileNameElectron}\",0,2.65,3.10,0.1,15.0,${NUM})" +root -b -l -q "/yourDir/hepmc_writer.C+(\"${fileNamePion}\",1,2.65,3.10,0.1,15.0,${NUM})" +root -b -l -q "/yourDir/hepmc_writer.C+(\"${fileNameKaon}\",2,2.65,3.10,0.1,15.0,${NUM})" +root -b -l -q "/yourDir/hepmc_writer.C+(\"${fileNameProton}\",3,2.65,3.10,0.1,15.0,${NUM})" simuNameElectron=pfRICH-epic_electron_${RUN}.root simuNamePion=pfRICH-epic_pion_${RUN}.root From bfcb664a7727738da0194bc9c4945cc5c1360dee Mon Sep 17 00:00:00 2001 From: Youqi Song Date: Mon, 6 May 2024 16:03:41 -0400 Subject: [PATCH 07/10] Updated binnings --- readTree.C | 525 +++++++++++++++++++++++++++-------------------------- 1 file changed, 267 insertions(+), 258 deletions(-) diff --git a/readTree.C b/readTree.C index 8f9c221..3115933 100644 --- a/readTree.C +++ b/readTree.C @@ -1,294 +1,303 @@ void readTree(const char *input, const char *output, const char *output_txt) { - cout << "running readTree.C" << endl; - - TFile *fa = new TFile(input); - - TTree *ta = (TTree *)fa->Get("T"); - - TFile *ofile = TFile::Open(output, "recreate"); - - Int_t partPDG, recoPDG, recoPhotons, recoCherenkovHitCount; - Double_t partMom, partEta, partTheta, partPhi, partVertX, partVertY, partVertZ, recoTrackCherenkovAngle; - Double_t zero = pow(10., -50.); - - ta->SetBranchAddress("partPDG", &partPDG); - ta->SetBranchAddress("partMom", &partMom); - ta->SetBranchAddress("partEta", &partEta); - ta->SetBranchAddress("partTheta", &partTheta); - ta->SetBranchAddress("partPhi", &partPhi); - ta->SetBranchAddress("partVertX", &partVertX); - ta->SetBranchAddress("partVertY", &partVertY); - ta->SetBranchAddress("partVertZ", &partVertZ); - ta->SetBranchAddress("recoPDG", &recoPDG); - ta->SetBranchAddress("recoPhotons", &recoPhotons); - ta->SetBranchAddress("recoTrackCherenkovAngle", &recoTrackCherenkovAngle); - - TH1D *h_p = new TH1D("p", "", 15, 0.1, 15.1); - TH1D *h_t = new TH1D("theta", "", 20, 2.7, 3.1); - TH1D *h_phi = new TH1D("phi", "", 360, -1.0 * TMath::Pi(), TMath::Pi()); - - static Float_t fee[15][20][360], fepi[15][20][360], fek[15][20][360], fep[15][20][360], feb[15][20][360]; - static Float_t fpie[15][20][360], fpipi[15][20][360], fpik[15][20][360], fpip[15][20][360], fpib[15][20][360]; - static Float_t fke[15][20][360], fkpi[15][20][360], fkk[15][20][360], fkp[15][20][360], fkb[15][20][360]; - static Float_t fpe[15][20][360], fppi[15][20][360], fpk[15][20][360], fpp[15][20][360], fpb[15][20][360]; - - for (int pbin = 0; pbin < 15; pbin++) - { - for (int tbin = 0; tbin < 20; tbin++) + cout << "running readTree.C" << endl; + + TFile *fa = new TFile(input); + + TTree *ta = (TTree *)fa->Get("T"); + + Int_t partPDG, recoPDG, recoPhotons, recoCherenkovHitCount; + Double_t partMom, partEta, partTheta, partPhi, partVertX, partVertY, partVertZ, recoTrackCherenkovAngle; + Double_t zero = pow(10., -50.); + + ta->SetBranchAddress("partPDG", &partPDG); + ta->SetBranchAddress("partMom", &partMom); + ta->SetBranchAddress("partEta", &partEta); + ta->SetBranchAddress("partTheta", &partTheta); + ta->SetBranchAddress("partPhi", &partPhi); + ta->SetBranchAddress("partVertX", &partVertX); + ta->SetBranchAddress("partVertY", &partVertY); + ta->SetBranchAddress("partVertZ", &partVertZ); + ta->SetBranchAddress("recoPDG", &recoPDG); + ta->SetBranchAddress("recoPhotons", &recoPhotons); + ta->SetBranchAddress("recoTrackCherenkovAngle", &recoTrackCherenkovAngle); + + const Int_t N_PBINS = 20; + Double_t pedges[N_PBINS + 1] = {0.4, 0.8, 1.2, 1.6, 2.0, 2.4, 2.8, 3.2, 3.6, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0}; + //TH1D *h_p = new TH1D("p", "", N_PBINS, 0.1, 15.0); + TH1D *h_p = new TH1D("p", "", N_PBINS, pedges); + TH1D *h_t = new TH1D("theta", "", 20, 2.65, 3.1); + TH1D *h_phi = new TH1D("phi", "", 120, -TMath::Pi(), TMath::Pi()); + + static Float_t fee[20][20][120], fepi[20][20][120], fek[20][20][120], fep[20][20][120], feb[20][20][120]; + static Float_t fpie[20][20][120], fpipi[20][20][120], fpik[20][20][120], fpip[20][20][120], fpib[20][20][120]; + static Float_t fke[20][20][120], fkpi[20][20][120], fkk[20][20][120], fkp[20][20][120], fkb[20][20][120]; + static Float_t fpe[20][20][120], fppi[20][20][120], fpk[20][20][120], fpp[20][20][120], fpb[20][20][120]; + + for (int pbin = 0; pbin < 20; pbin++) + { + for (int tbin = 0; tbin < 20; tbin++) + { + for (int phibin = 0; phibin < 120; phibin++) { - for (int phibin = 0; phibin < 360; phibin++) - { - fee[pbin][tbin][phibin] = fepi[pbin][tbin][phibin] = fek[pbin][tbin][phibin] = fep[pbin][tbin][phibin] = feb[pbin][tbin][phibin] = 0.; - fpie[pbin][tbin][phibin] = fpipi[pbin][tbin][phibin] = fpik[pbin][tbin][phibin] = fpip[pbin][tbin][phibin] = fpib[pbin][tbin][phibin] = 0.; - fke[pbin][tbin][phibin] = fkpi[pbin][tbin][phibin] = fkk[pbin][tbin][phibin] = fkp[pbin][tbin][phibin] = fkb[pbin][tbin][phibin] = 0.; - fpe[pbin][tbin][phibin] = fppi[pbin][tbin][phibin] = fpk[pbin][tbin][phibin] = fpp[pbin][tbin][phibin] = fpb[pbin][tbin][phibin] = 0.; - } + fee[pbin][tbin][phibin] = fepi[pbin][tbin][phibin] = fek[pbin][tbin][phibin] = fep[pbin][tbin][phibin] = feb[pbin][tbin][phibin] = 0.; + fpie[pbin][tbin][phibin] = fpipi[pbin][tbin][phibin] = fpik[pbin][tbin][phibin] = fpip[pbin][tbin][phibin] = fpib[pbin][tbin][phibin] = 0.; + fke[pbin][tbin][phibin] = fkpi[pbin][tbin][phibin] = fkk[pbin][tbin][phibin] = fkp[pbin][tbin][phibin] = fkb[pbin][tbin][phibin] = 0.; + fpe[pbin][tbin][phibin] = fppi[pbin][tbin][phibin] = fpk[pbin][tbin][phibin] = fpp[pbin][tbin][phibin] = fpb[pbin][tbin][phibin] = 0.; } - } + } + } - Int_t nEntries = (Int_t)ta->GetEntries(); - for (int i = 0; i < nEntries; i++) - { - if (i % 100000 == 0) - cout << "Processing Event " << i << endl; - ta->GetEntry(i); + Int_t nEntries = (Int_t)ta->GetEntries(); + for (int i = 0; i < nEntries; i++) + { + if (i % 1000000 == 0) + cout << "Processing Event " << i << endl; + ta->GetEntry(i); - Int_t pbin = h_p->FindBin(partMom) - 1; - Int_t tbin = h_t->FindBin(partTheta) - 1; - Int_t phibin = h_phi->FindBin(partPhi) - 1; + Int_t pbin = h_p->FindBin(partMom) - 1; + Int_t tbin = h_t->FindBin(partTheta) - 1; + Int_t phibin = h_phi->FindBin(partPhi) - 1; - // if (partMom < 1) cout << pbin << endl;// -> bin number starts from 1 + // if (partMom < 1) cout << pbin << endl;// -> bin number starts from 1 - if (partPDG == 11) // Look at Thrown Electrons + if (partPDG == 11) // Look at Thrown Electrons + { + if (recoPhotons < 3) { - if (recoPhotons < 4) - { - feb[pbin][tbin][phibin]++; - } - else - { - if (recoPDG == 11) - fee[pbin][tbin][phibin]++; - if (recoPDG == 211) - fepi[pbin][tbin][phibin]++; - if (recoPDG == 321) - fek[pbin][tbin][phibin]++; - if (recoPDG == 2212) - fep[pbin][tbin][phibin]++; - } + feb[pbin][tbin][phibin]++; } - - if (partPDG == 211) // Look at Thrown Pions + else { - if (recoPhotons < 4) - { - fpib[pbin][tbin][phibin]++; - } - else - { - if (recoPDG == 11) - fpie[pbin][tbin][phibin]++; - if (recoPDG == 211) - fpipi[pbin][tbin][phibin]++; - if (recoPDG == 321) - fpik[pbin][tbin][phibin]++; - if (recoPDG == 2212) - fpip[pbin][tbin][phibin]++; - } + if (recoPDG == 11) + fee[pbin][tbin][phibin]++; + if (recoPDG == 211) + fepi[pbin][tbin][phibin]++; + if (recoPDG == 321) + fek[pbin][tbin][phibin]++; + if (recoPDG == 2212) + fep[pbin][tbin][phibin]++; } + } - if (partPDG == 321) // Look at Thrown Kaons + if (partPDG == 211) // Look at Thrown Pions + { + if (recoPhotons < 3) { - if (recoPhotons < 4) - { - fkb[pbin][tbin][phibin]++; - } - else - { - if (recoPDG == 11) - fke[pbin][tbin][phibin]++; - if (recoPDG == 211) - fkpi[pbin][tbin][phibin]++; - if (recoPDG == 321) - fkk[pbin][tbin][phibin]++; - if (recoPDG == 2212) - fkp[pbin][tbin][phibin]++; - } + fpib[pbin][tbin][phibin]++; } - - if (partPDG == 2212) // Look at Thrown Protons + else { - if (recoPhotons < 4) - { - fpb[pbin][tbin][phibin]++; - } - else - { - if (recoPDG == 11) - fpe[pbin][tbin][phibin]++; - if (recoPDG == 211) - fppi[pbin][tbin][phibin]++; - if (recoPDG == 321) - fpk[pbin][tbin][phibin]++; - if (recoPDG == 2212) - fpp[pbin][tbin][phibin]++; - } + if (recoPDG == 11) + fpie[pbin][tbin][phibin]++; + if (recoPDG == 211) + fpipi[pbin][tbin][phibin]++; + if (recoPDG == 321) + fpik[pbin][tbin][phibin]++; + if (recoPDG == 2212) + fpip[pbin][tbin][phibin]++; } - } + } - TH2D *h[15][20][360]; - for (int pbin = 0; pbin < 15; pbin++) - { - for (int tbin = 0; tbin < 20; tbin++) + if (partPDG == 321) // Look at Thrown Kaons + { + if (recoPhotons < 3) { - for (int phibin = 0; phibin < 360; phibin++) - { - h[pbin][tbin][phibin] = new TH2D(Form("p%it%iphi%i", pbin, tbin, phibin), Form("p%it%iphi%i;reco;truth", pbin, tbin, phibin), 5, 0., 5., 4, 0., 4.); - } + fkb[pbin][tbin][phibin]++; } - } - - for (int pbin = 0; pbin < 15; pbin++) - { - for (int tbin = 0; tbin < 20; tbin++) + else { - for (int phibin = 0; phibin < 360; phibin++) - { - - Double_t fe_tot = fee[pbin][tbin][phibin] + fepi[pbin][tbin][phibin] + fek[pbin][tbin][phibin] + fep[pbin][tbin][phibin] + feb[pbin][tbin][phibin] + zero; - fee[pbin][tbin][phibin] /= fe_tot; - fepi[pbin][tbin][phibin] /= fe_tot; - fek[pbin][tbin][phibin] /= fe_tot; - fep[pbin][tbin][phibin] /= fe_tot; - feb[pbin][tbin][phibin] /= fe_tot; - - Double_t fpi_tot = fpie[pbin][tbin][phibin] + fpipi[pbin][tbin][phibin] + fpik[pbin][tbin][phibin] + fpip[pbin][tbin][phibin] + fpib[pbin][tbin][phibin] + zero; - fpie[pbin][tbin][phibin] /= fpi_tot; - fpipi[pbin][tbin][phibin] /= fpi_tot; - fpik[pbin][tbin][phibin] /= fpi_tot; - fpip[pbin][tbin][phibin] /= fpi_tot; - fpib[pbin][tbin][phibin] /= fpi_tot; - - Double_t fk_tot = fke[pbin][tbin][phibin] + fkpi[pbin][tbin][phibin] + fkk[pbin][tbin][phibin] + fkp[pbin][tbin][phibin] + fkb[pbin][tbin][phibin] + zero; - fke[pbin][tbin][phibin] /= fk_tot; - fkpi[pbin][tbin][phibin] /= fk_tot; - fkk[pbin][tbin][phibin] /= fk_tot; - fkp[pbin][tbin][phibin] /= fk_tot; - fkb[pbin][tbin][phibin] /= fk_tot; - - Double_t fp_tot = fpe[pbin][tbin][phibin] + fppi[pbin][tbin][phibin] + fpk[pbin][tbin][phibin] + fpp[pbin][tbin][phibin] + fpb[pbin][tbin][phibin] + zero; - fpe[pbin][tbin][phibin] /= fp_tot; - fppi[pbin][tbin][phibin] /= fp_tot; - fpk[pbin][tbin][phibin] /= fp_tot; - fpp[pbin][tbin][phibin] /= fp_tot; - fpb[pbin][tbin][phibin] /= fp_tot; - - h[pbin][tbin][phibin]->SetBinContent(1, 4, fee[pbin][tbin][phibin]); - h[pbin][tbin][phibin]->SetBinContent(2, 4, fepi[pbin][tbin][phibin]); - h[pbin][tbin][phibin]->SetBinContent(3, 4, fek[pbin][tbin][phibin]); - h[pbin][tbin][phibin]->SetBinContent(4, 4, fep[pbin][tbin][phibin]); - h[pbin][tbin][phibin]->SetBinContent(5, 4, feb[pbin][tbin][phibin]); - - h[pbin][tbin][phibin]->SetBinContent(1, 3, fpie[pbin][tbin][phibin]); - h[pbin][tbin][phibin]->SetBinContent(2, 3, fpipi[pbin][tbin][phibin]); - h[pbin][tbin][phibin]->SetBinContent(3, 3, fpik[pbin][tbin][phibin]); - h[pbin][tbin][phibin]->SetBinContent(4, 3, fpip[pbin][tbin][phibin]); - h[pbin][tbin][phibin]->SetBinContent(5, 3, fpib[pbin][tbin][phibin]); - - h[pbin][tbin][phibin]->SetBinContent(1, 2, fke[pbin][tbin][phibin]); - h[pbin][tbin][phibin]->SetBinContent(2, 2, fkpi[pbin][tbin][phibin]); - h[pbin][tbin][phibin]->SetBinContent(3, 2, fkk[pbin][tbin][phibin]); - h[pbin][tbin][phibin]->SetBinContent(4, 2, fkp[pbin][tbin][phibin]); - h[pbin][tbin][phibin]->SetBinContent(5, 2, fkb[pbin][tbin][phibin]); - - h[pbin][tbin][phibin]->SetBinContent(1, 1, fpe[pbin][tbin][phibin]); - h[pbin][tbin][phibin]->SetBinContent(2, 1, fppi[pbin][tbin][phibin]); - h[pbin][tbin][phibin]->SetBinContent(3, 1, fpk[pbin][tbin][phibin]); - h[pbin][tbin][phibin]->SetBinContent(4, 1, fpp[pbin][tbin][phibin]); - h[pbin][tbin][phibin]->SetBinContent(5, 1, fpb[pbin][tbin][phibin]); - - Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1); - Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1); - Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); - - if (pbin == 0 && tbin == 0 && phibin == 0) - { - cout << "truth PID, charge, p, theta, phi, prob(e), prob(pi), prob(k), prob(p), prob(fail)" << endl; - cout << "11 -1 " << p << " " << t << " " << phi << " " << fee[pbin][tbin][phibin] << " " << fepi[pbin][tbin][phibin] << " " << fek[pbin][tbin][phibin] << " " << fep[pbin][tbin][phibin] << " " << feb[pbin][tbin][phibin] << endl; - cout << "211 1 " << p << " " << t << " " << phi << " " << fpie[pbin][tbin][phibin] << " " << fpipi[pbin][tbin][phibin] << " " << fpik[pbin][tbin][phibin] << " " << fpip[pbin][tbin][phibin] << " " << fpib[pbin][tbin][phibin] << endl; - cout << "321 1 " << p << " " << t << " " << phi << " " << fke[pbin][tbin][phibin] << " " << fkpi[pbin][tbin][phibin] << " " << fkk[pbin][tbin][phibin] << " " << fkp[pbin][tbin][phibin] << " " << fkb[pbin][tbin][phibin] << endl; - cout << "2212 1 " << p << " " << t << " " << phi << " " << fpe[pbin][tbin][phibin] << " " << fppi[pbin][tbin][phibin] << " " << fpk[pbin][tbin][phibin] << " " << fpp[pbin][tbin][phibin] << " " << fpb[pbin][tbin][phibin] << endl; - } - } + if (recoPDG == 11) + fke[pbin][tbin][phibin]++; + if (recoPDG == 211) + fkpi[pbin][tbin][phibin]++; + if (recoPDG == 321) + fkk[pbin][tbin][phibin]++; + if (recoPDG == 2212) + fkp[pbin][tbin][phibin]++; } - } + } - FILE *filePointer; - filePointer = fopen(output_txt, "w"); - - for (int pbin = 0; pbin < 15; pbin++) - { - for (int tbin = 0; tbin < 20; tbin++) + if (partPDG == 2212) // Look at Thrown Protons + { + if (recoPhotons < 3) { - for (int phibin = 0; phibin < 360; phibin++) - { - Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1); - Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1); - Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); - - fprintf(filePointer, "11 -1 %f %f %f %f %f %f %f %f\n", p, t, phi, fee[pbin][tbin][phibin], fepi[pbin][tbin][phibin], fek[pbin][tbin][phibin], fep[pbin][tbin][phibin], feb[pbin][tbin][phibin]); - } + fpb[pbin][tbin][phibin]++; } - } - - for (int pbin = 0; pbin < 15; pbin++) - { - for (int tbin = 0; tbin < 20; tbin++) + else { - for (int phibin = 0; phibin < 360; phibin++) - { - Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1); - Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1); - Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); - - fprintf(filePointer, "211 1 %f %f %f %f %f %f %f %f\n", p, t, phi, fpie[pbin][tbin][phibin], fpipi[pbin][tbin][phibin], fpik[pbin][tbin][phibin], fpip[pbin][tbin][phibin], fpib[pbin][tbin][phibin]); - } + if (recoPDG == 11) + fpe[pbin][tbin][phibin]++; + if (recoPDG == 211) + fppi[pbin][tbin][phibin]++; + if (recoPDG == 321) + fpk[pbin][tbin][phibin]++; + if (recoPDG == 2212) + fpp[pbin][tbin][phibin]++; + } + } + } + + TFile *ofile = TFile::Open(output, "recreate"); + TH2D *h[20][20][120]; + for (int pbin = 0; pbin < 20; pbin++) + { + for (int tbin = 0; tbin < 20; tbin++) + { + for (int phibin = 0; phibin < 120; phibin++) + { + h[pbin][tbin][phibin] = new TH2D(Form("p%it%iphi%i", pbin, tbin, phibin), Form("p%it%iphi%i;reco;truth", pbin, tbin, phibin), 5, 0., 5., 4, 0., 4.); + } + } + } + + for (int pbin = 0; pbin < 20; pbin++) + { + for (int tbin = 0; tbin < 20; tbin++) + { + for (int phibin = 0; phibin < 120; phibin++) + { + + Double_t fe_tot = fee[pbin][tbin][phibin] + fepi[pbin][tbin][phibin] + fek[pbin][tbin][phibin] + fep[pbin][tbin][phibin] + feb[pbin][tbin][phibin] + zero; + fee[pbin][tbin][phibin] /= fe_tot; + fepi[pbin][tbin][phibin] /= fe_tot; + fek[pbin][tbin][phibin] /= fe_tot; + fep[pbin][tbin][phibin] /= fe_tot; + feb[pbin][tbin][phibin] /= fe_tot; + + Double_t fpi_tot = fpie[pbin][tbin][phibin] + fpipi[pbin][tbin][phibin] + fpik[pbin][tbin][phibin] + fpip[pbin][tbin][phibin] + fpib[pbin][tbin][phibin] + zero; + fpie[pbin][tbin][phibin] /= fpi_tot; + fpipi[pbin][tbin][phibin] /= fpi_tot; + fpik[pbin][tbin][phibin] /= fpi_tot; + fpip[pbin][tbin][phibin] /= fpi_tot; + fpib[pbin][tbin][phibin] /= fpi_tot; + + Double_t fk_tot = fke[pbin][tbin][phibin] + fkpi[pbin][tbin][phibin] + fkk[pbin][tbin][phibin] + fkp[pbin][tbin][phibin] + fkb[pbin][tbin][phibin] + zero; + fke[pbin][tbin][phibin] /= fk_tot; + fkpi[pbin][tbin][phibin] /= fk_tot; + fkk[pbin][tbin][phibin] /= fk_tot; + fkp[pbin][tbin][phibin] /= fk_tot; + fkb[pbin][tbin][phibin] /= fk_tot; + + Double_t fp_tot = fpe[pbin][tbin][phibin] + fppi[pbin][tbin][phibin] + fpk[pbin][tbin][phibin] + fpp[pbin][tbin][phibin] + fpb[pbin][tbin][phibin] + zero; + fpe[pbin][tbin][phibin] /= fp_tot; + fppi[pbin][tbin][phibin] /= fp_tot; + fpk[pbin][tbin][phibin] /= fp_tot; + fpp[pbin][tbin][phibin] /= fp_tot; + fpb[pbin][tbin][phibin] /= fp_tot; + + h[pbin][tbin][phibin]->SetBinContent(1, 4, fee[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(2, 4, fepi[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(3, 4, fek[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(4, 4, fep[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(5, 4, feb[pbin][tbin][phibin]); + + h[pbin][tbin][phibin]->SetBinContent(1, 3, fpie[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(2, 3, fpipi[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(3, 3, fpik[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(4, 3, fpip[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(5, 3, fpib[pbin][tbin][phibin]); + + h[pbin][tbin][phibin]->SetBinContent(1, 2, fke[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(2, 2, fkpi[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(3, 2, fkk[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(4, 2, fkp[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(5, 2, fkb[pbin][tbin][phibin]); + + h[pbin][tbin][phibin]->SetBinContent(1, 1, fpe[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(2, 1, fppi[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(3, 1, fpk[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(4, 1, fpp[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(5, 1, fpb[pbin][tbin][phibin]); + + Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1); + Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1); + Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); + + if (pbin == 0 && tbin == 0 && phibin == 0) + { + cout << "truth PID, charge, p, theta, phi, prob(e), prob(pi), prob(k), prob(p), prob(fail)" << endl; + cout << "11 -1 " << p << " " << t << " " << phi << " " << fee[pbin][tbin][phibin] << " " << fepi[pbin][tbin][phibin] << " " << fek[pbin][tbin][phibin] << " " << fep[pbin][tbin][phibin] << " " << feb[pbin][tbin][phibin] << endl; + cout << "211 1 " << p << " " << t << " " << phi << " " << fpie[pbin][tbin][phibin] << " " << fpipi[pbin][tbin][phibin] << " " << fpik[pbin][tbin][phibin] << " " << fpip[pbin][tbin][phibin] << " " << fpib[pbin][tbin][phibin] << endl; + cout << "321 1 " << p << " " << t << " " << phi << " " << fke[pbin][tbin][phibin] << " " << fkpi[pbin][tbin][phibin] << " " << fkk[pbin][tbin][phibin] << " " << fkp[pbin][tbin][phibin] << " " << fkb[pbin][tbin][phibin] << endl; + cout << "2212 1 " << p << " " << t << " " << phi << " " << fpe[pbin][tbin][phibin] << " " << fppi[pbin][tbin][phibin] << " " << fpk[pbin][tbin][phibin] << " " << fpp[pbin][tbin][phibin] << " " << fpb[pbin][tbin][phibin] << endl; + } } - } + } + } - for (int pbin = 0; pbin < 15; pbin++) - { - for (int tbin = 0; tbin < 20; tbin++) + FILE *filePointer; + filePointer = fopen(output_txt, "w"); + + for (int pbin = 0; pbin < 20; pbin++) + { + for (int tbin = 0; tbin < 20; tbin++) + { + for (int phibin = 0; phibin < 120; phibin++) { - for (int phibin = 0; phibin < 360; phibin++) - { - Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1); - Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1); - Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); - - fprintf(filePointer, "321 1 %f %f %f %f %f %f %f %f\n", p, t, phi, fke[pbin][tbin][phibin], fkpi[pbin][tbin][phibin], fkk[pbin][tbin][phibin], fkp[pbin][tbin][phibin], fkb[pbin][tbin][phibin]); - } + Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1); + Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1); + Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); + + fprintf(filePointer, "11 -1 %f %f %f %f %f %f %f %f\n", p, t, phi, fee[pbin][tbin][phibin], fepi[pbin][tbin][phibin], fek[pbin][tbin][phibin], fep[pbin][tbin][phibin], feb[pbin][tbin][phibin]); } - } + } + } + cout << "Filled out LUT for e" << endl; + + for (int pbin = 0; pbin < 20; pbin++) + { + for (int tbin = 0; tbin < 20; tbin++) + { + for (int phibin = 0; phibin < 120; phibin++) + { + Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1); + Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1); + Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); - for (int pbin = 0; pbin < 15; pbin++) - { - for (int tbin = 0; tbin < 20; tbin++) + fprintf(filePointer, "211 1 %f %f %f %f %f %f %f %f\n", p, t, phi, fpie[pbin][tbin][phibin], fpipi[pbin][tbin][phibin], fpik[pbin][tbin][phibin], fpip[pbin][tbin][phibin], fpib[pbin][tbin][phibin]); + } + } + } + cout << "Filled out LUT for pions" << endl; + + for (int pbin = 0; pbin < 20; pbin++) + { + for (int tbin = 0; tbin < 20; tbin++) + { + for (int phibin = 0; phibin < 120; phibin++) { - for (int phibin = 0; phibin < 360; phibin++) - { - Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1); - Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1); - Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); - - fprintf(filePointer, "2212 1 %f %f %f %f %f %f %f %f\n", p, t, phi, fpe[pbin][tbin][phibin], fppi[pbin][tbin][phibin], fpk[pbin][tbin][phibin], fpp[pbin][tbin][phibin], fpb[pbin][tbin][phibin]); - } + Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1); + Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1); + Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); + + fprintf(filePointer, "321 1 %f %f %f %f %f %f %f %f\n", p, t, phi, fke[pbin][tbin][phibin], fkpi[pbin][tbin][phibin], fkk[pbin][tbin][phibin], fkp[pbin][tbin][phibin], fkb[pbin][tbin][phibin]); } - } + } + } + cout << "Filled out LUT for kaons" << endl; + + for (int pbin = 0; pbin < 20; pbin++) + { + for (int tbin = 0; tbin < 20; tbin++) + { + for (int phibin = 0; phibin < 120; phibin++) + { + Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1); + Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1); + Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); + + fprintf(filePointer, "2212 1 %f %f %f %f %f %f %f %f\n", p, t, phi, fpe[pbin][tbin][phibin], fppi[pbin][tbin][phibin], fpk[pbin][tbin][phibin], fpp[pbin][tbin][phibin], fpb[pbin][tbin][phibin]); + } + } + } + cout << "Filled out LUT for protons" << endl; + + fclose(filePointer); - fclose(filePointer); - ofile->Write(); - ofile->Close(); + ofile->Write(); + cout << "Written to output file" << endl; + ofile->Close(); + cout << "Closed output file" << endl; } From f10a9ad6efad53699a86c786edbfb60634569de6 Mon Sep 17 00:00:00 2001 From: Youqi Song Date: Tue, 7 May 2024 01:16:07 -0400 Subject: [PATCH 08/10] Updated binnings. Take care of momentum underflow bins --- readTree.C | 542 ++++++++++++++++++++++++++--------------------------- 1 file changed, 270 insertions(+), 272 deletions(-) diff --git a/readTree.C b/readTree.C index 3115933..5e26b14 100644 --- a/readTree.C +++ b/readTree.C @@ -1,303 +1,301 @@ void readTree(const char *input, const char *output, const char *output_txt) { - cout << "running readTree.C" << endl; - - TFile *fa = new TFile(input); - - TTree *ta = (TTree *)fa->Get("T"); - - Int_t partPDG, recoPDG, recoPhotons, recoCherenkovHitCount; - Double_t partMom, partEta, partTheta, partPhi, partVertX, partVertY, partVertZ, recoTrackCherenkovAngle; - Double_t zero = pow(10., -50.); - - ta->SetBranchAddress("partPDG", &partPDG); - ta->SetBranchAddress("partMom", &partMom); - ta->SetBranchAddress("partEta", &partEta); - ta->SetBranchAddress("partTheta", &partTheta); - ta->SetBranchAddress("partPhi", &partPhi); - ta->SetBranchAddress("partVertX", &partVertX); - ta->SetBranchAddress("partVertY", &partVertY); - ta->SetBranchAddress("partVertZ", &partVertZ); - ta->SetBranchAddress("recoPDG", &recoPDG); - ta->SetBranchAddress("recoPhotons", &recoPhotons); - ta->SetBranchAddress("recoTrackCherenkovAngle", &recoTrackCherenkovAngle); - - const Int_t N_PBINS = 20; - Double_t pedges[N_PBINS + 1] = {0.4, 0.8, 1.2, 1.6, 2.0, 2.4, 2.8, 3.2, 3.6, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0}; - //TH1D *h_p = new TH1D("p", "", N_PBINS, 0.1, 15.0); - TH1D *h_p = new TH1D("p", "", N_PBINS, pedges); - TH1D *h_t = new TH1D("theta", "", 20, 2.65, 3.1); - TH1D *h_phi = new TH1D("phi", "", 120, -TMath::Pi(), TMath::Pi()); - - static Float_t fee[20][20][120], fepi[20][20][120], fek[20][20][120], fep[20][20][120], feb[20][20][120]; - static Float_t fpie[20][20][120], fpipi[20][20][120], fpik[20][20][120], fpip[20][20][120], fpib[20][20][120]; - static Float_t fke[20][20][120], fkpi[20][20][120], fkk[20][20][120], fkp[20][20][120], fkb[20][20][120]; - static Float_t fpe[20][20][120], fppi[20][20][120], fpk[20][20][120], fpp[20][20][120], fpb[20][20][120]; - - for (int pbin = 0; pbin < 20; pbin++) - { - for (int tbin = 0; tbin < 20; tbin++) - { - for (int phibin = 0; phibin < 120; phibin++) + cout << "running readTree.C" << endl; + + TFile *fa = new TFile(input); + + TTree *ta = (TTree *)fa->Get("T"); + + TFile *ofile = TFile::Open(output, "recreate"); + + Int_t partPDG, recoPDG, recoPhotons, recoCherenkovHitCount; + Double_t partMom, partEta, partTheta, partPhi, partVertX, partVertY, partVertZ, recoTrackCherenkovAngle; + Double_t zero = pow(10., -50.); + + ta->SetBranchAddress("partPDG", &partPDG); + ta->SetBranchAddress("partMom", &partMom); + ta->SetBranchAddress("partEta", &partEta); + ta->SetBranchAddress("partTheta", &partTheta); + ta->SetBranchAddress("partPhi", &partPhi); + ta->SetBranchAddress("partVertX", &partVertX); + ta->SetBranchAddress("partVertY", &partVertY); + ta->SetBranchAddress("partVertZ", &partVertZ); + ta->SetBranchAddress("recoPDG", &recoPDG); + ta->SetBranchAddress("recoPhotons", &recoPhotons); + ta->SetBranchAddress("recoTrackCherenkovAngle", &recoTrackCherenkovAngle); + + TH1D *h_p = new TH1D("p", "", 37, 0.4, 15.2); + TH1D *h_t = new TH1D("theta", "", 20, 2.65, 3.1); + TH1D *h_phi = new TH1D("phi", "", 120, -1.0 * TMath::Pi(), TMath::Pi()); + + static Float_t fee[37][20][120], fepi[37][20][120], fek[37][20][120], fep[37][20][120], feb[37][20][120]; + static Float_t fpie[37][20][120], fpipi[37][20][120], fpik[37][20][120], fpip[37][20][120], fpib[37][20][120]; + static Float_t fke[37][20][120], fkpi[37][20][120], fkk[37][20][120], fkp[37][20][120], fkb[37][20][120]; + static Float_t fpe[37][20][120], fppi[37][20][120], fpk[37][20][120], fpp[37][20][120], fpb[37][20][120]; + + for (int pbin = 0; pbin < 37; pbin++) + { + for (int tbin = 0; tbin < 20; tbin++) { - fee[pbin][tbin][phibin] = fepi[pbin][tbin][phibin] = fek[pbin][tbin][phibin] = fep[pbin][tbin][phibin] = feb[pbin][tbin][phibin] = 0.; - fpie[pbin][tbin][phibin] = fpipi[pbin][tbin][phibin] = fpik[pbin][tbin][phibin] = fpip[pbin][tbin][phibin] = fpib[pbin][tbin][phibin] = 0.; - fke[pbin][tbin][phibin] = fkpi[pbin][tbin][phibin] = fkk[pbin][tbin][phibin] = fkp[pbin][tbin][phibin] = fkb[pbin][tbin][phibin] = 0.; - fpe[pbin][tbin][phibin] = fppi[pbin][tbin][phibin] = fpk[pbin][tbin][phibin] = fpp[pbin][tbin][phibin] = fpb[pbin][tbin][phibin] = 0.; + for (int phibin = 0; phibin < 120; phibin++) + { + fee[pbin][tbin][phibin] = fepi[pbin][tbin][phibin] = fek[pbin][tbin][phibin] = fep[pbin][tbin][phibin] = feb[pbin][tbin][phibin] = 0.; + fpie[pbin][tbin][phibin] = fpipi[pbin][tbin][phibin] = fpik[pbin][tbin][phibin] = fpip[pbin][tbin][phibin] = fpib[pbin][tbin][phibin] = 0.; + fke[pbin][tbin][phibin] = fkpi[pbin][tbin][phibin] = fkk[pbin][tbin][phibin] = fkp[pbin][tbin][phibin] = fkb[pbin][tbin][phibin] = 0.; + fpe[pbin][tbin][phibin] = fppi[pbin][tbin][phibin] = fpk[pbin][tbin][phibin] = fpp[pbin][tbin][phibin] = fpb[pbin][tbin][phibin] = 0.; + } } - } - } - - Int_t nEntries = (Int_t)ta->GetEntries(); - for (int i = 0; i < nEntries; i++) - { - if (i % 1000000 == 0) - cout << "Processing Event " << i << endl; - ta->GetEntry(i); - - Int_t pbin = h_p->FindBin(partMom) - 1; - Int_t tbin = h_t->FindBin(partTheta) - 1; - Int_t phibin = h_phi->FindBin(partPhi) - 1; - - // if (partMom < 1) cout << pbin << endl;// -> bin number starts from 1 - - if (partPDG == 11) // Look at Thrown Electrons - { - if (recoPhotons < 3) - { - feb[pbin][tbin][phibin]++; - } - else + } + + Int_t nEntries = (Int_t)ta->GetEntries(); + for (int i = 0; i < nEntries; i++) + { + if (i % 100000 == 0) + cout << "Processing Event " << i << endl; + ta->GetEntry(i); + + Int_t pbin = h_p->FindBin(partMom) - 1; + Int_t tbin = h_t->FindBin(partTheta) - 1; + Int_t phibin = h_phi->FindBin(partPhi) - 1; + + // if (partMom < 1) cout << pbin << endl;// -> bin number starts from 1 + + if (pbin < 0) + continue; + if (tbin < 0) + continue; + if (phibin < 0) + continue + + if (partPDG == 11) // Look at Thrown Electrons { - if (recoPDG == 11) - fee[pbin][tbin][phibin]++; - if (recoPDG == 211) - fepi[pbin][tbin][phibin]++; - if (recoPDG == 321) - fek[pbin][tbin][phibin]++; - if (recoPDG == 2212) - fep[pbin][tbin][phibin]++; + if (recoPhotons < 3) + { + feb[pbin][tbin][phibin]++; + } + else + { + if (recoPDG == 11) + fee[pbin][tbin][phibin]++; + if (recoPDG == 211) + fepi[pbin][tbin][phibin]++; + if (recoPDG == 321) + fek[pbin][tbin][phibin]++; + if (recoPDG == 2212) + fep[pbin][tbin][phibin]++; + } } - } - if (partPDG == 211) // Look at Thrown Pions - { - if (recoPhotons < 3) + if (partPDG == 211) // Look at Thrown Pions { - fpib[pbin][tbin][phibin]++; + if (recoPhotons < 3) + { + fpib[pbin][tbin][phibin]++; + } + else + { + if (recoPDG == 11) + fpie[pbin][tbin][phibin]++; + if (recoPDG == 211) + fpipi[pbin][tbin][phibin]++; + if (recoPDG == 321) + fpik[pbin][tbin][phibin]++; + if (recoPDG == 2212) + fpip[pbin][tbin][phibin]++; + } } - else - { - if (recoPDG == 11) - fpie[pbin][tbin][phibin]++; - if (recoPDG == 211) - fpipi[pbin][tbin][phibin]++; - if (recoPDG == 321) - fpik[pbin][tbin][phibin]++; - if (recoPDG == 2212) - fpip[pbin][tbin][phibin]++; - } - } - if (partPDG == 321) // Look at Thrown Kaons - { - if (recoPhotons < 3) - { - fkb[pbin][tbin][phibin]++; - } - else + if (partPDG == 321) // Look at Thrown Kaons { - if (recoPDG == 11) - fke[pbin][tbin][phibin]++; - if (recoPDG == 211) - fkpi[pbin][tbin][phibin]++; - if (recoPDG == 321) - fkk[pbin][tbin][phibin]++; - if (recoPDG == 2212) - fkp[pbin][tbin][phibin]++; + if (recoPhotons < 3) + { + fkb[pbin][tbin][phibin]++; + } + else + { + if (recoPDG == 11) + fke[pbin][tbin][phibin]++; + if (recoPDG == 211) + fkpi[pbin][tbin][phibin]++; + if (recoPDG == 321) + fkk[pbin][tbin][phibin]++; + if (recoPDG == 2212) + fkp[pbin][tbin][phibin]++; + } } - } - if (partPDG == 2212) // Look at Thrown Protons - { - if (recoPhotons < 3) - { - fpb[pbin][tbin][phibin]++; - } - else + if (partPDG == 2212) // Look at Thrown Protons { - if (recoPDG == 11) - fpe[pbin][tbin][phibin]++; - if (recoPDG == 211) - fppi[pbin][tbin][phibin]++; - if (recoPDG == 321) - fpk[pbin][tbin][phibin]++; - if (recoPDG == 2212) - fpp[pbin][tbin][phibin]++; + if (recoPhotons < 3) + { + fpb[pbin][tbin][phibin]++; + } + else + { + if (recoPDG == 11) + fpe[pbin][tbin][phibin]++; + if (recoPDG == 211) + fppi[pbin][tbin][phibin]++; + if (recoPDG == 321) + fpk[pbin][tbin][phibin]++; + if (recoPDG == 2212) + fpp[pbin][tbin][phibin]++; + } } - } - } - - TFile *ofile = TFile::Open(output, "recreate"); - TH2D *h[20][20][120]; - for (int pbin = 0; pbin < 20; pbin++) - { - for (int tbin = 0; tbin < 20; tbin++) - { - for (int phibin = 0; phibin < 120; phibin++) + } + + TH2D *h[37][20][120]; + for (int pbin = 0; pbin < 37; pbin++) + { + for (int tbin = 0; tbin < 20; tbin++) { - h[pbin][tbin][phibin] = new TH2D(Form("p%it%iphi%i", pbin, tbin, phibin), Form("p%it%iphi%i;reco;truth", pbin, tbin, phibin), 5, 0., 5., 4, 0., 4.); + for (int phibin = 0; phibin < 120; phibin++) + { + h[pbin][tbin][phibin] = new TH2D(Form("p%it%iphi%i", pbin, tbin, phibin), Form("p%it%iphi%i;reco;truth", pbin, tbin, phibin), 5, 0., 5., 4, 0., 4.); + } } - } - } - - for (int pbin = 0; pbin < 20; pbin++) - { - for (int tbin = 0; tbin < 20; tbin++) - { - for (int phibin = 0; phibin < 120; phibin++) - { + } - Double_t fe_tot = fee[pbin][tbin][phibin] + fepi[pbin][tbin][phibin] + fek[pbin][tbin][phibin] + fep[pbin][tbin][phibin] + feb[pbin][tbin][phibin] + zero; - fee[pbin][tbin][phibin] /= fe_tot; - fepi[pbin][tbin][phibin] /= fe_tot; - fek[pbin][tbin][phibin] /= fe_tot; - fep[pbin][tbin][phibin] /= fe_tot; - feb[pbin][tbin][phibin] /= fe_tot; - - Double_t fpi_tot = fpie[pbin][tbin][phibin] + fpipi[pbin][tbin][phibin] + fpik[pbin][tbin][phibin] + fpip[pbin][tbin][phibin] + fpib[pbin][tbin][phibin] + zero; - fpie[pbin][tbin][phibin] /= fpi_tot; - fpipi[pbin][tbin][phibin] /= fpi_tot; - fpik[pbin][tbin][phibin] /= fpi_tot; - fpip[pbin][tbin][phibin] /= fpi_tot; - fpib[pbin][tbin][phibin] /= fpi_tot; - - Double_t fk_tot = fke[pbin][tbin][phibin] + fkpi[pbin][tbin][phibin] + fkk[pbin][tbin][phibin] + fkp[pbin][tbin][phibin] + fkb[pbin][tbin][phibin] + zero; - fke[pbin][tbin][phibin] /= fk_tot; - fkpi[pbin][tbin][phibin] /= fk_tot; - fkk[pbin][tbin][phibin] /= fk_tot; - fkp[pbin][tbin][phibin] /= fk_tot; - fkb[pbin][tbin][phibin] /= fk_tot; - - Double_t fp_tot = fpe[pbin][tbin][phibin] + fppi[pbin][tbin][phibin] + fpk[pbin][tbin][phibin] + fpp[pbin][tbin][phibin] + fpb[pbin][tbin][phibin] + zero; - fpe[pbin][tbin][phibin] /= fp_tot; - fppi[pbin][tbin][phibin] /= fp_tot; - fpk[pbin][tbin][phibin] /= fp_tot; - fpp[pbin][tbin][phibin] /= fp_tot; - fpb[pbin][tbin][phibin] /= fp_tot; - - h[pbin][tbin][phibin]->SetBinContent(1, 4, fee[pbin][tbin][phibin]); - h[pbin][tbin][phibin]->SetBinContent(2, 4, fepi[pbin][tbin][phibin]); - h[pbin][tbin][phibin]->SetBinContent(3, 4, fek[pbin][tbin][phibin]); - h[pbin][tbin][phibin]->SetBinContent(4, 4, fep[pbin][tbin][phibin]); - h[pbin][tbin][phibin]->SetBinContent(5, 4, feb[pbin][tbin][phibin]); - - h[pbin][tbin][phibin]->SetBinContent(1, 3, fpie[pbin][tbin][phibin]); - h[pbin][tbin][phibin]->SetBinContent(2, 3, fpipi[pbin][tbin][phibin]); - h[pbin][tbin][phibin]->SetBinContent(3, 3, fpik[pbin][tbin][phibin]); - h[pbin][tbin][phibin]->SetBinContent(4, 3, fpip[pbin][tbin][phibin]); - h[pbin][tbin][phibin]->SetBinContent(5, 3, fpib[pbin][tbin][phibin]); - - h[pbin][tbin][phibin]->SetBinContent(1, 2, fke[pbin][tbin][phibin]); - h[pbin][tbin][phibin]->SetBinContent(2, 2, fkpi[pbin][tbin][phibin]); - h[pbin][tbin][phibin]->SetBinContent(3, 2, fkk[pbin][tbin][phibin]); - h[pbin][tbin][phibin]->SetBinContent(4, 2, fkp[pbin][tbin][phibin]); - h[pbin][tbin][phibin]->SetBinContent(5, 2, fkb[pbin][tbin][phibin]); - - h[pbin][tbin][phibin]->SetBinContent(1, 1, fpe[pbin][tbin][phibin]); - h[pbin][tbin][phibin]->SetBinContent(2, 1, fppi[pbin][tbin][phibin]); - h[pbin][tbin][phibin]->SetBinContent(3, 1, fpk[pbin][tbin][phibin]); - h[pbin][tbin][phibin]->SetBinContent(4, 1, fpp[pbin][tbin][phibin]); - h[pbin][tbin][phibin]->SetBinContent(5, 1, fpb[pbin][tbin][phibin]); - - Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1); - Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1); - Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); - - if (pbin == 0 && tbin == 0 && phibin == 0) - { - cout << "truth PID, charge, p, theta, phi, prob(e), prob(pi), prob(k), prob(p), prob(fail)" << endl; - cout << "11 -1 " << p << " " << t << " " << phi << " " << fee[pbin][tbin][phibin] << " " << fepi[pbin][tbin][phibin] << " " << fek[pbin][tbin][phibin] << " " << fep[pbin][tbin][phibin] << " " << feb[pbin][tbin][phibin] << endl; - cout << "211 1 " << p << " " << t << " " << phi << " " << fpie[pbin][tbin][phibin] << " " << fpipi[pbin][tbin][phibin] << " " << fpik[pbin][tbin][phibin] << " " << fpip[pbin][tbin][phibin] << " " << fpib[pbin][tbin][phibin] << endl; - cout << "321 1 " << p << " " << t << " " << phi << " " << fke[pbin][tbin][phibin] << " " << fkpi[pbin][tbin][phibin] << " " << fkk[pbin][tbin][phibin] << " " << fkp[pbin][tbin][phibin] << " " << fkb[pbin][tbin][phibin] << endl; - cout << "2212 1 " << p << " " << t << " " << phi << " " << fpe[pbin][tbin][phibin] << " " << fppi[pbin][tbin][phibin] << " " << fpk[pbin][tbin][phibin] << " " << fpp[pbin][tbin][phibin] << " " << fpb[pbin][tbin][phibin] << endl; - } + for (int pbin = 0; pbin < 37; pbin++) + { + for (int tbin = 0; tbin < 20; tbin++) + { + for (int phibin = 0; phibin < 120; phibin++) + { + + Double_t fe_tot = fee[pbin][tbin][phibin] + fepi[pbin][tbin][phibin] + fek[pbin][tbin][phibin] + fep[pbin][tbin][phibin] + feb[pbin][tbin][phibin] + zero; + fee[pbin][tbin][phibin] /= fe_tot; + fepi[pbin][tbin][phibin] /= fe_tot; + fek[pbin][tbin][phibin] /= fe_tot; + fep[pbin][tbin][phibin] /= fe_tot; + feb[pbin][tbin][phibin] /= fe_tot; + + Double_t fpi_tot = fpie[pbin][tbin][phibin] + fpipi[pbin][tbin][phibin] + fpik[pbin][tbin][phibin] + fpip[pbin][tbin][phibin] + fpib[pbin][tbin][phibin] + zero; + fpie[pbin][tbin][phibin] /= fpi_tot; + fpipi[pbin][tbin][phibin] /= fpi_tot; + fpik[pbin][tbin][phibin] /= fpi_tot; + fpip[pbin][tbin][phibin] /= fpi_tot; + fpib[pbin][tbin][phibin] /= fpi_tot; + + Double_t fk_tot = fke[pbin][tbin][phibin] + fkpi[pbin][tbin][phibin] + fkk[pbin][tbin][phibin] + fkp[pbin][tbin][phibin] + fkb[pbin][tbin][phibin] + zero; + fke[pbin][tbin][phibin] /= fk_tot; + fkpi[pbin][tbin][phibin] /= fk_tot; + fkk[pbin][tbin][phibin] /= fk_tot; + fkp[pbin][tbin][phibin] /= fk_tot; + fkb[pbin][tbin][phibin] /= fk_tot; + + Double_t fp_tot = fpe[pbin][tbin][phibin] + fppi[pbin][tbin][phibin] + fpk[pbin][tbin][phibin] + fpp[pbin][tbin][phibin] + fpb[pbin][tbin][phibin] + zero; + fpe[pbin][tbin][phibin] /= fp_tot; + fppi[pbin][tbin][phibin] /= fp_tot; + fpk[pbin][tbin][phibin] /= fp_tot; + fpp[pbin][tbin][phibin] /= fp_tot; + fpb[pbin][tbin][phibin] /= fp_tot; + + h[pbin][tbin][phibin]->SetBinContent(1, 4, fee[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(2, 4, fepi[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(3, 4, fek[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(4, 4, fep[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(5, 4, feb[pbin][tbin][phibin]); + + h[pbin][tbin][phibin]->SetBinContent(1, 3, fpie[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(2, 3, fpipi[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(3, 3, fpik[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(4, 3, fpip[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(5, 3, fpib[pbin][tbin][phibin]); + + h[pbin][tbin][phibin]->SetBinContent(1, 2, fke[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(2, 2, fkpi[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(3, 2, fkk[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(4, 2, fkp[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(5, 2, fkb[pbin][tbin][phibin]); + + h[pbin][tbin][phibin]->SetBinContent(1, 1, fpe[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(2, 1, fppi[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(3, 1, fpk[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(4, 1, fpp[pbin][tbin][phibin]); + h[pbin][tbin][phibin]->SetBinContent(5, 1, fpb[pbin][tbin][phibin]); + + Double_t p = h_p->GetXaxis()->GetBinCenter(pbin + 1); + Double_t t = h_t->GetXaxis()->GetBinCenter(tbin + 1); + Double_t phi = h_phi->GetXaxis()->GetBinCenter(phibin + 1); + + if (pbin == 0 && tbin == 0 && phibin == 0) + { + cout << "truth PID, charge, p, theta, phi, prob(e), prob(pi), prob(k), prob(p), prob(fail)" << endl; + cout << "11 -1 " << p << " " << t << " " << phi << " " << fee[pbin][tbin][phibin] << " " << fepi[pbin][tbin][phibin] << " " << fek[pbin][tbin][phibin] << " " << fep[pbin][tbin][phibin] << " " << feb[pbin][tbin][phibin] << endl; + cout << "211 1 " << p << " " << t << " " << phi << " " << fpie[pbin][tbin][phibin] << " " << fpipi[pbin][tbin][phibin] << " " << fpik[pbin][tbin][phibin] << " " << fpip[pbin][tbin][phibin] << " " << fpib[pbin][tbin][phibin] << endl; + cout << "321 1 " << p << " " << t << " " << phi << " " << fke[pbin][tbin][phibin] << " " << fkpi[pbin][tbin][phibin] << " " << fkk[pbin][tbin][phibin] << " " << fkp[pbin][tbin][phibin] << " " << fkb[pbin][tbin][phibin] << endl; + cout << "2212 1 " << p << " " << t << " " << phi << " " << fpe[pbin][tbin][phibin] << " " << fppi[pbin][tbin][phibin] << " " << fpk[pbin][tbin][phibin] << " " << fpp[pbin][tbin][phibin] << " " << fpb[pbin][tbin][phibin] << endl; + } + } } - } - } + } - FILE *filePointer; - filePointer = fopen(output_txt, "w"); + FILE *filePointer; + filePointer = fopen(output_txt, "w"); - for (int pbin = 0; pbin < 20; pbin++) - { - for (int tbin = 0; tbin < 20; tbin++) - { - for (int phibin = 0; phibin < 120; phibin++) + for (int pbin = 0; pbin < 37; pbin++) + { + for (int tbin = 0; tbin < 20; tbin++) { - Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1); - Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1); - Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); - - fprintf(filePointer, "11 -1 %f %f %f %f %f %f %f %f\n", p, t, phi, fee[pbin][tbin][phibin], fepi[pbin][tbin][phibin], fek[pbin][tbin][phibin], fep[pbin][tbin][phibin], feb[pbin][tbin][phibin]); + for (int phibin = 0; phibin < 120; phibin++) + { + Double_t p = h_p->GetXaxis()->GetBinCenter(pbin + 1); + Double_t t = h_t->GetXaxis()->GetBinCenter(tbin + 1); + Double_t phi = h_phi->GetXaxis()->GetBinCenter(phibin + 1); + + fprintf(filePointer, "11 -1 %f %f %f %f %f %f %f %f\n", p, t, phi, fee[pbin][tbin][phibin], fepi[pbin][tbin][phibin], fek[pbin][tbin][phibin], fep[pbin][tbin][phibin], feb[pbin][tbin][phibin]); + } } - } - } - cout << "Filled out LUT for e" << endl; - - for (int pbin = 0; pbin < 20; pbin++) - { - for (int tbin = 0; tbin < 20; tbin++) - { - for (int phibin = 0; phibin < 120; phibin++) - { - Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1); - Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1); - Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); + } - fprintf(filePointer, "211 1 %f %f %f %f %f %f %f %f\n", p, t, phi, fpie[pbin][tbin][phibin], fpipi[pbin][tbin][phibin], fpik[pbin][tbin][phibin], fpip[pbin][tbin][phibin], fpib[pbin][tbin][phibin]); - } - } - } - cout << "Filled out LUT for pions" << endl; - - for (int pbin = 0; pbin < 20; pbin++) - { - for (int tbin = 0; tbin < 20; tbin++) - { - for (int phibin = 0; phibin < 120; phibin++) + for (int pbin = 0; pbin < 37; pbin++) + { + for (int tbin = 0; tbin < 20; tbin++) { - Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1); - Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1); - Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); - - fprintf(filePointer, "321 1 %f %f %f %f %f %f %f %f\n", p, t, phi, fke[pbin][tbin][phibin], fkpi[pbin][tbin][phibin], fkk[pbin][tbin][phibin], fkp[pbin][tbin][phibin], fkb[pbin][tbin][phibin]); + for (int phibin = 0; phibin < 120; phibin++) + { + Double_t p = h_p->GetXaxis()->GetBinCenter(pbin + 1); + Double_t t = h_t->GetXaxis()->GetBinCenter(tbin + 1); + Double_t phi = h_phi->GetXaxis()->GetBinCenter(phibin + 1); + + fprintf(filePointer, "211 1 %f %f %f %f %f %f %f %f\n", p, t, phi, fpie[pbin][tbin][phibin], fpipi[pbin][tbin][phibin], fpik[pbin][tbin][phibin], fpip[pbin][tbin][phibin], fpib[pbin][tbin][phibin]); + } } - } - } - cout << "Filled out LUT for kaons" << endl; - - for (int pbin = 0; pbin < 20; pbin++) - { - for (int tbin = 0; tbin < 20; tbin++) - { - for (int phibin = 0; phibin < 120; phibin++) - { - Double_t p = h_p->GetXaxis()->GetBinLowEdge(pbin + 1); - Double_t t = h_t->GetXaxis()->GetBinLowEdge(tbin + 1); - Double_t phi = h_phi->GetXaxis()->GetBinLowEdge(phibin + 1); + } - fprintf(filePointer, "2212 1 %f %f %f %f %f %f %f %f\n", p, t, phi, fpe[pbin][tbin][phibin], fppi[pbin][tbin][phibin], fpk[pbin][tbin][phibin], fpp[pbin][tbin][phibin], fpb[pbin][tbin][phibin]); + for (int pbin = 0; pbin < 37; pbin++) + { + for (int tbin = 0; tbin < 20; tbin++) + { + for (int phibin = 0; phibin < 120; phibin++) + { + Double_t p = h_p->GetXaxis()->GetBinCenter(pbin + 1); + Double_t t = h_t->GetXaxis()->GetBinCenter(tbin + 1); + Double_t phi = h_phi->GetXaxis()->GetBinCenter(phibin + 1); + + fprintf(filePointer, "321 1 %f %f %f %f %f %f %f %f\n", p, t, phi, fke[pbin][tbin][phibin], fkpi[pbin][tbin][phibin], fkk[pbin][tbin][phibin], fkp[pbin][tbin][phibin], fkb[pbin][tbin][phibin]); + } } - } - } - cout << "Filled out LUT for protons" << endl; + } - fclose(filePointer); + for (int pbin = 0; pbin < 37; pbin++) + { + for (int tbin = 0; tbin < 20; tbin++) + { + for (int phibin = 0; phibin < 120; phibin++) + { + Double_t p = h_p->GetXaxis()->GetBinCenter(pbin + 1); + Double_t t = h_t->GetXaxis()->GetBinCenter(tbin + 1); + Double_t phi = h_phi->GetXaxis()->GetBinCenter(phibin + 1); + + fprintf(filePointer, "2212 1 %f %f %f %f %f %f %f %f\n", p, t, phi, fpe[pbin][tbin][phibin], fppi[pbin][tbin][phibin], fpk[pbin][tbin][phibin], fpp[pbin][tbin][phibin], fpb[pbin][tbin][phibin]); + } + } + } - ofile->Write(); - cout << "Written to output file" << endl; - ofile->Close(); - cout << "Closed output file" << endl; -} + fclose(filePointer); + ofile->Write(); + ofile->Close(); +} \ No newline at end of file From 1ba29b3e2f8e7a324a871e7005fa61fa0dfeb8b7 Mon Sep 17 00:00:00 2001 From: Youqi Song Date: Tue, 7 May 2024 02:13:07 -0400 Subject: [PATCH 09/10] Added a missing semicolon --- readTree.C | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/readTree.C b/readTree.C index 5e26b14..cc377a0 100644 --- a/readTree.C +++ b/readTree.C @@ -51,7 +51,7 @@ void readTree(const char *input, const char *output, const char *output_txt) Int_t nEntries = (Int_t)ta->GetEntries(); for (int i = 0; i < nEntries; i++) { - if (i % 100000 == 0) + if (i % 1000000 == 0) cout << "Processing Event " << i << endl; ta->GetEntry(i); @@ -63,10 +63,10 @@ void readTree(const char *input, const char *output, const char *output_txt) if (pbin < 0) continue; - if (tbin < 0) - continue; - if (phibin < 0) - continue + //if (tbin < 0) + // continue; + //if (phibin < 0) + // continue; if (partPDG == 11) // Look at Thrown Electrons { @@ -298,4 +298,4 @@ void readTree(const char *input, const char *output, const char *output_txt) fclose(filePointer); ofile->Write(); ofile->Close(); -} \ No newline at end of file +} From c62cdae32d8cf67cb1ab368b12e11d316cf7db6b Mon Sep 17 00:00:00 2001 From: Youqi Song Date: Tue, 7 May 2024 02:13:34 -0400 Subject: [PATCH 10/10] Added a missing semicolon --- readTree.C | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/readTree.C b/readTree.C index cc377a0..dede7a6 100644 --- a/readTree.C +++ b/readTree.C @@ -63,10 +63,10 @@ void readTree(const char *input, const char *output, const char *output_txt) if (pbin < 0) continue; - //if (tbin < 0) - // continue; - //if (phibin < 0) - // continue; + if (tbin < 0) + continue; + if (phibin < 0) + continue; if (partPDG == 11) // Look at Thrown Electrons {