/** * \file CbmMuonAnalysis.cxx **/ #include "CbmMuonAnalysis.h" #include "FairRootManager.h" #include "CbmGlobalTrack.h" #include "CbmStsTrack.h" #include "CbmMuchTrack.h" #include "TCanvas.h" #include "TMath.h" #include "TLorentzVector.h" #include using std::cout; using std::endl; CbmMuonAnalysis::CbmMuonAnalysis(): FairTask("CbmMuonAnalysis"), fGlobalTracks(NULL), fStsTracks(NULL), fMuchTracks(NULL), fPrimVertex(NULL), fhStsParamFirstMom(NULL), fhMuchParamLastMom(NULL), fhMuchNofHitsInTrack(NULL), fhMuchParamLastZ(NULL), fhChiSqPrimary(NULL), fhInvMass(NULL), fKFFitter() { } CbmMuonAnalysis::~CbmMuonAnalysis() { } InitStatus CbmMuonAnalysis::Init() { ReadDataBranches(); CreateHistograms(); // Initialize track fit tool fKFFitter.Init(); return kSUCCESS; } void CbmMuonAnalysis::Exec( Option_t* option) { Int_t nofGlobalTracks = fGlobalTracks->GetEntriesFast(); // Loop over reconstructed global tracks. for (Int_t iGlobalTrack = 0; iGlobalTrack < nofGlobalTracks; iGlobalTrack++) { // Declare const pointer to CbmGlobalTrack. // Use static_cast to explicitly cast the object returned by TClonesArray::At(). const CbmGlobalTrack* globalTrack = static_cast(fGlobalTracks->At(iGlobalTrack)); // Get index of CbmStsTrack from CbmGlobalTrack. Int_t stsId = globalTrack->GetStsTrackIndex(); if (stsId != -1) { // Declare pointer to CbmStsTrack. // Use C-style type casting with braces (). CbmStsTrack* stsTrack = (CbmStsTrack*)fStsTracks->At(stsId); // Take the first STS track parameters. FairTrackParam* stsParamFirst = stsTrack->GetParamFirst(); // Fill the histogram with inverse Q/p parameter. fhStsParamFirstMom->Fill(1. / stsParamFirst->GetQp()); // Refit the track with muon mass fKFFitter.DoFit(stsTrack, 13); // Get chi squared vertex Double_t chi2Prim = fKFFitter.GetChiToVertex(stsTrack, fPrimVertex); fhChiSqPrimary->Fill(chi2Prim); } // GetIndex of CbmMuchTrack from CbmGlobalTrack. Int_t muchId = globalTrack->GetMuchTrackIndex(); if (muchId != -1) { const CbmMuchTrack* muchTrack = (const CbmMuchTrack*)fMuchTracks->At(muchId); // Take the last MUCH track parameter. const FairTrackParam* muchParamLast = muchTrack->GetParamLast(); // Fill the histogram with inverse Q/p parameter. fhMuchParamLastMom->Fill(1. / muchParamLast->GetQp()); // Fill histogram with the last Z position fhMuchParamLastZ->Fill(muchParamLast->GetZ()); // fill histogram with number of hits fhMuchNofHitsInTrack->Fill(muchTrack->GetNofHits()); } } // Invariant mass histogram for (Int_t igt = 0; igt < nofGlobalTracks - 1; igt++) { const CbmGlobalTrack* globalTrack1 = static_cast(fGlobalTracks->At(igt)); if (!IsMuon(globalTrack1)) continue; for (Int_t jgt = igt + 1; jgt < nofGlobalTracks; jgt++) { const CbmGlobalTrack* globalTrack2 = static_cast(fGlobalTracks->At(jgt)); if (!IsMuon(globalTrack2)) continue; Int_t stsId1 = globalTrack1->GetStsTrackIndex(); Int_t stsId2 = globalTrack2->GetStsTrackIndex(); if (stsId1 < 0 || stsId2 < 0) continue; CbmStsTrack* stsTrack1 = static_cast(fStsTracks->At(stsId1)); CbmStsTrack* stsTrack2 = static_cast(fStsTracks->At(stsId2)); // Fit tracks to the primary vertex FairTrackParam vtxParam1, vtxParam2; fKFFitter.FitToVertex(stsTrack1, fPrimVertex, &vtxParam1); fKFFitter.FitToVertex(stsTrack2, fPrimVertex, &vtxParam2); // check that charges are oposite if ( (vtxParam1.GetQp() > 0 && vtxParam2.GetQp() > 0) || (vtxParam1.GetQp() < 0 && vtxParam2.GetQp() < 0)) continue; TVector3 mom1, mom2; vtxParam1.Momentum(mom1); vtxParam2.Momentum(mom2); fhInvMass->Fill(CalculateInvMass(mom1, mom2)); } } static Int_t nofEvents = 0; cout << "-I- CbmMuonAnalysis::Exec: event: " << nofEvents++ << endl; } void CbmMuonAnalysis::Finish() { DrawHistograms(); WriteHistograms(); } void CbmMuonAnalysis::SetParContainers() { } void CbmMuonAnalysis::ReadDataBranches() { FairRootManager* ioman = FairRootManager::Instance(); if (NULL == ioman) { Fatal("Init","FairRootManager is not instantiated"); } fGlobalTracks = (TClonesArray*) ioman->GetObject("GlobalTrack"); if (NULL == fGlobalTracks) { Fatal("Init","No GlobalTrack array!"); } fStsTracks = (TClonesArray*) ioman->GetObject("StsTrack"); if (NULL == fStsTracks) { Fatal("Init","No StsTrack array!"); } fMuchTracks = (TClonesArray*) ioman->GetObject("MuchTrack"); if (NULL == fMuchTracks) { Fatal("Init","No MuchTrack array!"); } fPrimVertex = (CbmVertex*) ioman->GetObject("PrimaryVertex"); if (NULL == fPrimVertex) { Fatal("Init","No PrimaryVertex array!"); } } void CbmMuonAnalysis::CreateHistograms() { fhStsParamFirstMom = new TH1D("hStsParamFirstMom", "hStsParamFirstMom;P [GeV/c]; Entries", 100, -25., 25.); fhMuchParamLastMom = new TH1D("hMuchParamLastMom", "hMuchParamLastMom;P [GeV/c]; Entries", 100, -25., 25.); fhMuchNofHitsInTrack= new TH1D("hMuchNofHitsInTrack", "hMuchNofHitsInTrack;Number of hits; Entries", 25, 0., 25.); fhMuchParamLastZ = new TH1D("hMuchParamLastZ", "hMuchParamLastZ;Z [cm]; Entries", 600, 0, 600.); fhChiSqPrimary = new TH1D("hChiSqPrimary", "hChiSqPrimary;#chi^{2}_{primary}; Entries", 100, 0, 10.); fhInvMass = new TH1D("hInvMass", "hInvMass;M_{#mu+#mu-} [GeV/c]; Entries", 400, 0, 4.); } void CbmMuonAnalysis::DrawHistograms() { TCanvas* canvas = new TCanvas("muon_analysis_mom_distribution", "muon_analysis_mom_distribution", 1200, 600); canvas->Divide(2, 1); canvas->cd(1); fhStsParamFirstMom->Draw(); canvas->cd(2); fhMuchParamLastMom->Draw(); TCanvas* canvas2 = new TCanvas("muon_analysis_other_distributions", "muon_analysis_other_distributions", 1200, 1200); canvas2->Divide(2, 2); canvas2->cd(1); fhMuchNofHitsInTrack->Draw(); canvas2->cd(2); fhMuchParamLastZ->Draw(); canvas2->cd(3); fhChiSqPrimary->Draw(); TCanvas* canvas3 = new TCanvas("muon_analysis_inv_mass", "muon_analysis_inv_mass", 600, 600); canvas3->cd(1); fhInvMass->Draw(); } void CbmMuonAnalysis::WriteHistograms() { fhStsParamFirstMom->Write(); fhMuchParamLastMom->Write(); fhMuchNofHitsInTrack->Write(); fhMuchParamLastZ->Write(); fhChiSqPrimary->Write(); fhInvMass->Write(); } Bool_t CbmMuonAnalysis::IsMuon( const CbmGlobalTrack* track) { // First check STS track Bool_t stsOK = false; Int_t stsId = track->GetStsTrackIndex(); if (stsId >= 0) { CbmStsTrack* stsTrack = static_cast(fStsTracks->At(stsId)); // Refit the track with muon mass fKFFitter.DoFit(stsTrack, 13); // Get chi squared vertex Double_t chi2Prim = fKFFitter.GetChiToVertex(stsTrack, fPrimVertex); stsOK = chi2Prim < 3; // Cut on primary vertex } else { stsOK = false; // if no STS track in global track } // Check MUCH track Bool_t muchOK = false; Int_t muchId = track->GetMuchTrackIndex(); if (muchId >= 0) { const CbmMuchTrack* muchTrack = static_cast(fMuchTracks->At(muchId)); muchOK = muchTrack->GetNofHits() > 15; // Cut on number of hits in MUCH track } else { muchOK = false; } return stsOK && muchOK; } Double_t CbmMuonAnalysis::CalculateInvMass( const TVector3& momP, const TVector3& momN) { static const Double_t M2MU = 0.1057 * 0.1057; // squared muon mass Double_t energyP = TMath::Sqrt(momP.Mag2() + M2MU); TLorentzVector lorVecP(momP, energyP); Double_t energyN = TMath::Sqrt(momN.Mag2() + M2MU); TLorentzVector lorVecN(momN, energyN); Double_t anglePair = lorVecN.Angle(lorVecP.Vect()); return 2. * TMath::Sin(anglePair / 2.) * TMath::Sqrt(momN.Mag() * momP.Mag()); } ClassImp(CbmMuonAnalysis)