////////////////////////////////////////////////////////////////////////////// // // $Id: $ // // //_HADES_CLASS_DESCRIPTION ////////////////////////////////////////////////////////////////////////////// // // HDEAnaTask // // ////////////////////////////////////////////////////////////////////////////// #include "HDEAnaTask.h" #include "hades.h" #include "hruntimedb.h" #include "hcategory.h" #include "hevent.h" #include "hgeantrich.h" #include "hlinearcatiter.h" #include "hmatrixcatiter.h" #include "hparset.h" #include "hspectrometer.h" #include "richdef.h" #include "hparticledef.h" #include "hparticletracksorter.h" #include "hrichhit.h" #include "hrichhitsim.h" #include "hgeantkine.h" #include "hparticlecandsim.h" #include "hparticletool.h" #include "HDEAnaDrawHist.h" #include "HDEAnaUtils.h" #include "HDEAnaMyUtils.h" #include "TCanvas.h" #include "TH2D.h" #include "TPad.h" #include "TEllipse.h" #include "TRandom.h" #include "TLatex.h" #include #include using namespace std; using namespace HDEAna; ClassImp(HDEAnaTask) HDEAnaTask::HDEAnaTask(const string& outputFileName): fEventNum(0) { fOutputFileName = outputFileName; } HDEAnaTask::~HDEAnaTask() { } Bool_t HDEAnaTask::init() { fCatKine = gHades->getCurrentEvent()->getCategory(catGeantKine); if (NULL == fCatKine) { Error("init", "Initializatin of catGeantKine category failed, returning..."); return kFALSE; } fCatParticleCand = gHades->getCurrentEvent()->getCategory(catParticleCand); if (NULL == fCatParticleCand) { Error("init", "Initializatin of catParticleCand category failed, returning..."); return kFALSE; } fCatRichHit = gHades->getCurrentEvent()->getCategory(catRichHit); if (NULL == fCatRichHit) { Error("init", "Initialization of geant category for catRichHit failed, returning..."); return kFALSE; } HParticleTrackSorter sorter; sorter.setIgnoreInnerMDC(); // fSignalMBR =(7.89/3.5)*7.8e-8; fSignalMBR = 1.76e-7; fCuts.init(fCatParticleCand); initHist(); return kTRUE; } void HDEAnaTask::initHist() { fHM = new HDEAnaHistManager(); // Nof events fHM->Create1("fhNofEvents", "fhNofEvents;;Counter", 1, -.5, 0.5); // invariant mass histograms Int_t nofBinsMinv = 150; Double_t minMinv = 0.; Double_t maxMinv = 1500; // minv after Lepton cut fHM->Create1("fhMinvSignal", "fhMinvSignal;M_{ee} [MeV/c^{2}];Yield", nofBinsMinv, minMinv, maxMinv); fHM->Create1("fhMinvGamma", "fhMinvGamma;M_{ee} [MeV/c^{2}];Yield", nofBinsMinv, minMinv, maxMinv); fHM->Create1("fhMinvPi0", "fhMinvPi0;M_{ee} [MeV/c^{2}];Yield", nofBinsMinv, minMinv, maxMinv); fHM->Create1("fhMinvBg", "fhMinvBg;M_{ee} [MeV/c^{2}];Yield", nofBinsMinv, minMinv, maxMinv); // minv after momentum cut fHM->Create1("fhMinvMomCutSignal", "fhMinvMomCutSignal;M_{ee} [MeV/c^{2}];Yield", nofBinsMinv, minMinv, maxMinv); fHM->Create1("fhMinvMomCutGamma", "fhMinvMomCutGamma;M_{ee} [MeV/c^{2}];Yield", nofBinsMinv, minMinv, maxMinv); fHM->Create1("fhMinvMomCutPi0", "fhMinvMomCutPi0;M_{ee} [MeV/c^{2}];Yield", nofBinsMinv, minMinv, maxMinv); fHM->Create1("fhMinvMomCutBg", "fhMinvMomCutBg;M_{ee} [MeV/c^{2}];Yield", nofBinsMinv, minMinv, maxMinv); // minv after Opening angle cut fHM->Create1("fhMinvOACutSignal", "fhMinvOACutSignal;M_{ee} [MeV/c^{2}];Yield", nofBinsMinv, minMinv, maxMinv); fHM->Create1("fhMinvOACutGamma", "fhMinvoaCutGamma;M_{ee} [MeV/c^{2}];Yield", nofBinsMinv, minMinv, maxMinv); fHM->Create1("fhMinvOACutPi0", "fhMinvOACutPi0;M_{ee} [MeV/c^{2}];Yield", nofBinsMinv, minMinv, maxMinv); fHM->Create1("fhMinvOACutBg", "fhMinvOACutBg;M_{ee} [MeV/c^{2}];Yield", nofBinsMinv, minMinv, maxMinv); //Minus particle on x-axis: gamma e-, pi0 e-, pion-, other //Plus particle on y-axis: gamma e+, pi0 e+, pion+, other fHM->Create2("fhBgPairSourcesLeptonCut","fhBgPairSourcesLeptonCut;Minus Canidates;Plus Canidate;Yield", 4, -.5, 3.5, 4, -.5, 3.5); fHM->Create2("fhBgPairSourcesMomCut","fhBgPairSourcesMomCut;Minus Canidates;Plus Canidate;Yield", 4, -.5, 3.5, 4, -.5, 3.5); fHM->Create2("fhBgPairSourcesOACut","fhBgPairSourcesOACut;Minus Canidates;Plus Canidate;Yield", 4, -.5, 3.5, 4, -.5, 3.5); //Kine ID / Creation Mechanism fHM->Create1("fhKineID", "fhKineID;PID", 51, -.5, 50.5); fHM->Create1("fhKineMechanism", "fhKineMechanism; Creation Mechanism", 51, -.5, 50.5); fHM->Create1("fhParentFromGamma", "fhParentFromGamma;PID", 20, -.5, 19.5); // Kine parameters for signal, gamma and pi0 e+- fHM->Create1("fhMomKineSignalEl", "fhMomKineSignalEl;Momentum [MeV/c];Yield", 100, 0., 2000.); fHM->Create1("fhTransMomKineSignalEl", "fhTransMomKineSignalEl;Momentum [MeV/c];Yield", 100, 0., 2000.); fHM->Create1("fhMomKineGammaEl", "fhMomKineGammaEl;Momentum [MeV/c];Yield", 100, 0., 2000.); fHM->Create1("fhTransMomKineGammaEl", "fhTransMomKineGammaEl;Momentum [MeV/c];Yield", 100, 0., 2000.); fHM->Create1("fhMomKinePi0El", "fhMomKinePi0El;Momentum [MeV/c];Yield", 100, 0., 2000.); fHM->Create1("fhTransMomKinePi0El", "fhTransMomKinePi0El;Momentum [MeV/c];Yield", 100, 0., 2000.); fHM->Create1("fhOpenAngleKineSignalEl", "fhOpenAngleKineSignalEl;Opening angle [deg];Yield", 190, 0., 190.); fHM->Create1("fhOpenAngleKineGammaEl", "fhOpenAngleKineGammaEl;Opening angle [deg];Yield", 190, 0., 190.); fHM->Create1("fhOpenAngleKinePi0El", "fhOpenAngleKinePi0El;Opening angle [deg];Yield", 190, 0., 190.); fHM->Create1("fhNofKinePi0", "fhNofKinePi0;Nof particles/event;Yield", 100, -.5, 99.5); fHM->Create1("fhNofKinePi0Primary", "fhNofKinePi0Primary;Nof particles/event;Yield", 100, -.5, 99.5); // Number of candidates per event fHM->Create1("fhNofCandsAll", "fhNofCandsAll;Nof candidates/event;Yield", 100, -.5, 99.5); fHM->Create1("fhNofCandsGood", "fhNofCandsGood;Nof candidates/event;Yield", 100, -.5, 99.5); fHM->Create1("fhNofCandsLepton", "fhNofCandsLepton;Nof candidates/event;Yield", 100, -.5, 99.5); fHM->Create1("fhNofCandsGammaEl", "fhNofCandsGammaEl;Nof candidates/event;Yield", 100, -.5, 99.5); fHM->Create1("fhNofCandsPi0El", "fhNofCandsPi0El;Nof candidates/event;Yield", 100, -.5, 99.5); //Number of rings in RICH fHM->Create1("fhNofRingsAll", "fhNofRingsAll;Nof Rings/event;Yield", 15, -0.5, 14.5); // Number of hits in ring depended on Theta fHM->Create2("fhNofHitsInRingVsThetaSignalEl", "fhNofHitsInRingVsThetaSignalEl;Theta [deg];Nof hits in ring;Yield", 30, 0, 90, 25, -0.5, 24.5); fHM->Create2("fhNofHitsInRingVsThetaGammaEl", "fhNofHitsInRingVsThetaGammaEl;Theta [deg];Nof hits in ring;Yield", 30, 0, 90, 25, -0.5, 24.5); fHM->Create2("fhNofHitsInRingVsThetaPi0El", "fhNofHitsInRingVsThetaPi0El;Theta [deg];Nof hits in ring;Yield", 30, 0, 90, 25, -0.5, 24.5); fHM->Create2("fhNofHitsInRingVsThetaOther", "fhNofHitsInRingVsThetaOther;Theta [deg];Nof hits in ring;Yield", 30, 0, 90, 25, -0.5, 24.5); // Number of hits in ring depended on Theta starting with cand loop fHM->Create2("fhNofHitsInRingVsThetaSignalElCands", "fhNofHitsInRingVsThetaSignalElCands;Theta [deg];Nof hits in ring;Yield", 30, 0, 90, 25, -0.5, 24.5); fHM->Create2("fhNofHitsInRingVsThetaGammaElCands", "fhNofHitsInRingVsThetaGammaElCands;Theta [deg];Nof hits in ring;Yield", 30, 0, 90, 25, -0.5, 24.5); fHM->Create2("fhNofHitsInRingVsThetaPi0ElCands", "fhNofHitsInRingVsThetaPi0ElCands;Theta [deg];Nof hits in ring;Yield", 30, 0, 90, 25, -0.5, 24.5); fHM->Create2("fhNofHitsInRingVsThetaOtherCands", "fhNofHitsInRingVsThetaOtherCands;Theta [deg];Nof hits in ring;Yield", 30, 0, 90, 25, -0.5, 24.5); // Number of hits in ring depended on Theta starting with cand loop Lepton Cut fHM->Create2("fhNofHitsInRingVsThetaSignalElCandsLeptonCut", "fhNofHitsInRingVsThetaSignalElCandsLeptonCut;Theta [deg];Nof hits in ring;Yield", 30, 0, 90, 25, -0.5, 24.5); fHM->Create2("fhNofHitsInRingVsThetaGammaElCandsLeptonCut", "fhNofHitsInRingVsThetaGammaElCandsLeptonCut;Theta [deg];Nof hits in ring;Yield", 30, 0, 90, 25, -0.5, 24.5); fHM->Create2("fhNofHitsInRingVsThetaPi0ElCandsLeptonCut", "fhNofHitsInRingVsThetaPi0ElCandsLeptonCut;Theta [deg];Nof hits in ring;Yield", 30, 0, 90, 25, -0.5, 24.5); fHM->Create2("fhNofHitsInRingVsThetaOtherCandsLeptonCut", "fhNofHitsInRingVsThetaOtherCandsLeptonCut;Theta [deg];Nof hits in ring;Yield", 30, 0, 90, 25, -0.5, 24.5); // Number of hits in ring depended on Theta for shared rings fHM->Create2("fhNofHitsInRingVsThetaSignalElShared", "fhNofHitsInRingVsThetaSignalElShared;Theta [deg];Nof hits in ring;Yield", 30, 0, 90, 25, -0.5, 24.5); fHM->Create2("fhNofHitsInRingVsThetaGammaElShared", "fhNofHitsInRingVsThetaGammaElShared;Theta [deg];Nof hits in ring;Yield", 30, 0, 90, 25, -0.5, 24.5); fHM->Create2("fhNofHitsInRingVsThetaPi0ElShared", "fhNofHitsInRingVsThetaPi0ElShared;Theta [deg];Nof hits in ring;Yield", 30, 0, 90, 25, -0.5, 24.5); fHM->Create2("fhNofHitsInRingVsThetaOtherShared", "fhNofHitsInRingVsThetaOtherShared;Theta [deg];Nof hits in ring;Yield", 30, 0, 90, 25, -0.5, 24.5); // Number of hits in ring depended on Theta for shared rings fHM->Create2("fhNofHitsInRingVsThetaSignalElNotShared", "fhNofHitsInRingVsThetaSignalElNotShared;Theta [deg];Nof hits in ring;Yield", 30, 0, 90, 25, -0.5, 24.5); fHM->Create2("fhNofHitsInRingVsThetaGammaElNotShared", "fhNofHitsInRingVsThetaGammaElNotShared;Theta [deg];Nof hits in ring;Yield", 30, 0, 90, 25, -0.5, 24.5); fHM->Create2("fhNofHitsInRingVsThetaPi0ElNotShared", "fhNofHitsInRingVsThetaPi0ElNotShared;Theta [deg];Nof hits in ring;Yield", 30, 0, 90, 25, -0.5, 24.5); fHM->Create2("fhNofHitsInRingVsThetaOtherNotShared", "fhNofHitsInRingVsThetaOtherNotShared;Theta [deg];Nof hits in ring;Yield", 30, 0, 90, 25, -0.5, 24.5); // Number of hits in ring vs nof assigned candidates fHM->Create1("fhNofHitsInRingNotShared", "fhNofHitsInRingNotShared;Nof hits in ring;Yield", 25, -0.5, 24.5); fHM->Create1("fhNofHitsInRingNotSharedSignalEl", "fhNofHitsInRingNotSharedSignalEl;Nof hits in ring;Yield", 25, -0.5, 24.5); fHM->Create1("fhNofHitsInRingShared", "fhNofHitsInRingShared;Nof hits in ring;Yield", 25, -0.5, 24.5); fHM->Create1("fhNofHitsInRingNotAssigned", "fhNofHitsInRingNotAssigned;Nof hits in ring;Yield", 25, -0.5, 24.5); // Nof shared particle candidates per ring fHM->Create1("fhNofSharedCandidatesForRing", "fhNofSharedCandidatesForRing;Nof shared particle candidates/ring;Yield", 10, -0.5, 9.5); //Charge * momentum fHM->Create1("fhQPSignalEl", "fhQPSignalEl;Momentum*Charge [C*MeV/c];Yield", 200, -2000., 2000.); fHM->Create1("fhQPGammaEl", "fhQPGammaEl;Momentum*Charge [C*MeV/c];Yield", 200, -2000., 2000.); fHM->Create1("fhQPPi0El", "fhQPPi0El;Momentum*Charge [C*MeV/c];Yield", 200, -2000., 2000.); fHM->Create1("fhQPPion", "fhQPPion;Momentum*Charge [C*MeV/c];Yield", 200, -2000., 2000.); fHM->Create1("fhQPProton", "fhQPProton;Momentum*Charge [C*MeV/c];Yield", 200, -2000., 2000.); fHM->Create1("fhQPOther", "fhQPOther;Momentum*Charge [C*MeV/c];Yield", 200, -2000., 2000.); fHM->Create1("fhQPLeptonCutSignalEl", "fhQPLeptonCutSignalEl;Momentum*Charge [C*MeV/c];Yield", 200, -2000., 2000.); fHM->Create1("fhQPLeptonCutGammaEl", "fhQPLeptonCutGammaEl;Momentum*Charge [C*MeV/c];Yield", 200, -2000., 2000.); fHM->Create1("fhQPLeptonCutPi0El", "fhQPLeptonCutPi0El;Momentum*Charge [C*MeV/c];Yield", 200, -2000., 2000.); fHM->Create1("fhQPLeptonCutPion", "fhQPLeptonCutPion;Momentum*Charge [C*MeV/c];Yield", 200, -2000., 2000.); fHM->Create1("fhQPLeptonCutProton", "fhQPLeptonCutProton;Momentum*Charge [C*MeV/c];Yield", 200, -2000., 2000.); fHM->Create1("fhQPLeptonCutOther", "fhQPLeptonCutOther;Momentum*Charge [C*MeV/c];Yield", 200, -2000., 2000.); fHM->Create1("fhQPMomCutSignalEl", "fhQPMomCutSignalEl;Momentum*Charge [C*MeV/c];Yield", 200, -2000., 2000.); fHM->Create1("fhQPMomCutGammaEl", "fhQPMomCutGammaEl;Momentum*Charge [C*MeV/c];Yield", 200, -2000., 2000.); fHM->Create1("fhQPMomCutPi0El", "fhQPMomCutPi0El;Momentum*Charge [C*MeV/c];Yield", 200, -2000., 2000.); fHM->Create1("fhQPMomCutPion", "fhQPMomCutPion;Momentum*Charge [C*MeV/c];Yield", 200, -2000., 2000.); fHM->Create1("fhQPMomCutProton", "fhQPMomCutProton;Momentum*Charge [C*MeV/c];Yield", 200, -2000., 2000.); fHM->Create1("fhQPMomCutOther", "fhQPMomCutOther;Momentum*Charge [C*MeV/c];Yield", 200, -2000., 2000.); fHM->Create1("fhQPOACutSignalEl", "fhQPOACutSignalEl;Momentum*Charge [C*MeV/c];Yield", 200, -2000., 2000.); fHM->Create1("fhQPOACutGammaEl", "fhQPOACutGammaEl;Momentum*Charge [C*MeV/c];Yield", 200, -2000., 2000.); fHM->Create1("fhQPOACutPi0El", "fhQPOACutPi0El;Momentum*Charge [C*MeV/c];Yield", 200, -2000., 2000.); fHM->Create1("fhQPOACutPion", "fhQPOACutPion;Momentum*Charge [C*MeV/c];Yield", 200, -2000., 2000.); fHM->Create1("fhQPOACutProton", "fhQPOACutProton;Momentum*Charge [C*MeV/c];Yield", 200, -2000., 2000.); fHM->Create1("fhQPOACutOther", "fhQPOACutOther;Momentum*Charge [C*MeV/c];Yield", 200, -2000., 2000.); //Beta vs momentum fHM->Create2("fhBetaVsMomSignalEl", "fhBetaVsMomSignalEl;Momentum [MeV/c];#beta;Yield", 30., 0., 2000., 100, 0.7, 1.2); fHM->Create2("fhBetaVsMomGammaEl", "fhBetaVsMomGammaEl;Momentum [MeV/c];#beta;Yield", 30., 0., 2000., 100, 0.7, 1.2); fHM->Create2("fhBetaVsMomPi0El", "fhBetaVsMomPi0El;Momentum [MeV/c];#beta;Yield", 30., 0., 2000., 100, 0.7, 1.2); fHM->Create2("fhBetaVsMomPion", "fhBetaVsMomPion;Momentum [MeV/c];#beta;Yield", 30., 0., 2000., 100, 0.7, 1.2); fHM->Create2("fhBetaVsMomProton", "fhBetaVsMomProton;Momentum [MeV/c];#beta;Yield", 30., 0., 2000., 100, 0.7, 1.2); fHM->Create2("fhBetaVsMomOther", "fhBetaVsMomOther;Momentum [MeV/c];#beta;Yield", 30., 0., 2000., 100, 0.7, 1.2); fHM->Create2("fhBetaVsMomLeptonCutSignalEl", "fhBetaVsMomLeptonCutSignalEl;Momentum [MeV/c];#beta;Yield", 30., 0., 2000., 100, 0.7, 1.2); fHM->Create2("fhBetaVsMomLeptonCutGammaEl", "fhBetaVsMomLeptonCutGammaEl;Momentum [MeV/c];#beta;Yield", 30., 0., 2000., 100, 0.7, 1.2); fHM->Create2("fhBetaVsMomLeptonCutPi0El", "fhBetaVsMomLeptonCutPi0El;Momentum [MeV/c];#beta;Yield", 30., 0., 2000., 100, 0.7, 1.2); fHM->Create2("fhBetaVsMomLeptonCutPion", "fhBetaVsMomLeptonCutPion;Momentum [MeV/c];#beta;Yield", 30., 0., 2000., 100, 0.7, 1.2); fHM->Create2("fhBetaVsMomLeptonCutProton", "fhBetaVsMomLeptonCutProton;Momentum [MeV/c];#beta;Yield", 30., 0., 2000., 100, 0.7, 1.2); fHM->Create2("fhBetaVsMomLeptonCutOther", "fhBetaVsMomLeptonCutOther;Momentum [MeV/c];#beta;Yield", 30., 0., 2000., 100, 0.7, 1.2); fHM->Create2("fhBetaVsMomMomCutSignalEl", "fhBetaVsMomMomCutSignalEl;Momentum [MeV/c];#beta;Yield", 30., 0., 2000., 100, 0.7, 1.2); fHM->Create2("fhBetaVsMomMomCutGammaEl", "fhBetaVsMomMomCutGammaEl;Momentum [MeV/c];#beta;Yield", 30., 0., 2000., 100, 0.7, 1.2); fHM->Create2("fhBetaVsMomMomCutPi0El", "fhBetaVsMomMomCutPi0El;Momentum [MeV/c];#beta;Yield", 30., 0., 2000., 100, 0.7, 1.2); fHM->Create2("fhBetaVsMomMomCutPion", "fhBetaVsMomMomCutPion;Momentum [MeV/c];#beta;Yield", 30., 0., 2000., 100, 0.7, 1.2); fHM->Create2("fhBetaVsMomMomCutProton", "fhBetaVsMomMomCutProton;Momentum [MeV/c];#beta;Yield", 30., 0., 2000., 100, 0.7, 1.2); fHM->Create2("fhBetaVsMomMomCutOther", "fhBetaVsMomMomCutOther;Momentum [MeV/c];#beta;Yield", 30., 0., 2000., 100, 0.7, 1.2); fHM->Create2("fhBetaVsMomOACutSignalEl", "fhBetaVsMomOACutSignalEl;Momentum [MeV/c];#beta;Yield", 30., 0., 2000., 100, 0.7, 1.2); fHM->Create2("fhBetaVsMomOACutGammaEl", "fhBetaVsMomOACutGammaEl;Momentum [MeV/c];#beta;Yield", 30., 0., 2000., 100, 0.7, 1.2); fHM->Create2("fhBetaVsMomOACutPi0El", "fhBetaVsMomOACutPi0El;Momentum [MeV/c];#beta;Yield", 30., 0., 2000., 100, 0.7, 1.2); fHM->Create2("fhBetaVsMomOACutPion", "fhBetaVsMomOACutPion;Momentum [MeV/c];#beta;Yield", 30., 0., 2000., 100, 0.7, 1.2); fHM->Create2("fhBetaVsMomOACutProton", "fhBetaVsMomOACutProton;Momentum [MeV/c];#beta;Yield", 30., 0., 2000., 100, 0.7, 1.2); fHM->Create2("fhBetaVsMomOACutOther", "fhBetaVsMomOACutOther;Momentum [MeV/c];#beta;Yield", 30., 0., 2000., 100, 0.7, 1.2); //dEdX vs momentum fHM->Create2("fhdEdXVsMomSignalEl", "fhdEdXVsMomSignalEl;Momentum [MeV/c];dEdx [MeV/cm];Yield", 30., 0., 2000., 100, 0., 10.); fHM->Create2("fhdEdXVsMomGammaEl", "fhdEdXVsMomGammaEl;Momentum [MeV/c];dEdx [MeV/cm];Yield", 30., 0., 2000., 100, 0., 10.); fHM->Create2("fhdEdXVsMomPi0El", "fhdEdXVsMomPi0El;Momentum [MeV/c];dEdx [MeV/cm];Yield", 30., 0., 2000., 100, 0., 10.); fHM->Create2("fhdEdXVsMomPion", "fhdEdXVsMomPion;Momentum [MeV/c];dEdx [MeV/cm];Yield", 30., 0., 2000., 100, 0., 10.); fHM->Create2("fhdEdXVsMomProton", "fhdEdXVsMomProton;Momentum [MeV/c];dEdx [MeV/cm];Yield", 30., 0., 2000., 100, 0., 10.); fHM->Create2("fhdEdXVsMomOther", "fhdEdXVsMomOther;Momentum [MeV/c];dEdx [MeV/cm];Yield", 30., 0., 2000., 100, 0., 10.); fHM->Create2("fhdEdXVsMomLeptonCutSignalEl", "fhdEdXVsMomLeptonCutSignalEl;Momentum [MeV/c];dEdx [MeV/cm];Yield", 30., 0., 2000., 100, 0., 10.); fHM->Create2("fhdEdXVsMomLeptonCutGammaEl", "fhdEdXVsMomLeptonCutGammaEl;Momentum [MeV/c];dEdx [MeV/cm];Yield", 30., 0., 2000., 100, 0., 10.); fHM->Create2("fhdEdXVsMomLeptonCutPi0El", "fhdEdXVsMomLeptonCutPi0El;Momentum [MeV/c];dEdx [MeV/cm];Yield", 30., 0., 2000., 100, 0., 10.); fHM->Create2("fhdEdXVsMomLeptonCutPion", "fhdEdXVsMomLeptonCutPion;Momentum [MeV/c];dEdx [MeV/cm];Yield", 30., 0., 2000., 100, 0., 10.); fHM->Create2("fhdEdXVsMomLeptonCutProton", "fhdEdXVsMomLeptonCutProton;Momentum [MeV/c];dEdx [MeV/cm];Yield", 30., 0., 2000., 100, 0., 10.); fHM->Create2("fhdEdXVsMomLeptonCutOther", "fhdEdXVsMomLeptonCutOther;Momentum [MeV/c];dEdx [MeV/cm];Yield", 30., 0., 2000., 100, 0., 10.); fHM->Create2("fhdEdXVsMomMomCutSignalEl", "fhdEdXVsMomMomCutSignalEl;Momentum [MeV/c];dEdx [MeV/cm];Yield", 30., 0., 2000., 100, 0., 10.); fHM->Create2("fhdEdXVsMomMomCutGammaEl", "fhdEdXVsMomMomCutGammaEl;Momentum [MeV/c];dEdx [MeV/cm];Yield", 30., 0., 2000., 100, 0., 10.); fHM->Create2("fhdEdXVsMomMomCutPi0El", "fhdEdXVsMomMomCutPi0El;Momentum [MeV/c];dEdx [MeV/cm];Yield", 30., 0., 2000., 100, 0., 10.); fHM->Create2("fhdEdXVsMomMomCutPion", "fhdEdXVsMomMomCutPion;Momentum [MeV/c];dEdx [MeV/cm];Yield", 30., 0., 2000., 100, 0., 10.); fHM->Create2("fhdEdXVsMomMomCutProton", "fhdEdXVsMomMomCutProton;Momentum [MeV/c];dEdx [MeV/cm];Yield", 30., 0., 2000., 100, 0., 10.); fHM->Create2("fhdEdXVsMomMomCutOther", "fhdEdXVsMomMomCutOther;Momentum [MeV/c];dEdx [MeV/cm];Yield", 30., 0., 2000., 100, 0., 10.); fHM->Create2("fhdEdXVsMomOACutSignalEl", "fhdEdXVsMomOACutSignalEl;Momentum [MeV/c];dEdx [MeV/cm];Yield", 30., 0., 2000., 100, 0., 10.); fHM->Create2("fhdEdXVsMomOACutGammaEl", "fhdEdXVsMomOACutGammaEl;Momentum [MeV/c];dEdx [MeV/cm];Yield", 30., 0., 2000., 100, 0., 10.); fHM->Create2("fhdEdXVsMomOACutPi0El", "fhdEdXVsMomOACutPi0El;Momentum [MeV/c];dEdx [MeV/cm];Yield", 30., 0., 2000., 100, 0., 10.); fHM->Create2("fhdEdXVsMomOACutPion", "fhdEdXVsMomOACutPion;Momentum [MeV/c];dEdx [MeV/cm];Yield", 30., 0., 2000., 100, 0., 10.); fHM->Create2("fhdEdXVsMomOACutProton", "fhdEdXVsMomOACutProton;Momentum [MeV/c];dEdx [MeV/cm];Yield", 30., 0., 2000., 100, 0., 10.); fHM->Create2("fhdEdXVsMomOACutOther", "fhdEdXVsMomOACutOther;Momentum [MeV/c];dEdx [MeV/cm];Yield", 30., 0., 2000., 100, 0., 10.); // Geant PID of particle candidates before/aftercuts fHM->Create1("fhGeantPidAll", "fhGeantPidAll;GEANT PID", 51, -.5, 50.5); fHM->Create1("fhGeantPidLeptonCutAll", "fhGeantPidLeptonCutAll;GEANT PID", 51, -.5, 50.5); fHM->Create1("fhGeantPidMomCutAll", "fhGeantPidMomCutAll;GEANT PID", 51, -.5, 50.5); fHM->Create1("fhGeantPidOACutAll", "fhGeantPidOACutAll;GEANT PID", 51, -.5, 50.5); fHM->Create1("fhGeantPidSignal", "fhGeantPidSignal;GEANT PID", 51, -.5, 50.5); fHM->Create1("fhGeantPidLeptonCutSignal", "fhGeantPidLeptonCutSignal;GEANT PID", 51, -.5, 50.5); fHM->Create1("fhGeantPidMomCutSignal", "fhGeantPidMomCutSignal;GEANT PID", 51, -.5, 50.5); fHM->Create1("fhGeantPidOACutSignal", "fhGeantPidOACutSignal;GEANT PID", 51, -.5, 50.5); fHM->Create1("fhGeantPidBg", "fhGeantPidBg;GEANT PID", 51, -.5, 50.5); fHM->Create1("fhGeantPidLeptonCutBg", "fhGeantPidLeptonCutBg;GEANT PID", 51, -.5, 50.5); fHM->Create1("fhGeantPidMomCutBg", "fhGeantPidMomCutBg;GEANT PID", 51, -.5, 50.5); fHM->Create1("fhGeantPidOACutBg", "fhGeantPidOACutBg;GEANT PID", 51, -.5, 50.5); // particle candidate in RICH fHM->Create2("fhRingXYAll", "fhRingXYAll;X [mm];Y [mm]", 100, -600., 600., 100, -600., 600.); fHM->Create2("fhRingXYSignalEl", "fhRingXYSignalEl;X [mm];Y [mm]", 100, -600., 600., 100, -600., 600.); fHM->Create2("fhRingXYGammaEl", "fhRingXYGammaEl;X [mm];Y [mm]", 100, -600., 600., 100, -600., 600.); fHM->Create2("fhRingXYPi0El", "fhRingXYPi0El;X [mm];Y [mm]", 100, -600., 600., 100, -600., 600.); fHM->Create2("fhRingXYLeptonCut", "fhRingXYLeptonCut;X [mm];Y [mm]", 100, -600., 600., 100, -600., 600.); fHM->Create2("fhRadiusVsNofHitsAll", "fhRadiusVsNofHitsAll;Nof hits in ring;Radius [mm];Yield", 30, -0.5, 29.5, 70, 0., 35.); fHM->Create2("fhRadiusVsNofHitsSignalEl", "fhRadiusVsNofHitsSignalEl;Nof hits in ring;Radius [mm];Yield", 30, -0.5, 29.5, 70, 0., 35.); fHM->Create2("fhRadiusVsNofHitsGammaEl", "fhRadiusVsNofHitsGammaEl;Nof hits in ring;Radius [mm];Yield", 30, -0.5, 29.5, 70, 0., 35.); fHM->Create2("fhRadiusVsNofHitsPi0El", "fhRadiusVsNofHitsPi0El;Nof hits in ring;Radius [mm];Yield", 30,-0.5, 29.5, 70, 0., 35.); fHM->Create2("fhRadiusVsMomAll", "fhRadiusVsMomAll;Momentum [MeV/c];Radius [mm];Yield", 30, 0., 1500., 70, 0., 35.); fHM->Create2("fhRadiusVsMomSignalEl", "fhRadiusVsMomSignalEl;Momentum [MeV/c];Radius [mm];Yield", 30, 0., 1500., 70, 0., 35.); fHM->Create2("fhRadiusVsMomGammaEl", "fhRadiusVsMomGammaEl;Momentum [MeV/c];Radius [mm];Yield", 30, 0., 1500., 70, 0., 35.); fHM->Create2("fhRadiusVsMomPi0El", "fhRadiusVsMomPi0El;Momentum [MeV/c];Radius [mm];Yield", 30, 0., 1500., 70, 0., 35.); fHM->Create1("fhOpenAngleSignal", "fhOpenAngleSignal;Opening angle [deg];Yield", 190, 0., 190.); fHM->Create1("fhOpenAngleGamma", "fhOpenAngleGamma;Opening angle [deg];Yield", 190, 0., 190.); fHM->Create1("fhOpenAnglePi0", "fhOpenAnglePi0;Opening angle [deg];Yield", 190, 0., 190.); fHM->Create1("fhOpenAngleBg", "fhOpenAngleBg;Opening angle [deg];Yield", 190, 0., 190.); fHM->Create1("fhOpenAngleSignalPositronToBG", "fhOpenAngleSignalPositronToBG;Opening angle [deg];Yield", 190, 0., 190.); fHM->Create1("fhOpenAngleSignalPositronToBGMomCut", "fhOpenAngleSignalPositronToBGMomCut;Opening angle [deg];Yield", 190, 0., 190.); fHM->Create1("fhOpenAngleSignalPositronToBGLeptonCut", "fhOpenAngleSignalPositronToBGLeptonCut;Opening angle [deg];Yield", 190, 0., 190.); fHM->Create1("fhOpenAngleSignalPositronToBGOACut", "fhOpenAngleSignalPositronToBGOACut;Opening angle [deg];Yield", 190, 0., 190.); fHM->Create1("fhOpenAngleSignalElectronToBG", "fhOpenAngleSignalElectronToBG;Opening angle [deg];Yield", 190, 0., 190.); fHM->Create1("fhOpenAngleSignalElectronToBGMomCut", "fhOpenAngleSignalElectronToBGMomCut;Opening angle [deg];Yield", 190, 0., 190.); fHM->Create1("fhOpenAngleSignalElectronToBGLeptonCut", "fhOpenAngleSignalElectronToBGLeptonCut;Opening angle [deg];Yield", 190, 0., 190.); fHM->Create1("fhOpenAngleSignalElectronToBGOACut", "fhOpenAngleSignalElectronToBGOACut;Opening angle [deg];Yield", 190, 0., 190.); fHM->Create2("fhOpenAngleVsRingDistBg", "fhOpenAngleVsRingDistBg;Opening angle [deg];Ring distance [mm]", 150, 0., 150., 100, 0., 1000.); fHM->Create2("fhOpenAngleVsRingDistSignal", "fhOpenAngleVsRingDistSignal;Opening angle [deg];Ring distance [mm]", 190, 0., 190., 150, 0., 1500.); fHM->Create2("fhOpenAngleVsRingDistGamma", "fhOpenAngleVsRingDistGamma;Opening angle [deg];Ring distance [mm]", 40, 0., 20., 50, 0., 50.); fHM->Create2("fhOpenAngleVsRingDistPi0", "fhOpenAngleVsRingDistPi0;Opening angle [deg];Ring distance [mm]", 40, 0., 40., 150, 0., 150.); // Distance between 2 ring centers fHM->Create2("fhRingDistVsMomPairSignal","fhRingDistVsMomPairSignal;Momentum_{pair} [MeV/c];Distance_{Ring1-Ring2} [mm]", 30, 0., 2000., 100, 0., 2000.); fHM->Create2("fhRingDistVsMomPairGamma","fhRingDistVsMomPairGamma;Momentum_{pair} [MeV/c];Distance_{Ring1-Ring2} [mm]", 30, 0., 2000., 100, 0., 2000.); fHM->Create2("fhRingDistVsMomPairPi0","fhRingDistVsMomPairPi0;Momentum_{pair} [MeV/c];Distance_{Ring1-Ring2} [mm]", 30, 0., 2000., 100, 0., 2000.); fHM->Create2("fhRingDistVsMomPairBg","fhRingDistVsMomPairBg;Momentum_{pair} [MeV/c];Distance_{Ring1-Ring2} [mm]", 30, 0., 2000., 100, 0., 2000.); fHM->Create2("fhMinDistVsMinAngle", "fhMinDistVsMinAngle;Ring distance [mm];Opening angle [deg]", 150, 0., 150., 40, 0., 40.); // X axis : all, Pair, No pair // Y axis : pi0 e-, pi0 e+, gamma e-, gamma e+ fHM->Create2("fhPairRecoTopology", "fhPairRecoTopology;Pair type;Candidate;Yield", 3, -.5, 2.5, 4, -.5, 3.5); fHM->Create2("fhPairRecoTopologyNoLeptonCut", "fhPairRecoTopologyNoLeptonCut;Pair type;Candidate;Yield", 3, -.5, 2.5, 4, -.5, 3.5); } Bool_t HDEAnaTask::reinit() { return kTRUE; } Int_t HDEAnaTask::execute() { fHM->H1("fhNofEvents")->Fill(0); fEventNum++; cout << "HDEAnalysisTask::execute eventNum " << fEventNum << endl; fillParticleCandHists(); fillGeantPidHists(); fillQPHists(); fillRichHists(); fillOpenAngleToBackgroundHists(); fillNofCandsHists(); fillDecayHists(); fillNofRingsHists(); fillKineHists(); return 0; } void HDEAnaTask::fillNofCandsHists() { Int_t nofCands = fCatParticleCand->getEntries(); Int_t nofCandsLepton = 0; Int_t nofCandsPi0El = 0; Int_t nofCandsGammaEl = 0; Int_t nofCandsGood = 0; fHM->H1("fhNofCandsAll")->Fill(nofCands); for(Int_t i=0; i< nofCands; i++){ HParticleCandSim* cand = static_cast(fCatParticleCand->getObject(i)); Bool_t isLepton = fCuts.isLeptonOk(cand); Bool_t isGoodLepton = cand->isFlagBit(kIsAcceptedHitRICHMDC); if (isLepton) { nofCandsLepton++; if(isGoodLepton) { nofCandsGood++; } if(HDEAnaMyUtils::isMcPi0Electron(cand, fCatKine)) { nofCandsPi0El++; } if(HDEAnaMyUtils::isMcGammaElectron(cand, fCatKine)) { nofCandsGammaEl++; } } } fHM->H1("fhNofCandsLepton")->Fill(nofCandsLepton); fHM->H1("fhNofCandsPi0El")->Fill(nofCandsPi0El); fHM->H1("fhNofCandsGammaEl")->Fill(nofCandsGammaEl); fHM->H1("fhNofCandsGood")->Fill(nofCandsGood); } void HDEAnaTask::fillDecayHists() { Int_t nofCands = fCatParticleCand->getEntries(); for(Int_t i1 = 0; i1 < nofCands; i1++){ HParticleCandSim* cand1 = static_cast(fCatParticleCand->getObject(i1)); if ( !fCuts.isLeptonOk(cand1)) continue; // Y axis : pi0 e-, pi0 e+, gamma e-, gamma e+ Double_t yInd = -1.; bool isPi0El = HDEAnaMyUtils::isMcPi0Electron(cand1, fCatKine); bool isGammaEl = HDEAnaMyUtils::isMcGammaElectron(cand1, fCatKine); if (isPi0El && cand1->getCharge() < 0) { yInd = 0.; } else if (isPi0El && cand1->getCharge() > 0) { yInd = 1.; } else if (isGammaEl && cand1->getCharge() < 0) { yInd = 2.; } else if (isGammaEl && cand1->getCharge() > 0) { yInd = 3.; } if (yInd < 0.) continue; Bool_t hasPair = false; Bool_t hasPairNoLeptonCut = false; fHM->H2("fhPairRecoTopology")->Fill(0., yInd); fHM->H2("fhPairRecoTopologyNoLeptonCut")->Fill(0., yInd); for(Int_t i2 = 0; i2 < nofCands; i2++){ HParticleCandSim* cand2 = static_cast(fCatParticleCand->getObject(i2)); if (i1 == i2) continue; if ( fCuts.isLeptonOk(cand2) && (yInd == 0. || yInd == 1.) && HDEAnaMyUtils::isMcPi0Pair(cand1, cand2, fCatKine)) { hasPair = true; } if ( fCuts.isLeptonOk(cand2) && (yInd == 2. || yInd == 3.) && HDEAnaMyUtils::isMcGammaPair(cand1, cand2, fCatKine)) { hasPair = true; } // No lepton cut if ( (yInd == 0. || yInd == 1.) && HDEAnaMyUtils::isMcPi0Pair(cand1, cand2, fCatKine)) { hasPairNoLeptonCut = true; } if ( (yInd == 2. || yInd == 3.) && HDEAnaMyUtils::isMcGammaPair(cand1, cand2, fCatKine)) { hasPairNoLeptonCut = true; } } if (hasPair) { fHM->H2("fhPairRecoTopology")->Fill(1., yInd); } else { fHM->H2("fhPairRecoTopology")->Fill(2., yInd); } if (hasPairNoLeptonCut) { fHM->H2("fhPairRecoTopologyNoLeptonCut")->Fill(1., yInd); } else { fHM->H2("fhPairRecoTopologyNoLeptonCut")->Fill(2., yInd); } } } void HDEAnaTask::fillNofRingsHists() { Int_t nofRings = fCatRichHit->getEntries(); fHM->H1("fhNofRingsAll")->Fill(nofRings); for (Int_t iR = 0; iR < nofRings; iR++) { HRichHitSim* ring = static_cast(fCatRichHit->getObject(iR)); if (ring == NULL) continue; Int_t trackId = ring->track1; Int_t nofHits = ring->fRich700NofRichCals; Float_t theta = ring->getTheta(); Int_t nofSharedCands = HDEAnaMyUtils::getNofSharedCandidatesForRing(iR, fCatParticleCand); fHM->H1("fhNofSharedCandidatesForRing")->Fill(nofSharedCands); HGeantKine* kine = (HGeantKine*)fCatKine->getObject(trackId - 1); if (nofSharedCands == 1) { fHM->H1("fhNofHitsInRingNotShared")->Fill(nofHits); if (HDEAnaMyUtils::isMcSignalElectron(kine)) { fHM->H2("fhNofHitsInRingVsThetaSignalElNotShared")->Fill( theta, nofHits); fHM->H1("fhNofHitsInRingNotSharedSignalEl")->Fill(nofHits); } else if (HDEAnaMyUtils::isMcGammaElectron(kine, fCatKine)) { fHM->H2("fhNofHitsInRingVsThetaGammaElNotShared")->Fill(theta, nofHits); } else if (HDEAnaMyUtils::isMcPi0Electron(kine, fCatKine)) { fHM->H2("fhNofHitsInRingVsThetaPi0ElNotShared")->Fill(theta, nofHits); } else { fHM->H2("fhNofHitsInRingVsThetaOtherNotShared")->Fill(theta, nofHits); } } else if (nofSharedCands == 0){ fHM->H1("fhNofHitsInRingNotAssigned")->Fill(nofHits); } else { fHM->H1("fhNofHitsInRingShared")->Fill(nofHits); if (HDEAnaMyUtils::isMcSignalElectron(kine)) { fHM->H2("fhNofHitsInRingVsThetaSignalElShared")->Fill( theta, nofHits); } else if (HDEAnaMyUtils::isMcGammaElectron(kine, fCatKine)) { fHM->H2("fhNofHitsInRingVsThetaGammaElShared")->Fill(theta, nofHits); } else if (HDEAnaMyUtils::isMcPi0Electron(kine, fCatKine)) { fHM->H2("fhNofHitsInRingVsThetaPi0ElShared")->Fill(theta, nofHits); } else { fHM->H2("fhNofHitsInRingVsThetaOtherShared")->Fill(theta, nofHits); } } if (HDEAnaMyUtils::isMcSignalElectron(kine)) { fHM->H2("fhNofHitsInRingVsThetaSignalEl")->Fill( theta, nofHits); } else if (HDEAnaMyUtils::isMcGammaElectron(kine, fCatKine)) { fHM->H2("fhNofHitsInRingVsThetaGammaEl")->Fill(theta, nofHits); } else if (HDEAnaMyUtils::isMcPi0Electron(kine, fCatKine)) { fHM->H2("fhNofHitsInRingVsThetaPi0El")->Fill(theta, nofHits); } else { fHM->H2("fhNofHitsInRingVsThetaOther")->Fill(theta, nofHits); } } } void HDEAnaTask::fillGeantPidHists() { Int_t nofCands = fCatParticleCand->getEntries(); for( Int_t iC = 0; iC < nofCands; iC++){ HParticleCandSim* cand = static_cast(fCatParticleCand->getObject(iC)); Bool_t isLepton = fCuts.isLeptonOk(cand); Bool_t isMomCutOk = fCuts.isMomCutOk(cand); Bool_t isOACutOk = fCuts.isOACutOk(iC); Int_t geantPid = cand->getGeantPID(); fHM->H1("fhGeantPidAll")->Fill(geantPid); if (HDEAnaMyUtils::isMcSignalElectron(cand, fCatKine)) { fHM->H1("fhGeantPidSignal")->Fill(geantPid); } else { fHM->H1("fhGeantPidBg")->Fill(geantPid); } if (isLepton) { fHM->H1("fhGeantPidLeptonCutAll")->Fill(geantPid); if (HDEAnaMyUtils::isMcSignalElectron(cand, fCatKine)) { fHM->H1("fhGeantPidLeptonCutSignal")->Fill(geantPid); } else { fHM->H1("fhGeantPidLeptonCutBg")->Fill(geantPid); } } if (isLepton && isMomCutOk) { fHM->H1("fhGeantPidMomCutAll")->Fill(geantPid); if (HDEAnaMyUtils::isMcSignalElectron(cand, fCatKine)) { fHM->H1("fhGeantPidMomCutSignal")->Fill(geantPid); } else { fHM->H1("fhGeantPidMomCutBg")->Fill(geantPid); } } if (isLepton && isMomCutOk && isOACutOk) { fHM->H1("fhGeantPidOACutAll")->Fill(geantPid); if (HDEAnaMyUtils::isMcSignalElectron(cand, fCatKine)) { fHM->H1("fhGeantPidOACutSignal")->Fill(geantPid); } else { fHM->H1("fhGeantPidOACutBg")->Fill(geantPid); } } } } void HDEAnaTask::fillQPHists() { Int_t nofCands = fCatParticleCand->getEntries(); for( Int_t iC = 0; iC < nofCands; iC++){ HParticleCandSim* cand = static_cast(fCatParticleCand->getObject(iC)); Bool_t isLepton = fCuts.isLeptonOk(cand); Bool_t isMomCutOk = fCuts.isMomCutOk(cand); Bool_t isOACutOk = fCuts.isOACutOk(iC); Double_t dEdX = cand->getMdcdEdx(); Double_t mom = cand->getMomentum(); Double_t beta = cand->getBeta(); Double_t pq = cand->getCharge() * mom; if (HDEAnaMyUtils::isMcSignalElectron(cand, fCatKine)) { fHM->H1("fhQPSignalEl")->Fill(pq, fSignalMBR); fHM->H2("fhBetaVsMomSignalEl")->Fill(mom, beta, fSignalMBR); fHM->H2("fhdEdXVsMomSignalEl")->Fill(mom, dEdX, fSignalMBR); } else if (HDEAnaMyUtils::isMcGammaElectron(cand, fCatKine)) { fHM->H1("fhQPGammaEl")->Fill(pq); fHM->H2("fhBetaVsMomGammaEl")->Fill(mom, beta); fHM->H2("fhdEdXVsMomGammaEl")->Fill(mom, dEdX); } else if (HDEAnaMyUtils::isMcPi0Electron(cand, fCatKine)) { fHM->H1("fhQPPi0El")->Fill(pq); fHM->H2("fhBetaVsMomPi0El")->Fill(mom, beta); fHM->H2("fhdEdXVsMomPi0El")->Fill(mom, dEdX); } else if (HDEAnaMyUtils::isMcProton(cand)) { fHM->H1("fhQPProton")->Fill(pq); fHM->H2("fhBetaVsMomProton")->Fill(mom, beta); fHM->H2("fhdEdXVsMomProton")->Fill(mom, dEdX); } else if (HDEAnaMyUtils::isMcPion(cand)) { fHM->H1("fhQPPion")->Fill(pq); fHM->H2("fhBetaVsMomPion")->Fill(mom, beta); fHM->H2("fhdEdXVsMomPion")->Fill(mom, dEdX); } else { fHM->H1("fhQPOther")->Fill(pq); fHM->H2("fhBetaVsMomOther")->Fill(mom, beta); fHM->H2("fhdEdXVsMomOther")->Fill(mom, dEdX); } // Fill histogramms after Lepton cut if (isLepton) { if (HDEAnaMyUtils::isMcSignalElectron(cand, fCatKine)) { fHM->H1("fhQPLeptonCutSignalEl")->Fill(pq, fSignalMBR); fHM->H2("fhBetaVsMomLeptonCutSignalEl")->Fill(mom, beta, fSignalMBR); fHM->H2("fhdEdXVsMomLeptonCutSignalEl")->Fill(mom, dEdX, fSignalMBR); } else if (HDEAnaMyUtils::isMcGammaElectron(cand, fCatKine)) { fHM->H1("fhQPLeptonCutGammaEl")->Fill(pq); fHM->H2("fhBetaVsMomLeptonCutGammaEl")->Fill(mom, beta); fHM->H2("fhdEdXVsMomLeptonCutGammaEl")->Fill(mom, dEdX); } else if (HDEAnaMyUtils::isMcPi0Electron(cand, fCatKine)) { fHM->H1("fhQPLeptonCutPi0El")->Fill(pq); fHM->H2("fhBetaVsMomLeptonCutPi0El")->Fill(mom, beta); fHM->H2("fhdEdXVsMomLeptonCutPi0El")->Fill(mom, dEdX); } else if (HDEAnaMyUtils::isMcProton(cand)) { fHM->H1("fhQPLeptonCutProton")->Fill(pq); fHM->H2("fhBetaVsMomLeptonCutProton")->Fill(mom, beta); fHM->H2("fhdEdXVsMomLeptonCutProton")->Fill(mom, dEdX); } else if (HDEAnaMyUtils::isMcPion(cand)) { fHM->H1("fhQPLeptonCutPion")->Fill(pq); fHM->H2("fhBetaVsMomLeptonCutPion")->Fill(mom, beta); fHM->H2("fhdEdXVsMomLeptonCutPion")->Fill(mom, dEdX); } else { fHM->H1("fhQPLeptonCutOther")->Fill(pq); fHM->H2("fhBetaVsMomLeptonCutOther")->Fill(mom, beta); fHM->H2("fhdEdXVsMomLeptonCutOther")->Fill(mom, dEdX); } } // Fill histogramms after Lepton cut and Momentum cut if (isLepton && isMomCutOk) { if (HDEAnaMyUtils::isMcSignalElectron(cand, fCatKine)) { fHM->H1("fhQPMomCutSignalEl")->Fill(pq, fSignalMBR); fHM->H2("fhBetaVsMomMomCutSignalEl")->Fill(mom, beta, fSignalMBR); fHM->H2("fhdEdXVsMomMomCutSignalEl")->Fill(mom, dEdX, fSignalMBR); } else if (HDEAnaMyUtils::isMcGammaElectron(cand, fCatKine)) { fHM->H1("fhQPMomCutGammaEl")->Fill(pq); fHM->H2("fhBetaVsMomMomCutGammaEl")->Fill(mom, beta); fHM->H2("fhdEdXVsMomMomCutGammaEl")->Fill(mom, dEdX); } else if (HDEAnaMyUtils::isMcPi0Electron(cand, fCatKine)) { fHM->H1("fhQPMomCutPi0El")->Fill(pq); fHM->H2("fhBetaVsMomMomCutPi0El")->Fill(mom, beta); fHM->H2("fhdEdXVsMomMomCutPi0El")->Fill(mom, dEdX); } else if (HDEAnaMyUtils::isMcProton(cand)) { fHM->H1("fhQPMomCutProton")->Fill(pq); fHM->H2("fhBetaVsMomMomCutProton")->Fill(mom, beta); fHM->H2("fhdEdXVsMomMomCutProton")->Fill(mom, dEdX); } else if (HDEAnaMyUtils::isMcPion(cand)) { fHM->H1("fhQPMomCutPion")->Fill(pq); fHM->H2("fhBetaVsMomMomCutPion")->Fill(mom, beta); fHM->H2("fhdEdXVsMomMomCutPion")->Fill(mom, dEdX); } else { fHM->H1("fhQPMomCutOther")->Fill(pq); fHM->H2("fhBetaVsMomMomCutOther")->Fill(mom, beta); fHM->H2("fhdEdXVsMomMomCutOther")->Fill(mom, dEdX); } } // Fill histogramms after Lepton cut and Momentum cut and OpenAngle cut if (isLepton && isMomCutOk && isOACutOk) { if (HDEAnaMyUtils::isMcSignalElectron(cand, fCatKine)) { fHM->H1("fhQPOACutSignalEl")->Fill(pq, fSignalMBR); fHM->H2("fhBetaVsMomOACutSignalEl")->Fill(mom, beta, fSignalMBR); fHM->H2("fhdEdXVsMomOACutSignalEl")->Fill(mom, dEdX, fSignalMBR); } else if (HDEAnaMyUtils::isMcGammaElectron(cand, fCatKine)) { fHM->H1("fhQPOACutGammaEl")->Fill(pq); fHM->H2("fhBetaVsMomOACutGammaEl")->Fill(mom, beta); fHM->H2("fhdEdXVsMomOACutGammaEl")->Fill(mom, dEdX); } else if (HDEAnaMyUtils::isMcPi0Electron(cand, fCatKine)) { fHM->H1("fhQPOACutPi0El")->Fill(pq); fHM->H2("fhBetaVsMomOACutPi0El")->Fill(mom, beta); fHM->H2("fhdEdXVsMomOACutPi0El")->Fill(mom, dEdX); } else if (HDEAnaMyUtils::isMcProton(cand)) { fHM->H1("fhQPOACutProton")->Fill(pq); fHM->H2("fhBetaVsMomOACutProton")->Fill(mom, beta); fHM->H2("fhdEdXVsMomOACutProton")->Fill(mom, dEdX); } else if (HDEAnaMyUtils::isMcPion(cand)) { fHM->H1("fhQPOACutPion")->Fill(pq); fHM->H2("fhBetaVsMomOACutPion")->Fill(mom, beta); fHM->H2("fhdEdXVsMomOACutPion")->Fill(mom, dEdX); } else { fHM->H1("fhQPOACutOther")->Fill(pq); fHM->H2("fhBetaVsMomOACutOther")->Fill(mom, beta); fHM->H2("fhdEdXVsMomOACutOther")->Fill(mom, dEdX); } } } } void HDEAnaTask::fillRichHists() { Int_t nofCands = fCatParticleCand->getEntries(); for(Int_t i = 0; i< nofCands; i++){ HParticleCandSim* cand = static_cast(fCatParticleCand->getObject(i)); Bool_t isLepton = fCuts.isLeptonOk(cand); Int_t richIndex = cand->getRichInd(); if (richIndex >= 0) { HRichHitSim* ring = static_cast(fCatRichHit->getObject(richIndex)); if (ring != NULL) { Double_t mom = cand->getMomentum(); Double_t xC = ring->getRich700CircleCenterX(); Double_t yC = ring->getRich700CircleCenterY(); Double_t radius = ring->getRich700CircleRadius(); Double_t nofHits = ring->getRich700NofRichCals(); Float_t theta = ring->getTheta(); if (HDEAnaMyUtils::isMcSignalElectron(cand, fCatKine)) { fHM->H2("fhNofHitsInRingVsThetaSignalElCands")->Fill( theta, nofHits); } else if (HDEAnaMyUtils::isMcGammaElectron(cand, fCatKine)) { fHM->H2("fhNofHitsInRingVsThetaGammaElCands")->Fill(theta, nofHits); } else if (HDEAnaMyUtils::isMcPi0Electron(cand, fCatKine)) { fHM->H2("fhNofHitsInRingVsThetaPi0ElCands")->Fill(theta, nofHits); } else { fHM->H2("fhNofHitsInRingVsThetaOtherCands")->Fill(theta, nofHits); } if(isLepton) { if (HDEAnaMyUtils::isMcSignalElectron(cand, fCatKine)) { fHM->H2("fhNofHitsInRingVsThetaSignalElCandsLeptonCut")->Fill( theta, nofHits); } else if (HDEAnaMyUtils::isMcGammaElectron(cand, fCatKine)) { fHM->H2("fhNofHitsInRingVsThetaGammaElCandsLeptonCut")->Fill(theta, nofHits); } else if (HDEAnaMyUtils::isMcPi0Electron(cand, fCatKine)) { fHM->H2("fhNofHitsInRingVsThetaPi0ElCandsLeptonCut")->Fill(theta, nofHits); } else { fHM->H2("fhNofHitsInRingVsThetaOtherCandsLeptonCut")->Fill(theta, nofHits); } } fHM->H2("fhRingXYAll")->Fill(xC, yC); if (isLepton) { fHM->H2("fhRingXYLeptonCut")->Fill(xC, yC); } fHM->H1("fhRadiusVsNofHitsAll")->Fill(nofHits, radius); fHM->H1("fhRadiusVsMomAll")->Fill(mom, radius); if (HDEAnaMyUtils::isMcSignalElectron(cand, fCatKine)) { fHM->H2("fhRingXYSignalEl")->Fill(xC, yC); fHM->H1("fhRadiusVsNofHitsSignalEl")->Fill(nofHits, radius); fHM->H1("fhRadiusVsMomSignalEl")->Fill(mom, radius); } else if (HDEAnaMyUtils::isMcGammaElectron(cand, fCatKine)) { fHM->H2("fhRingXYGammaEl")->Fill(xC, yC); fHM->H1("fhRadiusVsNofHitsGammaEl")->Fill(nofHits, radius); fHM->H1("fhRadiusVsMomGammaEl")->Fill(mom, radius); } else if (HDEAnaMyUtils::isMcPi0Electron(cand, fCatKine)) { fHM->H2("fhRingXYPi0El")->Fill(xC, yC); fHM->H1("fhRadiusVsNofHitsPi0El")->Fill(nofHits, radius); fHM->H1("fhRadiusVsMomPi0El")->Fill(mom, radius); } } } } } void HDEAnaTask::fillOpenAngleToBackgroundHists() { Int_t nofCands = fCatParticleCand->getEntries(); // loop plus candidates for(Int_t iP = 0; iP < nofCands; iP++){ HParticleCandSim* candP = static_cast(fCatParticleCand->getObject(iP)); if (candP->getCharge() < 0) continue; candP->calc4vectorProperties(); // loop minus candidates for(Int_t iM = 0; iM < nofCands; iM++){ if (iP == iM) continue; HParticleCandSim* candM = static_cast(fCatParticleCand->getObject(iM)); if (candM->getCharge() > 0) continue; candM->calc4vectorProperties(); Float_t openingAngle = HParticleTool::getOpeningAngle(candP, candM); Bool_t isLepton = (fCuts.isLeptonOk(candP) && fCuts.isLeptonOk(candM)); Bool_t isMomCut = (fCuts.isMomCutOk(candP) && fCuts.isMomCutOk(candM)); Bool_t isOACut = (fCuts.isOACutOk(iP) && fCuts.isOACutOk(iM)); if (HDEAnaMyUtils::isMcSignalElectron(candP, fCatKine) && !(HDEAnaMyUtils::isMcSignalElectron(candM, fCatKine))) { fHM->H1("fhOpenAngleSignalPositronToBG")->Fill(openingAngle); if(isLepton) fHM->H1("fhOpenAngleSignalPositronToBGLeptonCut")->Fill(openingAngle); if(isLepton && isMomCut) fHM->H1("fhOpenAngleSignalPositronToBGMomCut")->Fill(openingAngle); if(isLepton && isMomCut && isOACut) fHM->H1("fhOpenAngleSignalPositronToBGOACut")->Fill(openingAngle); } if (HDEAnaMyUtils::isMcSignalElectron(candM, fCatKine) && !(HDEAnaMyUtils::isMcSignalElectron(candP, fCatKine))) { fHM->H1("fhOpenAngleSignalElectronToBG")->Fill(openingAngle); if(isLepton) fHM->H1("fhOpenAngleSignalElectronToBGLeptonCut")->Fill(openingAngle); if(isLepton && isMomCut) fHM->H1("fhOpenAngleSignalElectronToBGMomCut")->Fill(openingAngle); if(isLepton && isMomCut && isOACut) fHM->H1("fhOpenAngleSignalElectronToBGOACut")->Fill(openingAngle); } } } } void HDEAnaTask::fillParticleCandHists() { Int_t nofCands = fCatParticleCand->getEntries(); Double_t minDist = 9999999.; Double_t minAngle = 9999999.; // loop plus candidates for(Int_t iP = 0; iP < nofCands; iP++){ HParticleCandSim* candP = static_cast(fCatParticleCand->getObject(iP)); if (candP->getCharge() < 0) continue; if (!fCuts.isLeptonOk(candP)) continue; candP->calc4vectorProperties(); Int_t richIndexP = candP->getRichInd(); HRichHitSim* ringP = (richIndexP >= 0)?static_cast(fCatRichHit->getObject(richIndexP)):NULL; // loop minus candidates for(Int_t iM = 0; iM < nofCands; iM++){ if (iP == iM) continue; HParticleCandSim* candM = static_cast(fCatParticleCand->getObject(iM)); if (candM->getCharge() > 0) continue; if (!fCuts.isLeptonOk(candM)) continue; candM->calc4vectorProperties(); Double_t momM = candM->getMomentum(); Float_t openingAngle = HParticleTool::getOpeningAngle(candP, candM); TLorentzVector pair; pair = *candM + *candP; Double_t minv = pair.M(); Double_t momPair = pair.P(); Bool_t isMomCut = (fCuts.isMomCutOk(candP) && fCuts.isMomCutOk(candM)); Bool_t isOACut = (fCuts.isOACutOk(iP) && fCuts.isOACutOk(iM)); Bool_t isSignalPair = HDEAnaMyUtils::isMcSignalPair(candP, candM, fCatKine); Bool_t isGammaPair = HDEAnaMyUtils::isMcGammaPair(candP, candM, fCatKine); Bool_t isPi0Pair = HDEAnaMyUtils::isMcPi0Pair(candP, candM, fCatKine); Bool_t isBgPair = !isGammaPair && !isPi0Pair && !(HDEAnaMyUtils::isMcSignalElectron(candM, fCatKine) || HDEAnaMyUtils::isMcSignalElectron(candP, fCatKine)); Int_t richIndexM = candM->getRichInd(); Double_t ringDist = -1.; if (richIndexM >= 0) { HRichHitSim* ringM = static_cast(fCatRichHit->getObject(richIndexM)); if (ringM != NULL && ringP != NULL) { Double_t xCP = ringP->getRich700CircleCenterX(); Double_t yCP = ringP->getRich700CircleCenterY(); Double_t xCM = ringM->getRich700CircleCenterX(); Double_t yCM = ringM->getRich700CircleCenterY(); ringDist = sqrt((xCM-xCP)*(xCM-xCP)+(yCP-yCM)*(yCP-yCM)); } } if (ringDist >= 0. && minDist >= ringDist) { minDist = ringDist; if (minAngle > openingAngle) { minAngle = openingAngle; } } if (isSignalPair) { fHM->H1("fhOpenAngleSignal")->Fill(openingAngle, fSignalMBR); fHM->H1("fhMinvSignal")->Fill(minv, fSignalMBR); if (isMomCut) fHM->H1("fhMinvMomCutSignal")->Fill(minv, fSignalMBR); if (isMomCut && isOACut) fHM->H1("fhMinvOACutSignal")->Fill(minv, fSignalMBR); if (ringDist >= 0.){ fHM->H2("fhRingDistVsMomPairSignal")->Fill(momPair, ringDist, fSignalMBR); fHM->H2("fhOpenAngleVsRingDistSignal")->Fill(openingAngle, ringDist, fSignalMBR); } } if (isBgPair) { fHM->H1("fhOpenAngleBg")->Fill(openingAngle); fHM->H1("fhMinvBg")->Fill(minv); if (isMomCut) fHM->H1("fhMinvMomCutBg")->Fill(minv); if (isMomCut && isOACut) fHM->H1("fhMinvOACutBg")->Fill(minv); if (ringDist >= 0.){ fHM->H2("fhRingDistVsMomPairBg")->Fill(momPair, ringDist); fHM->H2("fhOpenAngleVsRingDistBg")->Fill(openingAngle, ringDist); } //Minus particle on x-axis: gamma e-, pi0 e-, pion-, other //Plus particle on y-axis: gamma e+, pi0 e+, pion+, other bool isDifferentParentId = (candM->getGeantParentTrackNum() != candP->getGeantParentTrackNum()); // candM bool isCandMGammaEl = HDEAnaMyUtils::isMcGammaElectron(candM, fCatKine); bool isCandMPi0El = HDEAnaMyUtils::isMcPi0Electron(candM, fCatKine); bool isCandMPion = HDEAnaMyUtils::isMcPion(candM); // candP bool isCandPGammaEl = HDEAnaMyUtils::isMcGammaElectron(candP, fCatKine); bool isCandPPi0El = HDEAnaMyUtils::isMcPi0Electron(candP, fCatKine); bool isCandPPion = HDEAnaMyUtils::isMcPion(candP); if ( isCandMGammaEl ) { if (isCandPGammaEl && isDifferentParentId){ fillPairSources(0, 0, true, isMomCut, isOACut); } else if ( isCandPPi0El ){ fillPairSources(0, 1, true, isMomCut, isOACut); } else if ( isCandPPion ){ fillPairSources(0, 2, true, isMomCut, isOACut); }else{ fillPairSources(0, 3, true, isMomCut, isOACut); } }else if ( isCandMPi0El ) { if (isCandPGammaEl ){ fillPairSources(1, 0, true, isMomCut, isOACut); } else if ( isCandPPi0El && isDifferentParentId ){ fillPairSources(1, 1, true, isMomCut, isOACut); } else if ( isCandPPion ){ fillPairSources(1, 2, true, isMomCut, isOACut); }else{ fillPairSources(1, 3, true, isMomCut, isOACut); } } else if ( isCandMPion ) { if (isCandPGammaEl ){ fillPairSources(2, 0, true, isMomCut, isOACut); } else if ( isCandPPi0El ){ fillPairSources(2, 1, true, isMomCut, isOACut); } else if ( isCandPPion ){ fillPairSources(2, 2, true, isMomCut, isOACut); }else{ fillPairSources(2, 3, true, isMomCut, isOACut); } } else { if (isCandPGammaEl ){ fillPairSources(3, 0, true, isMomCut, isOACut); } else if ( isCandPPi0El ){ fillPairSources(3, 1, true, isMomCut, isOACut); } else if ( isCandPPion ){ fillPairSources(3, 2, true, isMomCut, isOACut); }else{ fillPairSources(3, 3, true, isMomCut, isOACut); } } } //isBgPair if (isGammaPair) { fHM->H1("fhOpenAngleGamma")->Fill(openingAngle); fHM->H1("fhMinvGamma")->Fill(minv); if (isMomCut) fHM->H1("fhMinvMomCutGamma")->Fill(minv); if (isMomCut && isOACut) fHM->H1("fhMinvOACutGamma")->Fill(minv); if (ringDist >= 0.) { fHM->H2("fhRingDistVsMomPairGamma")->Fill(momPair, ringDist); fHM->H2("fhOpenAngleVsRingDistGamma")->Fill(openingAngle, ringDist); } } if (isPi0Pair) { fHM->H1("fhOpenAnglePi0")->Fill(openingAngle); fHM->H1("fhMinvPi0")->Fill(minv); if (isMomCut) fHM->H1("fhMinvMomCutPi0")->Fill(minv); if (isMomCut && isOACut) fHM->H1("fhMinvOACutPi0")->Fill(minv); if (ringDist >= 0.) { fHM->H2("fhRingDistVsMomPairPi0")->Fill(momPair, ringDist); fHM->H2("fhOpenAngleVsRingDistPi0")->Fill(openingAngle, ringDist); } } } } fHM->H2("fhMinDistVsMinAngle")->Fill(minDist, minAngle); } void HDEAnaTask::fillPairSources(double xBin, double yBin, bool isLeptonCut, bool isMomCut, bool isOACut) { if (isLeptonCut) fHM->H2("fhBgPairSourcesLeptonCut")->Fill(xBin, yBin); if (isLeptonCut && isMomCut) fHM->H2("fhBgPairSourcesMomCut")->Fill(xBin, yBin); if (isLeptonCut && isMomCut && isOACut) fHM->H2("fhBgPairSourcesOACut")->Fill(xBin, yBin); } void HDEAnaTask::fillKineHists() { Int_t nofKine = fCatKine->getEntries(); Int_t nofPi0 = 0; Int_t nofPi0Primary = 0; for(Int_t i=0; i< nofKine; i++){ HGeantKine* kine = (HGeantKine*)fCatKine->getObject(i); Double_t mom = kine->getTotalMomentum(); Float_t transmom = kine->getTransverseMomentum(); Int_t kineID = kine->getID(); Int_t mechanism = kine->getMechanism(); fHM->H1("fhKineID")->Fill(kineID); fHM->H1("fhKineMechanism")->Fill(mechanism); if (HDEAnaMyUtils::isMcSignalElectron(kine)) { fHM->H1("fhMomKineSignalEl")->Fill(mom, fSignalMBR); fHM->H1("fhTransMomKineSignalEl")->Fill(transmom, fSignalMBR); } if (HDEAnaMyUtils::isMcGammaElectron(kine, fCatKine)) { fHM->H1("fhMomKineGammaEl")->Fill(mom); fHM->H1("fhTransMomKineGammaEl")->Fill(transmom); } if (HDEAnaMyUtils::isMcPi0Electron(kine, fCatKine)) { fHM->H1("fhMomKinePi0El")->Fill(mom); fHM->H1("fhTransMomKinePi0El")->Fill(transmom); } if(HDEAnaMyUtils::isMcPi0(kine)) { nofPi0++; } if(HDEAnaMyUtils::isMcPi0Primary(kine)) { nofPi0Primary++; } if(kine->getID()==1) { if(kine->isPrimary()) { fHM->H1("fhParentFromGamma")->Fill(0); } Int_t parentInd = kine->getParentTrack() - 1; if (parentInd < 0) continue; HGeantKine* kineParent = static_cast(fCatKine->getObject(parentInd)); if (kineParent != NULL) { Int_t kineParentID = kineParent->getID(); fHM->H1("fhParentFromGamma")->Fill(kineParentID); } } } fHM->H1("fhNofKinePi0")->Fill(nofPi0); fHM->H1("fhNofKinePi0Primary")->Fill(nofPi0Primary); for(Int_t iP=0; iP < nofKine; iP++){ HGeantKine* kineP = static_cast(fCatKine->getObject(iP)); if (!kineP->getID()== 2 ) continue; for(Int_t iM=0; iM < nofKine; iM++){ HGeantKine* kineM = static_cast(fCatKine->getObject(iM)); if (!kineM->getID() == 3 ) continue; if (HDEAnaMyUtils::isMcSignalElectron(kineP) && HDEAnaMyUtils::isMcSignalElectron(kineM)) { Float_t openingAngle = HParticleTool::getOpeningAngle(kineP, kineM); fHM->H1("fhOpenAngleKineSignalEl")->Fill(openingAngle); } else if (HDEAnaMyUtils::isMcGammaElectron(kineP, fCatKine) && HDEAnaMyUtils::isMcGammaElectron(kineM, fCatKine) && kineP->getParentTrack() == kineM->getParentTrack()) { Float_t openingAngle = HParticleTool::getOpeningAngle(kineP, kineM); fHM->H1("fhOpenAngleKineGammaEl")->Fill(openingAngle); } else if (HDEAnaMyUtils::isMcPi0Electron(kineP, fCatKine) && HDEAnaMyUtils::isMcPi0Electron(kineM, fCatKine) && kineP->getParentTrack() == kineM->getParentTrack()) { Float_t openingAngle = HParticleTool::getOpeningAngle(kineP, kineM); fHM->H1("fhOpenAngleKinePi0El")->Fill(openingAngle); } } } } void HDEAnaTask::drawHist() { SetDefaultDrawStyle(); cout << "Number of analysed events:" << fHM->H1("fhNofEvents")->GetEntries() << endl; { TCanvas* c = fHM->CreateCanvas("dea_nof_cands", "dea_nof_cands", 1200, 600); c->Divide(2,1); c->cd(1); vector h1 = fHM->GetH1ByNames("fhNofCandsAll", "fhNofCandsLepton", "fhNofCandsGood", "fhNofCandsGammaEl", "fhNofCandsPi0El"); const string l[] = {"All", "Lepton ID", "Lepton Ring", "#gamma^{0} e^{#pm}", "#pi^{0} e^{#pm}"}; vector lM = HDEAnaMyUtils::getLegendLabelsWithMean(ToVector(l), h1); DrawH1(h1, lM, kLinear, kLog, true, 0.7, 0.80, 0.99, 0.99); DrawTextOnPad("All range", 0.2, 0.9, 0.75, 0.99); c->cd(2); Double_t zoomMin = -0.5; Double_t zoomMax = 9.5; vector h2; h2.push_back( (TH1D*)fHM->H1("fhNofCandsAll")->Clone() ); h2.push_back( (TH1D*)fHM->H1("fhNofCandsLepton")->Clone() ); h2.push_back( (TH1D*)fHM->H1("fhNofCandsGood")->Clone() ); h2.push_back( (TH1D*)fHM->H1("fhNofCandsGammaEl")->Clone() ); h2.push_back( (TH1D*)fHM->H1("fhNofCandsPi0El")->Clone() ); for (UInt_t i = 0; i < h2.size(); i++) { h2[i]->GetXaxis()->SetRangeUser(zoomMin, zoomMax); } DrawH1(h2, lM, kLinear, kLog, true, 0.7, 0.80, 0.99, 0.99); DrawTextOnPad("Zoom", 0.2, 0.9, 0.75, 0.99); } { TCanvas* c = fHM->CreateCanvas("dea_nof_rings","dea_nof_rings", 1000, 800); vector h; h.push_back(fHM->H1("fhNofRingsAll")); vector l; l.push_back("All ("+ HDEAna::NumberToString(fHM->H1("fhNofRingsAll")->GetMean()) + ")"); DrawH1(h,l, kLinear, kLinear, true, 0.8, 0.85, 0.99, 0.99); } { TCanvas* c = fHM->CreateCanvas("dea_nof_shared_particle_cands_for_ring", "dea_nof_shared_particle_cands_for_ring", 1000, 800); vector h; h.push_back(fHM->H1("fhNofSharedCandidatesForRing")); vector l; l.push_back("All ("+ HDEAna::NumberToString(fHM->H1("fhNofSharedCandidatesForRing")->GetMean()) + ")"); DrawH1(h,l, kLinear, kLinear, true, 0.8, 0.85, 0.99, 0.99); } { fHM->CreateCanvas("dea_nof_hits_in_ring", "dea_nof_hits_in_ring", 1000, 800); vector h = fHM->GetH2ProjectionYByNames("fhNofHitsInRingVsThetaSignalEl", "fhNofHitsInRingVsThetaGammaEl", "fhNofHitsInRingVsThetaPi0El", "fhNofHitsInRingVsThetaOther"); fHM->NormalizeToIntegral(h); const string l[] = {"Signal e^{#pm}", "#gamma e^{#pm}", "#pi^{0} e^{#pm}", "Other"}; vector lM = HDEAnaMyUtils::getLegendLabelsWithMean(ToVector(l), h); DrawH1(h, lM, kLinear, kLinear, true, 0.6, 0.8, 0.99, 0.99); } { fHM->CreateCanvas("dea_nof_hits_in_ring_cands", "dea_nof_hits_in_ring_cands", 1000, 800); vector h = fHM->GetH2ProjectionYByNames("fhNofHitsInRingVsThetaSignalElCands", "fhNofHitsInRingVsThetaGammaElCands", "fhNofHitsInRingVsThetaPi0ElCands", "fhNofHitsInRingVsThetaOtherCands"); fHM->NormalizeToIntegral(h); const string l[] = {"Signal e^{#pm}", "#gamma e^{#pm}", "#pi^{0} e^{#pm}", "Other"}; vector lM = HDEAnaMyUtils::getLegendLabelsWithMean(ToVector(l), h); DrawH1(h, lM, kLinear, kLinear, true, 0.6, 0.8, 0.99, 0.99); } { fHM->CreateCanvas("dea_nof_hits_in_ring_cands_lepton_cut", "dea_nof_hits_in_ring_cands_lepton_cut", 1000, 800); vector h = fHM->GetH2ProjectionYByNames("fhNofHitsInRingVsThetaSignalElCandsLeptonCut", "fhNofHitsInRingVsThetaGammaElCandsLeptonCut", "fhNofHitsInRingVsThetaPi0ElCandsLeptonCut", "fhNofHitsInRingVsThetaOtherCandsLeptonCut"); fHM->NormalizeToIntegral(h); const string l[] = {"Signal e^{#pm}", "#gamma e^{#pm}", "#pi^{0} e^{#pm}", "Other"}; vector lM = HDEAnaMyUtils::getLegendLabelsWithMean(ToVector(l), h); DrawH1(h, lM, kLinear, kLinear, true, 0.6, 0.8, 0.99, 0.99); } { TCanvas* c = fHM->CreateCanvas("dea_nof_hits_in_ring_vs_theta", "dea_nof_hits_in_ring_vs_theta", 1000, 800); c->Divide(2, 2); c->cd(1); DrawH2WithProfile(fHM->H2("fhNofHitsInRingVsThetaSignalEl")); DrawTextOnPad("Signal e^{#pm}", 0.6, 0.9, 0.75, 0.99); c->cd(2); DrawH2WithProfile(fHM->H2("fhNofHitsInRingVsThetaGammaEl")); DrawTextOnPad("#gamma e^{#pm}", 0.6, 0.9, 0.75, 0.99); c->cd(3); DrawH2WithProfile(fHM->H2("fhNofHitsInRingVsThetaPi0El")); DrawTextOnPad("#pi^{0} e^{#pm}", 0.6, 0.9, 0.75, 0.99); c->cd(4); DrawH2WithProfile(fHM->H2("fhNofHitsInRingVsThetaOther")); DrawTextOnPad("Other e^{#pm}", 0.6, 0.9, 0.75, 0.99); } { TCanvas* c = fHM->CreateCanvas("dea_nof_hits_in_ring_vs_theta_shared", "dea_nof_hits_in_ring_vs_theta_shared", 1000, 800); c->Divide(2, 2); c->cd(1); DrawH2WithProfile(fHM->H2("fhNofHitsInRingVsThetaSignalElShared")); DrawTextOnPad("Signal e^{#pm}", 0.6, 0.9, 0.75, 0.99); c->cd(2); DrawH2WithProfile(fHM->H2("fhNofHitsInRingVsThetaGammaElShared")); DrawTextOnPad("#gamma e^{#pm}", 0.6, 0.9, 0.75, 0.99); c->cd(3); DrawH2WithProfile(fHM->H2("fhNofHitsInRingVsThetaPi0ElShared")); DrawTextOnPad("#pi^{0} e^{#pm}", 0.6, 0.9, 0.75, 0.99); c->cd(4); DrawH2WithProfile(fHM->H2("fhNofHitsInRingVsThetaOtherShared")); DrawTextOnPad("Other e^{#pm}", 0.6, 0.9, 0.75, 0.99); } { TCanvas* c = fHM->CreateCanvas("dea_nof_hits_in_ring_vs_theta_notshared", "dea_nof_hits_in_ring_vs_theta_notshared", 1000, 800); c->Divide(2, 2); c->cd(1); DrawH2WithProfile(fHM->H2("fhNofHitsInRingVsThetaSignalElNotShared")); DrawTextOnPad("Signal e^{#pm}", 0.6, 0.9, 0.75, 0.99); c->cd(2); DrawH2WithProfile(fHM->H2("fhNofHitsInRingVsThetaGammaElNotShared")); DrawTextOnPad("#gamma e^{#pm}", 0.6, 0.9, 0.75, 0.99); c->cd(3); DrawH2WithProfile(fHM->H2("fhNofHitsInRingVsThetaPi0ElNotShared")); DrawTextOnPad("#pi^{0} e^{#pm}", 0.6, 0.9, 0.75, 0.99); c->cd(4); DrawH2WithProfile(fHM->H2("fhNofHitsInRingVsThetaOtherNotShared")); DrawTextOnPad("Other e^{#pm}", 0.6, 0.9, 0.75, 0.99); } { TCanvas* c = fHM->CreateCanvas("dea_nof_hits_in_ring_shared", "dea_nof_hits_in_ring_shared", 1000, 800); vector h = fHM->GetH1ByNames("fhNofHitsInRingNotShared", "fhNofHitsInRingNotSharedSignalEl", "fhNofHitsInRingNotAssigned", "fhNofHitsInRingShared"); fHM->NormalizeToIntegral(h); const string l[] = {"Not shared", "Not shared signal e^{#pm}", "Not assigned", "Shared"}; vector lM = HDEAnaMyUtils::getLegendLabelsWithMean(ToVector(l), h); DrawH1(h, lM, kLinear, kLinear, true, 0.6, 0.8, 0.99, 0.99); } { TCanvas* c = fHM->CreateCanvas("dea_geant_pid", "dea_geant_pid", 1500, 500); c->Divide(3,1); c->cd(1); vector h1 = fHM->GetH1ByNames("fhGeantPidAll", "fhGeantPidLeptonCutAll", "fhGeantPidMomCutAll", "fhGeantPidOACutAll"); const string l[] = {"No cuts", "+Lepton cut", "+Momentum cut", "+OAngle cut"}; DrawH1(h1, ToVector(l), kLinear, kLog, true, 0.75, 0.8, 0.99, 0.99); DrawTextOnPad("All", 0.2, 0.9, 0.75, 0.99); c->cd(2); vector h2 = fHM->GetH1ByNames("fhGeantPidSignal", "fhGeantPidLeptonCutSignal", "fhGeantPidMomCutSignal", "fhGeantPidOACutSignal"); DrawH1(h2, ToVector(l), kLinear, kLog, true, 0.75, 0.8, 0.99, 0.99); DrawTextOnPad("Signal", 0.2, 0.9, 0.75, 0.99); h2[0]->SetMinimum(1); c->cd(3); vector h3 = fHM->GetH1ByNames("fhGeantPidBg", "fhGeantPidLeptonCutBg", "fhGeantPidMomCutBg", "fhGeantPidOACutBg"); DrawH1(h3, ToVector(l), kLinear, kLog, true, 0.75, 0.8, 0.99, 0.99); DrawTextOnPad("BG", 0.2, 0.9, 0.75, 0.99); } { TCanvas* c = fHM->CreateCanvas("dea_qp", "dea_qp", 1500, 1000); c->Divide(3, 2); c->cd(1); vector h = fHM->GetH1ByNames("fhQPGammaEl", "fhQPPi0El", "fhQPPion", "fhQPProton", "fhQPOther"); const string l[] = {"#gamma e^{#pm}", "#pi^{0} e^{#pm}", "#pi^{#pm}", "protons", "Other"}; DrawH1(h, ToVector(l), kLinear, kLog, true, 0.85, 0.75, 0.99, 0.99); DrawTextOnPad("BG no cuts", 0.2, 0.9, 0.75, 0.99); c->cd(2); vector h2 = fHM->GetH1ByNames("fhQPLeptonCutGammaEl", "fhQPLeptonCutPi0El", "fhQPLeptonCutPion", "fhQPLeptonCutProton", "fhQPLeptonCutOther"); DrawH1(h2, ToVector(l), kLinear, kLog, true, 0.85, 0.75, 0.99, 0.99); DrawTextOnPad("BG Lepton cut", 0.2, 0.9, 0.75, 0.99); c->cd(3); vector h3 = fHM->GetH1ByNames("fhQPMomCutGammaEl", "fhQPMomCutPi0El", "fhQPMomCutPion", "fhQPMomCutProton", "fhQPMomCutOther"); DrawH1(h3, ToVector(l), kLinear, kLog, true, 0.85, 0.75, 0.99, 0.99); DrawTextOnPad("BG momentum cut", 0.2, 0.9, 0.75, 0.99); c->cd(4); vector h4 = fHM->GetH1ByNames("fhQPOACutGammaEl", "fhQPOACutPi0El", "fhQPOACutPion", "fhQPOACutProton", "fhQPOACutOther"); DrawH1(h4, ToVector(l), kLinear, kLog, true, 0.85, 0.75, 0.99, 0.99); DrawTextOnPad("BG OAngle cut", 0.2, 0.9, 0.75, 0.99); c->cd(5); vector h5 = fHM->GetH1ByNames("fhQPSignalEl", "fhQPLeptonCutSignalEl", "fhQPMomCutSignalEl", "fhQPOACutSignalEl"); const string l5[] = {"No cuts", "Lepton cut", "Mom cut", "OAngle cut"}; DrawH1(h5, ToVector(l5), kLinear, kLog, true, 0.75, 0.85, 0.99, 0.99); DrawTextOnPad("Signal e^{#pm}", 0.2, 0.9, 0.75, 0.99); } { TCanvas* c = fHM->CreateCanvas("dea_beta", "dea_beta", 1000, 1000); c->Divide(2, 2); c->cd(1); vector h = fHM->GetH2ProjectionYByNames("fhBetaVsMomSignalEl", "fhBetaVsMomGammaEl", "fhBetaVsMomPi0El", "fhBetaVsMomPion", "fhBetaVsMomProton", "fhBetaVsMomOther"); fHM->NormalizeToIntegral(h); const string l[] = {"Signal e^{#pm}", "#gamma e^{#pm}", "#pi^{0} e^{#pm}", "#pi^{#pm}", "protons", "Other"}; DrawH1(h, ToVector(l), kLinear, kLog, true, 0.85, 0.75, 0.99, 0.99); DrawTextOnPad("No cuts", 0.2, 0.9, 0.75, 0.99); c->cd(2); vector h2 = fHM->GetH2ProjectionYByNames("fhBetaVsMomLeptonCutSignalEl", "fhBetaVsMomLeptonCutGammaEl", "fhBetaVsMomLeptonCutPi0El", "fhBetaVsMomLeptonCutPion", "fhBetaVsMomLeptonCutProton", "fhBetaVsMomLeptonCutOther"); fHM->NormalizeToIntegral(h2); DrawH1(h2, ToVector(l), kLinear, kLog, true, 0.85, 0.75, 0.99, 0.99); DrawTextOnPad("Lepton cut", 0.2, 0.9, 0.75, 0.99); c->cd(3); vector h3 = fHM->GetH2ProjectionYByNames("fhBetaVsMomLeptonCutSignalEl", "fhBetaVsMomMomCutGammaEl", "fhBetaVsMomMomCutPi0El", "fhBetaVsMomMomCutPion", "fhBetaVsMomMomCutProton", "fhBetaVsMomMomCutOther"); fHM->NormalizeToIntegral(h3); DrawH1(h3, ToVector(l), kLinear, kLog, true, 0.85, 0.75, 0.99, 0.99); DrawTextOnPad("Momentum cut", 0.2, 0.9, 0.75, 0.99); c->cd(4); vector h4 = fHM->GetH2ProjectionYByNames("fhBetaVsMomLeptonCutSignalEl", "fhBetaVsMomOACutGammaEl", "fhBetaVsMomOACutPi0El", "fhBetaVsMomOACutPion", "fhBetaVsMomOACutProton", "fhBetaVsMomOACutOther"); fHM->NormalizeToIntegral(h4); DrawH1(h4, ToVector(l), kLinear, kLog, true, 0.85, 0.75, 0.99, 0.99); DrawTextOnPad("OAngle cut", 0.2, 0.9, 0.75, 0.99); } { TCanvas* c = fHM->CreateCanvas("dea_beta_vs_mom", "dea_beta_vs_mom", 1500, 1000); c->Divide(3, 2); c->cd(1); DrawH2(fHM->H2("fhBetaVsMomSignalEl")); DrawTextOnPad("Signal e^{#pm}", 0.2, 0.9, 0.75, 0.99); c->cd(2); DrawH2(fHM->H2("fhBetaVsMomGammaEl")); DrawTextOnPad("#gamma e^{#pm}", 0.2, 0.9, 0.75, 0.99); c->cd(3); DrawH2(fHM->H2("fhBetaVsMomPi0El")); DrawTextOnPad("#pi^{0} e^{#pm}", 0.2, 0.9, 0.75, 0.99); c->cd(4); DrawH2(fHM->H2("fhBetaVsMomPion")); DrawTextOnPad("#pi^{#pm}", 0.2, 0.9, 0.75, 0.99); c->cd(5); DrawH2(fHM->H2("fhBetaVsMomProton")); DrawTextOnPad("protons", 0.2, 0.9, 0.75, 0.99); c->cd(6); DrawH2(fHM->H2("fhBetaVsMomOther")); DrawTextOnPad("Other", 0.2, 0.9, 0.75, 0.99); } { TCanvas* c = fHM->CreateCanvas("dea_dedx", "dea_dedx", 1000, 1000); c->Divide(2, 2); c->cd(1); vector h = fHM->GetH2ProjectionYByNames("fhdEdXVsMomSignalEl", "fhdEdXVsMomGammaEl", "fhdEdXVsMomPi0El", "fhdEdXVsMomPion", "fhdEdXVsMomProton", "fhdEdXVsMomOther"); fHM->NormalizeToIntegral(h); const string l[] = {"Signal e^{#pm}", "#gamma e^{#pm}", "#pi^{0} e^{#pm}", "#pi^{#pm}", "protons", "Other"}; DrawH1(h, ToVector(l), kLinear, kLog, true, 0.85, 0.75, 0.99, 0.99); DrawTextOnPad("No cuts", 0.2, 0.9, 0.75, 0.99); c->cd(2); vector h2 = fHM->GetH2ProjectionYByNames("fhdEdXVsMomLeptonCutSignalEl", "fhdEdXVsMomLeptonCutGammaEl", "fhdEdXVsMomLeptonCutPi0El", "fhdEdXVsMomLeptonCutPion", "fhdEdXVsMomLeptonCutProton", "fhdEdXVsMomLeptonCutOther"); fHM->NormalizeToIntegral(h2); DrawH1(h2, ToVector(l), kLinear, kLog, true, 0.85, 0.75, 0.99, 0.99); DrawTextOnPad("Lepton cut", 0.2, 0.9, 0.75, 0.99); c->cd(3); vector h3 = fHM->GetH2ProjectionYByNames("fhdEdXVsMomLeptonCutSignalEl", "fhdEdXVsMomMomCutGammaEl", "fhdEdXVsMomMomCutPi0El", "fhdEdXVsMomMomCutPion", "fhdEdXVsMomMomCutProton", "fhdEdXVsMomMomCutOther"); fHM->NormalizeToIntegral(h3); DrawH1(h3, ToVector(l), kLinear, kLog, true, 0.85, 0.75, 0.99, 0.99); DrawTextOnPad("Momentum cut", 0.2, 0.9, 0.75, 0.99); c->cd(4); vector h4 = fHM->GetH2ProjectionYByNames("fhdEdXVsMomLeptonCutSignalEl", "fhdEdXVsMomOACutGammaEl", "fhdEdXVsMomOACutPi0El", "fhdEdXVsMomOACutPion", "fhdEdXVsMomOACutProton", "fhdEdXVsMomOACutOther"); fHM->NormalizeToIntegral(h4); DrawH1(h4, ToVector(l), kLinear, kLog, true, 0.85, 0.75, 0.99, 0.99); DrawTextOnPad("OAngle cut", 0.2, 0.9, 0.75, 0.99); } { TCanvas* c = fHM->CreateCanvas("dea_dedx_vs_mom", "dea_dedx_vs_mom", 1500, 1000); c->Divide(3, 2); c->cd(1); DrawH2(fHM->H2("fhdEdXVsMomSignalEl")); DrawTextOnPad("Signal e^{#pm}", 0.2, 0.9, 0.75, 0.99); c->cd(2); DrawH2(fHM->H2("fhdEdXVsMomGammaEl")); DrawTextOnPad("#gamma e^{#pm}", 0.2, 0.9, 0.75, 0.99); c->cd(3); DrawH2(fHM->H2("fhdEdXVsMomPi0El")); DrawTextOnPad("#pi^{0} e^{#pm}", 0.2, 0.9, 0.75, 0.99); c->cd(4); DrawH2(fHM->H2("fhdEdXVsMomPion")); DrawTextOnPad("#pi^{#pm}", 0.2, 0.9, 0.75, 0.99); c->cd(5); DrawH2(fHM->H2("fhdEdXVsMomProton")); DrawTextOnPad("protons", 0.2, 0.9, 0.75, 0.99); c->cd(6); DrawH2(fHM->H2("fhdEdXVsMomOther")); DrawTextOnPad("Other", 0.2, 0.9, 0.75, 0.99); } { TCanvas* c = fHM->CreateCanvas("dea_ring_xy", "dea_ring_xy", 1600, 1000); c->Divide(3, 2); c->cd(1); DrawH2(fHM->H2("fhRingXYAll")); DrawTextOnPad("All e^{#pm}", 0.6, 0.9, 0.75, 0.99); c->cd(2); DrawH2(fHM->H2("fhRingXYSignalEl")); DrawTextOnPad("Signal e^{#pm}", 0.6, 0.9, 0.75, 0.99); c->cd(3); DrawH2(fHM->H2("fhRingXYGammaEl")); DrawTextOnPad("#gamma e^{#pm}", 0.6, 0.9, 0.75, 0.99); c->cd(4); DrawH2(fHM->H2("fhRingXYPi0El")); DrawTextOnPad("#pi^{0} e^{#pm}", 0.6, 0.9, 0.75, 0.99); c->cd(5); DrawH2(fHM->H2("fhRingXYLeptonCut")); DrawTextOnPad("e^{#pm} Lepton Cut", 0.6, 0.9, 0.75, 0.99); } { fHM->CreateCanvas("dea_nof_kine_pi0","dea_nof_kine_pi0", 1000, 800); vector h = fHM->GetH1ByNames("fhNofKinePi0", "fhNofKinePi0Primary"); const string l[] = {"#pi^{0}_{all}", "#pi^{0}_{prim}"}; vector lM = HDEAnaMyUtils::getLegendLabelsWithMean(ToVector(l), h); DrawH1(h, lM, kLinear, kLinear, true, 0.75, 0.85, 0.99, 0.99); } { TCanvas* c = fHM->CreateCanvas("dea_kine_parent_from_gamma", "dea_kine_parent_from_gamma", 1000, 800); DrawH1(fHM->H1("fhParentFromGamma"), kLinear, kLog); } { TCanvas* c = fHM->CreateCanvas("dea_radius_vs_nof_hits", "dea_radius_vs_nof_hits", 1200, 1000); c->Divide(2, 2); c->cd(1); DrawH2WithProfile(fHM->H2("fhRadiusVsNofHitsAll")); DrawTextOnPad("All", 0.6, 0.9, 0.75, 0.99); c->cd(2); DrawH2WithProfile(fHM->H2("fhRadiusVsNofHitsSignalEl")); DrawTextOnPad("Signal e^{#pm}", 0.6, 0.9, 0.75, 0.99); c->cd(3); DrawH2WithProfile(fHM->H2("fhRadiusVsNofHitsGammaEl")); DrawTextOnPad("#gamma e^{#pm}", 0.6, 0.9, 0.75, 0.99); c->cd(4); DrawH2WithProfile(fHM->H2("fhRadiusVsNofHitsPi0El")); DrawTextOnPad("#pi^{0} e^{#pm}", 0.6, 0.9, 0.75, 0.99); } { TCanvas* c = fHM->CreateCanvas("dea_radius_vs_mom", "dea_radius_vs_mom", 1200, 1000); c->Divide(2, 2); c->cd(1); DrawH2WithProfile(fHM->H2("fhRadiusVsMomAll")); DrawTextOnPad("All", 0.6, 0.9, 0.75, 0.99); c->cd(2); DrawH2WithProfile(fHM->H2("fhRadiusVsMomSignalEl")); DrawTextOnPad("Signal e^{#pm}", 0.6, 0.9, 0.75, 0.99); c->cd(3); DrawH2WithProfile(fHM->H2("fhRadiusVsMomGammaEl")); DrawTextOnPad("#gamma e^{#pm}", 0.6, 0.9, 0.75, 0.99); c->cd(4); DrawH2WithProfile(fHM->H2("fhRadiusVsMomPi0El")); DrawTextOnPad("#pi^{0} e^{#pm}", 0.6, 0.9, 0.75, 0.99); } { TCanvas* c = fHM->CreateCanvas("dea_pair_reco_topology_2d", "dea_pair_reco_topology_2d", 1200, 600); c->Divide(2, 1); c->cd(1); fHM->H2("fhPairRecoTopology")->SetMarkerSize(2); DrawH2(fHM->H2("fhPairRecoTopology"), kLinear, kLinear, kLinear, "TEXT COLZ"); DrawTextOnPad("Lepton cut for 2nd candidate", 0.2, 0.9, 0.75, 0.99); c->cd(2); fHM->H2("fhPairRecoTopologyNoLeptonCut")->SetMarkerSize(2); DrawH2(fHM->H2("fhPairRecoTopologyNoLeptonCut"), kLinear, kLinear, kLinear, "TEXT COLZ"); DrawTextOnPad("No lepton cut for 2nd candidate", 0.2, 0.9, 0.75, 0.99); string labelsX[3] = {"All", "Pair", "No pair"}; for (Int_t i = 1; i <= 3; i++){ fHM->H2("fhPairRecoTopology")->GetXaxis()->SetBinLabel(i, labelsX[i-1].c_str()); fHM->H2("fhPairRecoTopologyNoLeptonCut")->GetXaxis()->SetBinLabel(i, labelsX[i-1].c_str()); } string labelsY[4] = { "#pi^{0} e^{-}", "#pi^{0} e^{+}", "#gamma e^{-}", "#gamma e^{+}"}; for (Int_t i = 1; i <= 4; i++){ fHM->H2("fhPairRecoTopology")->GetYaxis()->SetBinLabel(i, labelsY[i-1].c_str()); fHM->H2("fhPairRecoTopologyNoLeptonCut")->GetYaxis()->SetBinLabel(i, labelsY[i-1].c_str()); } } { TCanvas* c = fHM->CreateCanvas("dea_bg_pair_sources", "dea_bg_pair_sources", 1500, 500); c->Divide(3, 1); c->cd(1); fHM->H2("fhBgPairSourcesLeptonCut")->SetMarkerSize(2); DrawH2(fHM->H2("fhBgPairSourcesLeptonCut"), kLinear, kLinear, kLinear, "TEXT COLZ"); DrawTextOnPad("Lepton cut", 0.2, 0.9, 0.75, 0.99); c->cd(2); fHM->H2("fhBgPairSourcesMomCut")->SetMarkerSize(2); DrawH2(fHM->H2("fhBgPairSourcesMomCut"), kLinear, kLinear, kLinear, "TEXT COLZ"); DrawTextOnPad("Momentum cut", 0.2, 0.9, 0.75, 0.99); c->cd(3); fHM->H2("fhBgPairSourcesOACut")->SetMarkerSize(2); DrawH2(fHM->H2("fhBgPairSourcesOACut"), kLinear, kLinear, kLinear, "TEXT COLZ"); DrawTextOnPad("OAngle cut", 0.2, 0.9, 0.75, 0.99); //Minus particle on x-axis: gamma e-, pi0 e-, pion-, other //Plus particle on y-axis: gamma e+, pi0 e+, pion+, other string labelsX[4] = {"#gamma e^{-}", "#pi^{0} e^{-}", "#pi^{-}", "Other"}; for (Int_t i = 1; i <= 4; i++){ fHM->H2("fhBgPairSourcesLeptonCut")->GetXaxis()->SetBinLabel(i, labelsX[i-1].c_str()); fHM->H2("fhBgPairSourcesMomCut")->GetXaxis()->SetBinLabel(i, labelsX[i-1].c_str()); fHM->H2("fhBgPairSourcesOACut")->GetXaxis()->SetBinLabel(i, labelsX[i-1].c_str()); } string labelsY[4] = {"#gamma e^{+}", "#pi^{0} e^{+}", "#pi^{+}", "Other"}; for (Int_t i = 1; i <= 4; i++){ fHM->H2("fhBgPairSourcesLeptonCut")->GetYaxis()->SetBinLabel(i, labelsY[i-1].c_str()); fHM->H2("fhBgPairSourcesMomCut")->GetYaxis()->SetBinLabel(i, labelsY[i-1].c_str()); fHM->H2("fhBgPairSourcesOACut")->GetYaxis()->SetBinLabel(i, labelsY[i-1].c_str()); } } { fHM->CreateCanvas("dea_minv", "dea_minv", 1000, 800); vector h = fHM->GetH1ByNames("fhMinvSignal", "fhMinvGamma", "fhMinvPi0", "fhMinvBg"); const string l[] = {"Signal", "#gamma", "#pi^{0}", "BG"}; DrawH1(h, ToVector(l), kLinear, kLog); } { fHM->CreateCanvas("dea_minv_momcut", "dea_minv_momcut", 1000, 800); vector h = fHM->GetH1ByNames("fhMinvMomCutSignal", "fhMinvMomCutGamma", "fhMinvMomCutPi0", "fhMinvMomCutBg"); const string l[] = {"Signal", "#gamma", "#pi^{0}", "BG"}; DrawH1(h, ToVector(l), kLinear, kLog); } { fHM->CreateCanvas("dea_minv_oacut", "dea_minv_oacut", 1000, 800); vector h = fHM->GetH1ByNames("fhMinvOACutSignal", "fhMinvOACutGamma", "fhMinvOACutPi0", "fhMinvOACutBg"); const string l[] = {"Signal", "#gamma", "#pi^{0}", "BG"}; DrawH1(h, ToVector(l), kLinear, kLog); } { fHM->CreateCanvas("dea_minv_bg", "dea_minv_bg", 1000, 800); vector h; h.push_back((TH1D*)fHM->H1("fhMinvBg")->Clone()); h.push_back((TH1D*)fHM->H1("fhMinvMomCutBg")->Clone()); h.push_back((TH1D*)fHM->H1("fhMinvOACutBg")->Clone()); const string l[] = {"Lepton cut", "+Momentum cut", "+OAngle cut"}; DrawH1(h, ToVector(l)); gPad->SetLogy(true); } { TCanvas* c = fHM->CreateCanvas("dea_opening_angle_vs_ring_distance", "dea_opening_angle_vs_ring_distance", 1200, 1000); c->Divide(2, 2); c->cd(1); DrawH2(fHM->H2("fhOpenAngleVsRingDistBg")); DrawTextOnPad("BG", 0.6, 0.9, 0.75, 0.99); c->cd(2); DrawH2(fHM->H2("fhOpenAngleVsRingDistSignal")); DrawTextOnPad("Signal e^{#pm}", 0.6, 0.9, 0.75, 0.99); c->cd(3); DrawH2(fHM->H2("fhOpenAngleVsRingDistGamma")); DrawTextOnPad("#gamma e^{#pm}", 0.6, 0.9, 0.75, 0.99); c->cd(4); DrawH2(fHM->H2("fhOpenAngleVsRingDistPi0")); DrawTextOnPad("#pi^{0} e^{#pm}", 0.6, 0.9, 0.75, 0.99); } { fHM->CreateCanvas("dea_opening_angle", "dea_opening_angle", 1000, 800); vector h = fHM->GetH1ByNames("fhOpenAngleSignal", "fhOpenAngleGamma", "fhOpenAnglePi0", "fhOpenAngleBg"); fHM->NormalizeToIntegral(h); const string l[] = {"Signal e^{#pm}", "#gamma e^{#pm}", "#pi^{0} e^{#pm}", "BG"}; DrawH1(h, ToVector(l)); gPad->SetLogy(true); } { TCanvas* c = fHM->CreateCanvas("dea_opening_angle_signal_to_bg", "dea_opening_angle_signal_to_bg", 1400, 500); c->Divide(2, 1); c->cd(1); vector h = fHM->GetH1ByNames("fhOpenAngleSignalPositronToBG", "fhOpenAngleSignalPositronToBGLeptonCut", "fhOpenAngleSignalPositronToBGMomCut", "fhOpenAngleSignalPositronToBGOACut"); const string l[] = {"All Signal e^{+}", "+LeptonCut", "+MomCut", "+OACut"}; DrawH1(h, ToVector(l), kLinear, kLog); c->cd(2); vector h2 = fHM->GetH1ByNames("fhOpenAngleSignalElectronToBG", "fhOpenAngleSignalElectronToBGLeptonCut", "fhOpenAngleSignalElectronToBGMomCut", "fhOpenAngleSignalElectronToBGOACut"); const string l2[] = {"All Signal e^{-}", "+LeptonCut", "+MomCut" , "+OACut"}; DrawH1(h2, ToVector(l2), kLinear, kLog); } { fHM->CreateCanvas("dea_opening_angle_kine", "dea_opening_angle_kine", 1000, 800); vector h = fHM->GetH1ByNames("fhOpenAngleKineSignalEl", "fhOpenAngleKineGammaEl", "fhOpenAngleKinePi0El"); fHM->NormalizeToIntegral(h); const string l[] = {"Signal e^{#pm}", "#gamma e^{#pm}", "#pi^{0} e^{#pm}"}; DrawH1(h, ToVector(l), kLinear, kLog); } { fHM->CreateCanvas("dea_mom_kine", "dea_mom_kine", 1000, 800); vector h = fHM->GetH1ByNames("fhMomKineSignalEl", "fhMomKineGammaEl", "fhMomKinePi0El"); fHM->NormalizeToIntegral(h); const string l[] = {"Signal e^{#pm}", "#gamma e^{#pm}", "#pi^{0} e^{#pm}"}; DrawH1(h, ToVector(l), kLinear, kLog); } { fHM->CreateCanvas("dea_transmom_kine", "dea_transmom_kine", 1000, 800); vector h = fHM->GetH1ByNames("fhTransMomKineSignalEl", "fhTransMomKineGammaEl", "fhTransMomKinePi0El"); fHM->NormalizeToIntegral(h); const string l[] = {"Signal e^{#pm}", "#gamma e^{#pm}","#pi^{0} e^{#pm}"}; DrawH1(h, ToVector(l), kLinear, kLog); } { fHM->CreateCanvas("dea_min_distance_vs_opening_angle", "dea_min_distance_vs_opening_angle", 1000, 800); DrawH2(fHM->H2("fhMinDistVsMinAngle")); } { TCanvas* c = fHM->CreateCanvas("dea_kine_id", "dea_kine_id", 1000, 800); DrawH1(fHM->H1("fhKineID"), kLinear, kLog); } { TCanvas* c = fHM->CreateCanvas("dea_kine_ creation_mechanism", "dea_kine_ creation_mechanism", 800, 500); DrawH1(fHM->H1("fhKineMechanism"), kLinear, kLog); } { fHM->CreateCanvas("dea_ring_distance_pair", "dea_ring_distance_pair", 1000, 800); vector h = fHM->GetH2ProjectionYByNames("fhRingDistVsMomPairSignal", "fhRingDistVsMomPairGamma", "fhRingDistVsMomPairPi0", "fhRingDistVsMomPairBg"); fHM->NormalizeToIntegral(h); const string l[] = {"Signal", "#gamma", "#pi^{0}", "BG"}; DrawH1(h, ToVector(l), kLinear, kLog); } { TCanvas* c = fHM->CreateCanvas("dea_ring_distance_vs_momentum_pair", "dea_ring_distance_vs_momentum_pair", 1200, 1000); c->Divide(2,2); c->cd(1); DrawH2WithProfile(fHM->H2("fhRingDistVsMomPairSignal")); DrawTextOnPad("Signal e^{#pm}", 0.6, 0.9, 0.75, 0.99); c->cd(2); DrawH2WithProfile(fHM->H2("fhRingDistVsMomPairBg")); DrawTextOnPad("BG", 0.6, 0.9, 0.75, 0.99); c->cd(3); DrawH2WithProfile(fHM->H2("fhRingDistVsMomPairGamma")); DrawTextOnPad("#gamma e^{#pm}", 0.6, 0.9, 0.75, 0.99); c->cd(4); DrawH2WithProfile(fHM->H2("fhRingDistVsMomPairPi0")); DrawTextOnPad("#pi^{0} e^{#pm}", 0.6, 0.9, 0.75, 0.99); } } Bool_t HDEAnaTask::finalize() { drawHist(); fHM->SaveCanvasToImage(fOutputDir); TFile* outputFile = new TFile(fOutputFileName.c_str(), "RECREATE"); if(outputFile != NULL) { fHM->WriteToFile(); outputFile->Save(); outputFile->Close(); } return kTRUE; }