]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG3/vertexingHF/AliAnalysisTaskSECharmFraction.cxx
move macros to macros/
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSECharmFraction.cxx
index 5bfb588f49bb37087975a8b1882314743325087f..7ad7bd96cc095104ec5cb45be1a35e977ea6b2fb 100644 (file)
@@ -28,6 +28,9 @@
 #include <TDatabasePDG.h>
 #include <TMath.h>
 #include <TROOT.h>
+
+#include "AliAnalysisManager.h"
+#include "AliAODHandler.h"
 #include "AliAODEvent.h"
 #include "AliAODRecoDecayHF2Prong.h"
 #include "AliAODRecoDecayHF.h"
@@ -55,6 +58,7 @@ ClassImp(AliAnalysisTaskSECharmFraction)
     : AliAnalysisTaskSE(),
       fVHFloose(0),
       fVHFtight(0),
+      fReadMC(kFALSE),
       fmD0PDG(),
       fnbins(),
       fptbins(0),
@@ -92,6 +96,7 @@ ClassImp(AliAnalysisTaskSECharmFraction)
     : AliAnalysisTaskSE(name),
       fVHFloose(0),
       fVHFtight(0),
+      fReadMC(kFALSE),
       fmD0PDG(),
       fnbins(),
       fptbins(0),
@@ -155,6 +160,7 @@ AliAnalysisTaskSECharmFraction::AliAnalysisTaskSECharmFraction(const char *name,
   : AliAnalysisTaskSE(name),
     fVHFloose(0),
     fVHFtight(0),
+    fReadMC(kFALSE),
     fmD0PDG(),
     fnbins(),
     fptbins(0),
@@ -2538,13 +2544,25 @@ void AliAnalysisTaskSECharmFraction::UserExec(Option_t */*option*/)
   
   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 
-  // In case there is an AOD handler writing a standard AOD, use the AOD 
-  // event in memory rather than the input (ESD) event.
-  if (!aod && AODEvent() && IsStandardAOD()) aod = dynamic_cast<AliAODEvent*> (AODEvent());
+  TClonesArray *arrayD0toKpi=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();
+      arrayD0toKpi=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
+    }
+  } else {
+    arrayD0toKpi=(TClonesArray*)aod->GetList()->FindObject("D0toKpi");
+  }
 
-  // load D0->Kpi candidates                                                   
-  TClonesArray *arrayD0toKpi =
-    (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
   if(!arrayD0toKpi) {
     Printf("AliAnalysisTaskSECharmFraction::UserExec: D0toKpi branch not found!\n");
     return;
@@ -2552,36 +2570,39 @@ void AliAnalysisTaskSECharmFraction::UserExec(Option_t */*option*/)
   
   // AOD primary vertex
   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
-  
-  
-  // load MC particles
-  TClonesArray *arrayMC = 
-    (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
-  if(!arrayMC) {
-    Printf("AliAnalysisTaskSECharmFraction::UserExec: MC particles branch not found!\n");
-    return;
-  }
-  
-  // load MC header
-  AliAODMCHeader *aodmcHeader = 
-    (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
-  if(!aodmcHeader) {
-    Printf("AliAnalysisTaskSECharmFraction::UserExec: MC header branch not found!\n");
-    return;
-  }
-  
-  // MC primary vertex
+  TClonesArray *arrayMC=0;
+  AliAODMCHeader *aodmcHeader=0;
   Double_t vtxTrue[3];
-  aodmcHeader->GetVertex(vtxTrue);
   
+  if(fReadMC){
+    // load MC particles
+    arrayMC = 
+      (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
+    if(!arrayMC) {
+      Printf("AliAnalysisTaskSECharmFraction::UserExec: MC particles branch not found!\n");
+      return;
+    }
+    
+    // load MC header
+    aodmcHeader = 
+      (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
+    if(!aodmcHeader) {
+      Printf("AliAnalysisTaskSECharmFraction::UserExec: MC header branch not found!\n");
+      return;
+    }
+    
+    // MC primary vertex
+    aodmcHeader->GetVertex(vtxTrue);
+  }    
   if (!aod) {
     Printf("ERROR: aod not available");
     return;
   }
+  
   //histogram filled with 1 for every AOD
   fNentries->Fill(1);
   PostData(1,fNentries);
-
+  
 
   //Printf("There are %d tracks in this event", aod->GetNumberOfTracks());
   //  Int_t nTotHF=0,nTotDstar=0,nTot3Prong=0;
@@ -2595,7 +2616,7 @@ void AliAnalysisTaskSECharmFraction::UserExec(Option_t */*option*/)
   Double_t invMassD0,invMassD0bar,ptD0,massmumtrue;
   
  
-  AliAODRecoDecayHF *aodDMC=0x0;// to be used to create a fake true sec vertex
+  AliAODRecoDecayHF *aodDMC=0;// to be used to create a fake true sec vertex
   // make trkIDtoEntry register (temporary)
   Int_t trkIDtoEntry[100000];
   for(Int_t it=0;it<aod->GetNumberOfTracks();it++) {
@@ -2614,7 +2635,7 @@ void AliAnalysisTaskSECharmFraction::UserExec(Option_t */*option*/)
   //   cout<<"Number of D0->Kpi: "<<nD0toKpi<<endl;
   
   for (Int_t iD0toKpi = 0; iD0toKpi < nD0toKpi; iD0toKpi++) {
-    if(aodDMC!=0x0)delete aodDMC;
+    if(aodDMC)delete aodDMC;
       
     isPeakD0=kFALSE;
     isPeakD0bar=kFALSE;
@@ -2649,7 +2670,10 @@ void AliAnalysisTaskSECharmFraction::UserExec(Option_t */*option*/)
     if((isSideBandD0||isSideBandD0bar)&&!(isPeakD0||isPeakD0bar))isSideBand=kTRUE;// TEMPORARY, NOT DONE IN THE METHOD CALLED ABOVE ONLY FOR FURTHER SIDE BAND STUDY
   
     // INVESTIGATE SIGNAL TYPE : ACCESS TO MC INFORMATION
-    aodDMC=GetD0toKPiSignalType(d,arrayMC,signallevel,massmumtrue,vtxTrue);
+    if(fReadMC){
+      aodDMC=GetD0toKPiSignalType(d,arrayMC,signallevel,massmumtrue,vtxTrue);
+    }
+    else signallevel=0;
     if(!isinacceptance)signallevel=9;
     fSignalType->Fill(signallevel);
   
@@ -2676,30 +2700,30 @@ void AliAnalysisTaskSECharmFraction::UserExec(Option_t */*option*/)
     //            CANDIDATE VARIABLES   
 
     //NO CUTS Case: force okD0 and okD0bar = kTRUE
-    if(signallevel==1)FillHistos(d,flistNoCutsSignal,ptbin,kTRUE,kTRUE,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
+    if(signallevel==1||signallevel==0)FillHistos(d,flistNoCutsSignal,ptbin,kTRUE,kTRUE,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
     else if(signallevel==2)FillHistos(d,flistNoCutsFromDstar,ptbin,kTRUE,kTRUE,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
     else if(signallevel==3||signallevel==4)FillHistos(d,flistNoCutsFromB,ptbin,kTRUE,kTRUE,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
     else if(signallevel==-1||signallevel==7||signallevel==8||signallevel==10||signallevel==9)FillHistos(d,flistNoCutsBack,ptbin,kTRUE,kTRUE,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
     else if(signallevel==5||signallevel==6)FillHistos(d,flistNoCutsOther,ptbin,kTRUE,kTRUE,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
 
     //LOOSE CUTS Case
-    if(signallevel==1)FillHistos(d,flistLsCutsSignal,ptbin,okd0loose,okd0barloose,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
+    if(signallevel==1||signallevel==0)FillHistos(d,flistLsCutsSignal,ptbin,okd0loose,okd0barloose,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
     else if(signallevel==2)FillHistos(d,flistLsCutsFromDstar,ptbin,okd0loose,okd0barloose,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
     else if(signallevel==3||signallevel==4)FillHistos(d,flistLsCutsFromB,ptbin,okd0loose,okd0barloose,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
     else if(signallevel==-1||signallevel==7||signallevel==8||signallevel==10)FillHistos(d,flistLsCutsBack,ptbin,okd0loose,okd0barloose,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
     else if(signallevel==5||signallevel==6)FillHistos(d,flistLsCutsOther,ptbin,okd0loose,okd0barloose,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
 
     //TIGHT CUTS Case
-    if(signallevel==1)FillHistos(d,flistTghCutsSignal,ptbin,okd0tight,okd0bartight,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
+    if(signallevel==1||signallevel==0)FillHistos(d,flistTghCutsSignal,ptbin,okd0tight,okd0bartight,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
     else if(signallevel==2)FillHistos(d,flistTghCutsFromDstar,ptbin,okd0tight,okd0bartight,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
     else if(signallevel==3||signallevel==4)FillHistos(d,flistTghCutsFromB,ptbin,okd0tight,okd0bartight,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
     else if(signallevel==-1||signallevel==7||signallevel==8||signallevel==10)FillHistos(d,flistTghCutsBack,ptbin,okd0tight,okd0bartight,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
     else if(signallevel==5||signallevel==6)FillHistos(d,flistTghCutsOther,ptbin,okd0tight,okd0bartight,invMassD0,invMassD0bar,isPeakD0,isPeakD0bar,isSideBand,massmumtrue,aodDMC,vtxTrue);
     
     
-    if(aodDMC!=0x0){
+    if(aodDMC){
       delete aodDMC;
-      aodDMC=0x0;
+      aodDMC=0;
     }
     
     if(unsetvtx) d->UnsetOwnPrimaryVtx();
@@ -2822,9 +2846,9 @@ AliAODRecoDecayHF* AliAnalysisTaskSECharmFraction::GetD0toKPiSignalType(const Al
                         //                            11: end of the method without output
                         //                            12: different result than MatchToMC method
 
-  AliAODMCParticle *mum1=0x0;
-  AliAODMCParticle *b1=0x0,*b2=0x0;
-  AliAODMCParticle *grandmoth1=0x0;
+  AliAODMCParticle *mum1=0;
+  AliAODMCParticle *b1=0,*b2=0;
+  AliAODMCParticle *grandmoth1=0;
   massMumTrue=-1;
   
   Int_t pdgmum,dglabels[2],matchtoMC;
@@ -2832,8 +2856,8 @@ AliAODRecoDecayHF* AliAnalysisTaskSECharmFraction::GetD0toKPiSignalType(const Al
   // get daughter AOD tracks
   AliAODTrack *trk0 = (AliAODTrack*)d->GetDaughter(0);
   AliAODTrack *trk1 = (AliAODTrack*)d->GetDaughter(1);
-  AliAODRecoDecayHF *aodDMC=0x0;
-  if(trk0==0x0||trk1==0x0){
+  AliAODRecoDecayHF *aodDMC=0;
+  if(!trk0 || !trk1){
     AliDebug(2,"Delete tracks I AM \n");
   
     signaltype=-1;
@@ -2842,7 +2866,7 @@ AliAODRecoDecayHF* AliAnalysisTaskSECharmFraction::GetD0toKPiSignalType(const Al
   }
   dglabels[0]=trk0->GetLabel();
   dglabels[1]=trk1->GetLabel();
-  if(dglabels[0]==-1||dglabels[1]==-1){
+  if(dglabels[0]<0 || dglabels[1]<0){
     AliDebug(2,"HERE I AM \n");
 
     //fake tracks
@@ -2856,9 +2880,15 @@ AliAODRecoDecayHF* AliAnalysisTaskSECharmFraction::GetD0toKPiSignalType(const Al
   b1=(AliAODMCParticle*)arrayMC->At(trk0->GetLabel());
   b2=(AliAODMCParticle*)arrayMC->At(trk1->GetLabel());
   
-  if(b1->GetMother()==-1||b2->GetMother()==-1){
+  if(!b1 || !b2) {
+    //Tracks with no mother  ????? FAKE DECAY VERTEX
+    
+    signaltype=10;
+    return aodDMC;
+  }
+
+  if(b1->GetMother()<0||b2->GetMother()<0){
     //Tracks with no mother  ????? FAKE DECAY VERTEX
     
     signaltype=10;
     return aodDMC;
@@ -2880,7 +2910,7 @@ AliAODRecoDecayHF* AliAnalysisTaskSECharmFraction::GetD0toKPiSignalType(const Al
   matchtoMC=d->MatchToMC(421,arrayMC,2,pdgdaughters);
   aodDMC=ConstructFakeTrueSecVtx(b1,b2,mum1,primaryVtx);
  
-  if(aodDMC==0x0){
+  if(aodDMC){
     signaltype=10;
     return aodDMC;
   }
@@ -2915,7 +2945,7 @@ AliAODRecoDecayHF* AliAnalysisTaskSECharmFraction::GetD0toKPiSignalType(const Al
    
   }
   
-  if(mum1->GetMother()==-1){
+  if(mum1->GetMother()<0){
     // A particle coming from nothing
     signaltype=10;
     return aodDMC;
@@ -2941,7 +2971,7 @@ AliAODRecoDecayHF* AliAnalysisTaskSECharmFraction::GetD0toKPiSignalType(const Al
   */
   //     if(fcheckMCD0){//THIS CHECK IS NEEDED TO AVOID POSSIBLE FAILURE IN THE SECOND WHILE, FOR DEBUGGING
   while(TMath::Abs(grandmoth1->GetPdgCode())!=4&&TMath::Abs(grandmoth1->GetPdgCode())!=5){
-    if(grandmoth1->GetMother()==-1){
+    if(grandmoth1->GetMother()<0){
       //### THE FOLLOWING IN CASE OF DEBUGGING ##########à
       /*printf("mother=-1, pdgcode: %d \n",grandmoth1->GetPdgCode());
        Int_t son=grandmoth1->GetDaughter(0);
@@ -2993,15 +3023,15 @@ AliAODRecoDecayHF* AliAnalysisTaskSECharmFraction::GetD0toKPiSignalType(const Al
 AliAODRecoDecayHF* AliAnalysisTaskSECharmFraction::ConstructFakeTrueSecVtx(const AliAODMCParticle *b1, const AliAODMCParticle *b2, const AliAODMCParticle *mum,Double_t *primaryVtxTrue){
   // CONSTRUCT A FAKE TRUE SECONDARY VERTEX (aodDMC)  
   //!!!NOTE THAT ONLY ONE MOTHER IS CONSIDERED: THE METHOD REQUIRES THE DAUGHTERS COME FROM THE SAME MOTHER !!
-  if(b1==0x0||b2==0x0)return 0x0;
-  if(mum==0x0)return 0x0;
+  if(!b1 || !b2)return 0;
+  if(!mum)return 0;
   Double_t pD[3],xD[3],pXtrTrue[2],pYtrTrue[2],pZtrTrue[2],xtr1[3],xtr2[3];
   Int_t charge[2]={0,0};
   if(b1->Charge()==-1)charge[0]=1;
   else {
     if(b2->Charge()==-1){
       //printf("Same charges for prongs \n");
-      return 0x0;
+      return 0;
     }
     charge[1]=1;
   }
@@ -3010,7 +3040,7 @@ AliAODRecoDecayHF* AliAnalysisTaskSECharmFraction::ConstructFakeTrueSecVtx(const
   pYtrTrue[charge[0]]=b1->Py();
   pZtrTrue[charge[0]]=b1->Pz();
   if(!b1->XvYvZv(xtr1)){
-    return 0x0;
+    return 0;
   }
   
   pXtrTrue[charge[1]]=b2->Px();
@@ -3019,11 +3049,11 @@ AliAODRecoDecayHF* AliAnalysisTaskSECharmFraction::ConstructFakeTrueSecVtx(const
   
   if(!mum->PxPyPz(pD)){
     //printf("!D from B:Get momentum failed \n");
-    return 0x0;
+    return 0;
   }
   if(!mum->XvYvZv(xD)){
     //printf("!D from B:Get position failed \n");
-    return 0x0;
+    return 0;
   }
   /* ############ THIS HAPPENS FROM TIME TO TIME: NUMERIC PROBLEM KNOWN #################
      if(pXtrTrue[0]+pXtrTrue[1]!=pD[0]){
@@ -3031,7 +3061,7 @@ AliAODRecoDecayHF* AliAnalysisTaskSECharmFraction::ConstructFakeTrueSecVtx(const
   
   
   if(!b2->XvYvZv(xtr2)){
-    return 0x0;
+    return 0;
   }
   Double_t d0dummy[2]={0.,0.};//TEMPORARY : dummy d0 for AliAODRecoDeay constructor
   AliAODRecoDecayHF* aodDMC=new AliAODRecoDecayHF(primaryVtxTrue,xD,2,0,pXtrTrue,pYtrTrue,pZtrTrue,d0dummy);
@@ -3117,26 +3147,26 @@ Bool_t AliAnalysisTaskSECharmFraction::FillHistos(AliAODRecoDecayHF2Prong *d,TLi
     if(okD0)((TH1F*)list->FindObject(str.Data()))->Fill(invMassD0);
     if(okD0bar)((TH1F*)list->FindObject(str.Data()))->Fill(invMassD0bar);
   }
-  
-  if(massmumtrue>0.){
-    str="hMassTrue";
-    str.Append(namehist.Data());
-    ((TH1F*)list->FindObject(str.Data()))->Fill(massmumtrue);
-    
-    if(isPeakD0||isPeakD0bar){
+  if(fReadMC){
+    if(massmumtrue>0.){
       str="hMassTrue";
       str.Append(namehist.Data());
-      str.Append("PM");
-      ((TH1F*)list->FindObject(str.Data()))->Fill(massmumtrue);
-    }
-    if(isSideBand){
-      str="hMassTrue";
-      str.Append(namehist.Data());
-      str.Append("SB");
       ((TH1F*)list->FindObject(str.Data()))->Fill(massmumtrue);
+      
+      if(isPeakD0||isPeakD0bar){
+       str="hMassTrue";
+       str.Append(namehist.Data());
+       str.Append("PM");
+       ((TH1F*)list->FindObject(str.Data()))->Fill(massmumtrue);
+      }
+      if(isSideBand){
+       str="hMassTrue";
+       str.Append(namehist.Data());
+       str.Append("SB");
+       ((TH1F*)list->FindObject(str.Data()))->Fill(massmumtrue);
+      }
     }
   }
-  
   // ################ D0 Impact Parameter Histos #####################
   if((isPeakD0&&okD0)||(isPeakD0bar&&okD0bar)){
     str="hd0D0";
@@ -3151,7 +3181,7 @@ Bool_t AliAnalysisTaskSECharmFraction::FillHistos(AliAODRecoDecayHF2Prong *d,TLi
     ((TH1F*)list->FindObject(str.Data()))->Fill(d->ImpParXY()*10000.);
      
     
-    if(vtxTrue){
+    if(fReadMC && vtxTrue){
       str="hd0D0VtxTrue";
       str.Append(namehist.Data());
       str.Append("PM");
@@ -3164,7 +3194,7 @@ Bool_t AliAnalysisTaskSECharmFraction::FillHistos(AliAODRecoDecayHF2Prong *d,TLi
       ((TH1F*)list->FindObject(str.Data()))->Fill(d->AliAODRecoDecay::ImpParXY(vtxTrue)*10000.);
     }
     
-    if(aodDMC!=0x0){
+    if(fReadMC && aodDMC){
       aodDMC->Print("");
       aodDMC->ImpParXY();
       aodDMC->Print("");
@@ -3194,7 +3224,7 @@ Bool_t AliAnalysisTaskSECharmFraction::FillHistos(AliAODRecoDecayHF2Prong *d,TLi
     ((TH1F*)list->FindObject(str.Data()))->Fill(d->ImpParXY()*10000.);
     
     
-    if(vtxTrue){
+    if(fReadMC&&vtxTrue){
       str="hd0D0VtxTrue";
       str.Append(namehist.Data());
       str.Append("SB");
@@ -3208,7 +3238,7 @@ Bool_t AliAnalysisTaskSECharmFraction::FillHistos(AliAODRecoDecayHF2Prong *d,TLi
       
     }
     
-    if(aodDMC!=0x0){
+    if(fReadMC && aodDMC){
       str="hMCd0D0";
       str.Append(namehist.Data());
       str.Append("SB");
@@ -3228,7 +3258,8 @@ Bool_t AliAnalysisTaskSECharmFraction::FillHistos(AliAODRecoDecayHF2Prong *d,TLi
 
 
 void AliAnalysisTaskSECharmFraction::SetNPtBins(Int_t nbins,const Double_t *ptbins){
-  if((fptbins)!=0x0)delete fptbins;
+  // SET THE PT BINS
+  if(fptbins)delete fptbins;
   fnbins=nbins;fptbins=new Double_t[fnbins];
   memcpy(fptbins,ptbins,fnbins*sizeof(Double_t));
   return;