StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StEmbeddingQADraw.cxx
1 /****************************************************************************************************
2  * $Id: StEmbeddingQADraw.cxx,v 1.36 2019/07/10 05:45:34 zhux Exp $
3  * $Log: StEmbeddingQADraw.cxx,v $
4  * Revision 1.36 2019/07/10 05:45:34 zhux
5  * added option for btof pid for primary real tracks
6  *
7  * Revision 1.35 2015/04/19 01:46:59 zhux
8  * Wrong particle name (under sl64) is fixed in ::getParticleName().
9  *
10  * Revision 1.34 2013/03/20 01:08:39 huck
11  * use mu- to compare to pi-
12  *
13  * Revision 1.33 2011/04/26 20:27:45 hmasui
14  * Add gamma geantid check
15  *
16  * Revision 1.32 2011/04/12 03:01:04 hmasui
17  * Fix isMatchedPairOk() to properly process particles with decay daughters
18  *
19  * Revision 1.31 2011/04/07 02:29:57 hmasui
20  * Print (pseudo-)rapidity up to 2nd digit string
21  *
22  * Revision 1.30 2011/04/06 20:01:29 hmasui
23  * Put back Ncommon vs NHitFit histogram for all pt, and fix a bug for name of projected histograms for Ncommon vs NHitFit
24  *
25  * Revision 1.29 2011/04/05 04:26:42 hmasui
26  * Bug fix for printing trigger id's, and follow default ROOT color scheme
27  *
28  * Revision 1.28 2011/04/01 05:07:07 hmasui
29  * Added track selections in the 2nd page, Ncommon vs NhitFit (pt dependent), and 1/pt(RC)-1/pt(MC) vs pt plot
30  *
31  * Revision 1.27 2011/02/25 18:34:40 hmasui
32  * Add phi comparison between real and reconstructed primary (matched pairs) in embedding
33  *
34  * Revision 1.26 2011/02/11 03:44:53 hmasui
35  * Draw error messages in pdf if histogram is missing. Add error check for Ncommon histogram
36  *
37  * Revision 1.25 2011/01/31 21:33:53 hmasui
38  * Add setParentGeantId() function to allow the multiple decays
39  *
40  * Revision 1.24 2010/08/04 21:16:06 hmasui
41  * Replace geantid to nsigma cut for real tracks in legend
42  *
43  * Revision 1.23 2010/08/03 23:39:35 hmasui
44  * Fix a bug in SetMaximum() function for MC rapidity distribution
45  *
46  * Revision 1.22 2010/07/12 21:31:49 hmasui
47  * Use StEmbeddingQAUtilities::getParticleDefinition() instead of StParticleTable
48  *
49  * Revision 1.21 2010/06/28 21:29:41 hmasui
50  * Fixed the multiplicity x-axis range for all branches
51  *
52  * Revision 1.20 2010/06/22 16:31:21 hmasui
53  * Separate 2D and 1D QA for MC tracks. Add pol0 fit for MC eta, y and phi distributions.
54  *
55  * Revision 1.19 2010/06/10 14:51:25 hmasui
56  * Added legend for each page
57  *
58  * Revision 1.18 2010/05/27 16:29:00 hmasui
59  * Remove / character from particle name
60  *
61  * Revision 1.17 2010/05/14 19:51:45 hmasui
62  * Modify the text size for Nevts, MC particle name etc to fit the window
63  *
64  * Revision 1.16 2010/04/24 19:50:56 hmasui
65  * Optimize to draw run number list, and fix bugs for maximum of y-axis in several histograms
66  *
67  * Revision 1.15 2010/04/07 19:45:12 hmasui
68  * Use box option for dE/dx vs p to reduce the pdf file size
69  *
70  * Revision 1.14 2010/03/15 21:05:24 hmasui
71  * Separate MC vertices QA into 2 pages. Added constraint on z-vertex cut for vx(vy) vs vz histograms.
72  *
73  * Revision 1.13 2010/02/24 18:13:23 hmasui
74  * Modify color code explanation more explicitly. Comparison of phi distributions between reconstructed embedding with MC tracks
75  *
76  * Revision 1.12 2010/02/23 16:56:47 hmasui
77  * Add phi distributions QA (MC vs reconstructed)
78  *
79  * Revision 1.11 2010/02/19 18:07:45 hmasui
80  * Divide runlist into two different pads
81  *
82  * Revision 1.10 2010/02/16 02:14:03 hmasui
83  * Print PDF file only for all QA plots
84  *
85  * Revision 1.9 2010/01/28 21:51:24 hmasui
86  * Add Vx vs Vz and Vy vs Vz histograms in the event-wise QA
87  *
88  * Revision 1.8 2010/01/26 17:50:54 hmasui
89  * Fix geantidFound to match the correct geant id. Use number of accepted events, not number of all events without vertex cuts. Add QA for eventid, runid and # of particles per event
90  *
91  * Revision 1.7 2009/12/22 21:40:09 hmasui
92  * Add comments for functions and members
93  *
94  ****************************************************************************************************/
95 
96 #include <assert.h>
97 #include <string>
98 
99 #include "TCanvas.h"
100 #include "TError.h"
101 #include "TF1.h"
102 #include "TFile.h"
103 #include "TGraphErrors.h"
104 #include "TH1.h"
105 #include "TH2.h"
106 #include "TH3.h"
107 #include "TLatex.h"
108 #include "TLegend.h"
109 #include "TMath.h"
110 #include "TObject.h"
111 #include "TPaveText.h"
112 #include "TPDF.h"
113 #include "TProfile.h"
114 #include "TROOT.h"
115 #include "TStyle.h"
116 #include "TSystem.h"
117 #include "TLine.h"
118 
119 #include "StEmbeddingQADraw.h"
120 #include "StEmbeddingQAUtilities.h"
121 #include "StMessMgr.h"
122 #include "StParticleDefinition.hh"
123 
124 using namespace std ;
125 
126 ClassImp(StEmbeddingQADraw)
127 
128 
129  UInt_t StEmbeddingQADraw::mCanvasId = 0;
130 
131 //____________________________________________________________________________________________________
132 StEmbeddingQADraw::StEmbeddingQADraw(const TString embeddingFile, const TString realDataFile, const Int_t geantid,
133  const Bool_t isEmbeddingOnly)
134  : mIsEmbeddingOnly(isEmbeddingOnly), mIsGIFOn(kFALSE), mIsJPGOn(kFALSE), mIsEPSOn(kFALSE), mIsPSOn(kFALSE), mGeantId(geantid)
135 {
137  mCanvas = 0 ;
138  mMainPad = 0 ;
139  mPDF = 0 ;
140 
141  mInputEmbedding = 0 ;
142  mInputRealData = 0 ;
143 
145  mIsOpen = open(embeddingFile, realDataFile);
146 
148  if(mIsOpen){
149  TString fileName(embeddingFile);
150  fileName.Remove(fileName.Index(".root"), fileName.Length()); // remove .root
151 
152  for(Int_t i=0;i<3;i++){
153  const Int_t start = 0 ;
154  const Int_t stop = fileName.First("_") ;
155  const TString subString(fileName(start, stop));
156 
157  fileName.Remove(start, stop+1);
158 
159  if( i == 2 ){
160  mYear = subString.Atoi(); // Convert character to integer
161  mProduction = fileName ;
162  }
163  }
164  }
165  else{
167  mYear = -10 ;
168  mProduction = "";
169  }
170 
171  mInputParentGeantId = 0 ;
172 }
173 
174 //____________________________________________________________________________________________________
175 StEmbeddingQADraw::StEmbeddingQADraw(const TString embeddingFile, const TString realDataFile,
176  const Int_t year, const TString production, const Int_t geantid, const Bool_t isEmbeddingOnly)
177  : mIsEmbeddingOnly(isEmbeddingOnly), mIsGIFOn(kFALSE), mIsJPGOn(kFALSE), mIsEPSOn(kFALSE), mIsPSOn(kFALSE), mGeantId(geantid)
178 {
181 
183  mIsOpen = open(embeddingFile, realDataFile);
184 
186  mYear = year ;
187  mProduction = production ;
188 
189  mInputParentGeantId = 0 ;
190 }
191 
192 //____________________________________________________________________________________________________
193 StEmbeddingQADraw::~StEmbeddingQADraw()
194 {
196 
198  mMcGeantId.clear();
199  mDaughterGeantId.clear();
200  mParentGeantId.clear();
201  mParentParentGeantId.clear();
202 
204  if( mInputEmbedding || mInputEmbedding->IsOpen() ) mInputEmbedding->Close();
205  if( mInputRealData || mInputRealData->IsOpen() ) mInputRealData->Close();
206 }
207 
208 //____________________________________________________________________________________________________
209 void StEmbeddingQADraw::setParentGeantId(const Int_t parentgeantid)
210 {
211  mInputParentGeantId = parentgeantid ;
212  LOG_INFO << "StEmbeddingQADraw::setParentGeantId Set parent geantid = "
213  << mInputParentGeantId
214  << endm;
215 }
216 
217 //____________________________________________________________________________________________________
218 void StEmbeddingQADraw::setPtMax(const Double_t ptmax)
219 {
220  mPtMax = ptmax ;
221  LOG_INFO << "StEmbeddingQADraw::setPtMax() Set maximum pt = " << mPtMax
222  << " GeV/c to be drawn" << endm;
223 }
224 
225 
226 //____________________________________________________________________________________________________
227 Bool_t StEmbeddingQADraw::open(const TString embeddingFile, const TString realDataFile)
228 {
230 
231  // OPEN embedding file
232  mInputEmbedding = TFile::Open(embeddingFile);
233 
234  if(!mInputEmbedding || !mInputEmbedding->IsOpen() || mInputEmbedding->IsZombie() ){
235  Error("StEmbeddingQADraw", "can't find %s", embeddingFile.Data());
236  return kFALSE;
237  }
238  LOG_INFO << "OPEN input (embedding) : " << mInputEmbedding->GetName() << endm;
239  LOG_INFO << endm;
240 
242  if( mIsEmbeddingOnly ) return kTRUE ;
243 
244  // OPEN real data file
245  mInputRealData = TFile::Open(realDataFile);
246 
247  if(!mInputRealData || !mInputRealData->IsOpen() || mInputRealData->IsZombie() ){
248  Error("StEmbeddingQADraw", "can't find %s", realDataFile.Data());
249  return kFALSE;
250  }
251  LOG_INFO << "OPEN input (real data) : " << mInputRealData->GetName() << endm;
252  LOG_INFO << endm;
253 
254  return kTRUE ;
255 }
256 
257 
258 //____________________________________________________________________________________________________
259 void StEmbeddingQADraw::init()
260 {
262 
264  mPtMax = 5.0;
265 
267  mMcGeantId.clear();
268  mDaughterGeantId.clear();
269  mParentGeantId.clear();
270  mParentParentGeantId.clear();
271 
272  LOG_INFO << "#------------------------------------------------------------" << endm;
273  LOG_INFO << Form(" Year = %10d", mYear) << endm;
274  LOG_INFO << Form(" Production = %10s", mProduction.Data()) << endm;
275  LOG_INFO << Form(" Particle = %10s", getParticleName()) << endm;
276  LOG_INFO << "#------------------------------------------------------------" << endm;
277 
279  StEmbeddingQAUtilities::instance()->setStyle();
280 
282  checkInputGeantId() ;
283 
285  setDaughterGeantId();
286 
288  LOG_INFO << "Initialize canvas" << endm ;
289  mCanvas = new TCanvas(Form("c%d", mCanvasId), Form("c%d", mCanvasId++), 700, 800);
290  mCanvas->SetFillColor(1);
291  mCanvas->cd();
292 
294  mMainPad = new TPad("mainPad", "", 0.00, 0.00, 1.00, 0.90);
295  mMainPad->SetFillColor(10);
296  mMainPad->Draw();
297 
299  mPDF = new TPDF(Form("%sqa_%s.pdf", mOutputFigureDirectory.Data(), getBaseName()));
300  LOG_INFO << "Initialize PDF file : " << mPDF->GetName() << endm ;
301 
303  mCanvas->cd();
304  TPaveText* info = new TPaveText(0.00, 0.95, 0.63, 1.00);
305  info->SetBorderSize(1);
306  info->SetFillColor(kBlack);
307  info->SetTextFont(42);
308  info->SetTextColor(10);
309  info->SetTextSize(0.030);
310 
311  const Int_t nevents = getEntries() ;
312  TString neventsName(Form("%d", getEntries())) ;
313  if( nevents > 1000000 ) neventsName = Form("%1.2f M", (Double_t)nevents/1.0e+06);
314  else if ( nevents > 1000 ) neventsName = Form("%1.2f K", (Double_t)nevents/1.0e+03);
315 
316  info->AddText(Form("N_{events}=%s, MC=%s, %s, (%d)", neventsName.Data(), getParticleName(), mProduction.Data(), mYear));
317  info->Draw();
318 
320  mMainPad->cd();
321 
322  TLatex* title = new TLatex(0.1, 0.6, "Embedding QA");
323  title->SetTextSize(0.10);
324  title->Draw();
325 
326  TString qaName("embedding + real");
327  if ( mIsEmbeddingOnly ){
328  qaName = "embedding only";
329  }
330 
331  TLatex* qaTitle = new TLatex(0.12, 0.5, qaName);
332  qaTitle->SetTextSize(0.07);
333  qaTitle->Draw();
334 
337  TLine* lineEmbedding = new TLine(0, 0, 1, 1);
338  TLine* lineRealBlue = new TLine(0, 0, 1, 1);
339  TLine* lineRealBlack = new TLine(0, 0, 1, 1);
340  lineEmbedding->SetLineColor(kRed); lineEmbedding->SetLineWidth(2);
341  lineRealBlue ->SetLineColor(kBlue); lineRealBlue ->SetLineWidth(2);
342  lineRealBlack->SetLineColor(kBlack); lineRealBlack->SetLineWidth(2);
343 
344  TLegend* colorCode = new TLegend(0.05, 0.15, 0.8, 0.45);
345  colorCode->SetTextSize(0.035);
346  colorCode->SetFillColor(10);
347  colorCode->SetHeader("Color code for embedding and real data");
348  colorCode->AddEntry(lineRealBlack, "MC (black)", "L");
349  colorCode->AddEntry(lineEmbedding, "Reconstructed embedding tracks* (red)", "L");
350  colorCode->AddEntry(lineRealBlue, "Real** (blue)", "L");
351  colorCode->Draw();
352 
353  TLatex* note0 = new TLatex(0.1, 0.10, "* matched pairs or contaminated pairs");
354  TLatex* note1 = new TLatex(0.1, 0.05, "** black is also used, see legend for each pad");
355  note0->SetTextSize(0.03);
356  note1->SetTextSize(0.03);
357  note0->Draw();
358  note1->Draw();
359 
360  mCanvas->cd();
361  mCanvas->Update();
362 
363  mPDF->NewPage() ;
364 
365  // Print event and track selections
366  const StEmbeddingQAUtilities* utility = StEmbeddingQAUtilities::instance() ;
367  TPaveText* header = initCanvas("Event & track selections");
368 
369  // Event selections
370  TPaveText* eventSelections = new TPaveText(0.05, 0.60, 0.95, 0.92);
371  eventSelections->SetTextAlign(12);
372  eventSelections->SetBorderSize(1);
373  eventSelections->SetFillColor(10);
374  eventSelections->SetTextColor(kBlack);
375  eventSelections->SetTextSize(0.030);
376  eventSelections->SetTextFont(42);
377 
378  eventSelections->AddText("*** Event selections");
379  // z-vertex cut
380  eventSelections->AddText(Form(" z-vertex cut : |v_{z}| < %1.1f cm", utility->getZVertexCut()));
381  // trigger id's
382  const vector<UInt_t> triggerId(utility->getTriggerIdCut());
383  if ( !triggerId.empty() ) {
384  TString triggers("");
385  for(UInt_t i=0; i<triggerId.size(); i++) {
386  triggers += Form("%d", triggerId[i]);
387  if( i != triggerId.size() - 1 ) triggers += ", ";
388  }
389  eventSelections->AddText(Form(" trigger id cut : id = %s", triggers.Data()));
390  }
391  eventSelections->AddText("NOTE: Trigger id cut for real data has to be made manually in doEmbeddingQAMaker.C");
392  eventSelections->SetAllWith("***", "color", kRed);
393  eventSelections->SetAllWith("***", "font", 72);
394  eventSelections->SetAllWith("***", "size", 0.033);
395  eventSelections->SetAllWith("NOTE", "color", kBlue);
396  eventSelections->SetAllWith("NOTE", "font", 72);
397  eventSelections->SetAllWith("NOTE", "size", 0.020);
398  eventSelections->Draw();
399 
400  // Track selections
401  TPaveText* trackSelections = new TPaveText(0.05, 0.05, 0.95, 0.59);
402  trackSelections->SetTextAlign(12);
403  trackSelections->SetBorderSize(1);
404  trackSelections->SetFillColor(10);
405  trackSelections->SetTextColor(kBlack);
406  trackSelections->SetTextSize(0.030);
407  trackSelections->SetTextFont(42);
408  trackSelections->AddText("*** Track selections");
409  trackSelections->AddText(Form(" %1.1f < p_{T} < %1.1f GeV/c", utility->getPtMinCut(), utility->getPtMaxCut())); // pt cut
410  trackSelections->AddText(Form(" |#eta| < %1.2f", utility->getEtaCut()));
411  trackSelections->AddText(Form(" |y| < %1.2f", utility->getRapidityCut()));
412  trackSelections->AddText(Form(" nHitsFit > %3d", utility->getNHitCut()));
413  trackSelections->AddText(Form(" nHitsFit/nHitsPoss > %1.2f", utility->getNHitToNPossCut()));
414  trackSelections->AddText(Form(" global Dca < %1.1f cm", utility->getDcaCut()));
415  trackSelections->AddText(Form(" |n#sigma| < %1.1f, using %s", utility->getNSigmaCut(),utility->getBTofPid()?"BTOF":"TPC dE/dx"));
416  trackSelections->AddText("NOTE1: Rapidity cut for real data has to be made manually in doEmbeddingQAMaker.C");
417  trackSelections->AddText("NOTE2: Cut on its own variable is currently disabled, e.x. no dca cut for dca histogram");
418  trackSelections->SetAllWith("***", "color", kRed);
419  trackSelections->SetAllWith("***", "font", 72);
420  trackSelections->SetAllWith("***", "size", 0.033);
421  trackSelections->SetAllWith("NOTE", "color", kBlue);
422  trackSelections->SetAllWith("NOTE", "font", 72);
423  trackSelections->SetAllWith("NOTE", "size", 0.020);
424  trackSelections->Draw();
425 
426  mCanvas->cd();
427  mCanvas->Update();
428  mPDF->NewPage() ;
429 
430  delete header ;
431 }
432 
433 //____________________________________________________________________________________________________
434 void StEmbeddingQADraw::setOutputDirectory(const TString name)
435 {
437 
438  mOutputFigureDirectory = name ;
439 
441  if( gSystem->AccessPathName(name) == kTRUE ){ // 0 is true, i.e. directory exists
442  Error("setOutputDirectory", "Directory %s does not exist. Set current directory as the output location",name.Data());
443  mOutputFigureDirectory = "./";
444  }
445 
447  if ( mOutputFigureDirectory.Last('/') + 1 != mOutputFigureDirectory.Length() ){
448  LOG_INFO << endm;
449  LOG_INFO << "Put / at the end of output directory name" << endm;
450  mOutputFigureDirectory.Append("/");
451  }
452 
453  LOG_INFO << "Set output directory for figures : " << mOutputFigureDirectory << endm;
454 }
455 
456 //____________________________________________________________________________________________________
457 void StEmbeddingQADraw::print(const TCanvas& canvas, const TString name) const
458 {
460 
461  if( mIsPNGOn ) canvas.Print(Form("%s%s_%s.png", mOutputFigureDirectory.Data(), name.Data(), getBaseName()));
462  if( mIsGIFOn ) canvas.Print(Form("%s%s_%s.gif", mOutputFigureDirectory.Data(), name.Data(), getBaseName()));
463  if( mIsJPGOn ) canvas.Print(Form("%s%s_%s.jpg", mOutputFigureDirectory.Data(), name.Data(), getBaseName()));
464  if( mIsEPSOn ) canvas.Print(Form("%s%s_%s.eps", mOutputFigureDirectory.Data(), name.Data(), getBaseName()));
465  if( mIsPSOn ) canvas.Print(Form("%s%s_%s.ps", mOutputFigureDirectory.Data(), name.Data(), getBaseName()));
466 }
467 
468 //____________________________________________________________________________________________________
469 void StEmbeddingQADraw::checkInputGeantId()
470 {
473 
474  TH1* hGeantId = (TH1D*) getHistogram("hGeantId_0");
475 
477  if(!hGeantId){
478  mGeantId = -1 ;
479  return ;
480  }
481 
483  for(Int_t id=0; id<hGeantId->GetNbinsX(); id++){
484  const Bool_t isGeantIdOk = hGeantId->GetBinContent(id+1)>0.0;
485  if(!isGeantIdOk) continue ;
486 
487  const Int_t geantid = TMath::Nint(hGeantId->GetBinLowEdge(id+1));
488  mMcGeantId.push_back(geantid);
489  }
490 
491  const StEmbeddingQAUtilities* utility = StEmbeddingQAUtilities::instance() ;
492 
493  if( mMcGeantId.size() == 1 ){
498  if ( mGeantId == mMcGeantId[0] ) return ;
499  else{
500  LOG_INFO << "#----------------------------------------------------------------------------------------------------" << endm;
501  LOG_INFO << " Input geant id doesn't correspond to the geant id in the MC histogram" << endm ;
502  LOG_INFO << " Geantid in the MC histogram, geantid = " << mMcGeantId[0]
503  << ", particle name = " << utility->getParticleDefinition(mMcGeantId[0])->name().c_str() << endm;
504  LOG_INFO << " Do you want to proceed to the QA for geantid = " << mMcGeantId[0] << " ?" << endm;
505  LOG_INFO << " [yes=1/no=0]:" << endm;
506  UInt_t proceedQA = 0 ;
507  cin >> proceedQA ;
508  LOG_INFO << "#----------------------------------------------------------------------------------------------------" << endm;
509 
511  assert(proceedQA);
512 
514  mGeantId = mMcGeantId[0] ;
515  }
516  }
517  else{
522  for(vector<Int_t>::iterator iter = mMcGeantId.begin(); iter != mMcGeantId.end(); iter++){
523  const Int_t geantidFound = (*iter) ;
524  if ( mGeantId == geantidFound ) return ;
525  }
526 
528  LOG_INFO << "#----------------------------------------------------------------------------------------------------" << endm;
529  LOG_INFO << " Cannot find input geantid in the MC histogram" << endm;
530  LOG_INFO << " Available geantid's are listed below;" << endm;
531  LOG_INFO << " geantid particle name" << endm;
532  for(vector<Int_t>::iterator iter = mMcGeantId.begin(); iter != mMcGeantId.end(); iter++){
533  const Int_t geantidFound = (*iter) ;
534  LOG_INFO << Form(" %5d %10s", geantidFound,
535  utility->getParticleDefinition(geantidFound)->name().c_str()) << endm;
536  }
537  LOG_INFO << " Which geantid you want to do the QA ?" << endm;
538  LOG_INFO << " [Put the one geantid listed above. If you want to stop the program, please put 0]:" << endm;
539  cin >> mGeantId ;
540 
542  assert(mGeantId);
543 
544  LOG_INFO << " QA for geantid = " << mGeantId
545  << ", particle name = " << utility->getParticleDefinition(mGeantId)->name().c_str() << endm;
546  LOG_INFO << "#----------------------------------------------------------------------------------------------------" << endm;
547  }
548 
549 }
550 
551 //____________________________________________________________________________________________________
552 Bool_t StEmbeddingQADraw::isOpen() const
553 {
559 
560  if(!mIsOpen) LOG_INFO << "No input files found" << endm;
561  assert(mIsOpen);
562 
563  const Bool_t isGeantIdOk = mGeantId>0 ;
564  if(!isGeantIdOk) LOG_INFO << "Cannot find input geantid in the histogram or no histogram in the input" << endm;
565  assert(isGeantIdOk);
566 
567  return kTRUE ;
568 }
569 
570 //____________________________________________________________________________________________________
571 Bool_t StEmbeddingQADraw::isMatchedPairOk() const
572 {
574  // NOTE: If there are no matched pairs in the minimc due to the error
575  // this function won't work.
576  // The condition "!isDecay()" would give us the right answer
577 // TH1* hGeantId = (TH1D*) getHistogram("hGeantId_1");
578 // if(!hGeantId) return kFALSE ;
579 // return (hGeantId->GetEntries()>0) ;
580 
582  // NOTE: Except for electrons, pions, kaons and protons
583  // some of those are assigned as "unstable"
584  if ( StEmbeddingQAUtilities::instance()->isEPiKP(mGeantId) ) {
585  // stable
586  return kTRUE ;
587  }
588  else if ( StEmbeddingQAUtilities::instance()->isGamma(mGeantId) ) {
589  // unstable for gamma
590  return kFALSE ;
591  }
592  else{
593  return StEmbeddingQAUtilities::instance()->getParticleDefinition(mGeantId)->stable() ;
594  }
595 }
596 
597 
598 //____________________________________________________________________________________________________
599 void StEmbeddingQADraw::setDaughterGeantId()
600 {
602 
604  if(isMatchedPairOk()){
606  mDaughterGeantId.push_back(mGeantId);
607  // Set 0 for parent and parent-parent geantid
608  mParentGeantId.push_back(0);
609  mParentParentGeantId.push_back(0);
610  return;
611  }
612 
614  for(UInt_t id=0; id<1000; id++){
615  const StEmbeddingQAUtilities* utility = StEmbeddingQAUtilities::instance() ;
616  const Int_t contamCategoryId = utility->getCategoryId("CONTAM") ;
617 
618  Int_t n = 0 ;
619  if( mInputParentGeantId == 0 ) n = 1 ;
620  else n = 2 ;
621 
622  for(Int_t i=0; i<n; i++) {
623  Int_t daughterId = id ;
624  Int_t parentId = mGeantId ;
625  Int_t parentParentId = 0 ;
626  if( i == 1 ) {
627  daughterId = id ;
628  parentId = mInputParentGeantId ;
629  parentParentId = mGeantId ;
630  }
631 // hDca = (TH3D*) mInputEmbedding->Get(Form("hDca_%d_%d_%d_%d", contamCategoryId, 0, mGeantId, id));
632 // hDca = (TH3D*) mInputEmbedding->Get(Form("hDca_%d_%d_%d_%d", contamCategoryId, mGeantId, mInputParentGeantId, id));
633  TH3* hDca = (TH3D*) mInputEmbedding->Get(Form("hDca_%d_%d_%d_%d", contamCategoryId, parentParentId, parentId, daughterId)) ;
634 
635  if ( hDca ){
636  if( i == 0 ){
637  const Char_t* daughterName = utility->getParticleDefinition(daughterId)->name().c_str() ;
638  const Char_t* parentName = utility->getParticleDefinition(parentId)->name().c_str() ;
639 
640  LOG_INFO << Form("Find daughter %10s from parent %10s", daughterName, parentName) << endm;
641  }
642  else{
643  const Char_t* daughterName = utility->getParticleDefinition(daughterId)->name().c_str() ;
644  const Char_t* parentName = utility->getParticleDefinition(parentId)->name().c_str() ;
645  const Char_t* parentparentName = utility->getParticleDefinition(parentParentId)->name().c_str() ;
646 
647  LOG_INFO << Form("Find daughter %10s from parent %10s (from %10s)", daughterName, parentName, parentparentName) << endm;
648  }
649 
650  mDaughterGeantId.push_back(id);
651  mParentGeantId.push_back(parentId);
652  mParentParentGeantId.push_back(parentParentId);
653  }
654  }
655  }
656 
658  if ( mDaughterGeantId.empty() ) mDaughterGeantId.push_back(mGeantId);
659 }
660 
661 //____________________________________________________________________________________________________
662 Int_t StEmbeddingQADraw::getEntries() const
663 {
665 
666  TH1* hVz = (TH1D*) getHistogram("hVzAccepted");
667  if(!hVz) return 0;
668 
669  return (Int_t)( hVz->GetEntries() );
670 }
671 
672 //____________________________________________________________________________________________________
673 Bool_t StEmbeddingQADraw::isDecay() const
674 {
676 
679  if(isMatchedPairOk()) return kFALSE ;
680 
682  return ( mDaughterGeantId.size() > 1 ) ;
683 }
684 
685 //____________________________________________________________________________________________________
686 Int_t StEmbeddingQADraw::getGeantIdReal(const Int_t daughter) const
687 {
690 
691  if(mDaughterGeantId[daughter]==6) return 9; // if mu-, use pi-
692  return ( StEmbeddingQAUtilities::instance()->isEPiKP(mDaughterGeantId[daughter]) ) ? mDaughterGeantId[daughter] : 8 ;
693 }
694 
695 //____________________________________________________________________________________________________
696 Int_t StEmbeddingQADraw::getCategoryId(const Bool_t isEmbedding) const
697 {
703  const StEmbeddingQAUtilities* utility = StEmbeddingQAUtilities::instance() ;
704 
705  if ( isEmbedding ){
706  return (isDecay()) ? utility->getCategoryId("CONTAM") : utility->getCategoryId("MATCHED") ;
707  }
708  else{
709  return utility->getCategoryId("PRIMARY") ;
710  }
711 }
712 
713 //____________________________________________________________________________________________________
714 TObject* StEmbeddingQADraw::getHistogram(const TString name, const Bool_t isEmbedding) const
715 {
717 
718  if(!isOpen()) return 0;
719 
722  if( mIsEmbeddingOnly && !isEmbedding) return 0 ;
723 
724  TObject* obj = (isEmbedding) ? mInputEmbedding->Get(name) : mInputRealData->Get(name) ;
725  if ( !obj ){
726  Error("StEmbeddingQADraw::getHistogram", "Histogram %s doesn't exist", name.Data());
727  return 0;
728  }
729 
730  LOG_DEBUG << "StEmbeddingQADraw::getHistogram() get histogram = " << name << endm;
731 
732  return obj ;
733 }
734 
735 
736 //____________________________________________________________________________________________________
737 TObject* StEmbeddingQADraw::getHistogram(const TString name, const UInt_t daughter, const Bool_t isEmbedding,
738  const UInt_t parentparentid) const
739 {
744 
745  if ( daughter >= mDaughterGeantId.size() ){
746  Error("StEmbeddingQADraw::getHistogram", "Unknown daughter index, index=%3d", daughter);
747  return 0;
748  }
749 
750  return getHistogram(getHistogramName(name, daughter, isEmbedding, parentparentid), isEmbedding) ;
751 }
752 
753 //____________________________________________________________________________________________________
754 const Char_t* StEmbeddingQADraw::getHistogramName(const TString name, const UInt_t daughter, const Bool_t isEmbedding,
755  const UInt_t parentparentid) const
756 {
757  const Int_t category = getCategoryId(isEmbedding) ;
758 
759  if( isEmbedding ){
763 
764  if( isDecay() ) {
765 // if ( mParentGeantId == 0 ) return Form("%s_%d_%d_%d_%d", name.Data(), category, parentparentid, mGeantId, mDaughterGeantId[daughter]) ;
766 // else return Form("%s_%d_%d_%d_%d", name.Data(), category, mGeantId, mParentGeantId, mDaughterGeantId[daughter]) ;
767  return Form("%s_%d_%d_%d_%d", name.Data(), category, mParentParentGeantId[daughter], mParentGeantId[daughter], mDaughterGeantId[daughter]) ;
768  }
769  else {
770  return Form("%s_%d_%d", name.Data(), category, mGeantId) ;
771  }
772  }
773  else{
777  const Int_t geantidReal = getGeantIdReal(daughter) ;
778 
779  return Form("%s_%d_%d", name.Data(), category, geantidReal) ;
780  }
781 }
782 
783 //____________________________________________________________________________________________________
784 Double_t StEmbeddingQADraw::getNormalization(const TH1& h) const
785 {
788 
789  const Double_t ntrack = h.Integral() ;
790  if( ntrack == 0 ) return 1.0 ;
791 
792  const Double_t binWidth = h.GetBinWidth(1) ;
793  if( binWidth == 0.0 ) return 1.0 ;
794 
795  return 1.0/(ntrack*binWidth);
796 }
797 
798 //____________________________________________________________________________________________________
799 const Char_t* StEmbeddingQADraw::getBaseName() const
800 {
802 
803  return Form("%d_%s_%s", mYear, mProduction.Data(), getParticleName());
804 }
805 
806 //____________________________________________________________________________________________________
807 Bool_t StEmbeddingQADraw::drawStatistics(const Double_t x1, const Double_t y1, const Double_t x2, const Double_t y2,
808  const Double_t textSize) const
809 {
815 
816  TPaveText* statistics = new TPaveText(x1, y1, x2, y2);
817  statistics->SetTextFont(42);
818  statistics->SetTextSize(textSize);
819  statistics->SetBorderSize(1);
820  statistics->SetFillColor(10);
821  statistics->AddText(Form("N_{evts} = %d", getEntries()));
822  statistics->AddText(Form("Year: %d", mYear));
823  statistics->AddText(Form("Production: %s", mProduction.Data()));
824  statistics->AddText(Form("Particle: %s", getParticleName()));
825  statistics->Draw();
826 
827  return kTRUE ;
828 }
829 
830 //____________________________________________________________________________________________________
831 TPaveText* StEmbeddingQADraw::drawHeader(const TString description,
832  const Double_t x1, const Double_t y1, const Double_t x2, const Double_t y2, const Double_t textSize) const
833 {
835  TPaveText* header = new TPaveText(x1, y1, x2, y2);
836  header->SetBorderSize(1);
837  header->SetFillColor(kBlack);
838  header->SetTextColor(10);
839  header->SetTextSize(textSize);
840  header->SetTextFont(42);
841  header->AddText(description);
842  header->Draw();
843 
844  return header ; // need to be deleted after printing PDF
845 }
846 
847 //____________________________________________________________________________________________________
848 void StEmbeddingQADraw::drawLegend(const UInt_t id, TH1* hembed, TH1* hreal,
849  const Option_t* option, const Bool_t doSplit) const
850 {
851  TLegend* leg = new TLegend(0, 0.2, 1, 0.8);
852  leg->SetFillColor(10);
853 
854  // Use fixed font size
855  leg->SetTextFont(43);
856  leg->SetTextSize(15);
857 
858  leg->AddEntry(hembed, getEmbeddingParticleName(id, doSplit), option) ;
859  if(hreal) leg->AddEntry(hreal, getRealParticleName(id, doSplit), option) ;
860  leg->Draw();
861 }
862 
863 //____________________________________________________________________________________________________
864 void StEmbeddingQADraw::drawErrorMessages(const TString histogramName) const
865 {
866  const StEmbeddingQAUtilities* utility = StEmbeddingQAUtilities::instance() ;
867  const Int_t categoryId = getCategoryId() ;
868  const TString title(utility->getCategoryTitle(categoryId));
869 
870  TPaveText* error0 = new TPaveText(0.05, 0.60, 0.95, 0.92, "NDC");
871  error0->SetTextAlign(12);
872  error0->SetBorderSize(1);
873  error0->SetFillColor(10);
874  error0->SetTextColor(kBlack);
875  error0->SetTextSize(0.030);
876  error0->SetTextFont(42);
877  error0->AddText(Form("*** No histogram \"%s\" found", histogramName.Data()));
878  error0->AddText(Form("Please make sure that you have %s in the minimc.root.", title.Data()));
879  error0->AddText("Open minimc.root file and then check number of tracks.");
880  error0->AddText(Form("There are maybe no %s in the minimc.", title.Data()));
881  error0->AddText("Or you used older base QA codes.");
882  error0->AddText("In case you have finite number of tracks,");
883  error0->AddText("please also have a look at geantid.");
884  error0->SetAllWith("***", "color", kRed);
885  error0->SetAllWith("***", "font", 72);
886  error0->SetAllWith("***", "size", 0.033);
887 
888 
889  TPaveText* error1 = new TPaveText(0.05, 0.05, 0.95, 0.55, "NDC");
890  error1->SetTextAlign(12);
891  error1->SetBorderSize(1);
892  error1->SetFillColor(10);
893  error1->SetTextColor(kBlack);
894  error1->SetTextSize(0.025);
895  error1->SetTextFont(42);
896  error1->AddText("See example below how to check number of tracks and geantid");
897 
898  if ( isDecay() ) {
899  error1->AddText("[ROOT]> StMiniMcTree->Draw(\"mNContamPair\")");
900  error1->AddText("or");
901  error1->AddText("[ROOT]> StMiniMcTree->Scan(\"mNContamPair\")");
902  error1->AddText(Form("For geantid in %s", title.Data()));
903  error1->AddText("[ROOT]> StMiniMcTree->Draw(\"mContamPairs.mGeantId\")");
904  error1->AddText("or use \"Scan()\" function");
905  }
906  else {
907  error1->AddText("[ROOT]> StMiniMcTree->Draw(\"mNMatchedPair\")");
908  error1->AddText("or");
909  error1->AddText("[ROOT]> StMiniMcTree->Scan(\"mNMatchedPair\")");
910  error1->AddText(Form("For geantid in %s", title.Data()));
911  error1->AddText("[ROOT]> StMiniMcTree->Draw(\"mMatchedPairs.mGeantId\")");
912  error1->AddText("or use \"Scan()\" function");
913  }
914  error1->SetAllWith("ROOT", "color", kBlue);
915  error1->SetAllWith("ROOT", "font", 52);
916  error1->SetAllWith("ROOT", "size", 0.028);
917 
918  error0->Draw();
919  error1->Draw();
920 }
921 
922 //____________________________________________________________________________________________________
923 const Char_t* StEmbeddingQADraw::getParticleName(const Int_t geantid) const
924 {
926 
928  const Int_t id = (geantid<0) ? mGeantId : geantid ;
929  const StParticleDefinition* particle = StEmbeddingQAUtilities::instance()->getParticleDefinition(id) ;
930 
931  if(!particle){
932  Error("StEmbeddingQADraw::getParticleName", "Can't find particle geantid=%d in StParticleDefinition", id);
933  assert(0);
934  }
935 
941  static TString name;
942  name = particle->name().c_str();
943  while( name.Contains("/") ) name.Remove(name.Index("/"), 1); // Remove "/" from particle name
944  while( name.Contains("(") ) name.Remove(name.Index("("), 1); // Remove "(" from particle name
945  while( name.Contains(")") ) name.Remove(name.Index(")"), 1); // Remove "/" from particle name
946  while( name.Contains("*") ) name.Replace(name.Index("*"), 1, "star"); // Replace "*" to "star"
947 
948  return name.Data() ;
949 }
950 
951 //____________________________________________________________________________________________________
952 const Char_t* StEmbeddingQADraw::getMcParticleName() const
953 {
954  const StEmbeddingQAUtilities* utility = StEmbeddingQAUtilities::instance() ;
955  const Int_t categoryId = getCategoryId() ;
956  const TString name(utility->getCategoryName(categoryId));
957  const TString parent( (utility->isContaminated(name)) ? "Parent " : "" );
958 
959  return Form("%s%s (MC, geantid=%d)", parent.Data(), getParticleName(mGeantId), mGeantId);
960 }
961 
962 //____________________________________________________________________________________________________
963 const Char_t* StEmbeddingQADraw::getEmbeddingParticleName(const UInt_t id, const Bool_t doSplit) const
964 {
965  if(id>=mDaughterGeantId.size()){
966  Error("StEmbeddingQADraw::getEmbeddingParticleName", "Unknown daughter particle index, id=%3d", id);
967  return "";
968  }
969 
970  const StEmbeddingQAUtilities* utility = StEmbeddingQAUtilities::instance() ;
971  const Int_t categoryId = getCategoryId() ;
972  const TString name(utility->getCategoryName(categoryId));
973 
974 // const TString parent( (mParentGeantId==0) ? "" : Form(" (from %s)", getParticleName(mParentGeantId)) );
975 // const TString daughter( (utility->isContaminated(name)) ? "Daughter " : "" );
976 // const Int_t daughterId = mDaughterGeantId[id] ;
977  const Int_t parentParentId = mParentParentGeantId[id] ;
978  const Int_t parentId = mParentGeantId[id] ;
979  const Int_t daughterId = mDaughterGeantId[id] ;
980 
981  const TString daughter( (utility->isContaminated(name)) ? "Daughter " : "" );
982  TString parent("");
983  // For unstable particles
984  if ( parentId != 0 ) {
985  const TString parentName(getParticleName(parentId));
986 
987  if ( parentParentId == 0 ) {
988  parent = " (from " + parentName + ")";
989  }
990  else{
991  const TString parentParentName(getParticleName(parentParentId));
992  parent = " (from " + parentName + " #leftarrow " + parentParentName + ")";
993  }
994  }
995 
996  return (doSplit) ? Form("#splitline{%s%s%s}{(%s, geantid=%d)}", daughter.Data(), getParticleName(daughterId),
997  parent.Data(), utility->getCategoryName(categoryId).Data(), daughterId)
998  : Form("%s%s%s (%s, geantid=%d)", daughter.Data(), getParticleName(daughterId),
999  parent.Data(), utility->getCategoryName(categoryId).Data(), daughterId);
1000 }
1001 
1002 //____________________________________________________________________________________________________
1003 const Char_t* StEmbeddingQADraw::getRealParticleName(const UInt_t id, const Bool_t doSplit) const
1004 {
1005  // nsigma < 2 is hard-coded (2 sigma cut is implemented in StEmbeddingQATrack.cxx).
1006  const StEmbeddingQAUtilities* utility = StEmbeddingQAUtilities::instance() ;
1007 
1008  if(id>=mDaughterGeantId.size()){
1009  Error("StEmbeddingQADraw::getRealParticleName", "Unknown daughter particle index, id=%3d", id);
1010  return "";
1011  }
1012 
1013  const Int_t geantid = getGeantIdReal(id) ;
1014  return (doSplit) ? Form("#splitline{%s}{(%s, |n #sigma_{%s}|<2 %s)}", getParticleName(geantid),
1015  StEmbeddingQAUtilities::instance()->getCategoryName(getCategoryId(kFALSE)).Data(), getParticleName(geantid), utility->getBTofPid()?"BTOF":"TPC")
1016  : Form("%s (%s, |n #sigma_{%s}|<2 %s)", getParticleName(geantid),
1017  StEmbeddingQAUtilities::instance()->getCategoryName(getCategoryId(kFALSE)).Data(), getParticleName(geantid), utility->getBTofPid()?"BTOF":"TPC");
1018 }
1019 
1020 //____________________________________________________________________________________________________
1021 Bool_t StEmbeddingQADraw::drawEvent() const
1022 {
1024 
1025  LOG_INFO << endm;
1026  LOG_INFO << "----------------------------------------------------------------------------------------------------" << endm ;
1027  LOG_INFO << "QA for Event-wise informations ..." << endm;
1028 
1030  if(!isOpen()) return kFALSE ;
1031 
1032  //----------------------------------------------------------------------------------------------------
1033  // Event-wise informations
1034  //----------------------------------------------------------------------------------------------------
1035 
1036  //----------------------------------------------------------------------------------------------------
1037  // Vertices
1038  drawVertices() ;
1039 
1040  //----------------------------------------------------------------------------------------------------
1042  drawRunEventId() ;
1043 
1044  //----------------------------------------------------------------------------------------------------
1046  drawNumberOfParticles() ;
1047 
1048  LOG_INFO << "----------------------------------------------------------------------------------------------------" << endm ;
1049  LOG_INFO << endm;
1050 
1051  return kTRUE ;
1052 }
1053 
1054 //____________________________________________________________________________________________________
1055 Bool_t StEmbeddingQADraw::drawVertices() const
1056 {
1058  LOG_INFO << "QA for event-vertices ..." << endm;
1059 
1061  const Double_t vzMin = getVzAcceptedMinimum() ;
1062  const Double_t vzMax = getVzAcceptedMaximum() ;
1063  TPaveText* header = initCanvas(Form("Event vertices, offline cuts: %1.1f < v_{z} < %1.1f cm", vzMin, vzMax), 2, 2);
1064 
1066  mMainPad->cd(1);
1067  TH1* hVz = (TH1D*) getHistogram("hVz");
1068  if(!hVz) return kFALSE ;
1069 
1070  TH1* hVzAccepted = (TH1D*) getHistogram("hVzAccepted");
1071  if(!hVzAccepted) return kFALSE ;
1072 
1073  hVzAccepted->SetLineColor(kCyan);
1074  hVzAccepted->SetFillColor(kCyan);
1075  hVz->Draw();
1076  hVzAccepted->Draw("same");
1077  hVz->Draw("same");
1078 
1080  mMainPad->cd(2);
1081  TH2* hVyVx = (TH2D*) getHistogram("hVyVx");
1082  if(!hVyVx) return kFALSE ;
1083  hVyVx->Draw("colz");
1084 
1086  mMainPad->cd(3);
1087  TH2* hVxVz = (TH2D*) getHistogram("hVxVz");
1088  if(!hVxVz) return kFALSE ;
1089  hVxVz->SetAxisRange(vzMin, vzMax, "X");
1090  hVxVz->Draw("colz");
1091 
1093  mMainPad->cd(4);
1094  TH2* hVyVz = (TH2D*) getHistogram("hVyVz");
1095  if(!hVyVz) return kFALSE ;
1096  hVyVz->SetAxisRange(vzMin, vzMax, "X");
1097  hVyVz->Draw("colz");
1098 
1099  mCanvas->cd();
1100  mCanvas->Update();
1101  mPDF->NewPage() ;
1102 
1104  delete header ;
1105  header = initCanvas("Event vertices, #Deltav = v(Data) - v(MC)", 2, 2);
1106 
1108  for(Int_t i=0; i<3; i++){
1109  mMainPad->cd(i+1);
1110  TH1* hdV = 0 ;
1111  if( i == 0 ) hdV = (TH1D*) getHistogram("hdVx");
1112  if( i == 1 ) hdV = (TH1D*) getHistogram("hdVy");
1113  if( i == 2 ) hdV = (TH1D*) getHistogram("hdVz");
1114  if(!hdV) return kFALSE ;
1115 
1116  hdV->Draw();
1117  }
1118 
1119  mCanvas->cd();
1120  mCanvas->Update();
1121  mPDF->NewPage() ;
1122 
1123  //----------------------------------------------------------------------------------------------------
1124 
1125  delete header ;
1126 
1127  return kTRUE ;
1128 }
1129 
1130 //____________________________________________________________________________________________________
1131 Bool_t StEmbeddingQADraw::drawRunEventId() const
1132 {
1134  LOG_INFO << "QA for run and event id ..." << endm;
1135 
1136  TPaveText* header = initCanvas("Run and event id", 2, 2);
1137 
1138  mMainPad->cd(1);
1139  TH1* hRunNumber = (TH1D*) mInputEmbedding->Get("hRunNumber");
1140  if(!hRunNumber) return kFALSE ;
1141  hRunNumber->Draw();
1142 
1143  mMainPad->cd(2);
1144  TH1* hEventId = (TH1D*) mInputEmbedding->Get("hEventId");
1145  if(!hEventId) return kFALSE ;
1146  hEventId->Draw();
1147 
1149  TPaveText* runlist[2];
1150  for(Int_t i=0;i<2;i++){
1151  runlist[i] = new TPaveText(0.1, 0.05, 0.9, 0.95);
1152  runlist[i]->SetTextFont(42);
1153  runlist[i]->SetTextSize(0.04);
1154  runlist[i]->SetBorderSize(1);
1155  runlist[i]->SetFillColor(10);
1156  runlist[i]->AddText("Run id statistics");
1157  }
1158 
1159  // Count total number of runs
1160  Int_t runTotal = 0 ;
1161  for(Int_t irun=0; irun<hRunNumber->GetNbinsX(); irun++){
1162  if(hRunNumber->GetBinContent(irun+1)>0.0) runTotal++;
1163  }
1164 
1165  // If number of runs >= 10, divided into 2 pads
1166  Int_t runAccepted = 0 ;
1167  for(Int_t irun=0; irun<hRunNumber->GetNbinsX(); irun++){
1168  const Double_t count = hRunNumber->GetBinContent(irun+1);
1169  if(count==0.0) continue ;
1170 
1171  Int_t pad = 0;
1172  if( runTotal < 10 ) pad = 0 ; // always drawing left pad
1173  else{
1174  // Divide into 2 pads
1175  if( runAccepted < runTotal/2 ) pad = 0;
1176  else pad = 1 ;
1177  }
1178 
1179  const Int_t runid = StEmbeddingQAUtilities::instance()->getRunId((Int_t)hRunNumber->GetBinLowEdge(irun+1), mYear);
1180  runlist[pad]->AddText(Form("%10d %10d events", runid, (Int_t)count));
1181 
1182  runAccepted++ ;
1183  }
1184 
1185  mMainPad->cd(3) ;
1186  runlist[0]->Draw();
1187 
1188  mMainPad->cd(4) ;
1189  runlist[1]->Draw();
1190 
1191  mCanvas->cd();
1192  mCanvas->Update();
1193 
1194  mPDF->NewPage() ;
1195  //----------------------------------------------------------------------------------------------------
1196 
1197  delete header ;
1198 
1199  return kTRUE ;
1200 }
1201 
1202 //____________________________________________________________________________________________________
1203 Bool_t StEmbeddingQADraw::drawNumberOfParticles() const
1204 {
1206  LOG_INFO << "QA for multiplicity distributions (number of particles per event) ..." << endm;
1207 
1208  TPaveText* header = initCanvas("Multiplicity distribution", 2, 3);
1209 
1210  for(UInt_t ic=0; ic<StEmbeddingQAConst::mNEmbedding; ic++){
1211  mMainPad->cd(ic+1);
1212  mMainPad->GetPad(ic+1)->SetLogy();
1213 
1214  TH1* hNParticles = (TH1D*) mInputEmbedding->Get(Form("hNParticles_%d", ic));
1215  if(!hNParticles) return kFALSE ;
1216 
1217  TString title(hNParticles->GetTitle());
1218  title.Remove(0, title.Index(',')+2);
1219  hNParticles->SetTitle(title);
1220 
1221  // Set x-axis range
1222  if( hNParticles->GetMean() == 0 && hNParticles->GetRMS() == 0 ){
1223  // If no particles in the branch
1224  hNParticles->SetAxisRange(0, 10, "X");
1225 
1226  LOG_DEBUG << "StEmbeddingQADraw::drawNumberOfParticles no particles for " << title << endm ;
1227  }
1228  else if( hNParticles->GetRMS() == 0 ){
1229  // Fixed multiplicity should have RMS = 0
1230  const Double_t mean = hNParticles->GetMean() ;
1231  const Double_t xmin = (mean-10.0<0.0) ? 0.0 : mean-10.0 ;
1232  const Double_t xmax = mean + 10.0 ;
1233  hNParticles->SetAxisRange(xmin, xmax, "X");
1234 
1235  LOG_DEBUG << "StEmbeddingQADraw::drawNumberOfParticles constant number of particles = " << mean
1236  << " for " << title
1237  << endm;
1238  }
1239  else{
1240  // Varied multiplicity (like 5% of refmult)
1241 
1242  // Try to find the maximum multiplicity
1243  Double_t xmax = 0.0 ;
1244  for(Int_t i=hNParticles->GetNbinsX()-1; i!=0; i--){
1245  if( hNParticles->GetBinContent(i+1) != 0.0 ){
1246  // Find first bin, contains finite bin content
1247  xmax = hNParticles->GetBinCenter(i+1) ;
1248  break ;
1249  }
1250  }
1251  hNParticles->SetAxisRange(0, xmax+10, "X");
1252 
1253  LOG_DEBUG << "StEmbeddingQADraw::drawNumberOfParticles varied multiplicity, max range = " << xmax
1254  << " for " << title
1255  << endm;
1256  }
1257 
1258  hNParticles->Draw();
1259  }
1260 
1261  mCanvas->cd();
1262  mCanvas->Update();
1263 
1264  mPDF->NewPage() ;
1265 
1266  delete header ;
1267 
1268  return kTRUE ;
1269 }
1270 
1271 //____________________________________________________________________________________________________
1272 Bool_t StEmbeddingQADraw::drawMcTrack() const
1273 {
1275 
1276  LOG_INFO << "QA for MC tracks (pt, eta, y, phi) ..." << endm;
1277 
1279  if(!isOpen()) return kFALSE ;
1280 
1281  //----------------------------------------------------------------------------------------------------
1282  // Track-wise informations
1283  TPaveText* header = initCanvas("MC track QA (2D)", 2, 2);
1284 
1286  mMainPad->cd(1);
1287  TH2* hPtVsEta = (TH2D*) getHistogram(Form("hPtVsEta_0_%d", mGeantId));
1288  if(!hPtVsEta) return kFALSE;
1289 
1290  hPtVsEta->SetTitle("");
1291  hPtVsEta->Draw("colz");
1292  hPtVsEta->SetAxisRange(0, mPtMax, "Y");
1293 
1295  mMainPad->cd(2);
1296  TH2* hPtVsY = (TH2D*) getHistogram(Form("hPtVsY_0_%d", mGeantId));
1297  if(!hPtVsY) return kFALSE ;
1298 
1299  hPtVsY->SetTitle("");
1300  hPtVsY->Draw("colz");
1301  hPtVsY->SetAxisRange(0, mPtMax, "Y");
1302 
1304  mMainPad->cd(3);
1305  TH2* hPtVsPhi = (TH2D*) getHistogram(Form("hPtVsPhi_0_%d", mGeantId));
1306  if(!hPtVsPhi) return kFALSE ;
1307 
1308  hPtVsPhi->SetTitle("");
1309  hPtVsPhi->Draw("colz");
1310  hPtVsPhi->SetAxisRange(0, mPtMax, "Y");
1311 
1312  mCanvas->cd();
1313  mCanvas->Update();
1314 
1315  mPDF->NewPage() ;
1316 
1317  delete header ;
1318 
1319  // 1D projections
1320  header = initCanvas("MC track QA (1D)", 2, 2);
1321 
1323  TH1* hPt = (TH1D*) hPtVsPhi->ProjectionY("hPtMc");
1324  TH1* hEta = (TH1D*) hPtVsEta->ProjectionX("hEtaMc");
1325  TH1* hY = (TH1D*) hPtVsY->ProjectionX("hYMc");
1326  TH1* hPhi = (TH1D*) hPtVsPhi->ProjectionX("hPhiMc");
1327  hPt->SetTitleOffset(1.0, "X");
1328  hPt->SetXTitle("p_{T} (GeV/c)");
1329  hEta->SetXTitle("#eta");
1330  hY->SetXTitle("rapidity");
1331 
1332  hPt ->SetTitle("");
1333  hEta->SetTitle("");
1334  hY ->SetTitle("");
1335  hPhi->SetTitle("");
1336 
1337  hPt ->Sumw2() ;
1338  hEta->Sumw2() ;
1339  hY ->Sumw2() ;
1340  hPhi->Sumw2() ;
1341 
1342  hPt->SetMinimum(0.0);
1343  hPt->SetMaximum(hPt->GetMaximum()*1.2);
1344  hPt->SetAxisRange(0, mPtMax, "X");
1345  hEta->SetMinimum(0.0);
1346  hEta->SetMaximum(hEta->GetMaximum()*1.2);
1347  hY->SetMinimum(0.0);
1348  hY->SetMaximum(hY->GetMaximum()*1.2);
1349  hPhi->SetMinimum(0.0);
1350  hPhi->SetMaximum(hPhi->GetMaximum()*1.2);
1351 
1352  // 1. pT
1353  // 3. y
1354  // 3. phi
1355  // 4. eta
1356  gStyle->SetOptFit(1);
1357  for(Int_t ipad=0; ipad<4; ipad++){
1358  mMainPad->cd(ipad+1);
1359  TH1* hdraw = 0 ;
1360  if( ipad == 0 ) hdraw = hPt ;
1361  if( ipad == 1 ) hdraw = hEta ;
1362  if( ipad == 2 ) hdraw = hY ;
1363  if( ipad == 3 ) hdraw = hPhi ;
1364 
1365  hdraw->Draw();
1366 
1367  if( ipad != 0 ){
1368  // Fitting by pol0
1369  Double_t min = -1.0 ;
1370  Double_t max = 1.0 ;
1371  if( ipad == 3 ){ min = -TMath::Pi() ; max = TMath::Pi() ; } // phi
1372 
1373  TF1* pol0Fit = new TF1(Form("pol0Fit_%d", ipad), "pol0", min, max);
1374  pol0Fit->SetLineColor(kRed);
1375  pol0Fit->SetLineWidth(1);
1376  pol0Fit->SetLineStyle(2);
1377  pol0Fit->SetParameter(0, hdraw->GetMean(2));
1378  hdraw->Fit(pol0Fit, "rq0");
1379  pol0Fit->Draw("same");
1380  }
1381  }
1382 
1383  mCanvas->cd();
1384  mCanvas->Update();
1385 
1386  mPDF->NewPage() ;
1387 
1388  delete header ;
1389 
1390  gStyle->SetOptFit(0);
1391 
1392  return kTRUE ;
1393 }
1394 
1395 //____________________________________________________________________________________________________
1396 Bool_t StEmbeddingQADraw::drawTrack() const
1397 {
1399  LOG_INFO << endm ;
1400  LOG_INFO << "#----------------------------------------------------------------------------------------------------" << endm ;
1401  LOG_INFO << "QA for reconstructed tracks ..." << endm ;
1402 
1404  if(!isOpen()) return kFALSE ;
1405 
1406  //----------------------------------------------------------------------------------------------------
1407  // Track-wise informations (Comparison between MC and reconstructed particles)
1409  const Bool_t isGeantIdOk = drawGeantId();
1410 
1412  const Bool_t isPhiOk = drawPhi();
1413 
1415  const Bool_t isRapidityOk = drawRapidity();
1416 
1418  const Bool_t isMomentumOk = drawMomentum();
1419  const Bool_t isPtOk = drawPt();
1420 
1422  const Bool_t isdEdxOk = drawdEdx();
1423 
1425  const Bool_t isDcaOk = drawDca();
1426 
1428  const Bool_t isNHitOk = drawNHit();
1429 
1430  LOG_INFO << "#----------------------------------------------------------------------------------------------------" << endm ;
1431  LOG_INFO << endm ;
1432 
1433  return isGeantIdOk && isPhiOk && isRapidityOk && isMomentumOk && isPtOk && isdEdxOk && isDcaOk && isNHitOk ;
1434 }
1435 
1436 //____________________________________________________________________________________________________
1437 Bool_t StEmbeddingQADraw::drawGeantId() const
1438 {
1440 
1441  LOG_INFO << "QA for geant id ..." << endm;
1442 
1444  if(!isOpen()) return kFALSE ;
1445 
1449  TPaveText* header = initCanvas("Geant id", 1, 2);
1450 
1452  TPad* topPad = (TPad*) mMainPad->cd(1);
1453  topPad->Divide(2, 1);
1454  topPad->cd(1);
1455  TH1* hGeantIdMc = (TH1D*) getHistogram("hGeantId_0");
1456  if(!hGeantIdMc) return kFALSE ;
1457 
1459  const Int_t mcIdBin = hGeantIdMc->GetMaximumBin() ;
1460  const Double_t mcId = hGeantIdMc->GetBinCenter(mcIdBin);
1461  const Double_t mcIdMin = (mcId-10<0) ? 0 : mcId-10 ;
1462  const Double_t mcIdMax = mcId + 10 ;
1463 
1464  hGeantIdMc->SetTitle(Form("MC Geant id for %s", getParticleName()));
1465  hGeantIdMc->SetMaximum(hGeantIdMc->GetMaximum()*1.2);
1466  hGeantIdMc->SetAxisRange(mcIdMin, mcIdMax, "X");
1467  hGeantIdMc->Draw();
1468 
1470  TLine* hGeantIdMcExpected = new TLine(mGeantId+0.5, 0, mGeantId+0.5, hGeantIdMc->GetMaximum()) ;
1471  hGeantIdMcExpected->SetLineColor(kBlue);
1472  hGeantIdMcExpected->Draw();
1473 
1475  topPad->cd(2);
1476  TH1* hGeantIdReco = (TH1D*) getHistogram(Form("hGeantId_%d", getCategoryId()));
1477  if(!hGeantIdReco) return kFALSE ;
1478 
1479  hGeantIdReco->SetTitle("Reconstructed geant id");
1480  hGeantIdReco->SetMaximum(hGeantIdReco->GetMaximum()*1.2);
1481  hGeantIdReco->SetAxisRange(0, 50, "X"); // Max geant id is 50
1482  hGeantIdReco->Draw();
1483 
1485  TLine** hGeantIdRecoExpected = new TLine*[mDaughterGeantId.size()];
1486 
1487  for(UInt_t id=0; id<mDaughterGeantId.size(); id++){
1488  const Int_t geantid = mDaughterGeantId[id] ;
1489  hGeantIdRecoExpected[id] = new TLine(geantid+0.5, 0, geantid+0.5, hGeantIdReco->GetMaximum()) ;
1490  hGeantIdRecoExpected[id]->SetLineColor(kRed);
1491  hGeantIdRecoExpected[id]->SetLineStyle(id+1);
1492  hGeantIdRecoExpected[id]->Draw();
1493  }
1494 
1495  mMainPad->cd(2);
1496  TLegend* leg = new TLegend(0.1, 0.2, 0.9, 0.8);
1497  leg->SetTextSize(0.05);
1498  leg->SetFillColor(10);
1499  leg->SetHeader("Particle informations");
1500 
1501  leg->AddEntry(hGeantIdMcExpected, getMcParticleName(), "L");
1502 
1503  for(UInt_t id=0; id<mDaughterGeantId.size(); id++){
1504  leg->AddEntry( hGeantIdRecoExpected[id], getEmbeddingParticleName(id), "L");
1505  }
1506  leg->Draw();
1507 
1508  mCanvas->cd();
1509  mCanvas->Update();
1510 
1511  mPDF->NewPage() ;
1512 
1513  delete header ;
1514 
1515  return kTRUE ;
1516 }
1517 
1518 //____________________________________________________________________________________________________
1519 Bool_t StEmbeddingQADraw::drawPhi() const
1520 {
1522 
1524  if(!isOpen()) return kFALSE ;
1525 
1527  Bool_t isPhiEmbeddingVsMcOk = drawProjection2D("phi", kTRUE);
1528 
1530  Bool_t isPhiEmbeddingVsRealOk = drawProjection2D("phi", kFALSE);
1531 
1532  return isPhiEmbeddingVsMcOk && isPhiEmbeddingVsRealOk ;
1533 }
1534 
1535 //____________________________________________________________________________________________________
1536 Bool_t StEmbeddingQADraw::drawRapidity() const
1537 {
1539 
1541  if(!isOpen()) return kFALSE ;
1542 
1544  const Bool_t isEtaOk = drawProjection2D("eta");
1545 
1547  const Bool_t isYOk = drawProjection2D("y");
1548 
1549  return isEtaOk && isYOk ;
1550 }
1551 
1552 //____________________________________________________________________________________________________
1553 Bool_t StEmbeddingQADraw::drawMomentum() const
1554 {
1556 
1558  if(!isOpen()) return kFALSE ;
1559 
1561  for(UInt_t id=0; id<mDaughterGeantId.size(); id++){
1562  TPaveText* header = 0 ;
1563 
1564  TH2* hRecoPVsMcP = (TH2D*) getHistogram("hRecoPVsMcP", id, kTRUE);
1565  if(!hRecoPVsMcP){
1566  // Draw error messages
1567  header = initCanvas("Reconstructed momentum vs MC momentum");
1568  drawErrorMessages(getHistogramName("hRecoPVsMcP", id, kTRUE)) ;
1569 
1570  mCanvas->cd();
1571  mCanvas->Update();
1572  mPDF->NewPage();
1573 
1574  delete header ;
1575 
1576  return kFALSE ;
1577  }
1578 
1579  header = initCanvas("Reconstructed momentum vs MC momentum");
1580 
1581  hRecoPVsMcP->SetAxisRange(0, mPtMax, "X");
1582  hRecoPVsMcP->SetAxisRange(0, mPtMax, "Y");
1583 
1584  hRecoPVsMcP->SetTitle(getParticleName(mDaughterGeantId[id]));
1585 
1586  hRecoPVsMcP->Draw("colz");
1587 
1588  mCanvas->cd();
1589  mCanvas->Update();
1590 
1591  mPDF->NewPage() ;
1592 
1593  delete header ;
1594  }
1595 
1596  return drawProjection2D("momentum");
1597 }
1598 
1599 //____________________________________________________________________________________________________
1600 Bool_t StEmbeddingQADraw::drawPt() const
1601 {
1603 
1605  if(!isOpen()) return kFALSE ;
1606 
1607  // 1/pT(Gl) - 1/pT(MC) vs pT(MC) (Embedding only)
1608  for(UInt_t id=0; id<mDaughterGeantId.size(); id++){
1609  TPaveText* header = 0 ;
1610 
1611  TH2* hdInvPtVsPt = (TH2D*) getHistogram("hdInvPtVsPt", id, kTRUE);
1612  if(!hdInvPtVsPt){
1613  // Draw error messages
1614  header = initCanvas(Form("1/p_{T} (Gl) - 1/p_{T} (MC) vs p_{T} (MC) (%s)", getParticleName(mDaughterGeantId[id])));
1615  drawErrorMessages(getHistogramName("hdInvPtVsPt", id, kTRUE)) ;
1616 
1617  mCanvas->cd();
1618  mCanvas->Update();
1619  mPDF->NewPage();
1620 
1621  delete header ;
1622 
1623  return kFALSE ;
1624  }
1625 
1626  gStyle->SetPadRightMargin(0.15);
1627 
1628  header = initCanvas(Form("1/p_{T} (Gl) - 1/p_{T} (MC) vs p_{T} (MC) (%s)", getParticleName(mDaughterGeantId[id])), 1, 2);
1629  // 2D
1630  mMainPad->cd(1);
1631  hdInvPtVsPt->SetAxisRange(0, mPtMax, "X");
1632  hdInvPtVsPt->SetTitle("");
1633  hdInvPtVsPt->Draw("colz");
1634 
1635  // Profile
1636  mMainPad->cd(2);
1637  TProfile* pdInvPtVsPt = (TProfile*) hdInvPtVsPt->ProfileX(Form("p%s", hdInvPtVsPt->GetName()));
1638  pdInvPtVsPt->SetAxisRange(0, mPtMax, "X");
1639  pdInvPtVsPt->SetTitle("");
1640  pdInvPtVsPt->SetYTitle(hdInvPtVsPt->GetYaxis()->GetTitle());
1641  pdInvPtVsPt->SetLineWidth(2);
1642  pdInvPtVsPt->SetMinimum(hdInvPtVsPt->GetYaxis()->GetXmin());
1643  pdInvPtVsPt->SetMaximum(hdInvPtVsPt->GetYaxis()->GetXmax());
1644  pdInvPtVsPt->Draw();
1645 
1646  mCanvas->cd();
1647  mCanvas->Update();
1648 
1649  mPDF->NewPage() ;
1650 
1651  delete header ;
1652  }
1653 
1654  gStyle->SetPadRightMargin(0.05);
1655 
1656  return drawProjection2D("pt");
1657 }
1658 
1659 //____________________________________________________________________________________________________
1660 Bool_t StEmbeddingQADraw::drawProjection2D(const TString name, const Bool_t isMC) const
1661 {
1663 
1666  //
1676 
1677  TString nameLower(name);
1678  nameLower.ToLower();
1679 
1680  LOG_INFO << "QA for " << name << " ..." << endm;
1681 
1683  TString histoName("");
1684  Bool_t isProjectionX = kFALSE ;
1685  TString yTitle("");
1686  TString xTitle("");
1687 
1688  if( nameLower.Contains("pt") ){
1689  xTitle = "#eta";
1690  yTitle = "p_{T}";
1691  histoName = "hPtVsEta";
1692  isProjectionX = kFALSE ;
1693  }
1694  else if( nameLower.Contains("momentum") ){
1695  xTitle = "#eta";
1696  yTitle = "p";
1697  histoName = "hMomVsEta";
1698  isProjectionX = kFALSE ;
1699  }
1700  else if( nameLower.Contains("eta") || name.Contains("pseudorapidity") ){
1701  xTitle = "p_{T}";
1702  yTitle = "#eta";
1703  histoName = "hPtVsEta";
1704  isProjectionX = kTRUE ;
1705  }
1706  else if( nameLower.Contains("y") || name.Contains("rapidity") ){
1707  xTitle = "p_{T}";
1708  yTitle = "y";
1709  histoName = "hPtVsY";
1710  isProjectionX = kTRUE ;
1711  }
1712  else if( nameLower.Contains("phi")){
1713  xTitle = "p_{T}";
1714  yTitle = "#phi";
1715  histoName = "hPtVsPhi";
1716  isProjectionX = kTRUE ;
1717  }
1718  else{
1719  Error("DrawProjection2D", "Unknown variable, %s", name.Data());
1720  LOG_INFO << " Current implemented variables are" << endm;
1721  LOG_INFO << "-----------------------------------------------------------------" << endm;
1722  LOG_INFO << " Input variable" << endm;
1723  LOG_INFO << "-----------------------------------------------------------------" << endm;
1724  LOG_INFO << " pt p_{T}" << endm;
1725  LOG_INFO << " momentum p" << endm;
1726  LOG_INFO << " eta or pseudorapidity eta" << endm;
1727  LOG_INFO << " y or rapidity y" << endm;
1728  LOG_INFO << " phi phi" << endm;
1729  LOG_INFO << "-----------------------------------------------------------------" << endm;
1730  LOG_INFO << endm;
1731  LOG_INFO << "NOTE : Input is case insensitive" << endm;
1732  LOG_INFO << endm;
1733 
1734  return kFALSE ;
1735  }
1736 
1737  gStyle->SetPadRightMargin(0.05);
1738 
1739  // MC tracks
1740  TH2* h2DMc = (TH2D*) getHistogram(Form("%s_0_%d", histoName.Data(), mGeantId));
1741 
1742  for(UInt_t id=0; id<mDaughterGeantId.size(); id++){
1743  const TString headerTitle(Form("Projection of %s for each %s bin", yTitle.Data(), xTitle.Data()));
1744  TPaveText* header = 0 ;
1745 
1746  TH2* h2DEmbed = (TH2D*) getHistogram(histoName, id, kTRUE);
1747  if(!h2DEmbed){
1748  // Draw error message
1749  header = initCanvas(headerTitle);
1750  drawErrorMessages(getHistogramName(histoName, id, kTRUE)) ;
1751 
1752  mCanvas->cd();
1753  mCanvas->Update();
1754  mPDF->NewPage();
1755 
1756  delete header ;
1757 
1758  return kFALSE ;
1759  }
1760 
1761  TH2* h2DReal = (TH2D*) getHistogram(histoName, id, kFALSE);
1762 
1764  header = initCanvas(headerTitle);
1765 
1766  Int_t npad = 0 ;
1767  Int_t npadMax = 0 ;
1768  if( isProjectionX ){
1769  mMainPad->Divide(2, 3);
1770  npad = 5 ;
1771  npadMax = 6 ;
1772  }
1773  else{
1774  mMainPad->Divide(2, 4);
1775  npad = 6 ;
1776  npadMax = 8 ;
1777  }
1778 
1782  const Int_t nbins = (isProjectionX) ? 10 : 6 ;
1783  const Double_t binStep = 0.5 ;
1784  const Double_t binMin = (isProjectionX) ? 0.0 : -1.5 ;
1785 
1786  Int_t ipad = 1 ;
1787  for(Int_t ibin=0; ibin<nbins; ibin++){
1788  if( ipad%(npad+1) == 0 ){
1789  mCanvas->cd();
1790  mCanvas->Update();
1791  mPDF->NewPage();
1792 
1793  if(header) delete header ;
1794  header = initCanvas(headerTitle);
1795 
1796  if( isProjectionX ) mMainPad->Divide(2, 3);
1797  else mMainPad->Divide(2, 4);
1798  ipad = 1 ;
1799  }
1800 
1801  const Double_t xMin = (isProjectionX && ibin==0) ? 0.2 : binMin + binStep * ibin ;
1802  const Double_t xMax = (isProjectionX && ibin==0) ? 0.5 : binMin + binStep * (ibin+1.0) ;
1803  const Int_t xMinBin = (isProjectionX) ? h2DEmbed->GetYaxis()->FindBin(xMin) : h2DEmbed->GetXaxis()->FindBin(xMin) ;
1804  const Int_t xMaxBin = (isProjectionX) ? h2DEmbed->GetYaxis()->FindBin(xMax) - 1 : h2DEmbed->GetXaxis()->FindBin(xMax) - 1;
1805 
1806  TH1* hEmbed = 0;
1807  TH1* hReal = 0;
1808  TH1* hMc = 0;
1809  if( isProjectionX ){
1810  if( h2DEmbed ) hEmbed = (TH1D*) h2DEmbed->ProjectionX(Form("h%sEmbed_%d_%d", name.Data(), id, ibin), xMinBin, xMaxBin);
1811  if( h2DReal ) hReal = (TH1D*) h2DReal->ProjectionX(Form("h%sReal_%d_%d", name.Data(), id, ibin), xMinBin, xMaxBin);
1812  if( h2DMc ) hMc = (TH1D*) h2DMc->ProjectionX(Form("h%sMc_%d_%d", name.Data(), id, ibin), xMinBin, xMaxBin);
1813  }
1814  else{
1815  if( h2DEmbed ) hEmbed = (TH1D*) h2DEmbed->ProjectionY(Form("h%sEmbed_%d_%d", name.Data(), id, ibin), xMinBin, xMaxBin);
1816  if( h2DReal ) hReal = (TH1D*) h2DReal->ProjectionY(Form("h%sReal_%d_%d", name.Data(), id, ibin), xMinBin, xMaxBin);
1817  if( h2DMc ) hMc = (TH1D*) h2DMc->ProjectionY(Form("h%sMc_%d_%d", name.Data(), id, ibin), xMinBin, xMaxBin);
1818  }
1819 
1820  hEmbed->SetLineColor(kRed);
1821  if(hEmbed->GetSumw2N()==0) hEmbed->Sumw2();
1822  hEmbed->Scale( getNormalization(*hEmbed) );
1823 
1824  if( hReal ){
1825  hReal->SetLineColor(kBlue);
1826  if(hReal->GetSumw2N()==0) hReal->Sumw2();
1827  hReal->Scale( getNormalization(*hReal) );
1828  }
1829 
1830  if( hMc ){
1831  hMc->SetLineColor(kBlack);
1832  if(hMc->GetSumw2N()==0) hMc->Sumw2();
1833  hMc->Scale( getNormalization(*hMc) );
1834  }
1835 
1836  hEmbed->SetMinimum(0.0);
1837 
1838  if( isProjectionX ){
1839  hEmbed->SetTitle(Form("%1.1f < p_{T} < %1.1f (GeV/c)", xMin, xMax));
1840  }
1841  else{
1842  hEmbed->SetTitle(Form("%1.1f < #eta < %1.1f", xMin, xMax));
1843  hEmbed->SetAxisRange(0, mPtMax, "X");
1844  }
1845  hEmbed->SetYTitle(Form("(1/N_{trk})dN/d%s", yTitle.Data()));
1846 
1847  mMainPad->cd(ipad) ;
1848  hEmbed->Draw();
1849 
1851  if(isMC && !isDecay()){
1853  hEmbed->SetMaximum( TMath::Max(hMc->GetMaximum(), hEmbed->GetMaximum()) * 1.2 );
1854 
1855  hMc->Draw("same");
1856  }
1857  else{
1861  const Double_t yMax = (hReal) ? TMath::Max(hReal->GetMaximum(), hEmbed->GetMaximum()) : hEmbed->GetMaximum() ;
1862  hEmbed->SetMaximum( yMax * 1.2 );
1863 
1864  if(hReal) hReal->Draw("hsame");
1865  }
1866  hEmbed->Draw("same");
1867 
1868  // Draw legend in the last pad in each Canvas
1869  if ( ipad == 1 ){
1870  mMainPad->cd(npadMax);
1871  TLegend* leg = new TLegend(0.0, 0.2, 1.0, 0.8);
1872  leg->SetFillColor(10);
1873  leg->SetTextFont(43);
1874  leg->SetTextSize(15);
1875 
1876  leg->AddEntry(hEmbed, getEmbeddingParticleName(id, kTRUE), "L");
1877  if(isMC && !isDecay()){
1878  leg->AddEntry(hMc, getMcParticleName(), "L");
1879  }
1880  else{
1881  if(hReal) leg->AddEntry(hReal, getRealParticleName(id, kTRUE), "L");
1882  }
1883  leg->Draw();
1884  }
1885 
1886  ipad++;
1887  }
1888 
1889  mCanvas->cd();
1890  mCanvas->Update();
1891 
1892  mPDF->NewPage();
1893 
1894  if(header) delete header ;
1895  }
1896 
1897  gStyle->SetPadRightMargin(0.15);
1898 
1899  return kTRUE ;
1900 }
1901 
1902 //____________________________________________________________________________________________________
1903 Bool_t StEmbeddingQADraw::drawdEdx() const
1904 {
1906 
1908  if(!isOpen()) return kFALSE ;
1909 
1911  const Bool_t isMcOk = drawdEdxVsMomentum(kTRUE) ;
1912 
1914  const Bool_t isRecoOk = drawdEdxVsMomentum(kFALSE) ;
1915 
1916  return isMcOk && isRecoOk ;
1917 }
1918 
1919 //____________________________________________________________________________________________________
1920 Bool_t StEmbeddingQADraw::drawdEdxVsMomentum(const Bool_t isMcMomentum) const
1921 {
1925  //
1927 
1930  const TString momName = (isMcMomentum) ? "Mc" : "Reco" ;
1931 
1932  LOG_INFO << "QA for dE/dx (" << momName << " momentum ...)" << endm;
1933 
1934  gStyle->SetOptFit(0);
1935  gStyle->SetPadRightMargin(0.05);
1936 
1937  for(UInt_t id=0; id<mDaughterGeantId.size(); id++){
1938  //----------------------------------------------------------------------------------------------------
1939  // 2D
1940  //----------------------------------------------------------------------------------------------------
1941  TString headerTitle("");
1942  if( mIsEmbeddingOnly ){
1943  headerTitle = Form("dE/dx vs momentum (Embedding:%s)", getParticleName(mDaughterGeantId[id]) );
1944  }
1945  else{
1946  headerTitle = Form("dE/dx vs momentum (Embedding:%s, Real:%s)",
1947  getParticleName(mDaughterGeantId[id]), getParticleName(getGeantIdReal(id)) );
1948  }
1949  TPaveText* header = 0 ;
1950 
1951  // Embedding
1952  TH2* hdEdxVsMomEmbed = (TH2D*) getHistogram(Form("hdEdxVsMom%s", momName.Data()), id, kTRUE);
1953  if(!hdEdxVsMomEmbed){
1954  // Draw error messages
1955  header = initCanvas(headerTitle);
1956  drawErrorMessages(getHistogramName(Form("hdEdxVsMom%s", momName.Data()), id, kTRUE)) ;
1957 
1958  mCanvas->cd();
1959  mCanvas->Update();
1960  mPDF->NewPage();
1961 
1962  delete header ;
1963 
1964  return kFALSE ;
1965  }
1966 
1967  header = initCanvas(headerTitle, 1, 2);
1968 
1969  hdEdxVsMomEmbed->SetLineColor(kRed);
1970  hdEdxVsMomEmbed->SetMarkerColor(kRed);
1971 
1972  // Real data
1973  TH2* hdEdxVsMomReal = (TH2D*) getHistogram("hdEdxVsMomReco", id, kFALSE);
1974 // if(!hdEdxVsMomReal) return kFALSE ;
1975 
1976  if( hdEdxVsMomReal ){
1977  hdEdxVsMomReal->SetLineColor(kBlack);
1978  hdEdxVsMomReal->SetMarkerColor(kBlack);
1979  }
1980 
1981  // Real data (with PID)
1982  TH2* hdEdxVsMomPidReal = (TH2D*) getHistogram("hdEdxVsMomRecoPidCut", id, kFALSE);
1983 // if(!hdEdxVsMomPidReal) return kFALSE ;
1984 
1985  if( hdEdxVsMomPidReal ){
1986  hdEdxVsMomPidReal->SetLineColor(kBlue);
1987  hdEdxVsMomPidReal->SetMarkerColor(kBlue);
1988  }
1989 
1990  // dE/dx vs momentum
1991  mMainPad->cd(1);
1992  hdEdxVsMomEmbed->SetAxisRange(0, mPtMax, "X");
1993  hdEdxVsMomEmbed->SetXTitle(Form("%s momentum (GeV/c)", momName.Data()));
1994  hdEdxVsMomEmbed->SetTitle("");
1995  hdEdxVsMomEmbed->SetMinimum(0);
1996 
1997  hdEdxVsMomEmbed->Draw("box");
1998  if( hdEdxVsMomReal ) hdEdxVsMomReal->Draw("box same");
1999  if( hdEdxVsMomPidReal ) hdEdxVsMomPidReal->Draw("box same");
2000  hdEdxVsMomEmbed->Draw("box same");
2001 
2002  mMainPad->cd(2);
2003 
2004  const StEmbeddingQAUtilities* utility = StEmbeddingQAUtilities::instance() ;
2005  TLegend* leg = new TLegend(0.1, 0.25, 0.9, 0.9);
2006  leg->SetFillColor(10);
2007  leg->SetTextSize(0.05);
2008  leg->AddEntry( hdEdxVsMomEmbed, getEmbeddingParticleName(id), "L");
2009  if( hdEdxVsMomReal ) leg->AddEntry( hdEdxVsMomReal, "Real data", "L");
2010  if( hdEdxVsMomPidReal ) leg->AddEntry( hdEdxVsMomPidReal, Form("Real data with PID cut (#sigma<2) %s",utility->getBTofPid()?"BTOF":"TPC"), "L");
2011  leg->Draw();
2012 
2013  mCanvas->cd();
2014  mCanvas->Update();
2015 
2016  mPDF->NewPage();
2017 
2018  delete header ;
2019 
2020  //----------------------------------------------------------------------------------------------------
2023  //----------------------------------------------------------------------------------------------------
2024  headerTitle = "Projection of dE/dx for each p bin";
2025 // header = initCanvas(headerTitle, 3, 3);
2026  header = initCanvas(headerTitle, 2, 4);
2027  Int_t ipad = 1 ;
2028 // const Int_t npad = 8 ;
2029 // const Int_t npadMax = 9;
2030  const Int_t npad = 6 ;
2031  const Int_t npadMax = 8;
2032 
2033  TGraphErrors* gMeanVsMom[2]; // 0:embedding, 1:real data
2034  TGraphErrors* gSigmaVsMom[2]; // 0:embedding, 1:real data
2035 
2036  for(Int_t i=0; i<2; i++){
2037  gMeanVsMom[i] = new TGraphErrors();
2038  gSigmaVsMom[i] = new TGraphErrors();
2039  gMeanVsMom[i] ->SetMarkerStyle(24 - i*4);
2040  gSigmaVsMom[i]->SetMarkerStyle(24 - i*4);
2041  }
2042  gMeanVsMom[0] ->SetMarkerColor(kRed); gMeanVsMom[0] ->SetLineColor(kRed);
2043  gSigmaVsMom[0]->SetMarkerColor(kRed); gSigmaVsMom[0]->SetLineColor(kRed);
2044  gMeanVsMom[1] ->SetMarkerColor(kBlue); gMeanVsMom[1] ->SetLineColor(kBlue);
2045  gSigmaVsMom[1]->SetMarkerColor(kBlue); gSigmaVsMom[1]->SetLineColor(kBlue);
2046 
2047  const Int_t npt = 24 ;
2048  const Double_t ptBin = 0.2 ;
2049  for(Int_t ipt=0; ipt<npt; ipt++){
2050  if( ipad % (npad+1) == 0 ){
2051  mCanvas->cd();
2052  mCanvas->Update();
2053 
2054  mPDF->NewPage();
2055 
2056  if(header) delete header ;
2057  header = initCanvas(headerTitle, 2, 4);
2058  ipad = 1 ;
2059  }
2060 
2061  const Double_t ptMin = 0.2 + ipt*ptBin ;
2062  const Double_t ptMax = ptMin + ptBin ;
2063  const Int_t ptMinBin = hdEdxVsMomEmbed->GetXaxis()->FindBin(ptMin);
2064  const Int_t ptMaxBin = hdEdxVsMomEmbed->GetXaxis()->FindBin(ptMax-0.001);
2065  if( ptMinBin == ptMaxBin ){
2066  LOG_INFO << Form("%1.1f - %1.1f GeV/c : bin = (%4d, %4d)", ptMin, ptMax, ptMinBin, ptMaxBin) << endm;
2067  }
2068 
2069  // Projections
2070  TH1* hdEdxEmbed = (TH1D*) hdEdxVsMomEmbed->ProjectionY(Form("hdEdxEmbed%s_%d_%d", momName.Data(), id, ipt), ptMinBin, ptMaxBin);
2071  TH1* hdEdxReal = 0;
2072  if( hdEdxVsMomPidReal )
2073  hdEdxReal = (TH1D*) hdEdxVsMomPidReal->ProjectionY(Form("hdEdxReal%s_%d_%d", momName.Data(), id, ipt), ptMinBin, ptMaxBin);
2074 
2075  hdEdxEmbed->Sumw2() ;
2076  hdEdxEmbed->Scale( getNormalization(*hdEdxEmbed) ) ;
2077 
2078  if( hdEdxReal ){
2079  hdEdxReal->Sumw2() ;
2080  hdEdxReal->Scale( getNormalization(*hdEdxReal) ) ;
2081  }
2082 
2083  hdEdxEmbed->SetMinimum(0.0);
2084 
2085  // Use max(real, embed) if real exists
2086  // Use embed if not
2087  if( hdEdxReal ){
2088  hdEdxEmbed->SetMaximum( TMath::Max(hdEdxReal->GetMaximum(), hdEdxEmbed->GetMaximum())*1.2 );
2089  }
2090  else{
2091  hdEdxEmbed->SetMaximum( hdEdxEmbed->GetMaximum()*1.2 );
2092  }
2093 
2094  hdEdxEmbed->SetTitle(Form("%1.1f < %s p < %1.1f GeV/c", ptMin, momName.Data(), ptMax));
2095  hdEdxEmbed->SetYTitle("(1/N_{trk})dN/d(dE/dx)");
2096 
2097  mMainPad->cd(ipad);
2098  hdEdxEmbed->Draw("h");
2099  if( hdEdxReal ) hdEdxReal->Draw("hsame");
2100  hdEdxEmbed->Draw("hsame");
2101 
2102  //----------------------------------------------------------------------------------------------------
2103  // Extract mean and sigma (Oct/22/2009)
2104  //----------------------------------------------------------------------------------------------------
2105  TF1 fEmbed(Form("fEmbed%s_%d_%d", momName.Data(), id, ipt), "gaus", 0, 10);
2106  TF1 fReal(Form("fReal%s_%d_%d", momName.Data(), id, ipt), "gaus", 0, 10);
2107  fEmbed.SetLineColor(kRed);
2108  fReal.SetLineColor(kBlue);
2109  fEmbed.SetLineWidth(1);
2110  fReal.SetLineWidth(1);
2111 
2112  if( hdEdxReal ){
2113  hdEdxReal->Fit(fReal.GetName(), "rq0");
2114  fReal.Draw("same");
2115  }
2116  hdEdxEmbed->Fit(fEmbed.GetName(), "rq0");
2117  fEmbed.Draw("same");
2118 
2119  const Double_t pt = (ptMin+ptMax)/2.0 ;
2120  gMeanVsMom[0]->SetPoint(ipt, pt, fEmbed.GetParameter(1));
2121  gMeanVsMom[0]->SetPointError(ipt, 0.0, fEmbed.GetParError(1));
2122  gSigmaVsMom[0]->SetPoint(ipt, pt, fEmbed.GetParameter(2));
2123  gSigmaVsMom[0]->SetPointError(ipt, 0.0, fEmbed.GetParError(2));
2124 
2125  if( hdEdxReal ){
2126  gMeanVsMom[1]->SetPoint(ipt, pt, fReal.GetParameter(1));
2127  gMeanVsMom[1]->SetPointError(ipt, 0.0, fReal.GetParError(1));
2128  gSigmaVsMom[1]->SetPoint(ipt, pt, fReal.GetParameter(2));
2129  gSigmaVsMom[1]->SetPointError(ipt, 0.0, fReal.GetParError(2));
2130  }
2131  else{
2132  gMeanVsMom[1]->SetPoint(ipt, pt, -9999.);
2133  gMeanVsMom[1]->SetPointError(ipt, 0.0, 0.0);
2134  gSigmaVsMom[1]->SetPoint(ipt, pt, -9999.);
2135  gSigmaVsMom[1]->SetPointError(ipt, 0.0, 0.0);
2136  }
2137 
2138  // Draw legend in the last pad in each Canvas
2139  if( ipad == 1 ){
2140  mMainPad->cd(npadMax);
2141  drawLegend(id, hdEdxEmbed, hdEdxReal, "L", kTRUE) ;
2142 // TLegend* leg = new TLegend(0, 0.2, 1, 0.8);
2143 // leg->SetFillColor(10);
2144 // leg->SetTextFont(43);
2145 // leg->SetTextSize(15);
2146 // leg->AddEntry( hdEdxEmbed, getEmbeddingParticleName(id, kTRUE), "L");
2147 // if( hdEdxReal ) leg->AddEntry( hdEdxReal, getRealParticleName(id, kTRUE), "L");
2148 // leg->Draw();
2149  }
2150 
2151  ipad++;
2152  }// pt loop
2153 
2154  mCanvas->cd();
2155  mCanvas->Update();
2156 
2157  mPDF->NewPage();
2158 
2159  delete header ;
2160 
2161  //----------------------------------------------------------------------------------------------------
2163  //----------------------------------------------------------------------------------------------------
2164  headerTitle = "Mean/#sigma of dE/dx vs momentum";
2165  header = initCanvas(headerTitle, 1, 2);
2166 
2167  for(Int_t i=0; i<2; i++){
2168  mMainPad->cd(i+1);
2169 
2170  const Double_t ymax = (i==0) ? 7.2 : 1.4 ;
2171  TH1* frame = mMainPad->GetPad(i+1)->DrawFrame(0, 0.0, 5.0, ymax);
2172  frame->SetXTitle(Form("%s momentum (GeV/c)", momName.Data()));
2173  if( i == 0 ) frame->SetYTitle("Mean (keV/cm)");
2174  if( i == 1 ) frame->SetYTitle("#sigma (keV/cm)");
2175 
2176  TLegend* leg = new TLegend(0.23, 0.86, 0.92, 0.99);
2177  leg->SetBorderSize(1);
2178  leg->SetTextFont(43);
2179  leg->SetTextSize(15);
2180  leg->SetFillColor(10);
2181 
2182  if( i == 0 ){
2183  gMeanVsMom[0]->Draw("P");
2184  gMeanVsMom[1]->Draw("P");
2185 
2186  leg->AddEntry( gMeanVsMom[0], getEmbeddingParticleName(id), "P");
2187  leg->AddEntry( gMeanVsMom[1], getRealParticleName(id), "P");
2188  }
2189  if( i == 1 ){
2190  gSigmaVsMom[0]->Draw("P");
2191  gSigmaVsMom[1]->Draw("P");
2192 
2193  leg->AddEntry( gSigmaVsMom[0], getEmbeddingParticleName(id), "P");
2194  leg->AddEntry( gSigmaVsMom[1], getRealParticleName(id), "P");
2195  }
2196  leg->Draw();
2197  }
2198 
2199  mCanvas->cd();
2200  mCanvas->Update();
2201 
2202  mPDF->NewPage() ;
2203 
2204  delete header ;
2205  }
2206 
2207  return kTRUE ;
2208 }
2209 
2210 //____________________________________________________________________________________________________
2211 Bool_t StEmbeddingQADraw::drawDca() const
2212 {
2214 
2215  LOG_INFO << "QA for dca distributions ..." << endm;
2216 
2218  if(!isOpen()) return kFALSE ;
2219 
2220  return drawProjection3D("Dca");
2221 }
2222 
2223 //____________________________________________________________________________________________________
2224 Bool_t StEmbeddingQADraw::drawNHit() const
2225 {
2227 
2228  LOG_INFO << "QA for NHit distributions ..." << endm;
2229 
2231  if(!isOpen()) return kFALSE ;
2232 
2233  gStyle->SetPadRightMargin(0.17);
2234 
2236  for(UInt_t id=0; id<mDaughterGeantId.size(); id++){
2237  TH3* hNCommonHitVsNHit = (TH3D*) getHistogram("hNCommonHitVsNHit", id, kTRUE);
2238  if ( !hNCommonHitVsNHit ) continue ;
2239 
2240  TString headerTitle("");
2241  if( mIsEmbeddingOnly ){
2242  headerTitle = Form("N_{common} vs N_{hit} (Embedding:%s)", getParticleName(mDaughterGeantId[id]) );
2243  }
2244  else{
2245  headerTitle = Form("N_{common} vs N_{hit} (Embedding:%s, Real:%s)",
2246  getParticleName(mDaughterGeantId[id]), getParticleName(getGeantIdReal(id)) );
2247  }
2248 
2249  // pt integrated
2250  TPaveText* header = initCanvas(headerTitle);
2251 
2252  // Maximum pt = 5 GeV/c for ncommon hit histogram
2253  hNCommonHitVsNHit->SetAxisRange(0.2, 5.0);
2254  TH2* hNCommonHitVsNHit2D = (TH2D*) hNCommonHitVsNHit->Project3D("zy");
2255  hNCommonHitVsNHit2D->SetName(Form("hNCommonHitVsNHit2D_%d", id));
2256  hNCommonHitVsNHit2D->SetTitle(Form("%1.1f < p_{T} < %1.1f GeV/c", 0.2, 5.0));
2257  hNCommonHitVsNHit2D->Draw("colz");
2258 
2259  mCanvas->cd();
2260  mCanvas->Update();
2261  mPDF->NewPage();
2262  delete header ;
2263 
2264  // Slices for each pt bin (0.5 GeV/c step)
2265  if( mIsEmbeddingOnly ){
2266  headerTitle = Form("N_{common} vs N_{hit}, p_{T} dependence (Embedding:%s)", getParticleName(mDaughterGeantId[id]) );
2267  }
2268  else{
2269  headerTitle = Form("N_{common} vs N_{hit}, p_{T} dependence (Embedding:%s, Real:%s)",
2270  getParticleName(mDaughterGeantId[id]), getParticleName(getGeantIdReal(id)) );
2271  }
2272  for(Int_t jpt=0; jpt<2; jpt++) {
2273  TPaveText* header = initCanvas(headerTitle, 2, 3);
2274 
2275  const Int_t npt = hNCommonHitVsNHit->GetNbinsX()/2 ;
2276  for(Int_t ipt=0; ipt<npt; ipt++) {
2277  mMainPad->cd(ipt+1);
2278  const Int_t ptId = jpt*npt + ipt + 1;
2279  Double_t ptmin = hNCommonHitVsNHit->GetXaxis()->GetBinLowEdge(ptId) ;
2280  Double_t ptmax = hNCommonHitVsNHit->GetXaxis()->GetBinLowEdge(ptId+1) ;
2281  if( ptmin == 0.0 ) ptmin = 0.2 ;
2282 
2283  hNCommonHitVsNHit->SetAxisRange(ptmin, ptmax);
2284  hNCommonHitVsNHit2D = (TH2D*) hNCommonHitVsNHit->Project3D("zy");
2285  hNCommonHitVsNHit2D->SetName(Form("hNCommonHitVsNHit2D_%d_%d_%d", id, ipt, jpt));
2286  hNCommonHitVsNHit2D->SetTitle(Form("%1.1f < p_{T} < %1.1f GeV/c", ptmin, ptmax));
2287  hNCommonHitVsNHit2D->Draw("colz");
2288  }
2289  mCanvas->cd();
2290  mCanvas->Update();
2291  mPDF->NewPage();
2292  delete header ;
2293  }
2294  }
2295 
2296  gStyle->SetPadRightMargin(0.05);
2297 
2298  return drawProjection3D("NHit") ;
2299 }
2300 
2301 //____________________________________________________________________________________________________
2302 Bool_t StEmbeddingQADraw::drawProjection3D(const TString name) const
2303 {
2305 
2308 
2310  const Bool_t isBatch = gROOT->IsBatch() ;
2311  if( !isBatch ){
2312  LOG_INFO << "Enter batch mode ..." << endm;
2313  gROOT->SetBatch(kTRUE);
2314  }
2315 
2316  TString nameLower(name);
2317  nameLower.ToLower();
2318 
2319  TString headerTitle(Form("%s distribution for (p_{T}, #eta) slices", name.Data()));
2320  for(UInt_t id=0; id<mDaughterGeantId.size(); id++){
2321  TPaveText* header = 0 ;
2322 
2324  TH3* h3DEmbed = (TH3D*) getHistogram(Form("h%s", name.Data()), id, kTRUE);
2325  if(!h3DEmbed){
2326  // Draw error messages
2327  header = initCanvas(headerTitle);
2328  drawErrorMessages(getHistogramName(Form("h%s", name.Data()), id, kTRUE));
2329 
2330  mCanvas->cd();
2331  mCanvas->Update();
2332  mPDF->NewPage();
2333 
2334  delete header ;
2335 
2336  return kFALSE ;
2337  }
2338 
2339  TH3* h3DReal = (TH3D*) getHistogram(Form("h%s", name.Data()), id, kFALSE);
2340 
2341  const Int_t nPt = h3DEmbed->GetNbinsX() ;
2342  const Int_t nEta = h3DEmbed->GetNbinsY() ;
2343  const Double_t ptMin = h3DEmbed->GetXaxis()->GetXmin() ;
2344  const Double_t ptMax = h3DEmbed->GetXaxis()->GetXmax() ;
2345  const Double_t etaMin = h3DEmbed->GetYaxis()->GetXmin() ;
2346  const Double_t etaMax = h3DEmbed->GetYaxis()->GetXmax() ;
2347  const Double_t etaBin = (etaMax-etaMin)/(Double_t)nEta ;
2348  const Double_t ptBin = (ptMax-ptMin)/(Double_t)nPt ;
2349 
2350  for(Int_t ipt=0; ipt<nPt; ipt++){
2351  TString pt(Form("%1.1f < p_{T} < %1.1f (GeV/c)", ptMin+ipt*ptBin, ptMin+(ipt+1)*ptBin));
2352  if( ipt == 0 ) pt = Form("%1.1f < p_{T} < %1.1f (GeV/c)", 0.1, ptMin+(ipt+1)*ptBin);
2353 
2354  TPaveText* header = initCanvas(headerTitle, 2, 3);
2355 
2356  const Int_t npad = 5 ;
2357  const Int_t npadMax = 6 ;
2358  Int_t ipad = 1 ;
2359  for(Int_t ieta=0; ieta<nEta; ieta++){
2360  if( ipad % (npad+1) == 0 ){
2361  mCanvas->cd();
2362  mCanvas->Update();
2363  mPDF->NewPage();
2364 
2365  if(header) delete header ;
2366  header = initCanvas(headerTitle, 2, 3);
2367 
2368  ipad = 1 ;
2369  }
2370 
2371  TString eta(Form("%1.1f < #eta < %1.1f", etaMin+ieta*etaBin, etaMin+(ieta+1)*etaBin));
2372 
2373  TH1* hEmbed = (TH1D*) h3DEmbed->ProjectionZ(Form("h%sEmbed_%d_%d_%d", name.Data(), id, ipt, ieta), ipt+1, ipt+1, ieta+1, ieta+1);
2374  TH1* hReal = 0 ;
2375  if( h3DReal ){
2376  hReal = (TH1D*) h3DReal->ProjectionZ(Form("h%sReal_%d_%d_%d", name.Data(), id, ipt, ieta), ipt+1, ipt+1, ieta+1, ieta+1);
2377  hReal ->Sumw2();
2378  hReal ->Scale( getNormalization(*hReal) );
2379  hReal ->SetLineColor(kBlue);
2380  }
2381  hEmbed->Sumw2();
2382  hEmbed->Scale( getNormalization(*hEmbed) );
2383  hEmbed->SetLineColor(kRed);
2384  hEmbed->SetMinimum(0.0);
2385 
2386  // Set maximum
2387  // If real data exists, use max(real, embed)
2388  // If not, use embedding histogram
2389  if( hReal ){
2390  hEmbed->SetMaximum( TMath::Max(hReal->GetMaximum(), hEmbed->GetMaximum()) * 1.2 );
2391  }
2392  else{
2393  hEmbed->SetMaximum( hEmbed->GetMaximum() * 1.2 );
2394  }
2395 
2396  hEmbed->SetTitle(eta + ", " + pt);
2397  hEmbed->SetYTitle(Form("(1/N_{trk})dN/d%s", name.Data())) ;
2398 
2399  mMainPad->cd(ipad);
2400  hEmbed->Draw();
2401  if(hReal) hReal->Draw("hsame");
2402  hEmbed->Draw("same");
2403 
2404  // Draw legend in the last pad for each Canvas
2405  if( ipad == 1 ){
2406  mMainPad->cd(npadMax);
2407  drawLegend(id, hEmbed, hReal, "L", kTRUE);
2408  }
2409 
2410  ipad++;
2411  }// eta loop
2412  mCanvas->cd();
2413  mCanvas->Update();
2414  mPDF->NewPage();
2415 
2416  delete header ;
2417  }// pt loop
2418  }// daughter particle loop
2419 
2420 // mCanvas->cd();
2421 // mCanvas->Update();
2422 //
2423 // mPDF->NewPage() ;
2424 
2425  // Back to the original mode
2426  if( isBatch ) gROOT->SetBatch(kTRUE); // Stay batch mode if you've started the macro by batch mode
2427  else gROOT->SetBatch(kFALSE); // Exit batch mode if you've not
2428 
2429  return kTRUE ;
2430 }
2431 
2432 //____________________________________________________________________________________________________
2433 TPaveText* StEmbeddingQADraw::initCanvas(const TString headerTitle, const Int_t nx, const Int_t ny) const
2434 {
2435  mMainPad->Clear();
2436  mCanvas->cd();
2437 
2439  TPaveText* header = drawHeader(headerTitle);
2440 
2441  mMainPad->cd();
2442  if( nx != 0 && ny != 0 ) mMainPad->Divide(nx, ny);
2443 
2444  return header ;
2445 }
2446 
2447 //____________________________________________________________________________________________________
2448 Double_t StEmbeddingQADraw::getVzAcceptedMinimum() const
2449 {
2451 // TH1* hVzAccepted = (TH1D*) getHistogram("hVzAccepted");
2452 //
2453 // Double_t vzMin = -200.0 ;
2454 // for(Int_t ix=0; ix<hVzAccepted->GetNbinsX(); ix++){
2455 // const Double_t count = hVzAccepted->GetBinContent(ix+1);
2456 // if(count!=0){
2457 // vzMin = hVzAccepted->GetBinLowEdge(ix+1);
2458 // LOG_INFO << "Find minimum vz cut, v_z(min) = " << vzMin << " (cm) " << endm ;
2459 // break ;
2460 // }
2461 // }
2462 //
2463 // return vzMin ;
2464 
2465  // Use actual cut value
2466  return -StEmbeddingQAUtilities::instance()->getZVertexCut() ;
2467 }
2468 
2469 //____________________________________________________________________________________________________
2470 Double_t StEmbeddingQADraw::getVzAcceptedMaximum() const
2471 {
2473 // TH1* hVzAccepted = (TH1D*) getHistogram("hVzAccepted");
2474 //
2475 // Double_t vzMax = 200.0 ;
2476 // for(Int_t ix=hVzAccepted->GetNbinsX()-1; ix!=-1; ix--){
2477 // const Double_t count = hVzAccepted->GetBinContent(ix+1);
2478 // if(count!=0){
2479 // vzMax = hVzAccepted->GetBinLowEdge(ix+1) + hVzAccepted->GetBinWidth(ix+1);
2480 // LOG_INFO << "Find maximum vz cut, v_z(max) = " << vzMax << " (cm) " << endm ;
2481 // break ;
2482 // }
2483 // }
2484 //
2485 // return vzMax ;
2486 
2487  // Use actual cut value
2488  return StEmbeddingQAUtilities::instance()->getZVertexCut() ;
2489 }
2490 
2491 //____________________________________________________________________________________________________
2492 Bool_t StEmbeddingQADraw::draw() const
2493 {
2495 
2496  //====================================================================================================
2497  //
2498  // Start drawing
2499  //
2500  //====================================================================================================
2501 
2502  //----------------------------------------------------------------------------------------------------
2504  const Bool_t isEventOk = drawEvent();
2505  //----------------------------------------------------------------------------------------------------
2506 
2507  //----------------------------------------------------------------------------------------------------
2509  const Bool_t isMcTrackOk = drawMcTrack();
2510  //----------------------------------------------------------------------------------------------------
2511 
2512  //----------------------------------------------------------------------------------------------------
2514  const Bool_t isTrackOk = drawTrack();
2515  //----------------------------------------------------------------------------------------------------
2516 
2517  return isEventOk && isMcTrackOk && isTrackOk ;
2518 }
2519 
2520 //____________________________________________________________________________________________________
2521 Bool_t StEmbeddingQADraw::finish()
2522 {
2523  // Close PDF file
2524  LOG_INFO << "StEmbeddingQADraw::finish() Close PDF file : " << mPDF->GetName() << endm ;
2525 
2526  mMainPad->cd();
2527  mMainPad->Clear();
2528 
2529  TLatex* title = new TLatex(0.1, 0.6, "End of QA");
2530  title->SetTextSize(0.1);
2531  title->Draw();
2532 
2533  mCanvas->cd();
2534  mCanvas->Update();
2535 
2536  mPDF->Close() ;
2537 
2538  return kTRUE ;
2539 }
2540 
2541 void StEmbeddingQADraw::setPNGOn() { LOG_INFO << endm << "Print figures to PNG file" << endm << endm; mIsPNGOn = kTRUE ; }
2542 void StEmbeddingQADraw::setGIFOn() { LOG_INFO << endm << "Print figures to GIF file" << endm << endm; mIsGIFOn = kTRUE ; }
2543 void StEmbeddingQADraw::setJPGOn() { LOG_INFO << endm << "Print figures to JPG file" << endm << endm; mIsJPGOn = kTRUE ; }
2544 void StEmbeddingQADraw::setEPSOn() { LOG_INFO << endm << "Print figures to EPS file" << endm << endm; mIsEPSOn = kTRUE ; }
2545 void StEmbeddingQADraw::setPSOn() { LOG_INFO << endm << "Print figures to PS file" << endm << endm; mIsPSOn = kTRUE ; }
2546 
2547 
2548 
StEmbeddingQADraw::drawPhi
Bool_t drawPhi() const
Geant id.
Definition: StEmbeddingQADraw.cxx:1519
Error
Definition: FJcore.h:367
StEmbeddingQADraw::getParticleName
const Char_t * getParticleName(const Int_t geantid=-1) const
Get geant id.
Definition: StEmbeddingQADraw.cxx:923
StEmbeddingQADraw::drawPt
Bool_t drawPt() const
(pseudo-)rapidity in different eta bins
Definition: StEmbeddingQADraw.cxx:1600
StEmbeddingQADraw::drawDca
Bool_t drawDca() const
dE/dx (2D) and projections
Definition: StEmbeddingQADraw.cxx:2211
StEmbeddingQADraw::drawdEdx
Bool_t drawdEdx() const
momentum (|eta|&lt;2, 0.5 eta increment)
Definition: StEmbeddingQADraw.cxx:1903
StEmbeddingQADraw::drawNHit
Bool_t drawNHit() const
Dca vs (pt, eta)
Definition: StEmbeddingQADraw.cxx:2224
StEmbeddingQADraw::drawGeantId
Bool_t drawGeantId() const
Draw Reconstructed track histograms.
Definition: StEmbeddingQADraw.cxx:1437
StEmbeddingQADraw::setJPGOn
void setJPGOn()
Set true for gif file flag (default is false)
Definition: StEmbeddingQADraw.cxx:2543
StEmbeddingQADraw::setPNGOn
void setPNGOn()
Finish QA.
Definition: StEmbeddingQADraw.cxx:2541
StEmbeddingQAUtilities::getCategoryName
TString getCategoryName(const UInt_t id) const
Category from category id.
Definition: StEmbeddingQAUtilities.cxx:140
StEmbeddingQAUtilities::setStyle
void setStyle() const
Check the input geantid is gamma.
Definition: StEmbeddingQAUtilities.cxx:338
StEmbeddingQADraw::setEPSOn
void setEPSOn()
Set true for jpg file flag (default is false)
Definition: StEmbeddingQADraw.cxx:2544
StEmbeddingQAUtilities::isContaminated
Bool_t isContaminated(const TString name) const
Check whether the track is Ghost pair or not.
Definition: StEmbeddingQAUtilities.cxx:212
StEmbeddingQAUtilities::getParticleDefinition
StParticleDefinition * getParticleDefinition(const UInt_t geantid) const
runnumber = runid - (year - 2000 + 1) * 10^6
Definition: StEmbeddingQAUtilities.cxx:439
StEmbeddingQADraw::drawRapidity
Bool_t drawRapidity() const
Azimuthal angle (phi) distributions.
Definition: StEmbeddingQADraw.cxx:1536
StEmbeddingQADraw::drawTrack
Bool_t drawTrack() const
Draw MC track histograms.
Definition: StEmbeddingQADraw.cxx:1396
StEmbeddingQAUtilities::getPtMinCut
Float_t getPtMinCut() const
Get track and event selections.
Definition: StEmbeddingQAUtilities.cxx:469
StEmbeddingQAUtilities::getRunId
Int_t getRunId(const Int_t runnumber, const Int_t year) const
Set font, title and label styles.
Definition: StEmbeddingQAUtilities.h:223
StEmbeddingQADraw::finish
Bool_t finish()
NHit vs (pt, eta)
Definition: StEmbeddingQADraw.cxx:2521
StEmbeddingQAUtilities::instance
static StEmbeddingQAUtilities * instance()
Get instance.
Definition: StEmbeddingQAUtilities.cxx:115
StEmbeddingQADraw::setParentGeantId
void setParentGeantId(const Int_t parentgeantid)
Set parent geant id (default is 0)
Definition: StEmbeddingQADraw.cxx:209
StEmbeddingQADraw::setPSOn
void setPSOn()
Set true for eps file flag (default is false)
Definition: StEmbeddingQADraw.cxx:2545
StEmbeddingQADraw::drawMomentum
Bool_t drawMomentum() const
pt (|eta|&lt;2, 0.5 eta increment)
Definition: StEmbeddingQADraw.cxx:1553
StEmbeddingQADraw::setGIFOn
void setGIFOn()
Set true for png file flag (default is false)
Definition: StEmbeddingQADraw.cxx:2542
StEmbeddingQADraw::drawMcTrack
Bool_t drawMcTrack() const
Draw all event-wise histograms.
Definition: StEmbeddingQADraw.cxx:1272
StEmbeddingQAUtilities::getCategoryTitle
TString getCategoryTitle(const UInt_t id) const
Category name from category id.
Definition: StEmbeddingQAUtilities.cxx:152
StEmbeddingQADraw::drawEvent
Bool_t drawEvent() const
Draw all histograms (event and track histograms).
Definition: StEmbeddingQADraw.cxx:1021
StEmbeddingQAUtilities::getCategoryId
Int_t getCategoryId(const TString name) const
Category title from category id.
Definition: StEmbeddingQAUtilities.cxx:164
StEmbeddingQADraw::init
void init()
Initialization.
Definition: StEmbeddingQADraw.cxx:259
StEmbeddingQADraw::setOutputDirectory
void setOutputDirectory(const TString name="./")
Default is current directory.
Definition: StEmbeddingQADraw.cxx:434
StEmbeddingQADraw::StEmbeddingQADraw
StEmbeddingQADraw(const TString embeddingFile, const TString realDataFile, const Int_t geantid, const Bool_t isEmbeddingOnly=kFALSE)
Definition: StEmbeddingQADraw.cxx:132
Generated by   doxygen 1.8.5