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