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