StRoot: macros/qa/Plot_Embed.C Source File

00001 
00002 
00004 // $Id: Plot_Embed.C,v 1.1 2007/04/27 15:03:04 cristina Exp $
00005 // owner: Cristina
00006 //
00007 // What it does: This Macro reads two diferent scanned files and compare them. It takes the 3D histogram (nch,pt, nfit)
00008 // and makes the projection on the z-axis (nfit) and give a 1D  plot for Nfit for both files
00009 // with extension .gif - When running this macro, include 2 bins for Pt:pt1 and pt2, check the location for the imput  files and where
00010 // the output will be saved
00011 //Plot created are: Dca, Nfit, EnergyLoss, Efficiecy, Delta_Vertex (X, Y, Z), reference multilicity
00012 //              
00014 
00015 
00016 
00017 #ifndef __CINT__
00018 #include "TROOT.h"
00019 #include "TSystem.h"
00020 #include <iostream.h>
00021 #include "TH1.h"
00022 #include "TH2.h"
00023 #include "TH3.h"
00024 #include "TFile.h"
00025 #include "TTree.h"
00026 #include "TChain.h"
00027 #include "TTreeHelper.h"
00028 #endif
00029 
00030 void plot_Nfit(float pt1 = 0.2, float pt2 = 2.0, int nch1 = 0, int nch2 = 1000, Int_t id=6){
00031 
00032 gROOT->LoadMacro("~/mymacros/Utility.C");
00033 gROOT->LoadMacro("~/mymacros/UtilityHist.C");
00034 gStyle->SetOptStat(0);
00035 gStyle->SetOptDate(0);
00036 //gStyle->SetOptTitle(0);       
00037 gStyle->SetOptFit(0);   
00038 
00039 char* oFile1 = "~/outputs/Phi/out_scan_embed_mc_Phi1.root";
00040 char* oFile2 = "~/outputs/out_scan_embed_mc_KMinus.root";
00041 
00042 //Opening file 1
00043  f1 = new TFile(oFile1);
00044 
00045    TH3D* hDca = hMcDca->Clone(); hDca->SetName("hDca");
00046    TH3D* hNfit = hNfit3->Clone();hNfit->SetName("hNfit");
00047    TH2F* dEdx_ep = (TH2F *) ((TH2F *) f1->Get("dEdx"))->Clone();
00048 
00049  //Opening file 2
00050 
00051 f2 = new TFile(oFile2);
00052 
00053    TH3D* hDcaR = hMcDcaR->Clone();hDcaR->SetName("hDcaR");
00054    TH3D* hNfit = hNfit3->Clone();hNfit->SetName("hNfit");
00055    TH2F* dEdx_pi = (TH2F *) ((TH2F *) f2->Get("dEdx"))->Clone();
00056 
00057    // --------------------------------------Plotting DCA-----------------------
00058 
00059         TCanvas *c= new TCanvas("c","DCA",400, 400);
00060         
00061         c->SetGridx(0);
00062         c->SetGridy(0);
00063         c->SetLeftMargin(0.15);
00064         c->SetRightMargin(0.05);
00065         c->cd;
00066 
00067 //Bins for Multiplicity - for oFile1
00068         TH1D *hX =  (TH1D*)hDca->Project3D("X");
00069         Int_t bin_nch1 = hX->FindBin(nch1);
00070         Int_t bin_nch2 = hX->FindBin(nch2);
00071 
00072 //Bins for Pt - for oFile1
00073         TH1D *hY =  (TH1D*)hDca->Project3D("Y");
00074         Int_t bin_pt1 = hY->FindBin(pt1);
00075         Int_t bin_pt2 = hY->FindBin(pt2);
00076 
00077 //ProjectionZ (DCA) for oFile1 is done in nch bins and pt bins found above
00078 
00079         TH1D *hDca1 =(TH1D*) hDca->ProjectionZ("hDca1",bin_nch1, bin_nch2, bin_pt1, bin_pt2);
00080         
00081         hDca1->SetLineColor(2);// Data in Red - eplus
00082         hDca1->SetXTitle("Dca (cm)");
00083         hDca1->SetYTitle("Counts");
00084         hDca1->GetYaxis()->SetTitleOffset(1.8);
00085         hDca1->GetYaxis()->SetDecimals();
00086        
00087         
00088         float mscale = hDca->GetMaximum();
00089         
00090         
00091 //Bins for Multiplicity - for oFile2
00092         TH1D *hX_r =  (TH1D*)hDcaR->Project3D("X");
00093         Int_t bin_nch1_r = hX_r->FindBin(nch1);
00094         Int_t bin_nch2_r = hX_r->FindBin(nch2);
00095 
00096 //Bins for Pt -  for oFile2
00097         TH1D *hY_r =  (TH1D*)hDcaR->Project3D("Y");
00098         Int_t bin_pt1_r = hY_r->FindBin(pt1);
00099         Int_t bin_pt2_r = hY_r->FindBin(pt2);
00100 
00101 //ProjectionZ (DCA)  for oFile2
00102 
00103         TH1D *hDca1_r=(TH1D*)hDcaR->ProjectionZ("hDca1_r",bin_nch1_r, bin_nch2_r, bin_pt1_r, bin_pt2_r);
00104 
00105         Double_t Sum_Real = hDca1_r->GetSum();
00106         Double_t Sum_MC = hDca1->GetSum();
00107         hDca1_r->Scale(1./Sum_Real); //Scaling  Data
00108         hDca1->Scale(1./Sum_MC); 
00109         
00110         hDca1->Draw();
00111         hDca1_r->Draw("same")
00112 
00113         //-----------------------------------Plotting Nfit-------------------------------------
00114 
00115 
00116         TCanvas *c1= new TCanvas("c1","NFIT",400, 400);
00117         c1->SetGridx(0);
00118         c1->SetGridy(0);
00119         c1->SetLeftMargin(0.15);
00120         c1->SetRightMargin(0.05);
00121         c1->cd;
00122 
00123 
00124  //Bins for Multiplicity - oFile 1
00125         TH1D *hX =  (TH1D*)hNfit->Project3D("X");
00126         Int_t multL = hX->FindBin(nch1);
00127         Int_t multH = hX->FindBin(nch2);
00128 
00129 
00130 // Bins for Pt - oFile 1
00131         TH1D *hY =  (TH1D*)hNfit->Project3D("Y");
00132         Int_t bin_pt1 = hY->FindBin(pt1);
00133         Int_t bin_pt2 = hY->FindBin(pt2);
00134 
00135 
00136 // Projection Z(Nfit) for oFile 1
00137         TH1D *hNfitz =(TH1D*) hNfit->ProjectionZ("hNfitz",multL, multH, bin_pt1, bin_pt2);
00138 
00139         hNfitz->SetLineColor(2);
00140         hNfitz->SetXTitle("Nfit");
00141         hNfitz->SetYTitle("Counts");
00142         hNfitz->GetYaxis()->SetDecimals();
00143         hNfitz->GetYaxis()->SetTitleOffset(1.7);
00144         
00145         Double_t Sum_MC= hNfitz->GetSum();
00146 
00147 
00148 //Bins for Multiplicity - oFile 2
00149         
00150         TH1D *hX_r = (TH1D*) hNfitR->Project3D("X");
00151         Int_t multLr = hX_r->FindBin(nch1);
00152         Int_t multHr = hX_r->FindBin(nch2);
00153 
00154 
00155 //Bins for Pt -  oFile 2
00156         TH1D *hY_r =  (TH1D*)hNfitR->Project3D("Y");
00157         Int_t bin_pt1_r = hY_r->FindBin(pt1);
00158         Int_t bin_pt2_r = hY_r->FindBin(pt2);
00159 
00160 // Projection Z(Nfit) for oFile 2
00161 
00162         TH1D *hNfitz_r = (TH1D*)hNfitR->ProjectionZ("hNfitz_r",multLr, multHr, bin_pt1_r, bin_pt2_r);
00163 
00164         Double_t sum_Real = hNfitz_r->GetSum();
00165         hNfitz_r->Scale(1./sum_Real); //Scaling real Data
00166         hNfitz->Scale(1./Sum_MC); //Scaling MC Data
00167 
00168         //------------------- Plotting Global Variables-----------------------------
00169         // Delta Vertex position (dvx, dvY, dvZ and vx_vy)
00170         
00171         TCanvas *c2= new TCanvas("c2","Delta Vertex Position",400, 400);
00172         c2->SetGridx(0);
00173         c2->SetGridy(0);
00174         c2->cd;
00175         
00176         v_xy->Draw();
00177 
00178 
00179         TCanvas *c3= new TCanvas("c3","Delta Vertex Position",400, 400);
00180         c3->SetGridx(0);
00181         c3->SetGridy(0);
00182         c3->SetLeftMargin(0.15);
00183         c3->SetRightMargin(0.05);
00184         c3->cd;
00185         c3->Divide(3,1);
00186 
00187         c3_1->cd(); //3 plots in the same canvas
00188         dvx->Draw();
00189 
00190         c3_2->cd()
00191         dvy->Draw();
00192         
00193         c3_3->cd()
00194         dvz->Draw();
00195 
00196         //--------------Plotting Reference Multiplicity---------------------------
00197 
00198         TCanvas *c4= new TCanvas("c4","Reference Multiplicity",400, 400);
00199         c4->SetGridx(0);
00200         c4->SetGridy(0);
00201         c4->cd;
00202         
00203         hMult ->Fill(*nPos + *nNeg); // coming from MuDst ??
00204         hMult->Draw();
00205 
00206         //----------------------Energy Loss------------------
00207         //For Primary Tracks
00208         
00209         TCanvas *c5= new TCanvas("c5","Energy Loss",400, 400);
00210         c5->SetGridx(0);
00211         c5->SetGridy(0);
00212         c5->cd;
00213 
00214         hPtM_E_Pr ->Fill(PtPrMatched[itr],PtPrMatched[itr]-PtMcMatched[itr]);
00215         hPtM_E_Pr ->SetLineColor(9);
00216         hPtM_E_Pr ->Draw();
00217 
00218         //For Global Tracks
00219         TCanvas *c6= new TCanvas("c6","Energy Loss",400, 400);
00220         c6->SetGridx(0);
00221         c6->SetGridy(0);
00222         c6->cd;
00223 
00224         hPtM_E_Gl -> Fill (PtPrGlobal[itr],PtPrGlobal[itr]-PtMcGlobal[itr]);
00225         hPtM_E_Gl -> SetLineColor(8);
00226         hPtM_E_Gl -> Draw();
00227 
00228         //-------------------------Efficiency--------------------------
00229 
00230         TCanvas *c7= new TCanvas("c7","Energy Loss",400, 400);
00231         c7->SetGridx(0);
00232         c7->SetGridy(0);
00233         c7->cd;
00234 
00235         //for Primary / Embedded_MC
00236 
00237 
00238         //for Global / Embedded_MC
00239 
00240         //-------------------------MIPS----------------------
00241         
00242         TCanvas *c8= new TCanvas("c8","Energy Loss",400, 400);
00243         c8->SetGridx(0);
00244         c8->SetGridy(0);
00245         c8->Divide(2,1);
00246         c8_1->cd();     
00247         
00248         dEdx_ep->SetMarkerColor(2);//E-plus
00249         //dEdx_ep->Draw();
00250         dEdx_pi->ProfileX();
00251         //dEdx_pi->Draw("same");//Pi 
00252         dEdx_pi->SetXTitle("Momentum P (GeV/c)");
00253         dEdx_pi->SetYTitle("dE/dx");
00254         c8->cd;
00255         
00256                 
00257         //Bins for Pt - Pion
00258         TH1D *dEdx_pi_x = (TH1D*)dEdx_pi->ProjectionX();
00259         Int_t bin_pt1_r = dEdx_pi_x->FindBin(pt1);
00260         Int_t bin_pt2_r = dEdx_pi_x->FindBin(pt2);
00261 
00262         //ProjectionY (dEdx) - Pion
00263  
00264         TH1D *dEdx_pi_y = (TH1D*)dEdx_pi->ProjectionY("py",bin_pt1_r, bin_pt2_r);
00265         dEdx_pi_y->SetAxisRange(0.,4.,"X");
00266         dEdx_pi_y->SetXTitle("dE / dx");
00267         
00268         // Gaussian Fit for Pion Data
00269 
00270         Double_t par[6];
00271         g1    = new TF1("g1","gaus",0.8.,1.6);
00272         g2    = new TF1("g2","gaus",1.7,1.9.);
00273         PiFitPiE2 = new TF1("PiFitPiE2","([0]/[2]*exp(-0.5*pow((x-[1])/[2],2))+[3]/[5]*exp(-0.5*pow((x-[4])/[5],2)))/sqrt(6.283)");//Fitting pi Data with 2 Gaussians
00274 
00275                 PiFitPiE2->SetRange(0.0,4.0);
00276                 dEdx_pi_y->Fit(g1,"R");
00277                 dEdx_pi_y->Fit(g2,"R");
00278    
00279         g1->GetParameters(&par[0]);
00280         g2->GetParameters(&par[3]);
00281  
00282         PiFitPiE2->SetParameters(par);
00283         PiFitPiE2->SetLineWidth(2);
00284         PiFitPiE2->SetLineColor(3);
00285         dEdx_pi_y->Fit(PiFitPiE2,"R");
00286         keyLine(.5,.7,"dEdx Pi",1);
00287 
00288 
00289         //----------------PROJECTION Y  DATA-----------------------------------------
00290         
00291         c8_2->cd();
00292  
00293         //Bins for Pt -EPlus
00294 
00295         TH1D *dEdx_ep_x = (TH1D*)dEdx_ep->ProjectionX();
00296         Int_t bin_pt1 = dEdx_ep_x->FindBin(pt1);
00297         Int_t bin_pt2 = dEdx_ep_x->FindBin(pt2); 
00298 
00299 //ProjectionY (dEdx) - EPlus
00300 
00301         TH1D *dEdx_ep_y = (TH1D*)dEdx_ep->ProjectionY("py1",bin_pt1, bin_pt2);
00302         dEdx_ep_y->SetAxisRange(0.,4.,"X");
00303         dEdx_ep_y->Draw();
00304 
00305 // Gaussian Fit for ME_plus Data        
00306 
00307         g6  = new TF1("g6","gaus",.6,4.);
00308         g6->SetLineWidth(2);
00309         g6->SetLineColor(2);
00310         dEdx_ep_y->Fit(g6,"R");
00311         Double_t Mean = g6->GetParameters(&par[2]);
00312         
00313         cout<< "Mean is : "<<Mean<<endl;
00314         keyLine(.5,.8,"dEdx EPlus",2);
00315         keyLine(.5, .6, "Mean Value : Mean",2) ;
00316         
00317         cout<<"bin_pt1_r"<<bin_pt1_r<<endl;
00318         cout<<"bin_pt2_r"<<bin_pt2_r<<endl;
00319         cout<<"bin_pt1"<<bin_pt1<<endl;
00320         cout<<"bin_pt2"<<bin_pt2<<endl;
00321         
00322         
00323         //------------------------------------------------------------------------
00324 char title[100], gif[100];
00325 
00326  sprintf(gif,"~/Plots/dedx_%.2fpT%.2f_%d.jpg",pt1,pt2);
00327 
00328 dEdx_ep_y->SetTitle(title);
00329 dEdx_pi_y->SetTitle(title);
00330 dEdx_ep_y->Draw();
00331 
00332 
00333         hDca1->SetTitle(title);
00334         hDca->Write();
00335         hNfit->Write();
00336         hMult->Write();
00337 return;
00338 }