far_forward_dvcs/analysis/analyze_DVCS_eicrecon.cxx

0001 //-------------------------
0002 //
0003 // Simple analysis code to analyze EPIC simulation output
0004 //
0005 //
0006 // Author: Alex Jentsch
0007 //
0008 //
0009 //
0010 //------------------------
0011
0012
0013 #include <iostream>
0014 #include <fstream>
0015 #include <TH1D.h>
0016 #include <TH2D.h>
0017 #include <TFile.h>
0018 #include <TTreeReader.h>
0019 #include <TTreeReaderArray.h>
0020 #include <TString.h>
0021 #include <TVector3.h>
0022
0023 using namespace std;
0024
0025 #include "detectorResolution.hpp"
0026 #include "render.hpp"
0027
0028 void analyze_DVCS_eicrecon(TString fileList = "inputFileList.list", TString outputDir = "."){
0029
0030     cout << "Input FileList: " << fileList << endl;
0031     string fileName;
0032     TFile * inputRootFile;
0033     TTree * rootTree;
0034     TString outputFileName = outputDir + "/output.root";
0035     cout << "Output file: " << outputFileName << endl;
0036
0037
0038     ifstream fileListStream;
0039     fileListStream.open(fileList);
0040     if(!fileListStream) { cout << "NO_LIST_FILE " << fileList << endl; return;}
0041
0042     //--------------------------------------------------------------------------
0043
0044     //histograms -- only a few for now
0045
0046     //MC information
0047     TH1D* h_eta_MC = new TH1D("h_eta",";Pseudorapidity, #eta",100,-10.0,10.0);
0048     TH1D* h_px_MC = new TH1D("px_MC", ";p_{x} [GeV/c]", 100, -10.0, 10.0);
0049     TH1D* h_py_MC = new TH1D("py_MC", ";p_{y} [GeV/c]", 100, -10.0, 10.0);
0050     TH1D* h_pt_MC = new TH1D("pt_MC", ";p_{t} [GeV/c]", 150, 0.0, 2.0);
0051     TH1D* h_pz_MC = new TH1D("pz_MC", ";p_{z} [GeV/c]", 100, 0.0, 320.0);
0052     TH1D* h_theta_MC = new TH1D("theta_MC", ";#theta [mrad]", 100, 0.0, 25.0);
0053     TH1D* h_phi_MC = new TH1D("phi_MC", ";#phi [rad]", 100, -3.2, 3.2);
0054
0055     TH1D* h_pt_squared_MC = new TH1D("pt_squared_MC", "; Momentum Transfer, -t [GeV^{2}]", 100, 0.0, 2.0);
0056
0057     //TRUTH information
0058     TH1D* h_eta_TRUTH = new TH1D("h_eta_TRUTH",";Pseudorapidity, #eta",100,-10.0,10.0);
0059     TH1D* h_px_TRUTH = new TH1D("px_TRUTH", ";p_{x} [GeV/c]", 100, -10.0, 10.0);
0060     TH1D* h_py_TRUTH = new TH1D("py_TRUTH", ";p_{y} [GeV/c]", 100, -10.0, 10.0);
0061     TH1D* h_pt_TRUTH = new TH1D("pt_TRUTH", ";p_{t} [GeV/c]", 150, 0.0, 2.0);
0062     TH1D* h_pz_TRUTH = new TH1D("pz_TRUTH", ";p_{z} [GeV/c]", 100, 0.0, 320.0);
0063     TH1D* h_theta_TRUTH = new TH1D("theta_TRUTH", ";#theta [mrad]", 100, 0.0, 25.0);
0064     TH1D* h_phi_TRUTH = new TH1D("phi_TRUTH", ";#phi [rad]", 100, -3.2, 3.2);
0065
0066     TH1D* h_pt_squared_TRUTH = new TH1D("pt_squared_TRUTH", "; Momentum Transfer, -t [GeV^{2}]", 100, 0.0, 2.0);
0067
0068     //ACCEPTANCE ONLY information
0069     TH1D* h_eta_accep = new TH1D("h_eta_accep",";Pseudorapidity, #eta",100,-10.0,10.0);
0070     TH1D* h_px_accep = new TH1D("px_accep", ";p_{x} [GeV/c]", 100, -10.0, 10.0);
0071     TH1D* h_py_accep = new TH1D("py_accep", ";p_{y} [GeV/c]", 100, -10.0, 10.0);
0072     TH1D* h_pt_accep = new TH1D("pt_accep", ";p_{t} [GeV/c]", 150, 0.0, 2.0);
0073     TH1D* h_pz_accep = new TH1D("pz_accep", ";p_{z} [GeV/c]", 100, 0.0, 320.0);
0074     TH1D* h_theta_accep = new TH1D("theta_accep", ";#theta [mrad]", 100, 0.0, 25.0);
0075     TH1D* h_phi_accep = new TH1D("phi_accep", ";#phi [rad]", 100, -3.2, 3.2);
0076
0077     TH1D* h_pt_squared_accep = new TH1D("pt_squared_accep", "; Momentum Transfer, -t [GeV^{2}]", 100, 0.0, 2.0);
0078
0079     //Roman pots
0080     TH1D* h_px_RomanPots = new TH1D("px_RomanPots", ";p_{x} [GeV/c]", 100, -10.0, 10.0);
0081     TH1D* h_py_RomanPots = new TH1D("py_RomanPots", ";p_{y} [GeV/c]", 100, -10.0, 10.0);
0082     TH1D* h_pt_RomanPots = new TH1D("pt_RomanPots", ";p_{t} [GeV/c]", 150, 0.0, 2.0);
0083     TH1D* h_pz_RomanPots = new TH1D("pz_RomanPots", ";p_{z} [GeV/c]", 100, 0.0, 320.0);
0084     TH1D* h_pt_squared_RomanPots = new TH1D("pt_squared_RomanPots", "; Momentum Transfer, -t [GeV^{2}]", 100, 0.0, 2.0);
0085     TH2D* h_rp_GLOBAL_occupancy_map = new TH2D("Roman_pots_GLOBAL_occupancy_map", ";hit x [mm];hit y [mm]", 100, -1300, -900, 100, -70, 70);
0086     TH2D* h_rp_LOCAL_occupancy_map = new TH2D("Roman_pots_LOCAL_occupancy_map", ";hit x [mm];hit y [mm]", 100, -150, 150, 100, -80, 80);
0087     //100, -150, 150, 100, -70, 70);
0088
0089     //OMD
0090     TH1D* h_px_OMD = new TH1D("px_OMD", ";p_{x} [GeV/c]", 100, -10.0, 10.0);
0091     TH1D* h_py_OMD = new TH1D("py_OMD", ";p_{y} [GeV/c]", 100, -10.0, 10.0);
0092     TH1D* h_pt_OMD = new TH1D("pt_OMD", ";p_{t} [GeV/c]", 100, 0.0, 2.0);
0093     TH1D* h_pz_OMD = new TH1D("pz_OMD", ";p_{z} [GeV/c]", 100, 0.0, 320.0);
0094     TH2D* h_omd_occupancy_map = new TH2D("OMD_occupancy_map", ";hit x [mm];hit y [mm]", 100, -150, 150, 100, -70, 70);
0095
0096
0097     //B0 tracker hits
0098     TH2D* h_B0_occupancy_map_layer_0 = new TH2D("B0_occupancy_map_0", "B0_occupancy_map_0", 100, -400, 0, 100, -170, 170);
0099     TH2D* h_B0_occupancy_map_layer_1 = new TH2D("B0_occupancy_map_1", "B0_occupancy_map_1", 100, -400, 0, 100, -170, 170);
0100     TH2D* h_B0_occupancy_map_layer_2 = new TH2D("B0_occupancy_map_2", "B0_occupancy_map_2", 100, -400, 0, 100, -170, 170);
0101     TH2D* h_B0_occupancy_map_layer_3 = new TH2D("B0_occupancy_map_3", "B0_occupancy_map_3", 100, -400, 0, 100, -170, 170);
0102     TH1D* h_B0_hit_energy_deposit = new TH1D("B0_tracker_hit_energy_deposit", ";Deposited Energy [keV]", 100, 0.0, 500.0);
0103
0104     //B0 EMCAL clusters
0105     TH2D* h_B0_emcal_occupancy_map = new TH2D("B0_emcal_occupancy_map", ";hit x [mm];hit y [mm]", 100, -400, 0, 100, -170, 170);
0106     TH1D* h_B0_emcal_cluster_energy = new TH1D("B0_emcal_cluster_energy", ";Cluster Energy [GeV]", 100, 0.0, 100.0);
0107
0108     //Reconstructed tracks (for usage with B0 too!!)
0109     TH1D* h_px_reco_track = new TH1D("px_reco_track", ";p_{x} [GeV/c]", 100, -10.0, 10.0);
0110     TH1D* h_py_reco_track = new TH1D("py_reco_track", ";p_{y} [GeV/c]", 100, -10.0, 10.0);
0111     TH1D* h_pt_reco_track = new TH1D("pt_reco_track", ";p_{t} [GeV/c]", 150, 0.0, 2.0);
0112     TH1D* h_pz_reco_track = new TH1D("pz_reco_track", ";p_{z} [GeV/c]", 100, 0.0, 320.0);
0113     TH1D* h_pt_squared_B0 = new TH1D("pt_squared_B0", "; Momentum Transfer, -t [GeV^{2}]", 100, 0.0, 2.0);
0114
0115     //ZDC EMCAL clusters
0116     TH2D* h_ZDC_emcal_occupancy_map = new TH2D("ZDC_emcal_occupancy_map", "ZDC_emcal_occupancy_map", 100, -1150, -1050, 100, -60, 60);
0117     TH1D* h_ZDC_emcal_cluster_energy = new TH1D("ZDC_emcal_cluster_energy", ";Cluster Energy [GeV]", 100, 0.0, 100.0);
0118
0119     //B0 momentum resolution
0120
0121     TH1D* h_b0_pt_resolution = new TH1D("b0_pt_resolution", ";#Delta p_{T} [GeV/c]", 100, -2.0, 2.0);
0122     TH2D* h_b0_pt_resolution_percent = new TH2D("b0_deltaPt_over_pt_vs_pt", ";P_{T, MC} [GeV/c]; #Delta p_{T}/p_{T, MC} [percent/100]", 100, 0.0, 2.5, 100, -1.0, 1.0);
0123     TH2D* h_b0_p_resolution_percent = new TH2D("b0_deltaP_over_p_vs_p", ";Three-Momentum,  p_{MC} [GeV/c]; #Delta p/p_{MC} [percent/100]", 100, 0.0, 20.0, 100, -1.0, 1.0);
0124     TH2D* h_b0_px_resolution_percent = new TH2D("b0_deltaPx_over_px_vs_px", ";  P_{x, MC} [GeV/c]; #Delta p_{x}/p_{x, MC} [percent/100]", 100, -10.0, 10.0, 100, -1.0, 1.0);
0125     TH2D* h_b0_py_resolution_percent = new TH2D("b0_deltaPy_over_py_vs_py", ";  p_{y, MC} [GeV/c]; #Delta p_{y}/p_{y, MC} [percent/100]", 100, -10.0, 10.0, 100, -1.0, 1.0);
0126     TH2D* h_b0_pz_resolution_percent = new TH2D("b0_deltaPz_over_pz_vs_pz", ";  p_{z, MC} [GeV/c]; #Delta p_{z}/p_{z, MC} [percent/100]", 100, 0.0, 320.0, 100, -1.0, 1.0);
0127
0128     TH1D* h_b0_extracted_pt_resolution;
0129     TH1D* h_b0_extracted_p_resolution;
0130
0131     TH1D* h_extracted_px_resolution;
0132     TH1D* h_extracted_py_resolution;
0133     TH1D* h_extracted_pz_resolution;
0134
0135     //RP momentum resolution
0136
0137
0138     TH1D* h_RP_pt_resolution = new TH1D("RP_pt_resolution", ";#Delta p_{T} [GeV/c]", 100, -2.0, 2.0);
0139     TH2D* h_RP_pt_resolution_percent = new TH2D("RP_deltaPt_over_pt_vs_pt", ";P_{T, MC} [GeV/c]; #Delta p_{T}/p_{T, MC} [percent/100]", 100, 0.0, 2.5, 100, -1.0, 1.0);
0140     TH2D* h_RP_p_resolution_percent = new TH2D("RP_deltaP_over_p_vs_p", ";Three-Momentum,  p_{MC} [GeV/c]; #Delta p/p_{MC} [percent/100]", 100, 0.0, 320.0, 100, -1.0, 1.0);
0141     TH2D* h_RP_px_resolution_percent = new TH2D("RP_deltaPx_over_px_vs_px", ";  P_{x, MC} [GeV/c]; #Delta p_{x}/p_{x, MC} [percent/100]", 100, -10.0, 10.0, 100, -1.0, 1.0);
0142     TH2D* h_RP_py_resolution_percent = new TH2D("RP_deltaPy_over_py_vs_py", ";  p_{y, MC} [GeV/c]; #Delta p_{y}/p_{y, MC} [percent/100]", 100, -10.0, 10.0, 100, -1.0, 1.0);
0143     TH2D* h_RP_pz_resolution_percent = new TH2D("RP_deltaPz_over_pz_vs_pz", ";  p_{z, MC} [GeV/c]; #Delta p_{z}/p_{z, MC} [percent/100]", 100, 0.0, 320.0, 100, -1.0, 1.0);
0144
0145     TH1D* h_RP_extracted_pt_resolution;
0146     TH1D* h_RP_extracted_p_resolution;
0147
0148     TH1D* h_extracted_RP_px_resolution;
0149     TH1D* h_extracted_RP_py_resolution;
0150     TH1D* h_extracted_RP_pz_resolution;
0151
0152     //forward ECAL information
0153
0154     TH2D* h_forward_emcal_occupancy_map = new TH2D("forward_emcal_occupancy_map", "forward_emcal_occupancy_map", 100, -1000, 1000, 100, -1000, 1000);
0155     TH1D* h_num_forward_emcal_clusters = new TH1D("number_of_forward_EMCAL_clusters_per_event", "number_of_forward_EMCAL_clusters_per_event", 50, 0, 50);
0156
0157     //barrel ECAL information
0158
0159     //TH2D* h_forward_emcal_occupancy_map = new TH2D("forward_emcal_occupancy_map", "forward_emcal_occupancy_map", 100, -1000, 1000, 100, -1000, 1000);
0160     //TH1D* h_num_forward_emcal_clusters = new TH1D("number_of_forward_EMCAL_clusters_per_event", "number_of_forward_EMCAL_clusters_per_event", 50, 0, 50);
0161
0162
0163     int fileCounter = 0;
0164     int iEvent = 0;
0165
0166     while(getline(fileListStream, fileName) && iEvent < 150000){
0167
0168         TString tmp = fileName;
0169
0170         cout << "Input file " << fileCounter << ": " << fileName << endl;
0171
0172         inputRootFile = new TFile(tmp);
0173         if(inputRootFile->IsZombie()){ cout << "MISSING_ROOT_FILE"<< fileName << endl; continue;}
0174
0175         fileCounter++;
0176
0177         TTree * evtTree = (TTree*)inputRootFile->Get("events");
0178
0179         int numEvents = evtTree->GetEntries();
0180
0181         TTreeReader tree_reader(evtTree);       // !the tree reader
0182
0183         //MC particles
0184
0185         TTreeReaderArray<double> mc_px_array = {tree_reader, "MCParticles.momentum.x"};
0186         TTreeReaderArray<double> mc_py_array = {tree_reader, "MCParticles.momentum.y"};
0187         TTreeReaderArray<double> mc_pz_array = {tree_reader, "MCParticles.momentum.z"};
0188         TTreeReaderArray<double> mc_mass_array = {tree_reader, "MCParticles.mass"};
0189         TTreeReaderArray<int> mc_pdg_array = {tree_reader, "MCParticles.PDG"};
0190         TTreeReaderArray<int> mc_genStatus_array = {tree_reader, "MCParticles.generatorStatus"};
0191
0192         TTreeReaderArray<double> truth_px_array = {tree_reader, "MCParticlesHeadOnFrameNoBeamFX.momentum.x"};
0193         TTreeReaderArray<double> truth_py_array = {tree_reader, "MCParticlesHeadOnFrameNoBeamFX.momentum.y"};
0194         TTreeReaderArray<double> truth_pz_array = {tree_reader, "MCParticlesHeadOnFrameNoBeamFX.momentum.z"};
0195         TTreeReaderArray<double> truth_mass_array = {tree_reader, "MCParticlesHeadOnFrameNoBeamFX.mass"};
0196         TTreeReaderArray<int> truth_pdg_array = {tree_reader, "MCParticlesHeadOnFrameNoBeamFX.PDG"};
0197         TTreeReaderArray<int> truth_genStatus_array = {tree_reader, "MCParticlesHeadOnFrameNoBeamFX.generatorStatus"};
0198
0199
0200         //Roman pots -- momentum vector
0201         TTreeReaderArray<float> reco_RP_px = {tree_reader, "ForwardRomanPotRecParticles.momentum.x"};
0202         TTreeReaderArray<float> reco_RP_py = {tree_reader, "ForwardRomanPotRecParticles.momentum.y"};
0203         TTreeReaderArray<float> reco_RP_pz = {tree_reader, "ForwardRomanPotRecParticles.momentum.z"};
0204
0205         //Off-Momentum -- momentum vector
0206         TTreeReaderArray<float> reco_OMD_px = {tree_reader, "ForwardOffMRecParticles.momentum.x"};
0207         TTreeReaderArray<float> reco_OMD_py = {tree_reader, "ForwardOffMRecParticles.momentum.y"};
0208         TTreeReaderArray<float> reco_OMD_pz = {tree_reader, "ForwardOffMRecParticles.momentum.z"};
0209
0210         //hit locations (for debugging)
0211         TTreeReaderArray<float> global_hit_RP_x = {tree_reader, "ForwardRomanPotRecParticles.referencePoint.x"};
0212         TTreeReaderArray<float> global_hit_RP_y = {tree_reader, "ForwardRomanPotRecParticles.referencePoint.y"};
0213         TTreeReaderArray<float> global_hit_RP_z = {tree_reader, "ForwardRomanPotRecParticles.referencePoint.z"};
0214
0215         //hit locations (for debugging)
0216         TTreeReaderArray<float> global_hit_OMD_x = {tree_reader, "ForwardOffMRecParticles.referencePoint.x"};
0217         TTreeReaderArray<float> global_hit_OMD_y = {tree_reader, "ForwardOffMRecParticles.referencePoint.y"};
0218         TTreeReaderArray<float> global_hit_OMD_z = {tree_reader, "ForwardOffMRecParticles.referencePoint.z"};
0219
0220         //b0 tracker hits
0221         TTreeReaderArray<float> b0_hits_x = {tree_reader, "B0TrackerRecHits.position.x"};
0222         TTreeReaderArray<float> b0_hits_y = {tree_reader, "B0TrackerRecHits.position.y"};
0223         TTreeReaderArray<float> b0_hits_z = {tree_reader, "B0TrackerRecHits.position.z"};
0224         TTreeReaderArray<float> b0_MC_momentum_x = {tree_reader, "B0TrackerHits.momentum.x"};
0225         TTreeReaderArray<float> b0_MC_momentum_y = {tree_reader, "B0TrackerHits.momentum.y"};
0226         TTreeReaderArray<float> b0_MC_momentum_z = {tree_reader, "B0TrackerHits.momentum.z"};
0227         TTreeReaderArray<float> b0_hits_eDep = {tree_reader, "B0TrackerRecHits.edep"}; //deposited energy per hit
0228
0229
0230         //b0 EMCAL
0231         TTreeReaderArray<float> b0_cluster_x = {tree_reader, "B0ECalClusters.position.x"};
0232         TTreeReaderArray<float> b0_cluster_y = {tree_reader, "B0ECalClusters.position.y"};
0233         TTreeReaderArray<float> b0_cluster_z = {tree_reader, "B0ECalClusters.position.z"};
0234         TTreeReaderArray<float>  b0_cluster_energy = {tree_reader, "B0ECalClusters.energy"}; //deposited energy in cluster
0235
0236         //reco tracks (where b0 tracks live!!!)
0237         //TTreeReaderArray<float> reco_track_x = {tree_reader, "ReconstructedChargedParticles.momentum.x"};
0238         //TTreeReaderArray<float> reco_track_y = {tree_reader, "ReconstructedChargedParticles.momentum.y"};
0239         //TTreeReaderArray<float> reco_track_z = {tree_reader, "ReconstructedChargedParticles.momentum.z"};
0240         //TTreeReaderArray<int> reco_track_PDG = {tree_reader, "ReconstructedChargedParticles.PDG"};
0241
0242         TTreeReaderArray<float> reco_track_x = {tree_reader, "ReconstructedTruthSeededChargedParticles.momentum.x"};
0243         TTreeReaderArray<float> reco_track_y = {tree_reader, "ReconstructedTruthSeededChargedParticles.momentum.y"};
0244         TTreeReaderArray<float> reco_track_z = {tree_reader, "ReconstructedTruthSeededChargedParticles.momentum.z"};
0245         TTreeReaderArray<int> reco_track_PDG = {tree_reader, "ReconstructedTruthSeededChargedParticles.PDG"};
0246
0247         //forward EMCAL endcap
0248         TTreeReaderArray<float> pECAL_cluster_x = {tree_reader, "EcalEndcapPTruthClusters.position.x"};
0249         TTreeReaderArray<float> pECAL_cluster_y = {tree_reader, "EcalEndcapPTruthClusters.position.y"};
0250         TTreeReaderArray<float> pECAL_cluster_z = {tree_reader, "EcalEndcapPTruthClusters.position.z"};
0251         TTreeReaderArray<float>  pECAL_cluster_energy = {tree_reader, "EcalEndcapPTruthClusters.energy"};
0252
0253         //ZDC EMCAL
0254         //TTreeReaderArray<float> zdc_ecal_cluster_x = {tree_reader, "ZDCEcalClusters.position.x"};
0255         //TTreeReaderArray<float> zdc_ecal_cluster_y = {tree_reader, "ZDCEcalClusters.position.y"};
0256         //TTreeReaderArray<float> zdc_ecal_cluster_z = {tree_reader, "ZDCEcalClusters.position.z"};
0257         //TTreeReaderArray<float>  zdc_ecal_cluster_energy = {tree_reader, "ZDCEcalClusters.energy"}; //deposited energy in cluster
0258
0259
0260         cout << "file has " << evtTree->GetEntries() <<  " events..." << endl;
0261
0262         tree_reader.SetEntriesRange(0, evtTree->GetEntries());
0263         while (tree_reader.Next()) {
0264
0265             //cout << "Reading event: " << iEvent << endl;
0266
0267             //MCParticles
0268
0269             TVector3 mctrk;
0270             TVector3 rptrk;
0271             TVector3 final_mc_track;
0272
0273             double protonMomentum = 0.0;
0274
0275             double maxPt=-99.;
0276
0277             //loop over TRUTH particles where the afterburner effects + crossing angle are compeltely removed
0278             for(int imc = 0; imc < truth_px_array.GetSize(); imc++){
0279                 mctrk.SetXYZ(truth_px_array[imc], truth_py_array[imc], truth_pz_array[imc]);
0280                 final_mc_track.SetXYZ(-999, -999, -999);
0281
0282                 if(truth_pdg_array[imc] == 2212 && truth_genStatus_array[imc] == 1){ //only checking final-state protons here
0283
0284                     if(mctrk.Mag() < 0.0){ continue; }
0285
0286                     h_eta_TRUTH->Fill(mctrk.Eta());
0287                     h_px_TRUTH->Fill(mctrk.Px());
0288                     h_py_TRUTH->Fill(mctrk.Py());
0289                     h_pt_TRUTH->Fill(mctrk.Perp());
0290                     h_pz_TRUTH->Fill(mctrk.Pz());
0291                     h_theta_TRUTH->Fill(1000*mctrk.Theta());
0292                     h_phi_TRUTH->Fill(mctrk.Phi());
0293
0294                     h_pt_squared_TRUTH->Fill(mctrk.Perp()*mctrk.Perp());
0295
0296                     //if(reco_RP_px.GetSize() > 0){h_pt_RomanPots->Fill(mctrk.Perp());} //for acceptance only plot
0297
0298                     //final_mc_track = mctrk;
0299                     //break;
0300                 }
0301             }
0302
0303             //MC Particles loop
0304             for(int imc=0;imc<mc_px_array.GetSize();imc++){
0305                 mctrk.SetXYZ(mc_px_array[imc], mc_py_array[imc], mc_pz_array[imc]);
0306
0307                 //He4 -- 1000020040
0308                 if(mc_pdg_array[imc] == 2212 && mc_genStatus_array[imc] == 1 ){ //only checking for protons here -- change as desired
0309                     mctrk.RotateY(0.025);
0310
0311                     protonMomentum = mctrk.Mag();
0312
0313                     h_eta_MC->Fill(mctrk.Eta());
0314
0315                     h_px_MC->Fill(mctrk.Px());
0316                     h_py_MC->Fill(mctrk.Py());
0317                     h_pt_MC->Fill(mctrk.Perp());
0318                     h_pz_MC->Fill(mctrk.Pz());
0319                     h_theta_MC->Fill(mctrk.Theta()*1000);
0320                     h_phi_MC->Fill(mctrk.Phi());
0321
0322                     h_pt_squared_MC->Fill(mctrk.Perp()*mctrk.Perp());
0323
0324                     //if(reco_RP_px.GetSize() > 0){h_pt_RomanPots->Fill(mctrk.Perp());}//for acceptance only plot
0325                     final_mc_track = mctrk;
0326
0327                 }
0328
0329             }
0330
0331             //roman pots reco tracks
0332             for(int iRPPart = 0; iRPPart < reco_RP_px.GetSize(); iRPPart++){
0333
0334                 TVector3 prec_romanpots(reco_RP_px[iRPPart], reco_RP_py[iRPPart], reco_RP_pz[iRPPart]);
0335
0336                 h_px_RomanPots->Fill(prec_romanpots.Px());
0337                 h_py_RomanPots->Fill(prec_romanpots.Py());
0338                 h_pt_RomanPots->Fill(prec_romanpots.Perp());
0339                 //h_pt_RomanPots->Fill(final_mc_track.Perp());
0340                 h_pz_RomanPots->Fill(prec_romanpots.Pz());
0341                 h_pt_squared_RomanPots->Fill(prec_romanpots.Perp()*prec_romanpots.Perp());
0342
0343
0344                 h_eta_accep->Fill(final_mc_track.Eta());
0345                 h_px_accep->Fill(final_mc_track.Px());
0346                 h_py_accep->Fill(final_mc_track.Py());
0347                 h_pt_accep->Fill(final_mc_track.Perp());
0348                 h_pz_accep->Fill(final_mc_track.Pz());
0349                 h_theta_accep->Fill(final_mc_track.Theta()*1000);
0350                 h_phi_accep->Fill(final_mc_track.Phi());
0351                 h_pt_squared_accep->Fill(final_mc_track.Perp()*final_mc_track.Perp());
0352
0353                 double delPt = prec_romanpots.Perp() - final_mc_track.Perp();
0354                 double delP = prec_romanpots.Mag() - final_mc_track.Mag();
0355
0356                 double delPt_percent = delPt/final_mc_track.Perp();
0357                 double delP_percent = delP/final_mc_track.Mag();
0358
0359                 h_RP_pt_resolution->Fill(delPt_percent);
0360                 h_RP_pt_resolution_percent->Fill(final_mc_track.Perp(), delPt_percent);
0361                 h_RP_p_resolution_percent->Fill(final_mc_track.Mag(), delP_percent);
0362
0363                 double delPx = prec_romanpots.Px() - final_mc_track.Px();
0364                 double delPy = prec_romanpots.Py() - final_mc_track.Py();
0365                 double delPz = prec_romanpots.Pz() - final_mc_track.Pz();
0366
0367                 double delPx_percent = delPx/final_mc_track.Px();
0368                 double delPy_percent = delPy/final_mc_track.Py();
0369                 double delPz_percent = delPz/final_mc_track.Pz();
0370
0371                 h_RP_px_resolution_percent->Fill(final_mc_track.Px(), delPx_percent);
0372                 h_RP_py_resolution_percent->Fill(final_mc_track.Py(), delPy_percent);
0373                 h_RP_pz_resolution_percent->Fill(final_mc_track.Pz(), delPz_percent);
0374
0375                 if(global_hit_RP_z[iRPPart] < 32555.0){ //only the first sensor plane
0376
0377                     h_rp_GLOBAL_occupancy_map->Fill(global_hit_RP_x[iRPPart], global_hit_RP_y[iRPPart]);
0378                     h_rp_LOCAL_occupancy_map->Fill(global_hit_RP_x[iRPPart], global_hit_RP_y[iRPPart]);
0379                 }
0380             }
0381
0382             //OMD reco tracks
0383             for(int iOMDPart = 0; iOMDPart < reco_OMD_px.GetSize(); iOMDPart++){
0384
0385                 TVector3 prec_omd(reco_OMD_px[iOMDPart], reco_OMD_py[iOMDPart], reco_OMD_pz[iOMDPart]);
0386
0387                 h_px_OMD->Fill(prec_omd.Px());
0388                 h_py_OMD->Fill(prec_omd.Py());
0389                 h_pt_OMD->Fill(prec_omd.Perp());
0390                 h_pz_OMD->Fill(prec_omd.Pz());
0391
0392                 h_omd_occupancy_map->Fill(global_hit_OMD_x[iOMDPart], global_hit_OMD_y[iOMDPart]);
0393             }
0394
0395
0396
0397             double hit_x = -9999.;
0398             double hit_y = -9999.;
0399             double hit_z = -9999.;
0400             double hit_deposited_energy = -9999.;
0401
0402             for(int b0cluster = 0; b0cluster < b0_cluster_x.GetSize(); b0cluster++){
0403
0404                 hit_x = b0_cluster_x[b0cluster];
0405                 hit_y = b0_cluster_y[b0cluster];
0406                 hit_z = b0_cluster_z[b0cluster];
0407                 hit_deposited_energy = b0_cluster_energy[b0cluster]*1.246; //poor man's calibration constant, for now
0408
0409                 h_B0_emcal_occupancy_map->Fill(hit_x, hit_y);
0410                 h_B0_emcal_cluster_energy->Fill(hit_deposited_energy);
0411
0412             }
0413
0414
0415             //b0 tracker hits -- for debugging or external tracking
0416             for(int b0hit = 0; b0hit < b0_hits_x.GetSize(); b0hit++){
0417
0418                 hit_x = b0_hits_x[b0hit];
0419                 hit_y = b0_hits_y[b0hit];
0420                 hit_z = b0_hits_z[b0hit];
0421                 hit_deposited_energy = b0_hits_eDep[b0hit]*1e6; //convert GeV --> keV
0422
0423                 h_B0_hit_energy_deposit->Fill(hit_deposited_energy);
0424
0425                 //if(hit_deposited_energy < 10.0){ continue; } //threshold value -- 10 keV, arbitrary for now
0426
0427                 //ACLGAD layout
0428                 if(hit_z > 5700 && hit_z < 5990){ h_B0_occupancy_map_layer_0->Fill(hit_x, hit_y); }
0429                 if(hit_z > 6100 && hit_z < 6200){ h_B0_occupancy_map_layer_1->Fill(hit_x, hit_y); }
0430                 if(hit_z > 6400 && hit_z < 6500){ h_B0_occupancy_map_layer_2->Fill(hit_x, hit_y); }
0431                 if(hit_z > 6700 && hit_z < 6750){ h_B0_occupancy_map_layer_3->Fill(hit_x, hit_y); }
0432
0433             }
0434
0435
0436             //reconstructed tracks with ACTS -- used for B0
0437             //cout << "Event has " << reco_track_x.GetSize() << " reco tracks from ACTS..." << endl;
0438
0439             bool hasB0Track = false;
0440             bool hasHit1 = false;
0441             bool hasHit2 = false;
0442             bool hasHit3 = false;
0443             bool hasHit4 = false;
0444
0445             int b0Hit_layer_1 = 0;
0446
0447
0448             for(int iRecoTrk = 0; iRecoTrk < reco_track_x.GetSize(); iRecoTrk++){
0449
0450                 TVector3 prec_reco_tracks(reco_track_x[iRecoTrk], reco_track_y[iRecoTrk], reco_track_z[iRecoTrk]);
0451
0452                 prec_reco_tracks.RotateY(0.025); //remove crossing angle
0453
0454                 //looking for B0 protons, only -- eta cut is cheap way to do it without associations
0455                 if(reco_track_PDG[iRecoTrk] != 0 || prec_reco_tracks.Eta() < 4.0){ continue; }
0456
0457                 h_px_reco_track->Fill(prec_reco_tracks.Px());
0458                 h_py_reco_track->Fill(prec_reco_tracks.Py());
0459                 h_pt_reco_track->Fill(prec_reco_tracks.Perp());
0460                 h_pz_reco_track->Fill(prec_reco_tracks.Pz());
0461                 h_pt_squared_B0->Fill(prec_reco_tracks.Perp()*prec_reco_tracks.Perp());
0462
0463
0464                 double delPt = prec_reco_tracks.Perp() - final_mc_track.Perp();
0465                 double delP = prec_reco_tracks.Mag() - final_mc_track.Mag();
0466
0467                 double delPt_percent = delPt/final_mc_track.Perp();
0468                 double delP_percent = delP/final_mc_track.Mag();
0469
0470                 h_b0_pt_resolution->Fill(delPt_percent);
0471                 h_b0_pt_resolution_percent->Fill(final_mc_track.Perp(), delPt_percent);
0472                 h_b0_p_resolution_percent->Fill(final_mc_track.Mag(), delP_percent);
0473
0474                 double delPx = prec_reco_tracks.Px() - final_mc_track.Px();
0475                 double delPy = prec_reco_tracks.Py() - final_mc_track.Py();
0476                 double delPz = prec_reco_tracks.Pz() - final_mc_track.Pz();
0477
0478                 double delPx_percent = delPx/final_mc_track.Px();
0479                 double delPy_percent = delPy/final_mc_track.Py();
0480                 double delPz_percent = delPz/final_mc_track.Pz();
0481
0482                 h_b0_px_resolution_percent->Fill(final_mc_track.Px(), delPx_percent);
0483                 h_b0_py_resolution_percent->Fill(final_mc_track.Py(), delPy_percent);
0484                 h_b0_pz_resolution_percent->Fill(final_mc_track.Pz(), delPz_percent);
0485             }
0486
0487
0488
0489             h_num_forward_emcal_clusters->Fill(pECAL_cluster_x.GetSize());
0490
0491             for(int iForwardECALCluster = 0; iForwardECALCluster < pECAL_cluster_x.GetSize(); iForwardECALCluster++){
0492
0493                 h_forward_emcal_occupancy_map->Fill(pECAL_cluster_x[iForwardECALCluster], pECAL_cluster_y[iForwardECALCluster]);
0494
0495
0496             }
0497
0498             iEvent++;
0499
0500         }// event loop
0501
0502         inputRootFile->Close();
0503
0504     }// input file loop
0505
0506
0507     //Calculate detector resolutions here (comment out if not needed)
0508
0509
0510     h_b0_extracted_pt_resolution = extractResolution("b0_extracted_pt_resolution", h_b0_pt_resolution_percent, false);
0511     h_b0_extracted_pt_resolution->GetXaxis()->SetTitle("p_{T, MC} [GeV/c]");
0512     h_b0_extracted_pt_resolution->GetYaxis()->SetTitle("#Delta p_{T}/p_{T, MC} [percent/100]");
0513
0514     h_b0_extracted_p_resolution = extractResolution("b0_extracted_p_resolution", h_b0_p_resolution_percent, false);
0515     h_b0_extracted_p_resolution->GetXaxis()->SetTitle("p_{MC} [GeV/c]");
0516     h_b0_extracted_p_resolution->GetYaxis()->SetTitle("#Delta p/p_{MC} [percent/100]");
0517
0518     h_extracted_px_resolution = extractResolution("track_extracted_px_resolution", h_b0_px_resolution_percent, false);
0519     h_extracted_py_resolution = extractResolution("track_extracted_py_resolution", h_b0_py_resolution_percent, false);
0520     h_extracted_pz_resolution = extractResolution("track_extracted_pz_resolution", h_b0_pz_resolution_percent, false);
0521
0522     h_extracted_px_resolution->GetXaxis()->SetTitle("p_{x, MC} [GeV/c]");
0523     h_extracted_px_resolution->GetYaxis()->SetTitle("#Delta p_{x}/p_{x, MC} [percent/100]");
0524
0525     h_extracted_py_resolution->GetXaxis()->SetTitle("p_{y, MC} [GeV/c]");
0526     h_extracted_py_resolution->GetYaxis()->SetTitle("#Delta p_{y}/p_{y, MC} [percent/100]");
0527
0528     h_extracted_pz_resolution->GetXaxis()->SetTitle("p_{z, MC} [GeV/c]");
0529     h_extracted_pz_resolution->GetYaxis()->SetTitle("#Delta p_{z}/p_{z, MC} [percent/100]");
0530
0531     h_RP_extracted_pt_resolution = extractResolution("RP_extracted_pt_resolution", h_RP_pt_resolution_percent, false);
0532     h_RP_extracted_pt_resolution->GetXaxis()->SetTitle("p_{T, MC} [GeV/c]");
0533     h_RP_extracted_pt_resolution->GetYaxis()->SetTitle("#Delta p_{T}/p_{T, MC} [percent/100]");
0534
0535     h_RP_extracted_p_resolution = extractResolution("RP_extracted_p_resolution", h_RP_p_resolution_percent, false);
0536     h_RP_extracted_p_resolution->GetXaxis()->SetTitle("p_{MC} [GeV/c]");
0537     h_RP_extracted_p_resolution->GetYaxis()->SetTitle("#Delta p/p_{MC} [percent/100]");
0538
0539     h_extracted_RP_px_resolution = extractResolution("track_extracted_RP_px_resolution", h_RP_px_resolution_percent, false);
0540     h_extracted_RP_py_resolution = extractResolution("track_extracted_RP_py_resolution", h_RP_py_resolution_percent, false);
0541     h_extracted_RP_pz_resolution = extractResolution("track_extracted_RP_pz_resolution", h_RP_pz_resolution_percent, false);
0542
0543     h_extracted_RP_px_resolution->GetXaxis()->SetTitle("p_{x, MC} [GeV/c]");
0544     h_extracted_RP_px_resolution->GetYaxis()->SetTitle("#Delta p_{x}/p_{x, MC} [percent/100]");
0545
0546     h_extracted_RP_py_resolution->GetXaxis()->SetTitle("p_{y, MC} [GeV/c]");
0547     h_extracted_RP_py_resolution->GetYaxis()->SetTitle("#Delta p_{y}/p_{y, MC} [percent/100]");
0548
0549     h_extracted_RP_pz_resolution->GetXaxis()->SetTitle("p_{z, MC} [GeV/c]");
0550     h_extracted_RP_pz_resolution->GetYaxis()->SetTitle("#Delta p_{z}/p_{z, MC} [percent/100]");
0551
0552     //for quick checks - just print plots here
0553
0554     /*
0555
0556     TLegend * leg2 = new TLegend(0.55, 0.7, 0.9, 0.9);
0557     leg2->AddEntry(h_pt_squared_MC, "Burned MC, X'ing angle removed", "p");
0558     leg2->AddEntry(h_pt_squared_TRUTH, "Generator Truth", "p");
0559     leg2->AddEntry(h_pt_squared_RomanPots, "Roman Pots FULL reconstruction", "p");
0560     leg2->AddEntry(h_pt_squared_B0, "B0 FULL reconstruction", "p");
0561
0562
0563     TCanvas * cannyCompare2 = new TCanvas("can2", "can2", 600, 500);
0564
0565     //cannyCompare2->Divide(2,1);
0566
0567     cannyCompare2->cd(1)->SetLogy();
0568
0569     h_pt_squared_MC->SetLineColor(kBlack);
0570     h_pt_squared_MC->SetMarkerColor(kBlack);
0571     h_pt_squared_MC->SetMarkerStyle(20);
0572     h_pt_squared_MC->SetStats(0);
0573     h_pt_squared_MC->GetXaxis()->SetTitle("Momentum Transfer, -t [GeV^{2}]");
0574
0575     h_pt_squared_TRUTH->SetLineColor(kRed);
0576     h_pt_squared_TRUTH->SetMarkerColor(kRed);
0577     h_pt_squared_TRUTH->SetMarkerStyle(33);
0578     h_pt_squared_TRUTH->SetStats(0);
0579     h_pt_squared_TRUTH->GetXaxis()->SetTitle("Momentum Transfer, -t [GeV^{2}]");
0580
0581     h_pt_squared_RomanPots->SetLineColor(kGreen + 3);
0582     h_pt_squared_RomanPots->SetMarkerColor(kGreen + 3);
0583     h_pt_squared_RomanPots->SetMarkerStyle(23);
0584     h_pt_squared_RomanPots->SetStats(0);
0585     h_pt_squared_RomanPots->GetXaxis()->SetTitle("Momentum Transfer, -t [GeV^{2}]");
0586
0587     h_pt_squared_B0->SetLineColor(kBlue + 3);
0588     h_pt_squared_B0->SetMarkerColor(kBlue + 3);
0589     h_pt_squared_B0->SetMarkerStyle(23);
0590     h_pt_squared_B0->SetStats(0);
0591     h_pt_squared_B0->GetXaxis()->SetTitle("Momentum Transfer, -t [GeV^{2}]");
0592
0593     h_pt_squared_MC->Draw("EP");
0594     h_pt_squared_TRUTH->Draw("SAME EP");
0595     h_pt_squared_RomanPots->Draw("SAME EP");
0596     h_pt_squared_B0->Draw("SAME EP");
0597
0598     leg2->Draw("SAME");
0599
0600     TBox * RP_outline = new TBox(-128.0, -48.0, 128.0, 48.0);
0601     RP_outline->SetLineColor(kBlack);
0602     RP_outline->SetLineWidth(2);
0603     RP_outline->SetFillStyle(0);
0604
0605
0606     cannyCompare2->cd(2)->SetLogy();
0607
0608     h_pz_MC->SetLineColor(kBlack);
0609     h_pz_MC->SetMarkerColor(kBlack);
0610     h_pz_MC->SetMarkerStyle(20);
0611     h_pz_MC->SetStats(0);
0612     h_pz_MC->GetXaxis()->SetTitle("Longitudinal Momentum, p_{z} [GeV/c]");
0613
0614     h_pz_RomanPots->SetLineColor(kRed);
0615     h_pz_RomanPots->SetMarkerColor(kRed);
0616     h_pz_RomanPots->SetMarkerStyle(33);
0617     h_pz_RomanPots->SetStats(0);
0618     h_pz_RomanPots->GetXaxis()->SetTitle("Longitudinal Momentum, p_{z} [GeV/c]");
0619
0620     h_pz_reco_track->SetLineColor(kGreen + 3);
0621     h_pz_reco_track->SetMarkerColor(kGreen + 3);
0622     h_pz_reco_track->SetMarkerStyle(23);
0623     h_pz_reco_track->SetStats(0);
0624     h_pz_reco_track->GetXaxis()->SetTitle("Longitudinal Momentum, p_{z} [GeV/c]");
0625
0626     h_pz_MC->Draw("EP");
0627     h_pz_RomanPots->Draw("SAME EP");
0628     h_pz_reco_track->Draw("SAME EP");
0629
0630
0631     TString outputPtPlots = "pt_roman_pots_B0_comparisons_";
0632     outputPtPlots = outputPtPlots + outputName + date + run + fileType_PDF;
0633
0634     cannyCompare2->SaveAs(outputPtPlots);
0635
0636     //TCanvas * RP_occupancy_canvas = new TCanvas("canvas_rp", "canvas_rp", 500, 500);
0637
0638     //h_rp_occupancy_map->Draw("COLZ");
0639     //RP_outline->Draw("SAME");
0640
0641     */
0642
0643     SaveHistogramsFromList(*gDirectory->GetList(), outputDir + "/far_forward_dvcs_");
0644
0645     TFile * outputFile = new TFile(outputFileName, "RECREATE");
0646
0647     h_eta_MC->Write();
0648     h_px_MC->Write();
0649     h_py_MC->Write();
0650     h_pt_MC->Write();
0651     h_pz_MC->Write();
0652     h_theta_MC->Write();
0653     h_phi_MC->Write();
0654     h_pt_squared_MC->Write();
0655
0656     h_eta_TRUTH->Write();
0657     h_px_TRUTH->Write();
0658     h_py_TRUTH->Write();
0659     h_pt_TRUTH->Write();
0660     h_pz_TRUTH->Write();
0661     h_theta_TRUTH->Write();
0662     h_phi_TRUTH->Write();
0663     h_pt_squared_TRUTH->Write();
0664
0665     h_eta_accep->Write();
0666     h_px_accep->Write();
0667     h_py_accep->Write();
0668     h_pt_accep->Write();
0669     h_pz_accep->Write();
0670     h_theta_accep->Write();
0671     h_phi_accep->Write();
0672     h_pt_squared_accep->Write();
0673
0674     h_px_RomanPots->Write();
0675     h_py_RomanPots->Write();
0676     h_pt_RomanPots->Write();
0677     h_pz_RomanPots->Write();
0678     h_pt_squared_RomanPots->Write();
0679     h_rp_GLOBAL_occupancy_map->Write();
0680     h_rp_LOCAL_occupancy_map->Write();
0681
0682     h_px_OMD->Write();
0683     h_py_OMD->Write();
0684     h_pt_OMD->Write();
0685     h_pz_OMD->Write();
0686     h_omd_occupancy_map->Write();
0687
0688     h_B0_occupancy_map_layer_0->Write();
0689     h_B0_occupancy_map_layer_1->Write();
0690     h_B0_occupancy_map_layer_2->Write();
0691     h_B0_occupancy_map_layer_3->Write();
0692     h_B0_hit_energy_deposit->Write();
0693
0694     h_B0_emcal_occupancy_map->Write();
0695     h_B0_emcal_cluster_energy->Write();
0696
0697     h_px_reco_track->Write();
0698     h_py_reco_track->Write();
0699     h_pt_reco_track->Write();
0700     h_pz_reco_track->Write();
0701     h_pt_squared_B0->Write();
0702
0703     h_b0_pt_resolution->Write();
0704     h_b0_pt_resolution_percent->Write();
0705     h_b0_extracted_pt_resolution->Write();
0706     h_b0_p_resolution_percent->Write();
0707     h_b0_extracted_p_resolution->Write();
0708
0709     h_extracted_px_resolution->Write();
0710     h_extracted_py_resolution->Write();
0711     h_extracted_pz_resolution->Write();
0712
0713     h_RP_pt_resolution->Write();
0714     h_RP_pt_resolution_percent->Write();
0715     h_RP_extracted_pt_resolution->Write();
0716     h_RP_p_resolution_percent->Write();
0717     h_RP_extracted_p_resolution->Write();
0718
0719     h_extracted_RP_px_resolution->Write();
0720     h_extracted_RP_py_resolution->Write();
0721     h_extracted_RP_pz_resolution->Write();
0722
0723     h_ZDC_emcal_occupancy_map->Write();
0724     h_ZDC_emcal_cluster_energy->Write();
0725
0726     h_forward_emcal_occupancy_map->Write();
0727     h_num_forward_emcal_clusters->Write();
0728
0729     outputFile->Close();
0730
0731
0732
0733     return;
0734
0735 }
0736