PWGHF/vertexingHF/AliAnalysisTaskSEDvsMultiplicity.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 //*************************************************************************
  19 // Class AliAnalysisTaskSEDvsMultiplicity
  20 // AliAnalysisTaskSE for the D meson vs. multiplcity analysis
  21 // Authors: Renu Bala, Zaida Conesa del Valle, Francesco Prino
  22 /////////////////////////////////////////////////////////////
  23
  24 #include <TClonesArray.h>
  25 #include <TCanvas.h>
  26 #include <TList.h>
  27 #include <TString.h>
  28 #include <TDatabasePDG.h>
  29 #include <TH1F.h>
  30 #include <TH2F.h>
  31 #include <TH3F.h>
  32 #include <THnSparse.h>
  33 #include <TProfile.h>
  34 #include "AliAnalysisManager.h"
  35 #include "AliRDHFCuts.h"
  36 #include "AliRDHFCutsDplustoKpipi.h"
  37 #include "AliRDHFCutsDStartoKpipi.h"
  38 #include "AliRDHFCutsD0toKpi.h"
  39 #include "AliAODHandler.h"
  40 #include "AliAODEvent.h"
  41 #include "AliAODVertex.h"
  42 #include "AliAODTrack.h"
  43 #include "AliAODRecoDecayHF.h"
  44 #include "AliAODRecoCascadeHF.h"
  45 #include "AliAnalysisVertexingHF.h"
  46 #include "AliAnalysisTaskSE.h"
  47 #include "AliAnalysisTaskSEDvsMultiplicity.h"
  48 #include "AliNormalizationCounter.h"
  49 #include "AliVertexingHFUtils.h"
  50 ClassImp(AliAnalysisTaskSEDvsMultiplicity)
  51
  52
  53 //________________________________________________________________________
  54 AliAnalysisTaskSEDvsMultiplicity::AliAnalysisTaskSEDvsMultiplicity():
  55 AliAnalysisTaskSE(),
  56   fOutput(0),
  57   fListCuts(0),
  58   fOutputCounters(0),
  59   fHistNEvents(0),
  60   fPtVsMassVsMult(0),
  61   fPtVsMassVsMultNoPid(0),
  62   fPtVsMassVsMultUncorr(0),
  63   fPtVsMassVsMultPart(0),
  64   fPtVsMassVsMultAntiPart(0),
  65   fUpmasslimit(1.965),
  66   fLowmasslimit(1.765),
  67   fNMassBins(200),
  68   fRDCutsAnalysis(0),
  69   fCounter(0),
  70   fCounterU(0),
  71   fDoImpPar(kFALSE),
  72   fNImpParBins(400),
  73   fLowerImpPar(-2000.),
  74   fHigherImpPar(2000.),
  75   fReadMC(kFALSE),
  76   fMCOption(0),
  77   fUseBit(kTRUE),
  78   fRefMult(9.5),
  79   fPdgMeson(411)
  80 {
  81    // Default constructor
  82   for(Int_t i=0; i<5; i++) fHistMassPtImpPar[i]=0;
  83   for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=0;
  84 }
  85
  86 //________________________________________________________________________
  87 AliAnalysisTaskSEDvsMultiplicity::AliAnalysisTaskSEDvsMultiplicity(const char *name, Int_t pdgMeson,AliRDHFCuts *cuts):
  88   AliAnalysisTaskSE(name),
  89   fOutput(0),
  90   fListCuts(0),
  91   fOutputCounters(0),
  92   fHistNEvents(0),
  93   fPtVsMassVsMult(0),
  94   fPtVsMassVsMultNoPid(0),
  95   fPtVsMassVsMultUncorr(0),
  96   fPtVsMassVsMultPart(0),
  97   fPtVsMassVsMultAntiPart(0),
  98   fUpmasslimit(1.965),
  99   fLowmasslimit(1.765),
 100   fNMassBins(200),
 101   fRDCutsAnalysis(cuts),
 102   fCounter(0),
 103   fCounterU(0),
 104   fDoImpPar(kFALSE),
 105   fNImpParBins(400),
 106   fLowerImpPar(-2000.),
 107   fHigherImpPar(2000.),
 108   fReadMC(kFALSE),
 109   fMCOption(0),
 110   fUseBit(kTRUE),
 111   fRefMult(9.5),
 112   fPdgMeson(pdgMeson)
 113 {
 114   //
 115   // Standard constructor
 116   //
 117   for(Int_t i=0; i<5; i++) fHistMassPtImpPar[i]=0;
 118   for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=0;
 119   if(fPdgMeson==413){
 120     fNMassBins=200; // FIXME
 121     SetMassLimits(0.,0.2); // FIXME
 122   }else{
 123     fNMassBins=200;
 124     SetMassLimits(fPdgMeson,0.1);
 125   }
 126   // Default constructor
 127    // Output slot #1 writes into a TList container
 128   DefineOutput(1,TList::Class());  //My private output
 129   // Output slot #2 writes cut to private output
 130   DefineOutput(2,TList::Class());
 131   // Output slot #3 writes cut to private output
 132   DefineOutput(3,TList::Class());
 133 }
 134 //________________________________________________________________________
 135 AliAnalysisTaskSEDvsMultiplicity::~AliAnalysisTaskSEDvsMultiplicity()
 136 {
 137   //
 138   // Destructor
 139   //
 140   delete fOutput;
 141   delete fHistNEvents;
 142   delete fListCuts;
 143   delete fRDCutsAnalysis;
 144   delete fCounter;
 145   delete fCounterU;
 146   for(Int_t i=0; i<5; i++){
 147     delete fHistMassPtImpPar[i];
 148   }
 149 }
 150
 151 //_________________________________________________________________
 152 void  AliAnalysisTaskSEDvsMultiplicity::SetMassLimits(Double_t lowlimit, Double_t uplimit){
 153   // set invariant mass limits
 154   if(uplimit>lowlimit){
 155     fUpmasslimit = lowlimit;
 156     fLowmasslimit = uplimit;
 157   }else{
 158     AliError("Wrong mass limits: upper value should be larger than lower one");
 159   }
 160 }
 161 //_________________________________________________________________
 162 void  AliAnalysisTaskSEDvsMultiplicity::SetMassLimits(Int_t pdg, Double_t range){
 163   // set invariant mass limits
 164   Double_t mass=TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg))->Mass();
 165   SetMassLimits(mass-range,mass+range);
 166 }
 167 //________________________________________________________________________
 168 void AliAnalysisTaskSEDvsMultiplicity::Init(){
 169   //
 170   // Initialization
 171   //
 172   printf("AnalysisTaskSEDvsMultiplicity::Init() \n");
 173
 174   fListCuts=new TList();
 175
 176   if(fPdgMeson==411){
 177      AliRDHFCutsDplustoKpipi* copycut=new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fRDCutsAnalysis)));
 178     copycut->SetName("AnalysisCutsDplus");
 179     fListCuts->Add(copycut);
 180   }else if(fPdgMeson==421){
 181     AliRDHFCutsD0toKpi* copycut=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fRDCutsAnalysis)));
 182     copycut->SetName("AnalysisCutsDzero");
 183     fListCuts->Add(copycut);
 184   }else if(fPdgMeson==413){
 185      AliRDHFCutsDStartoKpipi* copycut=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fRDCutsAnalysis)));
 186     copycut->SetName("AnalysisCutsDStar");
 187     fListCuts->Add(copycut);
 188   }
 189   PostData(2,fListCuts);
 190
 191   return;
 192 }
 193
 194 //________________________________________________________________________
 195 void AliAnalysisTaskSEDvsMultiplicity::UserCreateOutputObjects()
 196 {
 197   // Create the output container
 198   //
 199   if(fDebug > 1) printf("AnalysisTaskSEDvsMultiplicity::UserCreateOutputObjects() \n");
 200
 201   // Several histograms are more conveniently managed in a TList
 202   fOutput = new TList();
 203   fOutput->SetOwner();
 204   fOutput->SetName("OutputHistos");
 205
 206
 207   TH1F *hspdmultCand = new TH1F("hspdmultCand", "Tracklets multiplicity for events with D candidates; Tracklets ; Entries",200,0.,200.);// Total multiplicity
 208   TH1F *hspdmultD = new TH1F("hspdmultD", "Tracklets multiplicity for events with D in mass region ; Tracklets ; Entries",200,0.,200.); //
 209   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
 210   TH2F *hNtrkvsVtxZ = new TH2F("hNtrkvsVtxZ","Ntracklet vs VtxZ; VtxZ;N_{tracklet};",300,-15,15,200,0,200.); //
 211   TH2F *hNtrkvsVtxZCorr = new TH2F("hNtrkvsVtxZCorr","Ntracklet vs VtxZ; VtxZ;N_{tracklet};",300,-15,15,200,0,200.); //
 212
 213
 214   hspdmultCand->Sumw2();
 215   hspdmultD->Sumw2();
 216
 217   fOutput->Add(hspdmultCand);
 218   fOutput->Add(hspdmultD);
 219   fOutput->Add(heta16vseta1);
 220   fOutput->Add(hNtrkvsVtxZ);
 221   fOutput->Add(hNtrkvsVtxZCorr);
 222
 223
 224   fHistNEvents = new TH1F("fHistNEvents", "number of events ",11,-0.5,10.5);
 225   fHistNEvents->GetXaxis()->SetBinLabel(1,"nEvents total");
 226   fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents with Z vertex");
 227   fHistNEvents->GetXaxis()->SetBinLabel(3,"nEvents selected");
 228   fHistNEvents->GetXaxis()->SetBinLabel(4,"Rejected due to trigger");
 229   fHistNEvents->GetXaxis()->SetBinLabel(5,"Rejected due to phys sel");
 230   fHistNEvents->GetXaxis()->SetBinLabel(6,"Rejected due to vertex cuts");
 231   fHistNEvents->GetXaxis()->SetBinLabel(7,"Rejected due to pileup");
 232   fHistNEvents->GetXaxis()->SetBinLabel(8,"Total no. of candidate");
 233   fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of cand wo bitmask");
 234   fHistNEvents->GetXaxis()->SetBinLabel(10,"D after cuts (No PID)");
 235   fHistNEvents->GetXaxis()->SetBinLabel(11,"D after cuts + PID)");
 236   fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
 237   fHistNEvents->Sumw2();
 238   fHistNEvents->SetMinimum(0);
 239   fOutput->Add(fHistNEvents);
 240
 241   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.,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
 242
 243   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.,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
 244
 245   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,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
 246
 247   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.,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
 248
 249   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.,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
 250
 251   fOutput->Add(fPtVsMassVsMult);
 252   fOutput->Add(fPtVsMassVsMultUncorr);
 253   fOutput->Add(fPtVsMassVsMultNoPid);
 254   fOutput->Add(fPtVsMassVsMultPart);
 255   fOutput->Add(fPtVsMassVsMultAntiPart);
 256
 257   if(fDoImpPar) CreateImpactParameterHistos();
 258
 259   fCounter = new AliNormalizationCounter("NormCounterCorrMult");
 260   fCounter->SetStudyMultiplicity(kTRUE,1.);
 261   fCounter->Init();
 262
 263   fCounterU = new AliNormalizationCounter("NormCounterUnCorrMult");
 264   fCounterU->SetStudyMultiplicity(kTRUE,1.);
 265   fCounterU->Init();
 266
 267   fOutputCounters = new TList();
 268   fOutputCounters->SetOwner();
 269   fOutputCounters->SetName("OutputCounters");
 270   fOutputCounters->Add(fCounter);
 271   fOutputCounters->Add(fCounterU);
 272
 273   PostData(1,fOutput);
 274   PostData(2,fListCuts);
 275   PostData(3,fOutputCounters);
 276
 277   return;
 278 }
 279
 280
 281
 282 //________________________________________________________________________
 283 void AliAnalysisTaskSEDvsMultiplicity::UserExec(Option_t */*option*/)
 284 {
 285   // Execute analysis for current event:
 286   // heavy flavor candidates association to MC truth
 287
 288   printf("UserExec\n");
 289
 290   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 291
 292   //  AliAODTracklets* tracklets = aod->GetTracklets();
 293   //Int_t ntracklets = tracklets->GetNumberOfTracklets();
 294
 295
 296   TClonesArray *arrayCand = 0;
 297   TString arrayName="";
 298   UInt_t pdgDau[3];
 299   Int_t nDau=0;
 300   Int_t selbit=0;
 301   if(fPdgMeson==411){
 302     arrayName="Charm3Prong";
 303     pdgDau[0]=211; pdgDau[1]=321; pdgDau[2]=211;
 304     nDau=3;
 305     selbit=AliRDHFCuts::kDplusCuts;
 306   }else if(fPdgMeson==421){
 307     arrayName="D0toKpi";
 308     pdgDau[0]=211; pdgDau[1]=321; pdgDau[2]=0;
 309     nDau=2;
 310     selbit=AliRDHFCuts::kD0toKpiCuts;
 311   }else if(fPdgMeson==413){
 312     arrayName="Dstar";
 313     pdgDau[0]=321; pdgDau[1]=211; pdgDau[2]=211;
 314     nDau=3;
 315     selbit=AliRDHFCuts::kDstarCuts;
 316   }
 317
 318   if(!aod && AODEvent() && IsStandardAOD()) {
 319     // In case there is an AOD handler writing a standard AOD, use the AOD
 320     // event in memory rather than the input (ESD) event.
 321     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 322      // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 323      // have to taken from the AOD event hold by the AliAODExtension
 324     AliAODHandler* aodHandler = (AliAODHandler*)
 325       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 326     if(aodHandler->GetExtensions()) {
 327       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 328       AliAODEvent *aodFromExt = ext->GetAOD();
 329       arrayCand=(TClonesArray*)aodFromExt->GetList()->FindObject(arrayName.Data());
 330     }
 331   } else if(aod) {
 332     arrayCand=(TClonesArray*)aod->GetList()->FindObject(arrayName.Data());
 333   }
 334
 335   if(!aod || !arrayCand) {
 336     printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: Charm3Prong branch not found!\n");
 337     return;
 338   }
 339
 340
 341
 342   // fix for temporary bug in ESDfilter
 343   // the AODs with null vertex pointer didn't pass the PhysSel
 344   if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return;
 345
 346   Int_t countTreta1=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.,1.);
 347   Int_t countTr=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.6,1.6);
 348
 349   fCounterU->StoreEvent(aod,fRDCutsAnalysis,fReadMC,countTreta1);
 350   fHistNEvents->Fill(0); // count event
 351
 352   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 353   TString primTitle = vtx1->GetTitle();
 354   Double_t countTreta1corr=countTreta1;
 355   if(vtx1 && vtx1->GetNContributors()>0){
 356     fHistNEvents->Fill(1);
 357     TProfile* estimatorAvg = GetEstimatorHistogram(aod);
 358     if(estimatorAvg){
 359       countTreta1corr=AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,countTreta1,vtx1->GetZ(),fRefMult);
 360     }
 361   }
 362
 363   fCounter->StoreEvent(aod,fRDCutsAnalysis,fReadMC,countTreta1corr);
 364
 365   Bool_t isEvSel=fRDCutsAnalysis->IsEventSelected(aod);
 366
 367   if(fRDCutsAnalysis->GetWhyRejection()==5) fHistNEvents->Fill(3);
 368   if(fRDCutsAnalysis->GetWhyRejection()==7) fHistNEvents->Fill(4);
 369   if(fRDCutsAnalysis->GetWhyRejection()==6)fHistNEvents->Fill(5);
 370   if(fRDCutsAnalysis->GetWhyRejection()==1)fHistNEvents->Fill(6);
 371
 372
 373   if(!isEvSel)return;
 374   fHistNEvents->Fill(2); // count events selected
 375
 376   ((TH2F*)(fOutput->FindObject("heta16vseta1")))->Fill(countTreta1,countTr);
 377   ((TH2F*)(fOutput->FindObject("heta1corrvseta1")))->Fill(countTreta1,countTreta1corr);
 378   ((TH2F*)(fOutput->FindObject("hNtrkvsVtxZ")))->Fill(vtx1->GetZ(),countTreta1);
 379   ((TH2F*)(fOutput->FindObject("hNtrkvsVtxZCorr")))->Fill(vtx1->GetZ(),countTreta1corr);
 380
 381
 382   TClonesArray *arrayMC=0;
 383   AliAODMCHeader *mcHeader=0;
 384
 385   // load MC particles
 386   if(fReadMC){
 387
 388     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 389     if(!arrayMC) {
 390       printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: MC particles branch not found!\n");
 391       return;
 392     }
 393     // load MC header
 394     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 395     if(!mcHeader) {
 396       printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: MC header branch not found!\n");
 397       return;
 398      }
 399   }
 400
 401   Int_t nCand = arrayCand->GetEntriesFast();
 402   Int_t nSelectedNoPID=0,nSelectedPID=0,nSelectedInMassPeak=0;
 403   Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
 404   Double_t mDplusPDG = TDatabasePDG::Instance()->GetParticle(411)->Mass();
 405   Double_t mDstarPDG = TDatabasePDG::Instance()->GetParticle(413)->Mass();
 406
 407   for (Int_t iCand = 0; iCand < nCand; iCand++) {
 408     AliAODRecoDecayHF *d = (AliAODRecoDecayHF*)arrayCand->UncheckedAt(iCand);
 409     fHistNEvents->Fill(7);
 410     if(fUseBit && !d->HasSelectionBit(selbit)){
 411       fHistNEvents->Fill(8);
 412       continue;
 413     }
 414
 415     Double_t ptCand = d->Pt();
 416     Double_t rapid=d->Y(fPdgMeson);
 417     Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
 418     Int_t passAllCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate,aod);
 419     Int_t passTopolCuts=fRDCutsAnalysis->GetIsSelectedCuts();
 420     if(passTopolCuts==0) continue;
 421     nSelectedNoPID++;
 422     fHistNEvents->Fill(9);
 423     if(passAllCuts){
 424       nSelectedPID++;
 425       fHistNEvents->Fill(10);
 426     }
 427     Bool_t isPrimary=kTRUE;
 428     Int_t labD=-1;
 429     Double_t trueImpParXY=9999.;
 430     Double_t impparXY=d->ImpParXY()*10000.;
 431     Double_t dlen=0.1; //FIXME
 432     Double_t mass[2];
 433     if(fPdgMeson==411){
 434       mass[0]=d->InvMass(nDau,pdgDau);
 435       mass[1]=-1.;
 436       if(TMath::Abs(mass[0]-mDplusPDG)<0.02) nSelectedInMassPeak++; //20 MeV for now... FIXME
 437     }else if(fPdgMeson==421){
 438       UInt_t pdgdaughtersD0[2]={211,321};//pi,K
 439       UInt_t pdgdaughtersD0bar[2]={321,211};//K,pi
 440       mass[0]=d->InvMass(2,pdgdaughtersD0);
 441       mass[1]=d->InvMass(2,pdgdaughtersD0bar);
 442       if(TMath::Abs(mass[0]-mD0PDG)<0.02 || TMath::Abs(mass[1]-mD0PDG)<0.02 ) nSelectedInMassPeak++; //20 MeV for now... FIXME
 443     }else if(fPdgMeson==413){
 444       // FIXME
 445       AliAODRecoCascadeHF* temp = (AliAODRecoCascadeHF*)d;
 446       mass[0]=temp->DeltaInvMass();
 447       mass[1]=-1.;
 448       if(TMath::Abs(mass[0]-(mDstarPDG-mD0PDG))<0.001) nSelectedInMassPeak++; //1 MeV for now... FIXME
 449     }
 450     for(Int_t iHyp=0; iHyp<2; iHyp++){
 451       if(mass[iHyp]<0.) continue; // for D+ and D* we have 1 mass hypothesis
 452       Double_t invMass=mass[iHyp];
 453       Double_t arrayForSparse[5]={invMass,ptCand,impparXY,dlen,countTreta1corr};
 454
 455       if(fReadMC){
 456
 457         labD = d->MatchToMC(fPdgMeson,arrayMC,nDau,(Int_t*)pdgDau);
 458         Bool_t fillHisto=fDoImpPar;
 459         if(labD>=0){
 460           AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
 461           Int_t code=partD->GetPdgCode();
 462           if(CheckOrigin(arrayMC,partD)==5) isPrimary=kFALSE;
 463           if(code<0 && iHyp==0) fillHisto=kFALSE;
 464           if(code>0 && iHyp==1) fillHisto=kFALSE;
 465           if(!isPrimary){
 466             if(fPdgMeson==411){
 467               trueImpParXY=AliVertexingHFUtils::GetTrueImpactParameterDplus(mcHeader,arrayMC,partD)*10000.;
 468             }else if(fPdgMeson==421){
 469               trueImpParXY=AliVertexingHFUtils::GetTrueImpactParameterDzero(mcHeader,arrayMC,partD)*10000.;
 470             }else if(fPdgMeson==413){
 471               trueImpParXY=0.; /// FIXME
 472             }
 473             Double_t arrayForSparseTrue[5]={invMass,ptCand,trueImpParXY,dlen,countTreta1corr};
 474             if(fillHisto && isFidAcc && passAllCuts){
 475               fHistMassPtImpPar[2]->Fill(arrayForSparse);
 476               fHistMassPtImpPar[3]->Fill(arrayForSparseTrue);
 477             }
 478           }else{
 479             if(fillHisto && isFidAcc && passAllCuts) fHistMassPtImpPar[1]->Fill(arrayForSparse);
 480           }
 481         }else{
 482           if(fillHisto && isFidAcc && passAllCuts)fHistMassPtImpPar[4]->Fill(arrayForSparse);
 483         }
 484         if(fPdgMeson==421){
 485           if(TMath::Abs(labD)==fPdgMeson && fMCOption==2) continue;
 486           if(TMath::Abs(labD)!=fPdgMeson && fMCOption==1) continue;
 487         }
 488       }
 489
 490       if(iHyp==0 && !(passTopolCuts&1)) continue; // candidate not passing as D0
 491       if(iHyp==1 && !(passTopolCuts&2)) continue; // candidate not passing as D0bar
 492
 493       if(isFidAcc){
 494         fPtVsMassVsMultNoPid->Fill(countTreta1corr,invMass,ptCand);
 495         if(passAllCuts){
 496           fPtVsMassVsMult->Fill(countTreta1corr,invMass,ptCand);
 497           fPtVsMassVsMultUncorr->Fill(countTreta1,invMass,ptCand);
 498           // Add separation between part antipart
 499           if(fPdgMeson==411){
 500             if(d->GetCharge()>0) fPtVsMassVsMultPart->Fill(countTreta1corr,invMass,ptCand);
 501             else fPtVsMassVsMultAntiPart->Fill(countTreta1corr,invMass,ptCand);
 502           }else if(fPdgMeson==421){
 503             if(passTopolCuts&1) fPtVsMassVsMultPart->Fill(countTreta1corr,invMass,ptCand);
 504             if(passTopolCuts&2) fPtVsMassVsMultAntiPart->Fill(countTreta1corr,invMass,ptCand);
 505           }else if(fPdgMeson==413){
 506             // FIXME ADD Dstar!!!!!!!!
 507           }
 508
 509
 510           if(fDoImpPar){
 511             fHistMassPtImpPar[0]->Fill(arrayForSparse);
 512           }
 513
 514         }
 515       }
 516
 517     }
 518   }
 519   fCounter->StoreCandidates(aod,nSelectedNoPID,kTRUE);
 520   fCounter->StoreCandidates(aod,nSelectedPID,kFALSE);
 521   if(nSelectedPID>0)  ((TH2F*)(fOutput->FindObject("hspdmultCand")))->Fill(countTreta1corr);
 522   if(nSelectedInMassPeak)  ((TH2F*)(fOutput->FindObject("hspdmultD")))->Fill(countTreta1corr);
 523
 524   PostData(1,fOutput);
 525   PostData(2,fListCuts);
 526   PostData(3,fOutput);
 527
 528
 529   return;
 530 }
 531 //________________________________________________________________________
 532 void AliAnalysisTaskSEDvsMultiplicity::CreateImpactParameterHistos(){
 533   // Histos for impact paramter study
 534   // mass . pt , impact parameter , decay length , multiplicity
 535
 536   Int_t nbins[5]={fNMassBins,200,fNImpParBins,50,100};
 537   Double_t xmin[5]={fLowmasslimit,0.,fLowerImpPar,0.,0.};
 538   Double_t xmax[5]={fUpmasslimit,20.,fHigherImpPar,1.,100.};
 539
 540   fHistMassPtImpPar[0]=new THnSparseF("hMassPtImpParAll",
 541                                         "Mass vs. pt vs.imppar - All",
 542                                         5,nbins,xmin,xmax);
 543   fHistMassPtImpPar[1]=new THnSparseF("hMassPtImpParPrompt",
 544                                         "Mass vs. pt vs.imppar - promptD",
 545                                         5,nbins,xmin,xmax);
 546   fHistMassPtImpPar[2]=new THnSparseF("hMassPtImpParBfeed",
 547                                         "Mass vs. pt vs.imppar - DfromB",
 548                                         5,nbins,xmin,xmax);
 549   fHistMassPtImpPar[3]=new THnSparseF("hMassPtImpParTrueBfeed",
 550                                         "Mass vs. pt vs.true imppar -DfromB",
 551                                         5,nbins,xmin,xmax);
 552   fHistMassPtImpPar[4]=new THnSparseF("hMassPtImpParBkg",
 553                                         "Mass vs. pt vs.imppar - backgr.",
 554                                         5,nbins,xmin,xmax);
 555   for(Int_t i=0; i<5;i++){
 556     fOutput->Add(fHistMassPtImpPar[i]);
 557   }
 558 }
 559
 560 //________________________________________________________________________
 561 void AliAnalysisTaskSEDvsMultiplicity::Terminate(Option_t */*option*/)
 562 {
 563   // Terminate analysis
 564   //
 565   if(fDebug > 1) printf("AnalysisTaskSEDvsMultiplicity: Terminate() \n");
 566
 567   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 568   if (!fOutput) {
 569     printf("ERROR: fOutput not available\n");
 570     return;
 571   }
 572   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
 573
 574
 575
 576   return;
 577 }
 578 //_________________________________________________________________________________________________
 579 Int_t AliAnalysisTaskSEDvsMultiplicity::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
 580   //
 581   // checking whether the mother of the particles come from a charm or a bottom quark
 582   //
 583
 584   Int_t pdgGranma = 0;
 585   Int_t mother = 0;
 586   mother = mcPartCandidate->GetMother();
 587   Int_t istep = 0;
 588   Int_t abspdgGranma =0;
 589   Bool_t isFromB=kFALSE;
 590   Bool_t isQuarkFound=kFALSE;
 591   while (mother >0 ){
 592     istep++;
 593     AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
 594     if (mcGranma){
 595       pdgGranma = mcGranma->GetPdgCode();
 596       abspdgGranma = TMath::Abs(pdgGranma);
 597       if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
 598         isFromB=kTRUE;
 599       }
 600       if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
 601       mother = mcGranma->GetMother();
 602     }else{
 603       AliError("Failed casting the mother particle!");
 604       break;
 605     }
 606   }
 607
 608   if(isFromB) return 5;
 609   else return 4;
 610 }
 611
 612
 613
 614 //____________________________________________________________________________
 615 TProfile* AliAnalysisTaskSEDvsMultiplicity::GetEstimatorHistogram(const AliVEvent* event){
 616   // Get Estimator Histogram from period event->GetRunNumber();
 617   //
 618   // If you select SPD tracklets in |eta|<1 you should use type == 1
 619   //
 620
 621   Int_t runNo  = event->GetRunNumber();
 622   Int_t period = -1;   // 0-LHC10b, 1-LHC10c, 2-LHC10d, 3-LHC10e
 623   if(runNo>114930 && runNo<117223) period = 0;
 624   if(runNo>119158 && runNo<120830) period = 1;
 625   if(runNo>122373 && runNo<126438) period = 2;
 626   if(runNo>127711 && runNo<130841) period = 3;
 627   if(period<0 || period>3) return 0;
 628
 629   return fMultEstimatorAvg[period];
 630 }
 631