tutorial-developing-benchmarks/files/uchannelrho.cxx

0001 #include <cmath>
0002 #include <fstream>
0003 #include <iostream>
0004 #include <string>
0005 #include <vector>
0006 #include <algorithm>
0007 #include <utility>
0008
0009 #include "ROOT/RDataFrame.hxx"
0010 #include <TH1D.h>
0011 #include <TFitResult.h>
0012 #include <TRandom3.h>
0013 #include <TCanvas.h>
0014 #include <TSystem.h>
0015 #include "TFile.h"
0016 #include "TLorentzVector.h"
0017 #include "TLorentzRotation.h"
0018 #include "TVector2.h"
0019 #include "TVector3.h"
0020
0021 #include "fmt/color.h"
0022 #include "fmt/core.h"
0023
0024 #include "edm4eic/InclusiveKinematicsData.h"
0025 #include "edm4eic/ReconstructedParticleData.h"
0026 #include "edm4hep/MCParticleData.h"
0027
0028 #define PI            3.1415926
0029 #define MASS_ELECTRON 0.00051
0030 #define MASS_PROTON   0.93827
0031 #define MASS_PION     0.13957
0032 #define MASS_KAON     0.493667
0033 #define MASS_AU197    183.45406466643374
0034
0035 int uchannelrho(TString rec_file="input.root", TString outputfile="output.root"){
0036     if (gSystem->AccessPathName(rec_file.Data()) != 0) {
0037        // File does not exist
0038        cout<<Form("File %s does not exist.", rec_file.Data())<<endl;
0039        return 0;
0040     }
0041
0042     // read our configuration
0043     auto tree = new TChain("events");
0044     TString name_of_input = (TString) rec_file;
0045     std::cout << "Input file = " << name_of_input << endl;
0046     tree->Add(name_of_input);
0047     TTreeReader tree_reader(tree);       // !the tree reader
0048
0049     //MC-level track attributes
0050     TTreeReaderArray<int>   mc_genStatus_array = {tree_reader, "MCParticles.generatorStatus"};
0051     TTreeReaderArray<float> mc_px_array = {tree_reader, "MCParticles.momentum.x"};
0052     TTreeReaderArray<float> mc_py_array = {tree_reader, "MCParticles.momentum.y"};
0053     TTreeReaderArray<float> mc_pz_array = {tree_reader, "MCParticles.momentum.z"};
0054     TTreeReaderArray<int>   mc_pdg_array= {tree_reader, "MCParticles.PDG"};
0055
0056     //Reco-level track attributes
0057     TTreeReaderArray<float> reco_px_array = {tree_reader, "ReconstructedChargedParticles.momentum.x"};
0058     TTreeReaderArray<float> reco_py_array = {tree_reader, "ReconstructedChargedParticles.momentum.y"};
0059     TTreeReaderArray<float> reco_pz_array = {tree_reader, "ReconstructedChargedParticles.momentum.z"};
0060     TTreeReaderArray<float> reco_charge_array = {tree_reader, "ReconstructedChargedParticles.charge"};
0061     TTreeReaderArray<int>   reco_type     = {tree_reader,"ReconstructedChargedParticles.type"};
0062
0063     TTreeReaderArray<unsigned int> rec_id = {tree_reader, "ReconstructedChargedParticleAssociations.recID"};
0064     TTreeReaderArray<unsigned int> sim_id = {tree_reader, "ReconstructedChargedParticleAssociations.simID"};
0065
0066     TString output_name_dir = outputfile;
0067     cout << "Output file = " << output_name_dir << endl;
0068     TFile* output = new TFile(output_name_dir,"RECREATE");
0069
0070     //Pion reconstruction efficiency histograms
0071     TProfile2D* h_effEtaPtPi  = new TProfile2D("h_effEtaPtPi",";#eta; p_{T}(GeV/c)",50,4,6,30,0,1.6);
0072     TProfile2D* h_effEtaPtPip = new TProfile2D("h_effEtaPtPip",";#eta; p_{T}(GeV/c)",50,4,6,30,0,1.6);
0073     TProfile2D* h_effEtaPtPim = new TProfile2D("h_effEtaPtPim",";#eta; p_{T}(GeV/c)",50,4,6,30,0,1.6);
0074     TProfile2D* h_effPhiEtaPi  = new TProfile2D("h_effPhiEtaPi",";#phi (rad);#eta",50,0,6.4,50,4,6);
0075     TProfile2D* h_effPhiEtaPip = new TProfile2D("h_effPhiEtaPip",";#phi (rad);#eta",50,0,6.4,50,4,6);
0076     TProfile2D* h_effPhiEtaPim = new TProfile2D("h_effPhiEtaPim",";#phi (rad);#eta",50,0,6.4,50,4,6);
0077
0078     //rho vector meson mass
0079     TH1D* h_VM_mass_MC = new TH1D("h_VM_mass_MC",";mass (GeV)",200,0,4);
0080     TH1D* h_VM_mass_MC_etacut = new TH1D("h_VM_mass_MC_etacut",";mass (GeV)",200,0,4);
0081     TH1D* h_VM_mass_REC = new TH1D("h_VM_mass_REC",";mass (GeV)",200,0,4);
0082     TH1D* h_VM_mass_REC_etacut = new TH1D("h_VM_mass_REC_etacut",";mass (GeV)",200,0,4);
0083     TH1D* h_VM_mass_REC_justpions = new TH1D("h_VM_mass_REC_justpions",";mass (GeV)",200,0,4);
0084
0085     cout<<"There are "<<tree->GetEntries()<<" events!!!!!!!"<<endl;
0086     tree_reader.SetEntriesRange(0, tree->GetEntries());
0087     while (tree_reader.Next()) {
0088
0089         TLorentzVector vmMC(0,0,0,0);
0090         TLorentzVector piplusMC(0,0,0,0);
0091         TLorentzVector piminusMC(0,0,0,0);
0092
0093         //MC level
0094         for(unsigned int imc=0;imc<mc_px_array.GetSize();imc++){
0095             TVector3 mctrk(mc_px_array[imc],mc_py_array[imc],mc_pz_array[imc]);
0096             if(mc_pdg_array[imc]!=11) mctrk.RotateY(0.025);//Rotate all non-electrons to hadron beam coordinate system
0097             if(mc_genStatus_array[imc]!=1) continue;
0098             if(mc_pdg_array[imc]==211 && mc_genStatus_array[imc]==1){
0099               piplusMC.SetVectM(mctrk,MASS_PION);
0100             }
0101           if(mc_pdg_array[imc]==-211 && mc_genStatus_array[imc]==1){
0102                 piminusMC.SetVectM(mctrk,MASS_PION);
0103             }
0104         }
0105
0106         vmMC=piplusMC+piminusMC;
0107         if(vmMC.E()!=0){
0108             h_VM_mass_MC->Fill(vmMC.M());
0109             if(piplusMC.Theta()>0.009  && piplusMC.Theta()<0.013 &&
0110                 piminusMC.Theta()>0.009 && piminusMC.Theta()<0.013 ) h_VM_mass_MC_etacut->Fill(vmMC.M());
0111         }
0112
0113         //4-vector for proton reconstructed as if it were a pi+
0114         TLorentzVector protonRECasifpion(0,0,0,0);
0115         //pion 4-vectors
0116         TLorentzVector piplusREC(0,0,0,0);
0117         TLorentzVector piminusREC(0,0,0,0);
0118         //rho reconstruction 4-vector
0119         TLorentzVector vmREC(0,0,0,0);
0120         //rho reconstruction from mis-identified proton 4-vector
0121         TLorentzVector vmRECfromproton(0,0,0,0);
0122
0123         bool isPiMinusFound = false;
0124         bool isPiPlusFound = false;
0125         bool isProtonFound = false;
0126
0127         //track loop
0128         int numpositivetracks = 0;
0129         int failed = 0;
0130         for(unsigned int itrk=0;itrk<reco_pz_array.GetSize();itrk++){
0131             TVector3 trk(reco_px_array[itrk],reco_py_array[itrk],reco_pz_array[itrk]);
0132
0133             //  Rotate in order to account for crossing angle
0134             //  and express coordinates in hadron beam pipe frame
0135             //  This is just a patch, not a final solution.
0136             trk.RotateY(0.025);
0137
0138             if(reco_type[itrk] == -1){
0139                 failed++;
0140                 continue;
0141             }
0142
0143
0144           if(reco_charge_array[itrk]>0){
0145                 numpositivetracks++;
0146               if ((sim_id[itrk - failed]==4 || sim_id[itrk - failed]==5) && reco_charge_array[itrk - failed]==1){
0147                 piplusREC.SetVectM(trk,MASS_PION);
0148                 isPiPlusFound=true;
0149               }
0150              if(sim_id[itrk - failed]==6){
0151                 protonRECasifpion.SetVectM(trk,MASS_PION);
0152                 isProtonFound=true;
0153              }
0154             }
0155           if(reco_charge_array[itrk]<0){
0156             piminusREC.SetVectM(trk,MASS_PION);
0157             if((sim_id[itrk - failed]==4 || sim_id[itrk - failed]==5) && reco_charge_array[itrk - failed]==-1)  isPiMinusFound=true;
0158           }
0159
0160         }
0161
0162         //4vector of VM;
0163         if(piplusREC.E()!=0. && piminusREC.E()!=0.){
0164             vmREC=piplusREC+piminusREC;
0165         }
0166         if(protonRECasifpion.E()!=0. && piminusREC.E()!=0.){
0167           vmRECfromproton=protonRECasifpion+piminusREC;
0168         }
0169
0170         //pion reconstruction efficiency
0171         double thispipeff = (isPiPlusFound) ? 1.0 : 0.0;
0172         double thispimeff = (isPiMinusFound) ? 1.0 : 0.0;
0173         h_effEtaPtPi ->Fill(piplusMC.Eta(), piplusMC.Pt(), thispipeff);
0174         h_effEtaPtPi ->Fill(piminusMC.Eta(),piminusMC.Pt(),thispimeff);
0175         h_effEtaPtPip->Fill(piplusMC.Eta(), piplusMC.Pt(), thispipeff);
0176         h_effEtaPtPim->Fill(piminusMC.Eta(),piminusMC.Pt(),thispimeff);
0177         //
0178         double thispipphi = piplusMC.Phi()>0  ? piplusMC.Phi()  : piplusMC.Phi()+6.2831853;
0179         double thispimphi = piminusMC.Phi()>0 ? piminusMC.Phi() : piminusMC.Phi()+6.2831853;
0180         if(piplusMC.Pt() >0.2) h_effPhiEtaPi ->Fill(thispipphi, piplusMC.Eta(), thispipeff);
0181         if(piminusMC.Pt()>0.2) h_effPhiEtaPi ->Fill(thispimphi,piminusMC.Eta(),thispimeff);
0182         if(piplusMC.Pt() >0.2) h_effPhiEtaPip->Fill(thispipphi, piplusMC.Eta(), thispipeff);
0183         if(piminusMC.Pt()>0.2) h_effPhiEtaPim->Fill(thispimphi,piminusMC.Eta(),thispimeff);
0184         //
0185
0186         //VM rec
0187         if(vmREC.E()==0 && vmRECfromproton.E()==0) continue;
0188         double rho_mass = vmREC.M();
0189         double rho_mass_fromproton = vmRECfromproton.M();
0190         h_VM_mass_REC->Fill(rho_mass);
0191         h_VM_mass_REC->Fill(rho_mass_fromproton);
0192         if(piplusMC.Theta()>0.009  && piplusMC.Theta()<0.013 &&
0193                           piminusMC.Theta()>0.009 && piminusMC.Theta()<0.013 ) h_VM_mass_REC_etacut->Fill(vmREC.M());
0194         if(isPiMinusFound && isPiPlusFound){
0195           h_VM_mass_REC_justpions->Fill(rho_mass);
0196         }
0197     }
0198     output->Write();
0199     output->Close();
0200
0201     return 0;
0202 }