From 4aae2a936ea5b32c49a993154bbf0f04563da8e8 Mon Sep 17 00:00:00 2001 From: snelling Date: Fri, 22 Mar 2013 10:21:23 +0000 Subject: [PATCH] patching flow cuts (Panos) --- .../Base/AliFlowAnalysisWithQCumulants.cxx | 6 +- .../Base/AliFlowEventSimpleMakerOnTheFly.cxx | 82 ++++++++++++++++++- .../Base/AliFlowEventSimpleMakerOnTheFly.h | 25 ++++++ PWG/FLOW/Tasks/AliAnalysisTaskFlowEvent.cxx | 2 +- PWG/FLOW/Tasks/AliFlowEvent.cxx | 21 +++-- PWG/FLOW/Tasks/AliFlowEventCuts.cxx | 10 +-- PWG/FLOW/Tasks/AliFlowEventCuts.h | 5 +- PWG/FLOW/Tasks/AliFlowTrackCuts.cxx | 1 + 8 files changed, 134 insertions(+), 18 deletions(-) diff --git a/PWG/FLOW/Base/AliFlowAnalysisWithQCumulants.cxx b/PWG/FLOW/Base/AliFlowAnalysisWithQCumulants.cxx index a1ce18a735b..120b3df10af 100644 --- a/PWG/FLOW/Base/AliFlowAnalysisWithQCumulants.cxx +++ b/PWG/FLOW/Base/AliFlowAnalysisWithQCumulants.cxx @@ -3595,9 +3595,11 @@ void AliFlowAnalysisWithQCumulants::CalculateIntFlowCorrelations() Double_t dMultiplicityBin = 0.; if(fMultiplicityIs==AliFlowCommonConstants::kRP) { - dMultiplicityBin = fNumberOfRPsEBE+0.5; + //Printf("RP multiplicity: %lf",fNumberOfRPsEBE); + dMultiplicityBin = fNumberOfRPsEBE+0.5; } else if(fMultiplicityIs==AliFlowCommonConstants::kExternal) { + //Printf("Reference multiplicity: %lf",fReferenceMultiplicityEBE); dMultiplicityBin = fReferenceMultiplicityEBE+0.5; } else if(fMultiplicityIs==AliFlowCommonConstants::kPOI) { @@ -15820,7 +15822,7 @@ void AliFlowAnalysisWithQCumulants::CalculateIntFlowSumOfEventWeights() dMultiplicityBin = fNumberOfRPsEBE+0.5; } else if(fMultiplicityIs==AliFlowCommonConstants::kExternal) { - dMultiplicityBin = fReferenceMultiplicityEBE+0.5; + dMultiplicityBin = fReferenceMultiplicityEBE+0.5; } else if(fMultiplicityIs==AliFlowCommonConstants::kPOI) { dMultiplicityBin = fNumberOfPOIsEBE+0.5; diff --git a/PWG/FLOW/Base/AliFlowEventSimpleMakerOnTheFly.cxx b/PWG/FLOW/Base/AliFlowEventSimpleMakerOnTheFly.cxx index 3cb69957718..7be4c3c589b 100644 --- a/PWG/FLOW/Base/AliFlowEventSimpleMakerOnTheFly.cxx +++ b/PWG/FLOW/Base/AliFlowEventSimpleMakerOnTheFly.cxx @@ -24,6 +24,7 @@ #include "Riostream.h" #include "TMath.h" #include "TF1.h" +#include "TH3.h" #include "TRandom3.h" #include "AliFlowEventSimpleMakerOnTheFly.h" #include "AliFlowEventSimple.h" @@ -68,7 +69,13 @@ fProbability1(0.), fPhiMin2(0.), fPhiMax2(0.), fProbability2(0.), -fPi(TMath::Pi()) +fPi(TMath::Pi()), +fSimulateDetectorEffects(kFALSE), +fNumberOfInefficientSectors(0), +fInefficiencyFactorInPhi(1.0), +fNumberOfDeadSectors(0), +fEfficiencyDropNearEtaEdges(kFALSE), +fEfficiencyMatrix(0) { // Constructor. @@ -127,8 +134,75 @@ void AliFlowEventSimpleMakerOnTheFly::Init() fPhiDistribution->SetParName(6,"Hexagonal Flow (v6)"); fPhiDistribution->SetParameter(6,fV6); + //==============Efficiency matrix==============// + if(fSimulateDetectorEffects) SetupEfficiencyMatrix(); + //==============Efficiency matrix==============// + } // end of void AliFlowEventSimpleMakerOnTheFly::Init() +//________________________________________________________________________ +void AliFlowEventSimpleMakerOnTheFly::SetupEfficiencyMatrix() { + //Setup the efficiency matrix + TH1F *hPt = new TH1F("hPt","",200,0.1,20.1); + TH1F *hEta = new TH1F("hEta","",20,-0.95,0.95); + TH1F *hPhi = new TH1F("hPhi","",72,0.,2.*TMath::Pi()); + fEfficiencyMatrix = new TH3F("fEfficiencyMatrix","", + hEta->GetNbinsX(), + hEta->GetXaxis()->GetXmin(), + hEta->GetXaxis()->GetXmax(), + hPt->GetNbinsX(), + hPt->GetXaxis()->GetXmin(), + hPt->GetXaxis()->GetXmax(), + hPhi->GetNbinsX(), + hPhi->GetXaxis()->GetXmin(), + hPhi->GetXaxis()->GetXmax()); + + //Efficiency in pt + Double_t epsilon[20] = {0.3,0.6,0.77,0.79,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80,0.80}; + for(Int_t i=1;i<=20;i++) { + hPt->SetBinContent(i,epsilon[i-1]); + hPt->SetBinError(i,0.01); + } + for(Int_t i=21;i<=200;i++) { + hPt->SetBinContent(i,epsilon[19]); + hPt->SetBinError(i,0.01); + } + + //Efficiency in eta + for(Int_t i=1;i<=hEta->GetNbinsX();i++) { + hEta->SetBinContent(i,1.0); + hEta->SetBinError(i,0.01); + } + if(fEfficiencyDropNearEtaEdges) { + hEta->SetBinContent(1,0.7); hEta->SetBinContent(2,0.8); + hEta->SetBinContent(3,0.9); + hEta->SetBinContent(18,0.9); hEta->SetBinContent(19,0.8); + hEta->SetBinContent(20,0.7); + } + + //Efficiency in phi + for(Int_t i=1;i<=hPhi->GetNbinsX();i++) { + hPhi->SetBinContent(i,1.0); + hPhi->SetBinError(i,0.01); + } + for(Int_t i=1;i<=fNumberOfInefficientSectors;i++) + hPhi->SetBinContent(hPhi->FindBin(hPhi->GetRandom()),fInefficiencyFactorInPhi); + for(Int_t i=1;i<=fNumberOfDeadSectors;i++) + hPhi->SetBinContent(hPhi->FindBin(hPhi->GetRandom()),0.0); + + //Fill the 3D efficiency map + for(Int_t iBinX = 1; iBinX<=fEfficiencyMatrix->GetXaxis()->GetNbins();iBinX++) { + //cout<<"==================================="<GetYaxis()->GetNbins();iBinY++) { + //cout<<"==================================="<GetZaxis()->GetNbins();iBinZ++) { + fEfficiencyMatrix->SetBinContent(iBinX,iBinY,iBinZ,hEta->GetBinContent(iBinX)*hPt->GetBinContent(iBinY)*hPhi->GetBinContent(iBinZ)); + //cout<<"Eta: "<GetBinCenter(iBinX)<<" - Pt: "<GetBinCenter(iBinY)<<" - Phi: "<GetBinCenter(iBinZ)<<" - "<GetBinContent(iBinX)<<" , "<GetBinContent(iBinY)<<" , "<GetBinContent(iBinZ)<<" - Efficiency: "<GetBinContent(iBinX)*hPt->GetBinContent(iBinY)*hPhi->GetBinContent(iBinZ)<SetCharge((gRandom->Integer(2)>0.5 ? 1 : -1)); // Check uniform acceptance: if(!fUniformAcceptance && !AcceptOrNot(pTrack)){continue;} + //Detector effects + if(fSimulateDetectorEffects) { + Double_t randomNumber = gRandom->Rndm(); + if(randomNumber > fEfficiencyMatrix->GetBinContent(fEfficiencyMatrix->FindBin(pTrack->Eta(),pTrack->Pt(),pTrack->Phi()))) + continue; + } // Checking the RP cuts: if(cutsRP->PassesCuts(pTrack)) { diff --git a/PWG/FLOW/Base/AliFlowEventSimpleMakerOnTheFly.h b/PWG/FLOW/Base/AliFlowEventSimpleMakerOnTheFly.h index c8a398afb80..6edff07f0ec 100644 --- a/PWG/FLOW/Base/AliFlowEventSimpleMakerOnTheFly.h +++ b/PWG/FLOW/Base/AliFlowEventSimpleMakerOnTheFly.h @@ -17,6 +17,7 @@ class TF1; class TRandom3; +class TH3F; class AliFlowEventSimple; class AliFlowTrackSimple; @@ -80,9 +81,25 @@ class AliFlowEventSimpleMakerOnTheFly{ Double_t GetSecondSectorPhiMax() const {return this->fPhiMax2;} void SetSecondSectorProbability(Double_t dProbability2) {this->fProbability2 = dProbability2;} Double_t GetSecondSectorProbability() const {return this->fProbability2;} + + //Acceptance - simulate detector effects/inefficiencies + void SimulateDetectorEffects() {fSimulateDetectorEffects = kTRUE;} + void SetNumberOfInefficientSectorsInPhi(Int_t numberOfInefficientSectors) { + fNumberOfInefficientSectors = numberOfInefficientSectors; + fInefficiencyFactorInPhi = 0.5;} + void SetInefficiencyFactor(Double_t gInefficiencyFactorInPhi) { + fInefficiencyFactorInPhi = gInefficiencyFactorInPhi;} + void SetNumberOfDeadSectorsInPhi(Int_t numberOfDeadSectors) { + fNumberOfDeadSectors = numberOfDeadSectors;} + void EnableEfficiencyDropNearEtaEdges() { + fEfficiencyDropNearEtaEdges = kTRUE;} + private: AliFlowEventSimpleMakerOnTheFly(const AliFlowEventSimpleMakerOnTheFly& anAnalysis); // copy constructor AliFlowEventSimpleMakerOnTheFly& operator=(const AliFlowEventSimpleMakerOnTheFly& anAnalysis); // assignment operator + + void SetupEfficiencyMatrix(); + Int_t fCount; // count number of events Int_t fMinMult; // uniformly sampled multiplicity is >= iMinMult Int_t fMaxMult; // uniformly sampled multiplicity is < iMaxMult @@ -116,6 +133,14 @@ class AliFlowEventSimpleMakerOnTheFly{ Double_t fProbability2; // particles emitted in fPhiMin2 < phi < fPhiMax2 are taken with probability fProbability2 Double_t fPi; // pi + //Simulate detector effects + Bool_t fSimulateDetectorEffects;//simulate detector effects in pT + Int_t fNumberOfInefficientSectors;//inefficient secotrs in phi + Double_t fInefficiencyFactorInPhi;//efficiency factor < 1 + Int_t fNumberOfDeadSectors;//number of dead sectors + Bool_t fEfficiencyDropNearEtaEdges;//efficiency drop in eta edges + TH3F *fEfficiencyMatrix; //efficiency matrix in eta-pt-phi + ClassDef(AliFlowEventSimpleMakerOnTheFly,1) // macro for rootcint }; diff --git a/PWG/FLOW/Tasks/AliAnalysisTaskFlowEvent.cxx b/PWG/FLOW/Tasks/AliAnalysisTaskFlowEvent.cxx index efb36b993b3..3643f8ebb20 100644 --- a/PWG/FLOW/Tasks/AliAnalysisTaskFlowEvent.cxx +++ b/PWG/FLOW/Tasks/AliAnalysisTaskFlowEvent.cxx @@ -334,7 +334,7 @@ void AliAnalysisTaskFlowEvent::UserExec(Option_t *) //fFlowEvent = new AliFlowEvent( fCutsRP, fCutsPOI ); // if (myESD) - fFlowEvent->SetReferenceMultiplicity(fCutsEvent->GetReferenceMultiplicity(InputEvent())); + fFlowEvent->SetReferenceMultiplicity(fCutsEvent->GetReferenceMultiplicity(InputEvent(),mcEvent)); if (mcEvent && mcEvent->GenEventHeader()) fFlowEvent->SetMCReactionPlaneAngle(mcEvent); } diff --git a/PWG/FLOW/Tasks/AliFlowEvent.cxx b/PWG/FLOW/Tasks/AliFlowEvent.cxx index c5c9f1620a4..eca9e523695 100644 --- a/PWG/FLOW/Tasks/AliFlowEvent.cxx +++ b/PWG/FLOW/Tasks/AliFlowEvent.cxx @@ -635,7 +635,8 @@ void AliFlowEvent::Fill( AliFlowTrackCuts* rpCuts, if (sourceRP==sourcePOI) { //loop over tracks - for (Int_t i=0; iGetNumberOfInputObjects(); i++) + Int_t numberOfInputObject = rpCuts->GetNumberOfInputObjects(); + for (Int_t i=0; iGetInputObject(i); @@ -667,8 +668,9 @@ void AliFlowEvent::Fill( AliFlowTrackCuts* rpCuts, //here we have two different sources of particles, so we fill //them independently //RP - for (Int_t i=0; iGetNumberOfInputObjects(); i++) - { + Int_t numberOfInputObject = rpCuts->GetNumberOfInputObjects(); + for (Int_t i=0; iGetInputObject(i); Bool_t rp = rpCuts->IsSelected(particle,i); if (!rp) continue; @@ -679,7 +681,8 @@ void AliFlowEvent::Fill( AliFlowTrackCuts* rpCuts, fNumberOfTracks++; } //POI - for (Int_t i=0; iGetNumberOfInputObjects(); i++) + numberOfInputObject = poiCuts->GetNumberOfInputObjects(); + for (Int_t i=0; iGetInputObject(i); Bool_t poi = poiCuts->IsSelected(particle,i); @@ -692,6 +695,7 @@ void AliFlowEvent::Fill( AliFlowTrackCuts* rpCuts, } } } + //----------------------------------------------------------------------- void AliFlowEvent::InsertTrack(AliFlowTrack *thisTrack) { // adds a flow track at the end of the container @@ -748,7 +752,8 @@ AliFlowEvent::AliFlowEvent( AliFlowTrackCuts* rpCuts, if (sourceRP==sourcePOI) { //loop over tracks - for (Int_t i=0; iGetNumberOfInputObjects(); i++) + Int_t numberOfInputObject = rpCuts->GetNumberOfInputObjects(); + for (Int_t i=0; iGetInputObject(i); @@ -783,7 +788,8 @@ AliFlowEvent::AliFlowEvent( AliFlowTrackCuts* rpCuts, //them independently AliFlowTrack* pTrack = NULL; //RP - for (Int_t i=0; iGetNumberOfInputObjects(); i++) + Int_t numberOfInputObject = rpCuts->GetNumberOfInputObjects(); + for (Int_t i=0; iGetInputObject(i); Bool_t rp = rpCuts->IsSelected(particle,i); @@ -794,7 +800,8 @@ AliFlowEvent::AliFlowEvent( AliFlowTrackCuts* rpCuts, AddTrack(pTrack); } //POI - for (Int_t i=0; iGetNumberOfInputObjects(); i++) + numberOfInputObject = poiCuts->GetNumberOfInputObjects(); + for (Int_t i=0; iGetInputObject(i); Bool_t poi = poiCuts->IsSelected(particle,i); diff --git a/PWG/FLOW/Tasks/AliFlowEventCuts.cxx b/PWG/FLOW/Tasks/AliFlowEventCuts.cxx index 3efa518eb54..f03b7a85388 100644 --- a/PWG/FLOW/Tasks/AliFlowEventCuts.cxx +++ b/PWG/FLOW/Tasks/AliFlowEventCuts.cxx @@ -368,7 +368,7 @@ Bool_t AliFlowEventCuts::PassesCuts(AliVEvent *event) if(fCutRefMult&&esdevent) { //reference multiplicity still to be defined - Double_t refMult = RefMult(event); + Double_t refMult = RefMult(event,0x0); if (refMult < fRefMultMin || refMult >= fRefMultMax ) { pass=kFALSE; @@ -461,7 +461,7 @@ AliFlowEventCuts* AliFlowEventCuts::StandardCuts() } //----------------------------------------------------------------------- -Int_t AliFlowEventCuts::RefMult(AliVEvent* event) +Int_t AliFlowEventCuts::RefMult(AliVEvent* event, AliMCEvent *mcEvent) { //calculate the reference multiplicity, if all fails return 0 AliESDVZERO* vzero = NULL; @@ -501,9 +501,9 @@ Int_t AliFlowEventCuts::RefMult(AliVEvent* event) } Int_t refmult=0; - fRefMultCuts->SetEvent(event); - for (Int_t i=0; iGetNumberOfInputObjects(); i++) - { + fRefMultCuts->SetEvent(event,mcEvent); + Int_t numberOfInputObjects = fRefMultCuts->GetNumberOfInputObjects(); + for (Int_t i=0; iIsSelected(fRefMultCuts->GetInputObject(i),i)) refmult++; } diff --git a/PWG/FLOW/Tasks/AliFlowEventCuts.h b/PWG/FLOW/Tasks/AliFlowEventCuts.h index 76fdb505a71..2d76801f556 100644 --- a/PWG/FLOW/Tasks/AliFlowEventCuts.h +++ b/PWG/FLOW/Tasks/AliFlowEventCuts.h @@ -14,6 +14,7 @@ #include "TNamed.h" class AliVEvent; +class AliMCEvent; class TBrowser; #include "TList.h" #include "TH1.h" @@ -74,9 +75,9 @@ class AliFlowEventCuts : public TNamed { TH1* QAbefore(Int_t i) {return static_cast(static_cast(fQA->At(0))->At(i));} TH1* QAafter(Int_t i) {return static_cast(static_cast(fQA->At(1))->At(i));} - Int_t RefMult(AliVEvent* event); + Int_t RefMult(AliVEvent* event, AliMCEvent *mcEvent = 0x0); //Int_t GetRefMult() {return fRefMult;} - Int_t GetReferenceMultiplicity(AliVEvent* event) {return RefMult(event);} + Int_t GetReferenceMultiplicity(AliVEvent* event, AliMCEvent *mcEvent) {return RefMult(event,mcEvent);} const char* CentrMethName(refMultMethod method) const; void SetCentralityPercentileRange(Float_t min, Float_t max){ fCentralityPercentileMin=min; fCentralityPercentileMax=max; diff --git a/PWG/FLOW/Tasks/AliFlowTrackCuts.cxx b/PWG/FLOW/Tasks/AliFlowTrackCuts.cxx index 159e39ef769..1ba66f3e344 100644 --- a/PWG/FLOW/Tasks/AliFlowTrackCuts.cxx +++ b/PWG/FLOW/Tasks/AliFlowTrackCuts.cxx @@ -485,6 +485,7 @@ Bool_t AliFlowTrackCuts::IsSelected(TObject* obj, Int_t id) if (fParamType==kMUON) return PassesMuonCuts(vparticle); // XZhang 20120604 return PassesCuts(vparticle); // XZhang 20120604 } // XZhang 20120604 + AliFlowTrackSimple* flowtrack = dynamic_cast(obj); if (flowtrack) return PassesCuts(flowtrack); AliMultiplicity* tracklets = dynamic_cast(obj); -- 2.43.0