Event Cut update + outlier centrality for flow
authorrbailhac <rbailhac@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 13 Feb 2013 13:10:39 +0000 (13:10 +0000)
committerrbailhac <rbailhac@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 13 Feb 2013 13:10:39 +0000 (13:10 +0000)
PWGHF/hfe/AliAnalysisTaskFlowTPCTOFEPSP.cxx
PWGHF/hfe/AliAnalysisTaskFlowTPCTOFEPSP.h
PWGHF/hfe/AliHFEcuts.cxx
PWGHF/hfe/AliHFEcuts.h
PWGHF/hfe/AliHFEextraEventCuts.cxx
PWGHF/hfe/AliHFEextraEventCuts.h

index 3af8721..a842e08 100644 (file)
@@ -147,8 +147,7 @@ AliAnalysisTaskFlowTPCTOFEPSP::AliAnalysisTaskFlowTPCTOFEPSP() :
   fHFEVZEROEventPlane(0x0),\r
   fHistEV(0),\r
   fHistPileUp(0),\r
-  fUpperPileUpCut(0),\r
-  fLowerPileUpCut(0),\r
+  fPileUpCut(kFALSE),\r
   fEventPlane(0x0),\r
   fEventPlaneaftersubtraction(0x0),\r
   fFractionContamination(0x0),\r
@@ -262,8 +261,7 @@ AliAnalysisTaskFlowTPCTOFEPSP:: AliAnalysisTaskFlowTPCTOFEPSP(const char *name)
   fHFEVZEROEventPlane(0x0),\r
   fHistEV(0),\r
   fHistPileUp(0),\r
-  fUpperPileUpCut(0),\r
-  fLowerPileUpCut(0),\r
+  fPileUpCut(kFALSE),\r
   fEventPlane(0x0),\r
   fEventPlaneaftersubtraction(0x0),\r
   fFractionContamination(0x0),\r
@@ -399,8 +397,7 @@ AliAnalysisTaskFlowTPCTOFEPSP::AliAnalysisTaskFlowTPCTOFEPSP(const AliAnalysisTa
   fHFEVZEROEventPlane(NULL),\r
   fHistEV(NULL),\r
   fHistPileUp(NULL),\r
-  fUpperPileUpCut(NULL),\r
-  fLowerPileUpCut(NULL),\r
+  fPileUpCut(kFALSE),\r
   fEventPlane(NULL),\r
   fEventPlaneaftersubtraction(NULL),\r
   fFractionContamination(NULL),\r
@@ -461,7 +458,6 @@ void AliAnalysisTaskFlowTPCTOFEPSP::Copy(TObject &o) const {
   // Copy into object o\r
   //\r
   AliAnalysisTaskFlowTPCTOFEPSP &target = dynamic_cast<AliAnalysisTaskFlowTPCTOFEPSP &>(o);\r
-  target.fListHist = fListHist;\r
   target.fAODAnalysis = fAODAnalysis;\r
   target.fUseFilterAOD = fUseFilterAOD;\r
   target.fApplyCut = fApplyCut;\r
@@ -498,6 +494,7 @@ void AliAnalysisTaskFlowTPCTOFEPSP::Copy(TObject &o) const {
   target.fMaxopening3D = fMaxopening3D;\r
   target.fMaxInvmass = fMaxInvmass;\r
   target.fSetMassConstraint =  fSetMassConstraint;\r
+  target.fAlgorithmMA = fAlgorithmMA;\r
   target.fCounterPoolBackground = fCounterPoolBackground;\r
   target.fDebugLevel = fDebugLevel;\r
   target.fMonitorEventPlane = fMonitorEventPlane;\r
@@ -511,58 +508,8 @@ void AliAnalysisTaskFlowTPCTOFEPSP::Copy(TObject &o) const {
   target.fHFECuts = fHFECuts;\r
   target.fRejectKinkMother = fRejectKinkMother;\r
   target.fPID = fPID;\r
-  target.fPIDTOFOnly = fPIDTOFOnly;\r
   target.fPIDqa = fPIDqa;\r
-  target.fflowEvent = fflowEvent;\r
-  target.fHFEBackgroundCuts = fHFEBackgroundCuts;\r
-  target.fPIDBackground = fPIDBackground;\r
-  target.fPIDBackgroundqa = fPIDBackgroundqa;\r
-  target.fAlgorithmMA = fAlgorithmMA;\r
-  target.fArraytrack = fArraytrack;\r
-  target.fCounterPoolBackground = fCounterPoolBackground;\r
   target.fHFEVZEROEventPlane = fHFEVZEROEventPlane;\r
-  target.fHistEV=fHistEV;      \r
-  target.fHistPileUp=fHistPileUp;  \r
-  target.fUpperPileUpCut=fUpperPileUpCut;        \r
-  target.fLowerPileUpCut=fLowerPileUpCut;        \r
-  target.fEventPlane=fEventPlane;     \r
-  target.fEventPlaneaftersubtraction=fEventPlaneaftersubtraction; \r
-  target.fFractionContamination=fFractionContamination;    \r
-  target.fContaminationv2=fContaminationv2;          \r
-  target.fCosSin2phiep=fCosSin2phiep;        \r
-  target.fCos2phie=fCos2phie;  \r
-  target.fSin2phie=fSin2phie;  \r
-  target.fCos2phiep=fCos2phiep; \r
-  target.fSin2phiep=fSin2phiep; \r
-  target.fSin2phiephiep=fSin2phiephiep;  \r
-  target.fCosResabc=fCosResabc; \r
-  target.fSinResabc=fSinResabc; \r
-  target.fProfileCosResab=fProfileCosResab; \r
-  target.fProfileCosResac=fProfileCosResac; \r
-  target.fProfileCosResbc=fProfileCosResbc; \r
-  target.fCosRes=fCosRes; \r
-  target.fSinRes=fSinRes; \r
-  target.fProfileCosRes=fProfileCosRes; \r
-  target.fTrackingCuts=fTrackingCuts; \r
-  target.fDeltaPhiMapsBeforePID=fDeltaPhiMapsBeforePID; \r
-  target.fCosPhiMapsBeforePID=fCosPhiMapsBeforePID; \r
-  target.fDeltaPhiMaps=fDeltaPhiMaps; \r
-  target.fDeltaPhiMapsContamination=fDeltaPhiMapsContamination; \r
-  target.fCosPhiMaps=fCosPhiMaps;         \r
-  target.fProfileCosPhiMaps=fProfileCosPhiMaps;  \r
-  target.fDeltaPhiMapsTaggedPhotonic=fDeltaPhiMapsTaggedPhotonic; \r
-  target.fDeltaPhiMapsTaggedNonPhotonic=fDeltaPhiMapsTaggedNonPhotonic; \r
-  target.fDeltaPhiMapsTaggedPhotonicLS=fDeltaPhiMapsTaggedPhotonicLS; \r
-  target.fMCSourceDeltaPhiMaps=fMCSourceDeltaPhiMaps; \r
-  target.fOppSignDeltaPhiMaps=fOppSignDeltaPhiMaps;  \r
-  target.fSameSignDeltaPhiMaps=fSameSignDeltaPhiMaps; \r
-  target.fOppSignAngle=fOppSignAngle;         \r
-  target.fSameSignAngle=fSameSignAngle;   \r
\r
-  for(Int_t k = 0; k < 10; k++) {\r
-    target.fBinCentralityLess[k] = fBinCentralityLess[k];\r
-  }     \r
-\r
   for(Int_t k = 0; k < 11; k++) {\r
     target.fContamination[k] = fContamination[k];\r
     target.fv2contamination[k] = fv2contamination[k];\r
@@ -807,6 +754,12 @@ void AliAnalysisTaskFlowTPCTOFEPSP::UserCreateOutputObjects()
   Double_t binLimCMore[nBinsCMore+1];\r
   for(Int_t i=0; i<=nBinsCMore; i++) binLimCMore[i]=(Double_t)minCMore + (maxCMore-minCMore)/nBinsCMore*(Double_t)i ;\r
 \r
+  Int_t nBinsCEvenMore = 100;\r
+  Double_t minCEvenMore = 0.0;\r
+  Double_t maxCEvenMore = 100.0;\r
+  Double_t binLimCEvenMore[nBinsCEvenMore+1];\r
+  for(Int_t i=0; i<=nBinsCEvenMore; i++) binLimCEvenMore[i]=(Double_t)minCEvenMore + (maxCEvenMore-minCEvenMore)/nBinsCEvenMore*(Double_t)i ;\r
+\r
   Int_t nBinsPhi = 8;\r
   Double_t minPhi = 0.0;\r
   Double_t maxPhi = TMath::Pi();\r
@@ -861,7 +814,7 @@ void AliAnalysisTaskFlowTPCTOFEPSP::UserCreateOutputObjects()
 \r
   Int_t nBinsMult = 100;\r
   Double_t minMult = 0.;\r
-  Double_t maxMult = 25000;\r
+  Double_t maxMult = 3000;\r
   Double_t binLimMult[nBinsMult+1];\r
   //for(Int_t i=0; i<=nBinsMult; i++) binLimMult[i]=TMath::Power((Double_t)minMult + (TMath::Sqrt(maxMult)-TMath::Sqrt(minMult))/nBinsMult*(Double_t)i,2);\r
   for(Int_t i=0; i<=nBinsMult; i++) binLimMult[i]=(Double_t)minMult + (maxMult-minMult)/nBinsMult*(Double_t)i;\r
@@ -885,12 +838,13 @@ void AliAnalysisTaskFlowTPCTOFEPSP::UserCreateOutputObjects()
   AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: histev");\r
 \r
   // V0 multiplicity vs # of tracks vs centraliy\r
-  const Int_t nDimPU=3;\r
-  Int_t nBinPU[nDimPU] = {nBinsMult,nBinsMult,nBinsCMore};\r
+  const Int_t nDimPU=4;\r
+  Int_t nBinPU[nDimPU] = {nBinsCEvenMore,nBinsCEvenMore,nBinsMult,nBinsMult};\r
   fHistPileUp = new THnSparseF("PileUp","PileUp",nDimPU,nBinPU);\r
-  fHistPileUp->SetBinEdges(0,binLimMult);\r
-  fHistPileUp->SetBinEdges(1,binLimMult);\r
-  fHistPileUp->SetBinEdges(2,binLimCMore);\r
+  fHistPileUp->SetBinEdges(0,binLimCEvenMore);\r
+  fHistPileUp->SetBinEdges(1,binLimCEvenMore);\r
+  fHistPileUp->SetBinEdges(2,binLimMult);\r
+  fHistPileUp->SetBinEdges(3,binLimMult);\r
   fHistPileUp->Sumw2();\r
   \r
   AliDebug(2,"AliAnalysisTaskFlowTPCTOFEPSP: eventplane");\r
@@ -1514,26 +1468,76 @@ void AliAnalysisTaskFlowTPCTOFEPSP::UserExec(Option_t */*option*/)
   ///////////////////////////////////////////////////////////\r
   // PileUpCut\r
   ///////////////////////////////////////////////////////////\r
\r
-  AliVVZERO* vzeroData=fInputEvent->GetVZEROData();\r
-  Double_t mult[3],multV0A(0),multV0C(0);\r
-  for(Int_t i=0; i<32; ++i) {\r
-    multV0A += vzeroData->GetMultiplicityV0A(i);\r
-    multV0C += vzeroData->GetMultiplicityV0C(i);\r
-  }\r
-  mult[0]=fInputEvent->GetNumberOfTracks();\r
-  mult[1]=multV0A+multV0C;\r
-  mult[2]=binctMore;\r
-  fHistPileUp->Fill(mult);\r
-\r
-  if(fUpperPileUpCut&&fLowerPileUpCut){\r
-    if((mult[0]<fLowerPileUpCut->Eval(mult[1])) || \r
-       (mult[0]>fUpperPileUpCut->Eval(mult[1]))){\r
-      AliDebug(2,"Does not pass the pileup cut");\r
-      PostData(1, fListHist);\r
+\r
+  Float_t multTPC(0.); // tpc mult estimate\r
+  Float_t multGlob(0.); // global multiplicity\r
+  const Int_t nGoodTracks = fInputEvent->GetNumberOfTracks();\r
+  for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill tpc mult\r
+    AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*>(fInputEvent->GetTrack(iTracks));\r
+    if (!trackAOD) continue;\r
+    if (!(trackAOD->TestFilterBit(1))) continue;\r
+    if ((trackAOD->Pt() < .2) || (trackAOD->Pt() > 5.0) || (TMath::Abs(trackAOD->Eta()) > .8) || (trackAOD->GetTPCNcls() < 70)  || (trackAOD->GetDetPid()->GetTPCsignal() < 10.0) || (trackAOD->Chi2perNDF() < 0.2)) continue;\r
+    multTPC++;\r
+  }\r
+  for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill global mult\r
+    AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*>(fInputEvent->GetTrack(iTracks));\r
+    if (!trackAOD) continue;\r
+    if (!(trackAOD->TestFilterBit(16))) continue;\r
+    if ((trackAOD->Pt() < .2) || (trackAOD->Pt() > 5.0) || (TMath::Abs(trackAOD->Eta()) > .8) || (trackAOD->GetTPCNcls() < 70) || (trackAOD->GetDetPid()->GetTPCsignal() < 10.0) || (trackAOD->Chi2perNDF() < 0.1)) continue;\r
+    Double_t b[2] = {-99., -99.};\r
+    Double_t bCov[3] = {-99., -99., -99.};\r
+    if (!(trackAOD->PropagateToDCA(fInputEvent->GetPrimaryVertex(), fInputEvent->GetMagneticField(), 100., b, bCov))) continue;\r
+    if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;\r
+    multGlob++;\r
+  } //track loop\r
+\r
+  Double_t pileup[4];\r
+  pileup[0]=fInputEvent->GetCentrality()->GetCentralityPercentile("V0M");\r
+  pileup[1]=fInputEvent->GetCentrality()->GetCentralityPercentile("TRK");\r
+  pileup[2]=multTPC;\r
+  pileup[3]=multGlob;\r
+  fHistPileUp->Fill(pileup);\r
+\r
+  if(fPileUpCut){\r
+    if (TMath::Abs(pileup[0]-pileup[1]) > 5) {\r
+      AliDebug(2,"Does not pass the centrality correlation cut");\r
+      return;\r
+    }\r
+    if(multTPC < (-36.81+1.48*multGlob) && multTPC > (63.03+1.78*multGlob)){\r
+      AliDebug(2,"Does not pass the multiplicity correlation cut");\r
       return;\r
     }\r
   }\r
\r
+  // AliVVZERO* vzeroData=fInputEvent->GetVZEROData();\r
+  // Double_t mult[3],multV0A(0),multV0C(0);\r
+  // for(Int_t i=0; i<32; ++i) {\r
+  //   multV0A += vzeroData->GetMultiplicityV0A(i);\r
+  //   multV0C += vzeroData->GetMultiplicityV0C(i);\r
+  // }\r
+\r
+  // int ntrk=0;\r
+  // for(Int_t k = 0; k < fInputEvent->GetNumberOfTracks(); k++){\r
+  //   AliVTrack *track = (AliVTrack *) fInputEvent->GetTrack(k);\r
+  //   if(!track) continue;\r
+  //   if(!(track->GetStatus()&AliVTrack::kITSrefit)) continue;\r
+  //   if(!(track->GetStatus()&AliVTrack::kTPCrefit)) continue;\r
+  //   ntrk++;\r
+  // }\r
+    \r
+  // mult[0]=fInputEvent->GetNumberOfTracks();\r
+  // mult[1]=multV0A+multV0C;\r
+  // mult[2]=binctMore;\r
+  // fHistPileUp->Fill(mult);\r
+\r
+  // if(fUpperPileUpCut&&fLowerPileUpCut){\r
+  //   if((mult[0]<fLowerPileUpCut->Eval(mult[1])) || \r
+  //      (mult[0]>fUpperPileUpCut->Eval(mult[1]))){\r
+  //     AliDebug(2,"Does not pass the pileup cut");\r
+  //     PostData(1, fListHist);\r
+  //     return;\r
+  //   }\r
+  // }\r
 \r
   ////////////////////////////////////  \r
   // First method event plane\r
@@ -1581,7 +1585,7 @@ void AliAnalysisTaskFlowTPCTOFEPSP::UserExec(Option_t */*option*/)
   else {\r
     \r
     Double_t qVx, qVy;  //TR: info\r
-    eventPlaneV0 = vEPa->CalculateVZEROEventPlane(fInputEvent,10,2,qVx,qVy);\r
+    eventPlaneV0 = TVector2::Phi_0_2pi(vEPa->CalculateVZEROEventPlane(fInputEvent,10,2,qVx,qVy));\r
     if(eventPlaneV0 > TMath::Pi()) eventPlaneV0 = eventPlaneV0 - TMath::Pi();\r
     qV0.Set(qVx,qVy);\r
     eventPlaneV0A = TVector2::Phi_0_2pi(vEPa->CalculateVZEROEventPlane(fInputEvent,8,2,qVx,qVy));\r
@@ -1596,8 +1600,8 @@ void AliAnalysisTaskFlowTPCTOFEPSP::UserExec(Option_t */*option*/)
     if(eventPlaneV0C<-900) return;\r
 \r
     eventPlaneV0=TVector2::Phi_0_2pi(eventPlaneV0);\r
-    eventPlaneV0A=TVector2::Phi_0_2pi(eventPlaneV0);\r
-    eventPlaneV0C=TVector2::Phi_0_2pi(eventPlaneV0);\r
+    eventPlaneV0A=TVector2::Phi_0_2pi(eventPlaneV0A);\r
+    eventPlaneV0C=TVector2::Phi_0_2pi(eventPlaneV0C);\r
   }\r
 \r
 \r
@@ -1784,6 +1788,7 @@ void AliAnalysisTaskFlowTPCTOFEPSP::UserExec(Option_t */*option*/)
     valuensparsea[1] = eventPlanesub1a;\r
     valuensparsea[2] = eventPlanesub2a;\r
   } \r
+  //printf("%f %f %f\n",valuensparsea[0],valuensparsea[1],valuensparsea[2]);\r
   fEventPlane->Fill(&valuensparsea[0]);\r
 \r
   // Fill\r
index 7f3d95c..274991e 100644 (file)
@@ -146,7 +146,7 @@ public:
   void  SetMaxopeningphi(Double_t maxOpeningphi) { fMaxopeningphi = maxOpeningphi; };\r
   void  SetAlgorithmMA(Bool_t algorithmMA) { fAlgorithmMA = algorithmMA; };\r
   void  SetMassConstraint(Bool_t massConstraint) { fSetMassConstraint = massConstraint; };\r
-  void  SetPileUpCuts(TF1 *up, TF1 *lo) { fUpperPileUpCut=up; fLowerPileUpCut=lo; }\r
+  void  SetPileUpCut(Bool_t cut=kTRUE) { fPileUpCut=cut; }\r
 \r
   Int_t    LookAtNonHFE(Int_t iTrack1, AliVTrack *track1, AliVEvent *fESD, AliMCEvent *mcEvent,Int_t binct,Double_t deltaphi,Int_t source,Int_t indexmother);\r
   \r
@@ -236,9 +236,8 @@ private:
   // Histos\r
   TH2D *fHistEV;               //! Number of events\r
   THnSparseF *fHistPileUp;     //! Pile up histogram\r
-  TF1 *fUpperPileUpCut;        //! Pile up cut\r
-  TF1 *fLowerPileUpCut;        //! Pile up cut\r
-  \r
+  Bool_t fPileUpCut;\r
+\r
   // A Event plane as function of phiepa, phiepb, phiepc, phiepd centrality \r
   // a V0A, b V0C, c TPC,\r
   THnSparseF *fEventPlane;     //! Event plane\r
index c9b8f9e..51ad106 100644 (file)
@@ -135,6 +135,7 @@ AliHFEcuts::AliHFEcuts():
   fUseMixedVertex(kTRUE),
   fUseSPDVertex(kFALSE),
   fUseCorrelationVertex(kFALSE),
+  fSPDVtxResolution(kFALSE),
   fIsIPSigmacut(kFALSE),
   fIsIPcharge(kFALSE),
   fIsIPOpp(kFALSE),
@@ -185,6 +186,7 @@ AliHFEcuts::AliHFEcuts(const Char_t *name, const Char_t *title):
   fUseMixedVertex(kTRUE),
   fUseSPDVertex(kFALSE),
   fUseCorrelationVertex(kFALSE),
+  fSPDVtxResolution(kFALSE),
   fIsIPSigmacut(kFALSE),
   fIsIPcharge(kFALSE),
   fIsIPOpp(kFALSE),
@@ -233,7 +235,8 @@ AliHFEcuts::AliHFEcuts(const AliHFEcuts &c):
   fITSpatternCut(c.fITSpatternCut),
   fUseMixedVertex(kTRUE),
   fUseSPDVertex(kFALSE),
-  fUseCorrelationVertex(kFALSE),
+  fUseCorrelationVertex(c.fUseCorrelationVertex),
+  fSPDVtxResolution(c.fSPDVtxResolution),
   fIsIPSigmacut(kFALSE),
   fIsIPcharge(kFALSE),
   fIsIPOpp(kFALSE),
@@ -291,6 +294,7 @@ void AliHFEcuts::Copy(TObject &c) const {
   target.fUseMixedVertex = fUseMixedVertex;
   target.fUseSPDVertex = fUseSPDVertex;
   target.fUseCorrelationVertex = fUseCorrelationVertex;
+  target.fSPDVtxResolution = fSPDVtxResolution;
   target.fIsIPSigmacut = fIsIPSigmacut;
   target.fIsIPcharge = fIsIPcharge;
   target.fIsIPOpp = fIsIPOpp;
@@ -497,6 +501,7 @@ void AliHFEcuts::SetEventCutList(Int_t istep){
     if(fUseSPDVertex) evRecCuts->SetUseSPDVertex();
     if(fUseMixedVertex) evRecCuts->SetUseMixedVertex();
     if(fUseCorrelationVertex) evRecCuts->SetCheckCorrelationSPDVtx();
+    if(fSPDVtxResolution) evRecCuts->SetCheckSPDResolution();
     evRecCuts->SetVertexZCut(-fVertexRangeZ, fVertexRangeZ);
     //evRecCuts->SetVertexNContributors(1,(Int_t)1.e9);
     if(IsQAOn()) evRecCuts->SetQAOn(fHistQA);
index d3bf672..c49e4f7 100644 (file)
@@ -168,6 +168,7 @@ class AliHFEcuts : public TNamed{
     inline void SetUseMixedVertex(Bool_t useMixedVertex);    
     inline void SetUseSPDVertex(Bool_t useSPDVertex);
     void SetUseCorrelationVertex() { fUseCorrelationVertex = kTRUE;};
+    void SetSPDVtxResolutionCut() {fSPDVtxResolution = kTRUE;}
     void SetFractionOfSharedTPCClusters(Double_t fractionOfSharedTPCClusters) {fFractionOfSharedTPCClusters = fractionOfSharedTPCClusters;};
     void SetMaxImpactParameterRpar(Bool_t maxImpactParameterRpar) { fMaxImpactParameterRpar = maxImpactParameterRpar; };
     
@@ -249,6 +250,7 @@ class AliHFEcuts : public TNamed{
     Bool_t   fUseMixedVertex;         // Use primary vertex from track if there otherwise SPD vertex
     Bool_t   fUseSPDVertex;           // Use primary SPD vertex 
     Bool_t   fUseCorrelationVertex;   // Use the correlation of the vertex in z
+    Bool_t   fSPDVtxResolution;       // Check resolution of the SPD vertex
     Float_t  fIPCutParams[4];         // Parameters of impact parameter cut parametrization
     Bool_t   fIsIPSigmacut;           // if IP cut or IP sigma cut 
     Bool_t   fIsIPcharge;             // if cut on IP * charge (cut using only positive side of distribution, to eliminate conversions)
@@ -265,7 +267,7 @@ class AliHFEcuts : public TNamed{
 
     Int_t fDebugLevel;            // Debug Level
     
-  ClassDef(AliHFEcuts, 4)         // Container for HFE cuts
+  ClassDef(AliHFEcuts, 5)         // Container for HFE cuts
 };
 
 //__________________________________________________________________
index d8fea67..ba4b000 100644 (file)
@@ -42,6 +42,7 @@ AliHFEextraEventCuts::AliHFEextraEventCuts() :
   fVtxMixed(0),
   fVtxSPD(0),
   fCheckCorrelationSPDVtx(0),
+  fVtxResolution(kFALSE),
   fBitMap(0x0)
 {
   //
@@ -59,6 +60,7 @@ AliHFEextraEventCuts::AliHFEextraEventCuts(Char_t* name, Char_t* title) :
   fVtxMixed(0),
   fVtxSPD(0),
   fCheckCorrelationSPDVtx(0),
+  fVtxResolution(kFALSE),
   fBitMap(0x0)
  {
   //
@@ -78,6 +80,7 @@ AliHFEextraEventCuts::AliHFEextraEventCuts(const AliHFEextraEventCuts& c) :
   fVtxMixed(c.fVtxMixed),
   fVtxSPD(c.fVtxSPD),
   fCheckCorrelationSPDVtx(c.fCheckCorrelationSPDVtx),
+  fVtxResolution(c.fVtxResolution),
   fBitMap(c.fBitMap)
 {
   //
@@ -140,6 +143,7 @@ AliHFEextraEventCuts& AliHFEextraEventCuts::operator=(const AliHFEextraEventCuts
     fVtxMixed=c.fVtxMixed;
     fVtxSPD=c.fVtxSPD;
     fCheckCorrelationSPDVtx=c.fCheckCorrelationSPDVtx;
+    fVtxResolution = c.fVtxResolution;
     fBitMap=c.fBitMap;
   }
 
@@ -181,147 +185,76 @@ void AliHFEextraEventCuts::SelectionBitMap(TObject* obj) {
   //
 
   //Check if the requested cuts are passed and return a bitmap
-  for(Int_t j=0;j<kNCuts;j++)fBitMap->SetBitNumber(j,kFALSE);
-
-
-
-  AliESDEvent* esd = dynamic_cast<AliESDEvent *>(obj);
-  AliAODEvent* aod = dynamic_cast<AliAODEvent *>(obj);
-
-
-  // ESD
-
-  if ( esd ) {
-    
-    //now start checking the cuts,
-    //first assume the event will be accepted: 
-    for(Int_t j=0;j<kNCuts;j++)fBitMap->SetBitNumber(j,kTRUE);
-    
-    
-    
-    if(fRequireVtxCuts){
-      const AliESDVertex* vtxESD = 0x0;
-      if      (fVtxMixed) {
-       vtxESD = esd->GetPrimaryVertexTracks() ;
-       if((!vtxESD) || (vtxESD->GetNContributors() <= 0)) {
-         vtxESD = esd->GetPrimaryVertexSPD() ;
-         if((!vtxESD) || (vtxESD->GetNContributors() <= 0)) {
-           for(Int_t j=1;j<kNCuts;j++)fBitMap->SetBitNumber(j,kFALSE); 
-           AliWarning("Cannot get vertex, skipping event");
-           return;
-         }
-       }
-      }
-      else if(fVtxSPD) {   
-       vtxESD = esd->GetPrimaryVertexSPD() ;
-       if(fCheckCorrelationSPDVtx) {
-         const AliESDVertex* vtxESDtr = esd->GetPrimaryVertexTracks();
-         if((vtxESD) && (vtxESD->GetNContributors() > 0) && (vtxESDtr) && (vtxESDtr->GetNContributors() > 0)) {
-           if(TMath::Abs(vtxESD->GetZv()-vtxESDtr->GetZv())>0.5) return;
-         }
-         else return;
-       }
-      }
-      else {
-       vtxESD = esd->GetPrimaryVertexTracks() ;
-      }
-      
-      if(!vtxESD){
-       for(Int_t j=1;j<kNCuts;j++)fBitMap->SetBitNumber(j,kFALSE); 
-       AliWarning("Cannot get vertex, skipping event");
-       return;
-      }
-      
-      // Pick up the position and uncertainties
-      Double_t vtxPos[3];
-      vtxPos[0] = vtxESD->GetXv();
-      vtxPos[1] = vtxESD->GetYv();
-      vtxPos[2] = vtxESD->GetZv();
-      
-      Int_t nCtrb = vtxESD->GetNContributors();
-      
-      // Apply the cut
-      
-      if (vtxPos[2]>fVtxZMax || vtxPos[2]<fVtxZMin)
-       fBitMap->SetBitNumber(0,kFALSE); 
-      if (nCtrb<fVtxNCtrbMin)
-       fBitMap->SetBitNumber(1,kFALSE);
-      
-    }  
-    return;
-    
-    
-  }
 
-  //
-
-  // AOD
-
-  if ( aod ) {
-    
-    //now start checking the cuts,
-    //first assume the event will be accepted: 
-    for(Int_t j=0;j<kNCuts;j++)fBitMap->SetBitNumber(j,kTRUE);
-    
-    
-    
-    if(fRequireVtxCuts){
-      const AliAODVertex* vtxAOD = 0x0;
-      if      (fVtxMixed) {
-       vtxAOD = aod->GetPrimaryVertex();
-       if((!vtxAOD) || (vtxAOD->GetNContributors() <= 0)) {
-         vtxAOD = aod->GetPrimaryVertexSPD() ;
-         if((!vtxAOD) || (vtxAOD->GetNContributors() <= 0)) {
-           for(Int_t j=1;j<kNCuts;j++)fBitMap->SetBitNumber(j,kFALSE); 
-           AliWarning("Cannot get vertex, skipping event");
-           return;
-         }
-       }
-      }
-      else if(fVtxSPD) {   
-       vtxAOD = aod->GetPrimaryVertexSPD() ;
-       if(fCheckCorrelationSPDVtx) {
-         const AliAODVertex* vtxAODtr = aod->GetPrimaryVertex();
-         if((vtxAOD) && (vtxAOD->GetNContributors() > 0) && (vtxAODtr) && (vtxAODtr->GetNContributors() > 0)) {
-           if(TMath::Abs(vtxAOD->GetZ()-vtxAODtr->GetZ())>0.5) return;
-         }
-         else return;
-       }
-      }
-      else {
-       vtxAOD = aod->GetPrimaryVertex() ;
-      }
-      
-      if(!vtxAOD){
-       for(Int_t j=1;j<kNCuts;j++)fBitMap->SetBitNumber(j,kFALSE); 
-       AliWarning("Cannot get vertex, skipping event");
-       return;
-      }
-      
-      // Pick up the position and uncertainties
-      Double_t vtxPos[3];
-      vtxPos[0] = vtxAOD->GetX();
-      vtxPos[1] = vtxAOD->GetY();
-      vtxPos[2] = vtxAOD->GetZ();
-      
-      Int_t nCtrb = vtxAOD->GetNContributors();
-      //printf("AliHFEextraCuts:: %d, %f, %f, %f\n",nCtrb,vtxPos[0],vtxPos[1],vtxPos[2]);
-      
-      // Apply the cut
-      
-      if (vtxPos[2]>fVtxZMax || vtxPos[2]<fVtxZMin)
-       fBitMap->SetBitNumber(0,kFALSE); 
-      if (nCtrb<fVtxNCtrbMin)
-       fBitMap->SetBitNumber(1,kFALSE);
-      
-    }  
+  AliVEvent *inputEvent = dynamic_cast<AliVEvent *>(obj);
+  if(!inputEvent){
+    AliDebug(1, "Not a virtual event");
+    for(Int_t j=0;j<kNCuts;j++)fBitMap->SetBitNumber(j,kFALSE);
     return;
-    
-    
   }
+  const AliVVertex *vtxTracks = GetPrimaryVertexTracks(inputEvent),
+                   *vtxSPD = GetPrimaryVertexSPD(inputEvent),
+                   *vtxPrim(NULL);
   
-  return;
-  
+  //first assume the event will be accepted: 
+  for(Int_t j=0;j<kNCuts;j++)fBitMap->SetBitNumber(j,kTRUE);
+
+  if(fVtxMixed){
+    // Use mixed vertex: Prefer vertex with tracks, in case not available use SPD vertex
+    if(vtxTracks && vtxTracks->GetNContributors() > 0) vtxPrim = vtxTracks;
+    else if(vtxSPD && vtxSPD->GetNContributors() > 0) vtxPrim = vtxSPD;
+  } else if(fVtxSPD){
+    if(vtxSPD && vtxSPD->GetNContributors () > 0) vtxPrim = vtxSPD;
+  } else {
+    if(vtxTracks && vtxTracks->GetNContributors() > 0) vtxPrim = vtxTracks;
+  }
+  if(!vtxPrim){
+    // No primary vertex: Reject event
+         for(Int_t j=1;j<kNCuts;j++)fBitMap->SetBitNumber(j,kFALSE); 
+         AliWarning("Cannot get vertex, skipping event");
+         return;
+  }
+  // Check standard vertex cuts using the primary vertex 
+  if(fVtxZMin > - 1000. && fVtxZMax < 1000.){
+    // Primary vertex z cut required
+    // Pick up the position and uncertainties
+    Double_t vtxPos[3];
+    vtxPos[0] = vtxPrim->GetX();
+    vtxPos[1] = vtxPrim->GetY();
+    vtxPos[2] = vtxPrim->GetZ();
+    if (vtxPos[2]>fVtxZMax || vtxPos[2]<fVtxZMin)
+           fBitMap->SetBitNumber(kVtxPosZ,kFALSE); 
+  }
+  if(fVtxNCtrbMin > 0){  
+    // cut required if the cut value is set to something larger than 0
+    // same effect as setting the min. cut value to 0
+    Int_t nCtrb = vtxPrim->GetNContributors();
+    if (nCtrb<fVtxNCtrbMin)
+           fBitMap->SetBitNumber(kVtxNCtrb,kFALSE);
+  }
+
+  // check vertex correlation cut
+  if(fCheckCorrelationSPDVtx){
+    if(vtxTracks && vtxTracks->GetNContributors() && vtxSPD && vtxSPD->GetNContributors()){
+      // Both vertices available, check correlation
+      if(TMath::Abs(vtxTracks->GetZ() - vtxSPD->GetZ()) >= 0.5) fBitMap->SetBitNumber(kCorrelation, kFALSE);
+    } else {
+      // No correlation available: set the cut to false
+      fBitMap->SetBitNumber(kCorrelation, kFALSE); 
+    }
+  }
+  // check cut on the SPD vertex resolution
+  if(fVtxResolution){
+    if(vtxSPD){
+      TString vtxTyp = vtxSPD->GetTitle();
+      Double_t cov[6]={0};
+      vtxSPD->GetCovarianceMatrix(cov);
+      Double_t zRes = TMath::Sqrt(cov[5]);
+      if (vtxTyp.Contains("vertexer:Z") && (zRes>0.25)) fBitMap->SetBitNumber(kResolution, kFALSE);
+    } else{
+      fBitMap->SetBitNumber(kResolution, kFALSE);
+    }
+  }
 }
 
 //_____________________________________________________________________________
@@ -332,35 +265,41 @@ void AliHFEextraEventCuts::FillHistograms(TObject* obj, Bool_t b)
   //
 
   if(!fIsQAOn) return;
-  // cast TObject into VParticle
-  AliESDEvent* esd = dynamic_cast<AliESDEvent *>(obj);
-  if (!esd ) return  ;
+  // cast TObject into VEvent
+  AliVEvent* inputEvent = dynamic_cast<AliESDEvent *>(obj);
+  if (!inputEvent) return;
 
   // index = 0: fill histograms before cuts
   // index = 1: fill histograms after cuts
   Int_t index = -1;
   index = ((b) ? 1 : 0);
 
+  // Obtain vertices
+  const AliVVertex *vtxTracks = GetPrimaryVertexTracks(inputEvent),
+                   *vtxSPD = GetPrimaryVertexSPD(inputEvent),
+                   *vtxPrim(NULL);
 
   //look at vertex parameters:
-  const AliESDVertex* vtxESD = 0x0;
-  if      (fVtxMixed) {
-    vtxESD = esd->GetPrimaryVertexTracks() ;
-    if((!vtxESD) || (vtxESD->GetNContributors() <= 0)) {
-      vtxESD = esd->GetPrimaryVertexSPD() ;
-      if((!vtxESD) || (vtxESD->GetNContributors() <= 0)) {
-       return;
-      }
-    }
+  if(fVtxMixed){
+    if(vtxTracks->GetNContributors() > 0) vtxPrim = vtxTracks;
+    else if(vtxSPD->GetNContributors()) vtxPrim = vtxSPD;
   }
   else {   
-    vtxESD = esd->GetPrimaryVertexTracks() ;
+    vtxPrim = vtxTracks; 
   }
-  if(!vtxESD)return;
+  if(!vtxPrim)return;
   // vertex position and uncertainties
-  fhQA[kVtxPosZ] [index]->Fill(vtxESD->GetZv());
-  fhQA[kVtxNCtrb][index]->Fill(vtxESD->GetNContributors());
-  
+  fhQA[kVtxPosZ] [index]->Fill(vtxPrim->GetZ());
+  fhQA[kVtxNCtrb][index]->Fill(vtxPrim->GetNContributors());
+  // SPD Vertex Position Correlation 
+  if(vtxTracks && vtxSPD){
+    fhQA[kCorrelation][index]->Fill(vtxTracks->GetZ()-vtxSPD->GetZ());
+  }
+  if(vtxSPD){
+    Double_t cov[6]={0};
+    vtxSPD->GetCovarianceMatrix(cov);
+    fhQA[kResolution][index]->Fill(TMath::Sqrt(cov[5]));
+  }
 }
 
 //____________________________________________________________________
@@ -405,7 +344,6 @@ void AliHFEextraEventCuts::SetHistogramBins(Int_t index, Int_t nbins, Double_t x
   }  
   
   // book QA histograms
-
   Char_t str[5];
   for (Int_t i=0; i<kNStepQA; i++) {
     if (i==0) snprintf(str,5," ");
@@ -413,6 +351,8 @@ void AliHFEextraEventCuts::SetHistogramBins(Int_t index, Int_t nbins, Double_t x
 
     fhQA[kVtxPosZ][i]  = new  TH1F(Form("%s_Vtx_Pos_Z%s",GetName(),str),               "",200,-50.,50.);
     fhQA[kVtxNCtrb][i] = new  TH1F(Form("%s_Vtx_N_Ctrb%s",GetName(),str),              "",1000,0.,1000.);
+    fhQA[kCorrelation][i] = new TH1F(Form("%s_SPDCorrelation_%s",GetName(),str), "",200, -10., 10);
+    fhQA[kResolution][i] = new TH1F(Form("%s_SPDResolution_%s",GetName(), str), "", 100, 0., 1.); 
  
     fhQA[kVtxPosZ][i]  ->SetXTitle("Vertex Position Z (cm)");
     fhQA[kVtxNCtrb][i] ->SetXTitle("Number of contributors");
@@ -435,3 +375,41 @@ void AliHFEextraEventCuts::AddQAHistograms(TList *qaList) {
        qaList->Add(fhQA[i][j]);
   }
 }
+
+//_____________________________________________________________________________
+const AliVVertex *AliHFEextraEventCuts::GetPrimaryVertexSPD(const AliVEvent * const inputEvent){
+  //
+  // Get SPD vertex from event
+  //
+  const AliVVertex *spdvtx(NULL);
+  const AliESDEvent *esd(NULL);
+  const AliAODEvent *aod(NULL);
+  if((esd = dynamic_cast<const AliESDEvent *>(inputEvent))){
+    spdvtx = esd->GetPrimaryVertexSPD();
+  } else if((aod = dynamic_cast<const AliAODEvent *>(inputEvent))){
+    spdvtx = aod->GetPrimaryVertexSPD();
+  }
+  return spdvtx;
+}
+
+//_____________________________________________________________________________
+const AliVVertex *AliHFEextraEventCuts::GetPrimaryVertexTracks(const AliVEvent *const inputEvent){
+  // 
+  // Get Primary Vertex from tracks
+  //
+  const AliVVertex *trkvtx(NULL);
+  const AliESDEvent *esd(NULL);
+  const AliAODEvent *aod(NULL);
+  if((esd = dynamic_cast<const AliESDEvent *>(inputEvent))){
+    trkvtx = esd->GetPrimaryVertexTracks();
+  } else if((aod = dynamic_cast<const AliAODEvent *>(inputEvent))){
+    const AliVVertex *vtxTmp = aod->GetPrimaryVertex();
+    // check whether the primary vertex is the vertex from tracks
+    TString vtxTtl = vtxTmp->GetTitle();
+    if(vtxTtl.Contains("VertexerTracks")){
+      trkvtx = vtxTmp;
+    }
+  }
+  return trkvtx;
+}
+
index 1fceb42..660b57a 100644 (file)
@@ -27,6 +27,8 @@
 #include "AliCFCutBase.h"
 class TH1F;
 class TBits;
+class AliVVertex;
+class AliVEvent;
 //_____________________________________________________________________________
 class AliHFEextraEventCuts: public AliCFCutBase 
 {
@@ -46,6 +48,7 @@ class AliHFEextraEventCuts: public AliCFCutBase
   void SetUseMixedVertex() {fVtxMixed=kTRUE; fVtxSPD=kFALSE;} //default is vertex from tracks
   void SetUseSPDVertex() {fVtxSPD=kTRUE; fVtxMixed=kFALSE;} //default is vertex from tracks
   void SetCheckCorrelationSPDVtx() {fCheckCorrelationSPDVtx=kTRUE;} // check the correlation between z of different vertices
+  void SetCheckSPDResolution() {fVtxResolution = kTRUE;} // check resolution on the SPD vertex
  
   Bool_t   GetRequireVtxCuts() const {return fRequireVtxCuts;} // cut value getter
   Double_t GetVertexZMax() const {return fVtxZMax;} // cut values getter
@@ -57,6 +60,8 @@ class AliHFEextraEventCuts: public AliCFCutBase
   void SetHistogramBins(Int_t index, Int_t nbins, Double_t xmin, Double_t xmax);
   enum{kVtxPosZ,
        kVtxNCtrb,
+       kCorrelation,
+       kResolution,
        kNCuts,
        kNStepQA=2
   };
@@ -67,6 +72,8 @@ class AliHFEextraEventCuts: public AliCFCutBase
   void DefineHistograms();             // books histograms 
   void Initialise();                   // sets everything to 0
   void FillHistograms(TObject* obj, Bool_t b);
+  const AliVVertex *GetPrimaryVertexSPD(const AliVEvent * const inputEvent);
+  const AliVVertex *GetPrimaryVertexTracks(const AliVEvent *const inputEvent);
 
   Bool_t fRequireVtxCuts ; //The type of trigger to be checked
   Double_t fVtxZMax ; //Z vertex position, maximum value
@@ -75,6 +82,7 @@ class AliHFEextraEventCuts: public AliCFCutBase
   Bool_t   fVtxMixed;  //Flag for use of mixed vertex (primary vertex with track, if not SPD vertex)
   Bool_t   fVtxSPD;    //Flag for use of SPD vertex 
   Bool_t   fCheckCorrelationSPDVtx;   //Check the correlation SPD, track vertex
+  Bool_t   fVtxResolution;            //Check vertex resolution cut
  
   TBits *fBitMap ; //cut mask