diff --git a/analysis/PWGDIL/dielectron/papaframework/scripts/Config_dilepton_testing.C b/analysis/PWGDIL/dielectron/papaframework/scripts/Config_dilepton_testing.C new file mode 100644 index 0000000000000000000000000000000000000000..fe6e83c6854871aed7719f40abfbd650adace57c --- /dev/null +++ b/analysis/PWGDIL/dielectron/papaframework/scripts/Config_dilepton_testing.C @@ -0,0 +1,1333 @@ +/// \file Config_dilepton_testing.C +// \brief A template task configuration macro with example and explanations +/// +/// PairAnalysis PAckage (PAPA) -- written by Julian Book +/// +/// +/// It configures a multi-task with 5 configurations for the study of the TRD performance +/// in the low and intermediate masss regions. +/// +/// UPDATES, NOTES: +/// - + +PairAnalysis* Config_Analysis(Int_t cutDefinition); + +//WEIGHTS: +// 1 - 4.5 GeV 0-5% centr +// 2 - 8 GeV 0-5% centr +// 3 - 12 GeV 0-5% centr +// 4 - 4.5 GeV mbias +// 5 - 8 GeV mbias +// 6 - 12 GeV mbias +// 7 - 3.42 GeV 0-5% centr +Int_t WEIGHTS = 3; + +Bool_t GMC = false; + +/// names of the tasks +TString names = ("ACC;REC;FULL;"); +enum +{ + kACCcfg, /// for PID setup + kRECcfg, /// for PID setup + kFullPIDcfg, /// for PID setup +}; + +////// OUTPUT +void InitHistograms(PairAnalysis* papa, Int_t cutDefinition); +void AddHitHistograms(PairAnalysis* papa, Int_t cutDefinition); +void AddTrackHistograms(PairAnalysis* papa, Int_t cutDefinition); +void AddTrackHistogramsReconstruction(PairAnalysisHistos* histos, Int_t cutDefinition); +void AddPairHistograms(PairAnalysis* papa, Int_t cutDefinition); +////// CUTS +void SetupEventCuts(AnalysisTaskMultiPairAnalysis* task); +void SetupTrackCuts(PairAnalysis* papa, Int_t cutDefinition); +void SetupTrackCutsMC(PairAnalysis* papa, Int_t cutDefinition); +void AddCutStepHistograms(PairAnalysis* papa, Int_t cutDefinition); +////// SETTINGS +void AddMCSignals(PairAnalysis* papa, Int_t cutDefinition); +////// MISC +TObjArray* arrNames = names.Tokenize(";"); +const Int_t nPapa = arrNames->GetEntries(); + +//______________________________________________________________________________________ +//AnalysisTaskMultiPairAnalysis *Config_dilepton_testing() +AnalysisTaskMultiPairAnalysis* Config_dilepton_testing(Int_t id = 1, TString TASK = "", Int_t scale = 0) +{ + /// + /// creation of one multi task + /// + // AnalysisTaskMultiPairAnalysis *task = new AnalysisTaskMultiPairAnalysis("testing",id); + // AnalysisTaskMultiPairAnalysis *task = new AnalysisTaskMultiPairAnalysis("testing",id,scale); + AnalysisTaskMultiPairAnalysis* task = new AnalysisTaskMultiPairAnalysis("testing", id, WEIGHTS); + // task->SetBeamEnergy(8.); //TODO: get internally from FairBaseParSet::GetBeamMom() + + /// apply event cuts + SetupEventCuts(task); + + /// add PAPA analysis with different cuts to the task + for (Int_t i = 0; i < nPapa; ++i) { + //activate configs + switch (i) { + case kACCcfg: + continue; + // case kRECcfg: continue; + // case kFullPIDcfg: continue; + default: + // Info(" Config_jbook_AA",Form("Configure PAPa-subtask %s",((TObjString*)arrNames->At(i))->GetName())); + break; + } + + /// load configuration + PairAnalysis* papa = Config_Analysis(i); + if (!papa) continue; + + /// add PAPA to the task + task->AddPairAnalysis(papa); + std::cout << "-I- : Added subtask " << papa->GetName() << std::endl; + } + std::cout << "-I- : Added task " << task->GetName() << std::endl; + return task; +} + +//______________________________________________________________________________________ +PairAnalysis* Config_Analysis(Int_t cutDefinition) +{ + /// + /// Setup the instance of PairAnalysis + /// + + /// task name + TString name = arrNames->At(cutDefinition)->GetName(); + printf(" Adding config %s \n", name.Data()); + + /// init managing object PairAnalysis with unique name,title + PairAnalysis* papa = new PairAnalysis(Form("%s", name.Data()), Form("%s", name.Data())); + papa->SetHasMC(kTRUE); /// TODO: set automatically + /// ~ type of analysis (leptonic, hadronic or semi-leptonic 2-particle decays are supported) + papa->SetLegPdg(-11, +11); /// default: dielectron + papa->SetRefitWithMassAssump(kTRUE); /// refit under legpdg-mass assumptions + //papa->SetUseKF(kTRUE); /// use KF particle for secondary vertexing + + /* vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv CUTS vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv */ + SetupTrackCuts(papa, cutDefinition); + if (!GMC) SetupTrackCutsMC(papa, cutDefinition); + + papa->SetPreFilterUnlikeOnly(); + + /// Monte Carlo Signals + AddMCSignals(papa, cutDefinition); + + /* vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv OUTPUT vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv */ + InitHistograms(papa, cutDefinition); + + /// some simple cut QA for events, tracks, OS-pairs + papa->SetCutQA(); + + return papa; +} + +//______________________________________________________________________________________ +void SetupEventCuts(AnalysisTaskMultiPairAnalysis* task) +{ + /// + /// Setup the event cuts + /// + + /// Cut using (a formula based on) variables which are defined and described in PairAnalysisVarManager.h + /// Cuts can be added by either excluding or including criteria (see PairAnalysisVarCuts.h) + /// formula function strings are listed here: http://root.cern.ch/root/html/TFormula.html#TFormula:Analyze + PairAnalysisVarCuts* eventCuts = new PairAnalysisVarCuts("vertex", "vertex"); + // eventCuts->AddCut(PairAnalysisVarManager::kNVtxContrib, 0.0, 800.); /// inclusion cut + // eventCuts->AddCut(PairAnalysisVarManager::kImpactParam, 0.0, 13.5); + // eventCuts->AddCut("VtxChi/VtxNDF", 6., 1.e3, kTRUE); /// 'kTRUE': exclusion cut based on formula + // eventCuts->AddCut("abs(ZvPrim)", 0., 10.); /// example of TMath in formula-cuts + + /// add cuts to the global event filter + task->SetEventFilter(eventCuts); + // papa->GetEventFilter().AddCuts(eventCuts); /// or to each config + + /// for debug purpose (recommended) + eventCuts->Print(); +} + +//______________________________________________________________________________________ +void SetupTrackCuts(PairAnalysis* papa, Int_t cutDefinition) +{ + /// + /// Setup the track cuts + /// + Bool_t hasMC = papa->GetHasMC(); + + + /* vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv TRACK QUALITY CUTS vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv */ + /// setup a cut group, in which cut logic can be set (kCompAND, kCompOR) + // PairAnalysisCutGroup *grpQualCuts = new PairAnalysisCutGroup("quality","quality",PairAnalysisCutGroup::kCompAND); + + /// mcPID cuts to reject heavy particle cocktail contributions + PairAnalysisVarCuts* mcPDGcuts = new PairAnalysisVarCuts("mcPDG", "mcPDG"); + mcPDGcuts->SetCutType(PairAnalysisVarCuts::ECutType::kAll); /// wheter 'kAll' or 'kAny' cut has to be fullfilled + mcPDGcuts->AddCut(PairAnalysisVarManager::kPdgCode, 1000010020. - 0.5, 1000010020. + 0.5, kTRUE); + mcPDGcuts->AddCut(PairAnalysisVarManager::kPdgCode, 1000010030. - 0.5, 1000010030. + 0.5, kTRUE); + mcPDGcuts->AddCut(PairAnalysisVarManager::kPdgCode, 1000020030. - 0.5, 1000020030. + 0.5, kTRUE); + mcPDGcuts->AddCut(PairAnalysisVarManager::kPdgCode, 1000020040. - 0.5, 1000020040. + 0.5, kTRUE); + mcPDGcuts->AddCut(PairAnalysisVarManager::kPdgCodeMother, 1000010020. - 0.5, 1000010020. + 0.5, kTRUE); + mcPDGcuts->AddCut(PairAnalysisVarManager::kPdgCodeMother, 1000010030. - 0.5, 1000010030. + 0.5, kTRUE); + mcPDGcuts->AddCut(PairAnalysisVarManager::kPdgCodeMother, 1000020030. - 0.5, 1000020030. + 0.5, kTRUE); + mcPDGcuts->AddCut(PairAnalysisVarManager::kPdgCodeMother, 1000020040. - 0.5, 1000020040. + 0.5, kTRUE); + + /// prefilter track cuts + PairAnalysisVarCuts* preCuts = new PairAnalysisVarCuts("preCuts", "preCuts"); + preCuts->SetCutType(PairAnalysisVarCuts::ECutType::kAll); /// wheter 'kAll' or 'kAny' cut has to be fullfilled + // preCuts->AddCut(PairAnalysisVarManager::kMvdFirstHitPosZ, 0., 7.5); /// a hit in 1st MVD layer + // preCuts->AddCut("MvdHits+StsHits", 3., 15.); + // preCuts->AddCut(PairAnalysisVarManager::kStsHits, 3., 15.); + // preCuts->AddCut(PairAnalysisVarManager::kStsHits, 3., 15.); + preCuts->AddCut(PairAnalysisVarManager::kChi2NDFtoVtx, 0., 3.); /// tracks pointing to the primary vertex + + PairAnalysisVarCuts* mvdCut = new PairAnalysisVarCuts("mvdCut", "mvdCut"); + mvdCut->SetCutType(PairAnalysisVarCuts::ECutType::kAll); /// wheter 'kAll' or 'kAny' cut has to be fullfilled + mvdCut->AddCut(PairAnalysisVarManager::kMvdHits, 3., 5.); + + /// acceptance cuts (applied after pair pre filter) + PairAnalysisVarCuts* accCuts = new PairAnalysisVarCuts("accRec", "accRec"); + accCuts->SetCutType(PairAnalysisVarCuts::ECutType::kAll); /// wheter 'kAll' or 'kAny' cut has to be fullfilled + //accCuts->AddCut(PairAnalysisVarManager::kMvdHitsMC, 0., 0.5, kTRUE); // MVD MC acceptance + accCuts->AddCut(PairAnalysisVarManager::kStsHitsMC, 1., 99.); // STS MC acceptance + accCuts->AddCut(PairAnalysisVarManager::kTrdHitsMC, 1., 99.); // TRD MC acceptance + // accCuts->AddCut(PairAnalysisVarManager::kRichhasProj, 0., 0.5, kTRUE); // RICH MC acceptance + // accCuts->AddCut(PairAnalysisVarManager::kPt, 0.05, 1e30); // NOTE: was 0.2 GeV/c + + + PairAnalysisVarCuts* pT = new PairAnalysisVarCuts("pT", "pT"); + pT->AddCut(PairAnalysisVarManager::kPt, 0.05, 1e30); // NOTE: was 0.2 GeV/c + + /// standard reconstruction cuts + PairAnalysisVarCuts* recSTS = new PairAnalysisVarCuts("recSTS", "recSTS"); + recSTS->SetCutType(PairAnalysisVarCuts::ECutType::kAll); + recSTS->AddCut(PairAnalysisVarManager::kStsHits, 3., 15.); + + PairAnalysisVarCuts* recMVD = new PairAnalysisVarCuts("recMVD", "recMVD"); + recMVD->SetCutType(PairAnalysisVarCuts::ECutType::kAll); + recMVD->AddCut(PairAnalysisVarManager::kMvdHits, 3., 5.); + + /// RICH acceptance cuts + PairAnalysisVarCuts* accRICH = new PairAnalysisVarCuts("accRICH", "accRICH"); + accRICH->SetCutType(PairAnalysisVarCuts::ECutType::kAll); + accRICH->AddCut(PairAnalysisVarManager::kRichhasProj, 0., 0.5, kTRUE); + + /// RICH reconstruction cuts + PairAnalysisVarCuts* recRICH = new PairAnalysisVarCuts("recRICH", "recRICH"); + recRICH->SetCutType(PairAnalysisVarCuts::ECutType::kAll); + recRICH->AddCut(PairAnalysisVarManager::kRichHits, 6., 100.); /// min+max inclusion for hits, 13=eff95% + // valid ring parameters + recRICH->AddCut(PairAnalysisVarManager::kRichAxisA, 0., 10.); + recRICH->AddCut(PairAnalysisVarManager::kRichAxisB, 0., 10.); + recRICH->AddCut(PairAnalysisVarManager::kRichDistance, 0., 999.); + recRICH->AddCut(PairAnalysisVarManager::kRichRadialPos, 0., 999.); + recRICH->AddCut(PairAnalysisVarManager::kRichRadialAngle, -TMath::TwoPi(), TMath::TwoPi()); + recRICH->AddCut(PairAnalysisVarManager::kRichPhi, -TMath::TwoPi(), TMath::TwoPi()); + // ellipse fit not working (chi2/ndf is nan) + recRICH->AddCut(PairAnalysisVarManager::kRichChi2NDF, 0., 100.); + + + /// TRD reconstruction cuts + PairAnalysisVarCuts* recTRD = new PairAnalysisVarCuts("recTRD", "recTRD"); + recTRD->SetCutType(PairAnalysisVarCuts::ECutType::kAll); + recTRD->AddCut(PairAnalysisVarManager::kTrdHits, 3., 10.); /// min+max requieremnt for hits + // recTRD->AddCut(PairAnalysisVarManager::kElossNew, 27., 35.); /// min+max requieremnt for hits + + /// TOF reconstruction cuts + PairAnalysisVarCuts* recTOF = new PairAnalysisVarCuts("recTOF", "recTOF"); + recTOF->SetCutType(PairAnalysisVarCuts::ECutType::kAll); + recTOF->AddCut(PairAnalysisVarManager::kTofHits, 1., 10.0); /// min+max requieremnt for hits + // recTOF->AddCut(PairAnalysisVarManager::kP, 0., 0.8); /// momentum selection + + PairAnalysisVarCuts* accPID = new PairAnalysisVarCuts("accPID", "accPID"); + accPID->SetCutType(PairAnalysisVarCuts::ECutType::kAll); + accPID->AddCut(PairAnalysisVarManager::kRichhasProj, 0., 0.5, kTRUE); + accPID->AddCut(PairAnalysisVarManager::kTrdHits, 2., 10.); + + /// TRD pid cuts - 2-dimensional + + TGraph* grTRD = NULL; // lower cut limt + TGraph* grMax = new TGraph(2); // upper cut limit + grMax->SetPoint(0, 0., 999.); + grMax->SetPoint(1, 999., 999.); + grMax->GetListOfFunctions()->Add(PairAnalysisHelper::GetFormula("varf", "P")); + + + Double_t x8[200] = { + 0.101342, 0.104063, 0.106857, 0.109725, 0.112671, 0.115695, 0.118801, 0.121991, 0.125266, 0.128628, 0.132082, + 0.135627, 0.139268, 0.143007, 0.146846, 0.150788, 0.154836, 0.158993, 0.163261, 0.167644, 0.172145, 0.176766, + 0.181511, 0.186384, 0.191388, 0.196526, 0.201802, 0.207219, 0.212782, 0.218494, 0.22436, 0.230383, 0.236568, + 0.242919, 0.24944, 0.256136, 0.263012, 0.270073, 0.277323, 0.284768, 0.292413, 0.300263, 0.308324, 0.316601, + 0.3251, 0.333828, 0.34279, 0.351992, 0.361441, 0.371144, 0.381108, 0.391339, 0.401845, 0.412633, 0.42371, + 0.435085, 0.446765, 0.458759, 0.471074, 0.483721, 0.496706, 0.510041, 0.523733, 0.537793, 0.55223, 0.567055, + 0.582278, 0.59791, 0.613961, 0.630443, 0.647368, 0.664747, 0.682592, 0.700917, 0.719734, 0.739055, 0.758896, + 0.779269, 0.800189, 0.82167, 0.843728, 0.866379, 0.889637, 0.91352, 0.938044, 0.963226, 0.989085, 1.01564, + 1.0429, 1.0709, 1.09965, 1.12917, 1.15948, 1.19061, 1.22257, 1.25539, 1.2891, 1.3237, 1.35924, + 1.39573, 1.4332, 1.47167, 1.51118, 1.55175, 1.59341, 1.63618, 1.68011, 1.72521, 1.77152, 1.81908, + 1.86792, 1.91806, 1.96955, 2.02243, 2.07672, 2.13247, 2.18972, 2.2485, 2.30886, 2.37085, 2.43449, + 2.49985, 2.56696, 2.63587, 2.70663, 2.77929, 2.85391, 2.93052, 3.00919, 3.08998, 3.17293, 3.25811, + 3.34557, 3.43539, 3.52761, 3.62231, 3.71956, 3.81941, 3.92194, 4.02723, 4.13534, 4.24636, 4.36036, + 4.47741, 4.59761, 4.72104, 4.84778, 4.97792, 5.11155, 5.24878, 5.38968, 5.53437, 5.68295, 5.83551, + 5.99217, 6.15303, 6.31821, 6.48783, 6.662, 6.84084, 7.02449, 7.21306, 7.4067, 7.60554, 7.80972, + 8.01937, 8.23466, 8.45572, 8.68272, 8.91581, 9.15516, 9.40094, 9.65331, 9.91246, 10.1786, 10.4518, + 10.7324, 11.0205, 11.3164, 11.6202, 11.9321, 12.2524, 12.5814, 12.9191, 13.2659, 13.6221, 13.9878, + 14.3633, 14.7489, 15.1448, 15.5514, 15.9689, 16.3976, 16.8378, 17.2898, 17.7539, 18.2306, 18.72, + 19.2225, 19.7386}; + Double_t y8[200] = { + 0.075, 0.075, 0.075, 0.075, 0.075, 0.075, 0.075, 0.075, 0.075, 0.075, 0.075, 0.075, 0.075, 0.075, 0.075, 0.075, + 0.075, 0.005, 0.075, 0.075, 0.075, 0.075, 0.075, 0.075, 0.075, 0.075, 0.075, 0.005, 0.005, 0.075, 0.305, 0.035, + 0.115, 0.125, 0.135, 0.115, 0.115, 0.075, 0.155, 0.115, 0.145, 0.105, 0.155, 0.155, 0.175, 0.165, 0.205, 0.185, + 0.145, 0.175, 0.165, 0.155, 0.165, 0.155, 0.175, 0.165, 0.185, 0.165, 0.165, 0.155, 0.175, 0.165, 0.175, 0.165, + 0.165, 0.165, 0.175, 0.185, 0.175, 0.175, 0.185, 0.175, 0.185, 0.195, 0.195, 0.195, 0.195, 0.185, 0.205, 0.195, + 0.195, 0.195, 0.195, 0.195, 0.185, 0.185, 0.185, 0.205, 0.215, 0.215, 0.215, 0.225, 0.225, 0.225, 0.225, 0.245, + 0.245, 0.255, 0.255, 0.255, 0.255, 0.265, 0.295, 0.305, 0.315, 0.345, 0.315, 0.355, 0.355, 0.365, 0.375, 0.375, + 0.405, 0.385, 0.405, 0.425, 0.435, 0.435, 0.435, 0.455, 0.445, 0.465, 0.445, 0.455, 0.455, 0.505, 0.495, 0.505, + 0.495, 0.525, 0.515, 0.535, 0.525, 0.495, 0.515, 0.515, 0.515, 0.495, 0.525, 0.515, 0.535, 0.475, 0.505, 0.485, + 0.505, 0.495, 0.515, 0.505, 0.485, 0.475, 0.255, 0.265, 0.295, 0.305, 0.315, 0.345, 0.315, 0.355, 0.355, 0.365, + 0.375, 0.375, 0.405, 0.385, 0.405, 0.425, 0.435, 0.435, 0.435, 0.455, 0.445, 0.465, 0.445, 0.455, 0.455, 0.505, + 0.495, 0.505, 0.495, 0.525, 0.515, 0.535, 0.525, 0.495, 0.515, 0.515, 0.515, 0.495, 0.525, 0.515, 0.535, 0.475, + 0.505, 0.485, 0.505, 0.495, 0.515, 0.505, 0.485, 0.475}; + + + Bool_t LH = kFALSE; + grTRD = new TGraph(200, x8, y8); + LH = kTRUE; + if (grTRD) grTRD->GetListOfFunctions()->Add(PairAnalysisHelper::GetFormula("varf", "P")); + + /// input to cut object + PairAnalysisObjectCuts* pidTRD2d = new PairAnalysisObjectCuts("pidTRD2d", "pidTRD2d"); + pidTRD2d->SetCutType(PairAnalysisObjectCuts::ECutType::kAll); /// wheter kAll or kAny cut has to be fullfilled + if (LH) pidTRD2d->AddCut(PairAnalysisVarManager::kTrdPidLikeEL, grTRD, grMax); + else + pidTRD2d->AddCut(PairAnalysisVarManager::kTrdPidANN, grTRD, grMax); + + // /// RICH pid cuts - 2-dimensional + Double_t xR[200] = { + 0.302653, 0.308006, 0.313454, 0.318998, 0.324641, 0.330383, 0.336226, 0.342173, 0.348225, 0.354384, 0.360653, + 0.367032, 0.373523, 0.38013, 0.386854, 0.393696, 0.400659, 0.407746, 0.414958, 0.422298, 0.429767, 0.437368, + 0.445104, 0.452977, 0.460989, 0.469143, 0.477441, 0.485885, 0.494479, 0.503225, 0.512126, 0.521184, 0.530403, + 0.539784, 0.549332, 0.559048, 0.568936, 0.578999, 0.58924, 0.599662, 0.610268, 0.621062, 0.632047, 0.643227, + 0.654604, 0.666182, 0.677965, 0.689956, 0.70216, 0.714579, 0.727218, 0.740081, 0.753171, 0.766493, 0.78005, + 0.793847, 0.807888, 0.822177, 0.83672, 0.851519, 0.86658, 0.881908, 0.897506, 0.913381, 0.929536, 0.945977, + 0.962709, 0.979737, 0.997066, 1.0147, 1.03265, 1.05091, 1.0695, 1.08842, 1.10767, 1.12726, 1.1472, + 1.16749, 1.18814, 1.20916, 1.23054, 1.25231, 1.27446, 1.297, 1.31994, 1.34329, 1.36705, 1.39123, + 1.41583, 1.44087, 1.46636, 1.4923, 1.51869, 1.54555, 1.57289, 1.60071, 1.62902, 1.65784, 1.68716, + 1.717, 1.74737, 1.77827, 1.80973, 1.84174, 1.87431, 1.90746, 1.9412, 1.97554, 2.01048, 2.04604, + 2.08223, 2.11906, 2.15654, 2.19468, 2.2335, 2.27301, 2.31321, 2.35412, 2.39576, 2.43814, 2.48126, + 2.52515, 2.56981, 2.61526, 2.66152, 2.7086, 2.7565, 2.80526, 2.85488, 2.90537, 2.95676, 3.00906, + 3.06228, 3.11645, 3.17157, 3.22766, 3.28475, 3.34285, 3.40198, 3.46215, 3.52339, 3.58571, 3.64913, + 3.71367, 3.77936, 3.8462, 3.91423, 3.98347, 4.05392, 4.12563, 4.1986, 4.27286, 4.34844, 4.42535, + 4.50362, 4.58328, 4.66434, 4.74684, 4.8308, 4.91625, 5.0032, 5.0917, 5.18176, 5.27341, 5.36668, + 5.4616, 5.5582, 5.65651, 5.75656, 5.85838, 5.962, 6.06745, 6.17477, 6.28399, 6.39513, 6.50825, + 6.62336, 6.74051, 6.85973, 6.98106, 7.10454, 7.2302, 7.35808, 7.48823, 7.62068, 7.75547, 7.89264, + 8.03224, 8.17431, 8.31889, 8.46603, 8.61577, 8.76817, 8.92325, 9.08108, 9.2417, 9.40516, 9.57152, + 9.74081, 9.9131}; + Double_t yR[200] = { + -1.17, -1.05, -1.11, -1.03, -1.11, -1.07, -1.05, -1.05, -1.03, -1.03, -1.01, -1.03, -0.97, -0.99, -1.01, -1.05, + -0.97, -0.99, -1.03, -0.99, -0.99, -0.91, -0.95, -0.87, -0.97, -0.93, -0.95, -0.89, -0.95, -0.95, -0.89, -0.85, + -0.91, -0.95, -0.91, -0.81, -0.87, -0.77, -0.77, -0.85, -0.79, -0.63, -0.69, -0.59, -0.65, -0.55, -0.77, -0.65, + -0.61, -0.63, -0.55, -0.49, -0.47, -0.63, -0.53, -0.51, -0.37, -0.27, -0.35, -0.41, -0.35, -0.33, -0.23, -0.25, + -0.23, -0.13, -0.11, -0.25, -0.19, -0.17, -0.13, -0.11, -0.17, -0.01, -0.01, 0.05, 0.05, 0.09, 0.09, 0.07, + 0.15, 0.13, 0.23, 0.19, 0.19, 0.25, 0.31, 0.27, 0.35, 0.33, 0.35, 0.33, 0.39, 0.39, 0.45, 0.45, + 0.45, 0.43, 0.47, 0.53, 0.53, 0.55, 0.57, 0.51, 0.57, 0.59, 0.61, 0.59, 0.63, 0.63, 0.65, 0.63, + 0.63, 0.67, 0.67, 0.67, 0.71, 0.69, 0.71, 0.71, 0.73, 0.71, 0.75, 0.75, 0.75, 0.75, 0.77, 0.75, + 0.79, 0.79, 0.79, 0.81, 0.79, 0.83, 0.83, 0.79, 0.81, 0.81, 0.83, 0.85, 0.83, 0.85, 0.85, 0.83, + 0.83, 0.85, 0.85, 0.85, 0.85, 0.85, 0.85, 0.87, 0.87, 0.85, 0.85, 0.85, 0.87, 0.87, 0.85, 0.87, + 0.87, 0.87, 0.85, 0.85, 0.87, 0.87, 0.87, 0.87, 0.87, 0.87, 0.87, 0.87, 0.89, 0.87, 0.85, 0.87, + 0.87, 0.87, 0.87, 0.87, 0.87, 0.85, 0.87, 0.87, 0.87, 0.87, 0.87, 0.89, 0.87, 0.87, 0.87, 0.87, + 0.89, 0.89, 0.89, 0.89, 0.87, 0.87, 0.89, 0.89}; + + + TGraph* grR = new TGraph(200, xR, yR); + grR->GetListOfFunctions()->Add(PairAnalysisHelper::GetFormula("varf", "P")); + + PairAnalysisObjectCuts* pidRICH2d = new PairAnalysisObjectCuts("pidRICH2d", "pidRICH2d"); + pidRICH2d->SetCutType(PairAnalysisObjectCuts::ECutType::kAll); /// wheter kAll or kAny cut has to be fullfilled + pidRICH2d->AddCut(PairAnalysisVarManager::kRichPidANN, grR, grMax); + + + Double_t mvP[8] = {0., 0.2, 0.4, 0.6, 0.8, 1., 1.2, 99.}; + Double_t mvD[8] = {0.4, 0.333, 0.2666, 0.2, 0.1333, 0.06666, 0., 0.}; + TGraph* grMVD = new TGraph(7, mvP, mvD); + grMVD->GetListOfFunctions()->Add(PairAnalysisHelper::GetFormula("varf", "P")); + TGraph* grMVDMax = new TGraph(2); + grMVDMax->SetPoint(0, 0., 999.); + grMVDMax->SetPoint(1, 100., 999.); + grMVDMax->GetListOfFunctions()->Add(PairAnalysisHelper::GetFormula("varf", "P")); + PairAnalysisObjectCuts* MVDClosest = new PairAnalysisObjectCuts("MVDClosest", "MVDClosest"); + MVDClosest->SetCutType(PairAnalysisObjectCuts::ECutType::kAll); /// wheter kAll or kAny cut has to be fullfilled + MVDClosest->AddCut(PairAnalysisVarManager::kMvdHitClosest, grMVD, grMVDMax); + + Double_t mvS[8] = {0., 0.2, 0.4, 0.6, 0.8, 1., 1.2, 99.}; + Double_t mvSD[8] = {1., 0.8333, 0.6666, 0.5, 0.3333, 0.1666, 0., 0.}; + TGraph* grSTS = new TGraph(7, mvS, mvSD); + grSTS->GetListOfFunctions()->Add(PairAnalysisHelper::GetFormula("varf", "P")); + TGraph* grSTSMax = new TGraph(2); + grSTSMax->SetPoint(0, 0., 999.); + grSTSMax->SetPoint(1, 100., 999.); + grSTSMax->GetListOfFunctions()->Add(PairAnalysisHelper::GetFormula("varf", "P")); + PairAnalysisObjectCuts* STSClosest = new PairAnalysisObjectCuts("STSClosest", "STSClosest"); + STSClosest->SetCutType(PairAnalysisObjectCuts::ECutType::kAll); /// wheter kAll or kAny cut has to be fullfilled + STSClosest->AddCut(PairAnalysisVarManager::kStsHitClosest, grSTS, grSTSMax); + + + Double_t mvdPO[12] = {0., 0.2, 0.4, 0.6, 0.8, 1., 1.2, 1.4, 1.6, 1.8, 2., 99.}; + Double_t mvdDO[12] = {0.035, 0.0315, 0.028, 0.0245, 0.021, 0.0175, 0.014, 0.0105, 0.007, 0.0035, 0., 0.}; + TGraph* grMVDOp = new TGraph(12, mvdPO, mvdDO); + grMVDOp->GetListOfFunctions()->Add(PairAnalysisHelper::GetFormula("varf", "MVDHitClosestMom")); + TGraph* grMVDMaxOp = new TGraph(2); + grMVDMaxOp->SetPoint(0, 0., 999.); + grMVDMaxOp->SetPoint(1, 100., 999.); + grMVDMaxOp->GetListOfFunctions()->Add(PairAnalysisHelper::GetFormula("varf", "MVDHitClosestMom")); + PairAnalysisObjectCuts* MVDClosestOpeningAngle = + new PairAnalysisObjectCuts("MVDClosestOpeningAngle", "MVDClosestOpeningAngle"); + MVDClosestOpeningAngle->SetCutType( + PairAnalysisObjectCuts::ECutType::kAll); /// wheter kAll or kAny cut has to be fullfilled + MVDClosestOpeningAngle->AddCut(PairAnalysisVarManager::kMvdHitClosestOpeningAngle, grMVDOp, grMVDMaxOp); + + Double_t stsPO[12] = {0., 0.2, 0.4, 0.6, 0.8, 1., 1.2, 1.4, 1.6, 1.8, 2., 99.}; + Double_t stsDO[12] = {0.035, 0.0315, 0.028, 0.0245, 0.021, 0.0175, 0.014, 0.0105, 0.007, 0.0035, 0., 0.}; + TGraph* grSTSOp = new TGraph(12, stsPO, stsDO); + grSTSOp->GetListOfFunctions()->Add(PairAnalysisHelper::GetFormula("varf", "STSHitClosestMom")); + TGraph* grSTSMaxOp = new TGraph(2); + grSTSMaxOp->SetPoint(0, 0., 999.); + grSTSMaxOp->SetPoint(1, 100., 999.); + grSTSMaxOp->GetListOfFunctions()->Add(PairAnalysisHelper::GetFormula("varf", "STSHitClosestMom")); + PairAnalysisObjectCuts* STSClosestOpeningAngle = + new PairAnalysisObjectCuts("STSClosestOpeningAngle", "STSClosestOpeningAngle"); + STSClosestOpeningAngle->SetCutType( + PairAnalysisObjectCuts::ECutType::kAll); /// wheter kAll or kAny cut has to be fullfilled + STSClosestOpeningAngle->AddCut(PairAnalysisVarManager::kStsHitClosestOpeningAngle, grSTSOp, grSTSMaxOp); + + // TOF Pid + PairAnalysisVarCuts* pidTOF = new PairAnalysisVarCuts("pidTOF", "pidTOF"); + pidTOF->AddCut(PairAnalysisVarManager::kTofPidDeltaBetaEL, -1.65 * 3.2e-03, +1.65 * 3.2e-03); //90% + // pidTOF->AddCut(PairAnalysisVarManager::kMassSq, -0.2, 0.2); + + PairAnalysisCutCombi* pidTOFavai = new PairAnalysisCutCombi("TOFPidAvai", "TOFPidAvai"); + pidTOFavai->AddCut(pidTOF, recTOF); + + /* vvvvvvvvvvvvvvvvvvvvvvvvvvvvvv PAIR PREFILTER CUTS vvvvvvvvvvvvvvvvvvvvvvvvvvvvv */ + /// example of gamma rejection cuts as prefilter in order to remove tracks from + /// the final track array. NOTE: inverted cut logic! + + /// rejection based on pair informations + /// leg cuts for pair prefilter + /// rejection based on pair informations + PairAnalysisVarCuts* gammaCuts = new PairAnalysisVarCuts("gammaCut", "gammaCut"); + gammaCuts->AddCut(PairAnalysisVarManager::kM, 0., 0.025); // exponential cut + + PairAnalysisVarCuts* prepairMvdCuts = new PairAnalysisVarCuts("prepairMVD", "prepairMVD"); + prepairMvdCuts->AddCut(PairAnalysisVarManager::kMvdHitDist, 0., 1.); // exponential cut + + PairAnalysisVarCuts* prepairStsCuts = new PairAnalysisVarCuts("prepairSTS", "prepairSTS"); + prepairStsCuts->AddCut(PairAnalysisVarManager::kStsHitDist, 0., 5.); // exponential cut + + + PairAnalysisVarCuts* preRecMvdCut = new PairAnalysisVarCuts("preRecMvdCut", "preRecMvdCut"); + preRecMvdCut->SetCutType(PairAnalysisVarCuts::ECutType::kAll); /// wheter 'kAll' or 'kAny' cut has to be fullfilled + preRecMvdCut->AddCut(PairAnalysisVarManager::kMvdHits, 3., 4.); + preRecMvdCut->AddCut(PairAnalysisVarManager::kStsHits, 3., 20.); + preRecMvdCut->AddCut(PairAnalysisVarManager::kRichHits, 6., 99.); + preRecMvdCut->AddCut(PairAnalysisVarManager::kTrdHits, 3., 5.); + + PairAnalysisVarCuts* preRecStsCut = new PairAnalysisVarCuts("preRecStsCut", "preRecStsCut"); + preRecStsCut->SetCutType(PairAnalysisVarCuts::ECutType::kAll); /// wheter 'kAll' or 'kAny' cut has to be fullfilled + preRecStsCut->AddCut(PairAnalysisVarManager::kStsHits, 3., 20.); + preRecStsCut->AddCut(PairAnalysisVarManager::kRichHits, 6., 99.); + preRecStsCut->AddCut(PairAnalysisVarManager::kTrdHits, 3., 5.); + + + /* vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv ACTIVATE CUTS vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv */ + /// activate the cut sets (order might be CPU timewise important) + switch (cutDefinition) { + case kACCcfg: + papa->GetTrackFilter().AddCuts(mcPDGcuts); + mcPDGcuts->Print(); + papa->GetTrackFilter().AddCuts(accCuts); + accCuts->Print(); + break; + + case kRECcfg: + papa->GetTrackFilter().AddCuts(mcPDGcuts); + mcPDGcuts->Print(); + papa->GetTrackFilter().AddCuts(preCuts); + preCuts->Print(); + papa->GetFinalTrackFilter().AddCuts(recMVD); + recMVD->Print(); + papa->GetFinalTrackFilter().AddCuts(recSTS); + recSTS->Print(); + papa->GetFinalTrackFilter().AddCuts(recRICH); + recRICH->Print(); + papa->GetFinalTrackFilter().AddCuts(recTRD); + recTRD->Print(); + papa->GetFinalTrackFilter().AddCuts(recTOF); + recTOF->Print(); + break; + + case kFullPIDcfg: + papa->GetTrackFilter().AddCuts(mcPDGcuts); + mcPDGcuts->Print(); + papa->GetTrackFilter().AddCuts(preCuts); + preCuts->Print(); + papa->GetTrackFilter().AddCuts(accCuts); + accCuts->Print(); + papa->GetTrackFilter().AddCuts(pT); + pT->Print(); + papa->GetTrackFilter().AddCuts(preRecMvdCut); + preRecMvdCut->Print(); + papa->GetPairPreFilter().AddCuts(gammaCuts); + gammaCuts->Print(); + papa->GetFinalTrackFilter().AddCuts(MVDClosestOpeningAngle); + MVDClosestOpeningAngle->Print(); + papa->GetFinalTrackFilter().AddCuts(recRICH); + recRICH->Print(); + papa->GetFinalTrackFilter().AddCuts(pidRICH2d); + pidRICH2d->Print(); + papa->GetFinalTrackFilter().AddCuts(recTRD); + recTRD->Print(); + papa->GetFinalTrackFilter().AddCuts(pidTRD2d); + pidTRD2d->Print(); + papa->GetFinalTrackFilter().AddCuts(pidTOFavai); + pidTOFavai->Print(); + break; + } +} + +//______________________________________________________________________________________ +void SetupTrackCutsMC(PairAnalysis* papa, Int_t cutDefinition) +{ + /// + /// Setup the track cuts based on MC information only + /// + if (!papa->GetHasMC()) return; + + /// NOTE: skip this + // return; + + /* vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv TRACK CUTS ON MCtruth vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv */ + /// example of adding acceptance cuts + PairAnalysisVarCuts* accCutsMC = new PairAnalysisVarCuts("accGen", "accGen"); + accCutsMC->SetCutType(PairAnalysisVarCuts::ECutType::kAll); /// wheter kAll or kAny cut has to be fullfilled + + /// example for config specific cuts + switch (cutDefinition) { + default: + accCutsMC->AddCut(PairAnalysisVarManager::kGeantId, 36.5, 37.5, kTRUE); + // accCutsMC->AddCut(PairAnalysisVarManager::kMvdHitsMC, 0., 0.5, kTRUE); // MVD MC acceptance + accCutsMC->AddCut(PairAnalysisVarManager::kStsHitsMC, 0., 0.5, kTRUE); // STS MC acceptance + // accCutsMC->AddCut(PairAnalysisVarManager::kRichHitsMC, 0., 0.5, kTRUE); // STS MC acceptance + accCutsMC->AddCut(PairAnalysisVarManager::kTrdHitsMC, 0., 0.5, kTRUE); // TRD MC acceptance + // accCutsMC->AddCut(PairAnalysisVarManager::kRichhasProj, 0., 0.5, kTRUE); // RICH MC acceptance + // accCutsMC->AddCut(PairAnalysisVarManager::kPtMC, 0.2, 1e30); + break; + } + /* vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv ACTIVATE CUTS vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv */ + // switch(cutDefinition) { + // default: + papa->GetTrackFilterMC().AddCuts(accCutsMC); + accCutsMC->Print(); + // } +} + +//______________________________________________________________________________________ +void InitHistograms(PairAnalysis* papa, Int_t cutDefinition) +{ + /// + /// Initialise the histograms + /// + /// NOTE: Histograms are added to a specific class type such as 'event,pair,track'. + /// Several different objects of x-dimensions can be added via 'AddHisogram', 'AddProfile', 'AddSparse'. + /// Please have a look at PairAnalysisHistos.h to understand all possible input arguments. + /// For histograms and profiles you can use formulas with full 'TMath'-support to calculate new variables: + /// checkout http://root.cern.ch/root/html/TFormula.html#TFormula:Analyze for formula function strings + /// The type of binning is provided by some 'PairAnalysisHelper' functions. + /// Some examples are given below.... + + Bool_t hasMC = papa->GetHasMC(); + + ///Setup histogram Manager + PairAnalysisHistos* histos = new PairAnalysisHistos(papa->GetName(), papa->GetTitle()); + papa->SetHistogramManager(histos); + histos->SetPrecision(PairAnalysisHistos::Eprecision::kFloat); + + /// Initialize superior histogram classes (reserved words) + histos->SetReservedWords("Hit;Track;Pair"); + + ////// PAIR HISTOS ///// + AddPairHistograms(papa, cutDefinition); + + ///// TRACK HISTOS ///// + AddTrackHistograms(papa, cutDefinition); + + // AddHitHistograms(papa, cutDefinition); + + if (!GMC) AddCutStepHistograms(papa, cutDefinition); + + + ////// DEBUG ////// + TIter nextClass(histos->GetHistogramList()); + THashList* l = 0; + while ((l = static_cast<THashList*>(nextClass()))) { + printf(" [D] HistogramManger: Class %s: Histograms: %04d \n", l->GetName(), l->GetEntries()); + } + +} //end: init histograms + + +//______________________________________________________________________________________ +void AddMCSignals(PairAnalysis* papa, Int_t cutDefinition) +{ + /// Do we have an MC handler? + if (!papa->GetHasMC()) return; + + + Double_t Momega = 0; + Double_t Mphi = 0; + Double_t Meta = 0; + Double_t Metap = 0; + Double_t Mrho = 0; + Double_t Mpi0 = 0; + Double_t Minmed = 0; + Double_t Mqgp = 0; + Double_t MDp = 0; + + if (WEIGHTS == 1) { + Momega = 0.270774; + Mphi = 0.022367; + Meta = 1.43997; + Metap = 0.0214123; + Mrho = 0.919404; + Minmed = 0.0100221; + Mqgp = 0.000483824; + MDp = 70.5; + } + if (WEIGHTS == 2) { + Momega = 2.28721; + Mphi = 0.311619; + Meta = 7.45497; + Metap = 0.262172; + Mrho = 7.57283; + Minmed = 0.0304706; + Mqgp = 0.000452941; + MDp = 168; + } + if (WEIGHTS == 3) { + Momega = 5.821; + Mphi = 1.018; + Meta = 14.8779; + Metap = 0.803; + Mrho = 19.0323; + Minmed = 0.041625; + Mqgp = 0.00061875; + MDp = 229.5; + } + if (WEIGHTS == 7) { + Momega = 0.270774; + Mphi = 0.022367; + Meta = 1.43997; + Metap = 0.0214123; + Mrho = 0.919404; + Minmed = 0.0100221; + Mqgp = 0.00; + MDp = 70.5; + } + + Double_t centr_mbias_ratio = 0.2613; + if (WEIGHTS == 4) { + Momega = 0.380235 * centr_mbias_ratio; + Mphi = 0.031409 * centr_mbias_ratio; + Meta = 2.02209 * centr_mbias_ratio; + Metap = 0.0300683 * centr_mbias_ratio; + Mrho = 1.29108 * centr_mbias_ratio; + Minmed = 0.0140735 * centr_mbias_ratio; + Mqgp = 0.000679412 * centr_mbias_ratio; + MDp = 99 * centr_mbias_ratio; + } + if (WEIGHTS == 5) { + Momega = 2.28721 * centr_mbias_ratio; + Mphi = 0.311619 * centr_mbias_ratio; + Meta = 7.45497 * centr_mbias_ratio; + Metap = 0.262172 * centr_mbias_ratio; + Mrho = 7.57283 * centr_mbias_ratio; + Minmed = 0.0304706 * centr_mbias_ratio; + Mqgp = 0.000452941 * centr_mbias_ratio; + MDp = 168 * centr_mbias_ratio; + } + if (WEIGHTS == 6) { + Momega = 5.821 * centr_mbias_ratio; + Mphi = 1.018 * centr_mbias_ratio; + Meta = 14.8779 * centr_mbias_ratio; + Metap = 0.803 * centr_mbias_ratio; + Mrho = 19.0323 * centr_mbias_ratio; + Minmed = 0.041625 * centr_mbias_ratio; + Mqgp = 0.00061875 * centr_mbias_ratio; + MDp = 229.5 * centr_mbias_ratio; + } + + ///// single particle signals ///// + TString fillMC = kFALSE; + PairAnalysisSignalMC* deltaele = new PairAnalysisSignalMC(PairAnalysisSignalMC::EDefinedSignal::kDeltaElectron); + deltaele->SetFillPureMCStep(fillMC); + + PairAnalysisSignalMC* eleGam = new PairAnalysisSignalMC("e^{#gamma}", "eGam"); + eleGam->SetFillPureMCStep(fillMC); + eleGam->SetLegPDGs(11, 0); + eleGam->SetCheckBothChargesLegs(kTRUE, kFALSE); + eleGam->SetMotherPDGs(22, 0); //0:default all + eleGam->SetMothersRelation(PairAnalysisSignalMC::EBranchRelation::kDifferent); + + PairAnalysisSignalMC* eleGamPi = new PairAnalysisSignalMC("e^{#gamma}Pi", "eGamPi"); + eleGamPi->SetFillPureMCStep(fillMC); + eleGamPi->SetLegPDGs(11, 0); + eleGamPi->SetCheckBothChargesLegs(kTRUE, kFALSE); + eleGamPi->SetMotherPDGs(111, 0); //0:default all + eleGamPi->SetMothersRelation(PairAnalysisSignalMC::EBranchRelation::kDifferent); + + PairAnalysisSignalMC* elePi = new PairAnalysisSignalMC("ePi", "ePi"); + elePi->SetFillPureMCStep(fillMC); + elePi->SetLegPDGs(11, 1); + elePi->SetIsSingleParticle(kTRUE); + elePi->SetCheckBothChargesLegs(kTRUE, kFALSE); + elePi->SetMotherPDGs(111, 0); //0:default all + + PairAnalysisSignalMC* conv = new PairAnalysisSignalMC(PairAnalysisSignalMC::EDefinedSignal::kConversion); + conv->SetFillPureMCStep(GMC); + + PairAnalysisSignalMC* conv2 = new PairAnalysisSignalMC("X #gamma e (comb.)", "X #gamma e (comb.)"); + conv2->SetFillPureMCStep(fillMC); + conv2->SetLegPDGs(11, -11); + conv2->SetMotherPDGs(22, 22); + conv2->SetGrandMotherPDGs(111, 111); + conv2->SetMothersRelation(PairAnalysisSignalMC::EBranchRelation::kSame); + conv2->SetCheckBothChargesLegs(kTRUE, kTRUE); + + + PairAnalysisSignalMC* ele = new PairAnalysisSignalMC(PairAnalysisSignalMC::EDefinedSignal::kPrimElectron); + ele->SetFillPureMCStep(GMC); + // ele->SetMothersRelation(PairAnalysisSignalMC::kSame); + + PairAnalysisSignalMC* e = new PairAnalysisSignalMC("e", "e"); + e->SetFillPureMCStep(fillMC); + e->SetIsSingleParticle(kTRUE); + e->SetLegPDGs(11, 1); + e->SetCheckBothChargesLegs(kTRUE, kFALSE); + // e->SetGEANTProcess(kPPrimary); + // e->SetInverse(); + + PairAnalysisSignalMC* pio = new PairAnalysisSignalMC(PairAnalysisSignalMC::EDefinedSignal::kPrimPion); + pio->SetFillPureMCStep(fillMC); + Double_t br = 1.0; // branching ratio + + // omega dalitz decays + PairAnalysisSignalMC* omegaDalitz = new PairAnalysisSignalMC(PairAnalysisSignalMC::EDefinedSignal::kOmegaDalitz); + omegaDalitz->SetFillPureMCStep(GMC); + br = 7.7e-04; + omegaDalitz->SetWeight(Momega * br); //HSD + + PairAnalysisSignalMC* etap = new PairAnalysisSignalMC("etap", "etap"); + etap->SetFillPureMCStep(GMC); + etap->SetIsDalitz(PairAnalysisSignalMC::EDalitz::kIsDalitz, 22); + etap->SetLegPDGs(-11, 11); + etap->SetCheckBothChargesLegs(kTRUE, kTRUE); + etap->SetMotherPDGs(331, 331); + etap->SetMothersRelation(PairAnalysisSignalMC::EBranchRelation::kSame); + br = 4.73e-4; + etap->SetWeight(Meta * br); + + PairAnalysisSignalMC* pibend = new PairAnalysisSignalMC("bendpi", "bendpi"); + pibend->SetFillPureMCStep(fillMC); + pibend->SetLegPDGs(-11, 11); + pibend->SetCheckBothChargesLegs(kTRUE, kTRUE); + pibend->SetMotherPDGs(111, 111); + // pibend->SetMothersRelation(PairAnalysisSignalMC::EBranchRelation::kSame); + pibend->SetIsSingleParticle(kTRUE); + pibend->SetGEANTProcess(kPMagneticFieldL); + + // omega + PairAnalysisSignalMC* omega = new PairAnalysisSignalMC(PairAnalysisSignalMC::EDefinedSignal::kOmega); + omega->SetFillPureMCStep(GMC); + br = 7.28e-04; + // omega->SetWeight(5.389 * br);//HSD + omega->SetWeight(Momega * br); //HSD + + // pi0 + PairAnalysisSignalMC* pi0 = new PairAnalysisSignalMC(PairAnalysisSignalMC::EDefinedSignal::kPi0); + pi0->SetFillPureMCStep(fillMC); + + + // pi0 -> gamma gamma + PairAnalysisSignalMC* pi0Gamma = new PairAnalysisSignalMC(PairAnalysisSignalMC::EDefinedSignal::kPi0Gamma); + pi0Gamma->SetFillPureMCStep(fillMC); + + // pi0 -> e+e- gamma + PairAnalysisSignalMC* pi0Dalitz = new PairAnalysisSignalMC(PairAnalysisSignalMC::EDefinedSignal::kPi0Dalitz); + // pi0Dalitz->SetFillPureMCStep(fillMC); + pi0Dalitz->SetFillPureMCStep(GMC); + pi0Dalitz->SetIsDalitz(PairAnalysisSignalMC::EDalitz::kWhoCares, 0); + + PairAnalysisSignalMC* eta = new PairAnalysisSignalMC(PairAnalysisSignalMC::EDefinedSignal::kEta); + eta->SetFillPureMCStep(fillMC); + + + // eta dalitz + PairAnalysisSignalMC* etaDalitz = new PairAnalysisSignalMC(PairAnalysisSignalMC::EDefinedSignal::kEtaDalitz); + // etaDalitz->SetFillPureMCStep(fillMC); + etaDalitz->SetFillPureMCStep(GMC); + etaDalitz->SetIsDalitz(PairAnalysisSignalMC::EDalitz::kWhoCares, 0); + br = 6.9e-3; + etaDalitz->SetWeight(Meta * br); + + /// activate mc signal + /// IMPORTANT NOTE: single particle and pair signals are sorted, BUT later pair signals + /// should be the most detailed ones if you apply weights + + PairAnalysisSignalMC* xx = new PairAnalysisSignalMC("xx (comb.)", "xx"); + xx->SetFillPureMCStep(GMC); + xx->SetLegPDGs(0, 0); + xx->SetMothersRelation(PairAnalysisSignalMC::EBranchRelation::kDifferent); + + PairAnalysisSignalMC* ee = new PairAnalysisSignalMC("ee (comb.)", "ee"); + ee->SetFillPureMCStep(fillMC); + ee->SetLegPDGs(11, -11); + /// eleele->SetCheckBothChargesLegs(kTRUE,kTRUE); + ee->SetMothersRelation(PairAnalysisSignalMC::EBranchRelation::kDifferent); + + PairAnalysisSignalMC* epi = new PairAnalysisSignalMC("epi (comb.)", "epi"); + epi->SetFillPureMCStep(fillMC); + epi->SetLegPDGs(11, 211); + epi->SetCheckBothChargesLegs(kTRUE, kTRUE); + epi->SetMothersRelation(PairAnalysisSignalMC::EBranchRelation::kDifferent); + + PairAnalysisSignalMC* pipi = new PairAnalysisSignalMC("pipi (comb.)", "pipi"); + pipi->SetFillPureMCStep(fillMC); + pipi->SetLegPDGs(211, 211); + pipi->SetCheckBothChargesLegs(kTRUE, kTRUE); + pipi->SetMothersRelation(PairAnalysisSignalMC::EBranchRelation::kDifferent); + + + PairAnalysisSignalMC* inmed = new PairAnalysisSignalMC("in-medium SF", "in-medium SF"); + inmed->SetFillPureMCStep(GMC); + inmed->SetLegPDGs(11, -11); + inmed->SetMotherPDGs(99009011, 99009011); //0:default all + // inmed->SetMotherPDGs(99009911,99009911); //0:default all + inmed->SetMothersRelation(PairAnalysisSignalMC::EBranchRelation::kSame); + inmed->SetGEANTProcess(kPPrimary); + br = 4.45e-02 / 2; + inmed->SetWeight(Minmed); //default + //inmed->SetWeight(Mqgp);//default + + + // phi + PairAnalysisSignalMC* phi = new PairAnalysisSignalMC(PairAnalysisSignalMC::EDefinedSignal::kPhi); + phi->SetFillPureMCStep(GMC); + br = 2.97e-04; + phi->SetWeight(Mphi * br); + + + PairAnalysisSignalMC* qgp = new PairAnalysisSignalMC("QGP rad.", "QGP rad."); + qgp->SetFillPureMCStep(GMC); + qgp->SetLegPDGs(11, -11); + qgp->SetMotherPDGs(99009111, 99009111); //0:default all + // qgp->SetMotherPDGs(99009911,99009911); //0:default all + qgp->SetMothersRelation(PairAnalysisSignalMC::EBranchRelation::kSame); + qgp->SetGEANTProcess(kPPrimary); + br = 1.15e-02 / 2; + qgp->SetWeight(Mqgp); //default + + PairAnalysisSignalMC* Dp = new PairAnalysisSignalMC("D+", "D+"); + Dp->SetFillPureMCStep(GMC); + Dp->SetLegPDGs(-11, 11); + Dp->SetCheckBothChargesLegs(kTRUE, kTRUE); + Dp->SetMotherPDGs(2214, 2214); + Dp->SetMothersRelation(PairAnalysisSignalMC::EBranchRelation::kSame); + br = 4.2e-5; + Dp->SetWeight(MDp * br); + + // rho0 + PairAnalysisSignalMC* rho0 = new PairAnalysisSignalMC(PairAnalysisSignalMC::EDefinedSignal::kRho0); + rho0->SetFillPureMCStep(GMC); + br = 4.72e-05; + rho0->SetWeight(Mrho * br); //HSD + + + switch (cutDefinition) { + default: + /// Single particles + papa->AddSignalMC(ele); + papa->AddSignalMC(conv); + papa->AddSignalMC(pio); + papa->AddSignalMC(xx); + papa->AddSignalMC(ee); + papa->AddSignalMC(epi); + papa->AddSignalMC(pipi); + + papa->AddSignalMC(omegaDalitz); + omegaDalitz->Print(); + papa->AddSignalMC(pi0Dalitz); + pi0Dalitz->Print(); + papa->AddSignalMC(omega); + omega->Print(); + } +} + +//______________________________________________________________________________________ +void AddTrackHistograms(PairAnalysis* papa, Int_t cutDefinition) +{ + // + // add track histograms + // + + /// skip histograms in case of internal train + if (!papa->DoEventProcess()) return; + + /// skip certain configs + + PairAnalysisHistos* histos = papa->GetHistoManager(); + + /// Add track classes - LEGS of the pairs + if (!papa->IsNoPairing()) { + /// loop over all pair types and add classes (pair types in PairAnalysis.h EPairType) + /// automatically skip pair types w.r.t. configured bgrd estimators + PairAnalysisMixingHandler* mix = papa->GetMixingHandler(); + for (Int_t i = 0; i < static_cast<Int_t>(PairAnalysis::EPairType::kPairTypes); i++) { + switch (i) { + case static_cast<Int_t>(PairAnalysis::EPairType::kSEPP): + case static_cast<Int_t>(PairAnalysis::EPairType::kSEMM): + if (!papa->DoProcessLS()) { + continue; + break; + } + case static_cast<Int_t>(PairAnalysis::EPairType::kMEPP): + case static_cast<Int_t>(PairAnalysis::EPairType::kMEMM): + if (!mix || mix->GetMixType() != PairAnalysisMixingHandler::EMixType::kOSonly) { + continue; + break; + } + case static_cast<Int_t>(PairAnalysis::EPairType::kMEMP): + if (!mix) { + continue; + break; + } + case static_cast<Int_t>(PairAnalysis::EPairType::kMEPM): + if (!mix || mix->GetMixType() != PairAnalysisMixingHandler::EMixType::kAll) { + continue; + break; + } + case static_cast<Int_t>(PairAnalysis::EPairType::kSEPMRot): + if (!papa->GetTrackRotator()) { + continue; + break; + } + } + // histos->AddClass(Form("Track.Legs.%s",PairAnalysis::PairClassName(i))); + } + } + + /// Add track classes - single tracks used for any pairing + /// loop over all leg types and add classes (leg types in PairAnalysis.h ELegType) + + for (Int_t i = 0; i < static_cast<Int_t>(PairAnalysis::ELegType::kLegTypes); ++i) + histos->AddClass(Form("Track.%s", PairAnalysis::TrackClassName(i))); + + + /// add MC signal (if any) histograms to pair class + if (papa->GetMCSignals()) { + for (Int_t i = 0; i < papa->GetMCSignals()->GetEntriesFast(); ++i) { + PairAnalysisSignalMC* sigMC = (PairAnalysisSignalMC*) papa->GetMCSignals()->At(i); + + /// selection + if (!sigMC) continue; + //if(!sigMC->IsSingleParticle()) continue; /// skip pair particle signals (no pairs) + TString sigMCname = sigMC->GetName(); + /// by hand switched off + if (sigMCname.EqualTo("eleGamPiOmega")) continue; + + /// mc truth - pair leg class + // if(!papa->IsNoPairing() && sigMC->GetFillPureMCStep()) histos->AddClass(Form("Track.Legs_%s_MCtruth",sigMCname.Data())); + + /// mc reconstructed - pair leg class + // if(!papa->IsNoPairing()) histos->AddClass(Form("Track.Legs_%s", sigMCname.Data())); + + /// single tracks (merged +-) + histos->AddClass(Form("Track.%s_%s", + PairAnalysis::PairClassName(static_cast<Int_t>(PairAnalysis::EPairType::kSEPM)), + sigMCname.Data())); + if (sigMC->GetFillPureMCStep() && GMC) + histos->AddClass(Form("Track.%s_%s_MCtruth", + PairAnalysis::PairClassName(static_cast<Int_t>(PairAnalysis::EPairType::kSEPM)), + sigMCname.Data())); + } + } + + + histos->AddProfile("Track", PairAnalysisHelper::MakeLogBinning(200, 0., 20.), PairAnalysisVarManager::kP, + PairAnalysisVarManager::kPRes, "I"); + histos->AddHistogram("Track", PairAnalysisHelper::MakeLogBinning(200, 0., 20.), PairAnalysisVarManager::kP, + PairAnalysisHelper::MakeLogBinning(200, 0., 0.2), PairAnalysisVarManager::kPRes); + + + /// define MC and REC histograms + AddTrackHistogramsReconstruction(histos, cutDefinition); +} + +void AddTrackHistogramsReconstruction(PairAnalysisHistos* histos, Int_t cutDefinition) +{ + // + // add track histograms + // + + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(49, -0.5, 48.5), PairAnalysisVarManager::kGeantId); + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(100, 0., 500.), PairAnalysisVarManager::kNTrk); + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(400, 0, 4.), PairAnalysisVarManager::kPt); + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(400, 0, 20.), PairAnalysisVarManager::kP); + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(300, -3., 6.), PairAnalysisVarManager::kYlab, + PairAnalysisHelper::MakeLinBinning(125, 0, 5.), PairAnalysisVarManager::kPt); + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(300, -3., 6.), PairAnalysisVarManager::kYlab, + PairAnalysisHelper::MakeLinBinning(125, 0, 5.), PairAnalysisVarManager::kPt, + PairAnalysisHelper::MakeLinBinning(5, -0.5, 4.5), PairAnalysisVarManager::kTrdHits); + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(300, -3.0, 6.0), PairAnalysisVarManager::kYlab); + TVectorD* loM2 = PairAnalysisHelper::MakeLinBinning(400, -0.25, 0.15); + TVectorD* hiM2 = PairAnalysisHelper::MakeLinBinning(550, 0.15, 11.15); + + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(100, -0.5, 9.5), + PairAnalysisVarManager::kChi2NDFtoVtx); + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(100, -0.5, 9.5), + PairAnalysisVarManager::kTrdChi2NDF); + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(17, -0.5, 16.5), PairAnalysisVarManager::kStsHits); + // histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(21,-0.5, 20.5), PairAnalysisVarManager::kStsMvdHits); + //histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(17,-0.5, 16.5), "MvdHits+StsHits"); + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(6, -0.5, 5.5), PairAnalysisVarManager::kMvdHits); + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(103, -2.5, 100.5), PairAnalysisVarManager::kZvMC, + PairAnalysisHelper::MakeLinBinning(200, -100, 100), PairAnalysisVarManager::kRvMC); + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(510, -0.5, 100.5), PairAnalysisVarManager::kZvMC, + PairAnalysisHelper::MakeLinBinning(6, -0.5, 5.5), PairAnalysisVarManager::kMvdHits); + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(7, -0.5, 6.5), PairAnalysisVarManager::kTrdHits); + + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(3, -0.5, 2.5), PairAnalysisVarManager::kTofHits); + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(46, 4.5, 50.5), PairAnalysisVarManager::kRichHits); + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(110, -0.05, 1.05), + PairAnalysisVarManager::kTrdPidLikeEL); + histos->AddHistogram("Track", PairAnalysisHelper::MakeLogBinning(200, 0.1, 20.), PairAnalysisVarManager::kTrdPin, + PairAnalysisHelper::MakeLinBinning(110, -0.05, 1.05), PairAnalysisVarManager::kTrdPidLikeEL); + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(200, -2.0, 2.0), + PairAnalysisVarManager::kRichPidANN); + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(75, 0., 15.), PairAnalysisVarManager::kP, + PairAnalysisHelper::MakeLinBinning(200, -2.0, 2.0), PairAnalysisVarManager::kRichPidANN); + + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(200, 0., 1.), + PairAnalysisVarManager::kMvdHitClosest); + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(100, 0., 10.), PairAnalysisVarManager::kP, + PairAnalysisHelper::MakeLinBinning(400, 0., 4.), PairAnalysisVarManager::kMvdHitClosest); + histos->AddHistogram( + "Track", PairAnalysisHelper::MakeLinBinning(100, 0., 10.), PairAnalysisVarManager::kMvdHitClosestMom, + PairAnalysisHelper::MakeLinBinning(200, 0., TMath::Pi() / 4), PairAnalysisVarManager::kMvdHitClosestOpeningAngle); + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(100, 0., 10.), PairAnalysisVarManager::kP, + PairAnalysisHelper::MakeLinBinning(200, 0., TMath::Pi() / 4), + PairAnalysisVarManager::kMvdHitClosestOpeningAngle); + + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(200, 0., 1.), + PairAnalysisVarManager::kStsHitClosest); + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(100, 0., 10.), PairAnalysisVarManager::kP, + PairAnalysisHelper::MakeLinBinning(400, 0., 4.), PairAnalysisVarManager::kStsHitClosest); + histos->AddHistogram( + "Track", PairAnalysisHelper::MakeLinBinning(100, 0., 10.), PairAnalysisVarManager::kStsHitClosestMom, + PairAnalysisHelper::MakeLinBinning(200, 0., TMath::Pi() / 4), PairAnalysisVarManager::kStsHitClosestOpeningAngle); + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(100, 0., 10.), PairAnalysisVarManager::kP, + PairAnalysisHelper::MakeLinBinning(200, 0., TMath::Pi() / 4), + PairAnalysisVarManager::kStsHitClosestOpeningAngle); + + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(200, 0., TMath::Pi() / 4), + PairAnalysisVarManager::kMvdHitClosestOpeningAngle); + + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(7, -0.5, 6.5), PairAnalysisVarManager::kTrdHitsMC); + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(7, -0.5, 6.5), PairAnalysisVarManager::kStsHitsMC); + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(7, -0.5, 6.5), PairAnalysisVarManager::kMvdHitsMC); + + histos->AddHistogram("Track", PairAnalysisHelper::MakePdgBinning(), PairAnalysisVarManager::kPdgCode); + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(49, -0.5, 48.5), PairAnalysisVarManager::kGeantId); + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(50, 0., 10.), PairAnalysisVarManager::kE); + + // histos->AddHistogram("Track", PairAnalysisHelper::MakePdgBinning(), PairAnalysisVarManager::kPdgCodeMother); + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(10, -0.5, 9.5), + PairAnalysisVarManager::kPdgCodeMother); + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(10, -0.5, 9.5), + PairAnalysisVarManager::kPdgCodeMother, PairAnalysisHelper::MakeLinBinning(12, -0.5, 11.5), + PairAnalysisVarManager::kGeantId); + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(10, -0.5, 9.5), + PairAnalysisVarManager::kPdgCodeMother, PairAnalysisHelper::MakeLinBinning(12, -0.5, 11.5), + PairAnalysisVarManager::kGeantId, PairAnalysisVarManager::kWeight); + + // histos->AddHistogram("Track",PairAnalysisHelper::MakeLinBinning(510,-0.5,100.), PairAnalysisVarManager::kZvMC, + // PairAnalysisHelper::MakeLinBinning(500,-10,100.), PairAnalysisVarManager::kRvMC); + + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(120, 0., 30.), + PairAnalysisVarManager::kMvdFirstHitPosZ); + + histos->AddHistogram("Track", PairAnalysisHelper::MakeLogBinning(100, 0., 10.), PairAnalysisVarManager::kP, + PairAnalysisHelper::MakeLinBinning(1000, -0.5, 0.5), PairAnalysisVarManager::kTofPidDeltaBetaEL); +} + + +//______________________________________________________________________________________ +void AddPairHistograms(PairAnalysis* papa, Int_t cutDefinition) +{ + /// + /// add pair histograms + /// + + /// skip if no pairing done + if (papa->IsNoPairing()) return; + + PairAnalysisHistos* histos = papa->GetHistoManager(); + + /// add histogram classes + /// loop over all pair types and add classes (pair types in PairAnalysis.h EPairType) + /// automatically skip pair types w.r.t. configured bgrd estimators + PairAnalysisMixingHandler* mix = papa->GetMixingHandler(); + for (Int_t i = 0; i < static_cast<Int_t>(PairAnalysis::EPairType::kPairTypes); i++) { + switch (i) { + case static_cast<Int_t>(PairAnalysis::EPairType::kSEPP): + case static_cast<Int_t>(PairAnalysis::EPairType::kSEMM): + if (!papa->DoProcessLS()) { + continue; + break; + } + case static_cast<Int_t>(PairAnalysis::EPairType::kMEPP): + case static_cast<Int_t>(PairAnalysis::EPairType::kMEMM): + if (!mix || mix->GetMixType() != PairAnalysisMixingHandler::EMixType::kOSonly) { + continue; + break; + } + case static_cast<Int_t>(PairAnalysis::EPairType::kMEMP): + if (!mix) { + continue; + break; + } + case static_cast<Int_t>(PairAnalysis::EPairType::kMEPM): + if (!mix || mix->GetMixType() != PairAnalysisMixingHandler::EMixType::kAll) { + continue; + break; + } + case static_cast<Int_t>(PairAnalysis::EPairType::kSEPMRot): + if (!papa->GetTrackRotator()) { + continue; + break; + } + } + histos->AddClass(Form("Pair.%s", PairAnalysis::PairClassName(i))); + } + + + /// mixing statistics + Int_t mixBins = 0; + if (mix) { + mixBins = mix->GetNumberOfBins(); + histos->AddHistogram("Pair", PairAnalysisHelper::MakeLinBinning(mixBins, 0., mixBins), + PairAnalysisVarManager::kMixingBin); + histos->AddHistogram("Pair", PairAnalysisHelper::MakeLinBinning(100, 0., 400.), + PairAnalysisVarManager::kNVtxContrib); + } + + ////// add MC signal histo classes + if (papa->GetMCSignals()) { + for (Int_t i = 0; i < papa->GetMCSignals()->GetEntriesFast(); ++i) { + PairAnalysisSignalMC* sigMC = (PairAnalysisSignalMC*) papa->GetMCSignals()->At(i); + if (!sigMC) continue; + if (sigMC->IsSingleParticle()) continue; /// skip pair particle signals (no pairs) + TString sigMCname = sigMC->GetName(); + histos->AddClass(Form("Pair_%s", sigMCname.Data())); + if (sigMC->GetFillPureMCStep() && GMC) histos->AddClass(Form("Pair_%s_MCtruth", sigMCname.Data())); + } + } + + + ///// define output objects for MC and REC ///// + histos->AddHistogram("Pair", PairAnalysisHelper::MakeLinBinning(300, 0., 3.), + PairAnalysisVarManager::kM); /// 20MeV bins, 5 GeV/c2 + histos->AddHistogram("Pair", PairAnalysisHelper::MakeLinBinning(100, 0., 10.), + PairAnalysisVarManager::kP); /// 20MeV bins, 5 GeV/2c + histos->AddHistogram("Pair", PairAnalysisHelper::MakeLinBinning(300, 0., 3.), PairAnalysisVarManager::kM, + PairAnalysisVarManager::kWeight); // 40MeV bins, 12GeV/c2 + // histos->AddHistogram("Pair", PairAnalysisHelper::MakeLinBinning(300,0.,3.), PairAnalysisVarManager::kM, + // PairAnalysisHelper::MakeLinBinning(1000,0.000000000001,0.0001),PairAnalysisVarManager::kWeight); // 40MeV bins, 12GeV/c2 + histos->AddHistogram("Pair", PairAnalysisHelper::MakeLinBinning(250, 0, 5.), PairAnalysisVarManager::kPt); + histos->AddHistogram("Pair", PairAnalysisHelper::MakeLinBinning(200, -2., 6.), PairAnalysisVarManager::kYlab); + histos->AddHistogram("Pair", PairAnalysisHelper::MakeLinBinning(200, -2., 6.), PairAnalysisVarManager::kYlab, + PairAnalysisHelper::MakeLinBinning(250, 0, 5.), PairAnalysisVarManager::kPt); + + + histos->AddHistogram("Pair", PairAnalysisHelper::MakeLinBinning(49, -0.5, 48.5), PairAnalysisVarManager::kGeantId); + histos->AddHistogram("Pair", PairAnalysisHelper::MakeLinBinning(10, -0.5, 9.5), + PairAnalysisVarManager::kPdgCodeMother); + // histos->AddHistogram("Pair", PairAnalysisHelper::MakeLinBinning(49,-0.5,48.5), PairAnalysisVarManager::kGeantIdMother); + + // histos->AddHistogram("Pair",PairAnalysisHelper::MakeLinBinning(900,0.,900.), PairAnalysisVarManager::kZvMC, + // PairAnalysisHelper::MakeLinBinning(1000,0.,5000.), PairAnalysisVarManager::kRvMC); + + // histos->AddHistogram("Pair",PairAnalysisHelper::MakeLinBinning(510,-0.5,100.5), PairAnalysisVarManager::kZvMC, + // PairAnalysisHelper::MakeLinBinning(100,0.,.5), PairAnalysisVarManager::kRvMC); + + histos->AddHistogram("Pair", PairAnalysisHelper::MakeLinBinning(21, -0.5, 20.5), + PairAnalysisVarManager::kStsMvdFirstDaughter, PairAnalysisHelper::MakeLinBinning(21, -0.5, 20.5), + PairAnalysisVarManager::kStsMvdSecondDaughter); + histos->AddHistogram( + "Pair", PairAnalysisHelper::MakeLinBinning(21, -0.5, 20.5), PairAnalysisVarManager::kStsMvdTrdFirstDaughter, + PairAnalysisHelper::MakeLinBinning(21, -0.5, 20.5), PairAnalysisVarManager::kStsMvdTrdSecondDaughter); + histos->AddHistogram( + "Pair", PairAnalysisHelper::MakeLinBinning(71, -0.5, 70.5), PairAnalysisVarManager::kStsMvdRichFirstDaughter, + PairAnalysisHelper::MakeLinBinning(71, -0.5, 70.5), PairAnalysisVarManager::kStsMvdRichSecondDaughter); + histos->AddHistogram("Pair", PairAnalysisHelper::MakeLinBinning(21, -0.5, 20.5), + PairAnalysisVarManager::kStsFirstDaughter, PairAnalysisHelper::MakeLinBinning(21, -0.5, 20.5), + PairAnalysisVarManager::kStsSecondDaughter); + histos->AddHistogram("Pair", PairAnalysisHelper::MakeLinBinning(21, -0.5, 20.5), + PairAnalysisVarManager::kMvdFirstDaughter, PairAnalysisHelper::MakeLinBinning(21, -0.5, 20.5), + PairAnalysisVarManager::kMvdSecondDaughter); + + histos->AddHistogram("Pair", PairAnalysisHelper::MakeLinBinning(200, 0., 10.), PairAnalysisVarManager::kMvdHitDist); + histos->AddHistogram("Pair", PairAnalysisHelper::MakeLinBinning(100, 0., 10.), PairAnalysisVarManager::kP, + PairAnalysisHelper::MakeLinBinning(200, 0., 10.), PairAnalysisVarManager::kMvdHitDist); + + histos->AddHistogram("Pair", PairAnalysisHelper::MakeLinBinning(400, 0., 40.), PairAnalysisVarManager::kStsHitDist); + histos->AddHistogram("Pair", PairAnalysisHelper::MakeLinBinning(100, 0., 10.), PairAnalysisVarManager::kP, + PairAnalysisHelper::MakeLinBinning(400, 0., 40.), PairAnalysisVarManager::kStsHitDist); + + + histos->AddHistogram("Pair", PairAnalysisHelper::MakeLinBinning(100, 0., TMath::Pi() / 2), + PairAnalysisVarManager::kOpeningAngle); + histos->AddHistogram("Pair", PairAnalysisHelper::MakeLinBinning(100, 0., 5.), PairAnalysisVarManager::kLegsP); + histos->AddHistogram("Pair", PairAnalysisHelper::MakeLinBinning(314, 0, 3.14), PairAnalysisVarManager::kPhivPair); + histos->AddHistogram("Pair", PairAnalysisHelper::MakeLinBinning(314, 0, 3.14), + PairAnalysisVarManager::kCosPointingAngleMC); + histos->AddHistogram("Pair", PairAnalysisHelper::MakeLinBinning(314, 0, 3.14), + PairAnalysisVarManager::kCosPointingAngle); + histos->AddHistogram("Pair", PairAnalysisHelper::MakeLinBinning(100, 0., 5.), PairAnalysisVarManager::kLegsP, + PairAnalysisHelper::MakeLinBinning(100, 0., 0.5), PairAnalysisVarManager::kOpeningAngle); + histos->AddHistogram("Pair", PairAnalysisHelper::MakeLinBinning(100, 0., 5.), PairAnalysisVarManager::kLegsP, + PairAnalysisHelper::MakeLinBinning(314, 0., 3.14), PairAnalysisVarManager::kPhivPair); + histos->AddHistogram("Pair", PairAnalysisHelper::MakeLinBinning(200, 0., TMath::Pi() / 2), + PairAnalysisVarManager::kOpeningAngle, PairAnalysisHelper::MakeLinBinning(314, 0., 3.14), + PairAnalysisVarManager::kPhivPair); + // histos->AddHistogram("Pair", PairAnalysisHelper::MakeLinBinning(200,-1.,1.), "0.035-LegsP*0.01-OpeningAngle"); +} + + +//______________________________________________________________________________________ +void AddHitHistograms(PairAnalysis* papa, Int_t cutDefinition) +{ + /// + /// add hit histograms + /// + + PairAnalysisHistos* histos = papa->GetHistoManager(); + + /// Add hit classes - TRD tracks (before pairing) + histos->AddClass("Hit.TRD"); + + + /// add MC signal (if any) histograms to hit class + if (papa->GetMCSignals()) { + for (Int_t i = 0; i < papa->GetMCSignals()->GetEntriesFast(); ++i) { + PairAnalysisSignalMC* sigMC = (PairAnalysisSignalMC*) papa->GetMCSignals()->At(i); + if (!sigMC) continue; + if (!sigMC->IsSingleParticle()) continue; /// skip pair particle signals (no pairs) + TString sigMCname = sigMC->GetName(); + /// by hand switched off + if (sigMCname.EqualTo("eleGamPiOmega")) continue; + + /// DET specific MC histograms + histos->AddClass(Form("Hit.TRD_%s", sigMCname.Data())); + + + histos->AddHistogram(Form("Hit.TRD_%s", sigMCname.Data()), PairAnalysisHelper::MakeLinBinning(400, 0., 1.e+2), + PairAnalysisVarManager::kEloss); + histos->AddHistogram(Form("Hit.TRD_%s", sigMCname.Data()), PairAnalysisHelper::MakeLinBinning(400, 0., 1.e+2), + PairAnalysisVarManager::kElossMC); + histos->AddHistogram(Form("Hit.TRD_%s", sigMCname.Data()), + PairAnalysisHelper::MakeArbitraryBinning("0.,0.25,0.5,1.,1.5,2.,3.,4.,5.,6.,7.,8.,9.,10."), + PairAnalysisVarManager::kTrdPin, PairAnalysisHelper::MakeLinBinning(400, 0., 1.e+2), + PairAnalysisVarManager::kEloss); + + histos->AddHistogram(Form("Hit.TRD_%s", sigMCname.Data()), PairAnalysisHelper::MakeLinBinning(16, -0.5, 15.5), + PairAnalysisVarManager::kTrdPads, PairAnalysisHelper::MakeLinBinning(10, -0.5, 9.5), + PairAnalysisVarManager::kLinksMC); + histos->AddProfile(Form("Hit.TRD_%s", sigMCname.Data()), PairAnalysisHelper::MakeLinBinning(16, -0.5, 15.5), + PairAnalysisVarManager::kTrdPads, PairAnalysisVarManager::kLinksMC, "I"); + + histos->AddHistogram(Form("Hit.TRD_%s", sigMCname.Data()), PairAnalysisHelper::MakeLinBinning(16, -0.5, 15.5), + PairAnalysisVarManager::kTrdPads); + histos->AddHistogram(Form("Hit.TRD_%s", sigMCname.Data()), PairAnalysisHelper::MakeLinBinning(7, -0.5, 6.5), + PairAnalysisVarManager::kTrdLayer); + } + } + + + /// histograms + histos->AddHistogram("Hit", PairAnalysisHelper::MakeLinBinning(200, 0., 5.e+1), PairAnalysisVarManager::kEloss); + histos->AddHistogram("Hit", + PairAnalysisHelper::MakeArbitraryBinning("0.,0.25,0.5,1.,1.5,2.,3.,4.,5.,6.,7.,8.,9.,10.,20."), + PairAnalysisVarManager::kTrdPin, PairAnalysisHelper::MakeLinBinning(800, 0., 8.e+1), + PairAnalysisVarManager::kEloss); + histos->AddHistogram( + "Hit", PairAnalysisHelper::MakeArbitraryBinning("0.,0.25,0.5,1.,1.5,2.,3.,4.,5.,6.,7.,8.,9.,10.,20."), + PairAnalysisVarManager::kP, PairAnalysisHelper::MakeLinBinning(800, 0., 8.e+1), PairAnalysisVarManager::kEloss); + + /// mc matching + histos->AddHistogram("Hit", PairAnalysisHelper::MakeLinBinning(200, 0., 5.e+1), PairAnalysisVarManager::kElossMC); + histos->AddHistogram("Hit", PairAnalysisHelper::MakeLogBinning(200, 0.1, 20.), PairAnalysisVarManager::kTrdPin, + PairAnalysisHelper::MakeLinBinning(300, 0., 3.e+1), PairAnalysisVarManager::kElossMC); + + histos->AddHistogram("Hit", PairAnalysisHelper::MakeLinBinning(400, -400., 400.), PairAnalysisVarManager::kPosX, + PairAnalysisHelper::MakeLinBinning(400, -400., 400.), PairAnalysisVarManager::kPosY); +} + + +//______________________________________________________________________________________ +void AddCutStepHistograms(PairAnalysis* papa, Int_t cutDefinition) +{ + // + // add QA histograms for each track cut applied in the analysis + // NOTE: this is rather CPU expencive (at least for MCtruth) since + // there signal check are reperformed + // + + /// skip certain configs + // //if(cutDefinition!=kRichTRDTOFcfg) return; + // if(cutDefinition!=kTrdcfg) return; + + PairAnalysisHistos* histos = new PairAnalysisHistos(); + + /// add MC signal (if any) histograms to pair class + if (papa->GetMCSignals()) { + for (Int_t i = 0; i < papa->GetMCSignals()->GetEntriesFast(); ++i) { + PairAnalysisSignalMC* sigMC = (PairAnalysisSignalMC*) papa->GetMCSignals()->At(i); + + /// selection + if (!sigMC) continue; + TString sigMCname = sigMC->GetName(); + // if(!sigMCname.Contains("PrimEle") && !sigMCname.Contains("PrimPio")) continue; /// skip pair particle signals (no pairs) + /// single tracks (merged +-) + histos->AddClass(Form("Track.%s_%s", + PairAnalysis::PairClassName(static_cast<Int_t>(PairAnalysis::EPairType::kSEPM)), + sigMCname.Data())); + // histos->AddClass(Form("Track.Legs.%s_%s",PairAnalysis::PairClassName(static_cast<Int_t>(PairAnalysis::EPairType::kSEPM)),sigMCname.Data())); + // histos->AddClass(Form("Pair_%s",sigMCname.Data())); + // if(sigMC->GetFillPureMCStep()) + // histos->AddClass(Form("Track.%s_%s_MCtruth",PairAnalysis::PairClassName(PairAnalysis::kSEPM),sigMCname.Data())); + } + } + + /// histograms + histos->AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(100, 0.0, TMath::Pi() / 4), + PairAnalysisVarManager::kTrdThetain); + histos->AddHistogram("Track", PairAnalysisHelper::MakeLogBinning(400, 0, 20.), PairAnalysisVarManager::kP); + histos->AddHistogram("Track", PairAnalysisHelper::MakeLogBinning(250, 0, 2.5), PairAnalysisVarManager::kM); + + papa->GetTrackFilter().AddHistos(histos); //prefilter histograms + papa->GetPairPreFilterLegs().AddHistos(histos); //prefilter histograsm + papa->GetFinalTrackFilter().AddHistos(histos); +} diff --git a/analysis/PWGDIL/dielectron/papaframework/scripts/ana.sh b/analysis/PWGDIL/dielectron/papaframework/scripts/ana.sh new file mode 100755 index 0000000000000000000000000000000000000000..ad001c2135c69f38be225274bb346dae0615cf4d --- /dev/null +++ b/analysis/PWGDIL/dielectron/papaframework/scripts/ana.sh @@ -0,0 +1,68 @@ +#!/bin/bash + +LOCATION=/lustre/nyx + +## choose cbm root installation +#source $LOCATION/cbm/users/$USER/CBMsoft/CbmRoot/trunk/build/config.sh +#source $LOCATION/cbm/users/$USER/CBMsoft/cbm-env.sh -n 4 + +## job content +export INDIR="$1" ## path to simreco directory +## number of events to process +NEVT=$2 +CONFIG=$3 ## analysis mode +GROUP=$4 ## grouping of files +TESTBIN="$5" +NAME=$6 + +export INFILE="$1" ## path to input file list for this jobs +export OUTDIR="$2"## output directory + +I=0 +n=0 +while [ "$n" -lt "$GROUP" ]; do + echo "file $n $8" + filelist[$n]=$7 + shift + let n=$n+1 +done +n=0 +while [ "$n" -lt "$GROUP" ]; do + echo "out $n $8" + outlist[$n]=$7 + shift + let n=$n+1 + +done + +echo "outlist ${outlist[*]}" +echo "filelist ${filelist[*]}" + +while [ "$I" -lt "$GROUP" ]; do + export OUTDIR="${outlist[$I]}" + export INFILE="${filelist[$I]}" ## path to simreco directory + + if [[ $CONFIG == "1" ]] ; then + if [[ $TESTBIN == "0" ]] ; then + echo "$OUTDIR/../run_common.C($NEVT)" + ${ROOTSYS}/bin/root -l -b -q "$OUTDIR/../run_testing.C($NEVT,kFALSE,$NAME)" + else + echo "$OUTDIR/../run_common.C($NEVT)" + ${ROOTSYS}/bin/root -l -b -q "$OUTDIR/../run_testing.C($NEVT,kTRUE,$NAME)" + fi + fi + if [[ $CONFIG == "2" ]] ; then + if [[ $TESTBIN == "0" ]] ; then + echo "$OUTDIR/../run_common_analysis.C($NEVT)" + ${ROOTSYS}/bin/root -l -b -q "$OUTDIR/../run_common_analysis.C($NEVT,kFALSE)" + else + echo "$OUTDIR/../run_common_analysis.C($NEVT)" + ${ROOTSYS}/bin/root -l -b -q "$OUTDIR/../run_common_analysis.C($NEVT,kTRUE)" + fi + fi + + ## cleanup + rm -v L1_histo.root + + let I=$I+1 +done diff --git a/analysis/PWGDIL/dielectron/papaframework/scripts/analysis.sh b/analysis/PWGDIL/dielectron/papaframework/scripts/analysis.sh index 6f55eb4fd86a350b882446f8a838791288f9c7b3..65fc9f820eea9c77925a872093e4329069ce0d55 100755 --- a/analysis/PWGDIL/dielectron/papaframework/scripts/analysis.sh +++ b/analysis/PWGDIL/dielectron/papaframework/scripts/analysis.sh @@ -3,17 +3,42 @@ LOCATION=/lustre/nyx ## choose cbm root installation -source $LOCATION/cbm/users/$USER/CBMsoft/cbm-env.sh -n +#source $LOCATION/cbm/users/$USER/CBMsoft/CbmRoot/trunk/build/config.sh +#source $LOCATION/cbm/users/$USER/CBMsoft/cbm-env.sh -n 4 ## job content export INDIR="$1" ## path to simreco directory export INFILE="$2" ## path to input file list for this jobs export OUTDIR="$3" ## output directory + ## number of events to process NEVT=$4 +CONFIG=$5 ## analysis mode +export TESTBIN="$6" +NAME=$7 -root -l -b -q "$OUTDIR/../run_analysis.C($NEVT)" +if [[ $CONFIG == "1" ]] ; then + if [[ $TESTBIN == "1" ]] ; then + echo "$OUTDIR/../run_testing.C($NEVT)" + ${ROOTSYS}/bin/root -l -b -q "$OUTDIR/../run_testing.C($NEVT,kTRUE,$NAME)" + fi + if [[ $TESTBIN == "0" ]] ; then + echo "$OUTDIR/../run_testing.C($NEVT)" + ${ROOTSYS}/bin/root -l -b -q "$OUTDIR/../run_testing.C($NEVT,kFALSE,$NAME)" + fi +fi +if [[ $CONFIG == "2" ]] ; then + if [[ $TESTBIN == "0" ]] ; then + echo "$OUTDIR/../run_common_analysis.C($NEVT)" + ${ROOTSYS}/bin/root -l -b -q "$OUTDIR/../run_common_analysis.C($NEVT,kFALSE)" + else + echo "$OUTDIR/../run_common_analysis.C($NEVT)" + cp -v "$OUTDIR/../../par.root" "$OUTDIR/par.root" + ${ROOTSYS}/bin/root -l -b -q "$OUTDIR/../run_common_analysis.C($NEVT,kTRUE)" + fi +fi +#${ROOTSYS}/bin/root -l -b -q "$OUTDIR/../run_analysis_old.C($NEVT)" ## cleanup rm -v L1_histo.root diff --git a/analysis/PWGDIL/dielectron/papaframework/scripts/inputfile_AuAu_8 b/analysis/PWGDIL/dielectron/papaframework/scripts/inputfile_AuAu_8 deleted file mode 100644 index f0616df3222c4ff6df826992b9a24db9d0c31c1f..0000000000000000000000000000000000000000 --- a/analysis/PWGDIL/dielectron/papaframework/scripts/inputfile_AuAu_8 +++ /dev/null @@ -1,18 +0,0 @@ -pro 197 79 -tar 197 79 - -nev 1000 -imp 0. - -plb 8 -tim 100 100 -rsd 1 - -f13 -#f14 -f15 -f16 -f19 -f20 - -xxx diff --git a/analysis/PWGDIL/dielectron/papaframework/scripts/inputfile_AuAu_8_010 b/analysis/PWGDIL/dielectron/papaframework/scripts/inputfile_AuAu_8_010 deleted file mode 100644 index d70e4e9d45a08a3ed3e0503601fa829b698dbc9d..0000000000000000000000000000000000000000 --- a/analysis/PWGDIL/dielectron/papaframework/scripts/inputfile_AuAu_8_010 +++ /dev/null @@ -1,18 +0,0 @@ -pro 197 79 -tar 197 79 - -nev 1000 -imp 0. 4.65 - -plb 8 -tim 100 100 -rsd 1 - -f13 -#f14 -f15 -f16 -f19 -f20 - -xxx diff --git a/analysis/PWGDIL/dielectron/papaframework/scripts/inputfile_pAu_30 b/analysis/PWGDIL/dielectron/papaframework/scripts/inputfile_pAu_30 deleted file mode 100644 index 80b610f9de5eada1410892e5a6a10fcd3cab45a5..0000000000000000000000000000000000000000 --- a/analysis/PWGDIL/dielectron/papaframework/scripts/inputfile_pAu_30 +++ /dev/null @@ -1,18 +0,0 @@ -pro 1 1 -tar 197 79 - -nev 1000 -imp 0. - -plb 30 -tim 100 100 -rsd 1 - -f13 -#f14 -f15 -f16 -f19 -f20 - -xxx diff --git a/analysis/PWGDIL/dielectron/papaframework/scripts/merge.C b/analysis/PWGDIL/dielectron/papaframework/scripts/merge.C index 4ec70adef259ca04933af3d51329448d8dd3726f..ee6599205b3bdb65aedd219f2c8b5548d327fb2c 100644 --- a/analysis/PWGDIL/dielectron/papaframework/scripts/merge.C +++ b/analysis/PWGDIL/dielectron/papaframework/scripts/merge.C @@ -1,4 +1,5 @@ -int merge(char* input_list = "list.txt", char* output_file = "merged.root") { +int merge(char* input_list = "list.txt", char* output_file = "merged.root") +{ /// set global debug level // gErrorIgnoreLevel = kInfo, kWarning, kError, kFatal; @@ -30,15 +31,14 @@ int merge(char* input_list = "list.txt", char* output_file = "merged.root") { /// was merge properly, you can savely ignore it (this results in job status FAILED) if (1) { /// selection/rejection of objects - FM->AddObjectNames( - "cbmout BranchList FileHeader cbmsim TimeBasedBranchList"); + FM->AddObjectNames("cbmout BranchList FileHeader cbmsim TimeBasedBranchList"); // Must add new merging flag on top of the the default ones - Int_t mode = TFileMerger::kAllIncremental - | TFileMerger::kSkipListed; /// skip ObjectNames + Int_t mode = TFileMerger::kAllIncremental | TFileMerger::kSkipListed; /// skip ObjectNames // TFileMerger::kAllIncremental | TFileMerger::kOnlyListed; /// merge only ObjectNames FM->PartialMerge(mode); - } else { + } + else { /// default merging FM->Merge(); } diff --git a/analysis/PWGDIL/dielectron/papaframework/scripts/run_analysis.C b/analysis/PWGDIL/dielectron/papaframework/scripts/run_analysis.C deleted file mode 100644 index 63af02d0db456bc3b79e38fe1e6860697e58728a..0000000000000000000000000000000000000000 --- a/analysis/PWGDIL/dielectron/papaframework/scripts/run_analysis.C +++ /dev/null @@ -1,216 +0,0 @@ -FairTask* addpapa(TString configPath, TString configFunct, TString configName) { - - TString myName = "run_analysis"; - std::cout << "-I- " << myName << ": Loading task " - << (configFunct + configName) << std::endl; - - TString cfgFile = configPath + configFunct + configName + ".C"; - // std::cout << "-I- " << myName << ": Loading macro " << cfgFile << std::endl; - gROOT->LoadMacro((configPath + configFunct + configName + ".C")); - - TString cfgFunct = configFunct + configName + "(\"" + configName + "\")"; - // std::cout << "-I- " << myName << ": process line " << cfgFunct << std::endl; - return (dynamic_cast<FairTask*>(gROOT->ProcessLine(cfgFunct))); -} - -void run_analysis(Int_t nEvents = 0) { - - // ----- Environment -------------------------------------------------- - TString myName = "run_analysis"; // this macro's name for screen output - TString srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory - // ------------------------------------------------------------------------ - - // ----- In- and output directory names ------------------------------- - TString inDir = gSystem->Getenv("INDIR"); - TString inFile = gSystem->Getenv("INFILE"); // input file list - TString outDir = gSystem->Getenv("OUTDIR"); - if (outDir.IsNull()) outFile = "."; - // ------------------------------------------------------------------------ - - // --- Logger settings ---------------------------------------------------- - gErrorIgnoreLevel = kFatal; //kInfo, kWarning, kError, kFatal; - // ------------------------------------------------------------------------ - - // ----- Load the geometry setup ------------------------------------- - /// load local copy of setup file (same as used for simulation and reconstruction) - std::cout << std::endl; - TString setupName = gSystem->GetFromPipe( - Form("echo $(basename $(ls %s/setup_*.C))", inDir.Data())); - setupName.ReplaceAll("setup_", ""); - setupName.ReplaceAll(".C", ""); - - TString setupFile = inDir + "/setup_" + setupName + ".C"; - TString setupFunct = "setup_"; - setupFunct = setupFunct + setupName + "()"; - std::cout << "-I- " << myName << ": Loading macro " << setupFile << std::endl; - gROOT->LoadMacro(setupFile); - gROOT->ProcessLine(setupFunct); - CbmSetup* setup = CbmSetup::Instance(); - // ------------------------------------------------------------------------ - - - // ----- Modify Cbm setup --------------------------------------------- - std::cout << std::endl; - // --- remove detector geometries - // setup->RemoveModule(kMvd); // remove mvd and its material budget - // setup->RemoveModule(kTrd); // e.g. for sts-tof-matching study - setup->RemoveModule(kPsd); // remove psd from setup - // --- change default geomerties - // setup->SetModule(kTrd, "v15d_1e", kTRUE); // 5 TRD layer - std::cout << "-I- " << myName << ": CbmSetup updated " << std::endl; - // ------------------------------------------------------------------------ - - - // ----- Run manager + I/O -------------------------------------------- - std::cout << std::endl; - FairRunAna* run = new FairRunAna(); - run->SetOutputFile(outDir + "/" + setupName + "_analysis.root"); - FairFileSource* src = NULL; - - /// stopwatch - TStopwatch timer; - timer.Start(); - - Int_t i = 0; - TString file; - // char filename[300]; - ifstream in(inFile); - TList* parList = new TList(); - /// loop over all file in list - while (in >> file) { - // TString file(filename);// + "/" + setupName; - file += "/"; - file += setupName; - - // mc sim file - if (!i) - src = new FairFileSource(file + "_mc.root"); - else - src->AddFile(file + "_mc.root"); - - // (skimmed) reco file - src->AddFriend(file + "_reco.root"); - - // parameter files - TObjString* parFile = new TObjString((file + "_params.root").Data()); - parList->Add(parFile); - // parIo1->open(parFile.GetString().Data()); - - i++; - } - - // add source to run - run->SetSource(src); - - // set parameter list - FairParRootFileIo* parIo1 = new FairParRootFileIo(); - parIo1->open(parList); - std::cout << "-I- " << myName << ": in/output added " << std::endl; - // ------------------------------------------------------------------------ - - - // ----- L1/KF tracking + PID ----------------------------------------- - std::cout << std::endl; - std::cout << "-I- " << myName << ": Loading tasks " << std::endl; - - //CbmKF is needed for Extrapolation - CbmKF* kf = new CbmKF(); - run->AddTask(kf); - std::cout << "-I- : Added task " << kf->GetName() << std::endl; - - CbmL1* l1 = new CbmL1(); - // --- Material budget file names - if (setup->IsActive(kMvd)) { - TString geoTag; - setup->GetGeoTag(kMvd, geoTag); - TString matFile = gSystem->Getenv("VMCWORKDIR"); - matFile = matFile + "/parameters/mvd/mvd_matbudget_" + geoTag + ".root"; - std::cout << "Using material budget file " << matFile << std::endl; - l1->SetMvdMaterialBudgetFileName(matFile.Data()); - } - if (setup->IsActive(kSts)) { - TString geoTag; - setup->GetGeoTag(kSts, geoTag); - TString matFile = gSystem->Getenv("VMCWORKDIR"); - matFile = matFile + "/parameters/sts/sts_matbudget_" + geoTag + ".root"; - std::cout << "Using material budget file " << matFile << std::endl; - l1->SetStsMaterialBudgetFileName(matFile.Data()); - } - run->AddTask(l1); - std::cout << "-I- : Added task " << l1->GetName() << std::endl; - - // --- TRD pid tasks - if (setup->IsActive(kTrd)) { - // ----------- TRD track Pid Ann ---------------------- - CbmTrdSetTracksPidANN* trdSetTracksPidAnnTask = - new CbmTrdSetTracksPidANN("Ann", "Ann"); - run->AddTask(trdSetTracksPidAnnTask); - std::cout << "-I- : Added task " << trdSetTracksPidAnnTask->GetName() - << std::endl; - // ---------------------------------------------------- - - // ----------- TRD track Pid Like ---------------------- - FairTask* getTRDli = dynamic_cast<FairTask*>( - gROOT->Macro("$VMCWORKDIR/macro/papa/Config_trd_PidLI.C")); - run->AddTask(getTRDli); - std::cout << "-I- : Added task " << getTRDli->GetName() << std::endl; - - // CbmTrdSetTracksPidLike* trdLI = new CbmTrdSetTracksPidLike("TRDLikelihood", "TRDLikelihood"); - // trdLI->SetUseMCInfo(kTRUE); - // trdLI->SetUseMomDependence(kFALSE); - // run->AddTask(trdLI); - // std::cout << "-I- : Added task " << trdLI->GetName() << std::endl; - // ------------------------------------------------------------------------ - } - - // ----- PAPA tasks --------------------------------------------------- - std::cout << std::endl; - std::cout << "-I- " << myName << ": Loading private tasks " << std::endl; - - // gSystem->Exec(Form("cp -u %s/*.C $OUTDIR/../.",configPath.Data())); - TString cfgPath = outDir + "/../"; - TString cfgFunc = "Config_jbook_"; - - // run->AddTask( addpapa(cfgPath,cfgFunc, "MC") ); - // run->AddTask( addpapa(cfgPath,cfgFunc, "URQMD") ); - // run->AddTask( addpapa(cfgPath,cfgFunc, "MVD") ); - // run->AddTask( addpapa(cfgPath,cfgFunc, "MUCH") ); - // run->AddTask( addpapa(cfgPath,cfgFunc, "TOF") ); - - // run->AddTask( addpapa(cfgPath,cfgFunc, "QA") ); - // run->AddTask( addpapa(cfgPath,cfgFunc, "StsQA") ); - - //run->AddTask( addpapa(cfgPath,cfgFunc, "AA") ); - // run->AddTask( addpapa(cfgPath,cfgFunc, "pA") ); - // run->AddTask( addpapa(cfgPath,cfgFunc, "RF") ); - - // run->AddTask( addtask(cfgPath,cfgFunc, "PID") ); - // run->AddTask( addtask(cfgPath,cfgFunc, "Perf") ); - // run->AddTask( addtask(cfgPath,cfgFunc, "Phi") ); - - // run->AddTask( addpapa(cfgPath,cfgFunc, "AAfast") ); - // run->AddTask( addpapa(cfgPath,cfgFunc, "RFfast") ); - // run->AddTask( addpapa(cfgPath,cfgFunc, "Test") ); - - // ------------------------------------------------------------------------ - - /// fair runtime database - FairRuntimeDb* rtdb = run->GetRuntimeDb(); - rtdb->setFirstInput(parIo1); - rtdb->setOutput(parIo1); - rtdb->saveOutput(); - - /// intialize and run - run->Init(); - run->Run(0, nEvents); - - timer.Stop(); - std::cout << "Macro finished succesfully." << std::endl; - std::cout << " Output file is " << (outDir + setupName + "_analysis.root") - << std::endl; - // std::cout << "Parameter file is " << parFile << std::endl; - std::cout << "Real time " << timer.RealTime() << " s, CPU time " - << timer.CpuTime() << " s" << std::endl; - std::cout << " Test passed" << std::endl; - std::cout << " All ok " << std::endl; -} diff --git a/analysis/PWGDIL/dielectron/papaframework/scripts/run_common_analysis.C b/analysis/PWGDIL/dielectron/papaframework/scripts/run_common_analysis.C new file mode 100644 index 0000000000000000000000000000000000000000..32537b614cc6fafe7e6d639e77366435be491657 --- /dev/null +++ b/analysis/PWGDIL/dielectron/papaframework/scripts/run_common_analysis.C @@ -0,0 +1,225 @@ +R__ADD_INCLUDE_PATH($PWD) + +#include "Config_dilepton_testing.C" +TString cfgFunc = "Config_dilepton_"; +TString configName = "testing"; + +// Includes needed for IDE +#if !defined(__CLING__) +#include "CbmDefs.h" +#include "CbmMCDataManager.h" +#include "CbmSetup.h" + +#include "FairSystemInfo.h" +#include <FairFileSource.h> +#include <FairMonitor.h> +#include <FairParAsciiFileIo.h> +#include <FairParRootFileIo.h> +#include <FairRunAna.h> +#include <FairRuntimeDb.h> + +#include <TStopwatch.h> +#endif + + +void run_common_analysis(Int_t nEvents = 0, Bool_t test = true, TString setup = "sis100_electron") +{ + + // ----- Environment -------------------------------------------------- + TString myName = "run_analysis"; // this macro's name for screen output + TString srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory + // ------------------------------------------------------------------------ + + // ----- In- and output directory names ------------------------------- + TString inDir = gSystem->Getenv("INDIR"); + TString inFile = gSystem->Getenv("INFILE"); // input file list + TString outDir = gSystem->Getenv("OUTDIR"); + // if(outDir.IsNull()) outFile = "."; + + // --- Logger settings ---------------------------------------------------- + gErrorIgnoreLevel = kFatal; //kInfo, kWarning, kError, kFatal; + // ------------------------------------------------------------------------ + + // ----- Load the geometry setup ------------------------------------- + /// load local copy of setup file (same as used for simulation and reconstruction) + std::cout << std::endl; + std::cout << "-I- " << myName << ": Loading setup " << setup << std::endl; + CbmSetup::Instance()->LoadSetup(setup); + + CbmSetup::Instance()->Print(); + + + // ----- run manager + I/O -------------------------------------------- + std::cout << std::endl; + FairRunAna* run = new FairRunAna(); + run->SetOutputFile(outDir + "/analysis.root"); + FairFileSource* src = NULL; + + /// stopwatch + TStopwatch timer; + timer.Start(); + + Int_t i = 0; + TString file; + // char filename[300]; + ifstream in(inFile); + TList* parList = new TList(); + TString parFile = ""; + + CbmMCDataManager* mcManager = new CbmMCDataManager("MCDataManager", 1); + + TString traFile = ""; + /// loop over all file in list + if (test) { + ifstream in(inFile); + while (in >> file) { + Int_t n = 1; + TString num = ""; + while (true) { + TString temp = (TString) file(file.Length() - n, file.Length() - 0); + // std::cout<<temp<<std::endl; + if (temp.Contains("/")) break; + num = temp; + n++; + } + + // // mc sim file + src = new FairFileSource(file + "/" + file(file.Length() - n + 1, file.Length() - 0) + ".tra.root"); + src->AddFriend(file + "/" + file(file.Length() - n + 1, file.Length() - 0) + ".event.raw.root"); + src->AddFriend(file + "/" + file(file.Length() - n + 1, file.Length() - 0) + ".rec.root"); + + parFile = file + "/" + file(file.Length() - n + 1, file.Length() - 0) + ".par.root"; + // parFile = "/lustre/cbm/users/ogolosov/mc/cbmsim/commonParFiles/"+file(60,file.Length()-n+1)+"/" + file(file.Length()-n+1,file.Length()-0)+".par.root"; + + traFile = file + "/" + file(file.Length() - n + 1, file.Length() - 0) + ".tra.root"; + + mcManager->AddFile(traFile); + + i++; + } + } + if (!test) { + ifstream in(outDir + "/" + inFile); + while (in >> file) { + Int_t n = 1; + TString num = ""; + while (true) { + TString temp = (TString) file(file.Length() - n, file.Length() - 0); + // std::cout<<temp<<std::endl; + if (temp.Contains("/")) break; + num = temp; + n++; + } + + // mc sim file + if (!i) src = new FairFileSource(file + "/" + file(file.Length() - n + 1, file.Length() - 0) + ".tra.root"); + else + src->AddFile(file + "/" + file(file.Length() - n + 1, file.Length() - 0) + ".tra.root"); + src->AddFriend(file + "/" + file(file.Length() - n + 1, file.Length() - 0) + ".event.raw.root"); + src->AddFriend(file + "/" + file(file.Length() - n + 1, file.Length() - 0) + ".rec.root"); + + parFile = file + "/" + file(file.Length() - n + 1, file.Length() - 0) + ".par.root"; + //parFile = "/lustre/cbm/users/ogolosov/mc/cbmsim/commonParFiles/"+file(60,file.Length()-n+1)+"/" + file(file.Length()-n+1,file.Length()-0)+".par.root"; + + traFile = file + "/" + file(file.Length() - n + 1, file.Length() - 0) + ".tra.root"; + + mcManager->AddFile(traFile); + + i++; + } + } + + // parList->Dump(); + // add source to run + run->SetSource(src); + + // // ------------------------------------------------------------------------ + + // ----- MCDataManager (legacy mode) ----------------------------------- + run->AddTask(mcManager); + // ------------------------------------------------------------------------ + + CbmMatchRecoToMC* match1 = new CbmMatchRecoToMC(); + run->AddTask(match1); + + + // ----- L1/KF tracking + PID ----------------------------------------- + std::cout << std::endl; + std::cout << "-I- " << myName << ": Loading tasks " << std::endl; + + //CbmKF is needed for Extrapolation + CbmKF* kf = new CbmKF(); + run->AddTask(kf); + std::cout << "-I- : Added task " << kf->GetName() << std::endl; + + CbmL1* l1 = new CbmL1(); + // --- Material budget file names + if (CbmSetup::Instance()->IsActive(ECbmModuleId::kMvd)) { + TString geoTag; + CbmSetup::Instance()->GetGeoTag(ECbmModuleId::kMvd, geoTag); + TString matFile = gSystem->Getenv("VMCWORKDIR"); + matFile = matFile + "/parameters/mvd/mvd_matbudget_" + geoTag + ".root"; + std::cout << "Using material budget file " << matFile << std::endl; + l1->SetMvdMaterialBudgetFileName(matFile.Data()); + } + if (CbmSetup::Instance()->IsActive(ECbmModuleId::kSts)) { + TString geoTag; + CbmSetup::Instance()->GetGeoTag(ECbmModuleId::kSts, geoTag); + TString matFile = gSystem->Getenv("VMCWORKDIR"); + matFile = matFile + "/parameters/sts/sts_matbudget_" + geoTag + ".root"; + // matFile = matFile + "/parameters/sts/sts_matbudget_v19a.root"; + // matFile = matFile + "/parameters/sts/sts_matbudget_v19i.root"; + std::cout << "Using material budget file " << matFile << std::endl; + l1->SetStsMaterialBudgetFileName(matFile.Data()); + // } + run->AddTask(l1); + std::cout << "-I- : Added task " << l1->GetName() << std::endl; + } + + // --- TRD pid tasks + if (CbmSetup::Instance()->IsActive(ECbmModuleId::kTrd)) { + CbmTrdSetTracksPidLike* trdLI = new CbmTrdSetTracksPidLike("TRDLikelihood", "TRDLikelihood"); + trdLI->SetUseMCInfo(kTRUE); + trdLI->SetUseMomDependence(kTRUE); + run->AddTask(trdLI); + std::cout << "-I- : Added task " << trdLI->GetName() << std::endl; + // ------------------------------------------------------------------------ + } + + // ----- PAPA tasks --------------------------------------------------- + std::cout << std::endl; + std::cout << "-I- " << myName << ": Loading private tasks " << std::endl; + + TString cfgPath = outDir + "/../"; + TString cfgFile = cfgPath + cfgFunc + configName + ".C"; //cfgPath + cfgFunc + configName + ".C"; + TString cfgFunct = cfgFunc + configName + "(\"" + configName + "\")"; + TString cfg = cfgFunc + configName + "()"; + gROOT->LoadMacro((cfgPath + cfgFunct + configName + ".C")); + FairTask* task = reinterpret_cast<FairTask*>(gROOT->ProcessLine(cfg)); + run->AddTask(task); + + // // ------------------------------------------------------------------------ + + // set parameter list + FairParRootFileIo* parIo1 = new FairParRootFileIo(); + // parIo1->open(parFile.Data(),"UPDATE"); + parIo1->open(parFile.Data(), "READ"); + // parIo1->open(parList); + FairRuntimeDb* rtdb = run->GetRuntimeDb(); + rtdb->setFirstInput(parIo1); + rtdb->setOutput(parIo1); + // rtdb->clearRunList(); + rtdb->saveOutput(); + + /// intialize and run + run->Init(); + run->Run(0, nEvents); + + timer.Stop(); + std::cout << "Macro finished succesfully." << std::endl; + std::cout << " Output file is " << (outDir + "analysis.root") << std::endl; + // std::cout << "Parameter file is " << parFile << std::endl; + std::cout << "Real time " << timer.RealTime() << " s, CPU time " << timer.CpuTime() << " s" << std::endl; + std::cout << " Test passed" << std::endl; + std::cout << " All ok " << std::endl; +} diff --git a/analysis/PWGDIL/dielectron/papaframework/scripts/run_reco_new.C b/analysis/PWGDIL/dielectron/papaframework/scripts/run_reco_new.C deleted file mode 100644 index 26c3d3bac53cd30cc25debdbc85a94a1569d89f5..0000000000000000000000000000000000000000 --- a/analysis/PWGDIL/dielectron/papaframework/scripts/run_reco_new.C +++ /dev/null @@ -1,310 +0,0 @@ -// -------------------------------------------------------------------------- -// -// Macro for reconstruction of simulated events with standard settings -// -// HitProducers in MVD, RICH, TRD, TOF, ECAL -// Digitizer and HitFinder in STS -// FAST MC for ECAL -// STS track finding and fitting (L1 / KF) -// TRD track finding and fitting (L1 / KF) -// RICH ring finding (ideal) and fitting -// Global track finding (ideal), rich assignment -// Primary vertex finding (ideal) -// Matching of reconstructed and MC tracks in STS, RICH and TRD -// -// V. Friese 24/02/2006 -// Version 24/04/2007 (V. Friese) -// -// -------------------------------------------------------------------------- - - -void run_reco_new(Int_t nEvents = 2, - const char* setupName = "sis100_electron", - const char* inputFile = "") { - - // ======================================================================== - // Adjust this part according to your requirements - - // --- Logger settings ---------------------------------------------------- - TString logLevel = "INFO"; - TString logVerbosity = "LOW"; - // ------------------------------------------------------------------------ - - - // ----- Environment -------------------------------------------------- - TString myName = "run_reco_new"; // this macro's name for screen output - TString srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory - // ------------------------------------------------------------------------ - - - // ----- In- and output file names ------------------------------------ - TString outDir = gSystem->Getenv("OUTDIR"); - if (outDir.IsNull()) outDir = "data/"; - TString inFile = outDir + setupName + "_mc.root"; // Input file (MC events) - TString outFile = outDir + setupName + "_reco.root"; // reco file - TString parFile = outDir + setupName + "_params.root"; // Parameter file - // ------------------------------------------------------------------------ - - - // ----- Remove old CTest runtime dependency file ---------------------- - TString depFile = Remove_CTest_Dependency_File(outDir, "run_reco", setupName); - // ------------------------------------------------------------------------ - - - // ----- Load the geometry setup ------------------------------------- - std::cout << std::endl; - TString setupFile = outDir + "/../setup_" + setupName + ".C"; - // TString setupFile = srcDir + "/geometry/setup/setup_" + setupName + ".C"; - TString setupFunct = "setup_"; - setupFunct = setupFunct + setupName + "()"; - std::cout << "-I- " << myName << ": Loading macro " << setupFile << std::endl; - gROOT->LoadMacro(setupFile); - gROOT->ProcessLine(setupFunct); - CbmSetup* setup = CbmSetup::Instance(); - // ------------------------------------------------------------------------ - - - // ----- Modify Cbm setup --------------------------------------------- - setup->SetActive(kPsd, kFALSE); // dont write psd points - // --- remove detector geometries - setup->RemoveModule(kPsd); // remove psd from setup - // setup->RemoveModule(kMvd); // remove mvd and its material budget - // setup->RemoveModule(kTrd); // e.g. for sts-tof-matching study - // --- change default geomerties - // setup->SetModule(kTrd, "v15d_1e", kTRUE); // 5 TRD layer - std::cout << "-I- " << myName << ": CbmSetup updated " << std::endl; - // ------------------------------------------------------------------------ - - - // ----- Parameter files as input to the runtime database ------------- - std::cout << std::endl; - std::cout << "-I- " << myName << ": Defining parameter files " << std::endl; - TList* parFileList = new TList(); - TString geoTag; - - // - TRD digitisation parameters - if (setup->GetGeoTag(kTrd, geoTag)) { - TObjString* trdFile = - new TObjString(srcDir + "/parameters/trd/trd_" + geoTag + ".digi.par"); - parFileList->Add(trdFile); - std::cout << "-I- " << myName << ": Using parameter file " - << trdFile->GetString() << std::endl; - } - - // - TOF digitisation parameters - if (setup->GetGeoTag(kTof, geoTag)) { - TObjString* tofFile = - new TObjString(srcDir + "/parameters/tof/tof_" + geoTag + ".digi.par"); - parFileList->Add(tofFile); - std::cout << "-I- " << myName << ": Using parameter file " - << tofFile->GetString() << std::endl; - TObjString* tofBdfFile = - new TObjString(srcDir + "/parameters/tof/tof_" + geoTag + ".digibdf.par"); - parFileList->Add(tofBdfFile); - std::cout << "-I- " << myName << ": Using parameter file " - << tofBdfFile->GetString() << std::endl; - } - // ------------------------------------------------------------------------ - - // In general, the following parts need not be touched - // ======================================================================== - - - // ----- Timer -------------------------------------------------------- - TStopwatch timer; - timer.Start(); - // ------------------------------------------------------------------------ - - - // ---- Debug option ------------------------------------------------- - gDebug = 0; - // ------------------------------------------------------------------------ - - - // ----- Input file --------------------------------------------------- - std::cout << std::endl; - TString defaultInputFile = "data/"; - defaultInputFile = defaultInputFile + setupName + "_test.mc.root"; - if (inFile.IsNull()) { // Not defined in the macro explicitly - if (strcmp(inputFile, "") == 0) { // not given as argument to the macro - inFile = defaultInputFile; - } else - inFile = inputFile; - } - std::cout << "-I- " << myName << ": Using input file " << inFile << std::endl; - // ------------------------------------------------------------------------ - - - // ----- FairRunAna --------------------------------------------------- - FairRunAna* run = new FairRunAna(); - run->SetInputFile(inFile); - run->SetOutputFile(outFile); - run->SetGenerateRunInfo(kTRUE); - Bool_t hasFairMonitor = Has_Fair_Monitor(); - if (hasFairMonitor) FairMonitor::GetMonitor()->EnableMonitor(kTRUE); - // ------------------------------------------------------------------------ - - - // ----- Logger settings ---------------------------------------------- - FairLogger* gLogger = FairLogger::GetLogger(); - gLogger->SetLogScreenLevel(logLevel.Data()); - gLogger->SetLogVerbosityLevel(logVerbosity.Data()); - // ------------------------------------------------------------------------ - - - // ----- Digitisers --------------------------------------------------- - std::cout << std::endl; - TString macroName = gSystem->Getenv("VMCWORKDIR"); - macroName += "/macro/run/modules/digitize.C"; - std::cout << "Loading macro " << macroName << std::endl; - gROOT->LoadMacro(macroName); - gROOT->ProcessLine("digitize()"); - // ------------------------------------------------------------------------ - - - // ----- Reconstruction tasks ----------------------------------------- - std::cout << std::endl; - macroName = srcDir + "/macro/run/modules/reconstruct.C"; - std::cout << "Loading macro " << macroName << std::endl; - gROOT->LoadMacro(macroName); - gROOT->ProcessLine("reconstruct()"); - // ------------------------------------------------------------------------ - - - // ----- Matching to Monte-Carlo tasks -------------------------------- - std::cout << std::endl; - CbmMCDataManager* mcManager = new CbmMCDataManager("MCManager", 1); - mcManager->AddFile(inFile); - run->AddTask(mcManager); - std::cout << "-I- : Added task " << mcManager->GetName() << std::endl; - CbmMatchRecoToMC* matchTask = new CbmMatchRecoToMC(); - run->AddTask(matchTask); - std::cout << "-I- : Added task " << matchTask->GetName() << std::endl; - // ------------------------------------------------------------------------ - - - // ----- PID tasks ---------------------------------------------------- - if (setup->IsActive(kTrd)) { - - // ----------- TRD track Pid Wkn ---------------------- - CbmTrdSetTracksPidWkn* trdSetTracksPidTask = - new CbmTrdSetTracksPidWkn("trdFindTracks", "trdFindTracks"); - run->AddTask(trdSetTracksPidTask); - std::cout << "-I- : Added task " << trdSetTracksPidTask->GetName() - << std::endl; - // ---------------------------------------------------- - - // ----------- TRD track Pid Ann ---------------------- - CbmTrdSetTracksPidANN* trdSetTracksPidAnnTask = - new CbmTrdSetTracksPidANN("Ann", "Ann"); - run->AddTask(trdSetTracksPidAnnTask); - std::cout << "-I- : Added task " << trdSetTracksPidAnnTask->GetName() - << std::endl; - // ---------------------------------------------------- - - // ----------- TRD track Pid Like ---------------------- - // Since in the newest version of this method depends on the global - // track the task has to move after the global tracking - CbmTrdSetTracksPidLike* trdSetTracksPidLikeTask = - new CbmTrdSetTracksPidLike("Likelihood", "Likelihood"); - run->AddTask(trdSetTracksPidLikeTask); - std::cout << "-I- : Added task " << trdSetTracksPidLikeTask->GetName() - << std::endl; - // ---------------------------------------------------- - } - // ------------------------------------------------------------------------ - - - // ----- control data persistency ------------------------------------- - std::cout << std::endl; - if (setup->IsActive(kMvd)) { - run->GetTask("MVD Digitiser")->SetOutputBranchPersistent("MvdDigi", kFALSE); - run->GetTask("MVD Digitiser") - ->SetOutputBranchPersistent("MvdDigiMatch", kFALSE); - } - if (setup->IsActive(kSts)) { - run->GetTask("StsDigitize")->SetOutputBranchPersistent("StsDigi", kFALSE); - run->GetTask("StsDigitize") - ->SetOutputBranchPersistent("StsDigiMatch", kFALSE); - } - if (setup->IsActive(kRich)) { - run->GetTask("CbmRichDigitizer") - ->SetOutputBranchPersistent("RichDigi", kFALSE); - } - if (setup->IsActive(kMuch)) { - run->GetTask("MuchDigitizeGem") - ->SetOutputBranchPersistent("MuchDigi", kFALSE); //gem - run->GetTask("MuchDigitizeGem") - ->SetOutputBranchPersistent("MuchDigiMatch", kFALSE); - run->GetTask("CbmMuchDigitizeStraws") - ->SetOutputBranchPersistent("MuchStrawDigi", kFALSE); //straw - run->GetTask("CbmMuchDigitizeStraws") - ->SetOutputBranchPersistent("MuchStrawDigiMatch", kFALSE); - } - if (setup->IsActive(kTrd)) { - run->GetTask("CbmTrdDigitizerPRF") - ->SetOutputBranchPersistent("TrdDigi", kFALSE); - run->GetTask("CbmTrdDigitizerPRF") - ->SetOutputBranchPersistent("TrdDigiMatch", kFALSE); - } - if (setup->IsActive(kTof)) { - run->GetTask("TOF Digitizer BDF") - ->SetOutputBranchPersistent("TofDigi", kFALSE); - run->GetTask("TOF Digitizer BDF") - ->SetOutputBranchPersistent("TofDigiMatchPoints", kFALSE); - run->GetTask("TOF Simple Clusterizer") - ->SetOutputBranchPersistent("TofDigiMatch", kFALSE); - matchTask->SetOutputBranchPersistent("TofDigiMatch", kFALSE); - } - std::cout << "-I- : Data persistency controlled " << std::endl; - // ------------------------------------------------------------------------ - - - // ----- Parameter database -------------------------------------------- - std::cout << std::endl << std::endl; - std::cout << "-I- " << myName << ": Set runtime DB" << std::endl; - FairRuntimeDb* rtdb = run->GetRuntimeDb(); - FairParRootFileIo* parIo1 = new FairParRootFileIo(); - FairParAsciiFileIo* parIo2 = new FairParAsciiFileIo(); - parIo1->open(parFile.Data()); - parIo2->open(parFileList, "in"); - rtdb->setFirstInput(parIo1); - rtdb->setSecondInput(parIo2); - rtdb->setOutput(parIo1); - rtdb->saveOutput(); - rtdb->print(); - // ------------------------------------------------------------------------ - - - // ----- Run initialisation ------------------------------------------- - std::cout << std::endl; - std::cout << "-I- " << myName << ": Initialise run" << std::endl; - run->Init(); - // ------------------------------------------------------------------------ - - - // ----- Start run ---------------------------------------------------- - std::cout << std::endl << std::endl; - std::cout << "-I- " << myName << ": Starting run" << std::endl; - run->Run(0, nEvents); - // ------------------------------------------------------------------------ - - - // ----- Finish ------------------------------------------------------- - timer.Stop(); - Double_t rtime = timer.RealTime(); - Double_t ctime = timer.CpuTime(); - std::cout << std::endl << std::endl; - std::cout << "Macro finished succesfully." << std::endl; - std::cout << "Output file is " << outFile << std::endl; - std::cout << "Parameter file is " << parFile << std::endl; - std::cout << "Real time " << rtime << " s, CPU time " << ctime << " s" - << std::endl; - std::cout << std::endl; - std::cout << " Test passed" << std::endl; - std::cout << " All ok " << std::endl; - // ------------------------------------------------------------------------ - - // Function needed for CTest runtime dependency - Generate_CTest_Dependency_File(depFile); -} diff --git a/analysis/PWGDIL/dielectron/papaframework/scripts/run_sim_new.C b/analysis/PWGDIL/dielectron/papaframework/scripts/run_sim_new.C deleted file mode 100644 index 4468cee64255a3a8934c5db7c29d8793d0ec131f..0000000000000000000000000000000000000000 --- a/analysis/PWGDIL/dielectron/papaframework/scripts/run_sim_new.C +++ /dev/null @@ -1,815 +0,0 @@ -// -------------------------------------------------------------------------- -// -// Macro for standard transport simulation using UrQMD input and GEANT3 -// -// V. Friese 22/02/2007 -// -// Version 2016-02-05 -// -// For the setup (geometry and field), predefined setups can be chosen -// by the second argument. A list of available setups is given below. -// The input file can be defined explicitly in this macro or by the -// third argument. If none of these options are chosen, a default -// input file distributed with the source code is selected. -// -// 2014-06-30 - DE - available setups from geometry/setup: -// 2014-06-30 - DE - sis100_hadron -// 2014-06-30 - DE - sis100_electron -// 2014-06-30 - DE - sis100_muon -// 2014-06-30 - DE - sis300_electron -// 2014-06-30 - DE - sis300_muon -// -// -------------------------------------------------------------------------- - - -void AddUrqmdGenerator(FairPrimaryGenerator* primGen, - Int_t idx, - TString inputFile); -void AddLMVMCocktail(FairPrimaryGenerator* primGen, - Int_t idx, - TString inputFile); -void AddRadiationCocktail(FairPrimaryGenerator* primGen, - Int_t idx, - TString inputFile); -void AddCharmoniaCocktail(FairPrimaryGenerator* primGen, - Int_t idx, - TString inputFile); -void AddFragmentsCocktail(FairPrimaryGenerator* primGen, Int_t idx); -void AddBoxGenerator(FairPrimaryGenerator* primGen, Int_t mult); - -void SetParticleDecays(); -void ConfigureMCStack(); - - -void run_sim_new(Int_t nEvents = 2, - const char* setupName = "sis100_electron", - const char* inputFile = "") { - - // ======================================================================== - // Adjust this part according to your requirements - - // --- Logger settings ---------------------------------------------------- - TString logLevel = "INFO"; - TString logVerbosity = "LOW"; - // ------------------------------------------------------------------------ - - - // ----- Environment -------------------------------------------------- - TString myName = "run_sim_new"; // this macro's name for screen output - TString srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory - // ------------------------------------------------------------------------ - - - // ----- In- and output file names ------------------------------------ - TString outDir = gSystem->Getenv("OUTDIR"); - if (outDir.IsNull()) outDir = "data/"; - TString geoFile = outDir + setupName + "_geofile_full.root"; - TString outFile = outDir + setupName + "_mc.root"; - TString parFile = outDir + setupName + "_params.root"; - // ------------------------------------------------------------------------ - - - // ----- Remove old CTest runtime dependency file --------------------- - TString depFile = Remove_CTest_Dependency_File(outDir, "run_sim", setupName); - // ------------------------------------------------------------------------ - - - // ----- Load the geometry setup ------------------------------------- - std::cout << std::endl; - TString setupFile = outDir + "/../setup_" + setupName + ".C"; - // TString setupFile = srcDir + "/geometry/setup/setup_" + setupName + ".C"; - TString setupFunct = "setup_"; - setupFunct = setupFunct + setupName + "()"; - std::cout << "-I- " << myName << ": Loading macro " << setupFile << std::endl; - gROOT->LoadMacro(setupFile); - gROOT->ProcessLine(setupFunct); - CbmSetup* setup = CbmSetup::Instance(); - // ------------------------------------------------------------------------ - - - // ----- Modify Cbm setup --------------------------------------------- - std::cout << std::endl; - setup->SetActive(kPsd, kFALSE); // dont write psd points - // --- remove detector geometries - setup->RemoveModule(kPsd); // remove psd from setup - // setup->RemoveModule(kMvd); // remove mvd and its material budget - // setup->RemoveModule(kTrd); // e.g. for sts-tof-matching study - // --- change default geomerties - // setup->SetModule(kTrd, "v15d_1e", kTRUE); // 5 TRD layer - std::cout << "-I- " << myName << ": CbmSetup updated " << std::endl; - // ------------------------------------------------------------------------ - - - // --- Define the target geometry ----------------------------------------- - // - // The target is not part of the setup, since one and the same setup can - // and will be used with different targets. - // The target is constructed as a tube in z direction with the specified - // diameter (in x and y) and thickness (in z). It will be placed at the - // specified position as daughter volume of the volume present there. It is - // in the responsibility of the user that no overlaps or extrusions are - // created by the placement of the target. - // - TString targetElement = "Gold"; - Double_t targetThickness = 0.0025; // full thickness in cm - Double_t targetDiameter = 2.5; // diameter in cm - Double_t targetPosX = 0.; // target x position in global c.s. [cm] - Double_t targetPosY = 0.; // target y position in global c.s. [cm] - Double_t targetPosZ = 0.; // target z position in global c.s. [cm] - Double_t targetRotY = 0.; // target rotation angle around the y axis [deg] - // ------------------------------------------------------------------------ - - - // --- Define the creation of the primary vertex ------------------------ - // - // By default, the primary vertex point is sampled from a Gaussian - // distribution in both x and y with the specified beam profile width, - // and from a flat distribution in z over the extension of the target. - // By setting the respective flags to kFALSE, the primary vertex will always - // at the (0., 0.) in x and y and in the z centre of the target, respectively. - // - Bool_t smearVertexXY = kTRUE; - Bool_t smearVertexZ = kTRUE; - Double_t beamWidthX = 1.; // Gaussian sigma of the beam profile in x [cm] - Double_t beamWidthY = 1.; // Gaussian sigma of the beam profile in y [cm] - // ------------------------------------------------------------------------ - - - // In general, the following parts need not be touched - // ======================================================================== - - - // ----- Timer -------------------------------------------------------- - TStopwatch timer; - timer.Start(); - // ------------------------------------------------------------------------ - - - // ---- Debug option ------------------------------------------------- - gDebug = 0; - // ------------------------------------------------------------------------ - - - // ----- Input file --------------------------------------------------- - std::cout << std::endl; - TString inFile = ""; // give here or as argument; otherwise default is taken - TString defaultInputFile = srcDir + "/input/urqmd.auau.10gev.centr.root"; - if (inFile.IsNull()) { // Not defined in the macro explicitly - if (strcmp(inputFile, "") == 0) { // not given as argument to the macro - inFile = defaultInputFile; - } else - inFile = inputFile; - } - // std::cout << "-I- " << myName << ": Using input file " << inFile << std::endl; - // ------------------------------------------------------------------------ - - - // ----- Create simulation run ---------------------------------------- - FairRunSim* run = new FairRunSim(); - run->SetName("TGeant3"); // Transport engine - run->SetOutputFile(outFile); // Output file - run->SetGenerateRunInfo(kTRUE); // Create FairRunInfo file - // ------------------------------------------------------------------------ - - - // ----- Logger settings ---------------------------------------------- - FairLogger* gLogger = FairLogger::GetLogger(); - gLogger->SetLogScreenLevel(logLevel.Data()); - gLogger->SetLogVerbosityLevel(logVerbosity.Data()); - // ------------------------------------------------------------------------ - - - // ----- Create media ------------------------------------------------- - std::cout << std::endl; - std::cout << "-I- " << myName << ": Setting media file" << std::endl; - run->SetMaterials("media.geo"); // Materials - // ------------------------------------------------------------------------ - - - // ----- Create and register modules ---------------------------------- - std::cout << std::endl; - TString macroName = gSystem->Getenv("VMCWORKDIR"); - macroName += "/macro/run/modules/registerSetup.C"; - std::cout << "Loading macro " << macroName << std::endl; - gROOT->LoadMacro(macroName); - gROOT->ProcessLine("registerSetup()"); - // ------------------------------------------------------------------------ - - - // ----- Create and register the target ------------------------------- - std::cout << std::endl; - std::cout << "-I- " << myName << ": Registering target" << std::endl; - CbmTarget* target = - new CbmTarget(targetElement.Data(), targetThickness, targetDiameter); - target->SetPosition(targetPosX, targetPosY, targetPosZ); - target->SetRotation(targetRotY); - target->Print(); - run->AddModule(target); - // ------------------------------------------------------------------------ - - // ----- Create magnetic field ---------------------------------------- - std::cout << std::endl; - std::cout << "-I- " << myName << ": Registering magnetic field" << std::endl; - CbmFieldMap* magField = CbmSetup::Instance()->CreateFieldMap(); - if (!magField) { - std::cout << "-E- " << myName << ": No valid field!"; - return; - } - run->SetField(magField); - // ------------------------------------------------------------------------ - - // ----- Create PrimaryGenerator -------------------------------------- - std::cout << std::endl; - std::cout << "-I- " << myName << ": Registering event generators" - << std::endl; - FairPrimaryGenerator* primGen = new FairPrimaryGenerator(); - // --- Uniform distribution of event plane angle - primGen->SetEventPlane(0., 2. * TMath::Pi()); - // --- Get target parameters - Double_t tX = 0.; - Double_t tY = 0.; - Double_t tZ = 0.; - Double_t tDz = 0.; - if (target) { - target->GetPosition(tX, tY, tZ); - tDz = target->GetThickness(); - } - primGen->SetTarget(tZ, tDz); - primGen->SetBeam(0., 0., beamWidthX, beamWidthY); - primGen->SmearGausVertexXY(smearVertexXY); - primGen->SmearVertexZ(smearVertexZ); - // - // TODO: Currently, there is no guaranteed consistency of the beam profile - // and the transversal target dimension, i.e., that the sampled primary - // vertex falls into the target volume. This would require changes - // in the FairPrimaryGenerator class. - // ------------------------------------------------------------------------ - - // --- Create particle cocktail - TObjArray* arr = outDir.Tokenize("/"); - Int_t idx = 1; - if (arr->GetEntriesFast() > 3) { - TString tmp = (static_cast<TObjString*>(arr->Last()))->GetString(); - idx = tmp.Atoi(); - } - delete arr; - - // NOTE: for random file selection set idx to -1 - idx = -1; - gRandom->SetSeed(0); - - TString location = "/lustre/nyx"; //"/hera"; //"/scratch" - - TString charmoniaFile = location + "/cbm/users/jbook/pluto/train/"; - /// NOTE: charmonia generator needs to be added before the urqmd, because of evtgen - //AddCharmoniaCocktail( primGen, idx, charmoniaFile ); - - /// AA urqmd - TString urqmdFile = - // location + "/cbm/prod/gen/urqmd/auau/8gev/mbias/urqmd.auau.8gev.mbias."; - // location + "/cbm/prod/gen/urqmd/auau/8gev/centr/urqmd.auau.8gev.centr."; - location + "/cbm/users/jbook/urqmd/auau/8gev/centr010/"; //selfmade UrQMD - /// pA urqmd - // location + "/cbm/users/jbook/urqmd/pau/30gev/centr/"; //selfmade UrQMD - // AddUrqmdGenerator( primGen, idx, urqmdFile ); - - TString plutoFile = location + "/cbm/users/jbook/pluto/train/"; - // location + "/cbm/users/ekrebs/pluto/jun15/auau/cktA/8gev/"; //for mumu - // AddLMVMCocktail( primGen, idx, plutoFile ); - - TString radiationFile = - location + "/cbm/users/jbook/pluto/inmed/out_rapp_pluto_ee_inmed_"; - // location + "/cbm/users/ekrebs/pluto/jun15/auau/cktRapp/8gev/"; // for mumu - //AddRadiationCocktail( primGen, idx, radiationFile ); - - // AddFragmentsCocktail( primGen, 1 ); - - AddBoxGenerator(primGen, 2); - - - run->SetGenerator(primGen); - // ------------------------------------------------------------------------ - - - // ----- Run initialisation ------------------------------------------- - std::cout << std::endl; - std::cout << "-I- " << myName << ": Initialise run" << std::endl; - run->Init(); // this inits gMC and the stack - // ------------------------------------------------------------------------ - - // ----- Particle decays in geant3 ------------------------------------ - SetParticleDecays(); - // ------------------------------------------------------------------------ - - // ----- Configure MC stack ------------------------------------------- - ConfigureMCStack(); - // ------------------------------------------------------------------------ - - - // ----- Parameter database --------------------------------------------- - std::cout << std::endl << std::endl; - std::cout << "-I- " << myName << ": Set runtime DB" << std::endl; - FairRuntimeDb* rtdb = run->GetRuntimeDb(); - CbmFieldPar* fieldPar = (CbmFieldPar*) rtdb->getContainer("CbmFieldPar"); - fieldPar->SetParameters(magField); - fieldPar->setChanged(); - fieldPar->setInputVersion(run->GetRunId(), 1); - Bool_t kParameterMerged = kTRUE; - FairParRootFileIo* rootIo = new FairParRootFileIo(kParameterMerged); - rootIo->open(parFile.Data()); - rtdb->setOutput(rootIo); - rtdb->saveOutput(); - rtdb->print(); - // ------------------------------------------------------------------------ - - - // ----- Start run ---------------------------------------------------- - std::cout << std::endl << std::endl; - std::cout << "-I- " << myName << ": Starting run" << std::endl; - run->Run(nEvents); - // ------------------------------------------------------------------------ - - - // ----- Finish ------------------------------------------------------- - run->CreateGeometryFile(geoFile); - timer.Stop(); - Double_t rtime = timer.RealTime(); - Double_t ctime = timer.CpuTime(); - std::cout << std::endl << std::endl; - std::cout << "Macro finished successfully." << std::endl; - std::cout << "Output file is " << outFile << std::endl; - std::cout << "Parameter file is " << parFile << std::endl; - std::cout << "Geometry file is " << geoFile << std::endl; - std::cout << "Real time " << rtime << " s, CPU time " << ctime << "s" - << std::endl - << std::endl; - // ------------------------------------------------------------------------ - - - // ----- Resource monitoring ------------------------------------------ - if (Has_Fair_Monitor()) { // FairRoot Version >= 15.11 - // Extract the maximal used memory an add is as Dart measurement - // This line is filtered by CTest and the value send to CDash - FairSystemInfo sysInfo; - Float_t maxMemory = sysInfo.GetMaxMemory(); - std::cout << "<DartMeasurement name=\"MaxMemory\" type=\"numeric/double\">"; - std::cout << maxMemory; - std::cout << "</DartMeasurement>" << std::endl; - - Float_t cpuUsage = ctime / rtime; - std::cout << "<DartMeasurement name=\"CpuLoad\" type=\"numeric/double\">"; - std::cout << cpuUsage; - std::cout << "</DartMeasurement>" << std::endl; - } - - std::cout << " Test passed" << std::endl; - std::cout << " All ok " << std::endl; - - // Function needed for CTest runtime dependency - Generate_CTest_Dependency_File(depFile); - // ------------------------------------------------------------------------ -} - -// ------------------------------------------------------------------------ -void AddUrqmdGenerator(FairPrimaryGenerator* primGen, - Int_t idx, - TString inputFile) { - /// UrQMD version 3.3 - 5k files with 1k events each - /// checkout https://lxcbmredmine01.gsi.de/projects/cbmroot/wiki/Mass_Production_at_GSI - /// - Bool_t isAA = inputFile.Contains("auau"); - Bool_t isOWN = inputFile.Contains("jbook"); - - Int_t module = (isAA ? (isOWN ? 10000 : 5000) : 10000); - Int_t idy = (idx >= 0 ? idx % module : gRandom->Integer(module)); - - TString file = ""; - if (isOWN) - file = inputFile + Form("/%05d", idy + 1) + "/urqmd.root"; - else if (isAA) - file = inputFile + Form("%05d", idy + 1) + ".root"; /// AA - else - file = inputFile + Form("%05d", idy + 1) + ".urqmd.f14"; /// pA - printf("urqmd file: %s \n", file.Data()); - - if (isAA || isOWN) { - CbmUnigenGenerator* uniGen = new CbmUnigenGenerator(file); - primGen->AddGenerator(uniGen); - } else { - FairUrqmdGenerator* urqmdGen = new FairUrqmdGenerator(file); - urqmdGen->SetEventPlane(0., 360.); - primGen->AddGenerator(urqmdGen); - } - return; -} - -// ------------------------------------------------------------------------ -void AddLMVMCocktail(FairPrimaryGenerator* primGen, - Int_t idx, - TString inputFile) { - /// - /// DIELECTRON-COCKTAIL: add particles via pluto generator - /// - - Int_t module = 10000; - /* - TString idy=""; - if(inputFile.Contains("ekrebs")) { - // dimuon - idy=Form("%04d",idx+1); - inputFile += "/rho0/mpmm/pluto.auau.8gev.rho0.mpmm." +idy + ".root"; - } else { - // dielectron - idy=Form("%05d",(idx%10000)+1); - inputFile += idy + "/pluto.auau.8gev.rho0.epem.root"; - } - printf("pluto vec.meson: %s \n",inputFile.Data()); - */ - - Int_t idy = (idx >= 0 ? idx % module : gRandom->Integer(module)); - TString file = - inputFile + Form("%05d", idy + 1) + "/pluto.auau.8gev.rho0.epem.root"; - printf("pluto vec.meson: %s \n", file.Data()); - CbmPlutoGenerator* plutoRho = new CbmPlutoGenerator(file); - // primGen->AddGenerator(plutoRho); // not added, included in inmed spectral func - - idy = (idx >= 0 ? idx % module : gRandom->Integer(module)); - file = inputFile + Form("%05d", idy + 1) + "/pluto.auau.8gev.rho0.epem.root"; - file.ReplaceAll("rho0", "omega"); - printf("pluto vec.meson: %s \n", file.Data()); - CbmPlutoGenerator* plutoOmega = new CbmPlutoGenerator(file); - primGen->AddGenerator(plutoOmega); - - idy = (idx >= 0 ? idx % module : gRandom->Integer(module)); - file = inputFile + Form("%05d", idy + 1) + "/pluto.auau.8gev.rho0.epem.root"; - file.ReplaceAll("rho0", "phi"); - printf("pluto vec.meson: %s \n", file.Data()); - CbmPlutoGenerator* plutoPhi = new CbmPlutoGenerator(file); - primGen->AddGenerator(plutoPhi); - - idy = (idx >= 0 ? idx % module : gRandom->Integer(module)); - file = inputFile + Form("%05d", idy + 1) + "/pluto.auau.8gev.rho0.epem.root"; - file.ReplaceAll("rho0", "omega"); - file.ReplaceAll("epem", "pi0epem"); - file.ReplaceAll("mpmm", "pi0mpmm"); - printf("pluto vec.meson: %s \n", file.Data()); - CbmPlutoGenerator* plutoOmegaDalitz = new CbmPlutoGenerator(file); - primGen->AddGenerator(plutoOmegaDalitz); - - // idy = (idx >= 0 ? idx%module : gRandom->Integer(module) ); - // file = inputFile + Form("%05d",idy+1) + "/pluto.auau.8gev.rho0.epem.root"; - // file.ReplaceAll("rho0","phi"); - // file.ReplaceAll("epem","kpkm"); - // file.ReplaceAll("mpmm","kpkm"); - // printf("pluto vec.meson: %s \n",file.Data()); - // CbmPlutoGenerator *plutoPhiKK= new CbmPlutoGenerator(file); - // primGen->AddGenerator(plutoPhiKK); - - return; -} - -// ------------------------------------------------------------------------ -void AddRadiationCocktail(FairPrimaryGenerator* primGen, - Int_t idx, - TString inputFile) { - /// - /// in-medium and qgp radiation - /// - - Int_t module = 500; - /* - TString idy=""; - if(inputFile.Contains("ekrebs")) { - // dimuon - idy=Form("%04d",(idx%5000)+1); - inputFile += "rapp.inmed/mpmm/pluto.auau.25gev.rapp.inmed.mpmm." + idy + ".root"; - } else { - // dielectron - idy=Form("%04d",(idx%500)+1); - inputFile += idy + ".root"; - } - */ - printf("pluto radiation: %s \n", inputFile.Data()); - - Int_t idy = (idx >= 0 ? idx % module : gRandom->Integer(module)); - TString file = inputFile + Form("%04d", idy + 1) + ".root"; - printf("pluto radiation: %s \n", file.Data()); - CbmPlutoGenerator* plutoRadInMed = new CbmPlutoGenerator(file); - plutoRadInMed->SetManualPDG(99009011); - primGen->AddGenerator(plutoRadInMed); - - if (!inputFile.Contains("ekrebs")) { - idy = (idx >= 0 ? idx % module : gRandom->Integer(module)); - file = inputFile + Form("%04d", idy + 1) + ".root"; - file.ReplaceAll("inmed", "qgp"); - printf("pluto radiation: %s \n", file.Data()); - CbmPlutoGenerator* plutoRadQGP = new CbmPlutoGenerator(file); - plutoRadQGP->SetManualPDG(99009111); - primGen->AddGenerator(plutoRadQGP); - } - return; -} - -// ------------------------------------------------------------------------ -void AddFragmentsCocktail(FairPrimaryGenerator* primGen, Int_t mult) { - /// - /// HeavyFragments-COCKTAIL: add particles via fair box generator - /// - - ///NOTE: the heavy fragements are not in the pdg table by default, so we have to add them: - TDatabasePDG* pdgDB = TDatabasePDG::Instance(); - // PDG nuclear states are 10-digit numbers - // 10LZZZAAAI e.g. deuteron is 1000010020 - // add IONS - const Int_t kion = 1000000000; - const Double_t khSlash = 1.0545726663e-27; - const Double_t kErg2Gev = 1 / 1.6021773349e-3; - const Double_t khShGev = khSlash * kErg2Gev; - const Double_t kYear2Sec = 3600 * 24 * 365.25; - // Done by default now from Pythia6 table - // check if already defined - - Int_t ionCode = kion + 10020; - if (!pdgDB->GetParticle(ionCode)) - pdgDB->AddParticle( - "Deuteron", "Deuteron", 1.875613, kTRUE, 0, 3, "Ion", ionCode); - pdgDB->AddAntiParticle("AntiDeuteron", -ionCode); - - ionCode = kion + 10030; - if (!pdgDB->GetParticle(ionCode)) - pdgDB->AddParticle("Triton", - "Triton", - 2.80925, - kFALSE, - khShGev / (12.33 * kYear2Sec), - 3, - "Ion", - ionCode); - pdgDB->AddAntiParticle("AntiTriton", -ionCode); - - ionCode = kion + 20030; - if (!pdgDB->GetParticle(ionCode)) - pdgDB->AddParticle("HE3", "HE3", 2.80923, kFALSE, 0, 6, "Ion", ionCode); - pdgDB->AddAntiParticle("AntiHE3", -ionCode); - - ionCode = kion + 20040; - if (!pdgDB->GetParticle(ionCode)) - pdgDB->AddParticle("Alpha", - "Alpha", - 3.727379, - kTRUE, - khShGev / (12.33 * kYear2Sec), - 6, - "Ion", - ionCode); - pdgDB->AddAntiParticle("AntiAlpha", -ionCode); - - // Use box generators for defined particles and per event multiplcities (pdg,mult) - FairBoxGenerator* boxGenDeu = new FairBoxGenerator(1000010020, mult); - boxGenDeu->SetPtRange(0., 3.); - boxGenDeu->SetPhiRange(0., 360.); - boxGenDeu->SetThetaRange(2.5, 25.); - boxGenDeu->SetCosTheta(); - boxGenDeu->SetDebug(kTRUE); - boxGenDeu->Init(); - primGen->AddGenerator(boxGenDeu); - - FairBoxGenerator* boxGenTri = new FairBoxGenerator(1000010030, mult); - boxGenTri->SetPtRange(0., 3.); - boxGenTri->SetPhiRange(0., 360.); - boxGenTri->SetThetaRange(2.5, 25.); - boxGenTri->SetCosTheta(); - boxGenTri->SetDebug(kTRUE); - boxGenTri->Init(); - primGen->AddGenerator(boxGenTri); - - FairBoxGenerator* boxGenHe3 = new FairBoxGenerator(1000020030, mult); - boxGenHe3->SetPtRange(0., 3.); - boxGenHe3->SetPhiRange(0., 360.); - boxGenHe3->SetThetaRange(2.5, 25.); - boxGenHe3->SetCosTheta(); - boxGenHe3->SetDebug(kTRUE); - boxGenHe3->Init(); - primGen->AddGenerator(boxGenHe3); - - FairBoxGenerator* boxGenAlp = new FairBoxGenerator(1000020040, mult); - boxGenAlp->SetPtRange(0., 3.); - boxGenAlp->SetPhiRange(0., 360.); - boxGenAlp->SetThetaRange(2.5, 25.); - boxGenAlp->SetCosTheta(); - boxGenAlp->SetDebug(kTRUE); - boxGenAlp->Init(); - primGen->AddGenerator(boxGenAlp); - - return; -} - -// ------------------------------------------------------------------------ -void AddCharmoniaCocktail(FairPrimaryGenerator* primGen, - Int_t idx, - TString inputFile) { - /// - /// CHARMONIA: add particles via pluto generator, decay them via TEvtGen - /// - - Int_t module = 10000; - - ///NOTE: this is for pAu collisisons simulations, either jpsi OR psi2S produced in the event - // TString idu=Form("%05d",idx+1); - // inputFile += idu + "/pluto.pau.30gev.jpsi.epem.root"; // OR "/pluto.pau.30gev.jpsi.root" - - Int_t idy = (idx >= 0 ? idx % module : gRandom->Integer(module)); - TString file = - inputFile + Form("%05d", idy + 1) + "/pluto.pau.30gev.jpsi.epem.root"; - printf("pluto vec.meson: %s \n", file.Data()); - CbmPlutoGenerator* plutoJpsi = new CbmPlutoGenerator(file); - if (idy % 2) primGen->AddGenerator(plutoJpsi); - - /// for radiative decay - CbmGenEvtGen* evtgen = new CbmGenEvtGen(); - evtgen->SetForceDecay(kBJpsiDiElectron); - evtgen->SetParticleSwitchedOff(CbmGenEvtGen::kJPsiPart); - if (idy % 2) primGen->AddGenerator(evtgen); - - /// NOTE: do NOT decay with evtgen - // idy = (idx >= 0 ? idx%module : gRandom->Integer(module) ); - file = inputFile + Form("%05d", idy + 1) + "/pluto.pau.30gev.jpsi.epem.root"; - file.ReplaceAll("jpsi", "psi2s"); - //inputFile.ReplaceAll("jpsi","psi2s.epem"); - printf("pluto vec.meson: %s \n", file.Data()); - CbmPlutoGenerator* plutoPsi2S = new CbmPlutoGenerator(file); - plutoPsi2S->SetManualPDG(100443); - if (!(idy % 2)) primGen->AddGenerator(plutoPsi2S); - - return; -} - -// ------------------------------------------------------------------------ -void AddBoxGenerator(FairPrimaryGenerator* primGen, Int_t mult) { - /// - /// Use box generators for defined particles and per event multiplcities (pdg,mult) - /// - - FairBoxGenerator* boxGenEle = new FairBoxGenerator(11, mult); - boxGenEle->SetPtRange(0., 3.); - boxGenEle->SetPhiRange(0., 360.); - boxGenEle->SetThetaRange(2.5, 25.); - boxGenEle->SetCosTheta(); - boxGenEle->SetDebug(kTRUE); - boxGenEle->Init(); - primGen->AddGenerator(boxGenEle); - - FairBoxGenerator* boxGenPiM = new FairBoxGenerator(-211, mult); - boxGenPiM->SetPtRange(0., 3.); - boxGenPiM->SetPhiRange(0., 360.); - boxGenPiM->SetThetaRange(2.5, 25.); - boxGenPiM->SetCosTheta(); - boxGenPiM->SetDebug(kTRUE); - boxGenPiM->Init(); - primGen->AddGenerator(boxGenPiM); - - FairBoxGenerator* boxGenPosi = new FairBoxGenerator(-11, mult); - boxGenPosi->SetPtRange(0., 3.); - boxGenPosi->SetPhiRange(0., 360.); - boxGenPosi->SetThetaRange(2.5, 25.); - boxGenPosi->SetCosTheta(); - boxGenPosi->SetDebug(kTRUE); - boxGenPosi->Init(); - primGen->AddGenerator(boxGenPosi); - - FairBoxGenerator* boxGenPiP = new FairBoxGenerator(+211, mult); - boxGenPiP->SetPtRange(0., 3.); - boxGenPiP->SetPhiRange(0., 360.); - boxGenPiP->SetThetaRange(2.5, 25.); - boxGenPiP->SetCosTheta(); - boxGenPiP->SetDebug(kTRUE); - boxGenPiP->Init(); - primGen->AddGenerator(boxGenPiP); - - return; -} - -// ------------------------------------------------------------------------ -void SetParticleDecays() { - /// - /// set and unset particle decays for geant3 - /// - TGeant3* gMC3 = (TGeant3*) gMC; - - /// switch off certain decays - gMC3->SetUserDecay(443); // Force the decay to be done w/external decayer - gMC3->SetUserDecay(100443); // Force the decay to be done w/external decayer - std::cout - << "----------------- switch OFF: jpsi decay by geant3 -----------------" - << std::endl; - - /// set branching ratios by hand - Float_t bratioEta[6]; - Int_t modeEta[6]; - - for (Int_t kz = 0; kz < 6; ++kz) { - bratioEta[kz] = 0.; - modeEta[kz] = 0; - } - - Int_t ipa = 17; - bratioEta[0] = 39.38; //2gamma - bratioEta[1] = 32.20; //3pi0 - bratioEta[2] = 22.70; //pi+pi-pi0 - bratioEta[3] = 4.69; //pi+pi-gamma - bratioEta[4] = 0.60; //e+e-gamma - bratioEta[5] = 4.4e-2; //pi02gamma - - modeEta[0] = 101; //2gamma - modeEta[1] = 70707; //3pi0 - modeEta[2] = 80907; //pi+pi-pi0 - modeEta[3] = 80901; //pi+pi-gamma - modeEta[4] = 30201; //e+e-gamma - modeEta[5] = 10107; //pi02gamma - gMC3->Gsdk(ipa, bratioEta, modeEta); - - Float_t bratioPi0[6]; - Int_t modePi0[6]; - - for (Int_t kz = 0; kz < 6; ++kz) { - bratioPi0[kz] = 0.; - modePi0[kz] = 0; - } - - ipa = 7; - bratioPi0[0] = 98.798; - bratioPi0[1] = 1.198; - - modePi0[0] = 101; - modePi0[1] = 30201; - - gMC3->Gsdk(ipa, bratioPi0, modePi0); - - /* - TDatabasePDG *pdg = TDatabasePDG::Instance(); - Int_t mother = 0; - Float_t br[10] = {0.}; - Int_t mode[10] = {0}; - - //////////// pi0 decay ////////// - mother = 111; - TParticlePDG *p1 = pdg->GetParticle(mother); - p1->Print(); //decay info - - // loop over all decay channels - TDecayChannel *dc=0x0; - for(Int_t i=0; i<10; i++) { br[i]=0.; mode[i]=0; } //reset - for(Int_t i=0; i<p1->NDecayChannels(); i++) { - dc = p1->DecayChannel(i); - br[i] = dc->BranchingRatio(); - mode[i] = 0; - // loop over all daughters - for(Int_t d=0; d<dc->NDaughters(); d++) { - mode[i] += pdg->ConvertPdgToGeant3(dc->DaughterPdgCode(d)) * TMath::Power(100,d); - } - } - for(Int_t i=0; i<10; i++) - Printf("i%d --> BR:%.3e \t mode:%d ",i,br[i],mode[i]); - gMC3->Gsdk(pdg->ConvertPdgToGeant3(mother), br, mode); - - - //////////// eta decay ////////// - mother = 221; - TParticlePDG *p1 = pdg->GetParticle(mother); - p1->Print(); //decay info - - // loop over all decay channels - TDecayChannel *dc=0x0; - for(Int_t i=0; i<10; i++) { br[i]=0.; mode[i]=0; } //reset - for(Int_t i=0; i<p1->NDecayChannels(); i++) { - dc = p1->DecayChannel(i); - br[i] = dc->BranchingRatio(); - mode[i] = 0; - // loop over all daughters - for(Int_t d=0; d<dc->NDaughters(); d++) { - mode[i] += pdg->ConvertPdgToGeant3(dc->DaughterPdgCode(d)) * TMath::Power(100,d); - } - } - for(Int_t i=0; i<10; i++) - Printf("i%d --> BR:%.3e \t mode:%d ",i,br[i],mode[i]); - gMC3->Gsdk(pdg->ConvertPdgToGeant3(mother), br, mode); - */ - - // set randomizer to gMC - TRandom3* rnd = new TRandom3(0); // computed via a TUUID object (unique!) - gMC->SetRandom(rnd); -} - -void ConfigureMCStack() { - /// - /// configure the cbm stack, apply/release cuts - /// - CbmStack* stack = dynamic_cast<CbmStack*>(gMC->GetStack()); - if (stack) { - // printf("[USER MODIFICATION]: set MC stack cuts by hand! \n"); - // stack->StoreSecondaries(kTRUE); // default: kTRUE - stack->SetMinPoints(0); // default: 1 - // stack->SetEnergyCut(0.00001); // default: 0. - // stack->StoreMothers(kTRUE); // default: kTRUE - } -} diff --git a/analysis/PWGDIL/dielectron/papaframework/scripts/run_testing.C b/analysis/PWGDIL/dielectron/papaframework/scripts/run_testing.C new file mode 100644 index 0000000000000000000000000000000000000000..7a59a250887761a063ef66f38ae90cdb3c98bbf2 --- /dev/null +++ b/analysis/PWGDIL/dielectron/papaframework/scripts/run_testing.C @@ -0,0 +1,193 @@ +R__ADD_INCLUDE_PATH($PWD) + +#include "Config_dilepton_testing.C" +TString cfgFunc = "Config_dilepton_"; +TString configName = "testing"; + +// Includes needed for IDE +#if !defined(__CLING__) +#include "CbmDefs.h" +#include "CbmMCDataManager.h" +#include "CbmSetup.h" + +#include "FairSystemInfo.h" +#include <FairFileSource.h> +#include <FairMonitor.h> +#include <FairParAsciiFileIo.h> +#include <FairParRootFileIo.h> +#include <FairRunAna.h> +#include <FairRuntimeDb.h> + +#include <TStopwatch.h> +#endif + + +void run_testing(Int_t nEvents = 0, Bool_t test = true, Int_t id = 1, TString TASK = "", + TString setup = "sis100_electron") +{ + + + // ----- Environment -------------------------------------------------- + TString myName = "run_analysis"; // this macro's name for screen output + TString srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory + // ------------------------------------------------------------------------ + + // ----- In- and output directory names ------------------------------- + TString inDir = gSystem->Getenv("INDIR"); + TString inFile = gSystem->Getenv("INFILE"); // input file list + TString outDir = gSystem->Getenv("OUTDIR"); + // if(outDir.IsNull()) outFile = "."; + + // --- Logger settings ---------------------------------------------------- + gErrorIgnoreLevel = kFatal; //kInfo, kWarning, kError, kFatal; + // ------------------------------------------------------------------------ + + + std::cout << std::endl; + std::cout << "-I- " << myName << ": Loading setup " << setup << std::endl; + CbmSetup::Instance()->LoadSetup(setup); + + // ----- run manager + I/O -------------------------------------------- + std::cout << std::endl; + FairRunAna* run = new FairRunAna(); + run->SetOutputFile(outDir + "/analysis.root"); + FairFileSource* src = NULL; + + /// stopwatch + TStopwatch timer; + timer.Start(); + + Int_t i = 0; + TString file; + // char filename[300]; + ifstream in(inFile); + TList* parList = new TList(); + TString parFile = ""; + + + CbmMCDataManager* mcManager = new CbmMCDataManager("MCDataManager", 1); + + TString traFile = ""; + /// loop over all file in list + if (test) { + ifstream in(inFile); + while (in >> file) { + // mc sim file + if (!i) src = new FairFileSource(file + "/tra.root"); + src->AddFriend(file + "/raw.root"); + src->AddFriend(file + "/rec.root"); + parFile = file + "/par.root"; + + mcManager->AddFile(traFile); + + i++; + } + } + if (!test) { + ifstream in(outDir + "/" + inFile); + while (in >> file) { + + if (!i) src = new FairFileSource(file + "/tra.root"); + src->AddFriend(file + "/raw.root"); + src->AddFriend(file + "/rec.root"); + parFile = file + "/par.root"; + // parFile = outDir + "/../../par.root"; + + mcManager->AddFile(traFile); + + i++; + } + } + + // parList->Dump(); + // add source to run + run->SetSource(src); + + // // ------------------------------------------------------------------------ + + // ----- MCDataManager (legacy mode) ----------------------------------- + run->AddTask(mcManager); + // ------------------------------------------------------------------------ + + // CbmMatchRecoToMC* match1 = new CbmMatchRecoToMC(); + // run->AddTask(match1); + + + // ----- L1/KF tracking + PID ----------------------------------------- + std::cout << std::endl; + std::cout << "-I- " << myName << ": Loading tasks " << std::endl; + + //CbmKF is needed for Extrapolation + CbmKF* kf = new CbmKF(); + run->AddTask(kf); + std::cout << "-I- : Added task " << kf->GetName() << std::endl; + + CbmL1* l1 = new CbmL1(); + // --- Material budget file names + if (CbmSetup::Instance()->IsActive(ECbmModuleId::kMvd)) { + TString geoTag; + CbmSetup::Instance()->GetGeoTag(ECbmModuleId::kMvd, geoTag); + TString matFile = gSystem->Getenv("VMCWORKDIR"); + matFile = matFile + "/parameters/mvd/mvd_matbudget_" + geoTag + ".root"; + std::cout << "Using material budget file " << matFile << std::endl; + l1->SetMvdMaterialBudgetFileName(matFile.Data()); + } + if (CbmSetup::Instance()->IsActive(ECbmModuleId::kSts)) { + TString geoTag; + CbmSetup::Instance()->GetGeoTag(ECbmModuleId::kSts, geoTag); + TString matFile = gSystem->Getenv("VMCWORKDIR"); + matFile = matFile + "/parameters/sts/sts_matbudget_" + geoTag + ".root"; + // matFile = matFile + "/parameters/sts/sts_matbudget_v19a.root"; + // matFile = matFile + "/parameters/sts/sts_matbudget_v19i.root"; + std::cout << "Using material budget file " << matFile << std::endl; + l1->SetStsMaterialBudgetFileName(matFile.Data()); + // } + run->AddTask(l1); + std::cout << "-I- : Added task " << l1->GetName() << std::endl; + } + + // --- TRD pid tasks + if (CbmSetup::Instance()->IsActive(ECbmModuleId::kTrd)) { + CbmTrdSetTracksPidLike* trdLI = new CbmTrdSetTracksPidLike("TRDLikelihood", "TRDLikelihood"); + trdLI->SetUseMCInfo(kTRUE); + trdLI->SetUseMomDependence(kTRUE); + run->AddTask(trdLI); + std::cout << "-I- : Added task " << trdLI->GetName() << std::endl; + // ------------------------------------------------------------------------ + } + + // ----- PAPA tasks --------------------------------------------------- + std::cout << std::endl; + std::cout << "-I- " << myName << ": Loading private tasks " << std::endl; + + TString cfgPath = outDir + "/../"; + TString cfgFile = cfgPath + cfgFunc + configName + ".C"; //cfgPath + cfgFunc + configName + ".C"; + TString cfgFunct = cfgFunc + configName + "(\"" + configName + "\")"; + TString cfg = cfgFunc + configName + "(" + std::to_string(id) + ",\"" + TASK + "\")"; + gROOT->LoadMacro((cfgPath + cfgFunct + configName + ".C")); + FairTask* task = reinterpret_cast<FairTask*>(gROOT->ProcessLine(cfg)); + run->AddTask(task); + + // // ------------------------------------------------------------------------ + // set parameter list + FairParRootFileIo* parIo1 = new FairParRootFileIo(); + // parIo1->open(parFile.Data(),"UPDATE"); + parIo1->open(parFile.Data(), "READ"); + // parIo1->open(parList); + FairRuntimeDb* rtdb = run->GetRuntimeDb(); + rtdb->setFirstInput(parIo1); + rtdb->setOutput(parIo1); + rtdb->saveOutput(); + + /// intialize and run + run->Init(); + run->Run(0, nEvents); + + timer.Stop(); + std::cout << "Macro finished succesfully." << std::endl; + std::cout << " Output file is " << (outDir + "analysis.root") << std::endl; + // std::cout << "Parameter file is " << parFile << std::endl; + std::cout << "Real time " << timer.RealTime() << " s, CPU time " << timer.CpuTime() << " s" << std::endl; + std::cout << " Test passed" << std::endl; + std::cout << " All ok " << std::endl; +} diff --git a/analysis/PWGDIL/dielectron/papaframework/scripts/simreco.sh b/analysis/PWGDIL/dielectron/papaframework/scripts/simreco.sh deleted file mode 100755 index e8cc8341d41f1504a122a52573bc3d379aab3b97..0000000000000000000000000000000000000000 --- a/analysis/PWGDIL/dielectron/papaframework/scripts/simreco.sh +++ /dev/null @@ -1,107 +0,0 @@ -#!/bin/bash - -LOCATION=/lustre/nyx - -## choose cbm root installation -source $LOCATION/cbm/users/$USER/CBMsoft/cbm-env.sh -n - -## go to directory with macros -export OUTDIR="$1/" -mkdir -p $OUTDIR -cd $OUTDIR -echo $PWD - -## copy rootrc -cp -v "../.rootrc" "." - -## get cbm setup macro forseen -export SETUP=$(basename $OUTDIR/../setup* .C) -SETUP=${SETUP/setup_/} -echo "cbm setup is: $SETUP" - -## number of events to process -NEVT=$2 - -## check files in output directory -echo "$OUTDIR" -ls -lhSra "$OUTDIR" - -## simulation -if [ -f "$1/run_sim_${SETUP}_ok" ] ; then - echo "mc file was produced already"; -else - echo "start mc generation"; - rm -v "$1/${SETUP}_mc.root" - rm -v "$1/${SETUP}_reco.root" - root -l -b -q "$PWD/../run_sim_new.C( $NEVT,\"$SETUP\")" - - # validate the root file - if [ -f /cvmfs/fairroot.gsi.de/cbm/bin/validate_cbmroot.sh ] ; then - export MACRODIR="/cvmfs/fairroot.gsi.de/cbm/macros" - /cvmfs/fairroot.gsi.de/cbm/bin/validate_cbmroot.sh "$1/${SETUP}_mc.root" $NEVT - result=$? - if [ "$result" = "0" ]; then - echo "mc production was succesfull!"; - else - echo "mc production was not succesfull!"; - rm -v "$1/run_sim_${SETUP}_ok" - rm -v "$1/${SETUP}_mc.root" - fi - fi - echo -fi - - -## reconstruction -if [ -f "$1/run_reco_${SETUP}_ok" ] ; then - echo "reco file was produced already"; -elif [ -f "$1/run_sim_${SETUP}_ok" ] ; then - echo "start reconstruction"; - rm -v "$1/${SETUP}_reco.root" - root -l -b -q "$PWD/../run_reco_new.C( $NEVT,\"$SETUP\")" - - # validate the root file - if [ -f /cvmfs/fairroot.gsi.de/cbm/bin/validate_cbmroot.sh ] ; then - export MACRODIR="/cvmfs/fairroot.gsi.de/cbm/macros" - /cvmfs/fairroot.gsi.de/cbm/bin/validate_cbmroot.sh "$1/${SETUP}_reco.root" $NEVT - result=$? - if [ "$result" = "0" ]; then - echo "reconstruction was succesfull!"; - else - echo "reconstruction was not succesfull!"; - rm -v "$1/run_reco_${SETUP}_ok" - rm -v "$1/${SETUP}_reco.root" - fi - fi - echo -else - echo "mc production was not succesfull!"; - rm -v "$1/${SETUP}_mc.root" -fi - - -## validation -if [ -f "$1/run_reco_${SETUP}_ok" ] ; then - echo "everything okay" - ls -lhSra "$1" - echo "clean up diretory" - rm -v "$1/.rootrc" - rm -v "$1/L1_histo.root" - rm -v "$1/TRhistos.root" - rm -v "$1/all_*.par" - rm -v "$1/gphysi.dat" - ls -lhSra "$1" - echo "everything done"; -else - echo "validation failed!" - rm -v "$1/${SETUP}_reco.root" - if [ -f "$1/run_sim_${SETUP}_ok" ] ; then - echo - ls -lhSra "$1" - echo "simulation okay"; - else - rm -v "$1/${SETUP}_mc.root" - fi -fi - - diff --git a/analysis/PWGDIL/dielectron/papaframework/scripts/submit_analysis.sh b/analysis/PWGDIL/dielectron/papaframework/scripts/submit_analysis.sh index f0a29b9b9e9411ab2a3003f93a4596c5da9bbbea..b958dd1f9d6a31012ac72d57e9c47a2c627f5616 100755 --- a/analysis/PWGDIL/dielectron/papaframework/scripts/submit_analysis.sh +++ b/analysis/PWGDIL/dielectron/papaframework/scripts/submit_analysis.sh @@ -1,85 +1,162 @@ #!/bin/bash -## test train or batch submission -TEST=1 - -## split level: defines number of files per jobs -split_level=10; +# indir need to be simulation directory with a subdirectory structure that contains the individual simulation runs (typically numbers) +indir=/home/$USER/sim/test/ +outname="test" # the output will be written in this directory under this directory name +CONFIG=1 # this one can be used to switch between data naming and configurations ; it is parsed to the next level script and changes the run_script that is called +train="test" # inside the output directory is a another directory for the specific analysis + +# CONFIG 1 - default dilepton testing config - sim names: tra.root/raw.root/rec.root +# CONFIG 2 - default dilepton testing config for common productions - sim names: 1.tra.root/1.event.raw.root/1.reco.root, etc + +out=$PWD/$outname/$train +## test train or batch submissio1 +TEST=1 # TEST 1 - calculates a testrun of one file locally ; TEST 0 - submits jobs to virgo +GROUP=10 # Number of files that is processed per job +DEBUG=0 # Single small job (50 events) that is send to the DEBUG partition to test that everything runs smoothly +useDEB=0 # For a very small number of very small jobs ; they can be processed on the DEBUG partition in rare cases ## number of events: 0=all NEVENT=0 +STOPJOB=20000 +MINJOB=0 -## unique train id -train="train0001" +if [[ $TEST == "1" ]] ; then + train="test" +fi ## storage element and directory LOCATION=/lustre/nyx DIR=sim +split_level=1; + ## test parameters if [[ $TEST == "1" ]] ; then train=${train/train/test}; + train="test" + GROUP=1; split_level=1; - NEVENT=100 + NEVENT=0 fi - -## input directory of simulated and reconstructed data -indir=$LOCATION/cbm/users/$USER/$DIR/sim_AA_UrQMD_eeCocktail_centr010_JUN16_25mum_4lay_wMVD - -## output directory -out=$indir/$train mkdir -p $out -## path for analysis configs -configPath="$LOCATION/cbm/users/$USER/papa-conf" - ## get CBM setup used in simulation and reconstruction SETUP=$(basename $indir/setup* .C) SETUP=${SETUP/setup_/} ## copy user config and analysis macro to train output directory ## for backup and documentation -cp -v "$configPath/"*.C "$out/." -cp -v "$PWD/run_analysis.C" "$out/." +if [[ $CONFIG == "1" ]] ; then + rm -v "$out/run_testing.C" + cp -v "$PWD/run_testing.C" "$out/." + cp -v "$PWD/Config_dilepton_testing.C" "$out/." +fi +if [[ $CONFIG == "2" ]] ; then + rm -v "$out/run_common_analysis.C" + cp -v "$PWD/run_common_analysis.C" "$out/." + cp -v "$PWD/Config_dilepton_testing.C" "$out/." +fi ## get run list and sort it -if grep -Fq "$indir/" runListSort.txt; then +#if grep -Fq "." runList.txt; then +if grep -Fq "." never.txt; then echo "sorted run list already there, only split the list"; + sort runList.txt > runListSort.txt else ## only collect file with validated reconstruction output if [[ $TEST == "1" ]] ; then - find $indir -type f -name "run_reco_${SETUP}_ok" -printf "%h\n" | head -n $split_level > runList.txt; + #find $indir -type f -name "run_transport_${SETUP}_ok" -printf "%h\n" | head -n $split_level > runList.txt; + if [[ $CONFIG == "1" ]] ; then + if [[ $GROUP == "1" ]] ; then + find $indir -type f -name "rec.root" -size +100M -printf "%h\n" | head -n $split_level > runList.txt; + else + find $indir -type f -name "rec.root" -size +100M -exec dirname {} \; > runList.txt; + fi + fi + if [[ $CONFIG == "2" ]] ; then + if [[ $GROUP == "1" ]] ; then + find $indir -type f -name "*.rec.root" -size +100M -printf "%h\n" | head -n $split_level > runList.txt; + else + find $indir -type f -name "*.rec.root" -size +100M -exec dirname {} \; > runList.txt; + fi + fi else - find $indir -type f -name "run_reco_${SETUP}_ok" -exec dirname {} \; > runList.txt; + if [[ $CONFIG == "1" ]] ; then + find $indir -type f -name "rec.root" -size +100M -exec dirname {} \; > runList.txt; + fi + if [[ $CONFIG == "2" ]] ; then + find $indir -type f -name "*.rec.root" -size +100M -exec dirname {} \; > runList.txt; + fi fi + echo ""; + echo ""; + echo ""; + echo "sorting"; + echo ""; + echo ""; + echo ""; sort runList.txt > runListSort.txt fi ## split file list by split_level into many files (1 list per job) split -a 4 -l $split_level runListSort.txt filelist_${train}_ -ListOfFilelist=`ls filelist_${train}_*` +ListOfFilelist=`ls filelist_*` # loop over all lists I=0 -for filelist in $ListOfFilelist ; do - let I=$I+1 +N=1 +COUNT=0 + +out1="" +file1="" +out2="" +file2="" +out3="" +file3="" + + +# echo $ListOfFilelist +for filelist in $ListOfFilelist ; do + if [[ "$I" -lt "$MINJOB" ]] ; then + let I=$I+1 + continue + fi + ## create output directory for job $I - outdir=$out/`printf "%04d/" $I` - echo "-$outdir-" - mkdir -p $outdir + outdir=$out/`printf "%05d/" $I` + mkdir -p $outdir + + if [[ "$GROUP" -gt "1" ]] ; then + outlist[$COUNT]=$outdir + file[$COUNT]=$filelist + let COUNT=$COUNT+1 + fi ## test analysis (abort after first job) if [[ $TEST == "1" ]] ; then - # clean up - rm $(ls filelist_${train}_* | grep -v aaaa) - rm -v runListSort.txt; - ./analysis.sh $indir $filelist $outdir $NEVENT # LOCAL run - rm -v filelist_${train}_* - break; - + if [[ $GROUP == "1" ]] ; then + ./analysis.sh $indir $filelist $outdir $NEVENT $CONFIG $TEST $I # LOCAL run + rm $(ls filelist_${train}_* | grep -v aaaa) + rm -v runListSort.txt; + + break + else + if [[ $N == "$GROUP" ]] ; then + echo ${outlist[*]} + echo ${file[*]} + ./ana.sh $indir $NEVENT $CONFIG $GROUP $TEST $I ${file[*]} ${outlist[*]} # LOCAL run + outdir="" + file="" + COUNT=0 + N=0 + rm -v runListSort.txt; + break + fi + fi else #### gsi batch environement for kronos (slurm) #### #### NOTE: possible partitions: sinfo -o "%9P %6g %10l %5w %5D %13C %N" @@ -90,32 +167,65 @@ for filelist in $ListOfFilelist ; do #### memory max 4096M, default 2048M ## configure batch job - ## copy filelist mv $filelist $outdir/. ## partition - resources="--mem=4000 --time=0-03:00:00" + resources="--mem=4000 --time=0-08:00:00" partition="--partition=main" -# resources="--mem=4000 --time=0-00:05:00" -# partition="--partition=debug" + if [[ $DEBUG == "1" ]] ; then + partition="--partition=debug" + NEVENT=50 + resources="--mem=4000 --time=0-00:20:00" + fi + + if [[ $useDEB == "1" ]] ; then + partition="--partition=debug" + resources="--mem=4000 --time=0-00:20:00" + fi + scriptdir="$PWD" workdir="--workdir ${outdir}" - - mail_type="--mail-type=FAIL" - mail_user="--mail-user=j.book" - - job_name="--job-name=ana_${train}_${I}" + workdir="--chdir ${outdir}" + + job_name="--job-name=${I}" log_output="--output=${outdir}/ana.slurm.out" log_error="--error=${outdir}/ana.slurm.err" ## submit job - command="${partition} ${resources} ${workdir} ${log_output} ${log_error} ${mail_type} ${mail_user} ${job_name} ${scriptdir}/analysis.sh $indir $filelist $outdir $NEVENT" - sbatch $command - break; - + if [[ $GROUP == "1" ]] ; then + command="${partition} ${resources} ${workdir} ${log_output} ${log_error} ${job_name} -- ${scriptdir}/analysis.sh $indir $filelist $outdir $NEVENT $CONFIG $TEST $I" + sbatch $command + if [[ $DEBUG == "1" ]] ; then + break + fi + else + if [[ $N == "$GROUP" ]] ; then + command="${partition} ${resources} ${workdir} ${log_output} ${log_error} ${job_name} -- ${scriptdir}/ana.sh $indir $NEVENT $CONFIG $GROUP $TEST $I ${file[*]} ${outlist[*]} " + echo $command + sbatch $command + outdir="" + file="" + COUNT=0 + N=0 + + if [[ $I == $STOPJOB ]] ; then + break; + fi + + if [[ $DEBUG == "1" ]] ; then + break + fi + fi + fi fi + if [[ $I == $STOPJOB ]] ; then + break; + fi + + let I=$I+1 + let N=$N+1 done diff --git a/analysis/PWGDIL/dielectron/papaframework/scripts/submit_merging.sh b/analysis/PWGDIL/dielectron/papaframework/scripts/submit_merging.sh index 9b442e28c93f32cdfd06d143b13a656c165a020d..acf4794c447e12defbedac3496d9bedd8785206d 100755 --- a/analysis/PWGDIL/dielectron/papaframework/scripts/submit_merging.sh +++ b/analysis/PWGDIL/dielectron/papaframework/scripts/submit_merging.sh @@ -2,78 +2,118 @@ ## split level: defines number of files per jobs split_level=10; +MERGEIT=$1 -## split step: iteration -splitstep=1 - -## unique train id -train="train0001" +TEST=0 ## storage element and directory LOCATION=/lustre/nyx DIR=sim -## output directory of analyzed data -indir=$LOCATION/cbm/users/$USER/$DIR/sim_AA_UrQMD_eeCocktail_centr010_JUN16_25mum_4lay_wMVD/$train - -## output directory -outdir=$indir/merge -mkdir -p $outdir - -## get run list and sort it -if [ -e $outdir/runList_$splitstep.txt ]; then - echo "run list already there, only split the list"; -elif [ $splitstep == "1" ]; then - # first iteration - find $indir/. -type f -size +100k -name "sis100_*analysis.root" > "$outdir/runList_$splitstep.txt"; -else - # iterations > 1 - laststep=$(expr $splitstep - 1) - find $outdir/. -type f -size +100k -name "merged_filelist_temp_$laststep*.root" > "$outdir/runList_$splitstep.txt"; +dirs=" +$PWD/test/test/ +" + +echo $dirs + +if [[ $MERGEIT == "1" ]] ; then + for direc in $dirs ; do + echo "" + echo "$direc/merge" + cd "$direc/merge" + hadd merged.root merged_* + rm merge1_* + rm filelist_temp_1a* + rm slurm-* + rm merged_* + rm -v ../merge.C + cd "-" + done + exit 0 fi -## clean up -rm $outdir/filelist_temp_* -rm "$outdir/merged_filelist_temp_${splitstep}*.root" -rm "$outdir/merge${splitstep}_*.slurm.*" - -## split file list by split_level into many files (1 list per job) -split -a 3 -l $split_level "$outdir/runList_$splitstep.txt" "$outdir/filelist_temp_$splitstep" -ListOfFilelist=`ls $outdir/filelist_temp_*` - - -# loop over all lists -I=0 -for filelist in $ListOfFilelist ; do - let I=$I+1 - - #### gsi batch environement for kronos (slurm) #### - #### NOTE: possible partitions: sinfo -o "%9P %6g %10l %5w %5D %13C %N" - #### and : scontrol show partition main - #### long: upto 7d (CPUs=3880) - #### main upto 8h (CPUs=7840) - #### debug: upto 20m (CPUs=200) - #### memory max 4096M, default 2048M - - ## configure batch job - tmpjobscriptname=$PWD/merge_${train}_${splitstep}_$I.sh - #### - echo "#! /bin/sh" > $tmpjobscriptname - echo "#SBATCH --output=$outdir/merge${splitstep}_${I}.slurm.out" >> $tmpjobscriptname - echo "#SBATCH --error=$outdir/merge${splitstep}_${I}.slurm.err" >> $tmpjobscriptname - echo "#SBATCH --job-name=m${filelist: -4}_${train}" >> $tmpjobscriptname - echo "#SBATCH --mail-type=FAIL" >> $tmpjobscriptname - echo "#SBATCH --mail-user=j.book@gsi.de" >> $tmpjobscriptname - echo "#SBATCH --mem-per-cpu=2000" >> $tmpjobscriptname - echo "#SBATCH --time=0-00:20:00" >> $tmpjobscriptname - echo "#SBATCH --partition=debug" >> $tmpjobscriptname - #### - echo "srun ./merging.sh $filelist" >> $tmpjobscriptname - - chmod u+x $tmpjobscriptname - - ## submit job - sbatch $tmpjobscriptname - +if [[ $MERGEIT == "2" ]] ; then + for direc in $dirs ; do + echo "" + echo "$direc/merge" + cd "$direc/merge" + hadd merged.root merged_* + rm merge1_* + rm filelist_temp_1a* + rm merged_* + rm -rv ../0* + rm -v ../merge.C + cd "-" + done + exit 0 +fi +for direc in $dirs ; do + indir=$direc + + cp -v merge.C $indir + ## output directory + outdir=$indir/merge + rm -rv $outdir + mkdir -p $outdir + + echo "- indir $indir" + echo "- outdir $outdir" + + ## get run list and sort it + if [ -e $outdir/runList_$splitstep.txt ]; then + echo "run list already there, only split the list"; + elif [ $splitstep == "1" ]; then + # first iteration + echo "- check -" + find $indir/. -type f -name "analysis.root" > "$outdir/runList_$splitstep.txt"; + else + # iterations > 1 + laststep=$(expr $splitstep - 1) + find $outdir/. -type f -size +100k -name "merged_filelist_temp_$laststep*.root" > "$outdir/runList_$splitstep.txt"; + fi + + ## clean up + rm $outdir/filelist_temp_* + rm "$outdir/merged_filelist_temp_${splitstep}*.root" + rm "$outdir/merge${splitstep}_*.slurm.*" + + ## split file list by split_level into many files (1 list per job) + split -a 3 -l $split_level "$outdir/runList_$splitstep.txt" "$outdir/filelist_temp_$splitstep" + ListOfFilelist=`ls $outdir/filelist_temp_*` + + + # loop over all lists + I=0 + for filelist in $ListOfFilelist ; do + let I=$I+1 + + #### gsi batch environement for kronos (slurm) #### + #### NOTE: possible partitions: sinfo -o "%9P %6g %10l %5w %5D %13C %N" + #### and : scontrol show partition main + #### long: upto 7d (CPUs=3880) + #### main upto 8h (CPUs=7840) + #### debug: upto 20m (CPUs=200) + #### memory max 4096M, default 2048M + + if [[ $TEST == "1" ]] ; then + ./merging.sh $filelist + break + fi + partition="--partition=debug" + resources="--mem=4000 --time=0-00:20:00" + scriptdir="$PWD" + workdir="--chdir ${outdir}" + + job_name="--job-name=${I}" + log_output="--output=${outdir}/ana.slurm.out" + log_error="--error=${outdir}/ana.slurm.err" + + command="${partition} ${resources} ${workdir} ${job_name} -- ${scriptdir}/merging.sh $filelist" + sbatch $command + + + done + rm ${PWD}/merge_${train}_* done + diff --git a/analysis/PWGDIL/dielectron/papaframework/scripts/submit_simreco.sh b/analysis/PWGDIL/dielectron/papaframework/scripts/submit_simreco.sh deleted file mode 100755 index ba670c6aa0a567c1f029cef17c05149597a85f1f..0000000000000000000000000000000000000000 --- a/analysis/PWGDIL/dielectron/papaframework/scripts/submit_simreco.sh +++ /dev/null @@ -1,89 +0,0 @@ -#!/bin/bash - -## test train or batch submission -TEST=1 - -## number of events: default=1000 -NEVENT=1000 - -## number of runs -NRUNS=2000 - -## storage element and directory -LOCATION=/lustre/nyx -DIR=sim - -## test parameters -if [[ $TEST == "1" ]] ; then - DIR="TEST${DIR}" - NEVENT=2 - NRUNS=1 -fi - - -## output directory for simulated and reconstructed data -outdir=$LOCATION/cbm/users/$USER/$DIR/sim_AA_UrQMD_eeCocktail_centr010_JUN16_25mum_4lay_wMVD -mkdir -p "$outdir" - -## source cbm root environment -source $LOCATION/cbm/users/$USER/CBMsoft/cbm-env.sh -n - -## choose CBM setup for simulation and reconstruction -export SETUP="sis100_electron" #"sis100_muon_jpsi" - -## copy cbm setup macro and simulation and reconstruction macros to output directory -## for backup and documentation -cp -nv "$VMCWORKDIR/macro/run/.rootrc" "$outdir/." -cp -nv "$VMCWORKDIR/geometry/setup/setup_${SETUP}.C" "$outdir/." -cp -nv "./run_sim_new.C" "$outdir/." -cp -nv "./run_reco_new.C" "$outdir/." - -## loop over all runs -I=0 -while [ "$I" -lt "$NRUNS" ]; do - - ## create output directory for job $I - out=$outdir/`printf "%05d" $I` - mkdir -p "$out" - - ## test analysis (abort after first job) - if [[ $TEST == "1" ]] ; then - - ./simreco.sh $out $NEVENT # LOCAL run - break; - - else - #### gsi batch environement for kronos (slurm) #### - #### NOTE: possible partitions: sinfo -o "%9P %6g %10l %5w %5D %13C %N" - #### and : scontrol show partition main - #### long: upto 7d (CPUs=3880) - #### main upto 8h (CPUs=7840) - #### debug: upto 20m (CPUs=200) - #### memory max 4096M, default 2048M - - ## configure batch job - - ## partition - resources="--mem=4000 --time=1-00:00:00" - partition="--partition=long" - - scriptdir="$PWD" - workdir="--workdir ${out}" - - mail_type="--mail-type=FAIL" - mail_user="--mail-user=j.book" - - job_name="--job-name=sr_${I}" - log_output="--output=${out}/simreco.slurm.out" - log_error="--error=${out}/simreco.slurm.err" - - ## submit job - command="${partition} ${resources} ${workdir} ${log_output} ${log_error} ${mail_type} ${mail_user} ${job_name} ${scriptdir}/simreco.sh $out $NEVENT" - sbatch $command - - fi - - let I=$I+1 - -done - diff --git a/analysis/PWGDIL/dielectron/papaframework/scripts/submit_urqmd.sh b/analysis/PWGDIL/dielectron/papaframework/scripts/submit_urqmd.sh deleted file mode 100755 index 208f9100225edba4445db21023f723e94c825360..0000000000000000000000000000000000000000 --- a/analysis/PWGDIL/dielectron/papaframework/scripts/submit_urqmd.sh +++ /dev/null @@ -1,86 +0,0 @@ -#!/bin/bash - -## test train or batch submission -TEST=1 - -## number of events: default=1000 -NEVENT=1000 - -## number of runs -NRUNS=10000 - -## storage element and directory -LOCATION=/lustre/nyx -DIR=urqmd - -## test parameters -if [[ $TEST == "1" ]] ; then - DIR="TEST${DIR}" - NEVENT=2 - NRUNS=1 -fi - - -## output directory for simulated and reconstructed data -outdir=$LOCATION/$USER/cbm/$DIR/auau/8gev/centr -mkdir -p "$outdir" - -## configuration file to output directory -## for backup and documentation -cp -nv "inputfile_AuAu_8" "$outdir/input" - -## set number of events per job -sed -i 's/nev 1000/nev '$NEVENT'/' "$outdir/input" - -## loop over all runs -I=0 -while [ "$I" -lt "$NRUNS" ]; do - - ## create output directory for job $I - out=$outdir/`printf "%05d" $I` - mkdir -p "$out" - - ## copy input file for each job and set random seed - cp "$outdir/input" "$out/." - sed -i 's/rsd 1/rsd '$I'/' "$out/input" - - ## test analysis (abort after first job) - if [[ $TEST == "1" ]] ; then - - ./urqmd.sh $out $NEVENT # LOCAL run - break; - - else - #### gsi batch environement for kronos (slurm) #### - #### NOTE: possible partitions: sinfo -o "%9P %6g %10l %5w %5D %13C %N" - #### and : scontrol show partition main - #### long: upto 7d (CPUs=3880) - #### main upto 8h (CPUs=7840) - #### debug: upto 20m (CPUs=200) - #### memory max 4096M, default 2048M - - ## configure batch job - tmpjobscriptname=$PWD/sr_$(basename $outdir)_$I.sh - #### - echo "#! /bin/sh" > $tmpjobscriptname - echo "#SBATCH --output=$out/urqmd.slurm.out" >> $tmpjobscriptname - echo "#SBATCH --error=$out/urqmd.slurm.err" >> $tmpjobscriptname - echo "#SBATCH --job-name=sr_${I}" >> $tmpjobscriptname - echo "#SBATCH --mail-type=END" >> $tmpjobscriptname - echo "#SBATCH --mail-user=j.book" >> $tmpjobscriptname - echo "#SBATCH --mem-per-cpu=4000" >> $tmpjobscriptname - echo "#SBATCH --time=1-05:00:00" >> $tmpjobscriptname - echo "#SBATCH --partition=long" >> $tmpjobscriptname - #### - echo "srun ./urqmd.sh $out $NEVENT" >> $tmpjobscriptname - - chmod u+x $tmpjobscriptname - - ## submit job - sbatch $tmpjobscriptname - fi - - let I=$I+1 - -done - diff --git a/analysis/PWGDIL/dielectron/papaframework/scripts/urqmd.sh b/analysis/PWGDIL/dielectron/papaframework/scripts/urqmd.sh deleted file mode 100755 index 94dcb2126776badad26d98cb82a6a3c7c8d1ba61..0000000000000000000000000000000000000000 --- a/analysis/PWGDIL/dielectron/papaframework/scripts/urqmd.sh +++ /dev/null @@ -1,66 +0,0 @@ -#!/bin/bash - -LOCATION=/lustre/nyx - -## choose cbm root installation -source $LOCATION/cbm/users/$USER/CBMsoft/cbm-env.sh -n - -## export variables needed by urqmd -export ftn09="$1/input" -export ftn13="/tmp/urqmd.f13" -export ftn14="$1/urqmd.f14" -export ftn15="/tmp/urqmd.f15" -export ftn16="/tmp/urqmd.f16" -export ftn19="/tmp/urqmd.f19" -export ftn20="/tmp/urqmd.f20" -echo "inputfile: $ftn09" -echo "outputfile: $ftn14" - -## number of events to process -NEVT=$2 - -## check files in output directory -echo "$1" -ls -lhSr "$1" - -## run urqmd -if [[ -f "$1/urqmd.f14" ]] || [[ -f "$1/urqmd.root" ]] ; then - echo "urqmd file was produced already"; -else - echo "start urqmd generation"; - rm -rfv "$1/urqmd.f14" - rm -rfv "$1/urqmd.root" - time urqmd.$(uname -i) -fi - - -## run unigen -if [ -f "$1/urqmd.root" ] ; then - echo "unigen file was produced already"; -elif [ -f "$1/urqmd.f14" ] ; then - echo "start unigen"; - time urqmd2u "$1/urqmd.f14" "$1/urqmd.root" $NEVT - echo -else - echo "unigen production was not succesfull!"; -fi - -## validation -if [ -f "$1/urqmd.root" ] ; then - echo - ls -lhSr "$1" - echo "clean up done"; - rm -v "$1/urqmd.f14" - rm -v "$1/input" - echo "everything done"; -else - echo "validation failed!" - rm -rfv "$1/urqmd.root" - if [ -f "$1/urqmd.f14" ] ; then - echo - ls -lhSr "$1" - echo "only urqmd production okay"; - fi -fi - -