Add event mixing option for background estimation in low pt D-meson task

diff --git a/PWGHF/vertexingHF/AliAnalysisTaskCombinHF.cxx b/PWGHF/vertexingHF/AliAnalysisTaskCombinHF.cxx
index 866116da9bc8a2ede45620e9ace8f9b728a1d513..3d60eaf3d8c03cf00ca8f4aabf456ede6c592e6b 100644 (file)
--- a/PWGHF/vertexingHF/AliAnalysisTaskCombinHF.cxx
+++ b/PWGHF/vertexingHF/AliAnalysisTaskCombinHF.cxx
@@ -69,6 +69,7 @@ AliAnalysisTaskCombinHF::AliAnalysisTaskCombinHF():
   fNormRotated(0x0),
   fDeltaMass(0x0),
   fDeltaMassFullAnalysis(0x0),
+  fMassVsPtVsYME(0x0),
   fFilterMask(BIT(4)),
   fTrackCutsAll(0x0),
   fTrackCutsPion(0x0),
@@ -98,7 +99,22 @@ AliAnalysisTaskCombinHF::AliAnalysisTaskCombinHF():
   fKeepNegID(kFALSE),
   fPIDselCaseZero(0),
   fBayesThresKaon(0.4),
-  fBayesThresPion(0.4)
+  fBayesThresPion(0.4),
+  fDoEventMixing(kTRUE),
+  fMaxNumberOfEventsForMixing(20),
+  fMaxzVertDistForMix(20.),
+  fMaxMultDiffForMix(9999.),
+  fUsePoolsZ(kFALSE),
+  fUsePoolsM(kFALSE),
+  fNzVertPools(1),
+  fzVertPoolLims(0x0),
+  fNMultPools(1),
+  fMultPoolLims(0x0),
+  fEventBuffer(0x0),
+  fVtxZ(0),
+  fMultiplicity(0),
+  fKaonTracks(0x0),
+  fPionTracks(0x0)
 {
   // default constructor
 }
@@ -129,6 +145,7 @@ AliAnalysisTaskCombinHF::AliAnalysisTaskCombinHF(Int_t meson, AliRDHFCuts* analy
   fNormRotated(0x0),
   fDeltaMass(0x0),
   fDeltaMassFullAnalysis(0x0),
+  fMassVsPtVsYME(0x0),
   fFilterMask(BIT(4)),
   fTrackCutsAll(0x0),
   fTrackCutsPion(0x0),
@@ -158,10 +175,25 @@ AliAnalysisTaskCombinHF::AliAnalysisTaskCombinHF(Int_t meson, AliRDHFCuts* analy
   fKeepNegID(kFALSE),
   fPIDselCaseZero(0),
   fBayesThresKaon(0.4),
-  fBayesThresPion(0.4)
+  fBayesThresPion(0.4),
+  fDoEventMixing(kTRUE),
+  fMaxNumberOfEventsForMixing(20),
+  fMaxzVertDistForMix(20.),
+  fMaxMultDiffForMix(9999.),
+  fUsePoolsZ(kFALSE),
+  fUsePoolsM(kFALSE),
+  fNzVertPools(1),
+  fzVertPoolLims(0x0),
+  fNMultPools(1),
+  fMultPoolLims(0x0),
+  fEventBuffer(0x0),
+  fVtxZ(0),
+  fMultiplicity(0),
+  fKaonTracks(0x0),
+  fPionTracks(0x0)
 {
   // standard constructor
-  
+
   DefineOutput(1,TList::Class());  //My private output
   DefineOutput(2,AliNormalizationCounter::Class());
 }
@@ -172,36 +204,68 @@ AliAnalysisTaskCombinHF::~AliAnalysisTaskCombinHF()
   //
   // Destructor
   //
+  if(!fOutput->IsOwner()){
+    delete fHistNEvents;
+    delete fHistTrackStatus;
+    delete fHistCheckOrigin;
+    delete fHistCheckOriginSel;
+    delete fHistCheckDecChan;
+    delete fHistCheckDecChanAcc;
+    delete fPtVsYGen;
+    delete fPtVsYGenLargeAcc;
+    delete fPtVsYGenLimAcc;
+    delete fPtVsYGenAcc;
+    delete fPtVsYReco;
+    delete fMassVsPtVsY;
+    delete fMassVsPtVsYLSpp;
+    delete fMassVsPtVsYLSmm;
+    delete fMassVsPtVsYRot;
+    delete fMassVsPtVsYSig;
+    delete fMassVsPtVsYRefl;
+    delete fMassVsPtVsYBkg;
+    delete fNSelected;
+    delete fNormRotated;
+    delete fDeltaMass;
+    delete fDeltaMassFullAnalysis;
+    delete fMassVsPtVsYME;
+  }
+
   delete fOutput;
-  delete fHistNEvents;
-  delete fHistTrackStatus;
-  delete fHistCheckOrigin;
-  delete fHistCheckOriginSel;
-  delete fHistCheckDecChan;
-  delete fHistCheckDecChanAcc;
-  delete fPtVsYGen;
-  delete fPtVsYGenLargeAcc;
-  delete fPtVsYGenLimAcc;
-  delete fPtVsYGenAcc;
-  delete fPtVsYReco;
-  delete fMassVsPtVsY;
-  delete fMassVsPtVsYLSpp;
-  delete fMassVsPtVsYLSmm;
-  delete fMassVsPtVsYRot;
-  delete fMassVsPtVsYSig;
-  delete fMassVsPtVsYRefl;
-  delete fMassVsPtVsYBkg;
-  delete fNSelected;
-  delete fNormRotated;
-  delete fDeltaMass;
   delete fCounter;
   delete fTrackCutsAll;
   delete fTrackCutsPion;
   delete fTrackCutsKaon;
   delete fPidHF;
   delete fAnalysisCuts;
+  fKaonTracks->Delete();
+  fPionTracks->Delete();
+  delete fKaonTracks;
+  delete fPionTracks;
+  delete fEventBuffer;
 }
 
+//________________________________________________________________________
+void AliAnalysisTaskCombinHF::ConfigureZVertPools(Int_t nPools, Double_t*  zVertLimits)
+{
+  // sets the pools for event mizing in zvertex
+  if(fzVertPoolLims) delete [] fzVertPoolLims;
+  fNzVertPools=nPools;
+  fzVertPoolLims = new Double_t[fNzVertPools+1];
+  for(Int_t ib=0; ib<nPools+1; ib++) fzVertPoolLims[ib]=zVertLimits[ib];
+  fUsePoolsZ=kTRUE;
+  return;  
+}
+//________________________________________________________________________
+void AliAnalysisTaskCombinHF::ConfigureMultiplicityPools(Int_t nPools, Double_t*  multLimits)
+{
+  // sets the pools for event mizing in zvertex
+  if(fMultPoolLims) delete [] fMultPoolLims;
+  fNMultPools=nPools;
+  fMultPoolLims = new Double_t[fNMultPools+1];
+  for(Int_t ib=0; ib<nPools+1; ib++) fMultPoolLims[ib]=multLimits[ib];
+  fUsePoolsM=kTRUE;
+  return;  
+}
 //________________________________________________________________________
 void AliAnalysisTaskCombinHF::UserCreateOutputObjects()
 {
@@ -353,10 +417,26 @@ void AliAnalysisTaskCombinHF::UserCreateOutputObjects()
   fDeltaMassFullAnalysis=new THnSparseF("fDeltaMassFullAnalysis","fDeltaMassFullAnalysis;inv mass (GeV/c);#Delta inv mass (GeV/c) ; p_{T}^{D} (GeV/c); #Delta p_{T} (GeV/c); daughter angle (2prongs) (rad);",5,binSparseDMassRot,edgeLowSparseDMassRot,edgeHighSparseDMassRot);
   fOutput->Add(fDeltaMassFullAnalysis);
   
+  fMassVsPtVsYME=new TH3F("hMassVsPtVsYME","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
+  fMassVsPtVsYME->Sumw2();
+  fMassVsPtVsYME->SetMinimum(0);
+  fOutput->Add(fMassVsPtVsYME);
+ 
+
   //Counter for Normalization
   fCounter = new AliNormalizationCounter("NormalizationCounter");
   fCounter->Init();
   
+  fKaonTracks = new TObjArray();
+  fPionTracks=new TObjArray();
+  fKaonTracks->SetOwner();
+  fPionTracks->SetOwner();
+  fEventBuffer = new TTree("EventBuffer", "Temporary buffer for event mixing");
+  fEventBuffer->Branch("zVertex", &fVtxZ);
+  fEventBuffer->Branch("multiplicity", &fMultiplicity);
+  fEventBuffer->Branch("karray", "TObjArray", &fKaonTracks);
+  fEventBuffer->Branch("parray", "TObjArray", &fPionTracks);
+
   PostData(1,fOutput);
   PostData(2,fCounter);
 }
@@ -423,7 +503,9 @@ void AliAnalysisTaskCombinHF::UserExec(Option_t */*option*/){
   
   
   Int_t ntracks=aod->GetNTracks();
-  
+  fVtxZ = aod->GetPrimaryVertex()->GetZ();
+  fMultiplicity = AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.,1.); 
+
   // select and flag tracks
   UChar_t* status = new UChar_t[ntracks];
   for(Int_t iTr=0; iTr<ntracks; iTr++){
@@ -472,10 +554,16 @@ void AliAnalysisTaskCombinHF::UserExec(Option_t */*option*/){
   Double_t tmpp[3];
   Double_t px[3],py[3],pz[3];
   Int_t dgLabels[3];
-  
+  fKaonTracks->Delete();
+  fPionTracks->Delete();
+ 
   for(Int_t iTr1=0; iTr1<ntracks; iTr1++){
     AliAODTrack* trK=aod->GetTrack(iTr1);
     if((status[iTr1] & 1)==0) continue;
+    if(fDoEventMixing){
+      if(status[iTr1] & 2) fKaonTracks->AddLast(new TLorentzVector(trK->Px(),trK->Py(),trK->Pz(),trK->Charge()));
+      if(status[iTr1] & 4) fPionTracks->AddLast(new TLorentzVector(trK->Px(),trK->Py(),trK->Pz(),trK->Charge()));
+    }
     if((status[iTr1] & 2)==0) continue;
     Int_t chargeK=trK->Charge();
     trK->GetPxPyPz(tmpp);
@@ -542,7 +630,8 @@ void AliAnalysisTaskCombinHF::UserExec(Option_t */*option*/){
   
   fCounter->StoreCandidates(aod,nFiltered,kTRUE);
   fCounter->StoreCandidates(aod,nSelected,kFALSE);
-  
+
+  if(fDoEventMixing) fEventBuffer->Fill();
   PostData(1,fOutput);
   PostData(2,fCounter);
   
@@ -714,6 +803,24 @@ Bool_t AliAnalysisTaskCombinHF::FillHistos(Int_t pdgD,Int_t nProngs, AliAODRecoD
   return accept;
   
 }
+//________________________________________________________________________
+void AliAnalysisTaskCombinHF::FillMEHistos(Int_t pdgD,Int_t nProngs, AliAODRecoDecay* tmpRD, Double_t* px, Double_t* py, Double_t* pz, UInt_t *pdgdau){
+  // Fill histos for candidates in MixedEvents
+    
+  tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
+  Double_t pt = tmpRD->Pt();
+  Double_t minv2 = tmpRD->InvMass2(nProngs,pdgdau);
+  Double_t mass=TMath::Sqrt(minv2);
+  
+  if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
+    Double_t rapid = tmpRD->Y(pdgD);
+    if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
+      fMassVsPtVsYME->Fill(mass,pt,rapid);
+    }
+  }
+  return;
+}
+
 //________________________________________________________________________
 Bool_t AliAnalysisTaskCombinHF::IsTrackSelected(AliAODTrack* track){
   // track selection cuts
@@ -837,7 +944,134 @@ Bool_t AliAnalysisTaskCombinHF::CheckAcceptance(TClonesArray* arrayMC,Int_t nPro
   }
   return kTRUE;
 }
+//_________________________________________________________________
+Bool_t AliAnalysisTaskCombinHF::CanBeMixed(Double_t zv1, Double_t zv2, Double_t mult1, Double_t mult2){
+  if(fUsePoolsZ && fzVertPoolLims){
+    Int_t theBin1=TMath::BinarySearch(fNzVertPools+1,fzVertPoolLims,zv1);
+    Int_t theBin2=TMath::BinarySearch(fNzVertPools+1,fzVertPoolLims,zv2);
+    if(theBin1!=theBin2) return kFALSE;
+    if(theBin1<0 || theBin2<0) return kFALSE;
+    if(theBin1>=fNzVertPools || theBin2>=fNzVertPools) return kFALSE;    
+  }else{
+    if(TMath::Abs(zv2-zv1)>fMaxzVertDistForMix) return kFALSE;
+  }
+
+  if(fUsePoolsM && fMultPoolLims){
+    Int_t theBin1=TMath::BinarySearch(fNMultPools+1,fMultPoolLims,mult1);
+    Int_t theBin2=TMath::BinarySearch(fNMultPools+1,fMultPoolLims,mult2);
+    if(theBin1!=theBin2) return kFALSE;
+    if(theBin1<0 || theBin2<0) return kFALSE;
+    if(theBin1>=fNMultPools || theBin2>=fNMultPools) return kFALSE;
+  }else{
+    if(TMath::Abs(mult2-mult1)>fMaxMultDiffForMix) return kFALSE;
+  }
+  return kTRUE;
+}
+
+//_________________________________________________________________
+void AliAnalysisTaskCombinHF::FinishTaskOutput()
+{
+  // perform mixed event analysis
+  if(!fDoEventMixing) return;
+
+  Int_t nEvents=fEventBuffer->GetEntries();
+  printf("Start Event Mixing of %d events\n",nEvents);
+
+  TObjArray* karray=0x0;
+  TObjArray* parray=0x0;
+  Double_t zVertex,mult;
+  fEventBuffer->SetBranchAddress("karray", &karray);
+  fEventBuffer->SetBranchAddress("parray", &parray);
+  fEventBuffer->SetBranchAddress("zVertex", &zVertex);
+  fEventBuffer->SetBranchAddress("multiplicity", &mult);
+
+  // dummy values of track impact parameter, needed to build an AliAODRecoDecay object
+  Double_t d02[2]={0.,0.};
+  Double_t d03[3]={0.,0.,0.};
+  AliAODRecoDecay* tmpRD2 = new AliAODRecoDecay(0x0,2,0,d02);
+  AliAODRecoDecay* tmpRD3 = new AliAODRecoDecay(0x0,3,1,d03);
+  UInt_t pdg0[2]={321,211};
+  UInt_t pdgp[3]={321,211,211};
+  Double_t px[3],py[3],pz[3];
+
+  for(Int_t iEv1=0; iEv1<nEvents; iEv1++){
+    fEventBuffer->GetEvent(iEv1);
+    Double_t zVertex1=zVertex;
+    Double_t mult1=mult;
+    TObjArray* karray1=new TObjArray(*karray);
+    for(Int_t iEv2=0; iEv2<fMaxNumberOfEventsForMixing; iEv2++){
+      Int_t iToMix=iEv1+iEv2+1;
+      if(iEv1>=(nEvents-fMaxNumberOfEventsForMixing)) iToMix=iEv1-iEv2-1;
+      if(iToMix<0) continue;
+      if(iToMix==iEv1) continue;
+      fEventBuffer->GetEvent(iToMix);
+      Double_t zVertex2=zVertex;
+      Double_t mult2=mult;
+      if(CanBeMixed(zVertex1,zVertex2,mult1,mult2)){
+       TObjArray* parray2=new TObjArray(*parray);
+       Int_t nKaons=karray1->GetEntries();
+       Int_t nPions=parray2->GetEntries();
+       Bool_t doThird=kTRUE;
+       Int_t iToMix3=iEv1+2*fMaxNumberOfEventsForMixing+1;
+       if(iToMix3>=nEvents) iToMix3=iEv1-2*fMaxNumberOfEventsForMixing-1;
+       if(iToMix3<0 || iToMix3==iEv1 || iToMix3==iToMix) doThird=kFALSE;
+       TObjArray* parray3=0x0;
+       Double_t zVertex3=-999.;
+       Double_t mult3=999999;
+       if(fMeson!=kDzero && doThird){
+         fEventBuffer->GetEvent(iToMix3);
+         zVertex3=zVertex;
+         mult3=mult;
+         if(CanBeMixed(zVertex1,zVertex3,mult1,mult3)){
+           parray3=new TObjArray(*parray);
+         }else{
+           doThird=kFALSE;
+         }
+       }
 
+       for(Int_t iTr1=0; iTr1<nKaons; iTr1++){
+         TLorentzVector* trK=(TLorentzVector*)karray1->At(iTr1);
+         Double_t chargeK=trK->T();
+         px[0] = trK->Px();
+         py[0] = trK->Py();
+         pz[0] = trK->Pz();
+         for(Int_t iTr2=0; iTr2<nPions; iTr2++){
+           TLorentzVector* trPi1=(TLorentzVector*)parray2->At(iTr2);
+           Double_t chargePi1=trPi1->T();
+           px[1] = trPi1->Px();
+           py[1] = trPi1->Py();
+           pz[1] = trPi1->Pz();
+           if(chargePi1*chargeK<0){
+             if(fMeson==kDzero){
+               FillMEHistos(421,2,tmpRD2,px,py,pz,pdg0);
+             }else{
+               if(doThird && parray3){
+                 Int_t nPions3=parray3->GetEntries();
+                 for(Int_t iTr3=iTr2+1; iTr3<nPions3; iTr3++){
+                   TLorentzVector* trPi2=(TLorentzVector*)parray3->At(iTr3);
+                   Double_t chargePi2=trPi2->T();
+                   px[2] = trPi2->Px();
+                   py[2] = trPi2->Py();
+                   pz[2] = trPi2->Pz();
+                   if(chargePi2*chargeK<0){
+                     FillMEHistos(411,3,tmpRD3,px,py,pz,pdgp);
+                   }
+                 }
+               }
+             }
+           }
+         }
+       }
+       delete parray3;
+       delete parray2;
+      }
+    }
+    delete karray1;
+  }
+  delete tmpRD2;
+  delete tmpRD3;
+
+}
 //_________________________________________________________________
 void AliAnalysisTaskCombinHF::Terminate(Option_t */*option*/)
 {

diff --git a/PWGHF/vertexingHF/AliAnalysisTaskCombinHF.h b/PWGHF/vertexingHF/AliAnalysisTaskCombinHF.h
index 192af438b632e703c7d3fe0aefca9603163fc927..6c9b7b42626cc5c339b0da6d9ec195e7946a8efa 100644 (file)
--- a/PWGHF/vertexingHF/AliAnalysisTaskCombinHF.h
+++ b/PWGHF/vertexingHF/AliAnalysisTaskCombinHF.h
@@ -34,8 +34,23 @@ public:
   virtual void LocalInit() {Init();}
   virtual void UserExec(Option_t *option);
   virtual void Terminate(Option_t *option);
+  virtual void        FinishTaskOutput();
   
   void SetReadMC(Bool_t read){fReadMC=read;}
+
+  void SetEventMixingOn(){fDoEventMixing=kTRUE;}
+  void SetEventMixingOff(){fDoEventMixing=kFALSE;}
+  void SetMaxNumberOfEventsForMixing(Int_t maxn){fMaxNumberOfEventsForMixing=maxn;}
+  void SetMaxzVertDistForMix(Double_t dist){
+    fUsePoolsZ=kFALSE; 
+    fMaxzVertDistForMix=dist;
+  }
+  void SetMaxMultDiffForMix(Double_t maxd){
+    fUsePoolsM=kFALSE; 
+    fMaxMultDiffForMix=maxd;
+  }
+  void ConfigureZVertPools(Int_t nPools, Double_t*  zVertLimits);
+  void ConfigureMultiplicityPools(Int_t nPools, Double_t*  multLimits);
   void SelectPromptD(){fPromptFeeddown=kPrompt;}
   void SelectFeeddownD(){fPromptFeeddown=kFeeddown;}
   void SelectPromptAndFeeddownD(){fPromptFeeddown=kBoth;}
@@ -89,9 +104,11 @@ public:
   
   Bool_t FillHistos(Int_t pdgD,Int_t nProngs, AliAODRecoDecay* tmpRD, Double_t* px, Double_t* py, Double_t* pz, UInt_t *pdgdau, TClonesArray *arrayMC, Int_t* dgLabels);
   void FillLSHistos(Int_t pdgD,Int_t nProngs, AliAODRecoDecay* tmpRD, Double_t* px, Double_t* py, Double_t* pz, UInt_t *pdgdau, Int_t charge);
+  void FillMEHistos(Int_t pdgD,Int_t nProngs, AliAODRecoDecay* tmpRD, Double_t* px, Double_t* py, Double_t* pz, UInt_t *pdgdau);
   void FillGenHistos(TClonesArray* arrayMC);
-  Bool_t CheckAcceptance(TClonesArray* arrayMC, Int_t nProng, Int_t *labDau);
-  
+  Bool_t CheckAcceptance(TClonesArray* arrayMC, Int_t nProng, Int_t *labDau); 
+  Bool_t CanBeMixed(Double_t zv1, Double_t zv2, Double_t mult1, Double_t mult2); 
+
   enum EMesonSpecies {kDzero, kDplus, kDstar, kDs};
   enum EPrompFd {kNone,kPrompt,kFeeddown,kBoth};
   enum EPIDstrategy {knSigma, kBayesianMaxProb, kBayesianThres};
@@ -124,6 +141,7 @@ private:
   TH1F *fNormRotated;      //! hist. rotated/selected D+
   TH1F *fDeltaMass;        //! hist. mass difference after rotations
   THnSparse *fDeltaMassFullAnalysis; //! hist. mass difference after rotations with more details
+  TH3F *fMassVsPtVsYME;   //! hist. of Y vs. Pt vs. Mass (mixedevents)
   
   UInt_t fFilterMask; // FilterMask
   AliESDtrackCuts* fTrackCutsAll; // track selection
@@ -160,8 +178,24 @@ private:
   Int_t    fPIDselCaseZero;  // flag to change PID strategy
   Double_t fBayesThresKaon;  // threshold for kaon identification via Bayesian PID
   Double_t fBayesThresPion;  // threshold for pion identification via Bayesian PID
-  
-  ClassDef(AliAnalysisTaskCombinHF,5); // D0D+ task from AOD tracks
+
+  Bool_t fDoEventMixing; // flag for event mixing
+  Int_t  fMaxNumberOfEventsForMixing; // maximum number of events to be used in event mixing
+  Double_t fMaxzVertDistForMix; // cut on zvertex distance for event mixing
+  Double_t fMaxMultDiffForMix; // cut on multiplicity difference for event mixing 
+  Bool_t fUsePoolsZ; // flag for using pools in z vertex
+  Bool_t fUsePoolsM; // flag for using pools in multiplicity
+  Int_t fNzVertPools; // number of pools in z vertex for event mixing
+  Double_t* fzVertPoolLims; //[fNzVertPools+1] limits of the pools in zVertex
+  Int_t fNMultPools; // number of pools in multiplicity for event mixing
+  Double_t* fMultPoolLims; //[fNMultPools+1] limits of the pools in multiplicity
+
+  TTree* fEventBuffer;    // structure for event mixing
+  Double_t fVtxZ;         // zVertex
+  Double_t fMultiplicity; // multiplicity
+  TObjArray* fKaonTracks; // array of kaon-compatible tracks (TLorentzVectors)
+  TObjArray* fPionTracks; // array of pion-compatible tracks (TLorentzVectors)  
+  ClassDef(AliAnalysisTaskCombinHF,6); // D0D+ task from AOD tracks
 };
 
 #endif