analysis_INTT/src/BeamspotReco.cxx

0001 #include "TF1.h"
0002 #include "TF2.h"
0003 #include "TFile.h"
0004 #include "TGraphErrors.h"
0005 #include "TH1F.h"
0006 #include "TProfile.h"
0007 #include "TTree.h"
0008
0009 #include <Eigen/Dense>
0010
0011 #include "Beamspot.h"
0012 #include "Hit.h"
0013 #include "Utilities.h"
0014
0015 Beamspot fit_PCAD0Phi0(std::vector<std::pair<std::shared_ptr<Hit>, std::shared_ptr<Hit>>> tracklets, Beamspot BS_init, double d0_cut)
0016 {
0017     Beamspot reco_BS;
0018
0019     int Ntracklets = 0;
0020
0021     std::cout << "processing " << tracklets.size() << " tracklets..." << std::endl;
0022
0023     for (auto &tracklet : tracklets)
0024     {
0025         double x1 = tracklet.first->posX();
0026         double y1 = tracklet.first->posY();
0027         double z1 = tracklet.first->posZ();
0028         double x2 = tracklet.second->posX();
0029         double y2 = tracklet.second->posY();
0030         double z2 = tracklet.second->posZ();
0031
0032         int size1 = tracklet.first->PhiSize();
0033         int size2 = tracklet.second->PhiSize();
0034
0035         // std::cout << "tracklet: (" << x1 << ", " << y1 << ", " << z1 <<") -> (" << x2 << ", " << y2 << ", " << z2 << ")" << std::endl;
0036
0037         // cluster position uncertainties
0038         double stripwidth = 0.0078; // cm
0039         double rphierr1 = stripwidth * size1;
0040         double rphierr2 = stripwidth * size2;
0041         double r1 = sqrt(x1 * x1 + y1 * y1);
0042         double r2 = sqrt(x2 * x2 + y2 * y2);
0043         double phi1 = atan2(y1, x1);
0044         double phi2 = atan2(y2, x2);
0045         double sigma2_x1 = pow(rphierr1 * cos(phi1), 2.);
0046         double sigma2_y1 = pow(rphierr1 * sin(phi1), 2.);
0047         double sigma2_x2 = pow(rphierr2 * cos(phi2), 2.);
0048         double sigma2_y2 = pow(rphierr2 * sin(phi2), 2.);
0049
0050         // phi0 (angle of momentum vector at PCA) and uncertainty
0051         double px = x2 - x1;
0052         double py = y2 - y1;
0053         double sigma2_px = sigma2_x2 + sigma2_x1;
0054         double sigma2_py = sigma2_y2 + sigma2_y1;
0055         double phi0 = atan2(py, px);
0056         double sigma2_phi0 = (sigma2_px * py * py + sigma2_py * px * px) / pow(px * px + py * py, 2.); // the aftermath of a bit of algebra
0057         // std::cout << "phi0: " << phi0 << " +- " << sqrt(sigma2_phi0) << std::endl;
0058
0059         // PCA and uncertainty
0060         double pca_x = x1 + ((BS_init.x - x1) * cos(phi0) + (BS_init.y - y1) * sin(phi0)) * cos(phi0);
0061         double pca_y = y1 + ((BS_init.x - x1) * cos(phi0) + (BS_init.y - y1) * sin(phi0)) * sin(phi0);
0062         // std::cout << "PCA: (" << pca_x << " +- " << sqrt(sigma2_pca_x) << ", " << pca_y << " +- " << sqrt(sigma2_pca_y) << std::endl;
0063
0064         double dca_origin = sqrt(pca_x * pca_x + pca_y * pca_y);
0065         double dzdr = (z2 - z1) / (r2 - r1);
0066         double z0 = z1 - dzdr * (r1 - dca_origin);
0067
0068         double d0 = sqrt(pow(pca_x - BS_init.x, 2.) + pow(pca_y - BS_init.y, 2.));
0069         if (d0 > d0_cut)
0070             continue;
0071         double phi_pca = atan2(pca_y - BS_init.y, pca_x - BS_init.x);
0072         double oppositephi_pca = phi_pca + M_PI;
0073         if (oppositephi_pca > M_PI)
0074             oppositephi_pca -= 2 * M_PI;
0075
0076         reco_BS.d0phi0dist->Fill(phi_pca, d0);
0077         reco_BS.d0phi0dist->Fill(oppositephi_pca, -d0);
0078
0079         reco_BS.pcadist->Fill(pca_x, pca_y);
0080         reco_BS.z0dist->Fill(z0);
0081
0082         Ntracklets++;
0083     }
0084
0085     TF1 z_f("gaus", "gaus", -60., 60.);
0086
0087     // Remove "background"
0088     reco_BS.d0phi0dist_bkgrm = (TH2F *)reco_BS.d0phi0dist->Clone("d0phi0dist_bkgrm");
0089     for (int i = 0; i < reco_BS.d0phi0dist->GetNbinsX(); i++)
0090     {
0091         TH1D *slice = reco_BS.d0phi0dist->ProjectionY("slice", i, i + 1);
0092         int maxbin = slice->GetMaximumBin();
0093         float max = slice->GetBinContent(maxbin);
0094         for (int j = 0; j < slice->GetNbinsX(); j++)
0095         {
0096             float content = slice->GetBinContent(j + 1);
0097             if (slice->GetBinContent(j + 1) < 0.995 * max)
0098                 reco_BS.d0phi0dist_bkgrm->SetBinContent(i + 1, j + 1, 0.);
0099         }
0100     }
0101
0102     reco_BS.d0phi0dist_bkgrm_prof = reco_BS.d0phi0dist_bkgrm->ProfileX("d0phi0dist_bkgrm_prof");
0103
0104     TGraphErrors *g = new TGraphErrors();
0105     for (int i = 0; i < reco_BS.d0phi0dist_bkgrm_prof->GetNbinsX(); i++)
0106     {
0107         double d0 = reco_BS.d0phi0dist_bkgrm_prof->GetBinContent(i + 1);
0108         double phi = -M_PI + i * (2. * M_PI) / 200.;
0109         g->SetPoint(i, phi, d0);
0110     }
0111
0112     TF1 d0cos("d0cos", "[0]*TMath::Cos(x-[1])", -M_PI, M_PI);
0113     g->Fit(&d0cos);
0114     g->Write("maxgraph_prermoutl");
0115
0116     vector<float> phi_2;
0117     vector<float> d0_2;
0118     TF1 d0cos_rmoutl("d0cos_rmoutl", "[0]*TMath::Cos(x-[1])", -M_PI, M_PI);
0119     for (int i = 0; i < g->GetN(); i++)
0120     {
0121         double x, y;
0122         g->GetPoint(i, x, y);
0123         double y_fit = d0cos.Eval(x);
0124         // cout << "x: " << x << " y: " << y << " y_fit: " << y_fit << "abs(y-y_fit): " << fabs(y - y_fit) << endl;
0125         if (fabs(y - y_fit) < 0.2)
0126         {
0127             phi_2.push_back(x);
0128             d0_2.push_back(y);
0129         }
0130     }
0131
0132     TGraphErrors *g2 = new TGraphErrors(phi_2.size(), &phi_2[0], &d0_2[0]);
0133     // fit the graph again
0134     g2->Fit(&d0cos_rmoutl);
0135     g2->Write("maxgraph_rmoutl");
0136
0137     reco_BS.z0dist->Fit(&z_f);
0138
0139     // build reconstructed beamspot: x and y first
0140     double BS_d0 = d0cos_rmoutl.GetParameter(0);
0141     double BS_phi0 = d0cos_rmoutl.GetParameter(1);
0142     double BS_sigmad0 = d0cos_rmoutl.GetParError(0);
0143     double BS_sigmaphi0 = d0cos_rmoutl.GetParError(1);
0144
0145     reco_BS.x = BS_d0 * cos(BS_phi0) + BS_init.x;
0146     reco_BS.y = BS_d0 * sin(BS_phi0) + BS_init.y;
0147     reco_BS.sigma_x = sqrt(pow(cos(BS_phi0), 2.) * pow(BS_sigmad0, 2.) + pow(BS_d0 * sin(BS_phi0), 2.) * pow(BS_sigmaphi0, 2.));
0148     reco_BS.sigma_y = sqrt(pow(sin(BS_phi0), 2.) * pow(BS_sigmad0, 2.) + pow(BS_d0 * cos(BS_phi0), 2.) * pow(BS_sigmaphi0, 2.));
0149
0150     reco_BS.Ntracklets = Ntracklets;
0151
0152     reco_BS.z = z_f.GetParameter("Mean");
0153     reco_BS.sigma_z = z_f.GetParameter("Sigma");
0154
0155     return reco_BS;
0156 }
0157
0158 int main(int argc, char **argv)
0159 {
0160     if (argc != 4)
0161     {
0162         std::cout << "Usage: ./BeamspotReco [infile] [outfile] [dphi_cut]" << std::endl;
0163         return 0;
0164     }
0165
0166     std::string infile = argv[1];
0167     std::string outfile = argv[2];
0168     double dphi_cut = atof(argv[3]);
0169
0170     TFile *inf = TFile::Open(infile.c_str());
0171     TTree *events = (TTree *)inf->Get("EventTree");
0172
0173     TFile *outf = new TFile(outfile.c_str(), "RECREATE");
0174     TTree *t = new TTree("reco_beamspot", "reco_beamspot");
0175
0176     uint64_t intt_bco;
0177     int Nclusters;
0178     std::vector<double> *clusters_x = 0;
0179     std::vector<double> *clusters_y = 0;
0180     std::vector<double> *clusters_z = 0;
0181     std::vector<int> *clusters_layer = 0;
0182     std::vector<int> *clusters_phisize = 0;
0183     std::vector<unsigned int> *clusters_adc = 0;
0184     bool InttBco_IsToBeRemoved;
0185
0186     events->SetBranchAddress("INTT_BCO", &intt_bco);
0187     if (events->GetListOfBranches()->FindObject("InttBco_IsToBeRemoved"))
0188     {
0189         events->SetBranchAddress("InttBco_IsToBeRemoved", &InttBco_IsToBeRemoved);
0190     }
0191     else
0192     {
0193         InttBco_IsToBeRemoved = false;
0194     }
0195     events->SetBranchAddress("NClus", &Nclusters);
0196     events->SetBranchAddress("ClusX", &clusters_x);
0197     events->SetBranchAddress("ClusY", &clusters_y);
0198     events->SetBranchAddress("ClusZ", &clusters_z);
0199     events->SetBranchAddress("ClusLayer", &clusters_layer);
0200     events->SetBranchAddress("ClusPhiSize", &clusters_phisize);
0201     events->SetBranchAddress("ClusAdc", &clusters_adc);
0202
0203     std::vector<std::pair<std::shared_ptr<Hit>, std::shared_ptr<Hit>>> tracklets;
0204     std::vector<uint64_t> bcos;
0205
0206     // assemble tracklets that pass dphi cut
0207     for (int i = 0; i < events->GetEntries(); i++)
0208     {
0209         std::vector<std::shared_ptr<Hit>> layer1_clusters;
0210         std::vector<std::shared_ptr<Hit>> layer2_clusters;
0211
0212         events->GetEntry(i);
0213
0214         bcos.push_back(intt_bco);
0215
0216         std::cout << "event " << i << " INTT_BCO " << intt_bco << std::endl;
0217
0218         if (InttBco_IsToBeRemoved == true)
0219             continue;
0220
0221         // use low-multiplicity events
0222         if (Nclusters < 20 || Nclusters > 350)
0223             continue;
0224
0225         for (int j = 0; j < Nclusters; j++)
0226         {
0227             if (clusters_adc->at(j) <= 35)
0228                 continue;
0229
0230             if (clusters_layer->at(j) == 3 || clusters_layer->at(j) == 4)
0231             {
0232                 layer1_clusters.push_back(std::make_shared<Hit>(clusters_x->at(j), clusters_y->at(j), clusters_z->at(j), 0., 0., 0., 0, clusters_phisize->at(j)));
0233             }
0234             if (clusters_layer->at(j) == 5 || clusters_layer->at(j) == 6)
0235             {
0236                 layer2_clusters.push_back(std::make_shared<Hit>(clusters_x->at(j), clusters_y->at(j), clusters_z->at(j), 0., 0., 0., 0, clusters_phisize->at(j)));
0237             }
0238         }
0239
0240         if (layer1_clusters.size() < 10 || layer2_clusters.size() < 10)
0241             continue;
0242
0243         for (auto &cluster1 : layer1_clusters)
0244         {
0245             for (auto &cluster2 : layer2_clusters)
0246             {
0247                 // if (fabs(cluster1->Phi() - cluster2->Phi()) < dphi_cut && fabs(cluster1->posZ() - cluster2->posZ()) < 5.) // bug in calculating delta phi (?);
0248                 if (deltaPhi(cluster1->Phi(), cluster2->Phi()) < dphi_cut)
0249                 {
0250                     tracklets.push_back(std::make_pair(cluster1, cluster2));
0251                 }
0252             }
0253         }
0254     }
0255
0256     Beamspot origin;
0257     // fit_PCAD0Phi0(std::vector<std::pair<std::shared_ptr<Hit>, std::shared_ptr<Hit>>> tracklets, Beamspot BS_init, double d0_cut)
0258     // Beamspot BS = fit_PCAD0Phi0(tracklets, origin, 100.);
0259     Beamspot BS = fit_PCAD0Phi0(tracklets, origin, 1.);
0260     // BS = fit_PCAD0Phi0(tracklets,BS,0.5);
0261     // BS = fit_PCAD0Phi0(tracklets,BS,0.2);
0262     // BS = fit_PCAD0Phi0(tracklets,BS,0.1);
0263
0264     // Get the minimum, median, and maximum values of bcos
0265     std::sort(bcos.begin(), bcos.end());
0266     uint64_t min_bco = bcos.front();
0267     uint64_t median_bco = (bcos.size() % 2 == 0) ? (bcos[bcos.size() / 2] + bcos[bcos.size() / 2 - 1]) / 2 : bcos[bcos.size() / 2];
0268     uint64_t max_bco = bcos.back();
0269
0270     t->Branch("x", &BS.x);
0271     t->Branch("y", &BS.y);
0272     t->Branch("z", &BS.z);
0273     t->Branch("sigma_x", &BS.sigma_x);
0274     t->Branch("sigma_y", &BS.sigma_y);
0275     t->Branch("sigma_z", &BS.sigma_z);
0276     t->Branch("xy_correlation", &BS.xy_correlation);
0277     t->Branch("Ntracklets", &BS.Ntracklets);
0278     t->Branch("intt_bco_min", &min_bco);
0279     t->Branch("intt_bco_median", &median_bco);
0280     t->Branch("intt_bco_max", &max_bco);
0281
0282     t->Fill();
0283     t->Write();
0284
0285     BS.d0phi0dist->Write("d0phi0dist", TObject::kOverwrite);
0286     BS.d0phi0dist_bkgrm->Write("d0phi0dist_bkgrm", TObject::kOverwrite);
0287     BS.d0phi0dist_bkgrm_prof->Write("d0phi0dist_bkgrm_prof", TObject::kOverwrite);
0288     BS.z0dist->Write("z0dist", TObject::kOverwrite);
0289     BS.pcadist->Write("pcadist", TObject::kOverwrite);
0290
0291     BS.identify();
0292
0293     outf->Close();
0294
0295     return 0;
0296 }