Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2026-07-16 08:12:45 

0001 #include "INTT_Calo_trkReco.h"
0002 #include "INTT_Calo_trkUser.h"
0003 
0004 void INTT_Calo_trkUser_Init(void) {
0005   hUSR_ntrk = new TH1F("hUSR_ntrk","number of Calo-INTT tracks",21,-0.5,20.5);
0006   hUSR_EmcIntt_pt  = new TH1F("hUSR_EmcIntt_pt","pt of Calo-INTT track",100,0.,10.);
0007   hUSR_EmcIntt_pt2 = new TH1F("hUSR_EmcIntt_pt2","pt of Calo-INTT track w Mvtx hit",100,0.,10.);
0008   hUSR_EmcIntt_pt3 = new TH1F("hUSR_EmcIntt_pt3","pt of Calo-INTT track w 3 Mvtx hit",100,0.,10.);
0009   hUSR_Intt_pt  = new TH1F("hUSR_Intt_pt","pt of INTT track",100,0.,10.);
0010   hUSR_Intt_pt2 = new TH1F("hUSR_Intt_pt2","pt of INTT track w Mvtx hit",100,0.,10.);
0011   hUSR_Intt_pt3 = new TH1F("hUSR_Intt_pt3","pt of INTT track w 3 Mvtx hit",100,0.,10.);
0012   hUSR_Rpt  = new TH1F("hUSR_Rpt","pt_Intt_mvtx/pt_EmcIntt",100,0.0,4.0);
0013   hUSR_Rpt2 = new TH1F("hUSR_Rpt2","pt_Intt_mvtx/pt_EmcIntt",100,0.0,4.0);
0014   hUSR_nEmc = new TH1F("hUSR_nEmc","# of Emc clusters",101,-0.5,100.5);
0015   hUSR_ntower = new TH1F("hUSR_ntower","# of towers in EMC cluster", 41.,-0.5,40.5);
0016   hUSR_tower_e = new TH1F("hUSR_tower_e","tower energy",100,0.0,5.0);
0017 
0018   hUSR_Intt_t = new TH1F("hUSR_Intt_t","Intt hit timing",11,-5.5,5.5);
0019   hUSR_Mvtx0_t = new TH1F("hUSR_Mvtx0_t","Mvtx0 hit timing",11,-5.5,5.5);
0020   hUSR_Mvtx1_t = new TH1F("hUSR_Mvtx1_t","Mvtx1 hit timing",11,-5.5,5.5);
0021   hUSR_Mvtx2_t = new TH1F("hUSR_Mvtx2_t","Mvtx2 hit timing",11,-5.5,5.5);
0022 
0023   hUSR_EIM_dr = new TH1F("hUSR_EIM_dr","rEmc - rIntt_mvtx",200,-50.,50);
0024   hUSR_EIM_dr2 = new TH1F("hUSR_EIM_dr2","rEmc - rIntt_mvtx",200,-50.,50);
0025   hUSR_EIM_dz = new TH1F("hUSR_EIM_dz","zEmc - zproj",200,-50.,50.);
0026   hUSR_EIM_dz2 = new TH1F("hUSR_EIM_dz2","zEmc - zproj",200,-50.,50.);
0027   hUSR_EIM_dz3 = new TH1F("hUSR_EIM_dz3","zEmc - zproj",200,-50.,50.);
0028 
0029   hUSR_Mee = new TH1F("hUSR_Mee","Mee",1000,0.0,4.0);
0030   hUSR_pe = new TH1F("hUSR_pe","pe(tagged by Mee<0.04)",100,0.,5.);
0031   hUSR_Ep = new TH1F("hUSR_Ep","E/p for conversion e+/e-",100,0.,2.0);
0032 
0033   hUSR_truth_pt = new TH1F("hUSR_truth_pt","pt (truth)",100,0.,5.0);
0034 
0035   hUSR_pt0m_tr_pt = new TH1F("hUSR_pt0m_tr_pt","pt_mvtx/pt_truth",100,0.0,2.0);
0036   hUSR_pt0m_tr_pt1 = new TH1F("hUSR_pt0m_tr_pt1","pt_mvtx/pt_truth (0.3<pt<0.5GeV/c)",100,0.0,2.0);
0037   hUSR_pt0m_tr_pt2 = new TH1F("hUSR_pt0m_tr_pt2","pt_mvtx/pt_truth (0.5<pt<0.7GeV/c)",100,0.0,2.0);
0038   hUSR_pt0m_tr_pt3 = new TH1F("hUSR_pt0m_tr_pt3","pt_mvtx/pt_truth (0.7<pt<1.0GeV/c)",100,0.0,2.0);
0039   hUSR_pt0m_tr_pt4 = new TH1F("hUSR_pt0m_tr_pt4","pt_mvtx/pt_truth (1.0<pt<1.5GeV/c)",100,0.0,2.0);
0040   hUSR_pt0m_tr_pt5 = new TH1F("hUSR_pt0m_tr_pt5","pt_mvtx/pt_truth (1.5<pt<2.0GeV/c)",100,0.0,2.0);
0041   hUSR_pt0m_tr_pt6 = new TH1F("hUSR_pt0m_tr_pt6","pt_mvtx/pt_truth (2.0<pt<3.0GeV/c)",100,0.0,2.0);
0042 
0043   hUSR_pt0e_tr_pt = new TH1F("hUSR_pt0e_tr_pt","pt_emc/pt_truth",100,0.0,2.0);
0044   hUSR_pt0e_tr_pt1 = new TH1F("hUSR_pt0e_tr_pt1","pt_emc/pt_truth (0.3<pt<0.5GeV/c)",100,0.0,2.0);
0045   hUSR_pt0e_tr_pt2 = new TH1F("hUSR_pt0e_tr_pt2","pt_emc/pt_truth (0.5<pt<0.7GeV/c)",100,0.0,2.0);
0046   hUSR_pt0e_tr_pt3 = new TH1F("hUSR_pt0e_tr_pt3","pt_emc/pt_truth (0.7<pt<1.0GeV/c)",100,0.0,2.0);
0047   hUSR_pt0e_tr_pt4 = new TH1F("hUSR_pt0e_tr_pt4","pt_emc/pt_truth (1.0<pt<1.5GeV/c)",100,0.0,2.0);
0048   hUSR_pt0e_tr_pt5 = new TH1F("hUSR_pt0e_tr_pt5","pt_emc/pt_truth (1.5<pt<2.0GeV/c)",100,0.0,2.0);
0049   hUSR_pt0e_tr_pt6 = new TH1F("hUSR_pt0e_tr_pt6","pt_emc/pt_truth (2.0<pt<3.0GeV/c)",100,0.0,2.0);
0050   
0051   nt_trk = new TNtuple("nt_trk","INTT track Ntuple","p_emc:p_mvtx:tot_e:emc_e:sign:ntrk:pt_mvtx");
0052   nt_sim_trk = new TNtuple("nt_sim_trk","Sim and Reco Ntuple","tr_pt:tr_phi0:tr_eta0:pt0e:phi_e:pt0m");
0053 }
0054 
0055 void INTT_Calo_trkUser_End(void) {
0056   hUSR_ntrk->Write();
0057   hUSR_Intt_pt->Write();
0058   hUSR_Intt_pt2->Write();
0059   hUSR_Intt_pt3->Write();
0060   hUSR_EmcIntt_pt->Write();
0061   hUSR_EmcIntt_pt2->Write();
0062   hUSR_EmcIntt_pt3->Write();
0063   hUSR_Rpt->Write();
0064   hUSR_Rpt2->Write();
0065   hUSR_nEmc->Write();
0066   hUSR_ntower->Write();
0067   hUSR_tower_e->Write();
0068 
0069   hUSR_Intt_t->Write();
0070   hUSR_Mvtx0_t->Write();
0071   hUSR_Mvtx1_t->Write();
0072   hUSR_Mvtx2_t->Write();
0073 
0074   hUSR_EIM_dr->Write();
0075   hUSR_EIM_dr2->Write();
0076   hUSR_EIM_dz->Write();
0077   hUSR_EIM_dz2->Write();
0078   hUSR_EIM_dz3->Write();
0079 
0080   hUSR_Mee->Write();
0081   hUSR_pe->Write();
0082   hUSR_Ep->Write();
0083 
0084   hUSR_truth_pt->Write();
0085 
0086   hUSR_pt0m_tr_pt->Write();
0087   hUSR_pt0m_tr_pt1->Write();
0088   hUSR_pt0m_tr_pt2->Write();
0089   hUSR_pt0m_tr_pt3->Write();
0090   hUSR_pt0m_tr_pt4->Write();
0091   hUSR_pt0m_tr_pt5->Write();
0092   hUSR_pt0m_tr_pt6->Write();
0093 
0094   hUSR_pt0e_tr_pt->Write();
0095   hUSR_pt0e_tr_pt1->Write();
0096   hUSR_pt0e_tr_pt2->Write();
0097   hUSR_pt0e_tr_pt3->Write();
0098   hUSR_pt0e_tr_pt4->Write();
0099   hUSR_pt0e_tr_pt5->Write();
0100   hUSR_pt0e_tr_pt6->Write();
0101   
0102   nt_trk->Write();
0103   nt_sim_trk->Write();
0104 }
0105 
0106 void INTT_Calo_trkUser(void) {
0107   int ntrk = vCaloInttMvtx.size();
0108   if(imode == 2) { cout << "ntrk="<<ntrk<<endl;}
0109   hUSR_ntrk->Fill(ntrk);
0110 
0111   int nEmc = emc_r->size();
0112   hUSR_nEmc->Fill(nEmc);
0113   if(nEmc>0) {
0114     for(int iemc=0;iemc<nEmc;iemc++) {
0115       vector<float> vtower_e= emc_tower_e->at(iemc);
0116       int ntower = vtower_e.size();
0117       hUSR_ntower->Fill(ntower);
0118       for(int itower=0;itower<ntower;itower++) {
0119     hUSR_tower_e->Fill(vtower_e.at(itower));
0120       }
0121     }
0122   }
0123   
0124   // positive tracks
0125   vector<float> vpxp;
0126   vector<float> vpyp;
0127   vector<float> vpzp;
0128   vector<float> vEp;  //Emc energy  
0129 
0130   // negative tracks
0131   vector<float> vpxm;
0132   vector<float> vpym;
0133   vector<float> vpzm;
0134   vector<float> vEm;  //Emc energy
0135   
0136   if(ntrk > 0) {
0137     /*
0138     cout << "i_event="<<i_event;
0139     cout << " ntrk="<<ntrk<<endl;
0140     */
0141 
0142     vpxp.clear();
0143     vpyp.clear();
0144     vpzp.clear();
0145     
0146     vpxm.clear();
0147     vpym.clear();
0148     vpzm.clear();
0149     
0150     for(int i=0;i<ntrk;i++) {
0151       CaloInttMvtx CIMtrack = vCaloInttMvtx.at(i);
0152       float pt_emc_intt  = CIMtrack.pt0e;
0153       float pt_intt_mvtx = CIMtrack.pt0m;
0154       float phi0m = CIMtrack.phi0m;
0155       float emc_e = CIMtrack.emc_e;
0156       float emcal_e = CIMtrack.emcal_e;
0157       float total_e = CIMtrack.total_e;
0158       float phi0 = CIMtrack.phi0e;
0159       float z0 = CIMtrack.z0e;
0160       int Intt0_t = vINTT0t.at(CIMtrack.iINTT0);
0161       int Intt1_t = vINTT1t.at(CIMtrack.iINTT1);
0162       hUSR_Intt_pt->Fill(pt_intt_mvtx);
0163       if(pt_emc_intt>0) hUSR_EmcIntt_pt->Fill(pt_emc_intt);
0164       hUSR_Intt_t->Fill(Intt0_t);
0165       hUSR_Intt_t->Fill(Intt1_t);
0166       
0167       int iMvtx0 = CIMtrack.iMvtx0;
0168       int iMvtx1 = CIMtrack.iMvtx1;
0169       int iMvtx2 = CIMtrack.iMvtx2;
0170       
0171       int nMvtxHits = 0;
0172       if( iMvtx2 >=0) nMvtxHits++;
0173       if( iMvtx1 >=0) nMvtxHits++;
0174       if( iMvtx0 >=0) nMvtxHits++;
0175 
0176       if(imode == 2) {
0177     cout << i << ": INTT0="<<CIMtrack.iINTT0<<" INTT1="<<CIMtrack.iINTT1;
0178     cout <<" Mvtx0="<<iMvtx0<<" Mvtx1="<<iMvtx1<<" iMvtx2= "<<iMvtx2<<endl;
0179       }
0180       if( nMvtxHits > 0) {
0181     hUSR_Intt_pt2->Fill(pt_intt_mvtx);
0182     if(pt_emc_intt>0) {
0183       hUSR_EmcIntt_pt2->Fill(pt_emc_intt);
0184       hUSR_Rpt->Fill(pt_intt_mvtx/pt_emc_intt);
0185       if(0.8 < pt_emc_intt && pt_emc_intt < 4.0) hUSR_Rpt2->Fill(pt_intt_mvtx/pt_emc_intt);
0186     }
0187       }//if( nMvtxHits)
0188 
0189       if( nMvtxHits >=3) {
0190     hUSR_Intt_pt3->Fill(pt_intt_mvtx);
0191     if(pt_emc_intt>0) hUSR_EmcIntt_pt3->Fill(pt_emc_intt);
0192       }//if(nMvtxHits>=3)
0193 
0194       if(nMvtxHits>=2) {
0195     // at least two Mvt layers has a hit for this track
0196     // thus an Intt-Mvtx track is formed and its track parameters
0197     // at the origin (x0m,y0m,z0m,pt0m,phi0m,pz0m) are well defined
0198     //
0199     int ntop = vTop_r.size();
0200     if(ntop>0) {
0201       for(int itop=0;itop<ntop;itop++) {
0202         if(vTop_emc_emc_e.at(itop) > 0) {// Emc Hit in the topocluster
0203           float emc_r   = vTop_emc_r.at(itop);
0204           float emc_phi = vTop_emc_phi.at(itop);
0205           float zemc = vTop_emc_z.at(itop);
0206           float xemc = emc_r*cos(emc_phi);
0207           float yemc = emc_r*sin(emc_phi);
0208           float x0m = CIMtrack.x0m;
0209           float y0m = CIMtrack.y0m;
0210           float z0m = CIMtrack.z0m;
0211           float pt0m = CIMtrack.pt0m;
0212           float phi0m = CIMtrack.phi0m;
0213           float pz0m  = CIMtrack.pz0m;
0214           float rc = CIMtrack.r_intt_mvtx;
0215           float xc = CIMtrack.xc_intt_mvtx;
0216           float yc = CIMtrack.yc_intt_mvtx;
0217           float rcEmc = sqrt((xemc-xc)*(xemc-xc)+(yemc-yc)*(yemc-yc));
0218           float dzdr = pz0m/pt0m;
0219           float zproj = z0m + dzdr*emc_r;
0220           hUSR_EIM_dr->Fill(rcEmc - rc);
0221           hUSR_EIM_dz->Fill(zemc - zproj);
0222           if(abs(rcEmc-rc)<20) hUSR_EIM_dz2->Fill(zemc - zproj);
0223           if(abs(zemc-zproj)<8) hUSR_EIM_dr2->Fill(rcEmc - rc);
0224           if(abs(rcEmc-rc)<20 && abs(zemc-zproj)<8) {
0225         hUSR_EIM_dz3->Fill(zemc - zproj);
0226           }
0227         }//if
0228       }//for
0229     }//if(ntop)
0230 
0231     if(pt_emc_intt>0 && pt_intt_mvtx > 0) {
0232       float pz0e = CIMtrack.pz0e;
0233       float pz0m = CIMtrack.pz0m;
0234       float px0m = pt_intt_mvtx*cos(phi0m);
0235       float py0m = pt_intt_mvtx*sin(phi0m);
0236       float p_emc  = sqrt(pt_emc_intt*pt_emc_intt+pz0e*pz0e);
0237       float p_mvtx = sqrt(pt_intt_mvtx*pt_intt_mvtx+pz0m*pz0m);
0238       int sign = CIMtrack.sign;
0239       nt_trk->Fill(p_emc,p_mvtx,total_e,emcal_e,sign,ntrk,pt_intt_mvtx);
0240     }//if(pt_emc_intt...)   
0241       }//if(nMvtxHits>=2)
0242 
0243       if(pt_intt_mvtx>0.2 && nMvtxHits>=2) {
0244     float pz0m = CIMtrack.pz0m;
0245     float px0m = pt_intt_mvtx*cos(phi0m);
0246     float py0m = pt_intt_mvtx*sin(phi0m);
0247     int sign = CIMtrack.sign;
0248     if(sign > 0) {
0249       vpxp.push_back(px0m);
0250       vpyp.push_back(py0m);
0251       vpzp.push_back(pz0m);
0252       vEp.push_back(emc_e);
0253     } else {
0254       vpxm.push_back(px0m);
0255       vpym.push_back(py0m);
0256       vpzm.push_back(pz0m);
0257       vEm.push_back(emc_e);
0258     }
0259       }//if(pt_intt_mvtx)
0260       
0261       if( iMvtx0 >= 0) hUSR_Mvtx0_t->Fill(vMvtx0t.at(iMvtx0));
0262       if( iMvtx1 >= 0) hUSR_Mvtx1_t->Fill(vMvtx1t.at(iMvtx1));
0263       if( iMvtx2 >= 0) hUSR_Mvtx2_t->Fill(vMvtx2t.at(iMvtx2));
0264     }// for(i)
0265 
0266     const float me = 0.000511;
0267     int Np = vpxp.size();
0268     int Nm = vpxm.size();
0269     if(Np>0 && Nm>0) {
0270       for(int ip=0;ip<Np;ip++) {
0271     float pxp=vpxp.at(ip);
0272     float pyp=vpyp.at(ip);
0273     float pzp=vpzp.at(ip);
0274     float pp = sqrt(pxp*pxp + pyp*pyp + pzp*pzp);
0275     float Ep = sqrt(pp*pp + me*me);
0276     for(int im=0;im<Nm;im++) {
0277       float pxm=vpxm.at(im);
0278       float pym=vpym.at(im);
0279       float pzm=vpzm.at(im);
0280       float pm = sqrt(pxm*pxm + pym*pym + pzm*pzm);
0281       float Em = sqrt(pm*pm + me*me);
0282 
0283       float Eee  = Ep+Em;
0284       float pxee = pxp+pxm;
0285       float pyee = pyp+pym;
0286       float pzee = pzp+pzm;
0287       float Mee = sqrt(Eee*Eee - pxee*pxee - pyee*pyee - pzee*pzee);
0288       hUSR_Mee->Fill(Mee);
0289 
0290       if( Mee < 0.04) {
0291         if(vEp.at(ip)>0) {
0292           hUSR_pe->Fill(pp);
0293           hUSR_Ep->Fill(vEp.at(ip)/pp);
0294         }
0295         if(vEm.at(im)>0) {
0296           hUSR_pe->Fill(pm);
0297           hUSR_Ep->Fill(vEm.at(im)/pm);
0298         }
0299       }
0300     }//for(im)
0301       }//for(ip)
0302     }
0303   }//if(ntrk)
0304 
0305   // Simulation 
0306   if(truthTree != nullptr) {
0307     int n_tr_trk = truth_pid->size();  //# of truth tracks
0308     //    cout <<"n_tr_trk="<<n_tr_trk<<endl;
0309     if(n_tr_trk>0) {
0310       for(int itr = 0; itr<n_tr_trk; itr++) {
0311     float tr_px=truth_px->at(itr);
0312     float tr_py=truth_py->at(itr);
0313     float tr_pz=truth_pz->at(itr);
0314     float tr_pt=sqrt(tr_px*tr_px + tr_py*tr_py);
0315     float tr_phi0 = atan2(tr_py,tr_px);
0316     float tr_eta0 = rz2eta(tr_pt,tr_pz);
0317     hUSR_truth_pt->Fill(tr_pt);
0318 
0319     if( n_tr_trk==1 && ntrk==1) {
0320       if(itr==0) {
0321         CaloInttMvtx CIMtrack = vCaloInttMvtx.at(0);
0322         float pt0m = CIMtrack.pt0m;
0323         float pt0e = CIMtrack.pt0e;
0324         float xemc = CIMtrack.xemc;
0325         float yemc = CIMtrack.yemc;
0326         float phi_e = atan2(yemc,xemc);
0327         
0328         nt_sim_trk->Fill(tr_pt,tr_phi0,tr_eta0,pt0e,phi_e,pt0m);
0329         if(tr_pt > 0.3) {
0330           hUSR_pt0m_tr_pt->Fill(pt0m/tr_pt);
0331           hUSR_pt0e_tr_pt->Fill(pt0e/tr_pt);
0332           if(tr_pt <0.5) {
0333         hUSR_pt0m_tr_pt1->Fill(pt0m/tr_pt);
0334         hUSR_pt0e_tr_pt1->Fill(pt0e/tr_pt);
0335           } else if(tr_pt<0.7) {
0336         hUSR_pt0m_tr_pt2->Fill(pt0m/tr_pt);
0337         hUSR_pt0e_tr_pt2->Fill(pt0e/tr_pt);
0338           } else if(tr_pt<1.0) {
0339         hUSR_pt0m_tr_pt3->Fill(pt0m/tr_pt);
0340         hUSR_pt0e_tr_pt3->Fill(pt0e/tr_pt);
0341           } else if(tr_pt<1.5) {
0342         hUSR_pt0m_tr_pt4->Fill(pt0m/tr_pt);
0343         hUSR_pt0e_tr_pt4->Fill(pt0e/tr_pt);
0344           } else if(tr_pt<2.0) {
0345         hUSR_pt0m_tr_pt5->Fill(pt0m/tr_pt);
0346         hUSR_pt0e_tr_pt5->Fill(pt0e/tr_pt);
0347           } else if(tr_pt<3.0) {
0348         hUSR_pt0m_tr_pt6->Fill(pt0m/tr_pt);
0349         hUSR_pt0e_tr_pt6->Fill(pt0e/tr_pt);
0350           }
0351         }
0352       }//if(itr)
0353     }//if(n_tr_trk,ntrk)
0354       }//for(itr)
0355     }//if(n_tr_trk)
0356   }//if(truthTree)
0357 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!