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/readTree.C b/readTree.C new file mode 100644 index 0000000..dede7a6 --- /dev/null +++ b/readTree.C @@ -0,0 +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"); + + 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++) + { + 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 (pbin < 0) + continue; + if (tbin < 0) + continue; + if (phibin < 0) + continue; + + if (partPDG == 11) // Look at Thrown Electrons + { + 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) + { + 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 < 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 (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[37][20][120]; + for (int pbin = 0; pbin < 37; 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 < 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"); + + 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, "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 < 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, "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 < 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]); + } + } + } + + 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]); + } + } + } + + fclose(filePointer); + ofile->Write(); + ofile->Close(); +} diff --git a/readTreeQA.C b/readTreeQA.C new file mode 100644 index 0000000..b34ca6e --- /dev/null +++ b/readTreeQA.C @@ -0,0 +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", "", 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++) + { + 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 && 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 (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) + { + 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 (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 && 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 (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++) + { + + 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"); + + 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); + + 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); + + 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); + + 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++) + { + 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(); +} diff --git a/runSimu.sh b/runSimu.sh new file mode 100755 index 0000000..fbe17f7 --- /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.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 +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..e1f2a0d --- /dev/null +++ b/scripts/reco-epic-LUT.C @@ -0,0 +1,109 @@ +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 partTheta; + Double_t partPhi; + Double_t partVertX; + Double_t partVertY; + Double_t partVertZ; + + Int_t recoPDG; + Int_t recoPhotons; + Double_t recoTrackCherenkovAngle; + + 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("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); + 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; + { + CherenkovEvent *event; + + // Event loop; + 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 + } + + tree->Write(); + ofile->Close(); +}