event_display/macro/analysis_tracking.cc

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