patching flow cuts (Panos)

diff --git a/PWG/FLOW/Base/AliFlowAnalysisWithQCumulants.cxx b/PWG/FLOW/Base/AliFlowAnalysisWithQCumulants.cxx
index a1ce18a735b03377dc236310bd185d54597c08c5..120b3df10afce5bfaab5725dc58b60b1b5865beb 100644 (file)
--- 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 3cb699577187a8aef885bb3ee594dd7a0c174705..7be4c3c589b5e3ee227b5a1f73a8c07a191a1626 100644 (file)
--- 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<<"==================================="<<endl;
+    for(Int_t iBinY = 1; iBinY<=fEfficiencyMatrix->GetYaxis()->GetNbins();iBinY++) {
+      //cout<<"==================================="<<endl;
+      for(Int_t iBinZ = 1; iBinZ<=fEfficiencyMatrix->GetZaxis()->GetNbins();iBinZ++) {
+       fEfficiencyMatrix->SetBinContent(iBinX,iBinY,iBinZ,hEta->GetBinContent(iBinX)*hPt->GetBinContent(iBinY)*hPhi->GetBinContent(iBinZ));
+       //cout<<"Eta: "<<hEta->GetBinCenter(iBinX)<<" - Pt: "<<hPt->GetBinCenter(iBinY)<<" - Phi: "<<hPhi->GetBinCenter(iBinZ)<<" - "<<hEta->GetBinContent(iBinX)<<" , "<<hPt->GetBinContent(iBinY)<<" , "<<hPhi->GetBinContent(iBinZ)<<" - Efficiency: "<<hEta->GetBinContent(iBinX)*hPt->GetBinContent(iBinY)*hPhi->GetBinContent(iBinZ)<<endl;
+      }
+    }
+  }
+}
+
 //====================================================================================================================
 
 Bool_t AliFlowEventSimpleMakerOnTheFly::AcceptOrNot(AliFlowTrackSimple *pTrack)
@@ -197,6 +271,12 @@ AliFlowEventSimple* AliFlowEventSimpleMakerOnTheFly::CreateEventOnTheFly(AliFlow
   pTrack->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 c8a398afb80eaca19e813532c5f8453cc950cb79..6edff07f0ec117cb701fabce740d3ef0bff4c1d7 100644 (file)
--- 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 efb36b993b33537343d8917cafe284079740619d..3643f8ebb20f18ef960f49747f4594d3bca04a24 100644 (file)
--- 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 c5c9f1620a489dc5fdfe7982e6674f01c85d77c8..eca9e523695910f83b3a23dda35af37ed4696061 100644 (file)
--- a/PWG/FLOW/Tasks/AliFlowEvent.cxx
+++ b/PWG/FLOW/Tasks/AliFlowEvent.cxx
@@ -635,7 +635,8 @@ void AliFlowEvent::Fill( AliFlowTrackCuts* rpCuts,
   if (sourceRP==sourcePOI)\r
   {\r
     //loop over tracks\r
-    for (Int_t i=0; i<rpCuts->GetNumberOfInputObjects(); i++)\r
+    Int_t numberOfInputObject = rpCuts->GetNumberOfInputObjects();\r
+    for (Int_t i=0; i<numberOfInputObject; i++)\r
     {\r
       //get input object (particle)\r
       TObject* particle = rpCuts->GetInputObject(i);\r
@@ -667,8 +668,9 @@ void AliFlowEvent::Fill( AliFlowTrackCuts* rpCuts,
     //here we have two different sources of particles, so we fill\r
     //them independently\r
     //RP\r
-    for (Int_t i=0; i<rpCuts->GetNumberOfInputObjects(); i++)\r
-    {\r
+    Int_t numberOfInputObject = rpCuts->GetNumberOfInputObjects();\r
+    for (Int_t i=0; i<numberOfInputObject; i++)\r
+      {\r
       TObject* particle = rpCuts->GetInputObject(i);\r
       Bool_t rp = rpCuts->IsSelected(particle,i);\r
       if (!rp) continue;\r
@@ -679,7 +681,8 @@ void AliFlowEvent::Fill( AliFlowTrackCuts* rpCuts,
       fNumberOfTracks++;\r
     }\r
     //POI\r
-    for (Int_t i=0; i<poiCuts->GetNumberOfInputObjects(); i++)\r
+    numberOfInputObject = poiCuts->GetNumberOfInputObjects();\r
+    for (Int_t i=0; i<numberOfInputObject; i++)\r
     {\r
       TObject* particle = poiCuts->GetInputObject(i);\r
       Bool_t poi = poiCuts->IsSelected(particle,i);\r
@@ -692,6 +695,7 @@ void AliFlowEvent::Fill( AliFlowTrackCuts* rpCuts,
     }\r
   }\r
 }\r
+\r
 //-----------------------------------------------------------------------\r
 void AliFlowEvent::InsertTrack(AliFlowTrack *thisTrack) {\r
   // adds a flow track at the end of the container\r
@@ -748,7 +752,8 @@ AliFlowEvent::AliFlowEvent( AliFlowTrackCuts* rpCuts,
   if (sourceRP==sourcePOI)\r
   {\r
     //loop over tracks\r
-    for (Int_t i=0; i<rpCuts->GetNumberOfInputObjects(); i++)\r
+    Int_t numberOfInputObject = rpCuts->GetNumberOfInputObjects();\r
+    for (Int_t i=0; i<numberOfInputObject; i++)\r
     {\r
       //get input object (particle)\r
       TObject* particle = rpCuts->GetInputObject(i);\r
@@ -783,7 +788,8 @@ AliFlowEvent::AliFlowEvent( AliFlowTrackCuts* rpCuts,
     //them independently\r
     AliFlowTrack* pTrack = NULL;\r
     //RP\r
-    for (Int_t i=0; i<rpCuts->GetNumberOfInputObjects(); i++)\r
+    Int_t numberOfInputObject = rpCuts->GetNumberOfInputObjects();\r
+    for (Int_t i=0; i<numberOfInputObject; i++)\r
     {\r
       TObject* particle = rpCuts->GetInputObject(i);\r
       Bool_t rp = rpCuts->IsSelected(particle,i);\r
@@ -794,7 +800,8 @@ AliFlowEvent::AliFlowEvent( AliFlowTrackCuts* rpCuts,
       AddTrack(pTrack);\r
     }\r
     //POI\r
-    for (Int_t i=0; i<poiCuts->GetNumberOfInputObjects(); i++)\r
+    numberOfInputObject = poiCuts->GetNumberOfInputObjects();\r
+    for (Int_t i=0; i<numberOfInputObject; i++)\r
     {\r
       TObject* particle = poiCuts->GetInputObject(i);\r
       Bool_t poi = poiCuts->IsSelected(particle,i);\r

diff --git a/PWG/FLOW/Tasks/AliFlowEventCuts.cxx b/PWG/FLOW/Tasks/AliFlowEventCuts.cxx
index 3efa518eb5472bc58b1f11d7091d9a563d48069f..f03b7a853886abc569a8e5186e9d744a8c9e3cac 100644 (file)
--- 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; i<fRefMultCuts->GetNumberOfInputObjects(); i++)
-  {
+  fRefMultCuts->SetEvent(event,mcEvent);
+  Int_t numberOfInputObjects = fRefMultCuts->GetNumberOfInputObjects();
+  for (Int_t i=0; i<numberOfInputObjects; i++) {
     if (fRefMultCuts->IsSelected(fRefMultCuts->GetInputObject(i),i))
       refmult++;
   }

diff --git a/PWG/FLOW/Tasks/AliFlowEventCuts.h b/PWG/FLOW/Tasks/AliFlowEventCuts.h
index 76fdb505a719ce6c5a631edf6b8603f4cce74ee2..2d76801f5566324ac664464ac41c2e7662f05f0f 100644 (file)
--- 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<TH1*>(static_cast<TList*>(fQA->At(0))->At(i));}
   TH1* QAafter(Int_t i) {return static_cast<TH1*>(static_cast<TList*>(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 159e39ef769037e8212d4cf2bdf531c763b42b82..1ba66f3e344b5e54a5da45a6cf48632ad52f2d72 100644 (file)
--- 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<AliFlowTrackSimple*>(obj);
   if (flowtrack) return PassesCuts(flowtrack);
   AliMultiplicity* tracklets = dynamic_cast<AliMultiplicity*>(obj);