StRoot: StEmbeddingUtilities/StEmbeddingQA.cxx Source File

00001 /****************************************************************************************************
00002  * $Id: StEmbeddingQA.cxx,v 1.22 2012/03/05 10:32:19 cpowell Exp $
00003  * $Log: StEmbeddingQA.cxx,v $
00004  * Revision 1.22  2012/03/05 10:32:19  cpowell
00005  * Functions added to cut on refMult
00006  *
00007  * Revision 1.21  2011/08/31 18:04:48  cpowell
00008  * Extended the range of hGeantId
00009  *
00010  * Revision 1.20  2011/04/01 05:05:49  hmasui
00011  * Track selections by StEmbeddingQAUtilities. Added 1/pt(RC)-1/pt(MC) vs pt, and pt dependent Ncommon vs NhitFit histograms
00012  *
00013  * Revision 1.19  2011/02/11 03:55:46  hmasui
00014  * Change geantid type to integer
00015  *
00016  * Revision 1.18  2011/01/31 21:32:12  hmasui
00017  * Modify histogram keys to TString to take into account parent geantid
00018  *
00019  * Revision 1.17  2011/01/14 23:46:16  hmasui
00020  * Add Ncommon hit cut for NHitFit histograms
00021  *
00022  * Revision 1.16  2011/01/12 21:36:29  hmasui
00023  * Add nHitsFit/nHitsPoss cut
00024  *
00025  * Revision 1.15  2010/11/01 03:10:23  hmasui
00026  * Modify geantid check for MC tracks in order to avoid non primary tracks
00027  *
00028  * Revision 1.14  2010/07/12 21:30:23  hmasui
00029  * Use StEmbeddingQAUtilities::getParticleDefinition(). Increase bin size, maximum for geantid histogram in order to cover geantid > 10k
00030  *
00031  * Revision 1.13  2010/06/28 17:33:47  hmasui
00032  * Added geant process = 6 (photon pair production) in contaminated pairs
00033  *
00034  * Revision 1.12  2010/05/14 19:50:12  hmasui
00035  * Add rapidity and trigger cuts.
00036  *
00037  * Revision 1.11  2010/04/24 20:21:21  hmasui
00038  * Add geant process check for contaminated pairs
00039  *
00040  * Revision 1.10  2010/03/15 21:03:56  hmasui
00041  * Modify binning for histograms of vertices
00042  *
00043  * Revision 1.9  2010/03/12 19:27:02  hmasui
00044  * Reduce bin size for event number histograms, and allow automatic bin extention
00045  *
00046  * Revision 1.8  2010/02/19 18:06:39  hmasui
00047  * Change the vertex range to +/-200 cm for vz histograms
00048  *
00049  * Revision 1.7  2010/02/16 02:13:34  hmasui
00050  * Add parent-parent geant id in the histogram name
00051  *
00052  * Revision 1.6  2010/02/12 16:24:13  hmasui
00053  * Extend the range of vz to +/-150cm for vx(vy) vs vz histograms
00054  *
00055  * Revision 1.5  2010/02/01 21:28:14  hmasui
00056  * Fix bugs for the binning of delta vx, vy, vz histograms
00057  *
00058  * Revision 1.4  2010/01/28 21:50:35  hmasui
00059  * Add Vx vs Vz and Vy vs Vz histograms.
00060  *
00061  * Revision 1.3  2010/01/26 18:07:33  hmasui
00062  * Change the binning for delta v_{x,y,z} from +/- 10 to +/- 1 cm
00063  *
00064  * Revision 1.2  2010/01/26 17:47:33  hmasui
00065  * Add histograms for eventid, runnumber and # of particles. Fix binning for delta pt vs pt
00066  *
00067  * Revision 1.1  2009/12/22 21:41:18  hmasui
00068  * Change class name from StEmbeddingQAMaker to StEmbeddingQA
00069  *
00070  ****************************************************************************************************/
00071 
00072 #include <algorithm>
00073 #include <fstream>
00074 #include <string>
00075 
00076 #include "TClonesArray.h"
00077 #include "TError.h"
00078 #include "TH1.h"
00079 #include "TH2.h"
00080 #include "TH3.h"
00081 #include "TFile.h"
00082 #include "TTree.h"
00083 #include "TObjArray.h"
00084 
00085 #include "StMessMgr.h"
00086 #include "StMiniMcEvent/StMiniMcEvent.h"
00087 #include "StMiniMcEvent/StTinyMcTrack.h"
00088 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
00089 #include "StMuDSTMaker/COMMON/StMuEvent.h"
00090 #include "StMuDSTMaker/COMMON/StMuTrack.h"
00091 #include "StParticleDefinition.hh"
00092 
00093 #include "StEmbeddingQA.h"
00094 #include "StEmbeddingQATrack.h"
00095 
00096 using namespace std ;
00097 
00098 ClassImp(StEmbeddingQA)
00099 
00100 //__________________________________________________________________________________________
00101 StEmbeddingQA::StEmbeddingQA()
00102   : mYear(2007), mProduction("P08ic"), mIsSimulation(kTRUE)
00103 {
00106 
00107   init() ;
00108 }
00109 
00110 //__________________________________________________________________________________________
00111 StEmbeddingQA::StEmbeddingQA(const Int_t year, const TString production, const Bool_t isSimulation)
00112   : mYear(year), mProduction(production), mIsSimulation(isSimulation)
00113 {
00115 
00116   init() ;
00117 }
00118 
00119 //__________________________________________________________________________________________
00120 StEmbeddingQA::~StEmbeddingQA()
00121 {
00122   clear();
00123 }
00124 
00125 //__________________________________________________________________________________________
00126 void StEmbeddingQA::clear()
00127 {
00129   LOG_DEBUG << "StEmbeddingQA::clear()" << endm;
00130 
00132   if ( mhRef ) delete mhRef ;
00133   if ( mhRefAccepted ) delete mhRefAccepted ;
00134   if ( mhVz ) delete mhVz ;
00135   if ( mhVzAccepted ) delete mhVzAccepted ;
00136   if ( mhVyVx ) delete mhVyVx ;
00137   if ( mhVxVz ) delete mhVxVz ;
00138   if ( mhVyVz ) delete mhVyVz ;
00139   if ( mhdVx ) delete mhdVx ;
00140   if ( mhdVy ) delete mhdVy ;
00141   if ( mhdVz ) delete mhdVz ;
00142   if ( mhEventId ) delete mhEventId ;
00143   if ( mhRunNumber ) delete mhRunNumber ;
00144 
00145   mGeantIdCollection.clear();
00146 
00148   for(UInt_t ic=0; ic<StEmbeddingQAConst::mNCategory; ic++){
00149     if ( mhGeantId[ic] ) delete mhGeantId[ic] ;
00150     if ( mhNParticles[ic] ) delete mhNParticles[ic] ;
00151 
00152     mGeantId[ic].clear();
00153     mhNHit[ic].clear();
00154     mhNCommonHitVsNHit[ic].clear();
00155     mhDca[ic].clear();
00156     mhPtVsEta[ic].clear();
00157     mhPtVsY[ic].clear();
00158     mhPtVsPhi[ic].clear();
00159     mhPtVsMom[ic].clear();
00160     mhdPtVsPt[ic].clear();
00161     mhdInvPtVsPt[ic].clear();
00162     mhMomVsEta[ic].clear();
00163     mhdEdxVsMomMc[ic].clear();
00164     mhdEdxVsMomMcPidCut[ic].clear();
00165     mhdEdxVsMomReco[ic].clear();
00166     mhdEdxVsMomRecoPidCut[ic].clear();
00167     mhRecoPVsMcP[ic].clear();
00168     mhNCommonHitVsNHit[ic].clear();
00169     mhEtaVsPhi[ic].clear();
00170     mhEtaVsVz[ic].clear();
00171     mhYVsVz[ic].clear();
00172   }
00173 }
00174 
00175 //__________________________________________________________________________________________
00176 void StEmbeddingQA::init()
00177 {
00179 
00180   LOG_INFO << endm;
00181   LOG_INFO << Form("   StEmbeddingQA::init() for year       : %10d", mYear) << endm;
00182   LOG_INFO << Form("   StEmbeddingQA::init() for production : %10s", mProduction.Data()) << endm;
00183   LOG_INFO << endm;
00184 
00185   mMuDstMaker = 0;
00186 
00187   mOutput = 0;
00188   mVz = -9999. ;
00189 
00190   mhRef = 0 ;
00191   mhRefAccepted = 0 ;
00192   mhVz = 0 ;
00193   mhVzAccepted = 0 ;
00194   mhVyVx = 0 ;
00195   mhVxVz = 0 ;
00196   mhVyVz = 0 ;
00197   mhdVx = 0 ;
00198   mhdVy = 0 ;
00199   mhdVz = 0 ;
00200   mhEventId = 0 ;
00201   mhRunNumber = 0 ;
00202 
00203   for(UInt_t ic=0; ic<StEmbeddingQAConst::mNCategory; ic++){
00204     mhNParticles[ic] = 0 ;
00205     mhGeantId[ic] = 0 ;
00206   }
00207 
00209   clear();
00210 }
00211 
00212 //__________________________________________________________________________________________
00213 void StEmbeddingQA::setRefMultMinCut(const Int_t refMultMin)
00214 {
00215   StEmbeddingQAUtilities::instance()->setRefMultMinCut(refMultMin) ;
00216 }
00217 
00218 //__________________________________________________________________________________________
00219 void StEmbeddingQA::setRefMultMaxCut(const Int_t refMultMax)
00220 {
00221   StEmbeddingQAUtilities::instance()->setRefMultMaxCut(refMultMax) ;
00222 }
00223 
00224 //__________________________________________________________________________________________
00225 void StEmbeddingQA::setZVertexCut(const Float_t vz)
00226 {
00227   StEmbeddingQAUtilities::instance()->setZVertexCut(vz) ;
00228 }
00229 
00230 //__________________________________________________________________________________________
00231 void StEmbeddingQA::addTriggerIdCut(const UInt_t id)
00232 {
00234 
00235   StEmbeddingQAUtilities::instance()->addTriggerIdCut(id) ;
00236 }
00237 
00238 //__________________________________________________________________________________________
00239 void StEmbeddingQA::setRapidityCut(const Float_t ycut)
00240 {
00242   StEmbeddingQAUtilities::instance()->setRapidityCut(ycut) ;
00243 }
00244 
00245 //__________________________________________________________________________________________
00246 Bool_t StEmbeddingQA::isRefMultOk(const StMiniMcEvent& mcevent) const
00247 {
00249   return StEmbeddingQAUtilities::instance()->isRefMultOk(mcevent.nUncorrectedPrimaries()) ;
00250 }
00251 
00252 //__________________________________________________________________________________________
00253 Bool_t StEmbeddingQA::isZVertexOk(const StMiniMcEvent& mcevent) const
00254 {
00256   return StEmbeddingQAUtilities::instance()->isZVertexOk(mcevent.vertexZ()) ;
00257 }
00258 
00259 //__________________________________________________________________________________________
00260 Bool_t StEmbeddingQA::isTriggerOk(StMuEvent* event) const
00261 {
00263   const vector<UInt_t> triggerId(StEmbeddingQAUtilities::instance()->getTriggerIdCut());
00264   if( triggerId.empty() ) return kTRUE ;
00265 
00267   for(UInt_t i=0; i<triggerId.size(); i++){
00268     if( event->triggerIdCollection().nominal().isTrigger( triggerId[i] ) ){
00269       LOG_DEBUG << "StEmbeddingQA::isTriggerOk  Trigger found: " << triggerId[i] << endm ;
00270       return kTRUE ;
00271     }
00272   }
00273 
00274   return kFALSE ;
00275 }
00276 
00277 //__________________________________________________________________________________________
00278 Bool_t StEmbeddingQA::book(const TString outputFileName)
00279 {
00281   LOG_DEBUG << "StEmbeddingQA::book()" << endm;
00282 
00283   TString fileName(outputFileName);
00284 
00288   if( fileName.IsWhitespace() ){
00289     const TString data = (mIsSimulation) ? "embedding" : "real" ;
00290     fileName = Form("qa_%s_%d_%s.root", data.Data(), mYear, mProduction.Data());
00291   }
00292 
00294   mOutput = TFile::Open(fileName, "recreate");
00295   mOutput->cd();
00296   LOG_INFO << "    OPEN " << mOutput->GetName() << endm;
00297 
00303   StEmbeddingQAUtilities* utility = StEmbeddingQAUtilities::instance() ;
00304 
00305   // Event-wise informations
00306   mhRef         = new TH1D("hRef", "refMult", 1000, 0, 1000);
00307   mhRefAccepted = new TH1D("hRefAccepted", "refMult", 1000, 0, 1000);
00308   mhVz         = new TH1D("hVz", "z-vertex", 400, -200, 200);
00309   mhVzAccepted = new TH1D("hVzAccepted", "z-vertex with z-vertex cut", 400, -200, 200);
00310   mhRef->SetXTitle("refMult");
00311   mhRefAccepted->SetXTitle("refMult");
00312   mhVz->SetXTitle("v_{z} (cm)");
00313   mhVzAccepted->SetXTitle("v_{z} (cm)");
00314 
00315   utility->setStyle(mhRef);
00316   utility->setStyle(mhRefAccepted);
00317   utility->setStyle(mhVz);
00318   utility->setStyle(mhVzAccepted);
00319 
00320   // Bins for vertices
00321   const Int_t vtxBin     = 800 ;
00322   const Double_t vtxMin  = -2.0 ;
00323   const Double_t vtxMax  =  2.0 ;
00324   const Double_t binSize = (vtxMax-vtxMin)/static_cast<Double_t>(vtxBin);
00325 
00326   // Shift vertices range by binSize/2 to locate v=0 at the middle of 0-th bin
00327   const Double_t vtxMinShift = vtxMin - binSize/2.0 ;
00328   const Double_t vtxMaxShift = vtxMax - binSize/2.0 ;
00329 
00330   mhVyVx = new TH2D("hVyVx", "v_{y} vs v_{x}", vtxBin, vtxMinShift, vtxMaxShift, vtxBin, vtxMinShift, vtxMaxShift);
00331   mhVxVz = new TH2D("hVxVz", "v_{x} vs v_{z}", 300, -150, 150, vtxBin, vtxMinShift, vtxMaxShift);
00332   mhVyVz = new TH2D("hVyVz", "v_{y} vs v_{z}", 300, -150, 150, vtxBin, vtxMinShift, vtxMaxShift);
00333   mhVyVx->SetXTitle("v_{x} (cm)"); mhVyVx->SetYTitle("v_{y} (cm)");
00334   mhVxVz->SetXTitle("v_{z} (cm)"); mhVxVz->SetYTitle("v_{x} (cm)");
00335   mhVyVz->SetXTitle("v_{z} (cm)"); mhVyVz->SetYTitle("v_{y} (cm)");
00336 
00337   utility->setStyle(mhVyVx);
00338   utility->setStyle(mhVxVz);
00339   utility->setStyle(mhVyVz);
00340 
00341   mhdVx = new TH1D("hdVx", "#Delta x = v_{x} - v_{x}(MC)",  vtxBin, vtxMinShift, vtxMaxShift);
00342   mhdVy = new TH1D("hdVy", "#Delta y = v_{y} - v_{y}(MC)",  vtxBin, vtxMinShift, vtxMaxShift);
00343   mhdVz = new TH1D("hdVz", "#Delta z = v_{z} - v_{z}(MC)",  vtxBin, vtxMinShift, vtxMaxShift);
00344   mhdVx->SetXTitle("#Deltav_{x} = v_{x} - v_{x}(MC) (cm)");
00345   mhdVy->SetXTitle("#Deltav_{y} = v_{y} - v_{y}(MC) (cm)");
00346   mhdVz->SetXTitle("#Deltav_{z} = v_{z} - v_{z}(MC) (cm)");
00347 
00348   utility->setStyle(mhdVx);
00349   utility->setStyle(mhdVy);
00350   utility->setStyle(mhdVz);
00351 
00352   mhEventId   = new TH1D("hEventId", "Event id", 1000, 0, 1000);
00353   mhRunNumber = new TH1D("hRunNumber", "Run id - (Year - 1999)#times10^{6}", 400000, 0, 400000);
00354   mhEventId->SetXTitle("Event id");
00355   mhRunNumber->SetXTitle("Run number");
00356 
00357   // Set bit to automatic bin extention
00358   mhEventId->SetBit(TH1::kCanRebin);
00359 
00360   utility->setStyle( mhEventId   ) ;
00361   utility->setStyle( mhRunNumber ) ;
00362 
00363   for(UInt_t ic=0; ic<StEmbeddingQAConst::mNCategory; ic++){
00364     mhNParticles[ic] = new TH1D(Form("hNParticles_%d", ic), 
00365         Form("Number of particles per event, %s", utility->getCategoryTitle(ic).Data()), 1000, 0, 1000);
00366     mhNParticles[ic]->SetXTitle("# of particles / event");
00367 
00368     utility->setStyle(mhNParticles[ic]);
00369 
00370     // Initialize geantid histogram. Increase the bin and maximum in order to cover id > 10k
00371     mhGeantId[ic] = new TH1D(Form("hGeantId_%d", ic), Form("Geantid, %s", utility->getCategoryTitle(ic).Data()), 100000, 0, 100000) ;
00372     mhGeantId[ic]->SetXTitle("Geantid");
00373 
00374     utility->setStyle(mhGeantId[ic]);
00375   }
00376 
00377   mOutput->GetList()->Sort();
00378   LOG_INFO << endm << endm << endm;
00379 
00380   return kStOk ;
00381 }
00382 
00383 //__________________________________________________________________________________________
00384 Bool_t StEmbeddingQA::make(const TString inputFileName, const Bool_t isSimulation)
00385 {
00387 
00389   if(!mOutput || !mOutput->IsOpen()){
00390     Error("StEmbeddingQA::Make", "Output file is not opened");
00391     return kStErr ;
00392   }
00393 
00394   if( isSimulation ){
00396     LOG_INFO << "------------------------------------------------------------------------------------" << endm;
00397     LOG_INFO << "            Fill embedding ..." << endm;
00398     LOG_INFO << "------------------------------------------------------------------------------------" << endm;
00399     fillEmbedding(inputFileName);
00400   }
00401   else{
00403     LOG_INFO << "------------------------------------------------------------------------------------" << endm;
00404     LOG_INFO << "            Fill real data ..." << endm;
00405     LOG_INFO << "------------------------------------------------------------------------------------" << endm;
00406     fillRealData(inputFileName);
00407   }
00408 
00409   return kStOk ;
00410 }
00411 
00412 //__________________________________________________________________________________________
00413 Bool_t StEmbeddingQA::fillEmbedding(const TString inputFileName)
00414 {
00416 
00418   TFile* file = TFile::Open(inputFileName);
00419   if( !file || !file->IsOpen() ){
00420     Error("StEmbeddingQA::fillEmbedding", "can't open %s", inputFileName.Data());
00421     return kStErr ;
00422   }
00423 
00425   const TString treeName("StMiniMcTree");
00426   TTree* tree = (TTree*) file->Get(treeName);
00427   if( !tree ){
00428     Error("StEmbeddingQA::fillEmbedding", "can't find %s tree", treeName.Data());
00429     return kStErr ;
00430   }
00431 
00433   StMiniMcEvent* mMiniMcEvent = new StMiniMcEvent();
00434   TBranch* b_MiniMcEvent = tree->GetBranch("StMiniMcEvent");
00435   b_MiniMcEvent->SetAddress(&mMiniMcEvent);
00436 
00437   // Tracks in minimc tree
00438   //     mMcTracks     = new TClonesArray("StTinyMcTrack",nMcTrack);
00439   //     mMatchedPairs = new TClonesArray("StMiniMcPair",nMatchedPair);
00440   //     mMergedPairs  = new TClonesArray("StMiniMcPair",nMergedPair);
00441   //     mSplitPairs   = new TClonesArray("StMiniMcPair",nSplitPair);
00442   //     
00443   //     mGhostPairs   = new TClonesArray("StMiniMcPair",nGhostPair);
00444   //     mContamPairs  = new TClonesArray("StContamPair",nContamPair); 
00445   //     mMatGlobPairs = new TClonesArray("StMiniMcPair",nMatGlobPair);
00446   //
00447   //  Track name definition in StMiniMcEvent.h
00448   //   enum Category { MC,MATCHED,MERGED,SPLIT,CONTAM,GHOST,MATGLOB};
00449 
00450   const Long64_t nentries = tree->GetEntries();
00451   LOG_INFO << Form("    OPEN %s,  # of event = %10d", inputFileName.Data(), nentries) << endm;
00452 
00454   for (Int_t ievent=0; ievent<nentries;ievent++) {
00455     tree->GetEntry(ievent);
00456 
00457     const Float_t vx   = mMiniMcEvent->vertexX() ;
00458     const Float_t vy   = mMiniMcEvent->vertexY() ;
00459     const Float_t vz   = mMiniMcEvent->vertexZ() ;
00460     const Float_t vxmc = mMiniMcEvent->mcVertexX() ;
00461     const Float_t vymc = mMiniMcEvent->mcVertexY() ;
00462     const Float_t vzmc = mMiniMcEvent->mcVertexZ() ;
00463     const Int_t nprimaries = mMiniMcEvent->nUncorrectedPrimaries(); 
00464 
00465     mhRef->Fill(nprimaries);
00466     mhVz->Fill(vz);
00467 
00468     mVz = vz ;
00469 
00471     if( !isZVertexOk(*mMiniMcEvent) ) continue ;
00472 
00474     if( !isRefMultOk(*mMiniMcEvent) ) continue ;
00475 
00477     mhRefAccepted->Fill(nprimaries);
00478     mhVzAccepted->Fill(vz);
00479     mhVyVx->Fill(vx, vy);
00480     mhVxVz->Fill(vz, vx);
00481     mhVyVz->Fill(vz, vy);
00482     mhdVx->Fill( vx - vxmc );
00483     mhdVy->Fill( vy - vymc );
00484     mhdVz->Fill( vz - vzmc );
00485 
00487     mhEventId->Fill( mMiniMcEvent->eventId() );
00488     mhRunNumber->Fill( StEmbeddingQAUtilities::instance()->getRunNumber(mMiniMcEvent->runId(), mYear) );
00489 
00491     for(UInt_t categoryid=0; categoryid<StEmbeddingQAConst::mNEmbedding; categoryid++){
00492       const Int_t nTrack = getNtrack(categoryid, *mMiniMcEvent) ;
00493       mhNParticles[categoryid]->Fill(nTrack);
00494 
00495       const Int_t nevents = (Int_t)mhVz->GetEntries();
00496       if( nevents % 100 == 0 && nTrack > 0 ){
00497         LOG_INFO << Form("####  accept/throw=%10d/%10d, category=%10s, ntrack=%10d",
00498             (Int_t)mhVzAccepted->GetEntries(), nevents, 
00499             StEmbeddingQAUtilities::instance()->getCategoryName(categoryid).Data(), nTrack)
00500           << endm;
00501       }
00502 
00504       for(Int_t itrk=0; itrk<nTrack; itrk++){
00505         fillEmbeddingTracks(*mMiniMcEvent, categoryid, itrk) ;
00506       } // Track loop
00507     }// Track category
00508 
00509   }// event loop
00510 
00512   delete mMiniMcEvent ;
00513 
00515   file->Close() ;
00516 
00517   return kStOk ;
00518 }
00519 
00520 //__________________________________________________________________________________________
00521 Bool_t StEmbeddingQA::fillRealData(const TString inputFileName)
00522 {
00524 
00526   Int_t ievent = 0;
00527   Int_t ieventAccept = 0;
00528   while( !mMuDstMaker->Make() ){ //read event
00529     ievent++;
00530     StMuEvent* muEvent = mMuDstMaker->muDst()->event() ;
00531     if(!muEvent) continue ;
00532 
00535     const Double_t vx = muEvent->primaryVertexPosition().x() ;
00536     const Double_t vy = muEvent->primaryVertexPosition().y() ;
00537     const Double_t vz = muEvent->primaryVertexPosition().z() ;
00538     const Int_t refMult = muEvent->refMult() ;
00539     mhRef->Fill(refMult);
00540 //    const Bool_t isVertexBad = TMath::Abs(vz) >= mVertexCut
00541     const Bool_t isVertexBad = !StEmbeddingQAUtilities::instance()->isZVertexOk(vz)
00542       || !StEmbeddingQAUtilities::instance()->isRefMultOk(refMult)
00543       || ( TMath::Abs(vx) < 1.0e-5 && TMath::Abs(vy) < 1.0e-5 && TMath::Abs(vz) < 1.0e-5 )
00544       || ( TMath::Abs(vx) > 1000 || TMath::Abs(vy) > 1000 || TMath::Abs(vz) > 1000 )
00545       ;
00546     if( isVertexBad ) continue ;
00547 
00549     if( !isTriggerOk(muEvent) ) continue ;
00550 
00551     mhRefAccepted->Fill(refMult);
00552     mhVz->Fill(vz);
00553     mhVyVx->Fill(vx, vy);
00554 
00555     mVz = vz ;
00556 
00558     for(UInt_t ic=0; ic<StEmbeddingQAConst::mNReal; ic++){
00559       const Int_t categoryid = ic + StEmbeddingQAConst::mNEmbedding ;
00560 
00561       if( ieventAccept % 100 == 0 ){
00562         LOG_INFO << Form("%85s ####  accept/throw=%10d/%10d, category=%10s",
00563             mMuDstMaker->GetFileName(), ieventAccept, ievent, StEmbeddingQAUtilities::instance()->getCategoryName(categoryid).Data())
00564           << endm;
00565       }
00566 
00567       const TObjArray* tracks = (ic==0) 
00568         ?  mMuDstMaker->muDst()->primaryTracks()   // Primary tracks
00569         :  mMuDstMaker->muDst()->globalTracks() ;  // Global tracks
00570 
00571       TObjArrayIter trackIterator(tracks);
00572       StMuTrack* track = 0;
00573       Int_t itrk = 0;
00574       while ( ( track = (StMuTrack*) trackIterator.Next() ) ){
00575         fillRealTracks(*track, categoryid, itrk);
00576         itrk++;
00577       }
00578     }
00579 
00580     ieventAccept++;
00581   }// event loop
00582 
00583   LOG_INFO << "End of real data" << endm;
00584   LOG_INFO << "Total # of events = " << ievent << endm;
00585 
00586   return kStOk ;
00587 }
00588 
00589 //__________________________________________________________________________________________
00590 Bool_t StEmbeddingQA::runRealData(const TString inputFileList)
00591 {
00593 
00595   if( mMuDstMaker ) delete mMuDstMaker ;
00596 
00597   LOG_INFO << "###    Read " << inputFileList << endm;
00598   mMuDstMaker  = new StMuDstMaker(0, 0, "", inputFileList, "MuDst", 1000);
00599 
00601   LOG_INFO << "StEmbeddingQA::runRealData()" << endm;
00602   LOG_INFO << "Add electrons, pions, kaons and protons for the real data QA" << endm;
00603   for(UInt_t ic=0; ic<StEmbeddingQAConst::mNReal; ic++){
00604     const Int_t categoryid = ic + StEmbeddingQAConst::mNEmbedding ;
00605 
00606     const Int_t parentid = 0 ;       // real tracks are assumed to be primary
00607     const Int_t parentparentid = 0 ; // real tracks are assumed to be primary
00608     const Int_t geantprocess = 0 ;
00609     expandHistograms(categoryid, 2, parentid, parentparentid, geantprocess);
00610     expandHistograms(categoryid, 3, parentid, parentparentid, geantprocess);
00611     expandHistograms(categoryid, 8, parentid, parentparentid, geantprocess);
00612     expandHistograms(categoryid, 9, parentid, parentparentid, geantprocess);
00613     expandHistograms(categoryid, 11, parentid, parentparentid, geantprocess);
00614     expandHistograms(categoryid, 12, parentid, parentparentid, geantprocess);
00615     expandHistograms(categoryid, 14, parentid, parentparentid, geantprocess);
00616     expandHistograms(categoryid, 15, parentid, parentparentid, geantprocess);
00617 
00618     // Make sure the input particle list
00619     for(vector<Int_t>::iterator iter = mGeantId[categoryid].begin(); iter != mGeantId[categoryid].end(); iter++){
00620       const Int_t geantid = (*iter) ;
00621       LOG_DEBUG << Form("  Input geant id = %10d,  name = %10s", geantid,
00622           StEmbeddingQAUtilities::instance()->getParticleDefinition(geantid)->name().c_str()) << endm ;
00623     }
00624   }
00625 
00626   make(inputFileList, kFALSE);
00627 
00628   return kStOk ;
00629 }
00630 
00631 //__________________________________________________________________________________________
00632 Bool_t StEmbeddingQA::runEmbedding(const TString inputFileList)
00633 {
00635 
00636   ifstream fEmbedding(inputFileList);
00637   if(!fEmbedding){
00638     Error("StEmbeddingQA::runEmbedding", "can't find %s", inputFileList.Data());
00639     return kFALSE ;
00640   }
00641   LOG_INFO << "###    Read " << inputFileList << endm;
00642 
00644   TString file ;
00645   while( fEmbedding >> file ){
00646     LOG_INFO << "####     Read file : " << file << endm;
00647     make(file, kTRUE);
00648   }
00649 
00650   return kStOk ;
00651 }
00652 
00653 //__________________________________________________________________________________________
00654 Bool_t StEmbeddingQA::run(const TString inputFileList)
00655 {
00657 
00658   LOG_INFO << "StEmbeddingQA::run()" << endm ;
00659 //  LOG_INFO << "  z-vertex cut is |vz| < " << mVertexCut << endm;
00660   StEmbeddingQAUtilities::instance()->PrintCuts() ;
00661 
00662   if( mIsSimulation ){
00663     // Embedding QA
00664     return runEmbedding(inputFileList) ;
00665   }
00666   else{
00667     // Real data QA
00668     return runRealData(inputFileList) ;
00669   }
00670 
00671   return kFALSE ;
00672 }
00673 
00674 //__________________________________________________________________________________________
00675 Bool_t StEmbeddingQA::end() const
00676 {
00678 
00679   LOG_INFO << "    End of StEmbeddingQA" << endm;
00680   LOG_INFO << "    Write output " << mOutput->GetName() << endm;
00681   mOutput->GetList()->Sort();
00682 
00683   mOutput->Write();
00684   mOutput->Close();
00685 
00686   return kStOk ;
00687 }
00688 
00689 //__________________________________________________________________________________________
00690 void StEmbeddingQA::fillEmbeddingTracks(const StMiniMcEvent& mcevent, const Int_t categoryid, const Int_t itrk)
00691 {
00693 
00695   StEmbeddingQATrack* miniTrack = getEmbeddingQATrack(mcevent, categoryid, itrk);
00696 
00698   if(!miniTrack) return ;
00699 
00701   fillHistograms(*miniTrack, categoryid);
00702 
00704   delete miniTrack ;
00705 }
00706 
00707 //____________________________________________________________________________________________________
00708 StEmbeddingQATrack* StEmbeddingQA::getEmbeddingQATrack(const StMiniMcEvent& mcevent, const Int_t categoryid, const Int_t itrk)
00709 {
00711   StEmbeddingQAUtilities* utility = StEmbeddingQAUtilities::instance();
00712 
00713   const TString name(utility->getCategoryName(categoryid));
00714   const Category category(utility->getCategory(categoryid));
00715 
00716   if ( utility->isMc(name) ){
00718     const StTinyMcTrack* track = (StTinyMcTrack*) mcevent.tracks(category)->At(itrk) ;
00719 
00721     if ( track->isPrimary() != 1 ){
00722       LOG_DEBUG << Form("StEmbeddingQA::getEmbeddingQATrack", "MC track with GeantID %3d is not a primary track", track->geantId()) << endm ;
00723       return 0;
00724     }
00725 
00726     return (new StEmbeddingQATrack(name, *track));
00727   }
00728   else if ( utility->isMatched(name) || utility->isGhost(name) || utility->isMatchedGlobal(name) ){
00730     StMiniMcPair* track = (StMiniMcPair*) mcevent.tracks(category)->At(itrk) ;
00731 
00733     if ( !utility->isGeantIdOk(track->geantId()) ){
00734       LOG_DEBUG << Form("StEmbeddingQA::getEmbeddingQATrack", "No geantid = %3d exists in StParticleTable", track->geantId()) << endm ;
00735       return 0;
00736     }
00737 
00738     return (new StEmbeddingQATrack(name, track));
00739   }
00740   else if ( utility->isContaminated(name) ){
00742     StContamPair* track = (StContamPair*) mcevent.tracks(category)->At(itrk) ;
00743 
00745     if ( !utility->isGeantIdOk(track->geantId()) ){
00746       LOG_DEBUG << Form("StEmbeddingQA::getEmbeddingQATrack", "No geantid = %3d exists in StParticleTable", track->geantId()) << endm ;
00747       return 0;
00748     }
00749 
00752     const Bool_t isGeantProcessOk = ( track->mGeantProcess == 5 || (track->parentGeantId() == 1 && track->mGeantProcess == 6) );
00753 
00754     if ( !isGeantProcessOk ){
00755       LOG_DEBUG << Form("StEmbeddingQA::getEmbeddingTrack()  geantprocess = %3d. Skip the track", track->mGeantProcess) << endm;
00756       return 0;
00757     }
00758 
00759     return (new StEmbeddingQATrack(name, track));
00760   }
00761   else{
00762     Warning("StEmbeddingQA::fillEmbeddingTracks", "Unknown category id, id=%3d", categoryid);
00763     return 0;
00764   }
00765 
00766   return 0;
00767 }
00768 
00769 //__________________________________________________________________________________________
00770 void StEmbeddingQA::fillRealTracks(const StMuTrack& track, const Int_t categoryid, const Int_t itrk)
00771 {
00773   for(vector<Int_t>::iterator iter = mGeantId[categoryid].begin(); iter != mGeantId[categoryid].end(); iter++){
00774     const Int_t geantid = (*iter) ;
00775     StEmbeddingQATrack miniTrack(StEmbeddingQAUtilities::instance()->getCategoryName(categoryid), track, geantid);
00776  
00777     fillHistograms(miniTrack, categoryid);
00778   }
00779 }
00780 
00781 //__________________________________________________________________________________________
00782 void StEmbeddingQA::fillHistograms(const StEmbeddingQATrack& track, const Int_t categoryid)
00783 {
00785 
00788 
00790   const Int_t geantid = track.getGeantId() ;
00791   if ( geantid < 0 ) return ;
00792 
00794   if( !track.isPtAndEtaOk() ) return ;
00795 
00797   StEmbeddingQAUtilities* utility = StEmbeddingQAUtilities::instance() ;
00798   const Double_t pt    = (utility->isReal(track.getName())) ? track.getPtRc()       : track.getPtMc() ;
00799   const Double_t mom   = (utility->isReal(track.getName())) ? track.getPRc()        : track.getPMc() ;
00800   const Double_t eta   = (utility->isReal(track.getName())) ? track.getEtaRc()      : track.getEtaMc() ;
00801   const Double_t y     = (utility->isReal(track.getName())) ? track.getRapidityRc() : track.getRapidityMc() ;
00802 
00803   // Reconstructed momentum
00804   const Double_t momRc = track.getPRc() ;
00805 
00809   expandHistograms(categoryid, geantid, track.getParentGeantId(), track.getParentParentGeantId(),
00810       track.getGeantProcess());
00811 
00812   // Fill geant id
00813   mhGeantId[categoryid]->Fill(track.getGeantId());
00814 
00815   // Rapidity cut, default rapidity cut is 10 (Use setRapidityCut(const Float_t ycut) to restrict y window)
00816 //  if(!track.isRapidityOk(mRapidityCut)) return ;
00817   if(!utility->isRapidityOk(y)) return ;
00818 
00819   const TString idcollection(getIdCollection(geantid, track.getParentGeantId(), track.getParentParentGeantId()));
00820 
00821   // dE/dx (no PID cut)
00822   //  - Add NHit cut (Nov/13/2009)
00823   //  - Add NHitFit/NHitPoss cut
00824   if( track.isDcaOk() && track.isNHitOk() && track.isNHitToNPossOk() ) {
00825     mhdEdxVsMomMc[categoryid][idcollection]->Fill(mom, track.getdEdxkeV());
00826     mhdEdxVsMomReco[categoryid][idcollection]->Fill(momRc, track.getdEdxkeV());
00827   }
00828 
00829   //  Oct/21/2009
00833 
00834   if ( !track.isNSigmaOk(geantid) ) return ;
00835 
00836   if( track.isDcaOk() ){
00837     // Fill NHit points
00838     //  Added common hit cuts
00839     if ( track.isCommonHitOk() ) {
00840       mhNHit[categoryid][idcollection]->Fill(pt, eta, track.getNHit());
00841     }
00842 
00843     if( track.isNHitOk() ){
00844       const Double_t phi = track.getPhi() ;
00845 
00846       // Fill Ncommon hits vs Nhits
00847       mhNCommonHitVsNHit[categoryid][idcollection]->Fill(pt, track.getNHit(), track.getNCommonHit());
00848 
00849       // NHitFit/NHitPoss cut
00850       if( track.isNHitToNPossOk() ) {
00851         // dE/dx (with PID cut)
00852         mhdEdxVsMomMcPidCut[categoryid][idcollection]->Fill(mom, track.getdEdxkeV());
00853         mhdEdxVsMomRecoPidCut[categoryid][idcollection]->Fill(momRc, track.getdEdxkeV());
00854   
00855         // Correlation between reconstructed and MC momentum
00856         mhRecoPVsMcP[categoryid][idcollection]->Fill(mom, momRc);
00857   
00858         // Pt, eta, phi
00859         mhPtVsEta[categoryid][idcollection]->Fill(eta, pt);
00860         mhPtVsY[categoryid][idcollection]->Fill(y, pt);
00861         mhPtVsPhi[categoryid][idcollection]->Fill(phi, pt);
00862         mhPtVsMom[categoryid][idcollection]->Fill(mom, pt);
00863         mhdPtVsPt[categoryid][idcollection]->Fill(pt, track.getPtRc()-pt);
00864         mhdInvPtVsPt[categoryid][idcollection]->Fill(pt, (1.0/track.getVectorGl().perp()-1.0/pt)*1000.);
00865         mhMomVsEta[categoryid][idcollection]->Fill(eta, mom);
00866   
00867         mhEtaVsPhi[categoryid][idcollection]->Fill(phi, eta);
00868         mhEtaVsVz[categoryid][idcollection]->Fill(mVz, eta);
00869         mhYVsVz[categoryid][idcollection]->Fill(mVz, y);
00870       }
00871     }
00872   }
00873 
00874   // Fill Dca
00875   if( track.isNHitOk() && track.isNHitToNPossOk() ){
00876     mhDca[categoryid][idcollection]->Fill(pt, eta, track.getDcaGl());
00877   }
00878 
00879   LOG_DEBUG << Form("     RC:(nfit, pt, eta, phi) = (%5d, %1.4f, %1.4f, %1.4f)  MC:(pt, eta) = (%1.4f, %1.4f)",
00880       track.getNHit(), track.getPtRc(), track.getEtaRc(), track.getPhi(), track.getPtMc(), track.getEtaMc()
00881       ) << endm;
00882 }
00883 
00884 //__________________________________________________________________________________________
00885 Bool_t StEmbeddingQA::pushBackGeantId(const Int_t categoryid, const Int_t geantid, const Int_t parentid,
00886     const Int_t parentparentid, const Int_t geantprocess)
00887 {
00889   StEmbeddingQAUtilities* utility = StEmbeddingQAUtilities::instance() ;
00890 
00891   const Bool_t isContaminated = utility->isContaminated(utility->getCategoryName(categoryid)) ;
00892 
00893   Bool_t isOk = kFALSE ;
00894   if ( mGeantId[categoryid].empty() ){
00896     LOG_INFO << "#----------------------------------------------------------------------------------------------------" << endm ;
00897     LOG_INFO << Form("StEmbeddingQA::pushBackGeantId()   mGeantId[%d] is empty", categoryid) << endm;
00898     if( isContaminated ){
00899       LOG_INFO << Form("StEmbeddingQA::pushBackGeantId()   Find a new geant id,  (geant, parent, parent-parent, process) = (%4d, %4d, %4d, %4d) (%10s)",
00900           geantid, parentid, parentparentid, geantprocess, utility->getCategoryName(categoryid).Data()) << endm;
00901     }
00902     else{
00903       LOG_INFO << Form("StEmbeddingQA::pushBackGeantId()   Find a new geant id,  geantid = %5d (%10s)", 
00904           geantid, utility->getCategoryName(categoryid).Data()) << endm;
00905     }
00906     LOG_INFO << "#----------------------------------------------------------------------------------------------------" << endm ;
00907     mGeantId[categoryid].push_back(geantid);
00908 
00909     isOk = kTRUE ;
00910   }
00911   else{
00915     const vector<Int_t>::iterator iter = find(mGeantId[categoryid].begin(), mGeantId[categoryid].end(), geantid) ;
00916  
00917     if ( iter != mGeantId[categoryid].end() ){
00919       isOk = kFALSE;
00920     }
00921     else{
00923       if( isContaminated ){
00924         LOG_INFO << Form("StEmbeddingQA::pushBackGeantId()   Find a new geant id,  (geant, parent, parent-parent, process) = (%4d, %4d, %4d, %4d) (%10s)",
00925             geantid, parentid, parentparentid, geantprocess, utility->getCategoryName(categoryid).Data()) << endm;
00926       }
00927       else{
00928         LOG_INFO << Form("StEmbeddingQA::pushBackGeantId()   Find a new geant id,  geantid = %5d (%10s)", 
00929             geantid, utility->getCategoryName(categoryid).Data()) << endm;
00930       }
00931       mGeantId[categoryid].push_back(geantid);
00932 
00933       isOk = kTRUE ;
00934     }
00935   }
00936 
00938   if( !utility->isContaminated(utility->getCategoryName(categoryid)) ) return isOk ;
00939 
00943 
00945   isOk = kFALSE ;
00946   if( parentid == 0 ){
00947     isOk = kFALSE ;
00948   }
00949   else{
00950     const TString idcollection(getIdCollection(geantid, parentid, parentparentid));
00951 
00952     if( mGeantIdCollection.empty() ){
00954       LOG_INFO << "#----------------------------------------------------------------------------------------------------" << endm ;
00955       LOG_INFO << "StEmbeddingQA::pushBackGeantId()  mGeantIdCollection is empty." << endm;
00956       LOG_INFO << Form("StEmbeddingQA::pushBackGeantId()   Push back             (geant, parent, parent-parent) = (%4d, %4d, %4d) (process=%4d)", 
00957           geantid, parentid, parentparentid, geantprocess) << endm;
00958       LOG_INFO << "#----------------------------------------------------------------------------------------------------" << endm ;
00959 
00960       mGeantIdCollection.push_back( idcollection );
00961       return kTRUE ;
00962     }
00963 
00965     const vector<TString>::iterator iterIdCollection = find(mGeantIdCollection.begin(), mGeantIdCollection.end(), idcollection) ;
00966     if ( iterIdCollection != mGeantIdCollection.end() ){
00968       isOk = kFALSE;
00969     }
00970     else{
00972       LOG_INFO << Form("StEmbeddingQA::pushBackGeantId()   Push back             (geant, parent, parent-parent) = (%4d, %4d, %4d) (process=%4d)", 
00973           geantid, parentid, parentparentid, geantprocess) << endm;
00974       mGeantIdCollection.push_back( idcollection );
00975 
00976       isOk = kTRUE ;
00977     }
00978   }// parentid != 0
00979 
00980   return isOk ;
00981 }
00982 
00983 
00984 //__________________________________________________________________________________________
00985 void StEmbeddingQA::expandHistograms(const Int_t categoryid, const Int_t geantid, const Int_t parentid,
00986     const Int_t parentparentid, const Int_t geantprocess)
00987 {
00990   if ( !pushBackGeantId(categoryid, geantid, parentid, parentparentid, geantprocess) ) return ;
00991 
00992   // Expand histograms
00993   mOutput->cd();
00994 
00995   const Int_t ptBin    = 100 ;
00996   const Float_t ptMin  = 0.0 ;
00997   const Float_t ptMax  = 10.0 ;
00998   const Int_t etaBin   = 100 ;
00999   const Float_t etaMin = -2.5 ;
01000   const Float_t etaMax =  2.5 ;
01001 
01003   const StEmbeddingQAUtilities* utility = StEmbeddingQAUtilities::instance() ;
01004   const Char_t* categoryTitle(utility->getCategoryTitle(categoryid).Data());
01005 
01007   const TString nameSuffix = (parentid>0) ? Form("_%d_%d_%d_%d", categoryid, parentparentid, parentid, geantid)
01008     : Form("_%d_%d", categoryid, geantid);
01009 
01011   TString CategoryAndGeantId = (parentid>0) ?  Form("%s (%s #rightarrow %s)", categoryTitle, 
01012       utility->getParticleDefinition(parentid)->name().c_str(), utility->getParticleDefinition(geantid)->name().c_str())
01013     : Form("%s (%s)", categoryTitle, utility->getParticleDefinition(geantid)->name().c_str());
01014 
01015   if ( parentparentid > 0 ){
01016     CategoryAndGeantId = Form("%s (%s #rightarrow %s #rightarrow %s)", categoryTitle,
01017         utility->getParticleDefinition(parentparentid)->name().c_str(),
01018         utility->getParticleDefinition(parentid)->name().c_str(),
01019         utility->getParticleDefinition(geantid)->name().c_str());
01020   }
01021 
01022   const TString categoryName(utility->getCategoryName(categoryid));
01023   const Bool_t isMc        = utility->isMc(categoryName);
01024   const Bool_t isEmbedding = utility->isEmbedding(categoryName);
01025 
01026   const TString idcollection(getIdCollection(geantid, parentid, parentparentid));
01027 
01028   // NHit vs eta vs MC pt
01029   TString title(Form("N_{fit} distribution (|dcaGl|<3cm), %s", CategoryAndGeantId.Data()));
01030   if( isMc ) title = Form("N_{fit} distribution, %s", CategoryAndGeantId.Data());
01031   TH3* hNhit = new TH3D(Form("hNHit%s", nameSuffix.Data()), title, 10, 0, 5, 10, -1.0, 1.0, 50, 0, 50) ;
01032   hNhit->SetXTitle("MC p_{T} (GeV/c)");
01033   hNhit->SetYTitle("#eta");
01034   hNhit->SetZTitle("N_{fit}");
01035   utility->setStyle(hNhit);
01036   mhNHit[categoryid].insert( std::make_pair(idcollection, hNhit) );
01037 
01038   // Ncommon hit vs Nfit
01039   title = Form("N_{common} hit vs N_{fit} (|dcaGl|<3cm), %s", CategoryAndGeantId.Data());
01040   if( isMc ) title = Form("N_{common} hit vs N_{fit}, %s", CategoryAndGeantId.Data());
01041   TH3* hNCommonHitVsNHit = new TH3D(Form("hNCommonHitVsNHit%s", nameSuffix.Data()), title, 10, 0, 5, 50, 0, 50, 50, 0, 50) ;
01042   hNCommonHitVsNHit->SetXTitle("p_{T} (GeV/c)");
01043   hNCommonHitVsNHit->SetYTitle("N_{fit}");
01044   hNCommonHitVsNHit->SetZTitle("N_{common}");
01045   utility->setStyle(hNCommonHitVsNHit);
01046   mhNCommonHitVsNHit[categoryid].insert( std::make_pair(idcollection, hNCommonHitVsNHit) );
01047 
01048   // Dca vs eta vs MC pt
01049   title = Form("Dca vs #eta vs MC p_{T} (N_{fit}#geq10), %s", CategoryAndGeantId.Data());
01050   if( isMc ) title = Form("Dca vs #eta vs MC p_{T}, %s", CategoryAndGeantId.Data());
01051   else if ( isEmbedding ) title = Form("Dca vs #eta vs MC p_{T} (N_{fit}#geq10 & N_{common}#geq10), %s", CategoryAndGeantId.Data());
01052 
01053   TH3* hDca = new TH3D(Form("hDca%s", nameSuffix.Data()), title, 10, 0, 5, 10, -1.0, 1.0, 100, 0, 3.0);
01054   hDca->SetXTitle("MC p_{T} (GeV/c)");
01055   hDca->SetYTitle("#eta");
01056   hDca->SetZTitle("Global dca (cm)");
01057   utility->setStyle(hDca);
01058   mhDca[categoryid].insert( std::make_pair(idcollection, hDca) );
01059 
01060   //--------------------------------------------------
01061   // Common cuts for the 
01062   //  - pt vs eta
01063   //  - pt vs y
01064   //  - pt vs phi
01065   //  - pt vs momentum
01066   //  - delta pt vs pt
01067   //  - momentum vs eta
01068   //  - dE/dx vs momentum
01069   //  - reco. momentum vs MC momentum
01070   //  - eta vs phi
01071   //  - eta vs vz
01072   //  - rapidity vs vz
01073   //
01074   //  - dE/dx vs momentum (with pid cuts for the real data)
01075   //--------------------------------------------------
01076   TString cut(" (N_{fit}#geq10 & |dcaGl|<3cm)"); // real data
01077   if( isMc ) cut = "";
01078   else if ( isEmbedding ) cut = "(N_{fit}#geq10 & N_{common}#geq10 & |dcaGl|<3cm)";
01079 
01080   const TString titleSuffix(Form("%s, %s", cut.Data(), CategoryAndGeantId.Data()));
01081 
01082   // pt vs eta
01083   TH2* hPtVsEta = new TH2D(Form("hPtVsEta%s", nameSuffix.Data()), Form("MC p_{T} vs #eta%s", titleSuffix.Data()),
01084       etaBin, etaMin, etaMax, ptBin, ptMin, ptMax);
01085   hPtVsEta->SetXTitle("#eta");
01086   hPtVsEta->SetYTitle("MC p_{T} (GeV/c)");
01087   utility->setStyle(hPtVsEta);
01088   mhPtVsEta[categoryid].insert( std::make_pair(idcollection, hPtVsEta) );
01089 
01090   // pt vs y
01091   TH2* hPtVsY = new TH2D(Form("hPtVsY%s", nameSuffix.Data()), Form("MC p_{T} vs y%s", titleSuffix.Data()),
01092       etaBin, etaMin, etaMax, ptBin, ptMin, ptMax);
01093   hPtVsY->SetXTitle("rapidity y");
01094   hPtVsY->SetYTitle("MC p_{T} (GeV/c)");
01095   utility->setStyle(hPtVsY);
01096   mhPtVsY[categoryid].insert( std::make_pair(idcollection, hPtVsY) );
01097 
01098   // pt vs phi
01099   TH2* hPtVsPhi = new TH2D(Form("hPtVsPhi%s", nameSuffix.Data()), Form("MC p_{T} vs #phi%s", titleSuffix.Data()),
01100       100, -TMath::Pi(), TMath::Pi(), ptBin, ptMin, ptMax);
01101   hPtVsPhi->SetXTitle("#phi (rad)");
01102   hPtVsPhi->SetYTitle("MC p_{T} (GeV/c)");
01103   utility->setStyle(hPtVsPhi);
01104   mhPtVsPhi[categoryid].insert( std::make_pair(idcollection, hPtVsPhi) );
01105 
01106   // pt vs momentum
01107   TH2* hPtVsMom = new TH2D(Form("hPtVsMom%s", nameSuffix.Data()), Form("MC p_{T} vs momentum%s", titleSuffix.Data()),
01108       ptBin, ptMin, ptMax, ptBin, ptMin, ptMax);
01109   hPtVsMom->SetXTitle("momentum (GeV/c)");
01110   hPtVsMom->SetYTitle("MC p_{T} (GeV/c)");
01111   utility->setStyle(hPtVsMom);
01112   mhPtVsMom[categoryid].insert( std::make_pair(idcollection, hPtVsMom) );
01113 
01114   // Delta pt vs pt
01115   TH2* hdPtVsPt = new TH2D(Form("hdPtVsPt%s", nameSuffix.Data()), Form("p_{T} - p_{T} (MC) vs p_{T}%s", titleSuffix.Data()),
01116       ptBin, ptMin, ptMax, 100, -5, 5);
01117   hdPtVsPt->SetXTitle("MC p_{T} (GeV/c)");
01118   hdPtVsPt->SetYTitle("reco. p_{T} - MC p_{T} (GeV/c)");
01119   utility->setStyle(hdPtVsPt);
01120   mhdPtVsPt[categoryid].insert( std::make_pair(idcollection, hdPtVsPt) );
01121 
01122   // Delta 1/pt vs pt
01123   TH2* hdInvPtVsPt = new TH2D(Form("hdInvPtVsPt%s", nameSuffix.Data()), Form("1/p_{T} (Gl) - 1/p_{T} (MC) vs p_{T}%s", titleSuffix.Data()),
01124       ptBin, ptMin, ptMax, 200, -50, 50);
01125   hdInvPtVsPt->SetXTitle("MC p_{T} (GeV/c)");
01126   hdInvPtVsPt->SetYTitle("Gl 1/p_{T} - MC 1/p_{T} (c/MeV)");
01127   utility->setStyle(hdInvPtVsPt);
01128   mhdInvPtVsPt[categoryid].insert( std::make_pair(idcollection, hdInvPtVsPt) );
01129 
01130   // momentum vs eta
01131   TH2* hMomVsEta = new TH2D(Form("hMomVsEta%s", nameSuffix.Data()), Form("Momentum vs #eta%s", titleSuffix.Data()),
01132       etaBin, etaMin, etaMax, ptBin, ptMin, ptMax);
01133   hMomVsEta->SetXTitle("#eta");
01134   hMomVsEta->SetYTitle("momentum (GeV/c)");
01135   utility->setStyle(hMomVsEta);
01136   mhMomVsEta[categoryid].insert( std::make_pair(idcollection, hMomVsEta) );
01137 
01138   // dE/dx vs MC momentum
01139   TH2* hdEdxVsMomMc = new TH2D(Form("hdEdxVsMomMc%s", nameSuffix.Data()), Form("dE/dx vs MC p%s", titleSuffix.Data()),
01140       ptBin, ptMin, ptMax, 100, 0, 10.0);
01141   hdEdxVsMomMc->SetXTitle("MC p (GeV/c)");
01142   hdEdxVsMomMc->SetYTitle("dE/dx (keV/cm)");
01143   utility->setStyle(hdEdxVsMomMc);
01144   mhdEdxVsMomMc[categoryid].insert( std::make_pair(idcollection, hdEdxVsMomMc) );
01145 
01146   // dE/dx vs reconstructed momentum
01147   TH2* hdEdxVsMomReco = new TH2D(Form("hdEdxVsMomReco%s", nameSuffix.Data()), Form("dE/dx vs Reconstructed p%s", titleSuffix.Data()),
01148       ptBin, ptMin, ptMax, 100, 0, 10.0);
01149   hdEdxVsMomReco->SetXTitle("Reconstructed p (GeV/c)");
01150   hdEdxVsMomReco->SetYTitle("dE/dx (keV/cm)");
01151   utility->setStyle(hdEdxVsMomReco);
01152   mhdEdxVsMomReco[categoryid].insert( std::make_pair(idcollection, hdEdxVsMomReco) );
01153 
01154   // Reconstructed momentum vs MC momentum
01155   TH2* hRecoPVsMcP = new TH2D(Form("hRecoPVsMcP%s", nameSuffix.Data()), Form("Reconstructed p vs MC p%s", titleSuffix.Data()),
01156       ptBin, ptMin, ptMax, ptBin, ptMin, ptMax);
01157   hRecoPVsMcP->SetXTitle("MC p (GeV/c)");
01158   hRecoPVsMcP->SetYTitle("Reconstructed p (GeV/c)");
01159   utility->setStyle(hRecoPVsMcP);
01160   mhRecoPVsMcP[categoryid].insert( std::make_pair(idcollection, hRecoPVsMcP) );
01161 
01162   // eta vs phi
01163   TH2* hEtaVsPhi = new TH2D(Form("hEtaVsPhi%s", nameSuffix.Data()), Form("#eta vs #phi%s", titleSuffix.Data()),
01164       100, -TMath::Pi(), TMath::Pi(), etaBin, etaMin, etaMax);
01165   hEtaVsPhi->SetXTitle("#phi (rad)");
01166   hEtaVsPhi->SetYTitle("#eta");
01167   utility->setStyle(hEtaVsPhi);
01168   mhRecoPVsMcP[categoryid].insert( std::make_pair(idcollection, hRecoPVsMcP) );
01169   mhEtaVsPhi[categoryid].insert( std::make_pair(idcollection, hEtaVsPhi) );
01170 
01171   // eta vs vz
01172   TH2* hEtaVsVz = new TH2D(Form("hEtaVsVz%s", nameSuffix.Data()), Form("#eta vs v_{z}%s", titleSuffix.Data()),
01173       200, -50, 50, 200, etaMin, etaMax);
01174   hEtaVsVz->SetXTitle("v_{z} (cm)");
01175   hEtaVsVz->SetYTitle("#eta");
01176   utility->setStyle(hEtaVsVz);
01177   mhEtaVsVz[categoryid].insert( std::make_pair(idcollection, hEtaVsVz) );
01178  
01179   // rapidity vs vz
01180   TH2* hYVsVz = new TH2D(Form("hYVsVz%s", nameSuffix.Data()), Form("rapidity y vs v_{z}%s", titleSuffix.Data()),
01181       200, -50, 50, 200, etaMin, etaMax);
01182   hYVsVz->SetXTitle("v_{z} (cm)");
01183   hYVsVz->SetYTitle("rapidity y");
01184   utility->setStyle(hYVsVz);
01185   mhYVsVz[categoryid].insert( std::make_pair(idcollection, hYVsVz) );
01186 
01187   // dE/dx vs MC momentum (with pid cut)
01188   TH2* hdEdxVsMomMcPidCut = new TH2D(Form("hdEdxVsMomMcPidCut%s", nameSuffix.Data()), Form("dE/dx vs MC p (with 2#sigma pid cut)%s", titleSuffix.Data()),
01189       ptBin, ptMin, ptMax, 100, 0, 10.0);
01190   hdEdxVsMomMcPidCut->SetXTitle("MC p (GeV/c)");
01191   hdEdxVsMomMcPidCut->SetYTitle("dE/dx (keV/cm)");
01192   utility->setStyle(hdEdxVsMomMcPidCut);
01193   mhdEdxVsMomMcPidCut[categoryid].insert( std::make_pair(idcollection, hdEdxVsMomMcPidCut) );
01194 
01195   // dE/dx vs Reconstructed momentum (with pid cut)
01196   TH2* hdEdxVsMomRecoPidCut = new TH2D(Form("hdEdxVsMomRecoPidCut%s", nameSuffix.Data()), 
01197       Form("dE/dx vs Reconstructed p (with 2#sigma pid cut)%s", titleSuffix.Data()),
01198       ptBin, ptMin, ptMax, 100, 0, 10.0);
01199   hdEdxVsMomRecoPidCut->SetXTitle("Reconstructed p (GeV/c)");
01200   hdEdxVsMomRecoPidCut->SetYTitle("dE/dx (keV/cm)");
01201   utility->setStyle(hdEdxVsMomRecoPidCut);
01202   mhdEdxVsMomRecoPidCut[categoryid].insert( std::make_pair(idcollection, hdEdxVsMomRecoPidCut) );
01203 }
01204 
01205 //__________________________________________________________________________________________
01206 Int_t StEmbeddingQA::getNtrack(const Int_t categoryid, const StMiniMcEvent& mcevent) const
01207 {
01208   switch ( categoryid ) {
01209     case 0: return mcevent.nMcTrack() ;
01210     case 1: return mcevent.nMatchedPair() ;
01211     case 2: return mcevent.nGhostPair() ;
01212     case 3: return mcevent.nContamPair() ;
01213     // NOTE:
01214     //   Currently, there is no function to get number of matched global paris in StMiniMcEvent
01215     //   Skip it
01216     case 4: return 0 ;
01217     default:
01218       Warning("getNtrack", "Unkown category id, id=%3d", categoryid);
01219       return 0 ;
01220   }
01221 
01222   return 0 ;
01223 }
01224 
01225 //__________________________________________________________________________________________
01226 TString StEmbeddingQA::getIdCollection(const Int_t geantid, const Int_t parentid, const Int_t parentparentid) const
01227 {
01228   return Form("%d_%d_%d", geantid, parentid, parentparentid) ;
01229 }
01230