diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmCand.h b/analysis/PWGDIL/dielectron/lmvm/LmvmCand.h index 575bef1faa82a7f8469d92246067ce291e28a178..b7a291f50fa9ef540f63e1083aa5535d9f3a0ef5 100644 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmCand.h +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmCand.h @@ -49,7 +49,7 @@ public: TVector3 fPosition; TVector3 fMomentum; double fMass = 0.; - double fMassSig = -1.; // mass of signal mother (remains '-1' if candidate is not signal electron) + double fMassSig = -1.; // mass of signal mother (remains '-1' if candidate is not signal electron) double fEnergy = 0.; double fRapidity = 0.; int fCharge = 0; @@ -58,14 +58,14 @@ public: double fChi2Rich = 0.; double fChi2Trd = 0.; double fChi2Tof = 0.; - double fLength = 0.; // length of according global track - double fTime = 0.; // - double fTofDist = -1.; // Hit-Track-Distance in ToF + double fLength = 0.; // length of according global track + double fTime = 0.; // + double fTofDist = -1.; // Hit-Track-Distance in ToF // To investigate misidentifications in single sub detectors bool fIsRichElectron = false; bool fIsTrdElectron = false; - bool fIsTofElectron = false; + bool fIsTofElectron = false; int fMcMotherId = -1; int fEventNumber = 0; @@ -76,8 +76,8 @@ public: int fStsInd = -1; int fRichInd = -1; int fTrdInd = -1; - int fTofInd = -1;// TODO: change to "fTofHitInd" - int fTofTrackInd = -1; + int fTofInd = -1; // TODO: change to "fTofHitInd" + int fTofTrackInd = -1; int fMcPdg = -1; int fGTrackInd = -1; double fRichAnn = 0.; @@ -89,7 +89,7 @@ public: // Cuts. If true then cut is passed bool fIsChi2Prim = false; bool fIsElectron = false; - bool fIsGammaCut = true; // Will be set to 'false' as soon as a partner with minv < 25 MeV is found + bool fIsGammaCut = true; // Will be set to 'false' as soon as a partner with minv < 25 MeV is found bool fIsMvd1Cut = false; bool fIsMvd2Cut = false; bool fIsTtCut = false; diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmDef.h b/analysis/PWGDIL/dielectron/lmvm/LmvmDef.h index 809283fd9aa7fe14534f12ca8c7d1d5d7ab10c82..f3da4ee19491a992377a4f633c35e39618d2a4b4 100644 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmDef.h +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmDef.h @@ -125,16 +125,11 @@ public: class LmvmLegend { public: - LmvmLegend() {} - LmvmLegend(TH1D* h, const char* name, Option_t* opt) - : fH(h) - , fName(name) - , fOpt(opt) - { - } - TH1D* fH = nullptr; + LmvmLegend() {} + LmvmLegend(TH1D* h, const char* name, Option_t* opt) : fH(h), fName(name), fOpt(opt) {} + TH1D* fH = nullptr; const char* fName = ""; - Option_t* fOpt = "l"; + Option_t* fOpt = "l"; }; #endif diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmDraw.cxx b/analysis/PWGDIL/dielectron/lmvm/LmvmDraw.cxx index a4b7d5173c9178bb631b441df1f0292af37fdd90..bb8e3b95d234ee942f3a7048d63a52dca1832836 100644 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmDraw.cxx +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmDraw.cxx @@ -85,7 +85,7 @@ void LmvmDraw::DrawHistFromFile(const string& fileName, const string& outputDir, DrawMismatchesAndGhosts(); DrawMvdCutQa(); DrawMvdAndStsHist(); - DrawAccRecVsMom(); // TODO: finish this method! + DrawAccRecVsMom(); // TODO: finish this method! DrawPmtXY(); DrawMinvBg(); // TODO: do not extra method DrawBetaMomSpectra(); @@ -127,7 +127,7 @@ void LmvmDraw::DrawBetaMomSpectra() TH2D* hPos = fH.H2Clone("hBetaMom_cands_plutoEl+"); TH2D* hEl = fH.H2Clone("hBetaMom_cands_plutoEl-"); TCanvas* c = fH.fHM.CreateCanvas("betaMom/", "betaMom/", 1600, 800); - c->Divide(2,1); + c->Divide(2, 1); c->cd(1); DrawH2(hEl, kLinear, kLinear, kLog, "colz"); c->cd(2); @@ -390,7 +390,8 @@ void LmvmDraw::Draw1DCut(const string& hist, const string& sigOption, double cut DrawCutEffH1(hist, sigOption); c->cd(3); - TH1D* sign = fH.CreateSignificanceH1(fH.H1(hist, ELmvmSrc::Signal), fH.H1(hist, ELmvmSrc::Bg), hist + "_significance", sigOption); + TH1D* sign = fH.CreateSignificanceH1(fH.H1(hist, ELmvmSrc::Signal), fH.H1(hist, ELmvmSrc::Bg), hist + "_significance", + sigOption); DrawH1(sign, kLinear, kLinear, "hist"); } @@ -798,21 +799,22 @@ void LmvmDraw::DrawAccRecVsMom() // Acceptance and reconstruction yields cs. momentum for various detector combinations for (const int& pdg : {11, 211, 2212, 321}) { vector<string> subNames {"mc", "acc", "recSts", "recStsRich", "recStsRichTrd", "recStsRichTrdTof"}; - vector<string> latex {"MC", "Acc", "Rec in STS", "Rec in STS-RICH", "Rec in STS-RICH-TRD", "Rec in STS-RICH-TRD-TOF"}; + vector<string> latex { + "MC", "Acc", "Rec in STS", "Rec in STS-RICH", "Rec in STS-RICH-TRD", "Rec in STS-RICH-TRD-TOF"}; vector<string> latexAll(latex.size()), latexPrim(latex.size()); string ptcl = (pdg == 11) ? "hEl" : (pdg == 211) ? "hPi" : (pdg == 2212) ? "hProton" : "hKaon"; - + vector<TH1*> histsAll, histsPrim; int i = 0; - + for (const string& subName : subNames) { - TH1D* hAll = fH.H1(ptcl + "Mom_all_" + subName); + TH1D* hAll = fH.H1(ptcl + "Mom_all_" + subName); hAll->SetMinimum(3e-6); hAll->SetMaximum(50); latexAll[i] = latex[i] + " (" + Cbm::NumberToString(hAll->GetEntries() / fNofEvents, 2) + "/ev.)"; histsAll.push_back(hAll); - TH1D* hPrim = fH.H1(ptcl + "Mom_prim_" + subName); + TH1D* hPrim = fH.H1(ptcl + "Mom_prim_" + subName); hPrim->SetMinimum(3e-6); hPrim->SetMaximum(50); latexPrim[i] = latex[i] + " (" + Cbm::NumberToString(hPrim->GetEntries() / fNofEvents, 2) + "/ev.)"; @@ -821,7 +823,7 @@ void LmvmDraw::DrawAccRecVsMom() } //if (pdg == 321) continue; TODO: with kaons? - double y1 = 0.17; //(pdg == 211) ? 0.20 : 0.74; + double y1 = 0.17; //(pdg == 211) ? 0.20 : 0.74; double y2 = 0.42; //(pdg == 211) ? 0.45 : 0.99; string cName = "AccRecMom/" + ptcl + "Mom"; fH.fHM.CreateCanvas(cName, cName, 900, 900); diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.cxx b/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.cxx index 9828b8866b331b16c2a9b34aef46b5580c6c190e..5967d6725af85351df3284dda7fbb9d96420bef4 100755 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.cxx +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.cxx @@ -15,13 +15,13 @@ #include "TEllipse.h" #include "TF1.h" #include "TFile.h" +#include "TGraph.h" #include "TH1.h" #include "TH1D.h" #include "TH2D.h" #include "TKey.h" -#include "TLine.h" -#include "TGraph.h" #include "TLatex.h" +#include "TLine.h" #include "TMath.h" #include "TStyle.h" #include "TSystem.h" @@ -78,7 +78,7 @@ void LmvmDrawAll::DrawHistFromFile(const string& fileInmed, const string& fileQg DrawMinvPtAll(); DrawMinvBgSourcesAll(); DrawSBgVsMinv(); - InvestigateMisid(); // TODO: move some procedures in here into seperate methods? + InvestigateMisid(); // TODO: move some procedures in here into seperate methods? DrawBetaMomSpectra(); DrawMomRecoPrecision(); DrawMomPluto(); @@ -135,17 +135,16 @@ void LmvmDrawAll::CreateMeanHistAll() CreateMeanHist<TH2D>("hIndexStsTof_" + ptcl, nofEvents);*/ CreateMeanHist<TH2D>("hRichRingTrackDist_gTracks_" + ptcl, nofEvents); CreateMeanHist<TH2D>("hMatchId_gTracks_" + ptcl, nofEvents); - CreateMeanHist<TH2D>("hTofM2Calc_gTracks_" + ptcl, nofEvents); + CreateMeanHist<TH2D>("hTofM2Calc_gTracks_" + ptcl, nofEvents); for (const string& cat : {"trueid", "misid", "truematch", "mismatch"}) { CreateMeanHist<TH2D>("hTofHitXY_" + cat + "_" + ptcl, nofEvents); - CreateMeanHist<TH2D>("hTofPointXY_" + cat + "_" + ptcl, nofEvents); - CreateMeanHist<TH1D>("hTofHitPointDist_" + cat + "_" + ptcl, nofEvents); + CreateMeanHist<TH2D>("hTofPointXY_" + cat + "_" + ptcl, nofEvents); + CreateMeanHist<TH1D>("hTofHitPointDist_" + cat + "_" + ptcl, nofEvents); CreateMeanHist<TH2D>("hTofTimeVsMom_gTracks_" + cat + "_" + ptcl, nofEvents); CreateMeanHist<TH2D>("hTofHitTrackDist_gTracks_" + cat + "_" + ptcl, nofEvents); - } - + for (const string& cat : {"all", "complete"}) { CreateMeanHist<TH1D>("hNofMismatches_gTracks_" + cat + "_" + ptcl, nofEvents); CreateMeanHist<TH1D>("hNofMismatchedTrackSegments_" + cat + "_" + ptcl, nofEvents); @@ -155,7 +154,7 @@ void LmvmDrawAll::CreateMeanHistAll() } } } - } // fGTrackNames + } // fGTrackNames // Candidate Loop for (const string& ptcl : fHMean.fCandNames) { @@ -173,7 +172,7 @@ void LmvmDrawAll::CreateMeanHistAll() CreateMeanHist<TH2D>(fHMean.GetName("hPtY_cands_" + ptcl, step), nofEvents); CreateMeanHist<TH2D>(fHMean.GetName("hTofM2_cands_" + ptcl, step), nofEvents); CreateMeanHist<TH2D>(fHMean.GetName("hRichRingTrackDist_cands_" + ptcl, step), nofEvents); - } // steps + } // steps for (const string& cat : {"gTracks", "cands"}) { CreateMeanHist<TH1D>("hBetaMom_" + cat + "_" + ptcl, nofEvents); @@ -182,15 +181,15 @@ void LmvmDrawAll::CreateMeanHistAll() for (const string det : {"sts", "rich", "trd", "tof"}) { CreateMeanHist<TH2D>("hChi2VsMom_" + det + "_" + ptcl, nofEvents); - CreateMeanHist<TH2D>("hTofTimeVsChi2_" + det + "_" + ptcl, nofEvents); + CreateMeanHist<TH2D>("hTofTimeVsChi2_" + det + "_" + ptcl, nofEvents); } - + for (const string det : {"StsRich", "StsTrd", "RichTrd"}) { - CreateMeanHist<TH2D>("hChi2Comb_" + det + "_" + ptcl, nofEvents); + CreateMeanHist<TH2D>("hChi2Comb_" + det + "_" + ptcl, nofEvents); } - } // fCandNames + } // fCandNames // Step Loop for (auto step : fHMean.fAnaSteps) { @@ -217,16 +216,16 @@ void LmvmDrawAll::CreateMeanHistAll() CreateMeanHist<TH1D>(fHMean.GetName("hVertexYZ_misidTof", step), nofEvents); CreateMeanHist<TH1D>(fHMean.GetName("hVertexXY_misidTof", step), nofEvents); CreateMeanHist<TH1D>(fHMean.GetName("hVertexRZ_misidTof", step), nofEvents); - } // steps + } // steps string hBgSrc = "hMinvBgSource2_elid_"; for (const string& pair : {"gg", "pipi", "pi0pi0", "oo", "gpi", "gpi0", "go", "pipi0", "pio", "pi0o"}) { string hPairFull = hBgSrc + pair; CreateMeanHist<TH1D>(hPairFull, nofEvents, fRebMinv); } - + for (const string& cat : {"true", "misid"}) { - for (size_t iP = 4; iP < fHMean.fCandNames.size(); iP++ ) { + for (size_t iP = 4; iP < fHMean.fCandNames.size(); iP++) { CreateMeanHist<TH1D>("hMom_" + fHMean.fCandNames[iP] + "_" + cat, nofEvents); CreateMeanHist<TH1D>("hPtY_" + fHMean.fCandNames[iP] + "_" + cat, nofEvents); CreateMeanHist<TH1D>("hTofM2_" + fHMean.fCandNames[iP] + "_" + cat, nofEvents); @@ -249,7 +248,10 @@ void LmvmDrawAll::CreateMeanHistAll() CreateMeanHist<TH2D>("hVertexGTrackRZ", nofEvents); } -TH1D* LmvmDrawAll::GetCocktailMinvH1(const string& name, ELmvmAnaStep step) { return GetCocktailMinv<TH1D>(name, step); } +TH1D* LmvmDrawAll::GetCocktailMinvH1(const string& name, ELmvmAnaStep step) +{ + return GetCocktailMinv<TH1D>(name, step); +} template<class T> T* LmvmDrawAll::GetCocktailMinv(const string& name, ELmvmAnaStep step) @@ -270,7 +272,8 @@ T* LmvmDrawAll::GetCocktailMinv(const string& name, ELmvmAnaStep step) sHist->Scale(1. / binWidth); } if (cocktail == nullptr) cocktail = sHist; - else cocktail->Add(sHist); + else + cocktail->Add(sHist); } cocktail->Add(sEta); cocktail->Add(sPi0); @@ -284,8 +287,8 @@ void LmvmDrawAll::DrawMinvScaleValues() TH1D* qgp = new TH1D("qgp", "qgp; M_{ee} [GeV/c^{2}]; Scale Value", 3400., 0., 3.4); for (int iB = 1; iB <= 3400; iB++) { double binCenter = inmed->GetBinCenter(iB); - double sInmed = LmvmUtils::GetMassScaleInmed(binCenter); - double sQgp = LmvmUtils::GetMassScaleQgp(binCenter); + double sInmed = LmvmUtils::GetMassScaleInmed(binCenter); + double sQgp = LmvmUtils::GetMassScaleQgp(binCenter); inmed->SetBinContent(iB, sInmed); qgp->SetBinContent(iB, sQgp); } @@ -298,19 +301,19 @@ void LmvmDrawAll::DrawMinvScaleValues() qgp2->GetXaxis()->SetRangeUser(1., 1.1); inmed2->GetYaxis()->SetRangeUser(3e-3, 5e-1); qgp2->GetYaxis()->SetRangeUser(3e-3, 5e-1); - - fHMean.fHM.CreateCanvas("MinvScaleValues2", "MinvScaleValues2", 2400, 800); - DrawH1({inmed2, qgp2}, {"inmed", "QGP"}, kLinear, kLog, true, 0.7, 0.7, 0.9, 0.9, "p"); + fHMean.fHM.CreateCanvas("MinvScaleValues2", "MinvScaleValues2", 2400, 800); + DrawH1({inmed2, qgp2}, {"inmed", "QGP"}, kLinear, kLog, true, 0.7, 0.7, 0.9, 0.9, "p"); } -void LmvmDrawAll::DrawSignificance(TH2D* hEl, TH2D* hBg, const string& name, double minZ, double maxZ, const string& option) +void LmvmDrawAll::DrawSignificance(TH2D* hEl, TH2D* hBg, const string& name, double minZ, double maxZ, + const string& option) { string hElProjName = name + "_yProjEl"; string hBgProjName = name + "_yProjBg"; - TH1D* el = hEl->ProjectionY(hElProjName.c_str(), 1, hEl->GetXaxis()->GetNbins(), ""); - TH1D* bg = hBg->ProjectionY(hBgProjName.c_str(), 1, hBg->GetXaxis()->GetNbins(), ""); - TH2D* rat = (TH2D*) hEl->Clone(); + TH1D* el = hEl->ProjectionY(hElProjName.c_str(), 1, hEl->GetXaxis()->GetNbins(), ""); + TH1D* bg = hBg->ProjectionY(hBgProjName.c_str(), 1, hBg->GetXaxis()->GetNbins(), ""); + TH2D* rat = (TH2D*) hEl->Clone(); rat->Divide(hBg); const string hName = name + "_sign"; @@ -348,10 +351,11 @@ void LmvmDrawAll::DrawSignificance(TH2D* hEl, TH2D* hBg, const string& name, dou } } DrawTextOnPad("max. at: " + Cbm::NumberToString(maxX, 1), 0.55, 0.2, 0.85, 0.5); - DrawTextOnPad("Significance", 0.25, 0.9, 0.75, 0.99); + DrawTextOnPad("Significance", 0.25, 0.9, 0.75, 0.99); } -void LmvmDrawAll::DrawSignificancesAll() // TODO: implement automatization to seperate electrons from not-electrons in gTracks and cands +void LmvmDrawAll:: + DrawSignificancesAll() // TODO: implement automatization to seperate electrons from not-electrons in gTracks and cands { // RICH ANN and Electron Likelihood { @@ -368,7 +372,7 @@ void LmvmDrawAll::DrawSignificancesAll() // TODO: implement automatization to s el->Add(fHMean.H2("hChi2VsMom_trd_plutoEl-")); TH2D* bg = fHMean.H2Clone("hChi2VsMom_trd_pion+"); - for (size_t iP = 2; iP < fHMean.fCandNames.size(); iP++) { // consider UrQMD-electrons as "BG" + for (size_t iP = 2; iP < fHMean.fCandNames.size(); iP++) { // consider UrQMD-electrons as "BG" string ptcl = fHMean.fCandNames[iP]; string hName1 = "hChi2VsMom_trd_" + ptcl; bg->Add(fHMean.H2(hName1.c_str())); @@ -385,10 +389,12 @@ void LmvmDrawAll::DrawSignificancesAll() // TODO: implement automatization to s for (size_t iP = 1; iP < fHMean.fCandNames.size(); iP++) { string hName = "hTofTimeVsMom_cands_" + fHMean.fCandNames[iP]; if (iP <= 3) el->Add(fHMean.H2(hName.c_str())); - else if (iP == 4) continue; - else bg->Add(fHMean.H2(hName.c_str())); + else if (iP == 4) + continue; + else + bg->Add(fHMean.H2(hName.c_str())); } - + DrawSignificance(el, bg, "hTofTimeVsMom", 1e-9, 1e-2, "right"); } @@ -399,11 +405,13 @@ void LmvmDrawAll::DrawSignificancesAll() // TODO: implement automatization to s for (size_t iP = 1; iP < fHMean.fCandNames.size(); iP++) { string hName = "hTofHitTrackDist_cands_" + fHMean.fCandNames[iP]; - if (iP <= 1) el->Add(fHMean.H2(hName.c_str())); // consider UrQMD-electrons as "BG" - else if (iP == 4) continue; - else bg->Add(fHMean.H2(hName.c_str())); + if (iP <= 1) el->Add(fHMean.H2(hName.c_str())); // consider UrQMD-electrons as "BG" + else if (iP == 4) + continue; + else + bg->Add(fHMean.H2(hName.c_str())); } - + DrawSignificance(el, bg, "hTofHitTrackDist", 1e-9, 1e-2, "right"); } } @@ -430,8 +438,8 @@ void LmvmDrawAll::DrawTofHitXY() //double min = 1e-7; //double max = 2e-3; double minD = 1e-7; - double maxD = 2e-2; - + double maxD = 2e-2; + /*TCanvas* c = fHMean.fHM.CreateCanvas("ToF/point-hit-dist_gTracks", "ToF/point-hit-dist_gTracks", 2400, 2400); // TODO: do these as loop for pile and not-pile histograms c->Divide(4, 4); int i = 1; @@ -475,28 +483,30 @@ void LmvmDrawAll::DrawTofHitXY() // draw only for particles in ToF pile that are identified as electrons in ToF double minH = 2e-6; double maxH = 1e-4; - - TCanvas* c3 = fHMean.fHM.CreateCanvas("ToF/HitXY/MisidsInTofPile/point-hit-dist_all", "ToF/HitXY/MisidsInTofPile/point-hit-dist_all", 2400, 1600); + + TCanvas* c3 = fHMean.fHM.CreateCanvas("ToF/HitXY/MisidsInTofPile/point-hit-dist_all", + "ToF/HitXY/MisidsInTofPile/point-hit-dist_all", 2400, 1600); c3->Divide(3, 2); int iC = 1; - for (size_t iP = 4; iP < fHMean.fCandNames.size(); iP++ ) { // only draw for not-electrons + for (size_t iP = 4; iP < fHMean.fCandNames.size(); iP++) { // only draw for not-electrons c3->cd(iC++); - string ptcl = fHMean.fCandNames[iP]; + string ptcl = fHMean.fCandNames[iP]; string hName = "hTofPileHitPointDist_" + ptcl; - TH1D* h = fHMean.H1Clone(hName.c_str()); + TH1D* h = fHMean.H1Clone(hName.c_str()); h->GetYaxis()->SetRangeUser(minD, maxD); h->Fit("gaus", "Q", "", 0., 2.); DrawH1(h, kLinear, kLog, "hist"); DrawTextOnPad(fHMean.fCandLatex[iP], 0.4, 0.9, 0.54, 0.99); - TF1* func = h->GetFunction("gaus"); + TF1* func = h->GetFunction("gaus"); double mean = (func != nullptr) ? func->GetParameter("Mean") : 0.; DrawTextOnPad("mean: " + Cbm::NumberToString(mean, 2), 0.2, 0.65, 0.6, 0.89); } - for (size_t iP = 4; iP < fHMean.fCandNames.size(); iP++ ) { // only draw for not-electrons + for (size_t iP = 4; iP < fHMean.fCandNames.size(); iP++) { // only draw for not-electrons string ptcl = fHMean.fCandNames[iP]; - TCanvas* c4 = fHMean.fHM.CreateCanvas("ToF/HitXY/MisidsInTofPile/" + ptcl, "ToF/HitXY/MisidsInTofPile/" + ptcl, 2400, 800); - c4->Divide(3,1); + TCanvas* c4 = + fHMean.fHM.CreateCanvas("ToF/HitXY/MisidsInTofPile/" + ptcl, "ToF/HitXY/MisidsInTofPile/" + ptcl, 2400, 800); + c4->Divide(3, 1); c4->cd(1); TH2D* hit = fHMean.H2Clone("hTofPileHitXY_" + ptcl); hit->GetZaxis()->SetRangeUser(minH, maxH); @@ -513,7 +523,7 @@ void LmvmDrawAll::DrawTofHitXY() dist->GetYaxis()->SetRangeUser(minD, maxD); DrawH1(dist, kLinear, kLog, "hist"); DrawTextOnPad("Distance Hit-Point", 0.3, 0.9, 0.7, 0.99); - TF1* func = dist->GetFunction("gaus"); + TF1* func = dist->GetFunction("gaus"); double mean = (func != nullptr) ? func->GetParameter("Mean") : 0.; DrawTextOnPad("mean: " + Cbm::NumberToString(mean, 1), 0.2, 0.65, 0.5, 0.89); } @@ -524,8 +534,8 @@ void LmvmDrawAll::DrawBetaMomSpectra() double min = 1e-5; double max = 10; for (const string& cat : {"gTracks", "cands"}) { - for (size_t i = 0; i < fHMean.fCandNames.size(); i++ ) { - string cName = "BetaMom/" + cat + "/" + fHMean.fCandNames[i]; + for (size_t i = 0; i < fHMean.fCandNames.size(); i++) { + string cName = "BetaMom/" + cat + "/" + fHMean.fCandNames[i]; string hName = "hBetaMom_" + cat + "_" + fHMean.fCandNames[i]; fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 800, 800); TH2D* h = fHMean.H2Clone(hName.c_str()); @@ -544,15 +554,15 @@ void LmvmDrawAll::DrawMomPluto() { for (const string& pi : {"Px", "Py", "Pz"}) { string cName = "MomDistPluto/" + pi; - TCanvas* c = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 2400, 1600); - c->Divide(3,2); + TCanvas* c = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 2400, 1600); + c->Divide(3, 2); int i = 1; for (ELmvmSignal signal : fHMean.fSignals) { c->cd(i++); string hName = "hMom" + pi; - TH1D* ep = H(signal)->H1Clone((hName + "+_signal_mc").c_str()); - TH1D* em = H(signal)->H1Clone((hName + "-_signal_mc").c_str()); - double bW = ep->GetBinWidth(1); + TH1D* ep = H(signal)->H1Clone((hName + "+_signal_mc").c_str()); + TH1D* em = H(signal)->H1Clone((hName + "-_signal_mc").c_str()); + double bW = ep->GetBinWidth(1); ep->Scale(1. / (H(signal)->H1("hEventNumber")->GetEntries() * bW)); em->Scale(1. / (H(signal)->H1("hEventNumber")->GetEntries() * bW)); ep->GetYaxis()->SetTitle("1/N dN/dP [GeV/c]^{-1}"); @@ -574,14 +584,15 @@ void LmvmDrawAll::DrawMomentum() TCanvas* c = fHMean.fHM.CreateCanvas("hMom/misid_vs_True", "hMom/misid_vs_True", 2400, 1600); c->Divide(3, 2); int iC = 1; - for (size_t iP = 4; iP < fHMean.fCandNames.size(); iP++ ) { // only draw for not-electrons + for (size_t iP = 4; iP < fHMean.fCandNames.size(); iP++) { // only draw for not-electrons c->cd(iC++); string ptcl = fHMean.fCandNames[iP]; TH1D* hMis = fHMean.H1Clone("hMom_" + ptcl + "_misid"); TH1D* hTrue = fHMean.H1Clone("hMom_" + ptcl + "_true"); hMis->GetYaxis()->SetRangeUser(min, max); hTrue->GetYaxis()->SetRangeUser(min, max); - DrawH1({hTrue, hMis}, {"not misidentified", "misidentified as electron"}, kLinear, kLog, true, 0.55, 0.8, 0.95, 0.91 ,"hist"); + DrawH1({hTrue, hMis}, {"not misidentified", "misidentified as electron"}, kLinear, kLog, true, 0.55, 0.8, 0.95, + 0.91, "hist"); DrawTextOnPad(fHMean.fCandLatex[iP], 0.4, 0.9, 0.54, 0.99); } } @@ -610,7 +621,8 @@ void LmvmDrawAll::DrawCutEffSignal() fHMean.fHM.CreateCanvas("CutEff_signal", "CutEff_signal", 900, 900); elid->GetYaxis()->SetTitle("Efficiency"); - DrawH1({elid, gamma, tt, pt}, {"El-ID", "Gamma cut", "TT cut", "P_{t} cut"}, kLinear, kLinear, true, 0.6, 0.75, 0.95, 0.91, "hist"); + DrawH1({elid, gamma, tt, pt}, {"El-ID", "Gamma cut", "TT cut", "P_{t} cut"}, kLinear, kLinear, true, 0.6, 0.75, 0.95, + 0.91, "hist"); DrawTextOnPad("Signal Efficiency", 0.3, 0.9, 0.7, 0.99); } @@ -632,7 +644,8 @@ void LmvmDrawAll::DrawPionSuppression() string text = fHMean.fAnaStepNames[static_cast<int>(step)]; fHMean.fHM.CreateCanvas("PionSuppression_" + text, "PionSuppression_" + text, 900, 900); elidMc->GetYaxis()->SetTitle("Pion Suppression"); - DrawH1({elidMc, gammaMc, ttMc, ptMc}, {"El-ID", "Gamma cut", "TT cut", "P_{t} cut"}, kLinear, kLog, true, 0.7, 0.8, 0.95, 0.91, "hist"); + DrawH1({elidMc, gammaMc, ttMc, ptMc}, {"El-ID", "Gamma cut", "TT cut", "P_{t} cut"}, kLinear, kLog, true, 0.7, 0.8, + 0.95, 0.91, "hist"); DrawTextOnPad("Pion Suppression (norm. to " + text + ")", 0.25, 0.88, 0.75, 0.99); } } @@ -646,20 +659,21 @@ void LmvmDrawAll::DrawTofM2() string hName = "hTofM2_gTracks_" + fHMean.fGTrackNames[iP]; hEl->Add(fHMean.H2(hName.c_str())); } - TH2D* hBg = fHMean.H2Clone("hTofM2_gTracks_"+ fHMean.fGTrackNames[4]); + TH2D* hBg = fHMean.H2Clone("hTofM2_gTracks_" + fHMean.fGTrackNames[4]); for (size_t iP = 5; iP < fHMean.fGTrackNames.size(); iP++) { string hName = "hTofM2_gTracks_" + fHMean.fGTrackNames[iP]; hBg->Add(fHMean.H2(hName.c_str())); } - - double minX = 0.; - double maxX = 2.5; + + double minX = 0.; + double maxX = 2.5; double minY = -0.05; - double maxY = 0.05; - double minZ = 1e-8; - double maxZ = 1; + double maxY = 0.05; + double minZ = 1e-8; + double maxZ = 1; vector<TLine*> lines {new TLine(0., 0.01, 1.3, 0.01), new TLine(1.3, 0.01, 2.5, 0.022)}; // set by hand - TCanvas* c = fHMean.fHM.CreateCanvas("ToF/Purity/gTracks", "ToF/Purity/gTracks", 1600, 800); // TODO: check: is this all gTracks? + TCanvas* c = fHMean.fHM.CreateCanvas("ToF/Purity/gTracks", "ToF/Purity/gTracks", 1600, + 800); // TODO: check: is this all gTracks? c->Divide(2, 1); c->cd(1); hEl->GetXaxis()->SetRangeUser(minX, maxX); @@ -672,7 +686,7 @@ void LmvmDrawAll::DrawTofM2() lines[i]->Draw(); } c->cd(2); - hBg->GetXaxis()->SetRangeUser(minX, maxX); + hBg->GetXaxis()->SetRangeUser(minX, maxX); hBg->GetYaxis()->SetRangeUser(minY, maxY); hBg->GetZaxis()->SetRangeUser(minZ, maxZ); DrawH2(hBg, kLinear, kLinear, kLog, "colz"); @@ -685,27 +699,27 @@ void LmvmDrawAll::DrawTofM2() // for candidates { - double minX = 0.; - double maxX = 2.5; + double minX = 0.; + double maxX = 2.5; double minY = -0.05; - double maxY = 0.05; - double minZ = 1e-8; - double maxZ = 10; + double maxY = 0.05; + double minZ = 1e-8; + double maxZ = 10; vector<TLine*> lines {new TLine(0., 0.01, 1.3, 0.01), new TLine(1.3, 0.01, 2.5, 0.022)}; // set by hand - + for (auto step : fHMean.fAnaSteps) { TH2D* hEl = fHMean.H2Clone(("hTofM2_cands_" + fHMean.fCandNames[0]), step); for (size_t iP = 1; iP < 4; iP++) { - string hName = "hTofM2_cands_" + fHMean.fCandNames[iP] + "_"+ fHMean.fAnaStepNames[static_cast<int>(step)]; + string hName = "hTofM2_cands_" + fHMean.fCandNames[iP] + "_" + fHMean.fAnaStepNames[static_cast<int>(step)]; hEl->Add(fHMean.H2(hName.c_str())); } TH2D* hBg = fHMean.H2Clone(("hTofM2_cands_" + fHMean.fCandNames[4]), step); for (size_t iP = 5; iP < fHMean.fCandNames.size(); iP++) { - string hName = "hTofM2_cands_" + fHMean.fCandNames[iP] + "_"+ fHMean.fAnaStepNames[static_cast<int>(step)]; + string hName = "hTofM2_cands_" + fHMean.fCandNames[iP] + "_" + fHMean.fAnaStepNames[static_cast<int>(step)]; hBg->Add(fHMean.H2(hName.c_str())); } string cName = "ToF/Purity/" + fHMean.fAnaStepNames[static_cast<int>(step)]; - TCanvas* c = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 1600, 800); + TCanvas* c = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 1600, 800); c->Divide(2, 1); c->cd(1); hEl->GetXaxis()->SetRangeUser(minX, maxX); @@ -734,36 +748,37 @@ void LmvmDrawAll::DrawTofM2() void LmvmDrawAll::DrawPtYAndTofM2Elid() { double x[200], y[200]; - int n = 200; - double m = 511e-6; - double m2 = m*m; - double p = 6; - double p2 = p*p; - for (int i=0; i < n; i++) { + int n = 200; + double m = 511e-6; + double m2 = m * m; + double p = 6; + double p2 = p * p; + for (int i = 0; i < n; i++) { x[i] = 0.02 * i; - double counter = std::sqrt(2 * m2 * TMath::Exp(2 * x[i]) - m2 * TMath::Exp(4 * x[i]) + 4 * p2 * TMath::Exp(2 * x[i]) - m2); - double denom = std::sqrt(1 + 2 * TMath::Exp(2 * x[i]) + TMath::Exp(4 * x[i])); - double r = counter / denom; - y[i] = r; + double counter = + std::sqrt(2 * m2 * TMath::Exp(2 * x[i]) - m2 * TMath::Exp(4 * x[i]) + 4 * p2 * TMath::Exp(2 * x[i]) - m2); + double denom = std::sqrt(1 + 2 * TMath::Exp(2 * x[i]) + TMath::Exp(4 * x[i])); + double r = counter / denom; + y[i] = r; } - auto pLine = new TGraph(n, x, y); // draw momentum line into Pt-Y histogram + auto pLine = new TGraph(n, x, y); // draw momentum line into Pt-Y histogram for (const string& cat : {"hPtY", "hTofM2"}) { // draw all global tracks seperate for diff. particles { - double min = 1e-8; - double max = 10; + double min = 1e-8; + double max = 10; string cName = cat + "/globalTracks/all"; - TCanvas* c = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 2400, 2400); + TCanvas* c = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 2400, 2400); c->Divide(4, 4); int i = 1; for (auto ptcl : fHMean.fGTrackNames) { c->cd(i++); string hName = cat + "_gTracks_" + ptcl; - TH2D* h = fHMean.H2Clone(hName.c_str()); + TH2D* h = fHMean.H2Clone(hName.c_str()); h->GetZaxis()->SetRangeUser(min, max); DrawH2(h, kLinear, kLinear, kLog, "colz"); - DrawTextOnPad(fHMean.fGTrackLatex[i-2], 0.25, 0.9, 0.75, 0.999); + DrawTextOnPad(fHMean.fGTrackLatex[i - 2], 0.25, 0.9, 0.75, 0.999); if (cat == "hPtY") pLine->Draw("C"); } } @@ -771,21 +786,21 @@ void LmvmDrawAll::DrawPtYAndTofM2Elid() // draw zoomed in TofM2 { if (cat == "hTofM2") { - double min = 1e-8; - double max = 10; + double min = 1e-8; + double max = 10; string cName = cat + "/globalTracks/all_zoom"; - TCanvas* c = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 2400, 2400); + TCanvas* c = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 2400, 2400); c->Divide(4, 4); int i = 1; for (auto ptcl : fHMean.fGTrackNames) { c->cd(i++); string hName = cat + "_gTracks_" + ptcl; - TH2D* h = fHMean.H2Clone(hName.c_str()); + TH2D* h = fHMean.H2Clone(hName.c_str()); h->GetXaxis()->SetRangeUser(0., 2.5); h->GetYaxis()->SetRangeUser(-0.05, 0.05); h->GetZaxis()->SetRangeUser(min, max); DrawH2(h, kLinear, kLinear, kLog, "colz"); - DrawTextOnPad(fHMean.fGTrackLatex[i-2], 0.25, 0.9, 0.75, 0.999); + DrawTextOnPad(fHMean.fGTrackLatex[i - 2], 0.25, 0.9, 0.75, 0.999); } } } @@ -793,42 +808,42 @@ void LmvmDrawAll::DrawPtYAndTofM2Elid() // draw misidentified and not-misidentified particles and ratio of both (all from global tracks) { double min = (cat == "hPtY") ? 1e-7 : 1e-8; - double max = (cat == "hPtY") ? 10 : 10; - for (size_t iP = 4; iP < fHMean.fCandNames.size(); iP++ ) { // only draw for not-electrons - string ptcl = fHMean.fCandNames[iP]; + double max = (cat == "hPtY") ? 10 : 10; + for (size_t iP = 4; iP < fHMean.fCandNames.size(); iP++) { // only draw for not-electrons + string ptcl = fHMean.fCandNames[iP]; string cName = cat + "/candidates/" + ptcl + "_misid"; - TCanvas* c = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 2400, 800); + TCanvas* c = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 2400, 800); c->Divide(3, 1); c->cd(1); TH2D* hTrue = fHMean.H2Clone(cat + "_" + ptcl + "_true"); - if (cat == "hTofM2") { + if (cat == "hTofM2") { hTrue->GetXaxis()->SetRangeUser(0., 2.5); - hTrue->GetYaxis()->SetRangeUser(-0.05, 0.05); - } + hTrue->GetYaxis()->SetRangeUser(-0.05, 0.05); + } hTrue->GetZaxis()->SetRangeUser(min, max); DrawH2(hTrue, kLinear, kLinear, kLog, "colz"); - DrawTextOnPad(fHMean.fCandLatex[iP] + " (not misidentified)" , 0.2, 0.9, 0.75, 0.99); + DrawTextOnPad(fHMean.fCandLatex[iP] + " (not misidentified)", 0.2, 0.9, 0.75, 0.99); if (cat == "hPtY") pLine->Draw("C"); c->cd(2); TH2D* hMis = fHMean.H2Clone(cat + "_" + ptcl + "_misid"); - if (cat == "hTofM2") { + if (cat == "hTofM2") { hMis->GetXaxis()->SetRangeUser(0., 2.5); - hMis->GetYaxis()->SetRangeUser(-0.05, 0.05); - } + hMis->GetYaxis()->SetRangeUser(-0.05, 0.05); + } hMis->GetZaxis()->SetRangeUser(min, max); DrawH2(hMis, kLinear, kLinear, kLog, "colz"); - DrawTextOnPad(fHMean.fCandLatex[iP] + " (misidentified as electron)" , 0.15, 0.9, 0.75, 0.99); + DrawTextOnPad(fHMean.fCandLatex[iP] + " (misidentified as electron)", 0.15, 0.9, 0.75, 0.99); if (cat == "hPtY") pLine->Draw("C"); - + c->cd(3); TH2D* hRat = fHMean.H2Clone(cat + "_" + ptcl + "_misid"); hRat->Divide(hTrue); - if (cat == "hTofM2") { + if (cat == "hTofM2") { hRat->GetXaxis()->SetRangeUser(0., 2.5); - hRat->GetYaxis()->SetRangeUser(-0.05, 0.05); - } - hRat->GetZaxis()->SetRangeUser(1e-5, 1); // TODO: check best min value + hRat->GetYaxis()->SetRangeUser(-0.05, 0.05); + } + hRat->GetZaxis()->SetRangeUser(1e-5, 1); // TODO: check best min value hRat->GetZaxis()->SetTitle("Ratio #misid./#not misid."); DrawH2(hRat, kLinear, kLinear, kLog, "colz"); DrawTextOnPad(fHMean.fCandLatex[iP] + ": Ratio #misid./#not misid.", 0.1, 0.9, 0.6, 0.99); @@ -840,9 +855,9 @@ void LmvmDrawAll::DrawPtYAndTofM2Elid() { double min = (cat == "hPtY") ? 2e-8 : 2e-8; double max = (cat == "hPtY") ? 5e-3 : 5e-3; - for (size_t iP = 4; iP < fHMean.fCandNames.size(); iP++ ) { - string ptcl = fHMean.fCandNames[iP]; - string cName = cat + "/candidates/" + ptcl + "_misid_steps"; + for (size_t iP = 4; iP < fHMean.fCandNames.size(); iP++) { + string ptcl = fHMean.fCandNames[iP]; + string cName = cat + "/candidates/" + ptcl + "_misid_steps"; TCanvas* cPty = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 1800, 1800); cPty->Divide(3, 3); int i = 1; @@ -850,28 +865,39 @@ void LmvmDrawAll::DrawPtYAndTofM2Elid() if (step < ELmvmAnaStep::ElId) continue; cPty->cd(i++); string hName = cat + "_cands_" + ptcl; - TH2D* h = fHMean.H2Clone(hName.c_str(), step); - if (cat == "hTofM2") { + TH2D* h = fHMean.H2Clone(hName.c_str(), step); + if (cat == "hTofM2") { h->GetXaxis()->SetRangeUser(0., 2.5); - h->GetYaxis()->SetRangeUser(-0.05, 0.05); - } - h->GetZaxis()->SetRangeUser(min, max); + h->GetYaxis()->SetRangeUser(-0.05, 0.05); + } + h->GetZaxis()->SetRangeUser(min, max); DrawH2(h, kLinear, kLinear, kLog, "colz"); fHMean.DrawAnaStepOnPad(step); if (cat == "hPtY") pLine->Draw("C"); } } } - } // cat: "hPtY", "hTofM2" + } // cat: "hPtY", "hTofM2" } void LmvmDrawAll::DrawTofPilePids() { // draw particles that occur in "ToF Pile" - vector<string> xLabel = {"e^{-}_{PLUTO}", "e^{+}_{PLUTO}", "e^{-}_{UrQMD_prim}", "e^{+}_{UrQMD_prim}", "e^{-}_{UrQMD_sec}", "e^{+}_{UrQMD_sec}", "#pi^{+}", "#pi^{-}", "p", "K^{+}", "K^{-}", "o."}; - double max = 3; - double min = 1e-5; - TCanvas* cTofPile = fHMean.fHM.CreateCanvas("ToF/TofPilePids", "ToF/TofPilePids", 1800, 1800); + vector<string> xLabel = {"e^{-}_{PLUTO}", + "e^{+}_{PLUTO}", + "e^{-}_{UrQMD_prim}", + "e^{+}_{UrQMD_prim}", + "e^{-}_{UrQMD_sec}", + "e^{+}_{UrQMD_sec}", + "#pi^{+}", + "#pi^{-}", + "p", + "K^{+}", + "K^{-}", + "o."}; + double max = 3; + double min = 1e-5; + TCanvas* cTofPile = fHMean.fHM.CreateCanvas("ToF/TofPilePids", "ToF/TofPilePids", 1800, 1800); cTofPile->Divide(4, 4); cTofPile->cd(1); @@ -912,7 +938,7 @@ void LmvmDrawAll::DrawPurity() for (size_t y = 1; y <= yLabel.size(); y++) { h->GetYaxis()->SetBinLabel(y, yLabel[y - 1].c_str()); } - double nEl = h->Integral(1, h->GetXaxis()->GetNbins(), 2, 2); // do not count PLUTO particles + double nEl = h->Integral(1, h->GetXaxis()->GetNbins(), 2, 2); // do not count PLUTO particles double purity = (nEl / h->Integral(1, h->GetXaxis()->GetNbins(), 2, h->GetYaxis()->GetNbins())) * 100; DrawTextOnPad("Purity: " + Cbm::NumberToString(purity, 1) + " %", 0.1, 0.9, 0.45, 0.99); } @@ -926,7 +952,8 @@ void LmvmDrawAll::DrawPurity() TH1D* purity = new TH1D("purity_Mom", "purity_Mom; P [GeV/c]; Purity [%]", nBins, xMin, xMax); for (int i = 1; i <= purity->GetNbinsX(); i++) { double nEl = fHMean.H2("hCandPdgVsMom", ELmvmAnaStep::ElId)->GetBinContent(i, 2); - double nAll = fHMean.H2("hCandPdgVsMom", ELmvmAnaStep::ElId)->Integral(i, i, 2, fHMean.H2("hCandPdgVsMom_elid")->GetYaxis()->GetNbins()); + double nAll = fHMean.H2("hCandPdgVsMom", ELmvmAnaStep::ElId) + ->Integral(i, i, 2, fHMean.H2("hCandPdgVsMom_elid")->GetYaxis()->GetNbins()); double val = (nAll != 0) ? 100 * nEl / nAll : 0.; purity->SetBinContent(i, val); } @@ -939,34 +966,34 @@ void LmvmDrawAll::DrawPurity() // Purity seperate for ID Detectors { vector<string> yLabel = {"#pi^{+}", "#pi^{-}", "p", "K^{+}", "K^{-}", "o."}; - double min = 1e-7; - double max = 10; + double min = 1e-7; + double max = 10; TH2D* rich1 = fHMean.H2Clone("hPdgVsMom_gTracks_rich_all"); TH2D* trd1 = fHMean.H2Clone("hPdgVsMom_gTracks_trd_all"); TH2D* tof1 = fHMean.H2Clone("hPdgVsMom_gTracks_tof_all"); TH2D* rich2 = fHMean.H2Clone("hPdgVsMom_gTracks_rich_complete"); TH2D* trd2 = fHMean.H2Clone("hPdgVsMom_gTracks_trd_complete"); - TH2D* tof2 = fHMean.H2Clone("hPdgVsMom_gTracks_tof_complete"); + TH2D* tof2 = fHMean.H2Clone("hPdgVsMom_gTracks_tof_complete"); rich1->GetZaxis()->SetRangeUser(min, max); - trd1 ->GetZaxis()->SetRangeUser(min, max); - tof1 ->GetZaxis()->SetRangeUser(min, max); + trd1->GetZaxis()->SetRangeUser(min, max); + tof1->GetZaxis()->SetRangeUser(min, max); rich2->GetZaxis()->SetRangeUser(min, max); - trd2 ->GetZaxis()->SetRangeUser(min, max); - tof2 ->GetZaxis()->SetRangeUser(min, max); + trd2->GetZaxis()->SetRangeUser(min, max); + tof2->GetZaxis()->SetRangeUser(min, max); for (size_t y = 1; y <= yLabel.size(); y++) { rich1->GetYaxis()->SetBinLabel(y, yLabel[y - 1].c_str()); - trd1 ->GetYaxis()->SetBinLabel(y, yLabel[y - 1].c_str()); - tof1 ->GetYaxis()->SetBinLabel(y, yLabel[y - 1].c_str()); + trd1->GetYaxis()->SetBinLabel(y, yLabel[y - 1].c_str()); + tof1->GetYaxis()->SetBinLabel(y, yLabel[y - 1].c_str()); rich2->GetYaxis()->SetBinLabel(y, yLabel[y - 1].c_str()); - trd2 ->GetYaxis()->SetBinLabel(y, yLabel[y - 1].c_str()); - tof2 ->GetYaxis()->SetBinLabel(y, yLabel[y - 1].c_str()); + trd2->GetYaxis()->SetBinLabel(y, yLabel[y - 1].c_str()); + tof2->GetYaxis()->SetBinLabel(y, yLabel[y - 1].c_str()); } TCanvas* c1 = fHMean.fHM.CreateCanvas("Misid/misid_gTracks_all", "Misid/misid_gTracks_all", 2400, 800); - c1->Divide(3,1); + c1->Divide(3, 1); c1->cd(1); DrawH2(rich1, kLinear, kLinear, kLog, "colz"); DrawTextOnPad("RICH", 0.25, 0.9, 0.75, 0.999); @@ -981,7 +1008,7 @@ void LmvmDrawAll::DrawPurity() DrawPurityHistText(tof1); TCanvas* c2 = fHMean.fHM.CreateCanvas("Misid/misid_gTracks_complete", "Misid/misid_gTracks_complete", 2400, 800); - c2->Divide(3,1); + c2->Divide(3, 1); c2->cd(1); DrawH2(rich2, kLinear, kLinear, kLog, "colz"); DrawTextOnPad("RICH", 0.25, 0.9, 0.75, 0.999); @@ -1002,22 +1029,22 @@ void LmvmDrawAll::DrawPurity() void LmvmDrawAll::DrawPurityHistText(TH2D* h) { int nX = h->GetXaxis()->GetNbins(); - + double xMin = 0.4; double xMax = 0.75; - double piP = h->Integral(1, nX, 1, 1); - double piM = h->Integral(1, nX, 2, 2); - double p = h->Integral(1, nX, 3, 3); - double KP = h->Integral(1, nX, 4, 4); - double KM = h->Integral(1, nX, 5, 5); - double o = h->Integral(1, nX, 6, 6); + double piP = h->Integral(1, nX, 1, 1); + double piM = h->Integral(1, nX, 2, 2); + double p = h->Integral(1, nX, 3, 3); + double KP = h->Integral(1, nX, 4, 4); + double KM = h->Integral(1, nX, 5, 5); + double o = h->Integral(1, nX, 6, 6); DrawTextOnPad("#pi^{+}: " + Cbm::NumberToString(piP, 1) + " /Ev", xMin, 0.12, xMax, 0.3); DrawTextOnPad("#pi^{-}: " + Cbm::NumberToString(piM, 1) + " /Ev", xMin, 0.28, xMax, 0.37); - DrawTextOnPad("p: " + Cbm::NumberToString(p, 1) + " /Ev", xMin, 0.43, xMax, 0.48); - DrawTextOnPad("K^{+}: " + Cbm::NumberToString(KP, 1) + " /Ev", xMin, 0.53, xMax, 0.65); - DrawTextOnPad("K^{-}: " + Cbm::NumberToString(KM, 1) + " /Ev", xMin, 0.66, xMax, 0.76); - DrawTextOnPad("o.: " + Cbm::NumberToString(o, 1) + " /Ev", xMin, 0.77, xMax, 0.9); + DrawTextOnPad("p: " + Cbm::NumberToString(p, 1) + " /Ev", xMin, 0.43, xMax, 0.48); + DrawTextOnPad("K^{+}: " + Cbm::NumberToString(KP, 1) + " /Ev", xMin, 0.53, xMax, 0.65); + DrawTextOnPad("K^{-}: " + Cbm::NumberToString(KM, 1) + " /Ev", xMin, 0.66, xMax, 0.76); + DrawTextOnPad("o.: " + Cbm::NumberToString(o, 1) + " /Ev", xMin, 0.77, xMax, 0.9); } void LmvmDrawAll::InvestigateMisid() @@ -1035,7 +1062,8 @@ void LmvmDrawAll::InvestigateMisid() } // RICH: Ring-Track Distance - fHMean.DrawAllGTracks(2, "hRichRingTrackDist/gTracks_all", "hRichRingTrackDist_gTracks", {""}, {""}, 1e-7, 1e-2); // TODO: move this anywhere else? + fHMean.DrawAllGTracks(2, "hRichRingTrackDist/gTracks_all", "hRichRingTrackDist_gTracks", {""}, {""}, 1e-7, + 1e-2); // TODO: move this anywhere else? fHMean.DrawAllCandsAndSteps(2, "hRichRingTrackDist/", "hRichRingTrackDist_cands", {""}, {""}, 1e-7, 1e-2); // ToF @@ -1045,29 +1073,32 @@ void LmvmDrawAll::InvestigateMisid() double minD = 3e-6; double maxD = 3e-1; for (const string& cat : {"trueid", "misid", "truematch", "mismatch"}) { - fHMean.DrawAllGTracks(2, "ToF/HitTrackDist/HitPointDist_" + cat, "hTofHitTrackDist_gTracks_" + cat, {""}, {""}, 1e-9, 1.); + fHMean.DrawAllGTracks(2, "ToF/HitTrackDist/HitPointDist_" + cat, "hTofHitTrackDist_gTracks_" + cat, {""}, {""}, + 1e-9, 1.); fHMean.DrawAllGTracks(2, "ToF/HitXY/PointXY_" + cat, "hTofPointXY_" + cat, {""}, {""}, minXY, maxXY); fHMean.DrawAllGTracks(2, "ToF/HitXY/HitXY_" + cat, "hTofHitXY_" + cat, {""}, {""}, minXY, maxXY); fHMean.DrawAllGTracks(1, "ToF/HitXY/HitPointDist_" + cat, "hTofHitPointDist_" + cat, {""}, {""}, minD, maxD); } fHMean.DrawAllCands(2, "ToF/Time/HitTrackDist_cands", "hTofHitTrackDist_cands", {""}, {""}, 1e-9, 1.); fHMean.DrawAllCands(2, "ToF/Time/TimeVsMom_", "hTofTimeVsMom_cands", {""}, {""}, 1e-9, 1e-2); - fHMean.DrawAllGTracks(2, "ToF/Time/IdAndMatch", "hTofM2Calc_gTracks", {""}, {"true-ID", "mis-ID", "truematch", "mismatch"}, 1e-9, 10); + fHMean.DrawAllGTracks(2, "ToF/Time/IdAndMatch", "hTofM2Calc_gTracks", {""}, + {"true-ID", "mis-ID", "truematch", "mismatch"}, 1e-9, 10); TCanvas* c1 = fHMean.fHM.CreateCanvas("ToF/HitXY/HitPointDist_all", "ToF/HitXY/HitPointDist_all", 2400, 2400); c1->Divide(4, 4); int iC = 1; for (size_t iP = 0; iP < fHMean.fGTrackNames.size(); iP++) { c1->cd(iC++); - string hName = (iP <= 3) ? "hTofHitPointDist_trueid_" + fHMean.fGTrackNames[iP] : "hTofHitPointDist_misid_" + fHMean.fGTrackNames[iP]; + string hName = (iP <= 3) ? "hTofHitPointDist_trueid_" + fHMean.fGTrackNames[iP] + : "hTofHitPointDist_misid_" + fHMean.fGTrackNames[iP]; string text = (iP <= 3) ? "true-ID" : "mis-ID"; - TH1D* h = fHMean.H1Clone(hName); + TH1D* h = fHMean.H1Clone(hName); h->GetYaxis()->SetRangeUser(minD, maxD); h->Fit("gaus", "Q", "", 0., 2.); DrawH1(h, kLinear, kLog, "hist"); DrawTextOnPad(text, 0.5, 0.7, 0.8, 0.9); DrawTextOnPad(fHMean.fGTrackLatex[iP], 0.25, 0.9, 0.75, 0.99); - TF1* func = h->GetFunction("gaus"); + TF1* func = h->GetFunction("gaus"); double mean = (func != nullptr) ? func->GetParameter("Mean") : 0.; DrawTextOnPad("mean: " + Cbm::NumberToString(mean, 2), 0.2, 0.65, 0.6, 0.89); } @@ -1089,7 +1120,8 @@ void LmvmDrawAll::InvestigateMisid() double min = 1e-7; for (auto step : fHMean.fAnaSteps) { string cName = "ToF/Vertex/tofMisid_" + fHMean.fAnaStepNames[static_cast<int>(step)]; - vector<TH2D*> xyz {fHMean.H2("hVertexXZ_misidTof", step), fHMean.H2("hVertexYZ_misidTof", step), fHMean.H2("hVertexXY_misidTof", step)}; + vector<TH2D*> xyz {fHMean.H2("hVertexXZ_misidTof", step), fHMean.H2("hVertexYZ_misidTof", step), + fHMean.H2("hVertexXY_misidTof", step)}; TCanvas* c = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 1800, 600); c->Divide(3, 1); @@ -1137,13 +1169,18 @@ void LmvmDrawAll::InvestigateMisid() double minMism = 1e-4; double maxMism = 100.; - vector<string> xLabelDet = {"STS_{all}", "RICH_{all}", "RICH_{mis}", "TRD_{all}", "TRD_{mis}", "ToF_{all}", "ToF_{mis}"}; - fHMean.DrawAllGTracks(1, "Mismatches/hNofMismatches_all", "hNofMismatches_gTracks_all", xLabelDet, {""}, minMism, maxMism); - fHMean.DrawAllGTracks(1, "Mismatches/hNofMismatches_complete", "hNofMismatches_gTracks_complete", xLabelDet, {""}, minMism, maxMism); + vector<string> xLabelDet = {"STS_{all}", "RICH_{all}", "RICH_{mis}", "TRD_{all}", + "TRD_{mis}", "ToF_{all}", "ToF_{mis}"}; + fHMean.DrawAllGTracks(1, "Mismatches/hNofMismatches_all", "hNofMismatches_gTracks_all", xLabelDet, {""}, minMism, + maxMism); + fHMean.DrawAllGTracks(1, "Mismatches/hNofMismatches_complete", "hNofMismatches_gTracks_complete", xLabelDet, {""}, + minMism, maxMism); vector<string> xLabelSeg = {"0", "1", "2", "3"}; - fHMean.DrawAllGTracks(1, "Mismatches/hNofMismatchedTrackSegments_all", "hNofMismatchedTrackSegments_all", xLabelSeg, {""}, minMism, maxMism); - fHMean.DrawAllGTracks(1, "Mismatches/hNofMismatchedTrackSegments_complete", "hNofMismatchedTrackSegments_complete", xLabelSeg, {""}, minMism, maxMism); + fHMean.DrawAllGTracks(1, "Mismatches/hNofMismatchedTrackSegments_all", "hNofMismatchedTrackSegments_all", xLabelSeg, + {""}, minMism, maxMism); + fHMean.DrawAllGTracks(1, "Mismatches/hNofMismatchedTrackSegments_complete", "hNofMismatchedTrackSegments_complete", + xLabelSeg, {""}, minMism, maxMism); } void LmvmDrawAll::DrawChi2() @@ -1154,7 +1191,7 @@ void LmvmDrawAll::DrawChi2() fHMean.DrawAllCands(2, "hChi2/Detectors/hChi2VsMom_" + det, "hChi2VsMom_" + det, {""}, {""}, min, max); fHMean.DrawAllCands(2, "hChi2/Detectors/hTofTimeVsChi2_" + det, "hTofTimeVsChi2_" + det, {""}, {""}, min, max); } - + for (const string det : {"StsRich", "StsTrd", "RichTrd"}) { fHMean.DrawAllCands(2, "hChi2/Detectors/hChi2Comb" + det, "hChi2Comb_" + det, {""}, {""}, min, max); } @@ -1210,7 +1247,7 @@ void LmvmDrawAll::DrawMinv(ELmvmAnaStep step) "in-medium #rho"); drawData.emplace_back(cocktail, -1, kRed + 2, 4, -1, "Cocktail"); drawData.emplace_back(cb, -1, kCyan + 2, 4, -1, "CB"); - + double min = std::numeric_limits<Double_t>::max(); double max = std::numeric_limits<Double_t>::min(); TH1D* h0 = nullptr; @@ -1242,7 +1279,7 @@ void LmvmDrawAll::DrawMomRecoPrecision() { for (const string& cat : {"", "_zoom", "_zoom2"}) { string cName = "hMom/precisionReco/mom_ratioRecOverMc" + cat; - TCanvas* c = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 1800, 1800); + TCanvas* c = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 1800, 1800); c->Divide(3, 3); int i = 1; for (auto step : fHMean.fAnaSteps) { @@ -1250,7 +1287,7 @@ void LmvmDrawAll::DrawMomRecoPrecision() vector<TH1*> hists; vector<string> legend; c->cd(i++); - for (size_t iP = 0; iP < fHMean.fCandNames.size(); iP++ ) { + for (size_t iP = 0; iP < fHMean.fCandNames.size(); iP++) { TH1D* h = fHMean.H1Clone("hMomRatio_cands_" + fHMean.fCandNames[iP], step); if (cat == "_zoom") h->GetXaxis()->SetRangeUser(0.9, 1.1); if (cat == "_zoom2") h->GetXaxis()->SetRangeUser(7., 7.5); @@ -1266,19 +1303,20 @@ void LmvmDrawAll::DrawMomRecoPrecision() // draw ratio P_rec/P_MC vs P { ELmvmAnaStep step = ELmvmAnaStep::ElId; - double min = 0.5; - double max = 1.5; - TCanvas* c = fHMean.fHM.CreateCanvas("hMom/precisionReco/mom_ratioRecOverMcVsMom", "hMom/precisionReco/mom_ratioRecOverMcVsMom", 2700, 1800); - c->Divide(4,3); + double min = 0.5; + double max = 1.5; + TCanvas* c = fHMean.fHM.CreateCanvas("hMom/precisionReco/mom_ratioRecOverMcVsMom", + "hMom/precisionReco/mom_ratioRecOverMcVsMom", 2700, 1800); + c->Divide(4, 3); int i = 1; for (auto ptcl : fHMean.fCandNames) { - int iL = i-1; + int iL = i - 1; c->cd(i++); string hName = "hMomRatioVsMom_cands_" + ptcl; - TH2D* h = fHMean.H2Clone(hName.c_str(), step); + TH2D* h = fHMean.H2Clone(hName.c_str(), step); h->GetYaxis()->SetRangeUser(min, max); DrawH2(h, kLinear, kLinear, kLog, "colz"); - string text = fHMean.fCandLatex[iL] + ", "+ fHMean.fAnaStepLatex[static_cast<int>(step)]; + string text = fHMean.fCandLatex[iL] + ", " + fHMean.fAnaStepLatex[static_cast<int>(step)]; DrawTextOnPad(text.c_str(), 0.25, 0.9, 0.75, 0.999); } } @@ -1287,24 +1325,24 @@ void LmvmDrawAll::DrawMomRecoPrecision() void LmvmDrawAll::DrawMinvOfficialStyle() { // variable binning - TH1D* h = fHMean.H1("hMinv", ELmvmSrc::Bg, ELmvmAnaStep::Mc); // template for values - double ChangeValue = 1.2; - int f = 2; // factor larger/smaller binsize - int binChange = h->GetXaxis()->FindBin(ChangeValue); - double bW = h->GetXaxis()->GetBinWidth(1); - const int nBins = (int) binChange + (h->GetNbinsX() - binChange)/f; - vector<double> binEdges(nBins+1, 0.); + TH1D* h = fHMean.H1("hMinv", ELmvmSrc::Bg, ELmvmAnaStep::Mc); // template for values + double ChangeValue = 1.2; + int f = 2; // factor larger/smaller binsize + int binChange = h->GetXaxis()->FindBin(ChangeValue); + double bW = h->GetXaxis()->GetBinWidth(1); + const int nBins = (int) binChange + (h->GetNbinsX() - binChange) / f; + vector<double> binEdges(nBins + 1, 0.); for (int iB = 1; iB <= nBins; iB++) { - binEdges[iB] = (iB < binChange) ? binEdges[iB-1] + bW : binEdges[iB-1] + f * bW; + binEdges[iB] = (iB < binChange) ? binEdges[iB - 1] + bW : binEdges[iB - 1] + f * bW; } - binEdges[nBins+1] = h->GetBinCenter(h->GetNbinsX()) + bW/f; + binEdges[nBins + 1] = h->GetBinCenter(h->GetNbinsX()) + bW / f; - int nofInmed = H(ELmvmSignal::Inmed) ->H1("hEventNumber")->GetEntries(); - int nofQgp = H(ELmvmSignal::Qgp) ->H1("hEventNumber")->GetEntries(); - int nofOmega = H(ELmvmSignal::Omega) ->H1("hEventNumber")->GetEntries(); + int nofInmed = H(ELmvmSignal::Inmed)->H1("hEventNumber")->GetEntries(); + int nofQgp = H(ELmvmSignal::Qgp)->H1("hEventNumber")->GetEntries(); + int nofOmega = H(ELmvmSignal::Omega)->H1("hEventNumber")->GetEntries(); int nofOmegaD = H(ELmvmSignal::OmegaD)->H1("hEventNumber")->GetEntries(); - int nofPhi = H(ELmvmSignal::Phi) ->H1("hEventNumber")->GetEntries(); + int nofPhi = H(ELmvmSignal::Phi)->H1("hEventNumber")->GetEntries(); for (auto step : fHMean.fAnaSteps) { if (step < ELmvmAnaStep::ElId) continue; @@ -1339,39 +1377,42 @@ void LmvmDrawAll::DrawMinvOfficialStyle() for (int iB = 1; iB <= h->GetNbinsX(); iB++) { if (iB < binChange) { - h_npm ->SetBinContent(iB, fHMean.H1("hMinvCombPM_sameEv", step)->GetBinContent(iB)); - h_bc ->SetBinContent(iB, fHMean.H1("hMinvCombBg", step)->GetBinContent(iB)); - h_coc ->SetBinContent(iB, GetCocktailMinvH1("hMinv", step)->GetBinContent(iB)); - h_eta ->SetBinContent(iB, fHMean.H1("hMinv", ELmvmSrc::Eta, step)->GetBinContent(iB)); - h_pi0 ->SetBinContent(iB, fHMean.H1("hMinv", ELmvmSrc::Pi0, step)->GetBinContent(iB)); - h_inmed ->SetBinContent(iB, hInmed0 ->GetBinContent(iB)); - h_qgp ->SetBinContent(iB, hQgp0 ->GetBinContent(iB)); - h_omega ->SetBinContent(iB, hOmega0 ->GetBinContent(iB)); + h_npm->SetBinContent(iB, fHMean.H1("hMinvCombPM_sameEv", step)->GetBinContent(iB)); + h_bc->SetBinContent(iB, fHMean.H1("hMinvCombBg", step)->GetBinContent(iB)); + h_coc->SetBinContent(iB, GetCocktailMinvH1("hMinv", step)->GetBinContent(iB)); + h_eta->SetBinContent(iB, fHMean.H1("hMinv", ELmvmSrc::Eta, step)->GetBinContent(iB)); + h_pi0->SetBinContent(iB, fHMean.H1("hMinv", ELmvmSrc::Pi0, step)->GetBinContent(iB)); + h_inmed->SetBinContent(iB, hInmed0->GetBinContent(iB)); + h_qgp->SetBinContent(iB, hQgp0->GetBinContent(iB)); + h_omega->SetBinContent(iB, hOmega0->GetBinContent(iB)); h_omegaD->SetBinContent(iB, hOmegaD0->GetBinContent(iB)); - h_phi ->SetBinContent(iB, hPhi0 ->GetBinContent(iB)); + h_phi->SetBinContent(iB, hPhi0->GetBinContent(iB)); } else { - int iB2 = (int) binChange + (iB-binChange)/f; - h_npm ->SetBinContent(iB2, h_npm->GetBinContent(iB2) + fHMean.H1("hMinvCombPM_sameEv", step)->GetBinContent(iB)/f); - h_bc ->SetBinContent(iB2, h_bc->GetBinContent(iB2) + fHMean.H1("hMinvCombBg", step)->GetBinContent(iB)/f); - h_coc ->SetBinContent(iB2, h_coc->GetBinContent(iB2) + GetCocktailMinvH1("hMinv", step)->GetBinContent(iB)/f); - h_eta ->SetBinContent(iB2, h_eta->GetBinContent(iB2) + fHMean.H1("hMinv", ELmvmSrc::Eta, step)->GetBinContent(iB)/f); - h_pi0 ->SetBinContent(iB2, h_pi0->GetBinContent(iB2) + fHMean.H1("hMinv", ELmvmSrc::Pi0, step)->GetBinContent(iB)/f); - h_inmed ->SetBinContent(iB2, h_inmed->GetBinContent(iB2) + hInmed0 ->GetBinContent(iB)/f); - h_qgp ->SetBinContent(iB2, h_qgp->GetBinContent(iB2) + hQgp0 ->GetBinContent(iB)/f); - h_omega ->SetBinContent(iB2, h_omega->GetBinContent(iB2) + hOmega0 ->GetBinContent(iB)/f); - h_omegaD->SetBinContent(iB2, h_omegaD->GetBinContent(iB2) + hOmegaD0->GetBinContent(iB)/f); - h_phi ->SetBinContent(iB2, h_phi->GetBinContent(iB2) + hPhi0 ->GetBinContent(iB)/f); + int iB2 = (int) binChange + (iB - binChange) / f; + h_npm->SetBinContent(iB2, + h_npm->GetBinContent(iB2) + fHMean.H1("hMinvCombPM_sameEv", step)->GetBinContent(iB) / f); + h_bc->SetBinContent(iB2, h_bc->GetBinContent(iB2) + fHMean.H1("hMinvCombBg", step)->GetBinContent(iB) / f); + h_coc->SetBinContent(iB2, h_coc->GetBinContent(iB2) + GetCocktailMinvH1("hMinv", step)->GetBinContent(iB) / f); + h_eta->SetBinContent(iB2, h_eta->GetBinContent(iB2) + + fHMean.H1("hMinv", ELmvmSrc::Eta, step)->GetBinContent(iB) / f); + h_pi0->SetBinContent(iB2, h_pi0->GetBinContent(iB2) + + fHMean.H1("hMinv", ELmvmSrc::Pi0, step)->GetBinContent(iB) / f); + h_inmed->SetBinContent(iB2, h_inmed->GetBinContent(iB2) + hInmed0->GetBinContent(iB) / f); + h_qgp->SetBinContent(iB2, h_qgp->GetBinContent(iB2) + hQgp0->GetBinContent(iB) / f); + h_omega->SetBinContent(iB2, h_omega->GetBinContent(iB2) + hOmega0->GetBinContent(iB) / f); + h_omegaD->SetBinContent(iB2, h_omegaD->GetBinContent(iB2) + hOmegaD0->GetBinContent(iB) / f); + h_phi->SetBinContent(iB2, h_phi->GetBinContent(iB2) + hPhi0->GetBinContent(iB) / f); } } // set errors int nEv = GetNofTotalEvents(); for (int iB = 1; iB <= h_npm->GetNbinsX(); iB++) { - double bW2 = h_npm->GetBinLowEdge(iB+1) - h_npm->GetBinLowEdge(iB); - //cout << "iB = " << iB << " | bW2 = " << bW2 << endl; + double bW2 = h_npm->GetBinLowEdge(iB + 1) - h_npm->GetBinLowEdge(iB); + //cout << "iB = " << iB << " | bW2 = " << bW2 << endl; h_npm->SetBinError(iB, std::sqrt(h_npm->GetBinContent(iB) * bW2 * nEv) / (bW2 * nEv)); - h_bc ->SetBinError(iB, std::sqrt(h_bc ->GetBinContent(iB) * bW2 * nEv) / (bW2 * nEv)); + h_bc->SetBinError(iB, std::sqrt(h_bc->GetBinContent(iB) * bW2 * nEv) / (bW2 * nEv)); } TH1D* h_sum = (TH1D*) h_eta->Clone(); @@ -1384,21 +1425,24 @@ void LmvmDrawAll::DrawMinvOfficialStyle() TH1D* h_sum2 = (TH1D*) h_sum->Clone(); - fHMean.SetOptH1(h_npm, "#font[52]{M}_{ee} [GeV/#font[52]{c}^{2}]", "1/N_{ev} dN/d#font[52]{M}_{ee} [GeV/#font[52]{c}^{2}]^{-1}", 510, 21, 1.3, kPink, "marker"); - fHMean.SetOptH1(h_bc, "#font[52]{M}_{ee} [GeV/#font[52]{c}^{2}]", "1/N_{ev} dN/d#font[52]{M}_{ee} [GeV/#font[52]{c}^{2}]^{-1}", 510, 22, 1.3, kBlue-2, "marker"); - fHMean.SetOptH1(h_sum, "#font[52]{M}_{ee} [GeV/#font[52]{c}^{2}]", "1/N_{ev} dN/d#font[52]{M}_{ee} [GeV/#font[52]{c}^{2}]^{-1}", 510, 20, 1.1, kBlack, "marker"); - - fHMean.SetOptH1(h_sum2, "", "",510, 1, 1, kBlack, "line"); - fHMean.SetOptH1(h_eta, "", "",510, 1, 1, kBlue+2, "line"); - fHMean.SetOptH1(h_pi0, "", "",510, 1, 1, kAzure-4, "line"); - fHMean.SetOptH1(h_inmed, "", "",510, 1, 1, kPink, "line"); - fHMean.SetOptH1(h_qgp, "", "",510, 1, 1, kOrange, "line"); - fHMean.SetOptH1(h_omega, "", "",510, 1, 1, kGreen+1, "line"); - fHMean.SetOptH1(h_omegaD, "", "",510, 1, 1, kGreen+3, "line"); - fHMean.SetOptH1(h_phi, "", "",510, 1, 1, kMagenta+2, "line"); + fHMean.SetOptH1(h_npm, "#font[52]{M}_{ee} [GeV/#font[52]{c}^{2}]", + "1/N_{ev} dN/d#font[52]{M}_{ee} [GeV/#font[52]{c}^{2}]^{-1}", 510, 21, 1.3, kPink, "marker"); + fHMean.SetOptH1(h_bc, "#font[52]{M}_{ee} [GeV/#font[52]{c}^{2}]", + "1/N_{ev} dN/d#font[52]{M}_{ee} [GeV/#font[52]{c}^{2}]^{-1}", 510, 22, 1.3, kBlue - 2, "marker"); + fHMean.SetOptH1(h_sum, "#font[52]{M}_{ee} [GeV/#font[52]{c}^{2}]", + "1/N_{ev} dN/d#font[52]{M}_{ee} [GeV/#font[52]{c}^{2}]^{-1}", 510, 20, 1.1, kBlack, "marker"); + + fHMean.SetOptH1(h_sum2, "", "", 510, 1, 1, kBlack, "line"); + fHMean.SetOptH1(h_eta, "", "", 510, 1, 1, kBlue + 2, "line"); + fHMean.SetOptH1(h_pi0, "", "", 510, 1, 1, kAzure - 4, "line"); + fHMean.SetOptH1(h_inmed, "", "", 510, 1, 1, kPink, "line"); + fHMean.SetOptH1(h_qgp, "", "", 510, 1, 1, kOrange, "line"); + fHMean.SetOptH1(h_omega, "", "", 510, 1, 1, kGreen + 1, "line"); + fHMean.SetOptH1(h_omegaD, "", "", 510, 1, 1, kGreen + 3, "line"); + fHMean.SetOptH1(h_phi, "", "", 510, 1, 1, kMagenta + 2, "line"); //h_npm->GetXaxis()->SetRangeUser(0.,2.); - h_npm->GetYaxis()->SetRangeUser(1e-9,10); + h_npm->GetYaxis()->SetRangeUser(1e-9, 10); string cName = "MinvOfficialStyle/minv_" + fHMean.fAnaStepNames[static_cast<int>(step)]; TCanvas* can = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 900, 918); @@ -1407,43 +1451,43 @@ void LmvmDrawAll::DrawMinvOfficialStyle() can->cd(); h_npm->Draw("peist"); - h_bc ->Draw("peistsame"); + h_bc->Draw("peistsame"); h_sum->Draw("phistsame"); - h_sum2 ->Draw("histsame"); - h_eta ->Draw("histsame"); - h_pi0 ->Draw("histsame"); - h_inmed ->Draw("histsame"); - h_qgp ->Draw("histsame"); - h_omega ->Draw("histsame"); + h_sum2->Draw("histsame"); + h_eta->Draw("histsame"); + h_pi0->Draw("histsame"); + h_inmed->Draw("histsame"); + h_qgp->Draw("histsame"); + h_omega->Draw("histsame"); h_omegaD->Draw("histsame"); - h_phi ->Draw("histsame"); + h_phi->Draw("histsame"); vector<LmvmLegend> legend3; legend3.emplace_back(h_npm, "N^{#pm}_{same}", "p"); - legend3.emplace_back(h_bc, "CB_{calc}", "p"); - legend3.emplace_back(h_sum, "MC Cocktail", "p"); - + legend3.emplace_back(h_bc, "CB_{calc}", "p"); + legend3.emplace_back(h_sum, "MC Cocktail", "p"); + vector<LmvmLegend> legend8; - legend8.emplace_back(h_sum2, "MC Cocktail", "l"); - legend8.emplace_back(h_eta, "#eta#rightarrow#gammae^{+}e^{-}", "l"); - legend8.emplace_back(h_pi0, "#pi^{0}#rightarrow#gammae^{+}e^{-}", "l"); + legend8.emplace_back(h_sum2, "MC Cocktail", "l"); + legend8.emplace_back(h_eta, "#eta#rightarrow#gammae^{+}e^{-}", "l"); + legend8.emplace_back(h_pi0, "#pi^{0}#rightarrow#gammae^{+}e^{-}", "l"); legend8.emplace_back(h_omegaD, "#omega#rightarrow#pi^{0}e^{+}e^{-}", "l"); - legend8.emplace_back(h_omega, "#omega#rightarrowe^{+}e^{-}", "l"); - legend8.emplace_back(h_phi, "#phi#rightarrowe^{+}e^{-}", "l"); - legend8.emplace_back(h_inmed, "Rapp in-medium SF", "l"); - legend8.emplace_back(h_qgp, "Rapp QGP", "l"); - + legend8.emplace_back(h_omega, "#omega#rightarrowe^{+}e^{-}", "l"); + legend8.emplace_back(h_phi, "#phi#rightarrowe^{+}e^{-}", "l"); + legend8.emplace_back(h_inmed, "Rapp in-medium SF", "l"); + legend8.emplace_back(h_qgp, "Rapp QGP", "l"); + fHMean.SetLegend(legend3, 0.035, 0.37, 0.64, 0.57, 0.77); fHMean.SetLegend(legend8, 0.025, 0.66, 0.57, 0.82, 0.87); - TLatex* tex = new TLatex(0.15, 2.5,"CBM Simulations"); - tex->SetTextColor(1); //kAzure+6); + TLatex* tex = new TLatex(0.15, 2.5, "CBM Simulations"); + tex->SetTextColor(1); //kAzure+6); tex->SetTextSize(0.035); tex->SetTextFont(42); tex->Draw(); - TLatex* tex2 = new TLatex(0.15, 0.7,"0 fm Au+Au, #sqrt{s_{NN}}=4.1 GeV"); - tex2->SetTextColor(1); //kAzure+6); + TLatex* tex2 = new TLatex(0.15, 0.7, "0 fm Au+Au, #sqrt{s_{NN}}=4.1 GeV"); + tex2->SetTextColor(1); //kAzure+6); tex2->SetTextSize(0.035); tex2->SetTextFont(42); tex2->Draw(); @@ -1460,7 +1504,7 @@ void LmvmDrawAll::DrawMinvOfficialStyle() delete h_phi; delete h_sum; delete h_sum2;*/ - } // steps + } // steps } void LmvmDrawAll::DrawMinvBgSourcesAll() @@ -1563,7 +1607,7 @@ void LmvmDrawAll::DrawMinvPtAll() void LmvmDrawAll::DrawMinvCombBgAndSignal() { - string folder = "CombinatorialBackground/"; + string folder = "CombinatorialBackground/"; string folderUAll = folder + "UrQMD/All/"; string folderUEl = folder + "UrQMD/Electrons/"; @@ -1571,14 +1615,16 @@ void LmvmDrawAll::DrawMinvCombBgAndSignal() { for (const string& cat : {"", "_urqmdEl", "_urqmdAll"}) { for (const string& ev : {"sameEv", "mixedEv"}) { - string cName = (cat == "_urqmdAll") ? folderUAll + "PairYields_" + ev :(cat == "_urqmdEl") ? folderUEl + "PairYields_" + ev : folder + "PairYields_" + ev; - TCanvas* c = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 1800, 1800); + string cName = (cat == "_urqmdAll") ? folderUAll + "PairYields_" + ev + : (cat == "_urqmdEl") ? folderUEl + "PairYields_" + ev + : folder + "PairYields_" + ev; + TCanvas* c = fHMean.fHM.CreateCanvas(cName.c_str(), cName.c_str(), 1800, 1800); c->Divide(3, 3); int i = 1; for (auto step : fHMean.fAnaSteps) { if (step < ELmvmAnaStep::ElId) continue; c->cd(i++); - TH1D* pp = fHMean.H1Clone("hMinvCombPP" + cat + "_"+ ev, step); + TH1D* pp = fHMean.H1Clone("hMinvCombPP" + cat + "_" + ev, step); TH1D* mm = fHMean.H1Clone("hMinvCombMM" + cat + "_" + ev, step); TH1D* pm = fHMean.H1Clone("hMinvCombPM" + cat + "_" + ev, step); pm->GetYaxis()->SetTitle("Yield"); @@ -1594,7 +1640,7 @@ void LmvmDrawAll::DrawMinvCombBgAndSignal() { for (const string& cat : {"", "_urqmdEl"}) { string cName = (cat == "_urqmdEl") ? folderUEl : folder; - TCanvas* c = fHMean.fHM.CreateCanvas(cName + "RatioMMPP_same", cName + "RatioMMPP_same", 1800, 1800); + TCanvas* c = fHMean.fHM.CreateCanvas(cName + "RatioMMPP_same", cName + "RatioMMPP_same", 1800, 1800); c->Divide(3, 3); int i = 1; for (auto step : fHMean.fAnaSteps) { @@ -1614,7 +1660,7 @@ void LmvmDrawAll::DrawMinvCombBgAndSignal() { for (const string& cat : {"", "_urqmdEl", "_urqmdAll"}) { string cName = (cat == "_urqmdAll") ? folderUAll : (cat == "_urqmdEl") ? folderUEl : folder; - TCanvas* c = fHMean.fHM.CreateCanvas(cName + "geomMean", cName + "geomMean", 1800, 1800); + TCanvas* c = fHMean.fHM.CreateCanvas(cName + "geomMean", cName + "geomMean", 1800, 1800); c->Divide(3, 3); int i = 1; for (auto step : fHMean.fAnaSteps) { @@ -1629,20 +1675,22 @@ void LmvmDrawAll::DrawMinvCombBgAndSignal() same->GetXaxis()->SetTitle("M_{ee} [GeV/c^{2}]"); same->SetMinimum(1e-8); mixed->SetMinimum(1e-8); - DrawH1({same, mixed}, {"geom. mean (same)", "geom. mean (mixed)"}, kLinear, kLog, true, 0.52, 0.8, 0.99, 0.99, "p"); + DrawH1({same, mixed}, {"geom. mean (same)", "geom. mean (mixed)"}, kLinear, kLog, true, 0.52, 0.8, 0.99, 0.99, + "p"); fHMean.DrawAnaStepOnPad(step); } } } // Draw k factor - { // draw all three k factors in one histo // TODO: keep this or next method or both? + { // draw all three k factors in one histo // TODO: keep this or next method or both? ELmvmAnaStep step = ELmvmAnaStep::ElId; fHMean.fHM.CreateCanvas(folder + "k_all", folder + "k_all", 1800, 1800); TH1D* kAll = fHMean.H1Clone("hMinvCombK", step); TH1D* kUAll = fHMean.H1Clone("hMinvCombK_urqmdAll", step); TH1D* kUEl = fHMean.H1Clone("hMinvCombK_urqmdEl", step); - DrawH1({kAll, kUAll, kUEl}, {"k (all)", "k (UrQMD all)", "k (UrQMD El"}, kLinear, kLinear, true, 0.65, 0.8, 0.99, 0.95, "p"); + DrawH1({kAll, kUAll, kUEl}, {"k (all)", "k (UrQMD all)", "k (UrQMD El"}, kLinear, kLinear, true, 0.65, 0.8, 0.99, + 0.95, "p"); kAll->GetYaxis()->SetTitle("k Factor"); kAll->GetYaxis()->SetRangeUser(0.8, 1.2); fHMean.DrawAnaStepOnPad(step); @@ -1654,13 +1702,13 @@ void LmvmDrawAll::DrawMinvCombBgAndSignal() for (const string& cat : {"", "_urqmdEl", "_urqmdAll"}) { string cName = (cat == "_urqmdAll") ? folderUAll : (cat == "_urqmdEl") ? folderUEl : folder; - TCanvas* c = fHMean.fHM.CreateCanvas(cName + "N+-VsBgCoc", cName + "N+-VsBgCoc", 1800, 900); + TCanvas* c = fHMean.fHM.CreateCanvas(cName + "N+-VsBgCoc", cName + "N+-VsBgCoc", 1800, 900); c->Divide(2, 1); - TH1D* npm = fHMean.H1Clone("hMinvCombPM" + cat + "_sameEv", step); - TH1D* cb = fHMean.H1Clone("hMinvCombBg" + cat, step); - TH1D* bg = fHMean.H1Clone("hMinv" + cat, ELmvmSrc::Bg, step); + TH1D* npm = fHMean.H1Clone("hMinvCombPM" + cat + "_sameEv", step); + TH1D* cb = fHMean.H1Clone("hMinvCombBg" + cat, step); + TH1D* bg = fHMean.H1Clone("hMinv" + cat, ELmvmSrc::Bg, step); TH1D* cocktail = GetCocktailMinvH1("hMinv" + cat, step); - TH1D* sbg = static_cast<TH1D*>(bg->Clone()); + TH1D* sbg = static_cast<TH1D*>(bg->Clone()); sbg->Add(cocktail); TH1D* ratS = (TH1D*) npm->Clone(); ratS->Divide(sbg); @@ -1668,13 +1716,15 @@ void LmvmDrawAll::DrawMinvCombBgAndSignal() ratB->Divide(bg); c->cd(1); npm->GetYaxis()->SetRangeUser(1e-6, 2e-2); - DrawH1({npm, sbg, cb, bg}, {"N^{+-}_{same}", "BG + Cocktail", "CB_{calc}", "Background"}, kLinear, kLog, true, 0.7, 0.7, 0.99, 0.95, "p"); + DrawH1({npm, sbg, cb, bg}, {"N^{+-}_{same}", "BG + Cocktail", "CB_{calc}", "Background"}, kLinear, kLog, true, + 0.7, 0.7, 0.99, 0.95, "p"); fHMean.DrawAnaStepOnPad(step); c->cd(2); ratS->GetYaxis()->SetTitle("Ratio"); ratS->GetYaxis()->SetRangeUser(0., 2.); ratB->GetYaxis()->SetRangeUser(0., 2.); - DrawH1({ratS, ratB}, {"N^{+-}_{same} / BG+Coc", "CB_{calc} / BG"}, kLinear, kLinear, true, 0.7, 0.8, 0.99, 0.95, "p"); + DrawH1({ratS, ratB}, {"N^{+-}_{same} / BG+Coc", "CB_{calc} / BG"}, kLinear, kLinear, true, 0.7, 0.8, 0.99, 0.95, + "p"); } } @@ -1685,14 +1735,14 @@ void LmvmDrawAll::DrawMinvCombBgAndSignal() TH1D* rsAll = fHMean.H1Clone("hMinvCombPM_sameEv", step); TH1D* rsUAll = fHMean.H1Clone("hMinvCombPM_urqmdAll_sameEv", step); TH1D* rsUEl = fHMean.H1Clone("hMinvCombPM_urqmdEl_sameEv", step); - + TH1D* cbgAll = GetCocktailMinvH1("hMinv", step); cbgAll->Add(fHMean.H1("hMinv", ELmvmSrc::Bg, step)); TH1D* cbgUAll = GetCocktailMinvH1("hMinv_urqmdAll", step); cbgUAll->Add(fHMean.H1("hMinv_urqmdAll", ELmvmSrc::Bg, step)); TH1D* cbgUEl = GetCocktailMinvH1("hMinv_urqmdEl", step); cbgUEl->Add(fHMean.H1("hMinv_urqmdEl", ELmvmSrc::Bg, step)); - + rsAll->Divide(cbgAll); rsUAll->Divide(cbgUAll); rsUEl->Divide(cbgUEl); @@ -1708,15 +1758,17 @@ void LmvmDrawAll::DrawMinvCombBgAndSignal() rsAll->GetYaxis()->SetTitle("N^{+-}_{same} / BG+Coc"); rsAll->GetYaxis()->SetRangeUser(0., 2.); rbAll->GetYaxis()->SetTitle("CB_{calc} / BG"); - rbAll->GetYaxis()->SetRangeUser(0., 2.); + rbAll->GetYaxis()->SetRangeUser(0., 2.); TCanvas* c = fHMean.fHM.CreateCanvas(folder + "CbVsMc_ratio", folder + "CbVsMc_ratio", 1800, 900); c->Divide(2, 1); c->cd(1); - DrawH1({rsAll, rsUAll, rsUEl}, {"UrQMD + PLUTO", "UrQMD", "UrQMD electrons"}, kLinear, kLinear, true, 0.55, 0.78, 0.99, 0.95, "hist"); + DrawH1({rsAll, rsUAll, rsUEl}, {"UrQMD + PLUTO", "UrQMD", "UrQMD electrons"}, kLinear, kLinear, true, 0.55, 0.78, + 0.99, 0.95, "hist"); fHMean.DrawAnaStepOnPad(step); c->cd(2); - DrawH1({rbAll, rbUAll, rbUEl}, {"UrQMD + PLUTO", "UrQMD", "UrQMD electrons"}, kLinear, kLinear, true, 0.55, 0.78, 0.99, 0.95, "hist"); + DrawH1({rbAll, rbUAll, rbUEl}, {"UrQMD + PLUTO", "UrQMD", "UrQMD electrons"}, kLinear, kLinear, true, 0.55, 0.78, + 0.99, 0.95, "hist"); fHMean.DrawAnaStepOnPad(step); } @@ -1739,7 +1791,7 @@ void LmvmDrawAll::DrawMinvCombBgAndSignal() { for (const string& cat : {"", "_urqmdEl", "_urqmdAll"}) { string cName = (cat == "_urqmdEl") ? folderUEl : (cat == "_urqmdAll") ? folderUAll : folder; - TCanvas* c = fHMean.fHM.CreateCanvas(cName + "CbVsInput_Npm", cName + "CbVsInput_Npm", 1800, 1800); + TCanvas* c = fHMean.fHM.CreateCanvas(cName + "CbVsInput_Npm", cName + "CbVsInput_Npm", 1800, 1800); c->Divide(3, 3); int i = 1; for (auto step : fHMean.fAnaSteps) { @@ -1752,8 +1804,8 @@ void LmvmDrawAll::DrawMinvCombBgAndSignal() pmSame->SetMaximum(1e-2); pmSame->SetMinimum(4e-9); DrawH1({pmSame, combBg, combSignalNpm, cocktail}, - {"N_{same}^{+-}", "B_{c}", "Signal (N_{same}^{+-} - B_{c})", "Cocktail"}, kLinear, kLog, true, 0.53, 0.75, - 0.99, 0.92, "p"); + {"N_{same}^{+-}", "B_{c}", "Signal (N_{same}^{+-} - B_{c})", "Cocktail"}, kLinear, kLog, true, 0.53, + 0.75, 0.99, 0.92, "p"); fHMean.DrawAnaStepOnPad(step); } } @@ -1763,7 +1815,7 @@ void LmvmDrawAll::DrawMinvCombBgAndSignal() { for (const string& cat : {"", "_urqmdEl", "_urqmdAll"}) { string cName = (cat == "_urqmdEl") ? folderUEl : (cat == "_urqmdAll") ? folderUAll : folder; - TCanvas* c = fHMean.fHM.CreateCanvas(cName + "CbVsInput_cocBg", cName + "CbVsInput_cocBg", 1800, 1800); + TCanvas* c = fHMean.fHM.CreateCanvas(cName + "CbVsInput_cocBg", cName + "CbVsInput_cocBg", 1800, 1800); c->Divide(3, 3); int i = 1; for (auto step : fHMean.fAnaSteps) { @@ -1775,9 +1827,9 @@ void LmvmDrawAll::DrawMinvCombBgAndSignal() TH1D* sBgSignal = fHMean.H1Clone("hMinvCombSignalMc" + cat, step); sbg2->SetMaximum(2e-2); sbg2->SetMinimum(4e-9); - TH1D* cocktail = GetCocktailMinvH1("hMinv" + cat, step); + TH1D* cocktail = GetCocktailMinvH1("hMinv" + cat, step); DrawH1({sbg2, combBg, sBgSignal, cocktail}, {"Cocktail + BG", "B_{C}", "Signal (Cock+BG - B_{C})", "Cocktail"}, - kLinear, kLog, true, 0.53, 0.72, 0.99, 0.92, "p"); + kLinear, kLog, true, 0.53, 0.72, 0.99, 0.92, "p"); fHMean.DrawAnaStepOnPad(step); } } diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.h b/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.h index ede49439e5b27848457b33737031cf41fa792066..1ec2314b6eeedc07515da5342b74007f7b04c5d1 100755 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.h +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmDrawAll.h @@ -69,7 +69,7 @@ private: /** * \brief Draw invariant mass spectra in official style. */ - void DrawMinvOfficialStyle(); + void DrawMinvOfficialStyle(); template<class T> void CreateMeanHist(const std::string& name, int nofEvents, int nofRebins = -1); // if nRebin = -1, no rebin @@ -118,7 +118,8 @@ private: void DrawGTrackVertex(); void DrawSignificancesAll(); - void DrawSignificance(TH2D* hEl, TH2D* hBg, const std::string& name, double minZ, double maxZ, const std::string& option); + void DrawSignificance(TH2D* hEl, TH2D* hBg, const std::string& name, double minZ, double maxZ, + const std::string& option); void DrawCutEffSignal(); @@ -148,7 +149,7 @@ private: */ void SaveCanvasToImage(); - double fZ = -44.; // z-position of target + double fZ = -44.; // z-position of target LmvmDrawAll(const LmvmDrawAll&); LmvmDrawAll operator=(const LmvmDrawAll&); diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmHist.cxx b/analysis/PWGDIL/dielectron/lmvm/LmvmHist.cxx index 957d44f9273751e198700ee828d033c99a1ba5f1..1621ee8dd25693d662676f0dd9510e054c765093 100755 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmHist.cxx +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmHist.cxx @@ -10,9 +10,9 @@ #include "Logger.h" -#include "TText.h" -#include "TStyle.h" #include "TLegend.h" +#include "TStyle.h" +#include "TText.h" #include "LmvmUtils.h" @@ -48,14 +48,24 @@ const vector<std::string> LmvmHist::fBgPairSrcLatex = {"#gamma-#gamma", "#pi^{0 // the following candidates and global track vectors are mainly set up to investigate misidentifications const vector<std::string> LmvmHist::fGTrackNames = { - "el+", "el+_prim", "el-", "el-_prim", "pion+", "pion+_prim", "pion-", "pion-_prim", - "proton", "proton_prim", "kaon+", "kaon+_prim", "kaon-", "kaon-_prim", "other"}; + "el+", "el+_prim", "el-", "el-_prim", "pion+", "pion+_prim", "pion-", "pion-_prim", + "proton", "proton_prim", "kaon+", "kaon+_prim", "kaon-", "kaon-_prim", "other"}; const vector<std::string> LmvmHist::fGTrackLatex = { - "e^{+}", "e^{+}_{prim}", "e^{-}", "e^{-}_{prim}", "#pi^{+}", "#pi^{+}_{prim}", "#pi^{-}", "#pi^{-}_{prim}", - "p", "p_{prim}", "K^{+}", "K^{+}_{prim}", "K^{-}", "K^{-}_{prim}", "other"}; - -const vector<std::string> LmvmHist::fCandNames = {"plutoEl+", "plutoEl-", "urqmdEl+", "urqmdEl-", "pion+", "pion-", "proton", "kaon+", "kaon-", "other"}; -const vector<std::string> LmvmHist::fCandLatex = {"e^{+}_{PLUTO}", "e^{-}_{PLUTO}", "e^{+}_{UrQMD}", "e^{-}_{UrQMD}", "#pi^{+}", "#pi^{-}", "p", "K^{+}", "K^{-}", "o."}; + "e^{+}", "e^{+}_{prim}", "e^{-}", "e^{-}_{prim}", "#pi^{+}", "#pi^{+}_{prim}", "#pi^{-}", "#pi^{-}_{prim}", + "p", "p_{prim}", "K^{+}", "K^{+}_{prim}", "K^{-}", "K^{-}_{prim}", "other"}; + +const vector<std::string> LmvmHist::fCandNames = {"plutoEl+", "plutoEl-", "urqmdEl+", "urqmdEl-", "pion+", + "pion-", "proton", "kaon+", "kaon-", "other"}; +const vector<std::string> LmvmHist::fCandLatex = {"e^{+}_{PLUTO}", + "e^{-}_{PLUTO}", + "e^{+}_{UrQMD}", + "e^{-}_{UrQMD}", + "#pi^{+}", + "#pi^{-}", + "p", + "K^{+}", + "K^{-}", + "o."}; LmvmHist::LmvmHist() {} @@ -154,7 +164,7 @@ void LmvmHist::FillH2(const string& name, ELmvmSrc src, double x, double y, doub { if (src == ELmvmSrc::Undefined) return; //double myWeight = (src == ELmvmSrc::Signal) ? wSignal : 1.; // TODO: delete commented lines - //H2(name, src)->Fill(x, y, myWeight); + //H2(name, src)->Fill(x, y, myWeight); H2(name, src)->Fill(x, y, wSignal); } @@ -169,7 +179,7 @@ void LmvmHist::FillH1(const string& name, ELmvmSrc src, ELmvmAnaStep step, doubl { if (src == ELmvmSrc::Undefined || step == ELmvmAnaStep::Undefined) return; //double myWeight = (src == ELmvmSrc::Signal) ? wSignal : 1.; // TODO: delete commented lines - //FillH1(GetName(name, src, step), x, myWeight); + //FillH1(GetName(name, src, step), x, myWeight); FillH1(GetName(name, src, step), x, wSignal); } @@ -205,7 +215,8 @@ string LmvmHist::GetName(const string& name, ELmvmSrc src, ELmvmAnaStep step) return GetName(GetName(name, src), step); } -void LmvmHist::SetOptH1(TH1D *hist, TString xAxisTitle=" ", TString yAxisTitle=" ", Int_t Ndevision=510, Int_t style=1, Float_t size=2, Int_t color=1, string opt) +void LmvmHist::SetOptH1(TH1D* hist, TString xAxisTitle = " ", TString yAxisTitle = " ", Int_t Ndevision = 510, + Int_t style = 1, Float_t size = 2, Int_t color = 1, string opt) { hist->GetXaxis()->SetTitle(xAxisTitle); hist->GetYaxis()->SetTitle(yAxisTitle); @@ -241,19 +252,20 @@ void LmvmHist::SetOptH1(TH1D *hist, TString xAxisTitle=" ", TString yAxisTitle=" else if (opt == "line") { hist->SetLineStyle(style); hist->SetLineColor(color); - hist->SetTitle(""); // TODO: with ""? + hist->SetTitle(""); // TODO: with ""? hist->SetLineWidth(size); } - else LOG(error) << "Option '" << opt << "' undefined. Choose 'marker' or 'line'." << std::endl; + else + LOG(error) << "Option '" << opt << "' undefined. Choose 'marker' or 'line'." << std::endl; } -void LmvmHist::SetOptCanvas(TCanvas *canvas) +void LmvmHist::SetOptCanvas(TCanvas* canvas) { gStyle->SetOptStat(0); gStyle->SetEndErrorSize(5); //gStyle->SetErrorX(0); // X errors of the data points set to be 0 - gStyle->SetLineStyleString(22,"80 18 12 18 12 12"); // special style for the line - gStyle->SetEndErrorSize(5); // define end width of error bars + gStyle->SetLineStyleString(22, "80 18 12 18 12 12"); // special style for the line + gStyle->SetEndErrorSize(5); // define end width of error bars gStyle->SetCanvasColor(10); gStyle->SetPadColor(10); canvas->SetLeftMargin(0.15); @@ -261,7 +273,7 @@ void LmvmHist::SetOptCanvas(TCanvas *canvas) canvas->SetTopMargin(0.05); canvas->SetBottomMargin(0.12); canvas->ToggleEventStatus(); - canvas->Range(-200,-10,1000,-2); + canvas->Range(-200, -10, 1000, -2); canvas->SetFillColor(0); canvas->SetBorderMode(0); canvas->SetBorderSize(0); @@ -273,7 +285,8 @@ void LmvmHist::SetOptCanvas(TCanvas *canvas) canvas->SetFrameBorderSize(0); } -void LmvmHist::SetLegend(vector<LmvmLegend> legendV, double textsize, Double_t x1, Double_t y1, Double_t x2, Double_t y2) +void LmvmHist::SetLegend(vector<LmvmLegend> legendV, double textsize, Double_t x1, Double_t y1, Double_t x2, + Double_t y2) { auto leg = new TLegend(x1, y1, x2, y2); leg->SetLineColor(10); @@ -375,31 +388,34 @@ TH2D* LmvmHist::CreateSignificanceH2(TH2D* signal, TH2D* bg, const string& name, return hsig; } -void LmvmHist::DrawAllGTracks(int dim, const string& cName, const string& hName, vector<string> xLabel, vector<string> yLabel, double min, double max) +void LmvmHist::DrawAllGTracks(int dim, const string& cName, const string& hName, vector<string> xLabel, + vector<string> yLabel, double min, double max) { TCanvas* c = fHM.CreateCanvas(cName, cName, 2400, 2400); c->Divide(4, 4); int i = 1; - for (const string& ptcl : fGTrackNames) { + for (const string& ptcl : fGTrackNames) { c->cd(i++); string hFullname = hName + "_" + ptcl; - DrawAll(dim, hFullname, fGTrackLatex[i-2], xLabel, yLabel, min, max); + DrawAll(dim, hFullname, fGTrackLatex[i - 2], xLabel, yLabel, min, max); } } -void LmvmHist::DrawAllCands(int dim, const string& cName, const string& hName, vector<string> xLabel, vector<string> yLabel, double min, double max) +void LmvmHist::DrawAllCands(int dim, const string& cName, const string& hName, vector<string> xLabel, + vector<string> yLabel, double min, double max) { TCanvas* c = fHM.CreateCanvas(cName, cName, 2400, 1800); c->Divide(4, 3); int i = 1; - for (const string& ptcl : fCandNames) { + for (const string& ptcl : fCandNames) { c->cd(i++); string hFullname = hName + "_" + ptcl; - DrawAll(dim, hFullname, fCandLatex[i-2], xLabel, yLabel, min, max); + DrawAll(dim, hFullname, fCandLatex[i - 2], xLabel, yLabel, min, max); } } -void LmvmHist::DrawAllCandsAndSteps(int dim, const string& cName, const string& hName, vector<string> xLabel, vector<string> yLabel, double min, double max) +void LmvmHist::DrawAllCandsAndSteps(int dim, const string& cName, const string& hName, vector<string> xLabel, + vector<string> yLabel, double min, double max) { for (auto step : fAnaSteps) { TCanvas* c = fHM.CreateCanvas(GetName(cName + "cands_all", step), GetName(cName + "cands_all", step), 2400, 1800); @@ -408,7 +424,7 @@ void LmvmHist::DrawAllCandsAndSteps(int dim, const string& cName, const string& for (const string& ptcl : fCandNames) { c->cd(i++); string hFullname = GetName(hName + "_" + ptcl, step); - DrawAll(dim, hFullname.c_str(), fCandLatex[i-2], xLabel, yLabel, min, max); + DrawAll(dim, hFullname.c_str(), fCandLatex[i - 2], xLabel, yLabel, min, max); } } @@ -424,7 +440,8 @@ void LmvmHist::DrawAllCandsAndSteps(int dim, const string& cName, const string& } } -void LmvmHist::DrawAll(int dim, const string& hFullname, const string& text, vector<string> xLabel, vector<string> yLabel, double min, double max) +void LmvmHist::DrawAll(int dim, const string& hFullname, const string& text, vector<string> xLabel, + vector<string> yLabel, double min, double max) { if (dim == 1) { TH1D* h = H1Clone(hFullname.c_str()); @@ -449,10 +466,11 @@ void LmvmHist::DrawAll(int dim, const string& hFullname, const string& text, vec for (size_t iL = 0; iL < yLabel.size(); iL++) { h->GetYaxis()->SetBinLabel(iL + 1, yLabel[iL].c_str()); } - } + } } - else LOG(error) << "LmvmHist::DrawAll: Choose dimension 1 or 2"; - + else + LOG(error) << "LmvmHist::DrawAll: Choose dimension 1 or 2"; + DrawTextOnPad(text, 0.25, 0.9, 0.75, 0.999); } diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmHist.h b/analysis/PWGDIL/dielectron/lmvm/LmvmHist.h index f36510336132fb251d5f834765948728cafc283f..6061ebfd03bbbde998e6c10f577804080ac66bef 100755 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmHist.h +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmHist.h @@ -49,7 +49,7 @@ public: const static int fNofCandNames = 10; const static std::vector<std::string> fCandNames; const static std::vector<std::string> fCandLatex; - + std::vector<std::string> CombineNames(const std::string& name, const std::vector<std::string>& subNames); @@ -134,21 +134,27 @@ public: std::string GetName(const std::string& name, ELmvmSrc src); std::string GetName(const std::string& name, ELmvmSrc src, ELmvmAnaStep step); - void SetOptH1(TH1D *hist, TString xAxisTitle, TString yAxisTitle, Int_t Ndevision, Int_t style, Float_t size, Int_t color, std::string opt = ""); // copied from Tetyanas macro - void SetOptCanvas(TCanvas *canvas); + void SetOptH1(TH1D* hist, TString xAxisTitle, TString yAxisTitle, Int_t Ndevision, Int_t style, Float_t size, + Int_t color, std::string opt = ""); // copied from Tetyanas macro + void SetOptCanvas(TCanvas* canvas); void SetLegend(std::vector<LmvmLegend>, double textsize, double x1, double y1, double x2, double y2); void Rebin(const std::string& name, int nGroup); // TODO: used? void Rebin(const std::string& name, const std::vector<std::string>& subNames, int nGroup); - void Rebin(const std::string& name, const std::vector<std::string>& subNames1, const std::vector<std::string>& subNames2, int nGroup); + void Rebin(const std::string& name, const std::vector<std::string>& subNames1, + const std::vector<std::string>& subNames2, int nGroup); TH1D* CreateSignificanceH1(TH1D* s, TH1D* bg, const std::string& name, const std::string& option); TH2D* CreateSignificanceH2(TH2D* signal, TH2D* bg, const std::string& name, const std::string& title); - void DrawAll(int dim, const std::string& hFullname, const std::string& padText, std::vector<std::string> xLabel, std::vector<std::string> yLabel, double min, double max); - void DrawAllGTracks(int dim, const std::string& cName, const std::string& hName, std::vector<std::string> xLabel, std::vector<std::string> yLabel, double min, double max); - void DrawAllCands(int dim, const std::string& cName, const std::string& hName, std::vector<std::string> xLabel, std::vector<std::string> yLabel, double min, double max); - void DrawAllCandsAndSteps(int dim, const std::string& cName, const std::string& hName, std::vector<std::string> xLabel, std::vector<std::string> yLabel, double min, double max); + void DrawAll(int dim, const std::string& hFullname, const std::string& padText, std::vector<std::string> xLabel, + std::vector<std::string> yLabel, double min, double max); + void DrawAllGTracks(int dim, const std::string& cName, const std::string& hName, std::vector<std::string> xLabel, + std::vector<std::string> yLabel, double min, double max); + void DrawAllCands(int dim, const std::string& cName, const std::string& hName, std::vector<std::string> xLabel, + std::vector<std::string> yLabel, double min, double max); + void DrawAllCandsAndSteps(int dim, const std::string& cName, const std::string& hName, + std::vector<std::string> xLabel, std::vector<std::string> yLabel, double min, double max); void WriteToFile(); diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmTask.cxx b/analysis/PWGDIL/dielectron/lmvm/LmvmTask.cxx index 149e86bac55bd8db6340d06540ddbde798483b02..ce2fa20be6d65f72605392ddbbc622f5bed2f098 100755 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmTask.cxx +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmTask.cxx @@ -14,16 +14,15 @@ #include "CbmRichPoint.h" #include "CbmRichRing.h" #include "CbmRichUtil.h" - #include "CbmStsHit.h" #include "CbmStsKFTrackFitter.h" #include "CbmStsTrack.h" -#include "CbmTrdHit.h" -#include "CbmTrdTrack.h" #include "CbmTofHit.h" #include "CbmTofPoint.h" #include "CbmTofTrack.h" #include "CbmTrackMatchNew.h" +#include "CbmTrdHit.h" +#include "CbmTrdTrack.h" #include "CbmVertex.h" #include "cbm/elid/CbmLitGlobalElectronId.h" #include "cbm/qa/mc/CbmLitMCTrackCreator.h" @@ -37,17 +36,17 @@ #include "TClonesArray.h" #include "TDatabasePDG.h" -#include "TMCProcess.h" #include "TFile.h" +#include "TMCProcess.h" #include "TRandom3.h" #include "TVector3.h" #include <sstream> #include <vector> +#include "LmvmHist.h" #include "LmvmSimParam.h" #include "LmvmUtils.h" -#include "LmvmHist.h" using namespace std; @@ -70,8 +69,10 @@ void LmvmTask::InitHists() fH.CreateH1("hMotherPdg", {"mc", "acc"}, "Pdg code", "Particles / Event", 7000, -3500., 3500.); fH.CreateH1("hCandPdg", fH.fAnaStepNames, "Pdg code", "Particles / Event", 7001, -3500., 3500.); - fH.CreateH2("hCandPdgVsMom", fH.fAnaStepNames, "P [GeV/c]", "Particle", "Yield/(Event * Bin)", 200, 0., 10., 6, 0., 6.); - fH.CreateH2("hCandElSrc", "Analysis step", "Mother of Electron Candidate", ax, fH.fNofAnaSteps, 0., fH.fNofAnaSteps, 8, 0., 8.); + fH.CreateH2("hCandPdgVsMom", fH.fAnaStepNames, "P [GeV/c]", "Particle", "Yield/(Event * Bin)", 200, 0., 10., 6, 0., + 6.); + fH.CreateH2("hCandElSrc", "Analysis step", "Mother of Electron Candidate", ax, fH.fNofAnaSteps, 0., fH.fNofAnaSteps, + 8, 0., 8.); fH.CreateH2("hBgPairPdg", fH.fAnaStepNames, "PDG of Candidate 1", "PDG of Candidate 2", ax, 8, 0., 8., 8, 0., 8.); fH.CreateH2("hPmtXY", fH.fSrcNames, "X [cm]", "Y [cm]", "Counter", 110, -110, 110, 200, -200, 200); @@ -79,15 +80,15 @@ void LmvmTask::InitHists() fH.CreateH2("hVertexGammaXZ", fH.fAnaStepNames, "Z [cm]", "X [cm]", ax, 200, -10., 190., 400, -130., 130.); fH.CreateH2("hVertexGammaYZ", fH.fAnaStepNames, "Z [cm]", "Y [cm]", ax, 200, -10., 190., 400, -130., 130.); fH.CreateH2("hVertexGammaXY", fH.fAnaStepNames, "X [cm]", "Y [cm]", ax, 400, -130., 130., 400, -130., 130.); - fH.CreateH2("hVertexGammaRZ", fH.fAnaStepNames, "Z [cm]", "#sqrt{X^{2}+Y^{2}} [cm]", ax, 110, fZ - 10., fZ + 100., 50, 0., - 50.); - + fH.CreateH2("hVertexGammaRZ", fH.fAnaStepNames, "Z [cm]", "#sqrt{X^{2}+Y^{2}} [cm]", ax, 110, fZ - 10., fZ + 100., 50, + 0., 50.); + // vertices of not-electron candidates that are misidentified as electrons in ToF fH.CreateH2("hVertexXZ_misidTof", fH.fAnaStepNames, "Z [cm]", "X [cm]", ax, 110, fZ - 10., fZ + 100., 100, -50., 50.); fH.CreateH2("hVertexYZ_misidTof", fH.fAnaStepNames, "Z [cm]", "Y [cm]", ax, 110, fZ - 10., fZ + 100., 100, -50., 50.); fH.CreateH2("hVertexXY_misidTof", fH.fAnaStepNames, "X [cm]", "Y [cm]", ax, 100, -50., 50., 100, -50., 50.); - fH.CreateH2("hVertexRZ_misidTof", fH.fAnaStepNames, "Z [cm]", "#sqrt{X^{2}+Y^{2}} [cm]", ax, 110, fZ - 10., fZ + 100., 100, 0., - 50.); + fH.CreateH2("hVertexRZ_misidTof", fH.fAnaStepNames, "Z [cm]", "#sqrt{X^{2}+Y^{2}} [cm]", ax, 110, fZ - 10., fZ + 100., + 100, 0., 50.); fH.CreateH1("hNofBgTracks", "Analysis step", "Tracks/event", fH.fNofAnaSteps, 0., fH.fNofAnaSteps); fH.CreateH1("hNofSignalTracks", "Analysis step", "Tracks/event", fH.fNofAnaSteps, 0., fH.fNofAnaSteps); fH.CreateH2("hBgSrcTracks", "Analysis step", "Candidate Source", ax, fH.fNofAnaSteps, 0., fH.fNofAnaSteps, 8, 0., 8.); @@ -96,13 +97,15 @@ void LmvmTask::InitHists() fH.CreateH1("hNofMismatches", {"all", "rich", "trd", "tof"}, "Analysis step", "Tracks/event", fH.fNofAnaSteps, 0., fH.fNofAnaSteps); fH.CreateH1("hNofMismatches_gTracks", {"all", "complete"}, fH.fGTrackNames, "Detector", ax, 7., 0., 7.); - fH.CreateH1("hNofMismatchedTrackSegments", {"all", "complete"}, fH.fGTrackNames, "nof mism. track segments", ax, 4., 0., 4.); //"all": also partial tracks; "complete": only tracks with entries in all detectors - fH.CreateH2("hMatchId_gTracks", fH.fGTrackNames, "Nof mismatched Track Segments", "Identification", ax, 4., 0., 4., 2., 0., 2.); - + fH.CreateH1("hNofMismatchedTrackSegments", {"all", "complete"}, fH.fGTrackNames, "nof mism. track segments", ax, 4., + 0., 4.); //"all": also partial tracks; "complete": only tracks with entries in all detectors + fH.CreateH2("hMatchId_gTracks", fH.fGTrackNames, "Nof mismatched Track Segments", "Identification", ax, 4., 0., 4., + 2., 0., 2.); + fH.CreateH1("hChi2_mismatch_all", {"rich", "trd", "tof"}, fH.fGTrackNames, "#chi^{2}", ax, 200, 0., 20.); fH.CreateH1("hChi2_truematch_all", {"rich", "trd", "tof"}, fH.fGTrackNames, "#chi^{2}", ax, 200, 0., 20.); fH.CreateH1("hChi2_mismatch_complete", {"rich", "trd", "tof"}, fH.fGTrackNames, "#chi^{2}", ax, 200, 0., 20.); - fH.CreateH1("hChi2_truematch_complete", {"rich", "trd", "tof"}, fH.fGTrackNames, "#chi^{2}", ax, 200, 0., 20.); + fH.CreateH1("hChi2_truematch_complete", {"rich", "trd", "tof"}, fH.fGTrackNames, "#chi^{2}", ax, 200, 0., 20.); fH.CreateH1("hNofGhosts", "Analysis step", "Tracks/event", fH.fNofAnaSteps, 0., fH.fNofAnaSteps); @@ -143,8 +146,10 @@ void LmvmTask::InitHists() fH.CreateH1("hMvdMcDist", {"1", "2"}, fH.fSrcNames, "Track-Hit distance [cm]", ax, 100, 0., 10.); fH.CreateH1("hMinv", fH.fSrcNames, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", axMinv, 250, 0., 2.5); - fH.CreateH1("hMinv_urqmdAll", fH.fSrcNames, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", axMinv, 250, 0., 2.5); // for UrQMD particles only - fH.CreateH1("hMinv_urqmdEl", fH.fSrcNames, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", axMinv, 250, 0., 2.5); // for UrQMD electrons only + fH.CreateH1("hMinv_urqmdAll", fH.fSrcNames, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", axMinv, 250, 0., + 2.5); // for UrQMD particles only + fH.CreateH1("hMinv_urqmdEl", fH.fSrcNames, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", axMinv, 250, 0., + 2.5); // for UrQMD electrons only // Pair Yield histograms for combinatorial BG calculation for (const string& comb : {"PM", "PP", "MM"}) { @@ -153,7 +158,7 @@ void LmvmTask::InitHists() fH.CreateH1(hName, {"sameEv", "mixedEv"}, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", axMinv, 250, 0., 2.5); } } - + fH.CreateH1("hMinvBgMatch", {"trueMatch", "trueMatchEl", "trueMatchNotEl", "mismatch"}, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", ax, 250, 0., 2.5); fH.CreateH1("hMinvBgSource", fH.fBgPairSrcNames, fH.fAnaStepNames, "M_{ee} [GeV/c^{2}]", axMinv, 250, 0., 2.5); @@ -188,15 +193,16 @@ void LmvmTask::InitHists() fH.CreateH1("hKaonMom", {"all", "prim"}, {"mc", "acc", "recSts", "recStsRich", "recStsRichTrd", "recStsRichTrdTof"}, "P [GeV/c]", ax, 100, 0., 10.); - fH.CreateH1("hPionSupp", {"pi", "idEl"},fH.fAnaStepNames, "P [GeV/c]", ax, 20, 0., 10.); + fH.CreateH1("hPionSupp", {"pi", "idEl"}, fH.fAnaStepNames, "P [GeV/c]", ax, 20, 0., 10.); fH.CreateH1("hMom_gTracks", fH.fGTrackNames, "P [GeV/c]", ax, 100, 0., 10.); fH.CreateH1("hMom_cands", fH.fCandNames, fH.fAnaStepNames, "P [GeV/c]", ax, 100, 0., 10.); fH.CreateH2("hPtY_gTracks", fH.fGTrackNames, "Rapidity", "P_{t} [GeV/c]", ax, 40, 0., 4., 20, 0., 2.); fH.CreateH2("hPtY_cands", fH.fCandNames, fH.fAnaStepNames, "Rapidity", "P_{t} [GeV/c]", ax, 40, 0., 4., 20, 0., 2.); fH.CreateH2("hTofM2_gTracks", fH.fGTrackNames, "P [GeV/c]", "m^{2} [GeV/c^{2}]^{2}", ax, 100, 0., 10., 500, -0.1, 1.); - fH.CreateH2("hTofM2_cands", fH.fCandNames, fH.fAnaStepNames, "P [GeV/c]", "m^{2} [GeV/c^{2}]^{2}", ax, 100, 0., 10., 500, -0.1, 1.); - + fH.CreateH2("hTofM2_cands", fH.fCandNames, fH.fAnaStepNames, "P [GeV/c]", "m^{2} [GeV/c^{2}]^{2}", ax, 100, 0., 10., + 500, -0.1, 1.); + fH.CreateH1("hTofPilePdgs_cands", fH.fAnaStepNames, "Particle", ax, 12, 0., 12.); fH.CreateH1("hTofPilePdgs_gTracks", "Particle", ax, 12, 0., 12.); fH.CreateH2("hTofPileHitXY", fH.fCandNames, "X [cm]", "Y [cm]", ax, 110., -550., 550., 90., -450., 450.); @@ -205,22 +211,30 @@ void LmvmTask::InitHists() fH.CreateH2("hTofPilePty_cands", fH.fCandNames, "Rapidity", "P_{t} [GeV/c]", ax, 40, 0., 4., 20, 0., 2.); // ToF Hits - fH.CreateH2("hTofPointXY", {"trueid", "misid", "truematch", "mismatch"}, fH.fGTrackNames, "X [cm]", "Y [cm]", ax, 110., -550., 550., 90., -450., 450.); // to see if maybe misid/mismatches stem from a certain region in ToF - fH.CreateH2("hTofHitXY", {"trueid", "misid", "truematch", "mismatch"}, fH.fGTrackNames, "X [cm]", "Y [cm]", ax, 110., -550., 550., 90., -450., 450.); - fH.CreateH1("hTofHitPointDist", {"trueid", "misid", "truematch", "mismatch"}, fH.fGTrackNames, "#sqrt{dX^{2} + dY^{2}} [cm]", ax, 400., 0., 20.); - fH.CreateH2("hTofHitTrackDist_gTracks", {"trueid", "misid", "truematch", "mismatch"}, fH.fGTrackNames, "P [GeV/c]", "Distance [cm]", ax, 100, 0., 10., 200., 0., 50.); - fH.CreateH2("hTofHitTrackDist_cands", fH.fCandNames, "P [GeV/c]", "Distance [cm]", ax, 100, 0., 10., 200., 0., 50.); + fH.CreateH2("hTofPointXY", {"trueid", "misid", "truematch", "mismatch"}, fH.fGTrackNames, "X [cm]", "Y [cm]", ax, + 110., -550., 550., 90., -450., + 450.); // to see if maybe misid/mismatches stem from a certain region in ToF + fH.CreateH2("hTofHitXY", {"trueid", "misid", "truematch", "mismatch"}, fH.fGTrackNames, "X [cm]", "Y [cm]", ax, 110., + -550., 550., 90., -450., 450.); + fH.CreateH1("hTofHitPointDist", {"trueid", "misid", "truematch", "mismatch"}, fH.fGTrackNames, + "#sqrt{dX^{2} + dY^{2}} [cm]", ax, 400., 0., 20.); + fH.CreateH2("hTofHitTrackDist_gTracks", {"trueid", "misid", "truematch", "mismatch"}, fH.fGTrackNames, "P [GeV/c]", + "Distance [cm]", ax, 100, 0., 10., 200., 0., 50.); + fH.CreateH2("hTofHitTrackDist_cands", fH.fCandNames, "P [GeV/c]", "Distance [cm]", ax, 100, 0., 10., 200., 0., 50.); fH.CreateH1("hMomRatio_cands", fH.fCandNames, fH.fAnaStepNames, "Ratio P_{rec}/P_{MC}", ax, 100, 0., 2.); - fH.CreateH2("hMomRatioVsMom_cands", fH.fCandNames, fH.fAnaStepNames, "P_{MC} [GeV/c]", "Ratio P_{rec}/P_{MC}", ax, 100, 0., 10., 200., 0., 2.); + fH.CreateH2("hMomRatioVsMom_cands", fH.fCandNames, fH.fAnaStepNames, "P_{MC} [GeV/c]", "Ratio P_{rec}/P_{MC}", ax, + 100, 0., 10., 200., 0., 2.); // compare misid. candidates (pions, proton, ...) with not misid. ones for (size_t iP = 4; iP < fH.fCandNames.size(); iP++) { - fH.CreateH1("hMom_" + fH.fCandNames[iP], {"true", "misid"}, "P [GeV/c]", ax, 100, 0., 10.); - fH.CreateH2("hPtY_" + fH.fCandNames[iP], {"true", "misid"}, "Rapidity", "P_{t} [GeV/c]", ax, 40, 0., 4., 40, 0., 2.); - fH.CreateH2("hTofM2_" + fH.fCandNames[iP], {"true", "misid"}, "P [GeV/c]", "m^{2} [GeV/c^{2}]^{2}", ax, 100, 0., 10., 500, -0.1, 1.); + fH.CreateH1("hMom_" + fH.fCandNames[iP], {"true", "misid"}, "P [GeV/c]", ax, 100, 0., 10.); + fH.CreateH2("hPtY_" + fH.fCandNames[iP], {"true", "misid"}, "Rapidity", "P_{t} [GeV/c]", ax, 40, 0., 4., 40, 0., + 2.); + fH.CreateH2("hTofM2_" + fH.fCandNames[iP], {"true", "misid"}, "P [GeV/c]", "m^{2} [GeV/c^{2}]^{2}", ax, 100, 0., + 10., 500, -0.1, 1.); } - + // beta-momentum spectra fH.CreateH2("hBetaMom", {"gTracks", "cands"}, fH.fCandNames, "P * Q ", "beta", ax, 200, -10., 10., 150., 0., 1.5); fH.CreateH2("hBetaMom_allGTracks", "P * Q ", "beta", ax, 200, -10., 10., 150., 0., 1.5); @@ -232,18 +246,26 @@ void LmvmTask::InitHists() fH.CreateH2("hIndexStsTrd", fH.fGTrackNames, "STS Index", "TRD Index", ax, 400., 0., 400., 400., 0., 400.); fH.CreateH2("hIndexStsTof", fH.fGTrackNames, "STS Index", "ToF Index", ax, 400., 0., 400., 800., 0., 800.);*/ - fH.CreateH2("hPdgVsMom_gTracks", {"rich", "trd", "tof"}, {"all", "complete"}, "P [GeV/c]", "Misid. Particle", "Yield/(Event * Bin)", 200, 0., 10., 6, 0., 6.); + fH.CreateH2("hPdgVsMom_gTracks", {"rich", "trd", "tof"}, {"all", "complete"}, "P [GeV/c]", "Misid. Particle", + "Yield/(Event * Bin)", 200, 0., 10., 6, 0., 6.); fH.CreateH2("hTofM2Calc_gTracks", fH.fGTrackNames, "Time", "", "Yield", 100., 5., 15., 4., 0., 4.); - fH.CreateH2("hTofTimeVsMom_gTracks", {"trueid", "misid", "truematch", "mismatch"}, fH.fGTrackNames, "P [GeV/c]", "Time (ct) [m]", "Yield", 100., 0., 10., 100., 5., 15.); - fH.CreateH2("hTofTimeVsMom_cands", fH.fCandNames, "P [GeV/c]", "Time (ct) [m]", "Yield", 100., 0., 10., 100., 5., 15.); - fH.CreateH2("hRichRingTrackDist_gTracks", fH.fGTrackNames, "P [GeV/c]", "D [cm]", "Yield", 100., 0., 10., 200., 0., 20.); - fH.CreateH2("hRichRingTrackDist_cands", fH.fCandNames, fH.fAnaStepNames, "P [GeV/c]", "D [cm]", "Yield", 100., 0., 10., 200., 0., 20.); + fH.CreateH2("hTofTimeVsMom_gTracks", {"trueid", "misid", "truematch", "mismatch"}, fH.fGTrackNames, "P [GeV/c]", + "Time (ct) [m]", "Yield", 100., 0., 10., 100., 5., 15.); + fH.CreateH2("hTofTimeVsMom_cands", fH.fCandNames, "P [GeV/c]", "Time (ct) [m]", "Yield", 100., 0., 10., 100., 5., + 15.); + fH.CreateH2("hRichRingTrackDist_gTracks", fH.fGTrackNames, "P [GeV/c]", "D [cm]", "Yield", 100., 0., 10., 200., 0., + 20.); + fH.CreateH2("hRichRingTrackDist_cands", fH.fCandNames, fH.fAnaStepNames, "P [GeV/c]", "D [cm]", "Yield", 100., 0., + 10., 200., 0., 20.); // check chi2 of candidates after El-ID step - fH.CreateH2("hChi2VsMom", {"sts", "rich", "trd", "tof"}, fH.fCandNames, "P [GeV/c]", "#chi^{2}", ax, 100., 0., 10., 200, 0., 20.); - fH.CreateH2("hTofTimeVsChi2", {"sts", "rich", "trd", "tof"}, fH.fCandNames, "#chi^{2}", "Time (ct) [m]", ax, 200., 0., 20., 200, 0., 20.); - fH.CreateH2("hChi2Comb", {"StsRich", "StsTrd", "StsTof", "RichTrd", "RichTof", "TrdTof"}, fH.fCandNames, "#chi^{2}_{1}", "#chi^{2}_{2}", ax, 200., 0., 20., 200, 0., 20.); + fH.CreateH2("hChi2VsMom", {"sts", "rich", "trd", "tof"}, fH.fCandNames, "P [GeV/c]", "#chi^{2}", ax, 100., 0., 10., + 200, 0., 20.); + fH.CreateH2("hTofTimeVsChi2", {"sts", "rich", "trd", "tof"}, fH.fCandNames, "#chi^{2}", "Time (ct) [m]", ax, 200., 0., + 20., 200, 0., 20.); + fH.CreateH2("hChi2Comb", {"StsRich", "StsTrd", "StsTof", "RichTrd", "RichTof", "TrdTof"}, fH.fCandNames, + "#chi^{2}_{1}", "#chi^{2}_{2}", ax, 200., 0., 20., 200, 0., 20.); } InitStatus LmvmTask::Init() @@ -273,8 +295,8 @@ InitStatus LmvmTask::Init() //fTofTracks = InitOrFatal<TClonesArray>("TofTracks"); // TODO: check this and next lines FairRootManager* fairRootMan = FairRootManager::Instance(); - fTofTracks = (TClonesArray*) fairRootMan->GetObject("TofTrack"); - if (fTofTracks == nullptr) { LOG(warning) << "LmvmTask::Init : no TofTrack array!"; } + fTofTracks = (TClonesArray*) fairRootMan->GetObject("TofTrack"); + if (fTofTracks == nullptr) { LOG(warning) << "LmvmTask::Init : no TofTrack array!"; } InitHists(); @@ -344,8 +366,10 @@ double LmvmTask::MinvScale(double minv) { double scale = -1.; if (fParticle == "inmed") scale = LmvmUtils::GetMassScaleInmed(minv); - else if (fParticle == "qgp") scale = LmvmUtils::GetMassScaleQgp(minv); - else scale = 1.; + else if (fParticle == "qgp") + scale = LmvmUtils::GetMassScaleQgp(minv); + else + scale = 1.; return scale; } @@ -383,7 +407,7 @@ void LmvmTask::DoMcTrack() string chargeStr = (mct->GetCharge() > 0) ? "+" : "-"; bool isAcc = IsMcTrackAccepted(i); double mom = mct->GetP(); - double w = (LmvmUtils::IsMcSignalEl(mct)) ? fW : 1.; + double w = (LmvmUtils::IsMcSignalEl(mct)) ? fW : 1.; FillMomHists(mct, nullptr, src, ELmvmAnaStep::Mc); if (isAcc) FillMomHists(mct, nullptr, src, ELmvmAnaStep::Acc); @@ -428,10 +452,10 @@ void LmvmTask::DoMcTrack() void LmvmTask::DoMcPair() { int nMcTracks = fMCTracks->GetEntries(); - for (int iMc1 = 0; iMc1 < nMcTracks; iMc1++) { // TODO: range until iMc1 < nMcTracks-1 ?? + for (int iMc1 = 0; iMc1 < nMcTracks; iMc1++) { // TODO: range until iMc1 < nMcTracks-1 ?? CbmMCTrack* mct1 = static_cast<CbmMCTrack*>(fMCTracks->At(iMc1)); ELmvmSrc src = LmvmUtils::GetMcSrc(mct1, fMCTracks); - + // To speed up: select only signal, eta and pi0 electrons if (!(src == ELmvmSrc::Signal || src == ELmvmSrc::Pi0 || src == ELmvmSrc::Eta)) continue; @@ -481,7 +505,7 @@ void LmvmTask::RichPmtXY() TVector3 v; mct->GetStartVertex(v); ELmvmSrc src = LmvmUtils::GetMcSrc(mct, fMCTracks); - double w = (LmvmUtils::IsMcSignalEl(mct)) ? fW : 1.; + double w = (LmvmUtils::IsMcSignalEl(mct)) ? fW : 1.; if (v.Z() < 2.) { fH.FillH2("hPmtXY", src, richHit->GetX(), richHit->GetY(), w); } } } @@ -501,8 +525,8 @@ void LmvmTask::AnalyseGlobalTracks() for (int i = 0; i < nofMcTracks; i++) { CbmMCTrack* mct = static_cast<CbmMCTrack*>(fMCTracks->At(i)); - int pdg = mct->GetPdgCode(); - bool isAccSts = (mct->GetNPoints(ECbmModuleId::kMvd) + mct->GetNPoints(ECbmModuleId::kSts) >= 4); + int pdg = mct->GetPdgCode(); + bool isAccSts = (mct->GetNPoints(ECbmModuleId::kMvd) + mct->GetNPoints(ECbmModuleId::kSts) >= 4); TVector3 vertex; mct->GetStartVertex(vertex); @@ -554,7 +578,7 @@ void LmvmTask::AnalyseGlobalTracks() if (isRich) fH.FillH1(hName + "Mom_all_recStsRich", p); if (isRich && isTrd) fH.FillH1(hName + "Mom_all_recStsRichTrd", p); if (isRich && isTrd && isTof) fH.FillH1(hName + "Mom_all_recStsRichTrdTof", p); - + if (isPrim) { fH.FillH1(hName + "Mom_prim_recSts", p); if (isRich) fH.FillH1(hName + "Mom_prim_recStsRich", p); @@ -570,34 +594,44 @@ void LmvmTask::AnalyseGlobalTracks() fH.FillH2("hIndexStsTrd_" + ptcl, stsInd, gTrack->GetTrdTrackIndex()); fH.FillH2("hIndexStsTof_" + ptcl, stsInd, gTrack->GetTofHitIndex());*/ - bool isTofEl = (CbmLitGlobalElectronId::GetInstance().IsTofElectron(i, p)) ? true : false; - bool isElectron = (CbmLitGlobalElectronId::GetInstance().IsRichElectron(i, p) && CbmLitGlobalElectronId::GetInstance().IsTrdElectron(i, p) && CbmLitGlobalElectronId::GetInstance().IsTofElectron(i, p)) ? true : false; - - if (gTrack != nullptr) { - CheckMismatches(gTrack, pdg, isElectron, ptcl, w); - } + bool isTofEl = (CbmLitGlobalElectronId::GetInstance().IsTofElectron(i, p)) ? true : false; + bool isElectron = (CbmLitGlobalElectronId::GetInstance().IsRichElectron(i, p) + && CbmLitGlobalElectronId::GetInstance().IsTrdElectron(i, p) + && CbmLitGlobalElectronId::GetInstance().IsTofElectron(i, p)) + ? true + : false; + + if (gTrack != nullptr) { CheckMismatches(gTrack, pdg, isElectron, ptcl, w); } if (std::abs(pdg) != 11) PidVsMom(gTrack, i, pdg, p); - Double_t richDist = CbmRichUtil::GetRingTrackDistance(i); + Double_t richDist = CbmRichUtil::GetRingTrackDistance(i); fH.FillH2("hRichRingTrackDist_gTracks_" + ptcl, p, richDist); // investigate misidentifications: compare misidentified candidates with same not-misidentified particles - // first check if global track has entries in all detectors - if (!IsInAllDets(gTrack)) continue; // Mind: all following actions are done only for fully rec. gTracks + // first check if global track has entries in all detectors + if (!IsInAllDets(gTrack)) continue; // Mind: all following actions are done only for fully rec. gTracks + + CheckTofIdentification(gTrack, ptcl, p, m2, pdg, isTofEl); - CheckTofIdentification(gTrack, ptcl, p, m2, pdg, isTofEl); - - fH.FillH1("hMom_gTracks_" + ptcl, p); - fH.FillH2("hPtY_gTracks_" + ptcl, rap, pt); + fH.FillH1("hMom_gTracks_" + ptcl, p); + fH.FillH2("hPtY_gTracks_" + ptcl, rap, pt); fH.FillH2("hTofM2_gTracks_" + ptcl, p, m2); - // check PIDs in "Tof pile" + // check PIDs in "Tof pile" if (p > 0.3 && p < 1. && m2 > -0.012 && m2 < 0.01) { bool isSignal = (mct != nullptr && mct->GetGeantProcessId() == kPPrimary && std::abs(mct->GetPdgCode()) == 11); - double pdgBin = (pdg == 11 && isSignal) ? 0.5 : (pdg == -11 && isSignal) ? 1.5 - : (pdg == 11 && !isSignal && isPrim) ? 2.5 : (pdg == -11 && !isSignal && isPrim) ? 3.5 - : (pdg == 11 && !isSignal && !isPrim) ? 4.5 : (pdg == -11 && !isSignal && !isPrim) ? 5.5 - : (pdg == 211) ? 6.5 : (pdg == -211) ? 7.5 : (pdg == 2212) ? 8.5 : (pdg == 321) ? 9.5 : (pdg == -321) ? 10.5 : 11.5; + double pdgBin = (pdg == 11 && isSignal) ? 0.5 + : (pdg == -11 && isSignal) ? 1.5 + : (pdg == 11 && !isSignal && isPrim) ? 2.5 + : (pdg == -11 && !isSignal && isPrim) ? 3.5 + : (pdg == 11 && !isSignal && !isPrim) ? 4.5 + : (pdg == -11 && !isSignal && !isPrim) ? 5.5 + : (pdg == 211) ? 6.5 + : (pdg == -211) ? 7.5 + : (pdg == 2212) ? 8.5 + : (pdg == 321) ? 9.5 + : (pdg == -321) ? 10.5 + : 11.5; fH.FillH1("hTofPilePdgs_gTracks", pdgBin); } @@ -607,7 +641,8 @@ void LmvmTask::AnalyseGlobalTracks() fH.FillH2("hPtY_" + ptcl2 + "_true", rap, pt); fH.FillH2("hTofM2_" + ptcl2 + "_true", p, m2); } - else if (isElectron && !(ptcl2 == "plutoEl+" || ptcl2 == "plutoEl-" || ptcl2 == "urqmdEl+" || ptcl2 == "urqmdEl-")) { + else if (isElectron + && !(ptcl2 == "plutoEl+" || ptcl2 == "plutoEl-" || ptcl2 == "urqmdEl+" || ptcl2 == "urqmdEl-")) { fH.FillH1("hMom_" + ptcl2 + "_misid", p); fH.FillH2("hPtY_" + ptcl2 + "_misid", rap, pt); fH.FillH2("hTofM2_" + ptcl2 + "_misid", p, m2); @@ -617,29 +652,32 @@ void LmvmTask::AnalyseGlobalTracks() } } -bool LmvmTask::IsInTofPile(double p, double m2) -{ - return (p > 0.3 && p < 1. && m2 > -0.012 && m2 < 0.01); -} +bool LmvmTask::IsInTofPile(double p, double m2) { return (p > 0.3 && p < 1. && m2 > -0.012 && m2 < 0.01); } void LmvmTask::PidVsMom(const CbmGlobalTrack* gTrack, int iGTrack, int pdg, double p) { bool isRichEl = CbmLitGlobalElectronId::GetInstance().IsRichElectron(iGTrack, p); bool isTrdEl = CbmLitGlobalElectronId::GetInstance().IsTrdElectron(iGTrack, p); bool isTofEl = CbmLitGlobalElectronId::GetInstance().IsTofElectron(iGTrack, p); - - double pdgBin = (pdg == 211) ? 0.5 : (pdg == -211) ? 1.5 : (pdg == 2212) ? 2.5 : (pdg == 321) ? 3.5 : (pdg == -321) ? 4.5 : 5.5; - + + double pdgBin = (pdg == 211) ? 0.5 + : (pdg == -211) ? 1.5 + : (pdg == 2212) ? 2.5 + : (pdg == 321) ? 3.5 + : (pdg == -321) ? 4.5 + : 5.5; + if (isRichEl) fH.FillH2("hPdgVsMom_gTracks_rich_all", p, pdgBin); - if (isTrdEl) fH.FillH2("hPdgVsMom_gTracks_trd_all", p, pdgBin); - if (isTofEl) fH.FillH2("hPdgVsMom_gTracks_tof_all", p, pdgBin); + if (isTrdEl) fH.FillH2("hPdgVsMom_gTracks_trd_all", p, pdgBin); + if (isTofEl) fH.FillH2("hPdgVsMom_gTracks_tof_all", p, pdgBin); if (isRichEl && IsInAllDets(gTrack)) fH.FillH2("hPdgVsMom_gTracks_rich_complete", p, pdgBin); - if (isTrdEl && IsInAllDets(gTrack)) fH.FillH2("hPdgVsMom_gTracks_trd_complete", p, pdgBin); - if (isTofEl && IsInAllDets(gTrack)) fH.FillH2("hPdgVsMom_gTracks_tof_complete", p, pdgBin); + if (isTrdEl && IsInAllDets(gTrack)) fH.FillH2("hPdgVsMom_gTracks_trd_complete", p, pdgBin); + if (isTofEl && IsInAllDets(gTrack)) fH.FillH2("hPdgVsMom_gTracks_tof_complete", p, pdgBin); } -void LmvmTask::CheckTofIdentification(const CbmGlobalTrack* gTrack, const string& gtString, double pMc, double m2, int pdg, bool isTofEl) +void LmvmTask::CheckTofIdentification(const CbmGlobalTrack* gTrack, const string& gtString, double pMc, double m2, + int pdg, bool isTofEl) { // Get STS ID int stsInd = gTrack->GetStsTrackIndex(); @@ -655,7 +693,7 @@ void LmvmTask::CheckTofIdentification(const CbmGlobalTrack* gTrack, const string CbmTofHit* tofHit = (CbmTofHit*) fTofHits->At(tofInd); if (tofHit == nullptr) return; CbmTofTrack* tofTrack = static_cast<CbmTofTrack*>(fTofTracks->At(gTrack->GetStsTrackIndex())); - if (tofTrack == nullptr) return; + if (tofTrack == nullptr) return; CbmMatch* tofHitMatch = static_cast<CbmMatch*>(fTofHitsMatches->At(tofInd)); if (tofHitMatch == nullptr) return; int tofPointIndex = tofHitMatch->GetMatchedLink().GetIndex(); @@ -666,7 +704,7 @@ void LmvmTask::CheckTofIdentification(const CbmGlobalTrack* gTrack, const string bool isTofMismatch = (tofMcTrackId == stsMcTrackId) ? false : true; bool isTofMisid = (std::abs(pdg) == 11 && isTofEl) ? false : (std::abs(pdg) != 11 && !isTofEl) ? false : true; - + double pointX = tofPoint->GetX(); double pointY = tofPoint->GetY(); double hitX = tofHit->GetX(); @@ -703,29 +741,34 @@ void LmvmTask::CheckTofIdentification(const CbmGlobalTrack* gTrack, const string // check ToF data for particles in ToF pile // TODO: needed or can be deleted? if (isTofEl && std::abs(pdg) != 11 && IsInTofPile(pMc, m2)) { - string candString = (pdg == 211) ? fH.fCandNames[4] : (pdg == -211) ? fH.fCandNames[5] : (pdg == 2212) ? fH.fCandNames[6] : (pdg == 321) ? fH.fCandNames[7] : (pdg == -321) ? fH.fCandNames[8] : fH.fCandNames[9]; + string candString = (pdg == 211) ? fH.fCandNames[4] + : (pdg == -211) ? fH.fCandNames[5] + : (pdg == 2212) ? fH.fCandNames[6] + : (pdg == 321) ? fH.fCandNames[7] + : (pdg == -321) ? fH.fCandNames[8] + : fH.fCandNames[9]; fH.FillH2("hTofPileHitXY_" + candString, hitX, hitY); fH.FillH2("hTofPilePointXY_" + candString, pointX, pointY); fH.FillH1("hTofPileHitPointDist_" + candString, dR); - } + } } void LmvmTask::BetaMom(const CbmMCTrack* mct, const CbmGlobalTrack* gTrack, const string& ptcl) { - Int_t tofInd = gTrack->GetTofHitIndex(); - if (tofInd < 0) return; - CbmTofHit* tofHit = (CbmTofHit*) fTofHits->At(tofInd); - if (NULL == tofHit) return; - double eventTime = FairRunAna::Instance()->GetEventHeader()->GetEventTime(); - Double_t noOffsetTime = tofHit->GetTime() - eventTime; - - double p = mct->GetP(); - double l = gTrack->GetLength(); - double t = 0.2998 * noOffsetTime; // time in ns -> transfrom to ct in m - double q = (mct->GetCharge() > 0) ? 1 : (mct->GetCharge() < 0) ? -1. : 0.; - double beta = l / (t * 100); // cm to m - fH.FillH2("hBetaMom_gTracks_" + ptcl, q * p, beta); - fH.FillH2("hBetaMom_allGTracks", q * p, beta); + Int_t tofInd = gTrack->GetTofHitIndex(); + if (tofInd < 0) return; + CbmTofHit* tofHit = (CbmTofHit*) fTofHits->At(tofInd); + if (NULL == tofHit) return; + double eventTime = FairRunAna::Instance()->GetEventHeader()->GetEventTime(); + Double_t noOffsetTime = tofHit->GetTime() - eventTime; + + double p = mct->GetP(); + double l = gTrack->GetLength(); + double t = 0.2998 * noOffsetTime; // time in ns -> transfrom to ct in m + double q = (mct->GetCharge() > 0) ? 1 : (mct->GetCharge() < 0) ? -1. : 0.; + double beta = l / (t * 100); // cm to m + fH.FillH2("hBetaMom_gTracks_" + ptcl, q * p, beta); + fH.FillH2("hBetaMom_allGTracks", q * p, beta); } void LmvmTask::CheckMismatches(const CbmGlobalTrack* gTrack, int pdg, bool isElectron, const string& ptcl, double w) @@ -739,7 +782,7 @@ void LmvmTask::CheckMismatches(const CbmGlobalTrack* gTrack, int pdg, bool isEle if (stsMcTrackId < 0) return; bool isFull = IsInAllDets(gTrack); - + // first calculate Chi2 CbmL1PFFitter fPFFitter; CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(stsInd)); @@ -753,21 +796,21 @@ void LmvmTask::CheckMismatches(const CbmGlobalTrack* gTrack, int pdg, bool isEle CbmTrackMatchNew* richMatch = nullptr; CbmTrackMatchNew* trdMatch = nullptr; - CbmMatch* tofMatch = nullptr; + CbmMatch* tofMatch = nullptr; int richRingInd = gTrack->GetRichRingIndex(); int trdTrackInd = gTrack->GetTrdTrackIndex(); int tofHitInd = gTrack->GetTofHitIndex(); if (richRingInd >= 0) richMatch = static_cast<CbmTrackMatchNew*>(fRichRingMatches->At(richRingInd)); - if (trdTrackInd >= 0) trdMatch = static_cast<CbmTrackMatchNew*>(fTrdTrackMatches->At(trdTrackInd)); - if (tofHitInd >= 0) tofMatch = static_cast<CbmMatch*>(fTofHitsMatches->At(tofHitInd)); + if (trdTrackInd >= 0) trdMatch = static_cast<CbmTrackMatchNew*>(fTrdTrackMatches->At(trdTrackInd)); + if (tofHitInd >= 0) tofMatch = static_cast<CbmMatch*>(fTofHitsMatches->At(tofHitInd)); int nofMismTrackSegmentsAll = 0; int nofMismTrackSegmentsComp = 0; fH.FillH1("hNofMismatches_gTracks_all_" + ptcl, 0.5, w); - if (isFull) fH.FillH1("hNofMismatches_gTracks_complete_" + ptcl, 0.5, w); + if (isFull) fH.FillH1("hNofMismatches_gTracks_complete_" + ptcl, 0.5, w); if (richMatch != nullptr) { int richMcTrackId = richMatch->GetMatchedLink().GetIndex(); @@ -779,8 +822,9 @@ void LmvmTask::CheckMismatches(const CbmGlobalTrack* gTrack, int pdg, bool isEle nofMismTrackSegmentsAll++; } if (isFull) { - if (richMcTrackId >= 0) fH.FillH1("hNofMismatches_gTracks_complete_" + ptcl, 1.5, w); - if (richMcTrackId >= 0 && richMcTrackId == stsMcTrackId) fH.FillH1("hChi2_truematch_complete_rich_" + ptcl, chi2, w); + if (richMcTrackId >= 0) fH.FillH1("hNofMismatches_gTracks_complete_" + ptcl, 1.5, w); + if (richMcTrackId >= 0 && richMcTrackId == stsMcTrackId) + fH.FillH1("hChi2_truematch_complete_rich_" + ptcl, chi2, w); if (richMcTrackId >= 0 && richMcTrackId != stsMcTrackId) { fH.FillH1("hChi2_mismatch_complete_rich_" + ptcl, chi2, w); fH.FillH1("hNofMismatches_gTracks_complete_" + ptcl, 2.5, w); @@ -790,7 +834,7 @@ void LmvmTask::CheckMismatches(const CbmGlobalTrack* gTrack, int pdg, bool isEle } if (trdMatch != nullptr) { - int trdMcTrackId = trdMatch->GetMatchedLink().GetIndex(); + int trdMcTrackId = trdMatch->GetMatchedLink().GetIndex(); if (trdMcTrackId >= 0) fH.FillH1("hNofMismatches_gTracks_all_" + ptcl, 3.5, w); if (trdMcTrackId >= 0 && trdMcTrackId == stsMcTrackId) fH.FillH1("hChi2_truematch_all_trd_" + ptcl, chi2, w); if (trdMcTrackId >= 0 && trdMcTrackId != stsMcTrackId) { @@ -805,10 +849,10 @@ void LmvmTask::CheckMismatches(const CbmGlobalTrack* gTrack, int pdg, bool isEle fH.FillH1("hNofMismatches_gTracks_complete_" + ptcl, 4.5, w); fH.FillH1("hChi2_mismatch_complete_trd_" + ptcl, chi2, w); nofMismTrackSegmentsComp++; - } + } } } - + if (tofMatch != nullptr && tofMatch->GetMatchedLink().GetIndex() >= 0) { FairMCPoint* tofPoint = static_cast<FairMCPoint*>(fTofPoints->At(tofMatch->GetMatchedLink().GetIndex())); if (tofPoint != nullptr) { @@ -822,12 +866,13 @@ void LmvmTask::CheckMismatches(const CbmGlobalTrack* gTrack, int pdg, bool isEle } if (isFull) { if (tofMcTrackId >= 0) fH.FillH1("hNofMismatches_gTracks_complete_" + ptcl, 5.5, w); - if (tofMcTrackId >= 0 && tofMcTrackId == stsMcTrackId) fH.FillH1("hChi2_truematch_complete_tof_" + ptcl, chi2, w); + if (tofMcTrackId >= 0 && tofMcTrackId == stsMcTrackId) + fH.FillH1("hChi2_truematch_complete_tof_" + ptcl, chi2, w); if (tofMcTrackId >= 0 && tofMcTrackId != stsMcTrackId) { fH.FillH1("hNofMismatches_gTracks_complete_" + ptcl, 6.5, w); fH.FillH1("hChi2_mismatch_complete_tof_" + ptcl, chi2, w); nofMismTrackSegmentsComp++; - } + } } } } @@ -835,9 +880,9 @@ void LmvmTask::CheckMismatches(const CbmGlobalTrack* gTrack, int pdg, bool isEle fH.FillH1("hNofMismatchedTrackSegments_all_" + ptcl, nofMismTrackSegmentsAll + 0.5, w); if (isFull) fH.FillH1("hNofMismatchedTrackSegments_complete_" + ptcl, nofMismTrackSegmentsComp + 0.5, w); - bool isTrueId = (std::abs(pdg) == 11 && isElectron) ? true : (std::abs(pdg) != 11 && !isElectron) ? true : false; - double binX = nofMismTrackSegmentsComp + 0.5; - double binY = (isTrueId == true) ? 0.5 : 1.5; + bool isTrueId = (std::abs(pdg) == 11 && isElectron) ? true : (std::abs(pdg) != 11 && !isElectron) ? true : false; + double binX = nofMismTrackSegmentsComp + 0.5; + double binY = (isTrueId == true) ? 0.5 : 1.5; fH.FillH2("hMatchId_gTracks_" + ptcl, binX, binY, w); } @@ -845,7 +890,7 @@ void LmvmTask::CheckMismatches(const CbmGlobalTrack* gTrack, int pdg, bool isEle bool LmvmTask::IsInAllDets(const CbmGlobalTrack* gTrack) { if (gTrack == nullptr) return false; - + int stsInd = gTrack->GetStsTrackIndex(); if (stsInd < 0) return false; CbmStsTrack* stsTrack = static_cast<CbmStsTrack*>(fStsTracks->At(stsInd)); @@ -858,7 +903,7 @@ bool LmvmTask::IsInAllDets(const CbmGlobalTrack* gTrack) CbmRichRing* richRing = static_cast<CbmRichRing*>(fRichRings->At(richInd)); CbmTrdTrack* trdTrack = static_cast<CbmTrdTrack*>(fTrdTracks->At(trdInd)); - CbmTofHit* tofHit = static_cast<CbmTofHit*>(fTofHits->At(tofInd)); + CbmTofHit* tofHit = static_cast<CbmTofHit*>(fTofHits->At(tofInd)); if (richRing == nullptr || trdTrack == nullptr || tofHit == nullptr) return false; return true; @@ -883,8 +928,8 @@ void LmvmTask::FillTopologyCands() cand.fRichInd = gTrack->GetRichRingIndex(); cand.fTrdInd = gTrack->GetTrdTrackIndex(); - cand.fTofInd = gTrack->GetTofHitIndex(); // TODO: is TofHitIndex; change this name everywhere - cand.fTofTrackInd = gTrack->GetStsTrackIndex(); // TODO: is this right?? + cand.fTofInd = gTrack->GetTofHitIndex(); // TODO: is TofHitIndex; change this name everywhere + cand.fTofTrackInd = gTrack->GetStsTrackIndex(); // TODO: is this right?? LmvmUtils::CalculateAndSetTrackParams(&cand, stsTrack, fKFVertex); cand.fIsChi2Prim = fCuts.IsChi2PrimaryOk(cand.fChi2Prim); @@ -1002,7 +1047,7 @@ void LmvmTask::FillCands() CbmRichRing* richRing = static_cast<CbmRichRing*>(fRichRings->At(cand.fRichInd)); CbmTrdTrack* trdTrack = static_cast<CbmTrdTrack*>(fTrdTracks->At(cand.fTrdInd)); - CbmTofHit* tofHit = static_cast<CbmTofHit*>(fTofHits->At(cand.fTofInd)); + CbmTofHit* tofHit = static_cast<CbmTofHit*>(fTofHits->At(cand.fTofInd)); CbmTofTrack* tofTrack = static_cast<CbmTofTrack*>(fTofTracks->At(cand.fTofTrackInd)); if (richRing == nullptr || trdTrack == nullptr || tofHit == nullptr || tofTrack == nullptr) continue; @@ -1035,17 +1080,19 @@ void LmvmTask::CombinatorialPairs() size_t nCand = fCandsTotal.size(); for (size_t iC1 = 0; iC1 < nCand - 1; iC1++) { const auto& cand1 = fCandsTotal[iC1]; - + for (size_t iC2 = iC1 + 1; iC2 < nCand; iC2++) { const auto& cand2 = fCandsTotal[iC2]; LmvmKinePar pRec = LmvmKinePar::Create(&cand1, &cand2); - double w = 1.; + double w = 1.; if (cand1.IsMcSignal()) w = fW * MinvScale(cand1.fMassSig); - else if (cand2.IsMcSignal()) w = fW * MinvScale(cand2.fMassSig); + else if (cand2.IsMcSignal()) + w = fW * MinvScale(cand2.fMassSig); bool isSameEvent = (cand1.fEventNumber == cand2.fEventNumber); // PLUTO electrons from two different events have to be scaled with w^2! - if (!isSameEvent && cand1.IsMcSignal() && cand2.IsMcSignal()) w = fW * fW * MinvScale(cand1.fMassSig) * MinvScale(cand2.fMassSig); + if (!isSameEvent && cand1.IsMcSignal() && cand2.IsMcSignal()) + w = fW * fW * MinvScale(cand1.fMassSig) * MinvScale(cand2.fMassSig); for (auto step : fH.fAnaSteps) { if (cand1.IsCutTill(step) && cand2.IsCutTill(step)) { @@ -1054,56 +1101,71 @@ void LmvmTask::CombinatorialPairs() if (cand1.fCharge * cand2.fCharge < 0) { if (isSameEvent) { fH.FillH1("hMinvCombPM_sameEv", step, pRec.fMinv, w); - cout << "LmvmTask::CombinatorialPairs(): e+e- same event (all). w = " << w << endl; //TODO: delete this line + cout << "LmvmTask::CombinatorialPairs(): e+e- same event (all). w = " << w + << endl; //TODO: delete this line } - else fH.FillH1("hMinvCombPM_mixedEv", step, pRec.fMinv, w); + else + fH.FillH1("hMinvCombPM_mixedEv", step, pRec.fMinv, w); } if (cand1.fCharge < 0 && cand2.fCharge < 0) { if (isSameEvent) { fH.FillH1("hMinvCombMM_sameEv", step, pRec.fMinv, w); - cout << "LmvmTask::CombinatorialPairs(): e-e- same event (all). w = " << w << endl; //TODO: delete this line + cout << "LmvmTask::CombinatorialPairs(): e-e- same event (all). w = " << w + << endl; //TODO: delete this line } - else fH.FillH1("hMinvCombMM_mixedEv", step, pRec.fMinv, w); + else + fH.FillH1("hMinvCombMM_mixedEv", step, pRec.fMinv, w); } if (cand1.fCharge > 0 && cand2.fCharge > 0) { if (isSameEvent) { fH.FillH1("hMinvCombPP_sameEv", step, pRec.fMinv, w); - cout << "LmvmTask::CombinatorialPairs(): e+e+ same event (all). w = " << w << endl; //TODO: delete this line - } - else fH.FillH1("hMinvCombPP_mixedEv", step, pRec.fMinv, w); + cout << "LmvmTask::CombinatorialPairs(): e+e+ same event (all). w = " << w + << endl; //TODO: delete this line + } + else + fH.FillH1("hMinvCombPP_mixedEv", step, pRec.fMinv, w); } - + // only UrQMD particles (also misidentified) if (!cand1.IsMcSignal() && !cand2.IsMcSignal()) { if (isSameEvent) { if (cand1.fCharge * cand2.fCharge < 0) fH.FillH1("hMinvCombPM_urqmdAll_sameEv", step, pRec.fMinv, w); - else if (cand1.fCharge < 0 && cand2.fCharge < 0) fH.FillH1("hMinvCombMM_urqmdAll_sameEv", step, pRec.fMinv, w); - else if (cand1.fCharge > 0 && cand2.fCharge > 0) fH.FillH1("hMinvCombPP_urqmdAll_sameEv", step, pRec.fMinv, w); + else if (cand1.fCharge < 0 && cand2.fCharge < 0) + fH.FillH1("hMinvCombMM_urqmdAll_sameEv", step, pRec.fMinv, w); + else if (cand1.fCharge > 0 && cand2.fCharge > 0) + fH.FillH1("hMinvCombPP_urqmdAll_sameEv", step, pRec.fMinv, w); } else { if (cand1.fCharge * cand2.fCharge < 0) fH.FillH1("hMinvCombPM_urqmdAll_mixedEv", step, pRec.fMinv, w); - else if (cand1.fCharge < 0 && cand2.fCharge < 0) fH.FillH1("hMinvCombMM_urqmdAll_mixedEv", step, pRec.fMinv, w); - else if (cand1.fCharge > 0 && cand2.fCharge > 0) fH.FillH1("hMinvCombPP_urqmdAll_mixedEv", step, pRec.fMinv, w); + else if (cand1.fCharge < 0 && cand2.fCharge < 0) + fH.FillH1("hMinvCombMM_urqmdAll_mixedEv", step, pRec.fMinv, w); + else if (cand1.fCharge > 0 && cand2.fCharge > 0) + fH.FillH1("hMinvCombPP_urqmdAll_mixedEv", step, pRec.fMinv, w); } } // only UrQMD electrons - if (!cand1.IsMcSignal() && !cand2.IsMcSignal() && std::abs(cand1.fMcPdg) == 11 && std::abs(cand2.fMcPdg) == 11) { // TODO: delete double brackets + if (!cand1.IsMcSignal() && !cand2.IsMcSignal() && std::abs(cand1.fMcPdg) == 11 + && std::abs(cand2.fMcPdg) == 11) { // TODO: delete double brackets if (isSameEvent) { if (cand1.fCharge * cand2.fCharge < 0) fH.FillH1("hMinvCombPM_urqmdEl_sameEv", step, pRec.fMinv, w); - else if (cand1.fCharge < 0 && cand2.fCharge < 0) fH.FillH1("hMinvCombMM_urqmdEl_sameEv", step, pRec.fMinv, w); - else if (cand1.fCharge > 0 && cand2.fCharge > 0) fH.FillH1("hMinvCombPP_urqmdEl_sameEv", step, pRec.fMinv, w); + else if (cand1.fCharge < 0 && cand2.fCharge < 0) + fH.FillH1("hMinvCombMM_urqmdEl_sameEv", step, pRec.fMinv, w); + else if (cand1.fCharge > 0 && cand2.fCharge > 0) + fH.FillH1("hMinvCombPP_urqmdEl_sameEv", step, pRec.fMinv, w); } else { if (cand1.fCharge * cand2.fCharge < 0) fH.FillH1("hMinvCombPM_urqmdEl_mixedEv", step, pRec.fMinv, w); - else if (cand1.fCharge < 0 && cand2.fCharge < 0) fH.FillH1("hMinvCombMM_urqmdEl_mixedEv", step, pRec.fMinv, w); - else if (cand1.fCharge > 0 && cand2.fCharge > 0) fH.FillH1("hMinvCombPP_urqmdEl_mixedEv", step, pRec.fMinv, w); + else if (cand1.fCharge < 0 && cand2.fCharge < 0) + fH.FillH1("hMinvCombMM_urqmdEl_mixedEv", step, pRec.fMinv, w); + else if (cand1.fCharge > 0 && cand2.fCharge > 0) + fH.FillH1("hMinvCombPP_urqmdEl_mixedEv", step, pRec.fMinv, w); } } - } // isCutTillStep - } // steps - } // cand2 - } // cand1 + } // isCutTillStep + } // steps + } // cand2 + } // cand1 } void LmvmTask::AssignMcToCands(vector<LmvmCand>& cands) @@ -1129,7 +1191,7 @@ void LmvmTask::AssignMcToCands(vector<LmvmCand>& cands) CbmMCTrack* mct2 = static_cast<CbmMCTrack*>(fMCTracks->At(iMc2)); if (mct2->GetMotherId() == cand.fMcMotherId && iMc2 != cand.fStsMcTrackId) { LmvmKinePar pKin = LmvmKinePar::Create(mct, mct2); - cand.fMassSig = pKin.fMinv; + cand.fMassSig = pKin.fMinv; } } } @@ -1181,7 +1243,7 @@ void LmvmTask::AssignMcToTopologyCands(vector<LmvmCand>& topoCands) void LmvmTask::PairSource(const LmvmCand& candP, const LmvmCand& candM, ELmvmAnaStep step, const LmvmKinePar& parRec) { ELmvmSrc src = LmvmUtils::GetMcPairSrc(candP, candM); - double w = (candP.IsMcSignal() || candM.IsMcSignal()) ? fW : 1.; + double w = (candP.IsMcSignal() || candM.IsMcSignal()) ? fW : 1.; fH.FillH1("hAnglePair", src, step, parRec.fAngle, w); if (src == ELmvmSrc::Bg) { @@ -1364,7 +1426,7 @@ void LmvmTask::FillPairHists(const LmvmCand& candP, const LmvmCand& candM, const fH.FillH1("hMinv", src, step, parRec.fMinv, w); if (!candP.IsMcSignal() && !candM.IsMcSignal()) fH.FillH1("hMinv_urqmdAll", src, step, parRec.fMinv, w); if (!candP.IsMcSignal() && !candM.IsMcSignal() && std::abs(candP.fMcPdg) == 11 && std::abs(candM.fMcPdg) == 11) - fH.FillH1("hMinv_urqmdEl", src, step, parRec.fMinv, w); + fH.FillH1("hMinv_urqmdEl", src, step, parRec.fMinv, w); fH.FillH2("hMinvPt", src, step, parRec.fMinv, parMc.fPt, w); @@ -1394,45 +1456,66 @@ string LmvmTask::GetPidString(const CbmMCTrack* mct, const LmvmCand* cand) return 0; } - int pdg = (mct != nullptr) ? mct->GetPdgCode() : cand->fMcPdg; - bool isMcSignal = (mct != nullptr) ? (mct->GetGeantProcessId() == kPPrimary && std::abs(mct->GetPdgCode()) == 11) : cand->IsMcSignal(); + int pdg = (mct != nullptr) ? mct->GetPdgCode() : cand->fMcPdg; + bool isMcSignal = (mct != nullptr) ? (mct->GetGeantProcessId() == kPPrimary && std::abs(mct->GetPdgCode()) == 11) + : cand->IsMcSignal(); string pidString = fH.fCandNames[9]; - if (isMcSignal && pdg == -11) pidString = fH.fCandNames[0]; - else if (isMcSignal && pdg == 11) pidString = fH.fCandNames[1]; - else if (!isMcSignal && pdg == -11) pidString = fH.fCandNames[2]; - else if (!isMcSignal && pdg == 11) pidString = fH.fCandNames[3]; - else if (pdg == 211) pidString = fH.fCandNames[4]; - else if (pdg == -211) pidString = fH.fCandNames[5]; - else if (pdg == 2212) pidString = fH.fCandNames[6]; - else if (pdg == 321) pidString = fH.fCandNames[7]; - else if (pdg == -321) pidString = fH.fCandNames[8]; - + if (isMcSignal && pdg == -11) pidString = fH.fCandNames[0]; + else if (isMcSignal && pdg == 11) + pidString = fH.fCandNames[1]; + else if (!isMcSignal && pdg == -11) + pidString = fH.fCandNames[2]; + else if (!isMcSignal && pdg == 11) + pidString = fH.fCandNames[3]; + else if (pdg == 211) + pidString = fH.fCandNames[4]; + else if (pdg == -211) + pidString = fH.fCandNames[5]; + else if (pdg == 2212) + pidString = fH.fCandNames[6]; + else if (pdg == 321) + pidString = fH.fCandNames[7]; + else if (pdg == -321) + pidString = fH.fCandNames[8]; + return pidString; } string LmvmTask::GetPidString(double vertexMag, int pdg) { string pidString = fH.fGTrackNames[14]; - bool isPrim = IsPrimary(vertexMag); + bool isPrim = IsPrimary(vertexMag); if (isPrim) { - if (pdg == -11) pidString = fH.fGTrackNames[1]; - else if (pdg == 11) pidString = fH.fGTrackNames[3]; - else if (pdg == 211) pidString = fH.fGTrackNames[5]; - else if (pdg == -211) pidString = fH.fGTrackNames[7]; - else if (pdg == 2212) pidString = fH.fGTrackNames[9]; - else if (pdg == 321) pidString = fH.fGTrackNames[11]; - else if (pdg == -321) pidString = fH.fGTrackNames[13]; + if (pdg == -11) pidString = fH.fGTrackNames[1]; + else if (pdg == 11) + pidString = fH.fGTrackNames[3]; + else if (pdg == 211) + pidString = fH.fGTrackNames[5]; + else if (pdg == -211) + pidString = fH.fGTrackNames[7]; + else if (pdg == 2212) + pidString = fH.fGTrackNames[9]; + else if (pdg == 321) + pidString = fH.fGTrackNames[11]; + else if (pdg == -321) + pidString = fH.fGTrackNames[13]; } else { - if (pdg == -11) pidString = fH.fGTrackNames[0]; - else if (pdg == 11) pidString = fH.fGTrackNames[2]; - else if (pdg == 211) pidString = fH.fGTrackNames[4]; - else if (pdg == -211) pidString = fH.fGTrackNames[6]; - else if (pdg == 2212) pidString = fH.fGTrackNames[8]; - else if (pdg == 321) pidString = fH.fGTrackNames[10]; - else if (pdg == -321) pidString = fH.fGTrackNames[12]; + if (pdg == -11) pidString = fH.fGTrackNames[0]; + else if (pdg == 11) + pidString = fH.fGTrackNames[2]; + else if (pdg == 211) + pidString = fH.fGTrackNames[4]; + else if (pdg == -211) + pidString = fH.fGTrackNames[6]; + else if (pdg == 2212) + pidString = fH.fGTrackNames[8]; + else if (pdg == 321) + pidString = fH.fGTrackNames[10]; + else if (pdg == -321) + pidString = fH.fGTrackNames[12]; } return pidString; @@ -1451,17 +1534,17 @@ void LmvmTask::AnalyseCandidates() // single candidates for (const auto& cand : fCands) { - double w = (cand.IsMcSignal()) ? fW : 1.; + double w = (cand.IsMcSignal()) ? fW : 1.; CbmMCTrack* mcTrack = nullptr; if (cand.fStsMcTrackId >= 0) mcTrack = static_cast<CbmMCTrack*>(fMCTracks->At(cand.fStsMcTrackId)); - int pdg = mcTrack->GetPdgCode(); - double mom = mcTrack->GetP(); + int pdg = mcTrack->GetPdgCode(); + double mom = mcTrack->GetP(); string pidString = GetPidString(nullptr, &cand); for (auto step : fH.fAnaSteps) { if (cand.IsCutTill(step)) { TrackSource(cand, step, pdg); - FillCandPidValues(mcTrack, cand, step); // if std::abs(pdg) != 11, this is misidentification + FillCandPidValues(mcTrack, cand, step); // if std::abs(pdg) != 11, this is misidentification if (mcTrack != nullptr) CheckTofId(mcTrack, cand, step, pdg); if (step >= ELmvmAnaStep::ElId && std::abs(pdg) == 211) fH.FillH1("hPionSupp_idEl", step, mom, w); @@ -1472,22 +1555,22 @@ void LmvmTask::AnalyseCandidates() } if (cand.IsCutTill(ELmvmAnaStep::ElId)) { - fH.FillH2("hChi2VsMom_sts_" + pidString, cand.fMomentum.Mag(), cand.fChi2Sts, w); + fH.FillH2("hChi2VsMom_sts_" + pidString, cand.fMomentum.Mag(), cand.fChi2Sts, w); fH.FillH2("hChi2VsMom_rich_" + pidString, cand.fMomentum.Mag(), cand.fChi2Rich, w); - fH.FillH2("hChi2VsMom_trd_" + pidString, cand.fMomentum.Mag(), cand.fChi2Trd, w); + fH.FillH2("hChi2VsMom_trd_" + pidString, cand.fMomentum.Mag(), cand.fChi2Trd, w); - fH.FillH2("hTofTimeVsChi2_sts_" + pidString, cand.fChi2Sts, cand.fTime, w); + fH.FillH2("hTofTimeVsChi2_sts_" + pidString, cand.fChi2Sts, cand.fTime, w); fH.FillH2("hTofTimeVsChi2_rich_" + pidString, cand.fChi2Rich, cand.fTime, w); - fH.FillH2("hTofTimeVsChi2_trd_" + pidString, cand.fChi2Trd, cand.fTime, w); + fH.FillH2("hTofTimeVsChi2_trd_" + pidString, cand.fChi2Trd, cand.fTime, w); fH.FillH2("hChi2Comb_StsRich_" + pidString, cand.fChi2Sts, cand.fChi2Rich, w); - fH.FillH2("hChi2Comb_StsTrd_" + pidString, cand.fChi2Sts, cand.fChi2Trd, w); + fH.FillH2("hChi2Comb_StsTrd_" + pidString, cand.fChi2Sts, cand.fChi2Trd, w); fH.FillH2("hChi2Comb_RichTrd_" + pidString, cand.fChi2Rich, cand.fChi2Trd, w); fH.FillH2("hTofTimeVsMom_cands_" + pidString, cand.fMomentum.Mag(), cand.fTime, w); fH.FillH2("hTofHitTrackDist_cands_" + pidString, cand.fMomentum.Mag(), cand.fTofDist, w); } - + DifferenceSignalAndBg(cand); // beta-momentum spectrum @@ -1496,9 +1579,9 @@ void LmvmTask::AnalyseCandidates() double m = cand.fMass; double r = p / m; double beta = r / sqrt(1. + r * r); - double q = (cand.fCharge > 0) ? 1 : (cand.fCharge < 0) ? -1. : 0.; + double q = (cand.fCharge > 0) ? 1 : (cand.fCharge < 0) ? -1. : 0.; fH.FillH2("hBetaMom_cands_" + ptcl, q * p, beta); - } // single candidates + } // single candidates // candidate pairs for (const auto& candP : fCands) { @@ -1526,18 +1609,27 @@ void LmvmTask::CheckTofId(const CbmMCTrack* mct, const LmvmCand& cand, ELmvmAnaS TVector3 v; mct->GetStartVertex(v); bool isPrim = IsPrimary(v.Mag()); - double pt = mct->GetPt(); - double rap = mct->GetRapidity(); + double pt = mct->GetPt(); + double rap = mct->GetRapidity(); // check PIDs in "Tof pile" if (cand.fMomentum.Mag() > 0.3 && cand.fMomentum.Mag() < 1. && cand.fMass2 > -0.012 && cand.fMass2 < 0.01) { string pidString = GetPidString(nullptr, &cand); - double pdgBin = (pdg == 11 && cand.IsMcSignal()) ? 0.5 : (pdg == -11 && cand.IsMcSignal()) ? 1.5 - : (pdg == 11 && !cand.IsMcSignal() && isPrim) ? 2.5 : (pdg == -11 && !cand.IsMcSignal() && isPrim) ? 3.5 - : (pdg == 11 && !cand.IsMcSignal() && !isPrim) ? 4.5 : (pdg == -11 && !cand.IsMcSignal() && !isPrim) ? 5.5 - : (pdg == 211) ? 6.5 : (pdg == -211) ? 7.5 : (pdg == 2212) ? 8.5 : (pdg == 321) ? 9.5 : (pdg == -321) ? 10.5 : 11.5; + double pdgBin = (pdg == 11 && cand.IsMcSignal()) ? 0.5 + : (pdg == -11 && cand.IsMcSignal()) ? 1.5 + : (pdg == 11 && !cand.IsMcSignal() && isPrim) ? 2.5 + : (pdg == -11 && !cand.IsMcSignal() && isPrim) ? 3.5 + : (pdg == 11 && !cand.IsMcSignal() && !isPrim) ? 4.5 + : (pdg == -11 && !cand.IsMcSignal() && !isPrim) ? 5.5 + : (pdg == 211) ? 6.5 + : (pdg == -211) ? 7.5 + : (pdg == 2212) ? 8.5 + : (pdg == 321) ? 9.5 + : (pdg == -321) ? 10.5 + : 11.5; fH.FillH1("hTofPilePdgs_cands", step, pdgBin); - if (step == ELmvmAnaStep::ElId && cand.IsCutTill(step)) fH.FillH2("hTofPilePty_cands_" + pidString, rap, pt); // TODO: check newly (25.10.22) added '&& IsTill(step)'!! + if (step == ELmvmAnaStep::ElId && cand.IsCutTill(step)) + fH.FillH2("hTofPilePty_cands_" + pidString, rap, pt); // TODO: check newly (25.10.22) added '&& IsTill(step)'!! } // check vertex of misidentified particles in ToF after electron ID // TODO: split this up into single contributions? @@ -1549,8 +1641,9 @@ void LmvmTask::CheckTofId(const CbmMCTrack* mct, const LmvmCand& cand, ELmvmAnaS } } -void LmvmTask::FillCandPidValues(const CbmMCTrack* mct, const LmvmCand& cand, ELmvmAnaStep step) { - double w = (cand.IsMcSignal()) ? fW : 1.; +void LmvmTask::FillCandPidValues(const CbmMCTrack* mct, const LmvmCand& cand, ELmvmAnaStep step) +{ + double w = (cand.IsMcSignal()) ? fW : 1.; string pidString = GetPidString(nullptr, &cand); fH.FillH1("hMom_cands_" + pidString, step, cand.fMomentum.Mag(), w); @@ -1558,7 +1651,7 @@ void LmvmTask::FillCandPidValues(const CbmMCTrack* mct, const LmvmCand& cand, EL fH.FillH2("hTofM2_cands_" + pidString, step, cand.fMomentum.Mag(), cand.fMass2, w); if (mct != nullptr) { - double rat = cand.fMomentum.Mag() / mct->GetP(); + double rat = cand.fMomentum.Mag() / mct->GetP(); string ptcl = GetPidString(nullptr, &cand); fH.FillH1("hMomRatio_cands_" + pidString, step, rat); fH.FillH2("hMomRatioVsMom_cands_" + ptcl, step, mct->GetP(), rat); @@ -1757,7 +1850,7 @@ void LmvmTask::DifferenceSignalAndBg(const LmvmCand& cand) } double mvd1r = sqrt(mvd1x * mvd1x + mvd1y * mvd1y); double mvd2r = sqrt(mvd2x * mvd2x + mvd2y * mvd2y); - fH.FillH1("hNofMvdHits", cand.fMcSrc, stsTrack->GetNofMvdHits(), fW); + fH.FillH1("hNofMvdHits", cand.fMcSrc, stsTrack->GetNofMvdHits(), fW); fH.FillH2("hMvdXY_1", cand.fMcSrc, mvd1x, mvd1y, fW); fH.FillH1("hMvdR_1", cand.fMcSrc, mvd1r, fW); fH.FillH2("hMvdXY_2", cand.fMcSrc, mvd2x, mvd2y, fW); diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmTask.h b/analysis/PWGDIL/dielectron/lmvm/LmvmTask.h index 8e8f8fcfbf6563d181c5fd4860720c40f1846815..263594b0951bc82143c3185acbec7545c861e9fe 100755 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmTask.h +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmTask.h @@ -5,6 +5,7 @@ #ifndef LMVM_TASK_H #define LMVM_TASK_H +#include "CbmGlobalTrack.h" #include "CbmKFVertex.h" #include "CbmStsKFTrackFitter.h" @@ -12,8 +13,6 @@ #include "FairRootManager.h" #include "FairTask.h" -#include "CbmGlobalTrack.h" - #include <fstream> #include <map> #include <string> @@ -123,7 +122,8 @@ public: void AnalyseGlobalTracks(); void CheckMismatches(const CbmGlobalTrack* gTrack, int pdg, bool isElectron, const std::string& ptcl, double weight); void BetaMom(const CbmMCTrack* mct, const CbmGlobalTrack* gTrack, const std::string& ptcl); - void CheckTofIdentification(const CbmGlobalTrack* gTrack, const std::string& pidString, double mom, double m2, int pdg, bool isTofEl); + void CheckTofIdentification(const CbmGlobalTrack* gTrack, const std::string& pidString, double mom, double m2, + int pdg, bool isTofEl); void PidVsMom(const CbmGlobalTrack* gTrack, int iGTrack, int pdg, double mom); void CheckGammaConvAndPi0(); @@ -217,7 +217,7 @@ private: std::map<int, int> fNofHitsInRingMap; // Number of hits in the MC RICH ring - double fZ = -44.; // z-position of target + double fZ = -44.; // z-position of target std::string fParticle = ""; diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmUtils.cxx b/analysis/PWGDIL/dielectron/lmvm/LmvmUtils.cxx index 149a5334bcfbb49d32e45811dfc20113e7ccdc79..b7f48d30b8beccf1617cae6c2ef95160536240de 100755 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmUtils.cxx +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmUtils.cxx @@ -4,14 +4,14 @@ #include "LmvmUtils.h" -#include "Logger.h" - #include "CbmKFVertex.h" #include "CbmL1PFFitter.h" #include "CbmMCTrack.h" #include "CbmStsTrack.h" #include "cbm/elid/CbmLitGlobalElectronId.h" +#include "Logger.h" + #include "TClonesArray.h" #include "TDatabasePDG.h" #include "TMCProcess.h" @@ -260,19 +260,23 @@ void LmvmUtils::IsElectron(int globalTrackIndex, double momentum, double momentu if (cand->fTofDist > momentum*slope + b) tofEl = false; }*/ - bool isMomCut = (momentumCut > 0.) ? (momentum < momentumCut) : true; + bool isMomCut = (momentumCut > 0.) ? (momentum < momentumCut) : true; cand->fIsElectron = (richEl && trdEl && tofEl && isMomCut); } -void LmvmUtils::IsRichElectron(int globalTrackIndex, double momentum, LmvmCand* cand) { +void LmvmUtils::IsRichElectron(int globalTrackIndex, double momentum, LmvmCand* cand) +{ cand->fIsRichElectron = CbmLitGlobalElectronId::GetInstance().IsRichElectron(globalTrackIndex, momentum); } -void LmvmUtils::IsTrdElectron(int globalTrackIndex, double momentum, LmvmCand* cand) { - cand->fIsTrdElectron = (momentum < 1.) ? true : CbmLitGlobalElectronId::GetInstance().IsTrdElectron(globalTrackIndex, momentum); +void LmvmUtils::IsTrdElectron(int globalTrackIndex, double momentum, LmvmCand* cand) +{ + cand->fIsTrdElectron = + (momentum < 1.) ? true : CbmLitGlobalElectronId::GetInstance().IsTrdElectron(globalTrackIndex, momentum); } -void LmvmUtils::IsTofElectron(int globalTrackIndex, double momentum, LmvmCand* cand) { +void LmvmUtils::IsTofElectron(int globalTrackIndex, double momentum, LmvmCand* cand) +{ cand->fIsTofElectron = CbmLitGlobalElectronId::GetInstance().IsTofElectron(globalTrackIndex, momentum); } @@ -303,18 +307,18 @@ string LmvmUtils::GetChargeStr(const CbmMCTrack* mct) double LmvmUtils::GetMassScaleInmed(double minv) // TODO: make these more elegant and delete cout's { - int nArray = sizeof(fMinvArray); - double weight = -1.; + int nArray = sizeof(fMinvArray); + double weight = -1.; if (minv < fMinvArray[0]) return fScaleArrayInmed[0]; else { for (int i = 1; i < nArray; i++) { if (fMinvArray[i] > minv) { - double dy = fScaleArrayInmed[i] - fScaleArrayInmed[i-1]; - double dx = fMinvArray[i] - fMinvArray[i-1]; - double slope = dy/dx; - double dLeft = minv - fMinvArray[i-1]; - weight = fScaleArrayInmed[i-1] + slope * dLeft; + double dy = fScaleArrayInmed[i] - fScaleArrayInmed[i - 1]; + double dx = fMinvArray[i] - fMinvArray[i - 1]; + double slope = dy / dx; + double dLeft = minv - fMinvArray[i - 1]; + weight = fScaleArrayInmed[i - 1] + slope * dLeft; return weight; } } @@ -324,18 +328,18 @@ double LmvmUtils::GetMassScaleInmed(double minv) // TODO: make these more elega double LmvmUtils::GetMassScaleQgp(double minv) // TODO: make these more elegant and delete cout's { - int nArray = sizeof(fMinvArray); - double weight = -1.; + int nArray = sizeof(fMinvArray); + double weight = -1.; if (minv < fMinvArray[0]) return fScaleArrayQgp[0]; else { for (int i = 1; i < nArray; i++) { if (fMinvArray[i] > minv) { - double dy = fScaleArrayQgp[i] - fScaleArrayQgp[i-1]; - double dx = fMinvArray[i] - fMinvArray[i-1]; - double slope = dy/dx; - double dLeft = minv - fMinvArray[i-1]; - weight = fScaleArrayQgp[i-1] + slope * dLeft; + double dy = fScaleArrayQgp[i] - fScaleArrayQgp[i - 1]; + double dx = fMinvArray[i] - fMinvArray[i - 1]; + double slope = dy / dx; + double dLeft = minv - fMinvArray[i - 1]; + weight = fScaleArrayQgp[i - 1] + slope * dLeft; return weight; } } diff --git a/analysis/PWGDIL/dielectron/lmvm/LmvmUtils.h b/analysis/PWGDIL/dielectron/lmvm/LmvmUtils.h index cd04c23d5272bc01da9c8f8b08d74052957e21db..7374453fa05fc199d02b96d3334d9cc39b2c8d36 100755 --- a/analysis/PWGDIL/dielectron/lmvm/LmvmUtils.h +++ b/analysis/PWGDIL/dielectron/lmvm/LmvmUtils.h @@ -99,71 +99,69 @@ public: static std::string GetChargeStr(const LmvmCand* cand); static std::string GetChargeStr(const CbmMCTrack* mct); - + static double GetMassScaleInmed(double minv); - static double GetMassScaleQgp(double minv); - - ClassDef(LmvmUtils, 1); + static double GetMassScaleQgp(double minv); -private: - - static constexpr double fMinvArray[170] = { - 0.0195, 0.0395, 0.0595, 0.0795, 0.0995, 0.1195, 0.1395, 0.1595, 0.1795, 0.1995, 0.2195, 0.2395, 0.2595, 0.2795, - 0.2995, 0.3195, 0.3395, 0.3595, 0.3795, 0.3995, 0.4195, 0.4395, 0.4595, 0.4795, 0.4995, 0.5195, 0.5395, 0.5595, - 0.5795, 0.5995, 0.6195, 0.6395, 0.6595, 0.6795, 0.6995, 0.7195, 0.7395, 0.7595, 0.7795, 0.7995, 0.8195, 0.8395, - 0.8595, 0.8795, 0.8995, 0.9195, 0.9395, 0.9595, 0.9795, 0.9995, 1.0195, 1.0395, 1.0595, 1.0795, 1.0995, 1.1195, - 1.1395, 1.1595, 1.1795, 1.1995, 1.2195, 1.2395, 1.2595, 1.2795, 1.2995, 1.3195, 1.3395, 1.3595, 1.3795, 1.3995, - 1.4195, 1.4395, 1.4595, 1.4795, 1.4995, 1.5195, 1.5395, 1.5595, 1.5795, 1.5995, 1.6195, 1.6395, 1.6595, 1.6795, - 1.6995, 1.7195, 1.7395, 1.7595, 1.7795, 1.7995, 1.8195, 1.8395, 1.8595, 1.8795, 1.8995, 1.9195, 1.9395, 1.9595, - 1.9795, 1.9995, 2.0195, 2.0395, 2.0595, 2.0795, 2.0995, 2.1195, 2.1395, 2.1595, 2.1795, 2.1995, 2.2195, 2.2395, - 2.2595, 2.2795, 2.2995, 2.3195, 2.3395, 2.3595, 2.3795, 2.3995, 2.4195, 2.4395, 2.4595, 2.4795, 2.4995, 2.5195, - 2.5395, 2.5595, 2.5795, 2.5995, 2.6195, 2.6395, 2.6595, 2.6795, 2.6995, 2.7195, 2.7395, 2.7595, 2.7795, 2.7995, - 2.8195, 2.8395, 2.8595, 2.8795, 2.8995, 2.9195, 2.9395, 2.9595, 2.9795, 2.9995, 3.0195, 3.0395, 3.0595, 3.0795, - 3.0995, 3.1195, 3.1395, 3.1595, 3.1795, 3.1995, 3.2195, 3.2395, 3.2595, 3.2795, 3.2995, 3.3195, 3.3395, 3.3595, - 3.3795, 3.3995}; - - static constexpr double fScaleArrayInmed[170] = { // for 12 AGeV - 41.706, 18.918, 11.465, 8.4388, 5.9176, 4.9025, 3.8087, 3.0387, 2.5856, - 2.1142, 1.7603, 1.5327, 1.28, 1.1579, 1.0367, 0.89355, 0.81317, 0.71582, - 0.65863, 0.59678, 0.53702, 0.45378, 0.41238, 0.37502, 0.33593, 0.28791, 0.26352, - 0.23939, 0.21167, 0.19479, 0.19204, 0.17492, 0.15811, 0.15479, 0.14935, 0.13803, - 0.1354, 0.11993, 0.1046, 0.08226, 0.073183, 0.055433, 0.043467, 0.033975, 0.028025, - 0.021504, 0.016863, 0.014108, 0.01094, 0.0088095, 0.007324, 0.0057162, 0.0046817, 0.0037459, - 0.0030017, 0.0024459, 0.0020671, 0.0016089, 0.0013754, 0.0011223, 0.00096256, 0.00081647, 0.00072656, - 0.00060776, 0.00051243, 0.00045705, 0.00039636, 0.00036259, 0.00033248, 0.0002953, 0.00027328, 0.00023776, - 0.00022163, 0.00019852, 0.000186, 0.00016846, 0.00015469, 0.00014169, 0.00013343, 0.00011594, 0.00010722, - 0.00010205, 9.1907e-05, 8.3718e-05, 7.5457e-05, 6.7192e-05, 6.2202e-05, 5.7372e-05, 4.8314e-05, 4.5502e-05, - 4.1334e-05, 3.7429e-05, 3.2131e-05, 3.0103e-05, 2.6125e-05, 2.3601e-05, 2.1167e-05, 1.94e-05, 1.7025e-05, - 1.5496e-05, 1.3704e-05, 1.1866e-05, 1.1135e-05, 9.8842e-06, 8.9101e-06, 7.9225e-06, 7.0706e-06, 6.3536e-06, - 5.3786e-06, 4.7179e-06, 4.2128e-06, 4.0015e-06, 3.4118e-06, 3.1864e-06, 2.734e-06, 2.3844e-06, 2.173e-06, - 1.8774e-06, 1.6468e-06, 1.501e-06, 1.3597e-06, 1.2113e-06, 1.0384e-06, 9.4105e-07, 8.4223e-07, 7.434e-07, - 6.5049e-07, 5.8824e-07, 5.3603e-07, 4.6756e-07, 4.1173e-07, 3.5872e-07, 3.2764e-07, 2.9889e-07, 2.5989e-07, - 2.219e-07, 1.9468e-07, 1.816e-07, 1.5707e-07, 1.3565e-07, 1.2619e-07, 1.0919e-07, 1.0071e-07, 8.4632e-08, - 7.6459e-08, 6.829e-08, 6.2046e-08, 5.5335e-08, 4.5937e-08, 4.2426e-08, 3.567e-08, 3.4051e-08, 2.9627e-08, - 2.5249e-08, 2.2767e-08, 2.1054e-08, 1.7873e-08, 1.574e-08, 1.3713e-08, 1.23e-08, 1.1045e-08, 9.5536e-09, - 8.5859e-09, 7.7217e-09, 6.9958e-09, 6.0992e-09, 5.3453e-09, 4.7659e-09, 4.3313e-09, 3.6575e-09}; - - static constexpr double fScaleArrayQgp[170] = { // for 12 AGeV - 39.496, 17.961, 11.024, 8.2093, 5.8331, 4.8995, 3.8612, 3.1258, 2.7006, - 2.2465, 1.908, 1.699, 1.4435, 1.3253, 1.2059, 1.049, 0.96753, 0.86685, - 0.81407, 0.75959, 0.70663, 0.61951, 0.58586, 0.55534, 0.51902, 0.46377, 0.4415, - 0.41412, 0.37414, 0.34883, 0.34494, 0.31141, 0.2762, 0.26331, 0.24693, 0.22286, - 0.21697, 0.1972, 0.1841, 0.16097, 0.16352, 0.14345, 0.13096, 0.11911, 0.11399, - 0.10111, 0.0913, 0.08764, 0.077745, 0.071417, 0.067561, 0.05987, 0.055543, 0.050193, - 0.045244, 0.04128, 0.03898, 0.03365, 0.031622, 0.028217, 0.026215, 0.023919, 0.022648, - 0.019915, 0.017524, 0.016145, 0.014357, 0.013362, 0.012368, 0.011036, 0.010198, 0.0088275, - 0.0081762, 0.0072697, 0.00675, 0.0060424, 0.0054788, 0.0049588, 0.0046174, 0.0039685, 0.00363, - 0.0034204, 0.0030534, 0.0027606, 0.0024723, 0.0021893, 0.0020174, 0.0018545, 0.0015584, 0.0014661, - 0.0013315, 0.0012065, 0.0010375, 0.00097456, 0.00084865, 0.00076982, 0.00069371, 0.00063931, 0.00056442, - 0.00051712, 0.00046054, 0.00040174, 0.00037996, 0.00034009, 0.00030921, 0.00027738, 0.00024981, 0.00022659, - 0.00019366, 0.00017153, 0.00015469, 0.00014841, 0.00012783, 0.00012061, 0.00010456, 9.2145e-05, 8.4856e-05, - 7.4087e-05, 6.5675e-05, 6.0496e-05, 5.5386e-05, 4.9865e-05, 4.3202e-05, 3.9571e-05, 3.5821e-05, 3.201e-05, - 2.8322e-05, 2.5886e-05, 2.384e-05, 2.1016e-05, 1.8703e-05, 1.6467e-05, 1.5199e-05, 1.4011e-05, 1.2311e-05, - 1.0621e-05, 9.4155e-06, 8.874e-06, 7.7548e-06, 6.7662e-06, 6.3589e-06, 5.5585e-06, 5.1791e-06, 4.3965e-06, - 4.012e-06, 3.6195e-06, 3.3215e-06, 2.9918e-06, 2.5084e-06, 2.3397e-06, 1.9865e-06, 1.915e-06, 1.6826e-06, - 1.448e-06, 1.3183e-06, 1.231e-06, 1.0551e-06, 9.3811e-07, 8.2511e-07, 7.4714e-07, 6.7735e-07, 5.9142e-07, - 5.3654e-07, 4.8709e-07, 4.4543e-07, 3.9199e-07, 3.4674e-07, 3.1203e-07, 2.862e-07, 2.4391e-07}; + ClassDef(LmvmUtils, 1); +private: + static constexpr double fMinvArray[170] = { + 0.0195, 0.0395, 0.0595, 0.0795, 0.0995, 0.1195, 0.1395, 0.1595, 0.1795, 0.1995, 0.2195, 0.2395, 0.2595, 0.2795, + 0.2995, 0.3195, 0.3395, 0.3595, 0.3795, 0.3995, 0.4195, 0.4395, 0.4595, 0.4795, 0.4995, 0.5195, 0.5395, 0.5595, + 0.5795, 0.5995, 0.6195, 0.6395, 0.6595, 0.6795, 0.6995, 0.7195, 0.7395, 0.7595, 0.7795, 0.7995, 0.8195, 0.8395, + 0.8595, 0.8795, 0.8995, 0.9195, 0.9395, 0.9595, 0.9795, 0.9995, 1.0195, 1.0395, 1.0595, 1.0795, 1.0995, 1.1195, + 1.1395, 1.1595, 1.1795, 1.1995, 1.2195, 1.2395, 1.2595, 1.2795, 1.2995, 1.3195, 1.3395, 1.3595, 1.3795, 1.3995, + 1.4195, 1.4395, 1.4595, 1.4795, 1.4995, 1.5195, 1.5395, 1.5595, 1.5795, 1.5995, 1.6195, 1.6395, 1.6595, 1.6795, + 1.6995, 1.7195, 1.7395, 1.7595, 1.7795, 1.7995, 1.8195, 1.8395, 1.8595, 1.8795, 1.8995, 1.9195, 1.9395, 1.9595, + 1.9795, 1.9995, 2.0195, 2.0395, 2.0595, 2.0795, 2.0995, 2.1195, 2.1395, 2.1595, 2.1795, 2.1995, 2.2195, 2.2395, + 2.2595, 2.2795, 2.2995, 2.3195, 2.3395, 2.3595, 2.3795, 2.3995, 2.4195, 2.4395, 2.4595, 2.4795, 2.4995, 2.5195, + 2.5395, 2.5595, 2.5795, 2.5995, 2.6195, 2.6395, 2.6595, 2.6795, 2.6995, 2.7195, 2.7395, 2.7595, 2.7795, 2.7995, + 2.8195, 2.8395, 2.8595, 2.8795, 2.8995, 2.9195, 2.9395, 2.9595, 2.9795, 2.9995, 3.0195, 3.0395, 3.0595, 3.0795, + 3.0995, 3.1195, 3.1395, 3.1595, 3.1795, 3.1995, 3.2195, 3.2395, 3.2595, 3.2795, 3.2995, 3.3195, 3.3395, 3.3595, + 3.3795, 3.3995}; + + static constexpr double fScaleArrayInmed[170] = { // for 12 AGeV + 41.706, 18.918, 11.465, 8.4388, 5.9176, 4.9025, 3.8087, 3.0387, 2.5856, + 2.1142, 1.7603, 1.5327, 1.28, 1.1579, 1.0367, 0.89355, 0.81317, 0.71582, + 0.65863, 0.59678, 0.53702, 0.45378, 0.41238, 0.37502, 0.33593, 0.28791, 0.26352, + 0.23939, 0.21167, 0.19479, 0.19204, 0.17492, 0.15811, 0.15479, 0.14935, 0.13803, + 0.1354, 0.11993, 0.1046, 0.08226, 0.073183, 0.055433, 0.043467, 0.033975, 0.028025, + 0.021504, 0.016863, 0.014108, 0.01094, 0.0088095, 0.007324, 0.0057162, 0.0046817, 0.0037459, + 0.0030017, 0.0024459, 0.0020671, 0.0016089, 0.0013754, 0.0011223, 0.00096256, 0.00081647, 0.00072656, + 0.00060776, 0.00051243, 0.00045705, 0.00039636, 0.00036259, 0.00033248, 0.0002953, 0.00027328, 0.00023776, + 0.00022163, 0.00019852, 0.000186, 0.00016846, 0.00015469, 0.00014169, 0.00013343, 0.00011594, 0.00010722, + 0.00010205, 9.1907e-05, 8.3718e-05, 7.5457e-05, 6.7192e-05, 6.2202e-05, 5.7372e-05, 4.8314e-05, 4.5502e-05, + 4.1334e-05, 3.7429e-05, 3.2131e-05, 3.0103e-05, 2.6125e-05, 2.3601e-05, 2.1167e-05, 1.94e-05, 1.7025e-05, + 1.5496e-05, 1.3704e-05, 1.1866e-05, 1.1135e-05, 9.8842e-06, 8.9101e-06, 7.9225e-06, 7.0706e-06, 6.3536e-06, + 5.3786e-06, 4.7179e-06, 4.2128e-06, 4.0015e-06, 3.4118e-06, 3.1864e-06, 2.734e-06, 2.3844e-06, 2.173e-06, + 1.8774e-06, 1.6468e-06, 1.501e-06, 1.3597e-06, 1.2113e-06, 1.0384e-06, 9.4105e-07, 8.4223e-07, 7.434e-07, + 6.5049e-07, 5.8824e-07, 5.3603e-07, 4.6756e-07, 4.1173e-07, 3.5872e-07, 3.2764e-07, 2.9889e-07, 2.5989e-07, + 2.219e-07, 1.9468e-07, 1.816e-07, 1.5707e-07, 1.3565e-07, 1.2619e-07, 1.0919e-07, 1.0071e-07, 8.4632e-08, + 7.6459e-08, 6.829e-08, 6.2046e-08, 5.5335e-08, 4.5937e-08, 4.2426e-08, 3.567e-08, 3.4051e-08, 2.9627e-08, + 2.5249e-08, 2.2767e-08, 2.1054e-08, 1.7873e-08, 1.574e-08, 1.3713e-08, 1.23e-08, 1.1045e-08, 9.5536e-09, + 8.5859e-09, 7.7217e-09, 6.9958e-09, 6.0992e-09, 5.3453e-09, 4.7659e-09, 4.3313e-09, 3.6575e-09}; + + static constexpr double fScaleArrayQgp[170] = { // for 12 AGeV + 39.496, 17.961, 11.024, 8.2093, 5.8331, 4.8995, 3.8612, 3.1258, 2.7006, + 2.2465, 1.908, 1.699, 1.4435, 1.3253, 1.2059, 1.049, 0.96753, 0.86685, + 0.81407, 0.75959, 0.70663, 0.61951, 0.58586, 0.55534, 0.51902, 0.46377, 0.4415, + 0.41412, 0.37414, 0.34883, 0.34494, 0.31141, 0.2762, 0.26331, 0.24693, 0.22286, + 0.21697, 0.1972, 0.1841, 0.16097, 0.16352, 0.14345, 0.13096, 0.11911, 0.11399, + 0.10111, 0.0913, 0.08764, 0.077745, 0.071417, 0.067561, 0.05987, 0.055543, 0.050193, + 0.045244, 0.04128, 0.03898, 0.03365, 0.031622, 0.028217, 0.026215, 0.023919, 0.022648, + 0.019915, 0.017524, 0.016145, 0.014357, 0.013362, 0.012368, 0.011036, 0.010198, 0.0088275, + 0.0081762, 0.0072697, 0.00675, 0.0060424, 0.0054788, 0.0049588, 0.0046174, 0.0039685, 0.00363, + 0.0034204, 0.0030534, 0.0027606, 0.0024723, 0.0021893, 0.0020174, 0.0018545, 0.0015584, 0.0014661, + 0.0013315, 0.0012065, 0.0010375, 0.00097456, 0.00084865, 0.00076982, 0.00069371, 0.00063931, 0.00056442, + 0.00051712, 0.00046054, 0.00040174, 0.00037996, 0.00034009, 0.00030921, 0.00027738, 0.00024981, 0.00022659, + 0.00019366, 0.00017153, 0.00015469, 0.00014841, 0.00012783, 0.00012061, 0.00010456, 9.2145e-05, 8.4856e-05, + 7.4087e-05, 6.5675e-05, 6.0496e-05, 5.5386e-05, 4.9865e-05, 4.3202e-05, 3.9571e-05, 3.5821e-05, 3.201e-05, + 2.8322e-05, 2.5886e-05, 2.384e-05, 2.1016e-05, 1.8703e-05, 1.6467e-05, 1.5199e-05, 1.4011e-05, 1.2311e-05, + 1.0621e-05, 9.4155e-06, 8.874e-06, 7.7548e-06, 6.7662e-06, 6.3589e-06, 5.5585e-06, 5.1791e-06, 4.3965e-06, + 4.012e-06, 3.6195e-06, 3.3215e-06, 2.9918e-06, 2.5084e-06, 2.3397e-06, 1.9865e-06, 1.915e-06, 1.6826e-06, + 1.448e-06, 1.3183e-06, 1.231e-06, 1.0551e-06, 9.3811e-07, 8.2511e-07, 7.4714e-07, 6.7735e-07, 5.9142e-07, + 5.3654e-07, 4.8709e-07, 4.4543e-07, 3.9199e-07, 3.4674e-07, 3.1203e-07, 2.862e-07, 2.4391e-07}; }; #endif