StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
plot_embed.C
1 //First run scan_embed.C to generate root file with all the histograms
2 // V. May 31 2007 - Cristina
3 //==================================================
4 //V2 Revised February 2009 - Cristina
5 // plots of dca and nfit (eta and pt) dependant simultaneously are being added
6 //============================================================
7 
8 
9 #ifndef __CINT__
10 #include "TROOT.h"
11 #include "TSystem.h"
12 #include <iostream.h>
13 #include "TH1.h"
14 #include "TH2.h"
15 #include "TH3.h"
16 #include "TFile.h"
17 #include "TTree.h"
18 #include "TChain.h"
19 #include "TTreeHelper.h"
20 #include "TText.h"
21 #include "TLatex.h"
22 #include "TAttLine.h"
23 #include "TCanvas.h"
24 #endif
25 
26 static Int_t iCanvas = 1 ;
27 //P08id_testXXXXXX_Phi_PiPlus_histo_040709.rootP08id_test_Phi_PiPlus_histo_031509.root"
28 
29 void plot_embed(
30  const Char_t* inputFileName = "~/Embed/P08id_testXXXXXX_Phi_PiPlus_histo_040709.root",
31  const Int_t id=11
32  )
33 {
34  Int_t plot = 0;
35  const TString PathSave = "/Embed/";
36 
37  //id -> is for the reconstructed track or Daughter
38 
39  gROOT->LoadMacro("StRoot/macros/embedding/Utility.C");
40  // gStyle->SetOptStat(0);
41  gStyle->SetOptFit(0);
42 
43  //Labels ->These will appear in plots and when saving them
44 
45  Char_t *Prod ="P08id";
46  Char_t *Particle = "Phi";//This is the particle embedded
47  Char_t *Daughter ="Pion";//reconstructing on this Daughter
48  Char_t *Date= "031509";
49 
50  char text1[80];
51  sprintf(text1,"P08id Phi");//this is going to show in all the histograms
52  TLatex t;
53  t.SetTextSize(0.07);
54  t.SetTextColor(1);
55  t.SetNDC();
56 
57  char title[100];
58  char gif[400];
59 
60  float mass2 = -1.0;
61  TString tag("K");
62  TString tagTex("K");
63 
64  if (id == 8) { tag = "Piplus"; tagTex = "#pi^{+}"; Char_t * Daughter = "Piplus" ;mass2 = 0.019;}
65  if (id == 9) { tag = "Piminus"; tagTex = "#pi^{-}"; Char_t * Daughter = "Piminus"; mass2 = 0.019;}
66  if (id == 11) { tag = "Kplus"; tagTex = "K^{+}"; Char_t * Daughter = "Kplus"; mass2 = 0.245;}
67  if (id == 12) { tag = "Kminus"; tagTex = "K^{-}"; Char_t * Daughter = "Kminus";mass2 = 0.245;}
68  if (id == 14) { tag = "Proton"; tagTex = "p"; Char_t * Daughter = "Proton";mass2 = 0.880;}
69  if (id == 15) { tag = "Pbar"; tagTex = "#bar{p}"; Char_t * Daughter = "Pbar"; mass2 = 0.880;}
70  if (id == 50) { tag = "Phi"; tagTex = "#phi"; mass2 = 1.020;}
71  if (id == 2) { tag = "Eplus"; tagTex = "e^{+}"; mass2 = 0.511;}
72  if (id == 1) { tag = "Dmeson"; tagTex = "D"; mass2 = 1.864;}
73 
74 
75  //==============================================================
76  //Cloning Histograms
77  //==============================================================
78 
79 
80  TFile* file = TFile::Open(inputFileName);
81  if(!file || !file->IsOpen()){
82  cout << "can't open " << inputFileName << endl;
83  return;
84  }
85 
86  TH3D *hDca1 = (TH3D*) file-> Get("hDca");//DCA
87  TH3D *hNfit1 = (TH3D*) file-> Get("hNfit");//Nfit
88 
89 
90  TH2D *hPtM_E1 = (TH2D*) file-> Get("hPtM_E");//Energy Loss
91  TH2F *dedx1 = (TH2F*) file-> Get("dedx");
92 
93  TH3D *hDcaG1 = (TH3D*) file-> Get("hDcaG");//DCA
94  TH3D *hNfitG1 = (TH3D*) file-> Get("hNfitG");//Nfit
95 
96  TH2F *dedxG1 = (TH2F*) file-> Get("dedxG");
97  TH1D *PhiG1 = (TH1D*) file-> Get("PhiG");
98 
99  TH2D *vxy1 = (TH2D*) file-> Get("vxy");
100 
101  TH1D *vz1 = (TH1D*) file-> Get("vz");
102  TH1D *dvx1 = (TH1D*) file-> Get("dvx");
103  TH1D *dvy1 = (TH1D*) file-> Get("dvy");
104  TH1D *dvz1 = (TH1D*) file-> Get("dvz");
105 
106  TH1D *PhiMc1 = (TH1D*) file-> Get("PhiMc");
107  TH1D *EtaMc1 = (TH1D*) file-> Get("EtaMc");
108  TH1D *YMc1 = (TH1D*) file-> Get("YMc");
109  TH1D *PtMc1 = (TH1D*) file-> Get("PtMc");
110  TH3D *PtYPhiMc1= (TH3D*) file-> Get("PtYPhiMc");
111 
112  TH1D *PhiM1 = (TH1D*) file-> Get("PhiM");
113  TH1D *EtaM1 = (TH1D*) file-> Get("EtaM");
114  TH1D *YM1 = (TH1D*) file-> Get("YM");
115  TH1D *PtM1 = (TH1D*) file-> Get("PtM");
116  TH3D *PtYPhiM1 = (TH3D*) file-> Get("PtYPhiM");
117 
118  TH2D *PtM_eff1 = (TH2D*) file->Get("PtM_E");//efficiency
119 
120  //if you have MuDst hist
121 
122  TH2D *dedx1R = (TH2D*) file-> Get("dedxR");
123 
124  //3D histos
125 
126  TH3D *hDcaM = (TH3D*) file-> Get("hPtEtaDcaM");//DCA
127  TH3D *hDcaG = (TH3D*) file-> Get("hPtEtaDcaG");//DCA
128  TH3D *hDcaR = (TH3D*) file-> Get("hPtEtaDcaR");
129 
130  TH3D *hNfitM = (TH3D*) file-> Get("hPtEtaNfitM");//Nfit
131  TH3D *hNfitR = (TH3D*) file-> Get("hPtEtaNfitR");//
132  TH3D *hNfitG = (TH3D*) file-> Get("hPtEtaNfitG");//Nfit
133 
134 
135  const int nch1 = 0;
136  const int nch2 = 1000;
137 
138  //The Following bins are for the 3X3 matrix of eta,pt d
139 
140  const int nEtaBin=3;
141  const Double_t eta[nEtaBin+1]= {0.2,0.5,0.8,1.0};
142 
143  const int nPtBin=3;
144  const Double_t pt[nPtBin+1] = {0.5,0.6,0.8,1.0};
145 
146 
147  //=============================================================
148  //Delta Vertex position
149  //=========================================================
150  /*
151  TCanvas *cdv= new TCanvas("cdv","Delta Vertex position",700, 500);
152  cdv->Divide(3, 1);
153 
154  cdv_1->cd();
155  // dvx1->Rebin(2);
156  dvx1->SetXTitle("#Delta VX");
157  dvx1->SetAxisRange(-1.5, 1.5, "X");
158  dvx1->Draw();
159 
160  cdv_2->cd();
161  dvy1->SetXTitle("#Delta VY");
162  dvy1->SetAxisRange(-1.5, 1.5, "X");
163  dvy1->Draw();
164 
165  cdv_3->cd();
166  dvz1->SetXTitle("#Delta VZ");
167  dvz1->SetAxisRange(-1.5, 1.5, "X");
168  dvz1->Draw();
169 
170 
171  t.DrawLatex(0.3, 0.9, text1);
172 
173  cdv->Update();
174 
175 if (plot)
176 {
177  sprintf(gif,"~/Embed/%s/%s/DeltaVertex_%s_%s_%s_%sXXXXXXXXX.gif",Prod, Particle, Prod, Particle,Daughter, Date);
178  // sprintf(gif,"~/Embed/Test_DeltaVertex_%s_%s_%s_%sXXXXXXXXX.gif", Prod, Particle,Daughter, Date);
179  // sprintf(gif,PathSave+"Test.gif");//, Prod, Particle,Daughter, Date);
180  cdv->SaveAs(gif);
181 }
182 
183 // return;
184 
185 
186  //=============================================================
187  //Vertex position
188  //=========================================================
189 
190 
191  TCanvas *cv= new TCanvas("cv","Vertex position",600, 400);
192  cv->Divide(2,1);
193 
194  cv_1->cd();
195  vz1->Rebin(2);
196  vz1->SetXTitle("Vertex Z");
197  vz1->Draw();
198 
199  cv_2->cd();
200 
201  vxy1->Draw("colz");
202  vxy1->SetAxisRange(-1.5, 1.5, "X");
203  vxy1->SetAxisRange(-1.5, 1.5, "Y");
204  vxy1->SetXTitle ("vertex X");
205  vxy1->SetYTitle ("vertex Y");
206 
207  //keyLine(.2, 1.05,Prod,1);
208  //keyLine(.2, 1.00,Particle,1);
209  t.DrawLatex(0.3, 0.9, text1);
210 
211  cv->Update();
212 
213  if (plot)
214  {
215  sprintf(gif,"~/Embed/%s/%s/Vertex_%s_%s_%s_%s.gif",Prod, Particle, Prod, Particle, Daughter,Date);
216  cv->SaveAs(gif);
217  }
218  // return;
219 
220  //=============================================================
221  //dEdx
222  //=========================================================
223 
224 
225  TCanvas *cdedx= new TCanvas("cdedx","dEdx vs P",600, 600);
226  // gStyle->SetPalette(1);
227 
228  cdedx->SetGridx(0);
229  cdedx->SetGridy(0);
230  cdedx->SetLeftMargin(0.15);
231  cdedx->SetRightMargin(0.05);
232  cdedx->cd();
233 
234  dedxG1->SetTitle("dEdx vs P ");
235  dedxG1->SetXTitle("Momentum (GeV/c)");
236  dedxG1->SetYTitle("dE/dx (x10^{-6})");
237 
238  dedxG1->SetMarkerStyle(1);
239  dedxG1->SetMarkerSize(0.2);
240  dedxG1->SetMarkerColor(4);
241  dedxG1->SetAxisRange(0, 8., "X");
242  dedxG1->SetAxisRange(0, 8., "Y");
243 
244 
245  dedx1->SetTitle("dEdx vs P ");
246  dedx1->SetXTitle("Momentum (GeV/c)");
247  dedx1->SetYTitle("dE/dx (x10^{-6})");
248  dedx1->SetAxisRange(0, 8., "X");
249  dedx1->SetAxisRange(0, 8., "Y");
250  dedx1->SetMarkerStyle(1);
251  dedx1->SetMarkerSize(2);
252  dedx1->SetMarkerColor(2);
253 
254  dedx1R->SetAxisRange(0, 8., "X");
255  dedx1R->SetAxisRange(0, 8., "Y");
256  dedx1R->SetMarkerStyle(1);
257  dedx1R->SetMarkerSize(0.2);
258  dedx1R->SetMarkerColor(1);
259 
260  //dedx1R->Draw();
261  dedxG1->Draw();//"colz");
262  dedx1->Draw("same");
263  // dedx1R->Draw("same");
264  t.SetTextSize(0.05);
265  t.DrawLatex(0.3, 0.9, text1);
266  //keyLine(.2, 1.05,Prod,1);
267  //keyLine(.2, 1.00,Particle,1);
268 
269  keySymbol(.4, .87, "Matched Tracks_"+tagTex,2,20, 0.03);
270  keySymbol(.4, .82, "Ghost Tracks",4,20,0.03);
271  // keySymbol(.4, .82, "MuDst Tracks_"+tagTex,1,20,0.03);
272 
273  cdedx->Update();
274  if(plot)
275  {
276 
277  // sprintf(gif,"~/Embed/%s/%s/Dedx_%s_%s_%s_%s.gif",Prod, Particle, Prod, Particle,Daughter, Date);
278  // cdedx->SaveAs(gif);
279  }
280  // return;
281 
282 
283  //=============================================================
284  //pt vs y of mc tracks
285  //=========================================================
286 
287 
288  TCanvas *c11= new TCanvas("c11","PT VS Y",400, 400);
289 
290  c11->SetGridx(0);
291  c11->SetGridy(0);
292  c11->cd();
293 
294  TH2D *hPt_Y = (TH2D*)PtYPhiMc1->Project3D("XY");
295  hPt_Y->SetTitle("pT Vs Y of MC tracks");
296  hPt_Y->SetYTitle ("pT (GeV/c)");
297  hPt_Y->SetXTitle ("Y");
298  hPt_Y->SetAxisRange(-1.2, 1.2, "X");
299  hPt_Y->SetAxisRange(0, 2.5, "Y");
300  hPt_Y->Draw();
301 
302  t.SetTextSize(0.05);
303  t.DrawLatex(0.3, 0.9, text1);
304 
305  c11->Update();
306  // sprintf(gif,"~/embedding/%s/%s/pT_Y_%s_%s_%s_%s.gif",Prod, Particle, Prod, Particle,Daughter, Date);
307  // c11->SaveAs(gif);
308 
309 
310  //=============================================================
311  //pt vs phi of mc tracks
312  //=========================================================
313 
314 
315 
316  TCanvas *c12= new TCanvas("c12","pT VS PHI",400, 400);
317 
318  c12->SetGridx(0);
319  c12->SetGridy(0);
320  c12->cd();
321 
322  TH2D *hPt_Phi = (TH2D*)PtYPhiMc1->Project3D("XZ");
323  hPt_Phi->SetTitle("pT Vs Phi of MC tracks");
324  hPt_Phi->SetYTitle ("p_{T} (GeV/c)");
325  hPt_Phi->SetXTitle ("Phi");
326  hPt_Phi->SetAxisRange(-3.4, 3.4, "X");
327  hPt_Phi->SetAxisRange(0, 2.5, "Y");
328  hPt_Phi->Draw();
329  t.SetTextSize(0.05);
330  t.DrawLatex(0.3, 0.9, text1);
331  // keySymbol(.2, 1.05,Prod,1);
332  // keySymbol(.4, 1.05,Particle,1);
333 
334  c12->Update();
335  // sprintf(gif,"~/embedding/%s/%s/pT_Phi_%s_%s_%s_%s.gif",Prod, Particle, Prod, Particle, Daughter, Date);
336  // c12->SaveAs(gif);
337  // return;
338 
339  //==================================================
340  //y vs phi of mc tracks
341  //===================================================
342  TCanvas *cyphi= new TCanvas("cyphi","Y VS PHI",400, 400);
343 
344  cyphi->SetGridx(0);
345  cyphi->SetGridy(0);
346  cyphi->cd();
347 
348  TH2D *hY_Phi = (TH2D*)PtYPhiMc1->Project3D("YZ");
349  hY_Phi->SetTitle("Y Vs Phi of MC tracks");
350  hY_Phi->SetYTitle ("Y (GeV/c)");
351  hY_Phi->SetXTitle ("Phi");
352  hY_Phi->SetAxisRange(-3.4, 3.4, "X");
353  hY_Phi->SetAxisRange(-1.2, 1.2, "Y");
354  hY_Phi->Draw();
355 
356  // keySymbol(.2, 1.05,Prod,1);
357  // keySymbol(.4, 1.05,Particle,1);
358  t.SetTextSize(0.05);
359  t.DrawLatex(0.3, 0.9, text1);
360  cyphi->Update();
361  if (plot)
362  {
363 
364  sprintf(gif,"~/embedding/%s/%s/Y_phi_%s_%s_%s_%s.gif",Prod, Particle, Prod, Particle, Daughter,Date);
365  cyphi->SaveAs(gif);
366  }
367 
368  // return;
369 
370  //=====================================================
371  //pt
372  //=========================================================
373 
374  TCanvas *cpt= new TCanvas("cpt","Pt",500, 500);
375 
376  cpt->SetGridx(0);
377  cpt->SetGridy(0);
378  cpt->cd();
379  gStyle->SetOptDate(0);
380 
381  PtM1->Rebin(2);
382  PtM1->Scale(1/2);
383  PtM1->Scale(PtM1->GetSum());
384 
385  PtMc1->Rebin(2);
386  PtMc1->Scale(1/2);
387  PtMc1->Scale(PtMc1->GetSum());
388 
389  PtMc1->SetLineColor(1);
390  PtMc1->SetTitle("p_{T}");
391  PtMc1->SetXTitle ("pT (GeV/c)");
392 
393  PtMc1->SetAxisRange(0.0, 8, "X");
394  PtMc1->SetAxisRange(0.0, 1.2*PtM1->GetMaximum(), "Y");
395  PtMc1->Draw();
396 
397  PtM1->SetLineColor(2);
398  PtM1->Draw("same");
399 
400  t.SetTextSize(0.05);
401  t.DrawLatex(0.3, 0.9, text1);
402 
403  keyLine(.4, 0.80,"MC tracks",1);
404  keyLine(.4, 0.75,"MatGlobal RC tracks "+tagTex,2);
405 
406  cpt->Update();
407 
408  if (plot)
409  {
410  sprintf(gif,"~/embedding/%s/%s/pT_%s_%s_%s_%s.gif",Prod, Particle, Prod, Particle,Daughter, Date);
411  cpt->SaveAs(gif);
412  }
413  // return;
414 
415  //======================================================
416  //phi
417  //==============================================================
418 
419 
420  TCanvas *cphi= new TCanvas("cphi","PHI",500, 500);
421  cphi->SetGridx(0);
422  cphi->SetGridy(0);
423 
424 
425  cphi->cd();
426 
427  PhiMc1->Rebin(2);
428  PhiMc1->SetLineColor(1);
429  PhiMc1->Scale(1/2);
430  PhiMc1->Scale(PhiMc1->GetSum());
431 
432  PhiMc1->SetTitle("Phi");
433  PhiMc1->SetXTitle ("Phi");
434  PhiMc1->Draw();
435 
436  //Reconstructed
437  PhiM1->Rebin(2);
438  PhiM1->Scale(1/2);
439  PhiM1->Scale(PhiM1->GetSum());
440 
441  PhiMc1->SetAxisRange(-3.3, 3.3, "X");
442  PhiMc1->SetAxisRange(0, 1.2*PhiMc1->GetMaximum(), "Y");
443 
444  PhiM1->SetLineColor(2);
445  PhiM1->Draw("same");
446 
447  PhiG1->Rebin(2);
448  PhiG1->Scale(1/2);
449 
450  PhiG1->Scale(PhiMc1->GetSum()/PhiG1->GetSum()*0.4);//this is just for dAu run 8
451  PhiG1->SetLineColor(3);
452  PhiG1->Draw("same");
453 
454  t.SetTextSize(0.05);
455  t.DrawLatex(0.3, 0.9, text1);
456 
457  keyLine(.3, 0.18,"MC tracks",1);
458  keyLine(.3, 0.13,"MatGlobal RC tracks "+tagTex,2);
459  keyLine(.3, 0.10,"Ghost tracks (rescaled)",3);
460 
461  cphi->Update();
462  if (plot)
463  {
464  sprintf(gif,"~/embedding/%s/%s/Phi_%s_%s_%s_%s.gif",Prod, Particle, Prod, Particle,Daughter, Date);
465  cphi->SaveAs(gif);
466  }
467  // return;
468 
469  //========================================================
470  //rapidity
471  //======================================================
472 
473  TCanvas *cy= new TCanvas("cy","Rapidity",500, 500);
474 
475  cy->SetGridx(0);
476  cy->SetGridy(0);
477  cy->cd();
478 
479  //embedded
480 
481 
482  YMc1->SetLineColor(1);
483  YMc1->SetTitle("Rapidity");
484  YMc1->SetXTitle ("y");
485  yMc1->SetAxisRange(-1.2, 1.2, "X");
486  YMc1->SetAxisRange(0, 1.4* YMc1->GetMaximum(), "Y");
487  YMc1->Draw();
488 
489  //Reco
490 
491  YM1->SetLineColor(2);
492  YM1->Draw("same");
493 
494 
495  t.SetTextSize(0.05);
496  t.DrawLatex(0.3, 0.9, text1);
497 
498  keyLine(.3, 0.21,"MC tracks",1);
499  keyLine(.3, 0.15,"MatGlobal RC tracks "+tagTex,2);
500 
501 
502  cy->Update();
503 
504  if (plot)
505  {
506  sprintf(gif,"~/embedding/%s/%s/Rapidity_%s_%s_%s_%s.gif",Prod, Particle, Prod, Particle,Daughter, Date);
507  cy->SaveAs(gif);
508  }
509  // return;
510 
511  */
512  //============================================================
513  //Dca vs Eta vs pt
514  //============================================================
515 
516  TCanvas *cde= new TCanvas("cde","DCA_Eta",900, 900);
517 
518  cde->Divide(nPtBin,nEtaBin);//etabin
519  cde->SetGridx(0);
520  cde->SetGridy(0);
521 
522  //Matched (Bins for Multiplicity)
523 
524  TH1D *hpt = (TH1D*)hDcaM->Project3D("X");
525 
526 
527  for (Int_t ipt=0; ipt <nPtBin ; ipt++)
528  {
529  Int_t bin_pt1 = hpt->FindBin(pt[ipt]);
530  Int_t bin_pt2 = hpt->FindBin(pt[ipt+1]);//this should be the same for both graphs (for 3 graphs)
531 
532 
533  TString nameM = "hDcaM";
534  TString nameR = "hDcaR";
535  TString nameG = "hDcaG";
536 
537  TH1D *hEtaM = (TH1D*)hDcaM->Project3D("Y");
538  TH1D *hEtaR = (TH1D*)hDcaR->Project3D("Y");
539  TH1D *hEtaG = (TH1D*)hDcaG->Project3D("Y");
540 
541 
542  for(int ieta=0; ieta<nEtaBin ; ieta++)
543  {
544  cde->cd(3*ipt+ieta+1);
545 
546  Double_t SumMC=0;
547  Double_t SumR=0;
548  Double_t SumG=0;
549 
550  Int_t bin_eta1 = hEtaM->FindBin(eta[ieta]);
551  Int_t bin_eta2 = hEtaM->FindBin(eta[ieta+1]);
552 
553  TH1D *hDcaMNew= (TH1D*)hDcaM->ProjectionZ(nameM+ieta,bin_pt1, bin_pt2, bin_eta1, bin_eta2);
554 
555  // hDcaMNew->Sumw2();
556 
557  // maxMc = hDcaMNew->GetMaximum();
558  //SumMC = hDcaMNew->GetSum();
559 
560  for(int kkk=1;kkk<hDcaMNew->GetNbinsX();kkk++)
561  {
562  SumMC += hDcaMNew->GetBinContent(kkk);
563  }
564 
565  hDcaMNew->Scale(1./SumMC);
566  hDcaMNew->SetLineColor(2);
567  hDcaMNew->SetXTitle("Dca (cm)");
568 
569  sprintf(title," %.2f GeV < pT < %.2f GeV - %.2f GeV/c < Eta < %.2f GeV/c", pt[ipt],pt[ipt+1], eta[ieta],eta[ieta+1]);
570  hDcaMNew->SetTitle(title);
571 
572  //----Now MuDSt
573 
574  Int_t bin_eta1R = hEtaR->FindBin(eta[ieta]);
575  Int_t bin_eta2R = hEtaR->FindBin(eta[ieta+1]);
576 
577  TH1D *hDcaRNew= (TH1D*)hDcaR->ProjectionZ(nameR+ieta,bin_pt1, bin_pt2, bin_eta1R, bin_eta2R);
578  // hDcaRNew->Sumw2();
579 
580  // SumR = hDcaRNew->GetSum();
581  for(int kkk=1;kkk<hDcaRNew->GetNbinsX();kkk++)
582  {
583  SumR += hDcaRNew->GetBinContent(kkk);
584 
585  }
586 
587  hDcaRNew->Scale(1./SumR);
588 
589 
590  //---Now Ghost
591 
592  Int_t bin_eta1G = hEtaG->FindBin(eta[ieta]);
593  Int_t bin_eta2G = hEtaG->FindBin(eta[ieta+1]);
594 
595  TH1D *hDcaGNew = (TH1D*)hDcaG->ProjectionZ(nameG+ieta, bin_pt1, bin_pt2, bin_eta1G, bin_eta2G);
596  // hDcaGNew->Sumw2();
597 
598  // SumG = hDcaGNew->GetSum();
599 
600  for(int kkk=1; kkk<hDcaGNew->GetNbinsX();kkk++)
601  {
602  SumG += hDcaGNew->GetBinContent(kkk);
603 
604  }
605 
606 
607  hDcaGNew->Scale(1./SumG);
608  hDcaGNew->SetLineColor(4);
609 
610  hDcaMNew->Draw();
611  hDcaRNew->Draw("esame");
612  // hDcaGNew->Draw("esame");
613 
614 
615  t.SetTextSize(0.045);
616  t.DrawLatex(0.4, 0.9, text1);
617 
618  t.SetTextSize(0.045);
619  t.SetTextColor(2);
620  t.DrawLatex(0.4, 0.8, "Matched Global "+tagTex);
621 
622  t.SetTextSize(0.045);
623  t.SetTextColor(1);
624  t.DrawLatex(0.4, 0.75, "MuDst");
625 
626  t.SetTextSize(0.045);
627  t.SetTextColor(1);
628  t.DrawLatex(0.4, 0.70, "Ghost Tracks");
629 
630  cde->Update();
631 
632  }
633 
634 
635  }
636  cde->Update();
637  // sprintf(gif," ~/embedding/%s/%s/dca_eta_pt_%s_%s_%s_%s.gif",Prod, Particle, Prod, Particle, Daughter,Date);
638  // cde->SaveAs(gif);
639  return;
640 
641  //===============================================
642  //nfit vs eta vs pt
643  //====================================================
644 
645  TCanvas *cn= new TCanvas("cn","NFIT",1000, 800);
646  cn->Divide(nPtBin,nEtaBin);
647  cn->SetGridx(0);
648  cn->SetGridy(0);
649 
650  //Bins for Multiplicity -Matched tracks
651 
652  TH1D *hpt = (TH1D*)hNfitM->Project3D("X");
653 
654  for (Int_t ipt=0; ipt <nPtBin ; ipt++)
655  {
656 
657  Int_t bin_pt1 = hpt->FindBin(pt[ipt]);
658  Int_t bin_pt2 = hpt->FindBin(pt[ipt+1]);//this should be the same for both graphs (for 3 graphs)
659 
660  TString nameM = "hNfitM";
661  TString nameR = "hNfitR";
662  TString nameG = "hNfitG";
663 
664  TH1D *hEtaM = (TH1D*)hNfitM->Project3D("Y");
665  TH1D *hEtaR = (TH1D*)hNfitR->Project3D("Y");
666  TH1D *hEtaG = (TH1D*)hNfitG->Project3D("Y");
667 
668 
669  for(int ieta=0; ieta<nEtaBin ; ieta++)
670 
671  {
672  cn->cd(3*ipt+ieta+1);
673 
674  Double_t SumM=0;
675  Double_t SumR=0;
676  Double_t SumG=0;
677 
678  Int_t bin_eta1 = hEtaM->FindBin(eta[ieta]);
679  Int_t bin_eta2 = hEtaM->FindBin(eta[ieta+1]);
680 
681 
682  TH1D *hNfitMNew= (TH1D*)hNfitM->ProjectionZ(nameM+ieta,bin_pt1, bin_pt2, bin_eta1, bin_eta2);
683  hNfitMNew->Sumw2();
684 
685  for(int kkk=1;kkk<hNfitMNew->GetNbinsX();kkk++)SumM += hNfitMNew->GetBinContent(kkk);
686  cout<<SumM<<endl;
687 
688  hNfitMNew->Scale(1./SumM);
689  hNfitMNew ->SetLineColor(2);
690  hNfitMNew->Draw();
691 
692  hNfitMNew->SetXTitle("Nfit");
693 
694  sprintf(title," %.2f GeV < pT < %.2f GeV - %.2f GeV/c < Eta < %.2f GeV/c", pt[ipt],pt[ipt+1], eta[ieta],eta[ieta+1]);
695  hNfitMNew->SetTitle(title);
696 
697  //----Now MuDSt
698 
699  Int_t bin_eta1R = hEtaR->FindBin(eta[ieta]);
700  Int_t bin_eta2R = hEtaR->FindBin(eta[ieta+1]);
701 
702  TH1D *hNfitRNew= (TH1D*)hNfitR->ProjectionZ(nameR+ieta,bin_pt1, bin_pt2, bin_eta1R, bin_eta2R);
703 
704  hNfitRNew->Sumw2();
705  SumR = hNfitRNew->GetSum();
706 
707  for(int kkk=1;kkk<hNfitR->GetNbinsX();kkk++)SumR += hNfitR->GetBinContent(kkk);
708  cout<<SumR<<endl;
709  hNfitRNew->Scale(1./SumR);
710  hNfitRNew->Draw("esame");
711 
712 
713  //Now Ghost
714 
715  Int_t bin_eta1G = hEtaG->FindBin(eta[ieta]);
716  Int_t bin_eta2G = hEtaG->FindBin(eta[ieta+1]);
717 
718 
719  TH1D *hNfitGNew = (TH1D*)hNfitG->ProjectionZ(nameG+ieta,bin_pt1, bin_pt2, bin_eta1, bin_eta2);
720  hNfitGNew->Sumw2();
721 
722  SumG = hNfitRNew->GetSum();
723  for(int kkk=1;kkk<hNfitGNew->GetNbinsX();kkk++)SumG += hNfitGNew->GetBinContent(kkk);
724  cout<<SumG<<endl;
725 
726  hNfitGNew->Scale(1./SumG);
727  hNfitGNew ->SetLineColor(4);
728  // hNfitGNew->Draw("same");
729 
730  t.SetTextSize(0.045);
731  t.DrawLatex(0.4, 0.9, text1);
732 
733  t.SetTextSize(0.045);
734  t.SetTextColor(2);
735  t.DrawLatex(0.4, 0.8, "Matched Global "+tagTex);
736 
737  t.SetTextSize(0.045);
738  t.SetTextColor(1);
739  t.DrawLatex(0.4, 0.75, "MuDst");
740 
741  cn->Update();
742  // keyLine(0.2, 0.80,"Ghost Tracks",4);
743  }
744 
745  }
746  cn->Update();
747  sprintf(gif,"~/embedding/%s/%s/Nfit_eta_pt_%s_%s_%s_%s.gif",Prod, Particle, Prod, Particle,Daughter, Date);
748  cn->SaveAs(gif);
749 
750  return;
751 }
752 
753 
754  //=====================================================================
755  //Following you can use them if need them.
756  // Inludes, MIPS, Efficiency, Energy Lost, Dca (just Pt dependant), and NFit ( just pt dependant)//
757 
758 /*
759  //==================================================================================
760  //MIPS (just for pions)
761  //=============================================================
762 
763  if (id==8 || id==9)
764  {
765 
766  TCanvas *c9= new TCanvas("c9","MIPS",500, 500);
767 
768  c9->SetGridx(0);
769  c9->SetGridy(0);
770  c9->SetLeftMargin(0.15);
771  c9->SetRightMargin(0.05);
772  c9->cd();
773 
774  double pt1 = 0.4;
775  double pt2 = 0.6;
776 
777  dedxG1 -> ProjectionX("rpx");
778 
779  int blG = rpx->FindBin(pt1);
780  int bhG = rpx->FindBin(pt2);
781 
782  cout<<blG<<endl;
783  cout<<bhG<<endl;
784 
785  dedxG1->ProjectionY("rpy",blG,bhG);
786  rpy->SetTitle("MIPS");
787  rpy->SetMarkerStyle(22);
788  // rpy->SetMarkerColor(2);
789  rpy->SetAxisRange(1.3, 4, "X");
790 
791  //dedxG1->Draw();
792 
793  dedx1->ProjectionX("mpx");
794 
795  int blm = mpx->FindBin(pt1);
796  int bhm = mpx->FindBin(pt2);
797 
798  cout<<blm<<endl;
799  cout<<bhm<<endl;
800 
801  dedx1->ProjectionY("mpy", blm,bhm);
802 
803  mpy->SetAxisRange(0.5, 6, "X");
804  mpy->SetMarkerStyle(22);
805  mpy->SetMarkerColor(2);
806 
807  float max_rpy = rpy->GetMaximum();
808  max_rpy /= 1.*mpy->GetMaximum();
809  mpy->Scale(max_rpy);
810 
811  cout<<"max_rpy is: "<<max_rpy<<endl;
812  cout<<"mpy is: "<<mpy<<endl;
813 
814  rpy->Sumw2();
815  mpy->Sumw2();
816 
817  rpy->Fit("gaus","","",1,4);
818  mpy->Fit("gaus","","", 1, 4);
819  mpy->GetFunction("gaus")->SetLineColor(2);
820 
821  rpy->SetAxisRange(0.5 ,6.0, "X");
822  mpy->Draw();
823  rpy->Draw("same");
824 
825  float mipMc = mpy->GetFunction("gaus")->GetParameter(1);//mean
826  float mipGhost = rpy->GetFunction("gaus")->GetParameter(1);
827 
828  float sigmaMc = mpy->GetFunction("gaus")->GetParameter(2);//mean
829  float sigmaGhost = rpy->GetFunction("gaus")->GetParameter(2);
830 
831  char label1[80];
832  char label2[80];
833  char label3[80];
834  char label4[80];
835 
836 
837  sprintf(label1,"mip MC %.3f",mipMc);
838  sprintf(label2,"mip Ghost %.3f",mipGhost);
839  sprintf(label3,"sigma MC %.3f",sigmaMc);
840  sprintf(label4,"sigma Ghost %.3f",sigmaGhost);
841 
842  keySymbol(.5, .9, label1,2,1);
843  keySymbol(.5, .85, label3,2,1);
844 
845  keySymbol(.5, .75, label2,1,1);
846  keySymbol(.5, .70, label4,1,1);
847 
848  keyLine(.2, 1.05,text1,1);
849 
850  char name[30];
851  sprintf(name,"%.2f GeV/c < Pt < %.2f GeV/c",pt1, pt2);
852  keySymbol(0.3, 0.65, name, 1, 1, 0.04);
853 
854 
855  c9->Update();
856 
857  }//close if pion
858 
859 
860  // pfx1->SetMarkerStyle(2,1);
861 
862 
863 
864  //================================================================
865  //DCA (pt)
866  //==================================================================
867  /*
868 
869  TCanvas *c= new TCanvas("c","DCA",400, 800);
870  c->Divide(1,nPtBin);
871  c->SetGridx(0);
872  c->SetGridy(0);
873 
874  //Matched (Bins for Multiplicity)
875 
876  TH1D *hX1 = (TH1D*)hDca1->Project3D("X");
877 
878  Int_t bin_nch1 = hX1->FindBin(nch1);
879  Int_t bin_nch2 = hX1->FindBin(nch2);//this should be the same for both graphs (for 3 graphs)
880 
881  cout<<bin_nch1<<" "<<bin_nch2<<endl;
882 
883  //Bins for Pt
884 
885  TString name1 = "hDca1";
886  TString namer1 = "hDcar1";
887 
888  TString namexx = "hDca";
889 
890 
891  TH1D *hY1 = (TH1D*)hDca1->Project3D("Y");
892  TH1D *hY1_r = (TH1D*)hDca1r->Project3D("Y");
893 
894  TH1D *hYG = (TH1D*)hDcaG1->Project3D("Y");
895 
896  for(Int_t i=0; i<nPtBin ; i++)
897  {
898  c->cd(i+1);
899 
900  Double_t Sum1_MC=0;
901  Double_t Sum1_Real=0;
902  Double_t Sum_MC=0;
903 
904  Int_t bin_ptl_1 = hY1->FindBin(pt[i]);
905  Int_t bin_pth_1 = hY1->FindBin(pt[i+1]);
906 
907  TH1D *hDcaNew1= (TH1D*)hDca1->ProjectionZ(name1+i,bin_nch1, bin_nch2, bin_ptl_1, bin_pth_1);
908  //Sum1_MC = hDcaNew1 ->GetSum();
909  for(int kkk=1;kkk<hDcaNew1->GetNbinsX();kkk++)Sum1_MC += hDcaNew1->GetBinContent(kkk);
910  cout<<Sum1_MC<<endl;
911 
912  hDcaNew1->Scale(1./Sum1_MC);
913  hDcaNew1 ->SetLineColor(2);
914  hDcaNew1->Draw();
915 
916  hDcaNew1->SetXTitle("Dca (cm)");
917 
918  //sprintf(title," %.2f GeV < pT < %.2f GeV, %d < nch < %d", pt[i],pt[i+1],nch1,nch2);
919  sprintf(title," %.2f GeV/c < pT < %.2f GeV/c", pt[i],pt[i+1]);
920  //hDcaNew1->SetTitle("Scale factor = 1.38");
921  hDcaNew1->SetTitle(title);
922 
923  //----Now MuDSt
924 
925  Int_t bin_ptrl_1r = hY1_r->FindBin(pt[i]);
926  Int_t bin_ptrh_1r = hY1_r->FindBin(pt[i+1]);
927 
928  TH1D *hDca_r1= (TH1D*)hDca1r->ProjectionZ(namer1+i,bin_nch1, bin_nch2, bin_ptrl_1r, bin_ptrh_1r);
929  //Sum1_Real = hDca_r1 ->GetSum();
930  for(int kkk=1;kkk<hDca_r1->GetNbinsX();kkk++)Sum1_Real += hDca_r1->GetBinContent(kkk);
931  cout<<Sum1_Real<<endl;
932 
933  hDca_r1->Scale(1./Sum1_Real);
934  //hDca_r1->Draw("same");
935 
936  //Now Previous Embedding
937 
938  Int_t bin_ptl = hYG->FindBin(pt[i]);
939  Int_t bin_pth = hYG->FindBin(pt[i+1]);
940 
941  TH1D *hDcaNew = (TH1D*)hDcaG1->ProjectionZ(namexx+i,bin_nch1, bin_nch2, bin_ptl, bin_pth);
942  //Sum_MC = hDcaNew ->GetSum();
943  for(int kkk=1;kkk<hDcaNew->GetNbinsX();kkk++)Sum_MC += hDcaNew->GetBinContent(kkk);
944  cout<<Sum_MC<<endl;
945 
946  hDcaNew->Scale(1./Sum_MC);
947  hDcaNew ->SetLineColor(4);
948  hDcaNew->Draw("same");
949 
950  //keySymbol(.4, .95,text1,1,1);
951  keyLine(0.4, 0.90,"MatGlobal Tracks",2);
952  //keyLine(0.4, 0.85,"MuDst",1);
953  keyLine(0.4, 0.80,"Ghost Tracks",4);
954  }
955 
956 
957  c->Update();
958  c->SaveAs("/home/mcsuarez/Embed/P08id/dca_f1.38.gif");
959 
960 
961 
962 
963  //===================================================
964  //NFIT
965  //=============================================================
966 
967  TCanvas *c1= new TCanvas("c1","NFIT",400, 800);
968  c1->Divide(1,nPtBin);
969  c1->SetGridx(0);
970  c1->SetGridy(0);
971 
972  //Bins for Multiplicity -Matched tracks
973 
974  TH1D *hX1 = (TH1D*)hNfit1->Project3D("X");
975 
976  Int_t bin_nch1 = hX1->FindBin(nch1);
977  Int_t bin_nch2 = hX1->FindBin(nch2);//this should be the same for both graphs (for 3 graphs)
978 
979  //Bins for Pt
980 
981  TString name_nfit1 = "hNfit1";
982  TString name_nfitr1 = "hNfitr1";
983  TString name_nfit = "hNfit";
984 
985  TH1D *hY1 = (TH1D*)hNfit1->Project3D("Y");
986  TH1D *hY1_r = (TH1D*)hNfit1r->Project3D("Y");
987  TH1D *hYG = (TH1D*)hNfitG1->Project3D("Y");
988 
989  for(Int_t i=0; i<nPtBin ; i++)
990 
991  {
992  c1->cd(i+1);
993 
994  Double_t Sum1_Nfit_MC = 0;
995  Double_t Sum1_Nfit_Real = 0;
996  Double_t Sum_Nfit_MC = 0;
997 
998  Int_t bin_ptl_1 = hY1->FindBin(pt[i]);
999  Int_t bin_pth_1 = hY1->FindBin(pt[i+1]);
1000 
1001  TH1D *hNfitNew1= (TH1D*)hNfit1->ProjectionZ(name_nfit1+i,bin_nch1, bin_nch2, bin_ptl_1, bin_pth_1);
1002  //Sum1_Nfit_MC = hNfitNew1 ->GetSum();
1003  for(int kkk=1;kkk<hNfitNew1->GetNbinsX();kkk++)Sum1_Nfit_MC += hNfitNew1->GetBinContent(kkk);
1004  cout<<Sum1_Nfit_MC<<endl;
1005 
1006  hNfitNew1->Scale(1./Sum1_Nfit_MC);
1007  hNfitNew1 ->SetLineColor(2);
1008  hNfitNew1->Draw();
1009 
1010  hNfitNew1->SetXTitle("Nfit");
1011 
1012  //sprintf(title," %.2f GeV < pT < %.2f GeV, %d < nch < %d", pt[i],pt[i+1],nch1,nch2);
1013  sprintf(title," %.2f GeV/c < pT < %.2f GeV/c", pt[i],pt[i+1],nch1,nch2);
1014  //hNfitNew1->SetTitle("Scale factor = 1.38");
1015  hNfitNew1->SetTitle(title);
1016 
1017  //----Now MuDSt
1018 
1019  Int_t bin_ptrl_1r = hY1_r->FindBin(pt[i]);
1020  Int_t bin_ptrh_1r = hY1_r->FindBin(pt[i+1]);
1021 
1022  TH1D *hNfit_r1= (TH1D*)hNfit1r->ProjectionZ(name_nfitr1+i,bin_nch1, bin_nch2, bin_ptrl_1r, bin_ptrh_1r);
1023  //Sum1_Nfit_Real = hNfit_r1 ->GetSum();
1024 
1025  for(int kkk=1;kkk<hNfit_r1->GetNbinsX();kkk++)Sum1_Nfit_Real += hNfit_r1->GetBinContent(kkk);
1026  cout<<Sum1_Nfit_Real<<endl;
1027  hNfit_r1->Scale(1./Sum1_Nfit_Real);
1028  //hNfit_r1->Draw("same");
1029 
1030  //Now Previous Embedding
1031 
1032  Int_t bin_ptl = hYG->FindBin(pt[i]);
1033  Int_t bin_pth = hYG->FindBin(pt[i+1]);
1034 
1035  TH1D *hNfitNew = (TH1D*)hNfitG1->ProjectionZ(name_nfit+i,bin_nch1, bin_nch2, bin_ptl, bin_pth);
1036  //Sum_Nfit_MC = hNfitNew ->GetSum();
1037  for(int kkk=1;kkk<hNfitNew->GetNbinsX();kkk++)Sum_Nfit_MC += hNfitNew->GetBinContent(kkk);
1038  cout<<Sum_Nfit_MC<<endl;
1039 
1040  hNfitNew->Scale(1./Sum_Nfit_MC);
1041  hNfitNew ->SetLineColor(4);
1042  hNfitNew->Draw("same");
1043 
1044  //keySymbol(.2, .95,text1,1,1);
1045  keyLine(0.2, 0.90,"MatGlobal Tracks",2);
1046  //keyLine(0.2, 0.85,"MuDst",1);
1047  keyLine(0.2, 0.80,"Ghost Tracks",4);
1048  }
1049 
1050  c1->Update();
1051  c1->SaveAs("/home/mcsuarez/Embed/P08id/nfit_f1.38.gif");
1052 
1053  return;
1054 
1055 
1056 
1057  // ===============================================================
1058  //Efficiency vs pT
1059  //this make sense just when reconstructed particels is the same as embedded (?)
1060 
1061  void PlotEfficiency(const TH1* hData, const TH1* hMC)
1062  {
1063  TCanvas *c1= new TCanvas(Form("c%d", iCanvas), Form("c%d", iCanvas));
1064  hData->Rebin(2);
1065  hMC->Rebin(2);
1066 
1067  TH1* hRatio = (TH1D*) hData->Clone();
1068  hRatio->Divide(hMC);
1069  hRatio->SetLineColor(1);
1070  hRatio->SetMarkerStyle(23);
1071  hRatio->SetMarkerColor(1);
1072  hRatio->SetXTitle ("p_{T} (GeV/c)");
1073  hRatio->SetAxisRange(0.0, 6.0, "X");
1074 
1075  hRatio->Draw();
1076 
1077  c1->cd();
1078  c1->Update();
1079  c1->SaveAs("eff.gif");
1080 
1081  // PtM2->Rebin(2);
1082  // PtMc2->Rebin(2);
1083  //
1084  // PtM2->Divide(PtMc2);
1085  // PtM2->SetLineColor(9);
1086  // PtM2->SetMarkerStyle(21);
1087  // PtM2->SetMarkerColor(9);
1088  // PtM2->Draw("same");
1089  //
1090  // PtM3->Rebin(2);
1091  // PtMc3->Rebin(2);
1092  //
1093  // PtM3->Divide(PtMc2);
1094  // PtM3->SetLineColor(2);
1095  // PtM3->SetMarkerStyle(22);
1096  // PtM3->SetMarkerColor(2);
1097  // PtM3->Draw("same");
1098 
1099 
1100  //============================================================
1101  // Vertex xx yy zz
1102  //Not too much sense
1103 
1104  TCanvas *cv1= new TCanvas("cv1","Vertex position",600, 400);
1105 
1106  cv1->Divide(3,1);
1107 
1108  cv1_1->cd();
1109 
1110  vxx1->SetXTitle("MonteCarlo Vertex X");
1111  vxx1->SetYTitle("Matched Vertex X");
1112 
1113  vxx1->Draw();
1114 
1115  cv1_2->cd();
1116 
1117  vyy1->SetXTitle("MonteCarlo Vertex Y");
1118  vyy1->SetYTitle("Matched Vertex Y");
1119  vyy1->Draw();
1120 
1121  cv1_3->cd();
1122 
1123  vzz1->SetXTitle("MonteCarlo Vertex Z");
1124  vzz1->SetYTitle("Matched Vertex Z");
1125  vzz1->Draw();
1126 
1127  cv1->Update();
1128  cv1->SaveAs("/home/mcsuarez/Embed/P08id/Vertex_xyz.gif");
1129 
1130 
1131 
1132  //========================================================================
1133  //Energy loss
1134  //=================================================================
1135 
1136  TCanvas *c7= new TCanvas("c7","Energy Loss",400, 400);
1137 
1138  c7->SetGridx(0);
1139  c7->SetGridy(0);
1140  c7->SetLeftMargin(0.20);
1141  c7->SetRightMargin(0.05);
1142  c7->cd();
1143 
1144  hPtM_E1->ProfileX("pfx");
1145  pfx->SetAxisRange(-0.01, 0.01, "Y");
1146  pfx->SetAxisRange(0, 2.5, "X");
1147  pfx->GetYaxis()->SetDecimals();
1148  pfx->SetMarkerStyle(23);
1149  pfx->SetMarkerSize(0.038);
1150  pfx->SetMarkerColor(4);
1151  pfx->SetLineColor(4);
1152  pfx->SetXTitle ("Pt-Reco");
1153  pfx->SetYTitle ("ptM - PtMc");
1154  pfx->SetTitleOffset(2,"Y");
1155 
1156  pfx->Draw();
1157 
1158 
1159  hPtM_E1->ProfileX("pfx1");
1160  pfx1->SetAxisRange(-0.007, 0.007, "Y");
1161  pfx1->GetYaxis()->SetDecimals();
1162  pfx1->SetLineColor(2);
1163  pfx1->SetMarkerSize(0.035);
1164  pfx1->SetXTitle ("Pt-Reco");
1165  pfx1->SetYTitle ("ptM - PtMc");
1166  pfx1->SetTitleOffset(2,"Y");
1167 
1168  pfx1->Draw("same");
1169 
1170 
1171  c7->Update();
1172  c7->SaveAs("/home/mcsuarez/Embed/P08id/eloss_f1.38.gif");
1173 
1174 
1175  //====================================================
1176  //Dedx (p08ic)
1177  //=======================================================
1178 
1179 
1180  TCanvas *c88= new TCanvas("c88","dEdx vs P (P08ic)", 600, 600);
1181 
1182  c88->SetGridx(0);
1183  c88->SetGridy(0);
1184  c88->SetLeftMargin(0.15);
1185  c88->SetRightMargin(0.05);
1186  c88->cd();
1187  gStyle->SetOptDate(0);
1188 
1189  dedx1R->SetTitle("dEdx vs P (P08ic)");
1190  dedx1R->SetXTitle("Momentum P (GeV/c)");
1191  dedx1R->SetYTitle("dE/dx (x10^{-6})");
1192  dedx1R->SetAxisRange(0, 5., "X");
1193  dedx1R->SetAxisRange(0, 8., "Y");
1194  dedx1R->SetMarkerColor(1);
1195  dedx1R->Draw();//"colz");
1196 
1197  c88->Update();
1198  c88->SaveAs("dedx_p08ic.gif");
1199  */