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   Double_t countTreta1corr=countTreta1;
 353   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 354   if(vtx1){
 355     if(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
 364   fCounter->StoreEvent(aod,fRDCutsAnalysis,fReadMC,countTreta1corr);
 365
 366   Bool_t isEvSel=fRDCutsAnalysis->IsEventSelected(aod);
 367
 368   if(fRDCutsAnalysis->GetWhyRejection()==5) fHistNEvents->Fill(3);
 369   if(fRDCutsAnalysis->GetWhyRejection()==7) fHistNEvents->Fill(4);
 370   if(fRDCutsAnalysis->GetWhyRejection()==6)fHistNEvents->Fill(5);
 371   if(fRDCutsAnalysis->GetWhyRejection()==1)fHistNEvents->Fill(6);
 372
 373
 374   if(!isEvSel)return;
 375   fHistNEvents->Fill(2); // count events selected
 376
 377   ((TH2F*)(fOutput->FindObject("heta16vseta1")))->Fill(countTreta1,countTr);
 378   ((TH2F*)(fOutput->FindObject("heta1corrvseta1")))->Fill(countTreta1,countTreta1corr);
 379   ((TH2F*)(fOutput->FindObject("hNtrkvsVtxZ")))->Fill(vtx1->GetZ(),countTreta1);
 380   ((TH2F*)(fOutput->FindObject("hNtrkvsVtxZCorr")))->Fill(vtx1->GetZ(),countTreta1corr);
 381
 382
 383   TClonesArray *arrayMC=0;
 384   AliAODMCHeader *mcHeader=0;
 385
 386   // load MC particles
 387   if(fReadMC){
 388
 389     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 390     if(!arrayMC) {
 391       printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: MC particles branch not found!\n");
 392       return;
 393     }
 394     // load MC header
 395     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 396     if(!mcHeader) {
 397       printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: MC header branch not found!\n");
 398       return;
 399      }
 400   }
 401
 402   Int_t nCand = arrayCand->GetEntriesFast();
 403   Int_t nSelectedNoPID=0,nSelectedPID=0,nSelectedInMassPeak=0;
 404   Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
 405   Double_t mDplusPDG = TDatabasePDG::Instance()->GetParticle(411)->Mass();
 406   Double_t mDstarPDG = TDatabasePDG::Instance()->GetParticle(413)->Mass();
 407
 408   for (Int_t iCand = 0; iCand < nCand; iCand++) {
 409     AliAODRecoDecayHF *d = (AliAODRecoDecayHF*)arrayCand->UncheckedAt(iCand);
 410     fHistNEvents->Fill(7);
 411     if(fUseBit && !d->HasSelectionBit(selbit)){
 412       fHistNEvents->Fill(8);
 413       continue;
 414     }
 415
 416     Double_t ptCand = d->Pt();
 417     Double_t rapid=d->Y(fPdgMeson);
 418     Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
 419     Int_t passAllCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate,aod);
 420     Int_t passTopolCuts=fRDCutsAnalysis->GetIsSelectedCuts();
 421     if(passTopolCuts==0) continue;
 422     nSelectedNoPID++;
 423     fHistNEvents->Fill(9);
 424     if(passAllCuts){
 425       nSelectedPID++;
 426       fHistNEvents->Fill(10);
 427     }
 428     Bool_t isPrimary=kTRUE;
 429     Int_t labD=-1;
 430     Double_t trueImpParXY=9999.;
 431     Double_t impparXY=d->ImpParXY()*10000.;
 432     Double_t dlen=0.1; //FIXME
 433     Double_t mass[2];
 434     if(fPdgMeson==411){
 435       mass[0]=d->InvMass(nDau,pdgDau);
 436       mass[1]=-1.;
 437       if(TMath::Abs(mass[0]-mDplusPDG)<0.02) nSelectedInMassPeak++; //20 MeV for now... FIXME
 438     }else if(fPdgMeson==421){
 439       UInt_t pdgdaughtersD0[2]={211,321};//pi,K
 440       UInt_t pdgdaughtersD0bar[2]={321,211};//K,pi
 441       mass[0]=d->InvMass(2,pdgdaughtersD0);
 442       mass[1]=d->InvMass(2,pdgdaughtersD0bar);
 443       if(TMath::Abs(mass[0]-mD0PDG)<0.02 || TMath::Abs(mass[1]-mD0PDG)<0.02 ) nSelectedInMassPeak++; //20 MeV for now... FIXME
 444     }else if(fPdgMeson==413){
 445       // FIXME
 446       AliAODRecoCascadeHF* temp = (AliAODRecoCascadeHF*)d;
 447       mass[0]=temp->DeltaInvMass();
 448       mass[1]=-1.;
 449       if(TMath::Abs(mass[0]-(mDstarPDG-mD0PDG))<0.001) nSelectedInMassPeak++; //1 MeV for now... FIXME
 450     }
 451     for(Int_t iHyp=0; iHyp<2; iHyp++){
 452       if(mass[iHyp]<0.) continue; // for D+ and D* we have 1 mass hypothesis
 453       Double_t invMass=mass[iHyp];
 454       Double_t arrayForSparse[5]={invMass,ptCand,impparXY,dlen,countTreta1corr};
 455
 456       if(fReadMC){
 457
 458         labD = d->MatchToMC(fPdgMeson,arrayMC,nDau,(Int_t*)pdgDau);
 459         Bool_t fillHisto=fDoImpPar;
 460         if(labD>=0){
 461           AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
 462           Int_t code=partD->GetPdgCode();
 463           if(CheckOrigin(arrayMC,partD)==5) isPrimary=kFALSE;
 464           if(code<0 && iHyp==0) fillHisto=kFALSE;
 465           if(code>0 && iHyp==1) fillHisto=kFALSE;
 466           if(!isPrimary){
 467             if(fPdgMeson==411){
 468               trueImpParXY=AliVertexingHFUtils::GetTrueImpactParameterDplus(mcHeader,arrayMC,partD)*10000.;
 469             }else if(fPdgMeson==421){
 470               trueImpParXY=AliVertexingHFUtils::GetTrueImpactParameterDzero(mcHeader,arrayMC,partD)*10000.;
 471             }else if(fPdgMeson==413){
 472               trueImpParXY=0.; /// FIXME
 473             }
 474             Double_t arrayForSparseTrue[5]={invMass,ptCand,trueImpParXY,dlen,countTreta1corr};
 475             if(fillHisto && isFidAcc && passAllCuts){
 476               fHistMassPtImpPar[2]->Fill(arrayForSparse);
 477               fHistMassPtImpPar[3]->Fill(arrayForSparseTrue);
 478             }
 479           }else{
 480             if(fillHisto && isFidAcc && passAllCuts) fHistMassPtImpPar[1]->Fill(arrayForSparse);
 481           }
 482         }else{
 483           if(fillHisto && isFidAcc && passAllCuts)fHistMassPtImpPar[4]->Fill(arrayForSparse);
 484         }
 485         if(fPdgMeson==421){
 486           if(TMath::Abs(labD)==fPdgMeson && fMCOption==2) continue;
 487           if(TMath::Abs(labD)!=fPdgMeson && fMCOption==1) continue;
 488         }
 489       }
 490
 491       if(iHyp==0 && !(passTopolCuts&1)) continue; // candidate not passing as D0
 492       if(iHyp==1 && !(passTopolCuts&2)) continue; // candidate not passing as D0bar
 493
 494       if(isFidAcc){
 495         fPtVsMassVsMultNoPid->Fill(countTreta1corr,invMass,ptCand);
 496         if(passAllCuts){
 497           fPtVsMassVsMult->Fill(countTreta1corr,invMass,ptCand);
 498           fPtVsMassVsMultUncorr->Fill(countTreta1,invMass,ptCand);
 499           // Add separation between part antipart
 500           if(fPdgMeson==411){
 501             if(d->GetCharge()>0) fPtVsMassVsMultPart->Fill(countTreta1corr,invMass,ptCand);
 502             else fPtVsMassVsMultAntiPart->Fill(countTreta1corr,invMass,ptCand);
 503           }else if(fPdgMeson==421){
 504             if(passTopolCuts&1) fPtVsMassVsMultPart->Fill(countTreta1corr,invMass,ptCand);
 505             if(passTopolCuts&2) fPtVsMassVsMultAntiPart->Fill(countTreta1corr,invMass,ptCand);
 506           }else if(fPdgMeson==413){
 507             // FIXME ADD Dstar!!!!!!!!
 508           }
 509
 510
 511           if(fDoImpPar){
 512             fHistMassPtImpPar[0]->Fill(arrayForSparse);
 513           }
 514
 515         }
 516       }
 517
 518     }
 519   }
 520   fCounter->StoreCandidates(aod,nSelectedNoPID,kTRUE);
 521   fCounter->StoreCandidates(aod,nSelectedPID,kFALSE);
 522   if(nSelectedPID>0)  ((TH2F*)(fOutput->FindObject("hspdmultCand")))->Fill(countTreta1corr);
 523   if(nSelectedInMassPeak)  ((TH2F*)(fOutput->FindObject("hspdmultD")))->Fill(countTreta1corr);
 524
 525   PostData(1,fOutput);
 526   PostData(2,fListCuts);
 527   PostData(3,fOutput);
 528
 529
 530   return;
 531 }
 532 //________________________________________________________________________
 533 void AliAnalysisTaskSEDvsMultiplicity::CreateImpactParameterHistos(){
 534   // Histos for impact paramter study
 535   // mass . pt , impact parameter , decay length , multiplicity
 536
 537   Int_t nbins[5]={fNMassBins,200,fNImpParBins,50,100};
 538   Double_t xmin[5]={fLowmasslimit,0.,fLowerImpPar,0.,0.};
 539   Double_t xmax[5]={fUpmasslimit,20.,fHigherImpPar,1.,100.};
 540
 541   fHistMassPtImpPar[0]=new THnSparseF("hMassPtImpParAll",
 542                                         "Mass vs. pt vs.imppar - All",
 543                                         5,nbins,xmin,xmax);
 544   fHistMassPtImpPar[1]=new THnSparseF("hMassPtImpParPrompt",
 545                                         "Mass vs. pt vs.imppar - promptD",
 546                                         5,nbins,xmin,xmax);
 547   fHistMassPtImpPar[2]=new THnSparseF("hMassPtImpParBfeed",
 548                                         "Mass vs. pt vs.imppar - DfromB",
 549                                         5,nbins,xmin,xmax);
 550   fHistMassPtImpPar[3]=new THnSparseF("hMassPtImpParTrueBfeed",
 551                                         "Mass vs. pt vs.true imppar -DfromB",
 552                                         5,nbins,xmin,xmax);
 553   fHistMassPtImpPar[4]=new THnSparseF("hMassPtImpParBkg",
 554                                         "Mass vs. pt vs.imppar - backgr.",
 555                                         5,nbins,xmin,xmax);
 556   for(Int_t i=0; i<5;i++){
 557     fOutput->Add(fHistMassPtImpPar[i]);
 558   }
 559 }
 560
 561 //________________________________________________________________________
 562 void AliAnalysisTaskSEDvsMultiplicity::Terminate(Option_t */*option*/)
 563 {
 564   // Terminate analysis
 565   //
 566   if(fDebug > 1) printf("AnalysisTaskSEDvsMultiplicity: Terminate() \n");
 567
 568   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 569   if (!fOutput) {
 570     printf("ERROR: fOutput not available\n");
 571     return;
 572   }
 573   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
 574
 575
 576
 577   return;
 578 }
 579 //_________________________________________________________________________________________________
 580 Int_t AliAnalysisTaskSEDvsMultiplicity::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
 581   //
 582   // checking whether the mother of the particles come from a charm or a bottom quark
 583   //
 584
 585   Int_t pdgGranma = 0;
 586   Int_t mother = 0;
 587   mother = mcPartCandidate->GetMother();
 588   Int_t istep = 0;
 589   Int_t abspdgGranma =0;
 590   Bool_t isFromB=kFALSE;
 591   Bool_t isQuarkFound=kFALSE;
 592   while (mother >0 ){
 593     istep++;
 594     AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
 595     if (mcGranma){
 596       pdgGranma = mcGranma->GetPdgCode();
 597       abspdgGranma = TMath::Abs(pdgGranma);
 598       if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
 599         isFromB=kTRUE;
 600       }
 601       if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
 602       mother = mcGranma->GetMother();
 603     }else{
 604       AliError("Failed casting the mother particle!");
 605       break;
 606     }
 607   }
 608
 609   if(isFromB) return 5;
 610   else return 4;
 611 }
 612
 613
 614
 615 //____________________________________________________________________________
 616 TProfile* AliAnalysisTaskSEDvsMultiplicity::GetEstimatorHistogram(const AliVEvent* event){
 617   // Get Estimator Histogram from period event->GetRunNumber();
 618   //
 619   // If you select SPD tracklets in |eta|<1 you should use type == 1
 620   //
 621
 622   Int_t runNo  = event->GetRunNumber();
 623   Int_t period = -1;   // 0-LHC10b, 1-LHC10c, 2-LHC10d, 3-LHC10e
 624   if(runNo>114930 && runNo<117223) period = 0;
 625   if(runNo>119158 && runNo<120830) period = 1;
 626   if(runNo>122373 && runNo<126438) period = 2;
 627   if(runNo>127711 && runNo<130841) period = 3;
 628   if(period<0 || period>3) return 0;
 629
 630   return fMultEstimatorAvg[period];
 631 }
 632