#include "StMuDstQAMaker.h" #include "StBTofHeader.h" #include "TBranch.h" ClassImp(StMuDstQAMaker) // // Set maximum file size to 1.9 GB (Root has a 2GB limit) // #define MAXFILESIZE 1900000000 StMuDstQAMaker::StMuDstQAMaker(StMuDstMaker *muMaker, const Char_t *oFileName) { std::cout << "StMuDstQAMaker::StMuDstQAMaker - Creating an instance..."; mMuDstMaker = muMaker; mAcceptTriggerFail = 0; mRefMultFail = 0; mAntiPileupFail = 0; mRankingFail = 0; mVtxNullFail = 0; mNPrimVtxFail = 0; mAcceptPrimVtxFail = 0; mEventIsGood = false; mNEventsIn = 0; mNEventsPassed = 0; mIsGoodTrack = false; mCurrentTrigger = 0; mTrackLoop = true; mTriggerCut = true; mFileName = oFileName; // // Initialize event cut variables // mPrimVtxZ[0] = -70.; mPrimVtxZ[1] = 70.; mPrimVtxR[0] = -1.; mPrimVtxR[1] = 10.; mPrimVtxVpdVzDiff[0] = -10.; mPrimVtxVpdVzDiff[1]= 10.; mPrimVtxXShift = 0.; mPrimVtxYShift = 0.; // // Initialize single-particle cut variables // mTrackP[0] = 0.1; mTrackP[1] = 2.; mTrackDca[0] = 0.; mTrackDca[1] = 5.; mTrackDcaGlobal[0] = 0.; mTrackDcaGlobal[1] = 5.; mTrackNHits[0] = 15; mTrackNHits[1] = 50; mTrackNHitsFit[0] = 15; mTrackNHitsFit[1] = 50; mTrackEta[0] = -1.1; mTrackEta[1] = 1.1; mTrackFlag[0] = 0; mTrackFlag[1] = 1000; /* // // Initialize TPC and TOF PID // mPionPionNSigma[0] = -3.; mPionPionNSigma[1] = 3.; mPionKaonNSigma[0] = 1.; mPionKaonNSigma[1] = -1.; */ mTTTPThreshold = 0.7; std::cout << "\t[DONE]" << std::endl; } //_________________ StMuDstQAMaker::~StMuDstQAMaker() { /* nothing to do */ } //_________________ Int_t StMuDstQAMaker::Init() { std::cout << "StMuDstQAMaker::Init - Initializing the maker" << std::endl; // // Create histograms and output file // mOutFile = new TFile(mFileName, "RECREATE"); // // General event distributions (before event cuts) // hR_BeforeCut = new TH1F("hR_BeforeCut", "R=#sqrt{v_{x}^{2} + v_{y}^{2}};R (cm);#", 160, 0., 0.8); hVx_BeforeCut = new TH1F("hVx_BeforeCut", ";v_{x} (cm);#", 320, -0.8, 0.8); hVy_BeforeCut = new TH1F("hVy_BeforeCut", ";v_{y} (cm);#", 320, -0.8, 0.8); hVz_BeforeCut = new TH1F("hVz_BeforeCut", ";v_{z} (cm);#", 400, -200., 200.); hVxVsVy_BeforeCut = new TH2F("hVxVsVy_BeforeCut", ";v_{x} (cm); v_{y} (cm)", 320, -0.8, 0.8, 320, -0.8, 0.8); hVxVsVz_BeforeCut = new TH2F("hVxVsVz_BeforeCut", ";v_{x} (cm); v_{z} (cm)", 600, -0.8, 0.8, 400, -200., 200.); hVyVsVz_BeforeCut = new TH2F("hVyVsVz_BeforeCut", ";v_{y} (cm); v_{z} (cm)", 320, -0.8, 0.8, 400, -200., 200.); hVpdVz2D_BeforeCut = new TH2F("hVpdVz2D_BeforeCut", ";vpd_{z} (cm); v_{z} (cm)", 400, -200., 200., 400, -200., 200.); hNPrimTr_BeforeCut = new TH1F("hNPrimTr_BeforeCut", ";N_{primary tracks};#", 150, 0, 150); hNGlobTr_BeforeCut = new TH1F("hNGlobTr_BeforeCut", ";N_{global tracks};#", 200, 0, 2000); hTofRefMult_BeforeCut = new TH1F("hTofRefMult_BeforeCut", ";TofRefMult;#", 100, 0, 1000); hTofRefMultVsRefMult_BeforeCut = new TH2F("hTofRefMultVsRefMult_BeforeCut", ";TofRefMult;RefMult", 500, 0., 1000, 500, 0., 1000); hVpd_BeforeCut = new TH1F("hVpd_BeforeCut", ";v_{z}^{vpd} (cm);#", 400, -200., 200.); hVpdVz_BeforeCut = new TH1F("hVpdVz_BeforeCut", ";v_{z}^{vpd}-v_{z}", 200, -20., 20.); hNPrimVtx_BeforeCut = new TH1F("hNPrimVtx_BeforeCut", ";N_{primary vertex};#", 15, 0, 15); hTriggerId_BeforePileupCut = new TH1I("hTriggerId_BeforePileupCut", ";trigger ID;#", 32, 0, 32); hTriggerId_BeforeAllCut = new TH1I("hTriggerId_BeforeAllCut", ";trigger ID;#", 32, 0, 32); hTriggerId_BeforeVertexCut = new TH1I("hTriggerId_BeforeVertexCut", ";trigger ID;#", 32, 0, 32); // // General event distributions (after event cuts) // hR_AfterCut = new TH1F("hR_AfterCut", "R=#sqrt{v_{x}^{2} + v_{y}^{2}};R (cm);#", 160, 0., 0.8); hVx_AfterCut = new TH1F("hVx_AfterCut", ";v_{x} (cm);#", 320, -0.8, 0.8); hVy_AfterCut = new TH1F("hVy_AfterCut", ";v_{y} (cm);#", 320, -0.8, 0.8); hVz_AfterCut = new TH1F("hVz_AfterCut", ";v_{z} (cm);#", 400, -200., 200.); hVxVsVy_AfterCut = new TH2F("hVxVsVy_AfterCut", ";v_{x} (cm); v_{y} (cm)", 320, -0.8, 0.8, 320, -0.8, 0.8); hVxVsVz_AfterCut = new TH2F("hVxVsVz_AfterCut", ";v_{x} (cm); v_{z} (cm)", 320, -0.8, 0.8, 400, -200., 200.); hVyVsVz_AfterCut = new TH2F("hVyVsVz_AfterCut", ";v_{y} (cm); v_{z} (cm)", 320, -0.8, 0.8, 400, -200., 200.); hVpdVz2D_AfterCut = new TH2F("hVpdVz2D_AfterCut", ";vpd_{z} (cm); v_{z} (cm)", 400, -200., 200., 400, -200., 200.); hNPrimTr_AfterCut = new TH1F("hNPrimTr_AfterCut", ";N_{primary tracks};#", 150, 0, 150); hNGlobTr_AfterCut = new TH1F("hNGlobTr_AfterCut", ";N_{global tracks};#", 200, 0, 2000); hTofRefMult_AfterCut = new TH1F("hTofRefMult_AfterCut", ";TofRefMult;#", 100, 0, 1000); hTofRefMultVsRefMult_AfterCut = new TH2F("hTofRefMultVsRefMult_AfterCut", ";TofRefMult;RefMult", 500, 0., 1000, 500, 0., 1000); hVpd_AfterCut = new TH1F("hVpd_AfterCut", ";v_{z}^{vpd} (cm);#", 400, -200., 200.); hVpdVz_AfterCut = new TH1F("hVpdVz_AfterCut", ";v_{z}^{vpd}-v_{z}", 200, -20., 20.); hNPrimVtx_AfterCut = new TH1F("hNPrimVtx_AfterCut", ";N_{primary vertex};#", 15, 0, 15); hTriggerId_AfterCut = new TH1I("hTriggerId_AfterCut", ";trigger ID;#", 32, 0, 32); // // General track distributions // if (mTrackLoop) { hNHitsVsMom = new TH2F("hNHitsVsMom", ";p (GeV/c);NHits", 400, -2.0, 2.0, 50, 0., 50.); pNHitsVsMom = new TProfile("pNHitsVsMom", ";p (GeV/c);", 400, -2.0, 2.0); hP = new TH1F("hP", ";p (GeV/c);#", 200, 0., 5.); hPt = new TH1F("hPt", ";p_{t} (GeV/c);#", 200, 0., 5.); hPx = new TH1F("hPx", ";p_{x} (GeV/c);#", 200, 0., 5.); hPy = new TH1F("hPy", ";p_{y} (GeV/c);#", 200, 0., 5.); hPz = new TH1F("hPz", ";p_{z} (GeV/c);#", 200, 0., 5.); hEta = new TH1F("hEta", ";#eta;#", 200, -2., 2.); hPAccept = new TH1F("hPAccept", "accepted tracks;p (GeV/c);#", 200, 0., 3.); hPtAccept = new TH1F("hPtAccept", "accepted tracks;p_{t} (GeV/c);#", 200, 0., 3.); hPxAccept = new TH1F("hPxAccept", "accepted tracks;p_{x} (GeV/c);#", 200, 0., 3.); hPyAccept = new TH1F("hPyAccept", "accepted tracks;p_{y} (GeV/c);#", 200, 0., 3.); hPzAccept = new TH1F("hPzAccept", "accepted tracks;p_{z} (GeV/c);#", 200, 0., 3.); hEtaAccept = new TH1F("hEtaAccept", "accepted tracks;#eta;#", 200, -2., 2.); hMassSqrVsPt = new TH2F("hMassSqrVsPt", ";p_{t}Q (GeV/c);m^{2} (GeV/c^{2})", 400, -2., 2., 1000, -0.1, 1.5); hDedxVsPt = new TH2F("hDedxVsPt", ";p_{t}Q (GeV/c);dE/dx (keV/cm)", 400, -2., 2., 400, 0., 20.); hInvBetaExpVsPt = new TH2F("hInvBetaExpVsPt", ";p_{t}Q (GeV/c);1/beta_{exp}", 400, -2., 2., 200, 0., 2.); hInvBetaThVsPt = new TH2F("hInvBetaThVsPt", ";p_{t}Q (GeV/c);1/beta_{exp}", 400, -2., 2., 200, 0., 2.); hTOF = new TH1F("hTOF", "time of flight;t ns;#", 100, 0., 50.); hDiffBeta = new TH1F("hDiffBeta", "difference in #beta between primary and global tracks;#beta_{primary} - #beta_{global};counts", 5000., -10., 10.); hDiffTOF = new TH1F("hDiffTOF", "difference in TOF between primary and global tracks;TOF_{primary} - TOF_{global};counts", 5000., -100., 100.); // Pion PID hNSigmaPionPion = new TH1F("hNSigmaPionPion", "#pi;n#sigma(#pi^{#pm});#", 100, -5., 5.); hNSigmaPionKaon = new TH1F("hNSigmaPionKaon", "#pi;n#sigma(#K^{#pm});#", 100, -5., 5.); hNSigmaPionProton = new TH1F("hNSigmaPionProton", "#pi;n#sigma(p);#", 100, -5., 5.); hNSigmaPionPionVsPt = new TH2F("hNSigmaPionPionVsPt", "#pi;p_{t}Q (GeV/c);n#sigma(#pi^{#pm})", 400, -2., 2., 1000, -5., 5.); hNSigmaPionKaonVsPt = new TH2F("hNSigmaPionKaonVsPt", "#pi;p_{t}Q (GeV/c);n#sigma(K^{#pm})", 400, -2., 2., 1000, -5., 5.); hNSigmaPionProtonVsPt = new TH2F("hNSigmaPionProtonVsPt", "#pi;p_{t}Q (GeV/c);n#sigma(p)", 400, -2., 2., 1000, -5., 5.); hDMassSqrVsPtPionTPC = new TH2F("hDMassSqrVsPtPionTPC", "TPC PID;p_{t} (GeV/c);dM(#pi^{#pm}) (GeV/c^{2})", 400, -2., 2., 500, -.5, .5); hDMassSqrVsPtPionTOF = new TH2F("hDMassSqrVsPtPionTOF", "TOF PID;p_{t} (GeV/c);dM(#pi^{#pm}) (GeV/c^{2})", 400, -2., 2., 500, -.5, .5); hInvBetaExpThPionTNT = new TH2F("hInvBetaExpThPionTNT", "TPC and TOF #pi^{#pm} PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}", 400, -2., 2., 500, -0.5, 0.5); hInvBetaExpThPionTPC = new TH2F("hInvBetaExpThPionTPC", "only TPC #pi^{#pm} PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}", 400, -2., 2., 500, -0.5, 0.5); hInvBetaExpThPionTOF = new TH2F("hInvBetaExpThPionTOF", "only TOF #pi^{#pm} PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}", 400, -2., 2., 500, -0.5, 0.5); hInvBetaExpThPionTTT = new TH2F("hInvBetaExpThPionTTT", Form("if p #leq %.2f then TPC, else TOF #pi^{#pm} PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}", mTTTPThreshold), 400, -2., 2., 500, -0.5, 0.5); hInvBetaExpThPionINT = new TH2F("hInvBetaExpThPionINT", "#pi^{#pm};p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}", 400, -2., 2., 500, -5., 5.); // Kaon PID hNSigmaKaonPion = new TH1F("hNSigmaKaonPion", "K;n#sigma(#pi^{#pm}};#", 100, -5., 5.); hNSigmaKaonKaon = new TH1F("hNSigmaKaonKaon", "K;n#sigma(#K^{#pm}};#", 100, -5., 5.); hNSigmaKaonProton = new TH1F("hNSigmaKaonProton", "K;n#sigma(p);#", 100, -5., 5.); hNSigmaKaonPionVsPt = new TH2F("hNSigmaKaonPionVsPt", "K;p_{t}Q (GeV/c);n#sigma(#pi^{#pm})", 400, -2., 2., 1000, -5., 5.); hNSigmaKaonKaonVsPt = new TH2F("hNSigmaKaonKaonVsPt", "K;p_{t}Q (GeV/c);n#sigma(K^{#pm})", 400, -2., 2., 1000, -5., 5.); hNSigmaKaonProtonVsPt = new TH2F("hNSigmaKaonProtonVsPt", "K;p_{t}Q (GeV/c);n#sigma(p)", 400, -2., 2., 1000, -5., 5.); hDMassSqrVsPtKaonTPC = new TH2F("hDMassSqrVsPtKaonTPC", "TPC PID;p_{t} (GeV/c);dM(K^{#pm}) (GeV/c^{2})", 400, -2., 2., 500, -.5, .5); hDMassSqrVsPtKaonTOF = new TH2F("hDMassSqrVsPtKaonTOF", "TOF PID;p_{t} (GeV/c);dM(K^{#pm}) (GeV/c^{2})", 400, -2., 2., 500, -.5, .5); hInvBetaExpThKaonTNT = new TH2F("hInvBetaExpThKaonTNT", "TPC and TOF K^{#pm} PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}", 400, -2., 2., 500, -0.5, 0.5); hInvBetaExpThKaonTPC = new TH2F("hInvBetaExpThKaonTPC", "only TPC K^{#pm} PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}", 400, -2., 2., 500, -0.5, 0.5); hInvBetaExpThKaonTOF = new TH2F("hInvBetaExpThKaonTOF", "only TOF K^{#pm} PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}", 400, -2., 2., 500, -0.5, 0.5); hInvBetaExpThKaonTTT = new TH2F("hInvBetaExpThKaonTTT", Form("if p #leq %.2f then TPC, else TOF K^{#pm} PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}", mTTTPThreshold), 400, -2., 2., 500, -0.5, 0.5); hInvBetaExpThKaonINT = new TH2F("hInvBetaExpThKaonINT", "K^{#pm};p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}", 400, -2., 2., 500, -5., 5.); // Proton PID hNSigmaProtonPion = new TH1F("hNSigmaProtonPion", "p;n#sigma(#pi^{#pm}};#", 100, -5., 5.); hNSigmaProtonKaon = new TH1F("hNSigmaProtonKaon", "p;n#sigma(#K^{#pm}};#", 100, -5., 5.); hNSigmaProtonProton = new TH1F("hNSigmaProtonProton", "p;n#sigma(p);#", 100, -5., 5.); hNSigmaProtonPionVsPt = new TH2F("hNSigmaProtonPionVsPt", "p;p_{t}Q (GeV/c);n#sigma(#pi^{#pm})", 400, -2., 2., 1000, -5., 5.); hNSigmaProtonKaonVsPt = new TH2F("hNSigmaProtonKaonVsPt", "p;p_{t}Q (GeV/c);n#sigma(K^{#pm})", 400, -2., 2., 1000, -5., 5.); hNSigmaProtonProtonVsPt = new TH2F("hNSigmaProtonProtonVsPt", "p;p_{t}Q (GeV/c);n#sigma(p)", 400, -2., 2., 1000, -5., 5.); hDMassSqrVsPtProtonTPC = new TH2F("hDMassSqrVsPtProtonTPC", "TPC PID;p_{t} (GeV/c);dM(p) (GeV/c^{2})", 400, -2., 2., 500, -.5, .5); hDMassSqrVsPtProtonTOF = new TH2F("hDMassSqrVsPtProtonTOF", "TOF PID;p_{t} (GeV/c);dM(p) (GeV/c^{2})", 400, -2., 2., 500, -.5, .5); hInvBetaExpThProtonTNT = new TH2F("hInvBetaExpThProtonTNT", "TPC and TOF p PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}", 400, -2., 2., 500, -0.5, 0.5); hInvBetaExpThProtonTPC = new TH2F("hInvBetaExpThProtonTPC", "only TPC p PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}", 400, -2., 2., 500, -0.5, 0.5); hInvBetaExpThProtonTOF = new TH2F("hInvBetaExpThProtonTOF", "only TOF p PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}", 400, -2., 2., 500, -0.5, 0.5); hInvBetaExpThProtonTTT = new TH2F("hInvBetaExpThProtonTTT", Form("if p #leq %.2f then TPC, else TOF p PID;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}", mTTTPThreshold), 400, -2., 2., 500, -0.5, 0.5); hInvBetaExpThProtonINT = new TH2F("hInvBetaExpThProtonINT", "p;p_{t}Q (GeV/c);1/beta_{exp} - 1/beta_{th}", 400, -2., 2., 500, -5., 5.); } // // Any charge // hRefMult_BeforeCut = new TH1F("hRefMult_BeforeCut", ";RefMult;#", 100, 0., 1000.); hRefMultAccept_BeforeCut = new TH1F("hRefMultAccept_BeforeCut", ";RefMultAccept;#", 100, 0., 1000.); hRefMult_AfterCut = new TH1F("hRefMult_AfterCut", ";RefMult;#", 100, 0., 1000.); hRefMultAccept_AfterCut = new TH1F("hRefMultAccept_AfterCut", ";RefMultAccept;#", 100, 0., 1000.); std::cout << "StMuDstQAMaker::Init - Initialization has been finished" << std::endl; return StMaker::Init(); } //________________ void StMuDstQAMaker::Clear(Option_t *option) { StMaker::Clear(); } //________________ Int_t StMuDstQAMaker::Make() { mNEventsIn++; mMuDst = NULL; mMuEvent = NULL; //MuDstMaker initialization if(!mMuDstMaker) { mMuDstMaker = (StMuDstMaker*)GetMaker("MuDst"); if(!mMuDstMaker) { LOG_ERROR << "StMuDstQAMaker::Make [ERROR] - Cannot find StMuDstMaker" << std::endl; return kStOk; } } //Obtaining MuDst mMuDst = mMuDstMaker->muDst(); if (!mMuDst) { mMuDst = (StMuDst*)GetInputDS("MuDst"); if (!mMuDst) { gMessMgr->Warning() << "StFemtoDstMaker::Make [WARNING] - No MuDst has been found" << endm; return kStOk; } } //Obtaining MuEvent mMuEvent = (StMuEvent*)mMuDst->event(); //Multiplicity cannot be negative unsigned short refMult = mMuEvent->refMult(); unsigned short refMultPos = 0; unsigned short refMultNeg = 0; if(refMult < 0) { mRefMultFail++; std::cout << "[kStOk] RefMult less than zero!\n"; return kStOk; } int mNVertices = mMuDst->numberOfPrimaryVertices(); int mNPrimTracks = mMuDst->numberOfPrimaryTracks(); int mNGlobTracks = mMuDst->numberOfGlobalTracks(); if (mNVertices == 0) { return kStOk; } //Some initializations of local variables mPrimVertex = NULL; StThreeVectorF mVertPosition; Float_t mVpdVz = 0.; Int_t mPrimVertIndex = -999; Float_t mRanking = -999.; //Clean index vectors CleanVariables(); //Vertex loop unsigned short refMultAccept = 0; unsigned short tofRefMult = 0; unsigned short refMultAcceptPos = 0; unsigned short refMultAcceptNeg = 0; // Not first primary vertex does not contain tracks with fast detectors (TOF) int iVert = 0; mEventIsGood = false; // // Fill histogram hTriggerId before all cuts // for(unsigned int iTrg=0; iTrgtriggerIdCollection().nominal().isTrigger(mTriggerIdCollection.at(iTrg))) { hTriggerId_BeforeAllCut->Fill(iTrg); } } mPrimVertex = mMuDst->primaryVertex(iVert); //Positive ranking only if(mPrimVertex->ranking() <= 0) { mRankingFail++; return kStOk; } //Not (0,0,0) position of the primary vertex if(mPrimVertex->position().x() == 0 && mPrimVertex->position().y() == 0 && mPrimVertex->position().z() == 0) { mVtxNullFail++; std::cout << "[kStOk] Vertex equals to zero!\n"; return kStOk; } //Reasonable amount of tracks if(mNPrimTracks < 0 || mNPrimTracks > 10000) { mNPrimVtxFail++; std::cout << "[kStOk] N primary tracks cut failed!\n"; return kStOk; } mPrimVertIndex = iVert; mRanking = mPrimVertex->ranking(); mVertPosition = mPrimVertex->position(); StBTofHeader *btofH = mMuDst->btofHeader(); if (!btofH) { std::cout << "StBTofHeader - NULL pointer!\n"; return kStOk; } mVpdVz = mMuDst->btofHeader()->vpdVz(); hRefMult_BeforeCut->Fill(refMult); hNPrimVtx_BeforeCut->Fill(mNVertices); hNPrimTr_BeforeCut->Fill(mNPrimTracks); hNGlobTr_BeforeCut->Fill(mNGlobTracks); // // Fill histogram hTriggerId before pile-up cuts // for(unsigned int iTrg=0; iTrgtriggerIdCollection().nominal().isTrigger(mTriggerIdCollection.at(iTrg))) { hTriggerId_BeforePileupCut->Fill(iTrg); } } // // Calculate the amount of the pile-up vertices // int mNVtxInBunchCross = 0; for(unsigned int i = 0; i < mMuDst->numberOfPrimaryVertices(); i++) { StMuPrimaryVertex *primVtx = mMuDst->primaryVertex(i); if(!primVtx || primVtx->ranking() <= 0.) continue; if(TMath::Abs(primVtx->position().x()) < 1e-5 && TMath::Abs(primVtx->position().y()) < 1e-5 && TMath::Abs(primVtx->position().z()) < 1e-5) continue; if(TMath::Abs(primVtx->position().z() - btofH->vpdVz()) < 5.0) mNVtxInBunchCross++; } // There should be only one primary vertex that caused the trigger if(mNVtxInBunchCross > 1) { mAntiPileupFail++; return kStOk; } // // Fill histogram hTriggerId before cuts // for(unsigned int iTrg=0; iTrg < mTriggerIdCollection.size(); iTrg++) { if(mMuEvent->triggerIdCollection().nominal().isTrigger(mTriggerIdCollection.at(iTrg))) { hTriggerId_BeforeVertexCut->Fill(iTrg); } } // // Trigger cut // if (mTriggerCut) { if(!AcceptTrigger(mMuEvent)) //If trigger is not found { mAcceptTriggerFail++; return kStOk; } } if( !AcceptPrimVtx(mVertPosition, mVpdVz) ) { mAcceptPrimVtxFail++; return kStOk; } mEventIsGood = true; // // Fill histogram hTriggerId after cuts // for(unsigned int iTrg=0; iTrgtriggerIdCollection().nominal().isTrigger(mTriggerIdCollection.at(iTrg))) { hTriggerId_AfterCut->Fill(iTrg); } } // //Loop over primary tracks // if (mTrackLoop) { for(int iTrk = 0; iTrk < mNPrimTracks; iTrk++) { mPrimTrack = mMuDst->primaryTracks(iTrk); int idxGlob = mPrimTrack->index2Global(); if (idxGlob < 0) continue; mGlobTrack = mMuDst->globalTracks(idxGlob); short charge = mPrimTrack->charge() > 0 ? 1 : -1; float eta = mPrimTrack->eta(); float pt = mPrimTrack->pt(); float p = mPrimTrack->p().mag(); float p2 = mPrimTrack->p().mag2(); float px = mPrimTrack->p().x(); float py = mPrimTrack->p().y(); float pz = mPrimTrack->p().z(); charge > 0 ? refMultPos++ : refMultNeg++; hP->Fill(p); hPt->Fill(pt); hPx->Fill(px); hPy->Fill(py); hPz->Fill(pz); hEta->Fill(eta); if( !AcceptTrack(mPrimTrack, iVert) ) continue; hPAccept->Fill(p); hPtAccept->Fill(pt); hPxAccept->Fill(px); hPyAccept->Fill(py); hPzAccept->Fill(pz); hEtaAccept->Fill(eta); refMultAccept++; charge > 0 ? refMultAcceptPos++ : refMultAcceptNeg++; float nSigmaPion = mPrimTrack->nSigmaPion(); float nSigmaKaon = mPrimTrack->nSigmaKaon(); float nSigmaProton = mPrimTrack->nSigmaProton(); float dedx = mPrimTrack->dEdx(); bool tofTrack = IsTofTrack(mGlobTrack, mPrimTrack); float betaExp; float massSqr; bool pionTPC = ((nSigmaPion >= mPionPionNSigma[0] || nSigmaPion <= mPionPionNSigma[1]) && (nSigmaKaon < mPionKaonNSigma[0] || nSigmaKaon > mPionKaonNSigma[1]) && (nSigmaProton < mPionProtonNSigma[0] || nSigmaProton > mPionProtonNSigma[1])); bool kaonTPC = ((nSigmaPion < mKaonPionNSigma[0] || nSigmaPion > mKaonPionNSigma[1]) && (nSigmaKaon >= mKaonKaonNSigma[0] || nSigmaKaon <= mKaonKaonNSigma[1]) && (nSigmaProton < mKaonProtonNSigma[0] || nSigmaProton > mKaonProtonNSigma[1])); bool protonTPC = ((nSigmaPion < mProtonPionNSigma[0] || nSigmaPion > mProtonPionNSigma[1]) && (nSigmaKaon < mProtonKaonNSigma[0] || nSigmaKaon > mProtonKaonNSigma[1]) && (nSigmaProton >= mProtonProtonNSigma[0] || nSigmaProton <= mProtonProtonNSigma[1])); bool pionTOF = false; bool kaonTOF = false; bool protonTOF = false; // // Pion TPC // if (pionTPC) { hNSigmaPionPion->Fill(nSigmaPion); hNSigmaPionKaon->Fill(nSigmaKaon); hNSigmaPionProton->Fill(nSigmaProton); hNSigmaPionPionVsPt->Fill(pt*charge, nSigmaPion); hNSigmaPionKaonVsPt->Fill(pt*charge, nSigmaKaon); hNSigmaPionProtonVsPt->Fill(pt*charge, nSigmaProton); } // // Kaon TPC // if (kaonTPC) { hNSigmaKaonPion->Fill(nSigmaPion); hNSigmaKaonKaon->Fill(nSigmaKaon); hNSigmaKaonProton->Fill(nSigmaProton); hNSigmaKaonPionVsPt->Fill(pt*charge, nSigmaPion); hNSigmaKaonKaonVsPt->Fill(pt*charge, nSigmaKaon); hNSigmaKaonProtonVsPt->Fill(pt*charge, nSigmaProton); } // // Proton TPC // if (protonTPC) { hNSigmaProtonPion->Fill(nSigmaPion); hNSigmaProtonKaon->Fill(nSigmaKaon); hNSigmaProtonProton->Fill(nSigmaProton); hNSigmaProtonPionVsPt->Fill(pt*charge, nSigmaPion); hNSigmaProtonKaonVsPt->Fill(pt*charge, nSigmaKaon); hNSigmaProtonProtonVsPt->Fill(pt*charge, nSigmaProton); } hDedxVsPt->Fill(pt*charge, dedx*1e6); if (tofTrack) { hTOF->Fill(mGlobTrack->btofPidTraits().timeOfFlight()); float betaThPion = sqrt(p2/(PION_MASS*PION_MASS + p2)); float betaThKaon = sqrt(p2/(KAON_MASS*KAON_MASS + p2)); float betaThProton = sqrt(p2/(PROTON_MASS*PROTON_MASS + p2)); tofRefMult++; betaExp = mGlobTrack->btofPidTraits().beta(); massSqr = p2*(1./(betaExp*betaExp) - 1.); hMassSqrVsPt->Fill(pt*charge, massSqr); pionTOF = massSqr >= mPionMass[0] && massSqr <= mPionMass[1]; kaonTOF = massSqr >= mKaonMass[0] && massSqr <= mKaonMass[1]; protonTOF = massSqr >= mProtonMass[0] && massSqr <= mProtonMass[1]; float dMpion = massSqr - PION_MASS*PION_MASS; float dMkaon = massSqr - KAON_MASS*KAON_MASS; float dMproton = massSqr - PROTON_MASS*PROTON_MASS; hInvBetaExpVsPt->Fill(pt*charge, 1./betaExp); // // TPC PID // if (pionTPC) { hInvBetaExpThPionTPC->Fill(pt*charge, 1./betaExp - 1./betaThPion); hDMassSqrVsPtPionTPC->Fill(pt*charge, dMpion); } else if (kaonTPC) { hInvBetaExpThKaonTPC->Fill(pt*charge, 1./betaExp - 1./betaThKaon); hDMassSqrVsPtKaonTPC->Fill(pt, dMkaon); } else if (protonTPC) { hInvBetaExpThProtonTPC->Fill(pt*charge, 1./betaExp - 1./betaThProton); hDMassSqrVsPtProtonTPC->Fill(pt, dMproton); } // // TOF PID // if (pionTOF) { hInvBetaExpThPionTOF->Fill(pt*charge, 1./betaExp - 1./betaThPion); hDMassSqrVsPtPionTOF->Fill(pt, dMpion); hInvBetaThVsPt->Fill(pt*charge, 1./betaThPion); hInvBetaExpThPionINT->Fill(pt*charge, 1./betaExp - 1./betaThPion); hInvBetaExpThPionINT->Fill(pt*charge, 1./betaExp - 1./betaThKaon); hInvBetaExpThPionINT->Fill(pt*charge, 1./betaExp - 1./betaThProton); } else if (kaonTOF) { hInvBetaExpThKaonTOF->Fill(pt*charge, 1./betaExp - 1./betaThKaon); hDMassSqrVsPtKaonTOF->Fill(pt, dMkaon); hInvBetaThVsPt->Fill(pt*charge, 1./betaThKaon); hInvBetaExpThKaonINT->Fill(pt*charge, 1./betaExp - 1./betaThPion); hInvBetaExpThKaonINT->Fill(pt*charge, 1./betaExp - 1./betaThKaon); hInvBetaExpThKaonINT->Fill(pt*charge, 1./betaExp - 1./betaThProton); } else if (protonTOF) { hInvBetaExpThProtonTOF->Fill(pt*charge, 1./betaExp - 1./betaThProton); hDMassSqrVsPtProtonTOF->Fill(pt, dMproton); hInvBetaThVsPt->Fill(pt*charge, 1./betaThProton); hInvBetaExpThProtonINT->Fill(pt*charge, 1./betaExp - 1./betaThPion); hInvBetaExpThProtonINT->Fill(pt*charge, 1./betaExp - 1./betaThKaon); hInvBetaExpThProtonINT->Fill(pt*charge, 1./betaExp - 1./betaThProton); } // // TNT PID // if (pionTOF && pionTPC) { hInvBetaExpThPionTNT->Fill(pt*charge, 1./betaExp - 1./betaThPion); } else if (kaonTOF && kaonTPC) { hInvBetaExpThKaonTNT->Fill(pt*charge, 1./betaExp - 1./betaThKaon); } else if (protonTOF && protonTPC) { hInvBetaExpThProtonTNT->Fill(pt*charge, 1./betaExp - 1./betaThProton); } // // TTT PID // if (p >= mTTTPThreshold) { if (pionTOF) { hInvBetaExpThPionTTT->Fill(pt*charge, 1./betaExp - 1./betaThPion); } else if (kaonTOF) { hInvBetaExpThKaonTTT->Fill(pt*charge, 1./betaExp - 1./betaThKaon); } else if (protonTOF) { hInvBetaExpThProtonTTT->Fill(pt*charge, 1./betaExp - 1./betaThProton); } } else { if (pionTPC) { hInvBetaExpThPionTTT->Fill(pt*charge, 1./betaExp - 1./betaThPion); } else if (kaonTPC) { hInvBetaExpThKaonTTT->Fill(pt*charge, 1./betaExp - 1./betaThKaon); } else if (protonTPC) { hInvBetaExpThProtonTTT->Fill(pt*charge, 1./betaExp - 1./betaThProton); } } } } //for(Int_t iTrk=0; iTrkFill(refMultAccept); hTofRefMult_AfterCut->Fill(tofRefMult); hTofRefMultVsRefMult_AfterCut->Fill(tofRefMult, refMultAccept); } return kStOk; } //_________________ Bool_t StMuDstQAMaker::AcceptTrigger(StMuEvent *muEvent) { Bool_t mIsGoodTrigger = false; if(mTriggerIdCollection.empty()) { mIsGoodTrigger = true; } else { for(UInt_t iTrg=0; iTrgtriggerIdCollection().nominal().isTrigger(mTriggerIdCollection.at(iTrg))) { mIsGoodTrigger = true; mCurrentTrigger = mTriggerIdCollection.at(iTrg); break; } } //for(UInt_t iTrg=0; iTrgFill(vtxR); hVx_BeforeCut->Fill(mVtxX); hVy_BeforeCut->Fill(mVtxY); hVz_BeforeCut->Fill(mVtxZ); hVxVsVy_BeforeCut->Fill(mVtxX, mVtxY); hVxVsVz_BeforeCut->Fill(mVtxX, mVtxZ); hVyVsVz_BeforeCut->Fill(mVtxY, mVtxZ); hVpd_BeforeCut->Fill(vpdVz); hVpdVz_BeforeCut->Fill(vpdDiff); hVpdVz2D_BeforeCut->Fill(vpdVz, mVtxZ); bool isGoodEvent = ( vtxPos.z() >= mPrimVtxZ[0] && vtxPos.z() <= mPrimVtxZ[1] && vtxR >= mPrimVtxR[0] && vtxR <= mPrimVtxR[1] && vpdDiff >= mPrimVtxVpdVzDiff[0] && vpdDiff <= mPrimVtxVpdVzDiff[1] ); if (isGoodEvent) { hR_AfterCut->Fill(vtxR); hVx_AfterCut->Fill(mVtxX); hVy_AfterCut->Fill(mVtxY); hVz_AfterCut->Fill(mVtxZ); hVxVsVy_AfterCut->Fill(mVtxX, mVtxY); hVxVsVz_AfterCut->Fill(mVtxX, mVtxZ); hVyVsVz_AfterCut->Fill(mVtxY, mVtxZ); hVpd_AfterCut->Fill(vpdVz); hVpdVz_AfterCut->Fill(vpdDiff); hVpdVz2D_AfterCut->Fill(vpdVz, mVtxZ); } return isGoodEvent; } //_________________ Bool_t StMuDstQAMaker::AcceptTrack(StMuTrack *trk, UShort_t vtxInd) { float mom = trk->momentum().mag(); float eta = trk->eta(); unsigned short nHitsFit = trk->nHitsFit(); float ratioH = (float)trk->nHitsFit()/(float)trk->nHitsPoss(); short flag = trk->flag(); float dca = trk->dcaGlobal(vtxInd).perp(); short charge = trk->charge(); bool isMomOk = ( mom >= mTrackP[0] && mom <= mTrackP[1] ); bool isEtaOk = ( eta >= mTrackEta[0] && eta <= mTrackEta[1] ); bool isNHitsFitOk = ( nHitsFit >= mTrackNHitsFit[0] && nHitsFit <= mTrackNHitsFit[1] && ratioH >= 0.51 ); bool isFlagOk = ( flag >= mTrackFlag[0] && flag <= mTrackFlag[1] ); // bool isDcaOk = ( dca >= mTrackDcaGlobal[0] && // dca <= mTrackDcaGlobal[1] ); bool isGoodTrack = isMomOk && isEtaOk && isNHitsFitOk && isFlagOk; if (isGoodTrack) { hNHitsVsMom->Fill(mom*charge, nHitsFit); pNHitsVsMom->Fill(mom*charge, nHitsFit); } return isGoodTrack; } //_________________ Int_t StMuDstQAMaker::Finish() { std::cout << "*************************************" << "\n" << "StMuDstQAMaker has been finished" << "\n" << "\t nEventsPassed : " << mNEventsPassed << "\n" << "\t nEventsProcessed: " << mNEventsIn << "\n" << "Accept trigger fail = " << mAcceptTriggerFail << "\n" << "RefMult fail = " << mRefMultFail << "\n" << "Anti-pileup fail = " << mAntiPileupFail << "\n" << "Ranking fail = " << mRankingFail << "\n" << "VtxNull fail = " << mVtxNullFail << "\n" << "NPrimVtx fail = " << mNPrimVtxFail << "\n" << "AcceptPrimVtx fail = " << mAcceptPrimVtxFail << "\n" << "*************************************" << "\n"; mOutFile->cd(); mOutFile->Write(); mOutFile->Close(); return kStOk; } //_________________ void StMuDstQAMaker::CleanVariables() { mEventIsGood = false; } //_________________ void StMuDstQAMaker::SetTriggerId(unsigned int id) { mTriggerIdCollection.push_back(id); } //_________________ void StMuDstQAMaker::SetMuDstMaker(StMuDstMaker *maker) { mMuDstMaker = maker; } //_________________ void StMuDstQAMaker::SetVtxZCut(float lo, float hi) { mPrimVtxZ[0] = lo; mPrimVtxZ[1] = hi; } //_________________ void StMuDstQAMaker::SetVtxRCut(float lo, float hi) { mPrimVtxR[0] = lo; mPrimVtxR[1] = hi; } //_________________ void StMuDstQAMaker::SetVtxShift(float xShift, float yShift) { mPrimVtxXShift = xShift; mPrimVtxYShift = yShift; } //_________________ void StMuDstQAMaker::SetVtxVpdVzDiffCut(float lo, float hi) { mPrimVtxVpdVzDiff[0] = lo; mPrimVtxVpdVzDiff[1] = hi; } //_________________ void StMuDstQAMaker::SetP(float lo, float hi) { mTrackP[0] = lo; mTrackP[1] = hi; } //_________________ void StMuDstQAMaker::SetPt(float lo, float hi) { mTrackPt[0] = lo; mTrackPt[1] = hi; } //_________________ void StMuDstQAMaker::SetTrackNHits(int lo, int hi) { mTrackNHits[0] = lo; mTrackNHits[1] = hi; } //_________________ void StMuDstQAMaker::SetTrackNHitsFit(int lo, int hi) { mTrackNHitsFit[0] = lo; mTrackNHitsFit[1] = hi; } //_________________ void StMuDstQAMaker::SetTrackEta(float lo, float hi) { mTrackEta[0] = lo; mTrackEta[1] = hi; } //_________________ void StMuDstQAMaker::SetTrackFlag(short lo, short hi) { mTrackFlag[0] = lo; mTrackFlag[1] = hi; } //_________________ void StMuDstQAMaker::SetTrackDCA(float lo, float hi) { mTrackDcaGlobal[0] = lo; mTrackDcaGlobal[1] = hi; } //_________________ Bool_t StMuDstQAMaker::IsTofTrack(StMuTrack *trkGlob, StMuTrack *trkPrim) { Bool_t isTofHit = false; float betaGlobal = trkGlob->btofPidTraits().beta(); float tofGlobal = trkGlob->btofPidTraits().timeOfFlight(); float betaPrimary = trkPrim->btofPidTraits().beta(); float tofPrimary = trkPrim->btofPidTraits().timeOfFlight(); hDiffBeta->Fill(betaPrimary - betaGlobal); hDiffTOF->Fill(tofPrimary - tofGlobal); if(betaGlobal > 0 && tofGlobal > 0) isTofHit = true; return isTofHit; } //_________________ void StMuDstQAMaker::SetPionPionNSigma(float lo, float hi) { mPionPionNSigma[0] = lo; mPionPionNSigma[1] = hi; } //_________________ void StMuDstQAMaker::SetPionKaonNSigma(float lo, float hi) { mPionKaonNSigma[0] = lo; mPionKaonNSigma[1] = hi; } //_________________ void StMuDstQAMaker::SetPionProtonNSigma(float lo, float hi) { mPionProtonNSigma[0] = lo; mPionProtonNSigma[1] = hi; } //_________________ void StMuDstQAMaker::SetKaonPionNSigma(float lo, float hi) { mKaonPionNSigma[0] = lo; mKaonPionNSigma[1] = hi; } //_________________ void StMuDstQAMaker::SetKaonKaonNSigma(float lo, float hi) { mKaonKaonNSigma[0] = lo; mKaonKaonNSigma[1] = hi; } //_________________ void StMuDstQAMaker::SetKaonProtonNSigma(float lo, float hi) { mKaonProtonNSigma[0] = lo; mKaonProtonNSigma[1] = hi; } //_________________ void StMuDstQAMaker::SetProtonPionNSigma(float lo, float hi) { mProtonPionNSigma[0] = lo; mProtonPionNSigma[1] = hi; } //_________________ void StMuDstQAMaker::SetProtonKaonNSigma(float lo, float hi) { mProtonKaonNSigma[0] = lo; mProtonKaonNSigma[1] = hi; } //_________________ void StMuDstQAMaker::SetProtonProtonNSigma(float lo, float hi) { mProtonProtonNSigma[0] = lo; mProtonProtonNSigma[1] = hi; } //_________________ void StMuDstQAMaker::SetPionMassSqr(float lo, float hi) { mPionMass[0] = lo; mPionMass[1] = hi; } //_________________ void StMuDstQAMaker::SetKaonMassSqr(float lo, float hi) { mKaonMass[0] = lo; mKaonMass[1] = hi; } //_________________ void StMuDstQAMaker::SetProtonMassSqr(float lo, float hi) { mProtonMass[0] = lo; mProtonMass[1] = hi; } //_________________ void StMuDstQAMaker::SetTTTPThreshold(float pThresh) { mTTTPThreshold = pThresh; } void StMuDstQAMaker::SetTracksLoop(bool val) { mTrackLoop = val; } void StMuDstQAMaker::SetTriggerCut(bool val) { mTriggerCut = val; }