Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:13:57 

0001 #define clustEvent_cc
0002 
0003 #include "clustEvent.hh"
0004 
0005 clustEvent::clustEvent()
0006 {
0007   colorset.push_back( kOrange - 3 );
0008   colorset.push_back( kOrange + 3 );
0009   colorset.push_back( kPink - 3 );
0010   colorset.push_back( kPink + 3 );
0011   colorset.push_back( kAzure - 3 );
0012   colorset.push_back( kAzure + 3 );
0013   colorset.push_back( kTeal - 3 );
0014   colorset.push_back( kTeal + 3 );
0015   colorset.push_back( kSpring - 3 );
0016   colorset.push_back( kSpring + 3 );
0017   colorset.push_back( kViolet - 3 );
0018   colorset.push_back( kViolet + 3 );
0019   
0020   // suggested by ChatGPT
0021   // (102, 255, 153)
0022   //   (255, 153, 102)
0023   //   (153, 102, 255)
0024   //   (255, 230, 102)
0025   //   (102, 204, 255)
0026 }
0027 
0028 void clustEvent::clear()
0029 {
0030   for (auto itr = vtrack.begin(); itr != vtrack.end(); ++itr)
0031     {
0032       delete *itr;
0033     }
0034   vtrack.clear();
0035 }
0036 
0037 void clustEvent::makeonetrack()
0038 {
0039   for (unsigned int i = 0; i < vtrack.size(); i++)
0040     {
0041       for (unsigned int j = i + 1; j < vtrack.size(); j++)
0042     {
0043       bool share_point = ((vtrack[i]->p1 == vtrack[j]->p1) || (vtrack[i]->p2 == vtrack[j]->p2));
0044       share_point = share_point || ((vtrack[i]->p1 == vtrack[j]->p2) || (vtrack[i]->p2 == vtrack[j]->p1));
0045       if (share_point)
0046         {
0047           if (fabs(vtrack[i]->dca_2d) > fabs(vtrack[j]->dca_2d))
0048         {
0049           vtrack[i]->is_deleted = true;
0050           continue;
0051         }
0052           else
0053         {
0054           vtrack[j]->is_deleted = true;
0055         }
0056         }
0057     }
0058     }
0059 }
0060 
0061 void clustEvent::dca_check(bool Debug)
0062 {
0063   for (unsigned int i = 0; i < vtrack.size(); i++)
0064     {
0065       if (vtrack[i]->dca_2d > dca2d_max)
0066     vtrack[i]->dca2d_range_out = true;
0067       if (vtrack[i]->dca_2d < dca2d_min)
0068     vtrack[i]->dca2d_range_out = true;
0069 
0070       if (vtrack[i]->dca[2] > dcaz_max)
0071     vtrack[i]->dcaz_range_out = true;
0072       if (vtrack[i]->dca[2] < dcaz_min)
0073     vtrack[i]->dcaz_range_out = true;
0074 
0075       if (vtrack[i]->dcaz_range_out || vtrack[i]->dca2d_range_out)
0076     vtrack[i]->dca_range_out = true;
0077 
0078       // // if (Debug)
0079       // {
0080       //     if (vtrack[i]->dca_mean[2] > vtrack[i]->p1.z() && vtrack[i]->dca_mean[2] < vtrack[i]->p2.z())
0081       //         vtrack[i]->dca_range_out = true;
0082       //     if (vtrack[i]->dca_mean[2] > vtrack[i]->p2.z() && vtrack[i]->dca_mean[2] < vtrack[i]->p1.z())
0083       //         vtrack[i]->dca_range_out = true;
0084       // }
0085     }
0086 }
0087 
0088 void clustEvent::truth_check()
0089 {
0090   for (unsigned int i = 0; i < vtruth.size(); i++)
0091     {
0092       double p = sqrt(vtruth[i]->m_truthpt * vtruth[i]->m_truthpt + vtruth[i]->m_truthpz * vtruth[i]->m_truthpz);
0093 
0094       if (vtruth[i]->m_status == 1 && fabs(vtruth[i]->m_trutheta) < 2 && p > (50 * 1e-3))
0095     {
0096       if (fabs(vtruth[i]->m_truthz_out) < 22.8)
0097         {
0098           vtruth[i]->is_intt = true;
0099         }
0100     }
0101     }
0102 }
0103 
0104 void clustEvent::cluster_check()
0105 {
0106   for (unsigned int itrt = 0; itrt < vtruth.size(); itrt++)
0107     {
0108       for (unsigned int icls = 0; icls < vclus.size(); icls++)
0109     {
0110       if (vclus[icls].layer < 2)
0111         continue;
0112 
0113       double d_phi = vclus[icls].getphi_clus() - vtruth[itrt]->m_truthphi;
0114       d_phi -= -2 * (M_PI * (int)(d_phi / (M_PI)));
0115 
0116       if (fabs(d_phi) < (0.0014 * 3))
0117         {
0118           vtruth[itrt]->have_cluster = true;
0119           break;
0120         }
0121     }
0122     }
0123 }
0124 
0125 void clustEvent::track_check()
0126 {
0127   for (unsigned int i = 0; i < vtruth.size(); i++)
0128     {
0129       for (unsigned int itrk = 0; itrk < vtrack.size(); itrk++)
0130     {
0131       double d_phi = vtrack[itrk]->getphi() - vtruth[i]->m_truthphi;
0132       d_phi -= -2 * (M_PI * (int)(d_phi / (M_PI)));
0133 
0134       if (fabs(d_phi) < (3.3e-4 * 3))
0135         {
0136           vtruth[i]->have_track = true;
0137           break;
0138         }
0139     }
0140     }
0141 }
0142 
0143 void clustEvent::charge_check()
0144 {
0145   for (unsigned int i = 0; i < vtruth.size(); i++)
0146     {
0147       if ((abs(vtruth[i]->m_truthpid) != 22 && abs(vtruth[i]->m_truthpid) != 2112 && abs(vtruth[i]->m_truthpid) != 130 && abs(vtruth[i]->m_truthpid) != 310 && abs(vtruth[i]->m_truthpid) != 111))
0148     vtruth[i]->is_charged = true;
0149     }
0150 }
0151 
0152 void clustEvent::draw_frame(int mode, bool is_preliminary )
0153 {
0154   TH1 *frame;
0155   string mag;
0156 
0157   if (mag_on)
0158     mag = "B-on";
0159   else
0160     mag = "B-off";
0161 
0162   string title;
0163   if( mode == 0 ) // x-y
0164     {
0165       if (run_no == 1)
0166     title = Form("x-y plane {MC(%s) event %d};x [cm];y [cm]", mag.c_str(), ievt);
0167       else
0168     title = Form("x-y plane {run%d(%s) event %d};x [cm];y [cm]", run_no, mag.c_str(), ievt);
0169       
0170       frame = gPad->DrawFrame(-15, -15, 15, 15, title.c_str());
0171       
0172       this->draw_signature( is_preliminary );
0173     }
0174   else if (mode == 1) // z-r
0175     {
0176       if (run_no == 1)
0177     title = Form("z-r plane {MC(%s) event %d};z [cm];r [cm]", mag.c_str(), ievt);
0178       else
0179     title = Form("z-r plane {run%d(%s) event %d};z [cm];r [cm]", run_no, mag.c_str(), ievt);
0180 
0181       //frame = gPad->DrawFrame(-25, -15, 25, 15, title.c_str());
0182       frame = gPad->DrawFrame(-30, -15, 30, 15, title.c_str());
0183       this->draw_date();
0184       
0185     }
0186 
0187   frame->GetXaxis()->CenterTitle();
0188   frame->GetYaxis()->CenterTitle();
0189   
0190 }
0191 
0192 string clustEvent::GetDate()
0193 {
0194   return "9/13/2024";
0195   
0196   TDatime dt;
0197   int year  = dt.GetYear();
0198   int month = dt.GetMonth();
0199   int day   = dt.GetDay();
0200 
0201   // format: mm/dd/yyyy
0202   std::stringstream ss;
0203   ss << month << "/" << day << "/" << year;
0204 
0205   return ss.str();
0206 }
0207 
0208 void clustEvent::draw_signature( bool is_preliminary = false )
0209 {
0210 
0211   ///////////////////////////////////////////////////////////////////////////////
0212   // Writting words in the canvas                                              //
0213   ///////////////////////////////////////////////////////////////////////////////
0214   TLatex* tex = new TLatex();
0215   double line_height = 0.05;
0216   double first_margin = 0.005;
0217   double top_margin = 0.1;  
0218   double pos_y = 1.0 - top_margin + first_margin;// - line_height;
0219   pos_y = 1.0 - line_height + 0.004 + 0.01;
0220 
0221   // sPHENIX Internal or sPHENIX Prelimnary
0222   // pos_y -= line_height - first_margin + 0.025;
0223   double pos_x = 0.175;
0224   if( is_preliminary == false )
0225     {
0226       tex->DrawLatexNDC( pos_x, pos_y, "#it{#bf{sPHENIX}} Internal" );
0227     }
0228   else
0229     {
0230       if( bco_intt == 396907563635 ||
0231       bco_intt == 396904684235 ||
0232       bco_intt == 396904268555 )
0233     tex->DrawLatexNDC( pos_x, pos_y, "#it{#bf{sPHENIX}} Preliminary" );
0234     }
0235 
0236     // p+p 200 GeV
0237   pos_y -= line_height;
0238   tex->DrawLatexNDC( pos_x, pos_y,
0239              ( string("Run-24 #it{p+p} 200 GeV Run" ) + to_string( run_no ) ).c_str()
0240              );
0241 
0242   // tex->DrawLatexNDC( 0.2, pos_y, "Run-24 #it{p+p} 200 GeV" ); // no run number version 
0243 
0244 }
0245 
0246 void clustEvent::draw_date()
0247 {
0248   TLatex* tex = new TLatex();
0249   double line_height = 0.05;
0250   double first_margin = 0.005;
0251   double top_margin = 0.1;  
0252   double pos_y = 1.0 - top_margin + first_margin;// - line_height;
0253   pos_y = 1.0 - line_height + 0.004 + 0.01;
0254 
0255   // top line
0256   if( run_no == 50889 )
0257     {
0258       //pos_y -= line_height;
0259       tex->DrawLatexNDC( 0.175, pos_y, "INTT Streaming readout" );
0260 
0261     }
0262 
0263   // second line
0264   pos_y -= line_height;
0265   // Date
0266   tex->DrawLatexNDC( 0.7, pos_y,
0267              string("#it{" + GetDate() + "}").c_str() );
0268 
0269   // tex->DrawLatexNDC( 0.175, pos_y,
0270   //             (string("BCO ") + to_string( bco_intt ) ).c_str()
0271   //             );
0272   tex->DrawLatexNDC( 0.175, pos_y, "Single crossing" );
0273 
0274   return;
0275 }
0276 
0277 void clustEvent::draw_intt(int mode)
0278 {
0279 
0280   const double kLayer_radii[4] = {7.1888, 7.800, 9.680, 10.330};
0281   if (mode == 0)
0282     {
0283       for (int i = 0; i < 4; i++)
0284     {
0285       auto circle = new TEllipse(0.0, 0.0, kLayer_radii[i]);
0286       circle->SetLineColorAlpha(kGray, 0.5);
0287       circle->SetLineWidth(2);
0288       circle->SetFillStyle(0);
0289       circle->Draw("same");
0290     }
0291     }
0292   else if (mode == 1)
0293     {
0294       for (int j = 0; j < 2; j++)
0295     {
0296       for (int i = 0; i < 4; i++)
0297         {
0298           TLine *line = new TLine(-22.8, (2 * j - 1) * kLayer_radii[i], 22.8, (2 * j - 1) * kLayer_radii[i]);
0299           line->SetLineColorAlpha(kGray, 0.5);
0300           line->SetLineWidth(2);
0301           line->Draw("same");
0302         }
0303     }
0304     }
0305 }
0306 
0307 void clustEvent::draw_tracklets(int mode = 0, bool does_overlay = false, int color = kBlack, bool dca_range_cut = false, bool is_deleted = false, bool reverse = false)
0308 {
0309 
0310   bool condition = !reverse;
0311   for (unsigned int i = 0; i < vtrack.size(); i++)
0312     {
0313       if (vtrack[i]->dca_range_out == condition && dca_range_cut == true)
0314     {
0315       continue;
0316     }
0317 
0318       if (vtrack[i]->is_deleted == condition && is_deleted == true)
0319     {
0320       continue;
0321     }
0322 
0323       TGraph *g_temp = new TGraph();
0324       g_temp->SetMarkerColor(color);
0325       g_temp->SetMarkerStyle(20);
0326       g_temp->SetMarkerSize(0.8);
0327 
0328       if (mode == 0)
0329     {
0330       g_temp->SetPoint(0, vtrack[i]->p1.x(), vtrack[i]->p1.y());
0331       g_temp->SetPoint(1, vtrack[i]->p2.x(), vtrack[i]->p2.y());
0332     }
0333       else if (mode == 1)
0334     {
0335       g_temp->SetPoint(0, vtrack[i]->p1.z(), vtrack[i]->getpointr(1));
0336       g_temp->SetPoint(1, vtrack[i]->p2.z(), vtrack[i]->getpointr(2));
0337     }
0338       g_temp->Draw("PLsame");
0339     }
0340 }
0341 
0342 void clustEvent::draw_trackline(int mode,
0343                 bool does_overlay,
0344                 //int color,
0345                 bool dca_range_cut,
0346                 bool is_deleted,
0347                 bool reverse)
0348 {
0349   bool condition = !reverse;    
0350   int color_ = 0;
0351 
0352   bool is_first = true;
0353   for (unsigned int i = 0; i < vtrack.size(); i++)
0354     {
0355 
0356       if (vtrack[i]->dca_range_out == condition && dca_range_cut == true)
0357     {
0358       continue;
0359     }
0360 
0361       if (vtrack[i]->is_deleted == condition && is_deleted == true)
0362     {
0363       continue;
0364     }
0365     
0366       color_ = color_ % colorset.size();
0367 
0368       TGraph *g_temp = new TGraph();
0369       g_temp->SetMarkerColor(colorset[color_]);
0370       // g_temp->SetLineColor(color);
0371       g_temp->SetMarkerStyle(20);
0372       g_temp->SetMarkerSize(0.8);
0373 
0374       if (mode == 0)
0375     {
0376       vtrack[i]->track_xy->SetLineColorAlpha(colorset[color_], 0.5);
0377       // vtrack[i]->track_xy->SetLineColorAlpha(color, 0.5);
0378       if (vtrack[i]->dca_mean[0] < vtrack[i]->p1.x())
0379         {
0380           vtrack[i]->track_xy->SetRange(vtrack[i]->dca_mean[0], 15);
0381         }
0382       else
0383         {
0384           vtrack[i]->track_xy->SetRange(-15, vtrack[i]->dca_mean[0]);
0385         }
0386 
0387       vtrack[i]->track_xy->Draw("Lsame");
0388       g_temp->SetPoint(0, vtrack[i]->dca_mean[0], vtrack[i]->dca_mean[1]);
0389       g_temp->SetPoint(1, vtrack[i]->p1.x(), vtrack[i]->p1.y());
0390       g_temp->SetPoint(2, vtrack[i]->p2.x(), vtrack[i]->p2.y());
0391     }
0392       else if (mode == 1)
0393     {
0394 
0395 
0396       if( isinf(vtrack[i]->track_rz->GetParameter(1)) )
0397         {
0398           TLine *line = nullptr;
0399           if (vtrack[i]->p1.y() > 0)
0400         line = new TLine(vtrack[i]->dca_mean[2], vtrack[i]->getpointr(3), vtrack[i]->dca_mean[2], 15);
0401           else
0402         line = new TLine(vtrack[i]->dca_mean[2], -15, vtrack[i]->dca_mean[2], vtrack[i]->getpointr(3));
0403 
0404           // line->SetLineColorAlpha(color, 0.5);
0405           line->SetLineColorAlpha(colorset[color_], 0.5);
0406           line->SetLineWidth(2);
0407           line->Draw("same");
0408         }
0409       else
0410         {
0411 
0412           vtrack[i]->track_rz->SetLineColorAlpha(colorset[color_], 0.5);
0413           // vtrack[i]->track_rz->SetLineColorAlpha(color, 0.5);
0414 
0415           int y = 0;
0416           if (vtrack[i]->p1.y() > 0)
0417         y = 15;
0418           else
0419         y = -15;
0420 
0421           double r1 = vtrack[i]->track_rz_inv->Eval(vtrack[i]->getpointr(3));
0422           double r2 = vtrack[i]->track_rz_inv->Eval(y);
0423 
0424           if (r1 <= r2)
0425         vtrack[i]->track_rz->SetRange(r1, r2);
0426           else
0427         vtrack[i]->track_rz->SetRange(r2, r1);
0428 
0429           vtrack[i]->track_rz->Draw("Lsame");
0430         }
0431 
0432       if( is_first ) // drawing graph only once is enough for z-r
0433         {
0434           g_temp->SetPoint(0, vtrack[i]->dca_mean[2], vtrack[i]->getpointr(3));
0435           is_first = false;
0436         }
0437 
0438       // draw hit strip
0439       vector < TLine * > zlines;
0440       for (int j = 0; j < 2; j++)
0441         {
0442           TLine *line = new TLine(vtrack[i]->zline[j][0],
0443                       vtrack[i]->getpointr(1 + j),
0444                       vtrack[i]->zline[j][1],
0445                       vtrack[i]->getpointr(1 + j) );
0446           line->SetLineColor(colorset[color_]);
0447           line->SetLineWidth(3);
0448           zlines.push_back(line);
0449         }
0450 
0451       for ( auto& line : zlines )
0452         line->Draw("same");
0453       
0454       // g_temp->SetPoint(1, vtrack[i]->p1.z(), vtrack[i]->getpointr(1));
0455       // g_temp->SetPoint(2, vtrack[i]->p2.z(), vtrack[i]->getpointr(2));
0456     }
0457     
0458       if( g_temp->GetN() > 0 )
0459     g_temp->Draw("Psame");
0460       
0461       color_++;
0462 
0463       if( color_ >= colorset.size() )
0464     color_ = 0;
0465       
0466     }
0467 
0468 }
0469 
0470 
0471 double clustEvent::rad_deg(double Rad)
0472 {
0473   Rad = Rad / M_PI * 180;
0474   if (Rad < 0)
0475     Rad += 360;
0476 
0477   return Rad;
0478 }
0479 
0480 void clustEvent::draw_curve2(int mode, double px, double py, double pz, double rad, double cx, double cy, int color, const vector<double> &V)
0481 {
0482   if (V.size() == 2)
0483     {
0484       // cout<<"cross : "<<V.at(0)<<", "<<V.at(1)<<endl;
0485       double phi1 = atan2(V.at(1) - cy, V.at(0) - cx);
0486       phi1 = rad_deg(phi1);
0487 
0488       double phic = atan2(dca_mean[1] - cy, dca_mean[0] - cx);
0489       phic = rad_deg(phic);
0490       double phi_min, phi_max;
0491 
0492       phi_min = (phi1 < phic) ? phi1 : phic;
0493       phi_max = (phi1 < phic) ? phic : phi1;
0494       TEllipse *circle = nullptr;
0495       // TEllipse *circle2 = nullptr;
0496       if ((phi_max - phi_min) > 180)
0497     circle = new TEllipse(cx, cy, rad, 0, -(360 - phi_max), phi_min);
0498       else
0499     circle = new TEllipse(cx, cy, rad, 0, phi_min, phi_max);
0500 
0501       if (mode == 0)
0502     {
0503       circle->SetLineColorAlpha(color, 0.5);
0504       circle->SetFillStyle(0);
0505       circle->SetLineWidth(2);
0506       circle->SetLineStyle(1);
0507       circle->SetNoEdges(1);
0508       circle->Draw("same");
0509     }
0510       else if (mode == 1)
0511     {
0512       circle->SetLineColorAlpha(color, 0.9);
0513       circle->SetFillStyle(0);
0514       circle->SetLineWidth(3);
0515       circle->SetLineStyle(9);
0516       circle->SetNoEdges(1);
0517       circle->Draw("same");
0518     }
0519     }
0520 
0521   // TEllipse *circle2 = new TEllipse(cx, cy, rad, 0, 0, 360);
0522   // circle2->SetLineColorAlpha(kGray, 0.2);
0523   // circle2->SetFillStyle(0);
0524   // circle2->SetLineWidth(2);
0525   // circle2->SetLineStyle(1);
0526   // circle2->SetNoEdges(1);
0527   // circle2->Draw("same");
0528 }
0529 
0530 void clustEvent::draw_trackcurve( int mode,
0531                   bool does_overlay,
0532                   bool dca_range_cut,
0533                   bool is_deleted,
0534                   bool reverse )
0535 {
0536   int color_ = 0;
0537   bool condition = !reverse;
0538 
0539   for (unsigned int i = 0; i < vtrack.size(); i++ )
0540     {
0541       if( color_ == colorset.size() )
0542     color_ = 0;
0543       
0544       if (vtrack[i]->dca_range_out == condition && dca_range_cut == true)
0545     {
0546       continue;
0547     }
0548 
0549       if (vtrack[i]->is_deleted == condition && is_deleted == true)
0550     {
0551       continue;
0552     }
0553       
0554       color_ = color_ % colorset.size();
0555 
0556       TGraph *g_temp = new TGraph();
0557       g_temp->SetMarkerColor(colorset[color_]);
0558       g_temp->SetMarkerStyle(20);
0559       g_temp->SetMarkerSize(0.8);
0560 
0561       if (mode == 0)
0562     {
0563       vector < double > cross_framepoint = get_crossframe(vtrack[i]->cx, vtrack[i]->cy, vtrack[i]->rad, vtrack[i]->p1.x(), vtrack[i]->p1.y());
0564       // cout<<"cross_framepoint.size() : "<<cross_framepoint.size()<<endl;
0565       // if(cross_framepoint.size() == 2)
0566       // cout<<"cross : "<<cross_framepoint.at(0)<<", "<<cross_framepoint.at(1)<<endl;
0567       draw_curve2(0,
0568               vtrack[i]->p_reco[0],
0569               vtrack[i]->p_reco[1],
0570               vtrack[i]->p_reco[2],
0571               vtrack[i]->rad,
0572               vtrack[i]->cx,
0573               vtrack[i]->cy,
0574               colorset[color_],
0575               cross_framepoint);
0576 
0577       // cout<<"p1 : "<<vtrack[i]->p1.x()<<"  "<<vtrack[i]->p1.y()<<endl;
0578 
0579       g_temp->AddPoint( vtrack[i]->dca_mean[0], vtrack[i]->dca_mean[1]);
0580       g_temp->AddPoint( vtrack[i]->p1.x(), vtrack[i]->p1.y());
0581       g_temp->AddPoint( vtrack[i]->p2.x(), vtrack[i]->p2.y());
0582 
0583       // cout << "cross_emc_size : " << vtrack[i]->cross_emc.size() << endl;
0584     }
0585 
0586       color_++;
0587       if( color_ >= colorset.size() )
0588     color_ = 0;
0589 
0590       g_temp->Draw("Psame");
0591     }
0592 }
0593 
0594 void clustEvent::draw_truthcurve(int mode = 0, bool does_overlay = false, int color = kBlack, bool only_far_cluster = false, bool only_far_track = false)
0595 {
0596   for (unsigned int i = 0; i < vtruth.size(); i++)
0597     {
0598       if (vtruth[i]->is_intt == false)
0599     continue;
0600 
0601       if (vtruth[i]->is_charged == false)
0602     continue;
0603 
0604       if (vtruth[i]->have_track == true && only_far_track)
0605     continue;
0606 
0607       if (vtruth[i]->have_cluster == true && only_far_cluster)
0608     continue;
0609 
0610       TGraph *g_temp = new TGraph();
0611       g_temp->SetMarkerColor(color);
0612       g_temp->SetLineColor(color);
0613       g_temp->SetMarkerStyle(20);
0614       g_temp->SetMarkerSize(0.8);
0615 
0616       if (mode == 0)
0617     {
0618       vector<double> cross_framepoint = get_crossframe(vtruth[i]->center[0], vtruth[i]->center[1], vtruth[i]->m_rad, vtruth[i]->m_truthpx, vtruth[i]->m_truthpy);
0619       draw_curve2(1, vtruth[i]->m_truthpx, vtruth[i]->m_truthpy, vtruth[i]->m_truthpz, vtruth[i]->m_rad, vtruth[i]->center[0], vtruth[i]->center[1], color, cross_framepoint);
0620 
0621       g_temp->SetPoint(0, vtruth[i]->m_xvtx, vtruth[i]->m_yvtx);
0622     }
0623 
0624       g_temp->Draw("Psame");
0625     }
0626 }
0627 
0628 void clustEvent::draw_truthline(int mode = 0, bool does_overlay = false, int color = kBlack, bool only_far_cluster = false, bool only_far_track = false)
0629 {
0630   for (unsigned int i = 0; i < vtruth.size(); i++)
0631     {
0632       if (vtruth[i]->is_intt == false)
0633     continue;
0634 
0635       if (vtruth[i]->is_charged == false)
0636     continue;
0637 
0638       if (vtruth[i]->have_track == true && only_far_track)
0639     continue;
0640 
0641       if (vtruth[i]->have_cluster == true && only_far_cluster)
0642     continue;
0643 
0644       TGraph *g_temp = new TGraph();
0645       g_temp->SetMarkerColor(color);
0646       g_temp->SetLineColor(color);
0647       g_temp->SetMarkerStyle(20);
0648       g_temp->SetMarkerSize(0.8);
0649 
0650       if (mode == 0)
0651     {
0652       vtruth[i]->truth_xy->SetLineColorAlpha(color, 0.8);
0653       vtruth[i]->truth_xy->SetLineStyle(9);
0654       vtruth[i]->truth_xy->SetLineWidth(3);
0655 
0656       if (vtruth[i]->m_truthpx > 0)
0657         {
0658           vtruth[i]->truth_xy->SetRange(vtruth[i]->m_xvtx, 15);
0659         }
0660       else
0661         {
0662           vtruth[i]->truth_xy->SetRange(-15, vtruth[i]->m_xvtx);
0663         }
0664       vtruth[i]->truth_xy->Draw("Lsame");
0665 
0666       g_temp->SetPoint(0, vtruth[i]->m_xvtx, vtruth[i]->m_yvtx);
0667     }
0668       else if (mode == 1)
0669     {
0670       vtruth[i]->truth_rz->SetLineColorAlpha(color, 0.8);
0671       vtruth[i]->truth_rz->SetLineStyle(9);
0672       vtruth[i]->truth_rz->SetLineWidth(3);
0673 
0674       int y = 0;
0675       if (vtruth[i]->m_truthpy > 0)
0676         y = 15;
0677       else
0678         y = -15;
0679       double r_min = vtruth[i]->truth_rz->GetX(vtruth[i]->getpointr(3));
0680       double r_max = vtruth[i]->truth_rz->GetX(y);
0681       if (r_min < r_max)
0682         vtruth[i]->truth_rz->SetRange(r_min, r_max);
0683       else
0684         vtruth[i]->truth_rz->SetRange(r_max, r_min);
0685 
0686       vtruth[i]->truth_rz->Draw("Lsame");
0687       g_temp->SetPoint(0, vtruth[i]->m_zvtx, vtruth[i]->getpointr(3));
0688     }
0689       g_temp->Draw("Psame");
0690     }
0691 }
0692 
0693 void clustEvent::draw_clusters(int mode, bool does_overlay, int color )
0694 {
0695   for (unsigned int i = 0; i < vclus.size(); i++)
0696     {
0697       TGraph *g_temp = new TGraph();
0698       g_temp->SetMarkerColorAlpha(color, 0.5);
0699       // g_temp->SetLineColor(color);
0700       g_temp->SetMarkerStyle(20);
0701       g_temp->SetMarkerSize(0.8);
0702       
0703       if (mode == 0)
0704     {
0705       g_temp->SetPoint(0, vclus[i].x, vclus[i].y);
0706       g_temp->Draw("psame");
0707     }
0708       else if (mode == 1)
0709     {
0710       TLine *line = new TLine( vclus[i].zline[0], vclus[i].r, vclus[i].zline[1], vclus[i].r );
0711       line->SetLineColorAlpha(color, 0.5);
0712       line->SetLineWidth(3);
0713       line->Draw("same");
0714       // g_temp->SetPoint(0, vclus[i].z, vclus[i].r);
0715     }
0716       // g_temp->Draw("psame");
0717     }
0718 }
0719 
0720 //////////////////////////////////
0721 //////////////////////////////////
0722 
0723 // vector<double> get_crossline(double x1, double y1, double r1)
0724 // {
0725 //     vector<double> cross_co;
0726 //     // ax + by + c = 0
0727 //     cross_co.push_back(2 * (cx - x1));                                                           // a
0728 //     cross_co.push_back(2 * (cy - y1));                                                           // b
0729 //     cross_co.push_back((x1 + cx) * (x1 - cx) + (y1 + cy) * (y1 - cy) + (rad + r1) * (rad - r1)); // c
0730 
0731 //     return cross_co;
0732 // }
0733 
0734 vector<double> clustEvent::get_crosspoint(const vector<double> &V, double cx, double cy, double rad, double p1_x, double p1_y) // ax + by + c = 0
0735 {
0736   vector<double> cross;
0737   double a = V[0];
0738   double b = V[1];
0739   double c = V[2];
0740 
0741   double d = fabs(a * cx + b * cy + c) / sqrt(a * a + b * b);
0742   // cout << "d ; " << d << " " << rad << endl;
0743 
0744   double theta = atan2(b, a);
0745 
0746   if (d > rad)
0747     {
0748       // cout << "d > rad" << endl;
0749     }
0750   else if (d == rad)
0751     {
0752       // cout << "d == rad" << endl;
0753 
0754       if (a * cx + b * cy + c > 0)
0755     theta += M_PI;
0756       cross.push_back(rad * cos(theta) + cx);
0757       cross.push_back(rad * sin(theta) + cy);
0758     }
0759   else
0760     {
0761       // cout << "d < rad" << endl;
0762       double alpha, beta, phi;
0763       phi = acos(d / rad);
0764       alpha = theta - phi;
0765       beta = theta + phi;
0766       if (a * cx + b * cy + c > 0)
0767     {
0768       alpha += M_PI;
0769       beta += M_PI;
0770     }
0771       // double temp_cross[2][2];
0772       vector<vector<double>> temp_cross = {{rad * cos(alpha) + cx, rad * cos(beta) + cx}, {rad * sin(alpha) + cy, rad * sin(beta) + cy}};
0773       // for (unsigned int j = 0; j < temp_cross.at(0).size(); j++)
0774       // {
0775       //     cout << "temp_cross : ";
0776       //     for (unsigned int i = 0; i < temp_cross.size(); i++)
0777       //         cout << temp_cross.at(j).at(i) << "  ";
0778       //     cout << endl;
0779       // }
0780       cross = get_closestpoint(temp_cross, p1_x, p1_y);
0781     }
0782   // cout << "cross size : " << cross.size() << endl;
0783   // for (int i = 0; i < cross.size(); i++)
0784   //     cout << cross[i] << endl;
0785   return cross;
0786 }
0787 
0788 vector<double> clustEvent::get_closestpoint(const vector<vector<double>> VV, double p1_x, double p1_y)
0789 {
0790   vector<double> closest;
0791   double pre_phi = 999;
0792   double phi_p1 = atan2(p1_y - dca_mean[1], p1_x - dca_mean[0]);
0793   for (unsigned int i = 0; i < VV.at(0).size(); i++)
0794     {
0795       // cout << "dca_mean : " << dca_mean[0] << "  " << dca_mean[1] << endl;
0796       // cout << "VV : " << VV[0][i] << "  " << VV[1][i] << endl;
0797       // cout << "pVV : " << VV[0][i] - dca_mean[0] << "  " << VV[1][i] - dca_mean[1] << endl;
0798       // cout << "p1  : " << p1.x() - dca_mean[0] << "  " << p1.y() - dca_mean[1] << endl;
0799       // cout << VV.at(i).at(0) << endl;
0800       double dot = (VV[0][i] - dca_mean[0]) * (p1_x - dca_mean[0]) + (VV[1][i] - dca_mean[1]) * (p1_y - dca_mean[1]);
0801       double temp_phi = atan2(VV[1][i] - dca_mean[1], VV[0][i] - dca_mean[0]);
0802       // cout << "dot : " << dot << endl;
0803       if (dot > 0)
0804     {
0805       if (fabs(temp_phi - phi_p1) < fabs(pre_phi - phi_p1))
0806         {
0807           closest.erase(closest.begin(), closest.end());
0808           // cout << "VV size : " << VV.size() << " " << VV.at(0).size() << endl;
0809           for (int j = 0; j < 2; j++)
0810         closest.push_back(VV[j][i]);
0811           pre_phi = temp_phi;
0812         }
0813     }
0814     }
0815   // closest = temp_closest;
0816   // for (unsigned int i = 0; i < VV.at(0).size(); i++)
0817   // {
0818   //     for (int j = 0; j < 2; j++)
0819   //         closest.push_back(VV[j][i]);
0820   // }
0821   // for (unsigned int i = 0; i < closest.size(); i++)
0822   //  // cout << closest[i] << endl;
0823 
0824   return closest;
0825 }
0826 
0827 vector<double> clustEvent::get_crossframe(double cx, double cy, double rad, double p1_x, double p1_y)
0828 {
0829   // cout << endl;
0830   // cout << "calc frame" << endl;
0831   double size = 15;
0832   vector<vector<double>> crossframe(2, vector<double>(0));
0833   for (int j = 0; j < 2; j++)
0834     {
0835       for (int i = 0; i < 2; i++)
0836     {
0837       double a_ = (j == 0) ? 1 : 0;
0838       double b_ = (j == 0) ? 0 : 1;
0839       double c_ = (2 * i - 1) * size;
0840       // cout << "kesu : " << a_ << " " << b_ << " " << c_ << endl;
0841       vector<double> co_ = {a_, b_, c_};
0842       vector<double> temp_crossframe = get_crosspoint(co_, cx, cy, rad, p1_x, p1_y);
0843       if (temp_crossframe.size() == 2)
0844         {
0845           // cout << "temp_crossframe : " << temp_crossframe[0] << "  " << temp_crossframe[1] << endl;
0846           crossframe.at(0).push_back(temp_crossframe[0]);
0847           crossframe.at(1).push_back(temp_crossframe[1]);
0848         }
0849     }
0850     }
0851   // cout << "crossframe" << endl;
0852   // cout << "size : " << crossframe.size() << "  " << crossframe.at(0).size() << endl;
0853   for (int j = 0; j < crossframe.at(0).size(); j++)
0854     {
0855       for (int i = 0; i < crossframe.size(); i++)
0856     {
0857       // cout << crossframe[i][j] << ", ";
0858     }
0859       // cout << endl;
0860     }
0861   vector<double> result;
0862   if (crossframe.at(0).size() != 0)
0863     result = get_closestpoint(crossframe, p1_x, p1_y);
0864   // cout << "closest frame" << endl;
0865   for (int i = 0; i < result.size(); i++)
0866     {
0867       // cout << result[i] << ", ";
0868     }
0869   // cout << endl;
0870   // cout << "calc end" << endl;
0871 
0872   return result;
0873 }
0874 
0875 int clustEvent::GetGoodTrackNum( bool dca_range_cut,
0876                  bool is_deleted,
0877                  bool reverse )
0878 
0879 {
0880   bool condition = !reverse;
0881   int counter = 0;
0882   for (unsigned int i = 0; i < vtrack.size(); i++ )
0883     {
0884       
0885       if (vtrack[i]->dca_range_out == condition && dca_range_cut == true)
0886     {
0887       continue;
0888     }
0889       
0890       if (vtrack[i]->is_deleted == condition && is_deleted == true)
0891     {
0892       continue;
0893     }
0894       
0895       counter++;
0896     }
0897   
0898   return counter;
0899 }
0900 
0901 double clustEvent::GetTrackLeftRightAsymmetry()
0902 {
0903   int num_right = 0, num_left = 0;
0904   for( auto& track : vtrack )
0905     {
0906       double phi = track->getphi_tracklet(); // in rad
0907 
0908       //          phi = pi/2
0909       //                   |
0910       //           up-left | up-right
0911       //                   |
0912       //                   |
0913       // phi = pi ---------+--------- phi = 0
0914       //                   |
0915       //                   |
0916       //         down-left | down-right
0917       //                   |
0918       //                   phi = 3 pi / 4
0919       //
0920       // LR asymmetry = (L - R) / (L + R)
0921       // UD asymmetry = (U - D) / (U + D)
0922       //
0923       if( 0 <= phi && phi < TMath::Pi()/2 )
0924     num_right++;
0925       else if( 3.0 * TMath::Pi()/2 < phi )
0926     num_right++;
0927       else
0928     num_left++;
0929     }
0930 
0931   if( vtrack.size() == 0 )
0932     return 9999;
0933 
0934   if( num_right == 0 )
0935     return -100;
0936   else if( num_left == 0 )
0937     return 100;
0938 
0939   
0940   return 1.0 * (num_left - num_right) / (num_left + num_right);
0941 }
0942 
0943 double clustEvent::GetTrackUpDownAsymmetry()
0944 {
0945   int num_up = 0, num_down = 0;
0946   for( auto& track : vtrack )
0947     {
0948       double phi = track->getphi_tracklet(); // in rad
0949 
0950       if( 0 <= phi && phi < TMath::Pi() )
0951     num_up++;
0952       else
0953     num_down++;
0954     }
0955 
0956   if( vtrack.size() == 0 )
0957     return 9999;
0958 
0959   if( num_up == 0 )
0960     return -100;
0961   else if( num_down == 0 )
0962     return 100;
0963 
0964   
0965   return 1.0 * (num_up - num_down) / (num_up + num_down);
0966 }
0967 
0968 
0969 void clustEvent::SetTrackInfoToCluster( ) // int i, bool flag, double theta, double phi )
0970 {
0971 
0972   for( auto& track : vtrack )
0973     {
0974       for( auto& index : track->cluster_ids )
0975     {
0976       vclus[index].is_associated = true; // flag;
0977       vclus[index].track_incoming_theta = track->track_theta;
0978       vclus[index].track_incoming_phi = 0.0; // not done yet
0979         
0980     }
0981     }
0982 
0983   return;
0984 }
0985 
0986 void clustEvent::Print()
0987 {
0988   cout << "+" << string( 100, '-' ) << "+" << endl;
0989   cout << "| " << "Event: " << ievt << endl;
0990   cout << "| " << "Magnet: " << mag_on << endl;
0991   cout << "| " << "#cluster: " << vclus.size() << endl;
0992   cout << "| " << "#truth track: " << vtruth.size() << endl;
0993   cout << "| " << "#track: " << vtrack.size() << endl;
0994   cout << "| " << "DCA: "
0995        << "x = " << dca_mean[0] << ", " 
0996        << "y = " << dca_mean[1] << ", " 
0997        << "z = " << dca_mean[2] << endl;
0998   cout << "+" << string( 100, '-' ) << "+" << endl;
0999 
1000 }
[ 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!