From c0cccc26b600afb57482a35f3334d6090c6fcf95 Mon Sep 17 00:00:00 2001 From: auras Date: Tue, 16 Dec 2014 17:44:36 +0100 Subject: [PATCH] MFT analysis class --- PWG/mftmuondep/CMakeLists.txt | 2 +- .../AliAnalysisTaskDimuonBackground.cxx | 525 ++++++++++++++++++ .../AliAnalysisTaskDimuonBackground.h | 80 +++ PWG/mftmuondep/charmonium/CMakeLists.txt | 41 +- .../PWGmftmuondepCharmoniumLinkDef.h | 11 + .../RunAnalysisTaskDimuonBackground.C | 297 ++++++++++ 6 files changed, 919 insertions(+), 37 deletions(-) create mode 100644 PWG/mftmuondep/charmonium/AliAnalysisTaskDimuonBackground.cxx create mode 100755 PWG/mftmuondep/charmonium/AliAnalysisTaskDimuonBackground.h create mode 100644 PWG/mftmuondep/charmonium/PWGmftmuondepCharmoniumLinkDef.h create mode 100755 PWG/mftmuondep/charmonium/RunAnalysisTaskDimuonBackground.C diff --git a/PWG/mftmuondep/CMakeLists.txt b/PWG/mftmuondep/CMakeLists.txt index ceda7e23e75..ac5ba0b857b 100644 --- a/PWG/mftmuondep/CMakeLists.txt +++ b/PWG/mftmuondep/CMakeLists.txt @@ -15,6 +15,6 @@ # Include the libraries add_subdirectory(charmonium) -add_subdirectory(openHF) +#add_subdirectory(openHF) message(STATUS "PWG mftmuondep enabled") diff --git a/PWG/mftmuondep/charmonium/AliAnalysisTaskDimuonBackground.cxx b/PWG/mftmuondep/charmonium/AliAnalysisTaskDimuonBackground.cxx new file mode 100644 index 00000000000..03a59597792 --- /dev/null +++ b/PWG/mftmuondep/charmonium/AliAnalysisTaskDimuonBackground.cxx @@ -0,0 +1,525 @@ +#include "TH1D.h" +#include "AliAODEvent.h" +#include "AliAODMCParticle.h" +#include "AliAODMCHeader.h" +#include "AliAnalysisTaskDimuonBackground.h" +#include "TDatabasePDG.h" +#include "AliMUONTrackParam.h" +#include "AliMUONTrackExtrap.h" +#include "AliAODDimuon.h" +#include "AliAODTrack.h" +#include "AliAODHeader.h" +#include "TClonesArray.h" +#include "AliMFTConstants.h" +#include "AliMFTAnalysisTools.h" +#include "TRandom.h" +#include "TList.h" + +ClassImp(AliAnalysisTaskDimuonBackground) + +//==================================================================================================================================================== + +AliAnalysisTaskDimuonBackground::AliAnalysisTaskDimuonBackground() : + AliAnalysisTaskSE(), + fVertexMode(0), + fMinTriggerMatch(0), + fSingleMuonMinEta(-9999), + fSingleMuonMaxEta(9999), + fSingleMuonMinPt(0), + fSingleMuonMaxChi2(9999), + fHistogramList(0), + fMainInputHandler(0), + fMixingInputHandler(0) +{ + + // Default constructor + + fHistSingleMuonsChi2VsPtVsRapidityBeforeCut = 0; + fHistSingleMuonsChi2VsPtVsRapidityAfterCut = 0; + + for (Int_t i=0; i<2; i++) fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity[i] = 0; + fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity_Resonances = 0; + + for (Int_t i=0; i<3; i++) { fPrimaryVertex[i]=0; fPrimaryVertexMixEv[i]=0; fPrimaryVertexTrue[i]=0; fPrimaryVertexMixEvTrue[i]=0; } + fVtxResolutionITS[0] = 5.e-4; + fVtxResolutionITS[1] = 5.e-4; + fVtxResolutionITS[2] = 4.e-4; + + fMassJpsi = TDatabasePDG::Instance()->GetParticle(443)->Mass(); + +} + +//==================================================================================================================================================== + +AliAnalysisTaskDimuonBackground::AliAnalysisTaskDimuonBackground(const Char_t *name) : + AliAnalysisTaskSE(name), + fVertexMode(0), + fMinTriggerMatch(0), + fSingleMuonMinEta(-9999), + fSingleMuonMaxEta(9999), + fSingleMuonMinPt(0), + fSingleMuonMaxChi2(9999), + fHistogramList(0), + fMainInputHandler(0), + fMixingInputHandler(0) +{ + + // Constructor + + fHistSingleMuonsChi2VsPtVsRapidityBeforeCut = 0; + fHistSingleMuonsChi2VsPtVsRapidityAfterCut = 0; + + for (Int_t i=0; i<2; i++) fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity[i] = 0; + fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity_Resonances = 0; + + for (Int_t i=0; i<3; i++) { fPrimaryVertex[i]=0; fPrimaryVertexMixEv[i]=0; fPrimaryVertexTrue[i]=0; fPrimaryVertexMixEvTrue[i]=0; } + fVtxResolutionITS[0] = 5.e-4; + fVtxResolutionITS[1] = 5.e-4; + fVtxResolutionITS[2] = 4.e-4; + + fMassJpsi = TDatabasePDG::Instance()->GetParticle(443)->Mass(); + + // Define input and output slots here + DefineOutput(1, TList::Class()); + +} + +//==================================================================================================================================================== + +void AliAnalysisTaskDimuonBackground::UserCreateOutputObjects() { + + // Called once + + fHistogramList = new TList(); + fHistogramList->SetOwner(kTRUE); + + // Max 100 dimensions + Int_t nBins[100] = {0}; + Double_t min[100] = {0}; + Double_t max[100] = {0}; + + TString tag[2] = {"SingleEvents", "MixedEvents"}; + + //------------------------------------------------------------------------------- + + nBins[0] = 100; + nBins[1] = 1000; + nBins[2] = 1000; + nBins[3] = 100; + nBins[4] = 100; + + min[0] = 0.; + min[1] = -10.; + min[2] = 0.; + min[3] = 0.; + min[4] = -4.5; + + max[0] = 1.; + max[1] = 10.; + max[2] = 10.; + max[3] = 10.; + max[4] = -2.0; + + for (Int_t i=0; i<2; i++) { + fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity[i] = new THnSparseD(Form("fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity_%s",tag[i].Data()), "", 5, nBins, min, max); + fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity[i] -> Sumw2(); + fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity[i] -> GetAxis(0)->SetTitle("PCA Quality"); + fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity[i] -> GetAxis(1)->SetTitle("t_{z} [ps]"); + fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity[i] -> GetAxis(2)->SetTitle("Mass [GeV/c^{2}]"); + fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity[i] -> GetAxis(3)->SetTitle("p_{T} [GeV/c]"); + fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity[i] -> GetAxis(4)->SetTitle("Rapidity"); + fHistogramList->Add(fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity[i]); + } + + fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity_Resonances = new THnSparseD("fHistPCAQualityVsPPDecayTimeVsMassVsPt_Resonances", "", 5, nBins, min, max); + fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity_Resonances -> Sumw2(); + fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity_Resonances -> GetAxis(0)->SetTitle("PCA Quality"); + fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity_Resonances -> GetAxis(1)->SetTitle("t_{z} [ps]"); + fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity_Resonances -> GetAxis(2)->SetTitle("Mass [GeV/c^{2}]"); + fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity_Resonances -> GetAxis(3)->SetTitle("p_{T} [GeV/c]"); + fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity_Resonances -> GetAxis(4)->SetTitle("Rapidity"); + fHistogramList->Add(fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity_Resonances); + + nBins[0] = 200; + nBins[1] = 100; + nBins[2] = 100; + + min[0] = 0.; + min[1] = 0.; + min[2] = -4.5; + + max[0] = 10.; + max[1] = 10.; + max[2] = -2.0; + + fHistSingleMuonsChi2VsPtVsRapidityBeforeCut = new THnSparseD("fHistSingleMuonsChi2VsPtVsRapidityBeforeCut", "", 3, nBins, min, max); + fHistSingleMuonsChi2VsPtVsRapidityBeforeCut -> Sumw2(); + fHistSingleMuonsChi2VsPtVsRapidityBeforeCut -> GetAxis(0)->SetTitle("p_{T} [GeV/c]"); + fHistSingleMuonsChi2VsPtVsRapidityBeforeCut -> GetAxis(1)->SetTitle("#chi^{2}"); + fHistSingleMuonsChi2VsPtVsRapidityBeforeCut -> GetAxis(2)->SetTitle("Rapidity"); + fHistogramList->Add(fHistSingleMuonsChi2VsPtVsRapidityBeforeCut); + + fHistSingleMuonsChi2VsPtVsRapidityAfterCut = new THnSparseD("fHistSingleMuonsChi2VsPtVsRapidityAfterCut", "", 3, nBins, min, max); + fHistSingleMuonsChi2VsPtVsRapidityAfterCut -> Sumw2(); + fHistSingleMuonsChi2VsPtVsRapidityAfterCut -> GetAxis(0)->SetTitle("p_{T} [GeV/c]"); + fHistSingleMuonsChi2VsPtVsRapidityAfterCut -> GetAxis(1)->SetTitle("#chi^{2}"); + fHistSingleMuonsChi2VsPtVsRapidityAfterCut -> GetAxis(2)->SetTitle("Rapidity"); + fHistogramList->Add(fHistSingleMuonsChi2VsPtVsRapidityAfterCut); + + //------------------------------------------------------------------------------- + + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + SetMainInputHandler(mgr); + if (fMainInputHandler) SetMixingInputHandler(fMainInputHandler); + + PostData(1, fHistogramList); + +} + +//==================================================================================================================================================== + +void AliAnalysisTaskDimuonBackground::UserExec(Option_t *) { + + // Main loop + // Called for each event + + AliAODEvent *aodEv = dynamic_cast(GetMainEvent()); + + if (!aodEv) return; + + if (!(AliMUONTrackExtrap::IsFieldON())) { + ((AliAODHeader*) aodEv->GetHeader())->InitMagneticField(); + AliMUONTrackExtrap::SetField(); + } + + // Getting primary vertex, either from the generation or from the reconstruction ------------------- + + AliAODMCHeader *mcHeader = (AliAODMCHeader*) (aodEv->GetList()->FindObject(AliAODMCHeader::StdBranchName())); + + if (fVertexMode == kGenerated) { + mcHeader->GetVertex(fPrimaryVertex); + for (Int_t i=0; i<3; i++) fPrimaryVertex[i] = gRandom->Gaus(fPrimaryVertex[i], fVtxResolutionITS[i]); + } + else if (fVertexMode == kReconstructed) { + aodEv->GetPrimaryVertex()->GetXYZ(fPrimaryVertex); + } + + TClonesArray *stackMC = (TClonesArray*) (aodEv->GetList()->FindObject(AliAODMCParticle::StdBranchName())); + + AliAODTrack *recMuon[2] = {0}; + AliAODMCParticle *mcMuon[2]={0}, *mcMother[2]={0}; + Double_t var[100] = {0}; + + //-------------------------------------------------------------------------------- + + for (Int_t iTrack=0; iTrackGetNumberOfTracks(); iTrack++) { + + recMuon[0] = (AliAODTrack*) aodEv->GetTrack(iTrack); + if (!(recMuon[0]->IsMuonGlobalTrack())) continue; + if (AliMFTAnalysisTools::IsTrackInjected(recMuon[0], mcHeader, stackMC)) continue; // Only HIJING tracks are considered to build the background + + var[0] = recMuon[0]->Pt(); + var[1] = recMuon[0]->Chi2perNDF(); + var[2] = recMuon[0]->Y(); + + fHistSingleMuonsChi2VsPtVsRapidityBeforeCut -> Fill(var); + + if (!IsSingleMuonCutPassed(recMuon[0])) continue; + + fHistSingleMuonsChi2VsPtVsRapidityAfterCut -> Fill(var); + + mcMuon[0] = NULL; + mcMother[0] = NULL; + + if (recMuon[0]->GetLabel()>=0) mcMuon[0] = (AliAODMCParticle*) stackMC->At(recMuon[0]->GetLabel()); + if (mcMuon[0] && mcMuon[0]->GetMother()>=0) mcMother[0] = (AliAODMCParticle*) stackMC->At(mcMuon[0]->GetMother()); + + for (Int_t jTrack=0; jTrackGetTrack(jTrack); + if (!(recMuon[1]->IsMuonGlobalTrack())) continue; + if (AliMFTAnalysisTools::IsTrackInjected(recMuon[1], mcHeader, stackMC)) continue; // Only HIJING tracks are considered to build the background + if (!IsSingleMuonCutPassed(recMuon[1])) continue; + + mcMuon[1] = NULL; + mcMother[1] = NULL; + + if (recMuon[1]->GetLabel()>=0) mcMuon[1] = (AliAODMCParticle*) stackMC->At(recMuon[1]->GetLabel()); + if (mcMuon[1] && mcMuon[1]->GetMother()>=0) mcMother[1] = (AliAODMCParticle*) stackMC->At(mcMuon[1]->GetMother()); + + TObjArray *dimuon = new TObjArray(); + dimuon -> Add(recMuon[0]); + dimuon -> Add(recMuon[1]); + + if (recMuon[0]->Charge() + recMuon[1]->Charge()) { + delete dimuon; + continue; + } + + Bool_t isDimuonResonance = kFALSE; + if (mcMother[0] && mcMother[1] && mcMuon[0]->GetMother()==mcMuon[1]->GetMother()) { + if (AliMFTAnalysisTools::IsPDGResonance(mcMother[0]->GetPdgCode())) isDimuonResonance = kTRUE; + } + + Double_t pca[3]={0}; + Double_t pcaQuality=0; + TLorentzVector kinem(0,0,0,0); + if (!AliMFTAnalysisTools::CalculatePCA(dimuon, pca, pcaQuality, kinem)) { + delete dimuon; + continue; + } + + Double_t ppDecayTime = AliMFTAnalysisTools::GetPseudoProperDecayTimeZ(fPrimaryVertex[2], pca[2], fMassJpsi, kinem.Pz()); + + var[0] = pcaQuality; + var[1] = ppDecayTime; + var[2] = kinem.M(); + var[3] = kinem.Pt(); + var[4] = kinem.Rapidity(); + + if (!isDimuonResonance) fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity[kSingleEvents] -> Fill(var); + else fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity_Resonances -> Fill(var); + + delete dimuon; + + } // end of loop on 2nd muon + + } // end of loop on 1st muon + + //-------------------------------------------------------------------------------- + + PostData(1, fHistogramList); + +} + +//==================================================================================================================================================== + +void AliAnalysisTaskDimuonBackground::UserExecMix(Option_t *) { + + // Performing Mixing + + if (!fMixingInputHandler) return; + + Int_t bufferSize = fMixingInputHandler->BufferSize(); + + AliAODEvent *aodEv = dynamic_cast(GetMainEvent()); + + if (!aodEv) return; + + if (!(AliMUONTrackExtrap::IsFieldON())) { + ((AliAODHeader*) aodEv->GetHeader())->InitMagneticField(); + AliMUONTrackExtrap::SetField(); + } + + // Getting primary vertex, either from the generation or from the reconstruction ------------------- + + AliAODMCHeader *mcHeader = (AliAODMCHeader*) (aodEv->GetList()->FindObject(AliAODMCHeader::StdBranchName())); + mcHeader->GetVertex(fPrimaryVertexTrue); + + if (fVertexMode == kGenerated) { + mcHeader->GetVertex(fPrimaryVertex); + for (Int_t i=0; i<3; i++) fPrimaryVertex[i] = gRandom->Gaus(fPrimaryVertex[i], fVtxResolutionITS[i]); + } + else if (fVertexMode == kReconstructed) { + aodEv->GetPrimaryVertex()->GetXYZ(fPrimaryVertex); + } + + TClonesArray *stackMC = (TClonesArray*) (aodEv->GetList()->FindObject(AliAODMCParticle::StdBranchName())); + + // ------------------------------------------------------------------------------------------------- + + AliAODTrack *recMuon[2] = {0}; + Double_t var[100] = {0}; + + for (Int_t iBuffer=0; iBuffer(GetMixedEvent(iBuffer)); + if (!aodEvMix) continue; + + // Getting primary vertex, either from the generation or from the reconstruction ------------------- + + AliAODMCHeader *mcHeaderMixEv = (AliAODMCHeader*) (aodEvMix->GetList()->FindObject(AliAODMCHeader::StdBranchName())); + mcHeaderMixEv->GetVertex(fPrimaryVertexMixEvTrue); + + if (fVertexMode == kGenerated) { + mcHeaderMixEv->GetVertex(fPrimaryVertexMixEv); + for (Int_t i=0; i<3; i++) fPrimaryVertexMixEv[i] = gRandom->Gaus(fPrimaryVertexMixEv[i], fVtxResolutionITS[i]); + } + else if (fVertexMode == kReconstructed) { + aodEvMix->GetPrimaryVertex()->GetXYZ(fPrimaryVertexMixEv); + } + + TClonesArray *stackMCMixEv = (TClonesArray*) (aodEvMix->GetList()->FindObject(AliAODMCParticle::StdBranchName())); + + // Loop over MUON+MFT muons of 1st and 2nd events + + //-------------------------------------------------------------------------------------------------- + + for (Int_t iTrack=0; iTrackGetNumberOfTracks(); iTrack++) { + + recMuon[0] = (AliAODTrack*) aodEv->GetTrack(iTrack); + if (!(recMuon[0]->IsMuonGlobalTrack())) continue; + if (AliMFTAnalysisTools::IsTrackInjected(recMuon[0], mcHeader, stackMC)) continue; // Only HIJING tracks are considered to build the background + if (!IsSingleMuonCutPassed(recMuon[0])) continue; + + for (Int_t jTrack=0; jTrackGetNumberOfTracks(); jTrack++) { + + recMuon[1] = (AliAODTrack*) aodEvMix->GetTrack(jTrack); + if (!(recMuon[1]->IsMuonGlobalTrack())) continue; + if (AliMFTAnalysisTools::IsTrackInjected(recMuon[1], mcHeaderMixEv, stackMCMixEv)) continue; // Only HIJING tracks are considered to build the background + if (!IsSingleMuonCutPassed(recMuon[1])) continue; + + // ----------- Translating muons to a common primary vertex (the origin). Preserve original tracks in case one wants to use them + + AliAODTrack *recMuonTranslated[2] = {0}; + for (Int_t i=0; i<2; i++) recMuonTranslated[i] = new AliAODTrack(*(recMuon[i])); + + if (!(AliMFTAnalysisTools::TranslateMuonToOrigin(recMuonTranslated[0], fPrimaryVertexTrue)) || + !(AliMFTAnalysisTools::TranslateMuonToOrigin(recMuonTranslated[1], fPrimaryVertexMixEvTrue))) { + printf("Error: tracks cannot be translated!!!\n"); + delete recMuonTranslated[0]; + delete recMuonTranslated[1]; + continue; + } + + TObjArray *dimuon = new TObjArray(); + dimuon -> Add(recMuonTranslated[0]); + dimuon -> Add(recMuonTranslated[1]); + + if (recMuonTranslated[0]->Charge() + recMuonTranslated[1]->Charge()) { + delete dimuon; + delete recMuonTranslated[0]; + delete recMuonTranslated[1]; + continue; + } + + Double_t pca[3]={0}; + Double_t pcaQuality=0; + TLorentzVector kinem(0,0,0,0); + if (!AliMFTAnalysisTools::CalculatePCA(dimuon, pca, pcaQuality, kinem)) { + delete dimuon; + delete recMuonTranslated[0]; + delete recMuonTranslated[1]; + continue; + } + + Double_t ppDecayTime = AliMFTAnalysisTools::GetPseudoProperDecayTimeZ(fPrimaryVertex[2]-fPrimaryVertexTrue[2], pca[2], fMassJpsi, kinem.Pz()); + + var[0] = pcaQuality; + var[1] = ppDecayTime; + var[2] = kinem.M(); + var[3] = kinem.Pt(); + var[4] = kinem.Rapidity(); + + fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity[kMixedEvents] -> Fill(var); + + delete dimuon; + delete recMuonTranslated[0]; + delete recMuonTranslated[1]; + + } // end of loop on 2nd muon + + } // end of loop on 1st muon + + //-------------------------------------------------------------------------------------------------- + + } // end of loop on events in buffer + + PostData(1, fHistogramList); + +} + +//==================================================================================================================================================== + +void AliAnalysisTaskDimuonBackground::Terminate(Option_t *) { + + // Draw result to the screen + // Called once at the end of the query + +} + +//==================================================================================================================================================== + +AliVEvent *AliAnalysisTaskDimuonBackground::GetMainEvent() { + + // Access to MainEvent + + AliMultiInputEventHandler *inEvHMainMulti = fMainInputHandler; + if (inEvHMainMulti) { + AliInputEventHandler *inEvMain = dynamic_cast(inEvHMainMulti->GetFirstInputEventHandler()); + if (inEvMain) return inEvMain->GetEvent(); + } + + return 0; + +} + +//==================================================================================================================================================== + +AliVEvent *AliAnalysisTaskDimuonBackground::GetMixedEvent(Int_t buffId) { + + // Access to Mixed event with buffer id + + AliMultiInputEventHandler *inEvHMain = fMainInputHandler; + + if (inEvHMain) { + + AliMixInputEventHandler *mixIH = fMixingInputHandler; + if (!mixIH) return 0; + + if (mixIH->CurrentBinIndex() < 0) { + AliDebug(AliLog::kDebug + 1, "Current event mixEH->CurrentEntry() == -1"); + return 0; + } + + AliMultiInputEventHandler *inEvHMixedCurrent = dynamic_cast(mixIH->InputEventHandler(buffId)); + if (!inEvHMixedCurrent) return 0; + + AliInputEventHandler *ihMixedCurrent = inEvHMixedCurrent->GetFirstInputEventHandler(); + if (ihMixedCurrent) return ihMixedCurrent->GetEvent(); + + } + + return 0; + +} + +//==================================================================================================================================================== + +AliMultiInputEventHandler *AliAnalysisTaskDimuonBackground::SetMainInputHandler(AliAnalysisManager *mgr) { + + // Sets main input handler + + if (!fMainInputHandler) fMainInputHandler = dynamic_cast(mgr->GetInputEventHandler()); + + return fMainInputHandler; + +} + +//==================================================================================================================================================== + +AliMixInputEventHandler *AliAnalysisTaskDimuonBackground::SetMixingInputHandler(AliMultiInputEventHandler *mainIH) { + + // Sets mixing input handler + + if (!fMixingInputHandler) fMixingInputHandler = dynamic_cast(mainIH->GetFirstMultiInputHandler()); + + return fMixingInputHandler; + +} + +//==================================================================================================================================================== + +Bool_t AliAnalysisTaskDimuonBackground::IsSingleMuonCutPassed(AliAODTrack *mu) { + + if (mu->GetMatchTrigger() < fMinTriggerMatch) return kFALSE; + if (mu->Eta()Eta()>fSingleMuonMaxEta) return kFALSE; + if (mu->Pt() < fSingleMuonMinPt) return kFALSE; + if (mu->Chi2perNDF() > fSingleMuonMaxChi2) return kFALSE; + + return kTRUE; + +} + +//==================================================================================================================================================== + diff --git a/PWG/mftmuondep/charmonium/AliAnalysisTaskDimuonBackground.h b/PWG/mftmuondep/charmonium/AliAnalysisTaskDimuonBackground.h new file mode 100755 index 00000000000..d4293bfb81d --- /dev/null +++ b/PWG/mftmuondep/charmonium/AliAnalysisTaskDimuonBackground.h @@ -0,0 +1,80 @@ +#ifndef AliAnalysisTaskDimuonBackground_H +#define AliAnalysisTaskDimuonBackground_H + +#include "AliAnalysisTaskSE.h" +#include "AliAnalysisManager.h" +#include "AliMultiInputEventHandler.h" +#include "AliMixInputEventHandler.h" +#include "TH1D.h" +#include "TList.h" +#include "AliAODTrack.h" +#include "THnSparse.h" + +//==================================================================================================================================================== + +class AliAnalysisTaskDimuonBackground : public AliAnalysisTaskSE { + +public: + + enum {kGenerated, kReconstructed}; + enum {kSingleEvents, kMixedEvents}; + + AliAnalysisTaskDimuonBackground(); + AliAnalysisTaskDimuonBackground(const char *name); + + virtual ~AliAnalysisTaskDimuonBackground() { if (fHistogramList) delete fHistogramList; } + + void SetVertexMode(Int_t vertexMode) { fVertexMode = vertexMode; } + void SetVtxResolutionITS(Double_t sigmaX, Double_t sigmaY, Double_t sigmaZ) { + fVtxResolutionITS[0] = sigmaX; + fVtxResolutionITS[1] = sigmaY; + fVtxResolutionITS[2] = sigmaZ; + } + + void SetMinTriggerMatch(Int_t minTriggerMatch) { fMinTriggerMatch = minTriggerMatch; } + + void SetSingleMuonMinPt(Double_t minPt) { fSingleMuonMinPt = minPt; } + void SetSingleMuonMinEta(Double_t minEta) { fSingleMuonMinEta = minEta; } + void SetSingleMuonMaxEta(Double_t maxEta) { fSingleMuonMaxEta = maxEta; } + void SetSingleMuonMaxChi2(Double_t maxChi2) { fSingleMuonMaxChi2 = maxChi2; } + + virtual void UserCreateOutputObjects(); + virtual void UserExec(Option_t *option); + virtual void UserExecMix(Option_t *); + virtual void Terminate(Option_t *); + + AliVEvent* GetMainEvent(); + AliVEvent* GetMixedEvent(Int_t buffId=0); + + AliMultiInputEventHandler* SetMainInputHandler(AliAnalysisManager *mgr); + AliMixInputEventHandler* SetMixingInputHandler(AliMultiInputEventHandler *mainIH); + + Bool_t IsSingleMuonCutPassed(AliAODTrack *mu); + +private: + + Double_t fMassJpsi; + + Double_t fPrimaryVertexTrue[3], fPrimaryVertexMixEvTrue[3], fPrimaryVertex[3], fPrimaryVertexMixEv[3], fVtxResolutionITS[3]; + Int_t fVertexMode; + + Int_t fMinTriggerMatch; + Double_t fSingleMuonMinEta, fSingleMuonMaxEta, fSingleMuonMinPt, fSingleMuonMaxChi2; + + TList *fHistogramList; + + THnSparse *fHistSingleMuonsChi2VsPtVsRapidityBeforeCut; //! + THnSparse *fHistSingleMuonsChi2VsPtVsRapidityAfterCut; //! + THnSparse *fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity[2]; //! + THnSparse *fHistPCAQualityVsPPDecayTimeVsMassVsPtVsRapidity_Resonances; //! + + AliMultiInputEventHandler *fMainInputHandler; //! tmp pointer to main input handler + AliMixInputEventHandler *fMixingInputHandler; //! tmp pointer to mixing input handler + + ClassDef(AliAnalysisTaskDimuonBackground, 1) // example of analysis + +}; + +//==================================================================================================================================================== + +#endif diff --git a/PWG/mftmuondep/charmonium/CMakeLists.txt b/PWG/mftmuondep/charmonium/CMakeLists.txt index ca37268ea2d..c6526d52c07 100644 --- a/PWG/mftmuondep/charmonium/CMakeLists.txt +++ b/PWG/mftmuondep/charmonium/CMakeLists.txt @@ -11,51 +11,20 @@ # * appear in the supporting documentation. The authors make no claims * # * about the suitability of this software for any purpose. It is * # * provided "as is" without express or implied warranty. * -# **************************************************************************/ +# ************************************************************************** -#Module -set(MODULE PWGSomething) +# Module +set (MODULE PWGmftmuondepCharmonium) # Module include folder include_directories(${AliRoot_SOURCE_DIR}/PWG/mftmuondep/charmonium) -# Additional includes - alphabetical order except ROOT -include_directories(${ROOT_INCLUDE_DIRS} - ) - # Sources - alphabetical order set(SRCS + AliAnalysisTaskDimuonBackground.cxx ) # Headers from sources string(REPLACE ".cxx" ".h" HDRS "${SRCS}") -# Generate the dictionary -# It will create G_ARG1.cxx and G_ARG1.h / ARG1 = function first argument -get_directory_property(incdirs INCLUDE_DIRECTORIES) -generate_dictionary("${MODULE}" "${MODULE}LinkDef.h" "${HDRS}" "${incdirs}") - -# Add a shared library -add_library(${MODULE} SHARED ${SRCS} G__${MODULE}.cxx) - -# Generate the ROOT map -# Dependecies -set(LIBDEPS) -generate_rootmap("${MODULE}" "${LIBDEPS}" "${CMAKE_CURRENT_SOURCE_DIR}/${MODULE}LinkDef.h") - -# Linking the library -target_link_libraries(${MODULE} ${LIBDEPS}) - -# Public include folders that will be propagated to the dependecies -target_include_directories(${MODULE} PUBLIC ${incdirs}) - -# System dependent: Modify the way the library is build -if(${CMAKE_SYSTEM} MATCHES Darwin) - set_target_properties(${MODULE} PROPERTIES LINK_FLAGS "-undefined dynamic_lookup") -endif(${CMAKE_SYSTEM} MATCHES Darwin) - -# Installation -install(TARGETS ${MODULE} - ARCHIVE DESTINATION lib - LIBRARY DESTINATION lib) -install(FILES ${HDRS} DESTINATION include) \ No newline at end of file +message(STATUS "PWG mftmuondep charmonium enabled") diff --git a/PWG/mftmuondep/charmonium/PWGmftmuondepCharmoniumLinkDef.h b/PWG/mftmuondep/charmonium/PWGmftmuondepCharmoniumLinkDef.h new file mode 100644 index 00000000000..ba5509f9501 --- /dev/null +++ b/PWG/mftmuondep/charmonium/PWGmftmuondepCharmoniumLinkDef.h @@ -0,0 +1,11 @@ +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class AliAnalysisTaskDimuonBackground.cxx + +#endif + + diff --git a/PWG/mftmuondep/charmonium/RunAnalysisTaskDimuonBackground.C b/PWG/mftmuondep/charmonium/RunAnalysisTaskDimuonBackground.C new file mode 100755 index 00000000000..dab98750cf6 --- /dev/null +++ b/PWG/mftmuondep/charmonium/RunAnalysisTaskDimuonBackground.C @@ -0,0 +1,297 @@ +//==================================================================================================================================================== + +Bool_t RunAnalysisTaskDimuonBackground(const Char_t *runType="local", // "grid" and "local" types have been tested + const Char_t *runMode="full") { + + // enum {kGenerated, kReconstructed}; + + LoadLibs(); + + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if (mgr) delete mgr; + mgr = new AliAnalysisManager("AM","Analysis Manager"); + + AliMultiInputEventHandler *handler = new AliMultiInputEventHandler(); + handler->AddInputEventHandler(new AliAODInputHandler()); + + //-------------------- Setup of the mixed event handler --------------- + + const Int_t mixNum = 20; // Number of times the main event is used for a mixing procedure + const Int_t bufferSize = 1; // Number of events involved in each mixing procedure with the main one (typically one) + + AliMixInputEventHandler *mixHandler = new AliMixInputEventHandler(bufferSize, mixNum); + mixHandler->SetInputHandlerForMixing(dynamic_cast handler); + AliMixEventPool *evPool = new AliMixEventPool(); + + AliMixEventCutObj *zVertex = new AliMixEventCutObj(AliMixEventCutObj::kZVertex, -10, 10, 1); + + evPool->AddCut(zVertex); + mixHandler->SetEventPool(evPool); + + mixHandler->SelectCollisionCandidates(AliVEvent::kAny); + handler->AddInputEventHandler(mixHandler); + + //--------------------------------------------------------------------- + + mgr->SetInputEventHandler(handler); + + if (!strcmp(runType,"grid")) mgr->SetGridHandler(CreateAlienHandler(runMode)); + + gROOT->LoadMacro("AliAnalysisTaskDimuonBackground.cxx++g"); + AliAnalysisTaskDimuonBackground *task = new AliAnalysisTaskDimuonBackground("AliAnalysisTaskDimuonBackground"); + + // in cm, taken from Fig. 7.4 of the ITS Upgrade TDR, in the hypothesis of ~80 tracks participating to the vtx + // task -> SetVtxResolutionITS(5.e-4, 5.e-4, 4.e-4); + + task -> SetVertexMode(AliAnalysisTaskDimuonBackground::kReconstructed); + + task -> SetMinTriggerMatch(1); + task -> SetSingleMuonMinPt(1.0); + task -> SetSingleMuonMinEta(-3.6); + task -> SetSingleMuonMaxEta(-2.5); + task -> SetSingleMuonMaxChi2(9999.); + + // create output container(s) + AliAnalysisDataContainer *histogramList = mgr->CreateContainer("list", TList::Class(), AliAnalysisManager::kOutputContainer, "AnalysisMFT_Output.root"); + + // connect input and output + mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer()); + mgr->ConnectOutput(task, 1, histogramList); + + // RUN ANALYSIS + + TStopwatch timer; + timer.Start(); + mgr->InitAnalysis(); + mgr->PrintStatus(); + + if (!strcmp(runType,"grid")) { + printf("Starting MFT analysis on the grid"); + mgr->StartAnalysis("grid"); + } + + else if (!strcmp(runType,"local")) { + printf("Starting MFT analysis locally"); + mgr->StartAnalysis("local", GetInputLocalData()); + } + + else AliError(Form("Specified run type %s is not supported", runType)); + + timer.Stop(); + timer.Print(); + return kTRUE; + +} + +//==================================================================================================================================================== + +AliAnalysisGrid* CreateAlienHandler(const Char_t *runMode) { + + // Set up the analysis plugin in case of a grid analysis + + TString rootVersion = "v5-34-08-6"; + TString alirootVersion = "vAN-20141128"; + + AliAnalysisAlien *analysisPlugin = new AliAnalysisAlien(); + if (!analysisPlugin) { Printf("Error : analysisPlugin is null !!!"); return kFALSE; } + analysisPlugin->SetAPIVersion("V1.1x"); + analysisPlugin->SetROOTVersion(rootVersion.Data()); + analysisPlugin->SetAliROOTVersion(alirootVersion.Data()); + analysisPlugin->SetExecutableCommand("aliroot -b -q"); + + // Overwrite all generated files, datasets and output results from a previous session + analysisPlugin->SetOverwriteMode(); + // Set the run mode (can be "full", "test", "offline", "submit" or "terminate") + analysisPlugin->SetRunMode(runMode); // VERY IMPORTANT + + analysisPlugin->SetAdditionalRootLibs("CORRFW"); + analysisPlugin->SetAdditionalRootLibs("PWGmuon"); + + // Method1: Location of Data and Working dir + // analysisPlugin->SetGridDataDir("/alice/sim/2014/LHC14j5_new_plus/"); + // analysisPlugin->SetDataPattern("AliAOD.root"); + // analysisPlugin->AddRunNumber(137161); + + // Method2: Use existing xml collection + + analysisPlugin->AddDataFile("../../collections/137844_1.xml"); + analysisPlugin->AddDataFile("../../collections/137848_1.xml"); + analysisPlugin->AddDataFile("../../collections/138190_1.xml"); + analysisPlugin->AddDataFile("../../collections/138192_1.xml"); + analysisPlugin->AddDataFile("../../collections/138197_1.xml"); + analysisPlugin->AddDataFile("../../collections/138201_1.xml"); + analysisPlugin->AddDataFile("../../collections/138225_1.xml"); + analysisPlugin->AddDataFile("../../collections/138275_1.xml"); + analysisPlugin->AddDataFile("../../collections/138364_1.xml"); + analysisPlugin->AddDataFile("../../collections/138396_1.xml"); + + // analysisPlugin->AddDataFile("137230_2.xml"); + // analysisPlugin->AddDataFile("137162_2.xml"); + // analysisPlugin->AddDataFile("137161_2.xml"); + // analysisPlugin->AddDataFile("137231_2.xml"); + // analysisPlugin->AddDataFile("137430_2.xml"); + // analysisPlugin->AddDataFile("137366_2.xml"); + // analysisPlugin->AddDataFile("137243_2.xml"); + // analysisPlugin->AddDataFile("137236_2.xml"); + // analysisPlugin->AddDataFile("137232_2.xml"); + // analysisPlugin->AddDataFile("137440_2.xml"); + // analysisPlugin->AddDataFile("137434_2.xml"); + // analysisPlugin->AddDataFile("137432_2.xml"); + // analysisPlugin->AddDataFile("137431_2.xml"); + // analysisPlugin->AddDataFile("137443_2.xml"); + // analysisPlugin->AddDataFile("137441_2.xml"); + // analysisPlugin->AddDataFile("137544_2.xml"); + // analysisPlugin->AddDataFile("137541_2.xml"); + // analysisPlugin->AddDataFile("137539_2.xml"); + // analysisPlugin->AddDataFile("137549_2.xml"); + // analysisPlugin->AddDataFile("137608_2.xml"); + // analysisPlugin->AddDataFile("137595_2.xml"); + // analysisPlugin->AddDataFile("137686_2.xml"); + // analysisPlugin->AddDataFile("137639_2.xml"); + // analysisPlugin->AddDataFile("137638_2.xml"); + // analysisPlugin->AddDataFile("137692_2.xml"); + // analysisPlugin->AddDataFile("137691_2.xml"); + // analysisPlugin->AddDataFile("137718_2.xml"); + // analysisPlugin->AddDataFile("137704_2.xml"); + // analysisPlugin->AddDataFile("137751_2.xml"); + // analysisPlugin->AddDataFile("137724_2.xml"); + // analysisPlugin->AddDataFile("137722_2.xml"); + // analysisPlugin->AddDataFile("137752_2.xml"); + // analysisPlugin->AddDataFile("137844_2.xml"); + // analysisPlugin->AddDataFile("138190_2.xml"); + // analysisPlugin->AddDataFile("137848_2.xml"); + // analysisPlugin->AddDataFile("138197_2.xml"); + // analysisPlugin->AddDataFile("138192_2.xml"); + // analysisPlugin->AddDataFile("138225_2.xml"); + // analysisPlugin->AddDataFile("138201_2.xml"); + // analysisPlugin->AddDataFile("138275_2.xml"); + // analysisPlugin->AddDataFile("138364_2.xml"); + // analysisPlugin->AddDataFile("138396_2.xml"); + + // analysisPlugin->AddDataFile("137161_3.xml"); + // analysisPlugin->AddDataFile("137162_3.xml"); + // analysisPlugin->AddDataFile("137230_3.xml"); + // analysisPlugin->AddDataFile("137231_3.xml"); + // analysisPlugin->AddDataFile("137232_3.xml"); + // analysisPlugin->AddDataFile("137236_3.xml"); + // analysisPlugin->AddDataFile("137243_3.xml"); + // analysisPlugin->AddDataFile("137366_3.xml"); + // analysisPlugin->AddDataFile("137430_3.xml"); + // analysisPlugin->AddDataFile("137431_3.xml"); + // analysisPlugin->AddDataFile("137432_3.xml"); + // analysisPlugin->AddDataFile("137434_3.xml"); + // analysisPlugin->AddDataFile("137440_3.xml"); + // analysisPlugin->AddDataFile("137441_3.xml"); + // analysisPlugin->AddDataFile("137443_3.xml"); + // analysisPlugin->AddDataFile("137539_3.xml"); + // analysisPlugin->AddDataFile("137541_3.xml"); + // analysisPlugin->AddDataFile("137544_3.xml"); + // analysisPlugin->AddDataFile("137549_3.xml"); + // analysisPlugin->AddDataFile("137595_3.xml"); + // analysisPlugin->AddDataFile("137608_3.xml"); + // analysisPlugin->AddDataFile("137638_3.xml"); + // analysisPlugin->AddDataFile("137639_3.xml"); + // analysisPlugin->AddDataFile("137686_3.xml"); + // analysisPlugin->AddDataFile("137691_3.xml"); + // analysisPlugin->AddDataFile("137692_3.xml"); + // analysisPlugin->AddDataFile("137704_3.xml"); + // analysisPlugin->AddDataFile("137718_3.xml"); + // analysisPlugin->AddDataFile("137722_3.xml"); + // analysisPlugin->AddDataFile("137724_3.xml"); + // analysisPlugin->AddDataFile("137751_3.xml"); + // analysisPlugin->AddDataFile("137752_3.xml"); + // analysisPlugin->AddDataFile("137844_3.xml"); + // analysisPlugin->AddDataFile("137848_3.xml"); + // analysisPlugin->AddDataFile("138190_3.xml"); + // analysisPlugin->AddDataFile("138192_3.xml"); + // analysisPlugin->AddDataFile("138197_3.xml"); + // analysisPlugin->AddDataFile("138201_3.xml"); + // analysisPlugin->AddDataFile("138225_3.xml"); + // analysisPlugin->AddDataFile("138275_3.xml"); + // analysisPlugin->AddDataFile("138364_3.xml"); + // analysisPlugin->AddDataFile("138396_3.xml"); + + analysisPlugin->SetGridWorkingDir("MFT/analysis/LHC14cj5_PbPb_0_5/charmonia/background"); + + // Declare alien output directory. Relative to working directory. + analysisPlugin->SetGridOutputDir("output"); // In this case will be $HOME/work/output + // Declare the analysis source files names separated by blancs. To be compiled runtime using ACLiC on the worker nodes. + analysisPlugin->SetAnalysisSource("AliAnalysisTaskDimuonBackground.cxx"); + // Declare all libraries (other than the default ones for the framework. These will be + // loaded by the generated analysis macro. Add all extra files (task .cxx/.h) here. + analysisPlugin->SetAdditionalLibs("libCORRFW.so libPWGmuon.so libEventMixing.so AliAnalysisTaskDimuonBackground.h AliAnalysisTaskDimuonBackground.cxx"); + + analysisPlugin->AddIncludePath("-I."); + + analysisPlugin->SetOutputToRunNo(); + + // Optionally set a name for the generated analysis macro (default MyAnalysis.C) + analysisPlugin->SetAnalysisMacro("MFTAnalysis.C"); + + // Optionally set maximum number of input files/subjob (default 100, put 0 to ignore) + analysisPlugin->SetSplitMaxInputFileNumber(100); + // Number of runs per master job + analysisPlugin->SetNrunsPerMaster(1); + + // Optionally modify the executable name (default analysis.sh) + analysisPlugin->SetExecutable("MFTAnalysis.sh"); + + // Optionally set number of failed jobs that will trigger killing waiting sub-jobs. + analysisPlugin->SetMaxInitFailed(5); + // Optionally resubmit threshold. + analysisPlugin->SetMasterResubmitThreshold(90); + // Optionally set time to live (default 30000 sec) + analysisPlugin->SetTTL(30000); + // Optionally set input format (default xml-single) + analysisPlugin->SetInputFormat("xml-single"); + // Optionally modify the name of the generated JDL (default analysis.jdl) + analysisPlugin->SetJDLName("MFTAnalysis.jdl"); + // Optionally modify job price (default 1) + analysisPlugin->SetPrice(1); + // Optionally modify split mode (default 'se') + analysisPlugin->SetSplitMode("se"); + + analysisPlugin->SetNtestFiles(5); + // analysisPlugin->SetMergeViaJDL(1); + analysisPlugin->SetOverwriteMode(kTRUE); + + return analysisPlugin; + +} + +//==================================================================================================================================================== + +TChain* GetInputLocalData() { + + // Set up the chain of input events in case of a local analysis + + TChain *chain = new TChain("aodTree"); + chain->Add("./test_03Nov2014_zVtx_Rndm/AliAOD.Muons.root"); + + return chain; + +} + +//==================================================================================================================================================== + +void LoadLibs() { + + gSystem->AddIncludePath("-I$ALICE_ROOT/include "); + + gSystem->Load("libTree.so") ; + gSystem->Load("libGeom.so") ; + gSystem->Load("libVMC.so") ; + gSystem->Load("libMinuit.so") ; + gSystem->Load("libPhysics.so") ; + gSystem->Load("libSTEERBase.so") ; + gSystem->Load("libESD.so") ; + gSystem->Load("libAOD.so") ; + gSystem->Load("libANALYSIS.so") ; + gSystem->Load("libOADB.so") ; + gSystem->Load("libANALYSISalice.so") ; + gSystem->Load("libCORRFW.so") ; + gSystem->Load("libPWGmuon.so") ; + +} + +//==================================================================================================================================================== -- 2.39.3