]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG4/JetTasks/AliAnalysisTaskJetServices.cxx
Change selection bit pattern to the same conventions as used in ANALYSIS/macros/...
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetServices.cxx
index aafbd5f44b1c5fea8758dfbffadae057833f084d..7e4017077c77fac85472edc9314fc0b95cee2aa1 100644 (file)
@@ -22,6 +22,7 @@
  
 #include <TROOT.h>
 #include <TRandom.h>
+#include <TString.h>
 #include <TSystem.h>
 #include <TInterpreter.h>
 #include <TChain.h>
@@ -34,7 +35,7 @@
 #include <TList.h>
 #include <TLorentzVector.h>
 #include <TClonesArray.h>
-#include  "TDatabasePDG.h"
+#include "TDatabasePDG.h"
 
 #include "AliAnalysisTaskJetServices.h"
 #include "AliAnalysisDataContainer.h"
@@ -59,7 +60,7 @@
 #include "AliGenCocktailEventHeader.h"
 #include "AliInputEventHandler.h"
 #include "AliPhysicsSelection.h"
-
+#include "AliTriggerAnalysis.h"
 
 #include "AliAnalysisHelperJetTasks.h"
 
@@ -67,20 +68,33 @@ ClassImp(AliAnalysisTaskJetServices)
 
 AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(): AliAnalysisTaskSE(),
   fUseAODInput(kFALSE),
+  fUsePhysicsSelection(kFALSE),
+  fMC(kFALSE),
+  fSelectionInfoESD(0),
+  fEventCutInfoESD(0),
   fAvgTrials(1),
-  fZVtxCut(8.),
-  fRealData(kFALSE),
-  fPhysicsSelection(0x0),
+  fVtxXMean(0),
+  fVtxYMean(0),
+  fVtxZMean(0),
+  fVtxRCut(1.),
+  fVtxZCut(8.),
+  fPtMinCosmic(5.),
+  fRIsolMinCosmic(3.),
+  fMaxCosmicAngle(0.01),
   fh1Xsec(0x0),
   fh1Trials(0x0),
   fh1PtHard(0x0),
   fh1PtHardTrials(0x0),
+  fh1SelectionInfoESD(0x0),
+  fh1EventCutInfoESD(0),
   fh2TriggerCount(0x0),
   fh2ESDTriggerCount(0x0),
   fh2TriggerVtx(0x0),
   fh2ESDTriggerVtx(0x0),
   fh2ESDTriggerRun(0x0),
   fh2VtxXY(0x0),
+  fh1NCosmicsPerEvent(0x0),
+  fTriggerAnalysis(0x0),
   fHistList(0x0)  
 {
   fRunRange[0] = fRunRange[1] = 0; 
@@ -89,20 +103,33 @@ AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(): AliAnalysisTaskSE(),
 AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(const char* name):
   AliAnalysisTaskSE(name),
   fUseAODInput(kFALSE),
+  fUsePhysicsSelection(kFALSE),
+  fMC(kFALSE),
+  fSelectionInfoESD(0),
+  fEventCutInfoESD(0),
   fAvgTrials(1),
-  fZVtxCut(8.),
-  fRealData(kFALSE),
-  fPhysicsSelection(0x0),
+  fVtxXMean(0),
+  fVtxYMean(0),
+  fVtxZMean(0),
+  fVtxRCut(1.),
+  fVtxZCut(8.),
+  fPtMinCosmic(5.),
+  fRIsolMinCosmic(3.),
+  fMaxCosmicAngle(0.01),
   fh1Xsec(0x0),
   fh1Trials(0x0),
   fh1PtHard(0x0),
   fh1PtHardTrials(0x0),
+  fh1SelectionInfoESD(0x0),
+  fh1EventCutInfoESD(0),
   fh2TriggerCount(0x0),
   fh2ESDTriggerCount(0x0),
   fh2TriggerVtx(0x0),
   fh2ESDTriggerVtx(0x0),
   fh2ESDTriggerRun(0x0),
   fh2VtxXY(0x0),
+  fh1NCosmicsPerEvent(0x0),
+  fTriggerAnalysis(0x0),
   fHistList(0x0)  
 {
   fRunRange[0] = fRunRange[1] = 0; 
@@ -157,13 +184,6 @@ void AliAnalysisTaskJetServices::UserCreateOutputObjects()
 
   Bool_t oldStatus = TH1::AddDirectoryStatus();
   TH1::AddDirectory(kFALSE);
-
-  if(!fPhysicsSelection)
-    fPhysicsSelection = new AliPhysicsSelection();
-  fPhysicsSelection->SetName("AliPhysicsSelection_outputlist"); // to prevent conflict with object that is automatically streamed back
-  //AliLog::SetClassDebugLevel("AliPhysicsSelection", AliLog::kDebug);
-  fHistList->Add(fPhysicsSelection);
-
   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
   fHistList->Add(fh1Xsec);
@@ -200,6 +220,13 @@ void AliAnalysisTaskJetServices::UserCreateOutputObjects()
   fHistList->Add(fh1PtHard);
   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
   fHistList->Add(fh1PtHardTrials);
+  fh1SelectionInfoESD = new TH1F("fh1SelectionInfoESD","Bit Masks that satisfy the selection info",
+                                AliAnalysisHelperJetTasks::kTotalSelections,0.5,AliAnalysisHelperJetTasks::kTotalSelections+0.5);
+  fHistList->Add(fh1SelectionInfoESD);
+
+  fh1EventCutInfoESD = new TH1F("fh1EventCutInfoESD","Bit Masks that for the events after each step of cuts",
+                               kTotalEventCuts,0.5,kTotalEventCuts+0.5);
+  fHistList->Add(fh1EventCutInfoESD);
 
   // 3 decisions, 0 trigger X, X + SPD vertex, X + SPD vertex in range  
   // 3 triggers BB BE/EB EE
@@ -220,6 +247,9 @@ void AliAnalysisTaskJetServices::UserCreateOutputObjects()
     if(hn)hn->Sumw2();
   }
 
+  fh1NCosmicsPerEvent = new TH1F("fh1NCosmicsPerEvent","Number of cosmic candidates per event",10,0.,10.);
+  fHistList->Add(fh1NCosmicsPerEvent),
+
 
   TH1::AddDirectory(oldStatus);
 }
@@ -229,8 +259,6 @@ void AliAnalysisTaskJetServices::Init()
   //
   // Initialization
   //
-
-  Printf(">>> AnalysisTaskJetServices::Init() debug level %d\n",fDebug);
   if (fDebug > 1) printf("AnalysisTaskJetServices::Init() \n");
 
 }
@@ -244,8 +272,12 @@ void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
  
   AliAODEvent *aod = 0;
   AliESDEvent *esd = 0;
-
+  
   AliAnalysisHelperJetTasks::Selected(kTRUE,kFALSE); // set slection to false
+  fSelectionInfoESD = 0; // reset
+  fEventCutInfoESD = 0; // reset
+  AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD); // set slection to false
+
 
   if(fUseAODInput){    
     aod = dynamic_cast<AliAODEvent*>(InputEvent());
@@ -265,25 +297,55 @@ void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
     esd = dynamic_cast<AliESDEvent*>(InputEvent());
   }
   
+  fSelectionInfoESD |= kNoEventCut;
+  fEventCutInfoESD |= kNoEventCut;
+
+  Bool_t esdVtxValid = false;
+  Bool_t esdVtxIn = false;
+  Bool_t aodVtxValid = false;
+  Bool_t aodVtxIn = false;
+
   if(esd){
+    // trigger analyisis
+    if(!fTriggerAnalysis){
+      fTriggerAnalysis = new AliTriggerAnalysis;
+      fTriggerAnalysis->SetAnalyzeMC(fMC);
+      fTriggerAnalysis->SetSPDGFOThreshhold(1);
+    }
+    //    fTriggerAnalysis->FillTriggerClasses(esd);
+    Bool_t v0A       = fTriggerAnalysis->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0A);
+    Bool_t v0C       = fTriggerAnalysis->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0C);
+    Bool_t v0ABG = fTriggerAnalysis->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0ABG);
+    Bool_t v0CBG = fTriggerAnalysis->IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0CBG);
+    Bool_t spdFO      = fTriggerAnalysis->SPDFiredChips(esd, 0);;
+    if(v0A)fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kV0A;
+    if(v0C)fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kV0C;
+    if(!(v0ABG||v0CBG))fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kNoV0BG;
+    if(spdFO)fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kSPDFO;
+
     Float_t run = (Float_t)esd->GetRunNumber();
     const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
+    esdVtxValid = 
+    esdVtxIn = IsVertexIn(vtxESD);
     Float_t zvtx = -999;
     Float_t xvtx = -999;
     Float_t yvtx = -999;
 
-    if(vtxESD->GetNContributors()>0){
+    if(esdVtxValid){
       zvtx = vtxESD->GetZ();
       yvtx = vtxESD->GetY();
       xvtx = vtxESD->GetX();
     }
 
+    // CKB this can be cleaned up a bit...
+    
     Int_t iTrig = -1;
     if(esd->GetFiredTriggerClasses().Contains("CINT1B")
        ||esd->GetFiredTriggerClasses().Contains("CSMBB")
        ||esd->GetFiredTriggerClasses().Contains("MB1")
        ||esd->GetFiredTriggerClasses().Contains("CINT6B")){
       iTrig = 0;
+      fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kBunchBunch;
     }
     else if(esd->GetFiredTriggerClasses().Contains("CINT1A")
            ||esd->GetFiredTriggerClasses().Contains("CSMBA")
@@ -293,94 +355,130 @@ void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
            ||esd->GetFiredTriggerClasses().Contains("CINT6C")){
       // empty bunch or bunch empty
       iTrig = 1;
+      fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kBunchEmpty;
     }
-    if(esd->GetFiredTriggerClasses().Contains("CINT1-E")
+    else if(esd->GetFiredTriggerClasses().Contains("CINT1-E")
        ||esd->GetFiredTriggerClasses().Contains("CINT6-E")){
       iTrig = 2;
+      fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kEmptyEmpty;
     }
 
     
     if(iTrig>=0){
       iTrig *= 3;
       fh2ESDTriggerRun->Fill(run,iTrig+1);
-      if(vtxESD->GetNContributors()>0){
+      if(vtxESD->GetNContributors()>2){
        fh2ESDTriggerRun->Fill(run,iTrig+2);
        fh2VtxXY->Fill(xvtx,yvtx);
       }
-      if(TMath::Abs(zvtx)<fZVtxCut&&TMath::Abs(xvtx)<0.5&&TMath::Abs(yvtx)<0.5)fh2ESDTriggerRun->Fill(run,iTrig+3);
+      xvtx -= fVtxXMean; 
+      yvtx -= fVtxYMean; 
+      zvtx -= fVtxZMean; 
+      Float_t r2 = xvtx *xvtx + yvtx *yvtx; 
+      if(TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut))fh2ESDTriggerRun->Fill(run,iTrig+3);
     }
     else{
       fh2ESDTriggerRun->Fill(run,0);
     }
+    // BKC
+  }
+  
+
+  // Apply additional constraints
+  Bool_t esdEventSelected = IsEventSelected(esd);
+  Bool_t esdEventPileUp = IsEventPileUp(esd);
+  Bool_t esdEventCosmic = IsEventCosmic(esd);
+
+  Bool_t aodEventSelected = IsEventSelected(aod);
 
+  Bool_t physicsSelection = ((fInputHandler->IsEventSelected())&AliVEvent::kMB);
+  fEventCutInfoESD |= kPhysicsSelectionCut; // other alreay set via IsEventSelected
+  fh1EventCutInfoESD->Fill(fEventCutInfoESD);
 
+  if(esdEventSelected) fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kVertexIn;
+  if(esdEventPileUp)   fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kIsPileUp;
+  if(esdEventCosmic)   fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kIsCosmic;
+  if(physicsSelection) fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kPhysicsSelection;
+
+
+  // here we have all selection information, fill histogram
+  for(unsigned int i = 1;i<(UInt_t)fh1SelectionInfoESD->GetNbinsX();i++){
+    if((i&fSelectionInfoESD)==i)fh1SelectionInfoESD->Fill(i);
   }
-  
-  // loop over all possible trigger and 
-  for(int it = AliAnalysisHelperJetTasks::kAcceptAll;it < AliAnalysisHelperJetTasks::kTrigger;it++){
-    Bool_t esdTrig = kFALSE;
-    if(esd){
+  AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD); 
+
+  if(esd&&aod&&fDebug){
+    if(esdEventSelected&&!aodEventSelected){
+      Printf("%s:%d Different Selection for ESD and AOD",(char*)__FILE__,__LINE__);
+      const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
+      const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
+      Printf("ESD Vtx %s %s %d",vtxESD->GetName(),vtxESD->GetTitle(),vtxESD->GetNContributors());
+      vtxESD->Print();
+      Printf("AOD Vtx %s %s %d",vtxAOD->GetName(),vtxAOD->GetTitle(),vtxAOD->GetNContributors());
+      vtxAOD->Print();
+    }
+  }
+
+  // loop over all possible triggers for esd 
+
+  if(esd){
+    const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
+    //      Printf(">> ESDvtx %s %s",vtxESD->GetName(),vtxESD->GetTitle());vtxESD->Print();
+    TString vtxName(vtxESD->GetName());
+    for(int it = AliAnalysisHelperJetTasks::kAcceptAll;it < AliAnalysisHelperJetTasks::kTrigger;it++){
+      Bool_t esdTrig = kFALSE;
       esdTrig = AliAnalysisHelperJetTasks::IsTriggerFired(esd,(AliAnalysisHelperJetTasks::Trigger)it);
       if(esdTrig)fh2ESDTriggerCount->Fill(it,kAllTriggered);
+      Bool_t cand = physicsSelection;
+      if(cand){
+       fh2ESDTriggerCount->Fill(it,kSelectedALICE); 
+      }
+      if(!fUsePhysicsSelection)cand =  AliAnalysisHelperJetTasks::IsTriggerFired(esd,AliAnalysisHelperJetTasks::kMB1);
+      if(vtxESD->GetNContributors()>2&&!vtxName.Contains("TPCVertex")){
+       if(esdTrig)fh2ESDTriggerCount->Fill(it,kTriggeredVertex);
+       Float_t zvtx = vtxESD->GetZ();
+       if(esdEventSelected&&esdTrig){
+         fh2ESDTriggerCount->Fill(it,kTriggeredVertexIn);
+         fh2ESDTriggerVtx->Fill(it,zvtx);
+       }
+       if(cand)fh2ESDTriggerCount->Fill(it,kSelectedALICEVertexValid);
+      }
+      if(cand&&esdEventSelected){
+       fh2ESDTriggerCount->Fill(it,kSelectedALICEVertexIn);
+       fh2ESDTriggerCount->Fill(it,kSelected);
+       AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
+      }
     }
-    Bool_t aodTrig = kFALSE;
-    if(aod){
+  }
+
+  if(aod){
+    const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
+    aodVtxValid = IsVertexValid(vtxAOD);
+    aodVtxIn = IsVertexIn(vtxAOD);
+
+    for(int it = AliAnalysisHelperJetTasks::kAcceptAll;it < AliAnalysisHelperJetTasks::kTrigger;it++){
+      Bool_t aodTrig = kFALSE;
       aodTrig = AliAnalysisHelperJetTasks::IsTriggerFired(aod,(AliAnalysisHelperJetTasks::Trigger)it);
+      Bool_t cand = aod->GetHeader()->GetOfflineTrigger()&AliVEvent::kMB;
       if(aodTrig)fh2TriggerCount->Fill(it,kAllTriggered);
-    }
-
-    // Check wether we have also an SPD vertex
-    
-    if(aod){
-      const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
-      //      Printf(">> AODvtx %s %s",vtxAOD->GetName(),vtxAOD->GetTitle());vtxAOD->Print();
-      
-      if(vtxAOD->GetNContributors()>0){
-       if(aodTrig)fh2TriggerCount->Fill(it,kTriggeredSPDVertex);
+      if(aodVtxValid){
+       if(aodTrig)fh2TriggerCount->Fill(it,kTriggeredVertex);
        Float_t zvtx = vtxAOD->GetZ();
-       Float_t yvtx = vtxAOD->GetY();
-       Float_t xvtx = vtxAOD->GetX();
-       fh2TriggerVtx->Fill(it,zvtx);
-       if(TMath::Abs(zvtx)<fZVtxCut&&aodTrig&&TMath::Abs(xvtx)<0.5&&TMath::Abs(yvtx)<0.5){
+       if(cand&&aodEventSelected){
+         fh2TriggerCount->Fill(it,kSelected);
+       }
+       if(aodTrig&&aodEventSelected){
+         fh2TriggerVtx->Fill(it,zvtx);
          fh2TriggerCount->Fill(it,kTriggeredVertexIn);
        }
-      }
-    }
-    if(esd){
-      const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
-      //      Printf(">> ESDvtx %s %s",vtxESD->GetName(),vtxESD->GetTitle());vtxESD->Print();
-      Bool_t cand = fPhysicsSelection->IsCollisionCandidate(esd);
-      if(cand) fh2ESDTriggerCount->Fill(it,kSelectedALICE);
-      if(vtxESD->GetNContributors()>0){
-       if(esdTrig)fh2ESDTriggerCount->Fill(it,kTriggeredSPDVertex);
-       Float_t zvtx = vtxESD->GetZ();
-       Float_t yvtx = vtxESD->GetY();
-       Float_t xvtx = vtxESD->GetX();
-       fh2ESDTriggerVtx->Fill(it,zvtx);
-       if(TMath::Abs(zvtx)<fZVtxCut&&esdTrig&&TMath::Abs(xvtx)<0.5&&TMath::Abs(yvtx)<0.5){
-         fh2ESDTriggerCount->Fill(it,kTriggeredVertexIn);
-         // here we select based on ESD info...
-         if(fRealData){
-           if(cand){
-             fh2ESDTriggerCount->Fill(it,kSelected);
-             AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
-           }
-         }
-         else{
-           fh2ESDTriggerCount->Fill(it,kSelected);
-           AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
-         }
+       if(cand){
+         fh2TriggerCount->Fill(it,kSelectedALICEVertexValid);
+         if(aodEventSelected)fh2TriggerCount->Fill(it,kSelectedALICEVertexIn);
        }
-
-
       }
-
     }
-
   }
 
-
-
   if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);
 
   
@@ -416,49 +514,196 @@ void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
   PostData(1, fHistList);
 }
 
+Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliESDEvent* esd){
+  if(!esd)return kFALSE;
+  const AliESDVertex *vtx = esd->GetPrimaryVertex();
+  return IsVertexIn(vtx); // vertex in calls vertex valid
+}
 
-void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
-{
-  // Terminate analysis
-  //
 
-  TDirectory* owd = gDirectory;
+Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliAODEvent* aod) const {
+  if(!aod)return kFALSE;
+  const AliAODVertex *vtx = aod->GetPrimaryVertex();
+  return IsVertexIn(vtx); // VertexIn calls VertexValid
+}
 
-  if (fDebug > 1) printf("AnalysisJetServices: Terminate() \n");
+Bool_t  AliAnalysisTaskJetServices::IsVertexValid ( const AliESDVertex* vtx) {
+
+  // We can check the type of the vertex by its name and title
+  // if vertexer failed title is empty (default c'tor) and ncontributors is 0
+  // vtx       name                  title
+  // Tracks    PrimaryVertex         VertexerTracksNoConstraint
+  // Tracks    PrimaryVertex         VertexerTracksWithConstraint
+  // TPC       TPCVertex             VertexerTracksNoConstraint
+  // TPC       TPCVertex             VertexerTracksWithConstraint
+  // SPD       SPDVertex             vertexer: 3D
+  // SPD       SPDVertex             vertexer: Z
+  
+  Int_t nCont = vtx->GetNContributors();
+
+  if(nCont>=1){
+    fEventCutInfoESD |= kContributorsCut1;    
+    if(nCont>=2){
+      fEventCutInfoESD |= kContributorsCut2;    
+      if(nCont>=3){
+       fEventCutInfoESD |= kContributorsCut3;    
+      }
+    }
+  }
+  
+  if(nCont<3)return kFALSE;
 
-  fHistList = dynamic_cast<TList*> (GetOutputData(1));
-  if (!fHistList)
-    Printf("ERROR: fHistList not available");
+  // do not want tpc only primary vertex
+  TString vtxName(vtx->GetName());
+  if(vtxName.Contains("TPCVertex")){
+    fEventCutInfoESD |= kVertexTPC;
+    return kFALSE;
+  }
+  if(vtxName.Contains("SPDVertex"))fEventCutInfoESD |= kVertexSPD;
+  if(vtxName.Contains("PrimaryVertex"))fEventCutInfoESD |= kVertexGlobal;
 
 
+  TString vtxTitle(vtx->GetTitle());
+  if(vtxTitle.Contains("vertexer: Z")){
+    if(vtx->GetDispersion()>0.02)return kFALSE;   
+  }
+  fEventCutInfoESD |= kSPDDispersionCut;
+  return kTRUE;
+}
 
-  AliAnalysisDataContainer *cont = GetOutputSlot(1)->GetContainer();
-  TString filename = cont->GetFileName();
-  TFile *f = NULL;
-  // Check first if the file is already opened
-  f = (TFile*)gROOT->GetListOfFiles()->FindObject(filename);
-  if (f) {
-    // Cd to file
-    f->cd();
-    // Check for a folder request
-    TString dir = cont->GetFolderName(); 
-    if (!dir.IsNull()) {
-      if (!f->GetDirectory(dir)) f->mkdir(dir);
-      f->cd(dir);
-    }
+
+Bool_t  AliAnalysisTaskJetServices::IsVertexValid ( const AliAODVertex* vtx) const {
+
+  // We can check the type of the vertex by its name and title
+  // if vertexer failed title is empty (default c'tor) and ncontributors is 0
+  // vtx       name                  title
+  // Tracks    PrimaryVertex         VertexerTracksNoConstraint
+  // TPC       TPCVertex             VertexerTracksNoConstraint
+  // SPD       SPDVertex             vertexer: 3D
+  // SPD       SPDVertex             vertexer: Z
+  
+  if(vtx->GetNContributors()<3)return kFALSE;
+  // do not want tpc only primary vertex
+  TString vtxName(vtx->GetName());
+  if(vtxName.Contains("TPCVertex"))return kFALSE;
+
+  // no dispersion yet...
+  /* 
+  TString vtxTitle(vtx->GetTitle());
+  if(vtxTitle.Contains("vertexer: Z")){
+    if(vtx->GetDispersion()>0.02)return kFALSE;
   }
+  */
+  return kTRUE;
+}
+
+
+Bool_t  AliAnalysisTaskJetServices::IsVertexIn (const AliESDVertex* vtx) {
+
+  if(!IsVertexValid(vtx))return kFALSE;
+
+  Float_t zvtx = vtx->GetZ();
+  Float_t yvtx = vtx->GetY();
+  Float_t xvtx = vtx->GetX();
+
+  xvtx -= fVtxXMean;
+  yvtx -= fVtxYMean;
+  zvtx -= fVtxZMean;
 
-  if (fHistList)
-  {
-    fPhysicsSelection = dynamic_cast<AliPhysicsSelection*> (fHistList->FindObject("AliPhysicsSelection_outputlist"));
+
+
+  if(TMath::Abs(zvtx)>fVtxZCut){
+    return kFALSE;
   }
+  fEventCutInfoESD |= kVertexZCut;  
+  Float_t r2   = yvtx*yvtx+xvtx*xvtx;      
+  if(r2>(fVtxRCut*fVtxRCut)){
+    return kFALSE;
+  }
+  fEventCutInfoESD |= kVertexRCut;  
+  return kTRUE;
+}
+
+
+Bool_t  AliAnalysisTaskJetServices::IsVertexIn (const AliAODVertex* vtx) const {
 
-  if (fPhysicsSelection)
-    {
-      fPhysicsSelection->SaveHistograms("physics_selection");
-      fPhysicsSelection->Print();
+  if(!IsVertexValid(vtx))return kFALSE;
+
+  Float_t zvtx = vtx->GetZ();
+  Float_t yvtx = vtx->GetY();
+  Float_t xvtx = vtx->GetX();
+
+  xvtx -= fVtxXMean;
+  yvtx -= fVtxYMean;
+  zvtx -= fVtxZMean;
+
+  Float_t r2   = yvtx*yvtx+xvtx*xvtx;      
+
+  Bool_t vertexIn = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
+  return vertexIn;
+}
+
+Bool_t AliAnalysisTaskJetServices::IsEventPileUp(const AliESDEvent* esd) const{
+  if(!esd)return kFALSE;
+  return esd->IsPileupFromSPD();
+}
+
+Bool_t AliAnalysisTaskJetServices::IsEventCosmic(const AliESDEvent* esd) const {
+  if(!esd)return kFALSE;
+  // add track cuts for which we look for cosmics...
+
+  Bool_t isCosmic = kFALSE;
+  Int_t nTracks = esd->GetNumberOfTracks();
+  Int_t nCosmicCandidates = 0;
+
+  for (Int_t iTrack1 = 0; iTrack1 < nTracks; iTrack1++) {
+    AliESDtrack* track1 = (AliESDtrack*)esd->GetTrack(iTrack1);
+    if (!track1)  continue;
+    UInt_t status1 = track1->GetStatus();
+    //If track is ITS stand alone track, skip the track
+    if (((status1 & AliESDtrack::kITSin) == 0 || (status1 & AliESDtrack::kTPCin))) continue;
+    if(track1->Pt()<fPtMinCosmic) continue;
+    //Start 2nd track loop to look for correlations
+    for (Int_t iTrack2 = iTrack1+1; iTrack2 < nTracks; iTrack2++) {
+      AliESDtrack* track2 = (AliESDtrack*)esd->GetTrack(iTrack2);
+      if(!track2) continue;
+      UInt_t status2 = track2->GetStatus();
+      //If track is ITS stand alone track, skip the track
+      if (((status2 & AliESDtrack::kITSin) == 0 || (status2 & AliESDtrack::kTPCin))) continue;
+      if(track2->Pt()<fPtMinCosmic) continue;
+      //Check if back-to-back
+      Double_t mom1[3],mom2[3];
+      track1->GetPxPyPz(mom1);
+      track2->GetPxPyPz(mom2);
+      TVector3 momv1(mom1[0],mom1[1],mom1[2]);
+      TVector3 momv2(mom2[0],mom2[1],mom2[2]);
+      Float_t theta = (float)(momv1.Phi()-momv2.Phi());
+      if(theta<-0.5*TMath::Pi()) theta+=2.*TMath::Pi();
+
+      Float_t deltaPhi = track1->Phi()-track2->Phi();
+      if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
+
+      Float_t rIsol = (float)(TMath::Sqrt( deltaPhi*deltaPhi+(track1->Eta()-track2->Eta())*(track1->Eta()-track2->Eta()) ));
+      if(rIsol<fRIsolMinCosmic) continue;
+
+      if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) {
+       nCosmicCandidates+=1;
+       isCosmic = kTRUE;
+      }
+      
     }
+  }
+
+  fh1NCosmicsPerEvent->Fill((float)nCosmicCandidates);
+
+  return isCosmic;
+}
 
-  owd->cd();
 
+
+
+void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
+{
+  // Terminate analysis
+  //
 }