diff --git a/analysis/PWGDIL/dimuon/CbmAnaDimuonAnalysis.cxx b/analysis/PWGDIL/dimuon/CbmAnaDimuonAnalysis.cxx index 96c9ad30d118ba585f15ba3873dfb912a332c576..c29c42569ac46d38d1d18ac7418b68d4aea87720 100644 --- a/analysis/PWGDIL/dimuon/CbmAnaDimuonAnalysis.cxx +++ b/analysis/PWGDIL/dimuon/CbmAnaDimuonAnalysis.cxx @@ -1,52 +1,52 @@ //---------------------------------------- // -// 2019 A. Senger a.senger@gsi.de +// 2019 A. Senger a.senger@gsi.de // //---------------------------------------- -#define PBINNING 100,0.,20. -#define THETABINNING 180,0.,180. -#define PTBINNING 100,0.,10. -#define YBINNING 100,-2.,6. -#define MBINNING 400,0.,4. +#define PBINNING 100, 0., 20. +#define THETABINNING 180, 0., 180. +#define PTBINNING 100, 0., 10. +#define YBINNING 100, -2., 6. +#define MBINNING 400, 0., 4. #include "CbmAnaDimuonAnalysis.h" -#include "CbmMCTrack.h" -#include "CbmStsTrack.h" -#include "CbmTrackMatchNew.h" -#include "CbmMuchTrack.h" -#include "CbmTrdTrack.h" -#include "CbmTrdHit.h" -#include "CbmTrackMatch.h" -#include "CbmMuchPixelHit.h" #include "CbmAnaMuonCandidate.h" #include "CbmGlobalTrack.h" +#include "CbmGlobalTrackFitterKF.h" +#include "CbmKFTrack.h" +#include "CbmKFVertex.h" +#include "CbmMCTrack.h" +#include "CbmMatch.h" #include "CbmMuchGeoScheme.h" +#include "CbmMuchPixelHit.h" +#include "CbmMuchTrack.h" +#include "CbmStsKFTrackFitter.h" +#include "CbmStsTrack.h" #include "CbmTofHit.h" #include "CbmTofPoint.h" -#include "CbmMatch.h" -#include "CbmVertex.h" -#include "CbmKFTrack.h" -#include "CbmStsKFTrackFitter.h" -#include "CbmGlobalTrackFitterKF.h" +#include "CbmTrackMatch.h" +#include "CbmTrackMatchNew.h" +#include "CbmTrdHit.h" +#include "CbmTrdTrack.h" #include "CbmTrdTrackFitterKF.h" -#include "CbmKFVertex.h" +#include "CbmVertex.h" #include "FairRootManager.h" #include "FairTrackParam.h" #include "TClonesArray.h" +#include "TDirectory.h" +#include "TFile.h" +#include "TH2D.h" +#include "TH3D.h" +#include "TLorentzVector.h" +#include "TMCProcess.h" #include "TMath.h" #include "TMultiLayerPerceptron.h" +#include "TProfile.h" #include "TSystem.h" #include "TTree.h" -#include "TLorentzVector.h" #include "TVector3.h" -#include "TH2D.h" -#include "TH3D.h" -#include "TProfile.h" -#include "TMCProcess.h" -#include "TFile.h" -#include "TDirectory.h" #include "PParticle.h" @@ -55,341 +55,488 @@ using std::vector; #include <sys/stat.h> -#include <iostream> #include <fstream> +#include <iostream> using namespace std; // ----- Default constructor ------------------------------------------- CbmAnaDimuonAnalysis::CbmAnaDimuonAnalysis(TString name, TString setup) - : FairTask("AnaDimuonAnalysis"), - fEvent(0), - fMCTracks(NULL), - fStsTracks(NULL), - fStsTrackMatches(NULL), - fMuchTracks(NULL), - fMuchTrackMatches(NULL), - fGlobalTracks(NULL), - fTrdTracks(NULL), - fTofHit(NULL), - fMuPlus(new TClonesArray("CbmAnaMuonCandidate",1)), - fMuMinus(new TClonesArray("CbmAnaMuonCandidate",1)), - fParticles(NULL), - fInputTree(NULL), - fPlutoFile(NULL), - fFitter(NULL), - fFitterTRD(NULL), - fFitterGlobal(NULL), - fVertex(NULL), - fChi2StsCut(2), - fChi2MuchCut(3), - fChi2VertexCut(3), - fAnnCut(-1), - fNeurons(0), - fSigmaTofCut(2), - fMass(0.105658), - fUseCuts(kTRUE), - fUseMC(kTRUE), - fNofMuchCut(11), - fNofStsCut(7), - fNofTrdCut(1), -// fFileName("histo.root"), -// fEffFileName("eff_histo.root"), - fGeoScheme(NULL) + : FairTask("AnaDimuonAnalysis") + , fEvent(0) + , fMCTracks(NULL) + , fStsTracks(NULL) + , fStsTrackMatches(NULL) + , fMuchTracks(NULL) + , fMuchTrackMatches(NULL) + , fGlobalTracks(NULL) + , fTrdTracks(NULL) + , fTofHit(NULL) + , fMuPlus(new TClonesArray("CbmAnaMuonCandidate", 1)) + , fMuMinus(new TClonesArray("CbmAnaMuonCandidate", 1)) + , fParticles(NULL) + , fInputTree(NULL) + , fPlutoFile(NULL) + , fFitter(NULL) + , fFitterTRD(NULL) + , fFitterGlobal(NULL) + , fVertex(NULL) + , fChi2StsCut(2) + , fChi2MuchCut(3) + , fChi2VertexCut(3) + , fAnnCut(-1) + , fNeurons(0) + , fSigmaTofCut(2) + , fMass(0.105658) + , fUseCuts(kTRUE) + , fUseMC(kTRUE) + , fNofMuchCut(11) + , fNofStsCut(7) + , fNofTrdCut(1) + , + // fFileName("histo.root"), + // fEffFileName("eff_histo.root"), + fGeoScheme(NULL) { - fPlutoFileName=name; - fSetupName =setup; + fPlutoFileName = name; + fSetupName = setup; } // ------------------------------------------------------------------------- // ----- SetParContainers ------------------------------------------------- -void CbmAnaDimuonAnalysis::SetParContainers() -{ -} +void CbmAnaDimuonAnalysis::SetParContainers() {} // ------------------------------------------------------------------------- // ----- Public method Init (abstract in base class) -------------------- -InitStatus CbmAnaDimuonAnalysis::Init() -{ +InitStatus CbmAnaDimuonAnalysis::Init() { // Get and check FairRootManager FairRootManager* fManager = FairRootManager::Instance(); - fMCTracks = (TClonesArray*) fManager->GetObject("MCTrack"); - if(nullptr == fMCTracks) LOG(FATAL) << "No MCTrack in input"; + fMCTracks = (TClonesArray*) fManager->GetObject("MCTrack"); + if (nullptr == fMCTracks) LOG(FATAL) << "No MCTrack in input"; - fStsTracks = (TClonesArray*) fManager->GetObject("StsTrack"); - if(nullptr == fStsTracks) LOG(FATAL) << "No StsTrack in input"; + fStsTracks = (TClonesArray*) fManager->GetObject("StsTrack"); + if (nullptr == fStsTracks) LOG(FATAL) << "No StsTrack in input"; - if(fUseMC)fStsTrackMatches = (TClonesArray*) fManager->GetObject("StsTrackMatch"); - if(nullptr == fStsTrackMatches && fUseMC) LOG(FATAL) << "No StsTrackMatch in input"; + if (fUseMC) + fStsTrackMatches = (TClonesArray*) fManager->GetObject("StsTrackMatch"); + if (nullptr == fStsTrackMatches && fUseMC) + LOG(FATAL) << "No StsTrackMatch in input"; - fMuchTracks = (TClonesArray*) fManager->GetObject("MuchTrack"); - if(nullptr == fMuchTracks) LOG(FATAL) << "No MuchTrack in input"; + fMuchTracks = (TClonesArray*) fManager->GetObject("MuchTrack"); + if (nullptr == fMuchTracks) LOG(FATAL) << "No MuchTrack in input"; - if(fUseMC)fMuchTrackMatches = (TClonesArray*) fManager->GetObject("MuchTrackMatch"); - if(nullptr == fMuchTrackMatches && fUseMC) LOG(FATAL) << "No MuchTrackMatch in input"; + if (fUseMC) + fMuchTrackMatches = (TClonesArray*) fManager->GetObject("MuchTrackMatch"); + if (nullptr == fMuchTrackMatches && fUseMC) + LOG(FATAL) << "No MuchTrackMatch in input"; - fGlobalTracks = (TClonesArray*) fManager->GetObject("GlobalTrack"); - if(nullptr == fGlobalTracks) LOG(FATAL) << "No GlobalTrack in input"; + fGlobalTracks = (TClonesArray*) fManager->GetObject("GlobalTrack"); + if (nullptr == fGlobalTracks) LOG(FATAL) << "No GlobalTrack in input"; - fTrdTracks = (TClonesArray*) fManager->GetObject("TrdTrack"); - if(nullptr == fTrdTracks) LOG(FATAL) << "No TrdTrack in input"; + fTrdTracks = (TClonesArray*) fManager->GetObject("TrdTrack"); + if (nullptr == fTrdTracks) LOG(FATAL) << "No TrdTrack in input"; - fTofHit = (TClonesArray*) fManager->GetObject("TofHit"); - if(nullptr == fTofHit) LOG(FATAL) << "No TofHit in input"; + fTofHit = (TClonesArray*) fManager->GetObject("TofHit"); + if (nullptr == fTofHit) LOG(FATAL) << "No TofHit in input"; fVertex = dynamic_cast<CbmVertex*>(fManager->GetObject("PrimaryVertex.")); - fEvent=0; - - fManager->Register("MuPlus" ,"Much",fMuPlus,kTRUE); - fManager->Register("MuMinus" ,"Much",fMuMinus,kTRUE); - + fEvent = 0; + + fManager->Register("MuPlus", "Much", fMuPlus, kTRUE); + fManager->Register("MuMinus", "Much", fMuMinus, kTRUE); + fGeoScheme = CbmMuchGeoScheme::Instance(); -// fGeoScheme->Init(fDigiFileName); -// fLastStationIndex = fGeoScheme->GetNStations()-1; -// fNLayers = 0; -// for (Int_t i=0;i<=fLastStationIndex;i++){ -// fNLayers+=fGeoScheme->GetLayerSides(i).size()/2; -// } + // fGeoScheme->Init(fDigiFileName); + // fLastStationIndex = fGeoScheme->GetNStations()-1; + // fNLayers = 0; + // for (Int_t i=0;i<=fLastStationIndex;i++){ + // fNLayers+=fGeoScheme->GetLayerSides(i).size()/2; + // } fFitter = new CbmStsKFTrackFitter(); fFitter->Init(); - + fFitterGlobal = new CbmGlobalTrackFitterKF(); fFitterGlobal->Init(); - + fFitterTRD = new CbmTrdTrackFitterKF(); fFitterTRD->Init(); fFitterTRD->SetPid(13); - - YPt_pluto = new TH2D("YPt_pluto", "PLUTO signal",YBINNING,PTBINNING); - - if(fPlutoFileName != ""){ - if(fPlutoFileName.Contains("root")){ - fParticles = new TClonesArray("PParticle",100); - fPlutoFile = new TFile(fPlutoFileName.Data()); - fInputTree = (TTree*)fPlutoFile->Get("data"); - fInputTree->SetBranchAddress("Particles", &fParticles); - for(int iEvent = 0;iEvent<fInputTree->GetEntries();iEvent++){ - fInputTree->GetEntry(iEvent); - Int_t NofPart = fParticles->GetEntriesFast(); - PParticle* Part1 = (PParticle*) fParticles->At(NofPart-2); - PParticle* Part2 = (PParticle*) fParticles->At(NofPart-1); - TLorentzVector mom1 = Part1->Vect4(); - TLorentzVector mom2 = Part2->Vect4(); - TLorentzVector Mom = mom1 + mom2; - YPt_pluto->Fill(Mom.Rapidity(),Mom.Pt(),1./(Double_t)fInputTree->GetEntries()); - }} + + YPt_pluto = new TH2D("YPt_pluto", "PLUTO signal", YBINNING, PTBINNING); + + if (fPlutoFileName != "") { + if (fPlutoFileName.Contains("root")) { + fParticles = new TClonesArray("PParticle", 100); + fPlutoFile = new TFile(fPlutoFileName.Data()); + fInputTree = (TTree*) fPlutoFile->Get("data"); + fInputTree->SetBranchAddress("Particles", &fParticles); + for (int iEvent = 0; iEvent < fInputTree->GetEntries(); iEvent++) { + fInputTree->GetEntry(iEvent); + Int_t NofPart = fParticles->GetEntriesFast(); + PParticle* Part1 = (PParticle*) fParticles->At(NofPart - 2); + PParticle* Part2 = (PParticle*) fParticles->At(NofPart - 1); + TLorentzVector mom1 = Part1->Vect4(); + TLorentzVector mom2 = Part2->Vect4(); + TLorentzVector Mom = mom1 + mom2; + YPt_pluto->Fill( + Mom.Rapidity(), Mom.Pt(), 1. / (Double_t) fInputTree->GetEntries()); + } + } } - + TString title1 = "STS accepted MC signal"; TString title2 = "STS+MUCH accepted MC signal"; TString title3 = "STS+MUCH+TRD accepted MC signal"; TString title4 = "STS+MUCH+TRD+TOF accepted MC signal"; - - YPt_StsAcc = new TH2D("YPt_StsAcc", title1,YBINNING,PTBINNING); + + YPt_StsAcc = new TH2D("YPt_StsAcc", title1, YBINNING, PTBINNING); YPt_StsAcc->GetXaxis()->SetTitle("Y"); YPt_StsAcc->GetYaxis()->SetTitle("P_{t} (GeV/c)"); - - YPt_StsMuchAcc = (TH2D*)YPt_StsAcc->Clone("YPt_StsMuchAcc"); YPt_StsMuchAcc->SetTitle(title2); - YPt_StsMuchTrdAcc = (TH2D*)YPt_StsAcc->Clone("YPt_StsMuchTrdAcc"); YPt_StsMuchTrdAcc->SetTitle(title3); - YPt_StsMuchTrdTofAcc= (TH2D*)YPt_StsAcc->Clone("YPt_StsMuchTrdTofAcc"); YPt_StsMuchTrdTofAcc->SetTitle(title4); - - YPtM = new TH3D("YPtM", "Reconstrcuted YPtM spectrum",YBINNING,PTBINNING,MBINNING); - - acc_P[0][0] = new TProfile("signal_accP_Sts",title1,PBINNING); + + YPt_StsMuchAcc = (TH2D*) YPt_StsAcc->Clone("YPt_StsMuchAcc"); + YPt_StsMuchAcc->SetTitle(title2); + YPt_StsMuchTrdAcc = (TH2D*) YPt_StsAcc->Clone("YPt_StsMuchTrdAcc"); + YPt_StsMuchTrdAcc->SetTitle(title3); + YPt_StsMuchTrdTofAcc = (TH2D*) YPt_StsAcc->Clone("YPt_StsMuchTrdTofAcc"); + YPt_StsMuchTrdTofAcc->SetTitle(title4); + + YPtM = new TH3D( + "YPtM", "Reconstrcuted YPtM spectrum", YBINNING, PTBINNING, MBINNING); + + acc_P[0][0] = new TProfile("signal_accP_Sts", title1, PBINNING); acc_P[0][0]->GetXaxis()->SetTitle("P (GeV/c)"); acc_P[0][0]->GetYaxis()->SetTitle("%"); - - acc_P[1][0] = (TProfile*)acc_P[0][0]->Clone("signal_accP_StsMuch"); acc_P[1][0]->SetTitle(title2); - acc_P[2][0] = (TProfile*)acc_P[0][0]->Clone("signal_accP_StsMuchTrd"); acc_P[2][0]->SetTitle(title3); - acc_P[3][0] = (TProfile*)acc_P[0][0]->Clone("signal_accP_StsMuchTrdTof");acc_P[3][0]->SetTitle(title4); - - acc_Theta[0][0] = new TProfile("signal_accTheta_Sts",title1,THETABINNING); + + acc_P[1][0] = (TProfile*) acc_P[0][0]->Clone("signal_accP_StsMuch"); + acc_P[1][0]->SetTitle(title2); + acc_P[2][0] = (TProfile*) acc_P[0][0]->Clone("signal_accP_StsMuchTrd"); + acc_P[2][0]->SetTitle(title3); + acc_P[3][0] = (TProfile*) acc_P[0][0]->Clone("signal_accP_StsMuchTrdTof"); + acc_P[3][0]->SetTitle(title4); + + acc_Theta[0][0] = new TProfile("signal_accTheta_Sts", title1, THETABINNING); acc_Theta[0][0]->GetXaxis()->SetTitle("#Theta (#circ)"); acc_Theta[0][0]->GetYaxis()->SetTitle("%"); - - acc_Theta[1][0] = (TProfile*)acc_Theta[0][0]->Clone("signal_accTheta_StsMuch"); acc_Theta[1][0]->SetTitle(title2); - acc_Theta[2][0] = (TProfile*)acc_Theta[0][0]->Clone("signal_accTheta_StsMuchTrd"); acc_Theta[2][0]->SetTitle(title3); - acc_Theta[3][0] = (TProfile*)acc_Theta[0][0]->Clone("signal_accTheta_StsMuchTrdTof");acc_Theta[3][0]->SetTitle(title4); - + + acc_Theta[1][0] = + (TProfile*) acc_Theta[0][0]->Clone("signal_accTheta_StsMuch"); + acc_Theta[1][0]->SetTitle(title2); + acc_Theta[2][0] = + (TProfile*) acc_Theta[0][0]->Clone("signal_accTheta_StsMuchTrd"); + acc_Theta[2][0]->SetTitle(title3); + acc_Theta[3][0] = + (TProfile*) acc_Theta[0][0]->Clone("signal_accTheta_StsMuchTrdTof"); + acc_Theta[3][0]->SetTitle(title4); + title1 = "STS accepted MC #mu+"; title2 = "STS+MUCH accepted MC #mu+"; title3 = "STS+MUCH+TRD accepted MC #mu+"; title4 = "STS+MUCH+TRD+TOF accepted MC #mu+"; - - acc_P[0][1] = (TProfile*)acc_P[0][0]->Clone("muPl_accP_Sts"); acc_P[0][1]->SetTitle(title1); - acc_P[1][1] = (TProfile*)acc_P[0][0]->Clone("muPl_accP_StsMuch"); acc_P[1][1]->SetTitle(title2); - acc_P[2][1] = (TProfile*)acc_P[0][0]->Clone("muPl_accP_StsMuchTrd"); acc_P[2][1]->SetTitle(title3); - acc_P[3][1] = (TProfile*)acc_P[0][0]->Clone("muPl_accP_StsMuchTrdTof");acc_P[3][1]->SetTitle(title4); - - acc_Theta[0][1] = (TProfile*)acc_Theta[0][0]->Clone("muPl_accTheta_Sts"); acc_Theta[0][1]->SetTitle(title1); - acc_Theta[1][1] = (TProfile*)acc_Theta[0][0]->Clone("muPl_accTheta_StsMuch"); acc_Theta[1][1]->SetTitle(title2); - acc_Theta[2][1] = (TProfile*)acc_Theta[0][0]->Clone("muPl_accTheta_StsMuchTrd"); acc_Theta[2][1]->SetTitle(title3); - acc_Theta[3][1] = (TProfile*)acc_Theta[0][0]->Clone("muPl_accTheta_StsMuchTrdTof");acc_Theta[3][1]->SetTitle(title4); - + + acc_P[0][1] = (TProfile*) acc_P[0][0]->Clone("muPl_accP_Sts"); + acc_P[0][1]->SetTitle(title1); + acc_P[1][1] = (TProfile*) acc_P[0][0]->Clone("muPl_accP_StsMuch"); + acc_P[1][1]->SetTitle(title2); + acc_P[2][1] = (TProfile*) acc_P[0][0]->Clone("muPl_accP_StsMuchTrd"); + acc_P[2][1]->SetTitle(title3); + acc_P[3][1] = (TProfile*) acc_P[0][0]->Clone("muPl_accP_StsMuchTrdTof"); + acc_P[3][1]->SetTitle(title4); + + acc_Theta[0][1] = (TProfile*) acc_Theta[0][0]->Clone("muPl_accTheta_Sts"); + acc_Theta[0][1]->SetTitle(title1); + acc_Theta[1][1] = (TProfile*) acc_Theta[0][0]->Clone("muPl_accTheta_StsMuch"); + acc_Theta[1][1]->SetTitle(title2); + acc_Theta[2][1] = + (TProfile*) acc_Theta[0][0]->Clone("muPl_accTheta_StsMuchTrd"); + acc_Theta[2][1]->SetTitle(title3); + acc_Theta[3][1] = + (TProfile*) acc_Theta[0][0]->Clone("muPl_accTheta_StsMuchTrdTof"); + acc_Theta[3][1]->SetTitle(title4); + title1 = "STS accepted MC #mu-"; title2 = "STS+MUCH accepted MC #mu-"; title3 = "STS+MUCH+TRD accepted MC #mu-"; title4 = "STS+MUCH+TRD+TOF accepted MC #mu-"; - - acc_P[0][2] = (TProfile*)acc_P[0][0]->Clone("muMn_accP_Sts"); acc_P[0][2]->SetTitle(title1); - acc_P[1][2] = (TProfile*)acc_P[0][0]->Clone("muMn_accP_StsMuch"); acc_P[1][2]->SetTitle(title2); - acc_P[2][2] = (TProfile*)acc_P[0][0]->Clone("muMn_accP_StsMuchTrd"); acc_P[2][2]->SetTitle(title3); - acc_P[3][2] = (TProfile*)acc_P[0][0]->Clone("muMn_accP_StsMuchTrdTof");acc_P[3][2]->SetTitle(title4); - - acc_Theta[0][2] = (TProfile*)acc_Theta[0][0]->Clone("muMn_accTheta_Sts"); acc_Theta[0][2]->SetTitle(title1); - acc_Theta[1][2] = (TProfile*)acc_Theta[0][0]->Clone("muMn_accTheta_StsMuch"); acc_Theta[1][2]->SetTitle(title2); - acc_Theta[2][2] = (TProfile*)acc_Theta[0][0]->Clone("muMn_accTheta_StsMuchTrd"); acc_Theta[2][2]->SetTitle(title3); - acc_Theta[3][2] = (TProfile*)acc_Theta[0][0]->Clone("muMn_accTheta_StsMuchTrdTof");acc_Theta[3][2]->SetTitle(title4); - - TString title0 = "reconstructed MC signal after primary vertex cut (PVC)"; // PVC: chi2 of STS track in target - title1 = "reconstructed MC signal after PVC and STS cuts (StsCs)"; // StsCs: chi2 of STS track and number of STS hits - title2 = "reconstructed MC signal after PVC+StsCs and MUCH cuts (MuchCs)"; // MuchCs: chi2 of MUCH track and number of MUCH hits - title3 = "reconstructed MC signal after PVC+StsCs+MuchCs and TRD cut (TrdC)"; // TrdC: number of TRD hits - title4 = "reconstructed MC signal after PVC+StsCs+MuchCs+TrdC and TOF cut"; // TOF: cut using mass distribution from time measurement - - effReco_P[0][0] = (TProfile*)acc_P[0][0]->Clone("signal_effRecoP_VtxSts"); effReco_P[0][0]->SetTitle(title1); - effReco_P[1][0] = (TProfile*)acc_P[0][0]->Clone("signal_effRecoP_VtxStsMuch"); effReco_P[1][0]->SetTitle(title2); - effReco_P[2][0] = (TProfile*)acc_P[0][0]->Clone("signal_effRecoP_VtxStsMuchTrd"); effReco_P[2][0]->SetTitle(title3); - effReco_P[3][0] = (TProfile*)acc_P[0][0]->Clone("signal_effRecoP_VtxStsMuchTrdTof");effReco_P[3][0]->SetTitle(title4); - - effReco_Theta[0][0] = (TProfile*)acc_Theta[0][0]->Clone("signal_effRecoTheta_VtxSts"); effReco_Theta[0][0]->SetTitle(title1); - effReco_Theta[1][0] = (TProfile*)acc_Theta[0][0]->Clone("signal_effRecoTheta_VtxStsMuch"); effReco_Theta[1][0]->SetTitle(title2); - effReco_Theta[2][0] = (TProfile*)acc_Theta[0][0]->Clone("signal_effRecoTheta_VtxStsMuchTrd"); effReco_Theta[2][0]->SetTitle(title3); - effReco_Theta[3][0] = (TProfile*)acc_Theta[0][0]->Clone("signal_effRecoTheta_VtxStsMuchTrdTof");effReco_Theta[3][0]->SetTitle(title4); - - eff4pi_P[0][0] = (TProfile*)acc_P[0][0]->Clone("signal_eff4piP_Vtx"); eff4pi_P[0][0]->SetTitle(title0); - eff4pi_P[1][0] = (TProfile*)acc_P[0][0]->Clone("signal_eff4piP_VtxSts"); eff4pi_P[1][0]->SetTitle(title1); - eff4pi_P[2][0] = (TProfile*)acc_P[0][0]->Clone("signal_eff4piP_VtxStsMuch"); eff4pi_P[2][0]->SetTitle(title2); - eff4pi_P[3][0] = (TProfile*)acc_P[0][0]->Clone("signal_eff4piP_VtxStsMuchTrd"); eff4pi_P[3][0]->SetTitle(title3); - eff4pi_P[4][0] = (TProfile*)acc_P[0][0]->Clone("signal_eff4piP_VtxStsMuchTrdTof");eff4pi_P[4][0]->SetTitle(title4); - - eff4pi_Theta[0][0] = (TProfile*)acc_Theta[0][0]->Clone("signal_eff4piTheta_Vtx"); eff4pi_Theta[0][0]->SetTitle(title0); - eff4pi_Theta[1][0] = (TProfile*)acc_Theta[0][0]->Clone("signal_eff4piTheta_VtxSts"); eff4pi_Theta[1][0]->SetTitle(title1); - eff4pi_Theta[2][0] = (TProfile*)acc_Theta[0][0]->Clone("signal_eff4piTheta_VtxStsMuch"); eff4pi_Theta[2][0]->SetTitle(title2); - eff4pi_Theta[3][0] = (TProfile*)acc_Theta[0][0]->Clone("signal_eff4piTheta_VtxStsMuchTrd"); eff4pi_Theta[3][0]->SetTitle(title3); - eff4pi_Theta[4][0] = (TProfile*)acc_Theta[0][0]->Clone("signal_eff4piTheta_VtxStsMuchTrdTof");eff4pi_Theta[4][0]->SetTitle(title4); - - YPt_VtxReco = (TH2D*)YPt_StsAcc->Clone("YPt_VtxReco"); YPt_VtxReco->SetTitle(title0); - YPt_VtxStsReco = (TH2D*)YPt_StsAcc->Clone("YPt_VtxStsReco"); YPt_VtxStsReco->SetTitle(title1); - YPt_VtxStsMuchReco = (TH2D*)YPt_StsAcc->Clone("YPt_VtxStsMuchReco"); YPt_VtxStsMuchReco->SetTitle(title2); - YPt_VtxStsMuchTrdReco = (TH2D*)YPt_StsAcc->Clone("YPt_VtxStsMuchTrdReco"); YPt_VtxStsMuchTrdReco->SetTitle(title3); - YPt_VtxStsMuchTrdTofReco= (TH2D*)YPt_StsAcc->Clone("YPt_VtxStsMuchTrdTofReco"); YPt_VtxStsMuchTrdTofReco->SetTitle(title4); - - title0 = "reconstructed MC #mu+ after primary vertex cut (PVC)"; - title1 = "reconstructed MC #mu+ after PVC and STS cuts (StsCs)"; - title2 = "reconstructed MC #mu+ after PVC+StsCs and MUCH cuts (MuchCs)"; - title3 = "reconstructed MC #mu+ after PVC+StsCs+MuchCs and TRD cut (TrdC)"; - title4 = "reconstructed MC #mu+ after PVC+StsCs+MuchCs+TrdC and TOF cut"; - - effReco_P[0][1] = (TProfile*)acc_P[0][0]->Clone("muPl_effRecoP_VtxSts"); effReco_P[0][1]->SetTitle(title1); - effReco_P[1][1] = (TProfile*)acc_P[0][0]->Clone("muPl_effRecoP_VtxStsMuch"); effReco_P[1][1]->SetTitle(title2); - effReco_P[2][1] = (TProfile*)acc_P[0][0]->Clone("muPl_effRecoP_VtxStsMuchTrd"); effReco_P[2][1]->SetTitle(title3); - effReco_P[3][1] = (TProfile*)acc_P[0][0]->Clone("muPl_effRecoP_VtxStsMuchTrdTof");effReco_P[3][1]->SetTitle(title4); - - effReco_Theta[0][1] = (TProfile*)acc_Theta[0][0]->Clone("muPl_effRecoTheta_VtxSts"); effReco_Theta[0][1]->SetTitle(title1); - effReco_Theta[1][1] = (TProfile*)acc_Theta[0][0]->Clone("muPl_effRecoTheta_VtxStsMuch"); effReco_Theta[1][1]->SetTitle(title2); - effReco_Theta[2][1] = (TProfile*)acc_Theta[0][0]->Clone("muPl_effRecoTheta_VtxStsMuchTrd"); effReco_Theta[2][1]->SetTitle(title3); - effReco_Theta[3][1] = (TProfile*)acc_Theta[0][0]->Clone("muPl_effRecoTheta_VtxStsMuchTrdTof");effReco_Theta[3][1]->SetTitle(title4); - - eff4pi_P[0][1] = (TProfile*)acc_P[0][0]->Clone("muPl_eff4piP_Vtx"); eff4pi_P[0][1]->SetTitle(title0); - eff4pi_P[1][1] = (TProfile*)acc_P[0][0]->Clone("muPl_eff4piP_VtxSts"); eff4pi_P[1][1]->SetTitle(title1); - eff4pi_P[2][1] = (TProfile*)acc_P[0][0]->Clone("muPl_eff4piP_VtxStsMuch"); eff4pi_P[2][1]->SetTitle(title2); - eff4pi_P[3][1] = (TProfile*)acc_P[0][0]->Clone("muPl_eff4piP_VtxStsMuchTrd"); eff4pi_P[3][1]->SetTitle(title3); - eff4pi_P[4][1] = (TProfile*)acc_P[0][0]->Clone("muPl_eff4piP_VtxStsMuchTrdTof");eff4pi_P[4][1]->SetTitle(title4); - - eff4pi_Theta[0][1] = (TProfile*)acc_Theta[0][0]->Clone("muPl_eff4piTheta_Vtx"); eff4pi_Theta[0][1]->SetTitle(title0); - eff4pi_Theta[1][1] = (TProfile*)acc_Theta[0][0]->Clone("muPl_eff4piTheta_VtxSts"); eff4pi_Theta[1][1]->SetTitle(title1); - eff4pi_Theta[2][1] = (TProfile*)acc_Theta[0][0]->Clone("muPl_eff4piTheta_VtxStsMuch"); eff4pi_Theta[2][1]->SetTitle(title2); - eff4pi_Theta[3][1] = (TProfile*)acc_Theta[0][0]->Clone("muPl_eff4piTheta_VtxStsMuchTrd"); eff4pi_Theta[3][1]->SetTitle(title3); - eff4pi_Theta[4][1] = (TProfile*)acc_Theta[0][0]->Clone("muPl_eff4piTheta_VtxStsMuchTrdTof");eff4pi_Theta[4][1]->SetTitle(title4); - - title0 = "reconstructed MC #mu- after primary vertex cut (PVC)"; - title1 = "reconstructed MC #mu- after PVC and STS cuts (StsCs)"; - title2 = "reconstructed MC #mu- after PVC+StsCs and MUCH cuts (MuchCs)"; - title3 = "reconstructed MC #mu- after PVC+StsCs+MuchCs and TRD cut (TrdC)"; - title4 = "reconstructed MC #mu- after PVC+StsCs+MuchCs+TrdC and TOF cut"; - - effReco_P[0][2] = (TProfile*)acc_P[0][0]->Clone("muMn_effRecoP_VtxSts"); effReco_P[0][2]->SetTitle(title1); - effReco_P[1][2] = (TProfile*)acc_P[0][0]->Clone("muMn_effRecoP_VtxStsMuch"); effReco_P[1][2]->SetTitle(title2); - effReco_P[2][2] = (TProfile*)acc_P[0][0]->Clone("muMn_effRecoP_VtxStsMuchTrd"); effReco_P[2][2]->SetTitle(title3); - effReco_P[3][2] = (TProfile*)acc_P[0][0]->Clone("muMn_effRecoP_VtxStsMuchTrdTof");effReco_P[3][2]->SetTitle(title4); - - effReco_Theta[0][2] = (TProfile*)acc_Theta[0][0]->Clone("muMn_effRecoTheta_VtxSts"); effReco_Theta[0][2]->SetTitle(title1); - effReco_Theta[1][2] = (TProfile*)acc_Theta[0][0]->Clone("muMn_effRecoTheta_VtxStsMuch"); effReco_Theta[1][2]->SetTitle(title2); - effReco_Theta[2][2] = (TProfile*)acc_Theta[0][0]->Clone("muMn_effRecoTheta_VtxStsMuchTrd"); effReco_Theta[2][2]->SetTitle(title3); - effReco_Theta[3][2] = (TProfile*)acc_Theta[0][0]->Clone("muMn_effRecoTheta_VtxStsMuchTrdTof");effReco_Theta[3][2]->SetTitle(title4); - - eff4pi_P[0][2] = (TProfile*)acc_P[0][0]->Clone("muMn_eff4piP_Vtx"); eff4pi_P[0][2]->SetTitle(title0); - eff4pi_P[1][2] = (TProfile*)acc_P[0][0]->Clone("muMn_eff4piP_VtxSts"); eff4pi_P[1][2]->SetTitle(title1); - eff4pi_P[2][2] = (TProfile*)acc_P[0][0]->Clone("muMn_eff4piP_VtxStsMuch"); eff4pi_P[2][2]->SetTitle(title2); - eff4pi_P[3][2] = (TProfile*)acc_P[0][0]->Clone("muMn_eff4piP_VtxStsMuchTrd"); eff4pi_P[3][2]->SetTitle(title3); - eff4pi_P[4][2] = (TProfile*)acc_P[0][0]->Clone("muMn_eff4piP_VtxStsMuchTrdTof");eff4pi_P[4][2]->SetTitle(title4); - - eff4pi_Theta[0][2] = (TProfile*)acc_Theta[0][0]->Clone("muMn_eff4piTheta_Vtx"); eff4pi_Theta[0][2]->SetTitle(title0); - eff4pi_Theta[1][2] = (TProfile*)acc_Theta[0][0]->Clone("muMn_eff4piTheta_VtxSts"); eff4pi_Theta[1][2]->SetTitle(title1); - eff4pi_Theta[2][2] = (TProfile*)acc_Theta[0][0]->Clone("muMn_eff4piTheta_VtxStsMuch"); eff4pi_Theta[2][2]->SetTitle(title2); - eff4pi_Theta[3][2] = (TProfile*)acc_Theta[0][0]->Clone("muMn_eff4piTheta_VtxStsMuchTrd"); eff4pi_Theta[3][2]->SetTitle(title3); - eff4pi_Theta[4][2] = (TProfile*)acc_Theta[0][0]->Clone("muMn_eff4piTheta_VtxStsMuchTrdTof");eff4pi_Theta[4][2]->SetTitle(title4); - - BgSup[0] = new TH1D("h0","all STS tracks",PBINNING); + + acc_P[0][2] = (TProfile*) acc_P[0][0]->Clone("muMn_accP_Sts"); + acc_P[0][2]->SetTitle(title1); + acc_P[1][2] = (TProfile*) acc_P[0][0]->Clone("muMn_accP_StsMuch"); + acc_P[1][2]->SetTitle(title2); + acc_P[2][2] = (TProfile*) acc_P[0][0]->Clone("muMn_accP_StsMuchTrd"); + acc_P[2][2]->SetTitle(title3); + acc_P[3][2] = (TProfile*) acc_P[0][0]->Clone("muMn_accP_StsMuchTrdTof"); + acc_P[3][2]->SetTitle(title4); + + acc_Theta[0][2] = (TProfile*) acc_Theta[0][0]->Clone("muMn_accTheta_Sts"); + acc_Theta[0][2]->SetTitle(title1); + acc_Theta[1][2] = (TProfile*) acc_Theta[0][0]->Clone("muMn_accTheta_StsMuch"); + acc_Theta[1][2]->SetTitle(title2); + acc_Theta[2][2] = + (TProfile*) acc_Theta[0][0]->Clone("muMn_accTheta_StsMuchTrd"); + acc_Theta[2][2]->SetTitle(title3); + acc_Theta[3][2] = + (TProfile*) acc_Theta[0][0]->Clone("muMn_accTheta_StsMuchTrdTof"); + acc_Theta[3][2]->SetTitle(title4); + + TString title0 = + "reconstructed MC signal after primary vertex cut (PVC)"; // PVC: chi2 of STS track in target + title1 = + "reconstructed MC signal after PVC and STS cuts (StsCs)"; // StsCs: chi2 of STS track and number of STS hits + title2 = + "reconstructed MC signal after PVC+StsCs and MUCH cuts (MuchCs)"; // MuchCs: chi2 of MUCH track and number of MUCH hits + title3 = + "reconstructed MC signal after PVC+StsCs+MuchCs and TRD cut (TrdC)"; // TrdC: number of TRD hits + title4 = + "reconstructed MC signal after PVC+StsCs+MuchCs+TrdC and TOF cut"; // TOF: cut using mass distribution from time measurement + + effReco_P[0][0] = (TProfile*) acc_P[0][0]->Clone("signal_effRecoP_VtxSts"); + effReco_P[0][0]->SetTitle(title1); + effReco_P[1][0] = + (TProfile*) acc_P[0][0]->Clone("signal_effRecoP_VtxStsMuch"); + effReco_P[1][0]->SetTitle(title2); + effReco_P[2][0] = + (TProfile*) acc_P[0][0]->Clone("signal_effRecoP_VtxStsMuchTrd"); + effReco_P[2][0]->SetTitle(title3); + effReco_P[3][0] = + (TProfile*) acc_P[0][0]->Clone("signal_effRecoP_VtxStsMuchTrdTof"); + effReco_P[3][0]->SetTitle(title4); + + effReco_Theta[0][0] = + (TProfile*) acc_Theta[0][0]->Clone("signal_effRecoTheta_VtxSts"); + effReco_Theta[0][0]->SetTitle(title1); + effReco_Theta[1][0] = + (TProfile*) acc_Theta[0][0]->Clone("signal_effRecoTheta_VtxStsMuch"); + effReco_Theta[1][0]->SetTitle(title2); + effReco_Theta[2][0] = + (TProfile*) acc_Theta[0][0]->Clone("signal_effRecoTheta_VtxStsMuchTrd"); + effReco_Theta[2][0]->SetTitle(title3); + effReco_Theta[3][0] = + (TProfile*) acc_Theta[0][0]->Clone("signal_effRecoTheta_VtxStsMuchTrdTof"); + effReco_Theta[3][0]->SetTitle(title4); + + eff4pi_P[0][0] = (TProfile*) acc_P[0][0]->Clone("signal_eff4piP_Vtx"); + eff4pi_P[0][0]->SetTitle(title0); + eff4pi_P[1][0] = (TProfile*) acc_P[0][0]->Clone("signal_eff4piP_VtxSts"); + eff4pi_P[1][0]->SetTitle(title1); + eff4pi_P[2][0] = (TProfile*) acc_P[0][0]->Clone("signal_eff4piP_VtxStsMuch"); + eff4pi_P[2][0]->SetTitle(title2); + eff4pi_P[3][0] = + (TProfile*) acc_P[0][0]->Clone("signal_eff4piP_VtxStsMuchTrd"); + eff4pi_P[3][0]->SetTitle(title3); + eff4pi_P[4][0] = + (TProfile*) acc_P[0][0]->Clone("signal_eff4piP_VtxStsMuchTrdTof"); + eff4pi_P[4][0]->SetTitle(title4); + + eff4pi_Theta[0][0] = + (TProfile*) acc_Theta[0][0]->Clone("signal_eff4piTheta_Vtx"); + eff4pi_Theta[0][0]->SetTitle(title0); + eff4pi_Theta[1][0] = + (TProfile*) acc_Theta[0][0]->Clone("signal_eff4piTheta_VtxSts"); + eff4pi_Theta[1][0]->SetTitle(title1); + eff4pi_Theta[2][0] = + (TProfile*) acc_Theta[0][0]->Clone("signal_eff4piTheta_VtxStsMuch"); + eff4pi_Theta[2][0]->SetTitle(title2); + eff4pi_Theta[3][0] = + (TProfile*) acc_Theta[0][0]->Clone("signal_eff4piTheta_VtxStsMuchTrd"); + eff4pi_Theta[3][0]->SetTitle(title3); + eff4pi_Theta[4][0] = + (TProfile*) acc_Theta[0][0]->Clone("signal_eff4piTheta_VtxStsMuchTrdTof"); + eff4pi_Theta[4][0]->SetTitle(title4); + + YPt_VtxReco = (TH2D*) YPt_StsAcc->Clone("YPt_VtxReco"); + YPt_VtxReco->SetTitle(title0); + YPt_VtxStsReco = (TH2D*) YPt_StsAcc->Clone("YPt_VtxStsReco"); + YPt_VtxStsReco->SetTitle(title1); + YPt_VtxStsMuchReco = (TH2D*) YPt_StsAcc->Clone("YPt_VtxStsMuchReco"); + YPt_VtxStsMuchReco->SetTitle(title2); + YPt_VtxStsMuchTrdReco = (TH2D*) YPt_StsAcc->Clone("YPt_VtxStsMuchTrdReco"); + YPt_VtxStsMuchTrdReco->SetTitle(title3); + YPt_VtxStsMuchTrdTofReco = + (TH2D*) YPt_StsAcc->Clone("YPt_VtxStsMuchTrdTofReco"); + YPt_VtxStsMuchTrdTofReco->SetTitle(title4); + + title0 = "reconstructed MC #mu+ after primary vertex cut (PVC)"; + title1 = "reconstructed MC #mu+ after PVC and STS cuts (StsCs)"; + title2 = "reconstructed MC #mu+ after PVC+StsCs and MUCH cuts (MuchCs)"; + title3 = "reconstructed MC #mu+ after PVC+StsCs+MuchCs and TRD cut (TrdC)"; + title4 = "reconstructed MC #mu+ after PVC+StsCs+MuchCs+TrdC and TOF cut"; + + effReco_P[0][1] = (TProfile*) acc_P[0][0]->Clone("muPl_effRecoP_VtxSts"); + effReco_P[0][1]->SetTitle(title1); + effReco_P[1][1] = (TProfile*) acc_P[0][0]->Clone("muPl_effRecoP_VtxStsMuch"); + effReco_P[1][1]->SetTitle(title2); + effReco_P[2][1] = + (TProfile*) acc_P[0][0]->Clone("muPl_effRecoP_VtxStsMuchTrd"); + effReco_P[2][1]->SetTitle(title3); + effReco_P[3][1] = + (TProfile*) acc_P[0][0]->Clone("muPl_effRecoP_VtxStsMuchTrdTof"); + effReco_P[3][1]->SetTitle(title4); + + effReco_Theta[0][1] = + (TProfile*) acc_Theta[0][0]->Clone("muPl_effRecoTheta_VtxSts"); + effReco_Theta[0][1]->SetTitle(title1); + effReco_Theta[1][1] = + (TProfile*) acc_Theta[0][0]->Clone("muPl_effRecoTheta_VtxStsMuch"); + effReco_Theta[1][1]->SetTitle(title2); + effReco_Theta[2][1] = + (TProfile*) acc_Theta[0][0]->Clone("muPl_effRecoTheta_VtxStsMuchTrd"); + effReco_Theta[2][1]->SetTitle(title3); + effReco_Theta[3][1] = + (TProfile*) acc_Theta[0][0]->Clone("muPl_effRecoTheta_VtxStsMuchTrdTof"); + effReco_Theta[3][1]->SetTitle(title4); + + eff4pi_P[0][1] = (TProfile*) acc_P[0][0]->Clone("muPl_eff4piP_Vtx"); + eff4pi_P[0][1]->SetTitle(title0); + eff4pi_P[1][1] = (TProfile*) acc_P[0][0]->Clone("muPl_eff4piP_VtxSts"); + eff4pi_P[1][1]->SetTitle(title1); + eff4pi_P[2][1] = (TProfile*) acc_P[0][0]->Clone("muPl_eff4piP_VtxStsMuch"); + eff4pi_P[2][1]->SetTitle(title2); + eff4pi_P[3][1] = (TProfile*) acc_P[0][0]->Clone("muPl_eff4piP_VtxStsMuchTrd"); + eff4pi_P[3][1]->SetTitle(title3); + eff4pi_P[4][1] = + (TProfile*) acc_P[0][0]->Clone("muPl_eff4piP_VtxStsMuchTrdTof"); + eff4pi_P[4][1]->SetTitle(title4); + + eff4pi_Theta[0][1] = + (TProfile*) acc_Theta[0][0]->Clone("muPl_eff4piTheta_Vtx"); + eff4pi_Theta[0][1]->SetTitle(title0); + eff4pi_Theta[1][1] = + (TProfile*) acc_Theta[0][0]->Clone("muPl_eff4piTheta_VtxSts"); + eff4pi_Theta[1][1]->SetTitle(title1); + eff4pi_Theta[2][1] = + (TProfile*) acc_Theta[0][0]->Clone("muPl_eff4piTheta_VtxStsMuch"); + eff4pi_Theta[2][1]->SetTitle(title2); + eff4pi_Theta[3][1] = + (TProfile*) acc_Theta[0][0]->Clone("muPl_eff4piTheta_VtxStsMuchTrd"); + eff4pi_Theta[3][1]->SetTitle(title3); + eff4pi_Theta[4][1] = + (TProfile*) acc_Theta[0][0]->Clone("muPl_eff4piTheta_VtxStsMuchTrdTof"); + eff4pi_Theta[4][1]->SetTitle(title4); + + title0 = "reconstructed MC #mu- after primary vertex cut (PVC)"; + title1 = "reconstructed MC #mu- after PVC and STS cuts (StsCs)"; + title2 = "reconstructed MC #mu- after PVC+StsCs and MUCH cuts (MuchCs)"; + title3 = "reconstructed MC #mu- after PVC+StsCs+MuchCs and TRD cut (TrdC)"; + title4 = "reconstructed MC #mu- after PVC+StsCs+MuchCs+TrdC and TOF cut"; + + effReco_P[0][2] = (TProfile*) acc_P[0][0]->Clone("muMn_effRecoP_VtxSts"); + effReco_P[0][2]->SetTitle(title1); + effReco_P[1][2] = (TProfile*) acc_P[0][0]->Clone("muMn_effRecoP_VtxStsMuch"); + effReco_P[1][2]->SetTitle(title2); + effReco_P[2][2] = + (TProfile*) acc_P[0][0]->Clone("muMn_effRecoP_VtxStsMuchTrd"); + effReco_P[2][2]->SetTitle(title3); + effReco_P[3][2] = + (TProfile*) acc_P[0][0]->Clone("muMn_effRecoP_VtxStsMuchTrdTof"); + effReco_P[3][2]->SetTitle(title4); + + effReco_Theta[0][2] = + (TProfile*) acc_Theta[0][0]->Clone("muMn_effRecoTheta_VtxSts"); + effReco_Theta[0][2]->SetTitle(title1); + effReco_Theta[1][2] = + (TProfile*) acc_Theta[0][0]->Clone("muMn_effRecoTheta_VtxStsMuch"); + effReco_Theta[1][2]->SetTitle(title2); + effReco_Theta[2][2] = + (TProfile*) acc_Theta[0][0]->Clone("muMn_effRecoTheta_VtxStsMuchTrd"); + effReco_Theta[2][2]->SetTitle(title3); + effReco_Theta[3][2] = + (TProfile*) acc_Theta[0][0]->Clone("muMn_effRecoTheta_VtxStsMuchTrdTof"); + effReco_Theta[3][2]->SetTitle(title4); + + eff4pi_P[0][2] = (TProfile*) acc_P[0][0]->Clone("muMn_eff4piP_Vtx"); + eff4pi_P[0][2]->SetTitle(title0); + eff4pi_P[1][2] = (TProfile*) acc_P[0][0]->Clone("muMn_eff4piP_VtxSts"); + eff4pi_P[1][2]->SetTitle(title1); + eff4pi_P[2][2] = (TProfile*) acc_P[0][0]->Clone("muMn_eff4piP_VtxStsMuch"); + eff4pi_P[2][2]->SetTitle(title2); + eff4pi_P[3][2] = (TProfile*) acc_P[0][0]->Clone("muMn_eff4piP_VtxStsMuchTrd"); + eff4pi_P[3][2]->SetTitle(title3); + eff4pi_P[4][2] = + (TProfile*) acc_P[0][0]->Clone("muMn_eff4piP_VtxStsMuchTrdTof"); + eff4pi_P[4][2]->SetTitle(title4); + + eff4pi_Theta[0][2] = + (TProfile*) acc_Theta[0][0]->Clone("muMn_eff4piTheta_Vtx"); + eff4pi_Theta[0][2]->SetTitle(title0); + eff4pi_Theta[1][2] = + (TProfile*) acc_Theta[0][0]->Clone("muMn_eff4piTheta_VtxSts"); + eff4pi_Theta[1][2]->SetTitle(title1); + eff4pi_Theta[2][2] = + (TProfile*) acc_Theta[0][0]->Clone("muMn_eff4piTheta_VtxStsMuch"); + eff4pi_Theta[2][2]->SetTitle(title2); + eff4pi_Theta[3][2] = + (TProfile*) acc_Theta[0][0]->Clone("muMn_eff4piTheta_VtxStsMuchTrd"); + eff4pi_Theta[3][2]->SetTitle(title3); + eff4pi_Theta[4][2] = + (TProfile*) acc_Theta[0][0]->Clone("muMn_eff4piTheta_VtxStsMuchTrdTof"); + eff4pi_Theta[4][2]->SetTitle(title4); + + BgSup[0] = new TH1D("h0", "all STS tracks", PBINNING); BgSup[0]->GetXaxis()->SetTitle("P (GeV/c)"); BgSup[0]->GetYaxis()->SetTitle("suppression"); - title0 = "reconstructed tracks after primary vertex cut (PVC)"; - title1 = "reconstructed tracks after PVC and STS cuts (StsCs)"; - title2 = "reconstructed tracks after PVC+StsCs and MUCH cuts (MuchCs)"; - title3 = "reconstructed tracks after PVC+StsCs+MuchCs and TRD cut (TrdC)"; - title4 = "reconstructed tracks after PVC+StsCs+MuchCs+TrdC and TOF cut"; - - BgSup[1] = (TH1D*)BgSup[0]->Clone("h1"); BgSup[1]->SetTitle(title0); - BgSup[2] = (TH1D*)BgSup[0]->Clone("h2"); BgSup[2]->SetTitle(title1); - BgSup[3] = (TH1D*)BgSup[0]->Clone("h3"); BgSup[3]->SetTitle(title2); - BgSup[4] = (TH1D*)BgSup[0]->Clone("h4"); BgSup[4]->SetTitle(title3); - BgSup[5] = (TH1D*)BgSup[0]->Clone("h5"); BgSup[5]->SetTitle(title4); - - TString dir = getenv("VMCWORKDIR"); - TString name= dir + "/parameters/much/TOF8gev_fitParam_sigma" - + std::to_string(fSigmaTofCut) - + "." + fSetupName - + ".root"; - - TFile *FF = new TFile(name); - - TTree *MinParamMu = (TTree*)FF->Get("MinParam"); - MinParamMu->SetBranchAddress("p0", &p0min); - MinParamMu->SetBranchAddress("p1", &p1min); - MinParamMu->SetBranchAddress("p2", &p2min); - MinParamMu->GetEntry(0); - - TTree *MaxParamMu = (TTree*)FF->Get("MaxParam"); - MaxParamMu->SetBranchAddress("p0", &p0max); - MaxParamMu->SetBranchAddress("p1", &p1max); - MaxParamMu->SetBranchAddress("p2", &p2max); - MaxParamMu->GetEntry(0); - - if(fNeurons>0)fFileAnnName = dir + "/parameters/much/muid_ann_" - + std::to_string(fNeurons)+ "_" + - + fSetupName + "_weights.txt"; - else fFileAnnName = dir + "/parameters/much/muid_ann_16_sis100_muon_lmvm_weights.txt"; - - FILE *file = fopen(fFileAnnName.Data(), "r"); - char buffer [100]; - - if(fgets(buffer, 100, file) == NULL){ + title0 = "reconstructed tracks after primary vertex cut (PVC)"; + title1 = "reconstructed tracks after PVC and STS cuts (StsCs)"; + title2 = "reconstructed tracks after PVC+StsCs and MUCH cuts (MuchCs)"; + title3 = "reconstructed tracks after PVC+StsCs+MuchCs and TRD cut (TrdC)"; + title4 = "reconstructed tracks after PVC+StsCs+MuchCs+TrdC and TOF cut"; + + BgSup[1] = (TH1D*) BgSup[0]->Clone("h1"); + BgSup[1]->SetTitle(title0); + BgSup[2] = (TH1D*) BgSup[0]->Clone("h2"); + BgSup[2]->SetTitle(title1); + BgSup[3] = (TH1D*) BgSup[0]->Clone("h3"); + BgSup[3]->SetTitle(title2); + BgSup[4] = (TH1D*) BgSup[0]->Clone("h4"); + BgSup[4]->SetTitle(title3); + BgSup[5] = (TH1D*) BgSup[0]->Clone("h5"); + BgSup[5]->SetTitle(title4); + + TString dir = getenv("VMCWORKDIR"); + TString name = dir + "/parameters/much/TOF8gev_fitParam_sigma" + + std::to_string(fSigmaTofCut) + "." + fSetupName + ".root"; + + TFile* FF = new TFile(name); + + TTree* MinParamMu = (TTree*) FF->Get("MinParam"); + MinParamMu->SetBranchAddress("p0", &p0min); + MinParamMu->SetBranchAddress("p1", &p1min); + MinParamMu->SetBranchAddress("p2", &p2min); + MinParamMu->GetEntry(0); + + TTree* MaxParamMu = (TTree*) FF->Get("MaxParam"); + MaxParamMu->SetBranchAddress("p0", &p0max); + MaxParamMu->SetBranchAddress("p1", &p1max); + MaxParamMu->SetBranchAddress("p2", &p2max); + MaxParamMu->GetEntry(0); + + if (fNeurons > 0) + fFileAnnName = dir + "/parameters/much/muid_ann_" + std::to_string(fNeurons) + + "_" + +fSetupName + "_weights.txt"; + else + fFileAnnName = + dir + "/parameters/much/muid_ann_16_sis100_muon_lmvm_weights.txt"; + + FILE* file = fopen(fFileAnnName.Data(), "r"); + char buffer[100]; + + if (fgets(buffer, 100, file) == NULL) { LOG(info) << "======================================================"; LOG(info) << "The ANN weights file " << fFileAnnName << " does not exist"; - - fFileAnnName = dir + "/parameters/much/muid_ann_16_sis100_muon_lmvm_weights.txt"; - - LOG(info) << "The default ANN weights file " << fFileAnnName << " will be used"; + + fFileAnnName = + dir + "/parameters/much/muid_ann_16_sis100_muon_lmvm_weights.txt"; + + LOG(info) << "The default ANN weights file " << fFileAnnName + << " will be used"; LOG(info) << "======================================================"; fNeurons = 16; } - if(fAnnCut>0)fUseCuts=kFALSE; + if (fAnnCut > 0) fUseCuts = kFALSE; return kSUCCESS; } @@ -397,464 +544,581 @@ InitStatus CbmAnaDimuonAnalysis::Init() // ----- Public method Exec -------------------------------------------- -void CbmAnaDimuonAnalysis::Exec(Option_t* /*opt*/){ - Int_t nMCTracks = fMCTracks->GetEntriesFast(); - Int_t nStsTracks = fStsTracks->GetEntriesFast(); - Int_t nMuchTracks = fMuchTracks->GetEntriesFast(); - Int_t nGlobalTracks = fGlobalTracks->GetEntriesFast(); - +void CbmAnaDimuonAnalysis::Exec(Option_t* /*opt*/) { + Int_t nMCTracks = fMCTracks->GetEntriesFast(); + Int_t nStsTracks = fStsTracks->GetEntriesFast(); + Int_t nMuchTracks = fMuchTracks->GetEntriesFast(); + Int_t nGlobalTracks = fGlobalTracks->GetEntriesFast(); + LOG(DEBUG) << "------------------------"; LOG(DEBUG) << GetName() << ": Event " << fEvent; - LOG(DEBUG) << "Number of tracks: MC - " << nMCTracks - << ", global - " << nGlobalTracks - << ", STS - " << nStsTracks - << ", MUCH - " << nMuchTracks; + LOG(DEBUG) << "Number of tracks: MC - " << nMCTracks << ", global - " + << nGlobalTracks << ", STS - " << nStsTracks << ", MUCH - " + << nMuchTracks; LOG(DEBUG) << "------------------------"; - + TLorentzVector pMC1, pMC2, M; - - struct CbmMuon{ - Bool_t Mu; - Bool_t Nsts; - Bool_t Nmuch; - Bool_t Ntrd; - Bool_t Ntof; - Bool_t Chi2V; - Bool_t Chi2sts; - Bool_t Chi2much; - CbmMuon() : Mu(kFALSE), Nsts(kFALSE), Nmuch(kFALSE), Ntrd(kFALSE), Ntof(kFALSE), - Chi2V(kFALSE), Chi2sts(kFALSE), Chi2much(kFALSE) {;} -}; - - struct CbmMuonMC{ - Bool_t Mu; - Bool_t Nsts; - Bool_t Nmuch; - Bool_t Ntrd; - Bool_t Ntof; - CbmMuonMC() : Mu(kFALSE), Nsts(kFALSE), Nmuch(kFALSE), Ntrd(kFALSE), Ntof(kFALSE) {;} -}; -//----------------- Sort MC tracks for acceptance histograms + + struct CbmMuon { + Bool_t Mu; + Bool_t Nsts; + Bool_t Nmuch; + Bool_t Ntrd; + Bool_t Ntof; + Bool_t Chi2V; + Bool_t Chi2sts; + Bool_t Chi2much; + CbmMuon() + : Mu(kFALSE) + , Nsts(kFALSE) + , Nmuch(kFALSE) + , Ntrd(kFALSE) + , Ntof(kFALSE) + , Chi2V(kFALSE) + , Chi2sts(kFALSE) + , Chi2much(kFALSE) { + ; + } + }; + + struct CbmMuonMC { + Bool_t Mu; + Bool_t Nsts; + Bool_t Nmuch; + Bool_t Ntrd; + Bool_t Ntof; + CbmMuonMC() + : Mu(kFALSE), Nsts(kFALSE), Nmuch(kFALSE), Ntrd(kFALSE), Ntof(kFALSE) { + ; + } + }; + //----------------- Sort MC tracks for acceptance histograms CbmMuon muPl_reco; CbmMuon muMn_reco; CbmMuon signal_reco; - + CbmMuonMC muPl_mc; CbmMuonMC muMn_mc; - CbmMuonMC signal_mc; - - for(Int_t iMCTrack=0;iMCTrack<nMCTracks;iMCTrack++){ + CbmMuonMC signal_mc; + + for (Int_t iMCTrack = 0; iMCTrack < nMCTracks; iMCTrack++) { CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTracks->At(iMCTrack); - if(mcTrack->GetGeantProcessId()!=kPPrimary)continue; - if(TMath::Abs(mcTrack->GetPdgCode())!=13)continue; - if(mcTrack->GetCharge()<0){ - muMn_mc.Mu=kTRUE; - if(mcTrack->GetNPoints(ECbmModuleId::kSts)>=fNofStsCut){ - muMn_mc.Nsts=kTRUE; - if(mcTrack->GetNPoints(ECbmModuleId::kMuch)>=fNofMuchCut){ - muMn_mc.Nmuch=kTRUE; - if(mcTrack->GetNPoints(ECbmModuleId::kTrd)>=fNofTrdCut){ - muMn_mc.Ntrd=kTRUE; - if(mcTrack->GetNPoints(ECbmModuleId::kTof)>=1)muMn_mc.Ntof=kTRUE; - } - } + if (mcTrack->GetGeantProcessId() != kPPrimary) continue; + if (TMath::Abs(mcTrack->GetPdgCode()) != 13) continue; + if (mcTrack->GetCharge() < 0) { + muMn_mc.Mu = kTRUE; + if (mcTrack->GetNPoints(ECbmModuleId::kSts) >= fNofStsCut) { + muMn_mc.Nsts = kTRUE; + if (mcTrack->GetNPoints(ECbmModuleId::kMuch) >= fNofMuchCut) { + muMn_mc.Nmuch = kTRUE; + if (mcTrack->GetNPoints(ECbmModuleId::kTrd) >= fNofTrdCut) { + muMn_mc.Ntrd = kTRUE; + if (mcTrack->GetNPoints(ECbmModuleId::kTof) >= 1) + muMn_mc.Ntof = kTRUE; + } + } } - }else{ - muPl_mc.Mu=kTRUE; - if(mcTrack->GetNPoints(ECbmModuleId::kSts)>=fNofStsCut){ - muPl_mc.Nsts=kTRUE; - if(mcTrack->GetNPoints(ECbmModuleId::kMuch)>=fNofMuchCut){ - muPl_mc.Nmuch=kTRUE; - if(mcTrack->GetNPoints(ECbmModuleId::kTrd)>=fNofTrdCut){ - muPl_mc.Ntrd=kTRUE; - if(mcTrack->GetNPoints(ECbmModuleId::kTof)>=1)muPl_mc.Ntof=kTRUE; - } - } - } - } + } else { + muPl_mc.Mu = kTRUE; + if (mcTrack->GetNPoints(ECbmModuleId::kSts) >= fNofStsCut) { + muPl_mc.Nsts = kTRUE; + if (mcTrack->GetNPoints(ECbmModuleId::kMuch) >= fNofMuchCut) { + muPl_mc.Nmuch = kTRUE; + if (mcTrack->GetNPoints(ECbmModuleId::kTrd) >= fNofTrdCut) { + muPl_mc.Ntrd = kTRUE; + if (mcTrack->GetNPoints(ECbmModuleId::kTof) >= 1) + muPl_mc.Ntof = kTRUE; + } + } + } + } } - if(muPl_mc.Mu && muMn_mc.Mu) signal_mc.Mu=kTRUE; - if(muPl_mc.Nsts && muMn_mc.Nsts) signal_mc.Nsts=kTRUE; - if(muPl_mc.Nmuch && muMn_mc.Nmuch)signal_mc.Nmuch=kTRUE; - if(muPl_mc.Ntrd && muMn_mc.Ntrd) signal_mc.Ntrd=kTRUE; - if(muPl_mc.Ntof && muMn_mc.Ntof) signal_mc.Ntof=kTRUE; -//----------------- + if (muPl_mc.Mu && muMn_mc.Mu) signal_mc.Mu = kTRUE; + if (muPl_mc.Nsts && muMn_mc.Nsts) signal_mc.Nsts = kTRUE; + if (muPl_mc.Nmuch && muMn_mc.Nmuch) signal_mc.Nmuch = kTRUE; + if (muPl_mc.Ntrd && muMn_mc.Ntrd) signal_mc.Ntrd = kTRUE; + if (muPl_mc.Ntof && muMn_mc.Ntof) signal_mc.Ntof = kTRUE; + //----------------- fMuPlus->Clear(); fMuMinus->Clear(); - - Int_t iMuPlus = 0; + + Int_t iMuPlus = 0; Int_t iMuMinus = 0; - - for (Int_t iTrack=0;iTrack<nGlobalTracks;iTrack++){ - - Bool_t analysis = kFALSE; - -//----------------- Global track parameters - + + for (Int_t iTrack = 0; iTrack < nGlobalTracks; iTrack++) { + + Bool_t analysis = kFALSE; + + //----------------- Global track parameters + CbmGlobalTrack* globalTrack = (CbmGlobalTrack*) fGlobalTracks->At(iTrack); -// Double_t chi2global = globalTrack->GetChi2()/globalTrack->GetNDF(); - + // Double_t chi2global = globalTrack->GetChi2()/globalTrack->GetNDF(); + Int_t iMuchTrack = globalTrack->GetMuchTrackIndex(); Int_t iStsTrack = globalTrack->GetStsTrackIndex(); Int_t iTrdTrack = globalTrack->GetTrdTrackIndex(); Int_t iTofHit = globalTrack->GetTofHitIndex(); - - if(iStsTrack<0) continue; -//----------------- STS track parameters - + if (iStsTrack < 0) continue; + + //----------------- STS track parameters + CbmStsTrack* stsTrack = (CbmStsTrack*) fStsTracks->At(iStsTrack); - if(!stsTrack)continue; - - Int_t nStsHits = stsTrack->GetNofHits(); + if (!stsTrack) continue; + + Int_t nStsHits = stsTrack->GetNofHits(); Double_t chi2vertex = fFitter->GetChiToVertex(stsTrack); - Double_t chi2sts = 1000; - if(stsTrack->GetNDF()!=0)chi2sts=stsTrack->GetChiSq()/stsTrack->GetNDF(); - + Double_t chi2sts = 1000; + if (stsTrack->GetNDF() != 0) + chi2sts = stsTrack->GetChiSq() / stsTrack->GetNDF(); + FairTrackParam par; TLorentzVector mom; TVector3 p; - - fFitter->Extrapolate(stsTrack,fVertex->GetZ(),&par); + + fFitter->Extrapolate(stsTrack, fVertex->GetZ(), &par); par.Momentum(p); - mom.SetVectM(p,fMass); - + mom.SetVectM(p, fMass); + Double_t momentum = mom.P(); - + BgSup[0]->Fill(momentum); - - Int_t q = par.GetQp()>0 ? 1 : -1; - -//----------------- MUCH track parameters - - Int_t nMuchHits = 0; + + Int_t q = par.GetQp() > 0 ? 1 : -1; + + //----------------- MUCH track parameters + + Int_t nMuchHits = 0; Double_t chi2much = 1000; - - if(iMuchTrack > -1){ + + if (iMuchTrack > -1) { CbmMuchTrack* muchTrack = (CbmMuchTrack*) fMuchTracks->At(iMuchTrack); - if(muchTrack){ + if (muchTrack) { nMuchHits = muchTrack->GetNofHits(); - if(muchTrack->GetNDF()!=0)chi2much=muchTrack->GetChiSq()/muchTrack->GetNDF(); + if (muchTrack->GetNDF() != 0) + chi2much = muchTrack->GetChiSq() / muchTrack->GetNDF(); } } - -//----------------- TRD track parameters - - Int_t nTrdHits = 0; - Double_t chi2trd = 1000; - - if (iTrdTrack > -1){ + + //----------------- TRD track parameters + + Int_t nTrdHits = 0; + Double_t chi2trd = 1000; + + if (iTrdTrack > -1) { CbmTrdTrack* trdTrack = (CbmTrdTrack*) fTrdTracks->At(iTrdTrack); - if(trdTrack){ - nTrdHits = trdTrack->GetNofHits(); - if(trdTrack->GetNDF()!=0)chi2trd=trdTrack->GetChiSq()/trdTrack->GetNDF(); // study of TRD chi2 ! + if (trdTrack) { + nTrdHits = trdTrack->GetNofHits(); + if (trdTrack->GetNDF() != 0) + chi2trd = + trdTrack->GetChiSq() / trdTrack->GetNDF(); // study of TRD chi2 ! } - } - -//----------------- TOF hit parameters - - Int_t nTofHits = 0; - Double_t mass = -1000; - - if(iTofHit > -1){ - CbmTofHit* th = (CbmTofHit*)fTofHit->At(iTofHit); - if(th){ - nTofHits = 1; - - Double_t time = th->GetTime(); - Double_t beta = globalTrack->GetLength()*0.01/(time*1e-9*TMath::C()); - - TVector3 momL; - - FairTrackParam* stpl = (FairTrackParam*)globalTrack->GetParamLast(); - stpl->Momentum(momL); - - if(beta!=0)mass = momL.Mag()*momL.Mag()*(1./beta/beta-1.); + } + + //----------------- TOF hit parameters + + Int_t nTofHits = 0; + Double_t mass = -1000; + + if (iTofHit > -1) { + CbmTofHit* th = (CbmTofHit*) fTofHit->At(iTofHit); + if (th) { + nTofHits = 1; + + Double_t time = th->GetTime(); + Double_t beta = + globalTrack->GetLength() * 0.01 / (time * 1e-9 * TMath::C()); + + TVector3 momL; + + FairTrackParam* stpl = (FairTrackParam*) globalTrack->GetParamLast(); + stpl->Momentum(momL); + + if (beta != 0) mass = momL.Mag() * momL.Mag() * (1. / beta / beta - 1.); } } -//----------------- STS MC matching - - Int_t isMu = 0; + //----------------- STS MC matching + + Int_t isMu = 0; Int_t stsPDG = 0; - - if(fUseMC){ - - CbmTrackMatchNew* stsMatch = (CbmTrackMatchNew*) fStsTrackMatches->At(iStsTrack); - - if(stsMatch){ - if (stsMatch->GetNofLinks() != 0){ - int stsMcTrackId = stsMatch->GetMatchedLink().GetIndex(); - CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTracks->At(stsMcTrackId); - if(mcTrack){ - int pdg = TMath::Abs(mcTrack->GetPdgCode()); - stsPDG = pdg; - if(mcTrack->GetGeantProcessId()==kPPrimary && pdg==13){ - isMu=1; - if(mcTrack->GetCharge()<0){ - muMn_reco.Mu=kTRUE; - - if(chi2vertex<=fChi2VertexCut){ - muMn_reco.Chi2V=kTRUE; - - if(nStsHits>=fNofStsCut && chi2sts<=fChi2StsCut){ - muMn_reco.Nsts=kTRUE; - - if(nMuchHits>=fNofMuchCut && chi2much<=fChi2MuchCut){ - muMn_reco.Nmuch=kTRUE; - - if(nTrdHits>=fNofTrdCut){ - muMn_reco.Ntrd=kTRUE; - - if(mass > (p0min+p1min*momentum+p2min*momentum*momentum) && - mass < (p0max+p1max*momentum+p2max*momentum*momentum) && fAnnCut < 0)muMn_reco.Ntof=kTRUE; - } - } - } - } - }else{ - muPl_reco.Mu=kTRUE; - if(chi2vertex<=fChi2VertexCut){ - muPl_reco.Chi2V=kTRUE; - - if(nStsHits>=fNofStsCut && chi2sts<=fChi2StsCut){ - muPl_reco.Nsts=kTRUE; - - if(nMuchHits>=fNofMuchCut && chi2much<=fChi2MuchCut){ - muPl_reco.Nmuch=kTRUE; - - if(nTrdHits>=fNofTrdCut){ - muPl_reco.Ntrd=kTRUE; - - if(mass > (p0min+p1min*momentum+p2min*momentum*momentum) && - mass < (p0max+p1max*momentum+p2max*momentum*momentum) && fAnnCut < 0)muPl_reco.Ntof=kTRUE; - } - } - } - } - } - } + + if (fUseMC) { + + CbmTrackMatchNew* stsMatch = + (CbmTrackMatchNew*) fStsTrackMatches->At(iStsTrack); + + if (stsMatch) { + if (stsMatch->GetNofLinks() != 0) { + int stsMcTrackId = stsMatch->GetMatchedLink().GetIndex(); + CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTracks->At(stsMcTrackId); + if (mcTrack) { + int pdg = TMath::Abs(mcTrack->GetPdgCode()); + stsPDG = pdg; + if (mcTrack->GetGeantProcessId() == kPPrimary && pdg == 13) { + isMu = 1; + if (mcTrack->GetCharge() < 0) { + muMn_reco.Mu = kTRUE; + + if (chi2vertex <= fChi2VertexCut) { + muMn_reco.Chi2V = kTRUE; + + if (nStsHits >= fNofStsCut && chi2sts <= fChi2StsCut) { + muMn_reco.Nsts = kTRUE; + + if (nMuchHits >= fNofMuchCut && chi2much <= fChi2MuchCut) { + muMn_reco.Nmuch = kTRUE; + + if (nTrdHits >= fNofTrdCut) { + muMn_reco.Ntrd = kTRUE; + + if (mass > (p0min + p1min * momentum + + p2min * momentum * momentum) + && mass < (p0max + p1max * momentum + + p2max * momentum * momentum) + && fAnnCut < 0) + muMn_reco.Ntof = kTRUE; + } + } + } + } + } else { + muPl_reco.Mu = kTRUE; + if (chi2vertex <= fChi2VertexCut) { + muPl_reco.Chi2V = kTRUE; + + if (nStsHits >= fNofStsCut && chi2sts <= fChi2StsCut) { + muPl_reco.Nsts = kTRUE; + + if (nMuchHits >= fNofMuchCut && chi2much <= fChi2MuchCut) { + muPl_reco.Nmuch = kTRUE; + + if (nTrdHits >= fNofTrdCut) { + muPl_reco.Ntrd = kTRUE; + + if (mass > (p0min + p1min * momentum + + p2min * momentum * momentum) + && mass < (p0max + p1max * momentum + + p2max * momentum * momentum) + && fAnnCut < 0) + muPl_reco.Ntof = kTRUE; + } + } + } + } + } + } + } + } } } - } - } - if(muPl_reco.Mu && muMn_reco.Mu) signal_reco.Mu=kTRUE; - if(muPl_reco.Chi2V && muMn_reco.Chi2V) signal_reco.Chi2V=kTRUE; - if(muPl_reco.Nsts && muMn_reco.Nsts) signal_reco.Nsts=kTRUE; - if(muPl_reco.Nmuch && muMn_reco.Nmuch) signal_reco.Nmuch=kTRUE; - if(muPl_reco.Ntrd && muMn_reco.Ntrd) signal_reco.Ntrd=kTRUE; - if(muPl_reco.Ntof && muMn_reco.Ntof && fAnnCut < 0) signal_reco.Ntof=kTRUE; - - if(chi2vertex<=fChi2VertexCut){ - BgSup[1]->Fill(momentum); - - if(nStsHits>=fNofStsCut && chi2sts<=fChi2StsCut){ - BgSup[2]->Fill(momentum); - - if(nMuchHits>=fNofMuchCut && chi2much<=fChi2MuchCut){ - BgSup[3]->Fill(momentum); - - if(nTrdHits>=fNofTrdCut){ - BgSup[4]->Fill(momentum); - - if(mass > (p0min+p1min*momentum+p2min*momentum*momentum) && - mass < (p0max+p1max*momentum+p2max*momentum*momentum)&& - fAnnCut < 0)BgSup[5]->Fill(momentum); - } - } - } - } - -//----------------- Muon track candidates + if (muPl_reco.Mu && muMn_reco.Mu) signal_reco.Mu = kTRUE; + if (muPl_reco.Chi2V && muMn_reco.Chi2V) signal_reco.Chi2V = kTRUE; + if (muPl_reco.Nsts && muMn_reco.Nsts) signal_reco.Nsts = kTRUE; + if (muPl_reco.Nmuch && muMn_reco.Nmuch) signal_reco.Nmuch = kTRUE; + if (muPl_reco.Ntrd && muMn_reco.Ntrd) signal_reco.Ntrd = kTRUE; + if (muPl_reco.Ntof && muMn_reco.Ntof && fAnnCut < 0) + signal_reco.Ntof = kTRUE; + + if (chi2vertex <= fChi2VertexCut) { + BgSup[1]->Fill(momentum); + + if (nStsHits >= fNofStsCut && chi2sts <= fChi2StsCut) { + BgSup[2]->Fill(momentum); + + if (nMuchHits >= fNofMuchCut && chi2much <= fChi2MuchCut) { + BgSup[3]->Fill(momentum); + + if (nTrdHits >= fNofTrdCut) { + BgSup[4]->Fill(momentum); + + if (mass > (p0min + p1min * momentum + p2min * momentum * momentum) + && mass + < (p0max + p1max * momentum + p2max * momentum * momentum) + && fAnnCut < 0) + BgSup[5]->Fill(momentum); + } + } + } + } + + //----------------- Muon track candidates Double_t id = -1; - - if(!fUseCuts && fAnnCut < 0)analysis=kTRUE; - else if(fUseCuts) - { - if(nStsHits >= fNofStsCut && - nMuchHits >= fNofMuchCut && - nTrdHits >= fNofTrdCut && - nTofHits >= 1 && - chi2vertex <= fChi2VertexCut && - chi2sts <= fChi2StsCut && - chi2much <= fChi2MuchCut && - (mass > (p0min+p1min*momentum+p2min*momentum*momentum) && - mass < (p0max+p1max*momentum+p2max*momentum*momentum))) - analysis=kTRUE; - }else if(!fUseCuts && fAnnCut > 0) { - id = CalculateAnnValue(mom.P(), mass, chi2vertex, chi2sts, chi2much, chi2trd, nStsHits, nMuchHits, nTrdHits, nTofHits); - if(id > fAnnCut){ - analysis=kTRUE; - BgSup[5]->Fill(momentum); - if(isMu==1 && q>0)muPl_reco.Ntof=kTRUE; - else if(isMu==1 && q<0)muMn_reco.Ntof=kTRUE; - if(muPl_reco.Ntof && muMn_reco.Ntof)signal_reco.Ntof=kTRUE; + + if (!fUseCuts && fAnnCut < 0) + analysis = kTRUE; + else if (fUseCuts) { + if (nStsHits >= fNofStsCut && nMuchHits >= fNofMuchCut + && nTrdHits >= fNofTrdCut && nTofHits >= 1 + && chi2vertex <= fChi2VertexCut && chi2sts <= fChi2StsCut + && chi2much <= fChi2MuchCut + && (mass > (p0min + p1min * momentum + p2min * momentum * momentum) + && mass + < (p0max + p1max * momentum + p2max * momentum * momentum))) + analysis = kTRUE; + } else if (!fUseCuts && fAnnCut > 0) { + id = CalculateAnnValue(mom.P(), + mass, + chi2vertex, + chi2sts, + chi2much, + chi2trd, + nStsHits, + nMuchHits, + nTrdHits, + nTofHits); + if (id > fAnnCut) { + analysis = kTRUE; + BgSup[5]->Fill(momentum); + if (isMu == 1 && q > 0) + muPl_reco.Ntof = kTRUE; + else if (isMu == 1 && q < 0) + muMn_reco.Ntof = kTRUE; + if (muPl_reco.Ntof && muMn_reco.Ntof) signal_reco.Ntof = kTRUE; } } - - if(!analysis)continue; - - CbmAnaMuonCandidate* mu; - if(q>0) - { - new((*fMuPlus)[iMuPlus++]) CbmAnaMuonCandidate(); - mu = (CbmAnaMuonCandidate*) (*fMuPlus)[iMuPlus-1]; - } - else - { - new((*fMuMinus)[iMuMinus++]) CbmAnaMuonCandidate(); - mu = (CbmAnaMuonCandidate*) (*fMuMinus)[iMuMinus-1]; - } - - mu->SetNStsHits(nStsHits); - mu->SetNMuchHits(nMuchHits); - mu->SetNTrdHits(nTrdHits); - mu->SetNTofHits(nTofHits); - - mu->SetChiToVertex(chi2vertex); - mu->SetChiSts(chi2sts); - mu->SetChiMuch(chi2much); - mu->SetChiTrd(chi2trd); - - mu->SetTrueMu(isMu); - mu->SetStsPdg(stsPDG); - mu->SetTofM(mass); - mu->SetSign(q); - mu->SetMomentum(mom); - - if(!fUseCuts && fAnnCut < 0)id = CalculateAnnValue(mu); - mu->SetAnnId(id); - } - -//----------------- Fill YPtM spectrum + + if (!analysis) continue; + + CbmAnaMuonCandidate* mu; + if (q > 0) { + new ((*fMuPlus)[iMuPlus++]) CbmAnaMuonCandidate(); + mu = (CbmAnaMuonCandidate*) (*fMuPlus)[iMuPlus - 1]; + } else { + new ((*fMuMinus)[iMuMinus++]) CbmAnaMuonCandidate(); + mu = (CbmAnaMuonCandidate*) (*fMuMinus)[iMuMinus - 1]; + } + + mu->SetNStsHits(nStsHits); + mu->SetNMuchHits(nMuchHits); + mu->SetNTrdHits(nTrdHits); + mu->SetNTofHits(nTofHits); + + mu->SetChiToVertex(chi2vertex); + mu->SetChiSts(chi2sts); + mu->SetChiMuch(chi2much); + mu->SetChiTrd(chi2trd); + + mu->SetTrueMu(isMu); + mu->SetStsPdg(stsPDG); + mu->SetTofM(mass); + mu->SetSign(q); + mu->SetMomentum(mom); + + if (!fUseCuts && fAnnCut < 0) id = CalculateAnnValue(mu); + mu->SetAnnId(id); + } + + //----------------- Fill YPtM spectrum int NofPlus = fMuPlus->GetEntriesFast(); int NofMinus = fMuMinus->GetEntriesFast(); - for (int iPart=0; iPart<NofPlus; iPart++){ - CbmAnaMuonCandidate *mu_pl = (CbmAnaMuonCandidate*)fMuPlus->At(iPart); - TLorentzVector *P_pl = mu_pl->GetMomentum(); - for (int jPart=0; jPart<NofMinus; jPart++){ - CbmAnaMuonCandidate *mu_mn = (CbmAnaMuonCandidate*)fMuMinus->At(jPart); - TLorentzVector *P_mn = mu_mn->GetMomentum(); - M = *P_pl + *P_mn; - YPtM->Fill(M.Rapidity(),M.Pt(),M.M()); + for (int iPart = 0; iPart < NofPlus; iPart++) { + CbmAnaMuonCandidate* mu_pl = (CbmAnaMuonCandidate*) fMuPlus->At(iPart); + TLorentzVector* P_pl = mu_pl->GetMomentum(); + for (int jPart = 0; jPart < NofMinus; jPart++) { + CbmAnaMuonCandidate* mu_mn = (CbmAnaMuonCandidate*) fMuMinus->At(jPart); + TLorentzVector* P_mn = mu_mn->GetMomentum(); + M = *P_pl + *P_mn; + YPtM->Fill(M.Rapidity(), M.Pt(), M.M()); } } -//----------------- Fill acceptance and efficiency spectra for MC signal muons and signal signal -//----------------- Fill YPt spectra for signal signal - - if(fPlutoFileName != ""){ + //----------------- Fill acceptance and efficiency spectra for MC signal muons and signal signal + //----------------- Fill YPt spectra for signal signal + + if (fPlutoFileName != "") { fInputTree->GetEntry(fEvent); - Int_t NofPart = fParticles->GetEntriesFast(); - PParticle* Part1 = (PParticle*) fParticles->At(NofPart-1); - PParticle* Part2 = (PParticle*) fParticles->At(NofPart-2); + Int_t NofPart = fParticles->GetEntriesFast(); + PParticle* Part1 = (PParticle*) fParticles->At(NofPart - 1); + PParticle* Part2 = (PParticle*) fParticles->At(NofPart - 2); TLorentzVector mom1 = Part1->Vect4(); TLorentzVector mom2 = Part2->Vect4(); - - if(muPl_mc.Mu){ - FillProfile(acc_P[0][1],mom2.P(),muPl_mc.Nsts); - FillProfile(acc_P[1][1],mom2.P(),muPl_mc.Nmuch); - FillProfile(acc_P[2][1],mom2.P(),muPl_mc.Ntrd); - FillProfile(acc_P[3][1],mom2.P(),muPl_mc.Ntof); - FillProfile(acc_Theta[0][1],mom2.Theta()*TMath::RadToDeg(),muPl_mc.Nsts); - FillProfile(acc_Theta[1][1],mom2.Theta()*TMath::RadToDeg(),muPl_mc.Nmuch); - FillProfile(acc_Theta[2][1],mom2.Theta()*TMath::RadToDeg(),muPl_mc.Ntrd); - FillProfile(acc_Theta[3][1],mom2.Theta()*TMath::RadToDeg(),muPl_mc.Ntof); - if(muPl_mc.Nsts) {FillProfile(effReco_P[0][1],mom2.P(),muPl_reco.Nsts); FillProfile(effReco_Theta[0][1],mom2.Theta()*TMath::RadToDeg(),muPl_reco.Nsts);} - if(muPl_mc.Nmuch){FillProfile(effReco_P[1][1],mom2.P(),muPl_reco.Nmuch);FillProfile(effReco_Theta[1][1],mom2.Theta()*TMath::RadToDeg(),muPl_reco.Nmuch);} - if(muPl_mc.Ntrd) {FillProfile(effReco_P[2][1],mom2.P(),muPl_reco.Ntrd); FillProfile(effReco_Theta[2][1],mom2.Theta()*TMath::RadToDeg(),muPl_reco.Ntrd);} - if(muPl_mc.Ntof) {FillProfile(effReco_P[3][1],mom2.P(),muPl_reco.Ntof); FillProfile(effReco_Theta[3][1],mom2.Theta()*TMath::RadToDeg(),muPl_reco.Ntof);} + + if (muPl_mc.Mu) { + FillProfile(acc_P[0][1], mom2.P(), muPl_mc.Nsts); + FillProfile(acc_P[1][1], mom2.P(), muPl_mc.Nmuch); + FillProfile(acc_P[2][1], mom2.P(), muPl_mc.Ntrd); + FillProfile(acc_P[3][1], mom2.P(), muPl_mc.Ntof); + FillProfile( + acc_Theta[0][1], mom2.Theta() * TMath::RadToDeg(), muPl_mc.Nsts); + FillProfile( + acc_Theta[1][1], mom2.Theta() * TMath::RadToDeg(), muPl_mc.Nmuch); + FillProfile( + acc_Theta[2][1], mom2.Theta() * TMath::RadToDeg(), muPl_mc.Ntrd); + FillProfile( + acc_Theta[3][1], mom2.Theta() * TMath::RadToDeg(), muPl_mc.Ntof); + if (muPl_mc.Nsts) { + FillProfile(effReco_P[0][1], mom2.P(), muPl_reco.Nsts); + FillProfile(effReco_Theta[0][1], + mom2.Theta() * TMath::RadToDeg(), + muPl_reco.Nsts); + } + if (muPl_mc.Nmuch) { + FillProfile(effReco_P[1][1], mom2.P(), muPl_reco.Nmuch); + FillProfile(effReco_Theta[1][1], + mom2.Theta() * TMath::RadToDeg(), + muPl_reco.Nmuch); + } + if (muPl_mc.Ntrd) { + FillProfile(effReco_P[2][1], mom2.P(), muPl_reco.Ntrd); + FillProfile(effReco_Theta[2][1], + mom2.Theta() * TMath::RadToDeg(), + muPl_reco.Ntrd); + } + if (muPl_mc.Ntof) { + FillProfile(effReco_P[3][1], mom2.P(), muPl_reco.Ntof); + FillProfile(effReco_Theta[3][1], + mom2.Theta() * TMath::RadToDeg(), + muPl_reco.Ntof); + } + } + + if (muPl_reco.Mu) { + FillProfile(eff4pi_P[0][1], mom2.P(), muPl_reco.Chi2V); + FillProfile(eff4pi_P[1][1], mom2.P(), muPl_reco.Nsts); + FillProfile(eff4pi_P[2][1], mom2.P(), muPl_reco.Nmuch); + FillProfile(eff4pi_P[3][1], mom2.P(), muPl_reco.Ntrd); + FillProfile(eff4pi_P[4][1], mom2.P(), muPl_reco.Ntof); + FillProfile( + eff4pi_Theta[0][1], mom2.Theta() * TMath::RadToDeg(), muPl_reco.Chi2V); + FillProfile( + eff4pi_Theta[1][1], mom2.Theta() * TMath::RadToDeg(), muPl_reco.Nsts); + FillProfile( + eff4pi_Theta[2][1], mom2.Theta() * TMath::RadToDeg(), muPl_reco.Nmuch); + FillProfile( + eff4pi_Theta[3][1], mom2.Theta() * TMath::RadToDeg(), muPl_reco.Ntrd); + FillProfile( + eff4pi_Theta[4][1], mom2.Theta() * TMath::RadToDeg(), muPl_reco.Ntof); } - - if(muPl_reco.Mu){ - FillProfile(eff4pi_P[0][1],mom2.P(),muPl_reco.Chi2V); - FillProfile(eff4pi_P[1][1],mom2.P(),muPl_reco.Nsts); - FillProfile(eff4pi_P[2][1],mom2.P(),muPl_reco.Nmuch); - FillProfile(eff4pi_P[3][1],mom2.P(),muPl_reco.Ntrd); - FillProfile(eff4pi_P[4][1],mom2.P(),muPl_reco.Ntof); - FillProfile(eff4pi_Theta[0][1],mom2.Theta()*TMath::RadToDeg(),muPl_reco.Chi2V); - FillProfile(eff4pi_Theta[1][1],mom2.Theta()*TMath::RadToDeg(),muPl_reco.Nsts); - FillProfile(eff4pi_Theta[2][1],mom2.Theta()*TMath::RadToDeg(),muPl_reco.Nmuch); - FillProfile(eff4pi_Theta[3][1],mom2.Theta()*TMath::RadToDeg(),muPl_reco.Ntrd); - FillProfile(eff4pi_Theta[4][1],mom2.Theta()*TMath::RadToDeg(),muPl_reco.Ntof); + + if (muMn_mc.Mu) { + FillProfile(acc_P[0][2], mom1.P(), muMn_mc.Nsts); + FillProfile(acc_P[1][2], mom1.P(), muMn_mc.Nmuch); + FillProfile(acc_P[2][2], mom1.P(), muMn_mc.Ntrd); + FillProfile(acc_P[3][2], mom1.P(), muMn_mc.Ntof); + FillProfile( + acc_Theta[0][2], mom1.Theta() * TMath::RadToDeg(), muMn_mc.Nsts); + FillProfile( + acc_Theta[1][2], mom1.Theta() * TMath::RadToDeg(), muMn_mc.Nmuch); + FillProfile( + acc_Theta[2][2], mom1.Theta() * TMath::RadToDeg(), muMn_mc.Ntrd); + FillProfile( + acc_Theta[3][2], mom1.Theta() * TMath::RadToDeg(), muMn_mc.Ntof); + if (muMn_mc.Nsts) { + FillProfile(effReco_P[0][2], mom1.P(), muMn_reco.Nsts); + FillProfile(effReco_Theta[0][2], + mom1.Theta() * TMath::RadToDeg(), + muMn_reco.Nsts); + } + if (muMn_mc.Nmuch) { + FillProfile(effReco_P[1][2], mom1.P(), muMn_reco.Nmuch); + FillProfile(effReco_Theta[1][2], + mom1.Theta() * TMath::RadToDeg(), + muMn_reco.Nmuch); + } + if (muMn_mc.Ntrd) { + FillProfile(effReco_P[2][2], mom1.P(), muMn_reco.Ntrd); + FillProfile(effReco_Theta[2][2], + mom1.Theta() * TMath::RadToDeg(), + muMn_reco.Ntrd); + } + if (muMn_mc.Ntof) { + FillProfile(effReco_P[3][2], mom1.P(), muMn_reco.Ntof); + FillProfile(effReco_Theta[3][2], + mom1.Theta() * TMath::RadToDeg(), + muMn_reco.Ntof); + } + } + + if (muMn_reco.Mu) { + FillProfile(eff4pi_P[0][2], mom1.P(), muMn_reco.Chi2V); + FillProfile(eff4pi_P[1][2], mom1.P(), muMn_reco.Nsts); + FillProfile(eff4pi_P[2][2], mom1.P(), muMn_reco.Nmuch); + FillProfile(eff4pi_P[3][2], mom1.P(), muMn_reco.Ntrd); + FillProfile(eff4pi_P[4][2], mom1.P(), muMn_reco.Ntof); + FillProfile( + eff4pi_Theta[0][2], mom1.Theta() * TMath::RadToDeg(), muMn_reco.Chi2V); + FillProfile( + eff4pi_Theta[1][2], mom1.Theta() * TMath::RadToDeg(), muMn_reco.Nsts); + FillProfile( + eff4pi_Theta[2][2], mom1.Theta() * TMath::RadToDeg(), muMn_reco.Nmuch); + FillProfile( + eff4pi_Theta[3][2], mom1.Theta() * TMath::RadToDeg(), muMn_reco.Ntrd); + FillProfile( + eff4pi_Theta[4][2], mom1.Theta() * TMath::RadToDeg(), muMn_reco.Ntof); } - - if(muMn_mc.Mu){ - FillProfile(acc_P[0][2],mom1.P(),muMn_mc.Nsts); - FillProfile(acc_P[1][2],mom1.P(),muMn_mc.Nmuch); - FillProfile(acc_P[2][2],mom1.P(),muMn_mc.Ntrd); - FillProfile(acc_P[3][2],mom1.P(),muMn_mc.Ntof); - FillProfile(acc_Theta[0][2],mom1.Theta()*TMath::RadToDeg(),muMn_mc.Nsts); - FillProfile(acc_Theta[1][2],mom1.Theta()*TMath::RadToDeg(),muMn_mc.Nmuch); - FillProfile(acc_Theta[2][2],mom1.Theta()*TMath::RadToDeg(),muMn_mc.Ntrd); - FillProfile(acc_Theta[3][2],mom1.Theta()*TMath::RadToDeg(),muMn_mc.Ntof); - if(muMn_mc.Nsts) {FillProfile(effReco_P[0][2],mom1.P(),muMn_reco.Nsts); FillProfile(effReco_Theta[0][2],mom1.Theta()*TMath::RadToDeg(),muMn_reco.Nsts);} - if(muMn_mc.Nmuch){FillProfile(effReco_P[1][2],mom1.P(),muMn_reco.Nmuch);FillProfile(effReco_Theta[1][2],mom1.Theta()*TMath::RadToDeg(),muMn_reco.Nmuch);} - if(muMn_mc.Ntrd) {FillProfile(effReco_P[2][2],mom1.P(),muMn_reco.Ntrd); FillProfile(effReco_Theta[2][2],mom1.Theta()*TMath::RadToDeg(),muMn_reco.Ntrd);} - if(muMn_mc.Ntof) {FillProfile(effReco_P[3][2],mom1.P(),muMn_reco.Ntof); FillProfile(effReco_Theta[3][2],mom1.Theta()*TMath::RadToDeg(),muMn_reco.Ntof);} - } - - if(muMn_reco.Mu){ - FillProfile(eff4pi_P[0][2],mom1.P(),muMn_reco.Chi2V); - FillProfile(eff4pi_P[1][2],mom1.P(),muMn_reco.Nsts); - FillProfile(eff4pi_P[2][2],mom1.P(),muMn_reco.Nmuch); - FillProfile(eff4pi_P[3][2],mom1.P(),muMn_reco.Ntrd); - FillProfile(eff4pi_P[4][2],mom1.P(),muMn_reco.Ntof); - FillProfile(eff4pi_Theta[0][2],mom1.Theta()*TMath::RadToDeg(),muMn_reco.Chi2V); - FillProfile(eff4pi_Theta[1][2],mom1.Theta()*TMath::RadToDeg(),muMn_reco.Nsts); - FillProfile(eff4pi_Theta[2][2],mom1.Theta()*TMath::RadToDeg(),muMn_reco.Nmuch); - FillProfile(eff4pi_Theta[3][2],mom1.Theta()*TMath::RadToDeg(),muMn_reco.Ntrd); - FillProfile(eff4pi_Theta[4][2],mom1.Theta()*TMath::RadToDeg(),muMn_reco.Ntof); - } - - TLorentzVector Mom = mom1 + mom2; - if(signal_mc.Mu){ - FillProfile(acc_P[0][0],Mom.P(),signal_mc.Nsts); - FillProfile(acc_P[1][0],Mom.P(),signal_mc.Nmuch); - FillProfile(acc_P[2][0],Mom.P(),signal_mc.Ntrd); - FillProfile(acc_P[3][0],Mom.P(),signal_mc.Ntof); - FillProfile(acc_Theta[0][0],Mom.Theta()*TMath::RadToDeg(),signal_mc.Nsts); - FillProfile(acc_Theta[1][0],Mom.Theta()*TMath::RadToDeg(),signal_mc.Nmuch); - FillProfile(acc_Theta[2][0],Mom.Theta()*TMath::RadToDeg(),signal_mc.Ntrd); - FillProfile(acc_Theta[3][0],Mom.Theta()*TMath::RadToDeg(),signal_mc.Ntof); - if(signal_mc.Nsts){ - FillProfile(effReco_P[0][0],Mom.P(),signal_reco.Nsts); FillProfile(effReco_Theta[0][0],Mom.Theta()*TMath::RadToDeg(),signal_reco.Nsts); - YPt_StsAcc->Fill(Mom.Rapidity(),Mom.Pt()); + + TLorentzVector Mom = mom1 + mom2; + if (signal_mc.Mu) { + FillProfile(acc_P[0][0], Mom.P(), signal_mc.Nsts); + FillProfile(acc_P[1][0], Mom.P(), signal_mc.Nmuch); + FillProfile(acc_P[2][0], Mom.P(), signal_mc.Ntrd); + FillProfile(acc_P[3][0], Mom.P(), signal_mc.Ntof); + FillProfile( + acc_Theta[0][0], Mom.Theta() * TMath::RadToDeg(), signal_mc.Nsts); + FillProfile( + acc_Theta[1][0], Mom.Theta() * TMath::RadToDeg(), signal_mc.Nmuch); + FillProfile( + acc_Theta[2][0], Mom.Theta() * TMath::RadToDeg(), signal_mc.Ntrd); + FillProfile( + acc_Theta[3][0], Mom.Theta() * TMath::RadToDeg(), signal_mc.Ntof); + if (signal_mc.Nsts) { + FillProfile(effReco_P[0][0], Mom.P(), signal_reco.Nsts); + FillProfile(effReco_Theta[0][0], + Mom.Theta() * TMath::RadToDeg(), + signal_reco.Nsts); + YPt_StsAcc->Fill(Mom.Rapidity(), Mom.Pt()); } - if(signal_mc.Nmuch){ - FillProfile(effReco_P[1][0],Mom.P(),signal_reco.Nmuch);FillProfile(effReco_Theta[1][0],Mom.Theta()*TMath::RadToDeg(),signal_reco.Nmuch); - YPt_StsMuchAcc->Fill(Mom.Rapidity(),Mom.Pt()); + if (signal_mc.Nmuch) { + FillProfile(effReco_P[1][0], Mom.P(), signal_reco.Nmuch); + FillProfile(effReco_Theta[1][0], + Mom.Theta() * TMath::RadToDeg(), + signal_reco.Nmuch); + YPt_StsMuchAcc->Fill(Mom.Rapidity(), Mom.Pt()); } - if(signal_mc.Ntrd){ - FillProfile(effReco_P[2][0],Mom.P(),signal_reco.Ntrd); FillProfile(effReco_Theta[2][0],Mom.Theta()*TMath::RadToDeg(),signal_reco.Ntrd); - YPt_StsMuchTrdAcc->Fill(Mom.Rapidity(),Mom.Pt()); + if (signal_mc.Ntrd) { + FillProfile(effReco_P[2][0], Mom.P(), signal_reco.Ntrd); + FillProfile(effReco_Theta[2][0], + Mom.Theta() * TMath::RadToDeg(), + signal_reco.Ntrd); + YPt_StsMuchTrdAcc->Fill(Mom.Rapidity(), Mom.Pt()); + } + if (signal_mc.Ntof) { + FillProfile(effReco_P[3][0], Mom.P(), signal_reco.Ntof); + FillProfile(effReco_Theta[3][0], + Mom.Theta() * TMath::RadToDeg(), + signal_reco.Ntof); + YPt_StsMuchTrdTofAcc->Fill(Mom.Rapidity(), Mom.Pt()); } - if(signal_mc.Ntof){ - FillProfile(effReco_P[3][0],Mom.P(),signal_reco.Ntof); FillProfile(effReco_Theta[3][0],Mom.Theta()*TMath::RadToDeg(),signal_reco.Ntof); - YPt_StsMuchTrdTofAcc->Fill(Mom.Rapidity(),Mom.Pt()); - } } - if(signal_reco.Mu){ - if(signal_reco.Chi2V)YPt_VtxReco->Fill(Mom.Rapidity(),Mom.Pt()); - if(signal_reco.Nsts) YPt_VtxStsReco->Fill(Mom.Rapidity(),Mom.Pt()); - if(signal_reco.Nmuch)YPt_VtxStsMuchReco->Fill(Mom.Rapidity(),Mom.Pt()); - if(signal_reco.Ntrd) YPt_VtxStsMuchTrdReco->Fill(Mom.Rapidity(),Mom.Pt()); - if(signal_reco.Ntof) YPt_VtxStsMuchTrdTofReco->Fill(Mom.Rapidity(),Mom.Pt()); + if (signal_reco.Mu) { + if (signal_reco.Chi2V) YPt_VtxReco->Fill(Mom.Rapidity(), Mom.Pt()); + if (signal_reco.Nsts) YPt_VtxStsReco->Fill(Mom.Rapidity(), Mom.Pt()); + if (signal_reco.Nmuch) YPt_VtxStsMuchReco->Fill(Mom.Rapidity(), Mom.Pt()); + if (signal_reco.Ntrd) + YPt_VtxStsMuchTrdReco->Fill(Mom.Rapidity(), Mom.Pt()); + if (signal_reco.Ntof) + YPt_VtxStsMuchTrdTofReco->Fill(Mom.Rapidity(), Mom.Pt()); } - FillProfile(eff4pi_P[0][0],Mom.P(),signal_reco.Chi2V); - FillProfile(eff4pi_P[1][0],Mom.P(),signal_reco.Nsts); - FillProfile(eff4pi_P[2][0],Mom.P(),signal_reco.Nmuch); - FillProfile(eff4pi_P[3][0],Mom.P(),signal_reco.Ntrd); - FillProfile(eff4pi_P[4][0],Mom.P(),signal_reco.Ntof); - FillProfile(eff4pi_Theta[0][0],Mom.Theta()*TMath::RadToDeg(),signal_reco.Chi2V); - FillProfile(eff4pi_Theta[1][0],Mom.Theta()*TMath::RadToDeg(),signal_reco.Nsts); - FillProfile(eff4pi_Theta[2][0],Mom.Theta()*TMath::RadToDeg(),signal_reco.Nmuch); - FillProfile(eff4pi_Theta[3][0],Mom.Theta()*TMath::RadToDeg(),signal_reco.Ntrd); - FillProfile(eff4pi_Theta[4][0],Mom.Theta()*TMath::RadToDeg(),signal_reco.Ntof); + FillProfile(eff4pi_P[0][0], Mom.P(), signal_reco.Chi2V); + FillProfile(eff4pi_P[1][0], Mom.P(), signal_reco.Nsts); + FillProfile(eff4pi_P[2][0], Mom.P(), signal_reco.Nmuch); + FillProfile(eff4pi_P[3][0], Mom.P(), signal_reco.Ntrd); + FillProfile(eff4pi_P[4][0], Mom.P(), signal_reco.Ntof); + FillProfile( + eff4pi_Theta[0][0], Mom.Theta() * TMath::RadToDeg(), signal_reco.Chi2V); + FillProfile( + eff4pi_Theta[1][0], Mom.Theta() * TMath::RadToDeg(), signal_reco.Nsts); + FillProfile( + eff4pi_Theta[2][0], Mom.Theta() * TMath::RadToDeg(), signal_reco.Nmuch); + FillProfile( + eff4pi_Theta[3][0], Mom.Theta() * TMath::RadToDeg(), signal_reco.Ntrd); + FillProfile( + eff4pi_Theta[4][0], Mom.Theta() * TMath::RadToDeg(), signal_reco.Ntof); } fEvent++; } @@ -862,252 +1126,296 @@ void CbmAnaDimuonAnalysis::Exec(Option_t* /*opt*/){ // used default ANN (16 neurons) weights generated for 8 GeV/c Au beam (UrQMD) and omega->µµ (PLUTO) // sis100_muon_lmvm setup // -Double_t CbmAnaDimuonAnalysis::CalculateAnnValue(CbmAnaMuonCandidate* mu){ - +Double_t CbmAnaDimuonAnalysis::CalculateAnnValue(CbmAnaMuonCandidate* mu) { + Double_t ID = -1; Double_t Chi2Vertex, Chi2STS, Chi2MUCH, Chi2TRD; Int_t NofSTS, NofMUCH, NofTRD, NofTOF, type; Double_t P, M; - TLorentzVector *PP = mu->GetMomentum(); - M = mu->GetTofM(); - P = PP->P(); - Chi2Vertex = mu->GetChiToVertex(); - Chi2STS = mu->GetChiSts(); - Chi2MUCH = mu->GetChiMuch(); - Chi2TRD = mu->GetChiTrd(); - NofSTS = mu->GetNStsHits(); - NofMUCH = mu->GetNMuchHits(); - NofTRD = mu->GetNTrdHits(); - NofTOF = mu->GetNTofHits(); - + TLorentzVector* PP = mu->GetMomentum(); + M = mu->GetTofM(); + P = PP->P(); + Chi2Vertex = mu->GetChiToVertex(); + Chi2STS = mu->GetChiSts(); + Chi2MUCH = mu->GetChiMuch(); + Chi2TRD = mu->GetChiTrd(); + NofSTS = mu->GetNStsHits(); + NofMUCH = mu->GetNMuchHits(); + NofTRD = mu->GetNTrdHits(); + NofTOF = mu->GetNTofHits(); + Double_t MUCHchi2 = 10.; - Double_t STSchi2 = 10.; + Double_t STSchi2 = 10.; Double_t Vchi2 = 10.; - - if((Chi2Vertex <= Vchi2 && Chi2Vertex >= 0) && - (Chi2STS <= STSchi2 && Chi2STS >= 0) && - (Chi2MUCH <= MUCHchi2 && Chi2MUCH >= 0) && - NofTOF > 0){ - Double_t params[9]; - - TTree *simu = new TTree("MonteCarlo", "Filtered Monte Carlo Events"); - simu->Branch("P", &P, "P/D"); - simu->Branch("M", &M, "M/D"); - simu->Branch("Chi2Vertex",&Chi2Vertex,"Chi2Vertex/D"); - simu->Branch("Chi2STS", &Chi2STS, "Chi2STS/D"); - simu->Branch("Chi2MUCH", &Chi2MUCH, "Chi2MUCH/D"); - simu->Branch("Chi2TRD", &Chi2TRD, "Chi2TRD/D"); - simu->Branch("NofSTS", &NofSTS, "NofSTS/I"); - simu->Branch("NofMUCH",&NofMUCH,"NofMUCH/I"); - simu->Branch("NofTRD", &NofTRD, "NofTRD/I"); - simu->Branch("type", &type, "type/I"); - - simu->Fill(); - - TMultiLayerPerceptron *mlp = new TMultiLayerPerceptron("@P,@M,@Chi2Vertex,@Chi2STS,@Chi2MUCH,@Chi2TRD,@NofSTS,@NofMUCH,@NofTRD:16:type",simu); - mlp->LoadWeights(fFileAnnName); - - params[0] = P; - params[1] = M; - params[2] = Chi2Vertex; - params[3] = Chi2STS; - params[4] = Chi2MUCH; - params[5] = Chi2TRD; - params[6] = (Double_t)NofSTS; - params[7] = (Double_t)NofMUCH; - params[8] = (Double_t)NofTRD; - ID = mlp->Evaluate(0, params); - - mlp->Delete(); - simu->Delete(); -} - + + if ((Chi2Vertex <= Vchi2 && Chi2Vertex >= 0) + && (Chi2STS <= STSchi2 && Chi2STS >= 0) + && (Chi2MUCH <= MUCHchi2 && Chi2MUCH >= 0) && NofTOF > 0) { + Double_t params[9]; + + TTree* simu = new TTree("MonteCarlo", "Filtered Monte Carlo Events"); + simu->Branch("P", &P, "P/D"); + simu->Branch("M", &M, "M/D"); + simu->Branch("Chi2Vertex", &Chi2Vertex, "Chi2Vertex/D"); + simu->Branch("Chi2STS", &Chi2STS, "Chi2STS/D"); + simu->Branch("Chi2MUCH", &Chi2MUCH, "Chi2MUCH/D"); + simu->Branch("Chi2TRD", &Chi2TRD, "Chi2TRD/D"); + simu->Branch("NofSTS", &NofSTS, "NofSTS/I"); + simu->Branch("NofMUCH", &NofMUCH, "NofMUCH/I"); + simu->Branch("NofTRD", &NofTRD, "NofTRD/I"); + simu->Branch("type", &type, "type/I"); + + simu->Fill(); + + TMultiLayerPerceptron* mlp = + new TMultiLayerPerceptron("@P,@M,@Chi2Vertex,@Chi2STS,@Chi2MUCH,@Chi2TRD," + "@NofSTS,@NofMUCH,@NofTRD:16:type", + simu); + mlp->LoadWeights(fFileAnnName); + + params[0] = P; + params[1] = M; + params[2] = Chi2Vertex; + params[3] = Chi2STS; + params[4] = Chi2MUCH; + params[5] = Chi2TRD; + params[6] = (Double_t) NofSTS; + params[7] = (Double_t) NofMUCH; + params[8] = (Double_t) NofTRD; + ID = mlp->Evaluate(0, params); + + mlp->Delete(); + simu->Delete(); + } + return ID; } // ------------------------------------------------------------------------- -Double_t CbmAnaDimuonAnalysis::CalculateAnnValue(Double_t P, Double_t M, Double_t Chi2Vertex, Double_t Chi2STS, Double_t Chi2MUCH, Double_t Chi2TRD, Int_t NofSTS, Int_t NofMUCH, Int_t NofTRD, Int_t NofTOF){ +Double_t CbmAnaDimuonAnalysis::CalculateAnnValue(Double_t P, + Double_t M, + Double_t Chi2Vertex, + Double_t Chi2STS, + Double_t Chi2MUCH, + Double_t Chi2TRD, + Int_t NofSTS, + Int_t NofMUCH, + Int_t NofTRD, + Int_t NofTOF) { Double_t ID = -1; Double_t MUCHchi2 = 10.; - Double_t STSchi2 = 10.; + Double_t STSchi2 = 10.; Double_t Vchi2 = 10.; - - if((Chi2Vertex <= Vchi2 && Chi2Vertex >= 0) && - (Chi2STS <= STSchi2 && Chi2STS >= 0) && - (Chi2MUCH <= MUCHchi2 && Chi2MUCH >= 0) && - NofTOF > 0){ - Int_t type; - Double_t params[9]; - - TTree *simu = new TTree("MonteCarlo", "Filtered Monte Carlo Events"); - simu->Branch("P", &P, "P/D"); - simu->Branch("M", &M, "M/D"); - simu->Branch("Chi2Vertex",&Chi2Vertex,"Chi2Vertex/D"); - simu->Branch("Chi2STS", &Chi2STS, "Chi2STS/D"); - simu->Branch("Chi2MUCH", &Chi2MUCH, "Chi2MUCH/D"); - simu->Branch("Chi2TRD", &Chi2TRD, "Chi2TRD/D"); - simu->Branch("NofSTS", &NofSTS, "NofSTS/I"); - simu->Branch("NofMUCH",&NofMUCH,"NofMUCH/I"); - simu->Branch("NofTRD", &NofTRD, "NofTRD/I"); - simu->Branch("type", &type, "type/I"); - - simu->Fill(); - - TString input; - input.Form("@P,@M,@Chi2Vertex,@Chi2STS,@Chi2MUCH,@Chi2TRD,@NofSTS,@NofMUCH,@NofTRD:%d:type",fNeurons); - TMultiLayerPerceptron *mlp = new TMultiLayerPerceptron(input.Data(),simu); - mlp->LoadWeights(fFileAnnName); - - params[0] = P; - params[1] = M; - params[2] = Chi2Vertex; - params[3] = Chi2STS; - params[4] = Chi2MUCH; - params[5] = Chi2TRD; - - params[6] = (Double_t)NofSTS; - params[7] = (Double_t)NofMUCH; - params[8] = (Double_t)NofTRD; - - ID = mlp->Evaluate(0, params); - - mlp->Delete(); - simu->Delete(); + + if ((Chi2Vertex <= Vchi2 && Chi2Vertex >= 0) + && (Chi2STS <= STSchi2 && Chi2STS >= 0) + && (Chi2MUCH <= MUCHchi2 && Chi2MUCH >= 0) && NofTOF > 0) { + Int_t type; + Double_t params[9]; + + TTree* simu = new TTree("MonteCarlo", "Filtered Monte Carlo Events"); + simu->Branch("P", &P, "P/D"); + simu->Branch("M", &M, "M/D"); + simu->Branch("Chi2Vertex", &Chi2Vertex, "Chi2Vertex/D"); + simu->Branch("Chi2STS", &Chi2STS, "Chi2STS/D"); + simu->Branch("Chi2MUCH", &Chi2MUCH, "Chi2MUCH/D"); + simu->Branch("Chi2TRD", &Chi2TRD, "Chi2TRD/D"); + simu->Branch("NofSTS", &NofSTS, "NofSTS/I"); + simu->Branch("NofMUCH", &NofMUCH, "NofMUCH/I"); + simu->Branch("NofTRD", &NofTRD, "NofTRD/I"); + simu->Branch("type", &type, "type/I"); + + simu->Fill(); + + TString input; + input.Form("@P,@M,@Chi2Vertex,@Chi2STS,@Chi2MUCH,@Chi2TRD,@NofSTS,@NofMUCH," + "@NofTRD:%d:type", + fNeurons); + TMultiLayerPerceptron* mlp = new TMultiLayerPerceptron(input.Data(), simu); + mlp->LoadWeights(fFileAnnName); + + params[0] = P; + params[1] = M; + params[2] = Chi2Vertex; + params[3] = Chi2STS; + params[4] = Chi2MUCH; + params[5] = Chi2TRD; + + params[6] = (Double_t) NofSTS; + params[7] = (Double_t) NofMUCH; + params[8] = (Double_t) NofTRD; + + ID = mlp->Evaluate(0, params); + + mlp->Delete(); + simu->Delete(); } return ID; - } // ------------------------------------------------------------------------- -void CbmAnaDimuonAnalysis::FillProfile(TProfile *profile, Double_t param, Bool_t trigger){ - if(trigger)profile->Fill(param,100);else profile->Fill(param,0); +void CbmAnaDimuonAnalysis::FillProfile(TProfile* profile, + Double_t param, + Bool_t trigger) { + if (trigger) + profile->Fill(param, 100); + else + profile->Fill(param, 0); } // ----- Public method Finish ------------------------------------------ -void CbmAnaDimuonAnalysis::Finish(){ - +void CbmAnaDimuonAnalysis::Finish() { + TString name; - if(fPlutoFileName != ""){ - - if(fAnnCut > 0)name = Form("YPt_histo_ANN_%1.2f.root", fAnnCut); - else name = "YPt_histo.root"; - - TFile* f = new TFile(name,"recreate"); - - YPt_StsAcc->Scale(1./(Double_t)fEvent); - YPt_StsMuchAcc->Scale(1./(Double_t)fEvent); - YPt_StsMuchTrdAcc->Scale(1./(Double_t)fEvent); - YPt_StsMuchTrdTofAcc->Scale(1./(Double_t)fEvent); - - YPt_VtxReco->Scale(1./(Double_t)fEvent); - YPt_VtxStsReco->Scale(1./(Double_t)fEvent); - YPt_VtxStsMuchReco->Scale(1./(Double_t)fEvent); - YPt_VtxStsMuchTrdReco->Scale(1./(Double_t)fEvent); - YPt_VtxStsMuchTrdTofReco->Scale(1./(Double_t)fEvent); - - YPt_pluto->Write(); - YPt_StsAcc->Write(); - YPt_StsMuchAcc->Write(); - YPt_StsMuchTrdAcc->Write(); - YPt_StsMuchTrdTofAcc->Write(); - - if(fAnnCut < 0){ - YPt_VtxReco->Write(); - YPt_VtxStsReco->Write(); - YPt_VtxStsMuchReco->Write(); - YPt_VtxStsMuchTrdReco->Write(); - } - YPt_VtxStsMuchTrdTofReco->Write(); + if (fPlutoFileName != "") { - f->Close(); + if (fAnnCut > 0) + name = Form("YPt_histo_ANN_%1.2f.root", fAnnCut); + else + name = "YPt_histo.root"; - if(fAnnCut > 0)name = Form("eff_histo_ANN_%1.2f.root", fAnnCut); - else name = "eff_histo.root"; - TFile* ff = new TFile(name,"recreate"); - - TDirectory *dir1 = ff->mkdir("mu Plus"); - dir1->cd(); - - TDirectory *dir1a = dir1->mkdir("accepted mu Plus"); - dir1a->cd(); - for(int j=0; j<4; j++)acc_P[j][1]->Write(); - for(int j=0; j<4; j++)acc_Theta[j][1]->Write(); - - dir1->cd(); - TDirectory *dir1b = dir1->mkdir("reconstructed mu Plus"); - dir1b->cd(); - if(fAnnCut < 0){ - for(int j=0; j<3; j++)effReco_P[j][1]->Write(); - for(int j=0; j<3; j++)effReco_Theta[j][1]->Write(); - for(int j=0; j<4; j++)eff4pi_P[j][1]->Write(); - for(int j=0; j<4; j++)eff4pi_Theta[j][1]->Write(); - } - effReco_P[3][1]->Write(); - effReco_Theta[3][1]->Write(); - eff4pi_P[4][1]->Write(); - eff4pi_Theta[4][1]->Write(); - - TDirectory *dir2 = ff->mkdir("mu Minus"); - dir2->cd(); - - TDirectory *dir2a = dir2->mkdir("accepted mu Minus"); - dir2a->cd(); - for(int j=0; j<4; j++)acc_P[j][2]->Write(); - for(int j=0; j<4; j++)acc_Theta[j][2]->Write(); - - dir2->cd(); - TDirectory *dir2b = dir2->mkdir("reconstructed mu Minus"); - dir2b->cd(); - if(fAnnCut < 0){ - for(int j=0; j<3; j++)effReco_P[j][2]->Write(); - for(int j=0; j<3; j++)effReco_Theta[j][2]->Write(); - for(int j=0; j<4; j++)eff4pi_P[j][2]->Write(); - for(int j=0; j<4; j++)eff4pi_Theta[j][2]->Write(); - } - effReco_P[3][2]->Write(); - effReco_Theta[3][2]->Write(); - eff4pi_P[4][2]->Write(); - eff4pi_Theta[4][2]->Write(); - - TDirectory *dir3 = ff->mkdir("signal"); - dir3->cd(); - - TDirectory *dir3a = dir3->mkdir("accepted signal"); - dir3a->cd(); - for(int j=0; j<4; j++)acc_P[j][0]->Write(); - for(int j=0; j<4; j++)acc_Theta[j][0]->Write(); - - dir3->cd(); - TDirectory *dir3b = dir3->mkdir("reconstructed signal"); - dir3b->cd(); - if(fAnnCut < 0){ - for(int j=0; j<3; j++)effReco_P[j][0]->Write(); - for(int j=0; j<3; j++)effReco_Theta[j][0]->Write(); - for(int j=0; j<4; j++)eff4pi_P[j][0]->Write(); - for(int j=0; j<4; j++)eff4pi_Theta[j][0]->Write(); - } - effReco_P[3][0]->Write(); - effReco_Theta[3][0]->Write(); - eff4pi_P[4][0]->Write(); - eff4pi_Theta[4][0]->Write(); - - ff->Close(); - // fPlutoFile->Close(); - }else{ - if(fAnnCut > 0)name = Form("sup_histo_ANN_%1.2f.root", fAnnCut); - else name = "sup_histo.root"; - TFile* f = new TFile(name,"recreate"); - - if(fAnnCut < 0)for(int j=0; j<5; j++)BgSup[j]->Write(); - BgSup[5]->Write(); - f->Close(); + TFile* f = new TFile(name, "recreate"); + + YPt_StsAcc->Scale(1. / (Double_t) fEvent); + YPt_StsMuchAcc->Scale(1. / (Double_t) fEvent); + YPt_StsMuchTrdAcc->Scale(1. / (Double_t) fEvent); + YPt_StsMuchTrdTofAcc->Scale(1. / (Double_t) fEvent); + + YPt_VtxReco->Scale(1. / (Double_t) fEvent); + YPt_VtxStsReco->Scale(1. / (Double_t) fEvent); + YPt_VtxStsMuchReco->Scale(1. / (Double_t) fEvent); + YPt_VtxStsMuchTrdReco->Scale(1. / (Double_t) fEvent); + YPt_VtxStsMuchTrdTofReco->Scale(1. / (Double_t) fEvent); + + YPt_pluto->Write(); + YPt_StsAcc->Write(); + YPt_StsMuchAcc->Write(); + YPt_StsMuchTrdAcc->Write(); + YPt_StsMuchTrdTofAcc->Write(); + + if (fAnnCut < 0) { + YPt_VtxReco->Write(); + YPt_VtxStsReco->Write(); + YPt_VtxStsMuchReco->Write(); + YPt_VtxStsMuchTrdReco->Write(); + } + YPt_VtxStsMuchTrdTofReco->Write(); + + f->Close(); + + if (fAnnCut > 0) + name = Form("eff_histo_ANN_%1.2f.root", fAnnCut); + else + name = "eff_histo.root"; + TFile* ff = new TFile(name, "recreate"); + + TDirectory* dir1 = ff->mkdir("mu Plus"); + dir1->cd(); + + TDirectory* dir1a = dir1->mkdir("accepted mu Plus"); + dir1a->cd(); + for (int j = 0; j < 4; j++) + acc_P[j][1]->Write(); + for (int j = 0; j < 4; j++) + acc_Theta[j][1]->Write(); + + dir1->cd(); + TDirectory* dir1b = dir1->mkdir("reconstructed mu Plus"); + dir1b->cd(); + if (fAnnCut < 0) { + for (int j = 0; j < 3; j++) + effReco_P[j][1]->Write(); + for (int j = 0; j < 3; j++) + effReco_Theta[j][1]->Write(); + for (int j = 0; j < 4; j++) + eff4pi_P[j][1]->Write(); + for (int j = 0; j < 4; j++) + eff4pi_Theta[j][1]->Write(); + } + effReco_P[3][1]->Write(); + effReco_Theta[3][1]->Write(); + eff4pi_P[4][1]->Write(); + eff4pi_Theta[4][1]->Write(); + + TDirectory* dir2 = ff->mkdir("mu Minus"); + dir2->cd(); + + TDirectory* dir2a = dir2->mkdir("accepted mu Minus"); + dir2a->cd(); + for (int j = 0; j < 4; j++) + acc_P[j][2]->Write(); + for (int j = 0; j < 4; j++) + acc_Theta[j][2]->Write(); + + dir2->cd(); + TDirectory* dir2b = dir2->mkdir("reconstructed mu Minus"); + dir2b->cd(); + if (fAnnCut < 0) { + for (int j = 0; j < 3; j++) + effReco_P[j][2]->Write(); + for (int j = 0; j < 3; j++) + effReco_Theta[j][2]->Write(); + for (int j = 0; j < 4; j++) + eff4pi_P[j][2]->Write(); + for (int j = 0; j < 4; j++) + eff4pi_Theta[j][2]->Write(); + } + effReco_P[3][2]->Write(); + effReco_Theta[3][2]->Write(); + eff4pi_P[4][2]->Write(); + eff4pi_Theta[4][2]->Write(); + + TDirectory* dir3 = ff->mkdir("signal"); + dir3->cd(); + + TDirectory* dir3a = dir3->mkdir("accepted signal"); + dir3a->cd(); + for (int j = 0; j < 4; j++) + acc_P[j][0]->Write(); + for (int j = 0; j < 4; j++) + acc_Theta[j][0]->Write(); + + dir3->cd(); + TDirectory* dir3b = dir3->mkdir("reconstructed signal"); + dir3b->cd(); + if (fAnnCut < 0) { + for (int j = 0; j < 3; j++) + effReco_P[j][0]->Write(); + for (int j = 0; j < 3; j++) + effReco_Theta[j][0]->Write(); + for (int j = 0; j < 4; j++) + eff4pi_P[j][0]->Write(); + for (int j = 0; j < 4; j++) + eff4pi_Theta[j][0]->Write(); + } + effReco_P[3][0]->Write(); + effReco_Theta[3][0]->Write(); + eff4pi_P[4][0]->Write(); + eff4pi_Theta[4][0]->Write(); + + ff->Close(); + // fPlutoFile->Close(); + } else { + if (fAnnCut > 0) + name = Form("sup_histo_ANN_%1.2f.root", fAnnCut); + else + name = "sup_histo.root"; + TFile* f = new TFile(name, "recreate"); + + if (fAnnCut < 0) + for (int j = 0; j < 5; j++) + BgSup[j]->Write(); + BgSup[5]->Write(); + f->Close(); } - - if(fAnnCut > 0)name = Form("YPtM_ANN_%1.2f.root", fAnnCut); - else name = "YPtM.root"; - TFile* f = new TFile(name,"recreate"); - - YPtM->Scale(1./(Double_t)fEvent); + + if (fAnnCut > 0) + name = Form("YPtM_ANN_%1.2f.root", fAnnCut); + else + name = "YPtM.root"; + TFile* f = new TFile(name, "recreate"); + + YPtM->Scale(1. / (Double_t) fEvent); YPtM->Write(); f->Close(); } diff --git a/analysis/PWGDIL/dimuon/CbmAnaDimuonAnalysis.h b/analysis/PWGDIL/dimuon/CbmAnaDimuonAnalysis.h index 50e9103375a879132a19a3773246d2e1dbd1edb1..fc8fd273ad0cb73f59e4906f35f7e6d32289fcc8 100644 --- a/analysis/PWGDIL/dimuon/CbmAnaDimuonAnalysis.h +++ b/analysis/PWGDIL/dimuon/CbmAnaDimuonAnalysis.h @@ -33,102 +33,116 @@ class TH3D; class TProfile; class TMultiLayerPerceptron; -class CbmAnaDimuonAnalysis : public FairTask{ +class CbmAnaDimuonAnalysis : public FairTask { public: CbmAnaDimuonAnalysis(TString name, TString setup); - - virtual ~CbmAnaDimuonAnalysis(){} + + virtual ~CbmAnaDimuonAnalysis() {} virtual InitStatus Init(); virtual void Exec(Option_t* opt); virtual void Finish(); virtual void SetParContainers(); - void SetChi2StsCut (Double_t cut) { fChi2StsCut = cut; } - void SetChi2MuchCut (Double_t cut) { fChi2MuchCut = cut; } - void SetChi2VertexCut(Double_t cut) { fChi2VertexCut = cut; } - - void SetNofMuchCut(Int_t cut) { fNofMuchCut = cut; } - void SetNofStsCut (Int_t cut) { fNofStsCut = cut; } - void SetNofTrdCut (Int_t cut) { fNofTrdCut = cut; } - void SetAnnCut (Double_t cut, Int_t neurons) { fAnnCut = cut; fNeurons = neurons;} - - void SetSigmaTofCut (Int_t cut) { fSigmaTofCut = cut; } - - void UseCuts(Bool_t cut) {fUseCuts = cut; } - void UseMC(Bool_t useMC) {fUseMC = useMC; } - -// void SetHistoFileName(TString name) {fFileName = name; } -// void SetEffFileName(TString name) {fEffFileName = name; } - + void SetChi2StsCut(Double_t cut) { fChi2StsCut = cut; } + void SetChi2MuchCut(Double_t cut) { fChi2MuchCut = cut; } + void SetChi2VertexCut(Double_t cut) { fChi2VertexCut = cut; } + + void SetNofMuchCut(Int_t cut) { fNofMuchCut = cut; } + void SetNofStsCut(Int_t cut) { fNofStsCut = cut; } + void SetNofTrdCut(Int_t cut) { fNofTrdCut = cut; } + void SetAnnCut(Double_t cut, Int_t neurons) { + fAnnCut = cut; + fNeurons = neurons; + } + + void SetSigmaTofCut(Int_t cut) { fSigmaTofCut = cut; } + + void UseCuts(Bool_t cut) { fUseCuts = cut; } + void UseMC(Bool_t useMC) { fUseMC = useMC; } + + // void SetHistoFileName(TString name) {fFileName = name; } + // void SetEffFileName(TString name) {fEffFileName = name; } + Double_t CalculateAnnValue(CbmAnaMuonCandidate* mu); - Double_t CalculateAnnValue(Double_t P, Double_t M, Double_t Chi2Vertex, Double_t Chi2STS, Double_t Chi2MUCH, Double_t Chi2TRD, Int_t NofSTS, Int_t NofMUCH, Int_t NofTRD, Int_t NofTOF); - - void FillProfile(TProfile *profile, Double_t param, Bool_t trigger); - + Double_t CalculateAnnValue(Double_t P, + Double_t M, + Double_t Chi2Vertex, + Double_t Chi2STS, + Double_t Chi2MUCH, + Double_t Chi2TRD, + Int_t NofSTS, + Int_t NofMUCH, + Int_t NofTRD, + Int_t NofTOF); + + void FillProfile(TProfile* profile, Double_t param, Bool_t trigger); + private: - Int_t fEvent; - TClonesArray* fMCTracks; - TClonesArray* fStsTracks; - TClonesArray* fStsTrackMatches; - TClonesArray* fMuchTracks; - TClonesArray* fMuchTrackMatches; - TClonesArray* fGlobalTracks; - TClonesArray* fTrdTracks; + Int_t fEvent; + TClonesArray* fMCTracks; + TClonesArray* fStsTracks; + TClonesArray* fStsTrackMatches; + TClonesArray* fMuchTracks; + TClonesArray* fMuchTrackMatches; + TClonesArray* fGlobalTracks; + TClonesArray* fTrdTracks; TClonesArray* fTofHit; - TClonesArray* fMuPlus; - TClonesArray* fMuMinus; + TClonesArray* fMuPlus; + TClonesArray* fMuMinus; TClonesArray* fParticles; TTree* fInputTree; TFile* fPlutoFile; - - CbmStsKFTrackFitter* fFitter; - CbmTrdTrackFitterKF* fFitterTRD; + + CbmStsKFTrackFitter* fFitter; + CbmTrdTrackFitterKF* fFitterTRD; CbmGlobalTrackFitterKF* fFitterGlobal; - + CbmVertex* fVertex; - Double_t fChi2StsCut; - Double_t fChi2MuchCut; + Double_t fChi2StsCut; + Double_t fChi2MuchCut; Double_t fChi2VertexCut; Double_t fAnnCut; Int_t fNeurons; Int_t fSigmaTofCut; - + Double_t fMass; - + Bool_t fUseCuts; Bool_t fUseMC; - - Int_t fNofMuchCut; - Int_t fNofStsCut; - Int_t fNofTrdCut; + + Int_t fNofMuchCut; + Int_t fNofStsCut; + Int_t fNofTrdCut; Double_t p0min, p1min, p2min; Double_t p0max, p1max, p2max; - + TString fFileAnnName; -// TString fEffFileName; -// TString fFileName; + // TString fEffFileName; + // TString fFileName; TString fPlutoFileName; TString fSetupName; - - CbmMuchGeoScheme* fGeoScheme; - - TH2D *YPt_pluto, *YPt_StsAcc, *YPt_StsMuchAcc, *YPt_StsMuchTrdAcc, *YPt_StsMuchTrdTofAcc; - TH2D *YPt_VtxReco, *YPt_VtxStsReco, *YPt_VtxStsMuchReco, *YPt_VtxStsMuchTrdReco, *YPt_VtxStsMuchTrdTofReco; - TH3D *YPtM; - - TProfile *acc_P[4][3], *acc_Theta[4][3]; + + CbmMuchGeoScheme* fGeoScheme; + + TH2D *YPt_pluto, *YPt_StsAcc, *YPt_StsMuchAcc, *YPt_StsMuchTrdAcc, + *YPt_StsMuchTrdTofAcc; + TH2D *YPt_VtxReco, *YPt_VtxStsReco, *YPt_VtxStsMuchReco, + *YPt_VtxStsMuchTrdReco, *YPt_VtxStsMuchTrdTofReco; + TH3D* YPtM; + + TProfile *acc_P[4][3], *acc_Theta[4][3]; TProfile *effReco_P[4][3], *effReco_Theta[4][3]; - TProfile *eff4pi_P[5][3], *eff4pi_Theta[5][3]; - - TH1D *BgSup[6]; - + TProfile *eff4pi_P[5][3], *eff4pi_Theta[5][3]; + + TH1D* BgSup[6]; + CbmAnaDimuonAnalysis(const CbmAnaDimuonAnalysis&); CbmAnaDimuonAnalysis operator=(const CbmAnaDimuonAnalysis&); - ClassDef(CbmAnaDimuonAnalysis,2); + ClassDef(CbmAnaDimuonAnalysis, 2); }; #endif diff --git a/analysis/PWGDIL/dimuon/CbmAnaMuonCandidate.cxx b/analysis/PWGDIL/dimuon/CbmAnaMuonCandidate.cxx index 88721399aceeab50c23ee4f89852681f09ca310e..e4061eb1761825dac6164d578b23551377318126 100644 --- a/analysis/PWGDIL/dimuon/CbmAnaMuonCandidate.cxx +++ b/analysis/PWGDIL/dimuon/CbmAnaMuonCandidate.cxx @@ -1,32 +1,28 @@ //---------------------------------------- // -// 2019 A. Senger a.senger@gsi.de +// 2019 A. Senger a.senger@gsi.de // //---------------------------------------- #include "CbmAnaMuonCandidate.h" CbmAnaMuonCandidate::CbmAnaMuonCandidate() - : TObject(), - track(), - fMom(), - fNStsHits(4), - fNMuchHits(4), - fNTrdHits(0), - fNTofHits(0), - fChiToVertex(100), - fSign(0), - fChiMuch(100), - fChiSts(100), - fChiTrd(100), - fChiGlobal(100), - fMu(0), - fPdg(13), - fId(-1), - fM(1e8) -{ - - -} + : TObject() + , track() + , fMom() + , fNStsHits(4) + , fNMuchHits(4) + , fNTrdHits(0) + , fNTofHits(0) + , fChiToVertex(100) + , fSign(0) + , fChiMuch(100) + , fChiSts(100) + , fChiTrd(100) + , fChiGlobal(100) + , fMu(0) + , fPdg(13) + , fId(-1) + , fM(1e8) {} ClassImp(CbmAnaMuonCandidate); diff --git a/analysis/PWGDIL/dimuon/CbmAnaMuonCandidate.h b/analysis/PWGDIL/dimuon/CbmAnaMuonCandidate.h index 620cf0aa2157fa423be89b67f3878b7c1e3d205b..2fdca2e81570107ef76a2d8f66e08df286805cbb 100644 --- a/analysis/PWGDIL/dimuon/CbmAnaMuonCandidate.h +++ b/analysis/PWGDIL/dimuon/CbmAnaMuonCandidate.h @@ -1,6 +1,6 @@ //---------------------------------------- // -// 2019 A. Senger a.senger@gsi.de +// 2019 A. Senger a.senger@gsi.de // //---------------------------------------- @@ -11,78 +11,78 @@ #include "TLorentzVector.h" #define NPLANES 31 -class CbmAnaMuonCandidate : public TObject{ - public: - CbmAnaMuonCandidate(); - virtual ~CbmAnaMuonCandidate(){}; - - void SetMomentum(TLorentzVector mom){fMom = TLorentzVector(mom);} - void SetSign(Int_t sign) {fSign = sign; } - - void SetNStsHits(Int_t nHits) {fNStsHits = nHits; } - void SetNMuchHits(Int_t nHits){fNMuchHits = nHits; } - void SetNTrdHits(Int_t nHits) {fNTrdHits = nHits; } - void SetNTofHits(Int_t nHits) {fNTofHits = nHits; } - - void SetChiToVertex(Double_t chi) {fChiToVertex = chi; } - - void SetChiMuch(Double_t chi) {fChiMuch = chi; } - void SetChiSts(Double_t chi) {fChiSts = chi; } - void SetChiTrd(Double_t chi) {fChiTrd = chi; } - void SetChiGlobal(Double_t chi) {fChiGlobal = chi; } - - void SetTrueMu(Int_t mu){fMu = mu;} - void SetStsPdg(Int_t pdg){fPdg = pdg;} - void SetAnnId(Double_t tID){fId = tID;} - - void SetTofM(Double_t m){fM = m;} - - TLorentzVector* GetMomentum() {return &fMom; } - Double_t GetSign() {return fSign;} - - Int_t GetNMuchHits(){ return fNMuchHits;} - Int_t GetNStsHits() { return fNStsHits; } - Int_t GetNTrdHits() { return fNTrdHits; } - Int_t GetNTofHits() { return fNTofHits; } - - Double_t GetChiToVertex() { return fChiToVertex; } - - Double_t GetChiMuch() { return fChiMuch; } - Double_t GetChiSts() { return fChiSts; } - Double_t GetChiTrd() { return fChiTrd; } - Double_t GetChiGlobal() { return fChiGlobal; } - - Int_t GetTrueMu(){ return fMu;} - Int_t GetStsPdg(){ return fPdg;} - Double_t GetAnnId(){ return fId;} - - Double_t GetTofM(){ return fM;} - - private: - CbmKFTrack track; - - TLorentzVector fMom; - - Int_t fNStsHits; - Int_t fNMuchHits; - Int_t fNTrdHits; - Int_t fNTofHits; - - Double_t fChiToVertex; - - Int_t fSign; - - Double_t fChiMuch; - Double_t fChiSts; - Double_t fChiTrd; - Double_t fChiGlobal; - - Int_t fMu; - Int_t fPdg; - Double_t fId; - Double_t fM; - - ClassDef(CbmAnaMuonCandidate,2); +class CbmAnaMuonCandidate : public TObject { +public: + CbmAnaMuonCandidate(); + virtual ~CbmAnaMuonCandidate() {}; + + void SetMomentum(TLorentzVector mom) { fMom = TLorentzVector(mom); } + void SetSign(Int_t sign) { fSign = sign; } + + void SetNStsHits(Int_t nHits) { fNStsHits = nHits; } + void SetNMuchHits(Int_t nHits) { fNMuchHits = nHits; } + void SetNTrdHits(Int_t nHits) { fNTrdHits = nHits; } + void SetNTofHits(Int_t nHits) { fNTofHits = nHits; } + + void SetChiToVertex(Double_t chi) { fChiToVertex = chi; } + + void SetChiMuch(Double_t chi) { fChiMuch = chi; } + void SetChiSts(Double_t chi) { fChiSts = chi; } + void SetChiTrd(Double_t chi) { fChiTrd = chi; } + void SetChiGlobal(Double_t chi) { fChiGlobal = chi; } + + void SetTrueMu(Int_t mu) { fMu = mu; } + void SetStsPdg(Int_t pdg) { fPdg = pdg; } + void SetAnnId(Double_t tID) { fId = tID; } + + void SetTofM(Double_t m) { fM = m; } + + TLorentzVector* GetMomentum() { return &fMom; } + Double_t GetSign() { return fSign; } + + Int_t GetNMuchHits() { return fNMuchHits; } + Int_t GetNStsHits() { return fNStsHits; } + Int_t GetNTrdHits() { return fNTrdHits; } + Int_t GetNTofHits() { return fNTofHits; } + + Double_t GetChiToVertex() { return fChiToVertex; } + + Double_t GetChiMuch() { return fChiMuch; } + Double_t GetChiSts() { return fChiSts; } + Double_t GetChiTrd() { return fChiTrd; } + Double_t GetChiGlobal() { return fChiGlobal; } + + Int_t GetTrueMu() { return fMu; } + Int_t GetStsPdg() { return fPdg; } + Double_t GetAnnId() { return fId; } + + Double_t GetTofM() { return fM; } + +private: + CbmKFTrack track; + + TLorentzVector fMom; + + Int_t fNStsHits; + Int_t fNMuchHits; + Int_t fNTrdHits; + Int_t fNTofHits; + + Double_t fChiToVertex; + + Int_t fSign; + + Double_t fChiMuch; + Double_t fChiSts; + Double_t fChiTrd; + Double_t fChiGlobal; + + Int_t fMu; + Int_t fPdg; + Double_t fId; + Double_t fM; + + ClassDef(CbmAnaMuonCandidate, 2); }; #endif diff --git a/macro/analysis/much/ANN.C b/macro/analysis/much/ANN.C index 4b07bd9502bfa85fd3cfa596f3e633a8e2fac6b0..0ff12c534259c5a508494b0781bbff2abf847885 100644 --- a/macro/analysis/much/ANN.C +++ b/macro/analysis/much/ANN.C @@ -8,102 +8,115 @@ // //--------------------------------------------------- -void ANN(Int_t energy=8, Int_t NofFiles=1000, - TString dir="../../much/data/", - TString setup="sis100_muon_lmvm", - Int_t ntrain=300, Int_t neurons=16) -{ - gStyle->SetCanvasColor(10); - gStyle->SetFrameFillColor(10); - gStyle->SetHistLineWidth(4); - gStyle->SetPadColor(10); - gStyle->SetStatColor(10); - gStyle->SetPalette(55); - gStyle->SetPaintTextFormat("2.1f"); - - Double_t Chi2Vertex, Chi2STS, Chi2MUCH, Chi2TRD; - Int_t NofSTS, NofMUCH, NofTRD; - Double_t P, M; - - TTree *Signal = new TTree("Signal","mu"); - Signal->Branch("P", &P, "P/D", 20000000); - Signal->Branch("M", &M, "M/D", 20000000); - Signal->Branch("Chi2Vertex",&Chi2Vertex,"Chi2Vertex/D",20000000); - Signal->Branch("Chi2STS", &Chi2STS, "Chi2STS/D", 20000000); - Signal->Branch("Chi2MUCH", &Chi2MUCH, "Chi2MUCH/D", 20000000); - Signal->Branch("Chi2TRD", &Chi2TRD, "Chi2TRD/D", 20000000); - Signal->Branch("NofSTS", &NofSTS, "NofSTS/I", 20000000); - Signal->Branch("NofMUCH",&NofMUCH,"NofMUCH/I",20000000); - Signal->Branch("NofTRD", &NofTRD, "NofTRD/I", 20000000); - - TTree *Bg = new TTree("Background","bg"); - Bg->Branch("P", &P, "P/D", 20000000); - Bg->Branch("M", &M, "M/D", 20000000); - Bg->Branch("Chi2Vertex",&Chi2Vertex,"Chi2Vertex/D",20000000); - Bg->Branch("Chi2STS", &Chi2STS, "Chi2STS/D", 20000000); - Bg->Branch("Chi2MUCH", &Chi2MUCH, "Chi2MUCH/D", 20000000); - Bg->Branch("Chi2TRD", &Chi2TRD, "Chi2TRD/D", 20000000); - Bg->Branch("NofSTS", &NofSTS, "NofSTS/I", 20000000); - Bg->Branch("NofMUCH",&NofMUCH,"NofMUCH/I",20000000); - Bg->Branch("NofTRD", &NofTRD, "NofTRD/I", 20000000); - - TString name; - - TClonesArray *MuPlus = new TClonesArray("CbmAnaMuonCandidate"); - TClonesArray *MuMinus = new TClonesArray("CbmAnaMuonCandidate"); - TTree *InputTree; - - Int_t MUCHhits= 9; - Int_t STShits = 4; - Int_t TRDhits = 1; - Int_t TOFhits = 1; - - Double_t MUCHchi2= 5.; - Double_t STSchi2 = 5.; - Double_t Vchi2 = 5.; - Double_t Mcut = 5.; - Double_t Pcut = 20.; - -//-------------------------------------------------------- - for(int i=0; i<2; i++){ - for(int k=1; k<NofFiles+1; k++) - { - if(i==0)name.Form("%s/%s/%dgev/omega/%d/muons.ana.root",dir.Data(),setup.Data(),energy,k); - else name.Form("%s/%s/%dgev/centr/%d/muons.ana.root",dir.Data(),setup.Data(),energy,k); - - TFile *f = new TFile(name); - if (f->IsZombie() || f->GetNkeys() < 1 || f->TestBit(TFile::kRecovered)){f->Close();continue;} - - if(k%100 == 0)cout<<"Input File "<<k<<" : "<<f->GetName()<<endl; - - InputTree = (TTree*)f->Get("cbmsim"); - InputTree->SetBranchAddress("MuPlus", &MuPlus); - InputTree->SetBranchAddress("MuMinus", &MuMinus); - - //----------------------------------------------------------- - - for(int iEvent=0; iEvent<InputTree->GetEntries(); iEvent++) - { - InputTree->GetEntry(iEvent); - - int NofPlus = MuPlus->GetEntriesFast(); - int NofMinus = MuMinus->GetEntriesFast(); - for (int iPart=0; iPart<NofPlus; iPart++){ - Bool_t trigger = kFALSE; - - CbmAnaMuonCandidate *mu_pl = (CbmAnaMuonCandidate*)MuPlus->At(iPart); - - if((i==0 && mu_pl->GetTrueMu()==1) || i==1)trigger=kTRUE; - if(!trigger)continue; +void ANN(Int_t energy = 8, + Int_t NofFiles = 1000, + TString dir = "../../much/data/", + TString setup = "sis100_muon_lmvm", + Int_t ntrain = 300, + Int_t neurons = 16) { + gStyle->SetCanvasColor(10); + gStyle->SetFrameFillColor(10); + gStyle->SetHistLineWidth(4); + gStyle->SetPadColor(10); + gStyle->SetStatColor(10); + gStyle->SetPalette(55); + gStyle->SetPaintTextFormat("2.1f"); - TLorentzVector *P_pl = mu_pl->GetMomentum(); - M = mu_pl->GetTofM(); - P = P_pl->P(); - - if(mu_pl->GetNStsHits() < STShits)continue; - if(mu_pl->GetNMuchHits()< MUCHhits)continue; - if(mu_pl->GetNTofHits() < 1)continue; - /* + Double_t Chi2Vertex, Chi2STS, Chi2MUCH, Chi2TRD; + Int_t NofSTS, NofMUCH, NofTRD; + Double_t P, M; + + TTree* Signal = new TTree("Signal", "mu"); + Signal->Branch("P", &P, "P/D", 20000000); + Signal->Branch("M", &M, "M/D", 20000000); + Signal->Branch("Chi2Vertex", &Chi2Vertex, "Chi2Vertex/D", 20000000); + Signal->Branch("Chi2STS", &Chi2STS, "Chi2STS/D", 20000000); + Signal->Branch("Chi2MUCH", &Chi2MUCH, "Chi2MUCH/D", 20000000); + Signal->Branch("Chi2TRD", &Chi2TRD, "Chi2TRD/D", 20000000); + Signal->Branch("NofSTS", &NofSTS, "NofSTS/I", 20000000); + Signal->Branch("NofMUCH", &NofMUCH, "NofMUCH/I", 20000000); + Signal->Branch("NofTRD", &NofTRD, "NofTRD/I", 20000000); + + TTree* Bg = new TTree("Background", "bg"); + Bg->Branch("P", &P, "P/D", 20000000); + Bg->Branch("M", &M, "M/D", 20000000); + Bg->Branch("Chi2Vertex", &Chi2Vertex, "Chi2Vertex/D", 20000000); + Bg->Branch("Chi2STS", &Chi2STS, "Chi2STS/D", 20000000); + Bg->Branch("Chi2MUCH", &Chi2MUCH, "Chi2MUCH/D", 20000000); + Bg->Branch("Chi2TRD", &Chi2TRD, "Chi2TRD/D", 20000000); + Bg->Branch("NofSTS", &NofSTS, "NofSTS/I", 20000000); + Bg->Branch("NofMUCH", &NofMUCH, "NofMUCH/I", 20000000); + Bg->Branch("NofTRD", &NofTRD, "NofTRD/I", 20000000); + + TString name; + + TClonesArray* MuPlus = new TClonesArray("CbmAnaMuonCandidate"); + TClonesArray* MuMinus = new TClonesArray("CbmAnaMuonCandidate"); + TTree* InputTree; + + Int_t MUCHhits = 9; + Int_t STShits = 4; + Int_t TRDhits = 1; + Int_t TOFhits = 1; + + Double_t MUCHchi2 = 5.; + Double_t STSchi2 = 5.; + Double_t Vchi2 = 5.; + Double_t Mcut = 5.; + Double_t Pcut = 20.; + + //-------------------------------------------------------- + for (int i = 0; i < 2; i++) { + for (int k = 1; k < NofFiles + 1; k++) { + if (i == 0) + name.Form("%s/%s/%dgev/omega/%d/muons.ana.root", + dir.Data(), + setup.Data(), + energy, + k); + else + name.Form("%s/%s/%dgev/centr/%d/muons.ana.root", + dir.Data(), + setup.Data(), + energy, + k); + + TFile* f = new TFile(name); + if (f->IsZombie() || f->GetNkeys() < 1 || f->TestBit(TFile::kRecovered)) { + f->Close(); + continue; + } + + if (k % 100 == 0) + cout << "Input File " << k << " : " << f->GetName() << endl; + + InputTree = (TTree*) f->Get("cbmsim"); + InputTree->SetBranchAddress("MuPlus", &MuPlus); + InputTree->SetBranchAddress("MuMinus", &MuMinus); + + //----------------------------------------------------------- + + for (int iEvent = 0; iEvent < InputTree->GetEntries(); iEvent++) { + InputTree->GetEntry(iEvent); + + int NofPlus = MuPlus->GetEntriesFast(); + int NofMinus = MuMinus->GetEntriesFast(); + for (int iPart = 0; iPart < NofPlus; iPart++) { + Bool_t trigger = kFALSE; + + CbmAnaMuonCandidate* mu_pl = (CbmAnaMuonCandidate*) MuPlus->At(iPart); + + if ((i == 0 && mu_pl->GetTrueMu() == 1) || i == 1) trigger = kTRUE; + if (!trigger) continue; + + TLorentzVector* P_pl = mu_pl->GetMomentum(); + M = mu_pl->GetTofM(); + P = P_pl->P(); + + if (mu_pl->GetNStsHits() < STShits) continue; + if (mu_pl->GetNMuchHits() < MUCHhits) continue; + if (mu_pl->GetNTofHits() < 1) continue; + /* if(M > Mcut || M < -10)continue; if(P > Pcut)continue; @@ -114,34 +127,37 @@ void ANN(Int_t energy=8, Int_t NofFiles=1000, if(mu_pl->GetNMuchHits()< 1)continue; if(M > Mcut)continue; */ - Chi2Vertex = mu_pl->GetChiToVertex(); - Chi2STS = mu_pl->GetChiSts(); - Chi2MUCH = mu_pl->GetChiMuch(); - Chi2TRD = mu_pl->GetChiTrd(); - - NofSTS = mu_pl->GetNStsHits(); - NofMUCH = mu_pl->GetNMuchHits(); - NofTRD = mu_pl->GetNTrdHits(); - if(i==0)Signal->Fill(); else Bg->Fill(); - - } - for (int jPart=0; jPart<NofMinus; jPart++){ - Bool_t trigger = kFALSE; - - CbmAnaMuonCandidate *mu_mn = (CbmAnaMuonCandidate*)MuMinus->At(jPart); - - if((i==0 && mu_mn->GetTrueMu()==1) || i==1)trigger=kTRUE; - if(!trigger)continue; + Chi2Vertex = mu_pl->GetChiToVertex(); + Chi2STS = mu_pl->GetChiSts(); + Chi2MUCH = mu_pl->GetChiMuch(); + Chi2TRD = mu_pl->GetChiTrd(); - TLorentzVector *P_mn = mu_mn->GetMomentum(); - M = mu_mn->GetTofM(); - P = P_mn->P(); - - if(mu_mn->GetNStsHits() < STShits)continue; - if(mu_mn->GetNMuchHits()< MUCHhits)continue; - if(mu_mn->GetNTofHits() < 1)continue; - - /* + NofSTS = mu_pl->GetNStsHits(); + NofMUCH = mu_pl->GetNMuchHits(); + NofTRD = mu_pl->GetNTrdHits(); + if (i == 0) + Signal->Fill(); + else + Bg->Fill(); + } + for (int jPart = 0; jPart < NofMinus; jPart++) { + Bool_t trigger = kFALSE; + + CbmAnaMuonCandidate* mu_mn = + (CbmAnaMuonCandidate*) MuMinus->At(jPart); + + if ((i == 0 && mu_mn->GetTrueMu() == 1) || i == 1) trigger = kTRUE; + if (!trigger) continue; + + TLorentzVector* P_mn = mu_mn->GetMomentum(); + M = mu_mn->GetTofM(); + P = P_mn->P(); + + if (mu_mn->GetNStsHits() < STShits) continue; + if (mu_mn->GetNMuchHits() < MUCHhits) continue; + if (mu_mn->GetNTofHits() < 1) continue; + + /* if(M > Mcut || M < -10)continue; if(P > Pcut)continue; @@ -152,125 +168,131 @@ void ANN(Int_t energy=8, Int_t NofFiles=1000, if(mu_mn->GetNMuchHits()< 1)continue; if(M > Mcut)continue; */ - Chi2Vertex = mu_mn->GetChiToVertex(); - Chi2STS = mu_mn->GetChiSts(); - Chi2MUCH = mu_mn->GetChiMuch(); - Chi2TRD = mu_mn->GetChiTrd(); - - NofSTS = mu_mn->GetNStsHits(); - NofMUCH = mu_mn->GetNMuchHits(); - NofTRD = mu_mn->GetNTrdHits(); - if(i==0)Signal->Fill(); else Bg->Fill(); - } - } - f->Close(); - }} - - Int_t type; - - TTree *simu = new TTree("MonteCarlo", "Filtered Monte Carlo Events"); - simu->Branch("P", &P, "P/D", 20000000); - simu->Branch("M", &M, "M/D", 20000000); - simu->Branch("Chi2Vertex",&Chi2Vertex,"Chi2Vertex/D", 20000000); - simu->Branch("Chi2STS", &Chi2STS, "Chi2STS/D", 20000000); - simu->Branch("Chi2MUCH", &Chi2MUCH, "Chi2MUCH/D", 20000000); - simu->Branch("Chi2TRD", &Chi2TRD, "Chi2TRD/D", 20000000); - simu->Branch("NofSTS", &NofSTS, "NofSTS/I", 20000000); - simu->Branch("NofMUCH",&NofMUCH,"NofMUCH/I", 20000000); - simu->Branch("NofTRD", &NofTRD, "NofTRD/I", 20000000); - simu->Branch("type", &type, "type/I", 20000000); - - type=1; // 1 = signal, 0 = background - Int_t i; - for (i = 0; i < Signal->GetEntries(); i++) { - Signal->GetEntry(i); - simu->Fill(); - } - - type = 0; - for (i = 0; i < Bg->GetEntries(); i++) { - Bg->GetEntry(i); - simu->Fill(); - } - - name.Form("@P,@M,@Chi2Vertex,@Chi2STS,@Chi2MUCH,@Chi2TRD,@NofSTS,@NofMUCH,@NofTRD:%d:type",neurons); - TMultiLayerPerceptron *mlp = new TMultiLayerPerceptron(name.Data(),simu,"Entry$%2"); - - mlp->Train(ntrain, "text,graph,update=10"); - name.Form("muid_ann_%d_%s_weights.txt",neurons,setup.Data()); - mlp->DumpWeights(name); - - TCanvas* mlpa_canvas = new TCanvas("mlpa_canvas","Network analysis"); - mlpa_canvas->Divide(2,2); - TMLPAnalyzer ana(mlp); - ana.GatherInformations(); - ana.CheckNetwork(); - mlpa_canvas->cd(1); - ana.DrawDInputs(); - mlpa_canvas->cd(2); - mlp->Draw(); - mlpa_canvas->cd(3); - ana.DrawNetwork(0,"type==1","type==0"); - mlpa_canvas->cd(4); - - TH1F *bg = new TH1F("bgh", "NN output", 50, -.5, 1.5); - TH1F *sig = new TH1F("sigh", "NN output", 50, -.5, 1.5); - bg->SetDirectory(0); - sig->SetDirectory(0); - - Double_t params[9]; - - for (i = 0; i < Bg->GetEntries(); i++) { - Bg->GetEntry(i); - params[0] = P; - params[1] = M; - params[2] = Chi2Vertex; - params[3] = Chi2STS; - params[4] = Chi2MUCH; - params[5] = Chi2TRD; - - params[6] = (Double_t)NofSTS; - params[7] = (Double_t)NofMUCH; - params[8] = (Double_t)NofTRD; - - bg->Fill(mlp->Evaluate(0, params)); - } - - for (i = 0; i < Signal->GetEntries(); i++) { - Signal->GetEntry(i); - params[0] = P; - params[1] = M; - params[2] = Chi2Vertex; - params[3] = Chi2STS; - params[4] = Chi2MUCH; - params[5] = Chi2TRD; - - params[6] = (Double_t)NofSTS; - params[7] = (Double_t)NofMUCH; - params[8] = (Double_t)NofTRD; - - sig->Fill(mlp->Evaluate(0,params)); - } - - bg->SetLineColor(kBlue); - bg->SetFillStyle(3008); bg->SetFillColor(kBlue); - sig->SetLineColor(kRed); - sig->SetFillStyle(3003); sig->SetFillColor(kRed); - bg->SetStats(0); - sig->SetStats(0); - bg->Draw(); - sig->Draw("same"); - TLegend *legend = new TLegend(.75, .80, .95, .95); - legend->AddEntry(bg, "Background"); - legend->AddEntry(sig, "Signal (#mu)"); - legend->Draw(); - mlpa_canvas->cd(0); - - name.Form("ann_%d_%s_canvas.root",neurons,setup.Data()); - TFile *FFF = new TFile(name,"recreate"); - mlpa_canvas->Write(); - FFF->Close(); -} + Chi2Vertex = mu_mn->GetChiToVertex(); + Chi2STS = mu_mn->GetChiSts(); + Chi2MUCH = mu_mn->GetChiMuch(); + Chi2TRD = mu_mn->GetChiTrd(); + + NofSTS = mu_mn->GetNStsHits(); + NofMUCH = mu_mn->GetNMuchHits(); + NofTRD = mu_mn->GetNTrdHits(); + if (i == 0) + Signal->Fill(); + else + Bg->Fill(); + } + } + f->Close(); + } + } + + Int_t type; + + TTree* simu = new TTree("MonteCarlo", "Filtered Monte Carlo Events"); + simu->Branch("P", &P, "P/D", 20000000); + simu->Branch("M", &M, "M/D", 20000000); + simu->Branch("Chi2Vertex", &Chi2Vertex, "Chi2Vertex/D", 20000000); + simu->Branch("Chi2STS", &Chi2STS, "Chi2STS/D", 20000000); + simu->Branch("Chi2MUCH", &Chi2MUCH, "Chi2MUCH/D", 20000000); + simu->Branch("Chi2TRD", &Chi2TRD, "Chi2TRD/D", 20000000); + simu->Branch("NofSTS", &NofSTS, "NofSTS/I", 20000000); + simu->Branch("NofMUCH", &NofMUCH, "NofMUCH/I", 20000000); + simu->Branch("NofTRD", &NofTRD, "NofTRD/I", 20000000); + simu->Branch("type", &type, "type/I", 20000000); + + type = 1; // 1 = signal, 0 = background + Int_t i; + for (i = 0; i < Signal->GetEntries(); i++) { + Signal->GetEntry(i); + simu->Fill(); + } + type = 0; + for (i = 0; i < Bg->GetEntries(); i++) { + Bg->GetEntry(i); + simu->Fill(); + } + name.Form("@P,@M,@Chi2Vertex,@Chi2STS,@Chi2MUCH,@Chi2TRD,@NofSTS,@NofMUCH,@" + "NofTRD:%d:type", + neurons); + TMultiLayerPerceptron* mlp = + new TMultiLayerPerceptron(name.Data(), simu, "Entry$%2"); + mlp->Train(ntrain, "text,graph,update=10"); + name.Form("muid_ann_%d_%s_weights.txt", neurons, setup.Data()); + mlp->DumpWeights(name); + + TCanvas* mlpa_canvas = new TCanvas("mlpa_canvas", "Network analysis"); + mlpa_canvas->Divide(2, 2); + TMLPAnalyzer ana(mlp); + ana.GatherInformations(); + ana.CheckNetwork(); + mlpa_canvas->cd(1); + ana.DrawDInputs(); + mlpa_canvas->cd(2); + mlp->Draw(); + mlpa_canvas->cd(3); + ana.DrawNetwork(0, "type==1", "type==0"); + mlpa_canvas->cd(4); + + TH1F* bg = new TH1F("bgh", "NN output", 50, -.5, 1.5); + TH1F* sig = new TH1F("sigh", "NN output", 50, -.5, 1.5); + bg->SetDirectory(0); + sig->SetDirectory(0); + + Double_t params[9]; + + for (i = 0; i < Bg->GetEntries(); i++) { + Bg->GetEntry(i); + params[0] = P; + params[1] = M; + params[2] = Chi2Vertex; + params[3] = Chi2STS; + params[4] = Chi2MUCH; + params[5] = Chi2TRD; + + params[6] = (Double_t) NofSTS; + params[7] = (Double_t) NofMUCH; + params[8] = (Double_t) NofTRD; + + bg->Fill(mlp->Evaluate(0, params)); + } + + for (i = 0; i < Signal->GetEntries(); i++) { + Signal->GetEntry(i); + params[0] = P; + params[1] = M; + params[2] = Chi2Vertex; + params[3] = Chi2STS; + params[4] = Chi2MUCH; + params[5] = Chi2TRD; + + params[6] = (Double_t) NofSTS; + params[7] = (Double_t) NofMUCH; + params[8] = (Double_t) NofTRD; + + sig->Fill(mlp->Evaluate(0, params)); + } + + bg->SetLineColor(kBlue); + bg->SetFillStyle(3008); + bg->SetFillColor(kBlue); + sig->SetLineColor(kRed); + sig->SetFillStyle(3003); + sig->SetFillColor(kRed); + bg->SetStats(0); + sig->SetStats(0); + bg->Draw(); + sig->Draw("same"); + TLegend* legend = new TLegend(.75, .80, .95, .95); + legend->AddEntry(bg, "Background"); + legend->AddEntry(sig, "Signal (#mu)"); + legend->Draw(); + mlpa_canvas->cd(0); + + name.Form("ann_%d_%s_canvas.root", neurons, setup.Data()); + TFile* FFF = new TFile(name, "recreate"); + mlpa_canvas->Write(); + FFF->Close(); +} diff --git a/macro/much/run_ana.C b/macro/much/run_ana.C index c2e8f0b87b84709f8a26ae1ebf6dc7dacc1f0145..d85a0bdb83df581adcf83091f0f9b9228e4ec264 100644 --- a/macro/much/run_ana.C +++ b/macro/much/run_ana.C @@ -4,32 +4,38 @@ // //--------------------------------------------------- -void run_ana(Int_t nEvents=1000, TString dataSet = "muons", - TString setup = "sis100_muon_lmvm", Bool_t useMC = kTRUE, TString pluto = "", Double_t ANN=-1) -{ - TString dir = ""; +void run_ana(Int_t nEvents = 1000, + TString dataSet = "muons", + TString setup = "sis100_muon_lmvm", + Bool_t useMC = kTRUE, + TString pluto = "", + Double_t ANN = -1) { + TString dir = ""; TString traFile = dir + dataSet + ".tra.root"; TString parFile = dir + dataSet + ".par.root"; TString recoFile = dir + dataSet + ".rec.root"; TString outFile; - if(ANN<0)outFile = dataSet + ".ana.root"; - else outFile = Form("%s.ana.ANN_%1.2f.root",dataSet.Data(), ANN); + if (ANN < 0) + outFile = dataSet + ".ana.root"; + else + outFile = Form("%s.ana.ANN_%1.2f.root", dataSet.Data(), ANN); - FairRunAna* run = new FairRunAna(); - run->SetInputFile(recoFile); - run->AddFriend(traFile); - run->SetOutputFile(outFile); - // run->SetGenerateRunInfo(kTRUE); - - // ----- Load the geometry setup ------------------------------------- + FairRunAna* run = new FairRunAna(); + run->SetInputFile(recoFile); + run->AddFriend(traFile); + run->SetOutputFile(outFile); + // run->SetGenerateRunInfo(kTRUE); + + // ----- Load the geometry setup ------------------------------------- // ----- Environment -------------------------------------------------- TString srcDir = gSystem->Getenv("VMCWORKDIR"); // top source directory // ------------------------------------------------------------------------ std::cout << std::endl; - TString setupFile = srcDir + "/geometry/setup/setup_" + setup + ".C"; + TString setupFile = srcDir + "/geometry/setup/setup_" + setup + ".C"; TString setupFunct = "setup_"; - setupFunct = setupFunct + setup + "()"; - std::cout << "-I- " << ": Loading macro " << setupFile << std::endl; + setupFunct = setupFunct + setup + "()"; + std::cout << "-I- " + << ": Loading macro " << setupFile << std::endl; gROOT->LoadMacro(setupFile); gROOT->ProcessLine(setupFunct); // You can modify the pre-defined setup by using @@ -43,7 +49,7 @@ void run_ana(Int_t nEvents=1000, TString dataSet = "muons", // ------------------------------------------------------------------------ CbmKF* kf = new CbmKF(); run->AddTask(kf); - + /* (VF) Not much sense in running L1 witout STS local reconstruction CbmL1* L1 = new CbmL1(); TString stsGeoTag; @@ -56,8 +62,8 @@ void run_ana(Int_t nEvents=1000, TString dataSet = "muons", } run->AddTask(L1); */ - - CbmAnaDimuonAnalysis* ana = new CbmAnaDimuonAnalysis(pluto,setup); + + CbmAnaDimuonAnalysis* ana = new CbmAnaDimuonAnalysis(pluto, setup); /* ana->SetChi2MuchCut(3.); ana->SetChi2StsCut(2.); @@ -73,8 +79,8 @@ void run_ana(Int_t nEvents=1000, TString dataSet = "muons", ana->UseMC(useMC); run->AddTask(ana); - // ----- Parameter database -------------------------------------------- - FairRuntimeDb* rtdb = run->GetRuntimeDb(); + // ----- Parameter database -------------------------------------------- + FairRuntimeDb* rtdb = run->GetRuntimeDb(); FairParRootFileIo* parIo1 = new FairParRootFileIo(); //FairParAsciiFileIo* parIo2 = new FairParAsciiFileIo(); parIo1->open(parFile.Data()); @@ -85,13 +91,13 @@ void run_ana(Int_t nEvents=1000, TString dataSet = "muons", rtdb->saveOutput(); // ------------------------------------------------------------------------ - // ----- Initialize and run -------------------------------------------- - run->Init(); - run->Run(0, nEvents); - // ------------------------------------------------------------------------ + // ----- Initialize and run -------------------------------------------- + run->Init(); + run->Run(0, nEvents); + // ------------------------------------------------------------------------ - cout << " Test passed" << endl; - cout << " All ok " << endl; + cout << " Test passed" << endl; + cout << " All ok " << endl; - // RemoveGeoManager(); + // RemoveGeoManager(); }