Add possibility to enable trigger patch QA in the jet preparation

[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskJetServices.cxx

diff --git a/PWGJE/AliAnalysisTaskJetServices.cxx b/PWGJE/AliAnalysisTaskJetServices.cxx
index 8d41a5c85c95ffa921da90909bf0cc927f53c229..fa939322071ec4b1609f91eec334746b71b90da4 100644 (file)
--- a/PWGJE/AliAnalysisTaskJetServices.cxx
+++ b/PWGJE/AliAnalysisTaskJetServices.cxx
@@ -93,6 +93,7 @@ AliAnalysisTaskJetServices::AliAnalysisTaskJetServices():
   fVtxZMean(0),
   fVtxRCut(1.),
   fVtxZCut(8.),
+  fVtxDeltaZCut(0.1),
   fPtMinCosmic(5.),
   fRIsolMinCosmic(3.),
   fMaxCosmicAngle(0.01),
@@ -108,14 +109,15 @@ AliAnalysisTaskJetServices::AliAnalysisTaskJetServices():
   fTriggerName(0x0),
   fh1Xsec(0x0),
   fh1Trials(0x0),
+  fh1AvgTrials(0x0),
   fh1PtHard(0x0),
   fh1PtHardTrials(0x0),
   fh1SelectionInfoESD(0x0),
   fh1EventCutInfoESD(0),
   fh1CentralityESD(0),
   fh1Centrality(0),
-  fh1ReducedTrigger(0),
   fh1RP(0),
+  fh2ReducedTrigger(0),
   fh2CentralityTriggerESD(0),
   fh2CentralityTrigger(0),
   fh2TriggerCount(0x0),
@@ -166,6 +168,7 @@ AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(const char* name):
   fVtxZMean(0),
   fVtxRCut(1.),
   fVtxZCut(8.),
+  fVtxDeltaZCut(0.1),
   fPtMinCosmic(5.),
   fRIsolMinCosmic(3.),
   fMaxCosmicAngle(0.01),
@@ -181,14 +184,15 @@ AliAnalysisTaskJetServices::AliAnalysisTaskJetServices(const char* name):
   fTriggerName(0x0),
   fh1Xsec(0x0),
   fh1Trials(0x0),
+  fh1AvgTrials(0x0),
   fh1PtHard(0x0),
   fh1PtHardTrials(0x0),
   fh1SelectionInfoESD(0x0),
   fh1EventCutInfoESD(0),
   fh1CentralityESD(0),
   fh1Centrality(0),
-  fh1ReducedTrigger(0),
   fh1RP(0),
+  fh2ReducedTrigger(0),
   fh2CentralityTriggerESD(0),
   fh2CentralityTrigger(0),
   fh2TriggerCount(0x0),
@@ -250,6 +254,7 @@ Bool_t AliAnalysisTaskJetServices::Notify()
     // construct a poor man average trials 
     Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
     if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
+    fh1Trials->Fill("#sum{avg ntrials}",ftrials); 
   }  
   return kTRUE;
 }
@@ -279,14 +284,18 @@ void AliAnalysisTaskJetServices::UserCreateOutputObjects()
   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
   fHistList->Add(fh1Trials);
 
-  const Int_t nBinPt = 100;
+  fh1AvgTrials = new TH1F("fh1AvgTrials","avg trials root file",1,0,1);
+  fh1AvgTrials->GetXaxis()->SetBinLabel(1,"#sum{avg ntrials}");
+  fHistList->Add(fh1AvgTrials);
+
+  const Int_t nBinPt = 125;
   Double_t binLimitsPt[nBinPt+1];
   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
     if(iPt == 0){
       binLimitsPt[iPt] = 0.0;
     }
     else {// 1.0
-      binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.5;
+      binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.;
     }
   }
   
@@ -300,8 +309,8 @@ void AliAnalysisTaskJetServices::UserCreateOutputObjects()
   fh1RP = new TH1F("fh1RP","RP;#Psi",440, -1.*TMath::Pi(), 2.*TMath::Pi());
   fHistList->Add(fh1RP);
 
-  fh1ReducedTrigger = new TH1F("fh1ReducedTrigger","red trigger;red trigger",1<<fNTrigger,-0.5,(1<<fNTrigger)-0.5);
-  fHistList->Add(fh1ReducedTrigger);
+  fh2ReducedTrigger = new TH2F("fh2ReducedTrigger","red trigger;red trigger;centrality",1<<fNTrigger,-0.5,(1<<fNTrigger)-0.5,100,-0.5,99.5);
+  fHistList->Add(fh2ReducedTrigger);
 
   fh2CentralityTriggerESD = new TH2F("fh2CentralityTriggerESD",";cent;trigger no",103,-1,102,fNTrigger,-0.5,fNTrigger-0.5);
   fHistList->Add(fh2CentralityTriggerESD);
@@ -310,8 +319,8 @@ void AliAnalysisTaskJetServices::UserCreateOutputObjects()
   fHistList->Add(fh2CentralityTrigger);
 
   for(int i = 0;i<fNTrigger;++i){
-    fh2CentralityTriggerESD->GetYaxis()->SetBinLabel(i+2,fTriggerName[i].Data());
-    fh2CentralityTrigger->GetYaxis()->SetBinLabel(i+2,fTriggerName[i].Data());
+    fh2CentralityTriggerESD->GetYaxis()->SetBinLabel(i+1,fTriggerName[i].Data());
+    fh2CentralityTrigger->GetYaxis()->SetBinLabel(i+1,fTriggerName[i].Data());
   }
 
 
@@ -494,7 +503,7 @@ void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
 
 
   fEventCutInfoESD |= kPhysicsSelectionCut; // other alreay set via IsEventSelected
-  fh1EventCutInfoESD->Fill(fEventCutInfoESD);
+  fh1EventCutInfoESD->Fill (fEventCutInfoESD);
 
   if(esdEventSelected) fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kVertexIn;
   if(esdEventPileUp)   fSelectionInfoESD |=  AliAnalysisHelperJetTasks::kIsPileUp;
@@ -526,7 +535,7 @@ void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
 
   Float_t cent = 0;
   if(fCollisionType==kPbPb){
-    if(aod)cent = aod->GetHeader()->GetCentrality();
+    if(aod)cent = ((AliVAODHeader*)aod->GetHeader())->GetCentrality();
     if(fDebug)Printf("%s:%d %3.3f",(char*)__FILE__,__LINE__,cent);
     if(cent<0)cent = 101;
   }
@@ -538,8 +547,9 @@ void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
 
   if(esd){
     const AliESDVertex *vtxESD = esd->GetPrimaryVertex();
+    const AliESDVertex *vtxESDSPD = esd->GetPrimaryVertexSPD();
     esdVtxValid = IsVertexValid(vtxESD);
-    esdVtxIn = IsVertexIn(vtxESD);
+    esdVtxIn = IsVertexIn(vtxESD,vtxESDSPD);
     if(aodH&&physicsSelection&&fFilterAODCollisions&&aod){
       if(fDebug)Printf("%s:%d Centrality %3.3f vtxin %d",(char*)__FILE__,__LINE__,cent,esdVtxIn);
       if(cent<=80&&esdVtxIn){
@@ -599,7 +609,7 @@ void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
            ir |= (1<<it);
          }
        }
-       fh1ReducedTrigger->Fill(ir);
+       fh2ReducedTrigger->Fill(ir,cent);
       }
     }
   }
@@ -608,13 +618,14 @@ void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
 
   if(aod){
     const AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
+    const AliAODVertex *vtxAODSPD = aod->GetPrimaryVertexSPD();
     aodVtxValid = IsVertexValid(vtxAOD);
-    aodVtxIn = IsVertexIn(vtxAOD);
+    aodVtxIn = IsVertexIn(vtxAOD,vtxAODSPD);
     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());
+    Bool_t cand = ((AliVAODHeader*)aod->GetHeader())->GetOfflineTrigger()&fPhysicsSelectionFlag;
+    if(fDebug)Printf("%s:%d AOD selection %d %d",(char*)__FILE__,__LINE__,cand,((AliVAODHeader*)aod->GetHeader())->GetOfflineTrigger());
     fh2TriggerCount->Fill(0.,kAllTriggered); 
     fh2TriggerCount->Fill(iCl,kAllTriggered); 
     if(cand){
@@ -660,7 +671,7 @@ void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
       TList recTracks;
       GetListOfTracks(&recTracks);
       CalculateReactionPlaneAngleVZERO(aod);
-      fRPAngle = aod->GetHeader()->GetEventplane();
+      fRPAngle = ((AliVAODHeader*)aod->GetHeader())->GetEventplane();
       fh1RP->Fill(fRPAngle);
       fh2RPCentrality->Fill(fCentrality,fRPAngle);
       fh2RPACentrality->Fill(fCentrality,fPsiVZEROA);
@@ -687,7 +698,7 @@ void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
   Double_t ptHard = 0; 
   Double_t nTrials = 1; // Trials for MC trigger 
 
-  fh1Trials->Fill("#sum{ntrials}",fAvgTrials); 
+  fh1AvgTrials->Fill("#sum{avg ntrials}",fAvgTrials); 
   AliMCEvent* mcEvent = MCEvent();
   //    AliStack *pStack = 0; 
   if(mcEvent){
@@ -772,7 +783,8 @@ void AliAnalysisTaskJetServices::UserExec(Option_t */*option*/)
 Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliESDEvent* esd){
   if(!esd)return kFALSE;
   const AliESDVertex *vtx = esd->GetPrimaryVertex();
-  return IsVertexIn(vtx); // vertex in calls vertex valid
+  const AliESDVertex *vtxSPD = esd->GetPrimaryVertex();
+  return IsVertexIn(vtx,vtxSPD); // vertex in calls vertex valid
 }
 
 AliAnalysisTaskJetServices::~AliAnalysisTaskJetServices(){
@@ -828,7 +840,8 @@ void AliAnalysisTaskJetServices::SetV0Centroids(TProfile *xa,TProfile *ya,TProfi
 Bool_t AliAnalysisTaskJetServices::IsEventSelected(const AliAODEvent* aod) const {
   if(!aod)return kFALSE;
   const AliAODVertex *vtx = aod->GetPrimaryVertex();
-  return IsVertexIn(vtx); // VertexIn calls VertexValid
+  const AliAODVertex *vtxSPD = aod->GetPrimaryVertexSPD();
+  return IsVertexIn(vtx,vtxSPD); // VertexIn calls VertexValid
 }
 
 Bool_t  AliAnalysisTaskJetServices::IsVertexValid ( const AliESDVertex* vtx) {
@@ -894,18 +907,19 @@ Bool_t  AliAnalysisTaskJetServices::IsVertexValid ( const AliAODVertex* vtx) con
     return kFALSE;
   }
 
-
+  TString vtxName(vtx->GetName());
   if(fDebug){
     Printf(" n contrib %d",vtx->GetNContributors());
+    Printf("vtxName: %s",vtxName.Data());
     vtx->Print();
   }
   
   //  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...
+  // no dispersion yet... 
   /* 
   TString vtxTitle(vtx->GetTitle());
   if(vtxTitle.Contains("vertexer: Z")){
@@ -916,13 +930,15 @@ Bool_t  AliAnalysisTaskJetServices::IsVertexValid ( const AliAODVertex* vtx) con
 }
 
 
-Bool_t  AliAnalysisTaskJetServices::IsVertexIn (const AliESDVertex* vtx) {
+Bool_t  AliAnalysisTaskJetServices::IsVertexIn (const AliESDVertex* vtx,const AliESDVertex* vtxSPD) {
 
   if(!IsVertexValid(vtx))return kFALSE;
 
   Float_t zvtx = vtx->GetZ();
   Float_t yvtx = vtx->GetY();
   Float_t xvtx = vtx->GetX();
+  Float_t deltaZ = zvtx - vtxSPD->GetZ();
+  if(fDebug)Printf("%s:%d deltaz %3.3f ",__FILE__,__LINE__,deltaZ);
 
   xvtx -= fVtxXMean;
   yvtx -= fVtxYMean;
@@ -939,17 +955,25 @@ Bool_t  AliAnalysisTaskJetServices::IsVertexIn (const AliESDVertex* vtx) {
     return kFALSE;
   }
   fEventCutInfoESD |= kVertexRCut;  
+
+
+
+
+  if(TMath::Abs(deltaZ)>fVtxDeltaZCut)return kFALSE;
+  fEventCutInfoESD |= kVertexDeltaZCut;  
   return kTRUE;
 }
 
 
-Bool_t  AliAnalysisTaskJetServices::IsVertexIn (const AliAODVertex* vtx) const {
+Bool_t  AliAnalysisTaskJetServices::IsVertexIn (const AliAODVertex* vtx,const AliAODVertex* vtxSPD) const {
 
   if(!IsVertexValid(vtx))return kFALSE;
 
   Float_t zvtx = vtx->GetZ();
   Float_t yvtx = vtx->GetY();
   Float_t xvtx = vtx->GetX();
+  Float_t deltaZ = zvtx - vtxSPD->GetZ();
+  if(fDebug)Printf("%s:%d deltaz %3.3f ",__FILE__,__LINE__,deltaZ);
 
   xvtx -= fVtxXMean;
   yvtx -= fVtxYMean;
@@ -958,6 +982,12 @@ Bool_t  AliAnalysisTaskJetServices::IsVertexIn (const AliAODVertex* vtx) const {
   Float_t r2   = yvtx*yvtx+xvtx*xvtx;      
 
   Bool_t vertexIn = TMath::Abs(zvtx)<fVtxZCut&&r2<(fVtxRCut*fVtxRCut);
+
+
+  if(TMath::Abs(deltaZ)>fVtxDeltaZCut)vertexIn = kFALSE;
+  if(fDebug)Printf("%s:%d vertexIn %d ",__FILE__,__LINE__,vertexIn);
+  
+
   return vertexIn;
 }
 
@@ -1039,7 +1069,7 @@ Int_t AliAnalysisTaskJetServices::GetEventClass(AliESDEvent *esd){
 Int_t AliAnalysisTaskJetServices::GetEventClass(AliAODEvent *aod){
 
   if(fCollisionType==kPbPb){
-    Float_t cent = aod->GetHeader()->GetCentrality();
+    Float_t cent = ((AliVAODHeader*)aod->GetHeader())->GetCentrality();
     if(cent>80||cent<0)return 5;
     if(cent>50)return 4;
     if(cent>30)return 3;
@@ -1196,7 +1226,8 @@ Int_t  AliAnalysisTaskJetServices::GetListOfTracks(TList *list){
     return iCount;
   }
   for(int it = 0;it < aod->GetNumberOfTracks();++it){
-    AliAODTrack *tr = aod->GetTrack(it);
+    AliAODTrack *tr = dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
+    if(!tr) AliFatal("Not a standard AOD");
     if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
     if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue;
     if(tr->Pt()<fMinTrackPt)continue;