macros/exclusive_analysis/jpsi_particle_plot.C

0001 // ----------------------------------------------------------------------------
0002 // TFiles and TTrees
0003 // ----------------------------------------------------------------------------
0004
0005 TFile *f_1 = new TFile("/sphenix/user/gregtom3/data/Summer2018/ExclusiveReco_studies/latest_run/10x100_1M_dvmp.root","READ");
0006
0007 TTree *t_1_exclusive = (TTree*)f_1->Get("event_exclusive");
0008 TTree *t_1_reco = (TTree*)f_1->Get("event_reco");
0009 TTree *t_1_truth = (TTree*)f_1->Get("event_truth");
0010
0011 // ----------------------------------------------------------------------------
0012 // Main Code
0013 // ----------------------------------------------------------------------------
0014
0015 int jpsi_particle_plot()
0016 {
0017   gROOT->SetBatch(kTRUE);
0018   // ----------------------------------------
0019   // Creating Histograms
0020   // ----------------------------------------
0021   const int nbins = 100;
0022   TH2F * h2_eta_phi_em_base = new TH2F("","",nbins,-4,4,nbins,-3.5,3.5);
0023   TH2F * h2_eta_phi_pr_base = new TH2F("","",nbins,-20,20,nbins,-3.5,3.5);
0024   TH2F * h2_mom_eta_em_base = new TH2F("","",nbins,-4,4,nbins,0,20);
0025   TH2F * h2_mom_eta_pr_base = new TH2F("","",nbins,-20,20,nbins,0,150);
0026   TH2F * h2_distance = new TH2F("","",nbins,-20,20,nbins,0,150);
0027   // ----------------------------------------
0028   // Setting Axes Titles
0029   // ----------------------------------------
0030   h2_eta_phi_em_base->GetXaxis()->SetTitle("#eta");
0031   h2_eta_phi_em_base->GetYaxis()->SetTitle("#phi");
0032   h2_eta_phi_pr_base->GetXaxis()->SetTitle("#eta");
0033   h2_eta_phi_pr_base->GetYaxis()->SetTitle("#phi");
0034
0035   h2_mom_eta_em_base->GetXaxis()->SetTitle("#eta");
0036   h2_mom_eta_em_base->GetYaxis()->SetTitle("Momentum [GeV]");
0037   h2_mom_eta_pr_base->GetXaxis()->SetTitle("#eta");
0038   h2_mom_eta_pr_base->GetYaxis()->SetTitle("Momentum [GeV]");
0039   // ----------------------------------------
0040   // Writing out types of hisotgrams
0041   // ----------------------------------------
0042   std::vector<TString> types_eta_phi_em_reco;
0043   std::vector<TString> types_eta_phi_em_true;
0044   std::vector<TString> types_eta_phi_pr;
0045   std::vector<TString> types_mom_eta_em_reco;
0046   std::vector<TString> types_mom_eta_em_true;
0047   std::vector<TString> types_mom_eta_pr;
0048   types_eta_phi_em_reco.push_back("reco_decay_positron_eta_phi");
0049   types_eta_phi_em_reco.push_back("reco_decay_electron_eta_phi");
0050   types_eta_phi_em_reco.push_back("reco_scattered_electron_eta_phi");
0051   types_mom_eta_em_reco.push_back("reco_decay_positron_mom_eta");
0052   types_mom_eta_em_reco.push_back("reco_decay_electron_mom_eta");
0053   types_mom_eta_em_reco.push_back("reco_scattered_electron_mom_eta");
0054   types_eta_phi_em_true.push_back("true_decay_positron_eta_phi");
0055   types_eta_phi_em_true.push_back("true_decay_electron_eta_phi");
0056   types_eta_phi_em_true.push_back("true_scattered_electron_eta_phi");
0057   types_eta_phi_pr.push_back("true_scattered_proton_eta_phi");
0058   types_mom_eta_em_true.push_back("true_decay_positron_mom_eta");
0059   types_mom_eta_em_true.push_back("true_decay_electron_mom_eta");
0060   types_mom_eta_em_true.push_back("true_scattered_electron_mom_eta");
0061   types_mom_eta_pr.push_back("true_scattered_proton_mom_eta");
0062   // ----------------------------------------
0063   // Filling histograms and saving them
0064   // ----------------------------------------
0065   for(unsigned idx = 0 ; idx < types_eta_phi_em_reco.size() ; idx++)
0066     {
0067       fillHist(h2_eta_phi_em_base, t_1_reco, t_1_truth, types_eta_phi_em_reco.at(idx), true);
0068       hist_to_png(h2_eta_phi_em_base,"10x100", types_eta_phi_em_reco.at(idx));
0069     }
0070
0071   for(unsigned idx = 0 ; idx < types_eta_phi_em_true.size() ; idx++)
0072     {
0073       fillHist(h2_eta_phi_em_base, t_1_reco, t_1_truth, types_eta_phi_em_true.at(idx), false);
0074       hist_to_png(h2_eta_phi_em_base,"10x100", types_eta_phi_em_true.at(idx));
0075     }
0076
0077   for(unsigned idx = 0 ; idx < types_eta_phi_pr.size() ; idx++)
0078     {
0079       fillHist(h2_eta_phi_pr_base, t_1_reco, t_1_truth, types_eta_phi_pr.at(idx), false);
0080       hist_to_png(h2_eta_phi_pr_base,"10x100", types_eta_phi_pr.at(idx));
0081     }
0082
0083   for(unsigned idx = 0 ; idx < types_mom_eta_em_reco.size() ; idx++)
0084     {
0085       fillHist(h2_mom_eta_em_base, t_1_reco, t_1_truth, types_mom_eta_em_reco.at(idx), true);
0086       hist_to_png(h2_mom_eta_em_base,"10x100", types_mom_eta_em_reco.at(idx));
0087     }
0088
0089   for(unsigned idx = 0 ; idx < types_mom_eta_em_true.size() ; idx++)
0090     {
0091       fillHist(h2_mom_eta_em_base, t_1_reco, t_1_truth, types_mom_eta_em_true.at(idx), false);
0092       hist_to_png(h2_mom_eta_em_base,"10x100", types_mom_eta_em_true.at(idx));
0093     }
0094
0095
0096   for(unsigned idx = 0 ; idx < types_mom_eta_pr.size() ; idx++)
0097     {
0098       fillHist(h2_mom_eta_pr_base, t_1_reco, t_1_truth, types_mom_eta_pr.at(idx), false);
0099       hist_to_png(h2_mom_eta_pr_base,"10x100", types_mom_eta_pr.at(idx));
0100     }
0101   return 0;
0102 }
0103
0104 void fillHist(TH2F * h, TTree *t_reco, TTree *t_truth, TString type_TString, bool reco)
0105 {
0106   // Create all vectors from tree
0107   std::vector<float> reco_eta;
0108   std::vector<float> true_eta;
0109   std::vector<float> reco_phi;
0110   std::vector<float> true_phi;
0111   std::vector<float> reco_ptotal;
0112   std::vector<float> true_ptotal;
0113   std::vector<int> reco_charge;
0114   std::vector<int> pid;
0115   std::vector<bool> reco_is_scattered_lepton;
0116   std::vector<bool> is_scattered_lepton;
0117   std::vector<float> reco_cluster_e;
0118
0119   // Point to these vectors
0120   std::vector<float>* reco_eta_pointer = &reco_eta;
0121   std::vector<float>* true_eta_pointer = &true_eta;
0122   std::vector<float>* reco_phi_pointer = &reco_phi;
0123   std::vector<float>* true_phi_pointer = &true_phi;
0124   std::vector<float>* reco_ptotal_pointer = &reco_ptotal;
0125   std::vector<float>* true_ptotal_pointer = &true_ptotal;
0126   std::vector<int>* reco_charge_pointer = &reco_charge;
0127   std::vector<int>* pid_pointer = &pid;
0128   std::vector<bool>* reco_is_scattered_lepton_pointer = &reco_is_scattered_lepton;
0129   std::vector<bool>* is_scattered_lepton_pointer = &is_scattered_lepton;
0130   std::vector<float>* reco_cluster_e_pointer = &reco_cluster_e;
0131
0132   // Empty out histogram
0133   h->Reset();
0134
0135   // Set branch addresses
0136   t_reco->SetBranchAddress("eta",&reco_eta_pointer);
0137   t_reco->SetBranchAddress("phi",&reco_phi_pointer);
0138   t_reco->SetBranchAddress("ptotal",&reco_ptotal_pointer);
0139   t_reco->SetBranchAddress("charge",&reco_charge_pointer);
0140   t_reco->SetBranchAddress("is_scattered_lepton",&reco_is_scattered_lepton_pointer);
0141   t_reco->SetBranchAddress("cluster_e",&reco_cluster_e_pointer);
0142
0143   t_truth->SetBranchAddress("geta",&true_eta_pointer);
0144   t_truth->SetBranchAddress("gphi",&true_phi_pointer);
0145   t_truth->SetBranchAddress("gptotal",&true_ptotal_pointer);
0146   t_truth->SetBranchAddress("pid",&pid_pointer);
0147   t_truth->SetBranchAddress("is_scattered_lepton",&is_scattered_lepton_pointer);
0148   // Are we filling in reco branches?
0149   if(reco)
0150     {
0151       // Get number of entries in reco tree
0152       Int_t nentries = Int_t(t_reco->GetEntries());
0153
0154       // Get comparison string for knowing what to fill
0155       char * type = type_TString.Data();
0156       int count = 0;
0157       // Loop over the tree
0158       for(Int_t entryInChain=0; entryInChain<nentries; entryInChain++)
0159     {
0160       Int_t entryInTree = t_reco->LoadTree(entryInChain);
0161       if (entryInTree < 0) break;
0162       t_reco->GetEntry(entryInChain);
0163       t_truth->GetEntry(entryInChain);
0164       if(entryInChain%1000==0)
0165         cout << "Processing event " << count << " of " << nentries << endl;
0166       count++;
0167       // Size of reco branch event
0168       unsigned size_reco = reco_eta.size();
0169
0170       // Loop over a single event
0171       for(unsigned i = 0 ; i < size_reco ; i++)
0172         {
0173           // Booleans for knowing scatter vs decay electrons
0174           bool is_scatter = reco_is_scattered_lepton.at(i);
0175           bool is_decay = !is_scatter;
0176
0177
0178           // Make sure scatter electron eta reco vs. true <0.2
0179           for(unsigned temp = 0 ; temp < true_eta.size() ; temp++)
0180         {
0181           bool is_true_scatter = is_scattered_lepton.at(temp);
0182           float eta_diff = abs(reco_eta.at(i)-true_eta.at(temp));
0183           if(is_true_scatter&&is_scatter&&eta_diff>0.2)
0184             {
0185               is_scatter = false;
0186             }
0187         }
0188
0189           // Ensure charge of electron and positron using truth info
0190
0191           for(unsigned temp = 0 ; temp < true_eta.size() ; temp++)
0192         {
0193           float phi_diff = abs(reco_phi.at(i)-true_phi.at(temp));
0194           float eta_diff = abs(reco_eta.at(i)-true_eta.at(temp));
0195           if(phi_diff<0.5&&eta_diff<0.2&&pid.at(temp)==11)
0196             reco_charge.at(i)==-1;
0197           else if(phi_diff<0.5&&eta_diff<0.2&&pid.at(temp)==-11)
0198             reco_charge.at(i)==1;
0199         }
0200
0201           // If the particle has an eta < -1.434, set the momentum of the particle to cluster energy
0202           if(reco_eta.at(i)<-1.434)
0203         reco_ptotal.at(i) = reco_cluster_e.at(i);
0204
0205
0206           // Use char * type to select what we fill
0207           if(strcmp(type,"reco_decay_positron_eta_phi")==0)
0208         {
0209           if(reco_charge.at(i)==1)
0210             h->Fill(reco_eta.at(i),reco_phi.at(i));
0211         }
0212           else if(strcmp(type,"reco_decay_electron_eta_phi")==0)
0213         {
0214           if(reco_charge.at(i)==-1&&is_decay)
0215             h->Fill(reco_eta.at(i),reco_phi.at(i));
0216
0217         }
0218           else if(strcmp(type,"reco_scattered_electron_eta_phi")==0)
0219         {
0220           if(reco_charge.at(i)==-1&&is_scatter)
0221             h->Fill(reco_eta.at(i),reco_phi.at(i));
0222         }
0223           else if(strcmp(type,"reco_decay_positron_mom_eta")==0)
0224         {
0225           if(reco_charge.at(i)==1)
0226             h->Fill(reco_eta.at(i),reco_ptotal.at(i));
0227         }
0228           else if(strcmp(type,"reco_decay_electron_mom_eta")==0)
0229         {
0230           if(reco_charge.at(i)==-1&&is_decay)
0231             h->Fill(reco_eta.at(i),reco_ptotal.at(i));
0232         }
0233           else if(strcmp(type,"reco_scattered_electron_mom_eta")==0)
0234         {
0235           if(reco_charge.at(i)==-1&&is_scatter)
0236             h->Fill(reco_eta.at(i),reco_ptotal.at(i));
0237         }
0238         }
0239     }
0240     }
0241   else
0242     {
0243       Int_t nentries = Int_t(t_truth->GetEntries());
0244       char * type = type_TString.Data();
0245       for(Int_t entryInChain=0; entryInChain<nentries; entryInChain++)
0246     {
0247       Int_t entryInTree = t_truth->LoadTree(entryInChain);
0248       if (entryInTree < 0) break;
0249       t_truth->GetEntry(entryInChain);
0250
0251       unsigned size_truth = true_eta.size();
0252
0253       for(unsigned i = 0 ; i < size_truth ; i++)
0254         {
0255
0256           bool is_scatter = is_scattered_lepton.at(i);
0257           bool is_decay = !is_scatter;
0258
0259           if(strcmp(type,"true_decay_positron_eta_phi")==0)
0260         {
0261           if(pid.at(i)==-11)
0262             h->Fill(true_eta.at(i),true_phi.at(i));
0263         }
0264           else if(strcmp(type,"true_decay_electron_eta_phi")==0)
0265         {
0266           if(pid.at(i)==11&&is_decay)
0267             h->Fill(true_eta.at(i),true_phi.at(i));
0268         }
0269           else if(strcmp(type,"true_scattered_electron_eta_phi")==0)
0270         {
0271           if(pid.at(i)==11&&is_scatter)
0272             h->Fill(true_eta.at(i),true_phi.at(i));
0273         }
0274           else if(strcmp(type,"true_decay_positron_mom_eta")==0)
0275         {
0276           if(pid.at(i)==-11&&is_decay)
0277             h->Fill(true_eta.at(i),true_ptotal.at(i));
0278         }
0279           else if(strcmp(type,"true_decay_electron_mom_eta")==0)
0280         {
0281           if(pid.at(i)==11&&is_decay)
0282             h->Fill(true_eta.at(i),true_ptotal.at(i));
0283         }
0284           else if(strcmp(type,"true_scattered_electron_mom_eta")==0)
0285         {
0286           if(pid.at(i)==11&&is_scatter)
0287             h->Fill(true_eta.at(i),true_ptotal.at(i));
0288         }
0289           else if(strcmp(type,"true_scattered_proton_mom_eta")==0)
0290         {
0291           if(pid.at(i)==2212)
0292             h->Fill(true_eta.at(i),true_ptotal.at(i));
0293         }
0294           else if(strcmp(type,"true_scattered_proton_eta_phi")==0)
0295         {
0296           if(pid.at(i)==2212)
0297             h->Fill(true_eta.at(i),true_phi.at(i));
0298         }
0299         } // end loop over event
0300     } // end loop over tree
0301     } // end if-else for selecting reco or truth filling
0302 } // end fillHist
0303
0304 void hist_to_png(TH2F * h, TString saveTitle, TString type2)
0305 {
0306   saveTitle = type2 + "_" + saveTitle+".png";
0307   TCanvas *cPNG = new TCanvas(saveTitle,"",700,500);
0308   TImage *img = TImage::Create();
0309   gStyle->SetOptStat(0);
0310   gPad->SetLogz();
0311   char * type = type2.Data();
0312   // Title the histogram //
0313   if(strcmp(type,"reco_decay_electron_eta_phi")==0)
0314     {
0315       h->SetTitle("#eta-#phi space : decay e^{-} reco");
0316     }
0317   else if(strcmp(type,"reco_decay_positron_eta_phi")==0)
0318     {
0319       h->SetTitle("#eta-#phi space : decay e^{+} reco");
0320     }
0321   else if(strcmp(type,"reco_scattered_electron_eta_phi")==0)
0322     {
0323       h->SetTitle("#eta-#phi space : scattered e^{-} reco");
0324     }
0325   else if(strcmp(type,"reco_decay_electron_mom_eta")==0)
0326     {
0327       h->SetTitle("Momentum vs Eta : decay e^{-} reco");
0328     }
0329   else if(strcmp(type,"reco_decay_positron_mom_eta")==0)
0330     {
0331       h->SetTitle("Momentum vs Eta : decay e^{+} reco");
0332     }
0333   else if(strcmp(type,"reco_scattered_electron_mom_eta")==0)
0334     {
0335       h->SetTitle("Momentum vs Eta : scattered e^{-} reco");
0336     }
0337   else if(strcmp(type,"true_decay_electron_eta_phi")==0)
0338     {
0339        h->SetTitle("#eta-#phi space : decay e^{-} truth");
0340     }
0341   else if(strcmp(type,"true_decay_positron_eta_phi")==0)
0342     {
0343       h->SetTitle("#eta-#phi space : decay e^{+} truth");
0344     }
0345   else if(strcmp(type,"true_scattered_electron_eta_phi")==0)
0346     {
0347       h->SetTitle("#eta-#phi space : scattered e^{-} truth");
0348     }
0349   else if(strcmp(type,"true_decay_electron_mom_eta")==0)
0350     {
0351       h->SetTitle("Momentum vs Eta : decay e^{-} truth");
0352     }
0353   else if(strcmp(type,"true_decay_positron_mom_eta")==0)
0354     {
0355       h->SetTitle("Momentum vs Eta : decay e^{+} truth");
0356     }
0357   else if(strcmp(type,"true_scattered_electron_mom_eta")==0)
0358     {
0359       h->SetTitle("Momentum vs Eta : scattered e^{-} truth");
0360     }
0361   else if(strcmp(type,"true_scattered_proton_mom_eta")==0)
0362     {
0363       h->SetTitle("Momentum vs Eta : scattered proton truth");
0364     }
0365   else if(strcmp(type,"true_scattered_proton_eta_phi")==0)
0366     {
0367       h->SetTitle("#eta-#phi space : scattered proton truth");
0368     }
0369
0370   h->Draw("colz");
0371
0372   img->FromPad(cPNG);
0373   img->WriteImage(saveTitle);
0374   delete img;
0375 }