analysis/runTrackingDiagnostics/plot_kshort.C

0001 #include <cmath>
0002 #include <TFile.h>
0003 #include <TNtuple.h>
0004 #include <TH2D.h>
0005 #include <TCut.h>
0006 #include <Eigen/Dense>
0007
0008
0009 void plot_kshort(std::string inputFile = "./kshort_G4sPHENIX.root", const std::string outputFile = "./kshort_plots.root")
0010 {
0011
0012
0013   bool plotBool = true; // Determines whether plots are displayed or not
0014   bool savePlot = true; // Determines whether histograms are saved to output file or not
0015   bool localVerbosity = true;
0016   bool analyzeMassSpectrum = true;
0017
0018   TFile fin(inputFile.c_str());
0019
0020   TNtuple *ntuple;
0021   fin.GetObject("ntp_reco_info",ntuple);
0022
0023
0024   /* Below are values in massRecoAnalysis ntuple 3/6/23
0025 "ntp_reco_info","decay_pairs","x1:y1:z1:px1:py1:pz1:dca3dxy1:dca3dz1:phi1:pca_rel1_x:pca_rel1_y:pca_rel1_z:eta1:charge1:tpcClusters_1:x2:y2:z2:px2:py2:pz2:dca3dxy2:dca3dz2:phi2:pca_rel2_x:pca_rel2_y:pca_rel2_z:eta2:charge2:tpcClusters_2:vertex_x:vertex_y:vertex_z:pair_dca:invariant_mass:invariant_pt:pathlength_x:pathlength_y:pathlength_z:pathlength:rapidity:pseudorapidity:projected_pos1_x:projected_pos1_y:projected_pos1_z:projected_pos2_x:projected_pos2_y:projected_pos2_z:projected_mom1_x:projected_mom1_y:projected_mom1_z:projected_mom2_x:projected_mom2_y:projected_mom2_z:projected_pca_rel1_x:projected_pca_rel1_y:projected_pca_rel1_z:projected_pca_rel2_x:projected_pca_rel2_y:projected_pca_rel2_z:projected_pair_dca:projected_pathlength_x:projected_pathlength_y:projected_pathlength_z:projected_pathlength:quality1:quality2:cosThetaReco"
0026   */
0027
0028
0029   float invariant_mass;
0030   float projected_pathlength_x;
0031   float projected_pathlength_y;
0032   float projected_pathlength_z;
0033   float projected_pathlength;
0034   float vertex_x;
0035   float vertex_y;
0036   float vertex_z;
0037   float projected_mom1_x;
0038   float projected_mom1_y;
0039   float projected_mom1_z;
0040   float projected_mom2_x;
0041   float projected_mom2_y;
0042   float projected_mom2_z;
0043   float projected_pair_dca;
0044   float invariant_pt;
0045   float quality1;
0046   float quality2;
0047   float dca3dxy1;
0048   float dca3dxy2;
0049   float dca3dz1;
0050   float dca3dz2;
0051
0052
0053   ntuple->SetBranchAddress("invariant_mass",&invariant_mass);
0054   ntuple->SetBranchAddress("projected_pathlength_x",&projected_pathlength_x);
0055   ntuple->SetBranchAddress("projected_pathlength_y",&projected_pathlength_y);
0056   ntuple->SetBranchAddress("projected_pathlength_z",&projected_pathlength_z);
0057   ntuple->SetBranchAddress("projected_pathlength",&projected_pathlength);
0058   ntuple->SetBranchAddress("vertex_x",&vertex_x);
0059   ntuple->SetBranchAddress("vertex_y",&vertex_y);
0060   ntuple->SetBranchAddress("vertex_z",&vertex_z);
0061   ntuple->SetBranchAddress("projected_mom1_x",&projected_mom1_x);
0062   ntuple->SetBranchAddress("projected_mom1_y",&projected_mom1_y);
0063   ntuple->SetBranchAddress("projected_mom1_z",&projected_mom1_z);
0064   ntuple->SetBranchAddress("projected_mom2_x",&projected_mom2_x);
0065   ntuple->SetBranchAddress("projected_mom2_y",&projected_mom2_y);
0066   ntuple->SetBranchAddress("projected_mom2_z",&projected_mom2_z);
0067   ntuple->SetBranchAddress("projected_pair_dca",&projected_pair_dca);
0068   ntuple->SetBranchAddress("invariant_pt",&invariant_pt);
0069   ntuple->SetBranchAddress("quality1",&quality1);
0070   ntuple->SetBranchAddress("quality2",&quality2);
0071   ntuple->SetBranchAddress("dca3dxy1",&dca3dxy1);
0072   ntuple->SetBranchAddress("dca3dxy2",&dca3dxy2);
0073   ntuple->SetBranchAddress("dca3dz1",&dca3dz1);
0074   ntuple->SetBranchAddress("dca3dz2",&dca3dz2);
0075
0076
0077   int entries = ntuple->GetEntries();
0078
0079   TH1D *invariant_mass_hist = new TH1D("invariant_mass_hist","Invariant Mass",1000,0.35,0.65);
0080   invariant_mass_hist->SetMinimum(0);
0081   TH2D *pathlengthpathlengthz = new TH2D("pathlengthpathlengthz","path v. pathlength_z",5000,-20.0,20.0,5000,0.0,10.0);
0082   TH2D *pathMass = new TH2D("pathMass","Invariant Mass vs. Path Length",5000,0.0,1.0,5000,0.0,10.0);
0083
0084
0085   int qual_cut = 10;
0086   float pathlength_cut = 0.09;
0087   double track_dca_cut = 0.0105;
0088
0089   for(int i=0; i<entries; ++i)
0090     {
0091       ntuple->GetEntry(i);
0092
0093       Eigen::Vector3f mom1(projected_mom1_x,projected_mom1_y,projected_mom1_z);
0094       Eigen::Vector3f mom2(projected_mom2_x,projected_mom2_y,projected_mom2_z);
0095       Eigen::Vector3f projected_momentum = mom1 + mom2;
0096       Eigen::Vector3f pathLength (projected_pathlength_x,projected_pathlength_y,projected_pathlength_z);
0097       float pathLengthMag = pathLength.norm();
0098
0099       Eigen::Vector3f vertex(vertex_x, vertex_y, vertex_z);
0100       float costheta = pathLength.dot(projected_momentum)/(projected_momentum.norm()*pathLength.norm());
0101       //float radius = sqrt(pow(projected_pathlength_x,2) + pow(projected_pathlength_y,2))
0102       /*
0103       //Aligned optimized cuts
0104       if(costheta < 0.9995) continue;
0105       //if(abs(projected_pathlength_x) < 0.1 || abs(projected_pathlength_y) < 0.1) continue;
0106       if(abs(projected_pathlength_x) < pathlength_cut || abs(projected_pathlength_y) < pathlength_cut) continue;
0107       if(abs(projected_pair_dca) > 0.035) continue;
0108       if(invariant_pt < 0.1) continue;
0109       if(quality1 > qual_cut || quality2 > qual_cut) continue;
0110       if(dca3dxy1 < track_dca_cut || dca3dxy2 < track_dca_cut || dca3dz1 < track_dca_cut || dca3dz2 < track_dca_cut) continue;
0111       */
0112       if(costheta < 0.99) continue;
0113       if(abs(projected_pathlength_x) < pathlength_cut || abs(projected_pathlength_y) < pathlength_cut) continue;
0114       if(abs(projected_pair_dca) > 0.055) continue;
0115       if(invariant_pt < 0.1) continue;
0116       if(quality1 > qual_cut || quality2 > qual_cut) continue;
0117       if(dca3dxy1 < track_dca_cut || dca3dxy2 < track_dca_cut || dca3dz1 < track_dca_cut || dca3dz2 < track_dca_cut) continue;
0118
0119       invariant_mass_hist->Fill(invariant_mass);
0120       pathlengthpathlengthz->Fill(projected_pathlength_z,pathLengthMag);
0121       pathMass->Fill(invariant_mass,pathLengthMag);
0122     }
0123
0124   TF1* f = new TF1("f","gaus(0)+[3]+[4]*x",0.45,0.55);
0125
0126   if(plotBool)
0127     {
0128       TCanvas *c1 = new TCanvas("c1","",10,10,800,800);
0129
0130       f->SetParameter(0,20);    //100 for aligned data
0131       f->SetParameter(1,0.497);
0132       f->SetParameter(2,0.020); //0.010 for aligned data
0133       f->SetParameter(3,0.0);
0134       f->SetParameter(4,0.0);
0135
0136       invariant_mass_hist->Fit("f","R L");//R intially for perfect geometry
0137       invariant_mass_hist->DrawCopy();
0138     }
0139
0140
0141   int binLow      = invariant_mass_hist->FindBin(0.48);
0142   int binHigh     = invariant_mass_hist->FindBin(0.52);
0143   double integral = invariant_mass_hist->Integral(binLow,binHigh);
0144
0145   TF1* fbackground  = new TF1("fbackground","[0]+[1]*x");
0146   fbackground->SetParameter(0,f->GetParameter(3));
0147   fbackground->SetParameter(1,f->GetParameter(4));
0148
0149   double backgroundIntegral    = fbackground->Integral(0.485,0.51);
0150   double foregroundIntegral    = f->Integral(0.485,0.51);
0151   double signal                = (foregroundIntegral - backgroundIntegral);
0152   double signalBackgroundRatio = signal/backgroundIntegral;
0153
0154
0155   TF1* fsignal = new TF1("fsignal","gaus(0)");
0156   fsignal->SetParameter(0,f->GetParameter(0));
0157   fsignal->SetParameter(1,f->GetParameter(1));
0158   fsignal->SetParameter(2,f->GetParameter(2));
0159   double signal_integral = fsignal->Integral(0.35,0.65);
0160   double binwidth = invariant_mass_hist->GetBinWidth(1);
0161
0162   if(plotBool)
0163     {
0164       TCanvas *c2 = new TCanvas("c2","",10,10,600,600); // path v path_z
0165       pathlengthpathlengthz->DrawCopy();
0166
0167       TCanvas *c3 = new TCanvas("c3","",10,10,600,600);
0168       pathMass->DrawCopy();
0169     }
0170
0171   if(localVerbosity)
0172     {
0173       std::cout << " Integral: "<< integral << std::endl;
0174       std::cout<< "signal: " << signal<<"  Signal to background ratio: " << signalBackgroundRatio << std::endl;
0175       std::cout << " Signal Integral: "<< signal_integral<< " yield: "<< signal_integral/binwidth<<std::endl;
0176     }
0177
0178   if(savePlot)
0179     {
0180       std::unique_ptr<TFile> kShortFile(TFile::Open(outputFile.c_str(), "recreate"));
0181       kShortFile->WriteObject(invariant_mass_hist,"invariant_mass");
0182       kShortFile->WriteObject(pathlengthpathlengthz,"pathlength_v_pathlengthz");
0183       kShortFile->WriteObject(pathMass,"pathlength_v_invariantMass");
0184     }
0185
0186 }