Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:12:55 

0001 // ------------------------------------------------------------------------//
0002 // Code Purpose
0003 // ------------------------------------------------------------------------//
0004 //   To show the difference between the track extrapolation code and truth
0005 //   track information
0006 // ------------------------------------------------------------------------//
0007 //
0008 // ------------------------------------------------------------------------//
0009 //
0010 //
0011 //
0012 //                                 Main Code 
0013 //
0014 //
0015 //
0016 // ------------------------------------------------------------------------//
0017 
0018 const int num_files=4;
0019 
0020 const int verbosity = 1;
0021 
0022 int track2cluster_plot()
0023 {
0024   gROOT->SetBatch(kTRUE);
0025   // ------------------------------------------------------------------------//
0026   //  Uploading Files and Cluster Trees
0027   // ------------------------------------------------------------------------//
0028   
0029   
0030   TFile *f_1 = new TFile("/sphenix/user/gregtom3/data/Summer2018/track2cluster_studies/e1DIS.root");
0031   TTree *t_truth_1 = (TTree*)f_1->Get("event_truth");
0032   TTree *t_cluster_1 = (TTree*)f_1->Get("event_cluster");
0033 
0034   TFile *f_2 = new TFile("/sphenix/user/gregtom3/data/Summer2018/track2cluster_studies/e2DIS.root");
0035   TTree *t_truth_2 = (TTree*)f_2->Get("event_truth");
0036   TTree *t_cluster_2 = (TTree*)f_2->Get("event_cluster");
0037 
0038   TFile *f_5 = new TFile("/sphenix/user/gregtom3/data/Summer2018/track2cluster_studies/e5DIS.root");
0039   TTree *t_truth_5 = (TTree*)f_5->Get("event_truth");
0040   TTree *t_cluster_5 = (TTree*)f_5->Get("event_cluster");
0041 
0042   TFile *f_10 = new TFile("/sphenix/user/gregtom3/data/Summer2018/track2cluster_studies/e10DIS.root");
0043   TTree *t_truth_10 = (TTree*)f_10->Get("event_truth");
0044   TTree *t_cluster_10 = (TTree*)f_10->Get("event_cluster");
0045 
0046   /*TFile *f_20 = new TFile("/sphenix/user/gregtom3/data/Summer2018/track2cluster_studies/e20DIS.root");
0047   TTree *t_truth_20 = (TTree*)f_20->Get("event_truth");
0048   TTree *t_cluster_20 = (TTree*)f_20->Get("event_cluster");*/
0049 
0050   // ------------------------------------------------------------------------//
0051   //  Creating Histograms (Phi, Eta Theta, momentum
0052   //                      posx, posy, posz, spatial distance,
0053   // ------------------------------------------------------------------------//
0054   const int nbins = 200;
0055   TH2F *h2_phi_base = new TH2F("","",nbins,-3.5,3.5,nbins,-3.5,3.5);
0056   TH2F *h2_eta_base = new TH2F("","",nbins,-10,10,nbins,-10,10);
0057   TH2F *h2_theta_base = new TH2F("","",nbins,-3.5,3.5,nbins,-3.5,3.5);
0058   TH2F *h2_mom_base = new TH2F("","",nbins,0,35,nbins,0,35);
0059   TH2F *h2_posx_base = new TH2F("","",nbins,-320,320,nbins,-320,320);
0060   TH2F *h2_posy_base = new TH2F("","",nbins,-320,320,nbins,-320,320);
0061   TH2F *h2_posz_base = new TH2F("","",nbins,-320,320,nbins,-320,320);
0062   TH1F *h1_spatial_base = new TH1F("","",nbins,0,30);
0063   // ------------------------------------------------------------------------//
0064   //  Creating Histogram Clones
0065   // ------------------------------------------------------------------------//
0066   TH2F *h2_phi = (TH2F*)h2_phi_base->Clone();
0067   TH2F *h2_eta = (TH2F*)h2_eta_base->Clone();
0068   TH2F *h2_theta = (TH2F*)h2_theta_base->Clone();
0069   TH2F *h2_mom = (TH2F*)h2_mom_base->Clone();
0070   TH2F *h2_posx = (TH2F*)h2_posx_base->Clone();
0071   TH2F *h2_posy = (TH2F*)h2_posy_base->Clone();
0072   TH2F *h2_posz = (TH2F*)h2_posz_base->Clone();
0073   TH1F *h1_spatial_cemc = (TH1F*)h1_spatial_base->Clone();
0074   TH1F *h1_spatial_eemc = (TH1F*)h1_spatial_base->Clone();
0075   TH1F *h1_spatial_femc = (TH1F*)h1_spatial_base->Clone();
0076   h1_spatial_cemc->SetLineColor(kRed);
0077   h1_spatial_eemc->SetLineColor(kBlue);
0078   h1_spatial_femc->SetLineColor(kGreen);
0079   // ------------------------------------------------------------------------//
0080   //  Which momentums would you like to iterate through?
0081   // ------------------------------------------------------------------------//
0082   std::vector<int> iter_p;
0083   iter_p.push_back(1);
0084   iter_p.push_back(2);
0085   iter_p.push_back(5);   // Only these momenta are currently available
0086   iter_p.push_back(10);
0087   iter_p.push_back(20);    
0088   // ------------------------------------------------------------------------//
0089   //  How many events per tree would you like to plot? nEvents < 0 = all evts
0090   // ------------------------------------------------------------------------//
0091   const int nEvents = 100000;
0092   // ------------------------------------------------------------------------//
0093   //  Analyze each tree
0094   // ------------------------------------------------------------------------//
0095   analyzeTree(h2_phi,h2_eta,h2_theta,h2_mom,h2_posx,h2_posy,h2_posz,h1_spatial_cemc,h1_spatial_eemc,h1_spatial_femc,t_truth_1,t_cluster_1,nEvents);
0096   printHists(h2_phi,h2_eta,h2_theta,h2_mom,h2_posx,h2_posy,h2_posz,h1_spatial_cemc,h1_spatial_eemc,h1_spatial_femc,"1GeV_electron.png");
0097 
0098   analyzeTree(h2_phi,h2_eta,h2_theta,h2_mom,h2_posx,h2_posy,h2_posz,h1_spatial_cemc,h1_spatial_eemc,h1_spatial_femc,t_truth_2,t_cluster_2,nEvents);
0099   printHists(h2_phi,h2_eta,h2_theta,h2_mom,h2_posx,h2_posy,h2_posz,h1_spatial_cemc,h1_spatial_eemc,h1_spatial_femc,"2GeV_electron.png");
0100 
0101   analyzeTree(h2_phi,h2_eta,h2_theta,h2_mom,h2_posx,h2_posy,h2_posz,h1_spatial_cemc,h1_spatial_eemc,h1_spatial_femc,t_truth_5,t_cluster_5,nEvents);
0102   printHists(h2_phi,h2_eta,h2_theta,h2_mom,h2_posx,h2_posy,h2_posz,h1_spatial_cemc,h1_spatial_eemc,h1_spatial_femc,"5GeV_electron.png");
0103 
0104   analyzeTree(h2_phi,h2_eta,h2_theta,h2_mom,h2_posx,h2_posy,h2_posz,h1_spatial_cemc,h1_spatial_eemc,h1_spatial_femc,t_truth_10,t_cluster_10,nEvents);
0105   printHists(h2_phi,h2_eta,h2_theta,h2_mom,h2_posx,h2_posy,h2_posz,h1_spatial_cemc,h1_spatial_eemc,h1_spatial_femc,"10GeV_electron.png");
0106   return 0;
0107 }
0108 void printHists(TH2F *h2_phi, TH2F* h2_eta, TH2F* h2_theta, TH2F* h2_mom, TH2F* h2_posx, TH2F* h2_posy, TH2F* h2_posz, TH1F* h1_spatial_cemc, TH1F *h1_spatial_eemc, TH1F* h1_spatial_femc,TString save)
0109 {
0110   hist_to_png(h2_phi,TString(TString("phi_")+save),"phi");
0111   hist_to_png(h2_eta,TString(TString("eta_")+save),"eta");
0112   hist_to_png(h2_theta,TString(TString("theta_")+save),"theta");
0113   hist_to_png(h2_mom,TString(TString("mom_")+save),"mom");
0114   hist_to_png(h2_posx,TString(TString("posx_")+save),"posx");
0115   hist_to_png(h2_posy,TString(TString("posy_")+save),"posy");
0116   hist_to_png(h2_posz,TString(TString("posz_")+save),"posz");
0117   hist_to_png_h1(h1_spatial_cemc,h1_spatial_eemc,h1_spatial_femc,TString(TString("spatial_")+save),"spatial");
0118 }
0119 void analyzeTree(TH2F *h2_phi, TH2F* h2_eta, TH2F* h2_theta, TH2F* h2_mom, TH2F* h2_posx, TH2F* h2_posy, TH2F* h2_posz, TH1F* h1_spatial_cemc, TH1F* h1_spatial_eemc, TH1F* h1_spatial_femc, TTree* t_truth, TTree* t_cluster, const int nEvents)
0120 {
0121   h2_phi->Reset();
0122   h2_eta->Reset();
0123   h2_theta->Reset();
0124   h2_mom->Reset();
0125   h2_posx->Reset();
0126   h2_posy->Reset();
0127   h2_posz->Reset();
0128   h1_spatial_cemc->Reset();
0129   h1_spatial_eemc->Reset();
0130   h1_spatial_femc->Reset();
0131       // -------------------------------//
0132       // Declare tree variables
0133       // -------------------------------//
0134       //   Track & Cluster
0135       std::vector<float> track_phi; std::vector<float> cluster_phi;
0136       std::vector<float> track_eta; std::vector<float> cluster_eta;
0137       std::vector<float> track_theta; std::vector<float> cluster_theta;
0138       std::vector<float> track_p; std::vector<float> cluster_p; //cluster energy
0139       std::vector<float> track_posx; std::vector<float> cluster_posx;
0140       std::vector<float> track_posy; std::vector<float> cluster_posy;
0141       std::vector<float> track_posz; std::vector<float> cluster_posz;
0142       //   Truth Particle Info
0143       std::vector<float> truth_phi; std::vector<float> truth_eta;
0144       std::vector<float> truth_theta; std::vector<float> truth_p;
0145       // -------------------------------//
0146       // Point to tree variables
0147       // -------------------------------//
0148       //   Track
0149       std::vector<float>* track_phi_pointer = &track_phi;
0150       std::vector<float>* track_eta_pointer = &track_eta;
0151       std::vector<float>* track_theta_pointer = &track_theta;
0152       std::vector<float>* track_p_pointer = &track_p;
0153       std::vector<float>* track_posx_pointer = &track_posx;
0154       std::vector<float>* track_posy_pointer = &track_posy;
0155       std::vector<float>* track_posz_pointer = &track_posz;
0156       //   Cluster
0157       std::vector<float>* cluster_phi_pointer = &cluster_phi;
0158       std::vector<float>* cluster_eta_pointer = &cluster_eta;
0159       std::vector<float>* cluster_theta_pointer = &cluster_theta;
0160       std::vector<float>* cluster_p_pointer = &cluster_p;
0161       std::vector<float>* cluster_posx_pointer = &cluster_posx;
0162       std::vector<float>* cluster_posy_pointer = &cluster_posy;
0163       std::vector<float>* cluster_posz_pointer = &cluster_posz;
0164       //   Truth Particle
0165       std::vector<float>* truth_phi_pointer = &truth_phi;
0166       std::vector<float>* truth_eta_pointer = &truth_eta;
0167       std::vector<float>* truth_theta_pointer = &truth_theta;
0168       std::vector<float>* truth_p_pointer = &truth_p;
0169       // -------------------------------//
0170       // Set tree branch addresses
0171       // -------------------------------//
0172       //   Track
0173       t_cluster->SetBranchAddress("em_track_phi2cluster",&track_phi_pointer);
0174       t_cluster->SetBranchAddress("em_track_eta2cluster",&track_eta_pointer);
0175       t_cluster->SetBranchAddress("em_track_theta2cluster",&track_theta_pointer);
0176       t_cluster->SetBranchAddress("em_track_ptotal",&track_p_pointer);
0177       t_cluster->SetBranchAddress("em_track_x2cluster",&track_posx_pointer);
0178       t_cluster->SetBranchAddress("em_track_y2cluster",&track_posy_pointer);
0179       t_cluster->SetBranchAddress("em_track_z2cluster",&track_posz_pointer);
0180       //   Cluster
0181       t_cluster->SetBranchAddress("em_cluster_phi",&cluster_phi_pointer);
0182       t_cluster->SetBranchAddress("em_cluster_eta",&cluster_eta_pointer);
0183       t_cluster->SetBranchAddress("em_cluster_theta",&cluster_theta_pointer);
0184       t_cluster->SetBranchAddress("em_cluster_e",&cluster_p_pointer);
0185       t_cluster->SetBranchAddress("em_cluster_posx",&cluster_posx_pointer);
0186       t_cluster->SetBranchAddress("em_cluster_posy",&cluster_posy_pointer);
0187       t_cluster->SetBranchAddress("em_cluster_posz",&cluster_posz_pointer);
0188       //   Truth Particle
0189       t_truth->SetBranchAddress("em_evtgen_phi",&truth_phi_pointer);
0190       t_truth->SetBranchAddress("em_evtgen_eta",&truth_eta_pointer);
0191       t_truth->SetBranchAddress("em_evtgen_theta",&truth_theta_pointer);
0192       t_truth->SetBranchAddress("em_evtgen_ptotal",&truth_p_pointer);
0193       // -------------------------------//
0194       // Set tree branch addresses
0195       // -------------------------------//
0196       Int_t nentries_truth =  Int_t(t_truth->GetEntries());
0197       Int_t nentries_cluster= Int_t(t_cluster->GetEntries());
0198       Int_t nentries=0;
0199       // Check if event_truth and event_cluster are equal
0200       
0201       if(nentries_truth!=nentries_cluster)
0202     {
0203       cout << "Warning: Entries in 'event_truth'("<<nentries_truth<<") does not equal Entries in 'event_cluster'("<<nentries_cluster<<"). Using the smaller of the two" << endl;
0204       if(nentries_truth>nentries_cluster) nentries = nentries_cluster;
0205     }
0206       else
0207     {
0208       nentries = nentries_truth;
0209     }
0210 
0211       // Check nEvents with nentries 
0212 
0213       if(nEvents>nentries)
0214     cout << "Warning : Number of Events specified to run is greater than entries in the tree, using all entries" << endl;
0215       else if(nEvents>0)
0216     nentries = nEvents;
0217 
0218       // Iterate through tree
0219       for(Int_t entryInChain=0; entryInChain<nentries; entryInChain++)
0220     {
0221       
0222       // Print out progress
0223       if(entryInChain%1000==0&&verbosity>0)
0224         cout << "Processing event " << entryInChain << " of " << nentries << endl;
0225 
0226       // Not sure what this line does
0227 
0228       Int_t entryInTree_truth = t_truth->LoadTree(entryInChain);
0229       if (entryInTree_truth < 0) break;
0230       Int_t entryInTree_cluster = t_cluster->LoadTree(entryInChain);
0231       if (entryInTree_cluster < 0) break;
0232       
0233       // Get entries (filling vectors)
0234       
0235       t_truth->GetEntry(entryInChain);
0236       t_cluster->GetEntry(entryInChain);
0237 
0238       // Iterate through all truth particles and see what we were able
0239       // to come up with (Fill histograms)
0240       for(unsigned k = 0; k<cluster_p.size(); k++)
0241         {
0242           if(cluster_eta.at(k)<-1.1)
0243         {
0244           h2_phi->Fill(track_phi.at(k),cluster_phi.at(k));
0245           h2_eta->Fill(track_eta.at(k),cluster_eta.at(k));
0246           h2_theta->Fill(track_theta.at(k),cluster_theta.at(k));
0247           h2_mom->Fill(track_p.at(k),cluster_p.at(k));
0248           h2_posx->Fill(track_posx.at(k),cluster_posx.at(k));
0249           h2_posy->Fill(track_posy.at(k),cluster_posy.at(k));
0250           h2_posz->Fill(track_posz.at(k),cluster_posz.at(k));
0251         }
0252           if(cluster_eta.at(k)>1)
0253         {
0254           h1_spatial_femc
0255             ->Fill(sqrt( (track_posx.at(k)-cluster_posx.at(k))*
0256                  (track_posx.at(k)-cluster_posx.at(k))
0257                                   + 
0258                  (track_posy.at(k)-cluster_posy.at(k))*
0259                  (track_posy.at(k)-cluster_posy.at(k))
0260                                   + 0*
0261                  (track_posz.at(k)-cluster_posz.at(k))*
0262                  (track_posz.at(k)-cluster_posz.at(k))));
0263         }
0264           else if(cluster_eta.at(k)<1&&cluster_eta.at(k)>-1.1)
0265         {
0266           h1_spatial_cemc
0267             ->Fill(sqrt( (track_posx.at(k)-cluster_posx.at(k))*
0268                  (track_posx.at(k)-cluster_posx.at(k))
0269                                   + 
0270                  (track_posy.at(k)-cluster_posy.at(k))*
0271                  (track_posy.at(k)-cluster_posy.at(k))
0272                                   + 
0273                  (track_posz.at(k)-cluster_posz.at(k))*
0274                  (track_posz.at(k)-cluster_posz.at(k))));
0275         }
0276           else
0277         {
0278           h1_spatial_eemc
0279             ->Fill(sqrt( (track_posx.at(k)-cluster_posx.at(k))*
0280                  (track_posx.at(k)-cluster_posx.at(k))
0281                                   + 
0282                  (track_posy.at(k)-cluster_posy.at(k))*
0283                  (track_posy.at(k)-cluster_posy.at(k))
0284                                   + 0*
0285                  (track_posz.at(k)-cluster_posz.at(k))*
0286                  (track_posz.at(k)-cluster_posz.at(k))));
0287           
0288         }
0289         }
0290     }
0291 }
0292 void openFiles()
0293 {
0294   
0295 }
0296 
0297 void openTrees()
0298 {
0299   
0300 }
0301 void BinLog(TH2F *h) 
0302 {
0303 
0304    TAxis *axis = h->GetXaxis(); 
0305    int bins = axis->GetNbins();
0306 
0307    Axis_t from = axis->GetXmin();
0308    Axis_t to = axis->GetXmax();
0309    Axis_t width = (to - from) / bins;
0310    Axis_t *new_bins = new Axis_t[bins + 1];
0311 
0312    for (int i = 0; i <= bins; i++) {
0313      new_bins[i] = TMath::Power(10, from + i * width);
0314    } 
0315    axis->Set(bins, new_bins); 
0316 
0317    TAxis *axis2 = h->GetYaxis(); 
0318    int bins2 = axis2->GetNbins();
0319 
0320    Axis_t from2 = axis2->GetXmin();
0321    Axis_t to2 = axis2->GetXmax();
0322    Axis_t width2 = (to2 - from2) / bins2;
0323    Axis_t *new_bins2 = new Axis_t[bins2 + 1];
0324 
0325    for (int i = 0; i <= bins2; i++) {
0326      new_bins2[i] = TMath::Power(10, from2 + i * width2);
0327    } 
0328    axis2->Set(bins2, new_bins2); 
0329 
0330    delete new_bins;
0331    delete new_bins2;
0332 }
0333 
0334 void hist_to_png_h1(TH1F * h_c, TH1F* h_e, TH1F* h_f, TString saveTitle, TString type_string)
0335 {
0336   TCanvas *cPNG = new TCanvas(saveTitle,"",700,500);
0337   TImage *img = TImage::Create();
0338   char * type = type_string.Data();
0339   gStyle->SetOptStat(0);
0340   gPad->SetLogy();
0341   auto legend = new TLegend(0.7,0.7,0.9,0.9);
0342   legend->SetTextSize(0.06);
0343   legend->AddEntry(h_c,"CEMC","l");
0344   legend->AddEntry(h_e,"EEMC","l");
0345   legend->AddEntry(h_f,"FEMC","l");
0346   
0347   if(strcmp(type,"spatial")==0)
0348     {  
0349       h_e->GetXaxis()->SetTitle("Extrapolation distance from cluster [cm]");
0350       h_e->GetYaxis()->SetTitle("Counts");
0351     }
0352   h_e->Draw();
0353   h_c->Draw("SAME");
0354   h_f->Draw("SAME");
0355   legend->Draw();
0356   img->FromPad(cPNG);
0357   img->WriteImage(saveTitle);
0358 }
0359 void hist_to_png(TH2F * h, TString saveTitle, TString type_string)
0360 {
0361   TCanvas *cPNG = new TCanvas(saveTitle,"",700,500);
0362   TImage *img = TImage::Create();
0363   char * type = type_string.Data();
0364   gPad->SetLogz();
0365   gStyle->SetOptStat(0);
0366   if(strcmp(type,"phi")==0)
0367     {
0368       h->GetXaxis()->SetTitle("Track Phi");
0369       h->GetYaxis()->SetTitle("Cluster Phi");
0370     }
0371   else if(strcmp(type,"eta")==0)
0372     {
0373       h->GetXaxis()->SetTitle("Track Eta");
0374       h->GetYaxis()->SetTitle("Cluster Eta"); 
0375     }
0376   else if(strcmp(type,"theta")==0)
0377     {
0378       h->GetXaxis()->SetTitle("Track Theta");
0379       h->GetYaxis()->SetTitle("Cluster Theta");
0380     }
0381   else if(strcmp(type,"mom")==0)
0382     {
0383       h->GetXaxis()->SetTitle("Track Momentum [GeV]");
0384       h->GetYaxis()->SetTitle("Cluster Momentum [GeV]");
0385     }
0386   else if(strcmp(type,"posx")==0)
0387     {
0388       h->GetXaxis()->SetTitle("Track Position X [cm]");
0389       h->GetYaxis()->SetTitle("Cluster Position X [cm]");
0390     }
0391   else if(strcmp(type,"posy")==0)
0392     {
0393       h->GetXaxis()->SetTitle("Track Position Y [cm]");
0394       h->GetYaxis()->SetTitle("Cluster Position Y [cm]");
0395     }
0396   else if(strcmp(type,"posz")==0)
0397     {  
0398       h->GetXaxis()->SetTitle("Track Position Z [cm]");
0399       h->GetYaxis()->SetTitle("Cluster Position Z [cm]");
0400     }
0401   h->Draw("colz9");
0402   img->FromPad(cPNG);
0403   img->WriteImage(saveTitle);
0404 
0405   delete img;
0406 }
[ 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!