From 6f719a5b2e598110e6c033b7d6769074e9de66b6 Mon Sep 17 00:00:00 2001 From: miweber Date: Wed, 7 Aug 2013 12:06:01 +0000 Subject: [PATCH] adding analysis option MCAODrec = AODs from MC with using MCtruth information for certain cuts (Alis Rodriguez Manso ) --- .../BalanceFunctions/AliAnalysisTaskBFPsi.cxx | 128 +++++++++++++++++- .../BalanceFunctions/AliAnalysisTaskBFPsi.h | 3 +- .../macros/AddTaskBalancePsiCentralityTrain.C | 16 +++ .../macros/configBalanceFunctionPsiAnalysis.C | 2 +- PWGCF/EBYE/macros/runBalanceFunctionPsi.C | 2 +- 5 files changed, 144 insertions(+), 7 deletions(-) diff --git a/PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.cxx b/PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.cxx index a3728433a02..72542863112 100755 --- a/PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.cxx +++ b/PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.cxx @@ -54,7 +54,8 @@ ClassImp(AliAnalysisTaskBFPsi) //________________________________________________________________________ AliAnalysisTaskBFPsi::AliAnalysisTaskBFPsi(const char *name) -: AliAnalysisTaskSE(name), +: AliAnalysisTaskSE(name), + fArrayMC(0), //+++++++++++++ fBalance(0), fRunShuffling(kFALSE), fShuffledBalance(0), @@ -70,7 +71,7 @@ AliAnalysisTaskBFPsi::AliAnalysisTaskBFPsi(const char *name) fHistListPIDQA(0), fHistEventStats(0), fHistCentStats(0), - fHistCentStatsUsed(0), //++++++++++++++++++++++++++ + fHistCentStatsUsed(0), fHistTriggerStats(0), fHistTrackStats(0), fHistVx(0), @@ -757,7 +758,7 @@ Double_t AliAnalysisTaskBFPsi::IsEventAccepted(AliVEvent *event){ //Centrality stuff if(fUseCentrality) { - if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD") { //centrality in AOD header + if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec") { //centrality in AOD header //+++++++++++ AliAODHeader *header = (AliAODHeader*) event->GetHeader(); if(header){ gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data()); @@ -961,7 +962,7 @@ Double_t AliAnalysisTaskBFPsi::GetRefMultiOrCentrality(AliVEvent *event){ TString gAnalysisLevel = fBalance->GetAnalysisLevel(); if(fEventClass == "Centrality"){ - if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD") { //centrality in AOD header + if(gAnalysisLevel == "AOD"|| gAnalysisLevel == "MCAOD" || gAnalysisLevel == "MCAODrec" ) { //centrality in AOD header //++++++++++++++ AliAODHeader *header = (AliAODHeader*) event->GetHeader(); if(header){ gCentrality = header->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data()); @@ -1332,6 +1333,125 @@ TObjArray* AliAnalysisTaskBFPsi::GetAcceptedTracks(AliVEvent *event, Double_t gC }//MCAOD //============================================================================================================== + //============================================================================================================== + else if(gAnalysisLevel == "MCAODrec") { + + /* fAOD = dynamic_cast(InputEvent()); + if (!fAOD) { + printf("ERROR: fAOD not available\n"); + return; + }*/ + + fArrayMC = dynamic_cast(event->FindListObject(AliAODMCParticle::StdBranchName())); + if (!fArrayMC) { + AliError("No array of MC particles found !!!"); + } + + AliMCEvent* mcEvent = MCEvent(); + if (!mcEvent) { + AliError("ERROR: Could not retrieve MC event"); + } + + for (Int_t iTracks = 0; iTracks < event->GetNumberOfTracks(); iTracks++) { + AliAODTrack* aodTrack = dynamic_cast(event->GetTrack(iTracks)); + if (!aodTrack) { + AliError(Form("Could not receive track %d", iTracks)); + continue; + } + + for(Int_t iTrackBit = 0; iTrackBit < 16; iTrackBit++){ + fHistTrackStats->Fill(iTrackBit,aodTrack->TestFilterBit(1<TestFilterBit(nAODtrackCutBit)) continue; + + vCharge = aodTrack->Charge(); + vEta = aodTrack->Eta(); + vY = aodTrack->Y(); + vPhi = aodTrack->Phi();// * TMath::RadToDeg(); + vPt = aodTrack->Pt(); + + Float_t dcaXY = aodTrack->DCA(); // this is the DCA from global track (not exactly what is cut on) + Float_t dcaZ = aodTrack->ZAtDCA(); // this is the DCA from global track (not exactly what is cut on) + + // Kinematics cuts from ESD track cuts + if( vPt < fPtMin || vPt > fPtMax) continue; + if( vEta < fEtaMin || vEta > fEtaMax) continue; + + // Extra DCA cuts (for systematic studies [!= -1]) + if( fDCAxyCut != -1 && fDCAzCut != -1){ + if(TMath::Sqrt((dcaXY*dcaXY)/(fDCAxyCut*fDCAxyCut)+(dcaZ*dcaZ)/(fDCAzCut*fDCAzCut)) > 1 ){ + continue; // 2D cut + } + } + + // Extra TPC cuts (for systematic studies [!= -1]) + if( fTPCchi2Cut != -1 && aodTrack->Chi2perNDF() > fTPCchi2Cut){ + continue; + } + if( fNClustersTPCCut != -1 && aodTrack->GetTPCNcls() < fNClustersTPCCut){ + continue; + } + + //Exclude resonances + if(fExcludeResonancesInMC) { + + Bool_t kExcludeParticle = kFALSE; + + Int_t label = TMath::Abs(aodTrack->GetLabel()); + AliAODMCParticle *AODmcTrack = (AliAODMCParticle*) fArrayMC->At(label); + + if (AODmcTrack){ + //if (AODmcTrack->IsPhysicalPrimary()){ + + Int_t gMotherIndex = AODmcTrack->GetMother(); + if(gMotherIndex != -1) { + AliAODMCParticle* motherTrack = dynamic_cast(mcEvent->GetTrack(gMotherIndex)); + if(motherTrack) { + Int_t pdgCodeOfMother = motherTrack->GetPdgCode(); + if(pdgCodeOfMother == 113 // rho0 + || pdgCodeOfMother == 213 || pdgCodeOfMother == -213 // rho+ + // || pdgCodeOfMother == 221 // eta + // || pdgCodeOfMother == 331 // eta' + // || pdgCodeOfMother == 223 // omega + // || pdgCodeOfMother == 333 // phi + || pdgCodeOfMother == 311 || pdgCodeOfMother == -311 // K0 + // || pdgCodeOfMother == 313 || pdgCodeOfMother == -313 // K0* + // || pdgCodeOfMother == 323 || pdgCodeOfMother == -323 // K+* + || pdgCodeOfMother == 3122 || pdgCodeOfMother == -3122 // Lambda + || pdgCodeOfMother == 111 // pi0 Dalitz + ) { + kExcludeParticle = kTRUE; + } + } + } + } + //Exclude from the analysis decay products of rho0, rho+, eta, eta' and phi + if(kExcludeParticle) continue; + } + + // fill QA histograms + fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls()); + fHistDCA->Fill(dcaZ,dcaXY); + fHistChi2->Fill(aodTrack->Chi2perNDF(),gCentrality); + fHistPt->Fill(vPt,gCentrality); + fHistEta->Fill(vEta,gCentrality); + fHistRapidity->Fill(vY,gCentrality); + if(vCharge > 0) fHistPhiPos->Fill(vPhi,gCentrality); + else if(vCharge < 0) fHistPhiNeg->Fill(vPhi,gCentrality); + fHistPhi->Fill(vPhi,gCentrality); + if(vCharge > 0) fHistEtaPhiPos->Fill(vEta,vPhi,gCentrality); + else if(vCharge < 0) fHistEtaPhiNeg->Fill(vEta,vPhi,gCentrality); + + //=======================================correction + Double_t correction = GetTrackbyTrackCorrectionMatrix(vEta, vPhi, vPt, vCharge, gCentrality); + //Printf("CORRECTIONminus: %.2f | Centrality %lf",correction,gCentrality); + + // add the track to the TObjArray + tracksAccepted->Add(new AliBFBasicParticle(vEta, vPhi, vPt, vCharge, correction)); + }//track loop + }//MCAODrec + //============================================================================================================== + else if(gAnalysisLevel == "ESD" || gAnalysisLevel == "MCESD") { // handling of TPC only tracks different in AOD and ESD AliESDtrack *trackTPC = NULL; diff --git a/PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.h b/PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.h index 716d3333951..6100b3d9d83 100755 --- a/PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.h +++ b/PWGCF/EBYE/BalanceFunctions/AliAnalysisTaskBFPsi.h @@ -172,7 +172,8 @@ class AliAnalysisTaskBFPsi : public AliAnalysisTaskSE { //===============================correction TObjArray* GetAcceptedTracks(AliVEvent* event, Double_t gCentrality, Double_t gReactionPlane); TObjArray* GetShuffledTracks(TObjArray* tracks, Double_t gCentrality); - + + TClonesArray* fArrayMC; //! AOD object //+++++++++++++++++++++ AliBalancePsi *fBalance; //BF object Bool_t fRunShuffling;//run shuffling or not AliBalancePsi *fShuffledBalance; //BF object (shuffled) diff --git a/PWGCF/EBYE/macros/AddTaskBalancePsiCentralityTrain.C b/PWGCF/EBYE/macros/AddTaskBalancePsiCentralityTrain.C index 44a66b86efb..8e813853766 100644 --- a/PWGCF/EBYE/macros/AddTaskBalancePsiCentralityTrain.C +++ b/PWGCF/EBYE/macros/AddTaskBalancePsiCentralityTrain.C @@ -104,6 +104,11 @@ AliAnalysisTaskBFPsi *AddTaskBalancePsiCentralityTrain(Double_t centrMin=0., if(gRunShuffling) bfs = GetBalanceFunctionObject("MCAOD",centralityEstimator,centrMin,centrMax,kTRUE,bResonancesCut,bHBTcut,bConversionCut,bMomentumDifferenceCut,fQCutMin,fArgEventClass,deltaEtaMax,bVertexBinning); if(gRunMixing) bfm = GetBalanceFunctionObject("MCAOD",centralityEstimator,centrMin,centrMax,kFALSE,bResonancesCut,bHBTcut,bConversionCut,bMomentumDifferenceCut,fQCutMin,fArgEventClass,deltaEtaMax,bVertexBinning); } + else if (analysisType=="MCAODrec"){ + bf = GetBalanceFunctionObject("MCAODrec",centralityEstimator,centrMin,centrMax,kFALSE,bResonancesCut,bHBTcut,bConversionCut,bMomentumDifferenceCut,fQCutMin,fArgEventClass,deltaEtaMax,bVertexBinning); + if(gRunShuffling) bfs = GetBalanceFunctionObject("MCAODrec",centralityEstimator,centrMin,centrMax,kTRUE,bResonancesCut,bHBTcut,bConversionCut,bMomentumDifferenceCut,fQCutMin,fArgEventClass,deltaEtaMax,bVertexBinning); + if(gRunMixing) bfm = GetBalanceFunctionObject("MCAODrec",centralityEstimator,centrMin,centrMax,kFALSE,bResonancesCut,bHBTcut,bConversionCut,bMomentumDifferenceCut,fQCutMin,fArgEventClass,deltaEtaMax,bVertexBinning); + } else{ ::Error("AddTaskBF", "analysis type NOT known."); return NULL; @@ -181,6 +186,17 @@ AliAnalysisTaskBFPsi *AddTaskBalancePsiCentralityTrain(Double_t centrMin=0., taskBF->SetAODtrackCutBit(AODfilterBit); taskBF->SetKinematicsCutsAOD(ptMin,ptMax,etaMin,etaMax); } + else if(analysisType == "MCAODrec") { //++++++++++++++++ + // pt and eta cut (pt_min, pt_max, eta_min, eta_max) + taskBF->SetAODtrackCutBit(AODfilterBit); + taskBF->SetKinematicsCutsAOD(ptMin,ptMax,etaMin,etaMax); + + // set extra DCA cuts (-1 no extra cut) + taskBF->SetExtraDCACutsAOD(DCAxy,DCAz); + + // set extra TPC chi2 / nr of clusters cut + taskBF->SetExtraTPCCutsAOD(maxTPCchi2, minNClustersTPC); + }//++++++++++++++++ // offline trigger selection (AliVEvent.h) // taskBF->UseOfflineTrigger(); // NOT used (selection is done with the AliAnalysisTaskSE::SelectCollisionCandidates()) diff --git a/PWGCF/EBYE/macros/configBalanceFunctionPsiAnalysis.C b/PWGCF/EBYE/macros/configBalanceFunctionPsiAnalysis.C index ea3343ff573..7c3d5a8f003 100644 --- a/PWGCF/EBYE/macros/configBalanceFunctionPsiAnalysis.C +++ b/PWGCF/EBYE/macros/configBalanceFunctionPsiAnalysis.C @@ -1,5 +1,5 @@ //__________________________________________________// -AliBalancePsi *GetBalanceFunctionObject(const char* analysisLevel = "AOD", +AliBalancePsi *GetBalanceFunctionObject(const char* analysisLevel = "MCAOD", // const char* centralityName = 0x0, Double_t centrMin = 0., Double_t centrMax = 100., diff --git a/PWGCF/EBYE/macros/runBalanceFunctionPsi.C b/PWGCF/EBYE/macros/runBalanceFunctionPsi.C index b07bf64be1f..d219e38c29a 100644 --- a/PWGCF/EBYE/macros/runBalanceFunctionPsi.C +++ b/PWGCF/EBYE/macros/runBalanceFunctionPsi.C @@ -133,7 +133,7 @@ void runBalanceFunctionPsi( }//local mode // input handler (ESD or AOD) - AliVEventHandler* inputH = NULL; + Aliveventhandler* inputH = NULL; if(!bAOD){ inputH = new AliESDInputHandler(); } -- 2.43.0