Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-12-16 09:18:01 

0001 #include <TFile.h>
0002 #include <TTree.h>
0003 #include <TArc.h>
0004 #include <TLine.h>
0005 #include <TGraph.h>
0006 #include <TGraph2D.h>
0007 #include <TH3.h>
0008 #include <TCanvas.h>
0009 #include <TVirtualPad.h>
0010 #include <TMath.h>
0011 #include <TStyle.h>
0012 
0013 #include <iostream>
0014 #include <vector>
0015 
0016 
0017 TH3F     *frame3d = nullptr;
0018 TGraph2D *gr2_xyz = nullptr;
0019 
0020 void calcTrajectoryXY(
0021         float trk_x, float trk_y, float trk_pt, float trk_phi, int trk_c,
0022         float& R, float& cx, float& cy, float& arcphi0, float& arcphi1)
0023 {
0024 
0025   const float B = 1.4; // Tesla
0026   const float C = 0.299792458; // m/ns
0027   /////////////////////////////
0028 
0029   float px = trk_pt*std::cos(trk_phi);
0030   float py = trk_pt*std::sin(trk_phi);
0031 
0032 
0033   R = trk_pt/(C*B) * 100.; // p=0.3BR in B field,  0.3 from C (speed of light)
0034                            // 100 to convert m to cm
0035 
0036   cx = R * ( py/trk_pt)*trk_c + trk_x;
0037   cy = R * (-px/trk_pt)*trk_c + trk_y;
0038 
0039   // phi angle calculation
0040   float rx = trk_x - cx; // circle origin to trackpoint
0041   float ry = trk_y - cy;
0042   float rv_x_pv = rx*py - ry*px; // direction of track relative to circle vector
0043 
0044   const float pi = TMath::Pi();
0045 
0046   //if(R>75.) 
0047   { 
0048     float arcphi = atan2(ry, rx);  // -pi ~ pi
0049 
0050     if(rv_x_pv<0.) { arcphi0 = arcphi-pi; arcphi1 = arcphi;}
0051     else           { arcphi0 = arcphi;    arcphi1 = arcphi+pi;}
0052   }
0053 }
0054 
0055 void drawDisplay(int eventnumber=10)
0056 {
0057  //std::string filepath = "/sphenix/user/jaein213/tracking/SiliconSeeding/MC/macro/ana/jpsi/ana_Allpart.root";
0058  //std::string filepath = "/sphenix/user/jaein213/tracking/SiliconSeeding/MC/macro/ana/jpsi/ana_all.root";
0059  //std::string filepath = "/sphenix/tg/tg01/commissioning/INTT/work/mahiro/SIliconCalo/MC/ana_e-/merged_10k.root";
0060  //std::string filepath = "/sphenix/tg/tg01/commissioning/INTT/work/mahiro/SIliconCalo/MC/ana/ana_0_1kevt.root";
0061  //std::string filepath = "/sphenix/tg/tg01/commissioning/INTT/work/mahiro/SIliconCalo/MC/ana/merged_100k.root";
0062  //std::string filepath = "/sphenix/tg/tg01/commissioning/INTT/work/hachiya/SiCaloTrack/data/ana_pythia_old/ana_addall.root";
0063  //std::string filepath = "/sphenix/tg/tg01/commissioning/INTT/work/hachiya/SiCaloTrack/data/ana/all_pythia.root";
0064  //std::string filepath = "/sphenix/tg/tg01/commissioning/INTT/work/hachiya/SiCaloTrack/data/ana/all_pythia_old.root";
0065  //std::string filepath = "../macro/ana_0.root";
0066  //std::string filepath = "/sphenix/user/hachiya/INTT/INTT/general_codes/hachiya/SiliconSeeding/macro/data/test_ana.root";
0067  std::string filepath = "/sphenix/tg/tg01/commissioning/INTT/work/hachiya/SiCaloTrack/data/ana_em/all_em.root";
0068 
0069 /*
0070   // Load ROOT libraries if not already loaded
0071   gSystem->Load("libTree");
0072   gSystem->Load("libGraf");
0073   gSystem->Load("libHist");
0074 */
0075 
0076   // Open the ROOT file
0077   TFile* file = TFile::Open(filepath.c_str(), "READ");
0078 
0079   // Get the trackTree and caloTree
0080   TTree* trackTree = (TTree*)file->Get("trackTree");
0081   TTree* caloTree  = (TTree*)file->Get("caloTree");
0082   //TTree* siClusTree= (TTree*)file->Get("SiClusTree");
0083   TTree* siClusTree= (TTree*)file->Get("SiClusAllTree");
0084   TTree* truthTree = (TTree*)file->Get("truthTree");
0085 
0086   if (!trackTree) {
0087     printf("Error: trackTree not found in file.\n");
0088     file->Close();
0089     return;
0090   }
0091   if (!siClusTree) {
0092     printf("Error: SiclusTree not found in file.\n");
0093     file->Close();
0094     return;
0095   }
0096   if (!caloTree) {
0097     printf("Error: caloTree not found in file.\n");
0098     file->Close();
0099     return;
0100   }
0101   if (!truthTree) {
0102     printf("No truthTree ound in file.\n");
0103   }
0104 
0105   // Set up branches for trackcluster_x, trackcluster_y, trackcluster_z, trackcluster_trackid, evt
0106   std::vector<float>* trackcluster_x = nullptr;
0107   std::vector<float>* trackcluster_y = nullptr;
0108   std::vector<float>* trackcluster_z = nullptr;
0109   std::vector<int>*   trackcluster_trackid = nullptr;
0110 
0111   std::vector<float>* track_x = nullptr;
0112   std::vector<float>* track_y = nullptr;
0113   std::vector<float>* track_z = nullptr;
0114   std::vector<float>* track_x_emc = nullptr;
0115   std::vector<float>* track_y_emc = nullptr;
0116   std::vector<float>* track_z_emc = nullptr;
0117   std::vector<float>* track_phi = nullptr;
0118   std::vector<float>* track_eta = nullptr;
0119   std::vector<float>* track_pt = nullptr;
0120   std::vector<int>*   track_charge = nullptr;
0121 
0122   int evt = -1;
0123   siClusTree->SetBranchAddress("Siclus_x", &trackcluster_x);
0124   siClusTree->SetBranchAddress("Siclus_y", &trackcluster_y);
0125   siClusTree->SetBranchAddress("Siclus_z", &trackcluster_z);
0126   siClusTree->SetBranchAddress("Siclus_trackid", &trackcluster_trackid);
0127 
0128   trackTree->SetBranchAddress("evt",    &evt);
0129   trackTree->SetBranchAddress("x0",     &track_x);
0130   trackTree->SetBranchAddress("y0",     &track_y);
0131   trackTree->SetBranchAddress("z0",     &track_z);
0132   trackTree->SetBranchAddress("phi0",   &track_phi);
0133   trackTree->SetBranchAddress("eta0",   &track_eta);
0134   trackTree->SetBranchAddress("pt0",    &track_pt);
0135   trackTree->SetBranchAddress("charge", &track_charge);
0136   // Add branch addresses for EMCAL state and track parameters
0137   trackTree->SetBranchAddress("x_proj_emc", &track_x_emc);
0138   trackTree->SetBranchAddress("y_proj_emc", &track_y_emc);
0139   trackTree->SetBranchAddress("z_proj_emc", &track_z_emc);
0140 
0141   // Set up branches for caloTree (calo_x, calo_y, calo_z, calo_evt)
0142   std::vector<float>* calo_x = nullptr;
0143   std::vector<float>* calo_y = nullptr;
0144   std::vector<float>* calo_z = nullptr;
0145   int calo_evt = -1;
0146   caloTree->SetBranchAddress("calo_evt", &calo_evt);
0147   caloTree->SetBranchAddress("x", &calo_x);
0148   caloTree->SetBranchAddress("y", &calo_y);
0149   caloTree->SetBranchAddress("z", &calo_z);
0150 
0151   // setup for truthTree
0152   std::vector<float>* truth_pt=0, *truth_eta=0, *truth_phi=0;
0153   std::vector<float>* truth_px=0, *truth_py=0, *truth_pz=0, *truth_e=0;
0154   std::vector<int>*  truth_pid=0;
0155   std::vector<int>*  truth_vtxid=0;
0156   std::vector<float>* truth_vtx_x=0, *truth_vtx_y=0, *truth_vtx_z=0; // indix vtxid, not trackid
0157   if(truthTree){
0158     truthTree->SetBranchAddress("truth_pid", &truth_pid);
0159     truthTree->SetBranchAddress("truth_pt", &truth_pt);
0160     truthTree->SetBranchAddress("truth_phi",&truth_phi);
0161     truthTree->SetBranchAddress("truth_px", &truth_px);
0162     truthTree->SetBranchAddress("truth_py", &truth_py);
0163     truthTree->SetBranchAddress("truth_pz", &truth_pz);
0164     truthTree->SetBranchAddress("truth_vtxid", &truth_vtxid);
0165     truthTree->SetBranchAddress("truth_vtx_x", &truth_vtx_x);
0166     truthTree->SetBranchAddress("truth_vtx_y", &truth_vtx_y);
0167     truthTree->SetBranchAddress("truth_vtx_z", &truth_vtx_z);
0168   }
0169  
0170   
0171 
0172   // Find the entry for the given eventnumber in trackTree
0173   Long64_t nentries = trackTree->GetEntries();
0174 
0175   bool found = false;
0176   std::vector<float> xs_track, ys_track;
0177   // For rz projection
0178   //std::vector<float> rz_trackcluster_z, rz_trackcluster_r;
0179   std::vector<float> rz_track_z, rz_track_r;
0180   std::vector<float> rz_track_z_emc, rz_track_r_emc;
0181   // For trajectory drawing
0182   //--std::vector<std::vector<float>> track_traj_xs;
0183   //--std::vector<std::vector<float>> track_traj_ys;
0184   //--std::vector<std::vector<float>> track_traj_zs;
0185   //--std::vector<std::vector<float>> track_traj_rs;
0186 
0187 
0188   TGraph* gr = NULL;
0189   TGraph* gr_rz = NULL;
0190   TGraph* gr_calo = NULL;
0191   TGraph* gr_rz_calo = NULL;
0192   TGraph* gr_phiz_calo = NULL;
0193   TGraph* gr_phiz_pemc = NULL;
0194   std::vector<TArc*> v_trk_xy;
0195   std::vector<TLine*> v_trk_zr;
0196 
0197   std::vector<TArc*> v_tr_xy;
0198   std::vector<TLine*> v_tr_zr;
0199 
0200   // frame constants
0201   const float si_cr[5] = {2.4, 3.2, 3.9, 7.4, 10.};
0202   const float si_cz[5] = {10,10,10,23.5, 23.5};
0203   const float calo_cr[1] = {93.5};
0204   const float calo_cz[1] = {120.};
0205   ///////////////////////////
0206   std::vector<TLine*>    v_line_rz;
0207   // Si RZ-frame
0208   v_line_rz.push_back(new TLine(-si_cz[0],  si_cr[0], si_cz[0],  si_cr[0]));
0209   v_line_rz.push_back(new TLine(-si_cz[1],  si_cr[1], si_cz[1],  si_cr[1]));
0210   v_line_rz.push_back(new TLine(-si_cz[2],  si_cr[2], si_cz[2],  si_cr[2]));
0211   v_line_rz.push_back(new TLine(-si_cz[3],  si_cr[3], si_cz[3],  si_cr[3]));
0212   v_line_rz.push_back(new TLine(-si_cz[4],  si_cr[4], si_cz[4],  si_cr[4]));
0213   v_line_rz.push_back(new TLine(-si_cz[0], -si_cr[0], si_cz[0], -si_cr[0]));
0214   v_line_rz.push_back(new TLine(-si_cz[1], -si_cr[1], si_cz[1], -si_cr[1]));
0215   v_line_rz.push_back(new TLine(-si_cz[2], -si_cr[2], si_cz[2], -si_cr[2]));
0216   v_line_rz.push_back(new TLine(-si_cz[3], -si_cr[3], si_cz[3], -si_cr[3]));
0217   v_line_rz.push_back(new TLine(-si_cz[4], -si_cr[4], si_cz[4], -si_cr[4]));
0218   // Calo RZ-frame
0219   v_line_rz.push_back(new TLine(-calo_cz[0], calo_cr[0], calo_cz[0], calo_cr[0]));
0220   v_line_rz.push_back(new TLine(-calo_cz[0],-calo_cr[0], calo_cz[0],-calo_cr[0]));
0221 
0222 
0223   std::vector<TEllipse*> v_elli_xy;
0224   // Si + calo X-Y frame
0225   v_elli_xy.push_back(new TEllipse(0,0, si_cr[0]));
0226   v_elli_xy.push_back(new TEllipse(0,0, si_cr[1]));
0227   v_elli_xy.push_back(new TEllipse(0,0, si_cr[2]));
0228   v_elli_xy.push_back(new TEllipse(0,0, si_cr[3]));
0229   v_elli_xy.push_back(new TEllipse(0,0, si_cr[4]));
0230   v_elli_xy.push_back(new TEllipse(0,0, calo_cr[0]));
0231 
0232   for (size_t il = 0; il < v_elli_xy.size(); ++il) {
0233     v_elli_xy[il]->SetFillStyle(0);
0234   }
0235 
0236   frame3d = new TH3F("frame3d","frame3d",1,-30, 30,1,-15, 15,1,-18, 18);  
0237 
0238   gStyle->SetOptStat(0);
0239 
0240   char input[256];
0241   TCanvas* can = new TCanvas("can", "Event Display XY", 1400, 700);
0242   can->Divide(2,1);
0243   TVirtualPad *p1 = can->cd(1);
0244   p1->Divide(2,2);
0245 
0246   ////////////////
0247   for (Long64_t i = 0; i < nentries; ++i) {
0248     if(truthTree) {
0249       truthTree->GetEntry(i);
0250       
0251       if(truth_vtx_z->size()>0 ) {
0252         std::cout<<"Truth Zvtx : "<<truth_vtx_z->at(0);
0253         if(fabs(truth_vtx_z->at(0))>5) {
0254           std::cout<<" : skipped"<<std::endl;
0255           continue;
0256         }
0257         std::cout<<" : good vtx"<<std::endl;
0258       }
0259       else {
0260         std::cout<<" : vtx size=0"<<std::endl;
0261       }
0262 
0263       // clear truth track vector
0264       for(size_t iarc=0; iarc<v_tr_xy.size(); iarc++){
0265         delete v_tr_xy[iarc];
0266         delete v_tr_zr[iarc];
0267       }
0268       v_tr_xy.clear();
0269       v_tr_zr.clear();
0270 
0271       bool ev_skip=false;
0272       std::cout<<"truth N: "<<truth_pid->size()<<std::endl;
0273       for (size_t k = 0; k < truth_pid->size(); ++k) {
0274 
0275         int   tr_vid = truth_vtxid->at(k)-1;
0276         float tr_x   = -999;
0277         float tr_y   = -999;
0278         float tr_z   = -999;
0279         if(0<=tr_vid&&tr_vid<truth_vtx_x->size()){
0280           tr_x = truth_vtx_x->at(tr_vid);
0281           tr_y = truth_vtx_y->at(tr_vid);
0282           tr_z = truth_vtx_z->at(tr_vid);
0283         }
0284         float& tr_px = truth_px->at(k);
0285         float& tr_py = truth_py->at(k);
0286         float& tr_pz = truth_pz->at(k);
0287 
0288         float tr_phi = truth_phi->at(k); //atan2(tr_py, tr_px);
0289         float tr_pt  = truth_pt->at(k); //sqrt(tr_px*tr_px + tr_py*tr_py);
0290         //float tr_eta = track_eta->at(k);
0291         int&   tr_pid = truth_pid->at(k);
0292         int    tr_c   = tr_pid<0 ? 1 : -1;
0293         float  tr_the = atan2(tr_pt, tr_pz);
0294 
0295         if(tr_pt>0.5) {
0296           cout<<"pt_excced : "<<tr_pt<<endl;
0297           ev_skip=true;
0298           break;
0299           //continue;
0300         }
0301 
0302         ////////////////////////
0303         // X-Y trajectory
0304         float tR, tcx, tcy, tarcphi0, tarcphi1;
0305         calcTrajectoryXY(tr_x, tr_y, tr_pt, tr_phi, tr_c,
0306                          tR, tcx, tcy, tarcphi0, tarcphi1);
0307 
0308 
0309          
0310         
0311         const float conv_ang = 180./TMath::Pi();
0312         float tarcphi0_ang = tarcphi0*conv_ang; 
0313         float tarcphi1_ang = tarcphi1*conv_ang; 
0314         
0315         //float phi_ang     = trk_phi*conv_ang; if(phi_ang<0) phi_ang += 360.;
0316         std::cout<<"px, py, R = "<<tr_pt<<" "<<tr_phi<<" "<<tR<<" "<<tcx<<" "<<tcy<<",   "
0317                  <<tarcphi0_ang<<" "<<tarcphi1_ang<<" "<<" "<<tr_x<<" "<<tr_y<<std::endl;
0318         std::cout<<" = "<<tr_vid<<" "<<tr_x<<" "<<tr_y<<" "<<tr_z<<std::endl;
0319 
0320 
0321 
0322         TArc *tr_xy = new TArc(tcx, tcy, tR, tarcphi0_ang, tarcphi1_ang); //, ang0*180/TMath::Pi(), ang1*180/TMath::Pi());
0323         tr_xy->SetFillStyle(20);
0324         tr_xy->SetLineStyle(5);
0325         tr_xy->SetLineWidth(1);
0326         tr_xy->SetLineColor(2);
0327         v_tr_xy.push_back(tr_xy);
0328         
0329         ////////////////////////
0330         // R-Z trajectory
0331         float r_end = 100.;
0332 
0333         float tr_r = std::sqrt(tr_x*tr_x+tr_y*tr_y);
0334 
0335         float slope = std::tan(tr_the);
0336         float z_end = (r_end-tr_r)/slope;
0337 
0338         //float rsign = (trk_yemc>0) ? 1 : -1;
0339         float rsign = (0<=tr_phi&&tr_phi<TMath::Pi()) ? 1 : -1;
0340 
0341         TLine *l = new TLine(tr_z, tr_r, z_end, rsign*r_end);
0342         l->SetLineStyle(5);
0343         l->SetLineColor(2);
0344         v_tr_zr.push_back(l);
0345       }
0346       if(ev_skip) continue;
0347     }
0348 
0349 
0350     trackTree->GetEntry(i);
0351     caloTree->GetEntry(i);
0352     siClusTree->GetEntry(i);
0353 
0354     
0355 
0356     ////////////////////////
0357     // clear vector
0358     for(size_t iarc=0; iarc<v_trk_xy.size(); iarc++){
0359       delete v_trk_xy[iarc];
0360       delete v_trk_zr[iarc];
0361     }
0362     v_trk_xy.clear();
0363     v_trk_zr.clear();
0364 
0365 
0366     ////////////////////////
0367     //--if (evt != eventnumber) continue;
0368 
0369     // SiSeed Hits
0370     std::vector<float> xs, ys, zs, rs; // cluster
0371     // Append all clusters for this event
0372     if (trackcluster_x && trackcluster_y && trackcluster_z) {
0373       for (size_t j = 0; j < trackcluster_x->size() && j < trackcluster_y->size() && j < trackcluster_z->size(); ++j) {
0374         xs.push_back(trackcluster_x->at(j));
0375         ys.push_back(trackcluster_y->at(j));
0376         zs.push_back(trackcluster_z->at(j));
0377 
0378         float sign = (trackcluster_y->at(j)>0)?1:-1;
0379         float r = sign* std::sqrt(trackcluster_x->at(j)*trackcluster_x->at(j) + trackcluster_y->at(j)*trackcluster_y->at(j));
0380         rs.push_back(r);
0381       }
0382     }
0383 
0384     // calo Clusters
0385     std::vector<float> xs_calo, ys_calo, zs_calo, rs_calo, phis_calo;
0386     if (calo_x && calo_y && calo_z) {
0387       for (size_t j = 0; j < calo_x->size() && j < calo_y->size() && j < calo_z->size(); ++j) {
0388         xs_calo.push_back(calo_x->at(j));
0389         ys_calo.push_back(calo_y->at(j));
0390         zs_calo.push_back(calo_z->at(j));
0391         phis_calo.push_back(atan2(calo_y->at(j), calo_x->at(j)));
0392 
0393         float sign = (calo_y->at(j)>0)?1:-1;
0394         float r = sign*std::sqrt(calo_x->at(j)*calo_x->at(j) + calo_y->at(j)*calo_y->at(j));
0395         rs_calo.push_back(r);
0396         
0397         //// Update rz_track_z_emc and rz_track_r_emc from caloTree
0398         //rz_track_z_emc.push_back(calo_z->at(j));
0399         //rz_track_r_emc.push_back(r);
0400       }
0401     }
0402 
0403     // track trajectory
0404     // Store track projection to EMC  for this event (one per track)
0405     std::vector<float> xs_track_pemc, ys_track_pemc, zs_track_pemc, phis_track_pemc;
0406     if (track_x_emc && track_y_emc && track_z_emc) {
0407       for (size_t j = 0; j < track_x_emc->size() && j < track_y_emc->size() && j < track_z_emc->size(); ++j) {
0408         xs_track_pemc.push_back(track_x_emc->at(j));
0409         ys_track_pemc.push_back(track_y_emc->at(j));
0410         zs_track_pemc.push_back(track_z_emc->at(j));
0411 
0412         float phi_pemc = atan2(track_y_emc->at(j), track_x_emc->at(j));
0413         phis_track_pemc.push_back( phi_pemc );
0414 
0415         if (calo_x && calo_y && calo_z) {
0416           for (size_t k = 0; k < calo_x->size() && k < calo_y->size(); ++k) {
0417             float phi_emc = atan2(calo_y->at(k), calo_x->at(k));
0418           }
0419         }
0420       }
0421     }
0422     
0423     // to show the track trajectories by using track parameters (track_x, track_y, track_z, pt, phi, eta, and field B=1.4T)
0424     std::vector<float> cxs_trackarc, cys_trackarc, crs_trackarc, crs_matchflag;
0425     if (track_x && track_y && track_z && track_phi && track_eta && track_pt) {
0426       size_t ntracks = std::min({track_x->size(), track_y->size(), track_z->size(), 
0427                                  track_pt->size()});
0428 
0429       for (size_t j = 0; j < ntracks; ++j) {
0430         float& trk_x   = track_x->at(j);
0431         float& trk_y   = track_y->at(j);
0432         float& trk_z   = track_z->at(j);
0433         float& trk_phi = track_phi->at(j);
0434         float& trk_eta = track_eta->at(j);
0435         float& trk_pt  = track_pt->at(j);
0436         int&   trk_c   = track_charge->at(j);
0437         float  trk_the = 2.*std::atan( std::exp(-trk_eta) );
0438 
0439         float& trk_yemc= (track_y_emc->at(j));
0440 
0441         ////////////////////////
0442         // X-Y trajectory
0443         float R, cx, cy, arcphi0, arcphi1;
0444         calcTrajectoryXY(trk_x, trk_y, trk_pt, trk_phi, trk_c,
0445                          R, cx, cy, arcphi0, arcphi1);
0446 
0447 
0448          
0449         
0450         const float conv_ang = 180./TMath::Pi();
0451         float arcphi0_ang = arcphi0*conv_ang; 
0452         float arcphi1_ang = arcphi1*conv_ang; 
0453         
0454         float phi_ang     = trk_phi*conv_ang; if(phi_ang<0) phi_ang += 360.;
0455         std::cout<<"px, py, R = "<<trk_pt<<" "<<trk_phi<<" "<<R<<" "<<cx<<" "<<cy<<",   "
0456                  <<phi_ang<<" "<<arcphi0_ang<<" "<<arcphi1_ang<<" "<<trk_x<<" "<<trk_y<<std::endl;
0457 
0458 
0459 
0460         TArc *trk_xy = new TArc(cx, cy, R, arcphi0_ang, arcphi1_ang); //, ang0*180/TMath::Pi(), ang1*180/TMath::Pi());
0461         trk_xy->SetFillStyle(20);
0462         trk_xy->SetLineWidth(1);
0463         trk_xy->SetLineColor(1);
0464         v_trk_xy.push_back(trk_xy);
0465 
0466         cxs_trackarc.push_back(cx);
0467         cys_trackarc.push_back(cy);
0468         crs_trackarc.push_back(R);
0469         
0470         ////////////////////////
0471         // R-Z trajectory
0472         float r_end = 100.;
0473 
0474         float trk_r = std::sqrt(trk_x*trk_x+trk_y*trk_y);
0475 
0476         float slope = std::tan(trk_the);
0477         float z_end = (r_end-trk_r)/slope;
0478 
0479         //float rsign = (trk_yemc>0) ? 1 : -1;
0480         float rsign = (0<=trk_phi&&trk_phi<TMath::Pi()) ? 1 : -1;
0481 
0482         TLine *l = new TLine(trk_z, trk_r, z_end, rsign*r_end);
0483         v_trk_zr.push_back(l);
0484 
0485         ////////////////////////////
0486         // check the association
0487         float& z_pemc   = track_z_emc->at(j);
0488         float  phi_pemc = atan2(track_y_emc->at(j), track_x_emc->at(j));
0489         phis_track_pemc.push_back( phi_pemc );
0490 
0491         uint flag=0;
0492         if (calo_x && calo_y && calo_z) {
0493           for (size_t k = 0; k < calo_x->size() && k < calo_y->size(); ++k) {
0494             float& z_emc   = calo_z->at(k);
0495             float  dz      = z_emc - z_pemc;
0496             float  phi_emc = atan2(calo_y->at(k), calo_x->at(k));
0497             float  dphi = phi_emc - phi_pemc;
0498 
0499             if( fabs(dz)<5){
0500               float resi = (trk_c * dphi) - (0.18*TMath::Power(trk_pt, -0.986));
0501               float resi2= (trk_c * dphi) - (0.67*TMath::Power(trk_pt, -0.986));
0502               if( fabs( (trk_c * dphi) - (0.18*TMath::Power(trk_pt, -0.986)) ) < 0.1 ) { flag|=1; }
0503               if( fabs( (trk_c * dphi) - (0.67*TMath::Power(trk_pt, -0.986)) ) < 0.1 ) { flag|=2; }
0504               std::cout<<"dphi : "<<dphi<<" "<<trk_pt<<" "<<resi<<" "<<resi2<<std::endl;
0505             }
0506           }
0507         }
0508         crs_matchflag.push_back(flag);
0509         std::cout<<"flag : "<<flag<<std::endl;
0510 
0511         if(flag>0) {
0512           trk_xy->SetLineColor(flag+5);
0513           l->SetLineColor(flag+5);
0514         }
0515       }
0516     }
0517 
0518 
0519     ///////////////////////
0520     // draw SiSeed Hits
0521     if(gr!=nullptr){ delete gr; }
0522     gr = new TGraph(xs.size(), &xs[0], &ys[0]);
0523     gr->SetTitle(Form("Track Clusters X vs Y for Event %d;X [cm];Y [cm]", eventnumber));
0524     gr->SetMarkerStyle(20);
0525     gr->SetMarkerSize(1.0);
0526     gr->SetMarkerColor(kBlue+1);
0527 
0528     if(gr_rz!=nullptr){ delete gr_rz; }
0529     gr_rz = new TGraph(zs.size(), &zs[0], &rs[0]);
0530     gr_rz->SetTitle(Form("Track Clusters R vs Z for Event %d;X [cm];Y [cm]", eventnumber));
0531     gr_rz->SetMarkerStyle(20);
0532     gr_rz->SetMarkerSize(1.0);
0533     gr_rz->SetMarkerColor(kBlue+1);
0534 
0535     if(gr2_xyz!=nullptr){ delete gr2_xyz; gr2_xyz = nullptr; }
0536     if(zs.size()>0){
0537       gr2_xyz = new TGraph2D(zs.size(), &zs[0], &xs[0], &ys[0]);
0538       gr2_xyz->SetTitle(Form("Track Clusters X, Y, Z for Event %d;X [cm];Y [cm]", eventnumber));
0539       gr2_xyz->SetMarkerStyle(20);
0540       gr2_xyz->SetMarkerSize(1.0);
0541       gr2_xyz->SetMarkerColor(kBlue+1);
0542     }
0543 
0544   
0545     // Draw EMCal hits if available
0546     if(gr_calo!=nullptr){ delete gr_calo; }
0547     gr_calo = new TGraph(xs_calo.size(), &xs_calo[0], &ys_calo[0]);
0548     gr_calo->SetMarkerStyle(20); // star marker
0549     gr_calo->SetMarkerSize(1.0);
0550     gr_calo->SetMarkerColor(kRed+1);
0551 
0552     if(gr_rz_calo!=nullptr){ delete gr_rz_calo; }
0553     gr_rz_calo = new TGraph(zs_calo.size(), &zs_calo[0], &rs_calo[0]);
0554     gr_rz_calo->SetMarkerStyle(20); // star marker
0555     gr_rz_calo->SetMarkerSize(1.0);
0556     gr_rz_calo->SetMarkerColor(kRed+1);
0557 
0558     if(gr_phiz_calo!=nullptr){ delete gr_phiz_calo; }
0559     gr_phiz_calo = new TGraph(zs_calo.size(), &zs_calo[0], &phis_calo[0]);
0560     gr_phiz_calo->SetMarkerStyle(20); // star marker
0561     gr_phiz_calo->SetMarkerSize(1.0);
0562     gr_phiz_calo->SetMarkerColor(kRed+1);
0563 
0564     if(gr_phiz_pemc!=nullptr){ delete gr_phiz_pemc; }
0565     gr_phiz_pemc = new TGraph(zs_track_pemc.size(), &zs_track_pemc[0], &phis_track_pemc[0]);
0566     gr_phiz_pemc->SetMarkerStyle(24); // star marker
0567     gr_phiz_pemc->SetMarkerSize(1.2);
0568     gr_phiz_pemc->SetMarkerColor(kBlue+1);
0569 
0570     //trajectory
0571 
0572     std::cout<<" Nsi : "<<xs.size()
0573              <<", Ncalo : "<<xs_calo.size()
0574              <<", Ntrk : "<<v_trk_xy.size()
0575              <<", Ntr : "<<v_tr_xy.size()
0576              <<std::endl;
0577 
0578     p1->cd(1);
0579     //frm1->Draw();
0580     gPad->DrawFrame(-15,-15,15,15);
0581     for (size_t il = 0; il < v_elli_xy.size()-1; ++il) {
0582       v_elli_xy[il]->Draw();
0583       //std::cout<<il<<" "<<v_elli_xy[il]->GetR1()<<std::endl;
0584     }
0585     if(gr->GetN()>0) gr->Draw("P");
0586     for(size_t iarc=0; iarc<v_trk_xy.size(); iarc++){
0587       v_trk_xy[iarc]->Draw("only");
0588     }
0589     for(size_t iarc=0; iarc<v_tr_xy.size(); iarc++){
0590       std::cout<<"aa"<<std::endl;
0591       v_tr_xy[iarc]->Draw("only");
0592     }
0593 
0594     p1->cd(2);
0595     //frm1->Draw();
0596     gPad->DrawFrame(-28,-15,28,15);
0597     for (size_t il = 0; il < v_line_rz.size()-2; ++il) {
0598       v_line_rz[il]->Draw();
0599     }
0600     if(gr->GetN()>0) gr_rz->Draw("P");
0601     for(size_t iarc=0; iarc<v_trk_zr.size(); iarc++){
0602       v_trk_zr[iarc]->Draw();
0603     }
0604     for(size_t iarc=0; iarc<v_tr_zr.size(); iarc++){
0605       v_tr_zr[iarc]->Draw();
0606     }
0607 
0608     p1->cd(3);
0609     //frm1->Draw();
0610     gPad->DrawFrame(-100,-100,100,100);
0611     for (size_t il = 0; il < v_elli_xy.size(); ++il) {
0612       v_elli_xy[il]->Draw();
0613     }
0614     if(gr->GetN()     >0) gr->Draw("P");
0615     if(gr_calo->GetN()>0) gr_calo->Draw("P");
0616     for(size_t iarc=0; iarc<v_trk_xy.size(); iarc++){
0617       v_trk_xy[iarc]->Draw("only");
0618     }
0619     for(size_t iarc=0; iarc<v_tr_xy.size(); iarc++){
0620       v_tr_xy[iarc]->Draw("only");
0621     }
0622 
0623     p1->cd(4);
0624     //frm1->Draw();
0625     gPad->DrawFrame(-150,-100,150,100);
0626     for (size_t il = 0; il < v_line_rz.size(); ++il) {
0627       v_line_rz[il]->Draw();
0628     }
0629     if(gr_rz->GetN()     >0) gr_rz->Draw("P");
0630     if(gr_rz_calo->GetN()>0) gr_rz_calo->Draw("P");
0631     for(size_t iarc=0; iarc<v_trk_zr.size(); iarc++){
0632       v_trk_zr[iarc]->Draw();
0633     }
0634     for(size_t iarc=0; iarc<v_tr_zr.size(); iarc++){
0635       v_tr_zr[iarc]->Draw();
0636     }
0637   
0638   
0639     can->cd(2);
0640     gPad->DrawFrame(-150, -3.2, 150, 3.2);
0641     gr_phiz_calo->Draw("P");
0642     gr_phiz_pemc->Draw("P");
0643     //--frame3d->Draw();
0644     //--if(gr2_xyz!=nullptr) {
0645     //--  //gr2_xyz->GetXaxis()->SetRangeUser(-30,30);
0646     //--  //gr2_xyz->GetYaxis()->SetRangeUser(-15,15);
0647     //--  //gr2_xyz->GetZaxis()->SetRangeUser(-15,15);
0648     //--  //gr2_xyz->Draw("Psame");
0649     //--  gr2_xyz->Draw("sameP");
0650     //--}
0651  
0652 
0653     can->Modified();
0654     can->Update();
0655 
0656     std::cin.getline(input, 256);
0657     if (input[0] == 'q' || input[0] == 'Q') {
0658       break;
0659     }
0660 
0661   }
0662 
0663 
0664   //can->Draw();
0665   //can->Update();
0666   file->Close();
0667 }
0668 
0669 /*
0670 void LoopEventDisplay()
0671 {
0672   int eventnumber = 0;
0673   char input[256];
0674   while (true) {
0675     DrawEventDisplay(eventnumber);
0676     printf("Press Enter to view next event, or 'q' then Enter to quit: ");
0677     std::cin.getline(input, 256);
0678     if (input[0] == 'q' || input[0] == 'Q') {
0679       break;
0680     }
0681     eventnumber++;
0682   }
0683 }
0684 */
[ 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!