]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
Cut class for the filtering of ESD V0s (Boris Hippolyte <hippolyt@in2p3.fr>)
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 15 May 2008 09:08:23 +0000 (09:08 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 15 May 2008 09:08:23 +0000 (09:08 +0000)
ANALYSIS/ANALYSISaliceLinkDef.h
ANALYSIS/AliESDv0Cuts.cxx [new file with mode: 0644]
ANALYSIS/AliESDv0Cuts.h [new file with mode: 0644]
ANALYSIS/libANALYSISalice.pkg

index a73ebfb8afc32fbb83afc852aee35e01acfa2eea..608a98ec14a5e558b0d60dd9caccb89f0af2659c 100644 (file)
@@ -9,6 +9,7 @@
 #pragma link C++ class  AliAnalysisTaskKineFilter+;
 #pragma link C++ class  AliKineTrackCuts+;
 #pragma link C++ class  AliESDtrackCuts+;
+#pragma link C++ class  AliESDv0Cuts+;
 #ifdef WITHXML
 #pragma link C++ class AliTagAnalysis+;
 #pragma link C++ class AliXMLCollection+;
diff --git a/ANALYSIS/AliESDv0Cuts.cxx b/ANALYSIS/AliESDv0Cuts.cxx
new file mode 100644 (file)
index 0000000..d316e9b
--- /dev/null
@@ -0,0 +1,842 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id: AliESDv0Cuts.cxx 24534 2008-03-16 22:22:11Z fca $ */
+
+#include "AliESDv0Cuts.h"
+
+#include <AliESDVertex.h>
+#include <AliESDtrack.h>
+#include <AliESDv0.h>
+#include <AliESD.h>
+#include <AliESDEvent.h>
+#include <AliLog.h>
+
+#include <TTree.h>
+#include <TCanvas.h>
+#include <TDirectory.h>
+
+//____________________________________________________________________
+ClassImp(AliESDv0Cuts)
+
+// Cut names
+const Char_t* AliESDv0Cuts::fgkCutNames[kNCuts] = {
+ "dca positive to pvtx",
+ "dca negative to pvtx",
+ "#Chi^{2}",
+ "dca v0 daughters",
+ "min decay radius",
+ "max decay radius",
+ "cosine pointing angle",
+ "on-the-fly status",
+ "dca v0 to pvtx",
+ "p",
+ "p_{T}",
+ "p_{x}",
+ "p_{y}",
+ "p_{z}"
+};
+
+//____________________________________________________________________
+AliESDv0Cuts::AliESDv0Cuts(const Char_t* name, const Char_t* title) : AliAnalysisCuts(name,title),
+  fCutMinDcaPosToVertex(0),
+  fCutMinDcaNegToVertex(0),
+  fCutMaxChi2(40),
+  fCutMaxDcaV0Daughters(0),
+  fCutMinRadius(0),
+  fCutMaxRadius(0),
+  fCutMinCosinePointingAngle(0),
+  fCutRequireOnFlyStatus(0),
+  fCutMaxDcaV0ToVertex(0),
+  fPMin(0),
+  fPMax(0),
+  fPtMin(0),
+  fPtMax(0),
+  fPxMin(0),
+  fPxMax(0),
+  fPyMin(0),
+  fPyMax(0),
+  fPzMin(0),
+  fPzMax(0),
+  fHistogramsOn(0),
+  fhCutStatistics(0),         
+  fhCutCorrelation(0)
+{
+  //
+  // constructor
+  //
+
+  Init();
+
+  //##############################################################################
+  // setting default cuts
+  SetMinDcaPosToVertex();
+  SetMinDcaNegToVertex();
+  SetMaxChi2();
+  SetMaxDcaV0Daughters();
+  SetMinRadius();
+  SetMaxRadius();
+  SetMinCosinePointingAngle();
+  SetRequireOnFlyStatus();
+  SetMaxDcaV0ToVertex();
+  SetPRange();
+  SetPtRange();
+  SetPxRange();
+  SetPyRange();
+  SetPzRange();
+
+  SetHistogramsOn();
+}
+
+//_____________________________________________________________________________
+AliESDv0Cuts::AliESDv0Cuts(const AliESDv0Cuts &c) : AliAnalysisCuts(c),
+  fCutMinDcaPosToVertex(0),
+  fCutMinDcaNegToVertex(0),
+  fCutMaxChi2(0),
+  fCutMaxDcaV0Daughters(0),
+  fCutMinRadius(0),
+  fCutMaxRadius(0),
+  fCutMinCosinePointingAngle(0),
+  fCutRequireOnFlyStatus(0),
+  fCutMaxDcaV0ToVertex(0),
+  fPMin(0),
+  fPMax(0),
+  fPtMin(0),
+  fPtMax(0),
+  fPxMin(0),
+  fPxMax(0),
+  fPyMin(0),
+  fPyMax(0),
+  fPzMin(0),
+  fPzMax(0),
+  fHistogramsOn(0),
+  fhCutStatistics(0),         
+  fhCutCorrelation(0)
+{
+  //
+  // copy constructor
+  //
+
+  ((AliESDv0Cuts &) c).Copy(*this);
+}
+
+AliESDv0Cuts::~AliESDv0Cuts()
+{
+  //
+  // destructor
+  //
+
+  for (Int_t i=0; i<2; i++) {
+    
+    if (fhDcaPosToVertex[i])
+      delete fhDcaPosToVertex[i];
+    if (fhDcaNegToVertex[i])
+      delete fhDcaNegToVertex[i];
+    if (fhChi2[i])
+      delete fhChi2[i];
+    if (fhDcaV0Daughters[i])
+      delete fhDcaV0Daughters[i];
+    if (fhRadius[i])
+      delete fhRadius[i];             
+    if (fhCosinePointingAngle[i])
+      delete fhCosinePointingAngle[i];
+    if (fhOnFlyStatus[i])
+    delete fhOnFlyStatus[i];
+    if (fhDcaV0ToVertex[i])
+      delete fhDcaV0ToVertex[i];
+    if (fhPt[i])
+      delete fhPt[i];
+  }
+
+  if (fhCutStatistics)
+    delete fhCutStatistics;             
+  if (fhCutCorrelation)
+    delete fhCutCorrelation;            
+}
+
+void AliESDv0Cuts::Init()
+{
+  //
+  // sets everything to zero
+  //
+  fCutMinDcaPosToVertex      = 0;
+  fCutMinDcaNegToVertex      = 0;
+  fCutMaxChi2                = 0;
+  fCutMaxDcaV0Daughters      = 0;
+  fCutMinRadius              = 0;
+  fCutMaxRadius              = 0;
+  fCutMinCosinePointingAngle = 0;
+  fCutRequireOnFlyStatus     = 0;
+  fCutMaxDcaV0ToVertex       = 0;
+
+  fPMin  = 0;
+  fPMax  = 0;
+  fPtMin = 0;
+  fPtMax = 0;
+  fPxMin = 0;
+  fPxMax = 0;
+  fPyMin = 0;
+  fPyMax = 0;
+  fPzMin = 0;
+  fPzMax = 0;
+
+  fHistogramsOn = kFALSE;
+
+  for (Int_t i=0; i<2; ++i)
+  {
+    fhDcaPosToVertex[i]      = 0;
+    fhDcaNegToVertex[i]      = 0;
+    fhChi2[i]                = 0;
+    fhDcaV0Daughters[i]      = 0;
+    fhRadius[i]              = 0;
+    fhCosinePointingAngle[i] = 0;
+    fhOnFlyStatus[i]         = 0;
+    fhDcaV0ToVertex[i]       = 0;
+    
+    fhPt[i]                  = 0;
+  }
+  fhCutStatistics = 0;
+  fhCutCorrelation = 0;
+}
+
+//_____________________________________________________________________________
+AliESDv0Cuts &AliESDv0Cuts::operator=(const AliESDv0Cuts &c)
+{
+  //
+  // Assignment operator
+  //
+
+  if (this != &c) ((AliESDv0Cuts &) c).Copy(*this);
+  return *this;
+}
+
+//_____________________________________________________________________________
+void AliESDv0Cuts::Copy(TObject &c) const
+{
+  //
+  // Copy function
+  //
+
+  AliESDv0Cuts& target = (AliESDv0Cuts &) c;
+
+  target.Init();
+
+  target.fCutMinDcaPosToVertex      = fCutMinDcaPosToVertex;
+  target.fCutMinDcaNegToVertex      = fCutMinDcaNegToVertex;
+  target.fCutMaxChi2                = fCutMaxChi2;
+  target.fCutMaxDcaV0Daughters      = fCutMaxDcaV0Daughters;
+  target.fCutMinRadius              = fCutMinRadius;
+  target.fCutMaxRadius              = fCutMaxRadius;
+  target.fCutMinCosinePointingAngle = fCutMinCosinePointingAngle;
+  target.fCutRequireOnFlyStatus     = fCutRequireOnFlyStatus;
+  target.fCutMaxDcaV0ToVertex       = fCutMaxDcaV0ToVertex;
+
+  target.fPMin  = fPMin;
+  target.fPMax  = fPMax;
+  target.fPtMin = fPtMin;
+  target.fPtMax = fPtMax;
+  target.fPxMin = fPxMin;
+  target.fPxMax = fPxMax;
+  target.fPyMin = fPyMin;
+  target.fPyMax = fPyMax;
+  target.fPzMin = fPzMin;
+  target.fPzMax = fPzMax;
+
+  target.fHistogramsOn = fHistogramsOn;
+
+  for (Int_t i=0; i<2; ++i)
+  {
+    if (fhDcaPosToVertex[i]) target.fhDcaPosToVertex[i] = (TH1F*) fhDcaPosToVertex[i]->Clone();
+    if (fhDcaNegToVertex[i]) target.fhDcaNegToVertex[i] = (TH1F*) fhDcaNegToVertex[i]->Clone();
+    if (fhChi2[i]) target.fhChi2[i] = (TH1F*) fhChi2[i]->Clone();
+    if (fhDcaV0Daughters[i]) target.fhDcaV0Daughters[i] = (TH1F*) fhDcaV0Daughters[i]->Clone();
+    if (fhRadius[i]) target.fhRadius[i] = (TH1F*) fhRadius[i]->Clone();
+    if (fhCosinePointingAngle[i]) target.fhCosinePointingAngle[i] = (TH1F*) fhCosinePointingAngle[i]->Clone();
+    if (fhOnFlyStatus[i]) target.fhOnFlyStatus[i] = (TH1F*) fhOnFlyStatus[i]->Clone();
+    if (fhDcaV0ToVertex[i]) target.fhDcaV0ToVertex[i] = (TH1F*) fhDcaV0ToVertex[i]->Clone();
+    
+    if (fhPt[i]) target.fhPt[i] = (TH1F*) fhPt[i]->Clone();
+  }
+  if (fhCutStatistics) target.fhCutStatistics = (TH1F*) fhCutStatistics->Clone();
+  if (fhCutCorrelation) target.fhCutCorrelation = (TH2F*) fhCutCorrelation->Clone();
+
+  TNamed::Copy(c);
+}
+
+//_____________________________________________________________________________
+Long64_t AliESDv0Cuts::Merge(TCollection* list) {
+  // Merge a list of AliESDv0Cuts objects with this (needed for PROOF)
+  // Returns the number of merged objects (including this)
+
+  if (!list)
+    return 0;
+  
+  if (list->IsEmpty())
+    return 1;
+
+  if (!fHistogramsOn)
+    return 0;
+
+  TIterator* iter = list->MakeIterator();
+  TObject* obj;
+
+
+  // collection of measured and generated histograms
+  Int_t count = 0;
+  while ((obj = iter->Next())) {
+
+    AliESDv0Cuts* entry = dynamic_cast<AliESDv0Cuts*>(obj);
+    if (entry == 0)
+      continue;
+
+    if (!entry->fHistogramsOn)
+      continue;
+    
+    for (Int_t i=0; i<2; i++) {
+      
+      fhDcaPosToVertex[i]     ->Add(entry->fhDcaPosToVertex[i]     );
+      fhDcaNegToVertex[i]     ->Add(entry->fhDcaNegToVertex[i]     );
+      fhChi2[i]               ->Add(entry->fhChi2[i]               );
+      fhDcaV0Daughters[i]     ->Add(entry->fhDcaV0Daughters[i]     );
+      fhRadius[i]             ->Add(entry->fhRadius[i]             );
+      fhCosinePointingAngle[i]->Add(entry->fhCosinePointingAngle[i]);
+      fhOnFlyStatus[i]        ->Add(entry->fhOnFlyStatus[i]        );
+      fhDcaV0ToVertex[i]      ->Add(entry->fhDcaV0ToVertex[i]      );
+
+      fhPt[i]                ->Add(entry->fhPt[i]);
+    }      
+    fhCutStatistics  ->Add(entry->fhCutStatistics);        
+    fhCutCorrelation ->Add(entry->fhCutCorrelation);      
+
+    count++;
+  }
+
+  return count+1;
+}
+
+void AliESDv0Cuts::EnableNeededBranches(TTree* tree)
+{
+  // enables the branches needed by AcceptV0, for a list see comment of AcceptV0
+
+  tree->SetBranchStatus("fV0s.fDcaV0Daughters", 1);
+  tree->SetBranchStatus("fV0s.fChi2V0", 1);
+  tree->SetBranchStatus("fV0s.fPos*", 1);
+  tree->SetBranchStatus("fV0s.fNmom*", 1);
+  tree->SetBranchStatus("fV0s.fPmom*", 1);
+  tree->SetBranchStatus("fV0s.fRr", 1);
+  tree->SetBranchStatus("fV0s.fPointAngle*", 1);
+  tree->SetBranchStatus("fV0s.fOnFlyStatus", 1);
+}
+
+//____________________________________________________________________
+Bool_t
+AliESDv0Cuts::IsSelected(TList* listObj) {
+  if(listObj->GetSize()!=4) return kFALSE;
+  AliESDv0           *esdV0     = (AliESDv0*)listObj->At(0);
+  AliESDtrack        *trackPos  = (AliESDtrack*)listObj->At(1);
+  AliESDtrack        *trackNeg  = (AliESDtrack*)listObj->At(2);
+  const AliESDVertex *esdVertex = (AliESDVertex*)listObj->At(3);
+  return AcceptV0(esdV0,trackPos,trackNeg,esdVertex);
+}
+
+//____________________________________________________________________
+Bool_t
+AliESDv0Cuts::AcceptV0(AliESDv0* esdV0, AliESDtrack* trackPos, AliESDtrack* trackNeg, const AliESDVertex* esdVertex) {
+  // 
+  // figure out if the v0s survives all the v0 cuts defined
+  //
+  // the different quality parameter and kinematic values are first
+  // retrieved from the v0. then it is found out what cuts the
+  // v0 did not survive and finally the cuts are imposed.
+
+  // this function needs the following branches... but this is not enough
+  // fV0s.fDcaV0Daughters
+  // fV0s.fChi2V0
+  // fV0s.fPos*
+  // fV0s.fNmom*
+  // fV0s.fPmom*
+  // fV0s.fRr
+  // fV0s.fPointAngle
+  // fV0s.fOnFlyStatus
+
+  Float_t  dcaPosToVertex = 0, dcaNegToVertex = 0;
+  Float_t  tdcaPosToVertex[2], tdcaNegToVertex[2];
+
+  if (trackPos) trackPos->GetImpactParameters(tdcaPosToVertex[0],tdcaPosToVertex[1]);
+  if (trackNeg) trackNeg->GetImpactParameters(tdcaNegToVertex[0],tdcaNegToVertex[1]);
+
+  dcaPosToVertex = TMath::Sqrt(tdcaPosToVertex[0]*tdcaPosToVertex[0]+tdcaPosToVertex[1]*tdcaPosToVertex[1]);
+  dcaNegToVertex = TMath::Sqrt(tdcaNegToVertex[0]*tdcaNegToVertex[0]+tdcaNegToVertex[1]*tdcaNegToVertex[1]);
+
+  UInt_t  status = esdV0->GetOnFlyStatus();
+  Float_t chi2  = esdV0->GetChi2V0();
+
+  Double_t dcaV0Daughters = esdV0->GetDcaV0Daughters();
+
+  Double_t vtxPosition[3]; esdVertex->GetXYZ(vtxPosition);
+  Double_t dcaV0ToVertex  = esdV0->GetD(vtxPosition[0],vtxPosition[1],vtxPosition[2]);
+
+  Double_t v0Position[3];
+  esdV0->GetXYZ(v0Position[0],v0Position[1],v0Position[2]);
+  Double_t radius = TMath::Sqrt(TMath::Power(v0Position[0],2) + TMath::Power(v0Position[1],2));
+  Double_t v0CosinePointingAngle = esdV0->GetV0CosineOfPointingAngle();
+
+  // getting the kinematic variables of the v0
+  Double_t p[3];
+  esdV0->GetPxPyPz(p[0],p[1],p[2]);
+  Float_t momentum = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2) + TMath::Power(p[2],2));
+  Float_t pt       = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2));
+
+  //########################################################################
+  // cut the v0?
+  
+  Bool_t cuts[kNCuts];
+  for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
+  
+  // v0 quality cuts
+  if (dcaPosToVertex < fCutMinDcaPosToVertex)
+    cuts[0]=kTRUE;
+  if (dcaNegToVertex < fCutMinDcaNegToVertex) 
+    cuts[1]=kTRUE;
+  if (chi2 > fCutMaxChi2) 
+    cuts[2]=kTRUE; 
+  if (dcaV0Daughters > fCutMaxDcaV0Daughters) 
+    cuts[3]=kTRUE;  
+  if (radius  < fCutMinRadius) 
+    cuts[4]=kTRUE;  
+  if (radius  > fCutMaxRadius) 
+    cuts[5]=kTRUE;  
+  if (v0CosinePointingAngle < fCutMinCosinePointingAngle)
+    cuts[6]=kTRUE;  
+  if (fCutRequireOnFlyStatus && !status)
+    cuts[7]=kTRUE;  
+  if (dcaV0ToVertex > fCutMaxDcaV0ToVertex)
+    cuts[8] = kTRUE;
+
+  // v0 kinematics cut
+  if((momentum < fPMin) || (momentum > fPMax)) 
+    cuts[9]=kTRUE;
+  if((pt < fPtMin) || (pt > fPtMax)) 
+    cuts[10] = kTRUE;
+  if((p[0] < fPxMin) || (p[0] > fPxMax)) 
+    cuts[11] = kTRUE;
+  if((p[1] < fPyMin) || (p[1] > fPyMax)) 
+    cuts[12] = kTRUE;
+  if((p[2] < fPzMin) || (p[2] > fPzMax))
+    cuts[13] = kTRUE;
+
+  Bool_t cut=kFALSE;
+  for (Int_t i=0; i<kNCuts; i++) 
+    if (cuts[i]) cut = kTRUE;
+  
+  //########################################################################
+  // filling histograms
+  if (fHistogramsOn) {
+    fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n v0s")));
+    
+    if (cut)
+      fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n cut v0s")));
+    
+    for (Int_t i=0; i<kNCuts; i++) {
+      if (cuts[i])
+       fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin(fgkCutNames[i])));
+      
+      for (Int_t j=i; j<kNCuts; j++) {
+       if (cuts[i] && cuts[j]) {
+         Float_t x = fhCutCorrelation->GetXaxis()->GetBinCenter(fhCutCorrelation->GetXaxis()->FindBin(fgkCutNames[i]));
+         Float_t y = fhCutCorrelation->GetYaxis()->GetBinCenter(fhCutCorrelation->GetYaxis()->FindBin(fgkCutNames[j]));
+         fhCutCorrelation->Fill(x,y);
+       }
+      }
+    }
+    
+    fhDcaPosToVertex[0]->Fill(dcaPosToVertex);
+    fhDcaNegToVertex[0]->Fill(dcaNegToVertex);
+    fhChi2[0]->Fill(chi2);
+    fhDcaV0Daughters[0]->Fill(dcaV0Daughters);
+    fhRadius[0]->Fill(radius);
+    fhCosinePointingAngle[0]->Fill(v0CosinePointingAngle);
+    fhOnFlyStatus[0]->Fill(status);
+    fhDcaV0ToVertex[0]->Fill(dcaV0ToVertex);
+    
+    fhPt[0]->Fill(pt);
+  }
+
+  //########################################################################
+  // cut the v0!
+  if (cut) return kFALSE;
+
+  //########################################################################
+  // filling histograms after cut
+  if (fHistogramsOn) {
+    fhDcaPosToVertex[1]->Fill(dcaPosToVertex);
+    fhDcaNegToVertex[1]->Fill(dcaNegToVertex);
+    fhChi2[1]->Fill(chi2);
+    fhDcaV0Daughters[1]->Fill(dcaV0Daughters);
+    fhRadius[1]->Fill(radius);
+    fhCosinePointingAngle[1]->Fill(v0CosinePointingAngle);
+    fhOnFlyStatus[1]->Fill(status);
+    fhDcaV0ToVertex[1]->Fill(dcaV0ToVertex);
+    
+    fhPt[1]->Fill(pt);
+  }
+
+  return kTRUE;
+}
+
+//____________________________________________________________________
+TObjArray* AliESDv0Cuts::GetAcceptedV0s(AliESD* esd)
+{
+  //
+  // returns an array of all v0s that pass the cuts
+  //
+
+  TObjArray* acceptedV0s = new TObjArray();
+  //  const AliESDVertex *spdVertex = esd->GetVertex();
+  const AliESDVertex *primaryVertex = esd->GetPrimaryVertex();
+  Int_t    lIndexTrackPos       = 0, lIndexTrackNeg       = 0;
+
+  // loop over esd v0s
+  for (Int_t iV0 = 0; iV0 < esd->GetNumberOfV0s(); iV0++) {
+    AliESDv0* v0 = esd->GetV0(iV0);
+
+    lIndexTrackPos = TMath::Abs(v0->GetPindex());
+    lIndexTrackNeg = TMath::Abs(v0->GetNindex());
+    AliESDtrack *trackPos = esd->GetTrack(lIndexTrackPos);
+    AliESDtrack *trackNeg = esd->GetTrack(lIndexTrackNeg);
+
+    if (AcceptV0(v0,trackPos,trackNeg,primaryVertex))
+      acceptedV0s->Add(v0);
+  }
+
+  return acceptedV0s;
+}
+
+//____________________________________________________________________
+Int_t AliESDv0Cuts::CountAcceptedV0s(AliESD* esd)
+{
+  //
+  // returns an the number of v0s that pass the cuts
+  //
+
+  Int_t count = 0;
+  //  const AliESDVertex *spdVertex = esd->GetVertex();
+  const AliESDVertex *primaryVertex = esd->GetPrimaryVertex();
+  Int_t    lIndexTrackPos       = 0, lIndexTrackNeg       = 0;
+
+  // loop over esd v0s
+  for (Int_t iV0 = 0; iV0 < esd->GetNumberOfV0s(); iV0++) {
+    AliESDv0* v0 = esd->GetV0(iV0);
+
+    lIndexTrackPos = TMath::Abs(v0->GetPindex());
+    lIndexTrackNeg = TMath::Abs(v0->GetNindex());
+    AliESDtrack *trackPos = esd->GetTrack(lIndexTrackPos);
+    AliESDtrack *trackNeg = esd->GetTrack(lIndexTrackNeg);
+
+    if (AcceptV0(v0,trackPos,trackNeg,primaryVertex))
+      count++;
+  }
+
+  return count;
+}
+
+//____________________________________________________________________
+TObjArray* AliESDv0Cuts::GetAcceptedV0s(AliESDEvent* esd)
+{
+  //
+  // returns an array of all v0s that pass the cuts
+  //
+
+  TObjArray* acceptedV0s = new TObjArray();
+  //  const AliESDVertex *spdVertex = esd->GetVertex();
+  const AliESDVertex *primaryVertex = esd->GetPrimaryVertex();
+  Int_t    lIndexTrackPos       = 0, lIndexTrackNeg       = 0;
+
+  // loop over esd v0s
+  for (Int_t iV0 = 0; iV0 < esd->GetNumberOfV0s(); iV0++) {
+    AliESDv0* v0 = esd->GetV0(iV0);
+
+    lIndexTrackPos = TMath::Abs(v0->GetPindex());
+    lIndexTrackNeg = TMath::Abs(v0->GetNindex());
+    AliESDtrack *trackPos = esd->GetTrack(lIndexTrackPos);
+    AliESDtrack *trackNeg = esd->GetTrack(lIndexTrackNeg);
+
+    if (AcceptV0(v0,trackPos,trackNeg,primaryVertex))
+      acceptedV0s->Add(v0);
+  }
+
+  return acceptedV0s;
+}
+
+//____________________________________________________________________
+Int_t AliESDv0Cuts::CountAcceptedV0s(AliESDEvent* esd)
+{
+  //
+  // returns an the number of v0s that pass the cuts
+  //
+
+  Int_t count = 0;
+  //  const AliESDVertex *spdVertex = esd->GetVertex();
+  const AliESDVertex *primaryVertex = esd->GetPrimaryVertex();
+  Int_t    lIndexTrackPos       = 0, lIndexTrackNeg       = 0;
+
+  // loop over esd v0s
+  for (Int_t iV0 = 0; iV0 < esd->GetNumberOfV0s(); iV0++) {
+    AliESDv0* v0 = esd->GetV0(iV0);
+
+    lIndexTrackPos = TMath::Abs(v0->GetPindex());
+    lIndexTrackNeg = TMath::Abs(v0->GetNindex());
+    AliESDtrack *trackPos = esd->GetTrack(lIndexTrackPos);
+    AliESDtrack *trackNeg = esd->GetTrack(lIndexTrackNeg);
+
+    if (AcceptV0(v0,trackPos,trackNeg,primaryVertex))
+      count++;
+  }
+
+  return count;
+}
+
+//____________________________________________________________________
+ void AliESDv0Cuts::DefineHistograms(Int_t color) {
+   // 
+   // diagnostics histograms are defined
+   // 
+
+   fHistogramsOn=kTRUE;
+   
+   Bool_t oldStatus = TH1::AddDirectoryStatus();
+   TH1::AddDirectory(kFALSE);
+   
+   //###################################################################################
+   // defining histograms
+
+   fhCutStatistics = new TH1F("cut_statistics","cut statistics",kNCuts+4,-0.5,kNCuts+3.5);
+
+   fhCutStatistics->GetXaxis()->SetBinLabel(1,"n v0s");
+   fhCutStatistics->GetXaxis()->SetBinLabel(2,"n cut v0s");
+
+   fhCutCorrelation = new TH2F("cut_correlation","cut correlation",kNCuts,-0.5,kNCuts-0.5,kNCuts,-0.5,kNCuts-0.5);;
+  
+   for (Int_t i=0; i<kNCuts; i++) {
+     fhCutStatistics->GetXaxis()->SetBinLabel(i+4,fgkCutNames[i]);
+     fhCutCorrelation->GetXaxis()->SetBinLabel(i+1,fgkCutNames[i]);
+     fhCutCorrelation->GetYaxis()->SetBinLabel(i+1,fgkCutNames[i]);
+   } 
+
+  fhCutStatistics  ->SetLineColor(color);
+  fhCutCorrelation ->SetLineColor(color);
+  fhCutStatistics  ->SetLineWidth(2);
+  fhCutCorrelation ->SetLineWidth(2);
+
+  Char_t str[256];
+  for (Int_t i=0; i<2; i++) {
+    if (i==0) sprintf(str," ");
+    else sprintf(str,"_cut");
+
+    fhDcaPosToVertex[i]      = new TH1F(Form("dcaPosToVertex%s",str),"",120,0,3);
+    fhDcaNegToVertex[i]      = new TH1F(Form("dcaNegToVertex%s",str),"",120,0,3);
+    fhChi2[i]                = new TH1F(Form("chi2%s",str),"",50,0,50);
+    fhDcaV0Daughters[i]      = new  TH1F(Form("dcaV0Daughters%s",str),"",200,0,5);
+    fhRadius[i]              = new  TH1F(Form("decayRadius%s",str),"",300,0,150);
+    fhCosinePointingAngle[i] = new  TH1F(Form("cosinePointingAngle%s",str),"",100,-1,1);
+    fhOnFlyStatus[i]         = new  TH1F(Form("onflystatus%s",str),"",5,0,5);
+    fhDcaV0ToVertex[i]       = new  TH1F(Form("dcaV0ToVertex%s",str),"",100,0,5);
+
+    fhPt[i]                  = new TH1F(Form("pt%s",str)     ,"p_{T} distribution;p_{T} (GeV/c)",500,0.0,100.0);
+    
+    fhDcaPosToVertex[i]->SetTitle("Dca of positive daughter to parent vertex");
+    fhDcaNegToVertex[i]->SetTitle("Dca of negative daughter to parent vertex");
+    fhChi2[i]->SetTitle("#Chi^{2} for v0");
+    fhDcaV0Daughters[i]->SetTitle("Dca between daughter tracks");
+    fhRadius[i]->SetTitle("Decay radius of the v0");
+    fhCosinePointingAngle[i]->SetTitle("Cosine of the Pointing Angle");
+    fhOnFlyStatus[i]->SetTitle("On-the-Fly Status");
+    fhDcaV0ToVertex[i]->SetTitle("Dca of v0 to parent vertex");
+
+    fhDcaPosToVertex[i]->SetLineColor(color);      fhDcaPosToVertex[i]->SetLineWidth(2);
+    fhDcaNegToVertex[i]->SetLineColor(color);      fhDcaNegToVertex[i]->SetLineWidth(2);
+    fhChi2[i]->SetLineColor(color);                fhChi2[i]->SetLineWidth(2);
+    fhDcaV0Daughters[i]->SetLineColor(color);      fhDcaV0Daughters[i]->SetLineWidth(2);
+    fhRadius[i]->SetLineColor(color);              fhRadius[i]->SetLineWidth(2);
+    fhCosinePointingAngle[i]->SetLineColor(color); fhCosinePointingAngle[i]->SetLineWidth(2);
+    fhOnFlyStatus[i]->SetLineColor(color);         fhOnFlyStatus[i]->SetLineWidth(2);
+    fhDcaV0ToVertex[i]->SetLineColor(color);       fhDcaV0ToVertex[i]->SetLineWidth(2); 
+  }
+
+  TH1::AddDirectory(oldStatus);
+}
+
+//____________________________________________________________________
+Bool_t AliESDv0Cuts::LoadHistograms(const Char_t* dir)
+{
+  //
+  // loads the histograms from a file
+  // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
+  //
+
+  if (!dir)
+    dir = GetName();
+
+  if (!gDirectory->cd(dir))
+    return kFALSE;
+
+  fhCutStatistics = dynamic_cast<TH1F*> (gDirectory->Get("cut_statistics"));
+  fhCutCorrelation = dynamic_cast<TH2F*> (gDirectory->Get("cut_correlation"));
+
+  Char_t str[5];
+  for (Int_t i=0; i<2; i++) {
+    if (i==0)
+    {
+      gDirectory->cd("before_cuts");
+      str[0] = 0;
+    }
+    else
+    {
+      gDirectory->cd("after_cuts");
+      sprintf(str,"_cut");
+    }
+
+    fhDcaPosToVertex[i]      = dynamic_cast<TH1F*> (gDirectory->Get(Form("dcaPosToVertex%s",str)     ));
+    fhDcaNegToVertex[i]      = dynamic_cast<TH1F*> (gDirectory->Get(Form("dcaNegToVertex%s",str)     ));
+    fhChi2[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("chi2%s",str)));
+    fhDcaV0Daughters[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("dcaV0Daughters%s",str)));
+    fhRadius[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("decayRadius%s",str)));
+    fhCosinePointingAngle[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("cosinepointingangle%s",str)));
+    fhOnFlyStatus[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("onflystatus%s",str)));
+    fhDcaV0ToVertex[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("dcaV0ToVertex%s",str)));
+
+    fhPt[i] = dynamic_cast<TH1F*> (gDirectory->Get(Form("pt%s",str)));
+
+    gDirectory->cd("../");
+  }
+
+  gDirectory->cd("..");
+
+  return kTRUE;
+}
+
+//____________________________________________________________________
+void AliESDv0Cuts::SaveHistograms(const Char_t* dir) {
+  //
+  // saves the histograms in a directory (dir)
+  //
+
+  if (!fHistogramsOn) {
+    AliDebug(0, "Histograms not on - cannot save histograms!!!");
+    return;
+  }
+
+  if (!dir)
+    dir = GetName();
+
+  gDirectory->mkdir(dir);
+  gDirectory->cd(dir);
+
+  gDirectory->mkdir("before_cuts");
+  gDirectory->mkdir("after_cuts");
+  fhCutStatistics->Write();
+  fhCutCorrelation->Write();
+
+  for (Int_t i=0; i<2; i++) {
+    if (i==0)
+      gDirectory->cd("before_cuts");
+    else
+      gDirectory->cd("after_cuts");
+
+    fhDcaPosToVertex[i]      ->Write();
+    fhDcaNegToVertex[i]      ->Write();
+    fhChi2[i]                ->Write();
+    fhDcaV0Daughters[i]      ->Write();
+    fhRadius[i]              ->Write();
+    fhCosinePointingAngle[i] ->Write();
+    fhOnFlyStatus[i]         ->Write();
+    fhDcaV0ToVertex[i]       ->Write();
+
+    fhPt[i]                  ->Write();
+    
+    gDirectory->cd("../");
+  }
+
+  gDirectory->cd("../");
+}
+
+//____________________________________________________________________
+void AliESDv0Cuts::DrawHistograms()
+{
+  // draws some histograms
+
+  TCanvas* canvas1 = new TCanvas(Form("%s_1", GetName()), "V0 Quality Results1", 800, 800);
+  canvas1->Divide(2, 2);
+
+  canvas1->cd(1);
+  fhDcaPosToVertex[0]->SetStats(kFALSE);
+  fhDcaPosToVertex[0]->Draw();
+
+  canvas1->cd(2);
+  fhChi2[0]->SetStats(kFALSE);
+  fhChi2[0]->Draw();
+
+  canvas1->cd(3);
+  fhDcaV0ToVertex[0]->SetStats(kFALSE);
+  fhDcaV0ToVertex[0]->Draw();
+
+  canvas1->SaveAs(Form("%s_%s.gif", GetName(), canvas1->GetName()));
+
+  TCanvas* canvas2 = new TCanvas(Form("%s_2", GetName()), "V0 Quality Results2", 1200, 800);
+  canvas2->Divide(3, 2);
+
+  canvas2->cd(1);
+  fhDcaV0Daughters[0]->SetStats(kFALSE);
+  gPad->SetLogy();
+  fhDcaV0Daughters[0]->Draw();
+
+  canvas2->cd(2);
+  fhRadius[0]->SetStats(kFALSE);
+  gPad->SetLogy();
+  fhRadius[0]->Draw();
+
+
+  canvas2->cd(4);
+  fhCosinePointingAngle[0]->SetStats(kFALSE);
+  gPad->SetLogy();
+  fhCosinePointingAngle[0]->Draw();
+
+  canvas2->cd(5);
+  fhOnFlyStatus[0]->SetStats(kFALSE);
+  gPad->SetLogy();
+  fhOnFlyStatus[0]->Draw();
+
+  canvas2->SaveAs(Form("%s_%s.gif", GetName(), canvas2->GetName()));
+
+  TCanvas* canvas3 = new TCanvas(Form("%s_4", GetName()), "V0 Quality Results3", 800, 500);
+  canvas3->Divide(2, 1);
+
+  canvas3->cd(1);
+  fhCutStatistics->SetStats(kFALSE);
+  fhCutStatistics->LabelsOption("v");
+  gPad->SetBottomMargin(0.3);
+  fhCutStatistics->Draw();
+
+  canvas3->cd(2);
+  fhCutCorrelation->SetStats(kFALSE);
+  fhCutCorrelation->LabelsOption("v");
+  gPad->SetBottomMargin(0.3);
+  gPad->SetLeftMargin(0.3);
+  fhCutCorrelation->Draw("COLZ");
+
+  canvas3->SaveAs(Form("%s_%s.gif", GetName(), canvas3->GetName()));
+}
+
diff --git a/ANALYSIS/AliESDv0Cuts.h b/ANALYSIS/AliESDv0Cuts.h
new file mode 100644 (file)
index 0000000..25bc86f
--- /dev/null
@@ -0,0 +1,121 @@
+//
+//  Class for handling of ESD v0 cuts.
+//
+//
+
+#ifndef ALIESDV0CUTS_H
+#define ALIESDV0CUTS_H
+
+#include <TF1.h>
+#include <TH2.h>
+#include "AliAnalysisCuts.h"
+
+class AliESD;
+class AliESDEvent;
+class AliESDVertex;
+class AliESDtrack;
+class AliESDv0;
+class AliLog;
+class TTree;
+
+class AliESDv0Cuts : public AliAnalysisCuts
+{
+public:
+  AliESDv0Cuts(const Char_t* name = "AliESDv0Cuts", const Char_t* title = "");
+  virtual ~AliESDv0Cuts();
+
+  Bool_t IsSelected(TObject* /*obj*/) {return kTRUE;}
+  Bool_t IsSelected(TList* listObj);
+  Bool_t IsSelected(TObject* obj1, TObject* obj2, TObject* obj3, TObject* obj4)
+  {return AcceptV0((AliESDv0*) obj1, (AliESDtrack*) obj2, (AliESDtrack*) obj3, (const AliESDVertex*) obj4);}
+  Bool_t AcceptV0(AliESDv0* esdV0, AliESDtrack* trackPos, AliESDtrack* trackNeg, const AliESDVertex* esdVertex);
+  TObjArray* GetAcceptedV0s(AliESD* esd);
+  Int_t CountAcceptedV0s(AliESD* esd);
+  TObjArray* GetAcceptedV0s(AliESDEvent* esd);
+  Int_t CountAcceptedV0s(AliESDEvent* esd);
+
+  virtual Long64_t Merge(TCollection* list);
+  virtual void Copy(TObject &c) const;
+  AliESDv0Cuts(const AliESDv0Cuts& pd);  // Copy Constructor
+  AliESDv0Cuts &operator=(const AliESDv0Cuts &c);
+
+  //######################################################
+  // v0 quality cut setters  
+  void SetMinDcaPosToVertex(Float_t min=-1)          {fCutMinDcaPosToVertex=min;}
+  void SetMinDcaNegToVertex(Float_t min=-1)          {fCutMinDcaNegToVertex=min;}
+  void SetMaxChi2(Float_t max=1e10)                  {fCutMaxChi2=max;}
+  void SetMaxDcaV0Daughters(Float_t max=1e10)        {fCutMaxDcaV0Daughters=max;}
+  void SetMinRadius(Float_t min=-1)                  {fCutMinRadius=min;}
+  void SetMaxRadius(Float_t max=1e10)                {fCutMaxRadius=max;}
+  void SetMinCosinePointingAngle(Float_t min=-1)     {fCutMinCosinePointingAngle=min;}
+  void SetRequireOnFlyStatus(Bool_t b=kFALSE)        {fCutRequireOnFlyStatus=b;}
+  void SetMaxDcaV0ToVertex(Float_t max=1e10)         {fCutMaxDcaV0ToVertex=max;}
+
+  // v0 kinematic cut setters
+  void SetPRange(Float_t r1=0, Float_t r2=1e10)      {fPMin=r1;   fPMax=r2;}
+  void SetPtRange(Float_t r1=0, Float_t r2=1e10)     {fPtMin=r1;  fPtMax=r2;}
+  void SetPxRange(Float_t r1=-1e10, Float_t r2=1e10) {fPxMin=r1;  fPxMax=r2;}
+  void SetPyRange(Float_t r1=-1e10, Float_t r2=1e10) {fPyMin=r1;  fPyMax=r2;}
+  void SetPzRange(Float_t r1=-1e10, Float_t r2=1e10) {fPzMin=r1;  fPzMax=r2;}
+
+  //######################################################
+  void SetHistogramsOn(Bool_t b=kFALSE) {fHistogramsOn = b;}
+  void DefineHistograms(Int_t color=1);
+  virtual Bool_t LoadHistograms(const Char_t* dir = 0);
+  void SaveHistograms(const Char_t* dir = 0);
+  void DrawHistograms();
+
+  static void EnableNeededBranches(TTree* tree);
+
+  // void SaveQualityCuts(Char_t* file)
+  // void LoadQualityCuts(Char_t* file)
+
+protected:
+  void Init(); // sets everything to 0
+
+  enum { kNCuts = 14 };
+
+  //######################################################
+  // esd v0 quality cuts
+  static const Char_t* fgkCutNames[kNCuts]; //! names of cuts (for internal use)
+
+  Float_t fCutMinDcaPosToVertex;      // min dca of the positive daughter to the primary vertex
+  Float_t fCutMinDcaNegToVertex;      // min dca of the negative daughter to the primary vertex
+  Float_t fCutMaxChi2;                // max chi2
+  Float_t fCutMaxDcaV0Daughters;      // max dca between the two v0 daughters
+  Float_t fCutMinRadius;              // min reconstruction radius (fiducial volume)
+  Float_t fCutMaxRadius;              // max reconstruction radius (fiducial volume)
+  Float_t fCutMinCosinePointingAngle; // min cosine of pointing angle
+  Bool_t  fCutRequireOnFlyStatus;     // require on fly status
+  Float_t fCutMaxDcaV0ToVertex;       // max dca of the v0 to the primary vertex
+
+  // v0 kinematics cuts
+  Float_t fPMin,   fPMax;             // definition of the range of the P
+  Float_t fPtMin,  fPtMax;            // definition of the range of the Pt
+  Float_t fPxMin,  fPxMax;            // definition of the range of the Px
+  Float_t fPyMin,  fPyMax;            // definition of the range of the Py
+  Float_t fPzMin,  fPzMax;            // definition of the range of the Pz
+
+  //######################################################
+  // diagnostics histograms
+  Bool_t fHistogramsOn;               // histograms on/off
+
+  TH1F* fhDcaPosToVertex[2];          //->
+  TH1F* fhDcaNegToVertex[2];          //->
+  TH1F* fhChi2[2];                    //->
+  TH1F* fhDcaV0Daughters[2];          //->
+  TH1F* fhRadius[2];                  //->
+  TH1F* fhCosinePointingAngle[2];     //->
+  TH1F* fhOnFlyStatus[2];             //->
+  TH1F* fhDcaV0ToVertex[2];           //->
+  
+  TH1F* fhPt[2];                      //-> pt of esd v0s
+
+  TH1F* fhCutStatistics;              //-> statistics of what cuts the v0s did not survive
+  TH2F* fhCutCorrelation;             //-> 2d statistics plot
+
+  ClassDef(AliESDv0Cuts, 1)
+};
+
+
+#endif
index 32b7152daea4a684a0091f4480ad4327586ccc6e..a8b9355b2e77270b1e50e2ca561323f594c37d26 100644 (file)
@@ -1,7 +1,7 @@
 #-*- Mode: Makefile -*-
 
 SRCS =         AliAnalysisTaskSE.cxx AliAnalysisTaskESDfilter.cxx AliAnalysisTaskKineFilter.cxx \
-       AliKineTrackCuts.cxx AliESDtrackCuts.cxx
+       AliKineTrackCuts.cxx AliESDtrackCuts.cxx AliESDv0Cuts.cxx
 
 CHECKALIEN = $(shell root-config --has-alien)
 ifeq (yes,$(CHECKALIEN))