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     fLowmasslimit = lowlimit;
 171     fUpmasslimit = 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     if(!isFidAcc) continue;
 436     Int_t passAllCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate,aod);
 437     Int_t passTopolCuts=fRDCutsAnalysis->GetIsSelectedCuts();
 438     if(passTopolCuts==0) continue;
 439     nSelectedNoPID++;
 440     fHistNEvents->Fill(9);
 441     if(passAllCuts){
 442       nSelectedPID++;
 443       fHistNEvents->Fill(10);
 444     }
 445     Bool_t isPrimary=kTRUE;
 446     Int_t labD=-1;
 447     Double_t trueImpParXY=9999.;
 448     Double_t impparXY=d->ImpParXY()*10000.;
 449     Double_t dlen=0.1; //FIXME
 450     Double_t mass[2];
 451     if(fPdgMeson==411){
 452       mass[0]=d->InvMass(nDau,pdgDau);
 453       mass[1]=-1.;
 454       if(TMath::Abs(mass[0]-mDplusPDG)<0.02) nSelectedInMassPeak++; //20 MeV for now... FIXME
 455     }else if(fPdgMeson==421){
 456       UInt_t pdgdaughtersD0[2]={211,321};//pi,K
 457       UInt_t pdgdaughtersD0bar[2]={321,211};//K,pi
 458       mass[0]=d->InvMass(2,pdgdaughtersD0);
 459       mass[1]=d->InvMass(2,pdgdaughtersD0bar);
 460       if(TMath::Abs(mass[0]-mD0PDG)<0.02 || TMath::Abs(mass[1]-mD0PDG)<0.02 ) nSelectedInMassPeak++; //20 MeV for now... FIXME
 461     }else if(fPdgMeson==413){
 462       // FIXME
 463       AliAODRecoCascadeHF* temp = (AliAODRecoCascadeHF*)d;
 464       mass[0]=temp->DeltaInvMass();
 465       mass[1]=-1.;
 466       if(TMath::Abs(mass[0]-(mDstarPDG-mD0PDG))<0.001) nSelectedInMassPeak++; //1 MeV for now... FIXME
 467     }
 468     for(Int_t iHyp=0; iHyp<2; iHyp++){
 469       if(mass[iHyp]<0.) continue; // for D+ and D* we have 1 mass hypothesis
 470       Double_t invMass=mass[iHyp];
 471       Double_t arrayForSparse[5]={invMass,ptCand,impparXY,dlen,countTreta1corr};
 472
 473       if(fReadMC){
 474
 475         labD = d->MatchToMC(fPdgMeson,arrayMC,nDau,(Int_t*)pdgDau);
 476         Bool_t fillHisto=fDoImpPar;
 477         if(labD>=0){
 478           AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
 479           Int_t code=partD->GetPdgCode();
 480           if(CheckOrigin(arrayMC,partD)==5) isPrimary=kFALSE;
 481           if(code<0 && iHyp==0) fillHisto=kFALSE;
 482           if(code>0 && iHyp==1) fillHisto=kFALSE;
 483           if(!isPrimary){
 484             if(fPdgMeson==411){
 485               trueImpParXY=AliVertexingHFUtils::GetTrueImpactParameterDplus(mcHeader,arrayMC,partD)*10000.;
 486             }else if(fPdgMeson==421){
 487               trueImpParXY=AliVertexingHFUtils::GetTrueImpactParameterDzero(mcHeader,arrayMC,partD)*10000.;
 488             }else if(fPdgMeson==413){
 489               trueImpParXY=0.; /// FIXME
 490             }
 491             Double_t arrayForSparseTrue[5]={invMass,ptCand,trueImpParXY,dlen,countTreta1corr};
 492             if(fillHisto && passAllCuts){
 493               fHistMassPtImpPar[2]->Fill(arrayForSparse);
 494               fHistMassPtImpPar[3]->Fill(arrayForSparseTrue);
 495             }
 496           }else{
 497             if(fillHisto && passAllCuts) fHistMassPtImpPar[1]->Fill(arrayForSparse);
 498           }
 499         }else{
 500           if(fillHisto && passAllCuts)fHistMassPtImpPar[4]->Fill(arrayForSparse);
 501         }
 502         if(fPdgMeson==421){
 503           if(TMath::Abs(labD)==fPdgMeson && fMCOption==2) continue;
 504           if(TMath::Abs(labD)!=fPdgMeson && fMCOption==1) continue;
 505         }
 506       }
 507
 508       if(fPdgMeson==421){
 509         if(iHyp==0 && !(passTopolCuts&1)) continue; // candidate not passing as D0
 510         if(iHyp==1 && !(passTopolCuts&2)) continue; // candidate not passing as D0bar
 511       }
 512
 513       fPtVsMassVsMultNoPid->Fill(countTreta1corr,invMass,ptCand);
 514       if(passAllCuts){
 515         fPtVsMassVsMult->Fill(countTreta1corr,invMass,ptCand);
 516         fPtVsMassVsMultUncorr->Fill(countTreta1,invMass,ptCand);
 517         // Add separation between part antipart
 518         if(fPdgMeson==411){
 519           if(d->GetCharge()>0) fPtVsMassVsMultPart->Fill(countTreta1corr,invMass,ptCand);
 520           else fPtVsMassVsMultAntiPart->Fill(countTreta1corr,invMass,ptCand);
 521         }else if(fPdgMeson==421){
 522           if(passTopolCuts&1) fPtVsMassVsMultPart->Fill(countTreta1corr,invMass,ptCand);
 523           if(passTopolCuts&2) fPtVsMassVsMultAntiPart->Fill(countTreta1corr,invMass,ptCand);
 524         }else if(fPdgMeson==413){
 525           // FIXME ADD Dstar!!!!!!!!
 526         }
 527
 528         if(fDoImpPar){
 529           fHistMassPtImpPar[0]->Fill(arrayForSparse);
 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,fOutputCounters);
 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
 591   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
 592   if(!fHistNEvents){
 593     printf("ERROR: fHistNEvents not available\n");
 594     return;
 595   }
 596   printf("Number of Analyzed Events = %d\n",(Int_t)fHistNEvents->GetBinContent(3));
 597
 598   return;
 599 }
 600 //_________________________________________________________________________________________________
 601 Int_t AliAnalysisTaskSEDvsMultiplicity::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
 602   //
 603   // checking whether the mother of the particles come from a charm or a bottom quark
 604   //
 605
 606   Int_t pdgGranma = 0;
 607   Int_t mother = 0;
 608   mother = mcPartCandidate->GetMother();
 609   Int_t istep = 0;
 610   Int_t abspdgGranma =0;
 611   Bool_t isFromB=kFALSE;
 612   Bool_t isQuarkFound=kFALSE;
 613   while (mother >0 ){
 614     istep++;
 615     AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
 616     if (mcGranma){
 617       pdgGranma = mcGranma->GetPdgCode();
 618       abspdgGranma = TMath::Abs(pdgGranma);
 619       if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
 620         isFromB=kTRUE;
 621       }
 622       if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
 623       mother = mcGranma->GetMother();
 624     }else{
 625       AliError("Failed casting the mother particle!");
 626       break;
 627     }
 628   }
 629
 630   if(isFromB) return 5;
 631   else return 4;
 632 }
 633
 634
 635
 636 //____________________________________________________________________________
 637 TProfile* AliAnalysisTaskSEDvsMultiplicity::GetEstimatorHistogram(const AliVEvent* event){
 638   // Get Estimator Histogram from period event->GetRunNumber();
 639   //
 640   // If you select SPD tracklets in |eta|<1 you should use type == 1
 641   //
 642
 643   Int_t runNo  = event->GetRunNumber();
 644   Int_t period = -1;   // 0-LHC10b, 1-LHC10c, 2-LHC10d, 3-LHC10e
 645   if(runNo>114930 && runNo<117223) period = 0;
 646   if(runNo>119158 && runNo<120830) period = 1;
 647   if(runNo>122373 && runNo<126438) period = 2;
 648   if(runNo>127711 && runNo<130841) period = 3;
 649   if(period<0 || period>3) return 0;
 650
 651   return fMultEstimatorAvg[period];
 652 }
 653