Back to home page

sPhenix code displayed by LXR
 
 

    

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

0001 #define Analysis_cc
0002 
0003 #include "Analysis.hh"
0004 
0005 /////////////////////////////////////////////////////////////////////////
0006 // Constractor                                                         //
0007 /////////////////////////////////////////////////////////////////////////
0008 Analysis::Analysis(int run_no, int fphx_bco, bool mag_on, bool debug, bool is_preliminary  ) :
0009   run_no_( run_no ),
0010   fphx_bco_( fphx_bco ),
0011   mag_on_( mag_on ),
0012   debug_( debug ),
0013   is_preliminary_( is_preliminary )
0014 {
0015 }
0016 
0017 /////////////////////////////////////////////////////////////////////////
0018 // private functions                                                    //
0019 /////////////////////////////////////////////////////////////////////////
0020 void Analysis::Init()
0021 {
0022 
0023   //  this->InitResetVariables();
0024   this->InitIO();
0025   this->InitHist();
0026   this->InitTree();
0027 }
0028 
0029 void Analysis::InitIO()
0030 {
0031   this->InitInput();
0032   this->InitOutput();
0033 }
0034 
0035 void Analysis::InitInput()
0036 {
0037 
0038   if( fphx_bco_ == -1 )
0039     fname_input_ = data_dir_ + Form("InttAna_run%d.root", run_no_ ); // with no any cut
0040   else
0041     fname_input_ = data_dir_  + Form("InttAna_run%d_FPHX_BCO_%d.root", run_no_, fphx_bco_ ); // FPHX BCO selection
0042   
0043   // if (mag_on)
0044   //   fname = Form("AnaTutorial_inttana_%d_mag.root", run_no_); // with no any cut
0045 
0046   if (is_geant_)
0047     {
0048       if (mag_on_)
0049         {
0050       fname_input_ = "/sphenix/tg/tg01/commissioning/INTT/work/tsujibata/tracking/F4AInttRead/macro/macros/detectors/sPHENIX/AnaTutorial_inttana_zvtxwidth10_ptmin0.root";
0051         }
0052 
0053       else
0054     {
0055       fname_input_ = "/sphenix/tg/tg01/commissioning/INTT/work/tsujibata/tracking/F4AInttRead/macro/macros/detectors/sPHENIX/AnaTutorial_inttana_vtx0.root";
0056     }
0057     }
0058 
0059   cout << fname_input_ << endl;
0060   f_input_ = TFile::Open( fname_input_.c_str() );
0061 
0062   // from input file
0063   ntp_clus_ = (TNtuple *)f_input_->Get( "ntp_clus" );
0064   ntp_clus_->SetBranchAddress( "nclus"  , &nclus_   );
0065   ntp_clus_->SetBranchAddress( "nclus2" , &nclus2_  );
0066   ntp_clus_->SetBranchAddress( "bco_full", &bco_full_   );
0067   ntp_clus_->SetBranchAddress( "evt"    , &evt_     );
0068   ntp_clus_->SetBranchAddress( "size"   , &size_    );
0069   ntp_clus_->SetBranchAddress( "adc"    , &adc_     );
0070   ntp_clus_->SetBranchAddress( "x"  , &x_       );
0071   ntp_clus_->SetBranchAddress( "y"  , &y_       );
0072   ntp_clus_->SetBranchAddress( "z"  , &z_       );
0073   ntp_clus_->SetBranchAddress( "lay"    , &lay_     );
0074   ntp_clus_->SetBranchAddress( "lad"    , &lad_     );
0075   ntp_clus_->SetBranchAddress( "sen"    , &sen_     );
0076   ntp_clus_->SetBranchAddress( "x_vtx"  , &x_vtx_   );
0077   ntp_clus_->SetBranchAddress( "y_vtx"  , &y_vtx_   );
0078   ntp_clus_->SetBranchAddress( "z_vtx"  , &z_vtx_   );
0079 
0080   ntp_evt_ = (TNtuple*)f_input_->Get( "ntp_evt" );
0081   ntp_evt_->SetBranchAddress( "zv", &vertex_z_ );
0082 
0083   bco_tree_ = (TTree*)f_input_->Get( "t_evt_bco" );
0084   bco_tree_->SetBranchAddress( "bco_intt", &bco_intt_ );
0085   
0086   hepmctree_ = (TTree *)f_input_->Get( "hepmctree" );
0087   hepmctree_->SetBranchAddress( "m_evt"     , &m_evt_       );
0088   hepmctree_->SetBranchAddress( "m_status"  , &m_status_        );
0089   hepmctree_->SetBranchAddress( "m_truthpx" , &m_truthpx_       );
0090   hepmctree_->SetBranchAddress( "m_truthpy" , &m_truthpy_       );
0091   hepmctree_->SetBranchAddress( "m_truthpz" , &m_truthpz_       );
0092   hepmctree_->SetBranchAddress( "m_truthpt" , &m_truthpt_       );
0093   hepmctree_->SetBranchAddress( "m_trutheta"    , &m_trutheta_      );
0094   hepmctree_->SetBranchAddress( "m_truthpid"    , &m_truthpid_      );
0095   hepmctree_->SetBranchAddress( "m_truththeta"  , &m_truththeta_    );
0096   hepmctree_->SetBranchAddress( "m_truthphi"    , &m_truthphi_      );
0097   hepmctree_->SetBranchAddress( "m_xvtx"    , &m_xvtx_      );
0098   hepmctree_->SetBranchAddress( "m_yvtx"    , &m_yvtx_      );
0099   hepmctree_->SetBranchAddress( "m_zvtx"    , &m_zvtx_      );
0100 
0101   no_mc_ = true;
0102   if( hepmctree_ == nullptr )
0103     no_mc_ = true;
0104   else if( hepmctree_->GetEntries() == 0 )    
0105     no_mc_ = true;
0106 
0107 }
0108 
0109 void Analysis::InitOutput()
0110 {
0111   string fname_output_basename = "tracking_run" + to_string( run_no_ );
0112   if( fphx_bco_ != -1 )
0113     fname_output_basename += "_fphx_bco_" + to_string( fphx_bco_ );
0114   
0115   fname_output_ = dir_output_ + fname_output_basename;  
0116   fname_output_ += ".root";
0117   
0118   f_output_ = new TFile(fname_output_, "recreate" );
0119   
0120   c_ = new TCanvas(fname_output_.ReplaceAll(".root", ".pdf" ), "c", 1200, 600);
0121   output_pdf_ = dir_output_ + fname_output_basename + ".pdf";
0122   output_good_pdf_ = dir_output_ + fname_output_basename + "_good.pdf";
0123 
0124   this->InitCanvas();
0125   c_->cd(1);
0126   c_->Print( (output_pdf_+"[").c_str() );
0127   //  c_->Print( (output_good_pdf_+"[").c_str() );
0128 
0129 }
0130 void Analysis::InitCanvas()
0131 {
0132   c_->Clear();
0133   
0134   c_->Divide(2, 1);
0135   double margin = 0.1;
0136   auto pad1 = c_->cd( 1 );
0137   pad1->SetRightMargin( margin );
0138   pad1->SetTopMargin( margin );
0139 
0140   auto pad2 = c_->cd( 2 );
0141   pad2->SetRightMargin( margin );
0142   pad2->SetTopMargin( margin );
0143 }
0144 
0145 void Analysis::InitHist()
0146 {
0147   
0148   h_dcaz_one_       = new TH1F("h_dcaz_one", "", 100, -200, 200);
0149   h_nclus_      = new TH1F("h_nclus", "", 100, 0, 100);
0150   h_dphi_nocut_     = new TH2F("h_dphi_nocut", "d_phi : phi_in;phi_in;d_phi", 200, -3.2, 3.2, 200, -0.1, 0.1);
0151   h_dcaz_sigma_one_ = new TH1F("h_dcaz_sigma_one", "h_dcaz_sigma_one;dca_z sigma of one event;entries", 100, 0, 15);
0152   h_xvtx_       = new TH1D("h_xvtx", "h_xvtx", 100, -0.01, 0.01);
0153   h_yvtx_       = new TH1D("h_yvtx", "h_yvtx", 100, -0.01, 0.01);
0154   h_zvtx_       = new TH1D("h_zvtx", "h_zvtx", 100, -20, 20);
0155   h_dcax_       = new TH1F("h_dcax", "h_dcax", 100, -0.4, 0.4);
0156   h_dcay_       = new TH1F("h_dcay", "h_dcay", 100, -0.4, 0.4);
0157   h_dcaz_       = new TH1F("h_dcaz", "h_dcaz;DCA_z[cm]", 60, -150, 150);
0158   h_dtheta_     = new TH1F("h_dtheta", "dtheta{inner - outer};dtheta;entries", 100, -3.2, 3.2);
0159   h_dphi_       = new TH1F("h_dphi", "#Delta #phi_{AB};#Delta #phi_{AB}", 100, -1, 1);
0160   h_dca2d_      = new TH1F("h_dca2d", "h_dca", 100, -10, 10);
0161 
0162   // for (int i = 0; i < 2; i++)
0163   //   h_zr_[i] = new TH2D(Form("h_zr_%d", i), "h_zr;z;r", 100, -20, 20, 100, -10, 10);
0164 
0165 }
0166 
0167 void Analysis::InitTree()
0168 {
0169   
0170   temp_tree_ = new TTree( "temp_tree", "temp_tree" );
0171   temp_tree_->Branch( "evt_temp_"   , &evt_temp_, "evt_temp_/I"     );
0172   temp_tree_->Branch( "d_phi_"      , &d_phi_               );
0173   temp_tree_->Branch( "d_theta_"    , &d_theta_             );
0174   temp_tree_->Branch( "phi_in_"     , &phi_in_              );
0175   temp_tree_->Branch( "dca_x_"      , &dca_x_               );
0176   temp_tree_->Branch( "dca_y_"      , &dca_y_               );
0177   temp_tree_->Branch( "dca_z_"      , &dca_z_               );
0178   temp_tree_->Branch( "dca_2d_"     , &dca_2d_              );
0179   temp_tree_->Branch( "zvtx_one_"   , &zvtx_one_, "zvtx_one_/D"     );
0180   temp_tree_->Branch( "zvtx_sigma_one_" , &zvtx_sigma_one_, "zvtx_sigma_one_/D" );
0181   
0182   clus_tree_ = new TTree( "clus_tree", "clus_tree" );
0183   clus_tree_->Branch( "evt_clus"        , &evt_clus_, "evt_clus/I"  );
0184   clus_tree_->Branch( "x_in"            , &x_in_            ); // x of clusters on the inner barrel
0185   clus_tree_->Branch( "y_in"            , &y_in_            ); // y of clusters on the inner barrel
0186   clus_tree_->Branch( "z_in"            , &z_in_            ); // z of clusters on the inner barrel
0187   clus_tree_->Branch( "r_in"            , &r_in_            ); // r of clusters on the inner barrel
0188   clus_tree_->Branch( "size_in"         , &size_in_         ); // size of inner cluster
0189   clus_tree_->Branch( "phi_in"          , &phi_in_          );
0190   clus_tree_->Branch( "theta_in"            , &theta_in_            );
0191   clus_tree_->Branch( "adc_in"          , &adc_in_          );
0192   clus_tree_->Branch( "is_associated_in"    , &is_associated_in_        );
0193   clus_tree_->Branch( "track_incoming_theta_in" , &track_incoming_theta_in_ );  
0194   clus_tree_->Branch( "x_out"           , &x_out_           ); // x of clusters on the inner barrel
0195   clus_tree_->Branch( "y_out"           , &y_out_           ); // y of clusters on the inner barrel
0196   clus_tree_->Branch( "z_out"           , &z_out_           ); // z of clusters on the inner barrel
0197   clus_tree_->Branch( "r_out"           , &r_out_           ); // r of clusters on the inner barrel
0198   clus_tree_->Branch( "size_out"        , &size_out_            ); // size of outer cluster
0199   clus_tree_->Branch( "phi_out"         , &phi_out_         );
0200   clus_tree_->Branch( "theta_out"       , &theta_out_           );
0201   clus_tree_->Branch( "adc_out"         , &adc_out_         );
0202   clus_tree_->Branch( "is_associated_out"   , &is_associated_out_       );
0203   clus_tree_->Branch( "track_incoming_theta_out", &track_incoming_theta_out_    );
0204   clus_tree_->Branch( "z_vertex"        , &z_vertex_, "z_vertex/D"  );
0205 
0206   track_tree_ = new TTree( "track_tree", "track_tree" );
0207   track_tree_->Branch( "evt_track"  , &evt_track_, "evt_track/I"    );
0208   track_tree_->Branch( "ntrack"     , &ntrack_, "ntrack/I"      );
0209   track_tree_->Branch( "phi_tracklets"  , &phi_tracklets_       );
0210   track_tree_->Branch( "slope_rz"   , &slope_rz_            );
0211   track_tree_->Branch( "theta_tracklets"    , &theta_tracklets_     );
0212   track_tree_->Branch( "phi_track"  , &phi_track_           );
0213   track_tree_->Branch( "theta_track"    , &theta_track_         );
0214   track_tree_->Branch( "z_tracklets_out"    , &z_tracklets_out_     );
0215   track_tree_->Branch( "dca_2d"     , &dca_2d_          );
0216   track_tree_->Branch( "dca_z"      , &dca_z_           );
0217   track_tree_->Branch( "is_deleted" , &is_deleted_          );
0218   track_tree_->Branch( "dca_range_out"  , &dca_range_out_       );
0219   track_tree_->Branch( "dcaz_range_out" , &dcaz_range_out_      );
0220   track_tree_->Branch( "dca2d_range_out"    , &dca2d_range_out_     );
0221   track_tree_->Branch( "pT"     , &pT_              );
0222   track_tree_->Branch( "px"     , &px_              );
0223   track_tree_->Branch( "py"     , &py_              );
0224   track_tree_->Branch( "pz"     , &pz_              );
0225   track_tree_->Branch( "px_truth"   , &px_truth_            );
0226   track_tree_->Branch( "py_truth"   , &py_truth_            );
0227   track_tree_->Branch( "pz_truth"   , &pz_truth_            );
0228   track_tree_->Branch( "pT_truth"   , &pT_truth_            );
0229   track_tree_->Branch( "charge"     , &charge_          );
0230   track_tree_->Branch( "rad"        , &rad_             );
0231   // track_tree_->Branch( "x_vertex"    , &x_vertex, "x_vertex/D"   );
0232   // track_tree_->Branch( "y_vertex"    , &y_vertex, "y_vertex/D"   );
0233   // track_tree_->Branch( "z_vertex"    , &z_vertex, "z_vertex/D"   );
0234   track_tree_->Branch( "vertex" , &vertex_ ); // , "z_vertex/D" );
0235   
0236   truth_tree_ = new TTree( "truth_tree", "truth_tree" );
0237   truth_tree_->Branch( "evt_truth"  , &evt_truth_, "evt_truth/I"    );
0238   truth_tree_->Branch( "ntruth_cut" , &ntruth_cut_, "ntruth_cut/I"  );
0239   truth_tree_->Branch( "truthpx"        , &truthpx_         );
0240   truth_tree_->Branch( "truthpy"        , &truthpy_         );
0241   truth_tree_->Branch( "truthpz"        , &truthpz_         );
0242   truth_tree_->Branch( "truthpt"        , &truthpt_         );
0243   truth_tree_->Branch( "trutheta"   , &trutheta_            );
0244   truth_tree_->Branch( "truththeta" , &truththeta_          );
0245   truth_tree_->Branch( "truthphi"   , &truthphi_            );
0246   truth_tree_->Branch( "truthxvtx"  , &truthxvtx_           );
0247   truth_tree_->Branch( "truthyvtx"  , &truthyvtx_           );
0248   truth_tree_->Branch( "truthzvtx"  , &truthzvtx_           );
0249   truth_tree_->Branch( "truthzout"  , &truthzout_           );
0250   truth_tree_->Branch( "truthpid"   , &truthpid_            );
0251   truth_tree_->Branch( "status"     , &status_          );
0252 
0253 }
0254 
0255 void Analysis::ResetVariables()
0256 {
0257 
0258   this->ResetTempTreeVariables();
0259   this->ResetClustTreeVariables();
0260   this->ResetTruthTreeVariables();
0261   this->ResetTrackTreeVariables();
0262   this->ResetHepMcTreeVariables();
0263 }
0264 
0265 void Analysis::ResetTempTreeVariables()
0266 {
0267   zvtx_one_     = zvtx_sigma_one_     = -9999;
0268 
0269   dca_x_.clear();
0270   dca_y_.clear();
0271   dca_z_.clear();
0272   dca_2d_.clear();
0273   d_phi_.clear();
0274   d_theta_.clear();
0275   phi_in_.clear();
0276 }
0277 
0278 void Analysis::ResetClustTreeVariables()
0279 {
0280   x_in_.clear();
0281   y_in_.clear();
0282   z_in_.clear();
0283   r_in_.clear();
0284   size_in_.clear();
0285   phi_in_.clear();
0286   theta_in_.clear();
0287   adc_in_.clear();
0288   is_associated_in_.clear();
0289   track_incoming_theta_in_.clear();
0290      
0291   x_out_.clear();
0292   y_out_.clear();
0293   z_out_.clear();
0294   r_out_.clear();
0295   size_out_.clear();
0296   phi_out_.clear();
0297   theta_out_.clear();
0298   adc_out_.clear();
0299   is_associated_out_.clear();
0300   track_incoming_theta_out_.clear();
0301 }
0302 
0303 void Analysis::ResetTruthTreeVariables()
0304 {
0305 
0306   evt_truth_ = ntruth_cut_ = -9999;
0307   truthpx_.clear();
0308   truthpy_.clear();
0309   truthpz_.clear();
0310   truthpt_.clear();
0311   trutheta_.clear();
0312   truththeta_.clear();
0313   truthphi_.clear();
0314   truthxvtx_.clear();
0315   truthyvtx_.clear();
0316   truthzvtx_.clear();
0317   truthzout_.clear();
0318   truthpid_.clear();
0319   status_.clear();
0320 }
0321 
0322 void Analysis::ResetTrackTreeVariables() // doesn't work at all!
0323 {
0324   evt_track_ =  ntrack_ = -9999;
0325   
0326   phi_tracklets_.clear();
0327   slope_rz_.clear();
0328   theta_tracklets_.clear();
0329   phi_track_.clear();
0330   theta_track_.clear();
0331   z_tracklets_out_.clear();
0332   dca_2d_.clear();
0333   dca_z_.clear();
0334   is_deleted_.clear();
0335   dca_range_out_.clear();
0336   dcaz_range_out_.clear();
0337   dca2d_range_out_.clear();
0338   pT_.clear();
0339   px_.clear();
0340   py_.clear();
0341   pz_.clear();
0342   px_truth_.clear();
0343   py_truth_.clear();
0344   pz_truth_.clear();
0345   pT_truth_.clear();
0346   charge_.clear();
0347   rad_.clear();
0348   vertex_ = TVector3();
0349 
0350 }
0351 
0352 void Analysis::ResetHepMcTreeVariables()
0353 {
0354   m_evt_ =  m_status_ =   m_truthpid_ = -9999;
0355   
0356   m_truthpx_ = m_truthpy_ = m_truthpz_ = m_truthpt_
0357     = m_trutheta_ = m_truththeta_ = m_truthphi_
0358     = m_xvtx_ = m_yvtx_ = m_zvtx_
0359     = -9999;
0360 
0361 }
0362 
0363 void Analysis::ProcessEvent()
0364 {
0365   
0366   int sum_event = 0;
0367   int ntrack_ = 0;
0368   for (int icls = 0; icls < ntp_clus_->GetEntries(); icls++)
0369     {
0370       cout << Form("  -----  event %d {event gene.}  -----  ", sum_event) << endl;
0371 
0372       clustEvent event{};
0373       cluster clust{};
0374 
0375       ntp_clus_->GetEntry(icls); // first cluster in one event
0376 
0377       clust.set(evt_, 0, x_, y_, z_, adc_, size_, lay_, x_vtx_, y_vtx_, z_vtx_ );
0378       event.vclus.push_back(clust);
0379 
0380       for (int j = 1; j < nclus_; j++) // second cluster in one event
0381         {
0382       ntp_clus_->GetEntry(++icls);
0383       cluster clust2{};
0384       clust2.set(evt_, 0, x_, y_, z_, adc_, size_, lay_, x_vtx_, y_vtx_, z_vtx_ );
0385       event.vclus.push_back(clust2);
0386         }
0387 
0388       dotracking(event, h_dphi_nocut_ );
0389       
0390       ntrack_ = event.vtrack.size();
0391       h_dcaz_one_->Reset();
0392       for (int itrk = 0; itrk < ntrack_; itrk++)
0393         {
0394       h_dcaz_one_->Fill(event.vtrack[itrk]->dca[2]);
0395         }
0396       zvtx_one_       = h_dcaz_one_->GetMean();
0397       zvtx_sigma_one_ = h_dcaz_one_->GetStdDev();
0398 
0399       evt_temp_ = sum_event;
0400       for (int itrk = 0; itrk < ntrack_; itrk++)
0401         {
0402       dca_x_.push_back( event.vtrack[itrk]->dca[0] );
0403       dca_y_.push_back( event.vtrack[itrk]->dca[1] );
0404       dca_z_.push_back( event.vtrack[itrk]->dca[2] );
0405       dca_2d_.push_back( event.vtrack[itrk]->dca_2d );
0406       d_phi_.push_back( event.vtrack[itrk]->d_phi );
0407       d_theta_.push_back( event.vtrack[itrk]->d_theta );
0408       phi_in_.push_back( event.vtrack[itrk]->phi_in );
0409         }
0410       
0411       temp_tree_->Fill();      
0412       this->ResetTempTreeVariables();
0413       
0414       event.clear();
0415       sum_event++;
0416       
0417       if (debug_ && sum_event > nloop_)
0418     break;
0419       
0420     } // end of for (int icls = 0; icls < ntp_clus_->GetEntries(); icls++)
0421   
0422   n_dotracking++;
0423   
0424   temp_tree_->Draw( "dca_x_>>h_dcax" );
0425   temp_tree_->Draw( "dca_y_>>h_dcay" );
0426   temp_tree_->Draw( "dca_z_>>h_dcaz" );
0427   temp_tree_->Draw( "d_theta_>>h_dtheta" );
0428   temp_tree_->Draw( "d_phi_>>h_dphi" );
0429   temp_tree_->Draw( "dca_2d_>>h_dca2d" );
0430 
0431   //vector < TH1F* > h_ = {h_dcax, h_dcay, h_dcaz, h_dphi, h_dtheta, h_dca2d};
0432 
0433   // the mean of DCA in all event
0434   double dcax_mean = h_dcax_->GetMean();
0435   double dcay_mean = h_dcay_->GetMean();
0436   if (!(is_geant_))
0437     {
0438       dcax_mean = -0.019;
0439       dcay_mean = 0.198;
0440     }
0441   double dcaz_mean = h_dcaz_->GetMean();
0442   double dca2d_mean = h_dca2d_->GetMean();
0443 
0444   double dcax_sigma = h_dcax_->GetStdDev();
0445   double dcay_sigma = h_dcay_->GetStdDev();
0446   double dcaz_sigma = h_dcaz_->GetStdDev();
0447   double dca2d_sigma = h_dca2d_->GetStdDev();
0448 
0449   dca2d_sigma *= sigma_;
0450   dcaz_sigma *= sigma_;
0451 
0452   double dca2d_max = dca2d_mean + dca2d_sigma;
0453   double dca2d_min = dca2d_mean - dca2d_sigma;
0454   double dcaz_max = dcaz_mean + dcaz_sigma;
0455   double dcaz_min = dcaz_mean - dcaz_sigma;
0456 
0457   // tracking
0458   int ihit = 0;
0459   int itruth = 0;
0460   int temp_evt = 0;
0461   for (int ievt = 0; ievt < sum_event; ievt++, ihit++, itruth++)
0462   //  for (int ievt = 0; ievt < 3; ievt++, ihit++, itruth++)
0463     {
0464       temp_tree_->GetEntry(ievt);
0465       ntp_evt_->GetEntry( ievt );
0466       bco_tree_->GetEntry( ievt );
0467       
0468       cout << "\r" << "event " << setw(10) << right << ievt;
0469       ntruth_ = 0;
0470       ntruth_cut_ = 0;
0471     
0472       clustEvent event;
0473       event.ievt = ievt;
0474       event.mag_on = mag_on_;
0475       event.run_no = run_no_;
0476       event.bco_intt = bco_intt_;
0477 
0478       if(ihit < ntp_clus_->GetEntries() )
0479         {
0480       ntp_clus_->GetEntry(ihit);
0481 
0482       if (!( no_mc_ ))
0483             {
0484           hepmctree_->GetEntry(itruth);
0485           while ( m_evt_ != evt_ )
0486                 {
0487           itruth++;
0488           hepmctree_->GetEntry(itruth);
0489                 }
0490           temp_evt = m_evt_;
0491             }
0492 
0493       cluster clust{};
0494       clust.set(evt_, 0, x_, y_, z_, adc_, size_, lay_, x_vtx_, y_vtx_, z_vtx_ );
0495       event.vclus.push_back(clust);
0496 
0497       for (int j = 1; j < nclus_; j++) // second~ cluster in one event
0498             {
0499           ntp_clus_->GetEntry(++ihit);
0500           cluster clust2{};
0501 
0502           clust2.set(evt_, 0, x_, y_, z_, adc_, size_, lay_, x_vtx_, y_vtx_, z_vtx_ );
0503           event.vclus.push_back(clust2);
0504             }
0505 
0506       if (!( no_mc_ ))
0507             {
0508           while (m_evt_ == temp_evt)
0509                 {
0510           ntruth_++;
0511           truth *tru = new truth();
0512           tru->set_truth(m_truthpx_, m_truthpy_, m_truthpz_, m_truthpt_,
0513                  m_status_, m_trutheta_, m_truthpid_,
0514                  m_truththeta_, m_truthphi_,
0515                  m_xvtx_, m_yvtx_, m_zvtx_);
0516           
0517           event.vtruth.push_back(tru);
0518           
0519           itruth++;
0520           if (itruth == hepmctree_->GetEntries())
0521             break;
0522 
0523           hepmctree_->GetEntry(itruth);
0524                 }
0525           itruth--;
0526             }
0527         }
0528 
0529       event.dca_mean[0] = dcax_mean;
0530       event.dca_mean[1] = dcay_mean;
0531       event.dca_mean[2] = zvtx_one_;
0532 
0533       // from dca_mean
0534       dotracking(event, h_dphi_nocut_ );
0535       // event.Print();
0536       // cout << zvtx_one_ << "\t" << z_vtx_ << endl;
0537       
0538       ntrack_ = event.vtrack.size();
0539       evt_clus_ = ievt;
0540       evt_track_ = ievt;
0541       evt_truth_ = ievt;
0542       h_dcaz_one_->Reset();
0543       for (int itrk = 0; itrk < ntrack_; itrk++) // from dca_mean
0544         {
0545       h_dcaz_one_->Fill(event.vtrack[itrk]->dca[2]);
0546         }
0547       
0548       double dcaz_mean_one = h_dcaz_one_->GetMean();    // one event
0549 
0550       // if( fabs(vertex_z_) < 300 )
0551       //    cout << dcaz_mean << "\t"
0552       //         << dcaz_mean_one << "\t"
0553       //         << vertex_z_ << endl;
0554       
0555       double dcaz_sigma_one = h_dcaz_one_->GetStdDev(); // one event
0556       h_dcaz_sigma_one_->Fill(dcaz_sigma_one);
0557       
0558       dcaz_sigma_one *= sigma_;
0559       double dcaz_max_one = dcaz_mean_one + dcaz_sigma;
0560       double dcaz_min_one = dcaz_mean_one - dcaz_sigma;
0561       zvtx_sigma_one_ *= sigma_;
0562       // double dcaz_max_one = zvtx_one_ + zvtx_sigma_one_;
0563       // double dcaz_min_one = zvtx_one_ - zvtx_sigma_one_;
0564 
0565       event.dca2d_max = dca2d_max;
0566       event.dca2d_min = dca2d_min;
0567       event.dcaz_max = dcaz_max_one;
0568       event.dcaz_min = dcaz_min_one;
0569 
0570       for (int itrk = 0; itrk < ntrack_; itrk++)
0571         {
0572       event.vtrack[itrk]->dca_mean[0] = dcax_mean;
0573       event.vtrack[itrk]->dca_mean[1] = dcay_mean;
0574       event.vtrack[itrk]->dca_mean[2] = dcaz_mean_one;
0575       event.vtrack[itrk]->calc_line();
0576       event.vtrack[itrk]->calc_inv();
0577       event.vtrack[itrk]->calc_pT();
0578 
0579         }
0580       
0581       event.SetTrackInfoToCluster(); // set track information to associated clusters
0582 
0583       if (!( no_mc_ ))
0584         {
0585       for (int itrt = 0; itrt < ntruth_; itrt++)
0586             {
0587           event.vtruth[itrt]->dca_mean[0] = dcax_mean;
0588           event.vtruth[itrt]->dca_mean[1] = dcay_mean;
0589           event.vtruth[itrt]->dca_mean[2] = dcaz_mean_one;
0590 
0591           event.vtruth[itrt]->calc_line();
0592           event.vtruth[itrt]->calc_center();
0593             }
0594         }
0595       h_nclus_->Fill(event.vclus.size());
0596 
0597       for (int iclus = 0; iclus < event.vclus.size(); iclus++)
0598         {
0599       event.vclus[iclus].dca_mean[0] = dcax_mean;
0600       event.vclus[iclus].dca_mean[1] = dcay_mean;
0601       event.vclus[iclus].dca_mean[2] = dcaz_mean_one;
0602       z_vertex_ = dcaz_mean_one;
0603 
0604       if (event.vclus[iclus].layer < 2)
0605             {
0606           x_in_.push_back               ( event.vclus[iclus].x              );
0607           y_in_.push_back               ( event.vclus[iclus].y              );
0608           z_in_.push_back               ( event.vclus[iclus].z              );
0609           r_in_.push_back               ( event.vclus[iclus].r              );
0610           size_in_.push_back            ( event.vclus[iclus].size           );
0611           phi_in_.push_back             ( event.vclus[iclus].getphi_clus()      );
0612           theta_in_.push_back           ( event.vclus[iclus].gettheta_clus()        );
0613           adc_in_.push_back             ( event.vclus[iclus].adc            );
0614           is_associated_in_.push_back       ( event.vclus[iclus].is_associated      );
0615           track_incoming_theta_in_.push_back    ( event.vclus[iclus].track_incoming_theta   );
0616           //event.vclus[iclus]_.print(                              );
0617           
0618             }
0619       if (1 < event.vclus[iclus].layer)
0620             {
0621           x_out_.push_back              ( event.vclus[iclus].x              );
0622           y_out_.push_back              ( event.vclus[iclus].y              );
0623           z_out_.push_back              ( event.vclus[iclus].z              );
0624           r_out_.push_back              ( event.vclus[iclus].r              );
0625           size_out_.push_back           ( event.vclus[iclus].size           );
0626           phi_out_.push_back            ( event.vclus[iclus].getphi_clus()      );
0627           theta_out_.push_back          ( event.vclus[iclus].gettheta_clus()        );
0628           adc_out_.push_back            ( event.vclus[iclus].adc            );
0629           is_associated_out_.push_back      ( event.vclus[iclus].is_associated      );
0630           track_incoming_theta_out_.push_back   ( event.vclus[iclus].track_incoming_theta   );
0631             }
0632         }
0633 
0634       clus_tree_->Fill();
0635       this->ResetClustTreeVariables();
0636 
0637       /*
0638       for (int itrk = 0; itrk < ntrack_; itrk++)
0639         {
0640       event.vtrack[itrk]->dca_mean[0] = dcax_mean;
0641       event.vtrack[itrk]->dca_mean[1] = dcay_mean;
0642       event.vtrack[itrk]->dca_mean[2] = dcaz_mean_one;
0643       event.vtrack[itrk]->calc_line();
0644       event.vtrack[itrk]->calc_inv();
0645       event.vtrack[itrk]->calc_pT();
0646 
0647       event.vtrack[itrk]->print();    
0648         }
0649       
0650       event.SetTrackInfoToCluster();
0651       */
0652       
0653       event.dca_check(debug_);
0654       event.makeonetrack();
0655 
0656       event.charge_check();
0657       event.truth_check();
0658 
0659       event.cluster_check();
0660       event.track_check();
0661       ntrack_ = 0;
0662       for (int itrk = 0; itrk < ntrack_; itrk++)
0663         {
0664       slope_rz_.push_back       (event.vtrack[itrk]->track_rz->GetParameter(1));
0665       phi_tracklets_.push_back  (event.vtrack[itrk]->getphi_tracklet());
0666       theta_tracklets_.push_back    (event.vtrack[itrk]->gettheta_tracklet());
0667       phi_track_.push_back      (event.vtrack[itrk]->getphi());
0668       theta_track_.push_back        (event.vtrack[itrk]->gettheta());
0669       dca_z_.push_back      (event.vtrack[itrk]->dca[2]);
0670       dca_2d_.push_back     (event.vtrack[itrk]->dca_2d);
0671       z_tracklets_out_.push_back    (event.vtrack[itrk]->p2.z());
0672       pT_.push_back         (event.vtrack[itrk]->pT);
0673       px_.push_back         (event.vtrack[itrk]->p_reco[0]);
0674       py_.push_back         (event.vtrack[itrk]->p_reco[1]);
0675       pz_.push_back         (event.vtrack[itrk]->p_reco[2]);
0676       rad_.push_back            (event.vtrack[itrk]->rad);
0677       charge_.push_back     (event.vtrack[itrk]->charge);
0678       if (!( no_mc_ ))
0679             {
0680           pT_truth_.push_back(event.vtruth[0]->m_truthpt);
0681           px_truth_.push_back(event.vtruth[0]->m_truthpx);
0682           py_truth_.push_back(event.vtruth[0]->m_truthpy);
0683           pz_truth_.push_back(event.vtruth[0]->m_truthpz);
0684             }
0685 
0686       is_deleted_.push_back(event.vtrack[itrk]->is_deleted);
0687       dca2d_range_out_.push_back(event.vtrack[itrk]->dca2d_range_out);
0688       dcaz_range_out_.push_back(event.vtrack[itrk]->dcaz_range_out);
0689       dca_range_out_.push_back(event.vtrack[itrk]->dca_range_out);
0690 
0691       // x_vertex_ = event.vtrack[itrk]->dca_mean[0];
0692       // y_vertex_ = event.vtrack[itrk]->dca_mean[1];
0693       z_vertex_ = event.vtrack[itrk]->dca_mean[2];
0694       vertex_.SetXYZ( event.vtrack[itrk]->dca_mean[0],
0695               event.vtrack[itrk]->dca_mean[1],
0696               event.vtrack[itrk]->dca_mean[2] );
0697 
0698       if (event.vtrack[itrk]->is_deleted == true)
0699         continue;
0700       if (event.vtrack[itrk]->dca_range_out == true)
0701         continue;
0702 
0703       ntrack_++;
0704 
0705       // h_xvtx->Fill(x_vertex_);
0706       // h_yvtx->Fill(y_vertex_);
0707       // h_zvtx->Fill(z_vertex_);
0708 
0709       h_xvtx_->Fill( vertex_.X() );
0710       h_yvtx_->Fill( vertex_.Y() );
0711       h_zvtx_->Fill( vertex_.Z() );
0712         }
0713 
0714       track_tree_->Fill();
0715       double temp = z_vertex_;
0716       this->ResetTrackTreeVariables();
0717       z_vertex_ = temp;
0718       
0719       if ( !( no_mc_ ) )
0720         {
0721       for (int itrt = 0; itrt < ntruth_; itrt++)
0722             {
0723 
0724           if (event.vtruth[itrt]->is_charged == false)
0725         continue;
0726 
0727           if (event.vtruth[itrt]->is_intt == false)
0728         continue;
0729           
0730           ntruth_cut_++;
0731 
0732           truthpx_.push_back    ( event.vtruth[itrt]->m_truthpx     );
0733           truthpy_.push_back    ( event.vtruth[itrt]->m_truthpy     );
0734           truthpz_.push_back    ( event.vtruth[itrt]->m_truthpz     );
0735           truthpt_.push_back    ( event.vtruth[itrt]->m_truthpt     );
0736           trutheta_.push_back   ( event.vtruth[itrt]->m_trutheta    );
0737           truththeta_.push_back ( event.vtruth[itrt]->m_truththeta  );
0738           truthphi_.push_back   ( event.vtruth[itrt]->m_truthphi    );
0739           truthpid_.push_back   ( event.vtruth[itrt]->m_truthpid    );
0740           status_.push_back     ( event.vtruth[itrt]->m_status      );
0741           truthxvtx_.push_back  ( event.vtruth[itrt]->m_xvtx        );
0742           truthyvtx_.push_back  ( event.vtruth[itrt]->m_yvtx        );
0743           truthzvtx_.push_back  ( event.vtruth[itrt]->m_zvtx        );
0744           truthzout_.push_back  ( event.vtruth[itrt]->m_truthz_out  );
0745             }
0746 
0747       truth_tree_->Fill();
0748       this->ResetTruthTreeVariables();
0749 
0750         }
0751 
0752       //      if( event.ievt > 39600 )
0753       this->ProcessEvent_Draw( event );
0754       
0755       event.clear();
0756       event_counter_++;
0757 
0758       this->ResetVariables(); // variables need to be reset at the end of event
0759       // if( ievt >20 )
0760       //    break;
0761       
0762       // if( page_counter_ > page_num_limit_ )
0763       //    break;
0764       
0765     } // End of  for (int ievt = 0; ievt < sum_event; ievt++, ihit++, itruth++)
0766 
0767 }
0768 
0769 
0770 void Analysis::ProcessEvent_Draw( clustEvent event )
0771 {
0772   // drawing
0773   // if (fabs(event.vclus.size() - 2 * ntrack) < 5 && /*debug_ &&*/ 10 < event.vclus.size() && event.vclus.size() < 50 /* && && is_geant_ &&*/ /*ievt < 500*/)
0774   // if (ievt < nloop && ievt == 56 /*&& ievt==56*/ /*&& fabs(event.dca_mean[2]) < 2. && ntrack < 20*/)
0775   
0776   {
0777     if( page_counter_ > page_num_limit_ )
0778       return;
0779 
0780     double vertex_range_z = 15;
0781     if( event.dca_mean[2] < -vertex_range_z || event.dca_mean[2] > vertex_range_z )
0782       return;
0783 
0784     int track_num_min = 5, track_num_max = 200;
0785     int cluster_num_tolerance = 5;
0786     if( event.vclus.size() < track_num_min * 2 || event.vclus.size() > track_num_max * 2 + cluster_num_tolerance )
0787       return;
0788     
0789     int good_track_num = event.GetGoodTrackNum( true, true, false ); // dca_range_cut, is_deleted, reverse
0790     if( good_track_num < track_num_min || track_num_max < good_track_num )
0791       return;
0792 
0793     double ratio = event.GetAssociatedClusterRatio();
0794     //    if( ratio < 0.8 || 0.95 < ratio )
0795     if( ratio < 0.7 || 0.95 < ratio )
0796       return;
0797 
0798     // double up_down_asymmetry = event.GetTrackUpDownAsymmetry();
0799     // if( up_down_asymmetry < -0.25 || 0.25 < up_down_asymmetry )
0800     //   return;
0801     
0802     // double left_right_asymmetry = event.GetTrackLeftRightAsymmetry();
0803     // // if( left_right_asymmetry < -0.25 || 0.25 < left_right_asymmetry )
0804     // //   return;
0805     
0806     cout << endl
0807      << "Event " << setw(10) << event_counter_ + 1 << " "
0808      << "Page " << setw(3) << page_counter_ + 1 << " "
0809      << event.vclus.size() << "\t"
0810      << good_track_num << "\t"
0811      << ratio << "\t"
0812      << event.dca_mean[2] << "\t"
0813      // << event.GetTrackUpDownAsymmetry() << "\t"
0814      // << event.GetTrackLeftRightAsymmetry()
0815      << endl;
0816     
0817     this->InitCanvas();
0818     c_->cd(1);
0819 
0820     event.draw_frame(0, is_preliminary_);
0821     event.draw_intt(0);
0822     event.draw_clusters(0, true, kBlack);
0823     
0824     if (mag_on_)
0825       {
0826     // if (!( no_mc_ ))
0827     //     event.draw_truthcurve(0, true, kGray + 1);
0828     // event.draw_trackcurve(0, true, TColor::GetColor(88, 171, 141), true, true, false);
0829     event.draw_trackcurve(0, true, true, true, false);
0830       }
0831     else
0832       {
0833     // if (!( no_mc_ ))
0834     //     event.draw_truthline(0, true, kGray + 1);
0835     // TColor::GetColor(232, 85, 98), // TColor::GetColor(88, 171, 141),
0836     event.draw_trackline(0, true, 
0837                  true, true, false);
0838     
0839     // event.draw_truthline(0, true, TColor::GetColor(232, 85, 98));
0840       }
0841     
0842     // event.draw_tracklets(0, true, kGray + 1, true, true);
0843 
0844     c_->cd(2);
0845     event.draw_frame(1, is_preliminary_);
0846     event.draw_intt(1);
0847 
0848     // event.draw_tracklets(1, false, kGray + 1, false, false);
0849     // event.draw_tracklets(1, true, TColor::GetColor(0, 118, 186), true, true);
0850     event.draw_clusters(1, true, kBlack);
0851     
0852     // if (!( no_mc_ ))
0853     //     event.draw_truthline(1, true, kGray + 1);
0854     event.draw_trackline(1, true, true, true, false);    
0855     
0856     if( page_counter_ < page_num_limit_ )
0857       {
0858     c_->Print( output_pdf_.c_str() );
0859 
0860     string output_individual = output_pdf_.substr( 0, output_pdf_.size() - 4 );
0861     output_individual += "_bco_" + to_string( event.bco_intt );
0862 
0863     if( is_preliminary_ == true )
0864       output_individual += "_preliminary.pdf";
0865     else
0866       output_individual += "_internal.pdf";
0867 
0868     cout << "individual output : " << output_individual << endl;
0869     c_->Print( output_individual.c_str() );
0870       }
0871 
0872     string canvas_name = string("canvas_") + to_string( event.ievt );
0873     f_output_->WriteTObject( c_, canvas_name.c_str() );
0874 
0875     page_counter_++;
0876   }
0877   
0878 }
0879 
0880 void Analysis::EndProcess()
0881 {
0882   // c_->Print(c_->GetName());
0883   c_->Print( (output_pdf_+"]").c_str() );
0884   //  c_->Print( (output_good_pdf_+"]").c_str() );
0885   f_output_->Write();
0886 
0887   cout << endl
0888        << page_counter_ << " pages of "
0889        << event_counter_ << " events "
0890        << "(" << 100.0 * page_counter_ / event_counter_ << "%) "
0891        << "were saved." << endl;
0892   
0893   vector < TH1F* > h_ = {h_dcax_, h_dcay_, h_dcaz_, h_dphi_, h_dtheta_, h_dca2d_};
0894   for (int ihst = 0; ihst < h_.size(); ihst++)
0895     {
0896       f_output_->WriteTObject(h_[ihst], h_[ihst]->GetName());
0897     }
0898 
0899   f_output_->Close();
0900   f_input_->Close();
0901 }
0902 
0903 /////////////////////////////////////////////////////////////////////////
0904 // public functions                                                    //
0905 /////////////////////////////////////////////////////////////////////////
0906 
0907 void Analysis::Begin()
0908 {
0909 
0910   this->Init();
0911   this->ProcessEvent();
0912   this->EndProcess();
0913 }
[ 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!