PWGJE/FlavourJetTasks/AliAnalysisTaskSEDmesonsFilterCJ.cxx

   1 // $Id$\r
   2 /**************************************************************************\r
   3  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *\r
   4  *                                                                        *\r
   5  * Author: The ALICE Off-line Project.                                    *\r
   6  * Contributors are mentioned in the code where appropriate.              *\r
   7  *                                                                        *\r
   8  * Permission to use, copy, modify and distribute this software and its   *\r
   9  * documentation strictly for non-commercial purposes is hereby granted   *\r
  10  * without fee, provided that the above copyright notice appears in all   *\r
  11  * copies and that both the copyright notice and this permission notice   *\r
  12  * appear in the supporting documentation. The authors make no claims     *\r
  13  * about the suitability of this software for any purpose. It is          *\r
  14  * provided "as is" without express or implied warranty.                  *\r
  15  **************************************************************************/\r
  16 //\r
  17 //\r
  18 // Task for selecting D mesons to be used as an input for D-Jet correlations\r
  19 //\r
  20 //-----------------------------------------------------------------------\r
  21 // Authors:\r
  22 // C. Bianchin (Utrecht University) chiara.bianchin@cern.ch\r
  23 // A.Grelli (Utrecht University) a.grelli@uu.nl\r
  24 // Xiaoming Zhang (LBNL)  XMZhang@lbl.gov\r
  25 //-----------------------------------------------------------------------\r
  26 \r
  27 #include <TDatabasePDG.h>\r
  28 #include <TParticle.h>\r
  29 #include <TVector3.h>\r
  30 #include "TROOT.h"\r
  31 #include <TH3F.h>\r
  32 \r
  33 #include "AliRDHFCutsDStartoKpipi.h"\r
  34 #include "AliRDHFCutsD0toKpi.h"\r
  35 #include "AliAODMCHeader.h"\r
  36 #include "AliAODHandler.h"\r
  37 #include "AliAnalysisManager.h"\r
  38 #include "AliLog.h"\r
  39 #include "AliAODVertex.h"\r
  40 #include "AliAODRecoDecay.h"\r
  41 #include "AliAODRecoCascadeHF.h"\r
  42 #include "AliAODRecoDecayHF2Prong.h"\r
  43 #include "AliESDtrack.h"\r
  44 #include "AliAODMCParticle.h"\r
  45 #include "AliAnalysisTaskSEDmesonsFilterCJ.h"\r
  46 \r
  47 ClassImp(AliAnalysisTaskSEDmesonsFilterCJ)\r
  48 \r
  49 //_____________________________________________________________________________\r
  50 AliAnalysisTaskSEDmesonsFilterCJ::AliAnalysisTaskSEDmesonsFilterCJ() :\r
  51   AliAnalysisTaskSE(),\r
  52   fUseMCInfo(kFALSE),\r
  53   fCandidateType(0),\r
  54   fCandidateName(""),\r
  55   fPDGmother(0),\r
  56   fNProngs(0),\r
  57   fBranchName(""),\r
  58   fOutput(0),\r
  59 //fOutputCandidates(0),\r
  60   fCuts(0),\r
  61   fMinMass(0.),\r
  62   fMaxMass(0.),\r
  63   fCandidateArray(0)\r
  64 //fIsSelectedArray(0)\r
  65 {\r
  66 //\r
  67 // Default constructor\r
  68 //\r
  69 \r
  70   for (Int_t i=4; i--;) fPDGdaughters[i] = 0;\r
  71   for (Int_t i=30; i--;) fSigmaD0[i] = 0.;\r
  72 }\r
  73 \r
  74 //_____________________________________________________________________________\r
  75 AliAnalysisTaskSEDmesonsFilterCJ::AliAnalysisTaskSEDmesonsFilterCJ(const char *name, AliRDHFCuts *cuts, ECandidateType candtype) :\r
  76   AliAnalysisTaskSE(name),\r
  77   fUseMCInfo(kFALSE),\r
  78   fCandidateType(candtype),\r
  79   fCandidateName(""),\r
  80   fPDGmother(0),\r
  81   fNProngs(0),\r
  82 //fPDGdaughters(),\r
  83   fBranchName(""),\r
  84   fOutput(0),\r
  85 //fOutputCandidates(0),\r
  86   fCuts(cuts),\r
  87   fMinMass(0.),\r
  88   fMaxMass(0.),\r
  89   fCandidateArray(0)\r
  90 //fIsSelectedArray(0)\r
  91 {\r
  92 //\r
  93 // Constructor. Initialization of Inputs and Outputs\r
  94 //\r
  95  \r
  96   Info("AliAnalysisTaskSEDmesonsFilterCJ","Calling Constructor");\r
  97 \r
  98   for (Int_t i=4; i--;) fPDGdaughters[i] = 0;\r
  99   for (Int_t i=30; i--;) fSigmaD0[i] = 0.;\r
 100 \r
 101   const Int_t nptbins = fCuts->GetNPtBins();\r
 102   Float_t defaultSigmaD013[13] = { 0.012, 0.012, 0.012, 0.015, 0.015, 0.018, 0.018, 0.020, 0.020, 0.030, 0.030, 0.037, 0.040 };\r
 103 \r
 104   switch (fCandidateType) {\r
 105   case 0 :\r
 106     fCandidateName = "D0";\r
 107     fPDGmother = 421;\r
 108     fNProngs = 2;\r
 109     fPDGdaughters[0] = 211;  // pi \r
 110     fPDGdaughters[1] = 321;  // K\r
 111     fPDGdaughters[2] = 0;    // empty\r
 112     fPDGdaughters[3] = 0;    // empty\r
 113     fBranchName = "D0toKpi";\r
 114     break;\r
 115   case 1 :\r
 116     fCandidateName = "Dstar";\r
 117     fPDGmother = 413;\r
 118     fNProngs = 3;\r
 119     fPDGdaughters[1] = 211; // pi soft\r
 120     fPDGdaughters[0] = 421; // D0\r
 121     fPDGdaughters[2] = 211; // pi fromD0\r
 122     fPDGdaughters[3] = 321; // K from D0\r
 123     fBranchName = "Dstar";\r
 124 \r
 125     if (nptbins<=13) {\r
 126       for (Int_t ipt=0; ipt<nptbins;ipt++) fSigmaD0[ipt] = defaultSigmaD013[ipt];\r
 127     } else {\r
 128       AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));\r
 129     }\r
 130     break;\r
 131   default :\r
 132     printf("%d not accepted!!\n",fCandidateType);\r
 133     break;\r
 134   }\r
 135 \r
 136   if (fCandidateType==kD0toKpi) SetMassLimits(0.15, fPDGmother);\r
 137   if (fCandidateType==kDstartoKpipi) SetMassLimits(0.015, fPDGmother);\r
 138 \r
 139   DefineOutput(1, TList::Class());       // histos\r
 140   DefineOutput(2, AliRDHFCuts::Class()); // my cuts\r
 141 }\r
 142 \r
 143 //_____________________________________________________________________________\r
 144 AliAnalysisTaskSEDmesonsFilterCJ::~AliAnalysisTaskSEDmesonsFilterCJ()\r
 145 {\r
 146 //\r
 147 // destructor\r
 148 //\r
 149 \r
 150   Info("~AliAnalysisTaskSEDmesonsFilterCJ","Calling Destructor");  \r
 151  \r
 152   if (fOutput) { delete fOutput; fOutput = 0; }\r
 153   if (fCuts)   { delete fCuts;   fCuts   = 0; }\r
 154   if (fCandidateArray)  { delete fCandidateArray;  fCandidateArray  = 0; }\r
 155 //if (fIsSelectedArray) { delete fIsSelectedArray; fIsSelectedArray = 0; }\r
 156 }\r
 157 \r
 158 //_____________________________________________________________________________\r
 159 void AliAnalysisTaskSEDmesonsFilterCJ::Init()\r
 160 {\r
 161 //\r
 162 // Initialization\r
 163 //\r
 164 \r
 165   if(fDebug>1) printf("AnalysisTaskSEDmesonsForJetCorrelations::Init() \n");\r
 166 \r
 167   switch (fCandidateType) {\r
 168     case 0: {\r
 169               AliRDHFCutsD0toKpi* copyfCutsDzero = new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));\r
 170               copyfCutsDzero->SetName("AnalysisCutsDzero");\r
 171               PostData(2, copyfCutsDzero);  // Post the data\r
 172             } break;\r
 173     case 1: {\r
 174               AliRDHFCutsDStartoKpipi* copyfCutsDstar = new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));\r
 175               copyfCutsDstar->SetName("AnalysisCutsDStar");\r
 176               PostData(2, copyfCutsDstar); // Post the cuts\r
 177             } break;\r
 178     default: return;\r
 179   }\r
 180 \r
 181   return;\r
 182 }\r
 183 \r
 184 //_____________________________________________________________________________\r
 185 void AliAnalysisTaskSEDmesonsFilterCJ::UserCreateOutputObjects()\r
 186 { \r
 187 //\r
 188 // output \r
 189 //\r
 190 \r
 191   Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());\r
 192   \r
 193 //OpenFile(1);\r
 194   fOutput = new TList(); fOutput->SetOwner();\r
 195   DefineHistoForAnalysis(); // define histograms\r
 196  \r
 197 //fOutputCandidates= new TList();\r
 198 //fOutputCandidates->SetOwner();\r
 199 \r
 200   fCandidateArray = new TClonesArray("AliAODRecoDecayHF",0);\r
 201   fCandidateArray->SetName(Form("fCandidateArray%s",fCandidateName.Data()));\r
 202 \r
 203 //fOutputCandidates->Add(fCandidateArray);\r
 204 //fIsSelectedArray = new TClonesArray("Int_t",0);\r
 205 //fIsSelectedArray->SetName(Form("fIsSelectedArray%s",suffix.Data()));\r
 206 //fOutputCandidates->Add(fIsSelectedArray);\r
 207 \r
 208   PostData(1, fOutput);\r
 209 //PostData(3, fOutputCandidates);\r
 210   return;\r
 211 }\r
 212 \r
 213 //_____________________________________________________________________________\r
 214 void AliAnalysisTaskSEDmesonsFilterCJ::UserExec(Option_t *)\r
 215 {\r
 216 //\r
 217 // user exec\r
 218 //\r
 219 \r
 220   // add cadidate branch\r
 221   fCandidateArray->Delete();\r
 222   if (!(InputEvent()->FindListObject(Form("fCandidateArray%s",fCandidateName.Data())))) InputEvent()->AddObject(fCandidateArray);\r
 223 \r
 224   // Load the event\r
 225   AliAODEvent *aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);\r
 226  \r
 227   TClonesArray *arrayDStartoD0pi = 0;\r
 228   if (!aodEvent && AODEvent() && IsStandardAOD()) {\r
 229 \r
 230     // In case there is an AOD handler writing a standard AOD, use the AOD \r
 231     // event in memory rather than the input (ESD) event.    \r
 232     aodEvent = dynamic_cast<AliAODEvent*>(AODEvent());\r
 233 \r
 234     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)\r
 235     // have to taken from the AOD event hold by the AliAODExtension\r
 236     AliAODHandler *aodHandler = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());\r
 237     if(aodHandler->GetExtensions()) {\r
 238       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");\r
 239       AliAODEvent *aodFromExt = ext->GetAOD();\r
 240       arrayDStartoD0pi = (TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());\r
 241     }\r
 242   } else {\r
 243     arrayDStartoD0pi = (TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data());\r
 244   }\r
 245 \r
 246   if (!arrayDStartoD0pi) {\r
 247     AliInfo(Form("Could not find array %s, skipping the event",fBranchName.Data()));\r
 248     return;\r
 249   } else {\r
 250     AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));   \r
 251   }\r
 252   \r
 253   TClonesArray* mcArray = 0x0;\r
 254   if (fUseMCInfo) {\r
 255     mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));\r
 256     if (!mcArray) {\r
 257       printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");\r
 258       return;\r
 259     }\r
 260   }\r
 261 \r
 262   //Histograms\r
 263   TH1I* hstat = (TH1I*)fOutput->FindObject("hstat");\r
 264   TH2F* hInvMassptD = (TH2F*)fOutput->FindObject("hInvMassptD");\r
 265 \r
 266   TH1F* hPtPion=0x0;\r
 267   if (fCandidateType==kDstartoKpipi) hPtPion = (TH1F*)fOutput->FindObject("hPtPion");\r
 268   hstat->Fill(0);\r
 269 \r
 270   // fix for temporary bug in ESDfilter \r
 271   // the AODs with null vertex pointer didn't pass the PhysSel\r
 272   if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;\r
 273 \r
 274   //Event selection\r
 275   Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);\r
 276 //TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();\r
 277   if (!iseventselected) return;\r
 278   hstat->Fill(1);\r
 279 \r
 280   const Int_t nD = arrayDStartoD0pi->GetEntriesFast();\r
 281   hstat->Fill(2, nD);\r
 282 \r
 283 //Int_t icountReco = 0;  // counters for efficiencies\r
 284 \r
 285   //D* and D0 prongs needed to MatchToMC method\r
 286   Int_t pdgDgDStartoD0pi[2] = { 421, 211 };  // D0,pi\r
 287   Int_t pdgDgD0toKpi[2] = { 321, 211 };      // K, pi\r
 288 \r
 289   //D0 from D0 bar\r
 290   Int_t pdgdaughtersD0[2] = { 211, 321 };     // pi,K \r
 291   Int_t pdgdaughtersD0bar[2] = { 321, 211 };  // K,pi \r
 292 \r
 293   Int_t iCand =0;\r
 294   Int_t isSelected = 0;\r
 295   AliAODRecoDecayHF *charmCand = 0;\r
 296 \r
 297   Int_t mcLabel = -9999;\r
 298   Int_t pdgMeson = 413;\r
 299   if (fCandidateType==kD0toKpi) pdgMeson = 421;\r
 300 \r
 301   for (Int_t icharm=0; icharm<nD; icharm++) {   //loop over D candidates\r
 302     charmCand = (AliAODRecoDecayHF*)arrayDStartoD0pi->At(icharm); // D candidates\r
 303     if (!charmCand) return;\r
 304 \r
 305 \r
 306     if (fUseMCInfo) { // Look in MC, try to simulate the z\r
 307       if (fCandidateType==kDstartoKpipi) {\r
 308         AliAODRecoCascadeHF *temp = (AliAODRecoCascadeHF*)charmCand;\r
 309         mcLabel = temp->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);\r
 310       }\r
 311 \r
 312       if (fCandidateType==kD0toKpi) mcLabel = charmCand->MatchToMC(421,mcArray,fNProngs,fPDGdaughters);\r
 313 \r
 314       if (mcLabel<=0) continue;        \r
 315     }\r
 316 \r
 317     // region of interest + cuts\r
 318     if (!fCuts->IsInFiducialAcceptance(charmCand->Pt(),charmCand->Y(pdgMeson))) continue;    \r
 319 \r
 320     isSelected = fCuts->IsSelected(charmCand,AliRDHFCuts::kAll,aodEvent); //selected\r
 321     if (!isSelected) continue;\r
 322 \r
 323     new ((*fCandidateArray)[iCand]) AliAODRecoDecayHF(*charmCand);\r
 324 //  new ((*fIsSelectedArray)[iCand]) Int_t(isSelected);\r
 325     iCand++;\r
 326 \r
 327     Double_t masses[2];\r
 328     Double_t ptD = charmCand->Pt();\r
 329     if (fCandidateType==kDstartoKpipi) { //D*->D0pi->Kpipi\r
 330 \r
 331       //softpion from D* decay\r
 332       AliAODRecoCascadeHF *temp = (AliAODRecoCascadeHF*)charmCand;\r
 333       AliAODTrack *track2 = (AliAODTrack*)temp->GetBachelor();  \r
 334       hstat->Fill(3);\r
 335 \r
 336       // select D* in the D0 window.\r
 337       // In the cut object window is loose to allow side bands\r
 338       Double_t mPDGD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();\r
 339 \r
 340       // retrieve the corresponding pt bin for the candidate\r
 341       // and set the expected D0 width (x3)\r
 342 //    static const Int_t n = fCuts->GetNPtBins();\r
 343       Int_t bin = fCuts->PtBin(ptD);\r
 344       if(bin<0 || bin>=fCuts->GetNPtBins()) {\r
 345         AliError(Form("Pt %.3f out of bounds",ptD));\r
 346         continue;\r
 347       }\r
 348 \r
 349       AliInfo(Form("Pt bin %d and sigma D0 %.4f",bin,fSigmaD0[bin]));\r
 350       if ((temp->InvMassD0()>=(mPDGD0-3.*fSigmaD0[bin])) && (temp->InvMassD0()<=(mPDGD0+3.*fSigmaD0[bin]))) {   \r
 351         masses[0] = temp->DeltaInvMass(); //D*\r
 352         masses[1] = 0.; //dummy for D*\r
 353 \r
 354         //D*  delta mass\r
 355         hInvMassptD->Fill(masses[0], ptD); // 2 D slice for pt bins\r
 356 \r
 357         // D* pt and soft pion pt for good candidates           \r
 358         Double_t mPDGDstar = TDatabasePDG::Instance()->GetParticle(413)->Mass();\r
 359         Double_t invmassDelta = temp->DeltaInvMass();\r
 360         if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))<0.0021) hPtPion->Fill(track2->Pt());\r
 361       }\r
 362     } //Dstar specific\r
 363 \r
 364     if (fCandidateType==kD0toKpi) { //D0->Kpi\r
 365 \r
 366       //needed quantities\r
 367       masses[0] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0); //D0\r
 368       masses[1] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0bar); //D0bar\r
 369       hstat->Fill(3);\r
 370 \r
 371       // mass vs pt\r
 372       if (isSelected==1 || isSelected==3) hInvMassptD->Fill(masses[0],ptD);\r
 373       if (isSelected>=2) hInvMassptD->Fill(masses[1],ptD);\r
 374     } //D0 specific\r
 375 \r
 376     charmCand = 0;  \r
 377   } // end of D cand loop\r
 378 \r
 379 //AliDebug(2, Form("Found %i Reco particles that are D*!!",icountReco));\r
 380 \r
 381 //PostData(1,fOutput);\r
 382 //PostData(3,fOutputCandidates);\r
 383 \r
 384   return;\r
 385 }\r
 386 \r
 387 //_____________________________________________________________________________\r
 388 void AliAnalysisTaskSEDmesonsFilterCJ::Terminate(Option_t*)\r
 389 {\r
 390 // The Terminate() function is the last function to be called during\r
 391 // a query. It always runs on the client, it can be used to present\r
 392 // the results graphically or save the results to file.\r
 393   \r
 394   Info("Terminate"," terminate");\r
 395   AliAnalysisTaskSE::Terminate();\r
 396 \r
 397   fOutput = dynamic_cast<TList*>(GetOutputData(1));\r
 398   if (!fOutput) {\r
 399     printf("ERROR: fOutput not available\n");\r
 400     return;\r
 401   }\r
 402 \r
 403   return;\r
 404 }\r
 405 \r
 406 //_____________________________________________________________________________\r
 407 void AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits(Double_t range, Int_t pdg)\r
 408 {\r
 409 //\r
 410 // AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits\r
 411 //\r
 412 \r
 413   Float_t mass = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg))->Mass();\r
 414 \r
 415   // compute the Delta mass for the D*\r
 416   if (fCandidateType==kDstartoKpipi) mass -= TDatabasePDG::Instance()->GetParticle(421)->Mass();\r
 417 //cout << "mass ---------------" << endl;\r
 418 \r
 419   fMinMass = mass - range;\r
 420   fMaxMass = mass + range;\r
 421 \r
 422   AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));\r
 423   if ((fMinMass<0.) || (fMaxMass<=0.) || (fMaxMass<fMinMass)) AliFatal("Wrong mass limits!\n");\r
 424 \r
 425   return;\r
 426 }\r
 427 \r
 428 //_____________________________________________________________________________\r
 429 void AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits(Double_t lowlimit, Double_t uplimit)\r
 430 {\r
 431 //\r
 432 // AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits\r
 433 //\r
 434 \r
 435   if (uplimit>lowlimit) {\r
 436     fMinMass = lowlimit;\r
 437     fMaxMass = uplimit;\r
 438   } else {\r
 439     printf("Error! Lower limit larger than upper limit!\n");\r
 440     fMinMass = uplimit - uplimit*0.2;\r
 441     fMaxMass = uplimit;\r
 442   }\r
 443 \r
 444   return;\r
 445 }\r
 446 \r
 447 //_____________________________________________________________________________\r
 448 Bool_t AliAnalysisTaskSEDmesonsFilterCJ::SetD0WidthForDStar(Int_t nptbins, Float_t *width)\r
 449 {\r
 450 //\r
 451 // AliAnalysisTaskSEDmesonsFilterCJ::SetD0WidthForDStar\r
 452 //\r
 453 \r
 454   if (nptbins>30) {\r
 455     AliInfo("Maximum number of bins allowed is 30!");\r
 456     return kFALSE;\r
 457   }\r
 458 \r
 459   if (!width) return kFALSE;\r
 460   for (Int_t ipt=0; ipt<nptbins; ipt++) fSigmaD0[ipt]=width[ipt];\r
 461 \r
 462   return kTRUE;\r
 463 }\r
 464 \r
 465 \r
 466 //_____________________________________________________________________________\r
 467 Bool_t AliAnalysisTaskSEDmesonsFilterCJ::DefineHistoForAnalysis()\r
 468 {\r
 469 //\r
 470 // AliAnalysisTaskSEDmesonsFilterCJ::DefineHistoForAnalysis\r
 471 //\r
 472 \r
 473   // Statistics \r
 474   TH1I* hstat = new TH1I("hstat","Statistics",6,-0.5,5.5);\r
 475   hstat->GetXaxis()->SetBinLabel(1, "N ev anal");\r
 476   hstat->GetXaxis()->SetBinLabel(2, "N ev sel");\r
 477   hstat->GetXaxis()->SetBinLabel(3, "N cand");\r
 478   hstat->GetXaxis()->SetBinLabel(4, "N cand sel cuts");\r
 479   hstat->GetXaxis()->SetBinLabel(5, "N jets");\r
 480   hstat->GetXaxis()->SetBinLabel(6, "N cand in jet");\r
 481 /*if(fUseMCInfo) {\r
 482     hstat->GetXaxis()->SetBinLabel(7,"N D");\r
 483     hstat->GetXaxis()->SetBinLabel(8,"N D in jet");\r
 484   }*/\r
 485   hstat->SetNdivisions(1);\r
 486   fOutput->Add(hstat);\r
 487 \r
 488   // Invariant mass related histograms\r
 489   const Int_t nbinsmass = 200;\r
 490   TH2F *hInvMass = new TH2F("hInvMassptD", "D invariant mass distribution", nbinsmass, fMinMass, fMaxMass, 100, 0., 50.);\r
 491   hInvMass->SetStats(kTRUE);\r
 492   hInvMass->GetXaxis()->SetTitle("mass (GeV/c)");\r
 493   hInvMass->GetYaxis()->SetTitle("p_{T} (GeV/c)");\r
 494   fOutput->Add(hInvMass);\r
 495   \r
 496   if (fCandidateType==kDstartoKpipi) {\r
 497     TH1F* hPtPion = new TH1F("hPtPion", "Primary pions candidates pt", 500, 0., 10.);\r
 498     hPtPion->SetStats(kTRUE);\r
 499     hPtPion->GetXaxis()->SetTitle("p_{T} (GeV/c)");\r
 500     hPtPion->GetYaxis()->SetTitle("entries");\r
 501     fOutput->Add(hPtPion);\r
 502   }\r
 503 \r
 504   return kTRUE; \r
 505 }\r