]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG4/JetTasks/AliAnalysisTaskJetServices.cxx
updates for differnt jet configs (including background), reduced number of UA1 Jets
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetServices.cxx
index 78b15e204b9e1a8742f02d9fe59ba2c796fd4cd6..64ace5ade4a13d85f9d55c499a8b631214fef5d6 100644 (file)
@@ -38,6 +38,7 @@
 #include "TDatabasePDG.h"
 
 #include "AliAnalysisTaskJetServices.h"
+#include "AliESDCentrality.h"
 #include "AliAnalysisDataContainer.h"
 #include "AliAnalysisDataSlot.h"
 #include "AliAnalysisManager.h"
 #include "AliGenCocktailEventHeader.h"
 #include "AliInputEventHandler.h"
 #include "AliPhysicsSelection.h"
-
+#include "AliTriggerAnalysis.h"
 
 #include "AliAnalysisHelperJetTasks.h"
 
 ClassImp(AliAnalysisTaskJetServices)
 
+AliAODHeader*  AliAnalysisTaskJetServices::fgAODHeader = NULL;
+TClonesArray*   AliAnalysisTaskJetServices::fgAODVertices = NULL;
+
 AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(): AliAnalysisTaskSE(),
   fUseAODInput(kFALSE),
   fUsePhysicsSelection(kFALSE),
-  fRealData(kFALSE),
+  fMC(kFALSE),
+  fFilterAODCollisions(kFALSE),
+  fPhysicsSelectionFlag(AliVEvent::kMB),
   fSelectionInfoESD(0),
+  fEventCutInfoESD(0),
   fAvgTrials(1),
   fVtxXMean(0),
   fVtxYMean(0),
@@ -80,11 +87,13 @@ AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(): AliAnalysisTaskSE(),
   fPtMinCosmic(5.),
   fRIsolMinCosmic(3.),
   fMaxCosmicAngle(0.01),
+  fNonStdFile(""),
   fh1Xsec(0x0),
   fh1Trials(0x0),
   fh1PtHard(0x0),
   fh1PtHardTrials(0x0),
   fh1SelectionInfoESD(0x0),
+  fh1EventCutInfoESD(0),
   fh2TriggerCount(0x0),
   fh2ESDTriggerCount(0x0),
   fh2TriggerVtx(0x0),
@@ -92,6 +101,7 @@ AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(): AliAnalysisTaskSE(),
   fh2ESDTriggerRun(0x0),
   fh2VtxXY(0x0),
   fh1NCosmicsPerEvent(0x0),
+  fTriggerAnalysis(0x0),
   fHistList(0x0)  
 {
   fRunRange[0] = fRunRange[1] = 0; 
@@ -101,8 +111,11 @@ AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(const char* name):
   AliAnalysisTaskSE(name),
   fUseAODInput(kFALSE),
   fUsePhysicsSelection(kFALSE),
-  fRealData(kFALSE),
+  fMC(kFALSE),
+  fFilterAODCollisions(kFALSE),
+  fPhysicsSelectionFlag(AliVEvent::kMB),
   fSelectionInfoESD(0),
+  fEventCutInfoESD(0),
   fAvgTrials(1),
   fVtxXMean(0),
   fVtxYMean(0),
@@ -112,11 +125,13 @@ AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(const char* name):
   fPtMinCosmic(5.),
   fRIsolMinCosmic(3.),
   fMaxCosmicAngle(0.01),
+  fNonStdFile(""),
   fh1Xsec(0x0),
   fh1Trials(0x0),
   fh1PtHard(0x0),
   fh1PtHardTrials(0x0),
   fh1SelectionInfoESD(0x0),
+  fh1EventCutInfoESD(0),
   fh2TriggerCount(0x0),
   fh2ESDTriggerCount(0x0),
   fh2TriggerVtx(0x0),
@@ -124,6 +139,7 @@ AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(const char* name):
   fh2ESDTriggerRun(0x0),
   fh2VtxXY(0x0),
   fh1NCosmicsPerEvent(0x0),
+  fTriggerAnalysis(0x0),
   fHistList(0x0)  
 {
   fRunRange[0] = fRunRange[1] = 0; 
@@ -202,11 +218,11 @@ void AliAnalysisTaskJetServices::UserCreateOutputObjects()
 
   fh2ESDTriggerCount = new TH2F("fh2ESDTriggerCount",";Trigger No.;constrained;Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,kConstraints,-0.5,kConstraints-0.5); 
   fHistList->Add(fh2ESDTriggerCount);
-
-  fh2TriggerVtx = new TH2F("fh2TriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.); 
+  const Int_t nBins = AliAnalysisHelperJetTasks::kTrigger*kConstraints;
+  fh2TriggerVtx = new TH2F("fh2TriggerVtx",";Constraint No. * (trig no+1);Vtx (cm);Count",nBins,-0.5,nBins-0.5,400,-20.,20.); 
   fHistList->Add(fh2TriggerVtx);
 
-  fh2ESDTriggerVtx = new TH2F("fh2ESDTriggerVtx",";Trigger No.;Vtx (cm);Count",AliAnalysisHelperJetTasks::kTrigger,-0.5,AliAnalysisHelperJetTasks::kTrigger-0.5,400,-20.,20.); 
+  fh2ESDTriggerVtx = new TH2F("fh2ESDTriggerVtx",";Constraint No.* (trg no+1);Vtx (cm);Count",nBins,-0.5,nBins-0.5,400,-20.,20.); 
   fHistList->Add(fh2ESDTriggerVtx);
   
 
@@ -216,12 +232,16 @@ void AliAnalysisTaskJetServices::UserCreateOutputObjects()
   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
 
-  fh2ESDTriggerRun = new TH2F("fh2ESDTriggerRun","Trigger vs run number:run;trigger",(Int_t)(1+fRunRange[1]-fRunRange[0]),fRunRange[0]-0.5,fRunRange[1]+0.5,10,-0.5,9.5);
+  fh2ESDTriggerRun = new TH2F("fh2ESDTriggerRun","Eventclass vs run number:run;trigger",(Int_t)(1+fRunRange[1]-fRunRange[0]),fRunRange[0]-0.5,fRunRange[1]+0.5,10,-0.5,9.5);
   fHistList->Add(fh2ESDTriggerRun);
 
   fh2VtxXY = new TH2F("fh2VtxXY","Beam Spot all INT triggered events;x (cm);y (cm)",160,-10,10,160,-10,10);
@@ -242,6 +262,19 @@ void AliAnalysisTaskJetServices::UserCreateOutputObjects()
 
 
   TH1::AddDirectory(oldStatus);
+
+  // Add an AOD branch for replication
+  if(fNonStdFile.Length()){
+     if (fDebug > 1) AliInfo("Replicating header");
+     fgAODHeader = new AliAODHeader;
+     AddAODBranch("AliAODHeader",&fgAODHeader,fNonStdFile.Data());
+     if (fDebug > 1) AliInfo("Replicating primary vertices");
+
+     fgAODVertices = new TClonesArray("AliAODVertex",3);
+     fgAODVertices->SetName("vertices");
+     AddAODBranch("TClonesArray",&fgAODVertices,fNonStdFile.Data());
+
+  }
 }
 
 void AliAnalysisTaskJetServices::Init()
@@ -263,11 +296,19 @@ void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
   AliAODEvent *aod = 0;
   AliESDEvent *esd = 0;
 
+
+  
   AliAnalysisHelperJetTasks::Selected(kTRUE,kFALSE); // set slection to false
+  AliAnalysisHelperJetTasks::EventClass(kTRUE,0);
   fSelectionInfoESD = 0; // reset
+  fEventCutInfoESD = 0; // reset
   AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD); // set slection to false
 
 
+  static AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
+   
+
+
   if(fUseAODInput){    
     aod = dynamic_cast<AliAODEvent*>(InputEvent());
     if(!aod){
@@ -285,76 +326,56 @@ void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
     }
     esd = dynamic_cast<AliESDEvent*>(InputEvent());
   }
-  
-  fSelectionInfoESD |= AliAnalysisHelperJetTasks::kNone;
+  if(aod&&fDebug>2){
+    aod->Print();
+    Printf("Vertices %d",aod->GetNumberOfVertices());
+    Printf("tracks %d",aod->GetNumberOfTracks());
+    Printf("jets %d",aod->GetNJets());
+  }
+  fSelectionInfoESD |= kNoEventCut;
+  fEventCutInfoESD |= kNoEventCut;
+
+  Bool_t esdVtxValid = false;
+  Bool_t esdVtxIn = false;
+  Bool_t aodVtxValid = false;
+  Bool_t aodVtxIn = false;
 
   if(esd){
-    Float_t run = (Float_t)esd->GetRunNumber();
-    const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
-    Float_t zvtx = -999;
-    Float_t xvtx = -999;
-    Float_t yvtx = -999;
-
-    if(vtxESD->GetNContributors()>2){
-      zvtx = vtxESD->GetZ();
-      yvtx = vtxESD->GetY();
-      xvtx = vtxESD->GetX();
+    // 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;
+  }
+  // Apply additional constraints
+  Bool_t esdEventSelected = IsEventSelected(esd);
+  Bool_t esdEventPileUp = IsEventPileUp(esd);
+  Bool_t esdEventCosmic = IsEventCosmic(esd);
 
-    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")
-           ||esd->GetFiredTriggerClasses().Contains("CINT6A")
-           ||esd->GetFiredTriggerClasses().Contains("CINT1C")
-           ||esd->GetFiredTriggerClasses().Contains("CSMBC")
-           ||esd->GetFiredTriggerClasses().Contains("CINT6C")){
-      // empty bunch or bunch empty
-      iTrig = 1;
-      fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kBunchEmpty;
-    }
-    else if(esd->GetFiredTriggerClasses().Contains("CINT1-E")
-       ||esd->GetFiredTriggerClasses().Contains("CINT6-E")){
-      iTrig = 2;
-      fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kEmptyEmpty;
-    }
+  Bool_t aodEventSelected = IsEventSelected(aod);
 
-    
-    if(iTrig>=0){
-      iTrig *= 3;
-      fh2ESDTriggerRun->Fill(run,iTrig+1);
-      if(vtxESD->GetNContributors()>2){
-       fh2ESDTriggerRun->Fill(run,iTrig+2);
-       fh2VtxXY->Fill(xvtx,yvtx);
-      }
-      xvtx -= fVtxXMean; 
-      yvtx -= fVtxYMean; 
-      zvtx -= fVtxZMean; 
-      Float_t r2 = xvtx *xvtx + yvtx *yvtx; 
-      if(TMath::Abs(zvtx)<fVtxZCut&&r2<fVtxRCut)fh2ESDTriggerRun->Fill(run,iTrig+3);
+  Bool_t physicsSelection = ((fInputHandler->IsEventSelected())&fPhysicsSelectionFlag);
+  if(aodH&&physicsSelection&&fFilterAODCollisions&&fgAODHeader){
+    if(fgAODHeader->GetCentrality()<=80){
+      aodH->SetFillAOD(kTRUE);
     }
-    else{
-      fh2ESDTriggerRun->Fill(run,0);
-    }
-
-
   }
-  
-
-  // Apply additional constraints
-  Bool_t esdEventSelected = IsEventSelectedESD(esd);
-  Bool_t esdEventPileUp = IsEventPileUpESD(esd);
-  Bool_t esdEventCosmic = IsEventCosmicESD(esd);
 
-  Bool_t aodEventSelected = IsEventSelectedAOD(aod);
 
-  Bool_t physicsSelection = fInputHandler->IsEventSelected();
+  fEventCutInfoESD |= kPhysicsSelectionCut; // other alreay set via IsEventSelected
+  fh1EventCutInfoESD->Fill(fEventCutInfoESD);
 
   if(esdEventSelected) fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kVertexIn;
   if(esdEventPileUp)   fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kIsPileUp;
@@ -366,7 +387,7 @@ void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
   for(unsigned int i = 1;i<(UInt_t)fh1SelectionInfoESD->GetNbinsX();i++){
     if((i&fSelectionInfoESD)==i)fh1SelectionInfoESD->Fill(i);
   }
-  AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD); // set slection to false
+  AliAnalysisHelperJetTasks::SelectInfo(kTRUE,fSelectionInfoESD); 
 
   if(esd&&aod&&fDebug){
     if(esdEventSelected&&!aodEventSelected){
@@ -384,50 +405,87 @@ void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
 
   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 = fInputHandler->IsEventSelected();
-      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);
-       }
+    esdVtxValid = IsVertexValid(vtxESD);
+    esdVtxIn = IsVertexIn(vtxESD);
+    Float_t zvtx = vtxESD->GetZ();
+    Int_t  iCl = GetEventClass(esd);
+    AliAnalysisHelperJetTasks::EventClass(kTRUE,iCl);
+    Bool_t cand = physicsSelection;
+
+    if(fDebug)Printf("%s:%d %d %d %d Icl %d",(char*)__FILE__,__LINE__,esdVtxValid,esdVtxIn,cand,iCl);
+
+    if(cand){
+      fh2ESDTriggerCount->Fill(0.,kSelectedALICE); 
+      fh2ESDTriggerCount->Fill(iCl,kSelectedALICE); 
+      fh2ESDTriggerVtx->Fill(kSelectedALICE*(iCl+1),zvtx);
+    }
+    //    if(!fUsePhysicsSelection)cand =  AliAnalysisHelperJetTasks::IsTriggerFired(esd,AliAnalysisHelperJetTasks::kMB1);
+    if(esdVtxValid){
+      fh2ESDTriggerCount->Fill(0.,kTriggeredVertex);
+      fh2ESDTriggerCount->Fill(iCl,kTriggeredVertex);
+      fh2ESDTriggerVtx->Fill(iCl,zvtx);
+      if(esdVtxIn){
+       fh2ESDTriggerCount->Fill(0.,kTriggeredVertexIn);
+       fh2ESDTriggerCount->Fill(iCl,kTriggeredVertexIn);
+       fh2ESDTriggerVtx->Fill(kTriggeredVertexIn*(iCl+1),zvtx);
       }
-      if(cand&&esdEventSelected){
-       fh2ESDTriggerCount->Fill(it,kSelectedALICEVertexIn);
-       fh2ESDTriggerCount->Fill(it,kSelected);
-       AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
+      if(cand){
+       fh2ESDTriggerCount->Fill(0.,kSelectedALICEVertexValid);
+       fh2ESDTriggerCount->Fill(iCl,kSelectedALICEVertexValid);
+       fh2ESDTriggerVtx->Fill(kSelectedALICEVertexValid*(iCl+1),zvtx);
       }
     }
+    if(cand&&esdVtxIn){
+      fh2ESDTriggerCount->Fill(0.,kSelectedALICEVertexIn);
+      fh2ESDTriggerCount->Fill(iCl,kSelectedALICEVertexIn);
+      fh2ESDTriggerVtx->Fill(kSelectedALICEVertexIn*(iCl+1),zvtx);
+      fh2ESDTriggerVtx->Fill(kSelected*(iCl+1),zvtx);
+      fh2ESDTriggerCount->Fill(iCl,kSelected);
+      fh2ESDTriggerCount->Fill(0.,kSelected);
+      AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
+    }
   }
 
+
+
   if(aod){
     const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
-    //      Printf(">> AODvtx %s %s",vtxAOD->GetName(),vtxAOD->GetTitle());vtxAOD->Print();
-    TString vtxTitle(vtxAOD->GetTitle());
-    for(int it = AliAnalysisHelperJetTasks::kAcceptAll;it < AliAnalysisHelperJetTasks::kTrigger;it++){
-      Bool_t aodTrig = kFALSE;
-      aodTrig = AliAnalysisHelperJetTasks::IsTriggerFired(aod,(AliAnalysisHelperJetTasks::Trigger)it);
-      if(aodTrig)fh2TriggerCount->Fill(it,kAllTriggered);
-      if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
-       if(aodTrig)fh2TriggerCount->Fill(it,kTriggeredVertex);
-       Float_t zvtx = vtxAOD->GetZ();
-       if(aodTrig&&aodEventSelected){
-         fh2TriggerVtx->Fill(it,zvtx);
-         fh2TriggerCount->Fill(it,kTriggeredVertexIn);
-       }
+    aodVtxValid = IsVertexValid(vtxAOD);
+    aodVtxIn = IsVertexIn(vtxAOD);
+    Float_t zvtx = vtxAOD->GetZ();
+    Int_t  iCl = GetEventClass(aod);
+    AliAnalysisHelperJetTasks::EventClass(kTRUE,iCl);
+    Bool_t cand = aod->GetHeader()->GetOfflineTrigger()&fPhysicsSelectionFlag;
+    if(fDebug)Printf("%s:%d AOD selection %d %d",(char*)__FILE__,__LINE__,cand,aod->GetHeader()->GetOfflineTrigger());
+    if(cand){
+      fh2TriggerCount->Fill(0.,kSelectedALICE); 
+      fh2TriggerCount->Fill(iCl,kSelectedALICE); 
+      fh2TriggerVtx->Fill(kSelectedALICE*(iCl+1),zvtx);
+    }
+    if(aodVtxValid){
+      fh2TriggerCount->Fill(0.,kTriggeredVertex);
+      fh2TriggerCount->Fill(iCl,kTriggeredVertex);
+      fh2TriggerVtx->Fill(iCl,zvtx);
+      if(aodVtxIn){
+       fh2TriggerCount->Fill(0.,kTriggeredVertexIn);
+       fh2TriggerCount->Fill(iCl,kTriggeredVertexIn);
+       fh2TriggerVtx->Fill(kTriggeredVertexIn*(iCl+1),zvtx);
+      }
+      if(cand){
+       fh2TriggerCount->Fill(0.,kSelectedALICEVertexValid);
+       fh2TriggerCount->Fill(iCl,kSelectedALICEVertexValid);
+       fh2TriggerVtx->Fill(kSelectedALICEVertexValid*(iCl+1),zvtx);
       }
     }
+    if(cand&&aodVtxIn){
+      fh2TriggerCount->Fill(0.,kSelectedALICEVertexIn);
+      fh2TriggerCount->Fill(iCl,kSelectedALICEVertexIn);
+      fh2TriggerVtx->Fill(kSelectedALICEVertexIn*(iCl+1),zvtx);
+      fh2TriggerVtx->Fill(kSelected*(iCl+1),zvtx);
+      fh2TriggerCount->Fill(iCl,kSelected);
+      fh2TriggerCount->Fill(0.,kSelected);
+      if(fUseAODInput)AliAnalysisHelperJetTasks::Selected(kTRUE,kTRUE);// select this event
+    }
   }
 
   if (fDebug > 1)printf(" AliAnalysisTaskJetServices: Analysing event # %5d\n", (Int_t) fEntry);
@@ -461,33 +519,188 @@ void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
 
   // trigger selection
   
+  // replication of 
+  if(fNonStdFile.Length()&&aod){
+    if (fgAODHeader){
+      *fgAODHeader =  *(dynamic_cast<AliAODHeader*>(aod->GetHeader()));
+    }
+    if(fgAODVertices){
+      Printf("%p",fgAODVertices);
+      Printf("%p",fgAODVertices->At(0));
+      fgAODVertices->Delete();
+      TClonesArray &vertices = *fgAODVertices;
+      const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
+      // we only use some basic information, 
+      
 
+      Double_t pos[3];
+      Double_t covVtx[6];
+
+      vtxAOD->GetXYZ(pos); // position                                                                
+      vtxAOD->GetCovMatrix(covVtx); //covariance matrix                                                      
+      Int_t jVertices = 0;
+      AliAODVertex * vtx = new(vertices[jVertices++])
+        AliAODVertex(pos, covVtx, vtxAOD->GetChi2perNDF(), NULL, -1, AliAODVertex::kPrimary);
+      vtx->SetName(vtxAOD->GetName());
+      vtx->SetTitle(vtxAOD->GetTitle());
+
+      TString vtitle = vtxAOD->GetTitle();
+      if (!vtitle.Contains("VertexerTracks"))
+       vtx->SetNContributors(vtxAOD->GetNContributors());
+
+      // Add SPD "main" vertex                                                                    
+      const AliAODVertex *vtxS = aod->GetPrimaryVertexSPD();
+      vtxS->GetXYZ(pos); // position
+      vtxS->GetCovMatrix(covVtx); //covariance matrix                                                      
+      AliAODVertex * mVSPD = new(vertices[jVertices++])
+       AliAODVertex(pos, covVtx, vtxS->GetChi2perNDF(), NULL, -1, AliAODVertex::kMainSPD);
+      mVSPD->SetName(vtxS->GetName());
+      mVSPD->SetTitle(vtxS->GetTitle());
+      mVSPD->SetNContributors(vtxS->GetNContributors());
+      
+      // Add tpc only vertex
+      if(esd){
+       const AliESDVertex *vtxT =  esd->GetPrimaryVertexTPC();
+       vtxT->GetXYZ(pos); // position                                                    
+       vtxT->GetCovMatrix(covVtx); //covariance matrix                                               
+       AliAODVertex * mVTPC = new(vertices[jVertices++])
+         AliAODVertex(pos, covVtx, vtxT->GetChi2toNDF(), NULL, -1, AliAODVertex::kMainTPC);
+       mVTPC->SetName(vtxT->GetName());
+       mVTPC->SetTitle(vtxT->GetTitle());
+       mVTPC->SetNContributors(vtxT->GetNContributors());
+      }
+    }
+  }
+  
   PostData(1, fHistList);
 }
 
-Bool_t AliAnalysisTaskJetServices::IsEventSelectedESD(AliESDEvent* esd){
+Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliESDEvent* esd){
   if(!esd)return kFALSE;
   const AliESDVertex *vtx = esd->GetPrimaryVertex();
+  return IsVertexIn(vtx); // vertex in calls vertex valid
+}
+
+AliAnalysisTaskJetServices::~AliAnalysisTaskJetServices(){
+  if(fgAODVertices){
+    fgAODVertices->Delete();
+    delete fgAODVertices;
+  }
+  if(fgAODHeader)delete fgAODHeader;
+}
+
+
+Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliAODEvent* aod) const {
+  if(!aod)return kFALSE;
+  const AliAODVertex *vtx = aod->GetPrimaryVertex();
+  return IsVertexIn(vtx); // VertexIn calls VertexValid
+}
+
+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;
 
+  // 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;
 
-  if(vtx->GetNContributors()<3)return kFALSE;
 
+  TString vtxTitle(vtx->GetTitle());
+  if(vtxTitle.Contains("vertexer: Z")){
+    if(vtx->GetDispersion()>0.02)return kFALSE;   
+  }
+  fEventCutInfoESD |= kSPDDispersionCut;
+  return kTRUE;
+}
+
+
+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();
 
-  // here we may fill the histos for run dependence....
+  xvtx -= fVtxXMean;
+  yvtx -= fVtxYMean;
+  zvtx -= fVtxZMean;
+
+
+
+  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(!IsVertexValid(vtx))return kFALSE;
+
+  Float_t zvtx = vtx->GetZ();
+  Float_t yvtx = vtx->GetY();
+  Float_t xvtx = vtx->GetX();
 
   xvtx -= fVtxXMean;
   yvtx -= fVtxYMean;
@@ -495,16 +708,16 @@ Bool_t AliAnalysisTaskJetServices::IsEventSelectedESD(AliESDEvent* esd){
 
   Float_t r2   = yvtx*yvtx+xvtx*xvtx;      
 
-  Bool_t eventSel = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
-  return eventSel;
+  Bool_t vertexIn = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
+  return vertexIn;
 }
 
-Bool_t AliAnalysisTaskJetServices::IsEventPileUpESD(AliESDEvent* esd){
+Bool_t AliAnalysisTaskJetServices::IsEventPileUp(const AliESDEvent* esd) const{
   if(!esd)return kFALSE;
   return esd->IsPileupFromSPD();
 }
 
-Bool_t AliAnalysisTaskJetServices::IsEventCosmicESD(AliESDEvent* esd){
+Bool_t AliAnalysisTaskJetServices::IsEventCosmic(const AliESDEvent* esd) const {
   if(!esd)return kFALSE;
   // add track cuts for which we look for cosmics...
 
@@ -515,11 +728,17 @@ Bool_t AliAnalysisTaskJetServices::IsEventCosmicESD(AliESDEvent* esd){
   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];
@@ -537,7 +756,7 @@ Bool_t AliAnalysisTaskJetServices::IsEventCosmicESD(AliESDEvent* esd){
       if(rIsol<fRIsolMinCosmic) continue;
 
       if(TMath::Abs(TMath::Pi()-theta)<fMaxCosmicAngle) {
-       nCosmicCandidates+=1.;
+       nCosmicCandidates+=1;
        isCosmic = kTRUE;
       }
       
@@ -550,45 +769,35 @@ Bool_t AliAnalysisTaskJetServices::IsEventCosmicESD(AliESDEvent* esd){
 }
 
 
-Bool_t AliAnalysisTaskJetServices::IsEventSelectedAOD(AliAODEvent* aod){
-  if(!aod)return kFALSE;
-  const AliAODVertex *vtx = aod->GetPrimaryVertex();
-  // 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
-  
+Int_t AliAnalysisTaskJetServices::GetEventClass(AliESDEvent *esd){
 
+  Float_t cent = 999;
+  if(esd->GetCentrality()){
+    cent = esd->GetCentrality()->GetCentralityPercentile("V0M");
+  }
+  if(cent>50)return 4;
+  if(cent>30)return 3;
+  if(cent>10)return 2;
+  return 1;
 
-  if(vtx->GetNContributors()<3)return kFALSE;
+}
 
-  // do not want tpc only primary vertex
-  TString vtxName(vtx->GetName());
-  if(vtxName.Contains("TPCVertex"))return kFALSE;
 
-  Float_t zvtx = vtx->GetZ();
-  Float_t yvtx = vtx->GetY();
-  Float_t xvtx = vtx->GetX();
+Int_t AliAnalysisTaskJetServices::GetEventClass(AliAODEvent *aod){
 
-  // here we may fill the histos for run dependence....
+  Float_t cent = aod->GetHeader()->GetCentrality();
 
-  xvtx -= fVtxXMean;
-  yvtx -= fVtxYMean;
-  zvtx -= fVtxZMean;
+  if(cent>50)return 4;
+  if(cent>30)return 3;
+  if(cent>10)return 2;
+  return 1;
 
-  Float_t r2   = yvtx*yvtx+xvtx*xvtx;      
-  Bool_t eventSel = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
-  return eventSel;
 }
 
 
-
-
 void AliAnalysisTaskJetServices::Terminate(Option_t */*option*/)
 {
   // Terminate analysis
   //
 }
+