]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG3/vertexingHF/AliAnalysisTaskSECompareHF.cxx
Merged tasks DStar and DStarSpectra (Alessandro, Yifei)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSECompareHF.cxx
index 5b061ae87c1673aa91a10e2496d7a472e26d7198..cbfdba42500f758e73bc8f01ca824fc221a41d64 100644 (file)
 #include <TList.h>
 #include <TH1F.h>
 
+#include "AliAnalysisManager.h"
+#include "AliAODHandler.h"
 #include "AliAODEvent.h"
 #include "AliAODVertex.h"
+#include "AliESDtrack.h"
 #include "AliAODTrack.h"
 #include "AliAODMCHeader.h"
 #include "AliAODMCParticle.h"
 #include "AliAODRecoDecayHF2Prong.h"
+#include "AliAODRecoDecayHF3Prong.h"
+#include "AliAODRecoDecayHF4Prong.h"
 #include "AliAODRecoCascadeHF.h"
 #include "AliAnalysisVertexingHF.h"
 #include "AliAnalysisTaskSE.h"
@@ -121,7 +126,8 @@ void AliAnalysisTaskSECompareHF::UserCreateOutputObjects()
   fHistNEvents->SetMinimum(0);
   fOutput->Add(fHistNEvents);
 
-  fNtupleCmp = new TNtuple("fNtupleCmp","Charm comparison","pdg:nprongs:VxRec:VxTrue:ErrVx:VyRec:VyTrue:ErrVy:VzRec:VzTrue:ErrVz:Chi2toNDF:PtRec:Mrec");
+  OpenFile(2); // 2 is the slot number of the ntuple
+  fNtupleCmp = new TNtuple("fNtupleCmp","Charm comparison","pdg:nprongs:VxRec:VxTrue:ErrVx:VyRec:VyTrue:ErrVy:VzRec:VzTrue:ErrVz:Chi2toNDF:PtRec:Mrec:CPta:Prodd0");
 
   return;
 }
@@ -134,37 +140,57 @@ void AliAnalysisTaskSECompareHF::UserExec(Option_t */*option*/)
 
   
   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
-  
-  fHistNEvents->Fill(0); // count event
-  // Post the data already here
-  PostData(1,fOutput);
 
-  // load HF vertices                     
-  TClonesArray *inputArrayVertices =
-    (TClonesArray*)aod->GetList()->FindObject("VerticesHF");
+  TClonesArray *inputArrayVertices = 0;
+  TClonesArray *inputArrayD0toKpi = 0;
+  TClonesArray *inputArrayDstar = 0;
+
+  if(!aod && AODEvent() && IsStandardAOD()) {
+    // In case there is an AOD handler writing a standard AOD, use the AOD 
+    // event in memory rather than the input (ESD) event.    
+    aod = dynamic_cast<AliAODEvent*> (AODEvent());
+    // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
+    // have to taken from the AOD event hold by the AliAODExtension
+    AliAODHandler* aodHandler = (AliAODHandler*) 
+      ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
+    if(aodHandler->GetExtensions()) {
+      AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
+      AliAODEvent *aodFromExt = ext->GetAOD();
+      // load HF vertices                
+      inputArrayVertices = (TClonesArray*)aodFromExt->GetList()->FindObject("VerticesHF");
+      // load D0->Kpi candidates
+      inputArrayD0toKpi = (TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
+      // load D*+ candidates                                                   
+      inputArrayDstar = (TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
+    }
+  } else {
+    // load HF vertices                
+    inputArrayVertices = (TClonesArray*)aod->GetList()->FindObject("VerticesHF");
+    // load D0->Kpi candidates                                                 
+    inputArrayD0toKpi = (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
+    // load D*+ candidates                                                   
+    inputArrayDstar = (TClonesArray*)aod->GetList()->FindObject("Dstar");
+  }
+
+
   if(!inputArrayVertices) {
     printf("AliAnalysisTaskSECompareHF::UserExec: Vertices branch not found!\n");
     return;
   }
-
-  // load D0->Kpi candidates                                                   
-  TClonesArray *inputArrayD0toKpi =
-    (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
   if(!inputArrayD0toKpi) {
     printf("AliAnalysisTaskSECompareHF::UserExec: D0toKpi branch not found!\n");
     return;
   }
-
-  // load D*+ candidates                                                   
-  TClonesArray *inputArrayDstar =
-    (TClonesArray*)aod->GetList()->FindObject("Dstar");
   if(!inputArrayDstar) {
     printf("AliAnalysisTaskSECompareHF::UserExec: Dstar branch not found!\n");
     return;
   }
   
 
+  fHistNEvents->Fill(0); // count event
+  // Post the data already here
+  PostData(1,fOutput);
+
   // AOD primary vertex
   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
   //vtx1->Print();
@@ -192,9 +218,13 @@ void AliAnalysisTaskSECompareHF::UserExec(Option_t */*option*/)
   AliAODRecoDecayHF3Prong *d3=0;
   AliAODRecoDecayHF4Prong *d4=0;
 
+  Int_t pdgDgD0toKpi[2]={321,211};
+  Int_t pdgDgDplustoKpipi[3]={321,211,211};
+  Int_t pdgDgD0toKpipipi[4]={321,211,211,211};
+
   // loop over vertices
   Int_t nVertices = inputArrayVertices->GetEntriesFast();
-  printf("Number of vertices: %d\n",nVertices);
+  if(fDebug>1) printf("Number of vertices: %d\n",nVertices);
 
   for (Int_t iVtx = 0; iVtx < nVertices; iVtx++) {
     AliAODVertex *vtx = (AliAODVertex*)inputArrayVertices->UncheckedAt(iVtx);
@@ -209,13 +239,22 @@ void AliAnalysisTaskSECompareHF::UserExec(Option_t */*option*/)
 
     // get parent AliAODRecoDecayHF
     nprongs= vtx->GetNDaughters();
+    // check that the daughters have kITSrefit and kTPCrefit
+    Bool_t allDgOK=kTRUE;
+    for(Int_t idg=0; idg<nprongs; idg++) {
+      AliAODTrack *track = (AliAODTrack*)vtx->GetDaughter(idg);
+      if(!(track->GetStatus()&AliESDtrack::kITSrefit)) allDgOK=kFALSE;
+      if(!(track->GetStatus()&AliESDtrack::kTPCrefit)) allDgOK=kFALSE;
+    }
+    if(!allDgOK) continue;
+
 
     switch(nprongs) {
     case 2: // look for D0->Kpi
       d2 = (AliAODRecoDecayHF2Prong*)vtx->GetParent();
       if(d2->IsLikeSign()) continue;
       if(d2->Charge() != 0) continue; // these are D* 
-      lab = d2->MatchToMC(421,mcArray);
+      lab = d2->MatchToMC(421,mcArray,2,pdgDgD0toKpi);
       if(lab>=0) {
        unsetvtx=kFALSE;
        if(!d2->GetOwnPrimaryVtx()) {
@@ -223,30 +262,32 @@ void AliAnalysisTaskSECompareHF::UserExec(Option_t */*option*/)
          unsetvtx=kTRUE;
        }
        okD0=0; okD0bar=0; 
-       //if(d2->SelectD0(fVHF->GetD0toKpiCuts(),okD0,okD0bar)) { // no cuts, otherwise we bias the resolution
-       AliAODMCParticle *dMC = (AliAODMCParticle*)mcArray->At(lab);
-       pdg = dMC->GetPdgCode();
-       invmass = (pdg==421 ? d2->InvMassD0() : d2->InvMassD0bar());
-       fHistMass->Fill(invmass);
-       // Post the data already here
-       PostData(1,fOutput);
-       // get a daughter for true pos of decay vertex
-       AliAODMCParticle *dg0MC = (AliAODMCParticle*)mcArray->At(dMC->GetDaughter(0));
-       dg0MC->XvYvZv(posTrue);
-       fNtupleCmp->Fill(pdg,nprongs,
-                        posRec[0],posTrue[0],errx,
-                        posRec[1],posTrue[1],erry,
-                        posRec[2],posTrue[2],errz,
-                        d2->GetReducedChi2(),d2->Pt(),invmass);
-       //PostData(2,fNtupleCmp);
-       //}
+       if(d2->SelectD0(fVHF->GetD0toKpiCuts(),okD0,okD0bar)) { // beware: cuts may bias the resolution!
+         AliAODMCParticle *dMC = (AliAODMCParticle*)mcArray->At(lab);
+         pdg = dMC->GetPdgCode();
+         invmass = (pdg==421 ? d2->InvMassD0() : d2->InvMassD0bar());
+         // get a daughter for true pos of decay vertex
+         AliAODMCParticle *dg0MC = (AliAODMCParticle*)mcArray->At(dMC->GetDaughter(0));
+         dg0MC->XvYvZv(posTrue);
+         fHistMass->Fill(invmass);
+         // Post the data already here
+         PostData(1,fOutput);
+         Float_t tmp[16]={(Float_t)pdg,(Float_t)nprongs,
+                          (Float_t)posRec[0],(Float_t)posTrue[0],(Float_t)errx,
+                          (Float_t)posRec[1],(Float_t)posTrue[1],(Float_t)erry,
+                          (Float_t)posRec[2],(Float_t)posTrue[2],(Float_t)errz,
+                          (Float_t)d2->GetReducedChi2(),(Float_t)d2->Pt(),(Float_t)invmass,
+                          (Float_t)d2->CosPointingAngle(),(Float_t)d2->Prodd0d0()};
+         fNtupleCmp->Fill(tmp);
+         PostData(2,fNtupleCmp);
+       }
        if(unsetvtx) d2->UnsetOwnPrimaryVtx();
       }
       break;
     case 3: // look for D+
       d3 = (AliAODRecoDecayHF3Prong*)vtx->GetParent();
       if(d3->IsLikeSign()) continue;
-      lab = d3->MatchToMC(411,mcArray);
+      lab = d3->MatchToMC(411,mcArray,3,pdgDgDplustoKpipi);
       if(lab>=0) {
        unsetvtx=kFALSE;
        if(!d3->GetOwnPrimaryVtx()) {
@@ -260,12 +301,14 @@ void AliAnalysisTaskSECompareHF::UserExec(Option_t */*option*/)
        // get a daughter for true pos of decay vertex
        AliAODMCParticle *dg0MC = (AliAODMCParticle*)mcArray->At(dMC->GetDaughter(0));
        dg0MC->XvYvZv(posTrue);
-       fNtupleCmp->Fill(pdg,nprongs,
-                        posRec[0],posTrue[0],errx,
-                        posRec[1],posTrue[1],erry,
-                        posRec[2],posTrue[2],errz,
-                        d3->GetReducedChi2(),d3->Pt(),invmass);
-       //PostData(2,fNtupleCmp);
+       Float_t tmp[16]={(Float_t)pdg,(Float_t)nprongs,
+                        (Float_t)posRec[0],(Float_t)posTrue[0],(Float_t)errx,
+                        (Float_t)posRec[1],(Float_t)posTrue[1],(Float_t)erry,
+                        (Float_t)posRec[2],(Float_t)posTrue[2],(Float_t)errz,
+                        (Float_t)d3->GetReducedChi2(),(Float_t)d3->Pt(),(Float_t)invmass,
+                        (Float_t)d3->CosPointingAngle(),(Float_t)(d3->Getd0Prong(0)*d3->Getd0Prong(1)*d3->Getd0Prong(2))};
+       fNtupleCmp->Fill(tmp);
+       PostData(2,fNtupleCmp);
        //}
        if(unsetvtx) d3->UnsetOwnPrimaryVtx();
       }
@@ -273,7 +316,7 @@ void AliAnalysisTaskSECompareHF::UserExec(Option_t */*option*/)
     case 4: // look for D0->Kpipipi
       d4 = (AliAODRecoDecayHF4Prong*)vtx->GetParent();
       if(d4->IsLikeSign()) continue;
-      lab = d4->MatchToMC(421,mcArray);
+      lab = d4->MatchToMC(421,mcArray,4,pdgDgD0toKpipipi);
       if(lab>=0) {
        unsetvtx=kFALSE;
        if(!d4->GetOwnPrimaryVtx()) {
@@ -289,12 +332,14 @@ void AliAnalysisTaskSECompareHF::UserExec(Option_t */*option*/)
        // get a daughter for true pos of decay vertex
        AliAODMCParticle *dg0MC = (AliAODMCParticle*)mcArray->At(dMC->GetDaughter(0));
        dg0MC->XvYvZv(posTrue);
-       fNtupleCmp->Fill(pdg,nprongs,
-                        posRec[0],posTrue[0],errx,
-                        posRec[1],posTrue[1],erry,
-                        posRec[2],posTrue[2],errz,
-                        d4->GetReducedChi2(),d4->Pt(),invmass);
-       //PostData(2,fNtupleCmp);
+       Float_t tmp[16]={(Float_t)pdg,(Float_t)nprongs,
+                        (Float_t)posRec[0],(Float_t)posTrue[0],(Float_t)errx,
+                        (Float_t)posRec[1],(Float_t)posTrue[1],(Float_t)erry,
+                        (Float_t)posRec[2],(Float_t)posTrue[2],(Float_t)errz,
+                        (Float_t)d4->GetReducedChi2(),(Float_t)d4->Pt(),(Float_t)invmass,
+                        (Float_t)d4->CosPointingAngle(),(Float_t)(d4->Getd0Prong(0)*d4->Getd0Prong(1)*d4->Getd0Prong(2)*d4->Getd0Prong(3))};
+       fNtupleCmp->Fill(tmp);
+       PostData(2,fNtupleCmp);
        //}
        if(unsetvtx) d4->UnsetOwnPrimaryVtx();
       }