Pile-up vertices added to AODEvent. (F. Prino)
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 9 Jul 2010 10:57:29 +0000 (10:57 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 9 Jul 2010 10:57:29 +0000 (10:57 +0000)
ANALYSIS/AliAnalysisTaskESDfilter.cxx
STEER/AliAODEvent.cxx
STEER/AliAODEvent.h
STEER/AliESDEvent.cxx
STEER/AliESDEvent.h
STEER/AliVEvent.h

index b061b04..323316f 100644 (file)
@@ -187,12 +187,17 @@ void AliAnalysisTaskESDfilter::ConvertESDtoAOD() {
     Float_t diamxy[2]={esd->GetDiamondX(),esd->GetDiamondY()};\r
     Float_t diamcov[3]; esd->GetDiamondCovXY(diamcov);\r
     header->SetDiamond(diamxy,diamcov);\r
+    header->SetDiamondZ(esd->GetDiamondZ(),esd->GetSigma2DiamondZ());\r
 //\r
 //\r
     Int_t nV0s      = esd->GetNumberOfV0s();\r
     Int_t nCascades = esd->GetNumberOfCascades();\r
     Int_t nKinks    = esd->GetNumberOfKinks();\r
     Int_t nVertices = nV0s + nCascades /*V0 wihtin cascade already counted*/+ nKinks + 1 /* = prim. vtx*/;\r
+    Int_t nPileSPDVertices=1+esd->GetNumberOfPileupVerticesSPD(); // also SPD main vertex\r
+    Int_t nPileTrkVertices=esd->GetNumberOfPileupVerticesTracks();\r
+    nVertices+=nPileSPDVertices;\r
+    nVertices+=nPileTrkVertices;\r
     Int_t nJets     = 0;\r
     Int_t nCaloClus = esd->GetNumberOfCaloClusters();\r
     Int_t nFmdClus  = 0;\r
@@ -280,6 +285,41 @@ void AliAnalysisTaskESDfilter::ConvertESDtoAOD() {
 \r
     if (fDebug > 0) primary->Print();\r
 \r
+\r
+    // Add SPD "main" vertex \r
+    const AliESDVertex *vtxS = esd->GetPrimaryVertexSPD();\r
+    vtxS->GetXYZ(pos); // position\r
+    vtxS->GetCovMatrix(covVtx); //covariance matrix\r
+    AliAODVertex * mVSPD = new(vertices[jVertices++])\r
+      AliAODVertex(pos, covVtx, vtxS->GetChi2toNDF(), NULL, -1, AliAODVertex::kMainSPD);\r
+    mVSPD->SetName(vtxS->GetName());\r
+    mVSPD->SetTitle(vtxS->GetTitle());\r
+    mVSPD->SetNContributors(vtxS->GetNContributors()); \r
+\r
+    // Add SPD pileup vertices\r
+    for(Int_t iV=0; iV<nPileSPDVertices-1; iV++){\r
+      const AliESDVertex *vtxP = esd->GetPileupVertexSPD(iV);\r
+      vtxP->GetXYZ(pos); // position\r
+      vtxP->GetCovMatrix(covVtx); //covariance matrix\r
+      AliAODVertex * pVSPD =  new(vertices[jVertices++])\r
+       AliAODVertex(pos, covVtx, vtxP->GetChi2toNDF(), NULL, -1, AliAODVertex::kPileupSPD);\r
+      pVSPD->SetName(vtxP->GetName());\r
+      pVSPD->SetTitle(vtxP->GetTitle());\r
+      pVSPD->SetNContributors(vtxP->GetNContributors()); \r
+    }\r
+\r
+    // Add TRK pileup vertices\r
+    for(Int_t iV=0; iV<nPileTrkVertices; iV++){\r
+      const AliESDVertex *vtxP = esd->GetPileupVertexTracks(iV);\r
+      vtxP->GetXYZ(pos); // position\r
+      vtxP->GetCovMatrix(covVtx); //covariance matrix\r
+      AliAODVertex * pVTRK = new(vertices[jVertices++])\r
+       AliAODVertex(pos, covVtx, vtxP->GetChi2toNDF(), NULL, -1, AliAODVertex::kPileupTracks);\r
+      pVTRK->SetName(vtxP->GetName());\r
+      pVTRK->SetTitle(vtxP->GetTitle());\r
+      pVTRK->SetNContributors(vtxP->GetNContributors());\r
+    }\r
+\r
     // Create vertices starting from the most complex objects\r
     Double_t chi2 = 0.;\r
     \r
index 4c4110c..4dbb5dd 100644 (file)
@@ -373,7 +373,7 @@ void AliAODEvent::ResetStd(Int_t trkArrSize,
   fVertices->Delete();
   if (vtxArrSize > fVertices->GetSize()) 
     fVertices->Expand(vtxArrSize);
-        
+
   fV0s->Delete();
   if (v0ArrSize > fV0s->GetSize()) 
     fV0s->Expand(v0ArrSize);
@@ -618,6 +618,118 @@ void AliAODEvent::ReadFromTree(TTree *tree, Option_t* opt /*= ""*/)
     fAODObjects->SetOwner(kTRUE);
   }
 }
+//______________________________________________________________________________
+Int_t  AliAODEvent::GetNumberOfPileupVerticesSPD() const{
+  // count number of SPD pileup vertices
+  Int_t nVertices=GetNumberOfVertices();
+  Int_t nPileupVertices=0;
+  for(Int_t iVert=0; iVert<nVertices; iVert++){
+    AliAODVertex *v=GetVertex(iVert);
+    if(v->GetType()==AliAODVertex::kPileupSPD) nPileupVertices++;
+  }
+  return nPileupVertices;
+}
+//______________________________________________________________________________
+Int_t  AliAODEvent::GetNumberOfPileupVerticesTracks() const{
+  // count number of track pileup vertices
+  Int_t nVertices=GetNumberOfVertices();
+  Int_t nPileupVertices=0;
+  for(Int_t iVert=0; iVert<nVertices; iVert++){
+    AliAODVertex *v=GetVertex(iVert);
+    if(v->GetType()==AliAODVertex::kPileupTracks) nPileupVertices++;
+  }
+  return nPileupVertices;
+}
+//______________________________________________________________________________
+AliAODVertex* AliAODEvent::GetPrimaryVertexSPD() const{
+  //
+  Int_t nVertices=GetNumberOfVertices();
+  for(Int_t iVert=0; iVert<nVertices; iVert++){
+    AliAODVertex *v=GetVertex(iVert);
+    if(v->GetType()==AliAODVertex::kMainSPD) return v;
+  }
+  return 0;
+}
+//______________________________________________________________________________
+AliAODVertex* AliAODEvent::GetPileupVertexSPD(Int_t iV) const{
+  //
+  Int_t nVertices=GetNumberOfVertices();
+  Int_t counter=0;
+  for(Int_t iVert=0; iVert<nVertices; iVert++){
+    AliAODVertex *v=GetVertex(iVert);
+    if(v->GetType()==AliAODVertex::kPileupSPD){
+      if(counter==iV) return v;
+      ++counter;
+    }
+  }
+  return 0;
+}
+//______________________________________________________________________________
+AliAODVertex* AliAODEvent::GetPileupVertexTracks(Int_t iV) const{
+  //
+  Int_t nVertices=GetNumberOfVertices();
+  Int_t counter=0;
+  for(Int_t iVert=0; iVert<nVertices; iVert++){
+    AliAODVertex *v=GetVertex(iVert);
+    if(v->GetType()==AliAODVertex::kPileupTracks){
+      if(counter==iV) return v;
+      ++counter;
+    }
+  }
+  return 0;
+}
+//______________________________________________________________________________
+Bool_t  AliAODEvent::IsPileupFromSPD(Int_t minContributors, 
+                                    Double_t minZdist, 
+                                    Double_t nSigmaZdist, 
+                                    Double_t nSigmaDiamXY, 
+                                    Double_t nSigmaDiamZ) const{
+  //
+  // This function checks if there was a pile up
+  // reconstructed with SPD
+  //
+  AliAODVertex *mainV=GetPrimaryVertexSPD();
+  if(!mainV) return kFALSE;
+  Int_t nc1=mainV->GetNContributors();
+  if(nc1<1) return kFALSE;
+  Int_t nPileVert=GetNumberOfPileupVerticesSPD();
+  if(nPileVert==0) return kFALSE;
+  Int_t nVertices=GetNumberOfVertices();
+  
+  for(Int_t iVert=0; iVert<nVertices; iVert++){
+    AliAODVertex *pv=GetVertex(iVert);
+    if(pv->GetType()!=AliAODVertex::kPileupSPD) continue;
+    Int_t nc2=pv->GetNContributors();
+    if(nc2>=minContributors){
+      Double_t z1=mainV->GetZ();
+      Double_t z2=pv->GetZ();
+      Double_t distZ=TMath::Abs(z2-z1);
+      Double_t distZdiam=TMath::Abs(z2-GetDiamondZ());
+      Double_t cutZdiam=nSigmaDiamZ*GetSigma2DiamondZ();
+      if(GetSigma2DiamondZ()<0.0001)cutZdiam=99999.; //protection for missing z diamond information
+      if(distZ>minZdist && distZdiam<cutZdiam){
+       Double_t x2=pv->GetX();
+       Double_t y2=pv->GetY();
+       Double_t distXdiam=TMath::Abs(x2-GetDiamondX());
+       Double_t distYdiam=TMath::Abs(y2-GetDiamondY());
+       Double_t cov1[6],cov2[6];       
+       mainV->GetCovarianceMatrix(cov1);
+       pv->GetCovarianceMatrix(cov2);
+       Double_t errxDist=TMath::Sqrt(cov2[0]+GetSigma2DiamondX());
+       Double_t erryDist=TMath::Sqrt(cov2[2]+GetSigma2DiamondY());
+       Double_t errzDist=TMath::Sqrt(cov1[5]+cov2[5]);
+       Double_t cutXdiam=nSigmaDiamXY*errxDist;
+       if(GetSigma2DiamondX()<0.0001)cutXdiam=99999.; //protection for missing diamond information
+       Double_t cutYdiam=nSigmaDiamXY*erryDist;
+       if(GetSigma2DiamondY()<0.0001)cutYdiam=99999.; //protection for missing diamond information
+       if( (distXdiam<cutXdiam) && (distYdiam<cutYdiam) && (distZ>nSigmaZdist*errzDist) ){
+         return kTRUE;
+       }
+      }
+    }
+  }
+  return kFALSE;
+}
 
 //______________________________________________________________________________
 void AliAODEvent::Print(Option_t *) const
index 3b99fb1..f54cdd3 100644 (file)
@@ -80,7 +80,7 @@ class AliAODEvent : public AliVEvent {
   void     SetMagneticField(Double_t mf){if (fHeader) fHeader->SetMagneticField(mf);}
   void     SetMuonMagFieldScale(Double_t mf){if (fHeader) fHeader->SetMuonMagFieldScale(mf);}
   void     SetDiamond(Float_t xy[2],Float_t cov[3]){if (fHeader) fHeader->SetDiamond(xy,cov);}
-
+  void     SetDiamondZ(Float_t z, Float_t sig2z){if (fHeader) fHeader->SetDiamondZ(z,sig2z);}
   Int_t    GetRunNumber() const {return fHeader ? fHeader->GetRunNumber() : -999;}
   UInt_t   GetPeriodNumber() const {return fHeader ? fHeader->GetPeriodNumber() : 0;}
   UInt_t   GetOrbitNumber() const {return fHeader ? fHeader->GetOrbitNumber() : 0;}
@@ -89,7 +89,11 @@ class AliAODEvent : public AliVEvent {
   Double_t GetMuonMagFieldScale() const {return fHeader ? fHeader->GetMuonMagFieldScale() : -999.;}
   Double_t GetDiamondX() const {return fHeader ? fHeader->GetDiamondX() : -999.;}
   Double_t GetDiamondY() const {return fHeader ? fHeader->GetDiamondY() : -999.;}
+  Double_t GetDiamondZ() const {return fHeader ? fHeader->GetDiamondZ() : -999.;}
   void     GetDiamondCovXY(Float_t cov[3]) const {cov[0]=-999.; if(fHeader) fHeader->GetDiamondCovXY(cov);}
+  Double_t GetSigma2DiamondX() const {return fHeader ? fHeader->GetSigma2DiamondX() : -999.;}
+  Double_t GetSigma2DiamondY() const {return fHeader ? fHeader->GetSigma2DiamondY() : -999.;}
+  Double_t GetSigma2DiamondZ() const {return fHeader ? fHeader->GetSigma2DiamondZ() : -999.;}
   
   void      SetEventType(UInt_t eventType){fHeader->SetEventType(eventType);}
   void      SetTriggerMask(ULong64_t n) {fHeader->SetTriggerMask(n);}
@@ -124,6 +128,16 @@ class AliAODEvent : public AliVEvent {
   
   // primary vertex
   virtual AliAODVertex *GetPrimaryVertex() const { return GetVertex(0); }
+  virtual AliAODVertex *GetPrimaryVertexSPD() const;
+
+  // -- Pileup vertices 
+  Int_t         GetNumberOfPileupVerticesTracks()   const;
+  Int_t         GetNumberOfPileupVerticesSPD()    const;
+  virtual AliAODVertex *GetPileupVertexSPD(Int_t iV=0) const;
+  virtual AliAODVertex *GetPileupVertexTracks(Int_t iV=0) const;
+  virtual Bool_t  IsPileupFromSPD(Int_t minContributors=3, Double_t minZdist=0.8, Double_t nSigmaZdist=3., Double_t nSigmaDiamXY=2., Double_t nSigmaDiamZ=5.) const;
+
+
 
   // V0
   TClonesArray *GetV0s()                 const { return fV0s; }
index 4f8f55e..5ef8893 100644 (file)
@@ -1603,64 +1603,54 @@ Bool_t    AliESDEvent::IsHLTTriggerFired(const char* name) const
   return kTRUE;
 }
 
-Bool_t  AliESDEvent::IsPileupFromSPD(Int_t ncont, Double_t distz, Double_t nSigmaDeltaZ, Double_t nSigmaXY, Int_t option) const {
+//______________________________________________________________________________
+Bool_t  AliESDEvent::IsPileupFromSPD(Int_t minContributors, 
+                                    Double_t minZdist, 
+                                    Double_t nSigmaZdist, 
+                                    Double_t nSigmaDiamXY, 
+                                    Double_t nSigmaDiamZ) const{
   //
   // This function checks if there was a pile up
   // reconstructed with SPD
   //
-  Double_t diamx= GetDiamondX();
-  Double_t diamsigma2x= GetSigma2DiamondX();
-  Double_t diamy= GetDiamondY();
-  Double_t diamsigma2y= GetSigma2DiamondY();
-  
-  Double_t sigmax= TMath::Sqrt(diamsigma2x);
-  Double_t sigmay= TMath::Sqrt(diamsigma2y);
-  
-  Double_t z1=fSPDVertex->GetZ();
   Int_t nc1=fSPDVertex->GetNContributors();
   if(nc1<1) return kFALSE;
   Int_t nPileVert=GetNumberOfPileupVerticesSPD();
   if(nPileVert==0) return kFALSE;
+  
   for(Int_t i=0; i<nPileVert;i++){
     const AliESDVertex* pv=GetPileupVertexSPD(i);
-    Double_t z2=pv->GetZ();
-    Double_t x2=pv->GetX();
-    Double_t y2=pv->GetY();
     Int_t nc2=pv->GetNContributors();
-    Double_t distanceZ=TMath::Abs(z2-z1);
-    Double_t distanceX=TMath::Abs(x2-diamx);
-    Double_t distanceY=TMath::Abs(y2-diamy);
-    Double_t errzDist=0.;
-    Double_t errxDist=0.;
-    Double_t erryDist=0.;
-    if(option==0){
-      Double_t ez1=fSPDVertex->GetZRes();
-      Double_t ez2=pv->GetZRes();
-      errzDist=TMath::Sqrt(ez1*ez1+ez2*ez2);
-    }else{  
-      Double_t resol1=-75.6+834.6/TMath::Sqrt((Double_t)nc1);
-      resol1/=10000.;
-      Double_t resol2=-75.6+834.6/TMath::Sqrt((Double_t)nc2);
-      resol2/=10000.;
-      errzDist=TMath::Sqrt(resol1*resol1+resol2*resol2); 
+    if(nc2>=minContributors){
+      Double_t z1=fSPDVertex->GetZ();
+      Double_t z2=pv->GetZ();
+      Double_t distZ=TMath::Abs(z2-z1);
+      Double_t distZdiam=TMath::Abs(z2-GetDiamondZ());
+      Double_t cutZdiam=nSigmaDiamZ*GetSigma2DiamondZ();
+      if(GetSigma2DiamondZ()<0.0001)cutZdiam=99999.; //protection for missing z diamond information
+      if(distZ>minZdist && distZdiam<cutZdiam){
+       Double_t x2=pv->GetX();
+       Double_t y2=pv->GetY();
+       Double_t distXdiam=TMath::Abs(x2-GetDiamondX());
+       Double_t distYdiam=TMath::Abs(y2-GetDiamondY());
+       Double_t cov1[6],cov2[6];       
+       fSPDVertex->GetCovarianceMatrix(cov1);
+       pv->GetCovarianceMatrix(cov2);
+       Double_t errxDist=TMath::Sqrt(cov2[0]+GetSigma2DiamondX());
+       Double_t erryDist=TMath::Sqrt(cov2[2]+GetSigma2DiamondY());
+       Double_t errzDist=TMath::Sqrt(cov1[5]+cov2[5]);
+       Double_t cutXdiam=nSigmaDiamXY*errxDist;
+       if(GetSigma2DiamondX()<0.0001)cutXdiam=99999.; //protection for missing diamond information
+       Double_t cutYdiam=nSigmaDiamXY*erryDist;
+       if(GetSigma2DiamondY()<0.0001)cutYdiam=99999.; //protection for missing diamond information
+       if( (distXdiam<cutXdiam) && (distYdiam<cutYdiam) && (distZ>nSigmaZdist*errzDist) ){
+         return kTRUE;
+       }
+      }
     }
-
-    
-      Double_t ex2 = pv->GetXRes();
-      Double_t ey2= pv->GetYRes();
-      errxDist=TMath::Sqrt(ex2*ex2+sigmax*sigmax); 
-      erryDist=TMath::Sqrt(ey2*ey2+sigmay*sigmay); 
-    
-      if(nc2>=ncont && distanceZ>nSigmaDeltaZ*errzDist && distanceX<nSigmaXY*errxDist && distanceY<nSigmaXY*erryDist && distanceZ>distz)
-       
-       return kTRUE;
-  }
-           
-           
-      
-      return kFALSE;
   }
-
+  return kFALSE;
+}
 
 //______________________________________________________________________________
 void AliESDEvent::EstimateMultiplicity(Int_t &tracklets, Int_t &trITSTPC, Int_t &trITSSApure, Double_t eta, Bool_t useDCAFlag,Bool_t useV0Flag) const
index a212270..abda314 100644 (file)
@@ -254,13 +254,17 @@ public:
     return (const AliESDVertex *)(fSPDPileupVertices?fSPDPileupVertices->UncheckedAt(i):0x0);
   }
   Char_t  AddPileupVertexSPD(const AliESDVertex *vtx);
-  Bool_t  IsPileupFromSPD(Int_t ncont=2, Double_t distz=1., Double_t nSigmaDeltaZ=3., Double_t nSigmaXY=2., Int_t option=0) const;
-
   const AliESDVertex *GetPileupVertexTracks(Int_t i) const {
     return (const AliESDVertex *)(fTrkPileupVertices?fTrkPileupVertices->UncheckedAt(i):0x0);
   }
   Char_t  AddPileupVertexTracks(const AliESDVertex *vtx);
 
+  virtual Bool_t  IsPileupFromSPD(Int_t minContributors=3, 
+                                 Double_t minZdist=0.8, 
+                                 Double_t nSigmaZdist=3., 
+                                 Double_t nSigmaDiamXY=2., 
+                                 Double_t nSigmaDiamZ=5.) const;
+  
   AliESDtrack *GetTrack(Int_t i) const {
     return (AliESDtrack *)(fTracks?fTracks->UncheckedAt(i):0x0);
   }
index 0e79e70..a1707ba 100644 (file)
@@ -91,7 +91,14 @@ public:
 
   // Primary vertex
   virtual const AliVVertex   *GetPrimaryVertex() const {return 0x0;}
-
+  virtual Bool_t IsPileupFromSPD(Int_t /*minContributors*/, 
+                                Double_t /*minZdist*/, 
+                                Double_t /*nSigmaZdist*/, 
+                                Double_t /*nSigmaDiamXY*/, 
+                                Double_t /*nSigmaDiamZ*/)
+                                const{
+    return kFALSE;
+  }
   //---------- end of new stuff