macros/exclusive_analysis/dvcs_plot.C

0001 // ------------------------------------------------------------------------//
0002 // Declaring Cut Structs
0003 // ------------------------------------------------------------------------//
0004 struct event_kinematics_cut{
0005   float min_Q2;
0006   float max_Q2;
0007   float min_y;
0008   float max_y;
0009   float min_t;
0010   float max_t;
0011   float min_angle_separation;
0012 };
0013
0014 struct lepton_cut{
0015   float min_eta;
0016   float max_eta;
0017   float min_mom;
0018   float max_mom;
0019   float min_theta;
0020   float max_theta;
0021 };
0022
0023 struct gamma_cut{
0024   float min_eta;
0025   float max_eta;
0026   float min_mom;
0027   float max_mom;
0028   float min_theta;
0029   float max_theta;
0030 };
0031
0032 struct hadron_cut{
0033   float min_eta;
0034   float max_eta;
0035   float min_mom;
0036   float max_mom;
0037   float min_theta;
0038   float max_theta;
0039 };
0040
0041 // ------------------------------------------------------------------------//
0042 // Standard Cuts
0043 // ------------------------------------------------------------------------//
0044 const event_kinematics_cut ev_standard_cut = {0,1000,0,100,0,100,0};
0045 const lepton_cut l_standard_cut = {-999,999,-999,999,-999,999};
0046 const hadron_cut h_standard_cut = {-999,999,-999,999,-999,999};
0047 const gamma_cut g_standard_cut = {-999,999,-999,999,-999,999};
0048 // ------------------------------------------------------------------------//
0049 // Relevant Arrays and Array Index Storing Variable
0050 // ------------------------------------------------------------------------//
0051 TH2F **h2array;
0052 event_kinematics_cut *ev_cutarray;
0053 lepton_cut *l_cutarray;
0054 gamma_cut *g_cutarray;
0055 hadron_cut *h_cutarray;
0056 TString *pngNames;
0057 int index=0;
0058 int arraySize;
0059 // ------------------------------------------------------------------------//
0060 // Other Global Variables
0061 // ------------------------------------------------------------------------//
0062 int nEvents;
0063 // ------------------------------------------------------------------------//
0064 // Creating Base Histograms for Types of Plots
0065 // ------------------------------------------------------------------------//
0066
0067 const int nbins=80;
0068 TH2F *h2_x_Q2_base = new TH2F("x_Q2","",nbins,-4,0,nbins,0,2);
0069 TH2F *h2_lepton_base = new TH2F("lepton_mom_eta","",nbins,-4,4,nbins,0,30);
0070 TH2F *h2_gamma_base = new TH2F("gamma_mom_eta","",nbins,-4,4,nbins,0,30);
0071 TH2F *h2_hadron_base = new TH2F("hadron_mom_eta","",nbins,-1,20,nbins,0,300);
0072 TH2F *h2_hadron_scatter_base = new TH2F("scatter","",nbins,0,1.5,nbins,0,0.006);
0073 TH2F *h2_electron_photon_angle_base = new TH2F("angle_separation","",nbins,0,3.2,nbins,0,1);
0074 TH2F *h2_nils_plot = new TH2F("nils_plot","",nbins,0,9,nbins,0,360);
0075
0076 // ------------------------------------------------------------------------//
0077 //
0078 //
0079 //
0080 //                                 Main Code
0081 //
0082 //
0083 //
0084 // ------------------------------------------------------------------------//
0085
0086 int dvcs_plot()
0087 {
0088   // ------------------------------------------------------------------------//
0089   //  Uploading File and Truth Tree
0090   // ------------------------------------------------------------------------//
0091
0092   TFile *f = new TFile("/sphenix/user/gregtom3/data/Summer2018/DVCS_studies/500k_10x250_sartre_dvcs.root","READ");
0093   TTree *t= (TTree*)f->Get("event_truth");
0094
0095   // ------------------------------------------------------------------------//
0096   // Creating Plotting Histograms
0097   // 1) Edit 'arraySize' to fit how many plots you will like to make
0098   // 2) Create your own event kinematics, lepton, hadron, or gamma cuts
0099   //    Or simply select the standard ones already created
0100   // 3) Add the desired histogram to the array using 'add_to_array'
0101   // ------------------------------------------------------------------------//
0102
0103   arraySize = 1;
0104
0105   h2array=malloc(arraySize * sizeof(TH2F));
0106   ev_cutarray=malloc(arraySize * sizeof(event_kinematics_cut));
0107   l_cutarray=malloc(arraySize * sizeof(lepton_cut));
0108   g_cutarray=malloc(arraySize * sizeof(gamma_cut));
0109   h_cutarray=malloc(arraySize * sizeof(hadron_cut));
0110   pngNames=malloc(arraySize * sizeof(TString));
0111
0112   // event_kinematics_cut *** = {float min_Q2, float max_Q2,
0113   //             float min_y , float max_y ,
0114   //             float min_t , float max_t ,
0115   //             float min_angle_separation}
0116
0117   // lepton_cut *** = {float min_eta, float max_eta,
0118   //                   float min_mom, float max_mom,
0119   //                   float min_theta, float max_theta}
0120
0121   // gamma_cut *** = {float min_eta, float max_eta,
0122   //                   float min_mom, float max_mom,
0123   //                   float min_theta, float max_theta}
0124
0125   // hadron_cut *** = {float min_eta, float max_eta,
0126   //                   float min_mom, float max_mom,
0127   //                   float min_theta, float max_theta}
0128
0129   gamma_cut EEMC = {-4,-1,-999,999,-999,999};
0130   gamma_cut CEMC = {-1,1,-999,999,-999,999};
0131   gamma_cut FEMC = {1,4,-999,999,-999,999};
0132
0133   // Adding to Histogram Plot Instructions
0134   // add_to_array(char *arg_1,  TString arg_2, TString arg_3
0135   //              event_kinematics_cut arg_4, lepton_cut arg_5, gamma_cut arg_6,                  hadron_cut arg_7)
0136   //
0137   // arg_1 : Type of Histogram to draw
0138   //       a) x_Q2
0139   //       b) gamma_mom_eta
0140   //       c) lepton_mom_eta
0141   //       d) proton_mom_eta
0142   //       e) scatter
0143   //       f) angle_separation
0144   //       g) nils_plot
0145   // arg_2 : After code is executed, histogram saves as png
0146   // arg_3 : Histogram Title
0147   // arg_4 : Cut on event kinematics and other special properties
0148   // arg_5 : Cut on the scattered electron's kinematics
0149   // arg_6 : Cut on the scattered photon's kinematics
0150   // arg_7 : Cut on the scattered proton's kinematics
0151   // **Note: Leaving any of arg_4 -> arg_7 empty defaults to standard cuts
0152   //         ev_standard_cut , l_standard_cut , g_standard_cut , h_standard_cut
0153
0154   add_to_array("x_Q2","/sphenix/user/gregtom3/data/Summer2018/DVCS_studies/png/milou_cuts/10x250_x_Q2_EEMC_plot.png","x-Q2 Distribution 10x250 DVCS EEMC (LoI cuts)",ev_standard_cut,l_standard_cut,EEMC,h_standard_cut);
0155
0156
0157
0158   // ------------------------------------------------------------------------//
0159   // Events to Run (If nEvents < 0, all events used)
0160   // ------------------------------------------------------------------------//
0161
0162   nEvents = -1;
0163
0164   // ------------------------------------------------------------------------//
0165   // Declare all relevant event variables
0166   // ------------------------------------------------------------------------//
0167
0168   Float_t Q2;
0169   Float_t W;
0170   Float_t x;
0171   Float_t y;
0172   Float_t T;
0173
0174   std::vector<float> energy;
0175   std::vector<float> eta;
0176   std::vector<float> pid;
0177   std::vector<float> theta;
0178   std::vector<float> phi;
0179   std::vector<float> xvector;
0180
0181   std::vector<float>* energy_pointer=&energy;
0182   std::vector<float>* eta_pointer=&eta;
0183   std::vector<float>* pid_pointer=&pid;
0184   std::vector<float>* theta_pointer=&theta;
0185   std::vector<float>* phi_pointer=&phi;
0186   std::vector<float>* xvector_pointer=&xvector;
0187
0188   t->SetBranchAddress("evtgen_Q2",&Q2);
0189   t->SetBranchAddress("evtgen_x",&x);
0190   t->SetBranchAddress("evtgen_y",&y);
0191   t->SetBranchAddress("evtgen_W",&W);
0192   t->SetBranchAddress("evtgen_t",&T);
0193   t->SetBranchAddress("em_evtgen_ptotal",&energy_pointer);
0194   t->SetBranchAddress("em_evtgen_eta",&eta_pointer);
0195   t->SetBranchAddress("em_evtgen_pid",&pid_pointer);
0196   t->SetBranchAddress("em_evtgen_theta",&theta_pointer);
0197   t->SetBranchAddress("em_evtgen_phi",&phi_pointer);
0198   t->SetBranchAddress("em_reco_x_e",&xvector_pointer);
0199
0200   // Momentum, Eta, Theta , Phi of Scattered Photon
0201   float gamma_eta;
0202   float gamma_mom;
0203   float gamma_theta;
0204   float gamma_phi;
0205
0206   // Momentum, Theta, Phi of Scattered Electron
0207   float lepton_eta;
0208   float lepton_mom;
0209   float lepton_theta;
0210   float lepton_phi;
0211
0212   // Momentum, Theta, Phi of Scattered Proton
0213   float hadron_eta;
0214   float hadron_mom;
0215   float hadron_theta;
0216   float hadron_phi;
0217
0218   // ------------------------------------------------------------------------//
0219   //  Loop through Tree
0220   // ------------------------------------------------------------------------//
0221
0222   Int_t nentries = Int_t(t->GetEntries());
0223   if(nEvents>nentries)
0224     cout << "Warning : Number of Events specified to run is greater than entries in the tree, using all entries" << endl;
0225   else if(nEvents>0)
0226     nentries = nEvents;
0227   for(Int_t entryInChain=0; entryInChain<nentries; entryInChain++)
0228     {
0229       if(entryInChain%1000==0)
0230     cout << "Processing event " << entryInChain << " of " << nentries << endl;
0231       Int_t entryInTree = t->LoadTree(entryInChain);
0232       if (entryInTree < 0) break;
0233       t->GetEntry(entryInChain);
0234       // Flip through all the particles in the event (only 3)
0235       for(unsigned i = 0; i<energy.size(); i++)
0236     {
0237       if(pid.at(i)==22) // Scattered photon
0238         {
0239           gamma_eta = eta.at(i);
0240           gamma_mom = energy.at(i);
0241           gamma_theta = theta.at(i);
0242           gamma_phi = phi.at(i);
0243         }
0244       if(pid.at(i)==11) // Scattered electron
0245         {
0246           lepton_eta = eta.at(i);
0247           lepton_mom = energy.at(i);
0248           lepton_theta = theta.at(i);
0249           lepton_phi = phi.at(i);
0250           x=xvector.at(i);
0251         }
0252       if(pid.at(i)==2212) // Scattered proton
0253         {
0254           hadron_eta = eta.at(i);
0255           hadron_mom = energy.at(i);
0256           hadron_theta = theta.at(i);
0257           hadron_phi = phi.at(i);
0258         }
0259     }
0260
0261       // Angle between Scattered Electron and Photon
0262       float angle_separation = get_angle_separation(lepton_phi, lepton_eta, gamma_phi, gamma_eta);
0263
0264       // Set 't' to abs for plotting
0265       if(T<0)
0266     T=-T;
0267       // Fill the histogram array
0268       for(int n = 0; n<arraySize; n++)
0269     {
0270       // This is really ugly, but we check each event parameters with the
0271       // cuts predefined by the plot.
0272       event_kinematics_cut myCut = ev_cutarray[n];
0273       lepton_cut lCut = l_cutarray[n];
0274       gamma_cut gCut = g_cutarray[n];
0275       hadron_cut hCut = h_cutarray[n];
0276       char *histTitle = h2array[n]->GetName();
0277       float sep = angle_separation;
0278
0279       bool passed_Q2 = (Q2 >= myCut.min_Q2 && Q2<= myCut.max_Q2);
0280       bool passed_y = (y >= myCut.min_y && y <= myCut.max_y);
0281       bool passed_t = (T >= myCut.min_t && T <= myCut.max_t);
0282       bool passed_separation = ( sep >= myCut.min_angle_separation );
0283       bool ev_passed = passed_Q2 && passed_y && passed_t && passed_separation;
0284
0285       bool passed_lepton_eta = (lepton_eta >= lCut.min_eta && lepton_eta <= lCut.max_eta);
0286       bool passed_lepton_mom = (lepton_mom >= lCut.min_mom && lepton_mom <= lCut.max_mom);
0287       bool passed_lepton_theta = (lepton_theta >= lCut.min_theta && lepton_theta <= lCut.max_theta);
0288       bool lepton_passed = passed_lepton_eta && passed_lepton_mom && passed_lepton_theta;
0289
0290       bool passed_gamma_eta = (gamma_eta >= gCut.min_eta && gamma_eta <= gCut.max_eta);
0291       bool passed_gamma_mom = (gamma_mom >= gCut.min_mom && gamma_mom <= gCut.max_mom);
0292       bool passed_gamma_theta = (gamma_theta >= gCut.min_theta && gamma_theta <= gCut.max_theta);
0293       bool gamma_passed = passed_gamma_eta && passed_gamma_mom && passed_gamma_theta;
0294
0295       bool passed_hadron_eta = (hadron_eta >= hCut.min_eta && hadron_eta <= hCut.max_eta);
0296       bool passed_hadron_mom = (hadron_mom >= hCut.min_mom && hadron_mom <= hCut.max_mom);
0297       bool passed_hadron_theta = (hadron_theta >= hCut.min_theta && hadron_theta <= hCut.max_theta);
0298       bool hadron_passed = passed_hadron_eta && passed_hadron_mom && passed_hadron_theta;
0299
0300       bool passed = ev_passed && lepton_passed && gamma_passed && hadron_passed;
0301       if(passed)
0302         {
0303           if(strcmp(histTitle,"x_Q2")==0)
0304         {
0305           h2array[n]->Fill(x,Q2);
0306         }
0307           else if(strcmp(histTitle,"lepton_mom_eta")==0)
0308         {
0309           h2array[n]->Fill(lepton_eta,lepton_mom);
0310         }
0311           else if(strcmp(histTitle,"gamma_mom_eta")==0)
0312         {
0313           h2array[n]->Fill(gamma_eta,gamma_mom);
0314         }
0315           else if(strcmp(histTitle,"hadron_mom_eta")==0)
0316         {
0317           h2array[n]->Fill(hadron_eta,hadron_mom);
0318         }
0319           else if(strcmp(histTitle,"scatter")==0)
0320         {
0321           h2array[n]->Fill(T,hadron_theta);
0322         }
0323           else if(strcmp(histTitle,"angle_separation")==0)
0324         {
0325           h2array[n]->Fill(lepton_theta,angle_separation);
0326         }
0327           else if(strcmp(histTitle,"nils_plot")==0)
0328         {
0329           float eta_diff = lepton_eta-gamma_eta;
0330           float phi_diff = lepton_phi-gamma_phi;
0331           if(eta_diff<0) eta_diff=-eta_diff;
0332           if(phi_diff<0) phi_diff=-phi_diff;
0333           h2array[n]->Fill(eta_diff,phi_diff*180/3.14159265);
0334         }
0335           else
0336         {
0337           cout << "Error: Histogram type '"<< histTitle << "' not found. Please use one of the specified types" << endl;
0338           return -1;
0339         }
0340         }
0341     }
0342     }
0343
0344   // ------------------------------------------------------------------------//
0345   //  Save PNG
0346   // ------------------------------------------------------------------------//
0347   for(int n = 0; n<arraySize; n++)
0348     {
0349       hist_to_png(h2array[n],pngNames[n]);
0350     }
0351 return 0;
0352
0353 }
0354
0355 float get_angle_separation(float e_phi, float e_theta, float g_phi, float g_theta)
0356 {
0357    float photon_x = cos(g_phi)*sin(g_theta);
0358    float photon_y = sin(g_phi)*sin(g_theta);
0359    float photon_z = cos(g_theta);
0360
0361    float electron_x = cos(e_phi)*sin(e_theta);
0362    float electron_y = sin(e_phi)*sin(e_theta);
0363    float electron_z = cos(e_theta);
0364
0365    float dot_prod = photon_x*electron_x+photon_y*electron_y+photon_z*electron_z;
0366
0367    if(dot_prod>=1)
0368      return 0;
0369    else if(dot_prod<=-1)
0370      return 3.14159265;
0371    else
0372      return acos(dot_prod);
0373 }
0374
0375 void hist_to_png(TH2F * h, TString saveTitle)
0376 {
0377   TCanvas *cPNG = new TCanvas(saveTitle,"",700,500);
0378   TImage *img = TImage::Create();
0379   gStyle->SetOptStat(0);
0380   if(strcmp(h->GetName(),"x_Q2")==0)
0381     {
0382       gPad->SetLogx();
0383       gPad->SetLogy();
0384       h->GetXaxis()->SetTitle("x");
0385       h->GetYaxis()->SetTitle("Q2 [GeV^{2}]");
0386     }
0387   else if(strcmp(h->GetName(),"gamma_mom_eta")==0)
0388     {
0389       h->GetXaxis()->SetTitle("Photon #eta");
0390       h->GetYaxis()->SetTitle("Photon Momentum [GeV]");
0391       gPad->SetLogz();
0392     }
0393   else if(strcmp(h->GetName(),"lepton_mom_eta")==0)
0394     {
0395       h->GetXaxis()->SetTitle("Electron #eta");
0396       h->GetYaxis()->SetTitle("Electron Momentum [GeV]");
0397       gPad->SetLogz();
0398     }
0399   else if(strcmp(h->GetName(),"hadron_mom_eta")==0)
0400     {
0401       h->GetXaxis()->SetTitle("Proton #eta");
0402       h->GetYaxis()->SetTitle("Proton Momentum [GeV]");
0403       gPad->SetLogz();
0404     }
0405   else if(strcmp(h->GetName(),"scatter")==0)
0406     {
0407       h->GetXaxis()->SetTitle("|t| [GeV^2]");
0408       h->GetYaxis()->SetTitle("Proton Scatter Angle [rad]");
0409     }
0410   else if(strcmp(h->GetName(),"angle_separation")==0)
0411     {
0412       h->GetXaxis()->SetTitle("Electron #theta");
0413       h->GetYaxis()->SetTitle("Electron Photon Angle Separation [rad]");
0414     }
0415   else if(strcmp(h->GetName(),"nils_plot")==0)
0416     {
0417       h->GetXaxis()->SetTitle("Electron Photon Eta Separation");
0418       h->GetYaxis()->SetTitle("Electron Photon Phi Separation [deg]");
0419     }
0420   h->Draw("colz9");
0421   img->FromPad(cPNG);
0422   img->WriteImage(saveTitle);
0423
0424   delete img;
0425 }
0426
0427 void add_to_array(char *type, TString pngName, TString title,
0428           event_kinematics_cut ev_c=ev_standard_cut,
0429           lepton_cut           l_c=l_standard_cut,
0430           gamma_cut            g_c=g_standard_cut,
0431           hadron_cut           h_c=h_standard_cut)
0432 {
0433   if(index+1>arraySize)
0434     {
0435       cout << "Warning: Set 'arraySize' to a value which matches the number of histograms you wish to plot" << endl;
0436       cout << "Not plotting Plot Index " << index << endl;
0437       exit();
0438     }
0439   else
0440     {
0441   TH2F *theH=NULL;
0442   if(strcmp(type,"x_Q2")==0)
0443     {
0444       theH = (TH2F*)h2_x_Q2_base->Clone();
0445       BinLog(theH);
0446     }
0447   else if(strcmp(type,"lepton_mom_eta")==0)
0448     {
0449       theH = (TH2F*)h2_lepton_base->Clone();
0450     }
0451   else if(strcmp(type,"gamma_mom_eta")==0)
0452     {
0453       theH = (TH2F*)h2_gamma_base->Clone();
0454     }
0455   else if(strcmp(type,"hadron_mom_eta")==0)
0456     {
0457       theH = (TH2F*)h2_hadron_base->Clone();
0458     }
0459   else if(strcmp(type,"scatter")==0)
0460     {
0461       theH = (TH2F*)h2_hadron_scatter_base->Clone();
0462     }
0463   else if(strcmp(type,"angle_separation")==0)
0464     {
0465       theH = (TH2F*)h2_electron_photon_angle_base->Clone();
0466     }
0467   else if(strcmp(type,"nils_plot")==0)
0468     {
0469       theH = (TH2F*)h2_nils_plot->Clone();
0470     }
0471   h2array[index] = theH;
0472   ev_cutarray[index] = ev_c;
0473   l_cutarray[index] = l_c;
0474   g_cutarray[index] = g_c;
0475   h_cutarray[index] = h_c;
0476   pngNames[index] = pngName;
0477   index++;
0478   theH->SetTitle(title);
0479     }
0480 }
0481
0482 void BinLog(TH2F *h)
0483 {
0484
0485    TAxis *axis = h->GetXaxis();
0486    int bins = axis->GetNbins();
0487
0488    Axis_t from = axis->GetXmin();
0489    Axis_t to = axis->GetXmax();
0490    Axis_t width = (to - from) / bins;
0491    Axis_t *new_bins = new Axis_t[bins + 1];
0492
0493    for (int i = 0; i <= bins; i++) {
0494      new_bins[i] = TMath::Power(10, from + i * width);
0495    }
0496    axis->Set(bins, new_bins);
0497
0498    TAxis *axis2 = h->GetYaxis();
0499    int bins2 = axis2->GetNbins();
0500
0501    Axis_t from2 = axis2->GetXmin();
0502    Axis_t to2 = axis2->GetXmax();
0503    Axis_t width2 = (to2 - from2) / bins2;
0504    Axis_t *new_bins2 = new Axis_t[bins2 + 1];
0505
0506    for (int i = 0; i <= bins2; i++) {
0507      new_bins2[i] = TMath::Power(10, from2 + i * width2);
0508    }
0509    axis2->Set(bins2, new_bins2);
0510
0511    delete new_bins;
0512    delete new_bins2;
0513 }