InttSeedTrackPerformance-Jingyu/src/JYInttSeedTrackPerformance.cxx

0001 #ifndef InttSeedTrackPerformance_cxx
0002 #define InttSeedTrackPerformance_cxx
0003
0004 #include "JYInttSeedTracking.h"
0005 #include "JYInttSeedTrackPerformance.h"
0006
0007 #include <TStyle.h>
0008 #include <TCanvas.h>
0009
0010 #include "SPHTracKuma.h"
0011
0012 bool bCaloClu = true;
0013 Int_t strockEvents = 1;
0014
0015 void InttSeedTrackPerformance::Loop(Int_t runNum)
0016 {
0017     if (fChain == 0) return;
0018
0019     // initial the histogram what we need
0020     HistInit();
0021
0022     Long64_t nentries = fChain->GetEntries();
0023     Long64_t nbytes = 0, nb = 0;
0024
0025     bool bTargetEV = false; // wether use the target event?
0026     std::vector<Int_t> vTargetEvents = {};
0027
0028     for (Long64_t jentry=(1000*runNum); jentry<(1000*runNum + 200); jentry++)
0029     {
0030         Long64_t tempJEntry = jentry; // store the jentry
0031         if(bTargetEV) jentry =  vTargetEvents.at(jentry); //用指定的runnumber
0032         pubEvNum = jentry;
0033
0034         Long64_t ientry = LoadTree(jentry);
0035         if (ientry < 0) break;
0036         nb = fChain->GetEntry(jentry);
0037         nbytes += nb;
0038
0039         if(jentry%100 == 0) std::cout << "===== event" << jentry << " =====" << std::endl; // cout to check current runnumber
0040
0041         // Extract the properties of primary particles, and store these data into std::vector<hitStruct> container that can be analyzed
0042         CheckPrimP(m_TruthParticle);
0043         // std::cout << "truth size = " << m_TruthParticle.size() << std::endl;
0044
0045         // get track subcluster information to create hitStruct(R,z,phi,eta), then store to the correspondence vector<hitStruct>
0046         ReadInttHitting(m_fMvtxHits, m_sMvtxHits, m_tMvtxHits, m_iInttHits, m_oInttHits, m_TpcHits);
0047
0048         // get calo information to create hitStruct(R,z,phi,eta), then store to the correspondence vector<hitStruct>
0049         if(bCaloClu) ReadCalCluHitting(m_emcalHits, m_iHCalHits, m_oHCalHits); // calo subcluster information into vector<hitStruct>
0050         else ReadCalHitting(m_emcalHits, m_iHCalHits, m_oHCalHits);            // elseif subtower information into vector<hitStruct>
0051
0052         // std::cout << "== event No." << jentry << "  =========================== "<< std::endl;
0053         // ShowEventInfo();
0054
0055         // vector<vector<hitStruct>>, vector of vector
0056         m_vTruthParticle.push_back(m_TruthParticle);
0057         m_viInttHits.push_back(m_iInttHits);
0058         m_voInttHits.push_back(m_oInttHits);
0059         m_vTpcHits.push_back(m_TpcHits);
0060         m_vemcalHits.push_back(m_emcalHits);
0061         m_viHCalHits.push_back(m_iHCalHits);
0062         m_voHCalHits.push_back(m_oHCalHits);
0063
0064         if(jentry%strockEvents == 0)
0065         {
0066             for(Int_t iEvent = 0; iEvent < strockEvents; iEvent++)
0067             {
0068                 m_iInttHits = m_viInttHits.at(iEvent);
0069                 m_oInttHits = m_voInttHits.at(iEvent);
0070                 m_TpcHits = m_vTpcHits.at(iEvent);
0071                 m_emcalHits = m_vemcalHits.at(iEvent);
0072
0073                 InttSeedTracking *TracKumaContainer = new InttSeedTracking(m_tracks, m_fMvtxHits, m_sMvtxHits, m_tMvtxHits, m_viInttHits.at(iEvent), m_voInttHits.at(iEvent), m_vemcalHits.at(iEvent), m_viHCalHits.at(iEvent), m_voHCalHits.at(iEvent));
0074
0075                 CaloEnergyQA(0, m_vemcalHits.at(iEvent));
0076                 CaloEnergyQA(1, m_viHCalHits.at(iEvent));
0077                 CaloEnergyQA(2, m_voHCalHits.at(iEvent));
0078
0079                 TrackQA(m_vTruthParticle.at(iEvent), m_vemcalHits.at(iEvent));
0080
0081                 m_bDecayEvent = false;
0082                 // ShowEventInfo();
0083             }
0084             AllResetValuesForEvent(); // clear then shrink_to_fit vector of vector
0085         }
0086         PartResetValuesForEvent(); // clear then shrink_to_fit vector of hitStruct
0087
0088         if(bTargetEV) jentry = tempJEntry;
0089     }
0090
0091     WrightHists();
0092
0093     std::cout << "std::vector<Int_t > vTargetEvents = {";
0094     for(Int_t i = 0; i < m_vTargetEvents.size(); i++) std::cout << m_vTargetEvents.at(i) << ", ";
0095     std::cout << "};" << std::endl;
0096 }
0097
0098
0099 void InttSeedTrackPerformance::HistInit()
0100 {
0101     m_HINTTHitMap = new TH2D( "m_HINTTHitMap", "Global Hit Map;x [cm];y [cm]",
0102                     2000, -100., 100., 2000, -100., 100. );
0103
0104     m_HGlobalHitMap = new TH2D( "m_HGlobalHitMap", "Global Hit Map;x [cm];y [cm]",
0105                     5000, -250., 250., 2500, 0., 250. );
0106
0107
0108     m_HMatchCalEVsHighestCalE = new TH1D( "m_HMatchCalEVsHighestCalE", \
0109        "EMCal Energy disparsion map;phi [rad];E [MeV]", \
0110        320, -1.6, 1.6);
0111
0112     m_HdPhiInttdPhiECal = new TH2D( "m_HdPhiInttdPhiECal", \
0113        "EMCal Energy disparsion map;d#it{#phi}_{INTT} [rad];d#it{#phi}_{ECal} [rad]", \
0114           320, -1.6, 1.6, 320, -1.6, 1.6);
0115     m_HTruthPtVsdPhiInttdPhiECal = new TH2D( "m_HTruthPtVsdPhiInttdPhiECal", \
0116        "truth pT vs dPhi (INTT to ECal);#it{p}_{T, truth} [cm];dPhi [rad]", \
0117           150, 0., 15.,  6, -3., 3. );
0118
0119     m_HTruthPtVsdPhiIIntt = new TH2D( "m_HTruthPtVsdPhiIIntt", \
0120        "truth pT vs dPhi (vertex to iINTT);#it{p}_{T, truth} [cm];dPhi [rad]",
0121             150, 0., 15., 6, -3., 3. );
0122     m_HTruthPtVsdPhiOIntt = new TH2D( "m_HTruthPtVsdPhiOIntt", \
0123        "truth pT vs dPhi (vertex to oINTT);#it{p}_{T, truth} [cm];dPhi [rad]",
0124             150, 0., 15.,  6, -3., 3. );
0125
0126
0127     // === pT resolution histograms  ==========
0128     m_HTruthPtVsSagittaPt_rough = new TH3D( "m_HTruthPtVsSagittaPt_rough", \
0129        "truth pT vs delta pt;#it{p}_{T} [GeV/#it{c}];#it{#delta}#it{p}_{T}/#it{p}_{T};#eta",
0130             150, 0., 15., 250, -2., 2., 20, -1., 1);
0131     m_HTruthPtVsSagittaPt_VtxIntt = new TH3D( "m_HTruthPtVsSagittaPt_VtxIntt", \
0132        "truth pT vs delta pt;#it{p}_{T} [GeV/#it{c}];#it{#delta}#it{p}_{T}/#it{p}_{T};#eta",
0133             150, 0., 15., 250, -2., 2., 20, -1., 1);
0134     m_HTruthPtVsSagittaPt_InttEmcal = new TH3D( "m_HTruthPtVsSagittaPt_InttEmcal", \
0135        "truth pT vs delta pt;#it{p}_{T} [GeV/#it{c}];#it{#delta}#it{p}_{T}/#it{p}_{T};#eta",
0136             150, 0., 15., 250, -2., 2., 20, -1., 1);
0137     m_HTruthPtVsSagittaPt_VtxInttEmcal = new TH3D( "m_HTruthPtVsSagittaPt_VtxInttEmcal", \
0138        "truth pT vs delta pt;#it{p}_{T} [GeV/#it{c}];#it{#delta}#it{p}_{T}/#it{p}_{T};#eta",
0139             150, 0., 15., 250, -2., 2., 20, -1., 1);
0140     m_HTruthPtVsSagittaPt_MvtxInttEmcal = new TH3D( "m_HTruthPtVsSagittaPt_MvtxInttEmcal", \
0141        "truth pT vs delta pt;#it{p}_{T} [GeV/#it{c}];#it{#delta}#it{p}_{T}/#it{p}_{T};#eta",
0142             150, 0., 15., 250, -2., 2., 20, -1., 1);
0143     m_HTruthPtVsSagittaPt_VtxMvtxInttEmcal = new TH3D( "m_HTruthPtVsSagittaPt_VtxMvtxInttEmcal", \
0144        "truth pT vs delta pt;#it{p}_{T} [GeV/#it{c}];#it{#delta}#it{p}_{T}/#it{p}_{T};#eta",
0145             150, 0., 15., 250, -2., 2., 20, -1., 1);
0146
0147     m_HTruthPtVsFitFuncPt_IInttOInttEmcal = new TH3D( "m_HTruthPtVsFitFuncPt_IInttOInttEmcal", \
0148        "truth pT vs delta pt;#it{p}_{T} [GeV/#it{c}];#it{#delta}#it{p}_{T}/#it{p}_{T};#eta",
0149             150, 0., 15., 250, -2., 2., 20, -1., 1);
0150
0151     m_HTruthPtVsFitFuncPt_VtxIInttEmcal = new TH3D( "m_HTruthPtVsFitFuncPt_VtxIInttEmcal", \
0152         "truth pT vs delta pt;#it{p}_{T} [GeV/#it{c}];#it{#delta}#it{p}_{T}/#it{p}_{T};#eta",
0153         150, 0., 15., 250, -2., 2., 20, -1., 1);
0154     m_HTruthPtVsFitFuncPt_VtxOInttEmcal = new TH3D( "m_HTruthPtVsFitFuncPt_VtxOInttEmcal", \
0155         "truth pT vs delta pt;#it{p}_{T} [GeV/#it{c}];#it{#delta}#it{p}_{T}/#it{p}_{T};#eta",
0156         150, 0., 15., 250, -2., 2., 20, -1., 1);
0157     m_HTruthPtVsFitFuncPt_VtxInttEmcal = new TH3D( "m_HTruthPtVsFitFuncPt_VtxInttEmcal", \
0158         "truth pT vs delta pt;#it{p}_{T} [GeV/#it{c}];#it{#delta}#it{p}_{T}/#it{p}_{T};#eta",
0159         150, 0., 15., 250, -2., 2., 20, -1., 1);
0160     m_HTruthPtVsFitFuncPt_MvtxIInttEmcal = new TH3D( "m_HTruthPtVsFitFuncPt_MVtxIInttEmcal", \
0161         "truth pT vs delta pt;#it{p}_{T} [GeV/#it{c}];#it{#delta}#it{p}_{T}/#it{p}_{T};#eta",
0162         150, 0., 15., 250, -2., 2., 20, -1., 1);
0163     m_HTruthPtVsFitFuncPt_MvtxOInttEmcal = new TH3D( "m_HTruthPtVsFitFuncPt_MVtxOInttEmcal", \
0164         "truth pT vs delta pt;#it{p}_{T} [GeV/#it{c}];#it{#delta}#it{p}_{T}/#it{p}_{T};#eta",
0165         150, 0., 15., 250, -2., 2., 20, -1., 1);
0166     m_HTruthPtVsFitFuncPt_MvtxInttEmcal = new TH3D( "m_HTruthPtVsFitFuncPt_MVtxInttEmcal", \
0167         "truth pT vs delta pt;#it{p}_{T} [GeV/#it{c}];#it{#delta}#it{p}_{T}/#it{p}_{T};#eta",
0168         150, 0., 15., 250, -2., 2., 20, -1., 1);
0169
0170
0171     // === theta resolution histograms  ==========
0172     m_HTruthPVsTheta_InttEmcal = new TH2D( "m_HTruthPVsTheta_InttEmcal", \
0173        "truth p vs delta pt;#it{p} [GeV/#it{c}];#it{#delta}#it{#theta}",
0174             150, 0., 15., 100, -0.5, 0.5 );
0175     m_HTruthPVsTheta_MvtxInttEmcal = new TH2D( "m_HTruthPVsTheta_MvtxInttEmcal", \
0176        "truth p vs delta p;#it{p} [GeV/#it{c}];#it{#delta}#it{#theta}",
0177             150, 0., 15., 100, -0.5, 0.5 );
0178
0179     // === phi resolution histograms  ==========
0180     m_HTruthPVsPhi_InttEmcal = new TH2D( "m_HTruthPVsPhi_InttEmcal", \
0181        "truth p vs delta pt;#it{p} [GeV/#it{c}];#it{#delta}#it{#phi}",
0182             150, 0., 15., 100, -0.5, 0.5 );
0183     m_HTruthPVsPhi_MvtxInttEmcal = new TH2D( "m_HTruthPVsPhi_MvtxInttEmcal", \
0184        "truth p vs delta p;#it{p} [GeV/#it{c}];#it{#delta}#it{#phi}",
0185             150, 0., 15., 100, -0.5, 0.5 );
0186
0187     // === dp histograms  ==========
0188     m_HTruthPVsRecoP_InttEmcal = new TH2D( "m_HTruthPVsRecoP_InttEmcal", \
0189        "truth p vs delta p;#it{p} [GeV/#it{c}];#it{#delta}#it{p}/#it{p}",
0190         150, 0., 15., 50, -2.5, 2.5 );
0191     m_HTruthPVsRecoP_MvtxInttEmcal = new TH2D( "m_HTruthPVsRecoP_MvtxInttEmcal", \
0192        "truth p vs delta p;#it{p} [GeV/#it{c}];#it{#delta}#it{p}/#it{p}",
0193         150, 0., 15., 50, -2.5, 2.5 );
0194     m_HTruthPVsRecoP_FitFunc = new TH2D( "m_HTruthPVsRecoP_FitFunc", \
0195        "truth p vs delta p;#it{p} [GeV/#it{c}];#it{#delta}#it{p}/#it{p}",
0196         150, 0., 15., 50, -2.5, 2.5 );
0197
0198     // === dE histograms  ==========
0199     m_HDE = new TH3D( "m_HDE", "truth E vs Reco E;#it{E} [MeV];#it{#delta}#it{E}/#it{E};#eta",
0200         200, 0., 20., 50, -2.5, 2.5, 20, -1., 1);
0201
0202
0203     // === E/p histograms  ==========
0204     m_HTruthPVsEOverP_InttEmcal = new TH2D( "m_HTruthPVsEOverP_InttEmcal", \
0205         "truth p vs delta pt;#it{p} [GeV/#it{c}];#it{E}/#it{p}",
0206             150, 0., 15., 60, 0., 3. );
0207     m_HTruthPVsEOverP_MvtxInttEmcal = new TH2D( "m_HTruthPVsEOverP_MvtxInttEmcal", \
0208         "truth p vs delta p;#it{p} [GeV/#it{c}];#it{E}/#it{p}",
0209             150, 0., 15., 60, 0., 3. );
0210     m_HTruthPVsEOverP_FitFunc = new TH2D( "m_HTruthPVsEOverP_FitFunc", \
0211         "truth p vs delta p;#it{p} [GeV/#it{c}];#it{E}/#it{p}",
0212             150, 0., 15., 60, 0., 3. );
0213
0214     // === dVtx histograms  ==========
0215     m_dVtxXY_InttEmcal = new TH2D( "m_dVtxXY_InttEmcal", "vertexReso;d#it{X}_{(reco - truth)} [cm];d#it{Y}_{(reco - truth)} [cm];count", 1000, -50., 50.,  1000, -50., 50.);
0216     m_dVtxXY_MvtxInttEmcal = new TH2D( "m_dVtxXY_MvtxInttEmcal", "vertexReso;d#it{X}_{(reco - truth)} [cm];d#it{Y}_{(reco - truth)} [cm];count", 1000, -50., 50., 1000, -50., 50.);
0217
0218     m_dVtxR_InttEmcal = new TH1D( "m_dVtxR_InttEmcal", "vertexReso;d#it{R}_{(reco - truth)} [cm];count", 500, 0., 50.);
0219     m_dVtxR_MvtxInttEmcal = new TH1D( "m_dVtxR_MvtxInttEmcal", "vertexReso;d#it{R}_{(reco - truth)} [cm];count", 500, 0., 50.);
0220
0221     m_dVtxZ_InttEmcal = new TH1D( "m_dVtxZ_InttEmcal", "vertexReso;d#it{Z}_{(reco - truth)} [cm];count", 1000, -50., 50.);
0222     m_dVtxZ_MvtxInttEmcal = new TH1D( "m_dVtxZ_MvtxInttEmcal", "vertexReso;d#it{Z}_{(reco - truth)} [cm];count", 1000, -50., 50.);
0223
0224     // == s == Calo QA Histograms ==================================
0225     m_HECalPhiVsE  = new TProfile2D( "m_HECalPhiVsE",  "EMCal Energy disparsion map;#phi [rad];#eta;E [MeV]", 320, -1.6, 1.6, 300, -1.5, 1.5);
0226     m_HIHCalPhiVsE = new TProfile2D( "m_HIHCalPhiVsE", "iHCal Energy disparsion map;#phi [rad];#eta;E [MeV]", 320, -1.6, 1.6, 300, -1.5, 1.5);
0227     m_HOHCalPhiVsE = new TProfile2D( "m_HOHCalPhiVsE", "oHCal Energy disparsion map;#phi [rad];#eta;E [MeV]", 320, -1.6, 1.6, 300, -1.5, 1.5);
0228     // == e == Calo QA  Histograms ==================================
0229
0230     // == s == Magnetic shift Histograms ==================================
0231     m_HDPhiVsDStraightVsTruPt = new TH3D("m_HDPhiVsDStraightVsTruPt", "dPhi/dl between vertex and inner INTT;d#it{l} [cm];d#it{phi} [rad];#it{p}_{T} [GeV/#it{c}];", 150, 0., 150.,  100, -0.05, 0.05, 150, 0., 15.);
0232
0233     m_HDPhiDStraight_VtxIIntt = new TH2D("m_HDPhiDStraight_VtxIIntt", "dPhi/dl between vertex and inner INTT;#it{p}_{T} [GeV/#it{c}];d#it{phi}/d#it{l} [rad/cm]", 150, 0., 15.,  100, -0.05, 0.05);
0234     m_HDPhiDStraight1D_VtxIIntt = new TH1D("m_HDPhiDStraight1D_VtxIIntt", "dPhi/dl between vertex and inner INTT;d#it{phi}/d#it{l} [rad/cm];count", 100, -0.05, 0.05);
0235
0236     m_HDPhiDStraight_IInttOIntt = new TH2D("m_HDPhiDStraight_IInttOIntt", "dPhi/dl between vertex and inner INTT;#it{p}_{T} [GeV/#it{c}];d#it{phi}/d#it{l} [rad/cm]", 150, 0., 15.,  100, -0.05, 0.05);
0237     m_HDPhiDStraight1D_IInttOIntt = new TH1D("m_HDPhiDStraight1D_IInttOIntt", "dPhi/dl between vertex and inner INTT;d#it{phi}/d#it{l} [rad/cm];count", 100, -0.05, 0.05);
0238
0239     m_HTruthPtVsDdPhiddStraight_IInttOInttEmcal = new TH2D("m_HTruthPtVsDdPhiddStraight_IInttOInttEmcal", "dPhi/dl between vertex and inner INTT;#it{p}_{T} [GeV/#it{c}];d#it{phi}/  d#it{l} [rad/cm]", 150, 0., 15., 100, -0.05, 0.05);
0240
0241     m_HDPhiDStraight_OInttEmcal = new TH2D("m_HDPhiDStraight_OInttEmcal", "dPhi/dl between vertex and inner INTT;#it{p}_{T} [GeV/#it{c}];d#it{phi}/d#it{l} [rad/cm]", 150, 0., 15.,  100, -0.05, 0.05);
0242     m_HDPhiDStraight1D_OInttEmcal = new TH1D("m_HDPhiDStraight1D_OInttEmcal", "dPhi/dl between vertex and inner INTT;d#it{phi}/d#it{l} [rad/cm];count", 100, -0.05, 0.05);
0243
0244
0245     m_HDLVsTruthPtVsEta_VtxIIntt = new TH3D( "m_HDLVsTruthPtVsEta_VtxIIntt", \
0246         "truth pT vs shift distanse (L);L [cm];#it{p}_{T} [GeV/#it{c}];#eta",
0247         1200, -10., 10., 1000, -50., 50., 20, -1., 1);
0248     m_HDL1D_VtxIIntt = new TH1D( "m_HDL1D_VtxIIntt", \
0249         "truth pT vs shift distanse (L);L [cm]; count", 200, -10., 10.);
0250     m_HDLVsTruthPtVsEta_IInttOIntt = new TH3D( "m_HDLVsTruthPtVsEta_IInttOIntt", \
0251         "truth pT vs shift distanse (L);L [cm];#it{p}_{T} [GeV/#it{c}];#eta",
0252         200, -10., 10., 1000, -50., 50., 20, -1., 1);
0253     m_HDL1D_IInttOIntt = new TH1D( "m_HDL1D_IInttOIntt", \
0254         "truth pT vs shift distanse (L);L [cm]; count", 200, -10., 10.);
0255     m_HDLVsTruthPtVsEta_OInttEmcal = new TH3D( "m_HDLVsTruthPtVsEta_OInttEmcal", \
0256         "truth pT vs shift distanse (L);L [cm];#it{p}_{T} [GeV/#it{c}];#eta",
0257         1000, -50., 50., 1000, -50., 50., 20, -1., 1);
0258     m_HDL1D_OInttEmcal = new TH1D( "m_HDL1D_OInttEmcal", \
0259         "truth pT vs shift distanse (L);L [cm]; count", 1000, -50., 50.);
0260
0261
0262     m_HDPhiVsTruthPtVsEta_VtxIIntt = new TH3D( "m_HDPhiVsTruthPtVsEta_VtxIIntt", \
0263         "truth pT vs shift distanse (L);#phi [rad];#it{p}_{T} [GeV/#it{c}];#eta",
0264         200, -1., 1., 150, 0., 15., 10, -1., 1);
0265     m_HDPhiVsTruthPtVsEta_IInttOIntt = new TH3D( "m_HDPhiVsTruthPtVsEta_IInttOIntt", \
0266         "truth pT vs shift distanse (L);#phi [rad];#it{p}_{T} [GeV/#it{c}];#eta",
0267         200, -1., 1., 150, 0., 15., 10, -1., 1);
0268     m_HDPhiVsTruthPtVsEta_OInttEmcal = new TH3D( "m_HDPhiVsTruthPtVsEta_OInttEmcal", \
0269         "truth pT vs shift distanse (L);#phi [rad];#it{p}_{T} [GeV/#it{c}];#eta",
0270         200, -1., 1., 150, 0., 15., 10, -1., 1);
0271
0272     m_HDPhiVsTruthPtVsEta_IInttOInttEmcal = new TH3D( "m_HDPhiVsTruthPtVsEta_IInttOInttEmcal", \
0273         "truth pT vs shift distanse (L);#phi [rad];#it{p}_{T} [GeV/#it{c}];#eta",
0274         200, -1., 1., 150, 0., 15., 10, -1., 1);
0275
0276     m_HDPhiVsTruthPtVsEta_VtxIInttEmcal = new TH3D( "m_HDPhiVsTruthPtVsEta_VtxIInttEmcal", \
0277         "truth pT vs shift distanse (L);#phi [rad];#it{p}_{T} [GeV/#it{c}];#eta",
0278         200, -1., 1., 150, 0., 15., 10, -1., 1);
0279     m_HDPhiVsTruthPtVsEta_VtxOInttEmcal = new TH3D( "m_HDPhiVsTruthPtVsEta_VtxOInttEmcal", \
0280         "truth pT vs shift distanse (L);#phi [rad];#it{p}_{T} [GeV/#it{c}];#eta",
0281         200, -1., 1., 150, 0., 15., 10, -1., 1);
0282     m_HDPhiVsTruthPtVsEta_VtxInttEmcal = new TH3D( "m_HDPhiVsTruthPtVsEta_VtxInttEmcal", \
0283         "truth pT vs shift distanse (L);#phi [rad];#it{p}_{T} [GeV/#it{c}];#eta",
0284         200, -1., 1., 150, 0., 15., 10, -1., 1);
0285
0286     m_HDPhiVsTruthPtVsEta_MvtxIInttEmcal = new TH3D( "m_HDPhiVsTruthPtVsEta_MvtxIInttEmcal", \
0287         "truth pT vs shift distanse (L);#phi [rad];#it{p}_{T} [GeV/#it{c}];#eta",
0288         200, -1., 1., 150, 0., 15., 10, -1., 1);
0289     m_HDPhiVsTruthPtVsEta_MvtxOInttEmcal = new TH3D( "m_HDPhiVsTruthPtVsEta_MvtxOInttEmcal", \
0290         "truth pT vs shift distanse (L);#phi [rad];#it{p}_{T} [GeV/#it{c}];#eta",
0291         200, -1., 1., 150, 0., 15., 10, -1., 1);
0292     m_HDPhiVsTruthPtVsEta_MvtxInttEmcal = new TH3D( "m_HDPhiVsTruthPtVsEta_MvtxInttEmcal", \
0293         "truth pT vs shift distanse (L);#phi [rad];#it{p}_{T} [GeV/#it{c}];#eta",
0294         200, -1., 1., 150, 0., 15., 10, -1., 1);
0295
0296
0297     m_HDPhiVsTruthPtProf_IInttOInttEmcal = new TProfile( "m_HDPhiVsTruthPtProf_IInttOInttEmcal", \
0298         "truth pT vs shift distanse (L);#phi [rad];#it{p}_{T} [GeV/#it{c}];", 200, -1., 1.);
0299
0300     m_HDPhiVsTruthPtProf_VtxIIntt = new TProfile( "m_HDPhiVsTruthPtProf_VtxIIntt", \
0301         "truth pT vs shift distanse (L);#phi [rad];#it{p}_{T} [GeV/#it{c}];", 200, -1., 1.);
0302     m_HDPhiVsTruthPtProf_IInttOIntt = new TProfile( "m_HDPhiVsTruthPtProf_IInttOIntt", \
0303         "truth pT vs shift distanse (L);#phi [rad];#it{p}_{T} [GeV/#it{c}]", 200, -1., 1.);
0304     m_HDPhiVsTruthPtProf_OInttEmcal = new TProfile( "m_HDPhiVsTruthPtProf_OInttEmcal", \
0305         "truth pT vs shift distanse (L);#phi [rad];#it{p}_{T} [GeV/#it{c}]", 200, -1., 1.);
0306
0307     m_HDPhiVsTruthPtProf_VtxIInttEmcal = new TProfile( "m_HDPhiVsTruthPtProf_VtxIInttEmcal", \
0308         "truth pT vs shift distanse (L);#phi [rad];#it{p}_{T} [GeV/#it{c}]", 200, -1., 1.);
0309     m_HDPhiVsTruthPtProf_VtxOInttEmcal = new TProfile( "m_HDPhiVsTruthPtProf_VtxOInttEmcal", \
0310         "truth pT vs shift distanse (L);#phi [rad];#it{p}_{T} [GeV/#it{c}]",200, -1., 1.);
0311     m_HDPhiVsTruthPtProf_VtxInttEmcal = new TProfile( "m_HDPhiVsTruthPtProf_VtxInttEmcal", \
0312         "truth pT vs shift distanse (L);#phi [rad];#it{p}_{T} [GeV/#it{c}]", 200, -1., 1.);
0313
0314     m_HDPhiVsTruthPtProf_MvtxIInttEmcal = new TProfile( "m_HDPhiVsTruthPtProf_MvtxIInttEmcal", \
0315         "truth pT vs shift distanse (L);#phi [rad];#it{p}_{T} [GeV/#it{c}]",200, -1., 1.);
0316     m_HDPhiVsTruthPtProf_MvtxOInttEmcal = new TProfile( "m_HDPhiVsTruthPtProf_MVtxOInttEmcal", \
0317         "truth pT vs shift distanse (L);#phi [rad];#it{p}_{T} [GeV/#it{c}]", 200, -1., 1.);
0318     m_HDPhiVsTruthPtProf_MvtxInttEmcal = new TProfile( "m_HDPhiVsTruthPtProf_MVtxInttEmcal", \
0319         "truth pT vs shift distanse (L);#phi [rad];#it{p}_{T} [GeV/#it{c}]", 200, -1., 1.);
0320
0321     m_HPtEfficiency = new TH1D( "m_HPtEfficiency", \
0322         "reconstrucntion efficinecy for pT;#it{p}_{T} [GeV/#it{c}];Efficiency [%]", 150, 0., 15);
0323     m_HTruTrackNum = new TH1D( "m_HTruTrackNum", \
0324         "reconstrucntion efficinecy for pT;#it{p}_{T} [GeV/#it{c}];Efficiency [%]", 150, 0., 15);
0325 }
0326
0327 // Extract the properties of primary particles, and store these data into std::vector<hitStruct> container that can be analyzed
0328 void InttSeedTrackPerformance::CheckPrimP(std::vector<hitStruct>& vTruthParticle)
0329 {
0330     // 检查粒子数量  number of TruthParticle
0331     Int_t numOfP = PrimaryG4P_PID->size();
0332     if(numOfP == 0) return;
0333
0334     // 提取粒子属性: 对每个粒子，提取属性, pPhi, pEta, pE, pPt
0335     for(Int_t iP = 0; iP < numOfP; iP++)
0336     {
0337         hitStruct truthParticle;
0338         Double_t pPhi = PrimaryG4P_Phi->at(iP);
0339         Double_t pEta = PrimaryG4P_Eta->at(iP);
0340         Double_t pE = PrimaryG4P_E->at(iP);
0341         Double_t pPt = PrimaryG4P_Pt->at(iP);
0342
0343         truthParticle.phi = pPhi;
0344         truthParticle.eta = pEta;
0345         truthParticle.pt = pPt;
0346         truthParticle.energy = pE;
0347
0348         vTruthParticle.push_back(truthParticle);
0349     }
0350 }
0351
0352 // get track subcluster information to create hitStruct(R,z,phi,eta), then store to the correspondence vector<hitStruct>
0353 void InttSeedTrackPerformance::ReadInttHitting(std::vector<hitStruct >& vFMvtxHits,\
0354    std::vector<hitStruct >& vSMvtxHits, std::vector<hitStruct >& vTMvtxHits,\
0355    std::vector<hitStruct >& vIInttHits, std::vector<hitStruct >& vOInttHits,\
0356    std::vector<hitStruct >& vTpcHits)
0357 {
0358     Int_t numOfClu = trk_system->size(); // 获取所有命中点的数量, trk_system - tree branch
0359     if(numOfClu == 0) return;
0360
0361     // 遍历所有命中点：
0362     for(Int_t iClu = 0; iClu < numOfClu; iClu++)
0363     {
0364         Int_t CluSysId = trk_system->at(iClu); // system - mvtx, intt, tpc
0365         Int_t CluLayId = trk_layer->at(iClu);  // layer  - 0-2 , 3-6,  7 ?
0366
0367         // get cluster xyz for calculate R Phi theta(eta)
0368         Double_t CluX = trk_X->at(iClu);
0369         Double_t CluY = trk_Y->at(iClu);
0370         Double_t CluZ = trk_Z->at(iClu);
0371
0372         Double_t CluR = sqrt(CluX*CluX + CluY*CluY);
0373         Double_t CluPhi = std::atan(CluY/CluX); // -pi/2 to pi/2, need to turn to -pi ~ pi (since 2,3rd phase cannot be get from atan)
0374         if((CluPhi < 0)&&(CluX < 0)) CluPhi += TMath::Pi();
0375         else if((CluPhi > 0)&&(CluX < 0)) CluPhi -= TMath::Pi();
0376
0377         Double_t CluTheta = std::atan(CluR/CluZ);// -pi - pi
0378         Double_t CluEta = (CluZ/std::abs(CluZ)) * (- log(std::abs(std::tan(CluTheta/2))));
0379
0380         // why need to xyz -> R phi eta -> xyz ?
0381         CluX = CluR*std::cos(CluPhi);
0382         CluY = CluR*std::sin(CluPhi);
0383         m_HINTTHitMap->Fill(CluX, CluY); // not only INTT? mvtx + intt + tpc
0384
0385         hitStruct cluHit;
0386         cluHit.r = CluR;
0387         cluHit.z = CluZ;
0388         cluHit.phi = CluPhi;
0389         cluHit.eta = CluEta;
0390
0391         if(CluSysId == 0) // mvtx
0392         {
0393             if(CluLayId == 0) vFMvtxHits.push_back(cluHit);
0394             else if(CluLayId == 1) vSMvtxHits.push_back(cluHit);
0395             else if(CluLayId == 2) vTMvtxHits.push_back(cluHit);
0396         }
0397         else if(CluSysId == 1) // intt
0398         {
0399             if((CluLayId == 3) || (CluLayId == 4)) vIInttHits.push_back(cluHit);
0400             else if((CluLayId == 5) || (CluLayId == 6)) vOInttHits.push_back(cluHit);
0401         }
0402         else // tpc
0403         {
0404             vTpcHits.push_back(cluHit);
0405         }
0406     }
0407 }
0408
0409 // get calo subtower information to create hitStruct(R,z,phi,eta), then store to the correspondence vector<hitStruct>
0410 void InttSeedTrackPerformance::ReadCalHitting(std::vector<hitStruct >& vEmcalHits,\
0411    std::vector<hitStruct >& vIHCalHits, std::vector<hitStruct >& vOHcalHits)
0412 {
0413     Int_t numOfCalTower = tower_Phi->size();
0414     if(numOfCalTower == 0) return;
0415
0416     for(Int_t iCalTower = 0; iCalTower < numOfCalTower; iCalTower++)
0417     {
0418         Int_t calId = tower_system->at(iCalTower);
0419
0420         Double_t caloX = tower_X->at(iCalTower);
0421         Double_t caloY = tower_Y->at(iCalTower);
0422         Double_t caloZ = tower_Z->at(iCalTower);
0423         Double_t caloPhi = tower_Phi->at(iCalTower);
0424
0425         Double_t caloEta = tower_Eta->at(iCalTower);
0426         Double_t caloE = tower_edep->at(iCalTower);
0427
0428         hitStruct caloHit;
0429         caloHit.r = std::sqrt(caloX*caloX + caloY*caloY);
0430         caloHit.z = caloZ;
0431         caloHit.phi = caloPhi;
0432         caloHit.eta = caloEta;
0433         caloHit.energy = caloE;
0434
0435         if(calId == 0) vEmcalHits.push_back(caloHit);
0436         else if(calId == 1) vIHCalHits.push_back(caloHit);
0437         else if(calId == 2) vOHcalHits.push_back(caloHit);
0438
0439         // if(caloE > 0.1){
0440         //    std::cout << "calE = " << caloE << ", bin(phi, eta) = " << tower_Phi_bin->at(iCalTower) << ", " << tower_Eta_bin->at(iCalTower) << std::endl;
0441         // }
0442     }
0443 }
0444
0445 // get emcal cluster information to create hitStruct(R,z,phi,eta), then store to the correspondence vector<hitStruct>
0446 void InttSeedTrackPerformance::ReadCalCluHitting(std::vector<hitStruct >& vEmcalHits,\
0447    std::vector<hitStruct >& vIHCalHits, std::vector<hitStruct >& vOHcalHits)
0448 {
0449     Int_t numOfCaloCluster = caloClus_Phi->size();
0450     if(numOfCaloCluster == 0) return;
0451
0452     for(Int_t iCaloCluster = 0; iCaloCluster < numOfCaloCluster; iCaloCluster++)
0453     {
0454         Int_t calId = caloClus_system->at(iCaloCluster); // cal system: 0-emcal, 1-ihcal, 2-ohcal
0455
0456         Double_t caloX = caloClus_X->at(iCaloCluster);
0457         Double_t caloY = caloClus_Y->at(iCaloCluster);
0458         Double_t caloZ = caloClus_Z->at(iCaloCluster);
0459
0460         Double_t caloPhi = caloClus_Phi->at(iCaloCluster);
0461         Double_t caloE = caloClus_edep->at(iCaloCluster);
0462
0463         hitStruct caloHit;
0464         caloHit.r = std::sqrt(caloX*caloX + caloY*caloY);
0465         caloHit.z = caloZ;
0466         caloHit.phi = caloPhi;
0467
0468         Double_t caloTheta = std::atan(caloHit.r/caloHit.z);
0469         Double_t caloEta = (caloZ/std::abs(caloZ)) * (- log(std::abs(std::tan(caloTheta/2))));
0470         caloHit.eta = caloEta;
0471         caloHit.energy = caloE;
0472
0473         if(calId == 0) vEmcalHits.push_back(caloHit);
0474         else if(calId == 1) vIHCalHits.push_back(caloHit);
0475         else if(calId == 2) vOHcalHits.push_back(caloHit);
0476     }
0477 }
0478
0479
0480 // == s == Track QA Functions  ###############################################
0481 void InttSeedTrackPerformance::TrackQA(std::vector<hitStruct> vTruthPs, std::vector<hitStruct> vEmcalHits)
0482 {
0483     std::vector<Int_t > matchiedRecoTrkId = {};
0484     Int_t numOfP = PrimaryG4P_PID->size();
0485     for(Int_t iP = 0; iP < numOfP; iP++)
0486     {
0487         m_HTruTrackNum->Fill(vTruthPs.at(iP).pt);
0488         if(m_tracks.size() == 0) continue;
0489
0490         Int_t trkId = TruRecoMatching(vTruthPs.at(iP), m_tracks, matchiedRecoTrkId);
0491         if(trkId == 9999) continue;
0492
0493         DeltaPtPerform(vTruthPs.at(iP), m_tracks.at(trkId));
0494         EstimateMagneticShift(vTruthPs.at(iP), m_tracks.at(trkId));
0495
0496         Double_t truE = vTruthPs.at(iP).energy;
0497         Double_t recoE = m_tracks.at(trkId).getTrackE();
0498         m_HDE->Fill(truE, (recoE - truE)/truE, vTruthPs.at(iP).eta);
0499         m_HPtEfficiency->Fill(vTruthPs.at(iP).pt);
0500     }
0501 }
0502
0503 Int_t InttSeedTrackPerformance::TruRecoMatching(hitStruct truthP, std::vector<tracKuma > vRecoTrk, std::vector<Int_t > vMatchiedRecoTrkId)
0504 {
0505     Int_t matchiedRecoTrkId = 9999;
0506     Double_t closestDPt = 9999.;
0507     for(Int_t iRecoTrk = 0; iRecoTrk < vRecoTrk.size(); iRecoTrk++)
0508     {
0509         Double_t tempTruTheta = 2*atan(std::exp(-truthP.eta));
0510         Double_t tempRecoTheta = vRecoTrk.at(iRecoTrk).getTrackTheta();
0511         Double_t tempRecoPhi = 0.;
0512         if(vRecoTrk.at(iRecoTrk).getHitIs(4)) tempRecoPhi = vRecoTrk.at(iRecoTrk).getHitPhi(4);
0513         else if(vRecoTrk.at(iRecoTrk).getHitIs(5)) tempRecoPhi = vRecoTrk.at(iRecoTrk).getHitPhi(5);
0514         // std::cout << "recoTheta, truTheta, dTheta = " << tempRecoTheta << ", " << tempTruTheta << ", " << std::abs(tempRecoTheta - tempTruTheta) << std::endl;
0515
0516         if((std::abs(tempRecoTheta - tempTruTheta) < 0.1))
0517         {
0518             // std::cout << "dPhi = " << std::abs(tempRecoPhi - truthP.phi) << std::endl;
0519             if((std::abs(tempRecoPhi - truthP.phi) < 0.3))
0520             {
0521                 // std::cout << "dPt = " << std::abs(vRecoTrk.at(iRecoTrk).getTrackPt() - truthP.pt) << std::endl;
0522                 if(closestDPt > std::abs(vRecoTrk.at(iRecoTrk).getTrackPt() - truthP.pt))
0523                 {
0524                     bool matchedIs = std::find( vMatchiedRecoTrkId.begin(),  vMatchiedRecoTrkId.end(), iRecoTrk) !=  vMatchiedRecoTrkId.end();
0525                     if(matchedIs) continue;
0526                     closestDPt = std::abs(vRecoTrk.at(iRecoTrk).getTrackPt() - truthP.pt);
0527                     matchiedRecoTrkId = matchedIs;
0528                 }
0529             }
0530         }
0531     }
0532     vMatchiedRecoTrkId.push_back(matchiedRecoTrkId);
0533
0534     return matchiedRecoTrkId;
0535 }
0536
0537 void InttSeedTrackPerformance::DeltaPtPerform(hitStruct truthP, tracKuma trk){
0538    InttSeedTracking truckF;
0539
0540    Double_t truthPt = truthP.pt;
0541    Double_t truthEta = truthP.eta;
0542
0543    Double_t sagittaR = 9999.;
0544    Double_t centerX = 9999.;
0545    Double_t centerY = 9999.;
0546
0547    Double_t HitsXY[3][2];
0548    truckF.Set3PointsXY(HitsXY, trk, 0);
0549    truckF.RoughEstiSagittaCenter3Point(sagittaR, centerX, centerY, HitsXY);
0550    Double_t sagittaPt = truckF.CalcSagittaPt(sagittaR);
0551    Double_t tempdPt = (sagittaPt - truthPt)/truthPt;
0552    m_HTruthPtVsSagittaPt_rough->Fill(truthPt, tempdPt, truthEta);
0553
0554    std::vector<Int_t > subDetIds_VtxIntt = {0, 4, 5};
0555    std::vector<Double_t > hitsR_VtxIntt = {};
0556    std::vector<Double_t > hitsPhi_VtxIntt = {};
0557    if(truckF.ReturnHitsRPhiVect(hitsR_VtxIntt, hitsPhi_VtxIntt, subDetIds_VtxIntt, trk)){
0558       Double_t tempHitOInttPhi = trk.getHitPhi(4);
0559       Double_t tempHitEmcalPhi = trk.getHitPhi(6);
0560       truckF.SagittaRByCircleFit(centerX, centerY, sagittaR, hitsR_VtxIntt, hitsPhi_VtxIntt, tempHitOInttPhi, tempHitEmcalPhi);
0561       Double_t recoPt = truckF.CalcSagittaPt(sagittaR);
0562       Double_t dPt = (sagittaPt - truthPt)/truthPt;
0563       m_HTruthPtVsSagittaPt_VtxIntt->Fill(truthPt, dPt, truthEta);
0564    }
0565
0566    // Double_t HitsXY2[3][2];
0567    // truckF.Set3PointsXY(HitsXY2, trk, 0);
0568    // truckF.RoughEstiSagittaCenter3Point(sagittaR, centerX, centerY, HitsXY);
0569
0570    std::vector<Int_t > subDetIds_InttEmcal = {4, 5, 6};
0571    std::vector<Double_t > hitsR_InttEmcal = {};
0572    std::vector<Double_t > hitsPhi_InttEmcal = {};
0573    if(truckF.ReturnHitsRPhiVect(hitsR_InttEmcal, hitsPhi_InttEmcal, subDetIds_InttEmcal, trk)){
0574       Double_t tempHitOInttPhi = trk.getHitPhi(4);
0575       Double_t tempHitEmcalPhi = trk.getHitPhi(6);
0576       truckF.SagittaRByCircleFit(centerX, centerY, sagittaR, hitsR_InttEmcal, hitsPhi_InttEmcal, tempHitOInttPhi, tempHitEmcalPhi);
0577       Double_t recoPt = truckF.CalcSagittaPt(sagittaR);
0578       Double_t dPt = (recoPt - truthPt)/truthPt;
0579       m_HTruthPtVsSagittaPt_InttEmcal->Fill(truthPt, dPt, truthEta);
0580
0581       TrackOtherPropertiesWTruth(truthP, trk, sagittaR, centerX, centerY, recoPt,\
0582             m_HTruthPVsTheta_InttEmcal, m_HTruthPVsPhi_InttEmcal, m_HTruthPVsRecoP_InttEmcal, m_HTruthPVsEOverP_InttEmcal,\
0583             m_dVtxXY_InttEmcal, m_dVtxR_InttEmcal, m_dVtxZ_InttEmcal);
0584    }
0585
0586    std::vector<Int_t > subDetIds_VtxInttEmcal = {0, 4, 5, 6};
0587    std::vector<Double_t > hitsR_VtxInttEmcal = {};
0588    std::vector<Double_t > hitsPhi_VtxInttEmcal = {};
0589    if(truckF.ReturnHitsRPhiVect(hitsR_VtxInttEmcal, hitsPhi_VtxInttEmcal, subDetIds_VtxInttEmcal, trk)){
0590       Double_t tempHitOInttPhi = trk.getHitPhi(4);
0591       Double_t tempHitEmcalPhi = trk.getHitPhi(6);
0592       truckF.SagittaRByCircleFit(centerX, centerY, sagittaR, hitsR_VtxInttEmcal, hitsPhi_VtxInttEmcal, tempHitOInttPhi, tempHitEmcalPhi);
0593       Double_t recoPt = truckF.CalcSagittaPt(sagittaR);
0594       Double_t dPt = (recoPt - truthPt)/truthPt;
0595       m_HTruthPtVsSagittaPt_VtxInttEmcal->Fill(truthPt, dPt, truthEta);
0596    }
0597
0598    std::vector<Int_t > subDetIds_MvtxInttEmcal = {1, 2, 3, 4, 5, 6};
0599    std::vector<Double_t > hitsR_MvtxInttEmcal = {};
0600    std::vector<Double_t > hitsPhi_MvtxInttEmcal = {};
0601    if(truckF.ReturnHitsRPhiVect(hitsR_MvtxInttEmcal, hitsPhi_MvtxInttEmcal, subDetIds_MvtxInttEmcal, trk)){
0602       Double_t tempHitOInttPhi = trk.getHitPhi(4);
0603       Double_t tempHitEmcalPhi = trk.getHitPhi(6);
0604       truckF.SagittaRByCircleFit(centerX, centerY, sagittaR, hitsR_MvtxInttEmcal, hitsPhi_MvtxInttEmcal, tempHitOInttPhi, tempHitEmcalPhi);
0605       Double_t recoPt = truckF.CalcSagittaPt(sagittaR);
0606       Double_t dPt = (recoPt - truthPt)/truthPt;
0607       m_HTruthPtVsSagittaPt_MvtxInttEmcal->Fill(truthPt, dPt, truthEta);
0608       // std::cout << "pubEvNum = " << pubEvNum << std::endl;
0609       // EventJudge(pubEvNum, dPt, -2., 2., false); //CheckumaDaYo!!!
0610
0611       TrackOtherPropertiesWTruth(truthP, trk, sagittaR, centerX, centerY, recoPt,\
0612          m_HTruthPVsTheta_MvtxInttEmcal, m_HTruthPVsPhi_MvtxInttEmcal,\
0613          m_HTruthPVsRecoP_MvtxInttEmcal, m_HTruthPVsEOverP_MvtxInttEmcal,\
0614          m_dVtxXY_MvtxInttEmcal, m_dVtxR_MvtxInttEmcal, m_dVtxZ_MvtxInttEmcal);
0615
0616       // std::cout << "dPt = " << dPt << std::endl;
0617       // CheckumaDaYo!! oooooooo
0618       // ShowTrackInfo(trk, dPt, centerX, centerY, sagittaR);
0619
0620    }
0621
0622    std::vector<Int_t > subDetIds_VtxMvtxInttEmcal = {0, 1, 2, 3, 4, 5, 6};
0623    std::vector<Double_t > hitsR_VtxMvtxInttEmcal = {};
0624    std::vector<Double_t > hitsPhi_VtxMvtxInttEmcal = {};
0625    if(truckF.ReturnHitsRPhiVect(hitsR_VtxMvtxInttEmcal, hitsPhi_VtxMvtxInttEmcal, subDetIds_VtxMvtxInttEmcal, trk)){
0626       Double_t tempHitOInttPhi = trk.getHitPhi(4);
0627       Double_t tempHitEmcalPhi = trk.getHitPhi(6);
0628       truckF.SagittaRByCircleFit(centerX, centerY, sagittaR, hitsR_VtxMvtxInttEmcal, hitsPhi_VtxMvtxInttEmcal, tempHitOInttPhi, tempHitEmcalPhi);
0629       Double_t recoPt = truckF.CalcSagittaPt(sagittaR);
0630       Double_t dPt = (recoPt - truthPt)/truthPt;
0631       m_HTruthPtVsSagittaPt_VtxMvtxInttEmcal->Fill(truthPt, dPt, truthEta);
0632    }
0633
0634    Double_t dPhiOInttEmcal = truckF.dPhiOInttEmcalEsti(trk);
0635    Double_t OriFunTrackPt = truckF.FitFunctionPt(dPhiOInttEmcal);
0636    Double_t dPtOriFun = (OriFunTrackPt - truthPt)/truthPt;
0637    m_HTruthPtVsFitFuncPt_IInttOInttEmcal->Fill(truthPt, dPtOriFun, truthEta);
0638    EventJudge(pubEvNum, dPtOriFun ,-2, -1, true);
0639    // std::cout << "111 dPt = " << dPtOriFun << std::endl;
0640    // CheckumaDaYo!! oooooooo
0641
0642    Double_t dPhiVtxIInttEmcal = truckF.dPhiVtxIInttEmcalEsti(trk);
0643    Double_t OriFunTrackPt_VtxIInttEmcal = truckF.FitFunctionPt_VtxIInttEmcal(dPhiOInttEmcal);
0644    Double_t dPtOriFun_VtxIInttEmcal = (OriFunTrackPt_VtxIInttEmcal - truthPt)/truthPt;
0645    m_HTruthPtVsFitFuncPt_VtxIInttEmcal->Fill(truthPt, dPtOriFun_VtxIInttEmcal, truthEta);
0646
0647    Double_t dPhiVtxOInttEmcal = truckF.dPhiVtxOInttEmcalEsti(trk);
0648    Double_t OriFunTrackPt_VtxOInttEmcal = truckF.FitFunctionPt_VtxOInttEmcal(dPhiVtxOInttEmcal);
0649    Double_t dPtOriFun_VtxOInttEmcal = (OriFunTrackPt_VtxOInttEmcal - truthPt)/truthPt;
0650    m_HTruthPtVsFitFuncPt_VtxOInttEmcal->Fill(truthPt, dPtOriFun_VtxOInttEmcal, truthEta);
0651
0652    Double_t dPhiVtxInttEmcal = truckF.dPhiVtxInttEmcalEsti(trk);
0653    Double_t OriFunTrackPt_VtxInttEmcal = truckF.FitFunctionPt_VtxInttEmcal(dPhiVtxInttEmcal);
0654    Double_t dPtOriFun_VtxInttEmcal = (OriFunTrackPt_VtxInttEmcal - truthPt)/truthPt;
0655    m_HTruthPtVsFitFuncPt_VtxInttEmcal->Fill(truthPt, dPtOriFun_VtxInttEmcal, truthEta);
0656
0657    Double_t dPhiMvtxIInttEmcal = truckF.dPhiMvtxIInttEmcalEsti(trk);
0658    Double_t OriFunTrackPt_MvtxIInttEmcal = truckF.FitFunctionPt_MVtxIInttEmcal(dPhiMvtxIInttEmcal);
0659    Double_t dPtOriFun_MvtxIInttEmcal = (OriFunTrackPt_MvtxIInttEmcal - truthPt)/truthPt;
0660    m_HTruthPtVsFitFuncPt_MvtxIInttEmcal->Fill(truthPt, dPtOriFun_MvtxIInttEmcal, truthEta);
0661
0662    Double_t dPhiMvtxOInttEmcal = truckF.dPhiMvtxOInttEmcalEsti(trk);
0663    Double_t OriFunTrackPt_MvtxOInttEmcal = truckF.FitFunctionPt_MVtxInttEmcal(dPhiMvtxOInttEmcal);
0664    Double_t dPtOriFun_MvtxOInttEmcal = (OriFunTrackPt_MvtxOInttEmcal - truthPt)/truthPt;
0665    m_HTruthPtVsFitFuncPt_MvtxOInttEmcal->Fill(truthPt, dPtOriFun_MvtxOInttEmcal, truthEta);
0666
0667    Double_t dPhiMvtxInttEmcal = truckF.dPhiMvtxInttEmcalEsti(trk);
0668    Double_t OriFunTrackPt_MvtxInttEmcal = truckF.FitFunctionPt_MVtxInttEmcal(dPhiMvtxInttEmcal);
0669    Double_t dPtOriFun_MvtxInttEmcal = (OriFunTrackPt_MvtxInttEmcal - truthPt)/truthPt;
0670    m_HTruthPtVsFitFuncPt_MvtxInttEmcal->Fill(truthPt, dPtOriFun_MvtxInttEmcal, truthEta);
0671 }
0672
0673
0674
0675 void InttSeedTrackPerformance::TrackOtherPropertiesWTruth(hitStruct truthP, tracKuma trk,\
0676    Double_t sagittaR, Double_t centerX, Double_t centerY, Double_t recoPt,\
0677    TH2D* hTruthPVsTheta, TH2D* hTruthPVsPhi, TH2D* hTruthPVsRecoP,TH2D* hTruthPVsEOverP,\
0678    TH2D* hDVtxXY, TH1D* hDVtxR, TH1D* hDVtxZ){
0679    InttSeedTracking truckF;
0680
0681    Double_t truthPt = truthP.pt;
0682    Double_t truthPhi = truthP.phi;
0683    Double_t truthTheta = 2*atan(std::exp(-truthP.eta));
0684    Double_t truVtxX = 0.; //TruthPV_trig_x
0685    Double_t truVtxY = 0.; //TruthPV_trig_y
0686    Double_t truVtxZ = 0.; //TruthPV_trig_y
0687
0688    Double_t recoTheta = truckF.EstimateRecoTheta(trk, 0);
0689
0690    Double_t trackE = trk.getTrackE();
0691    Double_t recoP = recoPt/sin(recoTheta);
0692
0693    Double_t vtxX = 0.;
0694    Double_t vtxY = 0.;
0695    //    truckF.EstiVertex(vtxX, vtxY, sagittaR, centerX, centerY);
0696    truckF.CrossLineCircle(vtxX, vtxY, centerX, centerY, sagittaR);
0697
0698    Double_t vtxR = std::sqrt(vtxX*vtxX + vtxY*vtxY);
0699    Double_t vtxPhi = atan(vtxY/vtxX);
0700    if((vtxPhi < 0)&&(vtxX < 0)) vtxPhi += TMath::Pi();
0701    else if((vtxPhi > 0)&&(vtxY < 0)) vtxPhi -= TMath::Pi();
0702
0703    Double_t trkPhi = -1/std::tan(vtxPhi);
0704    Double_t vtxZ = 0.;
0705    if(trk.getHitIs(6)){
0706       vtxZ = trk.getHitZ(6) - trk.getHitZ(6)/std::tan(recoTheta);
0707    }
0708
0709    Double_t dVtxX = 9999.;
0710    Double_t dVtxY = 9999.;
0711    Double_t dVtxR = 9999.;
0712    Double_t dVtxZ = 9999.;
0713    dVtxX = vtxX - truVtxX;
0714    dVtxY = vtxY - truVtxY;
0715    dVtxR = std::sqrt(dVtxX*dVtxX + dVtxY*dVtxY);
0716    dVtxZ = vtxZ - truVtxZ;
0717    // std::cout << "dvtx r = " << dVtxR << std::endl;
0718
0719    hTruthPVsTheta->Fill(truthPt, recoTheta - truthTheta);
0720    // hTruthPVsRecoP->Fill(recoP, (recoP - )/recoP);
0721    hTruthPVsEOverP->Fill(truthPt, trackE/recoP);
0722    hTruthPVsPhi->Fill(truthPt, trkPhi - truthPhi);
0723
0724    hDVtxXY->Fill(dVtxX, dVtxY);
0725    hDVtxR->Fill(dVtxR);
0726    hDVtxZ->Fill(dVtxZ);
0727 }
0728
0729 void InttSeedTrackPerformance::EstimateMagneticShift(hitStruct truthP, tracKuma trk){
0730    InttSeedTracking truckF;
0731
0732    Double_t truthPt  = truthP.pt;
0733    Double_t truthEta = truthP.eta;
0734    Double_t truthPhi = truthP.phi;
0735
0736    Double_t refPhi = truthPhi;
0737    refPhi = TMath::Pi()/2 - refPhi;
0738    Double_t tempTruPhi = truthPhi + refPhi;
0739
0740
0741    Double_t fMvtxHitR = trk.getHitR(1);
0742    Double_t fMvtxHitPhi = trk.getHitPhi(1) + refPhi;
0743    Double_t sMvtxHitR = trk.getHitR(2);
0744    Double_t sMvtxHitPhi = trk.getHitPhi(2) + refPhi;
0745    Double_t tMvtxHitR = trk.getHitR(3);
0746    Double_t tMvtxHitPhi = trk.getHitPhi(3) + refPhi;
0747    Double_t MvtxHitR = (fMvtxHitR + sMvtxHitR + tMvtxHitR)/3;
0748    Double_t MvtxHitPhi = (fMvtxHitPhi + sMvtxHitPhi + tMvtxHitPhi)/3;
0749
0750    Double_t iInttHitR = trk.getHitR(4);
0751    Double_t iInttHitPhi = trk.getHitPhi(4) + refPhi;
0752    Double_t oInttHitR = trk.getHitR(5);
0753    Double_t oInttHitPhi = trk.getHitPhi(5) + refPhi;
0754    Double_t InttHitR = (iInttHitR + oInttHitR)/2;
0755    Double_t InttHitPhi = (iInttHitPhi + oInttHitPhi)/2;
0756
0757    Double_t emcalHitR = trk.getHitR(6);
0758    Double_t emcalHitPhi = trk.getHitPhi(6) + refPhi;
0759
0760
0761    Double_t dStraight_VtxIIntt = iInttHitR*sin(iInttHitPhi);
0762    Double_t dShiftL_VtxIIntt   = iInttHitR*cos(iInttHitPhi);
0763
0764    Double_t slopeIInttOIntt = 0.;
0765    Double_t sectionIInttOIntt = 0.;
0766    truckF.CalcLineEq(slopeIInttOIntt, sectionIInttOIntt, \
0767       iInttHitR*cos(iInttHitPhi), iInttHitR*sin(iInttHitPhi),\
0768          oInttHitR*cos(oInttHitPhi), oInttHitR*sin(oInttHitPhi));
0769
0770    Double_t crossX1 = 0.;
0771    Double_t crossY1 = 0.;
0772    Double_t slopeIINTT = dStraight_VtxIIntt/dShiftL_VtxIIntt;
0773    Double_t sectionIINTT = - slopeIINTT*iInttHitR*cos(iInttHitPhi) + iInttHitR*sin(iInttHitPhi);
0774    truckF.CalcCrossPoint(crossX1, crossY1, slopeIINTT, sectionIINTT, \
0775       oInttHitR*cos(oInttHitPhi), oInttHitR*sin(oInttHitPhi));
0776    Double_t dShiftL_VtxOIntt   = oInttHitR * cos(oInttHitPhi);
0777    Double_t dShiftL_IInttOIntt = dShiftL_VtxOIntt - dShiftL_VtxIIntt;
0778    Double_t dStraight_IInttOIntt = \
0779    std::sqrt((crossX1 - iInttHitR*cos(iInttHitPhi))*(crossX1 - iInttHitR*cos(iInttHitPhi))\
0780             + (crossY1 - iInttHitR*sin(iInttHitPhi))*(crossY1 - iInttHitR*sin(iInttHitPhi)));
0781
0782
0783    Double_t slopeOInttEmcal = 0.;
0784    Double_t sectionOInttEmcal = 0.;
0785    truckF.CalcLineEq(slopeOInttEmcal, sectionOInttEmcal, oInttHitR*cos(oInttHitPhi), oInttHitR*sin(oInttHitPhi),\
0786       emcalHitR*cos(emcalHitPhi), emcalHitR*sin(emcalHitPhi));
0787    Double_t crossX2 = 0.;
0788    Double_t crossY2 = 0.;
0789    truckF.CalcCrossPoint(crossX2, crossY2, slopeIInttOIntt, sectionIInttOIntt, \
0790       emcalHitR*cos(emcalHitPhi), emcalHitR*sin(emcalHitPhi));
0791    Double_t dStraight_OInttEmcal \
0792       = std::sqrt((crossX2 - oInttHitR*cos(oInttHitPhi))*(crossX2 - oInttHitR * cos(oInttHitPhi))\
0793             + (crossY2 - oInttHitR*sin(oInttHitPhi))*(crossY2 - oInttHitR * sin(oInttHitPhi)));
0794
0795
0796    Double_t dShiftL_VtxEmcal   = emcalHitR*cos(emcalHitPhi);
0797    Double_t dStraight_VtxEmcal = emcalHitR*sin(emcalHitPhi);
0798    Double_t dShiftL_OInttEmcal = dShiftL_VtxEmcal - dShiftL_VtxOIntt;
0799
0800    Double_t dPhi_VtxIIntt = TMath::Pi()/2 - iInttHitPhi;
0801    Double_t temp_dPhi_IInttOIntt = atan(dStraight_IInttOIntt/dShiftL_IInttOIntt);
0802    if(temp_dPhi_IInttOIntt < 0) temp_dPhi_IInttOIntt += TMath::Pi();
0803    Double_t dPhi_IInttOIntt = TMath::Pi()/2 - temp_dPhi_IInttOIntt - dPhi_VtxIIntt;
0804    Double_t temp_dPhi_OInttEmcal = atan(dStraight_OInttEmcal/dShiftL_OInttEmcal);
0805    if(temp_dPhi_OInttEmcal < 0) temp_dPhi_OInttEmcal += TMath::Pi();
0806    Double_t dPhi_OInttEmcal = TMath::Pi()/2 - temp_dPhi_OInttEmcal - dPhi_VtxIIntt - dPhi_IInttOIntt;
0807
0808    Double_t dPhiDStraight_VtxIIntt = dPhi_VtxIIntt/dStraight_VtxIIntt;
0809    m_HDPhiDStraight_VtxIIntt->Fill(truthPt, dPhiDStraight_VtxIIntt);
0810    m_HDPhiDStraight1D_VtxIIntt->Fill(dPhiDStraight_VtxIIntt);
0811
0812    Double_t dPhiDStraight_IInttOIntt = dPhi_IInttOIntt/dStraight_IInttOIntt;
0813    m_HDPhiDStraight_IInttOIntt->Fill(truthPt,  dPhiDStraight_IInttOIntt);
0814    m_HDPhiDStraight1D_IInttOIntt->Fill(dPhiDStraight_IInttOIntt);
0815
0816    Double_t dPhiDStraight_OInttEmcal = dPhi_OInttEmcal/dStraight_OInttEmcal;
0817    m_HDPhiDStraight_OInttEmcal->Fill(truthPt, dPhiDStraight_OInttEmcal);
0818    m_HDPhiDStraight1D_OInttEmcal->Fill(dPhiDStraight_OInttEmcal);
0819
0820    Double_t ddPhiddStraight_IInttOInttEmcal = dPhiDStraight_OInttEmcal - dPhiDStraight_IInttOIntt;
0821    m_HTruthPtVsDdPhiddStraight_IInttOInttEmcal->Fill(truthPt, ddPhiddStraight_IInttOInttEmcal);
0822
0823    m_HDPhiVsDStraightVsTruPt->Fill(dStraight_VtxIIntt, dPhi_VtxIIntt, truthPt);
0824    m_HDPhiVsDStraightVsTruPt->Fill(dStraight_IInttOIntt, dPhi_IInttOIntt, truthPt);
0825    m_HDPhiVsDStraightVsTruPt->Fill(dPhiDStraight_OInttEmcal, dPhi_OInttEmcal, truthPt);
0826
0827    m_HDLVsTruthPtVsEta_VtxIIntt->Fill(dShiftL_VtxIIntt, truthPt, truthEta);
0828    m_HDL1D_VtxIIntt->Fill(dShiftL_VtxIIntt);
0829    m_HDLVsTruthPtVsEta_IInttOIntt->Fill(dShiftL_IInttOIntt, truthPt, truthEta);
0830    m_HDL1D_IInttOIntt->Fill(dShiftL_IInttOIntt);
0831    m_HDLVsTruthPtVsEta_OInttEmcal->Fill(dShiftL_OInttEmcal, truthPt, truthEta);
0832    m_HDL1D_OInttEmcal->Fill(dShiftL_OInttEmcal);
0833
0834
0835    Double_t tempMvtxX = MvtxHitR*std::cos(MvtxHitPhi);
0836    Double_t tempMvtxY = MvtxHitR*std::sin(MvtxHitPhi);
0837    Double_t tempIInttX = iInttHitR*std::cos(iInttHitPhi);
0838    Double_t tempIInttY = iInttHitR*std::sin(iInttHitPhi);
0839    Double_t tempOInttX = oInttHitR*std::cos(oInttHitPhi);
0840    Double_t tempOInttY = oInttHitR*std::sin(oInttHitPhi);
0841    Double_t tempInttX = InttHitR*std::cos(InttHitPhi);
0842    Double_t tempInttY = InttHitR*std::sin(InttHitPhi);
0843    Double_t tempEmcalX = emcalHitR*std::cos(emcalHitPhi);
0844    Double_t tempEmcalY = emcalHitR*std::sin(emcalHitPhi);
0845
0846    Double_t dPhiMvtxIIntt = std::atan((tempIInttY - tempMvtxY)/(tempIInttX - tempMvtxX));
0847    if(dPhiMvtxIIntt < 0) dPhiMvtxIIntt += TMath::Pi();
0848    Double_t dPhiVtxMvtxIIntt = dPhiMvtxIIntt - MvtxHitPhi;
0849
0850    Double_t dPhiMvtxOIntt = std::atan((tempOInttY - tempMvtxY)/(tempOInttX - tempMvtxX));
0851    if(dPhiMvtxOIntt < 0) dPhiMvtxOIntt += TMath::Pi();
0852    Double_t dPhiVtxMvtxOIntt = dPhiMvtxOIntt - MvtxHitPhi;
0853
0854    Double_t dPhiMvtxIntt = std::atan((tempInttY - tempMvtxY)/(tempInttX - tempMvtxX));
0855    if(dPhiMvtxIntt < 0) dPhiMvtxIntt += TMath::Pi();
0856    Double_t dPhiVtxMvtxIntt = dPhiMvtxIntt - MvtxHitPhi;
0857
0858    // Double_t tempHitsXY[3][2];
0859    // truckF.Set3PointsXY(tempHitsXY, trk, 3);
0860    // Double_t tempR = 0.;
0861    // Double_t tempXc = 0.;
0862    // Double_t tempYc = 0.;
0863    // truckF.RoughEstiSagittaCenter3Point(tempR, tempXc, tempYc, tempHitsXY);
0864    // Double_t crossXcal = 0.;
0865    // Double_t crossYcal = 0.;
0866    // Double_t calPhi = truckF.CrossCircleCircle(crossXcal, crossYcal, tempXc, tempYc, tempR, trk.getHitPhi(5));
0867
0868    // if(std::abs(trk.getHitPhi(6) - calPhi) > 0.3) bDecayEvent = true;
0869    // std::cout << "x, y = " << crossXcal << ",  " << crossYcal << std::endl;
0870    // std::cout << "dPhi = " << std::abs(trk.getHitPhi(6) - calPhi) << std::endl;
0871
0872    if(!m_bDecayEvent){
0873       Double_t dPhiIInttEmcal = std::atan((tempEmcalY - tempIInttY)/(tempEmcalX - tempIInttX));
0874       if(dPhiIInttEmcal < 0) dPhiIInttEmcal += TMath::Pi();
0875       Double_t dPhiVtxIInttEmcal = dPhiIInttEmcal - iInttHitPhi;
0876       Double_t dPhiMvtxIInttEmcal = dPhiIInttEmcal - dPhiMvtxIIntt;
0877
0878       Double_t dPhiOInttEmcal = std::atan((tempEmcalY - tempOInttY)/(tempEmcalX - tempOInttX));
0879       if(dPhiOInttEmcal < 0) dPhiOInttEmcal += TMath::Pi();
0880       Double_t dPhiVtxOInttEmcal = dPhiOInttEmcal - oInttHitPhi;
0881       Double_t dPhiMvtxOInttEmcal = dPhiOInttEmcal - dPhiMvtxOIntt;
0882
0883       Double_t dPhiInttEmcal = std::atan((tempEmcalY - tempInttY)/(tempEmcalX - tempInttX));
0884       if(dPhiInttEmcal < 0) dPhiInttEmcal += TMath::Pi();
0885       Double_t dPhiVtxInttEmcal = dPhiInttEmcal - InttHitPhi;
0886       Double_t dPhiMvtxInttEmcal = dPhiInttEmcal - dPhiMvtxIntt;
0887
0888       Double_t dPhiIInttOIntt = std::atan((tempOInttY - tempIInttY)/(tempOInttX - tempIInttX));
0889       if(dPhiIInttOIntt < 0) dPhiIInttOIntt += TMath::Pi();
0890       Double_t dPhi_IInttOInttEmcal = dPhiOInttEmcal - dPhiIInttOIntt;
0891
0892       m_HDPhiVsTruthPtVsEta_VtxIIntt->Fill(dPhi_VtxIIntt, truthPt, truthEta);
0893       m_HDPhiVsTruthPtVsEta_IInttOIntt->Fill(dPhi_IInttOIntt, truthPt, truthEta);
0894       m_HDPhiVsTruthPtVsEta_OInttEmcal->Fill(dPhi_OInttEmcal, truthPt, truthEta);
0895
0896       m_HDPhiVsTruthPtProf_VtxIIntt->Fill(dPhi_VtxIIntt, truthPt);
0897       m_HDPhiVsTruthPtProf_IInttOIntt->Fill(dPhi_IInttOIntt, truthPt);
0898       m_HDPhiVsTruthPtProf_OInttEmcal->Fill(dPhi_OInttEmcal, truthPt);
0899
0900       m_HDPhiVsTruthPtVsEta_IInttOInttEmcal->Fill(dPhi_IInttOInttEmcal, truthPt, truthEta);
0901       m_HDPhiVsTruthPtProf_IInttOInttEmcal->Fill(dPhi_IInttOInttEmcal, truthPt);
0902
0903       m_HDPhiVsTruthPtVsEta_VtxIInttEmcal->Fill(dPhiVtxIInttEmcal, truthPt, truthEta);
0904       m_HDPhiVsTruthPtVsEta_VtxOInttEmcal->Fill(dPhiVtxOInttEmcal, truthPt, truthEta);
0905       m_HDPhiVsTruthPtVsEta_VtxInttEmcal->Fill(dPhiVtxInttEmcal, truthPt, truthEta);
0906
0907       m_HDPhiVsTruthPtVsEta_MvtxIInttEmcal->Fill(dPhiMvtxIInttEmcal, truthPt, truthEta);
0908       m_HDPhiVsTruthPtVsEta_MvtxOInttEmcal->Fill(dPhiMvtxOInttEmcal, truthPt, truthEta);
0909       m_HDPhiVsTruthPtVsEta_MvtxInttEmcal->Fill(dPhiMvtxInttEmcal, truthPt, truthEta);
0910
0911       m_HDPhiVsTruthPtProf_VtxIInttEmcal->Fill(dPhiVtxIInttEmcal, truthPt);
0912       m_HDPhiVsTruthPtProf_VtxOInttEmcal->Fill(dPhiVtxOInttEmcal, truthPt);
0913       m_HDPhiVsTruthPtProf_VtxInttEmcal->Fill(dPhiVtxInttEmcal, truthPt);
0914
0915       m_HDPhiVsTruthPtProf_MvtxIInttEmcal->Fill(dPhiMvtxIInttEmcal, truthPt);
0916       m_HDPhiVsTruthPtProf_MvtxOInttEmcal->Fill(dPhiMvtxOInttEmcal, truthPt);
0917       m_HDPhiVsTruthPtProf_MvtxInttEmcal->Fill(dPhiMvtxInttEmcal, truthPt);
0918    }
0919 }
0920
0921 void InttSeedTrackPerformance::CaloEnergyQA(Int_t CalId, std::vector<hitStruct > calHits)
0922 {
0923     Int_t highestEID = 9999;
0924     Double_t highestE = 0.;
0925     Int_t numOfCalT = calHits.size();
0926     for(Int_t iCalTower = 0; iCalTower < numOfCalT; iCalTower++)
0927     {
0928         if(highestE < calHits.at(iCalTower).energy)
0929         {
0930             highestE = calHits.at(iCalTower).energy;
0931             highestEID = iCalTower;
0932         }
0933     }
0934     if(highestEID == 9999) return;
0935     Double_t refCalPhi = calHits.at(highestEID).phi;
0936     Double_t refCalEta = calHits.at(highestEID).eta;
0937     for(Int_t iCalT = 0; iCalT < numOfCalT; iCalT++)
0938     {
0939         Double_t phi = calHits.at(iCalT).phi - refCalPhi;
0940         Double_t eta = calHits.at(iCalT).eta - refCalEta;
0941         Double_t E = calHits.at(iCalT).energy;
0942         if(CalId == 0) m_HECalPhiVsE->Fill(phi, eta, E);
0943         else if(CalId == 1) m_HIHCalPhiVsE->Fill(phi, eta, E);
0944         else if(CalId == 2) m_HOHCalPhiVsE->Fill(phi, eta, E);
0945     }
0946 }
0947
0948
0949 bool InttSeedTrackPerformance::ParticleDecayCheck(tracKuma trk, std::vector<hitStruct > vEmcalHits){
0950    bool bIsDecay = false;
0951    Int_t numOfCal = vEmcalHits.size();
0952    for(Int_t iCal = 0; iCal < numOfCal; iCal++){
0953       if((vEmcalHits.at(iCal).energy > 0.1)\
0954          && (std::abs(vEmcalHits.at(iCal).phi - trk.getHitPhi(6)) > 0.04)){
0955          bIsDecay = true;
0956       }
0957    }
0958    return bIsDecay;
0959 }
0960 // == e == Track QA Functions  ###############################################
0961
0962
0963
0964 void InttSeedTrackPerformance::ResetValuesForEvent(){
0965    m_tracks.clear();
0966    m_tracks.shrink_to_fit();
0967
0968    m_TruthParticle.clear();
0969    m_TruthParticle.shrink_to_fit();
0970
0971    m_fMvtxHits.clear();
0972    m_sMvtxHits.clear();
0973    m_tMvtxHits.clear();
0974
0975    m_iInttHits.clear();
0976    m_oInttHits.clear();
0977
0978    m_TpcHits.clear();
0979
0980    m_emcalHits.clear();
0981    m_iHCalHits.clear();
0982    m_oHCalHits.clear();
0983
0984    // ==
0985    m_fMvtxHits.shrink_to_fit();
0986    m_sMvtxHits.shrink_to_fit();
0987    m_tMvtxHits.shrink_to_fit();
0988
0989    m_iInttHits.shrink_to_fit();
0990    m_oInttHits.shrink_to_fit();
0991
0992    m_TpcHits.shrink_to_fit();
0993
0994    m_emcalHits.shrink_to_fit();
0995    m_iHCalHits.shrink_to_fit();
0996    m_oHCalHits.shrink_to_fit();
0997
0998 }
0999
1000 void InttSeedTrackPerformance::AllResetValuesForEvent()
1001 {
1002     m_tracks.clear();
1003     m_tracks.shrink_to_fit();
1004
1005     m_vTruthParticle.clear();
1006     m_viInttHits.clear();
1007     m_voInttHits.clear();
1008     m_vTpcHits.clear();
1009     m_vemcalHits.clear();
1010     m_viHCalHits.clear();
1011     m_voHCalHits.clear();
1012
1013     m_vTruthParticle.shrink_to_fit();
1014     m_viInttHits.shrink_to_fit();
1015     m_voInttHits.shrink_to_fit();
1016     m_vTpcHits.shrink_to_fit();
1017     m_vemcalHits.shrink_to_fit();
1018     m_viHCalHits.shrink_to_fit();
1019     m_voHCalHits.shrink_to_fit();
1020
1021     m_TruthParticle.clear();
1022     m_TruthParticle.shrink_to_fit();
1023
1024     m_fMvtxHits.clear();
1025     m_sMvtxHits.clear();
1026     m_tMvtxHits.clear();
1027
1028     m_iInttHits.clear();
1029     m_oInttHits.clear();
1030
1031     m_TpcHits.clear();
1032
1033     m_emcalHits.clear();
1034     m_iHCalHits.clear();
1035     m_oHCalHits.clear();
1036
1037     // ==
1038     m_fMvtxHits.shrink_to_fit();
1039     m_sMvtxHits.shrink_to_fit();
1040     m_tMvtxHits.shrink_to_fit();
1041
1042     m_iInttHits.shrink_to_fit();
1043     m_oInttHits.shrink_to_fit();
1044
1045     m_TpcHits.shrink_to_fit();
1046
1047     m_emcalHits.shrink_to_fit();
1048     m_iHCalHits.shrink_to_fit();
1049     m_oHCalHits.shrink_to_fit();
1050 }
1051
1052 void InttSeedTrackPerformance::PartResetValuesForEvent()
1053 {
1054     m_TruthParticle.clear();
1055     m_TruthParticle.shrink_to_fit();
1056
1057     m_iInttHits.clear();
1058     m_oInttHits.clear();
1059
1060     m_TpcHits.clear();
1061
1062     m_emcalHits.clear();
1063     m_iHCalHits.clear();
1064     m_oHCalHits.clear();
1065
1066     // ==
1067     m_iInttHits.shrink_to_fit();
1068     m_oInttHits.shrink_to_fit();
1069
1070     m_TpcHits.shrink_to_fit();
1071
1072     m_emcalHits.shrink_to_fit();
1073     m_iHCalHits.shrink_to_fit();
1074     m_oHCalHits.shrink_to_fit();
1075 }
1076
1077
1078
1079 void InttSeedTrackPerformance::WrightHists()
1080 {
1081     oFile = new TFile(fOutput.c_str(), "recreate");
1082     oFile->cd();
1083     m_HINTTHitMap->Write();
1084
1085     m_HGlobalHitMap->Write();
1086
1087     m_HECalPhiVsE->Write();
1088     m_HIHCalPhiVsE->Write();
1089     m_HOHCalPhiVsE->Write();
1090
1091     m_HMatchCalEVsHighestCalE->Write();
1092
1093     m_HdPhiInttdPhiECal->Write();
1094     m_HTruthPtVsdPhiInttdPhiECal->Write();
1095
1096     m_HTruthPtVsSagittaPt_rough->Write();
1097     m_HTruthPtVsSagittaPt_VtxIntt->Write();
1098     m_HTruthPtVsSagittaPt_InttEmcal->Write();
1099     m_HTruthPtVsSagittaPt_VtxInttEmcal->Write();
1100     m_HTruthPtVsSagittaPt_MvtxInttEmcal->Write();
1101     m_HTruthPtVsSagittaPt_VtxMvtxInttEmcal->Write();
1102
1103     m_HTruthPtVsFitFuncPt_IInttOInttEmcal->Write();
1104
1105     m_HTruthPtVsFitFuncPt_VtxIInttEmcal->Write();
1106     m_HTruthPtVsFitFuncPt_VtxOInttEmcal->Write();
1107     m_HTruthPtVsFitFuncPt_VtxInttEmcal->Write();
1108
1109     m_HTruthPtVsFitFuncPt_MvtxIInttEmcal->Write();
1110     m_HTruthPtVsFitFuncPt_MvtxOInttEmcal->Write();
1111     m_HTruthPtVsFitFuncPt_MvtxInttEmcal->Write();
1112
1113     m_HTruthPVsRecoP_InttEmcal->Write();
1114     m_HTruthPVsRecoP_MvtxInttEmcal->Write();
1115     m_HTruthPVsRecoP_FitFunc->Write();
1116
1117     m_HTruthPVsEOverP_InttEmcal->Write();
1118     m_HTruthPVsEOverP_MvtxInttEmcal->Write();
1119     m_HTruthPVsEOverP_FitFunc->Write();
1120
1121     m_HTruthPVsTheta_InttEmcal->Write();
1122     m_HTruthPVsTheta_MvtxInttEmcal->Write();
1123
1124     m_HTruthPVsPhi_InttEmcal->Write();
1125     m_HTruthPVsPhi_MvtxInttEmcal->Write();
1126
1127
1128     m_dVtxXY_InttEmcal->Write();
1129     m_dVtxXY_MvtxInttEmcal->Write();
1130
1131     m_dVtxR_InttEmcal->Write();
1132     m_dVtxR_MvtxInttEmcal->Write();
1133
1134     m_dVtxZ_InttEmcal->Write();
1135     m_dVtxZ_MvtxInttEmcal->Write();
1136
1137     m_HDE->Write();
1138
1139     m_HDPhiVsDStraightVsTruPt->Write();
1140
1141     m_HDPhiDStraight_VtxIIntt->Write();
1142     m_HDPhiDStraight1D_VtxIIntt->Write();
1143     m_HDPhiDStraight_IInttOIntt->Write();
1144     m_HDPhiDStraight1D_IInttOIntt->Write();
1145     m_HDPhiDStraight_OInttEmcal->Write();
1146     m_HDPhiDStraight1D_OInttEmcal->Write();
1147     m_HTruthPtVsDdPhiddStraight_IInttOInttEmcal->Write();
1148
1149
1150     m_HDLVsTruthPtVsEta_VtxIIntt->Write();
1151     m_HDL1D_VtxIIntt->Write();
1152     m_HDLVsTruthPtVsEta_IInttOIntt->Write();
1153     m_HDL1D_IInttOIntt->Write();
1154     m_HDLVsTruthPtVsEta_OInttEmcal->Write();
1155     m_HDL1D_OInttEmcal->Write();
1156
1157     m_HDPhiVsTruthPtVsEta_VtxIIntt->Write();
1158     m_HDPhiVsTruthPtVsEta_IInttOIntt->Write();
1159     m_HDPhiVsTruthPtVsEta_OInttEmcal->Write();
1160
1161     m_HDPhiVsTruthPtVsEta_IInttOInttEmcal->Write();
1162
1163     m_HDPhiVsTruthPtVsEta_VtxIInttEmcal->Write();
1164     m_HDPhiVsTruthPtVsEta_VtxOInttEmcal->Write();
1165     m_HDPhiVsTruthPtVsEta_VtxInttEmcal->Write();
1166     m_HDPhiVsTruthPtVsEta_MvtxIInttEmcal->Write();
1167     m_HDPhiVsTruthPtVsEta_MvtxOInttEmcal->Write();
1168     m_HDPhiVsTruthPtVsEta_MvtxInttEmcal->Write();
1169     m_HDPhiVsTruthPtProf_VtxIIntt->Write();
1170
1171     m_HDPhiVsTruthPtProf_IInttOInttEmcal->Write();
1172     m_HDPhiVsTruthPtProf_IInttOIntt->Write();
1173     m_HDPhiVsTruthPtProf_OInttEmcal->Write();
1174     m_HDPhiVsTruthPtProf_VtxIInttEmcal->Write();
1175     m_HDPhiVsTruthPtProf_VtxOInttEmcal->Write();
1176     m_HDPhiVsTruthPtProf_VtxInttEmcal->Write();
1177     m_HDPhiVsTruthPtProf_MvtxIInttEmcal->Write();
1178     m_HDPhiVsTruthPtProf_MvtxOInttEmcal->Write();
1179     m_HDPhiVsTruthPtProf_MvtxInttEmcal->Write();
1180
1181     m_HPtEfficiency->Divide(m_HPtEfficiency, m_HTruTrackNum, 1.0, 1.0, "b");
1182     m_HPtEfficiency->Write();
1183     m_HTruTrackNum->Write();
1184
1185     oFile->Close();
1186
1187     std::cout << "root " << fOutput << std::endl;
1188 }
1189
1190 // MakeClass Original Functions  ###############################################
1191 InttSeedTrackPerformance::InttSeedTrackPerformance(TTree *tree, std::string fInputName, std::string fOutputName, Int_t runNum) : fChain(0)
1192 {
1193 // if parameter tree is not specified (or zero), connect the file
1194 // used to generate this class and read the Tree.
1195    fInput = fInputName;
1196    fOutput = fOutputName + std::to_string(runNum) + ".root";
1197
1198    if (tree == 0) {
1199       TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject(fInput.c_str());
1200       if (!f || !f->IsOpen()) {
1201          f = new TFile(fInput.c_str());
1202       }
1203       f->GetObject("tree",tree);
1204
1205    }
1206    Init(tree);
1207 }
1208
1209 InttSeedTrackPerformance::~InttSeedTrackPerformance()
1210 {
1211    if (!fChain) return;
1212    delete fChain->GetCurrentFile();
1213 }
1214
1215 Int_t InttSeedTrackPerformance::GetEntry(Long64_t entry)
1216 {
1217 // Read contents of entry.
1218    if (!fChain) return 0;
1219    return fChain->GetEntry(entry);
1220 }
1221
1222 Long64_t InttSeedTrackPerformance::LoadTree(Long64_t entry)
1223 {
1224 // Set the environment to read one entry
1225    if (!fChain) return -5;
1226    Long64_t centry = fChain->LoadTree(entry);
1227    if (centry < 0) return centry;
1228    if (fChain->GetTreeNumber() != fCurrent) {
1229       fCurrent = fChain->GetTreeNumber();
1230    }
1231    return centry;
1232 }
1233
1234 void InttSeedTrackPerformance::Init(TTree *tree)
1235 {
1236    // The Init() function is called when the selector needs to initialize
1237    // a new tree or chain. Typically here the branch addresses and branch
1238    // pointers of the tree will be set.
1239    // It is normally not necessary to make changes to the generated
1240    // code, but the routine can be extended by the user if needed.
1241    // Init() will be called many times when running on PROOF
1242    // (once per file to be processed).
1243
1244    // Set object pointer
1245    trk_system = 0;
1246    trk_layer = 0;
1247    trk_adc = 0;
1248    trk_X = 0;
1249    trk_Y = 0;
1250    trk_Z = 0;
1251    trk_size = 0;
1252    trk_phi_size = 0;
1253    trk_Z_size = 0;
1254
1255    tower_system = 0;
1256    tower_X = 0;
1257    tower_Y = 0;
1258    tower_Z = 0;
1259    tower_Eta = 0;
1260    tower_Phi = 0;
1261    tower_Eta_test = 0;
1262    tower_Phi_test = 0;
1263    tower_Eta_bin = 0;
1264    tower_Phi_bin = 0;
1265    tower_edep = 0;
1266
1267    caloClus_system = 0;
1268    caloClus_X = 0;
1269    caloClus_Y = 0;
1270    caloClus_Z = 0;
1271    caloClus_R = 0;
1272    caloClus_Phi = 0;
1273    caloClus_edep = 0;
1274
1275    PrimaryG4P_Pt = 0;
1276    PrimaryG4P_Eta = 0;
1277    PrimaryG4P_Phi = 0;
1278    PrimaryG4P_E = 0;
1279    PrimaryG4P_PID = 0;
1280    PrimaryG4P_isChargeHadron = 0;
1281    // Set branch addresses and branch pointers
1282    if (!tree) return;
1283    fChain = tree;
1284    fCurrent = -1;
1285    fChain->SetMakeClass(1);
1286
1287    fChain->SetBranchAddress("trk_NClus", &trk_NClus, &b_trk_NClus);
1288    fChain->SetBranchAddress("trk_system", &trk_system, &b_trk_system);
1289    fChain->SetBranchAddress("trk_layer", &trk_layer, &b_trk_layer);
1290    fChain->SetBranchAddress("trk_adc", &trk_adc, &b_trk_adc);
1291    fChain->SetBranchAddress("trk_X", &trk_X, &b_trk_X);
1292    fChain->SetBranchAddress("trk_Y", &trk_Y, &b_trk_Y);
1293    fChain->SetBranchAddress("trk_Z", &trk_Z, &b_trk_Z);
1294    fChain->SetBranchAddress("trk_size", &trk_size, &b_trk_size);
1295    fChain->SetBranchAddress("trk_phi_size", &trk_phi_size, &b_trk_phi_size);
1296    fChain->SetBranchAddress("trk_Z_size", &trk_Z_size, &b_trk_Z_size);
1297
1298    fChain->SetBranchAddress("nTowers", &nTowers, &b_nTowers);
1299    fChain->SetBranchAddress("tower_system", &tower_system, &b_tower_system);
1300    fChain->SetBranchAddress("tower_X", &tower_X, &b_tower_X);
1301    fChain->SetBranchAddress("tower_Y", &tower_Y, &b_tower_Y);
1302    fChain->SetBranchAddress("tower_Z", &tower_Z, &b_tower_Z);
1303    fChain->SetBranchAddress("tower_Eta", &tower_Eta, &b_tower_Eta);
1304    fChain->SetBranchAddress("tower_Phi", &tower_Phi, &b_tower_Phi);
1305    fChain->SetBranchAddress("tower_Eta_test", &tower_Eta_test, &b_tower_Eta_test);
1306    fChain->SetBranchAddress("tower_Phi_test", &tower_Phi_test, &b_tower_Phi_test);
1307    fChain->SetBranchAddress("tower_Eta_bin", &tower_Eta_bin, &b_tower_Eta_bin);
1308    fChain->SetBranchAddress("tower_Phi_bin", &tower_Phi_bin, &b_tower_Phi_bin);
1309    fChain->SetBranchAddress("tower_edep", &tower_edep, &b_tower_edep);
1310
1311    //checkumaDaYooo!!
1312    if(bCaloClu){
1313       fChain->SetBranchAddress("nCaloClus", &nCaloClus, &b_nCaloClus);
1314       fChain->SetBranchAddress("caloClus_system", &caloClus_system, &b_caloClus_system);
1315       fChain->SetBranchAddress("caloClus_X", &caloClus_X, &b_caloClus_X);
1316       fChain->SetBranchAddress("caloClus_Y", &caloClus_Y, &b_caloClus_Y);
1317       fChain->SetBranchAddress("caloClus_Z", &caloClus_Z, &b_caloClus_Z);
1318       fChain->SetBranchAddress("caloClus_R", &caloClus_R, &b_caloClus_R);
1319       fChain->SetBranchAddress("caloClus_Phi", &caloClus_Phi, &b_caloClus_Phi);
1320       fChain->SetBranchAddress("caloClus_edep", &caloClus_edep, &b_caloClus_edep);
1321    }
1322
1323
1324    fChain->SetBranchAddress("NTruthVtx", &NTruthVtx, &b_NTruthVtx);
1325    fChain->SetBranchAddress("TruthPV_trig_x", &TruthPV_trig_x, &b_TruthPV_trig_x);
1326    fChain->SetBranchAddress("TruthPV_trig_y", &TruthPV_trig_y, &b_TruthPV_trig_y);
1327    fChain->SetBranchAddress("TruthPV_trig_z", &TruthPV_trig_z, &b_TruthPV_trig_z);
1328    fChain->SetBranchAddress("NPrimaryG4P", &NPrimaryG4P, &b_NPrimaryG4P);
1329    fChain->SetBranchAddress("PrimaryG4P_Pt", &PrimaryG4P_Pt, &b_PrimaryG4P_Pt);
1330    fChain->SetBranchAddress("PrimaryG4P_Eta", &PrimaryG4P_Eta, &b_PrimaryG4P_Eta);
1331    fChain->SetBranchAddress("PrimaryG4P_Phi", &PrimaryG4P_Phi, &b_PrimaryG4P_Phi);
1332    fChain->SetBranchAddress("PrimaryG4P_E", &PrimaryG4P_E, &b_PrimaryG4P_E);
1333    fChain->SetBranchAddress("PrimaryG4P_PID", &PrimaryG4P_PID, &b_PrimaryG4P_PID);
1334    fChain->SetBranchAddress("PrimaryG4P_isChargeHadron", &PrimaryG4P_isChargeHadron, &b_PrimaryG4P_isChargeHadron);
1335
1336 }
1337 // MakeClass Original Functions  ###############################################
1338
1339
1340 void InttSeedTrackPerformance::DrawEventDisplay(Int_t eventId, std::vector<tracKuma > trks){
1341    if(qaEvent != 9999){
1342       TCanvas* c = new TCanvas("c","c: pT",800,850);
1343       c->cd();
1344       TPad* pad1 = new TPad("pad1", "pad1", 0.1, 0.1, 1, 1.0);
1345       pad1->SetBottomMargin(0.15);
1346       pad1->SetLeftMargin(0.15);
1347       pad1->SetRightMargin(0.05);
1348       pad1->SetTopMargin(0.05);
1349       pad1->Draw();
1350       pad1->cd();
1351
1352       TH2D* hHitMapED = new TH2D( "hHitMapED", "hHitMapED", 200, -100, 100, 200, -100, 100);
1353       std::vector<TArc* > trkLines;
1354       for(Int_t iTrk = 0; iTrk < trks.size(); iTrk++){
1355          tracKuma trk = trks.at(iTrk);
1356          for(Int_t iHit = 1; iHit < 7; iHit++){
1357             Double_t tempX = trk.getHitR(iHit)*std::cos(trk.getHitPhi(iHit));
1358             Double_t tempY = trk.getHitR(iHit)*std::sin(trk.getHitPhi(iHit));
1359             hHitMapED->Fill(tempX, tempY);
1360          }
1361
1362          Double_t tempR = trk.getTrackSagR();
1363          Double_t tempXc = trk.getTrackSagX();
1364          Double_t tempYc = trk.getTrackSagY();
1365
1366          Double_t edgeX1 = 0.;
1367          Double_t edgeY1 = 0.;
1368          Double_t edgePhi1 = std::atan((edgeY1 - tempYc)/(edgeX1 - tempXc));
1369          if((edgePhi1 < 0)&&((edgeX1 - tempXc) < 0)) edgePhi1 += TMath::Pi();
1370          else if((edgePhi1 > 0)&&((edgeX1 - tempXc) < 0))  edgePhi1 += TMath::Pi();
1371
1372          Double_t edgeX2 = trk.getHitR(6)*std::cos(trk.getHitPhi(6));
1373          Double_t edgeY2 = trk.getHitR(6)*std::sin(trk.getHitPhi(6));
1374          Double_t edgePhi2 = std::atan((edgeY2 - tempYc)/(edgeX2 - tempXc));
1375          if((edgePhi2 < 0)&&((edgeX2 - tempXc) < 0)) edgePhi2 += TMath::Pi();
1376          else if((edgePhi2 > 0)&&((edgeX2 - tempXc) < 0))  edgePhi2 += TMath::Pi();
1377
1378          Double_t diffPhi1 = std::abs(edgePhi2 - edgePhi1);
1379          Double_t diffPhi2 = std::abs((edgePhi2 - 2*TMath::Pi()) - edgePhi1);
1380          if(diffPhi2 < diffPhi1) edgePhi2 -= 2*TMath::Pi();
1381
1382          Double_t tempMinPhi = 0.;
1383          Double_t tempMaxPhi = 360.;
1384          if(edgePhi1 < edgePhi2){
1385             tempMinPhi = edgePhi1 - 10.;
1386             tempMaxPhi = edgePhi2 + 10.;
1387          }else{
1388             tempMinPhi = edgePhi2 - 10.;
1389             tempMaxPhi = edgePhi1 + 10.;
1390          }
1391
1392          auto arc = new TArc(tempXc, tempYc, tempR, tempMinPhi, tempMaxPhi);
1393          arc->SetFillStyle(0);
1394
1395          trkLines.push_back(arc);
1396       }
1397
1398       hHitMapED->SetFillColor(3);
1399       hHitMapED->Draw("BOX");
1400       for(Int_t iTrk = 0; iTrk < trkLines.size(); iTrk++){
1401          trkLines.at(iTrk)->Draw("only");
1402       }
1403
1404       c->SaveAs("hQaSagittaFit.root");
1405       c->Close();
1406       // oFile->cd();
1407       // hHitMap->Write();
1408       // fCircleQA->Write();
1409       qaEvent = 9999;
1410    }
1411 }
1412
1413 void InttSeedTrackPerformance::ShowEventInfo(){
1414    std::cout << "TH2D* hHitMapMvtx = new TH2D( \"hHitMapMvtx\", \"hHitMap\", 2000, -100, 100, 2000, -100, 100);" << std::endl;
1415    std::cout << "TH2D* hHitMapIntt = new TH2D( \"hHitMapIntt\", \"hHitMap\", 2000, -100, 100, 2000, -100, 100);" << std::endl;
1416    std::cout << "TH2D* hHitMapTpc = new TH2D( \"hHitMapTpc\", \"hHitMap\", 2000, -100, 100, 2000, -100, 100);" << std::endl;
1417    std::cout << "TH2D* hHitMapEmcal = new TH2D( \"hHitMapEmcal\", \"hHitMap\", 2000, -100, 100, 2000, -100, 100);" << std::endl;
1418    std::cout << "hHitMapMvtx->SetLineWidth(8);" << std::endl;
1419    std::cout << "hHitMapIntt->SetLineWidth(8);" << std::endl;
1420    std::cout << "hHitMapTpc->SetLineWidth(8);" << std::endl;
1421    std::cout << "hHitMapEmcal->SetLineWidth(8);" << std::endl;
1422    std::cout << "hHitMapMvtx->SetLineColor(860);" << std::endl;
1423    std::cout << "hHitMapIntt->SetLineColor(875);" << std::endl;
1424    std::cout << "hHitMapTpc->SetLineColor(870);" << std::endl;
1425    std::cout << "hHitMapEmcal->SetLineColor(807);" << std::endl;
1426
1427
1428    for(Int_t i=0; i<m_fMvtxHits.size(); i++){
1429       std::cout << "hHitMapMvtx->Fill(" << m_fMvtxHits.at(i).r*std::cos(m_fMvtxHits.at(i).phi) \
1430          << ", " << m_fMvtxHits.at(i).r*std::sin(m_fMvtxHits.at(i).phi) << ");" << std::endl;
1431    }
1432    for(Int_t i=0; i<m_sMvtxHits.size(); i++){
1433       std::cout << "hHitMapMvtx->Fill(" << m_sMvtxHits.at(i).r*std::cos(m_sMvtxHits.at(i).phi) \
1434          << ", " << m_sMvtxHits.at(i).r*std::sin(m_sMvtxHits.at(i).phi) << ");" << std::endl;
1435    }
1436    for(Int_t i=0; i<m_tMvtxHits.size(); i++){
1437       std::cout << "hHitMapMvtx->Fill(" << m_tMvtxHits.at(i).r*std::cos(m_tMvtxHits.at(i).phi) \
1438          << ", " << m_tMvtxHits.at(i).r*std::sin(m_tMvtxHits.at(i).phi) << ");" << std::endl;
1439    }
1440    for(Int_t i=0; i< m_iInttHits.size(); i++){
1441       std::cout << "hHitMapIntt->Fill(" << m_iInttHits.at(i).r*std::cos(m_iInttHits.at(i).phi) \
1442          << ", " << m_iInttHits.at(i).r*std::sin(m_iInttHits.at(i).phi) << ");" << std::endl;
1443    }
1444    for(Int_t i=0; i< m_oInttHits.size(); i++){
1445       std::cout << "hHitMapIntt->Fill(" << m_oInttHits.at(i).r*std::cos(m_oInttHits.at(i).phi) \
1446          << ", " << m_oInttHits.at(i).r*std::sin(m_oInttHits.at(i).phi) << ");" << std::endl;
1447    }
1448    for(Int_t i=0; i< m_TpcHits.size(); i++){
1449       std::cout << "hHitMapTpc->Fill(" << m_TpcHits.at(i).r*std::cos(m_TpcHits.at(i).phi) \
1450          << ", " << m_TpcHits.at(i).r*std::sin(m_TpcHits.at(i).phi) << ");" << std::endl;
1451    }
1452    for(Int_t i=0; i< m_emcalHits.size(); i++){
1453       if(m_emcalHits.at(i).energy > 0.1) std::cout << "hHitMapEmcal->Fill(" << m_emcalHits.at(i).r*std::cos(m_emcalHits.at(i).phi) \
1454          << ", " << m_emcalHits.at(i).r*std::sin(m_emcalHits.at(i).phi) << ");" << std::endl;
1455    }
1456
1457    std::cout << "hHitMapMvtx->Draw(\"BOX\");" << std::endl;
1458    std::cout << "hHitMapIntt->Draw(\"BOX same\");" << std::endl;
1459    std::cout << "hHitMapTpc->Draw(\"BOX same\");" << std::endl;
1460    std::cout << "hHitMapEmcal->Draw(\"BOX same\");" << std::endl;
1461 }
1462
1463
1464 void InttSeedTrackPerformance::ShowTrackInfo(tracKuma trk, Double_t dPt, Double_t centerX, Double_t centerY, Double_t sagittaR){
1465    std::cout << "222: dPt, R, Xc, Yc = " << dPt << ", " << sagittaR << ", " << centerX << ", " << centerY << std::endl;
1466    std::cout << "auto arc = new  TArc(" << centerX << ", " << centerY << ", " << sagittaR << ");"  << std::endl;
1467    std::cout << "arc->SetFillStyle(0);" << std::endl;
1468    std::cout << "arc->Draw(\"only\");" << std::endl;
1469
1470
1471    std::cout << "TH2D* hHitMapMvtxTrk = new TH2D( \"hHitMapMvtxTrk\", \"hHitMap\", 2000, -100, 100, 2000, -100, 100);" << std::endl;
1472    std::cout << "TH2D* hHitMapInttTrk = new TH2D( \"hHitMapInttTrk\", \"hHitMap\", 2000, -100, 100, 2000, -100, 100);" << std::endl;
1473    std::cout << "TH2D* hHitMapTpcTrk = new TH2D( \"hHitMapTpcTrk\", \"hHitMap\", 2000, -100, 100, 2000, -100, 100);" << std::endl;
1474    std::cout << "TH2D* hHitMapEmcalTrk = new TH2D( \"hHitMapEmcalTrk\", \"hHitMap\", 2000, -100, 100, 2000, -100, 100);" << std::endl;
1475    std::cout << "hHitMapMvtxTrk->SetLineWidth(3);" << std::endl;
1476    std::cout << "hHitMapInttTrk->SetLineWidth(3);" << std::endl;
1477    std::cout << "hHitMapTpcTrk->SetLineWidth(3);" << std::endl;
1478    std::cout << "hHitMapEmcalTrk->SetLineWidth(3);" << std::endl;
1479    std::cout << "hHitMapMvtxTrk->SetLineColor(632);" << std::endl;
1480    std::cout << "hHitMapInttTrk->SetLineColor(632);" << std::endl;
1481    std::cout << "hHitMapTpcTrk->SetLineColor(632);" << std::endl;
1482    std::cout << "hHitMapEmcalTrk->SetLineColor(632);" << std::endl;
1483
1484    std::cout << "hHitMapMvtxTrk->Fill(" << trk.getHitR(1)*std::cos(trk.getHitPhi(1)) \
1485             << ", " << trk.getHitR(1)*std::sin(trk.getHitPhi(1)) << ");" << std::endl;
1486    std::cout << "hHitMapMvtxTrk->Fill(" << trk.getHitR(2)*std::cos(trk.getHitPhi(2)) \
1487             << ", " << trk.getHitR(2)*std::sin(trk.getHitPhi(2)) << ");" << std::endl;
1488    std::cout << "hHitMapMvtxTrk->Fill(" << trk.getHitR(3)*std::cos(trk.getHitPhi(3)) \
1489             << ", " << trk.getHitR(3)*std::sin(trk.getHitPhi(3)) << ");" << std::endl;
1490    std::cout << "hHitMapInttTrk->Fill(" << trk.getHitR(4)*std::cos(trk.getHitPhi(4)) \
1491             << ", " << trk.getHitR(4)*std::sin(trk.getHitPhi(4)) << ");" << std::endl;
1492    std::cout << "hHitMapInttTrk->Fill(" << trk.getHitR(5)*std::cos(trk.getHitPhi(5)) \
1493             << ", " << trk.getHitR(5)*std::sin(trk.getHitPhi(5)) << ");" << std::endl;
1494    std::cout << "hHitMapEmcalTrk->Fill(" << trk.getHitR(6)*std::cos(trk.getHitPhi(6)) \
1495             << ", " << trk.getHitR(6)*std::sin(trk.getHitPhi(6)) << ");" << std::endl;
1496
1497    std::cout << "hHitMapMvtxTrk->Draw(\"BOX same\");" << std::endl;
1498    std::cout << "hHitMapInttTrk->Draw(\"BOX same\");" << std::endl;
1499    std::cout << "hHitMapTpcTrk->Draw(\"BOX same\");" << std::endl;
1500    std::cout << "hHitMapEmcalTrk->Draw(\"BOX same\");" << std::endl;
1501 }
1502
1503 void InttSeedTrackPerformance::EventJudge(Int_t eventNum, Double_t targetVal, Double_t minLim, Double_t maxLim, bool bIn){
1504    if(bIn){
1505       if((targetVal > minLim)&&(targetVal < maxLim)){
1506          m_vTargetEvents.push_back(eventNum);
1507       }
1508    }else{
1509       if((targetVal < minLim)||(targetVal > maxLim)){
1510          m_vTargetEvents.push_back(eventNum);
1511       }
1512    }
1513
1514 }
1515
1516
1517 void InttSeedTrackPerformance::ChecKuma(std::string checkNo){
1518     std::cout <<"CheeeeeeeeeeeeeeeeeeeeeeeeeeeeeeecKuma"<< checkNo << std::endl;
1519 }
1520
1521
1522 #endif // #ifdef InttSeedTrackPerformance_cxx
1523
1524
1525
1526
1527