Nrew task for D vs. multiplcity analysis
authorprino <prino@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 29 Mar 2012 23:17:00 +0000 (23:17 +0000)
committerprino <prino@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 29 Mar 2012 23:17:00 +0000 (23:17 +0000)
PWGHF/CMakelibPWGHFvertexingHF.pkg
PWGHF/PWGHFvertexingHFLinkDef.h
PWGHF/vertexingHF/AliAnalysisTaskSEDvsMultiplicity.cxx [new file with mode: 0644]
PWGHF/vertexingHF/AliAnalysisTaskSEDvsMultiplicity.h [new file with mode: 0644]

index a111d66..ec43124 100644 (file)
@@ -55,6 +55,7 @@ set ( SRCS
     vertexingHF/AliAnalysisTaskSELambdac.cxx 
     vertexingHF/AliAnalysisTaskSED0Mass.cxx 
     vertexingHF/AliAnalysisTaskSECharmFraction.cxx 
+    vertexingHF/AliAnalysisTaskSEDvsMultiplicity.cxx 
     vertexingHF/AliCFVertexingHF.cxx 
     vertexingHF/AliCFVertexingHF2Prong.cxx 
     vertexingHF/AliCFVertexingHF3Prong.cxx 
index a475026..e5fdbc6 100644 (file)
@@ -33,6 +33,7 @@
 #pragma link C++ class AliAnalysisTaskSELambdac+;
 #pragma link C++ class AliAnalysisTaskSED0Mass+;
 #pragma link C++ class AliAnalysisTaskSECharmFraction+;
+#pragma link C++ class AliAnalysisTaskSEDvsMultiplicity+;
 #pragma link C++ class AliCFVertexingHF+;
 #pragma link C++ class AliCFVertexingHF2Prong+;
 #pragma link C++ class AliCFVertexingHF3Prong+;
diff --git a/PWGHF/vertexingHF/AliAnalysisTaskSEDvsMultiplicity.cxx b/PWGHF/vertexingHF/AliAnalysisTaskSEDvsMultiplicity.cxx
new file mode 100644 (file)
index 0000000..1c17c8a
--- /dev/null
@@ -0,0 +1,633 @@
+/**************************************************************************
+ * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+//*************************************************************************
+// Class AliAnalysisTaskSEDvsMultiplicity
+// AliAnalysisTaskSE for the D meson vs. multiplcity analysis
+// Authors: Renu Bala, Zaida Conesa del Valle, Francesco Prino
+/////////////////////////////////////////////////////////////
+
+#include <TClonesArray.h>
+#include <TCanvas.h>
+#include <TList.h>
+#include <TString.h>
+#include <TDatabasePDG.h>
+#include <TH1F.h>
+#include <TH2F.h>     
+#include <TH3F.h>
+#include <THnSparse.h>
+#include <TProfile.h>  
+#include "AliAnalysisManager.h"
+#include "AliRDHFCuts.h"
+#include "AliRDHFCutsDplustoKpipi.h"
+#include "AliRDHFCutsDStartoKpipi.h"
+#include "AliRDHFCutsD0toKpi.h"
+#include "AliAODHandler.h"
+#include "AliAODEvent.h"
+#include "AliAODVertex.h"
+#include "AliAODTrack.h"
+#include "AliAODRecoDecayHF.h"
+#include "AliAODRecoCascadeHF.h"
+#include "AliAnalysisVertexingHF.h"
+#include "AliAnalysisTaskSE.h"
+#include "AliAnalysisTaskSEDvsMultiplicity.h"
+#include "AliNormalizationCounter.h"
+#include "AliVertexingHFUtils.h"
+ClassImp(AliAnalysisTaskSEDvsMultiplicity)
+
+
+//________________________________________________________________________
+AliAnalysisTaskSEDvsMultiplicity::AliAnalysisTaskSEDvsMultiplicity():
+AliAnalysisTaskSE(),
+  fOutput(0),
+  fListCuts(0),
+  fOutputCounters(0),
+  fHistNEvents(0),
+  fPtVsMassVsMult(0),
+  fPtVsMassVsMultNoPid(0),
+  fPtVsMassVsMultUncorr(0),
+  fPtVsMassVsMultPart(0),
+  fPtVsMassVsMultAntiPart(0),
+  fUpmasslimit(1.965),
+  fLowmasslimit(1.765),
+  fBinWidth(0.002),
+  fRDCutsAnalysis(0),
+  fCounter(0),
+  fCounterU(0),
+  fDoImpPar(kFALSE),
+  fNImpParBins(400),
+  fLowerImpPar(-2000.),
+  fHigherImpPar(2000.),
+  fReadMC(kFALSE),
+  fMCOption(0),
+  fUseBit(kTRUE),
+  fRefMult(9.5),
+  fPdgMeson(411)
+{
+   // Default constructor
+  for(Int_t i=0; i<5; i++) fHistMassPtImpPar[i]=0;
+  for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=0;
+}
+
+//________________________________________________________________________
+AliAnalysisTaskSEDvsMultiplicity::AliAnalysisTaskSEDvsMultiplicity(const char *name,AliRDHFCuts *cuts):
+  AliAnalysisTaskSE(name),
+  fOutput(0),
+  fListCuts(0),
+  fOutputCounters(0),
+  fHistNEvents(0),
+  fPtVsMassVsMult(0),
+  fPtVsMassVsMultNoPid(0),
+  fPtVsMassVsMultUncorr(0),
+  fPtVsMassVsMultPart(0),
+  fPtVsMassVsMultAntiPart(0),
+  fUpmasslimit(1.965),
+  fLowmasslimit(1.765),
+  fBinWidth(0.002),
+  fRDCutsAnalysis(cuts),
+  fCounter(0),
+  fCounterU(0),
+  fDoImpPar(kFALSE),
+  fNImpParBins(400),
+  fLowerImpPar(-2000.),
+  fHigherImpPar(2000.),
+  fReadMC(kFALSE),
+  fMCOption(0),
+  fUseBit(kTRUE),
+  fRefMult(9.5),
+  fPdgMeson(411)
+{
+  // 
+  // Standard constructor
+  //
+  for(Int_t i=0; i<5; i++) fHistMassPtImpPar[i]=0;
+  for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=0;
+
+  // Default constructor
+   // Output slot #1 writes into a TList container
+  DefineOutput(1,TList::Class());  //My private output
+  // Output slot #2 writes cut to private output
+  DefineOutput(2,TList::Class());
+  // Output slot #3 writes cut to private output
+  DefineOutput(3,TList::Class());
+ }
+//________________________________________________________________________
+AliAnalysisTaskSEDvsMultiplicity::~AliAnalysisTaskSEDvsMultiplicity()
+{
+  //
+  // Destructor
+  //
+  delete fOutput;
+  delete fHistNEvents;
+  delete fListCuts;
+  delete fRDCutsAnalysis;
+  delete fCounter;
+  delete fCounterU;
+  for(Int_t i=0; i<5; i++){
+    delete fHistMassPtImpPar[i];
+  }
+}  
+
+//_________________________________________________________________
+void  AliAnalysisTaskSEDvsMultiplicity::SetMassLimits(Float_t lowlimit, Float_t uplimit){
+  // set invariant mass limits
+  if(uplimit>lowlimit){
+    Float_t bw=GetBinWidth();
+    fUpmasslimit = lowlimit;
+    fLowmasslimit = uplimit;
+    SetBinWidth(bw);
+  }
+}
+//________________________________________________________________
+void AliAnalysisTaskSEDvsMultiplicity::SetBinWidth(Float_t w){
+  Float_t width=w;
+  Int_t nbins=(Int_t)((fUpmasslimit-fLowmasslimit)/width+0.5);
+  Int_t missingbins=4-nbins%4;
+  nbins=nbins+missingbins;
+  width=(fUpmasslimit-fLowmasslimit)/nbins;
+  if(missingbins!=0){
+    printf("AliAnalysisTaskSEDvsMultiplicity::SetBinWidth: W-bin width of %f will produce histograms not rebinnable by 4. New width set to %f\n",w,width);
+  }
+  else{
+    if(fDebug>1) printf("AliAnalysisTaskSEDvsMultiplicity::SetBinWidth: width set to %f\n",width);
+  }
+  fBinWidth=width;
+}
+//________________________________________________________________________
+void AliAnalysisTaskSEDvsMultiplicity::Init(){
+  //
+  // Initialization
+  //
+  printf("AnalysisTaskSEDvsMultiplicity::Init() \n");
+  
+  fListCuts=new TList();
+
+  if(fPdgMeson==411){
+     AliRDHFCutsDplustoKpipi* copycut=new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fRDCutsAnalysis)));
+    copycut->SetName("AnalysisCutsDplus");
+    fListCuts->Add(copycut);
+  }else if(fPdgMeson==421){
+    AliRDHFCutsD0toKpi* copycut=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fRDCutsAnalysis)));
+    copycut->SetName("AnalysisCutsDzero");
+    fListCuts->Add(copycut);
+  }else if(fPdgMeson==413){
+     AliRDHFCutsDStartoKpipi* copycut=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fRDCutsAnalysis)));
+    copycut->SetName("AnalysisCutsDStar");
+    fListCuts->Add(copycut);
+  }
+  PostData(2,fListCuts);
+  
+  return;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskSEDvsMultiplicity::UserCreateOutputObjects()
+{
+  // Create the output container
+  //
+  if(fDebug > 1) printf("AnalysisTaskSEDvsMultiplicity::UserCreateOutputObjects() \n");
+
+  // Several histograms are more conveniently managed in a TList
+  fOutput = new TList();
+  fOutput->SetOwner();
+  fOutput->SetName("OutputHistos");
+
+
+  TH1F *hspdmultCand = new TH1F("hspdmultCand", "Tracklets multiplicity for events with D candidates; Tracklets ; Entries",200,0.,200.);// Total multiplicity
+  TH1F *hspdmultD = new TH1F("hspdmultD", "Tracklets multiplicity for events with D in mass region ; Tracklets ; Entries",200,0.,200.); // 
+  TH2F *heta16vseta1 = new TH2F("heta16vseta1","Uncorrected Eta1.6 vs Eta1.0; Ntracklets #eta<1.0; Ntracklets #eta<1.6",200,-0.5,199.5,200,-0.5,199.5); //eta 1.6 vs eta 1.0 histogram 
+  TH2F *hNtrkvsVtxZ = new TH2F("hNtrkvsVtxZ","Ntracklet vs VtxZ; VtxZ;N_{tracklet};",300,-15,15,200,0,200.); // 
+  TH2F *hNtrkvsVtxZCorr = new TH2F("hNtrkvsVtxZCorr","Ntracklet vs VtxZ; VtxZ;N_{tracklet};",300,-15,15,200,0,200.); // 
+  
+
+  hspdmultCand->Sumw2();
+  hspdmultD->Sumw2();
+
+  fOutput->Add(hspdmultCand);
+  fOutput->Add(hspdmultD);
+  fOutput->Add(heta16vseta1);
+  fOutput->Add(hNtrkvsVtxZ); 
+  fOutput->Add(hNtrkvsVtxZCorr); 
+
+  
+  fHistNEvents = new TH1F("fHistNEvents", "number of events ",11,-0.5,10.5);
+  fHistNEvents->GetXaxis()->SetBinLabel(1,"nEvents total");
+  fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents with Z vertex");
+  fHistNEvents->GetXaxis()->SetBinLabel(3,"nEvents selected");
+  fHistNEvents->GetXaxis()->SetBinLabel(4,"Rejected due to trigger");
+  fHistNEvents->GetXaxis()->SetBinLabel(5,"Rejected due to phys sel");
+  fHistNEvents->GetXaxis()->SetBinLabel(6,"Rejected due to vertex cuts");
+  fHistNEvents->GetXaxis()->SetBinLabel(7,"Rejected due to pileup");
+  fHistNEvents->GetXaxis()->SetBinLabel(8,"Total no. of candidate");
+  fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of cand wo bitmask");
+  fHistNEvents->GetXaxis()->SetBinLabel(10,"D after cuts (No PID)");
+  fHistNEvents->GetXaxis()->SetBinLabel(11,"D after cuts + PID)"); 
+  fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);  
+  fHistNEvents->Sumw2();
+  fHistNEvents->SetMinimum(0);
+  fOutput->Add(fHistNEvents);
+
+  fPtVsMassVsMult=new TH3F("hPtVsMassvsMult", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",200,0.,200.,200,fLowmasslimit,fUpmasslimit,48,0.,24.);
+  fPtVsMassVsMultNoPid=new TH3F("hPtVsMassvsMultNoPid", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",200,0.,200.,200,fLowmasslimit,fUpmasslimit,48,0.,24.); 
+
+  fPtVsMassVsMultUncorr=new TH3F("hPtVsMassvsMultUncorr", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",200,-0.5,199.5,200,fLowmasslimit,fUpmasslimit,48,0.,24.);
+
+  fPtVsMassVsMultPart=new TH3F("hPtVsMassvsMultPart", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",200,0.,200.,200,fLowmasslimit,fUpmasslimit,48,0.,24.);
+
+  fPtVsMassVsMultAntiPart=new TH3F("hPtVsMassvsMultAntiPart", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",200,0.,200.,200,fLowmasslimit,fUpmasslimit,48,0.,24.);
+
+  fOutput->Add(fPtVsMassVsMult);
+  fOutput->Add(fPtVsMassVsMultUncorr);
+  fOutput->Add(fPtVsMassVsMultNoPid);
+  fOutput->Add(fPtVsMassVsMultPart);
+  fOutput->Add(fPtVsMassVsMultAntiPart);
+
+  if(fDoImpPar) CreateImpactParameterHistos();
+
+  fCounter = new AliNormalizationCounter("NormCounterCorrMult");
+  fCounter->SetStudyMultiplicity(kTRUE,1.);
+  fCounter->Init(); 
+
+  fCounterU = new AliNormalizationCounter("NormCounterUnCorrMult");
+  fCounterU->SetStudyMultiplicity(kTRUE,1.);
+  fCounterU->Init(); 
+
+  fOutputCounters = new TList();
+  fOutputCounters->SetOwner();
+  fOutputCounters->SetName("OutputCounters");
+  fOutputCounters->Add(fCounter);
+  fOutputCounters->Add(fCounterU);
+  
+  PostData(1,fOutput); 
+  PostData(2,fListCuts);
+  PostData(3,fOutputCounters);
+
+  return;
+}
+
+
+
+//________________________________________________________________________
+void AliAnalysisTaskSEDvsMultiplicity::UserExec(Option_t */*option*/)
+{
+  // Execute analysis for current event:
+  // heavy flavor candidates association to MC truth
+
+  printf("UserExec\n");
+
+  AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
+  
+  //  AliAODTracklets* tracklets = aod->GetTracklets();
+  //Int_t ntracklets = tracklets->GetNumberOfTracklets();
+  
+  TClonesArray *arrayCand = 0;
+  TString arrayName="";
+  UInt_t pdgDau[3];
+  Int_t nDau=0;
+  Int_t selbit=0;
+  if(fPdgMeson==411){
+    arrayName="Charm3Prong";
+    pdgDau[0]=211; pdgDau[1]=321; pdgDau[2]=211; 
+    nDau=3;
+    selbit=AliRDHFCuts::kDplusCuts;
+  }else if(fPdgMeson==421){
+    arrayName="D0toKpi";
+    pdgDau[0]=211; pdgDau[1]=321; pdgDau[2]=0;
+    nDau=2;
+    selbit=AliRDHFCuts::kD0toKpiCuts;
+  }else if(fPdgMeson==413){
+    arrayName="Dstar";
+    pdgDau[0]=321; pdgDau[1]=211; pdgDau[2]=211; 
+    nDau=3;
+    selbit=AliRDHFCuts::kDstarCuts;
+  }
+  printf("Get delta AOD\n");
+  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();
+      arrayCand=(TClonesArray*)aodFromExt->GetList()->FindObject(arrayName.Data());
+    }
+  } else if(aod) {
+    arrayCand=(TClonesArray*)aod->GetList()->FindObject(arrayName.Data());
+  }
+  printf("delta AOD OK\n");
+
+  if(!aod || !arrayCand) {
+    printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: Charm3Prong branch not found!\n");
+    return;
+  }
+
+
+
+  // fix for temporary bug in ESDfilter 
+  // the AODs with null vertex pointer didn't pass the PhysSel
+  if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return;
+
+  Int_t countTreta1=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.,1.);
+  Int_t countTr=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.6,1.6);
+
+  fCounterU->StoreEvent(aod,fRDCutsAnalysis,fReadMC,countTreta1);
+  fHistNEvents->Fill(0); // count event
+
+  AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
+  TString primTitle = vtx1->GetTitle();
+  Double_t countTreta1corr=countTreta1;
+  if(vtx1 && vtx1->GetNContributors()>0){    
+    fHistNEvents->Fill(1); 
+    TProfile* estimatorAvg = GetEstimatorHistogram(aod);
+    if(estimatorAvg){
+      countTreta1corr=AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,countTreta1,vtx1->GetZ(),fRefMult); 
+    }
+  }
+   
+  fCounter->StoreEvent(aod,fRDCutsAnalysis,fReadMC,countTreta1corr);
+
+  Bool_t isEvSel=fRDCutsAnalysis->IsEventSelected(aod);
+
+  if(fRDCutsAnalysis->GetWhyRejection()==5) fHistNEvents->Fill(3);
+  if(fRDCutsAnalysis->GetWhyRejection()==7) fHistNEvents->Fill(4); 
+  if(fRDCutsAnalysis->GetWhyRejection()==6)fHistNEvents->Fill(5);
+  if(fRDCutsAnalysis->GetWhyRejection()==1)fHistNEvents->Fill(6);
+
+  
+  if(!isEvSel)return;
+  fHistNEvents->Fill(2); // count events selected
+  ((TH2F*)(fOutput->FindObject("heta16vseta1")))->Fill(countTreta1,countTr);
+  ((TH2F*)(fOutput->FindObject("heta1corrvseta1")))->Fill(countTreta1,countTreta1corr);
+  ((TH2F*)(fOutput->FindObject("hNtrkvsVtxZ")))->Fill(vtx1->GetZ(),countTreta1);
+  ((TH2F*)(fOutput->FindObject("hNtrkvsVtxZCorr")))->Fill(vtx1->GetZ(),countTreta1corr);
+  
+  TClonesArray *arrayMC=0;
+  AliAODMCHeader *mcHeader=0;
+
+  // load MC particles
+  if(fReadMC){
+     
+    arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
+    if(!arrayMC) {
+      printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: MC particles branch not found!\n");
+      return;
+    }  
+    // load MC header
+    mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
+    if(!mcHeader) {
+      printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: MC header branch not found!\n");
+      return;
+     }
+  }
+  
+  Int_t nCand = arrayCand->GetEntriesFast(); 
+  Int_t nSelectedNoPID=0,nSelectedPID=0;
+
+  for (Int_t iCand = 0; iCand < nCand; iCand++) {
+    AliAODRecoDecayHF *d = (AliAODRecoDecayHF*)arrayCand->UncheckedAt(iCand);
+    fHistNEvents->Fill(7);
+    if(fUseBit && !d->HasSelectionBit(selbit)){
+      fHistNEvents->Fill(8);
+      continue;
+    }
+    
+    Double_t ptCand = d->Pt();
+    Double_t rapid=d->Y(fPdgMeson);
+    Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
+    Int_t passAllCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate,aod);
+    Int_t passTopolCuts=fRDCutsAnalysis->GetIsSelectedCuts();
+    if(passTopolCuts==0) continue;
+    nSelectedNoPID++;
+    fHistNEvents->Fill(9);
+    if(passAllCuts){
+      nSelectedPID++;
+      fHistNEvents->Fill(10);
+    }
+    Bool_t isPrimary=kTRUE;
+    Int_t labD=-1;
+    Double_t trueImpParXY=9999.;
+    Double_t impparXY=d->ImpParXY()*10000.;
+    Double_t dlen=0.1; //FIXME
+    Double_t mass[2];
+    if(fPdgMeson==411){
+      mass[0]=d->InvMass(nDau,pdgDau);
+      mass[1]=-1.;
+    }else if(fPdgMeson==421){
+      UInt_t pdgdaughtersD0[2]={211,321};//pi,K 
+      UInt_t pdgdaughtersD0bar[2]={321,211};//K,pi 
+      mass[0]=d->InvMass(2,pdgdaughtersD0);
+      mass[1]=d->InvMass(2,pdgdaughtersD0bar);
+    }else if(fPdgMeson==413){
+      // FIXME
+      AliAODRecoCascadeHF* temp = (AliAODRecoCascadeHF*)d;
+      mass[0]=temp->DeltaInvMass();
+      mass[1]=-1.;
+    }
+    for(Int_t iHyp=0; iHyp<2; iHyp++){
+      if(mass[iHyp]<0.) continue; // for D+ and D* we have 1 mass hypothesis
+      Double_t invMass=mass[iHyp];
+      Double_t arrayForSparse[5]={invMass,ptCand,impparXY,dlen,countTreta1corr};
+
+      if(fReadMC){
+       
+       labD = d->MatchToMC(fPdgMeson,arrayMC,nDau,(Int_t*)pdgDau);
+       Bool_t fillHisto=fDoImpPar;
+       if(labD>=0){
+         AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
+         Int_t code=partD->GetPdgCode();
+         if(CheckOrigin(arrayMC,partD)==5) isPrimary=kFALSE;
+         if(code<0 && iHyp==0) fillHisto=kFALSE;
+         if(code>0 && iHyp==1) fillHisto=kFALSE;
+         if(!isPrimary){
+           if(fPdgMeson==411){
+             trueImpParXY=AliVertexingHFUtils::GetTrueImpactParameterDplus(mcHeader,arrayMC,partD)*10000.;
+           }else if(fPdgMeson==421){
+             trueImpParXY=0.; /// FIXME
+           }else if(fPdgMeson==413){
+             trueImpParXY=0.; /// FIXME
+           }
+           Double_t arrayForSparseTrue[5]={invMass,ptCand,trueImpParXY,dlen,countTreta1corr};
+           if(fillHisto && isFidAcc && passAllCuts){
+             fHistMassPtImpPar[2]->Fill(arrayForSparse);
+             fHistMassPtImpPar[3]->Fill(arrayForSparseTrue);
+           }
+         }else{
+           if(fillHisto && isFidAcc && passAllCuts) fHistMassPtImpPar[1]->Fill(arrayForSparse);
+         }
+       }else{
+         if(fillHisto && isFidAcc && passAllCuts)fHistMassPtImpPar[4]->Fill(arrayForSparse);
+       }
+       if(fPdgMeson==421){
+         if(TMath::Abs(labD)==fPdgMeson && fMCOption==2) continue;
+         if(TMath::Abs(labD)!=fPdgMeson && fMCOption==1) continue;      
+       }
+      }
+      
+      if(iHyp==0 && !(passTopolCuts&1)) continue; // candidate not passing as D0
+      if(iHyp==1 && !(passTopolCuts&2)) continue; // candidate not passing as D0bar
+
+      if(isFidAcc){
+       fPtVsMassVsMultNoPid->Fill(countTreta1corr,invMass,ptCand);
+       if(passAllCuts){
+         fPtVsMassVsMult->Fill(countTreta1corr,invMass,ptCand);                      
+         fPtVsMassVsMultUncorr->Fill(countTreta1,invMass,ptCand);
+         // Add separation between part antipart
+         if(fPdgMeson==411){
+           if(d->GetCharge()>0) fPtVsMassVsMultPart->Fill(countTreta1corr,invMass,ptCand);
+           else fPtVsMassVsMultAntiPart->Fill(countTreta1corr,invMass,ptCand);
+         }else if(fPdgMeson==421){
+           if(passTopolCuts&1) fPtVsMassVsMultPart->Fill(countTreta1corr,invMass,ptCand);
+           if(passTopolCuts&2) fPtVsMassVsMultAntiPart->Fill(countTreta1corr,invMass,ptCand);
+         }else if(fPdgMeson==413){
+           // FIXME ADD Dstar!!!!!!!!
+         }
+      
+       
+         if(fDoImpPar){
+           fHistMassPtImpPar[0]->Fill(arrayForSparse);
+         }
+         
+       }
+      }
+
+    }
+  }
+  fCounter->StoreCandidates(aod,nSelectedNoPID,kTRUE);
+  fCounter->StoreCandidates(aod,nSelectedPID,kFALSE);
+  
+    
+  PostData(1,fOutput); 
+  PostData(2,fListCuts);
+  PostData(3,fOutput);
+  
+  
+  return;
+}
+//_________________________________________________________________
+Int_t AliAnalysisTaskSEDvsMultiplicity::GetNMassBins() const {
+  return (Int_t)((fUpmasslimit-fLowmasslimit)/fBinWidth+0.5);
+}
+//________________________________________________________________________
+void AliAnalysisTaskSEDvsMultiplicity::CreateImpactParameterHistos(){
+  // Histos for impact paramter study
+
+  Int_t nmassbins=GetNMassBins();
+  Int_t nbins[5]={nmassbins,200,fNImpParBins,50,100};
+  Double_t xmin[5]={fLowmasslimit,0.,fLowerImpPar,0.,-0.5};
+  Double_t xmax[5]={fUpmasslimit,20.,fHigherImpPar,1.,99.5};
+
+  fHistMassPtImpPar[0]=new THnSparseF("hMassPtImpParAll",
+                                       "Mass vs. pt vs.imppar - All",
+                                       5,nbins,xmin,xmax);
+  fHistMassPtImpPar[1]=new THnSparseF("hMassPtImpParPrompt",
+                                       "Mass vs. pt vs.imppar - promptD",
+                                       5,nbins,xmin,xmax);
+  fHistMassPtImpPar[2]=new THnSparseF("hMassPtImpParBfeed",
+                                       "Mass vs. pt vs.imppar - DfromB",
+                                       5,nbins,xmin,xmax);
+  fHistMassPtImpPar[3]=new THnSparseF("hMassPtImpParTrueBfeed",
+                                       "Mass vs. pt vs.true imppar -DfromB",
+                                       5,nbins,xmin,xmax);
+  fHistMassPtImpPar[4]=new THnSparseF("hMassPtImpParBkg",
+                                       "Mass vs. pt vs.imppar - backgr.",
+                                       5,nbins,xmin,xmax);
+  for(Int_t i=0; i<5;i++){
+    fOutput->Add(fHistMassPtImpPar[i]);
+  }
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskSEDvsMultiplicity::Terminate(Option_t */*option*/)
+{
+  // Terminate analysis
+  //
+  if(fDebug > 1) printf("AnalysisTaskSEDvsMultiplicity: Terminate() \n");
+
+  fOutput = dynamic_cast<TList*> (GetOutputData(1));
+  if (!fOutput) {     
+    printf("ERROR: fOutput not available\n");
+    return;
+  }
+  fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
+
+
+  return;
+}
+//_________________________________________________________________________________________________
+Int_t AliAnalysisTaskSEDvsMultiplicity::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {          
+  //
+  // checking whether the mother of the particles come from a charm or a bottom quark
+  //
+       
+  Int_t pdgGranma = 0;
+  Int_t mother = 0;
+  mother = mcPartCandidate->GetMother();
+  Int_t istep = 0;
+  Int_t abspdgGranma =0;
+  Bool_t isFromB=kFALSE;
+  Bool_t isQuarkFound=kFALSE;
+  while (mother >0 ){
+    istep++;
+    AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
+    if (mcGranma){
+      pdgGranma = mcGranma->GetPdgCode();
+      abspdgGranma = TMath::Abs(pdgGranma);
+      if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
+       isFromB=kTRUE;
+      }
+      if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
+      mother = mcGranma->GetMother();
+    }else{
+      AliError("Failed casting the mother particle!");
+      break;
+    }
+  }
+  
+  if(isFromB) return 5;
+  else return 4;
+}
+
+
+
+//____________________________________________________________________________
+TProfile* AliAnalysisTaskSEDvsMultiplicity::GetEstimatorHistogram(const AliVEvent* event){
+  // Get Estimator Histogram from period event->GetRunNumber();
+  //
+  // If you select SPD tracklets in |eta|<1 you should use type == 1
+  //
+
+  Int_t runNo  = event->GetRunNumber();
+  Int_t period = -1;   // 0-LHC10b, 1-LHC10c, 2-LHC10d, 3-LHC10e
+  if(runNo>114930 && runNo<117223) period = 0;
+  if(runNo>119158 && runNo<120830) period = 1;
+  if(runNo>122373 && runNo<126438) period = 2;
+  if(runNo>127711 && runNo<130841) period = 3;
+  if(period<0 || period>3) return 0;
+
+  return fMultEstimatorAvg[period];
+}
+
diff --git a/PWGHF/vertexingHF/AliAnalysisTaskSEDvsMultiplicity.h b/PWGHF/vertexingHF/AliAnalysisTaskSEDvsMultiplicity.h
new file mode 100644 (file)
index 0000000..85fe7d5
--- /dev/null
@@ -0,0 +1,133 @@
+#ifndef ALIANALYSISTASKSEDVSMULTIPLICITY_H
+#define ALIANALYSISTASKSEDVSMULTIPLICITY_H
+
+/* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */ 
+
+//*************************************************************************
+// Class AliAnalysisTaskSEDvsMultiplicity
+// AliAnalysisTaskSE for the D meson vs. multiplcity analysis
+// Authors: Renu Bala, Zaida Conesa del Valle, Francesco Prino
+//*************************************************************************
+
+#include <TROOT.h>
+#include <TSystem.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TH3F.h>
+#include <TArrayD.h>
+#include <TFile.h>
+#include <TRandom.h>
+#include <TProfile.h>
+#include "AliRDHFCutsDplustoKpipi.h"
+#include "AliAnalysisTaskSE.h"
+#include "AliAnalysisVertexingHF.h"
+#include "AliNormalizationCounter.h"
+#include "AliAODMCHeader.h"
+#include "AliAODMCParticle.h"
+#include "AliVertexingHFUtils.h"
+#include "AliVEvent.h"
+
+
+
+class AliAnalysisTaskSEDvsMultiplicity : public AliAnalysisTaskSE
+{
+ public:
+
+  AliAnalysisTaskSEDvsMultiplicity();
+  AliAnalysisTaskSEDvsMultiplicity(const char *name, AliRDHFCuts* cuts);
+  virtual ~AliAnalysisTaskSEDvsMultiplicity();
+
+  void SetReadMC(Bool_t readMC=kTRUE){fReadMC=readMC;}
+
+  void SetMassLimits(Float_t lowlimit, Float_t uplimit);
+  Float_t GetUpperMassLimit(){return fUpmasslimit;}
+  Float_t GetLowerMassLimit(){return fLowmasslimit;}
+  Int_t GetNMassBins() const;
+
+  void SetBinWidth(Float_t w);
+  Float_t GetBinWidth(){return fBinWidth;}
+
+  void SetImpactParameterBinning(Int_t nbins, Float_t dmin, Float_t dmax){
+    fNImpParBins=nbins;
+    fLowerImpPar=dmin;
+    fHigherImpPar=dmax;
+  }
+
+  void SetUseBit(Bool_t use=kTRUE){fUseBit=use;}
+  void SetDoImpactParameterHistos(Bool_t doImp=kTRUE){fDoImpPar=doImp;}
+
+
+  void SetMultiplVsZProfileLHC10b(TProfile* hprof){
+    fMultEstimatorAvg[0]=hprof;
+  }
+  void SetMultiplVsZProfileLHC10c(TProfile* hprof){
+    fMultEstimatorAvg[1]=hprof;
+  }
+  void SetMultiplVsZProfileLHC10d(TProfile* hprof){
+    fMultEstimatorAvg[2]=hprof;
+  }
+  void SetMultiplVsZProfileLHC10e(TProfile* hprof){
+    fMultEstimatorAvg[3]=hprof;
+  }
+
+
+  Int_t CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const;
+
+  // Implementation of interface methods
+  virtual void UserCreateOutputObjects();
+  virtual void Init();
+  virtual void LocalInit() {Init();}
+  virtual void UserExec(Option_t *option);
+  virtual void Terminate(Option_t *option);
+    
+ private:
+
+  AliAnalysisTaskSEDvsMultiplicity(const AliAnalysisTaskSEDvsMultiplicity &source);
+  AliAnalysisTaskSEDvsMultiplicity& operator=(const AliAnalysisTaskSEDvsMultiplicity& source); 
+
+  TProfile* GetEstimatorHistogram(const AliVEvent *event);
+  void CreateImpactParameterHistos();
+
+
+  TList  *fOutput; //! list send on output slot 1
+  TList  *fListCuts; //list of cuts
+  TList  *fOutputCounters; //! list send on output slot 3
+
+  TH1F *fHistNEvents; //!hist. for No. of events
+  TH3F *fPtVsMassVsMult;  //! hist. of Pt vs Mult vs. mass (
+  TH3F *fPtVsMassVsMultNoPid;  //! hist. of Pt vs Mult vs. mass (no pid)
+  TH3F *fPtVsMassVsMultUncorr;  //! hist. of Pt vs Mult vs. mass (raw mult)
+  TH3F *fPtVsMassVsMultPart;  //! hist. of Pt vs Mult vs. mass (particle)
+  TH3F *fPtVsMassVsMultAntiPart;  //! hist. of Pt vs Mult vs. mass (antiparticle)
+
+  THnSparseF *fHistMassPtImpPar[5];//! histograms for impact paramter studies
+
+  Float_t fUpmasslimit;  //upper inv mass limit for histos
+  Float_t fLowmasslimit; //lower inv mass limit for histos
+  Float_t fBinWidth;//width of one bin in output histos
+
+  AliRDHFCuts *fRDCutsAnalysis; // Cuts for Analysis
+  AliNormalizationCounter *fCounter;  //!Counter for normalization
+  AliNormalizationCounter *fCounterU; //!Counter for normalization, uncorr mult
+
+  Bool_t fDoImpPar;  //swicth for D impact parameter THnSparse
+  Int_t  fNImpParBins;   // nunber of bins in impact parameter histos
+  Float_t fLowerImpPar;  // lower limit in impact parameter (um)
+  Float_t fHigherImpPar; // higher limit in impact parameter (um)
+
+  Bool_t fReadMC;    //flag for access to MC
+  Int_t  fMCOption;  // 0=keep all cand, 1=keep only signal, 2= keep only back
+  Bool_t fUseBit;    // flag to use bitmask
+  
+  TProfile* fMultEstimatorAvg[4]; // TProfile with mult vs. Z per period
+  Double_t fRefMult;   // refrence multiplcity (period b)
+  Int_t fPdgMeson;   // pdg code of analyzed meson
+
+  
+   ClassDef(AliAnalysisTaskSEDvsMultiplicity,1); // AliAnalysisTaskSE for the MC association of heavy-flavour decay candidates
+};
+
+#endif