00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087 #include <assert.h>
00088 #include <string>
00089
00090 #include "TCanvas.h"
00091 #include "TError.h"
00092 #include "TF1.h"
00093 #include "TFile.h"
00094 #include "TGraphErrors.h"
00095 #include "TH1.h"
00096 #include "TH2.h"
00097 #include "TH3.h"
00098 #include "TLatex.h"
00099 #include "TLegend.h"
00100 #include "TMath.h"
00101 #include "TObject.h"
00102 #include "TPaveText.h"
00103 #include "TPDF.h"
00104 #include "TProfile.h"
00105 #include "TROOT.h"
00106 #include "TStyle.h"
00107 #include "TSystem.h"
00108
00109 #include "StEmbeddingQADraw.h"
00110 #include "StEmbeddingQAUtilities.h"
00111 #include "StMessMgr.h"
00112 #include "StParticleDefinition.hh"
00113
00114 using namespace std ;
00115
00116 ClassImp(StEmbeddingQADraw)
00117
00118
00119 UInt_t StEmbeddingQADraw::mCanvasId = 0;
00120
00121
00122 StEmbeddingQADraw::StEmbeddingQADraw(const TString embeddingFile, const TString realDataFile, const Int_t geantid,
00123 const Bool_t isEmbeddingOnly)
00124 : mIsEmbeddingOnly(isEmbeddingOnly), mIsGIFOn(kFALSE), mIsJPGOn(kFALSE), mIsEPSOn(kFALSE), mIsPSOn(kFALSE), mGeantId(geantid)
00125 {
00127 mCanvas = 0 ;
00128 mMainPad = 0 ;
00129 mPDF = 0 ;
00130
00131 mInputEmbedding = 0 ;
00132 mInputRealData = 0 ;
00133
00135 mIsOpen = open(embeddingFile, realDataFile);
00136
00138 if(mIsOpen){
00139 TString fileName(embeddingFile);
00140 fileName.Remove(fileName.Index(".root"), fileName.Length());
00141
00142 for(Int_t i=0;i<3;i++){
00143 const Int_t start = 0 ;
00144 const Int_t stop = fileName.First("_") ;
00145 const TString subString(fileName(start, stop));
00146
00147 fileName.Remove(start, stop+1);
00148
00149 if( i == 2 ){
00150 mYear = subString.Atoi();
00151 mProduction = fileName ;
00152 }
00153 }
00154 }
00155 else{
00157 mYear = -10 ;
00158 mProduction = "";
00159 }
00160
00161 mInputParentGeantId = 0 ;
00162 }
00163
00164
00165 StEmbeddingQADraw::StEmbeddingQADraw(const TString embeddingFile, const TString realDataFile,
00166 const Int_t year, const TString production, const Int_t geantid, const Bool_t isEmbeddingOnly)
00167 : mIsEmbeddingOnly(isEmbeddingOnly), mIsGIFOn(kFALSE), mIsJPGOn(kFALSE), mIsEPSOn(kFALSE), mIsPSOn(kFALSE), mGeantId(geantid)
00168 {
00171
00173 mIsOpen = open(embeddingFile, realDataFile);
00174
00176 mYear = year ;
00177 mProduction = production ;
00178
00179 mInputParentGeantId = 0 ;
00180 }
00181
00182
00183 StEmbeddingQADraw::~StEmbeddingQADraw()
00184 {
00186
00188 mMcGeantId.clear();
00189 mDaughterGeantId.clear();
00190 mParentGeantId.clear();
00191 mParentParentGeantId.clear();
00192
00194 if( mInputEmbedding || mInputEmbedding->IsOpen() ) mInputEmbedding->Close();
00195 if( mInputRealData || mInputRealData->IsOpen() ) mInputRealData->Close();
00196 }
00197
00198
00199 void StEmbeddingQADraw::setParentGeantId(const Int_t parentgeantid)
00200 {
00201 mInputParentGeantId = parentgeantid ;
00202 LOG_INFO << "StEmbeddingQADraw::setParentGeantId Set parent geantid = "
00203 << mInputParentGeantId
00204 << endm;
00205 }
00206
00207
00208 void StEmbeddingQADraw::setPtMax(const Double_t ptmax)
00209 {
00210 mPtMax = ptmax ;
00211 LOG_INFO << "StEmbeddingQADraw::setPtMax() Set maximum pt = " << mPtMax
00212 << " GeV/c to be drawn" << endm;
00213 }
00214
00215
00216
00217 Bool_t StEmbeddingQADraw::open(const TString embeddingFile, const TString realDataFile)
00218 {
00220
00221
00222 mInputEmbedding = TFile::Open(embeddingFile);
00223
00224 if(!mInputEmbedding || !mInputEmbedding->IsOpen() || mInputEmbedding->IsZombie() ){
00225 Error("StEmbeddingQADraw", "can't find %s", embeddingFile.Data());
00226 return kFALSE;
00227 }
00228 LOG_INFO << "OPEN input (embedding) : " << mInputEmbedding->GetName() << endm;
00229 LOG_INFO << endm;
00230
00232 if( mIsEmbeddingOnly ) return kTRUE ;
00233
00234
00235 mInputRealData = TFile::Open(realDataFile);
00236
00237 if(!mInputRealData || !mInputRealData->IsOpen() || mInputRealData->IsZombie() ){
00238 Error("StEmbeddingQADraw", "can't find %s", realDataFile.Data());
00239 return kFALSE;
00240 }
00241 LOG_INFO << "OPEN input (real data) : " << mInputRealData->GetName() << endm;
00242 LOG_INFO << endm;
00243
00244 return kTRUE ;
00245 }
00246
00247
00248
00249 void StEmbeddingQADraw::init()
00250 {
00252
00254 mPtMax = 5.0;
00255
00257 mMcGeantId.clear();
00258 mDaughterGeantId.clear();
00259 mParentGeantId.clear();
00260 mParentParentGeantId.clear();
00261
00262 LOG_INFO << "#------------------------------------------------------------" << endm;
00263 LOG_INFO << Form(" Year = %10d", mYear) << endm;
00264 LOG_INFO << Form(" Production = %10s", mProduction.Data()) << endm;
00265 LOG_INFO << Form(" Particle = %10s", getParticleName()) << endm;
00266 LOG_INFO << "#------------------------------------------------------------" << endm;
00267
00269 StEmbeddingQAUtilities::instance()->setStyle();
00270
00272 checkInputGeantId() ;
00273
00275 setDaughterGeantId();
00276
00278 LOG_INFO << "Initialize canvas" << endm ;
00279 mCanvas = new TCanvas(Form("c%d", mCanvasId), Form("c%d", mCanvasId++), 700, 800);
00280 mCanvas->SetFillColor(1);
00281 mCanvas->cd();
00282
00284 mMainPad = new TPad("mainPad", "", 0.00, 0.00, 1.00, 0.90);
00285 mMainPad->SetFillColor(10);
00286 mMainPad->Draw();
00287
00289 mPDF = new TPDF(Form("%sqa_%s.pdf", mOutputFigureDirectory.Data(), getBaseName()));
00290 LOG_INFO << "Initialize PDF file : " << mPDF->GetName() << endm ;
00291
00293 mCanvas->cd();
00294 TPaveText* info = new TPaveText(0.00, 0.95, 0.63, 1.00);
00295 info->SetBorderSize(1);
00296 info->SetFillColor(kBlack);
00297 info->SetTextFont(42);
00298 info->SetTextColor(10);
00299 info->SetTextSize(0.030);
00300
00301 const Int_t nevents = getEntries() ;
00302 TString neventsName(Form("%d", getEntries())) ;
00303 if( nevents > 1000000 ) neventsName = Form("%1.2f M", (Double_t)nevents/1.0e+06);
00304 else if ( nevents > 1000 ) neventsName = Form("%1.2f K", (Double_t)nevents/1.0e+03);
00305
00306 info->AddText(Form("N_{events}=%s, MC=%s, %s, (%d)", neventsName.Data(), getParticleName(), mProduction.Data(), mYear));
00307 info->Draw();
00308
00310 mMainPad->cd();
00311
00312 TLatex* title = new TLatex(0.1, 0.6, "Embedding QA");
00313 title->SetTextSize(0.10);
00314 title->Draw();
00315
00316 TString qaName("embedding + real");
00317 if ( mIsEmbeddingOnly ){
00318 qaName = "embedding only";
00319 }
00320
00321 TLatex* qaTitle = new TLatex(0.12, 0.5, qaName);
00322 qaTitle->SetTextSize(0.07);
00323 qaTitle->Draw();
00324
00327 TLine* lineEmbedding = new TLine(0, 0, 1, 1);
00328 TLine* lineRealBlue = new TLine(0, 0, 1, 1);
00329 TLine* lineRealBlack = new TLine(0, 0, 1, 1);
00330 lineEmbedding->SetLineColor(kRed); lineEmbedding->SetLineWidth(2);
00331 lineRealBlue ->SetLineColor(kBlue); lineRealBlue ->SetLineWidth(2);
00332 lineRealBlack->SetLineColor(kBlack); lineRealBlack->SetLineWidth(2);
00333
00334 TLegend* colorCode = new TLegend(0.05, 0.15, 0.8, 0.45);
00335 colorCode->SetTextSize(0.035);
00336 colorCode->SetFillColor(10);
00337 colorCode->SetHeader("Color code for embedding and real data");
00338 colorCode->AddEntry(lineRealBlack, "MC (black)", "L");
00339 colorCode->AddEntry(lineEmbedding, "Reconstructed embedding tracks* (red)", "L");
00340 colorCode->AddEntry(lineRealBlue, "Real** (blue)", "L");
00341 colorCode->Draw();
00342
00343 TLatex* note0 = new TLatex(0.1, 0.10, "* matched pairs or contaminated pairs");
00344 TLatex* note1 = new TLatex(0.1, 0.05, "** black is also used, see legend for each pad");
00345 note0->SetTextSize(0.03);
00346 note1->SetTextSize(0.03);
00347 note0->Draw();
00348 note1->Draw();
00349
00350 mCanvas->cd();
00351 mCanvas->Update();
00352
00353 mPDF->NewPage() ;
00354
00355
00356 const StEmbeddingQAUtilities* utility = StEmbeddingQAUtilities::instance() ;
00357 TPaveText* header = initCanvas("Event & track selections");
00358
00359
00360 TPaveText* eventSelections = new TPaveText(0.05, 0.60, 0.95, 0.92);
00361 eventSelections->SetTextAlign(12);
00362 eventSelections->SetBorderSize(1);
00363 eventSelections->SetFillColor(10);
00364 eventSelections->SetTextColor(kBlack);
00365 eventSelections->SetTextSize(0.030);
00366 eventSelections->SetTextFont(42);
00367
00368 eventSelections->AddText("*** Event selections");
00369
00370 eventSelections->AddText(Form(" z-vertex cut : |v_{z}| < %1.1f cm", utility->getZVertexCut()));
00371
00372 const vector<UInt_t> triggerId(utility->getTriggerIdCut());
00373 if ( !triggerId.empty() ) {
00374 TString triggers("");
00375 for(UInt_t i=0; i<triggerId.size(); i++) {
00376 triggers += Form("%d", triggerId[i]);
00377 if( i != triggerId.size() - 1 ) triggers += ", ";
00378 }
00379 eventSelections->AddText(Form(" trigger id cut : id = %s", triggers.Data()));
00380 }
00381 eventSelections->AddText("NOTE: Trigger id cut for real data has to be made manually in doEmbeddingQAMaker.C");
00382 eventSelections->SetAllWith("***", "color", kRed);
00383 eventSelections->SetAllWith("***", "font", 72);
00384 eventSelections->SetAllWith("***", "size", 0.033);
00385 eventSelections->SetAllWith("NOTE", "color", kBlue);
00386 eventSelections->SetAllWith("NOTE", "font", 72);
00387 eventSelections->SetAllWith("NOTE", "size", 0.020);
00388 eventSelections->Draw();
00389
00390
00391 TPaveText* trackSelections = new TPaveText(0.05, 0.05, 0.95, 0.59);
00392 trackSelections->SetTextAlign(12);
00393 trackSelections->SetBorderSize(1);
00394 trackSelections->SetFillColor(10);
00395 trackSelections->SetTextColor(kBlack);
00396 trackSelections->SetTextSize(0.030);
00397 trackSelections->SetTextFont(42);
00398 trackSelections->AddText("*** Track selections");
00399 trackSelections->AddText(Form(" %1.1f < p_{T} < %1.1f GeV/c", utility->getPtMinCut(), utility->getPtMaxCut()));
00400 trackSelections->AddText(Form(" |#eta| < %1.2f", utility->getEtaCut()));
00401 trackSelections->AddText(Form(" |y| < %1.2f", utility->getRapidityCut()));
00402 trackSelections->AddText(Form(" nHitsFit > %3d", utility->getNHitCut()));
00403 trackSelections->AddText(Form(" nHitsFit/nHitsPoss > %1.2f", utility->getNHitToNPossCut()));
00404 trackSelections->AddText(Form(" global Dca < %1.1f cm", utility->getDcaCut()));
00405 trackSelections->AddText(Form(" |n#sigma| < %1.1f", utility->getNSigmaCut()));
00406 trackSelections->AddText("NOTE1: Rapidity cut for real data has to be made manually in doEmbeddingQAMaker.C");
00407 trackSelections->AddText("NOTE2: Cut on its own variable is currently disabled, e.x. no dca cut for dca histogram");
00408 trackSelections->SetAllWith("***", "color", kRed);
00409 trackSelections->SetAllWith("***", "font", 72);
00410 trackSelections->SetAllWith("***", "size", 0.033);
00411 trackSelections->SetAllWith("NOTE", "color", kBlue);
00412 trackSelections->SetAllWith("NOTE", "font", 72);
00413 trackSelections->SetAllWith("NOTE", "size", 0.020);
00414 trackSelections->Draw();
00415
00416 mCanvas->cd();
00417 mCanvas->Update();
00418 mPDF->NewPage() ;
00419
00420 delete header ;
00421 }
00422
00423
00424 void StEmbeddingQADraw::setOutputDirectory(const TString name)
00425 {
00427
00428 mOutputFigureDirectory = name ;
00429
00431 if( gSystem->AccessPathName(name) == kTRUE ){
00432 Error("setOutputDirectory", "Directory %s does not exist. Set current directory as the output location");
00433 mOutputFigureDirectory = "./";
00434 }
00435
00437 if ( mOutputFigureDirectory.Last('/') + 1 != mOutputFigureDirectory.Length() ){
00438 LOG_INFO << endm;
00439 LOG_INFO << "Put / at the end of output directory name" << endm;
00440 mOutputFigureDirectory.Append("/");
00441 }
00442
00443 LOG_INFO << "Set output directory for figures : " << mOutputFigureDirectory << endm;
00444 }
00445
00446
00447 void StEmbeddingQADraw::print(const TCanvas& canvas, const TString name) const
00448 {
00450
00451 if( mIsPNGOn ) canvas.Print(Form("%s%s_%s.png", mOutputFigureDirectory.Data(), name.Data(), getBaseName()));
00452 if( mIsGIFOn ) canvas.Print(Form("%s%s_%s.gif", mOutputFigureDirectory.Data(), name.Data(), getBaseName()));
00453 if( mIsJPGOn ) canvas.Print(Form("%s%s_%s.jpg", mOutputFigureDirectory.Data(), name.Data(), getBaseName()));
00454 if( mIsEPSOn ) canvas.Print(Form("%s%s_%s.eps", mOutputFigureDirectory.Data(), name.Data(), getBaseName()));
00455 if( mIsPSOn ) canvas.Print(Form("%s%s_%s.ps", mOutputFigureDirectory.Data(), name.Data(), getBaseName()));
00456 }
00457
00458
00459 void StEmbeddingQADraw::checkInputGeantId()
00460 {
00463
00464 TH1* hGeantId = (TH1D*) getHistogram("hGeantId_0");
00465
00467 if(!hGeantId){
00468 mGeantId = -1 ;
00469 return ;
00470 }
00471
00473 for(Int_t id=0; id<hGeantId->GetNbinsX(); id++){
00474 const Bool_t isGeantIdOk = hGeantId->GetBinContent(id+1)>0.0;
00475 if(!isGeantIdOk) continue ;
00476
00477 const Int_t geantid = TMath::Nint(hGeantId->GetBinLowEdge(id+1));
00478 mMcGeantId.push_back(geantid);
00479 }
00480
00481 const StEmbeddingQAUtilities* utility = StEmbeddingQAUtilities::instance() ;
00482
00483 if( mMcGeantId.size() == 1 ){
00488 if ( mGeantId == mMcGeantId[0] ) return ;
00489 else{
00490 LOG_INFO << "#----------------------------------------------------------------------------------------------------" << endm;
00491 LOG_INFO << " Input geant id doesn't correspond to the geant id in the MC histogram" << endm ;
00492 LOG_INFO << " Geantid in the MC histogram, geantid = " << mMcGeantId[0]
00493 << ", particle name = " << utility->getParticleDefinition(mMcGeantId[0])->name().c_str() << endm;
00494 LOG_INFO << " Do you want to proceed to the QA for geantid = " << mMcGeantId[0] << " ?" << endm;
00495 LOG_INFO << " [yes=1/no=0]:" << endm;
00496 UInt_t proceedQA = 0 ;
00497 cin >> proceedQA ;
00498 LOG_INFO << "#----------------------------------------------------------------------------------------------------" << endm;
00499
00501 assert(proceedQA);
00502
00504 mGeantId = mMcGeantId[0] ;
00505 }
00506 }
00507 else{
00512 for(vector<Int_t>::iterator iter = mMcGeantId.begin(); iter != mMcGeantId.end(); iter++){
00513 const Int_t geantidFound = (*iter) ;
00514 if ( mGeantId == geantidFound ) return ;
00515 }
00516
00518 LOG_INFO << "#----------------------------------------------------------------------------------------------------" << endm;
00519 LOG_INFO << " Cannot find input geantid in the MC histogram" << endm;
00520 LOG_INFO << " Available geantid's are listed below;" << endm;
00521 LOG_INFO << " geantid particle name" << endm;
00522 for(vector<Int_t>::iterator iter = mMcGeantId.begin(); iter != mMcGeantId.end(); iter++){
00523 const Int_t geantidFound = (*iter) ;
00524 LOG_INFO << Form(" %5d %10s", geantidFound,
00525 utility->getParticleDefinition(geantidFound)->name().c_str()) << endm;
00526 }
00527 LOG_INFO << " Which geantid you want to do the QA ?" << endm;
00528 LOG_INFO << " [Put the one geantid listed above. If you want to stop the program, please put 0]:" << endm;
00529 cin >> mGeantId ;
00530
00532 assert(mGeantId);
00533
00534 LOG_INFO << " QA for geantid = " << mGeantId
00535 << ", particle name = " << utility->getParticleDefinition(mGeantId)->name().c_str() << endm;
00536 LOG_INFO << "#----------------------------------------------------------------------------------------------------" << endm;
00537 }
00538
00539 }
00540
00541
00542 Bool_t StEmbeddingQADraw::isOpen() const
00543 {
00549
00550 if(!mIsOpen) LOG_INFO << "No input files found" << endm;
00551 assert(mIsOpen);
00552
00553 const Bool_t isGeantIdOk = mGeantId>0 ;
00554 if(!isGeantIdOk) LOG_INFO << "Cannot find input geantid in the histogram or no histogram in the input" << endm;
00555 assert(isGeantIdOk);
00556
00557 return kTRUE ;
00558 }
00559
00560
00561 Bool_t StEmbeddingQADraw::isMatchedPairOk() const
00562 {
00564
00565
00566
00567
00568
00569
00570
00572
00573
00574 if ( StEmbeddingQAUtilities::instance()->isEPiKP(mGeantId) ) {
00575
00576 return kTRUE ;
00577 }
00578 else if ( StEmbeddingQAUtilities::instance()->isGamma(mGeantId) ) {
00579
00580 return kFALSE ;
00581 }
00582 else{
00583 return StEmbeddingQAUtilities::instance()->getParticleDefinition(mGeantId)->stable() ;
00584 }
00585 }
00586
00587
00588
00589 void StEmbeddingQADraw::setDaughterGeantId()
00590 {
00592
00594 if(isMatchedPairOk()){
00596 mDaughterGeantId.push_back(mGeantId);
00597
00598 mParentGeantId.push_back(0);
00599 mParentParentGeantId.push_back(0);
00600 return;
00601 }
00602
00604 for(UInt_t id=0; id<1000; id++){
00605 const StEmbeddingQAUtilities* utility = StEmbeddingQAUtilities::instance() ;
00606 const Int_t contamCategoryId = utility->getCategoryId("CONTAM") ;
00607
00608 Int_t n = 0 ;
00609 if( mInputParentGeantId == 0 ) n = 1 ;
00610 else n = 2 ;
00611
00612 for(Int_t i=0; i<n; i++) {
00613 Int_t daughterId = id ;
00614 Int_t parentId = mGeantId ;
00615 Int_t parentParentId = 0 ;
00616 if( i == 1 ) {
00617 daughterId = id ;
00618 parentId = mInputParentGeantId ;
00619 parentParentId = mGeantId ;
00620 }
00621
00622
00623 TH3* hDca = (TH3D*) mInputEmbedding->Get(Form("hDca_%d_%d_%d_%d", contamCategoryId, parentParentId, parentId, daughterId)) ;
00624
00625 if ( hDca ){
00626 if( i == 0 ){
00627 const Char_t* daughterName = utility->getParticleDefinition(daughterId)->name().c_str() ;
00628 const Char_t* parentName = utility->getParticleDefinition(parentId)->name().c_str() ;
00629
00630 LOG_INFO << Form("Find daughter %10s from parent %10s", daughterName, parentName) << endm;
00631 }
00632 else{
00633 const Char_t* daughterName = utility->getParticleDefinition(daughterId)->name().c_str() ;
00634 const Char_t* parentName = utility->getParticleDefinition(parentId)->name().c_str() ;
00635 const Char_t* parentparentName = utility->getParticleDefinition(parentParentId)->name().c_str() ;
00636
00637 LOG_INFO << Form("Find daughter %10s from parent %10s (from %10s)", daughterName, parentName, parentparentName) << endm;
00638 }
00639
00640 mDaughterGeantId.push_back(id);
00641 mParentGeantId.push_back(parentId);
00642 mParentParentGeantId.push_back(parentParentId);
00643 }
00644 }
00645 }
00646
00648 if ( mDaughterGeantId.empty() ) mDaughterGeantId.push_back(mGeantId);
00649 }
00650
00651
00652 Int_t StEmbeddingQADraw::getEntries() const
00653 {
00655
00656 TH1* hVz = (TH1D*) getHistogram("hVzAccepted");
00657 if(!hVz) return 0;
00658
00659 return (Int_t)( hVz->GetEntries() );
00660 }
00661
00662
00663 Bool_t StEmbeddingQADraw::isDecay() const
00664 {
00666
00669 if(isMatchedPairOk()) return kFALSE ;
00670
00672 return ( mDaughterGeantId.size() > 1 ) ;
00673 }
00674
00675
00676 Int_t StEmbeddingQADraw::getGeantIdReal(const Int_t daughter) const
00677 {
00680
00681 return ( StEmbeddingQAUtilities::instance()->isEPiKP(mDaughterGeantId[daughter]) ) ? mDaughterGeantId[daughter] : 8 ;
00682 }
00683
00684
00685 Int_t StEmbeddingQADraw::getCategoryId(const Bool_t isEmbedding) const
00686 {
00692 const StEmbeddingQAUtilities* utility = StEmbeddingQAUtilities::instance() ;
00693
00694 if ( isEmbedding ){
00695 return (isDecay()) ? utility->getCategoryId("CONTAM") : utility->getCategoryId("MATCHED") ;
00696 }
00697 else{
00698 return utility->getCategoryId("PRIMARY") ;
00699 }
00700 }
00701
00702
00703 TObject* StEmbeddingQADraw::getHistogram(const TString name, const Bool_t isEmbedding) const
00704 {
00706
00707 if(!isOpen()) return 0;
00708
00711 if( mIsEmbeddingOnly && !isEmbedding) return 0 ;
00712
00713 TObject* obj = (isEmbedding) ? mInputEmbedding->Get(name) : mInputRealData->Get(name) ;
00714 if ( !obj ){
00715 Error("StEmbeddingQADraw::getHistogram", "Histogram %s doesn't exist", name.Data());
00716 return 0;
00717 }
00718
00719 LOG_DEBUG << "StEmbeddingQADraw::getHistogram() get histogram = " << name << endm;
00720
00721 return obj ;
00722 }
00723
00724
00725
00726 TObject* StEmbeddingQADraw::getHistogram(const TString name, const UInt_t daughter, const Bool_t isEmbedding,
00727 const UInt_t parentparentid) const
00728 {
00733
00734 if ( daughter >= mDaughterGeantId.size() ){
00735 Error("StEmbeddingQADraw::getHistogram", "Unknown daughter index, index=%3d", daughter);
00736 return 0;
00737 }
00738
00739 return getHistogram(getHistogramName(name, daughter, isEmbedding, parentparentid), isEmbedding) ;
00740 }
00741
00742
00743 const Char_t* StEmbeddingQADraw::getHistogramName(const TString name, const UInt_t daughter, const Bool_t isEmbedding,
00744 const UInt_t parentparentid) const
00745 {
00746 const Int_t category = getCategoryId(isEmbedding) ;
00747
00748 if( isEmbedding ){
00752
00753 if( isDecay() ) {
00754
00755
00756 return Form("%s_%d_%d_%d_%d", name.Data(), category, mParentParentGeantId[daughter], mParentGeantId[daughter], mDaughterGeantId[daughter]) ;
00757 }
00758 else {
00759 return Form("%s_%d_%d", name.Data(), category, mGeantId) ;
00760 }
00761 }
00762 else{
00766 const Int_t geantidReal = getGeantIdReal(daughter) ;
00767
00768 return Form("%s_%d_%d", name.Data(), category, geantidReal) ;
00769 }
00770 }
00771
00772
00773 Double_t StEmbeddingQADraw::getNormalization(const TH1& h) const
00774 {
00777
00778 const Double_t ntrack = h.Integral() ;
00779 if( ntrack == 0 ) return 1.0 ;
00780
00781 const Double_t binWidth = h.GetBinWidth(1) ;
00782 if( binWidth == 0.0 ) return 1.0 ;
00783
00784 return 1.0/(ntrack*binWidth);
00785 }
00786
00787
00788 const Char_t* StEmbeddingQADraw::getBaseName() const
00789 {
00791
00792 return Form("%d_%s_%s", mYear, mProduction.Data(), getParticleName());
00793 }
00794
00795
00796 Bool_t StEmbeddingQADraw::drawStatistics(const Double_t x1, const Double_t y1, const Double_t x2, const Double_t y2,
00797 const Double_t textSize) const
00798 {
00804
00805 TPaveText* statistics = new TPaveText(x1, y1, x2, y2);
00806 statistics->SetTextFont(42);
00807 statistics->SetTextSize(textSize);
00808 statistics->SetBorderSize(1);
00809 statistics->SetFillColor(10);
00810 statistics->AddText(Form("N_{evts} = %d", getEntries()));
00811 statistics->AddText(Form("Year: %d", mYear));
00812 statistics->AddText(Form("Production: %s", mProduction.Data()));
00813 statistics->AddText(Form("Particle: %s", getParticleName()));
00814 statistics->Draw();
00815
00816 return kTRUE ;
00817 }
00818
00819
00820 TPaveText* StEmbeddingQADraw::drawHeader(const TString description,
00821 const Double_t x1, const Double_t y1, const Double_t x2, const Double_t y2, const Double_t textSize) const
00822 {
00824 TPaveText* header = new TPaveText(x1, y1, x2, y2);
00825 header->SetBorderSize(1);
00826 header->SetFillColor(kBlack);
00827 header->SetTextColor(10);
00828 header->SetTextSize(textSize);
00829 header->SetTextFont(42);
00830 header->AddText(description);
00831 header->Draw();
00832
00833 return header ;
00834 }
00835
00836
00837 void StEmbeddingQADraw::drawLegend(const UInt_t id, TH1* hembed, TH1* hreal,
00838 const Option_t* option, const Bool_t doSplit) const
00839 {
00840 TLegend* leg = new TLegend(0, 0.2, 1, 0.8);
00841 leg->SetFillColor(10);
00842
00843
00844 leg->SetTextFont(43);
00845 leg->SetTextSize(15);
00846
00847 leg->AddEntry(hembed, getEmbeddingParticleName(id, doSplit), option) ;
00848 if(hreal) leg->AddEntry(hreal, getRealParticleName(id, doSplit), option) ;
00849 leg->Draw();
00850 }
00851
00852
00853 void StEmbeddingQADraw::drawErrorMessages(const TString histogramName) const
00854 {
00855 const StEmbeddingQAUtilities* utility = StEmbeddingQAUtilities::instance() ;
00856 const Int_t categoryId = getCategoryId() ;
00857 const TString title(utility->getCategoryTitle(categoryId));
00858
00859 TPaveText* error0 = new TPaveText(0.05, 0.60, 0.95, 0.92, "NDC");
00860 error0->SetTextAlign(12);
00861 error0->SetBorderSize(1);
00862 error0->SetFillColor(10);
00863 error0->SetTextColor(kBlack);
00864 error0->SetTextSize(0.030);
00865 error0->SetTextFont(42);
00866 error0->AddText(Form("*** No histogram \"%s\" found", histogramName.Data()));
00867 error0->AddText(Form("Please make sure that you have %s in the minimc.root.", title.Data()));
00868 error0->AddText("Open minimc.root file and then check number of tracks.");
00869 error0->AddText(Form("There are maybe no %s in the minimc.", title.Data()));
00870 error0->AddText("Or you used older base QA codes.");
00871 error0->AddText("In case you have finite number of tracks,");
00872 error0->AddText("please also have a look at geantid.");
00873 error0->SetAllWith("***", "color", kRed);
00874 error0->SetAllWith("***", "font", 72);
00875 error0->SetAllWith("***", "size", 0.033);
00876
00877
00878 TPaveText* error1 = new TPaveText(0.05, 0.05, 0.95, 0.55, "NDC");
00879 error1->SetTextAlign(12);
00880 error1->SetBorderSize(1);
00881 error1->SetFillColor(10);
00882 error1->SetTextColor(kBlack);
00883 error1->SetTextSize(0.025);
00884 error1->SetTextFont(42);
00885 error1->AddText("See example below how to check number of tracks and geantid");
00886
00887 if ( isDecay() ) {
00888 error1->AddText("[ROOT]> StMiniMcTree->Draw(\"mNContamPair\")");
00889 error1->AddText("or");
00890 error1->AddText("[ROOT]> StMiniMcTree->Scan(\"mNContamPair\")");
00891 error1->AddText(Form("For geantid in %s", title.Data()));
00892 error1->AddText("[ROOT]> StMiniMcTree->Draw(\"mContamPairs.mGeantId\")");
00893 error1->AddText("or use \"Scan()\" function");
00894 }
00895 else {
00896 error1->AddText("[ROOT]> StMiniMcTree->Draw(\"mNMatchedPair\")");
00897 error1->AddText("or");
00898 error1->AddText("[ROOT]> StMiniMcTree->Scan(\"mNMatchedPair\")");
00899 error1->AddText(Form("For geantid in %s", title.Data()));
00900 error1->AddText("[ROOT]> StMiniMcTree->Draw(\"mMatchedPairs.mGeantId\")");
00901 error1->AddText("or use \"Scan()\" function");
00902 }
00903 error1->SetAllWith("ROOT", "color", kBlue);
00904 error1->SetAllWith("ROOT", "font", 52);
00905 error1->SetAllWith("ROOT", "size", 0.028);
00906
00907 error0->Draw();
00908 error1->Draw();
00909 }
00910
00911
00912 const Char_t* StEmbeddingQADraw::getParticleName(const Int_t geantid) const
00913 {
00915
00917 const Int_t id = (geantid<0) ? mGeantId : geantid ;
00918 const StParticleDefinition* particle = StEmbeddingQAUtilities::instance()->getParticleDefinition(id) ;
00919
00920 if(!particle){
00921 Error("StEmbeddingQADraw::getParticleName", "Can't find particle geantid=%d in StParticleDefinition", id);
00922 assert(0);
00923 }
00924
00930 TString name(particle->name().c_str());
00931 while( name.Contains("/") ) name.Remove(name.Index("/"), 1);
00932 while( name.Contains("(") ) name.Remove(name.Index("("), 1);
00933 while( name.Contains(")") ) name.Remove(name.Index(")"), 1);
00934 while( name.Contains("*") ) name.Replace(name.Index("*"), 1, "star");
00935
00936 return name.Data() ;
00937 }
00938
00939
00940 const Char_t* StEmbeddingQADraw::getMcParticleName() const
00941 {
00942 const StEmbeddingQAUtilities* utility = StEmbeddingQAUtilities::instance() ;
00943 const Int_t categoryId = getCategoryId() ;
00944 const TString name(utility->getCategoryName(categoryId));
00945 const TString parent( (utility->isContaminated(name)) ? "Parent " : "" );
00946
00947 return Form("%s%s (MC, geantid=%d)", parent.Data(), getParticleName(mGeantId), mGeantId);
00948 }
00949
00950
00951 const Char_t* StEmbeddingQADraw::getEmbeddingParticleName(const UInt_t id, const Bool_t doSplit) const
00952 {
00953 if(id>=mDaughterGeantId.size()){
00954 Error("StEmbeddingQADraw::getEmbeddingParticleName", "Unknown daughter particle index, id=%3d", id);
00955 return "";
00956 }
00957
00958 const StEmbeddingQAUtilities* utility = StEmbeddingQAUtilities::instance() ;
00959 const Int_t categoryId = getCategoryId() ;
00960 const TString name(utility->getCategoryName(categoryId));
00961
00962
00963
00964
00965 const Int_t parentParentId = mParentParentGeantId[id] ;
00966 const Int_t parentId = mParentGeantId[id] ;
00967 const Int_t daughterId = mDaughterGeantId[id] ;
00968
00969 const TString daughter( (utility->isContaminated(name)) ? "Daughter " : "" );
00970 TString parent("");
00971
00972 if ( parentId != 0 ) {
00973 const TString parentName(getParticleName(parentId));
00974
00975 if ( parentParentId == 0 ) {
00976 parent = " (from " + parentName + ")";
00977 }
00978 else{
00979 const TString parentParentName(getParticleName(parentParentId));
00980 parent = " (from " + parentName + " #leftarrow " + parentParentName + ")";
00981 }
00982 }
00983
00984 return (doSplit) ? Form("#splitline{%s%s%s}{(%s, geantid=%d)}", daughter.Data(), getParticleName(daughterId),
00985 parent.Data(), utility->getCategoryName(categoryId).Data(), daughterId)
00986 : Form("%s%s%s (%s, geantid=%d)", daughter.Data(), getParticleName(daughterId),
00987 parent.Data(), utility->getCategoryName(categoryId).Data(), daughterId);
00988 }
00989
00990
00991 const Char_t* StEmbeddingQADraw::getRealParticleName(const UInt_t id, const Bool_t doSplit) const
00992 {
00993
00994
00995 if(id>=mDaughterGeantId.size()){
00996 Error("StEmbeddingQADraw::getRealParticleName", "Unknown daughter particle index, id=%3d", id);
00997 return "";
00998 }
00999
01000 const Int_t geantid = getGeantIdReal(id) ;
01001 return (doSplit) ? Form("#splitline{%s}{(%s, |n #sigma_{%s}|<2)}", getParticleName(geantid),
01002 StEmbeddingQAUtilities::instance()->getCategoryName(getCategoryId(kFALSE)).Data(), getParticleName(geantid))
01003 : Form("%s (%s, |n #sigma_{%s}|<2)", getParticleName(geantid),
01004 StEmbeddingQAUtilities::instance()->getCategoryName(getCategoryId(kFALSE)).Data(), getParticleName(geantid));
01005 }
01006
01007
01008 Bool_t StEmbeddingQADraw::drawEvent() const
01009 {
01011
01012 LOG_INFO << endm;
01013 LOG_INFO << "----------------------------------------------------------------------------------------------------" << endm ;
01014 LOG_INFO << "QA for Event-wise informations ..." << endm;
01015
01017 if(!isOpen()) return kFALSE ;
01018
01019
01020
01021
01022
01023
01024
01025 drawVertices() ;
01026
01027
01029 drawRunEventId() ;
01030
01031
01033 drawNumberOfParticles() ;
01034
01035 LOG_INFO << "----------------------------------------------------------------------------------------------------" << endm ;
01036 LOG_INFO << endm;
01037
01038 return kTRUE ;
01039 }
01040
01041
01042 Bool_t StEmbeddingQADraw::drawVertices() const
01043 {
01045 LOG_INFO << "QA for event-vertices ..." << endm;
01046
01048 const Double_t vzMin = getVzAcceptedMinimum() ;
01049 const Double_t vzMax = getVzAcceptedMaximum() ;
01050 TPaveText* header = initCanvas(Form("Event vertices, offline cuts: %1.1f < v_{z} < %1.1f cm", vzMin, vzMax), 2, 2);
01051
01053 mMainPad->cd(1);
01054 TH1* hVz = (TH1D*) getHistogram("hVz");
01055 if(!hVz) return kFALSE ;
01056
01057 TH1* hVzAccepted = (TH1D*) getHistogram("hVzAccepted");
01058 if(!hVzAccepted) return kFALSE ;
01059
01060 hVzAccepted->SetLineColor(kCyan);
01061 hVzAccepted->SetFillColor(kCyan);
01062 hVz->Draw();
01063 hVzAccepted->Draw("same");
01064 hVz->Draw("same");
01065
01067 mMainPad->cd(2);
01068 TH2* hVyVx = (TH2D*) getHistogram("hVyVx");
01069 if(!hVyVx) return kFALSE ;
01070 hVyVx->Draw("colz");
01071
01073 mMainPad->cd(3);
01074 TH2* hVxVz = (TH2D*) getHistogram("hVxVz");
01075 if(!hVxVz) return kFALSE ;
01076 hVxVz->SetAxisRange(vzMin, vzMax, "X");
01077 hVxVz->Draw("colz");
01078
01080 mMainPad->cd(4);
01081 TH2* hVyVz = (TH2D*) getHistogram("hVyVz");
01082 if(!hVyVz) return kFALSE ;
01083 hVyVz->SetAxisRange(vzMin, vzMax, "X");
01084 hVyVz->Draw("colz");
01085
01086 mCanvas->cd();
01087 mCanvas->Update();
01088 mPDF->NewPage() ;
01089
01091 delete header ;
01092 header = initCanvas("Event vertices, #Deltav = v(Data) - v(MC)", 2, 2);
01093
01095 for(Int_t i=0; i<3; i++){
01096 mMainPad->cd(i+1);
01097 TH1* hdV = 0 ;
01098 if( i == 0 ) hdV = (TH1D*) getHistogram("hdVx");
01099 if( i == 1 ) hdV = (TH1D*) getHistogram("hdVy");
01100 if( i == 2 ) hdV = (TH1D*) getHistogram("hdVz");
01101 if(!hdV) return kFALSE ;
01102
01103 hdV->Draw();
01104 }
01105
01106 mCanvas->cd();
01107 mCanvas->Update();
01108 mPDF->NewPage() ;
01109
01110
01111
01112 delete header ;
01113
01114 return kTRUE ;
01115 }
01116
01117
01118 Bool_t StEmbeddingQADraw::drawRunEventId() const
01119 {
01121 LOG_INFO << "QA for run and event id ..." << endm;
01122
01123 TPaveText* header = initCanvas("Run and event id", 2, 2);
01124
01125 mMainPad->cd(1);
01126 TH1* hRunNumber = (TH1D*) mInputEmbedding->Get("hRunNumber");
01127 if(!hRunNumber) return kFALSE ;
01128 hRunNumber->Draw();
01129
01130 mMainPad->cd(2);
01131 TH1* hEventId = (TH1D*) mInputEmbedding->Get("hEventId");
01132 if(!hEventId) return kFALSE ;
01133 hEventId->Draw();
01134
01136 TPaveText* runlist[2];
01137 for(Int_t i=0;i<2;i++){
01138 runlist[i] = new TPaveText(0.1, 0.05, 0.9, 0.95);
01139 runlist[i]->SetTextFont(42);
01140 runlist[i]->SetTextSize(0.04);
01141 runlist[i]->SetBorderSize(1);
01142 runlist[i]->SetFillColor(10);
01143 runlist[i]->AddText("Run id statistics");
01144 }
01145
01146
01147 Int_t runTotal = 0 ;
01148 for(Int_t irun=0; irun<hRunNumber->GetNbinsX(); irun++){
01149 if(hRunNumber->GetBinContent(irun+1)>0.0) runTotal++;
01150 }
01151
01152
01153 Int_t runAccepted = 0 ;
01154 for(Int_t irun=0; irun<hRunNumber->GetNbinsX(); irun++){
01155 const Double_t count = hRunNumber->GetBinContent(irun+1);
01156 if(count==0.0) continue ;
01157
01158 Int_t pad = 0;
01159 if( runTotal < 10 ) pad = 0 ;
01160 else{
01161
01162 if( runAccepted < runTotal/2 ) pad = 0;
01163 else pad = 1 ;
01164 }
01165
01166 const Int_t runid = StEmbeddingQAUtilities::instance()->getRunId((Int_t)hRunNumber->GetBinLowEdge(irun+1), mYear);
01167 runlist[pad]->AddText(Form("%10d %10d events", runid, (Int_t)count));
01168
01169 runAccepted++ ;
01170 }
01171
01172 mMainPad->cd(3) ;
01173 runlist[0]->Draw();
01174
01175 mMainPad->cd(4) ;
01176 runlist[1]->Draw();
01177
01178 mCanvas->cd();
01179 mCanvas->Update();
01180
01181 mPDF->NewPage() ;
01182
01183
01184 delete header ;
01185
01186 return kTRUE ;
01187 }
01188
01189
01190 Bool_t StEmbeddingQADraw::drawNumberOfParticles() const
01191 {
01193 LOG_INFO << "QA for multiplicity distributions (number of particles per event) ..." << endm;
01194
01195 TPaveText* header = initCanvas("Multiplicity distribution", 2, 3);
01196
01197 for(UInt_t ic=0; ic<StEmbeddingQAConst::mNEmbedding; ic++){
01198 mMainPad->cd(ic+1);
01199 mMainPad->GetPad(ic+1)->SetLogy();
01200
01201 TH1* hNParticles = (TH1D*) mInputEmbedding->Get(Form("hNParticles_%d", ic));
01202 if(!hNParticles) return kFALSE ;
01203
01204 TString title(hNParticles->GetTitle());
01205 title.Remove(0, title.Index(',')+2);
01206 hNParticles->SetTitle(title);
01207
01208
01209 if( hNParticles->GetMean() == 0 && hNParticles->GetRMS() == 0 ){
01210
01211 hNParticles->SetAxisRange(0, 10, "X");
01212
01213 LOG_DEBUG << "StEmbeddingQADraw::drawNumberOfParticles no particles for " << title << endm ;
01214 }
01215 else if( hNParticles->GetRMS() == 0 ){
01216
01217 const Double_t mean = hNParticles->GetMean() ;
01218 const Double_t xmin = (mean-10.0<0.0) ? 0.0 : mean-10.0 ;
01219 const Double_t xmax = mean + 10.0 ;
01220 hNParticles->SetAxisRange(xmin, xmax, "X");
01221
01222 LOG_DEBUG << "StEmbeddingQADraw::drawNumberOfParticles constant number of particles = " << mean
01223 << " for " << title
01224 << endm;
01225 }
01226 else{
01227
01228
01229
01230 Double_t xmax = 0.0 ;
01231 for(Int_t i=hNParticles->GetNbinsX()-1; i!=0; i--){
01232 if( hNParticles->GetBinContent(i+1) != 0.0 ){
01233
01234 xmax = hNParticles->GetBinCenter(i+1) ;
01235 break ;
01236 }
01237 }
01238 hNParticles->SetAxisRange(0, xmax+10, "X");
01239
01240 LOG_DEBUG << "StEmbeddingQADraw::drawNumberOfParticles varied multiplicity, max range = " << xmax
01241 << " for " << title
01242 << endm;
01243 }
01244
01245 hNParticles->Draw();
01246 }
01247
01248 mCanvas->cd();
01249 mCanvas->Update();
01250
01251 mPDF->NewPage() ;
01252
01253 delete header ;
01254
01255 return kTRUE ;
01256 }
01257
01258
01259 Bool_t StEmbeddingQADraw::drawMcTrack() const
01260 {
01262
01263 LOG_INFO << "QA for MC tracks (pt, eta, y, phi) ..." << endm;
01264
01266 if(!isOpen()) return kFALSE ;
01267
01268
01269
01270 TPaveText* header = initCanvas("MC track QA (2D)", 2, 2);
01271
01273 mMainPad->cd(1);
01274 TH2* hPtVsEta = (TH2D*) getHistogram(Form("hPtVsEta_0_%d", mGeantId));
01275 if(!hPtVsEta) return kFALSE;
01276
01277 hPtVsEta->SetTitle("");
01278 hPtVsEta->Draw("colz");
01279 hPtVsEta->SetAxisRange(0, mPtMax, "Y");
01280
01282 mMainPad->cd(2);
01283 TH2* hPtVsY = (TH2D*) getHistogram(Form("hPtVsY_0_%d", mGeantId));
01284 if(!hPtVsY) return kFALSE ;
01285
01286 hPtVsY->SetTitle("");
01287 hPtVsY->Draw("colz");
01288 hPtVsY->SetAxisRange(0, mPtMax, "Y");
01289
01291 mMainPad->cd(3);
01292 TH2* hPtVsPhi = (TH2D*) getHistogram(Form("hPtVsPhi_0_%d", mGeantId));
01293 if(!hPtVsPhi) return kFALSE ;
01294
01295 hPtVsPhi->SetTitle("");
01296 hPtVsPhi->Draw("colz");
01297 hPtVsPhi->SetAxisRange(0, mPtMax, "Y");
01298
01299 mCanvas->cd();
01300 mCanvas->Update();
01301
01302 mPDF->NewPage() ;
01303
01304 delete header ;
01305
01306
01307 header = initCanvas("MC track QA (1D)", 2, 2);
01308
01310 TH1* hPt = (TH1D*) hPtVsPhi->ProjectionY("hPtMc");
01311 TH1* hEta = (TH1D*) hPtVsEta->ProjectionX("hEtaMc");
01312 TH1* hY = (TH1D*) hPtVsY->ProjectionX("hYMc");
01313 TH1* hPhi = (TH1D*) hPtVsPhi->ProjectionX("hPhiMc");
01314 hPt->SetTitleOffset(1.0, "X");
01315 hPt->SetXTitle("p_{T} (GeV/c)");
01316 hEta->SetXTitle("#eta");
01317 hY->SetXTitle("rapidity");
01318
01319 hPt ->SetTitle("");
01320 hEta->SetTitle("");
01321 hY ->SetTitle("");
01322 hPhi->SetTitle("");
01323
01324 hPt ->Sumw2() ;
01325 hEta->Sumw2() ;
01326 hY ->Sumw2() ;
01327 hPhi->Sumw2() ;
01328
01329 hPt->SetMinimum(0.0);
01330 hPt->SetMaximum(hPt->GetMaximum()*1.2);
01331 hPt->SetAxisRange(0, mPtMax, "X");
01332 hEta->SetMinimum(0.0);
01333 hEta->SetMaximum(hEta->GetMaximum()*1.2);
01334 hY->SetMinimum(0.0);
01335 hY->SetMaximum(hY->GetMaximum()*1.2);
01336 hPhi->SetMinimum(0.0);
01337 hPhi->SetMaximum(hPhi->GetMaximum()*1.2);
01338
01339
01340
01341
01342
01343 gStyle->SetOptFit(1);
01344 for(Int_t ipad=0; ipad<4; ipad++){
01345 mMainPad->cd(ipad+1);
01346 TH1* hdraw = 0 ;
01347 if( ipad == 0 ) hdraw = hPt ;
01348 if( ipad == 1 ) hdraw = hEta ;
01349 if( ipad == 2 ) hdraw = hY ;
01350 if( ipad == 3 ) hdraw = hPhi ;
01351
01352 hdraw->Draw();
01353
01354 if( ipad != 0 ){
01355
01356 Double_t min = -1.0 ;
01357 Double_t max = 1.0 ;
01358 if( ipad == 3 ){ min = -TMath::Pi() ; max = TMath::Pi() ; }
01359
01360 TF1* pol0Fit = new TF1(Form("pol0Fit_%d", ipad), "pol0", min, max);
01361 pol0Fit->SetLineColor(kRed);
01362 pol0Fit->SetLineWidth(1);
01363 pol0Fit->SetLineStyle(2);
01364 pol0Fit->SetParameter(0, hdraw->GetMean(2));
01365 hdraw->Fit(pol0Fit, "rq0");
01366 pol0Fit->Draw("same");
01367 }
01368 }
01369
01370 mCanvas->cd();
01371 mCanvas->Update();
01372
01373 mPDF->NewPage() ;
01374
01375 delete header ;
01376
01377 gStyle->SetOptFit(0);
01378
01379 return kTRUE ;
01380 }
01381
01382
01383 Bool_t StEmbeddingQADraw::drawTrack() const
01384 {
01386 LOG_INFO << endm ;
01387 LOG_INFO << "#----------------------------------------------------------------------------------------------------" << endm ;
01388 LOG_INFO << "QA for reconstructed tracks ..." << endm ;
01389
01391 if(!isOpen()) return kFALSE ;
01392
01393
01394
01396 const Bool_t isGeantIdOk = drawGeantId();
01397
01399 const Bool_t isPhiOk = drawPhi();
01400
01402 const Bool_t isRapidityOk = drawRapidity();
01403
01405 const Bool_t isMomentumOk = drawMomentum();
01406 const Bool_t isPtOk = drawPt();
01407
01409 const Bool_t isdEdxOk = drawdEdx();
01410
01412 const Bool_t isDcaOk = drawDca();
01413
01415 const Bool_t isNHitOk = drawNHit();
01416
01417 LOG_INFO << "#----------------------------------------------------------------------------------------------------" << endm ;
01418 LOG_INFO << endm ;
01419
01420 return isGeantIdOk && isPhiOk && isRapidityOk && isMomentumOk && isPtOk && isdEdxOk && isDcaOk && isNHitOk ;
01421 }
01422
01423
01424 Bool_t StEmbeddingQADraw::drawGeantId() const
01425 {
01427
01428 LOG_INFO << "QA for geant id ..." << endm;
01429
01431 if(!isOpen()) return kFALSE ;
01432
01436 TPaveText* header = initCanvas("Geant id", 1, 2);
01437
01439 TPad* topPad = (TPad*) mMainPad->cd(1);
01440 topPad->Divide(2, 1);
01441 topPad->cd(1);
01442 TH1* hGeantIdMc = (TH1D*) getHistogram("hGeantId_0");
01443 if(!hGeantIdMc) return kFALSE ;
01444
01446 const Int_t mcIdBin = hGeantIdMc->GetMaximumBin() ;
01447 const Double_t mcId = hGeantIdMc->GetBinCenter(mcIdBin);
01448 const Double_t mcIdMin = (mcId-10<0) ? 0 : mcId-10 ;
01449 const Double_t mcIdMax = mcId + 10 ;
01450
01451 hGeantIdMc->SetTitle(Form("MC Geant id for %s", getParticleName()));
01452 hGeantIdMc->SetMaximum(hGeantIdMc->GetMaximum()*1.2);
01453 hGeantIdMc->SetAxisRange(mcIdMin, mcIdMax, "X");
01454 hGeantIdMc->Draw();
01455
01457 TLine* hGeantIdMcExpected = new TLine(mGeantId+0.5, 0, mGeantId+0.5, hGeantIdMc->GetMaximum()) ;
01458 hGeantIdMcExpected->SetLineColor(kBlue);
01459 hGeantIdMcExpected->Draw();
01460
01462 topPad->cd(2);
01463 TH1* hGeantIdReco = (TH1D*) getHistogram(Form("hGeantId_%d", getCategoryId()));
01464 if(!hGeantIdReco) return kFALSE ;
01465
01466 hGeantIdReco->SetTitle("Reconstructed geant id");
01467 hGeantIdReco->SetMaximum(hGeantIdReco->GetMaximum()*1.2);
01468 hGeantIdReco->SetAxisRange(0, 50, "X");
01469 hGeantIdReco->Draw();
01470
01472 TLine** hGeantIdRecoExpected = new TLine*[mDaughterGeantId.size()];
01473
01474 for(UInt_t id=0; id<mDaughterGeantId.size(); id++){
01475 const Int_t geantid = mDaughterGeantId[id] ;
01476 hGeantIdRecoExpected[id] = new TLine(geantid+0.5, 0, geantid+0.5, hGeantIdReco->GetMaximum()) ;
01477 hGeantIdRecoExpected[id]->SetLineColor(kRed);
01478 hGeantIdRecoExpected[id]->SetLineStyle(id+1);
01479 hGeantIdRecoExpected[id]->Draw();
01480 }
01481
01482 mMainPad->cd(2);
01483 TLegend* leg = new TLegend(0.1, 0.2, 0.9, 0.8);
01484 leg->SetTextSize(0.05);
01485 leg->SetFillColor(10);
01486 leg->SetHeader("Particle informations");
01487
01488 leg->AddEntry(hGeantIdMcExpected, getMcParticleName(), "L");
01489
01490 for(UInt_t id=0; id<mDaughterGeantId.size(); id++){
01491 leg->AddEntry( hGeantIdRecoExpected[id], getEmbeddingParticleName(id), "L");
01492 }
01493 leg->Draw();
01494
01495 mCanvas->cd();
01496 mCanvas->Update();
01497
01498 mPDF->NewPage() ;
01499
01500 delete header ;
01501
01502 return kTRUE ;
01503 }
01504
01505
01506 Bool_t StEmbeddingQADraw::drawPhi() const
01507 {
01509
01511 if(!isOpen()) return kFALSE ;
01512
01514 Bool_t isPhiEmbeddingVsMcOk = drawProjection2D("phi", kTRUE);
01515
01517 Bool_t isPhiEmbeddingVsRealOk = drawProjection2D("phi", kFALSE);
01518
01519 return isPhiEmbeddingVsMcOk && isPhiEmbeddingVsRealOk ;
01520 }
01521
01522
01523 Bool_t StEmbeddingQADraw::drawRapidity() const
01524 {
01526
01528 if(!isOpen()) return kFALSE ;
01529
01531 const Bool_t isEtaOk = drawProjection2D("eta");
01532
01534 const Bool_t isYOk = drawProjection2D("y");
01535
01536 return isEtaOk && isYOk ;
01537 }
01538
01539
01540 Bool_t StEmbeddingQADraw::drawMomentum() const
01541 {
01543
01545 if(!isOpen()) return kFALSE ;
01546
01548 for(UInt_t id=0; id<mDaughterGeantId.size(); id++){
01549 TPaveText* header = 0 ;
01550
01551 TH2* hRecoPVsMcP = (TH2D*) getHistogram("hRecoPVsMcP", id, kTRUE);
01552 if(!hRecoPVsMcP){
01553
01554 header = initCanvas("Reconstructed momentum vs MC momentum");
01555 drawErrorMessages(getHistogramName("hRecoPVsMcP", id, kTRUE)) ;
01556
01557 mCanvas->cd();
01558 mCanvas->Update();
01559 mPDF->NewPage();
01560
01561 delete header ;
01562
01563 return kFALSE ;
01564 }
01565
01566 header = initCanvas("Reconstructed momentum vs MC momentum");
01567
01568 hRecoPVsMcP->SetAxisRange(0, mPtMax, "X");
01569 hRecoPVsMcP->SetAxisRange(0, mPtMax, "Y");
01570
01571 hRecoPVsMcP->SetTitle(getParticleName(mDaughterGeantId[id]));
01572
01573 hRecoPVsMcP->Draw("colz");
01574
01575 mCanvas->cd();
01576 mCanvas->Update();
01577
01578 mPDF->NewPage() ;
01579
01580 delete header ;
01581 }
01582
01583 return drawProjection2D("momentum");
01584 }
01585
01586
01587 Bool_t StEmbeddingQADraw::drawPt() const
01588 {
01590
01592 if(!isOpen()) return kFALSE ;
01593
01594
01595 for(UInt_t id=0; id<mDaughterGeantId.size(); id++){
01596 TPaveText* header = 0 ;
01597
01598 TH2* hdInvPtVsPt = (TH2D*) getHistogram("hdInvPtVsPt", id, kTRUE);
01599 if(!hdInvPtVsPt){
01600
01601 header = initCanvas(Form("1/p_{T} (Gl) - 1/p_{T} (MC) vs p_{T} (MC) (%s)", getParticleName(mDaughterGeantId[id])));
01602 drawErrorMessages(getHistogramName("hdInvPtVsPt", id, kTRUE)) ;
01603
01604 mCanvas->cd();
01605 mCanvas->Update();
01606 mPDF->NewPage();
01607
01608 delete header ;
01609
01610 return kFALSE ;
01611 }
01612
01613 gStyle->SetPadRightMargin(0.15);
01614
01615 header = initCanvas(Form("1/p_{T} (Gl) - 1/p_{T} (MC) vs p_{T} (MC) (%s)", getParticleName(mDaughterGeantId[id])), 1, 2);
01616
01617 mMainPad->cd(1);
01618 hdInvPtVsPt->SetAxisRange(0, mPtMax, "X");
01619 hdInvPtVsPt->SetTitle("");
01620 hdInvPtVsPt->Draw("colz");
01621
01622
01623 mMainPad->cd(2);
01624 TProfile* pdInvPtVsPt = (TProfile*) hdInvPtVsPt->ProfileX(Form("p%s", hdInvPtVsPt->GetName()));
01625 pdInvPtVsPt->SetAxisRange(0, mPtMax, "X");
01626 pdInvPtVsPt->SetTitle("");
01627 pdInvPtVsPt->SetYTitle(hdInvPtVsPt->GetYaxis()->GetTitle());
01628 pdInvPtVsPt->SetLineWidth(2);
01629 pdInvPtVsPt->SetMinimum(hdInvPtVsPt->GetYaxis()->GetXmin());
01630 pdInvPtVsPt->SetMaximum(hdInvPtVsPt->GetYaxis()->GetXmax());
01631 pdInvPtVsPt->Draw();
01632
01633 mCanvas->cd();
01634 mCanvas->Update();
01635
01636 mPDF->NewPage() ;
01637
01638 delete header ;
01639 }
01640
01641 gStyle->SetPadRightMargin(0.05);
01642
01643 return drawProjection2D("pt");
01644 }
01645
01646
01647 Bool_t StEmbeddingQADraw::drawProjection2D(const TString name, const Bool_t isMC) const
01648 {
01650
01653
01663
01664 TString nameLower(name);
01665 nameLower.ToLower();
01666
01667 LOG_INFO << "QA for " << name << " ..." << endm;
01668
01670 TString histoName("");
01671 Bool_t isProjectionX = kFALSE ;
01672 TString yTitle("");
01673 TString xTitle("");
01674
01675 if( nameLower.Contains("pt") ){
01676 xTitle = "#eta";
01677 yTitle = "p_{T}";
01678 histoName = "hPtVsEta";
01679 isProjectionX = kFALSE ;
01680 }
01681 else if( nameLower.Contains("momentum") ){
01682 xTitle = "#eta";
01683 yTitle = "p";
01684 histoName = "hMomVsEta";
01685 isProjectionX = kFALSE ;
01686 }
01687 else if( nameLower.Contains("eta") || name.Contains("pseudorapidity") ){
01688 xTitle = "p_{T}";
01689 yTitle = "#eta";
01690 histoName = "hPtVsEta";
01691 isProjectionX = kTRUE ;
01692 }
01693 else if( nameLower.Contains("y") || name.Contains("rapidity") ){
01694 xTitle = "p_{T}";
01695 yTitle = "y";
01696 histoName = "hPtVsY";
01697 isProjectionX = kTRUE ;
01698 }
01699 else if( nameLower.Contains("phi")){
01700 xTitle = "p_{T}";
01701 yTitle = "#phi";
01702 histoName = "hPtVsPhi";
01703 isProjectionX = kTRUE ;
01704 }
01705 else{
01706 Error("DrawProjection2D", "Unknown variable, %s", name.Data());
01707 LOG_INFO << " Current implemented variables are" << endm;
01708 LOG_INFO << "-----------------------------------------------------------------" << endm;
01709 LOG_INFO << " Input variable" << endm;
01710 LOG_INFO << "-----------------------------------------------------------------" << endm;
01711 LOG_INFO << " pt p_{T}" << endm;
01712 LOG_INFO << " momentum p" << endm;
01713 LOG_INFO << " eta or pseudorapidity eta" << endm;
01714 LOG_INFO << " y or rapidity y" << endm;
01715 LOG_INFO << " phi phi" << endm;
01716 LOG_INFO << "-----------------------------------------------------------------" << endm;
01717 LOG_INFO << endm;
01718 LOG_INFO << "NOTE : Input is case insensitive" << endm;
01719 LOG_INFO << endm;
01720
01721 return kFALSE ;
01722 }
01723
01724 gStyle->SetPadRightMargin(0.05);
01725
01726
01727 TH2* h2DMc = (TH2D*) getHistogram(Form("%s_0_%d", histoName.Data(), mGeantId));
01728
01729 for(UInt_t id=0; id<mDaughterGeantId.size(); id++){
01730 const TString headerTitle(Form("Projection of %s for each %s bin", yTitle.Data(), xTitle.Data()));
01731 TPaveText* header = 0 ;
01732
01733 TH2* h2DEmbed = (TH2D*) getHistogram(histoName, id, kTRUE);
01734 if(!h2DEmbed){
01735
01736 header = initCanvas(headerTitle);
01737 drawErrorMessages(getHistogramName(histoName, id, kTRUE)) ;
01738
01739 mCanvas->cd();
01740 mCanvas->Update();
01741 mPDF->NewPage();
01742
01743 delete header ;
01744
01745 return kFALSE ;
01746 }
01747
01748 TH2* h2DReal = (TH2D*) getHistogram(histoName, id, kFALSE);
01749
01751 header = initCanvas(headerTitle);
01752
01753 Int_t npad = 0 ;
01754 Int_t npadMax = 0 ;
01755 if( isProjectionX ){
01756 mMainPad->Divide(2, 3);
01757 npad = 5 ;
01758 npadMax = 6 ;
01759 }
01760 else{
01761 mMainPad->Divide(2, 4);
01762 npad = 6 ;
01763 npadMax = 8 ;
01764 }
01765
01769 const Int_t nbins = (isProjectionX) ? 10 : 6 ;
01770 const Double_t binStep = 0.5 ;
01771 const Double_t binMin = (isProjectionX) ? 0.0 : -1.5 ;
01772
01773 Int_t ipad = 1 ;
01774 for(Int_t ibin=0; ibin<nbins; ibin++){
01775 if( ipad%(npad+1) == 0 ){
01776 mCanvas->cd();
01777 mCanvas->Update();
01778 mPDF->NewPage();
01779
01780 if(header) delete header ;
01781 header = initCanvas(headerTitle);
01782
01783 if( isProjectionX ) mMainPad->Divide(2, 3);
01784 else mMainPad->Divide(2, 4);
01785 ipad = 1 ;
01786 }
01787
01788 const Double_t xMin = (isProjectionX && ibin==0) ? 0.2 : binMin + binStep * ibin ;
01789 const Double_t xMax = (isProjectionX && ibin==0) ? 0.5 : binMin + binStep * (ibin+1.0) ;
01790 const Int_t xMinBin = (isProjectionX) ? h2DEmbed->GetYaxis()->FindBin(xMin) : h2DEmbed->GetXaxis()->FindBin(xMin) ;
01791 const Int_t xMaxBin = (isProjectionX) ? h2DEmbed->GetYaxis()->FindBin(xMax) - 1 : h2DEmbed->GetXaxis()->FindBin(xMax) - 1;
01792
01793 TH1* hEmbed = 0;
01794 TH1* hReal = 0;
01795 TH1* hMc = 0;
01796 if( isProjectionX ){
01797 if( h2DEmbed ) hEmbed = (TH1D*) h2DEmbed->ProjectionX(Form("h%sEmbed_%d_%d", name.Data(), id, ibin), xMinBin, xMaxBin);
01798 if( h2DReal ) hReal = (TH1D*) h2DReal->ProjectionX(Form("h%sReal_%d_%d", name.Data(), id, ibin), xMinBin, xMaxBin);
01799 if( h2DMc ) hMc = (TH1D*) h2DMc->ProjectionX(Form("h%sMc_%d_%d", name.Data(), id, ibin), xMinBin, xMaxBin);
01800 }
01801 else{
01802 if( h2DEmbed ) hEmbed = (TH1D*) h2DEmbed->ProjectionY(Form("h%sEmbed_%d_%d", name.Data(), id, ibin), xMinBin, xMaxBin);
01803 if( h2DReal ) hReal = (TH1D*) h2DReal->ProjectionY(Form("h%sReal_%d_%d", name.Data(), id, ibin), xMinBin, xMaxBin);
01804 if( h2DMc ) hMc = (TH1D*) h2DMc->ProjectionY(Form("h%sMc_%d_%d", name.Data(), id, ibin), xMinBin, xMaxBin);
01805 }
01806
01807 hEmbed->SetLineColor(kRed);
01808 if(hEmbed->GetSumw2N()==0) hEmbed->Sumw2();
01809 hEmbed->Scale( getNormalization(*hEmbed) );
01810
01811 if( hReal ){
01812 hReal->SetLineColor(kBlue);
01813 if(hReal->GetSumw2N()==0) hReal->Sumw2();
01814 hReal->Scale( getNormalization(*hReal) );
01815 }
01816
01817 if( hMc ){
01818 hMc->SetLineColor(kBlack);
01819 if(hMc->GetSumw2N()==0) hMc->Sumw2();
01820 hMc->Scale( getNormalization(*hMc) );
01821 }
01822
01823 hEmbed->SetMinimum(0.0);
01824
01825 if( isProjectionX ){
01826 hEmbed->SetTitle(Form("%1.1f < p_{T} < %1.1f (GeV/c)", xMin, xMax));
01827 }
01828 else{
01829 hEmbed->SetTitle(Form("%1.1f < #eta < %1.1f", xMin, xMax));
01830 hEmbed->SetAxisRange(0, mPtMax, "X");
01831 }
01832 hEmbed->SetYTitle(Form("(1/N_{trk})dN/d%s", yTitle.Data()));
01833
01834 mMainPad->cd(ipad) ;
01835 hEmbed->Draw();
01836
01838 if(isMC && !isDecay()){
01840 hEmbed->SetMaximum( TMath::Max(hMc->GetMaximum(), hEmbed->GetMaximum()) * 1.2 );
01841
01842 hMc->Draw("same");
01843 }
01844 else{
01848 const Double_t yMax = (hReal) ? TMath::Max(hReal->GetMaximum(), hEmbed->GetMaximum()) : hEmbed->GetMaximum() ;
01849 hEmbed->SetMaximum( yMax * 1.2 );
01850
01851 if(hReal) hReal->Draw("hsame");
01852 }
01853 hEmbed->Draw("same");
01854
01855
01856 if ( ipad == 1 ){
01857 mMainPad->cd(npadMax);
01858 TLegend* leg = new TLegend(0.0, 0.2, 1.0, 0.8);
01859 leg->SetFillColor(10);
01860 leg->SetTextFont(43);
01861 leg->SetTextSize(15);
01862
01863 leg->AddEntry(hEmbed, getEmbeddingParticleName(id, kTRUE), "L");
01864 if(isMC && !isDecay()){
01865 leg->AddEntry(hMc, getMcParticleName(), "L");
01866 }
01867 else{
01868 if(hReal) leg->AddEntry(hReal, getRealParticleName(id, kTRUE), "L");
01869 }
01870 leg->Draw();
01871 }
01872
01873 ipad++;
01874 }
01875
01876 mCanvas->cd();
01877 mCanvas->Update();
01878
01879 mPDF->NewPage();
01880
01881 if(header) delete header ;
01882 }
01883
01884 gStyle->SetPadRightMargin(0.15);
01885
01886 return kTRUE ;
01887 }
01888
01889
01890 Bool_t StEmbeddingQADraw::drawdEdx() const
01891 {
01893
01895 if(!isOpen()) return kFALSE ;
01896
01898 const Bool_t isMcOk = drawdEdxVsMomentum(kTRUE) ;
01899
01901 const Bool_t isRecoOk = drawdEdxVsMomentum(kFALSE) ;
01902
01903 return isMcOk && isRecoOk ;
01904 }
01905
01906
01907 Bool_t StEmbeddingQADraw::drawdEdxVsMomentum(const Bool_t isMcMomentum) const
01908 {
01912
01914
01917 const TString momName = (isMcMomentum) ? "Mc" : "Reco" ;
01918
01919 LOG_INFO << "QA for dE/dx (" << momName << " momentum ...)" << endm;
01920
01921 gStyle->SetOptFit(0);
01922 gStyle->SetPadRightMargin(0.05);
01923
01924 for(UInt_t id=0; id<mDaughterGeantId.size(); id++){
01925
01926
01927
01928 TString headerTitle("");
01929 if( mIsEmbeddingOnly ){
01930 headerTitle = Form("dE/dx vs momentum (Embedding:%s)", getParticleName(mDaughterGeantId[id]) );
01931 }
01932 else{
01933 headerTitle = Form("dE/dx vs momentum (Embedding:%s, Real:%s)",
01934 getParticleName(mDaughterGeantId[id]), getParticleName(getGeantIdReal(id)) );
01935 }
01936 TPaveText* header = 0 ;
01937
01938
01939 TH2* hdEdxVsMomEmbed = (TH2D*) getHistogram(Form("hdEdxVsMom%s", momName.Data()), id, kTRUE);
01940 if(!hdEdxVsMomEmbed){
01941
01942 header = initCanvas(headerTitle);
01943 drawErrorMessages(getHistogramName(Form("hdEdxVsMom%s", momName.Data()), id, kTRUE)) ;
01944
01945 mCanvas->cd();
01946 mCanvas->Update();
01947 mPDF->NewPage();
01948
01949 delete header ;
01950
01951 return kFALSE ;
01952 }
01953
01954 header = initCanvas(headerTitle, 1, 2);
01955
01956 hdEdxVsMomEmbed->SetLineColor(kRed);
01957 hdEdxVsMomEmbed->SetMarkerColor(kRed);
01958
01959
01960 TH2* hdEdxVsMomReal = (TH2D*) getHistogram("hdEdxVsMomReco", id, kFALSE);
01961
01962
01963 if( hdEdxVsMomReal ){
01964 hdEdxVsMomReal->SetLineColor(kBlack);
01965 hdEdxVsMomReal->SetMarkerColor(kBlack);
01966 }
01967
01968
01969 TH2* hdEdxVsMomPidReal = (TH2D*) getHistogram("hdEdxVsMomRecoPidCut", id, kFALSE);
01970
01971
01972 if( hdEdxVsMomPidReal ){
01973 hdEdxVsMomPidReal->SetLineColor(kBlue);
01974 hdEdxVsMomPidReal->SetMarkerColor(kBlue);
01975 }
01976
01977
01978 mMainPad->cd(1);
01979 hdEdxVsMomEmbed->SetAxisRange(0, mPtMax, "X");
01980 hdEdxVsMomEmbed->SetXTitle(Form("%s momentum (GeV/c)", momName.Data()));
01981 hdEdxVsMomEmbed->SetTitle("");
01982 hdEdxVsMomEmbed->SetMinimum(0);
01983
01984 hdEdxVsMomEmbed->Draw("box");
01985 if( hdEdxVsMomReal ) hdEdxVsMomReal->Draw("box same");
01986 if( hdEdxVsMomPidReal ) hdEdxVsMomPidReal->Draw("box same");
01987 hdEdxVsMomEmbed->Draw("box same");
01988
01989 mMainPad->cd(2);
01990
01991 TLegend* leg = new TLegend(0.1, 0.25, 0.9, 0.9);
01992 leg->SetFillColor(10);
01993 leg->SetTextSize(0.05);
01994 leg->AddEntry( hdEdxVsMomEmbed, getEmbeddingParticleName(id), "L");
01995 if( hdEdxVsMomReal ) leg->AddEntry( hdEdxVsMomReal, "Real data", "L");
01996 if( hdEdxVsMomPidReal ) leg->AddEntry( hdEdxVsMomPidReal, "Real data with PID cut (#sigma<2)", "L");
01997 leg->Draw();
01998
01999 mCanvas->cd();
02000 mCanvas->Update();
02001
02002 mPDF->NewPage();
02003
02004 delete header ;
02005
02006
02009
02010 headerTitle = "Projection of dE/dx for each p bin";
02011
02012 header = initCanvas(headerTitle, 2, 4);
02013 Int_t ipad = 1 ;
02014
02015
02016 const Int_t npad = 6 ;
02017 const Int_t npadMax = 8;
02018
02019 TGraphErrors* gMeanVsMom[2];
02020 TGraphErrors* gSigmaVsMom[2];
02021
02022 for(Int_t i=0; i<2; i++){
02023 gMeanVsMom[i] = new TGraphErrors();
02024 gSigmaVsMom[i] = new TGraphErrors();
02025 gMeanVsMom[i] ->SetMarkerStyle(24 - i*4);
02026 gSigmaVsMom[i]->SetMarkerStyle(24 - i*4);
02027 }
02028 gMeanVsMom[0] ->SetMarkerColor(kRed); gMeanVsMom[0] ->SetLineColor(kRed);
02029 gSigmaVsMom[0]->SetMarkerColor(kRed); gSigmaVsMom[0]->SetLineColor(kRed);
02030 gMeanVsMom[1] ->SetMarkerColor(kBlue); gMeanVsMom[1] ->SetLineColor(kBlue);
02031 gSigmaVsMom[1]->SetMarkerColor(kBlue); gSigmaVsMom[1]->SetLineColor(kBlue);
02032
02033 const Int_t npt = 24 ;
02034 const Double_t ptBin = 0.2 ;
02035 for(Int_t ipt=0; ipt<npt; ipt++){
02036 if( ipad % (npad+1) == 0 ){
02037 mCanvas->cd();
02038 mCanvas->Update();
02039
02040 mPDF->NewPage();
02041
02042 if(header) delete header ;
02043 header = initCanvas(headerTitle, 2, 4);
02044 ipad = 1 ;
02045 }
02046
02047 const Double_t ptMin = 0.2 + ipt*ptBin ;
02048 const Double_t ptMax = ptMin + ptBin ;
02049 const Int_t ptMinBin = hdEdxVsMomEmbed->GetXaxis()->FindBin(ptMin);
02050 const Int_t ptMaxBin = hdEdxVsMomEmbed->GetXaxis()->FindBin(ptMax-0.001);
02051 if( ptMinBin == ptMaxBin ){
02052 LOG_INFO << Form("%1.1f - %1.1f GeV/c : bin = (%4d, %4d)", ptMin, ptMax, ptMinBin, ptMaxBin) << endm;
02053 }
02054
02055
02056 TH1* hdEdxEmbed = (TH1D*) hdEdxVsMomEmbed->ProjectionY(Form("hdEdxEmbed%s_%d_%d", momName.Data(), id, ipt), ptMinBin, ptMaxBin);
02057 TH1* hdEdxReal = 0;
02058 if( hdEdxVsMomPidReal )
02059 hdEdxReal = (TH1D*) hdEdxVsMomPidReal->ProjectionY(Form("hdEdxReal%s_%d_%d", momName.Data(), id, ipt), ptMinBin, ptMaxBin);
02060
02061 hdEdxEmbed->Sumw2() ;
02062 hdEdxEmbed->Scale( getNormalization(*hdEdxEmbed) ) ;
02063
02064 if( hdEdxReal ){
02065 hdEdxReal->Sumw2() ;
02066 hdEdxReal->Scale( getNormalization(*hdEdxReal) ) ;
02067 }
02068
02069 hdEdxEmbed->SetMinimum(0.0);
02070
02071
02072
02073 if( hdEdxReal ){
02074 hdEdxEmbed->SetMaximum( TMath::Max(hdEdxReal->GetMaximum(), hdEdxEmbed->GetMaximum())*1.2 );
02075 }
02076 else{
02077 hdEdxEmbed->SetMaximum( hdEdxEmbed->GetMaximum()*1.2 );
02078 }
02079
02080 hdEdxEmbed->SetTitle(Form("%1.1f < %s p < %1.1f GeV/c", ptMin, momName.Data(), ptMax));
02081 hdEdxEmbed->SetYTitle("(1/N_{trk})dN/d(dE/dx)");
02082
02083 mMainPad->cd(ipad);
02084 hdEdxEmbed->Draw("h");
02085 if( hdEdxReal ) hdEdxReal->Draw("hsame");
02086 hdEdxEmbed->Draw("hsame");
02087
02088
02089
02090
02091 TF1 fEmbed(Form("fEmbed%s_%d_%d", momName.Data(), id, ipt), "gaus", 0, 10);
02092 TF1 fReal(Form("fReal%s_%d_%d", momName.Data(), id, ipt), "gaus", 0, 10);
02093 fEmbed.SetLineColor(kRed);
02094 fReal.SetLineColor(kBlue);
02095 fEmbed.SetLineWidth(1);
02096 fReal.SetLineWidth(1);
02097
02098 if( hdEdxReal ){
02099 hdEdxReal->Fit(fReal.GetName(), "rq0");
02100 fReal.Draw("same");
02101 }
02102 hdEdxEmbed->Fit(fEmbed.GetName(), "rq0");
02103 fEmbed.Draw("same");
02104
02105 const Double_t pt = (ptMin+ptMax)/2.0 ;
02106 gMeanVsMom[0]->SetPoint(ipt, pt, fEmbed.GetParameter(1));
02107 gMeanVsMom[0]->SetPointError(ipt, 0.0, fEmbed.GetParError(1));
02108 gSigmaVsMom[0]->SetPoint(ipt, pt, fEmbed.GetParameter(2));
02109 gSigmaVsMom[0]->SetPointError(ipt, 0.0, fEmbed.GetParError(2));
02110
02111 if( hdEdxReal ){
02112 gMeanVsMom[1]->SetPoint(ipt, pt, fReal.GetParameter(1));
02113 gMeanVsMom[1]->SetPointError(ipt, 0.0, fReal.GetParError(1));
02114 gSigmaVsMom[1]->SetPoint(ipt, pt, fReal.GetParameter(2));
02115 gSigmaVsMom[1]->SetPointError(ipt, 0.0, fReal.GetParError(2));
02116 }
02117 else{
02118 gMeanVsMom[1]->SetPoint(ipt, pt, -9999.);
02119 gMeanVsMom[1]->SetPointError(ipt, 0.0, 0.0);
02120 gSigmaVsMom[1]->SetPoint(ipt, pt, -9999.);
02121 gSigmaVsMom[1]->SetPointError(ipt, 0.0, 0.0);
02122 }
02123
02124
02125 if( ipad == 1 ){
02126 mMainPad->cd(npadMax);
02127 drawLegend(id, hdEdxEmbed, hdEdxReal, "L", kTRUE) ;
02128
02129
02130
02131
02132
02133
02134
02135 }
02136
02137 ipad++;
02138 }
02139
02140 mCanvas->cd();
02141 mCanvas->Update();
02142
02143 mPDF->NewPage();
02144
02145 delete header ;
02146
02147
02149
02150 headerTitle = "Mean/#sigma of dE/dx vs momentum";
02151 header = initCanvas(headerTitle, 1, 2);
02152
02153 for(Int_t i=0; i<2; i++){
02154 mMainPad->cd(i+1);
02155
02156 const Double_t ymax = (i==0) ? 7.2 : 1.4 ;
02157 TH1* frame = mMainPad->GetPad(i+1)->DrawFrame(0, 0.0, 5.0, ymax);
02158 frame->SetXTitle(Form("%s momentum (GeV/c)", momName.Data()));
02159 if( i == 0 ) frame->SetYTitle("Mean (keV/cm)");
02160 if( i == 1 ) frame->SetYTitle("#sigma (keV/cm)");
02161
02162 TLegend* leg = new TLegend(0.23, 0.86, 0.92, 0.99);
02163 leg->SetBorderSize(1);
02164 leg->SetTextFont(43);
02165 leg->SetTextSize(15);
02166 leg->SetFillColor(10);
02167
02168 if( i == 0 ){
02169 gMeanVsMom[0]->Draw("P");
02170 gMeanVsMom[1]->Draw("P");
02171
02172 leg->AddEntry( gMeanVsMom[0], getEmbeddingParticleName(id), "P");
02173 leg->AddEntry( gMeanVsMom[1], getRealParticleName(id), "P");
02174 }
02175 if( i == 1 ){
02176 gSigmaVsMom[0]->Draw("P");
02177 gSigmaVsMom[1]->Draw("P");
02178
02179 leg->AddEntry( gSigmaVsMom[0], getEmbeddingParticleName(id), "P");
02180 leg->AddEntry( gSigmaVsMom[1], getRealParticleName(id), "P");
02181 }
02182 leg->Draw();
02183 }
02184
02185 mCanvas->cd();
02186 mCanvas->Update();
02187
02188 mPDF->NewPage() ;
02189
02190 delete header ;
02191 }
02192
02193 return kTRUE ;
02194 }
02195
02196
02197 Bool_t StEmbeddingQADraw::drawDca() const
02198 {
02200
02201 LOG_INFO << "QA for dca distributions ..." << endm;
02202
02204 if(!isOpen()) return kFALSE ;
02205
02206 return drawProjection3D("Dca");
02207 }
02208
02209
02210 Bool_t StEmbeddingQADraw::drawNHit() const
02211 {
02213
02214 LOG_INFO << "QA for NHit distributions ..." << endm;
02215
02217 if(!isOpen()) return kFALSE ;
02218
02219 gStyle->SetPadRightMargin(0.17);
02220
02222 for(UInt_t id=0; id<mDaughterGeantId.size(); id++){
02223 TH3* hNCommonHitVsNHit = (TH3D*) getHistogram("hNCommonHitVsNHit", id, kTRUE);
02224 if ( !hNCommonHitVsNHit ) continue ;
02225
02226 TString headerTitle("");
02227 if( mIsEmbeddingOnly ){
02228 headerTitle = Form("N_{common} vs N_{hit} (Embedding:%s)", getParticleName(mDaughterGeantId[id]) );
02229 }
02230 else{
02231 headerTitle = Form("N_{common} vs N_{hit} (Embedding:%s, Real:%s)",
02232 getParticleName(mDaughterGeantId[id]), getParticleName(getGeantIdReal(id)) );
02233 }
02234
02235
02236 TPaveText* header = initCanvas(headerTitle);
02237
02238
02239 hNCommonHitVsNHit->SetAxisRange(0.2, 5.0);
02240 TH2* hNCommonHitVsNHit2D = (TH2D*) hNCommonHitVsNHit->Project3D("zy");
02241 hNCommonHitVsNHit2D->SetName(Form("hNCommonHitVsNHit2D_%d", id));
02242 hNCommonHitVsNHit2D->SetTitle(Form("%1.1f < p_{T} < %1.1f GeV/c", 0.2, 5.0));
02243 hNCommonHitVsNHit2D->Draw("colz");
02244
02245 mCanvas->cd();
02246 mCanvas->Update();
02247 mPDF->NewPage();
02248 delete header ;
02249
02250
02251 if( mIsEmbeddingOnly ){
02252 headerTitle = Form("N_{common} vs N_{hit}, p_{T} dependence (Embedding:%s)", getParticleName(mDaughterGeantId[id]) );
02253 }
02254 else{
02255 headerTitle = Form("N_{common} vs N_{hit}, p_{T} dependence (Embedding:%s, Real:%s)",
02256 getParticleName(mDaughterGeantId[id]), getParticleName(getGeantIdReal(id)) );
02257 }
02258 for(Int_t jpt=0; jpt<2; jpt++) {
02259 TPaveText* header = initCanvas(headerTitle, 2, 3);
02260
02261 const Int_t npt = hNCommonHitVsNHit->GetNbinsX()/2 ;
02262 for(Int_t ipt=0; ipt<npt; ipt++) {
02263 mMainPad->cd(ipt+1);
02264 const Int_t ptId = jpt*npt + ipt + 1;
02265 Double_t ptmin = hNCommonHitVsNHit->GetXaxis()->GetBinLowEdge(ptId) ;
02266 Double_t ptmax = hNCommonHitVsNHit->GetXaxis()->GetBinLowEdge(ptId+1) ;
02267 if( ptmin == 0.0 ) ptmin = 0.2 ;
02268
02269 hNCommonHitVsNHit->SetAxisRange(ptmin, ptmax);
02270 hNCommonHitVsNHit2D = (TH2D*) hNCommonHitVsNHit->Project3D("zy");
02271 hNCommonHitVsNHit2D->SetName(Form("hNCommonHitVsNHit2D_%d_%d_%d", id, ipt, jpt));
02272 hNCommonHitVsNHit2D->SetTitle(Form("%1.1f < p_{T} < %1.1f GeV/c", ptmin, ptmax));
02273 hNCommonHitVsNHit2D->Draw("colz");
02274 }
02275 mCanvas->cd();
02276 mCanvas->Update();
02277 mPDF->NewPage();
02278 delete header ;
02279 }
02280 }
02281
02282 gStyle->SetPadRightMargin(0.05);
02283
02284 return drawProjection3D("NHit") ;
02285 }
02286
02287
02288 Bool_t StEmbeddingQADraw::drawProjection3D(const TString name) const
02289 {
02291
02294
02296 const Bool_t isBatch = gROOT->IsBatch() ;
02297 if( !isBatch ){
02298 LOG_INFO << "Enter batch mode ..." << endm;
02299 gROOT->SetBatch(kTRUE);
02300 }
02301
02302 TString nameLower(name);
02303 nameLower.ToLower();
02304
02305 TString headerTitle(Form("%s distribution for (p_{T}, #eta) slices", name.Data()));
02306 for(UInt_t id=0; id<mDaughterGeantId.size(); id++){
02307 TPaveText* header = 0 ;
02308
02310 TH3* h3DEmbed = (TH3D*) getHistogram(Form("h%s", name.Data()), id, kTRUE);
02311 if(!h3DEmbed){
02312
02313 header = initCanvas(headerTitle);
02314 drawErrorMessages(getHistogramName(Form("h%s", name.Data()), id, kTRUE));
02315
02316 mCanvas->cd();
02317 mCanvas->Update();
02318 mPDF->NewPage();
02319
02320 delete header ;
02321
02322 return kFALSE ;
02323 }
02324
02325 TH3* h3DReal = (TH3D*) getHistogram(Form("h%s", name.Data()), id, kFALSE);
02326
02327 const Int_t nPt = h3DEmbed->GetNbinsX() ;
02328 const Int_t nEta = h3DEmbed->GetNbinsY() ;
02329 const Double_t ptMin = h3DEmbed->GetXaxis()->GetXmin() ;
02330 const Double_t ptMax = h3DEmbed->GetXaxis()->GetXmax() ;
02331 const Double_t etaMin = h3DEmbed->GetYaxis()->GetXmin() ;
02332 const Double_t etaMax = h3DEmbed->GetYaxis()->GetXmax() ;
02333 const Double_t etaBin = (etaMax-etaMin)/(Double_t)nEta ;
02334 const Double_t ptBin = (ptMax-ptMin)/(Double_t)nPt ;
02335
02336 for(Int_t ipt=0; ipt<nPt; ipt++){
02337 TString pt(Form("%1.1f < p_{T} < %1.1f (GeV/c)", ptMin+ipt*ptBin, ptMin+(ipt+1)*ptBin));
02338 if( ipt == 0 ) pt = Form("%1.1f < p_{T} < %1.1f (GeV/c)", 0.1, ptMin+(ipt+1)*ptBin);
02339
02340 TPaveText* header = initCanvas(headerTitle, 2, 3);
02341
02342 const Int_t npad = 5 ;
02343 const Int_t npadMax = 6 ;
02344 Int_t ipad = 1 ;
02345 for(Int_t ieta=0; ieta<nEta; ieta++){
02346 if( ipad % (npad+1) == 0 ){
02347 mCanvas->cd();
02348 mCanvas->Update();
02349 mPDF->NewPage();
02350
02351 if(header) delete header ;
02352 header = initCanvas(headerTitle, 2, 3);
02353
02354 ipad = 1 ;
02355 }
02356
02357 TString eta(Form("%1.1f < #eta < %1.1f", etaMin+ieta*etaBin, etaMin+(ieta+1)*etaBin));
02358
02359 TH1* hEmbed = (TH1D*) h3DEmbed->ProjectionZ(Form("h%sEmbed_%d_%d_%d", name.Data(), id, ipt, ieta), ipt+1, ipt+1, ieta+1, ieta+1);
02360 TH1* hReal = 0 ;
02361 if( h3DReal ){
02362 hReal = (TH1D*) h3DReal->ProjectionZ(Form("h%sReal_%d_%d_%d", name.Data(), id, ipt, ieta), ipt+1, ipt+1, ieta+1, ieta+1);
02363 hReal ->Sumw2();
02364 hReal ->Scale( getNormalization(*hReal) );
02365 hReal ->SetLineColor(kBlue);
02366 }
02367 hEmbed->Sumw2();
02368 hEmbed->Scale( getNormalization(*hEmbed) );
02369 hEmbed->SetLineColor(kRed);
02370 hEmbed->SetMinimum(0.0);
02371
02372
02373
02374
02375 if( hReal ){
02376 hEmbed->SetMaximum( TMath::Max(hReal->GetMaximum(), hEmbed->GetMaximum()) * 1.2 );
02377 }
02378 else{
02379 hEmbed->SetMaximum( hEmbed->GetMaximum() * 1.2 );
02380 }
02381
02382 hEmbed->SetTitle(eta + ", " + pt);
02383 hEmbed->SetYTitle(Form("(1/N_{trk})dN/d%s", name.Data())) ;
02384
02385 mMainPad->cd(ipad);
02386 hEmbed->Draw();
02387 if(hReal) hReal->Draw("hsame");
02388 hEmbed->Draw("same");
02389
02390
02391 if( ipad == 1 ){
02392 mMainPad->cd(npadMax);
02393 drawLegend(id, hEmbed, hReal, "L", kTRUE);
02394 }
02395
02396 ipad++;
02397 }
02398 mCanvas->cd();
02399 mCanvas->Update();
02400 mPDF->NewPage();
02401
02402 delete header ;
02403 }
02404 }
02405
02406
02407
02408
02409
02410
02411
02412 if( isBatch ) gROOT->SetBatch(kTRUE);
02413 else gROOT->SetBatch(kFALSE);
02414
02415 return kTRUE ;
02416 }
02417
02418
02419 TPaveText* StEmbeddingQADraw::initCanvas(const TString headerTitle, const Int_t nx, const Int_t ny) const
02420 {
02421 mMainPad->Clear();
02422 mCanvas->cd();
02423
02425 TPaveText* header = drawHeader(headerTitle);
02426
02427 mMainPad->cd();
02428 if( nx != 0 && ny != 0 ) mMainPad->Divide(nx, ny);
02429
02430 return header ;
02431 }
02432
02433
02434 Double_t StEmbeddingQADraw::getVzAcceptedMinimum() const
02435 {
02437
02438
02439
02440
02441
02442
02443
02444
02445
02446
02447
02448
02449
02450
02451
02452 return -StEmbeddingQAUtilities::instance()->getZVertexCut() ;
02453 }
02454
02455
02456 Double_t StEmbeddingQADraw::getVzAcceptedMaximum() const
02457 {
02459
02460
02461
02462
02463
02464
02465
02466
02467
02468
02469
02470
02471
02472
02473
02474 return StEmbeddingQAUtilities::instance()->getZVertexCut() ;
02475 }
02476
02477
02478 Bool_t StEmbeddingQADraw::draw() const
02479 {
02481
02482
02483
02484
02485
02486
02487
02488
02490 const Bool_t isEventOk = drawEvent();
02491
02492
02493
02495 const Bool_t isMcTrackOk = drawMcTrack();
02496
02497
02498
02500 const Bool_t isTrackOk = drawTrack();
02501
02502
02503 return isEventOk && isMcTrackOk && isTrackOk ;
02504 }
02505
02506
02507 Bool_t StEmbeddingQADraw::finish()
02508 {
02509
02510 LOG_INFO << "StEmbeddingQADraw::finish() Close PDF file : " << mPDF->GetName() << endm ;
02511
02512 mMainPad->cd();
02513 mMainPad->Clear();
02514
02515 TLatex* title = new TLatex(0.1, 0.6, "End of QA");
02516 title->SetTextSize(0.1);
02517 title->Draw();
02518
02519 mCanvas->cd();
02520 mCanvas->Update();
02521
02522 mPDF->Close() ;
02523
02524 return kTRUE ;
02525 }
02526
02527 void StEmbeddingQADraw::setPNGOn() { LOG_INFO << endm << "Print figures to PNG file" << endm << endm; mIsPNGOn = kTRUE ; }
02528 void StEmbeddingQADraw::setGIFOn() { LOG_INFO << endm << "Print figures to GIF file" << endm << endm; mIsGIFOn = kTRUE ; }
02529 void StEmbeddingQADraw::setJPGOn() { LOG_INFO << endm << "Print figures to JPG file" << endm << endm; mIsJPGOn = kTRUE ; }
02530 void StEmbeddingQADraw::setEPSOn() { LOG_INFO << endm << "Print figures to EPS file" << endm << endm; mIsEPSOn = kTRUE ; }
02531 void StEmbeddingQADraw::setPSOn() { LOG_INFO << endm << "Print figures to PS file" << endm << endm; mIsPSOn = kTRUE ; }
02532
02533
02534
